Context Navigation

source: trunk/GSASIImath.py @ 974

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