source: trunk/GSASIIscriptable.py @ 3009

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

more documentation changes, removed (phase|histogram)_names, add occupancy

File size: 70.3 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 os.path as ospath
35import sys
36import cPickle
37import imp
38import copy
39import os
40import random as ran
41
42import numpy.ma as ma
43import scipy.interpolate as si
44import numpy as np
45import scipy as sp
46
47import GSASIIpath
48GSASIIpath.SetBinaryPath(False) # would rather have this in __name__ == '__main__' stanza
49import GSASIIobj as G2obj
50import GSASIIpwd as G2pwd
51import GSASIIstrMain as G2strMain
52import GSASIIspc as G2spc
53import GSASIIElem as G2elem
54
55# Delay imports to not slow down small scripts
56G2fil = None
57PwdrDataReaders = []
58PhaseReaders = []
59
60
61def LoadG2fil():
62    """Delay importing this module, it is slow"""
63    global G2fil
64    global PwdrDataReaders
65    global PhaseReaders
66    if G2fil is None:
67        import GSASIIfiles
68        G2fil = GSASIIfiles
69        PwdrDataReaders = G2fil.LoadImportRoutines("pwd", "Powder_Data")
70        PhaseReaders = G2fil.LoadImportRoutines("phase", "Phase")
71
72
73def LoadDictFromProjFile(ProjFile):
74    '''Read a GSAS-II project file and load items to dictionary
75    :param str ProjFile: GSAS-II project (name.gpx) full file name
76    :returns: Project,nameList, where
77
78      * Project (dict) is a representation of gpx file following the GSAS-II tree struture
79        for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict
80        data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)}
81      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
82
83    Example for fap.gpx::
84
85      Project = {                 #NB:dict order is not tree order
86        u'Phases':{'data':None,'fap':{phase dict}},
87        u'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc},
88        u'Rigid bodies':{'data': {rigid body dict}},
89        u'Covariance':{'data':{covariance data dict}},
90        u'Controls':{'data':{controls data dict}},
91        u'Notebook':{'data':[notebook list]},
92        u'Restraints':{'data':{restraint data dict}},
93        u'Constraints':{'data':{constraint data dict}}]
94        }
95      nameList = [                #NB: reproduces tree order
96        [u'Notebook',],
97        [u'Controls',],
98        [u'Covariance',],
99        [u'Constraints',],
100        [u'Restraints',],
101        [u'Rigid bodies',],
102        [u'PWDR FAP.XRA Bank 1',
103             u'Comments',
104             u'Limits',
105             u'Background',
106             u'Instrument Parameters',
107             u'Sample Parameters',
108             u'Peak List',
109             u'Index Peak List',
110             u'Unit Cells List',
111             u'Reflection Lists'],
112        [u'Phases', u'fap']
113        ]
114    '''
115    # Let IOError be thrown if file does not exist
116    # if not ospath.exists(ProjFile):
117    #     print ('\n*** Error attempt to open project file that does not exist:\n   '+
118    #         str(ProjFile))
119    #     return
120    file = open(ProjFile,'rb')
121    # print('loading from file: {}'.format(ProjFile))
122    Project = {}
123    nameList = []
124    try:
125        while True:
126            try:
127                data = cPickle.load(file)
128            except EOFError:
129                break
130            datum = data[0]
131            Project[datum[0]] = {'data':datum[1]}
132            nameList.append([datum[0],])
133            for datus in data[1:]:
134                Project[datum[0]][datus[0]] = datus[1]
135                nameList[-1].append(datus[0])
136        file.close()
137        # print('project load successful')
138    except:
139        raise IOError("Error reading file "+str(ProjFile)+". This is not a GSAS-II .gpx file")
140    finally:
141        file.close()
142    return Project,nameList
143
144def SaveDictToProjFile(Project,nameList,ProjFile):
145    '''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile
146
147    :param dict Project: representation of gpx file following the GSAS-II
148        tree structure as described for LoadDictFromProjFile
149    :param list nameList: names of main tree entries & subentries used to reconstruct project file
150    :param str ProjFile: full file name for output project.gpx file (including extension)
151    '''
152    file = open(ProjFile,'wb')
153    # print('save to file: {}'.format(ProjFile))
154    try:
155        for name in nameList:
156            data = []
157            item = Project[name[0]]
158            data.append([name[0],item['data']])
159            for item2 in name[1:]:
160                data.append([item2,item[item2]])
161            cPickle.dump(data,file,1)
162    finally:
163        file.close()
164    # print('project save successful')
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
379def make_empty_project(author=None, filename=None):
380    """Creates an dictionary in the style of GSASIIscriptable, for an empty
381    project.
382
383    If no author name or filename is supplied, 'no name' and
384    <current dir>/test_output.gpx are used , respectively.
385
386    Returns: project dictionary, name list
387
388    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
389    """
390    if not filename:
391        filename = os.path.join(os.getcwd(), 'test_output.gpx')
392    else:
393        filename = os.path.abspath(filename)
394    gsasii_version = str(GSASIIpath.GetVersionNumber())
395    LoadG2fil()
396    import matplotlib as mpl
397    python_library_versions = G2fil.get_python_versions([mpl, np, sp])
398
399    controls_data = dict(G2obj.DefaultControls)
400    controls_data['LastSavedAs'] = unicode(filename)
401    controls_data['LastSavedUsing'] = gsasii_version
402    controls_data['PythonVersions'] = python_library_versions
403    if author:
404        controls_data['Author'] = author
405
406    output = {'Constraints': {'data': {'HAP': [], 'Hist': [], 'Phase': [],
407                                       'Global': []}},
408              'Controls': {'data': controls_data},
409              u'Covariance': {'data': {}},
410              u'Notebook': {'data': ['']},
411              u'Restraints': {'data': {}},
412              u'Rigid bodies': {'data': {'RBIds': {'Residue': [], 'Vector': []},
413                                'Residue': {'AtInfo': {}},
414                                'Vector':  {'AtInfo': {}}}}}
415
416    names = [[u'Notebook'], [u'Controls'], [u'Covariance'],
417             [u'Constraints'], [u'Restraints'], [u'Rigid bodies']]
418
419    return output, names
420
421
422class G2ImportException(Exception):
423    pass
424
425
426def import_generic(filename, readerlist):
427    """Attempt to import a filename, using a list of reader objects.
428
429    Returns the first reader object which worked."""
430    # Translated from OnImportGeneric method in GSASII.py
431    primaryReaders, secondaryReaders = [], []
432    for reader in readerlist:
433        flag = reader.ExtensionValidator(filename)
434        if flag is None:
435            secondaryReaders.append(reader)
436        elif flag:
437            primaryReaders.append(reader)
438    if not secondaryReaders and not primaryReaders:
439        raise G2ImportException("Could not read file: ", filename)
440
441    with open(filename, 'Ur') as fp:
442        rd_list = []
443
444        for rd in primaryReaders + secondaryReaders:
445            # Initialize reader
446            rd.selections = []
447            rd.dnames = []
448            rd.ReInitialize()
449            # Rewind file
450            fp.seek(0)
451            rd.errors = ""
452            if not rd.ContentsValidator(fp):
453                # Report error
454                pass
455            if len(rd.selections) > 1:
456                # Select data?
457                # GSASII.py:543
458                raise G2ImportException("Not sure what data to select")
459
460            block = 0
461            rdbuffer = {}
462            fp.seek(0)
463            repeat = True
464            while repeat:
465                repeat = False
466                block += 1
467                rd.objname = os.path.basename(filename)
468                flag = rd.Reader(filename, fp, buffer=rdbuffer, blocknum=block)
469                if flag:
470                    # Omitting image loading special cases
471                    rd.readfilename = filename
472                    rd_list.append(copy.deepcopy(rd))
473                    repeat = rd.repeat
474                else:
475                    raise G2ImportException("{}.Reader() returned:".format(rd),
476                                            flag)
477
478            if rd_list:
479                if rd.warnings:
480                    print("Read warning by", rd.formatName, "reader:",
481                          rd.warnings, file=sys.stderr)
482                return rd_list
483    raise G2ImportException("No reader could read file: " + filename)
484
485
486def load_iprms(instfile, reader):
487    """Loads instrument parameters from a file, and edits the
488    given reader.
489
490    Returns a 2-tuple of (Iparm1, Iparm2) parameters
491    """
492    LoadG2fil()
493    ext = os.path.splitext(instfile)[1]
494
495    if ext.lower() == '.instprm':
496        # New GSAS File, load appropriately
497        with open(instfile) as f:
498            lines = f.readlines()
499        bank = reader.powderentry[2]
500        numbanks = reader.numbanks
501        iparms = G2fil.ReadPowderInstprm(lines, bank, numbanks, reader)
502
503        reader.instfile = instfile
504        reader.instmsg = 'GSAS-II file' + instfile
505        return iparms
506    elif ext.lower() not in ('.prm', '.inst', '.ins'):
507        raise ValueError('Expected .prm file, found: ', instfile)
508
509    # It's an old GSAS file, load appropriately
510    Iparm = {}
511    with open(instfile, 'Ur') as fp:
512        for line in fp:
513            if '#' in line:
514                continue
515            Iparm[line[:12]] = line[12:-1]
516    ibanks = int(Iparm.get('INS   BANK  ', '1').strip())
517    if ibanks == 1:
518        reader.instbank = 1
519        reader.powderentry[2] = 1
520        reader.instfile = instfile
521        reader.instmsg = instfile + ' bank ' + str(reader.instbank)
522        return G2fil.SetPowderInstParms(Iparm, reader)
523    # TODO handle >1 banks
524    raise NotImplementedError("Check GSASIIfiles.py: ReadPowderInstprm")
525
526def load_pwd_from_reader(reader, instprm, existingnames=[]):
527    """Loads powder data from a reader object, and assembles it into a GSASII data tree.
528
529    :returns: (name, tree) - 2-tuple of the histogram name (str), and data
530
531    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
532    """
533    HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_')
534    HistName = unicode(G2obj.MakeUniqueLabel(HistName, existingnames))
535
536    try:
537        Iparm1, Iparm2 = instprm
538    except ValueError:
539        Iparm1, Iparm2 = load_iprms(instprm, reader)
540
541    Ymin = np.min(reader.powderdata[1])
542    Ymax = np.max(reader.powderdata[1])
543    valuesdict = {'wtFactor': 1.0,
544                  'Dummy': False,
545                  'ranId': ran.randint(0, sys.maxint),
546                  'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax,
547                  'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax,
548                  'qPlot': False, 'dPlot': False, 'sqrtPlot': False,
549                  'Yminmax': [Ymin, Ymax]}
550    reader.Sample['ranId'] = valuesdict['ranId']
551
552    # Ending keys:
553    # [u'Reflection Lists',
554    #  u'Limits',
555    #  'data',
556    #  u'Index Peak List',
557    #  u'Comments',
558    #  u'Unit Cells List',
559    #  u'Sample Parameters',
560    #  u'Peak List',
561    #  u'Background',
562    #  u'Instrument Parameters']
563    Tmin = np.min(reader.powderdata[0])
564    Tmax = np.max(reader.powderdata[0])
565
566    default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0],
567                          {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}]
568
569    output_dict = {u'Reflection Lists': {},
570                   u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]),
571                   u'data': [valuesdict, reader.powderdata, HistName],
572                   u'Index Peak List': [[], []],
573                   u'Comments': reader.comments,
574                   u'Unit Cells List': [],
575                   u'Sample Parameters': reader.Sample,
576                   u'Peak List': {'peaks': [], 'sigDict': {}},
577                   u'Background': reader.pwdparms.get('Background', default_background),
578                   u'Instrument Parameters': [Iparm1, Iparm2],
579                   }
580
581    names = [u'Comments',
582             u'Limits',
583             u'Background',
584             u'Instrument Parameters',
585             u'Sample Parameters',
586             u'Peak List',
587             u'Index Peak List',
588             u'Unit Cells List',
589             u'Reflection Lists']
590
591    # TODO controls?? GSASII.py:1664-7
592
593    return HistName, [HistName] + names, output_dict
594
595
596def _deep_copy_into(from_, into):
597    """Helper function for reloading .gpx file. See G2Project.reload()
598
599    :author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
600    """
601    if isinstance(from_, dict) and isinstance(into, dict):
602        combined_keys = set(from_.keys()).union(into.keys())
603        for key in combined_keys:
604            if key in from_ and key in into:
605                both_dicts = (isinstance(from_[key], dict)
606                              and isinstance(into[key], dict))
607                both_lists = (isinstance(from_[key], list)
608                              and isinstance(into[key], list))
609                if both_dicts or both_lists:
610                    _deep_copy_into(from_[key], into[key])
611                else:
612                    into[key] = from_[key]
613            elif key in from_:
614                into[key] = from_[key]
615            else:  # key in into
616                del into[key]
617    elif isinstance(from_, list) and isinstance(into, list):
618        if len(from_) == len(into):
619            for i in xrange(len(from_)):
620                both_dicts = (isinstance(from_[i], dict)
621                              and isinstance(into[i], dict))
622                both_lists = (isinstance(from_[i], list)
623                              and isinstance(into[i], list))
624                if both_dicts or both_lists:
625                    _deep_copy_into(from_[i], into[i])
626                else:
627                    into[i] = from_[i]
628        else:
629            into[:] = from_
630
631
632class G2ObjectWrapper(object):
633    """Base class for all GSAS-II object wrappers.
634
635    The underlying GSAS-II format can be accessed as `wrapper.data`. A number
636    of overrides are implemented so that the wrapper behaves like a dictionary.
637
638    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
639    """
640    def __init__(self, datadict):
641        self.data = datadict
642
643    def __getitem__(self, key):
644        return self.data[key]
645
646    def __setitem__(self, key, value):
647        self.data[key] = value
648
649    def __contains__(self, key):
650        return key in self.data
651
652    def get(self, k, d=None):
653        return self.data.get(k, d)
654
655    def keys(self):
656        return self.data.keys()
657
658    def values(self):
659        return self.data.values()
660
661    def items(self):
662        return self.data.items()
663
664
665class G2Project(G2ObjectWrapper):
666    """Represents an entire GSAS-II project."""
667    def __init__(self, gpxfile=None, author=None, filename=None):
668        """Loads a GSAS-II project from a specified filename.
669
670        :param str gpxfile: Existing .gpx file to be loaded. If nonexistent,
671            creates an empty project.
672        :param str author: Author's name. Optional.
673        :param str filename: The filename the project should be saved to in
674            the future. If both filename and gpxfile are present, the project is
675            loaded from the gpxfile, then set to  save to filename in the future"""
676        if gpxfile is None:
677            filename = os.path.abspath(os.path.expanduser(filename))
678            self.filename = filename
679            self.data, self.names = make_empty_project(author=author, filename=filename)
680        elif isinstance(gpxfile, str):
681            # TODO set author, filename
682            self.filename = os.path.abspath(os.path.expanduser(gpxfile))
683            self.data, self.names = LoadDictFromProjFile(gpxfile)
684            self.index_ids()
685        else:
686            raise ValueError("Not sure what to do with gpxfile")
687
688    @classmethod
689    def from_dict_and_names(cls, gpxdict, names, filename=None):
690        """Creates a :class:`G2Project` directly from
691        a dictionary and a list of names. If in doubt, do not use this.
692
693        :returns: a :class:`G2Project`
694        """
695        out = cls()
696        if filename:
697            filename = os.path.abspath(os.path.expanduser(filename))
698            out.filename = filename
699            gpxdict['Controls']['data']['LastSavedAs'] = filename
700        else:
701            try:
702                out.filename = gpxdict['Controls']['data']['LastSavedAs']
703            except KeyError:
704                out.filename = None
705        out.data = gpxdict
706        out.names = names
707
708    def save(self, filename=None):
709        """Saves the project, either to the current filename, or to a new file.
710
711        Updates self.filename if a new filename provided"""
712        # TODO update LastSavedUsing ?
713        if filename:
714            filename = os.path.abspath(os.path.expanduser(filename))
715            self.data['Controls']['data']['LastSavedAs'] = filename
716            self.filename = filename
717        elif not self.filename:
718            raise AttributeError("No file name to save to")
719        SaveDictToProjFile(self.data, self.names, self.filename)
720
721    def add_powder_histogram(self, datafile, iparams):
722        """Loads a powder data histogram into the project.
723
724        Automatically checks for an instrument parameter file, or one can be
725        provided.
726
727        :param str datafile: The powder data file to read, a filename.
728        :param str iparams: The instrument parameters file, a filename.
729
730        :returns: A :class:`G2PwdrData` object representing
731            the histogram
732        """
733        LoadG2fil()
734        datafile = os.path.abspath(os.path.expanduser(datafile))
735        iparams = os.path.abspath(os.path.expanduser(iparams))
736        pwdrreaders = import_generic(datafile, PwdrDataReaders)
737        histname, new_names, pwdrdata = load_pwd_from_reader(
738                                          pwdrreaders[0], iparams,
739                                          [h.name for h in self.histograms()])
740        if histname in self.data:
741            print("Warning - redefining histogram", histname)
742        else:
743            if self.names[-1][0] == 'Phases':
744                self.names.insert(-1, new_names)
745            else:
746                self.names.append(new_names)
747        self.data[histname] = pwdrdata
748        return self.histogram(histname)
749
750    def add_phase(self, phasefile, phasename=None, histograms=[]):
751        """Loads a phase into the project from a .cif file
752
753        :param str phasefile: The CIF file from which to import the phase.
754        :param str phasename: The name of the new phase, or None for the default
755        :param list histograms: The names of the histograms to associate with
756            this phase
757
758        :returns: A :class:`G2Phase` object representing the
759            new phase.
760        """
761        LoadG2fil()
762        phasefile = os.path.abspath(os.path.expanduser(phasefile))
763
764        # TODO handle multiple phases in a file
765        phasereaders = import_generic(phasefile, PhaseReaders)
766        phasereader = phasereaders[0]
767
768        phasename = phasename or phasereader.Phase['General']['Name']
769        phaseNameList = [p.name for p in self.phases()]
770        phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
771        phasereader.Phase['General']['Name'] = phasename
772
773        if 'Phases' not in self.data:
774            self.data[u'Phases'] = { 'data': None }
775        assert phasename not in self.data['Phases'], "phase names should be unique"
776        self.data['Phases'][phasename] = phasereader.Phase
777
778        if phasereader.Constraints:
779            Constraints = self.data['Constraints']
780            for i in phasereader.Constraints:
781                if isinstance(i, dict):
782                    if '_Explain' not in Constraints:
783                        Constraints['_Explain'] = {}
784                    Constraints['_Explain'].update(i)
785                else:
786                    Constraints['Phase'].append(i)
787
788        data = self.data['Phases'][phasename]
789        generalData = data['General']
790        SGData = generalData['SGData']
791        NShkl = len(G2spc.MustrainNames(SGData))
792        NDij = len(G2spc.HStrainNames(SGData))
793        Super = generalData.get('Super', 0)
794        if Super:
795            SuperVec = np.array(generalData['SuperVec'][0])
796        else:
797            SuperVec = []
798        UseList = data['Histograms']
799
800        for histname in histograms:
801            if histname.startswith('HKLF '):
802                raise NotImplementedError("Does not support HKLF yet")
803            elif histname.startswith('PWDR '):
804                hist = self.histogram(histname)
805                hist['Reflection Lists'][generalData['Name']] = {}
806                UseList[histname] = SetDefaultDData('PWDR', histname, NShkl=NShkl, NDij=NDij)
807                for key, val in [('Use', True), ('LeBail', False),
808                                 ('newLeBail', True),
809                                 ('Babinet', {'BabA': [0.0, False],
810                                              'BabU': [0.0, False]})]:
811                    if key not in UseList[histname]:
812                        UseList[histname][key] = val
813            else:
814                raise NotImplementedError("Unexpected histogram" + histname)
815
816        for obj in self.names:
817            if obj[0] == 'Phases':
818                phasenames = obj
819                break
820        else:
821            phasenames = [u'Phases']
822            self.names.append(phasenames)
823        phasenames.append(unicode(phasename))
824
825        # TODO should it be self.filename, not phasefile?
826        SetupGeneral(data, os.path.dirname(phasefile))
827        self.index_ids()
828
829        return self.phase(phasename)
830
831    def reload(self):
832        """Reload self from self.filename"""
833        data, names = LoadDictFromProjFile(self.filename)
834        self.names = names
835        # Need to deep copy the new data file data into the current tree,
836        # so that any existing G2Phase, or G2PwdrData objects will still be
837        # valid
838        _deep_copy_into(from_=data, into=self.data)
839
840    def refine(self, newfile=None, printFile=None, makeBack=False):
841        # index_ids will automatically save the project
842        self.index_ids()
843        # TODO G2strMain does not properly use printFile
844        G2strMain.Refine(self.filename, makeBack=makeBack)
845        # Reload yourself
846        self.reload()
847
848    def histogram(self, histname):
849        """Returns the histogram named histname, or None if it does not exist.
850
851        :param histname: The name of the histogram (str), or ranId or index.
852        :returns: A :class:`G2PwdrData` object, or None if
853            the histogram does not exist
854
855        .. seealso::
856            :func:`~GSASIIscriptable.G2Project.histograms`
857            :func:`~GSASIIscriptable.G2Project.phase`
858            :func:`~GSASIIscriptable.G2Project.phases`
859            """
860        if histname in self.data:
861            return G2PwdrData(self.data[histname], self)
862        for key, val in G2obj.HistIdLookup.items():
863            name, ranId = val
864            # histname can be either ranId (key) or index (val)
865            if ranId == histname or key == str(histname):
866                return self.histogram(name)
867
868    def histograms(self):
869        """Return a list of all histograms, as
870        :class:`G2PwdrData` objects
871
872        .. seealso::
873            :func:`~GSASIIscriptable.G2Project.histograms`
874            :func:`~GSASIIscriptable.G2Project.phase`
875            :func:`~GSASIIscriptable.G2Project.phases`
876            """
877        output = []
878        for obj in self.names:
879            if len(obj) > 1 and obj[0] != u'Phases':
880                output.append(self.histogram(obj[0]))
881        return output
882
883    def phase(self, phasename):
884        """
885        Gives an object representing the specified phase in this project.
886
887        :param str phasename: The name of the desired phase. Either the name
888            (str), the phase's ranId, or the phase's index
889        :returns: A :class:`G2Phase` object
890        :raises: KeyError
891
892        .. seealso::
893            :func:`~GSASIIscriptable.G2Project.histograms`
894            :func:`~GSASIIscriptable.G2Project.phase`
895            :func:`~GSASIIscriptable.G2Project.phases`
896            """
897        phases = self.data['Phases']
898        if phasename in phases:
899            return G2Phase(phases[phasename], phasename, self)
900        for key, val in G2obj.PhaseIdLookup.items():
901            name, ranId = val
902            # phasename can be either ranId (key) or index (val)
903            if ranId == phasename or key == str(phasename):
904                return self.phase(name)
905
906    def phases(self):
907        """
908        Returns a list of all the phases in the project.
909
910        :returns: A list of :class:`G2Phase` objects
911
912        .. seealso::
913            :func:`~GSASIIscriptable.G2Project.histogram`
914            :func:`~GSASIIscriptable.G2Project.histograms`
915            :func:`~GSASIIscriptable.G2Project.phase`
916            """
917        for obj in self.names:
918            if obj[0] == 'Phases':
919                return [self.phase(p) for p in obj[1:]]
920        return []
921
922    def do_refinements(self, refinements, histogram='all', phase='all',
923                       outputnames=None):
924        """Conducts a series of refinements.
925
926        :param list refinements: A list of dictionaries defining refinements
927        :param str histogram: Name of histogram for refinements to be applied
928            to, or 'all'
929        :param str phase: Name of phase for refinements to be applied to, or
930            'all'
931        """
932        if outputnames:
933            if len(refinements) != len(outputnames):
934                raise ValueError("Should have same number of outuputs to"
935                                 "refinements")
936        else:
937            outputnames = [None for r in refinements]
938
939        for output, refinement in zip(outputnames, refinements):
940            self.set_refinement(refinement, histogram)
941            # Handle 'once' args - refinements that are disabled after this
942            # refinement
943            if 'once' in refinement:
944                temp = {'set': refinement['once']}
945                self.set_refinement(temp, histogram, phase)
946
947            if output:
948                self.save(output)
949
950            self.refine()  # newFile=output)
951
952            # Handle 'once' args - refinements that are disabled after this
953            # refinement
954            if 'once' in refinement:
955                temp = {'clear': refinement['once']}
956                self.set_refinement(temp, histogram, phase)
957
958    def set_refinement(self, refinement, histogram='all', phase='all'):
959        """Apply specified refinements to a given histogram(s) or phase(s).
960
961        Refinement parameters are categorize in three groups:
962
963        1. Histogram parameters
964        2. Phase parameters
965        3. Histogram-and-Phase (HAP) parameters
966
967        :param dict refinement: The refinements to be conducted
968        :param histogram: Either a name of a histogram (str), a list of
969            histogram names, or 'all' (default)
970        :param phase: Either a name of a phase (str), a list of phase names, or
971            'all' (default)
972
973        .. seealso::
974            :func:`~G2PwdrData.set_refinements`
975            :func:`~G2PwdrData.clear_refinements`
976            :func:`~G2Phase.set_refinements`
977            :func:`~G2Phase.clear_refinements`
978            :func:`~G2Phase.set_HAP_refinements`
979            :func:`~G2Phase.clear_HAP_refinements`"""
980
981        if histogram == 'all':
982            hists = self.histograms()
983        elif isinstance(histogram, str) or isinstance(histogram, int):
984            hists = [self.histogram(histogram)]
985        else:
986            hists = [self.histogram(name) for name in histogram]
987
988        if phase == 'all':
989            phases = self.phases()
990        elif isinstance(phase, str) or isinstance(phase, int):
991            phases = [self.phase(phase)]
992        else:
993            phases = [self.phase(name) for name in phase]
994
995
996        # TODO: HAP parameters:
997        #   Babinet
998        #   Extinction
999        #   HStrain
1000        #   Mustrain
1001        #   Pref. Ori
1002        #   Size
1003
1004        pwdr_set = {}
1005        phase_set = {}
1006        hap_set = {}
1007        for key, val in refinement.get('set', {}).items():
1008            # Apply refinement options
1009            if G2PwdrData.is_valid_refinement_key(key):
1010                pwdr_set[key] = val
1011            elif G2Phase.is_valid_refinement_key(key):
1012                phase_set[key] = val
1013            elif G2Phase.is_valid_HAP_refinement_key(key):
1014                hap_set[key] = val
1015            else:
1016                raise ValueError("Unknown refinement key", key)
1017
1018        for hist in hists:
1019            hist.set_refinements(pwdr_set)
1020        for phase in phases:
1021            phase.set_refinements(phase_set)
1022        for phase in phases:
1023            phase.set_HAP_refinements(hap_set, hists)
1024
1025        pwdr_clear = {}
1026        phase_clear = {}
1027        hap_clear = {}
1028        for key, val in refinement.get('clear', {}).items():
1029            # Clear refinement options
1030            if G2PwdrData.is_valid_refinement_key(key):
1031                pwdr_clear[key] = val
1032            elif G2Phase.is_valid_refinement_key(key):
1033                phase_clear[key] = val
1034            elif G2Phase.is_valid_HAP_refinement_key(key):
1035                hap_set[key] = val
1036            else:
1037                raise ValueError("Unknown refinement key", key)
1038
1039        for hist in hists:
1040            hist.clear_refinements(pwdr_clear)
1041        for phase in phases:
1042            phase.clear_refinements(phase_clear)
1043        for phase in phases:
1044            phase.clear_HAP_refinements(hap_clear, hists)
1045
1046    def index_ids(self):
1047        import GSASIIstrIO as G2strIO
1048        self.save()
1049        return G2strIO.GetUsedHistogramsAndPhases(self.filename)
1050
1051    def add_constraint_raw(self, cons_scope, constr):
1052        """Adds a constraint of type consType to the project.
1053        cons_scope should be one of "Hist", "Phase", "HAP", or "Global".
1054
1055        WARNING it does not check the constraint is well-constructed"""
1056        constrs = self['Constraints']['data']
1057        if 'Global' not in constrs:
1058            constrs['Global'] = []
1059        constrs[cons_scope].append(constr)
1060
1061    def hold_many(self, vars, type):
1062        """Apply holds for all the variables in vars, for constraint of a given type.
1063
1064        type is passed directly to add_constraint_raw as consType
1065
1066        :param list vars: A list of variables to hold. Either :class:`GSASIIobj.G2VarObj` objects,
1067            string variable specifiers, or arguments for :meth:`make_var_obj`
1068        :param str type: A string constraint type specifier. See
1069        :class:`~GSASIIscriptable.G2Project.add_constraint_raw`"""
1070        for var in vars:
1071            if isinstance(var, str):
1072                var = self.make_var_obj(var)
1073            elif not isinstance(var, G2obj.G2VarObj):
1074                var = self.make_var_obj(*var)
1075            self.add_constraint_raw(type, [[1.0, var], None, None, 'h'])
1076
1077    def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
1078                     reloadIdx=True):
1079        """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a")
1080        or individual names of phase, histogram, varname, and atomId.
1081
1082        Automatically converts string phase, hist, or atom names into the ID required
1083        by G2VarObj."""
1084
1085        if reloadIdx:
1086            self.index_ids()
1087
1088        # If string representation, short circuit
1089        if hist is None and varname is None and atomId is None:
1090            if isinstance(phase, str) and ':' in phase:
1091                return G2obj.G2VarObj(phase)
1092
1093        # Get phase index
1094        phaseObj = None
1095        if isinstance(phase, G2Phase):
1096            phaseObj = phase
1097            phase = G2obj.PhaseRanIdLookup[phase.ranId]
1098        elif phase in self['Phases']:
1099            phaseObj = self.phase(phase)
1100            phaseRanId = phaseObj.ranId
1101            phase = G2obj.PhaseRanIdLookup[phaseRanId]
1102
1103        # Get histogram index
1104        if isinstance(hist, G2PwdrData):
1105            hist = G2obj.HistRanIdLookup[hist.ranId]
1106        elif hist in self.data:
1107            histRanId = self.histogram(hist).ranId
1108            hist = G2obj.HistRanIdLookup[histRanId]
1109
1110        # Get atom index (if any)
1111        if isinstance(atomId, G2AtomRecord):
1112            atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
1113        elif phaseObj:
1114            atomObj = phaseObj.atom(atomId)
1115            if atomObj:
1116                atomRanId = atomObj.ranId
1117                atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
1118
1119        return G2obj.G2VarObj(phase, hist, varname, atomId)
1120
1121
1122class G2AtomRecord(G2ObjectWrapper):
1123    """Wrapper for an atom record. Has convenient accessors via @property.
1124
1125
1126    Available accessors: label, type, refinement_flags, coordinates,
1127        occupancy, ranId, id, adp_flag, uiso
1128
1129    Example:
1130
1131    >>> atom = some_phase.atom("O3")
1132    >>> # We can access the underlying data format
1133    >>> atom.data
1134    ['O3', 'O-2', '', ... ]
1135    >>> # We can also use wrapper accessors
1136    >>> atom.coordinates
1137    (0.33, 0.15, 0.5)
1138    >>> atom.refinement_flags
1139    u'FX'
1140    >>> atom.ranId
1141    4615973324315876477
1142    >>> atom.occupancy
1143    1.0
1144    """
1145    def __init__(self, data, indices, proj):
1146        self.data = data
1147        self.cx, self.ct, self.cs, self.cia = indices
1148        self.proj = proj
1149
1150    @property
1151    def label(self):
1152        return self.data[self.ct-1]
1153
1154    @property
1155    def type(self):
1156        return self.data[self.ct]
1157
1158    @property
1159    def refinement_flags(self):
1160        return self.data[self.ct+1]
1161
1162    @refinement_flags.setter
1163    def refinement_flags(self, other):
1164        # Automatically check it is a valid refinement
1165        for c in other:
1166            if c not in ' FXU':
1167                raise ValueError("Invalid atom refinement: ", other)
1168        self.data[self.ct+1] = unicode(other)
1169
1170    @property
1171    def coordinates(self):
1172        return tuple(self.data[self.cx:self.cx+3])
1173
1174    @property
1175    def occupancy(self):
1176        return self.data[self.cx+3]
1177
1178    @occupancy.setter
1179    def occupancy(self, val):
1180        self.data[self.cx+3] = float(val)
1181
1182    @property
1183    def ranId(self):
1184        return self.data[self.cia+8]
1185
1186    @property
1187    def adp_flag(self):
1188        # Either 'I' or 'A'
1189        return self.data[self.cia]
1190
1191    @property
1192    def uiso(self):
1193        if self.adp_flag == 'I':
1194            return self.data[self.cia+1]
1195        else:
1196            return self.data[self.cia+2:self.cia+8]
1197
1198
1199class G2PwdrData(G2ObjectWrapper):
1200    """Wraps a Powder Data Histogram."""
1201    def __init__(self, data, proj):
1202        self.data = data
1203        self.proj = proj
1204
1205    @staticmethod
1206    def is_valid_refinement_key(key):
1207        valid_keys = ['Limits', 'Sample Parameters', 'Background',
1208                      'Instrument Parameters']
1209        return key in valid_keys
1210
1211    @property
1212    def name(self):
1213        return self['data'][-1]
1214
1215    @property
1216    def ranId(self):
1217        return self['data'][0]['ranId']
1218
1219    @property
1220    def residuals(self):
1221        data = self['data'][0]
1222        return {key: data[key]
1223                for key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']}
1224
1225    @property
1226    def id(self):
1227        return G2obj.HistRanIdLookup[self.ranId]
1228
1229    def fit_fixed_points(self):
1230        """Attempts to apply a background fit to the fixed points currently specified."""
1231        def SetInstParms(Inst):
1232            dataType = Inst['Type'][0]
1233            insVary = []
1234            insNames = []
1235            insVals = []
1236            for parm in Inst:
1237                insNames.append(parm)
1238                insVals.append(Inst[parm][1])
1239                if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1240                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1241                        Inst[parm][2] = False
1242            instDict = dict(zip(insNames, insVals))
1243            instDict['X'] = max(instDict['X'], 0.01)
1244            instDict['Y'] = max(instDict['Y'], 0.01)
1245            if 'SH/L' in instDict:
1246                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
1247            return dataType, instDict, insVary
1248
1249        bgrnd = self['Background']
1250
1251        # Need our fixed points in order
1252        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
1253        X = [x for x, y in bgrnd[1]['FixedPoints']]
1254        Y = [y for x, y in bgrnd[1]['FixedPoints']]
1255
1256        limits = self['Limits'][1]
1257        if X[0] > limits[0]:
1258            X = [limits[0]] + X
1259            Y = [Y[0]] + Y
1260        if X[-1] < limits[1]:
1261            X += [limits[1]]
1262            Y += [Y[-1]]
1263
1264        # Some simple lookups
1265        controls = self.proj['Controls']['data']
1266        inst, inst2 = self['Instrument Parameters']
1267        pwddata = self['data'][1]
1268
1269        # Construct the data for background fitting
1270        xBeg = np.searchsorted(pwddata[0], limits[0])
1271        xFin = np.searchsorted(pwddata[0], limits[1])
1272        xdata = pwddata[0][xBeg:xFin]
1273        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
1274
1275        W = [1]*len(xdata)
1276        Z = [0]*len(xdata)
1277
1278        dataType, insDict, insVary = SetInstParms(inst)
1279        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1280
1281        # Do the fit
1282        data = np.array([xdata, ydata, W, Z, Z, Z])
1283        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
1284                        prevVaryList=bakVary, controls=controls)
1285
1286        # Post-fit
1287        parmDict = {}
1288        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1289        parmDict.update(bakDict)
1290        parmDict.update(insDict)
1291        pwddata[3][xBeg:xFin] *= 0
1292        pwddata[5][xBeg:xFin] *= 0
1293        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
1294
1295        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
1296        # TODO update background
1297        self.proj.save()
1298
1299    def y_calc(self):
1300        return self['data'][1][3]
1301
1302    def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
1303        try:
1304            import matplotlib.pyplot as plt
1305            data = self['data'][1]
1306            if Yobs:
1307                plt.plot(data[0], data[1], label='Yobs')
1308            if Ycalc:
1309                plt.plot(data[0], data[3], label='Ycalc')
1310            if Background:
1311                plt.plot(data[0], data[4], label='Background')
1312            if Residual:
1313                plt.plot(data[0], data[5], label="Residual")
1314        except ImportError:
1315            pass
1316
1317    def set_refinements(self, refs):
1318        """Sets the refinement parameter 'key' to the specification 'value'
1319
1320        :param dict refs: A dictionary of the parameters to be set.
1321
1322        :returns: None"""
1323        do_fit_fixed_points = False
1324        for key, value in refs.items():
1325            if key == 'Limits':
1326                old_limits = self.data['Limits'][1]
1327                new_limits = value
1328                if isinstance(new_limits, dict):
1329                    if 'low' in new_limits:
1330                        old_limits[0] = new_limits['low']
1331                    if 'high' in new_limits:
1332                        old_limits[1] = new_limits['high']
1333                else:
1334                    old_limits[0], old_limits[1] = new_limits
1335            elif key == 'Sample Parameters':
1336                sample = self.data['Sample Parameters']
1337                for sparam in value:
1338                    sample[sparam][1] = True
1339            elif key == 'Background':
1340                bkg, peaks = self.data['Background']
1341
1342                # If True or False, just set the refine parameter
1343                if value in (True, False):
1344                    bkg[1] = value
1345                    return
1346
1347                if 'type' in value:
1348                    bkg[0] = value['type']
1349                if 'refine' in value:
1350                    bkg[1] = value['refine']
1351                if 'no. coeffs' in value:
1352                    cur_coeffs = bkg[2]
1353                    n_coeffs = value['no. coeffs']
1354                    if n_coeffs > cur_coeffs:
1355                        for x in range(n_coeffs - cur_coeffs):
1356                            bkg.append(0.0)
1357                    else:
1358                        for _ in range(cur_coeffs - n_coeffs):
1359                            bkg.pop()
1360                    bkg[2] = n_coeffs
1361                if 'coeffs' in value:
1362                    bkg[3:] = value['coeffs']
1363                if 'FixedPoints' in value:
1364                    peaks['FixedPoints'] = [(float(a), float(b))
1365                                            for a, b in value['FixedPoints']]
1366                if value.get('fit fixed points', False):
1367                    do_fit_fixed_points = True
1368
1369            elif key == 'Instrument Parameters':
1370                instrument, secondary = self.data['Instrument Parameters']
1371                for iparam in value:
1372                    try:
1373                        instrument[iparam][2] = True
1374                    except IndexError:
1375                        raise ValueError("Invalid key:", iparam)
1376            else:
1377                raise ValueError("Unknown key:", key)
1378        # Fit fixed points after the fact - ensure they are after fixed points
1379        # are added
1380        if do_fit_fixed_points:
1381            # Background won't be fit if refinement flag not set
1382            orig = self['Background'][0][1]
1383            self['Background'][0][1] = True
1384            self.fit_fixed_points()
1385            # Restore the previous value
1386            self['Background'][0][1] = orig
1387
1388    def clear_refinements(self, refs):
1389        """Clears the refinement parameter 'key' and its associated value.
1390
1391        :param dict refs: A dictionary of parameters to clear."""
1392        for key, value in refs.items():
1393            if key == 'Limits':
1394                old_limits, cur_limits = self.data['Limits']
1395                cur_limits[0], cur_limits[1] = old_limits
1396            elif key == 'Sample Parameters':
1397                sample = self.data['Sample Parameters']
1398                for sparam in value:
1399                    sample[sparam][1] = False
1400            elif key == 'Background':
1401                bkg, peaks = self.data['Background']
1402
1403                # If True or False, just set the refine parameter
1404                if value in (True, False):
1405                    bkg[1] = False
1406                    return
1407
1408                bkg[1] = False
1409                if 'FixedPoints' in value:
1410                    if 'FixedPoints' in peaks:
1411                        del peaks['FixedPoints']
1412            elif key == 'Instrument Parameters':
1413                instrument, secondary = self.data['Instrument Parameters']
1414                for iparam in value:
1415                    instrument[iparam][2] = False
1416            else:
1417                raise ValueError("Unknown key:", key)
1418
1419
1420class G2Phase(G2ObjectWrapper):
1421    """A wrapper object around a given phase.
1422
1423    Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com)
1424    """
1425    def __init__(self, data, name, proj):
1426        self.data = data
1427        self.name = name
1428        self.proj = proj
1429
1430    @staticmethod
1431    def is_valid_refinement_key(key):
1432        valid_keys = ["Cell", "Atoms", "LeBail"]
1433        return key in valid_keys
1434
1435    @staticmethod
1436    def is_valid_HAP_refinement_key(key):
1437        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
1438                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
1439        return key in valid_keys
1440
1441    def atom(self, atomlabel):
1442        """Returns the atom specified by atomlabel, or None if it does not
1443        exist.
1444
1445        :param str atomlabel: The name of the atom (e.g. "O2")
1446        :returns: A :class:`G2AtomRecord` object
1447            representing the atom.
1448        """
1449        # Consult GSASIIobj.py for the meaning of this
1450        cx, ct, cs, cia = self.data['General']['AtomPtrs']
1451        ptrs = [cx, ct, cs, cia]
1452        atoms = self.data['Atoms']
1453        for atom in atoms:
1454            if atom[ct-1] == atomlabel:
1455                return G2AtomRecord(atom, ptrs, self.proj)
1456
1457    def atoms(self):
1458        """Returns a list of atoms present in the phase.
1459
1460        :returns: A list of :class:`G2AtomRecord` objects.
1461
1462        .. seealso::
1463            :func:`~GSASIIscriptable.G2Phase.atom`
1464            :class:`G2AtomRecord`
1465        """
1466        ptrs = self.data['General']['AtomPtrs']
1467        output = []
1468        atoms = self.data['Atoms']
1469        for atom in atoms:
1470            output.append(G2AtomRecord(atom, ptrs, self.proj))
1471        return output
1472
1473    def histograms(self):
1474        output = []
1475        for hname in self.data.get('Histograms', {}).keys():
1476            output.append(self.proj.histogram(hname))
1477        return output
1478
1479    @property
1480    def ranId(self):
1481        return self.data['ranId']
1482
1483    @property
1484    def id(self):
1485        return G2obj.PhaseRanIdLookup[self.ranId]
1486
1487    def cell_dict(self):
1488        """Returns a dictionary of the cell parameters, with keys:
1489            'a', 'b', 'c', 'alpha', 'beta', 'gamma', 'vol'
1490
1491        :returns: a dict"""
1492        cell = self['General']['Cell']
1493        return {'a': cell[1], 'b': cell[2], 'c': cell[3],
1494                'alpha': cell[4], 'beta': cell[5], 'gamma': cell[6],
1495                'vol': cell[7]}
1496
1497    def set_refinements(self, refs):
1498        """Sets the refinement parameter 'key' to the specification 'value'
1499
1500        :param dict refs: A dictionary of the parameters to be set.
1501
1502        :returns: None"""
1503        for key, value in refs.items():
1504            if key == "Cell":
1505                self.data['General']['Cell'][0] = True
1506            elif key == "Atoms":
1507                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1508
1509                for atomlabel, atomrefinement in value.items():
1510                    if atomlabel == 'all':
1511                        for atom in self.atoms():
1512                            atom.refinement_flags = atomrefinement
1513                    else:
1514                        atom = self.atom(atomlabel)
1515                        atom.refinement_flags = atomrefinement
1516            elif key == "LeBail":
1517                hists = self.data['Histograms']
1518                for hname, hoptions in hists.items():
1519                    if 'LeBail' not in hoptions:
1520                        hoptions['newLeBail'] = True
1521                    hoptions['LeBail'] = bool(value)
1522            else:
1523                raise ValueError("Unknown key:", key)
1524
1525    def clear_refinements(self, refs):
1526        """Clears a given set of parameters.
1527
1528        :param dict refs: The parameters to clear"""
1529        for key, value in refs.items():
1530            if key == "Cell":
1531                self.data['General']['Cell'][0] = False
1532            elif key == "Atoms":
1533                cx, ct, cs, cia = self.data['General']['AtomPtrs']
1534
1535                for atomlabel in value:
1536                    atom = self.atom(atomlabel)
1537                    # Set refinement to none
1538                    atom.refinement_flags = ' '
1539            elif key == "LeBail":
1540                hists = self.data['Histograms']
1541                for hname, hoptions in hists.items():
1542                    if 'LeBail' not in hoptions:
1543                        hoptions['newLeBail'] = True
1544                    hoptions['LeBail'] = False
1545            else:
1546                raise ValueError("Unknown key:", key)
1547
1548    def set_HAP_refinements(self, refs, histograms='all'):
1549        """Sets the given HAP refinement parameters between this phase and
1550        the given histograms
1551
1552        :param dict refs: A dictionary of the parameters to be set.
1553        :param histograms: Either 'all' (default) or a list of the histograms
1554            whose HAP parameters will be set with this phase. Histogram and phase
1555            must already be associated
1556
1557        :returns: None
1558        """
1559        if histograms == 'all':
1560            histograms = self['Histograms'].values()
1561        else:
1562            histograms = [self['Histograms'][h.name] for h in histograms]
1563
1564        for key, val in refs.items():
1565            for h in histograms:
1566                if key == 'Babinet':
1567                    try:
1568                        sets = list(val)
1569                    except ValueError:
1570                        sets = ['BabA', 'BabU']
1571                    for param in sets:
1572                        if param not in ['BabA', 'BabU']:
1573                            raise ValueError("Not sure what to do with" + param)
1574                        for hist in histograms:
1575                            hist['Babinet'][param][1] = True
1576                elif key == 'Extinction':
1577                    for h in histograms:
1578                        h['Extinction'][1] = True
1579                elif key == 'HStrain':
1580                    for h in histograms:
1581                        hist['HStrain'][1] = [True for p in hist['Hstrain'][0]]
1582                elif key == 'Mustrain':
1583                    for h in histograms:
1584                        mustrain = h['Mustrain']
1585                        newType = None
1586                        if isinstance(val, str):
1587                            if val in ['isotropic', 'uniaxial', 'generalized']:
1588                                newType = val
1589                            else:
1590                                raise ValueError("Not a Mustrain type: " + val)
1591                        elif isinstance(val, dict):
1592                            newType = val.get('type', None)
1593                            direction = val.get('direction', None)
1594
1595                        if newType:
1596                            mustrain[0] = newType
1597                            if newType == 'isotropic':
1598                                mustrain[2][0] = True
1599                                mustrain[5] = [False for p in mustrain[4]]
1600                            elif newType == 'uniaxial':
1601                                pass
1602                            else:  # newtype == 'generalized'
1603                                mustrain[2] = [False for p in mustrain[1]]
1604                        if direction:
1605                            direction = [int(n) for n in direction]
1606                            if len(direction) != 3:
1607                                raise ValueError("Expected hkl, found", direction)
1608                            mustrain[3] = direction
1609                elif key == 'Pref.Ori.':
1610                    for h in histograms:
1611                        h['Pref.Ori.'][2] = True
1612                elif key == 'Show':
1613                    for h in histograms:
1614                        h['Show'] = True
1615                elif key == 'Size':
1616                    raise NotImplementedError()
1617                elif key == 'Use':
1618                    for h in histograms:
1619                        h['Use'] = True
1620                elif key == 'Scale':
1621                    for h in histograms:
1622                        h['Scale'][1] = False
1623
1624    def clear_HAP_refinements(self, refs, histograms='all'):
1625        """Clears the given HAP refinement parameters between this phase and
1626        the given histograms
1627
1628        :param dict refs: A dictionary of the parameters to be cleared.
1629        :param histograms: Either 'all' (default) or a list of the histograms
1630            whose HAP parameters will be cleared with this phase. Histogram and
1631            phase must already be associated
1632
1633        :returns: None
1634        """
1635        if histograms == 'all':
1636            histograms = self['Histograms'].values()
1637        else:
1638            histograms = [self['Histograms'][h.name] for h in histograms]
1639
1640        for key, val in refs.items():
1641            for h in histograms:
1642                if key == 'Babinet':
1643                    try:
1644                        sets = list(val)
1645                    except ValueError:
1646                        sets = ['BabA', 'BabU']
1647                    for param in sets:
1648                        if param not in ['BabA', 'BabU']:
1649                            raise ValueError("Not sure what to do with" + param)
1650                        for hist in histograms:
1651                            hist['Babinet'][param][1] = False
1652                elif key == 'Extinction':
1653                    for h in histograms:
1654                        h['Extinction'][1] = False
1655                elif key == 'HStrain':
1656                    for h in histograms:
1657                        hist['HStrain'][1] = [False for p in hist['Hstrain'][0]]
1658                elif key == 'Mustrain':
1659                    for h in histograms:
1660                        mustrain = h['Mustrain']
1661                        mustrain[2] = [False for p in mustrain[1]]
1662                        mustrain[5] = [False for p in mustrain[4]]
1663                elif key == 'Pref.Ori.':
1664                    for h in histograms:
1665                        h['Pref.Ori.'][2] = False
1666                elif key == 'Show':
1667                    for h in histograms:
1668                        h['Show'] = False
1669                elif key == 'Size':
1670                    raise NotImplementedError()
1671                elif key == 'Use':
1672                    for h in histograms:
1673                        h['Use'] = False
1674                elif key == 'Scale':
1675                    for h in histograms:
1676                        h['Scale'][1] = False
1677
1678
1679##########################
1680# Command Line Interface #
1681##########################
1682
1683
1684def create(*args):
1685    """The create subcommand.
1686
1687    Should be passed all the command-line arguments after `create`"""
1688    proj = G2Project(filename=args[0])
1689
1690    isPhase = False
1691    isPowderData = False
1692    isInstPrms = False
1693    instPrms = None
1694
1695    # TODO how to associate phase with histogram?
1696    for arg in args[1:]:
1697        if arg == '--phases':
1698            isPhase = True
1699            isPowderData = False
1700            isInstPrms = False
1701        elif arg == '--powder':
1702            isPhase = False
1703            isPowderData = True
1704            isInstPrms = False
1705        # Instrument parameters must be specified before
1706        # any powder data files are passed
1707        elif arg == '--iparams':
1708            isPhase = False
1709            isPowderData = False
1710            isInstPrms = True
1711        elif isPhase:
1712            proj.add_phase(arg)
1713        elif isPowderData:
1714            proj.add_powder_histogram(arg, instPrms)
1715        elif isInstPrms:
1716            instPrms = arg
1717            isInstPrms = False
1718        else:
1719            print("Not sure what to do with: {}".format(arg))
1720
1721    proj.save()
1722
1723
1724def dump(*args):
1725    """The dump subcommand"""
1726    raw = True
1727    histograms = False
1728    phases = False
1729    filenames = []
1730    for arg in args:
1731        if arg in ['-h', '--histograms']:
1732            histograms = True
1733            raw = False
1734        elif arg in ['-p', '--phases']:
1735            phases = True
1736            raw = False
1737        elif arg in ['-r', '--raw']:
1738            raw = True
1739        else:
1740            filenames.append(arg)
1741    if raw:
1742        import IPython.lib.pretty as pretty
1743    for fname in filenames:
1744        if raw:
1745            proj, nameList = LoadDictFromProjFile(fname)
1746            print("file:", fname)
1747            print("names:", nameList)
1748            for key, val in proj.items():
1749                print(key, ":")
1750                pretty.pprint(val)
1751        else:
1752            proj = G2Project(fname)
1753            if histograms:
1754                hists = proj.histograms()
1755                for h in hists:
1756                    print(fname, "hist", h.id, h.name)
1757            if phases:
1758                phase_list = proj.phases()
1759                for p in phase_list:
1760                    print(fname, "phase", p.id, p.name)
1761
1762
1763def IPyBrowse(*args):
1764    """Load a .gpx file and then open a IPython shell to browse it
1765    """
1766    filename = []
1767    for arg in args:
1768        fname = arg
1769        proj, nameList = LoadDictFromProjFile(fname)
1770        msg = "\nfile {} loaded into proj (dict) with names in nameList".format(fname)
1771        GSASIIpath.IPyBreak_base(msg)
1772        break
1773
1774
1775def refine(*args):
1776    """The refine subcommand"""
1777    proj = G2Project(args[0])
1778    proj.refine()
1779
1780
1781def seqrefine(*args):
1782    """The seqrefine subcommand"""
1783    raise NotImplementedError("seqrefine is not yet implemented")
1784
1785
1786def export(*args):
1787    """The export subcommand"""
1788    # Export CIF or Structure or ...
1789    raise NotImplementedError("export is not yet implemented")
1790
1791
1792subcommands = {"create": create,
1793               "dump": dump,
1794               "refine": refine,
1795               "seqrefine": seqrefine,
1796               "export": export,
1797               "browse": IPyBrowse}
1798
1799
1800def main():
1801    '''TODO: the command-line options need some thought
1802    '''
1803    argv = sys.argv
1804    if len(argv) > 1 and argv[1] in subcommands:
1805        subcommands[argv[1]](*argv[2:])
1806    elif len(argv) == 1 or argv[1] in ('help', '--help', '-h'):
1807        # TODO print usage
1808        subcommand_names = ' | '.join(sorted(subcommands.keys()))
1809        print("USAGE: {} [ {} ] ...".format(argv[0], subcommand_names))
1810    else:
1811        print("Unknown subcommand: {}".format(argv[1]))
1812        print("Available subcommands:")
1813        for name in sorted(subcommands.keys()):
1814            print("\t{}".format(name))
1815        sys.exit(-1)
1816    # arg = sys.argv
1817    # print(arg)
1818    # if len(arg) > 1:
1819    #     GPXfile = arg[1]
1820    #     if not ospath.exists(GPXfile):
1821    #         print(u'ERROR - '+GPXfile+u" doesn't exist!")
1822    #         exit()
1823    #     Project,nameList = LoadDictFromProjFile(GPXfile)
1824    #     SaveDictToProjFile(Project,nameList,'testout.gpx')
1825    # else:
1826    #     print('ERROR - missing filename')
1827    #     exit()
1828    # print("Done")
1829
1830if __name__ == '__main__':
1831    main()
1832
1833    # from gpx_manipulatons.py
1834    # try:
1835    #     filename, authorname = sys.argv[1:3]
1836    #     proj, names = make_empty_project(authorname, filename)
1837    #     SaveDictToProjFile(proj, names, os.path.abspath(filename))
1838    # except ValueError:
1839    #     print("Usage: {} <filename> <author>".format(sys.argv[0]))
1840    #     sys.exit(-1)
1841
1842
1843    # from refinements.py
1844#      USAGE = """USAGE: {} datafile instparams phasefile projectname refinements
1845
1846# datafile:    Input powder data
1847# intparams:   Corresponding instrument parameter file
1848# phasefile:   Phase to refine against data
1849# projectname: Project file to be created, should end in .gpx
1850# refinements: JSON file of refinements to be executed
1851# """
1852#     try:
1853#         datafile, instprm, phasefile, projectname, refinements = sys.argv[1:]
1854#     except ValueError:
1855#         print(USAGE.format(sys.argv[0]))
1856#         sys.exit(-1)
1857
1858#     try:
1859#         with open(refinements) as f:
1860#             refinements = json.load(f)
1861#     except IOError:
1862#         print("No such refinements file: {}".format(refinements))
1863
1864#     print("Creating project file \"{}\"...".format(projectname))
1865#     proj = G2Project(filename=projectname)
1866#     # Add the histogram
1867#     hist = proj.add_powder_histogram(datafile, instprm)
1868#     # Add the phase, and associate it with the histogram
1869#     proj.add_phase(phasefile, histograms=[hist.name])
1870
1871#     proj.do_refinements(refinements['refinements'])
1872#     proj.save()
1873
1874
1875    # from gpx_dumper
1876    # import IPython.lib.pretty as pretty
1877    # proj, nameList = LoadDictFromProjFile(sys.argv[1])
1878    # print("names:", nameList)
1879    # for key, val in proj.items():
1880    #     print(key, ":", sep='')
1881    #     pretty.pprint(val)
Note: See TracBrowser for help on using the repository browser.