diff --git a/cmake/load_dependencies.cmake b/cmake/load_dependencies.cmake index bd67d89a..0fe4bb5d 100644 --- a/cmake/load_dependencies.cmake +++ b/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) diff --git a/conanfile.py b/conanfile.py index 662b4bb5..c7b01f44 100644 --- a/conanfile.py +++ b/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": diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index 78a0c3e6..db8d50e3 100644 --- a/doc/CMakeLists.txt +++ b/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) diff --git a/doc/latex_macros.tex b/doc/latex_macros.tex new file mode 100644 index 00000000..8b518e31 --- /dev/null +++ b/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}} diff --git a/doc/main.md b/doc/main.md index e7d63165..22e3e1d5 100644 --- a/doc/main.md +++ b/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" +

Theories

+ +- \subpage linear_reg "Linear regression" +

For Developers

- \subpage program_style_conventions "Programming styles and conventions" diff --git a/doc/references.bib b/doc/references.bib new file mode 100644 index 00000000..5f03a7bf --- /dev/null +++ b/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.} +} diff --git a/doc/theories/linear_reg.md b/doc/theories/linear_reg.md new file mode 100644 index 00000000..fc0c4645 --- /dev/null +++ b/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: + + +\f{equation}{ +y = f(x, a, b) = a \cdot x + b +\f} + + +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: + + +\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} + + +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: + + +\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} + + +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 + + +\f{equation}{ +z = \bar{y} - f(\bar{x}, a, b) +\label{eq:residual} +\f} + + +where @f$\bar{x}@f$ and @f$\bar{y}@f$ are both observed data, the minimization can then be expressed as: + + +\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} + + +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: + + +\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} + + +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: + + +\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} + + +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: + + +\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} + + +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: + + +\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} + + +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 + + +\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} + + +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. diff --git a/source/centipede/CMakeLists.txt b/source/centipede/CMakeLists.txt index 7d814919..a6d195df 100644 --- a/source/centipede/CMakeLists.txt +++ b/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) diff --git a/source/centipede/core/CMakeLists.txt b/source/centipede/core/CMakeLists.txt index bfe3f77b..db336665 100644 --- a/source/centipede/core/CMakeLists.txt +++ b/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 diff --git a/source/centipede/core/engines/base_engine.hpp b/source/centipede/core/engines/base_engine.hpp index 8c3bdba8..8849a771 100644 --- a/source/centipede/core/engines/base_engine.hpp +++ b/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& entry); + auto fill_data(this auto&& self, const Entry& 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 - void Base::fill_data(this auto&& self, const Entry& entry) + auto Base::fill_data(this auto&& self, const Entry& 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,13 +124,14 @@ 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 - auto Base::analyze(this auto&& self, double alpha) -> EnumError<> + auto Base::analyze(this auto&& self, double alpha) -> VoidError { return self.fit_local_pars() .and_then([&self]() -> EnumError> @@ -137,7 +139,7 @@ namespace centipede::core::engine .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(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 + template class Engine { }; diff --git a/source/centipede/core/engines/eigen_engine.hpp b/source/centipede/core/engines/eigen_engine.hpp index 41165042..4bd49d97 100644 --- a/source/centipede/core/engines/eigen_engine.hpp +++ b/source/centipede/core/engines/eigen_engine.hpp @@ -4,6 +4,7 @@ #include "centipede/core/engines/engine_types.hpp" #include "centipede/core/engines/result.hpp" #include "centipede/data/entry.hpp" +#include "centipede/util/common_definitions.hpp" #include "centipede/util/error_types.hpp" #include "centipede/util/return_types.hpp" #include @@ -20,14 +21,39 @@ #include #include +#ifdef HAS_LIBASSERT +#include "libassert/assert.hpp" +#else +#include +#endif + namespace centipede::core::engine { + /** * @brief Engine template specialization for Eigen library implementation. */ template - class Engine : public Base + class Engine : public Base { + private: + /** + * @brief RAII class to perform the malloc check for eigen matrix operations. + * + */ + class EigenMemGuard + { + public: + /** + * @brief Default constructor to enable the check. + */ + EigenMemGuard() { Eigen::internal::set_is_malloc_allowed(false); } + /** + * @brief Default destructor to disable the check. + */ + ~EigenMemGuard() { Eigen::internal::set_is_malloc_allowed(true); } + }; + public: /** * @brief Matrix and vector used to solve global parameter updates @@ -55,7 +81,13 @@ namespace centipede::core::engine */ static void solve(const Globals& globals, Result& result) { - assert(globals.factor_matrix.isApprox(globals.factor_matrix.transpose())); +#ifdef HAS_LIBASSERT + debug_assert(globals.factor_matrix.isApprox(globals.factor_matrix.transpose(), + static_cast(common::EIGEN_APPROX_PRECISION))); +#else + assert(globals.factor_matrix.isApprox(globals.factor_matrix.transpose(), + static_cast(common::EIGEN_APPROX_PRECISION))); +#endif if (globals.factor_matrix.isZero()) { @@ -83,14 +115,25 @@ namespace centipede::core::engine std::back_inserter(result.parameters)); result.error_status = ErrorCode::success; } - else - { - check_rank_deficit(globals, result); - } + // else + // { + check_rank_deficit(globals, result); + // } } void add_to_globals(Globals& globals) { + globals.factor_matrix.resize(globals_.factor_matrix.rows(), globals_.factor_matrix.cols()); + globals.rhs_vec.resize(globals_.rhs_vec.rows(), globals_.rhs_vec.cols()); +#ifdef HAS_LIBASSERT + debug_assert(globals.factor_matrix.cols() == globals_.factor_matrix.cols(), + "Check if column size matches for global factor matrix"); + debug_assert(globals.factor_matrix.rows() == globals_.factor_matrix.rows(), + "Check if row size matches for global factor matrix"); +#else + assert(globals.factor_matrix.cols() == globals_.factor_matrix.cols()); + assert(globals.factor_matrix.rows() == globals_.factor_matrix.rows()); +#endif globals.factor_matrix += globals_.factor_matrix.eval(); globals.rhs_vec += globals_.rhs_vec.eval(); } @@ -136,9 +179,9 @@ namespace centipede::core::engine { const auto& current_state = Base::get_current_state(); - const auto entrypoint_size = current_state.n_points; - const auto n_globals = current_state.n_globals; - const auto n_locals = current_state.n_locals; + const auto entrypoint_size = static_cast(current_state.n_points); + const auto n_globals = static_cast(current_state.n_globals); + const auto n_locals = static_cast(current_state.n_locals); // NOTE: resize may cause memory allocation. local_t_.resize(n_locals, entrypoint_size); @@ -163,8 +206,10 @@ namespace centipede::core::engine void fill_sigmas(const std::vector& data) { - std::ranges::copy(std::views::transform(data, [](DataType val) -> DataType { return 1. / (val * val); }), - sigmas_.begin()); + std::ranges::copy( + std::views::transform(data, + [](DataType val) -> DataType { return static_cast(1.) / (val * val); }), + sigmas_.begin()); } void fill_measurements(const std::vector& data) { std::ranges::copy(data, measurements_.begin()); } @@ -173,13 +218,18 @@ namespace centipede::core::engine { for (const auto& [point_idx, deriv] : data) { +#if HAS_LIBASSERT + debug_assert(point_idx < local_t_.cols()); + debug_assert(deriv.first < local_t_.rows()); +#else assert(point_idx < local_t_.cols()); assert(deriv.first < local_t_.rows()); +#endif local_t_(deriv.first, point_idx) = deriv.second; } } - void fill_global_derivs(const std::vector::Deriv>& data) + auto fill_global_derivs(const std::vector::Deriv>& data) -> VoidError { triplets_.clear(); @@ -190,12 +240,13 @@ namespace centipede::core::engine triplets_.emplace_back(deriv.first, point_idx, deriv.second); } global_t_.setFromSortedTriplets(triplets_.begin(), triplets_.end()); + return {}; } - auto fit_local_pars() -> EnumError<> + auto fit_local_pars() -> VoidError { // NOTE: Multiplications will trigger temporary object (memory allocation later during the assignment.) - Eigen::internal::set_is_malloc_allowed(false); + auto _ = EigenMemGuard{}; buffers_.local_weighted_t.noalias() = local_t_ * sigmas_.asDiagonal(); buffers_.local_weighted_square.noalias() = buffers_.local_weighted_t.lazyProduct(local_t_.transpose()); @@ -211,17 +262,15 @@ namespace centipede::core::engine buffers_.local_solutions.noalias() = buffers_.local_weighted_square_inv.lazyProduct(buffers_.local_weighted_t).lazyProduct(measurements_); - Eigen::internal::set_is_malloc_allowed(true); - return {}; } auto calculate_local_fit_chi_square() -> EnumError> { - Eigen::internal::set_is_malloc_allowed(false); + auto _ = EigenMemGuard{}; const auto entrypoint_size = Base::get_current_state().n_points; const auto local_size = buffers_.local_solutions.rows(); - const auto ndf = entrypoint_size - local_size; + const auto ndf = entrypoint_size - static_cast(local_size); if (ndf < 1) { @@ -230,11 +279,10 @@ namespace centipede::core::engine buffers_.residual_values.noalias() = measurements_ - (local_t_.transpose() * buffers_.local_solutions); const auto chi_square = buffers_.residual_values.dot(sigmas_.asDiagonal() * buffers_.residual_values); - Eigen::internal::set_is_malloc_allowed(true); return std::pair{ ndf, chi_square }; } - auto update_global_factor_matrix() -> EnumError<> + auto update_global_factor_matrix() -> VoidError { // TODO: Perform the production using index accessing. // NOTE: Seems that there is no way to prevent memory allocation with sparse matrices. @@ -251,13 +299,20 @@ namespace centipede::core::engine buffers_.local_weighted_square_inv_sparse * buffers_.global_local_weighted_t; buffers_.global_square_update += buffers_.global_weighted_square; - assert(buffers_.global_square_update.isApprox(buffers_.global_square_update.transpose())); +#ifdef HAS_LIBASSERT + debug_assert(buffers_.global_square_update.isApprox(buffers_.global_square_update.transpose(), + static_cast(common::EIGEN_APPROX_PRECISION))); +#else + assert(buffers_.global_square_update.isApprox(buffers_.global_square_update.transpose(), + static_cast(common::EIGEN_APPROX_PRECISION))); +#endif + globals_.factor_matrix += buffers_.global_square_update; // Eigen::internal::set_is_malloc_allowed(true); return {}; } - auto update_global_rhs_vector() -> EnumError<> + auto update_global_rhs_vector() -> VoidError { // TODO: Perform the production using index accessing. // NOTE: Seems that there is no way to prevent memory allocation with sparse matrices. @@ -274,9 +329,10 @@ namespace centipede::core::engine static void resize_globals(Globals& globals, std::size_t n_globals) { - globals.rhs_vec.resize(n_globals); + const auto num_of_globals = static_cast(n_globals); + globals.rhs_vec.resize(num_of_globals); globals.rhs_vec.setZero(); - globals.factor_matrix.resize(n_globals, n_globals); + globals.factor_matrix.resize(num_of_globals, num_of_globals); globals.factor_matrix.setZero(); } @@ -296,7 +352,7 @@ namespace centipede::core::engine } } const auto diagonal_values = (unitary_matrix * prob_mat * unitary_matrix.transpose()).diagonal().eval(); - for (const auto [idx, diagonal_val] : std::views::zip(std::views::iota(0), diagonal_values)) + for (const auto [idx, diagonal_val] : std::views::zip(std::views::iota(std::size_t{ 0 }), diagonal_values)) { if (std::abs(diagonal_val) > Eigen::NumTraits::dummy_precision()) { diff --git a/source/centipede/core/engines/engine_concept.hpp b/source/centipede/core/engines/engine_concept.hpp index a13c5658..299b7c42 100644 --- a/source/centipede/core/engines/engine_concept.hpp +++ b/source/centipede/core/engines/engine_concept.hpp @@ -14,7 +14,7 @@ namespace centipede::core::engine /** * @brief Concept used for core::engine::Master option. */ - template + template concept EngineLike = requires(Engine engine, Result& result, typename Engine::Globals& globals) { @@ -26,8 +26,8 @@ namespace centipede::core::engine { Engine::solve(globals, result) } -> std::same_as; { engine.add_to_globals(globals) } -> std::same_as; { engine.add_to_result(result) } -> std::same_as; - { engine.analyze(double{}) } -> std::same_as>; - { engine.fill_data(Entry{}) } -> std::same_as; + { engine.analyze(double{}) } -> std::same_as; + { engine.fill_data(Entry{}) } -> std::same_as; }; } // namespace centipede::core::engine diff --git a/source/centipede/core/engines/engine_types.hpp b/source/centipede/core/engines/engine_types.hpp index daddd5bf..3becb757 100644 --- a/source/centipede/core/engines/engine_types.hpp +++ b/source/centipede/core/engines/engine_types.hpp @@ -8,7 +8,7 @@ namespace centipede::core::engine * @brief Engine types. * */ - enum class MatrixEngineType : uint8_t + enum class MatrixEngine : uint8_t { eigen, xtensor, @@ -20,7 +20,7 @@ namespace centipede::core::engine */ struct MasterOpt { - MatrixEngineType engine_type = MatrixEngineType::eigen; + MatrixEngine engine_type = MatrixEngine::eigen; bool has_multi_slaves = false; //!< Use taskflow or not. }; @@ -28,5 +28,5 @@ namespace centipede::core::engine namespace centipede { - using EngineType = core::engine::MatrixEngineType; + using MatrixEngine = core::engine::MatrixEngine; } diff --git a/source/centipede/core/engines/master_engine.hpp b/source/centipede/core/engines/master_engine.hpp index bfd0ac01..658aa623 100644 --- a/source/centipede/core/engines/master_engine.hpp +++ b/source/centipede/core/engines/master_engine.hpp @@ -6,7 +6,7 @@ #include "centipede/core/engines/engine_types.hpp" #include "centipede/core/engines/result.hpp" #include "centipede/data/entry.hpp" -#include "centipede/data/entry_base.hpp" +#include "centipede/data/entrypoint.hpp" #include "centipede/util/error_types.hpp" #include "centipede/util/return_types.hpp" #include @@ -44,11 +44,11 @@ namespace centipede::core::engine */ struct State { - std::size_t point_index = 0; //!< Index used for each new entrypoint added. - Entry entry; //!< Data storing the current entry. + std::size_t next_point_index = 0; //!< Index used for each new entrypoint added. + Entry entry; //!< Data storing the current entry. }; - using Result = Result; + using ResultType = Result; using EngineImp = Engine; using DataTypeUsed = DataType; @@ -69,64 +69,55 @@ namespace centipede::core::engine * entry. */ template - [[nodiscard]] auto add_entrypoint(const EntryPoint& entry_point) -> EnumError<> + [[nodiscard]] auto add_entrypoint(const EntryPoint& entry_point) -> VoidError { - const auto n_locals = [&entry_point]() -> std::size_t - { - if constexpr (NLocals == internal::DYNAMIC_SIZE or NGlobals == internal::DYNAMIC_SIZE) - { - return entry_point.get_locals().size(); - } - else - { - return NLocals; - } - }(); - if (current_state_.entry.n_locals.has_value() and current_state_.entry.n_locals.value() != n_locals) - { - return std::unexpected{ ErrorCode::handler_incomp_n_locals }; - } - current_state_.entry.n_locals = n_locals; - - current_state_.entry.measurements.push_back(entry_point.get_measurement()); - current_state_.entry.sigmas.push_back(entry_point.get_sigma()); - std::ranges::copy(std::views::zip_transform( - [this](auto local_idx, auto deriv) -> Entry::Deriv - { return std::pair{ current_state_.point_index, std::pair{ local_idx, deriv } }; }, - std::views::iota(0), - entry_point.get_locals()), - std::back_inserter(current_state_.entry.local_derivs)); - std::ranges::copy(entry_point.get_globals() | - std::views::transform([this](const auto& deriv) -> Entry::Deriv - { return std::pair{ current_state_.point_index, deriv }; }), - std::back_inserter(current_state_.entry.global_derivs)); - std::ranges::sort(current_state_.entry.global_derivs, - [](const Entry::Deriv& left, const Entry::Deriv& right) -> bool - { - return left.first < right.first || - ((left.first == right.first) && (left.second.first < right.second.first)); - }); - - ++current_state_.point_index; - return {}; + return check_entrypoint_valid(entry_point) + .transform( + [this, &entry_point]() + { + current_state_.entry.n_locals = entry_point.get_n_locals(); + + current_state_.entry.measurements.push_back(entry_point.get_measurement()); + current_state_.entry.sigmas.push_back(entry_point.get_sigma()); + std::ranges::copy( + std::views::zip_transform( + [this](auto local_idx, auto deriv) -> Entry::Deriv + { return std::pair{ current_state_.next_point_index, std::pair{ local_idx, deriv } }; }, + std::views::iota(0), + entry_point.get_locals()), + std::back_inserter(current_state_.entry.local_derivs)); + std::ranges::copy( + entry_point.get_globals() | + std::views::transform([this](const auto& deriv) -> Entry::Deriv + { return std::pair{ current_state_.next_point_index, deriv }; }), + std::back_inserter(current_state_.entry.global_derivs)); + std::ranges::sort( + current_state_.entry.global_derivs, + [](const Entry::Deriv& left, const Entry::Deriv& right) -> bool + { + return left.first < right.first || + ((left.first == right.first) && (left.second.first < right.second.first)); + }); + ++current_state_.next_point_index; + }); } /** * @brief Fitting the current entry data. * - * This operation can be async. + * This operation can be async. The state is also reset to the default one. * @see ref */ - auto analyze() -> EnumError<> + auto analyze() -> VoidError { - engine_imp_.fill_data(current_state_.entry); - auto res = engine_imp_.analyze(config_.alpha); - reset_state(); - if (not res) - { - return std::unexpected{ res.error() }; - } - return {}; + return engine_imp_.fill_data(current_state_.entry) + .and_then( + [this]() + { + auto res = engine_imp_.analyze(config_.alpha); + reset_state(); + return res; + }); } /** @@ -135,17 +126,16 @@ namespace centipede::core::engine * @return name description * @see ref */ - auto solve() -> EnumError<> + auto solve() -> VoidError { engine_imp_.add_to_globals(globals_); engine_imp_.add_to_result(result_); - engine_imp_.solve(globals_, result_); + EngineImp::solve(globals_, result_); - return (result_.error_status == ErrorCode::success) ? EnumError<>{} - : std::unexpected{ result_.error_status }; + return (result_.error_status == ErrorCode::success) ? VoidError{} : std::unexpected{ result_.error_status }; } - [[nodiscard]] auto get_current_state() const -> const auto& { return current_state_; } + [[nodiscard]] auto get_current_state() const -> const State& { return current_state_; } [[nodiscard]] auto get_engine() const -> const auto& { return engine_imp_; } @@ -153,20 +143,39 @@ namespace centipede::core::engine private: Config config_; - Result result_; + ResultType result_; State current_state_; EngineImp engine_imp_{}; + + // TODO: represent globals as mdspan, instead of relying on eigen. EngineImp::Globals globals_{}; void reset_state() { - current_state_.point_index = 0; + current_state_.next_point_index = 0; current_state_.entry.global_derivs.clear(); current_state_.entry.local_derivs.clear(); current_state_.entry.measurements.clear(); current_state_.entry.sigmas.clear(); current_state_.entry.n_locals.reset(); } + + template + auto check_entrypoint_valid(const EntryPoint& entry_point) -> VoidError + { + const auto n_locals = entry_point.get_n_locals(); + if (current_state_.entry.n_locals.has_value() and current_state_.entry.n_locals.value() != n_locals) + { + return std::unexpected{ ErrorCode::handler_incomp_n_locals }; + } + const auto n_globals = config_.n_globals; + if (std::ranges::any_of(entry_point.get_globals(), + [n_globals](const auto& idx_value) { return idx_value.first >= n_globals; })) + { + return std::unexpected{ ErrorCode::analysis_global_idx_too_large }; + } + return {}; + } }; } // namespace centipede::core::engine diff --git a/source/centipede/core/engines/result.hpp b/source/centipede/core/engines/result.hpp index a9b45fb6..56d6fc76 100644 --- a/source/centipede/core/engines/result.hpp +++ b/source/centipede/core/engines/result.hpp @@ -48,8 +48,9 @@ struct std::formatter> static constexpr auto parse(std::format_parse_context& ctx) { return ctx.begin(); } static constexpr auto format(const Result& result, std::format_context& ctx) { - const auto percentage = - (result.n_entries == 0) ? 0. : result.n_entries_rejected / static_cast(result.n_entries) * 100.; + const auto percentage = (result.n_entries == 0) ? 0. + : static_cast(result.n_entries_rejected) / + static_cast(result.n_entries) * 100.; if (result.error_status == centipede::ErrorCode::success) { return std::format_to(ctx.out(), diff --git a/source/centipede/core/handler.cpp b/source/centipede/core/handler.cpp deleted file mode 100644 index e69de29b..00000000 diff --git a/source/centipede/core/handler.hpp b/source/centipede/core/handler.hpp index 5adc4bf4..09fa2b7c 100644 --- a/source/centipede/core/handler.hpp +++ b/source/centipede/core/handler.hpp @@ -2,11 +2,13 @@ #include "centipede/core/engines/engine_types.hpp" #include "centipede/core/engines/master_engine.hpp" -#include "centipede/data/entry.hpp" +#include "centipede/data/entrypoint.hpp" +#include "centipede/util/error_types.hpp" #include "centipede/util/return_types.hpp" #include +#include -namespace centipede::core +namespace centipede { /** @@ -14,7 +16,7 @@ namespace centipede::core * * This class should handle all inputs and configurations from users */ - template + template class Handler { public: @@ -24,41 +26,62 @@ namespace centipede::core */ struct Config { - std::size_t n_globals = 0; //!< Number of global parameters. + // std::size_t n_globals = 0; //!< Number of global parameters. }; - using EngineType = engine::Master; + using EngineType = core::engine::Master; /** * @brief Constructor * * This is the default constructor */ - explicit Handler(Config config = {}) + explicit Handler(std::size_t n_globals, Config config = {}) : config_{ config } - , engine_{ typename EngineType::Config{ .n_globals = config.n_globals } } + , engine_{ typename EngineType::Config{ .n_globals = n_globals } } { } - /** - * @brief Initialization of the instance - * - */ - [[nodiscard]] auto init() -> EnumError<> { return {}; } + // /** + // * @brief Initialization of the instance + // * + // */ + // [[nodiscard]] auto init() -> EnumError<> { + // engine_.init(); + // return {}; } template - [[nodiscard]] auto add_entrypoint(const EntryPoint& entry_point) -> EnumError<> + [[nodiscard]] auto add_entrypoint(const EntryPoint& entry_point) -> VoidError { return engine_.add_entrypoint(entry_point); } - auto analyze_current_entry() -> EnumError { return {}; }; + auto analyze_current_entry() -> EnumError + { + auto n_points = engine_.get_current_state().entry.measurements.size(); + + if (n_points == 0) + { + return std::unexpected{ ErrorCode::analysis_empty_entry }; + } + + return engine_.analyze().transform([n_points]() { return n_points; }); + }; [[nodiscard]] auto get_current_state() const -> const auto& { return engine_.get_current_state(); } + [[nodiscard]] auto get_current_entry_state() const -> const auto& + { + return engine_.get_engine().get_current_state(); + } + + [[nodiscard]] auto solve() -> VoidError { return engine_.solve(); } + + [[nodiscard]] auto get_result() const -> const auto& { return engine_.get_result(); } + private: Config config_; EngineType engine_; // std::vector<> }; -} // namespace centipede::core +} // namespace centipede diff --git a/source/centipede/data/entry.hpp b/source/centipede/data/entry.hpp index e2fa9309..f59f2535 100644 --- a/source/centipede/data/entry.hpp +++ b/source/centipede/data/entry.hpp @@ -1,170 +1,14 @@ #pragma once -#include "centipede/data/entry_base.hpp" -#include -#include +#include #include #include -#include #include #include #include namespace centipede { - /** - * @brief Structure of a entrypoint with dynamic sizes. - * - * Values of local and global derivative are store in `std::vector`. Values are added to the vector via the setters - * or #add_local and #add_global methods. It's highly recommended to reserve the memory first for the underlying - * data to avoid the memory reallocation. - */ - template - class EntryPoint : public internal::EntryPointBase - { - public: - using LocalDerivs = std::vector; //!< Data type of local deriv array. - using GlobalDerivs = std::vector>; //!< Data type of global deriv array. - - /** - * @brief Default constructor. - * - * All values are initialized to the default values. - * - */ - EntryPoint() = default; - - /** - * @brief Add a value to the local derivatives. - * @param value Local derivative value. - * @return Non-const reference to this object. - */ - constexpr auto add_local(std::floating_point auto value) -> auto& - { - locals_.push_back(static_cast(value)); - return *this; - } - - /** - * @brief Add an index-value pair to the global derivatives. - * @param index Global parameter ID (0-based indexing). - * @param value Global derivative value. - * @return Non-const reference to this object. - */ - constexpr auto add_global(std::integral auto index, std::floating_point auto value) -> auto& - { - globals_.emplace_back(static_cast(index), static_cast(value)); - return *this; - } - - /** - * @brief Set the capacity of the vector storing the local derivatives - * @param size New capacity value. - * @return Non-const reference to this object. - */ - constexpr auto reserve_locals(std::size_t size) -> auto& - { - locals_.reserve(size); - return *this; - } - - /** - * @brief Set the capacity of the vector storing the global derivatives - * @param size New capacity value. - * @return Non-const reference to this object. - */ - constexpr auto reserve_globals(std::size_t size) -> auto& - { - globals_.reserve(size); - return *this; - } - - /** - * @brief Default spaceship comparison. - */ - constexpr auto operator<=>(const EntryPoint& other) const = default; - - private: - friend EntryPointBase; - LocalDerivs locals_; //!< Local derivatives. - GlobalDerivs globals_; //!< Global label and derivatives pair. The label is using **0-based indexing**. - - template - constexpr void set_locals_imp(DataTypes... locals) - { - locals_.clear(); - (locals_.push_back(locals), ...); - } - - template - constexpr void set_globals_imp(DataTypes... globals) - { - globals_.clear(); - (globals_.emplace_back(static_cast(globals.first), static_cast(globals.second)), ...); - } - - constexpr void reset_derivs() - { - locals_.clear(); - globals_.clear(); - } - }; - - /** - * @brief Structure of a entrypoint with static sizes. - * - * Entrypoint contains the derivatives of local and global parameters, together with measurement and sigma values. - * The values are stored locally. - * @tparam NLocals Number of local parameters. - * @tparam NGlobals Number of global parameters. - */ - template - requires(NLocals < internal::DYNAMIC_SIZE and NGlobals < internal::DYNAMIC_SIZE) - class EntryPoint : public internal::EntryPointBase - { - public: - using LocalDerivs = std::array; //!< Data type of local deriv array. - using GlobalDerivs = std::array, NGlobals>; //!< Data type of global deriv array. - - /** - * @brief Default constructor. - * - * All values are initialized to the default values. - * - * @see Config - */ - EntryPoint() = default; - - /** - * @brief Default spaceship comparison. - */ - constexpr auto operator<=>(const EntryPoint& other) const = default; - - private: - friend EntryPointBase; - LocalDerivs locals_{}; //!< Local derivatives. - GlobalDerivs globals_{}; //!< Global label and derivatives pair. The label is using **0-based indexing**. - - template - constexpr void set_locals_imp(DataTypes... locals) - { - locals_ = std::array{ static_cast(locals)... }; - } - - template - constexpr void set_globals_imp(DataTypes... globals) - { - globals_ = - std::array{ std::pair{ static_cast(globals.first), static_cast(globals.second) }... }; - } - - constexpr void reset_derivs() - { - locals_ = std::array{}; - globals_ = std::array, NGlobals>{}; - } - }; - /** * @class Entry * @brief A structure to store the all entrypoints in a entry. @@ -188,23 +32,3 @@ namespace centipede }; }; // namespace centipede - -/** - * @brief Formatter for the class EntryPoint - */ -template -// NOLINTNEXTLINE (bugprone-std-namespace-modification) -struct std::formatter> -{ - static constexpr auto parse(std::format_parse_context& ctx) { return ctx.begin(); } - - static auto format(const centipede::EntryPoint& entry, std::format_context& ctx) - { - return std::format_to(ctx.out(), - "local derivatives: {}, global derivatives: {}, measurement: {}, sigma: {}", - entry.get_locals(), - entry.get_globals(), - entry.get_measurement(), - entry.get_sigma()); - } -}; diff --git a/source/centipede/data/entrypoint.hpp b/source/centipede/data/entrypoint.hpp new file mode 100644 index 00000000..c2ff6184 --- /dev/null +++ b/source/centipede/data/entrypoint.hpp @@ -0,0 +1,296 @@ +#pragma once + +#include "centipede/data/entrypoint_base.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace centipede +{ + /** + * @brief Structure of a entrypoint with dynamic sizes. + * + * Values of local and global derivative are store in `std::vector`. Values are added to the vector via the setters + * or #add_local and #add_global methods. It's highly recommended to reserve the memory first for the underlying + * data to avoid the memory reallocation. + */ + template + class EntryPoint : public internal::EntryPointBase + { + public: + using LocalDerivs = std::vector; //!< Data type of local deriv array. + using GlobalDerivs = std::vector>; //!< Data type of global deriv array. + + /** + * @brief Default constructor. + * + * All values are initialized to the default values. + * + */ + EntryPoint() = default; + + /** + * @brief Getter for local parameter sizes. + * @return Sizes of local parameters. + */ + [[nodiscard]] inline auto get_n_locals() const -> std::size_t { return locals_.size(); } + + /** + * @brief Getter for global parameter sizes. + * @return Sizes of global parameters. + */ + [[nodiscard]] inline auto get_n_globals() const -> std::size_t { return globals_.size(); } + + /** + * @brief Add a value to the local derivatives. + * @param value Local derivative value. + * @return Non-const reference to this object. + */ + constexpr auto add_local(std::floating_point auto value) -> auto& + { + locals_.push_back(static_cast(value)); + return *this; + } + + /** + * @brief Add an index-value pair to the global derivatives. + * @param index Global parameter ID (0-based indexing). + * @param value Global derivative value. + * @return Non-const reference to this object. + */ + constexpr auto add_global(std::integral auto index, std::floating_point auto value) -> auto& + { + globals_.emplace_back(static_cast(index), static_cast(value)); + return *this; + } + + using internal::EntryPointBase::set_globals; + using internal::EntryPointBase::set_locals; + + /** + * @brief Set global derivative values from a view. + * + * The view must reference a std::pair, which must satisfy EntryPointGlobalIdxPair concept. The global derivate + * values are emptied before being set by the new values. + * + * @tparam View Types of the view. + * @param view View object. + * @return Universal reference to the caller. + */ + template + requires EntryPointGlobalIdxPair> + constexpr auto set_globals(View view) -> auto&& + { + globals_.clear(); + std::ranges::copy(view, std::back_inserter(globals_)); + return *this; + } + + /** + * @brief Set local derivative values from a view. + * + * The view must reference a value of DataType. The local derivate values are emptied before being set by the + * new values. + * + * @tparam View Types of the view. + * @param view View object. + * @return Universal reference to the caller. + */ + template + requires std::same_as, float> + constexpr auto set_locals(View view) -> auto&& + { + locals_.clear(); + std::ranges::copy(view, std::back_inserter(locals_)); + return *this; + } + + /** + * @brief Set the capacity of the vector storing the local derivatives + * @param size New capacity value. + * @return Non-const reference to this object. + */ + constexpr auto reserve_locals(std::size_t size) -> auto& + { + locals_.reserve(size); + return *this; + } + + /** + * @brief Set the capacity of the vector storing the global derivatives + * @param size New capacity value. + * @return Non-const reference to this object. + */ + constexpr auto reserve_globals(std::size_t size) -> auto& + { + globals_.reserve(size); + return *this; + } + + /** + * @brief Default spaceship comparison. + */ + constexpr auto operator<=>(const EntryPoint& other) const = default; + + private: + friend EntryPointBase; + LocalDerivs locals_; //!< Local derivatives. + GlobalDerivs globals_; //!< Global label and derivatives pair. The label is using **0-based indexing**. + + template + constexpr void set_locals_imp(DataTypes... locals) + { + locals_.clear(); + (locals_.push_back(locals), ...); + } + + template + constexpr void set_globals_imp(DataTypes... globals) + { + globals_.clear(); + (globals_.emplace_back(static_cast(globals.first), static_cast(globals.second)), ...); + } + + constexpr void reset_derivs() + { + locals_.clear(); + globals_.clear(); + } + }; + + /** + * @brief Structure of a entrypoint with static sizes. + * + * Entrypoint contains the derivatives of local and global parameters, together with measurement and sigma values. + * The values are stored locally. + * @anchor NLocals_DY_Entrypoint + * @tparam NLocals Number of local parameters. + * @anchor NGlobals_DY_Entrypoint + * @tparam NGlobals Number of global parameters. + */ + template + requires(NLocals < internal::DYNAMIC_SIZE and NGlobals < internal::DYNAMIC_SIZE) + class EntryPoint : public internal::EntryPointBase + { + public: + using LocalDerivs = std::array; //!< Data type of local deriv array. + using GlobalDerivs = std::array, NGlobals>; //!< Data type of global deriv array. + + /** + * @brief Default constructor. + * + * All values are initialized to the default values. + * + * @see Config + */ + EntryPoint() = default; + + using internal::EntryPointBase::set_globals; + using internal::EntryPointBase::set_locals; + + /** + * @brief Getter for local parameter sizes. + * @return Sizes of local parameters. + */ + [[nodiscard]] inline auto get_n_locals() const -> std::size_t { return NLocals; } + + /** + * @brief Getter for global parameter sizes. + * @return Sizes of global parameters. + */ + [[nodiscard]] inline auto get_n_globals() const -> std::size_t { return NGlobals; } + + /** + * @brief Default spaceship comparison. + */ + constexpr auto operator<=>(const EntryPoint& other) const = default; + + /** + * @brief Set global derivative values from a view. + * + * The view must reference a std::pair, which must satisfy EntryPointGlobalIdxPair concept. The size of the + * input view is assumed to be larger or equal to @ref NGlobals_DY_Entrypoint "NGlobals". + * + * @tparam View Types of the view. + * @param view View object. + * @return Universal reference to the caller. + */ + template + requires EntryPointGlobalIdxPair> + constexpr auto set_globals(View view) -> auto&& + { + std::ranges::copy(view | std::views::take(NGlobals), globals_.begin()); + return *this; + } + + /** + * @brief Set local derivative values from a view. + * + * The view must reference a value of DataType. The size of the input view is assumed to be larger or equal to + * @ref NLocals_DY_Entrypoint "NLocals". + * + * @tparam View Types of the view. + * @param view View object. + * @return Universal reference to the caller. + */ + template + requires std::same_as, float> + constexpr auto set_locals(View view) -> auto&& + { + std::ranges::copy(view | std::views::take(NLocals), locals_.begin()); + return *this; + } + + private: + friend EntryPointBase; + LocalDerivs locals_{}; //!< Local derivatives. + GlobalDerivs globals_{}; //!< Global label and derivatives pair. The label is using **0-based indexing**. + + template + constexpr void set_locals_imp(DataTypes... locals) + { + locals_ = std::array{ static_cast(locals)... }; + } + + template + constexpr void set_globals_imp(DataTypes... globals) + { + globals_ = + std::array{ std::pair{ static_cast(globals.first), static_cast(globals.second) }... }; + } + + constexpr void reset_derivs() + { + locals_ = std::array{}; + globals_ = std::array, NGlobals>{}; + } + }; +}; // namespace centipede + +/** + * @brief Formatter for the class EntryPoint + */ +template +// NOLINTNEXTLINE (bugprone-std-namespace-modification) +struct std::formatter> +{ + static constexpr auto parse(std::format_parse_context& ctx) { return ctx.begin(); } + + static auto format(const centipede::EntryPoint& entry, std::format_context& ctx) + { + return std::format_to(ctx.out(), + "local derivatives: {}, global derivatives: {}, measurement: {}, sigma: {}", + entry.get_locals(), + entry.get_globals(), + entry.get_measurement(), + entry.get_sigma()); + } +}; diff --git a/source/centipede/data/entry_base.hpp b/source/centipede/data/entrypoint_base.hpp similarity index 99% rename from source/centipede/data/entry_base.hpp rename to source/centipede/data/entrypoint_base.hpp index d41e2b2b..d8522d58 100644 --- a/source/centipede/data/entry_base.hpp +++ b/source/centipede/data/entrypoint_base.hpp @@ -64,7 +64,6 @@ namespace centipede::internal public: using data_type = float; //!< Data type of deriv values. - // /** * @brief Reset all values to the default values. * diff --git a/source/centipede/reader/binary.hpp b/source/centipede/reader/binary.hpp index 2873318c..c7ff8701 100644 --- a/source/centipede/reader/binary.hpp +++ b/source/centipede/reader/binary.hpp @@ -1,6 +1,6 @@ #pragma once -#include "centipede/data/entry.hpp" +#include "centipede/data/entrypoint.hpp" #include "centipede/util/common_definitions.hpp" #include "centipede/util/error_types.hpp" #include "centipede/util/return_types.hpp" diff --git a/source/centipede/util/common_definitions.hpp b/source/centipede/util/common_definitions.hpp index 1854c256..6a4079ee 100644 --- a/source/centipede/util/common_definitions.hpp +++ b/source/centipede/util/common_definitions.hpp @@ -6,4 +6,7 @@ namespace centipede::common { constexpr auto DEFAULT_BUFFER_SIZE = std::size_t{ 10000 }; //!< Default maximum buffer size for binary readers/writers. + + constexpr auto EIGEN_APPROX_PRECISION = + 0.001; //!< Precision to compare whether two floating point values are equal. See Eigen::DenseBase::isApprox(). } // namespace centipede::common diff --git a/source/centipede/util/error_types.hpp b/source/centipede/util/error_types.hpp index 81010444..6a20b461 100644 --- a/source/centipede/util/error_types.hpp +++ b/source/centipede/util/error_types.hpp @@ -26,13 +26,15 @@ namespace centipede analysis_rank_deficit, //!< Rank deficit occurred during the analysis. analysis_zero_global_factor_mat, analysis_global_negative_definite, + analysis_global_idx_too_large, analysis_factor_matrix_zero, //!< Global factor matrix is zero matrix. - analysis_rhs_vector_zero, //!< Global right-hand-side vector is zero vector. - reader_file_fail_to_open, //!< Input file failed to be open. - reader_file_fail_to_read, //!< Input file failed to read - reader_uninitialized, //!< Reader is not initialized. - reader_buffer_overflow, //!< Buffer size is too small for a new entry occurs. See @ref reader::Binary. - reader_invalid_filename, //!< Filename is invalid or empty + analysis_empty_entry, + analysis_rhs_vector_zero, //!< Global right-hand-side vector is zero vector. + reader_file_fail_to_open, //!< Input file failed to be open. + reader_file_fail_to_read, //!< Input file failed to read + reader_uninitialized, //!< Reader is not initialized. + reader_buffer_overflow, //!< Buffer size is too small for a new entry occurs. See @ref reader::Binary. + reader_invalid_filename, //!< Filename is invalid or empty }; } // namespace centipede @@ -83,10 +85,15 @@ struct std::formatter return std::format_to(ctx.out(), "Rank deficit occurred during the analysis."); case analysis_zero_global_factor_mat: return std::format_to(ctx.out(), "Global factor matrix is zero."); + case analysis_global_idx_too_large: + return std::format_to( + ctx.out(), "Global index (0-based) is larger than or equal to the size of global parameters."); case analysis_factor_matrix_zero: return std::format_to(ctx.out(), "Global factor matrix is zero matrix."); case analysis_rhs_vector_zero: return std::format_to(ctx.out(), "Global right-hand-side vector is zero vector."); + case analysis_empty_entry: + return std::format_to(ctx.out(), "Current entry has no entry points."); case reader_file_fail_to_open: return std::format_to(ctx.out(), "Reader: Failed to open the file."); case reader_uninitialized: diff --git a/source/centipede/util/return_types.hpp b/source/centipede/util/return_types.hpp index 7ac79ee4..520da800 100644 --- a/source/centipede/util/return_types.hpp +++ b/source/centipede/util/return_types.hpp @@ -10,6 +10,9 @@ namespace centipede { + /** + * @brief Template alias for expected string results. + */ using StrError = std::expected; /** @@ -17,4 +20,9 @@ namespace centipede */ template using EnumError = std::expected; + + /** + * @brief Template alias for expected void results. + */ + using VoidError = std::expected; } // namespace centipede diff --git a/source/centipede/writer/binary.hpp b/source/centipede/writer/binary.hpp index d5550c85..a0a1f66e 100644 --- a/source/centipede/writer/binary.hpp +++ b/source/centipede/writer/binary.hpp @@ -1,6 +1,6 @@ #pragma once -#include "centipede/data/entry.hpp" +#include "centipede/data/entrypoint.hpp" #include "centipede/util/common_definitions.hpp" #include "centipede/util/error_types.hpp" #include "centipede/util/return_types.hpp" diff --git a/test/integration_tests/CMakeLists.txt b/test/integration_tests/CMakeLists.txt index b67ea9ee..bc860a0c 100644 --- a/test/integration_tests/CMakeLists.txt +++ b/test/integration_tests/CMakeLists.txt @@ -8,17 +8,14 @@ set(COMPILE_OPTIONS -fno-rtti ) +#===================================writer=================================== + add_executable(integration_test_writer test_writer.cpp) target_link_libraries(integration_test_writer PRIVATE centipede::centipede) target_compile_options(integration_test_writer PRIVATE ${COMPILE_OPTIONS}) -# if(ENABLE_COVERAGE) -# target_compile_options(integration_test_writer PRIVATE --coverage) -# target_link_options(integration_test_writer PRIVATE --coverage) -# endif() - add_test(NAME integration_test_writer_data_gen COMMAND integration_test_writer) set_tests_properties(integration_test_writer_data_gen PROPERTIES TIMEOUT 10) @@ -34,6 +31,8 @@ set_tests_properties( PROPERTIES DEPENDS integration_test_writer_data_gen ) +#===================================reader=============================== + add_executable(integration_test_reader test_reader.cpp) add_test(NAME integration_test_reader_data_read COMMAND integration_test_reader) @@ -46,3 +45,28 @@ set_tests_properties( integration_test_reader_data_read PROPERTIES DEPENDS integration_test_writer_data_gen ) +#===================================eigen core=============================== + +find_package(mps) + +if(NOT Boost_FOUND) + message(STATUS "Boost library is not used for the build!") +endif() + +if(mps_FOUND AND Boost_FOUND) + message(STATUS "Testing with simulated data from mps is enabled!") + add_executable(test_eigen test_eigen.cpp) + + target_link_libraries(test_eigen PRIVATE centipede::centipede mps::mps Boost::headers) + + target_compile_options(test_eigen PRIVATE ${COMPILE_OPTIONS}) + + target_compile_definitions(test_eigen PRIVATE BOOST_NO_EXCEPTIONS BOOST_THROW_EXCEPTION_HEADER_ONLY) + + add_test(NAME integration_test_eigen COMMAND test_eigen) + set_tests_properties(integration_test_eigen PROPERTIES TIMEOUT 10) +else() + message(STATUS "Testing with simulated data from mps is disabled!") +endif() + +#============================================================================== diff --git a/test/integration_tests/test_eigen.cpp b/test/integration_tests/test_eigen.cpp new file mode 100644 index 00000000..794600f4 --- /dev/null +++ b/test/integration_tests/test_eigen.cpp @@ -0,0 +1,105 @@ +#include "centipede/centipede.hpp" +#include "centipede/data/entrypoint.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ + enum class ResStatus + { + succeed, + fail, + }; +} // namespace + +namespace bh = boost::histogram; + +namespace boost +{ + + BOOST_NORETURN void throw_exception(std::exception const&) { std::abort(); } + + BOOST_NORETURN void throw_exception(std::exception const&, boost::source_location const&) { std::abort(); } + +} // namespace boost + +auto main() -> int +{ + const auto config = mps::MPSConfig{}; + + auto mps = mps::MPS{ config }; + mps.init(); + + const auto n_globals = config.detector.spec.num_modules * 2; + auto handler = centipede::Handler{ n_globals }; + + auto hist_residuals = bh::make_histogram(bh::axis::regular{ 20, 0., 17., "chi2" }); + auto hist_pvalues = bh::make_histogram(bh::axis::regular{ 20, 0., 1., "pvalues" }); + + auto entry_point_input = centipede::EntryPoint{}; + for (auto _ : std::views::iota(0, 1000)) + { + auto entry_span = mps.generate_one_entry_data(); + + for (const auto& entrypoint : entry_span) + { + entry_point_input.reset(); + entry_point_input.set_measurement(entrypoint.measurement) + // entry_point_input.set_measurement(0.) + .set_sigma(entrypoint.sigma) + .set_globals(entrypoint.globals | + std::views::transform([](const auto& globals) + { return std::pair{ globals.first - 1, globals.second }; })) + .set_locals(entrypoint.locals); + auto res = handler.add_entrypoint(entry_point_input); + if (not res) + { + std::println(stderr, "Error from adding the current point: {}", res.error()); + } + } + // std::println("Number of entrypoints in the current entry: {}", handler.get_current_state().next_point_index); + auto res = handler.analyze_current_entry(); + const auto& state = handler.get_current_entry_state(); + hist_residuals(state.chi2); + hist_pvalues(state.p_value); + if (not res) + { + if (res.error() != centipede::ErrorCode::analysis_empty_entry) + { + std::println(stderr, "Error from analyzing the current entry: {}", res.error()); + } + } + } + + std::println("Solving the equations"); + auto is_ok = handler.solve(); + + if (not is_ok) + { + std::println("Error from the solving: {}", is_ok.error()); + return EXIT_FAILURE; + } + + const auto& result = handler.get_result(); + + auto hist_ostream = std::ostringstream{}; + hist_ostream << hist_residuals << hist_pvalues; + + std::println("{}", hist_ostream.str()); + std::println("{}", result); + std::println("parameters: {}", result.parameters); + std::println("eigen values: {}", result.eigen_values); + + return EXIT_SUCCESS; +} diff --git a/test/integration_tests/test_writer.cpp b/test/integration_tests/test_writer.cpp index 25bcfa62..42692bca 100644 --- a/test/integration_tests/test_writer.cpp +++ b/test/integration_tests/test_writer.cpp @@ -1,4 +1,4 @@ -#include "centipede/data/entry.hpp" +#include "centipede/data/entrypoint.hpp" #include "centipede/writer/binary.hpp" #include // IWYU pragma: keep #include diff --git a/test/unit_tests/shared.hpp b/test/unit_tests/shared.hpp index 3623aa82..0de57799 100644 --- a/test/unit_tests/shared.hpp +++ b/test/unit_tests/shared.hpp @@ -1,6 +1,6 @@ #pragma once -#include "centipede/data/entry.hpp" +#include "centipede/data/entrypoint.hpp" #include #include #include diff --git a/test/unit_tests/test_base_engine.cpp b/test/unit_tests/test_base_engine.cpp index 320ea778..464d9f5a 100644 --- a/test/unit_tests/test_base_engine.cpp +++ b/test/unit_tests/test_base_engine.cpp @@ -27,10 +27,10 @@ namespace centipede::test MOCK_METHOD(void, fill_measurements, (const std::vector&), ()); MOCK_METHOD(void, fill_sigmas, (const std::vector&), ()); MOCK_METHOD(void, fill_local_derivs, (const std::vector::Deriv>&), ()); - MOCK_METHOD(void, fill_global_derivs, (const std::vector::Deriv>&), ()); + MOCK_METHOD(VoidError, fill_global_derivs, (const std::vector::Deriv>&), ()); // Called in analyze method. - MOCK_METHOD((EnumError<>), fit_local_pars, (), ()); + MOCK_METHOD((VoidError), fit_local_pars, (), ()); MOCK_METHOD((std::pair), calculate_local_fit_chi_square, (), ()); MOCK_METHOD((void), update_global_factor_matrix, (), ()); MOCK_METHOD((void), update_global_rhs_vector, (), ()); @@ -80,7 +80,7 @@ namespace centipede::test EXPECT_CALL(engine, fill_local_derivs(entry.local_derivs)).Times(1); EXPECT_CALL(engine, fill_global_derivs(entry.global_derivs)).Times(1); - engine.fill_data(entry); + EXPECT_TRUE_RES(engine.fill_data(entry)); EXPECT_EQ(engine.get_log().n_entries_read, 1); } @@ -95,7 +95,7 @@ namespace centipede::test EXPECT_CALL(engine, fill_local_derivs(entry.local_derivs)).Times(0); EXPECT_CALL(engine, fill_global_derivs(entry.global_derivs)).Times(0); - engine.fill_data(entry); + EXPECT_TRUE_RES(engine.fill_data(entry)); EXPECT_EQ(engine.get_log().n_entries_read, 0); } diff --git a/test/unit_tests/test_eigen_engine.cpp b/test/unit_tests/test_eigen_engine.cpp index 8c815ce8..48a3099f 100644 --- a/test/unit_tests/test_eigen_engine.cpp +++ b/test/unit_tests/test_eigen_engine.cpp @@ -13,7 +13,7 @@ namespace centipede::test TEST(eigen_engine, constructor) { constexpr auto n_global_pars = 10; - auto engine = core::engine::Engine{ n_global_pars }; + auto engine = core::engine::Engine{ n_global_pars }; const auto& factor_matrix = engine.get_global_factor_matrix(); EXPECT_EQ(factor_matrix.rows(), n_global_pars); EXPECT_EQ(factor_matrix.cols(), n_global_pars); @@ -24,7 +24,7 @@ namespace centipede::test TEST(eigen_engine, solve_rank_deficit) { - using EngineClass = core::engine::Engine; + using EngineClass = core::engine::Engine; auto result = Result{}; const auto globals = []() @@ -48,7 +48,7 @@ namespace centipede::test TEST(eigen_engine, solve_negative_definite) { - using EngineClass = core::engine::Engine; + using EngineClass = core::engine::Engine; auto result = Result{}; const auto globals = []() @@ -69,7 +69,7 @@ namespace centipede::test TEST(eigen_engine, solve_zero_factor_matrix) { - using EngineClass = core::engine::Engine; + using EngineClass = core::engine::Engine; auto result = Result{}; const auto globals = []() @@ -90,7 +90,7 @@ namespace centipede::test TEST(eigen_engine, solve_zero_rhs_vector) { - using EngineClass = core::engine::Engine; + using EngineClass = core::engine::Engine; auto result = Result{}; const auto globals = []() @@ -111,7 +111,7 @@ namespace centipede::test TEST(eigen_engine, solve) { - using EngineClass = core::engine::Engine; + using EngineClass = core::engine::Engine; auto result = Result{}; const auto globals = []() diff --git a/test/unit_tests/test_handler.cpp b/test/unit_tests/test_handler.cpp index ac7e32f6..89a1e065 100644 --- a/test/unit_tests/test_handler.cpp +++ b/test/unit_tests/test_handler.cpp @@ -6,41 +6,32 @@ namespace { constexpr auto DEFAULT_MAX_POINT = 30; - using centipede::core::Handler; + using centipede::Handler; } // namespace namespace centipede::test { - TEST(handler, constructor) { auto handler = Handler{}; } + TEST(handler, constructor) { auto handler = Handler{ DEFAULT_MAX_GLOBAL_ID }; } TEST(handler, constructor_float_eigen) { using HandlerType = Handler; using Config = HandlerType::Config; - auto handler = HandlerType{ Config{ .n_globals = DEFAULT_MAX_GLOBAL_ID } }; + auto handler = HandlerType{ DEFAULT_MAX_GLOBAL_ID }; } TEST(handler, constructor_double_eigen) { using HandlerType = Handler; using Config = HandlerType::Config; - auto handler = HandlerType{ Config{ .n_globals = DEFAULT_MAX_GLOBAL_ID } }; - } - - TEST(handler, init) - { - auto handler = Handler{}; - auto err = handler.init(); - EXPECT_TRUE(err.has_value()); + auto handler = HandlerType{ DEFAULT_MAX_GLOBAL_ID }; } // NOLINTBEGIN(readability-function-cognitive-complexity) TEST(handler, add_entrypoint) { - auto handler = Handler{}; - auto err = handler.init(); - EXPECT_TRUE(err.has_value()); + auto handler = Handler{ DEFAULT_MAX_GLOBAL_ID }; constexpr auto n_points = 100; const auto entrypoints = generate_random_entry_points(n_points); EXPECT_EQ(n_points, entrypoints.size()); @@ -62,9 +53,7 @@ namespace centipede::test TEST(handler, local_derivs_incomp_numbers) { - auto handler = Handler{}; - auto err = handler.init(); - EXPECT_TRUE(err.has_value()); + auto handler = Handler{ DEFAULT_MAX_GLOBAL_ID }; constexpr auto n_points = 10; const auto entrypoints = generate_random_entry_points(n_points); for (const auto& entry_point : entrypoints) diff --git a/test/unit_tests/test_master_engine.cpp b/test/unit_tests/test_master_engine.cpp index a19be642..46c51cb8 100644 --- a/test/unit_tests/test_master_engine.cpp +++ b/test/unit_tests/test_master_engine.cpp @@ -16,7 +16,7 @@ namespace centipede::core::engine { namespace { - struct Globals + struct GlobalsType { }; @@ -25,23 +25,23 @@ namespace centipede::core::engine { public: MOCK_METHOD(void, construct_with, (std::size_t n_globals), ()); - MOCK_METHOD(void, solve, (const Globals& globals, Result& result), ()); + MOCK_METHOD(void, solve, (const GlobalsType& globals, Result& result), ()); }; } // namespace /** @cond */ template - class Engine + class Engine { public: Engine(std::size_t n_globals) { mock_helper->construct_with(n_globals); } - using Globals = Globals; + using Globals = GlobalsType; static void solve(const Globals& globals, Result& result) { mock_helper->solve(globals, result); } MOCK_METHOD(void, add_to_globals, (Globals & globals), (const)); - MOCK_METHOD((void), fill_data, (const Entry& entry), (const)); + MOCK_METHOD((VoidError), fill_data, (const Entry& entry), (const)); MOCK_METHOD((void), add_to_result, (Result & result), (const)); MOCK_METHOD((EnumError<>), analyze, (double alpha), (const)); @@ -49,7 +49,7 @@ namespace centipede::core::engine }; template - MockHelper* Engine::mock_helper = nullptr; + MockHelper* Engine::mock_helper = nullptr; /** @endcond */ } // namespace centipede::core::engine @@ -60,13 +60,13 @@ namespace centipede::test class master_engine : public testing::Test { public: - using Master = engine::Master; + using Master = engine::Master; using Config = Master::Config; protected: using MockHelperType = core::engine::MockHelper; using EngineClass = Master::EngineImp; - using Result = Result; + using ResultType = Result; void SetUp() override { @@ -94,7 +94,7 @@ namespace centipede::test .set_globals(std::pair{ 9, 2.F }, std::pair{ 2, 3.F }) .set_locals(1.5F, 2.5F); const auto& state = master_->get_current_state(); - EXPECT_EQ(state.point_index, 0); + EXPECT_EQ(state.next_point_index, 0); EXPECT_FALSE(state.entry.n_locals); EXPECT_EQ(state.entry.measurements.size(), 0); EXPECT_EQ(state.entry.sigmas.size(), 0); @@ -103,7 +103,7 @@ namespace centipede::test EXPECT_TRUE_RES(master_->add_entrypoint(entrypoint)); - EXPECT_EQ(state.point_index, 1); + EXPECT_EQ(state.next_point_index, 1); ASSERT_TRUE(state.entry.n_locals); EXPECT_EQ(state.entry.n_locals.value(), 2); EXPECT_EQ(state.entry.measurements.size(), 1); @@ -158,7 +158,7 @@ namespace centipede::test EXPECT_TRUE_RES(master_->analyze()); const auto& state = master_->get_current_state(); - EXPECT_EQ(state.point_index, 0); + EXPECT_EQ(state.next_point_index, 0); EXPECT_FALSE(state.entry.n_locals); EXPECT_EQ(state.entry.measurements.size(), 0); EXPECT_EQ(state.entry.sigmas.size(), 0); @@ -184,7 +184,7 @@ namespace centipede::test EXPECT_CALL(*engine_class_, add_to_result(testing::_)).Times(1); EXPECT_CALL(*mock_helper_, solve(testing::_, testing::_)) .Times(1) - .WillOnce([](const auto&, Result& result) { result.error_status = ErrorCode::success; }); + .WillOnce([](const auto&, ResultType& result) { result.error_status = ErrorCode::success; }); EXPECT_TRUE_RES(master_->solve()); } @@ -195,7 +195,7 @@ namespace centipede::test EXPECT_CALL(*engine_class_, add_to_result(testing::_)).Times(1); EXPECT_CALL(*mock_helper_, solve(testing::_, testing::_)) .Times(1) - .WillOnce([](const auto&, Result& result) { result.error_status = ErrorCode::analysis_rank_deficit; }); + .WillOnce([](const auto&, ResultType& result) { result.error_status = ErrorCode::analysis_rank_deficit; }); auto res = master_->solve(); EXPECT_FALSE(res); diff --git a/util/mathjax-config.js b/util/mathjax-config.js index c8212065..50764287 100644 --- a/util/mathjax-config.js +++ b/util/mathjax-config.js @@ -1,5 +1,3 @@ -window.MathJax = { - tex: { - tags: "ams", - }, -}; +window.MathJax = window.MathJax || {}; +window.MathJax.tex = window.MathJax.tex || {}; +window.MathJax.tex.tags = "ams";