Context Navigation

source: trunk/GSASIImath.py @ 1926

Last change on this file since 1926 was 1926, checked in by vondreele, 8 years ago
Work on add Hydrogens
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 132.7 KB

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

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

Download in other formats: