source: trunk/GSASIImath.py @ 1363

Last change on this file since 1363 was 1363, checked in by vondreele, 8 years ago

fix to error in restraint processing

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 105.2 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2014-05-23 17:45:28 +0000 (Fri, 23 May 2014) $
5# $Author: vondreele $
6# $Revision: 1363 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1363 2014-05-23 17:45:28Z 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: 1363 $")
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    else:       #after last cycle or if zero cycles
155        M = func(x0,*args)
156        nfev += 1
157        Yvec,Amat = Hess(x0,*args)
158        Adiag = np.sqrt(np.diag(Amat))
159        Anorm = np.outer(Adiag,Adiag)
160        Amatlam = Amat*(One+Lam)/Anorm              #scale Amat to Marquardt array       
161    try:
162        Bmat = nl.inv(Amatlam)*(One+Lam)/Anorm      #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 LS'
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]+1)/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    '''default doc string
1999   
2000    :param type name: description
2001   
2002    :returns: type name: description
2003   
2004    '''
2005    tp = tand(pos/2.0)
2006    return ins['U']*tp**2+ins['V']*tp+ins['W']
2007   
2008def getCWsigDeriv(pos):
2009    '''default doc string
2010   
2011    :param type name: description
2012   
2013    :returns: type name: description
2014   
2015    '''
2016    tp = tand(pos/2.0)
2017    return tp**2,tp,1.0
2018   
2019def getCWgam(ins,pos):
2020    '''default doc string
2021   
2022    :param type name: description
2023   
2024    :returns: type name: description
2025   
2026    '''
2027    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
2028   
2029def getCWgamDeriv(pos):
2030    '''default doc string
2031   
2032    :param type name: description
2033   
2034    :returns: type name: description
2035   
2036    '''
2037    return 1./cosd(pos/2.0),tand(pos/2.0)
2038   
2039def getTOFsig(ins,dsp):
2040    '''default doc string
2041   
2042    :param type name: description
2043   
2044    :returns: type name: description
2045   
2046    '''
2047    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-q']*dsp
2048   
2049def getTOFsigDeriv(dsp):
2050    '''default doc string
2051   
2052    :param type name: description
2053   
2054    :returns: type name: description
2055   
2056    '''
2057    return 1.0,dsp**2,dsp
2058   
2059def getTOFgamma(ins,dsp):
2060    '''default doc string
2061   
2062    :param type name: description
2063   
2064    :returns: type name: description
2065   
2066    '''
2067    return ins['X']*dsp+ins['Y']*dsp**2
2068   
2069def getTOFgammaDeriv(dsp):
2070    '''default doc string
2071   
2072    :param type name: description
2073   
2074    :returns: type name: description
2075   
2076    '''
2077    return dsp,dsp**2
2078   
2079def getTOFbeta(ins,dsp):
2080    '''default doc string
2081   
2082    :param type name: description
2083   
2084    :returns: type name: description
2085   
2086    '''
2087    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp
2088   
2089def getTOFbetaDeriv(dsp):
2090    '''default doc string
2091   
2092    :param type name: description
2093   
2094    :returns: type name: description
2095   
2096    '''
2097    return 1.0,1./dsp**4,1./dsp
2098   
2099def getTOFalpha(ins,dsp):
2100    '''default doc string
2101   
2102    :param type name: description
2103   
2104    :returns: type name: description
2105   
2106    '''
2107    return ins['alpha']/dsp
2108   
2109def getTOFalphaDeriv(dsp):
2110    '''default doc string
2111   
2112    :param type name: description
2113   
2114    :returns: type name: description
2115   
2116    '''
2117    return 1./dsp
2118   
2119def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
2120    '''default doc string
2121   
2122    :param type name: description
2123   
2124    :returns: type name: description
2125   
2126    '''
2127    ind = 0
2128    if useFit:
2129        ind = 1
2130    ins = {}
2131    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
2132        for x in ['U','V','W','X','Y']:
2133            ins[x] = Parms[x][ind]
2134        if ifQ:                              #qplot - convert back to 2-theta
2135            pos = 2.0*asind(pos*wave/(4*math.pi))
2136        sig = getCWsig(ins,pos)
2137        gam = getCWgam(ins,pos)           
2138        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
2139    else:
2140        if ifQ:
2141            dsp = 2.*np.pi/pos
2142            pos = Parms['difC']*dsp
2143        else:
2144            dsp = pos/Parms['difC'][1]
2145        if 'Pdabc' in Parms2:
2146            for x in ['sig-0','sig-1','sig-q','X','Y']:
2147                ins[x] = Parms[x][ind]
2148            Pdabc = Parms2['Pdabc'].T
2149            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2150            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2151        else:
2152            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-q','X','Y']:
2153                ins[x] = Parms[x][ind]
2154            alp = getTOFalpha(ins,dsp)
2155            bet = getTOFbeta(ins,dsp)
2156        sig = getTOFsig(ins,dsp)
2157        gam = getTOFgamma(ins,dsp)
2158        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2159    return XY
2160   
2161################################################################################
2162##### MC/SA stuff
2163################################################################################
2164
2165#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2166# Original Author: Travis Oliphant 2002
2167# Bug-fixes in 2006 by Tim Leslie
2168
2169
2170import numpy
2171from numpy import asarray, tan, exp, ones, squeeze, sign, \
2172     all, log, sqrt, pi, shape, array, minimum, where
2173from numpy import random
2174
2175#__all__ = ['anneal']
2176
2177_double_min = numpy.finfo(float).min
2178_double_max = numpy.finfo(float).max
2179class base_schedule(object):
2180    def __init__(self):
2181        self.dwell = 20
2182        self.learn_rate = 0.5
2183        self.lower = -10
2184        self.upper = 10
2185        self.Ninit = 50
2186        self.accepted = 0
2187        self.tests = 0
2188        self.feval = 0
2189        self.k = 0
2190        self.T = None
2191
2192    def init(self, **options):
2193        self.__dict__.update(options)
2194        self.lower = asarray(self.lower)
2195        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2196        self.upper = asarray(self.upper)
2197        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2198        self.k = 0
2199        self.accepted = 0
2200        self.feval = 0
2201        self.tests = 0
2202
2203    def getstart_temp(self, best_state):
2204        """ Find a matching starting temperature and starting parameters vector
2205        i.e. find x0 such that func(x0) = T0.
2206
2207        Parameters
2208        ----------
2209        best_state : _state
2210            A _state object to store the function value and x0 found.
2211
2212        returns
2213        -------
2214        x0 : array
2215            The starting parameters vector.
2216        """
2217
2218        assert(not self.dims is None)
2219        lrange = self.lower
2220        urange = self.upper
2221        fmax = _double_min
2222        fmin = _double_max
2223        for _ in range(self.Ninit):
2224            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2225            fval = self.func(x0, *self.args)
2226            self.feval += 1
2227            if fval > fmax:
2228                fmax = fval
2229            if fval < fmin:
2230                fmin = fval
2231                best_state.cost = fval
2232                best_state.x = array(x0)
2233
2234        self.T0 = (fmax-fmin)*1.5
2235        return best_state.x
2236       
2237    def set_range(self,x0,frac):
2238        delrange = frac*(self.upper-self.lower)
2239        self.upper = x0+delrange
2240        self.lower = x0-delrange
2241
2242    def accept_test(self, dE):
2243        T = self.T
2244        self.tests += 1
2245        if dE < 0:
2246            self.accepted += 1
2247            return 1
2248        p = exp(-dE*1.0/self.boltzmann/T)
2249        if (p > random.uniform(0.0, 1.0)):
2250            self.accepted += 1
2251            return 1
2252        return 0
2253
2254    def update_guess(self, x0):
2255        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2256
2257    def update_temp(self, x0):
2258        pass
2259
2260
2261#  A schedule due to Lester Ingber modified to use bounds - OK
2262class fast_sa(base_schedule):
2263    def init(self, **options):
2264        self.__dict__.update(options)
2265        if self.m is None:
2266            self.m = 1.0
2267        if self.n is None:
2268            self.n = 1.0
2269        self.c = self.m * exp(-self.n * self.quench)
2270
2271    def update_guess(self, x0):
2272        x0 = asarray(x0)
2273        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2274        T = self.T
2275        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2276        xnew = xc*(self.upper - self.lower)+self.lower
2277        return xnew
2278#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2279#        xc = y*(self.upper - self.lower)
2280#        xnew = x0 + xc
2281#        return xnew
2282
2283    def update_temp(self):
2284        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2285        self.k += 1
2286        return
2287
2288class cauchy_sa(base_schedule):     #modified to use bounds - not good
2289    def update_guess(self, x0):
2290        x0 = asarray(x0)
2291        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2292        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2293        xnew = xc*(self.upper - self.lower)+self.lower
2294        return xnew
2295#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2296#        xc = self.learn_rate * self.T * tan(numbers)
2297#        xnew = x0 + xc
2298#        return xnew
2299
2300    def update_temp(self):
2301        self.T = self.T0/(1+self.k)
2302        self.k += 1
2303        return
2304
2305class boltzmann_sa(base_schedule):
2306#    def update_guess(self, x0):
2307#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2308#        x0 = asarray(x0)
2309#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2310#
2311#        xnew = x0 + xc*std*self.learn_rate
2312#        return xnew
2313
2314    def update_temp(self):
2315        self.k += 1
2316        self.T = self.T0 / log(self.k+1.0)
2317        return
2318
2319class log_sa(base_schedule):        #OK
2320
2321    def init(self,**options):
2322        self.__dict__.update(options)
2323       
2324    def update_guess(self,x0):     #same as default
2325        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2326       
2327    def update_temp(self):
2328        self.k += 1
2329        self.T = self.T0*self.slope**self.k
2330       
2331class _state(object):
2332    def __init__(self):
2333        self.x = None
2334        self.cost = None
2335
2336# TODO:
2337#     allow for general annealing temperature profile
2338#     in that case use update given by alpha and omega and
2339#     variation of all previous updates and temperature?
2340
2341# Simulated annealing
2342
2343def anneal(func, x0, args=(), schedule='fast', full_output=0,
2344           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2345           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2346           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
2347           ranRange=0.10,autoRan=False,dlg=None):
2348    """Minimize a function using simulated annealing.
2349
2350    Schedule is a schedule class implementing the annealing schedule.
2351    Available ones are 'fast', 'cauchy', 'boltzmann'
2352
2353    :param callable func: f(x, \*args)
2354        Function to be optimized.
2355    :param ndarray x0:
2356        Initial guess.
2357    :param tuple args:
2358        Extra parameters to `func`.
2359    :param base_schedule schedule:
2360        Annealing schedule to use (a class).
2361    :param bool full_output:
2362        Whether to return optional outputs.
2363    :param float T0:
2364        Initial Temperature (estimated as 1.2 times the largest
2365        cost-function deviation over random points in the range).
2366    :param float Tf:
2367        Final goal temperature.
2368    :param int maxeval:
2369        Maximum function evaluations.
2370    :param int maxaccept:
2371        Maximum changes to accept.
2372    :param int maxiter:
2373        Maximum cooling iterations.
2374    :param float learn_rate:
2375        Scale constant for adjusting guesses.
2376    :param float boltzmann:
2377        Boltzmann constant in acceptance test
2378        (increase for less stringent test at each temperature).
2379    :param float feps:
2380        Stopping relative error tolerance for the function value in
2381        last four coolings.
2382    :param float quench,m,n:
2383        Parameters to alter fast_sa schedule.
2384    :param float/ndarray lower,upper:
2385        Lower and upper bounds on `x`.
2386    :param int dwell:
2387        The number of times to search the space at each temperature.
2388    :param float slope:
2389        Parameter for log schedule
2390    :param bool ranStart=False:
2391        True for set 10% of ranges about x
2392
2393    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2394
2395     * xmin (ndarray): Point giving smallest value found.
2396     * Jmin (float): Minimum value of function found.
2397     * T (float): Final temperature.
2398     * feval (int): Number of function evaluations.
2399     * iters (int): Number of cooling iterations.
2400     * accept (int): Number of tests accepted.
2401     * retval (int): Flag indicating stopping condition:
2402
2403              *  0: Points no longer changing
2404              *  1: Cooled to final temperature
2405              *  2: Maximum function evaluations
2406              *  3: Maximum cooling iterations reached
2407              *  4: Maximum accepted query locations reached
2408              *  5: Final point not the minimum amongst encountered points
2409
2410    *Notes*:
2411    Simulated annealing is a random algorithm which uses no derivative
2412    information from the function being optimized. In practice it has
2413    been more useful in discrete optimization than continuous
2414    optimization, as there are usually better algorithms for continuous
2415    optimization problems.
2416
2417    Some experimentation by trying the difference temperature
2418    schedules and altering their parameters is likely required to
2419    obtain good performance.
2420
2421    The randomness in the algorithm comes from random sampling in numpy.
2422    To obtain the same results you can call numpy.random.seed with the
2423    same seed immediately before calling scipy.optimize.anneal.
2424
2425    We give a brief description of how the three temperature schedules
2426    generate new points and vary their temperature. Temperatures are
2427    only updated with iterations in the outer loop. The inner loop is
2428    over xrange(dwell), and new points are generated for
2429    every iteration in the inner loop. (Though whether the proposed
2430    new points are accepted is probabilistic.)
2431
2432    For readability, let d denote the dimension of the inputs to func.
2433    Also, let x_old denote the previous state, and k denote the
2434    iteration number of the outer loop. All other variables not
2435    defined below are input variables to scipy.optimize.anneal itself.
2436
2437    In the 'fast' schedule the updates are ::
2438
2439        u ~ Uniform(0, 1, size=d)
2440        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2441        xc = y * (upper - lower)
2442        x_new = x_old + xc
2443
2444        c = n * exp(-n * quench)
2445        T_new = T0 * exp(-c * k**quench)
2446
2447
2448    In the 'cauchy' schedule the updates are ::
2449
2450        u ~ Uniform(-pi/2, pi/2, size=d)
2451        xc = learn_rate * T * tan(u)
2452        x_new = x_old + xc
2453
2454        T_new = T0 / (1+k)
2455
2456    In the 'boltzmann' schedule the updates are ::
2457
2458        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2459        y ~ Normal(0, std, size=d)
2460        x_new = x_old + learn_rate * y
2461
2462        T_new = T0 / log(1+k)
2463
2464    """
2465    x0 = asarray(x0)
2466    lower = asarray(lower)
2467    upper = asarray(upper)
2468
2469    schedule = eval(schedule+'_sa()')
2470    #   initialize the schedule
2471    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2472                  learn_rate=learn_rate, lower=lower, upper=upper,
2473                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2474
2475    current_state, last_state, best_state = _state(), _state(), _state()
2476    if ranStart:
2477        schedule.set_range(x0,ranRange)
2478    if T0 is None:
2479        x0 = schedule.getstart_temp(best_state)
2480    else:
2481        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
2482        best_state.x = None
2483        best_state.cost = numpy.Inf
2484
2485    last_state.x = asarray(x0).copy()
2486    fval = func(x0,*args)
2487    schedule.feval += 1
2488    last_state.cost = fval
2489    if last_state.cost < best_state.cost:
2490        best_state.cost = fval
2491        best_state.x = asarray(x0).copy()
2492    schedule.T = schedule.T0
2493    fqueue = [100, 300, 500, 700]
2494    iters = 1
2495    keepGoing = True
2496    bestn = 0
2497    while keepGoing:
2498        retval = 0
2499        for n in xrange(dwell):
2500            current_state.x = schedule.update_guess(last_state.x)
2501            current_state.cost = func(current_state.x,*args)
2502            schedule.feval += 1
2503
2504            dE = current_state.cost - last_state.cost
2505            if schedule.accept_test(dE):
2506                last_state.x = current_state.x.copy()
2507                last_state.cost = current_state.cost
2508                if last_state.cost < best_state.cost:
2509                    best_state.x = last_state.x.copy()
2510                    best_state.cost = last_state.cost
2511                    bestn = n
2512                    if best_state.cost < 1.0 and autoRan:
2513                        schedule.set_range(x0,best_state.cost/2.)                       
2514        if dlg:
2515            GoOn = dlg.Update(min(100.,best_state.cost*100),
2516                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
2517                    'Best trial:',bestn,  \
2518                    'MC/SA Residual:',best_state.cost*100,'%', \
2519                    ))[0]
2520            if not GoOn:
2521                best_state.x = last_state.x.copy()
2522                best_state.cost = last_state.cost
2523                retval = 5
2524        schedule.update_temp()
2525        iters += 1
2526        # Stopping conditions
2527        # 0) last saved values of f from each cooling step
2528        #     are all very similar (effectively cooled)
2529        # 1) Tf is set and we are below it
2530        # 2) maxeval is set and we are past it
2531        # 3) maxiter is set and we are past it
2532        # 4) maxaccept is set and we are past it
2533        # 5) user canceled run via progress bar
2534
2535        fqueue.append(squeeze(last_state.cost))
2536        fqueue.pop(0)
2537        af = asarray(fqueue)*1.0
2538        if retval == 5:
2539            print ' User terminated run; incomplete MC/SA'
2540            keepGoing = False
2541            break
2542        if all(abs((af-af[0])/af[0]) < feps):
2543            retval = 0
2544            if abs(af[-1]-best_state.cost) > feps*10:
2545                retval = 5
2546#                print "Warning: Cooled to %f at %s but this is not" \
2547#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2548#                      + " the smallest point found."
2549            break
2550        if (Tf is not None) and (schedule.T < Tf):
2551            retval = 1
2552            break
2553        if (maxeval is not None) and (schedule.feval > maxeval):
2554            retval = 2
2555            break
2556        if (iters > maxiter):
2557            print "Warning: Maximum number of iterations exceeded."
2558            retval = 3
2559            break
2560        if (maxaccept is not None) and (schedule.accepted > maxaccept):
2561            retval = 4
2562            break
2563
2564    if full_output:
2565        return best_state.x, best_state.cost, schedule.T, \
2566               schedule.feval, iters, schedule.accepted, retval
2567    else:
2568        return best_state.x, retval
2569
2570def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
2571    outlist = []
2572    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
2573    for n in range(iCyc):
2574        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
2575        outlist.append(result[0])
2576        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
2577    out_q.put(outlist)
2578
2579def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
2580    import multiprocessing as mp
2581   
2582    nprocs = mp.cpu_count()
2583    out_q = mp.Queue()
2584    procs = []
2585    iCyc = np.zeros(nprocs)
2586    for i in range(nCyc):
2587        iCyc[i%nprocs] += 1
2588    for i in range(nprocs):
2589        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
2590        procs.append(p)
2591        p.start()
2592    resultlist = []
2593    for i in range(nprocs):
2594        resultlist += out_q.get()
2595    for p in procs:
2596        p.join()
2597    return resultlist
2598
2599def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
2600    '''default doc string
2601   
2602    :param type name: description
2603   
2604    :returns: type name: description
2605   
2606    '''
2607   
2608    twopi = 2.0*np.pi
2609    global tsum
2610    tsum = 0.
2611   
2612    def getMDparms(item,pfx,parmDict,varyList):
2613        parmDict[pfx+'MDaxis'] = item['axis']
2614        parmDict[pfx+'MDval'] = item['Coef'][0]
2615        if item['Coef'][1]:
2616            varyList += [pfx+'MDval',]
2617            limits = item['Coef'][2]
2618            lower.append(limits[0])
2619            upper.append(limits[1])
2620                       
2621    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
2622        parmDict[pfx+'Atype'] = item['atType']
2623        aTypes |= set([item['atType'],]) 
2624        pstr = ['Ax','Ay','Az']
2625        XYZ = [0,0,0]
2626        for i in range(3):
2627            name = pfx+pstr[i]
2628            parmDict[name] = item['Pos'][0][i]
2629            XYZ[i] = parmDict[name]
2630            if item['Pos'][1][i]:
2631                varyList += [name,]
2632                limits = item['Pos'][2][i]
2633                lower.append(limits[0])
2634                upper.append(limits[1])
2635        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
2636           
2637    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
2638        parmDict[mfx+'MolCent'] = item['MolCent']
2639        parmDict[mfx+'RBId'] = item['RBId']
2640        pstr = ['Px','Py','Pz']
2641        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
2642        for i in range(3):
2643            name = mfx+pstr[i]
2644            parmDict[name] = item['Pos'][0][i]
2645            if item['Pos'][1][i]:
2646                varyList += [name,]
2647                limits = item['Pos'][2][i]
2648                lower.append(limits[0])
2649                upper.append(limits[1])
2650        AV = item['Ori'][0]
2651        A = AV[0]
2652        V = AV[1:]
2653        for i in range(4):
2654            name = mfx+ostr[i]
2655            if i:
2656                parmDict[name] = V[i-1]
2657            else:
2658                parmDict[name] = A
2659            if item['Ovar'] == 'AV':
2660                varyList += [name,]
2661                limits = item['Ori'][2][i]
2662                lower.append(limits[0])
2663                upper.append(limits[1])
2664            elif item['Ovar'] == 'A' and not i:
2665                varyList += [name,]
2666                limits = item['Ori'][2][i]
2667                lower.append(limits[0])
2668                upper.append(limits[1])
2669        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
2670            for i in range(len(item['Tor'][0])):
2671                name = mfx+'Tor'+str(i)
2672                parmDict[name] = item['Tor'][0][i]
2673                if item['Tor'][1][i]:
2674                    varyList += [name,]
2675                    limits = item['Tor'][2][i]
2676                    lower.append(limits[0])
2677                    upper.append(limits[1])
2678        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
2679        aTypes |= set(atypes)
2680        atNo += len(atypes)
2681        return atNo
2682               
2683    def GetAtomM(Xdata,SGData):
2684        Mdata = []
2685        for xyz in Xdata:
2686            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
2687        return np.array(Mdata)
2688       
2689    def GetAtomT(RBdata,parmDict):
2690        'Needs a doc string'
2691        atNo = parmDict['atNo']
2692        nfixAt = parmDict['nfixAt']
2693        Tdata = atNo*[' ',]
2694        for iatm in range(nfixAt):
2695            parm = ':'+str(iatm)+':Atype'
2696            if parm in parmDict:
2697                Tdata[iatm] = aTypes.index(parmDict[parm])
2698        iatm = nfixAt
2699        for iObj in range(parmDict['nObj']):
2700            pfx = str(iObj)+':'
2701            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2702                if parmDict[pfx+'Type'] == 'Vector':
2703                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2704                    nAtm = len(RBRes['rbVect'][0])
2705                else:       #Residue
2706                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2707                    nAtm = len(RBRes['rbXYZ'])
2708                for i in range(nAtm):
2709                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
2710                    iatm += 1
2711            elif parmDict[pfx+'Type'] == 'Atom':
2712                atNo = parmDict[pfx+'atNo']
2713                parm = pfx+'Atype'              #remove extra ':'
2714                if parm in parmDict:
2715                    Tdata[atNo] = aTypes.index(parmDict[parm])
2716                iatm += 1
2717            else:
2718                continue        #skips March Dollase
2719        return Tdata
2720       
2721    def GetAtomX(RBdata,parmDict):
2722        'Needs a doc string'
2723        Bmat = parmDict['Bmat']
2724        atNo = parmDict['atNo']
2725        nfixAt = parmDict['nfixAt']
2726        Xdata = np.zeros((3,atNo))
2727        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
2728        for iatm in range(nfixAt):
2729            for key in keys:
2730                parm = ':'+str(iatm)+key
2731                if parm in parmDict:
2732                    keys[key][iatm] = parmDict[parm]
2733        iatm = nfixAt
2734        for iObj in range(parmDict['nObj']):
2735            pfx = str(iObj)+':'
2736            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2737                if parmDict[pfx+'Type'] == 'Vector':
2738                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
2739                    vecs = RBRes['rbVect']
2740                    mags = RBRes['VectMag']
2741                    Cart = np.zeros_like(vecs[0])
2742                    for vec,mag in zip(vecs,mags):
2743                        Cart += vec*mag
2744                elif parmDict[pfx+'Type'] == 'Residue':
2745                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
2746                    Cart = np.array(RBRes['rbXYZ'])
2747                    for itor,seq in enumerate(RBRes['rbSeq']):
2748                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
2749                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
2750                if parmDict[pfx+'MolCent'][1]:
2751                    Cart -= parmDict[pfx+'MolCent'][0]
2752                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2753                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
2754                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
2755                iatm += len(Cart)
2756            elif parmDict[pfx+'Type'] == 'Atom':
2757                atNo = parmDict[pfx+'atNo']
2758                for key in keys:
2759                    parm = pfx+key[1:]              #remove extra ':'
2760                    if parm in parmDict:
2761                        keys[key][atNo] = parmDict[parm]
2762                iatm += 1
2763            else:
2764                continue        #skips March Dollase
2765        return Xdata.T
2766       
2767    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
2768        allX = np.inner(Xdata,SGM)+SGT
2769        allT = np.repeat(Tdata,allX.shape[1])
2770        allM = np.repeat(Mdata,allX.shape[1])
2771        allX = np.reshape(allX,(-1,3))
2772        return allT,allM,allX
2773
2774    def getAllX(Xdata,SGM,SGT):
2775        allX = np.inner(Xdata,SGM)+SGT
2776        allX = np.reshape(allX,(-1,3))
2777        return allX
2778       
2779    def normQuaternions(RBdata,parmDict,varyList,values):
2780        for iObj in range(parmDict['nObj']):
2781            pfx = str(iObj)+':'
2782            if parmDict[pfx+'Type'] in ['Vector','Residue']:
2783                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
2784                A,V = Q2AVdeg(Qori)
2785                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
2786                    if i:
2787                        parmDict[pfx+name] = V[i-1]
2788                    else:
2789                        parmDict[pfx+name] = A
2790       
2791    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
2792        ''' Compute structure factors for all h,k,l for phase
2793        input:
2794            refList: [ref] where each ref = h,k,l,m,d,...
2795            rcov:   array[nref,nref] covariance terms between Fo^2 values
2796            ifInv:  bool True if centrosymmetric
2797            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
2798            RBdata: [dict] rigid body dictionary
2799            varyList: [list] names of varied parameters in MC/SA (not used here)           
2800            ParmDict: [dict] problem parameters
2801        puts result F^2 in each ref[5] in refList
2802        returns:
2803            delt-F*rcov*delt-F/sum(Fo^2)^2
2804        '''       
2805        global tsum
2806        t0 = time.time()
2807        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
2808        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
2809        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
2810        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
2811        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
2812        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
2813        refList[5] = Aterm
2814        if not ifInv:
2815            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
2816        if len(cosTable):        #apply MD correction
2817            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
2818        sumFcsq = np.sum(refList[5])
2819        scale = parmDict['sumFosq']/sumFcsq
2820        refList[5] *= scale
2821        refList[6] = refList[4]-refList[5]
2822        M = np.inner(refList[6],np.inner(rcov,refList[6]))
2823        tsum += (time.time()-t0)
2824        return M/np.sum(refList[4]**2)
2825
2826    sq8ln2 = np.sqrt(8*np.log(2))
2827    sq2pi = np.sqrt(2*np.pi)
2828    sq4pi = np.sqrt(4*np.pi)
2829    generalData = data['General']
2830    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2831    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
2832    SGData = generalData['SGData']
2833    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
2834    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
2835    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
2836    fixAtoms = data['Atoms']                       #if any
2837    cx,ct,cs = generalData['AtomPtrs'][:3]
2838    aTypes = set([])
2839    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
2840    varyList = []
2841    atNo = 0
2842    for atm in fixAtoms:
2843        pfx = ':'+str(atNo)+':'
2844        parmDict[pfx+'Atype'] = atm[ct]
2845        aTypes |= set([atm[ct],])
2846        pstr = ['Ax','Ay','Az']
2847        parmDict[pfx+'Amul'] = atm[cs+1]
2848        for i in range(3):
2849            name = pfx+pstr[i]
2850            parmDict[name] = atm[cx+i]
2851        atNo += 1
2852    parmDict['nfixAt'] = len(fixAtoms)       
2853    MCSA = generalData['MCSA controls']
2854    reflName = MCSA['Data source']
2855    phaseName = generalData['Name']
2856    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
2857    upper = []
2858    lower = []
2859    MDvec = np.zeros(3)
2860    for i,item in enumerate(MCSAObjs):
2861        mfx = str(i)+':'
2862        parmDict[mfx+'Type'] = item['Type']
2863        if item['Type'] == 'MD':
2864            getMDparms(item,mfx,parmDict,varyList)
2865            MDvec = np.array(item['axis'])
2866        elif item['Type'] == 'Atom':
2867            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
2868            parmDict[mfx+'atNo'] = atNo
2869            atNo += 1
2870        elif item['Type'] in ['Residue','Vector']:
2871            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
2872    parmDict['atNo'] = atNo                 #total no. of atoms
2873    parmDict['nObj'] = len(MCSAObjs)
2874    aTypes = list(aTypes)
2875    Tdata = GetAtomT(RBdata,parmDict)
2876    Xdata = GetAtomX(RBdata,parmDict)
2877    Mdata = GetAtomM(Xdata,SGData)
2878    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
2879    FFtables = G2el.GetFFtable(aTypes)
2880    refs = []
2881    allFF = []
2882    cosTable = []
2883    sumFosq = 0
2884    if 'PWDR' in reflName:
2885        for ref in reflData:
2886            h,k,l,m,d,pos,sig,gam,f = ref[:9]
2887            if d >= MCSA['dmin']:
2888                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
2889                SQ = 0.25/d**2
2890                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2891                refs.append([h,k,l,m,f*m,pos,sig])
2892                sumFosq += f*m
2893                Heqv = np.inner(np.array([h,k,l]),SGMT)
2894                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
2895        nRef = len(refs)
2896        cosTable = np.array(cosTable)**2
2897        rcov = np.zeros((nRef,nRef))
2898        for iref,refI in enumerate(refs):
2899            rcov[iref][iref] = 1./(sq4pi*refI[6])
2900            for jref,refJ in enumerate(refs[:iref]):
2901                t1 = refI[6]**2+refJ[6]**2
2902                t2 = (refJ[5]-refI[5])**2/(2.*t1)
2903                if t2 > 10.:
2904                    rcov[iref][jref] = 0.
2905                else:
2906                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
2907        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2908        Rdiag = np.sqrt(np.diag(rcov))
2909        Rnorm = np.outer(Rdiag,Rdiag)
2910        rcov /= Rnorm
2911    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
2912        vNames = []
2913        pfx = str(data['pId'])+'::PWLref:'
2914        for iref,refI in enumerate(reflData):           #Pawley reflection set
2915            h,k,l,m,d,v,f,s = refI
2916            if d >= MCSA['dmin'] and v:       #skip unrefined ones
2917                vNames.append(pfx+str(iref))
2918                SQ = 0.25/d**2
2919                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2920                refs.append([h,k,l,m,f*m,iref,0.])
2921                sumFosq += f*m
2922                Heqv = np.inner(np.array([h,k,l]),SGMT)
2923                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
2924        cosTable = np.array(cosTable)**2
2925        nRef = len(refs)
2926        if covData['freshCOV'] and generalData['doPawley'] and MCSA.get('newDmin',True):
2927            vList = covData['varyList']
2928            covMatrix = covData['covMatrix']
2929            rcov = getVCov(vNames,vList,covMatrix)
2930            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
2931            Rdiag = np.sqrt(np.diag(rcov))
2932            Rdiag = np.where(Rdiag,Rdiag,1.0)
2933            Rnorm = np.outer(Rdiag,Rdiag)
2934            rcov /= Rnorm
2935            MCSA['rcov'] = rcov
2936            covData['freshCOV'] = False
2937            MCSA['newDmin'] = False
2938        else:
2939            rcov = MCSA['rcov']
2940    elif 'HKLF' in reflName:
2941        for ref in reflData:
2942            [h,k,l,m,d],f = ref[:5],ref[6]
2943            if d >= MCSA['dmin']:
2944                SQ = 0.25/d**2
2945                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
2946                refs.append([h,k,l,m,f*m,0.,0.])
2947                sumFosq += f*m
2948        nRef = len(refs)
2949        rcov = np.identity(len(refs))
2950    allFF = np.array(allFF).T
2951    refs = np.array(refs).T
2952    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
2953    print ' Number of parameters varied: %d'%(len(varyList))
2954    parmDict['sumFosq'] = sumFosq
2955    x0 = [parmDict[val] for val in varyList]
2956    ifInv = SGData['SGInv']
2957    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
2958        schedule=MCSA['Algorithm'], full_output=True,
2959        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
2960        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
2961        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
2962        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
2963        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
2964    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
2965#    for ref in refs.T:
2966#        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])
2967#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
2968    Result = [False,False,results[1],results[2],]+list(results[0])
2969    Result.append(varyList)
2970    return Result,tsum
2971
2972       
2973################################################################################
2974##### Quaternion stuff
2975################################################################################
2976
2977def prodQQ(QA,QB):
2978    ''' Grassman quaternion product
2979        QA,QB quaternions; q=r+ai+bj+ck
2980    '''
2981    D = np.zeros(4)
2982    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
2983    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
2984    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
2985    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
2986   
2987#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
2988#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
2989   
2990    return D
2991   
2992def normQ(QA):
2993    ''' get length of quaternion & normalize it
2994        q=r+ai+bj+ck
2995    '''
2996    n = np.sqrt(np.sum(np.array(QA)**2))
2997    return QA/n
2998   
2999def invQ(Q):
3000    '''
3001        get inverse of quaternion
3002        q=r+ai+bj+ck; q* = r-ai-bj-ck
3003    '''
3004    return Q*np.array([1,-1,-1,-1])
3005   
3006def prodQVQ(Q,V):
3007    """
3008    compute the quaternion vector rotation qvq-1 = v'
3009    q=r+ai+bj+ck
3010    """
3011    T2 = Q[0]*Q[1]
3012    T3 = Q[0]*Q[2]
3013    T4 = Q[0]*Q[3]
3014    T5 = -Q[1]*Q[1]
3015    T6 = Q[1]*Q[2]
3016    T7 = Q[1]*Q[3]
3017    T8 = -Q[2]*Q[2]
3018    T9 = Q[2]*Q[3]
3019    T10 = -Q[3]*Q[3]
3020    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
3021    VP = 2.*np.inner(V,M)
3022    return VP+V
3023   
3024def Q2Mat(Q):
3025    ''' make rotation matrix from quaternion
3026        q=r+ai+bj+ck
3027    '''
3028    QN = normQ(Q)
3029    aa = QN[0]**2
3030    ab = QN[0]*QN[1]
3031    ac = QN[0]*QN[2]
3032    ad = QN[0]*QN[3]
3033    bb = QN[1]**2
3034    bc = QN[1]*QN[2]
3035    bd = QN[1]*QN[3]
3036    cc = QN[2]**2
3037    cd = QN[2]*QN[3]
3038    dd = QN[3]**2
3039    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
3040        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
3041        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
3042    return np.array(M)
3043   
3044def AV2Q(A,V):
3045    ''' convert angle (radians) & vector to quaternion
3046        q=r+ai+bj+ck
3047    '''
3048    Q = np.zeros(4)
3049    d = nl.norm(np.array(V))
3050    if d:
3051        V /= d
3052        if not A:       #==0.
3053            A = 2.*np.pi
3054        p = A/2.
3055        Q[0] = np.cos(p)
3056        Q[1:4] = V*np.sin(p)
3057    else:
3058        Q[3] = 1.
3059    return Q
3060   
3061def AVdeg2Q(A,V):
3062    ''' convert angle (degrees) & vector to quaternion
3063        q=r+ai+bj+ck
3064    '''
3065    Q = np.zeros(4)
3066    d = nl.norm(np.array(V))
3067    if not A:       #== 0.!
3068        A = 360.
3069    if d:
3070        V /= d
3071        p = A/2.
3072        Q[0] = cosd(p)
3073        Q[1:4] = V*sind(p)
3074    else:
3075        Q[3] = 1.
3076    return Q
3077   
3078def Q2AVdeg(Q):
3079    ''' convert quaternion to angle (degrees 0-360) & normalized vector
3080        q=r+ai+bj+ck
3081    '''
3082    A = 2.*acosd(Q[0])
3083    V = np.array(Q[1:])
3084    V /= sind(A/2.)
3085    return A,V
3086   
3087def Q2AV(Q):
3088    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
3089        q=r+ai+bj+ck
3090    '''
3091    A = 2.*np.arccos(Q[0])
3092    V = np.array(Q[1:])
3093    V /= np.sin(A/2.)
3094    return A,V
3095   
3096def randomQ(r0,r1,r2,r3):
3097    ''' create random quaternion from 4 random numbers in range (-1,1)
3098    '''
3099    sum = 0
3100    Q = np.array(4)
3101    Q[0] = r0
3102    sum += Q[0]**2
3103    Q[1] = np.sqrt(1.-sum)*r1
3104    sum += Q[1]**2
3105    Q[2] = np.sqrt(1.-sum)*r2
3106    sum += Q[2]**2
3107    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
3108    return Q
3109   
3110def randomAVdeg(r0,r1,r2,r3):
3111    ''' create random angle (deg),vector from 4 random number in range (-1,1)
3112    '''
3113    return Q2AVdeg(randomQ(r0,r1,r2,r3))
3114   
3115def makeQuat(A,B,C):
3116    ''' Make quaternion from rotation of A vector to B vector about C axis
3117
3118        :param np.array A,B,C: Cartesian 3-vectors
3119        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
3120    '''
3121
3122    V1 = np.cross(A,C)
3123    V2 = np.cross(B,C)
3124    if nl.norm(V1)*nl.norm(V2):
3125        V1 /= nl.norm(V1)
3126        V2 /= nl.norm(V2)
3127        V3 = np.cross(V1,V2)
3128    else:
3129        V3 = np.zeros(3)
3130    Q = np.array([0.,0.,0.,1.])
3131    D = 0.
3132    if nl.norm(V3):
3133        V3 /= nl.norm(V3)
3134        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
3135        D = np.arccos(D1)/2.0
3136        V1 = C-V3
3137        V2 = C+V3
3138        DM = nl.norm(V1)
3139        DP = nl.norm(V2)
3140        S = np.sin(D)
3141        Q[0] = np.cos(D)
3142        Q[1:] = V3*S
3143        D *= 2.
3144        if DM > DP:
3145            D *= -1.
3146    return Q,D
3147   
3148def annealtests():
3149    from numpy import cos
3150    # minimum expected at ~-0.195
3151    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
3152    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
3153    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
3154    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
3155
3156    # minimum expected at ~[-0.195, -0.1]
3157    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
3158    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')
3159    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')
3160    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')
3161
3162
3163if __name__ == '__main__':
3164    annealtests()
Note: See TracBrowser for help on using the repository browser.