source: trunk/GSASIIscriptable.py @ 4015

Last change on this file since 4015 was 4015, checked in by toby, 2 years ago

add setHAPvalues to scriptable; disable import G(r) because broken

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 224.9 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3########### SVN repository information ###################
4# $Date: 2019-06-02 18:26:06 +0000 (Sun, 02 Jun 2019) $
5# $Author: toby $
6# $Revision: 4015 $
7# $URL: trunk/GSASIIscriptable.py $
8# $Id: GSASIIscriptable.py 4015 2019-06-02 18:26:06Z toby $
9########### SVN repository information ###################
10#
11"""
12*GSASIIscriptable: Scripting Interface*
13=======================================
14
15Routines to use an increasing amount of GSAS-II's capabilities from scripts,
16without use of the graphical user interface (GUI). GSASIIscriptable can create and access
17GSAS-II project (.gpx) files and can directly perform image handling and refinements. 
18The module defines wrapper classes (inheriting from :class:`G2ObjectWrapper`) for a growing number
19of data tree items.
20
21GSASIIscriptable can be used in two ways. It offers a command-line mode
22(see :ref:`CommandlineInterface`) that
23provides access a number of features without writing Python scripts
24via shell/batch commands. The more powerful mode of GSASIIscriptable is
25use is through Python scripts that
26call the module's application interface (API), see API summary that follows or the :ref:`API`
27section.
28
29==================================================
30Application Interface (API) Summary
31==================================================
32This section of the documentation provides an overview to API, with full documentation
33in the :ref:`API` section. The typical API use will be with a Python script, such as this:
34
35.. code-block::  python
36
37    from __future__ import division, print_function
38    import os,sys
39    sys.path.insert(0,'/Users/toby/software/G2/GSASII') # needed to "find" GSAS-II modules
40    import GSASIIscriptable as G2sc
41    datadir = "/Users/Scratch/"
42    gpx = G2sc.G2Project(os.path.join(datadir,'test2.gpx'))
43    gpx.histogram(0).add_back_peak(4.5,30000,5000,0)
44    pardict = {'set': {'Sample Parameters': ['Absorption', 'Contrast', 'DisplaceX'],
45                       'Background': {'type': 'chebyschev', 'refine': True,
46                                      'peaks':[[0,True]]}}}
47    gpx.set_refinement(pardict)
48
49GSASIIscriptable can be used to setup and start sequential refinements. This script
50is used to take the single-dataset fit at the end of Step 1 of the
51`Sequential Refinement tutorial <https://subversion.xray.aps.anl.gov/pyGSAS/Tutorials/SeqRefine/SequentialTutorial.htm>`
52and turn off refinement flags, add histograms and setup the sequential fit, which is then run:
53
54.. code-block::  python
55
56    import os,sys,glob
57    sys.path.insert(0,'/Users/toby/software/G2/GSASII')
58    import GSASIIscriptable as G2sc
59    datadir = os.path.expanduser("~/Scratch/SeqTut2019Mar")
60    PathWrap = lambda fil: os.path.join(datadir,fil)
61    # load and rename project
62    gpx = G2sc.G2Project(PathWrap('7Konly.gpx'))
63    gpx.save(PathWrap('SeqRef.gpx'))
64    # turn off some variables; turn on Dijs
65    for p in gpx.phases():
66        p.set_refinements({"Cell": False})
67    gpx.phase(0).set_HAP_refinements(
68        {'Scale': False,
69         "Size": {'type':'isotropic', 'refine': False},
70         "Mustrain": {'type':'uniaxial', 'refine': False},
71         "HStrain":True,})
72    gpx.phase(1).set_HAP_refinements({'Scale': False})
73    gpx.histogram(0).clear_refinements({'Background':False,
74                     'Sample Parameters':['DisplaceX'],})
75    gpx.histogram(0).ref_back_peak(0,[])
76    gpx.phase(1).set_HAP_refinements({"HStrain":(1,1,1,0)})
77    for fil in sorted(glob.glob(PathWrap('*.fxye'))): # load in remaining fxye files
78        if '00' in fil: continue
79        gpx.add_powder_histogram(fil, PathWrap('OH_00.prm'), fmthint="GSAS powder",phases='all')
80    # copy HAP values, background, instrument params. & limits, not sample params.
81    gpx.copyHistParms(0,'all',['b','i','l'])
82    for p in gpx.phases(): p.copyHAPvalues(0,'all')
83    # setup and launch sequential fit
84    gpx.set_Controls('sequential',gpx.histograms())
85    gpx.set_Controls('cycles',10)
86    gpx.set_Controls('seqCopy',True)
87    gpx.refine() 
88
89Most functionality is provided via the objects and methods described in this section.
90
91---------------------
92:class:`G2Project`
93---------------------
94
95  All GSASIIscriptable scripts will need to create a :class:`G2Project` object
96  either for a new GSAS-II project or to read in an existing project (.gpx) file.
97  The most commonly used routines in this object are:
98
99.. tabularcolumns:: |l|p{3.5in}|
100
101==================================================    ===============================================================================================================
102method                                                Use
103==================================================    ===============================================================================================================
104:meth:`G2Project.save`                                Writes the current project to disk.
105
106:meth:`G2Project.add_powder_histogram`                Used to read in powder diffraction data into a project file.
107
108:meth:`G2Project.add_simulated_powder_histogram`      Defines a "dummy" powder diffraction data that will be simulated after a refinement step.
109
110:meth:`G2Project.add_image`                           Reads in an image into a project.
111
112:meth:`G2Project.add_phase`                           Adds a phase to a project
113
114:meth:`G2Project.add_PDF`                             Adds a PDF entry to a project (does not compute it)
115
116:meth:`G2Project.histograms`                          Provides a list of histograms in the current project, as :class:`G2PwdrData` objects
117
118:meth:`G2Project.phases`                              Provides a list of phases defined in the current project, as :class:`G2Phase` objects
119
120:meth:`G2Project.images`                              Provides a list of images in the current project, as :class:`G2Image` objects
121
122:meth:`G2Project.pdfs`                                Provides a list of PDFs in the current project, as :class:`G2PDF` objects
123
124:meth:`G2Project.seqref`                              Returns a :class:`G2SeqRefRes` object if there are Sequential Refinement results
125
126:meth:`G2Project.do_refinements`                      This is passed a list of dictionaries, where each dict defines a refinement step.
127                                                      Passing a list with a single empty dict initiates a refinement with the current
128                                                      parameters and flags. A refinement dict sets up a single refinement step
129                                                      (as described in :ref:`Project_dicts`). Also see :ref:`Refinement_recipe`.
130
131:meth:`G2Project.set_refinement`                      This is passed a single dict which is used to set parameters and flags.
132                                                      These actions can be performed also in :meth:`G2Project.do_refinements`.
133:meth:`G2Project.get_Variable`                        Retrieves the value and esd for a parameter
134:meth:`G2Project.get_Covariance`                      Retrieves values and covariance for a set of refined parameters
135:meth:`G2Project.set_Controls`                        Set overall GSAS-II control settings such as number of cycles and to set up a sequential
136                                                      fit. (Also see :meth:`G2Project.get_Controls` to read values.)
137==================================================    ===============================================================================================================
138
139---------------------
140:class:`G2Phase`
141---------------------
142
143  Another common object in GSASIIscriptable scripts is :class:`G2Phase`, used to encapsulate each phase in a project, with commonly used methods:
144
145.. tabularcolumns:: |l|p{3.5in}|
146
147==================================================    ===============================================================================================================
148method                                                Use
149==================================================    ===============================================================================================================
150:meth:`G2Phase.set_refinements`                       Provides a mechanism to set values and refinement flags for the phase. See :ref:`Phase_parameters_table`
151                                                      for more details. This information also can be supplied within a call to :meth:`G2Project.do_refinements`
152                                                      or :meth:`G2Project.set_refinement`.
153:meth:`G2Phase.clear_refinements`                     Unsets refinement flags for the phase.
154:meth:`G2Phase.set_HAP_refinements`                   Provides a mechanism to set values and refinement flags for parameters specific to both this phase and
155                                                      one of its histograms. See :ref:`HAP_parameters_table`. This information also can be supplied within
156                                                      a call to :meth:`G2Project.do_refinements` or :meth:`G2Project.set_refinement`.
157:meth:`G2Phase.clear_HAP_refinements`                 Clears refinement flags specific to both this phase and one of its histograms.
158:meth:`G2Phase.getHAPvalues`                          Returns values of parameters specific to both this phase and one of its histograms.
159:meth:`G2Phase.copyHAPvalues`                         Copies HAP settings between from one phase/histogram and to other histograms in same phase.
160:meth:`G2Phase.atoms`                                 Returns a list of atoms in the phase
161:meth:`G2Phase.atom`                                  Returns an atom from its label
162:meth:`G2Phase.histograms`                            Returns a list of histograms linked to the phase
163:meth:`G2Phase.get_cell`                              Returns unit cell parameters (also see :meth:`G2Phase.get_cell_and_esd`)
164:meth:`G2Phase.export_CIF`                            Writes a CIF for the phase
165==================================================    ===============================================================================================================
166
167---------------------
168:class:`G2PwdrData`
169---------------------
170
171  Another common object in GSASIIscriptable scripts is :class:`G2PwdrData`, which encapsulate each powder diffraction histogram in a project, with commonly used methods:
172
173.. tabularcolumns:: |l|p{3.5in}|
174
175==================================================    ===============================================================================================================
176method                                                Use
177==================================================    ===============================================================================================================
178:meth:`G2PwdrData.set_refinements`                    Provides a mechanism to set values and refinement flags for the powder histogram. See
179                                                      :ref:`Histogram_parameters_table` for details. 
180:meth:`G2PwdrData.clear_refinements`                  Unsets refinement flags for the the powder histogram.
181:meth:`G2PwdrData.residuals`                          Reports R-factors etc. for the the powder histogram (also see :meth:`G2PwdrData.get_wR`)
182:meth:`G2PwdrData.add_back_peak`                      Adds a background peak to the histogram. Also see :meth:`G2PwdrData.del_back_peak` and
183                                                      :meth:`G2PwdrData.ref_back_peak`.
184:meth:`G2PwdrData.fit_fixed_points`                   Fits background to the specified fixed points.
185:meth:`G2PwdrData.getdata`                            Provides access to the diffraction data associated with the histogram.
186:meth:`G2PwdrData.Export`                             Writes the diffraction data into a file
187==================================================    ===============================================================================================================
188
189---------------------
190:class:`G2Image`
191---------------------
192
193  When working with images, there will be a :class:`G2Image` object for each image (also see :meth:`G2Project.add_image`  and :meth:`G2Project.images`).
194
195.. tabularcolumns:: |l|p{3.5in}|
196
197==================================================    ===============================================================================================================
198method                                                Use
199==================================================    ===============================================================================================================
200:meth:`G2Image.Recalibrate`                           Invokes a recalibration fit starting from the current Image Controls calibration coefficients.
201:meth:`G2Image.Integrate`                             Invokes an image integration All parameters Image Controls will have previously been set.
202:meth:`G2Image.setControl`                            Set an Image Controls parameter in the current image.
203:meth:`G2Image.getControl`                            Return an Image Controls parameter in the current image.
204:meth:`G2Image.findControl`                           Get the names of Image Controls parameters.
205:meth:`G2Image.loadControls`                          Load controls from a .imctrl file (also see :meth:`G2Image.saveControls`).
206:meth:`G2Image.loadMasks`                             Load masks from a .immask file.
207:meth:`G2Image.setVary`                               Set a refinement flag for Image Controls parameter in the current image. (Also see :meth:`G2Image.getVary`)
208:meth:`G2Image.setCalibrant`                          Set a calibrant type (or show choices) for the current image.
209:meth:`G2Image.setControlFile`                        Set a image to be used as a background/dark/gain map image.
210==================================================    ===============================================================================================================
211
212
213---------------------
214:class:`G2PDF`
215---------------------
216
217  To work with PDF entries, object :class:`G2PDF`, encapsulates a PDF entry with methods:
218
219.. tabularcolumns:: |l|p{3.5in}|
220
221==================================================    ===============================================================================================================
222method                                                Use
223==================================================    ===============================================================================================================
224:meth:`G2PDF.export`                                   Used to write G(r), etc. as a file
225:meth:`G2PDF.calculate`                                Computes the PDF using parameters in the object
226:meth:`G2PDF.optimize`                                 Optimizes selected PDF parameters
227:meth:`G2PDF.set_background`                           Sets the histograms used for sample background, container, etc.
228:meth:`G2PDF.set_formula`                              Sets the chemical formula for the sample
229==================================================    ===============================================================================================================
230
231---------------------
232:class:`G2SeqRefRes`
233---------------------
234
235  To work with Sequential Refinement results, object :class:`G2SeqRefRes`, encapsulates the sequential refinement table with methods:
236
237.. tabularcolumns:: |l|p{3.5in}|
238
239==================================================    ===============================================================================================================
240method                                                Use
241==================================================    ===============================================================================================================
242:meth:`G2SeqRefRes.histograms`                         Provides a list of histograms used in the Sequential Refinement
243:meth:`G2SeqRefRes.get_cell_and_esd`                   Returns cell dimensions and standard uncertainies for a phase and histogram from the Sequential Refinement
244:meth:`G2SeqRefRes.get_Variable`                       Retrieves the value and esd for a parameter from a particular histogram in the Sequential Refinement
245:meth:`G2SeqRefRes.get_Covariance`                     Retrieves values and covariance for a set of refined parameters for a particular histogram
246==================================================    ===============================================================================================================
247
248----------------------
249:class:`G2AtomRecord`
250----------------------
251
252  When working with phases, :class:`G2AtomRecord` objects provide access to the contents of each atom in a phase. This provides access to "properties" that can be
253  used to get values of much of the atoms associated settings: label, type, refinement_flags, coordinates, occupancy, ranId, adp_flag, and uiso. In addition,
254  refinement_flags, occupancy and uiso can be used to set values. See the :class:`G2AtomRecord` docs and source code.
255
256.. _Refinement_dicts:
257
258=====================
259Refinement parameters
260=====================
261While scripts can be written that setup refinements by changing individual parameters
262through calls to the methods associated with objects that wrap each data tree item,
263many of these actions can be combined into fairly complex dict structures to conduct refinement
264steps. Use of these dicts is required with the :ref:`CommandlineInterface`. This section of the
265documentation describes these dicts.
266
267.. _Project_dicts:
268
269-----------------------------
270Project-level Parameter Dict
271-----------------------------
272
273As noted below (:ref:`Refinement_parameters_kinds`), there are three types of refinement parameters,
274which can be accessed individually by the objects that encapsulate individual phases and histograms
275but it will often be simplest to create a composite dictionary
276that is used at the project-level. A dict is created with keys
277"set" and "clear" that can be supplied to :meth:`G2Project.set_refinement`
278(or :meth:`G2Project.do_refinements`, see :ref:`Refinement_recipe` below) that will
279determine parameter values and will determine which parameters will be refined.
280
281The specific keys and subkeys that can be used are defined in tables
282:ref:`Histogram_parameters_table`, :ref:`Phase_parameters_table` and :ref:`HAP_parameters_table`.
283
284Note that optionally a list of histograms and/or phases can be supplied in the call to
285:meth:`G2Project.set_refinement`, but if not specified, the default is to use all defined
286phases and histograms.
287
288As an example:
289
290.. code-block::  python
291
292    pardict = {'set': { 'Limits': [0.8, 12.0],
293                       'Sample Parameters': ['Absorption', 'Contrast', 'DisplaceX'],
294                       'Background': {'type': 'chebyschev', 'refine': True,
295                                      'peaks':[[0,True],[1,1,1]] }},
296              'clear': {'Instrument Parameters': ['U', 'V', 'W']}}
297    my_project.set_refinement(pardict)
298   
299.. _Refinement_recipe:
300   
301------------------------
302Refinement recipe
303------------------------
304Building on the :ref:`Project_dicts`,
305it is possible to specify a sequence of refinement actions as a list of
306these dicts and supplying this list
307as an argument to :meth:`G2Project.do_refinements`.
308
309As an example, this code performs the same actions as in the example in the section above:
310
311.. code-block::  python
312   
313    pardict = {'set': { 'Limits': [0.8, 12.0],
314                       'Sample Parameters': ['Absorption', 'Contrast', 'DisplaceX'],
315                       'Background': {'type': 'chebyschev', 'refine': True}},
316              'clear': {'Instrument Parameters': ['U', 'V', 'W']}}
317    my_project.do_refinements([pardict])
318
319However, in addition to setting a number of parameters, this example will perform a refinement as well,
320after setting the parameters. More than one refinement can be performed by including more
321than one dict in the list.
322
323In this example, two refinement steps will be performed:
324
325.. code-block::  python
326
327    my_project.do_refinements([pardict,pardict1])
328
329
330The keys defined in the following table
331may be used in a dict supplied to :meth:`G2Project.do_refinements`. Note that keys ``histograms``
332and ``phases`` are used to limit actions to specific sets of parameters within the project.
333
334========== ============================================================================
335key         explanation
336========== ============================================================================
337set                    Specifies a dict with keys and subkeys as described in the
338                       :ref:`Refinement_parameters_fmt` section. Items listed here
339                       will be set to be refined.
340clear                  Specifies a dict, as above for set, except that parameters are
341                       cleared and thus will not be refined.
342once                   Specifies a dict as above for set, except that parameters are
343                       set for the next cycle of refinement and are cleared once the
344                       refinement step is completed.
345skip                   Normally, once parameters are processed with a set/clear/once
346                       action(s), a refinement is started. If skip is defined as True
347                       (or any other value) the refinement step is not performed.
348output                 If a file name is specified for output is will be used to save
349                       the current refinement.
350histograms             Should contain a list of histogram(s) to be used for the
351                       set/clear/once action(s) on :ref:`Histogram_parameters_table` or
352                       :ref:`HAP_parameters_table`. Note that this will be
353                       ignored for :ref:`Phase_parameters_table`. Histograms may be
354                       specified as a list of strings [('PWDR ...'),...], indices
355                       [0,1,2] or as list of objects [hist1, hist2].
356phases                 Should contain a list of phase(s) to be used for the
357                       set/clear/once action(s) on :ref:`Phase_parameters_table` or
358                       :ref:`HAP_parameters_table`. Note that this will be
359                       ignored for :ref:`Histogram_parameters_table`.
360                       Phases may be specified as a list of strings
361                       [('Phase name'),...], indices [0,1,2] or as list of objects
362                       [phase0, phase2].
363call                   Specifies a function to call after a refinement is completed.
364                       The value supplied can be the object (typically a function)
365                       that will be called or a string that will evaluate (in the
366                       namespace inside :meth:`G2Project.iter_refinements` where
367                       ``self`` references the project.)
368                       Nothing is called if this is not specified.
369callargs               Provides a list of arguments that will be passed to the function
370                       in call (if any). If call is defined and callargs is not, the
371                       current <tt>G2Project</tt> is passed as a single argument.
372========== ============================================================================
373
374An example that performs a series of refinement steps follows:
375
376.. code-block::  python
377
378    reflist = [
379            {"set": { "Limits": { "low": 0.7 },
380                      "Background": { "no. coeffs": 3,
381                                      "refine": True }}},
382            {"set": { "LeBail": True,
383                      "Cell": True }},
384            {"set": { "Sample Parameters": ["DisplaceX"]}},
385            {"set": { "Instrument Parameters": ["U", "V", "W", "X", "Y"]}},
386            {"set": { "Mustrain": { "type": "uniaxial",
387                                    "refine": "equatorial",
388                                    "direction": [0, 0, 1]}}},
389            {"set": { "Mustrain": { "type": "uniaxial",
390                                    "refine": "axial"}}},
391            {"clear": { "LeBail": True},
392             "set": { "Atoms": { "Mn": "X" }}},
393            {"set": { "Atoms": { "O1": "X", "O2": "X" }}},]
394    my_project.do_refinements(reflist)
395   
396
397In this example, a separate refinement step will be performed for each dict in the list. The keyword
398"skip" can be used to specify a dict that should not include a refinement.
399Note that in the second from last refinement step, parameters are both set and cleared.
400   
401.. _Refinement_parameters_kinds:
402
403----------------------------
404Refinement parameter types
405----------------------------
406
407Note that parameters and refinement flags used in GSAS-II fall into three classes:
408
409    * **Histogram**: There will be a set of these for each dataset loaded into a
410      project file. The parameters available depend on the type of histogram
411      (Bragg-Brentano, Single-Crystal, TOF,...). Typical Histogram parameters
412      include the overall scale factor, background, instrument and sample parameters;
413      see the :ref:`Histogram_parameters_table` table for a list of the histogram
414      parameters where access has been provided.
415     
416    * **Phase**: There will be a set of these for each phase loaded into a
417      project file. While some parameters are found in all types of phases,
418      others are only found in certain types (modulated, magnetic, protein...).
419      Typical phase parameters include unit cell lengths and atomic positions; see the
420      :ref:`Phase_parameters_table` table for a list of the phase     
421      parameters where access has been provided.
422     
423    * **Histogram-and-phase** (HAP): There is a set of these for every histogram
424      that is associated with each phase, so that if there are ``N`` phases and ``M``
425      histograms, there can be ``N*M`` total sets of "HAP" parameters sets (fewer if all
426      histograms are not linked to all phases.) Typical HAP parameters include the
427      phase fractions, sample microstrain and crystallite size broadening terms,
428      hydrostatic strain perturbations of the unit cell and preferred orientation
429      values.
430      See the :ref:`HAP_parameters_table` table for the HAP parameters where access has
431      been provided.
432
433.. _Refinement_parameters_fmt:
434
435=================================
436Specifying Refinement Parameters
437=================================
438
439Refinement parameter values and flags to turn refinement on and off are specified within dictionaries,
440where the details of these dicts are organized depends on the
441type of parameter (see :ref:`Refinement_parameters_kinds`), with a different set
442of keys (as described below) for each of the three types of parameters.
443
444.. _Histogram_parameters_table:
445
446--------------------
447Histogram parameters
448--------------------
449
450This table describes the dictionaries supplied to :func:`G2PwdrData.set_refinements`
451and :func:`G2PwdrData.clear_refinements`. As an example,
452
453.. code-block::  python
454
455   hist.set_refinements({"Background": {"no.coeffs": 3, "refine": True},
456                         "Sample Parameters": ["Scale"],
457                         "Limits": [10000, 40000]})
458
459With :meth:`G2Project.do_refinements`, these parameters should be placed inside a dict with a key
460``set``, ``clear``, or ``once``. Values will be set for all histograms, unless the ``histograms``
461key is used to define specific histograms. As an example:
462
463.. code-block::  python
464
465  gsas_proj.do_refinements([
466      {'set': {
467          'Background': {'no.coeffs': 3, 'refine': True},
468          'Sample Parameters': ['Scale'],
469          'Limits': [10000, 40000]},
470      'histograms': [1,2]}
471                            ])
472
473Note that below in the Instrument Parameters section,
474related profile parameters (such as U and V) are grouped together but
475separated by commas to save space in the table.
476
477.. tabularcolumns:: |l|l|p{3.5in}|
478
479===================== ====================  =================================================
480key                   subkey                explanation
481===================== ====================  =================================================
482Limits                                      The range of 2-theta (degrees) or TOF (in
483                                            microsec) range of values to use. Can
484                                            be either a dictionary of 'low' and/or 'high',
485                                            or a list of 2 items [low, high]
486\                     low                   Sets the low limit
487\                     high                  Sets the high limit
488
489Sample Parameters                           Should be provided as a **list** of subkeys
490                                            to set or clear, e.g. ['DisplaceX', 'Scale']
491\                     Absorption
492\                     Contrast
493\                     DisplaceX             Sample displacement along the X direction
494\                     DisplaceY             Sample displacement along the Y direction
495\                     Scale                 Histogram Scale factor
496
497Background                                  Sample background. Value will be a dict or
498                                            a boolean. If True or False, the refine
499                                            parameter for background is set to that.
500                                            Note that background peaks are not handled
501                                            via this; see
502                                            :meth:`G2PwdrData.ref_back_peak` instead.
503                                            When value is a dict,
504                                            supply any of the following keys:
505\                     type                  The background model, e.g. 'chebyschev'
506\                     refine                The value of the refine flag, boolean
507\                     no. coeffs            Number of coefficients to use, integer
508\                     coeffs                List of floats, literal values for background
509\                     FixedPoints           List of (2-theta, intensity) values for fixed points
510\                     fit fixed points      If True, triggers a fit to the fixed points to
511                                            be calculated. It is calculated when this key is
512                                            detected, regardless of calls to refine.
513                      peaks                 Specifies a set of flags for refining
514                                            background peaks as a nested list. There may
515                                            be an item for each defined background peak
516                                            (or fewer) and each item is a list with the flag
517                                            values for pos,int,sig & gam (fewer than 4 values
518                                            are allowed).
519
520Instrument Parameters                       As in Sample Paramters, provide as a **list** of
521                                            subkeys to
522                                            set or clear, e.g. ['X', 'Y', 'Zero', 'SH/L']
523\                     U, V, W               Gaussian peak profile terms
524\                     X, Y, Z               Lorentzian peak profile terms
525\                     alpha, beta-0,        TOF profile terms
526                      beta-1, beta-q,
527\                     sig-0, sig-1,         TOF profile terms
528                      sig-2, sig-q
529\                     difA, difB, difC      TOF Calibration constants
530\                     Zero                  Zero shift
531\                     SH/L                  Finger-Cox-Jephcoat low-angle peak asymmetry
532\                     Polariz.              Polarization parameter
533\                     Lam                   Lambda, the incident wavelength
534===================== ====================  =================================================
535
536.. _Phase_parameters_table:
537
538----------------
539Phase parameters
540----------------
541
542This table describes the dictionaries supplied to :func:`G2Phase.set_refinements`
543and :func:`G2Phase.clear_refinements`. With :meth:`G2Project.do_refinements`,
544these parameters should be placed inside a dict with a key
545``set``, ``clear``, or ``once``. Values will be set for all phases, unless the ``phases``
546key is used to define specific phase(s).
547
548
549.. tabularcolumns:: |l|p{4.5in}|
550
551======= ==========================================================
552key                   explanation
553======= ==========================================================
554Cell                  Whether or not to refine the unit cell.
555Atoms                 Dictionary of atoms and refinement flags.
556                      Each key should be an atom label, e.g.
557                      'O3', 'Mn5', and each value should be
558                      a string defining what values to refine.
559                      Values can be any combination of 'F'
560                      for fractional occupancy, 'X' for position,
561                      and 'U' for Debye-Waller factor
562LeBail                Enables LeBail intensity extraction.
563======= ==========================================================
564
565
566.. _HAP_parameters_table:
567
568
569Histogram-and-phase parameters
570------------------------------
571
572This table describes the dictionaries supplied to :func:`G2Phase.set_HAP_refinements`
573and :func:`G2Phase.clear_HAP_refinements`. When supplied to
574:meth:`G2Project.do_refinements`, these parameters should be placed inside a dict with a key
575``set``, ``clear``, or ``once``. Values will be set for all histograms used in each phase,
576unless the ``histograms`` and ``phases`` keys are used to define specific phases and histograms.
577
578.. tabularcolumns:: |l|l|p{3.5in}|
579
580=============  ==========  ============================================================
581key             subkey                 explanation
582=============  ==========  ============================================================
583Babinet                                Should be a **list** of the following
584                                       subkeys. If not, assumes both
585                                       BabA and BabU
586\               BabA
587\               BabU
588Extinction                             Boolean, True to refine.
589HStrain                                Boolean or list/tuple, True to refine all
590                                       appropriate D\ :sub:`ij` terms or False
591                                       to not refine any. If a list/tuple, will
592                                       be a set of True & False values for each
593                                       D\ :sub:`ij` term; number of items must
594                                       match number of terms.
595Mustrain
596\               type                   Mustrain model. One of 'isotropic',
597                                       'uniaxial', or 'generalized'. **Should always
598                                       be included when Mustrain is used.**
599\              direction               For uniaxial only. A list of three
600                                       integers,
601                                       the [hkl] direction of the axis.
602\               refine                 Usually boolean, set to True to refine.
603                                       or False to clear.
604                                       For uniaxial model, can specify a value
605                                       of 'axial' or 'equatorial' to set that flag
606                                       to True or a single
607                                       boolean sets both axial and equatorial.
608Size                                   
609\               type                   Size broadening model. One of 'isotropic',
610                                       'uniaxial', or 'ellipsoid'. **Should always
611                                       be specified when Size is used.**
612\              direction               For uniaxial only. A list of three
613                                       integers,
614                                       the [hkl] direction of the axis.
615\               refine                 Boolean, True to refine.
616\               value                  float, size value in microns
617Pref.Ori.                              Boolean, True to refine
618Show                                   Boolean, True to refine
619Use                                    Boolean, True to refine
620Scale                                  Phase fraction; Boolean, True to refine
621=============  ==========  ============================================================
622
623------------------------
624Histogram/Phase objects
625------------------------
626Each phase and powder histogram in a :class:`G2Project` object has an associated
627object. Parameters within each individual object can be turned on and off by calling
628:meth:`G2PwdrData.set_refinements` or :meth:`G2PwdrData.clear_refinements`
629for histogram parameters;
630:meth:`G2Phase.set_refinements` or :meth:`G2Phase.clear_refinements`
631for phase parameters; and :meth:`G2Phase.set_HAP_refinements` or
632:meth:`G2Phase.clear_HAP_refinements`. As an example, if some_histogram is a histogram object (of type :class:`G2PwdrData`), use this to set parameters in that histogram:
633
634.. code-block::  python
635
636    params = { 'Limits': [0.8, 12.0],
637               'Sample Parameters': ['Absorption', 'Contrast', 'DisplaceX'],
638               'Background': {'type': 'chebyschev', 'refine': True}}
639    some_histogram.set_refinements(params)
640
641Likewise to turn refinement flags on, use code such as this:
642
643.. code-block::  python
644
645    params = { 'Instrument Parameters': ['U', 'V', 'W']}
646    some_histogram.set_refinements(params)
647
648and to turn these refinement flags, off use this (Note that the
649``.clear_refinements()`` methods will usually will turn off refinement even
650if a refinement parameter is set in the dict to True.):
651
652.. code-block::  python
653
654    params = { 'Instrument Parameters': ['U', 'V', 'W']}
655    some_histogram.clear_refinements(params)
656
657For phase parameters, use code such as this:
658   
659.. code-block::  python
660
661    params = { 'LeBail': True, 'Cell': True,
662               'Atoms': { 'Mn1': 'X',
663                          'O3': 'XU',
664                          'V4': 'FXU'}}
665    some_histogram.set_refinements(params)
666
667
668and here is an example for HAP parameters:
669
670.. code-block::  python
671
672    params = { 'Babinet': 'BabA',
673               'Extinction': True,
674               'Mustrain': { 'type': 'uniaxial',
675                             'direction': [0, 0, 1],
676                             'refine': True}}
677    some_phase.set_HAP_refinements(params)
678
679Note that the parameters must match the object type and method (phase vs. histogram vs. HAP).
680
681.. _CommandlineInterface:
682
683=======================================
684GSASIIscriptable Command-line Interface
685=======================================
686
687The routines described above are intended to be called from a Python script, but an
688alternate way to access some of the same functionality is to
689invoke the ``GSASIIscriptable.py`` script from
690the command line usually from within a shell script or batch file. This
691will usually be done with a command such as::
692
693       python <path/>GSASIIscriptable.py <subcommand> <file.gpx> <options>
694
695The following subcommands are defined:
696
697        * create, see :func:`create`
698        * add, see :func:`add`
699        * dump, see :func:`dump`
700        * refine, see :func:`refine`
701        * export, :func:`export`
702        * browse, see :func:`IPyBrowse`
703       
704Run::
705
706   python GSASIIscriptable.py --help
707
708to show the available subcommands, and inspect each subcommand with
709`python GSASIIscriptable.py <subcommand> --help` or see the documentation for each of the above routines.
710
711.. _JsonFormat:
712
713-------------------------
714Parameters in JSON files
715-------------------------
716
717The refine command requires two inputs: an existing GSAS-II project (.gpx) file and
718a JSON format file
719(see `Introducing JSON <http://json.org/>`_) that contains a single dict.
720This dict may have two keys:
721
722refinements:
723  This defines the a set of refinement steps in a JSON representation of a
724  :ref:`Refinement_recipe` list.
725
726code:
727  This optionally defines Python code that will be executed after the project is loaded,
728  but before the refinement is started. This can be used to execute Python code to change
729  parameters that are not accessible via a :ref:`Refinement_recipe` dict (note that the
730  project object is accessed with variable ``proj``) or to define code that will be called
731  later (see key ``call`` in the :ref:`Refinement_recipe` section.)
732   
733JSON website: `Introducing JSON <http://json.org/>`_.
734
735.. _API:
736
737============================================================
738API: Complete Documentation
739============================================================
740
741The large number of classes and modules in this module are described below.
742A script will have one or more G2Project objects using :class:`G2Project` and then
743perform actions such as adding a histogram (method :meth:`G2Project.add_powder_histogram`),
744adding a phase (method :meth:`G2Project.add_phase`),
745or setting parameters and performing a refinement
746(method :meth:`G2Project.do_refinements`).
747
748To change settings within histograms, images and phases, one usually needs to use
749methods inside :class:`G2PwdrData`, :class:`G2Image` or :class:`G2Phase`.
750"""
751
752#============================================================================
753# adding a new object type
754# 1) add a new object class (e.g. G2PDF)
755# 2) add the wrapper into G2Project (e.g. _pdfs, pdf, pdfs)
756# 3) add a new method to add the object into a project (G2Project.add_PDF)
757# 4) add to documentation in section :class:`G2Project`
758# 5) add a new documentation section for the new class
759#============================================================================
760
761from __future__ import division, print_function
762import argparse
763import os.path as ospath
764import datetime as dt
765import sys
766import platform
767if '2' in platform.python_version_tuple()[0]:
768    import cPickle
769    strtypes = (str,unicode)
770else:
771    import pickle as cPickle
772    strtypes = (str,bytes)
773import imp
774import copy
775import os
776import random as ran
777
778import numpy.ma as ma
779import scipy.interpolate as si
780import numpy as np
781import scipy as sp
782
783import GSASIIpath
784GSASIIpath.SetBinaryPath(True)  # for now, this is needed before some of these modules can be imported
785import GSASIIobj as G2obj
786import GSASIIpwd as G2pwd
787import GSASIIstrMain as G2strMain
788#import GSASIIIO as G2IO
789import GSASIIstrIO as G2strIO
790import GSASIIspc as G2spc
791import GSASIIElem as G2elem
792import GSASIIfiles as G2fil
793import GSASIIimage as G2img
794
795# Delay imports to not slow down small scripts that don't need them
796Readers = {'Pwdr':[], 'Phase':[], 'Image':[]}
797'''Readers by reader type'''
798exportersByExtension = {}
799'''Specifies the list of extensions that are supported for Powder data export'''
800npsind = lambda x: np.sin(x*np.pi/180.)
801
802def LoadG2fil():
803    """Setup GSAS-II importers. Delay importing this module, it is slow"""
804    if len(Readers['Pwdr']) > 0: return
805    # initialize imports
806    Readers['Pwdr'] = G2fil.LoadImportRoutines("pwd", "Powder_Data")
807    Readers['Phase'] = G2fil.LoadImportRoutines("phase", "Phase")
808    Readers['Image'] = G2fil.LoadImportRoutines("img", "Image")
809
810    # initialize exports
811    for obj in G2fil.LoadExportRoutines(None):
812        try:
813            obj.Writer
814        except AttributeError:
815            continue
816        for typ in obj.exporttype:
817            if typ not in exportersByExtension:
818                exportersByExtension[typ] = {obj.extension:obj}
819            else:
820                exportersByExtension[typ][obj.extension] = obj
821
822def LoadDictFromProjFile(ProjFile):
823    '''Read a GSAS-II project file and load items to dictionary
824   
825    :param str ProjFile: GSAS-II project (name.gpx) full file name
826    :returns: Project,nameList, where
827
828      * Project (dict) is a representation of gpx file following the GSAS-II tree structure
829        for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict
830        data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)}
831      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
832
833    Example for fap.gpx::
834
835      Project = {                 #NB:dict order is not tree order
836        'Phases':{'data':None,'fap':{phase dict}},
837        'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc},
838        'Rigid bodies':{'data': {rigid body dict}},
839        'Covariance':{'data':{covariance data dict}},
840        'Controls':{'data':{controls data dict}},
841        'Notebook':{'data':[notebook list]},
842        'Restraints':{'data':{restraint data dict}},
843        'Constraints':{'data':{constraint data dict}}]
844        }
845      nameList = [                #NB: reproduces tree order
846        ['Notebook',],
847        ['Controls',],
848        ['Covariance',],
849        ['Constraints',],
850        ['Restraints',],
851        ['Rigid bodies',],
852        ['PWDR FAP.XRA Bank 1',
853             'Comments',
854             'Limits',
855             'Background',
856             'Instrument Parameters',
857             'Sample Parameters',
858             'Peak List',
859             'Index Peak List',
860             'Unit Cells List',
861             'Reflection Lists'],
862        ['Phases', 'fap']
863        ]
864    '''
865    # Let IOError be thrown if file does not exist
866    if not ospath.exists(ProjFile):
867        print ('\n*** Error attempt to open project file that does not exist: \n    {}'.
868                   format(ProjFile))
869        raise IOError('GPX file {} does not exist'.format(ProjFile))
870    try:
871        Project, nameList = G2strIO.GetFullGPX(ProjFile)
872    except Exception as msg:
873        raise IOError(msg)
874    return Project,nameList
875
876def SaveDictToProjFile(Project,nameList,ProjFile):
877    '''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile
878
879    :param dict Project: representation of gpx file following the GSAS-II
880        tree structure as described for LoadDictFromProjFile
881    :param list nameList: names of main tree entries & subentries used to reconstruct project file
882    :param str ProjFile: full file name for output project.gpx file (including extension)
883    '''
884    file = open(ProjFile,'wb')
885    try:
886        for name in nameList:
887            data = []
888            item = Project[name[0]]
889            data.append([name[0],item['data']])
890            for item2 in name[1:]:
891                data.append([item2,item[item2]])
892            cPickle.dump(data,file,1)
893    finally:
894        file.close()
895    print('gpx file saved as %s'%ProjFile)
896
897# def ImportPowder(reader,filename):
898#     '''Use a reader to import a powder diffraction data file
899
900#     :param str reader: a scriptable reader
901#     :param str filename: full name of powder data file; can be "multi-Bank" data
902
903#     :returns: list rdlist: list of reader objects containing powder data, one for each
904#         "Bank" of data encountered in file. Items in reader object of interest are:
905
906#           * rd.comments: list of str: comments found on powder file
907#           * rd.dnames: list of str: data nammes suitable for use in GSASII data tree NB: duplicated in all rd entries in rdlist
908#           * rd.powderdata: list of numpy arrays: pos,int,wt,zeros,zeros,zeros as needed for a PWDR entry in  GSASII data tree.
909#     '''
910#     rdfile,rdpath,descr = imp.find_module(reader)
911#     rdclass = imp.load_module(reader,rdfile,rdpath,descr)
912#     rd = rdclass.GSAS_ReaderClass()
913#     if not rd.scriptable:
914#         print(u'**** ERROR: '+reader+u' is not a scriptable reader')
915#         return None
916#     rdlist = []
917#     if rd.ContentsValidator(filename):
918#         repeat = True
919#         rdbuffer = {} # create temporary storage for file reader
920#         block = 0
921#         while repeat: # loop if the reader asks for another pass on the file
922#             block += 1
923#             repeat = False
924#             rd.objname = ospath.basename(filename)
925#             flag = rd.Reader(filename,None,buffer=rdbuffer,blocknum=block,)
926#             if flag:
927#                 rdlist.append(copy.deepcopy(rd)) # save the result before it is written over
928#                 if rd.repeat:
929#                     repeat = True
930#         return rdlist
931#     print(rd.errors)
932#     return None
933
934def SetDefaultDData(dType,histoName,NShkl=0,NDij=0):
935    '''Create an initial Histogram dictionary
936
937    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
938    '''
939    if dType in ['SXC','SNC']:
940        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
941            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
942            'Extinction':['Lorentzian','None', {'Tbar':0.1,'Cos2TM':0.955,
943            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}],
944            'Flack':[0.0,False]}
945    elif dType == 'SNT':
946        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
947            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
948            'Extinction':['Lorentzian','None', {
949            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}]}
950    elif 'P' in dType:
951        return {'Histogram':histoName,'Show':False,'Scale':[1.0,False],
952            'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{},[],0.1],
953            'Size':['isotropic',[1.,1.,1.],[False,False,False],[0,0,1],
954                [1.,1.,1.,0.,0.,0.],6*[False,]],
955            'Mustrain':['isotropic',[1000.0,1000.0,1.0],[False,False,False],[0,0,1],
956                NShkl*[0.01,],NShkl*[False,]],
957            'HStrain':[NDij*[0.0,],NDij*[False,]],
958            'Extinction':[0.0,False],'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}}
959
960
961def PreSetup(data):
962    '''Create part of an initial (empty) phase dictionary
963
964    from GSASIIphsGUI.py, near end of UpdatePhaseData
965
966    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
967    '''
968    if 'RBModels' not in data:
969        data['RBModels'] = {}
970    if 'MCSA' not in data:
971        data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[.8,1.2],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
972    if 'dict' in str(type(data['MCSA']['Results'])):
973        data['MCSA']['Results'] = []
974    if 'Modulated' not in data['General']:
975        data['General']['Modulated'] = False
976#    if 'modulated' in data['General']['Type']:
977#        data['General']['Modulated'] = True
978#        data['General']['Type'] = 'nuclear'
979
980
981def SetupGeneral(data, dirname):
982    """Helps initialize phase data.
983
984    From GSASIIphsGui.py, function of the same name. Minor changes for imports etc.
985
986    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
987    """
988    mapDefault = {'MapType':'','RefList':'','Resolution':0.5,'Show bonds':True,
989                'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False}
990    generalData = data['General']
991    atomData = data['Atoms']
992    generalData['AtomTypes'] = []
993    generalData['Isotopes'] = {}
994
995    if 'Isotope' not in generalData:
996        generalData['Isotope'] = {}
997    if 'Data plot type' not in generalData:
998        generalData['Data plot type'] = 'Mustrain'
999    if 'POhkl' not in generalData:
1000        generalData['POhkl'] = [0,0,1]
1001    if 'Map' not in generalData:
1002        generalData['Map'] = mapDefault.copy()
1003    if 'Flip' not in generalData:
1004        generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None',
1005            'k-factor':0.1,'k-Max':20.,}
1006    if 'testHKL' not in generalData['Flip']:
1007        generalData['Flip']['testHKL'] = [[0,0,2],[2,0,0],[1,1,1],[0,2,0],[1,2,3]]
1008    if 'doPawley' not in generalData:
1009        generalData['doPawley'] = False     #ToDo: change to ''
1010    if 'Pawley dmin' not in generalData:
1011        generalData['Pawley dmin'] = 1.0
1012    if 'Pawley neg wt' not in generalData:
1013        generalData['Pawley neg wt'] = 0.0
1014    if 'Algolrithm' in generalData.get('MCSA controls',{}) or \
1015        'MCSA controls' not in generalData:
1016        generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50],
1017        'dmin':2.0,'Algorithm':'log','Jump coeff':[0.95,0.5],'boltzmann':1.0,
1018        'fast parms':[1.0,1.0,1.0],'log slope':0.9,'Cycles':1,'Results':[],'newDmin':True}
1019    if 'AtomPtrs' not in generalData:
1020        generalData['AtomPtrs'] = [3,1,7,9]
1021        if generalData['Type'] == 'macromolecular':
1022            generalData['AtomPtrs'] = [6,4,10,12]
1023        elif generalData['Type'] == 'magnetic':
1024            generalData['AtomPtrs'] = [3,1,10,12]
1025    if generalData['Modulated']:
1026        generalData['Type'] = 'nuclear'
1027        if 'Super' not in generalData:
1028            generalData['Super'] = 1
1029            generalData['SuperVec'] = [[0,0,.1],False,4]
1030            generalData['SSGData'] = {}
1031        if '4DmapData' not in generalData:
1032            generalData['4DmapData'] = mapDefault.copy()
1033            generalData['4DmapData'].update({'MapType':'Fobs'})
1034    if 'Modulated' not in generalData:
1035        generalData['Modulated'] = False
1036    if 'HydIds' not in generalData:
1037        generalData['HydIds'] = {}
1038    cx,ct,cs,cia = generalData['AtomPtrs']
1039    generalData['NoAtoms'] = {}
1040    generalData['BondRadii'] = []
1041    generalData['AngleRadii'] = []
1042    generalData['vdWRadii'] = []
1043    generalData['AtomMass'] = []
1044    generalData['Color'] = []
1045    if generalData['Type'] == 'magnetic':
1046        generalData['MagDmin'] = generalData.get('MagDmin',1.0)
1047        landeg = generalData.get('Lande g',[])
1048    generalData['Mydir'] = dirname
1049    badList = {}
1050    for iat,atom in enumerate(atomData):
1051        atom[ct] = atom[ct].lower().capitalize()              #force to standard form
1052        if generalData['AtomTypes'].count(atom[ct]):
1053            generalData['NoAtoms'][atom[ct]] += atom[cx+3]*float(atom[cs+1])
1054        elif atom[ct] != 'UNK':
1055            Info = G2elem.GetAtomInfo(atom[ct])
1056            if not Info:
1057                if atom[ct] not in badList:
1058                    badList[atom[ct]] = 0
1059                badList[atom[ct]] += 1
1060                atom[ct] = 'UNK'
1061                continue
1062            atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo
1063            generalData['AtomTypes'].append(atom[ct])
1064            generalData['Z'] = Info['Z']
1065            generalData['Isotopes'][atom[ct]] = Info['Isotopes']
1066            generalData['BondRadii'].append(Info['Drad'])
1067            generalData['AngleRadii'].append(Info['Arad'])
1068            generalData['vdWRadii'].append(Info['Vdrad'])
1069            if atom[ct] in generalData['Isotope']:
1070                if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]:
1071                    isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
1072                    generalData['Isotope'][atom[ct]] = isotope
1073                generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][atom[ct]]]['Mass'])
1074            else:
1075                generalData['Isotope'][atom[ct]] = 'Nat. Abund.'
1076                if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]:
1077                    isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
1078                    generalData['Isotope'][atom[ct]] = isotope
1079                generalData['AtomMass'].append(Info['Mass'])
1080            generalData['NoAtoms'][atom[ct]] = atom[cx+3]*float(atom[cs+1])
1081            generalData['Color'].append(Info['Color'])
1082            if generalData['Type'] == 'magnetic':
1083                if len(landeg) < len(generalData['AtomTypes']):
1084                    landeg.append(2.0)
1085    if generalData['Type'] == 'magnetic':
1086        generalData['Lande g'] = landeg[:len(generalData['AtomTypes'])]
1087
1088    if badList:
1089        msg = 'Warning: element symbol(s) not found:'
1090        for key in badList:
1091            msg += '\n\t' + key
1092            if badList[key] > 1:
1093                msg += ' (' + str(badList[key]) + ' times)'
1094        raise G2ScriptException("Phase error:\n" + msg)
1095        # wx.MessageBox(msg,caption='Element symbol error')
1096    F000X = 0.
1097    F000N = 0.
1098    for i,elem in enumerate(generalData['AtomTypes']):
1099        F000X += generalData['NoAtoms'][elem]*generalData['Z']
1100        isotope = generalData['Isotope'][elem]
1101        F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope]['SL'][0]
1102    generalData['F000X'] = F000X
1103    generalData['F000N'] = F000N
1104    import GSASIImath as G2mth
1105    generalData['Mass'] = G2mth.getMass(generalData)
1106
1107
1108def make_empty_project(author=None, filename=None):
1109    """Creates an dictionary in the style of GSASIIscriptable, for an empty
1110    project.
1111
1112    If no author name or filename is supplied, 'no name' and
1113    <current dir>/test_output.gpx are used , respectively.
1114
1115    Returns: project dictionary, name list
1116
1117    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1118    """
1119    if not filename:
1120        filename = 'test_output.gpx'
1121    filename = os.path.abspath(filename)
1122    gsasii_version = str(GSASIIpath.GetVersionNumber())
1123    LoadG2fil()
1124    try:
1125        import matplotlib as mpl
1126        python_library_versions = G2fil.get_python_versions([mpl, np, sp])
1127    except ImportError:
1128        python_library_versions = G2fil.get_python_versions([np, sp])
1129
1130    controls_data = dict(G2obj.DefaultControls)
1131    controls_data['LastSavedAs'] = filename
1132    controls_data['LastSavedUsing'] = gsasii_version
1133    controls_data['PythonVersions'] = python_library_versions
1134    if author:
1135        controls_data['Author'] = author
1136
1137    output = {'Constraints': {'data': {'HAP': [], 'Hist': [], 'Phase': [],
1138                                       'Global': []}},
1139              'Controls': {'data': controls_data},
1140              u'Covariance': {'data': {}},
1141              u'Notebook': {'data': ['']},
1142              u'Restraints': {'data': {}},
1143              u'Rigid bodies': {'data': {'RBIds': {'Residue': [], 'Vector': []},
1144                                'Residue': {'AtInfo': {}},
1145                                'Vector':  {'AtInfo': {}}}}}
1146
1147    names = [[u'Notebook'], [u'Controls'], [u'Covariance'],
1148             [u'Constraints'], [u'Restraints'], [u'Rigid bodies']]
1149
1150    return output, names
1151
1152
1153class G2ImportException(Exception):
1154    pass
1155
1156class G2ScriptException(Exception):
1157    pass
1158
1159def import_generic(filename, readerlist, fmthint=None, bank=None):
1160    """Attempt to import a filename, using a list of reader objects.
1161
1162    Returns the first reader object which worked."""
1163    # Translated from OnImportGeneric method in GSASII.py
1164    primaryReaders, secondaryReaders = [], []
1165    for reader in readerlist:
1166        if fmthint is not None and fmthint not in reader.formatName: continue
1167        flag = reader.ExtensionValidator(filename)
1168        if flag is None:
1169            secondaryReaders.append(reader)
1170        elif flag:
1171            primaryReaders.append(reader)
1172    if not secondaryReaders and not primaryReaders:
1173        raise G2ImportException("Could not read file: ", filename)
1174
1175    with open(filename, 'Ur') as fp:
1176        rd_list = []
1177
1178        for rd in primaryReaders + secondaryReaders:
1179            # Initialize reader
1180            rd.selections = []
1181            if bank is None:
1182                rd.selections = []
1183            else:
1184                rd.selections = [bank-1]
1185            rd.dnames = []
1186            rd.ReInitialize()
1187            # Rewind file
1188            rd.errors = ""
1189            if not rd.ContentsValidator(filename):
1190                # Report error
1191                print("Warning: File {} has a validation error, continuing".format(filename))
1192            if len(rd.selections) > 1:
1193                raise G2ImportException("File {} has {} banks. Specify which bank to read with databank param."
1194                                .format(filename,len(rd.selections)))
1195
1196            block = 0
1197            rdbuffer = {}
1198            repeat = True
1199            while repeat:
1200                repeat = False
1201                block += 1
1202                rd.objname = os.path.basename(filename)
1203                try:
1204                    flag = rd.Reader(filename,buffer=rdbuffer, blocknum=block)
1205                except:
1206                    flag = False
1207                if flag:
1208                    # Omitting image loading special cases
1209                    rd.readfilename = filename
1210                    rd_list.append(copy.deepcopy(rd))
1211                    repeat = rd.repeat
1212                else:
1213                    if GSASIIpath.GetConfigValue('debug'): print("DBG_{} Reader failed to read {}".format(rd.formatName,filename))
1214            if rd_list:
1215                if rd.warnings:
1216                    print("Read warning by", rd.formatName, "reader:",
1217                          rd.warnings, file=sys.stderr)
1218                elif bank is None:
1219                    print("{} read by Reader {}"
1220                              .format(filename,rd.formatName))
1221                else:
1222                    print("{} block # {} read by Reader {}"
1223                              .format(filename,bank,rd.formatName))
1224                return rd_list
1225    raise G2ImportException("No reader could read file: " + filename)
1226
1227
1228def load_iprms(instfile, reader, bank=None):
1229    """Loads instrument parameters from a file, and edits the
1230    given reader.
1231
1232    Returns a 2-tuple of (Iparm1, Iparm2) parameters
1233    """
1234    LoadG2fil()
1235    ext = os.path.splitext(instfile)[1]
1236
1237    if ext.lower() == '.instprm':
1238        # New GSAS File, load appropriate bank
1239        with open(instfile) as f:
1240            lines = f.readlines()
1241        if bank is None: 
1242            bank = reader.powderentry[2] 
1243        numbanks = reader.numbanks
1244        iparms = G2fil.ReadPowderInstprm(lines, bank, numbanks, reader)
1245        reader.instfile = instfile
1246        reader.instmsg = '{} (G2 fmt) bank {}'.format(instfile,bank)
1247        return iparms
1248    elif ext.lower() not in ('.prm', '.inst', '.ins'):
1249        raise ValueError('Expected .prm file, found: ', instfile)
1250
1251    # It's an old GSAS file, load appropriately
1252    Iparm = {}
1253    with open(instfile, 'Ur') as fp:
1254        for line in fp:
1255            if '#' in line:
1256                continue
1257            Iparm[line[:12]] = line[12:-1]
1258    ibanks = int(Iparm.get('INS   BANK  ', '1').strip())
1259    if bank is not None:
1260        # pull out requested bank # bank from the data, and change the bank to 1
1261        Iparm,IparmC = {},Iparm
1262        for key in IparmC:
1263            if 'INS' not in key[:3]: continue   #skip around rubbish lines in some old iparm
1264            if key[4:6] == "  ":
1265                Iparm[key] = IparmC[key]
1266            elif int(key[4:6].strip()) == bank:
1267                Iparm[key[:4]+' 1'+key[6:]] = IparmC[key]           
1268        reader.instbank = bank
1269    elif ibanks == 1:
1270        reader.instbank = 1
1271    else: 
1272        raise G2ImportException("Instrument parameter file has {} banks, select one with instbank param."
1273                                    .format(ibanks))
1274    reader.powderentry[2] = 1
1275    reader.instfile = instfile
1276    reader.instmsg = '{} bank {}'.format(instfile,reader.instbank)
1277    return G2fil.SetPowderInstParms(Iparm, reader)
1278
1279def load_pwd_from_reader(reader, instprm, existingnames=[],bank=None):
1280    """Loads powder data from a reader object, and assembles it into a GSASII data tree.
1281
1282    :returns: (name, tree) - 2-tuple of the histogram name (str), and data
1283
1284    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1285    """
1286    HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_')
1287    HistName = G2obj.MakeUniqueLabel(HistName, existingnames)
1288
1289    try:
1290        Iparm1, Iparm2 = instprm
1291    except ValueError:
1292        Iparm1, Iparm2 = load_iprms(instprm, reader, bank=bank)
1293        print('Instrument parameters read:',reader.instmsg)
1294    Ymin = np.min(reader.powderdata[1])
1295    Ymax = np.max(reader.powderdata[1])
1296    valuesdict = {'wtFactor': 1.0,
1297                  'Dummy': False,
1298                  'ranId': ran.randint(0, sys.maxsize),
1299                  'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax,
1300                  'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax,
1301                  'Yminmax': [Ymin, Ymax]}
1302    reader.Sample['ranId'] = valuesdict['ranId']
1303    if 'T' in Iparm1['Type'][0]:
1304        if not reader.clockWd and reader.GSAS:
1305            reader.powderdata[0] *= 100.        #put back the CW centideg correction
1306
1307    # Ending keys:
1308    # [u'Reflection Lists',
1309    #  u'Limits',
1310    #  'data',
1311    #  u'Index Peak List',
1312    #  u'Comments',
1313    #  u'Unit Cells List',
1314    #  u'Sample Parameters',
1315    #  u'Peak List',
1316    #  u'Background',
1317    #  u'Instrument Parameters']
1318    Tmin = np.min(reader.powderdata[0])
1319    Tmax = np.max(reader.powderdata[0])
1320
1321    default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0],
1322                          {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}]
1323
1324    output_dict = {u'Reflection Lists': {},
1325                   u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]),
1326                   u'data': [valuesdict, reader.powderdata, HistName],
1327                   u'Index Peak List': [[], []],
1328                   u'Comments': reader.comments,
1329                   u'Unit Cells List': [],
1330                   u'Sample Parameters': reader.Sample,
1331                   u'Peak List': {'peaks': [], 'sigDict': {}},
1332                   u'Background': reader.pwdparms.get('Background', default_background),
1333                   u'Instrument Parameters': [Iparm1, Iparm2],
1334                   }
1335
1336    names = [u'Comments',
1337             u'Limits',
1338             u'Background',
1339             u'Instrument Parameters',
1340             u'Sample Parameters',
1341             u'Peak List',
1342             u'Index Peak List',
1343             u'Unit Cells List',
1344             u'Reflection Lists']
1345
1346    # TODO controls?? GSASII.py:1664-7
1347
1348    return HistName, [HistName] + names, output_dict
1349
1350def _deep_copy_into(from_, into):
1351    """Helper function for reloading .gpx file. See G2Project.reload()
1352
1353    :author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1354    """
1355    if isinstance(from_, dict) and isinstance(into, dict):
1356        combined_keys = set(from_.keys()).union(into.keys())
1357        for key in combined_keys:
1358            if key in from_ and key in into:
1359                both_dicts = (isinstance(from_[key], dict)
1360                              and isinstance(into[key], dict))
1361                both_lists = (isinstance(from_[key], list)
1362                              and isinstance(into[key], list))
1363                if both_dicts or both_lists:
1364                    _deep_copy_into(from_[key], into[key])
1365                else:
1366                    into[key] = from_[key]
1367            elif key in from_:
1368                into[key] = from_[key]
1369            else:  # key in into
1370                del into[key]
1371    elif isinstance(from_, list) and isinstance(into, list):
1372        if len(from_) == len(into):
1373            for i in range(len(from_)):
1374                both_dicts = (isinstance(from_[i], dict)
1375                              and isinstance(into[i], dict))
1376                both_lists = (isinstance(from_[i], list)
1377                              and isinstance(into[i], list))
1378                if both_dicts or both_lists:
1379                    _deep_copy_into(from_[i], into[i])
1380                else:
1381                    into[i] = from_[i]
1382        else:
1383            into[:] = from_
1384
1385def _getCorrImage(ImageReaderlist,proj,imageRef):
1386    '''Gets image & applies dark, background & flat background corrections.
1387    based on :func:`GSASIIimgGUI.GetImageZ`. Expected to be for internal
1388    use only.
1389
1390    :param list ImageReaderlist: list of Reader objects for images
1391    :param object ImageReaderlist: list of Reader objects for images
1392    :param imageRef: A reference to the desired image. Either the Image
1393      tree name (str), the image's index (int) or
1394      a image object (:class:`G2Image`)
1395
1396    :return: array sumImg: corrected image for background/dark/flat back
1397    '''
1398    ImgObj = proj.image(imageRef)
1399    Controls = ImgObj.data['Image Controls']
1400    formatName = Controls.get('formatName','')
1401    imagefile = ImgObj.data['data'][1]
1402    ImageTag = None # fix this for multiimage files
1403    sumImg = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
1404    if sumImg is None:
1405        return []
1406    darkImg = False
1407    if 'dark image' in Controls:
1408        darkImg,darkScale = Controls['dark image']
1409        if darkImg:
1410            dImgObj = proj.image(darkImg)
1411            formatName = dImgObj.data['Image Controls'].get('formatName','')
1412            imagefile = dImgObj.data['data'][1]
1413            ImageTag = None # fix this for multiimage files
1414            darkImage = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
1415            if darkImg is None:
1416                raise Exception('Error reading dark image {}'.format(imagefile))
1417            sumImg += np.array(darkImage*darkScale,dtype='int32')
1418    if 'background image' in Controls:
1419        backImg,backScale = Controls['background image']           
1420        if backImg:     #ignores any transmission effect in the background image
1421            bImgObj = proj.image(backImg)
1422            formatName = bImgObj.data['Image Controls'].get('formatName','')
1423            imagefile = bImgObj.data['data'][1]
1424            ImageTag = None # fix this for multiimage files
1425            backImage = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
1426            if backImage is None:
1427                raise Exception('Error reading background image {}'.format(imagefile))
1428            if darkImg:
1429                backImage += np.array(darkImage*darkScale/backScale,dtype='int32')
1430            else:
1431                sumImg += np.array(backImage*backScale,dtype='int32')
1432    if 'Gain map' in Controls:
1433        gainMap = Controls['Gain map']
1434        if gainMap:
1435            gImgObj = proj.image(gainMap)
1436            formatName = gImgObj.data['Image Controls'].get('formatName','')
1437            imagefile = gImgObj.data['data'][1]
1438            ImageTag = None # fix this for multiimage files
1439            GMimage = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
1440            if GMimage is None:
1441                raise Exception('Error reading Gain map image {}'.format(imagefile))
1442            sumImg = sumImg*GMimage/1000
1443    sumImg -= int(Controls.get('Flat Bkg',0))
1444    Imax = np.max(sumImg)
1445    Controls['range'] = [(0,Imax),[0,Imax]]
1446    return np.asarray(sumImg,dtype='int32')
1447
1448class G2ObjectWrapper(object):
1449    """Base class for all GSAS-II object wrappers.
1450
1451    The underlying GSAS-II format can be accessed as `wrapper.data`. A number
1452    of overrides are implemented so that the wrapper behaves like a dictionary.
1453
1454    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1455    """
1456    def __init__(self, datadict):
1457        self.data = datadict
1458
1459    def __getitem__(self, key):
1460        return self.data[key]
1461
1462    def __setitem__(self, key, value):
1463        self.data[key] = value
1464
1465    def __contains__(self, key):
1466        return key in self.data
1467
1468    def get(self, k, d=None):
1469        return self.data.get(k, d)
1470
1471    def keys(self):
1472        return self.data.keys()
1473
1474    def values(self):
1475        return self.data.values()
1476
1477    def items(self):
1478        return self.data.items()
1479
1480
1481class G2Project(G2ObjectWrapper):   
1482    """Represents an entire GSAS-II project.
1483
1484    :param str gpxfile: Existing .gpx file to be loaded. If nonexistent,
1485            creates an empty project.
1486    :param str author: Author's name (not yet implemented)
1487    :param str newgpx: The filename the project should be saved to in
1488            the future. If both newgpx and gpxfile are present, the project is
1489            loaded from the gpxfile, then when saved will be written to newgpx.
1490    :param str filename: Name to be used to save the project. Has same function as
1491            parameter newgpx (do not use both gpxfile and filename). Use of newgpx
1492            is preferred over filename.
1493
1494    There are two ways to initialize this object:
1495
1496    >>> # Load an existing project file
1497    >>> proj = G2Project('filename.gpx')
1498   
1499    >>> # Create a new project
1500    >>> proj = G2Project(newgpx='new_file.gpx')
1501   
1502    Histograms can be accessed easily.
1503
1504    >>> # By name
1505    >>> hist = proj.histogram('PWDR my-histogram-name')
1506   
1507    >>> # Or by index
1508    >>> hist = proj.histogram(0)
1509    >>> assert hist.id == 0
1510   
1511    >>> # Or by random id
1512    >>> assert hist == proj.histogram(hist.ranId)
1513
1514    Phases can be accessed the same way.
1515
1516    >>> phase = proj.phase('name of phase')
1517
1518    New data can also be loaded via :meth:`~G2Project.add_phase` and
1519    :meth:`~G2Project.add_powder_histogram`.
1520
1521    >>> hist = proj.add_powder_histogram('some_data_file.chi',
1522                                         'instrument_parameters.prm')
1523    >>> phase = proj.add_phase('my_phase.cif', histograms=[hist])
1524
1525    Parameters for Rietveld refinement can be turned on and off as well.
1526    See :meth:`~G2Project.set_refinement`, :meth:`~G2Project.clear_refinements`,
1527    :meth:`~G2Project.iter_refinements`, :meth:`~G2Project.do_refinements`.
1528    """
1529    def __init__(self, gpxfile=None, author=None, filename=None, newgpx=None):
1530        if filename is not None and newgpx is not None:
1531            raise G2ScriptException('Do not use filename and newgpx together')
1532        elif newgpx is not None:
1533            filename = newgpx
1534        if gpxfile is None:
1535            filename = os.path.abspath(os.path.expanduser(filename))
1536            self.filename = filename
1537            self.data, self.names = make_empty_project(author=author, filename=filename)
1538        elif isinstance(gpxfile, str): # TODO: create replacement for isinstance that checks if path exists
1539                                       # filename is valid, etc.
1540            if isinstance(filename, str): 
1541                self.filename = os.path.abspath(os.path.expanduser(filename)) # both are defined
1542            else: 
1543                self.filename = os.path.abspath(os.path.expanduser(gpxfile))
1544            # TODO set author
1545            self.data, self.names = LoadDictFromProjFile(gpxfile)
1546            self.update_ids()
1547        else:
1548            raise ValueError("Not sure what to do with gpxfile")
1549
1550    @classmethod
1551    def from_dict_and_names(cls, gpxdict, names, filename=None):
1552        """Creates a :class:`G2Project` directly from
1553        a dictionary and a list of names. If in doubt, do not use this.
1554
1555        :returns: a :class:`G2Project`
1556        """
1557        out = cls()
1558        if filename:
1559            filename = os.path.abspath(os.path.expanduser(filename))
1560            out.filename = filename
1561            gpxdict['Controls']['data']['LastSavedAs'] = filename
1562        else:
1563            try:
1564                out.filename = gpxdict['Controls']['data']['LastSavedAs']
1565            except KeyError:
1566                out.filename = None
1567        out.data = gpxdict
1568        out.names = names
1569
1570    def save(self, filename=None):
1571        """Saves the project, either to the current filename, or to a new file.
1572
1573        Updates self.filename if a new filename provided"""
1574        # TODO update LastSavedUsing ?
1575        if filename:
1576            filename = os.path.abspath(os.path.expanduser(filename))
1577            self.data['Controls']['data']['LastSavedAs'] = filename
1578            self.filename = filename
1579        elif not self.filename:
1580            raise AttributeError("No file name to save to")
1581        SaveDictToProjFile(self.data, self.names, self.filename)
1582
1583    def add_powder_histogram(self, datafile, iparams, phases=[], fmthint=None,
1584                                 databank=None, instbank=None):
1585        """Loads a powder data histogram into the project.
1586
1587        Automatically checks for an instrument parameter file, or one can be
1588        provided. Note that in unix fashion, "~" can be used to indicate the
1589        home directory (e.g. ~/G2data/data.fxye).
1590
1591        :param str datafile: The powder data file to read, a filename.
1592        :param str iparams: The instrument parameters file, a filename.
1593        :param list phases: A list of phases to link to the new histogram,
1594           phases can be references by object, name, rId or number.
1595           Alternately, use 'all' to link to all phases in the project.
1596        :param str fmthint: If specified, only importers where the format name
1597          (reader.formatName, as shown in Import menu) contains the
1598          supplied string will be tried as importers. If not specified, all
1599          importers consistent with the file extension will be tried
1600          (equivalent to "guess format" in menu).
1601        :param int databank: Specifies a dataset number to read, if file contains
1602          more than set of data. This should be 1 to read the first bank in
1603          the file (etc.) regardless of the number on the Bank line, etc.
1604          Default is None which means there should only be one dataset in the
1605          file.
1606        :param int instbank: Specifies an instrument parameter set to read, if
1607          the instrument parameter file contains more than set of parameters.
1608          This will match the INS # in an GSAS type file so it will typically
1609          be 1 to read the first parameter set in the file (etc.)
1610          Default is None which means there should only be one parameter set
1611          in the file.
1612
1613        :returns: A :class:`G2PwdrData` object representing
1614            the histogram
1615        """
1616        LoadG2fil()
1617        datafile = os.path.abspath(os.path.expanduser(datafile))
1618        iparams = os.path.abspath(os.path.expanduser(iparams))
1619        pwdrreaders = import_generic(datafile, Readers['Pwdr'],fmthint=fmthint,bank=databank)
1620        histname, new_names, pwdrdata = load_pwd_from_reader(
1621                                          pwdrreaders[0], iparams,
1622                                          [h.name for h in self.histograms()],bank=instbank)
1623        if histname in self.data:
1624            print("Warning - redefining histogram", histname)
1625        elif self.names[-1][0] == 'Phases':
1626            self.names.insert(-1, new_names)
1627        else:
1628            self.names.append(new_names)
1629        self.data[histname] = pwdrdata
1630        self.update_ids()
1631
1632        if phases == 'all':
1633            phases = self.phases()
1634        for phase in phases:
1635            phase = self.phase(phase)
1636            self.link_histogram_phase(histname, phase)
1637
1638        return self.histogram(histname)
1639
1640    def add_simulated_powder_histogram(self, histname, iparams, Tmin, Tmax, Tstep,
1641                                       wavelength=None, scale=None, phases=[]):
1642        """Loads a powder data histogram into the project.
1643
1644        Requires an instrument parameter file.
1645        Note that in unix fashion, "~" can be used to indicate the
1646        home directory (e.g. ~/G2data/data.prm). The instrument parameter file
1647        will determine if the histogram is x-ray, CW neutron, TOF, etc. as well
1648        as the instrument type.
1649
1650        :param str histname: A name for the histogram to be created.
1651        :param str iparams: The instrument parameters file, a filename.
1652        :param float Tmin: Minimum 2theta or TOF (ms) for dataset to be simulated
1653        :param float Tmax: Maximum 2theta or TOF (ms) for dataset to be simulated
1654        :param float Tstep: Step size in 2theta or TOF (ms) for dataset to be simulated       
1655        :param float wavelength: Wavelength for CW instruments, overriding the value
1656           in the instrument parameters file if specified.
1657        :param float scale: Histogram scale factor which multiplies the pattern. Note that
1658           simulated noise is added to the pattern, so that if the maximum intensity is
1659           small, the noise will mask the computed pattern. The scale
1660           needs to be a large number for CW neutrons.
1661           The default, None, provides a scale of 1 for x-rays and TOF; 10,000 for CW neutrons.
1662        :param list phases: Phases to link to the new histogram. Use proj.phases() to link to
1663           all defined phases.
1664
1665        :returns: A :class:`G2PwdrData` object representing the histogram
1666        """
1667        LoadG2fil()
1668        iparams = os.path.abspath(os.path.expanduser(iparams))
1669        if not os.path.exists(iparams):
1670            raise G2ScriptException("File does not exist:"+iparams)
1671        rd = G2obj.ImportPowderData( # Initialize a base class reader
1672            extensionlist=tuple(),
1673            strictExtension=False,
1674            formatName = 'Simulate dataset',
1675            longFormatName = 'Compute a simulated pattern')
1676        rd.powderentry[0] = '' # no filename
1677        rd.powderentry[2] = 1 # only one bank
1678        rd.comments.append('This is a dummy dataset for powder pattern simulation')
1679        rd.idstring = histname
1680        #Iparm1, Iparm2 = load_iprms(iparams, rd)
1681        if Tmax < Tmin:
1682            Tmin,Tmax = Tmax,Tmin
1683        Tstep = abs(Tstep)
1684        if 'TOF' in rd.idstring:
1685                N = (np.log(Tmax)-np.log(Tmin))/Tstep
1686                x = np.exp((np.arange(0,N))*Tstep+np.log(Tmin*1000.))
1687                N = len(x)
1688        else:           
1689                N = int((Tmax-Tmin)/Tstep)+1
1690                x = np.linspace(Tmin,Tmax,N,True)
1691                N = len(x)
1692        if N < 3:
1693            raise G2ScriptException("Error: Range is too small or step is too large, <3 points")
1694        rd.powderdata = [
1695            np.array(x), # x-axis values
1696            np.zeros_like(x), # powder pattern intensities
1697            np.ones_like(x), # 1/sig(intensity)^2 values (weights)
1698            np.zeros_like(x), # calc. intensities (zero)
1699            np.zeros_like(x), # calc. background (zero)
1700            np.zeros_like(x), # obs-calc profiles
1701            ]
1702        Tmin = rd.powderdata[0][0]
1703        Tmax = rd.powderdata[0][-1]
1704        histname, new_names, pwdrdata = load_pwd_from_reader(rd, iparams,
1705                                                            [h.name for h in self.histograms()])
1706        if histname in self.data:
1707            print("Warning - redefining histogram", histname)
1708        elif self.names[-1][0] == 'Phases':
1709            self.names.insert(-1, new_names)
1710        else:
1711            self.names.append(new_names)
1712        if scale is not None:
1713            pwdrdata['Sample Parameters']['Scale'][0] = scale
1714        elif pwdrdata['Instrument Parameters'][0]['Type'][0].startswith('PNC'):
1715            pwdrdata['Sample Parameters']['Scale'][0] = 10000.
1716        self.data[histname] = pwdrdata
1717        self.update_ids()
1718
1719        for phase in phases:
1720            phase = self.phase(phase)
1721            self.link_histogram_phase(histname, phase)
1722
1723        return self.histogram(histname)
1724   
1725    def add_phase(self, phasefile, phasename=None, histograms=[], fmthint=None):
1726        """Loads a phase into the project from a .cif file
1727
1728        :param str phasefile: The CIF file from which to import the phase.
1729        :param str phasename: The name of the new phase, or None for the default
1730        :param list histograms: The names of the histograms to associate with
1731            this phase. Use proj.histograms() to add to all histograms.
1732        :param str fmthint: If specified, only importers where the format name
1733          (reader.formatName, as shown in Import menu) contains the
1734          supplied string will be tried as importers. If not specified, all
1735          importers consistent with the file extension will be tried
1736          (equivalent to "guess format" in menu).
1737
1738        :returns: A :class:`G2Phase` object representing the
1739            new phase.
1740        """
1741        LoadG2fil()
1742        histograms = [self.histogram(h).name for h in histograms]
1743        phasefile = os.path.abspath(os.path.expanduser(phasefile))
1744
1745        # TODO handle multiple phases in a file
1746        phasereaders = import_generic(phasefile, Readers['Phase'], fmthint=fmthint)
1747        phasereader = phasereaders[0]
1748       
1749        phasename = phasename or phasereader.Phase['General']['Name']
1750        phaseNameList = [p.name for p in self.phases()]
1751        phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
1752        phasereader.Phase['General']['Name'] = phasename
1753
1754        if 'Phases' not in self.data:
1755            self.data[u'Phases'] = { 'data': None }
1756        assert phasename not in self.data['Phases'], "phase names should be unique"
1757        self.data['Phases'][phasename] = phasereader.Phase
1758
1759        if phasereader.Constraints:
1760            Constraints = self.data['Constraints']
1761            for i in phasereader.Constraints:
1762                if isinstance(i, dict):
1763                    if '_Explain' not in Constraints:
1764                        Constraints['_Explain'] = {}
1765                    Constraints['_Explain'].update(i)
1766                else:
1767                    Constraints['Phase'].append(i)
1768
1769        data = self.data['Phases'][phasename]
1770        generalData = data['General']
1771        SGData = generalData['SGData']
1772        NShkl = len(G2spc.MustrainNames(SGData))
1773        NDij = len(G2spc.HStrainNames(SGData))
1774        Super = generalData.get('Super', 0)
1775        if Super:
1776            SuperVec = np.array(generalData['SuperVec'][0])
1777        else:
1778            SuperVec = []
1779        UseList = data['Histograms']
1780
1781        for hist in histograms:
1782            self.link_histogram_phase(hist, phasename)
1783
1784        for obj in self.names:
1785            if obj[0] == 'Phases':
1786                phasenames = obj
1787                break
1788        else:
1789            phasenames = [u'Phases']
1790            self.names.append(phasenames)
1791        phasenames.append(phasename)
1792
1793        # TODO should it be self.filename, not phasefile?
1794        SetupGeneral(data, os.path.dirname(phasefile))
1795        self.index_ids()
1796
1797        self.update_ids()
1798        return self.phase(phasename)
1799
1800    def link_histogram_phase(self, histogram, phase):
1801        """Associates a given histogram and phase.
1802
1803        .. seealso::
1804
1805            :meth:`~G2Project.histogram`
1806            :meth:`~G2Project.phase`"""
1807        hist = self.histogram(histogram)
1808        phase = self.phase(phase)
1809
1810        generalData = phase['General']
1811
1812        if hist.name.startswith('HKLF '):
1813            raise NotImplementedError("HKLF not yet supported")
1814        elif hist.name.startswith('PWDR '):
1815            hist['Reflection Lists'][generalData['Name']] = {}
1816            UseList = phase['Histograms']
1817            SGData = generalData['SGData']
1818            NShkl = len(G2spc.MustrainNames(SGData))
1819            NDij = len(G2spc.HStrainNames(SGData))
1820            UseList[hist.name] = SetDefaultDData('PWDR', hist.name, NShkl=NShkl, NDij=NDij)
1821            UseList[hist.name]['hId'] = hist.id
1822            for key, val in [('Use', True), ('LeBail', False),
1823                             ('newLeBail', True),
1824                             ('Babinet', {'BabA': [0.0, False],
1825                                          'BabU': [0.0, False]})]:
1826                if key not in UseList[hist.name]:
1827                    UseList[hist.name][key] = val
1828        else:
1829            raise RuntimeError("Unexpected histogram" + hist.name)
1830
1831    def reload(self):
1832        """Reload self from self.filename"""
1833        data, names = LoadDictFromProjFile(self.filename)
1834        self.names = names
1835        # Need to deep copy the new data file data into the current tree,
1836        # so that any existing G2Phase, or G2PwdrData objects will still be
1837        # valid
1838        _deep_copy_into(from_=data, into=self.data)
1839
1840    def refine(self, newfile=None, printFile=None, makeBack=False):
1841        '''Invoke a refinement for the project. The project is written to
1842        the currently selected gpx file and then either a single or sequential refinement
1843        is performed depending on the setting of 'Seq Data' in Controls
1844        (set in :meth:`get_Controls`).
1845        '''
1846        seqSetting = self.data['Controls']['data'].get('Seq Data',[])
1847        if not seqSetting:
1848            self.index_ids()    # index_ids will automatically save the project
1849            # TODO: migrate to RefineCore G2strMain does not properly use printFile
1850            # G2strMain.RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
1851            #      calcControls,pawleyLookup,ifPrint,printFile,dlg)
1852            G2strMain.Refine(self.filename, makeBack=makeBack)
1853        else:
1854            self._seqrefine()
1855        self.reload() # get file from GPX
1856
1857    def _seqrefine(self):
1858        '''Perform a sequential refinement.
1859        '''
1860
1861        self.data['Controls']['data']['ShowCell'] = True
1862        # add to tree item to project, if not present
1863        if 'Sequential results' not in self.data:
1864            self.data['Sequential results'] = {'data':{}}
1865            self.names.append(['Sequential results'])
1866        self.index_ids()    # index_ids will automatically save the project
1867        #GSASIIpath.IPyBreak_base()
1868
1869        # check that constraints are OK
1870        errmsg, warnmsg = G2strIO.ReadCheckConstraints(self.filename)
1871        if errmsg:
1872            print('Refinement error',errmsg)
1873            raise Exception('Constraint error')
1874        if warnmsg:
1875            print(u'Conflict between refinment flag settings and constraints:\n'+
1876                  warnmsg+u'\nRefinement not possible')
1877            raise Exception('Constraint error')
1878        OK,Msg = G2strMain.SeqRefine(self.filename,None)
1879
1880    def histogram(self, histname):
1881        """Returns the histogram named histname, or None if it does not exist.
1882
1883        :param histname: The name of the histogram (str), or ranId or index.
1884        :returns: A :class:`G2PwdrData` object, or None if
1885            the histogram does not exist
1886
1887        .. seealso::
1888            :meth:`~G2Project.histograms`
1889            :meth:`~G2Project.phase`
1890            :meth:`~G2Project.phases`
1891        """
1892        if isinstance(histname, G2PwdrData):
1893            if histname.proj == self:
1894                return histname
1895            else:
1896                raise Exception('Histogram object (G2PwdrData) is not in current project')
1897        if histname in self.data:
1898            return G2PwdrData(self.data[histname], self, histname)
1899        try:
1900            # see if histname is an id or ranId
1901            histname = int(histname)
1902        except ValueError:
1903            return
1904
1905        for histogram in self.histograms():
1906            if histogram.id == histname or histogram.ranId == histname:
1907                return histogram
1908
1909    def histograms(self, typ=None):
1910        """Return a list of all histograms, as :class:`G2PwdrData` objects
1911
1912        For now this only finds Powder/Single Xtal histograms, since that is all that is
1913        currently implemented in this module.
1914
1915        :param ste typ: The prefix (type) the histogram such as 'PWDR '. If None
1916          (the default) all known histograms types are found.
1917        :returns: a list of objects
1918
1919        .. seealso::
1920            :meth:`~G2Project.histogram`
1921            :meth:`~G2Project.phase`
1922            :meth:`~G2Project.phases`
1923        """
1924        output = []
1925        # loop through each tree entry. If it is more than one level (more than one item in the
1926        # list of names). then it must be a histogram, unless labeled Phases or Restraints
1927        if typ is None:
1928            for obj in self.names:
1929                if obj[0].startswith('PWDR ') or obj[0].startswith('HKLF '): 
1930                    output.append(self.histogram(obj[0]))
1931        else:
1932            for obj in self.names:
1933                if len(obj) > 1 and obj[0].startswith(typ): 
1934                    output.append(self.histogram(obj[0]))
1935        return output
1936
1937    def phase(self, phasename):
1938        """
1939        Gives an object representing the specified phase in this project.
1940
1941        :param str phasename: A reference to the desired phase. Either the phase
1942            name (str), the phase's ranId, the phase's index (both int) or
1943            a phase object (:class:`G2Phase`)
1944        :returns: A :class:`G2Phase` object
1945        :raises: KeyError
1946
1947        .. seealso::
1948            :meth:`~G2Project.histograms`
1949            :meth:`~G2Project.phase`
1950            :meth:`~G2Project.phases`
1951            """
1952        if isinstance(phasename, G2Phase):
1953            if phasename.proj == self:
1954                return phasename
1955        phases = self.data['Phases']
1956        if phasename in phases:
1957            return G2Phase(phases[phasename], phasename, self)
1958
1959        try:
1960            # phasename should be phase index or ranId
1961            phasename = int(phasename)
1962        except ValueError:
1963            return
1964
1965        for phase in self.phases():
1966            if phase.id == phasename or phase.ranId == phasename:
1967                return phase
1968
1969    def phases(self):
1970        """
1971        Returns a list of all the phases in the project.
1972
1973        :returns: A list of :class:`G2Phase` objects
1974
1975        .. seealso::
1976            :meth:`~G2Project.histogram`
1977            :meth:`~G2Project.histograms`
1978            :meth:`~G2Project.phase`
1979            """
1980        for obj in self.names:
1981            if obj[0] == 'Phases':
1982                return [self.phase(p) for p in obj[1:]]
1983        return []
1984
1985    def _images(self):
1986        """Returns a list of all the phases in the project.
1987        """
1988        return [i[0] for i in self.names if i[0].startswith('IMG ')]
1989   
1990    def image(self, imageRef):
1991        """
1992        Gives an object representing the specified image in this project.
1993
1994        :param str imageRef: A reference to the desired image. Either the Image
1995          tree name (str), the image's index (int) or
1996          a image object (:class:`G2Image`)
1997        :returns: A :class:`G2Image` object
1998        :raises: KeyError
1999
2000        .. seealso::
2001            :meth:`~G2Project.images`
2002        """
2003        if isinstance(imageRef, G2Image):
2004            if imageRef.proj == self:
2005                return imageRef
2006            else:
2007                raise Exception("Image {} not in current selected project".format(imageRef.name))
2008        if imageRef in self._images():
2009            return G2Image(self.data[imageRef], imageRef, self)
2010
2011        try:
2012            # imageRef should be an index
2013            num = int(imageRef)
2014            imageRef = self._images()[num] 
2015            return G2Image(self.data[imageRef], imageRef, self)
2016        except ValueError:
2017            raise Exception("imageRef {} not an object, name or image index in current selected project"
2018                                .format(imageRef))
2019        except IndexError:
2020            raise Exception("imageRef {} out of range (max={}) in current selected project"
2021                                .format(imageRef,len(self._images())-1))
2022           
2023    def images(self):
2024        """
2025        Returns a list of all the images in the project.
2026
2027        :returns: A list of :class:`G2Image` objects
2028        """
2029        return [G2Image(self.data[i],i,self) for i in self._images()]
2030   
2031    def _pdfs(self):
2032        """Returns a list of all the PDF entries in the project.
2033        """
2034        return [i[0] for i in self.names if i[0].startswith('PDF ')]
2035   
2036    def pdf(self, pdfRef):
2037        """
2038        Gives an object representing the specified PDF entry in this project.
2039
2040        :param pdfRef: A reference to the desired image. Either the PDF
2041          tree name (str), the pdf's index (int) or
2042          a PDF object (:class:`G2PDF`)
2043        :returns: A :class:`G2PDF` object
2044        :raises: KeyError
2045
2046        .. seealso::
2047            :meth:`~G2Project.pdfs`
2048            :class:`~G2PDF`
2049        """
2050        if isinstance(pdfRef, G2PDF):
2051            if pdfRef.proj == self:
2052                return pdfRef
2053            else:
2054                raise Exception("PDF {} not in current selected project".format(pdfRef.name))
2055        if pdfRef in self._pdfs():
2056            return G2PDF(self.data[pdfRef], pdfRef, self)
2057
2058        try:
2059            # pdfRef should be an index
2060            num = int(pdfRef)
2061            pdfRef = self._pdfs()[num] 
2062            return G2PDF(self.data[pdfRef], pdfRef, self)
2063        except ValueError:
2064            raise Exception("pdfRef {} not an object, name or PDF index in current selected project"
2065                                .format(pdfRef))
2066        except IndexError:
2067            raise Exception("pdfRef {} out of range (max={}) in current selected project"
2068                                .format(pdfRef,len(self._images())-1))
2069    def pdfs(self):
2070        """
2071        Returns a list of all the PDFs in the project.
2072
2073        :returns: A list of :class:`G2PDF` objects
2074        """
2075        return [G2PDF(self.data[i],i,self) for i in self._pdfs()]
2076       
2077    def add_PDF(self, prmfile, histogram):
2078        '''Creates a PDF entry that can be used to compute a PDF.
2079        Note that this command places an entry in the project,
2080        but :meth:`G2PDF.calculate` must be used to actually perform
2081        the computation.
2082
2083        :param str datafile: The powder data file to read, a filename.
2084        :param histogram: A reference to a histogram,
2085          which can be reference by object, name, or number.
2086        :returns: A :class:`G2PDF` object for the PDF entry
2087        '''
2088       
2089        LoadG2fil()
2090        PDFname = 'PDF ' + self.histogram(histogram).name[4:]
2091        peaks = {'Limits':[1.,5.],'Background':[2,[0.,-0.2*np.pi],False],'Peaks':[]}
2092        Controls = {
2093            'Sample':{'Name':self.histogram(histogram).name,'Mult':1.0},
2094            'Sample Bkg.':{'Name':'','Mult':-1.0,'Refine':False},
2095            'Container':{'Name':'','Mult':-1.0,'Refine':False},
2096            'Container Bkg.':{'Name':'','Mult':-1.0},
2097            'ElList':{},
2098            'Geometry':'Cylinder','Diam':1.0,'Pack':0.50,'Form Vol':0.0,'Flat Bkg':0,
2099            'DetType':'Area detector','ObliqCoeff':0.2,'Ruland':0.025,'QScaleLim':[20,25],
2100            'Lorch':False,'BackRatio':0.0,'Rmax':100.,'noRing':False,'IofQmin':1.0,'Rmin':1.0,
2101            'I(Q)':[],'S(Q)':[],'F(Q)':[],'G(R)':[]}
2102
2103        fo = open(prmfile,'r')
2104        S = fo.readline()
2105        while S:
2106            if '#' in S:
2107                S = fo.readline()
2108                continue
2109            key,val = S.split(':',1)
2110            try:
2111                Controls[key] = eval(val)
2112            except:
2113                Controls[key] = val.strip()
2114            S = fo.readline()
2115        fo.close()
2116        Controls['Sample']['Name'] = self.histogram(histogram).name
2117        for i in 'Sample Bkg.','Container','Container Bkg.':
2118            Controls[i]['Name'] = ''
2119        PDFdict = {'data':None,'PDF Controls':Controls, 'PDF Peaks':peaks}
2120        self.names.append([PDFname]+['PDF Controls', 'PDF Peaks'])
2121        self.data[PDFname] = PDFdict
2122        print('Adding "{}" to project'.format(PDFname))
2123        return G2PDF(self.data[PDFname], PDFname, self)
2124
2125    def seqref(self):
2126        """
2127        Returns a sequential refinement results object, if present
2128
2129        :returns: A :class:`G2SeqRefRes` object or None if not present
2130        """
2131        if 'Sequential results' not in self.data: return
2132        return G2SeqRefRes(self.data['Sequential results']['data'], self)
2133   
2134    def update_ids(self):
2135        """Makes sure all phases and histograms have proper hId and pId"""
2136        # Translated from GetUsedHistogramsAndPhasesfromTree,
2137        #   GSASIIdataGUI.py:4107
2138        for i, h in enumerate(self.histograms()):
2139            h.id = i
2140        for i, p in enumerate(self.phases()):
2141            p.id = i
2142
2143    def do_refinements(self, refinements, histogram='all', phase='all',
2144                       outputnames=None, makeBack=False):
2145        """Conducts one or a series of refinements according to the
2146           input provided in parameter refinements. This is a wrapper
2147           around :meth:`iter_refinements`
2148
2149        :param list refinements: A list of dictionaries specifiying changes to be made to
2150            parameters before refinements are conducted.
2151            See the :ref:`Refinement_recipe` section for how this is defined.
2152        :param str histogram: Name of histogram for refinements to be applied
2153            to, or 'all'; note that this can be overridden for each refinement
2154            step via a "histograms" entry in the dict.
2155        :param str phase: Name of phase for refinements to be applied to, or
2156            'all'; note that this can be overridden for each refinement
2157            step via a "phases" entry in the dict.
2158        :param list outputnames: Provides a list of project (.gpx) file names
2159            to use for each refinement step (specifying None skips the save step).
2160            See :meth:`save`.
2161            Note that this can be overridden using an "output" entry in the dict.
2162        :param bool makeBack: determines if a backup ).bckX.gpx) file is made
2163            before a refinement is performed. The default is False.
2164           
2165        To perform a single refinement without changing any parameters, use this
2166        call:
2167
2168        .. code-block::  python
2169       
2170            my_project.do_refinements([])
2171        """
2172       
2173        for proj in self.iter_refinements(refinements, histogram, phase,
2174                                          outputnames, makeBack):
2175            pass
2176        return self
2177
2178    def iter_refinements(self, refinements, histogram='all', phase='all',
2179                         outputnames=None, makeBack=False):
2180        """Conducts a series of refinements, iteratively. Stops after every
2181        refinement and yields this project, to allow error checking or
2182        logging of intermediate results. Parameter use is the same as for
2183        :meth:`do_refinements` (which calls this method).
2184
2185        >>> def checked_refinements(proj):
2186        ...     for p in proj.iter_refinements(refs):
2187        ...         # Track intermediate results
2188        ...         log(p.histogram('0').residuals)
2189        ...         log(p.phase('0').get_cell())
2190        ...         # Check if parameter diverged, nonsense answer, or whatever
2191        ...         if is_something_wrong(p):
2192        ...             raise Exception("I need a human!")
2193
2194           
2195        """
2196        if outputnames:
2197            if len(refinements) != len(outputnames):
2198                raise ValueError("Should have same number of outputs to"
2199                                 "refinements")
2200        else:
2201            outputnames = [None for r in refinements]
2202
2203        for output, refinedict in zip(outputnames, refinements):
2204            if 'histograms' in refinedict:
2205                hist = refinedict['histograms']
2206            else:
2207                hist = histogram
2208            if 'phases' in refinedict:
2209                ph = refinedict['phases']
2210            else:
2211                ph = phase
2212            if 'output' in refinedict:
2213                output = refinedict['output']
2214            self.set_refinement(refinedict, hist, ph)
2215            # Handle 'once' args - refinements that are disabled after this
2216            # refinement
2217            if 'once' in refinedict:
2218                temp = {'set': refinedict['once']}
2219                self.set_refinement(temp, hist, ph)
2220
2221            if output:
2222                self.save(output)
2223
2224            if 'skip' not in refinedict:
2225                self.refine(makeBack=makeBack)
2226            yield self
2227
2228            # Handle 'once' args - refinements that are disabled after this
2229            # refinement
2230            if 'once' in refinedict:
2231                temp = {'clear': refinedict['once']}
2232                self.set_refinement(temp, hist, ph)
2233            if 'call' in refinedict:
2234                fxn = refinedict['call']
2235                if callable(fxn):
2236                    fxn(*refinedict.get('callargs',[self]))
2237                elif callable(eval(fxn)):
2238                    eval(fxn)(*refinedict.get('callargs',[self]))
2239                else:
2240                    raise G2ScriptException("Error: call value {} is not callable".format(fxn))
2241
2242    def set_refinement(self, refinement, histogram='all', phase='all'):
2243        """Apply specified refinements to a given histogram(s) or phase(s).
2244
2245        :param dict refinement: The refinements to be conducted
2246        :param histogram: Specifies either 'all' (default), a single histogram or
2247          a list of histograms. Histograms may be specified as histogram objects
2248          (see :class:`G2PwdrData`), the histogram name (str) or the index number (int)
2249          of the histogram in the project, numbered starting from 0.
2250          Omitting the parameter or the string 'all' indicates that parameters in
2251          all histograms should be set.
2252        :param phase: Specifies either 'all' (default), a single phase or
2253          a list of phases. Phases may be specified as phase objects
2254          (see :class:`G2Phase`), the phase name (str) or the index number (int)
2255          of the phase in the project, numbered starting from 0.
2256          Omitting the parameter or the string 'all' indicates that parameters in
2257          all phases should be set.
2258
2259        Note that refinement parameters are categorized as one of three types:
2260
2261        1. Histogram parameters
2262        2. Phase parameters
2263        3. Histogram-and-Phase (HAP) parameters
2264       
2265        .. seealso::
2266            :meth:`G2PwdrData.set_refinements`
2267            :meth:`G2PwdrData.clear_refinements`
2268            :meth:`G2Phase.set_refinements`
2269            :meth:`G2Phase.clear_refinements`
2270            :meth:`G2Phase.set_HAP_refinements`
2271            :meth:`G2Phase.clear_HAP_refinements`"""
2272
2273        if histogram == 'all':
2274            hists = self.histograms()
2275        elif isinstance(histogram, list) or isinstance(histogram, tuple):
2276            hists = []
2277            for h in histogram:
2278                if isinstance(h, str) or isinstance(h, int):
2279                    hists.append(self.histogram(h))
2280                else:
2281                    hists.append(h)
2282        elif isinstance(histogram, str) or isinstance(histogram, int):
2283            hists = [self.histogram(histogram)]
2284        else:
2285            hists = [histogram]
2286
2287        if phase == 'all':
2288            phases = self.phases()
2289        elif isinstance(phase, list) or isinstance(phase, tuple):
2290            phases = []
2291            for ph in phase:
2292                if isinstance(ph, str) or isinstance(ph, int):
2293                    phases.append(self.phase(ph))
2294                else:
2295                    phases.append(ph)
2296        elif isinstance(phase, str) or isinstance(phase, int):
2297            phases = [self.phase(phase)]
2298        else:
2299            phases = [phase]
2300
2301        pwdr_set = {}
2302        phase_set = {}
2303        hap_set = {}
2304        for key, val in refinement.get('set', {}).items():
2305            # Apply refinement options
2306            if G2PwdrData.is_valid_refinement_key(key):
2307                pwdr_set[key] = val
2308            elif G2Phase.is_valid_refinement_key(key):
2309                phase_set[key] = val
2310            elif G2Phase.is_valid_HAP_refinement_key(key):
2311                hap_set[key] = val
2312            else:
2313                raise ValueError("Unknown refinement key", key)
2314
2315        for hist in hists:
2316            hist.set_refinements(pwdr_set)
2317        for phase in phases:
2318            phase.set_refinements(phase_set)
2319        for phase in phases:
2320            phase.set_HAP_refinements(hap_set, hists)
2321
2322        pwdr_clear = {}
2323        phase_clear = {}
2324        hap_clear = {}
2325        for key, val in refinement.get('clear', {}).items():
2326            # Clear refinement options
2327            if G2PwdrData.is_valid_refinement_key(key):
2328                pwdr_clear[key] = val
2329            elif G2Phase.is_valid_refinement_key(key):
2330                phase_clear[key] = val
2331            elif G2Phase.is_valid_HAP_refinement_key(key):
2332                hap_set[key] = val
2333            else:
2334                raise ValueError("Unknown refinement key", key)
2335
2336        for hist in hists:
2337            hist.clear_refinements(pwdr_clear)
2338        for phase in phases:
2339            phase.clear_refinements(phase_clear)
2340        for phase in phases:
2341            phase.clear_HAP_refinements(hap_clear, hists)
2342
2343    def index_ids(self):
2344        self.save()
2345        return G2strIO.GetUsedHistogramsAndPhases(self.filename)
2346
2347    def add_constraint_raw(self, cons_scope, constr):
2348        """Adds a constraint of type consType to the project.
2349        cons_scope should be one of "Hist", "Phase", "HAP", or "Global".
2350
2351        WARNING it does not check the constraint is well-constructed"""
2352        constrs = self.data['Constraints']['data']
2353        if 'Global' not in constrs:
2354            constrs['Global'] = []
2355        constrs[cons_scope].append(constr)
2356
2357    def hold_many(self, vars, type):
2358        """Apply holds for all the variables in vars, for constraint of a given type.
2359
2360        type is passed directly to add_constraint_raw as consType
2361
2362        :param list vars: A list of variables to hold. Either :class:`GSASIIobj.G2VarObj` objects,
2363            string variable specifiers, or arguments for :meth:`make_var_obj`
2364        :param str type: A string constraint type specifier. See
2365            :class:`G2Project.add_constraint_raw`
2366
2367        """
2368        for var in vars:
2369            if isinstance(var, str):
2370                var = self.make_var_obj(var)
2371            elif not isinstance(var, G2obj.G2VarObj):
2372                var = self.make_var_obj(*var)
2373            self.add_constraint_raw(type, [[1.0, var], None, None, 'h'])
2374
2375    def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
2376                     reloadIdx=True):
2377        """Wrapper to create a G2VarObj. Takes either a string representation ("p:h:name:a")
2378        or individual names of phase, histogram, varname, and atomId.
2379
2380        Automatically converts string phase, hist, or atom names into the ID required
2381        by G2VarObj."""
2382
2383        if reloadIdx:
2384            self.index_ids()
2385
2386        # If string representation, short circuit
2387        if hist is None and varname is None and atomId is None:
2388            if isinstance(phase, str) and ':' in phase:
2389                return G2obj.G2VarObj(phase)
2390
2391        # Get phase index
2392        phaseObj = None
2393        if isinstance(phase, G2Phase):
2394            phaseObj = phase
2395            phase = G2obj.PhaseRanIdLookup[phase.ranId]
2396        elif phase in self.data['Phases']:
2397            phaseObj = self.phase(phase)
2398            phaseRanId = phaseObj.ranId
2399            phase = G2obj.PhaseRanIdLookup[phaseRanId]
2400
2401        # Get histogram index
2402        if isinstance(hist, G2PwdrData):
2403            hist = G2obj.HistRanIdLookup[hist.ranId]
2404        elif hist in self.data:
2405            histRanId = self.histogram(hist).ranId
2406            hist = G2obj.HistRanIdLookup[histRanId]
2407
2408        # Get atom index (if any)
2409        if isinstance(atomId, G2AtomRecord):
2410            atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
2411        elif phaseObj:
2412            atomObj = phaseObj.atom(atomId)
2413            if atomObj:
2414                atomRanId = atomObj.ranId
2415                atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
2416
2417        return G2obj.G2VarObj(phase, hist, varname, atomId)
2418
2419    def add_image(self, imagefile, fmthint=None, defaultImage=None):
2420        """Load an image into a project
2421
2422        :param str imagefile: The image file to read, a filename.
2423        :param str fmthint: If specified, only importers where the format name
2424          (reader.formatName, as shown in Import menu) contains the
2425          supplied string will be tried as importers. If not specified, all
2426          importers consistent with the file extension will be tried
2427          (equivalent to "guess format" in menu).
2428        :param str defaultImage: The name of an image to use as a default for
2429          setting parameters for the image file to read.
2430
2431        :returns: a list of G2Image object for the added image(s) [this routine
2432          has not yet been tested with files with multiple images...]
2433        """
2434        LoadG2fil()
2435        imagefile = os.path.abspath(os.path.expanduser(imagefile))
2436        readers = import_generic(imagefile, Readers['Image'], fmthint=fmthint)
2437        objlist = []
2438        for rd in readers:
2439            if rd.SciPy:        #was default read by scipy; needs 1 time fixes
2440                print('TODO: Image {} read by SciPy. Parameters likely need editing'.format(imagefile))
2441                #see this: G2IO.EditImageParms(self,rd.Data,rd.Comments,rd.Image,imagefile)
2442                rd.SciPy = False
2443            rd.readfilename = imagefile
2444            if rd.repeatcount == 1 and not rd.repeat: # skip image number if only one in set
2445                rd.Data['ImageTag'] = None
2446            else:
2447                rd.Data['ImageTag'] = rd.repeatcount
2448            rd.Data['formatName'] = rd.formatName
2449            if rd.sumfile:
2450                rd.readfilename = rd.sumfile
2451            # Load generic metadata, as configured
2452            G2fil.GetColumnMetadata(rd)
2453            # Code from G2IO.LoadImage2Tree(rd.readfilename,self,rd.Comments,rd.Data,rd.Npix,rd.Image)
2454            Imax = np.amax(rd.Image)
2455            ImgNames = [i[0] for i in self.names if i[0].startswith('IMG ')]
2456            TreeLbl = 'IMG '+os.path.basename(imagefile)
2457            ImageTag = rd.Data.get('ImageTag')
2458            if ImageTag:
2459                TreeLbl += ' #'+'%04d'%(ImageTag)
2460                imageInfo = (imagefile,ImageTag)
2461            else:
2462                imageInfo = imagefile
2463            TreeName = G2obj.MakeUniqueLabel(TreeLbl,ImgNames)
2464            # MT dict to contain image info
2465            ImgDict = {}
2466            ImgDict['data'] = [rd.Npix,imageInfo]
2467            ImgDict['Comments'] = rd.Comments
2468            if defaultImage:
2469                if isinstance(defaultImage, G2Image):
2470                    if defaultImage.proj == self:
2471                        defaultImage = self.data[defaultImage.name]['data']
2472                    else:
2473                        raise Exception("Image {} not in current selected project".format(defaultImage.name))
2474                elif defaultImage in self._images():
2475                    defaultImage = self.data[defaultImage]['data']
2476                else:
2477                    defaultImage = None
2478            Data = rd.Data
2479            if defaultImage:
2480                Data = copy.copy(defaultImage)
2481                Data['showLines'] = True
2482                Data['ring'] = []
2483                Data['rings'] = []
2484                Data['cutoff'] = 10.
2485                Data['pixLimit'] = 20
2486                Data['edgemin'] = 100000000
2487                Data['calibdmin'] = 0.5
2488                Data['calibskip'] = 0
2489                Data['ellipses'] = []
2490                Data['calibrant'] = ''
2491                Data['GonioAngles'] = [0.,0.,0.]
2492                Data['DetDepthRef'] = False
2493            else:
2494                Data['type'] = 'PWDR'
2495                Data['color'] = GSASIIpath.GetConfigValue('Contour_color','Paired')
2496                if 'tilt' not in Data:          #defaults if not preset in e.g. Bruker importer
2497                    Data['tilt'] = 0.0
2498                    Data['rotation'] = 0.0
2499                    Data['pixLimit'] = 20
2500                    Data['calibdmin'] = 0.5
2501                    Data['cutoff'] = 10.
2502                Data['showLines'] = False
2503                Data['calibskip'] = 0
2504                Data['ring'] = []
2505                Data['rings'] = []
2506                Data['edgemin'] = 100000000
2507                Data['ellipses'] = []
2508                Data['GonioAngles'] = [0.,0.,0.]
2509                Data['DetDepth'] = 0.
2510                Data['DetDepthRef'] = False
2511                Data['calibrant'] = ''
2512                Data['IOtth'] = [5.0,50.0]
2513                Data['LRazimuth'] = [0.,180.]
2514                Data['azmthOff'] = 0.0
2515                Data['outChannels'] = 2500
2516                Data['outAzimuths'] = 1
2517                Data['centerAzm'] = False
2518                Data['fullIntegrate'] = GSASIIpath.GetConfigValue('fullIntegrate',True)
2519                Data['setRings'] = False
2520                Data['background image'] = ['',-1.0]                           
2521                Data['dark image'] = ['',-1.0]
2522                Data['Flat Bkg'] = 0.0
2523                Data['Oblique'] = [0.5,False]
2524            Data['setDefault'] = False
2525            Data['range'] = [(0,Imax),[0,Imax]]
2526            ImgDict['Image Controls'] = Data
2527            ImgDict['Masks'] = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],
2528                                'Frames':[],'Thresholds':[(0,Imax),[0,Imax]]}
2529            ImgDict['Stress/Strain']  = {'Type':'True','d-zero':[],'Sample phi':0.0,
2530                                             'Sample z':0.0,'Sample load':0.0}
2531            self.names.append([TreeName]+['Comments','Image Controls','Masks','Stress/Strain'])
2532            self.data[TreeName] = ImgDict
2533            del rd.Image
2534            objlist.append(G2Image(self.data[TreeName], TreeName, self))
2535        return objlist
2536   
2537    def get_Controls(self, control):
2538        '''Return project controls settings
2539
2540        :param str control: the item to be returned. See below for allowed values.
2541        :returns: The value for the control.
2542
2543        Allowed values for parameter control:
2544
2545            * cycles: the maximum number of cycles (returns int)
2546            * sequential: the histograms used for a sequential refinement as a list
2547              of histogram names or an empty list when in non-sequential mode.
2548            * seqCopy: returns True or False. True indicates that results from
2549              each sequential fit are used as the starting point for the next
2550              histogram.
2551            * Anything else returns the value in the Controls dict, if present. An
2552              exception is raised if the control value is not present.
2553
2554        .. seealso::
2555            :meth:`set_Controls`
2556        '''
2557        if control == 'cycles':
2558            return self.data['Controls']['data']['max cyc']
2559        elif control == 'sequential':
2560            return self.data['Controls']['data']['Seq Data']
2561        elif control in self.data['Controls']['data']:
2562            return self.data['Controls']['data'][control]
2563        elif control == 'seqCopy':
2564            return self.data['Controls']['data']['Copy2Next']
2565        else:
2566            print('Defined Controls:',self.data['Controls']['data'].keys())
2567            raise Exception('{} is an invalid control value'.format(control))
2568       
2569    def set_Controls(self, control, value):
2570        '''Set project controls
2571
2572        :param str control: the item to be set. See below for allowed values.
2573        :param value: the value to be set.
2574
2575        Allowed values for parameter control:
2576
2577            * cycles: sets the maximum number of cycles (value must be int)
2578            * sequential: sets the histograms to be used for a sequential refinement.
2579              Use an empty list to turn off sequential fitting.
2580              The values in the list may be the name of the histogram (a str), or
2581              a ranId or index (int values), see :meth:`histogram`.
2582            * seqCopy: when True, the results from each sequential fit are used as
2583              the starting point for the next. After each fit is is set to False.
2584              Ignored for non-sequential fits.
2585
2586        .. seealso::
2587            :meth:`get_Controls`
2588        '''
2589        if control == 'cycles':
2590            self.data['Controls']['data']['max cyc'] = int(value)
2591        elif control == 'seqCopy':
2592            self.data['Controls']['data']['Copy2Next'] = bool(value)
2593        elif control == 'sequential':
2594            histlist = []
2595            for i,j in enumerate(value):
2596                h = self.histogram(j)
2597                if h:
2598                    histlist.append(h.name)
2599                else:
2600                    raise Exception('item #{} ({}) is an invalid histogram value'
2601                                        .format(i,j))
2602            self.data['Controls']['data']['Seq Data'] = histlist
2603        else:
2604            raise Exception('{} is an invalid control value'.format(control))
2605       
2606    def copyHistParms(self,sourcehist,targethistlist='all',modelist='all'):
2607        '''Copy histogram information from one histogram to others
2608
2609        :param sourcehist: is a histogram object (:class:`G2PwdrData`) or
2610            a histogram name or the index number of the histogram
2611        :param list targethistlist: a list of histograms where each item in the
2612            list can be a histogram object (:class:`G2PwdrData`),
2613            a histogram name or the index number of the histogram.
2614            if the string 'all' (default value), then all histograms in
2615            the project are used.
2616        :param list modelist: May be a list of sections to copy, which may
2617           include 'Background', 'Instrument Parameters', 'Limits' and
2618           'Sample Parameters' (items may be shortened to uniqueness and
2619           capitalization is ignored, so ['b','i','L','s'] will work.)
2620           The default value, 'all' causes the listed sections to
2621
2622        '''
2623        sections = ('Background','Instrument Parameters','Limits',
2624                        'Sample Parameters')
2625        hist_in = self.histogram(sourcehist)
2626        if not hist_in:
2627            raise Exception('{} is not a valid histogram'.format(sourcehist))
2628        if targethistlist == "all":
2629            targethistlist = self.histograms()
2630        if 'all' in modelist:
2631            copysections = sections
2632        else:
2633            copysections = set()
2634            for s in sections:
2635                for m in modelist:
2636                    if s.lower().startswith(m.lower()):
2637                        copysections.add(s)
2638        for h in targethistlist:
2639            hist_out = self.histogram(h)
2640            if not hist_out:
2641                raise Exception('{} is not a valid histogram'.format(h))
2642            for key in copysections: 
2643                hist_out[key] = copy.deepcopy(hist_in[key])
2644
2645    def get_VaryList(self):
2646        '''Returns a list of the refined variables in the
2647        last refinement cycle
2648
2649        :returns: a list of variables or None if no refinement has been
2650          performed.
2651        '''
2652        try:
2653            return self['Covariance']['data']['varyList']
2654        except:
2655            return
2656
2657    def get_ParmList(self):
2658        '''Returns a list of all the parameters defined in the
2659        last refinement cycle
2660
2661        :returns: a list of parameters or None if no refinement has been
2662          performed.
2663        '''
2664        try:
2665            return list(self['Covariance']['data']['parmDict'].keys())
2666        except:
2667            return
2668       
2669    def get_Variable(self,var):
2670        '''Returns the value and standard uncertainty (esd) for a variable
2671        parameters, as defined in the last refinement cycle
2672
2673        :param str var: a variable name of form '<p>:<h>:<name>', such as
2674          ':0:Scale'
2675        :returns: (value,esd) if the parameter is refined or
2676          (value, None) if the variable is in a constraint or is not
2677          refined or None if the parameter is not found.
2678        '''
2679        if var not in self['Covariance']['data']['parmDict']:
2680            return None
2681        val = self['Covariance']['data']['parmDict'][var]
2682        try:
2683            pos = self['Covariance']['data']['varyList'].index(var)
2684            esd = np.sqrt(self['Covariance']['data']['covMatrix'][pos,pos])
2685            return (val,esd)
2686        except ValueError:
2687            return (val,None)
2688
2689    def get_Covariance(self,varList):
2690        '''Returns the values and covariance matrix for a series of variable
2691        parameters. as defined in the last refinement cycle
2692
2693        :param tuple varList: a list of variable names of form '<p>:<h>:<name>'
2694        :returns: (valueList,CovMatrix) where valueList contains the (n) values
2695          in the same order as varList (also length n) and CovMatrix is a
2696          (n x n) matrix. If any variable name is not found in the varyList
2697          then None is returned.
2698
2699        Use this code, where sig provides standard uncertainties for
2700        parameters and where covArray provides the correlation between
2701        off-diagonal terms::
2702
2703            sig = np.sqrt(np.diag(covMatrix))
2704            xvar = np.outer(sig,np.ones_like(sig))
2705            covArray = np.divide(np.divide(covMatrix,xvar),xvar.T)
2706
2707        '''
2708        missing = [i for i in varList if i not in self['Covariance']['data']['varyList']]
2709        if missing:
2710            print('Variable(s) {} were not found in the varyList'.format(missing))
2711            return None
2712        vals = [self['Covariance']['data']['parmDict'][i] for i in varList]
2713        import GSASIImath as G2mth
2714        cov = G2mth.getVCov(varList,
2715                            self['Covariance']['data']['varyList'],
2716                            self['Covariance']['data']['covMatrix'])
2717        return (vals,cov)
2718
2719class G2AtomRecord(G2ObjectWrapper):
2720    """Wrapper for an atom record. Has convenient accessors via @property.
2721
2722    Available accessors: label, type, refinement_flags, coordinates,
2723        occupancy, ranId, id, adp_flag, uiso
2724
2725    Example:
2726
2727    >>> atom = some_phase.atom("O3")
2728    >>> # We can access the underlying data format
2729    >>> atom.data
2730    ['O3', 'O-2', '', ... ]
2731    >>> # We can also use wrapper accessors
2732    >>> atom.coordinates
2733    (0.33, 0.15, 0.5)
2734    >>> atom.refinement_flags
2735    u'FX'
2736    >>> atom.ranId
2737    4615973324315876477
2738    >>> atom.occupancy
2739    1.0
2740    """
2741    def __init__(self, data, indices, proj):
2742        self.data = data
2743        self.cx, self.ct, self.cs, self.cia = indices
2744        self.proj = proj
2745
2746    @property
2747    def label(self):
2748        '''Get the associated atom's label
2749        '''
2750        return self.data[self.ct-1]
2751
2752    @property
2753    def type(self):
2754        '''Get the associated atom's type
2755        '''
2756        return self.data[self.ct]
2757
2758    @property
2759    def refinement_flags(self):
2760        '''Get or set refinement flags for the associated atom
2761        '''
2762        return self.data[self.ct+1]
2763
2764    @refinement_flags.setter
2765    def refinement_flags(self, other):
2766        # Automatically check it is a valid refinement
2767        for c in other:
2768            if c not in ' FXU':
2769                raise ValueError("Invalid atom refinement: ", other)
2770        self.data[self.ct+1] = other
2771
2772    @property
2773    def coordinates(self):
2774        '''Get the associated atom's coordinates
2775        '''
2776        return tuple(self.data[self.cx:self.cx+3])
2777
2778    @property
2779    def occupancy(self):
2780        '''Get or set the associated atom's occupancy fraction
2781        '''
2782        return self.data[self.cx+3]
2783
2784    @occupancy.setter
2785    def occupancy(self, val):
2786        self.data[self.cx+3] = float(val)
2787
2788    @property
2789    def ranId(self):
2790        '''Get the associated atom's Random Id number
2791        '''
2792        return self.data[self.cia+8]
2793
2794    @property
2795    def adp_flag(self):
2796        '''Get the associated atom's iso/aniso setting, 'I' or 'A'
2797        '''
2798        # Either 'I' or 'A'
2799        return self.data[self.cia]
2800
2801    @property
2802    def uiso(self):
2803        '''Get or set the associated atom's Uiso or Uaniso value(s)
2804        '''
2805        if self.adp_flag == 'I':
2806            return self.data[self.cia+1]
2807        else:
2808            return self.data[self.cia+2:self.cia+8]
2809
2810    @uiso.setter
2811    def uiso(self, value):
2812        if self.adp_flag == 'I':
2813            self.data[self.cia+1] = float(value)
2814        else:
2815            assert len(value) == 6
2816            self.data[self.cia+2:self.cia+8] = [float(v) for v in value]
2817
2818class G2PwdrData(G2ObjectWrapper):
2819    """Wraps a Powder Data Histogram."""
2820    def __init__(self, data, proj, name):
2821        self.data = data
2822        self.proj = proj
2823        self.name = name
2824
2825    @staticmethod
2826    def is_valid_refinement_key(key):
2827        valid_keys = ['Limits', 'Sample Parameters', 'Background',
2828                      'Instrument Parameters']
2829        return key in valid_keys
2830
2831    #@property
2832    #def name(self):
2833    #    return self.data['data'][-1]
2834
2835    @property
2836    def ranId(self):
2837        return self.data['data'][0]['ranId']
2838
2839    @property
2840    def residuals(self):
2841        '''Provides a dictionary with with the R-factors for this histogram.
2842        Includes the weighted and unweighted profile terms (R, Rb, wR, wRb, wRmin)
2843        as well as the Bragg R-values for each phase (ph:H:Rf and ph:H:Rf^2).
2844        '''
2845        data = self.data['data'][0]
2846        return {key: data[key] for key in data
2847                if key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']
2848                   or ':' in key}
2849
2850    @property
2851    def InstrumentParameters(self):
2852        '''Provides a dictionary with with the Instrument Parameters
2853        for this histogram.
2854        '''
2855        return self.data['Instrument Parameters'][0]
2856
2857    @property
2858    def SampleParameters(self):
2859        '''Provides a dictionary with with the Sample Parameters
2860        for this histogram.
2861        '''
2862        return self.data['Sample Parameters']
2863
2864    @property
2865    def Background(self):
2866        '''Provides a list with with the Background parameters
2867        for this histogram.
2868
2869        :returns: list containing a list and dict with background values
2870        '''
2871        return self.data['Background']
2872
2873    def add_back_peak(self,pos,int,sig,gam,refflags=[]):
2874        '''Adds a background peak to the Background parameters
2875       
2876        :param float pos: position of peak, a 2theta or TOF value
2877        :param float int: integrated intensity of background peak, usually large
2878        :param float sig: Gaussian width of background peak, usually large
2879        :param float gam: Lorentzian width of background peak, usually unused (small)
2880        :param list refflags: a list of 1 to 4 boolean refinement flags for
2881            pos,int,sig & gam, respectively (use [0,1] to refine int only).
2882            Defaults to [] which means nothing is refined.
2883        '''
2884        if 'peaksList' not in self.Background[1]:
2885            self.Background[1]['peaksList'] = []
2886        flags = 4*[False]
2887        for i,f in enumerate(refflags):
2888            if i>3: break
2889            flags[i] = bool(f)
2890        bpk = []
2891        for i,j in zip((pos,int,sig,gam),flags):
2892            bpk += [float(i),j]
2893        self.Background[1]['peaksList'].append(bpk)
2894        self.Background[1]['nPeaks'] = len(self.Background[1]['peaksList'])
2895
2896    def del_back_peak(self,peaknum):
2897        '''Removes a background peak from the Background parameters
2898       
2899        :param int peaknum: the number of the peak (starting from 0)
2900        '''
2901        npks = self.Background[1].get('nPeaks',0)
2902        if peaknum >= npks:
2903            raise Exception('peak {} not found in histogram {}'.format(peaknum,self.name))
2904        del self.Background[1]['peaksList'][peaknum]
2905        self.Background[1]['nPeaks'] = len(self.Background[1]['peaksList'])
2906       
2907    def ref_back_peak(self,peaknum,refflags=[]):
2908        '''Sets refinement flag for a background peak
2909       
2910        :param int peaknum: the number of the peak (starting from 0)
2911        :param list refflags: a list of 1 to 4 boolean refinement flags for
2912            pos,int,sig & gam, respectively. If a flag is not specified
2913            it defaults to False (use [0,1] to refine int only).
2914            Defaults to [] which means nothing is refined.
2915        '''
2916        npks = self.Background[1].get('nPeaks',0)
2917        if peaknum >= npks:
2918            raise Exception('peak {} not found in histogram {}'.format(peaknum,self.name))
2919        flags = 4*[False]
2920        for i,f in enumerate(refflags):
2921            if i>3: break
2922            flags[i] = bool(f)
2923        for i,f in enumerate(flags):
2924            self.Background[1]['peaksList'][peaknum][2*i+1] = f
2925                   
2926    @property
2927    def id(self):
2928        self.proj.update_ids()
2929        return self.data['data'][0]['hId']
2930
2931    @id.setter
2932    def id(self, val):
2933        self.data['data'][0]['hId'] = val
2934
2935    def fit_fixed_points(self):
2936        """Attempts to apply a background fit to the fixed points currently specified."""
2937        def SetInstParms(Inst):
2938            dataType = Inst['Type'][0]
2939            insVary = []
2940            insNames = []
2941            insVals = []
2942            for parm in Inst:
2943                insNames.append(parm)
2944                insVals.append(Inst[parm][1])
2945                if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
2946                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
2947                        Inst[parm][2] = False
2948            instDict = dict(zip(insNames, insVals))
2949            instDict['X'] = max(instDict['X'], 0.01)
2950            instDict['Y'] = max(instDict['Y'], 0.01)
2951            if 'SH/L' in instDict:
2952                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
2953            return dataType, instDict, insVary
2954
2955        bgrnd = self.data['Background']
2956
2957        # Need our fixed points in order
2958        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
2959        X = [x for x, y in bgrnd[1]['FixedPoints']]
2960        Y = [y for x, y in bgrnd[1]['FixedPoints']]
2961
2962        limits = self.data['Limits'][1]
2963        if X[0] > limits[0]:
2964            X = [limits[0]] + X
2965            Y = [Y[0]] + Y
2966        if X[-1] < limits[1]:
2967            X += [limits[1]]
2968            Y += [Y[-1]]
2969
2970        # Some simple lookups
2971        controls = self.proj['Controls']['data']
2972        inst, inst2 = self.data['Instrument Parameters']
2973        pwddata = self.data['data'][1]
2974
2975        # Construct the data for background fitting
2976        xBeg = np.searchsorted(pwddata[0], limits[0])
2977        xFin = np.searchsorted(pwddata[0], limits[1])
2978        xdata = pwddata[0][xBeg:xFin]
2979        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
2980
2981        W = [1]*len(xdata)
2982        Z = [0]*len(xdata)
2983
2984        dataType, insDict, insVary = SetInstParms(inst)
2985        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
2986
2987        # Do the fit
2988        data = np.array([xdata, ydata, W, Z, Z, Z])
2989        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
2990                        prevVaryList=bakVary, controls=controls)
2991
2992        # Post-fit
2993        parmDict = {}
2994        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
2995        parmDict.update(bakDict)
2996        parmDict.update(insDict)
2997        pwddata[3][xBeg:xFin] *= 0
2998        pwddata[5][xBeg:xFin] *= 0
2999        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
3000
3001        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
3002        # TODO update background
3003        self.proj.save()
3004
3005    def getdata(self,datatype):
3006        '''Provides access to the histogram data of the selected data type
3007
3008        :param str datatype: must be one of the following values (case is ignored)
3009       
3010           * 'X': the 2theta or TOF values for the pattern
3011           * 'Yobs': the observed intensity values
3012           * 'Yweight': the weights for each data point (1/sigma**2)
3013           * 'Ycalc': the computed intensity values
3014           * 'Background': the computed background values
3015           * 'Residual': the difference between Yobs and Ycalc (obs-calc)
3016
3017        :returns: an numpy MaskedArray with data values of the requested type
3018       
3019        '''
3020        enums = ['x', 'yobs', 'yweight', 'ycalc', 'background', 'residual']
3021        if datatype.lower() not in enums:
3022            raise G2ScriptException("Invalid datatype = "+datatype+" must be one of "+str(enums))
3023        return self.data['data'][1][enums.index(datatype.lower())]
3024       
3025    def y_calc(self):
3026        return self.data['data'][1][3]
3027
3028    def Export(self,fileroot,extension):
3029        '''Write the histogram into a file. The path is specified by fileroot and
3030        extension.
3031       
3032        :param str fileroot: name of the file, optionally with a path (extension is
3033           ignored)
3034        :param str extension: includes '.', must match an extension in global
3035           exportersByExtension['powder'] or a Exception is raised.
3036        :returns: name of file that was written
3037        '''
3038        if extension not in exportersByExtension.get('powder',[]):
3039            raise G2ScriptException('No Writer for file type = "'+extension+'"')
3040        fil = os.path.abspath(os.path.splitext(fileroot)[0]+extension)
3041        obj = exportersByExtension['powder'][extension]
3042        obj.SetFromArray(hist=self.data,histname=self.name)
3043        obj.Writer(self.name,fil)
3044        return fil
3045   
3046    def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
3047        try:
3048            import matplotlib.pyplot as plt
3049        except ImportError:
3050            print('Unable to import matplotlib, skipping plot')
3051            return
3052        data = self.data['data'][1]
3053        if Yobs:
3054            plt.plot(data[0], data[1], label='Yobs')
3055        if Ycalc:
3056            plt.plot(data[0], data[3], label='Ycalc')
3057        if Background:
3058            plt.plot(data[0], data[4], label='Background')
3059        if Residual:
3060            plt.plot(data[0], data[5], label="Residual")
3061
3062    def get_wR(self):
3063        """returns the overall weighted profile R factor for a histogram
3064       
3065        :returns: a wR value as a percentage or None if not defined
3066        """
3067        return self['data'][0].get('wR')
3068
3069    def set_refinements(self, refs):
3070        """Sets the histogram refinement parameter 'key' to the specification 'value'.
3071
3072        :param dict refs: A dictionary of the parameters to be set. See the
3073                          :ref:`Histogram_parameters_table` table for a description of
3074                          what these dictionaries should be.
3075
3076        :returns: None
3077
3078        """
3079        do_fit_fixed_points = False
3080        for key, value in refs.items():
3081            if key == 'Limits':
3082                old_limits = self.data['Limits'][1]
3083                new_limits = value
3084                if isinstance(new_limits, dict):
3085                    if 'low' in new_limits:
3086                        old_limits[0] = new_limits['low']
3087                    if 'high' in new_limits:
3088                        old_limits[1] = new_limits['high']
3089                else:
3090                    old_limits[0], old_limits[1] = new_limits
3091            elif key == 'Sample Parameters':
3092                sample = self.data['Sample Parameters']
3093                for sparam in value:
3094                    if sparam not in sample:
3095                        raise ValueError("Unknown refinement parameter, "
3096                                         + str(sparam))
3097                    sample[sparam][1] = True
3098            elif key == 'Background':
3099                bkg, peaks = self.data['Background']
3100
3101                # If True or False, just set the refine parameter
3102                if value in (True, False):
3103                    bkg[1] = value
3104                    return
3105
3106                if 'type' in value:
3107                    bkg[0] = value['type']
3108                if 'refine' in value:
3109                    bkg[1] = value['refine']
3110                if 'no. coeffs' in value:
3111                    cur_coeffs = bkg[2]
3112                    n_coeffs = value['no. coeffs']
3113                    if n_coeffs > cur_coeffs:
3114                        for x in range(n_coeffs - cur_coeffs):
3115                            bkg.append(0.0)
3116                    else:
3117                        for _ in range(cur_coeffs - n_coeffs):
3118                            bkg.pop()
3119                    bkg[2] = n_coeffs
3120                if 'coeffs' in value:
3121                    bkg[3:] = value['coeffs']
3122                if 'FixedPoints' in value:
3123                    peaks['FixedPoints'] = [(float(a), float(b))
3124                                            for a, b in value['FixedPoints']]
3125                if value.get('fit fixed points', False):
3126                    do_fit_fixed_points = True
3127                if 'peaks' in value:
3128                    for i,flags in enumerate(value['peaks']):
3129                        self.ref_back_peak(i,flags)
3130
3131            elif key == 'Instrument Parameters':
3132                instrument, secondary = self.data['Instrument Parameters']
3133                for iparam in value:
3134                    try:
3135                        instrument[iparam][2] = True
3136                    except IndexError:
3137                        raise ValueError("Invalid key:", iparam)
3138            else:
3139                raise ValueError("Unknown key:", key)
3140        # Fit fixed points after the fact - ensure they are after fixed points
3141        # are added
3142        if do_fit_fixed_points:
3143            # Background won't be fit if refinement flag not set
3144            orig = self.data['Background'][0][1]
3145            self.data['Background'][0][1] = True
3146            self.fit_fixed_points()
3147            # Restore the previous value
3148            self.data['Background'][0][1] = orig
3149
3150    def clear_refinements(self, refs):
3151        """Clears the refinement parameter 'key' and its associated value.
3152
3153        :param dict refs: A dictionary of parameters to clear.
3154          See the :ref:`Histogram_parameters_table` table for what can be specified.
3155        """
3156        for key, value in refs.items():
3157            if key == 'Limits':
3158                old_limits, cur_limits = self.data['Limits']
3159                cur_limits[0], cur_limits[1] = old_limits
3160            elif key == 'Sample Parameters':
3161                sample = self.data['Sample Parameters']
3162                for sparam in value:
3163                    sample[sparam][1] = False
3164            elif key == 'Background':
3165                bkg, peaks = self.data['Background']
3166
3167                # If True or False, just set the refine parameter
3168                if value in (True, False):
3169                    bkg[1] = False
3170                    return
3171
3172                bkg[1] = False
3173                if 'FixedPoints' in value:
3174                    if 'FixedPoints' in peaks:
3175                        del peaks['FixedPoints']
3176                if 'peaks' in value:
3177                    for i in range(len(self.Background[1]['peaksList'])):
3178                        self.ref_back_peak(i,[])
3179            elif key == 'Instrument Parameters':
3180                instrument, secondary = self.data['Instrument Parameters']
3181                for iparam in value:
3182                    instrument[iparam][2] = False
3183            else:
3184                raise ValueError("Unknown key:", key)
3185
3186    def add_peak(self,area,dspace=None,Q=None,ttheta=None):
3187        '''Adds a single peak to the peak list
3188        :param float area: peak area
3189        :param float dspace: peak position as d-space (A)
3190        :param float Q: peak position as Q (A-1)
3191        :param float ttheta: peak position as 2Theta (deg)
3192
3193        Note: only one of the parameters dspace, Q or ttheta may be specified
3194        '''
3195        import GSASIIlattice as G2lat
3196        import GSASIImath as G2mth
3197        if (not dspace) + (not Q) + (not ttheta) != 2:
3198            print('add_peak error: too many or no peak position(s) specified')
3199            return
3200        pos = ttheta
3201        Parms,Parms2 = self.data['Instrument Parameters']
3202        if Q:
3203            pos = G2lat.Dsp2pos(Parms,2.*np.pi/Q)
3204        elif dspace:
3205            pos = G2lat.Dsp2pos(Parms,dspace)
3206        peaks = self.data['Peak List']
3207        peaks['sigDict'] = {}        #no longer valid
3208        peaks['peaks'].append(G2mth.setPeakparms(Parms,Parms2,pos,area))
3209
3210    def set_peakFlags(self,peaklist=None,area=None,pos=None,sig=None,gam=None):
3211        '''Set refinement flags for peaks
3212       
3213        :param list peaklist: a list of peaks to change flags. If None (default), changes
3214          are made to all peaks.
3215        :param bool area: Sets or clears the refinement flag for the peak area value.
3216          If None (the default), no change is made.
3217        :param bool pos: Sets or clears the refinement flag for the peak position value.
3218          If None (the default), no change is made.
3219        :param bool sig: Sets or clears the refinement flag for the peak sig (Gaussian width) value.
3220          If None (the default), no change is made.
3221        :param bool gam: Sets or clears the refinement flag for the peak sig (Lorentzian width) value.
3222          If None (the default), no change is made.
3223         
3224        Note that when peaks are first created the area flag is on and the other flags are
3225        initially off.
3226
3227        Example::
3228       
3229           set_peakFlags(sig=False,gam=True)
3230
3231        causes the sig refinement flag to be cleared and the gam flag to be set, in both cases for
3232        all peaks. The position and area flags are not changed from their previous values.
3233        '''
3234        peaks = self.data['Peak List']
3235        if peaklist is None:
3236            peaklist = range(len(peaks['peaks']))
3237        for i in peaklist:
3238            for var,j in [(area,3),(pos,1),(sig,5),(gam,7)]:
3239                if var is not None:
3240                    peaks['peaks'][i][j] = var
3241           
3242    def refine_peaks(self):
3243        '''Causes a refinement of peak position, background and instrument parameters
3244        '''
3245        import GSASIIpwd as G2pwd
3246        controls = self.proj.data.get('Controls',{})
3247        controls = controls.get('data',
3248                            {'deriv type':'analytic','min dM/M':0.001,}     #fill in defaults if needed
3249                            )
3250        peaks = self.data['Peak List']
3251        Parms,Parms2 = self.data['Instrument Parameters']
3252        background = self.data['Background']
3253        limits = self.data['Limits'][1]
3254        bxye = np.zeros(len(self.data['data'][1][1]))
3255        peaks['sigDict'] = G2pwd.DoPeakFit('LSQ',peaks['peaks'],background,limits,
3256                                           Parms,Parms2,self.data['data'][1],bxye,[],
3257                                           False,controls,None)[0]
3258
3259    def SaveProfile(self,filename):
3260        '''Writes a GSAS-II (new style) .instprm file
3261        '''
3262        data,Parms2 = self.data['Instrument Parameters']
3263        filename = os.path.splitext(filename)[0]+'.instprm'         # make sure extension is .instprm
3264        File = open(filename,'w')
3265        File.write("#GSAS-II instrument parameter file; do not add/delete items!\n")
3266        for item in data:
3267            File.write(item+':'+str(data[item][1])+'\n')
3268        File.close()
3269        print ('Instrument parameters saved to: '+filename)
3270
3271    def LoadProfile(self,filename,bank=0):
3272        '''Reads a GSAS-II (new style) .instprm file and overwrites the current
3273        parameters
3274
3275        :param str filename: instrument parameter file name, extension ignored if not
3276          .instprm
3277        :param int bank: bank number to read, defaults to zero
3278        '''
3279        filename = os.path.splitext(filename)[0]+'.instprm'         # make sure extension is .instprm
3280        File = open(filename,'r')
3281        S = File.readline()
3282        newItems = []
3283        newVals = []
3284        Found = False
3285        while S:
3286            if S[0] == '#':
3287                if Found:
3288                    break
3289                if 'Bank' in S:
3290                    if bank == int(S.split(':')[0].split()[1]):
3291                        S = File.readline()
3292                        continue
3293                    else:
3294                        S = File.readline()
3295                        while S and '#Bank' not in S:
3296                            S = File.readline()
3297                        continue
3298                else:   #a non #Bank file
3299                    S = File.readline()
3300                    continue
3301            Found = True
3302            [item,val] = S[:-1].split(':')
3303            newItems.append(item)
3304            try:
3305                newVals.append(float(val))
3306            except ValueError:
3307                newVals.append(val)                       
3308            S = File.readline()               
3309        File.close()
3310        LoadG2fil()
3311        self.data['Instrument Parameters'][0] = G2fil.makeInstDict(newItems,newVals,len(newVals)*[False,])
3312
3313       
3314class G2Phase(G2ObjectWrapper):
3315    """A wrapper object around a given phase.
3316
3317    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
3318    """
3319    def __init__(self, data, name, proj):
3320        self.data = data
3321        self.name = name
3322        self.proj = proj
3323
3324    @staticmethod
3325    def is_valid_refinement_key(key):
3326        valid_keys = ["Cell", "Atoms", "LeBail"]
3327        return key in valid_keys
3328
3329    @staticmethod
3330    def is_valid_HAP_refinement_key(key):
3331        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
3332                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
3333        return key in valid_keys
3334
3335    def atom(self, atomlabel):
3336        """Returns the atom specified by atomlabel, or None if it does not
3337        exist.
3338
3339        :param str atomlabel: The name of the atom (e.g. "O2")
3340        :returns: A :class:`G2AtomRecord` object
3341            representing the atom.
3342        """
3343        # Consult GSASIIobj.py for the meaning of this
3344        cx, ct, cs, cia = self.data['General']['AtomPtrs']
3345        ptrs = [cx, ct, cs, cia]
3346        atoms = self.data['Atoms']
3347        for atom in atoms:
3348            if atom[ct-1] == atomlabel:
3349                return G2AtomRecord(atom, ptrs, self.proj)
3350
3351    def atoms(self):
3352        """Returns a list of atoms present in the current phase.
3353
3354        :returns: A list of :class:`G2AtomRecord` objects.
3355
3356        .. seealso::
3357            :meth:`~G2Phase.atom`
3358            :class:`G2AtomRecord`
3359        """
3360        ptrs = self.data['General']['AtomPtrs']
3361        output = []
3362        atoms = self.data['Atoms']
3363        for atom in atoms:
3364            output.append(G2AtomRecord(atom, ptrs, self.proj))
3365        return output
3366
3367    def histograms(self):
3368        '''Returns a list of histogram names associated with the current phase ordered
3369        as they appear in the tree (see :meth:`G2Project.histograms`).
3370        '''
3371        return [i.name for i in self.proj.histograms() if i.name in self.data.get('Histograms', {})] 
3372
3373    @property
3374    def ranId(self):
3375        return self.data['ranId']
3376
3377    @property
3378    def id(self):
3379        return self.data['pId']
3380
3381    @id.setter
3382    def id(self, val):
3383        self.data['pId'] = val
3384
3385    def get_cell(self):
3386        """Returns a dictionary of the cell parameters, with keys:
3387            'length_a', 'length_b', 'length_c', 'angle_alpha', 'angle_beta', 'angle_gamma', 'volume'
3388
3389        :returns: a dict
3390
3391        .. seealso::
3392           :meth:`~G2Phase.get_cell_and_esd`
3393
3394        """
3395        cell = self.data['General']['Cell']
3396        return {'length_a': cell[1], 'length_b': cell[2], 'length_c': cell[3],
3397                'angle_alpha': cell[4], 'angle_beta': cell[5], 'angle_gamma': cell[6],
3398                'volume': cell[7]}
3399
3400    def get_cell_and_esd(self):
3401        """
3402        Returns a pair of dictionaries, the first representing the unit cell, the second
3403        representing the estimated standard deviations of the unit cell.
3404
3405        :returns: a tuple of two dictionaries
3406
3407        .. seealso::
3408           :meth:`~G2Phase.get_cell`
3409
3410        """
3411        # translated from GSASIIstrIO.ExportBaseclass.GetCell
3412        import GSASIIlattice as G2lat
3413        import GSASIImapvars as G2mv
3414        try:
3415            pfx = str(self.id) + '::'
3416            sgdata = self['General']['SGData']
3417            covDict = self.proj['Covariance']['data']
3418
3419            parmDict = dict(zip(covDict.get('varyList',[]),
3420                                covDict.get('variables',[])))
3421            sigDict = dict(zip(covDict.get('varyList',[]),
3422                               covDict.get('sig',[])))
3423
3424            if covDict.get('covMatrix') is not None:
3425                sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
3426                                                  covDict['varyList'],
3427                                                  parmDict))
3428
3429            A, sigA = G2strIO.cellFill(pfx, sgdata, parmDict, sigDict)
3430            cellSig = G2strIO.getCellEsd(pfx, sgdata, A, self.proj['Covariance']['data'])
3431            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
3432            cellDict, cellSigDict = {}, {}
3433            for i, key in enumerate(['length_a', 'length_b', 'length_c',
3434                                     'angle_alpha', 'angle_beta', 'angle_gamma',
3435                                     'volume']):
3436                cellDict[key] = cellList[i]
3437                cellSigDict[key] = cellSig[i]
3438            return cellDict, cellSigDict
3439        except KeyError:
3440            cell = self.get_cell()
3441            return cell, {key: 0.0 for key in cell}
3442
3443    def export_CIF(self, outputname, quickmode=True):
3444        """Write this phase to a .cif file named outputname
3445
3446        :param str outputname: The name of the .cif file to write to
3447        :param bool quickmode: Currently ignored. Carryover from exports.G2export_CIF"""
3448        # This code is all taken from exports/G2export_CIF.py
3449        # Functions copied have the same names
3450        import GSASIImath as G2mth
3451        import GSASIImapvars as G2mv
3452        from exports import G2export_CIF as cif
3453
3454        CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
3455        CIFname = os.path.splitext(self.proj.filename)[0]
3456        CIFname = os.path.split(CIFname)[1]
3457        CIFname = ''.join([c if ord(c) < 128 else ''
3458                           for c in CIFname.replace(' ', '_')])
3459        try:
3460            author = self.proj['Controls']['data'].get('Author','').strip()
3461        except KeyError:
3462            pass
3463        oneblock = True
3464
3465        covDict = self.proj['Covariance']['data']
3466        parmDict = dict(zip(covDict.get('varyList',[]),
3467                            covDict.get('variables',[])))
3468        sigDict = dict(zip(covDict.get('varyList',[]),
3469                           covDict.get('sig',[])))
3470
3471        if covDict.get('covMatrix') is not None:
3472            sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
3473                                              covDict['varyList'],
3474                                              parmDict))
3475
3476        with open(outputname, 'w') as fp:
3477            fp.write(' \n' + 70*'#' + '\n')
3478            cif.WriteCIFitem(fp, 'data_' + CIFname)
3479            # from exports.G2export_CIF.WritePhaseInfo
3480            cif.WriteCIFitem(fp, '\n# phase info for '+str(self.name) + ' follows')
3481            cif.WriteCIFitem(fp, '_pd_phase_name', self.name)
3482            # TODO get esds
3483            cellDict = self.get_cell()
3484            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
3485            names = ['length_a','length_b','length_c',
3486                     'angle_alpha','angle_beta ','angle_gamma',
3487                     'volume']
3488            for key, val in cellDict.items():
3489                cif.WriteCIFitem(fp, '_cell_' + key, G2mth.ValEsd(val))
3490
3491            cif.WriteCIFitem(fp, '_symmetry_cell_setting',
3492                         self.data['General']['SGData']['SGSys'])
3493
3494            spacegroup = self.data['General']['SGData']['SpGrp'].strip()
3495            # regularize capitalization and remove trailing H/R
3496            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
3497            cif.WriteCIFitem(fp, '_symmetry_space_group_name_H-M', spacegroup)
3498
3499            # generate symmetry operations including centering and center of symmetry
3500            SymOpList, offsetList, symOpList, G2oprList, G2opcodes = G2spc.AllOps(
3501                self.data['General']['SGData'])
3502            cif.WriteCIFitem(fp, 'loop_\n    _space_group_symop_id\n    _space_group_symop_operation_xyz')
3503            for i, op in enumerate(SymOpList,start=1):
3504                cif.WriteCIFitem(fp, '   {:3d}  {:}'.format(i,op.lower()))
3505
3506            # TODO skipped histograms, exports/G2export_CIF.py:880
3507
3508            # report atom params
3509            if self.data['General']['Type'] in ['nuclear','macromolecular']:        #this needs macromolecular variant, etc!
3510                cif.WriteAtomsNuclear(fp, self.data, self.name, parmDict, sigDict, [])
3511                # self._WriteAtomsNuclear(fp, parmDict, sigDict)
3512            else:
3513                raise G2ScriptException("no export for "+str(self.data['General']['Type'])+" coordinates implemented")
3514            # report cell contents
3515            cif.WriteComposition(fp, self.data, self.name, parmDict)
3516            if not quickmode and self.data['General']['Type'] == 'nuclear':      # report distances and angles
3517                # WriteDistances(fp,self.name,SymOpList,offsetList,symOpList,G2oprList)
3518                raise NotImplementedError("only quickmode currently supported")
3519            if 'Map' in self.data['General'] and 'minmax' in self.data['General']['Map']:
3520                cif.WriteCIFitem(fp,'\n# Difference density results')
3521                MinMax = self.data['General']['Map']['minmax']
3522                cif.WriteCIFitem(fp,'_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009))
3523                cif.WriteCIFitem(fp,'_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009))
3524
3525
3526    def set_refinements(self, refs):
3527        """Sets the phase refinement parameter 'key' to the specification 'value'
3528
3529        :param dict refs: A dictionary of the parameters to be set. See the
3530                          :ref:`Phase_parameters_table` table for a description of
3531                          this dictionary.
3532
3533        :returns: None"""
3534        for key, value in refs.items():
3535            if key == "Cell":
3536                self.data['General']['Cell'][0] = value
3537
3538            elif key == "Atoms":
3539                for atomlabel, atomrefinement in value.items():
3540                    if atomlabel == 'all':
3541                        for atom in self.atoms():
3542                            atom.refinement_flags = atomrefinement
3543                    else:
3544                        atom = self.atom(atomlabel)
3545                        if atom is None:
3546                            raise ValueError("No such atom: " + atomlabel)
3547                        atom.refinement_flags = atomrefinement
3548
3549            elif key == "LeBail":
3550                hists = self.data['Histograms']
3551                for hname, hoptions in hists.items():
3552                    if 'LeBail' not in hoptions:
3553                        hoptions['newLeBail'] = bool(True)
3554                    hoptions['LeBail'] = bool(value)
3555            else:
3556                raise ValueError("Unknown key:", key)
3557
3558    def clear_refinements(self, refs):
3559        """Clears a given set of parameters.
3560
3561        :param dict refs: The parameters to clear.
3562          See the :ref:`Phase_parameters_table` table for what can be specified.
3563        """
3564        for key, value in refs.items():
3565            if key == "Cell":
3566                self.data['General']['Cell'][0] = False
3567            elif key == "Atoms":
3568                cx, ct, cs, cia = self.data['General']['AtomPtrs']
3569
3570                for atomlabel in value:
3571                    atom = self.atom(atomlabel)
3572                    # Set refinement to none
3573                    atom.refinement_flags = ' '
3574            elif key == "LeBail":
3575                hists = self.data['Histograms']
3576                for hname, hoptions in hists.items():
3577                    if 'LeBail' not in hoptions:
3578                        hoptions['newLeBail'] = True
3579                    hoptions['LeBail'] = False
3580            else:
3581                raise ValueError("Unknown key:", key)
3582
3583    def set_HAP_refinements(self, refs, histograms='all'):
3584        """Sets the given HAP refinement parameters between the current phase and
3585        the specified histograms.
3586
3587        :param dict refs: A dictionary of the parameters to be set. See
3588                          the :ref:`HAP_parameters_table` table for a description of this
3589                          dictionary.
3590        :param histograms: Either 'all' (default) or a list of the histograms by index, name
3591            or object. The index number is relative to all histograms in the tree, not to
3592            those in the phase.
3593            Histograms not associated with the current phase will be ignored.
3594            whose HAP parameters will be set with this phase. Histogram and phase
3595            must already be associated.
3596        :returns: None
3597        """
3598        if not self.data.get('Histograms',[]):
3599            print("Error likely: Phase {} has no linked histograms".format(self.name))
3600            return
3601           
3602        if histograms == 'all':
3603            histograms = self.data['Histograms'].keys()
3604        else:
3605            histograms = [self._decodeHist(h) for h in histograms
3606                          if self._decodeHist(h) in self.data['Histograms']]   
3607        if not histograms:
3608            print("Skipping HAP set for phase {}, no selected histograms".format(self.name))
3609            return
3610        for key, val in refs.items():
3611            if key == 'Babinet':
3612                try:
3613                    sets = list(val)
3614                except ValueError:
3615                    sets = ['BabA', 'BabU']
3616                for param in sets:
3617                    if param not in ['BabA', 'BabU']:
3618                        raise ValueError("Not sure what to do with" + param)
3619                    for h in histograms:
3620                        self.data['Histograms'][h]['Babinet'][param][1] = True
3621            elif key == 'Extinction':
3622                for h in histograms:
3623                    self.data['Histograms'][h]['Extinction'][1] = bool(val)
3624            elif key == 'HStrain':
3625                if isinstance(val,list) or isinstance(val,tuple):
3626                    for h in histograms:
3627                        if len(self.data['Histograms'][h]['HStrain'][1]) != len(val):
3628                            raise Exception('Need {} HStrain terms for phase {} hist {}'
3629                                .format(len(self.data['Histograms'][h]['HStrain'][1]),self.name,h))
3630                        for i,v in enumerate(val):
3631                            self.data['Histograms'][h]['HStrain'][1][i] = bool(v)
3632                else:
3633                    for h in histograms:
3634                        self.data['Histograms'][h]['HStrain'][1] = [bool(val) for p in self.data['Histograms'][h]['HStrain'][1]]
3635            elif key == 'Mustrain':
3636                for h in histograms:
3637                    mustrain = self.data['Histograms'][h]['Mustrain']
3638                    newType = None
3639                    direction = None
3640                    if isinstance(val, strtypes):
3641                        if val in ['isotropic', 'uniaxial', 'generalized']:
3642                            newType = val
3643                        else:
3644                            raise ValueError("Not a Mustrain type: " + val)
3645                    elif isinstance(val, dict):
3646                        newType = val.get('type', None)
3647                        direction = val.get('direction', None)
3648
3649                    if newType:
3650                        mustrain[0] = newType
3651                        if newType == 'isotropic':
3652                            mustrain[2][0] = True == val.get('refine',False)
3653                            mustrain[5] = [False for p in mustrain[4]]
3654                        elif newType == 'uniaxial':
3655                            if 'refine' in val:
3656                                mustrain[2][0] = False
3657                                types = val['refine']
3658                                if isinstance(types, strtypes):
3659                                    types = [types]
3660                                elif isinstance(types, bool):
3661                                    mustrain[2][1] = types
3662                                    mustrain[2][2] = types
3663                                    types = []
3664                                else:
3665                                    raise ValueError("Not sure what to do with: "
3666                                                     + str(types))
3667                            else:
3668                                types = []
3669
3670                            for unitype in types:
3671                                if unitype == 'equatorial':
3672                                    mustrain[2][0] = True
3673                                elif unitype == 'axial':
3674                                    mustrain[2][1] = True
3675                                else:
3676                                    msg = 'Invalid uniaxial mustrain type'
3677                                    raise ValueError(msg + ': ' + unitype)
3678                        else:  # newtype == 'generalized'
3679                            mustrain[2] = [False for p in mustrain[1]]
3680                            if 'refine' in val:
3681                                mustrain[5] = [True == val['refine']]*len(mustrain[5])
3682
3683                    if direction:
3684                        if len(direction) != 3:
3685                            raise ValueError("Expected hkl, found", direction)
3686                        direction = [int(n) for n in direction]
3687                        mustrain[3] = direction
3688            elif key == 'Size':
3689                newSize = None
3690                if 'value' in val:
3691                    newSize = float(val['value'])
3692                for h in histograms:
3693                    size = self.data['Histograms'][h]['Size']
3694                    newType = None
3695                    direction = None
3696                    if isinstance(val, strtypes):
3697                        if val in ['isotropic', 'uniaxial', 'ellipsoidal']:
3698                            newType = val
3699                        else:
3700                            raise ValueError("Not a valid Size type: " + val)
3701                    elif isinstance(val, dict):
3702                        newType = val.get('type', None)
3703                        direction = val.get('direction', None)
3704
3705                    if newType:
3706                        size[0] = newType
3707                        refine = bool(val.get('refine'))
3708                        if newType == 'isotropic' and refine is not None:
3709                            size[2][0] = bool(refine)
3710                            if newSize: size[1][0] = newSize
3711                        elif newType == 'uniaxial' and refine is not None:
3712                            size[2][1] = bool(refine)
3713                            size[2][2] = bool(refine)
3714                            if newSize: size[1][1] = size[1][2] =newSize
3715                        elif newType == 'ellipsoidal' and refine is not None:
3716                            size[5] = [bool(refine) for p in size[5]]
3717                            if newSize: size[4] = [newSize for p in size[4]]
3718
3719                    if direction:
3720                        if len(direction) != 3:
3721                            raise ValueError("Expected hkl, found", direction)
3722                        direction = [int(n) for n in direction]
3723                        size[3] = direction
3724            elif key == 'Pref.Ori.':
3725                for h in histograms:
3726                    self.data['Histograms'][h]['Pref.Ori.'][2] = bool(val)
3727            elif key == 'Show':
3728                for h in histograms:
3729                    self.data['Histograms'][h]['Show'] = bool(val)
3730            elif key == 'Use':
3731                for h in histograms:
3732                    self.data['Histograms'][h]['Use'] = bool(val)
3733            elif key == 'Scale':
3734                for h in histograms:
3735                    self.data['Histograms'][h]['Scale'][1] = bool(val)
3736            else:
3737                print(u'Unknown HAP key: '+key)
3738
3739    def clear_HAP_refinements(self, refs, histograms='all'):
3740        """Clears the given HAP refinement parameters between this phase and
3741        the given histograms.
3742
3743        :param dict refs: A dictionary of the parameters to be cleared.
3744            See the the :ref:`HAP_parameters_table` table for what can be specified.
3745        :param histograms: Either 'all' (default) or a list of the histograms by index, name
3746            or object.
3747            The index number is relative to all histograms in the tree, not to
3748            those in the phase.
3749            Histograms not associated with the current phase will be ignored.
3750            whose HAP parameters will be set with this phase. Histogram and phase
3751            must already be associated
3752        :returns: None
3753        """
3754        if histograms == 'all':
3755            histograms = self.data['Histograms'].keys()
3756        else:
3757            histograms = [self._decodeHist(h) for h in histograms
3758                          if self._decodeHist(h) in self.data['Histograms']]   
3759
3760        for key, val in refs.items():
3761            for h in histograms:
3762                if key == 'Babinet':
3763                    try:
3764                        sets = list(val)
3765                    except ValueError:
3766                        sets = ['BabA', 'BabU']
3767                    for param in sets:
3768                        if param not in ['BabA', 'BabU']:
3769                            raise ValueError("Not sure what to do with" + param)
3770                        for h in histograms:
3771                            self.data['Histograms'][h]['Babinet'][param][1] = False
3772                elif key == 'Extinction':
3773                    for h in histograms:
3774                        self.data['Histograms'][h]['Extinction'][1] = False
3775                elif key == 'HStrain':
3776                    for h in histograms:
3777                        self.data['Histograms'][h]['HStrain'][1] = [False for p in self.data['Histograms'][h]['HStrain'][1]]
3778                elif key == 'Mustrain':
3779                    for h in histograms:
3780                        mustrain = self.data['Histograms'][h]['Mustrain']
3781                        mustrain[2] = [False for p in mustrain[2]]
3782                        mustrain[5] = [False for p in mustrain[4]]
3783                elif key == 'Pref.Ori.':
3784                    for h in histograms:
3785                        self.data['Histograms'][h]['Pref.Ori.'][2] = False
3786                elif key == 'Show':
3787                    for h in histograms:
3788                        self.data['Histograms'][h]['Show'] = False
3789                elif key == 'Size':
3790                    for h in histograms:
3791                        size = self.data['Histograms'][h]['Size']
3792                        size[2] = [False for p in size[2]]
3793                        size[5] = [False for p in size[5]]
3794                elif key == 'Use':
3795                    for h in histograms:
3796                        self.data['Histograms'][h]['Use'] = False
3797                elif key == 'Scale':
3798                    for h in histograms:
3799                        self.data['Histograms'][h]['Scale'][1] = False
3800                else:
3801                    print(u'Unknown HAP key: '+key)
3802
3803    def _decodeHist(self,hist):
3804        '''Convert a histogram reference to a histogram name string
3805        '''
3806        if isinstance(hist, G2PwdrData):
3807            return hist.name
3808        elif hist in self.data['Histograms']:
3809            return hist
3810        elif type(hist) is int:
3811            return self.proj.histograms()[hist].name
3812        else:
3813            raise G2ScriptException("Invalid histogram reference: "+str(hist))
3814       
3815    def getHAPvalues(self, histname):
3816        """Returns a dict with HAP values for the selected histogram
3817
3818        :param histogram: is a histogram object (:class:`G2PwdrData`) or
3819            a histogram name or the index number of the histogram.
3820            The index number is relative to all histograms in the tree, not to
3821            those in the phase.
3822        :returns: HAP value dict
3823        """
3824        return self.data['Histograms'][self._decodeHist(histname)]
3825
3826    def copyHAPvalues(self, sourcehist, targethistlist='all', skip=[], use=None):
3827        """Copies HAP parameters for one histogram to a list of other histograms.
3828        Use skip or use to select specific entries to be copied or not used.
3829
3830        :param sourcehist: is a histogram object (:class:`G2PwdrData`) or
3831            a histogram name or the index number of the histogram to copy
3832            parameters from.
3833            The index number is relative to all histograms in the tree, not to
3834            those in the phase.
3835        :param list targethistlist: a list of histograms where each item in the
3836            list can be a histogram object (:class:`G2PwdrData`),
3837            a histogram name or the index number of the histogram.
3838            If the string 'all' (default), then all histograms in the phase
3839            are used.
3840        :param list skip: items in the HAP dict that should not be
3841            copied. The default is an empty list, which causes all items
3842            to be copied. To see a list of items in the dict, use
3843            :meth:`getHAPvalues` or use an invalid item, such as '?'.
3844        :param list use: specifies the items in the HAP dict should be
3845            copied. The default is None, which causes all items
3846            to be copied.
3847
3848        examples::
3849
3850            ph0.copyHAPvalues(0,[1,2,3])
3851            ph0.copyHAPvalues(0,use=['HStrain','Size'])
3852
3853        The first example copies all HAP parameters from the first histogram to
3854        the second, third and fourth histograms (as listed in the project tree).
3855        The second example copies only the 'HStrain' (Dij parameters and
3856        refinement flags) and the 'Size' (crystallite size settings, parameters
3857        and refinement flags) from the first histogram to all histograms.
3858        """
3859        sourcehist = self._decodeHist(sourcehist)
3860        if targethistlist == 'all':
3861            targethistlist = self.histograms()
3862       
3863        copydict = copy.deepcopy(self.data['Histograms'][sourcehist])
3864        for item in skip:
3865            if item in list(copydict.keys()):
3866                del copydict[item]
3867            else:
3868                print('items in HAP dict are: {}'.format(
3869                    list(self.data['Histograms'][sourcehist])))
3870                raise Exception('HAP skip list entry {} invalid'.format(item))
3871        if use:
3872            for item in list(copydict.keys()):
3873                if item not in use:
3874                    del copydict[item]
3875
3876        print('Copying item(s) {} from histogram {}'.format(list(copydict.keys()),sourcehist))
3877        print(' to histogram(s) {}'.format([self._decodeHist(h) for h in targethistlist]))
3878        for h in targethistlist:
3879            h = self._decodeHist(h)
3880            if h not in self.data['Histograms']:
3881                print('Unexpected Warning: histogram {} not in phase {}'.format(h,self.name))
3882                continue
3883            self.data['Histograms'][h].update(copy.deepcopy(copydict))
3884
3885    def setHAPvalues(self, HAPdict, targethistlist='all', skip=[], use=None):
3886        """Copies HAP parameters for one histogram to a list of other histograms.
3887        Use skip or use to select specific entries to be copied or not used.
3888        Note that ``HStrain`` and sometimes ``Mustrain`` values can be specific to
3889        a Laue class and should be copied with care between phases of different
3890        symmetry. A "sanity check" on the number of Dij terms is made if ``HStrain``
3891        values are copied.
3892
3893        :param dict HAPdict: is a dict returned by :meth:`getHAPvalues` containing
3894            HAP parameters.
3895        :param list targethistlist: a list of histograms where each item in the
3896            list can be a histogram object (:class:`G2PwdrData`),
3897            a histogram name or the index number of the histogram.
3898            The index number is relative to all histograms in the tree, not to
3899            those in the phase.
3900            If the string 'all' (default), then all histograms in the phase
3901            are used.
3902        :param list skip: items in the HAP dict that should not be
3903            copied. The default is an empty list, which causes all items
3904            to be copied. To see a list of items in the dict, use
3905            :meth:`getHAPvalues` or use an invalid item, such as '?'.
3906        :param list use: specifies the items in the HAP dict should be
3907            copied. The default is None, which causes all items
3908            to be copied.
3909
3910        example::
3911
3912            HAPdict = ph0.getHAPvalues(0)
3913            ph1.setHAPvalues(HAPdict,use=['HStrain','Size'])
3914
3915        This copies the Dij (hydrostatic strain) HAP parameters and the
3916        crystallite size broadening terms from the first histogram in
3917        phase ``ph0`` to all histograms in phase ``ph1``.
3918        """
3919        if targethistlist == 'all':
3920            targethistlist = self.histograms()
3921        copydict = copy.deepcopy(HAPdict)
3922        for item in skip:
3923            if item in list(copydict.keys()):
3924                del copydict[item]
3925            else:
3926                print('items in HAP dict are: {}'.format(
3927                    list(self.data['Histograms'][sourcehist])))
3928                raise Exception('HAP skip list entry {} invalid'.format(item))
3929        if use:
3930            for item in list(copydict.keys()):
3931                if item not in use:
3932                    del copydict[item]
3933
3934        first = True
3935        for h in targethistlist:
3936            h = self._decodeHist(h)
3937            if h not in self.data['Histograms']:
3938                print('Warning: histogram {} not in phase {}'.format(h,self.name))
3939                continue
3940            if first:
3941                first = False
3942                if 'HStrain' in self.data['Histograms'][h] and 'HStrain' in copydict:
3943                    if len(copydict['HStrain'][0]) != len(self.data['Histograms'][h]['HStrain'][0]):
3944                        print('HStrain has differing numbers of terms. Input: {}, phase {}: {}'.
3945                                  format(len(copydict['HStrain'][0]),
3946                                        self.name,len(self.data['Histograms'][h]['HStrain'][0])))
3947                        raise Exception('HStrain has differing numbers of terms.')
3948            self.data['Histograms'][h].update(copy.deepcopy(copydict))           
3949        print('Copied item(s) {} from dict'.format(list(copydict.keys())))
3950        print(' to histogram(s) {}'.format([self._decodeHist(h) for h in targethistlist]))
3951
3952class G2SeqRefRes(G2ObjectWrapper):
3953    '''Wrapper for a Sequential Refinement Results tree entry, containing the
3954    results for a refinement
3955
3956    As an example::
3957
3958        from __future__ import division, print_function
3959        import os,sys
3960        sys.path.insert(0,'/Users/toby/software/G2/GSASII')
3961        PathWrap = lambda fil: os.path.join('/Users/toby/Scratch/SeqTut2019Mar',fil)
3962        import GSASIIscriptable as G2sc
3963        gpx = G2sc.G2Project(PathWrap('scr4.gpx'))
3964        seq = gpx.seqref()
3965        lbl = ('a','b','c','alpha','beta','gamma','Volume')
3966        for j,h in enumerate(seq.histograms()):
3967            cell,cellU,uniq = seq.get_cell_and_esd(1,h)
3968            print(h)
3969            print([cell[i] for i in list(uniq)+[6]])
3970            print([cellU[i] for i in list(uniq)+[6]])
3971            print('')
3972        print('printed',[lbl[i] for i in list(uniq)+[6]])
3973
3974    .. seealso::
3975        :meth:`G2Project.seqref`
3976    '''
3977    def __init__(self, data, proj):
3978        self.data = data
3979        self.proj = proj
3980        self.newCellDict = {} # dict with recp. cell tensor & Dij name
3981        # newAtomDict = {} # dict with atom positions; relative & absolute
3982        for name in self.data['histNames']:
3983            self.newCellDict.update(self.data[name].get('newCellDict',{}))
3984            #newAtomDict.update(self.data[name].get('newAtomDict',{}) # dict with atom positions; relative & absolute
3985        #ESDlookup = {self.newCellDict[item][0]:item for item in self.newCellDict}
3986        #Dlookup = {item:self.newCellDict[item][0] for item in self.newCellDict}
3987        # Possible error: the next might need to be data[histNames[0]]['varyList']
3988        #atomLookup = {newAtomDict[item][0]:item for item in newAtomDict if item in self.data['varyList']}
3989        #Dlookup.update({atomLookup[parm]:parm for parm in atomLookup}
3990        #ESDlookup.update({parm:atomLookup[parm] for parm in atomLookup})
3991
3992#    @property
3993    def histograms(self):
3994        '''returns a list of histograms in the squential fit
3995        '''
3996        return self.data['histNames']
3997
3998    def get_cell_and_esd(self,phase,hist):
3999        '''Returns a vector of cell lengths and esd values
4000
4001        :param phase: A phase, which may be specified as a phase object
4002          (see :class:`G2Phase`), the phase name (str) or the index number (int)
4003          of the phase in the project, numbered starting from 0.
4004        :param hist: Specify a histogram or using the histogram name (str)
4005          or the index number (int) of the histogram in the sequential
4006          refinement (not the project), numbered as in in the project tree
4007          starting from 0.
4008        :returns: cell,cellESD,uniqCellIndx where cell (list)
4009          with the unit cell parameters (a,b,c,alpha,beta,gamma,Volume);
4010          cellESD are the standard uncertainties on the 7 unit cell
4011          parameters; and uniqCellIndx is a tuple with indicies for the
4012          unique (non-symmetry determined) unit parameters (e.g.
4013          [0,2] for a,c in a tetragonal cell)
4014        '''
4015        def striphist(var,insChar=''):
4016            'strip a histogram number from a var name'
4017            sv = var.split(':')
4018            if len(sv) <= 1: return var
4019            if sv[1]:
4020                sv[1] = insChar
4021            return ':'.join(sv)
4022        import GSASIIlattice as G2lat
4023        import GSASIIstrIO as G2stIO
4024       
4025        uniqCellLookup = [
4026        [['m3','m3m'],(0,)],
4027        [['3R','3mR'],(0,3)],
4028        [['3','3m1','31m','6/m','6/mmm','4/m','4/mmm'],(0,2)],
4029        [['mmm'],(0,1,2)],
4030        [['2/m'+'a'],(0,1,2,3)],
4031        [['2/m'+'b'],(0,1,2,4)],
4032        [['2/m'+'c'],(0,1,2,5)],
4033        [['-1'],(0,1,2,3,4,5)],
4034        ]
4035
4036        seqData,histData = self.RefData(hist)
4037        hId = histData['data'][0]['hId']
4038        phasedict = self.proj.phase(phase).data
4039        pId = phasedict['pId']
4040        pfx = str(pId)+'::' # prefix for A values from phase
4041        phfx = '%d:%d:'%(pId,hId) # Dij prefix
4042        # get unit cell & symmetry for phase
4043        RecpCellTerms = G2lat.cell2A(phasedict['General']['Cell'][1:7])
4044        zeroDict = {pfx+'A'+str(i):0.0 for i in range(6)}
4045        SGdata = phasedict['General']['SGData']
4046        # determine the cell items not defined by symmetry
4047        laue = SGdata['SGLaue'][:]
4048        if laue == '2/m':
4049            laue += SGdata['SGUniq']
4050        for symlist,celllist in uniqCellLookup:
4051            if laue in symlist:
4052                uniqCellIndx = celllist
4053                break
4054        else: # should not happen
4055            uniqCellIndx = list(range(6))
4056        initialCell = {}
4057        for i in uniqCellIndx:
4058            initialCell[str(pId)+'::A'+str(i)] =  RecpCellTerms[i]
4059
4060        esdLookUp = {}
4061        dLookup = {}
4062        # note that varyList keys are p:h:Dij while newCellDict keys are p::Dij
4063        for nKey in seqData['newCellDict']:
4064            p = nKey.find('::')+1
4065            vKey = nKey[:p] + str(hId) + nKey[p:]
4066            if vKey in seqData['varyList']:
4067                esdLookUp[self.newCellDict[nKey][0]] = nKey
4068                dLookup[nKey] = self.newCellDict[nKey][0]
4069        covData = {'varyList': [dLookup.get(striphist(v),v) for v in seqData['varyList']],
4070                'covMatrix': seqData['covMatrix']}
4071        A = RecpCellTerms[:] # make copy of starting A values
4072
4073        for i,j in enumerate(('D11','D22','D33','D12','D13','D23')):
4074            var = pfx+'A'+str(i)
4075            Dvar = phfx+j
4076            # apply Dij value if non-zero
4077            if Dvar in seqData['parmDict']:
4078                A[i] += seqData['parmDict'][Dvar]
4079            # override with fit result if is Dij varied
4080            try:
4081                A[i] = seqData['newCellDict'][esdLookUp[var]][1] # get refined value
4082            except KeyError:
4083                pass
4084        Albls = [pfx+'A'+str(i) for i in range(6)]
4085        cellDict = dict(zip(Albls,A))
4086
4087       
4088        A,zeros = G2stIO.cellFill(pfx,SGdata,cellDict,zeroDict)
4089        # convert to direct cell
4090        c = G2lat.A2cell(A)
4091        vol = G2lat.calc_V(A)
4092        cE = G2stIO.getCellEsd(pfx,SGdata,A,covData)
4093        return list(c)+[vol],cE,uniqCellIndx
4094   
4095    def get_VaryList(self,hist):
4096        '''Returns a list of the refined variables in the
4097        last refinement cycle for the selected histogram
4098
4099        :param hist: Specify a histogram or using the histogram name (str)
4100          or the index number (int) of the histogram in the sequential
4101          refinement (not the project), numbered starting from 0.
4102        :returns: a list of variables or None if no refinement has been
4103          performed.
4104        '''
4105        try:
4106            seqData,histData = self.RefData(hist)
4107            return seqData['varyList']
4108        except:
4109            return
4110
4111    def get_ParmList(self,hist):
4112        '''Returns a list of all the parameters defined in the
4113        last refinement cycle for the selected histogram
4114
4115        :param hist: Specify a histogram or using the histogram name (str)
4116          or the index number (int) of the histogram in the sequential
4117          refinement (not the project), numbered as in the project tree
4118          starting from 0.
4119        :returns: a list of parameters or None if no refinement has been
4120          performed.
4121        '''
4122        try:
4123            seqData,histData = self.RefData(hist)
4124            return list(seqData['parmDict'].keys())
4125        except:
4126            return
4127
4128    def get_Variable(self,hist,var):
4129        '''Returns the value and standard uncertainty (esd) for a variable
4130        parameters, as defined for the selected histogram
4131        in the last sequential refinement cycle
4132
4133        :param hist: Specify a histogram or using the histogram name (str)
4134          or the index number (int) of the histogram in the sequential
4135          refinement (not the project), numbered as in the project tree
4136          starting from 0.
4137        :param str var: a variable name of form '<p>:<h>:<name>', such as
4138          ':0:Scale'
4139        :returns: (value,esd) if the parameter is refined or
4140          (value, None) if the variable is in a constraint or is not
4141          refined or None if the parameter is not found.
4142        '''
4143        try:
4144            seqData,histData = self.RefData(hist)
4145            val = seqData['parmDict'][var]
4146        except:
4147            return
4148        try:
4149            pos = seqData['varyList'].index(var)
4150            esd = np.sqrt(seqData['covMatrix'][pos,pos])
4151            return (val,esd)
4152        except ValueError:
4153            return (val,None)
4154
4155    def get_Covariance(self,hist,varList):
4156        '''Returns the values and covariance matrix for a series of variable
4157        parameters, as defined for the selected histogram
4158        in the last sequential refinement cycle
4159
4160        :param hist: Specify a histogram or using the histogram name (str)
4161          or the index number (int) of the histogram in the sequential
4162          refinement (not the project), numbered as in the project tree
4163          starting from 0.
4164        :param tuple varList: a list of variable names of form '<p>:<h>:<name>'
4165        :returns: (valueList,CovMatrix) where valueList contains the (n) values
4166          in the same order as varList (also length n) and CovMatrix is a
4167          (n x n) matrix. If any variable name is not found in the varyList
4168          then None is returned.
4169
4170        Use this code, where sig provides standard uncertainties for
4171        parameters and where covArray provides the correlation between
4172        off-diagonal terms::
4173
4174            sig = np.sqrt(np.diag(covMatrix))
4175            xvar = np.outer(sig,np.ones_like(sig))
4176            covArray = np.divide(np.divide(covMatrix,xvar),xvar.T)
4177
4178        '''
4179        try:
4180            seqData,histData = self.RefData(hist)
4181        except:
4182            print('Histogram {} not found in the sequential fit'.format(hist))
4183            return
4184        missing = [i for i in varList if i not in seqData['varyList']]
4185        if missing:
4186            print('Variable(s) {} were not found in the varyList'.format(missing))
4187            return None
4188        vals = [seqData['parmDict'][i] for i in varList]
4189        import GSASIImath as G2mth
4190        cov = G2mth.getVCov(varList,seqData['varyList'],seqData['covMatrix'])
4191        return (vals,cov)
4192   
4193    def RefData(self,hist):
4194        '''Provides access to the output from a particular histogram
4195
4196        :param hist: Specify a histogram or using the histogram name (str)
4197          or the index number (int) of the histogram in the sequential
4198          refinement (not the project), numbered as in the project tree
4199          starting from 0.
4200        :returns: a list of dicts where the first element has sequential
4201          refinement results and the second element has the contents of
4202          the histogram tree items.
4203        '''
4204        try:
4205            hist = self.data['histNames'][hist]
4206        except IndexError:
4207            raise Exception('Histogram #{} is out of range from the Sequential Refinement'
4208                                .format(hist))
4209        except TypeError:
4210            pass
4211        if hist not in self.data['histNames']:
4212            raise Exception('Histogram {} is not included in the Sequential Refinement'
4213                                .format(hist))
4214        return self.data[hist],self.proj.histogram(hist).data
4215
4216class G2PDF(G2ObjectWrapper):
4217    """Wrapper for a PDF tree entry, containing the information needed to
4218    compute a PDF and the S(Q), G(r) etc. after the computation is done.
4219    Note that in a GSASIIscriptable script, instances of G2PDF will be created by
4220    calls to :meth:`G2Project.add_PDF` or :meth:`G2Project.pdf`, not via calls
4221    to :meth:`G2PDF.__init__`.
4222
4223    Example use of :class:`G2PDF`::
4224
4225       gpx.add_PDF('250umSiO2.pdfprm',0)
4226       pdf.set_formula(['Si',1],['O',2])
4227       pdf.set_background('Container',1,-0.21)
4228       for i in range(5):
4229           if pdf.optimize(): break
4230       pdf.calculate()
4231       pdf.export(gpx.filename,'S(Q), pdfGUI')
4232       gpx.save('pdfcalc.gpx')
4233
4234    .. seealso::
4235        :meth:`G2Project.pdf`
4236        :meth:`G2Project.pdfs`
4237    """
4238    def __init__(self, data, name, proj):
4239        self.data = data
4240        self.name = name
4241        self.proj = proj
4242    def set_background(self,btype,histogram,mult=-1.,refine=False):
4243        '''Sets a histogram to be used as the 'Sample Background',
4244        the 'Container' or the 'Container Background.'
4245
4246        :param str btype: Type of background to set, must contain
4247          the string 'samp' for Sample Background', 'cont' and 'back'
4248          for the 'Container Background' or only 'cont' for the
4249          'Container'. Note that capitalization and extra characters
4250          are ignored, so the full strings (such as 'Sample
4251          Background' & 'Container Background') can be used.
4252       
4253        :param histogram: A reference to a histogram,
4254          which can be reference by object, name, or number.
4255        :param float mult: a multiplier for the histogram; defaults
4256          to -1.0
4257        :param bool refine: a flag to enable refinement (only
4258          implemented for 'Sample Background'); defaults to False
4259        '''
4260        if 'samp' in btype.lower():
4261            key = 'Sample Bkg.'
4262        elif 'cont' in btype.lower() and 'back' in btype.lower():
4263            key = 'Container Bkg.'
4264        elif 'cont' in btype.lower():
4265            key = 'Container'
4266        else:
4267            raise Exception('btype = {} is invalid'.format(btype))
4268        self.data['PDF Controls'][key]['Name'] = self.proj.histogram(histogram).name
4269        self.data['PDF Controls'][key]['Mult'] = mult
4270        self.data['PDF Controls'][key]['Refine'] = refine
4271
4272    def set_formula(self,*args):
4273        '''Set the chemical formula for the PDF computation.
4274        Use pdf.set_formula(['Si',1],['O',2]) for SiO2.
4275
4276        :param list item1: The element symbol and number of atoms in formula for first element
4277        :param list item2: The element symbol and number of atoms in formula for second element,...
4278
4279        repeat parameters as needed for all elements in the formula.
4280        '''
4281        powderHist = self.proj.histogram(self.data['PDF Controls']['Sample']['Name'])
4282        inst = powderHist.data['Instrument Parameters'][0]
4283        ElList = self.data['PDF Controls']['ElList']
4284        ElList.clear()
4285        sumVol = 0.
4286        for elem,mult in args:
4287            ElList[elem] = G2elem.GetElInfo(elem,inst)
4288            ElList[elem]['FormulaNo'] = mult
4289            Avol = (4.*np.pi/3.)*ElList[elem]['Drad']**3
4290            sumVol += Avol*ElList[elem]['FormulaNo']
4291        self.data['PDF Controls']['Form Vol'] = max(10.0,sumVol)
4292       
4293    def calculate(self):
4294        '''Compute the PDF using the current parameters. Results are set
4295        in the PDF object arrays (self.data['PDF Controls']['G(R)'] etc.)
4296        '''
4297        data = self.data['PDF Controls']
4298        xydata = {}
4299        for key in 'Sample Bkg.','Container Bkg.','Container','Sample':
4300            name = data[key]['Name'].strip()
4301            if name:
4302                xydata[key] = self.proj.histogram(name).data['data']
4303                if key != 'Sample': continue
4304                limits = self.proj.histogram(name).data['Limits'][1]
4305                inst = self.proj.histogram(name).data['Instrument Parameters'][0]
4306        G2pwd.CalcPDF(data,inst,limits,xydata)
4307        data['I(Q)'] = xydata['IofQ']
4308        data['S(Q)'] = xydata['SofQ']
4309        data['F(Q)'] = xydata['FofQ']
4310        data['G(R)'] = xydata['GofR']
4311
4312    def optimize(self,showFit=True,maxCycles=5):
4313        '''Optimize the low R portion of G(R) to minimize selected
4314        parameters. Note that this updates the parameters in the settings
4315        (self.data['PDF Controls']) but does not update the PDF object
4316        arrays (self.data['PDF Controls']['G(R)'] etc.) with the computed
4317        values, use :meth:`calculate` after a fit to do that.
4318
4319        :param bool showFit: if True (default) the optimized parameters
4320          are shown before and after the fit, as well as the RMS value
4321          in the minimized region.
4322        :param int maxCycles: the maximum number of least-squares cycles;
4323          defaults to 5.
4324        :returns: the result from the optimizer as True or False, depending
4325          on if the refinement converged.
4326        '''
4327        data = self.data['PDF Controls']
4328        xydata = {}
4329        for key in 'Sample Bkg.','Container Bkg.','Container','Sample':
4330            name = data[key]['Name'].strip()
4331            if name:
4332                xydata[key] = self.proj.histogram(name).data['data']
4333                if key != 'Sample': continue
4334                limits = self.proj.histogram(name).data['Limits'][1]
4335                inst = self.proj.histogram(name).data['Instrument Parameters'][0]
4336        res = G2pwd.OptimizePDF(data,xydata,limits,inst,showFit,maxCycles)
4337        return res['success']
4338       
4339    def export(self,fileroot,formats):
4340        '''Write out the PDF-related data (G(r), S(Q),...) into files
4341
4342        :param str fileroot: name of file(s) to be written. The extension
4343          will be ignored and set to .iq, .sq, .fq or .gr depending
4344          on the formats selected.
4345        :param str formats: string specifying the file format(s) to be written,
4346          should contain at least one of the following keywords:
4347          I(Q), S(Q), F(Q), G(r) and/or PDFgui (capitalization and
4348          punctuation is ignored). Note that G(r) and PDFgui should not
4349          be specifed together.
4350        '''
4351        PDFsaves = 5*[False]
4352        PDFentry = self.name
4353        name = self.data['PDF Controls']['Sample']['Name'].strip()
4354        limits = self.proj.histogram(name).data['Limits']
4355        inst = self.proj.histogram(name).data['Instrument Parameters'][0]
4356        for i,lbl in enumerate(['I(Q)', 'S(Q)', 'F(Q)', 'G(r)', 'PDFgui']):
4357            PDFsaves[i] = lbl.lower() in formats.lower()
4358        G2fil.PDFWrite(PDFentry,fileroot,PDFsaves,self.data['PDF Controls'],inst,limits)
4359
4360class G2Image(G2ObjectWrapper):
4361    """Wrapper for an IMG tree entry, containing an image and various metadata.
4362    Note that in a GSASIIscriptable script, instances of G2Image will be created by
4363    calls to :meth:`G2Project.add_image` or :meth:`G2Project.images`, not via calls
4364    to :meth:`G2Image.__init__`.
4365
4366    Example use of G2Image:
4367
4368    >>> gpx = G2sc.G2Project(filename='itest.gpx')
4369    >>> imlst = gpx.add_image(idata,fmthint="TIF")
4370    >>> imlst[0].loadControls('stdSettings.imctrl')
4371    >>> imlst[0].setCalibrant('Si    SRM640c')
4372    >>> imlst[0].loadMasks('stdMasks.immask')
4373    >>> imlst[0].Recalibrate()
4374    >>> imlst[0].setControl('outAzimuths',3)
4375    >>> pwdrList = imlst[0].Integrate()
4376
4377
4378    Sample script using an image:   
4379
4380.. code-block::  python
4381
4382    import os,sys
4383    sys.path.insert(0,'/Users/toby/software/G2/GSASII')
4384    import GSASIIscriptable as G2sc
4385    datadir = "/tmp/V4343_IntegrationError/"
4386    PathWrap = lambda fil: os.path.join(datadir,fil)
4387
4388    gpx = G2sc.G2Project(filename=PathWrap('inttest.gpx'))
4389    imlst = gpx.add_image(PathWrap('Si_free_dc800_1-00000.tif'),fmthint="TIF")
4390    imlst[0].loadControls(PathWrap('Si_free_dc800_1-00000.imctrl'))
4391    pwdrList = imlst[0].Integrate()
4392    gpx.save()
4393
4394    """
4395    # parameters in that can be accessed via setControl. This may need future attention
4396    ControlList = {
4397        'int': ['calibskip', 'pixLimit', 'edgemin', 'outChannels',
4398                    'outAzimuths'],
4399        'float': ['cutoff', 'setdist', 'wavelength', 'Flat Bkg',
4400                      'azmthOff', 'tilt', 'calibdmin', 'rotation',
4401                      'distance', 'DetDepth'],
4402        'bool': ['setRings', 'setDefault', 'centerAzm', 'fullIntegrate',
4403                     'DetDepthRef', 'showLines'],
4404        'str': ['SampleShape', 'binType', 'formatName', 'color',
4405                    'type', ],
4406        'list': ['GonioAngles', 'IOtth', 'LRazimuth', 'Oblique', 'PolaVal',
4407                   'SampleAbs', 'center', 'ellipses', 'linescan',
4408                    'pixelSize', 'range', 'ring', 'rings', 'size', ],
4409        'dict': ['varyList'],
4410        }
4411    '''Defines the items known to exist in the Image Controls tree section
4412    and the item's data types. A few are not included here
4413    ('background image', 'dark image', 'Gain map', and 'calibrant') because
4414    these items have special set routines,
4415    where references to entries are checked to make sure their values are
4416    correct.
4417    ''' 
4418       
4419    def __init__(self, data, name, proj):
4420        self.data = data
4421        self.name = name
4422        self.proj = proj
4423
4424    def setControl(self,arg,value):
4425        '''Set an Image Controls parameter in the current image.
4426        If the parameter is not found an exception is raised.
4427
4428        :param str arg: the name of a parameter (dict entry) in the
4429          image. The parameter must be found in :data:`ControlList`
4430          or an exception is raised.
4431        :param value: the value to set the parameter. The value is
4432          cast as the appropriate type from :data:`ControlList`.
4433        '''
4434        for typ in self.ControlList:
4435            if arg in self.ControlList[typ]: break
4436        else:
4437            print('Allowed args:\n',[nam for nam,typ in self.findControl('')])
4438            raise Exception('arg {} not defined in G2Image.setControl'
4439                                .format(arg))
4440        try:
4441            if typ == 'int':
4442                self.data['Image Controls'][arg] = int(value)
4443            elif typ == 'float':
4444                self.data['Image Controls'][arg] = float(value)
4445            elif typ == 'bool':
4446                self.data['Image Controls'][arg] = bool(value)
4447            elif typ == 'str':
4448                self.data['Image Controls'][arg] = str(value)
4449            elif typ == 'list':
4450                self.data['Image Controls'][arg] = list(value)
4451            elif typ == 'dict':
4452                self.data['Image Controls'][arg] = dict(value)
4453            else:
4454                raise Exception('Unknown type {} for arg {} in  G2Image.setControl'
4455                                    .format(typ,arg))
4456        except:
4457            raise Exception('Error formatting value {} as type {} for arg {} in  G2Image.setControl'
4458                                    .format(value,typ,arg))
4459
4460    def getControl(self,arg):
4461        '''Return an Image Controls parameter in the current image.
4462        If the parameter is not found an exception is raised.
4463
4464        :param str arg: the name of a parameter (dict entry) in the
4465          image.
4466        :returns: the value as a int, float, list,...
4467        '''
4468        if arg in self.data['Image Controls']:
4469            return self.data['Image Controls'][arg]
4470        print(self.findControl(''))
4471        raise Exception('arg {} not defined in G2Image.getControl'.format(arg))
4472
4473    def findControl(self,arg=''):
4474        '''Finds the Image Controls parameter(s) in the current image
4475        that match the string in arg. Default is '' which returns all
4476        parameters.
4477
4478            Example:
4479
4480            >>> findControl('calib')
4481            [['calibskip', 'int'], ['calibdmin', 'float'], ['calibrant', 'str']]
4482
4483        :param str arg: a string containing part of the name of a
4484          parameter (dict entry) in the image's Image Controls.
4485        :returns: a list of matching entries in form
4486          [['item','type'], ['item','type'],...] where each 'item' string
4487          contains the sting in arg.
4488        '''
4489        matchList = []
4490        for typ in self.ControlList:
4491            for item in self.ControlList[typ]:
4492                if arg in item:
4493                    matchList.append([item,typ])
4494        return matchList
4495
4496    def setCalibrant(self,calib):
4497        '''Set a calibrant for the current image
4498
4499        :param str calib: specifies a calibrant name which must be one of
4500          the entries in file ImageCalibrants.py. This is validated and
4501          an error provides a list of valid choices.
4502        '''
4503        import ImageCalibrants as calFile
4504        if calib in calFile.Calibrants.keys():
4505            self.data['Image Controls']['calibrant'] = calib
4506            return
4507        print('Calibrant {} is not valid. Valid calibrants'.format(calib))
4508        for i in calFile.Calibrants.keys():
4509            if i: print('\t"{}"'.format(i))
4510       
4511    def setControlFile(self,typ,imageRef,mult=None):
4512        '''Set a image to be used as a background/dark/gain map image
4513
4514        :param str typ: specifies image type, which must be one of:
4515           'background image', 'dark image', 'gain map'; N.B. only the first
4516           four characters must be specified and case is ignored.
4517        :param imageRef: A reference to the desired image. Either the Image
4518          tree name (str), the image's index (int) or
4519          a image object (:class:`G2Image`)
4520        :param float mult: a multiplier to be applied to the image (not used
4521          for 'Gain map'; required for 'background image', 'dark image'
4522        '''
4523        if 'back' in typ.lower():
4524            key = 'background image'
4525            mult = float(mult)
4526        elif 'dark' in typ.lower():
4527            key = 'dark image'
4528            mult = float(mult)
4529        elif 'gain' in typ.lower():
4530            #key = 'Gain map'
4531            if mult is not None:
4532                print('Ignoring multiplier for Gain map')
4533            mult = None
4534        else:
4535            raise Exception("Invalid typ {} for setControlFile".format(typ))
4536        imgNam = self.proj.image(imageRef).name
4537        if mult is None:
4538            self.data['Image Controls']['Gain map'] = imgNam
4539        else:
4540            self.data['Image Controls'][key] = [imgNam,mult]
4541
4542    def loadControls(self,filename=None,imgDict=None):
4543        '''load controls from a .imctrl file
4544
4545        :param str filename: specifies a file to be read, which should end
4546          with .imctrl (defaults to None, meaning parameters are input
4547          with imgDict.)
4548        :param dict imgDict: contains a set of image parameters (defaults to
4549          None, meaning parameters are input with filename.)
4550        '''
4551        if filename:
4552            File = open(filename,'r')
4553            Slines = File.readlines()
4554            File.close()
4555            G2fil.LoadControls(Slines,self.data['Image Controls'])
4556            print('file {} read into {}'.format(filename,self.name))
4557        elif imgDict:
4558            self.data['Image Controls'].update(imgDict)
4559            print('Image controls set')
4560        else:
4561            raise Exception("loadControls called without imgDict or filename specified")
4562
4563    def saveControls(self,filename):
4564        '''write current controls values to a .imctrl file
4565
4566        :param str filename: specifies a file to write, which should end
4567          with .imctrl
4568        '''
4569        G2fil.WriteControls(filename,self.data['Image Controls'])
4570        print('file {} written from {}'.format(filename,self.name))
4571
4572    def loadMasks(self,filename,ignoreThreshold=False):
4573        '''load masks from a .immask file
4574
4575        :param str filename: specifies a file to be read, which should end
4576          with .immask
4577        :param bool ignoreThreshold: If True, masks are loaded with
4578          threshold masks. Default is False which means any Thresholds
4579          in the file are ignored.
4580        '''
4581        G2fil.readMasks(filename,self.data['Masks'],ignoreThreshold)
4582        print('file {} read into {}'.format(filename,self.name))
4583       
4584    def getVary(self,*args):
4585        '''Return the refinement flag(s) for Image Controls parameter(s)
4586        in the current image.
4587        If the parameter is not found, an exception is raised.
4588
4589        :param str arg: the name of a refinement parameter in the
4590          varyList for the image. The name should be one of
4591          'dep', 'det-X', 'det-Y', 'dist', 'phi', 'tilt', or 'wave'
4592        :param str arg1: the name of a parameter (dict entry) as before,
4593          optional
4594
4595
4596        :returns: a list of bool value(s)
4597        '''
4598        res = []
4599        for arg in args:
4600            if arg in self.data['Image Controls']['varyList']:
4601                res.append(self.data['Image Controls']['varyList'][arg])
4602            else:
4603                raise Exception('arg {} not defined in G2Image.getVary'.format(arg))
4604        return res
4605   
4606    def setVary(self,arg,value):
4607        '''Set a refinement flag for Image Controls parameter in the
4608        current image.
4609        If the parameter is not found an exception is raised.
4610
4611        :param str arg: the name of a refinement parameter in the
4612          varyList for the image. The name should be one of
4613          'dep', 'det-X', 'det-Y', 'dist', 'phi', 'tilt', or 'wave'
4614        :param str arg: the name of a parameter (dict entry) in the
4615          image. The parameter must be found in :data:`ControlList`
4616          or an exception is raised.
4617        :param value: the value to set the parameter. The value is
4618          cast as the appropriate type from :data:`ControlList`.
4619        '''
4620        if arg in self.data['Image Controls']['varyList']:
4621            self.data['Image Controls']['varyList'][arg] = bool(value)
4622        else:
4623            raise Exception('arg {} not defined in G2Image.setVary'.format(arg))
4624
4625    def Recalibrate(self):
4626        '''Invokes a recalibration fit (same as Image Controls/Calibration/Recalibrate
4627        menu command). Note that for this to work properly, the calibration
4628        coefficients (center, wavelength, ditance & tilts) must be fairly close.
4629        This may produce a better result if run more than once.
4630        '''
4631        LoadG2fil()
4632        ImageZ = _getCorrImage(Readers['Image'],self.proj,self)
4633        G2img.ImageRecalibrate(None,ImageZ,self.data['Image Controls'],self.data['Masks'])
4634
4635    def Integrate(self,name=None):
4636        '''Invokes an image integration (same as Image Controls/Integration/Integrate
4637        menu command). All parameters will have previously been set with Image Controls
4638        so no input is needed here. Note that if integration is performed on an
4639        image more than once, histogram entries may be overwritten. Use the name
4640        parameter to prevent this if desired.
4641
4642        :param str name: base name for created histogram(s). If None (default),
4643          the histogram name is taken from the image name.
4644        :returns: a list of created histogram (:class:`G2PwdrData`) objects.
4645        '''
4646        blkSize = 256   #256 seems to be optimal; will break in polymask if >1024
4647        ImageZ = _getCorrImage(Readers['Image'],self.proj,self)
4648        # do integration
4649        ints,azms,Xvals,cancel = G2img.ImageIntegrate(ImageZ,self.data['Image Controls'],self.data['Masks'],blkSize=blkSize)
4650        # code from here on based on G2IO.SaveIntegration, but places results in the current
4651        # project rather than tree
4652        X = Xvals[:-1]
4653        N = len(X)
4654
4655        data = self.data['Image Controls']
4656        Comments = self.data['Comments']
4657        # make name from image, unless overridden
4658        if name:
4659            if not name.startswith(data['type']+' '):
4660                name = data['type']+' '+name
4661        else:
4662            name = self.name.replace('IMG ',data['type']+' ')
4663        if 'PWDR' in name:
4664            if 'target' in data:
4665                names = ['Type','Lam1','Lam2','I(L2)/I(L1)','Zero','Polariz.','U','V','W','X','Y','Z','SH/L','Azimuth'] 
4666                codes = [0 for i in range(14)]
4667            else:
4668                names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','Z','SH/L','Azimuth'] 
4669                codes = [0 for i in range(12)]
4670        elif 'SASD' in name:
4671            names = ['Type','Lam','Zero','Azimuth'] 
4672            codes = [0 for i in range(4)]
4673            X = 4.*np.pi*npsind(X/2.)/data['wavelength']    #convert to q
4674        Xminmax = [X[0],X[-1]]
4675        Azms = []
4676        dazm = 0.
4677        if data['fullIntegrate'] and data['outAzimuths'] == 1:
4678            Azms = [45.0,]                              #a poor man's average?
4679        else:
4680            for i,azm in enumerate(azms[:-1]):
4681                if azm > 360. and azms[i+1] > 360.:
4682                    Azms.append(G2img.meanAzm(azm%360.,azms[i+1]%360.))
4683                else:   
4684                    Azms.append(G2img.meanAzm(azm,azms[i+1]))
4685            dazm = np.min(np.abs(np.diff(azms)))/2.
4686        # pull out integration results and make histograms for each
4687        IntgOutList = []
4688        for i,azm in enumerate(azms[:-1]):
4689            Aname = name+" Azm= %.2f"%((azm+dazm)%360.)
4690            # MT dict to contain histogram
4691            HistDict = {}
4692            histItems = [Aname]
4693            Sample = G2obj.SetDefaultSample()       #set as Debye-Scherrer
4694            Sample['Gonio. radius'] = data['distance']
4695            Sample['Omega'] = data['GonioAngles'][0]
4696            Sample['Chi'] = data['GonioAngles'][1]
4697            Sample['Phi'] = data['GonioAngles'][2]
4698            Sample['Azimuth'] = (azm+dazm)%360.    #put here as bin center
4699            polariz = 0.99    #set default polarization for synchrotron radiation!
4700            for item in Comments:
4701                if 'polariz' in item:
4702                    try:
4703                        polariz = float(item.split('=')[1])
4704                    except:
4705                        polariz = 0.99
4706                for key in ('Temperature','Pressure','Time','FreePrm1','FreePrm2','FreePrm3','Omega',
4707                    'Chi','Phi'):
4708                    if key.lower() in item.lower():
4709                        try:
4710                            Sample[key] = float(item.split('=')[1])
4711                        except:
4712                            pass
4713                # if 'label_prm' in item.lower():
4714                #     for num in ('1','2','3'):
4715                #         if 'label_prm'+num in item.lower():
4716                #             Controls['FreePrm'+num] = item.split('=')[1].strip()
4717            if 'PWDR' in Aname:
4718                if 'target' in data:    #from lab x-ray 2D imaging data
4719                    waves = {'CuKa':[1.54051,1.54433],'TiKa':[2.74841,2.75207],'CrKa':[2.28962,2.29351],
4720                                 'FeKa':[1.93597,1.93991],'CoKa':[1.78892,1.79278],'MoKa':[0.70926,0.713543],
4721                                 'AgKa':[0.559363,0.563775]}
4722                    wave1,wave2 = waves[data['target']]
4723                    parms = ['PXC',wave1,wave2,0.5,0.0,polariz,290.,-40.,30.,6.,-14.,0.0,0.0001,Azms[i]]
4724                else:
4725                    parms = ['PXC',data['wavelength'],0.0,polariz,1.0,-0.10,0.4,0.30,1.0,0.0,0.0001,Azms[i]]
4726            elif 'SASD' in Aname:
4727                Sample['Trans'] = data['SampleAbs'][0]
4728                parms = ['LXC',data['wavelength'],0.0,Azms[i]]
4729            Y = ints[i]
4730            Ymin = np.min(Y)
4731            Ymax = np.max(Y)
4732            W = np.where(Y>0.,1./Y,1.e-6)                    #probably not true
4733            section = 'Comments'
4734            histItems += [section]
4735            HistDict[section] = Comments
4736            section = 'Limits'
4737            histItems += [section]
4738            HistDict[section] = copy.deepcopy([tuple(Xminmax),Xminmax])
4739            if 'PWDR' in Aname:
4740                section = 'Background'
4741                histItems += [section]
4742                HistDict[section] = [['chebyschev',1,3,1.0,0.0,0.0],
4743                    {'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}]
4744            inst = [dict(zip(names,zip(parms,parms,codes))),{}]
4745            for item in inst[0]:
4746                inst[0][item] = list(inst[0][item])
4747            section = 'Instrument Parameters'
4748            histItems += [section]
4749            HistDict[section] = inst
4750            if 'PWDR' in Aname:
4751                section = 'Sample Parameters'
4752                histItems += [section]
4753                HistDict[section] = Sample
4754                section = 'Peak List'
4755                histItems += [section]
4756                HistDict[section] = {'sigDict':{},'peaks':[]}
4757                section = 'Index Peak List'
4758                histItems += [section]
4759                HistDict[section] = [[],[]]
4760                section = 'Unit Cells List'
4761                histItems += [section]
4762                HistDict[section] = []
4763                section = 'Reflection Lists'
4764                histItems += [section]
4765                HistDict[section] = {}
4766            # elif 'SASD' in Aname:             
4767            #     section = 'Substances'
4768            #     histItems += [section]
4769            #     HistDict[section] = G2pdG.SetDefaultSubstances()  # this needs to be moved
4770            #     section = 'Sample Parameters'
4771            #     histItems += [section]
4772            #     HistDict[section] = Sample
4773            #     section = 'Models'
4774            #     histItems += [section]
4775            #     HistDict[section] = G2pdG.SetDefaultSASDModel() # this needs to be moved
4776            valuesdict = {
4777                'wtFactor':1.0,'Dummy':False,'ranId':ran.randint(0,sys.maxsize),'Offset':[0.0,0.0],'delOffset':0.02*Ymax,
4778                'refOffset':-0.1*Ymax,'refDelt':0.1*Ymax,'Yminmax':[Ymin,Ymax]}
4779            # if Aname is already in the project replace it
4780            for j in self.proj.names:
4781                if j[0] == Aname: 
4782                    print('Replacing "{}" in project'.format(Aname))
4783                    break
4784            else:
4785                print('Adding "{}" to project'.format(Aname))
4786                self.proj.names.append([Aname]+
4787                        [u'Comments',u'Limits',u'Background',u'Instrument Parameters',
4788                         u'Sample Parameters', u'Peak List', u'Index Peak List',
4789                         u'Unit Cells List', u'Reflection Lists'])
4790            HistDict['data'] = [valuesdict,
4791                    [np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]]
4792            self.proj.data[Aname] = HistDict
4793            IntgOutList.append(self.proj.histogram(Aname))
4794        return IntgOutList
4795
4796##########################
4797# Command Line Interface #
4798##########################
4799# Each of these takes an argparse.Namespace object as their argument,
4800# representing the parsed command-line arguments for the relevant subcommand.
4801# The argument specification for each is in the subcommands dictionary (see
4802# below)
4803
4804commandhelp={}
4805commandhelp["create"] = "creates a GSAS-II project, optionally adding histograms and/or phases"
4806def create(args):
4807    """Implements the create command-line subcommand. This creates a GSAS-II project, optionally adding histograms and/or phases::
4808
4809  usage: GSASIIscriptable.py create [-h] [-d HISTOGRAMS [HISTOGRAMS ...]]
4810                                  [-i IPARAMS [IPARAMS ...]]
4811                                  [-p PHASES [PHASES ...]]
4812                                  filename
4813                                 
4814positional arguments::
4815
4816  filename              the project file to create. should end in .gpx
4817
4818optional arguments::
4819
4820  -h, --help            show this help message and exit
4821  -d HISTOGRAMS [HISTOGRAMS ...], --histograms HISTOGRAMS [HISTOGRAMS ...]
4822                        list of datafiles to add as histograms
4823  -i IPARAMS [IPARAMS ...], --iparams IPARAMS [IPARAMS ...]
4824                        instrument parameter file, must be one for every
4825                        histogram
4826  -p PHASES [PHASES ...], --phases PHASES [PHASES ...]
4827                        list of phases to add. phases are automatically
4828                        associated with all histograms given.
4829
4830    """
4831    proj = G2Project(gpxfile=args.filename)
4832
4833    hist_objs = []
4834    if args.histograms:
4835        for h,i in zip(args.histograms,args.iparams):
4836            print("Adding histogram from",h,"with instparm ",i)
4837            hist_objs.append(proj.add_powder_histogram(h, i))
4838
4839    if args.phases: 
4840        for p in args.phases:
4841            print("Adding phase from",p)
4842            proj.add_phase(p, histograms=hist_objs)
4843        print('Linking phase(s) to histogram(s):')
4844        for h in hist_objs:
4845            print ('   '+h.name)
4846
4847    proj.save()
4848
4849commandhelp["add"] = "adds histograms and/or phases to GSAS-II project"
4850def add(args):
4851    """Implements the add command-line subcommand. This adds histograms and/or phases to GSAS-II project::
4852
4853  usage: GSASIIscriptable.py add [-h] [-d HISTOGRAMS [HISTOGRAMS ...]]
4854                               [-i IPARAMS [IPARAMS ...]]
4855                               [-hf HISTOGRAMFORMAT] [-p PHASES [PHASES ...]]
4856                               [-pf PHASEFORMAT] [-l HISTLIST [HISTLIST ...]]
4857                               filename
4858
4859
4860positional arguments::
4861
4862  filename              the project file to open. Should end in .gpx
4863
4864optional arguments::
4865
4866  -h, --help            show this help message and exit
4867  -d HISTOGRAMS [HISTOGRAMS ...], --histograms HISTOGRAMS [HISTOGRAMS ...]
4868                        list of datafiles to add as histograms
4869  -i IPARAMS [IPARAMS ...], --iparams IPARAMS [IPARAMS ...]
4870                        instrument parameter file, must be one for every
4871                        histogram
4872  -hf HISTOGRAMFORMAT, --histogramformat HISTOGRAMFORMAT
4873                        format hint for histogram import. Applies to all
4874                        histograms
4875  -p PHASES [PHASES ...], --phases PHASES [PHASES ...]
4876                        list of phases to add. phases are automatically
4877                        associated with all histograms given.
4878  -pf PHASEFORMAT, --phaseformat PHASEFORMAT
4879                        format hint for phase import. Applies to all phases.
4880                        Example: -pf CIF
4881  -l HISTLIST [HISTLIST ...], --histlist HISTLIST [HISTLIST ...]
4882                        list of histgram indices to associate with added
4883                        phases. If not specified, phases are associated with
4884                        all previously loaded histograms. Example: -l 2 3 4
4885   
4886    """
4887    proj = G2Project(args.filename)
4888
4889    if args.histograms:
4890        for h,i in zip(args.histograms,args.iparams):
4891            print("Adding histogram from",h,"with instparm ",i)
4892            proj.add_powder_histogram(h, i, fmthint=args.histogramformat)
4893
4894    if args.phases: 
4895        if not args.histlist:
4896            histlist = proj.histograms()
4897        else:
4898            histlist = [proj.histogram(i) for i in args.histlist]
4899
4900        for p in args.phases:
4901            print("Adding phase from",p)
4902            proj.add_phase(p, histograms=histlist, fmthint=args.phaseformat)
4903           
4904        if not args.histlist:
4905            print('Linking phase(s) to all histogram(s)')
4906        else:
4907            print('Linking phase(s) to histogram(s):')
4908            for h in histlist:
4909                print ('   '+h.name)
4910
4911    proj.save()
4912
4913
4914commandhelp["dump"] = "Shows the contents of a GSAS-II project"
4915def dump(args):
4916    """Implements the dump command-line subcommand, which shows the contents of a GSAS-II project::
4917
4918       usage: GSASIIscriptable.py dump [-h] [-d] [-p] [-r] files [files ...]
4919
4920positional arguments::
4921
4922  files
4923
4924optional arguments::
4925
4926  -h, --help        show this help message and exit
4927  -d, --histograms  list histograms in files, overrides --raw
4928  -p, --phases      list phases in files, overrides --raw
4929  -r, --raw         dump raw file contents, default
4930 
4931    """
4932    if not args.histograms and not args.phases:
4933        args.raw = True
4934    if args.raw:
4935        import IPython.lib.pretty as pretty
4936
4937    for fname in args.files:
4938        if args.raw:
4939            proj, nameList = LoadDictFromProjFile(fname)
4940            print("file:", fname)
4941            print("names:", nameList)
4942            for key, val in proj.items():
4943                print(key, ":")
4944                pretty.pprint(val)
4945        else:
4946            proj = G2Project(fname)
4947            if args.histograms:
4948                hists = proj.histograms()
4949                for h in hists:
4950                    print(fname, "hist", h.id, h.name)
4951            if args.phases:
4952                phase_list = proj.phases()
4953                for p in phase_list:
4954                    print(fname, "phase", p.id, p.name)
4955
4956
4957commandhelp["browse"] = "Load a GSAS-II project and then open a IPython shell to browse it"
4958def IPyBrowse(args):
4959    """Load a .gpx file and then open a IPython shell to browse it::
4960
4961  usage: GSASIIscriptable.py browse [-h] files [files ...]
4962
4963positional arguments::
4964
4965  files       list of files to browse
4966
4967optional arguments::
4968
4969  -h, --help  show this help message and exit
4970
4971    """
4972    for fname in args.files:
4973        proj, nameList = LoadDictFromProjFile(fname)
4974        msg = "\nfname {} loaded into proj (dict) with names in nameList".format(fname)
4975        GSASIIpath.IPyBreak_base(msg)
4976        break
4977
4978
4979commandhelp["refine"] = '''
4980Conducts refinements on GSAS-II projects according to a list of refinement
4981steps in a JSON dict
4982'''
4983def refine(args):
4984    """Implements the refine command-line subcommand:
4985    conducts refinements on GSAS-II projects according to a JSON refinement dict::
4986
4987        usage: GSASIIscriptable.py refine [-h] gpxfile [refinements]
4988
4989positional arguments::
4990
4991  gpxfile      the project file to refine
4992  refinements  json file of refinements to apply. if not present refines file
4993               as-is
4994
4995optional arguments::
4996
4997  -h, --help   show this help message and exit
4998 
4999    """
5000    proj = G2Project(args.gpxfile)
5001    if args.refinements is None:
5002        proj.refine()
5003    else:
5004        import json
5005        with open(args.refinements) as refs:
5006            refs = json.load(refs)
5007        if type(refs) is not dict:
5008            raise G2ScriptException("Error: JSON object must be a dict.")
5009        if "code" in refs:
5010            print("executing code:\n|  ",'\n|  '.join(refs['code']))
5011            exec('\n'.join(refs['code']))
5012        proj.do_refinements(refs['refinements'])
5013
5014
5015commandhelp["export"] = "Export phase as CIF"
5016def export(args):
5017    """Implements the export command-line subcommand: Exports phase as CIF::
5018
5019      usage: GSASIIscriptable.py export [-h] gpxfile phase exportfile
5020
5021positional arguments::
5022
5023  gpxfile     the project file from which to export
5024  phase       identifier of phase to export
5025  exportfile  the .cif file to export to
5026
5027optional arguments::
5028
5029  -h, --help  show this help message and exit
5030
5031    """
5032    proj = G2Project(args.gpxfile)
5033    phase = proj.phase(args.phase)
5034    phase.export_CIF(args.exportfile)
5035
5036
5037def _args_kwargs(*args, **kwargs):
5038    return args, kwargs
5039
5040# A dictionary of the name of each subcommand, and a tuple
5041# of its associated function and a list of its arguments
5042# The arguments are passed directly to the add_argument() method
5043# of an argparse.ArgumentParser
5044
5045subcommands = {"create":
5046               (create, [_args_kwargs('filename',
5047                                      help='the project file to create. should end in .gpx'),
5048
5049                         _args_kwargs('-d', '--histograms',
5050                                      nargs='+',
5051                                      help='list of datafiles to add as histograms'),
5052                                     
5053                         _args_kwargs('-i', '--iparams',
5054                                      nargs='+',
5055                                      help='instrument parameter file, must be one'
5056                                           ' for every histogram'
5057                                      ),
5058
5059                         _args_kwargs('-p', '--phases',
5060                                      nargs='+',
5061                                      help='list of phases to add. phases are '
5062                                           'automatically associated with all '
5063                                           'histograms given.')]),
5064               "add": (add, [_args_kwargs('filename',
5065                                      help='the project file to open. Should end in .gpx'),
5066
5067                         _args_kwargs('-d', '--histograms',
5068                                      nargs='+',
5069                                      help='list of datafiles to add as histograms'),
5070                                     
5071                         _args_kwargs('-i', '--iparams',
5072                                      nargs='+',
5073                                      help='instrument parameter file, must be one'
5074                                           ' for every histogram'
5075                                      ),
5076                                     
5077                         _args_kwargs('-hf', '--histogramformat',
5078                                      help='format hint for histogram import. Applies to all'
5079                                           ' histograms'
5080                                      ),
5081
5082                         _args_kwargs('-p', '--phases',
5083                                      nargs='+',
5084                                      help='list of phases to add. phases are '
5085                                           'automatically associated with all '
5086                                           'histograms given.'),
5087
5088                         _args_kwargs('-pf', '--phaseformat',
5089                                      help='format hint for phase import. Applies to all'
5090                                           ' phases. Example: -pf CIF'
5091                                      ),
5092                                     
5093                         _args_kwargs('-l', '--histlist',
5094                                      nargs='+',
5095                                      help='list of histgram indices to associate with added'
5096                                           ' phases. If not specified, phases are'
5097                                           ' associated with all previously loaded'
5098                                           ' histograms. Example: -l 2 3 4')]),
5099                                           
5100               "dump": (dump, [_args_kwargs('-d', '--histograms',
5101                                     action='store_true',
5102                                     help='list histograms in files, overrides --raw'),
5103
5104                               _args_kwargs('-p', '--phases',
5105                                            action='store_true',
5106                                            help='list phases in files, overrides --raw'),
5107
5108                               _args_kwargs('-r', '--raw',
5109                                      action='store_true', help='dump raw file contents, default'),
5110
5111                               _args_kwargs('files', nargs='+')]),
5112
5113               "refine":
5114               (refine, [_args_kwargs('gpxfile', help='the project file to refine'),
5115                         _args_kwargs('refinements',
5116                                      help='JSON file of refinements to apply. if not present'
5117                                           ' refines file as-is',
5118                                      default=None,
5119                                      nargs='?')]),
5120
5121               "export": (export, [_args_kwargs('gpxfile',
5122                                                help='the project file from which to export'),
5123                                   _args_kwargs('phase', help='identifier of phase to export'),
5124                                   _args_kwargs('exportfile', help='the .cif file to export to')]),
5125               "browse": (IPyBrowse, [_args_kwargs('files', nargs='+',
5126                                                   help='list of files to browse')])}
5127
5128
5129def main():
5130    '''The command-line interface for calling GSASIIscriptable as a shell command,
5131    where it is expected to be called as::
5132
5133       python GSASIIscriptable.py <subcommand> <file.gpx> <options>
5134
5135    The following subcommands are defined:
5136
5137        * create, see :func:`create`
5138        * add, see :func:`add`
5139        * dump, see :func:`dump`
5140        * refine, see :func:`refine`
5141        * export, :func:`export`
5142        * browse, see :func:`IPyBrowse`
5143
5144    .. seealso::
5145        :func:`create`
5146        :func:`add`
5147        :func:`dump`
5148        :func:`refine`
5149        :func:`export`
5150        :func:`IPyBrowse`
5151    '''
5152    parser = argparse.ArgumentParser(description=
5153        "Use of "+os.path.split(__file__)[1]+" Allows GSAS-II actions from command line."
5154        )
5155    subs = parser.add_subparsers()
5156
5157    # Create all of the specified subparsers
5158    for name, (func, args) in subcommands.items():
5159        new_parser = subs.add_parser(name,help=commandhelp.get(name),
5160                                     description='Command "'+name+'" '+commandhelp.get(name))
5161        for listargs, kwds in args:
5162            new_parser.add_argument(*listargs, **kwds)
5163        new_parser.set_defaults(func=func)
5164
5165    # Parse and trigger subcommand
5166    result = parser.parse_args()
5167    result.func(result)
5168
5169if __name__ == '__main__':
5170    #fname='/tmp/corundum-template.gpx'
5171    #prj = G2Project(fname)
5172    main()
Note: See TracBrowser for help on using the repository browser.