Context Navigation

source: trunk/GSASIImath.py @ 989

Last change on this file since 989 was 989, checked in by toby, 10 years ago
add constraint derivs for single xtal; allow cellFill to work without esds; update docs; fix controls copy bug
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 96.8 KB

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