Add SPDEGP

Project.toml

@@ -4,18 +4,22 @@ authors = ["Elias Carvalho <[email protected]> and contributors"]
 version = "0.1.0"
 
 [deps]
+Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 CpuId = "adafc99b-e345-5852-983c-f28acb93d879"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 GeoStatsBase = "323cb8eb-fbf6-51c0-afd0-f8fba70507b2"
 GeoTables = "e502b557-6362-48c1-8219-d30d308dcdb0"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
 CpuId = "0.3"
 GeoStatsBase = "0.39"
 GeoTables = "1.9"
 Meshes = "0.35"
 ProgressMeter = "1.9"
+Tables = "1"
 julia = "1.9"

src/GeoStatsProcesses.jl

@@ -1,15 +1,26 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
 module GeoStatsProcesses
 
-using GeoTables
 using Meshes
-using Random
+using GeoTables
 
+using CpuId
+using Tables
 using ProgressMeter
+
+using Random
 using Distributed
-using CpuId
+using LinearAlgebra
 
-using GeoStatsBase: Ensemble
+using GeoStatsBase: Ensemble, integrate
+using Bessels: gamma
 
 include("interface.jl")
+include("spde.jl")
+
+export SPDEGP
 
 end

src/interface.jl

@@ -1,3 +1,7 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
 """
     GeoStatsProcess
 

src/spde.jl

@@ -0,0 +1,79 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+@kwdef struct SPDEGP <: GeoStatsProcess
+  sill::Float64 = 1.0
+  range::Float64 = 1.0
+end
+
+function randprep(::AbstractRNG, process::SPDEGP, setup::RandSetup)
+  isnothing(setup.geotable) || @error "conditional simulation is not implemented"
+
+  # retrieve sill and range
+  σ = process.sill
+  𝓁 = process.range
+
+  @assert σ > 0 "sill must be positive"
+  @assert 𝓁 > 0 "range must be positive"
+
+  # retrieve problem info
+  𝒟 = setup.domain
+  d = paramdim(𝒟)
+
+  # Beltrami-Laplace discretization
+  B = laplacematrix(𝒟)
+  M = measurematrix(𝒟)
+  Δ = inv(M) * B
+
+  # result of preprocessing
+  pairs = map(setup.varnames) do var
+    # LHS of SPDE (κ² - Δ)Z = τW
+    α = 2one(σ + 𝓁)
+    ν = α - d / 2
+    κ = 1 / 𝓁
+    A = κ^2 * I - Δ
+
+    # covariance structure
+    τ² = σ^2 * κ^(2ν) * (4π)^(d / 2) * gamma(α) / gamma(ν)
+    Q = A'A / τ²
+
+    # factorization
+    F = cholesky(Array(Q))
+    L = inv(Array(F.U))
+
+    # save preprocessed inputs for variable
+    var => (; L)
+  end
+
+  Dict(pairs)
+end
+
+function randsingle(rng::AbstractRNG, ::SPDEGP, setup::RandSetup, prep)
+  # retrieve problem info
+  𝒟 = setup.domain
+  n = nvertices(𝒟)
+
+  # simulation at vertices
+  varreal = map(setup.vartypes, setup.varnames) do V, var
+    # unpack preprocessed parameters
+    L = prep[var].L
+
+    # perform simulation
+    w = randn(rng, V, n)
+    z = L * w
+
+    var => z
+  end
+
+  # vertex table
+  vtable = (; varreal...)
+
+  # change of support
+  vdata = georef(vtable, 𝒟)
+  edata = integrate(vdata, setup.varnames...)
+
+  # columns of element table
+  cols = Tables.columns(values(edata))
+  Dict(var => Tables.getcolumn(cols, var) for var in setup.varnames)
+end