Enforce rate units in MatrixExponentialTransiogram

JuliaEarth · Oct 29, 2024 · a56628a · a56628a
1 parent 0203c72
commit a56628a
Show file tree

Hide file tree

Showing 3 changed files with 7 additions and 2 deletions.
diff --git a/src/theoretical/transiogram/matrixexponential.jl b/src/theoretical/transiogram/matrixexponential.jl
@@ -36,7 +36,7 @@ struct MatrixExponentialTransiogram{R<:StaticMatrix,B<:MetricBall} <: Transiogra
 end
 
 function MatrixExponentialTransiogram(ball::MetricBall, rate::AbstractMatrix)
-  srate = SMatrix{size(rate)...}(rate)
+  srate = SMatrix{size(rate)...}(float(asinvlen.(rate)))
   MatrixExponentialTransiogram{typeof(srate),typeof(ball)}(srate, ball)
 end
 

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

diff --git a/test/theoretical/transiogram.jl b/test/theoretical/transiogram.jl
@@ -19,7 +19,7 @@
   R = A ./ sum(A, dims=2)
   t = MatrixExponentialTransiogram(R)
   @test t isa MatrixExponentialTransiogram
-  @test range(t) == maximum(1 ./ -diag(R))
+  @test range(t) == maximum(1 ./ -diag(R)) * u"m"
 
   # invalid transition rate matrix
   A = rand(3, 2)