source: trunk/GSASIImath.py @ 956

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

snapshot of CIF dev

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