fpParametricCurve.hpp

Go to the documentation of this file.

/***************************************************************************************************
!                                ____________________  ___   ________ __
!                               / ____/  _/_  __/ __ \/   | / ____/ //_/
!                              / /_   / /  / / / /_/ / /| |/ /   / ,<
!                             / __/ _/ /  / / / ____/ ___ / /___/ /| |
!                            /_/   /___/ /_/ /_/   /_/  |_\____/_/ |_|
!
!                                     A Curve Fitting Package
!
!   fpParametricCurve.hpp (class fpParametricCurve)
!> @brief Standalone C++ wrapper for fitpack_parametric_curve (no fortran-arrays dependency)
!
!   @author Federico Perini
!   @date   2026-08-27
!
!   References :
!     - C. De Boor, "On calculating with b-splines", J Approx Theory 6 (1972) 50-62
!     - M. G. Cox, "The numerical evaluation of b-splines", J Inst Maths Applics 10 (1972) 134-149
!     - P. Dierckx, "Curve and surface fitting with splines", Monographs on numerical analysis,
!                    Oxford university press, 1993.
!
! **************************************************************************************************/
 
#ifndef FPPARAMETRICCURVE_HPP_INCLUDED
#define FPPARAMETRICCURVE_HPP_INCLUDED
 
#include "fitpack_config.h"
 
#if HAVE_FXARRAY
#include "fxArrays.hpp"
#endif
 
#include "fitpack_parametric_curve_c.h"
#include "fpFitter.hpp"
#include <string>
#include <vector>
#include <cstdint>
#include <memory>
#include <stdexcept>
#include <variant>
#include <optional>
#include <array>
#include "fpPoint.hpp"
 
static_assert(sizeof(fitpack_parametric_curve_c) == sizeof(fitpack_fitter_c),
    "C descriptor layout mismatch: fitpack_parametric_curve_c vs fitpack_fitter_c");

class fpParametricCurve : public fpFitter {
public:
    // ===========================================================================================
    // Constructors and Destructor
    // ===========================================================================================

    fpParametricCurve() : fpFitter(NoAlloc{}) {
        fitpack_parametric_curve_c_allocate(as<fitpack_parametric_curve_c>(), nullptr);
    }

    ~fpParametricCurve() override {
        fitpack_parametric_curve_c_destroy(as<fitpack_parametric_curve_c>(), nullptr);
    }

    fpParametricCurve(const fpParametricCurve& other) : fpFitter(NoAlloc{}) {
        *as<fitpack_parametric_curve_c>() = fitpack_parametric_curve_c_null;
        fitpack_parametric_curve_c_copy(as<fitpack_parametric_curve_c>(), other.as<fitpack_parametric_curve_c>(), false, nullptr);
    }

    fpParametricCurve& operator=(const fpParametricCurve& other) {
        if (this != &other) {
            fitpack_parametric_curve_c_destroy(as<fitpack_parametric_curve_c>(), nullptr);
            fitpack_parametric_curve_c_copy(as<fitpack_parametric_curve_c>(), other.as<fitpack_parametric_curve_c>(), false, nullptr);
        }
        return *this;
    }

    fpParametricCurve(fpParametricCurve&& other) noexcept : fpFitter(NoAlloc{}) {
        *as<fitpack_parametric_curve_c>() = fitpack_parametric_curve_c_null;
        fitpack_parametric_curve_c_move_alloc(as<fitpack_parametric_curve_c>(), other.as<fitpack_parametric_curve_c>(), nullptr);
    }

    fpParametricCurve& operator=(fpParametricCurve&& other) noexcept {
        if (this != &other) {
            fitpack_parametric_curve_c_destroy(as<fitpack_parametric_curve_c>(), nullptr);
            fitpack_parametric_curve_c_move_alloc(as<fitpack_parametric_curve_c>(), other.as<fitpack_parametric_curve_c>(), nullptr);
        }
        return *this;
    }

