Context Navigation

source: trunk/GSASIImath.py @ 1055

Last change on this file since 1055 was 1055, checked in by vondreele, 10 years ago
working MC/SA
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 98.5 KB

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

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

Download in other formats: