ScatteredInterpolator.hpp

Go to the documentation of this file.

#pragma once
 
#include "janus/core/JanusConcepts.hpp"
#include "janus/core/JanusError.hpp"
#include "janus/core/JanusTypes.hpp"
#include "janus/math/Interpolate.hpp"
#include "janus/math/Linalg.hpp"
#include <cmath>
#include <vector>
 
namespace janus {
 
// ============================================================================
// RBF Kernel Types
// ============================================================================

enum class RBFKernel {
    ThinPlateSpline, 
    Multiquadric,    
    Gaussian,        
    Linear,          
    Cubic            
};
 
// ============================================================================
// RBF Kernel Evaluation
// ============================================================================
 
namespace detail {

constexpr double rbf_regularization = 1e-10;
constexpr double rbf_zero_threshold = 1e-15;

inline double rbf_phi(double r, RBFKernel kernel, double epsilon = 1.0) {
    switch (kernel) {
    case RBFKernel::ThinPlateSpline:
        // φ(r) = r² log(r), with φ(0) = 0
        return (r > rbf_zero_threshold) ? r * r * std::log(r) : 0.0;
    case RBFKernel::Multiquadric:
        return std::sqrt(1.0 + (epsilon * r) * (epsilon * r));
    case RBFKernel::Gaussian:
        return std::exp(-(epsilon * r) * (epsilon * r));
    case RBFKernel::Linear:
        return r;
    case RBFKernel::Cubic:
        return r * r * r;
    default:
        return r * r * std::log(r + rbf_zero_threshold);
    }
}

inline double euclidean_distance(const double *p1, const double *p2, int dims) {
    double sum = 0.0;
    for (int d = 0; d < dims; ++d) {
        double diff = p1[d] - p2[d];
        sum += diff * diff;
    }
    return std::sqrt(sum);
}
 
} // namespace detail
 
// ============================================================================
// Scattered Interpolator Class
// ============================================================================

class ScatteredInterpolator {
  private:
    Interpolator m_gridded;
    int m_dims = 0;
    double m_reconstruction_error = 0.0;
    bool m_valid = false;
 
    // Store original data for error computation
    NumericMatrix m_original_points;
    NumericVector m_original_values;
 
  public:
    ScatteredInterpolator() = default;

    ScatteredInterpolator(const NumericMatrix &points, const NumericVector &values,
                          int grid_resolution = 20, RBFKernel kernel = RBFKernel::ThinPlateSpline,
                          double epsilon = 1.0,
                          InterpolationMethod method = InterpolationMethod::Linear) {
        if (points.rows() == 0 || points.cols() == 0) {
            throw InterpolationError("ScatteredInterpolator: points cannot be empty");
        }
        if (points.rows() != values.size()) {
            throw InterpolationError(
                "ScatteredInterpolator: points rows must equal values size. Got " +
                std::to_string(points.rows()) + " points, " + std::to_string(values.size()) +
                " values");
        }
        if (grid_resolution < 2) {
            throw InterpolationError("ScatteredInterpolator: grid_resolution must be >= 2");
        }
 
        m_dims = static_cast<int>(points.cols());
        m_original_points = points;
        m_original_values = values;
 
        // Build uniform grid spanning data range
        std::vector<NumericVector> grid(m_dims);
        for (int d = 0; d < m_dims; ++d) {
            double min_val = points.col(d).minCoeff();
            double max_val = points.col(d).maxCoeff();
            // Add small padding to avoid edge issues
            double padding = 0.01 * (max_val - min_val);
            if (padding < detail::rbf_regularization)
                padding = 0.1; // Handle constant dimension
            min_val -= padding;
            max_val += padding;
 
            grid[d] = NumericVector::LinSpaced(grid_resolution, min_val, max_val);
        }
 
        build_interpolator(points, values, grid, kernel, epsilon, method);
    }

    ScatteredInterpolator(const NumericMatrix &points, const NumericVector &values,
                          const std::vector<NumericVector> &grid_points,
                          RBFKernel kernel = RBFKernel::ThinPlateSpline, double epsilon = 1.0,
                          InterpolationMethod method = InterpolationMethod::Linear) {
        if (points.rows() == 0 || points.cols() == 0) {
            throw InterpolationError("ScatteredInterpolator: points cannot be empty");
        }
        if (points.rows() != values.size()) {
            throw InterpolationError("ScatteredInterpolator: points rows must equal values size");
        }
        if (static_cast<int>(grid_points.size()) != points.cols()) {
            throw InterpolationError(
                "ScatteredInterpolator: grid_points dimensions must match points columns");
        }
 
        m_dims = static_cast<int>(points.cols());
        m_original_points = points;
        m_original_values = values;
 
        build_interpolator(points, values, grid_points, kernel, epsilon, method);
    }