    explicit fpParametricCurve(fitpack_parametric_curve_c& c_wrapper, bool move = false) : fpFitter(NoAlloc{}) {
        if (move) {
            fitpack_parametric_curve_c_move_alloc(as<fitpack_parametric_curve_c>(), &c_wrapper, nullptr);
        } else {
            *as<fitpack_parametric_curve_c>() = fitpack_parametric_curve_c_null;
            fitpack_parametric_curve_c_copy(as<fitpack_parametric_curve_c>(), &c_wrapper, false, nullptr);
        }
    }

    explicit fpParametricCurve(fitpack_parametric_curve_c& c_wrapper, ViewTag) : fpFitter(NoAlloc{}) {
        *as<fitpack_parametric_curve_c>() = c_wrapper;
        as<fitpack_parametric_curve_c>()->is_pointer = true;
    }
 
    // ===========================================================================================
    // Utility Methods
    // ===========================================================================================

    bool is_allocated() const override {
        return as<fitpack_parametric_curve_c>()->cptr != nullptr;
    }

    bool is_pointer() const override {
        return as<fitpack_parametric_curve_c>()->is_pointer;
    }

    const char* c_type_name() const override {
        return fitpack_parametric_curve_c_c_type_name(*as<fitpack_parametric_curve_c>());
    }

    const char* cpp_type_name() const override {
        return "fpParametricCurve";
    }

    fitpack_parametric_curve_c& c_handle() {
        return *as<fitpack_parametric_curve_c>();
    }

    const fitpack_parametric_curve_c& c_handle() const {
        return *as<fitpack_parametric_curve_c>();
    }

    static fpParametricCurve make_view() {
        return fpParametricCurve(NoAlloc{});
    }

    fpFitter& as_parent() {
        return static_cast<fpFitter&>(*this);
    }
    const fpFitter& as_parent() const {
        return static_cast<const fpFitter&>(*this);
    }
 
    // ===========================================================================================
    // Method Wrappers (standalone — no fxArray dependency)
    // ===========================================================================================

    virtual void new_points(int32_t x_n1, int32_t x_n2, std::vector<double>& x, double* u = nullptr, double* w = nullptr) {
        fitpack_parametric_curve_c_new_points(as<fitpack_parametric_curve_c>(), x_n1, x_n2, x.data(), u, w);
    }
 
 
    virtual void set_default_parameters() {
        fitpack_parametric_curve_c_set_default_parameters(as<fitpack_parametric_curve_c>());
    }

    virtual int32_t new_fit(int32_t x_n1, int32_t x_n2, std::vector<double>& x, double* u = nullptr, double* w = nullptr, double* smoothing = nullptr, int32_t* order = nullptr) {
        return fitpack_parametric_curve_c_new_fit(as<fitpack_parametric_curve_c>(), x_n1, x_n2, x.data(), u, w, smoothing, order);
    }
 
 
    virtual int32_t fit(double* smoothing = nullptr, int32_t* order = nullptr, bool* keep_knots = nullptr) {
        return fitpack_parametric_curve_c_fit(as<fitpack_parametric_curve_c>(), smoothing, order, keep_knots);
    }
 
    virtual int32_t interpolate(int32_t* order = nullptr, bool* reset_knots = nullptr) {
        return fitpack_parametric_curve_c_interpolate(as<fitpack_parametric_curve_c>(), order, reset_knots);
    }
 
    virtual int32_t least_squares(double* smoothing = nullptr, bool* reset_knots = nullptr) {
        return fitpack_parametric_curve_c_least_squares(as<fitpack_parametric_curve_c>(), smoothing, reset_knots);
    }

    virtual std::vector<double> eval_one(double u, int32_t n_result, int32_t* ierr = nullptr) {
        std::vector<double> result(n_result);
        fitpack_parametric_curve_c_eval_one(as<fitpack_parametric_curve_c>(), u, ierr, result.data(), n_result);
        return result;
    }
 

