source: trunk/GSASIIscriptable.py @ 3223

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

fix histogram export file naming

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