diff --git a/src/GeoStatsFunctions.jl b/src/GeoStatsFunctions.jl index 0061269..ed70633 100644 --- a/src/GeoStatsFunctions.jl +++ b/src/GeoStatsFunctions.jl @@ -97,7 +97,6 @@ export # theoretical transiogram Transiogram, ExponentialTransiogram, - ratematrix, # fitting algorithms VariogramFitAlgo, diff --git a/src/matrices.jl b/src/matrices.jl index e29b1e7..5feb775 100644 --- a/src/matrices.jl +++ b/src/matrices.jl @@ -2,10 +2,6 @@ # Licensed under the MIT License. See LICENSE in the project root. # ------------------------------------------------------------------ -# ---------- -# EMPIRICAL -# ---------- - """ ratematrix(data, var; [parameters]) @@ -102,46 +98,6 @@ function _countmatrix(dom, vals, minlag) C end -# ------------ -# THEORETICAL -# ------------ - -""" - baseratematrix(l, p) - -Transition rate matrix from mean lengths `l` and proportions `p` -that assumes random transitions based on relative proportions. - -The transition rate for the pair `i -> j` is given by `-1 / l[i]` -if `i == j` and by `(p[j] / (1 - p[i])) / l[i]` otherwise. - -## References - -* Carle et al 1998. [Conditional Simulation of Hydrofacies Architecture: - A Transition Probability/Markov Approach](https://doi.org/10.2110/sepmcheg.01.147) -""" -function baseratematrix(l, p) - nₗ = length(l) - nₚ = length(p) - - # sanity checks - if nₗ != nₚ - throw(ArgumentError("lengths and proportions must have the same length")) - end - if !all(pᵢ -> 0 ≤ pᵢ ≤ 1, p) - throw(ArgumentError("proportions must be in interval [0, 1]")) - end - - # Eq. 17 and 18 of Carle et al 1998. - map(Iterators.product(1:nₗ, 1:nₗ)) do (i, j) - if i == j - -1 / l[i] - else - (p[j] / (1 - p[i])) / l[i] - end - end -end - # ----------------- # HELPER FUNCTIONS # ----------------- diff --git a/src/theoretical/transiogram.jl b/src/theoretical/transiogram.jl index 0e7b28e..c46d76d 100644 --- a/src/theoretical/transiogram.jl +++ b/src/theoretical/transiogram.jl @@ -20,12 +20,6 @@ abstract type Transiogram <: GeoStatsFunction end Base.range(t::Transiogram) = maximum(ranges(t)) -# ---------------- -# TRANSIOGRAM API -# ---------------- - -ranges(t::Transiogram) = 1 ./ -diag(ratematrix(t)) - # ---------------- # IMPLEMENTATIONS # ---------------- diff --git a/src/theoretical/transiogram/exponential.jl b/src/theoretical/transiogram/exponential.jl index beecabc..ede9843 100644 --- a/src/theoretical/transiogram/exponential.jl +++ b/src/theoretical/transiogram/exponential.jl @@ -3,15 +3,26 @@ # ------------------------------------------------------------------ """ - ExponentialTransiogram(rate, [ball]) + ExponentialTransiogram(rate; levels=l) + ExponentialTransiogram(ball, rate; levels=l) An exponential transiogram with transition `rate` matrix. -Optionally, specify a metric `ball` to model anisotropy. +Optionally, specify a metric `ball` to model anisotropy, +and the `levels` or categories. + + ExponentialTransiogram(lengths, proportions; levels=l) + ExponentialTransiogram(ball, lengths, proportions; levels=l) + +Alternatively, build transition rate matrix from mean `lengths` +and relative `proportions`. ## References * Carle, S.F. & Fogg, G.E. 1996. [Transition probability-based indicator geostatistics](https://link.springer.com/article/10.1007/BF02083656) + +* Carle et al 1998. [Conditional Simulation of Hydrofacies Architecture: + A Transition Probability/Markov Approach](https://doi.org/10.2110/sepmcheg.01.147) """ struct ExponentialTransiogram{R<:StaticMatrix,L<:AbstractVector,B<:MetricBall} <: Transiogram rate::R @@ -39,13 +50,43 @@ function ExponentialTransiogram(rate::AbstractMatrix; levels=1:size(rate, 1)) ExponentialTransiogram(ball, rate; levels) end -""" - ratematrix(t) +ExponentialTransiogram(ball::MetricBall, lens::AbstractVector, props::AbstractVector; levels=1:length(lens)) = + ExponentialTransiogram(ball, baseratematrix(lens, props); levels) -Return the transition rate matrix of the transiogram `t`. -""" -ratematrix(t::ExponentialTransiogram) = t.rate +ExponentialTransiogram(lens::AbstractVector, props::AbstractVector; levels=1:length(lens)) = + ExponentialTransiogram(baseratematrix(lens, props); levels) + +ranges(t::Transiogram) = 1 ./ -diag(t.rate) levels(t::ExponentialTransiogram) = t.levs (t::ExponentialTransiogram)(h) = exp(h * t.rate) + +# ----------------- +# HELPER FUNCTIONS +# ----------------- + +function baseratematrix(l, p) + nₗ = length(l) + nₚ = length(p) + + # sanity checks + if nₗ != nₚ + throw(ArgumentError("lengths and proportions must have the same length")) + end + if !all(pᵢ -> 0 ≤ pᵢ ≤ 1, p) + throw(ArgumentError("proportions must be in interval [0, 1]")) + end + if !(sum(p) ≈ 1) + throw(ArgumentError("proportions must add up to unit")) + end + + # Eq. 17 and 18 of Carle et al 1998. + map(Iterators.product(1:nₗ, 1:nₗ)) do (i, j) + if i == j + -1 / l[i] + else + (p[j] / (1 - p[i])) / l[i] + end + end +end diff --git a/test/matrices.jl b/test/matrices.jl index f49c13f..8ddc6e3 100644 --- a/test/matrices.jl +++ b/test/matrices.jl @@ -44,13 +44,4 @@ R = GeoStatsFunctions.ratematrix(gtb, "FACIES") @test size(R) == (15, 15) testratematrix(R) - - # base transition rate matrix - R = GeoStatsFunctions.baseratematrix([1.0, 2.0, 3.0]u"m", [0.2, 0.5, 0.3]) - @test R == - [ - -1/1.0 0.5 / (1 - 0.2)/1.0 0.3 / (1 - 0.2)/1.0 - 0.2 / (1 - 0.5)/2.0 -1/2.0 0.3 / (1 - 0.5)/2.0 - 0.2 / (1 - 0.3)/3.0 0.5 / (1 - 0.3)/3.0 -1/3.0 - ] * u"m^-1" end diff --git a/test/transiogram.jl b/test/transiogram.jl index e583bbe..bc6ed5e 100644 --- a/test/transiogram.jl +++ b/test/transiogram.jl @@ -1,9 +1,24 @@ @testset "Transiogram" begin + # base transition rate matrix + R = GeoStatsFunctions.baseratematrix([1.0, 2.0, 3.0]u"m", [0.2, 0.5, 0.3]) + @test R == + [ + -1/1.0 0.5 / (1 - 0.2)/1.0 0.3 / (1 - 0.2)/1.0 + 0.2 / (1 - 0.5)/2.0 -1/2.0 0.3 / (1 - 0.5)/2.0 + 0.2 / (1 - 0.3)/3.0 0.5 / (1 - 0.3)/3.0 -1/3.0 + ] * u"m^-1" + + # corresponding exponential transiogram + t = ExponentialTransiogram([1.0, 2.0, 3.0]u"m", [0.2, 0.5, 0.3]) + @test t isa ExponentialTransiogram + @test GeoStatsFunctions.ranges(t) == [1.0, 2.0, 3.0]u"m" + @test range(t) == 3.0u"m" + + # random transition rate matrix A = rand(3, 3) R = A ./ sum(A, dims=2) t = ExponentialTransiogram(R) @test t isa ExponentialTransiogram - @test ratematrix(t) isa StaticMatrix @test range(t) == maximum(1 ./ -diag(R)) # invalid transition rate matrix