Signature	Description
enum class VariogramModel : unsigned char { spherical = 1, // Bounded - Most common default. Linear near origin. exponential = 2, // Bounded - Gradual decay. Common in hydrology. gaussian = 3, // Bounded - Parabolic near origin - very smooth fields. matern = 4, // Bounded - Generalizes exponential/Gaussian via smooth. power = 5, // Unbounded - Non-stationary / fractal processes. linear = 6, // Unbounded - Special case of power with exponent = 1. logarithmic = 7, // Unbounded - Intrinsic random functions of order 1. };	The theoretical variogram/covariance model used to characterize spatial (or spatio-temporal) correlation as a function of separation distance h.
template<typename T = double> struct KrigingParams { using value_type = T; // Which theoretical variogram/covariance family to use. // VariogramModel model { VariogramModel::spherical }; // Nugget (c0): discontinuity at h=0 -- measurement error / micro-scale // variability. Must be >= 0. // value_type nugget { 0 }; // Partial sill (c): structural variance above the nugget. Total sill is // (nugget + sill). Ignored by power/linear/logarithmic (unbounded). // Must be >= 0. // value_type sill { 1 }; // Range (a): practical correlation distance for bounded models. Must be // > 0. For Matérn this acts as the length-scale ρ. // value_type range { 1 }; // Matérn smoothness parameter ν (nu). Common values: 0.5 (== exponential) // 1.5, 2.5, and the limit ν -> infinity (== Gaussian). Only used when // model == matern. Must be > 0. // value_type matern_smoothness { 1.5 }; // Exponent α for the power model, 0 < α < 2. α == 1 reduces to the // linear model (also selectable directly via model == linear). // Only used when model == power. // value_type power_exponent { 1.0 }; // Small positive jitter added to the diagonal of the kriging system // before factorization, purely for numerical stability (handles // near-duplicate coordinates / ill-conditioning). This is independent of // the statistical nugget effect above. // value_type ridge { 1e-10 }; // If true, perform Ordinary Kriging (unknown constant mean, estimated // via the Lagrange-multiplier-augmented system). If false, perform // Simple Kriging assuming a known mean (simple_mean below) with zero-mean // residuals -- this drops the Lagrange row/column and is cheaper. // bool ordinary { true }; // Used only when ordinary == false. // value_type simple_mean { 0 }; // Geometric anisotropy: stretch distances along a preferred direction // before evaluating the variogram. anisotropy_ratio is range_minor / // range_major (1 == isotropic). anisotropy_angle is the direction of the // major axis in radians, measured counter-clockwise from the x-axis. // Only meaningful for 2D coordinates. // value_type anisotropy_ratio { 1.0 }; value_type anisotropy_angle { 0.0 }; // If > 0, switches to local (moving-window) Kriging: only the // max_neighbors closest known points are used to predict each unknown // location, turning an O(N^3) global solve into many small solves -- // essential for large N. 0 means use all points (global Kriging). // std::size_t max_neighbors { 0 }; inline void sanity_check() const { #ifdef HMDF_SANITY_EXCEPTIONS if (nugget < 0) throw DataFrameError("KrigingParams: nugget must be >= 0"); if (sill < 0) throw DataFrameError("KrigingParams: sill must be >= 0"); if (range <= 0) throw DataFrameError("KrigingParams: range must be > 0"); if (model == VariogramModel::matern && matern_smoothness <= 0) throw DataFrameError("KrigingParams: matern_smoothness must be > 0"); if (model == VariogramModel::power && (power_exponent <= 0 \|\| power_exponent >= 2)) throw DataFrameError("KrigingParams: power_exponent must be in (0, 2)"); if (anisotropy_ratio <= 0 \|\| anisotropy_ratio > 1.0) throw DataFrameError("KrigingParams: anisotropy_ratio must be in (0, 1]"); #endif // HMDF_SANITY_EXCEPTIONS } };	All tunable hyper-parameters of a Kriging model are grouped into this one structure, decoupling "what the model is" from "how it is computed". The same KrigingVisitor implementation dispatches on params_.model at runtime, so switching models is a one-field change with no recompilation of logic.

Signature

Description

enum class  VariogramModel : unsigned char  {

    spherical = 1,    // Bounded - Most common default. Linear near origin.
    exponential = 2,  // Bounded - Gradual decay. Common in hydrology.
    gaussian = 3,     // Bounded - Parabolic near origin - very smooth fields.
    matern = 4,       // Bounded - Generalizes exponential/Gaussian via smooth.
    power = 5,        // Unbounded - Non-stationary / fractal processes.
    linear = 6,       // Unbounded - Special case of power with exponent = 1.
    logarithmic = 7,  // Unbounded - Intrinsic random functions of order 1.
};

The theoretical variogram/covariance model used to characterize spatial (or spatio-temporal) correlation as a function of separation distance h.

template<typename T = double>
struct  KrigingParams  {

    using value_type = T;

    // Which theoretical variogram/covariance family to use.
    //
    VariogramModel  model { VariogramModel::spherical };

    // Nugget (c0): discontinuity at h=0 -- measurement error / micro-scale
    // variability. Must be >= 0.
    //
    value_type      nugget { 0 };

    // Partial sill (c): structural variance above the nugget. Total sill is
    // (nugget + sill). Ignored by power/linear/logarithmic (unbounded).
    // Must be >= 0.
    //
    value_type      sill { 1 };

    // Range (a): practical correlation distance for bounded models. Must be
    // > 0. For Matérn this acts as the length-scale ρ.
    //
    value_type      range { 1 };

    // Matérn smoothness parameter ν (nu). Common values: 0.5 (== exponential)
    // 1.5, 2.5, and the limit ν -> infinity (== Gaussian). Only used when
    // model == matern. Must be > 0.
    //
    value_type      matern_smoothness { 1.5 };

    // Exponent α for the power model, 0 < α < 2.  α == 1 reduces to the
    // linear model (also selectable directly via model == linear).
    // Only used when model == power.
    //
    value_type      power_exponent { 1.0 };

    // Small positive jitter added to the diagonal of the kriging system
    // before factorization, purely for numerical stability (handles
    // near-duplicate coordinates / ill-conditioning). This is independent of
    // the statistical nugget effect above.
    //
    value_type      ridge { 1e-10 };

    // If true, perform Ordinary Kriging (unknown constant mean, estimated
    // via the Lagrange-multiplier-augmented system). If false, perform
    // Simple Kriging assuming a known mean (simple_mean below) with zero-mean
    // residuals -- this drops the Lagrange row/column and is cheaper.
    //
    bool            ordinary { true };

    // Used only when ordinary == false.
    //
    value_type      simple_mean { 0 };

    // Geometric anisotropy: stretch distances along a preferred direction
    // before evaluating the variogram. anisotropy_ratio is range_minor /
    // range_major (1 == isotropic). anisotropy_angle is the direction of the
    // major axis in radians, measured counter-clockwise from the x-axis.
    // Only meaningful for 2D coordinates.
    //
    value_type      anisotropy_ratio { 1.0 };
    value_type      anisotropy_angle { 0.0 };

    // If > 0, switches to local (moving-window) Kriging: only the
    // max_neighbors closest known points are used to predict each unknown
    // location, turning an O(N^3) global solve into many small solves --
    // essential for large N. 0 means use all points (global Kriging).
    //
    std::size_t     max_neighbors { 0 };

    inline void sanity_check() const  {

#ifdef HMDF_SANITY_EXCEPTIONS
        if (nugget < 0)
            throw DataFrameError("KrigingParams: nugget must be >= 0");
        if (sill < 0)
            throw DataFrameError("KrigingParams: sill must be >= 0");
        if (range <= 0)
            throw DataFrameError("KrigingParams: range must be > 0");
        if (model == VariogramModel::matern && matern_smoothness <= 0)
            throw DataFrameError("KrigingParams: matern_smoothness must be > 0");
        if (model == VariogramModel::power && (power_exponent <= 0 || power_exponent >= 2))
            throw DataFrameError("KrigingParams: power_exponent must be in (0, 2)");
        if (anisotropy_ratio <= 0 || anisotropy_ratio > 1.0)
            throw DataFrameError("KrigingParams: anisotropy_ratio must be in (0, 1]");
#endif // HMDF_SANITY_EXCEPTIONS
    }
};

All tunable hyper-parameters of a Kriging model are grouped into this one structure, decoupling "what the model is" from "how it is computed". The same KrigingVisitor implementation dispatches on params_.model at runtime, so switching models is a one-field change with no recompilation of logic.

Signature	Description	Parameters
#include <DataFrame/DataFrameMLVisitors.h> template<std::size_t DIM, typename T, typename I = unsigned long, std::size_t A = 0> struct KrigingVisitor; // ------------------------------------- template<std::size_t DIM, typename T, typename I = unsigned long, std::size_t A = 0> using krig_v = KrigingVisitor<DIM, T, I, A>;	This is a "single action visitor", meaning it is passed the whole data vector in one call and you must use the single_act_visit() interface. Kriging algorithm performs spatial interpolation (a.k.a. Gaussian Process regression with a fixed/known covariance kernel) over scattered data. Kriging models are advanced geostatistical interpolation techniques used to predict unknown values at unsampled locations based on a set of known, scattered data points. Unlike simple averaging, Kriging calculates weights using both the distance between points and the overall spatial arrangement of the data, while also estimating the uncertainty of the predictions. Kriging sits at the intersection of interpolation and probabilistic regression — it's better categorized as a spatial interpolation / Gaussian Process regression algorithm. A fitting algorithm finds a compact parametric model (a line, a curve, coefficients) that summarizes the data, accepting residuals. It reduces the data to a small set of parameters. Kriging does the opposite: it retains all training points, and every prediction is a weighted combination of all observed values. There are no "fitted parameters" — nugget, sill, range are supplied by the user, not estimated from the data. The output is not a model summary but an interpolated surface. The visitor takes two parallel columns as inputs: Coordinates — a column of std::array/std::vector values, one per data point. Each element of the array is one spatial dimension, so DIM=1 is a 1D transect, DIM=2 is a 2D scatter (easting/northing, lon/lat, x/y), DIM=3 is volumetric, and so on. The dimensionality is a compile-time template parameter. Observations — a parallel column of scalar T values, one per data point. These are the measured quantities at each coordinate — the "z" values you want to interpolate. So the mental model is: given N scattered points each with a known location and a known value, fit the Kriging system so you can later call predict() at any arbitrary location and get back a value plus its uncertainty. This visitor has the following methods to get results: // Point estimate + Kriging variance (the latter quantifies prediction // uncertainty -- near-zero close to known points, growing with // distance / sparsity). // struct Estimate { value_type value { 0 }; value_type variance { 0 }; }; using result_type = vec_t<Estimate>; // Predict the value (and its Kriging variance) at an arbitrary query // point. Safe to call repeatedly and concurrently (read-only) once the // visitor has been invoked via operator() / fit(). // Estimate predict(const coord_type &query_pt) const; // Convenience batch API -- predicts every point in [begin, end) and // returns results in the same order. Each prediction only re-solves // against the cached factorization (or local neighborhood), so this is // O(m * n) for m query points rather than O(m * n^3). // template<typename QV> result_type predict_many (const QV &begin, const QV &end) const; // Leave-one-out cross validation: refits N times, each time omitting // one observation, and reports the predicted-vs-actual residual. Useful // for comparing variogram models/parameters objectively. This is O(N) // refits of an O(N^3) system, so intended for moderate N or via // max_neighbors-restricted local Kriging. // vec_t<value_type> loo_cross_validate () const; explicit KrigingVisitor(const KrigingParams ¶ms = { }); : params_(params) params: Necessary parameters as explained above.	DIM: Dimensionality of coordinates. T: Column data type. I: Index type. A: Memory alignment boundary for vectors. Default is system default alignment

static void test_KrigingVisitor()  {

    std::cout << "\nTesting KrigingVisitor{ } ..." << std::endl;

    ULDataFrame df;

    try  {
        df.read("FORD.csv", io_format::csv2);
    }
    catch (const DataFrameError &ex)  {
        std::cout << ex.what() << std::endl;
        ::exit(-1);
    }

    // Duplicate points
    //
    std::vector<std::array<double, 2>>  coords { { 0.0, 0.0 }, { 0.0, 0.0 }, { 5.0, 5.0 }, { 1.0, 1.0 } };
    std::vector<double>                 obs { 1.0, 1.2, 4.0, 1.5 };

    df.load_column("COORDS 1", std::move(coords), nan_policy::dont_pad_with_nans);
    df.load_column("OBSERV 1", std::move(obs), nan_policy::dont_pad_with_nans);

    KrigingParams<double>  params;

    params.model = VariogramModel::spherical;
    params.range = 3.0;
    params.ridge = 1e-8;

    KrigingVisitor<2, double>   kv1 { params };

    df.single_act_visit<std::array<double, 2>, double>("COORDS 1", "OBSERV 1", kv1);

    const auto  estimate1 { kv1.predict({ 0.0, 0.0 }) };
    const auto  leave_one_out1 { kv1.loo_cross_validate() };

    assert(std::abs(estimate1.value - 1.1) < 0.001);
    assert(estimate1.variance < 0.00000001);
    assert(leave_one_out1.size() == 4);
    assert(std::abs(leave_one_out1[0] - -0.2) < 0.001);
    assert(std::abs(leave_one_out1[1] - 0.2) < 0.001);
    assert(std::abs(leave_one_out1[3] - -0.549356) < 0.000001);

    // Filtering NaN
    //
    constexpr double   qnan = std::numeric_limits<double>::quiet_NaN();

    std::vector<std::array<double, 2>>  coords2 { { 0.0, 0.0 }, { 1.0, 1.0 }, { qnan, 2.0 }, { 3.0, 3.0 }, { 4.0, 4.0 } };
    std::vector<double>                 obs2 { 1.0, 2.0, 3.0, qnan, 5.0 };

    df.load_column("COORDS 2", std::move(coords2), nan_policy::dont_pad_with_nans);
    df.load_column("OBSERV 2", std::move(obs2), nan_policy::dont_pad_with_nans);

    KrigingParams<double>   params2;

    params2.model = VariogramModel::exponential;
    params2.range = 2.0;

    KrigingVisitor<2, double>   kv2 { params2 };

    df.single_act_visit<std::array<double, 2>, double>("COORDS 2", "OBSERV 2", kv2);

    const auto  estimate2 { kv2.predict({ 0.0, 0.0 }) };
    const auto  leave_one_out2 { kv2.loo_cross_validate() };

    assert(std::abs(estimate2.value - 1.0) < 0.001);
    assert(estimate2.variance < 0.000000001);
    assert(leave_one_out2.size() == 3);
    assert(std::abs(leave_one_out2[0] - -1.7604) < 0.0001);
    assert(std::abs(leave_one_out2[1] - -0.206722) < 0.000001);
    assert(std::abs(leave_one_out2[2] - 3.44006) < 0.00001);

    // 3D
    //
    std::vector<std::array<double, 3>>  coords3;
    std::vector<double>                 obs3;

    for (int i = 0; i < 30; ++i)  {
        const double    x = (i % 3) * 2.0;
        const double    y = ((i / 3) % 3) * 2.0;
        const double    z = (i / 9) * 2.0;

        coords3.push_back({ x, y, z });
        obs3.push_back(x + 2.0 * y - z);
    }
    df.load_column("COORDS 3", std::move(coords3), nan_policy::dont_pad_with_nans);
    df.load_column("OBSERV 3", std::move(obs3), nan_policy::dont_pad_with_nans);

    KrigingParams<double>   params3;

    params3.model = VariogramModel::matern;
    params3.matern_smoothness = 2.5;
    params3.range = 4.0;

    KrigingVisitor<3, double>   kv3 { params3 };

    df.single_act_visit<std::array<double, 3>, double>("COORDS 3", "OBSERV 3", kv3);

    const auto  estimate3 { kv3.predict({ 1.0, 1.0, 1.0 }) };
    const auto  leave_one_out3 { kv3.loo_cross_validate() };

    assert(std::abs(estimate3.value - 1.60219) < 0.00001);
    assert(std::abs(estimate3.variance - 0.018453) < 0.000001);
    assert(leave_one_out3.size() == 30);
    assert(std::abs(leave_one_out3[0] - -0.276349) < 0.000001);
    assert(std::abs(leave_one_out3[10] - 0.030186) < 0.000001);
    assert(std::abs(leave_one_out3[29] - -0.800075) < 0.000001);

    // Generic Bessel-K path
    //
    std::vector<std::array<double, 1>>  coords4 { { 0.0 }, { 1.0 }, { 2.0 }, { 3.0 }, { 4.0 }, { 5.0 } };
    std::vector<double>                 obs4 { 0.0, 1.0, 0.5, 1.5, 1.0, 2.0 };

    df.load_column("COORDS 4", std::move(coords4), nan_policy::dont_pad_with_nans);
    df.load_column("OBSERV 4", std::move(obs4), nan_policy::dont_pad_with_nans);
    for (double nu : { 0.7, 1.0, 1.8, 3.3 })  {
        KrigingParams<double>  params;

        params.model = VariogramModel::matern;
        params.matern_smoothness = nu;
        params.range = 2.0;

        KrigingVisitor<1, double>     kv { params };

        df.single_act_visit<std::array<double, 1>, double>("COORDS 4", "OBSERV 4", kv);

        const auto  estimate = kv.predict({ 2.5 });
        const auto  leave_one_out { kv.loo_cross_validate() };

        assert(leave_one_out.size() == 6);
        assert(std::abs(estimate.value - 1.0) < 0.001);
        if (nu == 0.7)  {
            assert(std::abs(estimate.variance - 0.38837) < 0.00001);
            assert(std::abs(leave_one_out[0] - -1.08785) < 0.00001);
            assert(std::abs(leave_one_out[3] - 0.663442) < 0.000001);
            assert(std::abs(leave_one_out[5] - 1.08785) < 0.00001);
        }
        else if (nu == 1.0)  {
            assert(std::abs(estimate.variance - 0.072554) < 0.000001);
            assert(std::abs(leave_one_out[0] - -1.25466) < 0.00001);
            assert(std::abs(leave_one_out[3] - 0.930086) < 0.000001);
            assert(std::abs(leave_one_out[5] - 1.25466) < 0.00001);
        }
        else if (nu == 1.8)  {
            assert(std::abs(estimate.variance - 0.593665) < 0.000001);
            assert(std::abs(leave_one_out[0] - -1.11348) < 0.00001);
            assert(std::abs(leave_one_out[3] - 0.682445) < 0.000001);
            assert(std::abs(leave_one_out[5] - 1.11348) < 0.00001);
        }
        else if (nu == 3.3)  {
            assert(std::abs(estimate.variance - 0.872824) < 0.000001);
            assert(std::abs(leave_one_out[0] - -1.02868) < 0.00001);
            assert(std::abs(leave_one_out[3] - 0.443655) < 0.000001);
            assert(std::abs(leave_one_out[5] - 1.02868) < 0.00001);
        }
    }
}