    virtual std::vector<double> eval_many(std::vector<double>& u, int32_t n_result, int32_t* ierr = nullptr) {
        int32_t n = static_cast<int32_t>(u.size());
        std::vector<double> result(n_result);
        fitpack_parametric_curve_c_eval_many(as<fitpack_parametric_curve_c>(), n, u.data(), ierr, result.data(), n_result);
        return result;
    }
 

    virtual std::vector<double> dfdx(double u, int32_t order, int32_t n_result, int32_t* ierr = nullptr) {
        std::vector<double> result(n_result);
        fitpack_parametric_curve_c_curve_derivative(as<fitpack_parametric_curve_c>(), u, order, ierr, result.data(), n_result);
        return result;
    }
 

    virtual std::vector<double> dfdx(std::vector<double>& u, int32_t order, int32_t n_result, int32_t* ierr = nullptr) {
        int32_t n = static_cast<int32_t>(u.size());
        std::vector<double> result(n_result);
        fitpack_parametric_curve_c_curve_derivatives(as<fitpack_parametric_curve_c>(), n, u.data(), order, ierr, result.data(), n_result);
        return result;
    }
 

    virtual std::vector<double> dfdx_all(double u, int32_t n_result, int32_t* ierr = nullptr) {
        std::vector<double> result(n_result);
        fitpack_parametric_curve_c_curve_all_derivatives(as<fitpack_parametric_curve_c>(), u, ierr, result.data(), n_result);
        return result;
    }
 
 
    int32_t comm_size() const override {
        return fitpack_parametric_curve_c_comm_size(as<fitpack_parametric_curve_c>());
    }

    void comm_pack(std::vector<double>& buffer) const override {
        int32_t n = static_cast<int32_t>(buffer.size());
        fitpack_parametric_curve_c_comm_pack(as<fitpack_parametric_curve_c>(), n, buffer.data());
    }
 

    void comm_expand(std::vector<double>& buffer) override {
        int32_t n = static_cast<int32_t>(buffer.size());
        fitpack_parametric_curve_c_comm_expand(as<fitpack_parametric_curve_c>(), n, buffer.data());
    }
 
 
    double mse() const override {
        return fitpack_parametric_curve_c_mse(as<fitpack_parametric_curve_c>());
    }
 
    int32_t core_comm_size() const override {
        return fitpack_parametric_curve_c_core_comm_size(as<fitpack_parametric_curve_c>());
    }

    void core_comm_pack(std::vector<double>& buffer) const override {
        int32_t n = static_cast<int32_t>(buffer.size());
        fitpack_parametric_curve_c_core_comm_pack(as<fitpack_parametric_curve_c>(), n, buffer.data());
    }
 

    void core_comm_expand(std::vector<double>& buffer) override {
        int32_t n = static_cast<int32_t>(buffer.size());
        fitpack_parametric_curve_c_core_comm_expand(as<fitpack_parametric_curve_c>(), n, buffer.data());
    }
 
 
    void destroy_base() override {
        fitpack_parametric_curve_c_destroy_base(as<fitpack_parametric_curve_c>());
    }
 
    // ===========================================================================================
    // Component Array Accessors
    // ===========================================================================================

    std::vector<double> x_vector() const {
        double* raw = nullptr;
        int64_t extents[2] = {0, 0};
        fitpack_parametric_curve_c_getcomp_x(as<fitpack_parametric_curve_c>(), &raw, extents);
        const int64_t total = extents[0] * extents[1];
        if (raw == nullptr || total <= 0) return {};
        const double* first = reinterpret_cast<const double*>(raw);
        return std::vector<double>(first, first + total);
    }

