source: trunk/GSASIImath.py @ 954

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

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