Context Navigation

source: trunk/GSASIImath.py @ 1780

Last change on this file since 1780 was 1780, checked in by vondreele, 8 years ago
force collapse of all children when main tree item is collapsed begin texture fitting from seq. refinements fix bugs in plotting if all SH coeff = 0 set background on Fcsq when delt = 3sig & 10 sig
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 123.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: