From e273be1be2aa1af46964d560b4ba5d6b6493a831 Mon Sep 17 00:00:00 2001 From: Yanzhao Wang Date: Fri, 8 May 2026 17:37:00 +0200 Subject: [PATCH 1/9] add integration test for eigen engine --- source/centipede/core/engines/master_engine.hpp | 4 ++-- test/integration_tests/CMakeLists.txt | 8 ++++++++ test/integration_tests/test_eigen.cpp | 12 ++++++++++++ 3 files changed, 22 insertions(+), 2 deletions(-) create mode 100644 test/integration_tests/test_eigen.cpp diff --git a/source/centipede/core/engines/master_engine.hpp b/source/centipede/core/engines/master_engine.hpp index bfd0ac01..86571db3 100644 --- a/source/centipede/core/engines/master_engine.hpp +++ b/source/centipede/core/engines/master_engine.hpp @@ -48,7 +48,7 @@ namespace centipede::core::engine Entry entry; //!< Data storing the current entry. }; - using Result = Result; + using ResultType = Result; using EngineImp = Engine; using DataTypeUsed = DataType; @@ -153,7 +153,7 @@ namespace centipede::core::engine private: Config config_; - Result result_; + ResultType result_; State current_state_; EngineImp engine_imp_{}; EngineImp::Globals globals_{}; diff --git a/test/integration_tests/CMakeLists.txt b/test/integration_tests/CMakeLists.txt index b67ea9ee..c12aa85a 100644 --- a/test/integration_tests/CMakeLists.txt +++ b/test/integration_tests/CMakeLists.txt @@ -14,11 +14,19 @@ target_link_libraries(integration_test_writer PRIVATE centipede::centipede) target_compile_options(integration_test_writer PRIVATE ${COMPILE_OPTIONS}) +add_executable(test_eigen test_eigen.cpp) + +target_link_libraries(test_eigen PRIVATE centipede::centipede) + +target_compile_options(test_eigen 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_eigen COMMAND test_eigen) +set_tests_properties(integration_test_eigen PROPERTIES TIMEOUT 10) + add_test(NAME integration_test_writer_data_gen COMMAND integration_test_writer) set_tests_properties(integration_test_writer_data_gen PROPERTIES TIMEOUT 10) diff --git a/test/integration_tests/test_eigen.cpp b/test/integration_tests/test_eigen.cpp new file mode 100644 index 00000000..82a30323 --- /dev/null +++ b/test/integration_tests/test_eigen.cpp @@ -0,0 +1,12 @@ +#include "centipede/centipede.hpp" +#include +#include + +void entrypoints_generator() {} + +auto main() -> int +{ + auto entries = std::vector>{}; + + return EXIT_SUCCESS; +} From 327d0bbfef660afa0419ad0ebb4906bae4c17316 Mon Sep 17 00:00:00 2001 From: Yanzhao Wang Date: Fri, 8 May 2026 17:49:07 +0200 Subject: [PATCH 2/9] fix warnings --- test/unit_tests/test_master_engine.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/unit_tests/test_master_engine.cpp b/test/unit_tests/test_master_engine.cpp index a19be642..58dc3e68 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,7 +25,7 @@ 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 @@ -36,7 +36,7 @@ namespace centipede::core::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); } @@ -66,7 +66,7 @@ namespace centipede::test protected: using MockHelperType = core::engine::MockHelper; using EngineClass = Master::EngineImp; - using Result = Result; + using ResultType = Result; void SetUp() override { @@ -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); From c5c8dda23c21acf14986c5329482122600cf748b Mon Sep 17 00:00:00 2001 From: Yanzhao Wang Date: Fri, 8 May 2026 17:54:47 +0200 Subject: [PATCH 3/9] switch to CMakeConfigDeps --- conanfile.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/conanfile.py b/conanfile.py index 662b4bb5..4d2f56a6 100644 --- a/conanfile.py +++ b/conanfile.py @@ -8,7 +8,7 @@ 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 From e592568ebf30edb5ebc62419738b39c42e3c6fd3 Mon Sep 17 00:00:00 2001 From: Yanzhao Wang Date: Sat, 16 May 2026 15:42:42 +0200 Subject: [PATCH 4/9] Fix compile errors --- source/centipede/CMakeLists.txt | 4 ++- source/centipede/core/CMakeLists.txt | 2 -- source/centipede/core/engines/base_engine.hpp | 2 +- .../centipede/core/engines/eigen_engine.hpp | 13 +++---- .../centipede/core/engines/engine_concept.hpp | 2 +- .../centipede/core/engines/engine_types.hpp | 6 ++-- source/centipede/core/handler.cpp | 0 source/centipede/core/handler.hpp | 8 ++--- source/centipede/util/return_types.hpp | 8 +++++ test/integration_tests/CMakeLists.txt | 36 ++++++++++++------- test/integration_tests/test_eigen.cpp | 7 +++- test/unit_tests/test_eigen_engine.cpp | 12 +++---- test/unit_tests/test_handler.cpp | 2 +- test/unit_tests/test_master_engine.cpp | 6 ++-- 14 files changed, 66 insertions(+), 42 deletions(-) delete mode 100644 source/centipede/core/handler.cpp diff --git a/source/centipede/CMakeLists.txt b/source/centipede/CMakeLists.txt index 7d814919..e71158f6 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,8 @@ target_sources( FILES centipede.hpp ) +target_link_libraries(core PUBLIC Eigen3::Eigen GSL::gsl) + 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..2146ea33 100644 --- a/source/centipede/core/engines/base_engine.hpp +++ b/source/centipede/core/engines/base_engine.hpp @@ -178,7 +178,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..176215a9 100644 --- a/source/centipede/core/engines/eigen_engine.hpp +++ b/source/centipede/core/engines/eigen_engine.hpp @@ -26,7 +26,7 @@ namespace centipede::core::engine * @brief Engine template specialization for Eigen library implementation. */ template - class Engine : public Base + class Engine : public Base { public: /** @@ -192,7 +192,7 @@ namespace centipede::core::engine global_t_.setFromSortedTriplets(triplets_.begin(), triplets_.end()); } - 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); @@ -234,7 +234,7 @@ namespace centipede::core::engine 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. @@ -257,7 +257,7 @@ namespace centipede::core::engine 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 +274,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(); } diff --git a/source/centipede/core/engines/engine_concept.hpp b/source/centipede/core/engines/engine_concept.hpp index a13c5658..a8e8bbdb 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) { 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/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..1a99c8ba 100644 --- a/source/centipede/core/handler.hpp +++ b/source/centipede/core/handler.hpp @@ -6,7 +6,7 @@ #include "centipede/util/return_types.hpp" #include -namespace centipede::core +namespace centipede { /** @@ -14,7 +14,7 @@ namespace centipede::core * * This class should handle all inputs and configurations from users */ - template + template class Handler { public: @@ -26,7 +26,7 @@ namespace centipede::core { std::size_t n_globals = 0; //!< Number of global parameters. }; - using EngineType = engine::Master; + using EngineType = core::engine::Master; /** * @brief Constructor @@ -61,4 +61,4 @@ namespace centipede::core // std::vector<> }; -} // namespace centipede::core +} // namespace centipede 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/test/integration_tests/CMakeLists.txt b/test/integration_tests/CMakeLists.txt index c12aa85a..e934f0d2 100644 --- a/test/integration_tests/CMakeLists.txt +++ b/test/integration_tests/CMakeLists.txt @@ -8,25 +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}) -add_executable(test_eigen test_eigen.cpp) - -target_link_libraries(test_eigen PRIVATE centipede::centipede) - -target_compile_options(test_eigen 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_eigen COMMAND test_eigen) -set_tests_properties(integration_test_eigen PROPERTIES TIMEOUT 10) - add_test(NAME integration_test_writer_data_gen COMMAND integration_test_writer) set_tests_properties(integration_test_writer_data_gen PROPERTIES TIMEOUT 10) @@ -42,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) @@ -54,3 +45,22 @@ set_tests_properties( integration_test_reader_data_read PROPERTIES DEPENDS integration_test_writer_data_gen ) +#===================================eigen core=============================== + +find_package(mps) + +if(mps_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) + + target_compile_options(test_eigen PRIVATE ${COMPILE_OPTIONS}) + + 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 index 82a30323..8081e64a 100644 --- a/test/integration_tests/test_eigen.cpp +++ b/test/integration_tests/test_eigen.cpp @@ -2,11 +2,16 @@ #include #include -void entrypoints_generator() {} +namespace +{ + // void entrypoints_generator() {} +} // namespace auto main() -> int { auto entries = std::vector>{}; + auto handler = centipede::Handler{}; + return EXIT_SUCCESS; } 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..fdeb0cd6 100644 --- a/test/unit_tests/test_handler.cpp +++ b/test/unit_tests/test_handler.cpp @@ -6,7 +6,7 @@ namespace { constexpr auto DEFAULT_MAX_POINT = 30; - using centipede::core::Handler; + using centipede::Handler; } // namespace diff --git a/test/unit_tests/test_master_engine.cpp b/test/unit_tests/test_master_engine.cpp index 58dc3e68..b8eade34 100644 --- a/test/unit_tests/test_master_engine.cpp +++ b/test/unit_tests/test_master_engine.cpp @@ -32,7 +32,7 @@ namespace centipede::core::engine /** @cond */ template - class Engine + class Engine { public: Engine(std::size_t n_globals) { mock_helper->construct_with(n_globals); } @@ -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,7 +60,7 @@ namespace centipede::test class master_engine : public testing::Test { public: - using Master = engine::Master; + using Master = engine::Master; using Config = Master::Config; protected: From 4e4f1a2b7e5768f9186faa2c5b6f87ac3baa96f4 Mon Sep 17 00:00:00 2001 From: Yanzhao Wang Date: Tue, 19 May 2026 23:30:00 +0200 Subject: [PATCH 5/9] add entrypoint.hpp --- source/centipede/core/engines/base_engine.hpp | 2 +- .../centipede/core/engines/eigen_engine.hpp | 16 +- .../centipede/core/engines/master_engine.hpp | 14 +- source/centipede/core/engines/result.hpp | 5 +- source/centipede/core/handler.hpp | 25 +- source/centipede/data/entry.hpp | 178 +----------- source/centipede/data/entrypoint.hpp | 274 ++++++++++++++++++ .../{entry_base.hpp => entrypoint_base.hpp} | 0 source/centipede/writer/binary.hpp | 2 +- test/integration_tests/test_eigen.cpp | 45 ++- test/integration_tests/test_writer.cpp | 2 +- test/unit_tests/shared.hpp | 2 +- test/unit_tests/test_handler.cpp | 11 - 13 files changed, 359 insertions(+), 217 deletions(-) create mode 100644 source/centipede/data/entrypoint.hpp rename source/centipede/data/{entry_base.hpp => entrypoint_base.hpp} (100%) diff --git a/source/centipede/core/engines/base_engine.hpp b/source/centipede/core/engines/base_engine.hpp index 2146ea33..1be1b3ad 100644 --- a/source/centipede/core/engines/base_engine.hpp +++ b/source/centipede/core/engines/base_engine.hpp @@ -137,7 +137,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; diff --git a/source/centipede/core/engines/eigen_engine.hpp b/source/centipede/core/engines/eigen_engine.hpp index 176215a9..39a6a3c4 100644 --- a/source/centipede/core/engines/eigen_engine.hpp +++ b/source/centipede/core/engines/eigen_engine.hpp @@ -136,9 +136,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 +163,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()); } @@ -221,7 +223,7 @@ namespace centipede::core::engine Eigen::internal::set_is_malloc_allowed(false); 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) { @@ -297,7 +299,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/master_engine.hpp b/source/centipede/core/engines/master_engine.hpp index 86571db3..0b37098f 100644 --- a/source/centipede/core/engines/master_engine.hpp +++ b/source/centipede/core/engines/master_engine.hpp @@ -6,7 +6,8 @@ #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/data/entrypoint_base.hpp" #include "centipede/util/error_types.hpp" #include "centipede/util/return_types.hpp" #include @@ -114,7 +115,7 @@ namespace centipede::core::engine /** * @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<> @@ -135,17 +136,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_; } 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.hpp b/source/centipede/core/handler.hpp index 1a99c8ba..8674ce07 100644 --- a/source/centipede/core/handler.hpp +++ b/source/centipede/core/handler.hpp @@ -2,7 +2,7 @@ #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/return_types.hpp" #include @@ -39,11 +39,13 @@ namespace centipede { } - /** - * @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<> @@ -51,10 +53,19 @@ namespace centipede 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(); + + return engine_.analyze().transform([n_points]() { return n_points; }); + }; [[nodiscard]] auto get_current_state() const -> const auto& { return 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_; 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..43782c71 --- /dev/null +++ b/source/centipede/data/entrypoint.hpp @@ -0,0 +1,274 @@ +#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 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 self The dynamic reference to the caller. + * @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 self The dynamic reference to the caller. + * @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. + * @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; + + using internal::EntryPointBase::set_globals; + using internal::EntryPointBase::set_locals; + + /** + * @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 #NGlobals. + * + * @tparam View Types of the view. + * @param self The dynamic reference to the caller. + * @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 + * #NLocals. + * + * @tparam View Types of the view. + * @param self The dynamic reference to the caller. + * @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 100% rename from source/centipede/data/entry_base.hpp rename to source/centipede/data/entrypoint_base.hpp 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/test_eigen.cpp b/test/integration_tests/test_eigen.cpp index 8081e64a..62f57a3d 100644 --- a/test/integration_tests/test_eigen.cpp +++ b/test/integration_tests/test_eigen.cpp @@ -1,5 +1,10 @@ #include "centipede/centipede.hpp" +#include "centipede/data/entrypoint.hpp" #include +#include +#include +#include +#include #include namespace @@ -9,9 +14,45 @@ namespace auto main() -> int { - auto entries = std::vector>{}; + const auto config = mps::MPSConfig{}; - auto handler = centipede::Handler{}; + auto mps = mps::MPS{ config }; + mps.init(); + + auto handler = centipede::Handler{}; + + 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) + .set_sigma(entrypoint.sigma) + .set_globals(entrypoint.globals) + .set_locals(entrypoint.locals); + } + auto res = handler.analyze_current_entry(); + if (not res) + { + std::println("Error: {}", res.error()); + return EXIT_FAILURE; + } + } + + auto is_ok = handler.solve(); + + if (not is_ok) + { + std::println("Error: {}", is_ok.error()); + return EXIT_FAILURE; + } + + const auto& result = handler.get_result(); + + std::println("{}", result); 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_handler.cpp b/test/unit_tests/test_handler.cpp index fdeb0cd6..9b80afd7 100644 --- a/test/unit_tests/test_handler.cpp +++ b/test/unit_tests/test_handler.cpp @@ -28,19 +28,10 @@ namespace centipede::test 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()); - } - // NOLINTBEGIN(readability-function-cognitive-complexity) TEST(handler, add_entrypoint) { auto handler = Handler{}; - auto err = handler.init(); - EXPECT_TRUE(err.has_value()); constexpr auto n_points = 100; const auto entrypoints = generate_random_entry_points(n_points); EXPECT_EQ(n_points, entrypoints.size()); @@ -63,8 +54,6 @@ namespace centipede::test TEST(handler, local_derivs_incomp_numbers) { auto handler = Handler{}; - auto err = handler.init(); - EXPECT_TRUE(err.has_value()); constexpr auto n_points = 10; const auto entrypoints = generate_random_entry_points(n_points); for (const auto& entry_point : entrypoints) From 54716f3ac3fc1e7aae999276d19cc740104565c0 Mon Sep 17 00:00:00 2001 From: Yanzhao Wang Date: Wed, 20 May 2026 00:02:17 +0200 Subject: [PATCH 6/9] fix rebase --- source/centipede/reader/binary.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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" From 9e02a765ac2ceda22388ac0a0a8ce2de95045705 Mon Sep 17 00:00:00 2001 From: Yanzhao Wang Date: Fri, 5 Jun 2026 16:30:18 +0200 Subject: [PATCH 7/9] add libassert and fix errors --- cmake/load_dependencies.cmake | 1 + conanfile.py | 1 + source/centipede/CMakeLists.txt | 6 + source/centipede/core/engines/base_engine.hpp | 14 ++- .../centipede/core/engines/eigen_engine.hpp | 61 ++++++++-- .../centipede/core/engines/engine_concept.hpp | 4 +- .../centipede/core/engines/master_engine.hpp | 113 ++++++++++-------- source/centipede/core/handler.hpp | 15 ++- source/centipede/data/entrypoint.hpp | 28 +++++ source/centipede/util/common_definitions.hpp | 3 + source/centipede/util/error_types.hpp | 19 ++- test/integration_tests/test_eigen.cpp | 27 ++++- test/unit_tests/test_base_engine.cpp | 8 +- test/unit_tests/test_handler.cpp | 10 +- test/unit_tests/test_master_engine.cpp | 8 +- 15 files changed, 221 insertions(+), 97 deletions(-) diff --git a/cmake/load_dependencies.cmake b/cmake/load_dependencies.cmake index bd67d89a..e4f68e79 100644 --- a/cmake/load_dependencies.cmake +++ b/cmake/load_dependencies.cmake @@ -3,6 +3,7 @@ find_package(magic_enum REQUIRED CONFIG) find_package(CLI11 REQUIRED CONFIG) find_package(Eigen3 REQUIRED CONFIG) find_package(GSL REQUIRED) +find_package(libassert) if(ENABLE_TEST) find_package(GTest CONFIG REQUIRED) diff --git a/conanfile.py b/conanfile.py index 4d2f56a6..c7b01f44 100644 --- a/conanfile.py +++ b/conanfile.py @@ -15,6 +15,7 @@ def requirements(self): 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/source/centipede/CMakeLists.txt b/source/centipede/CMakeLists.txt index e71158f6..a6d195df 100644 --- a/source/centipede/CMakeLists.txt +++ b/source/centipede/CMakeLists.txt @@ -24,6 +24,12 @@ target_sources( 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/engines/base_engine.hpp b/source/centipede/core/engines/base_engine.hpp index 1be1b3ad..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> diff --git a/source/centipede/core/engines/eigen_engine.hpp b/source/centipede/core/engines/eigen_engine.hpp index 39a6a3c4..3bef76e8 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,8 +21,25 @@ #include #include +#ifdef HAS_LIBASSERT +#include "libassert/assert.hpp" +#else +#include +#endif + namespace centipede::core::engine { + + namespace + { + class EigenMemGuard + { + public: + EigenMemGuard() { Eigen::internal::set_is_malloc_allowed(false); } + ~EigenMemGuard() { Eigen::internal::set_is_malloc_allowed(true); } + }; + } // namespace + /** * @brief Engine template specialization for Eigen library implementation. */ @@ -55,7 +73,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()) { @@ -91,6 +115,17 @@ namespace centipede::core::engine 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(); } @@ -175,13 +210,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(); @@ -192,12 +232,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() -> 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()); @@ -213,14 +254,12 @@ 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 - static_cast(local_size); @@ -232,7 +271,6 @@ 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 }; } @@ -253,7 +291,14 @@ 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 {}; diff --git a/source/centipede/core/engines/engine_concept.hpp b/source/centipede/core/engines/engine_concept.hpp index a8e8bbdb..299b7c42 100644 --- a/source/centipede/core/engines/engine_concept.hpp +++ b/source/centipede/core/engines/engine_concept.hpp @@ -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/master_engine.hpp b/source/centipede/core/engines/master_engine.hpp index 0b37098f..51157649 100644 --- a/source/centipede/core/engines/master_engine.hpp +++ b/source/centipede/core/engines/master_engine.hpp @@ -7,7 +7,6 @@ #include "centipede/core/engines/result.hpp" #include "centipede/data/entry.hpp" #include "centipede/data/entrypoint.hpp" -#include "centipede/data/entrypoint_base.hpp" #include "centipede/util/error_types.hpp" #include "centipede/util/return_types.hpp" #include @@ -45,8 +44,8 @@ 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 ResultType = Result; @@ -70,46 +69,37 @@ 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; + }); } /** @@ -118,16 +108,16 @@ namespace centipede::core::engine * 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; + }); } /** @@ -156,17 +146,36 @@ namespace centipede::core::engine 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/handler.hpp b/source/centipede/core/handler.hpp index 8674ce07..483c9c55 100644 --- a/source/centipede/core/handler.hpp +++ b/source/centipede/core/handler.hpp @@ -3,8 +3,10 @@ #include "centipede/core/engines/engine_types.hpp" #include "centipede/core/engines/master_engine.hpp" #include "centipede/data/entrypoint.hpp" +#include "centipede/util/error_types.hpp" #include "centipede/util/return_types.hpp" #include +#include namespace centipede { @@ -24,7 +26,7 @@ namespace centipede */ struct Config { - std::size_t n_globals = 0; //!< Number of global parameters. + // std::size_t n_globals = 0; //!< Number of global parameters. }; using EngineType = core::engine::Master; @@ -33,9 +35,9 @@ namespace centipede * * 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 } } { } @@ -48,7 +50,7 @@ namespace centipede // 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); } @@ -57,6 +59,11 @@ namespace centipede { 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; }); }; diff --git a/source/centipede/data/entrypoint.hpp b/source/centipede/data/entrypoint.hpp index 43782c71..e20785be 100644 --- a/source/centipede/data/entrypoint.hpp +++ b/source/centipede/data/entrypoint.hpp @@ -37,6 +37,20 @@ namespace centipede */ EntryPoint() = default; + /** + * @brief Getter for local parameter sizes. + * @param self The dynamic reference to the caller. + * @return Sizes of local parameters. + */ + [[nodiscard]] inline auto get_n_locals() const -> std::size_t { return locals_.size(); } + + /** + * @brief Getter for global parameter sizes. + * @param self The dynamic reference to the caller. + * @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. @@ -184,6 +198,20 @@ namespace centipede using internal::EntryPointBase::set_globals; using internal::EntryPointBase::set_locals; + /** + * @brief Getter for local parameter sizes. + * @param self The dynamic reference to the caller. + * @return Sizes of local parameters. + */ + [[nodiscard]] inline auto get_n_locals() const -> std::size_t { return NLocals; } + + /** + * @brief Getter for global parameter sizes. + * @param self The dynamic reference to the caller. + * @return Sizes of global parameters. + */ + [[nodiscard]] inline auto get_n_globals() const -> std::size_t { return NGlobals; } + /** * @brief Default spaceship comparison. */ 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/test/integration_tests/test_eigen.cpp b/test/integration_tests/test_eigen.cpp index 62f57a3d..f1334f54 100644 --- a/test/integration_tests/test_eigen.cpp +++ b/test/integration_tests/test_eigen.cpp @@ -1,11 +1,14 @@ #include "centipede/centipede.hpp" #include "centipede/data/entrypoint.hpp" +#include +#include +#include #include #include #include #include #include -#include +#include namespace { @@ -19,7 +22,8 @@ auto main() -> int auto mps = mps::MPS{ config }; mps.init(); - auto handler = centipede::Handler{}; + const auto n_globals = config.detector.spec.num_modules * 2; + auto handler = centipede::Handler{ n_globals }; auto entry_point_input = centipede::EntryPoint{}; for (auto _ : std::views::iota(0, 1000)) @@ -31,22 +35,33 @@ auto main() -> int entry_point_input.reset(); entry_point_input.set_measurement(entrypoint.measurement) .set_sigma(entrypoint.sigma) - .set_globals(entrypoint.globals) + .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(); if (not res) { - std::println("Error: {}", res.error()); - return EXIT_FAILURE; + 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: {}", is_ok.error()); + std::println("Error from the solving: {}", is_ok.error()); return EXIT_FAILURE; } 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_handler.cpp b/test/unit_tests/test_handler.cpp index 9b80afd7..89a1e065 100644 --- a/test/unit_tests/test_handler.cpp +++ b/test/unit_tests/test_handler.cpp @@ -12,26 +12,26 @@ 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 } }; + auto handler = HandlerType{ DEFAULT_MAX_GLOBAL_ID }; } // NOLINTBEGIN(readability-function-cognitive-complexity) TEST(handler, add_entrypoint) { - auto handler = Handler{}; + 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()); @@ -53,7 +53,7 @@ namespace centipede::test TEST(handler, local_derivs_incomp_numbers) { - auto handler = Handler{}; + 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 b8eade34..46c51cb8 100644 --- a/test/unit_tests/test_master_engine.cpp +++ b/test/unit_tests/test_master_engine.cpp @@ -41,7 +41,7 @@ namespace centipede::core::engine 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)); @@ -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); From 5c3005b267dc4dad3e1caa9a42d322a69a04d053 Mon Sep 17 00:00:00 2001 From: Yanzhao Wang Date: Thu, 18 Jun 2026 12:24:32 +0200 Subject: [PATCH 8/9] add more docs --- cmake/load_dependencies.cmake | 1 + doc/CMakeLists.txt | 3 + doc/latex_macros.tex | 2 + doc/main.md | 4 ++ doc/references.bib | 7 +++ doc/theories/linear_reg.md | 55 +++++++++++++++++++ .../centipede/core/engines/eigen_engine.hpp | 32 +++++++---- .../centipede/core/engines/master_engine.hpp | 2 +- source/centipede/core/handler.hpp | 5 ++ source/centipede/data/entrypoint.hpp | 14 ++--- source/centipede/data/entrypoint_base.hpp | 1 - test/integration_tests/CMakeLists.txt | 10 +++- test/integration_tests/test_eigen.cpp | 36 +++++++++++- util/mathjax-config.js | 8 +-- 14 files changed, 147 insertions(+), 33 deletions(-) create mode 100644 doc/latex_macros.tex create mode 100644 doc/references.bib create mode 100644 doc/theories/linear_reg.md diff --git a/cmake/load_dependencies.cmake b/cmake/load_dependencies.cmake index e4f68e79..0fe4bb5d 100644 --- a/cmake/load_dependencies.cmake +++ b/cmake/load_dependencies.cmake @@ -4,6 +4,7 @@ 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/doc/CMakeLists.txt b/doc/CMakeLists.txt index 78a0c3e6..5dd01700 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -31,6 +31,9 @@ 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") +set(DOXYGEN_FORMULA_MACROFILE "${PROJECT_SOURCE_DIR}/doc/latex_macros.tex") +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..8722fba4 --- /dev/null +++ b/doc/references.bib @@ -0,0 +1,7 @@ +@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}", + note = "[Online; accessed 17-June-2026]" + } diff --git a/doc/theories/linear_reg.md b/doc/theories/linear_reg.md new file mode 100644 index 00000000..3f8522cc --- /dev/null +++ b/doc/theories/linear_reg.md @@ -0,0 +1,55 @@ +# Linear regression {#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. + +## 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}} +\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) +\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} diff --git a/source/centipede/core/engines/eigen_engine.hpp b/source/centipede/core/engines/eigen_engine.hpp index 3bef76e8..4bd49d97 100644 --- a/source/centipede/core/engines/eigen_engine.hpp +++ b/source/centipede/core/engines/eigen_engine.hpp @@ -30,22 +30,30 @@ namespace centipede::core::engine { - namespace + /** + * @brief Engine template specialization for Eigen library implementation. + */ + template + 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); } }; - } // namespace - /** - * @brief Engine template specialization for Eigen library implementation. - */ - template - class Engine : public Base - { public: /** * @brief Matrix and vector used to solve global parameter updates @@ -107,10 +115,10 @@ 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) diff --git a/source/centipede/core/engines/master_engine.hpp b/source/centipede/core/engines/master_engine.hpp index 51157649..658aa623 100644 --- a/source/centipede/core/engines/master_engine.hpp +++ b/source/centipede/core/engines/master_engine.hpp @@ -147,7 +147,7 @@ namespace centipede::core::engine State current_state_; EngineImp engine_imp_{}; - //TODO: represent globals as mdspan, instead of relying on eigen. + // TODO: represent globals as mdspan, instead of relying on eigen. EngineImp::Globals globals_{}; void reset_state() diff --git a/source/centipede/core/handler.hpp b/source/centipede/core/handler.hpp index 483c9c55..09fa2b7c 100644 --- a/source/centipede/core/handler.hpp +++ b/source/centipede/core/handler.hpp @@ -69,6 +69,11 @@ namespace centipede [[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(); } diff --git a/source/centipede/data/entrypoint.hpp b/source/centipede/data/entrypoint.hpp index e20785be..c2ff6184 100644 --- a/source/centipede/data/entrypoint.hpp +++ b/source/centipede/data/entrypoint.hpp @@ -39,14 +39,12 @@ namespace centipede /** * @brief Getter for local parameter sizes. - * @param self The dynamic reference to the caller. * @return Sizes of local parameters. */ [[nodiscard]] inline auto get_n_locals() const -> std::size_t { return locals_.size(); } /** * @brief Getter for global parameter sizes. - * @param self The dynamic reference to the caller. * @return Sizes of global parameters. */ [[nodiscard]] inline auto get_n_globals() const -> std::size_t { return globals_.size(); } @@ -84,7 +82,6 @@ namespace centipede * values are emptied before being set by the new values. * * @tparam View Types of the view. - * @param self The dynamic reference to the caller. * @param view View object. * @return Universal reference to the caller. */ @@ -104,7 +101,6 @@ namespace centipede * new values. * * @tparam View Types of the view. - * @param self The dynamic reference to the caller. * @param view View object. * @return Universal reference to the caller. */ @@ -175,7 +171,9 @@ namespace centipede * * 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 @@ -200,14 +198,12 @@ namespace centipede /** * @brief Getter for local parameter sizes. - * @param self The dynamic reference to the caller. * @return Sizes of local parameters. */ [[nodiscard]] inline auto get_n_locals() const -> std::size_t { return NLocals; } /** * @brief Getter for global parameter sizes. - * @param self The dynamic reference to the caller. * @return Sizes of global parameters. */ [[nodiscard]] inline auto get_n_globals() const -> std::size_t { return NGlobals; } @@ -221,10 +217,9 @@ namespace centipede * @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 #NGlobals. + * input view is assumed to be larger or equal to @ref NGlobals_DY_Entrypoint "NGlobals". * * @tparam View Types of the view. - * @param self The dynamic reference to the caller. * @param view View object. * @return Universal reference to the caller. */ @@ -240,10 +235,9 @@ namespace centipede * @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 - * #NLocals. + * @ref NLocals_DY_Entrypoint "NLocals". * * @tparam View Types of the view. - * @param self The dynamic reference to the caller. * @param view View object. * @return Universal reference to the caller. */ diff --git a/source/centipede/data/entrypoint_base.hpp b/source/centipede/data/entrypoint_base.hpp index d41e2b2b..d8522d58 100644 --- a/source/centipede/data/entrypoint_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/test/integration_tests/CMakeLists.txt b/test/integration_tests/CMakeLists.txt index e934f0d2..bc860a0c 100644 --- a/test/integration_tests/CMakeLists.txt +++ b/test/integration_tests/CMakeLists.txt @@ -49,14 +49,20 @@ set_tests_properties( find_package(mps) -if(mps_FOUND) +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) + 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() diff --git a/test/integration_tests/test_eigen.cpp b/test/integration_tests/test_eigen.cpp index f1334f54..794600f4 100644 --- a/test/integration_tests/test_eigen.cpp +++ b/test/integration_tests/test_eigen.cpp @@ -1,20 +1,39 @@ #include "centipede/centipede.hpp" #include "centipede/data/entrypoint.hpp" #include -#include +#include +#include +#include +#include #include #include #include #include #include #include +#include #include namespace { - // void entrypoints_generator() {} + 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{}; @@ -25,6 +44,9 @@ auto main() -> int 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)) { @@ -34,6 +56,7 @@ auto main() -> int { 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) @@ -47,6 +70,9 @@ auto main() -> int } // 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) @@ -67,7 +93,13 @@ auto main() -> int 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/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"; From 5dab47bd94cd1e1fe729f5b9584c93e6be6b1129 Mon Sep 17 00:00:00 2001 From: Yanzhao Wang Date: Fri, 19 Jun 2026 11:31:06 +0200 Subject: [PATCH 9/9] add more docs --- doc/CMakeLists.txt | 3 +- doc/references.bib | 14 +++++++- doc/theories/linear_reg.md | 68 ++++++++++++++++++++++++++++++++++++-- 3 files changed, 80 insertions(+), 5 deletions(-) diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index 5dd01700..db8d50e3 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -31,8 +31,9 @@ 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") +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") diff --git a/doc/references.bib b/doc/references.bib index 8722fba4..5f03a7bf 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -3,5 +3,17 @@ @misc{ wiki:GausNewton 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}", - note = "[Online; accessed 17-June-2026]" } + +@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 index 3f8522cc..fc0c4645 100644 --- a/doc/theories/linear_reg.md +++ b/doc/theories/linear_reg.md @@ -1,55 +1,117 @@ -# Linear regression {#linear_reg} +# 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. -## Objective functions +## 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 +### 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.