Context Navigation

source: trunk/GSASIImath.py @ 1977

Last change on this file since 1977 was 1977, checked in by vondreele, 7 years ago
More modulation structure factor stuff, begin SSF derivatives
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 144.7 KB

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

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

Download in other formats: