source: trunk/GSASIImath.py @ 1674

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

bad response for extended structures & assemble molecule - now fixed

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