source: trunk/GSASIImath.py @ 485

Last change on this file since 485 was 485, checked in by vondreele, 11 years ago
File size: 17.1 KB
Line 
1#GSASIImath - major mathematics routines
2########### SVN repository information ###################
3# $Date: 2012-01-13 11:48:53 -0600 (Fri, 13 Jan 2012) $
4# $Author: vondreele & toby $
5# $Revision: 451 $
6# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIImath.py $
7# $Id: GSASIImath.py 451 2012-01-13 17:48:53Z vondreele $
8########### SVN repository information ###################
9import sys
10import os
11import os.path as ospath
12import numpy as np
13import numpy.linalg as nl
14import cPickle
15import time
16import math
17import GSASIIpath
18import GSASIIspc as G2spc
19import scipy.optimize as so
20import scipy.linalg as sl
21
22sind = lambda x: np.sin(x*np.pi/180.)
23cosd = lambda x: np.cos(x*np.pi/180.)
24tand = lambda x: np.tan(x*np.pi/180.)
25asind = lambda x: 180.*np.arcsin(x)/np.pi
26acosd = lambda x: 180.*np.arccos(x)/np.pi
27atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
28
29def HessianLSQ(func,x0,Hess,args=(),ftol=1.49012e-8,xtol=1.49012e-8, maxcyc=0):
30   
31    """
32    Minimize the sum of squares of a set of equations.
33
34    ::
35   
36                    Nobs
37        x = arg min(sum(func(y)**2,axis=0))
38                    y=0
39
40    Parameters
41    ----------
42    func : callable
43        should take at least one (possibly length N vector) argument and
44        returns M floating point numbers.
45    x0 : ndarray
46        The starting estimate for the minimization of length N
47    Hess : callable
48        A required function or method to compute the weighted vector and Hessian for func.
49        It must be a symmetric NxN array
50    args : tuple
51        Any extra arguments to func are placed in this tuple.
52    ftol : float
53        Relative error desired in the sum of squares.
54    xtol : float
55        Relative error desired in the approximate solution.
56    maxcyc : int
57        The maximum number of cycles of refinement to execute, if -1 refine
58        until other limits are met (ftol, xtol)
59
60    Returns
61    -------
62    x : ndarray
63        The solution (or the result of the last iteration for an unsuccessful
64        call).
65    cov_x : ndarray
66        Uses the fjac and ipvt optional outputs to construct an
67        estimate of the jacobian around the solution.  ``None`` if a
68        singular matrix encountered (indicates very flat curvature in
69        some direction).  This matrix must be multiplied by the
70        residual standard deviation to get the covariance of the
71        parameter estimates -- see curve_fit.
72    infodict : dict
73        a dictionary of optional outputs with the key s::
74
75            - 'fvec' : the function evaluated at the output
76
77
78    Notes
79    -----
80
81    """
82               
83    x0 = np.array(x0, ndmin=1)      #might be redundant?
84    n = len(x0)
85    if type(args) != type(()):
86        args = (args,)
87       
88    icycle = 0
89    One = np.ones((n,n))
90    lam = 0.001
91    lamMax = lam
92    nfev = 0
93    while icycle <= maxcyc:
94        lamMax = max(lamMax,lam)
95        M = func(x0,*args)
96        nfev += 1
97        chisq0 = np.sum(M**2)
98        Yvec,Amat = Hess(x0,*args)
99        Adiag = np.sqrt(np.diag(Amat))
100        if 0.0 in Adiag:                #hard singularity in matrix
101            psing = list(np.where(Adiag == 0.)[0])
102            return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
103        Anorm = np.outer(Adiag,Adiag)
104        Yvec /= Adiag
105        Amat /= Anorm       
106        while True:
107            Lam = np.eye(Amat.shape[0])*lam
108            Amatlam = Amat*(One+Lam)
109            try:
110                Xvec = nl.solve(Amatlam,Yvec)
111            except LinAlgError:
112                psing = list(np.where(np.diag(nl.gr(Amatlam)[1]) < 1.e-14)[0])
113                return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
114            Xvec /= Adiag
115            M2 = func(x0+Xvec,*args)
116            nfev += 1
117            chisq1 = np.sum(M2**2)
118            if chisq1 > chisq0:
119                lam *= 10.
120            else:
121                x0 += Xvec
122                lam /= 10.
123                break
124        if (chisq0-chisq1)/chisq0 < ftol:
125            break
126        icycle += 1
127    M = func(x0,*args)
128    nfev += 1
129    Yvec,Amat = Hess(x0,*args)
130    try:
131        Bmat = nl.inv(Amat)
132        return [x0,Bmat,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':[]}]
133    except LinAlgError:
134        psing = list(np.where(np.diag(nl.gr(Amat)[1]) < 1.e-14)[0])
135        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}] 
136   
137def getVCov(varyNames,varyList,covMatrix):
138    vcov = np.zeros((len(varyNames),len(varyNames)))
139    for i1,name1 in enumerate(varyNames):
140        for i2,name2 in enumerate(varyNames):
141            try:
142                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
143            except ValueError:
144                vcov[i1][i2] = 0.0
145    return vcov
146   
147def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
148   
149    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
150        TxT = inv*(np.inner(M,Tx)+T)+C+U
151        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
152       
153    inv = Top/abs(Top)
154    cent = abs(Top)/100
155    op = abs(Top)%100-1
156    M,T = SGData['SGOps'][op]
157    C = SGData['SGCen'][cent]
158    dx = .00001
159    deriv = np.zeros(6)
160    for i in [0,1,2]:
161        Oxyz[i] += dx
162        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
163        Oxyz[i] -= 2*dx
164        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
165        Oxyz[i] += dx
166        Txyz[i] += dx
167        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
168        Txyz[i] -= 2*dx
169        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
170        Txyz[i] += dx
171    return deriv
172   
173def getAngSig(VA,VB,Amat,SGData,covData={}):
174   
175    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
176        TxT = inv*(np.inner(M,Tx)+T)+C
177        TxT = G2spc.MoveToUnitCell(TxT)+U
178        return np.inner(Amat,(TxT-Ox))
179       
180    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
181        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
182        VecA /= np.sqrt(np.sum(VecA**2))
183        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
184        VecB /= np.sqrt(np.sum(VecB**2))
185        edge = VecB-VecA
186        edge = np.sum(edge**2)
187        angle = (2.-edge)/2.
188        angle = max(angle,-1.)
189        return acosd(angle)
190       
191    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
192    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
193    invA = invB = 1
194    invA = TopA/abs(TopA)
195    invB = TopB/abs(TopB)
196    centA = abs(TopA)/100
197    centB = abs(TopB)/100
198    opA = abs(TopA)%100-1
199    opB = abs(TopB)%100-1
200    MA,TA = SGData['SGOps'][opA]
201    MB,TB = SGData['SGOps'][opB]
202    CA = SGData['SGCen'][centA]
203    CB = SGData['SGCen'][centB]
204    if 'covMatrix' in covData:
205        covMatrix = covData['covMatrix']
206        varyList = covData['varyList']
207        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
208        dx = .00001
209        dadx = np.zeros(9)
210        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
211        for i in [0,1,2]:
212            OxA[i] += dx
213            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
214            OxA[i] -= 2*dx
215            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
216            OxA[i] += dx
217           
218            TxA[i] += dx
219            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
220            TxA[i] -= 2*dx
221            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
222            TxA[i] += dx
223           
224            TxB[i] += dx
225            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
226            TxB[i] -= 2*dx
227            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
228            TxB[i] += dx
229           
230        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
231        if sigAng < 0.01:
232            sigAng = 0.0
233        return Ang,sigAng
234    else:
235        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
236
237def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
238
239    def calcDist(Atoms,SyOps,Amat):
240        XYZ = []
241        for i,atom in enumerate(Atoms):
242            Inv,M,T,C,U = SyOps[i]
243            XYZ.append(np.array(atom[1:4]))
244            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
245            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
246        V1 = XYZ[1]-XYZ[0]
247        return np.sqrt(np.sum(V1**2))
248       
249    Inv = []
250    SyOps = []
251    names = []
252    for i,atom in enumerate(Oatoms):
253        names += atom[-1]
254        Op,unit = Atoms[i][-1]
255        inv = Op/abs(Op)
256        m,t = SGData['SGOps'][abs(Op)%100-1]
257        c = SGData['SGCen'][abs(Op)/100]
258        SyOps.append([inv,m,t,c,unit])
259    Dist = calcDist(Oatoms,SyOps,Amat)
260   
261    sig = -0.001
262    if 'covMatrix' in covData:
263        parmNames = []
264        dx = .00001
265        dadx = np.zeros(6)
266        for i in range(6):
267            ia = i/3
268            ix = i%3
269            Oatoms[ia][ix+1] += dx
270            a0 = calcDist(Oatoms,SyOps,Amat)
271            Oatoms[ia][ix+1] -= 2*dx
272            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
273        covMatrix = covData['covMatrix']
274        varyList = covData['varyList']
275        DistVcov = getVCov(names,varyList,covMatrix)
276        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
277        if sig < 0.001:
278            sig = -0.001
279   
280    return Dist,sig
281
282def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
283       
284    def calcAngle(Atoms,SyOps,Amat):
285        XYZ = []
286        for i,atom in enumerate(Atoms):
287            Inv,M,T,C,U = SyOps[i]
288            XYZ.append(np.array(atom[1:4]))
289            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
290            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
291        V1 = XYZ[1]-XYZ[0]
292        V1 /= np.sqrt(np.sum(V1**2))
293        V2 = XYZ[1]-XYZ[2]
294        V2 /= np.sqrt(np.sum(V2**2))
295        V3 = V2-V1
296        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
297        return acosd(cang)
298
299    Inv = []
300    SyOps = []
301    names = []
302    for i,atom in enumerate(Oatoms):
303        names += atom[-1]
304        Op,unit = Atoms[i][-1]
305        inv = Op/abs(Op)
306        m,t = SGData['SGOps'][abs(Op)%100-1]
307        c = SGData['SGCen'][abs(Op)/100]
308        SyOps.append([inv,m,t,c,unit])
309    Angle = calcAngle(Oatoms,SyOps,Amat)
310   
311    sig = -0.01
312    if 'covMatrix' in covData:
313        parmNames = []
314        dx = .00001
315        dadx = np.zeros(9)
316        for i in range(9):
317            ia = i/3
318            ix = i%3
319            Oatoms[ia][ix+1] += dx
320            a0 = calcAngle(Oatoms,SyOps,Amat)
321            Oatoms[ia][ix+1] -= 2*dx
322            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
323        covMatrix = covData['covMatrix']
324        varyList = covData['varyList']
325        AngVcov = getVCov(names,varyList,covMatrix)
326        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
327        if sig < 0.01:
328            sig = -0.01
329   
330    return Angle,sig
331
332def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
333   
334    def calcTorsion(Atoms,SyOps,Amat):
335       
336        XYZ = []
337        for i,atom in enumerate(Atoms):
338            Inv,M,T,C,U = SyOps[i]
339            XYZ.append(np.array(atom[1:4]))
340            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
341            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
342        V1 = XYZ[1]-XYZ[0]
343        V2 = XYZ[2]-XYZ[1]
344        V3 = XYZ[3]-XYZ[2]
345        V1 /= np.sqrt(np.sum(V1**2))
346        V2 /= np.sqrt(np.sum(V2**2))
347        V3 /= np.sqrt(np.sum(V3**2))
348        M = np.array([V1,V2,V3])
349        D = nl.det(M)
350        Ang = 1.0
351        P12 = np.dot(V1,V2)
352        P13 = np.dot(V1,V3)
353        P23 = np.dot(V2,V3)
354        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
355        return Tors
356           
357    Inv = []
358    SyOps = []
359    names = []
360    for i,atom in enumerate(Oatoms):
361        names += atom[-1]
362        Op,unit = Atoms[i][-1]
363        inv = Op/abs(Op)
364        m,t = SGData['SGOps'][abs(Op)%100-1]
365        c = SGData['SGCen'][abs(Op)/100]
366        SyOps.append([inv,m,t,c,unit])
367    Tors = calcTorsion(Oatoms,SyOps,Amat)
368   
369    sig = -0.01
370    if 'covMatrix' in covData:
371        parmNames = []
372        dx = .00001
373        dadx = np.zeros(12)
374        for i in range(12):
375            ia = i/3
376            ix = i%3
377            Oatoms[ia][ix+1] += dx
378            a0 = calcTorsion(Oatoms,SyOps,Amat)
379            Oatoms[ia][ix+1] -= 2*dx
380            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
381        covMatrix = covData['covMatrix']
382        varyList = covData['varyList']
383        TorVcov = getVCov(names,varyList,covMatrix)
384        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
385        if sig < 0.01:
386            sig = -0.01
387   
388    return Tors,sig
389       
390def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
391   
392    def calcDist(Atoms,SyOps,Amat):
393        XYZ = []
394        for i,atom in enumerate(Atoms):
395            Inv,M,T,C,U = SyOps[i]
396            XYZ.append(np.array(atom[1:4]))
397            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
398            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
399        V1 = XYZ[1]-XYZ[0]
400        return np.sqrt(np.sum(V1**2))
401       
402    def calcAngle(Atoms,SyOps,Amat):
403        XYZ = []
404        for i,atom in enumerate(Atoms):
405            Inv,M,T,C,U = SyOps[i]
406            XYZ.append(np.array(atom[1:4]))
407            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
408            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
409        V1 = XYZ[1]-XYZ[0]
410        V1 /= np.sqrt(np.sum(V1**2))
411        V2 = XYZ[1]-XYZ[2]
412        V2 /= np.sqrt(np.sum(V2**2))
413        V3 = V2-V1
414        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
415        return acosd(cang)
416
417    def calcTorsion(Atoms,SyOps,Amat):
418       
419        XYZ = []
420        for i,atom in enumerate(Atoms):
421            Inv,M,T,C,U = SyOps[i]
422            XYZ.append(np.array(atom[1:4]))
423            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
424            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
425        V1 = XYZ[1]-XYZ[0]
426        V2 = XYZ[2]-XYZ[1]
427        V3 = XYZ[3]-XYZ[2]
428        V1 /= np.sqrt(np.sum(V1**2))
429        V2 /= np.sqrt(np.sum(V2**2))
430        V3 /= np.sqrt(np.sum(V3**2))
431        M = np.array([V1,V2,V3])
432        D = nl.det(M)
433        Ang = 1.0
434        P12 = np.dot(V1,V2)
435        P13 = np.dot(V1,V3)
436        P23 = np.dot(V2,V3)
437        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
438        return Tors
439           
440    Inv = []
441    SyOps = []
442    names = []
443    for i,atom in enumerate(Oatoms):
444        names += atom[-1]
445        Op,unit = Atoms[i][-1]
446        inv = Op/abs(Op)
447        m,t = SGData['SGOps'][abs(Op)%100-1]
448        c = SGData['SGCen'][abs(Op)/100]
449        SyOps.append([inv,m,t,c,unit])
450    M = len(Oatoms)
451    if M == 2:
452        Val = calcDist(Oatoms,SyOps,Amat)
453    elif M == 3:
454        Val = calcAngle(Oatoms,SyOps,Amat)
455    else:
456        Val = calcTorsion(Oatoms,SyOps,Amat)
457   
458    sigVals = [-0.001,-0.01,-0.01]
459    sig = sigVals[M-3]
460    if 'covMatrix' in covData:
461        parmNames = []
462        dx = .00001
463        N = M*3
464        dadx = np.zeros(N)
465        for i in range(N):
466            ia = i/3
467            ix = i%3
468            Oatoms[ia][ix+1] += dx
469            if M == 2:
470                a0 = calcDist(Oatoms,SyOps,Amat)
471            elif M == 3:
472                a0 = calcAngle(Oatoms,SyOps,Amat)
473            else:
474                a0 = calcTorsion(Oatoms,SyOps,Amat)
475            Oatoms[ia][ix+1] -= 2*dx
476            if M == 2:
477                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
478            elif M == 3:
479                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
480            else:
481                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
482        covMatrix = covData['covMatrix']
483        varyList = covData['varyList']
484        Vcov = getVCov(names,varyList,covMatrix)
485        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
486        if sig < sigVals[M-3]:
487            sig = sigVals[M-3]
488   
489    return Val,sig
490       
491   
492def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
493    # returns value(esd) string; nTZ=True for no trailing zeros
494    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
495    #get the 2 significant digits in the esd
496    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
497    #get the number of digits to represent them
498    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
499   
500    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
501    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
502    if esd > 0:
503        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
504        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
505    elif esd < 0:
506         return str(round(value,mdec(esd)-1))
507    else:
508        text = str("%f"%(value))
509        if nTZ:
510            return text.rstrip('0')
511        else:
512            return text
513
514   
Note: See TracBrowser for help on using the repository browser.