Context Navigation

source: trunk/GSASIImath.py @ 2060

Last change on this file since 2060 was 2060, checked in by vondreele, 7 years ago
revise 64bit fortran libraries for linux
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 150.9 KB

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

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

Download in other formats: