diff --git a/src/specmicp_common/CMakeLists.txt b/src/specmicp_common/CMakeLists.txt index 309f9a1..ea4f390 100644 --- a/src/specmicp_common/CMakeLists.txt +++ b/src/specmicp_common/CMakeLists.txt @@ -1,203 +1,204 @@ # Main # ==== # var to store the files set(specmicp_common_srcs "" CACHE INTERNAL "specmicp_common files" FORCE) set(specmicp_common_headers "" CACHE INTERNAL "specmicp_common headers" FORCE ) # macro to add to vars macro(add_to_main_srcs_list LIST_NAME) set(tmp "") foreach (src ${${LIST_NAME}}) list(APPEND tmp ${CMAKE_CURRENT_SOURCE_DIR}/${src}) endforeach(src) set(specmicp_common_srcs "${specmicp_common_srcs};${tmp}" CACHE INTERNAL "specmicp common files" FORCE) endmacro(add_to_main_srcs_list) macro(add_to_main_headers_list LIST_NAME) set( tmp "") foreach(header ${${LIST_NAME}}) LIST(APPEND tmp ${CMAKE_CURRENT_SOURCE_DIR}/${header}) endforeach(header) set(specmicp_common_headers "${specmicp_common_headers};${tmp}" CACHE INTERNAL "headers for specmicp common" FORCE) endmacro(add_to_main_headers_list) # set( specmicp_common_main_srcs dateandtime.cpp filesystem.cpp log.cpp moving_average.cpp openmp_support.cpp string_algorithms.cpp timer.cpp ) add_to_main_srcs_list( specmicp_common_main_srcs ) set( specmicp_common_main_headers cached_vector.hpp compat.hpp dateandtime.hpp filesystem.hpp log_impl.hpp log.hpp macros.hpp moving_average.hpp openmp_support.hpp options_handler.hpp perfs_handler.hpp pimpl_ptr.hpp range_iterator.hpp range_iterator.inl return_code.hpp scope_guard.hpp string_algorithms.hpp timer.hpp types.hpp ) add_to_main_headers_list( specmicp_common_main_headers ) INSTALL(FILES ${specmicp_common_headers} DESTINATION ${INCLUDE_INSTALL_DIR}/specmicp_common ) add_subdirectory(cli) add_subdirectory(eigen) add_subdirectory(io) add_subdirectory(micpsolver) add_subdirectory(odeint) add_subdirectory(physics) add_subdirectory(plugins) add_subdirectory(sparse_solvers) +add_subdirectory(simplex) set_property(DIRECTORY APPEND PROPERTY CMAKE_CONFIGURE_DEPENDS config.h.in) add_custom_target(config_h SOURCES config.h.in) configure_file(config.h.in ${CMAKE_CURRENT_BINARY_DIR}/config.h ) INSTALL(FILES ${CMAKE_CURRENT_BINARY_DIR}/config.h DESTINATION ${INCLUDE_INSTALL_DIR}/specmicp_common) # build library # ============= spc_set_pgo_flag(${specmicp_common_srcs}) add_library(objspecmicp_common OBJECT ${specmicp_common_srcs} ${specmicp_common_headers} ) set_property(TARGET objspecmicp_common PROPERTY POSITION_INDEPENDENT_CODE 1) target_include_directories(objspecmicp_common PUBLIC ${EIGEN3_INCLUDE_DIR} ${HDF5_INCLUDE_DIRS}) target_include_directories(objspecmicp_common PUBLIC ${YAML_INCLUDE_DIRS}) target_include_directories(objspecmicp_common PUBLIC $) if(EIGEN3_UNSUPPORTED_FOUND) set_property(TARGET objspecmicp_common APPEND PROPERTY COMPILE_DEFINITIONS EIGEN3_UNSUPPORTED_FOUND) target_include_directories(objspecmicp_common PRIVATE ${EIGEN3_UNSUPPORTED_INCLUDE_DIR}) # added after endif(EIGEN3_UNSUPPORTED_FOUND) set_target_properties( objspecmicp_common PROPERTIES CXX_STANDARD ${SPECMICP_CXX_STANDARD} CXX_STANDARD_REQUIRED ON ) add_library(specmicp_common SHARED $) spc_set_objlib_optimization_flags(objspecmicp_common) spc_set_visibility_flag(objspecmicp_common) if (UNIX) LIST(APPEND specmicp_common_link_libraries dl) else() message(FATAL_ERROR "Plugin system only for POSIX at this time !") endif() install(TARGETS specmicp_common EXPORT ${SPECMICP_TARGET} LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR} ) # static libraries # ---------------- if(SPECMICP_BUILD_STATIC) add_library(specmicp_common_static STATIC $ ) install(TARGETS specmicp_common_static EXPORT ${SPECMICP_TARGET} ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR} ) else() add_library(specmicp_common_static EXCLUDE_FROM_ALL STATIC $ ) endif() set_target_properties(specmicp_common_static PROPERTIES OUTPUT_NAME specmicp_common) # link libraries and include directories # --------------------------------------- foreach(libtarget specmicp_common specmicp_common_static) target_link_libraries(${libtarget} PRIVATE ${specmicp_common_link_libraries}) target_link_libraries(${libtarget} PUBLIC PkgConfig::YAML) target_link_libraries(${libtarget} PUBLIC HDF5::hdf5) target_link_libraries(${libtarget} INTERFACE Eigen::eigen) spc_target_link_profile(${libtarget}) #if(SPECMICP_USE_BLAS) # target_link_libraries(${libtarget} PUBLIC ${EIGEN3_EXTRA_LINKFLAGS}) #endif() get_property(eigen_include_dirs TARGET Eigen::eigen PROPERTY INTERFACE_INCLUDE_DIRECTORIES) get_property(hdf5_include_dirs TARGET HDF5::hdf5 PROPERTY INTERFACE_INCLUDE_DIRECTORIES) get_property(yaml_include_dirs TARGET PkgConfig::YAML PROPERTY INTERFACE_INCLUDE_DIRECTORIES) target_include_directories(${libtarget} PUBLIC $) target_include_directories(${libtarget} PUBLIC $ $ ) target_include_directories(${libtarget} INTERFACE ${eigen_include_dirs} ) target_include_directories(${libtarget} INTERFACE ${hdf5_include_dirs} ) target_include_directories(${libtarget} INTERFACE ${yaml_include_dirs} ) # if(EIGEN3_UNSUPPORTED_FOUND) set_property(TARGET ${libtarget} APPEND PROPERTY COMPILE_DEFINITIONS EIGEN3_UNSUPPORTED_FOUND) set_property(TARGET ${libtarget} APPEND PROPERTY INTERFACE_COMPILE_DEFINITIONS EIGEN3_UNSUPPORTED_FOUND) target_include_directories(${libtarget} INTERFACE ${EIGEN3_UNSUPPORTED_INCLUDE_DIR}) # added after endif(EIGEN3_UNSUPPORTED_FOUND) spc_set_cxx_version_interface(${libtarget}) spc_set_target_optimization_flags(${libtarget}) endforeach(libtarget) diff --git a/src/specmicp_common/simplex/CMakeLists.txt b/src/specmicp_common/simplex/CMakeLists.txt new file mode 100644 index 0000000..e606345 --- /dev/null +++ b/src/specmicp_common/simplex/CMakeLists.txt @@ -0,0 +1,11 @@ +set(specmicp_common_simplex_headers + simplex_solver.hpp + simplex_solver.inl + simplex_program.hpp +) + +add_to_main_headers_list(specmicp_common_simplex_headers) + +install(FILES ${simplex_headers} + DESTINATION ${INCLUDE_INSTALL_DIR}/specmicp_common/simplex +) diff --git a/src/specmicp_common/simplex/simplex_program.hpp b/src/specmicp_common/simplex/simplex_program.hpp new file mode 100644 index 0000000..b95f156 --- /dev/null +++ b/src/specmicp_common/simplex/simplex_program.hpp @@ -0,0 +1,126 @@ +/* ============================================================================= + + Copyright (c) 2014-2017 F. Georget Princeton University + Copyright (c) 2017-2019 F. Georget EPFL + All rights reserved. + + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + 3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * + +============================================================================= */ + +//! \file simplex/simplex_program.hpp +//! \brief Base class for a simplex program + +#ifndef SPECMICP_COMMON_SIMPLEX_PROGRAM_HPP +#define SPECMICP_COMMON_SIMPLEX_PROGRAM_HPP + +#include "specmicp_common/types.hpp" + +namespace specmicp { +namespace simplex { + +//! \brief Base class for a simplex program +//! +//! Uses the CRTP. +//! +//! \f[ \mathrm{min} c^T \cdot x \f] +//! \f[ \mathrm{s.t.} A\cdot x = b +template +class SimplexProgram { +public: + using BasisType = Eigen::Matrix; // always the case for now + using MatrixType = Matrix; + + //! \brief Number of equalities constraints + index_t nb_equalities() const { + return derived()->nb_equalities_impl(); + } + //! \brief Number of variables + index_t nb_variables() const { + return derived()->nb_variables_impl(); + } + //! \brief Number of non zero variables + index_t size_B() const { + return nb_equalities(); + } + //! \brief Number of zero variables + index_t size_N() const { + return nb_variables()-size_B(); + } + //! \brief Fill the starting guess + //! + //! As B assumed to be the identity value, set the value to b + void get_starting_guess(Eigen::Ref m_xB) { + derived()->get_starting_guess_impl(m_xB); + } + //! \brief Coefficients of the objective functions + void get_cB(Eigen::Ref cB) { + derived()->get_cB_impl(cB); + } + //! \brief compute the \f$s_N\f$ vector + //! + //! \f[ c_N - N^T \cdot \lambda \f] + void compute_sN( + const Eigen::Ref& Nbasis, + const Eigen::Ref& lambda, + Eigen::Ref sN + ) + { + derived()->compute_sN_impl(Nbasis, lambda, sN); + } + + //! \brief Return \f$A(:,q)\f$ + void fill_Aq( + const Eigen::Ref& Nbasis, + Eigen::Ref Aq, + index_t q + ) + { + derived()->fill_Aq_impl(Nbasis, Aq, q); + } + //! \brief Return the coefficient of the objective function corresponding to q + scalar_t cN_q( + const Eigen::Ref& Nbasis, + index_t q) { + return derived()->cN_q_impl(Nbasis, q); + } +private: + //! \brief Return the derived class + Derived* derived() { + return static_cast(this); + } + //! \brief Return a const pointer to the derived class + const Derived* const derived() const { + return static_cast(this); + } +}; + +} +} + +#endif // SPECMICP_COMMON_SIMPLEX_PROGRAM_HPP + diff --git a/src/specmicp_common/simplex/simplex_solver.hpp b/src/specmicp_common/simplex/simplex_solver.hpp new file mode 100644 index 0000000..880b06f --- /dev/null +++ b/src/specmicp_common/simplex/simplex_solver.hpp @@ -0,0 +1,135 @@ +/* ============================================================================= + + Copyright (c) 2014-2017 F. Georget Princeton University + Copyright (c) 2017-2019 F. Georget EPFL + All rights reserved. + + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + 3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * + +============================================================================= */ + +#ifndef SPECMICP_COMMON_SIMPLEX_SOLVER_HPP +#define SPECMICP_COMMON_SIMPLEX_SOLVER_HPP + +//! \file simplex/simplex_solver.hpp +//! \brief Declaration of the simplex solver + +#include "specmicp_common/types.hpp" +#include + + +namespace specmicp { +//! \namespace specmicp::simplex +//! \brief A simplex solver +namespace simplex { + +//! \brief The return code of the simplex algorithm +enum class SimplexSolverReturnCode +{ + Error = -10, //!< Generic error + Unbounded = -1, //!< Problem is unbounded + NotConverged = 0, //!< Solver has not converged yet + Success = 1 //!< Solver has converged +}; + +//! \brief The Simplex solver +//! +//! Using the Revised Simplex algorithm +//! Ref: Nocedal & Wright (2006) +template +class SimplexSolver { + +using MatrixType = Matrix; // has to be for now +using BasisType = Eigen::Matrix;// + +public: + SimplexSolver(SimplexProgram& program): + m_program(program) + {} + + //! \brief Solve the linear system + //! + //! Return SimplexSolverReturnCode::Success if converged + SimplexSolverReturnCode solve(); + //! \brief Solve one step of the simplex method + //! + //! Return SimplexSolverReturnCode::Success if converged + SimplexSolverReturnCode one_step(); + //! \brief Return the solution + //! + //! \param[out] val vector of the variables + //! \return the value of the objective function + scalar_t get_solution(Vector& val); + +private: + //! \brief Initialize the data structures and fill values + void initialize(); + //! \brief The pricing step, how to improve the minimization + void pricing_step(); + //! \brief Check convergence + SimplexSolverReturnCode check_objective(); + //! \brief Select the variable to enter the basis + index_t select_q(); + //! \brief Needed to select which varible to exit + void solve_for_d(index_t q); + //! \brief Check if the problem is unbounded + SimplexSolverReturnCode check_unboundness(); + //! \brief Select the variable to exit the basis + std::pair select_p(); + //! \brief Update the variables vector + void update_xB(scalar_t xq); + //! \brief Modify the basis + void modify_basis(index_t p, index_t q, scalar_t xq); + //! \brief Start the simplex algorithm + void compute_start_step(); + + SimplexProgram& m_program; //!< Reference to the program + BasisType m_B; //!< Basis of non-zero variables + BasisType m_N; //!< Basis of zero variables + + Matrix m_mB; //!< The non-zero part of the equalities + + Vector m_xB; //!< The non-zero variables + + Vector m_sN; //!< Gradient of the objective function + Vector m_cB; //!< Coefficient of non-zero term the objective functions + + Vector m_db; //!< Column of the equalities to switch in + Vector m_d; //!< Vector used to find the varaible to set to zero + + Eigen::PartialPivLU m_LU; //!< LU decomposition + +}; + +} +} + +#endif // SPECMICP_COMMON_SIMPLEX_SOLVER_HPP + +#ifndef SPECMICP_COMMON_SIMPLEX_SOLVER_INL +#include "simplex_solver.inl" +#endif diff --git a/src/specmicp_common/simplex/simplex_solver.inl b/src/specmicp_common/simplex/simplex_solver.inl new file mode 100644 index 0000000..60e037b --- /dev/null +++ b/src/specmicp_common/simplex/simplex_solver.inl @@ -0,0 +1,208 @@ +/* ============================================================================= + + Copyright (c) 2014-2017 F. Georget Princeton University + Copyright (c) 2017-2019 F. Georget EPFL + All rights reserved. + + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + + 3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * + +============================================================================= */ + +#ifndef SPECMICP_COMMON_SIMPLEX_SOLVER_INL +#define SPECMICP_COMMON_SIMPLEX_SOLVER_INL + +// for syntax coloring +#ifndef SPECMICP_COMMON_SIMPLEX_SOLVER_HPP +#include "simplex_solver.hpp" +#endif //SPECMICP_COMMON_SIMPLEX_SOLVER_HPP + +#include "specmicp_common/eigen/vector_checker.hpp" + +#include + +namespace specmicp { +namespace simplex { + +template +void SimplexSolver::initialize() +{ + auto m = m_program.size_B(); + auto p = m_program.size_N(); + + m_B.setZero(m); + for (auto i=0; i +SimplexSolverReturnCode SimplexSolver::solve() +{ + initialize(); + SimplexSolverReturnCode ret = SimplexSolverReturnCode::NotConverged; + index_t cnt = 0; + while (ret == SimplexSolverReturnCode::NotConverged and cnt<1000) { + ret = one_step(); + ++cnt; + } + return ret; +} + +template +SimplexSolverReturnCode SimplexSolver::one_step() +{ + compute_start_step(); + pricing_step(); + SimplexSolverReturnCode ret = check_objective(); + if (ret != SimplexSolverReturnCode::NotConverged) { + return ret; + } + index_t q = select_q(); + solve_for_d(q); + ret = check_unboundness(); + if (ret != SimplexSolverReturnCode::NotConverged) { + return ret; + } + std::pair pandxq = select_p(); + update_xB(pandxq.second); + modify_basis(pandxq.first, q, pandxq.second); + return SimplexSolverReturnCode::NotConverged; +} + + +template +void SimplexSolver::compute_start_step() +{ + m_LU.compute(m_mB); +} + +template +void SimplexSolver::pricing_step() +{ + Vector lamb = m_LU.transpose().solve(m_cB); + m_program.compute_sN(m_N, lamb, m_sN); + +} + +template +SimplexSolverReturnCode SimplexSolver::check_objective() +{ + if (vector_checker::all(m_sN) >= 0.0) { + return SimplexSolverReturnCode::Success; + } + return SimplexSolverReturnCode::NotConverged; +} + + + +template +index_t SimplexSolver::select_q() +{ + index_t q; + m_sN.minCoeff(&q); + return q; + +} + +template +void SimplexSolver::solve_for_d(index_t q) +{ + m_program.fill_Aq(m_N, m_db, q); + + m_d = m_LU.solve(m_db); +} + +template +SimplexSolverReturnCode SimplexSolver::check_unboundness() +{ + if (vector_checker::all(m_d) <= 0.0) { + return SimplexSolverReturnCode::Unbounded; + } + return SimplexSolverReturnCode::NotConverged; +} + +template +std::pair SimplexSolver::select_p() +{ + index_t p; + scalar_t x_q = (m_xB.array()/m_d.array()).minCoeff(&p); + return std::pair(p, x_q); +} + +template +void SimplexSolver::update_xB(scalar_t xq) +{ + m_xB -= m_d*xq; +} + +template +void SimplexSolver::modify_basis(index_t p, index_t q, scalar_t xq) +{ + // B+ + m_mB.col(p) = m_db; + //update cB + m_cB(p) = m_program.cN_q(m_N, q); + + auto tmp = m_B(p); + m_B(p)=m_N(q); + m_N(q) = tmp; + + m_xB(p) = xq; +} + +template +scalar_t SimplexSolver::get_solution(Vector& x) +{ + x = Vector::Zero(m_program.nb_variables()); + for (index_t i=0; i + +using namespace specmicp; +using namespace specmicp::simplex; + +class TestSimplexProgram: + public SimplexProgram +{ +friend class SimplexProgram; + +public: + TestSimplexProgram( + Vector& atotal_concentration, + Vector& ac_component, + Vector& ac_secondary, + Matrix& anu + ): + m_nb_equalities(atotal_concentration.rows()), + m_nb_variables(ac_component.rows()+ac_secondary.rows()), + total_concentration(atotal_concentration), + c_component(ac_component), + c_secondary(ac_secondary), + nu(anu) + {} +protected: + index_t nb_equalities_impl() const{ + return m_nb_equalities; + } + index_t nb_variables_impl() const { + return m_nb_variables; + } + + void get_starting_guess_impl(Eigen::Ref m_xB) { + m_xB = total_concentration; + } + + void get_cB_impl(Eigen::Ref cB) { + cB = c_component; + } + + void compute_sN_impl( + const Eigen::Ref& Nbasis, + const Eigen::Ref& lambda, + Eigen::Ref sN + ) { + for (index_t k=0; k& Nbasis, + Eigen::Ref Aq, + index_t q + ) + { + auto r = Nbasis(q); + if (is_component(r)) { + Aq.setZero(); + Aq(r) = 1.0; + } else { + const auto j = index_secondary(Nbasis(q)); + for (index_t i=0; i& Nbasis, + index_t q) { + auto k = Nbasis(q); + if (is_component(k)) { + return c_component(k); + } else { + auto j = index_secondary(k); + return c_secondary(j); + } + } + + bool is_component(index_t i) { + return i < nb_equalities_impl(); + } + bool is_secondary(index_t i) { + return (not is_component(i)); + } + index_t index_secondary(index_t i) { + return i-nb_equalities_impl(); + } + +private: + index_t m_nb_equalities; + index_t m_nb_variables; + + Vector total_concentration; + Vector c_component; + Vector c_secondary; + + Matrix nu; +}; + + +TEST_CASE("Simplex solver 3x4", "[simplex]") { + + // A B C + // AB AC A2BC A3B2C + + Vector total_concentration(3); + Vector c_component(3); + Vector c_secondary(4); + Matrix nu(4,3); + + total_concentration << 3,2,2; + c_component << -2,-5,-2; + c_secondary << -10, -6, -20, -10; + + nu << 1, 1, 0, + 1, 0, 1, + 2, 1, 1, + 3, 2, 1; + + + SECTION("Program test") { + auto prog = TestSimplexProgram(total_concentration, c_component, c_secondary, nu); + CHECK(prog.nb_equalities() == 3); + CHECK(prog.nb_variables() == 7); + + Vector start_guess(3); + prog.get_starting_guess(start_guess); + CHECK((total_concentration - start_guess).norm() <= 1e-8); + + Vector cB(3); + prog.get_cB(cB); + CHECK((c_component-cB).norm() <= 1e-8); + + } + + SECTION("Solver test") { + auto prog = TestSimplexProgram(total_concentration, c_component, c_secondary, nu); + auto solver = SimplexSolver(prog); + + auto ret = solver.solve(); + CHECK(ret == SimplexSolverReturnCode::Success); + + Vector x; + auto c = solver.get_solution(x); + std::cout << "Objective : " << c << std::endl; + + std::cout << "Vars : \n" << x << std::endl; + } + + +}