From eb41cb32cadc513d104e1f81576ab6704625eede Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Thu, 23 May 2024 15:55:09 -0300 Subject: [PATCH] Update to Meshes.jl v0.43 --- .github/workflows/FormatPR.yml | 2 +- Project.toml | 4 +-- src/algorithms/ballsearch.jl | 9 ++++-- src/algorithms/fullsearch.jl | 9 ++++-- src/empirical.jl | 5 ++-- src/empirical/partition.jl | 12 ++++---- src/empirical/variogram.jl | 10 +++---- src/empirical/varioplane.jl | 10 +++---- src/fitting.jl | 8 +++--- src/sampling.jl | 4 +-- src/utils.jl | 16 +++++++++++ src/variogram.jl | 6 ++-- src/variogram/circular.jl | 7 +++-- src/variogram/cubic.jl | 8 ++++-- src/variogram/exponential.jl | 5 ++-- src/variogram/gaussian.jl | 5 ++-- src/variogram/matern.jl | 9 +++--- src/variogram/nugget.jl | 4 +-- src/variogram/pentaspherical.jl | 8 ++++-- src/variogram/power.jl | 9 +++--- src/variogram/sinehole.jl | 9 +++--- src/variogram/spherical.jl | 8 ++++-- test/covariance.jl | 36 +++++++++++------------ test/empirical.jl | 16 +++++------ test/fitting.jl | 10 +++---- test/nesting.jl | 14 ++++----- test/variogram.jl | 51 ++++++++++++++++----------------- 27 files changed, 163 insertions(+), 131 deletions(-) diff --git a/.github/workflows/FormatPR.yml b/.github/workflows/FormatPR.yml index f66bbcd..a0a1cd9 100644 --- a/.github/workflows/FormatPR.yml +++ b/.github/workflows/FormatPR.yml @@ -2,7 +2,7 @@ name: FormatPR on: push: branches: - - master + - main jobs: build: runs-on: ubuntu-latest diff --git a/Project.toml b/Project.toml index 2a3dba0..ec47d82 100644 --- a/Project.toml +++ b/Project.toml @@ -29,11 +29,11 @@ GeoStatsFunctionsMakieExt = "Makie" [compat] Bessels = "0.2" Distances = "0.10" -GeoTables = "1.19" +GeoTables = "1.21" InteractiveUtils = "1.9" LinearAlgebra = "1.9" Makie = "0.20" -Meshes = "0.42" +Meshes = "0.43" NearestNeighbors = "0.4" Optim = "1.7" Printf = "1.9" diff --git a/src/algorithms/ballsearch.jl b/src/algorithms/ballsearch.jl index 4529632..9621a03 100644 --- a/src/algorithms/ballsearch.jl +++ b/src/algorithms/ballsearch.jl @@ -3,17 +3,20 @@ # ------------------------------------------------------------------ """ - BallSearchAccum(maxlag, nlags, distance) + BallSearchAccum(nlags, maxlag, distance) Accumulate pairs of points in geospatial data with nearest neighbors inside metric ball. """ -struct BallSearchAccum{T,D} <: VariogramAccumAlgo +struct BallSearchAccum{ℒ<:Len,D} <: VariogramAccumAlgo nlags::Int - maxlag::T + maxlag::ℒ distance::D + BallSearchAccum(nlags, maxlag::ℒ, distance::D) where {ℒ<:Len,D} = new{float(ℒ),D}(nlags, maxlag, distance) end +BallSearchAccum(nlags, maxlag, distance) = BallSearchAccum(nlags, addunit(maxlag, u"m"), distance) + function neighfun(algo::BallSearchAccum, pset) ball = MetricBall(algo.maxlag, algo.distance) searcher = BallSearch(pset, ball) diff --git a/src/algorithms/fullsearch.jl b/src/algorithms/fullsearch.jl index 782ab4c..495414b 100644 --- a/src/algorithms/fullsearch.jl +++ b/src/algorithms/fullsearch.jl @@ -3,17 +3,20 @@ # ------------------------------------------------------------------ """ - FullSearchAccum(maxlag, nlags, distance) + FullSearchAccum(nlags, maxlag, distance) Accumulate pairs of points in geospatial data with exhaustive (or full) search. """ -struct FullSearchAccum{T,D} <: VariogramAccumAlgo +struct FullSearchAccum{ℒ<:Len,D} <: VariogramAccumAlgo nlags::Int - maxlag::T + maxlag::ℒ distance::D + FullSearchAccum(nlags, maxlag::ℒ, distance::D) where {ℒ<:Len,D} = new{float(ℒ),D}(nlags, maxlag, distance) end +FullSearchAccum(nlags, maxlag, distance) = FullSearchAccum(nlags, addunit(maxlag, u"m"), distance) + neighfun(::FullSearchAccum, pset) = j -> (j + 1):nelements(pset) skipfun(::FullSearchAccum) = (i, j) -> false diff --git a/src/empirical.jl b/src/empirical.jl index 1bd166c..1b0f175 100644 --- a/src/empirical.jl +++ b/src/empirical.jl @@ -60,7 +60,8 @@ function accumulate(data, var₁, var₂, estim::VariogramEstimator, algo::Vario V = returntype(estim, z₁, z₂) # lag sums and counts - Σx = zeros(nlags) + ℒ = Meshes.lentype(𝒫) + Σx = zeros(ℒ, nlags) Σy = zeros(V, nlags) ns = zeros(Int, nlags) @@ -78,7 +79,7 @@ function accumulate(data, var₁, var₂, estim::VariogramEstimator, algo::Vario z₂ᵢ = z₂[i] # evaluate geospatial lag - h = evaluate(distance, coordinates(pᵢ), coordinates(pⱼ)) + h = uevaluate(distance, pᵢ, pⱼ) # early exit if out of range exit(h) && continue diff --git a/src/empirical/partition.jl b/src/empirical/partition.jl index 3403f84..4be0898 100644 --- a/src/empirical/partition.jl +++ b/src/empirical/partition.jl @@ -26,29 +26,29 @@ function EmpiricalVariogram(partition::Partition, var₁, var₂=var₁; kwargs. end """ - DirectionalVariogram(direction, data, var₁, var₂=var₁; dtol=1e-6, [parameters]) + DirectionalVariogram(direction, data, var₁, var₂=var₁; dtol=1e-6u"m", [parameters]) Computes the empirical (cross-)variogram for the variables `var₁` and `var₂` stored in -geospatial `data` along a given `direction` with band tolerance `dtol`. +geospatial `data` along a given `direction` with band tolerance `dtol` in length units. Optionally, forward `parameters` for the underlying [`EmpiricalVariogram`](@ref). """ -function DirectionalVariogram(dir, data::AbstractGeoTable, var₁, var₂=var₁; dtol=1e-6, kwargs...) +function DirectionalVariogram(dir, data::AbstractGeoTable, var₁, var₂=var₁; dtol=1e-6u"m", kwargs...) rng = MersenneTwister(123) Π = partition(rng, data, DirectionPartition(dir; tol=dtol)) EmpiricalVariogram(Π, var₁, var₂; kwargs...) end """ - PlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6, [parameters]) + PlanarVariogram(normal, data, var₁, var₂=var₁; ntol=1e-6u"m", [parameters]) Computes the empirical (cross-)variogram for the variables `var₁` and `var₂` stored in geospatial `data` along a plane perpendicular to a `normal` direction with plane -tolerance `ntol`. +tolerance `ntol` in length units. Optionally, forward `parameters` for the underlying [`EmpiricalVariogram`](@ref). """ -function PlanarVariogram(normal, data::AbstractGeoTable, var₁, var₂=var₁; ntol=1e-6, kwargs...) +function PlanarVariogram(normal, data::AbstractGeoTable, var₁, var₂=var₁; ntol=1e-6u"m", kwargs...) rng = MersenneTwister(123) Π = partition(rng, data, PlanePartition(normal; tol=ntol)) EmpiricalVariogram(Π, var₁, var₂; kwargs...) diff --git a/src/empirical/variogram.jl b/src/empirical/variogram.jl index bfb7840..5e43441 100644 --- a/src/empirical/variogram.jl +++ b/src/empirical/variogram.jl @@ -12,7 +12,7 @@ geospatial `data`. ## Parameters * nlags - number of lags (default to `20`) - * maxlag - maximum lag (default to 1/10 of maximum lag of data) + * maxlag - maximum lag in length units (default to 1/10 of maximum lag of data) * distance - custom distance function (default to `Euclidean` distance) * estimator - variogram estimator (default to `:matheron` estimator) * algorithm - accumulation algorithm (default to `:ball`) @@ -44,8 +44,8 @@ See also: [`DirectionalVariogram`](@ref), [`PlanarVariogram`](@ref), * Hoffimann, J and Zadrozny, B. 2019. [Efficient variography with partition variograms] (https://www.sciencedirect.com/science/article/pii/S0098300419302936) """ -struct EmpiricalVariogram{V,D,E} - abscissa::Vector{Float64} +struct EmpiricalVariogram{ℒ<:Len,V,D,E} + abscissa::Vector{ℒ} ordinate::Vector{V} counts::Vector{Int} distance::D @@ -73,7 +73,7 @@ function EmpiricalVariogram( # sanity checks @assert nelem > 1 "variogram requires at least 2 elements" @assert nlags > 0 "number of lags must be positive" - @assert maxlag > 0 "maximum lag distance must be positive" + @assert maxlag > zero(maxlag) "maximum lag distance must be positive" @assert estimator ∈ (:matheron, :cressie) "invalid variogram estimator" @assert algorithm ∈ (:full, :ball) "invalid accumulation algorithm" @@ -82,7 +82,7 @@ function EmpiricalVariogram( # ball search with NearestNeighbors.jl requires AbstractFloat and MinkowskiMetric # https://github.com/KristofferC/NearestNeighbors.jl/issues/13 - isfloat = coordtype(𝒟) <: AbstractFloat + isfloat = Unitful.numtype(Meshes.lentype(𝒟)) <: AbstractFloat isminkowski = distance isa MinkowskiMetric # warn users requesting :ball option with invalid parameters diff --git a/src/empirical/varioplane.jl b/src/empirical/varioplane.jl index bf42837..8645f57 100644 --- a/src/empirical/varioplane.jl +++ b/src/empirical/varioplane.jl @@ -5,14 +5,14 @@ """ EmpiricalVarioplane(data, var₁, var₂=var₁; normal=spheredir(0,0), - nangs=50, ptol=0.5, dtol=0.5, + nangs=50, ptol=0.5u"m", dtol=0.5u"m", [parameters]) Given a `normal` direction, estimate the (cross-)variogram of variables `var₁` and `var₂` along all directions in the corresponding plane of variation. -Optionally, specify the tolerance `ptol` for the plane partition, the tolerance -`dtol` for the direction partition, the number of angles `nangs` in the plane, +Optionally, specify the tolerance `ptol` in length units for the plane partition, the tolerance +`dtol` in length units for the direction partition, the number of angles `nangs` in the plane, and forward the `parameters` to the underlying [`EmpiricalVariogram`](@ref). """ struct EmpiricalVarioplane{T,V} @@ -26,8 +26,8 @@ function EmpiricalVarioplane( var₂=var₁; normal=spheredir(0, 0), nangs=50, - ptol=0.5, - dtol=0.5, + ptol=0.5u"m", + dtol=0.5u"m", kwargs... ) # sanity checks diff --git a/src/fitting.jl b/src/fitting.jl index 18822a8..498a897 100644 --- a/src/fitting.jl +++ b/src/fitting.jl @@ -135,7 +135,7 @@ function fit_impl( # evaluate weights f = algo.weightfun - w = isnothing(f) ? n / sum(n) : map(f, x) + w = isnothing(f) ? n / sum(n) : map(xᵢ -> f(ustrip(xᵢ)), x) # objective function function J(θ) @@ -160,15 +160,15 @@ function fit_impl( rₒ = isnothing(range) ? rmax / 3 : range sₒ = isnothing(sill) ? 0.95 * smax : sill nₒ = isnothing(nugget) ? 1e-6 : nugget - θₒ = [rₒ, sₒ, nₒ] + θₒ = [ustrip(rₒ), sₒ, nₒ] # box constraints δ = 1e-8 rₗ, rᵤ = isnothing(range) ? (0.0, rmax) : (range - δ, range + δ) sₗ, sᵤ = isnothing(sill) ? (0.0, smax) : (sill - δ, sill + δ) nₗ, nᵤ = isnothing(nugget) ? (0.0, nmax) : (nugget - δ, nugget + δ) - l = [rₗ, sₗ, nₗ] - u = [rᵤ, sᵤ, nᵤ] + l = [ustrip(rₗ), sₗ, nₗ] + u = [ustrip(rᵤ), sᵤ, nᵤ] # solve optimization problem sol = Optim.optimize(θ -> J(θ) + λ * L(θ), l, u, θₒ) diff --git a/src/sampling.jl b/src/sampling.jl index 6c7d3ea..8df7ee8 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -40,8 +40,8 @@ end function _spacing(γ::Variogram, g::Geometry) r = range(γ) s = sides(boundingbox(g)) - l = minimum(filter(>(0), s)) - r > 0 ? min(r, l) / 3 : l / 3 + l = minimum(filter(>(zero(eltype(s))), s)) + r > zero(r) ? min(r, l) / 3 : l / 3 end function _dims(γ::Variogram, g::Geometry) diff --git a/src/utils.jl b/src/utils.jl index 9d0b990..994cb32 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -2,6 +2,22 @@ # Licensed under the MIT License. See LICENSE in the project root. # ------------------------------------------------------------------ +const Len{T} = Quantity{T,u"𝐋"} + +addunit(x::Number, u) = x * u +addunit(::Quantity, _) = throw(ArgumentError("invalid units, please check the documentation")) + +unitless(xs::Quantity...) = ustrip.(promote(xs...)) + +function uevaluate(distance, p₁, p₂) + u₁ = unit(Meshes.lentype(p₁)) + u₂ = unit(Meshes.lentype(p₂)) + u = Unitful.promote_unit(u₁, u₂) + v₁ = ustrip.(u, to(p₁)) + v₂ = ustrip.(u, to(p₂)) + evaluate(distance, v₁, v₂) * u +end + """ spheredir(θ, φ) diff --git a/src/variogram.jl b/src/variogram.jl index de09bec..8b5d36f 100644 --- a/src/variogram.jl +++ b/src/variogram.jl @@ -105,9 +105,7 @@ Evaluate the variogram at points `u` and `v`. """ function (γ::Variogram)(u::Point, v::Point) d = metric(γ.ball) - x = coordinates(u) - y = coordinates(v) - h = evaluate(d, x, y) + h = uevaluate(d, u, v) γ(h) end @@ -139,6 +137,8 @@ function (γ::Variogram)(U::Geometry, V::Geometry) mean(γ(u, v) for u in us, v in vs) end +(γ::Variogram)(h) = γ(addunit(h, u"m")) + """ pairwise(γ, domain) diff --git a/src/variogram/circular.jl b/src/variogram/circular.jl index 9b0ab5d..d09bec3 100644 --- a/src/variogram/circular.jl +++ b/src/variogram/circular.jl @@ -19,10 +19,11 @@ variotype(::CircularVariogram) = CircularVariogram isstationary(::Type{<:CircularVariogram}) = true -function (γ::CircularVariogram)(h) +function (γ::CircularVariogram)(h::Len) r = radius(γ.ball) s = γ.sill n = γ.nugget - v = h ≤ r ? 1 - (2 / π) * acos(h / r) + (2h / (π * r)) * sqrt(1 - (h^2 / r^2)) : one(h) - (s - n) * v + (h > zero(h)) * n + h′, r′ = unitless(h, r) + v = h′ ≤ r′ ? 1 - (2 / π) * acos(h′ / r′) + (2h′ / (π * r′)) * sqrt(1 - (h′^2 / r′^2)) : one(h′) + (s - n) * v + (h′ > zero(h′)) * n end diff --git a/src/variogram/cubic.jl b/src/variogram/cubic.jl index 06e8421..1598e79 100644 --- a/src/variogram/cubic.jl +++ b/src/variogram/cubic.jl @@ -23,17 +23,19 @@ variotype(::CubicVariogram) = CubicVariogram isstationary(::Type{<:CubicVariogram}) = true -function (γ::CubicVariogram)(h::T) where {T} +function (γ::CubicVariogram)(h::Len) r = radius(γ.ball) s = γ.sill n = γ.nugget + h′, r′ = unitless(h, r) # constants + T = typeof(h′) c1 = T(35) / T(4) c2 = T(7) / T(2) c3 = T(3) / T(4) - s1 = 7 * (h / r)^2 - c1 * (h / r)^3 + c2 * (h / r)^5 - c3 * (h / r)^7 + s1 = 7 * (h′ / r′)^2 - c1 * (h′ / r′)^3 + c2 * (h′ / r′)^5 - c3 * (h′ / r′)^7 s2 = T(1) - (h < r) * (s - n) * s1 + (h ≥ r) * (s - n) * s2 + (h > 0) * n + (h′ < r′) * (s - n) * s1 + (h′ ≥ r′) * (s - n) * s2 + (h′ > 0) * n end diff --git a/src/variogram/exponential.jl b/src/variogram/exponential.jl index ec57bba..4e99449 100644 --- a/src/variogram/exponential.jl +++ b/src/variogram/exponential.jl @@ -24,9 +24,10 @@ variotype(::ExponentialVariogram) = ExponentialVariogram isstationary(::Type{<:ExponentialVariogram}) = true -function (γ::ExponentialVariogram)(h) +function (γ::ExponentialVariogram)(h::Len) r = radius(γ.ball) s = γ.sill n = γ.nugget - (s - n) * (1 - exp(-3(h / r))) + (h > 0) * n + h′, r′ = unitless(h, r) + (s - n) * (1 - exp(-3(h′ / r′))) + (h′ > 0) * n end diff --git a/src/variogram/gaussian.jl b/src/variogram/gaussian.jl index 2d1c110..41811d0 100644 --- a/src/variogram/gaussian.jl +++ b/src/variogram/gaussian.jl @@ -23,11 +23,12 @@ variotype(::GaussianVariogram) = GaussianVariogram isstationary(::Type{<:GaussianVariogram}) = true -function (γ::GaussianVariogram)(h) +function (γ::GaussianVariogram)(h::Len) # add small eps to nugget # for numerical stability r = radius(γ.ball) s = γ.sill n = γ.nugget + typeof(s)(1e-6) - (s - n) * (1 - exp(-3(h / r)^2)) + (h > 0) * n + h′, r′ = unitless(h, r) + (s - n) * (1 - exp(-3(h′ / r′)^2)) + (h′ > 0) * n end diff --git a/src/variogram/matern.jl b/src/variogram/matern.jl index 7ddf81c..cdeb5b7 100644 --- a/src/variogram/matern.jl +++ b/src/variogram/matern.jl @@ -26,19 +26,20 @@ variotype(::MaternVariogram) = MaternVariogram isstationary(::Type{<:MaternVariogram}) = true -function (γ::MaternVariogram)(h::T) where {T} +function (γ::MaternVariogram)(h::Len) r = radius(γ.ball) s = γ.sill n = γ.nugget ν = γ.order + h′, r′ = unitless(h, r) # shift lag by machine precision to # avoid explosion at the origin - h′ = h + eps(T) + h′ = h′ + eps(typeof(h′)) - δ = √(2ν) * 3(h′ / r) + δ = √(2ν) * 3(h′ / r′) Β = besselk(ν, δ) Γ = gamma(ν) - (s - n) * (1 - 2^(1 - ν) / Γ * δ^ν * Β) + (h > 0) * n + (s - n) * (1 - 2^(1 - ν) / Γ * δ^ν * Β) + (h′ > 0) * n end diff --git a/src/variogram/nugget.jl b/src/variogram/nugget.jl index 02464f5..373e029 100644 --- a/src/variogram/nugget.jl +++ b/src/variogram/nugget.jl @@ -21,10 +21,10 @@ isisotropic(::NuggetEffect) = true sill(γ::NuggetEffect) = γ.nugget -Base.range(::NuggetEffect{T}) where {T} = zero(T) +Base.range(::NuggetEffect) = 0u"m" scale(γ::NuggetEffect, ::Real) = γ -(γ::NuggetEffect)(h) = (h > 0) * γ.nugget +(γ::NuggetEffect)(h::Len) = (h > zero(h)) * γ.nugget (γ::NuggetEffect)(u::Point, v::Point) = ifelse(u == v, zero(γ.nugget), γ.nugget) diff --git a/src/variogram/pentaspherical.jl b/src/variogram/pentaspherical.jl index a423a57..0816dc8 100644 --- a/src/variogram/pentaspherical.jl +++ b/src/variogram/pentaspherical.jl @@ -24,17 +24,19 @@ variotype(::PentasphericalVariogram) = PentasphericalVariogram isstationary(::Type{<:PentasphericalVariogram}) = true -function (γ::PentasphericalVariogram)(h::T) where {T} +function (γ::PentasphericalVariogram)(h::Len) r = radius(γ.ball) s = γ.sill n = γ.nugget + h′, r′ = unitless(h, r) # constants + T = typeof(h′) c1 = T(15) / T(8) c2 = T(5) / T(4) c3 = T(3) / T(8) - s1 = c1 * (h / r) - c2 * (h / r)^3 + c3 * (h / r)^5 + s1 = c1 * (h′ / r′) - c2 * (h′ / r′)^3 + c3 * (h′ / r′)^5 s2 = T(1) - (h < r) * (s - n) * s1 + (h ≥ r) * (s - n) * s2 + (h > 0) * n + (h′ < r′) * (s - n) * s1 + (h′ ≥ r′) * (s - n) * s2 + (h′ > 0) * n end diff --git a/src/variogram/power.jl b/src/variogram/power.jl index 3ff349e..103095c 100644 --- a/src/variogram/power.jl +++ b/src/variogram/power.jl @@ -21,17 +21,16 @@ isstationary(::Type{<:PowerVariogram}) = false isisotropic(::PowerVariogram) = true -function (γ::PowerVariogram)(h) +function (γ::PowerVariogram)(h::Len) s = γ.scaling a = γ.exponent n = γ.nugget - s * h^a + (h > 0) * n + h′ = ustrip(h) + s * h′^a + (h′ > 0) * n end function (γ::PowerVariogram)(u::Point, v::Point) d = Euclidean() - x = coordinates(u) - y = coordinates(v) - h = evaluate(d, x, y) + h = uevaluate(d, u, v) γ(h) end diff --git a/src/variogram/sinehole.jl b/src/variogram/sinehole.jl index 5664abe..e66b0c7 100644 --- a/src/variogram/sinehole.jl +++ b/src/variogram/sinehole.jl @@ -23,15 +23,16 @@ variotype(::SineHoleVariogram) = SineHoleVariogram isstationary(::Type{<:SineHoleVariogram}) = true -function (γ::SineHoleVariogram)(h::T) where {T} +function (γ::SineHoleVariogram)(h::Len) r = radius(γ.ball) s = γ.sill n = γ.nugget + h′, r′ = unitless(h, r) # shift lag by machine precision to # avoid explosion at the origin - h′ = h + eps(T) - c = T(π) + h′ = h′ + eps(typeof(h′)) + c = oftype(h′, π) - (s - n) * (1 - sin(c * h′ / r) / (c * h′ / r)) + (h′ > 0) * n + (s - n) * (1 - sin(c * h′ / r′) / (c * h′ / r′)) + (h′ > 0) * n end diff --git a/src/variogram/spherical.jl b/src/variogram/spherical.jl index b7bbdd4..0b6c16e 100644 --- a/src/variogram/spherical.jl +++ b/src/variogram/spherical.jl @@ -24,16 +24,18 @@ variotype(::SphericalVariogram) = SphericalVariogram isstationary(::Type{<:SphericalVariogram}) = true -function (γ::SphericalVariogram)(h::T) where {T} +function (γ::SphericalVariogram)(h::Len) r = radius(γ.ball) s = γ.sill n = γ.nugget + h′, r′ = unitless(h, r) # constants + T = typeof(h′) c1 = T(3) / T(2) c2 = T(1) / T(2) - s1 = c1 * (h / r) - c2 * (h / r)^3 + s1 = c1 * (h′ / r′) - c2 * (h′ / r′)^3 s2 = T(1) - (h < r) * (s - n) * s1 + (h ≥ r) * (s - n) * s2 + (h > 0) * n + (h′ < r′) * (s - n) * s1 + (h′ ≥ r′) * (s - n) * s2 + (h′ > 0) * n end diff --git a/test/covariance.jl b/test/covariance.jl index 94ae706..f3c6baa 100644 --- a/test/covariance.jl +++ b/test/covariance.jl @@ -1,5 +1,5 @@ @testset "Covariance" begin - x, y = rand(Point3), rand(Point3) + x, y = rand(Point{3}), rand(Point{3}) for (CovType, VarioType) in [ (CircularCovariance, CircularVariogram), (CubicCovariance, CubicVariogram), @@ -46,80 +46,80 @@ @test sill(cov) == 1.0 @test nugget(cov) == 0.0 @test metricball(cov) == MetricBall(1.0) - @test range(cov) == 1.0 + @test range(cov) == 1.0u"m" @test GeoStatsFunctions.scale(cov, 2) == GaussianCovariance(range=2.0, sill=1.0, nugget=0.0) # shows cov = CircularCovariance() - @test sprint(show, cov) == "CircularCovariance(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, cov) == "CircularCovariance(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), cov) == """ CircularCovariance ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" cov = CubicCovariance() - @test sprint(show, cov) == "CubicCovariance(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, cov) == "CubicCovariance(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), cov) == """ CubicCovariance ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" cov = ExponentialCovariance() - @test sprint(show, cov) == "ExponentialCovariance(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, cov) == "ExponentialCovariance(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), cov) == """ ExponentialCovariance ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" cov = GaussianCovariance() - @test sprint(show, cov) == "GaussianCovariance(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, cov) == "GaussianCovariance(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), cov) == """ GaussianCovariance ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" cov = MaternCovariance() - @test sprint(show, cov) == "MaternCovariance(sill: 1.0, nugget: 0.0, order: 1.0, range: 1.0, distance: Euclidean)" + @test sprint(show, cov) == "MaternCovariance(sill: 1.0, nugget: 0.0, order: 1.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), cov) == """ MaternCovariance ├─ sill: 1.0 ├─ nugget: 0.0 ├─ order: 1.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" cov = PentasphericalCovariance() - @test sprint(show, cov) == "PentasphericalCovariance(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, cov) == "PentasphericalCovariance(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), cov) == """ PentasphericalCovariance ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" cov = SineHoleCovariance() - @test sprint(show, cov) == "SineHoleCovariance(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, cov) == "SineHoleCovariance(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), cov) == """ SineHoleCovariance ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" cov = SphericalCovariance() - @test sprint(show, cov) == "SphericalCovariance(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, cov) == "SphericalCovariance(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), cov) == """ SphericalCovariance ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" end diff --git a/test/empirical.jl b/test/empirical.jl index 98cac06..bbd0164 100644 --- a/test/empirical.jl +++ b/test/empirical.jl @@ -4,7 +4,7 @@ sdata = georef((z=ones(3),), Matrix(1.0I, 3, 3)) γ = EmpiricalVariogram(sdata, :z, nlags=2, maxlag=2.0) x, y, n = values(γ) - @test x ≈ [1 / 2, √2] + @test x ≈ [1 / 2, √2] * u"m" @test y[2] == 0.0 @test n == [0, 3] @@ -20,7 +20,7 @@ sdata = georef((z=ones(3),), Matrix(1I, 3, 3)) γ = EmpiricalVariogram(sdata, :z, nlags=2, maxlag=2, algorithm=:full) x, y, n = values(γ) - @test x ≈ [1 / 2, √2] + @test x ≈ [1 / 2, √2] * u"m" @test y[2] == 0.0 @test n == [0, 3] @@ -30,7 +30,7 @@ 𝒟 = georef((z=z,), X) γ = EmpiricalVariogram(𝒟, :z, maxlag=1.0, nlags=5) x, y, n = values(γ) - @test x == [0.1, 0.3, 0.5, 0.7, 0.9] + @test x == [0.1, 0.3, 0.5, 0.7, 0.9] * u"m" @test all(iszero.(n)) # accumulation algorithms give the same result @@ -51,10 +51,10 @@ d = georef((z=rand(rng, 100, 100),)) γ = EmpiricalVariogram(d, :z) @test sprint(show, γ) == - "EmpiricalVariogram(abscissa: [0.353553, ..., 13.8426], ordinate: [0.0, ..., 0.0841137], distance: Euclidean(0.0), estimator: MatheronEstimator(), npairs: 2790126)" + "EmpiricalVariogram(abscissa: [0.353553 m, ..., 13.8426 m], ordinate: [0.0, ..., 0.0841137], distance: Euclidean(0.0), estimator: MatheronEstimator(), npairs: 2790126)" @test sprint(show, MIME"text/plain"(), γ) == """ EmpiricalVariogram - ├─ abscissa: [0.353553, 1.20607, 2.0, ..., 12.2868, 13.1058, 13.8426] + ├─ abscissa: [0.353553 m, 1.20607 m, 2.0 m, ..., 12.2868 m, 13.1058 m, 13.8426 m] ├─ ordinate: [0.0, 0.084454, 0.0849279, ..., 0.0841661, 0.0841251, 0.0841137] ├─ distance: Euclidean(0.0) ├─ estimator: MatheronEstimator() @@ -64,7 +64,7 @@ data = georef((z=rand(Composition{3}, 100),), rand(2, 100)) γ = EmpiricalVariogram(data, :z, maxlag=1.0, algorithm=:full) x, y, n = values(γ) - @test all(≥(0), x) + @test all(≥(0u"m"), x) @test all(≥(0), y) @test all(>(0), n) @@ -72,7 +72,7 @@ data = georef((z=[1 * u"K" for i in 1:100],), rand(2, 100)) γ = EmpiricalVariogram(data, :z, nlags=20) x, y, n = values(γ) - @test all(≥(0), x) + @test all(≥(0u"m"), x) @test y == fill(0.0 * u"K^2", 20) # Matheron's vs Cressie's estimator @@ -91,7 +91,7 @@ data = georef((; Z=img)) γ = EmpiricalVariogram(data, "Z", maxlag=50.0) x, y, n = values(γ) - @test all(≥(0), x) + @test all(≥(0u"m"), x) @test all(>(0.8), y[11:end]) @test all(≥(0), n) end diff --git a/test/fitting.jl b/test/fitting.jl index da7e761..a8d2ba7 100644 --- a/test/fitting.jl +++ b/test/fitting.jl @@ -15,28 +15,28 @@ # fix parameters γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0) - @test isapprox(range(γ), 12.0, atol=1e-3) + @test isapprox(range(γ), 12.0u"m", atol=1e-3u"m") γ = GeoStatsFunctions.fit(GaussianVariogram, g, sill=0.07) @test isapprox(sill(γ), 0.07, atol=1e-3) γ = GeoStatsFunctions.fit(GaussianVariogram, g, nugget=0.05) @test isapprox(nugget(γ), 0.05, atol=1e-3) γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0, sill=0.07) - @test isapprox(range(γ), 12.0, atol=1e-3) + @test isapprox(range(γ), 12.0u"m", atol=1e-3u"m") @test isapprox(sill(γ), 0.07, atol=1e-3) γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0, nugget=0.05) - @test isapprox(range(γ), 12.0, atol=1e-3) + @test isapprox(range(γ), 12.0u"m", atol=1e-3u"m") @test isapprox(nugget(γ), 0.05, atol=1e-3) γ = GeoStatsFunctions.fit(GaussianVariogram, g, sill=0.07, nugget=0.05) @test isapprox(sill(γ), 0.07, atol=1e-3) @test isapprox(nugget(γ), 0.05, atol=1e-3) γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0, sill=0.07, nugget=0.05) - @test isapprox(range(γ), 12.0, atol=1e-3) + @test isapprox(range(γ), 12.0u"m", atol=1e-3u"m") @test isapprox(sill(γ), 0.07, atol=1e-3) @test isapprox(nugget(γ), 0.05, atol=1e-3) # fix maximum parameters γ = GeoStatsFunctions.fit(GaussianVariogram, g, maxrange=5.0) - @test isapprox(range(γ), 5.0, atol=1e-3) + @test isapprox(range(γ), 5.0u"m", atol=1e-3u"m") γ = GeoStatsFunctions.fit(GaussianVariogram, g, maxsill=0.04) @test isapprox(sill(γ), 0.04, atol=1e-3) γ = GeoStatsFunctions.fit(GaussianVariogram, g, maxnugget=0.004) diff --git a/test/nesting.jl b/test/nesting.jl index 8386c89..e2d0922 100644 --- a/test/nesting.jl +++ b/test/nesting.jl @@ -3,7 +3,7 @@ γ = NuggetEffect(0.2) + GaussianVariogram(nugget=0.1, sill=0.8, range=50.0) @test nugget(γ) ≈ 0.3 @test sill(γ) ≈ 1.0 - @test range(γ) ≈ 50.0 + @test range(γ) ≈ 50.0u"m" @test (@elapsed sill(γ)) < 1e-6 @test (@elapsed nugget(γ)) < 1e-6 @test (@allocated sill(γ)) < 32 @@ -11,7 +11,7 @@ γ = 2.0 * NuggetEffect(0.2) @test nugget(γ) ≈ 0.4 @test sill(γ) ≈ 0.4 - @test range(γ) ≈ 0.0 + @test range(γ) ≈ 0.0u"m" @test (@elapsed sill(γ)) < 1e-5 @test (@elapsed nugget(γ)) < 1e-5 @test (@allocated sill(γ)) < 32 @@ -43,15 +43,15 @@ # result type is defined for nested models # see https://github.com/JuliaEarth/GeoStats.jl/issues/121 γ = GaussianVariogram() + ExponentialVariogram() - @test GeoStatsFunctions.returntype(γ, rand(Point3), rand(Point3)) == Float64 + @test GeoStatsFunctions.returntype(γ, Point(0.0, 0.0, 0.0), Point(0.0, 0.0, 0.0)) == Float64 γ = GaussianVariogram(sill=1.0f0, range=1.0f0, nugget=0.1f0) - @test GeoStatsFunctions.returntype(γ, rand(Point3f), rand(Point3f)) == Float32 + @test GeoStatsFunctions.returntype(γ, Point(0f0, 0f0, 0f0), Point(0f0, 0f0, 0f0)) == Float32 # nested model with matrix coefficients C₁ = [1.0 0.5; 0.5 2.0] C₂ = [3.0 0.0; 0.0 3.0] γ = C₁ * GaussianVariogram(range=1.0) + C₂ * SphericalVariogram(range=2.0) - @test range(γ) ≈ 2.0 + @test range(γ) ≈ 2.0u"m" @test sill(γ) ≈ C₁ .+ C₂ @test γ(10.0) ≈ sill(γ) @test γ(Point(10.0, 0.0), Point(0.0, 0.0)) ≈ sill(γ) @@ -60,7 +60,7 @@ # nested model with matrix coefficients C = [1.0 0.0; 0.0 1.0] γ = C * GaussianVariogram() + C * ExponentialVariogram() + C * CubicVariogram() - @test range(γ) ≈ 1.0 + @test range(γ) ≈ 1.0u"m" @test sill(γ) ≈ [3.0 0.0; 0.0 3.0] @test γ(10.0) ≈ sill(γ) @test γ(Point(10.0, 0.0), Point(0.0, 0.0)) ≈ sill(γ) @@ -69,7 +69,7 @@ # test constructor explicitly γ = NestedVariogram((1.0, 2.0), (ExponentialVariogram(), SphericalVariogram())) @test sill(γ) == 3.0 - @test range(γ) == 1.0 + @test range(γ) == 1.0u"m" @test nugget(γ) == 0.0 @test (@elapsed sill(γ)) < 1e-5 @test (@elapsed nugget(γ)) < 1e-5 diff --git a/test/variogram.jl b/test/variogram.jl index c8649ce..bdce539 100644 --- a/test/variogram.jl +++ b/test/variogram.jl @@ -1,7 +1,7 @@ @testset "Variogram" begin rng = StableRNG(123) h = range(0, stop=10, length=50) - x, y = rand(rng, Point3), rand(rng, Point3) + x, y = rand(rng, Point{3}), rand(rng, Point{3}) # stationary variogram models γs = [ @@ -80,7 +80,7 @@ γ = NuggetEffect(nugget=0.2) @test nugget(γ) == 0.2 @test sill(γ) == 0.2 - @test range(γ) == 0.0 + @test range(γ) == 0u"m" # ill-conditioned models and nugget regularization # see https://github.com/JuliaEarth/GeoStats.jl/issues/29 @@ -104,7 +104,7 @@ # sill and nugget in single precision for G in [GaussianVariogram, SphericalVariogram, ExponentialVariogram, MaternVariogram] γ = G(sill=1.0f0) - @test typeof(range(γ)) == Float64 + @test typeof(range(γ)) == typeof(1.0u"m") @test typeof(sill(γ)) == Float32 @test typeof(nugget(γ)) == Float32 end @@ -173,70 +173,69 @@ end # scale - for γ in [ - GaussianVariogram(), - SphericalVariogram(), - ExponentialVariogram() - ] + for γ in [GaussianVariogram(), SphericalVariogram(), ExponentialVariogram()] g = GeoStatsFunctions.scale(γ, 2) - @test range(g) == 2 + @test range(g) == 2.0u"m" end # scale with NestedVariogram γ = GaussianVariogram(range=2.0) + ExponentialVariogram(range=3.0) g = GeoStatsFunctions.scale(γ, 2) - @test range(g) == 6 + @test range(g) == 6.0u"m" # scale doesn't affect NuggetEffect γ = NuggetEffect() g = GeoStatsFunctions.scale(γ, 2) @test g == γ + # error: h must have length units + @test_throws ArgumentError γ(0.0u"s") + # shows γ = CircularVariogram() - @test sprint(show, γ) == "CircularVariogram(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, γ) == "CircularVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), γ) == """ CircularVariogram ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" γ = CubicVariogram() - @test sprint(show, γ) == "CubicVariogram(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, γ) == "CubicVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), γ) == """ CubicVariogram ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" γ = ExponentialVariogram() - @test sprint(show, γ) == "ExponentialVariogram(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, γ) == "ExponentialVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), γ) == """ ExponentialVariogram ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" γ = GaussianVariogram() - @test sprint(show, γ) == "GaussianVariogram(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, γ) == "GaussianVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), γ) == """ GaussianVariogram ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" γ = MaternVariogram() - @test sprint(show, γ) == "MaternVariogram(sill: 1.0, nugget: 0.0, order: 1.0, range: 1.0, distance: Euclidean)" + @test sprint(show, γ) == "MaternVariogram(sill: 1.0, nugget: 0.0, order: 1.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), γ) == """ MaternVariogram ├─ sill: 1.0 ├─ nugget: 0.0 ├─ order: 1.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" γ = NuggetEffect() @@ -246,12 +245,12 @@ └─ nugget: 1.0""" γ = PentasphericalVariogram() - @test sprint(show, γ) == "PentasphericalVariogram(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, γ) == "PentasphericalVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), γ) == """ PentasphericalVariogram ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" γ = PowerVariogram() @@ -263,20 +262,20 @@ └─ exponent: 1.0""" γ = SineHoleVariogram() - @test sprint(show, γ) == "SineHoleVariogram(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, γ) == "SineHoleVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), γ) == """ SineHoleVariogram ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" γ = SphericalVariogram() - @test sprint(show, γ) == "SphericalVariogram(sill: 1.0, nugget: 0.0, range: 1.0, distance: Euclidean)" + @test sprint(show, γ) == "SphericalVariogram(sill: 1.0, nugget: 0.0, range: 1.0 m, distance: Euclidean)" @test sprint(show, MIME"text/plain"(), γ) == """ SphericalVariogram ├─ sill: 1.0 ├─ nugget: 0.0 - ├─ range: 1.0 + ├─ range: 1.0 m └─ distance: Euclidean""" end