Context Navigation

source: trunk/GSASIImath.py @ 1045

Last change on this file since 1045 was 1045, checked in by vondreele, 10 years ago
fix Fourier map problems - didn't like Fo<0. use ErrorDialog? for some errors.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 96.3 KB

Line
1	# -- coding: utf-8 --
2	#GSASIImath - major mathematics routines
3	########### SVN repository information ###################
4	# $Date: 2013-08-29 17:52:36 +0000 (Thu, 29 Aug 2013) $
5	# $Author: vondreele $
6	# $Revision: 1045 $
7	# $URL: trunk/GSASIImath.py $
8	# $Id: GSASIImath.py 1045 2013-08-29 17:52:36Z 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: 1045 $")
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	if valoff+esdoff < 0:
1301	valoff = esdoff = 0
1302	out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
1303	elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
1304	out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
1305	else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
1306	if abs(value) > 0:
1307	extra = -math.log10(abs(value))
1308	else:
1309	extra = 0
1310	if extra > 0: extra += 1
1311	out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
1312	if esd > 0:
1313	out += ("({:d})").format(intesd) # add the esd
1314	elif nTZ and '.' in out:
1315	out = out.rstrip('0') # strip zeros to right of decimal
1316	out = out.rstrip('.') # and decimal place when not needed
1317	if valoff != 0:
1318	out += ("e{:d}").format(valoff) # add an exponent, when needed
1319	return out
1320
1321	################################################################################
1322	##### Fourier & charge flip stuff
1323	################################################################################
1324
1325	def adjHKLmax(SGData,Hmax):
1326	'''default doc string
1327
1328	:param type name: description
1329
1330	:returns: type name: description
1331
1332	'''
1333	if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1334	Hmax[0] = ((Hmax[0]+3)/6)*6
1335	Hmax[1] = ((Hmax[1]+3)/6)*6
1336	Hmax[2] = ((Hmax[2]+1)/4)*4
1337	else:
1338	Hmax[0] = ((Hmax[0]+2)/4)*4
1339	Hmax[1] = ((Hmax[1]+2)/4)*4
1340	Hmax[2] = ((Hmax[2]+1)/4)*4
1341
1342	def OmitMap(data,reflData):
1343	'''default doc string
1344
1345	:param type name: description
1346
1347	:returns: type name: description
1348
1349	'''
1350	generalData = data['General']
1351	if not generalData['Map']['MapType']:
1352	print '**** ERROR - Fourier map not defined'
1353	return
1354	mapData = generalData['Map']
1355	dmin = mapData['Resolution']
1356	SGData = generalData['SGData']
1357	cell = generalData['Cell'][1:8]
1358	A = G2lat.cell2A(cell[:6])
1359	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1360	adjHKLmax(SGData,Hmax)
1361	Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1362	time0 = time.time()
1363	for ref in reflData:
1364	if ref[4] >= dmin:
1365	Fosq,Fcsq,ph = ref[8:11]
1366	for i,hkl in enumerate(ref[11]):
1367	hkl = np.asarray(hkl,dtype='i')
1368	dp = 360.*ref[12][i]
1369	a = cosd(ph+dp)
1370	b = sind(ph+dp)
1371	phasep = complex(a,b)
1372	phasem = complex(a,-b)
1373	F = np.sqrt(Fosq)
1374	h,k,l = hkl+Hmax
1375	Fhkl[h,k,l] = F*phasep
1376	h,k,l = -hkl+Hmax
1377	Fhkl[h,k,l] = F*phasem
1378	rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1379	print 'NB: this is just an Fobs map for now - under development'
1380	print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1381	print rho.shape
1382	mapData['rho'] = np.real(rho)
1383	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1384	return mapData
1385
1386	def FourierMap(data,reflData):
1387	'''default doc string
1388
1389	:param type name: description
1390
1391	:returns: type name: description
1392
1393	'''
1394	generalData = data['General']
1395	if not generalData['Map']['MapType']:
1396	print '**** ERROR - Fourier map not defined'
1397	return
1398	mapData = generalData['Map']
1399	dmin = mapData['Resolution']
1400	SGData = generalData['SGData']
1401	cell = generalData['Cell'][1:8]
1402	A = G2lat.cell2A(cell[:6])
1403	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1404	adjHKLmax(SGData,Hmax)
1405	Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1406	# Fhkl[0,0,0] = generalData['F000X']
1407	time0 = time.time()
1408	for ref in reflData:
1409	if ref[4] >= dmin:
1410	Fosq,Fcsq,ph = ref[8:11]
1411	for i,hkl in enumerate(ref[11]):
1412	hkl = np.asarray(hkl,dtype='i')
1413	dp = 360.*ref[12][i]
1414	a = cosd(ph+dp)
1415	b = sind(ph+dp)
1416	phasep = complex(a,b)
1417	phasem = complex(a,-b)
1418	if 'Fobs' in mapData['MapType']:
1419	F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1420	h,k,l = hkl+Hmax
1421	Fhkl[h,k,l] = F*phasep
1422	h,k,l = -hkl+Hmax
1423	Fhkl[h,k,l] = F*phasem
1424	elif 'Fcalc' in mapData['MapType']:
1425	F = np.sqrt(Fcsq)
1426	h,k,l = hkl+Hmax
1427	Fhkl[h,k,l] = F*phasep
1428	h,k,l = -hkl+Hmax
1429	Fhkl[h,k,l] = F*phasem
1430	elif 'delt-F' in mapData['MapType']:
1431	dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1432	h,k,l = hkl+Hmax
1433	Fhkl[h,k,l] = dF*phasep
1434	h,k,l = -hkl+Hmax
1435	Fhkl[h,k,l] = dF*phasem
1436	elif '2*Fo-Fc' in mapData['MapType']:
1437	F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1438	h,k,l = hkl+Hmax
1439	Fhkl[h,k,l] = F*phasep
1440	h,k,l = -hkl+Hmax
1441	Fhkl[h,k,l] = F*phasem
1442	elif 'Patterson' in mapData['MapType']:
1443	h,k,l = hkl+Hmax
1444	Fhkl[h,k,l] = complex(Fosq,0.)
1445	h,k,l = -hkl+Hmax
1446	Fhkl[h,k,l] = complex(Fosq,0.)
1447	rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1448	print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1449	mapData['rho'] = np.real(rho)
1450	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1451	return mapData
1452
1453	# map printing for testing purposes
1454	def printRho(SGLaue,rho,rhoMax):
1455	'''default doc string
1456
1457	:param type name: description
1458
1459	:returns: type name: description
1460
1461	'''
1462	dim = len(rho.shape)
1463	if dim == 2:
1464	ix,jy = rho.shape
1465	for j in range(jy):
1466	line = ''
1467	if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1468	line += (jy-j)*' '
1469	for i in range(ix):
1470	r = int(100*rho[i,j]/rhoMax)
1471	line += '%4d'%(r)
1472	print line+'\n'
1473	else:
1474	ix,jy,kz = rho.shape
1475	for k in range(kz):
1476	print 'k = ',k
1477	for j in range(jy):
1478	line = ''
1479	if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1480	line += (jy-j)*' '
1481	for i in range(ix):
1482	r = int(100*rho[i,j,k]/rhoMax)
1483	line += '%4d'%(r)
1484	print line+'\n'
1485	## keep this
1486
1487	def findOffset(SGData,A,Fhkl):
1488	'''default doc string
1489
1490	:param type name: description
1491
1492	:returns: type name: description
1493
1494	'''
1495	if SGData['SpGrp'] == 'P 1':
1496	return [0,0,0]
1497	hklShape = Fhkl.shape
1498	hklHalf = np.array(hklShape)/2
1499	sortHKL = np.argsort(Fhkl.flatten())
1500	Fdict = {}
1501	for hkl in sortHKL:
1502	HKL = np.unravel_index(hkl,hklShape)
1503	F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1504	if F == 0.:
1505	break
1506	Fdict['%.6f'%(np.absolute(F))] = hkl
1507	Flist = np.flipud(np.sort(Fdict.keys()))
1508	F = str(1.e6)
1509	i = 0
1510	DH = []
1511	Dphi = []
1512	Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1513	for F in Flist:
1514	hkl = np.unravel_index(Fdict[F],hklShape)
1515	iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1516	Uniq = np.array(Uniq,dtype='i')
1517	Phi = np.array(Phi)
1518	Uniq = np.concatenate((Uniq,-Uniq))+hklHalf # put in Friedel pairs & make as index to Farray
1519	Phi = np.concatenate((Phi,-Phi)) # and their phase shifts
1520	Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1521	ang0 = np.angle(Fh0,deg=True)/360.
1522	for H,phi in zip(Uniq,Phi)[1:]:
1523	ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
1524	dH = H-hkl
1525	dang = ang-ang0
1526	if np.any(np.abs(dH)-Hmax > 0): #keep low order DHs
1527	continue
1528	DH.append(dH)
1529	Dphi.append((dang+.5) % 1.0)
1530	if i > 20 or len(DH) > 30:
1531	break
1532	i += 1
1533	DH = np.array(DH)
1534	print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1535	Dphi = np.array(Dphi)
1536	steps = np.array(hklShape)
1537	X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1538	XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1539	Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
1540	Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1541	hist,bins = np.histogram(Mmap,bins=1000)
1542	# for i,item in enumerate(hist[:10]):
1543	# print item,bins[i]
1544	chisq = np.min(Mmap)
1545	DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1546	print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1547	# print (np.dot(DX,DH.T)+.5)%1.-Dphi
1548	return DX
1549
1550	def ChargeFlip(data,reflData,pgbar):
1551	'''default doc string
1552
1553	:param type name: description
1554
1555	:returns: type name: description
1556
1557	'''
1558	generalData = data['General']
1559	mapData = generalData['Map']
1560	flipData = generalData['Flip']
1561	FFtable = {}
1562	if 'None' not in flipData['Norm element']:
1563	normElem = flipData['Norm element'].upper()
1564	FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1565	for ff in FFs:
1566	if ff['Symbol'] == normElem:
1567	FFtable.update(ff)
1568	dmin = flipData['Resolution']
1569	SGData = generalData['SGData']
1570	cell = generalData['Cell'][1:8]
1571	A = G2lat.cell2A(cell[:6])
1572	Vol = cell[6]
1573	Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1574	adjHKLmax(SGData,Hmax)
1575	Ehkl = np.zeros(shape=2*Hmax,dtype='c16') #2X64bits per complex no.
1576	time0 = time.time()
1577	for ref in reflData:
1578	dsp = ref[4]
1579	if dsp >= dmin:
1580	ff = 0.1Vol #est. no. atoms for ~10A*3/atom
1581	if FFtable:
1582	SQ = 0.25/dsp**2
1583	ff *= G2el.ScatFac(FFtable,SQ)[0]
1584	if ref[8] > 0.: #use only +ve Fobs**2
1585	E = np.sqrt(ref[8])/ff
1586	else:
1587	E = 0.
1588	ph = ref[10]
1589	ph = rn.uniform(0.,360.)
1590	for i,hkl in enumerate(ref[11]):
1591	hkl = np.asarray(hkl,dtype='i')
1592	dp = 360.*ref[12][i]
1593	a = cosd(ph+dp)
1594	b = sind(ph+dp)
1595	phasep = complex(a,b)
1596	phasem = complex(a,-b)
1597	h,k,l = hkl+Hmax
1598	Ehkl[h,k,l] = E*phasep
1599	h,k,l = -hkl+Hmax #Friedel pair refl.
1600	Ehkl[h,k,l] = E*phasem
1601	# Ehkl[Hmax] = 0.00001 #this to preserve F[0,0,0]
1602	CEhkl = copy.copy(Ehkl)
1603	MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1604	Emask = ma.getmask(MEhkl)
1605	sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1606	Ncyc = 0
1607	old = np.seterr(all='raise')
1608	while True:
1609	CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1610	CEsig = np.std(CErho)
1611	CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1612	CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho) #solves U atom problem! make 20. adjustible
1613	CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1614	CFhkl = np.where(CFhkl,CFhkl,1.0) #avoid divide by zero
1615	phase = CFhkl/np.absolute(CFhkl)
1616	CEhkl = np.absolute(Ehkl)*phase
1617	Ncyc += 1
1618	sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1619	DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1620	Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1621	if Rcf < 5.:
1622	break
1623	GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1624	if not GoOn or Ncyc > 10000:
1625	break
1626	np.seterr(**old)
1627	print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1628	CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1629	print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1630	roll = findOffset(SGData,A,CEhkl)
1631
1632	mapData['Rcf'] = Rcf
1633	mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1634	mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1635	mapData['rollMap'] = [0,0,0]
1636	return mapData
1637
1638	def SearchMap(data):
1639	'''Does a search of a density map for peaks meeting the criterion of peak
1640	height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
1641	mapData is data['General']['mapData']; the map is also in mapData.
1642
1643	:param data: the phase data structure
1644
1645	:returns: (peaks,mags,dzeros) where
1646
1647	* peaks : ndarray
1648	x,y,z positions of the peaks found in the map
1649	* mags : ndarray
1650	the magnitudes of the peaks
1651	* dzeros : ndarray
1652	the distance of the peaks from the unit cell origin
1653
1654	'''
1655	rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1656
1657	norm = 1./(np.sqrt(3.)np.sqrt(2.np.pi)**3)
1658
1659	def noDuplicate(xyz,peaks,Amat):
1660	XYZ = np.inner(Amat,xyz)
1661	if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1662	print ' Peak',xyz,' <0.5A from another peak'
1663	return False
1664	return True
1665
1666	def fixSpecialPos(xyz,SGData,Amat):
1667	equivs = G2spc.GenAtom(xyz,SGData,Move=True)
1668	X = []
1669	xyzs = [equiv[0] for equiv in equivs]
1670	for x in xyzs:
1671	if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
1672	X.append(x)
1673	if len(X) > 1:
1674	return np.average(X,axis=0)
1675	else:
1676	return xyz
1677
1678	def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
1679	Mag,x0,y0,z0,sig = parms
1680	z = -((x0-rX)2+(y0-rY)2+(z0-rZ)*2)/(2.sig**2)
1681	# return normMagnp.exp(z)/(sigres*3) #not slower but some faults in LS
1682	return normMag(1.+z+z*2/2.)/(sigres**3)
1683
1684	def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
1685	Mag,x0,y0,z0,sig = parms
1686	M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1687	return M
1688
1689	def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
1690	Mag,x0,y0,z0,sig = parms
1691	dMdv = np.zeros(([5,]+list(rX.shape)))
1692	delt = .01
1693	for i in range(5):
1694	parms[i] -= delt
1695	rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1696	parms[i] += 2.*delt
1697	rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1698	parms[i] -= delt
1699	dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
1700	rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1701	Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
1702	dMdv = np.reshape(dMdv,(5,rX.size))
1703	Hess = np.inner(dMdv,dMdv)
1704
1705	return Vec,Hess
1706
1707	generalData = data['General']
1708	phaseName = generalData['Name']
1709	SGData = generalData['SGData']
1710	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1711	drawingData = data['Drawing']
1712	peaks = []
1713	mags = []
1714	dzeros = []
1715	try:
1716	mapData = generalData['Map']
1717	contLevel = mapData['cutOff']*mapData['rhoMax']/100.
1718	rho = copy.copy(mapData['rho']) #don't mess up original
1719	mapHalf = np.array(rho.shape)/2
1720	res = mapData['Resolution']
1721	incre = np.array(rho.shape,dtype=np.float)
1722	step = max(1.0,1./res)+1
1723	steps = np.array(3*[step,])
1724	except KeyError:
1725	print '**** ERROR - Fourier map not defined'
1726	return peaks,mags
1727	rhoMask = ma.array(rho,mask=(rho<contLevel))
1728	indices = (-1,0,1)
1729	rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1730	for roll in rolls:
1731	if np.any(roll):
1732	rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
1733	indx = np.transpose(rhoMask.nonzero())
1734	peaks = indx/incre
1735	mags = rhoMask[rhoMask.nonzero()]
1736	for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
1737	rho = rollMap(rho,ind)
1738	rMM = mapHalf-steps
1739	rMP = mapHalf+steps+1
1740	rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1741	peakInt = np.sum(rhoPeak)res*3
1742	rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1743	x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0] #magnitude, position & width(sig)
1744	result = HessianLSQ(peakFunc,x0,Hess=peakHess,
1745	args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
1746	x1 = result[0]
1747	if not np.any(x1 < 0):
1748	mag = x1[0]
1749	peak = (np.array(x1[1:4])-ind)/incre
1750	peak = fixSpecialPos(peak,SGData,Amat)
1751	rho = rollMap(rho,-ind)
1752	dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
1753	return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1754
1755	def sortArray(data,pos,reverse=False):
1756	'''data is a list of items
1757	sort by pos in list; reverse if True
1758	'''
1759	T = []
1760	for i,M in enumerate(data):
1761	T.append((M[pos],i))
1762	D = dict(zip(T,data))
1763	T.sort()
1764	if reverse:
1765	T.reverse()
1766	X = []
1767	for key in T:
1768	X.append(D[key])
1769	return X
1770
1771	def PeaksEquiv(data,Ind):
1772	'''Find the equivalent map peaks for those selected. Works on the
1773	contents of data['Map Peaks'].
1774
1775	:param data: the phase data structure
1776	:param list Ind: list of selected peak indices
1777	:returns: augmented list of peaks including those related by symmetry to the
1778	ones in Ind
1779
1780	'''
1781	def Duplicate(xyz,peaks,Amat):
1782	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1783	return True
1784	return False
1785
1786	generalData = data['General']
1787	cell = generalData['Cell'][1:7]
1788	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1789	A = G2lat.cell2A(cell)
1790	SGData = generalData['SGData']
1791	mapPeaks = data['Map Peaks']
1792	XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
1793	Indx = {}
1794	for ind in Ind:
1795	xyz = np.array(mapPeaks[ind][1:4])
1796	xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1797	# for x in xyzs: print x
1798	for jnd,xyz in enumerate(XYZ):
1799	Indx[jnd] = Duplicate(xyz,xyzs,Amat)
1800	Ind = []
1801	for ind in Indx:
1802	if Indx[ind]:
1803	Ind.append(ind)
1804	return Ind
1805
1806	def PeaksUnique(data,Ind):
1807	'''Finds the symmetry unique set of peaks from those selected. Works on the
1808	contents of data['Map Peaks'].
1809
1810	:param data: the phase data structure
1811	:param list Ind: list of selected peak indices
1812	:returns: the list of symmetry unique peaks from among those given in Ind
1813
1814	'''
1815	# XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
1816
1817	def noDuplicate(xyz,peaks,Amat):
1818	if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1819	return False
1820	return True
1821
1822	generalData = data['General']
1823	cell = generalData['Cell'][1:7]
1824	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1825	A = G2lat.cell2A(cell)
1826	SGData = generalData['SGData']
1827	mapPeaks = data['Map Peaks']
1828	Indx = {}
1829	XYZ = {}
1830	for ind in Ind:
1831	XYZ[ind] = np.array(mapPeaks[ind][1:4])
1832	Indx[ind] = True
1833	for ind in Ind:
1834	if Indx[ind]:
1835	xyz = XYZ[ind]
1836	for jnd in Ind:
1837	if ind != jnd and Indx[jnd]:
1838	Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
1839	xyzs = np.array([equiv[0] for equiv in Equiv])
1840	Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
1841	Ind = []
1842	for ind in Indx:
1843	if Indx[ind]:
1844	Ind.append(ind)
1845	return Ind
1846
1847	################################################################################
1848	##### single peak fitting profile fxn stuff
1849	################################################################################
1850
1851	def getCWsig(ins,pos):
1852	'''default doc string
1853
1854	:param type name: description
1855
1856	:returns: type name: description
1857
1858	'''
1859	tp = tand(pos/2.0)
1860	return ins['U']tp2+ins['V']tp+ins['W']
1861
1862	def getCWsigDeriv(pos):
1863	'''default doc string
1864
1865	:param type name: description
1866
1867	:returns: type name: description
1868
1869	'''
1870	tp = tand(pos/2.0)
1871	return tp**2,tp,1.0
1872
1873	def getCWgam(ins,pos):
1874	'''default doc string
1875
1876	:param type name: description
1877
1878	:returns: type name: description
1879
1880	'''
1881	return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
1882
1883	def getCWgamDeriv(pos):
1884	'''default doc string
1885
1886	:param type name: description
1887
1888	:returns: type name: description
1889
1890	'''
1891	return 1./cosd(pos/2.0),tand(pos/2.0)
1892
1893	def getTOFsig(ins,dsp):
1894	'''default doc string
1895
1896	:param type name: description
1897
1898	:returns: type name: description
1899
1900	'''
1901	return ins['sig-0']+ins['sig-1']dsp2+ins['sig-q']dsp
1902
1903	def getTOFsigDeriv(dsp):
1904	'''default doc string
1905
1906	:param type name: description
1907
1908	:returns: type name: description
1909
1910	'''
1911	return 1.0,dsp**2,dsp
1912
1913	def getTOFgamma(ins,dsp):
1914	'''default doc string
1915
1916	:param type name: description
1917
1918	:returns: type name: description
1919
1920	'''
1921	return ins['X']dsp+ins['Y']dsp**2
1922
1923	def getTOFgammaDeriv(dsp):
1924	'''default doc string
1925
1926	:param type name: description
1927
1928	:returns: type name: description
1929
1930	'''
1931	return dsp,dsp**2
1932
1933	def getTOFbeta(ins,dsp):
1934	'''default doc string
1935
1936	:param type name: description
1937
1938	:returns: type name: description
1939
1940	'''
1941	return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp
1942
1943	def getTOFbetaDeriv(dsp):
1944	'''default doc string
1945
1946	:param type name: description
1947
1948	:returns: type name: description
1949
1950	'''
1951	return 1.0,1./dsp**4,1./dsp
1952
1953	def getTOFalpha(ins,dsp):
1954	'''default doc string
1955
1956	:param type name: description
1957
1958	:returns: type name: description
1959
1960	'''
1961	return ins['alpha']/dsp
1962
1963	def getTOFalphaDeriv(dsp):
1964	'''default doc string
1965
1966	:param type name: description
1967
1968	:returns: type name: description
1969
1970	'''
1971	return 1./dsp
1972
1973	def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
1974	'''default doc string
1975
1976	:param type name: description
1977
1978	:returns: type name: description
1979
1980	'''
1981	ind = 0
1982	if useFit:
1983	ind = 1
1984	ins = {}
1985	if 'C' in Parms['Type'][0]: #CW data - TOF later in an elif
1986	for x in ['U','V','W','X','Y']:
1987	ins[x] = Parms[x][ind]
1988	if ifQ: #qplot - convert back to 2-theta
1989	pos = 2.0asind(poswave/(4*math.pi))
1990	sig = getCWsig(ins,pos)
1991	gam = getCWgam(ins,pos)
1992	XY = [pos,0, mag,1, sig,0, gam,0] #default refine intensity 1st
1993	else:
1994	if ifQ:
1995	dsp = 2.*np.pi/pos
1996	pos = Parms['difC']*dsp
1997	else:
1998	dsp = pos/Parms['difC'][1]
1999	if 'Pdabc' in Parms2:
2000	for x in ['sig-0','sig-1','sig-q','X','Y']:
2001	ins[x] = Parms[x][ind]
2002	Pdabc = Parms2['Pdabc'].T
2003	alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2004	bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2005	else:
2006	for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','X','Y']:
2007	ins[x] = Parms[x][ind]
2008	alp = getTOFalpha(ins,dsp)
2009	bet = getTOFbeta(ins,dsp)
2010	sig = getTOFsig(ins,dsp)
2011	gam = getTOFgamma(ins,dsp)
2012	XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2013	return XY
2014
2015	################################################################################
2016	##### MC/SA stuff
2017	################################################################################
2018
2019	#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2020	# Original Author: Travis Oliphant 2002
2021	# Bug-fixes in 2006 by Tim Leslie
2022
2023
2024	import numpy
2025	from numpy import asarray, tan, exp, ones, squeeze, sign, \
2026	all, log, sqrt, pi, shape, array, minimum, where
2027	from numpy import random
2028
2029	#__all__ = ['anneal']
2030
2031	_double_min = numpy.finfo(float).min
2032	_double_max = numpy.finfo(float).max
2033	class base_schedule(object):
2034	def __init__(self):
2035	self.dwell = 20
2036	self.learn_rate = 0.5
2037	self.lower = -10
2038	self.upper = 10
2039	self.Ninit = 50
2040	self.accepted = 0
2041	self.tests = 0
2042	self.feval = 0
2043	self.k = 0
2044	self.T = None
2045
2046	def init(self, **options):
2047	self.__dict__.update(options)
2048	self.lower = asarray(self.lower)
2049	self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2050	self.upper = asarray(self.upper)
2051	self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2052	self.k = 0
2053	self.accepted = 0
2054	self.feval = 0
2055	self.tests = 0
2056
2057	def getstart_temp(self, best_state):
2058	""" Find a matching starting temperature and starting parameters vector
2059	i.e. find x0 such that func(x0) = T0.
2060
2061	Parameters
2062	----------
2063	best_state : _state
2064	A _state object to store the function value and x0 found.
2065
2066	returns
2067	-------
2068	x0 : array
2069	The starting parameters vector.
2070	"""
2071
2072	assert(not self.dims is None)
2073	lrange = self.lower
2074	urange = self.upper
2075	fmax = _double_min
2076	fmin = _double_max
2077	for _ in range(self.Ninit):
2078	x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2079	fval = self.func(x0, *self.args)
2080	self.feval += 1
2081	if fval > fmax:
2082	fmax = fval
2083	if fval < fmin:
2084	fmin = fval
2085	best_state.cost = fval
2086	best_state.x = array(x0)
2087
2088	self.T0 = (fmax-fmin)*1.5
2089	return best_state.x
2090
2091	def accept_test(self, dE):
2092	T = self.T
2093	self.tests += 1
2094	if dE < 0:
2095	self.accepted += 1
2096	return 1
2097	p = exp(-dE*1.0/self.boltzmann/T)
2098	if (p > random.uniform(0.0, 1.0)):
2099	self.accepted += 1
2100	return 1
2101	return 0
2102
2103	def update_guess(self, x0):
2104	return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2105
2106	def update_temp(self, x0):
2107	pass
2108
2109
2110	# A schedule due to Lester Ingber modified to use bounds - OK
2111	class fast_sa(base_schedule):
2112	def init(self, **options):
2113	self.__dict__.update(options)
2114	if self.m is None:
2115	self.m = 1.0
2116	if self.n is None:
2117	self.n = 1.0
2118	self.c = self.m * exp(-self.n * self.quench)
2119
2120	def update_guess(self, x0):
2121	x0 = asarray(x0)
2122	u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2123	T = self.T
2124	xc = (sign(u-0.5)T((1+1.0/T)*abs(2u-1)-1.0)+1.0)/2.0
2125	xnew = xc*(self.upper - self.lower)+self.lower
2126	return xnew
2127	# y = sign(u-0.5)T((1+1.0/T)*abs(2u-1)-1.0)
2128	# xc = y*(self.upper - self.lower)
2129	# xnew = x0 + xc
2130	# return xnew
2131
2132	def update_temp(self):
2133	self.T = self.T0exp(-self.c self.k**(self.quench))
2134	self.k += 1
2135	return
2136
2137	class cauchy_sa(base_schedule): #modified to use bounds - not good
2138	def update_guess(self, x0):
2139	x0 = asarray(x0)
2140	numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2141	xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2142	xnew = xc*(self.upper - self.lower)+self.lower
2143	return xnew
2144	# numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2145	# xc = self.learn_rate * self.T * tan(numbers)
2146	# xnew = x0 + xc
2147	# return xnew
2148
2149	def update_temp(self):
2150	self.T = self.T0/(1+self.k)
2151	self.k += 1
2152	return
2153
2154	class boltzmann_sa(base_schedule):
2155	# def update_guess(self, x0):
2156	# std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2157	# x0 = asarray(x0)
2158	# xc = squeeze(random.normal(0, 1.0, size=self.dims))
2159	#
2160	# xnew = x0 + xcstdself.learn_rate
2161	# return xnew
2162
2163	def update_temp(self):
2164	self.k += 1
2165	self.T = self.T0 / log(self.k+1.0)
2166	return
2167
2168	class log_sa(base_schedule): #OK
2169
2170	def init(self,**options):
2171	self.__dict__.update(options)
2172
2173	# def update_guess(self,x0):
2174	# return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2175
2176	def update_temp(self):
2177	self.k += 1
2178	self.T = self.T0self.slope*self.k
2179
2180	class Tremayne_sa(base_schedule): #needs fixing for two steps
2181
2182	def init(self,**options):
2183	self.__dict__.update(options)
2184
2185	# def update_guess(self,x0):
2186	# x0 = np.asarray(x0)
2187	# u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
2188	# xnew = x0+u
2189	# return xnew
2190
2191	def update_temp(self):
2192	self.k += 1
2193	self.T = self.T0self.slope*k
2194
2195	class _state(object):
2196	def __init__(self):
2197	self.x = None
2198	self.cost = None
2199
2200	# TODO:
2201	# allow for general annealing temperature profile
2202	# in that case use update given by alpha and omega and
2203	# variation of all previous updates and temperature?
2204
2205	# Simulated annealing
2206
2207	def anneal(func, x0, args=(), schedule='fast', full_output=0,
2208	T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2209	boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2210	lower=-100, upper=100, dwell=50, slope=0.9,ranStart=True,dlg=None):
2211	"""Minimize a function using simulated annealing.
2212
2213	Schedule is a schedule class implementing the annealing schedule.
2214	Available ones are 'fast', 'cauchy', 'boltzmann'
2215
2216	:param callable func: f(x, \*args)
2217	Function to be optimized.
2218	:param ndarray x0:
2219	Initial guess.
2220	:param tuple args:
2221	Extra parameters to `func`.
2222	:param base_schedule schedule:
2223	Annealing schedule to use (a class).
2224	:param bool full_output:
2225	Whether to return optional outputs.
2226	:param float T0:
2227	Initial Temperature (estimated as 1.2 times the largest
2228	cost-function deviation over random points in the range).
2229	:param float Tf:
2230	Final goal temperature.
2231	:param int maxeval:
2232	Maximum function evaluations.
2233	:param int maxaccept:
2234	Maximum changes to accept.
2235	:param int maxiter:
2236	Maximum cooling iterations.
2237	:param float learn_rate:
2238	Scale constant for adjusting guesses.
2239	:param float boltzmann:
2240	Boltzmann constant in acceptance test
2241	(increase for less stringent test at each temperature).
2242	:param float feps:
2243	Stopping relative error tolerance for the function value in
2244	last four coolings.
2245	:param float quench,m,n:
2246	Parameters to alter fast_sa schedule.
2247	:param float/ndarray lower,upper:
2248	Lower and upper bounds on `x`.
2249	:param int dwell:
2250	The number of times to search the space at each temperature.
2251	:param float slope:
2252	Parameter for log schedule
2253	:param bool ranStart=True:
2254	False for fixed point start
2255
2256	:returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2257
2258	* xmin (ndarray): Point giving smallest value found.
2259	* Jmin (float): Minimum value of function found.
2260	* T (float): Final temperature.
2261	* feval (int): Number of function evaluations.
2262	* iters (int): Number of cooling iterations.
2263	* accept (int): Number of tests accepted.
2264	* retval (int): Flag indicating stopping condition:
2265
2266	* 0: Points no longer changing
2267	* 1: Cooled to final temperature
2268	* 2: Maximum function evaluations
2269	* 3: Maximum cooling iterations reached
2270	* 4: Maximum accepted query locations reached
2271	* 5: Final point not the minimum amongst encountered points
2272
2273	Notes:
2274	Simulated annealing is a random algorithm which uses no derivative
2275	information from the function being optimized. In practice it has
2276	been more useful in discrete optimization than continuous
2277	optimization, as there are usually better algorithms for continuous
2278	optimization problems.
2279
2280	Some experimentation by trying the difference temperature
2281	schedules and altering their parameters is likely required to
2282	obtain good performance.
2283
2284	The randomness in the algorithm comes from random sampling in numpy.
2285	To obtain the same results you can call numpy.random.seed with the
2286	same seed immediately before calling scipy.optimize.anneal.
2287
2288	We give a brief description of how the three temperature schedules
2289	generate new points and vary their temperature. Temperatures are
2290	only updated with iterations in the outer loop. The inner loop is
2291	over xrange(dwell), and new points are generated for
2292	every iteration in the inner loop. (Though whether the proposed
2293	new points are accepted is probabilistic.)
2294
2295	For readability, let d denote the dimension of the inputs to func.
2296	Also, let x_old denote the previous state, and k denote the
2297	iteration number of the outer loop. All other variables not
2298	defined below are input variables to scipy.optimize.anneal itself.
2299
2300	In the 'fast' schedule the updates are ::
2301
2302	u ~ Uniform(0, 1, size=d)
2303	y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2304	xc = y * (upper - lower)
2305	x_new = x_old + xc
2306
2307	c = n * exp(-n * quench)
2308	T_new = T0 * exp(-c * k**quench)
2309
2310
2311	In the 'cauchy' schedule the updates are ::
2312
2313	u ~ Uniform(-pi/2, pi/2, size=d)
2314	xc = learn_rate * T * tan(u)
2315	x_new = x_old + xc
2316
2317	T_new = T0 / (1+k)
2318
2319	In the 'boltzmann' schedule the updates are ::
2320
2321	std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2322	y ~ Normal(0, std, size=d)
2323	x_new = x_old + learn_rate * y
2324
2325	T_new = T0 / log(1+k)
2326
2327	"""
2328	x0 = asarray(x0)
2329	lower = asarray(lower)
2330	upper = asarray(upper)
2331
2332	schedule = eval(schedule+'_sa()')
2333	# initialize the schedule
2334	schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2335	learn_rate=learn_rate, lower=lower, upper=upper,
2336	m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2337
2338	current_state, last_state, best_state = _state(), _state(), _state()
2339	if T0 is None:
2340	x0 = schedule.getstart_temp(best_state)
2341	else:
2342	if ranStart:
2343	x0 = random.uniform(size=len(x0))*(upper-lower) + lower #comment to avoid random start
2344	best_state.x = None
2345	best_state.cost = numpy.Inf
2346
2347	last_state.x = asarray(x0).copy()
2348	fval = func(x0,*args)
2349	schedule.feval += 1
2350	last_state.cost = fval
2351	if last_state.cost < best_state.cost:
2352	best_state.cost = fval
2353	best_state.x = asarray(x0).copy()
2354	schedule.T = schedule.T0
2355	fqueue = [100, 300, 500, 700]
2356	iters = 1
2357	keepGoing = True
2358	while keepGoing:
2359	retval = 0
2360	for n in xrange(dwell):
2361	current_state.x = schedule.update_guess(last_state.x)
2362	current_state.cost = func(current_state.x,*args)
2363	schedule.feval += 1
2364
2365	dE = current_state.cost - last_state.cost
2366	if schedule.accept_test(dE):
2367	last_state.x = current_state.x.copy()
2368	last_state.cost = current_state.cost
2369	if last_state.cost < best_state.cost:
2370	best_state.x = last_state.x.copy()
2371	best_state.cost = last_state.cost
2372	if dlg:
2373	GoOn = dlg.Update(best_state.cost*100,
2374	newmsg='%s%8.5f\n%s%8.4f%s'%('Temperature =',schedule.T,'MC/SA Residual =',best_state.cost*100,'%'))[0]
2375	if not GoOn:
2376	best_state.x = last_state.x.copy()
2377	best_state.cost = last_state.cost
2378	retval = 5
2379	schedule.update_temp()
2380	iters += 1
2381	# Stopping conditions
2382	# 0) last saved values of f from each cooling step
2383	# are all very similar (effectively cooled)
2384	# 1) Tf is set and we are below it
2385	# 2) maxeval is set and we are past it
2386	# 3) maxiter is set and we are past it
2387	# 4) maxaccept is set and we are past it
2388	# 5) user canceled run via progress bar
2389
2390	fqueue.append(squeeze(last_state.cost))
2391	fqueue.pop(0)
2392	af = asarray(fqueue)*1.0
2393	if retval == 5:
2394	print ' User terminated run; incomplete MC/SA'
2395	keepGoing = False
2396	break
2397	if all(abs((af-af[0])/af[0]) < feps):
2398	retval = 0
2399	if abs(af[-1]-best_state.cost) > feps*10:
2400	retval = 5
2401	# print "Warning: Cooled to %f at %s but this is not" \
2402	# % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2403	# + " the smallest point found."
2404	break
2405	if (Tf is not None) and (schedule.T < Tf):
2406	retval = 1
2407	break
2408	if (maxeval is not None) and (schedule.feval > maxeval):
2409	retval = 2
2410	break
2411	if (iters > maxiter):
2412	print "Warning: Maximum number of iterations exceeded."
2413	retval = 3
2414	break
2415	if (maxaccept is not None) and (schedule.accepted > maxaccept):
2416	retval = 4
2417	break
2418
2419	if full_output:
2420	return best_state.x, best_state.cost, schedule.T, \
2421	schedule.feval, iters, schedule.accepted, retval
2422	else:
2423	return best_state.x, retval
2424
2425	def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q):
2426	outlist = []
2427	for n in range(iCyc):
2428	result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
2429	outlist.append(result[0])
2430	print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
2431	out_q.put(outlist)
2432
2433	def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
2434	import multiprocessing as mp
2435
2436	nprocs = mp.cpu_count()
2437	out_q = mp.Queue()
2438	procs = []
2439	iCyc = np.zeros(nprocs)
2440	for i in range(nCyc):
2441	iCyc[i%nprocs] += 1
2442	for i in range(nprocs):
2443	p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q))
2444	procs.append(p)
2445	p.start()
2446	resultlist = []
2447	for i in range(nprocs):
2448	resultlist += out_q.get()
2449	for p in procs:
2450	p.join()
2451	return resultlist
2452
2453	def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2454	'''default doc string
2455
2456	:param type name: description
2457
2458	:returns: type name: description
2459
2460	'''
2461
2462	twopi = 2.0*np.pi
2463	global tsum
2464	tsum = 0.
2465
2466	def getMDparms(item,pfx,parmDict,varyList):
2467	parmDict[pfx+'MDaxis'] = item['axis']
2468	parmDict[pfx+'MDval'] = item['Coef'][0]
2469	if item['Coef'][1]:
2470	varyList += [pfx+'MDval',]
2471	limits = item['Coef'][2]
2472	lower.append(limits[0])
2473	upper.append(limits[1])
2474
2475	def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2476	parmDict[pfx+'Atype'] = item['atType']
2477	aTypes \|= set([item['atType'],])
2478	pstr = ['Ax','Ay','Az']
2479	XYZ = [0,0,0]
2480	for i in range(3):
2481	name = pfx+pstr[i]
2482	parmDict[name] = item['Pos'][0][i]
2483	XYZ[i] = parmDict[name]
2484	if item['Pos'][1][i]:
2485	varyList += [name,]
2486	limits = item['Pos'][2][i]
2487	lower.append(limits[0])
2488	upper.append(limits[1])
2489	parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2490
2491	def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2492	parmDict[mfx+'MolCent'] = item['MolCent']
2493	parmDict[mfx+'RBId'] = item['RBId']
2494	pstr = ['Px','Py','Pz']
2495	ostr = ['Qa','Qi','Qj','Qk']
2496	for i in range(3):
2497	name = mfx+pstr[i]
2498	parmDict[name] = item['Pos'][0][i]
2499	if item['Pos'][1][i]:
2500	varyList += [name,]
2501	limits = item['Pos'][2][i]
2502	lower.append(limits[0])
2503	upper.append(limits[1])
2504	for i in range(4):
2505	name = mfx+ostr[i]
2506	parmDict[name] = item['Ori'][0][i]
2507	if item['Ovar'] == 'AV' and i:
2508	varyList += [name,]
2509	limits = item['Ori'][2][i]
2510	lower.append(limits[0])
2511	upper.append(limits[1])
2512	elif item['Ovar'] == 'A' and not i:
2513	varyList += [name,]
2514	limits = item['Ori'][2][i]
2515	lower.append(limits[0])
2516	upper.append(limits[1])
2517	if 'Tor' in item: #'Tor' not there for 'Vector' RBs
2518	for i in range(len(item['Tor'][0])):
2519	name = mfx+'Tor'+str(i)
2520	parmDict[name] = item['Tor'][0][i]
2521	if item['Tor'][1][i]:
2522	varyList += [name,]
2523	limits = item['Tor'][2][i]
2524	lower.append(limits[0])
2525	upper.append(limits[1])
2526	atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2527	aTypes \|= set(atypes)
2528	atNo += len(atypes)
2529	return atNo
2530
2531	def GetAtomM(Xdata,SGData):
2532	Mdata = []
2533	for xyz in Xdata:
2534	Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2535	return np.array(Mdata)
2536
2537	def GetAtomTX(RBdata,parmDict):
2538	'Needs a doc string'
2539	Bmat = parmDict['Bmat']
2540	atNo = parmDict['atNo']
2541	nfixAt = parmDict['nfixAt']
2542	Tdata = atNo*[' ',]
2543	Xdata = np.zeros((3,atNo))
2544	keys = {':Atype':Tdata,':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
2545	for iatm in range(nfixAt):
2546	for key in keys:
2547	parm = ':'+str(iatm)+key
2548	if parm in parmDict:
2549	if key == ':Atype':
2550	keys[key][iatm] = aTypes.index(parmDict[parm])
2551	else:
2552	keys[key][iatm] = parmDict[parm]
2553	iatm = nfixAt
2554	for iObj in range(parmDict['nObj']):
2555	pfx = str(iObj)+':'
2556	if parmDict[pfx+'Type'] in ['Vector','Residue']:
2557	if parmDict[pfx+'Type'] == 'Vector':
2558	RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2559	vecs = RBRes['rbVect']
2560	mags = RBRes['VectMag']
2561	Cart = np.zeros_like(vecs[0])
2562	for vec,mag in zip(vecs,mags):
2563	Cart += vec*mag
2564	elif parmDict[pfx+'Type'] == 'Residue':
2565	RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2566	Cart = np.array(RBRes['rbXYZ'])
2567	for itor,seq in enumerate(RBRes['rbSeq']):
2568	tName = pfx+'Tor'+str(itor)
2569	QuatA = AVdeg2Q(parmDict[tName],Cart[seq[0]]-Cart[seq[1]])
2570	for ride in seq[3]:
2571	Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
2572	if parmDict[pfx+'MolCent'][1]:
2573	Cart -= parmDict[pfx+'MolCent'][0]
2574	Qori = normQ(np.array([parmDict[pfx+'Qa'],parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']]))
2575	Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
2576	for i,x in enumerate(Cart):
2577	X = np.inner(Bmat,prodQVQ(Qori,x))+Pos
2578	for j in range(3):
2579	Xdata[j][iatm] = X[j]
2580	Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
2581	iatm += 1
2582	elif parmDict[pfx+'Type'] == 'Atom':
2583	atNo = parmDict[pfx+'atNo']
2584	for key in keys:
2585	parm = pfx+key[1:] #remove extra ':'
2586	if parm in parmDict:
2587	if key == ':Atype':
2588	keys[key][atNo] = aTypes.index(parmDict[parm])
2589	else:
2590	keys[key][atNo] = parmDict[parm]
2591	iatm += 1
2592	else:
2593	continue #skips March Dollase
2594	return Tdata,Xdata.T
2595
2596
2597	def mcsaCalc(values,refList,rcov,ifInv,RBdata,varyList,parmDict):
2598	''' Compute structure factors for all h,k,l for phase
2599	input:
2600	refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
2601	ParmDict:
2602	puts result F^2 in each ref[8] in refList
2603	'''
2604	def mcsaMDSFcalc(mul,FFs,Uniq,Phi):
2605	''' Calls fortran routine'''
2606	Icalc = pyd.pymcsamdsfcalc(ifInv,len(Tdata),Tdata,Mdata,Xdata.flatten(),
2607	MDval,MDaxis,Gmat,mul,len(FFs),FFs,len(Uniq),Uniq.flatten(),Phi)
2608	return Icalc
2609	global tsum
2610	parmDict.update(dict(zip(varyList,values)))
2611	Tdata,Xdata = GetAtomTX(RBdata,parmDict)
2612	Mdata = parmDict['Mdata']
2613	MDval = parmDict['0:MDval']
2614	MDaxis = parmDict['0:MDaxis']
2615	Gmat = parmDict['Gmat']
2616	Srefs = np.zeros(len(refList))
2617	sumFcsq = 0.
2618	for refl in refList:
2619	t0 = time.time()
2620	refl[5] = mcsaMDSFcalc(refl[3],refl[7],refl[8],refl[9])
2621	tsum += (time.time()-t0)
2622	sumFcsq += refl[5]
2623	scale = (parmDict['sumFosq']/sumFcsq)
2624	for iref,refl in enumerate(refList):
2625	refl[5] *= scale
2626	Srefs[iref] = refl[4]-refl[5]
2627	M = np.inner(Srefs,np.inner(rcov,Srefs))
2628	return M/parmDict['sumFosq']**2
2629
2630	sq8ln2 = np.sqrt(8*np.log(2))
2631	sq2pi = np.sqrt(2*np.pi)
2632	sq4pi = np.sqrt(4*np.pi)
2633	generalData = data['General']
2634	Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2635	Gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])[0]
2636	SGData = generalData['SGData']
2637	fixAtoms = data['Atoms'] #if any
2638	cx,ct,cs = generalData['AtomPtrs'][:3]
2639	aTypes = set([])
2640	parmDict = {'Bmat':Bmat,'Gmat':Gmat}
2641	varyList = []
2642	atNo = 0
2643	for atm in fixAtoms:
2644	pfx = ':'+str(atNo)+':'
2645	parmDict[pfx+'Atype'] = atm[ct]
2646	aTypes \|= set([atm[ct],])
2647	pstr = ['Ax','Ay','Az']
2648	parmDict[pfx+'Amul'] = atm[cs+1]
2649	for i in range(3):
2650	name = pfx+pstr[i]
2651	parmDict[name] = atm[cx+i]
2652	atNo += 1
2653	parmDict['nfixAt'] = len(fixAtoms)
2654	MCSA = generalData['MCSA controls']
2655	reflName = MCSA['Data source']
2656	phaseName = generalData['Name']
2657	MCSAObjs = data['MCSA']['Models'] #list of MCSA models
2658	upper = []
2659	lower = []
2660	for i,item in enumerate(MCSAObjs):
2661	mfx = str(i)+':'
2662	parmDict[mfx+'Type'] = item['Type']
2663	if item['Type'] == 'MD':
2664	getMDparms(item,mfx,parmDict,varyList)
2665	elif item['Type'] == 'Atom':
2666	getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
2667	parmDict[mfx+'atNo'] = atNo
2668	atNo += 1
2669	elif item['Type'] in ['Residue','Vector']:
2670	atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
2671	parmDict['atNo'] = atNo #total no. of atoms
2672	parmDict['nObj'] = len(MCSAObjs)
2673	aTypes = list(aTypes)
2674	Tdata,Xdata = GetAtomTX(RBdata,parmDict)
2675	parmDict['Mdata'] = GetAtomM(Xdata,SGData)
2676	FFtables = G2el.GetFFtable(aTypes)
2677	refs = []
2678	sumFosq = 0
2679	if 'PWDR' in reflName:
2680	for ref in reflData:
2681	h,k,l,m,d,pos,sig,gam,f = ref[:9]
2682	if d >= MCSA['dmin']:
2683	sig = G2pwd.getgamFW(sig,gam)/sq8ln2 #--> sig from FWHM
2684	SQ = 0.25/d**2
2685	Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2686	FFs = G2el.getFFvalues(FFtables,SQ,True)
2687	refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi])
2688	sumFosq += f*m
2689	nRef = len(refs)
2690	rcov = np.zeros((nRef,nRef))
2691	for iref,refI in enumerate(refs):
2692	rcov[iref][iref] = 1./(sq4pi*refI[6])
2693	for jref,refJ in enumerate(refs[:iref]):
2694	t1 = refI[6]2+refJ[6]2
2695	t2 = (refJ[5]-refI[5])*2/(2.t1)
2696	if t2 > 10.:
2697	rcov[iref][jref] = 0.
2698	else:
2699	rcov[iref][jref] = 1./(sq2pinp.sqrt(t1)np.exp(t2))
2700	rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2701	Rdiag = np.sqrt(np.diag(rcov))
2702	Rnorm = np.outer(Rdiag,Rdiag)
2703	rcov /= Rnorm
2704	elif 'Pawley' in reflName:
2705	covMatrix = covData['covMatrix']
2706	vList = covData['varyList']
2707	for iref,refI in enumerate(reflData):
2708	h,k,l,m,d,v,f,s = refI
2709	if d >= MCSA['dmin'] and v: #skip unrefined ones
2710	SQ = 0.25/d**2
2711	Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2712	FFs = G2el.getFFvalues(FFtables,SQ,True)
2713	refs.append([h,k,l,m,f*m,iref,0.,FFs,Uniq,phi])
2714	sumFosq += f*m
2715	nRef = len(refs)
2716	pfx = str(data['pId'])+'::PWLref:'
2717	rcov = np.zeros((nRef,nRef))
2718	for iref,refI in enumerate(refs):
2719	I = refI[5]
2720	nameI = pfx+str(I)
2721	if nameI in vList:
2722	Iindx = vList.index(nameI)
2723	rcov[iref][iref] = covMatrix[Iindx][Iindx]
2724	for jref,refJ in enumerate(refs[:iref]):
2725	J = refJ[5]
2726	nameJ = pfx+str(J)
2727	try:
2728	rcov[iref][jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)]
2729	except ValueError:
2730	rcov[iref][jref] = rcov[iref][jref-1]
2731	else:
2732	rcov[iref] = rcov[iref-1]
2733	rcov[iref][iref] = rcov[iref-1][iref-1]
2734	rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2735	Rdiag = np.sqrt(np.diag(rcov))
2736	Rnorm = np.outer(Rdiag,Rdiag)
2737	rcov /= Rnorm
2738	elif 'HKLF' in reflName:
2739	for ref in reflData:
2740	[h,k,l,m,d],f = ref[:5],ref[6]
2741	if d >= MCSA['dmin']:
2742	SQ = 0.25/d**2
2743	Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2744	FFs = G2el.getFFvalues(FFtables,SQ,True)
2745	refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
2746	sumFosq += f*m
2747	nRef = len(refs)
2748	rcov = np.identity(len(refs))
2749	print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
2750	print ' Number of parameters varied: %d'%(len(varyList))
2751	parmDict['sumFosq'] = sumFosq
2752	x0 = [parmDict[val] for val in varyList]
2753	ifInv = SGData['SGInv']
2754	results = anneal(mcsaCalc,x0,args=(refs,rcov,ifInv,RBdata,varyList,parmDict),
2755	schedule=MCSA['Algorithm'], full_output=True,
2756	T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
2757	boltzmann=MCSA['boltzmann'], learn_rate=0.5,
2758	quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
2759	lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',True),dlg=pgbar)
2760	Result = [False,False,results[1],results[2],]+list(results[0])
2761	Result.append(varyList)
2762	return Result,tsum
2763
2764
2765	################################################################################
2766	##### Quaternion stuff
2767	################################################################################
2768
2769	def prodQQ(QA,QB):
2770	''' Grassman quaternion product
2771	QA,QB quaternions; q=r+ai+bj+ck
2772	'''
2773	D = np.zeros(4)
2774	D[0] = QA[0]QB[0]-QA[1]QB[1]-QA[2]QB[2]-QA[3]QB[3]
2775	D[1] = QA[0]QB[1]+QA[1]QB[0]+QA[2]QB[3]-QA[3]QB[2]
2776	D[2] = QA[0]QB[2]-QA[1]QB[3]+QA[2]QB[0]+QA[3]QB[1]
2777	D[3] = QA[0]QB[3]+QA[1]QB[2]-QA[2]QB[1]+QA[3]QB[0]
2778
2779	# D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
2780	# D[1:] = QA[0]QB[1:]+QB[0]QA[1:]+np.cross(QA[1:],QB[1:])
2781
2782	return D
2783
2784	def normQ(QA):
2785	''' get length of quaternion & normalize it
2786	q=r+ai+bj+ck
2787	'''
2788	n = np.sqrt(np.sum(np.array(QA)**2))
2789	return QA/n
2790
2791	def invQ(Q):
2792	'''
2793	get inverse of quaternion
2794	q=r+ai+bj+ck; q* = r-ai-bj-ck
2795	'''
2796	return Q*np.array([1,-1,-1,-1])
2797
2798	def prodQVQ(Q,V):
2799	"""
2800	compute the quaternion vector rotation qvq-1 = v'
2801	q=r+ai+bj+ck
2802	"""
2803	VP = np.zeros(3)
2804	T2 = Q[0]*Q[1]
2805	T3 = Q[0]*Q[2]
2806	T4 = Q[0]*Q[3]
2807	T5 = -Q[1]*Q[1]
2808	T6 = Q[1]*Q[2]
2809	T7 = Q[1]*Q[3]
2810	T8 = -Q[2]*Q[2]
2811	T9 = Q[2]*Q[3]
2812	T10 = -Q[3]*Q[3]
2813	VP[0] = 2.((T8+T10)V[0]+(T6-T4)V[1]+(T3+T7)V[2])+V[0]
2814	VP[1] = 2.((T4+T6)V[0]+(T5+T10)V[1]+(T9-T2)V[2])+V[1]
2815	VP[2] = 2.((T7-T3)V[0]+(T2+T9)V[1]+(T5+T8)V[2])+V[2]
2816	return VP
2817
2818	def Q2Mat(Q):
2819	''' make rotation matrix from quaternion
2820	q=r+ai+bj+ck
2821	'''
2822	QN = normQ(Q)
2823	aa = QN[0]**2
2824	ab = QN[0]*QN[1]
2825	ac = QN[0]*QN[2]
2826	ad = QN[0]*QN[3]
2827	bb = QN[1]**2
2828	bc = QN[1]*QN[2]
2829	bd = QN[1]*QN[3]
2830	cc = QN[2]**2
2831	cd = QN[2]*QN[3]
2832	dd = QN[3]**2
2833	M = [[aa+bb-cc-dd, 2.(bc-ad), 2.(ac+bd)],
2834	[2(ad+bc), aa-bb+cc-dd, 2.(cd-ab)],
2835	[2(bd-ac), 2.(ab+cd), aa-bb-cc+dd]]
2836	return np.array(M)
2837
2838	def AV2Q(A,V):
2839	''' convert angle (radians) & vector to quaternion
2840	q=r+ai+bj+ck
2841	'''
2842	Q = np.zeros(4)
2843	d = np.sqrt(np.sum(np.array(V)**2))
2844	if d:
2845	V /= d
2846	else:
2847	V = np.array([0.,0.,1.])
2848	p = A/2.
2849	Q[0] = np.cos(p)
2850	Q[1:4] = V*np.sin(p)
2851	return Q
2852
2853	def AVdeg2Q(A,V):
2854	''' convert angle (degrees) & vector to quaternion
2855	q=r+ai+bj+ck
2856	'''
2857	Q = np.zeros(4)
2858	d = np.sqrt(np.sum(np.array(V)**2))
2859	if d:
2860	V /= d
2861	else:
2862	V = np.array([0.,0.,1.])
2863	p = A/2.
2864	Q[0] = cosd(p)
2865	Q[1:4] = V*sind(p)
2866	return Q
2867
2868	def Q2AVdeg(Q):
2869	''' convert quaternion to angle (degrees 0-360) & normalized vector
2870	q=r+ai+bj+ck
2871	'''
2872	A = 2.*acosd(Q[0])
2873	V = np.array(Q[1:])
2874	d = np.sqrt(np.sum(V**2))
2875	if d:
2876	V /= d
2877	else:
2878	A = 0.
2879	V = np.array([0.,0.,0.])
2880	return A,V
2881
2882	def Q2AV(Q):
2883	''' convert quaternion to angle (radians 0-2pi) & normalized vector
2884	q=r+ai+bj+ck
2885	'''
2886	A = 2.*np.arccos(Q[0])
2887	V = np.array(Q[1:])
2888	d = np.sqrt(np.sum(V**2))
2889	if d:
2890	V /= d
2891	else:
2892	A = 0.
2893	V = np.array([0.,0.,0.])
2894	return A,V
2895
2896	def makeQuat(A,B,C):
2897	''' Make quaternion from rotation of A vector to B vector about C axis
2898
2899	:param np.array A,B,C: Cartesian 3-vectors
2900	:returns: quaternion & rotation angle in radians q=r+ai+bj+ck
2901	'''
2902
2903	V1 = np.cross(A,C)
2904	V2 = np.cross(B,C)
2905	if nl.norm(V1)*nl.norm(V2):
2906	V1 /= nl.norm(V1)
2907	V2 /= nl.norm(V2)
2908	V3 = np.cross(V1,V2)
2909	else:
2910	V3 = np.zeros(3)
2911	Q = np.array([0.,0.,0.,1.])
2912	D = 0.
2913	if nl.norm(V3):
2914	V3 /= nl.norm(V3)
2915	D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
2916	D = np.arccos(D1)/2.0
2917	V1 = C-V3
2918	V2 = C+V3
2919	DM = nl.norm(V1)
2920	DP = nl.norm(V2)
2921	S = np.sin(D)
2922	Q[0] = np.cos(D)
2923	Q[1:] = V3*S
2924	D *= 2.
2925	if DM > DP:
2926	D *= -1.
2927	return Q,D
2928
2929	if __name__ == "__main__":
2930	from numpy import cos
2931	# minimum expected at ~-0.195
2932	func = lambda x: cos(14.5x-0.3) + (x+0.2)x
2933	print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
2934	print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
2935	print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
2936
2937	# minimum expected at ~[-0.195, -0.1]
2938	func = lambda x: cos(14.5x[0]-0.3) + (x[1]+0.2)x[1] + (x[0]+0.2)*x[0]
2939	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')
2940	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')
2941	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: