Context Navigation

source: trunk/GSASIImath.py @ 2278

Last change on this file since 2278 was 2275, checked in by vondreele, 9 years ago
add warning for high Levenberg-Marquardt lambda.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 161.0 KB

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

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

Download in other formats: