Multivariate SPDE¶

import gstlearn as gl
import scipy as sc
import numpy as np
import scipy.sparse.linalg

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))

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
        
    

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)

simref = pop.evalSimulate(u)
simnew = pops.evalSimulate(u)