source: trunk/GSASIIscriptable.py @ 4119

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

add reflection output in scripting

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