kriging_rs/variogram/
models.rs

1use crate::Real;
2use crate::error::KrigingError;
3
4/// Matérn semivariance: γ(h) = nugget + partial_sill * (1 - (2^(1-ν)/Γ(ν)) * x^ν * K_ν(x)) with x = h√(2ν)/range.
5fn matern_semivariance(d: Real, nugget: Real, partial_sill: Real, range: Real, nu: Real) -> Real {
6    if d <= 0.0 {
7        return nugget;
8    }
9    let nu_f64 = nu as f64;
10    let x_f64 = (d as f64) * (2.0 * nu_f64).sqrt() / (range as f64);
11    if x_f64 <= 0.0 {
12        return nugget;
13    }
14    let (_i_nu, k_nu) = puruspe::bessel::Inu_Knu(nu_f64, x_f64);
15    let gamma_nu = puruspe::gamma::gamma(nu_f64);
16    let factor = (2.0_f64).powf(1.0 - nu_f64) / gamma_nu * x_f64.powf(nu_f64) * k_nu;
17    let correlation = factor.clamp(0.0, 1.0);
18    nugget + partial_sill * (1.0 - (correlation as Real))
19}
20
21/// Parametric variogram family used to construct a [`VariogramModel`].
22#[derive(Debug, Clone, Copy, PartialEq, Eq)]
23pub enum VariogramType {
24    Spherical,
25    Exponential,
26    Gaussian,
27    Cubic,
28    /// Power-law family; requires shape parameter `alpha` in (0, 2].
29    Stable,
30    /// Matérn family; requires smoothness parameter `nu` > 0.
31    Matern,
32}
33
34/// Parametric variogram (nugget, sill, range, and optional shape).
35///
36/// Build with [`VariogramModel::new`] from a [`VariogramType`]; use with
37/// [`OrdinaryKrigingModel::new`](crate::OrdinaryKrigingModel::new) or
38/// [`BinomialKrigingModel::new`](crate::BinomialKrigingModel::new).
39#[derive(Debug, Clone, Copy, PartialEq)]
40pub enum VariogramModel {
41    Spherical {
42        nugget: Real,
43        sill: Real,
44        range: Real,
45    },
46    Exponential {
47        nugget: Real,
48        sill: Real,
49        range: Real,
50    },
51    Gaussian {
52        nugget: Real,
53        sill: Real,
54        range: Real,
55    },
56    Cubic {
57        nugget: Real,
58        sill: Real,
59        range: Real,
60    },
61    Stable {
62        nugget: Real,
63        sill: Real,
64        range: Real,
65        alpha: Real,
66    },
67    Matern {
68        nugget: Real,
69        sill: Real,
70        range: Real,
71        nu: Real,
72    },
73}
74
75impl VariogramModel {
76    /// Constructs a variogram model with validated parameters: `nugget >= 0`, `sill > nugget`, `range > 0`, all finite.
77    pub fn new(
78        nugget: Real,
79        sill: Real,
80        range: Real,
81        model_type: VariogramType,
82    ) -> Result<Self, KrigingError> {
83        if !nugget.is_finite() || nugget < 0.0 {
84            return Err(KrigingError::FittingError(
85                "nugget must be finite and non-negative".to_string(),
86            ));
87        }
88        if !sill.is_finite() || sill <= nugget {
89            return Err(KrigingError::FittingError(
90                "sill must be finite and greater than nugget".to_string(),
91            ));
92        }
93        if !range.is_finite() || range <= 0.0 {
94            return Err(KrigingError::FittingError(
95                "range must be finite and positive".to_string(),
96            ));
97        }
98        Ok(match model_type {
99            VariogramType::Spherical => VariogramModel::Spherical {
100                nugget,
101                sill,
102                range,
103            },
104            VariogramType::Exponential => VariogramModel::Exponential {
105                nugget,
106                sill,
107                range,
108            },
109            VariogramType::Gaussian => VariogramModel::Gaussian {
110                nugget,
111                sill,
112                range,
113            },
114            VariogramType::Cubic => VariogramModel::Cubic {
115                nugget,
116                sill,
117                range,
118            },
119            VariogramType::Stable => VariogramModel::Stable {
120                nugget,
121                sill,
122                range,
123                alpha: 1.0,
124            },
125            VariogramType::Matern => VariogramModel::Matern {
126                nugget,
127                sill,
128                range,
129                nu: 0.5,
130            },
131        })
132    }
133
134    /// Constructs a variogram model with an explicit shape parameter for Stable (alpha) or Matérn (nu).
135    /// For other model types, `shape` is ignored. Stable: alpha in (0, 2]. Matérn: nu > 0.
136    pub fn new_with_shape(
137        nugget: Real,
138        sill: Real,
139        range: Real,
140        model_type: VariogramType,
141        shape: Real,
142    ) -> Result<Self, KrigingError> {
143        if !nugget.is_finite() || nugget < 0.0 {
144            return Err(KrigingError::FittingError(
145                "nugget must be finite and non-negative".to_string(),
146            ));
147        }
148        if !sill.is_finite() || sill <= nugget {
149            return Err(KrigingError::FittingError(
150                "sill must be finite and greater than nugget".to_string(),
151            ));
152        }
153        if !range.is_finite() || range <= 0.0 {
154            return Err(KrigingError::FittingError(
155                "range must be finite and positive".to_string(),
156            ));
157        }
158        Ok(match model_type {
159            VariogramType::Spherical => VariogramModel::Spherical {
160                nugget,
161                sill,
162                range,
163            },
164            VariogramType::Exponential => VariogramModel::Exponential {
165                nugget,
166                sill,
167                range,
168            },
169            VariogramType::Gaussian => VariogramModel::Gaussian {
170                nugget,
171                sill,
172                range,
173            },
174            VariogramType::Cubic => VariogramModel::Cubic {
175                nugget,
176                sill,
177                range,
178            },
179            VariogramType::Stable => {
180                if !shape.is_finite() || shape <= 0.0 || shape > 2.0 {
181                    return Err(KrigingError::FittingError(
182                        "Stable shape (alpha) must be in (0, 2]".to_string(),
183                    ));
184                }
185                VariogramModel::Stable {
186                    nugget,
187                    sill,
188                    range,
189                    alpha: shape,
190                }
191            }
192            VariogramType::Matern => {
193                if !shape.is_finite() || shape <= 0.0 {
194                    return Err(KrigingError::FittingError(
195                        "Matérn shape (nu) must be positive".to_string(),
196                    ));
197                }
198                VariogramModel::Matern {
199                    nugget,
200                    sill,
201                    range,
202                    nu: shape,
203                }
204            }
205        })
206    }
207
208    pub fn variogram_type(&self) -> VariogramType {
209        match self {
210            Self::Spherical { .. } => VariogramType::Spherical,
211            Self::Exponential { .. } => VariogramType::Exponential,
212            Self::Gaussian { .. } => VariogramType::Gaussian,
213            Self::Cubic { .. } => VariogramType::Cubic,
214            Self::Stable { .. } => VariogramType::Stable,
215            Self::Matern { .. } => VariogramType::Matern,
216        }
217    }
218
219    pub fn params(&self) -> (Real, Real, Real) {
220        match self {
221            Self::Spherical {
222                nugget,
223                sill,
224                range,
225            }
226            | Self::Exponential {
227                nugget,
228                sill,
229                range,
230            }
231            | Self::Gaussian {
232                nugget,
233                sill,
234                range,
235            }
236            | Self::Cubic {
237                nugget,
238                sill,
239                range,
240            }
241            | Self::Stable {
242                nugget,
243                sill,
244                range,
245                ..
246            }
247            | Self::Matern {
248                nugget,
249                sill,
250                range,
251                ..
252            } => (*nugget, *sill, *range),
253        }
254    }
255
256    /// Shape parameter for Stable (alpha) or Matérn (nu). Returns `None` for 3-parameter models.
257    pub fn shape(&self) -> Option<Real> {
258        match self {
259            Self::Stable { alpha, .. } => Some(*alpha),
260            Self::Matern { nu, .. } => Some(*nu),
261            _ => None,
262        }
263    }
264
265    /// Semivariance at `distance`. Assumes `distance >= 0` (e.g. from haversine); clamps negative input to 0.
266    pub fn semivariance(&self, distance: Real) -> Real {
267        let d = distance.max(0.0);
268        let (nugget, sill, range) = self.params();
269        let partial_sill = sill - nugget;
270        let r = range.max(Real::EPSILON);
271
272        match self {
273            Self::Spherical { .. } => {
274                if d >= range {
275                    sill
276                } else {
277                    let x = d / r;
278                    nugget + partial_sill * (1.5 * x - 0.5 * x.powi(3))
279                }
280            }
281            Self::Exponential { .. } => nugget + partial_sill * (1.0 - (-3.0 * d / r).exp()),
282            Self::Gaussian { .. } => {
283                nugget + partial_sill * (1.0 - (-3.0 * (d * d) / (r * r)).exp())
284            }
285            Self::Cubic { .. } => {
286                if d >= range {
287                    sill
288                } else {
289                    let x = d / r;
290                    let poly = 7.0 * x * x - 8.5 * x.powi(3) + 3.5 * x.powi(5) - 0.5 * x.powi(7);
291                    nugget + partial_sill * poly
292                }
293            }
294            Self::Stable { alpha, .. } => {
295                let x = (d / r).powf(*alpha);
296                nugget + partial_sill * (1.0 - (-x).exp())
297            }
298            Self::Matern { nu, .. } => matern_semivariance(d, nugget, partial_sill, r, *nu),
299        }
300    }
301
302    pub fn covariance(&self, distance: Real) -> Real {
303        let (_, sill, _) = self.params();
304        sill - self.semivariance(distance)
305    }
306}
307
308#[cfg(test)]
309mod tests {
310    use super::*;
311    use approx::assert_relative_eq;
312
313    #[test]
314    fn spherical_hits_sill_at_range() {
315        let model = VariogramModel::new(0.1, 1.0, 10.0, VariogramType::Spherical).unwrap();
316        assert_relative_eq!(model.semivariance(10.0), 1.0, epsilon = 1e-6);
317        assert_relative_eq!(model.semivariance(20.0), 1.0, epsilon = 1e-6);
318    }
319
320    #[test]
321    fn exponential_and_gaussian_start_at_nugget() {
322        let exp = VariogramModel::new(0.2, 1.2, 15.0, VariogramType::Exponential).unwrap();
323        let gauss = VariogramModel::new(0.2, 1.2, 15.0, VariogramType::Gaussian).unwrap();
324        assert_relative_eq!(exp.semivariance(0.0), 0.2, epsilon = 1e-6);
325        assert_relative_eq!(gauss.semivariance(0.0), 0.2, epsilon = 1e-6);
326    }
327
328    #[test]
329    fn covariance_complements_semivariance() {
330        let model = VariogramModel::new(0.1, 1.0, 5.0, VariogramType::Exponential).unwrap();
331        let d = 2.2;
332        assert_relative_eq!(
333            model.covariance(d) + model.semivariance(d),
334            1.0,
335            epsilon = 1e-5
336        );
337    }
338
339    #[test]
340    fn cubic_hits_sill_at_range() {
341        let model = VariogramModel::new(0.1, 1.0, 10.0, VariogramType::Cubic).unwrap();
342        assert_relative_eq!(model.semivariance(10.0), 1.0, epsilon = 1e-5);
343        assert_relative_eq!(model.semivariance(20.0), 1.0, epsilon = 1e-5);
344    }
345
346    #[test]
347    fn cubic_stable_matern_start_at_nugget() {
348        let cubic = VariogramModel::new(0.2, 1.2, 15.0, VariogramType::Cubic).unwrap();
349        let stable = VariogramModel::new(0.2, 1.2, 15.0, VariogramType::Stable).unwrap();
350        let matern = VariogramModel::new(0.2, 1.2, 15.0, VariogramType::Matern).unwrap();
351        assert_relative_eq!(cubic.semivariance(0.0), 0.2, epsilon = 1e-6);
352        assert_relative_eq!(stable.semivariance(0.0), 0.2, epsilon = 1e-6);
353        assert_relative_eq!(matern.semivariance(0.0), 0.2, epsilon = 1e-6);
354    }
355
356    #[test]
357    fn stable_with_alpha_one_increases_to_sill() {
358        let stable =
359            VariogramModel::new_with_shape(0.1, 2.0, 10.0, VariogramType::Stable, 1.0).unwrap();
360        let (nugget, sill, _) = stable.params();
361        assert_relative_eq!(stable.semivariance(0.0), nugget, epsilon = 1e-6);
362        let mut prev = nugget;
363        for d in [1.0, 3.0, 5.0, 10.0, 20.0] {
364            let g = stable.semivariance(d);
365            assert!(
366                g >= prev && g <= sill,
367                "stable semivariance should increase toward sill"
368            );
369            prev = g;
370        }
371        assert_relative_eq!(stable.semivariance(100.0), sill, epsilon = 1e-4);
372    }
373
374    #[test]
375    fn shape_returns_none_for_three_param_models() {
376        let m = VariogramModel::new(0.0, 1.0, 1.0, VariogramType::Cubic).unwrap();
377        assert_eq!(m.shape(), None);
378    }
379
380    #[test]
381    fn shape_returns_some_for_stable_and_matern() {
382        let stable =
383            VariogramModel::new_with_shape(0.0, 1.0, 1.0, VariogramType::Stable, 1.5).unwrap();
384        let matern =
385            VariogramModel::new_with_shape(0.0, 1.0, 1.0, VariogramType::Matern, 2.5).unwrap();
386        assert_relative_eq!(stable.shape().unwrap(), 1.5, epsilon = 1e-6);
387        assert_relative_eq!(matern.shape().unwrap(), 2.5, epsilon = 1e-6);
388    }
389}
kriging_rs/variogram/models.rs

kriging_rs/variogram/
models.rs
