source: trunk/GSASIIscriptable.py @ 3209

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

Add simulation to scriptable & add to tutorial; add sphinx index to ctrlsGUI; remove path from RunGSASII.bat

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