Context Navigation

source: trunk/GSASIIscriptable.py @ 3009

Last change on this file since 3009 was 3009, checked in by odonnell, 5 years ago
more documentation changes, removed (phase\|histogram)_names, add occupancy
File size: 70.3 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: