source: trunk/GSASIImath.py @ 1683

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

further fixes to FindMolecule? - remove SymOp? stuff; not needed & caused memory errors

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