source: trunk/GSASIImath.py @ 530

Last change on this file since 530 was 530, checked in by vondreele, 11 years ago

small fix to HessianLSQ about # of cycles
finish map peak search (no GUI/plot output yet)
File size: 17.2 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 nl.LinAlgError:
134        psing = []
135        if maxcyc:
136            psing = list(np.where(np.diag(nl.qr(Amat)[1]) < 1.e-14)[0])
137        return [x0,None,{'num cyc':icycle,'fvec':M,'nfev':nfev,'lamMax':lamMax,'psing':psing}]
138   
139def getVCov(varyNames,varyList,covMatrix):
140    vcov = np.zeros((len(varyNames),len(varyNames)))
141    for i1,name1 in enumerate(varyNames):
142        for i2,name2 in enumerate(varyNames):
143            try:
144                vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)]
145            except ValueError:
146                vcov[i1][i2] = 0.0
147    return vcov
148   
149def getDistDerv(Oxyz,Txyz,Amat,Tunit,Top,SGData):
150   
151    def calcDist(Ox,Tx,U,inv,C,M,T,Amat):
152        TxT = inv*(np.inner(M,Tx)+T)+C+U
153        return np.sqrt(np.sum(np.inner(Amat,(TxT-Ox))**2))
154       
155    inv = Top/abs(Top)
156    cent = abs(Top)/100
157    op = abs(Top)%100-1
158    M,T = SGData['SGOps'][op]
159    C = SGData['SGCen'][cent]
160    dx = .00001
161    deriv = np.zeros(6)
162    for i in [0,1,2]:
163        Oxyz[i] += dx
164        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
165        Oxyz[i] -= 2*dx
166        deriv[i] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
167        Oxyz[i] += dx
168        Txyz[i] += dx
169        d0 = calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)
170        Txyz[i] -= 2*dx
171        deriv[i+3] = (calcDist(Oxyz,Txyz,Tunit,inv,C,M,T,Amat)-d0)/(2.*dx)
172        Txyz[i] += dx
173    return deriv
174   
175def getAngSig(VA,VB,Amat,SGData,covData={}):
176   
177    def calcVec(Ox,Tx,U,inv,C,M,T,Amat):
178        TxT = inv*(np.inner(M,Tx)+T)+C
179        TxT = G2spc.MoveToUnitCell(TxT)+U
180        return np.inner(Amat,(TxT-Ox))
181       
182    def calcAngle(Ox,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat):
183        VecA = calcVec(Ox,TxA,unitA,invA,CA,MA,TA,Amat)
184        VecA /= np.sqrt(np.sum(VecA**2))
185        VecB = calcVec(Ox,TxB,unitB,invB,CB,MB,TB,Amat)
186        VecB /= np.sqrt(np.sum(VecB**2))
187        edge = VecB-VecA
188        edge = np.sum(edge**2)
189        angle = (2.-edge)/2.
190        angle = max(angle,-1.)
191        return acosd(angle)
192       
193    OxAN,OxA,TxAN,TxA,unitA,TopA = VA
194    OxBN,OxB,TxBN,TxB,unitB,TopB = VB
195    invA = invB = 1
196    invA = TopA/abs(TopA)
197    invB = TopB/abs(TopB)
198    centA = abs(TopA)/100
199    centB = abs(TopB)/100
200    opA = abs(TopA)%100-1
201    opB = abs(TopB)%100-1
202    MA,TA = SGData['SGOps'][opA]
203    MB,TB = SGData['SGOps'][opB]
204    CA = SGData['SGCen'][centA]
205    CB = SGData['SGCen'][centB]
206    if 'covMatrix' in covData:
207        covMatrix = covData['covMatrix']
208        varyList = covData['varyList']
209        AngVcov = getVCov(OxAN+TxAN+TxBN,varyList,covMatrix)
210        dx = .00001
211        dadx = np.zeros(9)
212        Ang = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
213        for i in [0,1,2]:
214            OxA[i] += dx
215            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
216            OxA[i] -= 2*dx
217            dadx[i] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
218            OxA[i] += dx
219           
220            TxA[i] += dx
221            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
222            TxA[i] -= 2*dx
223            dadx[i+3] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
224            TxA[i] += dx
225           
226            TxB[i] += dx
227            a0 = calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)
228            TxB[i] -= 2*dx
229            dadx[i+6] = (calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat)-a0)/dx
230            TxB[i] += dx
231           
232        sigAng = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
233        if sigAng < 0.01:
234            sigAng = 0.0
235        return Ang,sigAng
236    else:
237        return calcAngle(OxA,TxA,TxB,unitA,unitB,invA,CA,MA,TA,invB,CB,MB,TB,Amat),0.0
238
239def GetDistSig(Oatoms,Atoms,Amat,SGData,covData={}):
240
241    def calcDist(Atoms,SyOps,Amat):
242        XYZ = []
243        for i,atom in enumerate(Atoms):
244            Inv,M,T,C,U = SyOps[i]
245            XYZ.append(np.array(atom[1:4]))
246            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
247            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
248        V1 = XYZ[1]-XYZ[0]
249        return np.sqrt(np.sum(V1**2))
250       
251    Inv = []
252    SyOps = []
253    names = []
254    for i,atom in enumerate(Oatoms):
255        names += atom[-1]
256        Op,unit = Atoms[i][-1]
257        inv = Op/abs(Op)
258        m,t = SGData['SGOps'][abs(Op)%100-1]
259        c = SGData['SGCen'][abs(Op)/100]
260        SyOps.append([inv,m,t,c,unit])
261    Dist = calcDist(Oatoms,SyOps,Amat)
262   
263    sig = -0.001
264    if 'covMatrix' in covData:
265        parmNames = []
266        dx = .00001
267        dadx = np.zeros(6)
268        for i in range(6):
269            ia = i/3
270            ix = i%3
271            Oatoms[ia][ix+1] += dx
272            a0 = calcDist(Oatoms,SyOps,Amat)
273            Oatoms[ia][ix+1] -= 2*dx
274            dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)
275        covMatrix = covData['covMatrix']
276        varyList = covData['varyList']
277        DistVcov = getVCov(names,varyList,covMatrix)
278        sig = np.sqrt(np.inner(dadx,np.inner(DistVcov,dadx)))
279        if sig < 0.001:
280            sig = -0.001
281   
282    return Dist,sig
283
284def GetAngleSig(Oatoms,Atoms,Amat,SGData,covData={}):
285       
286    def calcAngle(Atoms,SyOps,Amat):
287        XYZ = []
288        for i,atom in enumerate(Atoms):
289            Inv,M,T,C,U = SyOps[i]
290            XYZ.append(np.array(atom[1:4]))
291            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
292            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
293        V1 = XYZ[1]-XYZ[0]
294        V1 /= np.sqrt(np.sum(V1**2))
295        V2 = XYZ[1]-XYZ[2]
296        V2 /= np.sqrt(np.sum(V2**2))
297        V3 = V2-V1
298        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
299        return acosd(cang)
300
301    Inv = []
302    SyOps = []
303    names = []
304    for i,atom in enumerate(Oatoms):
305        names += atom[-1]
306        Op,unit = Atoms[i][-1]
307        inv = Op/abs(Op)
308        m,t = SGData['SGOps'][abs(Op)%100-1]
309        c = SGData['SGCen'][abs(Op)/100]
310        SyOps.append([inv,m,t,c,unit])
311    Angle = calcAngle(Oatoms,SyOps,Amat)
312   
313    sig = -0.01
314    if 'covMatrix' in covData:
315        parmNames = []
316        dx = .00001
317        dadx = np.zeros(9)
318        for i in range(9):
319            ia = i/3
320            ix = i%3
321            Oatoms[ia][ix+1] += dx
322            a0 = calcAngle(Oatoms,SyOps,Amat)
323            Oatoms[ia][ix+1] -= 2*dx
324            dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)
325        covMatrix = covData['covMatrix']
326        varyList = covData['varyList']
327        AngVcov = getVCov(names,varyList,covMatrix)
328        sig = np.sqrt(np.inner(dadx,np.inner(AngVcov,dadx)))
329        if sig < 0.01:
330            sig = -0.01
331   
332    return Angle,sig
333
334def GetTorsionSig(Oatoms,Atoms,Amat,SGData,covData={}):
335   
336    def calcTorsion(Atoms,SyOps,Amat):
337       
338        XYZ = []
339        for i,atom in enumerate(Atoms):
340            Inv,M,T,C,U = SyOps[i]
341            XYZ.append(np.array(atom[1:4]))
342            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
343            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
344        V1 = XYZ[1]-XYZ[0]
345        V2 = XYZ[2]-XYZ[1]
346        V3 = XYZ[3]-XYZ[2]
347        V1 /= np.sqrt(np.sum(V1**2))
348        V2 /= np.sqrt(np.sum(V2**2))
349        V3 /= np.sqrt(np.sum(V3**2))
350        M = np.array([V1,V2,V3])
351        D = nl.det(M)
352        Ang = 1.0
353        P12 = np.dot(V1,V2)
354        P13 = np.dot(V1,V3)
355        P23 = np.dot(V2,V3)
356        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
357        return Tors
358           
359    Inv = []
360    SyOps = []
361    names = []
362    for i,atom in enumerate(Oatoms):
363        names += atom[-1]
364        Op,unit = Atoms[i][-1]
365        inv = Op/abs(Op)
366        m,t = SGData['SGOps'][abs(Op)%100-1]
367        c = SGData['SGCen'][abs(Op)/100]
368        SyOps.append([inv,m,t,c,unit])
369    Tors = calcTorsion(Oatoms,SyOps,Amat)
370   
371    sig = -0.01
372    if 'covMatrix' in covData:
373        parmNames = []
374        dx = .00001
375        dadx = np.zeros(12)
376        for i in range(12):
377            ia = i/3
378            ix = i%3
379            Oatoms[ia][ix+1] += dx
380            a0 = calcTorsion(Oatoms,SyOps,Amat)
381            Oatoms[ia][ix+1] -= 2*dx
382            dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
383        covMatrix = covData['covMatrix']
384        varyList = covData['varyList']
385        TorVcov = getVCov(names,varyList,covMatrix)
386        sig = np.sqrt(np.inner(dadx,np.inner(TorVcov,dadx)))
387        if sig < 0.01:
388            sig = -0.01
389   
390    return Tors,sig
391       
392def GetDATSig(Oatoms,Atoms,Amat,SGData,covData={}):
393   
394    def calcDist(Atoms,SyOps,Amat):
395        XYZ = []
396        for i,atom in enumerate(Atoms):
397            Inv,M,T,C,U = SyOps[i]
398            XYZ.append(np.array(atom[1:4]))
399            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
400            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
401        V1 = XYZ[1]-XYZ[0]
402        return np.sqrt(np.sum(V1**2))
403       
404    def calcAngle(Atoms,SyOps,Amat):
405        XYZ = []
406        for i,atom in enumerate(Atoms):
407            Inv,M,T,C,U = SyOps[i]
408            XYZ.append(np.array(atom[1:4]))
409            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
410            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
411        V1 = XYZ[1]-XYZ[0]
412        V1 /= np.sqrt(np.sum(V1**2))
413        V2 = XYZ[1]-XYZ[2]
414        V2 /= np.sqrt(np.sum(V2**2))
415        V3 = V2-V1
416        cang = min(1.,max((2.-np.sum(V3**2))/2.,-1.))
417        return acosd(cang)
418
419    def calcTorsion(Atoms,SyOps,Amat):
420       
421        XYZ = []
422        for i,atom in enumerate(Atoms):
423            Inv,M,T,C,U = SyOps[i]
424            XYZ.append(np.array(atom[1:4]))
425            XYZ[-1] = Inv*(np.inner(M,np.array(XYZ[-1]))+T)+C+U
426            XYZ[-1] = np.inner(Amat,XYZ[-1]).T
427        V1 = XYZ[1]-XYZ[0]
428        V2 = XYZ[2]-XYZ[1]
429        V3 = XYZ[3]-XYZ[2]
430        V1 /= np.sqrt(np.sum(V1**2))
431        V2 /= np.sqrt(np.sum(V2**2))
432        V3 /= np.sqrt(np.sum(V3**2))
433        M = np.array([V1,V2,V3])
434        D = nl.det(M)
435        Ang = 1.0
436        P12 = np.dot(V1,V2)
437        P13 = np.dot(V1,V3)
438        P23 = np.dot(V2,V3)
439        Tors = acosd((P12*P23-P13)/(np.sqrt(1.-P12**2)*np.sqrt(1.-P23**2)))*D/abs(D)
440        return Tors
441           
442    Inv = []
443    SyOps = []
444    names = []
445    for i,atom in enumerate(Oatoms):
446        names += atom[-1]
447        Op,unit = Atoms[i][-1]
448        inv = Op/abs(Op)
449        m,t = SGData['SGOps'][abs(Op)%100-1]
450        c = SGData['SGCen'][abs(Op)/100]
451        SyOps.append([inv,m,t,c,unit])
452    M = len(Oatoms)
453    if M == 2:
454        Val = calcDist(Oatoms,SyOps,Amat)
455    elif M == 3:
456        Val = calcAngle(Oatoms,SyOps,Amat)
457    else:
458        Val = calcTorsion(Oatoms,SyOps,Amat)
459   
460    sigVals = [-0.001,-0.01,-0.01]
461    sig = sigVals[M-3]
462    if 'covMatrix' in covData:
463        parmNames = []
464        dx = .00001
465        N = M*3
466        dadx = np.zeros(N)
467        for i in range(N):
468            ia = i/3
469            ix = i%3
470            Oatoms[ia][ix+1] += dx
471            if M == 2:
472                a0 = calcDist(Oatoms,SyOps,Amat)
473            elif M == 3:
474                a0 = calcAngle(Oatoms,SyOps,Amat)
475            else:
476                a0 = calcTorsion(Oatoms,SyOps,Amat)
477            Oatoms[ia][ix+1] -= 2*dx
478            if M == 2:
479                dadx[i] = (calcDist(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
480            elif M == 3:
481                dadx[i] = (calcAngle(Oatoms,SyOps,Amat)-a0)/(2.*dx)               
482            else:
483                dadx[i] = (calcTorsion(Oatoms,SyOps,Amat)-a0)/(2.*dx)
484        covMatrix = covData['covMatrix']
485        varyList = covData['varyList']
486        Vcov = getVCov(names,varyList,covMatrix)
487        sig = np.sqrt(np.inner(dadx,np.inner(Vcov,dadx)))
488        if sig < sigVals[M-3]:
489            sig = sigVals[M-3]
490   
491    return Val,sig
492       
493   
494def ValEsd(value,esd=0,nTZ=False):                  #NOT complete - don't use
495    # returns value(esd) string; nTZ=True for no trailing zeros
496    # use esd < 0 for level of precision shown e.g. esd=-0.01 gives 2 places beyond decimal
497    #get the 2 significant digits in the esd
498    edig = lambda esd: int(round(10**(math.log10(esd) % 1+1)))
499    #get the number of digits to represent them
500    epl = lambda esd: 2+int(1.545-math.log10(10*edig(esd)))
501   
502    mdec = lambda esd: -int(round(math.log10(abs(esd))))+1
503    ndec = lambda esd: int(1.545-math.log10(abs(esd)))
504    if esd > 0:
505        fmt = '"%.'+str(ndec(esd))+'f(%d)"'
506        return str(fmt%(value,int(round(esd*10**(mdec(esd)))))).strip('"')
507    elif esd < 0:
508         return str(round(value,mdec(esd)-1))
509    else:
510        text = str("%f"%(value))
511        if nTZ:
512            return text.rstrip('0')
513        else:
514            return text
515
516   
Note: See TracBrowser for help on using the repository browser.