source: trunk/GSASIImath.py @ 1780

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

force collapse of all children when main tree item is collapsed
begin texture fitting from seq. refinements
fix bugs in plotting if all SH coeff = 0
set background on Fcsq when delt = 3sig & 10 sig

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 123.4 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIImath - major mathematics routines
3########### SVN repository information ###################
4# $Date: 2015-04-07 22:14:33 +0000 (Tue, 07 Apr 2015) $
5# $Author: vondreele $
6# $Revision: 1780 $
7# $URL: trunk/GSASIImath.py $
8# $Id: GSASIImath.py 1780 2015-04-07 22:14:33Z 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: 1780 $")
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##### Texture fitting stuff
1592################################################################################
1593
1594def FitTexture(General,SamAngs,refData):
1595    print 'fit texture for '+General["Name"]
1596    SGData = General['SGData']
1597    Texture = General['SH Texture']
1598    if not Texture['Order']:
1599        return 'No spherical harmonics coefficients'
1600       
1601       
1602    return None
1603   
1604################################################################################
1605##### Fourier & charge flip stuff
1606################################################################################
1607
1608def adjHKLmax(SGData,Hmax):
1609    '''default doc string
1610   
1611    :param type name: description
1612   
1613    :returns: type name: description
1614   
1615    '''
1616    if SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm']:
1617        Hmax[0] = int(math.ceil(Hmax[0]/6.))*6
1618        Hmax[1] = int(math.ceil(Hmax[1]/6.))*6
1619        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
1620    else:
1621        Hmax[0] = int(math.ceil(Hmax[0]/4.))*4
1622        Hmax[1] = int(math.ceil(Hmax[1]/4.))*4
1623        Hmax[2] = int(math.ceil(Hmax[2]/4.))*4
1624
1625def OmitMap(data,reflDict,pgbar=None):
1626    '''default doc string
1627   
1628    :param type name: description
1629   
1630    :returns: type name: description
1631   
1632    '''
1633    generalData = data['General']
1634    if not generalData['Map']['MapType']:
1635        print '**** ERROR - Fourier map not defined'
1636        return
1637    mapData = generalData['Map']
1638    dmin = mapData['Resolution']
1639    SGData = generalData['SGData']
1640    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1641    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1642    cell = generalData['Cell'][1:8]       
1643    A = G2lat.cell2A(cell[:6])
1644    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1645    adjHKLmax(SGData,Hmax)
1646    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1647    time0 = time.time()
1648    for iref,ref in enumerate(reflDict['RefList']):
1649        if ref[4] >= dmin:
1650            Fosq,Fcsq,ph = ref[8:11]
1651            Uniq = np.inner(ref[:3],SGMT)
1652            Phi = np.inner(ref[:3],SGT)
1653            for i,hkl in enumerate(Uniq):        #uses uniq
1654                hkl = np.asarray(hkl,dtype='i')
1655                dp = 360.*Phi[i]                #and phi
1656                a = cosd(ph+dp)
1657                b = sind(ph+dp)
1658                phasep = complex(a,b)
1659                phasem = complex(a,-b)
1660                Fo = np.sqrt(Fosq)
1661                if '2Fo-Fc' in mapData['MapType']:
1662                    F = 2.*np.sqrt(Fosq)-np.sqrt(Fcsq)
1663                else:
1664                    F = np.sqrt(Fosq)
1665                h,k,l = hkl+Hmax
1666                Fhkl[h,k,l] = F*phasep
1667                h,k,l = -hkl+Hmax
1668                Fhkl[h,k,l] = F*phasem
1669    rho0 = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1670    M = np.mgrid[0:4,0:4,0:4]
1671    blkIds = np.array(zip(M[0].flatten(),M[1].flatten(),M[2].flatten()))
1672    iBeg = blkIds*rho0.shape/4
1673    iFin = (blkIds+1)*rho0.shape/4
1674    rho_omit = np.zeros_like(rho0)
1675    nBlk = 0
1676    for iB,iF in zip(iBeg,iFin):
1677        rho1 = np.copy(rho0)
1678        rho1[iB[0]:iF[0],iB[1]:iF[1],iB[2]:iF[2]] = 0.
1679        Fnew = fft.ifftshift(fft.ifftn(rho1))
1680        Fnew = np.where(Fnew,Fnew,1.0)           #avoid divide by zero
1681        phase = Fnew/np.absolute(Fnew)
1682        OFhkl = np.absolute(Fhkl)*phase
1683        rho1 = np.real(fft.fftn(fft.fftshift(OFhkl)))*(1.+0j)
1684        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]])
1685        nBlk += 1
1686        pgbar.Update(nBlk)
1687    mapData['rho'] = np.real(rho_omit)/cell[6]
1688    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1689    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
1690    print 'Omit map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1691    return mapData
1692   
1693def Fourier4DMap(data,reflDict):
1694    '''default doc string
1695   
1696    :param type name: description
1697   
1698    :returns: type name: description
1699   
1700    '''
1701    generalData = data['General']
1702    mapData = generalData['4DmapData']
1703    dmin = mapData['Resolution']
1704    SGData = generalData['SGData']
1705    SSGData = generalData['SSGData']
1706    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
1707    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
1708    cell = generalData['Cell'][1:8]       
1709    A = G2lat.cell2A(cell[:6])
1710    maxM = generalData['SuperVec'][2]
1711    Hmax = G2lat.getHKLmax(dmin,SGData,A)+[maxM,]
1712    adjHKLmax(SGData,Hmax)
1713    Hmax = np.asarray(Hmax,dtype='i')+1
1714    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1715    time0 = time.time()
1716    for iref,ref in enumerate(reflDict['RefList']):
1717        if ref[5] > dmin:
1718            Fosq,Fcsq,ph = ref[9:12]
1719            Uniq = np.inner(ref[:4],SSGMT)
1720            Phi = np.inner(ref[:4],SSGT)
1721            for i,hkl in enumerate(Uniq):        #uses uniq
1722                hkl = np.asarray(hkl,dtype='i')
1723                dp = 360.*Phi[i]                #and phi
1724                a = cosd(ph+dp)
1725                b = sind(ph+dp)
1726                phasep = complex(a,b)
1727                phasem = complex(a,-b)
1728                if 'Fobs' in mapData['MapType']:
1729                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1730                    h,k,l,m = hkl+Hmax
1731                    Fhkl[h,k,l,m] = F*phasep
1732                    h,k,l,m = -hkl+Hmax
1733                    Fhkl[h,k,l,m] = F*phasem
1734                elif 'delt-F' in mapData['MapType']:
1735                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1736                    h,k,l,m = hkl+Hmax
1737                    Fhkl[h,k,l,m] = dF*phasep
1738                    h,k,l,m = -hkl+Hmax
1739                    Fhkl[h,k,l,m] = dF*phasem
1740    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1741    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1742    mapData['Type'] = reflDict['Type']
1743    mapData['rho'] = np.real(rho)
1744    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1745    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
1746    return mapData
1747
1748def FourierMap(data,reflDict):
1749    '''default doc string
1750   
1751    :param type name: description
1752   
1753    :returns: type name: description
1754   
1755    '''
1756    generalData = data['General']
1757    mapData = generalData['Map']
1758    dmin = mapData['Resolution']
1759    SGData = generalData['SGData']
1760    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1761    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1762    cell = generalData['Cell'][1:8]       
1763    A = G2lat.cell2A(cell[:6])
1764    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1765    adjHKLmax(SGData,Hmax)
1766    Fhkl = np.zeros(shape=2*Hmax,dtype='c16')
1767#    Fhkl[0,0,0] = generalData['F000X']
1768    time0 = time.time()
1769    for iref,ref in enumerate(reflDict['RefList']):
1770        if ref[4] > dmin:
1771            Fosq,Fcsq,ph = ref[8:11]
1772            Uniq = np.inner(ref[:3],SGMT)
1773            Phi = np.inner(ref[:3],SGT)
1774            for i,hkl in enumerate(Uniq):        #uses uniq
1775                hkl = np.asarray(hkl,dtype='i')
1776                dp = 360.*Phi[i]                #and phi
1777                a = cosd(ph+dp)
1778                b = sind(ph+dp)
1779                phasep = complex(a,b)
1780                phasem = complex(a,-b)
1781                if 'Fobs' in mapData['MapType']:
1782                    F = np.where(Fosq>0.,np.sqrt(Fosq),0.)
1783                    h,k,l = hkl+Hmax
1784                    Fhkl[h,k,l] = F*phasep
1785                    h,k,l = -hkl+Hmax
1786                    Fhkl[h,k,l] = F*phasem
1787                elif 'Fcalc' in mapData['MapType']:
1788                    F = np.sqrt(Fcsq)
1789                    h,k,l = hkl+Hmax
1790                    Fhkl[h,k,l] = F*phasep
1791                    h,k,l = -hkl+Hmax
1792                    Fhkl[h,k,l] = F*phasem
1793                elif 'delt-F' in mapData['MapType']:
1794                    dF = np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1795                    h,k,l = hkl+Hmax
1796                    Fhkl[h,k,l] = dF*phasep
1797                    h,k,l = -hkl+Hmax
1798                    Fhkl[h,k,l] = dF*phasem
1799                elif '2*Fo-Fc' in mapData['MapType']:
1800                    F = 2.*np.where(Fosq>0.,np.sqrt(Fosq),0.)-np.sqrt(Fcsq)
1801                    h,k,l = hkl+Hmax
1802                    Fhkl[h,k,l] = F*phasep
1803                    h,k,l = -hkl+Hmax
1804                    Fhkl[h,k,l] = F*phasem
1805                elif 'Patterson' in mapData['MapType']:
1806                    h,k,l = hkl+Hmax
1807                    Fhkl[h,k,l] = complex(Fosq,0.)
1808                    h,k,l = -hkl+Hmax
1809                    Fhkl[h,k,l] = complex(Fosq,0.)
1810    rho = fft.fftn(fft.fftshift(Fhkl))/cell[6]
1811    print 'Fourier map time: %.4f'%(time.time()-time0),'no. elements: %d'%(Fhkl.size)
1812    mapData['Type'] = reflDict['Type']
1813    mapData['rho'] = np.real(rho)
1814    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
1815    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
1816    return mapData
1817   
1818# map printing for testing purposes
1819def printRho(SGLaue,rho,rhoMax):                         
1820    '''default doc string
1821   
1822    :param type name: description
1823   
1824    :returns: type name: description
1825   
1826    '''
1827    dim = len(rho.shape)
1828    if dim == 2:
1829        ix,jy = rho.shape
1830        for j in range(jy):
1831            line = ''
1832            if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1833                line += (jy-j)*'  '
1834            for i in range(ix):
1835                r = int(100*rho[i,j]/rhoMax)
1836                line += '%4d'%(r)
1837            print line+'\n'
1838    else:
1839        ix,jy,kz = rho.shape
1840        for k in range(kz):
1841            print 'k = ',k
1842            for j in range(jy):
1843                line = ''
1844                if SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1845                    line += (jy-j)*'  '
1846                for i in range(ix):
1847                    r = int(100*rho[i,j,k]/rhoMax)
1848                    line += '%4d'%(r)
1849                print line+'\n'
1850## keep this
1851               
1852def findOffset(SGData,A,Fhkl):   
1853    '''default doc string
1854   
1855    :param type name: description
1856   
1857    :returns: type name: description
1858   
1859    '''
1860    if SGData['SpGrp'] == 'P 1':
1861        return [0,0,0]   
1862    hklShape = Fhkl.shape
1863    hklHalf = np.array(hklShape)/2
1864    sortHKL = np.argsort(Fhkl.flatten())
1865    Fdict = {}
1866    for hkl in sortHKL:
1867        HKL = np.unravel_index(hkl,hklShape)
1868        F = Fhkl[HKL[0]][HKL[1]][HKL[2]]
1869        if F == 0.:
1870            break
1871        Fdict['%.6f'%(np.absolute(F))] = hkl
1872    Flist = np.flipud(np.sort(Fdict.keys()))
1873    F = str(1.e6)
1874    i = 0
1875    DH = []
1876    Dphi = []
1877    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1878    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1879    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
1880    for F in Flist:
1881        hkl = np.unravel_index(Fdict[F],hklShape)
1882        if np.any(np.abs(hkl-hklHalf)-Hmax > 0):
1883            continue
1884        iabsnt,mulp,Uniq,Phi = G2spc.GenHKLf(list(hkl-hklHalf),SGData)
1885        Uniq = np.array(Uniq,dtype='i')
1886        Phi = np.array(Phi)
1887        Uniq = np.concatenate((Uniq,-Uniq))+hklHalf         # put in Friedel pairs & make as index to Farray
1888        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
1889        Fh0 = Fhkl[hkl[0],hkl[1],hkl[2]]
1890        ang0 = np.angle(Fh0,deg=True)/360.
1891        for H,phi in zip(Uniq,Phi)[1:]:
1892            ang = (np.angle(Fhkl[H[0],H[1],H[2]],deg=True)/360.-phi)
1893            dH = H-hkl
1894            dang = ang-ang0
1895            DH.append(dH)
1896            Dphi.append((dang+.5) % 1.0)
1897        if i > 20 or len(DH) > 30:
1898            break
1899        i += 1
1900    DH = np.array(DH)
1901    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
1902    Dphi = np.array(Dphi)
1903    steps = np.array(hklShape)
1904    X,Y,Z = np.mgrid[0:1:1./steps[0],0:1:1./steps[1],0:1:1./steps[2]]
1905    XYZ = np.array(zip(X.flatten(),Y.flatten(),Z.flatten()))
1906    Dang = (np.dot(XYZ,DH.T)+.5)%1.-Dphi
1907    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
1908    hist,bins = np.histogram(Mmap,bins=1000)
1909#    for i,item in enumerate(hist[:10]):
1910#        print item,bins[i]
1911    chisq = np.min(Mmap)
1912    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
1913    print ' map offset chi**2: %.3f, map offset: %d %d %d'%(chisq,DX[0],DX[1],DX[2])
1914#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
1915    return DX
1916   
1917def ChargeFlip(data,reflDict,pgbar):
1918    '''default doc string
1919   
1920    :param type name: description
1921   
1922    :returns: type name: description
1923   
1924    '''
1925    generalData = data['General']
1926    mapData = generalData['Map']
1927    flipData = generalData['Flip']
1928    FFtable = {}
1929    if 'None' not in flipData['Norm element']:
1930        normElem = flipData['Norm element'].upper()
1931        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
1932        for ff in FFs:
1933            if ff['Symbol'] == normElem:
1934                FFtable.update(ff)
1935    dmin = flipData['Resolution']
1936    SGData = generalData['SGData']
1937    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
1938    SGT = np.array([ops[1] for ops in SGData['SGOps']])
1939    cell = generalData['Cell'][1:8]       
1940    A = G2lat.cell2A(cell[:6])
1941    Vol = cell[6]
1942    im = 0
1943    if generalData['Type'] in ['modulated','magnetic',]:
1944        im = 1
1945    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A),dtype='i')+1
1946    adjHKLmax(SGData,Hmax)
1947    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
1948    time0 = time.time()
1949    for iref,ref in enumerate(reflDict['RefList']):
1950        dsp = ref[4+im]
1951        if im and ref[3]:   #skip super lattice reflections - result is 3D projection
1952            continue
1953        if dsp > dmin:
1954            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
1955            if FFtable:
1956                SQ = 0.25/dsp**2
1957                ff *= G2el.ScatFac(FFtable,SQ)[0]
1958            if ref[8+im] > 0.:         #use only +ve Fobs**2
1959                E = np.sqrt(ref[8+im])/ff
1960            else:
1961                E = 0.
1962            ph = ref[10]
1963            ph = rn.uniform(0.,360.)
1964            Uniq = np.inner(ref[:3],SGMT)
1965            Phi = np.inner(ref[:3],SGT)
1966            for i,hkl in enumerate(Uniq):        #uses uniq
1967                hkl = np.asarray(hkl,dtype='i')
1968                dp = 360.*Phi[i]                #and phi
1969                a = cosd(ph+dp)
1970                b = sind(ph+dp)
1971                phasep = complex(a,b)
1972                phasem = complex(a,-b)
1973                h,k,l = hkl+Hmax
1974                Ehkl[h,k,l] = E*phasep
1975                h,k,l = -hkl+Hmax       #Friedel pair refl.
1976                Ehkl[h,k,l] = E*phasem
1977#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
1978    CEhkl = copy.copy(Ehkl)
1979    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
1980    Emask = ma.getmask(MEhkl)
1981    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
1982    Ncyc = 0
1983    old = np.seterr(all='raise')
1984    while True:       
1985        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
1986        CEsig = np.std(CErho)
1987        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
1988        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
1989        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
1990        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
1991        phase = CFhkl/np.absolute(CFhkl)
1992        CEhkl = np.absolute(Ehkl)*phase
1993        Ncyc += 1
1994        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
1995        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
1996        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
1997        if Rcf < 5.:
1998            break
1999        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2000        if not GoOn or Ncyc > 10000:
2001            break
2002    np.seterr(**old)
2003    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2004    CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))
2005    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2006    roll = findOffset(SGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2007       
2008    mapData['Rcf'] = Rcf
2009    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2010    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2011    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2012    mapData['Type'] = reflDict['Type']
2013    return mapData
2014   
2015def findSSOffset(SGData,SSGData,A,Fhklm):   
2016    '''default doc string
2017   
2018    :param type name: description
2019   
2020    :returns: type name: description
2021   
2022    '''
2023    if SGData['SpGrp'] == 'P 1':
2024        return [0,0,0,0]   
2025    hklmShape = Fhklm.shape
2026    hklmHalf = np.array(hklmShape)/2
2027    sortHKLM = np.argsort(Fhklm.flatten())
2028    Fdict = {}
2029    for hklm in sortHKLM:
2030        HKLM = np.unravel_index(hklm,hklmShape)
2031        F = Fhklm[HKLM[0]][HKLM[1]][HKLM[2]][HKLM[3]]
2032        if F == 0.:
2033            break
2034        Fdict['%.6f'%(np.absolute(F))] = hklm
2035    Flist = np.flipud(np.sort(Fdict.keys()))
2036    F = str(1.e6)
2037    i = 0
2038    DH = []
2039    Dphi = []
2040    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2041    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2042    Hmax = 2*np.asarray(G2lat.getHKLmax(3.5,SGData,A),dtype='i')
2043    for F in Flist:
2044        hklm = np.unravel_index(Fdict[F],hklmShape)
2045        if np.any(np.abs(hklm-hklmHalf)[:3]-Hmax > 0):
2046            continue
2047        Uniq = np.inner(hklm-hklmHalf,SSGMT)
2048        Phi = np.inner(hklm-hklmHalf,SSGT)
2049        Uniq = np.concatenate((Uniq,-Uniq))+hklmHalf         # put in Friedel pairs & make as index to Farray
2050        Phi = np.concatenate((Phi,-Phi))                      # and their phase shifts
2051        Fh0 = Fhklm[hklm[0],hklm[1],hklm[2],hklm[3]]
2052        ang0 = np.angle(Fh0,deg=True)/360.
2053        for H,phi in zip(Uniq,Phi)[1:]:
2054            ang = (np.angle(Fhklm[H[0],H[1],H[2],H[3]],deg=True)/360.-phi)
2055            dH = H-hklm
2056            dang = ang-ang0
2057            DH.append(dH)
2058            Dphi.append((dang+.5) % 1.0)
2059        if i > 20 or len(DH) > 30:
2060            break
2061        i += 1
2062    DH = np.array(DH)
2063    print ' map offset no.of terms: %d from %d reflections'%(len(DH),len(Flist))
2064    Dphi = np.array(Dphi)
2065    steps = np.array(hklmShape)
2066    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]]
2067    XYZT = np.array(zip(X.flatten(),Y.flatten(),Z.flatten(),T.flatten()))
2068    Dang = (np.dot(XYZT,DH.T)+.5)%1.-Dphi
2069    Mmap = np.reshape(np.sum((Dang)**2,axis=1),newshape=steps)/len(DH)
2070    hist,bins = np.histogram(Mmap,bins=1000)
2071#    for i,item in enumerate(hist[:10]):
2072#        print item,bins[i]
2073    chisq = np.min(Mmap)
2074    DX = -np.array(np.unravel_index(np.argmin(Mmap),Mmap.shape))
2075    print ' map offset chi**2: %.3f, map offset: %d %d %d %d'%(chisq,DX[0],DX[1],DX[2],DX[3])
2076#    print (np.dot(DX,DH.T)+.5)%1.-Dphi
2077    return DX
2078   
2079def SSChargeFlip(data,reflDict,pgbar):
2080    '''default doc string
2081   
2082    :param type name: description
2083   
2084    :returns: type name: description
2085   
2086    '''
2087    generalData = data['General']
2088    mapData = generalData['Map']
2089    map4DData = {}
2090    flipData = generalData['Flip']
2091    FFtable = {}
2092    if 'None' not in flipData['Norm element']:
2093        normElem = flipData['Norm element'].upper()
2094        FFs = G2el.GetFormFactorCoeff(normElem.split('+')[0].split('-')[0])
2095        for ff in FFs:
2096            if ff['Symbol'] == normElem:
2097                FFtable.update(ff)
2098    dmin = flipData['Resolution']
2099    SGData = generalData['SGData']
2100    SSGData = generalData['SSGData']
2101    SGMT = np.array([ops[0].T for ops in SGData['SGOps']])
2102    SGT = np.array([ops[1] for ops in SGData['SGOps']])
2103    SSGMT = np.array([ops[0].T for ops in SSGData['SSGOps']])
2104    SSGT = np.array([ops[1] for ops in SSGData['SSGOps']])
2105    cell = generalData['Cell'][1:8]       
2106    A = G2lat.cell2A(cell[:6])
2107    Vol = cell[6]
2108    maxM = generalData['SuperVec'][2]
2109    Hmax = np.asarray(G2lat.getHKLmax(dmin,SGData,A)+[maxM,],dtype='i')+1
2110    adjHKLmax(SGData,Hmax)
2111    Ehkl = np.zeros(shape=2*Hmax,dtype='c16')       #2X64bits per complex no.
2112    time0 = time.time()
2113    for iref,ref in enumerate(reflDict['RefList']):
2114        dsp = ref[5]
2115        if dsp > dmin:
2116            ff = 0.1*Vol    #est. no. atoms for ~10A**3/atom
2117            if FFtable:
2118                SQ = 0.25/dsp**2
2119                ff *= G2el.ScatFac(FFtable,SQ)[0]
2120            if ref[9] > 0.:         #use only +ve Fobs**2
2121                E = np.sqrt(ref[9])/ff
2122            else:
2123                E = 0.
2124            ph = ref[11]
2125            ph = rn.uniform(0.,360.)
2126            Uniq = np.inner(ref[:4],SSGMT)
2127            Phi = np.inner(ref[:4],SSGT)
2128            for i,hklm in enumerate(Uniq):        #uses uniq
2129                hklm = np.asarray(hklm,dtype='i')
2130                dp = 360.*Phi[i]                #and phi
2131                a = cosd(ph+dp)
2132                b = sind(ph+dp)
2133                phasep = complex(a,b)
2134                phasem = complex(a,-b)
2135                h,k,l,m = hklm+Hmax
2136                Ehkl[h,k,l,m] = E*phasep
2137                h,k,l,m = -hklm+Hmax       #Friedel pair refl.
2138                Ehkl[h,k,l,m] = E*phasem
2139#    Ehkl[Hmax] = 0.00001           #this to preserve F[0,0,0]
2140    CEhkl = copy.copy(Ehkl)
2141    MEhkl = ma.array(Ehkl,mask=(Ehkl==0.0))
2142    Emask = ma.getmask(MEhkl)
2143    sumE = np.sum(ma.array(np.absolute(CEhkl),mask=Emask))
2144    Ncyc = 0
2145    old = np.seterr(all='raise')
2146    while True:       
2147        CErho = np.real(fft.fftn(fft.fftshift(CEhkl)))*(1.+0j)
2148        CEsig = np.std(CErho)
2149        CFrho = np.where(np.real(CErho) >= flipData['k-factor']*CEsig,CErho,-CErho)
2150        CFrho = np.where(np.real(CErho) <= flipData['k-Max']*CEsig,CFrho,-CFrho)      #solves U atom problem!
2151        CFhkl = fft.ifftshift(fft.ifftn(CFrho))
2152        CFhkl = np.where(CFhkl,CFhkl,1.0)           #avoid divide by zero
2153        phase = CFhkl/np.absolute(CFhkl)
2154        CEhkl = np.absolute(Ehkl)*phase
2155        Ncyc += 1
2156        sumCF = np.sum(ma.array(np.absolute(CFhkl),mask=Emask))
2157        DEhkl = np.absolute(np.absolute(Ehkl)/sumE-np.absolute(CFhkl)/sumCF)
2158        Rcf = min(100.,np.sum(ma.array(DEhkl,mask=Emask)*100.))
2159        if Rcf < 5.:
2160            break
2161        GoOn = pgbar.Update(Rcf,newmsg='%s%8.3f%s\n%s %d'%('Residual Rcf =',Rcf,'%','No.cycles = ',Ncyc))[0]
2162        if not GoOn or Ncyc > 10000:
2163            break
2164    np.seterr(**old)
2165    print ' Charge flip time: %.4f'%(time.time()-time0),'no. elements: %d'%(Ehkl.size)
2166    CErho = np.real(fft.fftn(fft.fftshift(CEhkl[:,:,:,maxM+1])))
2167    SSrho = np.real(fft.fftn(fft.fftshift(CEhkl)))
2168    print ' No.cycles = ',Ncyc,'Residual Rcf =%8.3f%s'%(Rcf,'%')+' Map size:',CErho.shape
2169    roll = findSSOffset(SGData,SSGData,A,CEhkl)               #CEhkl needs to be just the observed set, not the full set!
2170
2171    mapData['Rcf'] = Rcf
2172    mapData['rho'] = np.roll(np.roll(np.roll(CErho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2173    mapData['rhoMax'] = max(np.max(mapData['rho']),-np.min(mapData['rho']))
2174    mapData['minmax'] = [np.max(mapData['rho']),np.min(mapData['rho'])]
2175    mapData['Type'] = reflDict['Type']
2176
2177    map4DData['Rcf'] = Rcf
2178    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))
2179    map4DData['rhoMax'] = max(np.max(map4DData['rho']),-np.min(map4DData['rho']))
2180    map4DData['minmax'] = [np.max(map4DData['rho']),np.min(map4DData['rho'])]
2181    map4DData['Type'] = reflDict['Type']
2182    return mapData,map4DData
2183   
2184def SearchMap(generalData,drawingData,Neg=False):
2185    '''Does a search of a density map for peaks meeting the criterion of peak
2186    height is greater than mapData['cutOff']/100 of mapData['rhoMax'] where
2187    mapData is data['General']['mapData']; the map is also in mapData.
2188
2189    :param generalData: the phase data structure; includes the map
2190    :param drawingData: the drawing data structure
2191    :param Neg:  if True then search for negative peaks (i.e. H-atoms & neutron data)
2192
2193    :returns: (peaks,mags,dzeros) where
2194
2195        * peaks : ndarray
2196          x,y,z positions of the peaks found in the map
2197        * mags : ndarray
2198          the magnitudes of the peaks
2199        * dzeros : ndarray
2200          the distance of the peaks from  the unit cell origin
2201        * dcent : ndarray
2202          the distance of the peaks from  the unit cell center
2203
2204    '''       
2205    rollMap = lambda rho,roll: np.roll(np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1),roll[2],axis=2)
2206   
2207    norm = 1./(np.sqrt(3.)*np.sqrt(2.*np.pi)**3)
2208   
2209#    def noDuplicate(xyz,peaks,Amat):
2210#        XYZ = np.inner(Amat,xyz)
2211#        if True in [np.allclose(XYZ,np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2212#            print ' Peak',xyz,' <0.5A from another peak'
2213#            return False
2214#        return True
2215#                           
2216    def fixSpecialPos(xyz,SGData,Amat):
2217        equivs = G2spc.GenAtom(xyz,SGData,Move=True)
2218        X = []
2219        xyzs = [equiv[0] for equiv in equivs]
2220        for x in xyzs:
2221            if np.sqrt(np.sum(np.inner(Amat,xyz-x)**2,axis=0))<0.5:
2222                X.append(x)
2223        if len(X) > 1:
2224            return np.average(X,axis=0)
2225        else:
2226            return xyz
2227       
2228    def rhoCalc(parms,rX,rY,rZ,res,SGLaue):
2229        Mag,x0,y0,z0,sig = parms
2230        z = -((x0-rX)**2+(y0-rY)**2+(z0-rZ)**2)/(2.*sig**2)
2231#        return norm*Mag*np.exp(z)/(sig*res**3)     #not slower but some faults in LS
2232        return norm*Mag*(1.+z+z**2/2.)/(sig*res**3)
2233       
2234    def peakFunc(parms,rX,rY,rZ,rho,res,SGLaue):
2235        Mag,x0,y0,z0,sig = parms
2236        M = rho-rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2237        return M
2238       
2239    def peakHess(parms,rX,rY,rZ,rho,res,SGLaue):
2240        Mag,x0,y0,z0,sig = parms
2241        dMdv = np.zeros(([5,]+list(rX.shape)))
2242        delt = .01
2243        for i in range(5):
2244            parms[i] -= delt
2245            rhoCm = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2246            parms[i] += 2.*delt
2247            rhoCp = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2248            parms[i] -= delt
2249            dMdv[i] = (rhoCp-rhoCm)/(2.*delt)
2250        rhoC = rhoCalc(parms,rX,rY,rZ,res,SGLaue)
2251        Vec = np.sum(np.sum(np.sum(dMdv*(rho-rhoC),axis=3),axis=2),axis=1)
2252        dMdv = np.reshape(dMdv,(5,rX.size))
2253        Hess = np.inner(dMdv,dMdv)
2254       
2255        return Vec,Hess
2256       
2257    phaseName = generalData['Name']
2258    SGData = generalData['SGData']
2259    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2260    peaks = []
2261    mags = []
2262    dzeros = []
2263    dcent = []
2264    try:
2265        mapData = generalData['Map']
2266        contLevel = mapData['cutOff']*mapData['rhoMax']/100.
2267        if Neg:
2268            rho = -copy.copy(mapData['rho'])     #flip +/-
2269        else:
2270            rho = copy.copy(mapData['rho'])     #don't mess up original
2271        mapHalf = np.array(rho.shape)/2
2272        res = mapData['Resolution']
2273        incre = np.array(rho.shape,dtype=np.float)
2274        step = max(1.0,1./res)+1
2275        steps = np.array(3*[step,])
2276    except KeyError:
2277        print '**** ERROR - Fourier map not defined'
2278        return peaks,mags
2279    rhoMask = ma.array(rho,mask=(rho<contLevel))
2280    indices = (-1,0,1)
2281    rolls = np.array([[h,k,l] for h in indices for k in indices for l in indices])
2282    for roll in rolls:
2283        if np.any(roll):
2284            rhoMask = ma.array(rhoMask,mask=(rhoMask-rollMap(rho,roll)<=0.))
2285    indx = np.transpose(rhoMask.nonzero())
2286    peaks = indx/incre
2287    mags = rhoMask[rhoMask.nonzero()]
2288    for i,[ind,peak,mag] in enumerate(zip(indx,peaks,mags)):
2289        rho = rollMap(rho,ind)
2290        rMM = mapHalf-steps
2291        rMP = mapHalf+steps+1
2292        rhoPeak = rho[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2293        peakInt = np.sum(rhoPeak)*res**3
2294        rX,rY,rZ = np.mgrid[rMM[0]:rMP[0],rMM[1]:rMP[1],rMM[2]:rMP[2]]
2295        x0 = [peakInt,mapHalf[0],mapHalf[1],mapHalf[2],2.0]          #magnitude, position & width(sig)
2296        result = HessianLSQ(peakFunc,x0,Hess=peakHess,
2297            args=(rX,rY,rZ,rhoPeak,res,SGData['SGLaue']),ftol=.01,maxcyc=10)
2298        x1 = result[0]
2299        if not np.any(x1 < 0):
2300            mag = x1[0]
2301            peak = (np.array(x1[1:4])-ind)/incre
2302        peak = fixSpecialPos(peak,SGData,Amat)
2303        rho = rollMap(rho,-ind)
2304    cent = np.ones(3)*.5     
2305    dzeros = np.sqrt(np.sum(np.inner(Amat,peaks)**2,axis=0))
2306    dcent = np.sqrt(np.sum(np.inner(Amat,peaks-cent)**2,axis=0))
2307    if Neg:     #want negative magnitudes for negative peaks
2308        return np.array(peaks),-np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
2309    else:
2310        return np.array(peaks),np.array([mags,]).T,np.array([dzeros,]).T,np.array([dcent,]).T
2311   
2312def sortArray(data,pos,reverse=False):
2313    '''data is a list of items
2314    sort by pos in list; reverse if True
2315    '''
2316    T = []
2317    for i,M in enumerate(data):
2318        try:
2319            T.append((M[pos],i))
2320        except IndexError:
2321            return data
2322    D = dict(zip(T,data))
2323    T.sort()
2324    if reverse:
2325        T.reverse()
2326    X = []
2327    for key in T:
2328        X.append(D[key])
2329    return X
2330
2331def PeaksEquiv(data,Ind):
2332    '''Find the equivalent map peaks for those selected. Works on the
2333    contents of data['Map Peaks'].
2334
2335    :param data: the phase data structure
2336    :param list Ind: list of selected peak indices
2337    :returns: augmented list of peaks including those related by symmetry to the
2338      ones in Ind
2339
2340    '''       
2341    def Duplicate(xyz,peaks,Amat):
2342        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2343            return True
2344        return False
2345                           
2346    generalData = data['General']
2347    cell = generalData['Cell'][1:7]
2348    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2349    A = G2lat.cell2A(cell)
2350    SGData = generalData['SGData']
2351    mapPeaks = data['Map Peaks']
2352    XYZ = np.array([xyz[1:4] for xyz in mapPeaks])
2353    Indx = {}
2354    for ind in Ind:
2355        xyz = np.array(mapPeaks[ind][1:4])
2356        xyzs = np.array([equiv[0] for equiv in G2spc.GenAtom(xyz,SGData,Move=True)])
2357        for jnd,xyz in enumerate(XYZ):       
2358            Indx[jnd] = Duplicate(xyz,xyzs,Amat)
2359    Ind = []
2360    for ind in Indx:
2361        if Indx[ind]:
2362            Ind.append(ind)
2363    return Ind
2364               
2365def PeaksUnique(data,Ind):
2366    '''Finds the symmetry unique set of peaks from those selected. Works on the
2367    contents of data['Map Peaks'].
2368
2369    :param data: the phase data structure
2370    :param list Ind: list of selected peak indices
2371    :returns: the list of symmetry unique peaks from among those given in Ind
2372
2373    '''       
2374#    XYZE = np.array([[equiv[0] for equiv in G2spc.GenAtom(xyz[1:4],SGData,Move=True)] for xyz in mapPeaks]) #keep this!!
2375
2376    def noDuplicate(xyz,peaks,Amat):
2377        if True in [np.allclose(np.inner(Amat,xyz),np.inner(Amat,peak),atol=0.5) for peak in peaks]:
2378            return False
2379        return True
2380                           
2381    generalData = data['General']
2382    cell = generalData['Cell'][1:7]
2383    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
2384    A = G2lat.cell2A(cell)
2385    SGData = generalData['SGData']
2386    mapPeaks = data['Map Peaks']
2387    Indx = {}
2388    XYZ = {}
2389    for ind in Ind:
2390        XYZ[ind] = np.array(mapPeaks[ind][1:4])
2391        Indx[ind] = True
2392    for ind in Ind:
2393        if Indx[ind]:
2394            xyz = XYZ[ind]
2395            for jnd in Ind:
2396                if ind != jnd and Indx[jnd]:                       
2397                    Equiv = G2spc.GenAtom(XYZ[jnd],SGData,Move=True)
2398                    xyzs = np.array([equiv[0] for equiv in Equiv])
2399                    Indx[jnd] = noDuplicate(xyz,xyzs,Amat)
2400    Ind = []
2401    for ind in Indx:
2402        if Indx[ind]:
2403            Ind.append(ind)
2404    return Ind
2405   
2406################################################################################
2407##### single peak fitting profile fxn stuff
2408################################################################################
2409
2410def getCWsig(ins,pos):
2411    '''get CW peak profile sigma
2412   
2413    :param dict ins: instrument parameters with at least 'U', 'V', & 'W'
2414      as values only
2415    :param float pos: 2-theta of peak
2416    :returns: float getCWsig: peak sigma
2417   
2418    '''
2419    tp = tand(pos/2.0)
2420    return ins['U']*tp**2+ins['V']*tp+ins['W']
2421   
2422def getCWsigDeriv(pos):
2423    '''get derivatives of CW peak profile sigma wrt U,V, & W
2424   
2425    :param float pos: 2-theta of peak
2426   
2427    :returns: list getCWsigDeriv: d(sig)/dU, d(sig)/dV & d(sig)/dW
2428   
2429    '''
2430    tp = tand(pos/2.0)
2431    return tp**2,tp,1.0
2432   
2433def getCWgam(ins,pos):
2434    '''get CW peak profile gamma
2435   
2436    :param dict ins: instrument parameters with at least 'X' & 'Y'
2437      as values only
2438    :param float pos: 2-theta of peak
2439    :returns: float getCWgam: peak gamma
2440   
2441    '''
2442    return ins['X']/cosd(pos/2.0)+ins['Y']*tand(pos/2.0)
2443   
2444def getCWgamDeriv(pos):
2445    '''get derivatives of CW peak profile gamma wrt X & Y
2446   
2447    :param float pos: 2-theta of peak
2448   
2449    :returns: list getCWgamDeriv: d(gam)/dX & d(gam)/dY
2450   
2451    '''
2452    return 1./cosd(pos/2.0),tand(pos/2.0)
2453   
2454def getTOFsig(ins,dsp):
2455    '''get TOF peak profile sigma
2456   
2457    :param dict ins: instrument parameters with at least 'sig-0', 'sig-1' & 'sig-q'
2458      as values only
2459    :param float dsp: d-spacing of peak
2460   
2461    :returns: float getTOFsig: peak sigma
2462   
2463    '''
2464    return ins['sig-0']+ins['sig-1']*dsp**2+ins['sig-2']*dsp**4+ins['sig-q']/dsp**2
2465   
2466def getTOFsigDeriv(dsp):
2467    '''get derivatives of TOF peak profile gamma wrt sig-0, sig-1, & sig-q
2468   
2469    :param float dsp: d-spacing of peak
2470   
2471    :returns: list getTOFsigDeriv: d(sig0/d(sig-0), d(sig)/d(sig-1) & d(sig)/d(sig-q)
2472   
2473    '''
2474    return 1.0,dsp**2,dsp**4,1./dsp**2
2475   
2476def getTOFgamma(ins,dsp):
2477    '''get TOF peak profile gamma
2478   
2479    :param dict ins: instrument parameters with at least 'X' & 'Y'
2480      as values only
2481    :param float dsp: d-spacing of peak
2482   
2483    :returns: float getTOFgamma: peak gamma
2484   
2485    '''
2486    return ins['X']*dsp+ins['Y']*dsp**2
2487   
2488def getTOFgammaDeriv(dsp):
2489    '''get derivatives of TOF peak profile gamma wrt X & Y
2490   
2491    :param float dsp: d-spacing of peak
2492   
2493    :returns: list getTOFgammaDeriv: d(gam)/dX & d(gam)/dY
2494   
2495    '''
2496    return dsp,dsp**2
2497   
2498def getTOFbeta(ins,dsp):
2499    '''get TOF peak profile beta
2500   
2501    :param dict ins: instrument parameters with at least 'beat-0', 'beta-1' & 'beta-q'
2502      as values only
2503    :param float dsp: d-spacing of peak
2504   
2505    :returns: float getTOFbeta: peak beat
2506   
2507    '''
2508    return ins['beta-0']+ins['beta-1']/dsp**4+ins['beta-q']/dsp**2
2509   
2510def getTOFbetaDeriv(dsp):
2511    '''get derivatives of TOF peak profile beta wrt beta-0, beta-1, & beat-q
2512   
2513    :param float dsp: d-spacing of peak
2514   
2515    :returns: list getTOFbetaDeriv: d(beta)/d(beat-0), d(beta)/d(beta-1) & d(beta)/d(beta-q)
2516   
2517    '''
2518    return 1.0,1./dsp**4,1./dsp**2
2519   
2520def getTOFalpha(ins,dsp):
2521    '''get TOF peak profile alpha
2522   
2523    :param dict ins: instrument parameters with at least 'alpha'
2524      as values only
2525    :param float dsp: d-spacing of peak
2526   
2527    :returns: flaot getTOFalpha: peak alpha
2528   
2529    '''
2530    return ins['alpha']/dsp
2531   
2532def getTOFalphaDeriv(dsp):
2533    '''get derivatives of TOF peak profile beta wrt alpha
2534   
2535    :param float dsp: d-spacing of peak
2536   
2537    :returns: float getTOFalphaDeriv: d(alp)/d(alpha)
2538   
2539    '''
2540    return 1./dsp
2541   
2542def setPeakparms(Parms,Parms2,pos,mag,ifQ=False,useFit=False):
2543    '''set starting peak parameters for single peak fits from plot selection or auto selection
2544   
2545    :param dict Parms: instrument parameters dictionary
2546    :param dict Parms2: table lookup for TOF profile coefficients
2547    :param float pos: peak position in 2-theta, TOF or Q (ifQ=True)
2548    :param float mag: peak top magnitude from pick
2549    :param bool ifQ: True if pos in Q
2550    :param bool useFit: True if use fitted CW Parms values (not defaults)
2551   
2552    :returns: list XY: peak list entry:
2553        for CW: [pos,0,mag,1,sig,0,gam,0]
2554        for TOF: [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2555        NB: mag refinement set by default, all others off
2556   
2557    '''
2558    ind = 0
2559    if useFit:
2560        ind = 1
2561    ins = {}
2562    if 'C' in Parms['Type'][0]:                            #CW data - TOF later in an elif
2563        for x in ['U','V','W','X','Y']:
2564            ins[x] = Parms[x][ind]
2565        if ifQ:                              #qplot - convert back to 2-theta
2566            pos = 2.0*asind(pos*wave/(4*math.pi))
2567        sig = getCWsig(ins,pos)
2568        gam = getCWgam(ins,pos)           
2569        XY = [pos,0, mag,1, sig,0, gam,0]       #default refine intensity 1st
2570    else:
2571        if ifQ:
2572            dsp = 2.*np.pi/pos
2573            pos = Parms['difC']*dsp
2574        else:
2575            dsp = pos/Parms['difC'][1]
2576        if 'Pdabc' in Parms2:
2577            for x in ['sig-0','sig-1','sig-2','sig-q','X','Y']:
2578                ins[x] = Parms[x][ind]
2579            Pdabc = Parms2['Pdabc'].T
2580            alp = np.interp(dsp,Pdabc[0],Pdabc[1])
2581            bet = np.interp(dsp,Pdabc[0],Pdabc[2])
2582        else:
2583            for x in ['alpha','beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','X','Y']:
2584                ins[x] = Parms[x][ind]
2585            alp = getTOFalpha(ins,dsp)
2586            bet = getTOFbeta(ins,dsp)
2587        sig = getTOFsig(ins,dsp)
2588        gam = getTOFgamma(ins,dsp)
2589        XY = [pos,0,mag,1,alp,0,bet,0,sig,0,gam,0]
2590    return XY
2591   
2592################################################################################
2593##### MC/SA stuff
2594################################################################################
2595
2596#scipy/optimize/anneal.py code modified by R. Von Dreele 2013
2597# Original Author: Travis Oliphant 2002
2598# Bug-fixes in 2006 by Tim Leslie
2599
2600
2601import numpy
2602from numpy import asarray, tan, exp, ones, squeeze, sign, \
2603     all, log, sqrt, pi, shape, array, minimum, where
2604from numpy import random
2605
2606#__all__ = ['anneal']
2607
2608_double_min = numpy.finfo(float).min
2609_double_max = numpy.finfo(float).max
2610class base_schedule(object):
2611    def __init__(self):
2612        self.dwell = 20
2613        self.learn_rate = 0.5
2614        self.lower = -10
2615        self.upper = 10
2616        self.Ninit = 50
2617        self.accepted = 0
2618        self.tests = 0
2619        self.feval = 0
2620        self.k = 0
2621        self.T = None
2622
2623    def init(self, **options):
2624        self.__dict__.update(options)
2625        self.lower = asarray(self.lower)
2626        self.lower = where(self.lower == numpy.NINF, -_double_max, self.lower)
2627        self.upper = asarray(self.upper)
2628        self.upper = where(self.upper == numpy.PINF, _double_max, self.upper)
2629        self.k = 0
2630        self.accepted = 0
2631        self.feval = 0
2632        self.tests = 0
2633
2634    def getstart_temp(self, best_state):
2635        """ Find a matching starting temperature and starting parameters vector
2636        i.e. find x0 such that func(x0) = T0.
2637
2638        Parameters
2639        ----------
2640        best_state : _state
2641            A _state object to store the function value and x0 found.
2642
2643        returns
2644        -------
2645        x0 : array
2646            The starting parameters vector.
2647        """
2648
2649        assert(not self.dims is None)
2650        lrange = self.lower
2651        urange = self.upper
2652        fmax = _double_min
2653        fmin = _double_max
2654        for _ in range(self.Ninit):
2655            x0 = random.uniform(size=self.dims)*(urange-lrange) + lrange
2656            fval = self.func(x0, *self.args)
2657            self.feval += 1
2658            if fval > fmax:
2659                fmax = fval
2660            if fval < fmin:
2661                fmin = fval
2662                best_state.cost = fval
2663                best_state.x = array(x0)
2664
2665        self.T0 = (fmax-fmin)*1.5
2666        return best_state.x
2667       
2668    def set_range(self,x0,frac):
2669        delrange = frac*(self.upper-self.lower)
2670        self.upper = x0+delrange
2671        self.lower = x0-delrange
2672
2673    def accept_test(self, dE):
2674        T = self.T
2675        self.tests += 1
2676        if dE < 0:
2677            self.accepted += 1
2678            return 1
2679        p = exp(-dE*1.0/self.boltzmann/T)
2680        if (p > random.uniform(0.0, 1.0)):
2681            self.accepted += 1
2682            return 1
2683        return 0
2684
2685    def update_guess(self, x0):
2686        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2687
2688    def update_temp(self, x0):
2689        pass
2690
2691
2692#  A schedule due to Lester Ingber modified to use bounds - OK
2693class fast_sa(base_schedule):
2694    def init(self, **options):
2695        self.__dict__.update(options)
2696        if self.m is None:
2697            self.m = 1.0
2698        if self.n is None:
2699            self.n = 1.0
2700        self.c = self.m * exp(-self.n * self.quench)
2701
2702    def update_guess(self, x0):
2703        x0 = asarray(x0)
2704        u = squeeze(random.uniform(0.0, 1.0, size=self.dims))
2705        T = self.T
2706        xc = (sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)+1.0)/2.0
2707        xnew = xc*(self.upper - self.lower)+self.lower
2708        return xnew
2709#        y = sign(u-0.5)*T*((1+1.0/T)**abs(2*u-1)-1.0)
2710#        xc = y*(self.upper - self.lower)
2711#        xnew = x0 + xc
2712#        return xnew
2713
2714    def update_temp(self):
2715        self.T = self.T0*exp(-self.c * self.k**(self.quench))
2716        self.k += 1
2717        return
2718
2719class cauchy_sa(base_schedule):     #modified to use bounds - not good
2720    def update_guess(self, x0):
2721        x0 = asarray(x0)
2722        numbers = squeeze(random.uniform(-pi/4, pi/4, size=self.dims))
2723        xc = (1.+(self.learn_rate * self.T * tan(numbers))%1.)
2724        xnew = xc*(self.upper - self.lower)+self.lower
2725        return xnew
2726#        numbers = squeeze(random.uniform(-pi/2, pi/2, size=self.dims))
2727#        xc = self.learn_rate * self.T * tan(numbers)
2728#        xnew = x0 + xc
2729#        return xnew
2730
2731    def update_temp(self):
2732        self.T = self.T0/(1+self.k)
2733        self.k += 1
2734        return
2735
2736class boltzmann_sa(base_schedule):
2737#    def update_guess(self, x0):
2738#        std = minimum(sqrt(self.T)*ones(self.dims), (self.upper-self.lower)/3.0/self.learn_rate)
2739#        x0 = asarray(x0)
2740#        xc = squeeze(random.normal(0, 1.0, size=self.dims))
2741#
2742#        xnew = x0 + xc*std*self.learn_rate
2743#        return xnew
2744
2745    def update_temp(self):
2746        self.k += 1
2747        self.T = self.T0 / log(self.k+1.0)
2748        return
2749
2750class log_sa(base_schedule):        #OK
2751
2752    def init(self,**options):
2753        self.__dict__.update(options)
2754       
2755    def update_guess(self,x0):     #same as default
2756        return np.squeeze(np.random.uniform(0.,1.,size=self.dims))*(self.upper-self.lower)+self.lower
2757       
2758    def update_temp(self):
2759        self.k += 1
2760        self.T = self.T0*self.slope**self.k
2761       
2762class _state(object):
2763    def __init__(self):
2764        self.x = None
2765        self.cost = None
2766
2767# TODO:
2768#     allow for general annealing temperature profile
2769#     in that case use update given by alpha and omega and
2770#     variation of all previous updates and temperature?
2771
2772# Simulated annealing
2773
2774def anneal(func, x0, args=(), schedule='fast', full_output=0,
2775           T0=None, Tf=1e-12, maxeval=None, maxaccept=None, maxiter=400,
2776           boltzmann=1.0, learn_rate=0.5, feps=1e-6, quench=1.0, m=1.0, n=1.0,
2777           lower=-100, upper=100, dwell=50, slope=0.9,ranStart=False,
2778           ranRange=0.10,autoRan=False,dlg=None):
2779    """Minimize a function using simulated annealing.
2780
2781    Schedule is a schedule class implementing the annealing schedule.
2782    Available ones are 'fast', 'cauchy', 'boltzmann'
2783
2784    :param callable func: f(x, \*args)
2785        Function to be optimized.
2786    :param ndarray x0:
2787        Initial guess.
2788    :param tuple args:
2789        Extra parameters to `func`.
2790    :param base_schedule schedule:
2791        Annealing schedule to use (a class).
2792    :param bool full_output:
2793        Whether to return optional outputs.
2794    :param float T0:
2795        Initial Temperature (estimated as 1.2 times the largest
2796        cost-function deviation over random points in the range).
2797    :param float Tf:
2798        Final goal temperature.
2799    :param int maxeval:
2800        Maximum function evaluations.
2801    :param int maxaccept:
2802        Maximum changes to accept.
2803    :param int maxiter:
2804        Maximum cooling iterations.
2805    :param float learn_rate:
2806        Scale constant for adjusting guesses.
2807    :param float boltzmann:
2808        Boltzmann constant in acceptance test
2809        (increase for less stringent test at each temperature).
2810    :param float feps:
2811        Stopping relative error tolerance for the function value in
2812        last four coolings.
2813    :param float quench,m,n:
2814        Parameters to alter fast_sa schedule.
2815    :param float/ndarray lower,upper:
2816        Lower and upper bounds on `x`.
2817    :param int dwell:
2818        The number of times to search the space at each temperature.
2819    :param float slope:
2820        Parameter for log schedule
2821    :param bool ranStart=False:
2822        True for set 10% of ranges about x
2823
2824    :returns: (xmin, Jmin, T, feval, iters, accept, retval) where
2825
2826     * xmin (ndarray): Point giving smallest value found.
2827     * Jmin (float): Minimum value of function found.
2828     * T (float): Final temperature.
2829     * feval (int): Number of function evaluations.
2830     * iters (int): Number of cooling iterations.
2831     * accept (int): Number of tests accepted.
2832     * retval (int): Flag indicating stopping condition:
2833
2834              *  0: Points no longer changing
2835              *  1: Cooled to final temperature
2836              *  2: Maximum function evaluations
2837              *  3: Maximum cooling iterations reached
2838              *  4: Maximum accepted query locations reached
2839              *  5: Final point not the minimum amongst encountered points
2840
2841    *Notes*:
2842    Simulated annealing is a random algorithm which uses no derivative
2843    information from the function being optimized. In practice it has
2844    been more useful in discrete optimization than continuous
2845    optimization, as there are usually better algorithms for continuous
2846    optimization problems.
2847
2848    Some experimentation by trying the difference temperature
2849    schedules and altering their parameters is likely required to
2850    obtain good performance.
2851
2852    The randomness in the algorithm comes from random sampling in numpy.
2853    To obtain the same results you can call numpy.random.seed with the
2854    same seed immediately before calling scipy.optimize.anneal.
2855
2856    We give a brief description of how the three temperature schedules
2857    generate new points and vary their temperature. Temperatures are
2858    only updated with iterations in the outer loop. The inner loop is
2859    over xrange(dwell), and new points are generated for
2860    every iteration in the inner loop. (Though whether the proposed
2861    new points are accepted is probabilistic.)
2862
2863    For readability, let d denote the dimension of the inputs to func.
2864    Also, let x_old denote the previous state, and k denote the
2865    iteration number of the outer loop. All other variables not
2866    defined below are input variables to scipy.optimize.anneal itself.
2867
2868    In the 'fast' schedule the updates are ::
2869
2870        u ~ Uniform(0, 1, size=d)
2871        y = sgn(u - 0.5) * T * ((1+ 1/T)**abs(2u-1) -1.0)
2872        xc = y * (upper - lower)
2873        x_new = x_old + xc
2874
2875        c = n * exp(-n * quench)
2876        T_new = T0 * exp(-c * k**quench)
2877
2878
2879    In the 'cauchy' schedule the updates are ::
2880
2881        u ~ Uniform(-pi/2, pi/2, size=d)
2882        xc = learn_rate * T * tan(u)
2883        x_new = x_old + xc
2884
2885        T_new = T0 / (1+k)
2886
2887    In the 'boltzmann' schedule the updates are ::
2888
2889        std = minimum( sqrt(T) * ones(d), (upper-lower) / (3*learn_rate) )
2890        y ~ Normal(0, std, size=d)
2891        x_new = x_old + learn_rate * y
2892
2893        T_new = T0 / log(1+k)
2894
2895    """
2896    x0 = asarray(x0)
2897    lower = asarray(lower)
2898    upper = asarray(upper)
2899
2900    schedule = eval(schedule+'_sa()')
2901    #   initialize the schedule
2902    schedule.init(dims=shape(x0),func=func,args=args,boltzmann=boltzmann,T0=T0,
2903                  learn_rate=learn_rate, lower=lower, upper=upper,
2904                  m=m, n=n, quench=quench, dwell=dwell, slope=slope)
2905
2906    current_state, last_state, best_state = _state(), _state(), _state()
2907    if ranStart:
2908        schedule.set_range(x0,ranRange)
2909    if T0 is None:
2910        x0 = schedule.getstart_temp(best_state)
2911    else:
2912        x0 = random.uniform(size=len(x0))*(upper-lower) + lower
2913        best_state.x = None
2914        best_state.cost = numpy.Inf
2915
2916    last_state.x = asarray(x0).copy()
2917    fval = func(x0,*args)
2918    schedule.feval += 1
2919    last_state.cost = fval
2920    if last_state.cost < best_state.cost:
2921        best_state.cost = fval
2922        best_state.x = asarray(x0).copy()
2923    schedule.T = schedule.T0
2924    fqueue = [100, 300, 500, 700]
2925    iters = 1
2926    keepGoing = True
2927    bestn = 0
2928    while keepGoing:
2929        retval = 0
2930        for n in xrange(dwell):
2931            current_state.x = schedule.update_guess(last_state.x)
2932            current_state.cost = func(current_state.x,*args)
2933            schedule.feval += 1
2934
2935            dE = current_state.cost - last_state.cost
2936            if schedule.accept_test(dE):
2937                last_state.x = current_state.x.copy()
2938                last_state.cost = current_state.cost
2939                if last_state.cost < best_state.cost:
2940                    best_state.x = last_state.x.copy()
2941                    best_state.cost = last_state.cost
2942                    bestn = n
2943                    if best_state.cost < 1.0 and autoRan:
2944                        schedule.set_range(x0,best_state.cost/2.)                       
2945        if dlg:
2946            GoOn = dlg.Update(min(100.,best_state.cost*100),
2947                newmsg='%s%8.5f, %s%d\n%s%8.4f%s'%('Temperature =',schedule.T, \
2948                    'Best trial:',bestn,  \
2949                    'MC/SA Residual:',best_state.cost*100,'%', \
2950                    ))[0]
2951            if not GoOn:
2952                best_state.x = last_state.x.copy()
2953                best_state.cost = last_state.cost
2954                retval = 5
2955        schedule.update_temp()
2956        iters += 1
2957        # Stopping conditions
2958        # 0) last saved values of f from each cooling step
2959        #     are all very similar (effectively cooled)
2960        # 1) Tf is set and we are below it
2961        # 2) maxeval is set and we are past it
2962        # 3) maxiter is set and we are past it
2963        # 4) maxaccept is set and we are past it
2964        # 5) user canceled run via progress bar
2965
2966        fqueue.append(squeeze(last_state.cost))
2967        fqueue.pop(0)
2968        af = asarray(fqueue)*1.0
2969        if retval == 5:
2970            print ' User terminated run; incomplete MC/SA'
2971            keepGoing = False
2972            break
2973        if all(abs((af-af[0])/af[0]) < feps):
2974            retval = 0
2975            if abs(af[-1]-best_state.cost) > feps*10:
2976                retval = 5
2977#                print "Warning: Cooled to %f at %s but this is not" \
2978#                      % (squeeze(last_state.cost), str(squeeze(last_state.x))) \
2979#                      + " the smallest point found."
2980            break
2981        if (Tf is not None) and (schedule.T < Tf):
2982            retval = 1
2983            break
2984        if (maxeval is not None) and (schedule.feval > maxeval):
2985            retval = 2
2986            break
2987        if (iters > maxiter):
2988            print "Warning: Maximum number of iterations exceeded."
2989            retval = 3
2990            break
2991        if (maxaccept is not None) and (schedule.accepted > maxaccept):
2992            retval = 4
2993            break
2994
2995    if full_output:
2996        return best_state.x, best_state.cost, schedule.T, \
2997               schedule.feval, iters, schedule.accepted, retval
2998    else:
2999        return best_state.x, retval
3000
3001def worker(iCyc,data,RBdata,reflType,reflData,covData,out_q,nprocess=-1):
3002    outlist = []
3003    random.seed(int(time.time())%100000+nprocess)   #make sure each process has a different random start
3004    for n in range(iCyc):
3005        result = mcsaSearch(data,RBdata,reflType,reflData,covData,None)
3006        outlist.append(result[0])
3007        print ' MC/SA residual: %.3f%% structure factor time: %.3f'%(100*result[0][2],result[1])
3008    out_q.put(outlist)
3009
3010def MPmcsaSearch(nCyc,data,RBdata,reflType,reflData,covData):
3011    import multiprocessing as mp
3012   
3013    nprocs = mp.cpu_count()
3014    out_q = mp.Queue()
3015    procs = []
3016    iCyc = np.zeros(nprocs)
3017    for i in range(nCyc):
3018        iCyc[i%nprocs] += 1
3019    for i in range(nprocs):
3020        p = mp.Process(target=worker,args=(int(iCyc[i]),data,RBdata,reflType,reflData,covData,out_q,i))
3021        procs.append(p)
3022        p.start()
3023    resultlist = []
3024    for i in range(nprocs):
3025        resultlist += out_q.get()
3026    for p in procs:
3027        p.join()
3028    return resultlist
3029
3030def mcsaSearch(data,RBdata,reflType,reflData,covData,pgbar):
3031    '''default doc string
3032   
3033    :param type name: description
3034   
3035    :returns: type name: description
3036   
3037    '''
3038   
3039    global tsum
3040    tsum = 0.
3041   
3042    def getMDparms(item,pfx,parmDict,varyList):
3043        parmDict[pfx+'MDaxis'] = item['axis']
3044        parmDict[pfx+'MDval'] = item['Coef'][0]
3045        if item['Coef'][1]:
3046            varyList += [pfx+'MDval',]
3047            limits = item['Coef'][2]
3048            lower.append(limits[0])
3049            upper.append(limits[1])
3050                       
3051    def getAtomparms(item,pfx,aTypes,SGData,parmDict,varyList):
3052        parmDict[pfx+'Atype'] = item['atType']
3053        aTypes |= set([item['atType'],]) 
3054        pstr = ['Ax','Ay','Az']
3055        XYZ = [0,0,0]
3056        for i in range(3):
3057            name = pfx+pstr[i]
3058            parmDict[name] = item['Pos'][0][i]
3059            XYZ[i] = parmDict[name]
3060            if item['Pos'][1][i]:
3061                varyList += [name,]
3062                limits = item['Pos'][2][i]
3063                lower.append(limits[0])
3064                upper.append(limits[1])
3065        parmDict[pfx+'Amul'] = len(G2spc.GenAtom(XYZ,SGData))
3066           
3067    def getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList):
3068        parmDict[mfx+'MolCent'] = item['MolCent']
3069        parmDict[mfx+'RBId'] = item['RBId']
3070        pstr = ['Px','Py','Pz']
3071        ostr = ['Qa','Qi','Qj','Qk']    #angle,vector not quaternion
3072        for i in range(3):
3073            name = mfx+pstr[i]
3074            parmDict[name] = item['Pos'][0][i]
3075            if item['Pos'][1][i]:
3076                varyList += [name,]
3077                limits = item['Pos'][2][i]
3078                lower.append(limits[0])
3079                upper.append(limits[1])
3080        AV = item['Ori'][0]
3081        A = AV[0]
3082        V = AV[1:]
3083        for i in range(4):
3084            name = mfx+ostr[i]
3085            if i:
3086                parmDict[name] = V[i-1]
3087            else:
3088                parmDict[name] = A
3089            if item['Ovar'] == 'AV':
3090                varyList += [name,]
3091                limits = item['Ori'][2][i]
3092                lower.append(limits[0])
3093                upper.append(limits[1])
3094            elif item['Ovar'] == 'A' and not i:
3095                varyList += [name,]
3096                limits = item['Ori'][2][i]
3097                lower.append(limits[0])
3098                upper.append(limits[1])
3099        if 'Tor' in item:      #'Tor' not there for 'Vector' RBs
3100            for i in range(len(item['Tor'][0])):
3101                name = mfx+'Tor'+str(i)
3102                parmDict[name] = item['Tor'][0][i]
3103                if item['Tor'][1][i]:
3104                    varyList += [name,]
3105                    limits = item['Tor'][2][i]
3106                    lower.append(limits[0])
3107                    upper.append(limits[1])
3108        atypes = RBdata[item['Type']][item['RBId']]['rbTypes']
3109        aTypes |= set(atypes)
3110        atNo += len(atypes)
3111        return atNo
3112               
3113    def GetAtomM(Xdata,SGData):
3114        Mdata = []
3115        for xyz in Xdata:
3116            Mdata.append(float(len(G2spc.GenAtom(xyz,SGData))))
3117        return np.array(Mdata)
3118       
3119    def GetAtomT(RBdata,parmDict):
3120        'Needs a doc string'
3121        atNo = parmDict['atNo']
3122        nfixAt = parmDict['nfixAt']
3123        Tdata = atNo*[' ',]
3124        for iatm in range(nfixAt):
3125            parm = ':'+str(iatm)+':Atype'
3126            if parm in parmDict:
3127                Tdata[iatm] = aTypes.index(parmDict[parm])
3128        iatm = nfixAt
3129        for iObj in range(parmDict['nObj']):
3130            pfx = str(iObj)+':'
3131            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3132                if parmDict[pfx+'Type'] == 'Vector':
3133                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3134                    nAtm = len(RBRes['rbVect'][0])
3135                else:       #Residue
3136                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3137                    nAtm = len(RBRes['rbXYZ'])
3138                for i in range(nAtm):
3139                    Tdata[iatm] = aTypes.index(RBRes['rbTypes'][i])
3140                    iatm += 1
3141            elif parmDict[pfx+'Type'] == 'Atom':
3142                atNo = parmDict[pfx+'atNo']
3143                parm = pfx+'Atype'              #remove extra ':'
3144                if parm in parmDict:
3145                    Tdata[atNo] = aTypes.index(parmDict[parm])
3146                iatm += 1
3147            else:
3148                continue        #skips March Dollase
3149        return Tdata
3150       
3151    def GetAtomX(RBdata,parmDict):
3152        'Needs a doc string'
3153        Bmat = parmDict['Bmat']
3154        atNo = parmDict['atNo']
3155        nfixAt = parmDict['nfixAt']
3156        Xdata = np.zeros((3,atNo))
3157        keys = {':Ax':Xdata[0],':Ay':Xdata[1],':Az':Xdata[2]}
3158        for iatm in range(nfixAt):
3159            for key in keys:
3160                parm = ':'+str(iatm)+key
3161                if parm in parmDict:
3162                    keys[key][iatm] = parmDict[parm]
3163        iatm = nfixAt
3164        for iObj in range(parmDict['nObj']):
3165            pfx = str(iObj)+':'
3166            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3167                if parmDict[pfx+'Type'] == 'Vector':
3168                    RBRes = RBdata['Vector'][parmDict[pfx+'RBId']]
3169                    vecs = RBRes['rbVect']
3170                    mags = RBRes['VectMag']
3171                    Cart = np.zeros_like(vecs[0])
3172                    for vec,mag in zip(vecs,mags):
3173                        Cart += vec*mag
3174                elif parmDict[pfx+'Type'] == 'Residue':
3175                    RBRes = RBdata['Residue'][parmDict[pfx+'RBId']]
3176                    Cart = np.array(RBRes['rbXYZ'])
3177                    for itor,seq in enumerate(RBRes['rbSeq']):
3178                        QuatA = AVdeg2Q(parmDict[pfx+'Tor'+str(itor)],Cart[seq[0]]-Cart[seq[1]])
3179                        Cart[seq[3]] = prodQVQ(QuatA,Cart[seq[3]]-Cart[seq[1]])+Cart[seq[1]]
3180                if parmDict[pfx+'MolCent'][1]:
3181                    Cart -= parmDict[pfx+'MolCent'][0]
3182                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3183                Pos = np.array([parmDict[pfx+'Px'],parmDict[pfx+'Py'],parmDict[pfx+'Pz']])
3184                Xdata.T[iatm:iatm+len(Cart)] = np.inner(Bmat,prodQVQ(Qori,Cart)).T+Pos
3185                iatm += len(Cart)
3186            elif parmDict[pfx+'Type'] == 'Atom':
3187                atNo = parmDict[pfx+'atNo']
3188                for key in keys:
3189                    parm = pfx+key[1:]              #remove extra ':'
3190                    if parm in parmDict:
3191                        keys[key][atNo] = parmDict[parm]
3192                iatm += 1
3193            else:
3194                continue        #skips March Dollase
3195        return Xdata.T
3196       
3197    def getAllTX(Tdata,Mdata,Xdata,SGM,SGT):
3198        allX = np.inner(Xdata,SGM)+SGT
3199        allT = np.repeat(Tdata,allX.shape[1])
3200        allM = np.repeat(Mdata,allX.shape[1])
3201        allX = np.reshape(allX,(-1,3))
3202        return allT,allM,allX
3203
3204    def getAllX(Xdata,SGM,SGT):
3205        allX = np.inner(Xdata,SGM)+SGT
3206        allX = np.reshape(allX,(-1,3))
3207        return allX
3208       
3209    def normQuaternions(RBdata,parmDict,varyList,values):
3210        for iObj in range(parmDict['nObj']):
3211            pfx = str(iObj)+':'
3212            if parmDict[pfx+'Type'] in ['Vector','Residue']:
3213                Qori = AVdeg2Q(parmDict[pfx+'Qa'],[parmDict[pfx+'Qi'],parmDict[pfx+'Qj'],parmDict[pfx+'Qk']])
3214                A,V = Q2AVdeg(Qori)
3215                for i,name in enumerate(['Qa','Qi','Qj','Qk']):
3216                    if i:
3217                        parmDict[pfx+name] = V[i-1]
3218                    else:
3219                        parmDict[pfx+name] = A
3220       
3221    def mcsaCalc(values,refList,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict):
3222        ''' Compute structure factors for all h,k,l for phase
3223        input:
3224            refList: [ref] where each ref = h,k,l,m,d,...
3225            rcov:   array[nref,nref] covariance terms between Fo^2 values
3226            ifInv:  bool True if centrosymmetric
3227            allFF: array[nref,natoms] each value is mult*FF(H)/max(mult)
3228            RBdata: [dict] rigid body dictionary
3229            varyList: [list] names of varied parameters in MC/SA (not used here)           
3230            ParmDict: [dict] problem parameters
3231        puts result F^2 in each ref[5] in refList
3232        returns:
3233            delt-F*rcov*delt-F/sum(Fo^2)^2
3234        '''       
3235        global tsum
3236        t0 = time.time()
3237        parmDict.update(dict(zip(varyList,values)))             #update parameter tables
3238        Xdata = GetAtomX(RBdata,parmDict)                       #get new atom coords from RB
3239        allX = getAllX(Xdata,SGM,SGT)                           #fill unit cell - dups. OK
3240        MDval = parmDict['0:MDval']                             #get March-Dollase coeff
3241        HX2pi = 2.*np.pi*np.inner(allX,refList[:3].T)           #form 2piHX for every H,X pair
3242        Aterm = refList[3]*np.sum(allFF*np.cos(HX2pi),axis=0)**2    #compute real part for all H
3243        refList[5] = Aterm
3244        if not ifInv:
3245            refList[5] += refList[3]*np.sum(allFF*np.sin(HX2pi),axis=0)**2    #imaginary part for all H
3246        if len(cosTable):        #apply MD correction
3247            refList[5] *= np.sum(np.sqrt((MDval/(cosTable*(MDval**3-1.)+1.))**3),axis=1)/cosTable.shape[1]
3248        sumFcsq = np.sum(refList[5])
3249        scale = parmDict['sumFosq']/sumFcsq
3250        refList[5] *= scale
3251        refList[6] = refList[4]-refList[5]
3252        M = np.inner(refList[6],np.inner(rcov,refList[6]))
3253        tsum += (time.time()-t0)
3254        return M/np.sum(refList[4]**2)
3255
3256    sq8ln2 = np.sqrt(8*np.log(2))
3257    sq2pi = np.sqrt(2*np.pi)
3258    sq4pi = np.sqrt(4*np.pi)
3259    generalData = data['General']
3260    Amat,Bmat = G2lat.cell2AB(generalData['Cell'][1:7])
3261    Gmat,gmat = G2lat.cell2Gmat(generalData['Cell'][1:7])
3262    SGData = generalData['SGData']
3263    SGM = np.array([SGData['SGOps'][i][0] for i in range(len(SGData['SGOps']))])
3264    SGMT = np.array([SGData['SGOps'][i][0].T for i in range(len(SGData['SGOps']))])
3265    SGT = np.array([SGData['SGOps'][i][1] for i in range(len(SGData['SGOps']))])
3266    fixAtoms = data['Atoms']                       #if any
3267    cx,ct,cs = generalData['AtomPtrs'][:3]
3268    aTypes = set([])
3269    parmDict = {'Bmat':Bmat,'Gmat':Gmat}
3270    varyList = []
3271    atNo = 0
3272    for atm in fixAtoms:
3273        pfx = ':'+str(atNo)+':'
3274        parmDict[pfx+'Atype'] = atm[ct]
3275        aTypes |= set([atm[ct],])
3276        pstr = ['Ax','Ay','Az']
3277        parmDict[pfx+'Amul'] = atm[cs+1]
3278        for i in range(3):
3279            name = pfx+pstr[i]
3280            parmDict[name] = atm[cx+i]
3281        atNo += 1
3282    parmDict['nfixAt'] = len(fixAtoms)       
3283    MCSA = generalData['MCSA controls']
3284    reflName = MCSA['Data source']
3285    phaseName = generalData['Name']
3286    MCSAObjs = data['MCSA']['Models']               #list of MCSA models
3287    upper = []
3288    lower = []
3289    MDvec = np.zeros(3)
3290    for i,item in enumerate(MCSAObjs):
3291        mfx = str(i)+':'
3292        parmDict[mfx+'Type'] = item['Type']
3293        if item['Type'] == 'MD':
3294            getMDparms(item,mfx,parmDict,varyList)
3295            MDvec = np.array(item['axis'])
3296        elif item['Type'] == 'Atom':
3297            getAtomparms(item,mfx,aTypes,SGData,parmDict,varyList)
3298            parmDict[mfx+'atNo'] = atNo
3299            atNo += 1
3300        elif item['Type'] in ['Residue','Vector']:
3301            atNo = getRBparms(item,mfx,aTypes,RBdata,SGData,atNo,parmDict,varyList)
3302    parmDict['atNo'] = atNo                 #total no. of atoms
3303    parmDict['nObj'] = len(MCSAObjs)
3304    aTypes = list(aTypes)
3305    Tdata = GetAtomT(RBdata,parmDict)
3306    Xdata = GetAtomX(RBdata,parmDict)
3307    Mdata = GetAtomM(Xdata,SGData)
3308    allT,allM = getAllTX(Tdata,Mdata,Xdata,SGM,SGT)[:2]
3309    FFtables = G2el.GetFFtable(aTypes)
3310    refs = []
3311    allFF = []
3312    cosTable = []
3313    sumFosq = 0
3314    if 'PWDR' in reflName:
3315        for ref in reflData:
3316            h,k,l,m,d,pos,sig,gam,f = ref[:9]
3317            if d >= MCSA['dmin']:
3318                sig = G2pwd.getgamFW(sig,gam)/sq8ln2        #--> sig from FWHM
3319                SQ = 0.25/d**2
3320                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3321                refs.append([h,k,l,m,f*m,pos,sig])
3322                sumFosq += f*m
3323                Heqv = np.inner(np.array([h,k,l]),SGMT)
3324                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3325        nRef = len(refs)
3326        cosTable = np.array(cosTable)**2
3327        rcov = np.zeros((nRef,nRef))
3328        for iref,refI in enumerate(refs):
3329            rcov[iref][iref] = 1./(sq4pi*refI[6])
3330            for jref,refJ in enumerate(refs[:iref]):
3331                t1 = refI[6]**2+refJ[6]**2
3332                t2 = (refJ[5]-refI[5])**2/(2.*t1)
3333                if t2 > 10.:
3334                    rcov[iref][jref] = 0.
3335                else:
3336                    rcov[iref][jref] = 1./(sq2pi*np.sqrt(t1)*np.exp(t2))
3337        rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3338        Rdiag = np.sqrt(np.diag(rcov))
3339        Rnorm = np.outer(Rdiag,Rdiag)
3340        rcov /= Rnorm
3341    elif 'Pawley' in reflName:  #need a bail out if Pawley cov matrix doesn't exist.
3342        vNames = []
3343        pfx = str(data['pId'])+'::PWLref:'
3344        for iref,refI in enumerate(reflData):           #Pawley reflection set
3345            h,k,l,m,d,v,f,s = refI
3346            if d >= MCSA['dmin'] and v:       #skip unrefined ones
3347                vNames.append(pfx+str(iref))
3348                SQ = 0.25/d**2
3349                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3350                refs.append([h,k,l,m,f*m,iref,0.])
3351                sumFosq += f*m
3352                Heqv = np.inner(np.array([h,k,l]),SGMT)
3353                cosTable.append(G2lat.CosAngle(Heqv,MDvec,Gmat))
3354        cosTable = np.array(cosTable)**2
3355        nRef = len(refs)
3356#        if generalData['doPawley'] and (covData['freshCOV'] or  MCSA['newDmin']):
3357        if covData['freshCOV'] or  MCSA['newDmin']:
3358            vList = covData['varyList']
3359            covMatrix = covData['covMatrix']
3360            rcov = getVCov(vNames,vList,covMatrix)
3361            rcov += (rcov.T-np.diagflat(np.diagonal(rcov)))
3362            Rdiag = np.sqrt(np.diag(rcov))
3363            Rdiag = np.where(Rdiag,Rdiag,1.0)
3364            Rnorm = np.outer(Rdiag,Rdiag)
3365            rcov /= Rnorm
3366            MCSA['rcov'] = rcov
3367            covData['freshCOV'] = False
3368            MCSA['newDmin'] = False
3369        else:
3370            rcov = MCSA['rcov']
3371    elif 'HKLF' in reflName:
3372        for ref in reflData:
3373            [h,k,l,m,d],f = ref[:5],ref[6]
3374            if d >= MCSA['dmin']:
3375                SQ = 0.25/d**2
3376                allFF.append(allM*[G2el.getFFvalues(FFtables,SQ,True)[i] for i in allT]/np.max(allM))
3377                refs.append([h,k,l,m,f*m,0.,0.])
3378                sumFosq += f*m
3379        nRef = len(refs)
3380        rcov = np.identity(len(refs))
3381    allFF = np.array(allFF).T
3382    refs = np.array(refs).T
3383    print ' Minimum d-spacing used: %.2f No. reflections used: %d'%(MCSA['dmin'],nRef)
3384    print ' Number of parameters varied: %d'%(len(varyList))
3385    parmDict['sumFosq'] = sumFosq
3386    x0 = [parmDict[val] for val in varyList]
3387    ifInv = SGData['SGInv']
3388    # consider replacing anneal with scipy.optimize.basinhopping
3389    results = anneal(mcsaCalc,x0,args=(refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict),
3390        schedule=MCSA['Algorithm'], full_output=True,
3391        T0=MCSA['Annealing'][0], Tf=MCSA['Annealing'][1],dwell=MCSA['Annealing'][2],
3392        boltzmann=MCSA['boltzmann'], learn_rate=0.5, 
3393        quench=MCSA['fast parms'][0], m=MCSA['fast parms'][1], n=MCSA['fast parms'][2],
3394        lower=lower, upper=upper, slope=MCSA['log slope'],ranStart=MCSA.get('ranStart',False),
3395        ranRange=MCSA.get('ranRange',0.10),autoRan=MCSA.get('autoRan',False),dlg=pgbar)
3396    M = mcsaCalc(results[0],refs,rcov,cosTable,ifInv,allFF,RBdata,varyList,parmDict)
3397#    for ref in refs.T:
3398#        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])
3399#    print np.sqrt((np.sum(refs[6]**2)/np.sum(refs[4]**2)))
3400    Result = [False,False,results[1],results[2],]+list(results[0])
3401    Result.append(varyList)
3402    return Result,tsum
3403
3404       
3405################################################################################
3406##### Quaternion stuff
3407################################################################################
3408
3409def prodQQ(QA,QB):
3410    ''' Grassman quaternion product
3411        QA,QB quaternions; q=r+ai+bj+ck
3412    '''
3413    D = np.zeros(4)
3414    D[0] = QA[0]*QB[0]-QA[1]*QB[1]-QA[2]*QB[2]-QA[3]*QB[3]
3415    D[1] = QA[0]*QB[1]+QA[1]*QB[0]+QA[2]*QB[3]-QA[3]*QB[2]
3416    D[2] = QA[0]*QB[2]-QA[1]*QB[3]+QA[2]*QB[0]+QA[3]*QB[1]
3417    D[3] = QA[0]*QB[3]+QA[1]*QB[2]-QA[2]*QB[1]+QA[3]*QB[0]
3418   
3419#    D[0] = QA[0]*QB[0]-np.dot(QA[1:],QB[1:])
3420#    D[1:] = QA[0]*QB[1:]+QB[0]*QA[1:]+np.cross(QA[1:],QB[1:])
3421   
3422    return D
3423   
3424def normQ(QA):
3425    ''' get length of quaternion & normalize it
3426        q=r+ai+bj+ck
3427    '''
3428    n = np.sqrt(np.sum(np.array(QA)**2))
3429    return QA/n
3430   
3431def invQ(Q):
3432    '''
3433        get inverse of quaternion
3434        q=r+ai+bj+ck; q* = r-ai-bj-ck
3435    '''
3436    return Q*np.array([1,-1,-1,-1])
3437   
3438def prodQVQ(Q,V):
3439    """
3440    compute the quaternion vector rotation qvq-1 = v'
3441    q=r+ai+bj+ck
3442    """
3443    T2 = Q[0]*Q[1]
3444    T3 = Q[0]*Q[2]
3445    T4 = Q[0]*Q[3]
3446    T5 = -Q[1]*Q[1]
3447    T6 = Q[1]*Q[2]
3448    T7 = Q[1]*Q[3]
3449    T8 = -Q[2]*Q[2]
3450    T9 = Q[2]*Q[3]
3451    T10 = -Q[3]*Q[3]
3452    M = np.array([[T8+T10,T6-T4,T3+T7],[T4+T6,T5+T10,T9-T2],[T7-T3,T2+T9,T5+T8]])
3453    VP = 2.*np.inner(V,M)
3454    return VP+V
3455   
3456def Q2Mat(Q):
3457    ''' make rotation matrix from quaternion
3458        q=r+ai+bj+ck
3459    '''
3460    QN = normQ(Q)
3461    aa = QN[0]**2
3462    ab = QN[0]*QN[1]
3463    ac = QN[0]*QN[2]
3464    ad = QN[0]*QN[3]
3465    bb = QN[1]**2
3466    bc = QN[1]*QN[2]
3467    bd = QN[1]*QN[3]
3468    cc = QN[2]**2
3469    cd = QN[2]*QN[3]
3470    dd = QN[3]**2
3471    M = [[aa+bb-cc-dd, 2.*(bc-ad),  2.*(ac+bd)],
3472        [2*(ad+bc),   aa-bb+cc-dd,  2.*(cd-ab)],
3473        [2*(bd-ac),    2.*(ab+cd), aa-bb-cc+dd]]
3474    return np.array(M)
3475   
3476def AV2Q(A,V):
3477    ''' convert angle (radians) & vector to quaternion
3478        q=r+ai+bj+ck
3479    '''
3480    Q = np.zeros(4)
3481    d = nl.norm(np.array(V))
3482    if d:
3483        V /= d
3484        if not A:       #==0.
3485            A = 2.*np.pi
3486        p = A/2.
3487        Q[0] = np.cos(p)
3488        Q[1:4] = V*np.sin(p)
3489    else:
3490        Q[3] = 1.
3491    return Q
3492   
3493def AVdeg2Q(A,V):
3494    ''' convert angle (degrees) & vector to quaternion
3495        q=r+ai+bj+ck
3496    '''
3497    Q = np.zeros(4)
3498    d = nl.norm(np.array(V))
3499    if not A:       #== 0.!
3500        A = 360.
3501    if d:
3502        V /= d
3503        p = A/2.
3504        Q[0] = cosd(p)
3505        Q[1:4] = V*sind(p)
3506    else:
3507        Q[3] = 1.
3508    return Q
3509   
3510def Q2AVdeg(Q):
3511    ''' convert quaternion to angle (degrees 0-360) & normalized vector
3512        q=r+ai+bj+ck
3513    '''
3514    A = 2.*acosd(Q[0])
3515    V = np.array(Q[1:])
3516    V /= sind(A/2.)
3517    return A,V
3518   
3519def Q2AV(Q):
3520    ''' convert quaternion to angle (radians 0-2pi) & normalized vector
3521        q=r+ai+bj+ck
3522    '''
3523    A = 2.*np.arccos(Q[0])
3524    V = np.array(Q[1:])
3525    V /= np.sin(A/2.)
3526    return A,V
3527   
3528def randomQ(r0,r1,r2,r3):
3529    ''' create random quaternion from 4 random numbers in range (-1,1)
3530    '''
3531    sum = 0
3532    Q = np.array(4)
3533    Q[0] = r0
3534    sum += Q[0]**2
3535    Q[1] = np.sqrt(1.-sum)*r1
3536    sum += Q[1]**2
3537    Q[2] = np.sqrt(1.-sum)*r2
3538    sum += Q[2]**2
3539    Q[3] = np.sqrt(1.-sum)*np.where(r3<0.,-1.,1.)
3540    return Q
3541   
3542def randomAVdeg(r0,r1,r2,r3):
3543    ''' create random angle (deg),vector from 4 random number in range (-1,1)
3544    '''
3545    return Q2AVdeg(randomQ(r0,r1,r2,r3))
3546   
3547def makeQuat(A,B,C):
3548    ''' Make quaternion from rotation of A vector to B vector about C axis
3549
3550        :param np.array A,B,C: Cartesian 3-vectors
3551        :returns: quaternion & rotation angle in radians q=r+ai+bj+ck
3552    '''
3553
3554    V1 = np.cross(A,C)
3555    V2 = np.cross(B,C)
3556    if nl.norm(V1)*nl.norm(V2):
3557        V1 /= nl.norm(V1)
3558        V2 /= nl.norm(V2)
3559        V3 = np.cross(V1,V2)
3560    else:
3561        V3 = np.zeros(3)
3562    Q = np.array([0.,0.,0.,1.])
3563    D = 0.
3564    if nl.norm(V3):
3565        V3 /= nl.norm(V3)
3566        D1 = min(1.0,max(-1.0,np.vdot(V1,V2)))
3567        D = np.arccos(D1)/2.0
3568        V1 = C-V3
3569        V2 = C+V3
3570        DM = nl.norm(V1)
3571        DP = nl.norm(V2)
3572        S = np.sin(D)
3573        Q[0] = np.cos(D)
3574        Q[1:] = V3*S
3575        D *= 2.
3576        if DM > DP:
3577            D *= -1.
3578    return Q,D
3579   
3580def annealtests():
3581    from numpy import cos
3582    # minimum expected at ~-0.195
3583    func = lambda x: cos(14.5*x-0.3) + (x+0.2)*x
3584    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='cauchy')
3585    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='fast')
3586    print anneal(func,1.0,full_output=1,upper=3.0,lower=-3.0,feps=1e-4,maxiter=2000,schedule='boltzmann')
3587
3588    # minimum expected at ~[-0.195, -0.1]
3589    func = lambda x: cos(14.5*x[0]-0.3) + (x[1]+0.2)*x[1] + (x[0]+0.2)*x[0]
3590    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')
3591    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')
3592    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')
3593
3594
3595if __name__ == '__main__':
3596    annealtests()
Note: See TracBrowser for help on using the repository browser.