alglib_wrapper.cpp

/*
 * alglib_wrapper.cpp - C wrapper for ALGLIB spline functions
 *
 * Provides C-callable interface to ALGLIB's spline1d functions for use in
 * swirl2d_ipopt.c. Supports cubic, Akima, monotone, Hermite, and Catmull-Rom spline types.
 */

#include "alglib-cpp/src/interpolation.h"
#include <cstdlib>
#include <cstring>
#include <cmath>

/* Spline type constants (must match AlglibSplineType in swirl2d_ipopt.c) */
#define ALGLIB_CUBIC      0
#define ALGLIB_AKIMA      1
#define ALGLIB_MONOTONE   2
#define ALGLIB_HERMITE    3
#define ALGLIB_CATMULLROM 4

/* Global boundary condition parameters for cubic spline */
static int g_alglib_bc_type_left = 2;    /* 0=parabolic, 1=1st deriv, 2=2nd deriv */
static int g_alglib_bc_type_right = 2;
static double g_alglib_bc_value_left = 0.0;
static double g_alglib_bc_value_right = 0.0;

/* Wrapper structure to hold ALGLIB spline object */
struct alglib_spline_wrapper {
    alglib::spline1dinterpolant spline;
};

extern "C" {

/*
 * Set boundary condition parameters for cubic spline
 * bc_type: 0=parabolically-terminated, 1=first derivative, 2=second derivative
 */
void alglib_set_bc(int bc_type_left, double bc_value_left,
                   int bc_type_right, double bc_value_right)
{
    g_alglib_bc_type_left = bc_type_left;
    g_alglib_bc_value_left = bc_value_left;
    g_alglib_bc_type_right = bc_type_right;
    g_alglib_bc_value_right = bc_value_right;
}

/*
 * Estimate derivatives for Hermite spline using finite differences
 * Uses 3-point formula at interior, one-sided at endpoints
 */
static void estimate_derivatives(const double *x, const double *y, int n, double *d)
{
    if (n < 2) {
        d[0] = 0.0;
        return;
    }

    /* Left endpoint: forward difference */
    d[0] = (y[1] - y[0]) / (x[1] - x[0]);

    /* Interior points: central difference */
    for (int i = 1; i < n - 1; i++) {
        double h1 = x[i] - x[i-1];
        double h2 = x[i+1] - x[i];
        /* Weighted average for non-uniform spacing */
        d[i] = (h2 * (y[i] - y[i-1]) / h1 + h1 * (y[i+1] - y[i]) / h2) / (h1 + h2);
    }

    /* Right endpoint: backward difference */
    d[n-1] = (y[n-1] - y[n-2]) / (x[n-1] - x[n-2]);
}

/*
 * Build an ALGLIB spline from data points
 *
 * Parameters:
 *   x, y      - input data arrays (must be sorted by x)
 *   n         - number of data points
 *   spline_type - 0=cubic, 1=akima, 2=monotone, 3=hermite, 4=catmullrom
 *   spline_ptr  - output pointer to allocated spline object
 */
void alglib_spline_build(const double *x, const double *y, int n,
                         int spline_type, void **spline_ptr)
{
    /* Allocate wrapper structure */
    alglib_spline_wrapper *wrapper = new alglib_spline_wrapper();

    /* Convert C arrays to ALGLIB arrays */
    alglib::real_1d_array xa, ya;
    xa.setcontent(n, x);
    ya.setcontent(n, y);

    /* Build the appropriate spline type */
    switch (spline_type) {
        case ALGLIB_CUBIC:
            /* Cubic spline with configurable boundary conditions */
            alglib::spline1dbuildcubic(xa, ya, n,
                                       g_alglib_bc_type_left, g_alglib_bc_value_left,
                                       g_alglib_bc_type_right, g_alglib_bc_value_right,
                                       wrapper->spline);
            break;

        case ALGLIB_AKIMA:
            /* Akima spline - local, avoids overshoot */
            alglib::spline1dbuildakima(xa, ya, n, wrapper->spline);
            break;

        case ALGLIB_MONOTONE:
            /* Monotone spline - shape-preserving */
            alglib::spline1dbuildmonotone(xa, ya, n, wrapper->spline);
            break;

        case ALGLIB_HERMITE:
            {
                /* Hermite spline with estimated derivatives */
                double *d = new double[n];
                estimate_derivatives(x, y, n, d);
                alglib::real_1d_array da;
                da.setcontent(n, d);
                alglib::spline1dbuildhermite(xa, ya, da, n, wrapper->spline);
                delete[] d;
            }
            break;

        case ALGLIB_CATMULLROM:
            /* Catmull-Rom spline - smooth interpolation through all points */
            /* boundtype=0: periodic, boundtype=1: parabolic at ends */
            alglib::spline1dbuildcatmullrom(xa, ya, n, 1, 0.0, wrapper->spline);
            break;

        default:
            /* Default to cubic with natural boundary conditions */
            alglib::spline1dbuildcubic(xa, ya, wrapper->spline);
            break;
    }

    *spline_ptr = wrapper;
}

/*
 * Evaluate ALGLIB spline at a single point
 */
double alglib_spline_eval(void *spline_ptr, double x)
{
    alglib_spline_wrapper *wrapper = (alglib_spline_wrapper *)spline_ptr;
    return alglib::spline1dcalc(wrapper->spline, x);
}

/*
 * Free ALGLIB spline object
 */
void alglib_spline_free(void *spline_ptr)
{
    alglib_spline_wrapper *wrapper = (alglib_spline_wrapper *)spline_ptr;
    delete wrapper;
}

} /* extern "C" */