Context Navigation

source: trunk/GSASIImath.py @ 1003

Last change on this file since 1003 was 1003, checked in by vondreele, 10 years ago
put Hessian eds calc back where they were (maybe a bit cleaner now).
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 96.1 KB

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