Add integration test for eigen engine and more documentation

cmake/load_dependencies.cmake

@@ -3,6 +3,8 @@ find_package(magic_enum REQUIRED CONFIG)
 find_package(CLI11 REQUIRED CONFIG)
 find_package(Eigen3 REQUIRED CONFIG)
 find_package(GSL REQUIRED)
+find_package(libassert)
+find_package(Boost)
 
 if(ENABLE_TEST)
     find_package(GTest CONFIG REQUIRED)

conanfile.py

@@ -8,13 +8,14 @@
 
 class CompressorRecipe(ConanFile):
     settings = "os", "compiler", "build_type", "arch"
-    generators = "CMakeDeps"
+    generators = "CMakeConfigDeps"
 
     def requirements(self):
         self.requires("spdlog/1.16.0", options={"use_std_fmt": True, "no_exceptions": True})  # type: ignore
         self.requires("magic_enum/0.9.7")   # type: ignore
         self.requires("cli11/2.6.0")        # type: ignore
         self.requires("eigen/5.0.1")        # type: ignore
+        self.requires("libassert/2.2.1")    # type: ignore
 
         # Conditions on cmake variables set from cmake/project_options
         if os.environ["CMAKE_ENABLE_TEST"] == "ON":

doc/CMakeLists.txt

@@ -31,6 +31,10 @@ set(DOXYGEN_CLASS_GRAPH NO)
 set(DOXYGEN_COLLABORATION_GRAPH NO)
 set(DOXYGEN_JAVADOC_AUTOBRIEF YES)
 set(DOXYGEN_USE_MATHJAX YES)
+set(DOXYGEN_MATHJAX_EXTENSIONS "ams" "boldsymbol")
+set(DOXYGEN_FORMULA_MACROFILE "${PROJECT_SOURCE_DIR}/doc/latex_macros.tex")
+set(DOXYGEN_LATEX_BIB_STYLE "unsrtnat")
+set(DOXYGEN_CITE_BIB_FILES "${PROJECT_SOURCE_DIR}/doc/references.bib")
 set(DOXYGEN_MATHJAX_VERSION "MathJax_3")
 set(DOXYGEN_MATHJAX_CODEFILE "${PROJECT_SOURCE_DIR}/util/mathjax-config.js")
 set(DOXYGEN_RECURSIVE YES)

doc/latex_macros.tex

