source: trunk/GSASIImath.py @ 948

Last change on this file since 948 was 948, checked in by toby, 9 years ago

update/reformat comments

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 81.0 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2013-06-11 15:47:27 +0000 (Tue, 11 Jun 2013) $
5# $Author: toby $
6# $Revision: 948 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 948 2013-06-11 15:47:27Z toby $
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: 948 $")
30import GSASIIElem as G2el
31import GSASIIlattice as G2lat
32import GSASIIspc as G2spc
33import numpy.fft as fft
34
35sind = lambda x: np.sin(x*np.pi/180.)
36cosd = lambda x: np.cos(x*np.pi/180.)
37tand = lambda x: np.tan(x*np.pi/180.)
38asind = lambda x: 180.*np.arcsin(x)/np.pi
39acosd = lambda x: 180.*np.arccos(x)/np.pi
40atand = lambda x: 180.*np.arctan(x)/np.pi
41atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
42
43def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
44   
45    """
46    Minimize the sum of squares of a set of equations.
47
48    ::
49   
50                    Nobs
51        x = arg min(sum(func(y)**2,axis=0))
52                    y=0
53
54    :param function func: callable method or function
55        should take at least one (possibly length N vector) argument and
56        returns M floating point numbers.
57    :param np.ndarray x0: The starting estimate for the minimization of length N
58    :param function Hess: callable method or function
59        A required function or method to compute the weighted vector and Hessian for func.
60        It must be a symmetric NxN array
61    :param tuple args: Any extra arguments to func are placed in this tuple.
62    :param float ftol: Relative error desired in the sum of squares.
63    :param float xtol: Relative error desired in the approximate solution.
64    :param int maxcyc: The maximum number of cycles of refinement to execute, if -1 refine
65        until other limits are met (ftol, xtol)
66
67    :Returns: (x,cov_x,infodict) where
68
69      * x : ndarray
70        The solution (or the result of the last iteration for an unsuccessful
71        call).
72      * cov_x : ndarray
73        Uses the fjac and ipvt optional outputs to construct an
74        estimate of the jacobian around the solution.  ``None`` if a
75        singular matrix encountered (indicates very flat curvature in
76        some direction).  This matrix must be multiplied by the
77        residual standard deviation to get the covariance of the
78        parameter estimates -- see curve_fit.
79      * infodict : dict
80        a dictionary of optional outputs with the keys:
81
82         * 'fvec' : the function evaluated at the output
83         * 'num cyc':
84         * 'nfev':
85         * 'lamMax':
86         * 'psing':
87           
88    """
89               
90    x0 = np.array(x0, ndmin=1)      #might be redundant?
91    n = len(x0)
92    if type(args) != type(()):
93        args = (args,)
94       
95    icycle = 0
96    One = np.ones((n,n))
97    lam = 0.001
98    lamMax = lam
99    nfev = 0
100    while icycle < maxcyc:
101        lamMax = max(lamMax,lam)
102        M = func(x0,*args)
103        nfev += 1
104        chisq0 = np.sum(M**2)
105        Yvec,Amat = Hess(x0,*args)
106        Adiag = np.sqrt(np.diag(Amat))
107        psing = np.where(np.abs(Adiag) < 1.e-14,True,False)
108        if np.any(psing):                #hard singularity in matrix
109            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
110        Anorm = np.outer(Adiag,Adiag)
111        Yvec /= Adiag
112        Amat /= Anorm       
113        while True:
114            Lam = np.eye(Amat.shape[0])*lam
115            Amatlam = Amat*(One+Lam)
116            try:
117                Xvec = nl.solve(Amatlam,Yvec)
118            except nl.LinAlgError:
119                print 'ouch #1'
120                psing = list(np.where(np.diag(nl.qr(Amatlam)[1]) < 1.e-14)[0])
121                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
122            Xvec /= Adiag
123            M2 = func(x0+Xvec,*args)
124            nfev += 1
125            chisq1 = np.sum(M2**2)
126            if chisq1 > chisq0:
127                lam *= 10.
128            else:
129                x0 += Xvec
130                lam /= 10.
131                break
132        if (chisq0-chisq1)/chisq0 < ftol:
133            break
134        icycle += 1
135    M = func(x0,*args)
136    nfev += 1
137    Yvec,Amat = Hess(x0,*args)
138    try:
139        Bmat = nl.inv(Amat)
140        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
141    except nl.LinAlgError:
142        print 'ouch #2 linear algebra error in LS'
143        psing = []
144        if maxcyc:
145            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
146        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
147   
148def getVCov(varyNames,varyList,covMatrix):
149    'Needs a doc string'
150    vcov = np.zeros((len(varyNames),len(varyNames)))
151    for i1,name1 in enumerate(varyNames):
152        for i2,name2 in enumerate(varyNames):
153            try:
154                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
155            except ValueError:
156                vcov[i1][i2] = 0.0
157    return vcov
158
159def FindAtomIndexByIDs(atomData,IDs,Draw=True):
160    'Needs a doc string'
161    indx = []
162    for i,atom in enumerate(atomData):
163        if Draw and atom[-3] in IDs:
164            indx.append(i)
165        elif atom[-1] in IDs:
166            indx.append(i)
167    return indx
168
169def FillAtomLookUp(atomData):
170    'Needs a doc string'
171    atomLookUp = {}
172    for iatm,atom in enumerate(atomData):
173        atomLookUp[atom[-1]] = iatm
174    return atomLookUp
175
176def GetAtomsById(atomData,atomLookUp,IdList):
177    'Needs a doc string'
178    atoms = []
179    for id in IdList:
180        atoms.append(atomData[atomLookUp[id]])
181    return atoms
182   
183def GetAtomItemsById(atomData,atomLookUp,IdList,itemLoc,numItems=1):
184    'Needs a doc string'
185    Items = []
186    if not isinstance(IdList,list):
187        IdList = [IdList,]
188    for id in IdList:
189        if numItems == 1:
190            Items.append(atomData[atomLookUp[id]][itemLoc])
191        else:
192            Items.append(atomData[atomLookUp[id]][itemLoc:itemLoc+numItems])
193    return Items
194   
195def GetAtomCoordsByID(pId,parmDict,AtLookup,indx):
196    'Needs a doc string'
197    pfx = [str(pId)+'::A'+i+':' for i in ['x','y','z']]
198    dpfx = [str(pId)+'::dA'+i+':' for i in ['x','y','z']]
199    XYZ = []
200    for ind in indx:
201        names = [pfx[i]+str(AtLookup[ind]) for i in range(3)]
202        dnames = [dpfx[i]+str(AtLookup[ind]) for i in range(3)]
203        XYZ.append([parmDict[name]+parmDict[dname] for name,dname in zip(names,dnames)])
204    return XYZ
205
206def AtomUij2TLS(atomData,atPtrs,Amat,Bmat,rbObj):   #unfinished & not used
207    'Needs a doc string; unfinished & not used'
208    for atom in atomData:
209        XYZ = np.inner(Amat,atom[cx:cx+3])
210        if atom[cia] == 'A':
211            UIJ = atom[cia+2:cia+8]
212               
213def TLS2Uij(xyz,g,Amat,rbObj):    #not used anywhere, but could be?
214    'Needs a doc string; not used anywhere'
215    TLStype,TLS = rbObj['ThermalMotion'][:2]
216    Tmat = np.zeros((3,3))
217    Lmat = np.zeros((3,3))
218    Smat = np.zeros((3,3))
219    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
220        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
221    if 'T' in TLStype:
222        Tmat = G2lat.U6toUij(TLS[:6])
223    if 'L' in TLStype:
224        Lmat = G2lat.U6toUij(TLS[6:12])
225    if 'S' in TLStype:
226        Smat = np.array([[TLS[18],TLS[12],TLS[13]],[TLS[14],TLS[19],TLS[15]],[TLS[16],TLS[17],0] ])
227    XYZ = np.inner(Amat,xyz)
228    Axyz = np.array([[ 0,XYZ[2],-XYZ[1]], [-XYZ[2],0,XYZ[0]], [XYZ[1],-XYZ[0],0]] )
229    Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
230    beta = np.inner(np.inner(g,Umat),g)
231    return G2lat.UijtoU6(beta)*gvec   
232       
233def AtomTLS2UIJ(atomData,atPtrs,Amat,rbObj):    #not used anywhere, but could be?
234    'Needs a doc string; not used anywhere'
235    cx,ct,cs,cia = atPtrs
236    TLStype,TLS = rbObj['ThermalMotion'][:2]
237    Tmat = np.zeros((3,3))
238    Lmat = np.zeros((3,3))
239    Smat = np.zeros((3,3))
240    G,g = G2lat.A2Gmat(Amat)
241    gvec = 1./np.sqrt(np.array([g[0][0],g[1][1],g[2][2],g[0][1],g[0][2],g[1][2]]))
242    if 'T' in TLStype:
243        Tmat = G2lat.U6toUij(TLS[:6])
244    if 'L' in TLStype:
245        Lmat = G2lat.U6toUij(TLS[6:12])
246    if 'S' in TLStype:
247        Smat = np.array([ [TLS[18],TLS[12],TLS[13]], [TLS[14],TLS[19],TLS[15]], [TLS[16],TLS[17],0] ])
248    for atom in atomData:
249        XYZ = np.inner(Amat,atom[cx:cx+3])
250        Axyz = np.array([ 0,XYZ[2],-XYZ[1], -XYZ[2],0,XYZ[0], XYZ[1],-XYZ[0],0],ndmin=2 )
251        if 'U' in TSLtype:
252            atom[cia+1] = TLS[0]
253            atom[cia] = 'I'
254        else:
255            atom[cia] = 'A'
256            Umat = Tmat+np.inner(Axyz,Smat)+np.inner(Smat.T,Axyz.T)+np.inner(np.inner(Axyz,Lmat),Axyz.T)
257            beta = np.inner(np.inner(g,Umat),g)
258            atom[cia+2:cia+8] = G2spc.U2Uij(beta/gvec)
259
260def GetXYZDist(xyz,XYZ,Amat):
261    ''' gets distance from position xyz to all XYZ, xyz & XYZ are np.array
262        and are in crystal coordinates; Amat is crystal to Cart matrix
263    '''
264    return np.sqrt(np.sum(np.inner(Amat,XYZ-xyz)**2,axis=0))
265
266def getAtomXYZ(atoms,cx):
267    'Needs a doc string; not used anywhere'
268    XYZ = []
269    for atom in atoms:
270        XYZ.append(atom[cx:cx+3])
271    return np.array(XYZ)
272
273def UpdateRBXYZ(Bmat,RBObj,RBData,RBType):
274    ''' Returns crystal coordinates for atoms described by RBObj
275    '''
276    RBRes = RBData[RBType][RBObj['RBId']]
277    if RBType == 'Vector':
278        vecs = RBRes['rbVect']
279        mags = RBRes['VectMag']
280        Cart = np.zeros_like(vecs[0])
281        for vec,mag in zip(vecs,mags):
282            Cart += vec*mag
283    elif RBType == 'Residue':
284        Cart = np.array(RBRes['rbXYZ'])
285        for tor,seq in zip(RBObj['Torsions'],RBRes['rbSeq']):
286            QuatA = AVdeg2Q(tor[0],Cart[seq[0]]-Cart[seq[1]])
287            for ride in seq[3]:
288                Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
289    XYZ = np.zeros_like(Cart)
290    for i,xyz in enumerate(Cart):
291        X = prodQVQ(RBObj['Orient'][0],xyz)
292        XYZ[i] = np.inner(Bmat,X)+RBObj['Orig'][0]
293    return XYZ,Cart
294
295def UpdateMCSAxyz(Bmat,MCSA):
296    'Needs a doc string'
297    xyz = []
298    atTypes = []
299    iatm = 0
300    for model in MCSA['Models'][1:]:        #skip the MD model
301        if model['Type'] == 'Atom':
302            xyz.append(model['Pos'][0])
303            atTypes.append(model['atType'])
304            iatm += 1
305        else:
306            RBRes = MCSA['rbData'][model['Type']][model['RBId']]
307            Pos = np.array(model['Pos'][0])
308            Qori = np.array(model['Ori'][0])
309            if model['Type'] == 'Vector':
310                vecs = RBRes['rbVect']
311                mags = RBRes['VectMag']
312                Cart = np.zeros_like(vecs[0])
313                for vec,mag in zip(vecs,mags):
314                    Cart += vec*mag
315            elif model['Type'] == 'Residue':
316                Cart = np.array(RBRes['rbXYZ'])
317                for itor,seq in enumerate(RBRes['rbSeq']):
318                    QuatA = AVdeg2Q(model['Tor'][0][itor],Cart[seq[0]]-Cart[seq[1]])
319                    for ride in seq[3]:
320                        Cart[ride] = prodQVQ(QuatA,Cart[ride]-Cart[seq[1]])+Cart[seq[1]]
321            if model['MolCent'][1]:
322                Cart -= model['MolCent'][0]
323            for i,x in enumerate(Cart):
324                xyz.append(np.inner(Bmat,prodQVQ(Qori,x))+Pos)
325                atType = RBRes['rbTypes'][i]
326                atTypes.append(atType)
327                iatm += 1
328    return np.array(xyz),atTypes
329   
330def SetMolCent(model,RBData):
331    'Needs a doc string'
332    rideList = []
333    RBRes = RBData[model['Type']][model['RBId']]
334    if model['Type'] == 'Vector':
335        vecs = RBRes['rbVect']
336        mags = RBRes['VectMag']
337        Cart = np.zeros_like(vecs[0])
338        for vec,mag in zip(vecs,mags):
339            Cart += vec*mag
340    elif model['Type'] == 'Residue':
341        Cart = np.array(RBRes['rbXYZ'])
342        for seq in RBRes['rbSeq']:
343            rideList += seq[3]
344    centList = set(range(len(Cart)))-set(rideList)
345    cent = np.zeros(3)
346    for i in centList:
347        cent += Cart[i]
348    model['MolCent'][0] = cent/len(centList)   
349   
350def UpdateRBUIJ(Bmat,Cart,RBObj):
351    ''' Returns atom I/A, Uiso or UIJ for atoms at XYZ as described by RBObj
352    '''
353    TLStype,TLS = RBObj['ThermalMotion'][:2]
354    T = np.zeros(6)
355    L = np.zeros(6)
356    S = np.zeros(8)
357    if 'T' in TLStype:
358        T = TLS[:6]
359    if 'L' in TLStype:
360        L = np.array(TLS[6:12])*(np.pi/180.)**2
361    if 'S' in TLStype:
362        S = np.array(TLS[12:])*(np.pi/180.)
363    g = nl.inv(np.inner(Bmat,Bmat))
364    gvec = np.sqrt(np.array([g[0][0]**2,g[1][1]**2,g[2][2]**2,
365        g[0][0]*g[1][1],g[0][0]*g[2][2],g[1][1]*g[2][2]]))
366    Uout = []
367    Q = RBObj['Orient'][0]
368    for X in Cart:
369        X = prodQVQ(Q,X)
370        if 'U' in TLStype:
371            Uout.append(['I',TLS[0],0,0,0,0,0,0])
372        elif not 'N' in TLStype:
373            U = [0,0,0,0,0,0]
374            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])
375            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])
376            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])
377            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]+  \
378                S[4]*X[0]-S[5]*X[1]-(S[6]+S[7])*X[2]
379            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]+  \
380                S[3]*X[2]-S[2]*X[0]+S[6]*X[1]
381            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]+  \
382                S[0]*X[1]-S[1]*X[2]+S[7]*X[0]
383            Umat = G2lat.U6toUij(U)
384            beta = np.inner(np.inner(Bmat.T,Umat),Bmat)
385            Uout.append(['A',0.0,]+list(G2lat.UijtoU6(beta)*gvec))
386        else:
387            Uout.append(['N',])
388    return Uout
389   
390def GetSHCoeff(pId,parmDict,SHkeys):
391    'Needs a doc string'
392    SHCoeff = {}
393    for shkey in SHkeys:
394        shname = str(pId)+'::'+shkey
395        SHCoeff[shkey] = parmDict[shname]
396    return SHCoeff
397       
398def getMass(generalData):
399    'Computes mass of unit cell contents'
400    mass = 0.
401    for i,elem in enumerate(generalData['AtomTypes']):
402        mass += generalData['NoAtoms'][elem]*generalData['AtomMass'][i]
403    return mass   
404
405def getDensity(generalData):
406    'Computes density of unit cell contents'   
407    mass = getMass(generalData)
408    Volume = generalData['Cell'][7]
409    density = mass/(0.6022137*Volume)
410    return density,Volume/mass
411   
412def getWave(Parms):
413    try:
414        return Parms['Lam'][1]
415    except KeyError:
416        return Parms['Lam1'][1]
417   
418################################################################################
419##### distance, angle, planes, torsion stuff stuff
420################################################################################
421
422def getSyXYZ(XYZ,ops,SGData):
423    'Needs a doc string'
424    XYZout = np.zeros_like(XYZ)
425    for i,[xyz,op] in enumerate(zip(XYZ,ops)):
426        if op == '1':
427            XYZout[i] = xyz
428        else:
429            oprs = op.split('+')
430            unit = [0,0,0]
431            if oprs[1]:
432                unit = np.array(list(eval(oprs[1])))
433            syop =int(oprs[0])
434            inv = syop/abs(syop)
435            syop *= inv
436            cent = syop/100
437            syop %= 100
438            syop -= 1
439            M,T = SGData['SGOps'][syop]
440            XYZout[i] = (np.inner(M,xyz)+T)*inv+SGData['SGCen'][cent]+unit
441    return XYZout
442   
443def getRestDist(XYZ,Amat):
444    'Needs a doc string'
445    return np.sqrt(np.sum(np.inner(Amat,(XYZ[1]-XYZ[0]))**2))
446   
447def getRestDeriv(Func,XYZ,Amat,ops,SGData):
448    'Needs a doc string'
449    deriv = np.zeros((len(XYZ),3))
450    dx = 0.00001
451    for j,xyz in enumerate(XYZ):
452        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
453            XYZ[j] -= x
454            d1 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
455            XYZ[j] += 2*x
456            d2 = Func(getSyXYZ(XYZ,ops,SGData),Amat)
457            XYZ[j] -= x
458            deriv[j][i] = (d1-d2)/(2*dx)
459    return deriv.flatten()
460
461def getRestAngle(XYZ,Amat):
462    'Needs a doc string'
463   
464    def calcVec(Ox,Tx,Amat):
465        return np.inner(Amat,(Tx-Ox))
466
467    VecA = calcVec(XYZ[1],XYZ[0],Amat)
468    VecA /= np.sqrt(np.sum(VecA**2))
469    VecB = calcVec(XYZ[1],XYZ[2],Amat)
470    VecB /= np.sqrt(np.sum(VecB**2))
471    edge = VecB-VecA
472    edge = np.sum(edge**2)
473    angle = (2.-edge)/2.
474    angle = max(angle,-1.)
475    return acosd(angle)
476   
477def getRestPlane(XYZ,Amat):
478    'Needs a doc string'
479    sumXYZ = np.zeros(3)
480    for xyz in XYZ:
481        sumXYZ += xyz
482    sumXYZ /= len(XYZ)
483    XYZ = np.array(XYZ)-sumXYZ
484    XYZ = np.inner(Amat,XYZ).T
485    Zmat = np.zeros((3,3))
486    for i,xyz in enumerate(XYZ):
487        Zmat += np.outer(xyz.T,xyz)
488    Evec,Emat = nl.eig(Zmat)
489    Evec = np.sqrt(Evec)/(len(XYZ)-3)
490    Order = np.argsort(Evec)
491    return Evec[Order[0]]
492   
493def getRestChiral(XYZ,Amat):   
494    'Needs a doc string'
495    VecA = np.empty((3,3))   
496    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
497    VecA[1] = np.inner(XYZ[2]-XYZ[0],Amat)
498    VecA[2] = np.inner(XYZ[3]-XYZ[0],Amat)
499    return nl.det(VecA)
500   
501def getRestTorsion(XYZ,Amat):
502    'Needs a doc string'
503    VecA = np.empty((3,3))
504    VecA[0] = np.inner(XYZ[1]-XYZ[0],Amat)
505    VecA[1] = np.inner(XYZ[2]-XYZ[1],Amat)
506    VecA[2] = np.inner(XYZ[3]-XYZ[2],Amat)
507    D = nl.det(VecA)
508    Mag = np.sqrt(np.sum(VecA*VecA,axis=1))
509    P12 = np.sum(VecA[0]*VecA[1])/(Mag[0]*Mag[1])
510    P13 = np.sum(VecA[0]*VecA[2])/(Mag[0]*Mag[2])
511    P23 = np.sum(VecA[1]*VecA[2])/(Mag[1]*Mag[2])
512    Ang = 1.0
513    if abs(P12) < 1.0 and abs(P23) < 1.0:
514        Ang = (P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2))
515    TOR = (acosd(Ang)*D/abs(D)+720.)%360.
516    return TOR
517   
518def calcTorsionEnergy(TOR,Coeff=[]):
519    'Needs a doc string'
520    sum = 0.
521    if len(Coeff):
522        cof = np.reshape(Coeff,(3,3)).T
523        delt = TOR-cof[1]
524        delt = np.where(delt<-180.,delt+360.,delt)
525        delt = np.where(delt>180.,delt-360.,delt)
526        term = -cof[2]*delt**2
527        val = cof[0]*np.exp(term/1000.0)
528        pMax = cof[0][np.argmin(val)]
529        Eval = np.sum(val)
530        sum = Eval-pMax
531    return sum,Eval
532
533def getTorsionDeriv(XYZ,Amat,Coeff):
534    'Needs a doc string'
535    deriv = np.zeros((len(XYZ),3))
536    dx = 0.00001
537    for j,xyz in enumerate(XYZ):
538        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
539            XYZ[j] -= x
540            tor = getRestTorsion(XYZ,Amat)
541            p,d1 = calcTorsionEnergy(tor,Coeff)
542            XYZ[j] += 2*x
543            tor = getRestTorsion(XYZ,Amat)
544            p,d2 = calcTorsionEnergy(tor,Coeff)           
545            XYZ[j] -= x
546            deriv[j][i] = (d2-d1)/(2*dx)
547    return deriv.flatten()
548
549def getRestRama(XYZ,Amat):
550    'Needs a doc string'
551    phi = getRestTorsion(XYZ[:5],Amat)
552    psi = getRestTorsion(XYZ[1:],Amat)
553    return phi,psi
554   
555def calcRamaEnergy(phi,psi,Coeff=[]):
556    'Needs a doc string'
557    sum = 0.
558    if len(Coeff):
559        cof = Coeff.T
560        dPhi = phi-cof[1]
561        dPhi = np.where(dPhi<-180.,dPhi+360.,dPhi)
562        dPhi = np.where(dPhi>180.,dPhi-360.,dPhi)
563        dPsi = psi-cof[2]
564        dPsi = np.where(dPsi<-180.,dPsi+360.,dPsi)
565        dPsi = np.where(dPsi>180.,dPsi-360.,dPsi)
566        val = -cof[3]*dPhi**2-cof[4]*dPsi**2-2.0*cof[5]*dPhi*dPsi
567        val = cof[0]*np.exp(val/1000.)
568        pMax = cof[0][np.argmin(val)]
569        Eval = np.sum(val)
570        sum = Eval-pMax
571    return sum,Eval
572
573def getRamaDeriv(XYZ,Amat,Coeff):
574    'Needs a doc string'
575    deriv = np.zeros((len(XYZ),3))
576    dx = 0.00001
577    for j,xyz in enumerate(XYZ):
578        for i,x in enumerate(np.array([[dx,0,0],[0,dx,0],[0,0,dx]])):
579            XYZ[j] -= x
580            phi,psi = getRestRama(XYZ,Amat)
581            p,d1 = calcRamaEnergy(phi,psi,Coeff)
582            XYZ[j] += 2*x
583            phi,psi = getRestRama(XYZ,Amat)
584            p,d2 = calcRamaEnergy(phi,psi,Coeff)
585            XYZ[j] -= x
586            deriv[j][i] = (d2-d1)/(2*dx)
587    return deriv.flatten()
588
589def getRestPolefig(ODFln,SamSym,Grid):
590    'Needs a doc string'
591    X,Y = np.meshgrid(np.linspace(1.,-1.,Grid),np.linspace(-1.,1.,Grid))
592    R,P = np.sqrt(X**2+Y**2).flatten(),atan2d(Y,X).flatten()
593    R = np.where(R <= 1.,2.*atand(R),0.0)
594    Z = np.zeros_like(R)
595    Z = G2lat.polfcal(ODFln,SamSym,R,P)
596    Z = np.reshape(Z,(Grid,Grid))
597    return np.reshape(R,(Grid,Grid)),np.reshape(P,(Grid,Grid)),Z
598
599def getRestPolefigDerv(HKL,Grid,SHCoeff):
600    'Needs a doc string'
601    pass
602       
603def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
604    'Needs a doc string'
605    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
606        TxT = inv*(np.inner(M,Tx)+T)+C+U
607        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
608       
609    inv = Top/abs(Top)
610    cent = abs(Top)/100
611    op = abs(Top)%100-1
612    M,T = SGData['SGOps'][op]
613    C = SGData['SGCen'][cent]
614    dx = .00001
615    deriv = np.zeros(6)
616    for i in [0,1,2]:
617        Oxyz[i] -= dx
618        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
619        Oxyz[i] += 2*dx
620        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
621        Oxyz[i] -= dx
622        Txyz[i] -= dx
623        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
624        Txyz[i] += 2*dx
625        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
626        Txyz[i] -= dx
627    return deriv
628   
629def getAngSig(VA,VB,Amat,SGData,covData={}):
630    'Needs a doc string'   
631    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
632        TxT = inv*(np.inner(M,Tx)+T)+C
633        TxT = G2spc.MoveToUnitCell(TxT)+U
634        return np.inner(Amat,(TxT-Ox))
635       
636    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
637        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
638        VecA /= np.sqrt(np.sum(VecA**2))
639        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
640        VecB /= np.sqrt(np.sum(VecB**2))
641        edge = VecB-VecA
642        edge = np.sum(edge**2)
643        angle = (2.-edge)/2.
644        angle = max(angle,-1.)
645        return acosd(angle)
646       
647    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
648    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
649    invA = invB = 1
650    invA = TopA/abs(TopA)
651    invB = TopB/abs(TopB)
652    centA = abs(TopA)/100
653    centB = abs(TopB)/100
654    opA = abs(TopA)%100-1
655    opB = abs(TopB)%100-1
656    MA,TA = SGData['SGOps'][opA]
657    MB,TB = SGData['SGOps'][opB]
658    CA = SGData['SGCen'][centA]
659    CB = SGData['SGCen'][centB]
660    if 'covMatrix' in covData:
661        covMatrix = covData['covMatrix']
662        varyList = covData['varyList']
663        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
664        dx = .00001
665        dadx = np.zeros(9)
666        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
667        for i in [0,1,2]:
668            OxA[i] -= dx
669            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
670            OxA[i] += 2*dx
671            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
672            OxA[i] -= dx
673           
674            TxA[i] -= dx
675            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
676            TxA[i] += 2*dx
677            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
678            TxA[i] -= dx
679           
680            TxB[i] -= dx
681            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
682            TxB[i] += 2*dx
683            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/(2*dx)
684            TxB[i] -= dx
685           
686        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
687        if sigAng < 0.01:
688            sigAng = 0.0
689        return Ang,sigAng
690    else:
691        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
692
693def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
694    'Needs a doc string'
695    def calcDist(Atoms,SyOps,Amat):
696        XYZ = []
697        for i,atom in enumerate(Atoms):
698            Inv,M,T,C,U = SyOps[i]
699            XYZ.append(np.array(atom[1:4]))
700            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
701            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
702        V1 = XYZ[1]-XYZ[0]
703        return np.sqrt(np.sum(V1**2))
704       
705    Inv = []
706    SyOps = []
707    names = []
708    for i,atom in enumerate(Oatoms):
709        names += atom[-1]
710        Op,unit = Atoms[i][-1]
711        inv = Op/abs(Op)
712        m,t = SGData['SGOps'][abs(Op)%100-1]
713        c = SGData['SGCen'][abs(Op)/100]
714        SyOps.append([inv,m,t,c,unit])
715    Dist = calcDist(Oatoms,SyOps,Amat)
716   
717    sig = -0.001
718    if 'covMatrix' in covData:
719        parmNames = []
720        dx = .00001
721        dadx = np.zeros(6)
722        for i in range(6):
723            ia = i/3
724            ix = i%3
725            Oatoms[ia][ix+1] += dx
726            a0 = calcDist(Oatoms,SyOps,Amat)
727            Oatoms[ia][ix+1] -= 2*dx
728            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
729        covMatrix = covData['covMatrix']
730        varyList = covData['varyList']
731        DistVcov = getVCov(names,varyList,covMatrix)
732        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
733        if sig < 0.001:
734            sig = -0.001
735   
736    return Dist,sig
737
738def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
739    'Needs a doc string'       
740    def calcAngle(Atoms,SyOps,Amat):
741        XYZ = []
742        for i,atom in enumerate(Atoms):
743            Inv,M,T,C,U = SyOps[i]
744            XYZ.append(np.array(atom[1:4]))
745            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
746            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
747        V1 = XYZ[1]-XYZ[0]
748        V1 /= np.sqrt(np.sum(V1**2))
749        V2 = XYZ[1]-XYZ[2]
750        V2 /= np.sqrt(np.sum(V2**2))
751        V3 = V2-V1
752        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
753        return acosd(cang)
754
755    Inv = []
756    SyOps = []
757    names = []
758    for i,atom in enumerate(Oatoms):
759        names += atom[-1]
760        Op,unit = Atoms[i][-1]
761        inv = Op/abs(Op)
762        m,t = SGData['SGOps'][abs(Op)%100-1]
763        c = SGData['SGCen'][abs(Op)/100]
764        SyOps.append([inv,m,t,c,unit])
765    Angle = calcAngle(Oatoms,SyOps,Amat)
766   
767    sig = -0.01
768    if 'covMatrix' in covData:
769        parmNames = []
770        dx = .00001
771        dadx = np.zeros(9)
772        for i in range(9):
773            ia = i/3
774            ix = i%3
775            Oatoms[ia][ix+1] += dx
776            a0 = calcAngle(Oatoms,SyOps,Amat)
777            Oatoms[ia][ix+1] -= 2*dx
778            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
779        covMatrix = covData['covMatrix']
780        varyList = covData['varyList']
781        AngVcov = getVCov(names,varyList,covMatrix)
782        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
783        if sig < 0.01:
784            sig = -0.01
785   
786    return Angle,sig
787
788def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
789    'Needs a doc string'
790    def calcTorsion(Atoms,SyOps,Amat):
791       
792        XYZ = []
793        for i,atom in enumerate(Atoms):
794            Inv,M,T,C,U = SyOps[i]
795            XYZ.append(np.array(atom[1:4]))
796            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
797            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
798        V1 = XYZ[1]-XYZ[0]
799        V2 = XYZ[2]-XYZ[1]
800        V3 = XYZ[3]-XYZ[2]
801        V1 /= np.sqrt(np.sum(V1**2))
802        V2 /= np.sqrt(np.sum(V2**2))
803        V3 /= np.sqrt(np.sum(V3**2))
804        M = np.array([V1,V2,V3])
805        D = nl.det(M)
806        Ang = 1.0
807        P12 = np.dot(V1,V2)
808        P13 = np.dot(V1,V3)
809        P23 = np.dot(V2,V3)
810        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
811        return Tors
812           
813    Inv = []
814    SyOps = []
815    names = []
816    for i,atom in enumerate(Oatoms):
817        names += atom[-1]
818        Op,unit = Atoms[i][-1]
819        inv = Op/abs(Op)
820        m,t = SGData['SGOps'][abs(Op)%100-1]
821        c = SGData['SGCen'][abs(Op)/100]
822        SyOps.append([inv,m,t,c,unit])
823    Tors = calcTorsion(Oatoms,SyOps,Amat)
824   
825    sig = -0.01
826    if 'covMatrix' in covData:
827        parmNames = []
828        dx = .00001
829        dadx = np.zeros(12)
830        for i in range(12):
831            ia = i/3
832            ix = i%3
833            Oatoms[ia][ix+1] -= dx
834            a0 = calcTorsion(Oatoms,SyOps,Amat)
835            Oatoms[ia][ix+1] += 2*dx
836            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
837            Oatoms[ia][ix+1] -= dx           
838        covMatrix = covData['covMatrix']
839        varyList = covData['varyList']
840        TorVcov = getVCov(names,varyList,covMatrix)
841        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
842        if sig < 0.01:
843            sig = -0.01
844   
845    return Tors,sig
846       
847def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
848    'Needs a doc string'   
849    def calcDist(Atoms,SyOps,Amat):
850        XYZ = []
851        for i,atom in enumerate(Atoms):
852            Inv,M,T,C,U = SyOps[i]
853            XYZ.append(np.array(atom[1:4]))
854            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
855            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
856        V1 = XYZ[1]-XYZ[0]
857        return np.sqrt(np.sum(V1**2))
858       
859    def calcAngle(Atoms,SyOps,Amat):
860        XYZ = []
861        for i,atom in enumerate(Atoms):
862            Inv,M,T,C,U = SyOps[i]
863            XYZ.append(np.array(atom[1:4]))
864            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
865            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
866        V1 = XYZ[1]-XYZ[0]
867        V1 /= np.sqrt(np.sum(V1**2))
868        V2 = XYZ[1]-XYZ[2]
869        V2 /= np.sqrt(np.sum(V2**2))
870        V3 = V2-V1
871        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
872        return acosd(cang)
873
874    def calcTorsion(Atoms,SyOps,Amat):
875       
876        XYZ = []
877        for i,atom in enumerate(Atoms):
878            Inv,M,T,C,U = SyOps[i]
879            XYZ.append(np.array(atom[1:4]))
880            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
881            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
882        V1 = XYZ[1]-XYZ[0]
883        V2 = XYZ[2]-XYZ[1]
884        V3 = XYZ[3]-XYZ[2]
885        V1 /= np.sqrt(np.sum(V1**2))
886        V2 /= np.sqrt(np.sum(V2**2))
887        V3 /= np.sqrt(np.sum(V3**2))
888        M = np.array([V1,V2,V3])
889        D = nl.det(M)
890        Ang = 1.0
891        P12 = np.dot(V1,V2)
892        P13 = np.dot(V1,V3)
893        P23 = np.dot(V2,V3)
894        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
895        return Tors
896           
897    Inv = []
898    SyOps = []
899    names = []
900    for i,atom in enumerate(Oatoms):
901        names += atom[-1]
902        Op,unit = Atoms[i][-1]
903        inv = Op/abs(Op)
904        m,t = SGData['SGOps'][abs(Op)%100-1]
905        c = SGData['SGCen'][abs(Op)/100]
906        SyOps.append([inv,m,t,c,unit])
907    M = len(Oatoms)
908    if M == 2:
909        Val = calcDist(Oatoms,SyOps,Amat)
910    elif M == 3:
911        Val = calcAngle(Oatoms,SyOps,Amat)
912    else:
913        Val = calcTorsion(Oatoms,SyOps,Amat)
914   
915    sigVals = [-0.001,-0.01,-0.01]
916    sig = sigVals[M-3]
917    if 'covMatrix' in covData:
918        parmNames = []
919        dx = .00001
920        N = M*3
921        dadx = np.zeros(N)
922        for i in range(N):
923            ia = i/3
924            ix = i%3
925            Oatoms[ia][ix+1] += dx
926            if M == 2:
927                a0 = calcDist(Oatoms,SyOps,Amat)
928            elif M == 3:
929                a0 = calcAngle(Oatoms,SyOps,Amat)
930            else:
931                a0 = calcTorsion(Oatoms,SyOps,Amat)
932            Oatoms[ia][ix+1] -= 2*dx
933            if M == 2:
934                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
935            elif M == 3:
936                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
937            else:
938                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
939        covMatrix = covData['covMatrix']
940        varyList = covData['varyList']
941        Vcov = getVCov(names,varyList,covMatrix)
942        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
943        if sig < sigVals[M-3]:
944            sig = sigVals[M-3]
945   
946    return Val,sig
947       
948def ValEsd(value,esd=0,nTZ=False):
949    '''Format a floating point number with a given level of precision or
950    with in crystallographic format with a "esd", as value(esd). If esd is
951    negative the number is formatted with the level of significant figures
952    appropriate if abs(esd) were the esd, but the esd is not included.
953    if the esd is zero, approximately 6 significant figures are printed.
954    nTZ=True causes "extra" zeros to be removed after the decimal place.
955    for example:
956
957      * "1.235(3)" for value=1.2346 & esd=0.003
958      * "1.235(3)e4" for value=12346. & esd=30
959      * "1.235(3)e6" for value=0.12346e7 & esd=3000
960      * "1.235" for value=1.2346 & esd=-0.003
961      * "1.240" for value=1.2395 & esd=-0.003
962      * "1.24" for value=1.2395 & esd=-0.003 with nTZ=True
963      * "1.23460" for value=1.2346 & esd=0.0
964
965    :param float value: number to be formatted
966    :param float esd: uncertainty or if esd < 0, specifies level of
967      precision to be shown e.g. esd=-0.01 gives 2 places beyond decimal
968    :param bool nTZ: True to remove trailing zeros (default is False)
969    :returns: value(esd) or value as a string
970
971    '''
972    # Note: this routine is Python 3 compatible -- I think
973    if esd != 0:
974        # transform the esd to a one or two digit integer
975        l = math.log10(abs(esd)) % 1
976        # cut off of 19 1.9==>(19) but 1.95==>(2) N.B. log10(1.95) = 0.2900...
977        if l < 0.290034611362518: l+= 1.       
978        intesd = int(round(10**l)) # esd as integer
979        # determine the number of digits offset for the esd
980        esdoff = int(round(math.log10(intesd*1./abs(esd))))
981    else:
982        esdoff = 5
983    valoff = 0
984    if esdoff < 0 or abs(value) > 1.0e6 or abs(value) < 1.0e-4: # use scientific notation
985        # where the digit offset is to the left of the decimal place or where too many
986        # digits are needed
987        if abs(value) > 1:
988            valoff = int(math.log10(abs(value)))
989        else:
990            valoff = int(math.log10(abs(value))-0.9999999)
991    if esd != 0:
992        out = ("{:."+str(valoff+esdoff)+"f}").format(value/10**valoff) # format the value
993    elif valoff != 0: # esd = 0; exponential notation ==> esdoff decimal places
994        out = ("{:."+str(esdoff)+"f}").format(value/10**valoff) # format the value
995    else: # esd = 0; non-exponential notation ==> esdoff+1 significant digits
996        extra = -math.log10(abs(value))
997        if extra > 0: extra += 1
998        out = ("{:."+str(max(0,esdoff+int(extra)))+"f}").format(value) # format the value
999    if esd > 0:
1000        out += ("({:d})").format(intesd)  # add the esd
1001    elif nTZ and '.' in out:
1002        out = out.rstrip('0')  # strip digits to right of decimal
1003    if valoff != 0:
1004        out += ("e{:d}").format(valoff) # add an exponent, when needed
1005    return out
1006
1007################################################################################
1008##### Fourier & charge flip stuff
1009################################################################################
1010
1011def adjHKLmax(SGData,Hmax):
1012    'Needs a doc string'       
1013    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1014        Hmax[0] = ((Hmax[0]+3)/6)*6
1015        Hmax[1] = ((Hmax[1]+3)/6)*6
1016        Hmax[2] = ((Hmax[2]+1)/4)*4
1017    else:
1018        Hmax[0] = ((Hmax[0]+2)/4)*4
1019        Hmax[1] = ((Hmax[1]+2)/4)*4
1020        Hmax[2] = ((Hmax[2]+1)/4)*4
1021
1022def OmitMap(data,reflData):
1023    'Needs a doc string'       
1024    generalData = data['General']
1025    if not generalData['Map']['MapType']:
1026        print '**** ERROR - Fourier map not defined'
1027        return
1028    mapData = generalData['Map']
1029    dmin = mapData['Resolution']
1030    SGData = generalData['SGData']
1031    cell = generalData['Cell'][1:8]       
1032    A = G2lat.cell2A(cell[:6])
1033    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1034    adjHKLmax(SGData,Hmax)
1035    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1036    time0 = time.time()
1037    for ref in reflData:
1038        if ref[4] >= dmin:
1039            Fosq,Fcsq,ph = ref[8:11]
1040            for i,hkl in enumerate(ref[11]):
1041                hkl = np.asarray(hkl,dtype='i')
1042                dp = 360.*ref[12][i]
1043                a = cosd(ph+dp)
1044                b = sind(ph+dp)
1045                phasep = complex(a,b)
1046                phasem = complex(a,-b)
1047                F = np.sqrt(Fosq)
1048                h,k,l = hkl+Hmax
1049                Fhkl[h,k,l] = F*phasep
1050                h,k,l = -hkl+Hmax
1051                Fhkl[h,k,l] = F*phasem
1052    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1053    print 'NB: this is just an Fobs map for now - under development'
1054    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1055    print rho.shape
1056    mapData['rho'] = np.real(rho)
1057    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1058    return mapData
1059
1060def FourierMap(data,reflData):
1061    'Needs a doc string'           
1062    generalData = data['General']
1063    if not generalData['Map']['MapType']:
1064        print '**** ERROR - Fourier map not defined'
1065        return
1066    mapData = generalData['Map']
1067    dmin = mapData['Resolution']
1068    SGData = generalData['SGData']
1069    cell = generalData['Cell'][1:8]       
1070    A = G2lat.cell2A(cell[:6])
1071    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1072    adjHKLmax(SGData,Hmax)
1073    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1074#    Fhkl[0,0,0] = generalData['F000X']
1075    time0 = time.time()
1076    for ref in reflData:
1077        if ref[4] >= dmin:
1078            Fosq,Fcsq,ph = ref[8:11]
1079            for i,hkl in enumerate(ref[11]):
1080                hkl = np.asarray(hkl,dtype='i')
1081                dp = 360.*ref[12][i]
1082                a = cosd(ph+dp)
1083                b = sind(ph+dp)
1084                phasep = complex(a,b)
1085                phasem = complex(a,-b)
1086                if 'Fobs' in mapData['MapType']:
1087                    F = np.sqrt(Fosq)
1088                    h,k,l = hkl+Hmax
1089                    Fhkl[h,k,l] = F*phasep
1090                    h,k,l = -hkl+Hmax
1091                    Fhkl[h,k,l] = F*phasem
1092                elif 'Fcalc' in mapData['MapType']:
1093                    F = np.sqrt(Fcsq)
1094                    h,k,l = hkl+Hmax
1095                    Fhkl[h,k,l] = F*phasep
1096                    h,k,l = -hkl+Hmax
1097                    Fhkl[h,k,l] = F*phasem
1098                elif 'delt-F' in mapData['MapType']:
1099                    dF = np.sqrt(Fosq)-np.sqrt(Fcsq)
1100                    h,k,l = hkl+Hmax
1101                    Fhkl[h,k,l] = dF*phasep
1102                    h,k,l = -hkl+Hmax
1103                    Fhkl[h,k,l] = dF*phasem
1104                elif '2*Fo-Fc' in mapData['MapType']:
1105                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1106                    h,k,l = hkl+Hmax
1107                    Fhkl[h,k,l] = F*phasep
1108                    h,k,l = -hkl+Hmax
1109                    Fhkl[h,k,l] = F*phasem
1110                elif 'Patterson' in mapData['MapType']:
1111                    h,k,l = hkl+Hmax
1112                    Fhkl[h,k,l] = complex(Fosq,0.)
1113                    h,k,l = -hkl+Hmax
1114                    Fhkl[h,k,l] = complex(Fosq,0.)
1115    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1116    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1117    mapData['rho'] = np.real(rho)
1118    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1119    return mapData
1120   
1121# map printing for testing purposes
1122def printRho(SGLaue,rho,rhoMax):                         
1123    'Needs a doc string'       
1124    dim = len(rho.shape)
1125    if dim == 2:
1126        ix,jy = rho.shape
1127        for j in range(jy):
1128            line = ''
1129            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1130                line += (jy-j)*'  '
1131            for i in range(ix):
1132                r = int(100*rho[i,j]/rhoMax)
1133                line += '%4d'%(r)
1134            print line+'\n'
1135    else:
1136        ix,jy,kz = rho.shape
1137        for k in range(kz):
1138            print 'k = ',k
1139            for j in range(jy):
1140                line = ''
1141                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1142                    line += (jy-j)*'  '
1143                for i in range(ix):
1144                    r = int(100*rho[i,j,k]/rhoMax)
1145                    line += '%4d'%(r)
1146                print line+'\n'
1147## keep this
1148               
1149def findOffset(SGData,A,Fhkl):   
1150    'Needs a doc string'       
1151    if SGData['SpGrp'] == 'P 1':
1152        return [0,0,0]   
1153    hklShape = Fhkl.shape
1154    hklHalf = np.array(hklShape)/2
1155    sortHKL = np.argsort(Fhkl.flatten())
1156    Fdict = {}
1157    for hkl in sortHKL:
1158        HKL = np.unravel_index(hkl,hklShape)
1159        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1160        if F == 0.:
1161            break
1162        Fdict['%.6f'%(np.absolute(F))] = hkl
1163    Flist = np.flipud(np.sort(Fdict.keys()))
1164    F = str(1.e6)
1165    i = 0
1166    DH = []
1167    Dphi = []
1168    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1169    for F in Flist:
1170        hkl = np.unravel_index(Fdict[F],hklShape)
1171        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1172        Uniq = np.array(Uniq,dtype='i')
1173        Phi = np.array(Phi)
1174        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
1175        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
1176        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1177        ang0 = np.angle(Fh0,deg=True)/360.
1178        for H,phi in zip(Uniq,Phi)[1:]:
1179            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
1180            dH = H-hkl
1181            dang = ang-ang0
1182            if np.any(np.abs(dH)-Hmax > 0):    #keep low order DHs
1183                continue
1184            DH.append(dH)
1185            Dphi.append((dang+.5) % 1.0)
1186        if i > 20 or len(DH) > 30:
1187            break
1188        i += 1
1189    DH = np.array(DH)
1190    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1191    Dphi = np.array(Dphi)
1192    steps = np.array(hklShape)
1193    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1194    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1195    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
1196    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1197    hist,bins = np.histogram(Mmap,bins=1000)
1198#    for i,item in enumerate(hist[:10]):
1199#        print item,bins[i]
1200    chisq = np.min(Mmap)
1201    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1202    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1203#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
1204    return DX
1205   
1206def ChargeFlip(data,reflData,pgbar):
1207    'Needs a doc string'       
1208    generalData = data['General']
1209    mapData = generalData['Map']
1210    flipData = generalData['Flip']
1211    FFtable = {}
1212    if 'None' not in flipData['Norm element']:
1213        normElem = flipData['Norm element'].upper()
1214        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1215        for ff in FFs:
1216            if ff['Symbol'] == normElem:
1217                FFtable.update(ff)
1218    dmin = flipData['Resolution']
1219    SGData = generalData['SGData']
1220    cell = generalData['Cell'][1:8]       
1221    A = G2lat.cell2A(cell[:6])
1222    Vol = cell[6]
1223    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1224    adjHKLmax(SGData,Hmax)
1225    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
1226    time0 = time.time()
1227    for ref in reflData:
1228        dsp = ref[4]
1229        if dsp >= dmin:
1230            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
1231            if FFtable:
1232                SQ = 0.25/dsp**2
1233                ff *= G2el.ScatFac(FFtable,SQ)[0]
1234            if ref[8] > 0.:         #use only +ve Fobs**2
1235                E = np.sqrt(ref[8])/ff
1236            else:
1237                E = 0.
1238            ph = ref[10]
1239            ph = rn.uniform(0.,360.)
1240            for i,hkl in enumerate(ref[11]):
1241                hkl = np.asarray(hkl,dtype='i')
1242                dp = 360.*ref[12][i]
1243                a = cosd(ph+dp)
1244                b = sind(ph+dp)
1245                phasep = complex(a,b)
1246                phasem = complex(a,-b)
1247                h,k,l = hkl+Hmax
1248                Ehkl[h,k,l] = E*phasep
1249                h,k,l = -hkl+Hmax       #Friedel pair refl.
1250                Ehkl[h,k,l] = E*phasem
1251#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
1252    CEhkl = copy.copy(Ehkl)
1253    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1254    Emask = ma.getmask(MEhkl)
1255    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1256    Ncyc = 0
1257    old = np.seterr(all='raise')
1258    while True:       
1259        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1260        CEsig = np.std(CErho)
1261        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1262        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem! make 20. adjustible
1263        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1264        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
1265        phase = CFhkl/np.absolute(CFhkl)
1266        CEhkl = np.absolute(Ehkl)*phase
1267        Ncyc += 1
1268        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1269        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1270        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1271        if Rcf < 5.:
1272            break
1273        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
1274        if not GoOn or Ncyc > 10000:
1275            break
1276    np.seterr(**old)
1277    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
1278    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
1279    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
1280    roll = findOffset(SGData,A,CEhkl)
1281       
1282    mapData['Rcf'] = Rcf
1283    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1284    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1285    mapData['rollMap'] = [0,0,0]
1286    return mapData
1287   
1288def SearchMap(data):
1289    'Needs a doc string'       
1290    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
1291   
1292    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
1293   
1294    def noDuplicate(xyz,peaks,Amat):
1295        XYZ = np.inner(Amat,xyz)
1296        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1297            print ' Peak',xyz,' <0.5A from another peak'
1298            return False
1299        return True
1300                           
1301    def fixSpecialPos(xyz,SGData,Amat):
1302        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
1303        X = []
1304        xyzs = [equiv[0] for equiv in equivs]
1305        for x in xyzs:
1306            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
1307                X.append(x)
1308        if len(X) > 1:
1309            return np.average(X,axis=0)
1310        else:
1311            return xyz
1312       
1313    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
1314        Mag,x0,y0,z0,sig = parms
1315        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
1316#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
1317        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
1318       
1319    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
1320        Mag,x0,y0,z0,sig = parms
1321        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1322        return M
1323       
1324    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
1325        Mag,x0,y0,z0,sig = parms
1326        dMdv = np.zeros(([5,]+list(rX.shape)))
1327        delt = .01
1328        for i in range(5):
1329            parms[i] -= delt
1330            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1331            parms[i] += 2.*delt
1332            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1333            parms[i] -= delt
1334            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
1335        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
1336        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
1337        dMdv = np.reshape(dMdv,(5,rX.size))
1338        Hess = np.inner(dMdv,dMdv)
1339       
1340        return Vec,Hess
1341       
1342    generalData = data['General']
1343    phaseName = generalData['Name']
1344    SGData = generalData['SGData']
1345    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1346    drawingData = data['Drawing']
1347    peaks = []
1348    mags = []
1349    dzeros = []
1350    try:
1351        mapData = generalData['Map']
1352        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
1353        rho = copy.copy(mapData['rho'])     #don't mess up original
1354        mapHalf = np.array(rho.shape)/2
1355        res = mapData['Resolution']
1356        incre = np.array(rho.shape,dtype=np.float)
1357        step = max(1.0,1./res)+1
1358        steps = np.array(3*[step,])
1359    except KeyError:
1360        print '**** ERROR - Fourier map not defined'
1361        return peaks,mags
1362    rhoMask = ma.array(rho,mask=(rho<contLevel))
1363    indices = (-1,0,1)
1364    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1365    for roll in rolls:
1366        if np.any(roll):
1367            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
1368    indx = np.transpose(rhoMask.nonzero())
1369    peaks = indx/incre
1370    mags = rhoMask[rhoMask.nonzero()]
1371    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
1372        rho = rollMap(rho,ind)
1373        rMM = mapHalf-steps
1374        rMP = mapHalf+steps+1
1375        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1376        peakInt = np.sum(rhoPeak)*res**3
1377        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
1378        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
1379        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
1380            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
1381        x1 = result[0]
1382        if not np.any(x1 < 0):
1383            mag = x1[0]
1384            peak = (np.array(x1[1:4])-ind)/incre
1385        peak = fixSpecialPos(peak,SGData,Amat)
1386        rho = rollMap(rho,-ind)       
1387    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
1388    return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T
1389   
1390def sortArray(data,pos,reverse=False):
1391    '''data is a list of items
1392    sort by pos in list; reverse if True
1393    '''
1394    T = []
1395    for i,M in enumerate(data):
1396        T.append((M[pos],i))
1397    D = dict(zip(T,data))
1398    T.sort()
1399    if reverse:
1400        T.reverse()
1401    X = []
1402    for key in T:
1403        X.append(D[key])
1404    return X
1405
1406def PeaksEquiv(data,Ind):
1407    'Needs a doc string'       
1408    def Duplicate(xyz,peaks,Amat):
1409        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1410            return True
1411        return False
1412                           
1413    generalData = data['General']
1414    cell = generalData['Cell'][1:7]
1415    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1416    A = G2lat.cell2A(cell)
1417    SGData = generalData['SGData']
1418    mapPeaks = data['Map Peaks']
1419    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
1420    Indx = {}
1421    for ind in Ind:
1422        xyz = np.array(mapPeaks[ind][1:4])
1423        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
1424#        for x in xyzs: print x
1425        for jnd,xyz in enumerate(XYZ):       
1426            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
1427    Ind = []
1428    for ind in Indx:
1429        if Indx[ind]:
1430            Ind.append(ind)
1431    return Ind
1432               
1433def PeaksUnique(data,Ind):
1434    'Needs a doc string'       
1435#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
1436
1437    def noDuplicate(xyz,peaks,Amat):
1438        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
1439            return False
1440        return True
1441                           
1442    generalData = data['General']
1443    cell = generalData['Cell'][1:7]
1444    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
1445    A = G2lat.cell2A(cell)
1446    SGData = generalData['SGData']
1447    mapPeaks = data['Map Peaks']
1448    Indx = {}
1449    XYZ = {}
1450    for ind in Ind:
1451        XYZ[ind] = np.array(mapPeaks[ind][1:4])
1452        Indx[ind] = True
1453    for ind in Ind:
1454        if Indx[ind]:
1455            xyz = XYZ[ind]
1456            for jnd in Ind:
1457                if ind != jnd and Indx[jnd]:                       
1458                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
1459                    xyzs = np.array([equiv[0] for equiv in Equiv])
1460                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
1461    Ind = []
1462    for ind in Indx:
1463        if Indx[ind]:
1464            Ind.append(ind)
1465    return Ind
1466   
1467def getCWsig(ins,pos):
1468    tp = tand(pos/2.0)
1469    return ins['U']*tp**2+ins['V']*tp+ins['W']
1470   
1471def getCWsigDeriv(pos):
1472    tp = tand(pos/2.0)
1473    return tp**2,tp,1.0
1474   
1475def getCWgam(ins,pos):
1476    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
1477   
1478def getCWgamDeriv(pos):
1479    return 1./cosd(pos/2.0),tand(pos/2.0)
1480   
1481def getTOFsig(ins,dsp):
1482    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-q']*dsp
1483   
1484def getTOFsigDeriv(dsp):
1485    return 1.0,dsp**2,dsp
1486   
1487def getTOFgamma(ins,dsp):
1488    return ins['X']*dsp+ins['Y']*dsp**2
1489   
1490def getTOFgammaDeriv(dsp):
1491    return dsp,dsp**2
1492   
1493def getTOFbeta(ins,dsp):
1494    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp
1495   
1496def getTOFbetaDeriv(dsp):
1497    return 1.0,1./dsp**4,1./dsp
1498   
1499def getTOFalpha(ins,dsp):
1500    return ins['alpha']/dsp
1501   
1502def getTOFalphaDeriv(dsp):
1503    return 1./dsp
1504   
1505def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
1506    'Needs a doc string'       
1507    ind = 0
1508    if useFit:
1509        ind = 1
1510    ins = {}
1511    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
1512        for x in ['U','V','W','X','Y']:
1513            ins[x] = Parms[x][ind]
1514        if ifQ:                              #qplot - convert back to 2-theta
1515            pos = 2.0*asind(pos*wave/(4*math.pi))
1516        sig = getCWsig(ins,pos)
1517        gam = getCWgam(ins,pos)           
1518        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
1519    else:
1520        if ifQ:
1521            dsp = 2.*np.pi/pos
1522            pos = Parms['difC']*dsp
1523        else:
1524            dsp = pos/Parms['difC'][1]
1525        if 'Pdabc' in Parms2:
1526            for x in ['sig-0','sig-1','sig-q','X','Y']:
1527                ins[x] = Parms[x][ind]
1528            Pdabc = Parms2['Pdabc'].T
1529            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1530            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1531        else:
1532            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','X','Y']:
1533                ins[x] = Parms[x][ind]
1534            alp = getTOFalpha(ins,dsp)
1535            bet = getTOFbeta(ins,dsp)
1536        sig = getTOFsig(ins,dsp)
1537        gam = getTOFgamma(ins,dsp)
1538        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
1539    return XY
1540   
1541################################################################################
1542##### MC/SA stuff
1543################################################################################
1544
1545#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
1546# Original Author: Travis Oliphant 2002
1547# Bug-fixes in 2006 by Tim Leslie
1548
1549
1550import numpy
1551from numpy import asarray, tan, exp, ones, squeeze, sign, \
1552     all, log, sqrt, pi, shape, array, minimum, where
1553from numpy import random
1554
1555__all__ = ['anneal']
1556
1557_double_min = numpy.finfo(float).min
1558_double_max = numpy.finfo(float).max
1559class base_schedule(object):
1560    def __init__(self):
1561        self.dwell = 20
1562        self.learn_rate = 0.5
1563        self.lower = -10
1564        self.upper = 10
1565        self.Ninit = 50
1566        self.accepted = 0
1567        self.tests = 0
1568        self.feval = 0
1569        self.k = 0
1570        self.T = None
1571
1572    def init(self, **options):
1573        self.__dict__.update(options)
1574        self.lower = asarray(self.lower)
1575        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
1576        self.upper = asarray(self.upper)
1577        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
1578        self.k = 0
1579        self.accepted = 0
1580        self.feval = 0
1581        self.tests = 0
1582
1583    def getstart_temp(self, best_state):
1584        """ Find a matching starting temperature and starting parameters vector
1585        i.e. find x0 such that func(x0) = T0.
1586
1587        Parameters
1588        ----------
1589        best_state : _state
1590            A _state object to store the function value and x0 found.
1591
1592        Returns
1593        -------
1594        x0 : array
1595            The starting parameters vector.
1596        """
1597
1598        assert(not self.dims is None)
1599        lrange = self.lower
1600        urange = self.upper
1601        fmax = _double_min
1602        fmin = _double_max
1603        for _ in range(self.Ninit):
1604            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
1605            fval = self.func(x0, *self.args)
1606            self.feval += 1
1607            if fval > fmax:
1608                fmax = fval
1609            if fval < fmin:
1610                fmin = fval
1611                best_state.cost = fval
1612                best_state.x = array(x0)
1613
1614        self.T0 = (fmax-fmin)*1.5
1615        return best_state.x
1616
1617    def accept_test(self, dE):
1618        T = self.T
1619        self.tests += 1
1620        if dE < 0:
1621            self.accepted += 1
1622            return 1
1623        p = exp(-dE*1.0/self.boltzmann/T)
1624        if (p > random.uniform(0.0, 1.0)):
1625            self.accepted += 1
1626            return 1
1627        return 0
1628
1629    def update_guess(self, x0):
1630        pass
1631
1632    def update_temp(self, x0):
1633        pass
1634
1635
1636#  A schedule due to Lester Ingber
1637class fast_sa(base_schedule):
1638    def init(self, **options):
1639        self.__dict__.update(options)
1640        if self.m is None:
1641            self.m = 1.0
1642        if self.n is None:
1643            self.n = 1.0
1644        self.c = self.m * exp(-self.n * self.quench)
1645
1646    def update_guess(self, x0):
1647        x0 = asarray(x0)
1648        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
1649        T = self.T
1650        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
1651        xc = y*(self.upper - self.lower)
1652        xnew = x0 + xc
1653        return xnew
1654
1655    def update_temp(self):
1656        self.T = self.T0*exp(-self.c * self.k**(self.quench))
1657        self.k += 1
1658        return
1659
1660class cauchy_sa(base_schedule):
1661    def update_guess(self, x0):
1662        x0 = asarray(x0)
1663        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
1664        xc = self.learn_rate * self.T * tan(numbers)
1665        xnew = x0 + xc
1666        return xnew
1667
1668    def update_temp(self):
1669        self.T = self.T0/(1+self.k)
1670        self.k += 1
1671        return
1672
1673class boltzmann_sa(base_schedule):
1674    def update_guess(self, x0):
1675        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
1676        x0 = asarray(x0)
1677        xc = squeeze(random.normal(0, 1.0, size=self.dims))
1678
1679        xnew = x0 + xc*std*self.learn_rate
1680        return xnew
1681
1682    def update_temp(self):
1683        self.k += 1
1684        self.T = self.T0 / log(self.k+1.0)
1685        return
1686
1687class log_sa(base_schedule):
1688
1689    def init(self,**options):
1690        self.__dict__.update(options)
1691       
1692    def update_guess(self,x0):
1693        x0 = np.asarray(x0)
1694        u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
1695        xnew = x0+u
1696        return xnew
1697       
1698    def update_temp(self):
1699        self.k += 1
1700        self.T = self.T0*self.slope**k
1701       
1702class Tremayne_sa(base_schedule):   #needs fixing for two steps
1703
1704    def init(self,**options):
1705        self.__dict__.update(options)
1706
1707    def update_guess(self,x0):
1708        x0 = np.asarray(x0)
1709        u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
1710        xnew = x0+u
1711        return xnew       
1712   
1713    def update_temp(self):
1714        self.k += 1
1715        self.T = self.T0*self.slope**k
1716   
1717class _state(object):
1718    def __init__(self):
1719        self.x = None
1720        self.cost = None
1721
1722# TODO:
1723#     allow for general annealing temperature profile
1724#     in that case use update given by alpha and omega and
1725#     variation of all previous updates and temperature?
1726
1727# Simulated annealing
1728
1729def anneal(func, x0, args=(), schedule='fast', full_output=0,
1730           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
1731           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
1732           lower=-100, upper=100, dwell=50, slope=0.9):
1733    """Minimize a function using simulated annealing.
1734
1735    Schedule is a schedule class implementing the annealing schedule.
1736    Available ones are 'fast', 'cauchy', 'boltzmann'
1737
1738    :param callable func: f(x, \*args)
1739        Function to be optimized.
1740    :param ndarray x0:
1741        Initial guess.
1742    :param tuple args:
1743        Extra parameters to `func`.
1744    :param base_schedule schedule:
1745        Annealing schedule to use (a class).
1746    :param bool full_output:
1747        Whether to return optional outputs.
1748    :param float T0:
1749        Initial Temperature (estimated as 1.2 times the largest
1750        cost-function deviation over random points in the range).
1751    :param float Tf:
1752        Final goal temperature.
1753    :param int maxeval:
1754        Maximum function evaluations.
1755    :param int maxaccept:
1756        Maximum changes to accept.
1757    :param int maxiter:
1758        Maximum cooling iterations.
1759    :param float learn_rate:
1760        Scale constant for adjusting guesses.
1761    :param float boltzmann:
1762        Boltzmann constant in acceptance test
1763        (increase for less stringent test at each temperature).
1764    :param float feps:
1765        Stopping relative error tolerance for the function value in
1766        last four coolings.
1767    :param float quench,m,n:
1768        Parameters to alter fast_sa schedule.
1769    :param float/ndarray lower,upper:
1770        Lower and upper bounds on `x`.
1771    :param int dwell:
1772        The number of times to search the space at each temperature.
1773    :param float slope:
1774        Parameter for log schedule
1775
1776    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
1777
1778     * xmin (ndarray): Point giving smallest value found.
1779     * Jmin (float): Minimum value of function found.
1780     * T (float): Final temperature.
1781     * feval (int): Number of function evaluations.
1782     * iters (int): Number of cooling iterations.
1783     * accept (int): Number of tests accepted.
1784     * retval (int): Flag indicating stopping condition:
1785
1786              *  0: Points no longer changing
1787              *  1: Cooled to final temperature
1788              *  2: Maximum function evaluations
1789              *  3: Maximum cooling iterations reached
1790              *  4: Maximum accepted query locations reached
1791              *  5: Final point not the minimum amongst encountered points
1792
1793    *Notes*:
1794    Simulated annealing is a random algorithm which uses no derivative
1795    information from the function being optimized. In practice it has
1796    been more useful in discrete optimization than continuous
1797    optimization, as there are usually better algorithms for continuous
1798    optimization problems.
1799
1800    Some experimentation by trying the difference temperature
1801    schedules and altering their parameters is likely required to
1802    obtain good performance.
1803
1804    The randomness in the algorithm comes from random sampling in numpy.
1805    To obtain the same results you can call numpy.random.seed with the
1806    same seed immediately before calling scipy.optimize.anneal.
1807
1808    We give a brief description of how the three temperature schedules
1809    generate new points and vary their temperature. Temperatures are
1810    only updated with iterations in the outer loop. The inner loop is
1811    over loop over xrange(dwell), and new points are generated for
1812    every iteration in the inner loop. (Though whether the proposed
1813    new points are accepted is probabilistic.)
1814
1815    For readability, let d denote the dimension of the inputs to func.
1816    Also, let x_old denote the previous state, and k denote the
1817    iteration number of the outer loop. All other variables not
1818    defined below are input variables to scipy.optimize.anneal itself.
1819
1820    In the 'fast' schedule the updates are ::
1821
1822        u ~ Uniform(0, 1, size=d)
1823        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
1824        xc = y * (upper - lower)
1825        x_new = x_old + xc
1826
1827        c = n * exp(-n * quench)
1828        T_new = T0 * exp(-c * k**quench)
1829
1830
1831    In the 'cauchy' schedule the updates are ::
1832
1833        u ~ Uniform(-pi/2, pi/2, size=d)
1834        xc = learn_rate * T * tan(u)
1835        x_new = x_old + xc
1836
1837        T_new = T0 / (1+k)
1838
1839    In the 'boltzmann' schedule the updates are ::
1840
1841        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
1842        y ~ Normal(0, std, size=d)
1843        x_new = x_old + learn_rate * y
1844
1845        T_new = T0 / log(1+k)
1846
1847    """
1848    x0 = asarray(x0)
1849    lower = asarray(lower)
1850    upper = asarray(upper)
1851
1852    schedule = eval(schedule+'_sa()')
1853    #   initialize the schedule
1854    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
1855                  learn_rate=learn_rate, lower=lower, upper=upper,
1856                  m=m, n=n, quench=quench, dwell=dwell)
1857
1858    current_state, last_state, best_state = _state(), _state(), _state()
1859    if T0 is None:
1860        x0 = schedule.getstart_temp(best_state)
1861    else:
1862        best_state.x = None
1863        best_state.cost = numpy.Inf
1864
1865    last_state.x = asarray(x0).copy()
1866    fval = func(x0,*args)
1867    schedule.feval += 1
1868    last_state.cost = fval
1869    if last_state.cost < best_state.cost:
1870        best_state.cost = fval
1871        best_state.x = asarray(x0).copy()
1872    schedule.T = schedule.T0
1873    fqueue = [100, 300, 500, 700]
1874    iters = 0
1875    while 1:
1876        for n in xrange(dwell):
1877            current_state.x = schedule.update_guess(last_state.x)
1878            current_state.cost = func(current_state.x,*args)
1879            schedule.feval += 1
1880
1881            dE = current_state.cost - last_state.cost
1882            if schedule.accept_test(dE):
1883                last_state.x = current_state.x.copy()
1884                last_state.cost = current_state.cost
1885                if last_state.cost < best_state.cost:
1886                    best_state.x = last_state.x.copy()
1887                    best_state.cost = last_state.cost
1888        schedule.update_temp()
1889        iters += 1
1890        # Stopping conditions
1891        # 0) last saved values of f from each cooling step
1892        #     are all very similar (effectively cooled)
1893        # 1) Tf is set and we are below it
1894        # 2) maxeval is set and we are past it
1895        # 3) maxiter is set and we are past it
1896        # 4) maxaccept is set and we are past it
1897
1898        fqueue.append(squeeze(last_state.cost))
1899        fqueue.pop(0)
1900        af = asarray(fqueue)*1.0
1901        if all(abs((af-af[0])/af[0]) < feps):
1902            retval = 0
1903            if abs(af[-1]-best_state.cost) > feps*10:
1904                retval = 5
1905                print "Warning: Cooled to %f at %s but this is not" \
1906                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
1907                      + " the smallest point found."
1908            break
1909        if (Tf is not None) and (schedule.T < Tf):
1910            retval = 1
1911            break
1912        if (maxeval is not None) and (schedule.feval > maxeval):
1913            retval = 2
1914            break
1915        if (iters > maxiter):
1916            print "Warning: Maximum number of iterations exceeded."
1917            retval = 3
1918            break
1919        if (maxaccept is not None) and (schedule.accepted > maxaccept):
1920            retval = 4
1921            break
1922
1923    if full_output:
1924        return best_state.x, best_state.cost, schedule.T, \
1925               schedule.feval, iters, schedule.accepted, retval
1926    else:
1927        return best_state.x, retval
1928
1929
1930def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
1931    gamFW = lambda s,g: math.exp(math.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.)
1932   
1933    def getMDparms(item,pfx,parmDict,varyList):
1934        parmDict[pfx+'MDaxis'] = item['axis']
1935        parmDict[pfx+'MDval'] = item['Coef'][0]
1936        if item['Coef'][1]:
1937            varyList += [pfx+'MDval',]
1938            limits = item['Coef'][2]
1939            lower.append(limits[0])
1940            upper.append(limits[1])
1941           
1942    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
1943        parmDict[pfx+'Atype'] = item['atType']
1944        aTypes |= set([item['atType'],]) 
1945        pstr = ['Ax','Ay','Az']
1946        XYZ = [0,0,0]
1947        for i in range(3):
1948            name = pfx+pstr[i]
1949            parmDict[name] = item['Pos'][0][i]
1950            XYZ[i] = parmDict[name]
1951            if item['Pos'][1][i]:
1952                varyList += [name,]
1953                limits = item['Pos'][2][i]
1954                lower.append(limits[0])
1955                upper.append(limits[1])
1956        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
1957           
1958    def getRBparms(item,mfx,aTypes,RBdata,atNo,parmDict,varyList):
1959        parmDict[mfx+'MolCent'] = item['MolCent']
1960        parmDict[mfx+'RBId'] = item['RBId']
1961        pstr = ['Px','Py','Pz']
1962        ostr = ['Qa','Qi','Qj','Qk']
1963        for i in range(3):
1964            name = mfx+pstr[i]
1965            parmDict[name] = item['Pos'][0][i]
1966            if item['Pos'][1][i]:
1967                varyList += [name,]
1968                limits = item['Pos'][2][i]
1969                lower.append(limits[0])
1970                upper.append(limits[1])
1971        for i in range(4):
1972            name = mfx+ostr[i]
1973            parmDict[name] = item['Ori'][0][i]
1974            if item['Ovar'] == 'AV' and i:
1975                varyList += [name,]
1976                limits = item['Ori'][2][i]
1977                lower.append(limits[0])
1978                upper.append(limits[1])
1979            elif item['Ovar'] == 'A' and not i:
1980                varyList += [name,]
1981                limits = item['Ori'][2][i]
1982                lower.append(limits[0])
1983                upper.append(limits[1])
1984        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
1985            for i in range(len(item['Tor'][0])):
1986                name = mfx+'Tor'+str(i)
1987                parmDict[name] = item['Tor'][0][i]
1988                if item['Tor'][1][i]:
1989                    varyList += [name,]
1990                    limits = item['Tor'][2][i]
1991                    lower.append(limits[0])
1992                    upper.append(limits[1])
1993        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
1994        aTypes |= set(atypes)
1995        atNo += len(atypes)
1996        return atNo
1997
1998    def GetAtomTMX(pfx,RBdata,parmDict):
1999        'Needs a doc string'
2000        atNo = parmDIct['atNo']
2001        nfixAt = parmDict['nfixAt']
2002        Tdata = atNo*[' ',]
2003        Mdata = np.zeros(atNo)
2004        Fdata = np.zeros(atNo)
2005        Xdata = np.zeros((3,atNo))
2006        keys = {'Atype:':Tdata,'Amul:':Mdata,
2007            'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2]}
2008        nObjs = parmDict['nObjs']
2009        for iatm in range(nfixAt):
2010            for key in keys:
2011                parm = pfx+key+str(iatm)
2012                if parm in parmDict:
2013                    keys[key][iatm] = parmDict[parm]
2014        return Tdata,Mdata,Xdata
2015   
2016#    def mcsaSfCalc(refList,RBdata,G,SGData,parmDict):
2017#        ''' Compute structure factors for all h,k,l for phase
2018#        input:
2019#            refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
2020#            G:      reciprocal metric tensor
2021#            SGData: space group info. dictionary output from SpcGroup
2022#            ParmDict:
2023#        puts result F^2 in each ref[8] in refList
2024#        '''       
2025#        twopi = 2.0*np.pi
2026#        twopisq = 2.0*np.pi**2
2027#        ast = np.sqrt(np.diag(G))
2028#        Mast = twopisq*np.multiply.outer(ast,ast)
2029#        Tdata,Mdata,Fdata,Xdata = GetAtomFX(pfx,calcControls,parmDict)
2030#        FF = np.zeros(len(Tdata))
2031#        if 'N' in parmDict[hfx+'Type']:
2032#            FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
2033#        else:
2034#            FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2035#            FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2036#        Uij = np.array(G2lat.U6toUij(Uijdata))
2037#        bij = Mast*Uij.T
2038#        for refl in refList:
2039#            fbs = np.array([0,0])
2040#            H = refl[:3]
2041#            SQ = 1./(2.*refl[4])**2
2042#            SQfactor = 4.0*SQ*twopisq
2043#            Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
2044#            if not len(refl[-1]):                #no form factors
2045#                if 'N' in parmDict[hfx+'Type']:
2046#                    refl[-1] = getBLvalues(BLtables)
2047#                else:       #'X'
2048#                    refl[-1] = getFFvalues(FFtables,SQ)
2049#            for i,El in enumerate(Tdata):
2050#                FF[i] = refl[-1][El]           
2051#            Uniq = refl[11]
2052#            phi = refl[12]
2053#            phase = twopi*(np.inner(Uniq,(Xdata.T))+phi[:,np.newaxis])
2054#            sinp = np.sin(phase)
2055#            cosp = np.cos(phase)
2056#            occ = Mdata*Fdata/len(Uniq)
2057#            fa = np.array([(FF+FP-Bab)*occ*cosp,-FPP*occ*sinp])
2058#            fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
2059#            if not SGData['SGInv']:
2060#                fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
2061#                fbs = np.sum(np.sum(fb,axis=1),axis=1)
2062#            fasq = fas**2
2063#            fbsq = fbs**2        #imaginary
2064#            refl[9] = np.sum(fasq)+np.sum(fbsq)
2065#            refl[10] = atan2d(fbs[0],fas[0])
2066   
2067    sq8ln2 = np.sqrt(8*np.log(2))
2068    sq2pi = np.sqrt(2*np.pi)
2069    sq4pi = np.sqrt(4*np.pi)
2070    generalData = data['General']
2071    SGData = generalData['SGData']
2072    fixAtoms = data['Atoms']                       #if any
2073    cx,ct,cs = generalData['AtomPtrs'][:3]
2074    aTypes = set([])
2075    parmDict = {}
2076    varyList = []
2077    atNo = 0
2078    for atm in fixAtoms:
2079        pfx = ':'+str(atNo)+':'
2080        parmDict[pfx+'Atype'] = atm[ct]
2081        aTypes |= set([atm[ct],])
2082        pstr = ['Ax','Ay','Az']
2083        parmDict[pfx+'Amul'] = atm[cs+1]
2084        for i in range(3):
2085            name = pfx+pstr[i]
2086            parmDict[name] = atm[cx+i]
2087        atNo += 1
2088    parmDict['nfixAt'] = len(fixAtoms)       
2089    mcsaControls = generalData['MCSA controls']
2090    reflName = mcsaControls['Data source']
2091    phaseName = generalData['Name']
2092    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
2093    upper = []
2094    lower = []
2095    for i,item in enumerate(MCSAObjs):
2096        mfx = str(i)+':'
2097        parmDict[mfx+'Type'] = item['Type']
2098        if item['Type'] == 'MD':
2099            getMDparms(item,mfx,parmDict,varyList)
2100        elif item['Type'] == 'Atom':
2101            pfx = mfx+str(atNo)+':'
2102            getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList)
2103            atNo += 1
2104        elif item['Type'] in ['Residue','Vector']:
2105            pfx = mfx+':'
2106            atNo = getRBparms(item,pfx,aTypes,RBdata,atNo,parmDict,varyList)
2107    parmDict['atNo'] = atNo                 #total no. of atoms
2108    parmDict['nObj'] = len(MCSAObjs)
2109    FFtables = G2el.GetFFtable(aTypes)
2110    refs = []
2111    if 'PWDR' in reflName:
2112        for ref in reflData:
2113            h,k,l,m,d,pos,sig,gam,f = ref[:9]
2114            if d >= mcsaControls['dmin']:
2115                sig = gamFW(sig,gam)/sq8ln2        #--> sig from FWHM
2116                SQ = 0.25/d**2
2117                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2118                FFs = G2el.getFFvalues(FFtables,SQ)
2119                refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi])
2120        nRef = len(refs)
2121        rcov = np.zeros((nRef,nRef))
2122        for iref,refI in enumerate(refs):
2123            rcov[iref][iref] = 1./(sq4pi*refI[6])
2124            for jref,refJ in enumerate(refs[:iref]):
2125                t1 = refI[6]**2+refJ[6]**2
2126                t2 = (refJ[5]-refI[5])**2/(2.*t1)
2127                if t2 > 10.:
2128                    rcov[iref][jref] = 0.
2129                else:
2130                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
2131        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2132        Rdiag = np.sqrt(np.diag(rcov))
2133        Rnorm = np.outer(Rdiag,Rdiag)
2134        rcov /= Rnorm
2135    elif 'Pawley' in reflName:
2136        covMatrix = covData['covMatrix']
2137        vList = covData['varyList']
2138        for iref,refI in enumerate(reflData):
2139            h,k,l,m,d,v,f,s = refI
2140            if d >= mcsaControls['dmin'] and v:       #skip unrefined ones
2141                SQ = 0.25/d**2
2142                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2143                FFs = G2el.getFFvalues(FFtables,SQ)
2144                refs.append([h,k,l,m,f*m,iref,0.,FFs,Uniq,phi])
2145        nRef = len(refs)
2146        pfx = str(data['pId'])+'::PWLref:'
2147        rcov = np.zeros((nRef,nRef))       
2148        for iref,refI in enumerate(refs):
2149            I = refI[5]
2150            nameI = pfx+str(I)
2151            if nameI in vList:
2152                Iindx = vList.index(nameI)
2153                rcov[iref][iref] = covMatrix[Iindx][Iindx]
2154                for jref,refJ in enumerate(refs[:iref]):
2155                    J = refJ[5]
2156                    nameJ = pfx+str(J)
2157                    try:
2158                        rcov[iref][jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)]
2159                    except ValueError:
2160                        rcov[iref][jref] = rcov[iref][jref-1]
2161            else:
2162                rcov[iref] = rcov[iref-1]
2163                rcov[iref][iref] = rcov[iref-1][iref-1]
2164        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2165        Rdiag = np.sqrt(np.diag(rcov))
2166        Rnorm = np.outer(Rdiag,Rdiag)
2167        rcov /= Rnorm
2168    elif 'HKLF' in reflName:
2169        for ref in reflData:
2170            [h,k,l,m,d],f = ref[:5],ref[6]
2171            if d >= mcsaControls['dmin']:
2172                SQ = 0.25/d**2
2173                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2174                FFs = G2el.getFFvalues(FFtables,SQ)
2175                refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
2176        rcov = np.identity(len(refs))
2177       
2178    for parm in parmDict:
2179        print parm,parmDict[parm] 
2180                       
2181#    XYZ,aTypes = UpdateMCSAxyz(Bmat,MCSA)       
2182           
2183#    generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50,1.e-6],
2184#    'dmin':2.0,'Algorithm':'fast','Jump coeff':[0.95,0.5],'nRuns':1,'boltzmann':1.0,
2185#    'fast parms':[1.0,1.0,1.0],'log slope':0.9}
2186
2187    return {}
2188
2189       
2190################################################################################
2191##### Quaternion stuff
2192################################################################################
2193
2194def prodQQ(QA,QB):
2195    ''' Grassman quaternion product
2196        QA,QB quaternions; q=r+ai+bj+ck
2197    '''
2198    D = np.zeros(4)
2199    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
2200    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
2201    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
2202    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
2203    return D
2204   
2205def normQ(QA):
2206    ''' get length of quaternion & normalize it
2207        q=r+ai+bj+ck
2208    '''
2209    n = np.sqrt(np.sum(np.array(QA)**2))
2210    return QA/n
2211   
2212def invQ(Q):
2213    '''
2214        get inverse of quaternion
2215        q=r+ai+bj+ck; q* = r-ai-bj-ck
2216    '''
2217    return Q*np.array([1,-1,-1,-1])
2218   
2219def prodQVQ(Q,V):
2220    """
2221    compute the quaternion vector rotation qvq-1 = v'
2222    q=r+ai+bj+ck
2223    """
2224    VP = np.zeros(3)
2225    T2 = Q[0]*Q[1]
2226    T3 = Q[0]*Q[2]
2227    T4 = Q[0]*Q[3]
2228    T5 = -Q[1]*Q[1]
2229    T6 = Q[1]*Q[2]
2230    T7 = Q[1]*Q[3]
2231    T8 = -Q[2]*Q[2]
2232    T9 = Q[2]*Q[3]
2233    T10 = -Q[3]*Q[3]
2234    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
2235    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
2236    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
2237    return VP   
2238   
2239def Q2Mat(Q):
2240    ''' make rotation matrix from quaternion
2241        q=r+ai+bj+ck
2242    '''
2243    QN = normQ(Q)
2244    aa = QN[0]**2
2245    ab = QN[0]*QN[1]
2246    ac = QN[0]*QN[2]
2247    ad = QN[0]*QN[3]
2248    bb = QN[1]**2
2249    bc = QN[1]*QN[2]
2250    bd = QN[1]*QN[3]
2251    cc = QN[2]**2
2252    cd = QN[2]*QN[3]
2253    dd = QN[3]**2
2254    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
2255        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
2256        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
2257    return np.array(M)
2258   
2259def AV2Q(A,V):
2260    ''' convert angle (radians) & vector to quaternion
2261        q=r+ai+bj+ck
2262    '''
2263    Q = np.zeros(4)
2264    d = np.sqrt(np.sum(np.array(V)**2))
2265    if d:
2266        V /= d
2267    else:
2268        V = np.array([0.,0.,1.])
2269    p = A/2.
2270    Q[0] = np.cos(p)
2271    Q[1:4] = V*np.sin(p)
2272    return Q
2273   
2274def AVdeg2Q(A,V):
2275    ''' convert angle (degrees) & vector to quaternion
2276        q=r+ai+bj+ck
2277    '''
2278    Q = np.zeros(4)
2279    d = np.sqrt(np.sum(np.array(V)**2))
2280    if d:
2281        V /= d
2282    else:
2283        V = np.array([0.,0.,1.])
2284    p = A/2.
2285    Q[0] = cosd(p)
2286    Q[1:4] = V*sind(p)
2287    return Q
2288   
2289def Q2AVdeg(Q):
2290    ''' convert quaternion to angle (degrees 0-360) & normalized vector
2291        q=r+ai+bj+ck
2292    '''
2293    A = 2.*acosd(Q[0])
2294    V = np.array(Q[1:])
2295    d = np.sqrt(np.sum(V**2))
2296    if d:
2297        V /= d
2298    else:
2299        A = 0.
2300        V = np.array([0.,0.,0.])
2301    return A,V
2302   
2303def Q2AV(Q):
2304    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
2305        q=r+ai+bj+ck
2306    '''
2307    A = 2.*np.arccos(Q[0])
2308    V = np.array(Q[1:])
2309    d = np.sqrt(np.sum(V**2))
2310    if d:
2311        V /= d
2312    else:
2313        A = 0.
2314        V = np.array([0.,0.,0.])
2315    return A,V
2316   
2317def makeQuat(A,B,C):
2318    ''' Make quaternion from rotation of A vector to B vector about C axis
2319
2320        :param np.array A,B,C: Cartesian 3-vectors
2321        :Returns: quaternion & rotation angle in radians q=r+ai+bj+ck
2322    '''
2323
2324    V1 = np.cross(A,C)
2325    V2 = np.cross(B,C)
2326    if nl.norm(V1)*nl.norm(V2):
2327        V1 /= nl.norm(V1)
2328        V2 /= nl.norm(V2)
2329        V3 = np.cross(V1,V2)
2330    else:
2331        V3 = np.zeros(3)
2332    Q = np.array([0.,0.,0.,1.])
2333    D = 0.
2334    if nl.norm(V3):
2335        V3 /= nl.norm(V3)
2336        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
2337        D = np.arccos(D1)/2.0
2338        V1 = C-V3
2339        V2 = C+V3
2340        DM = nl.norm(V1)
2341        DP = nl.norm(V2)
2342        S = np.sin(D)
2343        Q[0] = np.cos(D)
2344        Q[1:] = V3*S
2345        D *= 2.
2346        if DM > DP:
2347            D *= -1.
2348    return Q,D
2349   
2350if __name__ == "__main__":
2351    from numpy import cos
2352    # minimum expected at ~-0.195
2353    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
2354    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
2355    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
2356    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
2357
2358    # minimum expected at ~[-0.195, -0.1]
2359    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
2360    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')
2361    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')
2362    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')
Note: See TracBrowser for help on using the repository browser.