371 lines
12 KiB
C++
371 lines
12 KiB
C++
#pragma once
|
|
|
|
#include <expected>
|
|
#include <limits>
|
|
#include <memory>
|
|
#include <ranges>
|
|
#include <vector>
|
|
|
|
namespace mcc::fitpack
|
|
{
|
|
|
|
extern "C" {
|
|
void curfit(int* iopt,
|
|
int* m,
|
|
double* x,
|
|
double* y,
|
|
double* w,
|
|
double* xb,
|
|
double* xe,
|
|
int* k,
|
|
double* s,
|
|
int* nest,
|
|
int* n,
|
|
double* t,
|
|
double* c,
|
|
double* fp,
|
|
double* wrk,
|
|
int* lwrk,
|
|
int* iwrk,
|
|
int* ier);
|
|
|
|
void splev(double* t, int* n, double* c, int* k, double* x, double* y, int* m, int* e, int* ier);
|
|
|
|
void splder(double* t, int* n, double* c, int* k, int* nu, double* x, double* y, int* m, int* e, double* wrk, int* ier);
|
|
|
|
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);
|
|
}
|
|
|
|
|
|
|
|
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_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>,
|
|
"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!");
|
|
}
|
|
|
|
auto 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;
|
|
}
|
|
|
|
const auto nt = std::ranges::size(tetha_knots);
|
|
const auto 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;
|
|
|
|
const int lwrk1 = 185 + 52 * v + 10 * u + 14 * uv + 8 * (u - 1) * v2 + 8 * m;
|
|
const int lwrk2 = 48 + 21 * v + 7 * uv + 4 * (u - 1) * v2;
|
|
const int kwrk = m + uv;
|
|
|
|
std::unique_ptr<double> wrk1{new double(lwrk1)};
|
|
std::unique_ptr<double> wrk2{new double(lwrk2)};
|
|
std::unique_ptr<int> iwrk{new int(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);
|
|
}
|
|
|
|
const int iopt = 0;
|
|
|
|
if constexpr (std::ranges::contiguous_range<WeightT>) {
|
|
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.get(), &lwrk1,
|
|
wrk2.get(), &lwrk2, iwrk.get(), &kwrk, &res);
|
|
} else {
|
|
std::vector<double> weight_vec(m, weight);
|
|
|
|
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,
|
|
const TKnotT& tetha_knots,
|
|
const 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(tetha, phi, func, weight, s, nt, np, tetha_knots, phi_knots, coeffs, resi2_sum, eps);
|
|
}
|
|
|
|
/*
|
|
|
|
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,
|
|
const TKnotT& tetha_knots,
|
|
const 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>,
|
|
"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!");
|
|
}
|
|
|
|
auto 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;
|
|
|
|
const auto nt = std::ranges::size(tetha_knots);
|
|
const auto np = std::ranges::size(phi_knots);
|
|
|
|
// number of knots must be >= 8 for the qubic b-splines
|
|
if (std::min(nt, np) < 8) {
|
|
return 10;
|
|
}
|
|
|
|
if constexpr (std::ranges::contiguous_range<WeightT>) {
|
|
m = std::min(m, std::ranges::size(weight));
|
|
if (m < 2) {
|
|
return 10;
|
|
}
|
|
}
|
|
|
|
// compute working arrays sizes according to sphere.f
|
|
const int u = nt - 7;
|
|
const int v = np - 7;
|
|
const int uv = u * v;
|
|
const int v2 = v * v;
|
|
|
|
const int lwrk1 = 185 + 52 * v + 10 * u + 14 * uv + 8 * (u - 1) * v2 + 8 * m;
|
|
const int lwrk2 = 48 + 21 * v + 7 * uv + 4 * (u - 1) * v2;
|
|
const int kwrk = m + uv;
|
|
|
|
std::unique_ptr<double> wrk1{new double(lwrk1)};
|
|
std::unique_ptr<double> wrk2{new double(lwrk2)};
|
|
std::unique_ptr<int> iwrk{new int(kwrk)};
|
|
|
|
const int n_coeffs = (nt - 4) * (np - 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);
|
|
}
|
|
|
|
const int iopt = -1;
|
|
const double s = 100.0;
|
|
|
|
if constexpr (std::ranges::contiguous_range<WeightT>) {
|
|
res = sphere(&iopt, &m, std::ranges::data(tetha), std::ranges::data(phi), std::ranges::data(func),
|
|
std::ranges::data(weight), &s, &nt, &np, &eps, &nt, std::ranges::data(tetha_knots), &np,
|
|
std::ranges::data(phi_knots), std::ranges::data(coeffs), &resi2_sum, wrk1.get(), &lwrk1,
|
|
wrk2.get(), &lwrk2, iwrk.get(), &kwrk, &res);
|
|
} else {
|
|
std::vector<double> weight_vec(m, weight);
|
|
|
|
res = sphere(&iopt, &m, std::ranges::data(tetha), std::ranges::data(phi), std::ranges::data(func),
|
|
weight_vec.data(), &s, &nt, &np, &eps, &nt, std::ranges::data(tetha_knots), &np,
|
|
std::ranges::data(phi_knots), std::ranges::data(coeffs), &resi2_sum, wrk1.get(), &lwrk1,
|
|
wrk2.get(), &lwrk2, iwrk.get(), &kwrk, &res);
|
|
}
|
|
|
|
|
|
return res;
|
|
}
|
|
|
|
*/
|
|
|
|
} // namespace mcc::fitpack
|