source: trunk/GSASIIscriptable.py @ 3082

Last change on this file since 3082 was 3069, checked in by odonnell, 8 years ago

phase and histogram access works without G2obj
File size: 82.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.
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        try:
906            # see if histname is an id or ranId
907            histname = int(histname)
908        except ValueError:
909            return
910
911        for histogram in self.histograms():
912            if histogram.id == histname or histogram.ranId == histname:
913                return histogram
914
915    def histograms(self):
916        """Return a list of all histograms, as
917        :class:`G2PwdrData` objects
918
919        .. seealso::
920            :func:`~GSASIIscriptable.G2Project.histograms`
921            :func:`~GSASIIscriptable.G2Project.phase`
922            :func:`~GSASIIscriptable.G2Project.phases`
923            """
924        output = []
925        for obj in self.names:
926            if len(obj) > 1 and obj[0] != u'Phases':
927                output.append(self.histogram(obj[0]))
928        return output
929
930    def phase(self, phasename):
931        """
932        Gives an object representing the specified phase in this project.
933
934        :param str phasename: The name of the desired phase. Either the name
935            (str), the phase's ranId, or the phase's index
936        :returns: A :class:`G2Phase` object
937        :raises: KeyError
938
939        .. seealso::
940            :func:`~GSASIIscriptable.G2Project.histograms`
941            :func:`~GSASIIscriptable.G2Project.phase`
942            :func:`~GSASIIscriptable.G2Project.phases`
943            """
944        phases = self.data['Phases']
945        if phasename in phases:
946            return G2Phase(phases[phasename], phasename, self)
947
948        try:
949            # phasename should be phase index or ranId
950            phasename = int(phasename)
951        except ValueError:
952            return
953
954        for phase in self.phases():
955            if phase.id == phasename or phase.ranId == phasename:
956                return phase
957
958    def phases(self):
959        """
960        Returns a list of all the phases in the project.
961
962        :returns: A list of :class:`G2Phase` objects
963
964        .. seealso::
965            :func:`~GSASIIscriptable.G2Project.histogram`
966            :func:`~GSASIIscriptable.G2Project.histograms`
967            :func:`~GSASIIscriptable.G2Project.phase`
968            """
969        for obj in self.names:
970            if obj[0] == 'Phases':
971                return [self.phase(p) for p in obj[1:]]
972        return []
973
974    def update_ids(self):
975        """Makes sure all phases and histograms have proper hId and pId"""
976        # Translated from GetUsedHistogramsAndPhasesfromTree,
977        #   GSASIIdataGUI.py:4107
978        for i, h in enumerate(self.histograms()):
979            h.id = i
980        for i, p in enumerate(self.phases()):
981            p.id = i
982
983    def do_refinements(self, refinements, histogram='all', phase='all',
984                       outputnames=None, makeBack=False):
985        """Conducts a series of refinements. Wrapper around iter_refinements
986
987        :param list refinements: A list of dictionaries defining refinements
988        :param str histogram: Name of histogram for refinements to be applied
989            to, or 'all'
990        :param str phase: Name of phase for refinements to be applied to, or
991            'all'
992        """
993        for proj in self.iter_refinements(refinements, histogram, phase,
994                                          outputnames, makeBack):
995            pass
996        return self
997
998    def iter_refinements(self, refinements, histogram='all', phase='all',
999                         outputnames=None, makeBack=False):
1000        """Conducts a series of refinements, iteratively. Stops after every
1001        refinement and yields this project, to allow error checking or
1002        logging of intermediate results.
1003
1004        >>> def checked_refinements(proj):
1005        ...     for p in proj.iter_refinements(refs):
1006        ...         # Track intermediate results
1007        ...         log(p.histogram('0').residuals)
1008        ...         log(p.phase('0').get_cell())
1009        ...         # Check if parameter diverged, nonsense answer, or whatever
1010        ...         if is_something_wrong(p):
1011        ...             raise Exception("I need a human!")
1012
1013        :param list refinements: A list of dictionaries defining refinements
1014        :param str histogram: Name of histogram for refinements to be applied
1015            to, or 'all'
1016        :param str phase: Name of phase for refinements to be applied to, or
1017            'all'
1018        """
1019        if outputnames:
1020            if len(refinements) != len(outputnames):
1021                raise ValueError("Should have same number of outuputs to"
1022                                 "refinements")
1023        else:
1024            outputnames = [None for r in refinements]
1025
1026        for output, refinement in zip(outputnames, refinements):
1027            self.set_refinement(refinement, histogram)
1028            # Handle 'once' args - refinements that are disabled after this
1029            # refinement
1030            if 'once' in refinement:
1031                temp = {'set': refinement['once']}
1032                self.set_refinement(temp, histogram, phase)
1033
1034            if output:
1035                self.save(output)
1036
1037            self.refine(makeBack=makeBack)
1038            yield self
1039
1040            # Handle 'once' args - refinements that are disabled after this
1041            # refinement
1042            if 'once' in refinement:
1043                temp = {'clear': refinement['once']}
1044                self.set_refinement(temp, histogram, phase)
1045
1046    def set_refinement(self, refinement, histogram='all', phase='all'):
1047        """Apply specified refinements to a given histogram(s) or phase(s).
1048
1049        Refinement parameters are categorize in three groups:
1050
1051        1. Histogram parameters
1052        2. Phase parameters
1053        3. Histogram-and-Phase (HAP) parameters
1054
1055        :param dict refinement: The refinements to be conducted
1056        :param histogram: Either a name of a histogram (str), a list of
1057            histogram names, or 'all' (default)
1058        :param phase: Either a name of a phase (str), a list of phase names, or
1059            'all' (default)
1060
1061        .. seealso::
1062            :func:`~G2PwdrData.set_refinements`
1063            :func:`~G2PwdrData.clear_refinements`
1064            :func:`~G2Phase.set_refinements`
1065            :func:`~G2Phase.clear_refinements`
1066            :func:`~G2Phase.set_HAP_refinements`
1067            :func:`~G2Phase.clear_HAP_refinements`"""
1068
1069        if histogram == 'all':
1070            hists = self.histograms()
1071        elif isinstance(histogram, str) or isinstance(histogram, int):
1072            hists = [self.histogram(histogram)]
1073        else:
1074            hists = [self.histogram(name) for name in histogram]
1075
1076        if phase == 'all':
1077            phases = self.phases()
1078        elif isinstance(phase, str) or isinstance(phase, int):
1079            phases = [self.phase(phase)]
1080        else:
1081            phases = [self.phase(name) for name in phase]
1082
1083        # TODO: HAP parameters:
1084        #   Babinet
1085        #   Extinction
1086        #   HStrain
1087        #   Mustrain
1088        #   Pref. Ori
1089        #   Size
1090
1091        pwdr_set = {}
1092        phase_set = {}
1093        hap_set = {}
1094        for key, val in refinement.get('set', {}).items():
1095            # Apply refinement options
1096            if G2PwdrData.is_valid_refinement_key(key):
1097                pwdr_set[key] = val
1098            elif G2Phase.is_valid_refinement_key(key):
1099                phase_set[key] = val
1100            elif G2Phase.is_valid_HAP_refinement_key(key):
1101                hap_set[key] = val
1102            else:
1103                raise ValueError("Unknown refinement key", key)
1104
1105        for hist in hists:
1106            hist.set_refinements(pwdr_set)
1107        for phase in phases:
1108            phase.set_refinements(phase_set)
1109        for phase in phases:
1110            phase.set_HAP_refinements(hap_set, hists)
1111
1112        pwdr_clear = {}
1113        phase_clear = {}
1114        hap_clear = {}
1115        for key, val in refinement.get('clear', {}).items():
1116            # Clear refinement options
1117            if G2PwdrData.is_valid_refinement_key(key):
1118                pwdr_clear[key] = val
1119            elif G2Phase.is_valid_refinement_key(key):
1120                phase_clear[key] = val
1121            elif G2Phase.is_valid_HAP_refinement_key(key):
1122                hap_set[key] = val
1123            else:
1124                raise ValueError("Unknown refinement key", key)
1125
1126        for hist in hists:
1127            hist.clear_refinements(pwdr_clear)
1128        for phase in phases:
1129            phase.clear_refinements(phase_clear)
1130        for phase in phases:
1131            phase.clear_HAP_refinements(hap_clear, hists)
1132
1133    def index_ids(self):
1134        import GSASIIstrIO as G2strIO
1135        self.save()
1136        return G2strIO.GetUsedHistogramsAndPhases(self.filename)
1137
1138    def add_constraint_raw(self, cons_scope, constr):
1139        """Adds a constraint of type consType to the project.
1140        cons_scope should be one of "Hist", "Phase", "HAP", or "Global".
1141
1142        WARNING it does not check the constraint is well-constructed"""
1143        constrs = self.data['Constraints']['data']
1144        if 'Global' not in constrs:
1145            constrs['Global'] = []
1146        constrs[cons_scope].append(constr)
1147
1148    def hold_many(self, vars, type):
1149        """Apply holds for all the variables in vars, for constraint of a given type.
1150
1151        type is passed directly to add_constraint_raw as consType
1152
1153        :param list vars: A list of variables to hold. Either :class:`GSASIIobj.G2VarObj` objects,
1154            string variable specifiers, or arguments for :meth:`make_var_obj`
1155        :param str type: A string constraint type specifier. See
1156        :class:`~GSASIIscriptable.G2Project.add_constraint_raw`"""
1157        for var in vars:
1158            if isinstance(var, str):
1159                var = self.make_var_obj(var)
1160            elif not isinstance(var, G2obj.G2VarObj):
1161                var = self.make_var_obj(*var)
1162            self.add_constraint_raw(type, [[1.0, var], None, None, 'h'])
1163
1164    def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
1165                     reloadIdx=True):
1166        """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a")
1167        or individual names of phase, histogram, varname, and atomId.
1168
1169        Automatically converts string phase, hist, or atom names into the ID required
1170        by G2VarObj."""
1171
1172        if reloadIdx:
1173            self.index_ids()
1174
1175        # If string representation, short circuit
1176        if hist is None and varname is None and atomId is None:
1177            if isinstance(phase, str) and ':' in phase:
1178                return G2obj.G2VarObj(phase)
1179
1180        # Get phase index
1181        phaseObj = None
1182        if isinstance(phase, G2Phase):
1183            phaseObj = phase
1184            phase = G2obj.PhaseRanIdLookup[phase.ranId]
1185        elif phase in self.data['Phases']:
1186            phaseObj = self.phase(phase)
1187            phaseRanId = phaseObj.ranId
1188            phase = G2obj.PhaseRanIdLookup[phaseRanId]
1189
1190        # Get histogram index
1191        if isinstance(hist, G2PwdrData):
1192            hist = G2obj.HistRanIdLookup[hist.ranId]
1193        elif hist in self.data:
1194            histRanId = self.histogram(hist).ranId
1195            hist = G2obj.HistRanIdLookup[histRanId]
1196
1197        # Get atom index (if any)
1198        if isinstance(atomId, G2AtomRecord):
1199            atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
1200        elif phaseObj:
1201            atomObj = phaseObj.atom(atomId)
1202            if atomObj:
1203                atomRanId = atomObj.ranId
1204                atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
1205
1206        return G2obj.G2VarObj(phase, hist, varname, atomId)
1207
1208
1209class G2AtomRecord(G2ObjectWrapper):
1210    """Wrapper for an atom record. Has convenient accessors via @property.
1211
1212
1213    Available accessors: label, type, refinement_flags, coordinates,
1214        occupancy, ranId, id, adp_flag, uiso
1215
1216    Example:
1217
1218    >>> atom = some_phase.atom("O3")
1219    >>> # We can access the underlying data format
1220    >>> atom.data
1221    ['O3', 'O-2', '', ... ]
1222    >>> # We can also use wrapper accessors
1223    >>> atom.coordinates
1224    (0.33, 0.15, 0.5)
1225    >>> atom.refinement_flags
1226    u'FX'
1227    >>> atom.ranId
1228    4615973324315876477
1229    >>> atom.occupancy
1230    1.0
1231    """
1232    def __init__(self, data, indices, proj):
1233        self.data = data
1234        self.cx, self.ct, self.cs, self.cia = indices
1235        self.proj = proj
1236
1237    @property
1238    def label(self):
1239        return self.data[self.ct-1]
1240
1241    @property
1242    def type(self):
1243        return self.data[self.ct]
1244
1245    @property
1246    def refinement_flags(self):
1247        return self.data[self.ct+1]
1248
1249    @refinement_flags.setter
1250    def refinement_flags(self, other):
1251        # Automatically check it is a valid refinement
1252        for c in other:
1253            if c not in ' FXU':
1254                raise ValueError("Invalid atom refinement: ", other)
1255        self.data[self.ct+1] = unicode(other)
1256
1257    @property
1258    def coordinates(self):
1259        return tuple(self.data[self.cx:self.cx+3])
1260
1261    @property
1262    def occupancy(self):
1263        return self.data[self.cx+3]
1264
1265    @occupancy.setter
1266    def occupancy(self, val):
1267        self.data[self.cx+3] = float(val)
1268
1269    @property
1270    def ranId(self):
1271        return self.data[self.cia+8]
1272
1273    @property
1274    def adp_flag(self):
1275        # Either 'I' or 'A'
1276        return self.data[self.cia]
1277
1278    @property
1279    def uiso(self):
1280        if self.adp_flag == 'I':
1281            return self.data[self.cia+1]
1282        else:
1283            return self.data[self.cia+2:self.cia+8]
1284
1285    @uiso.setter
1286    def uiso(self, value):
1287        if self.adp_flag == 'I':
1288            self.data[self.cia+1] = float(value)
1289        else:
1290            assert len(value) == 6
1291            self.data[self.cia+2:self.cia+8] = [float(v) for v in value]
1292
1293
1294class G2PwdrData(G2ObjectWrapper):
1295    """Wraps a Powder Data Histogram."""
1296    def __init__(self, data, proj):
1297        self.data = data
1298        self.proj = proj
1299
1300    @staticmethod
1301    def is_valid_refinement_key(key):
1302        valid_keys = ['Limits', 'Sample Parameters', 'Background',
1303                      'Instrument Parameters']
1304        return key in valid_keys
1305
1306    @property
1307    def name(self):
1308        return self.data['data'][-1]
1309
1310    @property
1311    def ranId(self):
1312        return self.data['data'][0]['ranId']
1313
1314    @property
1315    def residuals(self):
1316        data = self.data['data'][0]
1317        return {key: data[key]
1318                for key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']}
1319
1320    @property
1321    def id(self):
1322        self.proj.update_ids()
1323        return self.data['data'][0]['hId']
1324
1325    @id.setter
1326    def id(self, val):
1327        self.data['data'][0]['hId'] = val
1328
1329    def fit_fixed_points(self):
1330        """Attempts to apply a background fit to the fixed points currently specified."""
1331        def SetInstParms(Inst):
1332            dataType = Inst['Type'][0]
1333            insVary = []
1334            insNames = []
1335            insVals = []
1336            for parm in Inst:
1337                insNames.append(parm)
1338                insVals.append(Inst[parm][1])
1339                if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1340                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1341                        Inst[parm][2] = False
1342            instDict = dict(zip(insNames, insVals))
1343            instDict['X'] = max(instDict['X'], 0.01)
1344            instDict['Y'] = max(instDict['Y'], 0.01)
1345            if 'SH/L' in instDict:
1346                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
1347            return dataType, instDict, insVary
1348
1349        bgrnd = self.data['Background']
1350
1351        # Need our fixed points in order
1352        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
1353        X = [x for x, y in bgrnd[1]['FixedPoints']]
1354        Y = [y for x, y in bgrnd[1]['FixedPoints']]
1355
1356        limits = self.data['Limits'][1]
1357        if X[0] > limits[0]:
1358            X = [limits[0]] + X
1359            Y = [Y[0]] + Y
1360        if X[-1] < limits[1]:
1361            X += [limits[1]]
1362            Y += [Y[-1]]
1363
1364        # Some simple lookups
1365        controls = self.proj['Controls']['data']
1366        inst, inst2 = self.data['Instrument Parameters']
1367        pwddata = self.data['data'][1]
1368
1369        # Construct the data for background fitting
1370        xBeg = np.searchsorted(pwddata[0], limits[0])
1371        xFin = np.searchsorted(pwddata[0], limits[1])
1372        xdata = pwddata[0][xBeg:xFin]
1373        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
1374
1375        W = [1]*len(xdata)
1376        Z = [0]*len(xdata)
1377
1378        dataType, insDict, insVary = SetInstParms(inst)
1379        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1380
1381        # Do the fit
1382        data = np.array([xdata, ydata, W, Z, Z, Z])
1383        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
1384                        prevVaryList=bakVary, controls=controls)
1385
1386        # Post-fit
1387        parmDict = {}
1388        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1389        parmDict.update(bakDict)
1390        parmDict.update(insDict)
1391        pwddata[3][xBeg:xFin] *= 0
1392        pwddata[5][xBeg:xFin] *= 0
1393        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
1394
1395        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
1396        # TODO update background
1397        self.proj.save()
1398
1399    def y_calc(self):
1400        return self.data['data'][1][3]
1401
1402    def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
1403        try:
1404            import matplotlib.pyplot as plt
1405            data = self.data['data'][1]
1406            if Yobs:
1407                plt.plot(data[0], data[1], label='Yobs')
1408            if Ycalc:
1409                plt.plot(data[0], data[3], label='Ycalc')
1410            if Background:
1411                plt.plot(data[0], data[4], label='Background')
1412            if Residual:
1413                plt.plot(data[0], data[5], label="Residual")
1414        except ImportError:
1415            pass
1416
1417    def set_refinements(self, refs):
1418        """Sets the refinement parameter 'key' to the specification 'value'
1419
1420        :param dict refs: A dictionary of the parameters to be set.
1421
1422        :returns: None"""
1423        do_fit_fixed_points = False
1424        for key, value in refs.items():
1425            if key == 'Limits':
1426                old_limits = self.data['Limits'][1]
1427                new_limits = value
1428                if isinstance(new_limits, dict):
1429                    if 'low' in new_limits:
1430                        old_limits[0] = new_limits['low']
1431                    if 'high' in new_limits:
1432                        old_limits[1] = new_limits['high']
1433                else:
1434                    old_limits[0], old_limits[1] = new_limits
1435            elif key == 'Sample Parameters':
1436                sample = self.data['Sample Parameters']
1437                for sparam in value:
1438                    sample[sparam][1] = True
1439            elif key == 'Background':
1440                bkg, peaks = self.data['Background']
1441
1442                # If True or False, just set the refine parameter
1443                if value in (True, False):
1444                    bkg[1] = value
1445                    return
1446
1447                if 'type' in value:
1448                    bkg[0] = value['type']
1449                if 'refine' in value:
1450                    bkg[1] = value['refine']
1451                if 'no. coeffs' in value:
1452                    cur_coeffs = bkg[2]
1453                    n_coeffs = value['no. coeffs']
1454                    if n_coeffs > cur_coeffs:
1455                        for x in range(n_coeffs - cur_coeffs):
1456                            bkg.append(0.0)
1457                    else:
1458                        for _ in range(cur_coeffs - n_coeffs):
1459                            bkg.pop()
1460                    bkg[2] = n_coeffs
1461                if 'coeffs' in value:
1462                    bkg[3:] = value['coeffs']
1463                if 'FixedPoints' in value:
1464                    peaks['FixedPoints'] = [(float(a), float(b))
1465                                            for a, b in value['FixedPoints']]
1466                if value.get('fit fixed points', False):
1467                    do_fit_fixed_points = True
1468
1469            elif key == 'Instrument Parameters':
1470                instrument, secondary = self.data['Instrument Parameters']
1471                for iparam in value:
1472                    try:
1473                        instrument[iparam][2] = True
1474                    except IndexError:
1475                        raise ValueError("Invalid key:", iparam)
1476            else:
1477                raise ValueError("Unknown key:", key)
1478        # Fit fixed points after the fact - ensure they are after fixed points
1479        # are added
1480        if do_fit_fixed_points:
1481            # Background won't be fit if refinement flag not set
1482            orig = self.data['Background'][0][1]
1483            self.data['Background'][0][1] = True
1484            self.fit_fixed_points()
1485            # Restore the previous value
1486            self.data['Background'][0][1] = orig
1487
1488    def clear_refinements(self, refs):
1489        """Clears the refinement parameter 'key' and its associated value.
1490
1491        :param dict refs: A dictionary of parameters to clear."""
1492        for key, value in refs.items():
1493            if key == 'Limits':
1494                old_limits, cur_limits = self.data['Limits']
1495                cur_limits[0], cur_limits[1] = old_limits
1496            elif key == 'Sample Parameters':
1497                sample = self.data['Sample Parameters']
1498                for sparam in value:
1499                    sample[sparam][1] = False
1500            elif key == 'Background':
1501                bkg, peaks = self.data['Background']
1502
1503                # If True or False, just set the refine parameter
1504                if value in (True, False):
1505                    bkg[1] = False
1506                    return
1507
1508                bkg[1] = False
1509                if 'FixedPoints' in value:
1510                    if 'FixedPoints' in peaks:
1511                        del peaks['FixedPoints']
1512            elif key == 'Instrument Parameters':
1513                instrument, secondary = self.data['Instrument Parameters']
1514                for iparam in value:
1515                    instrument[iparam][2] = False
1516            else:
1517                raise ValueError("Unknown key:", key)
1518
1519
1520class G2Phase(G2ObjectWrapper):
1521    """A wrapper object around a given phase.
1522
1523    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1524    """
1525    def __init__(self, data, name, proj):
1526        self.data = data
1527        self.name = name
1528        self.proj = proj
1529
1530    @staticmethod
1531    def is_valid_refinement_key(key):
1532        valid_keys = ["Cell", "Atoms", "LeBail"]
1533        return key in valid_keys
1534
1535    @staticmethod
1536    def is_valid_HAP_refinement_key(key):
1537        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
1538                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
1539        return key in valid_keys
1540
1541    def atom(self, atomlabel):
1542        """Returns the atom specified by atomlabel, or None if it does not
1543        exist.
1544
1545        :param str atomlabel: The name of the atom (e.g. "O2")
1546        :returns: A :class:`G2AtomRecord` object
1547            representing the atom.
1548        """
1549        # Consult GSASIIobj.py for the meaning of this
1550        cx, ct, cs, cia = self.data['General']['AtomPtrs']
1551        ptrs = [cx, ct, cs, cia]
1552        atoms = self.data['Atoms']
1553        for atom in atoms:
1554            if atom[ct-1] == atomlabel:
1555                return G2AtomRecord(atom, ptrs, self.proj)
1556
1557    def atoms(self):
1558        """Returns a list of atoms present in the phase.
1559
1560        :returns: A list of :class:`G2AtomRecord` objects.
1561
1562        .. seealso::
1563            :func:`~GSASIIscriptable.G2Phase.atom`
1564            :class:`G2AtomRecord`
1565        """
1566        ptrs = self.data['General']['AtomPtrs']
1567        output = []
1568        atoms = self.data['Atoms']
1569        for atom in atoms:
1570            output.append(G2AtomRecord(atom, ptrs, self.proj))
1571        return output
1572
1573    def histograms(self):
1574        output = []
1575        for hname in self.data.get('Histograms', {}).keys():
1576            output.append(self.proj.histogram(hname))
1577        return output
1578
1579    @property
1580    def ranId(self):
1581        return self.data['ranId']
1582
1583    @property
1584    def id(self):
1585        return self.data['pId']
1586
1587    @id.setter
1588    def id(self, val):
1589        self.data['pId'] = val
1590
1591    def get_cell(self):
1592        """Returns a dictionary of the cell parameters, with keys:
1593            'length_a', 'length_b', 'length_c', 'angle_alpha', 'angle_beta', 'angle_gamma', 'volume'
1594
1595        :returns: a dict
1596
1597        .. seealso::
1598           :func:`~G2Phase.get_cell_and_esd`
1599        """
1600        cell = self.data['General']['Cell']
1601        return {'length_a': cell[1], 'length_b': cell[2], 'length_c': cell[3],
1602                'angle_alpha': cell[4], 'angle_beta': cell[5], 'angle_gamma': cell[6],
1603                'volume': cell[7]}
1604
1605    def get_cell_and_esd(self):
1606        """
1607        Returns a pair of dictionaries, the first representing the unit cell, the second
1608        representing the estimated standard deviations of the unit cell.
1609
1610        :returns: a tuple of two dictionaries
1611        .. seealso::
1612           :func:`~G2Phase.get_cell`
1613        """
1614        # translated from GSASIIstrIO.ExportBaseclass.GetCell
1615        import GSASIIstrIO as G2stIO
1616        import GSASIIlattice as G2lat
1617        import GSASIImapvars as G2mv
1618        try:
1619            pfx = str(self.id) + '::'
1620            sgdata = self['General']['SGData']
1621            covDict = self.proj['Covariance']['data']
1622
1623            parmDict = dict(zip(covDict.get('varyList',[]),
1624                                covDict.get('variables',[])))
1625            sigDict = dict(zip(covDict.get('varyList',[]),
1626                               covDict.get('sig',[])))
1627
1628            if covDict.get('covMatrix') is not None:
1629                sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
1630                                                  covDict['varyList'],
1631                                                  parmDict))
1632
1633            A, sigA = G2stIO.cellFill(pfx, sgdata, parmDict, sigDict)
1634            cellSig = G2stIO.getCellEsd(pfx, sgdata, A, self.proj['Covariance']['data'])
1635            cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),)
1636            cellDict, cellSigDict = {}, {}
1637            for i, key in enumerate(['length_a', 'length_b', 'length_c',
1638                                     'angle_alpha', 'angle_beta', 'angle_gamma',
1639                                     'volume']):
1640                cellDict[key] = cellList[i]
1641                cellSigDict[key] = cellSig[i]
1642            return cellDict, cellSigDict
1643        except KeyError:
1644            cell = self.get_cell()
1645            return cell, {key: 0.0 for key in cell}
1646
1647    def export_CIF(self, outputname, quickmode=True):
1648        """Write this phase to a .cif file named outputname
1649
1650        :param str outputname: The name of the .cif file to write to
1651        :param bool quickmode: Currently ignored. Carryover from exports.G2export_CIF"""
1652        # This code is all taken from exports/G2export_CIF.py
1653        # Functions copied have the same names
1654        import GSASIImath as G2mth
1655        import GSASIImapvars as G2mv
1656        from exports import G2export_CIF as cif
1657
1658        CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M")
1659        CIFname = os.path.splitext(self.proj.filename)[0]
1660        CIFname = os.path.split(CIFname)[1]
1661        CIFname = ''.join([c if ord(c) < 128 else ''
1662                           for c in CIFname.replace(' ', '_')])
1663        try:
1664            author = self.proj['Controls']['data'].get('Author','').strip()
1665        except KeyError:
1666            pass
1667        oneblock = True
1668
1669        covDict = self.proj['Covariance']['data']
1670        parmDict = dict(zip(covDict.get('varyList',[]),
1671                            covDict.get('variables',[])))
1672        sigDict = dict(zip(covDict.get('varyList',[]),
1673                           covDict.get('sig',[])))
1674
1675        if covDict.get('covMatrix') is not None:
1676            sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'],
1677                                              covDict['varyList'],
1678                                              parmDict))
1679
1680        with open(outputname, 'w') as fp:
1681            fp.write(' \n' + 70*'#' + '\n')
1682            cif.WriteCIFitem(fp, 'data_' + CIFname)
1683            # from exports.G2export_CIF.WritePhaseInfo
1684            cif.WriteCIFitem(fp, '\n# phase info for '+str(self.name) + ' follows')
1685            cif.WriteCIFitem(fp, '_pd_phase_name', self.name)
1686            # TODO get esds
1687            cellDict = self.get_cell()
1688            defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma
1689            names = ['length_a','length_b','length_c',
1690                     'angle_alpha','angle_beta ','angle_gamma',
1691                     'volume']
1692            for key, val in cellDict.items():
1693                cif.WriteCIFitem(fp, '_cell_' + key, G2mth.ValEsd(val))
1694
1695            cif.WriteCIFitem(fp, '_symmetry_cell_setting',
1696                         self.data['General']['SGData']['SGSys'])
1697
1698            spacegroup = self.data['General']['SGData']['SpGrp'].strip()
1699            # regularize capitalization and remove trailing H/R
1700            spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ')
1701            cif.WriteCIFitem(fp, '_symmetry_space_group_name_H-M', spacegroup)
1702
1703            # generate symmetry operations including centering and center of symmetry
1704            SymOpList, offsetList, symOpList, G2oprList, G2opcodes = G2spc.AllOps(
1705                self.data['General']['SGData'])
1706            cif.WriteCIFitem(fp, 'loop_\n    _space_group_symop_id\n    _space_group_symop_operation_xyz')
1707            for i, op in enumerate(SymOpList,start=1):
1708                cif.WriteCIFitem(fp, '   {:3d}  {:}'.format(i,op.lower()))
1709
1710            # TODO skipped histograms, exports/G2export_CIF.py:880
1711
1712            # report atom params
1713            if self.data['General']['Type'] in ['nuclear','macromolecular']:        #this needs macromolecular variant, etc!
1714                cif.WriteAtomsNuclear(fp, self.data, self.name, parmDict, sigDict, [])
1715                # self._WriteAtomsNuclear(fp, parmDict, sigDict)
1716            else:
1717                raise Exception("no export for "+str(self.data['General']['Type'])+" coordinates implemented")
1718            # report cell contents
1719            cif.WriteComposition(fp, self.data, self.name, parmDict)
1720            if not quickmode and self.data['General']['Type'] == 'nuclear':      # report distances and angles
1721                # WriteDistances(fp,self.name,SymOpList,offsetList,symOpList,G2oprList)
1722                raise NotImplementedError("only quickmode currently supported")
1723            if 'Map' in self.data['General'] and 'minmax' in self.data['General']['Map']:
1724                cif.WriteCIFitem(fp,'\n# Difference density results')
1725                MinMax = self.data['General']['Map']['minmax']
1726                cif.WriteCIFitem(fp,'_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009))
1727                cif.WriteCIFitem(fp,'_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009))
1728
1729
1730    def set_refinements(self, refs):
1731        """Sets the refinement parameter 'key' to the specification 'value'
1732
1733        :param dict refs: A dictionary of the parameters to be set.
1734
1735        :returns: None"""
1736        for key, value in refs.items():
1737            if key == "Cell":
1738                self.data['General']['Cell'][0] = True
1739            elif key == "Atoms":
1740                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1741
1742                for atomlabel, atomrefinement in value.items():
1743                    if atomlabel == 'all':
1744                        for atom in self.atoms():
1745                            atom.refinement_flags = atomrefinement
1746                    else:
1747                        atom = self.atom(atomlabel)
1748                        atom.refinement_flags = atomrefinement
1749            elif key == "LeBail":
1750                hists = self.data['Histograms']
1751                for hname, hoptions in hists.items():
1752                    if 'LeBail' not in hoptions:
1753                        hoptions['newLeBail'] = True
1754                    hoptions['LeBail'] = bool(value)
1755            else:
1756                raise ValueError("Unknown key:", key)
1757
1758    def clear_refinements(self, refs):
1759        """Clears a given set of parameters.
1760
1761        :param dict refs: The parameters to clear"""
1762        for key, value in refs.items():
1763            if key == "Cell":
1764                self.data['General']['Cell'][0] = False
1765            elif key == "Atoms":
1766                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1767
1768                for atomlabel in value:
1769                    atom = self.atom(atomlabel)
1770                    # Set refinement to none
1771                    atom.refinement_flags = ' '
1772            elif key == "LeBail":
1773                hists = self.data['Histograms']
1774                for hname, hoptions in hists.items():
1775                    if 'LeBail' not in hoptions:
1776                        hoptions['newLeBail'] = True
1777                    hoptions['LeBail'] = False
1778            else:
1779                raise ValueError("Unknown key:", key)
1780
1781    def set_HAP_refinements(self, refs, histograms='all'):
1782        """Sets the given HAP refinement parameters between this phase and
1783        the given histograms
1784
1785        :param dict refs: A dictionary of the parameters to be set.
1786        :param histograms: Either 'all' (default) or a list of the histograms
1787            whose HAP parameters will be set with this phase. Histogram and phase
1788            must already be associated
1789
1790        :returns: None
1791        """
1792        if histograms == 'all':
1793            histograms = self.data['Histograms'].values()
1794        else:
1795            histograms = [self.data['Histograms'][h.name] for h in histograms]
1796
1797        for key, val in refs.items():
1798            for h in histograms:
1799                if key == 'Babinet':
1800                    try:
1801                        sets = list(val)
1802                    except ValueError:
1803                        sets = ['BabA', 'BabU']
1804                    for param in sets:
1805                        if param not in ['BabA', 'BabU']:
1806                            raise ValueError("Not sure what to do with" + param)
1807                        for hist in histograms:
1808                            hist['Babinet'][param][1] = True
1809                elif key == 'Extinction':
1810                    for h in histograms:
1811                        h['Extinction'][1] = True
1812                elif key == 'HStrain':
1813                    for h in histograms:
1814                        hist['HStrain'][1] = [True for p in hist['Hstrain'][0]]
1815                elif key == 'Mustrain':
1816                    for h in histograms:
1817                        mustrain = h['Mustrain']
1818                        newType = None
1819                        if isinstance(val, str):
1820                            if val in ['isotropic', 'uniaxial', 'generalized']:
1821                                newType = val
1822                            else:
1823                                raise ValueError("Not a Mustrain type: " + val)
1824                        elif isinstance(val, dict):
1825                            newType = val.get('type', None)
1826                            direction = val.get('direction', None)
1827
1828                        if newType:
1829                            mustrain[0] = newType
1830                            if newType == 'isotropic':
1831                                mustrain[2][0] = True
1832                                mustrain[5] = [False for p in mustrain[4]]
1833                            elif newType == 'uniaxial':
1834                                pass
1835                            else:  # newtype == 'generalized'
1836                                mustrain[2] = [False for p in mustrain[1]]
1837                        if direction:
1838                            direction = [int(n) for n in direction]
1839                            if len(direction) != 3:
1840                                raise ValueError("Expected hkl, found", direction)
1841                            mustrain[3] = direction
1842                elif key == 'Pref.Ori.':
1843                    for h in histograms:
1844                        h['Pref.Ori.'][2] = True
1845                elif key == 'Show':
1846                    for h in histograms:
1847                        h['Show'] = True
1848                elif key == 'Size':
1849                    raise NotImplementedError()
1850                elif key == 'Use':
1851                    for h in histograms:
1852                        h['Use'] = True
1853                elif key == 'Scale':
1854                    for h in histograms:
1855                        h['Scale'][1] = False
1856
1857    def clear_HAP_refinements(self, refs, histograms='all'):
1858        """Clears the given HAP refinement parameters between this phase and
1859        the given histograms
1860
1861        :param dict refs: A dictionary of the parameters to be cleared.
1862        :param histograms: Either 'all' (default) or a list of the histograms
1863            whose HAP parameters will be cleared with this phase. Histogram and
1864            phase must already be associated
1865
1866        :returns: None
1867        """
1868        if histograms == 'all':
1869            histograms = self.data['Histograms'].values()
1870        else:
1871            histograms = [self.data['Histograms'][h.name] for h in histograms]
1872
1873        for key, val in refs.items():
1874            for h in histograms:
1875                if key == 'Babinet':
1876                    try:
1877                        sets = list(val)
1878                    except ValueError:
1879                        sets = ['BabA', 'BabU']
1880                    for param in sets:
1881                        if param not in ['BabA', 'BabU']:
1882                            raise ValueError("Not sure what to do with" + param)
1883                        for hist in histograms:
1884                            hist['Babinet'][param][1] = False
1885                elif key == 'Extinction':
1886                    for h in histograms:
1887                        h['Extinction'][1] = False
1888                elif key == 'HStrain':
1889                    for h in histograms:
1890                        hist['HStrain'][1] = [False for p in hist['Hstrain'][0]]
1891                elif key == 'Mustrain':
1892                    for h in histograms:
1893                        mustrain = h['Mustrain']
1894                        mustrain[2] = [False for p in mustrain[1]]
1895                        mustrain[5] = [False for p in mustrain[4]]
1896                elif key == 'Pref.Ori.':
1897                    for h in histograms:
1898                        h['Pref.Ori.'][2] = False
1899                elif key == 'Show':
1900                    for h in histograms:
1901                        h['Show'] = False
1902                elif key == 'Size':
1903                    raise NotImplementedError()
1904                elif key == 'Use':
1905                    for h in histograms:
1906                        h['Use'] = False
1907                elif key == 'Scale':
1908                    for h in histograms:
1909                        h['Scale'][1] = False
1910
1911
1912##########################
1913# Command Line Interface #
1914##########################
1915# Each of these takes an argparse.Namespace object as their argument,
1916# representing the parsed command-line arguments for the relevant subcommand.
1917# The argument specification for each is in the subcommands dictionary (see
1918# below)
1919
1920
1921def create(args):
1922    """The create subcommand."""
1923    proj = G2Project(filename=args.filename)
1924
1925    hist_objs = []
1926    for h in args.histograms:
1927        hist_objs.append(proj.add_powder_histogram(h, args.iparams))
1928
1929    for p in args.phases:
1930        proj.add_phase(p, histograms=hist_objs)
1931
1932    proj.save()
1933
1934
1935def dump(args):
1936    """The dump subcommand"""
1937    if not args.histograms and not args.phases:
1938        args.raw = True
1939    if args.raw:
1940        import IPython.lib.pretty as pretty
1941
1942    for fname in args.files:
1943        if args.raw:
1944            proj, nameList = LoadDictFromProjFile(fname)
1945            print("file:", fname)
1946            print("names:", nameList)
1947            for key, val in proj.items():
1948                print(key, ":")
1949                pretty.pprint(val)
1950        else:
1951            proj = G2Project(fname)
1952            if args.histograms:
1953                hists = proj.histograms()
1954                for h in hists:
1955                    print(fname, "hist", h.id, h.name)
1956            if args.phases:
1957                phase_list = proj.phases()
1958                for p in phase_list:
1959                    print(fname, "phase", p.id, p.name)
1960
1961
1962def IPyBrowse(args):
1963    """Load a .gpx file and then open a IPython shell to browse it
1964    """
1965    for fname in args.files:
1966        proj, nameList = LoadDictFromProjFile(fname)
1967        msg = "\nfname {} loaded into proj (dict) with names in nameList".format(fname)
1968        GSASIIpath.IPyBreak_base(msg)
1969        break
1970
1971
1972def refine(args):
1973    """The refine subcommand"""
1974    proj = G2Project(args.gpxfile)
1975    if args.refinements is None:
1976        proj.refine()
1977    else:
1978        import json
1979        with open(args.refinements) as refs:
1980            refs = json.load(refs)
1981        proj.do_refinements(refs['refinements'])
1982
1983
1984def seqrefine(args):
1985    """The seqrefine subcommand"""
1986    raise NotImplementedError("seqrefine is not yet implemented")
1987
1988
1989def export(args):
1990    """The export subcommand"""
1991    # Export phase as CIF to args.exportfile
1992    proj = G2Project(args.gpxfile)
1993    phase = proj.phase(args.phase)
1994    phase.export_CIF(args.exportfile)
1995
1996
1997def _args_kwargs(*args, **kwargs):
1998    return args, kwargs
1999
2000# A dictionary of the name of each subcommand, and a tuple
2001# of its associated function and a list of its arguments
2002# The arguments are passed directly to the add_argument() method
2003# of an argparse.ArgumentParser
2004subcommands = {"create":
2005               (create, [_args_kwargs('filename',
2006                                      help='the project file to create. should end in .gpx'),
2007
2008                         _args_kwargs('-i', '--iparams',
2009                                      help='instrument parameter file'),
2010
2011                         _args_kwargs('-d', '--histograms',
2012                                      nargs='+',
2013                                      help='list of datafiles to add as histograms'),
2014
2015                         _args_kwargs('-p', '--phases',
2016                                      nargs='+',
2017                                      help='list of phases to add. phases are '
2018                                           'automatically associated with all '
2019                                           'histograms given.')]),
2020
2021               "dump": (dump, [_args_kwargs('-d', '--histograms',
2022                                     action='store_true',
2023                                     help='list histograms in files, overrides --raw'),
2024
2025                               _args_kwargs('-p', '--phases',
2026                                            action='store_true',
2027                                            help='list phases in files, overrides --raw'),
2028
2029                               _args_kwargs('-r', '--raw',
2030                                      action='store_true', help='dump raw file contents, default'),
2031
2032                               _args_kwargs('files', nargs='+')]),
2033
2034               "refine":
2035               (refine, [_args_kwargs('gpxfile', help='the project file to refine'),
2036                         _args_kwargs('refinements',
2037                                      help='json file of refinements to apply. if not present'
2038                                           ' refines file as-is',
2039                                      default=None,
2040                                      nargs='?')]),
2041
2042               "seqrefine": (seqrefine, []),
2043               "export": (export, [_args_kwargs('gpxfile',
2044                                                help='the project file from which to export'),
2045                                   _args_kwargs('phase', help='identifier of phase to export'),
2046                                   _args_kwargs('exportfile', help='the .cif file to export to')]),
2047               "browse": (IPyBrowse, [_args_kwargs('files', nargs='+',
2048                                                   help='list of files to browse')])}
2049
2050
2051def main():
2052    '''The command line interface for GSASIIscriptable.
2053
2054    Executes one of the following subcommands:
2055
2056        * :func:`create`
2057        * :func:`dump`
2058        * :func:`refine`
2059        * :func:`seqrefine`
2060        * :func:`export`
2061        * :func:`browse`
2062
2063    .. seealso::
2064        :func:`create`
2065        :func:`dump`
2066        :func:`refine`
2067        :func:`seqrefine`
2068        :func:`export`
2069        :func:`browse`
2070    '''
2071    parser = argparse.ArgumentParser()
2072    subs = parser.add_subparsers()
2073
2074    # Create all of the specified subparsers
2075    for name, (func, args) in subcommands.items():
2076        new_parser = subs.add_parser(name)
2077        for listargs, kwds in args:
2078            new_parser.add_argument(*listargs, **kwds)
2079        new_parser.set_defaults(func=func)
2080
2081    # Parse and trigger subcommand
2082    result = parser.parse_args()
2083    result.func(result)
2084
2085    # argv = sys.argv
2086    # if len(argv) > 1 and argv[1] in subcommands:
2087    #     subcommands[argv[1]](*argv[2:])
2088    # elif len(argv) == 1 or argv[1] in ('help', '--help', '-h'):
2089    #     # TODO print usage
2090    #     subcommand_names = ' | '.join(sorted(subcommands.keys()))
2091    #     print("USAGE: {} [ {} ] ...".format(argv[0], subcommand_names))
2092    # else:
2093    #     print("Unknown subcommand: {}".format(argv[1]))
2094    #     print("Available subcommands:")
2095    #     for name in sorted(subcommands.keys()):
2096    #         print("\t{}".format(name))
2097    #     sys.exit(-1)
2098    # sys.exit(0)
2099
2100if __name__ == '__main__':
2101    main()
2102
2103
2104    # from refinements.py
2105#      USAGE = """USAGE: {} datafile instparams phasefile projectname refinements
2106
2107# datafile:    Input powder data
2108# intparams:   Corresponding instrument parameter file
2109# phasefile:   Phase to refine against data
2110# projectname: Project file to be created, should end in .gpx
2111# refinements: JSON file of refinements to be executed
2112# """
2113#     try:
2114#         datafile, instprm, phasefile, projectname, refinements = sys.argv[1:]
2115#     except ValueError:
2116#         print(USAGE.format(sys.argv[0]))
2117#         sys.exit(-1)
2118
2119#     try:
2120#         with open(refinements) as f:
2121#             refinements = json.load(f)
2122#     except IOError:
2123#         print("No such refinements file: {}".format(refinements))
2124
2125#     print("Creating project file \"{}\"...".format(projectname))
2126#     proj = G2Project(filename=projectname)
2127#     # Add the histogram
2128#     hist = proj.add_powder_histogram(datafile, instprm)
2129#     # Add the phase, and associate it with the histogram
2130#     proj.add_phase(phasefile, histograms=[hist.name])
2131
2132#     proj.do_refinements(refinements['refinements'])
2133#     proj.save()
2134
2135
2136    # from gpx_dumper
2137    # import IPython.lib.pretty as pretty
2138    # proj, nameList = LoadDictFromProjFile(sys.argv[1])
2139    # print("names:", nameList)
2140    # for key, val in proj.items():
2141    #     print(key, ":", sep='')
2142    #     pretty.pprint(val)
Note: See TracBrowser for help on using the repository browser.