source: trunk/GSASIIscriptable.py @ 3130

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

updates for G2Scriptable and associated tutorial

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