diff --git a/CMakeLists.txt b/CMakeLists.txt index bcfa1dc..ecc8063 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,6 +27,12 @@ option( ON ) +option( + USE_EIGEN3 + "Use of Eigen3 library as dependency for implemetation of PCM construction" + ON +) + option(USE_BSPLINE_PCM "Use of FITPACK bivariate splines for PCM" ON) option( @@ -42,6 +48,27 @@ find_package(Threads REQUIRED) include(FetchContent) include(ExternalProject) +# ******** Eigen3 LIBRARY ******* + +if(USE_EIGEN3) + find_package(Eigen3 CONFIG) + if(TARGET Eigen3::Eigen) + message(STATUS "Eigen3 library was found!") + else() + message(STATUS \tfetch Eigen3 ...) + FetchContent_Declare( + eigen + GIT_REPOSITORY "https://gitlab.com/libeigen/eigen.git" + GIT_SHALLOW TRUE + GIT_PROGRESS TRUE + OVERRIDE_FIND_PACKAGE + ) + FetchContent_MakeAvailable(eigen) + + find_package(Eigen3 REQUIRED CONFIG) + endif() +endif() + # ******* SPDLOG LIBRARY ******* if(USE_SPDLOG) @@ -250,6 +277,10 @@ if(USE_BSPLINE_PCM) add_subdirectory(fitpack) endif() +if(USE_EIGEN3) + list(APPEND MCC_SRC include/mcc/mcc_pcm_construct.h) +endif() + add_library(${PROJECT_NAME} INTERFACE ${MCC_SRC}) target_compile_features(${PROJECT_NAME} INTERFACE cxx_std_23) target_include_directories( @@ -272,6 +303,10 @@ if(USE_BSPLINE_PCM) ) endif() +if(USE_EIGEN3) + target_link_libraries(${PROJECT_NAME} INTERFACE Eigen3::Eigen) +endif() + if(USE_ASIO) target_link_libraries(${PROJECT_NAME} INTERFACE PkgConfig::ASIOLIB) endif() diff --git a/include/mcc/mcc_pcm_construct.h b/include/mcc/mcc_pcm_construct.h index 645df3d..1697822 100644 --- a/include/mcc/mcc_pcm_construct.h +++ b/include/mcc/mcc_pcm_construct.h @@ -117,8 +117,11 @@ public: MccDefaultPCMType pcm_type; MccError error{}; // final model computation error - std::vector model_colon, model_colat; // fitting model values - std::vector colon_res, colat_res; // target - model + size_t max_iter{100}; // max number of iterations for robust linear regression method + size_t final_iter{}; // number of final iteration + + std::vector model_colon{}, model_colat{}; // fitting model values + std::vector colon_res{}, colat_res{}; // target - model #ifdef USE_BSPLINE_PCM int bspline_fit_err{}; // bivariate B-spline fitting exit code (see FITPACK) @@ -229,93 +232,29 @@ public: } // robust linear regression with Tukey's loss function + } else { + if (_table.size() < min_data_size) { + result.error = MccDefaultPCMConstructorErrorCode::ERROR_NOT_ENOUGH_DATA; + return result; + } } #ifdef USE_BSPLINE_PCM - if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE || - pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE) { - if (pcm_data.bspline.knotsX.size() < 2 || pcm_data.bspline.knotsY.size() < 2) { - return MccDefaultPCMConstructorErrorCode::ERROR_INVALID_KNOTS_NUMBER; - } + if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) { + return bsplineFitting(pcm_data); + } else if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE) { + // the fitting for geometrical coefficients is already done above so + // one must fit residuals by bivariate B-splines - double resi2x, resi2y; // fitting residuals - - std::vector tx(pcm_data.bspline.knotsX.size() + 6), ty(pcm_data.bspline.knotsY.size() + 6); - - size_t Ncoeffs = (tx.size() - 4) * (ty.size() - 4); - - pcm_data.bspline.coeffsX.resize(Ncoeffs); - pcm_data.bspline.coeffsY.resize(Ncoeffs); - - if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) { - if (_table.size() < min_data_size) { - result.error = MccDefaultPCMConstructorErrorCode::ERROR_NOT_ENOUGH_DATA; - return result; - } - - // here both direct and inverse coefficients will be calculated - pcm_data.inverseBspline.coeffsX.resize(Ncoeffs); - pcm_data.inverseBspline.coeffsY.resize(Ncoeffs); - - // direct (celestial = encoder + pcm) - result.bspline_fit_err = bsplines::fitpack_sphere_fit( - _table.hw_colat, _table.hw_colon, _table.colon_res, 1.0, pcm_data.bspline.knotsY, - pcm_data.bspline.knotsX, pcm_data.bspline.coeffsX, resi2x); - if (result.bspline_fit_err > 0) { - result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; - return result; - } - - result.bspline_fit_err = bsplines::fitpack_sphere_fit( - _table.hw_colat, _table.hw_colon, _table.colat_res, 1.0, pcm_data.bspline.knotsY, - pcm_data.bspline.knotsX, pcm_data.bspline.coeffsY, resi2y); - if (result.bspline_fit_err > 0) { - result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; - return result; - } - - - // inverse (encoder = celestial + pcm) - std::vector colon_res = _table.colon_res; - std::vector colat_res = _table.colat_res; - for (size_t i = 0; i < colat_res.size(); ++i) { - colon_res[i] = -colon_res[i]; - colat_res[i] = -colat_res[i]; - } - - result.bspline_fit_err = bsplines::fitpack_sphere_fit( - _table.target_colon, _table.target_colat, colon_res, 1.0, pcm_data.bspline.knotsY, - pcm_data.bspline.knotsX, pcm_data.bspline.inverseCoeffsX, resi2x); - if (result.bspline_fit_err > 0) { - result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; - return result; - } - - result.bspline_fit_err = bsplines::fitpack_sphere_fit( - _table.target_colon, _table.target_colat, colat_res, 1.0, pcm_data.bspline.knotsY, - pcm_data.bspline.knotsX, pcm_data.bspline.inverseCoeffsY, resi2y); - if (result.bspline_fit_err > 0) { - result.error = MccDefaultPCMConstructorErrorCode::ERROR_BSPLINE_FIT; - return result; - } - } else { // geometry + B-spline - // the fitting for geometrical coefficients is already done above so - // one must fit residuals by bivariate B-splines - - std::vector xres(numberOfPoints()), yres(numberOfPoints()); - // for (size_t i = 0; i < _table.size(); ++i) { - // xres = _table[i].target_colon; - // yres = _table[i].target_colat; - // } - for (size_t i = 0; i < numberOfPoints(); ++i) { - xres = _table.target_colon[i] - result.; - yres = _table.target_colat[i]; - } + std::vector xres(numberOfPoints()), yres(numberOfPoints()); + for (size_t i = 0; i < numberOfPoints(); ++i) { + xres = _table.target_colon[i] - result.model_colon; + yres = _table.target_colat[i] - result.model_colat; } } - } #endif + } protected: // std::vector _table; @@ -327,10 +266,12 @@ protected: compute_result_t result{.pcm_type = pcm_data.type, .error = MccDefaultPCMConstructorErrorCode::ERROR_OK}; if (pcm_data.bspline.knotsX.size() < 2 || pcm_data.bspline.knotsY.size() < 2) { - return MccDefaultPCMConstructorErrorCode::ERROR_INVALID_KNOTS_NUMBER; + result.error = MccDefaultPCMConstructorErrorCode::ERROR_INVALID_KNOTS_NUMBER; + + return result; } - double resi2x, resi2y; // fitting residuals + double resi2x, resi2y; // sum of fitting residuals squares std::vector tx(pcm_data.bspline.knotsX.size() + 6), ty(pcm_data.bspline.knotsY.size() + 6); @@ -432,7 +373,85 @@ protected: } - void robustLinearRegress() {} + compute_result_t robustLinearRegress(MccDefaultPCM::pcm_data_t& pcm_data, size_t max_iter = 100) + { + auto tukey_weight = [](Eigen::VectorXd const& resi, double tukey_c, Eigen::VectorXd& tukey_weight) { + Eigen::ArrayXd u = (resi / tukey_c).cwiseAbs(); + auto mask = (u < 1.0); + + Eigen::VectorXd r = mask.select(1.0 - u * u, 0.0); + + tukey_weight = r * r; + }; + + auto gastwirth_est = [](Eigen::VectorXd const& resi) { + // Gastwirth’s estimation: 0.3*Q(1/3) + 0.4*Q(1/2) + 0.3*Q(2/3) + + Eigen::ArrayXd arr = resi.cwiseAbs(); + + std::nth_element(arr.begin(), arr.begin() + 2 * arr.size() / 3, arr.end()); + + // add 1.0E-8 to avoid 0.0 + return (0.3 * arr(arr.size() / 3) + 0.4 * arr(arr.size() / 2) + 0.3 * arr(2 * arr.size() / 3)) / 0.6745 + + 1.0E-8; + }; + + compute_result_t result{ + .pcm_type = pcm_data.type, .error = MccDefaultPCMConstructorErrorCode::ERROR_OK, .max_iter = max_iter}; + + const size_t n_coeffs = pcm_data.type == MOUNT_TYPE == MccMountType::FORK_TYPE ? 9 : 8; + + Eigen::MatrixXd X(2 * numberOfPoints(), n_coeffs); + Eigen::VectorXd y(2 * numberOfPoints()), beta(n_coeffs); + + Eigen::RowVectorXd row(n_coeffs); + double tgY, cosY, rcosY, cosX, sinX, cosPhi = std::cos(pcm_data.siteLatitude), + sinPhi = std::sin(pcm_data.siteLatitude); + + // system matrix and right-hand vector + for (size_t i = 0; i < numberOfPoints(); ++i) { + cosY = std::cos(_table.target_colat[i]); + rcosY = 1.0 / cosY; + tgY = std::tan(_table.target_colat[i]); + cosX = std::cos(_table.target_colon[i]); + sinX = std::sin(_table.target_colon[i]); + + + // [X0, Y0, collim, nonperp, misal1, misal2, tubefl, DEcaxfl, forkfl] or + // [X0, Y0, collim, nonperp, misal1, misal2, tubefl, DEcaxfl] + if constexpr (MOUNT_TYPE == MccMountType::FORK_TYPE) { + row << 1.0, 0.0, rcosY, tgY, -cosX * tgY, sinX * tgY, cosPhi * sinX * rcosY, + -(cosPhi * cosX + sinPhi * tgY), 0.0; + } else { + row << 1.0, 0.0, rcosY, tgY, -cosX * tgY, sinX * tgY, cosPhi * sinX * rcosY, + -(cosPhi * cosX + sinPhi * tgY); + } + + X.row(i) = row; + + if constexpr (MOUNT_TYPE == MccMountType::FORK_TYPE) { + row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * std::sin(_table.target_colat[i]) - sinPhi * cosY, + 0.0, 1.0 / cosX; + } else { + row << 0.0, 1.0, 0.0, 0.0, sinX, cosX, cosPhi * cosX * std::sin(_table.target_colat[i]) - sinPhi * cosY, + 0.0; + } + + X.row(i + numberOfPoints()) = row; + + y(i) = _table.colon_res[i]; + y(i + numberOfPoints()) = _table.colat_res[i]; + } + + Eigen::VectorXd resi = y; + Eigen::VectorXd weights(2 * numberOfPoints()); + + double scale; + + for (result.final_iter = 1; result.final_iter <= max_iter; ++result.final_iter) { + scale = gastwirth_est(resi); + } + } }; } // namespace mcc::impl \ No newline at end of file