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