    std::array<int64_t, 2> x_shape() const {
        double* raw = nullptr;
        int64_t extents[2] = {0, 0};
        fitpack_parametric_curve_c_getcomp_x(as<fitpack_parametric_curve_c>(), &raw, extents);
        return { extents[0], extents[1] };
    }
 
#if HAVE_FXARRAY
    fxArray<double> x() const {
        double* raw = nullptr;
        int64_t extents[2] = {0, 0};
        fitpack_parametric_curve_c_getcomp_x(as<fitpack_parametric_curve_c>(), &raw, extents);
        FX_SIZE bounds[4];   // (lower, upper) per dimension, Fortran lbound 1
        for (int k = 0; k < 2; ++k) {
            bounds[2 * k] = 1;
            bounds[2 * k + 1] = static_cast<FX_SIZE>(extents[k]);
        }
        array_c descr = array_c_null;
        array_c_from_ptr(&descr, "x", static_cast<void*>(raw),
                         getCFITypeFlag<double>(),
                         static_cast<FX_SIZE>(sizeof(double)),
                         static_cast<FX_RANK>(2), bounds);
        return fxArray<double>(descr);
    }
#endif // HAVE_FXARRAY

    std::vector<double> u_vector() const {
        double* raw = nullptr;
        int64_t extents[1] = {0};
        fitpack_parametric_curve_c_getcomp_u(as<fitpack_parametric_curve_c>(), &raw, extents);
        const int64_t total = extents[0];
        if (raw == nullptr || total <= 0) return {};
        const double* first = reinterpret_cast<const double*>(raw);
        return std::vector<double>(first, first + total);
    }
 
#if HAVE_FXARRAY
    fxArray<double> u() const {
        double* raw = nullptr;
        int64_t extents[1] = {0};
        fitpack_parametric_curve_c_getcomp_u(as<fitpack_parametric_curve_c>(), &raw, extents);
        FX_SIZE bounds[2];   // (lower, upper) per dimension, Fortran lbound 1
        for (int k = 0; k < 1; ++k) {
            bounds[2 * k] = 1;
            bounds[2 * k + 1] = static_cast<FX_SIZE>(extents[k]);
        }
        array_c descr = array_c_null;
        array_c_from_ptr(&descr, "u", static_cast<void*>(raw),
                         getCFITypeFlag<double>(),
                         static_cast<FX_SIZE>(sizeof(double)),
                         static_cast<FX_RANK>(1), bounds);
        return fxArray<double>(descr);
    }
#endif // HAVE_FXARRAY

    std::vector<double> sp_vector() const {
        double* raw = nullptr;
        int64_t extents[1] = {0};
        fitpack_parametric_curve_c_getcomp_sp(as<fitpack_parametric_curve_c>(), &raw, extents);
        const int64_t total = extents[0];
        if (raw == nullptr || total <= 0) return {};
        const double* first = reinterpret_cast<const double*>(raw);
        return std::vector<double>(first, first + total);
    }
 
#if HAVE_FXARRAY
    fxArray<double> sp() const {
        double* raw = nullptr;
        int64_t extents[1] = {0};
        fitpack_parametric_curve_c_getcomp_sp(as<fitpack_parametric_curve_c>(), &raw, extents);
        FX_SIZE bounds[2];   // (lower, upper) per dimension, Fortran lbound 1
        for (int k = 0; k < 1; ++k) {
            bounds[2 * k] = 1;
            bounds[2 * k + 1] = static_cast<FX_SIZE>(extents[k]);
        }
        array_c descr = array_c_null;
        array_c_from_ptr(&descr, "sp", static_cast<void*>(raw),
                         getCFITypeFlag<double>(),
                         static_cast<FX_SIZE>(sizeof(double)),
                         static_cast<FX_RANK>(1), bounds);
        return fxArray<double>(descr);
    }
#endif // HAVE_FXARRAY

