Context Navigation

source: trunk/GSASIIstrMain.py @ 4857

Last change on this file since 4857 was 4857, checked in by vondreele, 2 years ago
Add new button to Covariance page - Show last shift/esd bar plot Add PlotNamedFloatBarGraph? to G2plot to do it. Fix a missing default on get in ReportProblems? in G2strMain
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 57.6 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIstrMain: main structure routine
4	---------------------------------------
5
6	'''
7	########### SVN repository information ###################
8	# $Date: 2021-03-19 17:56:44 +0000 (Fri, 19 Mar 2021) $
9	# $Author: vondreele $
10	# $Revision: 4857 $
11	# $URL: trunk/GSASIIstrMain.py $
12	# $Id: GSASIIstrMain.py 4857 2021-03-19 17:56:44Z vondreele $
13	########### SVN repository information ###################
14	from __future__ import division, print_function
15	import platform
16	import sys
17	import os.path as ospath
18	import time
19	import math
20	import copy
21	if '2' in platform.python_version_tuple()[0]:
22	import cPickle
23	else:
24	import pickle as cPickle
25	import numpy as np
26	import numpy.linalg as nl
27	import scipy.optimize as so
28	import GSASIIpath
29	GSASIIpath.SetBinaryPath()
30	GSASIIpath.SetVersionNumber("$Revision: 4857 $")
31	import GSASIIlattice as G2lat
32	import GSASIIspc as G2spc
33	import GSASIImapvars as G2mv
34	import GSASIImath as G2mth
35	import GSASIIstrIO as G2stIO
36	import GSASIIstrMath as G2stMth
37	import GSASIIobj as G2obj
38	import GSASIIfiles as G2fil
39	import GSASIIElem as G2elem
40	import atmdata
41
42	sind = lambda x: np.sin(x*np.pi/180.)
43	cosd = lambda x: np.cos(x*np.pi/180.)
44	tand = lambda x: np.tan(x*np.pi/180.)
45	asind = lambda x: 180.*np.arcsin(x)/np.pi
46	acosd = lambda x: 180.*np.arccos(x)/np.pi
47	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
48
49	ateln2 = 8.0*math.log(2.0)
50	DEBUG = True
51
52	def ReportProblems(result,Rvals,varyList):
53	'''Create a message based results from the refinement
54	'''
55	#report on SVD 0's and highly correlated variables
56	msg = ''
57	# process singular variables; all vars go to console, first 10 to
58	# dialog window
59	psing = result[2].get('psing',[])
60	if len(psing):
61	if msg: msg += '\n'
62	m = 'Error: {} Parameter(s) dropped:'.format(len(psing))
63	msg += m
64	G2fil.G2Print(m, mode='warn')
65	m = ''
66	for i,val in enumerate(psing):
67	if i == 0:
68	msg += '\n{}'.format(varyList[val])
69	m = ' {}'.format(varyList[val])
70	else:
71	if len(m) > 70:
72	G2fil.G2Print(m, mode='warn')
73	m = ' '
74	else:
75	m += ', '
76	m += '{}'.format(varyList[val])
77	if i == 10:
78	msg += ', {}... see console for full list'.format(varyList[val])
79	elif i > 10:
80	pass
81	else:
82	msg += ', {}'.format(varyList[val])
83	if m: G2fil.G2Print(m, mode='warn')
84	SVD0 = result[2].get('SVD0',0)
85	if SVD0 == 1:
86	msg += 'Warning: Soft (SVD) singularity in the Hessian'
87	elif SVD0 > 0:
88	msg += 'Warning: {} soft (SVD) Hessian singularities'.format(SVD0)
89	SVDsing = result[2].get('SVDsing',[])
90	if len(SVDsing):
91	if msg: msg += '\n'
92	m = 'SVD problem(s) likely from:'
93	msg += m
94	G2fil.G2Print(m, mode='warn')
95	m = ''
96	for i,val in enumerate(SVDsing):
97	if i == 0:
98	msg += '\n{}'.format(varyList[val])
99	m = ' {}'.format(varyList[val])
100	else:
101	if len(m) > 70:
102	G2fil.G2Print(m, mode='warn')
103	m = ' '
104	else:
105	m += ', '
106	m += '{}'.format(varyList[val])
107	if i == 10:
108	msg += ', {}... see console for full list'.format(varyList[val])
109	elif i > 10:
110	pass
111	else:
112	msg += ', {}'.format(varyList[val])
113	if m: G2fil.G2Print(m, mode='warn')
114	#report on highly correlated variables
115	Hcorr = result[2].get('Hcorr',[])
116	if len(Hcorr) > 0:
117	if msg: msg += '\n'
118	m = 'Note highly correlated parameters:'
119	G2fil.G2Print(m, mode='warn')
120	msg += m
121	elif SVD0 > 0:
122	if msg: msg += '\n'
123	m = 'Check covariance matrix for parameter correlation'
124	G2fil.G2Print(m, mode='warn')
125	msg += m
126	for i,(v1,v2,corr) in enumerate(Hcorr):
127	if corr > .95:
128	stars = '**'
129	else:
130	stars = ' '
131	m = ' {} {} and {} (@{:.2f}%)'.format(
132	stars,varyList[v1],varyList[v2],100.*corr)
133	G2fil.G2Print(m, mode='warn')
134	if i == 5:
135	msg += '\n' + m
136	msg += '\n ... check console for more'
137	elif i < 5:
138	msg += '\n' + m
139	if msg:
140	if 'msg' not in Rvals: Rvals['msg'] = ''
141	Rvals['msg'] += msg
142
143	def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
144	calcControls,pawleyLookup,ifSeq,printFile,dlg,refPlotUpdate=None):
145	'''Core optimization routines, shared between SeqRefine and Refine
146
147	:returns: 5-tuple of ifOk (bool), Rvals (dict), result, covMatrix, sig
148	'''
149	#patch (added Oct 2020) convert variable names for parm limits to G2VarObj
150	import GSASIIscriptable as G2sc
151	G2sc.patchControls(Controls)
152	# end patch
153	# print 'current',varyList
154	# for item in parmDict: print item,parmDict[item] ######### show dict just before refinement
155	ifPrint = True
156	if ifSeq:
157	ifPrint = False
158	Rvals = {}
159	chisq0 = None
160	while True:
161	G2mv.Map2Dict(parmDict,varyList)
162	begin = time.time()
163	values = np.array(G2stMth.Dict2Values(parmDict, varyList))
164	if np.any(np.isnan(values)):
165	raise G2obj.G2Exception('ERROR - nan found in LS parameters - use Calculate/View LS parms to locate')
166	# test code to compute GOF and save for external repeat
167	#args = ([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)
168	#print '* before fit chi2',np.sum(G2stMth.errRefine(values,args)*2)
169	#fl = open('beforeFit.cpickle','wb')
170	#cPickle.dump(values,fl,1)
171	#cPickle.dump(args[:-1],fl,1)
172	#fl.close()
173	Ftol = Controls['min dM/M']
174	Xtol = Controls['SVDtol']
175	Factor = Controls['shift factor']
176	if 'Jacobian' in Controls['deriv type']:
177	result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
178	ftol=Ftol,col_deriv=True,factor=Factor,
179	args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
180	ncyc = int(result[2]['nfev']/2)
181	if refPlotUpdate is not None: refPlotUpdate(Histograms) # update plot after completion
182	elif 'analytic Hessian' in Controls['deriv type']:
183	Lamda = Controls.get('Marquardt',-3)
184	maxCyc = Controls['max cyc']
185	result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
186	args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
187	refPlotUpdate=refPlotUpdate)
188	ncyc = result[2]['num cyc']+1
189	Rvals['lamMax'] = result[2]['lamMax']
190	if 'Ouch#4' in result[2]:
191	Rvals['Aborted'] = True
192	if 'msg' in result[2]:
193	Rvals['msg'] = result[2]['msg']
194	Controls['Marquardt'] = -3 #reset to default
195	if 'chisq0' in result[2] and chisq0 is None:
196	chisq0 = result[2]['chisq0']
197	elif 'Hessian SVD' in Controls['deriv type']:
198	maxCyc = Controls['max cyc']
199	result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
200	args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
201	refPlotUpdate=refPlotUpdate)
202	if result[1] is None:
203	IfOK = False
204	covMatrix = []
205	sig = len(varyList)*[None,]
206	break
207	ncyc = result[2]['num cyc']+1
208	if 'chisq0' in result[2] and chisq0 is None:
209	chisq0 = result[2]['chisq0']
210	else: #'numeric'
211	result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
212	args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
213	ncyc = 1
214	if len(varyList):
215	ncyc = int(result[2]['nfev']/len(varyList))
216	if refPlotUpdate is not None: refPlotUpdate(Histograms) # update plot
217	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
218	# for item in table: print item,table[item] #useful debug - are things shifting?
219	runtime = time.time()-begin
220	Rvals['SVD0'] = result[2].get('SVD0',0)
221	Rvals['converged'] = result[2].get('Converged')
222	Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
223	Rvals['chisq'] = np.sum(result[2]['fvec']**2)
224	G2stMth.Values2Dict(parmDict, varyList, result[0])
225	G2mv.Dict2Map(parmDict,varyList)
226	Rvals['Nobs'] = Histograms['Nobs']
227	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100. #to %
228	Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
229	printFile.write(' Number of function calls: %d No. of observations: %d No. of parameters: %d User rejected: %d Sp. gp. extinct: %d\n'% \
230	(result[2]['nfev'],Histograms['Nobs'],len(varyList),Histograms['Nrej'],Histograms['Next']))
231	if ncyc:
232	printFile.write(' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles\n'%(runtime,runtime/ncyc,ncyc))
233	printFile.write(' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f\n'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']))
234	sig = len(varyList)*[None,]
235	if 'None' in str(type(result[1])) and ifSeq: #this bails out of a sequential refinement on singular matrix
236	IfOK = False
237	covMatrix = []
238	G2fil.G2Print ('Warning: ** Refinement failed - singular matrix **')
239	if 'Hessian' in Controls['deriv type']:
240	num = len(varyList)-1
241	# BHT -- I am not sure if this works correctly:
242	for i,val in enumerate(np.flipud(result[2]['psing'])):
243	if val:
244	G2fil.G2Print('Bad parameter: '+varyList[num-i],mode='warn')
245	else:
246	Ipvt = result[2]['ipvt']
247	for i,ipvt in enumerate(Ipvt):
248	if not np.sum(result[2]['fjac'],axis=1)[i]:
249	G2fil.G2Print('Bad parameter: '+varyList[ipvt-1],mode='warn')
250	break
251	IfOK = True
252	try:
253	covMatrix = result[1]Rvals['GOF']*2
254	sig = np.sqrt(np.diag(covMatrix))
255	Lastshft = result[0]-values
256	Rvals['Max shft/sig'] = np.max(np.nan_to_num(Lastshft/sig))
257	if np.any(np.isnan(sig)) or not sig.shape:
258	G2fil.G2Print ('* Least squares aborted - some invalid esds possible *',mode='error')
259	# report on refinement issues. Result in Rvals['msg']
260	ReportProblems(result,Rvals,varyList)
261	break #refinement succeeded - finish up!
262	except TypeError:
263	# if we get here, no result[1] (covar matrix) was returned or other calc error: refinement failed
264	IfOK = False
265	if not len(varyList):
266	covMatrix = []
267	break
268	if 'Hessian' in Controls['deriv type']:
269	SVD0 = result[2].get('SVD0')
270	if SVD0 == -1:
271	G2fil.G2Print ('** Refinement failed - singular matrix **',mode='error')
272	elif SVD0 == -2:
273	G2fil.G2Print ('** Refinement failed - other problem **',mode='error')
274	elif SVD0 > 0:
275	G2fil.G2Print ('** Refinement failed with {} SVD singularities **'.format(SVD0),mode='error')
276	else:
277	G2fil.G2Print ('** Refinement failed **',mode='error')
278	if result[1] is None:
279	IfOK = False
280	covMatrix = []
281	sig = len(varyList)*[None,]
282	# report on highly correlated variables
283	ReportProblems(result,Rvals,varyList)
284	# process singular variables
285	if dlg: break # refining interactively
286	else:
287	G2fil.G2Print ('** Refinement failed - singular matrix **',mode='error')
288	Ipvt = result[2]['ipvt']
289	for i,ipvt in enumerate(Ipvt):
290	if not np.sum(result[2]['fjac'],axis=1)[i]:
291	G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
292	del(varyList[ipvt-1])
293	break
294	if IfOK:
295	if CheckLeBail(Phases): # only needed for LeBail extraction
296	G2stMth.errRefine([],[Histograms,Phases,restraintDict,rigidbodyDict],
297	parmDict,[],calcControls,pawleyLookup,dlg)
298	G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
299	if chisq0 is not None:
300	Rvals['GOF0'] = np.sqrt(chisq0/(Histograms['Nobs']-len(varyList)))
301	return IfOK,Rvals,result,covMatrix,sig,Lastshft
302
303	def Refine(GPXfile,dlg=None,makeBack=True,refPlotUpdate=None):
304	'''Global refinement -- refines to minimize against all histograms.
305	This can be called in one of three ways, from :meth:`GSASIIdataGUI.GSASII.OnRefine` in an
306	interactive refinement, where dlg will be a wx.ProgressDialog, or non-interactively from
307	:meth:`GSASIIscriptable.G2Project.refine` or from :func:`main`, where dlg will be None.
308	'''
309	import GSASIImpsubs as G2mp
310	G2mp.InitMP()
311	import pytexture as ptx
312	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
313
314	printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
315	G2stIO.ShowBanner(printFile)
316	varyList = []
317	parmDict = {}
318	G2mv.InitVars()
319	Controls = G2stIO.GetControls(GPXfile)
320	G2stIO.ShowControls(Controls,printFile)
321	calcControls = {}
322	calcControls.update(Controls)
323	constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
324	restraintDict = G2stIO.GetRestraints(GPXfile)
325	Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
326	if not Phases:
327	G2fil.G2Print (' * ERROR - you have no phases to refine! *')
328	G2fil.G2Print (' * Refine aborted *')
329	return False,{'msg':'No phases'}
330	if not Histograms:
331	G2fil.G2Print (' * ERROR - you have no data to refine with! *')
332	G2fil.G2Print (' * Refine aborted *')
333	return False,{'msg':'No data'}
334	rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
335	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
336	rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
337	(Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
338	maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
339	calcControls['atomIndx'] = atomIndx
340	calcControls['Natoms'] = Natoms
341	calcControls['FFtables'] = FFtables
342	calcControls['BLtables'] = BLtables
343	calcControls['MFtables'] = MFtables
344	calcControls['maxSSwave'] = maxSSwave
345	hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
346	TwConstr,TwFixed = G2stIO.makeTwinFrConstr(Phases,Histograms,hapVary)
347	constrDict += TwConstr
348	fixedList += TwFixed
349	calcControls.update(controlDict)
350	histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
351	calcControls.update(controlDict)
352	varyList = rbVary+phaseVary+hapVary+histVary
353	parmDict.update(rbDict)
354	parmDict.update(phaseDict)
355	parmDict.update(hapDict)
356	parmDict.update(histDict)
357	G2stIO.GetFprime(calcControls,Histograms)
358	# do constraint processing
359	varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
360	msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
361	if msg:
362	return False,{'msg':'Unable to interpret multiplier(s): '+msg}
363	try:
364	G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict)
365	#print(G2mv.VarRemapShow(varyList))
366	#print('DependentVars',G2mv.GetDependentVars())
367	#print('IndependentVars',G2mv.GetIndependentVars())
368	except G2mv.ConstraintException:
369	G2fil.G2Print (' * ERROR - your constraints are internally inconsistent *')
370	#errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
371	#print 'Errors',errmsg
372	#if warnmsg: print 'Warnings',warnmsg
373	return False,{'msg':' Constraint error'}
374	# print G2mv.VarRemapShow(varyList)
375
376	# remove frozen vars from refinement
377	if 'parmFrozen' not in Controls:
378	Controls['parmFrozen'] = {}
379	if 'FrozenList' not in Controls['parmFrozen']:
380	Controls['parmFrozen']['FrozenList'] = []
381	parmFrozenList = Controls['parmFrozen']['FrozenList']
382	frozenList = [i for i in varyList if i in parmFrozenList]
383	if len(frozenList) != 0:
384	varyList = [i for i in varyList if i not in parmFrozenList]
385	G2fil.G2Print(
386	'Frozen refined variables (due to exceeding limits)\n\t:{}'
387	.format(frozenList))
388
389	ifSeq = False
390	printFile.write('\n Refinement results:\n')
391	printFile.write(135*'-'+'\n')
392	Rvals = {}
393	try:
394	covData = {}
395	IfOK,Rvals,result,covMatrix,sig,Lastshft = RefineCore(Controls,Histograms,Phases,restraintDict,
396	rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
397	refPlotUpdate=refPlotUpdate)
398	if IfOK:
399	sigDict = dict(zip(varyList,sig))
400	newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
401	newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
402	covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
403	'varyListStart':varyListStart,'Lastshft':Lastshft,
404	'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
405	'newCellDict':newCellDict,'freshCOV':True}
406	# add the uncertainties into the esd dictionary (sigDict)
407	sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
408	# check for variables outside their allowed range, reset and freeze them
409	frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
410	G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
411	G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
412	G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
413	G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
414	G2stIO.PrintISOmodes(printFile,Phases,parmDict,sigDict)
415	G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,pFile=printFile)
416	G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls,pFile=printFile)
417	if len(frozen):
418	if 'msg' in Rvals:
419	Rvals['msg'] += '\n'
420	else:
421	Rvals['msg'] = ''
422	msg = ('Warning: {} variable(s) refined outside limits and were frozen ({} total frozen)'
423	.format(len(frozen),len(parmFrozenList))
424	)
425	G2fil.G2Print(msg)
426	Rvals['msg'] += msg
427	elif len(parmFrozenList):
428	if 'msg' in Rvals:
429	Rvals['msg'] += '\n'
430	else:
431	Rvals['msg'] = ''
432	msg = ('Note: a total of {} variable(s) are frozen due to refining outside limits'
433	.format(len(parmFrozenList))
434	)
435	G2fil.G2Print('Note: ',msg)
436	Rvals['msg'] += msg
437	G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,parmFrozenList,makeBack)
438	printFile.close()
439	G2fil.G2Print (' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
440	G2fil.G2Print (' *** Refinement successful ***')
441	else:
442	G2fil.G2Print ('****ERROR - Refinement failed',mode='error')
443	if 'msg' in Rvals:
444	G2fil.G2Print ('Note refinement problem:',mode='warn')
445	G2fil.G2Print (Rvals['msg'],mode='warn')
446	raise G2obj.G2Exception('** ERROR: Refinement failed **')
447	except G2obj.G2RefineCancel as Msg:
448	printFile.close()
449	G2fil.G2Print (' *** Refinement stopped ***')
450	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
451	if 'msg' in Rvals:
452	Rvals['msg'] += '\n'
453	Rvals['msg'] += Msg.msg
454	if not dlg:
455	G2fil.G2Print ('Note refinement problem:',mode='warn')
456	G2fil.G2Print (Rvals['msg'],mode='warn')
457	else:
458	Rvals['msg'] = Msg.msg
459	return False,Rvals
460	# except G2obj.G2Exception as Msg: # cell metric error, others?
461	except Exception as Msg: # cell metric error, others?
462	if GSASIIpath.GetConfigValue('debug'):
463	import traceback
464	print(traceback.format_exc())
465	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
466	printFile.close()
467	G2fil.G2Print (' *** Refinement error ***')
468	if 'msg' in Rvals:
469	Rvals['msg'] += '\n\n'
470	Rvals['msg'] += Msg.msg
471	if not dlg:
472	G2fil.G2Print ('Note refinement problem:',mode='warn')
473	G2fil.G2Print (Rvals['msg'],mode='warn')
474	else:
475	Rvals['msg'] = Msg.msg
476	return False,Rvals
477
478	#for testing purposes, create a file for testderiv
479	if GSASIIpath.GetConfigValue('debug'): # and IfOK:
480	#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
481	fl = open(ospath.splitext(GPXfile)[0]+'.testDeriv','wb')
482	cPickle.dump(result[0],fl,1)
483	cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
484	cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
485	cPickle.dump(parmDict,fl,1)
486	cPickle.dump(varyList,fl,1)
487	cPickle.dump(calcControls,fl,1)
488	cPickle.dump(pawleyLookup,fl,1)
489	fl.close()
490	if dlg:
491	return True,Rvals
492	elif 'msg' in Rvals:
493	G2fil.G2Print ('Reported from refinement:',mode='warn')
494	G2fil.G2Print (Rvals['msg'],mode='warn')
495
496	def CheckLeBail(Phases):
497	'''Check if there is a LeBail extraction in any histogram
498
499	:returns: True if there is at least one LeBail flag turned on, False otherwise
500	'''
501	for key in Phases:
502	phase = Phases[key]
503	for h in phase['Histograms']:
504	#phase['Histograms'][h]
505	if not phase['Histograms'][h]['Use']: continue
506	try:
507	if phase['Histograms'][h]['LeBail']:
508	return True
509	except KeyError: #HKLF & old gpx files
510	pass
511	return False
512
513	def DoLeBail(GPXfile,dlg=None,cycles=3,refPlotUpdate=None):
514	'''Fit LeBail intensities without changes to any other refined parameters.
515	This is a stripped-down version of :func:`Refine` that does not perform
516	any refinement cycles
517	'''
518	import GSASIImpsubs as G2mp
519	G2mp.InitMP()
520	import pytexture as ptx
521	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
522
523	#varyList = []
524	parmDict = {}
525	Controls = G2stIO.GetControls(GPXfile)
526	calcControls = {}
527	calcControls.update(Controls)
528	constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
529	restraintDict = {}
530	Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
531	if not Phases:
532	G2fil.G2Print (' * ERROR - you have no phases to refine! *')
533	return False,{'msg':'No phases'}
534	if not Histograms:
535	G2fil.G2Print (' * ERROR - you have no data to refine with! *')
536	return False,{'msg':'No data'}
537	if not CheckLeBail(Phases):
538	msg = 'Warning: There are no histograms with LeBail extraction enabled'
539	G2fil.G2Print ('* '+msg+' *')
540	return False,{'msg': msg}
541	rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
542	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
543	rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
544	(Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
545	maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,Print=False)
546	calcControls['atomIndx'] = atomIndx
547	calcControls['Natoms'] = Natoms
548	calcControls['FFtables'] = FFtables
549	calcControls['BLtables'] = BLtables
550	calcControls['MFtables'] = MFtables
551	calcControls['maxSSwave'] = maxSSwave
552	hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,Print=False)
553	calcControls.update(controlDict)
554	histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,Print=False)
555	calcControls.update(controlDict)
556	parmDict.update(rbDict)
557	parmDict.update(phaseDict)
558	parmDict.update(hapDict)
559	parmDict.update(histDict)
560	G2stIO.GetFprime(calcControls,Histograms)
561	try:
562	for i in range(cycles):
563	M = G2stMth.errRefine([],[Histograms,Phases,restraintDict,rigidbodyDict],parmDict,[],calcControls,pawleyLookup,dlg)
564	G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
565	if refPlotUpdate is not None: refPlotUpdate(Histograms,i)
566	Rvals = {}
567	Rvals['chisq'] = np.sum(M**2)
568	Rvals['Nobs'] = Histograms['Nobs']
569	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100. #to %
570	Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs'])) # no variables
571
572	covData = {'variables':0,'varyList':[],'sig':[],'Rvals':Rvals,
573	'varyListStart':[],
574	'covMatrix':np.zeros([0,0]),'title':GPXfile,
575	#'newAtomDict':newAtomDict,'newCellDict':newCellDict,
576	'freshCOV':True}
577
578	G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,[],True)
579	G2fil.G2Print (' *** LeBail fit completed ***')
580	return True,Rvals
581	except Exception as Msg:
582	G2fil.G2Print (' *** LeBail fit error ***')
583	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
584	if GSASIIpath.GetConfigValue('debug'):
585	import traceback
586	print(traceback.format_exc())
587	return False,{'msg':Msg.msg}
588
589	def phaseCheck(phaseVary,Phases,histogram):
590	'''
591	Removes unused parameters from phase varylist if phase not in histogram
592	for seq refinement removes vars in "Fix FXU" and "FixedSeqVars" here
593	'''
594	NewVary = []
595	for phase in Phases:
596	if histogram not in Phases[phase]['Histograms']: continue
597	if Phases[phase]['Histograms'][histogram]['Use']:
598	pId = Phases[phase]['pId']
599	newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
600	FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
601	if 'F' in FixVals:
602	newVary = [item for item in newVary if not 'Afrac' in item]
603	if 'X' in FixVals:
604	newVary = [item for item in newVary if not 'dA' in item]
605	if 'U' in FixVals:
606	newVary = [item for item in newVary if not 'AU' in item]
607	if 'M' in FixVals:
608	newVary = [item for item in newVary if not 'AM' in item]
609	removeVars = Phases[phase]['Histograms'][histogram].get('FixedSeqVars',[])
610	newVary = [item for item in newVary if item not in removeVars]
611	NewVary += newVary
612	return NewVary
613
614	def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
615	'''Perform a sequential refinement -- cycles through all selected histgrams,
616	one at a time
617	'''
618	import GSASIImpsubs as G2mp
619	G2mp.InitMP()
620	import pytexture as ptx
621	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
622	msgs = {}
623	printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
624	G2fil.G2Print ('Starting Sequential Refinement')
625	G2stIO.ShowBanner(printFile)
626	Controls = G2stIO.GetControls(GPXfile)
627	preFrozenCount = 0
628	for h in Controls['parmFrozen']:
629	if h == 'FrozenList': continue
630	preFrozenCount += len(Controls['parmFrozen'][h])
631	G2stIO.ShowControls(Controls,printFile,SeqRef=True,
632	preFrozenCount=preFrozenCount)
633	restraintDict = G2stIO.GetRestraints(GPXfile)
634	Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
635	if not Phases:
636	G2fil.G2Print (' * ERROR - you have no phases to refine! *')
637	G2fil.G2Print (' * Refine aborted *')
638	return False,'No phases'
639	if not Histograms:
640	G2fil.G2Print (' * ERROR - you have no data to refine with! *')
641	G2fil.G2Print (' * Refine aborted *')
642	return False,'No data'
643	rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
644	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
645	rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
646	G2mv.InitVars()
647	(Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
648	maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,
649	Print=False,pFile=printFile,seqRef=True)
650	for item in phaseVary:
651	if '::A0' in item:
652	G2fil.G2Print ('** WARNING - lattice parameters should not be refined in a sequential refinement **')
653	G2fil.G2Print ('** instead use the Dij parameters for each powder histogram **')
654	return False,'Lattice parameter refinement error - see console message'
655	if '::C(' in item:
656	G2fil.G2Print ('** WARNING - phase texture parameters should not be refined in a sequential refinement **')
657	G2fil.G2Print ('** instead use the C(L,N) parameters for each powder histogram **')
658	return False,'Phase texture refinement error - see console message'
659	if 'Seq Data' in Controls:
660	histNames = Controls['Seq Data']
661	else: # patch from before Controls['Seq Data'] was implemented?
662	histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
663	if Controls.get('Reverse Seq'):
664	histNames.reverse()
665	SeqResult = G2stIO.GetSeqResult(GPXfile)
666	# SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
667	Histo = {}
668	NewparmDict = {}
669	G2stIO.SetupSeqSavePhases(GPXfile)
670	msgs['steepestNum'] = 0
671	msgs['maxshift/sigma'] = []
672	for ihst,histogram in enumerate(histNames):
673	if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
674	G2fil.G2Print('\nRefining with '+str(histogram))
675	G2mv.InitVars()
676	(Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
677	FFtables,BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
678	Phases,restraintDict,rbIds,
679	Print=False,pFile=printFile,seqRef=True)
680	ifPrint = False
681	if dlg:
682	dlg.SetTitle('Residual for histogram '+str(ihst))
683	calcControls = {}
684	calcControls['atomIndx'] = atomIndx
685	calcControls['Natoms'] = Natoms
686	calcControls['FFtables'] = FFtables
687	calcControls['BLtables'] = BLtables
688	calcControls['MFtables'] = MFtables
689	calcControls['maxSSwave'] = maxSSwave
690	if histogram not in Histograms:
691	G2fil.G2Print("Error: not found!")
692	raise G2obj.G2Exception("refining with invalid histogram {}".
693	format(histogram))
694	hId = Histograms[histogram]['hId']
695	redphaseVary = phaseCheck(phaseVary,Phases,histogram)
696	Histo = {histogram:Histograms[histogram],}
697	hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
698	calcControls.update(controlDict)
699	histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
700	calcControls.update(controlDict)
701	varyList = rbVary+redphaseVary+hapVary+histVary
702	# if not ihst:
703	# save the initial vary list, but without histogram numbers on parameters
704	saveVaryList = varyList[:]
705	for i,item in enumerate(saveVaryList):
706	items = item.split(':')
707	if items[1]:
708	items[1] = ''
709	item = ':'.join(items)
710	saveVaryList[i] = item
711	if not ihst:
712	SeqResult['varyList'] = saveVaryList
713	else:
714	SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
715	parmDict = {}
716	parmDict.update(rbDict)
717	parmDict.update(phaseDict)
718	parmDict.update(hapDict)
719	parmDict.update(histDict)
720	if Controls['Copy2Next']: # update with parms from last histogram
721	#parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
722	for parm in NewparmDict:
723	if parm in parmDict:
724	parmDict[parm] = NewparmDict[parm]
725	elif histogram in SeqResult: # update phase from last seq ref
726	NewparmDict = SeqResult[histogram].get('parmDict',{})
727	for parm in NewparmDict:
728	if '::' in parm and parm in parmDict:
729	parmDict[parm] = NewparmDict[parm]
730
731	G2stIO.GetFprime(calcControls,Histo)
732	# do constraint processing
733	#reload(G2mv) # debug
734	constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
735	varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
736	msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
737	if msg:
738	return False,'Unable to interpret multiplier(s): '+msg
739	try:
740	groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict,SeqHist=hId)
741	# if GSASIIpath.GetConfigValue('debug'): print("DBG_"+
742	# G2mv.VarRemapShow(varyList,True))
743	# print('DependentVars',G2mv.GetDependentVars())
744	# print('IndependentVars',G2mv.GetIndependentVars())
745	constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
746	except G2mv.ConstraintException:
747	G2fil.G2Print (' * ERROR - your constraints are internally inconsistent *')
748	#errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
749	#print 'Errors',errmsg
750	#if warnmsg: print 'Warnings',warnmsg
751	return False,' Constraint error'
752	#print G2mv.VarRemapShow(varyList)
753	if not ihst:
754	# first histogram to refine against
755	firstVaryList = []
756	for item in varyList:
757	items = item.split(':')
758	if items[1]:
759	items[1] = ''
760	item = ':'.join(items)
761	firstVaryList.append(item)
762	newVaryList = firstVaryList
763	else:
764	newVaryList = []
765	for item in varyList:
766	items = item.split(':')
767	if items[1]:
768	items[1] = ''
769	item = ':'.join(items)
770	newVaryList.append(item)
771	if newVaryList != firstVaryList and Controls['Copy2Next']:
772	# variable lists are expected to match between sequential refinements when Copy2Next is on
773	#print '**** ERROR - variable list for this histogram does not match previous'
774	#print ' Copy of variables is not possible'
775	#print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
776	combined = list(set(firstVaryList+newVaryList))
777	c = [var for var in combined if var not in newVaryList]
778	p = [var for var in combined if var not in firstVaryList]
779	G2fil.G2Print('* Variables change *')
780	for typ,vars in [('Removed',c),('Added',p)]:
781	line = ' '+typ+': '
782	if vars:
783	for var in vars:
784	if len(line) > 70:
785	G2fil.G2Print(line)
786	line = ' '
787	line += var + ', '
788	else:
789	line += 'none, '
790	G2fil.G2Print(line[:-2])
791	firstVaryList = newVaryList
792
793	ifSeq = True
794	printFile.write('\n Refinement results for histogram id {}: {}\n'
795	.format(hId,histogram))
796	printFile.write(135*'-'+'\n')
797	# remove frozen vars
798	if 'parmFrozen' not in Controls:
799	Controls['parmFrozen'] = {}
800	if histogram not in Controls['parmFrozen']:
801	Controls['parmFrozen'][histogram] = []
802	parmFrozenList = Controls['parmFrozen'][histogram]
803	frozenList = [i for i in varyList if i in parmFrozenList]
804	if len(frozenList) != 0:
805	varyList = [i for i in varyList if i not in parmFrozenList]
806	s = ''
807	for a in frozenList:
808	if s:
809	s+= ', '
810	s += a
811	printFile.write(
812	' The following refined variables have previously been frozen due to exceeding limits:\n\t{}\n'
813	.format(s))
814	try:
815	IfOK,Rvals,result,covMatrix,sig,Lastshft = RefineCore(Controls,Histo,Phases,restraintDict,
816	rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
817	refPlotUpdate=refPlotUpdate)
818	try:
819	shft = '%.4f'% Rvals['Max shft/sig']
820	except:
821	shft = '?'
822	G2fil.G2Print (' wR = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f, last delta chi = %.4f, last shft/sig = %s'%(
823	Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'],shft))
824	if Rvals.get('lamMax',0) >= 10.:
825	msgs['steepestNum'] += 1
826	if Rvals.get('Max shft/sig'):
827	msgs['maxshift/sigma'].append(Rvals['Max shft/sig'])
828	# add the uncertainties into the esd dictionary (sigDict)
829	if not IfOK:
830	G2fil.G2Print('***** Sequential refinement failed at histogram '+histogram,mode='warn')
831	break
832	sigDict = dict(zip(varyList,sig))
833	# the uncertainties for dependent constrained parms into the esd dict
834	sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
835	# check for variables outside their allowed range, reset and freeze them
836	frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
837	msg = None
838	if len(frozen) > 0:
839	msg = ('Hist {}: {} variables were outside limits and were frozen (now {} frozen total)'
840	.format(ihst,len(frozen),len(parmFrozenList)))
841	G2fil.G2Print(msg)
842	msg = (' {} variables were outside limits and were frozen (now {} frozen total)'
843	.format(len(frozen),len(parmFrozenList)))
844	for p in frozen:
845	if p not in varyList:
846	print('Frozen Warning: {} not in varyList. This should not happen!'.format(p))
847	continue
848	i = varyList.index(p)
849	result[0][i] = parmDict[p]
850	sig[i] = -0.1
851	# a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
852	depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
853	newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
854	newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
855	histRefData = {
856	'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
857	'varyListStart':varyListStart,'Lastshft':Lastshft,
858	'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
859	'newCellDict':newCellDict,'depParmDict':depParmDict,
860	'constraintInfo':constraintInfo,
861	'parmDict':parmDict,
862	}
863	SeqResult[histogram] = histRefData
864	G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
865	# G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
866	G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
867	G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile,seq=True)
868	G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,Controls['parmFrozen'])
869	if msg:
870	printFile.write(msg+'\n')
871	NewparmDict = {}
872	# make dict of varied parameters in current histogram, renamed to
873	# next histogram, for use in next refinement.
874	if Controls['Copy2Next'] and ihst < len(histNames)-1:
875	hId = Histo[histogram]['hId'] # current histogram
876	nexthId = Histograms[histNames[ihst+1]]['hId']
877	for parm in set(list(varyList)+list(varyListStart)):
878	items = parm.split(':')
879	if len(items) < 3:
880	continue
881	if str(hId) in items[1]:
882	items[1] = str(nexthId)
883	newparm = ':'.join(items)
884	NewparmDict[newparm] = parmDict[parm]
885	else:
886	if items[2].startswith('dA'): parm = parm.replace(':dA',':A')
887	NewparmDict[parm] = parmDict[parm]
888
889	except G2obj.G2RefineCancel as Msg:
890	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
891	printFile.close()
892	G2fil.G2Print (' *** Refinement stopped ***')
893	return False,Msg.msg
894	except G2obj.G2Exception as Msg: # cell metric error, others?
895	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
896	printFile.close()
897	G2fil.G2Print (' *** Refinement error ***')
898	return False,Msg.msg
899	if GSASIIpath.GetConfigValue('Show_timing'):
900	t2 = time.time()
901	G2fil.G2Print("Fit step time {:.2f} sec.".format(t2-t1))
902	t1 = t2
903	SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
904	try:
905	G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
906	except Exception as msg:
907	print('Error reading Sequential results\n',str(msg))
908	if GSASIIpath.GetConfigValue('debug'):
909	import traceback
910	print(traceback.format_exc())
911	postFrozenCount = 0
912	for h in Controls['parmFrozen']:
913	if h == 'FrozenList': continue
914	postFrozenCount += len(Controls['parmFrozen'][h])
915	if postFrozenCount:
916	msgs['Frozen'] = 'Ending refinement with {} Frozen variables ({} added now)\n'.format(postFrozenCount,postFrozenCount-preFrozenCount)
917	printFile.write('\n'+msgs['Frozen'])
918	printFile.close()
919	G2fil.G2Print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
920	G2fil.G2Print (' *** Sequential refinement successful ***')
921	return True,msgs
922
923	def dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList):
924	'''Find variables in the parameters dict that are outside the ranges
925	(in parmMinDict and parmMaxDict) and set them to the limits values.
926	Add any such variables into the list of frozen variable
927	(parmFrozenList). Returns a list of newly frozen variables, if any.
928	'''
929	parmMinDict = Controls.get('parmMinDict',{})
930	parmMaxDict = Controls.get('parmMaxDict',{})
931	freeze = []
932	if parmMinDict or parmMaxDict:
933	for name in varyList:
934	if name not in parmDict: continue
935	n,val = G2obj.prmLookup(name,parmMinDict)
936	if n is not None:
937	if parmDict[name] < parmMinDict[n]:
938	parmDict[name] = parmMinDict[n]
939	sigDict[name] = 0.0
940	freeze.append(name)
941	continue
942	n,val = G2obj.prmLookup(name,parmMaxDict)
943	if n is not None:
944	if parmDict[name] > parmMaxDict[n]:
945	parmDict[name] = parmMaxDict[n]
946	sigDict[name] = 0.0
947	freeze.append(name)
948	continue
949	for v in freeze:
950	if v not in parmFrozenList:
951	parmFrozenList.append(v)
952	return freeze
953
954	def RetDistAngle(DisAglCtls,DisAglData,dlg=None):
955	'''Compute and return distances and angles
956
957	:param dict DisAglCtls: contains distance/angle radii usually defined using
958	:func:`GSASIIctrlGUI.DisAglDialog`
959	:param dict DisAglData: contains phase data:
960	Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
961	for distance/angle origins and atoms to be used as targets.
962	Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
963
964	:returns: AtomLabels,DistArray,AngArray where:
965
966	AtomLabels is a dict of atom labels, keys are the atom number
967
968	DistArray is a dict keyed by the origin atom number where the value is a list
969	of distance entries. The value for each distance is a list containing:
970
971	0) the target atom number (int);
972	1) the unit cell offsets added to x,y & z (tuple of int values)
973	2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
974	3) an interatomic distance in A (float)
975	4) an uncertainty on the distance in A or 0.0 (float)
976
977	AngArray is a dict keyed by the origin (central) atom number where
978	the value is a list of
979	angle entries. The value for each angle entry consists of three values:
980
981	0) a distance item reference for one neighbor (int)
982	1) a distance item reference for a second neighbor (int)
983	2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
984
985	The AngArray distance reference items refer directly to the index of the items in the
986	DistArray item for the list of distances for the central atom.
987	'''
988	import numpy.ma as ma
989
990	SGData = DisAglData['SGData']
991	Cell = DisAglData['Cell']
992
993	Amat,Bmat = G2lat.cell2AB(Cell[:6])
994	covData = {}
995	if len(DisAglData.get('covData',{})):
996	covData = DisAglData['covData']
997	covMatrix = covData['covMatrix']
998	varyList = covData['varyList']
999	pfx = str(DisAglData['pId'])+'::'
1000
1001	Factor = DisAglCtls['Factors']
1002	Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
1003	indices = (-2,-1,0,1,2)
1004	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1005	origAtoms = DisAglData['OrigAtoms']
1006	targAtoms = DisAglData['TargAtoms']
1007	AtomLabels = {}
1008	for Oatom in origAtoms:
1009	AtomLabels[Oatom[0]] = Oatom[1]
1010	for Oatom in targAtoms:
1011	AtomLabels[Oatom[0]] = Oatom[1]
1012	DistArray = {}
1013	AngArray = {}
1014	for iO,Oatom in enumerate(origAtoms):
1015	DistArray[Oatom[0]] = []
1016	AngArray[Oatom[0]] = []
1017	OxyzNames = ''
1018	IndBlist = []
1019	Dist = []
1020	Vect = []
1021	VectA = []
1022	angles = []
1023	for Tatom in targAtoms:
1024	Xvcov = []
1025	TxyzNames = ''
1026	if len(DisAglData.get('covData',{})):
1027	OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
1028	TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
1029	Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
1030	result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
1031	BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
1032	AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
1033	for [Txyz,Top,Tunit,Spn] in result:
1034	Dx = (Txyz-np.array(Oatom[3:6]))+Units
1035	dx = np.inner(Amat,Dx)
1036	dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
1037	IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
1038	if np.any(IndB):
1039	for indb in IndB:
1040	for i in range(len(indb)):
1041	if str(dx.T[indb][i]) not in IndBlist:
1042	IndBlist.append(str(dx.T[indb][i]))
1043	unit = Units[indb][i]
1044	tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
1045	sig = 0.0
1046	if len(Xvcov):
1047	pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
1048	sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
1049	Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
1050	if (Dist[-1][-2]-AsumR) <= 0.:
1051	Vect.append(dx.T[indb][i]/Dist[-1][-2])
1052	VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
1053	else:
1054	Vect.append([0.,0.,0.])
1055	VectA.append([])
1056	if dlg is not None:
1057	dlg.Update(iO,newmsg='Atoms done=%d'%(iO))
1058	for D in Dist:
1059	DistArray[Oatom[0]].append(D[1:])
1060	Vect = np.array(Vect)
1061	angles = np.zeros((len(Vect),len(Vect)))
1062	angsig = np.zeros((len(Vect),len(Vect)))
1063	for i,veca in enumerate(Vect):
1064	if np.any(veca):
1065	for j,vecb in enumerate(Vect):
1066	if np.any(vecb):
1067	angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
1068	if i <= j: continue
1069	AngArray[Oatom[0]].append((i,j,
1070	G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
1071	return AtomLabels,DistArray,AngArray
1072
1073	def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
1074	'''Print distances and angles
1075
1076	:param dict DisAglCtls: contains distance/angle radii usually defined using
1077	:func:`GSASIIctrlGUI.DisAglDialog`
1078	:param dict DisAglData: contains phase data:
1079	Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
1080	for distance/angle origins and atoms to be used as targets.
1081	Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
1082	:param file out: file object for output. Defaults to sys.stdout.
1083	'''
1084	def MyPrint(s):
1085	out.write(s+'\n')
1086	# print(s,file=out) # use in Python 3
1087
1088	def ShowBanner(name):
1089	MyPrint(80'')
1090	MyPrint(' Interatomic Distances and Angles for phase '+name)
1091	MyPrint((80'')+'\n')
1092
1093	ShowBanner(DisAglCtls['Name'])
1094	SGData = DisAglData['SGData']
1095	SGtext,SGtable = G2spc.SGPrint(SGData)
1096	for line in SGtext: MyPrint(line)
1097	if len(SGtable) > 1:
1098	for i,item in enumerate(SGtable[::2]):
1099	if 2*i+1 == len(SGtable):
1100	line = ' %s'%(item.ljust(30))
1101	else:
1102	line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
1103	MyPrint(line)
1104	else:
1105	MyPrint(' ( 1) %s'%(SGtable[0])) #triclinic case
1106	Cell = DisAglData['Cell']
1107
1108	Amat,Bmat = G2lat.cell2AB(Cell[:6])
1109	covData = {}
1110	if len(DisAglData.get('covData',{})):
1111	covData = DisAglData['covData']
1112	pfx = str(DisAglData['pId'])+'::'
1113	A = G2lat.cell2A(Cell[:6])
1114	cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
1115	names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
1116	valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
1117	line = '\n Unit cell:'
1118	for name,vals in zip(names,valEsd):
1119	line += name+vals
1120	MyPrint(line)
1121	else:
1122	MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
1123	' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
1124	('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
1125
1126	AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
1127	origAtoms = DisAglData['OrigAtoms']
1128	for Oatom in origAtoms:
1129	i = Oatom[0]
1130	Dist = DistArray[i]
1131	nDist = len(Dist)
1132	angles = np.zeros((nDist,nDist))
1133	angsig = np.zeros((nDist,nDist))
1134	for k,j,tup in AngArray[i]:
1135	angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
1136	line = ''
1137	for i,x in enumerate(Oatom[3:6]):
1138	line += ('%12.5f'%x).rstrip('0')
1139	MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
1140	MyPrint(80'')
1141	line = ''
1142	for dist in Dist[:-1]:
1143	line += '%12s'%(AtomLabels[dist[0]].center(12))
1144	MyPrint(' To cell +(sym. op.) dist. '+line.rstrip())
1145	BVS = {}
1146	BVdat = {}
1147	Otyp = G2elem.FixValence(Oatom[2]).split('+')[0].split('-')[0]
1148	BVox = [BV for BV in atmdata.BVSoxid[Otyp] if '+' in BV]
1149	if len(BVox):
1150	BVS = {BV:0.0 for BV in BVox}
1151	BVdat = {BV:dict(zip(['O','F','Cl'],atmdata.BVScoeff[BV])) for BV in BVox}
1152	pvline = 'Bond Valence sums for: '
1153	for i,dist in enumerate(Dist):
1154	line = ''
1155	for j,angle in enumerate(angles[i][0:i]):
1156	sig = angsig[i][j]
1157	if angle:
1158	if sig:
1159	line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
1160	else:
1161	val = '%.3f'%(angle)
1162	line += '%12s'%(val.center(12))
1163	else:
1164	line += 12*' '
1165	if dist[4]: #sig exists!
1166	val = G2mth.ValEsd(dist[3],dist[4])
1167	else:
1168	val = '%8.4f'%(dist[3])
1169	if len(BVox):
1170	Tatm = G2elem.FixValence(DisAglData['TargAtoms'][dist[0]][2]).split('-')[0]
1171	if Tatm in ['O','F','Cl']:
1172	for BV in BVox:
1173	BVS[BV] += np.exp((BVdat[BV][Tatm]-dist[3])/0.37)
1174	tunit = '[%2d%2d%2d]'% dist[1]
1175	MyPrint((' %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
1176	if len(BVox):
1177	MyPrint(80'')
1178	for BV in BVox:
1179	pvline += ' %s: %.2f '%(BV,BVS[BV])
1180	MyPrint(pvline)
1181
1182	def DisAglTor(DATData):
1183	'Needs a doc string'
1184	SGData = DATData['SGData']
1185	Cell = DATData['Cell']
1186
1187	Amat,Bmat = G2lat.cell2AB(Cell[:6])
1188	covData = {}
1189	pfx = ''
1190	if 'covData' in DATData:
1191	covData = DATData['covData']
1192	pfx = str(DATData['pId'])+'::'
1193	Datoms = []
1194	Oatoms = []
1195	for i,atom in enumerate(DATData['Datoms']):
1196	symop = atom[-1].split('+')
1197	if len(symop) == 1:
1198	symop.append('0,0,0')
1199	symop[0] = int(symop[0])
1200	symop[1] = eval(symop[1])
1201	atom.append(symop)
1202	Datoms.append(atom)
1203	oatom = DATData['Oatoms'][i]
1204	names = ['','','']
1205	if pfx:
1206	names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
1207	oatom += [names,]
1208	Oatoms.append(oatom)
1209	atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
1210	if DATData['Natoms'] == 4: #torsion
1211	Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1212	G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
1213	G2fil.G2Print (' NB: Atom sequence determined by selection order')
1214	return # done with torsion
1215	elif DATData['Natoms'] == 3: #angle
1216	Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1217	G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
1218	G2fil.G2Print (' NB: Atom sequence determined by selection order')
1219	else: #2 atoms - distance
1220	Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1221	G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
1222
1223	def BestPlane(PlaneData):
1224	'Needs a doc string'
1225
1226	def ShowBanner(name):
1227	G2fil.G2Print (80'')
1228	G2fil.G2Print (' Best plane result for phase '+name)
1229	G2fil.G2Print (80'','\n')
1230
1231	ShowBanner(PlaneData['Name'])
1232
1233	Cell = PlaneData['Cell']
1234	Amat,Bmat = G2lat.cell2AB(Cell[:6])
1235	Atoms = PlaneData['Atoms']
1236	sumXYZ = np.zeros(3)
1237	XYZ = []
1238	Natoms = len(Atoms)
1239	for atom in Atoms:
1240	xyz = np.array(atom[3:6])
1241	XYZ.append(xyz)
1242	sumXYZ += xyz
1243	sumXYZ /= Natoms
1244	XYZ = np.array(XYZ)-sumXYZ
1245	XYZ = np.inner(Amat,XYZ).T
1246	Zmat = np.zeros((3,3))
1247	for i,xyz in enumerate(XYZ):
1248	Zmat += np.outer(xyz.T,xyz)
1249	G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
1250	Evec,Emat = nl.eig(Zmat)
1251	Evec = np.sqrt(Evec)/(Natoms-3)
1252	Order = np.argsort(Evec)
1253	XYZ = np.inner(XYZ,Emat.T).T
1254	XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
1255	G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
1256	G2fil.G2Print (' Name X Y Z')
1257	for i,xyz in enumerate(XYZ):
1258	G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
1259	G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
1260
1261	def main():
1262	'Called to run a refinement when this module is executed '
1263	starttime = time.time()
1264	arg = sys.argv
1265	if len(arg) > 1:
1266	GPXfile = arg[1]
1267	if not ospath.exists(GPXfile):
1268	G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
1269	exit()
1270	else:
1271	G2fil.G2Print ('ERROR - missing filename')
1272	exit()
1273	# TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
1274	Refine(GPXfile,None)
1275	G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
1276
1277	if __name__ == '__main__':
1278	GSASIIpath.InvokeDebugOpts()
1279	main()

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

Download in other formats: