Multivariate SPDE¶
In [1]:
import gstlearn as gl
import scipy as sc
import numpy as np
import scipy.sparse.linalg
In [2]:
multistruct = True
### Stationary case
mesh1 = gl.MeshETurbo([10,10])
r1 = 0.9
s11 = 2
s21 = 1
s121 = r1 * np.sqrt(s11*s21)
modelMulti = gl.Model.createFromParam(gl.ECov.MATERN,param=1,sills = [s11,s121,s121,s21],range = 5)
mesh2 = gl.MeshETurbo([20,20])
r2 = -0.
s12 = .01
s22 = .04
s122 = r2 * np.sqrt(s12*s22)
modelMulti2 = gl.Model.createFromParam(gl.ECov.MATERN,param=1,sills = [s12,s122,s122,s22],range=5)
if multistruct:
modelMulti.addCov(modelMulti2.getCova(0))
In [3]:
class PrecisionOpMulti:
def __init__(self,model,meshes):
self.sills = []
self.invsill = []
self.cholsill = []
self.nvertex = []
self.SigmaMult = []
self.invSigmaMult = []
self.Qop = []
self.temp = []
self.ncovar = model.getCovaNumber()
for i in range(self.ncovar):
cova = model.getCova(i)
self.sills += [cova.getSill().toTL()]
self.invsill += [np.linalg.inv(self.sills[i])]
self.cholsill += [np.linalg.cholesky(self.sills[i])]
self.nvertex += [meshes[i].getNApices()]
modelMono = gl.Model.createFromParam(cova.getType(),
param = cova.getParam(),
range = cova.getRange())
self.Qop += [gl.PrecisionOp(meshes[i],modelMono)]
self.temp += [gl.VectorDouble(np.zeros(shape=self.nvertex[i]))]
# Quantité utilisées pour coder la ref. Inutile de le faire en C++, c'est juste pour tester
QopCs = gl.PrecisionOpCs(meshes[i],modelMono)
Qmat = QopCs.getQ().toTL()
Sigma = sc.sparse.linalg.inv(Qmat).todense()
self.SigmaMult += [np.kron(self.sills[i],Sigma)]
self.invSigmaMult += [np.linalg.inv(self.SigmaMult[i])]
self.nvar = self.invsill[0].shape[0]
self.nvertextot = np.sum(self.nvertex)
self.sizetot = self.nvertextot * self.nvar
# Fonction pour évaluer les refs. Inutile de la coder en C++
def evalRefDirect(self,u):
return np.array(self.invSigmaMult@u).reshape(-1)
# Fonction pour évaluer les refs. Inutile de la coder en C++
def evalRefSigma(self,u):
return np.array(self.SigmaMult@u).reshape(-1)
# Fonction pour évaluer les refs. Inutile de la coder en C++
def evalQinvByQmhTimesQmhTr(self,gauss):
v = gl.VectorDouble(self.nvertex)
tempsimu = np.zeros(shape=self.nvertex * self.nvar)
for i in range(self.nvar):
temp = gl.VectorDouble(gauss[(i*self.nvertex):((i+1)*self.nvertex)])
self.Qop.getShiftOp().prodLambda(temp,v,gl.EPowerPT.MINUSONE)
tempsimu[(i*self.nvertex):((i+1)*self.nvertex)] = v.getVector()
w = self.evalSimulateTranspose(tempsimu)
for i in range(self.nvar):
temp = gl.VectorDouble(w[(i*self.nvertex):((i+1)*self.nvertex)])
self.Qop.getShiftOp().prodLambda(temp,v,gl.EPowerPT.ONE)
tempsimu[(i*self.nvertex):((i+1)*self.nvertex)] = v.getVector()
w = self.evalSimulate(tempsimu)
return w
# Privilégier un proto evalDirect(const Vector inv&, Vector outv&)
def evalDirect(self,inv):
outv = np.zeros_like(inv)
s=0
for ii in range(self.ncovar):
for i in range(self.nvar):
self.Qop[ii].evalDirect(inv[(s+self.nvertex[ii]*i):(s+self.nvertex[ii]*(i+1))],self.temp[ii])
for j in range(self.nvar):
outv[(s+self.nvertex[ii]*j):(s+self.nvertex[ii]*(j+1))]+=self.invsill[ii][i,j] * \
np.array(self.temp[ii].getVector())
s+= self.nvertex[ii] * self.nvar
return outv
def evalSimulate(self,gauss):
outv = np.zeros_like(gauss)
s=0
iadx = 0
for ii in range(self.ncovar):
nv = self.nvertex[ii]
for i in range(self.nvar):
self.Qop[ii].evalSimulate(
gl.VectorDouble(gauss[iadx:(iadx+nv)]),
self.temp[ii])
iady = s
for j in range(0,self.nvar):
outv[iady:(iady+nv)]+=self.cholsill[ii][j,i] *\
np.array(self.temp[ii].getVector())
iady+=nv
iadx += nv
s = iadx
return outv
def getCholSill(self):
return self.cholsill
#Ne pas coder en C++, c'est juste pour un test
def evalSimulateTranspose(self,gauss):
outv = np.zeros_like(gauss)
for i in range(self.nvar):
self.Qop.evalSimulate(gl.VectorDouble(gauss[(self.nvertex*i):(self.nvertex*(i+1))]),self.temp)
for j in range(0,i+1):
outv[(self.nvertex*j):(self.nvertex*(j+1))]+=self.cholsill[j,i] *\
np.array(self.temp.getVector())
return outv
In [4]:
pops = gl.PrecisionOpMulti(modelMulti)
pops.addMesh(mesh1)
if multistruct:
pops.addMesh(mesh2)
pops.display()
Number of Variables = 2 Number of Covariances = 2 Number of Meshes = 2 Vector dimension = 1000 Indices of Matérn Covariance = 0 1 Dimensions of the Meshes = 100 400 Class is Valid for operations!
In [5]:
meshes = [mesh1]
if multistruct:
meshes += [mesh2]
pop = PrecisionOpMulti(modelMulti, [mesh1,mesh2])
nvar = pop.nvar
u = np.random.normal(size=pop.sizetot )
res = pops.evalDirect(u)
ref = pop.evalDirect(u)
In [6]:
print("Error in the operator product =",np.max(np.abs(res-ref)))
Error in the operator product = 0.0
In [7]:
import matplotlib.pyplot as plt
ind = range(pop.nvertextot*pop.nvar)
ax = plt.scatter(res[ind],ref[ind])
In [8]:
resnew = pops.evalDirect(u)
np.max(np.abs(resnew-ref))
Out[8]:
0.0
In [9]:
simref = pop.evalSimulate(u)
simnew = pops.evalSimulate(u)
In [10]:
print("Error in the simulation operator =",np.max(np.abs(simnew-simref)))
Error in the simulation operator = 0.0
In [11]:
import matplotlib.pyplot as plt
ind = range(pop.nvertextot,pop.sizetot)
ind = range(pop.sizetot)
ax = plt.scatter(simnew[ind],simref[ind])