From eb41cb32cadc513d104e1f81576ab6704625eede Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Thu, 23 May 2024 15:55:09 -0300 Subject: [PATCH 1/8] 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 From 81d081ee5e76ac342f53844caf412575118280ae Mon Sep 17 00:00:00 2001 From: Elias Carvalho <73039601+eliascarv@users.noreply.github.com> Date: Thu, 23 May 2024 17:51:39 -0300 Subject: [PATCH 2/8] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Júlio Hoffimann --- src/variogram/matern.jl | 2 +- src/variogram/sinehole.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/variogram/matern.jl b/src/variogram/matern.jl index cdeb5b7..adac2f0 100644 --- a/src/variogram/matern.jl +++ b/src/variogram/matern.jl @@ -35,7 +35,7 @@ function (γ::MaternVariogram)(h::Len) # shift lag by machine precision to # avoid explosion at the origin - h′ = h′ + eps(typeof(h′)) + h′ += eps(typeof(h′)) δ = √(2ν) * 3(h′ / r′) Β = besselk(ν, δ) diff --git a/src/variogram/sinehole.jl b/src/variogram/sinehole.jl index e66b0c7..8763066 100644 --- a/src/variogram/sinehole.jl +++ b/src/variogram/sinehole.jl @@ -31,7 +31,7 @@ function (γ::SineHoleVariogram)(h::Len) # shift lag by machine precision to # avoid explosion at the origin - h′ = h′ + eps(typeof(h′)) + h′ += eps(typeof(h′)) c = oftype(h′, π) (s - n) * (1 - sin(c * h′ / r′) / (c * h′ / r′)) + (h′ > 0) * n From 45ac9715ea1e352c9057a6a4be82567ec62ea904 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Thu, 23 May 2024 18:47:15 -0300 Subject: [PATCH 3/8] Apply suggestions --- src/fitting.jl | 19 ++++++++++--------- src/sampling.jl | 2 +- src/variogram.jl | 2 -- test/empirical.jl | 16 ++++++++++++++-- test/fitting.jl | 10 +++++----- test/nesting.jl | 4 ++-- test/variogram.jl | 15 ++++++--------- 7 files changed, 38 insertions(+), 30 deletions(-) diff --git a/src/fitting.jl b/src/fitting.jl index 498a897..6cb304f 100644 --- a/src/fitting.jl +++ b/src/fitting.jl @@ -130,12 +130,13 @@ function fit_impl( n = n[n .> 0] # strip units if necessary - 𝓊 = unit(first(y)) + ux = unit(first(x)) + uy = unit(first(y)) y = ustrip.(y) # evaluate weights f = algo.weightfun - w = isnothing(f) ? n / sum(n) : map(xᵢ -> f(ustrip(xᵢ)), x) + w = isnothing(f) ? n / sum(n) : map(xᵢ -> ustrip(f(xᵢ)), x) # objective function function J(θ) @@ -160,15 +161,15 @@ function fit_impl( rₒ = isnothing(range) ? rmax / 3 : range sₒ = isnothing(sill) ? 0.95 * smax : sill nₒ = isnothing(nugget) ? 1e-6 : nugget - θₒ = [ustrip(rₒ), sₒ, nₒ] + θₒ = [ustrip(ux, 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 = [ustrip(rₗ), sₗ, nₗ] - u = [ustrip(rᵤ), sᵤ, nᵤ] + rₗ, rᵤ = isnothing(range) ? (zero(rmax), rmax) : (range - δ * unit(range), range + δ * unit(range)) + sₗ, sᵤ = isnothing(sill) ? (zero(smax), smax) : (sill - δ, sill + δ) + nₗ, nᵤ = isnothing(nugget) ? (zero(nmax), nmax) : (nugget - δ, nugget + δ) + l = [ustrip(ux, rₗ), sₗ, nₗ] + u = [ustrip(ux, rᵤ), sᵤ, nᵤ] # solve optimization problem sol = Optim.optimize(θ -> J(θ) + λ * L(θ), l, u, θₒ) @@ -176,7 +177,7 @@ function fit_impl( θ = Optim.minimizer(sol) # optimal variogram (with units) - γ = V(ball(θ[1]), sill=θ[2] * 𝓊, nugget=θ[3] * 𝓊) + γ = V(ball(θ[1]), sill=θ[2] * uy, nugget=θ[3] * uy) γ, ϵ end diff --git a/src/sampling.jl b/src/sampling.jl index 8df7ee8..b4e42ac 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -40,7 +40,7 @@ end function _spacing(γ::Variogram, g::Geometry) r = range(γ) s = sides(boundingbox(g)) - l = minimum(filter(>(zero(eltype(s))), s)) + l = minimum(filter(sᵢ -> sᵢ > zero(sᵢ), s)) r > zero(r) ? min(r, l) / 3 : l / 3 end diff --git a/src/variogram.jl b/src/variogram.jl index 8b5d36f..ce3501c 100644 --- a/src/variogram.jl +++ b/src/variogram.jl @@ -137,8 +137,6 @@ 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/test/empirical.jl b/test/empirical.jl index bbd0164..fdd47f8 100644 --- a/test/empirical.jl +++ b/test/empirical.jl @@ -101,8 +101,20 @@ data = georef((z=img,)) γ = EmpiricalVarioplane(data, :z, maxlag=50.0) @test sprint(show, γ) == "EmpiricalVarioplane" - @test sprint(show, MIME"text/plain"(), γ) == - "EmpiricalVarioplane\n N° pairs\n └─0.00° → 372500\n └─3.67° → 304782\n └─7.35° → 298306\n └─11.02° → 297432\n └─14.69° → 297243\n ⋮\n └─165.31° → 293643\n └─168.98° → 295850\n └─172.65° → 296931\n └─176.33° → 306528\n └─180.00° → 372500" + @test sprint(show, MIME"text/plain"(), γ) == """ + EmpiricalVarioplane + N° pairs + └─0.00° → 372500 + └─3.67° → 304782 + └─7.35° → 298306 + └─11.02° → 297432 + └─14.69° → 297243 + ⋮ + └─165.31° → 293643 + └─168.98° → 295850 + └─172.65° → 296931 + └─176.33° → 306528 + └─180.00° → 372500""" end @testset "Directional" begin diff --git a/test/fitting.jl b/test/fitting.jl index a8d2ba7..878b9ed 100644 --- a/test/fitting.jl +++ b/test/fitting.jl @@ -14,28 +14,28 @@ @test isapprox(sill(γ₄), 0.054, atol=1e-3) # fix parameters - γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0) + γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0u"m") @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) + γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0u"m", sill=0.07) @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) + γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0u"m", nugget=0.05) @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) + γ = GeoStatsFunctions.fit(GaussianVariogram, g, range=12.0u"m", sill=0.07, nugget=0.05) @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) + γ = GeoStatsFunctions.fit(GaussianVariogram, g, maxrange=5.0u"m") @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) diff --git a/test/nesting.jl b/test/nesting.jl index e2d0922..8a59b4c 100644 --- a/test/nesting.jl +++ b/test/nesting.jl @@ -53,7 +53,7 @@ γ = C₁ * GaussianVariogram(range=1.0) + C₂ * SphericalVariogram(range=2.0) @test range(γ) ≈ 2.0u"m" @test sill(γ) ≈ C₁ .+ C₂ - @test γ(10.0) ≈ sill(γ) + @test γ(10.0u"m") ≈ sill(γ) @test γ(Point(10.0, 0.0), Point(0.0, 0.0)) ≈ sill(γ) @test isstationary(γ) @@ -62,7 +62,7 @@ γ = C * GaussianVariogram() + C * ExponentialVariogram() + C * CubicVariogram() @test range(γ) ≈ 1.0u"m" @test sill(γ) ≈ [3.0 0.0; 0.0 3.0] - @test γ(10.0) ≈ sill(γ) + @test γ(10.0u"m") ≈ sill(γ) @test γ(Point(10.0, 0.0), Point(0.0, 0.0)) ≈ sill(γ) @test isstationary(γ) diff --git a/test/variogram.jl b/test/variogram.jl index bdce539..a907aa9 100644 --- a/test/variogram.jl +++ b/test/variogram.jl @@ -1,6 +1,6 @@ @testset "Variogram" begin rng = StableRNG(123) - h = range(0, stop=10, length=50) + h = range(0, stop=10, length=50) * u"m" x, y = rand(rng, Point{3}), rand(rng, Point{3}) # stationary variogram models @@ -61,12 +61,12 @@ # some variograms are non-decreasing for γ in (γnd ∪ [sum(γnd)]) - @test all(γ.(h) .≤ γ.(h .+ 1)) + @test all(γ.(h) .≤ γ.(h .+ 1u"m")) end # variograms are valid at the origin for γ in (γs ∪ γn ∪ γnd) - @test !isnan(γ(0.0)) && !isinf(γ(0.0)) + @test !isnan(γ(0.0u"m")) && !isinf(γ(0.0u"m")) end # practical ranges @@ -97,8 +97,8 @@ # nugget regularization γ = GaussianVariogram(range=20.0, nugget=0.1) C = sill(γ) .- GeoStatsFunctions.pairwise(γ, pset) - @test γ(0) == 0 - @test γ(1e-6) > 0 + @test γ(0u"m") == 0 + @test γ(1e-6u"m") > 0 @test cond(C) < 100.0 # sill and nugget in single precision @@ -125,7 +125,7 @@ # unitful non-stationary types γn = [PowerVariogram(scaling=1.0u"K^2")] for γ in γs - @test unit(γ(1.0)) == u"K^2" + @test unit(γ(1.0u"m")) == u"K^2" end 𝒟 = PointSet(Matrix(1.0I, 3, 3)) @@ -188,9 +188,6 @@ 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 m, distance: Euclidean)" From cb1c74fe6b35ffa6506fa14e84944bf7fa6d6806 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Thu, 23 May 2024 19:05:11 -0300 Subject: [PATCH 4/8] Apply suggestions --- src/fitting.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/fitting.jl b/src/fitting.jl index 6cb304f..28c1ad0 100644 --- a/src/fitting.jl +++ b/src/fitting.jl @@ -134,6 +134,10 @@ function fit_impl( uy = unit(first(y)) y = ustrip.(y) + # strip units of `range` and `maxrange` + range′ = isnothing(range) ? range : ustrip(ux, range) + maxrange′ = isnothing(maxrange) ? maxrange : ustrip(ux, maxrange) + # evaluate weights f = algo.weightfun w = isnothing(f) ? n / sum(n) : map(xᵢ -> ustrip(f(xᵢ)), x) @@ -151,25 +155,25 @@ function fit_impl( λ = sum(yᵢ -> yᵢ^2, y) # maximum range, sill and nugget - xmax = maximum(x) + xmax = ustrip(maximum(x)) ymax = maximum(y) - rmax = isnothing(maxrange) ? xmax : maxrange + rmax = isnothing(maxrange′) ? xmax : maxrange′ smax = isnothing(maxsill) ? ymax : maxsill nmax = isnothing(maxnugget) ? ymax : maxnugget # initial guess - rₒ = isnothing(range) ? rmax / 3 : range + rₒ = isnothing(range′) ? rmax / 3 : range′ sₒ = isnothing(sill) ? 0.95 * smax : sill nₒ = isnothing(nugget) ? 1e-6 : nugget - θₒ = [ustrip(ux, rₒ), sₒ, nₒ] + θₒ = [rₒ, sₒ, nₒ] # box constraints δ = 1e-8 - rₗ, rᵤ = isnothing(range) ? (zero(rmax), rmax) : (range - δ * unit(range), range + δ * unit(range)) + rₗ, rᵤ = isnothing(range′) ? (zero(rmax), rmax) : (range′ - δ, range′ + δ) sₗ, sᵤ = isnothing(sill) ? (zero(smax), smax) : (sill - δ, sill + δ) nₗ, nᵤ = isnothing(nugget) ? (zero(nmax), nmax) : (nugget - δ, nugget + δ) - l = [ustrip(ux, rₗ), sₗ, nₗ] - u = [ustrip(ux, rᵤ), sᵤ, nᵤ] + l = [rₗ, sₗ, nₗ] + u = [rᵤ, sᵤ, nᵤ] # solve optimization problem sol = Optim.optimize(θ -> J(θ) + λ * L(θ), l, u, θₒ) From 602335ac01319df1bfef92c731c1586a361180f6 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Fri, 24 May 2024 08:50:55 -0300 Subject: [PATCH 5/8] Apply suggestions --- src/fitting.jl | 32 ++++++++++++++++++-------------- src/utils.jl | 6 +++++- src/variogram/circular.jl | 2 +- src/variogram/cubic.jl | 2 +- src/variogram/exponential.jl | 2 +- src/variogram/gaussian.jl | 2 +- src/variogram/matern.jl | 2 +- src/variogram/nugget.jl | 2 +- src/variogram/pentaspherical.jl | 2 +- src/variogram/power.jl | 2 +- src/variogram/sinehole.jl | 2 +- src/variogram/spherical.jl | 2 +- test/nesting.jl | 4 ++-- test/variogram.jl | 12 ++++++------ 14 files changed, 41 insertions(+), 33 deletions(-) diff --git a/src/fitting.jl b/src/fitting.jl index 28c1ad0..40f7039 100644 --- a/src/fitting.jl +++ b/src/fitting.jl @@ -130,13 +130,17 @@ function fit_impl( n = n[n .> 0] # strip units if necessary - ux = unit(first(x)) - uy = unit(first(y)) - y = ustrip.(y) + ux = unit(eltype(x)) + uy = unit(eltype(y)) + y′ = ustrip.(y) - # strip units of `range` and `maxrange` + # strip units of kwargs range′ = isnothing(range) ? range : ustrip(ux, range) - maxrange′ = isnothing(maxrange) ? maxrange : ustrip(ux, maxrange) + sill′ = isnothing(sill) ? sill : ustrip(uy, sill) + nugget′ = isnothing(nugget) ? nugget : ustrip(uy, nugget) + maxrange′ = isnothing(maxrange) ? maxrange : ustrip(ux, maxrange) + maxsill′ = isnothing(maxsill) ? maxsill : ustrip(uy, maxsill) + maxnugget′ = isnothing(maxnugget) ? maxnugget : ustrip(uy, maxnugget) # evaluate weights f = algo.weightfun @@ -145,33 +149,33 @@ function fit_impl( # objective function function J(θ) γ = V(ball(θ[1]), sill=θ[2], nugget=θ[3]) - sum(i -> w[i] * (γ(x[i]) - y[i])^2, eachindex(w, x, y)) + sum(i -> w[i] * (γ(x[i]) - y′[i])^2, eachindex(w, x, y′)) end # linear constraint (sill ≥ nugget) L(θ) = θ[2] ≥ θ[3] ? 0.0 : θ[3] - θ[2] # penalty for linear constraint (J + λL) - λ = sum(yᵢ -> yᵢ^2, y) + λ = sum(yᵢ -> yᵢ^2, y′) # maximum range, sill and nugget xmax = ustrip(maximum(x)) - ymax = maximum(y) + ymax = maximum(y′) rmax = isnothing(maxrange′) ? xmax : maxrange′ - smax = isnothing(maxsill) ? ymax : maxsill - nmax = isnothing(maxnugget) ? ymax : maxnugget + smax = isnothing(maxsill′) ? ymax : maxsill′ + nmax = isnothing(maxnugget′) ? ymax : maxnugget′ # initial guess rₒ = isnothing(range′) ? rmax / 3 : range′ - sₒ = isnothing(sill) ? 0.95 * smax : sill - nₒ = isnothing(nugget) ? 1e-6 : nugget + sₒ = isnothing(sill′) ? 0.95 * smax : sill′ + nₒ = isnothing(nugget′) ? 1e-6 : nugget′ θₒ = [rₒ, sₒ, nₒ] # box constraints δ = 1e-8 rₗ, rᵤ = isnothing(range′) ? (zero(rmax), rmax) : (range′ - δ, range′ + δ) - sₗ, sᵤ = isnothing(sill) ? (zero(smax), smax) : (sill - δ, sill + δ) - nₗ, nᵤ = isnothing(nugget) ? (zero(nmax), nmax) : (nugget - δ, nugget + δ) + sₗ, sᵤ = isnothing(sill′) ? (zero(smax), smax) : (sill′ - δ, sill′ + δ) + nₗ, nᵤ = isnothing(nugget′) ? (zero(nmax), nmax) : (nugget′ - δ, nugget′ + δ) l = [rₗ, sₗ, nₗ] u = [rᵤ, sᵤ, nᵤ] diff --git a/src/utils.jl b/src/utils.jl index 994cb32..625c703 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -7,7 +7,11 @@ 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...)) +unitless(a::Number, b::Quantity) = a, ustrip(b) +function unitless(a::Quantity, b::Quantity) + u = Unitful.promote_unit(unit(a), unit(b)) + ustrip(u, a), ustrip(u, b) +end function uevaluate(distance, p₁, p₂) u₁ = unit(Meshes.lentype(p₁)) diff --git a/src/variogram/circular.jl b/src/variogram/circular.jl index d09bec3..e884f99 100644 --- a/src/variogram/circular.jl +++ b/src/variogram/circular.jl @@ -19,7 +19,7 @@ variotype(::CircularVariogram) = CircularVariogram isstationary(::Type{<:CircularVariogram}) = true -function (γ::CircularVariogram)(h::Len) +function (γ::CircularVariogram)(h) r = radius(γ.ball) s = γ.sill n = γ.nugget diff --git a/src/variogram/cubic.jl b/src/variogram/cubic.jl index 1598e79..d7cba8e 100644 --- a/src/variogram/cubic.jl +++ b/src/variogram/cubic.jl @@ -23,7 +23,7 @@ variotype(::CubicVariogram) = CubicVariogram isstationary(::Type{<:CubicVariogram}) = true -function (γ::CubicVariogram)(h::Len) +function (γ::CubicVariogram)(h) r = radius(γ.ball) s = γ.sill n = γ.nugget diff --git a/src/variogram/exponential.jl b/src/variogram/exponential.jl index 4e99449..4623a3b 100644 --- a/src/variogram/exponential.jl +++ b/src/variogram/exponential.jl @@ -24,7 +24,7 @@ variotype(::ExponentialVariogram) = ExponentialVariogram isstationary(::Type{<:ExponentialVariogram}) = true -function (γ::ExponentialVariogram)(h::Len) +function (γ::ExponentialVariogram)(h) r = radius(γ.ball) s = γ.sill n = γ.nugget diff --git a/src/variogram/gaussian.jl b/src/variogram/gaussian.jl index 41811d0..7294f71 100644 --- a/src/variogram/gaussian.jl +++ b/src/variogram/gaussian.jl @@ -23,7 +23,7 @@ variotype(::GaussianVariogram) = GaussianVariogram isstationary(::Type{<:GaussianVariogram}) = true -function (γ::GaussianVariogram)(h::Len) +function (γ::GaussianVariogram)(h) # add small eps to nugget # for numerical stability r = radius(γ.ball) diff --git a/src/variogram/matern.jl b/src/variogram/matern.jl index adac2f0..5e5117f 100644 --- a/src/variogram/matern.jl +++ b/src/variogram/matern.jl @@ -26,7 +26,7 @@ variotype(::MaternVariogram) = MaternVariogram isstationary(::Type{<:MaternVariogram}) = true -function (γ::MaternVariogram)(h::Len) +function (γ::MaternVariogram)(h) r = radius(γ.ball) s = γ.sill n = γ.nugget diff --git a/src/variogram/nugget.jl b/src/variogram/nugget.jl index 373e029..39bfe84 100644 --- a/src/variogram/nugget.jl +++ b/src/variogram/nugget.jl @@ -25,6 +25,6 @@ Base.range(::NuggetEffect) = 0u"m" scale(γ::NuggetEffect, ::Real) = γ -(γ::NuggetEffect)(h::Len) = (h > zero(h)) * γ.nugget +(γ::NuggetEffect)(h) = (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 0816dc8..f8958fb 100644 --- a/src/variogram/pentaspherical.jl +++ b/src/variogram/pentaspherical.jl @@ -24,7 +24,7 @@ variotype(::PentasphericalVariogram) = PentasphericalVariogram isstationary(::Type{<:PentasphericalVariogram}) = true -function (γ::PentasphericalVariogram)(h::Len) +function (γ::PentasphericalVariogram)(h) r = radius(γ.ball) s = γ.sill n = γ.nugget diff --git a/src/variogram/power.jl b/src/variogram/power.jl index 103095c..d0651b0 100644 --- a/src/variogram/power.jl +++ b/src/variogram/power.jl @@ -21,7 +21,7 @@ isstationary(::Type{<:PowerVariogram}) = false isisotropic(::PowerVariogram) = true -function (γ::PowerVariogram)(h::Len) +function (γ::PowerVariogram)(h) s = γ.scaling a = γ.exponent n = γ.nugget diff --git a/src/variogram/sinehole.jl b/src/variogram/sinehole.jl index 8763066..c5cd6ee 100644 --- a/src/variogram/sinehole.jl +++ b/src/variogram/sinehole.jl @@ -23,7 +23,7 @@ variotype(::SineHoleVariogram) = SineHoleVariogram isstationary(::Type{<:SineHoleVariogram}) = true -function (γ::SineHoleVariogram)(h::Len) +function (γ::SineHoleVariogram)(h) r = radius(γ.ball) s = γ.sill n = γ.nugget diff --git a/src/variogram/spherical.jl b/src/variogram/spherical.jl index 0b6c16e..98d8cb6 100644 --- a/src/variogram/spherical.jl +++ b/src/variogram/spherical.jl @@ -24,7 +24,7 @@ variotype(::SphericalVariogram) = SphericalVariogram isstationary(::Type{<:SphericalVariogram}) = true -function (γ::SphericalVariogram)(h::Len) +function (γ::SphericalVariogram)(h) r = radius(γ.ball) s = γ.sill n = γ.nugget diff --git a/test/nesting.jl b/test/nesting.jl index 8a59b4c..e2d0922 100644 --- a/test/nesting.jl +++ b/test/nesting.jl @@ -53,7 +53,7 @@ γ = C₁ * GaussianVariogram(range=1.0) + C₂ * SphericalVariogram(range=2.0) @test range(γ) ≈ 2.0u"m" @test sill(γ) ≈ C₁ .+ C₂ - @test γ(10.0u"m") ≈ sill(γ) + @test γ(10.0) ≈ sill(γ) @test γ(Point(10.0, 0.0), Point(0.0, 0.0)) ≈ sill(γ) @test isstationary(γ) @@ -62,7 +62,7 @@ γ = C * GaussianVariogram() + C * ExponentialVariogram() + C * CubicVariogram() @test range(γ) ≈ 1.0u"m" @test sill(γ) ≈ [3.0 0.0; 0.0 3.0] - @test γ(10.0u"m") ≈ sill(γ) + @test γ(10.0) ≈ sill(γ) @test γ(Point(10.0, 0.0), Point(0.0, 0.0)) ≈ sill(γ) @test isstationary(γ) diff --git a/test/variogram.jl b/test/variogram.jl index a907aa9..ab05cc7 100644 --- a/test/variogram.jl +++ b/test/variogram.jl @@ -1,6 +1,6 @@ @testset "Variogram" begin rng = StableRNG(123) - h = range(0, stop=10, length=50) * u"m" + h = range(0, stop=10, length=50) x, y = rand(rng, Point{3}), rand(rng, Point{3}) # stationary variogram models @@ -61,12 +61,12 @@ # some variograms are non-decreasing for γ in (γnd ∪ [sum(γnd)]) - @test all(γ.(h) .≤ γ.(h .+ 1u"m")) + @test all(γ.(h) .≤ γ.(h .+ 1)) end # variograms are valid at the origin for γ in (γs ∪ γn ∪ γnd) - @test !isnan(γ(0.0u"m")) && !isinf(γ(0.0u"m")) + @test !isnan(γ(0.0)) && !isinf(γ(0.0)) end # practical ranges @@ -97,8 +97,8 @@ # nugget regularization γ = GaussianVariogram(range=20.0, nugget=0.1) C = sill(γ) .- GeoStatsFunctions.pairwise(γ, pset) - @test γ(0u"m") == 0 - @test γ(1e-6u"m") > 0 + @test γ(0) == 0 + @test γ(1e-6) > 0 @test cond(C) < 100.0 # sill and nugget in single precision @@ -125,7 +125,7 @@ # unitful non-stationary types γn = [PowerVariogram(scaling=1.0u"K^2")] for γ in γs - @test unit(γ(1.0u"m")) == u"K^2" + @test unit(γ(1.0)) == u"K^2" end 𝒟 = PointSet(Matrix(1.0I, 3, 3)) From a65885505fbf1a1483893c99b32664b64aab93b7 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Fri, 24 May 2024 08:59:21 -0300 Subject: [PATCH 6/8] Apply suggestions --- src/fitting.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/fitting.jl b/src/fitting.jl index 40f7039..1a3270a 100644 --- a/src/fitting.jl +++ b/src/fitting.jl @@ -132,6 +132,7 @@ function fit_impl( # strip units if necessary ux = unit(eltype(x)) uy = unit(eltype(y)) + x′ = ustrip.(x) y′ = ustrip.(y) # strip units of kwargs @@ -149,7 +150,7 @@ function fit_impl( # objective function function J(θ) γ = V(ball(θ[1]), sill=θ[2], nugget=θ[3]) - sum(i -> w[i] * (γ(x[i]) - y′[i])^2, eachindex(w, x, y′)) + sum(i -> w[i] * (γ(x′[i]) - y′[i])^2, eachindex(w, x′, y′)) end # linear constraint (sill ≥ nugget) @@ -159,7 +160,7 @@ function fit_impl( λ = sum(yᵢ -> yᵢ^2, y′) # maximum range, sill and nugget - xmax = ustrip(maximum(x)) + xmax = maximum(x′) ymax = maximum(y′) rmax = isnothing(maxrange′) ? xmax : maxrange′ smax = isnothing(maxsill′) ? ymax : maxsill′ From 688966bab6cdf4462cd112e622f561844f656182 Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Fri, 24 May 2024 09:03:32 -0300 Subject: [PATCH 7/8] Apply suggestions --- src/fitting.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fitting.jl b/src/fitting.jl index 1a3270a..3af3708 100644 --- a/src/fitting.jl +++ b/src/fitting.jl @@ -129,7 +129,7 @@ function fit_impl( y = y[n .> 0] n = n[n .> 0] - # strip units if necessary + # strip units of coordinates ux = unit(eltype(x)) uy = unit(eltype(y)) x′ = ustrip.(x) From 09e21356c9f469a102d50a038c055832286b4eaf Mon Sep 17 00:00:00 2001 From: Elias Carvalho Date: Fri, 24 May 2024 09:20:55 -0300 Subject: [PATCH 8/8] Add tests --- src/fitting.jl | 2 +- test/fitting.jl | 12 ++++++++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/src/fitting.jl b/src/fitting.jl index 3af3708..4cef8e4 100644 --- a/src/fitting.jl +++ b/src/fitting.jl @@ -186,7 +186,7 @@ function fit_impl( θ = Optim.minimizer(sol) # optimal variogram (with units) - γ = V(ball(θ[1]), sill=θ[2] * uy, nugget=θ[3] * uy) + γ = V(ball(θ[1] * ux), sill=θ[2] * uy, nugget=θ[3] * uy) γ, ϵ end diff --git a/test/fitting.jl b/test/fitting.jl index 878b9ed..d978935 100644 --- a/test/fitting.jl +++ b/test/fitting.jl @@ -63,4 +63,16 @@ γ = GeoStatsFunctions.fit(Variogram, g) @test unit(sill(γ)) == u"K^2" @test unit(nugget(γ)) == u"K^2" + γ = GeoStatsFunctions.fit(Variogram, g, range=12.0u"m") + @test isapprox(range(γ), 12.0u"m", atol=1e-3u"m") + γ = GeoStatsFunctions.fit(Variogram, g, sill=0.06u"K^2") + @test isapprox(sill(γ), 0.06u"K^2", atol=1e-3u"K^2") + γ = GeoStatsFunctions.fit(Variogram, g, nugget=0.02u"K^2") + @test isapprox(nugget(γ), 0.02u"K^2", atol=1e-3u"K^2") + γ = GeoStatsFunctions.fit(Variogram, g, maxrange=6.0u"m") + @test isapprox(range(γ), 6.0u"m", atol=1e-3u"m") + γ = GeoStatsFunctions.fit(Variogram, g, sill=0.03u"K^2") + @test isapprox(sill(γ), 0.03u"K^2", atol=1e-3u"K^2") + γ = GeoStatsFunctions.fit(Variogram, g, nugget=0.01u"K^2") + @test isapprox(nugget(γ), 0.01u"K^2", atol=1e-3u"K^2") end