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. Run `python GSASIIscriptable.py --help` to |
---|
17 | show the available subcommands, and inspect each subcommand with |
---|
18 | `python GSASIIscriptable.py <subcommand> --help`. |
---|
19 | |
---|
20 | This file specifies several wrapper classes around GSAS-II data representations. |
---|
21 | They all inherit from :class:`G2ObjectWrapper`. The chief class is :class:`G2Project`, |
---|
22 | which represents an entire GSAS-II project and provides several methods to access |
---|
23 | phases, powder histograms, and execute Rietveld refinements. |
---|
24 | |
---|
25 | .. _Refinement_parameters_kinds: |
---|
26 | |
---|
27 | ===================== |
---|
28 | Refinement parameters |
---|
29 | ===================== |
---|
30 | |
---|
31 | There are three classes of refinement parameters: |
---|
32 | |
---|
33 | * Histogram. Turned on and off through :func:`G2PwdrData.set_refinements` |
---|
34 | and :func:`G2PwdrData.clear_refinements` |
---|
35 | * Phase. Turned on and off through :func:`G2Phase.set_refinements` |
---|
36 | and :func:`G2Phase.clear_refinements` |
---|
37 | * Histogram-and-phase (HAP). Turned on and off through |
---|
38 | :func:`G2Phase.set_HAP_refinements` and :func:`G2Phase.clear_HAP_refinements` |
---|
39 | |
---|
40 | |
---|
41 | ============================ |
---|
42 | Refinement specifiers format |
---|
43 | ============================ |
---|
44 | |
---|
45 | Refinement parameters are specified as dictionaries, supplied to any of the functions |
---|
46 | named in :ref:`Refinement_parameters_kinds`. Each method accepts a different set |
---|
47 | of keys, described below for each of the three parameter classes. |
---|
48 | |
---|
49 | .. _Histogram_parameters_table: |
---|
50 | |
---|
51 | -------------------- |
---|
52 | Histogram parameters |
---|
53 | -------------------- |
---|
54 | |
---|
55 | This table describes the dictionaries supplied to :func:`G2PwdrData.set_refinements` |
---|
56 | and :func:`G2PwdrData.clear_refinements`. |
---|
57 | |
---|
58 | Example: |
---|
59 | |
---|
60 | .. code-block:: python |
---|
61 | |
---|
62 | params = {'set': { 'Limits': [0.8, 12.0], |
---|
63 | 'Sample Parameters': ['Absorption', 'Contrast', 'DisplaceX'], |
---|
64 | 'Background': {'type': 'chebyschev', 'refine': True}}, |
---|
65 | 'clear': {'Instrument Parameters': ['U', 'V', 'W']}} |
---|
66 | some_histogram.set_refinements(params['set']) |
---|
67 | some_histogram.clear_refinements(params['clear']) |
---|
68 | |
---|
69 | .. tabularcolumns:: |l|p|p| |
---|
70 | |
---|
71 | ===================== ==================== ==================== |
---|
72 | key subkey explanation |
---|
73 | ===================== ==================== ==================== |
---|
74 | Limits The 2-theta range of values to consider. Can |
---|
75 | be either a dictionary of 'low' and/or 'high', |
---|
76 | or a list of 2 items [low, high] |
---|
77 | \ low Sets the low limit |
---|
78 | \ high Sets the high limit |
---|
79 | Sample Parameters Should be provided as a **list** of subkeys |
---|
80 | to set or clear, e.g. ['DisplaceX', 'Scale'] |
---|
81 | \ Absorption |
---|
82 | \ Contrast |
---|
83 | \ DisplaceX Sample displacement along the X direction |
---|
84 | \ DisplaceY Sample displacement along the Y direction |
---|
85 | \ Scale Histogram Scale factor |
---|
86 | Background Sample background. If value is a boolean, |
---|
87 | the background's 'refine' parameter is set |
---|
88 | to the given boolean. Usually should be a |
---|
89 | dictionary with any of the following keys: |
---|
90 | \ type The background model, e.g. 'chebyschev' |
---|
91 | \ refine The value of the refine flag, boolean |
---|
92 | \ no. coeffs Number of coefficients to use, integer |
---|
93 | \ coeffs List of floats, literal values for background |
---|
94 | \ FixedPoints List of (2-theta, intensity) values for fixed points |
---|
95 | \ fit fixed points If True, triggers a fit to the fixed points to be calculated. It is calculated when this key is detected, regardless of calls to refine. |
---|
96 | Instrument Parameters As in Sample Paramters, Should be provided as a **list** of subkeys to |
---|
97 | set or clear, e.g. ['X', 'Y', 'Zero', 'SH/L'] |
---|
98 | \ U, V, W All separate keys. Gaussian peak profile terms |
---|
99 | \ X, Y Separate keys. Lorentzian peak profile terms |
---|
100 | \ Zero Zero shift |
---|
101 | \ SH/L |
---|
102 | \ Polariz. Polarization parameter |
---|
103 | \ Lam Lambda, the incident wavelength |
---|
104 | ===================== ==================== ==================== |
---|
105 | |
---|
106 | |
---|
107 | .. _Phase_parameters_table: |
---|
108 | |
---|
109 | ---------------- |
---|
110 | Phase parameters |
---|
111 | ---------------- |
---|
112 | |
---|
113 | This table describes the dictionaries supplied to :func:`G2Phase.set_refinements` |
---|
114 | and :func:`G2Phase.clear_refinements`. |
---|
115 | |
---|
116 | Example: |
---|
117 | |
---|
118 | .. code-block:: python |
---|
119 | |
---|
120 | params = { 'LeBail': True, 'Cell': True, |
---|
121 | 'Atoms': { 'Mn1': 'X', |
---|
122 | 'O3': 'XU', |
---|
123 | 'V4': 'FXU'}} |
---|
124 | some_histogram.set_refinements(params) |
---|
125 | |
---|
126 | .. tabularcolumns:: |l|p| |
---|
127 | |
---|
128 | ===================== ==================== |
---|
129 | key explanation |
---|
130 | ===================== ==================== |
---|
131 | Cell Whether or not to refine the unit cell. |
---|
132 | Atoms Dictionary of atoms and refinement flags. |
---|
133 | Each key should be an atom label, e.g. |
---|
134 | 'O3', 'Mn5', and each value should be |
---|
135 | a string defining what values to refine. |
---|
136 | Values can be any combination of 'F' |
---|
137 | for fractional occupancy, 'X' for position, |
---|
138 | and 'U' for Debye-Waller factor |
---|
139 | LeBail Enables LeBail intensity extraction. |
---|
140 | ===================== ==================== |
---|
141 | |
---|
142 | .. _HAP_parameters_table: |
---|
143 | |
---|
144 | ------------------------------ |
---|
145 | Histogram-and-phase parameters |
---|
146 | ------------------------------ |
---|
147 | |
---|
148 | This table describes the dictionaries supplied to :func:`G2Phase.set_HAP_refinements` |
---|
149 | and :func:`G2Phase.clear_HAP_refinements`. |
---|
150 | |
---|
151 | Example: |
---|
152 | |
---|
153 | .. code-block:: python |
---|
154 | |
---|
155 | params = { 'Babinet': 'BabA', |
---|
156 | 'Extinction': True, |
---|
157 | 'Mustrain': { 'type': 'uniaxial', |
---|
158 | 'direction': [0, 0, 1], |
---|
159 | 'refine': True}} |
---|
160 | some_phase.set_HAP_refinements(params) |
---|
161 | |
---|
162 | .. tabularcolumns:: |l|p|p| |
---|
163 | |
---|
164 | ===================== ===================== ==================== |
---|
165 | key subkey explanation |
---|
166 | ===================== ===================== ==================== |
---|
167 | Babinet Should be a **list** of the following |
---|
168 | subkeys. If not, assumes both |
---|
169 | BabA and BabU |
---|
170 | \ BabA |
---|
171 | \ BabU |
---|
172 | Extinction Should be boolean, whether or not to |
---|
173 | refine. |
---|
174 | HStrain Should be boolean, whether or not to |
---|
175 | refine. |
---|
176 | Mustrain |
---|
177 | \ type Mustrain model. One of 'isotropic', |
---|
178 | 'uniaxial', or 'generalized' |
---|
179 | \ direction For uniaxial. A list of three integers, |
---|
180 | the [hkl] direction of the axis. |
---|
181 | \ refine Usually boolean, whether or not to refine. |
---|
182 | When in doubt, set it to true. |
---|
183 | For uniaxial model, can specify list |
---|
184 | of 'axial' or 'equatorial'. If boolean |
---|
185 | given sets both axial and equatorial. |
---|
186 | Pref.Ori. Boolean, True to refine |
---|
187 | Show Boolean, True to refine |
---|
188 | Size Not implemented |
---|
189 | Use Boolean, True to refine |
---|
190 | Scale Phase fraction; Boolean, True to refine |
---|
191 | ===================== ===================== ==================== |
---|
192 | |
---|
193 | ============================ |
---|
194 | Scriptable API |
---|
195 | ============================ |
---|
196 | """ |
---|
197 | from __future__ import division, print_function # needed? |
---|
198 | import argparse |
---|
199 | import os.path as ospath |
---|
200 | import datetime as dt |
---|
201 | import sys |
---|
202 | import cPickle |
---|
203 | import imp |
---|
204 | import copy |
---|
205 | import os |
---|
206 | import random as ran |
---|
207 | |
---|
208 | import numpy.ma as ma |
---|
209 | import scipy.interpolate as si |
---|
210 | import numpy as np |
---|
211 | import scipy as sp |
---|
212 | |
---|
213 | import GSASIIpath |
---|
214 | GSASIIpath.SetBinaryPath(False) # would rather have this in __name__ == '__main__' stanza |
---|
215 | import GSASIIobj as G2obj |
---|
216 | import GSASIIpwd as G2pwd |
---|
217 | import GSASIIstrMain as G2strMain |
---|
218 | import GSASIIspc as G2spc |
---|
219 | import GSASIIElem as G2elem |
---|
220 | |
---|
221 | # Delay imports to not slow down small scripts |
---|
222 | G2fil = None |
---|
223 | PwdrDataReaders = [] |
---|
224 | PhaseReaders = [] |
---|
225 | |
---|
226 | |
---|
227 | def LoadG2fil(): |
---|
228 | """Delay importing this module, it is slow""" |
---|
229 | global G2fil |
---|
230 | global PwdrDataReaders |
---|
231 | global PhaseReaders |
---|
232 | if G2fil is None: |
---|
233 | import GSASIIfiles |
---|
234 | G2fil = GSASIIfiles |
---|
235 | PwdrDataReaders = G2fil.LoadImportRoutines("pwd", "Powder_Data") |
---|
236 | PhaseReaders = G2fil.LoadImportRoutines("phase", "Phase") |
---|
237 | |
---|
238 | |
---|
239 | def LoadDictFromProjFile(ProjFile): |
---|
240 | '''Read a GSAS-II project file and load items to dictionary |
---|
241 | :param str ProjFile: GSAS-II project (name.gpx) full file name |
---|
242 | :returns: Project,nameList, where |
---|
243 | |
---|
244 | * Project (dict) is a representation of gpx file following the GSAS-II tree struture |
---|
245 | for each item: key = tree name (e.g. 'Controls','Restraints',etc.), data is dict |
---|
246 | data dict = {'data':item data whch may be list, dict or None,'subitems':subdata (if any)} |
---|
247 | * nameList (list) has names of main tree entries & subentries used to reconstruct project file |
---|
248 | |
---|
249 | Example for fap.gpx:: |
---|
250 | |
---|
251 | Project = { #NB:dict order is not tree order |
---|
252 | u'Phases':{'data':None,'fap':{phase dict}}, |
---|
253 | u'PWDR FAP.XRA Bank 1':{'data':[histogram data list],'Comments':comments,'Limits':limits, etc}, |
---|
254 | u'Rigid bodies':{'data': {rigid body dict}}, |
---|
255 | u'Covariance':{'data':{covariance data dict}}, |
---|
256 | u'Controls':{'data':{controls data dict}}, |
---|
257 | u'Notebook':{'data':[notebook list]}, |
---|
258 | u'Restraints':{'data':{restraint data dict}}, |
---|
259 | u'Constraints':{'data':{constraint data dict}}] |
---|
260 | } |
---|
261 | nameList = [ #NB: reproduces tree order |
---|
262 | [u'Notebook',], |
---|
263 | [u'Controls',], |
---|
264 | [u'Covariance',], |
---|
265 | [u'Constraints',], |
---|
266 | [u'Restraints',], |
---|
267 | [u'Rigid bodies',], |
---|
268 | [u'PWDR FAP.XRA Bank 1', |
---|
269 | u'Comments', |
---|
270 | u'Limits', |
---|
271 | u'Background', |
---|
272 | u'Instrument Parameters', |
---|
273 | u'Sample Parameters', |
---|
274 | u'Peak List', |
---|
275 | u'Index Peak List', |
---|
276 | u'Unit Cells List', |
---|
277 | u'Reflection Lists'], |
---|
278 | [u'Phases', u'fap'] |
---|
279 | ] |
---|
280 | ''' |
---|
281 | # Let IOError be thrown if file does not exist |
---|
282 | # if not ospath.exists(ProjFile): |
---|
283 | # print ('\n*** Error attempt to open project file that does not exist:\n '+ |
---|
284 | # str(ProjFile)) |
---|
285 | # return |
---|
286 | file = open(ProjFile,'rb') |
---|
287 | # print('loading from file: {}'.format(ProjFile)) |
---|
288 | Project = {} |
---|
289 | nameList = [] |
---|
290 | try: |
---|
291 | while True: |
---|
292 | try: |
---|
293 | data = cPickle.load(file) |
---|
294 | except EOFError: |
---|
295 | break |
---|
296 | datum = data[0] |
---|
297 | Project[datum[0]] = {'data':datum[1]} |
---|
298 | nameList.append([datum[0],]) |
---|
299 | for datus in data[1:]: |
---|
300 | Project[datum[0]][datus[0]] = datus[1] |
---|
301 | nameList[-1].append(datus[0]) |
---|
302 | file.close() |
---|
303 | # print('project load successful') |
---|
304 | except: |
---|
305 | raise IOError("Error reading file "+str(ProjFile)+". This is not a GSAS-II .gpx file") |
---|
306 | finally: |
---|
307 | file.close() |
---|
308 | return Project,nameList |
---|
309 | |
---|
310 | def SaveDictToProjFile(Project,nameList,ProjFile): |
---|
311 | '''Save a GSAS-II project file from dictionary/nameList created by LoadDictFromProjFile |
---|
312 | |
---|
313 | :param dict Project: representation of gpx file following the GSAS-II |
---|
314 | tree structure as described for LoadDictFromProjFile |
---|
315 | :param list nameList: names of main tree entries & subentries used to reconstruct project file |
---|
316 | :param str ProjFile: full file name for output project.gpx file (including extension) |
---|
317 | ''' |
---|
318 | file = open(ProjFile,'wb') |
---|
319 | try: |
---|
320 | for name in nameList: |
---|
321 | data = [] |
---|
322 | item = Project[name[0]] |
---|
323 | data.append([name[0],item['data']]) |
---|
324 | for item2 in name[1:]: |
---|
325 | data.append([item2,item[item2]]) |
---|
326 | cPickle.dump(data,file,1) |
---|
327 | finally: |
---|
328 | file.close() |
---|
329 | |
---|
330 | def ImportPowder(reader,filename): |
---|
331 | '''Use a reader to import a powder diffraction data file |
---|
332 | |
---|
333 | :param str reader: a scriptable reader |
---|
334 | :param str filename: full name of powder data file; can be "multi-Bank" data |
---|
335 | |
---|
336 | :returns list rdlist: list of reader objects containing powder data, one for each |
---|
337 | "Bank" of data encountered in file. Items in reader object of interest are: |
---|
338 | |
---|
339 | * rd.comments: list of str: comments found on powder file |
---|
340 | * rd.dnames: list of str: data nammes suitable for use in GSASII data tree NB: duplicated in all rd entries in rdlist |
---|
341 | * rd.powderdata: list of numpy arrays: pos,int,wt,zeros,zeros,zeros as needed for a PWDR entry in GSASII data tree. |
---|
342 | ''' |
---|
343 | rdfile,rdpath,descr = imp.find_module(reader) |
---|
344 | rdclass = imp.load_module(reader,rdfile,rdpath,descr) |
---|
345 | rd = rdclass.GSAS_ReaderClass() |
---|
346 | if not rd.scriptable: |
---|
347 | print(u'**** ERROR: '+reader+u' is not a scriptable reader') |
---|
348 | return None |
---|
349 | fl = open(filename,'rb') |
---|
350 | rdlist = [] |
---|
351 | if rd.ContentsValidator(fl): |
---|
352 | fl.seek(0) |
---|
353 | repeat = True |
---|
354 | rdbuffer = {} # create temporary storage for file reader |
---|
355 | block = 0 |
---|
356 | while repeat: # loop if the reader asks for another pass on the file |
---|
357 | block += 1 |
---|
358 | repeat = False |
---|
359 | rd.objname = ospath.basename(filename) |
---|
360 | flag = rd.Reader(filename,fl,None,buffer=rdbuffer,blocknum=block,) |
---|
361 | if flag: |
---|
362 | rdlist.append(copy.deepcopy(rd)) # save the result before it is written over |
---|
363 | if rd.repeat: |
---|
364 | repeat = True |
---|
365 | return rdlist |
---|
366 | print(rd.errors) |
---|
367 | return None |
---|
368 | |
---|
369 | def SetDefaultDData(dType,histoName,NShkl=0,NDij=0): |
---|
370 | '''Create an initial Histogram dictionary |
---|
371 | |
---|
372 | Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com) |
---|
373 | ''' |
---|
374 | if dType in ['SXC','SNC']: |
---|
375 | return {'Histogram':histoName,'Show':False,'Scale':[1.0,True], |
---|
376 | 'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}, |
---|
377 | 'Extinction':['Lorentzian','None', {'Tbar':0.1,'Cos2TM':0.955, |
---|
378 | 'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}], |
---|
379 | 'Flack':[0.0,False]} |
---|
380 | elif dType == 'SNT': |
---|
381 | return {'Histogram':histoName,'Show':False,'Scale':[1.0,True], |
---|
382 | 'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}, |
---|
383 | 'Extinction':['Lorentzian','None', { |
---|
384 | 'Eg':[1.e-10,False],'Es':[1.e-10,False],'Ep':[1.e-10,False]}]} |
---|
385 | elif 'P' in dType: |
---|
386 | return {'Histogram':histoName,'Show':False,'Scale':[1.0,False], |
---|
387 | 'Pref.Ori.':['MD',1.0,False,[0,0,1],0,{},[],0.1], |
---|
388 | 'Size':['isotropic',[1.,1.,1.],[False,False,False],[0,0,1], |
---|
389 | [1.,1.,1.,0.,0.,0.],6*[False,]], |
---|
390 | 'Mustrain':['isotropic',[1000.0,1000.0,1.0],[False,False,False],[0,0,1], |
---|
391 | NShkl*[0.01,],NShkl*[False,]], |
---|
392 | 'HStrain':[NDij*[0.0,],NDij*[False,]], |
---|
393 | 'Extinction':[0.0,False],'Babinet':{'BabA':[0.0,False],'BabU':[0.0,False]}} |
---|
394 | |
---|
395 | |
---|
396 | def PreSetup(data): |
---|
397 | '''Create part of an initial (empty) phase dictionary |
---|
398 | |
---|
399 | from GSASIIphsGUI.py, near end of UpdatePhaseData |
---|
400 | |
---|
401 | Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com) |
---|
402 | ''' |
---|
403 | if 'RBModels' not in data: |
---|
404 | data['RBModels'] = {} |
---|
405 | if 'MCSA' not in data: |
---|
406 | data['MCSA'] = {'Models':[{'Type':'MD','Coef':[1.0,False,[.8,1.2],],'axis':[0,0,1]}],'Results':[],'AtInfo':{}} |
---|
407 | if 'dict' in str(type(data['MCSA']['Results'])): |
---|
408 | data['MCSA']['Results'] = [] |
---|
409 | if 'Modulated' not in data['General']: |
---|
410 | data['General']['Modulated'] = False |
---|
411 | if 'modulated' in data['General']['Type']: |
---|
412 | data['General']['Modulated'] = True |
---|
413 | data['General']['Type'] = 'nuclear' |
---|
414 | |
---|
415 | |
---|
416 | def SetupGeneral(data, dirname): |
---|
417 | """Helps initialize phase data. |
---|
418 | |
---|
419 | From GSASIIphsGui.py, function of the same name. Minor changes for imports etc. |
---|
420 | |
---|
421 | Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com) |
---|
422 | """ |
---|
423 | mapDefault = {'MapType':'','RefList':'','Resolution':0.5,'Show bonds':True, |
---|
424 | 'rho':[],'rhoMax':0.,'mapSize':10.0,'cutOff':50.,'Flip':False} |
---|
425 | generalData = data['General'] |
---|
426 | atomData = data['Atoms'] |
---|
427 | generalData['AtomTypes'] = [] |
---|
428 | generalData['Isotopes'] = {} |
---|
429 | |
---|
430 | if 'Isotope' not in generalData: |
---|
431 | generalData['Isotope'] = {} |
---|
432 | if 'Data plot type' not in generalData: |
---|
433 | generalData['Data plot type'] = 'Mustrain' |
---|
434 | if 'POhkl' not in generalData: |
---|
435 | generalData['POhkl'] = [0,0,1] |
---|
436 | if 'Map' not in generalData: |
---|
437 | generalData['Map'] = mapDefault.copy() |
---|
438 | if 'Flip' not in generalData: |
---|
439 | generalData['Flip'] = {'RefList':'','Resolution':0.5,'Norm element':'None', |
---|
440 | 'k-factor':0.1,'k-Max':20.,} |
---|
441 | if 'testHKL' not in generalData['Flip']: |
---|
442 | generalData['Flip']['testHKL'] = [[0,0,2],[2,0,0],[1,1,1],[0,2,0],[1,2,3]] |
---|
443 | if 'doPawley' not in generalData: |
---|
444 | generalData['doPawley'] = False #ToDo: change to '' |
---|
445 | if 'Pawley dmin' not in generalData: |
---|
446 | generalData['Pawley dmin'] = 1.0 |
---|
447 | if 'Pawley neg wt' not in generalData: |
---|
448 | generalData['Pawley neg wt'] = 0.0 |
---|
449 | if 'Algolrithm' in generalData.get('MCSA controls',{}) or \ |
---|
450 | 'MCSA controls' not in generalData: |
---|
451 | generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50], |
---|
452 | 'dmin':2.0,'Algorithm':'log','Jump coeff':[0.95,0.5],'boltzmann':1.0, |
---|
453 | 'fast parms':[1.0,1.0,1.0],'log slope':0.9,'Cycles':1,'Results':[],'newDmin':True} |
---|
454 | if 'AtomPtrs' not in generalData: |
---|
455 | generalData['AtomPtrs'] = [3,1,7,9] |
---|
456 | if generalData['Type'] == 'macromolecular': |
---|
457 | generalData['AtomPtrs'] = [6,4,10,12] |
---|
458 | elif generalData['Type'] == 'magnetic': |
---|
459 | generalData['AtomPtrs'] = [3,1,10,12] |
---|
460 | if generalData['Type'] in ['modulated',]: |
---|
461 | generalData['Modulated'] = True |
---|
462 | generalData['Type'] = 'nuclear' |
---|
463 | if 'Super' not in generalData: |
---|
464 | generalData['Super'] = 1 |
---|
465 | generalData['SuperVec'] = [[0,0,.1],False,4] |
---|
466 | generalData['SSGData'] = {} |
---|
467 | if '4DmapData' not in generalData: |
---|
468 | generalData['4DmapData'] = mapDefault.copy() |
---|
469 | generalData['4DmapData'].update({'MapType':'Fobs'}) |
---|
470 | if 'Modulated' not in generalData: |
---|
471 | generalData['Modulated'] = False |
---|
472 | if 'HydIds' not in generalData: |
---|
473 | generalData['HydIds'] = {} |
---|
474 | cx,ct,cs,cia = generalData['AtomPtrs'] |
---|
475 | generalData['NoAtoms'] = {} |
---|
476 | generalData['BondRadii'] = [] |
---|
477 | generalData['AngleRadii'] = [] |
---|
478 | generalData['vdWRadii'] = [] |
---|
479 | generalData['AtomMass'] = [] |
---|
480 | generalData['Color'] = [] |
---|
481 | if generalData['Type'] == 'magnetic': |
---|
482 | generalData['MagDmin'] = generalData.get('MagDmin',1.0) |
---|
483 | landeg = generalData.get('Lande g',[]) |
---|
484 | generalData['Mydir'] = dirname |
---|
485 | badList = {} |
---|
486 | for iat,atom in enumerate(atomData): |
---|
487 | atom[ct] = atom[ct].lower().capitalize() #force to standard form |
---|
488 | if generalData['AtomTypes'].count(atom[ct]): |
---|
489 | generalData['NoAtoms'][atom[ct]] += atom[cx+3]*float(atom[cs+1]) |
---|
490 | elif atom[ct] != 'UNK': |
---|
491 | Info = G2elem.GetAtomInfo(atom[ct]) |
---|
492 | if not Info: |
---|
493 | if atom[ct] not in badList: |
---|
494 | badList[atom[ct]] = 0 |
---|
495 | badList[atom[ct]] += 1 |
---|
496 | atom[ct] = 'UNK' |
---|
497 | continue |
---|
498 | atom[ct] = Info['Symbol'] # N.B. symbol might be changed by GetAtomInfo |
---|
499 | generalData['AtomTypes'].append(atom[ct]) |
---|
500 | generalData['Z'] = Info['Z'] |
---|
501 | generalData['Isotopes'][atom[ct]] = Info['Isotopes'] |
---|
502 | generalData['BondRadii'].append(Info['Drad']) |
---|
503 | generalData['AngleRadii'].append(Info['Arad']) |
---|
504 | generalData['vdWRadii'].append(Info['Vdrad']) |
---|
505 | if atom[ct] in generalData['Isotope']: |
---|
506 | if generalData['Isotope'][atom[ct]] not in generalData['Isotopes'][atom[ct]]: |
---|
507 | isotope = generalData['Isotopes'][atom[ct]].keys()[-1] |
---|
508 | generalData['Isotope'][atom[ct]] = isotope |
---|
509 | generalData['AtomMass'].append(Info['Isotopes'][generalData['Isotope'][atom[ct]]]['Mass']) |
---|
510 | else: |
---|
511 | generalData['Isotope'][atom[ct]] = 'Nat. Abund.' |
---|
512 | if 'Nat. Abund.' not in generalData['Isotopes'][atom[ct]]: |
---|
513 | isotope = generalData['Isotopes'][atom[ct]].keys()[-1] |
---|
514 | generalData['Isotope'][atom[ct]] = isotope |
---|
515 | generalData['AtomMass'].append(Info['Mass']) |
---|
516 | generalData['NoAtoms'][atom[ct]] = atom[cx+3]*float(atom[cs+1]) |
---|
517 | generalData['Color'].append(Info['Color']) |
---|
518 | if generalData['Type'] == 'magnetic': |
---|
519 | if len(landeg) < len(generalData['AtomTypes']): |
---|
520 | landeg.append(2.0) |
---|
521 | if generalData['Type'] == 'magnetic': |
---|
522 | generalData['Lande g'] = landeg[:len(generalData['AtomTypes'])] |
---|
523 | |
---|
524 | if badList: |
---|
525 | msg = 'Warning: element symbol(s) not found:' |
---|
526 | for key in badList: |
---|
527 | msg += '\n\t' + key |
---|
528 | if badList[key] > 1: |
---|
529 | msg += ' (' + str(badList[key]) + ' times)' |
---|
530 | raise Exception("Phase error:\n" + msg) |
---|
531 | # wx.MessageBox(msg,caption='Element symbol error') |
---|
532 | F000X = 0. |
---|
533 | F000N = 0. |
---|
534 | for i,elem in enumerate(generalData['AtomTypes']): |
---|
535 | F000X += generalData['NoAtoms'][elem]*generalData['Z'] |
---|
536 | isotope = generalData['Isotope'][elem] |
---|
537 | F000N += generalData['NoAtoms'][elem]*generalData['Isotopes'][elem][isotope]['SL'][0] |
---|
538 | generalData['F000X'] = F000X |
---|
539 | generalData['F000N'] = F000N |
---|
540 | import GSASIImath as G2mth |
---|
541 | generalData['Mass'] = G2mth.getMass(generalData) |
---|
542 | |
---|
543 | |
---|
544 | def make_empty_project(author=None, filename=None): |
---|
545 | """Creates an dictionary in the style of GSASIIscriptable, for an empty |
---|
546 | project. |
---|
547 | |
---|
548 | If no author name or filename is supplied, 'no name' and |
---|
549 | <current dir>/test_output.gpx are used , respectively. |
---|
550 | |
---|
551 | Returns: project dictionary, name list |
---|
552 | |
---|
553 | Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com) |
---|
554 | """ |
---|
555 | if not filename: |
---|
556 | filename = os.path.join(os.getcwd(), 'test_output.gpx') |
---|
557 | else: |
---|
558 | filename = os.path.abspath(filename) |
---|
559 | gsasii_version = str(GSASIIpath.GetVersionNumber()) |
---|
560 | LoadG2fil() |
---|
561 | import matplotlib as mpl |
---|
562 | python_library_versions = G2fil.get_python_versions([mpl, np, sp]) |
---|
563 | |
---|
564 | controls_data = dict(G2obj.DefaultControls) |
---|
565 | controls_data['LastSavedAs'] = unicode(filename) |
---|
566 | controls_data['LastSavedUsing'] = gsasii_version |
---|
567 | controls_data['PythonVersions'] = python_library_versions |
---|
568 | if author: |
---|
569 | controls_data['Author'] = author |
---|
570 | |
---|
571 | output = {'Constraints': {'data': {'HAP': [], 'Hist': [], 'Phase': [], |
---|
572 | 'Global': []}}, |
---|
573 | 'Controls': {'data': controls_data}, |
---|
574 | u'Covariance': {'data': {}}, |
---|
575 | u'Notebook': {'data': ['']}, |
---|
576 | u'Restraints': {'data': {}}, |
---|
577 | u'Rigid bodies': {'data': {'RBIds': {'Residue': [], 'Vector': []}, |
---|
578 | 'Residue': {'AtInfo': {}}, |
---|
579 | 'Vector': {'AtInfo': {}}}}} |
---|
580 | |
---|
581 | names = [[u'Notebook'], [u'Controls'], [u'Covariance'], |
---|
582 | [u'Constraints'], [u'Restraints'], [u'Rigid bodies']] |
---|
583 | |
---|
584 | return output, names |
---|
585 | |
---|
586 | |
---|
587 | class G2ImportException(Exception): |
---|
588 | pass |
---|
589 | |
---|
590 | |
---|
591 | def import_generic(filename, readerlist): |
---|
592 | """Attempt to import a filename, using a list of reader objects. |
---|
593 | |
---|
594 | Returns the first reader object which worked.""" |
---|
595 | # Translated from OnImportGeneric method in GSASII.py |
---|
596 | primaryReaders, secondaryReaders = [], [] |
---|
597 | for reader in readerlist: |
---|
598 | flag = reader.ExtensionValidator(filename) |
---|
599 | if flag is None: |
---|
600 | secondaryReaders.append(reader) |
---|
601 | elif flag: |
---|
602 | primaryReaders.append(reader) |
---|
603 | if not secondaryReaders and not primaryReaders: |
---|
604 | raise G2ImportException("Could not read file: ", filename) |
---|
605 | |
---|
606 | with open(filename, 'Ur') as fp: |
---|
607 | rd_list = [] |
---|
608 | |
---|
609 | for rd in primaryReaders + secondaryReaders: |
---|
610 | # Initialize reader |
---|
611 | rd.selections = [] |
---|
612 | rd.dnames = [] |
---|
613 | rd.ReInitialize() |
---|
614 | # Rewind file |
---|
615 | fp.seek(0) |
---|
616 | rd.errors = "" |
---|
617 | if not rd.ContentsValidator(fp): |
---|
618 | # Report error |
---|
619 | pass |
---|
620 | if len(rd.selections) > 1: |
---|
621 | # Select data? |
---|
622 | # GSASII.py:543 |
---|
623 | raise G2ImportException("Not sure what data to select") |
---|
624 | |
---|
625 | block = 0 |
---|
626 | rdbuffer = {} |
---|
627 | fp.seek(0) |
---|
628 | repeat = True |
---|
629 | while repeat: |
---|
630 | repeat = False |
---|
631 | block += 1 |
---|
632 | rd.objname = os.path.basename(filename) |
---|
633 | try: |
---|
634 | flag = rd.Reader(filename, fp, buffer=rdbuffer, blocknum=block) |
---|
635 | except: |
---|
636 | flag = False |
---|
637 | if flag: |
---|
638 | # Omitting image loading special cases |
---|
639 | rd.readfilename = filename |
---|
640 | rd_list.append(copy.deepcopy(rd)) |
---|
641 | repeat = rd.repeat |
---|
642 | else: |
---|
643 | if GSASIIpath.GetConfigValue('debug'): print("{} Reader failed to read {}".format(rd.formatName,filename)) |
---|
644 | if rd_list: |
---|
645 | if rd.warnings: |
---|
646 | print("Read warning by", rd.formatName, "reader:", |
---|
647 | rd.warnings, file=sys.stderr) |
---|
648 | else: |
---|
649 | print("{} read by Reader {}\n".format(filename,rd.formatName)) |
---|
650 | return rd_list |
---|
651 | raise G2ImportException("No reader could read file: " + filename) |
---|
652 | |
---|
653 | |
---|
654 | def load_iprms(instfile, reader): |
---|
655 | """Loads instrument parameters from a file, and edits the |
---|
656 | given reader. |
---|
657 | |
---|
658 | Returns a 2-tuple of (Iparm1, Iparm2) parameters |
---|
659 | """ |
---|
660 | LoadG2fil() |
---|
661 | ext = os.path.splitext(instfile)[1] |
---|
662 | |
---|
663 | if ext.lower() == '.instprm': |
---|
664 | # New GSAS File, load appropriately |
---|
665 | with open(instfile) as f: |
---|
666 | lines = f.readlines() |
---|
667 | bank = reader.powderentry[2] |
---|
668 | numbanks = reader.numbanks |
---|
669 | iparms = G2fil.ReadPowderInstprm(lines, bank, numbanks, reader) |
---|
670 | |
---|
671 | reader.instfile = instfile |
---|
672 | reader.instmsg = 'GSAS-II file' + instfile |
---|
673 | return iparms |
---|
674 | elif ext.lower() not in ('.prm', '.inst', '.ins'): |
---|
675 | raise ValueError('Expected .prm file, found: ', instfile) |
---|
676 | |
---|
677 | # It's an old GSAS file, load appropriately |
---|
678 | Iparm = {} |
---|
679 | with open(instfile, 'Ur') as fp: |
---|
680 | for line in fp: |
---|
681 | if '#' in line: |
---|
682 | continue |
---|
683 | Iparm[line[:12]] = line[12:-1] |
---|
684 | ibanks = int(Iparm.get('INS BANK ', '1').strip()) |
---|
685 | if ibanks == 1: |
---|
686 | reader.instbank = 1 |
---|
687 | reader.powderentry[2] = 1 |
---|
688 | reader.instfile = instfile |
---|
689 | reader.instmsg = instfile + ' bank ' + str(reader.instbank) |
---|
690 | return G2fil.SetPowderInstParms(Iparm, reader) |
---|
691 | # TODO handle >1 banks |
---|
692 | raise NotImplementedError("Check GSASIIfiles.py: ReadPowderInstprm") |
---|
693 | |
---|
694 | def load_pwd_from_reader(reader, instprm, existingnames=[]): |
---|
695 | """Loads powder data from a reader object, and assembles it into a GSASII data tree. |
---|
696 | |
---|
697 | :returns: (name, tree) - 2-tuple of the histogram name (str), and data |
---|
698 | |
---|
699 | Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com) |
---|
700 | """ |
---|
701 | HistName = 'PWDR ' + G2obj.StripUnicode(reader.idstring, '_') |
---|
702 | HistName = unicode(G2obj.MakeUniqueLabel(HistName, existingnames)) |
---|
703 | |
---|
704 | try: |
---|
705 | Iparm1, Iparm2 = instprm |
---|
706 | except ValueError: |
---|
707 | Iparm1, Iparm2 = load_iprms(instprm, reader) |
---|
708 | |
---|
709 | Ymin = np.min(reader.powderdata[1]) |
---|
710 | Ymax = np.max(reader.powderdata[1]) |
---|
711 | valuesdict = {'wtFactor': 1.0, |
---|
712 | 'Dummy': False, |
---|
713 | 'ranId': ran.randint(0, sys.maxint), |
---|
714 | 'Offset': [0.0, 0.0], 'delOffset': 0.02*Ymax, |
---|
715 | 'refOffset': -0.1*Ymax, 'refDelt': 0.1*Ymax, |
---|
716 | 'qPlot': False, 'dPlot': False, 'sqrtPlot': False, |
---|
717 | 'Yminmax': [Ymin, Ymax]} |
---|
718 | reader.Sample['ranId'] = valuesdict['ranId'] |
---|
719 | |
---|
720 | # Ending keys: |
---|
721 | # [u'Reflection Lists', |
---|
722 | # u'Limits', |
---|
723 | # 'data', |
---|
724 | # u'Index Peak List', |
---|
725 | # u'Comments', |
---|
726 | # u'Unit Cells List', |
---|
727 | # u'Sample Parameters', |
---|
728 | # u'Peak List', |
---|
729 | # u'Background', |
---|
730 | # u'Instrument Parameters'] |
---|
731 | Tmin = np.min(reader.powderdata[0]) |
---|
732 | Tmax = np.max(reader.powderdata[0]) |
---|
733 | |
---|
734 | default_background = [['chebyschev', False, 3, 1.0, 0.0, 0.0], |
---|
735 | {'nDebye': 0, 'debyeTerms': [], 'nPeaks': 0, 'peaksList': []}] |
---|
736 | |
---|
737 | output_dict = {u'Reflection Lists': {}, |
---|
738 | u'Limits': reader.pwdparms.get('Limits', [(Tmin, Tmax), [Tmin, Tmax]]), |
---|
739 | u'data': [valuesdict, reader.powderdata, HistName], |
---|
740 | u'Index Peak List': [[], []], |
---|
741 | u'Comments': reader.comments, |
---|
742 | u'Unit Cells List': [], |
---|
743 | u'Sample Parameters': reader.Sample, |
---|
744 | u'Peak List': {'peaks': [], 'sigDict': {}}, |
---|
745 | u'Background': reader.pwdparms.get('Background', default_background), |
---|
746 | u'Instrument Parameters': [Iparm1, Iparm2], |
---|
747 | } |
---|
748 | |
---|
749 | names = [u'Comments', |
---|
750 | u'Limits', |
---|
751 | u'Background', |
---|
752 | u'Instrument Parameters', |
---|
753 | u'Sample Parameters', |
---|
754 | u'Peak List', |
---|
755 | u'Index Peak List', |
---|
756 | u'Unit Cells List', |
---|
757 | u'Reflection Lists'] |
---|
758 | |
---|
759 | # TODO controls?? GSASII.py:1664-7 |
---|
760 | |
---|
761 | return HistName, [HistName] + names, output_dict |
---|
762 | |
---|
763 | |
---|
764 | def _deep_copy_into(from_, into): |
---|
765 | """Helper function for reloading .gpx file. See G2Project.reload() |
---|
766 | |
---|
767 | :author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com) |
---|
768 | """ |
---|
769 | if isinstance(from_, dict) and isinstance(into, dict): |
---|
770 | combined_keys = set(from_.keys()).union(into.keys()) |
---|
771 | for key in combined_keys: |
---|
772 | if key in from_ and key in into: |
---|
773 | both_dicts = (isinstance(from_[key], dict) |
---|
774 | and isinstance(into[key], dict)) |
---|
775 | both_lists = (isinstance(from_[key], list) |
---|
776 | and isinstance(into[key], list)) |
---|
777 | if both_dicts or both_lists: |
---|
778 | _deep_copy_into(from_[key], into[key]) |
---|
779 | else: |
---|
780 | into[key] = from_[key] |
---|
781 | elif key in from_: |
---|
782 | into[key] = from_[key] |
---|
783 | else: # key in into |
---|
784 | del into[key] |
---|
785 | elif isinstance(from_, list) and isinstance(into, list): |
---|
786 | if len(from_) == len(into): |
---|
787 | for i in xrange(len(from_)): |
---|
788 | both_dicts = (isinstance(from_[i], dict) |
---|
789 | and isinstance(into[i], dict)) |
---|
790 | both_lists = (isinstance(from_[i], list) |
---|
791 | and isinstance(into[i], list)) |
---|
792 | if both_dicts or both_lists: |
---|
793 | _deep_copy_into(from_[i], into[i]) |
---|
794 | else: |
---|
795 | into[i] = from_[i] |
---|
796 | else: |
---|
797 | into[:] = from_ |
---|
798 | |
---|
799 | |
---|
800 | class G2ObjectWrapper(object): |
---|
801 | """Base class for all GSAS-II object wrappers. |
---|
802 | |
---|
803 | The underlying GSAS-II format can be accessed as `wrapper.data`. A number |
---|
804 | of overrides are implemented so that the wrapper behaves like a dictionary. |
---|
805 | |
---|
806 | Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com) |
---|
807 | """ |
---|
808 | def __init__(self, datadict): |
---|
809 | self.data = datadict |
---|
810 | |
---|
811 | def __getitem__(self, key): |
---|
812 | return self.data[key] |
---|
813 | |
---|
814 | def __setitem__(self, key, value): |
---|
815 | self.data[key] = value |
---|
816 | |
---|
817 | def __contains__(self, key): |
---|
818 | return key in self.data |
---|
819 | |
---|
820 | def get(self, k, d=None): |
---|
821 | return self.data.get(k, d) |
---|
822 | |
---|
823 | def keys(self): |
---|
824 | return self.data.keys() |
---|
825 | |
---|
826 | def values(self): |
---|
827 | return self.data.values() |
---|
828 | |
---|
829 | def items(self): |
---|
830 | return self.data.items() |
---|
831 | |
---|
832 | |
---|
833 | class G2Project(G2ObjectWrapper): |
---|
834 | """ |
---|
835 | Represents an entire GSAS-II project. |
---|
836 | |
---|
837 | There are two ways to initialize it: |
---|
838 | |
---|
839 | >>> # Load an existing project file |
---|
840 | >>> proj = G2Project('filename.gpx') |
---|
841 | >>> # Create a new project |
---|
842 | >>> proj = G2Project(filename='new_file.gpx') |
---|
843 | >>> # Specify an author |
---|
844 | >>> proj = G2Project(author='Dr. So-And-So', filename='my_data.gpx') |
---|
845 | |
---|
846 | Histograms can be accessed easily. |
---|
847 | |
---|
848 | >>> # By name |
---|
849 | >>> hist = proj.histogram('PWDR my-histogram-name') |
---|
850 | >>> # Or by index |
---|
851 | >>> hist = proj.histogram(0) |
---|
852 | >>> assert hist.id == 0 |
---|
853 | >>> # Or by random id |
---|
854 | >>> assert hist == proj.histogram(hist.ranId) |
---|
855 | |
---|
856 | Phases can be accessed the same way. |
---|
857 | |
---|
858 | >>> phase = proj.phase('name of phase') |
---|
859 | |
---|
860 | New data can also be loaded via :func:`~G2Project.add_phase` and |
---|
861 | :func:`~G2Project.add_powder_histogram`. |
---|
862 | |
---|
863 | >>> hist = proj.add_powder_histogram('some_data_file.chi', |
---|
864 | 'instrument_parameters.prm') |
---|
865 | >>> phase = proj.add_phase('my_phase.cif', histograms=[hist]) |
---|
866 | |
---|
867 | Parameters for Rietveld refinement can be turned on and off as well. |
---|
868 | See :func:`~G2Project.set_refinement`, :func:`~G2Project.refine`, |
---|
869 | :func:`~G2Project.iter_refinements`, :func:`~G2Project.do_refinements`. |
---|
870 | """ |
---|
871 | def __init__(self, gpxfile=None, author=None, filename=None): |
---|
872 | """Loads a GSAS-II project from a specified filename. |
---|
873 | |
---|
874 | :param str gpxfile: Existing .gpx file to be loaded. If nonexistent, |
---|
875 | creates an empty project. |
---|
876 | :param str author: Author's name. Optional. |
---|
877 | :param str filename: The filename the project should be saved to in |
---|
878 | the future. If both filename and gpxfile are present, the project is |
---|
879 | loaded from the gpxfile, then set to save to filename in the future""" |
---|
880 | if gpxfile is None: |
---|
881 | filename = os.path.abspath(os.path.expanduser(filename)) |
---|
882 | self.filename = filename |
---|
883 | self.data, self.names = make_empty_project(author=author, filename=filename) |
---|
884 | elif isinstance(gpxfile, str): |
---|
885 | # TODO set author, filename |
---|
886 | self.filename = os.path.abspath(os.path.expanduser(gpxfile)) |
---|
887 | self.data, self.names = LoadDictFromProjFile(gpxfile) |
---|
888 | self.update_ids() |
---|
889 | else: |
---|
890 | raise ValueError("Not sure what to do with gpxfile") |
---|
891 | |
---|
892 | @classmethod |
---|
893 | def from_dict_and_names(cls, gpxdict, names, filename=None): |
---|
894 | """Creates a :class:`G2Project` directly from |
---|
895 | a dictionary and a list of names. If in doubt, do not use this. |
---|
896 | |
---|
897 | :returns: a :class:`G2Project` |
---|
898 | """ |
---|
899 | out = cls() |
---|
900 | if filename: |
---|
901 | filename = os.path.abspath(os.path.expanduser(filename)) |
---|
902 | out.filename = filename |
---|
903 | gpxdict['Controls']['data']['LastSavedAs'] = filename |
---|
904 | else: |
---|
905 | try: |
---|
906 | out.filename = gpxdict['Controls']['data']['LastSavedAs'] |
---|
907 | except KeyError: |
---|
908 | out.filename = None |
---|
909 | out.data = gpxdict |
---|
910 | out.names = names |
---|
911 | |
---|
912 | def save(self, filename=None): |
---|
913 | """Saves the project, either to the current filename, or to a new file. |
---|
914 | |
---|
915 | Updates self.filename if a new filename provided""" |
---|
916 | # TODO update LastSavedUsing ? |
---|
917 | if filename: |
---|
918 | filename = os.path.abspath(os.path.expanduser(filename)) |
---|
919 | self.data['Controls']['data']['LastSavedAs'] = filename |
---|
920 | self.filename = filename |
---|
921 | elif not self.filename: |
---|
922 | raise AttributeError("No file name to save to") |
---|
923 | SaveDictToProjFile(self.data, self.names, self.filename) |
---|
924 | |
---|
925 | def add_powder_histogram(self, datafile, iparams, phases=[]): |
---|
926 | """Loads a powder data histogram into the project. |
---|
927 | |
---|
928 | Automatically checks for an instrument parameter file, or one can be |
---|
929 | provided. Note that in unix fashion, "~" can be used to indicate the |
---|
930 | home directory (e.g. ~/G2data/data.fxye). |
---|
931 | |
---|
932 | :param str datafile: The powder data file to read, a filename. |
---|
933 | :param str iparams: The instrument parameters file, a filename. |
---|
934 | :param list phases: Phases to link to the new histogram |
---|
935 | |
---|
936 | :returns: A :class:`G2PwdrData` object representing |
---|
937 | the histogram |
---|
938 | """ |
---|
939 | LoadG2fil() |
---|
940 | datafile = os.path.abspath(os.path.expanduser(datafile)) |
---|
941 | iparams = os.path.abspath(os.path.expanduser(iparams)) |
---|
942 | pwdrreaders = import_generic(datafile, PwdrDataReaders) |
---|
943 | histname, new_names, pwdrdata = load_pwd_from_reader( |
---|
944 | pwdrreaders[0], iparams, |
---|
945 | [h.name for h in self.histograms()]) |
---|
946 | if histname in self.data: |
---|
947 | print("Warning - redefining histogram", histname) |
---|
948 | elif self.names[-1][0] == 'Phases': |
---|
949 | self.names.insert(-1, new_names) |
---|
950 | else: |
---|
951 | self.names.append(new_names) |
---|
952 | self.data[histname] = pwdrdata |
---|
953 | self.update_ids() |
---|
954 | |
---|
955 | for phase in phases: |
---|
956 | phase = self.phase(phase) |
---|
957 | self.link_histogram_phase(histname, phase) |
---|
958 | |
---|
959 | return self.histogram(histname) |
---|
960 | |
---|
961 | def add_phase(self, phasefile, phasename=None, histograms=[]): |
---|
962 | """Loads a phase into the project from a .cif file |
---|
963 | |
---|
964 | :param str phasefile: The CIF file from which to import the phase. |
---|
965 | :param str phasename: The name of the new phase, or None for the default |
---|
966 | :param list histograms: The names of the histograms to associate with |
---|
967 | this phase |
---|
968 | |
---|
969 | :returns: A :class:`G2Phase` object representing the |
---|
970 | new phase. |
---|
971 | """ |
---|
972 | LoadG2fil() |
---|
973 | histograms = [self.histogram(h).name for h in histograms] |
---|
974 | phasefile = os.path.abspath(os.path.expanduser(phasefile)) |
---|
975 | |
---|
976 | # TODO handle multiple phases in a file |
---|
977 | phasereaders = import_generic(phasefile, PhaseReaders) |
---|
978 | phasereader = phasereaders[0] |
---|
979 | |
---|
980 | phasename = phasename or phasereader.Phase['General']['Name'] |
---|
981 | phaseNameList = [p.name for p in self.phases()] |
---|
982 | phasename = G2obj.MakeUniqueLabel(phasename, phaseNameList) |
---|
983 | phasereader.Phase['General']['Name'] = phasename |
---|
984 | |
---|
985 | if 'Phases' not in self.data: |
---|
986 | self.data[u'Phases'] = { 'data': None } |
---|
987 | assert phasename not in self.data['Phases'], "phase names should be unique" |
---|
988 | self.data['Phases'][phasename] = phasereader.Phase |
---|
989 | |
---|
990 | if phasereader.Constraints: |
---|
991 | Constraints = self.data['Constraints'] |
---|
992 | for i in phasereader.Constraints: |
---|
993 | if isinstance(i, dict): |
---|
994 | if '_Explain' not in Constraints: |
---|
995 | Constraints['_Explain'] = {} |
---|
996 | Constraints['_Explain'].update(i) |
---|
997 | else: |
---|
998 | Constraints['Phase'].append(i) |
---|
999 | |
---|
1000 | data = self.data['Phases'][phasename] |
---|
1001 | generalData = data['General'] |
---|
1002 | SGData = generalData['SGData'] |
---|
1003 | NShkl = len(G2spc.MustrainNames(SGData)) |
---|
1004 | NDij = len(G2spc.HStrainNames(SGData)) |
---|
1005 | Super = generalData.get('Super', 0) |
---|
1006 | if Super: |
---|
1007 | SuperVec = np.array(generalData['SuperVec'][0]) |
---|
1008 | else: |
---|
1009 | SuperVec = [] |
---|
1010 | UseList = data['Histograms'] |
---|
1011 | |
---|
1012 | for hist in histograms: |
---|
1013 | self.link_histogram_phase(hist, phasename) |
---|
1014 | |
---|
1015 | for obj in self.names: |
---|
1016 | if obj[0] == 'Phases': |
---|
1017 | phasenames = obj |
---|
1018 | break |
---|
1019 | else: |
---|
1020 | phasenames = [u'Phases'] |
---|
1021 | self.names.append(phasenames) |
---|
1022 | phasenames.append(unicode(phasename)) |
---|
1023 | |
---|
1024 | # TODO should it be self.filename, not phasefile? |
---|
1025 | SetupGeneral(data, os.path.dirname(phasefile)) |
---|
1026 | self.index_ids() |
---|
1027 | |
---|
1028 | self.update_ids() |
---|
1029 | return self.phase(phasename) |
---|
1030 | |
---|
1031 | def link_histogram_phase(self, histogram, phase): |
---|
1032 | """Associates a given histogram and phase. |
---|
1033 | |
---|
1034 | .. seealso:: |
---|
1035 | |
---|
1036 | :func:`~G2Project.histogram` |
---|
1037 | :func:`~G2Project.phase`""" |
---|
1038 | hist = self.histogram(histogram) |
---|
1039 | phase = self.phase(phase) |
---|
1040 | |
---|
1041 | generalData = phase['General'] |
---|
1042 | |
---|
1043 | if hist.name.startswith('HKLF '): |
---|
1044 | raise NotImplementedError("HKLF not yet supported") |
---|
1045 | elif hist.name.startswith('PWDR '): |
---|
1046 | hist['Reflection Lists'][generalData['Name']] = {} |
---|
1047 | UseList = phase['Histograms'] |
---|
1048 | SGData = generalData['SGData'] |
---|
1049 | NShkl = len(G2spc.MustrainNames(SGData)) |
---|
1050 | NDij = len(G2spc.HStrainNames(SGData)) |
---|
1051 | UseList[hist.name] = SetDefaultDData('PWDR', hist.name, NShkl=NShkl, NDij=NDij) |
---|
1052 | UseList[hist.name]['hId'] = hist.id |
---|
1053 | for key, val in [('Use', True), ('LeBail', False), |
---|
1054 | ('newLeBail', True), |
---|
1055 | ('Babinet', {'BabA': [0.0, False], |
---|
1056 | 'BabU': [0.0, False]})]: |
---|
1057 | if key not in UseList[hist.name]: |
---|
1058 | UseList[hist.name][key] = val |
---|
1059 | else: |
---|
1060 | raise RuntimeError("Unexpected histogram" + hist.name) |
---|
1061 | |
---|
1062 | |
---|
1063 | def reload(self): |
---|
1064 | """Reload self from self.filename""" |
---|
1065 | data, names = LoadDictFromProjFile(self.filename) |
---|
1066 | self.names = names |
---|
1067 | # Need to deep copy the new data file data into the current tree, |
---|
1068 | # so that any existing G2Phase, or G2PwdrData objects will still be |
---|
1069 | # valid |
---|
1070 | _deep_copy_into(from_=data, into=self.data) |
---|
1071 | |
---|
1072 | def refine(self, newfile=None, printFile=None, makeBack=False): |
---|
1073 | # TODO migrate to RefineCore |
---|
1074 | # G2strMain.RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList, |
---|
1075 | # calcControls,pawleyLookup,ifPrint,printFile,dlg) |
---|
1076 | # index_ids will automatically save the project |
---|
1077 | self.index_ids() |
---|
1078 | # TODO G2strMain does not properly use printFile |
---|
1079 | G2strMain.Refine(self.filename, makeBack=makeBack) |
---|
1080 | # Reload yourself |
---|
1081 | self.reload() |
---|
1082 | |
---|
1083 | def histogram(self, histname): |
---|
1084 | """Returns the histogram named histname, or None if it does not exist. |
---|
1085 | |
---|
1086 | :param histname: The name of the histogram (str), or ranId or index. |
---|
1087 | :returns: A :class:`G2PwdrData` object, or None if |
---|
1088 | the histogram does not exist |
---|
1089 | |
---|
1090 | .. seealso:: |
---|
1091 | :func:`~GSASIIscriptable.G2Project.histograms` |
---|
1092 | :func:`~GSASIIscriptable.G2Project.phase` |
---|
1093 | :func:`~GSASIIscriptable.G2Project.phases` |
---|
1094 | """ |
---|
1095 | if isinstance(histname, G2PwdrData): |
---|
1096 | if histname.proj == self: |
---|
1097 | return histname |
---|
1098 | if histname in self.data: |
---|
1099 | return G2PwdrData(self.data[histname], self) |
---|
1100 | try: |
---|
1101 | # see if histname is an id or ranId |
---|
1102 | histname = int(histname) |
---|
1103 | except ValueError: |
---|
1104 | return |
---|
1105 | |
---|
1106 | for histogram in self.histograms(): |
---|
1107 | if histogram.id == histname or histogram.ranId == histname: |
---|
1108 | return histogram |
---|
1109 | |
---|
1110 | def histograms(self): |
---|
1111 | """Return a list of all histograms, as |
---|
1112 | :class:`G2PwdrData` objects |
---|
1113 | |
---|
1114 | .. seealso:: |
---|
1115 | :func:`~GSASIIscriptable.G2Project.histograms` |
---|
1116 | :func:`~GSASIIscriptable.G2Project.phase` |
---|
1117 | :func:`~GSASIIscriptable.G2Project.phases` |
---|
1118 | """ |
---|
1119 | output = [] |
---|
1120 | for obj in self.names: |
---|
1121 | if len(obj) > 1 and obj[0] != u'Phases': |
---|
1122 | output.append(self.histogram(obj[0])) |
---|
1123 | return output |
---|
1124 | |
---|
1125 | def phase(self, phasename): |
---|
1126 | """ |
---|
1127 | Gives an object representing the specified phase in this project. |
---|
1128 | |
---|
1129 | :param str phasename: The name of the desired phase. Either the name |
---|
1130 | (str), the phase's ranId, or the phase's index |
---|
1131 | :returns: A :class:`G2Phase` object |
---|
1132 | :raises: KeyError |
---|
1133 | |
---|
1134 | .. seealso:: |
---|
1135 | :func:`~GSASIIscriptable.G2Project.histograms` |
---|
1136 | :func:`~GSASIIscriptable.G2Project.phase` |
---|
1137 | :func:`~GSASIIscriptable.G2Project.phases` |
---|
1138 | """ |
---|
1139 | phases = self.data['Phases'] |
---|
1140 | if phasename in phases: |
---|
1141 | return G2Phase(phases[phasename], phasename, self) |
---|
1142 | |
---|
1143 | try: |
---|
1144 | # phasename should be phase index or ranId |
---|
1145 | phasename = int(phasename) |
---|
1146 | except ValueError: |
---|
1147 | return |
---|
1148 | |
---|
1149 | for phase in self.phases(): |
---|
1150 | if phase.id == phasename or phase.ranId == phasename: |
---|
1151 | return phase |
---|
1152 | |
---|
1153 | def phases(self): |
---|
1154 | """ |
---|
1155 | Returns a list of all the phases in the project. |
---|
1156 | |
---|
1157 | :returns: A list of :class:`G2Phase` objects |
---|
1158 | |
---|
1159 | .. seealso:: |
---|
1160 | :func:`~GSASIIscriptable.G2Project.histogram` |
---|
1161 | :func:`~GSASIIscriptable.G2Project.histograms` |
---|
1162 | :func:`~GSASIIscriptable.G2Project.phase` |
---|
1163 | """ |
---|
1164 | for obj in self.names: |
---|
1165 | if obj[0] == 'Phases': |
---|
1166 | return [self.phase(p) for p in obj[1:]] |
---|
1167 | return [] |
---|
1168 | |
---|
1169 | def update_ids(self): |
---|
1170 | """Makes sure all phases and histograms have proper hId and pId""" |
---|
1171 | # Translated from GetUsedHistogramsAndPhasesfromTree, |
---|
1172 | # GSASIIdataGUI.py:4107 |
---|
1173 | for i, h in enumerate(self.histograms()): |
---|
1174 | h.id = i |
---|
1175 | for i, p in enumerate(self.phases()): |
---|
1176 | p.id = i |
---|
1177 | |
---|
1178 | def do_refinements(self, refinements, histogram='all', phase='all', |
---|
1179 | outputnames=None, makeBack=False): |
---|
1180 | """Conducts a series of refinements. Wrapper around iter_refinements |
---|
1181 | |
---|
1182 | :param list refinements: A list of dictionaries defining refinements |
---|
1183 | :param str histogram: Name of histogram for refinements to be applied |
---|
1184 | to, or 'all' |
---|
1185 | :param str phase: Name of phase for refinements to be applied to, or |
---|
1186 | 'all' |
---|
1187 | """ |
---|
1188 | for proj in self.iter_refinements(refinements, histogram, phase, |
---|
1189 | outputnames, makeBack): |
---|
1190 | pass |
---|
1191 | return self |
---|
1192 | |
---|
1193 | def iter_refinements(self, refinements, histogram='all', phase='all', |
---|
1194 | outputnames=None, makeBack=False): |
---|
1195 | """Conducts a series of refinements, iteratively. Stops after every |
---|
1196 | refinement and yields this project, to allow error checking or |
---|
1197 | logging of intermediate results. |
---|
1198 | |
---|
1199 | >>> def checked_refinements(proj): |
---|
1200 | ... for p in proj.iter_refinements(refs): |
---|
1201 | ... # Track intermediate results |
---|
1202 | ... log(p.histogram('0').residuals) |
---|
1203 | ... log(p.phase('0').get_cell()) |
---|
1204 | ... # Check if parameter diverged, nonsense answer, or whatever |
---|
1205 | ... if is_something_wrong(p): |
---|
1206 | ... raise Exception("I need a human!") |
---|
1207 | |
---|
1208 | :param list refinements: A list of dictionaries defining refinements |
---|
1209 | :param str histogram: Name of histogram for refinements to be applied |
---|
1210 | to, or 'all' |
---|
1211 | :param str phase: Name of phase for refinements to be applied to, or |
---|
1212 | 'all' |
---|
1213 | """ |
---|
1214 | if outputnames: |
---|
1215 | if len(refinements) != len(outputnames): |
---|
1216 | raise ValueError("Should have same number of outuputs to" |
---|
1217 | "refinements") |
---|
1218 | else: |
---|
1219 | outputnames = [None for r in refinements] |
---|
1220 | |
---|
1221 | for output, refinement in zip(outputnames, refinements): |
---|
1222 | self.set_refinement(refinement, histogram) |
---|
1223 | # Handle 'once' args - refinements that are disabled after this |
---|
1224 | # refinement |
---|
1225 | if 'once' in refinement: |
---|
1226 | temp = {'set': refinement['once']} |
---|
1227 | self.set_refinement(temp, histogram, phase) |
---|
1228 | |
---|
1229 | if output: |
---|
1230 | self.save(output) |
---|
1231 | |
---|
1232 | self.refine(makeBack=makeBack) |
---|
1233 | yield self |
---|
1234 | |
---|
1235 | # Handle 'once' args - refinements that are disabled after this |
---|
1236 | # refinement |
---|
1237 | if 'once' in refinement: |
---|
1238 | temp = {'clear': refinement['once']} |
---|
1239 | self.set_refinement(temp, histogram, phase) |
---|
1240 | |
---|
1241 | def set_refinement(self, refinement, histogram='all', phase='all'): |
---|
1242 | """Apply specified refinements to a given histogram(s) or phase(s). |
---|
1243 | |
---|
1244 | Refinement parameters are categorize in three groups: |
---|
1245 | |
---|
1246 | 1. Histogram parameters |
---|
1247 | 2. Phase parameters |
---|
1248 | 3. Histogram-and-Phase (HAP) parameters |
---|
1249 | |
---|
1250 | :param dict refinement: The refinements to be conducted |
---|
1251 | :param histogram: Either a name of a histogram (str), a list of |
---|
1252 | histogram names, or 'all' (default) |
---|
1253 | :param phase: Either a name of a phase (str), a list of phase names, or |
---|
1254 | 'all' (default) |
---|
1255 | |
---|
1256 | .. seealso:: |
---|
1257 | :func:`~G2PwdrData.set_refinements` |
---|
1258 | :func:`~G2PwdrData.clear_refinements` |
---|
1259 | :func:`~G2Phase.set_refinements` |
---|
1260 | :func:`~G2Phase.clear_refinements` |
---|
1261 | :func:`~G2Phase.set_HAP_refinements` |
---|
1262 | :func:`~G2Phase.clear_HAP_refinements`""" |
---|
1263 | |
---|
1264 | if histogram == 'all': |
---|
1265 | hists = self.histograms() |
---|
1266 | elif isinstance(histogram, str) or isinstance(histogram, int): |
---|
1267 | hists = [self.histogram(histogram)] |
---|
1268 | else: |
---|
1269 | hists = [self.histogram(name) for name in histogram] |
---|
1270 | |
---|
1271 | if phase == 'all': |
---|
1272 | phases = self.phases() |
---|
1273 | elif isinstance(phase, str) or isinstance(phase, int): |
---|
1274 | phases = [self.phase(phase)] |
---|
1275 | else: |
---|
1276 | phases = [self.phase(name) for name in phase] |
---|
1277 | |
---|
1278 | # TODO: HAP parameters: |
---|
1279 | # Babinet |
---|
1280 | # Extinction |
---|
1281 | # HStrain |
---|
1282 | # Mustrain |
---|
1283 | # Pref. Ori |
---|
1284 | # Size |
---|
1285 | |
---|
1286 | pwdr_set = {} |
---|
1287 | phase_set = {} |
---|
1288 | hap_set = {} |
---|
1289 | for key, val in refinement.get('set', {}).items(): |
---|
1290 | # Apply refinement options |
---|
1291 | if G2PwdrData.is_valid_refinement_key(key): |
---|
1292 | pwdr_set[key] = val |
---|
1293 | elif G2Phase.is_valid_refinement_key(key): |
---|
1294 | phase_set[key] = val |
---|
1295 | elif G2Phase.is_valid_HAP_refinement_key(key): |
---|
1296 | hap_set[key] = val |
---|
1297 | else: |
---|
1298 | raise ValueError("Unknown refinement key", key) |
---|
1299 | |
---|
1300 | for hist in hists: |
---|
1301 | hist.set_refinements(pwdr_set) |
---|
1302 | for phase in phases: |
---|
1303 | phase.set_refinements(phase_set) |
---|
1304 | for phase in phases: |
---|
1305 | phase.set_HAP_refinements(hap_set, hists) |
---|
1306 | |
---|
1307 | pwdr_clear = {} |
---|
1308 | phase_clear = {} |
---|
1309 | hap_clear = {} |
---|
1310 | for key, val in refinement.get('clear', {}).items(): |
---|
1311 | # Clear refinement options |
---|
1312 | if G2PwdrData.is_valid_refinement_key(key): |
---|
1313 | pwdr_clear[key] = val |
---|
1314 | elif G2Phase.is_valid_refinement_key(key): |
---|
1315 | phase_clear[key] = val |
---|
1316 | elif G2Phase.is_valid_HAP_refinement_key(key): |
---|
1317 | hap_set[key] = val |
---|
1318 | else: |
---|
1319 | raise ValueError("Unknown refinement key", key) |
---|
1320 | |
---|
1321 | for hist in hists: |
---|
1322 | hist.clear_refinements(pwdr_clear) |
---|
1323 | for phase in phases: |
---|
1324 | phase.clear_refinements(phase_clear) |
---|
1325 | for phase in phases: |
---|
1326 | phase.clear_HAP_refinements(hap_clear, hists) |
---|
1327 | |
---|
1328 | def index_ids(self): |
---|
1329 | import GSASIIstrIO as G2strIO |
---|
1330 | self.save() |
---|
1331 | return G2strIO.GetUsedHistogramsAndPhases(self.filename) |
---|
1332 | |
---|
1333 | def add_constraint_raw(self, cons_scope, constr): |
---|
1334 | """Adds a constraint of type consType to the project. |
---|
1335 | cons_scope should be one of "Hist", "Phase", "HAP", or "Global". |
---|
1336 | |
---|
1337 | WARNING it does not check the constraint is well-constructed""" |
---|
1338 | constrs = self.data['Constraints']['data'] |
---|
1339 | if 'Global' not in constrs: |
---|
1340 | constrs['Global'] = [] |
---|
1341 | constrs[cons_scope].append(constr) |
---|
1342 | |
---|
1343 | def hold_many(self, vars, type): |
---|
1344 | """Apply holds for all the variables in vars, for constraint of a given type. |
---|
1345 | |
---|
1346 | type is passed directly to add_constraint_raw as consType |
---|
1347 | |
---|
1348 | :param list vars: A list of variables to hold. Either :class:`GSASIIobj.G2VarObj` objects, |
---|
1349 | string variable specifiers, or arguments for :meth:`make_var_obj` |
---|
1350 | :param str type: A string constraint type specifier. See |
---|
1351 | :class:`~GSASIIscriptable.G2Project.add_constraint_raw` |
---|
1352 | |
---|
1353 | """ |
---|
1354 | for var in vars: |
---|
1355 | if isinstance(var, str): |
---|
1356 | var = self.make_var_obj(var) |
---|
1357 | elif not isinstance(var, G2obj.G2VarObj): |
---|
1358 | var = self.make_var_obj(*var) |
---|
1359 | self.add_constraint_raw(type, [[1.0, var], None, None, 'h']) |
---|
1360 | |
---|
1361 | def make_var_obj(self, phase=None, hist=None, varname=None, atomId=None, |
---|
1362 | reloadIdx=True): |
---|
1363 | """Wrapper to create a G2VarObj. Takes either a string representaiton ("p:h:name:a") |
---|
1364 | or individual names of phase, histogram, varname, and atomId. |
---|
1365 | |
---|
1366 | Automatically converts string phase, hist, or atom names into the ID required |
---|
1367 | by G2VarObj.""" |
---|
1368 | |
---|
1369 | if reloadIdx: |
---|
1370 | self.index_ids() |
---|
1371 | |
---|
1372 | # If string representation, short circuit |
---|
1373 | if hist is None and varname is None and atomId is None: |
---|
1374 | if isinstance(phase, str) and ':' in phase: |
---|
1375 | return G2obj.G2VarObj(phase) |
---|
1376 | |
---|
1377 | # Get phase index |
---|
1378 | phaseObj = None |
---|
1379 | if isinstance(phase, G2Phase): |
---|
1380 | phaseObj = phase |
---|
1381 | phase = G2obj.PhaseRanIdLookup[phase.ranId] |
---|
1382 | elif phase in self.data['Phases']: |
---|
1383 | phaseObj = self.phase(phase) |
---|
1384 | phaseRanId = phaseObj.ranId |
---|
1385 | phase = G2obj.PhaseRanIdLookup[phaseRanId] |
---|
1386 | |
---|
1387 | # Get histogram index |
---|
1388 | if isinstance(hist, G2PwdrData): |
---|
1389 | hist = G2obj.HistRanIdLookup[hist.ranId] |
---|
1390 | elif hist in self.data: |
---|
1391 | histRanId = self.histogram(hist).ranId |
---|
1392 | hist = G2obj.HistRanIdLookup[histRanId] |
---|
1393 | |
---|
1394 | # Get atom index (if any) |
---|
1395 | if isinstance(atomId, G2AtomRecord): |
---|
1396 | atomId = G2obj.AtomRanIdLookup[phase][atomId.ranId] |
---|
1397 | elif phaseObj: |
---|
1398 | atomObj = phaseObj.atom(atomId) |
---|
1399 | if atomObj: |
---|
1400 | atomRanId = atomObj.ranId |
---|
1401 | atomId = G2obj.AtomRanIdLookup[phase][atomRanId] |
---|
1402 | |
---|
1403 | return G2obj.G2VarObj(phase, hist, varname, atomId) |
---|
1404 | |
---|
1405 | |
---|
1406 | class G2AtomRecord(G2ObjectWrapper): |
---|
1407 | """Wrapper for an atom record. Has convenient accessors via @property. |
---|
1408 | |
---|
1409 | |
---|
1410 | Available accessors: label, type, refinement_flags, coordinates, |
---|
1411 | occupancy, ranId, id, adp_flag, uiso |
---|
1412 | |
---|
1413 | Example: |
---|
1414 | |
---|
1415 | >>> atom = some_phase.atom("O3") |
---|
1416 | >>> # We can access the underlying data format |
---|
1417 | >>> atom.data |
---|
1418 | ['O3', 'O-2', '', ... ] |
---|
1419 | >>> # We can also use wrapper accessors |
---|
1420 | >>> atom.coordinates |
---|
1421 | (0.33, 0.15, 0.5) |
---|
1422 | >>> atom.refinement_flags |
---|
1423 | u'FX' |
---|
1424 | >>> atom.ranId |
---|
1425 | 4615973324315876477 |
---|
1426 | >>> atom.occupancy |
---|
1427 | 1.0 |
---|
1428 | """ |
---|
1429 | def __init__(self, data, indices, proj): |
---|
1430 | self.data = data |
---|
1431 | self.cx, self.ct, self.cs, self.cia = indices |
---|
1432 | self.proj = proj |
---|
1433 | |
---|
1434 | @property |
---|
1435 | def label(self): |
---|
1436 | return self.data[self.ct-1] |
---|
1437 | |
---|
1438 | @property |
---|
1439 | def type(self): |
---|
1440 | return self.data[self.ct] |
---|
1441 | |
---|
1442 | @property |
---|
1443 | def refinement_flags(self): |
---|
1444 | return self.data[self.ct+1] |
---|
1445 | |
---|
1446 | @refinement_flags.setter |
---|
1447 | def refinement_flags(self, other): |
---|
1448 | # Automatically check it is a valid refinement |
---|
1449 | for c in other: |
---|
1450 | if c not in ' FXU': |
---|
1451 | raise ValueError("Invalid atom refinement: ", other) |
---|
1452 | self.data[self.ct+1] = unicode(other) |
---|
1453 | |
---|
1454 | @property |
---|
1455 | def coordinates(self): |
---|
1456 | return tuple(self.data[self.cx:self.cx+3]) |
---|
1457 | |
---|
1458 | @property |
---|
1459 | def occupancy(self): |
---|
1460 | return self.data[self.cx+3] |
---|
1461 | |
---|
1462 | @occupancy.setter |
---|
1463 | def occupancy(self, val): |
---|
1464 | self.data[self.cx+3] = float(val) |
---|
1465 | |
---|
1466 | @property |
---|
1467 | def ranId(self): |
---|
1468 | return self.data[self.cia+8] |
---|
1469 | |
---|
1470 | @property |
---|
1471 | def adp_flag(self): |
---|
1472 | # Either 'I' or 'A' |
---|
1473 | return self.data[self.cia] |
---|
1474 | |
---|
1475 | @property |
---|
1476 | def uiso(self): |
---|
1477 | if self.adp_flag == 'I': |
---|
1478 | return self.data[self.cia+1] |
---|
1479 | else: |
---|
1480 | return self.data[self.cia+2:self.cia+8] |
---|
1481 | |
---|
1482 | @uiso.setter |
---|
1483 | def uiso(self, value): |
---|
1484 | if self.adp_flag == 'I': |
---|
1485 | self.data[self.cia+1] = float(value) |
---|
1486 | else: |
---|
1487 | assert len(value) == 6 |
---|
1488 | self.data[self.cia+2:self.cia+8] = [float(v) for v in value] |
---|
1489 | |
---|
1490 | |
---|
1491 | class G2PwdrData(G2ObjectWrapper): |
---|
1492 | """Wraps a Powder Data Histogram.""" |
---|
1493 | def __init__(self, data, proj): |
---|
1494 | self.data = data |
---|
1495 | self.proj = proj |
---|
1496 | |
---|
1497 | @staticmethod |
---|
1498 | def is_valid_refinement_key(key): |
---|
1499 | valid_keys = ['Limits', 'Sample Parameters', 'Background', |
---|
1500 | 'Instrument Parameters'] |
---|
1501 | return key in valid_keys |
---|
1502 | |
---|
1503 | @property |
---|
1504 | def name(self): |
---|
1505 | return self.data['data'][-1] |
---|
1506 | |
---|
1507 | @property |
---|
1508 | def ranId(self): |
---|
1509 | return self.data['data'][0]['ranId'] |
---|
1510 | |
---|
1511 | @property |
---|
1512 | def residuals(self): |
---|
1513 | data = self.data['data'][0] |
---|
1514 | return {key: data[key] |
---|
1515 | for key in ['R', 'Rb', 'wR', 'wRb', 'wRmin']} |
---|
1516 | |
---|
1517 | @property |
---|
1518 | def id(self): |
---|
1519 | self.proj.update_ids() |
---|
1520 | return self.data['data'][0]['hId'] |
---|
1521 | |
---|
1522 | @id.setter |
---|
1523 | def id(self, val): |
---|
1524 | self.data['data'][0]['hId'] = val |
---|
1525 | |
---|
1526 | def fit_fixed_points(self): |
---|
1527 | """Attempts to apply a background fit to the fixed points currently specified.""" |
---|
1528 | def SetInstParms(Inst): |
---|
1529 | dataType = Inst['Type'][0] |
---|
1530 | insVary = [] |
---|
1531 | insNames = [] |
---|
1532 | insVals = [] |
---|
1533 | for parm in Inst: |
---|
1534 | insNames.append(parm) |
---|
1535 | insVals.append(Inst[parm][1]) |
---|
1536 | if parm in ['U','V','W','X','Y','SH/L','I(L2)/I(L1)','alpha', |
---|
1537 | 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',] and Inst[parm][2]: |
---|
1538 | Inst[parm][2] = False |
---|
1539 | instDict = dict(zip(insNames, insVals)) |
---|
1540 | instDict['X'] = max(instDict['X'], 0.01) |
---|
1541 | instDict['Y'] = max(instDict['Y'], 0.01) |
---|
1542 | if 'SH/L' in instDict: |
---|
1543 | instDict['SH/L'] = max(instDict['SH/L'], 0.002) |
---|
1544 | return dataType, instDict, insVary |
---|
1545 | |
---|
1546 | bgrnd = self.data['Background'] |
---|
1547 | |
---|
1548 | # Need our fixed points in order |
---|
1549 | bgrnd[1]['FixedPoints'].sort(key=lambda pair: pair[0]) |
---|
1550 | X = [x for x, y in bgrnd[1]['FixedPoints']] |
---|
1551 | Y = [y for x, y in bgrnd[1]['FixedPoints']] |
---|
1552 | |
---|
1553 | limits = self.data['Limits'][1] |
---|
1554 | if X[0] > limits[0]: |
---|
1555 | X = [limits[0]] + X |
---|
1556 | Y = [Y[0]] + Y |
---|
1557 | if X[-1] < limits[1]: |
---|
1558 | X += [limits[1]] |
---|
1559 | Y += [Y[-1]] |
---|
1560 | |
---|
1561 | # Some simple lookups |
---|
1562 | controls = self.proj['Controls']['data'] |
---|
1563 | inst, inst2 = self.data['Instrument Parameters'] |
---|
1564 | pwddata = self.data['data'][1] |
---|
1565 | |
---|
1566 | # Construct the data for background fitting |
---|
1567 | xBeg = np.searchsorted(pwddata[0], limits[0]) |
---|
1568 | xFin = np.searchsorted(pwddata[0], limits[1]) |
---|
1569 | xdata = pwddata[0][xBeg:xFin] |
---|
1570 | ydata = si.interp1d(X,Y)(ma.getdata(xdata)) |
---|
1571 | |
---|
1572 | W = [1]*len(xdata) |
---|
1573 | Z = [0]*len(xdata) |
---|
1574 | |
---|
1575 | dataType, insDict, insVary = SetInstParms(inst) |
---|
1576 | bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd) |
---|
1577 | |
---|
1578 | # Do the fit |
---|
1579 | data = np.array([xdata, ydata, W, Z, Z, Z]) |
---|
1580 | G2pwd.DoPeakFit('LSQ', [], bgrnd, limits, inst, inst2, data, |
---|
1581 | prevVaryList=bakVary, controls=controls) |
---|
1582 | |
---|
1583 | # Post-fit |
---|
1584 | parmDict = {} |
---|
1585 | bakType, bakDict, bakVary = G2pwd.SetBackgroundParms(bgrnd) |
---|
1586 | parmDict.update(bakDict) |
---|
1587 | parmDict.update(insDict) |
---|
1588 | pwddata[3][xBeg:xFin] *= 0 |
---|
1589 | pwddata[5][xBeg:xFin] *= 0 |
---|
1590 | pwddata[4][xBeg:xFin] = G2pwd.getBackground('', parmDict, bakType, dataType, xdata)[0] |
---|
1591 | |
---|
1592 | # TODO adjust pwddata? GSASIIpwdGUI.py:1041 |
---|
1593 | # TODO update background |
---|
1594 | self.proj.save() |
---|
1595 | |
---|
1596 | def y_calc(self): |
---|
1597 | return self.data['data'][1][3] |
---|
1598 | |
---|
1599 | def plot(self, Yobs=True, Ycalc=True, Background=True, Residual=True): |
---|
1600 | try: |
---|
1601 | import matplotlib.pyplot as plt |
---|
1602 | data = self.data['data'][1] |
---|
1603 | if Yobs: |
---|
1604 | plt.plot(data[0], data[1], label='Yobs') |
---|
1605 | if Ycalc: |
---|
1606 | plt.plot(data[0], data[3], label='Ycalc') |
---|
1607 | if Background: |
---|
1608 | plt.plot(data[0], data[4], label='Background') |
---|
1609 | if Residual: |
---|
1610 | plt.plot(data[0], data[5], label="Residual") |
---|
1611 | except ImportError: |
---|
1612 | pass |
---|
1613 | |
---|
1614 | def set_refinements(self, refs): |
---|
1615 | """Sets the refinement parameter 'key' to the specification 'value' |
---|
1616 | |
---|
1617 | :param dict refs: A dictionary of the parameters to be set. See |
---|
1618 | :ref:`Histogram_parameters_table` for a description of |
---|
1619 | what these dictionaries should be. |
---|
1620 | |
---|
1621 | :returns: None |
---|
1622 | |
---|
1623 | """ |
---|
1624 | do_fit_fixed_points = False |
---|
1625 | for key, value in refs.items(): |
---|
1626 | if key == 'Limits': |
---|
1627 | old_limits = self.data['Limits'][1] |
---|
1628 | new_limits = value |
---|
1629 | if isinstance(new_limits, dict): |
---|
1630 | if 'low' in new_limits: |
---|
1631 | old_limits[0] = new_limits['low'] |
---|
1632 | if 'high' in new_limits: |
---|
1633 | old_limits[1] = new_limits['high'] |
---|
1634 | else: |
---|
1635 | old_limits[0], old_limits[1] = new_limits |
---|
1636 | elif key == 'Sample Parameters': |
---|
1637 | sample = self.data['Sample Parameters'] |
---|
1638 | for sparam in value: |
---|
1639 | if sparam not in sample: |
---|
1640 | raise ValueError("Unknown refinement parameter, " |
---|
1641 | + str(sparam)) |
---|
1642 | sample[sparam][1] = True |
---|
1643 | elif key == 'Background': |
---|
1644 | bkg, peaks = self.data['Background'] |
---|
1645 | |
---|
1646 | # If True or False, just set the refine parameter |
---|
1647 | if value in (True, False): |
---|
1648 | bkg[1] = value |
---|
1649 | return |
---|
1650 | |
---|
1651 | if 'type' in value: |
---|
1652 | bkg[0] = value['type'] |
---|
1653 | if 'refine' in value: |
---|
1654 | bkg[1] = value['refine'] |
---|
1655 | if 'no. coeffs' in value: |
---|
1656 | cur_coeffs = bkg[2] |
---|
1657 | n_coeffs = value['no. coeffs'] |
---|
1658 | if n_coeffs > cur_coeffs: |
---|
1659 | for x in range(n_coeffs - cur_coeffs): |
---|
1660 | bkg.append(0.0) |
---|
1661 | else: |
---|
1662 | for _ in range(cur_coeffs - n_coeffs): |
---|
1663 | bkg.pop() |
---|
1664 | bkg[2] = n_coeffs |
---|
1665 | if 'coeffs' in value: |
---|
1666 | bkg[3:] = value['coeffs'] |
---|
1667 | if 'FixedPoints' in value: |
---|
1668 | peaks['FixedPoints'] = [(float(a), float(b)) |
---|
1669 | for a, b in value['FixedPoints']] |
---|
1670 | if value.get('fit fixed points', False): |
---|
1671 | do_fit_fixed_points = True |
---|
1672 | |
---|
1673 | elif key == 'Instrument Parameters': |
---|
1674 | instrument, secondary = self.data['Instrument Parameters'] |
---|
1675 | for iparam in value: |
---|
1676 | try: |
---|
1677 | instrument[iparam][2] = True |
---|
1678 | except IndexError: |
---|
1679 | raise ValueError("Invalid key:", iparam) |
---|
1680 | else: |
---|
1681 | raise ValueError("Unknown key:", key) |
---|
1682 | # Fit fixed points after the fact - ensure they are after fixed points |
---|
1683 | # are added |
---|
1684 | if do_fit_fixed_points: |
---|
1685 | # Background won't be fit if refinement flag not set |
---|
1686 | orig = self.data['Background'][0][1] |
---|
1687 | self.data['Background'][0][1] = True |
---|
1688 | self.fit_fixed_points() |
---|
1689 | # Restore the previous value |
---|
1690 | self.data['Background'][0][1] = orig |
---|
1691 | |
---|
1692 | def clear_refinements(self, refs): |
---|
1693 | """Clears the refinement parameter 'key' and its associated value. |
---|
1694 | |
---|
1695 | :param dict refs: A dictionary of parameters to clear.""" |
---|
1696 | for key, value in refs.items(): |
---|
1697 | if key == 'Limits': |
---|
1698 | old_limits, cur_limits = self.data['Limits'] |
---|
1699 | cur_limits[0], cur_limits[1] = old_limits |
---|
1700 | elif key == 'Sample Parameters': |
---|
1701 | sample = self.data['Sample Parameters'] |
---|
1702 | for sparam in value: |
---|
1703 | sample[sparam][1] = False |
---|
1704 | elif key == 'Background': |
---|
1705 | bkg, peaks = self.data['Background'] |
---|
1706 | |
---|
1707 | # If True or False, just set the refine parameter |
---|
1708 | if value in (True, False): |
---|
1709 | bkg[1] = False |
---|
1710 | return |
---|
1711 | |
---|
1712 | bkg[1] = False |
---|
1713 | if 'FixedPoints' in value: |
---|
1714 | if 'FixedPoints' in peaks: |
---|
1715 | del peaks['FixedPoints'] |
---|
1716 | elif key == 'Instrument Parameters': |
---|
1717 | instrument, secondary = self.data['Instrument Parameters'] |
---|
1718 | for iparam in value: |
---|
1719 | instrument[iparam][2] = False |
---|
1720 | else: |
---|
1721 | raise ValueError("Unknown key:", key) |
---|
1722 | |
---|
1723 | |
---|
1724 | class G2Phase(G2ObjectWrapper): |
---|
1725 | """A wrapper object around a given phase. |
---|
1726 | |
---|
1727 | Author: Jackson O'Donnell (jacksonhodonnell .at. gmail.com) |
---|
1728 | """ |
---|
1729 | def __init__(self, data, name, proj): |
---|
1730 | self.data = data |
---|
1731 | self.name = name |
---|
1732 | self.proj = proj |
---|
1733 | |
---|
1734 | @staticmethod |
---|
1735 | def is_valid_refinement_key(key): |
---|
1736 | valid_keys = ["Cell", "Atoms", "LeBail"] |
---|
1737 | return key in valid_keys |
---|
1738 | |
---|
1739 | @staticmethod |
---|
1740 | def is_valid_HAP_refinement_key(key): |
---|
1741 | valid_keys = ["Babinet", "Extinction", "HStrain", "Mustrain", |
---|
1742 | "Pref.Ori.", "Show", "Size", "Use", "Scale"] |
---|
1743 | return key in valid_keys |
---|
1744 | |
---|
1745 | def atom(self, atomlabel): |
---|
1746 | """Returns the atom specified by atomlabel, or None if it does not |
---|
1747 | exist. |
---|
1748 | |
---|
1749 | :param str atomlabel: The name of the atom (e.g. "O2") |
---|
1750 | :returns: A :class:`G2AtomRecord` object |
---|
1751 | representing the atom. |
---|
1752 | """ |
---|
1753 | # Consult GSASIIobj.py for the meaning of this |
---|
1754 | cx, ct, cs, cia = self.data['General']['AtomPtrs'] |
---|
1755 | ptrs = [cx, ct, cs, cia] |
---|
1756 | atoms = self.data['Atoms'] |
---|
1757 | for atom in atoms: |
---|
1758 | if atom[ct-1] == atomlabel: |
---|
1759 | return G2AtomRecord(atom, ptrs, self.proj) |
---|
1760 | |
---|
1761 | def atoms(self): |
---|
1762 | """Returns a list of atoms present in the phase. |
---|
1763 | |
---|
1764 | :returns: A list of :class:`G2AtomRecord` objects. |
---|
1765 | |
---|
1766 | .. seealso:: |
---|
1767 | :func:`~GSASIIscriptable.G2Phase.atom` |
---|
1768 | :class:`G2AtomRecord` |
---|
1769 | """ |
---|
1770 | ptrs = self.data['General']['AtomPtrs'] |
---|
1771 | output = [] |
---|
1772 | atoms = self.data['Atoms'] |
---|
1773 | for atom in atoms: |
---|
1774 | output.append(G2AtomRecord(atom, ptrs, self.proj)) |
---|
1775 | return output |
---|
1776 | |
---|
1777 | def histograms(self): |
---|
1778 | output = [] |
---|
1779 | for hname in self.data.get('Histograms', {}).keys(): |
---|
1780 | output.append(self.proj.histogram(hname)) |
---|
1781 | return output |
---|
1782 | |
---|
1783 | @property |
---|
1784 | def ranId(self): |
---|
1785 | return self.data['ranId'] |
---|
1786 | |
---|
1787 | @property |
---|
1788 | def id(self): |
---|
1789 | return self.data['pId'] |
---|
1790 | |
---|
1791 | @id.setter |
---|
1792 | def id(self, val): |
---|
1793 | self.data['pId'] = val |
---|
1794 | |
---|
1795 | def get_cell(self): |
---|
1796 | """Returns a dictionary of the cell parameters, with keys: |
---|
1797 | 'length_a', 'length_b', 'length_c', 'angle_alpha', 'angle_beta', 'angle_gamma', 'volume' |
---|
1798 | |
---|
1799 | :returns: a dict |
---|
1800 | |
---|
1801 | .. seealso:: |
---|
1802 | :func:`~G2Phase.get_cell_and_esd` |
---|
1803 | |
---|
1804 | """ |
---|
1805 | cell = self.data['General']['Cell'] |
---|
1806 | return {'length_a': cell[1], 'length_b': cell[2], 'length_c': cell[3], |
---|
1807 | 'angle_alpha': cell[4], 'angle_beta': cell[5], 'angle_gamma': cell[6], |
---|
1808 | 'volume': cell[7]} |
---|
1809 | |
---|
1810 | def get_cell_and_esd(self): |
---|
1811 | """ |
---|
1812 | Returns a pair of dictionaries, the first representing the unit cell, the second |
---|
1813 | representing the estimated standard deviations of the unit cell. |
---|
1814 | |
---|
1815 | :returns: a tuple of two dictionaries |
---|
1816 | |
---|
1817 | .. seealso:: |
---|
1818 | :func:`~G2Phase.get_cell` |
---|
1819 | |
---|
1820 | """ |
---|
1821 | # translated from GSASIIstrIO.ExportBaseclass.GetCell |
---|
1822 | import GSASIIstrIO as G2stIO |
---|
1823 | import GSASIIlattice as G2lat |
---|
1824 | import GSASIImapvars as G2mv |
---|
1825 | try: |
---|
1826 | pfx = str(self.id) + '::' |
---|
1827 | sgdata = self['General']['SGData'] |
---|
1828 | covDict = self.proj['Covariance']['data'] |
---|
1829 | |
---|
1830 | parmDict = dict(zip(covDict.get('varyList',[]), |
---|
1831 | covDict.get('variables',[]))) |
---|
1832 | sigDict = dict(zip(covDict.get('varyList',[]), |
---|
1833 | covDict.get('sig',[]))) |
---|
1834 | |
---|
1835 | if covDict.get('covMatrix') is not None: |
---|
1836 | sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'], |
---|
1837 | covDict['varyList'], |
---|
1838 | parmDict)) |
---|
1839 | |
---|
1840 | A, sigA = G2stIO.cellFill(pfx, sgdata, parmDict, sigDict) |
---|
1841 | cellSig = G2stIO.getCellEsd(pfx, sgdata, A, self.proj['Covariance']['data']) |
---|
1842 | cellList = G2lat.A2cell(A) + (G2lat.calc_V(A),) |
---|
1843 | cellDict, cellSigDict = {}, {} |
---|
1844 | for i, key in enumerate(['length_a', 'length_b', 'length_c', |
---|
1845 | 'angle_alpha', 'angle_beta', 'angle_gamma', |
---|
1846 | 'volume']): |
---|
1847 | cellDict[key] = cellList[i] |
---|
1848 | cellSigDict[key] = cellSig[i] |
---|
1849 | return cellDict, cellSigDict |
---|
1850 | except KeyError: |
---|
1851 | cell = self.get_cell() |
---|
1852 | return cell, {key: 0.0 for key in cell} |
---|
1853 | |
---|
1854 | def export_CIF(self, outputname, quickmode=True): |
---|
1855 | """Write this phase to a .cif file named outputname |
---|
1856 | |
---|
1857 | :param str outputname: The name of the .cif file to write to |
---|
1858 | :param bool quickmode: Currently ignored. Carryover from exports.G2export_CIF""" |
---|
1859 | # This code is all taken from exports/G2export_CIF.py |
---|
1860 | # Functions copied have the same names |
---|
1861 | import GSASIImath as G2mth |
---|
1862 | import GSASIImapvars as G2mv |
---|
1863 | from exports import G2export_CIF as cif |
---|
1864 | |
---|
1865 | CIFdate = dt.datetime.strftime(dt.datetime.now(),"%Y-%m-%dT%H:%M") |
---|
1866 | CIFname = os.path.splitext(self.proj.filename)[0] |
---|
1867 | CIFname = os.path.split(CIFname)[1] |
---|
1868 | CIFname = ''.join([c if ord(c) < 128 else '' |
---|
1869 | for c in CIFname.replace(' ', '_')]) |
---|
1870 | try: |
---|
1871 | author = self.proj['Controls']['data'].get('Author','').strip() |
---|
1872 | except KeyError: |
---|
1873 | pass |
---|
1874 | oneblock = True |
---|
1875 | |
---|
1876 | covDict = self.proj['Covariance']['data'] |
---|
1877 | parmDict = dict(zip(covDict.get('varyList',[]), |
---|
1878 | covDict.get('variables',[]))) |
---|
1879 | sigDict = dict(zip(covDict.get('varyList',[]), |
---|
1880 | covDict.get('sig',[]))) |
---|
1881 | |
---|
1882 | if covDict.get('covMatrix') is not None: |
---|
1883 | sigDict.update(G2mv.ComputeDepESD(covDict['covMatrix'], |
---|
1884 | covDict['varyList'], |
---|
1885 | parmDict)) |
---|
1886 | |
---|
1887 | with open(outputname, 'w') as fp: |
---|
1888 | fp.write(' \n' + 70*'#' + '\n') |
---|
1889 | cif.WriteCIFitem(fp, 'data_' + CIFname) |
---|
1890 | # from exports.G2export_CIF.WritePhaseInfo |
---|
1891 | cif.WriteCIFitem(fp, '\n# phase info for '+str(self.name) + ' follows') |
---|
1892 | cif.WriteCIFitem(fp, '_pd_phase_name', self.name) |
---|
1893 | # TODO get esds |
---|
1894 | cellDict = self.get_cell() |
---|
1895 | defsigL = 3*[-0.00001] + 3*[-0.001] + [-0.01] # significance to use when no sigma |
---|
1896 | names = ['length_a','length_b','length_c', |
---|
1897 | 'angle_alpha','angle_beta ','angle_gamma', |
---|
1898 | 'volume'] |
---|
1899 | for key, val in cellDict.items(): |
---|
1900 | cif.WriteCIFitem(fp, '_cell_' + key, G2mth.ValEsd(val)) |
---|
1901 | |
---|
1902 | cif.WriteCIFitem(fp, '_symmetry_cell_setting', |
---|
1903 | self.data['General']['SGData']['SGSys']) |
---|
1904 | |
---|
1905 | spacegroup = self.data['General']['SGData']['SpGrp'].strip() |
---|
1906 | # regularize capitalization and remove trailing H/R |
---|
1907 | spacegroup = spacegroup[0].upper() + spacegroup[1:].lower().rstrip('rh ') |
---|
1908 | cif.WriteCIFitem(fp, '_symmetry_space_group_name_H-M', spacegroup) |
---|
1909 | |
---|
1910 | # generate symmetry operations including centering and center of symmetry |
---|
1911 | SymOpList, offsetList, symOpList, G2oprList, G2opcodes = G2spc.AllOps( |
---|
1912 | self.data['General']['SGData']) |
---|
1913 | cif.WriteCIFitem(fp, 'loop_\n _space_group_symop_id\n _space_group_symop_operation_xyz') |
---|
1914 | for i, op in enumerate(SymOpList,start=1): |
---|
1915 | cif.WriteCIFitem(fp, ' {:3d} {:}'.format(i,op.lower())) |
---|
1916 | |
---|
1917 | # TODO skipped histograms, exports/G2export_CIF.py:880 |
---|
1918 | |
---|
1919 | # report atom params |
---|
1920 | if self.data['General']['Type'] in ['nuclear','macromolecular']: #this needs macromolecular variant, etc! |
---|
1921 | cif.WriteAtomsNuclear(fp, self.data, self.name, parmDict, sigDict, []) |
---|
1922 | # self._WriteAtomsNuclear(fp, parmDict, sigDict) |
---|
1923 | else: |
---|
1924 | raise Exception("no export for "+str(self.data['General']['Type'])+" coordinates implemented") |
---|
1925 | # report cell contents |
---|
1926 | cif.WriteComposition(fp, self.data, self.name, parmDict) |
---|
1927 | if not quickmode and self.data['General']['Type'] == 'nuclear': # report distances and angles |
---|
1928 | # WriteDistances(fp,self.name,SymOpList,offsetList,symOpList,G2oprList) |
---|
1929 | raise NotImplementedError("only quickmode currently supported") |
---|
1930 | if 'Map' in self.data['General'] and 'minmax' in self.data['General']['Map']: |
---|
1931 | cif.WriteCIFitem(fp,'\n# Difference density results') |
---|
1932 | MinMax = self.data['General']['Map']['minmax'] |
---|
1933 | cif.WriteCIFitem(fp,'_refine_diff_density_max',G2mth.ValEsd(MinMax[0],-0.009)) |
---|
1934 | cif.WriteCIFitem(fp,'_refine_diff_density_min',G2mth.ValEsd(MinMax[1],-0.009)) |
---|
1935 | |
---|
1936 | |
---|
1937 | def set_refinements(self, refs): |
---|
1938 | """Sets the refinement parameter 'key' to the specification 'value' |
---|
1939 | |
---|
1940 | :param dict refs: A dictionary of the parameters to be set. See |
---|
1941 | :ref:`Phase_parameters_table` for a description of |
---|
1942 | this dictionary. |
---|
1943 | |
---|
1944 | :returns: None""" |
---|
1945 | for key, value in refs.items(): |
---|
1946 | if key == "Cell": |
---|
1947 | self.data['General']['Cell'][0] = value |
---|
1948 | |
---|
1949 | elif key == "Atoms": |
---|
1950 | for atomlabel, atomrefinement in value.items(): |
---|
1951 | if atomlabel == 'all': |
---|
1952 | for atom in self.atoms(): |
---|
1953 | atom.refinement_flags = atomrefinement |
---|
1954 | else: |
---|
1955 | atom = self.atom(atomlabel) |
---|
1956 | if atom is None: |
---|
1957 | raise ValueError("No such atom: " + atomlabel) |
---|
1958 | atom.refinement_flags = atomrefinement |
---|
1959 | |
---|
1960 | elif key == "LeBail": |
---|
1961 | hists = self.data['Histograms'] |
---|
1962 | for hname, hoptions in hists.items(): |
---|
1963 | if 'LeBail' not in hoptions: |
---|
1964 | hoptions['newLeBail'] = bool(True) |
---|
1965 | hoptions['LeBail'] = bool(value) |
---|
1966 | else: |
---|
1967 | raise ValueError("Unknown key:", key) |
---|
1968 | |
---|
1969 | def clear_refinements(self, refs): |
---|
1970 | """Clears a given set of parameters. |
---|
1971 | |
---|
1972 | :param dict refs: The parameters to clear""" |
---|
1973 | for key, value in refs.items(): |
---|
1974 | if key == "Cell": |
---|
1975 | self.data['General']['Cell'][0] = False |
---|
1976 | elif key == "Atoms": |
---|
1977 | cx, ct, cs, cia = self.data['General']['AtomPtrs'] |
---|
1978 | |
---|
1979 | for atomlabel in value: |
---|
1980 | atom = self.atom(atomlabel) |
---|
1981 | # Set refinement to none |
---|
1982 | atom.refinement_flags = ' ' |
---|
1983 | elif key == "LeBail": |
---|
1984 | hists = self.data['Histograms'] |
---|
1985 | for hname, hoptions in hists.items(): |
---|
1986 | if 'LeBail' not in hoptions: |
---|
1987 | hoptions['newLeBail'] = True |
---|
1988 | hoptions['LeBail'] = False |
---|
1989 | else: |
---|
1990 | raise ValueError("Unknown key:", key) |
---|
1991 | |
---|
1992 | def set_HAP_refinements(self, refs, histograms='all'): |
---|
1993 | """Sets the given HAP refinement parameters between this phase and |
---|
1994 | the given histograms |
---|
1995 | |
---|
1996 | :param dict refs: A dictionary of the parameters to be set. See |
---|
1997 | :ref:`HAP_parameters_table` for a description of this |
---|
1998 | dictionary. |
---|
1999 | :param histograms: Either 'all' (default) or a list of the histograms |
---|
2000 | whose HAP parameters will be set with this phase. Histogram and phase |
---|
2001 | must already be associated |
---|
2002 | |
---|
2003 | :returns: None |
---|
2004 | """ |
---|
2005 | if histograms == 'all': |
---|
2006 | histograms = self.data['Histograms'].values() |
---|
2007 | else: |
---|
2008 | histograms = [self.data['Histograms'][h.name] for h in histograms] |
---|
2009 | |
---|
2010 | for key, val in refs.items(): |
---|
2011 | for h in histograms: |
---|
2012 | if key == 'Babinet': |
---|
2013 | try: |
---|
2014 | sets = list(val) |
---|
2015 | except ValueError: |
---|
2016 | sets = ['BabA', 'BabU'] |
---|
2017 | for param in sets: |
---|
2018 | if param not in ['BabA', 'BabU']: |
---|
2019 | raise ValueError("Not sure what to do with" + param) |
---|
2020 | for hist in histograms: |
---|
2021 | hist['Babinet'][param][1] = True |
---|
2022 | elif key == 'Extinction': |
---|
2023 | for h in histograms: |
---|
2024 | h['Extinction'][1] = bool(val) |
---|
2025 | elif key == 'HStrain': |
---|
2026 | for h in histograms: |
---|
2027 | hist['HStrain'][1] = [bool(val) for p in hist['Hstrain'][0]] |
---|
2028 | elif key == 'Mustrain': |
---|
2029 | for h in histograms: |
---|
2030 | mustrain = h['Mustrain'] |
---|
2031 | newType = None |
---|
2032 | if isinstance(val, (unicode, str)): |
---|
2033 | if val in ['isotropic', 'uniaxial', 'generalized']: |
---|
2034 | newType = val |
---|
2035 | else: |
---|
2036 | raise ValueError("Not a Mustrain type: " + val) |
---|
2037 | elif isinstance(val, dict): |
---|
2038 | newType = val.get('type', None) |
---|
2039 | direction = val.get('direction', None) |
---|
2040 | |
---|
2041 | if newType: |
---|
2042 | mustrain[0] = newType |
---|
2043 | if newType == 'isotropic': |
---|
2044 | mustrain[2][0] = True |
---|
2045 | mustrain[5] = [False for p in mustrain[4]] |
---|
2046 | elif newType == 'uniaxial': |
---|
2047 | if 'refine' in val: |
---|
2048 | types = val['refine'] |
---|
2049 | if isinstance(types, (unicode, str)): |
---|
2050 | types = [types] |
---|
2051 | elif isinstance(types, bool): |
---|
2052 | mustrain[2][0] = types |
---|
2053 | mustrain[2][1] = types |
---|
2054 | types = [] |
---|
2055 | else: |
---|
2056 | raise ValueError("Not sure what to do with: " |
---|
2057 | + str(types)) |
---|
2058 | else: |
---|
2059 | types = [] |
---|
2060 | |
---|
2061 | for unitype in types: |
---|
2062 | if unitype == 'equatorial': |
---|
2063 | mustrain[2][0] = True |
---|
2064 | elif unitype == 'axial': |
---|
2065 | mustrain[2][1] = True |
---|
2066 | else: |
---|
2067 | msg = 'Invalid uniaxial mustrain type' |
---|
2068 | raise ValueError(msg + ': ' + unitype) |
---|
2069 | else: # newtype == 'generalized' |
---|
2070 | mustrain[2] = [False for p in mustrain[1]] |
---|
2071 | |
---|
2072 | if direction: |
---|
2073 | direction = [int(n) for n in direction] |
---|
2074 | if len(direction) != 3: |
---|
2075 | raise ValueError("Expected hkl, found", direction) |
---|
2076 | mustrain[3] = direction |
---|
2077 | elif key == 'Pref.Ori.': |
---|
2078 | for h in histograms: |
---|
2079 | h['Pref.Ori.'][2] = bool(val) |
---|
2080 | elif key == 'Show': |
---|
2081 | for h in histograms: |
---|
2082 | h['Show'] = bool(val) |
---|
2083 | elif key == 'Size': |
---|
2084 | # TODO |
---|
2085 | raise NotImplementedError() |
---|
2086 | elif key == 'Use': |
---|
2087 | for h in histograms: |
---|
2088 | h['Use'] = bool(val) |
---|
2089 | elif key == 'Scale': |
---|
2090 | for h in histograms: |
---|
2091 | h['Scale'][1] = bool(val) |
---|
2092 | |
---|
2093 | def clear_HAP_refinements(self, refs, histograms='all'): |
---|
2094 | """Clears the given HAP refinement parameters between this phase and |
---|
2095 | the given histograms |
---|
2096 | |
---|
2097 | :param dict refs: A dictionary of the parameters to be cleared. |
---|
2098 | :param histograms: Either 'all' (default) or a list of the histograms |
---|
2099 | whose HAP parameters will be cleared with this phase. Histogram and |
---|
2100 | phase must already be associated |
---|
2101 | |
---|
2102 | :returns: None |
---|
2103 | """ |
---|
2104 | if histograms == 'all': |
---|
2105 | histograms = self.data['Histograms'].values() |
---|
2106 | else: |
---|
2107 | histograms = [self.data['Histograms'][h.name] for h in histograms] |
---|
2108 | |
---|
2109 | for key, val in refs.items(): |
---|
2110 | for h in histograms: |
---|
2111 | if key == 'Babinet': |
---|
2112 | try: |
---|
2113 | sets = list(val) |
---|
2114 | except ValueError: |
---|
2115 | sets = ['BabA', 'BabU'] |
---|
2116 | for param in sets: |
---|
2117 | if param not in ['BabA', 'BabU']: |
---|
2118 | raise ValueError("Not sure what to do with" + param) |
---|
2119 | for hist in histograms: |
---|
2120 | hist['Babinet'][param][1] = False |
---|
2121 | elif key == 'Extinction': |
---|
2122 | for h in histograms: |
---|
2123 | h['Extinction'][1] = False |
---|
2124 | elif key == 'HStrain': |
---|
2125 | for h in histograms: |
---|
2126 | hist['HStrain'][1] = [False for p in hist['Hstrain'][0]] |
---|
2127 | elif key == 'Mustrain': |
---|
2128 | for h in histograms: |
---|
2129 | mustrain = h['Mustrain'] |
---|
2130 | mustrain[2] = [False for p in mustrain[1]] |
---|
2131 | mustrain[5] = [False for p in mustrain[4]] |
---|
2132 | elif key == 'Pref.Ori.': |
---|
2133 | for h in histograms: |
---|
2134 | h['Pref.Ori.'][2] = False |
---|
2135 | elif key == 'Show': |
---|
2136 | for h in histograms: |
---|
2137 | h['Show'] = False |
---|
2138 | elif key == 'Size': |
---|
2139 | raise NotImplementedError() |
---|
2140 | elif key == 'Use': |
---|
2141 | for h in histograms: |
---|
2142 | h['Use'] = False |
---|
2143 | elif key == 'Scale': |
---|
2144 | for h in histograms: |
---|
2145 | h['Scale'][1] = False |
---|
2146 | |
---|
2147 | |
---|
2148 | ########################## |
---|
2149 | # Command Line Interface # |
---|
2150 | ########################## |
---|
2151 | # Each of these takes an argparse.Namespace object as their argument, |
---|
2152 | # representing the parsed command-line arguments for the relevant subcommand. |
---|
2153 | # The argument specification for each is in the subcommands dictionary (see |
---|
2154 | # below) |
---|
2155 | |
---|
2156 | |
---|
2157 | def create(args): |
---|
2158 | """The create subcommand.""" |
---|
2159 | proj = G2Project(filename=args.filename) |
---|
2160 | |
---|
2161 | hist_objs = [] |
---|
2162 | for h in args.histograms: |
---|
2163 | hist_objs.append(proj.add_powder_histogram(h, args.iparams)) |
---|
2164 | |
---|
2165 | for p in args.phases: |
---|
2166 | proj.add_phase(p, histograms=hist_objs) |
---|
2167 | |
---|
2168 | proj.save() |
---|
2169 | |
---|
2170 | |
---|
2171 | def dump(args): |
---|
2172 | """The dump subcommand""" |
---|
2173 | if not args.histograms and not args.phases: |
---|
2174 | args.raw = True |
---|
2175 | if args.raw: |
---|
2176 | import IPython.lib.pretty as pretty |
---|
2177 | |
---|
2178 | for fname in args.files: |
---|
2179 | if args.raw: |
---|
2180 | proj, nameList = LoadDictFromProjFile(fname) |
---|
2181 | print("file:", fname) |
---|
2182 | print("names:", nameList) |
---|
2183 | for key, val in proj.items(): |
---|
2184 | print(key, ":") |
---|
2185 | pretty.pprint(val) |
---|
2186 | else: |
---|
2187 | proj = G2Project(fname) |
---|
2188 | if args.histograms: |
---|
2189 | hists = proj.histograms() |
---|
2190 | for h in hists: |
---|
2191 | print(fname, "hist", h.id, h.name) |
---|
2192 | if args.phases: |
---|
2193 | phase_list = proj.phases() |
---|
2194 | for p in phase_list: |
---|
2195 | print(fname, "phase", p.id, p.name) |
---|
2196 | |
---|
2197 | |
---|
2198 | def IPyBrowse(args): |
---|
2199 | """Load a .gpx file and then open a IPython shell to browse it |
---|
2200 | """ |
---|
2201 | for fname in args.files: |
---|
2202 | proj, nameList = LoadDictFromProjFile(fname) |
---|
2203 | msg = "\nfname {} loaded into proj (dict) with names in nameList".format(fname) |
---|
2204 | GSASIIpath.IPyBreak_base(msg) |
---|
2205 | break |
---|
2206 | |
---|
2207 | |
---|
2208 | def refine(args): |
---|
2209 | """The refine subcommand""" |
---|
2210 | proj = G2Project(args.gpxfile) |
---|
2211 | if args.refinements is None: |
---|
2212 | proj.refine() |
---|
2213 | else: |
---|
2214 | import json |
---|
2215 | with open(args.refinements) as refs: |
---|
2216 | refs = json.load(refs) |
---|
2217 | proj.do_refinements(refs['refinements']) |
---|
2218 | |
---|
2219 | |
---|
2220 | def seqrefine(args): |
---|
2221 | """The seqrefine subcommand""" |
---|
2222 | raise NotImplementedError("seqrefine is not yet implemented") |
---|
2223 | |
---|
2224 | |
---|
2225 | def export(args): |
---|
2226 | """The export subcommand""" |
---|
2227 | # Export phase as CIF to args.exportfile |
---|
2228 | proj = G2Project(args.gpxfile) |
---|
2229 | phase = proj.phase(args.phase) |
---|
2230 | phase.export_CIF(args.exportfile) |
---|
2231 | |
---|
2232 | |
---|
2233 | def _args_kwargs(*args, **kwargs): |
---|
2234 | return args, kwargs |
---|
2235 | |
---|
2236 | # A dictionary of the name of each subcommand, and a tuple |
---|
2237 | # of its associated function and a list of its arguments |
---|
2238 | # The arguments are passed directly to the add_argument() method |
---|
2239 | # of an argparse.ArgumentParser |
---|
2240 | subcommands = {"create": |
---|
2241 | (create, [_args_kwargs('filename', |
---|
2242 | help='the project file to create. should end in .gpx'), |
---|
2243 | |
---|
2244 | _args_kwargs('-i', '--iparams', |
---|
2245 | help='instrument parameter file'), |
---|
2246 | |
---|
2247 | _args_kwargs('-d', '--histograms', |
---|
2248 | nargs='+', |
---|
2249 | help='list of datafiles to add as histograms'), |
---|
2250 | |
---|
2251 | _args_kwargs('-p', '--phases', |
---|
2252 | nargs='+', |
---|
2253 | help='list of phases to add. phases are ' |
---|
2254 | 'automatically associated with all ' |
---|
2255 | 'histograms given.')]), |
---|
2256 | |
---|
2257 | "dump": (dump, [_args_kwargs('-d', '--histograms', |
---|
2258 | action='store_true', |
---|
2259 | help='list histograms in files, overrides --raw'), |
---|
2260 | |
---|
2261 | _args_kwargs('-p', '--phases', |
---|
2262 | action='store_true', |
---|
2263 | help='list phases in files, overrides --raw'), |
---|
2264 | |
---|
2265 | _args_kwargs('-r', '--raw', |
---|
2266 | action='store_true', help='dump raw file contents, default'), |
---|
2267 | |
---|
2268 | _args_kwargs('files', nargs='+')]), |
---|
2269 | |
---|
2270 | "refine": |
---|
2271 | (refine, [_args_kwargs('gpxfile', help='the project file to refine'), |
---|
2272 | _args_kwargs('refinements', |
---|
2273 | help='json file of refinements to apply. if not present' |
---|
2274 | ' refines file as-is', |
---|
2275 | default=None, |
---|
2276 | nargs='?')]), |
---|
2277 | |
---|
2278 | "seqrefine": (seqrefine, []), |
---|
2279 | "export": (export, [_args_kwargs('gpxfile', |
---|
2280 | help='the project file from which to export'), |
---|
2281 | _args_kwargs('phase', help='identifier of phase to export'), |
---|
2282 | _args_kwargs('exportfile', help='the .cif file to export to')]), |
---|
2283 | "browse": (IPyBrowse, [_args_kwargs('files', nargs='+', |
---|
2284 | help='list of files to browse')])} |
---|
2285 | |
---|
2286 | |
---|
2287 | def main(): |
---|
2288 | '''The command line interface for GSASIIscriptable. |
---|
2289 | |
---|
2290 | Executes one of the following subcommands: |
---|
2291 | |
---|
2292 | * :func:`create` |
---|
2293 | * :func:`dump` |
---|
2294 | * :func:`refine` |
---|
2295 | * :func:`seqrefine` |
---|
2296 | * :func:`export` |
---|
2297 | * :func:`browse` |
---|
2298 | |
---|
2299 | .. seealso:: |
---|
2300 | :func:`create` |
---|
2301 | :func:`dump` |
---|
2302 | :func:`refine` |
---|
2303 | :func:`seqrefine` |
---|
2304 | :func:`export` |
---|
2305 | :func:`browse` |
---|
2306 | ''' |
---|
2307 | parser = argparse.ArgumentParser() |
---|
2308 | subs = parser.add_subparsers() |
---|
2309 | |
---|
2310 | # Create all of the specified subparsers |
---|
2311 | for name, (func, args) in subcommands.items(): |
---|
2312 | new_parser = subs.add_parser(name) |
---|
2313 | for listargs, kwds in args: |
---|
2314 | new_parser.add_argument(*listargs, **kwds) |
---|
2315 | new_parser.set_defaults(func=func) |
---|
2316 | |
---|
2317 | # Parse and trigger subcommand |
---|
2318 | result = parser.parse_args() |
---|
2319 | result.func(result) |
---|
2320 | |
---|
2321 | if __name__ == '__main__': |
---|
2322 | main() |
---|