Context Navigation

source: trunk/GSASIImath.py @ 1056

Last change on this file since 1056 was 1056, checked in by vondreele, 10 years ago
restore multiprocessing MCSA
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 98.1 KB

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