source: branch/2frame/GSASIIscriptable.py @ 2974

Last change on this file since 2974 was 2974, checked in by odonnell, 6 years ago

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