Source code for PyMca5.PyMcaMath.mva.Lanczos

## /************************************************************************
##
##   Copyright
##   Alessandro MIRONE
##   mirone@esrf.fr
##
##   Copyright 2002  by European Synchrotron Radiation Facility, Grenoble,
##                   France
##
##                                ----------
##
##                            All Rights Reserved
##
##                                ----------
##
## Permission to use, copy, modify, and distribute this software and its
## documentation for any purpose and without fee is hereby granted,
## provided that the above copyright notice appear in all copies and that
## both that copyright notice and this permission notice appear in
## supporting documentation, and that the names of European Synchrotron
## Radiation Facility or ESRF or SCISOFT not be used in advertising or
## publicity pertaining to distribution of the software without specific,
## written prior permission.
##
## EUROPEAN SYNCHROTRON RADIATION FACILITY DISCLAIMS ALL WARRANTIES WITH
## REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL EUROPEAN SYNCHROTRON
## RADIATION FACILITY OR ESRF BE LIABLE FOR ANY SPECIAL, INDIRECT OR
## CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE,
## DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
## TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
## PERFORMANCE OF THIS SOFTWARE.
##
## **************************************************************************/
__author__ = "A. Mirone - ESRF SciSoft Group"
__license__ = "MIT"
import math
sparsamodulo=0
PARALLEL=0

import numpy
import numpy.linalg
try:
    import numpy.core._dotblas as dotblas
except ImportError:
    dotblas = numpy
import random

######################################
## interfaccia per sparse
# __class__
# static    .getClass4Vect()
# object    .gohersch()
# object    .goherschMax()
# object    .goherschMin()

######################################
## interfaccia per vettori
# __class__ === class4vect
# class Vector(dim)
# object.set_value(n,v)
# object.set_all_random(v)

