-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* Add SPDEGP * Update compat
- Loading branch information
Showing
4 changed files
with
103 additions
and
4 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -4,18 +4,23 @@ 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] | ||
Bessels = "0.2" | ||
CpuId = "0.3" | ||
GeoStatsBase = "0.39" | ||
GeoTables = "1.9" | ||
Meshes = "0.35" | ||
ProgressMeter = "1.9" | ||
Tables = "1.11" | ||
julia = "1.9" |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |