Spatio Temporal with gstlearn

In this tutorial, we study the family of SPDEs :

$$\left(\frac{\partial}{\partial t} + c\left[(1-\nabla H. \nabla)^{\alpha} + v.\nabla\right]\right)Z(s,t)=\sqrt{c}W_T(t)\otimes X_S(s)$$

where

• $c>0$ is a time scale parameter
• $H$ is an anisotropic matrix
• $v$ is a velocity vector
• $W_T$ is a temporal white-noise

• $X_S$ is a spatially structured noise, solution of

  $$(1-\nabla H_S \nabla)X_S = \mathcal{W}_S$$

import numpy as np
import gstlearn as gl
from scipy.special import gamma

sudo apt-get install python-scipy libsuitesparse-dev git clone https://github.com/scikit-sparse/scikit-sparse.git cd scikit-sparse python3 setup.py install

from sksparse.cholmod import cholesky
import scipy as sc
from scipy.sparse import diags

scale = 4.
kappa = 1./scale
kappa2 = kappa**2
dt = .1
dx = 1. 
c = 40
sqc = np.sqrt(c)
sqdt = np.sqrt(dt)

nx = [100,100]
m = gl.Model.createFromParam(gl.ECov.BESSEL_K,range=1,param=2,flagRange=False)
mesh = gl.MeshETurbo(nx ,[dx,dx])
S = gl.ShiftOpCs(mesh,m)
Smat = S.getSToTriplet().toTL()
TildeC12 = diags(np.sqrt(S.getTildeC()))
invTildeC12 = diags(1./np.sqrt(S.getTildeC()))
G = TildeC12 @ Smat @ TildeC12
M = TildeC12 @ TildeC12
K = (kappa2 * M + G) @ (kappa2 * M + G) @ (kappa2 * M + G)
param = 3
P = M +  c*dt * K
cholP = cholesky(P)

def evalInvA(x) :
    return  cholP.solve_A(M @ x)

def evalInvB(x) :
    return sqc * sqdt * cholP.solve_A(TildeC12 @ x)

# TODO: Reason for the warning ??!
# https://stackoverflow.com/a/48273717

/tmp/ipykernel_729697/2651008202.py:13: CholmodTypeConversionWarning: converting matrix of class csr_matrix to CSC format
  cholP = cholesky(P)

np.random.seed(123)
x = np.random.normal(size=Smat.shape[0])

for i in range(100):
    u = np.random.normal(size=Smat.shape[0])
    x = evalInvA(x) + evalInvB(u)

anim = True

if anim :
    %matplotlib notebook
    import matplotlib.pyplot as plt
    from matplotlib.animation import FuncAnimation
    #x = np.random.normal(size=Smat.shape[0])
    xtot = np.zeros(shape = [Smat.shape[0],50])
    xtot[:,0] = x
    fig, ax = plt.subplots(figsize=(2,2))
    ln = plt.imshow(x.reshape(nx), 'BrBG')

def init():
    ln.set_data(x.reshape(nx))
    return ln

def update(frame):
    u = np.random.normal(size=Smat.shape[0])
    xtot[:,frame] = evalInvA(xtot[:,frame-1]) + evalInvB(u)
    ln.set_data(xtot[:,frame].reshape(nx))
    return ln

if anim :
    ani = FuncAnimation(fig, update, frames = 50,
                        init_func=init, blit=False,interval=10)
    p = plt.show()

%matplotlib inline

model = gl.Model.createFromParam(gl.ECov.BESSEL_K,range = scale,param = param,flagRange=False)

rangev = 3.5 * model.getCova(0).getRange()
h = np.linspace(0,rangev,100)
points = [gl.SpacePoint([0.,i]) for i in h]
cov = [model.eval(points[0],i) for i in points]

def correc(kappa,alpha,d=2):
    return gamma(alpha-d/2)/ (gamma(alpha)*(4*np.pi)**(d/2)*kappa**(2*alpha-d))
    

d = 2
param = 3
alpha = param+d/2
hmax = 3.5 * model.getCova(0).getRange()
N = 2**8


ind = np.concatenate([np.arange(int(N/2),N-1,2),np.arange(0,int(N/2),2)])

a= np.pi * (N-1) / hmax
v = np.linspace(-1.,1.,N)
u = a/2 * v
deltau=a/(N-1)

xi = np.meshgrid(u,u)
normxi = np.array([i**2 + j**2 for i in u for j in u]).reshape((len(u),len(u)))

grxi = c * (kappa2 + normxi)**alpha * correc(kappa,param+d/2,d)

nres = 200
time = np.arange(0,nres,1) 
result = np.zeros(shape=[len(ind),nres])
for k in time :
    fourier = np.exp(-k *dt * np.abs(grxi))/(  1./c * (2*np.pi)**d*2*grxi)
    B = np.real(np.fft.fftn(fourier,norm="backward"))
    result[:,k] = B[0,:][ind]
    

X= np.pi * v[np.arange(0,N,2)] /deltau 
points = [gl.SpacePoint([0.,i]) for i in X]
cov = [model.eval(gl.SpacePoint([0.,0]),i)/2. for i in points]
 