    std::vector<double> w_vector() const {
        double* raw = nullptr;
        int64_t extents[1] = {0};
        fitpack_parametric_curve_c_getcomp_w(as<fitpack_parametric_curve_c>(), &raw, extents);
        const int64_t total = extents[0];
        if (raw == nullptr || total <= 0) return {};
        const double* first = reinterpret_cast<const double*>(raw);
        return std::vector<double>(first, first + total);
    }
 
#if HAVE_FXARRAY
    fxArray<double> w() const {
        double* raw = nullptr;
        int64_t extents[1] = {0};
        fitpack_parametric_curve_c_getcomp_w(as<fitpack_parametric_curve_c>(), &raw, extents);
        FX_SIZE bounds[2];   // (lower, upper) per dimension, Fortran lbound 1
        for (int k = 0; k < 1; ++k) {
            bounds[2 * k] = 1;
            bounds[2 * k + 1] = static_cast<FX_SIZE>(extents[k]);
        }
        array_c descr = array_c_null;
        array_c_from_ptr(&descr, "w", static_cast<void*>(raw),
                         getCFITypeFlag<double>(),
                         static_cast<FX_SIZE>(sizeof(double)),
                         static_cast<FX_RANK>(1), bounds);
        return fxArray<double>(descr);
    }
#endif // HAVE_FXARRAY

    std::vector<double> dd_vector() const {
        double* raw = nullptr;
        int64_t extents[2] = {0, 0};
        fitpack_parametric_curve_c_getcomp_dd(as<fitpack_parametric_curve_c>(), &raw, extents);
        const int64_t total = extents[0] * extents[1];
        if (raw == nullptr || total <= 0) return {};
        const double* first = reinterpret_cast<const double*>(raw);
        return std::vector<double>(first, first + total);
    }

    std::array<int64_t, 2> dd_shape() const {
        double* raw = nullptr;
        int64_t extents[2] = {0, 0};
        fitpack_parametric_curve_c_getcomp_dd(as<fitpack_parametric_curve_c>(), &raw, extents);
        return { extents[0], extents[1] };
    }
 
#if HAVE_FXARRAY
    fxArray<double> dd() const {
        double* raw = nullptr;
        int64_t extents[2] = {0, 0};
        fitpack_parametric_curve_c_getcomp_dd(as<fitpack_parametric_curve_c>(), &raw, extents);
        FX_SIZE bounds[4];   // (lower, upper) per dimension, Fortran lbound 1
        for (int k = 0; k < 2; ++k) {
            bounds[2 * k] = 1;
            bounds[2 * k + 1] = static_cast<FX_SIZE>(extents[k]);
        }
        array_c descr = array_c_null;
        array_c_from_ptr(&descr, "dd", static_cast<void*>(raw),
                         getCFITypeFlag<double>(),
                         static_cast<FX_SIZE>(sizeof(double)),
                         static_cast<FX_RANK>(2), bounds);
        return fxArray<double>(descr);
    }
#endif // HAVE_FXARRAY

    std::vector<double> t_vector() const {
        double* raw = nullptr;
        int64_t extents[1] = {0};
        fitpack_parametric_curve_c_getcomp_t(as<fitpack_parametric_curve_c>(), &raw, extents);
        const int64_t total = extents[0];
        if (raw == nullptr || total <= 0) return {};
        const double* first = reinterpret_cast<const double*>(raw);
        return std::vector<double>(first, first + total);
    }
 
#if HAVE_FXARRAY
    fxArray<double> t() const {
        double* raw = nullptr;
        int64_t extents[1] = {0};
        fitpack_parametric_curve_c_getcomp_t(as<fitpack_parametric_curve_c>(), &raw, extents);
        FX_SIZE bounds[2];   // (lower, upper) per dimension, Fortran lbound 1
        for (int k = 0; k < 1; ++k) {
            bounds[2 * k] = 1;
            bounds[2 * k + 1] = static_cast<FX_SIZE>(extents[k]);
        }
        array_c descr = array_c_null;
        array_c_from_ptr(&descr, "t", static_cast<void*>(raw),
                         getCFITypeFlag<double>(),
                         static_cast<FX_SIZE>(sizeof(double)),
                         static_cast<FX_RANK>(1), bounds);
        return fxArray<double>(descr);
    }
#endif // HAVE_FXARRAY
 
