Enforce units in PiecewiseLinearTransiogram

JuliaEarth · Oct 29, 2024 · a802c58 · a802c58
1 parent 8798694
commit a802c58
Show file tree

Hide file tree

Showing 3 changed files with 12 additions and 2 deletions.
diff --git a/src/theoretical/transiogram/piecewiselinear.jl b/src/theoretical/transiogram/piecewiselinear.jl
@@ -32,7 +32,7 @@ function PiecewiseLinearTransiogram(ball::MetricBall, abscissas::AbstractVector,
   m = size(ordinates, 1)
 
   # abscissa vector
-  x = abscissas
+  x = float(aslen.(abscissas))
 
   # ordinate matrices
   Y = map(1:n) do k

diff --git a/src/utils.jl b/src/utils.jl
@@ -5,6 +5,10 @@
 const Len{T} = Quantity{T,u"𝐋"}
 const InvLen{T} = Quantity{T,u"𝐋^-1"}
 
+aslen(x::Number) = x * u"m"
+aslen(x::Len) = x
+aslen(::Quantity) = throw(ArgumentError("invalid units, please check the documentation"))
+
 asinvlen(x::Number) = x * u"m^-1"
 asinvlen(x::InvLen) = x
 asinvlen(::Quantity) = throw(ArgumentError("invalid units, please check the documentation"))

diff --git a/test/theoretical/transiogram.jl b/test/theoretical/transiogram.jl
@@ -28,5 +28,11 @@
 end
 
 @testset "PiecewiseLinearTransiogram" begin
-  # TODO
+  csv = CSV.File(joinpath(datadir, "facies5.csv"))
+  gtb = georef(csv, ("X", "Y", "Z"))
+  t = EmpiricalTransiogram(gtb, "FACIES", maxlag=20, nlags=20)
+  τ = PiecewiseLinearTransiogram(t.abscissas, t.ordinates)
+  @test τ(0.0u"m") == I(5)
+  @test all(x -> 0 < x < 1, τ(5.0u"m"))
+  @test all(allequal, eachcol(τ(100.0u"m")))
 end