Context Navigation

source: trunk/GSASIImath.py @ 2079

Last change on this file since 2079 was 2079, checked in by vondreele, 7 years ago
new importer for hk6 files from 15ID
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 154.8 KB

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

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

Download in other formats: