Implement GeoStatsFunctions.scale for all functions

src/theoretical.jl

@@ -58,6 +58,14 @@ Return the maximum effective range of the geostatistical function `f`.
 """
 Base.range(f::GeoStatsFunction) = maximum(radii(metricball(f)))
 
+"""
+    scale(f, s)
+
+Scale the metric ball of the geostatistical function `f`
+with a strictly positive scaling factor `s`.
+"""
+function scale end
+
 """
     nvariates(f)
 

src/theoretical/composite.jl

@@ -21,6 +21,10 @@ struct CompositeFunction{CS,FS} <: GeoStatsFunction
   end
 end
 
+# ---------------------
+# GEOSTATSFUNCTION API
+# ---------------------
+
 isstationary(f::CompositeFunction) = all(isstationary, f.fs)
 
 isisotropic(f::CompositeFunction) = all(isisotropic, f.fs)
@@ -29,17 +33,21 @@ issymmetric(f::CompositeFunction) = all(issymmetric, f.cs) && all(issymmetric, f
 
 isbanded(f::CompositeFunction) = all(isbanded, f.fs)
 
-sill(f::CompositeFunction) = sum(f.cs .* map(sill, f.fs))
-
-nugget(f::CompositeFunction) = sum(f.cs .* map(nugget, f.fs))
-
 metricball(f::CompositeFunction) = metricball(argmax(fᵢ -> range(fᵢ), f.fs))
 
 Base.range(f::CompositeFunction) = maximum(range(fᵢ) for fᵢ in f.fs)
 
+scale(f::CompositeFunction, s::Real) = sum(f.cs .* map(fᵢ -> scale(fᵢ, s), f.fs))
+
 nvariates(f::CompositeFunction) = maximum(size(cᵢ, 1) for cᵢ in f.cs)
 
-scale(f::CompositeFunction, s::Real) = CompositeFunction(f.cs, map(fᵢ -> scale(fᵢ, s), f.fs))
+# -------------------------
+# VARIOGRAM/COVARIANCE API
+# -------------------------
+
+sill(f::CompositeFunction) = sum(f.cs .* map(sill, f.fs))
+
+nugget(f::CompositeFunction) = sum(f.cs .* map(nugget, f.fs))
 
 function structures(f::CompositeFunction)
   ks, fs = f.cs, f.fs
@@ -62,6 +70,10 @@ function structures(f::CompositeFunction)
   ucₒ, ucs, fs
 end
 
+# -----------
+# EVALUATION
+# -----------
+
 (f::CompositeFunction)(h) = sum(f.cs .* map(fᵢ -> fᵢ(h), f.fs))
 (f::CompositeFunction)(u::Point, v::Point) = sum(f.cs .* map(fᵢ -> fᵢ(u, v), f.fs))
 

src/theoretical/covariance.jl

@@ -21,6 +21,8 @@ isbanded(::Type{<:Covariance}) = true
 
 metricball(cov::Covariance) = metricball(cov.γ)
 
+scale(cov::C, s::Real) where {C<:Covariance} = C(scale(cov.γ, s))
+
 nvariates(::Type{<:Covariance}) = 1
 
 (cov::Covariance)(h) = sill(cov.γ) - cov.γ(h)
@@ -65,14 +67,6 @@ function structures(cov::C) where {C<:Covariance}
   cₒ, c, C.(γ)
 end
 
-"""
-    scale(cov, s)
-
-Scale metric ball of covariance `cov` with strictly
-positive scaling factor `s`.
-"""
-scale(cov::C, s::Real) where {C<:Covariance} = C(scale(cov.γ, s))
-
 # -----------
 # IO METHODS
 # -----------

src/theoretical/transiogram.jl

@@ -24,6 +24,11 @@ issymmetric(::Type{<:Transiogram}) = false
 
 isbanded(::Type{<:Transiogram}) = true
 
+function scale(t::Transiogram, s::Real)
+  T = constructor(t)
+  T(s * metricball(t); proportions=proportions(t))
+end
+
 nvariates(t::Transiogram) = length(proportions(t))
 
 # ----------------

src/theoretical/transiogram/matrixexponential.jl

@@ -71,6 +71,13 @@ MatrixExponentialTransiogram(ball::MetricBall; lengths=(1.0, 1.0), proportions=(
 MatrixExponentialTransiogram(; range=1.0, ranges=nothing, rotation=I, lengths=(1.0, 1.0), proportions=(0.5, 0.5)) =
   MatrixExponentialTransiogram(rangeball(range, ranges, rotation), baseratematrix(lengths, proportions))
 
+constructor(::MatrixExponentialTransiogram) = MatrixExponentialTransiogram
+
+function scale(t::MatrixExponentialTransiogram, s::Real)
+  T = constructor(t)
+  T(s * metricball(t), t.rate)
+end
+
 meanlengths(t::MatrixExponentialTransiogram) = Tuple(1 ./ -diag(t.rate))
 
 proportions(t::MatrixExponentialTransiogram) = Tuple(normalize(diag(t(100range(t))), 1))

src/theoretical/transiogram/piecewiselinear.jl

@@ -63,6 +63,13 @@ function PiecewiseLinearTransiogram(abscissas::AbstractVector, ordinates::Abstra
   PiecewiseLinearTransiogram(ball, abscissas, ordinates)
 end
 
+constructor(::PiecewiseLinearTransiogram) = PiecewiseLinearTransiogram
+
+function scale(t::PiecewiseLinearTransiogram, s::Real)
+  T = constructor(t)
+  T(s * metricball(t), t.abscissas, t.ordinates, t.ordinfinity)
+end
+
 proportions(t::PiecewiseLinearTransiogram) = Tuple(diag(t.ordinfinity))
 
 function (t::PiecewiseLinearTransiogram)(h)

src/theoretical/variogram.jl

@@ -17,6 +17,11 @@ issymmetric(::Type{<:Variogram}) = true
 
 isbanded(::Type{<:Variogram}) = false
 
+function scale(γ::Variogram, s::Real)
+  V = constructor(γ)
+  V(s * metricball(γ); sill=sill(γ), nugget=nugget(γ))
+end
+
 nvariates(::Type{<:Variogram}) = 1
 
 """
@@ -60,17 +65,6 @@ function structures(γ::Variogram)
   ucₒ, (uc,), (γ,)
 end
 
-"""
-    scale(γ, s)
-
-Scale metric ball of variogram `γ` with strictly
-positive scaling factor `s`.
-"""
-function scale(γ::Variogram, s::Real)
-  V = constructor(γ)
-  V(s * metricball(γ); sill=sill(γ), nugget=nugget(γ))
-end
-
 # leverage symmetry of variograms
 function pairwise!(Γ, γ::Variogram, domain)
   _, (_, n, k) = matrixparams(γ, domain, domain)

test/theoretical/transiogram.jl

@@ -18,6 +18,12 @@
   # check bandness
   @test all(isbanded, ts)
 
+  # scale metric ball
+  for t in ts
+    t′ = GeoStatsFunctions.scale(t, 2)
+    @test range(t′) == 2.0u"m"
+  end
+
   # number of variates
   @test all(nvariates.(ts) .== 2)
 
@@ -70,6 +76,7 @@ end
   # corresponding exponential transiogram
   t = MatrixExponentialTransiogram(lengths=(1.0u"m", 2.0u"m", 3.0u"m"), proportions=(0.2, 0.5, 0.3))
   @test t isa MatrixExponentialTransiogram
+  @test range(GeoStatsFunctions.scale(t, 2)) == 2u"m"
   @test nvariates(t) == 3
   @test meanlengths(t) == (1.0u"m", 2.0u"m", 3.0u"m")
   @test all(proportions(t) .≈ (0.12403100775193801, 0.38759689922480606, 0.48837209302325607))
@@ -102,6 +109,7 @@ end
   gtb = georef(csv, ("X", "Y", "Z"))
   t = EmpiricalTransiogram(gtb, "FACIES", maxlag=20, nlags=20)
   τ = PiecewiseLinearTransiogram(t.abscissas, t.ordinates)
+  @test range(GeoStatsFunctions.scale(τ, 2)) == 2u"m"
   @test nvariates(τ) == 5
   @test τ(0.0u"m") == I(5)
   @test all(x -> 0 < x < 1, τ(5.0u"m"))

test/theoretical/variogram.jl

@@ -216,7 +216,7 @@
     @test GeoStatsFunctions.constructor(γ)() == γ
   end
 
-  # scale
+  # scale metric ball
   for γ in [GaussianVariogram(), SphericalVariogram(), ExponentialVariogram()]
     g = GeoStatsFunctions.scale(γ, 2)
     @test range(g) == 2.0u"m"