source: trunk/GSASIIscriptable.py @ 3007

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

G2Project.refine does not back up, refactored phase/histogram access, HAP params

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
18
19Look at :class:`G2Project` to start.
20"""
21from __future__ import division, print_function # needed?
22import os.path as ospath
23import sys
24import cPickle
25import imp
26import copy
27import os
28import random as ran
29
30import numpy.ma as ma
31import scipy.interpolate as si
32import numpy as np
33import scipy as sp
34
35import GSASIIpath
36GSASIIpath.SetBinaryPath(False) # would rather have this in __name__ == '__main__' stanza
37import GSASIIobj as G2obj
38import GSASIIpwd as G2pwd
39import GSASIIstrMain as G2strMain
40import GSASIIspc as G2spc
41import GSASIIElem as G2elem
42
43# Delay imports to not slow down small scripts
44G2fil = None
45PwdrDataReaders = []
46PhaseReaders = []
47
48
49def LoadG2fil():
50    """Delay importing this module, it is slow"""
51    global G2fil
52    global PwdrDataReaders
53    global PhaseReaders
54    if G2fil is None:
55        import GSASIIfiles
56        G2fil = GSASIIfiles
57        PwdrDataReaders = G2fil.LoadImportRoutines("pwd", "Powder_Data")
58        PhaseReaders = G2fil.LoadImportRoutines("phase", "Phase")
59
60
61def LoadDictFromProjFile(ProjFile):
62    '''Read a GSAS-II project file and load items to dictionary
63    :param str ProjFile: GSAS-II project (name.gpx) full file name
64    :returns: Project,nameList, where
65
66      * Project (dict) is a representation of gpx file following the GSAS-II tree struture
67        for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict
68        data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)}
69      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
70
71    Example for fap.gpx::
72
73      Project = {                 #NB:dict order is not tree order
74        u'Phases':{'data':None,'fap':{phase dict}},
75        u'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc},
76        u'Rigid bodies':{'data': {rigid body dict}},
77        u'Covariance':{'data':{covariance data dict}},
78        u'Controls':{'data':{controls data dict}},
79        u'Notebook':{'data':[notebook list]},
80        u'Restraints':{'data':{restraint data dict}},
81        u'Constraints':{'data':{constraint data dict}}]
82        }
83      nameList = [                #NB: reproduces tree order
84        [u'Notebook',],
85        [u'Controls',],
86        [u'Covariance',],
87        [u'Constraints',],
88        [u'Restraints',],
89        [u'Rigid bodies',],
90        [u'PWDR FAP.XRA Bank 1',
91             u'Comments',
92             u'Limits',
93             u'Background',
94             u'Instrument Parameters',
95             u'Sample Parameters',
96             u'Peak List',
97             u'Index Peak List',
98             u'Unit Cells List',
99             u'Reflection Lists'],
100        [u'Phases', u'fap']
101        ]
102    '''
103    # Let IOError be thrown if file does not exist
104    # if not ospath.exists(ProjFile):
105    #     print ('\n*** Error attempt to open project file that does not exist:\n   '+
106    #         str(ProjFile))
107    #     return
108    file = open(ProjFile,'rb')
109    # print('loading from file: {}'.format(ProjFile))
110    Project = {}
111    nameList = []
112    try:
113        while True:
114            try:
115                data = cPickle.load(file)
116            except EOFError:
117                break
118            datum = data[0]
119            Project[datum[0]] = {'data':datum[1]}
120            nameList.append([datum[0],])
121            for datus in data[1:]:
122                Project[datum[0]][datus[0]] = datus[1]
123                nameList[-1].append(datus[0])
124        file.close()
125        # print('project load successful')
126    except:
127        raise IOError("Error reading file "+str(ProjFile)+". This is not a GSAS-II .gpx file")
128    finally:
129        file.close()
130    return Project,nameList
131
132def SaveDictToProjFile(Project,nameList,ProjFile):
133    '''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile
134
135    :param dict Project: representation of gpx file following the GSAS-II
136      tree structure as described for LoadDictFromProjFile
137    :param list nameList: names of main tree entries & subentries used to reconstruct project file
138    :param str ProjFile: full file name for output project.gpx file (including extension)
139    '''
140    file = open(ProjFile,'wb')
141    # print('save to file: {}'.format(ProjFile))
142    try:
143        for name in nameList:
144            data = []
145            item = Project[name[0]]
146            data.append([name[0],item['data']])
147            for item2 in name[1:]:
148                data.append([item2,item[item2]])
149            cPickle.dump(data,file,1)
150    finally:
151        file.close()
152    # print('project save successful')
153
154def ImportPowder(reader,filename):
155    '''Use a reader to import a powder diffraction data file
156
157    :param str reader: a scriptable reader
158    :param str filename: full name of powder data file; can be "multi-Bank" data
159
160    :returns list rdlist: list of reader objects containing powder data, one for each
161      "Bank" of data encountered in file. Items in reader object of interest are:
162
163        * rd.comments: list of str: comments found on powder file
164        * rd.dnames: list of str: data nammes suitable for use in GSASII data tree NB: duplicated in all rd entries in rdlist
165        * rd.powderdata: list of numpy arrays: pos,int,wt,zeros,zeros,zeros as needed for a PWDR entry in  GSASII data tree.
166    '''
167    rdfile,rdpath,descr = imp.find_module(reader)
168    rdclass = imp.load_module(reader,rdfile,rdpath,descr)
169    rd = rdclass.GSAS_ReaderClass()
170    if not rd.scriptable:
171        print(u'**** ERROR: '+reader+u' is not a scriptable reader')
172        return None
173    fl = open(filename,'rb')
174    rdlist = []
175    if rd.ContentsValidator(fl):
176        fl.seek(0)
177        repeat = True
178        rdbuffer = {} # create temporary storage for file reader
179        block = 0
180        while repeat: # loop if the reader asks for another pass on the file
181            block += 1
182            repeat = False
183            rd.objname = ospath.basename(filename)
184            flag = rd.Reader(filename,fl,None,buffer=rdbuffer,blocknum=block,)
185            if flag:
186                rdlist.append(copy.deepcopy(rd)) # save the result before it is written over
187                if rd.repeat:
188                    repeat = True
189        return rdlist
190    print(rd.errors)
191    return None
192
193def SetDefaultDData(dType,histoName,NShkl=0,NDij=0):
194    '''Create an initial Histogram dictionary
195
196    Author: Jackson O'Donnell jacksonhodonnell@gmail.com
197    '''
198    if dType in ['SXC','SNC']:
199        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
200            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
201            'Extinction':['Lorentzian','None', {'Tbar':0.1,'Cos2TM':0.955,
202            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}],
203            'Flack':[0.0,False]}
204    elif dType == 'SNT':
205        return {'Histogram':histoName,'Show':False,'Scale':[1.0,True],
206            'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]},
207            'Extinction':['Lorentzian','None', {
208            'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}]}
209    elif 'P' in dType:
210        return {'Histogram':histoName,'Show':False,'Scale':[1.0,False],
211            'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{},[],0.1],
212            'Size':['isotropic',[1.,1.,1.],[False,False,False],[0,0,1],
213                [1.,1.,1.,0.,0.,0.],6*[False,]],
214            'Mustrain':['isotropic',[1000.0,1000.0,1.0],[False,False,False],[0,0,1],
215                NShkl*[0.01,],NShkl*[False,]],
216            'HStrain':[NDij*[0.0,],NDij*[False,]],
217            'Extinction':[0.0,False],'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}}
218
219
220def PreSetup(data):
221    '''Create part of an initial (empty) phase dictionary
222
223    from GSASIIphsGUI.py, near end of UpdatePhaseData
224    Author: Jackson O'Donnell jacksonhodonnell@gmail.com
225    '''
226    if 'RBModels' not in data:
227        data['RBModels'] = {}
228    if 'MCSA' not in data:
229        data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[.8,1.2],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}}
230    if 'dict' in str(type(data['MCSA']['Results'])):
231        data['MCSA']['Results'] = []
232    if 'Modulated' not in data['General']:
233        data['General']['Modulated'] = False
234    if 'modulated' in data['General']['Type']:
235        data['General']['Modulated'] = True
236        data['General']['Type'] = 'nuclear'
237
238
239def SetupGeneral(data, dirname):
240    """Helps initialize phase data.
241
242    From GSASIIphsGui.py, function of the same name. Minor changes for imports etc.
243    Author: Jackson O'Donnell jacksonhodonnell@gmail.com
244    """
245    mapDefault = {'MapType':'','RefList':'','Resolution':0.5,'Show bonds':True,
246                'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False}
247    generalData = data['General']
248    atomData = data['Atoms']
249    generalData['AtomTypes'] = []
250    generalData['Isotopes'] = {}
251
252    if 'Isotope' not in generalData:
253        generalData['Isotope'] = {}
254    if 'Data plot type' not in generalData:
255        generalData['Data plot type'] = 'Mustrain'
256    if 'POhkl' not in generalData:
257        generalData['POhkl'] = [0,0,1]
258    if 'Map' not in generalData:
259        generalData['Map'] = mapDefault.copy()
260    if 'Flip' not in generalData:
261        generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None',
262            'k-factor':0.1,'k-Max':20.,}
263    if 'testHKL' not in generalData['Flip']:
264        generalData['Flip']['testHKL'] = [[0,0,2],[2,0,0],[1,1,1],[0,2,0],[1,2,3]]
265    if 'doPawley' not in generalData:
266        generalData['doPawley'] = False     #ToDo: change to ''
267    if 'Pawley dmin' not in generalData:
268        generalData['Pawley dmin'] = 1.0
269    if 'Pawley neg wt' not in generalData:
270        generalData['Pawley neg wt'] = 0.0
271    if 'Algolrithm' in generalData.get('MCSA controls',{}) or \
272        'MCSA controls' not in generalData:
273        generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50],
274        'dmin':2.0,'Algorithm':'log','Jump coeff':[0.95,0.5],'boltzmann':1.0,
275        'fast parms':[1.0,1.0,1.0],'log slope':0.9,'Cycles':1,'Results':[],'newDmin':True}
276    if 'AtomPtrs' not in generalData:
277        generalData['AtomPtrs'] = [3,1,7,9]
278        if generalData['Type'] == 'macromolecular':
279            generalData['AtomPtrs'] = [6,4,10,12]
280        elif generalData['Type'] == 'magnetic':
281            generalData['AtomPtrs'] = [3,1,10,12]
282    if generalData['Type'] in ['modulated',]:
283        generalData['Modulated'] = True
284        generalData['Type'] = 'nuclear'
285        if 'Super' not in generalData:
286            generalData['Super'] = 1
287            generalData['SuperVec'] = [[0,0,.1],False,4]
288            generalData['SSGData'] = {}
289        if '4DmapData' not in generalData:
290            generalData['4DmapData'] = mapDefault.copy()
291            generalData['4DmapData'].update({'MapType':'Fobs'})
292    if 'Modulated' not in generalData:
293        generalData['Modulated'] = False
294    if 'HydIds' not in generalData:
295        generalData['HydIds'] = {}
296    cx,ct,cs,cia = generalData['AtomPtrs']
297    generalData['NoAtoms'] = {}
298    generalData['BondRadii'] = []
299    generalData['AngleRadii'] = []
300    generalData['vdWRadii'] = []
301    generalData['AtomMass'] = []
302    generalData['Color'] = []
303    if generalData['Type'] == 'magnetic':
304        generalData['MagDmin'] = generalData.get('MagDmin',1.0)
305        landeg = generalData.get('Lande g',[])
306    generalData['Mydir'] = dirname
307    badList = {}
308    for iat,atom in enumerate(atomData):
309        atom[ct] = atom[ct].lower().capitalize()              #force to standard form
310        if generalData['AtomTypes'].count(atom[ct]):
311            generalData['NoAtoms'][atom[ct]] += atom[cx+3]*float(atom[cs+1])
312        elif atom[ct] != 'UNK':
313            Info = G2elem.GetAtomInfo(atom[ct])
314            if not Info:
315                if atom[ct] not in badList:
316                    badList[atom[ct]] = 0
317                badList[atom[ct]] += 1
318                atom[ct] = 'UNK'
319                continue
320            atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo
321            generalData['AtomTypes'].append(atom[ct])
322            generalData['Z'] = Info['Z']
323            generalData['Isotopes'][atom[ct]] = Info['Isotopes']
324            generalData['BondRadii'].append(Info['Drad'])
325            generalData['AngleRadii'].append(Info['Arad'])
326            generalData['vdWRadii'].append(Info['Vdrad'])
327            if atom[ct] in generalData['Isotope']:
328                if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]:
329                    isotope = generalData['Isotopes'][atom[ct]].keys()[-1]
330                    generalData['Isotope'][atom[ct]] = isotope
331                generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][atom[ct]]]['Mass'])
332            else:
333                generalData['Isotope'][atom[ct]] = 'Nat. Abund.'
334                if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]:
335                    isotope = generalData['Isotopes'][atom[ct]].keys()[-1]
336                    generalData['Isotope'][atom[ct]] = isotope
337                generalData['AtomMass'].append(Info['Mass'])
338            generalData['NoAtoms'][atom[ct]] = atom[cx+3]*float(atom[cs+1])
339            generalData['Color'].append(Info['Color'])
340            if generalData['Type'] == 'magnetic':
341                if len(landeg) < len(generalData['AtomTypes']):
342                    landeg.append(2.0)
343    if generalData['Type'] == 'magnetic':
344        generalData['Lande g'] = landeg[:len(generalData['AtomTypes'])]
345
346    if badList:
347        msg = 'Warning: element symbol(s) not found:'
348        for key in badList:
349            msg += '\n\t' + key
350            if badList[key] > 1:
351                msg += ' (' + str(badList[key]) + ' times)'
352        raise Exception("Phase error:\n" + msg)
353        # wx.MessageBox(msg,caption='Element symbol error')
354    F000X = 0.
355    F000N = 0.
356    for i,elem in enumerate(generalData['AtomTypes']):
357        F000X += generalData['NoAtoms'][elem]*generalData['Z']
358        isotope = generalData['Isotope'][elem]
359        F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope]['SL'][0]
360    generalData['F000X'] = F000X
361    generalData['F000N'] = F000N
362    import GSASIImath as G2mth
363    generalData['Mass'] = G2mth.getMass(generalData)
364
365def make_empty_project(author=None, filename=None):
366    """Creates an dictionary in the style of GSASIIscriptable, for an empty
367    project.
368
369    If no author name or filename is supplied, 'no name' and
370    <current dir>/test_output.gpx are used , respectively.
371
372    Returns: project dictionary, name list
373    Author: Jackson O'Donnell jacksonhodonnell@gmail.com
374    """
375    if not filename:
376        filename = os.path.join(os.getcwd(), 'test_output.gpx')
377    else:
378        filename = os.path.abspath(filename)
379    gsasii_version = str(GSASIIpath.GetVersionNumber())
380    LoadG2fil()
381    import matplotlib as mpl
382    python_library_versions = G2fil.get_python_versions([mpl, np, sp])
383
384    controls_data = dict(G2obj.DefaultControls)
385    controls_data['LastSavedAs'] = unicode(filename)
386    controls_data['LastSavedUsing'] = gsasii_version
387    controls_data['PythonVersions'] = python_library_versions
388    if author:
389        controls_data['Author'] = author
390
391    output = {'Constraints': {'data': {'HAP': [], 'Hist': [], 'Phase': [],
392                                       'Global': []}},
393              'Controls': {'data': controls_data},
394              u'Covariance': {'data': {}},
395              u'Notebook': {'data': ['']},
396              u'Restraints': {'data': {}},
397              u'Rigid bodies': {'data': {'RBIds': {'Residue': [], 'Vector': []},
398                                'Residue': {'AtInfo': {}},
399                                'Vector':  {'AtInfo': {}}}}}
400
401    names = [[u'Notebook'], [u'Controls'], [u'Covariance'],
402             [u'Constraints'], [u'Restraints'], [u'Rigid bodies']]
403
404    return output, names
405
406
407class G2ImportException(Exception):
408    pass
409
410
411def import_generic(filename, readerlist):
412    """Attempt to import a filename, using a list of reader objects.
413
414    Returns the first reader object which worked."""
415    # Translated from OnImportGeneric method in GSASII.py
416    primaryReaders, secondaryReaders = [], []
417    for reader in readerlist:
418        flag = reader.ExtensionValidator(filename)
419        if flag is None:
420            secondaryReaders.append(reader)
421        elif flag:
422            primaryReaders.append(reader)
423    if not secondaryReaders and not primaryReaders:
424        raise G2ImportException("Could not read file: ", filename)
425
426    with open(filename, 'Ur') as fp:
427        rd_list = []
428
429        for rd in primaryReaders + secondaryReaders:
430            # Initialize reader
431            rd.selections = []
432            rd.dnames = []
433            rd.ReInitialize()
434            # Rewind file
435            fp.seek(0)
436            rd.errors = ""
437            if not rd.ContentsValidator(fp):
438                # Report error
439                pass
440            if len(rd.selections) > 1:
441                # Select data?
442                # GSASII.py:543
443                raise G2ImportException("Not sure what data to select")
444
445            block = 0
446            rdbuffer = {}
447            fp.seek(0)
448            repeat = True
449            while repeat:
450                repeat = False
451                block += 1
452                rd.objname = os.path.basename(filename)
453                flag = rd.Reader(filename, fp, buffer=rdbuffer, blocknum=block)
454                if flag:
455                    # Omitting image loading special cases
456                    rd.readfilename = filename
457                    rd_list.append(copy.deepcopy(rd))
458                    repeat = rd.repeat
459                else:
460                    raise G2ImportException("{}.Reader() returned:".format(rd),
461                                            flag)
462
463            if rd_list:
464                if rd.warnings:
465                    print("Read warning by", rd.formatName, "reader:",
466                          rd.warnings, file=sys.stderr)
467                return rd_list
468    raise G2ImportException("No reader could read file: " + filename)
469
470
471def load_iprms(instfile, reader):
472    """Loads instrument parameters from a file, and edits the
473    given reader.
474
475    Returns a 2-tuple of (Iparm1, Iparm2) parameters
476    """
477    LoadG2fil()
478    ext = os.path.splitext(instfile)[1]
479
480    if ext.lower() == '.instprm':
481        # New GSAS File, load appropriately
482        with open(instfile) as f:
483            lines = f.readlines()
484        bank = reader.powderentry[2]
485        numbanks = reader.numbanks
486        iparms = G2fil.ReadPowderInstprm(lines, bank, numbanks, reader)
487
488        reader.instfile = instfile
489        reader.instmsg = 'GSAS-II file' + instfile
490        return iparms
491    elif ext.lower() not in ('.prm', '.inst', '.ins'):
492        raise ValueError('Expected .prm file, found: ', instfile)
493
494    # It's an old GSAS file, load appropriately
495    Iparm = {}
496    with open(instfile, 'Ur') as fp:
497        for line in fp:
498            if '#' in line:
499                continue
500            Iparm[line[:12]] = line[12:-1]
501    ibanks = int(Iparm.get('INS   BANK  ', '1').strip())
502    if ibanks == 1:
503        reader.instbank = 1
504        reader.powderentry[2] = 1
505        reader.instfile = instfile
506        reader.instmsg = instfile + ' bank ' + str(reader.instbank)
507        return G2fil.SetPowderInstParms(Iparm, reader)
508    # TODO handle >1 banks
509    raise NotImplementedError("Check GSASIIfiles.py: ReadPowderInstprm")
510
511def load_pwd_from_reader(reader, instprm, existingnames=[]):
512    """Loads powder data from a reader object, and assembles it into a GSASII data tree.
513
514    Returns: (name, tree) - 2-tuple of the histogram name (str), and data
515    Author: Jackson O'Donnell jacksonhodonnell@gmail.com
516    """
517    HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_')
518    HistName = unicode(G2obj.MakeUniqueLabel(HistName, existingnames))
519
520    try:
521        Iparm1, Iparm2 = instprm
522    except ValueError:
523        Iparm1, Iparm2 = load_iprms(instprm, reader)
524
525    Ymin = np.min(reader.powderdata[1])
526    Ymax = np.max(reader.powderdata[1])
527    valuesdict = {'wtFactor': 1.0,
528                  'Dummy': False,
529                  'ranId': ran.randint(0, sys.maxint),
530                  'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax,
531                  'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax,
532                  'qPlot': False, 'dPlot': False, 'sqrtPlot': False,
533                  'Yminmax': [Ymin, Ymax]}
534    reader.Sample['ranId'] = valuesdict['ranId']
535
536    # Ending keys:
537    # [u'Reflection Lists',
538    #  u'Limits',
539    #  'data',
540    #  u'Index Peak List',
541    #  u'Comments',
542    #  u'Unit Cells List',
543    #  u'Sample Parameters',
544    #  u'Peak List',
545    #  u'Background',
546    #  u'Instrument Parameters']
547    Tmin = np.min(reader.powderdata[0])
548    Tmax = np.max(reader.powderdata[0])
549
550    default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0],
551                          {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}]
552
553    output_dict = {u'Reflection Lists': {},
554                   u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]),
555                   u'data': [valuesdict, reader.powderdata, HistName],
556                   u'Index Peak List': [[], []],
557                   u'Comments': reader.comments,
558                   u'Unit Cells List': [],
559                   u'Sample Parameters': reader.Sample,
560                   u'Peak List': {'peaks': [], 'sigDict': {}},
561                   u'Background': reader.pwdparms.get('Background', default_background),
562                   u'Instrument Parameters': [Iparm1, Iparm2],
563                   }
564
565    names = [u'Comments',
566             u'Limits',
567             u'Background',
568             u'Instrument Parameters',
569             u'Sample Parameters',
570             u'Peak List',
571             u'Index Peak List',
572             u'Unit Cells List',
573             u'Reflection Lists']
574
575    # TODO controls?? GSASII.py:1664-7
576
577    return HistName, [HistName] + names, output_dict
578
579
580def _deep_copy_into(from_, into):
581    """Helper function for reloading .gpx file. See G2Project.reload()
582    Author: Jackson O'Donnell jacksonhodonnell@gmail.com
583    """
584    if isinstance(from_, dict) and isinstance(into, dict):
585        combined_keys = set(from_.keys()).union(into.keys())
586        for key in combined_keys:
587            if key in from_ and key in into:
588                both_dicts = (isinstance(from_[key], dict)
589                              and isinstance(into[key], dict))
590                both_lists = (isinstance(from_[key], list)
591                              and isinstance(into[key], list))
592                if both_dicts or both_lists:
593                    _deep_copy_into(from_[key], into[key])
594                else:
595                    into[key] = from_[key]
596            elif key in from_:
597                into[key] = from_[key]
598            else:  # key in into
599                del into[key]
600    elif isinstance(from_, list) and isinstance(into, list):
601        if len(from_) == len(into):
602            for i in xrange(len(from_)):
603                both_dicts = (isinstance(from_[i], dict)
604                              and isinstance(into[i], dict))
605                both_lists = (isinstance(from_[i], list)
606                              and isinstance(into[i], list))
607                if both_dicts or both_lists:
608                    _deep_copy_into(from_[i], into[i])
609                else:
610                    into[i] = from_[i]
611        else:
612            into[:] = from_
613
614
615class G2ObjectWrapper(object):
616    """Base class for all GSAS-II object wrappers.
617
618    The underlying GSAS-II format can be accessed as `wrapper.data`. A number
619    of overrides are implemented so that the wrapper behaves like a dictionary.
620
621    Author: Jackson O'Donnell jacksonhodonnell@gmail.com
622    """
623    def __init__(self, datadict):
624        self.data = datadict
625
626    def __getitem__(self, key):
627        return self.data[key]
628
629    def __setitem__(self, key, value):
630        self.data[key] = value
631
632    def __contains__(self, key):
633        return key in self.data
634
635    def get(self, k, d=None):
636        return self.data.get(k, d)
637
638    def keys(self):
639        return self.data.keys()
640
641    def values(self):
642        return self.data.values()
643
644    def items(self):
645        return self.data.items()
646
647
648class G2Project(G2ObjectWrapper):
649    """Represents an entire GSAS-II project."""
650    def __init__(self, gpxfile=None, author=None, filename=None):
651        """Loads a GSAS-II project from a specified filename.
652
653        :param str gpxfile: Existing .gpx file to be loaded. If nonexistent,
654        creates an empty project.
655        :param str author: Author's name. Optional.
656        :param str filename: The filename the project should be saved to in
657        the future. If both filename and gpxfile are present, the project is
658        loaded from the gpxfile, then set to  save to filename in the future"""
659        if gpxfile is None:
660            filename = os.path.abspath(os.path.expanduser(filename))
661            self.filename = filename
662            self.data, self.names = make_empty_project(author=author, filename=filename)
663        elif isinstance(gpxfile, str):
664            # TODO set author, filename
665            self.filename = os.path.abspath(os.path.expanduser(gpxfile))
666            self.data, self.names = LoadDictFromProjFile(gpxfile)
667            self.index_ids()
668        else:
669            raise ValueError("Not sure what to do with gpxfile")
670
671    @classmethod
672    def from_dict_and_names(cls, gpxdict, names, filename=None):
673        """Creates a :class:`G2Project` directly from
674        a dictionary and a list of names. If in doubt, do not use this.
675
676        :returns: a :class:`G2Project`
677        """
678        out = cls()
679        if filename:
680            filename = os.path.abspath(os.path.expanduser(filename))
681            out.filename = filename
682            gpxdict['Controls']['data']['LastSavedAs'] = filename
683        else:
684            try:
685                out.filename = gpxdict['Controls']['data']['LastSavedAs']
686            except KeyError:
687                out.filename = None
688        out.data = gpxdict
689        out.names = names
690
691    def save(self, filename=None):
692        """Saves the project, either to the current filename, or to a new file.
693
694        Updates self.filename if a new filename provided"""
695        # TODO update LastSavedUsing ?
696        if filename:
697            filename = os.path.abspath(os.path.expanduser(filename))
698            self.data['Controls']['data']['LastSavedAs'] = filename
699            self.filename = filename
700        elif not self.filename:
701            raise AttributeError("No file name to save to")
702        SaveDictToProjFile(self.data, self.names, self.filename)
703
704    def add_powder_histogram(self, datafile, iparams):
705        """Loads a powder data histogram into the project.
706
707        Automatically checks for an instrument parameter file, or one can be
708        provided.
709
710        :param str datafile: The powder data file to read, a filename.
711        :param str iparams: The instrument parameters file, a filename.
712
713        :returns: A :class:`G2PwdrData` object representing
714        the histogram
715        """
716        LoadG2fil()
717        datafile = os.path.abspath(os.path.expanduser(datafile))
718        iparams = os.path.abspath(os.path.expanduser(iparams))
719        pwdrreaders = import_generic(datafile, PwdrDataReaders)
720        histname, new_names, pwdrdata = load_pwd_from_reader(
721                                          pwdrreaders[0], iparams,
722                                          [h.name for h in self.histograms()])
723        if histname in self.data:
724            print("Warning - redefining histogram", histname)
725        else:
726            if self.names[-1][0] == 'Phases':
727                self.names.insert(-1, new_names)
728            else:
729                self.names.append(new_names)
730        self.data[histname] = pwdrdata
731        return self.histogram(histname)
732
733    def add_phase(self, phasefile, phasename=None, histograms=[]):
734        """Loads a phase into the project from a .cif file
735
736        :param str phasefile: The CIF file from which to import the phase.
737        :param str phasename: The name of the new phase, or None for the default
738        :param list histograms: The names of the histograms to associate with
739        this phase
740
741        :returns: A :class:`G2Phase` object representing the
742        new phase.
743        """
744        LoadG2fil()
745        phasefile = os.path.abspath(os.path.expanduser(phasefile))
746
747        # TODO handle multiple phases in a file
748        phasereaders = import_generic(phasefile, PhaseReaders)
749        phasereader = phasereaders[0]
750
751        phasename = phasename or phasereader.Phase['General']['Name']
752        phaseNameList = [p.name for p in self.phases()]
753        phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList)
754        phasereader.Phase['General']['Name'] = phasename
755
756        if 'Phases' not in self.data:
757            self.data[u'Phases'] = { 'data': None }
758        assert phasename not in self.data['Phases'], "phase names should be unique"
759        self.data['Phases'][phasename] = phasereader.Phase
760
761        if phasereader.Constraints:
762            Constraints = self.data['Constraints']
763            for i in phasereader.Constraints:
764                if isinstance(i, dict):
765                    if '_Explain' not in Constraints:
766                        Constraints['_Explain'] = {}
767                    Constraints['_Explain'].update(i)
768                else:
769                    Constraints['Phase'].append(i)
770
771        data = self.data['Phases'][phasename]
772        generalData = data['General']
773        SGData = generalData['SGData']
774        NShkl = len(G2spc.MustrainNames(SGData))
775        NDij = len(G2spc.HStrainNames(SGData))
776        Super = generalData.get('Super', 0)
777        if Super:
778            SuperVec = np.array(generalData['SuperVec'][0])
779        else:
780            SuperVec = []
781        UseList = data['Histograms']
782
783        for histname in histograms:
784            if histname.startswith('HKLF '):
785                raise NotImplementedError("Does not support HKLF yet")
786            elif histname.startswith('PWDR '):
787                hist = self.histogram(histname)
788                hist['Reflection Lists'][generalData['Name']] = {}
789                UseList[histname] = SetDefaultDData('PWDR', histname, NShkl=NShkl, NDij=NDij)
790                for key, val in [('Use', True), ('LeBail', False),
791                                 ('newLeBail', True),
792                                 ('Babinet', {'BabA': [0.0, False],
793                                              'BabU': [0.0, False]})]:
794                    if key not in UseList[histname]:
795                        UseList[histname][key] = val
796            else:
797                raise NotImplementedError("Unexpected histogram" + histname)
798
799        for obj in self.names:
800            if obj[0] == 'Phases':
801                phasenames = obj
802                break
803        else:
804            phasenames = [u'Phases']
805            self.names.append(phasenames)
806        phasenames.append(unicode(phasename))
807
808        # TODO should it be self.filename, not phasefile?
809        SetupGeneral(data, os.path.dirname(phasefile))
810        self.index_ids()
811
812        return self.phase(phasename)
813
814    def reload(self):
815        """Reload self from self.filename"""
816        data, names = LoadDictFromProjFile(self.filename)
817        self.names = names
818        # Need to deep copy the new data file data into the current tree,
819        # so that any existing G2Phase, or G2PwdrData objects will still be
820        # valid
821        _deep_copy_into(from_=data, into=self.data)
822
823    def refine(self, newfile=None, printFile=None, makeBack=False):
824        # index_ids will automatically save the project
825        self.index_ids()
826        # TODO G2strMain does not properly use printFile
827        G2strMain.Refine(self.filename, makeBack=makeBack)
828        # Reload yourself
829        self.reload()
830
831    def histogram_names(self):
832        """Gives a list of the names of each histogram in the project.
833
834        :returns: list of strings
835
836        .. seealso::
837            :func:`~GSASIIscriptable.G2Project.histogram`
838            :func:`~GSASIIscriptable.G2Project.histograms`
839            :func:`~GSASIIscriptable.G2Project.phase`
840            :func:`~GSASIIscriptable.G2Project.phases`
841            :func:`~GSASIIscriptable.G2Project.phase_names`
842            """
843        output = []
844        for obj in self.names:
845            if len(obj) > 1 and obj[0] != u'Phases':
846                output.append(obj[0])
847        return output
848
849    def histogram(self, histname):
850        """Returns the histogram named histname, or None if it does not exist.
851
852        :param histname: The name of the histogram (str), or ranId or index.
853        :returns: A :class:`G2PwdrData` object, or None if
854        the histogram does not exist
855
856        .. seealso::
857            :func:`~GSASIIscriptable.G2Project.histograms`
858            :func:`~GSASIIscriptable.G2Project.phase`
859            :func:`~GSASIIscriptable.G2Project.phases`
860            """
861        if histname in self.data:
862            return G2PwdrData(self.data[histname], self)
863        histRanId = None
864        for key, val in G2obj.HistIdLookup.items():
865            name, ranId = val
866            # histname can be either ranId (key) or index (val)
867            if ranId == histname or key == str(histname):
868                return self.histogram(name)
869
870    def histograms(self):
871        """Return a list of all histograms, as
872        :class:`G2PwdrData` objects
873
874        .. seealso::
875            :func:`~GSASIIscriptable.G2Project.histograms`
876            :func:`~GSASIIscriptable.G2Project.phase`
877            :func:`~GSASIIscriptable.G2Project.phases`
878            """
879        output = []
880        for obj in self.names:
881            if len(obj) > 1 and obj[0] != u'Phases':
882                output.append(self.histogram(obj[0]))
883        return output
884
885    def phase_names(self):
886        """Gives an list of the names of each phase in the project.
887
888        :returns: A list of strings
889
890        .. seealso::
891            :func:`~GSASIIscriptable.G2Project.histogram`
892            :func:`~GSASIIscriptable.G2Project.histograms`
893            :func:`~GSASIIscriptable.G2Project.phase`
894            :func:`~GSASIIscriptable.G2Project.phases`
895            """
896        for obj in self.names:
897            if obj[0] == 'Phases':
898                return obj[1:]
899        return []
900
901    def phase(self, phasename):
902        """
903        Gives an object representing the specified phase in this project.
904
905        :param str phasename: The name of the desired phase
906        :returns: A :class:`G2Phase` object
907        :raises: KeyError
908
909        .. seealso::
910            :func:`~GSASIIscriptable.G2Project.histograms`
911            :func:`~GSASIIscriptable.G2Project.phase`
912            :func:`~GSASIIscriptable.G2Project.phases`
913            """
914        phases = self.data['Phases']
915        if phasename in phases:
916            return G2Phase(phases[phasename], phasename, self)
917        phaseRanId = None
918        for key, val in G2obj.PhaseIdLookup.items():
919            name, ranId = val
920            # phasename can be either ranId (key) or index (val)
921            if ranId == phasename or key == str(phasename):
922                return self.phase(name)
923
924    def phases(self):
925        """
926        Returns a list of all the phases in the project.
927
928        :returns: A :class:`G2Phase`
929
930        .. seealso::
931            :func:`~GSASIIscriptable.G2Project.histogram`
932            :func:`~GSASIIscriptable.G2Project.histograms`
933            :func:`~GSASIIscriptable.G2Project.phase`
934            """
935        for obj in self.names:
936            if obj[0] == 'Phases':
937                return [self.phase(p) for p in obj[1:]]
938        return []
939
940    def do_refinements(self, refinements, histogram='all', phase='all',
941                       outputnames=None):
942        """Conducts a series of refinements.
943
944        :param list refinements: A list of dictionaries defining refinements
945        :param str histogram: Name of histogram for refinements to be applied
946        to, or 'all'
947        :param str phase: Name of phase for refinements to be applied to, or
948        'all'
949        """
950        if outputnames:
951            if len(refinements) != len(outputnames):
952                raise ValueError("Should have same number of outuputs to"
953                                 "refinements")
954        else:
955            outputnames = [None for r in refinements]
956
957        for output, refinement in zip(outputnames, refinements):
958            self.set_refinement(refinement, histogram)
959            # Handle 'once' args - refinements that are disabled after this
960            # refinement
961            if 'once' in refinement:
962                temp = {'set': refinement['once']}
963                self.set_refinement(temp, histogram, phase)
964
965            if output:
966                self.save(output)
967
968            self.refine()  # newFile=output)
969
970            # Handle 'once' args - refinements that are disabled after this
971            # refinement
972            if 'once' in refinement:
973                temp = {'clear': refinement['once']}
974                self.set_refinement(temp, histogram, phase)
975
976    def set_refinement(self, refinement, histogram='all', phase='all'):
977        """Apply specified refinements to a given histogram(s) or phase(s).
978
979        Refinement parameters are categorize in three groups:
980
981        1. Histogram parameters
982        2. Phase parameters
983        3. Histogram-and-Phase (HAP) parameters
984
985        :param dict refinement: The refinements to be conducted
986        :param histogram: Either a name of a histogram (str), a list of
987        histogram names, or 'all' (default)
988        :param phase: Either a name of a phase (str), a list of phase names, or
989        'all' (default)
990
991        .. seealso::
992            :func:`~G2PwdrData.set_refinements`
993            :func:`~G2PwdrData.clear_refinements`
994            :func:`~G2Phase.set_refinements`
995            :func:`~G2Phase.clear_refinements`
996            :func:`~G2Phase.set_HAP_refinements`
997            :func:`~G2Phase.clear_HAP_refinements`"""
998
999        if histogram == 'all':
1000            hists = self.histograms()
1001        elif isinstance(histogram, str) or isinstance(histogram, int):
1002            hists = [self.histogram(histogram)]
1003        else:
1004            hists = [self.histogram(name) for name in histogram]
1005
1006        if phase == 'all':
1007            phases = self.phases()
1008        elif isinstance(phase, str) or isinstance(phase, int):
1009            phases = [self.phase(phase)]
1010        else:
1011            phases = [self.phase(name) for name in phase]
1012
1013
1014        # TODO: HAP parameters:
1015        #   Babinet
1016        #   Extinction
1017        #   HStrain
1018        #   Mustrain
1019        #   Pref. Ori
1020        #   Size
1021
1022        pwdr_set = {}
1023        phase_set = {}
1024        hap_set = {}
1025        for key, val in refinement.get('set', {}).items():
1026            # Apply refinement options
1027            if G2PwdrData.is_valid_refinement_key(key):
1028                pwdr_set[key] = val
1029            elif G2Phase.is_valid_refinement_key(key):
1030                phase_set[key] = val
1031            elif G2Phase.is_valid_HAP_refinement_key(key):
1032                hap_set[key] = val
1033            else:
1034                raise ValueError("Unknown refinement key", key)
1035
1036        for hist in hists:
1037            hist.set_refinements(pwdr_set)
1038        for phase in phases:
1039            phase.set_refinements(phase_set)
1040        for phase in phases:
1041            phase.set_HAP_refinements(hap_set, hists)
1042
1043        pwdr_clear = {}
1044        phase_clear = {}
1045        hap_clear = {}
1046        for key, val in refinement.get('clear', {}).items():
1047            # Clear refinement options
1048            if G2PwdrData.is_valid_refinement_key(key):
1049                pwdr_clear[key] = val
1050            elif G2Phase.is_valid_refinement_key(key):
1051                phase_clear[key] = val
1052            elif G2Phase.is_valid_HAP_refinement_key(key):
1053                hap_set[key] = val
1054            else:
1055                raise ValueError("Unknown refinement key", key)
1056
1057        for hist in hists:
1058            hist.clear_refinements(pwdr_clear)
1059        for phase in phases:
1060            phase.clear_refinements(phase_clear)
1061        for phase in phases:
1062            phase.clear_HAP_refinements(hap_clear, hists)
1063
1064    def index_ids(self):
1065        import GSASIIstrIO as G2strIO
1066        self.save()
1067        return G2strIO.GetUsedHistogramsAndPhases(self.filename)
1068
1069    def add_constraint_raw(self, cons_scope, constr):
1070        """Adds a constraint of type consType to the project.
1071        cons_scope should be one of "Hist", "Phase", "HAP", or "Global".
1072
1073        WARNING it does not check the constraint is well-constructed"""
1074        constrs = self['Constraints']['data']
1075        if 'Global' not in constrs:
1076            constrs['Global'] = []
1077        constrs[cons_scope].append(constr)
1078
1079    def hold_many(self, vars, type):
1080        """Apply holds for all the variables in vars, for constraint of a given type.
1081
1082        type is passed directly to add_constraint_raw as consType
1083
1084        :param list vars: A list of variables to hold. Either G2obj.G2VarObj objects,
1085        string variable specifiers, or arguments for :meth:`make_var_obj`
1086        :param str type: A string constraint type specifier. See
1087        :class:`~GSASIIscriptable.G2Project.add_constraint_raw`"""
1088        for var in vars:
1089            if isinstance(var, str):
1090                var = self.make_var_obj(var)
1091            elif not isinstance(var, G2obj.G2VarObj):
1092                var = self.make_var_obj(*var)
1093            self.add_constraint_raw(type, [[1.0, var], None, None, 'h'])
1094
1095    def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None,
1096                     reloadIdx=True):
1097        """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a")
1098        or individual names of phase, histogram, varname, and atomId.
1099
1100        Automatically converts string phase, hist, or atom names into the ID required
1101        by G2VarObj."""
1102
1103        if reloadIdx:
1104            self.index_ids()
1105
1106        # If string representation, short circuit
1107        if hist is None and varname is None and atomId is None:
1108            if isinstance(phase, str) and ':' in phase:
1109                return G2obj.G2VarObj(phase)
1110
1111        # Get phase index
1112        phaseObj = None
1113        if isinstance(phase, G2Phase):
1114            phaseObj = phase
1115            phase = G2obj.PhaseRanIdLookup[phase.ranId]
1116        elif phase in self['Phases']:
1117            phaseObj = self.phase(phase)
1118            phaseRanId = phaseObj.ranId
1119            phase = G2obj.PhaseRanIdLookup[phaseRanId]
1120
1121        # Get histogram index
1122        if isinstance(hist, G2PwdrData):
1123            hist = G2obj.HistRanIdLookup[hist.ranId]
1124        elif hist in self.data:
1125            histRanId = self.histogram(hist).ranId
1126            hist = G2obj.HistRanIdLookup[histRanId]
1127
1128        # Get atom index (if any)
1129        if isinstance(atomId, G2AtomRecord):
1130            atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId]
1131        elif phaseObj:
1132            atomObj = phaseObj.atom(atomId)
1133            if atomObj:
1134                atomRanId = atomObj.ranId
1135                atomId = G2obj.AtomRanIdLookup[phase][atomRanId]
1136
1137        return G2obj.G2VarObj(phase, hist, varname, atomId)
1138
1139
1140class G2AtomRecord(G2ObjectWrapper):
1141    """Wrapper for an atom record. Has convenient accessors via @property.
1142
1143    Example:
1144
1145    >>> atom = some_phase.atom("O3")
1146    >>> # We can access the underlying data format
1147    >>> atom.data
1148    ['O3', 'O-2', '', ... ]
1149    >>> # We can also use wrapper accessors
1150    >>> atom.coordinates
1151    (0.33, 0.15, 0.5)
1152    >>> atom.refinement_flags
1153    u'FX'
1154    >>> atom.ranId
1155    4615973324315876477
1156    """
1157    def __init__(self, data, indices, proj):
1158        self.data = data
1159        self.cx, self.ct, self.cs, self.cia = indices
1160        self.proj = proj
1161
1162    @property
1163    def label(self):
1164        return self.data[self.ct-1]
1165
1166    @property
1167    def type(self):
1168        return self.data[self.ct]
1169
1170    @property
1171    def refinement_flags(self):
1172        return self.data[self.ct+1]
1173
1174    @refinement_flags.setter
1175    def refinement_flags(self, other):
1176        # Automatically check it is a valid refinement
1177        for c in other:
1178            if c not in ' FXU':
1179                raise ValueError("Invalid atom refinement: ", other)
1180        self.data[self.ct+1] = unicode(other)
1181
1182    @property
1183    def coordinates(self):
1184        return tuple(self.data[self.cx:self.cx+3])
1185
1186    @property
1187    def ranId(self):
1188        return self.data[self.cia+8]
1189
1190    @property
1191    def adp_flag(self):
1192        # Either 'I' or 'A'
1193        return self.data[self.cia]
1194
1195    @property
1196    def uiso(self):
1197        if self.adp_flag == 'I':
1198            return self.data[self.cia+1]
1199        else:
1200            return self.data[self.cia+2:self.cia+8]
1201
1202
1203class G2PwdrData(G2ObjectWrapper):
1204    """Wraps a Powder Data Histogram."""
1205    def __init__(self, data, proj):
1206        self.data = data
1207        self.proj = proj
1208
1209    @staticmethod
1210    def is_valid_refinement_key(key):
1211        valid_keys = ['Limits', 'Sample Parameters', 'Background',
1212                      'Instrument Parameters']
1213        return key in valid_keys
1214
1215    @property
1216    def name(self):
1217        return self['data'][-1]
1218
1219    @property
1220    def ranId(self):
1221        return self['data'][0]['ranId']
1222
1223    @property
1224    def residuals(self):
1225        data = self['data'][0]
1226        return {key: data[key]
1227                for key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']}
1228
1229    @property
1230    def id(self):
1231        return G2obj.HistRanIdLookup[self.ranId]
1232
1233    def fit_fixed_points(self):
1234        """Attempts to apply a background fit to the fixed points currently specified."""
1235        def SetInstParms(Inst):
1236            dataType = Inst['Type'][0]
1237            insVary = []
1238            insNames = []
1239            insVals = []
1240            for parm in Inst:
1241                insNames.append(parm)
1242                insVals.append(Inst[parm][1])
1243                if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha',
1244                    'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]:
1245                        Inst[parm][2] = False
1246            instDict = dict(zip(insNames, insVals))
1247            instDict['X'] = max(instDict['X'], 0.01)
1248            instDict['Y'] = max(instDict['Y'], 0.01)
1249            if 'SH/L' in instDict:
1250                instDict['SH/L'] = max(instDict['SH/L'], 0.002)
1251            return dataType, instDict, insVary
1252
1253        bgrnd = self['Background']
1254
1255        # Need our fixed points in order
1256        bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0])
1257        X = [x for x, y in bgrnd[1]['FixedPoints']]
1258        Y = [y for x, y in bgrnd[1]['FixedPoints']]
1259
1260        limits = self['Limits'][1]
1261        if X[0] > limits[0]:
1262            X = [limits[0]] + X
1263            Y = [Y[0]] + Y
1264        if X[-1] < limits[1]:
1265            X += [limits[1]]
1266            Y += [Y[-1]]
1267
1268        # Some simple lookups
1269        controls = self.proj['Controls']['data']
1270        inst, inst2 = self['Instrument Parameters']
1271        pwddata = self['data'][1]
1272
1273        # Construct the data for background fitting
1274        xBeg = np.searchsorted(pwddata[0], limits[0])
1275        xFin = np.searchsorted(pwddata[0], limits[1])
1276        xdata = pwddata[0][xBeg:xFin]
1277        ydata = si.interp1d(X,Y)(ma.getdata(xdata))
1278
1279        W = [1]*len(xdata)
1280        Z = [0]*len(xdata)
1281
1282        dataType, insDict, insVary = SetInstParms(inst)
1283        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1284
1285        # Do the fit
1286        data = np.array([xdata, ydata, W, Z, Z, Z])
1287        G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data,
1288                        prevVaryList=bakVary, controls=controls)
1289
1290        # Post-fit
1291        parmDict = {}
1292        bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd)
1293        parmDict.update(bakDict)
1294        parmDict.update(insDict)
1295        pwddata[3][xBeg:xFin] *= 0
1296        pwddata[5][xBeg:xFin] *= 0
1297        pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0]
1298
1299        # TODO adjust pwddata? GSASIIpwdGUI.py:1041
1300        # TODO update background
1301        self.proj.save()
1302
1303    def y_calc(self):
1304        return self['data'][1][3]
1305
1306    def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True):
1307        try:
1308            import matplotlib.pyplot as plt
1309            data = self['data'][1]
1310            if Yobs:
1311                plt.plot(data[0], data[1], label='Yobs')
1312            if Ycalc:
1313                plt.plot(data[0], data[3], label='Ycalc')
1314            if Background:
1315                plt.plot(data[0], data[4], label='Background')
1316            if Residual:
1317                plt.plot(data[0], data[5], label="Residual")
1318        except ImportError:
1319            pass
1320
1321    def set_refinements(self, refs):
1322        """Sets the refinement parameter 'key' to the specification 'value'
1323
1324        :param dict refs: A dictionary of the parameters to be set.
1325
1326        :returns: None"""
1327        do_fit_fixed_points = False
1328        for key, value in refs.items():
1329            if key == 'Limits':
1330                old_limits = self.data['Limits'][1]
1331                new_limits = value
1332                if isinstance(new_limits, dict):
1333                    if 'low' in new_limits:
1334                        old_limits[0] = new_limits['low']
1335                    if 'high' in new_limits:
1336                        old_limits[1] = new_limits['high']
1337                else:
1338                    old_limits[0], old_limits[1] = new_limits
1339            elif key == 'Sample Parameters':
1340                sample = self.data['Sample Parameters']
1341                for sparam in value:
1342                    sample[sparam][1] = True
1343            elif key == 'Background':
1344                bkg, peaks = self.data['Background']
1345
1346                # If True or False, just set the refine parameter
1347                if value in (True, False):
1348                    bkg[1] = value
1349                    return
1350
1351                if 'type' in value:
1352                    bkg[0] = value['type']
1353                if 'refine' in value:
1354                    bkg[1] = value['refine']
1355                if 'no. coeffs' in value:
1356                    cur_coeffs = bkg[2]
1357                    n_coeffs = value['no. coeffs']
1358                    if n_coeffs > cur_coeffs:
1359                        for x in range(n_coeffs - cur_coeffs):
1360                            bkg.append(0.0)
1361                    else:
1362                        for _ in range(cur_coeffs - n_coeffs):
1363                            bkg.pop()
1364                    bkg[2] = n_coeffs
1365                if 'coeffs' in value:
1366                    bkg[3:] = value['coeffs']
1367                if 'FixedPoints' in value:
1368                    peaks['FixedPoints'] = [(float(a), float(b))
1369                                            for a, b in value['FixedPoints']]
1370                if value.get('fit fixed points', False):
1371                    do_fit_fixed_points = True
1372
1373            elif key == 'Instrument Parameters':
1374                instrument, secondary = self.data['Instrument Parameters']
1375                for iparam in value:
1376                    try:
1377                        instrument[iparam][2] = True
1378                    except IndexError:
1379                        raise ValueError("Invalid key:", iparam)
1380            else:
1381                raise ValueError("Unknown key:", key)
1382        # Fit fixed points after the fact - ensure they are after fixed points
1383        # are added
1384        if do_fit_fixed_points:
1385            # Background won't be fit if refinement flag not set
1386            orig = self['Background'][0][1]
1387            self['Background'][0][1] = True
1388            self.fit_fixed_points()
1389            # Restore the previous value
1390            self['Background'][0][1] = orig
1391
1392    def clear_refinements(self, refs):
1393        """Clears the refinement parameter 'key' and its associated value.
1394
1395        :param dict refs: A dictionary of parameters to clear."""
1396        for key, value in refs.items():
1397            if key == 'Limits':
1398                old_limits, cur_limits = self.data['Limits']
1399                cur_limits[0], cur_limits[1] = old_limits
1400            elif key == 'Sample Parameters':
1401                sample = self.data['Sample Parameters']
1402                for sparam in value:
1403                    sample[sparam][1] = False
1404            elif key == 'Background':
1405                bkg, peaks = self.data['Background']
1406
1407                # If True or False, just set the refine parameter
1408                if value in (True, False):
1409                    bkg[1] = False
1410                    return
1411
1412                bkg[1] = False
1413                if 'FixedPoints' in value:
1414                    if 'FixedPoints' in peaks:
1415                        del peaks['FixedPoints']
1416            elif key == 'Instrument Parameters':
1417                instrument, secondary = self.data['Instrument Parameters']
1418                for iparam in value:
1419                    instrument[iparam][2] = False
1420            else:
1421                raise ValueError("Unknown key:", key)
1422
1423
1424class G2Phase(G2ObjectWrapper):
1425    """A wrapper object around a given phase.
1426    Author: Jackson O'Donnell jacksonhodonnell@gmail.com
1427    """
1428    def __init__(self, data, name, proj):
1429        self.data = data
1430        self.name = name
1431        self.proj = proj
1432
1433    @staticmethod
1434    def is_valid_refinement_key(key):
1435        valid_keys = ["Cell", "Atoms", "LeBail"]
1436        return key in valid_keys
1437
1438    @staticmethod
1439    def is_valid_HAP_refinement_key(key):
1440        valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain",
1441                      "Pref.Ori.", "Show", "Size", "Use", "Scale"]
1442        return key in valid_keys
1443
1444    def atom(self, atomlabel):
1445        """Returns the atom specified by atomlabel, or None if it does not
1446        exist.
1447
1448        :param str atomlabel: The name of the atom (e.g. "O2")
1449        :returns: A :class:`G2AtomRecord` object
1450        representing the atom.
1451        """
1452        # Consult GSASIIobj.py for the meaning of this
1453        cx, ct, cs, cia = self.data['General']['AtomPtrs']
1454        ptrs = [cx, ct, cs, cia]
1455        atoms = self.data['Atoms']
1456        for atom in atoms:
1457            if atom[ct-1] == atomlabel:
1458                return G2AtomRecord(atom, ptrs, self.proj)
1459
1460    def atoms(self):
1461        """Returns a list of atoms present in the phase.
1462
1463        :returns: A list of :class:`G2AtomRecord` objects. See
1464        also :func:`~GSASIIscriptable.G2Phase.atom`
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.