source: trunk/GSASIIscriptable.py @ 3026

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

add uiso setter

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