Add pentaspherical and sine hole variogram models

docs/src/theoretical_variograms.md

@@ -101,6 +101,16 @@ SphericalVariogram
 CubicVariogram
 ```
 
+## Pentaspherical
+
+```math
+\gamma(h) = (s - n) \left[\left(\frac{15}{8}\left(\frac{h}{r}\right) - \frac{5}{4}\left(\frac{h}{r}\right)^3 + \frac{3}{8}\left(\frac{h}{r}\right)^5\right) \cdot \1_{(0,r)}(h) + \1_{[r,\infty)}(h)\right] + n \cdot \1_{(0,\infty)}(h)
+```
+
+```@docs
+PentasphericalVariogram
+```
+
 ## Power
 
 ```math
@@ -111,6 +121,16 @@ CubicVariogram
 PowerVariogram
 ```
 
+## Sine hole
+
+```math
+\gamma(h) = (s - n) \left[1 - \frac{\sin(\pi h / r)}{\pi h / r}\right] + n \cdot \1_{(0,\infty)}(h)
+```
+
+```@docs
+SineHoleVariogram
+```
+
 ## Composite
 
 ```@docs

src/GeoStats.jl

@@ -53,7 +53,9 @@ export
   MaternVariogram,
   SphericalVariogram,
   CubicVariogram,
+  PentasphericalVariogram,
   PowerVariogram,
+  SineHoleVariogram,
   CompositeVariogram,
 
   # estimators

src/theoretical_variograms.jl

@@ -135,6 +135,29 @@ end
 (γ::CubicVariogram)(x, y) = γ(γ.distance(x, y))
 isstationary(::CubicVariogram) = true
 
+"""
+    PentasphericalVariogram
+
+A pentaspherical variogram with sill `s`, range `r` and nugget `n`.
+Optionally, use a custom distance `d`.
+"""
+@with_kw struct PentasphericalVariogram{T<:Real,D<:AbstractDistance} <: AbstractVariogram
+  sill::T   = 1.
+  range::T  = 1.
+  nugget::T = 0.
+  distance::D = EuclideanDistance()
+end
+(γ::PentasphericalVariogram)(h) = begin
+  s = γ.sill
+  r = γ.range
+  n = γ.nugget
+
+  (h .< r) .* (s - n) .* ((15/8)*(h/r) - (5/4)*(h/r).^3 + (3/8)*(h/r).^5) +
+  (h .≥ r) .* (s - n) + n
+end
+(γ::PentasphericalVariogram)(x, y) = γ(γ.distance(x, y))
+isstationary(::PentasphericalVariogram) = true
+
 """
     PowerVariogram(scaling=s, exponent=a, nugget=n, distance=d)
 
@@ -150,6 +173,22 @@ end
 (γ::PowerVariogram)(h) = γ.scaling*h.^γ.exponent + γ.nugget
 (γ::PowerVariogram)(x, y) = γ(γ.distance(x, y))
 
+"""
+    SineHoleVariogram(sill=s, range=r, nugget=n, distance=d)
+
+A sine hole variogram with sill `s`, range `r` and nugget `n`.
+Optionally, use a custom distance `d`.
+"""
+@with_kw struct SineHoleVariogram{T<:Real,D<:AbstractDistance} <: AbstractVariogram
+  sill::T   = 1.
+  range::T  = 1.
+  nugget::T = 0.
+  distance::D = EuclideanDistance()
+end
+(γ::SineHoleVariogram)(h) = (γ.sill - γ.nugget) * (1 - sin.(π*h/γ.range)./(π*h/γ.range)) + γ.nugget
+(γ::SineHoleVariogram)(x, y) = γ(γ.distance(x, y))
+isstationary(::SineHoleVariogram) = true
+
 """
     CompositeVariogram(γ₁, γ₂, ..., γₙ)
 

test/theoretical_variograms.jl

@@ -5,32 +5,44 @@
   # stationary variogram models
   γs = [GaussianVariogram(), ExponentialVariogram(),
         MaternVariogram(), SphericalVariogram(),
-        SphericalVariogram(range=2.), CubicVariogram()]
+        SphericalVariogram(range=2.), CubicVariogram(),
+        PentasphericalVariogram(), SineHoleVariogram()]
 
   # non-stationary variogram models
-  γn = [PowerVariogram()]
+  γn = [PowerVariogram(), PowerVariogram(exponent=.4)]
+
+  # non-decreasing variogram models
+  γnd = [GaussianVariogram(), ExponentialVariogram(),
+         MaternVariogram(), SphericalVariogram(),
+         SphericalVariogram(range=2.), CubicVariogram(),
+         PentasphericalVariogram(), PowerVariogram()]
 
   # check stationarity
   @test all(isstationary(γ) for γ ∈ γs)
-
-  # check non-stationarity
   @test all(!isstationary(γ) for γ ∈ γn)
 
-  for γ in [γs..., CompositeVariogram(γs...)]
-    # variograms are increasing
-    @test all(γ(h) .≤ γ(h+1))
-
-    # variograms are symmetric
+  # variograms are symmetric under Euclidean distance
+  for γ ∈ (γs ∪ γn ∪ γnd ∪ [CompositeVariogram(γs..., γn..., γnd...)])
     @test γ(x, y) ≈ γ(y, x)
   end
 
+  # some variograms are non-decreasing
+  for γ ∈ (γnd ∪ [CompositeVariogram(γnd...)])
+    @test all(γ(h) .≤ γ(h+1))
+  end
+
   if ismaintainer || istravis
     @testset "Plot recipe" begin
       function plot_variograms(fname)
-        plt = plot()
+        plt1 = plot()
         for γ ∈ γs
-          plot!(plt, γ, maxlag=3.)
+          plot!(plt1, γ, maxlag=3.)
+        end
+        plt2 = plot()
+        for γ ∈ γn
+          plot!(plt2, γ, maxlag=3.)
         end
+        plot(plt1, plt2, size=(1000,400))
         png(fname)
       end
       refimg = joinpath(datadir,"TheoreticalVariograms.png")