Context Navigation

source: trunk/GSASIImath.py @ 1097

Last change on this file since 1097 was 1097, checked in by vondreele, 10 years ago
mods to restraint GUI & LS output for restraints correct torsion & Ramachandran penalties
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 100.1 KB

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