...

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()

include/mcc/mcc_pcm_construct.h

@@ -117,8 +117,11 @@ public:
         MccDefaultPCMType pcm_type;
         MccError error{};  // final model computation error
 
-        std::vector<double> model_colon, model_colat;  // fitting model values
+        size_t max_iter{100};  // max number of iterations for robust linear regression method
-        std::vector<double> colon_res, colat_res;      // target - model
+        size_t final_iter{};   // number of final iteration
+
+        std::vector<double> model_colon{}, model_colat{};  // fitting model values
+        std::vector<double> 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 ||
+        if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_BSPLINE) {
-            pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE) {
+            return bsplineFitting(pcm_data);
-            if (pcm_data.bspline.knotsX.size() < 2 || pcm_data.bspline.knotsY.size() < 2) {
+        } else if (pcm_data.type == MccDefaultPCMType::PCM_TYPE_GEOMETRY_BSPLINE) {
-                return MccDefaultPCMConstructorErrorCode::ERROR_INVALID_KNOTS_NUMBER;
+            // 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<double> xres(numberOfPoints()), yres(numberOfPoints());
-
+            for (size_t i = 0; i < numberOfPoints(); ++i) {
-            std::vector<double> tx(pcm_data.bspline.knotsX.size() + 6), ty(pcm_data.bspline.knotsY.size() + 6);
+                xres = _table.target_colon[i] - result.model_colon;
-
+                yres = _table.target_colat[i] - result.model_colat;
-            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<double> colon_res = _table.colon_res;
-                std::vector<double> 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<double> 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];
-                }
             }
         }
-    }
 #endif
+    }
 
 protected:
     // std::vector<table_elem_t> _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<double> 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<MOUNT_TYPE>::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