Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update FieldProcess constructors #14

Merged
merged 9 commits into from
Nov 1, 2023
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 9 additions & 4 deletions src/field/gaussian.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,9 @@
# ------------------------------------------------------------------

"""
GaussianProcess([parameters])
GaussianProcess(variogram, mean)
GaussianProcess(variogram; mean)
GaussianProcess(; variogram, mean)
eliascarv marked this conversation as resolved.
Show resolved Hide resolved

Gaussian process with given variogram and mean.

Expand All @@ -13,11 +15,14 @@ Gaussian process with given variogram and mean.
* `mean` - Mean field value (default to `0.0`)

"""
@kwdef struct GaussianProcess{V,T} <: FieldProcess
variogram::V = GaussianVariogram()
mean::T = 0.0
struct GaussianProcess{V,T} <: FieldProcess
variogram::V
mean::T
end

GaussianProcess(variogram; mean=0.0) = GaussianProcess(variogram, mean)
GaussianProcess(; variogram=GaussianVariogram(), mean=0.0) = GaussianProcess(variogram, mean)
eliascarv marked this conversation as resolved.
Show resolved Hide resolved

#---------
# METHODS
#---------
Expand Down
2 changes: 1 addition & 1 deletion src/field/gaussian/fft.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
# ------------------------------------------------------------------

"""
FFTMethod([paramaters])
FFTMethod(; [paramaters])

The FFT Gaussian simulation method introduced by Gutjahr 1997.
The covariance function is perturbed in the frequency domain
Expand Down
2 changes: 1 addition & 1 deletion src/field/gaussian/lu.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
# ------------------------------------------------------------------

"""
LUMethod([paramaters])
LUMethod(; [paramaters])

The LU Gaussian simulation method introduced by Alabert 1987.
The full covariance matrix is built to include all locations
Expand Down
4 changes: 2 additions & 2 deletions src/field/gaussian/seq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
# ------------------------------------------------------------------

"""
SEQMethod([paramaters])
SEQMethod(; [paramaters])

The sequential simulation method introduced by Gomez-Hernandez 1993.
It traverses all locations of the geospatial domain according to a path,
Expand Down Expand Up @@ -49,7 +49,7 @@ function randprep(::AbstractRNG, process::GaussianProcess, method::SEQMethod, se
(; path, minneighbors, maxneighbors, neighborhood, distance, init) = method
probmodel = GeoStatsModels.SimpleKriging(variogram, mean)
marginal = Normal(mean, √sill(variogram))
SequentialProcess(; probmodel, marginal, path, minneighbors, maxneighbors, neighborhood, distance, init)
SequentialProcess(probmodel, marginal; path, minneighbors, maxneighbors, neighborhood, distance, init)
end

function randsingle(rng::AbstractRNG, ::GaussianProcess, ::SEQMethod, setup::RandSetup, prep)
Expand Down
13 changes: 9 additions & 4 deletions src/field/lindgren.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,9 @@
# ------------------------------------------------------------------

"""
LindgrenProcess([paramaters])
LindgrenProcess(range, sill)
LindgrenProcess(range; sill)
LindgrenProcess(; range, sill)

