source: trunk/GSASIIscriptable.py @ 3015

Last change on this file since 3015 was 3015, checked in by odonnell, 4 years ago

refactor do_refinements

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