source: trunk/GSASIIscriptable.py @ 3205

Last change on this file since 3205 was 3205, checked in by toby, 4 years ago

expand scripting capabilities and documentation

File size: 116.6 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3########### SVN repository information ###################
4# $Date: 2017-04-12 15:12:45 -0500 (Wed, 12 Apr 2017) $
5# $Author: vondreele $
6# $Revision: 2777 $
7# $URL: https://subversion.xray.aps.anl.gov/pyGSAS/trunk/GSASIIscriptable.py $
8# $Id: GSASIIIO.py 2777 2017-04-12 20:12:45Z vondreele $
9########### SVN repository information ###################
10"""
11*GSASIIscriptable: Scripting Interface*
12=======================================
13
14Routines for reading, writing, modifying and creating GSAS-II project (.gpx) files.
15This file specifies several wrapper classes around GSAS-II data representations.
16They all inherit from :class:`G2ObjectWrapper`. The chief class is :class:`G2Project`,
17which represents an entire GSAS-II project and provides several methods to access
18phases, powder histograms, and execute Rietveld refinements. These routines are
19typically called by using the :ref:`CommandlineInterface` to access a number of features or
20the :ref:`API`, which allows much more versatile access.
21
22.. _Refinement_parameters_kinds:
23
24=====================
25Refinement parameters
26=====================
27
28Note that parameters and refinement flags used in GSAS-II fall into three classes:
29
30    * **Histogram**: There will be a set of these for each dataset loaded into a
31      project file. The parameters available depend on the type of histogram
32      (Bragg-Brentano, Single-Crystal, TOF,...). Typical Histogram parameters
33      include the overall scale factor, background, instrument and sample parameters;
34      see the :ref:`Histogram_parameters_table` table for a list of the histogram
35      parameters where access has been provided.
36     
37    * **Phase**: There will be a set of these for each phase loaded into a
38      project file. While some parameters are found in all types of phases,
39      others are only found in certain types (modulated, magnetic, protein...).
40      Typical phase parameters include unit cell lengths and atomic positions; see the
41      :ref:`Phase_parameters_table` table for a list of the phase     
42      parameters where access has been provided.
43     
44    * **Histogram-and-phase** (HAP): There is a set of these for every histogram
45      that is associated with each phase, so that if there are ``N`` phases and ``M``
46      histograms, there can be ``N*M`` total sets of "HAP" parameters sets (fewer if all
47      histograms are not linked to all phases.) Typical HAP parameters include the
48      phase fractions, sample microstrain and crystallite size broadening terms,
49      hydrostatic strain pertibations of the unit cell and preferred orientation
50      values.
51      See the :ref:`HAP_parameters_table` table for the HAP parameters where access has
52      been provided.
53
54There are several ways to set parameters using different objects, as described below.
55
56------------------------
57Histogram/Phase objects
58------------------------
59Each phase and powder histogram in a :class:`G2Project` object has an associated
60object. Parameters within each individual object can be turned on and off by calling
61:meth:`G2PwdrData.set_refinements` or :meth:`G2PwdrData.clear_refinements`
62for histogram parameters;
63:meth:`G2Phase.set_refinements` or :meth:`G2Phase.clear_refinements`
64for phase parameters; and :meth:`G2Phase.set_HAP_refinements` or
65: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:
66
67.. code-block::  python
68
69    params = { 'Limits': [0.8, 12.0],
70               'Sample Parameters': ['Absorption', 'Contrast', 'DisplaceX'],
71               'Background': {'type': 'chebyschev', 'refine': True}}
72    some_histogram.set_refinements(params)
73
74Likewise to turn refinement flags on, use code such as this:
75
76.. code-block::  python
77
78    params = { 'Instrument Parameters': ['U', 'V', 'W']}
79    some_histogram.set_refinements(params)
80
81and to turn these refinement flags, off use this (Note that the
82``.clear_refinements()`` methods will usually will turn off refinement even
83if a refinement parameter is set in the dict to True.):
84
85.. code-block::  python
86
87    params = { 'Instrument Parameters': ['U', 'V', 'W']}
88    some_histogram.clear_refinements(params)
89
90For phase parameters, use code such as this:
91   
92.. code-block::  python
93
94    params = { 'LeBail': True, 'Cell': True,
95               'Atoms': { 'Mn1': 'X',
96                          'O3': 'XU',
97                          'V4': 'FXU'}}
98    some_histogram.set_refinements(params)
99
100
101and here is an example for HAP parameters:
102
103.. code-block::  python
104
105    params = { 'Babinet': 'BabA',
106               'Extinction': True,
107               'Mustrain': { 'type': 'uniaxial',
108                             'direction': [0, 0, 1],
109                             'refine': True}}
110    some_phase.set_HAP_refinements(params)
111
112Note that the parameters must match the object type and method (phase vs. histogram vs. HAP).
113
114.. _Project_objects:
115
116-----------------------------
117Project-level Parameter Dict
118-----------------------------
119It is also possible to create a composite dictionary
120that references all three of the above types of refinement parameters.
121In this case dictionaries are nested with keys at the outer level, such as
122"set" and "clear" which determine function is used with function
123:meth:`G2Project.set_refinement`.
124
125Note that optionally a list of histograms and/or phases can be supplied to
126:meth:`G2Project.set_refinement` or :meth:`G2Project.do_refinements`,
127where the default is to use all phases and histograms, but more commonly for
128:meth:`G2Project.do_refinements` will be to define the "histograms" and "phases"
129items within individual dictionaries and these will override the call arguments.
130
131
132As an example:
133
134.. code-block::  python
135
136    pardict = {'set': { 'Limits': [0.8, 12.0],
137                       'Sample Parameters': ['Absorption', 'Contrast', 'DisplaceX'],
138                       'Background': {'type': 'chebyschev', 'refine': True}},
139              'clear': {'Instrument Parameters': ['U', 'V', 'W']}}
140    my_project.set_refinement(pardict)
141
142.. _Refinement_recipe:
143   
144------------------------
145Refinement recipe
146------------------------
147Finally, it is possible to specify a sequence of refinement actions as a list of
148dicts that contain :ref:`Project_objects` objects. This list is
149supplied as an argument to :meth:`G2Project.do_refinements`.
150These dicts in this list are each like those described in the
151:ref:`Project_objects` section,
152except that additional keys, as are described in the table below may be used.
153
154========== ============================================================================
155key         explanation
156========== ============================================================================
157set                    Specifies a dict with keys and subkeys as described in the
158                       :ref:`Refinement_parameters_fmt` section. Items listed here
159                       will be set to be refined.
160clear                  Specifies a dict as above for set, except that parameters are
161                       cleared and thus will not be refined.
162once                   Specifies a dict as above for set, except that parameters are
163                       set for the next cycle of refinement and are cleared once the
164                       refinement step is completed.
165skip                   Normally, once parameters are processed with a set/clear/once
166                       action(s), a refinement is started. If skip is defined as True
167                       (or any other value) the refinement step is not performed.
168output                 If a file name is specified for output is will be used for
169                       the current refinement.
170histograms             Should contain a list of histogram(s) to be used for the
171                       set/clear/once action(s) on :ref:`Histogram_parameters_table` or
172                       :ref:`HAP_parameters_table`. Note that this will be
173                       ignored for :ref:`Phase_parameters_table`. Histograms may be
174                       specified as a list of strings [('PWDR ...'),...], indices
175                       [0,1,2] or as list of objects [hist1, hist2].
176phases                 Should contain a list of phase(s) to be used for the
177                       set/clear/once action(s) on :ref:`Phase_parameters_table` or
178                       :ref:`HAP_parameters_table`. Note that this will be
179                       ignored for :ref:`Histogram_parameters_table`.
180                       Phases may be specified as a list of strings
181                       [('Phase name'),...], indices [0,1,2] or as list of objects
182                       [phase0, phase2].
183call                   Specifies a function to call after a refinement is completed.
184                       No function is called if this is not specified.
185callargs               Provides a list of arguments that will be passed to the function
186                       in call (if any). If call is defined and callargs is not, the
187                       current <tt>G2Project</tt> is passed as a single argument.
188========== ============================================================================
189
190An example follows:
191
192.. code-block::  python
193
194    reflist = [
195            {"set": { "Limits": { "low": 0.7 },
196                      "Background": { "no. coeffs": 3,
197                                      "refine": True }}},
198            {"set": { "LeBail": True,
199                      "Cell": True }},
200            {"set": { "Sample Parameters": ["DisplaceX"]}},
201            {"set": { "Instrument Parameters": ["U", "V", "W", "X", "Y"]}},
202            {"set": { "Mustrain": { "type": "uniaxial",
203                                    "refine": "equatorial",
204                                    "direction": [0, 0, 1]}}},
205            {"set": { "Mustrain": { "type": "uniaxial",
206                                    "refine": "axial"}}},
207            {"clear": { "LeBail": True},
208             "set": { "Atoms": { "Mn": "X" }}},
209            {"set": { "Atoms": { "O1": "X", "O2": "X" }}},]
210    my_project.do_refinements(reflist)
211   
212
213In this example, the list contains a set of dicts, each defined as before.
214A separate refinement step will be performed for each element in the list unless
215"skip" is included.
216Note that in the second from last refinement step, parameters are both set and cleared.
217
218.. _Refinement_parameters_fmt:
219
220============================
221Refinement specifiers format
222============================
223
224Refinement parameters are specified as dictionaries, supplied to any of the functions
225named in :ref:`Refinement_parameters_kinds`. Each method accepts a different set
226of keys, described below for each of the three parameter classes.
227
228.. _Histogram_parameters_table:
229
230--------------------
231Histogram parameters
232--------------------
233
234This table describes the dictionaries supplied to :func:`G2PwdrData.set_refinements`
235and :func:`G2PwdrData.clear_refinements`.
236
237.. tabularcolumns:: |l|l|p{3.5in}|
238
239===================== ====================  =================================================
240key                   subkey                explanation
241===================== ====================  =================================================
242Limits                                      The 2-theta range of values to consider. Can
243                                            be either a dictionary of 'low' and/or 'high',
244                                            or a list of 2 items [low, high]
245\                     low                   Sets the low limit
246\                     high                  Sets the high limit
247Sample Parameters                           Should be provided as a **list** of subkeys
248                                            to set or clear, e.g. ['DisplaceX', 'Scale']
249\                     Absorption
250\                     Contrast
251\                     DisplaceX             Sample displacement along the X direction
252\                     DisplaceY             Sample displacement along the Y direction
253\                     Scale                 Histogram Scale factor
254Background                                  Sample background. If value is a boolean,
255                                            the background's 'refine' parameter is set
256                                            to the given boolean. Usually should be a
257                                            dictionary with any of the following keys:
258\                     type                  The background model, e.g. 'chebyschev'
259\                     refine                The value of the refine flag, boolean
260\                     no. coeffs            Number of coefficients to use, integer
261\                     coeffs                List of floats, literal values for background
262\                     FixedPoints           List of (2-theta, intensity) values for fixed points
263\                     fit fixed points      If True, triggers a fit to the fixed points to be calculated. It is calculated when this key is detected, regardless of calls to refine.
264Instrument Parameters                       As in Sample Paramters, Should be provided as a **list** of subkeys to
265                                            set or clear, e.g. ['X', 'Y', 'Zero', 'SH/L']
266\                     U, V, W               All separate keys. Gaussian peak profile terms
267\                     X, Y                  Separate keys. Lorentzian peak profile terms
268\                     Zero                  Zero shift
269\                     SH/L
270\                     Polariz.              Polarization parameter
271\                     Lam                   Lambda, the incident wavelength
272===================== ====================  =================================================
273
274.. _Phase_parameters_table:
275
276----------------
277Phase parameters
278----------------
279
280This table describes the dictionaries supplied to :func:`G2Phase.set_refinements`
281and :func:`G2Phase.clear_refinements`.
282
283.. tabularcolumns:: |l|p{4.5in}|
284
285======= ==========================================================
286key                   explanation
287======= ==========================================================
288Cell                  Whether or not to refine the unit cell.
289Atoms                 Dictionary of atoms and refinement flags.
290                      Each key should be an atom label, e.g.
291                      'O3', 'Mn5', and each value should be
292                      a string defining what values to refine.
293                      Values can be any combination of 'F'
294                      for fractional occupancy, 'X' for position,
295                      and 'U' for Debye-Waller factor
296LeBail                Enables LeBail intensity extraction.
297======= ==========================================================
298
299
300.. _HAP_parameters_table:
301
302------------------------------
303Histogram-and-phase parameters
304------------------------------
305
306This table describes the dictionaries supplied to :func:`G2Phase.set_HAP_refinements`
307and :func:`G2Phase.clear_HAP_refinements`.
308
309.. tabularcolumns:: |l|l|p{3.5in}|
310
311=============  ==========  ============================================================
312key             subkey                 explanation
313=============  ==========  ============================================================
314Babinet                                Should be a **list** of the following
315                                       subkeys. If not, assumes both
316                                       BabA and BabU
317\               BabA
318\               BabU
319Extinction                             Boolean, True to refine.
320HStrain                                Boolean, True to refine all appropriate
321                                       $D_ij$ terms.
322Mustrain
323\               type                   Mustrain model. One of 'isotropic',
324                                       'uniaxial', or 'generalized'. Should always
325                                       be specified.
326\              direction               For uniaxial only. A list of three
327                                       integers,
328                                       the [hkl] direction of the axis.
329\               refine                 Usually boolean, set to True to refine.
330                                       When in doubt, set it to true.
331                                       For uniaxial model, can specify list
332                                       of 'axial' or 'equatorial' or a single
333                                       boolean sets both axial and equatorial.
334Size                                   Not yet implemented
335\               type                   Size broadening model. One of 'isotropic',
336                                       'uniaxial', or 'ellipsoid'. Should always
337                                       be specified.
338\              direction               For uniaxial only. A list of three
339                                       integers,
340                                       the [hkl] direction of the axis.
341\               refine                 A boolean, True to refine.
342Pref.Ori.                              Boolean, True to refine
343Show                                   Boolean, True to refine
344Use                                    Boolean, True to refine
345Scale                                  Phase fraction; Boolean, True to refine
346=============  ==========  ============================================================
347
348.. _CommandlineInterface:
349
350=======================================
351GSASIIscriptable Command-line Interface
352=======================================
353
354One way to access these routines is by calling this script
355via a command line interface as a shell command, where it is expected to be called as::
356
357       python GSASIIscriptable.py <subcommand> <file.gpx> <options>
358
359    The following subcommands are defined:
360
361        * create, see :func:`create`
362        * add, see :func:`add`
363        * dump, see :func:`dump`
364        * refine, see :func:`refine`
365        * seqrefine, see :func:`seqrefine`
366        * export, :func:`export`
367        * browse, see :func:`IPyBrowse`
368       
369Run::
370
371   python GSASIIscriptable.py --help
372
373to show the available subcommands, and inspect each subcommand with
374`python GSASIIscriptable.py <subcommand> --help` or see the documentation for each of the above routines.
375
376.. _JsonFormat:
377
378-------------------------
379Parameters in JSON files
380-------------------------
381
382The refine command requires two inputs: an existing GSAS-II project (.gpx) file and
383a JSON format file
384(see `Introducing JSON <http://json.org/>`_) that contains a single dict.
385This dict may have two keys:
386
387refinements:
388  This defines the a set of refinement steps in a JSON representation of a
389  :ref:`Refinement_recipe` list.
390
391code:
392  This optionally defines Python code that will be executed after the project is loaded,
393  but before the refinement is started. This can be used to execute Python code to change
394  parameters that are not accessible via a :ref:`Refinement_recipe` dict (note that the
395  project object is accessed with variable ``proj``) or to define code that will be called
396  later (see key ``call`` in the :ref:`Refinement_recipe` section.)
397   
398JSON website: `Introducing JSON <http://json.org/>`_.
399
400.. _API:
401
402===================================
403GSASIIscriptable Application Layer
404===================================
405
406This module provides a large number of classes and modules, as described below.
407Most commonly a script will create a G2Project object using :class:`G2Project` and then
408perform actions such as
409adding a histogram (method :meth:`G2Project.add_powder_histogram`),
410adding a phase (method :meth:`G2Project.add_phase`),
411or setting parameters and performing a refinement
412(method :meth:`G2Project.do_refinements`).
413
414In some cases, it may be easier or more options may be available by direct access to
415methods inside :class:`G2PwdrData` or :class:`G2Phase`
416
417---------------------------------------
418GSASIIscriptable Classes and functions
419---------------------------------------
420"""
421from __future__ import division, print_function
422import argparse
423import os.path as ospath
424import datetime as dt
425import sys
426import platform
427if '2' in platform.python_version_tuple()[0]:
428    import cPickle
429    strtypes = (str,unicode)
430else:
431    import _pickle as cPickle
432    strtypes = (str,bytes)
433import imp
434import copy
435import os
436import random as ran
437
438import numpy.ma as ma
439import scipy.interpolate as si
440import numpy as np
441import scipy as sp
442
443import GSASIIpath
444GSASIIpath.SetBinaryPath(True)  # for now, this is needed before some of these modules can be imported
445import GSASIIobj as G2obj
446import GSASIIpwd as G2pwd
447import GSASIIstrMain as G2strMain
448import GSASIIspc as G2spc
449import GSASIIElem as G2elem
450
451
452# Delay imports to not slow down small scripts
453G2fil = None
454PwdrDataReaders = []
455PhaseReaders = []
456
457def LoadG2fil():
458    """Delay importing this module, it is slow"""
459    global G2fil
460    global PwdrDataReaders
461    global PhaseReaders
462    if G2fil is None:
463        import GSASIIfiles
464        G2fil = GSASIIfiles
465        PwdrDataReaders = G2fil.LoadImportRoutines("pwd", "Powder_Data")
466        PhaseReaders = G2fil.LoadImportRoutines("phase", "Phase")
467
468
469def LoadDictFromProjFile(ProjFile):
470    '''Read a GSAS-II project file and load items to dictionary
471   
472    :param str ProjFile: GSAS-II project (name.gpx) full file name
473    :returns: Project,nameList, where
474
475      * Project (dict) is a representation of gpx file following the GSAS-II tree structure
476        for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict
477        data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)}
478      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
479
480    Example for fap.gpx::
481
482      Project = {                 #NB:dict order is not tree order
483        u'Phases':{'data':None,'fap':{phase dict}},
484        u'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc},
485        u'Rigid bodies':{'data': {rigid body dict}},
486        u'Covariance':{'data':{covariance data dict}},
487        u'Controls':{'data':{controls data dict}},
488        u'Notebook':{'data':[notebook list]},
489        u'Restraints':{'data':{restraint data dict}},
490        u'Constraints':{'data':{constraint data dict}}]
491        }
492      nameList = [                #NB: reproduces tree order
493        [u'Notebook',],
494        [u'Controls',],
495        [u'Covariance',],
496        [u'Constraints',],
497        [u'Restraints',],
498        [u'Rigid bodies',],
499        [u'PWDR FAP.XRA Bank 1',
500             u'Comments',
501             u'Limits',
502             u'Background',
503             u'Instrument Parameters',
504             u'Sample Parameters',
505             u'Peak List',
506             u'Index Peak List',
507             u'Unit Cells List',
508             u'Reflection Lists'],
509        [u'Phases', u'fap']
510        ]
511    '''
512    # Let IOError be thrown if file does not exist
513    # if not ospath.exists(ProjFile):
514    #     print ('\n*** Error attempt to open project file that does not exist:\n   '+
515    #         str(ProjFile))
516    #     return
517    file = open(ProjFile,'rb')
518    # print('loading from file: {}'.format(ProjFile))
519    Project = {}
520    nameList = []
521    try:
522        while True:
523            try:
524                data = cPickle.load(file)
525            except EOFError:
526                break
527            datum = data[0]
528            Project[datum[0]] = {'data':datum[1]}
529            nameList.append([datum[0],])
530            for datus in data[1:]:
531                Project[datum[0]][datus[0]] = datus[1]
532                nameList[-1].append(datus[0])
533        file.close()
534        # print('project load successful')
535    except:
536        raise IOError("Error reading file "+str(ProjFile)+". This is not a GSAS-II .gpx file")
537    finally:
538        file.close()
539    return Project,nameList
540
541def SaveDictToProjFile(Project,nameList,ProjFile):
542    '''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile
543
544    :param dict Project: representation of gpx file following the GSAS-II
545        tree structure as described for LoadDictFromProjFile
546    :param list nameList: names of main tree entries & subentries used to reconstruct project file
547    :param str ProjFile: full file name for output project.gpx file (including extension)
548    '''
549    file = open(ProjFile,'wb')
550    try:
551        for name in nameList:
552            data = []
553            item = Project[name[0]]
554            data.append([name[0],item['data']])
555            for item2 in name[1:]:
556                data.append([item2,item[item2]])
557            cPickle.dump(data,file,1)
558    finally:
559        file.close()
560
561def ImportPowder(reader,filename):
562    '''Use a reader to import a powder diffraction data file
563
564    :param str reader: a scriptable reader
565    :param str filename: full name of powder data file; can be "multi-Bank" data
566
567    :returns: list rdlist: list of reader objects containing powder data, one for each
568        "Bank" of data encountered in file. Items in reader object of interest are:
569
570          * rd.comments: list of str: comments found on powder file
571          * rd.dnames: list of str: data nammes suitable for use in GSASII data tree NB: duplicated in all rd entries in rdlist
572          * rd.powderdata: list of numpy arrays: pos,int,wt,zeros,zeros,zeros as needed for a PWDR entry in  GSASII data tree.
573    '''
574    rdfile,rdpath,descr = imp.find_module(reader)
575    rdclass = imp.load_module(reader,rdfile,rdpath,descr)
576    rd = rdclass.GSAS_ReaderClass()
577    if not rd.scriptable:
578        print(u'**** ERROR: '+reader+u' is not a scriptable reader')
579        return None
580    rdlist = []
581    if rd.ContentsValidator(filename):
582        repeat = True
583        rdbuffer = {} # create temporary storage for file reader
584        block = 0
585        while repeat: # loop if the reader asks for another pass on the file
586            block += 1
587            repeat = False
588            rd.objname = ospath.basename(filename)
589            flag = rd.Reader(filename,None,buffer=rdbuffer,blocknum=block,)
590            if flag:
591                rdlist.append(copy.deepcopy(rd)) # save the result before it is written over
592                if rd.repeat:
593                    repeat = True
594        return rdlist
595    print(rd.errors)
596    return None
597
598def SetDefaultDData(dType,histoName,NShkl=0,NDij=0):
599    '''Create an initial Histogram dictionary
600
601    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
602    '''
603    if dType in ['SXC','SNC']:
604        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
605            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
606            'Extinction':['Lorentzian','None', {'Tbar':0.1,'Cos2TM':0.955,
607            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}],
608            'Flack':[0.0,False]}
609    elif dType == 'SNT':
610        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
611            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
612            'Extinction':['Lorentzian','None', {
613            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}]}
614    elif 'P' in dType:
615        return {'Histogram':histoName,'Show':False,'Scale':[1.0,False],
616            'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{},[],0.1],
617            'Size':['isotropic',[1.,1.,1.],[False,False,False],[0,0,1],
618                [1.,1.,1.,0.,0.,0.],6*[False,]],
619            'Mustrain':['isotropic',[1000.0,1000.0,1.0],[False,False,False],[0,0,1],
620                NShkl*[0.01,],NShkl*[False,]],
621            'HStrain':[NDij*[0.0,],NDij*[False,]],
622            'Extinction':[0.0,False],'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}}
623
624
625def PreSetup(data):
626    '''Create part of an initial (empty) phase dictionary
627
628    from GSASIIphsGUI.py, near end of UpdatePhaseData
629
630    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
631    '''
632    if 'RBModels' not in data:
633        data['RBModels'] = {}
634    if 'MCSA' not in data:
635        data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[.8,1.2],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
636    if 'dict' in str(type(data['MCSA']['Results'])):
637        data['MCSA']['Results'] = []
638    if 'Modulated' not in data['General']:
639        data['General']['Modulated'] = False
640    if 'modulated' in data['General']['Type']:
641        data['General']['Modulated'] = True
642        data['General']['Type'] = 'nuclear'
643
644
645def SetupGeneral(data, dirname):
646    """Helps initialize phase data.
647
648    From GSASIIphsGui.py, function of the same name. Minor changes for imports etc.
649
650    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
651    """
652    mapDefault = {'MapType':'','RefList':'','Resolution':0.5,'Show bonds':True,
653                'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False}
654    generalData = data['General']
655    atomData = data['Atoms']
656    generalData['AtomTypes'] = []
657    generalData['Isotopes'] = {}
658
659    if 'Isotope' not in generalData:
660        generalData['Isotope'] = {}
661    if 'Data plot type' not in generalData:
662        generalData['Data plot type'] = 'Mustrain'
663    if 'POhkl' not in generalData:
664        generalData['POhkl'] = [0,0,1]
665    if 'Map' not in generalData:
666        generalData['Map'] = mapDefault.copy()
667    if 'Flip' not in generalData:
668        generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None',
669            'k-factor':0.1,'k-Max':20.,}
670    if 'testHKL' not in generalData['Flip']:
671        generalData['Flip']['testHKL'] = [[0,0,2],[2,0,0],[1,1,1],[0,2,0],[1,2,3]]
672    if 'doPawley' not in generalData:
673        generalData['doPawley'] = False     #ToDo: change to ''
674    if 'Pawley dmin' not in generalData:
675        generalData['Pawley dmin'] = 1.0
676    if 'Pawley neg wt' not in generalData:
677        generalData['Pawley neg wt'] = 0.0
678    if 'Algolrithm' in generalData.get('MCSA controls',{}) or \
679        'MCSA controls' not in generalData:
680        generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50],
681        'dmin':2.0,'Algorithm':'log','Jump coeff':[0.95,0.5],'boltzmann':1.0,
682        'fast parms':[1.0,1.0,1.0],'log slope':0.9,'Cycles':1,'Results':[],'newDmin':True}
683    if 'AtomPtrs' not in generalData:
684        generalData['AtomPtrs'] = [3,1,7,9]
685        if generalData['Type'] == 'macromolecular':
686            generalData['AtomPtrs'] = [6,4,10,12]
687        elif generalData['Type'] == 'magnetic':
688            generalData['AtomPtrs'] = [3,1,10,12]
689    if generalData['Type'] in ['modulated',]:
690        generalData['Modulated'] = True
691        generalData['Type'] = 'nuclear'
692        if 'Super' not in generalData:
693            generalData['Super'] = 1
694            generalData['SuperVec'] = [[0,0,.1],False,4]
695            generalData['SSGData'] = {}
696        if '4DmapData' not in generalData:
697            generalData['4DmapData'] = mapDefault.copy()
698            generalData['4DmapData'].update({'MapType':'Fobs'})
699    if 'Modulated' not in generalData:
700        generalData['Modulated'] = False
701    if 'HydIds' not in generalData:
702        generalData['HydIds'] = {}
703    cx,ct,cs,cia = generalData['AtomPtrs']
704    generalData['NoAtoms'] = {}
705    generalData['BondRadii'] = []
706    generalData['AngleRadii'] = []
707    generalData['vdWRadii'] = []
708    generalData['AtomMass'] = []
709    generalData['Color'] = []
710    if generalData['Type'] == 'magnetic':
711        generalData['MagDmin'] = generalData.get('MagDmin',1.0)
712        landeg = generalData.get('Lande g',[])
713    generalData['Mydir'] = dirname
714    badList = {}
715    for iat,atom in enumerate(atomData):
716        atom[ct] = atom[ct].lower().capitalize()              #force to standard form
717        if generalData['AtomTypes'].count(atom[ct]):
718            generalData['NoAtoms'][atom[ct]] += atom[cx+3]*float(atom[cs+1])
719        elif atom[ct] != 'UNK':
720            Info = G2elem.GetAtomInfo(atom[ct])
721            if not Info:
722                if atom[ct] not in badList:
723                    badList[atom[ct]] = 0
724                badList[atom[ct]] += 1
725                atom[ct] = 'UNK'
726                continue
727            atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo
728            generalData['AtomTypes'].append(atom[ct])
729            generalData['Z'] = Info['Z']
730            generalData['Isotopes'][atom[ct]] = Info['Isotopes']
731            generalData['BondRadii'].append(Info['Drad'])
732            generalData['AngleRadii'].append(Info['Arad'])
733            generalData['vdWRadii'].append(Info['Vdrad'])
734            if atom[ct] in generalData['Isotope']:
735                if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]:
736                    isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
737                    generalData['Isotope'][atom[ct]] = isotope
738                generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][atom[ct]]]['Mass'])
739            else:
740                generalData['Isotope'][atom[ct]] = 'Nat. Abund.'
741                if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]:
742                    isotope = list(generalData['Isotopes'][atom[ct]].keys())[-1]
743                    generalData['Isotope'][atom[ct]] = isotope
744                generalData['AtomMass'].append(Info['Mass'])
745            generalData['NoAtoms'][atom[ct]] = atom[cx+3]*float(atom[cs+1])
746            generalData['Color'].append(Info['Color'])
747            if generalData['Type'] == 'magnetic':
748                if len(landeg) < len(generalData['AtomTypes']):
749                    landeg.append(2.0)
750    if generalData['Type'] == 'magnetic':
751        generalData['Lande g'] = landeg[:len(generalData['AtomTypes'])]
752
753    if badList:
754        msg = 'Warning: element symbol(s) not found:'
755        for key in badList:
756            msg += '\n\t' + key
757            if badList[key] > 1:
758                msg += ' (' + str(badList[key]) + ' times)'
759        raise G2ScriptException("Phase error:\n" + msg)
760        # wx.MessageBox(msg,caption='Element symbol error')
761    F000X = 0.
762    F000N = 0.
763    for i,elem in enumerate(generalData['AtomTypes']):
764        F000X += generalData['NoAtoms'][elem]*generalData['Z']
765        isotope = generalData['Isotope'][elem]
766        F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope]['SL'][0]
767    generalData['F000X'] = F000X
768    generalData['F000N'] = F000N
769    import GSASIImath as G2mth
770    generalData['Mass'] = G2mth.getMass(generalData)
771
772
773def make_empty_project(author=None, filename=None):
774    """Creates an dictionary in the style of GSASIIscriptable, for an empty
775    project.
776
777    If no author name or filename is supplied, 'no name' and
778    <current dir>/test_output.gpx are used , respectively.
779
780    Returns: project dictionary, name list
781
782    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
783    """
784    if not filename:
785        filename = 'test_output.gpx'
786    filename = os.path.abspath(filename)
787    gsasii_version = str(GSASIIpath.GetVersionNumber())
788    LoadG2fil()
789    import matplotlib as mpl
790    python_library_versions = G2fil.get_python_versions([mpl, np, sp])
791
792    controls_data = dict(G2obj.DefaultControls)
793    controls_data['LastSavedAs'] = filename
794    controls_data['LastSavedUsing'] = gsasii_version
795    controls_data['PythonVersions'] = python_library_versions
796    if author:
797        controls_data['Author'] = author
798
799    output = {'Constraints': {'data': {'HAP': [], 'Hist': [], 'Phase': [],
800                                       'Global': []}},
801              'Controls': {'data': controls_data},
802              u'Covariance': {'data': {}},
803              u'Notebook': {'data': ['']},
804              u'Restraints': {'data': {}},
805              u'Rigid bodies': {'data': {'RBIds': {'Residue': [], 'Vector': []},
806                                'Residue': {'AtInfo': {}},
807                                'Vector':  {'AtInfo': {}}}}}
808
809    names = [[u'Notebook'], [u'Controls'], [u'Covariance'],
810             [u'Constraints'], [u'Restraints'], [u'Rigid bodies']]
811
812    return output, names
813
814
815class G2ImportException(Exception):
816    pass
817
818class G2ScriptException(Exception):
819    pass
820
821def import_generic(filename, readerlist, fmthint=None):
822    """Attempt to import a filename, using a list of reader objects.
823
824    Returns the first reader object which worked."""
825    # Translated from OnImportGeneric method in GSASII.py
826    primaryReaders, secondaryReaders = [], []
827    for reader in readerlist:
828        if fmthint is not None and fmthint not in reader.formatName: continue
829        flag = reader.ExtensionValidator(filename)
830        if flag is None:
831            secondaryReaders.append(reader)
832        elif flag:
833            primaryReaders.append(reader)
834    if not secondaryReaders and not primaryReaders:
835        raise G2ImportException("Could not read file: ", filename)
836
837    with open(filename, 'Ur') as fp:
838        rd_list = []
839
840        for rd in primaryReaders + secondaryReaders:
841            # Initialize reader
842            rd.selections = []
843            rd.dnames = []
844            rd.ReInitialize()
845            # Rewind file
846            rd.errors = ""
847            if not rd.ContentsValidator(filename):
848                # Report error
849                pass
850            if len(rd.selections) > 1:
851                # Select data?
852                # GSASII.py:543
853                raise G2ImportException("Not sure what data to select")
854
855            block = 0
856            rdbuffer = {}
857            repeat = True
858            while repeat:
859                repeat = False
860                block += 1
861                rd.objname = os.path.basename(filename)
862                try:
863                    flag = rd.Reader(filename,buffer=rdbuffer, blocknum=block)
864                except:
865                    flag = False
866                if flag:
867                    # Omitting image loading special cases
868                    rd.readfilename = filename
869                    rd_list.append(copy.deepcopy(rd))
870                    repeat = rd.repeat
871                else:
872                    if GSASIIpath.GetConfigValue('debug'): print("{} Reader failed to read {}".format(rd.formatName,filename))
873            if rd_list:
874                if rd.warnings:
875                    print("Read warning by", rd.formatName, "reader:",
876                          rd.warnings, file=sys.stderr)
877                else:
878                    print("{} read by Reader {}\n".format(filename,rd.formatName))                   
879                return rd_list
880    raise G2ImportException("No reader could read file: " + filename)
881
882
883def load_iprms(instfile, reader):
884    """Loads instrument parameters from a file, and edits the
885    given reader.
886
887    Returns a 2-tuple of (Iparm1, Iparm2) parameters
888    """
889    LoadG2fil()
890    ext = os.path.splitext(instfile)[1]
891
892    if ext.lower() == '.instprm':
893        # New GSAS File, load appropriately
894        with open(instfile) as f:
895            lines = f.readlines()
896        bank = reader.powderentry[2]
897        numbanks = reader.numbanks
898        iparms = G2fil.ReadPowderInstprm(lines, bank, numbanks, reader)
899
900        reader.instfile = instfile
901        reader.instmsg = 'GSAS-II file' + instfile
902        return iparms
903    elif ext.lower() not in ('.prm', '.inst', '.ins'):
904        raise ValueError('Expected .prm file, found: ', instfile)
905
906    # It's an old GSAS file, load appropriately
907    Iparm = {}
908    with open(instfile, 'Ur') as fp:
909        for line in fp:
910            if '#' in line:
911                continue
912            Iparm[line[:12]] = line[12:-1]
913    ibanks = int(Iparm.get('INS   BANK  ', '1').strip())
914    if ibanks == 1:
915        reader.instbank = 1
916        reader.powderentry[2] = 1
917        reader.instfile = instfile
918        reader.instmsg = instfile + ' bank ' + str(reader.instbank)
919        return G2fil.SetPowderInstParms(Iparm, reader)
920    # TODO handle >1 banks
921    raise NotImplementedError("Check GSASIIfiles.py: ReadPowderInstprm")
922
923def load_pwd_from_reader(reader, instprm, existingnames=[]):
924    """Loads powder data from a reader object, and assembles it into a GSASII data tree.
925
926    :returns: (name, tree) - 2-tuple of the histogram name (str), and data
927
928    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
929    """
930    HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_')
931    HistName = G2obj.MakeUniqueLabel(HistName, existingnames)
932
933    try:
934        Iparm1, Iparm2 = instprm
935    except ValueError:
936        Iparm1, Iparm2 = load_iprms(instprm, reader)
937
938    Ymin = np.min(reader.powderdata[1])
939    Ymax = np.max(reader.powderdata[1])
940    valuesdict = {'wtFactor': 1.0,
941                  'Dummy': False,
942                  'ranId': ran.randint(0, sys.maxsize),
943                  'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax,
944                  'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax,
945                  'qPlot': False, 'dPlot': False, 'sqrtPlot': False,
946                  'Yminmax': [Ymin, Ymax]}
947    reader.Sample['ranId'] = valuesdict['ranId']
948
949    # Ending keys:
950    # [u'Reflection Lists',
951    #  u'Limits',
952    #  'data',
953    #  u'Index Peak List',
954    #  u'Comments',
955    #  u'Unit Cells List',
956    #  u'Sample Parameters',
957    #  u'Peak List',
958    #  u'Background',
959    #  u'Instrument Parameters']
960    Tmin = np.min(reader.powderdata[0])
961    Tmax = np.max(reader.powderdata[0])
962
963    default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0],
964                          {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}]
965
966    output_dict = {u'Reflection Lists': {},
967                   u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]),
968                   u'data': [valuesdict, reader.powderdata, HistName],
969                   u'Index Peak List': [[], []],
970                   u'Comments': reader.comments,
971                   u'Unit Cells List': [],
972                   u'Sample Parameters': reader.Sample,
973                   u'Peak List': {'peaks': [], 'sigDict': {}},
974                   u'Background': reader.pwdparms.get('Background', default_background),
975                   u'Instrument Parameters': [Iparm1, Iparm2],
976                   }
977
978    names = [u'Comments',
979             u'Limits',
980             u'Background',
981             u'Instrument Parameters',
982             u'Sample Parameters',
983             u'Peak List',
984             u'Index Peak List',
985             u'Unit Cells List',
986             u'Reflection Lists']
987
988    # TODO controls?? GSASII.py:1664-7
989
990    return HistName, [HistName] + names, output_dict
991
992
993def _deep_copy_into(from_, into):
994    """Helper function for reloading .gpx file. See G2Project.reload()
995
996    :author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
997    """
998    if isinstance(from_, dict) and isinstance(into, dict):
999        combined_keys = set(from_.keys()).union(into.keys())
1000        for key in combined_keys:
1001            if key in from_ and key in into:
1002                both_dicts = (isinstance(from_[key], dict)
1003                              and isinstance(into[key], dict))
1004                both_lists = (isinstance(from_[key], list)
1005                              and isinstance(into[key], list))
1006                if both_dicts or both_lists:
1007                    _deep_copy_into(from_[key], into[key])
1008                else:
1009                    into[key] = from_[key]
1010            elif key in from_:
1011                into[key] = from_[key]
1012            else:  # key in into
1013                del into[key]
1014    elif isinstance(from_, list) and isinstance(into, list):
1015        if len(from_) == len(into):
1016            for i in range(len(from_)):
1017                both_dicts = (isinstance(from_[i], dict)
1018                              and isinstance(into[i], dict))
1019                both_lists = (isinstance(from_[i], list)
1020                              and isinstance(into[i], list))
1021                if both_dicts or both_lists:
1022                    _deep_copy_into(from_[i], into[i])
1023                else:
1024                    into[i] = from_[i]
1025        else:
1026            into[:] = from_
1027
1028
1029class G2ObjectWrapper(object):
1030    """Base class for all GSAS-II object wrappers.
1031
1032    The underlying GSAS-II format can be accessed as `wrapper.data`. A number
1033    of overrides are implemented so that the wrapper behaves like a dictionary.
1034
1035    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1036    """
1037    def __init__(self, datadict):
1038        self.data = datadict
1039
1040    def __getitem__(self, key):
1041        return self.data[key]
1042
1043    def __setitem__(self, key, value):
1044        self.data[key] = value
1045
1046    def __contains__(self, key):
1047        return key in self.data
1048
1049    def get(self, k, d=None):
1050        return self.data.get(k, d)
1051
1052    def keys(self):
1053        return self.data.keys()
1054
1055    def values(self):
1056        return self.data.values()
1057
1058    def items(self):
1059        return self.data.items()
1060
1061
1062class G2Project(G2ObjectWrapper):
1063    """
1064    Represents an entire GSAS-II project.
1065
1066    There are two ways to initialize it:
1067
1068    >>> # Load an existing project file
1069    >>> proj = G2Project('filename.gpx')
1070    >>> # Create a new project
1071    >>> proj = G2Project(filename='new_file.gpx')
1072    >>> # Specify an author
1073    >>> proj = G2Project(author='Dr. So-And-So', filename='my_data.gpx')
1074
1075    Histograms can be accessed easily.
1076
1077    >>> # By name
1078    >>> hist = proj.histogram('PWDR my-histogram-name')
1079    >>> # Or by index
1080    >>> hist = proj.histogram(0)
1081    >>> assert hist.id == 0
1082    >>> # Or by random id
1083    >>> assert hist == proj.histogram(hist.ranId)
1084
1085    Phases can be accessed the same way.
1086
1087    >>> phase = proj.phase('name of phase')
1088
1089    New data can also be loaded via :meth:`~G2Project.add_phase` and
1090    :meth:`~G2Project.add_powder_histogram`.
1091
1092    >>> hist = proj.add_powder_histogram('some_data_file.chi',
1093                                         'instrument_parameters.prm')
1094    >>> phase = proj.add_phase('my_phase.cif', histograms=[hist])
1095
1096    Parameters for Rietveld refinement can be turned on and off as well.
1097    See :meth:`~G2Project.set_refinement`, :meth:`~G2Project.refine`,
1098    :meth:`~G2Project.iter_refinements`, :meth:`~G2Project.do_refinements`.
1099    """
1100    def __init__(self, gpxfile=None, author=None, filename=None):
1101        """Loads a GSAS-II project from a specified filename.
1102
1103        :param str gpxfile: Existing .gpx file to be loaded. If nonexistent,
1104            creates an empty project.
1105        :param str author: Author's name. Optional.
1106        :param str filename: The filename the project should be saved to in
1107            the future. If both filename and gpxfile are present, the project is
1108            loaded from the gpxfile, then set to  save to filename in the future"""
1109        if gpxfile is None:
1110            filename = os.path.abspath(os.path.expanduser(filename))
1111            self.filename = filename
1112            self.data, self.names = make_empty_project(author=author, filename=filename)
1113        elif isinstance(gpxfile, str):
1114            # TODO set author, filename
1115            self.filename = os.path.abspath(os.path.expanduser(gpxfile))
1116            self.data, self.names = LoadDictFromProjFile(gpxfile)
1117            self.update_ids()
1118        else:
1119            raise ValueError("Not sure what to do with gpxfile")
1120
1121    @classmethod
1122    def from_dict_and_names(cls, gpxdict, names, filename=None):
1123        """Creates a :class:`G2Project` directly from
1124        a dictionary and a list of names. If in doubt, do not use this.
1125
1126        :returns: a :class:`G2Project`
1127        """
1128        out = cls()
1129        if filename:
1130            filename = os.path.abspath(os.path.expanduser(filename))
1131            out.filename = filename
1132            gpxdict['Controls']['data']['LastSavedAs'] = filename
1133        else:
1134            try:
1135                out.filename = gpxdict['Controls']['data']['LastSavedAs']
1136            except KeyError:
1137                out.filename = None
1138        out.data = gpxdict
1139        out.names = names
1140
1141    def save(self, filename=None):
1142        """Saves the project, either to the current filename, or to a new file.
1143
1144        Updates self.filename if a new filename provided"""
1145        # TODO update LastSavedUsing ?
1146        if filename:
1147            filename = os.path.abspath(os.path.expanduser(filename))
1148            self.data['Controls']['data']['LastSavedAs'] = filename
1149            self.filename = filename
1150        elif not self.filename:
1151            raise AttributeError("No file name to save to")
1152        SaveDictToProjFile(self.data, self.names, self.filename)
1153
1154    def add_powder_histogram(self, datafile, iparams, phases=[], fmthint=None):
1155        """Loads a powder data histogram into the project.
1156
1157        Automatically checks for an instrument parameter file, or one can be
1158        provided. Note that in unix fashion, "~" can be used to indicate the
1159        home directory (e.g. ~/G2data/data.fxye).
1160
1161        :param str datafile: The powder data file to read, a filename.
1162        :param str iparams: The instrument parameters file, a filename.
1163        :param list phases: Phases to link to the new histogram
1164        :param str fmthint: If specified, only importers where the format name
1165          (reader.formatName, as shown in Import menu) containing the
1166          supplied string will be tried as importers. If not specified, all
1167          importers consistent with the file extension will be tried
1168          (equivalent to "guess format" in menu).
1169
1170        :returns: A :class:`G2PwdrData` object representing
1171            the histogram
1172        """
1173        LoadG2fil()
1174        datafile = os.path.abspath(os.path.expanduser(datafile))
1175        iparams = os.path.abspath(os.path.expanduser(iparams))
1176        pwdrreaders = import_generic(datafile, PwdrDataReaders,fmthint=fmthint)
1177        histname, new_names, pwdrdata = load_pwd_from_reader(
1178                                          pwdrreaders[0], iparams,
1179                                          [h.name for h in self.histograms()])
1180        if histname in self.data:
1181            print("Warning - redefining histogram", histname)
1182        elif self.names[-1][0] == 'Phases':
1183            self.names.insert(-1, new_names)
1184        else:
1185            self.names.append(new_names)
1186        self.data[histname] = pwdrdata
1187        self.update_ids()
1188
1189        for phase in phases:
1190            phase = self.phase(phase)
1191            self.link_histogram_phase(histname, phase)
1192
1193        return self.histogram(histname)
1194
1195    def add_phase(self, phasefile, phasename=None, histograms=[], fmthint=None):
1196        """Loads a phase into the project from a .cif file
1197
1198        :param str phasefile: The CIF file from which to import the phase.
1199        :param str phasename: The name of the new phase, or None for the default
1200        :param list histograms: The names of the histograms to associate with
1201            this phase
1202        :param str fmthint: If specified, only importers where the format name
1203          (reader.formatName, as shown in Import menu) containing the
1204          supplied string will be tried as importers. If not specified, all
1205          importers consistent with the file extension will be tried
1206          (equivalent to "guess format" in menu).
1207
1208        :returns: A :class:`G2Phase` object representing the
1209            new phase.
1210        """
1211        LoadG2fil()
1212        histograms = [self.histogram(h).name for h in histograms]
1213        phasefile = os.path.abspath(os.path.expanduser(phasefile))
1214
1215        # TODO handle multiple phases in a file
1216        phasereaders = import_generic(phasefile, PhaseReaders, fmthint=fmthint)
1217        phasereader = phasereaders[0]
1218       
1219        phasename = phasename or phasereader.Phase['General']['Name']
1220        phaseNameList = [p.name for p in self.phases()]
1221        phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
1222        phasereader.Phase['General']['Name'] = phasename
1223
1224        if 'Phases' not in self.data:
1225            self.data[u'Phases'] = { 'data': None }
1226        assert phasename not in self.data['Phases'], "phase names should be unique"
1227        self.data['Phases'][phasename] = phasereader.Phase
1228
1229        if phasereader.Constraints:
1230            Constraints = self.data['Constraints']
1231            for i in phasereader.Constraints:
1232                if isinstance(i, dict):
1233                    if '_Explain' not in Constraints:
1234                        Constraints['_Explain'] = {}
1235                    Constraints['_Explain'].update(i)
1236                else:
1237                    Constraints['Phase'].append(i)
1238
1239        data = self.data['Phases'][phasename]
1240        generalData = data['General']
1241        SGData = generalData['SGData']
1242        NShkl = len(G2spc.MustrainNames(SGData))
1243        NDij = len(G2spc.HStrainNames(SGData))
1244        Super = generalData.get('Super', 0)
1245        if Super:
1246            SuperVec = np.array(generalData['SuperVec'][0])
1247        else:
1248            SuperVec = []
1249        UseList = data['Histograms']
1250
1251        for hist in histograms:
1252            self.link_histogram_phase(hist, phasename)
1253
1254        for obj in self.names:
1255            if obj[0] == 'Phases':
1256                phasenames = obj
1257                break
1258        else:
1259            phasenames = [u'Phases']
1260            self.names.append(phasenames)
1261        phasenames.append(phasename)
1262
1263        # TODO should it be self.filename, not phasefile?
1264        SetupGeneral(data, os.path.dirname(phasefile))
1265        self.index_ids()
1266
1267        self.update_ids()
1268        return self.phase(phasename)
1269
1270    def link_histogram_phase(self, histogram, phase):
1271        """Associates a given histogram and phase.
1272
1273        .. seealso::
1274
1275            :meth:`G2Project.histogram`
1276            :meth:`G2Project.phase`"""
1277        hist = self.histogram(histogram)
1278        phase = self.phase(phase)
1279
1280        generalData = phase['General']
1281
1282        if hist.name.startswith('HKLF '):
1283            raise NotImplementedError("HKLF not yet supported")
1284        elif hist.name.startswith('PWDR '):
1285            hist['Reflection Lists'][generalData['Name']] = {}
1286            UseList = phase['Histograms']
1287            SGData = generalData['SGData']
1288            NShkl = len(G2spc.MustrainNames(SGData))
1289            NDij = len(G2spc.HStrainNames(SGData))
1290            UseList[hist.name] = SetDefaultDData('PWDR', hist.name, NShkl=NShkl, NDij=NDij)
1291            UseList[hist.name]['hId'] = hist.id
1292            for key, val in [('Use', True), ('LeBail', False),
1293                             ('newLeBail', True),
1294                             ('Babinet', {'BabA': [0.0, False],
1295                                          'BabU': [0.0, False]})]:
1296                if key not in UseList[hist.name]:
1297                    UseList[hist.name][key] = val
1298        else:
1299            raise RuntimeError("Unexpected histogram" + hist.name)
1300
1301
1302    def reload(self):
1303        """Reload self from self.filename"""
1304        data, names = LoadDictFromProjFile(self.filename)
1305        self.names = names
1306        # Need to deep copy the new data file data into the current tree,
1307        # so that any existing G2Phase, or G2PwdrData objects will still be
1308        # valid
1309        _deep_copy_into(from_=data, into=self.data)
1310
1311    def refine(self, newfile=None, printFile=None, makeBack=False):
1312        # TODO migrate to RefineCore
1313        # G2strMain.RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
1314        #      calcControls,pawleyLookup,ifPrint,printFile,dlg)
1315        # index_ids will automatically save the project
1316        self.index_ids()
1317        # TODO G2strMain does not properly use printFile
1318        G2strMain.Refine(self.filename, makeBack=makeBack)
1319        # Reload yourself
1320        self.reload()
1321
1322    def histogram(self, histname):
1323        """Returns the histogram named histname, or None if it does not exist.
1324
1325        :param histname: The name of the histogram (str), or ranId or index.
1326        :returns: A :class:`G2PwdrData` object, or None if
1327            the histogram does not exist
1328
1329        .. seealso::
1330            :meth:`G2Project.histograms`
1331            :meth:`G2Project.phase`
1332            :meth:`G2Project.phases`
1333            """
1334        if isinstance(histname, G2PwdrData):
1335            if histname.proj == self:
1336                return histname
1337        if histname in self.data:
1338            return G2PwdrData(self.data[histname], self)
1339        try:
1340            # see if histname is an id or ranId
1341            histname = int(histname)
1342        except ValueError:
1343            return
1344
1345        for histogram in self.histograms():
1346            if histogram.id == histname or histogram.ranId == histname:
1347                return histogram
1348
1349    def histograms(self):
1350        """Return a list of all histograms, as
1351        :class:`G2PwdrData` objects
1352
1353        .. seealso::
1354            :meth:`G2Project.histograms`
1355            :meth:`G2Project.phase`
1356            :meth:`G2Project.phases`
1357            """
1358        output = []
1359        for obj in self.names:
1360            if len(obj) > 1 and obj[0] != u'Phases':
1361                output.append(self.histogram(obj[0]))
1362        return output
1363
1364    def phase(self, phasename):
1365        """
1366        Gives an object representing the specified phase in this project.
1367
1368        :param str phasename: The name of the desired phase. Either the name
1369            (str), the phase's ranId, or the phase's index
1370        :returns: A :class:`G2Phase` object
1371        :raises: KeyError
1372
1373        .. seealso::
1374            :meth:`G2Project.histograms`
1375            :meth:`G2Project.phase`
1376            :meth:`G2Project.phases`
1377            """
1378        phases = self.data['Phases']
1379        if phasename in phases:
1380            return G2Phase(phases[phasename], phasename, self)
1381
1382        try:
1383            # phasename should be phase index or ranId
1384            phasename = int(phasename)
1385        except ValueError:
1386            return
1387
1388        for phase in self.phases():
1389            if phase.id == phasename or phase.ranId == phasename:
1390                return phase
1391
1392    def phases(self):
1393        """
1394        Returns a list of all the phases in the project.
1395
1396        :returns: A list of :class:`G2Phase` objects
1397
1398        .. seealso::
1399            :meth:`G2Project.histogram`
1400            :meth:`G2Project.histograms`
1401            :meth:`G2Project.phase`
1402            """
1403        for obj in self.names:
1404            if obj[0] == 'Phases':
1405                return [self.phase(p) for p in obj[1:]]
1406        return []
1407
1408    def update_ids(self):
1409        """Makes sure all phases and histograms have proper hId and pId"""
1410        # Translated from GetUsedHistogramsAndPhasesfromTree,
1411        #   GSASIIdataGUI.py:4107
1412        for i, h in enumerate(self.histograms()):
1413            h.id = i
1414        for i, p in enumerate(self.phases()):
1415            p.id = i
1416
1417    def do_refinements(self, refinements, histogram='all', phase='all',
1418                       outputnames=None, makeBack=False):
1419        """Conducts one or a series of refinements according to the
1420           input provided in parameter refinements. This is a wrapper
1421           around :meth:`iter_refinements`
1422
1423        :param list refinements: A list of dictionaries specifiying changes to be made to
1424            parameters before refinements are conducted.
1425            See the :ref:`Refinement_recipe` section for how this is defined.
1426        :param str histogram: Name of histogram for refinements to be applied
1427            to, or 'all'; note that this can be overridden for each refinement
1428            step via a "histograms" entry in the dict.
1429        :param str phase: Name of phase for refinements to be applied to, or
1430            'all'; note that this can be overridden for each refinement
1431            step via a "phases" entry in the dict.
1432        :param list outputnames: Provides a list of project (.gpx) file names
1433            to use for each refinement step (specifying None skips the save step).
1434            See :meth:`save`.
1435            Note that this can be overridden using an "output" entry in the dict.
1436        :param bool makeBack: determines if a backup ).bckX.gpx) file is made
1437            before a refinement is performed. The default is False.
1438           
1439        To perform a single refinement without changing any parameters, use this
1440        call:
1441
1442        .. code-block::  python
1443       
1444            my_project.do_refinements([])
1445        """
1446       
1447        for proj in self.iter_refinements(refinements, histogram, phase,
1448                                          outputnames, makeBack):
1449            pass
1450        return self
1451
1452    def iter_refinements(self, refinements, histogram='all', phase='all',
1453                         outputnames=None, makeBack=False):
1454        """Conducts a series of refinements, iteratively. Stops after every
1455        refinement and yields this project, to allow error checking or
1456        logging of intermediate results. Parameter use is the same as for
1457        :meth:`do_refinements` (which calls this method).
1458
1459        >>> def checked_refinements(proj):
1460        ...     for p in proj.iter_refinements(refs):
1461        ...         # Track intermediate results
1462        ...         log(p.histogram('0').residuals)
1463        ...         log(p.phase('0').get_cell())
1464        ...         # Check if parameter diverged, nonsense answer, or whatever
1465        ...         if is_something_wrong(p):
1466        ...             raise Exception("I need a human!")
1467
1468           
1469        """
1470        if outputnames:
1471            if len(refinements) != len(outputnames):
1472                raise ValueError("Should have same number of outputs to"
1473                                 "refinements")
1474        else:
1475            outputnames = [None for r in refinements]
1476
1477        for output, refinedict in zip(outputnames, refinements):
1478            if 'histograms' in refinedict:
1479                hist = refinedict['histograms']
1480            else:
1481                hist = histogram
1482            if 'phases' in refinedict:
1483                ph = refinedict['phases']
1484            else:
1485                ph = phase
1486            if 'output' in refinedict:
1487                output = refinedict['output']
1488            self.set_refinement(refinedict, hist, ph)
1489            # Handle 'once' args - refinements that are disabled after this
1490            # refinement
1491            if 'once' in refinedict:
1492                temp = {'set': refinedict['once']}
1493                self.set_refinement(temp, hist, ph)
1494
1495            if output:
1496                self.save(output)
1497
1498            if 'skip' not in refinedict:
1499                self.refine(makeBack=makeBack)
1500            yield self
1501
1502            # Handle 'once' args - refinements that are disabled after this
1503            # refinement
1504            if 'once' in refinedict:
1505                temp = {'clear': refinedict['once']}
1506                self.set_refinement(temp, hist, ph)
1507            if 'call' in refinedict:
1508                fxn = refinedict['call']
1509                if callable(fxn):
1510                    fxn(*refinedict.get('callargs',[self]))
1511                elif callable(eval(fxn)):
1512                    eval(fxn)(*refinedict.get('callargs',[self]))
1513                else:
1514                    raise G2ScriptException("Error: call value {} is not callable".format(fxn))
1515
1516    def set_refinement(self, refinement, histogram='all', phase='all'):
1517        """Apply specified refinements to a given histogram(s) or phase(s).
1518
1519        :param dict refinement: The refinements to be conducted
1520        :param histogram: Specifies either 'all' (default), a single histogram or
1521          a list of histograms. Histograms may be specified as histogram objects
1522          (see :class:`G2PwdrData`), the histogram name (str) or the index number (int)
1523          of the histogram in the project, numbered starting from 0.
1524          Omitting the parameter or the string 'all' indicates that parameters in
1525          all histograms should be set.
1526        :param phase: Specifies either 'all' (default), a single phase or
1527          a list of phases. Phases may be specified as phase objects
1528          (see :class:`G2Phase`), the phase name (str) or the index number (int)
1529          of the phase in the project, numbered starting from 0.
1530          Omitting the parameter or the string 'all' indicates that parameters in
1531          all phases should be set.
1532
1533        Note that refinement parameters are categorized as one of three types:
1534
1535        1. Histogram parameters
1536        2. Phase parameters
1537        3. Histogram-and-Phase (HAP) parameters
1538       
1539        .. seealso::
1540            :meth:`G2PwdrData.set_refinements`
1541            :meth:`G2PwdrData.clear_refinements`
1542            :meth:`G2Phase.set_refinements`
1543            :meth:`G2Phase.clear_refinements`
1544            :meth:`G2Phase.set_HAP_refinements`
1545            :meth:`G2Phase.clear_HAP_refinements`"""
1546
1547        if histogram == 'all':
1548            hists = self.histograms()
1549        elif isinstance(histogram, list) or isinstance(histogram, tuple):
1550            hists = []
1551            for h in histogram:
1552                if isinstance(h, str) or isinstance(h, int):
1553                    hists.append(self.histogram(h))
1554                else:
1555                    hists.append(h)
1556        elif isinstance(histogram, str) or isinstance(histogram, int):
1557            hists = [self.histogram(histogram)]
1558        else:
1559            hists = [histogram]
1560
1561        if phase == 'all':
1562            phases = self.phases()
1563        elif isinstance(phase, list) or isinstance(phase, tuple):
1564            phases = []
1565            for ph in phase:
1566                if isinstance(ph, str) or isinstance(ph, int):
1567                    phases.append(self.phase(ph))
1568                else:
1569                    phases.append(ph)
1570        elif isinstance(phase, str) or isinstance(phase, int):
1571            phases = [self.phase(phase)]
1572        else:
1573            phases = [phase]
1574
1575        # TODO: HAP parameters:
1576        #   Babinet
1577        #   Extinction
1578        #   HStrain
1579        #   Mustrain
1580        #   Pref. Ori
1581        #   Size
1582
1583        pwdr_set = {}
1584        phase_set = {}
1585        hap_set = {}
1586        for key, val in refinement.get('set', {}).items():
1587            # Apply refinement options
1588            if G2PwdrData.is_valid_refinement_key(key):
1589                pwdr_set[key] = val
1590            elif G2Phase.is_valid_refinement_key(key):
1591                phase_set[key] = val
1592            elif G2Phase.is_valid_HAP_refinement_key(key):
1593                hap_set[key] = val
1594            else:
1595                raise ValueError("Unknown refinement key", key)
1596
1597        for hist in hists:
1598            hist.set_refinements(pwdr_set)
1599        for phase in phases:
1600            phase.set_refinements(phase_set)
1601        for phase in phases:
1602            phase.set_HAP_refinements(hap_set, hists)
1603
1604        pwdr_clear = {}
1605        phase_clear = {}
1606        hap_clear = {}
1607        for key, val in refinement.get('clear', {}).items():
1608            # Clear refinement options
1609            if G2PwdrData.is_valid_refinement_key(key):
1610                pwdr_clear[key] = val
1611            elif G2Phase.is_valid_refinement_key(key):
1612                phase_clear[key] = val
1613            elif G2Phase.is_valid_HAP_refinement_key(key):
1614                hap_set[key] = val
1615            else:
1616                raise ValueError("Unknown refinement key", key)
1617
1618        for hist in hists:
1619            hist.clear_refinements(pwdr_clear)
1620        for phase in phases:
1621            phase.clear_refinements(phase_clear)
1622        for phase in phases:
1623            phase.clear_HAP_refinements(hap_clear, hists)
1624
1625    def index_ids(self):
1626        import GSASIIstrIO as G2strIO
1627        self.save()
1628        return G2strIO.GetUsedHistogramsAndPhases(self.filename)
1629
1630    def add_constraint_raw(self, cons_scope, constr):
1631        """Adds a constraint of type consType to the project.
1632        cons_scope should be one of "Hist", "Phase", "HAP", or "Global".
1633
1634        WARNING it does not check the constraint is well-constructed"""
1635        constrs = self.data['Constraints']['data']
1636        if 'Global' not in constrs:
1637            constrs['Global'] = []
1638        constrs[cons_scope].append(constr)
1639
1640    def hold_many(self, vars, type):
1641        """Apply holds for all the variables in vars, for constraint of a given type.
1642
1643        type is passed directly to add_constraint_raw as consType
1644
1645        :param list vars: A list of variables to hold. Either :class:`GSASIIobj.G2VarObj` objects,
1646            string variable specifiers, or arguments for :meth:`make_var_obj`
1647        :param str type: A string constraint type specifier. See
1648            :class:`G2Project.add_constraint_raw`
1649
1650        """
1651        for var in vars:
1652            if isinstance(var, str):
1653                var = self.make_var_obj(var)
1654            elif not isinstance(var, G2obj.G2VarObj):
1655                var = self.make_var_obj(*var)
1656            self.add_constraint_raw(type, [[1.0, var], None, None, 'h'])
1657
1658    def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
1659                     reloadIdx=True):
1660        """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a")
1661        or individual names of phase, histogram, varname, and atomId.
1662
1663        Automatically converts string phase, hist, or atom names into the ID required
1664        by G2VarObj."""
1665
1666        if reloadIdx:
1667            self.index_ids()
1668
1669        # If string representation, short circuit
1670        if hist is None and varname is None and atomId is None:
1671            if isinstance(phase, str) and ':' in phase:
1672                return G2obj.G2VarObj(phase)
1673
1674        # Get phase index
1675        phaseObj = None
1676        if isinstance(phase, G2Phase):
1677            phaseObj = phase
1678            phase = G2obj.PhaseRanIdLookup[phase.ranId]
1679        elif phase in self.data['Phases']:
1680            phaseObj = self.phase(phase)
1681            phaseRanId = phaseObj.ranId
1682            phase = G2obj.PhaseRanIdLookup[phaseRanId]
1683
1684        # Get histogram index
1685        if isinstance(hist, G2PwdrData):
1686            hist = G2obj.HistRanIdLookup[hist.ranId]
1687        elif hist in self.data:
1688            histRanId = self.histogram(hist).ranId
1689            hist = G2obj.HistRanIdLookup[histRanId]
1690
1691        # Get atom index (if any)
1692        if isinstance(atomId, G2AtomRecord):
1693            atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
1694        elif phaseObj:
1695            atomObj = phaseObj.atom(atomId)
1696            if atomObj:
1697                atomRanId = atomObj.ranId
1698                atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
1699
1700        return G2obj.G2VarObj(phase, hist, varname, atomId)
1701
1702
1703class G2AtomRecord(G2ObjectWrapper):
1704    """Wrapper for an atom record. Has convenient accessors via @property.
1705
1706
1707    Available accessors: label, type, refinement_flags, coordinates,
1708        occupancy, ranId, id, adp_flag, uiso
1709
1710    Example:
1711
1712    >>> atom = some_phase.atom("O3")
1713    >>> # We can access the underlying data format
1714    >>> atom.data
1715    ['O3', 'O-2', '', ... ]
1716    >>> # We can also use wrapper accessors
1717    >>> atom.coordinates
1718    (0.33, 0.15, 0.5)
1719    >>> atom.refinement_flags
1720    u'FX'
1721    >>> atom.ranId
1722    4615973324315876477
1723    >>> atom.occupancy
1724    1.0
1725    """
1726    def __init__(self, data, indices, proj):
1727        self.data = data
1728        self.cx, self.ct, self.cs, self.cia = indices
1729        self.proj = proj
1730
1731    @property
1732    def label(self):
1733        return self.data[self.ct-1]
1734
1735    @property
1736    def type(self):
1737        return self.data[self.ct]
1738
1739    @property
1740    def refinement_flags(self):
1741        return self.data[self.ct+1]
1742
1743    @refinement_flags.setter
1744    def refinement_flags(self, other):
1745        # Automatically check it is a valid refinement
1746        for c in other:
1747            if c not in ' FXU':
1748                raise ValueError("Invalid atom refinement: ", other)
1749        self.data[self.ct+1] = other
1750
1751    @property
1752    def coordinates(self):
1753        return tuple(self.data[self.cx:self.cx+3])
1754
1755    @property
1756    def occupancy(self):
1757        return self.data[self.cx+3]
1758
1759    @occupancy.setter
1760    def occupancy(self, val):
1761        self.data[self.cx+3] = float(val)
1762
1763    @property
1764    def ranId(self):
1765        return self.data[self.cia+8]
1766
1767    @property
1768    def adp_flag(self):
1769        # Either 'I' or 'A'
1770        return self.data[self.cia]
1771
1772    @property
1773    def uiso(self):
1774        if self.adp_flag == 'I':
1775            return self.data[self.cia+1]
1776        else:
1777            return self.data[self.cia+2:self.cia+8]
1778
1779    @uiso.setter
1780    def uiso(self, value):
1781        if self.adp_flag == 'I':
1782            self.data[self.cia+1] = float(value)
1783        else:
1784            assert len(value) == 6
1785            self.data[self.cia+2:self.cia+8] = [float(v) for v in value]
1786
1787
1788class G2PwdrData(G2ObjectWrapper):
1789    """Wraps a Powder Data Histogram."""
1790    def __init__(self, data, proj):
1791        self.data = data
1792        self.proj = proj
1793
1794    @staticmethod
1795    def is_valid_refinement_key(key):
1796        valid_keys = ['Limits', 'Sample Parameters', 'Background',
1797                      'Instrument Parameters']
1798        return key in valid_keys
1799
1800    @property
1801    def name(self):
1802        return self.data['data'][-1]
1803
1804    @property
1805    def ranId(self):
1806        return self.data['data'][0]['ranId']
1807
1808    @property
1809    def residuals(self):
1810        data = self.data['data'][0]
1811        return {key: data[key]
1812                for key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']}
1813
1814    @property
1815    def id(self):
1816        self.proj.update_ids()
1817        return self.data['data'][0]['hId']
1818
1819    @id.setter
1820    def id(self, val):
1821        self.data['data'][0]['hId'] = val
1822
1823    def fit_fixed_points(self):
1824        """Attempts to apply a background fit to the fixed points currently specified."""
1825        def SetInstParms(Inst):
1826            dataType = Inst['Type'][0]
1827            insVary = []
1828            insNames = []
1829            insVals = []
1830            for parm in Inst:
1831                insNames.append(parm)
1832                insVals.append(Inst[parm][1])
1833                if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha',
1834                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1835                        Inst[parm][2] = False
1836            instDict = dict(zip(insNames, insVals))
1837            instDict['X'] = max(instDict['X'], 0.01)
1838            instDict['Y'] = max(instDict['Y'], 0.01)
1839            if 'SH/L' in instDict:
1840                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
1841            return dataType, instDict, insVary
1842
1843        bgrnd = self.data['Background']
1844
1845        # Need our fixed points in order
1846        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
1847        X = [x for x, y in bgrnd[1]['FixedPoints']]
1848        Y = [y for x, y in bgrnd[1]['FixedPoints']]
1849
1850        limits = self.data['Limits'][1]
1851        if X[0] > limits[0]:
1852            X = [limits[0]] + X
1853            Y = [Y[0]] + Y
1854        if X[-1] < limits[1]:
1855            X += [limits[1]]
1856            Y += [Y[-1]]
1857
1858        # Some simple lookups
1859        controls = self.proj['Controls']['data']
1860        inst, inst2 = self.data['Instrument Parameters']
1861        pwddata = self.data['data'][1]
1862
1863        # Construct the data for background fitting
1864        xBeg = np.searchsorted(pwddata[0], limits[0])
1865        xFin = np.searchsorted(pwddata[0], limits[1])
1866        xdata = pwddata[0][xBeg:xFin]
1867        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
1868
1869        W = [1]*len(xdata)
1870        Z = [0]*len(xdata)
1871
1872        dataType, insDict, insVary = SetInstParms(inst)
1873        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1874
1875        # Do the fit
1876        data = np.array([xdata, ydata, W, Z, Z, Z])
1877        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
1878                        prevVaryList=bakVary, controls=controls)
1879
1880        # Post-fit
1881        parmDict = {}
1882        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1883        parmDict.update(bakDict)
1884        parmDict.update(insDict)
1885        pwddata[3][xBeg:xFin] *= 0
1886        pwddata[5][xBeg:xFin] *= 0
1887        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
1888
1889        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
1890        # TODO update background
1891        self.proj.save()
1892
1893    def y_calc(self):
1894        return self.data['data'][1][3]
1895
1896    def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
1897        try:
1898            import matplotlib.pyplot as plt
1899            data = self.data['data'][1]
1900            if Yobs:
1901                plt.plot(data[0], data[1], label='Yobs')
1902            if Ycalc:
1903                plt.plot(data[0], data[3], label='Ycalc')
1904            if Background:
1905                plt.plot(data[0], data[4], label='Background')
1906            if Residual:
1907                plt.plot(data[0], data[5], label="Residual")
1908        except ImportError:
1909            pass
1910
1911    def get_wR(self):
1912        """returns the overall weighted profile R factor for a histogram
1913       
1914        :returns: a wR value as a percentage or None if not defined
1915        """
1916        return self['data'][0].get('wR')
1917
1918    def set_refinements(self, refs):
1919        """Sets the refinement parameter 'key' to the specification 'value'
1920
1921        :param dict refs: A dictionary of the parameters to be set. See
1922                          :ref:`Histogram_parameters_table` for a description of
1923                          what these dictionaries should be.
1924
1925        :returns: None
1926
1927        """
1928        do_fit_fixed_points = False
1929        for key, value in refs.items():
1930            if key == 'Limits':
1931                old_limits = self.data['Limits'][1]
1932                new_limits = value
1933                if isinstance(new_limits, dict):
1934                    if 'low' in new_limits:
1935                        old_limits[0] = new_limits['low']
1936                    if 'high' in new_limits:
1937                        old_limits[1] = new_limits['high']
1938                else:
1939                    old_limits[0], old_limits[1] = new_limits
1940            elif key == 'Sample Parameters':
1941                sample = self.data['Sample Parameters']
1942                for sparam in value:
1943                    if sparam not in sample:
1944                        raise ValueError("Unknown refinement parameter, "
1945                                         + str(sparam))
1946                    sample[sparam][1] = True
1947            elif key == 'Background':
1948                bkg, peaks = self.data['Background']
1949
1950                # If True or False, just set the refine parameter
1951                if value in (True, False):
1952                    bkg[1] = value
1953                    return
1954
1955                if 'type' in value:
1956                    bkg[0] = value['type']
1957                if 'refine' in value:
1958                    bkg[1] = value['refine']
1959                if 'no. coeffs' in value:
1960                    cur_coeffs = bkg[2]
1961                    n_coeffs = value['no. coeffs']
1962                    if n_coeffs > cur_coeffs:
1963                        for x in range(n_coeffs - cur_coeffs):
1964                            bkg.append(0.0)
1965                    else:
1966                        for _ in range(cur_coeffs - n_coeffs):
1967                            bkg.pop()
1968                    bkg[2] = n_coeffs
1969                if 'coeffs' in value:
1970                    bkg[3:] = value['coeffs']
1971                if 'FixedPoints' in value:
1972                    peaks['FixedPoints'] = [(float(a), float(b))
1973                                            for a, b in value['FixedPoints']]
1974                if value.get('fit fixed points', False):
1975                    do_fit_fixed_points = True
1976
1977            elif key == 'Instrument Parameters':
1978                instrument, secondary = self.data['Instrument Parameters']
1979                for iparam in value:
1980                    try:
1981                        instrument[iparam][2] = True
1982                    except IndexError:
1983                        raise ValueError("Invalid key:", iparam)
1984            else:
1985                raise ValueError("Unknown key:", key)
1986        # Fit fixed points after the fact - ensure they are after fixed points
1987        # are added
1988        if do_fit_fixed_points:
1989            # Background won't be fit if refinement flag not set
1990            orig = self.data['Background'][0][1]
1991            self.data['Background'][0][1] = True
1992            self.fit_fixed_points()
1993            # Restore the previous value
1994            self.data['Background'][0][1] = orig
1995
1996    def clear_refinements(self, refs):
1997        """Clears the refinement parameter 'key' and its associated value.
1998
1999        :param dict refs: A dictionary of parameters to clear."""
2000        for key, value in refs.items():
2001            if key == 'Limits':
2002                old_limits, cur_limits = self.data['Limits']
2003                cur_limits[0], cur_limits[1] = old_limits
2004            elif key == 'Sample Parameters':
2005                sample = self.data['Sample Parameters']
2006                for sparam in value:
2007                    sample[sparam][1] = False
2008            elif key == 'Background':
2009                bkg, peaks = self.data['Background']
2010
2011                # If True or False, just set the refine parameter
2012                if value in (True, False):
2013                    bkg[1] = False
2014                    return
2015
2016                bkg[1] = False
2017                if 'FixedPoints' in value:
2018                    if 'FixedPoints' in peaks:
2019                        del peaks['FixedPoints']
2020            elif key == 'Instrument Parameters':
2021                instrument, secondary = self.data['Instrument Parameters']
2022                for iparam in value:
2023                    instrument[iparam][2] = False
2024            else:
2025                raise ValueError("Unknown key:", key)
2026
2027
2028class G2Phase(G2ObjectWrapper):
2029    """A wrapper object around a given phase.
2030
2031    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
2032    """
2033    def __init__(self, data, name, proj):
2034        self.data = data
2035        self.name = name
2036        self.proj = proj
2037
2038    @staticmethod
2039    def is_valid_refinement_key(key):
2040        valid_keys = ["Cell", "Atoms", "LeBail"]
2041        return key in valid_keys
2042
2043    @staticmethod
2044    def is_valid_HAP_refinement_key(key):
2045        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
2046                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
2047        return key in valid_keys
2048
2049    def atom(self, atomlabel):
2050        """Returns the atom specified by atomlabel, or None if it does not
2051        exist.
2052
2053        :param str atomlabel: The name of the atom (e.g. "O2")
2054        :returns: A :class:`G2AtomRecord` object
2055            representing the atom.
2056        """
2057        # Consult GSASIIobj.py for the meaning of this
2058        cx, ct, cs, cia = self.data['General']['AtomPtrs']
2059        ptrs = [cx, ct, cs, cia]
2060        atoms = self.data['Atoms']
2061        for atom in atoms:
2062            if atom[ct-1] == atomlabel:
2063                return G2AtomRecord(atom, ptrs, self.proj)
2064
2065    def atoms(self):
2066        """Returns a list of atoms present in the phase.
2067
2068        :returns: A list of :class:`G2AtomRecord` objects.
2069
2070        .. seealso::
2071            :meth:`G2Phase.atom`
2072            :class:`G2AtomRecord`
2073        """
2074        ptrs = self.data['General']['AtomPtrs']
2075        output = []
2076        atoms = self.data['Atoms']
2077        for atom in atoms:
2078            output.append(G2AtomRecord(atom, ptrs, self.proj))
2079        return output
2080
2081    def histograms(self):
2082        output = []
2083        for hname in self.data.get('Histograms', {}).keys():
2084            output.append(self.proj.histogram(hname))
2085        return output
2086
2087    @property
2088    def ranId(self):
2089        return self.data['ranId']
2090
2091    @property
2092    def id(self):
2093        return self.data['pId']
2094
2095    @id.setter
2096    def id(self, val):
2097        self.data['pId'] = val
2098
2099    def get_cell(self):
2100        """Returns a dictionary of the cell parameters, with keys:
2101            'length_a', 'length_b', 'length_c', 'angle_alpha', 'angle_beta', 'angle_gamma', 'volume'
2102
2103        :returns: a dict
2104
2105        .. seealso::
2106           :meth:`G2Phase.get_cell_and_esd`
2107
2108        """
2109        cell = self.data['General']['Cell']
2110        return {'length_a': cell[1], 'length_b': cell[2], 'length_c': cell[3],
2111                'angle_alpha': cell[4], 'angle_beta': cell[5], 'angle_gamma': cell[6],
2112                'volume': cell[7]}
2113
2114    def get_cell_and_esd(self):
2115        """
2116        Returns a pair of dictionaries, the first representing the unit cell, the second
2117        representing the estimated standard deviations of the unit cell.
2118
2119        :returns: a tuple of two dictionaries
2120
2121        .. seealso::
2122           :meth:`G2Phase.get_cell`
2123
2124        """
2125        # translated from GSASIIstrIO.ExportBaseclass.GetCell
2126        import GSASIIstrIO as G2stIO
2127        import GSASIIlattice as G2lat
2128        import GSASIImapvars as G2mv
2129        try:
2130            pfx = str(self.id) + '::'
2131            sgdata = self['General']['SGData']
2132            covDict = self.proj['Covariance']['data']
2133
2134            parmDict = dict(zip(covDict.get('varyList',[]),
2135                                covDict.get('variables',[])))
2136            sigDict = dict(zip(covDict.get('varyList',[]),
2137                               covDict.get('sig',[])))
2138
2139            if covDict.get('covMatrix') is not None:
2140                sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
2141                                                  covDict['varyList'],
2142                                                  parmDict))
2143
2144            A, sigA = G2stIO.cellFill(pfx, sgdata, parmDict, sigDict)
2145            cellSig = G2stIO.getCellEsd(pfx, sgdata, A, self.proj['Covariance']['data'])
2146            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
2147            cellDict, cellSigDict = {}, {}
2148            for i, key in enumerate(['length_a', 'length_b', 'length_c',
2149                                     'angle_alpha', 'angle_beta', 'angle_gamma',
2150                                     'volume']):
2151                cellDict[key] = cellList[i]
2152                cellSigDict[key] = cellSig[i]
2153            return cellDict, cellSigDict
2154        except KeyError:
2155            cell = self.get_cell()
2156            return cell, {key: 0.0 for key in cell}
2157
2158    def export_CIF(self, outputname, quickmode=True):
2159        """Write this phase to a .cif file named outputname
2160
2161        :param str outputname: The name of the .cif file to write to
2162        :param bool quickmode: Currently ignored. Carryover from exports.G2export_CIF"""
2163        # This code is all taken from exports/G2export_CIF.py
2164        # Functions copied have the same names
2165        import GSASIImath as G2mth
2166        import GSASIImapvars as G2mv
2167        from exports import G2export_CIF as cif
2168
2169        CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
2170        CIFname = os.path.splitext(self.proj.filename)[0]
2171        CIFname = os.path.split(CIFname)[1]
2172        CIFname = ''.join([c if ord(c) < 128 else ''
2173                           for c in CIFname.replace(' ', '_')])
2174        try:
2175            author = self.proj['Controls']['data'].get('Author','').strip()
2176        except KeyError:
2177            pass
2178        oneblock = True
2179
2180        covDict = self.proj['Covariance']['data']
2181        parmDict = dict(zip(covDict.get('varyList',[]),
2182                            covDict.get('variables',[])))
2183        sigDict = dict(zip(covDict.get('varyList',[]),
2184                           covDict.get('sig',[])))
2185
2186        if covDict.get('covMatrix') is not None:
2187            sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
2188                                              covDict['varyList'],
2189                                              parmDict))
2190
2191        with open(outputname, 'w') as fp:
2192            fp.write(' \n' + 70*'#' + '\n')
2193            cif.WriteCIFitem(fp, 'data_' + CIFname)
2194            # from exports.G2export_CIF.WritePhaseInfo
2195            cif.WriteCIFitem(fp, '\n# phase info for '+str(self.name) + ' follows')
2196            cif.WriteCIFitem(fp, '_pd_phase_name', self.name)
2197            # TODO get esds
2198            cellDict = self.get_cell()
2199            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
2200            names = ['length_a','length_b','length_c',
2201                     'angle_alpha','angle_beta ','angle_gamma',
2202                     'volume']
2203            for key, val in cellDict.items():
2204                cif.WriteCIFitem(fp, '_cell_' + key, G2mth.ValEsd(val))
2205
2206            cif.WriteCIFitem(fp, '_symmetry_cell_setting',
2207                         self.data['General']['SGData']['SGSys'])
2208
2209            spacegroup = self.data['General']['SGData']['SpGrp'].strip()
2210            # regularize capitalization and remove trailing H/R
2211            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
2212            cif.WriteCIFitem(fp, '_symmetry_space_group_name_H-M', spacegroup)
2213
2214            # generate symmetry operations including centering and center of symmetry
2215            SymOpList, offsetList, symOpList, G2oprList, G2opcodes = G2spc.AllOps(
2216                self.data['General']['SGData'])
2217            cif.WriteCIFitem(fp, 'loop_\n    _space_group_symop_id\n    _space_group_symop_operation_xyz')
2218            for i, op in enumerate(SymOpList,start=1):
2219                cif.WriteCIFitem(fp, '   {:3d}  {:}'.format(i,op.lower()))
2220
2221            # TODO skipped histograms, exports/G2export_CIF.py:880
2222
2223            # report atom params
2224            if self.data['General']['Type'] in ['nuclear','macromolecular']:        #this needs macromolecular variant, etc!
2225                cif.WriteAtomsNuclear(fp, self.data, self.name, parmDict, sigDict, [])
2226                # self._WriteAtomsNuclear(fp, parmDict, sigDict)
2227            else:
2228                raise G2ScriptException("no export for "+str(self.data['General']['Type'])+" coordinates implemented")
2229            # report cell contents
2230            cif.WriteComposition(fp, self.data, self.name, parmDict)
2231            if not quickmode and self.data['General']['Type'] == 'nuclear':      # report distances and angles
2232                # WriteDistances(fp,self.name,SymOpList,offsetList,symOpList,G2oprList)
2233                raise NotImplementedError("only quickmode currently supported")
2234            if 'Map' in self.data['General'] and 'minmax' in self.data['General']['Map']:
2235                cif.WriteCIFitem(fp,'\n# Difference density results')
2236                MinMax = self.data['General']['Map']['minmax']
2237                cif.WriteCIFitem(fp,'_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009))
2238                cif.WriteCIFitem(fp,'_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009))
2239
2240
2241    def set_refinements(self, refs):
2242        """Sets the refinement parameter 'key' to the specification 'value'
2243
2244        :param dict refs: A dictionary of the parameters to be set. See
2245                          :ref:`Phase_parameters_table` for a description of
2246                          this dictionary.
2247
2248        :returns: None"""
2249        for key, value in refs.items():
2250            if key == "Cell":
2251                self.data['General']['Cell'][0] = value
2252
2253            elif key == "Atoms":
2254                for atomlabel, atomrefinement in value.items():
2255                    if atomlabel == 'all':
2256                        for atom in self.atoms():
2257                            atom.refinement_flags = atomrefinement
2258                    else:
2259                        atom = self.atom(atomlabel)
2260                        if atom is None:
2261                            raise ValueError("No such atom: " + atomlabel)
2262                        atom.refinement_flags = atomrefinement
2263
2264            elif key == "LeBail":
2265                hists = self.data['Histograms']
2266                for hname, hoptions in hists.items():
2267                    if 'LeBail' not in hoptions:
2268                        hoptions['newLeBail'] = bool(True)
2269                    hoptions['LeBail'] = bool(value)
2270            else:
2271                raise ValueError("Unknown key:", key)
2272
2273    def clear_refinements(self, refs):
2274        """Clears a given set of parameters.
2275
2276        :param dict refs: The parameters to clear"""
2277        for key, value in refs.items():
2278            if key == "Cell":
2279                self.data['General']['Cell'][0] = False
2280            elif key == "Atoms":
2281                cx, ct, cs, cia = self.data['General']['AtomPtrs']
2282
2283                for atomlabel in value:
2284                    atom = self.atom(atomlabel)
2285                    # Set refinement to none
2286                    atom.refinement_flags = ' '
2287            elif key == "LeBail":
2288                hists = self.data['Histograms']
2289                for hname, hoptions in hists.items():
2290                    if 'LeBail' not in hoptions:
2291                        hoptions['newLeBail'] = True
2292                    hoptions['LeBail'] = False
2293            else:
2294                raise ValueError("Unknown key:", key)
2295
2296    def set_HAP_refinements(self, refs, histograms='all'):
2297        """Sets the given HAP refinement parameters between this phase and
2298        the given histograms
2299
2300        :param dict refs: A dictionary of the parameters to be set. See
2301                          :ref:`HAP_parameters_table` for a description of this
2302                          dictionary.
2303        :param histograms: Either 'all' (default) or a list of the histograms
2304            whose HAP parameters will be set with this phase. Histogram and phase
2305            must already be associated
2306
2307        :returns: None
2308        """
2309        if histograms == 'all':
2310            histograms = self.data['Histograms'].values()
2311        else:
2312            histograms = [self.data['Histograms'][h.name] for h in histograms]
2313
2314        for key, val in refs.items():
2315            for h in histograms:
2316                if key == 'Babinet':
2317                    try:
2318                        sets = list(val)
2319                    except ValueError:
2320                        sets = ['BabA', 'BabU']
2321                    for param in sets:
2322                        if param not in ['BabA', 'BabU']:
2323                            raise ValueError("Not sure what to do with" + param)
2324                        for hist in histograms:
2325                            hist['Babinet'][param][1] = True
2326                elif key == 'Extinction':
2327                    for h in histograms:
2328                        h['Extinction'][1] = bool(val)
2329                elif key == 'HStrain':
2330                    for h in histograms:
2331                        h['HStrain'][1] = [bool(val) for p in h['HStrain'][1]]
2332                elif key == 'Mustrain':
2333                    for h in histograms:
2334                        mustrain = h['Mustrain']
2335                        newType = None
2336                        direction = None
2337                        if isinstance(val, strtypes):
2338                            if val in ['isotropic', 'uniaxial', 'generalized']:
2339                                newType = val
2340                            else:
2341                                raise ValueError("Not a Mustrain type: " + val)
2342                        elif isinstance(val, dict):
2343                            newType = val.get('type', None)
2344                            direction = val.get('direction', None)
2345
2346                        if newType:
2347                            mustrain[0] = newType
2348                            if newType == 'isotropic':
2349                                mustrain[2][0] = True
2350                                mustrain[5] = [False for p in mustrain[4]]
2351                            elif newType == 'uniaxial':
2352                                if 'refine' in val:
2353                                    types = val['refine']
2354                                    if isinstance(types, strtypes):
2355                                        types = [types]
2356                                    elif isinstance(types, bool):
2357                                        mustrain[2][0] = types
2358                                        mustrain[2][1] = types
2359                                        types = []
2360                                    else:
2361                                        raise ValueError("Not sure what to do with: "
2362                                                         + str(types))
2363                                else:
2364                                    types = []
2365
2366                                for unitype in types:
2367                                    if unitype == 'equatorial':
2368                                        mustrain[2][0] = True
2369                                    elif unitype == 'axial':
2370                                        mustrain[2][1] = True
2371                                    else:
2372                                        msg = 'Invalid uniaxial mustrain type'
2373                                        raise ValueError(msg + ': ' + unitype)
2374                            else:  # newtype == 'generalized'
2375                                mustrain[2] = [False for p in mustrain[1]]
2376
2377                        if direction:
2378                            if len(direction) != 3:
2379                                raise ValueError("Expected hkl, found", direction)
2380                            direction = [int(n) for n in direction]
2381                            mustrain[3] = direction
2382                elif key == 'Size':
2383                    for h in histograms:
2384                        size = h['Size']
2385                        newType = None
2386                        direction = None
2387                        if isinstance(val, strtypes):
2388                            if val in ['isotropic', 'uniaxial', 'ellipsoidal']:
2389                                newType = val
2390                            else:
2391                                raise ValueError("Not a valid Size type: " + val)
2392                        elif isinstance(val, dict):
2393                            newType = val.get('type', None)
2394                            direction = val.get('direction', None)
2395
2396                        if newType:
2397                            size[0] = newType
2398                            refine = val.get('refine')
2399                            if newType == 'isotropic' and refine is not None:
2400                                size[2][0] = bool(refine)
2401                            elif newType == 'uniaxial' and refine is not None:
2402                                size[2][1] = bool(refine)
2403                                size[2][2] = bool(refine)
2404                            elif newType == 'ellipsoidal' and refine is not None:
2405                                size[5] = [bool(refine) for p in size[5]]
2406
2407                        if direction:
2408                            if len(direction) != 3:
2409                                raise ValueError("Expected hkl, found", direction)
2410                            direction = [int(n) for n in direction]
2411                            size[3] = direction
2412                elif key == 'Pref.Ori.':
2413                    for h in histograms:
2414                        h['Pref.Ori.'][2] = bool(val)
2415                elif key == 'Show':
2416                    for h in histograms:
2417                        h['Show'] = bool(val)
2418                elif key == 'Use':
2419                    for h in histograms:
2420                        h['Use'] = bool(val)
2421                elif key == 'Scale':
2422                    for h in histograms:
2423                        h['Scale'][1] = bool(val)
2424                else:
2425                    print(u'Unknown HAP key: '+key)
2426
2427    def clear_HAP_refinements(self, refs, histograms='all'):
2428        """Clears the given HAP refinement parameters between this phase and
2429        the given histograms
2430
2431        :param dict refs: A dictionary of the parameters to be cleared.
2432        :param histograms: Either 'all' (default) or a list of the histograms
2433            whose HAP parameters will be cleared with this phase. Histogram and
2434            phase must already be associated
2435
2436        :returns: None
2437        """
2438        if histograms == 'all':
2439            histograms = self.data['Histograms'].values()
2440        else:
2441            histograms = [self.data['Histograms'][h.name] for h in histograms]
2442
2443        for key, val in refs.items():
2444            for h in histograms:
2445                if key == 'Babinet':
2446                    try:
2447                        sets = list(val)
2448                    except ValueError:
2449                        sets = ['BabA', 'BabU']
2450                    for param in sets:
2451                        if param not in ['BabA', 'BabU']:
2452                            raise ValueError("Not sure what to do with" + param)
2453                        for hist in histograms:
2454                            hist['Babinet'][param][1] = False
2455                elif key == 'Extinction':
2456                    for h in histograms:
2457                        h['Extinction'][1] = False
2458                elif key == 'HStrain':
2459                    for h in histograms:
2460                        h['HStrain'][1] = [False for p in h['HStrain'][1]]
2461                elif key == 'Mustrain':
2462                    for h in histograms:
2463                        mustrain = h['Mustrain']
2464                        mustrain[2] = [False for p in mustrain[2]]
2465                        mustrain[5] = [False for p in mustrain[4]]
2466                elif key == 'Pref.Ori.':
2467                    for h in histograms:
2468                        h['Pref.Ori.'][2] = False
2469                elif key == 'Show':
2470                    for h in histograms:
2471                        h['Show'] = False
2472                elif key == 'Size':
2473                    for h in histograms:
2474                        size = h['Size']
2475                        size[2] = [False for p in size[2]]
2476                        size[5] = [False for p in size[5]]
2477                elif key == 'Use':
2478                    for h in histograms:
2479                        h['Use'] = False
2480                elif key == 'Scale':
2481                    for h in histograms:
2482                        h['Scale'][1] = False
2483                else:
2484                    print(u'Unknown HAP key: '+key)
2485
2486
2487##########################
2488# Command Line Interface #
2489##########################
2490# Each of these takes an argparse.Namespace object as their argument,
2491# representing the parsed command-line arguments for the relevant subcommand.
2492# The argument specification for each is in the subcommands dictionary (see
2493# below)
2494
2495commandhelp={}
2496commandhelp["create"] = "creates a GSAS-II project, optionally adding histograms and/or phases"
2497def create(args):
2498    """The create subcommand. This creates a GSAS-II project, optionally adding histograms and/or phases::
2499
2500  usage: GSASIIscriptable.py create [-h] [-d HISTOGRAMS [HISTOGRAMS ...]]
2501                                  [-i IPARAMS [IPARAMS ...]]
2502                                  [-p PHASES [PHASES ...]]
2503                                  filename
2504                                 
2505positional arguments::
2506
2507  filename              the project file to create. should end in .gpx
2508
2509optional arguments::
2510
2511  -h, --help            show this help message and exit
2512  -d HISTOGRAMS [HISTOGRAMS ...], --histograms HISTOGRAMS [HISTOGRAMS ...]
2513                        list of datafiles to add as histograms
2514  -i IPARAMS [IPARAMS ...], --iparams IPARAMS [IPARAMS ...]
2515                        instrument parameter file, must be one for every
2516                        histogram
2517  -p PHASES [PHASES ...], --phases PHASES [PHASES ...]
2518                        list of phases to add. phases are automatically
2519                        associated with all histograms given.
2520
2521    """
2522    proj = G2Project(filename=args.filename)
2523
2524    hist_objs = []
2525    if args.histograms:
2526        for h,i in zip(args.histograms,args.iparams):
2527            print("Adding histogram from",h,"with instparm ",i)
2528            hist_objs.append(proj.add_powder_histogram(h, i))
2529
2530    if args.phases: 
2531        for p in args.phases:
2532            print("Adding phase from",p)
2533            proj.add_phase(p, histograms=hist_objs)
2534        print('Linking phase(s) to histogram(s):')
2535        for h in hist_objs:
2536            print ('   '+h.name)
2537
2538    proj.save()
2539
2540commandhelp["add"] = "adds histograms and/or phases to GSAS-II project"
2541def add(args):
2542    """The add subcommand. This adds histograms and/or phases to GSAS-II project::
2543
2544  usage: GSASIIscriptable.py add [-h] [-d HISTOGRAMS [HISTOGRAMS ...]]
2545                               [-i IPARAMS [IPARAMS ...]]
2546                               [-hf HISTOGRAMFORMAT] [-p PHASES [PHASES ...]]
2547                               [-pf PHASEFORMAT] [-l HISTLIST [HISTLIST ...]]
2548                               filename
2549
2550
2551positional arguments::
2552
2553  filename              the project file to open. Should end in .gpx
2554
2555optional arguments::
2556
2557  -h, --help            show this help message and exit
2558  -d HISTOGRAMS [HISTOGRAMS ...], --histograms HISTOGRAMS [HISTOGRAMS ...]
2559                        list of datafiles to add as histograms
2560  -i IPARAMS [IPARAMS ...], --iparams IPARAMS [IPARAMS ...]
2561                        instrument parameter file, must be one for every
2562                        histogram
2563  -hf HISTOGRAMFORMAT, --histogramformat HISTOGRAMFORMAT
2564                        format hint for histogram import. Applies to all
2565                        histograms
2566  -p PHASES [PHASES ...], --phases PHASES [PHASES ...]
2567                        list of phases to add. phases are automatically
2568                        associated with all histograms given.
2569  -pf PHASEFORMAT, --phaseformat PHASEFORMAT
2570                        format hint for phase import. Applies to all phases.
2571                        Example: -pf CIF
2572  -l HISTLIST [HISTLIST ...], --histlist HISTLIST [HISTLIST ...]
2573                        list of histgram indices to associate with added
2574                        phases. If not specified, phases are associated with
2575                        all previously loaded histograms. Example: -l 2 3 4
2576   
2577    """
2578    proj = G2Project(args.filename)
2579
2580    if args.histograms:
2581        for h,i in zip(args.histograms,args.iparams):
2582            print("Adding histogram from",h,"with instparm ",i)
2583            proj.add_powder_histogram(h, i, fmthint=args.histogramformat)
2584
2585    if args.phases: 
2586        if not args.histlist:
2587            histlist = proj.histograms()
2588        else:
2589            histlist = [proj.histogram(i) for i in args.histlist]
2590
2591        for p in args.phases:
2592            print("Adding phase from",p)
2593            proj.add_phase(p, histograms=histlist, fmthint=args.phaseformat)
2594           
2595        if not args.histlist:
2596            print('Linking phase(s) to all histogram(s)')
2597        else:
2598            print('Linking phase(s) to histogram(s):')
2599            for h in histlist:
2600                print ('   '+h.name)
2601
2602    proj.save()
2603
2604
2605commandhelp["dump"] = "Shows the contents of a GSAS-II project"
2606def dump(args):
2607    """The dump subcommand shows the contents of a GSAS-II project::
2608
2609       usage: GSASIIscriptable.py dump [-h] [-d] [-p] [-r] files [files ...]
2610
2611positional arguments::
2612
2613  files
2614
2615optional arguments::
2616
2617  -h, --help        show this help message and exit
2618  -d, --histograms  list histograms in files, overrides --raw
2619  -p, --phases      list phases in files, overrides --raw
2620  -r, --raw         dump raw file contents, default
2621 
2622    """
2623    if not args.histograms and not args.phases:
2624        args.raw = True
2625    if args.raw:
2626        import IPython.lib.pretty as pretty
2627
2628    for fname in args.files:
2629        if args.raw:
2630            proj, nameList = LoadDictFromProjFile(fname)
2631            print("file:", fname)
2632            print("names:", nameList)
2633            for key, val in proj.items():
2634                print(key, ":")
2635                pretty.pprint(val)
2636        else:
2637            proj = G2Project(fname)
2638            if args.histograms:
2639                hists = proj.histograms()
2640                for h in hists:
2641                    print(fname, "hist", h.id, h.name)
2642            if args.phases:
2643                phase_list = proj.phases()
2644                for p in phase_list:
2645                    print(fname, "phase", p.id, p.name)
2646
2647
2648commandhelp["browse"] = "Load a GSAS-II project and then open a IPython shell to browse it"
2649def IPyBrowse(args):
2650    """Load a .gpx file and then open a IPython shell to browse it::
2651
2652  usage: GSASIIscriptable.py browse [-h] files [files ...]
2653
2654positional arguments::
2655
2656  files       list of files to browse
2657
2658optional arguments::
2659
2660  -h, --help  show this help message and exit
2661
2662    """
2663    for fname in args.files:
2664        proj, nameList = LoadDictFromProjFile(fname)
2665        msg = "\nfname {} loaded into proj (dict) with names in nameList".format(fname)
2666        GSASIIpath.IPyBreak_base(msg)
2667        break
2668
2669
2670commandhelp["refine"] = '''
2671Conducts refinements on GSAS-II projects according to a list of refinement
2672steps in a JSON dict
2673'''
2674def refine(args):
2675    """Conducts refinements on GSAS-II projects according to a JSON refinement dict::
2676   
2677  usage: GSASIIscriptable.py refine [-h] gpxfile [refinements]
2678
2679positional arguments::
2680
2681  gpxfile      the project file to refine
2682  refinements  json file of refinements to apply. if not present refines file
2683               as-is
2684
2685optional arguments::
2686
2687  -h, --help   show this help message and exit
2688 
2689    """
2690    proj = G2Project(args.gpxfile)
2691    if args.refinements is None:
2692        proj.refine()
2693    else:
2694        import json
2695        with open(args.refinements) as refs:
2696            refs = json.load(refs)
2697        if type(refs) is not dict:
2698            raise G2ScriptException("Error: JSON object must be a dict.")
2699        if "code" in refs:
2700            print("executing code:\n|  ",'\n|  '.join(refs['code']))
2701            exec('\n'.join(refs['code']))
2702        proj.do_refinements(refs['refinements'])
2703
2704
2705commandhelp["seqrefine"] = "Not implemented. Placeholder for eventual sequential refinement implementation"
2706def seqrefine(args):
2707    """The seqrefine subcommand"""
2708    raise NotImplementedError("seqrefine is not yet implemented")
2709
2710
2711commandhelp["export"] = "Export phase as CIF"
2712def export(args):
2713    """The export subcommand: Exports phase as CIF::
2714
2715      usage: GSASIIscriptable.py export [-h] gpxfile phase exportfile
2716
2717positional arguments::
2718
2719  gpxfile     the project file from which to export
2720  phase       identifier of phase to export
2721  exportfile  the .cif file to export to
2722
2723optional arguments::
2724
2725  -h, --help  show this help message and exit
2726
2727    """
2728    proj = G2Project(args.gpxfile)
2729    phase = proj.phase(args.phase)
2730    phase.export_CIF(args.exportfile)
2731
2732
2733def _args_kwargs(*args, **kwargs):
2734    return args, kwargs
2735
2736# A dictionary of the name of each subcommand, and a tuple
2737# of its associated function and a list of its arguments
2738# The arguments are passed directly to the add_argument() method
2739# of an argparse.ArgumentParser
2740
2741subcommands = {"create":
2742               (create, [_args_kwargs('filename',
2743                                      help='the project file to create. should end in .gpx'),
2744
2745                         _args_kwargs('-d', '--histograms',
2746                                      nargs='+',
2747                                      help='list of datafiles to add as histograms'),
2748                                     
2749                         _args_kwargs('-i', '--iparams',
2750                                      nargs='+',
2751                                      help='instrument parameter file, must be one'
2752                                           ' for every histogram'
2753                                      ),
2754
2755                         _args_kwargs('-p', '--phases',
2756                                      nargs='+',
2757                                      help='list of phases to add. phases are '
2758                                           'automatically associated with all '
2759                                           'histograms given.')]),
2760               "add": (add, [_args_kwargs('filename',
2761                                      help='the project file to open. Should end in .gpx'),
2762
2763                         _args_kwargs('-d', '--histograms',
2764                                      nargs='+',
2765                                      help='list of datafiles to add as histograms'),
2766                                     
2767                         _args_kwargs('-i', '--iparams',
2768                                      nargs='+',
2769                                      help='instrument parameter file, must be one'
2770                                           ' for every histogram'
2771                                      ),
2772                                     
2773                         _args_kwargs('-hf', '--histogramformat',
2774                                      help='format hint for histogram import. Applies to all'
2775                                           ' histograms'
2776                                      ),
2777
2778                         _args_kwargs('-p', '--phases',
2779                                      nargs='+',
2780                                      help='list of phases to add. phases are '
2781                                           'automatically associated with all '
2782                                           'histograms given.'),
2783
2784                         _args_kwargs('-pf', '--phaseformat',
2785                                      help='format hint for phase import. Applies to all'
2786                                           ' phases. Example: -pf CIF'
2787                                      ),
2788                                     
2789                         _args_kwargs('-l', '--histlist',
2790                                      nargs='+',
2791                                      help='list of histgram indices to associate with added'
2792                                           ' phases. If not specified, phases are'
2793                                           ' associated with all previously loaded'
2794                                           ' histograms. Example: -l 2 3 4')]),
2795                                           
2796               "dump": (dump, [_args_kwargs('-d', '--histograms',
2797                                     action='store_true',
2798                                     help='list histograms in files, overrides --raw'),
2799
2800                               _args_kwargs('-p', '--phases',
2801                                            action='store_true',
2802                                            help='list phases in files, overrides --raw'),
2803
2804                               _args_kwargs('-r', '--raw',
2805                                      action='store_true', help='dump raw file contents, default'),
2806
2807                               _args_kwargs('files', nargs='+')]),
2808
2809               "refine":
2810               (refine, [_args_kwargs('gpxfile', help='the project file to refine'),
2811                         _args_kwargs('refinements',
2812                                      help='JSON file of refinements to apply. if not present'
2813                                           ' refines file as-is',
2814                                      default=None,
2815                                      nargs='?')]),
2816
2817               "seqrefine": (seqrefine, []),
2818               "export": (export, [_args_kwargs('gpxfile',
2819                                                help='the project file from which to export'),
2820                                   _args_kwargs('phase', help='identifier of phase to export'),
2821                                   _args_kwargs('exportfile', help='the .cif file to export to')]),
2822               "browse": (IPyBrowse, [_args_kwargs('files', nargs='+',
2823                                                   help='list of files to browse')])}
2824
2825
2826def main():
2827    '''The command-line interface for calling GSASIIscriptable as a shell command,
2828    where it is expected to be called as::
2829
2830       python GSASIIscriptable.py <subcommand> <file.gpx> <options>
2831
2832    The following subcommands are defined:
2833
2834        * create, see :func:`create`
2835        * add, see :func:`add`
2836        * dump, see :func:`dump`
2837        * refine, see :func:`refine`
2838        * seqrefine, see :func:`seqrefine`
2839        * export, :func:`export`
2840        * browse, see :func:`IPyBrowse`
2841
2842    .. seealso::
2843        :func:`create`
2844        :func:`add`
2845        :func:`dump`
2846        :func:`refine`
2847        :func:`seqrefine`
2848        :func:`export`
2849        :func:`IPyBrowse`
2850    '''
2851    parser = argparse.ArgumentParser(description=
2852        "Use of "+os.path.split(__file__)[1]+" Allows GSAS-II actions from command line."
2853        )
2854    subs = parser.add_subparsers()
2855
2856    # Create all of the specified subparsers
2857    for name, (func, args) in subcommands.items():
2858        new_parser = subs.add_parser(name,help=commandhelp.get(name),
2859                                     description='Command "'+name+'" '+commandhelp.get(name))
2860        for listargs, kwds in args:
2861            new_parser.add_argument(*listargs, **kwds)
2862        new_parser.set_defaults(func=func)
2863
2864    # Parse and trigger subcommand
2865    result = parser.parse_args()
2866    result.func(result)
2867
2868if __name__ == '__main__':
2869    main()
Note: See TracBrowser for help on using the repository browser.