Context Navigation

source: trunk/GSASIImath.py @ 956

Last change on this file since 956 was 956, checked in by toby, 10 years ago
snapshot of CIF dev
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 93.7 KB

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