#pragma once

#include <limits>
#include <ranges>
#include <vector>

#include <FortranCInterface.h>

extern "C" {
void surfit(int* iopt,
            int* m,
            double* x,
            double* y,
            double* z,
            double* w,
            double* xb,
            double* xe,
            double* yb,
            double* ye,
            int* kx,
            int* ky,
            double* s,
            int* nxest,
            int* nyest,
            int* nmax,
            double* eps,
            int* nx,
            double* tx,
            int* ny,
            double* ty,
            double* c,
            double* fp,
            double* wrk1,
            int* lwrk1,
            double* wrk2,
            int* lwrk2,
            int* iwrk,
            int* kwrk,
            int* ier);

void bispev(double* tx,
            int* nx,
            double* ty,
            int* ny,
            double* c,
            int* kx,
            int* ky,
            double* x,
            int* mx,
            double* y,
            int* my,
            double* z,
            double* wrk,
            int* lwrk,
            int* iwrk,
            int* kwrk,
            int* ier);

void parder(double* tx,
            int* nx,
            double* ty,
            int* ny,
            double* c,
            int* kx,
            int* ky,
            int* nux,
            int* nuy,
            double* x,
            int* mx,
            double* y,
            int* my,
            double* z,
            double* wrk,
            int* lwrk,
            int* iwrk,
            int* kwrk,
            int* ier);


void sphere(int* iopt,
            int* m,
            double* teta,
            double* phi,
            double* r,
            double* w,
            double* s,
            int* ntest,
            int* npest,
            double* eps,
            int* nt,
            double* tt,
            int* np,
            double* tp,
            double* c,
            double* fp,
            double* wrk1,
            int* lwrk1,
            double* wrk2,
            int* lwrk2,
            int* iwrk,
            int* kwrk,
            int* ier);
}



