source: trunk/GSASIImath.py @ 1799

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

trap missing data error
numpy array sph.harm. texture fit function & derivative routines - speedier
change pole figure titles

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