    // ===========================================================================================
    // Scalar Property Accessors
    // ===========================================================================================
 
    int32_t& m() {
        return *fitpack_parametric_curve_c_ref_m(as<fitpack_parametric_curve_c>());
    }
    const int32_t& m() const {
        return *fitpack_parametric_curve_c_ref_m(as<fitpack_parametric_curve_c>());
    }
 
    int32_t& idim() {
        return *fitpack_parametric_curve_c_ref_idim(as<fitpack_parametric_curve_c>());
    }
    const int32_t& idim() const {
        return *fitpack_parametric_curve_c_ref_idim(as<fitpack_parametric_curve_c>());
    }
 
    bool has_params() const {
        return fitpack_parametric_curve_c_get_has_params(as<fitpack_parametric_curve_c>());
    }
    void set_has_params(bool value) {
        fitpack_parametric_curve_c_set_has_params(as<fitpack_parametric_curve_c>(), value);
    }
 
    int32_t& order() {
        return *fitpack_parametric_curve_c_ref_order(as<fitpack_parametric_curve_c>());
    }
    const int32_t& order() const {
        return *fitpack_parametric_curve_c_ref_order(as<fitpack_parametric_curve_c>());
    }
 
    double& ubegin() {
        return *fitpack_parametric_curve_c_ref_ubegin(as<fitpack_parametric_curve_c>());
    }
    const double& ubegin() const {
        return *fitpack_parametric_curve_c_ref_ubegin(as<fitpack_parametric_curve_c>());
    }
 
    double& uend() {
        return *fitpack_parametric_curve_c_ref_uend(as<fitpack_parametric_curve_c>());
    }
    const double& uend() const {
        return *fitpack_parametric_curve_c_ref_uend(as<fitpack_parametric_curve_c>());
    }
 
    int32_t& nest() {
        return *fitpack_parametric_curve_c_ref_nest(as<fitpack_parametric_curve_c>());
    }
    const int32_t& nest() const {
        return *fitpack_parametric_curve_c_ref_nest(as<fitpack_parametric_curve_c>());
    }
 
    int32_t& knots() {
        return *fitpack_parametric_curve_c_ref_knots(as<fitpack_parametric_curve_c>());
    }
    const int32_t& knots() const {
        return *fitpack_parametric_curve_c_ref_knots(as<fitpack_parametric_curve_c>());
    }
 
    // ===========================================================================================
    // fpPoint<dim> overloads — extra_methods/fpParametricPoints.hpp (hand-maintained)
    //
    // The dimension is a template argument, so a point is a fixed-size POD rather than a heap
    // allocation, and a std::vector<fpPoint<dim>> is bit-identical to the Fortran x(dim,m) it
    // feeds. Generalizes the idea of PR #61 (@illionj) to the whole parametric family.
    //
    // Everything here forwards to a virtual raw method, so fpClosedCurve and fpConstrainedCurve
    // inherit these overloads and still reach their own Fortran routines through dispatch.
    // Methods that read an existing fit check idim() against dim and report
    // FITPACK_INPUT_ERROR on a mismatch rather than reading past the result buffer.
    // ===========================================================================================

    FP_SIZE degree() const { return order(); }

    FP_SIZE ndim() const { return idim(); }

    template <FP_SIZE dim>
    FP_FLAG new_fit(const std::vector<fpPoint<dim>>& x, FP_REAL smoothing = 1000.0, FP_SIZE order = 3)
    {
        std::vector<FP_REAL> xw = fpPointFlatten<dim>(x);
        return new_fit(dim, static_cast<FP_SIZE>(x.size()), xw, nullptr, nullptr, &smoothing, &order);
    }

