source: trunk/GSASIImath.py @ 1474

Last change on this file since 1474 was 1474, checked in by vondreele, 7 years ago

1st MC/SA tutorial
various MC/SA fixes
fix to background peak fitting for CW & TOF

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