source: trunk/GSASIIscriptable.py @ 3834

Last change on this file since 3834 was 3834, checked in by toby, 3 years ago

set idstring for simulations as needed; thanx to Thomas Proffen

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