source: trunk/GSASIIscriptable.py @ 3290

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

scripting doc updates; magnification Py3 updates

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