source: trunk/GSASIIscriptable.py @ 3213

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

fix EXP import with extra chars on line

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