namespace mcc::fitpack
{

template <int IOPT = 0,
          std::ranges::contiguous_range TethaT,
          std::ranges::contiguous_range PhiT,
          std::ranges::contiguous_range FuncT,
          typename WeightT,
          std::ranges::contiguous_range TKnotT,
          std::ranges::contiguous_range PKnotT,
          std::ranges::contiguous_range CoeffT>
int fitpack_sphere_smooth(const TethaT& tetha,
                          const PhiT& phi,
                          const FuncT& func,
                          const WeightT& weight,
                          double s_par,
                          int& ntest,
                          int& npest,
                          TKnotT& tetha_knots,
                          PKnotT& phi_knots,
                          CoeffT& coeffs,
                          double& resi2_sum,
                          double eps = std::numeric_limits<double>::epsilon())
    requires((std::ranges::contiguous_range<WeightT> || std::convertible_to<WeightT, double>) &&
             std::ranges::output_range<CoeffT, double>)
{
    static_assert(std::same_as<std::ranges::range_value_t<TethaT>, double> &&
                      std::same_as<std::ranges::range_value_t<PhiT>, double> &&
                      std::same_as<std::ranges::range_value_t<FuncT>, double> &&
                      std::same_as<std::ranges::range_value_t<TKnotT>, double> &&
                      std::same_as<std::ranges::range_value_t<PKnotT>, double> &&
                      std::same_as<std::ranges::range_value_t<CoeffT>, double>,
                  "Input ranges elements type must be double!");

    if constexpr (std::ranges::contiguous_range<WeightT>) {
        static_assert(std::same_as<std::ranges::range_value_t<WeightT>, double>,
                      "Input ranges elements type must be double!");
    }


    static_assert(IOPT != -1 || IOPT != 0 || IOPT != 1, "Invalid IOPT template parameter value!");


    int m = std::min(std::min(std::ranges::size(tetha), std::ranges::size(phi)), std::ranges::size(func));
    if (m < 2) {
        return 10;
    }

    int res = 0;

    /*   do checking here to avoid possibly unnecessary memory allocations  */

    // number of knots must be >= 8 for the qubic b-splines
    if (ntest < 8 || npest < 8) {
        return 10;
    }

    if constexpr (std::ranges::contiguous_range<WeightT>) {
        m = std::min(m, std::ranges::size(weight));
        if (m < 2) {
            return 10;
        }
    }

    if (s_par <= 0.0) {
        return 10;
    }

    if (eps <= 0.0 || eps >= 1.0) {
        return 10;
    }

    int nt = std::ranges::size(tetha_knots);
    int np = std::ranges::size(phi_knots);

    if (nt < ntest) {
        std::ranges::fill_n(std::back_inserter(tetha_knots), ntest - nt, 0.0);
    }
    if (np < npest) {
        std::ranges::fill_n(std::back_inserter(phi_knots), npest - np, 0.0);
    }


    // compute working arrays sizes according to sphere.f
    const int u = ntest - 7;
    const int v = npest - 7;
    const int uv = u * v;
    const int v2 = v * v;

    int lwrk1 = 185 + 52 * v + 10 * u + 14 * uv + 8 * (u - 1) * v2 + 8 * m;
    int lwrk2 = 48 + 21 * v + 7 * uv + 4 * (u - 1) * v2;
    int kwrk = m + uv;

    std::vector<double> wrk1(lwrk1);
    std::vector<double> wrk2(lwrk2);
    std::vector<int> iwrk(kwrk);

    const int n_coeffs = (ntest - 4) * (npest - 4);
    if (std::ranges::size(coeffs) < n_coeffs) {  // resize
        std::ranges::fill_n(std::back_inserter(coeffs), n_coeffs - std::ranges::size(coeffs), 0.0);
    }

    int iopt = IOPT;

    auto tetha_ptr = const_cast<double*>(std::ranges::data(tetha));
    auto phi_ptr = const_cast<double*>(std::ranges::data(phi));
    auto func_ptr = const_cast<double*>(std::ranges::data(func));
    auto tetha_knots_ptr = const_cast<double*>(std::ranges::data(tetha_knots));
    auto phi_knots_ptr = const_cast<double*>(std::ranges::data(phi_knots));

    if constexpr (std::ranges::contiguous_range<WeightT>) {
        auto weight_ptr = const_cast<double*>(std::ranges::data(weight));

        sphere(&iopt, &m, tetha_ptr, phi_ptr, func_ptr, weight_ptr, &s_par, &nt, &np, &eps, &ntest, tetha_knots_ptr,
               &npest, phi_knots_ptr, std::ranges::data(coeffs), &resi2_sum, wrk1.data(), &lwrk1, wrk2.data(), &lwrk2,
               iwrk.data(), &kwrk, &res);
        // sphere(&iopt, &m, std::ranges::data(tetha), std::ranges::data(phi), std::ranges::data(func),
        //        std::ranges::data(weight), &s_par, &ntest, &npest, &eps, &ntest, std::ranges::data(tetha_knots),
        //        &npest, std::ranges::data(phi_knots), std::ranges::data(coeffs), &resi2_sum, wrk1.data(), &lwrk1,
        //        wrk2.data(), &lwrk2, iwrk.data(), &kwrk, &res);
    } else {
        std::vector<double> weight_vec(m, weight);

        sphere(&iopt, &m, tetha_ptr, phi_ptr, func_ptr, weight_vec.data(), &s_par, &nt, &np, &eps, &ntest,
               tetha_knots_ptr, &npest, phi_knots_ptr, std::ranges::data(coeffs), &resi2_sum, wrk1.data(), &lwrk1,
               wrk2.data(), &lwrk2, iwrk.data(), &kwrk, &res);
        // res = sphere(&iopt, &m, std::ranges::data(tetha), std::ranges::data(phi), std::ranges::data(func),
        //              weight_vec.data(), &s_par, &ntest, &npest, &eps, &ntest, std::ranges::data(tetha_knots), &npest,
        //              std::ranges::data(phi_knots), std::ranges::data(coeffs), &resi2_sum, wrk1.get(), &lwrk1,
        //              wrk2.get(), &lwrk2, iwrk.get(), &kwrk, &res);
    }


    return res;
}

/*
    least-squares fitting (iopt=-1)
*/
template <std::ranges::contiguous_range TethaT,
          std::ranges::contiguous_range PhiT,
          std::ranges::contiguous_range FuncT,
          typename WeightT,
          std::ranges::contiguous_range TKnotT,
          std::ranges::contiguous_range PKnotT,
          std::ranges::contiguous_range CoeffT>
int fitpack_sphere_fit(const TethaT& tetha,
                       const PhiT& phi,
                       const FuncT& func,
                       const WeightT& weight,
                       TKnotT& tetha_knots,
                       PKnotT& phi_knots,
                       CoeffT& coeffs,
                       double& resi2_sum,
                       double eps = std::numeric_limits<double>::epsilon())
{
    const double s = 100.0;
    int nt = std::ranges::size(tetha_knots);
    int np = std::ranges::size(phi_knots);

    return fitpack_sphere_smooth<-1>(tetha, phi, func, weight, s, nt, np, tetha_knots, phi_knots, coeffs, resi2_sum,
                                     eps);
}


template <std::ranges::contiguous_range TXT,
          std::ranges::contiguous_range TYT,
          std::ranges::contiguous_range XT,
          std::ranges::contiguous_range YT,
          std::ranges::contiguous_range CoeffT,
          std::ranges::contiguous_range FuncT>
int fitpack_eval_spl2d(const TXT& tx,
                       const TYT& ty,
                       const CoeffT& coeffs,
                       const XT& x,
                       const YT& y,
                       FuncT& func,
                       int kx = 3,
                       int ky = 3)
{
    static_assert(std::same_as<std::ranges::range_value_t<TXT>, double> &&
                      std::same_as<std::ranges::range_value_t<TYT>, double> &&
                      std::same_as<std::ranges::range_value_t<XT>, double> &&
                      std::same_as<std::ranges::range_value_t<YT>, double> &&
                      std::same_as<std::ranges::range_value_t<CoeffT>, double> &&
                      std::same_as<std::ranges::range_value_t<FuncT>, double>,
                  "Input ranges elements type must be double!");

    if (kx < 0 || ky < 0) {
        return 10;
    }


    int ier = 0;

    int ntx = std::ranges::size(tx);
    int nty = std::ranges::size(ty);

    auto n_coeffs = (ntx - kx - 1) * (nty - ky - 1);
    if (std::ranges::size(coeffs) < n_coeffs) {
        return 10;
    }

    int mx = std::ranges::size(x), my = std::ranges::size(y);
    int N = mx * my;

    if (std::ranges::size(func) < N) {
        std::ranges::fill_n(std::back_inserter(func), N - std::ranges::size(func), 0.0);
    }

    // compute sizes of working arrays according to bispev.f
    int lwrk = mx * (kx + 1) + my * (ky + 1);
    std::vector<double> wrk(lwrk);

    int kwrk = mx + my;
    std::vector<int> iwrk(kwrk);

    auto tx_ptr = const_cast<double*>(std::ranges::data(tx));
    auto ty_ptr = const_cast<double*>(std::ranges::data(ty));
    auto coeffs_ptr = const_cast<double*>(std::ranges::data(coeffs));
    auto x_ptr = const_cast<double*>(std::ranges::data(x));
    auto y_ptr = const_cast<double*>(std::ranges::data(y));

    bispev(tx_ptr, &ntx, ty_ptr, &nty, coeffs_ptr, &kx, &ky, x_ptr, &mx, y_ptr, &my, std::ranges::data(func),
           wrk.data(), &lwrk, iwrk.data(), &kwrk, &ier);

    return ier;
}


// scalar x,y version
template <std::ranges::contiguous_range TXT,
          std::ranges::contiguous_range TYT,
          typename XT,
          typename YT,
          std::ranges::contiguous_range CoeffT,
          typename FuncT>
int fitpack_eval_spl2d(const TXT& tx,
                       const TYT& ty,
                       const CoeffT& coeffs,
                       const XT& x,
                       const YT& y,
                       FuncT& func,
                       int kx = 3,
                       int ky = 3)
{
    static_assert(std::same_as<std::ranges::range_value_t<TXT>, double> &&
                      std::same_as<std::ranges::range_value_t<TYT>, double> &&
                      std::same_as<std::ranges::range_value_t<CoeffT>, double>,
                  "Input ranges elements type must be double!");

    static_assert(
        std::convertible_to<XT, double> && std::convertible_to<YT, double> && std::convertible_to<FuncT, double>,
        "XT, YT and FuncT types must be a scalar convertible to double!");

    auto xv = std::vector<double>(1, x);
    auto yv = std::vector<double>(1, y);
    auto fv = std::vector<double>(1, func);

    return fitpack_eval_spl2d(tx, ty, coeffs, xv, yv, fv, kx, ky);
}

}  // namespace mcc::fitpack