source: trunk/GSASIIscriptable.py @ 3129

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

bug fixes

File size: 98.6 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 phase and powder histogram in a :class:`G2Project` object has an associated
66object. Parameters 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 get_wR(self):
1744        """returns the overall weighted profile R factor for a histogram
1745       
1746        :returns: a wR value as a percentage or None if not defined
1747        """
1748        return self['data'][0].get('wR')
1749
1750    def set_refinements(self, refs):
1751        """Sets the refinement parameter 'key' to the specification 'value'
1752
1753        :param dict refs: A dictionary of the parameters to be set. See
1754                          :ref:`Histogram_parameters_table` for a description of
1755                          what these dictionaries should be.
1756
1757        :returns: None
1758
1759        """
1760        do_fit_fixed_points = False
1761        for key, value in refs.items():
1762            if key == 'Limits':
1763                old_limits = self.data['Limits'][1]
1764                new_limits = value
1765                if isinstance(new_limits, dict):
1766                    if 'low' in new_limits:
1767                        old_limits[0] = new_limits['low']
1768                    if 'high' in new_limits:
1769                        old_limits[1] = new_limits['high']
1770                else:
1771                    old_limits[0], old_limits[1] = new_limits
1772            elif key == 'Sample Parameters':
1773                sample = self.data['Sample Parameters']
1774                for sparam in value:
1775                    if sparam not in sample:
1776                        raise ValueError("Unknown refinement parameter, "
1777                                         + str(sparam))
1778                    sample[sparam][1] = True
1779            elif key == 'Background':
1780                bkg, peaks = self.data['Background']
1781
1782                # If True or False, just set the refine parameter
1783                if value in (True, False):
1784                    bkg[1] = value
1785                    return
1786
1787                if 'type' in value:
1788                    bkg[0] = value['type']
1789                if 'refine' in value:
1790                    bkg[1] = value['refine']
1791                if 'no. coeffs' in value:
1792                    cur_coeffs = bkg[2]
1793                    n_coeffs = value['no. coeffs']
1794                    if n_coeffs > cur_coeffs:
1795                        for x in range(n_coeffs - cur_coeffs):
1796                            bkg.append(0.0)
1797                    else:
1798                        for _ in range(cur_coeffs - n_coeffs):
1799                            bkg.pop()
1800                    bkg[2] = n_coeffs
1801                if 'coeffs' in value:
1802                    bkg[3:] = value['coeffs']
1803                if 'FixedPoints' in value:
1804                    peaks['FixedPoints'] = [(float(a), float(b))
1805                                            for a, b in value['FixedPoints']]
1806                if value.get('fit fixed points', False):
1807                    do_fit_fixed_points = True
1808
1809            elif key == 'Instrument Parameters':
1810                instrument, secondary = self.data['Instrument Parameters']
1811                for iparam in value:
1812                    try:
1813                        instrument[iparam][2] = True
1814                    except IndexError:
1815                        raise ValueError("Invalid key:", iparam)
1816            else:
1817                raise ValueError("Unknown key:", key)
1818        # Fit fixed points after the fact - ensure they are after fixed points
1819        # are added
1820        if do_fit_fixed_points:
1821            # Background won't be fit if refinement flag not set
1822            orig = self.data['Background'][0][1]
1823            self.data['Background'][0][1] = True
1824            self.fit_fixed_points()
1825            # Restore the previous value
1826            self.data['Background'][0][1] = orig
1827
1828    def clear_refinements(self, refs):
1829        """Clears the refinement parameter 'key' and its associated value.
1830
1831        :param dict refs: A dictionary of parameters to clear."""
1832        for key, value in refs.items():
1833            if key == 'Limits':
1834                old_limits, cur_limits = self.data['Limits']
1835                cur_limits[0], cur_limits[1] = old_limits
1836            elif key == 'Sample Parameters':
1837                sample = self.data['Sample Parameters']
1838                for sparam in value:
1839                    sample[sparam][1] = False
1840            elif key == 'Background':
1841                bkg, peaks = self.data['Background']
1842
1843                # If True or False, just set the refine parameter
1844                if value in (True, False):
1845                    bkg[1] = False
1846                    return
1847
1848                bkg[1] = False
1849                if 'FixedPoints' in value:
1850                    if 'FixedPoints' in peaks:
1851                        del peaks['FixedPoints']
1852            elif key == 'Instrument Parameters':
1853                instrument, secondary = self.data['Instrument Parameters']
1854                for iparam in value:
1855                    instrument[iparam][2] = False
1856            else:
1857                raise ValueError("Unknown key:", key)
1858
1859
1860class G2Phase(G2ObjectWrapper):
1861    """A wrapper object around a given phase.
1862
1863    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1864    """
1865    def __init__(self, data, name, proj):
1866        self.data = data
1867        self.name = name
1868        self.proj = proj
1869
1870    @staticmethod
1871    def is_valid_refinement_key(key):
1872        valid_keys = ["Cell", "Atoms", "LeBail"]
1873        return key in valid_keys
1874
1875    @staticmethod
1876    def is_valid_HAP_refinement_key(key):
1877        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
1878                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
1879        return key in valid_keys
1880
1881    def atom(self, atomlabel):
1882        """Returns the atom specified by atomlabel, or None if it does not
1883        exist.
1884
1885        :param str atomlabel: The name of the atom (e.g. "O2")
1886        :returns: A :class:`G2AtomRecord` object
1887            representing the atom.
1888        """
1889        # Consult GSASIIobj.py for the meaning of this
1890        cx, ct, cs, cia = self.data['General']['AtomPtrs']
1891        ptrs = [cx, ct, cs, cia]
1892        atoms = self.data['Atoms']
1893        for atom in atoms:
1894            if atom[ct-1] == atomlabel:
1895                return G2AtomRecord(atom, ptrs, self.proj)
1896
1897    def atoms(self):
1898        """Returns a list of atoms present in the phase.
1899
1900        :returns: A list of :class:`G2AtomRecord` objects.
1901
1902        .. seealso::
1903            :meth:`G2Phase.atom`
1904            :class:`G2AtomRecord`
1905        """
1906        ptrs = self.data['General']['AtomPtrs']
1907        output = []
1908        atoms = self.data['Atoms']
1909        for atom in atoms:
1910            output.append(G2AtomRecord(atom, ptrs, self.proj))
1911        return output
1912
1913    def histograms(self):
1914        output = []
1915        for hname in self.data.get('Histograms', {}).keys():
1916            output.append(self.proj.histogram(hname))
1917        return output
1918
1919    @property
1920    def ranId(self):
1921        return self.data['ranId']
1922
1923    @property
1924    def id(self):
1925        return self.data['pId']
1926
1927    @id.setter
1928    def id(self, val):
1929        self.data['pId'] = val
1930
1931    def get_cell(self):
1932        """Returns a dictionary of the cell parameters, with keys:
1933            'length_a', 'length_b', 'length_c', 'angle_alpha', 'angle_beta', 'angle_gamma', 'volume'
1934
1935        :returns: a dict
1936
1937        .. seealso::
1938           :meth:`G2Phase.get_cell_and_esd`
1939
1940        """
1941        cell = self.data['General']['Cell']
1942        return {'length_a': cell[1], 'length_b': cell[2], 'length_c': cell[3],
1943                'angle_alpha': cell[4], 'angle_beta': cell[5], 'angle_gamma': cell[6],
1944                'volume': cell[7]}
1945
1946    def get_cell_and_esd(self):
1947        """
1948        Returns a pair of dictionaries, the first representing the unit cell, the second
1949        representing the estimated standard deviations of the unit cell.
1950
1951        :returns: a tuple of two dictionaries
1952
1953        .. seealso::
1954           :meth:`G2Phase.get_cell`
1955
1956        """
1957        # translated from GSASIIstrIO.ExportBaseclass.GetCell
1958        import GSASIIstrIO as G2stIO
1959        import GSASIIlattice as G2lat
1960        import GSASIImapvars as G2mv
1961        try:
1962            pfx = str(self.id) + '::'
1963            sgdata = self['General']['SGData']
1964            covDict = self.proj['Covariance']['data']
1965
1966            parmDict = dict(zip(covDict.get('varyList',[]),
1967                                covDict.get('variables',[])))
1968            sigDict = dict(zip(covDict.get('varyList',[]),
1969                               covDict.get('sig',[])))
1970
1971            if covDict.get('covMatrix') is not None:
1972                sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
1973                                                  covDict['varyList'],
1974                                                  parmDict))
1975
1976            A, sigA = G2stIO.cellFill(pfx, sgdata, parmDict, sigDict)
1977            cellSig = G2stIO.getCellEsd(pfx, sgdata, A, self.proj['Covariance']['data'])
1978            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
1979            cellDict, cellSigDict = {}, {}
1980            for i, key in enumerate(['length_a', 'length_b', 'length_c',
1981                                     'angle_alpha', 'angle_beta', 'angle_gamma',
1982                                     'volume']):
1983                cellDict[key] = cellList[i]
1984                cellSigDict[key] = cellSig[i]
1985            return cellDict, cellSigDict
1986        except KeyError:
1987            cell = self.get_cell()
1988            return cell, {key: 0.0 for key in cell}
1989
1990    def export_CIF(self, outputname, quickmode=True):
1991        """Write this phase to a .cif file named outputname
1992
1993        :param str outputname: The name of the .cif file to write to
1994        :param bool quickmode: Currently ignored. Carryover from exports.G2export_CIF"""
1995        # This code is all taken from exports/G2export_CIF.py
1996        # Functions copied have the same names
1997        import GSASIImath as G2mth
1998        import GSASIImapvars as G2mv
1999        from exports import G2export_CIF as cif
2000
2001        CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
2002        CIFname = os.path.splitext(self.proj.filename)[0]
2003        CIFname = os.path.split(CIFname)[1]
2004        CIFname = ''.join([c if ord(c) < 128 else ''
2005                           for c in CIFname.replace(' ', '_')])
2006        try:
2007            author = self.proj['Controls']['data'].get('Author','').strip()
2008        except KeyError:
2009            pass
2010        oneblock = True
2011
2012        covDict = self.proj['Covariance']['data']
2013        parmDict = dict(zip(covDict.get('varyList',[]),
2014                            covDict.get('variables',[])))
2015        sigDict = dict(zip(covDict.get('varyList',[]),
2016                           covDict.get('sig',[])))
2017
2018        if covDict.get('covMatrix') is not None:
2019            sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
2020                                              covDict['varyList'],
2021                                              parmDict))
2022
2023        with open(outputname, 'w') as fp:
2024            fp.write(' \n' + 70*'#' + '\n')
2025            cif.WriteCIFitem(fp, 'data_' + CIFname)
2026            # from exports.G2export_CIF.WritePhaseInfo
2027            cif.WriteCIFitem(fp, '\n# phase info for '+str(self.name) + ' follows')
2028            cif.WriteCIFitem(fp, '_pd_phase_name', self.name)
2029            # TODO get esds
2030            cellDict = self.get_cell()
2031            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
2032            names = ['length_a','length_b','length_c',
2033                     'angle_alpha','angle_beta ','angle_gamma',
2034                     'volume']
2035            for key, val in cellDict.items():
2036                cif.WriteCIFitem(fp, '_cell_' + key, G2mth.ValEsd(val))
2037
2038            cif.WriteCIFitem(fp, '_symmetry_cell_setting',
2039                         self.data['General']['SGData']['SGSys'])
2040
2041            spacegroup = self.data['General']['SGData']['SpGrp'].strip()
2042            # regularize capitalization and remove trailing H/R
2043            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
2044            cif.WriteCIFitem(fp, '_symmetry_space_group_name_H-M', spacegroup)
2045
2046            # generate symmetry operations including centering and center of symmetry
2047            SymOpList, offsetList, symOpList, G2oprList, G2opcodes = G2spc.AllOps(
2048                self.data['General']['SGData'])
2049            cif.WriteCIFitem(fp, 'loop_\n    _space_group_symop_id\n    _space_group_symop_operation_xyz')
2050            for i, op in enumerate(SymOpList,start=1):
2051                cif.WriteCIFitem(fp, '   {:3d}  {:}'.format(i,op.lower()))
2052
2053            # TODO skipped histograms, exports/G2export_CIF.py:880
2054
2055            # report atom params
2056            if self.data['General']['Type'] in ['nuclear','macromolecular']:        #this needs macromolecular variant, etc!
2057                cif.WriteAtomsNuclear(fp, self.data, self.name, parmDict, sigDict, [])
2058                # self._WriteAtomsNuclear(fp, parmDict, sigDict)
2059            else:
2060                raise Exception("no export for "+str(self.data['General']['Type'])+" coordinates implemented")
2061            # report cell contents
2062            cif.WriteComposition(fp, self.data, self.name, parmDict)
2063            if not quickmode and self.data['General']['Type'] == 'nuclear':      # report distances and angles
2064                # WriteDistances(fp,self.name,SymOpList,offsetList,symOpList,G2oprList)
2065                raise NotImplementedError("only quickmode currently supported")
2066            if 'Map' in self.data['General'] and 'minmax' in self.data['General']['Map']:
2067                cif.WriteCIFitem(fp,'\n# Difference density results')
2068                MinMax = self.data['General']['Map']['minmax']
2069                cif.WriteCIFitem(fp,'_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009))
2070                cif.WriteCIFitem(fp,'_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009))
2071
2072
2073    def set_refinements(self, refs):
2074        """Sets the refinement parameter 'key' to the specification 'value'
2075
2076        :param dict refs: A dictionary of the parameters to be set. See
2077                          :ref:`Phase_parameters_table` for a description of
2078                          this dictionary.
2079
2080        :returns: None"""
2081        for key, value in refs.items():
2082            if key == "Cell":
2083                self.data['General']['Cell'][0] = value
2084
2085            elif key == "Atoms":
2086                for atomlabel, atomrefinement in value.items():
2087                    if atomlabel == 'all':
2088                        for atom in self.atoms():
2089                            atom.refinement_flags = atomrefinement
2090                    else:
2091                        atom = self.atom(atomlabel)
2092                        if atom is None:
2093                            raise ValueError("No such atom: " + atomlabel)
2094                        atom.refinement_flags = atomrefinement
2095
2096            elif key == "LeBail":
2097                hists = self.data['Histograms']
2098                for hname, hoptions in hists.items():
2099                    if 'LeBail' not in hoptions:
2100                        hoptions['newLeBail'] = bool(True)
2101                    hoptions['LeBail'] = bool(value)
2102            else:
2103                raise ValueError("Unknown key:", key)
2104
2105    def clear_refinements(self, refs):
2106        """Clears a given set of parameters.
2107
2108        :param dict refs: The parameters to clear"""
2109        for key, value in refs.items():
2110            if key == "Cell":
2111                self.data['General']['Cell'][0] = False
2112            elif key == "Atoms":
2113                cx, ct, cs, cia = self.data['General']['AtomPtrs']
2114
2115                for atomlabel in value:
2116                    atom = self.atom(atomlabel)
2117                    # Set refinement to none
2118                    atom.refinement_flags = ' '
2119            elif key == "LeBail":
2120                hists = self.data['Histograms']
2121                for hname, hoptions in hists.items():
2122                    if 'LeBail' not in hoptions:
2123                        hoptions['newLeBail'] = True
2124                    hoptions['LeBail'] = False
2125            else:
2126                raise ValueError("Unknown key:", key)
2127
2128    def set_HAP_refinements(self, refs, histograms='all'):
2129        """Sets the given HAP refinement parameters between this phase and
2130        the given histograms
2131
2132        :param dict refs: A dictionary of the parameters to be set. See
2133                          :ref:`HAP_parameters_table` for a description of this
2134                          dictionary.
2135        :param histograms: Either 'all' (default) or a list of the histograms
2136            whose HAP parameters will be set with this phase. Histogram and phase
2137            must already be associated
2138
2139        :returns: None
2140        """
2141        if histograms == 'all':
2142            histograms = self.data['Histograms'].values()
2143        else:
2144            histograms = [self.data['Histograms'][h.name] for h in histograms]
2145
2146        for key, val in refs.items():
2147            for h in histograms:
2148                if key == 'Babinet':
2149                    try:
2150                        sets = list(val)
2151                    except ValueError:
2152                        sets = ['BabA', 'BabU']
2153                    for param in sets:
2154                        if param not in ['BabA', 'BabU']:
2155                            raise ValueError("Not sure what to do with" + param)
2156                        for hist in histograms:
2157                            hist['Babinet'][param][1] = True
2158                elif key == 'Extinction':
2159                    for h in histograms:
2160                        h['Extinction'][1] = bool(val)
2161                elif key == 'HStrain':
2162                    for h in histograms:
2163                        h['HStrain'][1] = [bool(val) for p in h['HStrain'][1]]
2164                elif key == 'Mustrain':
2165                    for h in histograms:
2166                        mustrain = h['Mustrain']
2167                        newType = None
2168                        if isinstance(val, (unicode, str)):
2169                            if val in ['isotropic', 'uniaxial', 'generalized']:
2170                                newType = val
2171                            else:
2172                                raise ValueError("Not a Mustrain type: " + val)
2173                        elif isinstance(val, dict):
2174                            newType = val.get('type', None)
2175                            direction = val.get('direction', None)
2176
2177                        if newType:
2178                            mustrain[0] = newType
2179                            if newType == 'isotropic':
2180                                mustrain[2][0] = True
2181                                mustrain[5] = [False for p in mustrain[4]]
2182                            elif newType == 'uniaxial':
2183                                if 'refine' in val:
2184                                    types = val['refine']
2185                                    if isinstance(types, (unicode, str)):
2186                                        types = [types]
2187                                    elif isinstance(types, bool):
2188                                        mustrain[2][0] = types
2189                                        mustrain[2][1] = types
2190                                        types = []
2191                                    else:
2192                                        raise ValueError("Not sure what to do with: "
2193                                                         + str(types))
2194                                else:
2195                                    types = []
2196
2197                                for unitype in types:
2198                                    if unitype == 'equatorial':
2199                                        mustrain[2][0] = True
2200                                    elif unitype == 'axial':
2201                                        mustrain[2][1] = True
2202                                    else:
2203                                        msg = 'Invalid uniaxial mustrain type'
2204                                        raise ValueError(msg + ': ' + unitype)
2205                            else:  # newtype == 'generalized'
2206                                mustrain[2] = [False for p in mustrain[1]]
2207
2208                        if direction:
2209                            direction = [int(n) for n in direction]
2210                            if len(direction) != 3:
2211                                raise ValueError("Expected hkl, found", direction)
2212                            mustrain[3] = direction
2213                elif key == 'Pref.Ori.':
2214                    for h in histograms:
2215                        h['Pref.Ori.'][2] = bool(val)
2216                elif key == 'Show':
2217                    for h in histograms:
2218                        h['Show'] = bool(val)
2219                elif key == 'Size':
2220                    for h in histograms:
2221                        if h['Size'][0] == 'isotropic':
2222                            h['Size'][2][0] = bool(val)
2223                        elif h['Size'][0] == 'uniaxial':
2224                            h['Size'][2][1] = bool(val)
2225                            h['Size'][2][2] = bool(val)
2226                        else:   # TODO
2227                            raise NotImplementedError()
2228                elif key == 'Use':
2229                    for h in histograms:
2230                        h['Use'] = bool(val)
2231                elif key == 'Scale':
2232                    for h in histograms:
2233                        h['Scale'][1] = bool(val)
2234                else:
2235                    print(u'Unknown HAP key: '+key)
2236
2237    def clear_HAP_refinements(self, refs, histograms='all'):
2238        """Clears the given HAP refinement parameters between this phase and
2239        the given histograms
2240
2241        :param dict refs: A dictionary of the parameters to be cleared.
2242        :param histograms: Either 'all' (default) or a list of the histograms
2243            whose HAP parameters will be cleared with this phase. Histogram and
2244            phase must already be associated
2245
2246        :returns: None
2247        """
2248        if histograms == 'all':
2249            histograms = self.data['Histograms'].values()
2250        else:
2251            histograms = [self.data['Histograms'][h.name] for h in histograms]
2252
2253        for key, val in refs.items():
2254            for h in histograms:
2255                if key == 'Babinet':
2256                    try:
2257                        sets = list(val)
2258                    except ValueError:
2259                        sets = ['BabA', 'BabU']
2260                    for param in sets:
2261                        if param not in ['BabA', 'BabU']:
2262                            raise ValueError("Not sure what to do with" + param)
2263                        for hist in histograms:
2264                            hist['Babinet'][param][1] = False
2265                elif key == 'Extinction':
2266                    for h in histograms:
2267                        h['Extinction'][1] = False
2268                elif key == 'HStrain':
2269                    for h in histograms:
2270                        h['HStrain'][1] = [False for p in h['HStrain'][1]]
2271                elif key == 'Mustrain':
2272                    for h in histograms:
2273                        mustrain = h['Mustrain']
2274                        mustrain[2] = [False for p in mustrain[1]]
2275                        mustrain[5] = [False for p in mustrain[4]]
2276                elif key == 'Pref.Ori.':
2277                    for h in histograms:
2278                        h['Pref.Ori.'][2] = False
2279                elif key == 'Show':
2280                    for h in histograms:
2281                        h['Show'] = False
2282                elif key == 'Size':
2283                    raise NotImplementedError()
2284                elif key == 'Use':
2285                    for h in histograms:
2286                        h['Use'] = False
2287                elif key == 'Scale':
2288                    for h in histograms:
2289                        h['Scale'][1] = False
2290                else:
2291                    print(u'Unknown HAP key: '+key)
2292
2293
2294##########################
2295# Command Line Interface #
2296##########################
2297# Each of these takes an argparse.Namespace object as their argument,
2298# representing the parsed command-line arguments for the relevant subcommand.
2299# The argument specification for each is in the subcommands dictionary (see
2300# below)
2301
2302
2303def create(args):
2304    """The create subcommand."""
2305    proj = G2Project(filename=args.filename)
2306
2307    hist_objs = []
2308    for h in args.histograms:
2309        hist_objs.append(proj.add_powder_histogram(h, args.iparams))
2310
2311    for p in args.phases:
2312        proj.add_phase(p, histograms=hist_objs)
2313
2314    proj.save()
2315
2316
2317def dump(args):
2318    """The dump subcommand. This is intended to be called by the :func:`main` routine.
2319    This is typically called by invoking this script with a subcommand::
2320
2321       python GSASIIscriptable.py dump <file.gpx>
2322    """
2323    if not args.histograms and not args.phases:
2324        args.raw = True
2325    if args.raw:
2326        import IPython.lib.pretty as pretty
2327
2328    for fname in args.files:
2329        if args.raw:
2330            proj, nameList = LoadDictFromProjFile(fname)
2331            print("file:", fname)
2332            print("names:", nameList)
2333            for key, val in proj.items():
2334                print(key, ":")
2335                pretty.pprint(val)
2336        else:
2337            proj = G2Project(fname)
2338            if args.histograms:
2339                hists = proj.histograms()
2340                for h in hists:
2341                    print(fname, "hist", h.id, h.name)
2342            if args.phases:
2343                phase_list = proj.phases()
2344                for p in phase_list:
2345                    print(fname, "phase", p.id, p.name)
2346
2347
2348def IPyBrowse(args):
2349    """Load a .gpx file and then open a IPython shell to browse it
2350    """
2351    for fname in args.files:
2352        proj, nameList = LoadDictFromProjFile(fname)
2353        msg = "\nfname {} loaded into proj (dict) with names in nameList".format(fname)
2354        GSASIIpath.IPyBreak_base(msg)
2355        break
2356
2357
2358def refine(args):
2359    """The refine subcommand"""
2360    proj = G2Project(args.gpxfile)
2361    if args.refinements is None:
2362        proj.refine()
2363    else:
2364        import json
2365        with open(args.refinements) as refs:
2366            refs = json.load(refs)
2367        proj.do_refinements(refs['refinements'])
2368
2369
2370def seqrefine(args):
2371    """The seqrefine subcommand"""
2372    raise NotImplementedError("seqrefine is not yet implemented")
2373
2374
2375def export(args):
2376    """The export subcommand"""
2377    # Export phase as CIF to args.exportfile
2378    proj = G2Project(args.gpxfile)
2379    phase = proj.phase(args.phase)
2380    phase.export_CIF(args.exportfile)
2381
2382
2383def _args_kwargs(*args, **kwargs):
2384    return args, kwargs
2385
2386# A dictionary of the name of each subcommand, and a tuple
2387# of its associated function and a list of its arguments
2388# The arguments are passed directly to the add_argument() method
2389# of an argparse.ArgumentParser
2390subcommands = {"create":
2391               (create, [_args_kwargs('filename',
2392                                      help='the project file to create. should end in .gpx'),
2393
2394                         _args_kwargs('-i', '--iparams',
2395                                      help='instrument parameter file'),
2396
2397                         _args_kwargs('-d', '--histograms',
2398                                      nargs='+',
2399                                      help='list of datafiles to add as histograms'),
2400
2401                         _args_kwargs('-p', '--phases',
2402                                      nargs='+',
2403                                      help='list of phases to add. phases are '
2404                                           'automatically associated with all '
2405                                           'histograms given.')]),
2406
2407               "dump": (dump, [_args_kwargs('-d', '--histograms',
2408                                     action='store_true',
2409                                     help='list histograms in files, overrides --raw'),
2410
2411                               _args_kwargs('-p', '--phases',
2412                                            action='store_true',
2413                                            help='list phases in files, overrides --raw'),
2414
2415                               _args_kwargs('-r', '--raw',
2416                                      action='store_true', help='dump raw file contents, default'),
2417
2418                               _args_kwargs('files', nargs='+')]),
2419
2420               "refine":
2421               (refine, [_args_kwargs('gpxfile', help='the project file to refine'),
2422                         _args_kwargs('refinements',
2423                                      help='json file of refinements to apply. if not present'
2424                                           ' refines file as-is',
2425                                      default=None,
2426                                      nargs='?')]),
2427
2428               "seqrefine": (seqrefine, []),
2429               "export": (export, [_args_kwargs('gpxfile',
2430                                                help='the project file from which to export'),
2431                                   _args_kwargs('phase', help='identifier of phase to export'),
2432                                   _args_kwargs('exportfile', help='the .cif file to export to')]),
2433               "browse": (IPyBrowse, [_args_kwargs('files', nargs='+',
2434                                                   help='list of files to browse')])}
2435
2436
2437def main():
2438    '''The command line interface for GSASIIscriptable.
2439
2440    Executes one of the following subcommands:
2441
2442        * :func:`create`
2443        * :func:`dump`
2444        * :func:`refine`
2445        * :func:`seqrefine`
2446        * :func:`export`
2447        * browse (:func:`IPyBrowse`)
2448
2449    These commands are typically called by invoking this script with a subcommand,
2450    for example::
2451
2452       python GSASIIscriptable.py dump <file.gpx>
2453       
2454
2455    .. seealso::
2456        :func:`create`
2457        :func:`dump`
2458        :func:`refine`
2459        :func:`seqrefine`
2460        :func:`export`
2461        :func:`IPyBrowse`
2462    '''
2463    parser = argparse.ArgumentParser()
2464    subs = parser.add_subparsers()
2465
2466    # Create all of the specified subparsers
2467    for name, (func, args) in subcommands.items():
2468        new_parser = subs.add_parser(name)
2469        for listargs, kwds in args:
2470            new_parser.add_argument(*listargs, **kwds)
2471        new_parser.set_defaults(func=func)
2472
2473    # Parse and trigger subcommand
2474    result = parser.parse_args()
2475    result.func(result)
2476
2477if __name__ == '__main__':
2478    main()
Note: See TracBrowser for help on using the repository browser.