The SPDE process introduced by Lindgren 2011. It relies on a
discretization of the Laplace-Beltrami operator on meshes and
Expand All @@ -20,11 +22,14 @@ is adequate for highly curved domains (e.g. surfaces).
Gaussian Markov random fields: the stochastic partial differential
equation approach](https://rss.onlinelibrary.wiley.com/doi/10.1111/j.1467-9868.2011.00777.x)
"""
@kwdef struct LindgrenProcess <: FieldProcess
range::Float64 = 1.0
sill::Float64 = 1.0
struct LindgrenProcess <: FieldProcess
range::Float64
sill::Float64
end

LindgrenProcess(range; sill=1.0) = LindgrenProcess(range, sill)
LindgrenProcess(; range=1.0, sill=1.0) = LindgrenProcess(range, sill)

function randprep(::AbstractRNG, process::LindgrenProcess, ::DefaultRandMethod, setup::RandSetup)
isnothing(setup.geotable) || @error "conditional simulation is not implemented"

Expand Down
27 changes: 19 additions & 8 deletions src/field/quilting.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
# ------------------------------------------------------------------

"""
QuiltingProcess([paramaters])
QuiltingProcess(trainimg, tilesize; [paramaters])

Image quilting process as described in Hoffimann et al. 2017.

Expand All @@ -28,13 +28,24 @@ Image quilting process as described in Hoffimann et al. 2017.
* Hoffimann et al 2017. *Stochastic simulation by image quilting of process-based geological models.*
* Hoffimann et al 2015. *Geostatistical modeling of evolving landscapes by means of image quilting.*
"""
@kwdef struct QuiltingProcess{TR,TS,O,P,IN,S,T,I} <: FieldProcess
struct QuiltingProcess{TR,TS,O,P,IN,S,T,I} <: FieldProcess
trainimg::TR
tilesize::TS
overlap::O = nothing
path::P = :raster
inactive::IN = nothing
soft::S = nothing
tol::T = 0.1
init::I = NearestInit()
overlap::O
path::P
inactive::IN
soft::S
tol::T
init::I
end

QuiltingProcess(
trainimg,
tilesize;
overlap=nothing,
path=:raster,
inactive=nothing,
soft=nothing,
tol=0.1,
init=NearestInit()
) = QuiltingProcess(trainimg, tilesize, overlap, path, inactive, soft, tol, init)
33 changes: 24 additions & 9 deletions src/field/sequential.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,8 @@
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------


"""
SequentialProcess([paramaters])
SequentialProcess(probmodel, marginal; [paramaters])

A sequential simulation process.

Expand All @@ -15,26 +14,42 @@ and in case there are none, use a `marginal` distribution.

## Parameters

### Required

* `probmodel` - Conditional distribution model from GeoStatsModels.jl
* `marginal` - Marginal distribution from Distributions.jl

### Optional

* `path` - Simulation path (default to `LinearPath()`)
* `minneighbors` - Minimum number of neighbors (default to `1`)
* `maxneighbors` - Maximum number of neighbors (default to `10`)
* `neighborhood` - Search neighborhood (default to `nothing`)
* `distance` - Distance used to find nearest neighbors (default to `Euclidean()`)
* `init` - Data initialization method (default to `NearestInit()`)
"""
@kwdef struct SequentialProcess{M,MD,P,N,D,I} <: FieldProcess
struct SequentialProcess{M,MD,P,N,D,I} <: FieldProcess
probmodel::M
marginal::MD
path::P = LinearPath()
minneighbors::Int = 1
maxneighbors::Int = 10
neighborhood::N = nothing
distance::D = Euclidean()
init::I = NearestInit()
path::P
minneighbors::Int
maxneighbors::Int
neighborhood::N
distance::D
init::I
end

SequentialProcess(
probmodel,
marginal;
path=LinearPath(),
minneighbors=1,
maxneighbors=10,
neighborhood=nothing,
distance=Euclidean(),
init=NearestInit()
) = SequentialProcess(probmodel, marginal, path, minneighbors, maxneighbors, neighborhood, distance, init)

function randprep(::AbstractRNG, process::SequentialProcess, ::DefaultRandMethod, setup::RandSetup)
# retrieve domain info
domain = setup.domain
Expand Down
22 changes: 15 additions & 7 deletions src/field/strata.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,13 +3,18 @@
# ------------------------------------------------------------------

"""
StrataProcess([paramaters])
StrataProcess(environment; [paramaters])

A process to simulate 3D stratigraphic models.

## Parameters

### Required

* `environment` - Geological environment

### Optional

* `state` - Initial geological state
* `stack` - Stacking scheme (:erosional or :depositional)
* `nepochs` - Number of epochs (default to 10)
Expand All @@ -21,11 +26,14 @@ A process to simulate 3D stratigraphic models.
Hoffimann 2018. *Morphodynamic analysis and statistical
synthesis of geormorphic data.*
"""
@kwdef struct StrataProcess{E,S,ST,N,FB,FT} <: FieldProcess
struct StrataProcess{E,S,ST,N,FB,FT} <: FieldProcess
environment::E
state::S = nothing
stack::ST = :erosional
nepochs::N = 10
fillbase::FB = NaN
filltop::FT = NaN
state::S
stack::ST
nepochs::N
fillbase::FB
filltop::FT
end

StrataProcess(environment; state=nothing, stack=:erosional, nepochs=10, fillbase=NaN, filltop=NaN) =
StrataProcess(environment, state, stack, nepochs, fillbase, filltop)
14 changes: 8 additions & 6 deletions src/field/turing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
# ------------------------------------------------------------------

"""
TuringProcess([paramaters])
TuringProcess(; [paramaters])

Turing process as described in Turing 1952.

Expand All @@ -18,9 +18,11 @@ Turing process as described in Turing 1952.

Turing 1952. *The chemical basis of morphogenesis.*
"""
@kwdef struct TuringProcess{P,B,E,I} <: FieldProcess
params::P = nothing
blur::B = nothing
edge::E = nothing
iter::I = 100
struct TuringProcess{P,B,E,I} <: FieldProcess
params::P
blur::B
edge::E
iter::I
end

TuringProcess(; params=nothing, blur=nothing, edge=nothing, iter=100) = TuringProcess(params, blur, edge, iter)
30 changes: 15 additions & 15 deletions test/processes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,22 +22,22 @@
# isotropic simulation
Random.seed!(2019)
dom = CartesianGrid(100, 100)
process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
process = GaussianProcess(GaussianVariogram(range=10.0))
method = FFTMethod()
sims = rand(process, dom, [:z => Float64], 3, method)

# anisotropic simulation
Random.seed!(2019)
dom = CartesianGrid(100, 100)
process = GaussianProcess(variogram=GaussianVariogram(MetricBall((20.0, 5.0))))
process = GaussianProcess(GaussianVariogram(MetricBall((20.0, 5.0))))
method = FFTMethod()
sims = rand(process, dom, [:z => Float64], 3, method)

# simulation on view of grid
Random.seed!(2022)
grid = CartesianGrid(100, 100)
vgrid = view(grid, 1:5000)
process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
process = GaussianProcess(GaussianVariogram(range=10.0))
method = FFTMethod()
sim = rand(process, vgrid, [:z => Float64], method)
@test domain(sim) == vgrid
Expand All @@ -49,7 +49,7 @@
coords = [(25.0, 25.0), (50.0, 75.0), (75.0, 50.0)]
samples = georef(table, coords)
sdomain = CartesianGrid(100, 100)
process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
process = GaussianProcess(GaussianVariogram(range=10.0))
method = FFTMethod(maxneighbors=3)
sim = rand(process, sdomain, samples, method)
end
Expand All @@ -59,7 +59,7 @@
𝒟 = CartesianGrid((100, 100), (0.5, 0.5), (1.0, 1.0))
N = 3

process = GaussianProcess(variogram=SphericalVariogram(range=35.0))
process = GaussianProcess(SphericalVariogram(range=35.0))
method = SEQMethod(neighborhood=MetricBall(30.0), maxneighbors=3)

Random.seed!(2017)
Expand All @@ -82,15 +82,15 @@
# conditional simulation
# ----------------------
rng = MersenneTwister(123)
process = GaussianProcess(variogram=SphericalVariogram(range=10.0))
process = GaussianProcess(SphericalVariogram(range=10.0))
method = LUMethod()
sims = rand(rng, process, 𝒟, 𝒮, 2, method)

# ------------------------
# unconditional simulation
# ------------------------
rng = MersenneTwister(123)
process = GaussianProcess(variogram=SphericalVariogram(range=10.0))
process = GaussianProcess(SphericalVariogram(range=10.0))
method = LUMethod()
sims = rand(rng, process, 𝒟, [:z => Float64], 2, method)

Expand All @@ -99,7 +99,7 @@
# -------------
𝒟 = CartesianGrid(500)
rng = MersenneTwister(123)
process = GaussianProcess(variogram=(SphericalVariogram(range=10.0), GaussianVariogram(range=10.0)))
process = GaussianProcess((SphericalVariogram(range=10.0), GaussianVariogram(range=10.0)))
method = LUMethod(correlation=0.95)
sim = rand(rng, process, 𝒟, [:a => Float64, :b => Float64], method)

Expand All @@ -108,7 +108,7 @@
# -----------
𝒟 = CartesianGrid(100, 100)
rng = MersenneTwister(123)
process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
process = GaussianProcess(GaussianVariogram(range=10.0))
method = LUMethod()
sims = rand(rng, process, 𝒟, [:z => Float64], 3, method)

Expand All @@ -118,7 +118,7 @@
𝒟 = CartesianGrid(100, 100)
rng = MersenneTwister(123)
ball = MetricBall((20.0, 5.0))
process = GaussianProcess(variogram=GaussianVariogram(ball))
process = GaussianProcess(GaussianVariogram(ball))
method = LUMethod()
sims = rand(rng, process, 𝒟, [:z => Float64], 3, method)

Expand All @@ -127,19 +127,19 @@
# ---------------------
𝒟 = CartesianGrid(100)
rng = MersenneTwister(123)
process = GaussianProcess(variogram=SphericalVariogram(range=10.0))
process = GaussianProcess(SphericalVariogram(range=10.0))
method1 = LUMethod(factorization=lu)
method2 = LUMethod(factorization=cholesky)
sim1 = rand(rng, process, 𝒟, 𝒮, 2, method1)
sim2 = rand(rng, process, 𝒟, 𝒮, 2, method2)

# throws
𝒟 = CartesianGrid(100, 100)
process = GaussianProcess(variogram=GaussianVariogram(range=10.0))
process = GaussianProcess(GaussianVariogram(range=10.0))
method = LUMethod()
# only 1 or 2 variables can be simulated simultaneously
@test_throws AssertionError rand(process, 𝒟, [:a => Float64, :b => Float64, :c => Float64], method)
process = GaussianProcess(variogram=(GaussianVariogram(range=10.0),))
process = GaussianProcess((GaussianVariogram(range=10.0),))
# the number of parameters must be equal to the number of variables
@test_throws AssertionError rand(process, 𝒟, [:a => Float64, :b => Float64], method)
end
Expand All @@ -152,7 +152,7 @@
rng = MersenneTwister(2017)
trainimg = geostatsimage("Strebelle")
inactive = [CartesianIndex(i, j) for i in 1:30 for j in 1:30]
process = QuiltingProcess(trainimg=trainimg, tilesize=(30, 30), inactive=inactive)
process = QuiltingProcess(trainimg, (30, 30), inactive=inactive)

sims = rand(rng, process, sdomain, sdata, 3)
@test length(sims) == 3
Expand All @@ -172,7 +172,7 @@
proc = SmoothingProcess()
env = Environment(rng, [proc, proc], [0.5 0.5; 0.5 0.5], ExponentialDuration(rng, 1.0))
sdomain = CartesianGrid(50, 50, 20)
sims = rand(StrataProcess(environment=env), sdomain, [:z => Float64], 3)
sims = rand(StrataProcess(env), sdomain, [:z => Float64], 3)
@test length(sims) == 3
@test size(domain(sims[1])) == (50, 50, 20)
end
Expand Down
Loading