source: trunk/GSASIIscriptable.py @ 3207

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

New command-line tutorial

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