    template <FP_SIZE dim>
    FP_FLAG new_fit(const std::vector<fpPoint<dim>>& x, const std::vector<FP_REAL>& u,
                    FP_REAL smoothing = 1000.0, FP_SIZE order = 3)
    {
        std::vector<FP_REAL> xw = fpPointFlatten<dim>(x), uw(u);
        return new_fit(dim, static_cast<FP_SIZE>(x.size()), xw, uw.data(), nullptr, &smoothing, &order);
    }

    template <FP_SIZE dim>
    FP_FLAG new_fit(const std::vector<fpPoint<dim>>& x, const std::vector<FP_REAL>& u,
                    const std::vector<FP_REAL>& w, FP_REAL smoothing = 1000.0, FP_SIZE order = 3)
    {
        std::vector<FP_REAL> xw = fpPointFlatten<dim>(x), uw(u), ww(w);
        return new_fit(dim, static_cast<FP_SIZE>(x.size()), xw, uw.data(), ww.data(), &smoothing, &order);
    }

    FP_FLAG fit(FP_SIZE order)                    { return fit(nullptr, &order, nullptr); }
    FP_FLAG fit(FP_REAL smoothing)                { return fit(&smoothing, nullptr, nullptr); }
    FP_FLAG fit(FP_REAL smoothing, FP_SIZE order) { return fit(&smoothing, &order, nullptr); }

    template <FP_SIZE dim>
    fpPoint<dim> eval(FP_REAL u, FP_FLAG* ierr = nullptr)
    {
        fpPoint<dim> y{};
        if (idim() != dim) { if (ierr) *ierr = FITPACK_INPUT_ERROR; return y; }
 
        FP_FLAG ierr0 = FITPACK_OK;
        const std::vector<FP_REAL> flat = eval_one(u, dim, &ierr0);
        if (ierr) *ierr = ierr0;
        for (FP_SIZE j = 0; j < dim; ++j) y[j] = flat[static_cast<std::size_t>(j)];
        return y;
    }

    template <FP_SIZE dim>
    std::vector<fpPoint<dim>> eval(const std::vector<FP_REAL>& u, FP_FLAG* ierr = nullptr)
    {
        if (idim() != dim) { if (ierr) *ierr = FITPACK_INPUT_ERROR; return {}; }
 
        std::vector<FP_REAL> uw(u);
        FP_FLAG ierr0 = FITPACK_OK;
        const std::vector<FP_REAL> flat = eval_many(uw, dim * static_cast<FP_SIZE>(u.size()), &ierr0);
        if (ierr) *ierr = ierr0;
        return fpPointGather<dim>(flat);
    }

    template <FP_SIZE dim>
    fpPoint<dim> ddu(FP_REAL u, FP_SIZE order, FP_FLAG* ierr = nullptr)
    {
        fpPoint<dim> d{};
        if (idim() != dim) { if (ierr) *ierr = FITPACK_INPUT_ERROR; return d; }
 
        FP_FLAG ierr0 = FITPACK_OK;
        const std::vector<FP_REAL> flat = dfdx(u, order, dim, &ierr0);
        if (ierr) *ierr = ierr0;
        for (FP_SIZE j = 0; j < dim; ++j) d[j] = flat[static_cast<std::size_t>(j)];
        return d;
    }

    template <FP_SIZE dim>
    std::vector<fpPoint<dim>> ddu_all(FP_REAL u, FP_FLAG* ierr = nullptr)
    {
        if (idim() != dim) { if (ierr) *ierr = FITPACK_INPUT_ERROR; return {}; }
 
        FP_FLAG ierr0 = FITPACK_OK;
        const std::vector<FP_REAL> flat = dfdx_all(u, dim * (degree() + 1), &ierr0);
        if (ierr) *ierr = ierr0;
        return fpPointGather<dim>(flat);
    }
 
protected:
    explicit fpParametricCurve(NoAlloc tag) : fpFitter(tag) {}
 
};
 
#endif /* FPPARAMETRICCURVE_HPP_INCLUDED */