source: trunk/GSASIImath.py @ 1625

Last change on this file since 1625 was 1625, checked in by vondreele, 7 years ago

add a parameter to result from G2stIO.GetPhaseData?
add modulation functions to G2Math
add modulation names to G2obj
implement various wave types for modulations
plot position modulation wave function on contour plots
implement shift of modulation plot by +/-/0 keys
temporarily default G2spc.GetSSfxuinel to 1,-1 site symm.
move maxSSwave dict out of parmDict - now in controlDict
implement import of Sawtooth parms from J2K files.

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