source: trunk/GSASIIscriptable.py @ 3067

Last change on this file since 3067 was 3067, checked in by odonnell, 5 years ago

Remove some index_id calls

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