source: trunk/GSASIIscriptable.py @ 3822

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

fixing a broken doc reference?

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