Context Navigation

source: trunk/GSASIImath.py @ 1282

Last change on this file since 1282 was 1282, checked in by toby, 9 years ago
sequential refinement updates
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 105.2 KB

Line
1	# -- coding: utf-8 --
2	#GSASIImath - major mathematics routines
3	########### SVN repository information ###################
4	# $Date: 2014-04-18 15:22:29 +0000 (Fri, 18 Apr 2014) $
5	# $Author: toby $
6	# $Revision: 1282 $
7	# $URL: trunk/GSASIImath.py $
8	# $Id: GSASIImath.py 1282 2014-04-18 15:22:29Z toby $
9	########### SVN repository information ###################
10	'''
11	GSASIImath: computation module
12	================================
13
14	Routines for least-squares minimization and other stuff
15
16	'''
17	import sys
18	import os
19	import os.path as ospath
20	import random as rn
21	import numpy as np
22	import numpy.linalg as nl
23	import numpy.ma as ma
24	import cPickle
25	import time
26	import math
27	import copy
28	import GSASIIpath
29	GSASIIpath.SetVersionNumber("$Revision: 1282 $")
30	import GSASIIElem as G2el
31	import GSASIIlattice as G2lat
32	import GSASIIspc as G2spc
33	import GSASIIpwd as G2pwd
34	import numpy.fft as fft
35	import pypowder as pyd
36
37	sind = lambda x: np.sin(x*np.pi/180.)
38	cosd = lambda x: np.cos(x*np.pi/180.)
39	tand = lambda x: np.tan(x*np.pi/180.)
40	asind = lambda x: 180.*np.arcsin(x)/np.pi
41	acosd = lambda x: 180.*np.arccos(x)/np.pi
42	atand = lambda x: 180.*np.arctan(x)/np.pi
43	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
44
45	def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0,Print=False):
46
47	"""
48	Minimize the sum of squares of a function (:math:`f`) evaluated on a series of
49	values (y): :math:`\sum_{y=0}^{N_{obs}} f(y,{args})`
50
51	::
52
53	Nobs
54	x = arg min(sum(func(y)**2,axis=0))
55	y=0
56
57	:param function func: callable method or function
58	should take at least one (possibly length N vector) argument and
59	returns M floating point numbers.
60	:param np.ndarray x0: The starting estimate for the minimization of length N
61	:param function Hess: callable method or function
62	A required function or method to compute the weighted vector and Hessian for func.
63	It must be a symmetric NxN array
64	:param tuple args: Any extra arguments to func are placed in this tuple.
65	:param float ftol: Relative error desired in the sum of squares.
66	:param float xtol: Relative error desired in the approximate solution.
67	:param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
68	until other limits are met (ftol, xtol)
69	:param bool Print: True for printing results (residuals & times) by cycle
70
71	:returns: (x,cov_x,infodict) where
72
73	* x : ndarray
74	The solution (or the result of the last iteration for an unsuccessful
75	call).
76	* cov_x : ndarray
77	Uses the fjac and ipvt optional outputs to construct an
78	estimate of the jacobian around the solution. ``None`` if a
79	singular matrix encountered (indicates very flat curvature in
80	some direction). This matrix must be multiplied by the
81	residual standard deviation to get the covariance of the
82	parameter estimates -- see curve_fit.
83	* infodict : dict
84	a dictionary of optional outputs with the keys:
85
86	* 'fvec' : the function evaluated at the output
87	* 'num cyc':
88	* 'nfev':
89	* 'lamMax':
90	* 'psing':
91
92	"""
93
94	ifConverged = False
95	deltaChi2 = -10.
96	x0 = np.array(x0, ndmin=1) #might be redundant?
97	n = len(x0)
98	if type(args) != type(()):
99	args = (args,)
100
101	icycle = 0
102	One = np.ones((n,n))
103	lam = 0.001
104	lamMax = lam
105	nfev = 0
106	if Print:
107	print ' Hessian refinement on %d variables:'%(n)
108	Lam = np.zeros((n,n))
109	while icycle < maxcyc:
110	time0 = time.time()
111	M = func(x0,*args)
112	nfev += 1
113	chisq0 = np.sum(M**2)
114	Yvec,Amat = Hess(x0,*args)
115	Adiag = np.sqrt(np.diag(Amat))
116	psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
117	if np.any(psing): #hard singularity in matrix
118	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
119	Anorm = np.outer(Adiag,Adiag)
120	Yvec /= Adiag
121	Amat /= Anorm
122	while True:
123	Lam = np.eye(Amat.shape[0])*lam
124	Amatlam = Amat*(One+Lam)
125	try:
126	Xvec = nl.solve(Amatlam,Yvec)
127	except nl.LinAlgError:
128	print 'ouch #1'
129	psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
130	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
131	Xvec /= Adiag
132	M2 = func(x0+Xvec,*args)
133	nfev += 1
134	chisq1 = np.sum(M2**2)
135	if chisq1 > chisq0:
136	lam *= 10.
137	else:
138	x0 += Xvec
139	lam /= 10.
140	break
141	if lam > 10.e5:
142	print 'ouch #3 chisq1 ',chisq1,' stuck > chisq0 ',chisq0
143	break
144	lamMax = max(lamMax,lam)
145	deltaChi2 = (chisq0-chisq1)/chisq0
146	if Print:
147	print ' Cycle: %d, Time: %.2fs, Chi**2: %.3g, Lambda: %.3g, Delta: %.3g'%(
148	icycle,time.time()-time0,chisq1,lam,deltaChi2)
149	if deltaChi2 < ftol:
150	ifConverged = True
151	if Print: print "converged"
152	break
153	icycle += 1
154	else: #after last cycle or if zero cycles
155	M = func(x0,*args)
156	nfev += 1
157	Yvec,Amat = Hess(x0,*args)
158	Adiag = np.sqrt(np.diag(Amat))
159	Anorm = np.outer(Adiag,Adiag)
160	Amatlam = Amat*(One+Lam)/Anorm #scale Amat to Marquardt array
161	try:
162	Bmat = nl.inv(Amatlam)*(One+Lam)/Anorm #rescale Bmat to Marquardt array
163	return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[], 'Converged': ifConverged, 'DelChi2':deltaChi2}]
164	except nl.LinAlgError:
165	print 'ouch #2 linear algebra error in LS'
166	psing = []
167	if maxcyc:
168	psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
169	return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
170
171	def getVCov(varyNames,varyList,covMatrix):
172	'''obtain variance-covariance terms for a set of variables. NB: the varyList
173	and covMatrix were saved by the last least squares refinement so they must match.
174
175	:param list varyNames: variable names to find v-cov matric for
176	:param list varyList: full list of all variables in v-cov matrix
177	:param nparray covMatrix: full variance-covariance matrix from the last
178	least squares refinement
179
180	:returns: nparray vcov: variance-covariance matrix for the variables given
181	in varyNames
182
183	'''
184	vcov = np.zeros((len(varyNames),len(varyNames)))
185	for i1,name1 in enumerate(varyNames):
186	for i2,name2 in enumerate(varyNames):
187	try:
188	vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
189	except ValueError:
190	vcov[i1][i2] = 0.0
191	return vcov
192
193	def FindAtomIndexByIDs(atomData,IDs,Draw=True):
194	'''finds the set of atom array indices for a list of atom IDs. Will search
195	either the Atom table or the drawAtom table.
196
197	:param list atomData: Atom or drawAtom table containting coordinates, etc.
198	:param list IDs: atom IDs to be found
199	:param bool Draw: True if drawAtom table to be searched; False if Atom table
200	is searched
201
202	:returns: list indx: atom (or drawAtom) indices
203
204	'''
205	indx = []
206	for i,atom in enumerate(atomData):
207	if Draw and atom[-3] in IDs:
208	indx.append(i)
209	elif atom[-1] in IDs:
210	indx.append(i)
211	return indx
212
213	def FillAtomLookUp(atomData):
214	'''create a dictionary of atom indexes with atom IDs as keys
215
216	:param list atomData: Atom table to be used
217
218	:returns: dict atomLookUp: dictionary of atom indexes with atom IDs as keys
219
220	'''
221	atomLookUp = {}
222	for iatm,atom in enumerate(atomData):
223	atomLookUp[atom[-1]] = iatm
224	return atomLookUp
225
226	def GetAtomsById(atomData,atomLookUp,IdList):
227	'''gets a list of atoms from Atom table that match a set of atom IDs
228
229	:param list atomData: Atom table to be used
230	:param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
231	:param list IdList: atom IDs to be found
232
233	:returns: list atoms: list of atoms found
234
235	'''
236	atoms = []
237	for id in IdList:
238	atoms.append(atomData[atomLookUp[id]])
239	return atoms
240
241	def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
242	'''gets atom parameters for atoms using atom IDs
243
244	:param list atomData: Atom table to be used
245	:param dict atomLookUp: dictionary of atom indexes with atom IDs as keys
246	:param list IdList: atom IDs to be found
247	:param int itemLoc: pointer to desired 1st item in an atom table entry
248	:param int numItems: number of items to be retrieved
249
250	:returns: type name: description
251
252	'''
253	Items = []
254	if not isinstance(IdList,list):
255	IdList = [IdList,]
256	for id in IdList:
257	if numItems == 1:
258	Items.append(atomData[atomLookUp[id]][itemLoc])
259	else:
260	Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
261	return Items
262
263	def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
264	'''default doc string
265
266	:param type name: description
267
268	:returns: type name: description
269
270	'''
271	pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
272	dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
273	XYZ = []
274	for ind in indx:
275	names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
276	dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
277	XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
278	return XYZ
279
280	def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj): #unfinished & not used
281	'''default doc string
282
283	:param type name: description
284
285	:returns: type name: description
286
287	'''
288	for atom in atomData:
289	XYZ = np.inner(Amat,atom[cx:cx+3])
290	if atom[cia] == 'A':
291	UIJ = atom[cia+2:cia+8]
292
293	def TLS2Uij(xyz,g,Amat,rbObj): #not used anywhere, but could be?
294	'''default doc string
295
296	:param type name: description
297
298	:returns: type name: description
299
300	'''
301	TLStype,TLS = rbObj['ThermalMotion'][:2]
302	Tmat = np.zeros((3,3))
303	Lmat = np.zeros((3,3))
304	Smat = np.zeros((3,3))
305	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
306	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
307	if 'T' in TLStype:
308	Tmat = G2lat.U6toUij(TLS[:6])
309	if 'L' in TLStype:
310	Lmat = G2lat.U6toUij(TLS[6:12])
311	if 'S' in TLStype:
312	Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
313	XYZ = np.inner(Amat,xyz)
314	Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
315	Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
316	beta = np.inner(np.inner(g,Umat),g)
317	return G2lat.UijtoU6(beta)*gvec
318
319	def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj): #not used anywhere, but could be?
320	'''default doc string
321
322	:param type name: description
323
324	:returns: type name: description
325
326	'''
327	cx,ct,cs,cia = atPtrs
328	TLStype,TLS = rbObj['ThermalMotion'][:2]
329	Tmat = np.zeros((3,3))
330	Lmat = np.zeros((3,3))
331	Smat = np.zeros((3,3))
332	G,g = G2lat.A2Gmat(Amat)
333	gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
334	if 'T' in TLStype:
335	Tmat = G2lat.U6toUij(TLS[:6])
336	if 'L' in TLStype:
337	Lmat = G2lat.U6toUij(TLS[6:12])
338	if 'S' in TLStype:
339	Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
340	for atom in atomData:
341	XYZ = np.inner(Amat,atom[cx:cx+3])
342	Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
343	if 'U' in TSLtype:
344	atom[cia+1] = TLS[0]
345	atom[cia] = 'I'
346	else:
347	atom[cia] = 'A'
348	Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
349	beta = np.inner(np.inner(g,Umat),g)
350	atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
351
352	def GetXYZDist(xyz,XYZ,Amat):
353	'''gets distance from position xyz to all XYZ, xyz & XYZ are np.array
354	and are in crystal coordinates; Amat is crystal to Cart matrix
355
356	:param type name: description
357
358	:returns: type name: description
359
360	'''
361	return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
362
363	def getAtomXYZ(atoms,cx):
364	'''default doc string
365
366	:param type name: description
367
368	:returns: type name: description
369
370	'''
371	XYZ = []
372	for atom in atoms:
373	XYZ.append(atom[cx:cx+3])
374	return np.array(XYZ)
375
376	def RotateRBXYZ(Bmat,Cart,oriQ):
377	'''rotate & transform cartesian coordinates to crystallographic ones
378	no translation applied. To be used for numerical derivatives
379
380	:param type name: description
381
382	:returns: type name: description
383
384	'''
385	''' returns crystal coordinates for atoms described by RBObj
386	'''
387	XYZ = np.zeros_like(Cart)
388	for i,xyz in enumerate(Cart):
389	XYZ[i] = np.inner(Bmat,prodQVQ(oriQ,xyz))
390	return XYZ
391
392	def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
393	'''default doc string
394
395	:param type name: description
396
397	:returns: type name: description
398
399	'''
400	''' returns crystal coordinates for atoms described by RBObj
401	'''
402	RBRes = RBData[RBType][RBObj['RBId']]
403	if RBType == 'Vector':
404	vecs = RBRes['rbVect']
405	mags = RBRes['VectMag']
406	Cart = np.zeros_like(vecs[0])
407	for vec,mag in zip(vecs,mags):
408	Cart += vec*mag
409	elif RBType == 'Residue':
410	Cart = np.array(RBRes['rbXYZ'])
411	for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
412	QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
413	Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
414	XYZ = np.zeros_like(Cart)
415	for i,xyz in enumerate(Cart):
416	XYZ[i] = np.inner(Bmat,prodQVQ(RBObj['Orient'][0],xyz))+RBObj['Orig'][0]
417	return XYZ,Cart
418
419	def UpdateMCSAxyz(Bmat,MCSA):
420	'''default doc string
421
422	:param type name: description
423
424	:returns: type name: description
425
426	'''
427	xyz = []
428	atTypes = []
429	iatm = 0
430	for model in MCSA['Models'][1:]: #skip the MD model
431	if model['Type'] == 'Atom':
432	xyz.append(model['Pos'][0])
433	atTypes.append(model['atType'])
434	iatm += 1
435	else:
436	RBRes = MCSA['rbData'][model['Type']][model['RBId']]
437	Pos = np.array(model['Pos'][0])
438	Ori = np.array(model['Ori'][0])
439	Qori = AVdeg2Q(Ori[0],Ori[1:])
440	if model['Type'] == 'Vector':
441	vecs = RBRes['rbVect']
442	mags = RBRes['VectMag']
443	Cart = np.zeros_like(vecs[0])
444	for vec,mag in zip(vecs,mags):
445	Cart += vec*mag
446	elif model['Type'] == 'Residue':
447	Cart = np.array(RBRes['rbXYZ'])
448	for itor,seq in enumerate(RBRes['rbSeq']):
449	QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
450	Cart[seq[3]] = prodQVQ(QuatA,(Cart[seq[3]]-Cart[seq[1]]))+Cart[seq[1]]
451	if model['MolCent'][1]:
452	Cart -= model['MolCent'][0]
453	for i,x in enumerate(Cart):
454	xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
455	atType = RBRes['rbTypes'][i]
456	atTypes.append(atType)
457	iatm += 1
458	return np.array(xyz),atTypes
459
460	def SetMolCent(model,RBData):
461	'''default doc string
462
463	:param type name: description
464
465	:returns: type name: description
466
467	'''
468	rideList = []
469	RBRes = RBData[model['Type']][model['RBId']]
470	if model['Type'] == 'Vector':
471	vecs = RBRes['rbVect']
472	mags = RBRes['VectMag']
473	Cart = np.zeros_like(vecs[0])
474	for vec,mag in zip(vecs,mags):
475	Cart += vec*mag
476	elif model['Type'] == 'Residue':
477	Cart = np.array(RBRes['rbXYZ'])
478	for seq in RBRes['rbSeq']:
479	rideList += seq[3]
480	centList = set(range(len(Cart)))-set(rideList)
481	cent = np.zeros(3)
482	for i in centList:
483	cent += Cart[i]
484	model['MolCent'][0] = cent/len(centList)
485
486	def UpdateRBUIJ(Bmat,Cart,RBObj):
487	'''default doc string
488
489	:param type name: description
490
491	:returns: type name: description
492
493	'''
494	''' returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
495	'''
496	TLStype,TLS = RBObj['ThermalMotion'][:2]
497	T = np.zeros(6)
498	L = np.zeros(6)
499	S = np.zeros(8)
500	if 'T' in TLStype:
501	T = TLS[:6]
502	if 'L' in TLStype:
503	L = np.array(TLS[6:12])(np.pi/180.)*2
504	if 'S' in TLStype:
505	S = np.array(TLS[12:])*(np.pi/180.)
506	g = nl.inv(np.inner(Bmat,Bmat))
507	gvec = np.sqrt(np.array([g[0][0]2,g[1][1]2,g[2][2]**2,
508	g[0][0]g[1][1],g[0][0]g[2][2],g[1][1]*g[2][2]]))
509	Uout = []
510	Q = RBObj['Orient'][0]
511	for X in Cart:
512	X = prodQVQ(Q,X)
513	if 'U' in TLStype:
514	Uout.append(['I',TLS[0],0,0,0,0,0,0])
515	elif not 'N' in TLStype:
516	U = [0,0,0,0,0,0]
517	U[0] = T[0]+L[1]X[2]2+L[2]X[1]*2-2.0L[5]X[1]X[2]+2.0(S[2]X[2]-S[4]*X[1])
518	U[1] = T[1]+L[0]X[2]2+L[2]X[0]*2-2.0L[4]X[0]X[2]+2.0(S[5]X[0]-S[0]*X[2])
519	U[2] = T[2]+L[1]X[0]2+L[0]X[1]*2-2.0L[3]X[1]X[0]+2.0(S[1]X[1]-S[3]*X[0])
520	U[3] = T[3]+L[4]X[1]X[2]+L[5]X[0]X[2]-L[3]X[2]2-L[2]X[0]*X[1]+ \
521	S[4]X[0]-S[5]X[1]-(S[6]+S[7])*X[2]
522	U[4] = T[4]+L[3]X[1]X[2]+L[5]X[0]X[1]-L[4]X[1]2-L[1]X[0]*X[2]+ \
523	S[3]X[2]-S[2]X[0]+S[6]*X[1]
524	U[5] = T[5]+L[3]X[0]X[2]+L[4]X[0]X[1]-L[5]X[0]2-L[0]X[2]*X[1]+ \
525	S[0]X[1]-S[1]X[2]+S[7]*X[0]
526	Umat = G2lat.U6toUij(U)
527	beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
528	Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
529	else:
530	Uout.append(['N',])
531	return Uout
532
533	def GetSHCoeff(pId,parmDict,SHkeys):
534	'''default doc string
535
536	:param type name: description
537
538	:returns: type name: description
539
540	'''
541	SHCoeff = {}
542	for shkey in SHkeys:
543	shname = str(pId)+'::'+shkey
544	SHCoeff[shkey] = parmDict[shname]
545	return SHCoeff
546
547	def getMass(generalData):
548	'''Computes mass of unit cell contents
549
550	:param dict generalData: The General dictionary in Phase
551
552	:returns: float mass: Crystal unit cell mass in AMU.
553
554	'''
555	mass = 0.
556	for i,elem in enumerate(generalData['AtomTypes']):
557	mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
558	return mass
559
560	def getDensity(generalData):
561	'''calculate crystal structure density
562
563	:param dict generalData: The General dictionary in Phase
564
565	:returns: float density: crystal density in gm/cm^3
566
567	'''
568	mass = getMass(generalData)
569	Volume = generalData['Cell'][7]
570	density = mass/(0.6022137*Volume)
571	return density,Volume/mass
572
573	def getWave(Parms):
574	'''returns wavelength from Instrument parameters dictionary
575
576	:param dict Parms: Instrument parameters;
577	must contain:
578	Lam: single wavelength
579	or
580	Lam1: Ka1 radiation wavelength
581
582	:returns: float wave: wavelength
583
584	'''
585	try:
586	return Parms['Lam'][1]
587	except KeyError:
588	return Parms['Lam1'][1]
589
590	def El2Mass(Elements):
591	'''compute molecular weight from Elements
592
593	:param dict Elements: elements in molecular formula;
594	each must contain
595	Num: number of atoms in formula
596	Mass: at. wt.
597
598	:returns: float mass: molecular weight.
599
600	'''
601	mass = 0
602	for El in Elements:
603	mass += Elements[El]['Num']*Elements[El]['Mass']
604	return mass
605
606	def Den2Vol(Elements,density):
607	'''converts density to molecular volume
608
609	:param dict Elements: elements in molecular formula;
610	each must contain
611	Num: number of atoms in formula
612	Mass: at. wt.
613	:param float density: material density in gm/cm^3
614
615	:returns: float volume: molecular volume in A^3
616
617	'''
618	return El2Mass(Elements)/(density*0.6022137)
619
620	def Vol2Den(Elements,volume):
621	'''converts volume to density
622
623	:param dict Elements: elements in molecular formula;
624	each must contain
625	Num: number of atoms in formula
626	Mass: at. wt.
627	:param float volume: molecular volume in A^3
628
629	:returns: float density: material density in gm/cm^3
630
631	'''
632	return El2Mass(Elements)/(volume*0.6022137)
633
634	def El2EstVol(Elements):
635	'''Estimate volume from molecular formula; assumes atom volume = 10A^3
636
637	:param dict Elements: elements in molecular formula;
638	each must contain
639	Num: number of atoms in formula
640
641	:returns: float volume: estimate of molecular volume in A^3
642
643	'''
644	vol = 0
645	for El in Elements:
646	vol += 10.*Elements[El]['Num']
647	return vol
648
649	def XScattDen(Elements,vol,wave=0.):
650	'''Estimate X-ray scattering density from molecular formula & volume;
651	ignores valence, but includes anomalous effects
652
653	:param dict Elements: elements in molecular formula;
654	each element must contain
655	Num: number of atoms in formula
656	Z: atomic number
657	:param float vol: molecular volume in A^3
658	:param float wave: optional wavelength in A
659
660	:returns: float rho: scattering density in 10^10cm^-2;
661	if wave > 0 the includes f' contribution
662	:returns: float mu: if wave>0 absorption coeff in cm^-1 ; otherwise 0
663
664	'''
665	rho = 0
666	mu = 0
667	if wave:
668	Xanom = XAnomAbs(Elements,wave)
669	for El in Elements:
670	f0 = Elements[El]['Z']
671	if wave:
672	f0 += Xanom[El][0]
673	mu += Xanom[El][2]*Elements[El]['Num']
674	rho += Elements[El]['Num']*f0
675	return 28.179rho/vol,0.1mu/vol
676
677	def wavekE(wavekE):
678	'''Convert wavelength to energy & vise versa
679
680	:param float waveKe:wavelength in A or energy in kE
681
682	:returns float waveKe:the other one
683
684	'''
685	return 12.397639/wavekE
686
687	def XAnomAbs(Elements,wave):
688	kE = wavekE(wave)
689	Xanom = {}
690	for El in Elements:
691	Orbs = G2el.GetXsectionCoeff(El)
692	Xanom[El] = G2el.FPcalc(Orbs, kE)
693	return Xanom
694
695
696	################################################################################
697	##### distance, angle, planes, torsion stuff stuff
698	################################################################################
699
700	def getSyXYZ(XYZ,ops,SGData):
701	'''default doc string
702
703	:param type name: description
704
705	:returns: type name: description
706
707	'''
708	XYZout = np.zeros_like(XYZ)
709	for i,[xyz,op] in enumerate(zip(XYZ,ops)):
710	if op == '1':
711	XYZout[i] = xyz
712	else:
713	oprs = op.split('+')
714	unit = [0,0,0]
715	if oprs[1]:
716	unit = np.array(list(eval(oprs[1])))
717	syop =int(oprs[0])
718	inv = syop/abs(syop)
719	syop *= inv
720	cent = syop/100
721	syop %= 100
722	syop -= 1
723	M,T = SGData['SGOps'][syop]
724	XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
725	return XYZout
726
727	def getRestDist(XYZ,Amat):
728	'''default doc string
729
730	:param type name: description
731
732	:returns: type name: description
733
734	'''
735	return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
736
737	def getRestDeriv(Func,XYZ,Amat,ops,SGData):
738	'''default doc string
739
740	:param type name: description
741
742	:returns: type name: description
743
744	'''
745	deriv = np.zeros((len(XYZ),3))
746	dx = 0.00001
747	for j,xyz in enumerate(XYZ):
748	for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
749	XYZ[j] -= x
750	d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
751	XYZ[j] += 2*x
752	d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
753	XYZ[j] -= x
754	deriv[j][i] = (d1-d2)/(2*dx)
755	return deriv.flatten()
756
757	def getRestAngle(XYZ,Amat):
758	'''default doc string
759
760	:param type name: description
761
762	:returns: type name: description
763
764	'''
765
766	def calcVec(Ox,Tx,Amat):
767	return np.inner(Amat,(Tx-Ox))
768
769	VecA = calcVec(XYZ[1],XYZ[0],Amat)
770	VecA /= np.sqrt(np.sum(VecA**2))
771	VecB = calcVec(XYZ[1],XYZ[2],Amat)
772	VecB /= np.sqrt(np.sum(VecB**2))
773	edge = VecB-VecA
774	edge = np.sum(edge**2)
775	angle = (2.-edge)/2.
776	angle = max(angle,-1.)
777	return acosd(angle)
778
779	def getRestPlane(XYZ,Amat):
780	'''default doc string
781
782	:param type name: description
783
784	:returns: type name: description
785
786	'''
787	sumXYZ = np.zeros(3)
788	for xyz in XYZ:
789	sumXYZ += xyz
790	sumXYZ /= len(XYZ)
791	XYZ = np.array(XYZ)-sumXYZ
792	XYZ = np.inner(Amat,XYZ).T
793	Zmat = np.zeros((3,3))
794	for i,xyz in enumerate(XYZ):
795	Zmat += np.outer(xyz.T,xyz)
796	Evec,Emat = nl.eig(Zmat)
797	Evec = np.sqrt(Evec)/(len(XYZ)-3)
798	Order = np.argsort(Evec)
799	return Evec[Order[0]]
800
801	def getRestChiral(XYZ,Amat):
802	'''default doc string
803
804	:param type name: description
805
806	:returns: type name: description
807
808	'''
809	VecA = np.empty((3,3))
810	VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
811	VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
812	VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
813	return nl.det(VecA)
814
815	def getRestTorsion(XYZ,Amat):
816	'''default doc string
817
818	:param type name: description
819
820	:returns: type name: description
821
822	'''
823	VecA = np.empty((3,3))
824	VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
825	VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
826	VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
827	D = nl.det(VecA)
828	Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
829	P12 = np.sum(VecA[0]VecA[1])/(Mag[0]Mag[1])
830	P13 = np.sum(VecA[0]VecA[2])/(Mag[0]Mag[2])
831	P23 = np.sum(VecA[1]VecA[2])/(Mag[1]Mag[2])
832	Ang = 1.0
833	if abs(P12) < 1.0 and abs(P23) < 1.0:
834	Ang = (P12P23-P13)/(np.sqrt(1.-P122)np.sqrt(1.-P23**2))
835	TOR = (acosd(Ang)*D/abs(D)+720.)%360.
836	return TOR
837
838	def calcTorsionEnergy(TOR,Coeff=[]):
839	'''default doc string
840
841	:param type name: description
842
843	:returns: type name: description
844
845	'''
846	sum = 0.
847	if len(Coeff):
848	cof = np.reshape(Coeff,(3,3)).T
849	delt = TOR-cof[1]
850	delt = np.where(delt<-180.,delt+360.,delt)
851	delt = np.where(delt>180.,delt-360.,delt)
852	term = -cof[2]delt*2
853	val = cof[0]*np.exp(term/1000.0)
854	pMax = cof[0][np.argmin(val)]
855	Eval = np.sum(val)
856	sum = Eval-pMax
857	return sum,Eval
858
859	def getTorsionDeriv(XYZ,Amat,Coeff):
860	'''default doc string
861
862	:param type name: description
863
864	:returns: type name: description
865
866	'''
867	deriv = np.zeros((len(XYZ),3))
868	dx = 0.00001
869	for j,xyz in enumerate(XYZ):
870	for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
871	XYZ[j] -= x
872	tor = getRestTorsion(XYZ,Amat)
873	p1,d1 = calcTorsionEnergy(tor,Coeff)
874	XYZ[j] += 2*x
875	tor = getRestTorsion(XYZ,Amat)
876	p2,d2 = calcTorsionEnergy(tor,Coeff)
877	XYZ[j] -= x
878	deriv[j][i] = (p2-p1)/(2*dx)
879	return deriv.flatten()
880
881	def getRestRama(XYZ,Amat):
882	'''Computes a pair of torsion angles in a 5 atom string
883
884	:param nparray XYZ: crystallographic coordinates of 5 atoms
885	:param nparray Amat: crystal to cartesian transformation matrix
886
887	:returns: list (phi,psi) two torsion angles in degrees
888
889	'''
890	phi = getRestTorsion(XYZ[:5],Amat)
891	psi = getRestTorsion(XYZ[1:],Amat)
892	return phi,psi
893
894	def calcRamaEnergy(phi,psi,Coeff=[]):
895	'''Computes pseudo potential energy from a pair of torsion angles and a
896	numerical description of the potential energy surface. Used to create
897	penalty function in LS refinement:
898	:math:`Eval(\\phi,\\psi) = C[0]*exp(-V/1000)`
899
900	where :math:`V = -C[3] * (\\phi-C[1])^2 - C[4](\\psi-C[2])^2 - 2(\\phi-C[1])*(\\psi-C[2])`
901
902	:param float phi: first torsion angle (:math:`\\phi`)
903	:param float psi: second torsion angle (:math:`\\psi`)
904	:param list Coeff: pseudo potential coefficients
905
906	:returns: list (sum,Eval): pseudo-potential difference from minimum & value;
907	sum is used for penalty function.
908
909	'''
910	sum = 0.
911	Eval = 0.
912	if len(Coeff):
913	cof = Coeff.T
914	dPhi = phi-cof[1]
915	dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
916	dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
917	dPsi = psi-cof[2]
918	dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
919	dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
920	val = -cof[3]dPhi2-cof[4]dPsi*2-2.0cof[5]dPhidPsi
921	val = cof[0]*np.exp(val/1000.)
922	pMax = cof[0][np.argmin(val)]
923	Eval = np.sum(val)
924	sum = Eval-pMax
925	return sum,Eval
926
927	def getRamaDeriv(XYZ,Amat,Coeff):
928	'''Computes numerical derivatives of torsion angle pair pseudo potential
929	with respect of crystallographic atom coordinates of the 5 atom sequence
930
931	:param nparray XYZ: crystallographic coordinates of 5 atoms
932	:param nparray Amat: crystal to cartesian transformation matrix
933	:param list Coeff: pseudo potential coefficients
934
935	:returns: list (deriv) derivatives of pseudopotential with respect to 5 atom
936	crystallographic xyz coordinates.
937
938	'''
939	deriv = np.zeros((len(XYZ),3))
940	dx = 0.00001
941	for j,xyz in enumerate(XYZ):
942	for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
943	XYZ[j] -= x
944	phi,psi = getRestRama(XYZ,Amat)
945	p1,d1 = calcRamaEnergy(phi,psi,Coeff)
946	XYZ[j] += 2*x
947	phi,psi = getRestRama(XYZ,Amat)
948	p2,d2 = calcRamaEnergy(phi,psi,Coeff)
949	XYZ[j] -= x
950	deriv[j][i] = (p2-p1)/(2*dx)
951	return deriv.flatten()
952
953	def getRestPolefig(ODFln,SamSym,Grid):
954	'''default doc string
955
956	:param type name: description
957
958	:returns: type name: description
959
960	'''
961	X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
962	R,P = np.sqrt(X2+Y2).flatten(),atan2d(Y,X).flatten()
963	R = np.where(R <= 1.,2.*atand(R),0.0)
964	Z = np.zeros_like(R)
965	Z = G2lat.polfcal(ODFln,SamSym,R,P)
966	Z = np.reshape(Z,(Grid,Grid))
967	return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
968
969	def getRestPolefigDerv(HKL,Grid,SHCoeff):
970	'''default doc string
971
972	:param type name: description
973
974	:returns: type name: description
975
976	'''
977	pass
978
979	def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
980	'''default doc string
981
982	:param type name: description
983
984	:returns: type name: description
985
986	'''
987	def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
988	TxT = inv*(np.inner(M,Tx)+T)+C+U
989	return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
990
991	inv = Top/abs(Top)
992	cent = abs(Top)/100
993	op = abs(Top)%100-1
994	M,T = SGData['SGOps'][op]
995	C = SGData['SGCen'][cent]
996	dx = .00001
997	deriv = np.zeros(6)
998	for i in [0,1,2]:
999	Oxyz[i] -= dx
1000	d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1001	Oxyz[i] += 2*dx
1002	deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1003	Oxyz[i] -= dx
1004	Txyz[i] -= dx
1005	d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
1006	Txyz[i] += 2*dx
1007	deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
1008	Txyz[i] -= dx
1009	return deriv
1010
1011	def getAngSig(VA,VB,Amat,SGData,covData={}):
1012	'''default doc string
1013
1014	:param type name: description
1015
1016	:returns: type name: description
1017
1018	'''
1019	def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
1020	TxT = inv*(np.inner(M,Tx)+T)+C
1021	TxT = G2spc.MoveToUnitCell(TxT)+U
1022	return np.inner(Amat,(TxT-Ox))
1023
1024	def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
1025	VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
1026	VecA /= np.sqrt(np.sum(VecA**2))
1027	VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
1028	VecB /= np.sqrt(np.sum(VecB**2))
1029	edge = VecB-VecA
1030	edge = np.sum(edge**2)
1031	angle = (2.-edge)/2.
1032	angle = max(angle,-1.)
1033	return acosd(angle)
1034
1035	OxAN,OxA,TxAN,TxA,unitA,TopA = VA
1036	OxBN,OxB,TxBN,TxB,unitB,TopB = VB
1037	invA = invB = 1
1038	invA = TopA/abs(TopA)
1039	invB = TopB/abs(TopB)
1040	centA = abs(TopA)/100
1041	centB = abs(TopB)/100
1042	opA = abs(TopA)%100-1
1043	opB = abs(TopB)%100-1
1044	MA,TA = SGData['SGOps'][opA]
1045	MB,TB = SGData['SGOps'][opB]
1046	CA = SGData['SGCen'][centA]
1047	CB = SGData['SGCen'][centB]
1048	if 'covMatrix' in covData:
1049	covMatrix = covData['covMatrix']
1050	varyList = covData['varyList']
1051	AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
1052	dx = .00001
1053	dadx = np.zeros(9)
1054	Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1055	for i in [0,1,2]:
1056	OxA[i] -= dx
1057	a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1058	OxA[i] += 2*dx
1059	dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1060	OxA[i] -= dx
1061
1062	TxA[i] -= dx
1063	a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1064	TxA[i] += 2*dx
1065	dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1066	TxA[i] -= dx
1067
1068	TxB[i] -= dx
1069	a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
1070	TxB[i] += 2*dx
1071	dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
1072	TxB[i] -= dx
1073
1074	sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1075	if sigAng < 0.01:
1076	sigAng = 0.0
1077	return Ang,sigAng
1078	else:
1079	return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
1080
1081	def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
1082	'''default doc string
1083
1084	:param type name: description
1085
1086	:returns: type name: description
1087
1088	'''
1089	def calcDist(Atoms,SyOps,Amat):
1090	XYZ = []
1091	for i,atom in enumerate(Atoms):
1092	Inv,M,T,C,U = SyOps[i]
1093	XYZ.append(np.array(atom[1:4]))
1094	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1095	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1096	V1 = XYZ[1]-XYZ[0]
1097	return np.sqrt(np.sum(V1**2))
1098
1099	Inv = []
1100	SyOps = []
1101	names = []
1102	for i,atom in enumerate(Oatoms):
1103	names += atom[-1]
1104	Op,unit = Atoms[i][-1]
1105	inv = Op/abs(Op)
1106	m,t = SGData['SGOps'][abs(Op)%100-1]
1107	c = SGData['SGCen'][abs(Op)/100]
1108	SyOps.append([inv,m,t,c,unit])
1109	Dist = calcDist(Oatoms,SyOps,Amat)
1110
1111	sig = -0.001
1112	if 'covMatrix' in covData:
1113	parmNames = []
1114	dx = .00001
1115	dadx = np.zeros(6)
1116	for i in range(6):
1117	ia = i/3
1118	ix = i%3
1119	Oatoms[ia][ix+1] += dx
1120	a0 = calcDist(Oatoms,SyOps,Amat)
1121	Oatoms[ia][ix+1] -= 2*dx
1122	dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1123	covMatrix = covData['covMatrix']
1124	varyList = covData['varyList']
1125	DistVcov = getVCov(names,varyList,covMatrix)
1126	sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
1127	if sig < 0.001:
1128	sig = -0.001
1129
1130	return Dist,sig
1131
1132	def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
1133	'''default doc string
1134
1135	:param type name: description
1136
1137	:returns: type name: description
1138
1139	'''
1140
1141	def calcAngle(Atoms,SyOps,Amat):
1142	XYZ = []
1143	for i,atom in enumerate(Atoms):
1144	Inv,M,T,C,U = SyOps[i]
1145	XYZ.append(np.array(atom[1:4]))
1146	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1147	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1148	V1 = XYZ[1]-XYZ[0]
1149	V1 /= np.sqrt(np.sum(V1**2))
1150	V2 = XYZ[1]-XYZ[2]
1151	V2 /= np.sqrt(np.sum(V2**2))
1152	V3 = V2-V1
1153	cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1154	return acosd(cang)
1155
1156	Inv = []
1157	SyOps = []
1158	names = []
1159	for i,atom in enumerate(Oatoms):
1160	names += atom[-1]
1161	Op,unit = Atoms[i][-1]
1162	inv = Op/abs(Op)
1163	m,t = SGData['SGOps'][abs(Op)%100-1]
1164	c = SGData['SGCen'][abs(Op)/100]
1165	SyOps.append([inv,m,t,c,unit])
1166	Angle = calcAngle(Oatoms,SyOps,Amat)
1167
1168	sig = -0.01
1169	if 'covMatrix' in covData:
1170	parmNames = []
1171	dx = .00001
1172	dadx = np.zeros(9)
1173	for i in range(9):
1174	ia = i/3
1175	ix = i%3
1176	Oatoms[ia][ix+1] += dx
1177	a0 = calcAngle(Oatoms,SyOps,Amat)
1178	Oatoms[ia][ix+1] -= 2*dx
1179	dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1180	covMatrix = covData['covMatrix']
1181	varyList = covData['varyList']
1182	AngVcov = getVCov(names,varyList,covMatrix)
1183	sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
1184	if sig < 0.01:
1185	sig = -0.01
1186
1187	return Angle,sig
1188
1189	def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
1190	'''default doc string
1191
1192	:param type name: description
1193
1194	:returns: type name: description
1195
1196	'''
1197
1198	def calcTorsion(Atoms,SyOps,Amat):
1199
1200	XYZ = []
1201	for i,atom in enumerate(Atoms):
1202	Inv,M,T,C,U = SyOps[i]
1203	XYZ.append(np.array(atom[1:4]))
1204	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1205	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1206	V1 = XYZ[1]-XYZ[0]
1207	V2 = XYZ[2]-XYZ[1]
1208	V3 = XYZ[3]-XYZ[2]
1209	V1 /= np.sqrt(np.sum(V1**2))
1210	V2 /= np.sqrt(np.sum(V2**2))
1211	V3 /= np.sqrt(np.sum(V3**2))
1212	M = np.array([V1,V2,V3])
1213	D = nl.det(M)
1214	Ang = 1.0
1215	P12 = np.dot(V1,V2)
1216	P13 = np.dot(V1,V3)
1217	P23 = np.dot(V2,V3)
1218	Tors = acosd((P12P23-P13)/(np.sqrt(1.-P122)np.sqrt(1.-P23*2)))D/abs(D)
1219	return Tors
1220
1221	Inv = []
1222	SyOps = []
1223	names = []
1224	for i,atom in enumerate(Oatoms):
1225	names += atom[-1]
1226	Op,unit = Atoms[i][-1]
1227	inv = Op/abs(Op)
1228	m,t = SGData['SGOps'][abs(Op)%100-1]
1229	c = SGData['SGCen'][abs(Op)/100]
1230	SyOps.append([inv,m,t,c,unit])
1231	Tors = calcTorsion(Oatoms,SyOps,Amat)
1232
1233	sig = -0.01
1234	if 'covMatrix' in covData:
1235	parmNames = []
1236	dx = .00001
1237	dadx = np.zeros(12)
1238	for i in range(12):
1239	ia = i/3
1240	ix = i%3
1241	Oatoms[ia][ix+1] -= dx
1242	a0 = calcTorsion(Oatoms,SyOps,Amat)
1243	Oatoms[ia][ix+1] += 2*dx
1244	dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1245	Oatoms[ia][ix+1] -= dx
1246	covMatrix = covData['covMatrix']
1247	varyList = covData['varyList']
1248	TorVcov = getVCov(names,varyList,covMatrix)
1249	sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
1250	if sig < 0.01:
1251	sig = -0.01
1252
1253	return Tors,sig
1254
1255	def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
1256	'''default doc string
1257
1258	:param type name: description
1259
1260	:returns: type name: description
1261
1262	'''
1263
1264	def calcDist(Atoms,SyOps,Amat):
1265	XYZ = []
1266	for i,atom in enumerate(Atoms):
1267	Inv,M,T,C,U = SyOps[i]
1268	XYZ.append(np.array(atom[1:4]))
1269	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1270	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1271	V1 = XYZ[1]-XYZ[0]
1272	return np.sqrt(np.sum(V1**2))
1273
1274	def calcAngle(Atoms,SyOps,Amat):
1275	XYZ = []
1276	for i,atom in enumerate(Atoms):
1277	Inv,M,T,C,U = SyOps[i]
1278	XYZ.append(np.array(atom[1:4]))
1279	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1280	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1281	V1 = XYZ[1]-XYZ[0]
1282	V1 /= np.sqrt(np.sum(V1**2))
1283	V2 = XYZ[1]-XYZ[2]
1284	V2 /= np.sqrt(np.sum(V2**2))
1285	V3 = V2-V1
1286	cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
1287	return acosd(cang)
1288
1289	def calcTorsion(Atoms,SyOps,Amat):
1290
1291	XYZ = []
1292	for i,atom in enumerate(Atoms):
1293	Inv,M,T,C,U = SyOps[i]
1294	XYZ.append(np.array(atom[1:4]))
1295	XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
1296	XYZ[-1] = np.inner(Amat,XYZ[-1]).T
1297	V1 = XYZ[1]-XYZ[0]
1298	V2 = XYZ[2]-XYZ[1]
1299	V3 = XYZ[3]-XYZ[2]
1300	V1 /= np.sqrt(np.sum(V1**2))
1301	V2 /= np.sqrt(np.sum(V2**2))
1302	V3 /= np.sqrt(np.sum(V3**2))
1303	M = np.array([V1,V2,V3])
1304	D = nl.det(M)
1305	Ang = 1.0
1306	P12 = np.dot(V1,V2)
1307	P13 = np.dot(V1,V3)
1308	P23 = np.dot(V2,V3)
1309	Tors = acosd((P12P23-P13)/(np.sqrt(1.-P122)np.sqrt(1.-P23*2)))D/abs(D)
1310	return Tors
1311
1312	Inv = []
1313	SyOps = []
1314	names = []
1315	for i,atom in enumerate(Oatoms):
1316	names += atom[-1]
1317	Op,unit = Atoms[i][-1]
1318	inv = Op/abs(Op)
1319	m,t = SGData['SGOps'][abs(Op)%100-1]
1320	c = SGData['SGCen'][abs(Op)/100]
1321	SyOps.append([inv,m,t,c,unit])
1322	M = len(Oatoms)
1323	if M == 2:
1324	Val = calcDist(Oatoms,SyOps,Amat)
1325	elif M == 3:
1326	Val = calcAngle(Oatoms,SyOps,Amat)
1327	else:
1328	Val = calcTorsion(Oatoms,SyOps,Amat)
1329
1330	sigVals = [-0.001,-0.01,-0.01]
1331	sig = sigVals[M-3]
1332	if 'covMatrix' in covData:
1333	parmNames = []
1334	dx = .00001
1335	N = M*3
1336	dadx = np.zeros(N)
1337	for i in range(N):
1338	ia = i/3
1339	ix = i%3
1340	Oatoms[ia][ix+1] += dx
1341	if M == 2:
1342	a0 = calcDist(Oatoms,SyOps,Amat)
1343	elif M == 3:
1344	a0 = calcAngle(Oatoms,SyOps,Amat)
1345	else:
1346	a0 = calcTorsion(Oatoms,SyOps,Amat)
1347	Oatoms[ia][ix+1] -= 2*dx
1348	if M == 2:
1349	dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1350	elif M == 3:
1351	dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1352	else:
1353	dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
1354	covMatrix = covData['covMatrix']
1355	varyList = covData['varyList']
1356	Vcov = getVCov(names,varyList,covMatrix)
1357	sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
1358	if sig < sigVals[M-3]:
1359	sig = sigVals[M-3]
1360
1361	return Val,sig
1362
1363	def ValEsd(value,esd=0,nTZ=False):
1364	'''Format a floating point number with a given level of precision or
1365	with in crystallographic format with a "esd", as value(esd). If esd is
1366	negative the number is formatted with the level of significant figures
1367	appropriate if abs(esd) were the esd, but the esd is not included.
1368	if the esd is zero, approximately 6 significant figures are printed.
1369	nTZ=True causes "extra" zeros to be removed after the decimal place.
1370	for example:
1371
1372	* "1.235(3)" for value=1.2346 & esd=0.003
1373	* "1.235(3)e4" for value=12346. & esd=30
1374	* "1.235(3)e6" for value=0.12346e7 & esd=3000
1375	* "1.235" for value=1.2346 & esd=-0.003
1376	* "1.240" for value=1.2395 & esd=-0.003
1377	* "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
1378	* "1.23460" for value=1.2346 & esd=0.0
1379
1380	:param float value: number to be formatted
1381	:param float esd: uncertainty or if esd < 0, specifies level of
1382	precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
1383	:param bool nTZ: True to remove trailing zeros (default is False)
1384	:returns: value(esd) or value as a string
1385
1386	'''
1387	# Note: this routine is Python 3 compatible -- I think
1388	if math.isnan(value): # invalid value, bail out
1389	return '?'
1390	if math.isnan(esd): # invalid esd, treat as zero
1391	esd = 0
1392	esdoff = 5
1393	elif esd != 0:
1394	# transform the esd to a one or two digit integer
1395	l = math.log10(abs(esd)) % 1
1396	# cut off of 19 1.9==>(19) but 1.95==>(2) N.B. log10(1.95) = 0.2900...
1397	if l < 0.290034611362518: l+= 1.
1398	intesd = int(round(10**l)) # esd as integer
1399	# determine the number of digits offset for the esd
1400	esdoff = int(round(math.log10(intesd*1./abs(esd))))
1401	else:
1402	esdoff = 5
1403	valoff = 0
1404	if abs(value) < abs(esdoff): # value is effectively zero
1405	pass
1406	elif esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
1407	# where the digit offset is to the left of the decimal place or where too many
1408	# digits are needed
1409	if abs(value) > 1:
1410	valoff = int(math.log10(abs(value)))
1411	elif abs(value) > 0:
1412	valoff = int(math.log10(abs(value))-0.9999999)
1413	else:
1414	valoff = 0
1415	if esd != 0:
1416	if valoff+esdoff < 0:
1417	valoff = esdoff = 0
1418	out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
1419	elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
1420	out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
1421	else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
1422	if abs(value) > 0:
1423	extra = -math.log10(abs(value))
1424	else:
1425	extra = 0
1426	if extra > 0: extra += 1
1427	out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
1428	if esd > 0:
1429	out += ("({:d})").format(intesd) # add the esd
1430	elif nTZ and '.' in out:
1431	out = out.rstrip('0') # strip zeros to right of decimal
1432	out = out.rstrip('.') # and decimal place when not needed
1433	if valoff != 0:
1434	out += ("e{:d}").format(valoff) # add an exponent, when needed
1435	return out
1436
1437	################################################################################
1438	##### Fourier & charge flip stuff
1439	################################################################################
1440
1441	def adjHKLmax(SGData,Hmax):
1442	'''default doc string
1443
1444	:param type name: description
1445
1446	:returns: type name: description
1447
1448	'''
1449	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1450	Hmax[0] = ((Hmax[0]+3)/6)*6
1451	Hmax[1] = ((Hmax[1]+3)/6)*6
1452	Hmax[2] = ((Hmax[2]+1)/4)*4
1453	else:
1454	Hmax[0] = ((Hmax[0]+2)/4)*4
1455	Hmax[1] = ((Hmax[1]+2)/4)*4
1456	Hmax[2] = ((Hmax[2]+1)/4)*4
1457
1458	def OmitMap(data,reflDict,pgbar=None):
1459	'''default doc string
1460
1461	:param type name: description
1462
1463	:returns: type name: description
1464
1465	'''
1466	generalData = data['General']
1467	if not generalData['Map']['MapType']:
1468	print '**** ERROR - Fourier map not defined'
1469	return
1470	mapData = generalData['Map']
1471	dmin = mapData['Resolution']
1472	SGData = generalData['SGData']
1473	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1474	SGT = np.array([ops[1] for ops in SGData['SGOps']])
1475	cell = generalData['Cell'][1:8]
1476	A = G2lat.cell2A(cell[:6])
1477	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1478	adjHKLmax(SGData,Hmax)
1479	Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1480	time0 = time.time()
1481	for iref,ref in enumerate(reflDict['RefList']):
1482	if ref[4] >= dmin:
1483	Fosq,Fcsq,ph = ref[8:11]
1484	Uniq = np.inner(ref[:3],SGMT)
1485	Phi = np.inner(ref[:3],SGT)
1486	for i,hkl in enumerate(Uniq): #uses uniq
1487	hkl = np.asarray(hkl,dtype='i')
1488	dp = 360.*Phi[i] #and phi
1489	a = cosd(ph+dp)
1490	b = sind(ph+dp)
1491	phasep = complex(a,b)
1492	phasem = complex(a,-b)
1493	Fo = np.sqrt(Fosq)
1494	if '2Fo-Fc' in mapData['MapType']:
1495	F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1496	else:
1497	F = np.sqrt(Fosq)
1498	h,k,l = hkl+Hmax
1499	Fhkl[h,k,l] = F*phasep
1500	h,k,l = -hkl+Hmax
1501	Fhkl[h,k,l] = F*phasem
1502	rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1503	M = np.mgrid[0:4,0:4,0:4]
1504	blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
1505	iBeg = blkIds*rho0.shape/4
1506	iFin = (blkIds+1)*rho0.shape/4
1507	rho_omit = np.zeros_like(rho0)
1508	nBlk = 0
1509	for iB,iF in zip(iBeg,iFin):
1510	rho1 = np.copy(rho0)
1511	rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
1512	Fnew = fft.ifftshift(fft.ifftn(rho1))
1513	Fnew = np.where(Fnew,Fnew,1.0) #avoid divide by zero
1514	phase = Fnew/np.absolute(Fnew)
1515	OFhkl = np.absolute(Fhkl)*phase
1516	rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
1517	rho_omit[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = np.copy(rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]])
1518	nBlk += 1
1519	pgbar.Update(nBlk)
1520	mapData['rho'] = np.real(rho_omit)/cell[6]
1521	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1522	print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1523	return mapData
1524
1525	def FourierMap(data,reflDict):
1526	'''default doc string
1527
1528	:param type name: description
1529
1530	:returns: type name: description
1531
1532	'''
1533	generalData = data['General']
1534	if not generalData['Map']['MapType']:
1535	print '**** ERROR - Fourier map not defined'
1536	return
1537	mapData = generalData['Map']
1538	dmin = mapData['Resolution']
1539	SGData = generalData['SGData']
1540	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1541	SGT = np.array([ops[1] for ops in SGData['SGOps']])
1542	cell = generalData['Cell'][1:8]
1543	A = G2lat.cell2A(cell[:6])
1544	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1545	adjHKLmax(SGData,Hmax)
1546	Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1547	# Fhkl[0,0,0] = generalData['F000X']
1548	time0 = time.time()
1549	for iref,ref in enumerate(reflDict['RefList']):
1550	if ref[4] >= dmin:
1551	Fosq,Fcsq,ph = ref[8:11]
1552	Uniq = np.inner(ref[:3],SGMT)
1553	Phi = np.inner(ref[:3],SGT)
1554	for i,hkl in enumerate(Uniq): #uses uniq
1555	hkl = np.asarray(hkl,dtype='i')
1556	dp = 360.*Phi[i] #and phi
1557	a = cosd(ph+dp)
1558	b = sind(ph+dp)
1559	phasep = complex(a,b)
1560	phasem = complex(a,-b)
1561	if 'Fobs' in mapData['MapType']:
1562	F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1563	h,k,l = hkl+Hmax
1564	Fhkl[h,k,l] = F*phasep
1565	h,k,l = -hkl+Hmax
1566	Fhkl[h,k,l] = F*phasem
1567	elif 'Fcalc' in mapData['MapType']:
1568	F = np.sqrt(Fcsq)
1569	h,k,l = hkl+Hmax
1570	Fhkl[h,k,l] = F*phasep
1571	h,k,l = -hkl+Hmax
1572	Fhkl[h,k,l] = F*phasem
1573	elif 'delt-F' in mapData['MapType']:
1574	dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1575	h,k,l = hkl+Hmax
1576	Fhkl[h,k,l] = dF*phasep
1577	h,k,l = -hkl+Hmax
1578	Fhkl[h,k,l] = dF*phasem
1579	elif '2*Fo-Fc' in mapData['MapType']:
1580	F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1581	h,k,l = hkl+Hmax
1582	Fhkl[h,k,l] = F*phasep
1583	h,k,l = -hkl+Hmax
1584	Fhkl[h,k,l] = F*phasem
1585	elif 'Patterson' in mapData['MapType']:
1586	h,k,l = hkl+Hmax
1587	Fhkl[h,k,l] = complex(Fosq,0.)
1588	h,k,l = -hkl+Hmax
1589	Fhkl[h,k,l] = complex(Fosq,0.)
1590	rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1591	print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1592	mapData['rho'] = np.real(rho)
1593	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1594	return mapData
1595
1596	# map printing for testing purposes
1597	def printRho(SGLaue,rho,rhoMax):
1598	'''default doc string
1599
1600	:param type name: description
1601
1602	:returns: type name: description
1603
1604	'''
1605	dim = len(rho.shape)
1606	if dim == 2:
1607	ix,jy = rho.shape
1608	for j in range(jy):
1609	line = ''
1610	if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1611	line += (jy-j)*' '
1612	for i in range(ix):
1613	r = int(100*rho[i,j]/rhoMax)
1614	line += '%4d'%(r)
1615	print line+'\n'
1616	else:
1617	ix,jy,kz = rho.shape
1618	for k in range(kz):
1619	print 'k = ',k
1620	for j in range(jy):
1621	line = ''
1622	if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1623	line += (jy-j)*' '
1624	for i in range(ix):
1625	r = int(100*rho[i,j,k]/rhoMax)
1626	line += '%4d'%(r)
1627	print line+'\n'
1628	## keep this
1629
1630	def findOffset(SGData,A,Fhkl):
1631	'''default doc string
1632
1633	:param type name: description
1634
1635	:returns: type name: description
1636
1637	'''
1638	if SGData['SpGrp'] == 'P 1':
1639	return [0,0,0]
1640	hklShape = Fhkl.shape
1641	hklHalf = np.array(hklShape)/2
1642	sortHKL = np.argsort(Fhkl.flatten())
1643	Fdict = {}
1644	for hkl in sortHKL:
1645	HKL = np.unravel_index(hkl,hklShape)
1646	F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1647	if F == 0.:
1648	break
1649	Fdict['%.6f'%(np.absolute(F))] = hkl
1650	Flist = np.flipud(np.sort(Fdict.keys()))
1651	F = str(1.e6)
1652	i = 0
1653	DH = []
1654	Dphi = []
1655	Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1656	for F in Flist:
1657	hkl = np.unravel_index(Fdict[F],hklShape)
1658	iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1659	Uniq = np.array(Uniq,dtype='i')
1660	Phi = np.array(Phi)
1661	Uniq = np.concatenate((Uniq,-Uniq))+hklHalf # put in Friedel pairs & make as index to Farray
1662	Phi = np.concatenate((Phi,-Phi)) # and their phase shifts
1663	Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1664	ang0 = np.angle(Fh0,deg=True)/360.
1665	for H,phi in zip(Uniq,Phi)[1:]:
1666	ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
1667	dH = H-hkl
1668	dang = ang-ang0
1669	if np.any(np.abs(dH)- Hmax > 0): #keep low order DHs
1670	continue
1671	DH.append(dH)
1672	Dphi.append((dang+.5) % 1.0)
1673	if i > 20 or len(DH) > 30:
1674	break
1675	i += 1
1676	DH = np.array(DH)
1677	print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1678	Dphi = np.array(Dphi)
1679	steps = np.array(hklShape)
1680	X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1681	XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1682	Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
1683	Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1684	hist,bins = np.histogram(Mmap,bins=1000)
1685	# for i,item in enumerate(hist[:10]):
1686	# print item,bins[i]
1687	chisq = np.min(Mmap)
1688	DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1689	print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1690	# print (np.dot(DX,DH.T)+.5)%1.-Dphi
1691	return DX
1692
1693	def ChargeFlip(data,reflDict,pgbar):
1694	'''default doc string
1695
1696	:param type name: description
1697
1698	:returns: type name: description
1699
1700	'''
1701	generalData = data['General']
1702	mapData = generalData['Map']
1703	flipData = generalData['Flip']
1704	FFtable = {}
1705	if 'None' not in flipData['Norm element']:
1706	normElem = flipData['Norm element'].upper()
1707	FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1708	for ff in FFs:
1709	if ff['Symbol'] == normElem:
1710	FFtable.update(ff)
1711	dmin = flipData['Resolution']
1712	SGData = generalData['SGData']
1713	SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1714	SGT = np.array([ops[1] for ops in SGData['SGOps']])
1715	cell = generalData['Cell'][1:8]
1716	A = G2lat.cell2A(cell[:6])
1717	Vol = cell[6]
1718	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1719	adjHKLmax(SGData,Hmax)
1720	Ehkl = np.zeros(shape=2*Hmax,dtype='c16') #2X64bits per complex no.
1721	time0 = time.time()
1722	for iref,ref in enumerate(reflDict['RefList']):
1723	dsp = ref[4]
1724	if dsp >= dmin:
1725	ff = 0.1Vol #est. no. atoms for ~10A*3/atom
1726	if FFtable:
1727	SQ = 0.25/dsp**2
1728	ff *= G2el.ScatFac(FFtable,SQ)[0]
1729	if ref[8] > 0.: #use only +ve Fobs**2
1730	E = np.sqrt(ref[8])/ff
1731	else:
1732	E = 0.
1733	ph = ref[10]
1734	ph = rn.uniform(0.,360.)
1735	Uniq = np.inner(ref[:3],SGMT)
1736	Phi = np.inner(ref[:3],SGT)
1737	for i,hkl in enumerate(Uniq): #uses uniq
1738	hkl = np.asarray(hkl,dtype='i')
1739	dp = 360.*Phi[i] #and phi
1740	a = cosd(ph+dp)
1741	b = sind(ph+dp)
1742	phasep = complex(a,b)
1743	phasem = complex(a,-b)
1744	h,k,l = hkl+Hmax
1745	Ehkl[h,k,l] = E*phasep
1746	h,k,l = -hkl+Hmax #Friedel pair refl.
1747	Ehkl[h,k,l] = E*phasem
1748	# Ehkl[Hmax] = 0.00001 #this to preserve F[0,0,0]
1749	CEhkl = copy.copy(Ehkl)
1750	MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1751	Emask = ma.getmask(MEhkl)
1752	sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1753	Ncyc = 0
1754	old = np.seterr(all='raise')
1755	while True:
1756	CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1757	CEsig = np.std(CErho)
1758	CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1759	CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho) #solves U atom problem!
1760	CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1761	CFhkl = np.where(CFhkl,CFhkl,1.0) #avoid divide by zero
1762	phase = CFhkl/np.absolute(CFhkl)
1763	CEhkl = np.absolute(Ehkl)*phase
1764	Ncyc += 1
1765	sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1766	DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1767	Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1768	if Rcf < 5.:
1769	break
1770	GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1771	if not GoOn or Ncyc > 10000:
1772	break
1773	np.seterr(**old)
1774	print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1775	CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1776	print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1777	roll = findOffset(SGData,A,CEhkl) #CEhkl needs to be just the observed set, not the full set!
1778
1779	mapData['Rcf'] = Rcf
1780	mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1781	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1782	mapData['rollMap'] = [0,0,0]
1783	return mapData
1784
1785	def SearchMap(data):
1786	'''Does a search of a density map for peaks meeting the criterion of peak
1787	height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
1788	mapData is data['General']['mapData']; the map is also in mapData.
1789
1790	:param data: the phase data structure
1791
1792	:returns: (peaks,mags,dzeros) where
1793
1794	* peaks : ndarray
1795	x,y,z positions of the peaks found in the map
1796	* mags : ndarray
1797	the magnitudes of the peaks
1798	* dzeros : ndarray
1799	the distance of the peaks from the unit cell origin
1800
1801	'''
1802	rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1803
1804	norm = 1./(np.sqrt(3.)np.sqrt(2.np.pi)**3)
1805
1806	# def noDuplicate(xyz,peaks,Amat):
1807	# XYZ = np.inner(Amat,xyz)
1808	# if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1809	# print ' Peak',xyz,' <0.5A from another peak'
1810	# return False
1811	# return True
1812	#
1813	def fixSpecialPos(xyz,SGData,Amat):
1814	equivs = G2spc.GenAtom(xyz,SGData,Move=True)
1815	X = []
1816	xyzs = [equiv[0] for equiv in equivs]
1817	for x in xyzs:
1818	if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
1819	X.append(x)
1820	if len(X) > 1:
1821	return np.average(X,axis=0)
1822	else:
1823	return xyz
1824
1825	def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
1826	Mag,x0,y0,z0,sig = parms
1827	z = -((x0-rX)2+(y0-rY)2+(z0-rZ)*2)/(2.sig**2)
1828	# return normMagnp.exp(z)/(sigres*3) #not slower but some faults in LS
1829	return normMag(1.+z+z*2/2.)/(sigres**3)
1830
1831	def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
1832	Mag,x0,y0,z0,sig = parms
1833	M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1834	return M
1835
1836	def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
1837	Mag,x0,y0,z0,sig = parms
1838	dMdv = np.zeros(([5,]+list(rX.shape)))
1839	delt = .01
1840	for i in range(5):
1841	parms[i] -= delt
1842	rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1843	parms[i] += 2.*delt
1844	rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1845	parms[i] -= delt
1846	dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
1847	rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1848	Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
1849	dMdv = np.reshape(dMdv,(5,rX.size))
1850	Hess = np.inner(dMdv,dMdv)
1851
1852	return Vec,Hess
1853
1854	generalData = data['General']
1855	phaseName = generalData['Name']
1856	SGData = generalData['SGData']
1857	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1858	drawingData = data['Drawing']
1859	peaks = []
1860	mags = []
1861	dzeros = []
1862	try:
1863	mapData = generalData['Map']
1864	contLevel = mapData['cutOff']*mapData['rhoMax']/100.
1865	rho = copy.copy(mapData['rho']) #don't mess up original
1866	mapHalf = np.array(rho.shape)/2
1867	res = mapData['Resolution']
1868	incre = np.array(rho.shape,dtype=np.float)
1869	step = max(1.0,1./res)+1
1870	steps = np.array(3*[step,])
1871	except KeyError:
1872	print '**** ERROR - Fourier map not defined'
1873	return peaks,mags
1874	rhoMask = ma.array(rho,mask=(rho<contLevel))
1875	indices = (-1,0,1)
1876	rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1877	for roll in rolls:
1878	if np.any(roll):
1879	rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
1880	indx = np.transpose(rhoMask.nonzero())
1881	peaks = indx/incre
1882	mags = rhoMask[rhoMask.nonzero()]
1883	for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
1884	rho = rollMap(rho,ind)
1885	rMM = mapHalf-steps
1886	rMP = mapHalf+steps+1
1887	rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1888	peakInt = np.sum(rhoPeak)res*3
1889	rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1890	x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0] #magnitude, position & width(sig)
1891	result = HessianLSQ(peakFunc,x0,Hess=peakHess,
1892	args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
1893	x1 = result[0]
1894	if not np.any(x1 < 0):
1895	mag = x1[0]
1896	peak = (np.array(x1[1:4])-ind)/incre
1897	peak = fixSpecialPos(peak,SGData,Amat)
1898	rho = rollMap(rho,-ind)
1899	dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
1900	return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1901
1902	def sortArray(data,pos,reverse=False):
1903	'''data is a list of items
1904	sort by pos in list; reverse if True
1905	'''
1906	T = []
1907	for i,M in enumerate(data):
1908	T.append((M[pos],i))
1909	D = dict(zip(T,data))
1910	T.sort()
1911	if reverse:
1912	T.reverse()
1913	X = []
1914	for key in T:
1915	X.append(D[key])
1916	return X
1917
1918	def PeaksEquiv(data,Ind):
1919	'''Find the equivalent map peaks for those selected. Works on the
1920	contents of data['Map Peaks'].
1921
1922	:param data: the phase data structure
1923	:param list Ind: list of selected peak indices
1924	:returns: augmented list of peaks including those related by symmetry to the
1925	ones in Ind
1926
1927	'''
1928	def Duplicate(xyz,peaks,Amat):
1929	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1930	return True
1931	return False
1932
1933	generalData = data['General']
1934	cell = generalData['Cell'][1:7]
1935	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1936	A = G2lat.cell2A(cell)
1937	SGData = generalData['SGData']
1938	mapPeaks = data['Map Peaks']
1939	XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
1940	Indx = {}
1941	for ind in Ind:
1942	xyz = np.array(mapPeaks[ind][1:4])
1943	xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1944	for jnd,xyz in enumerate(XYZ):
1945	Indx[jnd] = Duplicate(xyz,xyzs,Amat)
1946	Ind = []
1947	for ind in Indx:
1948	if Indx[ind]:
1949	Ind.append(ind)
1950	return Ind
1951
1952	def PeaksUnique(data,Ind):
1953	'''Finds the symmetry unique set of peaks from those selected. Works on the
1954	contents of data['Map Peaks'].
1955
1956	:param data: the phase data structure
1957	:param list Ind: list of selected peak indices
1958	:returns: the list of symmetry unique peaks from among those given in Ind
1959
1960	'''
1961	# XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
1962
1963	def noDuplicate(xyz,peaks,Amat):
1964	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=1.0) for peak in peaks]:
1965	return False
1966	return True
1967
1968	generalData = data['General']
1969	cell = generalData['Cell'][1:7]
1970	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1971	A = G2lat.cell2A(cell)
1972	SGData = generalData['SGData']
1973	mapPeaks = data['Map Peaks']
1974	Indx = {}
1975	XYZ = {}
1976	for ind in Ind:
1977	XYZ[ind] = np.array(mapPeaks[ind][1:4])
1978	Indx[ind] = True
1979	for ind in Ind:
1980	if Indx[ind]:
1981	xyz = XYZ[ind]
1982	for jnd in Ind:
1983	if ind != jnd and Indx[jnd]:
1984	Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
1985	xyzs = np.array([equiv[0] for equiv in Equiv])
1986	Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
1987	Ind = []
1988	for ind in Indx:
1989	if Indx[ind]:
1990	Ind.append(ind)
1991	return Ind
1992
1993	################################################################################
1994	##### single peak fitting profile fxn stuff
1995	################################################################################
1996
1997	def getCWsig(ins,pos):
1998	'''default doc string
1999
2000	:param type name: description
2001
2002	:returns: type name: description
2003
2004	'''
2005	tp = tand(pos/2.0)
2006	return ins['U']tp2+ins['V']tp+ins['W']
2007
2008	def getCWsigDeriv(pos):
2009	'''default doc string
2010
2011	:param type name: description
2012
2013	:returns: type name: description
2014
2015	'''
2016	tp = tand(pos/2.0)
2017	return tp**2,tp,1.0
2018
2019	def getCWgam(ins,pos):
2020	'''default doc string
2021
2022	:param type name: description
2023
2024	:returns: type name: description
2025
2026	'''
2027	return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
2028
2029	def getCWgamDeriv(pos):
2030	'''default doc string
2031
2032	:param type name: description
2033
2034	:returns: type name: description
2035
2036	'''
2037	return 1./cosd(pos/2.0),tand(pos/2.0)
2038
2039	def getTOFsig(ins,dsp):
2040	'''default doc string
2041
2042	:param type name: description
2043
2044	:returns: type name: description
2045
2046	'''
2047	return ins['sig-0']+ins['sig-1']dsp2+ins['sig-q']dsp
2048
2049	def getTOFsigDeriv(dsp):
2050	'''default doc string
2051
2052	:param type name: description
2053
2054	:returns: type name: description
2055
2056	'''
2057	return 1.0,dsp**2,dsp
2058
2059	def getTOFgamma(ins,dsp):
2060	'''default doc string
2061
2062	:param type name: description
2063
2064	:returns: type name: description
2065
2066	'''
2067	return ins['X']dsp+ins['Y']dsp**2
2068
2069	def getTOFgammaDeriv(dsp):
2070	'''default doc string
2071
2072	:param type name: description
2073
2074	:returns: type name: description
2075
2076	'''
2077	return dsp,dsp**2
2078
2079	def getTOFbeta(ins,dsp):
2080	'''default doc string
2081
2082	:param type name: description
2083
2084	:returns: type name: description
2085
2086	'''
2087	return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp
2088
2089	def getTOFbetaDeriv(dsp):
2090	'''default doc string
2091
2092	:param type name: description
2093
2094	:returns: type name: description
2095
2096	'''
2097	return 1.0,1./dsp**4,1./dsp
2098
2099	def getTOFalpha(ins,dsp):
2100	'''default doc string
2101
2102	:param type name: description
2103
2104	:returns: type name: description
2105
2106	'''
2107	return ins['alpha']/dsp
2108
2109	def getTOFalphaDeriv(dsp):
2110	'''default doc string
2111
2112	:param type name: description
2113
2114	:returns: type name: description
2115
2116	'''
2117	return 1./dsp
2118
2119	def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
2120	'''default doc string
2121
2122	:param type name: description
2123
2124	:returns: type name: description
2125
2126	'''
2127	ind = 0
2128	if useFit:
2129	ind = 1
2130	ins = {}
2131	if 'C' in Parms['Type'][0]: #CW data - TOF later in an elif
2132	for x in ['U','V','W','X','Y']:
2133	ins[x] = Parms[x][ind]
2134	if ifQ: #qplot - convert back to 2-theta
2135	pos = 2.0asind(poswave/(4*math.pi))
2136	sig = getCWsig(ins,pos)
2137	gam = getCWgam(ins,pos)
2138	XY = [pos,0, mag,1, sig,0, gam,0] #default refine intensity 1st
2139	else:
2140	if ifQ:
2141	dsp = 2.*np.pi/pos
2142	pos = Parms['difC']*dsp
2143	else:
2144	dsp = pos/Parms['difC'][1]
2145	if 'Pdabc' in Parms2:
2146	for x in ['sig-0','sig-1','sig-q','X','Y']:
2147	ins[x] = Parms[x][ind]
2148	Pdabc = Parms2['Pdabc'].T
2149	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2150	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2151	else:
2152	for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','X','Y']:
2153	ins[x] = Parms[x][ind]
2154	alp = getTOFalpha(ins,dsp)
2155	bet = getTOFbeta(ins,dsp)
2156	sig = getTOFsig(ins,dsp)
2157	gam = getTOFgamma(ins,dsp)
2158	XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2159	return XY
2160
2161	################################################################################
2162	##### MC/SA stuff
2163	################################################################################
2164
2165	#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2166	# Original Author: Travis Oliphant 2002
2167	# Bug-fixes in 2006 by Tim Leslie
2168
2169
2170	import numpy
2171	from numpy import asarray, tan, exp, ones, squeeze, sign, \
2172	all, log, sqrt, pi, shape, array, minimum, where
2173	from numpy import random
2174
2175	#__all__ = ['anneal']
2176
2177	_double_min = numpy.finfo(float).min
2178	_double_max = numpy.finfo(float).max
2179	class base_schedule(object):
2180	def __init__(self):
2181	self.dwell = 20
2182	self.learn_rate = 0.5
2183	self.lower = -10
2184	self.upper = 10
2185	self.Ninit = 50
2186	self.accepted = 0
2187	self.tests = 0
2188	self.feval = 0
2189	self.k = 0
2190	self.T = None
2191
2192	def init(self, **options):
2193	self.__dict__.update(options)
2194	self.lower = asarray(self.lower)
2195	self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2196	self.upper = asarray(self.upper)
2197	self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2198	self.k = 0
2199	self.accepted = 0
2200	self.feval = 0
2201	self.tests = 0
2202
2203	def getstart_temp(self, best_state):
2204	""" Find a matching starting temperature and starting parameters vector
2205	i.e. find x0 such that func(x0) = T0.
2206
2207	Parameters
2208	----------
2209	best_state : _state
2210	A _state object to store the function value and x0 found.
2211
2212	returns
2213	-------
2214	x0 : array
2215	The starting parameters vector.
2216	"""
2217
2218	assert(not self.dims is None)
2219	lrange = self.lower
2220	urange = self.upper
2221	fmax = _double_min
2222	fmin = _double_max
2223	for _ in range(self.Ninit):
2224	x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2225	fval = self.func(x0, *self.args)
2226	self.feval += 1
2227	if fval > fmax:
2228	fmax = fval
2229	if fval < fmin:
2230	fmin = fval
2231	best_state.cost = fval
2232	best_state.x = array(x0)
2233
2234	self.T0 = (fmax-fmin)*1.5
2235	return best_state.x
2236
2237	def set_range(self,x0,frac):
2238	delrange = frac*(self.upper-self.lower)
2239	self.upper = x0+delrange
2240	self.lower = x0-delrange
2241
2242	def accept_test(self, dE):
2243	T = self.T
2244	self.tests += 1
2245	if dE < 0:
2246	self.accepted += 1
2247	return 1
2248	p = exp(-dE*1.0/self.boltzmann/T)
2249	if (p > random.uniform(0.0, 1.0)):
2250	self.accepted += 1
2251	return 1
2252	return 0
2253
2254	def update_guess(self, x0):
2255	return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2256
2257	def update_temp(self, x0):
2258	pass
2259
2260
2261	# A schedule due to Lester Ingber modified to use bounds - OK
2262	class fast_sa(base_schedule):
2263	def init(self, **options):
2264	self.__dict__.update(options)
2265	if self.m is None:
2266	self.m = 1.0
2267	if self.n is None:
2268	self.n = 1.0
2269	self.c = self.m * exp(-self.n * self.quench)
2270
2271	def update_guess(self, x0):
2272	x0 = asarray(x0)
2273	u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2274	T = self.T
2275	xc = (sign(u-0.5)T((1+1.0/T)*abs(2u-1)-1.0)+1.0)/2.0
2276	xnew = xc*(self.upper - self.lower)+self.lower
2277	return xnew
2278	# y = sign(u-0.5)T((1+1.0/T)*abs(2u-1)-1.0)
2279	# xc = y*(self.upper - self.lower)
2280	# xnew = x0 + xc
2281	# return xnew
2282
2283	def update_temp(self):
2284	self.T = self.T0exp(-self.c self.k**(self.quench))
2285	self.k += 1
2286	return
2287
2288	class cauchy_sa(base_schedule): #modified to use bounds - not good
2289	def update_guess(self, x0):
2290	x0 = asarray(x0)
2291	numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2292	xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2293	xnew = xc*(self.upper - self.lower)+self.lower
2294	return xnew
2295	# numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2296	# xc = self.learn_rate * self.T * tan(numbers)
2297	# xnew = x0 + xc
2298	# return xnew
2299
2300	def update_temp(self):
2301	self.T = self.T0/(1+self.k)
2302	self.k += 1
2303	return
2304
2305	class boltzmann_sa(base_schedule):
2306	# def update_guess(self, x0):
2307	# std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2308	# x0 = asarray(x0)
2309	# xc = squeeze(random.normal(0, 1.0, size=self.dims))
2310	#
2311	# xnew = x0 + xcstdself.learn_rate
2312	# return xnew
2313
2314	def update_temp(self):
2315	self.k += 1
2316	self.T = self.T0 / log(self.k+1.0)
2317	return
2318
2319	class log_sa(base_schedule): #OK
2320
2321	def init(self,**options):
2322	self.__dict__.update(options)
2323
2324	def update_guess(self,x0): #same as default
2325	return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2326
2327	def update_temp(self):
2328	self.k += 1
2329	self.T = self.T0self.slope*self.k
2330
2331	class _state(object):
2332	def __init__(self):
2333	self.x = None
2334	self.cost = None
2335
2336	# TODO:
2337	# allow for general annealing temperature profile
2338	# in that case use update given by alpha and omega and
2339	# variation of all previous updates and temperature?
2340
2341	# Simulated annealing
2342
2343	def anneal(func, x0, args=(), schedule='fast', full_output=0,
2344	T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2345	boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2346	lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
2347	ranRange=0.10,autoRan=False,dlg=None):
2348	"""Minimize a function using simulated annealing.
2349
2350	Schedule is a schedule class implementing the annealing schedule.
2351	Available ones are 'fast', 'cauchy', 'boltzmann'
2352
2353	:param callable func: f(x, \*args)
2354	Function to be optimized.
2355	:param ndarray x0:
2356	Initial guess.
2357	:param tuple args:
2358	Extra parameters to `func`.
2359	:param base_schedule schedule:
2360	Annealing schedule to use (a class).
2361	:param bool full_output:
2362	Whether to return optional outputs.
2363	:param float T0:
2364	Initial Temperature (estimated as 1.2 times the largest
2365	cost-function deviation over random points in the range).
2366	:param float Tf:
2367	Final goal temperature.
2368	:param int maxeval:
2369	Maximum function evaluations.
2370	:param int maxaccept:
2371	Maximum changes to accept.
2372	:param int maxiter:
2373	Maximum cooling iterations.
2374	:param float learn_rate:
2375	Scale constant for adjusting guesses.
2376	:param float boltzmann:
2377	Boltzmann constant in acceptance test
2378	(increase for less stringent test at each temperature).
2379	:param float feps:
2380	Stopping relative error tolerance for the function value in
2381	last four coolings.
2382	:param float quench,m,n:
2383	Parameters to alter fast_sa schedule.
2384	:param float/ndarray lower,upper:
2385	Lower and upper bounds on `x`.
2386	:param int dwell:
2387	The number of times to search the space at each temperature.
2388	:param float slope:
2389	Parameter for log schedule
2390	:param bool ranStart=False:
2391	True for set 10% of ranges about x
2392
2393	:returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2394
2395	* xmin (ndarray): Point giving smallest value found.
2396	* Jmin (float): Minimum value of function found.
2397	* T (float): Final temperature.
2398	* feval (int): Number of function evaluations.
2399	* iters (int): Number of cooling iterations.
2400	* accept (int): Number of tests accepted.
2401	* retval (int): Flag indicating stopping condition:
2402
2403	* 0: Points no longer changing
2404	* 1: Cooled to final temperature
2405	* 2: Maximum function evaluations
2406	* 3: Maximum cooling iterations reached
2407	* 4: Maximum accepted query locations reached
2408	* 5: Final point not the minimum amongst encountered points
2409
2410	Notes:
2411	Simulated annealing is a random algorithm which uses no derivative
2412	information from the function being optimized. In practice it has
2413	been more useful in discrete optimization than continuous
2414	optimization, as there are usually better algorithms for continuous
2415	optimization problems.
2416
2417	Some experimentation by trying the difference temperature
2418	schedules and altering their parameters is likely required to
2419	obtain good performance.
2420
2421	The randomness in the algorithm comes from random sampling in numpy.
2422	To obtain the same results you can call numpy.random.seed with the
2423	same seed immediately before calling scipy.optimize.anneal.
2424
2425	We give a brief description of how the three temperature schedules
2426	generate new points and vary their temperature. Temperatures are
2427	only updated with iterations in the outer loop. The inner loop is
2428	over xrange(dwell), and new points are generated for
2429	every iteration in the inner loop. (Though whether the proposed
2430	new points are accepted is probabilistic.)
2431
2432	For readability, let d denote the dimension of the inputs to func.
2433	Also, let x_old denote the previous state, and k denote the
2434	iteration number of the outer loop. All other variables not
2435	defined below are input variables to scipy.optimize.anneal itself.
2436
2437	In the 'fast' schedule the updates are ::
2438
2439	u ~ Uniform(0, 1, size=d)
2440	y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2441	xc = y * (upper - lower)
2442	x_new = x_old + xc
2443
2444	c = n * exp(-n * quench)
2445	T_new = T0 * exp(-c * k**quench)
2446
2447
2448	In the 'cauchy' schedule the updates are ::
2449
2450	u ~ Uniform(-pi/2, pi/2, size=d)
2451	xc = learn_rate * T * tan(u)
2452	x_new = x_old + xc
2453
2454	T_new = T0 / (1+k)
2455
2456	In the 'boltzmann' schedule the updates are ::
2457
2458	std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2459	y ~ Normal(0, std, size=d)
2460	x_new = x_old + learn_rate * y
2461
2462	T_new = T0 / log(1+k)
2463
2464	"""
2465	x0 = asarray(x0)
2466	lower = asarray(lower)
2467	upper = asarray(upper)
2468
2469	schedule = eval(schedule+'_sa()')
2470	# initialize the schedule
2471	schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2472	learn_rate=learn_rate, lower=lower, upper=upper,
2473	m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2474
2475	current_state, last_state, best_state = _state(), _state(), _state()
2476	if ranStart:
2477	schedule.set_range(x0,ranRange)
2478	if T0 is None:
2479	x0 = schedule.getstart_temp(best_state)
2480	else:
2481	x0 = random.uniform(size=len(x0))*(upper-lower) + lower
2482	best_state.x = None
2483	best_state.cost = numpy.Inf
2484
2485	last_state.x = asarray(x0).copy()
2486	fval = func(x0,*args)
2487	schedule.feval += 1
2488	last_state.cost = fval
2489	if last_state.cost < best_state.cost:
2490	best_state.cost = fval
2491	best_state.x = asarray(x0).copy()
2492	schedule.T = schedule.T0
2493	fqueue = [100, 300, 500, 700]
2494	iters = 1
2495	keepGoing = True
2496	bestn = 0
2497	while keepGoing:
2498	retval = 0
2499	for n in xrange(dwell):
2500	current_state.x = schedule.update_guess(last_state.x)
2501	current_state.cost = func(current_state.x,*args)
2502	schedule.feval += 1
2503
2504	dE = current_state.cost - last_state.cost
2505	if schedule.accept_test(dE):
2506	last_state.x = current_state.x.copy()
2507	last_state.cost = current_state.cost
2508	if last_state.cost < best_state.cost:
2509	best_state.x = last_state.x.copy()
2510	best_state.cost = last_state.cost
2511	bestn = n
2512	if best_state.cost < 1.0 and autoRan:
2513	schedule.set_range(x0,best_state.cost/2.)
2514	if dlg:
2515	GoOn = dlg.Update(min(100.,best_state.cost*100),
2516	newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
2517	'Best trial:',bestn, \
2518	'MC/SA Residual:',best_state.cost*100,'%', \
2519	))[0]
2520	if not GoOn:
2521	best_state.x = last_state.x.copy()
2522	best_state.cost = last_state.cost
2523	retval = 5
2524	schedule.update_temp()
2525	iters += 1
2526	# Stopping conditions
2527	# 0) last saved values of f from each cooling step
2528	# are all very similar (effectively cooled)
2529	# 1) Tf is set and we are below it
2530	# 2) maxeval is set and we are past it
2531	# 3) maxiter is set and we are past it
2532	# 4) maxaccept is set and we are past it
2533	# 5) user canceled run via progress bar
2534
2535	fqueue.append(squeeze(last_state.cost))
2536	fqueue.pop(0)
2537	af = asarray(fqueue)*1.0
2538	if retval == 5:
2539	print ' User terminated run; incomplete MC/SA'
2540	keepGoing = False
2541	break
2542	if all(abs((af-af[0])/af[0]) < feps):
2543	retval = 0
2544	if abs(af[-1]-best_state.cost) > feps*10:
2545	retval = 5
2546	# print "Warning: Cooled to %f at %s but this is not" \
2547	# % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2548	# + " the smallest point found."
2549	break
2550	if (Tf is not None) and (schedule.T < Tf):
2551	retval = 1
2552	break
2553	if (maxeval is not None) and (schedule.feval > maxeval):
2554	retval = 2
2555	break
2556	if (iters > maxiter):
2557	print "Warning: Maximum number of iterations exceeded."
2558	retval = 3
2559	break
2560	if (maxaccept is not None) and (schedule.accepted > maxaccept):
2561	retval = 4
2562	break
2563
2564	if full_output:
2565	return best_state.x, best_state.cost, schedule.T, \
2566	schedule.feval, iters, schedule.accepted, retval
2567	else:
2568	return best_state.x, retval
2569
2570	def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
2571	outlist = []
2572	random.seed(int(time.time())%100000+nprocess) #make sure each process has a different random start
2573	for n in range(iCyc):
2574	result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
2575	outlist.append(result[0])
2576	print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
2577	out_q.put(outlist)
2578
2579	def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
2580	import multiprocessing as mp
2581
2582	nprocs = mp.cpu_count()
2583	out_q = mp.Queue()
2584	procs = []
2585	iCyc = np.zeros(nprocs)
2586	for i in range(nCyc):
2587	iCyc[i%nprocs] += 1
2588	for i in range(nprocs):
2589	p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
2590	procs.append(p)
2591	p.start()
2592	resultlist = []
2593	for i in range(nprocs):
2594	resultlist += out_q.get()
2595	for p in procs:
2596	p.join()
2597	return resultlist
2598
2599	def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2600	'''default doc string
2601
2602	:param type name: description
2603
2604	:returns: type name: description
2605
2606	'''
2607
2608	twopi = 2.0*np.pi
2609	global tsum
2610	tsum = 0.
2611
2612	def getMDparms(item,pfx,parmDict,varyList):
2613	parmDict[pfx+'MDaxis'] = item['axis']
2614	parmDict[pfx+'MDval'] = item['Coef'][0]
2615	if item['Coef'][1]:
2616	varyList += [pfx+'MDval',]
2617	limits = item['Coef'][2]
2618	lower.append(limits[0])
2619	upper.append(limits[1])
2620
2621	def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2622	parmDict[pfx+'Atype'] = item['atType']
2623	aTypes \|= set([item['atType'],])
2624	pstr = ['Ax','Ay','Az']
2625	XYZ = [0,0,0]
2626	for i in range(3):
2627	name = pfx+pstr[i]
2628	parmDict[name] = item['Pos'][0][i]
2629	XYZ[i] = parmDict[name]
2630	if item['Pos'][1][i]:
2631	varyList += [name,]
2632	limits = item['Pos'][2][i]
2633	lower.append(limits[0])
2634	upper.append(limits[1])
2635	parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2636
2637	def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2638	parmDict[mfx+'MolCent'] = item['MolCent']
2639	parmDict[mfx+'RBId'] = item['RBId']
2640	pstr = ['Px','Py','Pz']
2641	ostr = ['Qa','Qi','Qj','Qk'] #angle,vector not quaternion
2642	for i in range(3):
2643	name = mfx+pstr[i]
2644	parmDict[name] = item['Pos'][0][i]
2645	if item['Pos'][1][i]:
2646	varyList += [name,]
2647	limits = item['Pos'][2][i]
2648	lower.append(limits[0])
2649	upper.append(limits[1])
2650	AV = item['Ori'][0]
2651	A = AV[0]
2652	V = AV[1:]
2653	for i in range(4):
2654	name = mfx+ostr[i]
2655	if i:
2656	parmDict[name] = V[i-1]
2657	else:
2658	parmDict[name] = A
2659	if item['Ovar'] == 'AV':
2660	varyList += [name,]
2661	limits = item['Ori'][2][i]
2662	lower.append(limits[0])
2663	upper.append(limits[1])
2664	elif item['Ovar'] == 'A' and not i:
2665	varyList += [name,]
2666	limits = item['Ori'][2][i]
2667	lower.append(limits[0])
2668	upper.append(limits[1])
2669	if 'Tor' in item: #'Tor' not there for 'Vector' RBs
2670	for i in range(len(item['Tor'][0])):
2671	name = mfx+'Tor'+str(i)
2672	parmDict[name] = item['Tor'][0][i]
2673	if item['Tor'][1][i]:
2674	varyList += [name,]
2675	limits = item['Tor'][2][i]
2676	lower.append(limits[0])
2677	upper.append(limits[1])
2678	atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2679	aTypes \|= set(atypes)
2680	atNo += len(atypes)
2681	return atNo
2682
2683	def GetAtomM(Xdata,SGData):
2684	Mdata = []
2685	for xyz in Xdata:
2686	Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2687	return np.array(Mdata)
2688
2689	def GetAtomT(RBdata,parmDict):
2690	'Needs a doc string'
2691	atNo = parmDict['atNo']
2692	nfixAt = parmDict['nfixAt']
2693	Tdata = atNo*[' ',]
2694	for iatm in range(nfixAt):
2695	parm = ':'+str(iatm)+':Atype'
2696	if parm in parmDict:
2697	Tdata[iatm] = aTypes.index(parmDict[parm])
2698	iatm = nfixAt
2699	for iObj in range(parmDict['nObj']):
2700	pfx = str(iObj)+':'
2701	if parmDict[pfx+'Type'] in ['Vector','Residue']:
2702	if parmDict[pfx+'Type'] == 'Vector':
2703	RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2704	nAtm = len(RBRes['rbVect'][0])
2705	else: #Residue
2706	RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2707	nAtm = len(RBRes['rbXYZ'])
2708	for i in range(nAtm):
2709	Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
2710	iatm += 1
2711	elif parmDict[pfx+'Type'] == 'Atom':
2712	atNo = parmDict[pfx+'atNo']
2713	parm = pfx+'Atype' #remove extra ':'
2714	if parm in parmDict:
2715	Tdata[atNo] = aTypes.index(parmDict[parm])
2716	iatm += 1
2717	else:
2718	continue #skips March Dollase
2719	return Tdata
2720
2721	def GetAtomX(RBdata,parmDict):
2722	'Needs a doc string'
2723	Bmat = parmDict['Bmat']
2724	atNo = parmDict['atNo']
2725	nfixAt = parmDict['nfixAt']
2726	Xdata = np.zeros((3,atNo))
2727	keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
2728	for iatm in range(nfixAt):
2729	for key in keys:
2730	parm = ':'+str(iatm)+key
2731	if parm in parmDict:
2732	keys[key][iatm] = parmDict[parm]
2733	iatm = nfixAt
2734	for iObj in range(parmDict['nObj']):
2735	pfx = str(iObj)+':'
2736	if parmDict[pfx+'Type'] in ['Vector','Residue']:
2737	if parmDict[pfx+'Type'] == 'Vector':
2738	RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2739	vecs = RBRes['rbVect']
2740	mags = RBRes['VectMag']
2741	Cart = np.zeros_like(vecs[0])
2742	for vec,mag in zip(vecs,mags):
2743	Cart += vec*mag
2744	elif parmDict[pfx+'Type'] == 'Residue':
2745	RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2746	Cart = np.array(RBRes['rbXYZ'])
2747	for itor,seq in enumerate(RBRes['rbSeq']):
2748	QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
2749	Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
2750	if parmDict[pfx+'MolCent'][1]:
2751	Cart -= parmDict[pfx+'MolCent'][0]
2752	Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2753	Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
2754	Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
2755	iatm += len(Cart)
2756	elif parmDict[pfx+'Type'] == 'Atom':
2757	atNo = parmDict[pfx+'atNo']
2758	for key in keys:
2759	parm = pfx+key[1:] #remove extra ':'
2760	if parm in parmDict:
2761	keys[key][atNo] = parmDict[parm]
2762	iatm += 1
2763	else:
2764	continue #skips March Dollase
2765	return Xdata.T
2766
2767	def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
2768	allX = np.inner(Xdata,SGM)+SGT
2769	allT = np.repeat(Tdata,allX.shape[1])
2770	allM = np.repeat(Mdata,allX.shape[1])
2771	allX = np.reshape(allX,(-1,3))
2772	return allT,allM,allX
2773
2774	def getAllX(Xdata,SGM,SGT):
2775	allX = np.inner(Xdata,SGM)+SGT
2776	allX = np.reshape(allX,(-1,3))
2777	return allX
2778
2779	def normQuaternions(RBdata,parmDict,varyList,values):
2780	for iObj in range(parmDict['nObj']):
2781	pfx = str(iObj)+':'
2782	if parmDict[pfx+'Type'] in ['Vector','Residue']:
2783	Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2784	A,V = Q2AVdeg(Qori)
2785	for i,name in enumerate(['Qa','Qi','Qj','Qk']):
2786	if i:
2787	parmDict[pfx+name] = V[i-1]
2788	else:
2789	parmDict[pfx+name] = A
2790
2791	def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
2792	''' Compute structure factors for all h,k,l for phase
2793	input:
2794	refList: [ref] where each ref = h,k,l,m,d,...
2795	rcov: array[nref,nref] covariance terms between Fo^2 values
2796	ifInv: bool True if centrosymmetric
2797	allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
2798	RBdata: [dict] rigid body dictionary
2799	varyList: [list] names of varied parameters in MC/SA (not used here)
2800	ParmDict: [dict] problem parameters
2801	puts result F^2 in each ref[5] in refList
2802	returns:
2803	delt-Frcovdelt-F/sum(Fo^2)^2
2804	'''
2805	global tsum
2806	t0 = time.time()
2807	parmDict.update(dict(zip(varyList,values))) #update parameter tables
2808	Xdata = GetAtomX(RBdata,parmDict) #get new atom coords from RB
2809	allX = getAllX(Xdata,SGM,SGT) #fill unit cell - dups. OK
2810	MDval = parmDict['0:MDval'] #get March-Dollase coeff
2811	HX2pi = 2.np.pinp.inner(allX,refList[:3].T) #form 2piHX for every H,X pair
2812	Aterm = refList[3]np.sum(allFFnp.cos(HX2pi),axis=0)**2 #compute real part for all H
2813	refList[5] = Aterm
2814	if not ifInv:
2815	refList[5] += refList[3]np.sum(allFFnp.sin(HX2pi),axis=0)**2 #imaginary part for all H
2816	if len(cosTable): #apply MD correction
2817	refList[5] = np.sum(np.sqrt((MDval/(cosTable(MDval3-1.)+1.))3),axis=1)/cosTable.shape[1]
2818	sumFcsq = np.sum(refList[5])
2819	scale = parmDict['sumFosq']/sumFcsq
2820	refList[5] *= scale
2821	refList[6] = refList[4]-refList[5]
2822	M = np.inner(refList[6],np.inner(rcov,refList[6]))
2823	tsum += (time.time()-t0)
2824	return M/np.sum(refList[4]**2)
2825
2826	sq8ln2 = np.sqrt(8*np.log(2))
2827	sq2pi = np.sqrt(2*np.pi)
2828	sq4pi = np.sqrt(4*np.pi)
2829	generalData = data['General']
2830	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2831	Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
2832	SGData = generalData['SGData']
2833	SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
2834	SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
2835	SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
2836	fixAtoms = data['Atoms'] #if any
2837	cx,ct,cs = generalData['AtomPtrs'][:3]
2838	aTypes = set([])
2839	parmDict = {'Bmat':Bmat,'Gmat':Gmat}
2840	varyList = []
2841	atNo = 0
2842	for atm in fixAtoms:
2843	pfx = ':'+str(atNo)+':'
2844	parmDict[pfx+'Atype'] = atm[ct]
2845	aTypes \|= set([atm[ct],])
2846	pstr = ['Ax','Ay','Az']
2847	parmDict[pfx+'Amul'] = atm[cs+1]
2848	for i in range(3):
2849	name = pfx+pstr[i]
2850	parmDict[name] = atm[cx+i]
2851	atNo += 1
2852	parmDict['nfixAt'] = len(fixAtoms)
2853	MCSA = generalData['MCSA controls']
2854	reflName = MCSA['Data source']
2855	phaseName = generalData['Name']
2856	MCSAObjs = data['MCSA']['Models'] #list of MCSA models
2857	upper = []
2858	lower = []
2859	MDvec = np.zeros(3)
2860	for i,item in enumerate(MCSAObjs):
2861	mfx = str(i)+':'
2862	parmDict[mfx+'Type'] = item['Type']
2863	if item['Type'] == 'MD':
2864	getMDparms(item,mfx,parmDict,varyList)
2865	MDvec = np.array(item['axis'])
2866	elif item['Type'] == 'Atom':
2867	getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
2868	parmDict[mfx+'atNo'] = atNo
2869	atNo += 1
2870	elif item['Type'] in ['Residue','Vector']:
2871	atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
2872	parmDict['atNo'] = atNo #total no. of atoms
2873	parmDict['nObj'] = len(MCSAObjs)
2874	aTypes = list(aTypes)
2875	Tdata = GetAtomT(RBdata,parmDict)
2876	Xdata = GetAtomX(RBdata,parmDict)
2877	Mdata = GetAtomM(Xdata,SGData)
2878	allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
2879	FFtables = G2el.GetFFtable(aTypes)
2880	refs = []
2881	allFF = []
2882	cosTable = []
2883	sumFosq = 0
2884	if 'PWDR' in reflName:
2885	for ref in reflData:
2886	h,k,l,m,d,pos,sig,gam,f = ref[:9]
2887	if d >= MCSA['dmin']:
2888	sig = G2pwd.getgamFW(sig,gam)/sq8ln2 #--> sig from FWHM
2889	SQ = 0.25/d**2
2890	allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2891	refs.append([h,k,l,m,f*m,pos,sig])
2892	sumFosq += f*m
2893	Heqv = np.inner(np.array([h,k,l]),SGMT)
2894	cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
2895	nRef = len(refs)
2896	cosTable = np.array(cosTable)**2
2897	rcov = np.zeros((nRef,nRef))
2898	for iref,refI in enumerate(refs):
2899	rcov[iref][iref] = 1./(sq4pi*refI[6])
2900	for jref,refJ in enumerate(refs[:iref]):
2901	t1 = refI[6]2+refJ[6]2
2902	t2 = (refJ[5]-refI[5])*2/(2.t1)
2903	if t2 > 10.:
2904	rcov[iref][jref] = 0.
2905	else:
2906	rcov[iref][jref] = 1./(sq2pinp.sqrt(t1)np.exp(t2))
2907	rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2908	Rdiag = np.sqrt(np.diag(rcov))
2909	Rnorm = np.outer(Rdiag,Rdiag)
2910	rcov /= Rnorm
2911	elif 'Pawley' in reflName: #need a bail out if Pawley cov matrix doesn't exist.
2912	vNames = []
2913	pfx = str(data['pId'])+'::PWLref:'
2914	for iref,refI in enumerate(reflData): #Pawley reflection set
2915	h,k,l,m,d,v,f,s = refI
2916	if d >= MCSA['dmin'] and v: #skip unrefined ones
2917	vNames.append(pfx+str(iref))
2918	SQ = 0.25/d**2
2919	allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2920	refs.append([h,k,l,m,f*m,iref,0.])
2921	sumFosq += f*m
2922	Heqv = np.inner(np.array([h,k,l]),SGMT)
2923	cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
2924	cosTable = np.array(cosTable)**2
2925	nRef = len(refs)
2926	if covData['freshCOV'] and generalData['doPawley'] and MCSA.get('newDmin',True):
2927	vList = covData['varyList']
2928	covMatrix = covData['covMatrix']
2929	rcov = getVCov(vNames,vList,covMatrix)
2930	rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2931	Rdiag = np.sqrt(np.diag(rcov))
2932	Rdiag = np.where(Rdiag,Rdiag,1.0)
2933	Rnorm = np.outer(Rdiag,Rdiag)
2934	rcov /= Rnorm
2935	MCSA['rcov'] = rcov
2936	covData['freshCOV'] = False
2937	MCSA['newDmin'] = False
2938	else:
2939	rcov = MCSA['rcov']
2940	elif 'HKLF' in reflName:
2941	for ref in reflData:
2942	[h,k,l,m,d],f = ref[:5],ref[6]
2943	if d >= MCSA['dmin']:
2944	SQ = 0.25/d**2
2945	allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2946	refs.append([h,k,l,m,f*m,0.,0.])
2947	sumFosq += f*m
2948	nRef = len(refs)
2949	rcov = np.identity(len(refs))
2950	allFF = np.array(allFF).T
2951	refs = np.array(refs).T
2952	print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
2953	print ' Number of parameters varied: %d'%(len(varyList))
2954	parmDict['sumFosq'] = sumFosq
2955	x0 = [parmDict[val] for val in varyList]
2956	ifInv = SGData['SGInv']
2957	results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
2958	schedule=MCSA['Algorithm'], full_output=True,
2959	T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
2960	boltzmann=MCSA['boltzmann'], learn_rate=0.5,
2961	quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
2962	lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
2963	ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
2964	M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
2965	# for ref in refs.T:
2966	# print ' %4d %4d %4d %10.3f %10.3f %10.3f'%(int(ref[0]),int(ref[1]),int(ref[2]),ref[4],ref[5],ref[6])
2967	# print np.sqrt((np.sum(refs[6]2)/np.sum(refs[4]2)))
2968	Result = [False,False,results[1],results[2],]+list(results[0])
2969	Result.append(varyList)
2970	return Result,tsum
2971
2972
2973	################################################################################
2974	##### Quaternion stuff
2975	################################################################################
2976
2977	def prodQQ(QA,QB):
2978	''' Grassman quaternion product
2979	QA,QB quaternions; q=r+ai+bj+ck
2980	'''
2981	D = np.zeros(4)
2982	D[0] = QA[0]QB[0]-QA[1]QB[1]-QA[2]QB[2]-QA[3]QB[3]
2983	D[1] = QA[0]QB[1]+QA[1]QB[0]+QA[2]QB[3]-QA[3]QB[2]
2984	D[2] = QA[0]QB[2]-QA[1]QB[3]+QA[2]QB[0]+QA[3]QB[1]
2985	D[3] = QA[0]QB[3]+QA[1]QB[2]-QA[2]QB[1]+QA[3]QB[0]
2986
2987	# D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
2988	# D[1:] = QA[0]QB[1:]+QB[0]QA[1:]+np.cross(QA[1:],QB[1:])
2989
2990	return D
2991
2992	def normQ(QA):
2993	''' get length of quaternion & normalize it
2994	q=r+ai+bj+ck
2995	'''
2996	n = np.sqrt(np.sum(np.array(QA)**2))
2997	return QA/n
2998
2999	def invQ(Q):
3000	'''
3001	get inverse of quaternion
3002	q=r+ai+bj+ck; q* = r-ai-bj-ck
3003	'''
3004	return Q*np.array([1,-1,-1,-1])
3005
3006	def prodQVQ(Q,V):
3007	"""
3008	compute the quaternion vector rotation qvq-1 = v'
3009	q=r+ai+bj+ck
3010	"""
3011	T2 = Q[0]*Q[1]
3012	T3 = Q[0]*Q[2]
3013	T4 = Q[0]*Q[3]
3014	T5 = -Q[1]*Q[1]
3015	T6 = Q[1]*Q[2]
3016	T7 = Q[1]*Q[3]
3017	T8 = -Q[2]*Q[2]
3018	T9 = Q[2]*Q[3]
3019	T10 = -Q[3]*Q[3]
3020	M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
3021	VP = 2.*np.inner(V,M)
3022	return VP+V
3023
3024	def Q2Mat(Q):
3025	''' make rotation matrix from quaternion
3026	q=r+ai+bj+ck
3027	'''
3028	QN = normQ(Q)
3029	aa = QN[0]**2
3030	ab = QN[0]*QN[1]
3031	ac = QN[0]*QN[2]
3032	ad = QN[0]*QN[3]
3033	bb = QN[1]**2
3034	bc = QN[1]*QN[2]
3035	bd = QN[1]*QN[3]
3036	cc = QN[2]**2
3037	cd = QN[2]*QN[3]
3038	dd = QN[3]**2
3039	M = [[aa+bb-cc-dd, 2.(bc-ad), 2.(ac+bd)],
3040	[2(ad+bc), aa-bb+cc-dd, 2.(cd-ab)],
3041	[2(bd-ac), 2.(ab+cd), aa-bb-cc+dd]]
3042	return np.array(M)
3043
3044	def AV2Q(A,V):
3045	''' convert angle (radians) & vector to quaternion
3046	q=r+ai+bj+ck
3047	'''
3048	Q = np.zeros(4)
3049	d = nl.norm(np.array(V))
3050	if d:
3051	V /= d
3052	if not A: #==0.
3053	A = 2.*np.pi
3054	p = A/2.
3055	Q[0] = np.cos(p)
3056	Q[1:4] = V*np.sin(p)
3057	else:
3058	Q[3] = 1.
3059	return Q
3060
3061	def AVdeg2Q(A,V):
3062	''' convert angle (degrees) & vector to quaternion
3063	q=r+ai+bj+ck
3064	'''
3065	Q = np.zeros(4)
3066	d = nl.norm(np.array(V))
3067	if not A: #== 0.!
3068	A = 360.
3069	if d:
3070	V /= d
3071	p = A/2.
3072	Q[0] = cosd(p)
3073	Q[1:4] = V*sind(p)
3074	else:
3075	Q[3] = 1.
3076	return Q
3077
3078	def Q2AVdeg(Q):
3079	''' convert quaternion to angle (degrees 0-360) & normalized vector
3080	q=r+ai+bj+ck
3081	'''
3082	A = 2.*acosd(Q[0])
3083	V = np.array(Q[1:])
3084	V /= sind(A/2.)
3085	return A,V
3086
3087	def Q2AV(Q):
3088	''' convert quaternion to angle (radians 0-2pi) & normalized vector
3089	q=r+ai+bj+ck
3090	'''
3091	A = 2.*np.arccos(Q[0])
3092	V = np.array(Q[1:])
3093	V /= np.sin(A/2.)
3094	return A,V
3095
3096	def randomQ(r0,r1,r2,r3):
3097	''' create random quaternion from 4 random numbers in range (-1,1)
3098	'''
3099	sum = 0
3100	Q = np.array(4)
3101	Q[0] = r0
3102	sum += Q[0]**2
3103	Q[1] = np.sqrt(1.-sum)*r1
3104	sum += Q[1]**2
3105	Q[2] = np.sqrt(1.-sum)*r2
3106	sum += Q[2]**2
3107	Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
3108	return Q
3109
3110	def randomAVdeg(r0,r1,r2,r3):
3111	''' create random angle (deg),vector from 4 random number in range (-1,1)
3112	'''
3113	return Q2AVdeg(randomQ(r0,r1,r2,r3))
3114
3115	def makeQuat(A,B,C):
3116	''' Make quaternion from rotation of A vector to B vector about C axis
3117
3118	:param np.array A,B,C: Cartesian 3-vectors
3119	:returns: quaternion & rotation angle in radians q=r+ai+bj+ck
3120	'''
3121
3122	V1 = np.cross(A,C)
3123	V2 = np.cross(B,C)
3124	if nl.norm(V1)*nl.norm(V2):
3125	V1 /= nl.norm(V1)
3126	V2 /= nl.norm(V2)
3127	V3 = np.cross(V1,V2)
3128	else:
3129	V3 = np.zeros(3)
3130	Q = np.array([0.,0.,0.,1.])
3131	D = 0.
3132	if nl.norm(V3):
3133	V3 /= nl.norm(V3)
3134	D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
3135	D = np.arccos(D1)/2.0
3136	V1 = C-V3
3137	V2 = C+V3
3138	DM = nl.norm(V1)
3139	DP = nl.norm(V2)
3140	S = np.sin(D)
3141	Q[0] = np.cos(D)
3142	Q[1:] = V3*S
3143	D *= 2.
3144	if DM > DP:
3145	D *= -1.
3146	return Q,D
3147
3148	def annealtests():
3149	from numpy import cos
3150	# minimum expected at ~-0.195
3151	func = lambda x: cos(14.5x-0.3) + (x+0.2)x
3152	print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
3153	print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
3154	print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
3155
3156	# minimum expected at ~[-0.195, -0.1]
3157	func = lambda x: cos(14.5x[0]-0.3) + (x[1]+0.2)x[1] + (x[0]+0.2)*x[0]
3158	print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='cauchy')
3159	print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='fast')
3160	print anneal(func,[1.0, 1.0],full_output=1,upper=[3.0, 3.0],lower=[-3.0, -3.0],feps=1e-4,maxiter=2000,schedule='boltzmann')
3161
3162
3163	if __name__ == '__main__':
3164	annealtests()

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

Download in other formats: