Program Listing for File Interpolation.h

Program Listing for File Interpolation.h#

Return to documentation for file (include/Karana/Math/Interpolation.h)

/*
 * Copyright (c) 2024-2026 Karana Dynamics Pty Ltd. All rights reserved.
 *
 * NOTICE TO USER:
 *
 * This source code and/or documentation (the "Licensed Materials") is
 * the confidential and proprietary information of Karana Dynamics Inc.
 * Use of these Licensed Materials is governed by the terms and conditions
 * of a separate software license agreement between Karana Dynamics and the
 * Licensee ("License Agreement"). Unless expressly permitted under that
 * agreement, any reproduction, modification, distribution, or disclosure
 * of the Licensed Materials, in whole or in part, to any third party
 * without the prior written consent of Karana Dynamics is strictly prohibited.
 *
 * THE LICENSED MATERIALS ARE PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND.
 * KARANA DYNAMICS DISCLAIMS ALL WARRANTIES, EXPRESS OR IMPLIED, INCLUDING
 * BUT NOT LIMITED TO WARRANTIES OF MERCHANTABILITY, NON-INFRINGEMENT, AND
 * FITNESS FOR A PARTICULAR PURPOSE.
 *
 * IN NO EVENT SHALL KARANA DYNAMICS BE LIABLE FOR ANY DAMAGES WHATSOEVER,
 * INCLUDING BUT NOT LIMITED TO LOSS OF PROFITS, DATA, OR USE, EVEN IF
 * ADVISED OF THE POSSIBILITY OF SUCH DAMAGES, WHETHER IN CONTRACT, TORT,
 * OR OTHERWISE ARISING OUT OF OR IN CONNECTION WITH THE LICENSED MATERIALS.
 *
 * U.S. Government End Users: The Licensed Materials are a "commercial item"
 * as defined at 48 C.F.R. 2.101, and are provided to the U.S. Government
 * only as a commercial end item under the terms of this license.
 *
 * Any use of the Licensed Materials in individual or commercial software must
 * include, in the user documentation and internal source code comments,
 * this Notice, Disclaimer, and U.S. Government Use Provision.
 */

/** @file
 * @brief Various interpolation methods
 */

#pragma once

#include "Karana/Math/Defs.h"
#include <optional>

namespace Karana::Math {
    /** Base class for interpolation methods */
    class BaseInterpolator {
      public:
        /// @brief Constructor.
        BaseInterpolator() = default;

        /// @brief Copy constructor.
        BaseInterpolator(const BaseInterpolator &) = default;

        /// @brief BaseInterpolator virtual destructor
        virtual ~BaseInterpolator() = default;

        /**
         * @brief Sample the dependent variable
         * @param x A value for the independent variable
         * @return The corresponding value for the dependent variable
         */
        virtual double operator()(double x) const = 0;

        /**
         * @brief Sample the dependent variable at multiple values
         * @param x A vector containing values for the independent variable
         * @param out A matching sized vector to contain the result
         */
        virtual void operator()(ConstVecSlice x, VecSlice out) const;

        /**
         * @brief Sample the dependent variable at multiple values
         * @param x A vector containing values for the independent variable
         * @return A corresponding vector of dependent variable values
         */
        virtual Vec operator()(ConstVecSlice x) const;
    };

    /** Constant value implementation of BaseInterpolator
     *
     * This just interpolates to the same constant value everywhere,
     * which may be useful when a BaseInterpolator is expected but only
     * a constant value is needed.
     */
    class ConstantInterpolator : public BaseInterpolator {
      public:
        /**
         * @brief ConstantInterpolator constructor
         * @param constant The constant value to use
         */
        ConstantInterpolator(double constant);

        /**
         * @brief Copy constructor.
         */
        ConstantInterpolator(const ConstantInterpolator &) = default;

        /**
         * @brief Copy assignment operator.
         *
         * @return A reference the recently assigned to ConstantInterpolator.
         */
        ConstantInterpolator &operator=(const ConstantInterpolator &) = default;

        double operator()(double x) const override;
        Vec operator()(ConstVecSlice x) const override;
        void operator()(ConstVecSlice x, VecSlice out) const override;

      private:
        double _constant;
    };

    /** Nearest neighbor interpolation method */
    class NearestNeighborInterpolator : public BaseInterpolator {
      public:
        /**
         * @brief NearestNeighborInterpolator constructor
         * @param indep Monotonically increasing independent variable values
         * @param dep Corresponding dependent variable values
         */
        NearestNeighborInterpolator(ConstVecSlice indep, ConstVecSlice dep);

        /**
         * @brief Copy constructor.
         */
        NearestNeighborInterpolator(const NearestNeighborInterpolator &) = default;

        /**
         * @brief Copy assignment operator.
         *
         * @return A reference the recently assigned to ConstantInterpolator.
         */
        NearestNeighborInterpolator &operator=(const NearestNeighborInterpolator &) = default;

        using BaseInterpolator::operator(); // Bring in overloads we aren't overriding

        /**
         * @brief Call operator to do interpolation.
         *
         * @param x The value to interpolate at.
         * @return The interpolated value.
         */
        double operator()(double x) const override;

      private:
        /**
         * @brief Get the first position the element could be inserted
         *        without changing the ordering. The element currently at that
         *        position is >= to x.
         *
         * @param x The element to find the nearest index to.
         * @return The first index the element could be inserted
         *         into without changing the ordering.
         */
        int _findNearestIndex(double x) const;

        Vec _indep;
        Vec _dep;
    };

    /** Piecewise linear interpolation method */
    class LinearInterpolator : public BaseInterpolator {
      public:
        /** How to extrapolate when sampling past the given data points */
        enum class Extrapolation {
            CONSTANT, ///< Use the value of the boundary data point
            LINEAR,   ///< Linearly extend from the boundary point
        };

