Context Navigation

source: trunk/GSASIImath.py @ 1063

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

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