@@ -0,0 +1,2 @@
+\newcommand{\argmin}[1]{\mathop{\mathrm{arg\,min}}\limits_{#1}}
+\newcommand{\argmax}[1]{\mathop{\mathrm{arg\,max}}\limits_{#1}}

doc/main.md

@@ -11,6 +11,10 @@ This program includes a header-only library for the direct integration of the al
 - \subpage installation "Installation"
 - \subpage basic_usages "Basic usages"
 
+<h2> Theories </h2>
+
+- \subpage linear_reg "Linear regression"
+
 <h2> For Developers</h2>
 
 - \subpage program_style_conventions "Programming styles and conventions"

doc/references.bib

@@ -0,0 +1,19 @@
+@misc{ wiki:GausNewton,
+    author = "{Wikipedia contributors}",
+    title = "Gauss–Newton algorithm --- {Wikipedia}{,} The Free Encyclopedia",
+    year = "2025",
+    howpublished = "\url{https://en.wikipedia.org/w/index.php?title=Gauss%E2%80%93Newton_algorithm&oldid=1295135960}",
+  }
+
+@techreport{boggs1989orthogonal,
+  author      = {Boggs, Paul T. and Rogers, Janet E.},
+  title       = {Orthogonal Distance Regression},
+  institution = {National Institute of Standards and Technology},
+  type        = {NISTIR},
+  number      = {89-4197},
+  address     = {Gaithersburg, MD},
+  year        = {1989},
+  month       = nov,
+  url         = {https://www.mechanicalkern.com/static/odr_ams.pdf},
+  note        = {Revised July 1990. Equation (1.3) cited.}
+}

doc/theories/linear_reg.md

@@ -0,0 +1,117 @@
+# Theoretical background {#linear_reg}
+
+[TOC]
+
+The millepede algorithm is fundamentally a linear regression model, which minimizes the _objective function_ with respect to parameters in a linear equation, based on a set of available data points. In this page, the traditional linear regression is explained, followed by a revised version better suited for millepede algorithm and the real data characteristics.
+
+## Simple linear regression
+
+### Objective functions
+
+In the simplest case of linear regression on a calibration equation:
+
+<!-- prettier-ignore-start -->
+\f{equation}{
+y = f(x, a, b) = a \cdot x + b
+\f}
+<!-- prettier-ignore-end -->
+
+where @f$(x, y)@f$ are available data points with parameters @f$a@f$ and @f$b@f$ that need to be optimized. The first-principle approach of parameter optimization is through the so-called _Maximum Likelihood Estimation_, which assumes that @f$y@f$ values follow the normal distribution and the probability to get such @f$y@f$ values being observed is at its maximum when parameters are opitimized. Thus, the probability of the observed data sets @f$(x_i, y_i)@f$ can be express as:
+
+<!-- prettier-ignore-start -->
+\f{equation}{
+Prob(\mathbf{x}, \mathbf{y}) \sim \prod^{n}_{i}\exp{ -\frac{\left(y_i - f(x_i, a, b)\right)^2}{2\sigma_i^2}}
+\label{eq:MLE}
+\f}
+<!-- prettier-ignore-end -->
+
+where @f$\mathbf{x}@f$ and @f$\mathbf{y}@f$ represents observed data values @f$(x_0, x_1, \ldots, x_n)@f$ and @f$(y_0, y_1, \ldots, y_n)@f$ respectively, and @f$\sigma_i@f$ represents the @f$y@f$ error of @f$i@f$-th data point. Maximizing the formula above is equally maximizing its logarithmic transformed formula:
+
+<!-- prettier-ignore-start -->
+\f{align}{
+\argmax{a, b} \prod^{n}_{i}\exp{ -\frac{\left(y_i - f(x_i, a, b)\right)^2}{2\sigma_i^2}} &= \argmax{a, b} \log\left(\prod^{n}_{i}\exp{ -\frac{\left(y_i - f(x_i, a, b)\right)^2}{2\sigma_i^2}}\right) \notag \\
+&= \argmax{a, b} \sum^{n}_{i} -\frac{\left(y_i - f(x_i, a, b)\right)^2}{2\sigma_i^2} \notag \\
+&= \argmin{a, b} \sum^{n}_{i}\frac{\left(y_i - f(x_i, a, b)\right)^2}{2\sigma_i^2}
+\f}
+<!-- prettier-ignore-end -->
+
+Thus, a probability maximization turns into the minimization of a function, which is defined as the **objective function**. By defining another term, called **residual**, denoted as
+
+<!-- prettier-ignore-start -->
+\f{equation}{
+z = \bar{y} - f(\bar{x}, a, b)
+\label{eq:residual}
+\f}
+<!-- prettier-ignore-end -->
+
+where @f$\bar{x}@f$ and @f$\bar{y}@f$ are both observed data, the minimization can then be expressed as:
+
+<!-- prettier-ignore-start -->
+\f{equation}{
+\hat{a}, \hat{b} = \argmin{a, b} \sum^{n}_{i}\frac{z(x_i, y_i, a, b)^2}{2\sigma_i^2}
+\f}
+<!-- prettier-ignore-end -->
+
+In general cases where @f$m@f$ parameters need to be optimized, denoted as @f$\mathbf{p} = (p_0, p_1, \ldots, p_m)@f$, the _objective function_ is expressed as:
+
+<!-- prettier-ignore-start -->
+\f{equation}{
+\mathcal{F}(\mathbf{p}) = \sum^{n}_{i}\frac{\left(y_i - f(x_i, \mathbf{p})\right)^2}{2\sigma_i^2}
+\f}
+<!-- prettier-ignore-end -->
+
+where @f$f(x_i, \mathbf{p})@f$ can be any functions, including non-linear ones in parameters @f$\mathcal{p}@f$.
+
+### Minimization process
+
+The minimization process in this project is using so-called _Gaussian Newton's Method_ \cite wiki:GausNewton, which basically assumes the linear relations of the optimization parameters in the calibration equation:
+
+<!-- prettier-ignore-start -->
+\f{equation}{
+f(x, \mathbf{p}) = f(x, \mathbf{p}) \bigg\rvert_{\mathbf{p} = \mathbf{p}_\text{init}} + \nabla_{\mathbf{p}} f(x, \mathbf{p}) \bigg\rvert_{\mathbf{p} = \mathbf{p}_\text{init}} \cdot(\mathbf{p} - \mathbf{p}_\text{init})
+\f}
+<!-- prettier-ignore-end -->
+
+Thus, the iteration update on the parameters @f$\mathbf{p}@f$, which is derived from the traditional _Newton's Method_, can be further simplified as:
+
+<!-- prettier-ignore-start -->
+\f{flalign}{
+   & & \nabla_{\mathbf{p}}^2 \mathcal{F} \, \delta \mathbf{p} &= \nabla_{\mathbf{p}} \mathcal{F} \notag & \\
+   &\implies & \nabla_{\mathbf{p}}^2 \left(\sum^{n}_{i}\frac{\left(y_i - f(x_i, \mathbf{p})\right)^2}{2\sigma_i^2} \right) \, \delta\mathbf{p} &= \nabla_{\mathbf{p}} \left( \sum^{n}_{i}\frac{\left(y_i - f(x_i, \mathbf{p})\right)^2}{2\sigma_i^2} \right) \notag & \\
+   &\implies & \left(\sum^{n}_i \frac{\nabla_{\mathbf{p}} f(x_i, \mathbf{p}) \nabla_{\mathbf{p}}^{\dagger} f(x_i, \mathbf{p})}{\sigma_i^2} \right) \, \delta\mathbf{p} &= - \sum^{n}_i \frac{z(x_i, y_i, \mathbf{p}) \nabla_{\mathbf{p}} f(x_i, \mathbf{p})}{\sigma_i^2} &\\
+\f}
+<!-- prettier-ignore-end -->
+
+where @f$z(x_i, y_i, \mathbf{p})@f$ is defined from equation @f$\eqref{eq:residual}@f$. The left-hand side matrix, @f$\nabla_{\mathbf{p}}^2 \mathcal{F}@f$ is called _Hessian matrix_, which is simplified to be the external product of the gradient of the calibration function@f$f(x, \mathbf{p})@f$ in the Gaussian Newton's method. There are four prerequisites for a successful minimization:
+
+1. Initial values of all parameters @f$\mathbf{p}@f$ must be available.
+2. Error values @f$\sigma_i@f$ must not be zero.
+3. Hessian matrix on the left-hand side must not have rank-deficit and must be invertible.
+4. Hessian matrix on the left-hand side must be positive definite.
+
+Fortunately, a majority of the calibration equations concerning real-world detectors already fulfilled these prerequisites, except the third one, which can be easily fixed by fixing some parameters (see section). Nevertheless, there are two more issues when applying this minimization to the real experimental data:
+
+* Both @f$x@f$ and @f$y@f$ have error values. Ignoring @f$x@f$ error values could cause inaccuracies on the optimized parameter values.
+* Number of parameters could grow very large when the calibration involves with data from lots of tracks. At some points, trying to invert the Hessian matrix becomes impractical and very expensive.
+
+Thus, in the next sessions, several methods addressing these issues will be explained in detail.
+
+## Regression with both x and y errors
+
+To introduce @f$x@f$ errors both in the objective function and minimization step, some changes are required when calculating the maximum probability in equation @f$\eqref{eq:MLE}@f$. Assuming both @f$x@f$ and @f$y@f$ values follow the Gaussian distribution and they are independent, the probability for the observed data points can be expressed as:
+
+<!-- prettier-ignore-start -->
+\f{equation}{
+    Prob(\mathbf{x}, \mathbf{y}) \sim \prod^{n}_{i} \exp{ -\frac{\left(y_i - f(\mu_i, \mathbf{p})\right)^2}{2(\sigma^y_i)^2}} \exp{ -\frac{\left(x_i - \mu_i\right)^2}{2(\sigma^x_i)^2}}
+\f}
+<!-- prettier-ignore-end -->
+
+where @f$\mu_i@f$ represents the true @f$x@f$ values and @f$f(\mu_i, \mathbf{p})@f$ the true @f$y@f$ values. Following a similar derivation, the objective function can then be expressed as
+
+<!-- prettier-ignore-start -->
+\f{equation}{
+    \mathcal{F}(\mathbf{p}, \boldsymbol{\mu}) = \sum^{n}_{i}\left(\frac{\left(y_i - f(\mu_i, \mathbf{p})\right)^2}{2(\sigma^y_i)^2} + \frac{\left(x_i - \mu_i\right)^2}{2(\sigma^x_i)^2} \right)
+\f}
+<!-- prettier-ignore-end -->
+
+It can be seen here that the @f$x@f$ true values @f$\boldsymbol{\mu}@f$ are treated as additional optimization parameters and its form is exactly the same as the objective function derived in the _Orthogonal Distance Regression_ (ODR)\cite boggs1989orthogonal.

source/centipede/CMakeLists.txt

@@ -1,6 +1,6 @@
 add_library(core SHARED)
 add_library(centipede::centipede ALIAS core)
-target_compile_features(core PRIVATE cxx_std_26)
+target_compile_features(core PUBLIC cxx_std_26)
 target_compile_options(
     core
     PRIVATE
@@ -22,6 +22,14 @@ target_sources(
             FILES centipede.hpp
 )
 
+target_link_libraries(core PUBLIC Eigen3::Eigen GSL::gsl)
+
+if(libassert_FOUND)
+    message(STATUS "Libassert is available")
+    target_link_libraries(core PUBLIC libassert::assert)
+    target_compile_definitions(core PUBLIC HAS_LIBASSERT LIBASSERT_LOWERCASE)
+endif()
+
 add_subdirectory(core)
 add_subdirectory(util)
 add_subdirectory(writer)

source/centipede/core/CMakeLists.txt

@@ -1,6 +1,5 @@
 target_sources(
     core
-    PRIVATE handler.cpp
     PUBLIC
         FILE_SET publicHeaders
             TYPE HEADERS
@@ -10,7 +9,6 @@ target_sources(
                 engines/master_engine.hpp
                 handler.hpp
 )
-target_link_libraries(core PUBLIC Eigen3::Eigen GSL::gsl)
 target_compile_definitions(
     core
     PUBLIC

source/centipede/core/engines/base_engine.hpp

@@ -65,9 +65,10 @@ namespace centipede::core::engine
          *
          * @param self Reference to the caller object.
          * @param entry Entry data.
+         * @return Possible error.
          * @see #Entry
          */
-        void fill_data(this auto&& self, const Entry<DataType>& entry);
+        auto fill_data(this auto&& self, const Entry<DataType>& entry) -> VoidError;
 
         /**
          * @brief Analyze the data from the current entry.
@@ -83,7 +84,7 @@ namespace centipede::core::engine
          * @param alpha Significance level to reject the current entry.
          * @return An error value if error occurs.
          */
-        auto analyze(this auto&& self, double alpha) -> EnumError<>;
+        auto analyze(this auto&& self, double alpha) -> VoidError;
 
         [[nodiscard]] auto get_current_state() const -> const auto& { return state_; }
         [[nodiscard]] auto get_log() const -> const auto& { return log_; }
@@ -108,11 +109,11 @@ namespace centipede::core::engine
     };
 
     template <typename DataType>
-    void Base<DataType>::fill_data(this auto&& self, const Entry<DataType>& entry)
+    auto Base<DataType>::fill_data(this auto&& self, const Entry<DataType>& entry) -> VoidError
     {
         if (not entry.n_locals)
         {
-            return;
+            return {};
         }
         self.state_.n_points = entry.measurements.size();
         assert(self.state_.n_points == entry.sigmas.size());
@@ -123,21 +124,22 @@ namespace centipede::core::engine
         self.fill_measurements(entry.measurements);
         self.fill_sigmas(entry.sigmas);
         self.fill_local_derivs(entry.local_derivs);
-        self.fill_global_derivs(entry.global_derivs);
+        auto res = self.fill_global_derivs(entry.global_derivs);
 
         ++self.log_.n_entries_read;
+        return res;
     }
 
     template <typename DataType>
-    auto Base<DataType>::analyze(this auto&& self, double alpha) -> EnumError<>
+    auto Base<DataType>::analyze(this auto&& self, double alpha) -> VoidError
     {
         return self.fit_local_pars()
             .and_then([&self]() -> EnumError<std::pair<std::size_t, double>>
                       { return self.calculate_local_fit_chi_square(); })
             .and_then(
                 [&self, alpha](const auto& ndf_chi2) -> EnumError<>
                 {
-                    const auto p_value = gsl_cdf_chisq_Q(ndf_chi2.second, ndf_chi2.first);
+                    const auto p_value = gsl_cdf_chisq_Q(ndf_chi2.second, static_cast<double>(ndf_chi2.first));
 
                     self.state_.p_value = p_value;
                     self.state_.chi2 = ndf_chi2.second;
@@ -178,7 +180,7 @@ namespace centipede::core::engine
     /**
      * @brief Empty base engine class. The real implementation is defined in its specialization.
      */
-    template <MatrixEngineType engine_type, typename DataType>
+    template <MatrixEngine engine_type, typename DataType>
     class Engine
     {
     };