source: trunk/GSASIImath.py @ 1929

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

fix sign error in getRho 8-pt interpolation

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