source: trunk/GSASIImath.py @ 458

Last change on this file since 458 was 458, checked in by vondreele, 10 years ago

finish dist angle calculations

File size: 9.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 $
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
151        TxT = G2spc.MoveToUnitCell(TxT)+U
152        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
153       
154    inv = 1
155    if Top < 0:
156        inv = -1
157    cent = abs(Top)/100
158    op = abs(Top)%100-1
159    M,T = SGData['SGOps'][op]
160    C = SGData['SGCen'][cent]
161    dx = .00001
162    deriv = np.zeros(6)
163    for i in [0,1,2]:
164        Oxyz[i] += dx
165        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
166        Oxyz[i] -= 2*dx
167        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
168        Oxyz[i] += dx
169        Txyz[i] += dx
170        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
171        Txyz[i] -= 2*dx
172        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
173        Txyz[i] += dx
174    return deriv
175   
176def getAngSig(VA,VB,Amat,SGData,covData={}):
177   
178    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
179        TxT = inv*(np.inner(M,Tx)+T)+C
180        TxT = G2spc.MoveToUnitCell(TxT)+U
181        return np.inner(Amat,(TxT-Ox))
182       
183    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
184        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
185        VecA /= np.sqrt(np.sum(VecA**2))
186        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
187        VecB /= np.sqrt(np.sum(VecB**2))
188        edge = VecB-VecA
189        edge = np.sum(edge**2)
190        angle = (2.-edge)/2.
191        angle = max(angle,-1.)
192        return acosd(angle)
193       
194    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
195    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
196    invA = invB = 1
197    if TopA < 0:
198        invA = -1
199    if TopB < 0:
200        invB = -1
201    centA = abs(TopA)/100
202    centB = abs(TopB)/100
203    opA = abs(TopA)%100-1
204    opB = abs(TopB)%100-1
205    MA,TA = SGData['SGOps'][opA]
206    MB,TB = SGData['SGOps'][opB]
207    CA = SGData['SGCen'][centA]
208    CB = SGData['SGCen'][centB]
209    if 'covMatrix' in covData:
210        covMatrix = covData['covMatrix']
211        varyList = covData['varyList']
212        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
213        dx = .00001
214        dadx = np.zeros(9)
215        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
216        for i in [0,1,2]:
217            OxA[i] += dx
218            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
219            OxA[i] -= 2*dx
220            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
221            OxA[i] += dx
222           
223            TxA[i] += dx
224            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
225            TxA[i] -= 2*dx
226            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
227            TxA[i] += dx
228           
229            TxB[i] += dx
230            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
231            TxB[i] -= 2*dx
232            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
233            TxB[i] += dx
234           
235        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
236        if sigAng < 0.01:
237            sigAng = 0.0
238        return Ang,sigAng
239    else:
240        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
241       
242   
243def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
244    # returns value(esd) string; nTZ=True for no trailing zeros
245    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
246    #get the 2 significant digits in the esd
247    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
248    #get the number of digits to represent them
249    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
250   
251    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
252    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
253    if esd > 0:
254        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
255        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
256    elif esd < 0:
257         return str(round(value,mdec(esd)))
258    else:
259        text = str("%f"%(value))
260        if nTZ:
261            return text.rstrip('0')
262        else:
263            return text
264
265   
Note: See TracBrowser for help on using the repository browser.