    ScatteredInterpolator(const NumericVector &x, const NumericVector &y, int grid_resolution = 50,
                          RBFKernel kernel = RBFKernel::ThinPlateSpline) {
        if (x.size() != y.size()) {
            throw InterpolationError("ScatteredInterpolator: x and y must have same size");
        }
        if (x.size() < 2) {
            throw InterpolationError("ScatteredInterpolator: need at least 2 points");
        }
 
        m_dims = 1;
        m_original_points.resize(x.size(), 1);
        m_original_points.col(0) = x;
        m_original_values = y;
 
        // Build 1D grid
        double min_val = x.minCoeff();
        double max_val = x.maxCoeff();
        double padding = 0.01 * (max_val - min_val);
        if (padding < detail::rbf_regularization)
            padding = 0.1;
 
        std::vector<NumericVector> grid(1);
        grid[0] = NumericVector::LinSpaced(grid_resolution, min_val - padding, max_val + padding);
 
        build_interpolator(m_original_points, y, grid, kernel, 1.0, InterpolationMethod::Linear);
    }
 
    // ========================================================================
    // Query Methods
    // ========================================================================

    template <JanusScalar Scalar> Scalar operator()(const JanusVector<Scalar> &query) const {
        if (!m_valid) {
            throw InterpolationError("ScatteredInterpolator: not initialized");
        }
        return m_gridded(query);
    }

    template <JanusScalar Scalar> Scalar operator()(const Scalar &query) const {
        if (!m_valid) {
            throw InterpolationError("ScatteredInterpolator: not initialized");
        }
        if (m_dims != 1) {
            throw InterpolationError("ScatteredInterpolator: scalar query only valid for 1D");
        }
        return m_gridded(query);
    }
 
    // ========================================================================
    // Properties
    // ========================================================================

    int dims() const { return m_dims; }

    bool valid() const { return m_valid; }

    const Interpolator &gridded() const { return m_gridded; }

    double reconstruction_error() const { return m_reconstruction_error; }
 
  private:
    void build_interpolator(const NumericMatrix &points, const NumericVector &values,
                            const std::vector<NumericVector> &grid, RBFKernel kernel,
                            double epsilon, InterpolationMethod method) {
        int n_points = static_cast<int>(points.rows());
        int n_dims = static_cast<int>(points.cols());
 
        // ====================================================================
        // Step 1: Build RBF interpolation matrix and solve for weights
        // ====================================================================
        // For thin plate spline in 2D, we need polynomial terms for uniqueness
        // General form: f(x) = Σ w_i φ(||x - x_i||) + p(x)
        // We'll use the simpler pure RBF form for now
 
        NumericMatrix Phi(n_points, n_points);
        for (int i = 0; i < n_points; ++i) {
            for (int j = 0; j < n_points; ++j) {
                double r =
                    detail::euclidean_distance(points.row(i).data(), points.row(j).data(), n_dims);
                Phi(i, j) = detail::rbf_phi(r, kernel, epsilon);
            }
        }
 
        // Add small regularization for numerical stability
        Phi.diagonal().array() += detail::rbf_regularization;
 
        // Solve Phi * weights = values
        NumericVector weights = solve(Phi, values);
 
        // ====================================================================
        // Step 2: Evaluate RBF at grid points
        // ====================================================================
        // Compute total grid size
        int grid_size = 1;
        for (int d = 0; d < n_dims; ++d) {
            grid_size *= static_cast<int>(grid[d].size());
        }
 
        // Flatten grid values (Fortran order for Interpolator)
        NumericVector grid_values(grid_size);
 
        // Iterate over all grid points
        std::vector<int> dims_sizes(n_dims);
        for (int d = 0; d < n_dims; ++d) {
            dims_sizes[d] = static_cast<int>(grid[d].size());
        }
 
        for (int flat_idx = 0; flat_idx < grid_size; ++flat_idx) {
            // Convert flat index to multi-index (Fortran order)
            std::vector<int> multi_idx(n_dims);
            int remaining = flat_idx;
            for (int d = 0; d < n_dims; ++d) {
                multi_idx[d] = remaining % dims_sizes[d];
                remaining /= dims_sizes[d];
            }
 
            // Get grid point coordinates
            std::vector<double> grid_pt(n_dims);
            for (int d = 0; d < n_dims; ++d) {
                grid_pt[d] = grid[d](multi_idx[d]);
            }
 
            // Evaluate RBF sum at this grid point
            double val = 0.0;
            for (int i = 0; i < n_points; ++i) {
                double r = detail::euclidean_distance(grid_pt.data(), points.row(i).data(), n_dims);
                val += weights(i) * detail::rbf_phi(r, kernel, epsilon);
            }
            grid_values(flat_idx) = val;
        }
 
        // ====================================================================
        // Step 3: Create gridded interpolator
        // ====================================================================
        m_gridded = Interpolator(grid, grid_values, method);
        m_valid = true;
 
        // ====================================================================
        // Step 4: Compute reconstruction error
        // ====================================================================
        compute_reconstruction_error();
    }

    void compute_reconstruction_error() {
        if (!m_valid)
            return;
 
        double sum_sq_error = 0.0;
        int n = static_cast<int>(m_original_values.size());
 
        for (int i = 0; i < n; ++i) {
            NumericVector query = m_original_points.row(i).transpose();
            double predicted = m_gridded(query);
            double error = predicted - m_original_values(i);
            sum_sq_error += error * error;
        }
 
        m_reconstruction_error = std::sqrt(sum_sq_error / n);
    }
};
 
} // namespace janus