From cfe1fc8fab88cde4a521cbed6213071dce9a2436 Mon Sep 17 00:00:00 2001 From: "Timur A. Fatkhullin" Date: Wed, 7 May 2025 00:01:14 +0300 Subject: [PATCH] ... --- cxx/CMakeLists.txt | 5 + cxx/fitpack/CMakeLists.txt | 9 +- cxx/fitpack/fitpack.h | 215 ++++++++++++++++--------------------- cxx/tests/fitpack_test.cpp | 108 +++++++++++++++++++ 4 files changed, 209 insertions(+), 128 deletions(-) create mode 100644 cxx/tests/fitpack_test.cpp diff --git a/cxx/CMakeLists.txt b/cxx/CMakeLists.txt index 1e9f4e6..4a57beb 100644 --- a/cxx/CMakeLists.txt +++ b/cxx/CMakeLists.txt @@ -142,4 +142,9 @@ if (WITH_TESTS) set(ASTROM_TEST_APP astrom_test) add_executable(${ASTROM_TEST_APP} tests/astrom_test.cpp) target_link_libraries(${ASTROM_TEST_APP} ERFA_LIB) + + set(FITPACK_TEST_APP fitpack_test) + add_executable(${FITPACK_TEST_APP} tests/fitpack_test.cpp) + target_link_libraries(${FITPACK_TEST_APP} fitpack) + # target_include_directories(${FITPACK_TEST_APP} PRIVATE fitpack) endif() diff --git a/cxx/fitpack/CMakeLists.txt b/cxx/fitpack/CMakeLists.txt index 5820206..af78251 100644 --- a/cxx/fitpack/CMakeLists.txt +++ b/cxx/fitpack/CMakeLists.txt @@ -20,11 +20,12 @@ enable_language(Fortran CXX) include(FortranCInterface) FortranCInterface_HEADER(FortranCInterface.h MACRO_NAMESPACE "FC_" - SYMBOL_NAMESPACE "fp_" - SYMBOLS ${func_str} + # SYMBOL_NAMESPACE "fp_" + SYMBOL_NAMESPACE "" + # SYMBOLS ${func_str} + SYMBOLS ${func_name} ) +FortranCInterface_VERIFY(CXX) add_library(fitpack STATIC ${src_files} fitpack.h) -set(FITPACK_INCLUDE_DIR ${CMAKE_CURRENT_BINARY_DIR}) - diff --git a/cxx/fitpack/fitpack.h b/cxx/fitpack/fitpack.h index cf62cde..3d52618 100644 --- a/cxx/fitpack/fitpack.h +++ b/cxx/fitpack/fitpack.h @@ -1,38 +1,13 @@ #pragma once -#include #include #include #include #include -namespace mcc::fitpack -{ +#include "/home/timur/PROGRAMS/C++/build-mountcontrol-Desktop-Debug/cxx/fitpack/FortranCInterface.h" 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, @@ -130,7 +105,11 @@ void sphere(int* iopt, -template , double> && std::same_as, double> && std::same_as, double> && - std::same_as, double>, + std::same_as, double> && + std::same_as, double>, "Input ranges elements type must be double!"); if constexpr (std::ranges::contiguous_range) { @@ -164,7 +144,11 @@ int fitpack_sphere_smooth(const TethaT& tetha, "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)); + + 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; } @@ -210,33 +194,47 @@ int fitpack_sphere_smooth(const TethaT& tetha, 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; + 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::unique_ptr wrk1{new double(lwrk1)}; - std::unique_ptr wrk2{new double(lwrk2)}; - std::unique_ptr iwrk{new int(kwrk)}; + std::vector wrk1(lwrk1); + std::vector wrk2(lwrk2); + std::vector 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); } - const int iopt = 0; + int iopt = IOPT; + + auto tetha_ptr = const_cast(std::ranges::data(tetha)); + auto phi_ptr = const_cast(std::ranges::data(phi)); + auto func_ptr = const_cast(std::ranges::data(func)); + auto tetha_knots_ptr = const_cast(std::ranges::data(tetha_knots)); + auto phi_knots_ptr = const_cast(std::ranges::data(phi_knots)); if constexpr (std::ranges::contiguous_range) { - 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); + auto weight_ptr = const_cast(std::ranges::data(weight)); + + sphere(&iopt, &m, tetha_ptr, phi_ptr, func_ptr, weight_ptr, &s_par, &ntest, &npest, &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 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); + sphere(&iopt, &m, tetha_ptr, phi_ptr, func_ptr, weight_vec.data(), &s_par, &ntest, &npest, &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); } @@ -257,8 +255,8 @@ int fitpack_sphere_fit(const TethaT& tetha, const PhiT& phi, const FuncT& func, const WeightT& weight, - const TKnotT& tetha_knots, - const PKnotT& phi_knots, + TKnotT& tetha_knots, + PKnotT& phi_knots, CoeffT& coeffs, double& resi2_sum, double eps = std::numeric_limits::epsilon()) @@ -267,104 +265,73 @@ int fitpack_sphere_fit(const TethaT& tetha, 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); + return fitpack_sphere_smooth<-1>(tetha, phi, func, weight, s, nt, np, tetha_knots, phi_knots, coeffs, resi2_sum, + eps); } -/* -template -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::epsilon()) - requires((std::ranges::contiguous_range || std::convertible_to) && - std::ranges::output_range) +template +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, double> && - std::same_as, double> && - std::same_as, double> && - std::same_as, double> && - std::same_as, double>, + static_assert(std::same_as, double> && + std::same_as, double> && + std::same_as, double> && + std::same_as, double> && + std::same_as, double> && + std::same_as, double>, "Input ranges elements type must be double!"); - if constexpr (std::ranges::contiguous_range) { - static_assert(std::same_as, 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) { + if (kx < 0 || ky < 0) { return 10; } - int res = 0; - const auto nt = std::ranges::size(tetha_knots); - const auto np = std::ranges::size(phi_knots); + int ier = 0; - // number of knots must be >= 8 for the qubic b-splines - if (std::min(nt, np) < 8) { + 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; } - if constexpr (std::ranges::contiguous_range) { - m = std::min(m, std::ranges::size(weight)); - if (m < 2) { - 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 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; + // compute sizes of working arrays according to bispev.f + int lwrk = mx * (kx + 1) + my * (ky + 1); + std::vector wrk(lwrk); - 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; + int kwrk = mx + my; + std::vector iwrk(kwrk); - std::unique_ptr wrk1{new double(lwrk1)}; - std::unique_ptr wrk2{new double(lwrk2)}; - std::unique_ptr iwrk{new int(kwrk)}; + auto tx_ptr = const_cast(std::ranges::data(tx)); + auto ty_ptr = const_cast(std::ranges::data(ty)); + auto coeffs_ptr = const_cast(std::ranges::data(coeffs)); + auto x_ptr = const_cast(std::ranges::data(x)); + auto y_ptr = const_cast(std::ranges::data(y)); - 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); - } + 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); - const int iopt = -1; - const double s = 100.0; - - if constexpr (std::ranges::contiguous_range) { - 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 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; + return ier; } -*/ - } // namespace mcc::fitpack diff --git a/cxx/tests/fitpack_test.cpp b/cxx/tests/fitpack_test.cpp new file mode 100644 index 0000000..af5ea69 --- /dev/null +++ b/cxx/tests/fitpack_test.cpp @@ -0,0 +1,108 @@ +#include +#include +#include +#include +#include + +// #include +#include "../fitpack/fitpack.h" + +int main() +{ + // size_t nt = 3, np = 5, N = nt * np, i = 1; + size_t nt = 10, np = 20, N = nt * np, i = 1; + double ts = std::numbers::pi / (nt + 1); + double ps = 2.0 * std::numbers::pi / (np + 1); + + std::vector tetha(N), phi(N), func(N); + + auto gen_func = [](double st, size_t& idx) { + double v = st * idx; + ++idx; + + return v; + }; + + auto print_func = [](const auto& r, std::string_view name) { + std::cout << name << ": "; + for (auto& el : r) { + std::cout << el << " "; + } + std::cout << "\n"; + }; + + + // std::ranges::generate(tetha, std::bind(gen_func, ts, i)); + + // i = 1; + // std::ranges::generate(phi, std::bind(gen_func, ps, i)); + + size_t k = 1; + i = 0; + for (size_t j = 0; j < nt; ++j) { + std::ranges::fill_n(tetha.begin() + i * np, np, ts * (i + 1)); + std::ranges::generate(phi | std::views::drop(i * np) | std::views::take(np), std::bind(gen_func, ps, k)); + ++i; + k = 1; + } + + std::uniform_real_distribution distr{-0.1, 0.1}; + std::random_device device; + std::mt19937 engine{device()}; + + std::ranges::generate(func, [ii = 0, &distr, &engine, &tetha, &phi]() mutable { + double v = (5.0 + tetha[ii]) * 1.3 + (3.0 + phi[ii]) * 3.1 + distr(engine); + ++ii; + return v; + }); + + + size_t ntk = 3, npk = 6, nf = (ntk + 4) * (npk + 4); + std::vector tk(ntk + 8), pk(npk + 8), coeffs(nf); + ts = std::numbers::pi / (ntk + 1); + ps = 2.0 * std::numbers::pi / (npk + 1); + + i = 1; + std::ranges::generate(tk | std::views::drop(4) | std::views::take(ntk), std::bind(gen_func, ts, i)); + i = 1; + std::ranges::generate(pk | std::views::drop(4) | std::views::take(npk), std::bind(gen_func, ps, i)); + + + double rs = 0.0; + int ec = mcc::fitpack::fitpack_sphere_fit(tetha, phi, func, 1.0, tk, pk, coeffs, rs); + + std::cout << "FIT EC = " << ec << "\n"; + std::cout << "FIT RESI = " << rs << "\n"; + + print_func(coeffs, "coeffs"); + + // print_func(tetha, "tetha"); + // print_func(phi, "phi"); + print_func(tk, "tetha_knots"); + print_func(func, "func"); + + std::cout << "\n\n"; + + k = 1; + ts = std::numbers::pi / (nt + 1); + std::ranges::generate_n(tetha.begin(), nt, std::bind(gen_func, ts, k)); + + std::vector f_func; + tetha.resize(nt); + phi.resize(np); + + ec = mcc::fitpack::fitpack_eval_spl2d(tk, pk, coeffs, tetha, phi, f_func); + // ec = mcc::fitpack::fitpack_eval_spl2d(pk, tk, coeffs, phi, tetha, f_func); + std::cout << "EVAL EC = " << ec << "\n"; + print_func(f_func, "func"); + + std::cout << "\n\n"; + + for (size_t l = 0; l < f_func.size(); ++l) { + auto r = f_func[l] - func[l]; + std::cout << r << " "; + } + std::cout << "\n"; + + return 0; +}