source: trunk/GSASIImath.py @ 943

Last change on this file since 943 was 943, checked in by vondreele, 10 years ago

more work on MC/SA
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 80.1 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2013-06-05 19:38:00 +0000 (Wed, 05 Jun 2013) $
5# $Author: vondreele $
6# $Revision: 943 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 943 2013-06-05 19:38:00Z 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: 943 $")
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 setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
1468    'Needs a doc string'       
1469    ind = 0
1470    if useFit:
1471        ind = 1
1472    ins = {}
1473    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
1474        for x in ['U','V','W','X','Y']:
1475            ins[x] = Parms[x][ind]
1476        if ifQ:                              #qplot - convert back to 2-theta
1477            pos = 2.0*asind(pos*wave/(4*math.pi))
1478        sig = ins['U']*tand(pos/2.0)**2+ins['V']*tand(pos/2.0)+ins['W']
1479        gam = ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)           
1480        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
1481    else:
1482        if ifQ:
1483            dsp = 2.*np.pi/pos
1484            pos = Parms['difC']*dsp
1485        else:
1486            dsp = pos/Parms['difC'][1]
1487        if 'Pdabc' in Parms2:
1488            for x in ['sig-0','sig-1','X','Y']:
1489                ins[x] = Parms[x][ind]
1490            Pdabc = Parms2['Pdabc'].T
1491            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
1492            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
1493        else:
1494            for x in ['alpha','beta-0','beta-1','sig-0','sig-1','X','Y']:
1495                ins[x] = Parms[x][ind]
1496            alp = ins['alpha']/dsp
1497            bet = ins['beta-0']+ins['beta-1']/dsp**4
1498        sig = ins['sig-0']+ins['sig-1']*dsp**2
1499        gam = ins['X']*dsp+ins['Y']*dsp**2
1500        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
1501    return XY
1502   
1503################################################################################
1504##### MC/SA stuff
1505################################################################################
1506
1507#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
1508# Original Author: Travis Oliphant 2002
1509# Bug-fixes in 2006 by Tim Leslie
1510
1511
1512import numpy
1513from numpy import asarray, tan, exp, ones, squeeze, sign, \
1514     all, log, sqrt, pi, shape, array, minimum, where
1515from numpy import random
1516
1517__all__ = ['anneal']
1518
1519_double_min = numpy.finfo(float).min
1520_double_max = numpy.finfo(float).max
1521class base_schedule(object):
1522    def __init__(self):
1523        self.dwell = 20
1524        self.learn_rate = 0.5
1525        self.lower = -10
1526        self.upper = 10
1527        self.Ninit = 50
1528        self.accepted = 0
1529        self.tests = 0
1530        self.feval = 0
1531        self.k = 0
1532        self.T = None
1533
1534    def init(self, **options):
1535        self.__dict__.update(options)
1536        self.lower = asarray(self.lower)
1537        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
1538        self.upper = asarray(self.upper)
1539        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
1540        self.k = 0
1541        self.accepted = 0
1542        self.feval = 0
1543        self.tests = 0
1544
1545    def getstart_temp(self, best_state):
1546        """ Find a matching starting temperature and starting parameters vector
1547        i.e. find x0 such that func(x0) = T0.
1548
1549        Parameters
1550        ----------
1551        best_state : _state
1552            A _state object to store the function value and x0 found.
1553
1554        Returns
1555        -------
1556        x0 : array
1557            The starting parameters vector.
1558        """
1559
1560        assert(not self.dims is None)
1561        lrange = self.lower
1562        urange = self.upper
1563        fmax = _double_min
1564        fmin = _double_max
1565        for _ in range(self.Ninit):
1566            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
1567            fval = self.func(x0, *self.args)
1568            self.feval += 1
1569            if fval > fmax:
1570                fmax = fval
1571            if fval < fmin:
1572                fmin = fval
1573                best_state.cost = fval
1574                best_state.x = array(x0)
1575
1576        self.T0 = (fmax-fmin)*1.5
1577        return best_state.x
1578
1579    def accept_test(self, dE):
1580        T = self.T
1581        self.tests += 1
1582        if dE < 0:
1583            self.accepted += 1
1584            return 1
1585        p = exp(-dE*1.0/self.boltzmann/T)
1586        if (p > random.uniform(0.0, 1.0)):
1587            self.accepted += 1
1588            return 1
1589        return 0
1590
1591    def update_guess(self, x0):
1592        pass
1593
1594    def update_temp(self, x0):
1595        pass
1596
1597
1598#  A schedule due to Lester Ingber
1599class fast_sa(base_schedule):
1600    def init(self, **options):
1601        self.__dict__.update(options)
1602        if self.m is None:
1603            self.m = 1.0
1604        if self.n is None:
1605            self.n = 1.0
1606        self.c = self.m * exp(-self.n * self.quench)
1607
1608    def update_guess(self, x0):
1609        x0 = asarray(x0)
1610        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
1611        T = self.T
1612        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
1613        xc = y*(self.upper - self.lower)
1614        xnew = x0 + xc
1615        return xnew
1616
1617    def update_temp(self):
1618        self.T = self.T0*exp(-self.c * self.k**(self.quench))
1619        self.k += 1
1620        return
1621
1622class cauchy_sa(base_schedule):
1623    def update_guess(self, x0):
1624        x0 = asarray(x0)
1625        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
1626        xc = self.learn_rate * self.T * tan(numbers)
1627        xnew = x0 + xc
1628        return xnew
1629
1630    def update_temp(self):
1631        self.T = self.T0/(1+self.k)
1632        self.k += 1
1633        return
1634
1635class boltzmann_sa(base_schedule):
1636    def update_guess(self, x0):
1637        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
1638        x0 = asarray(x0)
1639        xc = squeeze(random.normal(0, 1.0, size=self.dims))
1640
1641        xnew = x0 + xc*std*self.learn_rate
1642        return xnew
1643
1644    def update_temp(self):
1645        self.k += 1
1646        self.T = self.T0 / log(self.k+1.0)
1647        return
1648
1649class log_sa(base_schedule):
1650
1651    def init(self,**options):
1652        self.__dict__.update(options)
1653       
1654    def update_guess(self,x0):
1655        x0 = np.asarray(x0)
1656        u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
1657        xnew = x0+u
1658        return xnew
1659       
1660    def update_temp(self):
1661        self.k += 1
1662        self.T = self.T0*self.slope**k
1663       
1664class Tremayne_sa(base_schedule):   #needs fixing for two steps
1665
1666    def init(self,**options):
1667        self.__dict__.update(options)
1668
1669    def update_guess(self,x0):
1670        x0 = np.asarray(x0)
1671        u = np.squeeze(np.random.uniform(0.,1.,size=self.dims))
1672        xnew = x0+u
1673        return xnew       
1674   
1675    def update_temp(self):
1676        self.k += 1
1677        self.T = self.T0*self.slope**k
1678   
1679class _state(object):
1680    def __init__(self):
1681        self.x = None
1682        self.cost = None
1683
1684# TODO:
1685#     allow for general annealing temperature profile
1686#     in that case use update given by alpha and omega and
1687#     variation of all previous updates and temperature?
1688
1689# Simulated annealing
1690
1691def anneal(func, x0, args=(), schedule='fast', full_output=0,
1692           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
1693           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
1694           lower=-100, upper=100, dwell=50, slope=0.9):
1695    """Minimize a function using simulated annealing.
1696
1697    Schedule is a schedule class implementing the annealing schedule.
1698    Available ones are 'fast', 'cauchy', 'boltzmann'
1699
1700    Parameters
1701    ----------
1702    func : callable f(x, *args)
1703        Function to be optimized.
1704    x0 : ndarray
1705        Initial guess.
1706    args : tuple
1707        Extra parameters to `func`.
1708    schedule : base_schedule
1709        Annealing schedule to use (a class).
1710    full_output : bool
1711        Whether to return optional outputs.
1712    T0 : float
1713        Initial Temperature (estimated as 1.2 times the largest
1714        cost-function deviation over random points in the range).
1715    Tf : float
1716        Final goal temperature.
1717    maxeval : int
1718        Maximum function evaluations.
1719    maxaccept : int
1720        Maximum changes to accept.
1721    maxiter : int
1722        Maximum cooling iterations.
1723    learn_rate : float
1724        Scale constant for adjusting guesses.
1725    boltzmann : float
1726        Boltzmann constant in acceptance test
1727        (increase for less stringent test at each temperature).
1728    feps : float
1729        Stopping relative error tolerance for the function value in
1730        last four coolings.
1731    quench, m, n : float
1732        Parameters to alter fast_sa schedule.
1733    lower, upper : float or ndarray
1734        Lower and upper bounds on `x`.
1735    dwell : int
1736        The number of times to search the space at each temperature.
1737    slope : float
1738        Parameter for log schedule
1739
1740    Returns
1741    -------
1742    xmin : ndarray
1743        Point giving smallest value found.
1744    Jmin : float
1745        Minimum value of function found.
1746    T : float
1747        Final temperature.
1748    feval : int
1749        Number of function evaluations.
1750    iters : int
1751        Number of cooling iterations.
1752    accept : int
1753        Number of tests accepted.
1754    retval : int
1755        Flag indicating stopping condition::
1756
1757                0 : Points no longer changing
1758                1 : Cooled to final temperature
1759                2 : Maximum function evaluations
1760                3 : Maximum cooling iterations reached
1761                4 : Maximum accepted query locations reached
1762                5 : Final point not the minimum amongst encountered points
1763
1764    Notes
1765    -----
1766    Simulated annealing is a random algorithm which uses no derivative
1767    information from the function being optimized. In practice it has
1768    been more useful in discrete optimization than continuous
1769    optimization, as there are usually better algorithms for continuous
1770    optimization problems.
1771
1772    Some experimentation by trying the difference temperature
1773    schedules and altering their parameters is likely required to
1774    obtain good performance.
1775
1776    The randomness in the algorithm comes from random sampling in numpy.
1777    To obtain the same results you can call numpy.random.seed with the
1778    same seed immediately before calling scipy.optimize.anneal.
1779
1780    We give a brief description of how the three temperature schedules
1781    generate new points and vary their temperature. Temperatures are
1782    only updated with iterations in the outer loop. The inner loop is
1783    over loop over xrange(dwell), and new points are generated for
1784    every iteration in the inner loop. (Though whether the proposed
1785    new points are accepted is probabilistic.)
1786
1787    For readability, let d denote the dimension of the inputs to func.
1788    Also, let x_old denote the previous state, and k denote the
1789    iteration number of the outer loop. All other variables not
1790    defined below are input variables to scipy.optimize.anneal itself.
1791
1792    In the 'fast' schedule the updates are ::
1793
1794        u ~ Uniform(0, 1, size=d)
1795        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
1796        xc = y * (upper - lower)
1797        x_new = x_old + xc
1798
1799        c = n * exp(-n * quench)
1800        T_new = T0 * exp(-c * k**quench)
1801
1802
1803    In the 'cauchy' schedule the updates are ::
1804
1805        u ~ Uniform(-pi/2, pi/2, size=d)
1806        xc = learn_rate * T * tan(u)
1807        x_new = x_old + xc
1808
1809        T_new = T0 / (1+k)
1810
1811    In the 'boltzmann' schedule the updates are ::
1812
1813        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
1814        y ~ Normal(0, std, size=d)
1815        x_new = x_old + learn_rate * y
1816
1817        T_new = T0 / log(1+k)
1818
1819    """
1820    x0 = asarray(x0)
1821    lower = asarray(lower)
1822    upper = asarray(upper)
1823
1824    schedule = eval(schedule+'_sa()')
1825    #   initialize the schedule
1826    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
1827                  learn_rate=learn_rate, lower=lower, upper=upper,
1828                  m=m, n=n, quench=quench, dwell=dwell)
1829
1830    current_state, last_state, best_state = _state(), _state(), _state()
1831    if T0 is None:
1832        x0 = schedule.getstart_temp(best_state)
1833    else:
1834        best_state.x = None
1835        best_state.cost = numpy.Inf
1836
1837    last_state.x = asarray(x0).copy()
1838    fval = func(x0,*args)
1839    schedule.feval += 1
1840    last_state.cost = fval
1841    if last_state.cost < best_state.cost:
1842        best_state.cost = fval
1843        best_state.x = asarray(x0).copy()
1844    schedule.T = schedule.T0
1845    fqueue = [100, 300, 500, 700]
1846    iters = 0
1847    while 1:
1848        for n in xrange(dwell):
1849            current_state.x = schedule.update_guess(last_state.x)
1850            current_state.cost = func(current_state.x,*args)
1851            schedule.feval += 1
1852
1853            dE = current_state.cost - last_state.cost
1854            if schedule.accept_test(dE):
1855                last_state.x = current_state.x.copy()
1856                last_state.cost = current_state.cost
1857                if last_state.cost < best_state.cost:
1858                    best_state.x = last_state.x.copy()
1859                    best_state.cost = last_state.cost
1860        schedule.update_temp()
1861        iters += 1
1862        # Stopping conditions
1863        # 0) last saved values of f from each cooling step
1864        #     are all very similar (effectively cooled)
1865        # 1) Tf is set and we are below it
1866        # 2) maxeval is set and we are past it
1867        # 3) maxiter is set and we are past it
1868        # 4) maxaccept is set and we are past it
1869
1870        fqueue.append(squeeze(last_state.cost))
1871        fqueue.pop(0)
1872        af = asarray(fqueue)*1.0
1873        if all(abs((af-af[0])/af[0]) < feps):
1874            retval = 0
1875            if abs(af[-1]-best_state.cost) > feps*10:
1876                retval = 5
1877                print "Warning: Cooled to %f at %s but this is not" \
1878                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
1879                      + " the smallest point found."
1880            break
1881        if (Tf is not None) and (schedule.T < Tf):
1882            retval = 1
1883            break
1884        if (maxeval is not None) and (schedule.feval > maxeval):
1885            retval = 2
1886            break
1887        if (iters > maxiter):
1888            print "Warning: Maximum number of iterations exceeded."
1889            retval = 3
1890            break
1891        if (maxaccept is not None) and (schedule.accepted > maxaccept):
1892            retval = 4
1893            break
1894
1895    if full_output:
1896        return best_state.x, best_state.cost, schedule.T, \
1897               schedule.feval, iters, schedule.accepted, retval
1898    else:
1899        return best_state.x, retval
1900
1901
1902def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
1903    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.)
1904   
1905    def getMDparms(item,pfx,parmDict,varyList):
1906        parmDict[pfx+'MDaxis'] = item['axis']
1907        parmDict[pfx+'MDval'] = item['Coef'][0]
1908        if item['Coef'][1]:
1909            varyList += [pfx+'MDval',]
1910            limits = item['Coef'][2]
1911            lower.append(limits[0])
1912            upper.append(limits[1])
1913           
1914    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
1915        parmDict[pfx+'Atype'] = item['atType']
1916        aTypes |= set([item['atType'],]) 
1917        pstr = ['Ax','Ay','Az']
1918        XYZ = [0,0,0]
1919        for i in range(3):
1920            name = pfx+pstr[i]
1921            parmDict[name] = item['Pos'][0][i]
1922            XYZ[i] = parmDict[name]
1923            if item['Pos'][1][i]:
1924                varyList += [name,]
1925                limits = item['Pos'][2][i]
1926                lower.append(limits[0])
1927                upper.append(limits[1])
1928        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
1929           
1930    def getRBparms(item,mfx,aTypes,RBdata,atNo,parmDict,varyList):
1931        parmDict[mfx+'MolCent'] = item['MolCent']
1932        parmDict[mfx+'RBId'] = item['RBId']
1933        pstr = ['Px','Py','Pz']
1934        ostr = ['Qa','Qi','Qj','Qk']
1935        for i in range(3):
1936            name = mfx+pstr[i]
1937            parmDict[name] = item['Pos'][0][i]
1938            if item['Pos'][1][i]:
1939                varyList += [name,]
1940                limits = item['Pos'][2][i]
1941                lower.append(limits[0])
1942                upper.append(limits[1])
1943        for i in range(4):
1944            name = mfx+ostr[i]
1945            parmDict[name] = item['Ori'][0][i]
1946            if item['Ovar'] == 'AV' and i:
1947                varyList += [name,]
1948                limits = item['Ori'][2][i]
1949                lower.append(limits[0])
1950                upper.append(limits[1])
1951            elif item['Ovar'] == 'A' and not i:
1952                varyList += [name,]
1953                limits = item['Ori'][2][i]
1954                lower.append(limits[0])
1955                upper.append(limits[1])
1956        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
1957            for i in range(len(item['Tor'][0])):
1958                name = mfx+'Tor'+str(i)
1959                parmDict[name] = item['Tor'][0][i]
1960                if item['Tor'][1][i]:
1961                    varyList += [name,]
1962                    limits = item['Tor'][2][i]
1963                    lower.append(limits[0])
1964                    upper.append(limits[1])
1965        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
1966        aTypes |= set(atypes)
1967        atNo += len(atypes)
1968        return atNo
1969
1970    def GetAtomTMX(pfx,RBdata,parmDict):
1971        'Needs a doc string'
1972        atNo = parmDIct['atNo']
1973        nfixAt = parmDict['nfixAt']
1974        Tdata = atNo*[' ',]
1975        Mdata = np.zeros(atNo)
1976        Fdata = np.zeros(atNo)
1977        Xdata = np.zeros((3,atNo))
1978        keys = {'Atype:':Tdata,'Amul:':Mdata,
1979            'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2]}
1980        nObjs = parmDict['nObjs']
1981        for iatm in range(nfixAt):
1982            for key in keys:
1983                parm = pfx+key+str(iatm)
1984                if parm in parmDict:
1985                    keys[key][iatm] = parmDict[parm]
1986        return Tdata,Mdata,Xdata
1987   
1988#    def mcsaSfCalc(refList,RBdata,G,SGData,parmDict):
1989#        ''' Compute structure factors for all h,k,l for phase
1990#        input:
1991#            refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv]
1992#            G:      reciprocal metric tensor
1993#            SGData: space group info. dictionary output from SpcGroup
1994#            ParmDict:
1995#        puts result F^2 in each ref[8] in refList
1996#        '''       
1997#        twopi = 2.0*np.pi
1998#        twopisq = 2.0*np.pi**2
1999#        ast = np.sqrt(np.diag(G))
2000#        Mast = twopisq*np.multiply.outer(ast,ast)
2001#        Tdata,Mdata,Fdata,Xdata = GetAtomFX(pfx,calcControls,parmDict)
2002#        FF = np.zeros(len(Tdata))
2003#        if 'N' in parmDict[hfx+'Type']:
2004#            FP,FPP = G2el.BlenRes(Tdata,BLtables,parmDict[hfx+'Lam'])
2005#        else:
2006#            FP = np.array([FFtables[El][hfx+'FP'] for El in Tdata])
2007#            FPP = np.array([FFtables[El][hfx+'FPP'] for El in Tdata])
2008#        Uij = np.array(G2lat.U6toUij(Uijdata))
2009#        bij = Mast*Uij.T
2010#        for refl in refList:
2011#            fbs = np.array([0,0])
2012#            H = refl[:3]
2013#            SQ = 1./(2.*refl[4])**2
2014#            SQfactor = 4.0*SQ*twopisq
2015#            Bab = parmDict[phfx+'BabA']*np.exp(-parmDict[phfx+'BabU']*SQfactor)
2016#            if not len(refl[-1]):                #no form factors
2017#                if 'N' in parmDict[hfx+'Type']:
2018#                    refl[-1] = getBLvalues(BLtables)
2019#                else:       #'X'
2020#                    refl[-1] = getFFvalues(FFtables,SQ)
2021#            for i,El in enumerate(Tdata):
2022#                FF[i] = refl[-1][El]           
2023#            Uniq = refl[11]
2024#            phi = refl[12]
2025#            phase = twopi*(np.inner(Uniq,(Xdata.T))+phi[:,np.newaxis])
2026#            sinp = np.sin(phase)
2027#            cosp = np.cos(phase)
2028#            occ = Mdata*Fdata/len(Uniq)
2029#            fa = np.array([(FF+FP-Bab)*occ*cosp,-FPP*occ*sinp])
2030#            fas = np.sum(np.sum(fa,axis=1),axis=1)        #real
2031#            if not SGData['SGInv']:
2032#                fb = np.array([(FF+FP-Bab)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr])
2033#                fbs = np.sum(np.sum(fb,axis=1),axis=1)
2034#            fasq = fas**2
2035#            fbsq = fbs**2        #imaginary
2036#            refl[9] = np.sum(fasq)+np.sum(fbsq)
2037#            refl[10] = atan2d(fbs[0],fas[0])
2038   
2039    sq8ln2 = np.sqrt(8*np.log(2))
2040    sq2pi = np.sqrt(2*np.pi)
2041    sq4pi = np.sqrt(4*np.pi)
2042    generalData = data['General']
2043    SGData = generalData['SGData']
2044    fixAtoms = data['Atoms']                       #if any
2045    cx,ct,cs = generalData['AtomPtrs'][:3]
2046    aTypes = set([])
2047    parmDict = {}
2048    varyList = []
2049    atNo = 0
2050    for atm in fixAtoms:
2051        pfx = ':'+str(atNo)+':'
2052        parmDict[pfx+'Atype'] = atm[ct]
2053        aTypes |= set([atm[ct],])
2054        pstr = ['Ax','Ay','Az']
2055        parmDict[pfx+'Amul'] = atm[cs+1]
2056        for i in range(3):
2057            name = pfx+pstr[i]
2058            parmDict[name] = atm[cx+i]
2059        atNo += 1
2060    parmDict['nfixAt'] = len(fixAtoms)       
2061    mcsaControls = generalData['MCSA controls']
2062    reflName = mcsaControls['Data source']
2063    phaseName = generalData['Name']
2064    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
2065    upper = []
2066    lower = []
2067    for i,item in enumerate(MCSAObjs):
2068        mfx = str(i)+':'
2069        parmDict[mfx+'Type'] = item['Type']
2070        if item['Type'] == 'MD':
2071            getMDparms(item,mfx,parmDict,varyList)
2072        elif item['Type'] == 'Atom':
2073            pfx = mfx+str(atNo)+':'
2074            getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList)
2075            atNo += 1
2076        elif item['Type'] in ['Residue','Vector']:
2077            pfx = mfx+':'
2078            atNo = getRBparms(item,pfx,aTypes,RBdata,atNo,parmDict,varyList)
2079    parmDict['atNo'] = atNo                 #total no. of atoms
2080    parmDict['nObj'] = len(MCSAObjs)
2081    FFtables = G2el.GetFFtable(aTypes)
2082    refs = []
2083    if 'PWDR' in reflName:
2084        for ref in reflData:
2085            h,k,l,m,d,pos,sig,gam,f = ref[:9]
2086            if d >= mcsaControls['dmin']:
2087                sig = gamFW(sig,gam)/sq8ln2        #--> sig from FWHM
2088                SQ = 0.25/d**2
2089                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2090                FFs = G2el.getFFvalues(FFtables,SQ)
2091                refs.append([h,k,l,m,f*m,pos,sig,FFs,Uniq,phi])
2092        nRef = len(refs)
2093        rcov = np.zeros((nRef,nRef))
2094        for iref,refI in enumerate(refs):
2095            rcov[iref][iref] = 1./(sq4pi*refI[6])
2096            for jref,refJ in enumerate(refs[:iref]):
2097                t1 = refI[6]**2+refJ[6]**2
2098                t2 = (refJ[5]-refI[5])**2/(2.*t1)
2099                if t2 > 10.:
2100                    rcov[iref][jref] = 0.
2101                else:
2102                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
2103        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2104        Rdiag = np.sqrt(np.diag(rcov))
2105        Rnorm = np.outer(Rdiag,Rdiag)
2106        rcov /= Rnorm
2107    elif 'Pawley' in reflName:
2108        covMatrix = covData['covMatrix']
2109        vList = covData['varyList']
2110        for iref,refI in enumerate(reflData):
2111            h,k,l,m,d,v,f,s = refI
2112            if d >= mcsaControls['dmin'] and v:       #skip unrefined ones
2113                SQ = 0.25/d**2
2114                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2115                FFs = G2el.getFFvalues(FFtables,SQ)
2116                refs.append([h,k,l,m,f*m,iref,0.,FFs,Uniq,phi])
2117        nRef = len(refs)
2118        pfx = str(data['pId'])+'::PWLref:'
2119        rcov = np.zeros((nRef,nRef))       
2120        for iref,refI in enumerate(refs):
2121            I = refI[5]
2122            nameI = pfx+str(I)
2123            if nameI in vList:
2124                Iindx = vList.index(nameI)
2125                rcov[iref][iref] = covMatrix[Iindx][Iindx]
2126                for jref,refJ in enumerate(refs[:iref]):
2127                    J = refJ[5]
2128                    nameJ = pfx+str(J)
2129                    try:
2130                        rcov[iref][jref] = covMatrix[vList.index(nameI)][vList.index(nameJ)]
2131                    except ValueError:
2132                        rcov[iref][jref] = rcov[iref][jref-1]
2133            else:
2134                rcov[iref] = rcov[iref-1]
2135                rcov[iref][iref] = rcov[iref-1][iref-1]
2136        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2137        Rdiag = np.sqrt(np.diag(rcov))
2138        Rnorm = np.outer(Rdiag,Rdiag)
2139        rcov /= Rnorm
2140    elif 'HKLF' in reflName:
2141        for ref in reflData:
2142            [h,k,l,m,d],f = ref[:5],ref[6]
2143            if d >= mcsaControls['dmin']:
2144                SQ = 0.25/d**2
2145                Uniq,phi = G2spc.GenHKLf([h,k,l],SGData)[2:]
2146                FFs = G2el.getFFvalues(FFtables,SQ)
2147                refs.append([h,k,l,m,f*m,0.,0.,FFs,Uniq,phi])
2148        rcov = np.identity(len(refs))
2149       
2150    for parm in parmDict:
2151        print parm,parmDict[parm] 
2152                       
2153#    XYZ,aTypes = UpdateMCSAxyz(Bmat,MCSA)       
2154           
2155#    generalData['MCSA controls'] = {'Data source':'','Annealing':[50.,0.001,50,1.e-6],
2156#    'dmin':2.0,'Algorithm':'fast','Jump coeff':[0.95,0.5],'nRuns':1,'boltzmann':1.0,
2157#    'fast parms':[1.0,1.0,1.0],'log slope':0.9}
2158
2159    return {}
2160
2161       
2162################################################################################
2163##### Quaternion stuff
2164################################################################################
2165
2166def prodQQ(QA,QB):
2167    ''' Grassman quaternion product
2168        QA,QB quaternions; q=r+ai+bj+ck
2169    '''
2170    D = np.zeros(4)
2171    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
2172    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
2173    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
2174    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
2175    return D
2176   
2177def normQ(QA):
2178    ''' get length of quaternion & normalize it
2179        q=r+ai+bj+ck
2180    '''
2181    n = np.sqrt(np.sum(np.array(QA)**2))
2182    return QA/n
2183   
2184def invQ(Q):
2185    '''
2186        get inverse of quaternion
2187        q=r+ai+bj+ck; q* = r-ai-bj-ck
2188    '''
2189    return Q*np.array([1,-1,-1,-1])
2190   
2191def prodQVQ(Q,V):
2192    """
2193    compute the quaternion vector rotation qvq-1 = v'
2194    q=r+ai+bj+ck
2195    """
2196    VP = np.zeros(3)
2197    T2 = Q[0]*Q[1]
2198    T3 = Q[0]*Q[2]
2199    T4 = Q[0]*Q[3]
2200    T5 = -Q[1]*Q[1]
2201    T6 = Q[1]*Q[2]
2202    T7 = Q[1]*Q[3]
2203    T8 = -Q[2]*Q[2]
2204    T9 = Q[2]*Q[3]
2205    T10 = -Q[3]*Q[3]
2206    VP[0] = 2.*((T8+T10)*V[0]+(T6-T4)*V[1]+(T3+T7)*V[2])+V[0]
2207    VP[1] = 2.*((T4+T6)*V[0]+(T5+T10)*V[1]+(T9-T2)*V[2])+V[1]
2208    VP[2] = 2.*((T7-T3)*V[0]+(T2+T9)*V[1]+(T5+T8)*V[2])+V[2] 
2209    return VP   
2210   
2211def Q2Mat(Q):
2212    ''' make rotation matrix from quaternion
2213        q=r+ai+bj+ck
2214    '''
2215    QN = normQ(Q)
2216    aa = QN[0]**2
2217    ab = QN[0]*QN[1]
2218    ac = QN[0]*QN[2]
2219    ad = QN[0]*QN[3]
2220    bb = QN[1]**2
2221    bc = QN[1]*QN[2]
2222    bd = QN[1]*QN[3]
2223    cc = QN[2]**2
2224    cd = QN[2]*QN[3]
2225    dd = QN[3]**2
2226    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
2227        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
2228        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
2229    return np.array(M)
2230   
2231def AV2Q(A,V):
2232    ''' convert angle (radians) & vector to quaternion
2233        q=r+ai+bj+ck
2234    '''
2235    Q = np.zeros(4)
2236    d = np.sqrt(np.sum(np.array(V)**2))
2237    if d:
2238        V /= d
2239    else:
2240        V = np.array([0.,0.,1.])
2241    p = A/2.
2242    Q[0] = np.cos(p)
2243    Q[1:4] = V*np.sin(p)
2244    return Q
2245   
2246def AVdeg2Q(A,V):
2247    ''' convert angle (degrees) & vector to quaternion
2248        q=r+ai+bj+ck
2249    '''
2250    Q = np.zeros(4)
2251    d = np.sqrt(np.sum(np.array(V)**2))
2252    if d:
2253        V /= d
2254    else:
2255        V = np.array([0.,0.,1.])
2256    p = A/2.
2257    Q[0] = cosd(p)
2258    Q[1:4] = V*sind(p)
2259    return Q
2260   
2261def Q2AVdeg(Q):
2262    ''' convert quaternion to angle (degrees 0-360) & normalized vector
2263        q=r+ai+bj+ck
2264    '''
2265    A = 2.*acosd(Q[0])
2266    V = np.array(Q[1:])
2267    d = np.sqrt(np.sum(V**2))
2268    if d:
2269        V /= d
2270    else:
2271        A = 0.
2272        V = np.array([0.,0.,0.])
2273    return A,V
2274   
2275def Q2AV(Q):
2276    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
2277        q=r+ai+bj+ck
2278    '''
2279    A = 2.*np.arccos(Q[0])
2280    V = np.array(Q[1:])
2281    d = np.sqrt(np.sum(V**2))
2282    if d:
2283        V /= d
2284    else:
2285        A = 0.
2286        V = np.array([0.,0.,0.])
2287    return A,V
2288   
2289def makeQuat(A,B,C):
2290    ''' Make quaternion from rotation of A vector to B vector about C axis
2291
2292        :param np.array A,B,C: Cartesian 3-vectors
2293        :Returns: quaternion & rotation angle in radians q=r+ai+bj+ck
2294    '''
2295
2296    V1 = np.cross(A,C)
2297    V2 = np.cross(B,C)
2298    if nl.norm(V1)*nl.norm(V2):
2299        V1 /= nl.norm(V1)
2300        V2 /= nl.norm(V2)
2301        V3 = np.cross(V1,V2)
2302    else:
2303        V3 = np.zeros(3)
2304    Q = np.array([0.,0.,0.,1.])
2305    D = 0.
2306    if nl.norm(V3):
2307        V3 /= nl.norm(V3)
2308        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
2309        D = np.arccos(D1)/2.0
2310        V1 = C-V3
2311        V2 = C+V3
2312        DM = nl.norm(V1)
2313        DP = nl.norm(V2)
2314        S = np.sin(D)
2315        Q[0] = np.cos(D)
2316        Q[1:] = V3*S
2317        D *= 2.
2318        if DM > DP:
2319            D *= -1.
2320    return Q,D
2321   
2322if __name__ == "__main__":
2323    from numpy import cos
2324    # minimum expected at ~-0.195
2325    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
2326    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
2327    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
2328    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
2329
2330    # minimum expected at ~[-0.195, -0.1]
2331    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
2332    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')
2333    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')
2334    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.