[docs]class LanczosNumericMatrix(object): tipo=numpy.float64 def __init__(self, mR): self.mR=mR self.dim= self.mR[0].shape[0] self.shift=0.0
[docs] def Moltiplica(self,res,v): if( len(self.mR)==1 ): res.vr[:]+=dotblas.dot(self.mR[0],v.vr) else: res.vr[:]+=dotblas.dot(self.mR[0],dotblas.dot(self.mR[1],v.vr)) if( self.shift !=0.0): res.add_from_vect_with_fact(v,self.shift)
[docs] def trasforma(self, fattore, addendo): self.shift+=addendo if( fattore!=1): self.mR[0][:]=self.mR[0]*numpy.array([fattore], self.tipo)
[docs] def getClass4Vect(self): return LanczosNumericVector
[docs]class LanczosNumericVector(object): tipo=numpy.float64 def __init__(self, *dim): if( dim!=(0,) ): self.vr=numpy.zeros(dim,self.tipo) else: pass def __getitem__(self, i): res=LanczosNumericVector(0) res.vr=self.vr[i] return res
[docs] def mat_mult(self, evect , q): self.vr[:evect.shape[0]]=dotblas.dot(evect.astype(self.tipo),q.vr[:evect.shape[1]])
[docs] def dividebyarray(self,prec): self.vr[:]=self.vr/prec
[docs] def multbyarray(self,prec): self.vr[:]=self.vr*prec
[docs] def len(self): return len(self.vr)
[docs] def copy(self): duma=numpy.array(self.vr) res=numpy.Vector(0) res.vr=duma return res
[docs] def copy_to_a_from_b(self,b): self.vr[:]=b.vr
[docs] def set_value(self, n, val): self.vr[n]=val
[docs] def set_to_zero(self): self.vr[:]=0
[docs] def set_to_one(self): self.vr[:]=1
[docs] def set_all_random(self, v): self.vr[:]=[random.random() for i in range(len(self.vr)) ]
[docs] def scalare(self,b ): resR =numpy.sum(self.vr*b.vr) return resR
[docs] def sqrtscalare(self,b): assert( self is b) return numpy.sqrt( self.scalare(b) )
[docs] def normalizzaauto(self): norma = self.sqrtscalare(self) self.vr[:]=self.vr/norma
[docs] def normalizza(self,norma): self.vr.normalizza(norma)
[docs] def rescale(self,fact): if(fact==0.0): self.set_to_zero() else: norma=1.0/fact self.normalizza(norma)
[docs] def add_from_vect(self, b): self.vr[:]=self.vr+b.vr
[docs] def add_from_vect_with_fact(self, b, fact): self.vr[:]=self.vr+numpy.array([fact],self.tipo)*b.vr
[docs]def REAL(a): if( type(a)==type(1.0+1.0j) ): return a.real else: return a
[docs]def Real(x): if( x.__class__ == complex): return x.real else: return x
[docs]class Lanczos: dump_count=0; countdumpab=0 def __init__(self, sparse, metrica=None, tol=1.0e-15): self.matrice=sparse self.metrica=metrica self.class4sparse = sparse.__class__ try: self.class4vect = self.class4sparse.getClass4Vect() except: self.class4vect = sparse.getClass4Vect() self.nsteps = 0 self.dim=self.matrice.dim self.q = None self.alpha = None self.beta = None self.omega=None self.tol = tol self.maxIt = 50 self.old = False
[docs] def diagoCustom(self, minDim=5, shift=None):
if shift is None: self.matrice.gohersch() shift = - self.matrice.goherschMax() self.cerca(minDim, shift) for ne in range(len(self.eval)): self.eval[ne] = self.eval[ne] - shift
[docs] def passeggia(self, k, m, start, gram=0, NT=4):
if k<0 or m>self.nsteps: print("k = ",k," m = ",m," nsteps = ",self.nsteps) print("Something wrong in passe k<0 or m>nsteps") raise ValueError(\ "Lanczos. Something wrong in passe k<0 or m>nsteps") sn = math.sqrt(float(self.dim)) eu = 1.1e-16 eusn = eu*sn if k==0: self.class4vect.copy_to_a_from_b(self.q[0],start) if self.metrica is None: self.q[0].normalizzaauto() else: self.tmp4met.set_to_zero() self.metrica.Moltiplica( self.tmp4met , self.q[0] ) tmpfat = math.sqrt(REAL(self.class4vect.scalare(self.q[0] , self.tmp4met))) self.q[0].normalizza(tmpfat) # self.q[0].dumptofile("q_0") p= self.class4vect(self.dim) for i in range(k,m): p.set_to_zero() # self.q[i].dumptofile("qbef_"+str(self.dump_count) ) self.matrice.Moltiplica(p,self.q[i]) if self.metrica is not None: self.tmp4met.set_to_zero() self.metrica.Moltiplica(self.tmp4met , self.q[i]) self.alpha[i] = REAL(self.class4vect.scalare(p , self.tmp4met)) else: self.alpha[i] = REAL(self.class4vect.scalare(p , self.q[i])) # p.dumptofile("p"+str(self.dump_count) ) self.dump_count+=1 if i==k: p.add_from_vect_with_fact( self.q[k] , -self.alpha[k] ) for l in range(k): p.add_from_vect_with_fact( self.q[l] , -self.beta[l] ) else: # self.q[i].dumptofile("q_i"+str(self.dump_count) ) # self.q[i-1].dumptofile("q_im1"+str(self.dump_count) ) p.add_from_vect_with_fact(self.q[i] , -self.alpha[i] ) p.add_from_vect_with_fact(self.q[i-1] , -self.beta[i-1] ) # p.dumptofile("pp"+str(self.dump_count) ) self.dump_count+=1 if self.metrica is not None: self.tmp4met.set_to_zero() self.metrica.Moltiplica(self.tmp4met , p) self.beta[i] = math.sqrt(REAL(self.class4vect.scalare(p , self.tmp4met))) else: self.beta[i]=self.class4vect.sqrtscalare(p,p) self.omega[i,i]=1. max0 = 0.0 if self.beta[i] != 0: #and not scipy.isnan(self.beta[i]): for j in range(i+1): self.omega[i+1,j] = eusn if j<k: add = 2 * eusn + abs(self.alpha[j]-self.alpha[i]) \ * abs(self.omega[i,j]) if i!=k: add += self.beta[j]*abs(self.omega[i,k]) if i>0 and j!=i-1: add += self.beta[i-1]*abs(self.omega[i-1,j]) self.omega[i+1,j] += add / self.beta[i] elif j==k : add = 2 * eusn + abs(self.alpha[j]-self.alpha[i]) \ * abs(self.omega[i,j]) for w in range(k): add += self.beta[w]*abs(self.omega[i,w]) if i!=k+1: add += self.beta[k]*abs(self.omega[i,k+1]) if i>0 and i!=k+1: add += self.beta[i-1]*abs(self.omega[i-1,k]) self.omega[i+1,j] += add / self.beta[i] elif j<i : add = 2 * eusn + abs(self.alpha[j]-self.alpha[i]) \ * abs(self.omega[i,j]) if i!=j+1: add += self.beta[j]*abs(self.omega[i,j+1]) if i>0 and j>0: add += self.beta[j-1]*abs(self.omega[i-1,j-1]) if i>0 and i!=j+1: add += self.beta[i-1]*abs(self.omega[i-1,j]) self.omega[i+1,j] += add / self.beta[i] else: add = eusn if i>0: add += self.beta[i-1]*abs(self.omega[i,i-1] ) self.omega[i+1,j] += add / self.beta[i] self.omega[j,i+1] = self.omega[i+1,j] max0 += self.omega[i+1,j]**2 if self.beta[i]==0 or max0>eu or gram : ## if self.beta[i]==0 or max0>0. : if i>0: #print " GRAMO self.beta[i]==0 or max0>eu and i>0", i," ", self.dump_count self.GramSchmidt(self.q[i],i, NT) if self.metrica is None: self.q[i].normalizzaauto() else: self.tmp4met.set_to_zero() self.metrica.Moltiplica( self.tmp4met , self.q[i] ) tmpfat = math.sqrt(REAL(self.class4vect.scalare(self.q[i] , self.tmp4met))) self.q[i].normalizza(tmpfat) p.set_to_zero() self.matrice.Moltiplica(p, self.q[i]) if self.metrica is not None: self.tmp4met.set_to_zero() self.metrica.Moltiplica(self.tmp4met , self.q[i]) self.alpha[i] = REAL(self.class4vect.scalare(p , self.tmp4met)) else: self.alpha[i] = REAL(self.class4vect.scalare(p , self.q[i])) #print " GRAMO bis ", i ## print self.alpha[:20] ## raise " OK " if self.metrica is None: self.GramSchmidt(p,i+1,NT) self.beta[i] = self.class4vect.sqrtscalare(p,p) p.normalizzaauto() else: # p.add_from_vect_with_fact( self.q[i] , -self.alpha[i] ) # p.add_from_vect_with_fact( self.q[i-1] , -self.beta[i-1] ) self.GramSchmidt(p,i+1,NT) self.tmp4met.set_to_zero() self.metrica.Moltiplica(self.tmp4met , p) tmpfat = math.sqrt(REAL(self.class4vect.scalare(p , self.tmp4met))) self.beta[i] = tmpfat p.normalizza(tmpfat) if i>0: condition = eu * \ math.sqrt(self.dim * (self.alpha[i]**2+\ self.beta[i-1]**2)) else: condition = eu * math.sqrt(self.dim * (self.alpha[i]**2)) if self.beta[i]< condition: self.beta[i]=0. p.set_all_random(1.0) self.GramSchmidt(p,i+1) if self.metrica is None: p.normalizzaauto() else: self.tmp4met.set_to_zero() self.metrica.Moltiplica( self.tmp4met , p ) tmpfat = math.sqrt(REAL(self.class4vect.scalare(p , self.tmp4met))) p.normalizza(tmpfat) for l in range(i) : self.omega[i,l]=self.omega[l,i]=eusn for l in range(i+1): self.omega[i+1,l]=self.omega[l,i+1]=eusn else: if self.metrica is None: p.normalizzaauto() else: self.tmp4met.set_to_zero() self.metrica.Moltiplica( self.tmp4met , p ) tmpfat = math.sqrt(REAL(self.class4vect.scalare(p , self.tmp4met))) p.normalizza(tmpfat) # p.dumptofile("normprelude"+str(self.dump_count)) self.class4vect.copy_to_a_from_b(self.q[i+1],p)
[docs] def GramSchmidt(self, vect , n, NT=4): for h in range(NT): if self.metrica is None : for i in range(n): s=self.class4vect.scalare( self.q[i], vect) vect.add_from_vect_with_fact(self.q[i],-s) else: self.tmp4met.set_to_zero() self.metrica.Moltiplica(self.tmp4met, vect) for i in range(n): s=self.class4vect.scalare( self.q[i], self.tmp4met) vect.add_from_vect_with_fact(self.q[i],-s)
[docs] def allocaMemory(self): self.alpha = numpy.zeros(self.nsteps,numpy.float64) self.beta = numpy.zeros(self.nsteps,numpy.float64) self.omega = numpy.zeros((self.nsteps+1,self.nsteps+1),numpy.float64) self.evect = numpy.zeros((self.nsteps,self.nsteps),numpy.float64) self.eval = numpy.zeros((self.nsteps),numpy.float64) self.oldalpha = numpy.zeros((self.nsteps),numpy.float64) self.q=self.class4vect(self.nsteps+1,self.dim) if self.metrica is not None: self.tmp4met = self.class4vect( self.dim )
[docs] def cerca(self, nd, shift):
self.shift=shift self.matrice.trasforma(1.0, shift) m = min(4*nd, self.dim) self.nsteps = m self.allocaMemory() vect_init = self.class4vect(self.dim) # vect_init.set_value(0,1.0) vect_init.set_all_random(1.0) # vect_init.set_to_one() k=0 nc=0 self.passeggia(k,m,vect_init) while nc<nd : #print " DIAGONALIZZAZIONE " self.diago(k,m) nc = self.converged(m) if k and not numpy.sometrue(abs(self.beta[:k])>self.tol) : break if (nc+2*nd) >= m: k=m-1 else: k=nc+2*nd #print "KKKKKKKKK " , k self.ricipolla(k,m) self.countdumpab+=1 #print " k,m , dim", k, m, self.dim # return m # sentinell self.passeggia(k,m,self.q[m]) if m==self.dim: return m else: return k
[docs] def converged(self, m):
for j in range(1,m): for i in range(m-j): if abs(self.eval[i]) < abs(self.eval[i+1]): self.eval[i],self.eval[i+1]=self.eval[i+1],self.eval[i] self.evect[i:i+2]=self.evect[i+1],self.evect[i] # print "eval", self.eval[0:m] res = 0 if self.old: while res<m: #print self.oldalpha[res] if (abs(self.eval[res]-self.oldalpha[res])\ /abs(self.oldalpha[res])>self.tol): break res+=1 else: self.old=True self.oldalpha[0:m]=self.eval[0:m] return res
[docs] def ricipolla(self, k,m):
#print "k,m",k,m self.alpha[0:k]= self.eval[0:k] self.beta[0:k] = self.beta[m-1]*self.evect[0:k,m-1] #print "beta beta",self.beta,self.beta[m-1],self.evect[0:k,m-1] # E = self.evect[0:k,0:m].copy() a=self.class4vect(k, self.dim ) #print " LANCIO mat MULT " self.class4vect.mat_mult(a, self.evect[0:k,0:m] , self.q) for i in range(k): a[i].normalizzaauto() self.class4vect.copy_to_a_from_b( self.q[i] , a[i] ) a=None self.class4vect.copy_to_a_from_b(self.q[k] , self.q[m] ) o = self.omega[0:m,0:m].copy() o = dotblas.dot(o, numpy.transpose(self.evect)) for i in range(k): self.omega[i,k]=self.omega[k,i]=o[i,k] o = dotblas.dot(self.evect,o) self.omega[0:k,0:k]=o[0:k,0:k]
[docs] def diago(self, k, m): mat = numpy.zeros([m,m], numpy.float) mat.shape=[m*m] mat[0:m*m:m+1] = self.alpha mat[k*m+k+1:m*m:m+1] =self.beta[k:m-1] mat[(k+1)*m+k:m*m:m+1] =self.beta[k:m-1] mat.shape=[m,m] mat[ k ,0:k ] = self.beta[:k] mat[ 0:k, k ] = self.beta[:k] self.eval,self.evect = numpy.linalg.eigh(mat)
[docs]def solveEigenSystem( S_base , nsearchedeigen, shift=None, metrica=None, tol=1.0e-15): lnczs=Lanczos( S_base , metrica=metrica, tol=tol) lnczs.diagoCustom(minDim=nsearchedeigen, shift=shift) return lnczs.eval, lnczs.q