source: trunk/GSASIIscriptable.py @ 4053

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

complete multiprocessing & caching of maps; results not yet tested

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