Refactor LUSIM and improve docstrings

Author: juliohm

src/simulation/field.jl

@@ -210,7 +210,7 @@ function defaultsimulation(process::GaussianProcess, domain; data=nothing)
   end
 end
 
-defaultsimulation(::IndicatorProcess, domain; data=nothing) = SEQSIM(path=RandomPath())
+defaultsimulation(::IndicatorProcess, domain; data=nothing) = SEQSIM()
 
 # ----------------
 # IMPLEMENTATIONS

src/simulation/field/gaussian/lu.jl

@@ -7,9 +7,6 @@ function preprocess(::AbstractRNG, process::GaussianProcess, method::LUSIM, init
   f = process.func
   μ = process.mean
 
-  # method options
-  factorization = method.factorization
-
   # sanity checks
   isvalid(f) = isstationary(f) && issymmetric(f) && isbanded(f)
   if !isvalid(f)
@@ -59,18 +56,18 @@ function preprocess(::AbstractRNG, process::GaussianProcess, method::LUSIM, init
 
     if isempty(dinds)
       d₂ = zero(eltype(z₁))
-      L₂₂ = factorization(Symmetric(C₂₂)).L
+      L₂₂ = cholesky(Symmetric(C₂₂)).L
     else
       # covariance beween data locations
       C₁₁ = _pairwise(cov, ddom)
       C₁₂ = _pairwise(cov, ddom, sdom)
 
-      L₁₁ = factorization(Symmetric(C₁₁)).L
+      L₁₁ = cholesky(Symmetric(C₁₁)).L
       B₁₂ = L₁₁ \ C₁₂
       A₂₁ = transpose(B₁₂)
 
       d₂ = A₂₁ * (L₁₁ \ z₁)
-      L₂₂ = factorization(Symmetric(C₂₂ - A₂₁ * B₁₂)).L
+      L₂₂ = cholesky(Symmetric(C₂₂ - A₂₁ * B₁₂)).L
     end
 
     (; var, z₁, μ₁, d₂, L₂₂, dinds, sinds)

src/simulation/field/methods.jl

@@ -10,18 +10,14 @@ A method for simulating field processes.
 abstract type FieldSimulationMethod end
 
 """
-    LUSIM(; [options])
+    LUSIM()
 
-The LU simulation method introduced by Alabert 1987.
+The LU simulation method introduced by Alabert 1987 and Davis 1987.
 
 The full covariance matrix is built to include all locations
 of the data, and samples from the multivariate Gaussian are
 drawn via a lower-upper (LU) matrix factorization.
 
-## Options
-
-`factorization` - Factorization function from LinearAlgebra (default to `cholesky`)
-
 ## References
 
 * Alabert 1987. [The practice of fast conditional simulations
@@ -44,9 +40,7 @@ factorize the full covariance.
 For larger domains (e.g. 3D grids), other methods are preferred
 such as [`SEQSIM`](@ref) and [`FFTSIM`](@ref).
 """
-@kwdef struct LUSIM{F} <: FieldSimulationMethod
-  factorization::F = cholesky
-end
+struct LUSIM <: FieldSimulationMethod end
 
 """
     SEQSIM(; [options])
@@ -82,9 +76,9 @@ the nearest neighbors are used without additional constraints.
 
 ### Notes
 
-This method is very sensitive to the neighbor search options.
-Care must be taken to make sure that enough neighbors are used
-in the underlying geostatistical model.
+This method is very sensitive to neighbor search options and
+simulation path. Care must be taken to make sure that enough
+neighbors are used in the underlying geostatistical model.
 """
 @kwdef struct SEQSIM{P,N,D} <: FieldSimulationMethod
   path::P = LinearPath()
@@ -104,6 +98,9 @@ after a fast Fourier transform. White noise is added to the
 phase of the spectrum, and a realization is produced by an
 inverse Fourier transform.
 
+Data conditioning is currently performed with Kriging, which
+accepts the following neighbor search options.
+
 ## Options
 
 * `minneighbors` - Minimum number of neighbors (default to `1`)

test/field.jl

@@ -42,13 +42,6 @@
       @test eltype(real.field1) <: Float64
       @test eltype(real.field2) <: Float64
 
-      # custom factorization
-      rng = StableRNG(123)
-      grid = CartesianGrid(100)
-      proc = GaussianProcess(SphericalCovariance(range=10.0))
-      real = rand(rng, proc, grid; method=LUSIM(factorization=lu), data=data)
-      @test eltype(real.Z) <: Float64
-
       # 2D example
       rng = StableRNG(123)
       grid = CartesianGrid(100, 100)
@@ -61,7 +54,7 @@
       grid = CartesianGrid((100, 100), (0.5, 0.5), (1.0, 1.0))
       data = georef((; Z=[1.0, 0.0, 1.0]), [(25.0, 25.0), (50.0, 75.0), (75.0, 50.0)])
 
-      method = SEQSIM(neighborhood=MetricBall(35.0), maxneighbors=3)
+      method = SEQSIM(maxneighbors=3)
 
       # unconditional simulation
       rng = StableRNG(2017)
@@ -99,13 +92,6 @@
       data = georef((; Z=[1.0, 0.0, 1.0]), [(25.0, 25.0), (50.0, 75.0), (75.0, 50.0)])
       real = rand(rng, proc, grid, method=FFTSIM(maxneighbors=3), data=data)
     end
-
-    proc = GaussianProcess(GaussianVariogram())
-    @test sprint(show, proc) == "GaussianProcess(func: GaussianVariogram(range: 1.0 m, sill: 1.0, nugget: 0.0), mean: 0.0)"
-    @test sprint(show, MIME("text/plain"), proc) == """
-    GaussianProcess
-    ├─ func: GaussianVariogram(range: 1.0 m, sill: 1.0, nugget: 0.0)
-    └─ mean: 0.0"""
   end
 
   @testset "IndicatorProcess" begin
@@ -174,13 +160,6 @@
     rng = StableRNG(2024)
     real = rand(rng, proc, mesh; data)
     @test isapprox(sum(real.Z) / length(real.Z), 0.0, atol=1e-3)
-
-    proc = LindgrenProcess()
-    @test sprint(show, proc) == "LindgrenProcess(range: 1.0 m, sill: 1.0)"
-    @test sprint(show, MIME("text/plain"), proc) == """
-    LindgrenProcess
-    ├─ range: 1.0 m
-    └─ sill: 1.0"""
   end
 
   @testset "QuiltingProcess" begin