source: trunk/GSASIIscriptable.py @ 3840

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

fix scriptable read of TOF; fraction editor for wx 3.0

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 177.9 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3########### SVN repository information ###################
4# $Date: 2019-03-07 00:18:01 +0000 (Thu, 07 Mar 2019) $
5# $Author: toby $
6# $Revision: 3840 $
7# $URL: trunk/GSASIIscriptable.py $
8# $Id: GSASIIscriptable.py 3840 2019-03-07 00:18:01Z 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    if 'T' in Iparm1['Type'][0]:
1196        if not reader.clockWd and reader.GSAS:
1197            reader.powderdata[0] *= 100.        #put back the CW centideg correction
1198
1199    # Ending keys:
1200    # [u'Reflection Lists',
1201    #  u'Limits',
1202    #  'data',
1203    #  u'Index Peak List',
1204    #  u'Comments',
1205    #  u'Unit Cells List',
1206    #  u'Sample Parameters',
1207    #  u'Peak List',
1208    #  u'Background',
1209    #  u'Instrument Parameters']
1210    Tmin = np.min(reader.powderdata[0])
1211    Tmax = np.max(reader.powderdata[0])
1212
1213    default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0],
1214                          {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}]
1215
1216    output_dict = {u'Reflection Lists': {},
1217                   u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]),
1218                   u'data': [valuesdict, reader.powderdata, HistName],
1219                   u'Index Peak List': [[], []],
1220                   u'Comments': reader.comments,
1221                   u'Unit Cells List': [],
1222                   u'Sample Parameters': reader.Sample,
1223                   u'Peak List': {'peaks': [], 'sigDict': {}},
1224                   u'Background': reader.pwdparms.get('Background', default_background),
1225                   u'Instrument Parameters': [Iparm1, Iparm2],
1226                   }
1227
1228    names = [u'Comments',
1229             u'Limits',
1230             u'Background',
1231             u'Instrument Parameters',
1232             u'Sample Parameters',
1233             u'Peak List',
1234             u'Index Peak List',
1235             u'Unit Cells List',
1236             u'Reflection Lists']
1237
1238    # TODO controls?? GSASII.py:1664-7
1239
1240    return HistName, [HistName] + names, output_dict
1241
1242
1243def _deep_copy_into(from_, into):
1244    """Helper function for reloading .gpx file. See G2Project.reload()
1245
1246    :author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1247    """
1248    if isinstance(from_, dict) and isinstance(into, dict):
1249        combined_keys = set(from_.keys()).union(into.keys())
1250        for key in combined_keys:
1251            if key in from_ and key in into:
1252                both_dicts = (isinstance(from_[key], dict)
1253                              and isinstance(into[key], dict))
1254                both_lists = (isinstance(from_[key], list)
1255                              and isinstance(into[key], list))
1256                if both_dicts or both_lists:
1257                    _deep_copy_into(from_[key], into[key])
1258                else:
1259                    into[key] = from_[key]
1260            elif key in from_:
1261                into[key] = from_[key]
1262            else:  # key in into
1263                del into[key]
1264    elif isinstance(from_, list) and isinstance(into, list):
1265        if len(from_) == len(into):
1266            for i in range(len(from_)):
1267                both_dicts = (isinstance(from_[i], dict)
1268                              and isinstance(into[i], dict))
1269                both_lists = (isinstance(from_[i], list)
1270                              and isinstance(into[i], list))
1271                if both_dicts or both_lists:
1272                    _deep_copy_into(from_[i], into[i])
1273                else:
1274                    into[i] = from_[i]
1275        else:
1276            into[:] = from_
1277
1278def GetCorrImage(ImageReaderlist,proj,imageRef):
1279    '''Gets image & applies dark, background & flat background corrections.
1280    based on :func:`GSASIIimgGUI.GetImageZ`. Likely for internal use only.
1281
1282    :param list ImageReaderlist: list of Reader objects for images
1283    :param object ImageReaderlist: list of Reader objects for images
1284    :param imageRef: A reference to the desired image. Either the Image
1285      tree name (str), the image's index (int) or
1286      a image object (:class:`G2Image`)
1287
1288    :return: array sumImg: corrected image for background/dark/flat back
1289    '''
1290    ImgObj = proj.image(imageRef)
1291    Controls = ImgObj.data['Image Controls']
1292    formatName = Controls.get('formatName','')
1293    imagefile = ImgObj.data['data'][1]
1294    ImageTag = None # fix this for multiimage files
1295    sumImg = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
1296    if sumImg is None:
1297        return []
1298    darkImg = False
1299    if 'dark image' in Controls:
1300        darkImg,darkScale = Controls['dark image']
1301        if darkImg:
1302            dImgObj = proj.image(darkImg)
1303            formatName = dImgObj.data['Image Controls'].get('formatName','')
1304            imagefile = dImgObj.data['data'][1]
1305            ImageTag = None # fix this for multiimage files
1306            darkImg = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
1307            if darkImg is None:
1308                raise Exception('Error reading dark image {}'.format(imagefile))
1309            sumImg += np.array(darkImage*darkScale,dtype='int32')
1310    if 'background image' in Controls:
1311        backImg,backScale = Controls['background image']           
1312        if backImg:     #ignores any transmission effect in the background image
1313            bImgObj = proj.image(backImg)
1314            formatName = bImgObj.data['Image Controls'].get('formatName','')
1315            imagefile = bImgObj.data['data'][1]
1316            ImageTag = None # fix this for multiimage files
1317            backImg = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
1318            if backImage is None:
1319                raise Exception('Error reading background image {}'.format(imagefile))
1320            if darkImg:
1321                backImage += np.array(darkImage*darkScale/backScale,dtype='int32')
1322            else:
1323                sumImg += np.array(backImage*backScale,dtype='int32')
1324    if 'Gain map' in Controls:
1325        gainMap = Controls['Gain map']
1326        if gainMap:
1327            gImgObj = proj.image(gainMap)
1328            formatName = gImgObj.data['Image Controls'].get('formatName','')
1329            imagefile = gImgObj.data['data'][1]
1330            ImageTag = None # fix this for multiimage files
1331            GMimage = G2fil.RereadImageData(ImageReaderlist,imagefile,ImageTag=ImageTag,FormatName=formatName)
1332            if GMimage is None:
1333                raise Exception('Error reading Gain map image {}'.format(imagefile))
1334            sumImg = sumImg*GMimage/1000
1335    sumImg -= int(Controls.get('Flat Bkg',0))
1336    Imax = np.max(sumImg)
1337    Controls['range'] = [(0,Imax),[0,Imax]]
1338    return np.asarray(sumImg,dtype='int32')
1339
1340class G2ObjectWrapper(object):
1341    """Base class for all GSAS-II object wrappers.
1342
1343    The underlying GSAS-II format can be accessed as `wrapper.data`. A number
1344    of overrides are implemented so that the wrapper behaves like a dictionary.
1345
1346    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1347    """
1348    def __init__(self, datadict):
1349        self.data = datadict
1350
1351    def __getitem__(self, key):
1352        return self.data[key]
1353
1354    def __setitem__(self, key, value):
1355        self.data[key] = value
1356
1357    def __contains__(self, key):
1358        return key in self.data
1359
1360    def get(self, k, d=None):
1361        return self.data.get(k, d)
1362
1363    def keys(self):
1364        return self.data.keys()
1365
1366    def values(self):
1367        return self.data.values()
1368
1369    def items(self):
1370        return self.data.items()
1371
1372
1373class G2Project(G2ObjectWrapper):   
1374    """Represents an entire GSAS-II project.
1375
1376    :param str gpxfile: Existing .gpx file to be loaded. If nonexistent,
1377            creates an empty project.
1378    :param str author: Author's name (not yet implemented)
1379    :param str newgpx: The filename the project should be saved to in
1380            the future. If both newgpx and gpxfile are present, the project is
1381            loaded from the gpxfile, then when saved will be written to newgpx.
1382    :param str filename: Name to be used to save the project. Has same function as
1383            parameter newgpx (do not use both gpxfile and filename). Use of newgpx
1384            is preferred over filename.
1385
1386    There are two ways to initialize this object:
1387
1388    >>> # Load an existing project file
1389    >>> proj = G2Project('filename.gpx')
1390   
1391    >>> # Create a new project
1392    >>> proj = G2Project(newgpx='new_file.gpx')
1393   
1394    Histograms can be accessed easily.
1395
1396    >>> # By name
1397    >>> hist = proj.histogram('PWDR my-histogram-name')
1398   
1399    >>> # Or by index
1400    >>> hist = proj.histogram(0)
1401    >>> assert hist.id == 0
1402   
1403    >>> # Or by random id
1404    >>> assert hist == proj.histogram(hist.ranId)
1405
1406    Phases can be accessed the same way.
1407
1408    >>> phase = proj.phase('name of phase')
1409
1410    New data can also be loaded via :meth:`~G2Project.add_phase` and
1411    :meth:`~G2Project.add_powder_histogram`.
1412
1413    >>> hist = proj.add_powder_histogram('some_data_file.chi',
1414                                         'instrument_parameters.prm')
1415    >>> phase = proj.add_phase('my_phase.cif', histograms=[hist])
1416
1417    Parameters for Rietveld refinement can be turned on and off as well.
1418    See :meth:`~G2Project.set_refinement`, :meth:`~G2Project.clear_refinements`,
1419    :meth:`~G2Project.iter_refinements`, :meth:`~G2Project.do_refinements`.
1420    """
1421    def __init__(self, gpxfile=None, author=None, filename=None, newgpx=None):
1422        if filename is not None and newgpx is not None:
1423            raise G2ScriptException('Do not use filename and newgpx together')
1424        elif newgpx is not None:
1425            filename = newgpx
1426        if gpxfile is None:
1427            filename = os.path.abspath(os.path.expanduser(filename))
1428            self.filename = filename
1429            self.data, self.names = make_empty_project(author=author, filename=filename)
1430        elif isinstance(gpxfile, str): # TODO: create replacement for isinstance that checks if path exists
1431                                       # filename is valid, etc.
1432            if isinstance(filename, str): 
1433                self.filename = os.path.abspath(os.path.expanduser(filename)) # both are defined
1434            else: 
1435                self.filename = os.path.abspath(os.path.expanduser(gpxfile))
1436            # TODO set author
1437            self.data, self.names = LoadDictFromProjFile(gpxfile)
1438            self.update_ids()
1439        else:
1440            raise ValueError("Not sure what to do with gpxfile")
1441
1442    @classmethod
1443    def from_dict_and_names(cls, gpxdict, names, filename=None):
1444        """Creates a :class:`G2Project` directly from
1445        a dictionary and a list of names. If in doubt, do not use this.
1446
1447        :returns: a :class:`G2Project`
1448        """
1449        out = cls()
1450        if filename:
1451            filename = os.path.abspath(os.path.expanduser(filename))
1452            out.filename = filename
1453            gpxdict['Controls']['data']['LastSavedAs'] = filename
1454        else:
1455            try:
1456                out.filename = gpxdict['Controls']['data']['LastSavedAs']
1457            except KeyError:
1458                out.filename = None
1459        out.data = gpxdict
1460        out.names = names
1461
1462    def save(self, filename=None):
1463        """Saves the project, either to the current filename, or to a new file.
1464
1465        Updates self.filename if a new filename provided"""
1466        # TODO update LastSavedUsing ?
1467        if filename:
1468            filename = os.path.abspath(os.path.expanduser(filename))
1469            self.data['Controls']['data']['LastSavedAs'] = filename
1470            self.filename = filename
1471        elif not self.filename:
1472            raise AttributeError("No file name to save to")
1473        SaveDictToProjFile(self.data, self.names, self.filename)
1474
1475    def add_powder_histogram(self, datafile, iparams, phases=[], fmthint=None,
1476                                 databank=None, instbank=None):
1477        """Loads a powder data histogram into the project.
1478
1479        Automatically checks for an instrument parameter file, or one can be
1480        provided. Note that in unix fashion, "~" can be used to indicate the
1481        home directory (e.g. ~/G2data/data.fxye).
1482
1483        :param str datafile: The powder data file to read, a filename.
1484        :param str iparams: The instrument parameters file, a filename.
1485        :param list phases: Phases to link to the new histogram
1486        :param str fmthint: If specified, only importers where the format name
1487          (reader.formatName, as shown in Import menu) contains the
1488          supplied string will be tried as importers. If not specified, all
1489          importers consistent with the file extension will be tried
1490          (equivalent to "guess format" in menu).
1491        :param int databank: Specifies a dataset number to read, if file contains
1492          more than set of data. This should be 1 to read the first bank in
1493          the file (etc.) regardless of the number on the Bank line, etc.
1494          Default is None which means there should only be one dataset in the
1495          file.
1496        :param int instbank: Specifies an instrument parameter set to read, if
1497          the instrument parameter file contains more than set of parameters.
1498          This will match the INS # in an GSAS type file so it will typically
1499          be 1 to read the first parameter set in the file (etc.)
1500          Default is None which means there should only be one parameter set
1501          in the file.
1502
1503        :returns: A :class:`G2PwdrData` object representing
1504            the histogram
1505        """
1506        LoadG2fil()
1507        datafile = os.path.abspath(os.path.expanduser(datafile))
1508        iparams = os.path.abspath(os.path.expanduser(iparams))
1509        pwdrreaders = import_generic(datafile, Readers['Pwdr'],fmthint=fmthint,bank=databank)
1510        histname, new_names, pwdrdata = load_pwd_from_reader(
1511                                          pwdrreaders[0], iparams,
1512                                          [h.name for h in self.histograms()],bank=instbank)
1513        if histname in self.data:
1514            print("Warning - redefining histogram", histname)
1515        elif self.names[-1][0] == 'Phases':
1516            self.names.insert(-1, new_names)
1517        else:
1518            self.names.append(new_names)
1519        self.data[histname] = pwdrdata
1520        self.update_ids()
1521
1522        for phase in phases:
1523            phase = self.phase(phase)
1524            self.link_histogram_phase(histname, phase)
1525
1526        return self.histogram(histname)
1527
1528    def add_simulated_powder_histogram(self, histname, iparams, Tmin, Tmax, Tstep,
1529                                       wavelength=None, scale=None, phases=[]):
1530        """Loads a powder data histogram into the project.
1531
1532        Requires an instrument parameter file.
1533        Note that in unix fashion, "~" can be used to indicate the
1534        home directory (e.g. ~/G2data/data.prm). The instrument parameter file
1535        will determine if the histogram is x-ray, CW neutron, TOF, etc. as well
1536        as the instrument type.
1537
1538        :param str histname: A name for the histogram to be created.
1539        :param str iparams: The instrument parameters file, a filename.
1540        :param float Tmin: Minimum 2theta or TOF (ms) for dataset to be simulated
1541        :param float Tmax: Maximum 2theta or TOF (ms) for dataset to be simulated
1542        :param float Tstep: Step size in 2theta or TOF (ms) for dataset to be simulated       
1543        :param float wavelength: Wavelength for CW instruments, overriding the value
1544           in the instrument parameters file if specified.
1545        :param float scale: Histogram scale factor which multiplies the pattern. Note that
1546           simulated noise is added to the pattern, so that if the maximum intensity is
1547           small, the noise will mask the computed pattern. The scale
1548           needs to be a large number for CW neutrons.
1549           The default, None, provides a scale of 1 for x-rays and TOF; 10,000 for CW neutrons.
1550        :param list phases: Phases to link to the new histogram. Use proj.phases() to link to
1551           all defined phases.
1552
1553        :returns: A :class:`G2PwdrData` object representing the histogram
1554        """
1555        LoadG2fil()
1556        iparams = os.path.abspath(os.path.expanduser(iparams))
1557        if not os.path.exists(iparams):
1558            raise G2ScriptException("File does not exist:"+iparams)
1559        rd = G2obj.ImportPowderData( # Initialize a base class reader
1560            extensionlist=tuple(),
1561            strictExtension=False,
1562            formatName = 'Simulate dataset',
1563            longFormatName = 'Compute a simulated pattern')
1564        rd.powderentry[0] = '' # no filename
1565        rd.powderentry[2] = 1 # only one bank
1566        rd.comments.append('This is a dummy dataset for powder pattern simulation')
1567        rd.idstring = histname
1568        #Iparm1, Iparm2 = load_iprms(iparams, rd)
1569        if Tmax < Tmin:
1570            Tmin,Tmax = Tmax,Tmin
1571        Tstep = abs(Tstep)
1572        if 'TOF' in rd.idstring:
1573                N = (np.log(Tmax)-np.log(Tmin))/Tstep
1574                x = np.exp((np.arange(0,N))*Tstep+np.log(Tmin*1000.))
1575                N = len(x)
1576        else:           
1577                N = int((Tmax-Tmin)/Tstep)+1
1578                x = np.linspace(Tmin,Tmax,N,True)
1579                N = len(x)
1580        if N < 3:
1581            raise G2ScriptException("Error: Range is too small or step is too large, <3 points")
1582        rd.powderdata = [
1583            np.array(x), # x-axis values
1584            np.zeros_like(x), # powder pattern intensities
1585            np.ones_like(x), # 1/sig(intensity)^2 values (weights)
1586            np.zeros_like(x), # calc. intensities (zero)
1587            np.zeros_like(x), # calc. background (zero)
1588            np.zeros_like(x), # obs-calc profiles
1589            ]
1590        Tmin = rd.powderdata[0][0]
1591        Tmax = rd.powderdata[0][-1]
1592        histname, new_names, pwdrdata = load_pwd_from_reader(rd, iparams,
1593                                                            [h.name for h in self.histograms()])
1594        if histname in self.data:
1595            print("Warning - redefining histogram", histname)
1596        elif self.names[-1][0] == 'Phases':
1597            self.names.insert(-1, new_names)
1598        else:
1599            self.names.append(new_names)
1600        if scale is not None:
1601            pwdrdata['Sample Parameters']['Scale'][0] = scale
1602        elif pwdrdata['Instrument Parameters'][0]['Type'][0].startswith('PNC'):
1603            pwdrdata['Sample Parameters']['Scale'][0] = 10000.
1604        self.data[histname] = pwdrdata
1605        self.update_ids()
1606
1607        for phase in phases:
1608            phase = self.phase(phase)
1609            self.link_histogram_phase(histname, phase)
1610
1611        return self.histogram(histname)
1612   
1613    def add_phase(self, phasefile, phasename=None, histograms=[], fmthint=None):
1614        """Loads a phase into the project from a .cif file
1615
1616        :param str phasefile: The CIF file from which to import the phase.
1617        :param str phasename: The name of the new phase, or None for the default
1618        :param list histograms: The names of the histograms to associate with
1619            this phase. Use proj.histograms() to add to all histograms.
1620        :param str fmthint: If specified, only importers where the format name
1621          (reader.formatName, as shown in Import menu) contains the
1622          supplied string will be tried as importers. If not specified, all
1623          importers consistent with the file extension will be tried
1624          (equivalent to "guess format" in menu).
1625
1626        :returns: A :class:`G2Phase` object representing the
1627            new phase.
1628        """
1629        LoadG2fil()
1630        histograms = [self.histogram(h).name for h in histograms]
1631        phasefile = os.path.abspath(os.path.expanduser(phasefile))
1632
1633        # TODO handle multiple phases in a file
1634        phasereaders = import_generic(phasefile, Readers['Phase'], fmthint=fmthint)
1635        phasereader = phasereaders[0]
1636       
1637        phasename = phasename or phasereader.Phase['General']['Name']
1638        phaseNameList = [p.name for p in self.phases()]
1639        phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
1640        phasereader.Phase['General']['Name'] = phasename
1641
1642        if 'Phases' not in self.data:
1643            self.data[u'Phases'] = { 'data': None }
1644        assert phasename not in self.data['Phases'], "phase names should be unique"
1645        self.data['Phases'][phasename] = phasereader.Phase
1646
1647        if phasereader.Constraints:
1648            Constraints = self.data['Constraints']
1649            for i in phasereader.Constraints:
1650                if isinstance(i, dict):
1651                    if '_Explain' not in Constraints:
1652                        Constraints['_Explain'] = {}
1653                    Constraints['_Explain'].update(i)
1654                else:
1655                    Constraints['Phase'].append(i)
1656
1657        data = self.data['Phases'][phasename]
1658        generalData = data['General']
1659        SGData = generalData['SGData']
1660        NShkl = len(G2spc.MustrainNames(SGData))
1661        NDij = len(G2spc.HStrainNames(SGData))
1662        Super = generalData.get('Super', 0)
1663        if Super:
1664            SuperVec = np.array(generalData['SuperVec'][0])
1665        else:
1666            SuperVec = []
1667        UseList = data['Histograms']
1668
1669        for hist in histograms:
1670            self.link_histogram_phase(hist, phasename)
1671
1672        for obj in self.names:
1673            if obj[0] == 'Phases':
1674                phasenames = obj
1675                break
1676        else:
1677            phasenames = [u'Phases']
1678            self.names.append(phasenames)
1679        phasenames.append(phasename)
1680
1681        # TODO should it be self.filename, not phasefile?
1682        SetupGeneral(data, os.path.dirname(phasefile))
1683        self.index_ids()
1684
1685        self.update_ids()
1686        return self.phase(phasename)
1687
1688    def link_histogram_phase(self, histogram, phase):
1689        """Associates a given histogram and phase.
1690
1691        .. seealso::
1692
1693            :meth:`G2Project.histogram`
1694            :meth:`G2Project.phase`"""
1695        hist = self.histogram(histogram)
1696        phase = self.phase(phase)
1697
1698        generalData = phase['General']
1699
1700        if hist.name.startswith('HKLF '):
1701            raise NotImplementedError("HKLF not yet supported")
1702        elif hist.name.startswith('PWDR '):
1703            hist['Reflection Lists'][generalData['Name']] = {}
1704            UseList = phase['Histograms']
1705            SGData = generalData['SGData']
1706            NShkl = len(G2spc.MustrainNames(SGData))
1707            NDij = len(G2spc.HStrainNames(SGData))
1708            UseList[hist.name] = SetDefaultDData('PWDR', hist.name, NShkl=NShkl, NDij=NDij)
1709            UseList[hist.name]['hId'] = hist.id
1710            for key, val in [('Use', True), ('LeBail', False),
1711                             ('newLeBail', True),
1712                             ('Babinet', {'BabA': [0.0, False],
1713                                          'BabU': [0.0, False]})]:
1714                if key not in UseList[hist.name]:
1715                    UseList[hist.name][key] = val
1716        else:
1717            raise RuntimeError("Unexpected histogram" + hist.name)
1718
1719
1720    def reload(self):
1721        """Reload self from self.filename"""
1722        data, names = LoadDictFromProjFile(self.filename)
1723        self.names = names
1724        # Need to deep copy the new data file data into the current tree,
1725        # so that any existing G2Phase, or G2PwdrData objects will still be
1726        # valid
1727        _deep_copy_into(from_=data, into=self.data)
1728
1729    def refine(self, newfile=None, printFile=None, makeBack=False):
1730        # TODO migrate to RefineCore
1731        # G2strMain.RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
1732        #      calcControls,pawleyLookup,ifPrint,printFile,dlg)
1733        # index_ids will automatically save the project
1734        self.index_ids()
1735        # TODO G2strMain does not properly use printFile
1736        G2strMain.Refine(self.filename, makeBack=makeBack)
1737        # Reload yourself
1738        self.reload()
1739
1740    def histogram(self, histname):
1741        """Returns the histogram named histname, or None if it does not exist.
1742
1743        :param histname: The name of the histogram (str), or ranId or index.
1744        :returns: A :class:`G2PwdrData` object, or None if
1745            the histogram does not exist
1746
1747        .. seealso::
1748            :meth:`G2Project.histograms`
1749            :meth:`G2Project.phase`
1750            :meth:`G2Project.phases`
1751        """
1752        if isinstance(histname, G2PwdrData):
1753            if histname.proj == self:
1754                return histname
1755        if histname in self.data:
1756            return G2PwdrData(self.data[histname], self)
1757        try:
1758            # see if histname is an id or ranId
1759            histname = int(histname)
1760        except ValueError:
1761            return
1762
1763        for histogram in self.histograms():
1764            if histogram.id == histname or histogram.ranId == histname:
1765                return histogram
1766
1767    def histograms(self, typ=None):
1768        """Return a list of all histograms, as :class:`G2PwdrData` objects
1769
1770        For now this only finds Powder/Single Xtal histograms, since that is all that is
1771        currently implemented in this module.
1772
1773        :param ste typ: The prefix (type) the histogram such as 'PWDR '. If None
1774          (the default) all known histograms types are found.
1775        :returns: a list of objects
1776
1777        .. seealso::
1778            :meth:`G2Project.histogram`
1779            :meth:`G2Project.phase`
1780            :meth:`G2Project.phases`
1781        """
1782        output = []
1783        # loop through each tree entry. If it is more than one level (more than one item in the
1784        # list of names). then it must be a histogram, unless labeled Phases or Restraints
1785        if typ is None:
1786            for obj in self.names:
1787                if obj[0].startswith('PWDR ') or obj[0].startswith('HKLF '): 
1788                    output.append(self.histogram(obj[0]))
1789        else:
1790            for obj in self.names:
1791                if len(obj) > 1 and obj[0].startswith(typ): 
1792                    output.append(self.histogram(obj[0]))
1793        return output
1794
1795    def phase(self, phasename):
1796        """
1797        Gives an object representing the specified phase in this project.
1798
1799        :param str phasename: A reference to the desired phase. Either the phase
1800            name (str), the phase's ranId, the phase's index (both int) or
1801            a phase object (:class:`G2Phase`)
1802        :returns: A :class:`G2Phase` object
1803        :raises: KeyError
1804
1805        .. seealso::
1806            :meth:`G2Project.histograms`
1807            :meth:`G2Project.phase`
1808            :meth:`G2Project.phases`
1809            """
1810        if isinstance(phasename, G2Phase):
1811            if phasename.proj == self:
1812                return phasename
1813        phases = self.data['Phases']
1814        if phasename in phases:
1815            return G2Phase(phases[phasename], phasename, self)
1816
1817        try:
1818            # phasename should be phase index or ranId
1819            phasename = int(phasename)
1820        except ValueError:
1821            return
1822
1823        for phase in self.phases():
1824            if phase.id == phasename or phase.ranId == phasename:
1825                return phase
1826
1827    def phases(self):
1828        """
1829        Returns a list of all the phases in the project.
1830
1831        :returns: A list of :class:`G2Phase` objects
1832
1833        .. seealso::
1834            :meth:`G2Project.histogram`
1835            :meth:`G2Project.histograms`
1836            :meth:`G2Project.phase`
1837            """
1838        for obj in self.names:
1839            if obj[0] == 'Phases':
1840                return [self.phase(p) for p in obj[1:]]
1841        return []
1842
1843    def _images(self):
1844        """Returns a list of all the phases in the project.
1845        """
1846        return [i[0] for i in self.names if i[0].startswith('IMG ')]
1847   
1848    def image(self, imageRef):
1849        """
1850        Gives an object representing the specified image in this project.
1851
1852        :param str imageRef: A reference to the desired image. Either the Image
1853          tree name (str), the image's index (int) or
1854          a image object (:class:`G2Image`)
1855        :returns: A :class:`G2Image` object
1856        :raises: KeyError
1857
1858        .. seealso::
1859            :meth:`G2Project.images`
1860            """
1861        if isinstance(imageRef, G2Image):
1862            if imageRef.proj == self:
1863                return imageRef
1864            else:
1865                raise Exception("Image {} not in current selected project".format(imageRef.name))
1866        if imageRef in self._images():
1867            return G2Image(self.data[imageRef], imageRef, self)
1868
1869        try:
1870            # imageRef should be an index
1871            num = int(imageRef)
1872            imageRef = self._images()[num] 
1873            return G2Image(self.data[imageRef], imageRef, self)
1874        except ValueError:
1875            raise Exception("imageRef {} not an object, name or image index in current selected project"
1876                                .format(imageRef))
1877        except IndexError:
1878            raise Exception("imageRef {} out of range (max={}) in current selected project"
1879                                .format(imageRef,len(self._images())-1))
1880           
1881    def images(self):
1882        """
1883        Returns a list of all the images in the project.
1884
1885        :returns: A list of :class:`G2Image` objects
1886        """
1887        return [G2Image(self.data[i],i,self) for i in self._images()]
1888   
1889    def update_ids(self):
1890        """Makes sure all phases and histograms have proper hId and pId"""
1891        # Translated from GetUsedHistogramsAndPhasesfromTree,
1892        #   GSASIIdataGUI.py:4107
1893        for i, h in enumerate(self.histograms()):
1894            h.id = i
1895        for i, p in enumerate(self.phases()):
1896            p.id = i
1897
1898    def do_refinements(self, refinements, histogram='all', phase='all',
1899                       outputnames=None, makeBack=False):
1900        """Conducts one or a series of refinements according to the
1901           input provided in parameter refinements. This is a wrapper
1902           around :meth:`iter_refinements`
1903
1904        :param list refinements: A list of dictionaries specifiying changes to be made to
1905            parameters before refinements are conducted.
1906            See the :ref:`Refinement_recipe` section for how this is defined.
1907        :param str histogram: Name of histogram for refinements to be applied
1908            to, or 'all'; note that this can be overridden for each refinement
1909            step via a "histograms" entry in the dict.
1910        :param str phase: Name of phase for refinements to be applied to, or
1911            'all'; note that this can be overridden for each refinement
1912            step via a "phases" entry in the dict.
1913        :param list outputnames: Provides a list of project (.gpx) file names
1914            to use for each refinement step (specifying None skips the save step).
1915            See :meth:`save`.
1916            Note that this can be overridden using an "output" entry in the dict.
1917        :param bool makeBack: determines if a backup ).bckX.gpx) file is made
1918            before a refinement is performed. The default is False.
1919           
1920        To perform a single refinement without changing any parameters, use this
1921        call:
1922
1923        .. code-block::  python
1924       
1925            my_project.do_refinements([])
1926        """
1927       
1928        for proj in self.iter_refinements(refinements, histogram, phase,
1929                                          outputnames, makeBack):
1930            pass
1931        return self
1932
1933    def iter_refinements(self, refinements, histogram='all', phase='all',
1934                         outputnames=None, makeBack=False):
1935        """Conducts a series of refinements, iteratively. Stops after every
1936        refinement and yields this project, to allow error checking or
1937        logging of intermediate results. Parameter use is the same as for
1938        :meth:`do_refinements` (which calls this method).
1939
1940        >>> def checked_refinements(proj):
1941        ...     for p in proj.iter_refinements(refs):
1942        ...         # Track intermediate results
1943        ...         log(p.histogram('0').residuals)
1944        ...         log(p.phase('0').get_cell())
1945        ...         # Check if parameter diverged, nonsense answer, or whatever
1946        ...         if is_something_wrong(p):
1947        ...             raise Exception("I need a human!")
1948
1949           
1950        """
1951        if outputnames:
1952            if len(refinements) != len(outputnames):
1953                raise ValueError("Should have same number of outputs to"
1954                                 "refinements")
1955        else:
1956            outputnames = [None for r in refinements]
1957
1958        for output, refinedict in zip(outputnames, refinements):
1959            if 'histograms' in refinedict:
1960                hist = refinedict['histograms']
1961            else:
1962                hist = histogram
1963            if 'phases' in refinedict:
1964                ph = refinedict['phases']
1965            else:
1966                ph = phase
1967            if 'output' in refinedict:
1968                output = refinedict['output']
1969            self.set_refinement(refinedict, hist, ph)
1970            # Handle 'once' args - refinements that are disabled after this
1971            # refinement
1972            if 'once' in refinedict:
1973                temp = {'set': refinedict['once']}
1974                self.set_refinement(temp, hist, ph)
1975
1976            if output:
1977                self.save(output)
1978
1979            if 'skip' not in refinedict:
1980                self.refine(makeBack=makeBack)
1981            yield self
1982
1983            # Handle 'once' args - refinements that are disabled after this
1984            # refinement
1985            if 'once' in refinedict:
1986                temp = {'clear': refinedict['once']}
1987                self.set_refinement(temp, hist, ph)
1988            if 'call' in refinedict:
1989                fxn = refinedict['call']
1990                if callable(fxn):
1991                    fxn(*refinedict.get('callargs',[self]))
1992                elif callable(eval(fxn)):
1993                    eval(fxn)(*refinedict.get('callargs',[self]))
1994                else:
1995                    raise G2ScriptException("Error: call value {} is not callable".format(fxn))
1996
1997    def set_refinement(self, refinement, histogram='all', phase='all'):
1998        """Apply specified refinements to a given histogram(s) or phase(s).
1999
2000        :param dict refinement: The refinements to be conducted
2001        :param histogram: Specifies either 'all' (default), a single histogram or
2002          a list of histograms. Histograms may be specified as histogram objects
2003          (see :class:`G2PwdrData`), the histogram name (str) or the index number (int)
2004          of the histogram in the project, numbered starting from 0.
2005          Omitting the parameter or the string 'all' indicates that parameters in
2006          all histograms should be set.
2007        :param phase: Specifies either 'all' (default), a single phase or
2008          a list of phases. Phases may be specified as phase objects
2009          (see :class:`G2Phase`), the phase name (str) or the index number (int)
2010          of the phase in the project, numbered starting from 0.
2011          Omitting the parameter or the string 'all' indicates that parameters in
2012          all phases should be set.
2013
2014        Note that refinement parameters are categorized as one of three types:
2015
2016        1. Histogram parameters
2017        2. Phase parameters
2018        3. Histogram-and-Phase (HAP) parameters
2019       
2020        .. seealso::
2021            :meth:`G2PwdrData.set_refinements`
2022            :meth:`G2PwdrData.clear_refinements`
2023            :meth:`G2Phase.set_refinements`
2024            :meth:`G2Phase.clear_refinements`
2025            :meth:`G2Phase.set_HAP_refinements`
2026            :meth:`G2Phase.clear_HAP_refinements`"""
2027
2028        if histogram == 'all':
2029            hists = self.histograms()
2030        elif isinstance(histogram, list) or isinstance(histogram, tuple):
2031            hists = []
2032            for h in histogram:
2033                if isinstance(h, str) or isinstance(h, int):
2034                    hists.append(self.histogram(h))
2035                else:
2036                    hists.append(h)
2037        elif isinstance(histogram, str) or isinstance(histogram, int):
2038            hists = [self.histogram(histogram)]
2039        else:
2040            hists = [histogram]
2041
2042        if phase == 'all':
2043            phases = self.phases()
2044        elif isinstance(phase, list) or isinstance(phase, tuple):
2045            phases = []
2046            for ph in phase:
2047                if isinstance(ph, str) or isinstance(ph, int):
2048                    phases.append(self.phase(ph))
2049                else:
2050                    phases.append(ph)
2051        elif isinstance(phase, str) or isinstance(phase, int):
2052            phases = [self.phase(phase)]
2053        else:
2054            phases = [phase]
2055
2056        # TODO: HAP parameters:
2057        #   Babinet
2058        #   Extinction
2059        #   HStrain
2060        #   Mustrain
2061        #   Pref. Ori
2062        #   Size
2063
2064        pwdr_set = {}
2065        phase_set = {}
2066        hap_set = {}
2067        for key, val in refinement.get('set', {}).items():
2068            # Apply refinement options
2069            if G2PwdrData.is_valid_refinement_key(key):
2070                pwdr_set[key] = val
2071            elif G2Phase.is_valid_refinement_key(key):
2072                phase_set[key] = val
2073            elif G2Phase.is_valid_HAP_refinement_key(key):
2074                hap_set[key] = val
2075            else:
2076                raise ValueError("Unknown refinement key", key)
2077
2078        for hist in hists:
2079            hist.set_refinements(pwdr_set)
2080        for phase in phases:
2081            phase.set_refinements(phase_set)
2082        for phase in phases:
2083            phase.set_HAP_refinements(hap_set, hists)
2084
2085        pwdr_clear = {}
2086        phase_clear = {}
2087        hap_clear = {}
2088        for key, val in refinement.get('clear', {}).items():
2089            # Clear refinement options
2090            if G2PwdrData.is_valid_refinement_key(key):
2091                pwdr_clear[key] = val
2092            elif G2Phase.is_valid_refinement_key(key):
2093                phase_clear[key] = val
2094            elif G2Phase.is_valid_HAP_refinement_key(key):
2095                hap_set[key] = val
2096            else:
2097                raise ValueError("Unknown refinement key", key)
2098
2099        for hist in hists:
2100            hist.clear_refinements(pwdr_clear)
2101        for phase in phases:
2102            phase.clear_refinements(phase_clear)
2103        for phase in phases:
2104            phase.clear_HAP_refinements(hap_clear, hists)
2105
2106    def index_ids(self):
2107        import GSASIIstrIO as G2strIO
2108        self.save()
2109        return G2strIO.GetUsedHistogramsAndPhases(self.filename)
2110
2111    def add_constraint_raw(self, cons_scope, constr):
2112        """Adds a constraint of type consType to the project.
2113        cons_scope should be one of "Hist", "Phase", "HAP", or "Global".
2114
2115        WARNING it does not check the constraint is well-constructed"""
2116        constrs = self.data['Constraints']['data']
2117        if 'Global' not in constrs:
2118            constrs['Global'] = []
2119        constrs[cons_scope].append(constr)
2120
2121    def hold_many(self, vars, type):
2122        """Apply holds for all the variables in vars, for constraint of a given type.
2123
2124        type is passed directly to add_constraint_raw as consType
2125
2126        :param list vars: A list of variables to hold. Either :class:`GSASIIobj.G2VarObj` objects,
2127            string variable specifiers, or arguments for :meth:`make_var_obj`
2128        :param str type: A string constraint type specifier. See
2129            :class:`G2Project.add_constraint_raw`
2130
2131        """
2132        for var in vars:
2133            if isinstance(var, str):
2134                var = self.make_var_obj(var)
2135            elif not isinstance(var, G2obj.G2VarObj):
2136                var = self.make_var_obj(*var)
2137            self.add_constraint_raw(type, [[1.0, var], None, None, 'h'])
2138
2139    def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
2140                     reloadIdx=True):
2141        """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a")
2142        or individual names of phase, histogram, varname, and atomId.
2143
2144        Automatically converts string phase, hist, or atom names into the ID required
2145        by G2VarObj."""
2146
2147        if reloadIdx:
2148            self.index_ids()
2149
2150        # If string representation, short circuit
2151        if hist is None and varname is None and atomId is None:
2152            if isinstance(phase, str) and ':' in phase:
2153                return G2obj.G2VarObj(phase)
2154
2155        # Get phase index
2156        phaseObj = None
2157        if isinstance(phase, G2Phase):
2158            phaseObj = phase
2159            phase = G2obj.PhaseRanIdLookup[phase.ranId]
2160        elif phase in self.data['Phases']:
2161            phaseObj = self.phase(phase)
2162            phaseRanId = phaseObj.ranId
2163            phase = G2obj.PhaseRanIdLookup[phaseRanId]
2164
2165        # Get histogram index
2166        if isinstance(hist, G2PwdrData):
2167            hist = G2obj.HistRanIdLookup[hist.ranId]
2168        elif hist in self.data:
2169            histRanId = self.histogram(hist).ranId
2170            hist = G2obj.HistRanIdLookup[histRanId]
2171
2172        # Get atom index (if any)
2173        if isinstance(atomId, G2AtomRecord):
2174            atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
2175        elif phaseObj:
2176            atomObj = phaseObj.atom(atomId)
2177            if atomObj:
2178                atomRanId = atomObj.ranId
2179                atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
2180
2181        return G2obj.G2VarObj(phase, hist, varname, atomId)
2182
2183    def add_image(self, imagefile, fmthint=None, defaultImage=None):
2184        """Load an image into a project
2185
2186        :param str imagefile: The image file to read, a filename.
2187        :param str fmthint: If specified, only importers where the format name
2188          (reader.formatName, as shown in Import menu) contains the
2189          supplied string will be tried as importers. If not specified, all
2190          importers consistent with the file extension will be tried
2191          (equivalent to "guess format" in menu).
2192        :param str defaultImage: The name of an image to use as a default for
2193          setting parameters for the image file to read.
2194
2195        :returns: a list of G2Image object for the added image(s) [this routine
2196          has not yet been tested with files with multiple images...]
2197        """
2198        LoadG2fil()
2199        imagefile = os.path.abspath(os.path.expanduser(imagefile))
2200        readers = import_generic(imagefile, Readers['Image'], fmthint=fmthint)
2201        objlist = []
2202        for rd in readers:
2203            if rd.SciPy:        #was default read by scipy; needs 1 time fixes
2204                print('TODO: Image {} read by SciPy. Parameters likely need editing'.format(imagefile))
2205                #see this: G2IO.EditImageParms(self,rd.Data,rd.Comments,rd.Image,imagefile)
2206                rd.SciPy = False
2207            rd.readfilename = imagefile
2208            if rd.repeatcount == 1 and not rd.repeat: # skip image number if only one in set
2209                rd.Data['ImageTag'] = None
2210            else:
2211                rd.Data['ImageTag'] = rd.repeatcount
2212            rd.Data['formatName'] = rd.formatName
2213            if rd.sumfile:
2214                rd.readfilename = rd.sumfile
2215            # Load generic metadata, as configured
2216            G2fil.GetColumnMetadata(rd)
2217            # Code from G2IO.LoadImage2Tree(rd.readfilename,self,rd.Comments,rd.Data,rd.Npix,rd.Image)
2218            Imax = np.amax(rd.Image)
2219            ImgNames = [i[0] for i in self.names if i[0].startswith('IMG ')]
2220            TreeLbl = 'IMG '+os.path.basename(imagefile)
2221            ImageTag = rd.Data.get('ImageTag')
2222            if ImageTag:
2223                TreeLbl += ' #'+'%04d'%(ImageTag)
2224                imageInfo = (imagefile,ImageTag)
2225            else:
2226                imageInfo = imagefile
2227            TreeName = G2obj.MakeUniqueLabel(TreeLbl,ImgNames)
2228            # MT dict to contain image info
2229            ImgDict = {}
2230            ImgDict['data'] = [rd.Npix,imageInfo]
2231            ImgDict['Comments'] = rd.Comments
2232            if defaultImage:
2233                if isinstance(defaultImage, G2Image):
2234                    if defaultImage.proj == self:
2235                        defaultImage = self.data[defaultImage.name]['data']
2236                    else:
2237                        raise Exception("Image {} not in current selected project".format(defaultImage.name))
2238                elif defaultImage in self._images():
2239                    defaultImage = self.data[defaultImage]['data']
2240                else:
2241                    defaultImage = None
2242            Data = rd.Data
2243            if defaultImage:
2244                Data = copy.copy(defaultImage)
2245                Data['showLines'] = True
2246                Data['ring'] = []
2247                Data['rings'] = []
2248                Data['cutoff'] = 10.
2249                Data['pixLimit'] = 20
2250                Data['edgemin'] = 100000000
2251                Data['calibdmin'] = 0.5
2252                Data['calibskip'] = 0
2253                Data['ellipses'] = []
2254                Data['calibrant'] = ''
2255                Data['GonioAngles'] = [0.,0.,0.]
2256                Data['DetDepthRef'] = False
2257            else:
2258                Data['type'] = 'PWDR'
2259                Data['color'] = GSASIIpath.GetConfigValue('Contour_color','Paired')
2260                if 'tilt' not in Data:          #defaults if not preset in e.g. Bruker importer
2261                    Data['tilt'] = 0.0
2262                    Data['rotation'] = 0.0
2263                    Data['pixLimit'] = 20
2264                    Data['calibdmin'] = 0.5
2265                    Data['cutoff'] = 10.
2266                Data['showLines'] = False
2267                Data['calibskip'] = 0
2268                Data['ring'] = []
2269                Data['rings'] = []
2270                Data['edgemin'] = 100000000
2271                Data['ellipses'] = []
2272                Data['GonioAngles'] = [0.,0.,0.]
2273                Data['DetDepth'] = 0.
2274                Data['DetDepthRef'] = False
2275                Data['calibrant'] = ''
2276                Data['IOtth'] = [5.0,50.0]
2277                Data['LRazimuth'] = [0.,180.]
2278                Data['azmthOff'] = 0.0
2279                Data['outChannels'] = 2250
2280                Data['outAzimuths'] = 1
2281                Data['centerAzm'] = False
2282                Data['fullIntegrate'] = GSASIIpath.GetConfigValue('fullIntegrate',True)
2283                Data['setRings'] = False
2284                Data['background image'] = ['',-1.0]                           
2285                Data['dark image'] = ['',-1.0]
2286                Data['Flat Bkg'] = 0.0
2287                Data['Oblique'] = [0.5,False]
2288            Data['setDefault'] = False
2289            Data['range'] = [(0,Imax),[0,Imax]]
2290            ImgDict['Image Controls'] = Data
2291            ImgDict['Masks'] = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],
2292                                'Frames':[],'Thresholds':[(0,Imax),[0,Imax]]}
2293            ImgDict['Stress/Strain']  = {'Type':'True','d-zero':[],'Sample phi':0.0,
2294                                             'Sample z':0.0,'Sample load':0.0}
2295            self.names.append([TreeName]+['Comments','Image Controls','Masks','Stress/Strain'])
2296            self.data[TreeName] = ImgDict
2297            del rd.Image
2298            objlist.append(G2Image(self.data[TreeName], TreeName, self))
2299        return objlist
2300
2301class G2AtomRecord(G2ObjectWrapper):
2302    """Wrapper for an atom record. Has convenient accessors via @property.
2303
2304
2305    Available accessors: label, type, refinement_flags, coordinates,
2306        occupancy, ranId, id, adp_flag, uiso
2307
2308    Example:
2309
2310    >>> atom = some_phase.atom("O3")
2311    >>> # We can access the underlying data format
2312    >>> atom.data
2313    ['O3', 'O-2', '', ... ]
2314    >>> # We can also use wrapper accessors
2315    >>> atom.coordinates
2316    (0.33, 0.15, 0.5)
2317    >>> atom.refinement_flags
2318    u'FX'
2319    >>> atom.ranId
2320    4615973324315876477
2321    >>> atom.occupancy
2322    1.0
2323    """
2324    def __init__(self, data, indices, proj):
2325        self.data = data
2326        self.cx, self.ct, self.cs, self.cia = indices
2327        self.proj = proj
2328
2329    @property
2330    def label(self):
2331        '''Get the associated atom's label
2332        '''
2333        return self.data[self.ct-1]
2334
2335    @property
2336    def type(self):
2337        '''Get the associated atom's type
2338        '''
2339        return self.data[self.ct]
2340
2341    @property
2342    def refinement_flags(self):
2343        '''Get or set refinement flags for the associated atom
2344        '''
2345        return self.data[self.ct+1]
2346
2347    @refinement_flags.setter
2348    def refinement_flags(self, other):
2349        # Automatically check it is a valid refinement
2350        for c in other:
2351            if c not in ' FXU':
2352                raise ValueError("Invalid atom refinement: ", other)
2353        self.data[self.ct+1] = other
2354
2355    @property
2356    def coordinates(self):
2357        '''Get the associated atom's coordinates
2358        '''
2359        return tuple(self.data[self.cx:self.cx+3])
2360
2361    @property
2362    def occupancy(self):
2363        '''Get or set the associated atom's occupancy fraction
2364        '''
2365        return self.data[self.cx+3]
2366
2367    @occupancy.setter
2368    def occupancy(self, val):
2369        self.data[self.cx+3] = float(val)
2370
2371    @property
2372    def ranId(self):
2373        '''Get the associated atom's Random Id number
2374        '''
2375        return self.data[self.cia+8]
2376
2377    @property
2378    def adp_flag(self):
2379        '''Get the associated atom's iso/aniso setting, 'I' or 'A'
2380        '''
2381        # Either 'I' or 'A'
2382        return self.data[self.cia]
2383
2384    @property
2385    def uiso(self):
2386        '''Get or set the associated atom's Uiso or Uaniso value(s)
2387        '''
2388        if self.adp_flag == 'I':
2389            return self.data[self.cia+1]
2390        else:
2391            return self.data[self.cia+2:self.cia+8]
2392
2393    @uiso.setter
2394    def uiso(self, value):
2395        if self.adp_flag == 'I':
2396            self.data[self.cia+1] = float(value)
2397        else:
2398            assert len(value) == 6
2399            self.data[self.cia+2:self.cia+8] = [float(v) for v in value]
2400
2401
2402class G2PwdrData(G2ObjectWrapper):
2403    """Wraps a Powder Data Histogram."""
2404    def __init__(self, data, proj):
2405        self.data = data
2406        self.proj = proj
2407
2408    @staticmethod
2409    def is_valid_refinement_key(key):
2410        valid_keys = ['Limits', 'Sample Parameters', 'Background',
2411                      'Instrument Parameters']
2412        return key in valid_keys
2413
2414    @property
2415    def name(self):
2416        return self.data['data'][-1]
2417
2418    @property
2419    def ranId(self):
2420        return self.data['data'][0]['ranId']
2421
2422    @property
2423    def residuals(self):
2424        '''Provides a dictionary with with the R-factors for this histogram.
2425        Includes the weighted and unweighted profile terms (R, Rb, wR, wRb, wRmin)
2426        as well as the Bragg R-values for each phase (ph:H:Rf and ph:H:Rf^2).
2427        '''
2428        data = self.data['data'][0]
2429        return {key: data[key] for key in data
2430                if key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']
2431                   or ':' in key}
2432
2433    @property
2434    def InstrumentParameters(self):
2435        '''Provides a dictionary with with the Instrument Parameters
2436        for this histogram.
2437        '''
2438        return self.data['Instrument Parameters'][0]
2439
2440    @property
2441    def SampleParameters(self):
2442        '''Provides a dictionary with with the Sample Parameters
2443        for this histogram.
2444        '''
2445        return self.data['Sample Parameters']
2446
2447    @property
2448    def Background(self):
2449        '''Provides a list with with the Background parameters
2450        for this histogram.
2451
2452        :returns: list containing a list and dict with background values
2453        '''
2454        return self.data['Background']
2455
2456    def add_back_peak(self,pos,int,sig,gam,refflags=[]):
2457        '''Adds a background peak to the Background parameters
2458       
2459        :param float pos: position of peak, a 2theta or TOF value
2460        :param float int: integrated intensity of background peak, usually large
2461        :param float sig: Gaussian width of background peak, usually large
2462        :param float gam: Lorentzian width of background peak, usually unused (small)
2463        :param list refflags: a list of 1 to 4 boolean refinement flags for
2464            pos,int,sig & gam, respectively (use [0,1] to refine int only).
2465            Defaults to [] which means nothing is refined.
2466        '''
2467        if 'peaksList' not in self.Background[1]:
2468            self.Background[1]['peaksList'] = []
2469        flags = 4*[False]
2470        for i,f in enumerate(refflags):
2471            if i>3: break
2472            flags[i] = bool(f)
2473        bpk = []
2474        for i,j in zip((pos,int,sig,gam),flags):
2475            bpk += [float(i),j]
2476        self.Background[1]['peaksList'].append(bpk)
2477        self.Background[1]['nPeaks'] = len(self.Background[1]['peaksList'])
2478
2479    def del_back_peak(self,peaknum):
2480        '''Removes a background peak from the Background parameters
2481       
2482        :param int peaknum: the number of the peak (starting from 0)
2483        '''
2484        npks = self.Background[1].get('nPeaks',0)
2485        if peaknum >= npks:
2486            raise Exception('peak {} not found in histogram {}'.format(peaknum,self.name))
2487        del self.Background[1]['peaksList'][peaknum]
2488        self.Background[1]['nPeaks'] = len(self.Background[1]['peaksList'])
2489       
2490    def ref_back_peak(self,peaknum,refflags=[]):
2491        '''Sets refinement flag for a background peak
2492       
2493        :param int peaknum: the number of the peak (starting from 0)
2494        :param list refflags: a list of 1 to 4 boolean refinement flags for
2495            pos,int,sig & gam, respectively. If a flag is not specified
2496            it defaults to False (use [0,1] to refine int only).
2497            Defaults to [] which means nothing is refined.
2498        '''
2499        npks = self.Background[1].get('nPeaks',0)
2500        if peaknum >= npks:
2501            raise Exception('peak {} not found in histogram {}'.format(peaknum,self.name))
2502        flags = 4*[False]
2503        for i,f in enumerate(refflags):
2504            if i>3: break
2505            flags[i] = bool(f)
2506        for i,f in enumerate(flags):
2507            self.Background[1]['peaksList'][peaknum][2*i+1] = f
2508                   
2509    @property
2510    def id(self):
2511        self.proj.update_ids()
2512        return self.data['data'][0]['hId']
2513
2514    @id.setter
2515    def id(self, val):
2516        self.data['data'][0]['hId'] = val
2517
2518    def fit_fixed_points(self):
2519        """Attempts to apply a background fit to the fixed points currently specified."""
2520        def SetInstParms(Inst):
2521            dataType = Inst['Type'][0]
2522            insVary = []
2523            insNames = []
2524            insVals = []
2525            for parm in Inst:
2526                insNames.append(parm)
2527                insVals.append(Inst[parm][1])
2528                if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
2529                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
2530                        Inst[parm][2] = False
2531            instDict = dict(zip(insNames, insVals))
2532            instDict['X'] = max(instDict['X'], 0.01)
2533            instDict['Y'] = max(instDict['Y'], 0.01)
2534            if 'SH/L' in instDict:
2535                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
2536            return dataType, instDict, insVary
2537
2538        bgrnd = self.data['Background']
2539
2540        # Need our fixed points in order
2541        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
2542        X = [x for x, y in bgrnd[1]['FixedPoints']]
2543        Y = [y for x, y in bgrnd[1]['FixedPoints']]
2544
2545        limits = self.data['Limits'][1]
2546        if X[0] > limits[0]:
2547            X = [limits[0]] + X
2548            Y = [Y[0]] + Y
2549        if X[-1] < limits[1]:
2550            X += [limits[1]]
2551            Y += [Y[-1]]
2552
2553        # Some simple lookups
2554        controls = self.proj['Controls']['data']
2555        inst, inst2 = self.data['Instrument Parameters']
2556        pwddata = self.data['data'][1]
2557
2558        # Construct the data for background fitting
2559        xBeg = np.searchsorted(pwddata[0], limits[0])
2560        xFin = np.searchsorted(pwddata[0], limits[1])
2561        xdata = pwddata[0][xBeg:xFin]
2562        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
2563
2564        W = [1]*len(xdata)
2565        Z = [0]*len(xdata)
2566
2567        dataType, insDict, insVary = SetInstParms(inst)
2568        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
2569
2570        # Do the fit
2571        data = np.array([xdata, ydata, W, Z, Z, Z])
2572        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
2573                        prevVaryList=bakVary, controls=controls)
2574
2575        # Post-fit
2576        parmDict = {}
2577        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
2578        parmDict.update(bakDict)
2579        parmDict.update(insDict)
2580        pwddata[3][xBeg:xFin] *= 0
2581        pwddata[5][xBeg:xFin] *= 0
2582        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
2583
2584        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
2585        # TODO update background
2586        self.proj.save()
2587
2588    def getdata(self,datatype):
2589        '''Provides access to the histogram data of the selected data type
2590
2591        :param str datatype: must be one of the following values (case is ignored)
2592       
2593           * 'X': the 2theta or TOF values for the pattern
2594           * 'Yobs': the observed intensity values
2595           * 'Yweight': the weights for each data point (1/sigma**2)
2596           * 'Ycalc': the computed intensity values
2597           * 'Background': the computed background values
2598           * 'Residual': the difference between Yobs and Ycalc (obs-calc)
2599
2600        :returns: an numpy MaskedArray with data values of the requested type
2601       
2602        '''
2603        enums = ['x', 'yobs', 'yweight', 'ycalc', 'background', 'residual']
2604        if datatype.lower() not in enums:
2605            raise G2ScriptException("Invalid datatype = "+datatype+" must be one of "+str(enums))
2606        return self.data['data'][1][enums.index(datatype.lower())]
2607       
2608    def y_calc(self):
2609        return self.data['data'][1][3]
2610
2611    def Export(self,fileroot,extension):
2612        '''Write the histogram into a file. The path is specified by fileroot and
2613        extension.
2614       
2615        :param str fileroot: name of the file, optionally with a path (extension is
2616           ignored)
2617        :param str extension: includes '.', must match an extension in global
2618           exportersByExtension['powder'] or a Exception is raised.
2619        :returns: name of file that was written
2620        '''
2621        if extension not in exportersByExtension.get('powder',[]):
2622            raise G2ScriptException('No Writer for file type = "'+extension+'"')
2623        fil = os.path.abspath(os.path.splitext(fileroot)[0]+extension)
2624        obj = exportersByExtension['powder'][extension]
2625        obj.SetFromArray(hist=self.data,histname=self.name)
2626        obj.Writer(self.name,fil)
2627           
2628    def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
2629        try:
2630            import matplotlib.pyplot as plt
2631            data = self.data['data'][1]
2632            if Yobs:
2633                plt.plot(data[0], data[1], label='Yobs')
2634            if Ycalc:
2635                plt.plot(data[0], data[3], label='Ycalc')
2636            if Background:
2637                plt.plot(data[0], data[4], label='Background')
2638            if Residual:
2639                plt.plot(data[0], data[5], label="Residual")
2640        except ImportError:
2641            pass
2642
2643    def get_wR(self):
2644        """returns the overall weighted profile R factor for a histogram
2645       
2646        :returns: a wR value as a percentage or None if not defined
2647        """
2648        return self['data'][0].get('wR')
2649
2650    def set_refinements(self, refs):
2651        """Sets the histogram refinement parameter 'key' to the specification 'value'
2652
2653        :param dict refs: A dictionary of the parameters to be set. See
2654                          :ref:`Histogram_parameters_table` for a description of
2655                          what these dictionaries should be.
2656
2657        :returns: None
2658
2659        """
2660        do_fit_fixed_points = False
2661        for key, value in refs.items():
2662            if key == 'Limits':
2663                old_limits = self.data['Limits'][1]
2664                new_limits = value
2665                if isinstance(new_limits, dict):
2666                    if 'low' in new_limits:
2667                        old_limits[0] = new_limits['low']
2668                    if 'high' in new_limits:
2669                        old_limits[1] = new_limits['high']
2670                else:
2671                    old_limits[0], old_limits[1] = new_limits
2672            elif key == 'Sample Parameters':
2673                sample = self.data['Sample Parameters']
2674                for sparam in value:
2675                    if sparam not in sample:
2676                        raise ValueError("Unknown refinement parameter, "
2677                                         + str(sparam))
2678                    sample[sparam][1] = True
2679            elif key == 'Background':
2680                bkg, peaks = self.data['Background']
2681
2682                # If True or False, just set the refine parameter
2683                if value in (True, False):
2684                    bkg[1] = value
2685                    return
2686
2687                if 'type' in value:
2688                    bkg[0] = value['type']
2689                if 'refine' in value:
2690                    bkg[1] = value['refine']
2691                if 'no. coeffs' in value:
2692                    cur_coeffs = bkg[2]
2693                    n_coeffs = value['no. coeffs']
2694                    if n_coeffs > cur_coeffs:
2695                        for x in range(n_coeffs - cur_coeffs):
2696                            bkg.append(0.0)
2697                    else:
2698                        for _ in range(cur_coeffs - n_coeffs):
2699                            bkg.pop()
2700                    bkg[2] = n_coeffs
2701                if 'coeffs' in value:
2702                    bkg[3:] = value['coeffs']
2703                if 'FixedPoints' in value:
2704                    peaks['FixedPoints'] = [(float(a), float(b))
2705                                            for a, b in value['FixedPoints']]
2706                if value.get('fit fixed points', False):
2707                    do_fit_fixed_points = True
2708                if 'peaks' in value:
2709                    for i,flags in enumerate(value['peaks']):
2710                        self.ref_back_peak(i,flags)
2711
2712            elif key == 'Instrument Parameters':
2713                instrument, secondary = self.data['Instrument Parameters']
2714                for iparam in value:
2715                    try:
2716                        instrument[iparam][2] = True
2717                    except IndexError:
2718                        raise ValueError("Invalid key:", iparam)
2719            else:
2720                raise ValueError("Unknown key:", key)
2721        # Fit fixed points after the fact - ensure they are after fixed points
2722        # are added
2723        if do_fit_fixed_points:
2724            # Background won't be fit if refinement flag not set
2725            orig = self.data['Background'][0][1]
2726            self.data['Background'][0][1] = True
2727            self.fit_fixed_points()
2728            # Restore the previous value
2729            self.data['Background'][0][1] = orig
2730
2731    def clear_refinements(self, refs):
2732        """Clears the refinement parameter 'key' and its associated value.
2733
2734        :param dict refs: A dictionary of parameters to clear."""
2735        for key, value in refs.items():
2736            if key == 'Limits':
2737                old_limits, cur_limits = self.data['Limits']
2738                cur_limits[0], cur_limits[1] = old_limits
2739            elif key == 'Sample Parameters':
2740                sample = self.data['Sample Parameters']
2741                for sparam in value:
2742                    sample[sparam][1] = False
2743            elif key == 'Background':
2744                bkg, peaks = self.data['Background']
2745
2746                # If True or False, just set the refine parameter
2747                if value in (True, False):
2748                    bkg[1] = False
2749                    return
2750
2751                bkg[1] = False
2752                if 'FixedPoints' in value:
2753                    if 'FixedPoints' in peaks:
2754                        del peaks['FixedPoints']
2755                if 'peaks' in value:
2756                    for i in range(len(self.Background[1]['peaksList'])):
2757                        self.ref_back_peak(i,[])
2758            elif key == 'Instrument Parameters':
2759                instrument, secondary = self.data['Instrument Parameters']
2760                for iparam in value:
2761                    instrument[iparam][2] = False
2762            else:
2763                raise ValueError("Unknown key:", key)
2764
2765    def add_peak(self,area,dspace=None,Q=None,ttheta=None):
2766        '''Adds a single peak to the peak list
2767        :param float area: peak area
2768        :param float dspace: peak position as d-space (A)
2769        :param float Q: peak position as Q (A-1)
2770        :param float ttheta: peak position as 2Theta (deg)
2771
2772        Note: only one of the parameters dspace, Q or ttheta may be specified
2773        '''
2774        import GSASIIlattice as G2lat
2775        import GSASIImath as G2mth
2776        if (not dspace) + (not Q) + (not ttheta) != 2:
2777            print('add_peak error: too many or no peak position(s) specified')
2778            return
2779        pos = ttheta
2780        Parms,Parms2 = self.data['Instrument Parameters']
2781        if Q:
2782            pos = G2lat.Dsp2pos(Parms,2.*np.pi/Q)
2783        elif dspace:
2784            pos = G2lat.Dsp2pos(Parms,dspace)
2785        peaks = self.data['Peak List']
2786        peaks['sigDict'] = {}        #no longer valid
2787        peaks['peaks'].append(G2mth.setPeakparms(Parms,Parms2,pos,area))
2788
2789    def set_peakFlags(self,peaklist=None,area=None,pos=None,sig=None,gam=None):
2790        '''Set refinement flags for peaks
2791       
2792        :param list peaklist: a list of peaks to change flags. If None (default), changes
2793          are made to all peaks.
2794        :param bool area: Sets or clears the refinement flag for the peak area value.
2795          If None (the default), no change is made.
2796        :param bool pos: Sets or clears the refinement flag for the peak position value.
2797          If None (the default), no change is made.
2798        :param bool sig: Sets or clears the refinement flag for the peak sig (Gaussian width) value.
2799          If None (the default), no change is made.
2800        :param bool gam: Sets or clears the refinement flag for the peak sig (Lorentzian width) value.
2801          If None (the default), no change is made.
2802         
2803        Note that when peaks are first created the area flag is on and the other flags are
2804        initially off.
2805
2806        Example::
2807       
2808           set_peakFlags(sig=False,gam=True)
2809
2810        causes the sig refinement flag to be cleared and the gam flag to be set, in both cases for
2811        all peaks. The position and area flags are not changed from their previous values.
2812        '''
2813        peaks = self.data['Peak List']
2814        if peaklist is None:
2815            peaklist = range(len(peaks['peaks']))
2816        for i in peaklist:
2817            for var,j in [(area,3),(pos,1),(sig,5),(gam,7)]:
2818                if var is not None:
2819                    peaks['peaks'][i][j] = var
2820           
2821    def refine_peaks(self):
2822        '''Causes a refinement of peak position, background and instrument parameters
2823        '''
2824        import GSASIIpwd as G2pwd
2825        controls = self.proj.data.get('Controls',{})
2826        controls = controls.get('data',
2827                            {'deriv type':'analytic','min dM/M':0.001,}     #fill in defaults if needed
2828                            )
2829        peaks = self.data['Peak List']
2830        Parms,Parms2 = self.data['Instrument Parameters']
2831        background = self.data['Background']
2832        limits = self.data['Limits'][1]
2833        bxye = np.zeros(len(self.data['data'][1][1]))
2834        peaks['sigDict'] = G2pwd.DoPeakFit('LSQ',peaks['peaks'],background,limits,
2835                                           Parms,Parms2,self.data['data'][1],bxye,[],
2836                                           False,controls,None)[0]
2837
2838    def SaveProfile(self,filename):
2839        '''Writes a GSAS-II (new style) .instprm file
2840        '''
2841        data,Parms2 = self.data['Instrument Parameters']
2842        filename = os.path.splitext(filename)[0]+'.instprm'         # make sure extension is .instprm
2843        File = open(filename,'w')
2844        File.write("#GSAS-II instrument parameter file; do not add/delete items!\n")
2845        for item in data:
2846            File.write(item+':'+str(data[item][1])+'\n')
2847        File.close()
2848        print ('Instrument parameters saved to: '+filename)
2849
2850    def LoadProfile(self,filename):
2851        '''Reads a GSAS-II (new style) .instprm file and overwrites the current parameters
2852        '''
2853        filename = os.path.splitext(filename)[0]+'.instprm'         # make sure extension is .instprm
2854        File = open(filename,'r')
2855        S = File.readline()
2856        newItems = []
2857        newVals = []
2858        Found = False
2859        while S:
2860            if S[0] == '#':
2861                if Found:
2862                    break
2863                if 'Bank' in S:
2864                    if bank == int(S.split(':')[0].split()[1]):
2865                        S = File.readline()
2866                        continue
2867                    else:
2868                        S = File.readline()
2869                        while S and '#Bank' not in S:
2870                            S = File.readline()
2871                        continue
2872                else:   #a non #Bank file
2873                    S = File.readline()
2874                    continue
2875            Found = True
2876            [item,val] = S[:-1].split(':')
2877            newItems.append(item)
2878            try:
2879                newVals.append(float(val))
2880            except ValueError:
2881                newVals.append(val)                       
2882            S = File.readline()               
2883        File.close()
2884        LoadG2fil()
2885        self.data['Instrument Parameters'][0] = G2fil.makeInstDict(newItems,newVals,len(newVals)*[False,])
2886
2887       
2888class G2Phase(G2ObjectWrapper):
2889    """A wrapper object around a given phase.
2890
2891    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
2892    """
2893    def __init__(self, data, name, proj):
2894        self.data = data
2895        self.name = name
2896        self.proj = proj
2897
2898    @staticmethod
2899    def is_valid_refinement_key(key):
2900        valid_keys = ["Cell", "Atoms", "LeBail"]
2901        return key in valid_keys
2902
2903    @staticmethod
2904    def is_valid_HAP_refinement_key(key):
2905        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
2906                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
2907        return key in valid_keys
2908
2909    def atom(self, atomlabel):
2910        """Returns the atom specified by atomlabel, or None if it does not
2911        exist.
2912
2913        :param str atomlabel: The name of the atom (e.g. "O2")
2914        :returns: A :class:`G2AtomRecord` object
2915            representing the atom.
2916        """
2917        # Consult GSASIIobj.py for the meaning of this
2918        cx, ct, cs, cia = self.data['General']['AtomPtrs']
2919        ptrs = [cx, ct, cs, cia]
2920        atoms = self.data['Atoms']
2921        for atom in atoms:
2922            if atom[ct-1] == atomlabel:
2923                return G2AtomRecord(atom, ptrs, self.proj)
2924
2925    def atoms(self):
2926        """Returns a list of atoms present in the phase.
2927
2928        :returns: A list of :class:`G2AtomRecord` objects.
2929
2930        .. seealso::
2931            :meth:`G2Phase.atom`
2932            :class:`G2AtomRecord`
2933        """
2934        ptrs = self.data['General']['AtomPtrs']
2935        output = []
2936        atoms = self.data['Atoms']
2937        for atom in atoms:
2938            output.append(G2AtomRecord(atom, ptrs, self.proj))
2939        return output
2940
2941    def histograms(self):
2942        output = []
2943        for hname in self.data.get('Histograms', {}).keys():
2944            output.append(self.proj.histogram(hname))
2945        return output
2946
2947    @property
2948    def ranId(self):
2949        return self.data['ranId']
2950
2951    @property
2952    def id(self):
2953        return self.data['pId']
2954
2955    @id.setter
2956    def id(self, val):
2957        self.data['pId'] = val
2958
2959    def get_cell(self):
2960        """Returns a dictionary of the cell parameters, with keys:
2961            'length_a', 'length_b', 'length_c', 'angle_alpha', 'angle_beta', 'angle_gamma', 'volume'
2962
2963        :returns: a dict
2964
2965        .. seealso::
2966           :meth:`G2Phase.get_cell_and_esd`
2967
2968        """
2969        cell = self.data['General']['Cell']
2970        return {'length_a': cell[1], 'length_b': cell[2], 'length_c': cell[3],
2971                'angle_alpha': cell[4], 'angle_beta': cell[5], 'angle_gamma': cell[6],
2972                'volume': cell[7]}
2973
2974    def get_cell_and_esd(self):
2975        """
2976        Returns a pair of dictionaries, the first representing the unit cell, the second
2977        representing the estimated standard deviations of the unit cell.
2978
2979        :returns: a tuple of two dictionaries
2980
2981        .. seealso::
2982           :meth:`G2Phase.get_cell`
2983
2984        """
2985        # translated from GSASIIstrIO.ExportBaseclass.GetCell
2986        import GSASIIstrIO as G2stIO
2987        import GSASIIlattice as G2lat
2988        import GSASIImapvars as G2mv
2989        try:
2990            pfx = str(self.id) + '::'
2991            sgdata = self['General']['SGData']
2992            covDict = self.proj['Covariance']['data']
2993
2994            parmDict = dict(zip(covDict.get('varyList',[]),
2995                                covDict.get('variables',[])))
2996            sigDict = dict(zip(covDict.get('varyList',[]),
2997                               covDict.get('sig',[])))
2998
2999            if covDict.get('covMatrix') is not None:
3000                sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
3001                                                  covDict['varyList'],
3002                                                  parmDict))
3003
3004            A, sigA = G2stIO.cellFill(pfx, sgdata, parmDict, sigDict)
3005            cellSig = G2stIO.getCellEsd(pfx, sgdata, A, self.proj['Covariance']['data'])
3006            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
3007            cellDict, cellSigDict = {}, {}
3008            for i, key in enumerate(['length_a', 'length_b', 'length_c',
3009                                     'angle_alpha', 'angle_beta', 'angle_gamma',
3010                                     'volume']):
3011                cellDict[key] = cellList[i]
3012                cellSigDict[key] = cellSig[i]
3013            return cellDict, cellSigDict
3014        except KeyError:
3015            cell = self.get_cell()
3016            return cell, {key: 0.0 for key in cell}
3017
3018    def export_CIF(self, outputname, quickmode=True):
3019        """Write this phase to a .cif file named outputname
3020
3021        :param str outputname: The name of the .cif file to write to
3022        :param bool quickmode: Currently ignored. Carryover from exports.G2export_CIF"""
3023        # This code is all taken from exports/G2export_CIF.py
3024        # Functions copied have the same names
3025        import GSASIImath as G2mth
3026        import GSASIImapvars as G2mv
3027        from exports import G2export_CIF as cif
3028
3029        CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
3030        CIFname = os.path.splitext(self.proj.filename)[0]
3031        CIFname = os.path.split(CIFname)[1]
3032        CIFname = ''.join([c if ord(c) < 128 else ''
3033                           for c in CIFname.replace(' ', '_')])
3034        try:
3035            author = self.proj['Controls']['data'].get('Author','').strip()
3036        except KeyError:
3037            pass
3038        oneblock = True
3039
3040        covDict = self.proj['Covariance']['data']
3041        parmDict = dict(zip(covDict.get('varyList',[]),
3042                            covDict.get('variables',[])))
3043        sigDict = dict(zip(covDict.get('varyList',[]),
3044                           covDict.get('sig',[])))
3045
3046        if covDict.get('covMatrix') is not None:
3047            sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
3048                                              covDict['varyList'],
3049                                              parmDict))
3050
3051        with open(outputname, 'w') as fp:
3052            fp.write(' \n' + 70*'#' + '\n')
3053            cif.WriteCIFitem(fp, 'data_' + CIFname)
3054            # from exports.G2export_CIF.WritePhaseInfo
3055            cif.WriteCIFitem(fp, '\n# phase info for '+str(self.name) + ' follows')
3056            cif.WriteCIFitem(fp, '_pd_phase_name', self.name)
3057            # TODO get esds
3058            cellDict = self.get_cell()
3059            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
3060            names = ['length_a','length_b','length_c',
3061                     'angle_alpha','angle_beta ','angle_gamma',
3062                     'volume']
3063            for key, val in cellDict.items():
3064                cif.WriteCIFitem(fp, '_cell_' + key, G2mth.ValEsd(val))
3065
3066            cif.WriteCIFitem(fp, '_symmetry_cell_setting',
3067                         self.data['General']['SGData']['SGSys'])
3068
3069            spacegroup = self.data['General']['SGData']['SpGrp'].strip()
3070            # regularize capitalization and remove trailing H/R
3071            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
3072            cif.WriteCIFitem(fp, '_symmetry_space_group_name_H-M', spacegroup)
3073
3074            # generate symmetry operations including centering and center of symmetry
3075            SymOpList, offsetList, symOpList, G2oprList, G2opcodes = G2spc.AllOps(
3076                self.data['General']['SGData'])
3077            cif.WriteCIFitem(fp, 'loop_\n    _space_group_symop_id\n    _space_group_symop_operation_xyz')
3078            for i, op in enumerate(SymOpList,start=1):
3079                cif.WriteCIFitem(fp, '   {:3d}  {:}'.format(i,op.lower()))
3080
3081            # TODO skipped histograms, exports/G2export_CIF.py:880
3082
3083            # report atom params
3084            if self.data['General']['Type'] in ['nuclear','macromolecular']:        #this needs macromolecular variant, etc!
3085                cif.WriteAtomsNuclear(fp, self.data, self.name, parmDict, sigDict, [])
3086                # self._WriteAtomsNuclear(fp, parmDict, sigDict)
3087            else:
3088                raise G2ScriptException("no export for "+str(self.data['General']['Type'])+" coordinates implemented")
3089            # report cell contents
3090            cif.WriteComposition(fp, self.data, self.name, parmDict)
3091            if not quickmode and self.data['General']['Type'] == 'nuclear':      # report distances and angles
3092                # WriteDistances(fp,self.name,SymOpList,offsetList,symOpList,G2oprList)
3093                raise NotImplementedError("only quickmode currently supported")
3094            if 'Map' in self.data['General'] and 'minmax' in self.data['General']['Map']:
3095                cif.WriteCIFitem(fp,'\n# Difference density results')
3096                MinMax = self.data['General']['Map']['minmax']
3097                cif.WriteCIFitem(fp,'_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009))
3098                cif.WriteCIFitem(fp,'_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009))
3099
3100
3101    def set_refinements(self, refs):
3102        """Sets the refinement parameter 'key' to the specification 'value'
3103
3104        :param dict refs: A dictionary of the parameters to be set. See
3105                          :ref:`Phase_parameters_table` for a description of
3106                          this dictionary.
3107
3108        :returns: None"""
3109        for key, value in refs.items():
3110            if key == "Cell":
3111                self.data['General']['Cell'][0] = value
3112
3113            elif key == "Atoms":
3114                for atomlabel, atomrefinement in value.items():
3115                    if atomlabel == 'all':
3116                        for atom in self.atoms():
3117                            atom.refinement_flags = atomrefinement
3118                    else:
3119                        atom = self.atom(atomlabel)
3120                        if atom is None:
3121                            raise ValueError("No such atom: " + atomlabel)
3122                        atom.refinement_flags = atomrefinement
3123
3124            elif key == "LeBail":
3125                hists = self.data['Histograms']
3126                for hname, hoptions in hists.items():
3127                    if 'LeBail' not in hoptions:
3128                        hoptions['newLeBail'] = bool(True)
3129                    hoptions['LeBail'] = bool(value)
3130            else:
3131                raise ValueError("Unknown key:", key)
3132
3133    def clear_refinements(self, refs):
3134        """Clears a given set of parameters.
3135
3136        :param dict refs: The parameters to clear"""
3137        for key, value in refs.items():
3138            if key == "Cell":
3139                self.data['General']['Cell'][0] = False
3140            elif key == "Atoms":
3141                cx, ct, cs, cia = self.data['General']['AtomPtrs']
3142
3143                for atomlabel in value:
3144                    atom = self.atom(atomlabel)
3145                    # Set refinement to none
3146                    atom.refinement_flags = ' '
3147            elif key == "LeBail":
3148                hists = self.data['Histograms']
3149                for hname, hoptions in hists.items():
3150                    if 'LeBail' not in hoptions:
3151                        hoptions['newLeBail'] = True
3152                    hoptions['LeBail'] = False
3153            else:
3154                raise ValueError("Unknown key:", key)
3155
3156    def set_HAP_refinements(self, refs, histograms='all'):
3157        """Sets the given HAP refinement parameters between this phase and
3158        the given histograms
3159
3160        :param dict refs: A dictionary of the parameters to be set. See
3161                          :ref:`HAP_parameters_table` for a description of this
3162                          dictionary.
3163        :param histograms: Either 'all' (default) or a list of the histograms
3164            whose HAP parameters will be set with this phase. Histogram and phase
3165            must already be associated
3166
3167        :returns: None
3168        """
3169        if not self.data['Histograms']:
3170            print("Error likely: Phase {} has no linked histograms".format(self.name))
3171            return
3172           
3173        if histograms == 'all':
3174            histograms = self.data['Histograms'].values()
3175        else:
3176            histograms = [self.data['Histograms'][h.name] for h in histograms
3177                          if h.name in self.data['Histograms']]
3178
3179        if not histograms:
3180            print("Skipping HAP set for phase {}, no selected histograms".format(self.name))
3181            return
3182        for key, val in refs.items():
3183            for h in histograms:
3184                if key == 'Babinet':
3185                    try:
3186                        sets = list(val)
3187                    except ValueError:
3188                        sets = ['BabA', 'BabU']
3189                    for param in sets:
3190                        if param not in ['BabA', 'BabU']:
3191                            raise ValueError("Not sure what to do with" + param)
3192                        for hist in histograms:
3193                            hist['Babinet'][param][1] = True
3194                elif key == 'Extinction':
3195                    for h in histograms:
3196                        h['Extinction'][1] = bool(val)
3197                elif key == 'HStrain':
3198                    for h in histograms:
3199                        h['HStrain'][1] = [bool(val) for p in h['HStrain'][1]]
3200                elif key == 'Mustrain':
3201                    for h in histograms:
3202                        mustrain = h['Mustrain']
3203                        newType = None
3204                        direction = None
3205                        if isinstance(val, strtypes):
3206                            if val in ['isotropic', 'uniaxial', 'generalized']:
3207                                newType = val
3208                            else:
3209                                raise ValueError("Not a Mustrain type: " + val)
3210                        elif isinstance(val, dict):
3211                            newType = val.get('type', None)
3212                            direction = val.get('direction', None)
3213
3214                        if newType:
3215                            mustrain[0] = newType
3216                            if newType == 'isotropic':
3217                                mustrain[2][0] = True == val.get('refine',False)
3218                                mustrain[5] = [False for p in mustrain[4]]
3219                            elif newType == 'uniaxial':
3220                                if 'refine' in val:
3221                                    mustrain[2][0] = False
3222                                    types = val['refine']
3223                                    if isinstance(types, strtypes):
3224                                        types = [types]
3225                                    elif isinstance(types, bool):
3226                                        mustrain[2][1] = types
3227                                        mustrain[2][2] = types
3228                                        types = []
3229                                    else:
3230                                        raise ValueError("Not sure what to do with: "
3231                                                         + str(types))
3232                                else:
3233                                    types = []
3234
3235                                for unitype in types:
3236                                    if unitype == 'equatorial':
3237                                        mustrain[2][0] = True
3238                                    elif unitype == 'axial':
3239                                        mustrain[2][1] = True
3240                                    else:
3241                                        msg = 'Invalid uniaxial mustrain type'
3242                                        raise ValueError(msg + ': ' + unitype)
3243                            else:  # newtype == 'generalized'
3244                                mustrain[2] = [False for p in mustrain[1]]
3245                                if 'refine' in val:
3246                                    mustrain[5] = [True == val['refine']]*len(mustrain[5])
3247
3248                        if direction:
3249                            if len(direction) != 3:
3250                                raise ValueError("Expected hkl, found", direction)
3251                            direction = [int(n) for n in direction]
3252                            mustrain[3] = direction
3253                elif key == 'Size':
3254                    newSize = None
3255                    if 'value' in val:
3256                        newSize = float(val['value'])
3257                    for h in histograms:
3258                        size = h['Size']
3259                        newType = None
3260                        direction = None
3261                        if isinstance(val, strtypes):
3262                            if val in ['isotropic', 'uniaxial', 'ellipsoidal']:
3263                                newType = val
3264                            else:
3265                                raise ValueError("Not a valid Size type: " + val)
3266                        elif isinstance(val, dict):
3267                            newType = val.get('type', None)
3268                            direction = val.get('direction', None)
3269
3270                        if newType:
3271                            size[0] = newType
3272                            refine = True == val.get('refine')
3273                            if newType == 'isotropic' and refine is not None:
3274                                size[2][0] = bool(refine)
3275                                if newSize: size[1][0] = newSize
3276                            elif newType == 'uniaxial' and refine is not None:
3277                                size[2][1] = bool(refine)
3278                                size[2][2] = bool(refine)
3279                                if newSize: size[1][1] = size[1][2] =newSize
3280                            elif newType == 'ellipsoidal' and refine is not None:
3281                                size[5] = [bool(refine) for p in size[5]]
3282                                if newSize: size[4] = [newSize for p in size[4]]
3283
3284                        if direction:
3285                            if len(direction) != 3:
3286                                raise ValueError("Expected hkl, found", direction)
3287                            direction = [int(n) for n in direction]
3288                            size[3] = direction
3289                elif key == 'Pref.Ori.':
3290                    for h in histograms:
3291                        h['Pref.Ori.'][2] = bool(val)
3292                elif key == 'Show':
3293                    for h in histograms:
3294                        h['Show'] = bool(val)
3295                elif key == 'Use':
3296                    for h in histograms:
3297                        h['Use'] = bool(val)
3298                elif key == 'Scale':
3299                    for h in histograms:
3300                        h['Scale'][1] = bool(val)
3301                else:
3302                    print(u'Unknown HAP key: '+key)
3303
3304    def getHAPvalues(self, histname):
3305        """Returns a dict with HAP values for the selected histogram
3306
3307        :param histogram: is a histogram object (:class:`G2PwdrData`) or
3308            a histogram name or the index number of the histogram
3309
3310        :returns: HAP value dict
3311        """
3312        if isinstance(histname, G2PwdrData):
3313            histname = histname.name
3314        elif histname in self.data['Histograms']:
3315            pass
3316        elif type(histname) is int:
3317            histname = self.proj.histograms()[histname].name
3318        else:
3319            raise G2ScriptException("Invalid histogram reference: "+str(histname))
3320        return self.data['Histograms'][histname]
3321                   
3322    def clear_HAP_refinements(self, refs, histograms='all'):
3323        """Clears the given HAP refinement parameters between this phase and
3324        the given histograms
3325
3326        :param dict refs: A dictionary of the parameters to be cleared.
3327        :param histograms: Either 'all' (default) or a list of the histograms
3328            whose HAP parameters will be cleared with this phase. Histogram and
3329            phase must already be associated
3330
3331        :returns: None
3332        """
3333        if histograms == 'all':
3334            histograms = self.data['Histograms'].values()
3335        else:
3336            histograms = [self.data['Histograms'][h.name] for h in histograms
3337                          if h.name in self.data['Histograms']]
3338
3339        for key, val in refs.items():
3340            for h in histograms:
3341                if key == 'Babinet':
3342                    try:
3343                        sets = list(val)
3344                    except ValueError:
3345                        sets = ['BabA', 'BabU']
3346                    for param in sets:
3347                        if param not in ['BabA', 'BabU']:
3348                            raise ValueError("Not sure what to do with" + param)
3349                        for hist in histograms:
3350                            hist['Babinet'][param][1] = False
3351                elif key == 'Extinction':
3352                    for h in histograms:
3353                        h['Extinction'][1] = False
3354                elif key == 'HStrain':
3355                    for h in histograms:
3356                        h['HStrain'][1] = [False for p in h['HStrain'][1]]
3357                elif key == 'Mustrain':
3358                    for h in histograms:
3359                        mustrain = h['Mustrain']
3360                        mustrain[2] = [False for p in mustrain[2]]
3361                        mustrain[5] = [False for p in mustrain[4]]
3362                elif key == 'Pref.Ori.':
3363                    for h in histograms:
3364                        h['Pref.Ori.'][2] = False
3365                elif key == 'Show':
3366                    for h in histograms:
3367                        h['Show'] = False
3368                elif key == 'Size':
3369                    for h in histograms:
3370                        size = h['Size']
3371                        size[2] = [False for p in size[2]]
3372                        size[5] = [False for p in size[5]]
3373                elif key == 'Use':
3374                    for h in histograms:
3375                        h['Use'] = False
3376                elif key == 'Scale':
3377                    for h in histograms:
3378                        h['Scale'][1] = False
3379                else:
3380                    print(u'Unknown HAP key: '+key)
3381
3382class G2Image(G2ObjectWrapper):
3383    """Wrapper for an IMG tree entry, containing an image and various metadata.
3384
3385    Example:
3386
3387    >>> gpx = G2sc.G2Project(filename='itest.gpx')
3388    >>> img2 = gpx.add_image(idata,fmthint="TIF")
3389    >>> img2[0].loadControls('stdSettings.imctrl')
3390    >>> img2[0].setCalibrant('Si    SRM640c')
3391    >>> img2[0].loadMasks('stdMasks.immask')
3392    >>> img2[0].Recalibrate()
3393    >>> img2[0].setControl('outAzimuths',3)
3394    >>> pwdrList = img2[0].Integrate()
3395
3396    """
3397    # parameters in that can be accessed via setControl
3398    ControlList = {
3399        'int': ['calibskip', 'pixLimit', 'edgemin', 'outChannels',
3400                    'outAzimuths'],
3401        'float': ['cutoff', 'setdist', 'wavelength', 'Flat Bkg',
3402                      'azmthOff', 'tilt', 'calibdmin', 'rotation',
3403                      'distance', 'DetDepth'],
3404        'bool': ['setRings', 'setDefault', 'centerAzm', 'fullIntegrate',
3405                     'DetDepthRef', 'showLines'],
3406        'str': ['SampleShape', 'binType', 'formatName', 'color',
3407                    'type', ],
3408        'list': ['GonioAngles', 'IOtth', 'LRazimuth', 'Oblique', 'PolaVal',
3409                   'SampleAbs', 'center', 'ellipses', 'linescan',
3410                    'pixelSize', 'range', 'ring', 'rings', 'size', ],
3411        'dict': ['varyList'],
3412        }
3413    '''Defines the items known to exist in the Image Controls tree section
3414    and the item's data types. A few are not included
3415    ('background image', 'dark image', 'Gain map', and 'calibrant') because
3416    these items have special set routines,
3417    where references to entries are checked to make sure their values are
3418    correct.
3419    ''' 
3420    # this may need future attention
3421       
3422    def __init__(self, data, name, proj):
3423        self.data = data
3424        self.name = name
3425        self.proj = proj
3426
3427    def setControl(self,arg,value):
3428        '''Set an Image Controls parameter in the current image.
3429        If the parameter is not found an exception is raised.
3430
3431        :param str arg: the name of a parameter (dict entry) in the
3432          image. The parameter must be found in :data:`ControlList`
3433          or an exception is raised.
3434        :param value: the value to set the parameter. The value is
3435          cast as the appropriate type from :data:`ControlList`.
3436        '''
3437        for typ in self.ControlList:
3438            if arg in self.ControlList[typ]: break
3439        else:
3440            print('Allowed args:\n',[nam for nam,typ in self.findControl('')])
3441            raise Exception('arg {} not defined in G2Image.setControl'
3442                                .format(arg))
3443        try:
3444            if typ == 'int':
3445                self.data['Image Controls'][arg] = int(value)
3446            elif typ == 'float':
3447                self.data['Image Controls'][arg] = float(value)
3448            elif typ == 'bool':
3449                self.data['Image Controls'][arg] = bool(value)
3450            elif typ == 'str':
3451                self.data['Image Controls'][arg] = str(value)
3452            elif typ == 'list':
3453                self.data['Image Controls'][arg] = list(value)
3454            elif typ == 'dict':
3455                self.data['Image Controls'][arg] = dict(value)
3456            else:
3457                raise Exception('Unknown type {} for arg {} in  G2Image.setControl'
3458                                    .format(typ,arg))
3459        except:
3460            raise Exception('Error formatting value {} as type {} for arg {} in  G2Image.setControl'
3461                                    .format(value,typ,arg))
3462
3463    def getControl(self,arg):
3464        '''Return an Image Controls parameter in the current image.
3465        If the parameter is not found an exception is raised.
3466
3467        :param str arg: the name of a parameter (dict entry) in the
3468          image.
3469        :returns: the value as a int, float, list,...
3470        '''
3471        if arg in self.data['Image Controls']:
3472            return self.data['Image Controls'][arg]
3473        print(self.findControl(''))
3474        raise Exception('arg {} not defined in G2Image.getControl'.format(arg))
3475
3476    def findControl(self,arg=''):
3477        '''Finds the Image Controls parameter(s) in the current image
3478        that match the string in arg. Default is '' which returns all
3479        parameters.
3480
3481            Example:
3482
3483            >>> findControl('calib')
3484            [['calibskip', 'int'], ['calibdmin', 'float'], ['calibrant', 'str']]
3485
3486        :param str arg: a string containing part of the name of a
3487          parameter (dict entry) in the image's Image Controls.
3488        :returns: a list of matching entries in form
3489          [['item','type'], ['item','type'],...] where each 'item' string
3490          contains the sting in arg.
3491        '''
3492        matchList = []
3493        for typ in self.ControlList:
3494            for item in self.ControlList[typ]:
3495                if arg in item:
3496                    matchList.append([item,typ])
3497        return matchList
3498
3499    def setCalibrant(self,calib):
3500        '''Set a calibrant for the current image
3501
3502        :param str calib: specifies a calibrant name which must be one of
3503          the entries in file ImageCalibrants.py. This is validated and
3504          an error provides a list of valid choices.
3505        '''
3506        import ImageCalibrants as calFile
3507        if calib in calFile.Calibrants.keys():
3508            self.data['Image Controls']['calibrant'] = calib
3509            return
3510        print('Calibrant {} is not valid. Valid calibrants'.format(calib))
3511        for i in calFile.Calibrants.keys():
3512            if i: print('\t"{}"'.format(i))
3513       
3514    def setControlFile(self,typ,imageRef,mult=None):
3515        '''Set a image to be used as a background/dark/gain map image
3516
3517        :param str typ: specifies image type, which must be one of:
3518           'background image', 'dark image', 'gain map'; N.B. only the first
3519           four characters must be specified and case is ignored.
3520        :param imageRef: A reference to the desired image. Either the Image
3521          tree name (str), the image's index (int) or
3522          a image object (:class:`G2Image`)
3523        :param float mult: a multiplier to be applied to the image (not used
3524          for 'Gain map'; required for 'background image', 'dark image'
3525        '''
3526        if 'back' in typ.lower():
3527            key = 'background image'
3528            mult = float(mult)
3529        elif 'dark' in typ.lower():
3530            key = 'dark image'
3531            mult = float(mult)
3532        elif 'gain' in typ.lower():
3533            #key = 'Gain map'
3534            if mult is not None:
3535                print('Ignoring multiplier for Gain map')
3536            mult = None
3537        else:
3538            raise Exception("Invalid typ {} for setControlFile".format(typ))
3539        imgNam = self.proj.image(imageRef).name
3540        if mult is None:
3541            self.data['Image Controls']['Gain map'] = imgNam
3542        else:
3543            self.data['Image Controls'][key] = [imgNam,mult]
3544
3545    def loadControls(self,filename):
3546        '''load controls from a .imctrl file
3547
3548        :param str filename: specifies a file to be read, which should end
3549          with .imctrl
3550        '''
3551        File = open(filename,'r')
3552        Slines = File.readlines()
3553        File.close()
3554        G2fil.LoadControls(Slines,self.data['Image Controls'])
3555        print('file {} read into {}'.format(filename,self.name))
3556
3557    def saveControls(self,filename):
3558        '''write current controls values to a .imctrl file
3559
3560        :param str filename: specifies a file to write, which should end
3561          with .imctrl
3562        '''
3563        G2fil.WriteControls(filename,self.data['Image Controls'])
3564        print('file {} written from {}'.format(filename,self.name))
3565
3566    def loadMasks(self,filename,ignoreThreshold=False):
3567        '''load masks from a .immask file
3568
3569        :param str filename: specifies a file to be read, which should end
3570          with .immask
3571        :param bool ignoreThreshold: If True, masks are loaded with
3572          threshold masks. Default is False which means any Thresholds
3573          in the file are ignored.
3574        '''
3575        G2fil.readMasks(filename,self.data['Masks'],ignoreThreshold)
3576        print('file {} read into {}'.format(filename,self.name))
3577       
3578    def getVary(self,*args):
3579        '''Return the refinement flag(s) for Image Controls parameter(s)
3580        in the current image.
3581        If the parameter is not found, an exception is raised.
3582
3583        :param str arg: the name of a refinement parameter in the
3584          varyList for the image. The name should be one of
3585          'dep', 'det-X', 'det-Y', 'dist', 'phi', 'tilt', or 'wave'
3586        :param str arg1: the name of a parameter (dict entry) as before,
3587          optional
3588
3589
3590        :returns: a list of bool value(s)
3591        '''
3592        res = []
3593        for arg in args:
3594            if arg in self.data['Image Controls']['varyList']:
3595                res.append(self.data['Image Controls']['varyList'][arg])
3596            else:
3597                raise Exception('arg {} not defined in G2Image.getVary'.format(arg))
3598        return res
3599   
3600    def setVary(self,arg,value):
3601        '''Set a refinement flag for Image Controls parameter in the
3602        current image.
3603        If the parameter is not found an exception is raised.
3604
3605        :param str arg: the name of a refinement parameter in the
3606          varyList for the image. The name should be one of
3607          'dep', 'det-X', 'det-Y', 'dist', 'phi', 'tilt', or 'wave'
3608        :param str arg: the name of a parameter (dict entry) in the
3609          image. The parameter must be found in :data:`ControlList`
3610          or an exception is raised.
3611        :param value: the value to set the parameter. The value is
3612          cast as the appropriate type from :data:`ControlList`.
3613        '''
3614        if arg in self.data['Image Controls']['varyList']:
3615            self.data['Image Controls']['varyList'][arg] = bool(value)
3616        else:
3617            raise Exception('arg {} not defined in G2Image.setVary'.format(arg))
3618
3619    def Recalibrate(self):
3620        '''Invokes a recalibration fit (same as Image Controls/Calibration/Recalibrate
3621        menu command). Note that for this to work properly, the calibration
3622        coefficients (center, wavelength, ditance & tilts) must be fairly close.
3623        This may produce a better result if run more than once.
3624        '''
3625        LoadG2fil()
3626        ImageZ = GetCorrImage(Readers['Image'],self.proj,self)
3627        G2img.ImageRecalibrate(None,ImageZ,self.data['Image Controls'],self.data['Masks'])
3628
3629    def Integrate(self,name=None):
3630        '''Invokes an image integration (same as Image Controls/Integration/Integrate
3631        menu command). All parameters will have previously been set with Image Controls
3632        so no input is needed here. Note that if integration is performed on an
3633        image more than once, histogram entries may be overwritten. Use the name
3634        parameter to prevent this if desired.
3635
3636        :param str name: base name for created histogram(s). If None (default),
3637          the histogram name is taken from the image name.
3638        :returns: a list of created histogram (:class:`G2PwdrData`) objects.
3639        '''
3640        ImageZ = GetCorrImage(Readers['Image'],self.proj,self)
3641        # do integration
3642        ints,azms,Xvals,cancel = G2img.ImageIntegrate(ImageZ,self.data['Image Controls'],self.data['Masks'])
3643        # code from here on based on G2IO.SaveIntegration, but places results in the current
3644        # project rather than tree
3645        X = Xvals[:-1]
3646        N = len(X)
3647
3648        data = self.data['Image Controls']
3649        Comments = self.data['Comments']
3650        # make name from image, unless overridden
3651        if name:
3652            if not name.startswith(data['type']+' '):
3653                name = data['type']+' '+name
3654        else:
3655            name = self.name.replace('IMG ',data['type']+' ')
3656        if 'PWDR' in name:
3657            if 'target' in data:
3658                names = ['Type','Lam1','Lam2','I(L2)/I(L1)','Zero','Polariz.','U','V','W','X','Y','Z','SH/L','Azimuth'] 
3659                codes = [0 for i in range(14)]
3660            else:
3661                names = ['Type','Lam','Zero','Polariz.','U','V','W','X','Y','Z','SH/L','Azimuth'] 
3662                codes = [0 for i in range(12)]
3663        elif 'SASD' in name:
3664            names = ['Type','Lam','Zero','Azimuth'] 
3665            codes = [0 for i in range(4)]
3666            X = 4.*np.pi*npsind(X/2.)/data['wavelength']    #convert to q
3667        Xminmax = [X[0],X[-1]]
3668        Azms = []
3669        dazm = 0.
3670        if data['fullIntegrate'] and data['outAzimuths'] == 1:
3671            Azms = [45.0,]                              #a poor man's average?
3672        else:
3673            for i,azm in enumerate(azms[:-1]):
3674                if azm > 360. and azms[i+1] > 360.:
3675                    Azms.append(G2img.meanAzm(azm%360.,azms[i+1]%360.))
3676                else:   
3677                    Azms.append(G2img.meanAzm(azm,azms[i+1]))
3678            dazm = np.min(np.abs(np.diff(azms)))/2.
3679        # pull out integration results and make histograms for each
3680        IntgOutList = []
3681        for i,azm in enumerate(azms[:-1]):
3682            Aname = name+" Azm= %.2f"%((azm+dazm)%360.)
3683            # MT dict to contain histogram
3684            HistDict = {}
3685            histItems = [Aname]
3686            Sample = G2obj.SetDefaultSample()       #set as Debye-Scherrer
3687            Sample['Gonio. radius'] = data['distance']
3688            Sample['Omega'] = data['GonioAngles'][0]
3689            Sample['Chi'] = data['GonioAngles'][1]
3690            Sample['Phi'] = data['GonioAngles'][2]
3691            Sample['Azimuth'] = (azm+dazm)%360.    #put here as bin center
3692            polariz = 0.99    #set default polarization for synchrotron radiation!
3693            for item in Comments:
3694                if 'polariz' in item:
3695                    try:
3696                        polariz = float(item.split('=')[1])
3697                    except:
3698                        polariz = 0.99
3699                for key in ('Temperature','Pressure','Time','FreePrm1','FreePrm2','FreePrm3','Omega',
3700                    'Chi','Phi'):
3701                    if key.lower() in item.lower():
3702                        try:
3703                            Sample[key] = float(item.split('=')[1])
3704                        except:
3705                            pass
3706                if 'label_prm' in item.lower():
3707                    for num in ('1','2','3'):
3708                        if 'label_prm'+num in item.lower():
3709                            Controls['FreePrm'+num] = item.split('=')[1].strip()
3710            if 'PWDR' in Aname:
3711                if 'target' in data:    #from lab x-ray 2D imaging data
3712                    wave1,wave2 = waves[data['target']]
3713                    parms = ['PXC',wave1,wave2,0.5,0.0,polariz,290.,-40.,30.,6.,-14.,0.0,0.0001,Azms[i]]
3714                else:
3715                    parms = ['PXC',data['wavelength'],0.0,polariz,1.0,-0.10,0.4,0.30,1.0,0.0,0.0001,Azms[i]]
3716            elif 'SASD' in Aname:
3717                Sample['Trans'] = data['SampleAbs'][0]
3718                parms = ['LXC',data['wavelength'],0.0,Azms[i]]
3719            Y = ints[i]
3720            Ymin = np.min(Y)
3721            Ymax = np.max(Y)
3722            W = np.where(Y>0.,1./Y,1.e-6)                    #probably not true
3723            section = 'Comments'
3724            histItems += [section]
3725            HistDict[section] = Comments
3726            section = 'Limits'
3727            histItems += [section]
3728            HistDict[section] = copy.deepcopy([tuple(Xminmax),Xminmax])
3729            if 'PWDR' in Aname:
3730                section = 'Background'
3731                histItems += [section]
3732                HistDict[section] = [['chebyschev',1,3,1.0,0.0,0.0],
3733                    {'nDebye':0,'debyeTerms':[],'nPeaks':0,'peaksList':[]}]
3734            inst = [dict(zip(names,zip(parms,parms,codes))),{}]
3735            for item in inst[0]:
3736                inst[0][item] = list(inst[0][item])
3737            section = 'Instrument Parameters'
3738            histItems += [section]
3739            HistDict[section] = inst
3740            if 'PWDR' in Aname:
3741                section = 'Sample Parameters'
3742                histItems += [section]
3743                HistDict[section] = Sample
3744                section = 'Peak List'
3745                histItems += [section]
3746                HistDict[section] = {'sigDict':{},'peaks':[]}
3747                section = 'Index Peak List'
3748                histItems += [section]
3749                HistDict[section] = [[],[]]
3750                section = 'Unit Cells List'
3751                histItems += [section]
3752                HistDict[section] = []
3753                section = 'Reflection Lists'
3754                histItems += [section]
3755                HistDict[section] = {}
3756            elif 'SASD' in Aname:             
3757                section = 'Substances'
3758                histItems += [section]
3759                HistDict[section] = G2pdG.SetDefaultSubstances()  # this needs to be moved
3760                section = 'Sample Parameters'
3761                histItems += [section]
3762                HistDict[section] = Sample
3763                section = 'Models'
3764                histItems += [section]
3765                HistDict[section] = G2pdG.SetDefaultSASDModel() # this needs to be moved
3766            valuesdict = {
3767                'wtFactor':1.0,'Dummy':False,'ranId':ran.randint(0,sys.maxsize),'Offset':[0.0,0.0],'delOffset':0.02*Ymax,
3768                'refOffset':-0.1*Ymax,'refDelt':0.1*Ymax,'Yminmax':[Ymin,Ymax]}
3769            # if Aname is already in the project replace it
3770            for j in self.proj.names:
3771                if j[0] == Aname: 
3772                    print('Replacing "{}" in project'.format(Aname))
3773                    break
3774            else:
3775                print('Adding "{}" to project'.format(Aname))
3776                self.proj.names.append([Aname]+
3777                        [u'Comments',u'Limits',u'Background',u'Instrument Parameters',
3778                         u'Sample Parameters', u'Peak List', u'Index Peak List',
3779                         u'Unit Cells List', u'Reflection Lists'])
3780            HistDict['data'] = [valuesdict,
3781                    [np.array(X),np.array(Y),np.array(W),np.zeros(N),np.zeros(N),np.zeros(N)]]
3782            self.proj.data[Aname] = HistDict
3783            IntgOutList.append(self.proj.histogram(Aname))
3784        return IntgOutList
3785
3786##########################
3787# Command Line Interface #
3788##########################
3789# Each of these takes an argparse.Namespace object as their argument,
3790# representing the parsed command-line arguments for the relevant subcommand.
3791# The argument specification for each is in the subcommands dictionary (see
3792# below)
3793
3794commandhelp={}
3795commandhelp["create"] = "creates a GSAS-II project, optionally adding histograms and/or phases"
3796def create(args):
3797    """Implements the create command-line subcommand. This creates a GSAS-II project, optionally adding histograms and/or phases::
3798
3799  usage: GSASIIscriptable.py create [-h] [-d HISTOGRAMS [HISTOGRAMS ...]]
3800                                  [-i IPARAMS [IPARAMS ...]]
3801                                  [-p PHASES [PHASES ...]]
3802                                  filename
3803                                 
3804positional arguments::
3805
3806  filename              the project file to create. should end in .gpx
3807
3808optional arguments::
3809
3810  -h, --help            show this help message and exit
3811  -d HISTOGRAMS [HISTOGRAMS ...], --histograms HISTOGRAMS [HISTOGRAMS ...]
3812                        list of datafiles to add as histograms
3813  -i IPARAMS [IPARAMS ...], --iparams IPARAMS [IPARAMS ...]
3814                        instrument parameter file, must be one for every
3815                        histogram
3816  -p PHASES [PHASES ...], --phases PHASES [PHASES ...]
3817                        list of phases to add. phases are automatically
3818                        associated with all histograms given.
3819
3820    """
3821    proj = G2Project(gpxname=args.filename)
3822
3823    hist_objs = []
3824    if args.histograms:
3825        for h,i in zip(args.histograms,args.iparams):
3826            print("Adding histogram from",h,"with instparm ",i)
3827            hist_objs.append(proj.add_powder_histogram(h, i))
3828
3829    if args.phases: 
3830        for p in args.phases:
3831            print("Adding phase from",p)
3832            proj.add_phase(p, histograms=hist_objs)
3833        print('Linking phase(s) to histogram(s):')
3834        for h in hist_objs:
3835            print ('   '+h.name)
3836
3837    proj.save()
3838
3839commandhelp["add"] = "adds histograms and/or phases to GSAS-II project"
3840def add(args):
3841    """Implements the add command-line subcommand. This adds histograms and/or phases to GSAS-II project::
3842
3843  usage: GSASIIscriptable.py add [-h] [-d HISTOGRAMS [HISTOGRAMS ...]]
3844                               [-i IPARAMS [IPARAMS ...]]
3845                               [-hf HISTOGRAMFORMAT] [-p PHASES [PHASES ...]]
3846                               [-pf PHASEFORMAT] [-l HISTLIST [HISTLIST ...]]
3847                               filename
3848
3849
3850positional arguments::
3851
3852  filename              the project file to open. Should end in .gpx
3853
3854optional arguments::
3855
3856  -h, --help            show this help message and exit
3857  -d HISTOGRAMS [HISTOGRAMS ...], --histograms HISTOGRAMS [HISTOGRAMS ...]
3858                        list of datafiles to add as histograms
3859  -i IPARAMS [IPARAMS ...], --iparams IPARAMS [IPARAMS ...]
3860                        instrument parameter file, must be one for every
3861                        histogram
3862  -hf HISTOGRAMFORMAT, --histogramformat HISTOGRAMFORMAT
3863                        format hint for histogram import. Applies to all
3864                        histograms
3865  -p PHASES [PHASES ...], --phases PHASES [PHASES ...]
3866                        list of phases to add. phases are automatically
3867                        associated with all histograms given.
3868  -pf PHASEFORMAT, --phaseformat PHASEFORMAT
3869                        format hint for phase import. Applies to all phases.
3870                        Example: -pf CIF
3871  -l HISTLIST [HISTLIST ...], --histlist HISTLIST [HISTLIST ...]
3872                        list of histgram indices to associate with added
3873                        phases. If not specified, phases are associated with
3874                        all previously loaded histograms. Example: -l 2 3 4
3875   
3876    """
3877    proj = G2Project(args.filename)
3878
3879    if args.histograms:
3880        for h,i in zip(args.histograms,args.iparams):
3881            print("Adding histogram from",h,"with instparm ",i)
3882            proj.add_powder_histogram(h, i, fmthint=args.histogramformat)
3883
3884    if args.phases: 
3885        if not args.histlist:
3886            histlist = proj.histograms()
3887        else:
3888            histlist = [proj.histogram(i) for i in args.histlist]
3889
3890        for p in args.phases:
3891            print("Adding phase from",p)
3892            proj.add_phase(p, histograms=histlist, fmthint=args.phaseformat)
3893           
3894        if not args.histlist:
3895            print('Linking phase(s) to all histogram(s)')
3896        else:
3897            print('Linking phase(s) to histogram(s):')
3898            for h in histlist:
3899                print ('   '+h.name)
3900
3901    proj.save()
3902
3903
3904commandhelp["dump"] = "Shows the contents of a GSAS-II project"
3905def dump(args):
3906    """Implements the dump command-line subcommand, which shows the contents of a GSAS-II project::
3907
3908       usage: GSASIIscriptable.py dump [-h] [-d] [-p] [-r] files [files ...]
3909
3910positional arguments::
3911
3912  files
3913
3914optional arguments::
3915
3916  -h, --help        show this help message and exit
3917  -d, --histograms  list histograms in files, overrides --raw
3918  -p, --phases      list phases in files, overrides --raw
3919  -r, --raw         dump raw file contents, default
3920 
3921    """
3922    if not args.histograms and not args.phases:
3923        args.raw = True
3924    if args.raw:
3925        import IPython.lib.pretty as pretty
3926
3927    for fname in args.files:
3928        if args.raw:
3929            proj, nameList = LoadDictFromProjFile(fname)
3930            print("file:", fname)
3931            print("names:", nameList)
3932            for key, val in proj.items():
3933                print(key, ":")
3934                pretty.pprint(val)
3935        else:
3936            proj = G2Project(fname)
3937            if args.histograms:
3938                hists = proj.histograms()
3939                for h in hists:
3940                    print(fname, "hist", h.id, h.name)
3941            if args.phases:
3942                phase_list = proj.phases()
3943                for p in phase_list:
3944                    print(fname, "phase", p.id, p.name)
3945
3946
3947commandhelp["browse"] = "Load a GSAS-II project and then open a IPython shell to browse it"
3948def IPyBrowse(args):
3949    """Load a .gpx file and then open a IPython shell to browse it::
3950
3951  usage: GSASIIscriptable.py browse [-h] files [files ...]
3952
3953positional arguments::
3954
3955  files       list of files to browse
3956
3957optional arguments::
3958
3959  -h, --help  show this help message and exit
3960
3961    """
3962    for fname in args.files:
3963        proj, nameList = LoadDictFromProjFile(fname)
3964        msg = "\nfname {} loaded into proj (dict) with names in nameList".format(fname)
3965        GSASIIpath.IPyBreak_base(msg)
3966        break
3967
3968
3969commandhelp["refine"] = '''
3970Conducts refinements on GSAS-II projects according to a list of refinement
3971steps in a JSON dict
3972'''
3973def refine(args):
3974    """Implements the refine command-line subcommand:
3975    conducts refinements on GSAS-II projects according to a JSON refinement dict::
3976
3977        usage: GSASIIscriptable.py refine [-h] gpxfile [refinements]
3978
3979positional arguments::
3980
3981  gpxfile      the project file to refine
3982  refinements  json file of refinements to apply. if not present refines file
3983               as-is
3984
3985optional arguments::
3986
3987  -h, --help   show this help message and exit
3988 
3989    """
3990    proj = G2Project(args.gpxfile)
3991    if args.refinements is None:
3992        proj.refine()
3993    else:
3994        import json
3995        with open(args.refinements) as refs:
3996            refs = json.load(refs)
3997        if type(refs) is not dict:
3998            raise G2ScriptException("Error: JSON object must be a dict.")
3999        if "code" in refs:
4000            print("executing code:\n|  ",'\n|  '.join(refs['code']))
4001            exec('\n'.join(refs['code']))
4002        proj.do_refinements(refs['refinements'])
4003
4004
4005commandhelp["seqrefine"] = "Not implemented. Placeholder for eventual sequential refinement implementation"
4006def seqrefine(args):
4007    """Future implementation for the seqrefine command-line subcommand """
4008    raise NotImplementedError("seqrefine is not yet implemented")
4009
4010
4011commandhelp["export"] = "Export phase as CIF"
4012def export(args):
4013    """Implements the export command-line subcommand: Exports phase as CIF::
4014
4015      usage: GSASIIscriptable.py export [-h] gpxfile phase exportfile
4016
4017positional arguments::
4018
4019  gpxfile     the project file from which to export
4020  phase       identifier of phase to export
4021  exportfile  the .cif file to export to
4022
4023optional arguments::
4024
4025  -h, --help  show this help message and exit
4026
4027    """
4028    proj = G2Project(args.gpxfile)
4029    phase = proj.phase(args.phase)
4030    phase.export_CIF(args.exportfile)
4031
4032
4033def _args_kwargs(*args, **kwargs):
4034    return args, kwargs
4035
4036# A dictionary of the name of each subcommand, and a tuple
4037# of its associated function and a list of its arguments
4038# The arguments are passed directly to the add_argument() method
4039# of an argparse.ArgumentParser
4040
4041subcommands = {"create":
4042               (create, [_args_kwargs('filename',
4043                                      help='the project file to create. should end in .gpx'),
4044
4045                         _args_kwargs('-d', '--histograms',
4046                                      nargs='+',
4047                                      help='list of datafiles to add as histograms'),
4048                                     
4049                         _args_kwargs('-i', '--iparams',
4050                                      nargs='+',
4051                                      help='instrument parameter file, must be one'
4052                                           ' for every histogram'
4053                                      ),
4054
4055                         _args_kwargs('-p', '--phases',
4056                                      nargs='+',
4057                                      help='list of phases to add. phases are '
4058                                           'automatically associated with all '
4059                                           'histograms given.')]),
4060               "add": (add, [_args_kwargs('filename',
4061                                      help='the project file to open. Should end in .gpx'),
4062
4063                         _args_kwargs('-d', '--histograms',
4064                                      nargs='+',
4065                                      help='list of datafiles to add as histograms'),
4066                                     
4067                         _args_kwargs('-i', '--iparams',
4068                                      nargs='+',
4069                                      help='instrument parameter file, must be one'
4070                                           ' for every histogram'
4071                                      ),
4072                                     
4073                         _args_kwargs('-hf', '--histogramformat',
4074                                      help='format hint for histogram import. Applies to all'
4075                                           ' histograms'
4076                                      ),
4077
4078                         _args_kwargs('-p', '--phases',
4079                                      nargs='+',
4080                                      help='list of phases to add. phases are '
4081                                           'automatically associated with all '
4082                                           'histograms given.'),
4083
4084                         _args_kwargs('-pf', '--phaseformat',
4085                                      help='format hint for phase import. Applies to all'
4086                                           ' phases. Example: -pf CIF'
4087                                      ),
4088                                     
4089                         _args_kwargs('-l', '--histlist',
4090                                      nargs='+',
4091                                      help='list of histgram indices to associate with added'
4092                                           ' phases. If not specified, phases are'
4093                                           ' associated with all previously loaded'
4094                                           ' histograms. Example: -l 2 3 4')]),
4095                                           
4096               "dump": (dump, [_args_kwargs('-d', '--histograms',
4097                                     action='store_true',
4098                                     help='list histograms in files, overrides --raw'),
4099
4100                               _args_kwargs('-p', '--phases',
4101                                            action='store_true',
4102                                            help='list phases in files, overrides --raw'),
4103
4104                               _args_kwargs('-r', '--raw',
4105                                      action='store_true', help='dump raw file contents, default'),
4106
4107                               _args_kwargs('files', nargs='+')]),
4108
4109               "refine":
4110               (refine, [_args_kwargs('gpxfile', help='the project file to refine'),
4111                         _args_kwargs('refinements',
4112                                      help='JSON file of refinements to apply. if not present'
4113                                           ' refines file as-is',
4114                                      default=None,
4115                                      nargs='?')]),
4116
4117               "seqrefine": (seqrefine, []),
4118               "export": (export, [_args_kwargs('gpxfile',
4119                                                help='the project file from which to export'),
4120                                   _args_kwargs('phase', help='identifier of phase to export'),
4121                                   _args_kwargs('exportfile', help='the .cif file to export to')]),
4122               "browse": (IPyBrowse, [_args_kwargs('files', nargs='+',
4123                                                   help='list of files to browse')])}
4124
4125
4126def main():
4127    '''The command-line interface for calling GSASIIscriptable as a shell command,
4128    where it is expected to be called as::
4129
4130       python GSASIIscriptable.py <subcommand> <file.gpx> <options>
4131
4132    The following subcommands are defined:
4133
4134        * create, see :func:`create`
4135        * add, see :func:`add`
4136        * dump, see :func:`dump`
4137        * refine, see :func:`refine`
4138        * seqrefine, see :func:`seqrefine`
4139        * export, :func:`export`
4140        * browse, see :func:`IPyBrowse`
4141
4142    .. seealso::
4143        :func:`create`
4144        :func:`add`
4145        :func:`dump`
4146        :func:`refine`
4147        :func:`seqrefine`
4148        :func:`export`
4149        :func:`IPyBrowse`
4150    '''
4151    parser = argparse.ArgumentParser(description=
4152        "Use of "+os.path.split(__file__)[1]+" Allows GSAS-II actions from command line."
4153        )
4154    subs = parser.add_subparsers()
4155
4156    # Create all of the specified subparsers
4157    for name, (func, args) in subcommands.items():
4158        new_parser = subs.add_parser(name,help=commandhelp.get(name),
4159                                     description='Command "'+name+'" '+commandhelp.get(name))
4160        for listargs, kwds in args:
4161            new_parser.add_argument(*listargs, **kwds)
4162        new_parser.set_defaults(func=func)
4163
4164    # Parse and trigger subcommand
4165    result = parser.parse_args()
4166    result.func(result)
4167
4168if __name__ == '__main__':
4169    #fname='/tmp/corundum-template.gpx'
4170    #prj = G2Project(fname)
4171    main()
Note: See TracBrowser for help on using the repository browser.