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

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

More sphinx docstrings, start of command line interface, wrappers have project reference

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