source: trunk/GSASIIscriptable.py @ 3123

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

do update when load of binaries fails; improve docs and bug fixes in G2scriptable

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