Refactor transiogram

src/GeoStatsFunctions.jl

@@ -97,7 +97,6 @@ export
   # theoretical transiogram
   Transiogram,
   ExponentialTransiogram,
-  ratematrix,
 
   # fitting algorithms
   VariogramFitAlgo,

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
 # -----------------

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
 # ----------------

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

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

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