source: trunk/GSASIIscriptable.py @ 3208

Last change on this file since 3208 was 3208, checked in by vondreele, 4 years ago

fix transformation bug abc -> abc* needed transform of matrix - G2phsGUI line 373
work on modullated magnetic cif import; fix some magnetic & incommensurate bugs - more to do.

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