source: trunk/GSASIImath.py @ 945

Last change on this file since 945 was 945, checked in by vondreele, 9 years ago

some mods for POWGEN data for single peak fitting
unify profile parm calcs into G2mth

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 80.9 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2013-06-11 14:58:44 +0000 (Tue, 11 Jun 2013) $
5# $Author: vondreele $
6# $Revision: 945 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 945 2013-06-11 14:58:44Z 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: 945 $")
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    Parameters
1739    ----------
1740    func : callable f(x, *args)
1741        Function to be optimized.
1742    x0 : ndarray
1743        Initial guess.
1744    args : tuple
1745        Extra parameters to `func`.
1746    schedule : base_schedule
1747        Annealing schedule to use (a class).
1748    full_output : bool
1749        Whether to return optional outputs.
1750    T0 : float
1751        Initial Temperature (estimated as 1.2 times the largest
1752        cost-function deviation over random points in the range).
1753    Tf : float
1754        Final goal temperature.
1755    maxeval : int
1756        Maximum function evaluations.
1757    maxaccept : int
1758        Maximum changes to accept.
1759    maxiter : int
1760        Maximum cooling iterations.
1761    learn_rate : float
1762        Scale constant for adjusting guesses.
1763    boltzmann : float
1764        Boltzmann constant in acceptance test
1765        (increase for less stringent test at each temperature).
1766    feps : float
1767        Stopping relative error tolerance for the function value in
1768        last four coolings.
1769    quench, m, n : float
1770        Parameters to alter fast_sa schedule.
1771    lower, upper : float or ndarray
1772        Lower and upper bounds on `x`.
1773    dwell : int
1774        The number of times to search the space at each temperature.
1775    slope : float
1776        Parameter for log schedule
1777
1778    Returns
1779    -------
1780    xmin : ndarray
1781        Point giving smallest value found.
1782    Jmin : float
1783        Minimum value of function found.
1784    T : float
1785        Final temperature.
1786    feval : int
1787        Number of function evaluations.
1788    iters : int
1789        Number of cooling iterations.
1790    accept : int
1791        Number of tests accepted.
1792    retval : int
1793        Flag indicating stopping condition::
1794
1795                0 : Points no longer changing
1796                1 : Cooled to final temperature
1797                2 : Maximum function evaluations
1798                3 : Maximum cooling iterations reached
1799                4 : Maximum accepted query locations reached
1800                5 : Final point not the minimum amongst encountered points
1801
1802    Notes
1803    -----
1804    Simulated annealing is a random algorithm which uses no derivative
1805    information from the function being optimized. In practice it has
1806    been more useful in discrete optimization than continuous
1807    optimization, as there are usually better algorithms for continuous
1808    optimization problems.
1809
1810    Some experimentation by trying the difference temperature
1811    schedules and altering their parameters is likely required to
1812    obtain good performance.
1813
1814    The randomness in the algorithm comes from random sampling in numpy.
1815    To obtain the same results you can call numpy.random.seed with the
1816    same seed immediately before calling scipy.optimize.anneal.
1817
1818    We give a brief description of how the three temperature schedules
1819    generate new points and vary their temperature. Temperatures are
1820    only updated with iterations in the outer loop. The inner loop is
1821    over loop over xrange(dwell), and new points are generated for
1822    every iteration in the inner loop. (Though whether the proposed
1823    new points are accepted is probabilistic.)
1824
1825    For readability, let d denote the dimension of the inputs to func.
1826    Also, let x_old denote the previous state, and k denote the
1827    iteration number of the outer loop. All other variables not
1828    defined below are input variables to scipy.optimize.anneal itself.
1829
1830    In the 'fast' schedule the updates are ::
1831
1832        u ~ Uniform(0, 1, size=d)
1833        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
1834        xc = y * (upper - lower)
1835        x_new = x_old + xc
1836
1837        c = n * exp(-n * quench)
1838        T_new = T0 * exp(-c * k**quench)
1839
1840
1841    In the 'cauchy' schedule the updates are ::
1842
1843        u ~ Uniform(-pi/2, pi/2, size=d)
1844        xc = learn_rate * T * tan(u)
1845        x_new = x_old + xc
1846
1847        T_new = T0 / (1+k)
1848
1849    In the 'boltzmann' schedule the updates are ::
1850
1851        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
1852        y ~ Normal(0, std, size=d)
1853        x_new = x_old + learn_rate * y
1854
1855        T_new = T0 / log(1+k)
1856
1857    """
1858    x0 = asarray(x0)
1859    lower = asarray(lower)
1860    upper = asarray(upper)
1861
1862    schedule = eval(schedule+'_sa()')
1863    #   initialize the schedule
1864    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
1865                  learn_rate=learn_rate, lower=lower, upper=upper,
1866                  m=m, n=n, quench=quench, dwell=dwell)
1867
1868    current_state, last_state, best_state = _state(), _state(), _state()
1869    if T0 is None:
1870        x0 = schedule.getstart_temp(best_state)
1871    else:
1872        best_state.x = None
1873        best_state.cost = numpy.Inf
1874
1875    last_state.x = asarray(x0).copy()
1876    fval = func(x0,*args)
1877    schedule.feval += 1
1878    last_state.cost = fval
1879    if last_state.cost < best_state.cost:
1880        best_state.cost = fval
1881        best_state.x = asarray(x0).copy()
1882    schedule.T = schedule.T0
1883    fqueue = [100, 300, 500, 700]
1884    iters = 0
1885    while 1:
1886        for n in xrange(dwell):
1887            current_state.x = schedule.update_guess(last_state.x)
1888            current_state.cost = func(current_state.x,*args)
1889            schedule.feval += 1
1890
1891            dE = current_state.cost - last_state.cost
1892            if schedule.accept_test(dE):
1893                last_state.x = current_state.x.copy()
1894                last_state.cost = current_state.cost
1895                if last_state.cost < best_state.cost:
1896                    best_state.x = last_state.x.copy()
1897                    best_state.cost = last_state.cost
1898        schedule.update_temp()
1899        iters += 1
1900        # Stopping conditions
1901        # 0) last saved values of f from each cooling step
1902        #     are all very similar (effectively cooled)
1903        # 1) Tf is set and we are below it
1904        # 2) maxeval is set and we are past it
1905        # 3) maxiter is set and we are past it
1906        # 4) maxaccept is set and we are past it
1907
1908        fqueue.append(squeeze(last_state.cost))
1909        fqueue.pop(0)
1910        af = asarray(fqueue)*1.0
1911        if all(abs((af-af[0])/af[0]) < feps):
1912            retval = 0
1913            if abs(af[-1]-best_state.cost) > feps*10:
1914                retval = 5
1915                print "Warning: Cooled to %f at %s but this is not" \
1916                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
1917                      + " the smallest point found."
1918            break
1919        if (Tf is not None) and (schedule.T < Tf):
1920            retval = 1
1921            break
1922        if (maxeval is not None) and (schedule.feval > maxeval):
1923            retval = 2
1924            break
1925        if (iters > maxiter):
1926            print "Warning: Maximum number of iterations exceeded."
1927            retval = 3
1928            break
1929        if (maxaccept is not None) and (schedule.accepted > maxaccept):
1930            retval = 4
1931            break
1932
1933    if full_output:
1934        return best_state.x, best_state.cost, schedule.T, \
1935               schedule.feval, iters, schedule.accepted, retval
1936    else:
1937        return best_state.x, retval
1938
1939
1940def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
1941    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.)
1942   
1943    def getMDparms(item,pfx,parmDict,varyList):
1944        parmDict[pfx+'MDaxis'] = item['axis']
1945        parmDict[pfx+'MDval'] = item['Coef'][0]
1946        if item['Coef'][1]:
1947            varyList += [pfx+'MDval',]
1948            limits = item['Coef'][2]
1949            lower.append(limits[0])
1950            upper.append(limits[1])
1951           
1952    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
1953        parmDict[pfx+'Atype'] = item['atType']
1954        aTypes |= set([item['atType'],]) 
1955        pstr = ['Ax','Ay','Az']
1956        XYZ = [0,0,0]
1957        for i in range(3):
1958            name = pfx+pstr[i]
1959            parmDict[name] = item['Pos'][0][i]
1960            XYZ[i] = parmDict[name]
1961            if item['Pos'][1][i]:
1962                varyList += [name,]
1963                limits = item['Pos'][2][i]
1964                lower.append(limits[0])
1965                upper.append(limits[1])
1966        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
1967           
1968    def getRBparms(item,mfx,aTypes,RBdata,atNo,parmDict,varyList):
1969        parmDict[mfx+'MolCent'] = item['MolCent']
1970        parmDict[mfx+'RBId'] = item['RBId']
1971        pstr = ['Px','Py','Pz']
1972        ostr = ['Qa','Qi','Qj','Qk']
1973        for i in range(3):
1974            name = mfx+pstr[i]
1975            parmDict[name] = item['Pos'][0][i]
1976            if item['Pos'][1][i]:
1977                varyList += [name,]
1978                limits = item['Pos'][2][i]
1979                lower.append(limits[0])
1980                upper.append(limits[1])
1981        for i in range(4):
1982            name = mfx+ostr[i]
1983            parmDict[name] = item['Ori'][0][i]
1984            if item['Ovar'] == 'AV' and i:
1985                varyList += [name,]
1986                limits = item['Ori'][2][i]
1987                lower.append(limits[0])
1988                upper.append(limits[1])
1989            elif item['Ovar'] == 'A' and not i:
1990                varyList += [name,]
1991                limits = item['Ori'][2][i]
1992                lower.append(limits[0])
1993                upper.append(limits[1])
1994        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
1995            for i in range(len(item['Tor'][0])):
1996                name = mfx+'Tor'+str(i)
1997                parmDict[name] = item['Tor'][0][i]
1998                if item['Tor'][1][i]:
1999                    varyList += [name,]
2000                    limits = item['Tor'][2][i]
2001                    lower.append(limits[0])
2002                    upper.append(limits[1])
2003        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2004        aTypes |= set(atypes)
2005        atNo += len(atypes)
2006        return atNo
2007
2008    def GetAtomTMX(pfx,RBdata,parmDict):
2009        'Needs a doc string'
2010        atNo = parmDIct['atNo']
2011        nfixAt = parmDict['nfixAt']
2012        Tdata = atNo*[' ',]
2013        Mdata = np.zeros(atNo)
2014        Fdata = np.zeros(atNo)
2015        Xdata = np.zeros((3,atNo))
2016        keys = {'Atype:':Tdata,'Amul:':Mdata,
2017            'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2]}
2018        nObjs = parmDict['nObjs']
2019        for iatm in range(nfixAt):
2020            for key in keys:
2021                parm = pfx+key+str(iatm)
2022                if parm in parmDict:
2023                    keys[key][iatm] = parmDict[parm]
2024        return Tdata,Mdata,Xdata
2025   
2026#    def mcsaSfCalc(refList,RBdata,G,SGData,parmDict):
2027#        ''' Compute structure factors for all h,k,l for phase
2028#        input:
2029#            refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
2030#            G:      reciprocal metric tensor
2031#            SGData: space group info. dictionary output from SpcGroup
2032#            ParmDict:
2033#        puts result F^2 in each ref[8] in refList
2034#        '''       
2035#        twopi = 2.0*np.pi
2036#        twopisq = 2.0*np.pi**2
2037#        ast = np.sqrt(np.diag(G))
2038#        Mast = twopisq*np.multiply.outer(ast,ast)
2039#        Tdata,Mdata,Fdata,Xdata = GetAtomFX(pfx,calcControls,parmDict)
2040#        FF = np.zeros(len(Tdata))
2041#        if 'N' in parmDict[hfx+'Type']:
2042#            FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
2043#        else:
2044#            FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2045#            FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2046#        Uij = np.array(G2lat.U6toUij(Uijdata))
2047#        bij = Mast*Uij.T
2048#        for refl in refList:
2049#            fbs = np.array([0,0])
2050#            H = refl[:3]
2051#            SQ = 1./(2.*refl[4])**2
2052#            SQfactor = 4.0*SQ*twopisq
2053#            Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
2054#            if not len(refl[-1]):                #no form factors
2055#                if 'N' in parmDict[hfx+'Type']:
2056#                    refl[-1] = getBLvalues(BLtables)
2057#                else:       #'X'
2058#                    refl[-1] = getFFvalues(FFtables,SQ)
2059#            for i,El in enumerate(Tdata):
2060#                FF[i] = refl[-1][El]           
2061#            Uniq = refl[11]
2062#            phi = refl[12]
2063#            phase = twopi*(np.inner(Uniq,(Xdata.T))+phi[:,np.newaxis])
2064#            sinp = np.sin(phase)
2065#            cosp = np.cos(phase)
2066#            occ = Mdata*Fdata/len(Uniq)
2067#            fa = np.array([(FF+FP-Bab)*occ*cosp,-FPP*occ*sinp])
2068#            fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
2069#            if not SGData['SGInv']:
2070#                fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
2071#                fbs = np.sum(np.sum(fb,axis=1),axis=1)
2072#            fasq = fas**2
2073#            fbsq = fbs**2        #imaginary
2074#            refl[9] = np.sum(fasq)+np.sum(fbsq)
2075#            refl[10] = atan2d(fbs[0],fas[0])
2076   
2077    sq8ln2 = np.sqrt(8*np.log(2))
2078    sq2pi = np.sqrt(2*np.pi)
2079    sq4pi = np.sqrt(4*np.pi)
2080    generalData = data['General']
2081    SGData = generalData['SGData']
2082    fixAtoms = data['Atoms']                       #if any
2083    cx,ct,cs = generalData['AtomPtrs'][:3]
2084    aTypes = set([])
2085    parmDict = {}
2086    varyList = []
2087    atNo = 0
2088    for atm in fixAtoms:
2089        pfx = ':'+str(atNo)+':'
2090        parmDict[pfx+'Atype'] = atm[ct]
2091        aTypes |= set([atm[ct],])
2092        pstr = ['Ax','Ay','Az']
2093        parmDict[pfx+'Amul'] = atm[cs+1]
2094        for i in range(3):
2095            name = pfx+pstr[i]
2096            parmDict[name] = atm[cx+i]
2097        atNo += 1
2098    parmDict['nfixAt'] = len(fixAtoms)       
2099    mcsaControls = generalData['MCSA controls']
2100    reflName = mcsaControls['Data source']
2101    phaseName = generalData['Name']
2102    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
2103    upper = []
2104    lower = []
2105    for i,item in enumerate(MCSAObjs):
2106        mfx = str(i)+':'
2107        parmDict[mfx+'Type'] = item['Type']
2108        if item['Type'] == 'MD':
2109            getMDparms(item,mfx,parmDict,varyList)
2110        elif item['Type'] == 'Atom':
2111            pfx = mfx+str(atNo)+':'
2112            getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList)
2113            atNo += 1
2114        elif item['Type'] in ['Residue','Vector']:
2115            pfx = mfx+':'
2116            atNo = getRBparms(item,pfx,aTypes,RBdata,atNo,parmDict,varyList)
2117    parmDict['atNo'] = atNo                 #total no. of atoms
2118    parmDict['nObj'] = len(MCSAObjs)
2119    FFtables = G2el.GetFFtable(aTypes)
2120    refs = []
2121    if 'PWDR' in reflName:
2122        for ref in reflData:
2123            h,k,l,m,d,pos,sig,gam,f = ref[:9]
2124            if d >= mcsaControls['dmin']:
2125                sig = gamFW(sig,gam)/sq8ln2        #--> sig from FWHM
2126                SQ = 0.25/d**2
2127                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2128                FFs = G2el.getFFvalues(FFtables,SQ)
2129                refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi])
2130        nRef = len(refs)
2131        rcov = np.zeros((nRef,nRef))
2132        for iref,refI in enumerate(refs):
2133            rcov[iref][iref] = 1./(sq4pi*refI[6])
2134            for jref,refJ in enumerate(refs[:iref]):
2135                t1 = refI[6]**2+refJ[6]**2
2136                t2 = (refJ[5]-refI[5])**2/(2.*t1)
2137                if t2 > 10.:
2138                    rcov[iref][jref] = 0.
2139                else:
2140                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
2141        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2142        Rdiag = np.sqrt(np.diag(rcov))
2143        Rnorm = np.outer(Rdiag,Rdiag)
2144        rcov /= Rnorm
2145    elif 'Pawley' in reflName:
2146        covMatrix = covData['covMatrix']
2147        vList = covData['varyList']
2148        for iref,refI in enumerate(reflData):
2149            h,k,l,m,d,v,f,s = refI
2150            if d >= mcsaControls['dmin'] and v:       #skip unrefined ones
2151                SQ = 0.25/d**2
2152                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2153                FFs = G2el.getFFvalues(FFtables,SQ)
2154                refs.append([h,k,l,m,f*m,iref,0.,FFs,Uniq,phi])
2155        nRef = len(refs)
2156        pfx = str(data['pId'])+'::PWLref:'
2157        rcov = np.zeros((nRef,nRef))       
2158        for iref,refI in enumerate(refs):
2159            I = refI[5]
2160            nameI = pfx+str(I)
2161            if nameI in vList:
2162                Iindx = vList.index(nameI)
2163                rcov[iref][iref] = covMatrix[Iindx][Iindx]
2164                for jref,refJ in enumerate(refs[:iref]):
2165                    J = refJ[5]
2166                    nameJ = pfx+str(J)
2167                    try:
2168                        rcov[iref][jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)]
2169                    except ValueError:
2170                        rcov[iref][jref] = rcov[iref][jref-1]
2171            else:
2172                rcov[iref] = rcov[iref-1]
2173                rcov[iref][iref] = rcov[iref-1][iref-1]
2174        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2175        Rdiag = np.sqrt(np.diag(rcov))
2176        Rnorm = np.outer(Rdiag,Rdiag)
2177        rcov /= Rnorm
2178    elif 'HKLF' in reflName:
2179        for ref in reflData:
2180            [h,k,l,m,d],f = ref[:5],ref[6]
2181            if d >= mcsaControls['dmin']:
2182                SQ = 0.25/d**2
2183                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2184                FFs = G2el.getFFvalues(FFtables,SQ)
2185                refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
2186        rcov = np.identity(len(refs))
2187       
2188    for parm in parmDict:
2189        print parm,parmDict[parm] 
2190                       
2191#    XYZ,aTypes = UpdateMCSAxyz(Bmat,MCSA)       
2192           
2193#    generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50,1.e-6],
2194#    'dmin':2.0,'Algorithm':'fast','Jump coeff':[0.95,0.5],'nRuns':1,'boltzmann':1.0,
2195#    'fast parms':[1.0,1.0,1.0],'log slope':0.9}
2196
2197    return {}
2198
2199       
2200################################################################################
2201##### Quaternion stuff
2202################################################################################
2203
2204def prodQQ(QA,QB):
2205    ''' Grassman quaternion product
2206        QA,QB quaternions; q=r+ai+bj+ck
2207    '''
2208    D = np.zeros(4)
2209    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
2210    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
2211    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
2212    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
2213    return D
2214   
2215def normQ(QA):
2216    ''' get length of quaternion & normalize it
2217        q=r+ai+bj+ck
2218    '''
2219    n = np.sqrt(np.sum(np.array(QA)**2))
2220    return QA/n
2221   
2222def invQ(Q):
2223    '''
2224        get inverse of quaternion
2225        q=r+ai+bj+ck; q* = r-ai-bj-ck
2226    '''
2227    return Q*np.array([1,-1,-1,-1])
2228   
2229def prodQVQ(Q,V):
2230    """
2231    compute the quaternion vector rotation qvq-1 = v'
2232    q=r+ai+bj+ck
2233    """
2234    VP = np.zeros(3)
2235    T2 = Q[0]*Q[1]
2236    T3 = Q[0]*Q[2]
2237    T4 = Q[0]*Q[3]
2238    T5 = -Q[1]*Q[1]
2239    T6 = Q[1]*Q[2]
2240    T7 = Q[1]*Q[3]
2241    T8 = -Q[2]*Q[2]
2242    T9 = Q[2]*Q[3]
2243    T10 = -Q[3]*Q[3]
2244    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
2245    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
2246    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
2247    return VP   
2248   
2249def Q2Mat(Q):
2250    ''' make rotation matrix from quaternion
2251        q=r+ai+bj+ck
2252    '''
2253    QN = normQ(Q)
2254    aa = QN[0]**2
2255    ab = QN[0]*QN[1]
2256    ac = QN[0]*QN[2]
2257    ad = QN[0]*QN[3]
2258    bb = QN[1]**2
2259    bc = QN[1]*QN[2]
2260    bd = QN[1]*QN[3]
2261    cc = QN[2]**2
2262    cd = QN[2]*QN[3]
2263    dd = QN[3]**2
2264    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
2265        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
2266        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
2267    return np.array(M)
2268   
2269def AV2Q(A,V):
2270    ''' convert angle (radians) & vector to quaternion
2271        q=r+ai+bj+ck
2272    '''
2273    Q = np.zeros(4)
2274    d = np.sqrt(np.sum(np.array(V)**2))
2275    if d:
2276        V /= d
2277    else:
2278        V = np.array([0.,0.,1.])
2279    p = A/2.
2280    Q[0] = np.cos(p)
2281    Q[1:4] = V*np.sin(p)
2282    return Q
2283   
2284def AVdeg2Q(A,V):
2285    ''' convert angle (degrees) & vector to quaternion
2286        q=r+ai+bj+ck
2287    '''
2288    Q = np.zeros(4)
2289    d = np.sqrt(np.sum(np.array(V)**2))
2290    if d:
2291        V /= d
2292    else:
2293        V = np.array([0.,0.,1.])
2294    p = A/2.
2295    Q[0] = cosd(p)
2296    Q[1:4] = V*sind(p)
2297    return Q
2298   
2299def Q2AVdeg(Q):
2300    ''' convert quaternion to angle (degrees 0-360) & normalized vector
2301        q=r+ai+bj+ck
2302    '''
2303    A = 2.*acosd(Q[0])
2304    V = np.array(Q[1:])
2305    d = np.sqrt(np.sum(V**2))
2306    if d:
2307        V /= d
2308    else:
2309        A = 0.
2310        V = np.array([0.,0.,0.])
2311    return A,V
2312   
2313def Q2AV(Q):
2314    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
2315        q=r+ai+bj+ck
2316    '''
2317    A = 2.*np.arccos(Q[0])
2318    V = np.array(Q[1:])
2319    d = np.sqrt(np.sum(V**2))
2320    if d:
2321        V /= d
2322    else:
2323        A = 0.
2324        V = np.array([0.,0.,0.])
2325    return A,V
2326   
2327def makeQuat(A,B,C):
2328    ''' Make quaternion from rotation of A vector to B vector about C axis
2329
2330        :param np.array A,B,C: Cartesian 3-vectors
2331        :Returns: quaternion & rotation angle in radians q=r+ai+bj+ck
2332    '''
2333
2334    V1 = np.cross(A,C)
2335    V2 = np.cross(B,C)
2336    if nl.norm(V1)*nl.norm(V2):
2337        V1 /= nl.norm(V1)
2338        V2 /= nl.norm(V2)
2339        V3 = np.cross(V1,V2)
2340    else:
2341        V3 = np.zeros(3)
2342    Q = np.array([0.,0.,0.,1.])
2343    D = 0.
2344    if nl.norm(V3):
2345        V3 /= nl.norm(V3)
2346        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
2347        D = np.arccos(D1)/2.0
2348        V1 = C-V3
2349        V2 = C+V3
2350        DM = nl.norm(V1)
2351        DP = nl.norm(V2)
2352        S = np.sin(D)
2353        Q[0] = np.cos(D)
2354        Q[1:] = V3*S
2355        D *= 2.
2356        if DM > DP:
2357            D *= -1.
2358    return Q,D
2359   
2360if __name__ == "__main__":
2361    from numpy import cos
2362    # minimum expected at ~-0.195
2363    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
2364    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
2365    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
2366    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
2367
2368    # minimum expected at ~[-0.195, -0.1]
2369    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
2370    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')
2371    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')
2372    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.