        /**
         * @brief LinearInterpolator constructor
         * @param indep Monotonically increasing independent variable values
         * @param dep Corresponding dependent variable values
         * @param extrapolation How to extrapolate outside the data points.
         *  If omitted, extrapolated values are NaN.
         */
        LinearInterpolator(ConstVecSlice indep,
                           ConstVecSlice dep,
                           std::optional<Extrapolation> extrapolation = std::nullopt);

        /**
         * @brief Copy constructor.
         */
        LinearInterpolator(const LinearInterpolator &) = default;

        /**
         * @brief Copy assignment operator.
         *
         * @return A reference the recently assigned to ConstantInterpolator.
         */

        LinearInterpolator &operator=(const LinearInterpolator &) = default;

        using BaseInterpolator::operator(); // Bring in overloads we aren't overriding

        /**
         * @brief Call operator to do interpolation.
         *
         * @param x The value to interpolate at.
         * @return The interpolated value.
         */
        double operator()(double x) const override;

      private:
        double _extendLowerBound(double x) const;
        double _extendUpperBound(double x) const;

        Vec _indep;
        Vec _dep;
        Vec _slopes;

        // How to extrapolate outside the given data points
        std::optional<Extrapolation> _extrapolation;
    };

    /** Akima spline interpolation method
     *
     * Uses a piecewise-cubic spline well-suited to cases where the
     * underlying data has a rapidly varying second derivative.
     * */
    class AkimaSplineInterpolator : public BaseInterpolator {
      public:
        /** How to extrapolate when sampling past the given data points */
        enum class Extrapolation {
            CONSTANT, ///< Use the value of the boundary data point
            LINEAR,   ///< Linearly extend from the boundary point
            CUBIC,    ///< Extend the boundary cubic spline segment
        };

        /**
         * @brief AkimaSplineInterpolator constructor
         * @param indep Strictly increasing independent variable values
         * @param dep Corresponding dependent variable values
         * @param extrapolation How to extrapolate outside the data points.
         *  If omitted, extrapolated values are NaN.
         */
        AkimaSplineInterpolator(ConstVecSlice indep,
                                ConstVecSlice dep,
                                std::optional<Extrapolation> extrapolation = std::nullopt);

        /**
         * @brief Copy constructor.
         */
        AkimaSplineInterpolator(const AkimaSplineInterpolator &);

        /**
         * @brief Copy assignment operator.
         *
         * @return A reference the recently assigned to ConstantInterpolator.
         */
        AkimaSplineInterpolator &operator=(const AkimaSplineInterpolator &);

        using BaseInterpolator::operator(); // Bring in overloads we aren't overriding

        /**
         * @brief Call operator to do interpolation.
         *
         * @param x The value to interpolate at.
         * @return The interpolated value.
         */
        double operator()(double x) const override;

      private:
        double _extendLowerBound(double x) const;
        double _extendUpperBound(double x) const;

        Vec _indep;
        Vec _dep;

        // Slope of the spline at each point.
        // Called `s` in some formulations.
        Vec _spline_slopes;

        // Each spline segment's constant terms
        VecSlice _a;

        // Coefficient for each spline segment's linear term
        VecSlice _b;

        // Coefficient for each spline segment's quadratic term
        Vec _c;

        // Coefficient for each spline segment's cubic term
        Vec _d;

        // How to extrapolate outside the given data points
        std::optional<Extrapolation> _extrapolation;
    };

    /**
     * @class CubicHermiteInterpolator
     * @brief Interpolate using a cubic Hermite spline.
     */
    class CubicHermiteInterpolator : public BaseInterpolator {
      public:
        /** How to extrapolate when sampling past the given data points */
        enum class Extrapolation {
            CONSTANT, ///< Use the value of the boundary data point
            LINEAR,   ///< Linearly extend from the boundary point
            CUBIC,    ///< Extend the boundary cubic spline segment
        };

        /**
         * @brief CubicHermiteInterpolator constructor
         * @param indep Strictly increasing independent variable values
         * @param dep Corresponding dependent variable values
         * @param extrapolation How to extrapolate outside the data points.
         *  If omitted, extrapolated values are NaN.
         * @param deriv Derivative at each indep value; computed using
         *  finite difference if omitted.
         */
        CubicHermiteInterpolator(ConstVecSlice indep,
                                 ConstVecSlice dep,
                                 std::optional<Extrapolation> extrapolation = std::nullopt,
                                 std::optional<ConstVecSlice> deriv = std::nullopt);

        /**
         * @brief Copy constructor.
         */
        CubicHermiteInterpolator(const CubicHermiteInterpolator &) = default;

        /**
         * @brief Copy assignment operator.
         *
         * @return A reference the recently assigned to ConstantInterpolator.
         */
        CubicHermiteInterpolator &operator=(const CubicHermiteInterpolator &) = default;

        using BaseInterpolator::operator(); // Bring in overloads we aren't overriding

        /**
         * @brief Call operator to do interpolation.
         *
         * @param x The value to interpolate at.
         * @return The interpolated value.
         */
        double operator()(double x) const override;

        /**
         * @brief Derivatives at each indep value
         * @return A vector containing the derivatives
         */
        ConstVecSlice derivatives() const;

      private:
        void _calculateDerivatives();
        double _interpolate(size_t seg, double x) const;
        double _extendLowerBound(double x) const;
        double _extendUpperBound(double x) const;

        Vec _indep;
        Vec _dep;
        std::optional<Extrapolation> _extrapolation;
        Vec _deriv;
    };
} // namespace Karana::Math