plt.plot(X,cov)
for i in range(10):
    covZ=result[:,i * 10]*deltau**d
    plt.plot(X, covZ)


lim = 1.5 * model.getCova(0).getRange()
plt.xlim([-lim,lim])
a = plt.axvline(x=0.)

i = plt.imshow(result)

v

array([-1.        , -0.99215686, -0.98431373, -0.97647059, -0.96862745,
       -0.96078431, -0.95294118, -0.94509804, -0.9372549 , -0.92941176,
       -0.92156863, -0.91372549, -0.90588235, -0.89803922, -0.89019608,
       -0.88235294, -0.8745098 , -0.86666667, -0.85882353, -0.85098039,
       -0.84313725, -0.83529412, -0.82745098, -0.81960784, -0.81176471,
       -0.80392157, -0.79607843, -0.78823529, -0.78039216, -0.77254902,
       -0.76470588, -0.75686275, -0.74901961, -0.74117647, -0.73333333,
       -0.7254902 , -0.71764706, -0.70980392, -0.70196078, -0.69411765,
       -0.68627451, -0.67843137, -0.67058824, -0.6627451 , -0.65490196,
       -0.64705882, -0.63921569, -0.63137255, -0.62352941, -0.61568627,
       -0.60784314, -0.6       , -0.59215686, -0.58431373, -0.57647059,
       -0.56862745, -0.56078431, -0.55294118, -0.54509804, -0.5372549 ,
       -0.52941176, -0.52156863, -0.51372549, -0.50588235, -0.49803922,
       -0.49019608, -0.48235294, -0.4745098 , -0.46666667, -0.45882353,
       -0.45098039, -0.44313725, -0.43529412, -0.42745098, -0.41960784,
       -0.41176471, -0.40392157, -0.39607843, -0.38823529, -0.38039216,
       -0.37254902, -0.36470588, -0.35686275, -0.34901961, -0.34117647,
       -0.33333333, -0.3254902 , -0.31764706, -0.30980392, -0.30196078,
       -0.29411765, -0.28627451, -0.27843137, -0.27058824, -0.2627451 ,
       -0.25490196, -0.24705882, -0.23921569, -0.23137255, -0.22352941,
       -0.21568627, -0.20784314, -0.2       , -0.19215686, -0.18431373,
       -0.17647059, -0.16862745, -0.16078431, -0.15294118, -0.14509804,
       -0.1372549 , -0.12941176, -0.12156863, -0.11372549, -0.10588235,
       -0.09803922, -0.09019608, -0.08235294, -0.0745098 , -0.06666667,
       -0.05882353, -0.05098039, -0.04313725, -0.03529412, -0.02745098,
       -0.01960784, -0.01176471, -0.00392157,  0.00392157,  0.01176471,
        0.01960784,  0.02745098,  0.03529412,  0.04313725,  0.05098039,
        0.05882353,  0.06666667,  0.0745098 ,  0.08235294,  0.09019608,
        0.09803922,  0.10588235,  0.11372549,  0.12156863,  0.12941176,
        0.1372549 ,  0.14509804,  0.15294118,  0.16078431,  0.16862745,
        0.17647059,  0.18431373,  0.19215686,  0.2       ,  0.20784314,
        0.21568627,  0.22352941,  0.23137255,  0.23921569,  0.24705882,
        0.25490196,  0.2627451 ,  0.27058824,  0.27843137,  0.28627451,
        0.29411765,  0.30196078,  0.30980392,  0.31764706,  0.3254902 ,
        0.33333333,  0.34117647,  0.34901961,  0.35686275,  0.36470588,
        0.37254902,  0.38039216,  0.38823529,  0.39607843,  0.40392157,
        0.41176471,  0.41960784,  0.42745098,  0.43529412,  0.44313725,
        0.45098039,  0.45882353,  0.46666667,  0.4745098 ,  0.48235294,
        0.49019608,  0.49803922,  0.50588235,  0.51372549,  0.52156863,
        0.52941176,  0.5372549 ,  0.54509804,  0.55294118,  0.56078431,
        0.56862745,  0.57647059,  0.58431373,  0.59215686,  0.6       ,
        0.60784314,  0.61568627,  0.62352941,  0.63137255,  0.63921569,
        0.64705882,  0.65490196,  0.6627451 ,  0.67058824,  0.67843137,
        0.68627451,  0.69411765,  0.70196078,  0.70980392,  0.71764706,
        0.7254902 ,  0.73333333,  0.74117647,  0.74901961,  0.75686275,
        0.76470588,  0.77254902,  0.78039216,  0.78823529,  0.79607843,
        0.80392157,  0.81176471,  0.81960784,  0.82745098,  0.83529412,
        0.84313725,  0.85098039,  0.85882353,  0.86666667,  0.8745098 ,
        0.88235294,  0.89019608,  0.89803922,  0.90588235,  0.91372549,
        0.92156863,  0.92941176,  0.9372549 ,  0.94509804,  0.95294118,
        0.96078431,  0.96862745,  0.97647059,  0.98431373,  0.99215686,
        1.        ])

result[:,0].max()

357.4611315408981