diff --git a/.github/workflows/macos-unit.yml b/.github/workflows/macos-unit.yml index c493f89..ed26c85 100644 --- a/.github/workflows/macos-unit.yml +++ b/.github/workflows/macos-unit.yml @@ -3,8 +3,8 @@ name: macOS unit on: push: branches: - - master - #- develop + - main + - develop pull_request: branches: - '**' diff --git a/.github/workflows/ubuntu-unit.yml b/.github/workflows/ubuntu-unit.yml index 3e1543e..654ad3b 100644 --- a/.github/workflows/ubuntu-unit.yml +++ b/.github/workflows/ubuntu-unit.yml @@ -3,8 +3,8 @@ name: Ubuntu unit on: push: branches: - - master - #- develop + - main + - develop pull_request: branches: - '**' @@ -30,24 +30,24 @@ jobs: - CC: gcc-13 CXX: g++-13 compiler: gcc-13 g++-13 - - CC: clang-12 - CXX: clang++-12 - compiler: clang-12 - - CC: clang-13 - CXX: clang++-13 - compiler: clang-13 - - CC: clang-14 - CXX: clang++-14 - compiler: clang-14 + # - CC: clang-12 -> error: call to consteval function 'std::chrono::hh_mm_ss::_S_fractional_width' is not a constant expression + # CXX: clang++-12 + # compiler: clang-12 + # - CC: clang-13 + # CXX: clang++-13 + # compiler: clang-13 + # - CC: clang-14 + # CXX: clang++-14 + # compiler: clang-14 - CC: clang-15 CXX: clang++-15 compiler: clang-15 - - CC: clang-16 - CXX: clang++-16 - compiler: clang-16 - - CC: clang-17 - CXX: clang++-17 - compiler: clang-17 + # - CC: clang-16 -> we don't have these as package. + # CXX: clang++-16 + # compiler: clang-16 + # - CC: clang-17 + # CXX: clang++-17 + # compiler: clang-17 steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/windows-unit.yml b/.github/workflows/windows-unit.yml index f9545c1..6dc4c06 100644 --- a/.github/workflows/windows-unit.yml +++ b/.github/workflows/windows-unit.yml @@ -3,8 +3,8 @@ name: Windows unit on: push: branches: - - master - #- develop + - main + - develop pull_request: branches: - '**' diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e4b232..1b25410 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,6 @@ This changelog contains a non-exhaustive list of new features and notable bug-fi ## New features * OSQP support is removed. -* Simplex method is added. * HiGHS solver is added. diff --git a/CMakeLists.txt b/CMakeLists.txt index 6caddb5..befd911 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,11 +26,11 @@ project(DTWC++ VERSION "0.0.4" set(DTWC_ENABLE_GUROBI ON) +include(cmake/Dependencies.cmake) # Include external projects include(cmake/FindGUROBI.cmake) # include(cmake/ProjectOptions.cmake) include(cmake/StandardProjectSettings.cmake) include(cmake/PreventInSourceBuilds.cmake) -include(cmake/Dependencies.cmake) # Include external projects dtwc_setup_dependencies() @@ -43,7 +43,7 @@ target_compile_features(project_options INTERFACE cxx_std_20) enable_testing() -add_subdirectory(dtwc) +add_subdirectory(dtwc bin) diff --git a/cmake/Coverage.cmake b/cmake/Coverage.cmake index 29b43ad..90710b8 100644 --- a/cmake/Coverage.cmake +++ b/cmake/Coverage.cmake @@ -1,9 +1,9 @@ option (DTWC_ENABLE_COVERAGE "Enable coverage reporting for GCC or Clang" OFF) # Setup macro for coverage testing for GCC or Clang -macro(add_executable_with_coverage_and_test TARGET_NAME CPP_FILE_NAME) - add_executable(${TARGET_NAME} ${CPP_FILE_NAME}) +macro(add_executable_with_coverage_and_test TARGET_NAME) + add_executable(${TARGET_NAME} "${TARGET_NAME}.cpp") target_link_libraries(${TARGET_NAME} PRIVATE dtwc++ Catch2::Catch2WithMain) - add_test(NAME ${TARGET_NAME} COMMAND ${TARGET_NAME}) + add_test(NAME ${TARGET_NAME} COMMAND ${TARGET_NAME} WORKING_DIRECTORY bin) if (DTWC_ENABLE_COVERAGE) if (${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU" OR ${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang") message (STATUS "Configuring with coverage") diff --git a/cmake/Dependencies.cmake b/cmake/Dependencies.cmake index bda02cf..91a5a0a 100644 --- a/cmake/Dependencies.cmake +++ b/cmake/Dependencies.cmake @@ -1,5 +1,5 @@ # This cmake file is to add external dependency projects. -# Adapted from https://github.com/cpp-best-practices/cmake_template/tree/main +# Adapted from https://github.com/cpp-best-practices/cmake_tecomplate/tree/main include(cmake/CPM.cmake) # Done as a function so that updates to variables like @@ -10,35 +10,37 @@ function(dtwc_setup_dependencies) # For each dependency, see if it's # already been provided to us by a parent project - if(NOT TARGET range-v3) # Range-v3 library: - CPMAddPackage( - NAME range-v3 - URL "https://github.com/ericniebler/range-v3/archive/refs/tags/0.12.0.tar.gz" - DOWNLOAD_ONLY YES - ) + # if(NOT TARGET range-v3) # Range-v3 library: + # CPMAddPackage( + # NAME range-v3 + # URL "https://github.com/ericniebler/range-v3/archive/refs/tags/0.12.0.tar.gz" + # DOWNLOAD_ONLY YES + # ) - if(range-v3_ADDED) - add_library(range-v3 INTERFACE IMPORTED) - target_include_directories(range-v3 INTERFACE ${range-v3_SOURCE_DIR}/include) - endif() - endif() + # if(range-v3_ADDED) + # add_library(range-v3 INTERFACE IMPORTED) + # target_include_directories(range-v3 INTERFACE ${range-v3_SOURCE_DIR}/include) + # endif() + # endif() - if(NOT TARGET eigen) # Eigen library: - CPMAddPackage( - NAME eigen - URL "https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz" - DOWNLOAD_ONLY YES - ) + # if(NOT TARGET eigen) # Eigen library: + # CPMAddPackage( + # NAME eigen + # URL "https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz" + # DOWNLOAD_ONLY YES + # ) - add_library(eigen INTERFACE IMPORTED) - target_include_directories(eigen SYSTEM INTERFACE ${eigen_SOURCE_DIR}) - endif() + # add_library(eigen INTERFACE IMPORTED) + # target_include_directories(eigen SYSTEM INTERFACE ${eigen_SOURCE_DIR}) + # endif() if(NOT TARGET Catch2::Catch2WithMain) # Catch2 library: CPMAddPackage( NAME Catch2 - URL "https://github.com/catchorg/Catch2/archive/refs/tags/v3.4.0.tar.gz" + URL "https://github.com/catchorg/Catch2/archive/refs/tags/v3.3.2.tar.gz" + OPTIONS + "CATCH_INSTALL_DOCS OFF" "CATCH_INSTALL_EXTRAS OFF" "CATCH_BUILD_TESTING OFF" ) endif() @@ -50,17 +52,15 @@ function(dtwc_setup_dependencies) ) endif() + # HiGHS library: if(NOT TARGET highs::highs)# HiGHS library: - include(FetchContent) - - FetchContent_Declare( - highs - GIT_REPOSITORY "https://github.com/ERGO-Code/HiGHS.git" - GIT_TAG "bazel" - ) - set(FAST_BUILD ON CACHE INTERNAL "Fast Build") - - FetchContent_MakeAvailable(highs) - + CPMAddPackage( + NAME highs + URL "https://github.com/ERGO-Code/HiGHS/archive/refs/tags/v1.6.0.tar.gz" + SYSTEM + EXCLUDE_FROM_ALL + OPTIONS + "FAST_BUILD ON" "BUILD_TESTING OFF" "BUILD_EXAMPLES OFF") endif() + endfunction() diff --git a/cmake/StandardProjectSettings.cmake b/cmake/StandardProjectSettings.cmake index fe1682a..2bc79b4 100644 --- a/cmake/StandardProjectSettings.cmake +++ b/cmake/StandardProjectSettings.cmake @@ -48,7 +48,7 @@ run_vcvarsall() message(STATUS "Host system: ${CMAKE_HOST_SYSTEM}") message(STATUS "Build type: ${CMAKE_BUILD_TYPE}") -# set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_DEBUG ${CMAKE_CURRENT_SOURCE_DIR}/bin/Debug) -# set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_CURRENT_SOURCE_DIR}/bin/Release) -# set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELWITHDEBINFO ${CMAKE_CURRENT_SOURCE_DIR}/bin/RelWithDebInfo) -# set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_MINSIZEREL ${CMAKE_CURRENT_SOURCE_DIR}/bin/MinSizeRel) \ No newline at end of file +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_DEBUG ${CMAKE_CURRENT_SOURCE_DIR}/bin/Debug) +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELEASE ${CMAKE_CURRENT_SOURCE_DIR}/bin/Release) +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_RELWITHDEBINFO ${CMAKE_CURRENT_SOURCE_DIR}/bin/RelWithDebInfo) +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY_MINSIZEREL ${CMAKE_CURRENT_SOURCE_DIR}/bin/MinSizeRel) \ No newline at end of file diff --git a/develop/TODO.md b/develop/TODO.md index 04bc822..e8b4178 100644 --- a/develop/TODO.md +++ b/develop/TODO.md @@ -29,8 +29,8 @@ - [ ] dtwclust in R - [ ] Encapsulating Data and related functions in one folder. - [ ] Open-source solver addition. - - [ ] Simplex is added. - - [ ] HiGHS is added. + - [ ] Simplex is added. -> temporarily removed. + - [x] HiGHS is added. - [x] Exploration of totally unimodular matrices. -> Not totally unimodular. - [ ] Creating DTW objects taking distance/band as a policy-based design - [x] Doxygen website? diff --git a/dtwc/CMakeLists.txt b/dtwc/CMakeLists.txt index ee70edb..95b114a 100644 --- a/dtwc/CMakeLists.txt +++ b/dtwc/CMakeLists.txt @@ -1,4 +1,3 @@ -add_subdirectory(solver) add_library(dtwc++ STATIC) @@ -16,9 +15,8 @@ target_link_libraries(dtwc++ PRIVATE project_warnings project_options - solvers - eigen - range-v3 + # eigen + # range-v3 fmt highs::highs ) @@ -33,7 +31,7 @@ if(DTWC_ENABLE_GUROBI) endif() -target_include_directories(dtwc++ PRIVATE ${HIGHS_INCLUDE_DIRS} ${range-v3_INCLUDE_DIRS}) +target_include_directories(dtwc++ PRIVATE ${range-v3_INCLUDE_DIRS} ) target_compile_definitions(dtwc++ PRIVATE ROOT_FOLDER="${PROJECT_SOURCE_DIR}") add_executable(dtwc_main diff --git a/dtwc/Problem.cpp b/dtwc/Problem.cpp index 251b0df..3571a39 100644 --- a/dtwc/Problem.cpp +++ b/dtwc/Problem.cpp @@ -159,12 +159,6 @@ void Problem::init_Kmeanspp() void Problem::cluster_by_MIP() { MIP_clustering_byHiGHS(*this); - // MIP_clustering_byDenseSimplex(*this); - // MIP_clustering_bySparseSimplex(*this); - - // if (settings::is_relaxed) - // MIP_clustering_byGurobi_relaxed(*this); - // else // MIP_clustering_byGurobi(*this); } diff --git a/dtwc/fileOperations.hpp b/dtwc/fileOperations.hpp index 2f7b1de..3380f23 100644 --- a/dtwc/fileOperations.hpp +++ b/dtwc/fileOperations.hpp @@ -21,6 +21,8 @@ #include // for pair #include // for vector #include +#include +#include namespace dtwc { diff --git a/dtwc/main.cpp b/dtwc/main.cpp index f4b3c8f..79d6715 100644 --- a/dtwc/main.cpp +++ b/dtwc/main.cpp @@ -4,21 +4,8 @@ #include - -#include "solver/EqualityConstraints.hpp" -#include "solver/test.hpp" - int main() { - // auto [eq, c] = dtwc::solver::get_prob_small(); - - // dtwc::solver::SparseSimplex prob_small(eq, c); - - // prob_small.gomoryAlgorithm(); - // auto [solution_small, copt_small] = prob_small.getResults(); - - // fmt::println("Solution: {} and Copt = [{}]\n", solution_small, copt_small); - // Here are some examples. You can either take the contents of example functions into main or modify and run them. dtwc::Clock clk; // Create a clock object @@ -78,7 +65,7 @@ int main() // std::cout << "Finished all tasks " << clk << "\n"; - dtwc::benchmarks::run_all(); + // dtwc::benchmarks::run_all(); std::cout << "Finished all tasks " << clk << "\n"; } diff --git a/dtwc/mip.cpp b/dtwc/mip.cpp index df8e77a..1b76c49 100644 --- a/dtwc/mip.cpp +++ b/dtwc/mip.cpp @@ -8,13 +8,11 @@ */ #include "Data.hpp" // for Data -#include "sparse_util.hpp" // for Triplet, RowMajor +#include "types/types.hpp" // for Triplet, RowMajor #include "mip.hpp" #include "Problem.hpp" #include "settings.hpp" #include "timing.hpp" -#include "solver/Simplex.hpp" -#include "solver/SparseSimplex.hpp" #include #include @@ -54,37 +52,6 @@ void extract_solution(Problem &prob, auto &solution) } } -template -void MIP_clustering_bySimplex(Problem &prob) -{ - std::cout << "Simplex is being called!" << std::endl; - dtwc::Clock clk; // Create a clock object - - prob.clear_clusters(); - - thread_local auto simplexSolver = T(prob); - - std::cout << "Problem formulation finished in " << clk << '\n'; - simplexSolver.gomoryAlgorithm(); - std::cout << "Problem solution finished in " << clk << '\n'; - - auto [solution, copt] = simplexSolver.getResults(); - - fmt::println("Solution: {} and Copt = [{}]\n", solution, copt); - - extract_solution(prob, solution); -} - -void MIP_clustering_bySparseSimplex(Problem &prob) -{ - MIP_clustering_bySimplex(prob); -} - -void MIP_clustering_byDenseSimplex(Problem &prob) -{ - MIP_clustering_bySimplex(prob); -} - void MIP_clustering_byHiGHS(Problem &prob) { std::cout << "HiGS is being called!" << std::endl; diff --git a/dtwc/mip.hpp b/dtwc/mip.hpp index 04518d4..be06819 100644 --- a/dtwc/mip.hpp +++ b/dtwc/mip.hpp @@ -13,10 +13,6 @@ namespace dtwc { class Problem; void MIP_clustering_byGurobi(Problem &prob); -void MIP_clustering_byGurobi_relaxed(Problem &prob); - -void MIP_clustering_bySparseSimplex(Problem &prob); -void MIP_clustering_byDenseSimplex(Problem &prob); void MIP_clustering_byHiGHS(Problem &prob); diff --git a/dtwc/mip_Gurobi.cpp b/dtwc/mip_Gurobi.cpp index 3c6506c..b035c55 100644 --- a/dtwc/mip_Gurobi.cpp +++ b/dtwc/mip_Gurobi.cpp @@ -108,107 +108,4 @@ void MIP_clustering_byGurobi(Problem &prob) #endif } - -void MIP_clustering_byGurobi_relaxed(Problem &prob) -{ -#ifdef DTWC_ENABLE_GUROBI - std::cout << "Gurobi-relaxed has been called!" << std::endl; - - const auto Nb = prob.data.size(); - const auto Nc = prob.cluster_size(); - - prob.clear_clusters(); - - try { - GRBEnv env = GRBEnv(); - GRBModel model = GRBModel(env); - - // Create variables - std::unique_ptr isCluster{ model.addVars(Nb, GRB_CONTINUOUS) }; - std::unique_ptr w{ model.addVars(Nb * Nb, GRB_CONTINUOUS) }; - - for (int i{ 0 }; i < Nb; i++) { - GRBLinExpr lhs = 0; - for (int j{ 0 }; j < Nb; j++) { - lhs += w[j + i * Nb]; - model.addConstr(w[j + i * Nb] <= 1); // For relaxed version. - model.addConstr(w[j + i * Nb] >= 0); // For relaxed version. - } - model.addConstr(lhs, '=', 1.0); - } - - - for (int j{ 0 }; j < Nb; j++) - for (int i{ 0 }; i < Nb; i++) - model.addConstr(w[i + j * Nb] <= isCluster[i]); - - { - GRBLinExpr lhs = 0; - for (int i{ 0 }; i < Nb; i++) { - lhs += isCluster[i]; - model.addConstr(isCluster[i] <= 1); // For relaxed version. - model.addConstr(isCluster[i] >= 0); // For relaxed version. - } - - model.addConstr(lhs == Nc); // There should be Nc clusters. - } - - // Set objective - GRBLinExpr obj = 0; - for (int j{ 0 }; j < Nb; j++) - for (int i{ 0 }; i < Nb; i++) - obj += w[i + j * Nb] * prob.distByInd_scaled(i, j); // Some scaling! - - model.setObjective(obj, GRB_MINIMIZE); - std::cout << "Finished setting up the MILP problem." << std::endl; - - model.set(GRB_IntParam_Threads, 3); // Set to dual simplex? - - // model.set(GRB_IntParam_Method, 1); // Set to dual simplex? - // model.set(GRB_IntParam_NumericFocus, 3); // Much numerics - - // model.set(GRB_IntParam_Presolve, 2); - - model.optimize(); - - std::vector test; - for (int i{ 0 }; i < Nb * Nb; i++) - test.push_back(w[i].get(GRB_DoubleAttr_X)); - - for (int i{ 0 }; i < Nb; i++) - if (isCluster[i].get(GRB_DoubleAttr_X) > 0.9) - prob.centroids_ind.push_back(i); - else if (isCluster[i].get(GRB_DoubleAttr_X) > 0.1) { - std::cerr << "Cluster " << i << " has value of " << isCluster[i].get(GRB_DoubleAttr_X) << " which should not happen for turtley unimodular matrices!\n"; - throw 10000; // #TODO more meaningful error codes? - } - - prob.clusters_ind.resize(Nb); - - int i_cluster = 0; - for (auto i : prob.centroids_ind) { - prob.cluster_members.emplace_back(); - for (int j{ 0 }; j < Nb; j++) - if (w[i + j * Nb].get(GRB_DoubleAttr_X) > 0.9) { - prob.clusters_ind[j] = i_cluster; - prob.cluster_members.back().push_back(j); - } else if (w[i + j * Nb].get(GRB_DoubleAttr_X) > 0.1) { - std::cerr << "Weight " << i + j * Nb << " has value of " << w[i + j * Nb].get(GRB_DoubleAttr_X) << " which should not happen for turtley unimodular matrices!\n"; - throw 10000; // #TODO more meaningful error codes? - } - - i_cluster++; - } - - } catch (GRBException &e) { - std::cout << "Error code = " << e.getErrorCode() << std::endl - << e.getMessage() << std::endl; - } catch (...) { - std::cout << "Unknown Exception during Gurobi relaxed, try regular Gurobi optimisation" << std::endl; - MIP_clustering_byGurobi(prob); - } - -#endif -} - } // namespace dtwc \ No newline at end of file diff --git a/dtwc/solver/CMakeLists.txt b/dtwc/solver/CMakeLists.txt deleted file mode 100644 index 9ed702f..0000000 --- a/dtwc/solver/CMakeLists.txt +++ /dev/null @@ -1,27 +0,0 @@ -cmake_minimum_required(VERSION 3.21) -add_library(solvers OBJECT) - -target_sources(solvers - PRIVATE - Simplex.cpp - SparseSimplex.cpp - SimplexTable.cpp - SimplexRowTable.cpp - SimplexFlatRowTable.cpp - SimplexFlatTable.cpp - - PUBLIC - solver_util.hpp - Simplex.hpp - EqualityConstraints.hpp - SparseSimplex.hpp - EqualityConstraints.hpp - SimplexTable.hpp - SimplexRowTable.hpp - SimplexFlatRowTable.hpp - SimplexFlatTable.hpp - - ) - -target_include_directories(solvers PUBLIC .) -target_link_libraries(solvers PRIVATE eigen fmt range-v3 highs::highs) \ No newline at end of file diff --git a/dtwc/solver/EqualityConstraints.hpp b/dtwc/solver/EqualityConstraints.hpp deleted file mode 100644 index 866225d..0000000 --- a/dtwc/solver/EqualityConstraints.hpp +++ /dev/null @@ -1,85 +0,0 @@ -/* - * SimplexSolver.hpp - * - * Helper class for sparse equality constraint. - - * Created on: 22 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#pragma once - -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "SparseMatrix.hpp" - -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -namespace dtwc::solver { - -struct EqualityConstraints -{ - SparseMatrix A; - std::vector b; - - EqualityConstraints() = default; - EqualityConstraints(int m_, int n_) : A(m_, n_), b(m_, 0.0) {} -}; - -EqualityConstraints inline defaultConstraints(int Nb, int Nc) -{ - const auto Neq = Nb + 1; - const auto Nineq = Nb * (Nb - 1); - const auto Nconstraints = Neq + Nineq; - - const auto Nvar_original = Nb * Nb; - const auto N_slack = Nineq; - const auto Nvar = Nvar_original + N_slack; // x1--xN^2 + s_slack - - auto eq = EqualityConstraints(Nconstraints, Nvar); - // Create b matrix: - eq.b[0] = Nc; - for (int i = 0; i < Nb; ++i) - eq.b[i + 1] = 1; - - // Create A matrix: - for (int i = 0; i < Nineq; i++) - eq.A(Neq + i, Nvar_original + i) = 1; - - for (int i = 0; i < Nb; ++i) { - eq.A(0, i * (Nb + 1)) = 1.0; // Sum of diagonals is Nc - - for (int j = 0; j < Nb; j++) - eq.A(1 + j, Nb * i + j) = 1; // Every element belongs to one cluster. - - // --------------- - int shift = 0; - for (int j = 0; j < Nb; j++) { - const int block_begin_row = Nb + 1 + (Nb - 1) * i; - const int block_begin_col = Nb * i; - if (i == j) { - for (int k = 0; k < (Nb - 1); k++) - eq.A(block_begin_row + k, block_begin_col + j) = -1; - shift = 1; - } else - eq.A(block_begin_row + j - shift, block_begin_col + j) = 1; - } - } - - return eq; -}; - -} // namespace dtwc::solver \ No newline at end of file diff --git a/dtwc/solver/Simplex.cpp b/dtwc/solver/Simplex.cpp deleted file mode 100644 index 63a4b02..0000000 --- a/dtwc/solver/Simplex.cpp +++ /dev/null @@ -1,383 +0,0 @@ -/* - * SimplexSolver.cpp - * - * LP solution - - * Created on: 22 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#include "Simplex.hpp" -#include "SimplexTable.hpp" - -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "../Problem.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - - -namespace dtwc::solver { - -int getRow(const MatrixType &tableau, int index) -{ - if (index < 0 || index >= (tableau.cols() - 1)) // Check if index is in the valid range - throw std::runtime_error(fmt::format("The index of the variable ({}) must be between 0 and {}", index, tableau.cols() - 2)); - - int rowIndex = -1; // Using -1 to represent None - // So this is checking there is one and only one 1.0 element! -> #TODO change with something keeping book of basic variables in future. - for (int j = 0; j < tableau.rows(); ++j) - if (!isAround(tableau(j, index), 0)) // The entry is non zero - { - if (rowIndex == -1 && isAround(tableau(j, index), 1.0)) // The entry is one, and the index has not been found yet. - rowIndex = j; - else - return -1; - } - - return rowIndex; -} - -void pivoting(MatrixType &tableau, int p, int q) -{ - // Make p,q element one and eliminate all other nonzero elements in that column by basic row operations. - const double thepivot = tableau(q, p); - if (isAround(thepivot, 0.0)) - throw std::runtime_error(fmt::format("The pivot is too close to zero: {}", thepivot)); - - tableau.row(q) /= thepivot; // Make (p,q) one. - - for (int i = 0; i < tableau.rows(); ++i) - if (i != q) - tableau.row(i) -= tableau(i, p) * tableau.row(q); -} - -std::tuple simplexTableau(const MatrixType &tableau) -{ - const int mtab = tableau.rows(), ntab = tableau.cols(); - const int m = mtab - 1, n = ntab - 1; - - // Find the first negative cost, if there are none, then table is optimal. - int p = -1; - for (int i = 0; i < n; i++) // Reduced cost. - if (tableau(Eigen::last, i) < -epsilon) { - p = i; - break; - } - - if (p == -1) // The tableau is optimal - return std::make_tuple(-1, -1, true, true); - - double vkTolerance = 1e-10; - // Calculate the maximum step that can be done along the basic direction d[p] - double minStep = std::numeric_limits::infinity(); - int q{ -1 }; // row of min step. -1 means unbounded. - - for (int k = 0; k < m; k++) - if (tableau(k, p) > vkTolerance) { - const auto step = tableau(k, Eigen::last) / tableau(k, p); // tableau(k,p) = minus(d) - if (step < minStep) { - minStep = step; - q = k; - } - } - - if (q == -1) // The tableau is unbounded - return std::make_tuple(-1, -1, false, false); - - return std::make_tuple(p, q, false, true); -} - -std::pair simplexAlgorithmTableau(MatrixType &tableau) -{ - size_t iter{}; - double duration_table{}, duration_pivoting{}; - dtwc::Clock clk; - while (true) { - auto [colPivot, rowPivot, optimal, bounded] = simplexTableau(tableau); - duration_table += clk.duration(); - if (optimal) return { true, false }; - if (!bounded) return { false, true }; - duration_pivoting += clk.duration(); - - - pivoting(tableau, colPivot, rowPivot); - - if (iter % 100 == 0) { - std::cout << "Iteration " << iter << " is finished!\n"; - std::cout << "Duration table: "; - dtwc::Clock::print_duration(std::cout, duration_table); - std::cout << "Duration pivoting: "; - dtwc::Clock::print_duration(std::cout, duration_pivoting); - - duration_table = 0; - duration_pivoting = 0; - - - } - - iter++; - } -} - -MatrixType createTableau(const MatrixType &A, const VectorXd &b, VectorXd &c) -{ - const int m = A.rows(), n = A.cols(); - MatrixType tableau(m + 1, n + m + 1); - - tableau.block(0, 0, m, n) = A; - tableau.block(0, n, m, m) = MatrixType::Identity(m, m); - tableau.block(0, n + m, m, 1) = b; - - // Set the first n columns of the last row - for (int k = 0; k < n; ++k) - tableau.row(m)(k) = -tableau.block(0, k, m, 1).sum(); - - // Set columns n through n+m of the last row to 0.0 - tableau.block(m, n, 1, m).setZero(); - - // Set the last element of the last row - tableau(m, n + m) = -tableau.block(0, n + m, m, 1).sum(); - - return tableau; -} - -std::tuple simplex(MatrixType &A, VectorXd &b, VectorXd &c) -{ - // bool unbounded{}, optimal{}; - const int m = A.rows(), n = A.cols(); - - if (b.rows() != m) - throw std::runtime_error(fmt::format("Incompatible sizes: A is {}x{}, b is of length {}, and should be {}", m, n, b.rows(), m)); - - if (c.size() != n) - throw std::runtime_error(fmt::format("Incompatible sizes: A is {}x{}, c is of length {}, and should be {}", m, n, c.size(), n)); - - // Ensure all elements of b are non-negative - for (int i = 0; i < m; ++i) - if (b(i) < 0) { - A.row(i) = -A.row(i); - b(i) = -b(i); - } - - MatrixType phaseOneTableau = createTableau(A, b, c); - - auto [optimal, unbounded] = simplexAlgorithmTableau(phaseOneTableau); - - if (unbounded) return { phaseOneTableau, true, false }; - - const auto isInfeasible = phaseOneTableau(Eigen::last, Eigen::last) < -epsilon; - if (isInfeasible) return { phaseOneTableau, false, true }; // Infeasible problem - - std::vector basicRows(m + n); - std::set basicIndices; - std::set tobeCleaned; - - for (int k = 0; k < m + n; ++k) { - basicRows[k] = getRow(phaseOneTableau, k); - if (basicRows[k] != -1) { - basicIndices.insert(k); - if (k >= n) tobeCleaned.insert(k); // If k>= n are basic indices then clean them. - } - } - - while (!tobeCleaned.empty()) { - int auxiliaryColumn = *tobeCleaned.begin(); - tobeCleaned.erase(tobeCleaned.begin()); - - int rowpivotIndex = basicRows[auxiliaryColumn]; - VectorXd rowpivot = phaseOneTableau.row(rowpivotIndex); - - std::vector originalNonbasic; - for (int i = 0; i < n; ++i) - if (basicIndices.find(i) == basicIndices.end()) - originalNonbasic.push_back(i); - - int colpivot = -1; - double maxVal = std::numeric_limits::epsilon(); - - for (int col : originalNonbasic) - if (std::abs(rowpivot(col)) > maxVal) { - maxVal = std::abs(rowpivot(col)); - colpivot = col; - } - - if (colpivot != -1) { - pivoting(phaseOneTableau, colpivot, rowpivotIndex); - basicRows[colpivot] = rowpivotIndex; - basicRows[auxiliaryColumn] = -1; - } else { - phaseOneTableau.conservativeResize(phaseOneTableau.rows() - 1, phaseOneTableau.cols()); - for (int k = 0; k < m + n; ++k) - basicRows[k] = getRow(phaseOneTableau, k); - } - } - - MatrixType leftPart = phaseOneTableau.leftCols(n); - MatrixType rightPart = phaseOneTableau.rightCols(phaseOneTableau.cols() - n - m); - - MatrixType phaseTwoTableau(leftPart.rows(), leftPart.cols() + rightPart.cols()); - phaseTwoTableau << leftPart, rightPart; - - // Reset basicRows for startPhaseTwo -> phaseTwoTableau - basicRows.resize(n); - for (int k = 0; k < n; ++k) - basicRows[k] = getRow(phaseTwoTableau, k); - - // Calculate last row - std::vector basicIndicesList, nonbasicIndicesList; - for (int i = 0; i < n; ++i) { - if (basicRows[i] != -1) - basicIndicesList.push_back(i); - else - nonbasicIndicesList.push_back(i); - } - - for (int k : nonbasicIndicesList) { - double sumVal = 0.0; - for (int j : basicIndicesList) { - sumVal += c(j) * phaseTwoTableau(basicRows[j], k); - } - - phaseTwoTableau(phaseTwoTableau.rows() - 1, k) = c(k) - sumVal; - } - - double lastRowSum = 0.0; - for (int j : basicIndicesList) - lastRowSum += c(j) * phaseTwoTableau(basicRows[j], phaseTwoTableau.cols() - 1); - - phaseTwoTableau(phaseTwoTableau.rows() - 1, phaseTwoTableau.cols() - 1) = -lastRowSum; - - // Phase II - auto [optimal2, unbounded2] = simplexAlgorithmTableau(phaseTwoTableau); // it was startPhaseTwo - return { phaseTwoTableau, unbounded2, false }; -} - -std::pair, double> Simplex::getResults() const -{ - int n = tableau.cols() - 1; // finalTableau - - std::vector solution(n); - for (int k = 0; k < n; ++k) { - const auto basicRowNo = getRow(tableau, k); - solution[k] = (basicRowNo != -1) ? tableau(basicRowNo, Eigen::last) : 0.0; - } - - return { solution, -tableau(Eigen::last, Eigen::last) }; // Solution and optimal cost -} - - -void Simplex::gomory() -{ - fmt::println("=================================="); - fmt::println("Problem with {} variables and {} constraints", A.cols(), A.rows()); - - bool unbounded{}, infeasible{}; - std::tie(tableau, unbounded, infeasible) = simplex(A, b, c); // Assuming simplex is defined - - if (unbounded) { - fmt::println("Unbounded problem"); - nGomory = 0; - return; - } - - if (infeasible) { - fmt::println("Infeasible problem"); - nGomory = 0; - return; - } - - auto [solution, copt] = getResults(); // Assuming getResults is defined - if constexpr (settings::debug_Simplex) - fmt::println("Solution: {} and Copt = [{}]\n", solution, copt); - - std::vector fractionalMask; - - for (int i = 0; i < (tableau.rows() - 1); i++) - if (isFractional(tableau(i, Eigen::last))) // Scan last column - fractionalMask.push_back(i); - - nGomory = fractionalMask.size(); - fmt::println("Number of Gomory cuts: {}", nGomory); - - MatrixType gamma = tableau(fractionalMask, Eigen::seqN(0, tableau.cols() - 1)).array().floor(); - VectorXd bplus = tableau(fractionalMask, tableau.cols() - 1).array().floor(); - - MatrixType newA(A.rows() + nGomory, A.cols() + nGomory); - VectorXd newb(A.rows() + nGomory); - VectorXd newc(A.cols() + nGomory); - - newA << A, MatrixType::Zero(A.rows(), nGomory), - gamma, MatrixType::Identity(nGomory, nGomory); - - newb << b, bplus; - - newc << c, VectorXd::Zero(nGomory); - - A = newA; - b = newb; - c = newc; -} - -Simplex::Simplex(Problem &prob) -{ - const auto Nb = prob.data.size(); - const auto Nc = prob.cluster_size(); - - const auto Neq = Nb + 1; - const auto Nineq = Nb * (Nb - 1); - const auto Nconstraints = Neq + Nineq; - - const auto Nvar_original = Nb * Nb; - const auto N_slack = Nineq; - const auto Nvar = Nvar_original + N_slack; // x1--xN^2 + s_slack - - A = MatrixType::Zero(Nconstraints, Nvar); - b = VectorXd::Zero(Nconstraints); - c = VectorXd::Zero(Nvar); - - // Create A matrix: - A.bottomRightCorner(N_slack, N_slack) = MatrixType::Identity(N_slack, N_slack); - - for (int i = 0; i < Nb; ++i) { - A(0, i * (Nb + 1)) = 1.0; // Sum of diagonals is Nc - A.block(1, Nb * i, Nb, Nb) = MatrixType::Identity(Nb, Nb); // Every element belongs to one cluster. - - // --------------- - int shift = 0; - for (int j = 0; j < Nb; j++) { - const int block_begin_row = Nb + 1 + (Nb - 1) * i; - const int block_begin_col = Nb * i; - if (i == j) { - A.block(block_begin_row, block_begin_col + j, Nb - 1, 1) = -1 * MatrixType::Ones(Nb - 1, 1); - shift = 1; - } else - A(block_begin_row + j - shift, block_begin_col + j) = 1; - } - } - - // Create b matrix: - b(0) = Nc; - for (int i = 0; i < Nb; ++i) - b(i + 1) = 1; - - for (size_t j{ 0 }; j < Nb; j++) - for (size_t i{ 0 }; i < Nb; i++) - c[i + j * Nb] = prob.distByInd_scaled(i, j); -} - -} // namespace dtwc::solver diff --git a/dtwc/solver/Simplex.hpp b/dtwc/solver/Simplex.hpp deleted file mode 100644 index 9087763..0000000 --- a/dtwc/solver/Simplex.hpp +++ /dev/null @@ -1,76 +0,0 @@ -/* - * Simplex.hpp - * - * LP solution - - * Created on: 22 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#pragma once - -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - - -namespace dtwc { -class Problem; -} - -namespace dtwc::solver { - -using Eigen::MatrixXd, Eigen::VectorXd; -using MatrixTypeRowMajor = Eigen::Matrix; -using MatrixType = MatrixTypeRowMajor; // Faster than column-major - - -class Simplex -{ - // -1 means None. - size_t N{ 0 }, Nc{ 0 }; - int mtab{}, ntab{}; - MatrixType tableau; - MatrixType A; // A*x = b - VectorXd b, c; // c*x = cost. - int nGomory{ -1 }; - -public: - Simplex(MatrixType A_, VectorXd b_, VectorXd c_) : A(A_), b(b_), c(c_) {} - Simplex() = default; - Simplex(Problem &prob); - - void gomory(); - void gomoryAlgorithm() - { - while (nGomory != 0) gomory(); - } - - std::pair, double> getResults() const; -}; - - -int getRow(const MatrixType &tableau, int index); -void inline pivoting(MatrixType &tableau, int p, int q); -std::tuple inline simplexTableau(const MatrixType &tableau); -std::pair inline simplexAlgorithmTableau(MatrixType &input_tableau); -MatrixType inline createTableau(const MatrixType &A, const VectorXd &b, VectorXd &c); -std::tuple inline simplex(MatrixType &A, VectorXd &b, VectorXd &c); - -} // namespace dtwc::solver diff --git a/dtwc/solver/SimplexFlatRowTable.cpp b/dtwc/solver/SimplexFlatRowTable.cpp deleted file mode 100644 index 69d9f6b..0000000 --- a/dtwc/solver/SimplexFlatRowTable.cpp +++ /dev/null @@ -1,282 +0,0 @@ -/* - * SimplexFlatRowTable.cpp - * - * Sparse implementation of a Simplex table. - - * Created on: 29 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#include "SimplexFlatRowTable.hpp" -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "EqualityConstraints.hpp" -#include -#include -#include "../timing.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - - -namespace dtwc::solver { - - -int SimplexFlatRowTable::getRow(int col) const -{ - if (col < 0 || col >= (ntab - 1)) // Check if index is in the valid range - throw std::runtime_error(fmt::format("The index of the variable ({}) must be between 0 and {}", col, ntab - 2)); - - int rowIndex = -1; // Using -1 to represent None - // So this is checking there is one and only one 1.0 element! -> #TODO change with something keeping book of basic variables in future. - - if (!isAround(reducedCosts[col], 0.0)) return rowIndex; - - for (int row = 0; row < innerTable.size(); row++) { - for (const auto [key, val] : innerTable[row]) - if (key == col) { // The entry is non-zero - if (rowIndex == -1 && isAround(val, 1.0)) // The entry is one, and the index has not been found yet. - rowIndex = row; - else - return -1; - } else if (key > col) - break; - } - - return rowIndex; -} - -int SimplexFlatRowTable::findMostNegativeCost() const -{ - // Returns -1 if table is optimal, otherwise the index of the most negative reduced cost. - int p = -1; - double mostNegativeValue = -epsilon; // start with a small negative threshold value - - for (int i = 0; i < reducedCosts.size(); i++) { - if (reducedCosts[i] < mostNegativeValue) { - mostNegativeValue = reducedCosts[i]; - p = i; - } - } - - return p; -} - -int SimplexFlatRowTable::findNegativeCost() const -{ - // Returns -1 if table is optimal otherwise first negative reduced cost. - // Find the first negative cost, if there are none, then table is optimal. - int p = -1; - for (int i = 0; i < reducedCosts.size(); i++) // Reduced cost. - if (reducedCosts[i] < -epsilon) { - p = i; - break; - } - - return p; -} - -int SimplexFlatRowTable::findMinStep(int p) -{ - pivotRows.clear(); // Prep for findMinStep; - pivotRows.resize(innerTable.size(), { -1, -1 }); - // Calculate the maximum step that can be done along the basic direction d[p] - using StepIndexPair = std::pair; - - - StepIndexPair initial = { std::numeric_limits::infinity(), -1 }; // row of min step. -1 means unbounded. - - - auto result = std::transform_reduce( - std::execution::par_unseq, // Use parallel execution policy - innerTable.begin(), - innerTable.end(), - initial, - [](const StepIndexPair &a, const StepIndexPair &b) -> StepIndexPair { - if (isAround(a.first, b.first, 1e-12)) - return (a.second < b.second) ? a : b; - else - return (a.first < b.first) ? a : b; - }, - [p, this](const auto &row) -> StepIndexPair { - const int rowIndex = (&row - &innerTable[0]); - double minus_d = -1; - for (int colIndex{}; colIndex < row.size(); colIndex++) - if (row[colIndex].index == p) { - pivotRows[rowIndex] = { colIndex, row[colIndex].value }; - minus_d = row[colIndex].value; - break; - } else if (row[colIndex].index > p) - break; - - - if (minus_d > epsilon) - return { rhs[rowIndex] / minus_d, rowIndex }; - - return { std::numeric_limits::infinity(), -1 }; - }); - - return result.second; -} - -std::tuple SimplexFlatRowTable::simplexTableau() -{ - // Find the first negative cost, if there are none, then table is optimal. - std::uniform_real_distribution dist(0.0, 1.0); - // Find the first negative cost, if there are none, then table is optimal. - // const int p = (dist(randGenerator) < 0.5) ? findNegativeCost() : findMostNegativeCost(); - const int p = findNegativeCost(); - if (p == -1) // The table is optimal - return std::make_tuple(-1, -1, true, true); - - const int q = findMinStep(p); - - if (q == -1) // The table is unbounded - return std::make_tuple(-1, -1, false, false); - - return std::make_tuple(p, q, false, true); -} - -void SimplexFlatRowTable::pivoting(int p, int q) -{ - constexpr double tol_sparsity = epsilon / 1000; - // p column, q = row. - // Make p,q element one and eliminate all other nonzero elements in that column by basic row operations. - const double thepivot = pivotRows[q].value; - if (isAround(thepivot, 0.0)) - throw std::runtime_error(fmt::format("The pivot is too close to zero: {}", thepivot)); - - auto &pivotRow = innerTable[q]; - - const auto reducedcost_p = reducedCosts[p]; - rhs[q] /= thepivot; - negativeObjective -= rhs[q] * reducedcost_p; - - for (auto &[key, val] : pivotRow) // Make the pivot row normalised. - { - val /= thepivot; - reducedCosts[key] -= reducedcost_p * val; // Remove from last row. - } - - auto oneTask = [this, thepivot, p, q, &pivotRow](auto &rowNow) { - if (&rowNow == &pivotRow) return; // Do not process pivot row. - - const int rowIndex = &rowNow - &innerTable[0]; - - if (pivotRows[rowIndex].index == -1) return; - - const auto p_val = pivotRows[rowIndex].value; - - rhs[rowIndex] = std::max(rhs[rowIndex] - p_val * rhs[q], 0.0); - - - const auto N_now = rowNow.size(); - const auto N_piv = pivotRow.size(); - - thread_local std::vector temporary; - temporary.clear(); - temporary.reserve(N_now); - - size_t i_now{}, i_piv{}; - - while (i_now != N_now && i_piv != N_piv) { - const auto &eNow = rowNow[i_now]; - const auto &ePiv = pivotRow[i_piv]; - - if (eNow.index < ePiv.index) { - temporary.push_back(eNow); - ++i_now; - } else if (eNow.index == ePiv.index) { - const auto new_value = eNow.value - p_val * ePiv.value; - if (!isAround(new_value, 0.0, tol_sparsity)) - temporary.emplace_back(Element{ eNow.index, new_value }); - - ++i_now; - ++i_piv; - } else { - const auto new_value = -p_val * ePiv.value; - if (!isAround(new_value, 0.0, tol_sparsity)) - temporary.emplace_back(ePiv.index, new_value); - ++i_piv; - } - } - - while (i_now != N_now) - temporary.push_back(rowNow[i_now++]); - - while (i_piv != N_piv) { - const auto &ePiv = pivotRow[i_piv++]; - const auto new_value = -p_val * ePiv.value; - if (!isAround(new_value, 0.0, tol_sparsity)) - temporary.emplace_back(ePiv.index, new_value); - } - - std::swap(temporary, rowNow); - }; - - std::for_each(std::execution::par_unseq, innerTable.begin(), innerTable.end(), oneTask); -} - - -std::pair SimplexFlatRowTable::simplexAlgorithmTableau() -{ - size_t iter{}; - double duration_table{}, duration_pivoting{}; - // static std::vector colNumbers; - // colNumbers.reserve(rhs.size()); - - - while (true) { - dtwc::Clock clk; - // colNumbers.clear(); - auto [colPivot, rowPivot, optimal, bounded] = simplexTableau(); - duration_table += clk.duration(); - - if (optimal) return { true, false }; - if (!bounded) return { false, true }; - - clk = dtwc::Clock{}; - - pivoting(colPivot, rowPivot); - duration_pivoting += clk.duration(); - - - if (iter % 500 == 0) { - std::cout << "Iteration " << iter << " is finished!\n"; - std::cout << "Duration table: "; - dtwc::Clock::print_duration(std::cout, duration_table); - std::cout << "Duration pivoting: "; - dtwc::Clock::print_duration(std::cout, duration_pivoting); - - duration_table = 0; - duration_pivoting = 0; - - size_t innerSize = 0; - for (auto &map : innerTable) { - innerSize += map.size(); - } - - std::cout << "Inner size per row: " << (double)innerSize / innerTable.size() - << " per " << innerTable.size() << '\n'; - - std::cout << "Minimum reduced cost: " << *std::min_element(reducedCosts.begin(), reducedCosts.end()) << '\n'; - std::cout << "Minimum rhs: " << *std::min_element(rhs.begin(), rhs.end()) << '\n'; - } - - iter++; - } -} - -} // namespace dtwc::solver diff --git a/dtwc/solver/SimplexFlatRowTable.hpp b/dtwc/solver/SimplexFlatRowTable.hpp deleted file mode 100644 index 05a8f77..0000000 --- a/dtwc/solver/SimplexFlatRowTable.hpp +++ /dev/null @@ -1,143 +0,0 @@ -/* - * SimplexFlatRowTable.hpp - * - * Sparse implementation of a Simplex table but this time each row. - * Using vectors for a flat map instead of std::map to avoid allocations. - * Created on: 30 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#pragma once - -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "sparse_util.hpp" -#include "EqualityConstraints.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include -#include - - -namespace dtwc::solver { - - -class SimplexFlatRowTable -{ - // Table is mtab x ntab - int mtab{}, ntab{}; - // Inner table is m x n - std::vector> innerTable; // Each is a row - std::vector reducedCosts, rhs; - double negativeObjective{}; - - std::vector pivotRows; // Literal Index of pivotRow / value of pivot col there. - -public: - SimplexFlatRowTable() = default; - SimplexFlatRowTable(int mtab_, int ntab_) : mtab{ mtab_ }, ntab{ ntab_ }, innerTable(mtab_ - 1), - reducedCosts(ntab_ - 1), rhs(mtab_ - 1) {} - - int rows() const { return mtab; } - int cols() const { return ntab; } - - void clear() - { - rhs.clear(); - reducedCosts.clear(); - - for (auto &map : innerTable) // Should not remove allocated memory for maps. - map.clear(); - } - - void removeColumns(int a, int b) - { - // Removes columns [a, b) - for (auto &row : innerTable) { - auto itBegin = std::lower_bound(row.begin(), row.end(), a, [](const Element &elem, int idx) { - return elem.index < idx; - }); - auto itEnd = std::upper_bound(itBegin, row.end(), b, [](int idx, const Element &elem) { - return idx < elem.index; - }); - - row.erase(itBegin, itEnd); - } - - reducedCosts.erase(reducedCosts.begin() + a, reducedCosts.begin() + b); - ntab -= (b - a); - } - - - void createPhaseOneTableau(const EqualityConstraints &eq) - { - const int m = eq.A.rows(), n = eq.A.cols(); - mtab = m + 1; - ntab = m + n + 1; - clear(); // Clear the things. - innerTable.resize(m); // Adding m auxillary variables. - rhs = eq.b; - reducedCosts.resize(n + m, 0.0); // Set columns n through n+m of the last row to 0.0 - - for (const auto [key, val] : eq.A.data) { - const auto [i, j] = key; - innerTable[i].push_back(Element{ j, val }); - reducedCosts[j] -= val; // Set the first n columns of the last row - } - - for (int i = 0; i < m; i++) - innerTable[i].emplace_back(n + i, 1.0); - - negativeObjective = -std::reduce(rhs.begin(), rhs.end()); - - std::for_each(innerTable.begin(), innerTable.end(), [](auto &elemVec) { - std::sort(elemVec.begin(), elemVec.end(), CompElementIndices{}); - }); - } - - int getRow(int col) const; - double getObjective() const { return -negativeObjective; } - double getRHS(int k) const { return rhs[k]; } - double getReducedCost(int k) const { return reducedCosts[k]; } - - double getValue(int k) const - { - const auto basicRowNo = getRow(k); - return (basicRowNo != -1) ? rhs[basicRowNo] : 0.0; - } - - double inner(int i, int j) - { - for (const auto [key, val] : innerTable[i]) - if (key == j) - return val; - else if (key > j) - return 0.0; - - return 0.0; - } - - double &setReducedCost(int k) { return reducedCosts[k]; } - void setNegativeObjective(double val) { negativeObjective = val; } - - int findNegativeCost() const; - int findMostNegativeCost() const; - int findMinStep(int p); - - std::tuple simplexTableau(); - - void pivoting(int p, int q); - - std::pair simplexAlgorithmTableau(); -}; - -} // namespace dtwc::solver diff --git a/dtwc/solver/SimplexFlatTable.cpp b/dtwc/solver/SimplexFlatTable.cpp deleted file mode 100644 index aede0b2..0000000 --- a/dtwc/solver/SimplexFlatTable.cpp +++ /dev/null @@ -1,273 +0,0 @@ -/* - * SimplexFlatTable.cpp - * - * Sparse implementation of a Simplex table. - - * Created on: 29 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#include "SimplexFlatTable.hpp" -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "EqualityConstraints.hpp" -#include -#include -#include "../timing.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - - -namespace dtwc::solver { - - -int SimplexFlatTable::getRow(int col) const -{ - if (col < 0 || col >= (ntab - 1)) // Check if index is in the valid range - throw std::runtime_error(fmt::format("The index of the variable ({}) must be between 0 and {}", col, ntab - 2)); - - // So this is checking there is one and only one 1.0 element! -> #TODO change with something keeping book of basic variables in future. - const bool isBasic = isAround(reducedCosts[col], 0.0) && innerTable[col].size() == 1 && isAround(innerTable[col][0].value, 1); - return isBasic ? innerTable[col][0].index : -1; // Using -1 to represent None -} - -int SimplexFlatTable::findNegativeCost() -{ - // Returns -1 if table is optimal otherwise first negative reduced cost. - // Find the first negative cost, if there are none, then table is optimal. - int p = -1; - for (int i = 0; i < reducedCosts.size(); i++) // Reduced cost. - { - const auto i_search = i % reducedCosts.size(); - if (reducedCosts[i_search] < -epsilon) { - p = i_search; - break; - } - } - - - return p; -} - -int SimplexFlatTable::findMostNegativeCost() const -{ - // Returns -1 if table is optimal, otherwise the index of the most negative reduced cost. - int p = -1; - double mostNegativeValue = -epsilon; // start with a small negative threshold value - - for (int i = 0; i < reducedCosts.size(); i++) { - if (reducedCosts[i] < mostNegativeValue) { - mostNegativeValue = reducedCosts[i]; - p = i; - } - } - - return p; -} - -int SimplexFlatTable::findMinStep(int p) const -{ - // Calculate the maximum step that can be done along the basic direction d[p] - double minStep = std::numeric_limits::infinity(); - int q{ -1 }; // row of min step. -1 means unbounded. - - for (auto [k, minus_d] : innerTable[p]) - if (minus_d > epsilon) { - const auto step = rhs[k] / minus_d; - if (step < minStep || minStep < 0) { - minStep = step; - q = k; - if (minStep == 0) return q; - } - } - - return q; -} - -std::tuple SimplexFlatTable::simplexTableau() -{ - std::uniform_real_distribution dist(0.0, 1.0); - // Find the first negative cost, if there are none, then table is optimal. - const int p = (dist(randGenerator) < 0.5) ? findNegativeCost() : findMostNegativeCost(); - - if (p == -1) // The table is optimal - return std::make_tuple(-1, -1, true, true); - - const int q = findMinStep(p); - - if (q == -1) // The table is unbounded - return std::make_tuple(-1, -1, false, false); - - return std::make_tuple(p, q, false, true); -} - -void SimplexFlatTable::pivoting(int p, int q) -{ - constexpr double tol_sparsity = 1e-15; - // p column, q = row. - // Make p,q element one and eliminate all other nonzero elements in that column by basic row operations. - const double thepivot = inner(q, p); - if (isAround(thepivot, 0.0, tol_sparsity)) { - - std::cout << "The pivot is too close to zero for (" << q << "," << p << ')' << std::endl; - for (auto elem : innerTable[p]) - std::cout << "(" << elem.index << ',' << elem.value << ")\n"; - - std::cout << std::endl; - - throw std::runtime_error(fmt::format("The pivot is too close to zero: {}", thepivot)); - } - - auto &pivotCol = innerTable[p]; - - auto oneTask = [this, thepivot, p, q, &pivotCol](auto &colNow) { - // if (!std::is_sorted(colNow.begin(), colNow.end(), CompElementIndices{})) { - // std::cout << "Series is not sorted" << std::endl; - // throw "Given column is not sorted!"; - // } - - if (&colNow == &pivotCol) return; // Do not process pivot row. - - const int colIndex = &colNow - &innerTable[0]; - - auto it_q = std::lower_bound(colNow.begin(), colNow.end(), q, [](const Element &elem, int idx) { - return elem.index < idx; - }); - - if (it_q == colNow.end() || it_q->index != q) return; // We don't have that index. - - (it_q->value) /= thepivot; // Normalise by pivot; - - const auto q_val = (it_q->value); - - reducedCosts[colIndex] -= reducedCosts[p] * q_val; // Remove from last row. - - // if (reducedCosts[colIndex] < -1e20) { - // std::cout << "Reduced cost is too low! " << reducedCosts[colIndex] << std::endl; - // throw "Reduced cost is too low!!!!\n"; - // } - - auto N_now = colNow.size(); - auto N_piv = pivotCol.size(); - - thread_local std::vector temporary; - temporary.clear(); - size_t i_now{}, i_piv{}; - - while (i_now != N_now && i_piv != N_piv) { - const auto &eNow = colNow[i_now]; - const auto &ePiv = pivotCol[i_piv]; - - if (ePiv.index == q) - ++i_piv; - else if (eNow.index < ePiv.index) { - // if (!isAround(eNow.value, 0.0, 10 * tol_sparsity)) - temporary.push_back(eNow); - ++i_now; - } else if (eNow.index == ePiv.index) { - const auto new_value = eNow.value - q_val * ePiv.value; - // if (!isAround(new_value, 0.0, 10 * tol_sparsity)) - temporary.emplace_back(Element{ eNow.index, new_value }); - ++i_now; - ++i_piv; - } else { - const auto new_value = -q_val * ePiv.value; - // if (!isAround(new_value, 0.0, 10 * tol_sparsity)) - temporary.emplace_back(ePiv.index, new_value); - ++i_piv; - } - } - - while (i_now != N_now) - // if (!isAround(colNow[i_now].value)) - temporary.push_back(colNow[i_now++]); - - while (i_piv != N_piv) { - const auto &ePiv = pivotCol[i_piv++]; - const auto new_value = -q_val * ePiv.value; - // if (!isAround(new_value, 0.0, 10 * tol_sparsity)) - temporary.emplace_back(ePiv.index, new_value); - } - - std::swap(temporary, colNow); - }; - - std::for_each(std::execution::par_unseq, innerTable.begin(), innerTable.end(), oneTask); - - rhs[q] /= thepivot; - // We always have RHS. - for (auto [key, val] : innerTable[p]) - if (key != q) { - const auto new_value = rhs[key] - val * rhs[q]; - rhs[key] = new_value; // std::max(new_value, 0.0); - } - - - negativeObjective -= rhs[q] * reducedCosts[p]; - - // Deal with the pivot column now. - pivotCol.clear(); - pivotCol.emplace_back(q, 1); - reducedCosts[p] = 0; -} - - -std::pair SimplexFlatTable::simplexAlgorithmTableau() -{ - size_t iter{}; - double duration_table{}, duration_pivoting{}; - while (true) { - dtwc::Clock clk; - auto [colPivot, rowPivot, optimal, bounded] = simplexTableau(); - duration_table += clk.duration(); - - if (optimal) return { true, false }; - if (!bounded) return { false, true }; - - clk = dtwc::Clock{}; - - pivoting(colPivot, rowPivot); - duration_pivoting += clk.duration(); - - - if (iter % 500 == 0) { - std::cout << "Iteration " << iter << " is finished!" << std::endl; - std::cout << "Duration table: "; - dtwc::Clock::print_duration(std::cout, duration_table); - std::cout << "Duration pivoting: "; - dtwc::Clock::print_duration(std::cout, duration_pivoting); - - duration_table = 0; - duration_pivoting = 0; - - size_t innerSize = 0; - for (auto &map : innerTable) { - innerSize += map.size(); - } - - std::cout << "Inner size per row: " << (double)innerSize / innerTable.size() - << " per " << innerTable.size() << '\n'; - - - std::cout << "Minimum reduced cost: " << *std::min_element(reducedCosts.begin(), reducedCosts.end()) << '\n'; - std::cout << "Minimum rhs: " << *std::min_element(rhs.begin(), rhs.end()) << '\n'; - } - - iter++; - } -} - -} // namespace dtwc::solver diff --git a/dtwc/solver/SimplexFlatTable.hpp b/dtwc/solver/SimplexFlatTable.hpp deleted file mode 100644 index 8fda1c9..0000000 --- a/dtwc/solver/SimplexFlatTable.hpp +++ /dev/null @@ -1,138 +0,0 @@ -/* - * SimplexFlatTable.hpp - * - * Sparse implementation of a Simplex table but this time each row. - * Using vectors for a flat map instead of std::map to avoid allocations. - * Created on: 30 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#pragma once - -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "sparse_util.hpp" -#include "EqualityConstraints.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include -#include - - -namespace dtwc::solver { - - -class SimplexFlatTable -{ - // Table is mtab x ntab - int mtab{}, ntab{}; - // Inner table is m x n - std::vector> innerTable; // Each is a column - std::vector reducedCosts, rhs; - double negativeObjective{}; - - int iLastPivot{ -1 }; // Last pivoted column - -public: - SimplexFlatTable() = default; - SimplexFlatTable(int mtab_, int ntab_) : mtab{ mtab_ }, ntab{ ntab_ }, innerTable(ntab_ - 1), - reducedCosts(ntab_ - 1), rhs(mtab_ - 1) {} - - int rows() const { return mtab; } - int cols() const { return ntab; } - - void clear() - { - rhs.clear(); - reducedCosts.clear(); - - for (auto &map : innerTable) // Should not remove allocated memory for maps. - map.clear(); - } - - void removeColumns(int a, int b) - { - // Removes columns [a, b) - innerTable.erase(innerTable.begin() + a, innerTable.begin() + b); - reducedCosts.erase(reducedCosts.begin() + a, reducedCosts.begin() + b); - ntab -= (b - a); - } - - - void createPhaseOneTableau(const EqualityConstraints &eq) - { - const int m = eq.A.rows(), n = eq.A.cols(); - mtab = m + 1; - ntab = m + n + 1; - clear(); // Clear the things. - innerTable.resize(m + n); // Adding m auxillary variables. - rhs = eq.b; - reducedCosts.resize(n + m, 0.0); // Set columns n through n+m of the last row to 0.0 - - for (const auto [key, val] : eq.A.data) { - const auto [i, j] = key; - innerTable[j].emplace_back(i, val); - reducedCosts[j] -= val; // Set the first n columns of the last row - } - - for (int i = 0; i < m; i++) - innerTable[n + i].emplace_back(i, 1.0); - - negativeObjective = -std::reduce(rhs.begin(), rhs.end()); - - std::for_each(innerTable.begin(), innerTable.end(), [](auto &elemVec) { - std::sort(elemVec.begin(), elemVec.end(), CompElementIndices{}); - } - - ); - } - - int getRow(int col) const; - - double getObjective() const { return -negativeObjective; } - - double getValue(int k) const - { - const auto basicRowNo = getRow(k); - return (basicRowNo != -1) ? rhs[basicRowNo] : 0.0; - } - - double getRHS(int k) const { return rhs[k]; } - - double getReducedCost(int k) const { return reducedCosts[k]; } - - double inner(int j, int i) - { - const auto &vec = innerTable[i]; - auto it = std::lower_bound(vec.begin(), vec.end(), j, [](const Element &elem, int idx) { - return elem.index < idx; - }); - - if (it != vec.end() && it->index == j) - return it->value; - - return 0.0; - } - - double &setReducedCost(int k) { return reducedCosts[k]; } - void setNegativeObjective(double val) { negativeObjective = val; } - - int findNegativeCost(); - int findMinStep(int p) const; - int findMostNegativeCost() const; - - void pivoting(int p, int q); - std::tuple simplexTableau(); - std::pair simplexAlgorithmTableau(); -}; - -} // namespace dtwc::solver diff --git a/dtwc/solver/SimplexRowTable.cpp b/dtwc/solver/SimplexRowTable.cpp deleted file mode 100644 index 7d05d61..0000000 --- a/dtwc/solver/SimplexRowTable.cpp +++ /dev/null @@ -1,231 +0,0 @@ -/* - * SimplexRowTable.cpp - * - * Sparse implementation of a Simplex table. - - * Created on: 29 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#include "SimplexRowTable.hpp" -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "EqualityConstraints.hpp" -#include -#include -#include "../timing.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - - -namespace dtwc::solver { - - -int SimplexRowTable::getRow(int col) const -{ - if (col < 0 || col >= (ntab - 1)) // Check if index is in the valid range - throw std::runtime_error(fmt::format("The index of the variable ({}) must be between 0 and {}", col, ntab - 2)); - - int rowIndex = -1; // Using -1 to represent None - // So this is checking there is one and only one 1.0 element! -> #TODO change with something keeping book of basic variables in future. - - if (!isAround(reducedCosts[col], 0.0)) return -1; - - for (int row = 0; row < innerTable.size(); row++) { - auto it = innerTable[row].find(col); - if (it != innerTable[row].end()) { - if (!isAround(it->second, 0)) // The entry is non zero - { - if (rowIndex == -1 && isAround(it->second, 1.0)) // The entry is one, and the index has not been found yet. - rowIndex = row; - else - return -1; - } - } - } - - return rowIndex; -} - -int SimplexRowTable::findNegativeCost() -{ - // Returns -1 if table is optimal otherwise first negative reduced cost. - // Find the first negative cost, if there are none, then table is optimal. - int p = -1; - for (int i = 0; i < reducedCosts.size(); i++) // Reduced cost. - if (reducedCosts[i] < -epsilon) { - p = i; - break; - } - - return p; -} - -int SimplexRowTable::findMinStep(int p) -{ - - // Calculate the maximum step that can be done along the basic direction d[p] - - constexpr double tol = 1e-10; - using StepIndexPair = std::pair; - StepIndexPair initial = { std::numeric_limits::infinity(), -1 }; // row of min step. -1 means unbounded. - - auto result = std::transform_reduce( - std::execution::par_unseq, // Use parallel execution policy - innerTable.begin(), - innerTable.end(), - initial, - [](const StepIndexPair &a, const StepIndexPair &b) -> StepIndexPair { - return (a.first < b.first) ? a : b; - }, - [p, tol, this](const auto &row) -> StepIndexPair { - auto it = row.find(p); - if (it != row.end()) { - auto [k, minus_d] = *it; - if (minus_d > tol) { - int rowIndex = std::distance(innerTable.begin(), std::find(innerTable.begin(), innerTable.end(), row)); - return { rhs[rowIndex] / minus_d, rowIndex }; - } - } - return { std::numeric_limits::infinity(), -1 }; - }); - - return result.second; -} - -std::tuple SimplexRowTable::simplexTableau() -{ - // Find the first negative cost, if there are none, then table is optimal. - const int p = findNegativeCost(); - - if (p == -1) // The table is optimal - return std::make_tuple(-1, -1, true, true); - - const int q = findMinStep(p); - - if (q == -1) // The table is unbounded - return std::make_tuple(-1, -1, false, false); - - return std::make_tuple(p, q, false, true); -} - -void SimplexRowTable::pivoting(int p, int q) -{ - // p column, q = row. - // Make p,q element one and eliminate all other nonzero elements in that column by basic row operations. - const double thepivot = innerTable[q][p]; - if (isAround(thepivot, 0.0)) - throw std::runtime_error(fmt::format("The pivot is too close to zero: {}", thepivot)); - - auto &pivotRow = innerTable[q]; - - const auto reducedcost_p = reducedCosts[p]; - rhs[q] /= thepivot; - negativeObjective -= rhs[q] * reducedcost_p; - - for (auto &[key, val] : pivotRow) // Make the pivot row normalised. - { - val /= thepivot; - reducedCosts[key] -= reducedcost_p * val; // Remove from last row. - } - - auto oneTask = [this, thepivot, p, q, pivotRow](int i) { - if (q == i) return; // Do not process pivot row. - - auto &rowNow = innerTable[i]; - - auto it_p = rowNow.find(p); - if (it_p == rowNow.end()) return; - - const auto p_val = it_p->second; - auto it_now = rowNow.begin(); - auto it_pivot = pivotRow.begin(); - - rhs[i] -= p_val * rhs[q]; - - while (it_now != rowNow.end() && it_pivot != pivotRow.end()) { - const auto [key_piv, val_piv] = (*it_pivot); - auto &[key_now, val_now] = (*it_now); - - if (key_now < key_piv) - ++it_now; - else if (key_now == key_piv) { - val_now -= p_val * val_piv; - ++it_now; - ++it_pivot; - } else { - rowNow.insert(it_now, { key_piv, -(p_val * val_piv) }); - ++it_pivot; - } - } - - while (it_pivot != pivotRow.end()) { - rowNow[it_pivot->first] = -p_val * (it_pivot->second); - ++it_pivot; - } - - - std::erase_if(rowNow, [](const auto &item) { - auto const &[key, value] = item; - return isAround(value, 0.0); // Remove zero elements to make it compact. - }); - }; - - dtwc::run(oneTask, innerTable.size()); -} - - -std::pair SimplexRowTable::simplexAlgorithmTableau() -{ - size_t iter{}; - double duration_table{}, duration_pivoting{}; - while (true) { - dtwc::Clock clk; - auto [colPivot, rowPivot, optimal, bounded] = simplexTableau(); - duration_table += clk.duration(); - - if (optimal) return { true, false }; - if (!bounded) return { false, true }; - - clk = dtwc::Clock{}; - - pivoting(colPivot, rowPivot); - duration_pivoting += clk.duration(); - - - if (iter % 100 == 0) { - std::cout << "Iteration " << iter << " is finished!" << std::endl; - std::cout << "Duration table: "; - dtwc::Clock::print_duration(std::cout, duration_table); - std::cout << "Duration pivoting: "; - dtwc::Clock::print_duration(std::cout, duration_pivoting); - - duration_table = 0; - duration_pivoting = 0; - - size_t innerSize = 0; - for (auto &map : innerTable) { - innerSize += map.size(); - } - - std::cout << "Inner size: " << innerSize << '\n'; - } - - iter++; - } -} - -} // namespace dtwc::solver diff --git a/dtwc/solver/SimplexRowTable.hpp b/dtwc/solver/SimplexRowTable.hpp deleted file mode 100644 index 07830b0..0000000 --- a/dtwc/solver/SimplexRowTable.hpp +++ /dev/null @@ -1,141 +0,0 @@ -/* - * SimplexRowTable.hpp - * - * Sparse implementation of a Simplex table but this time each row. - - * Created on: 22 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#pragma once - -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "EqualityConstraints.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include - - -namespace dtwc::solver { - -class SimplexRowTable -{ - // Table is mtab x ntab - // Inner table is m x n - int mtab{}, ntab{}; - std::vector> innerTable; // Each is a row - std::vector reducedCosts, rhs; - std::vector rowIndices; - double negativeObjective{}; - -public: - SimplexRowTable() = default; - SimplexRowTable(int mtab_, int ntab_) : mtab{ mtab_ }, ntab{ ntab_ }, innerTable(mtab - 1), - reducedCosts(ntab_ - 1), rhs(mtab_ - 1), rowIndices(ntab_ - 1) {} - - int rows() const { return mtab; } - int cols() const { return ntab; } - - void clear() - { - rowIndices.clear(); - rhs.clear(); - reducedCosts.clear(); - - for (auto &map : innerTable) // Should not remove allocated memory for maps. - map.clear(); - } - - void removeColumns(int a, int b) - { - // Removes columns [a, b) - for (auto &row : innerTable) { - auto itBegin = row.lower_bound(a); - auto itEnd = row.upper_bound(b); - row.erase(itBegin, itEnd); - } - - reducedCosts.erase(reducedCosts.begin() + a, reducedCosts.begin() + b); - ntab -= (b - a); - } - - - void createPhaseOneTableau(const EqualityConstraints &eq) - { - // Table size (m + 1, m + n + 1) - // table.block(0, 0, m, n) = A; - // table.block(0, n, m, m) = SimplexRowTable::Identity(m, m); - // table.block(0, n + m, m, 1) = b; - - // // Set the first n columns of the last row - // for (int k = 0; k < n; ++k) - // table.row(m)(k) = -table.block(0, k, m, 1).sum(); - - // // Set columns n through n+m of the last row to 0.0 - // table.block(m, n, 1, m).setZero(); - - // // Set the last element of the last row - // table(m, n + m) = -table.block(0, n + m, m, 1).sum(); - - const int m = eq.A.rows(), n = eq.A.cols(); - mtab = m + 1; - ntab = m + n + 1; - clear(); // Clear the things. - innerTable.resize(m); // Adding m auxillary variables. - rhs = eq.b; - reducedCosts.resize(n + m, 0.0); // Set columns n through n+m of the last row to 0.0 - - for (const auto [key, val] : eq.A.data) { - const auto [i, j] = key; - innerTable[i][j] = val; - reducedCosts[j] -= val; // Set the first n columns of the last row - } - - for (int i = 0; i < m; i++) - innerTable[i][n + i] = 1.0; - - negativeObjective = -std::reduce(rhs.begin(), rhs.end()); - } - - int getRow(int col) const; - double getObjective() const { return -negativeObjective; } - - double getValue(int k) const - { - const auto basicRowNo = getRow(k); - return (basicRowNo != -1) ? rhs[basicRowNo] : 0.0; - } - - double getRHS(int k) const { return rhs[k]; } - - double getReducedCost(int k) const { return reducedCosts[k]; } - - double inner(int i, int j) { return innerTable[i][j]; } - - - double &setReducedCost(int k) { return reducedCosts[k]; } - void setNegativeObjective(double val) { negativeObjective = val; } - - int findNegativeCost(); - - int findMinStep(int p); - - std::tuple simplexTableau(); - - void pivoting(int p, int q); - - std::pair simplexAlgorithmTableau(); -}; - - -} // namespace dtwc::solver diff --git a/dtwc/solver/SimplexTable.cpp b/dtwc/solver/SimplexTable.cpp deleted file mode 100644 index ce805f8..0000000 --- a/dtwc/solver/SimplexTable.cpp +++ /dev/null @@ -1,201 +0,0 @@ -/* - * SimplexTable.cpp - * - * Sparse implementation of a Simplex table. - - * Created on: 29 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#include "SimplexTable.hpp" -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "EqualityConstraints.hpp" -#include -#include -#include "../timing.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - - -namespace dtwc::solver { - - -int SimplexTable::getRow(int col) const -{ - if (col < 0 || col >= (ntab - 1)) // Check if index is in the valid range - throw std::runtime_error(fmt::format("The index of the variable ({}) must be between 0 and {}", col, ntab - 2)); - - int rowIndex = -1; // Using -1 to represent None - if (!isAround(reducedCosts[col], 0.0)) return rowIndex; - - // So this is checking there is one and only one 1.0 element! -> #TODO change with something keeping book of basic variables in future. - for (auto [key, value] : innerTable[col]) - if (!isAround(value, 0)) // The entry is non zero - { - if (rowIndex == -1 && isAround(value, 1.0)) // The entry is one, and the index has not been found yet. - rowIndex = key; - else - return -1; - } - - return rowIndex; -} - - -int SimplexTable::findNegativeCost() -{ - // Returns -1 if table is optimal otherwise first negative reduced cost. - // Find the first negative cost, if there are none, then table is optimal. - int p = -1; - for (int i = 0; i < reducedCosts.size(); i++) // Reduced cost. - if (reducedCosts[i] < -epsilon) { - p = i; - break; - } - - return p; -} - -int SimplexTable::findMinStep(int p) -{ - constexpr double tol = 1e-10; - // Calculate the maximum step that can be done along the basic direction d[p] - double minStep = std::numeric_limits::infinity(); - int q{ -1 }; // row of min step. -1 means unbounded. - - for (auto [k, minus_d] : innerTable[p]) - if (minus_d > tol) { - const auto step = rhs[k] / minus_d; - if (step < minStep) { - minStep = step; - q = k; - } - } - - return q; -} - -std::tuple SimplexTable::simplexTableau() -{ - // Find the first negative cost, if there are none, then table is optimal. - const int p = findNegativeCost(); - - if (p == -1) // The table is optimal - return std::make_tuple(-1, -1, true, true); - - const int q = findMinStep(p); - - if (q == -1) // The table is unbounded - return std::make_tuple(-1, -1, false, false); - - return std::make_tuple(p, q, false, true); -} - - -void SimplexTable::pivoting(int p, int q) -{ - // p column, q = row. - // Make p,q element one and eliminate all other nonzero elements in that column by basic row operations. - const double thepivot = innerTable[p][q]; - if (isAround(thepivot, 0.0)) - throw std::runtime_error(fmt::format("The pivot is too close to zero: {}", thepivot)); - - - auto oneTask = [this, thepivot, p, q](int i) { - if (i == p) return; // Dont delete the pivot column yet. - - auto currentCol_q = innerTable[i].find(q); - - if (currentCol_q != innerTable[i].end()) // We have a row like that! - { - (currentCol_q->second) /= thepivot; // Normalise by pivot; - - reducedCosts[i] -= reducedCosts[p] * (currentCol_q->second); // Remove from last row. - - for (auto [key, val] : innerTable[p]) // pivot Column - if (key != q) { - auto it_now = innerTable[i].find(key); - if (it_now != innerTable[i].end()) // If that row exists only then subtract otherwise equate. - (it_now->second) -= val * (currentCol_q->second); - else - innerTable[i][key] = -val * (currentCol_q->second); - } - - std::erase_if(innerTable[i], [](const auto &item) { - auto const &[key, value] = item; - return isAround(value, 0.0); // Remove zero elements to make it compact. - }); - } - }; - - dtwc::run(oneTask, innerTable.size()); - - rhs[q] /= thepivot; - // We always have RHS. - for (auto [key, val] : innerTable[p]) - if (key != q) - rhs[key] -= val * rhs[q]; - - negativeObjective -= rhs[q] * reducedCosts[p]; - - // Deal with the pivot column now. - innerTable[p].clear(); - innerTable[p][q] = 1; - reducedCosts[p] = 0; -} - - -std::pair SimplexTable::simplexAlgorithmTableau() -{ - size_t iter{}; - double duration_table{}, duration_pivoting{}; - while (true) { - dtwc::Clock clk; - auto [colPivot, rowPivot, optimal, bounded] = simplexTableau(); - duration_table += clk.duration(); - - if (optimal) return { true, false }; - if (!bounded) return { false, true }; - - clk = dtwc::Clock{}; - - pivoting(colPivot, rowPivot); - duration_pivoting += clk.duration(); - - - if (iter % 100 == 0) { - std::cout << "Iteration " << iter << " is finished!" << std::endl; - std::cout << "Duration table: "; - dtwc::Clock::print_duration(std::cout, duration_table); - std::cout << "Duration pivoting: "; - dtwc::Clock::print_duration(std::cout, duration_pivoting); - - duration_table = 0; - duration_pivoting = 0; - - size_t innerSize = 0; - for (auto &map : innerTable) { - innerSize += map.size(); - } - - std::cout << "Inner size: " << innerSize << '\n'; - } - - iter++; - } -} - -} // namespace dtwc::solver diff --git a/dtwc/solver/SimplexTable.hpp b/dtwc/solver/SimplexTable.hpp deleted file mode 100644 index abc8387..0000000 --- a/dtwc/solver/SimplexTable.hpp +++ /dev/null @@ -1,122 +0,0 @@ -/* - * SimplexTable.hpp - * - * Sparse implementation of a Simplex table. - - * Created on: 22 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#pragma once - -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "EqualityConstraints.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include - - -namespace dtwc::solver { - -class SimplexTable -{ - // Table is mtab x ntab - // Inner table is m x n - int mtab{}, ntab{}; - std::vector> innerTable; // Each is a column - std::vector reducedCosts, rhs; - std::vector rowIndices; - double negativeObjective{}; - -public: - SimplexTable() = default; - SimplexTable(int mtab_, int ntab_) : mtab{ mtab_ }, ntab{ ntab_ }, innerTable(ntab_ - 1), - reducedCosts(ntab_ - 1), rhs(mtab_ - 1), rowIndices(ntab_ - 1) {} - - int rows() const { return mtab; } - int cols() const { return ntab; } - - void clear() - { - rowIndices.clear(); - rhs.clear(); - reducedCosts.clear(); - - for (auto &map : innerTable) // Should not remove allocated memory for maps. - map.clear(); - } - - void removeColumns(int a, int b) - { - // Removes columns [a, b) - innerTable.erase(innerTable.begin() + a, innerTable.begin() + b); - reducedCosts.erase(reducedCosts.begin() + a, reducedCosts.begin() + b); - ntab -= (b - a); - } - - - void createPhaseOneTableau(const EqualityConstraints &eq) - { - const int m = eq.A.rows(), n = eq.A.cols(); - mtab = m + 1; - ntab = m + n + 1; - clear(); // Clear the things. - innerTable.resize(m + n); // Adding m auxillary variables. - rhs = eq.b; - reducedCosts.resize(n + m, 0.0); // Set columns n through n+m of the last row to 0.0 - - for (const auto [key, val] : eq.A.data) { - const auto [i, j] = key; - innerTable[j][i] = val; - reducedCosts[j] -= val; // Set the first n columns of the last row - } - - for (int i = 0; i < m; i++) - innerTable[n + i][i] = 1.0; - - negativeObjective = -std::reduce(rhs.begin(), rhs.end()); - } - - int getRow(int col) const; - - double getObjective() const { return -negativeObjective; } - - double getValue(int k) const - { - const auto basicRowNo = getRow(k); - return (basicRowNo != -1) ? rhs[basicRowNo] : 0.0; - } - - double getRHS(int k) const { return rhs[k]; } - - double getReducedCost(int k) const { return reducedCosts[k]; } - - double inner(int i, int j) { return innerTable[j][i]; } - - - double &setReducedCost(int k) { return reducedCosts[k]; } - void setNegativeObjective(double val) { negativeObjective = val; } - - int findNegativeCost(); - - int findMinStep(int p); - - std::tuple simplexTableau(); - - void pivoting(int p, int q); - - std::pair simplexAlgorithmTableau(); -}; - - -} // namespace dtwc::solver diff --git a/dtwc/solver/SparseMatrix.hpp b/dtwc/solver/SparseMatrix.hpp deleted file mode 100644 index 7badc8b..0000000 --- a/dtwc/solver/SparseMatrix.hpp +++ /dev/null @@ -1,72 +0,0 @@ -/* - * SparseMatrix.hpp - * - * A map based implementation of a Sparse Matrix. - - * Created on: 22 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - -#pragma once - -#include "sparse_util.hpp" - -#include -#include -#include -#include -#include - -namespace dtwc::solver { - -template -struct SparseMatrix -{ - std::map data; - int m{}, n{}; // rows and columns - - SparseMatrix() = default; - SparseMatrix(int m_, int n_) : m{ m_ }, n{ n_ } {} - - double operator()(int x, int y) const - { - auto it = data.find(Coordinate{ x, y }); - return (it != data.end()) ? (it->second) : 0.0; - } - - double &operator()(int x, int y) { return data[Coordinate{ x, y }]; } - - void compress() - { - std::erase_if(data, [](const auto &item) { - auto const &[key, value] = item; - return isAround(value, 0.0); - }); - } - - void print() const - { - for (int i = 0; i < m; i++) { - for (int j = 0; j < n; j++) { - if (j != 0) std::cout << ','; - std::cout << std::setw(3) << (*this)(i, j); - } - std::cout << '\n'; - } - } - - int rows() const { return m; } - int cols() const { return n; } - - void expand(int m_, int n_) - { - m = m_; - n = n_; - } - - auto row_begin(int i) { return data.lower_bound({ i, 0 }); } - auto row_end(int i) { return data.lower_bound({ i + 1, 0 }); } - - // std::span> row_range(int i) { return std::span>{ row_begin(i), row_end(i) }; } -}; -} // namespace dtwc::solver \ No newline at end of file diff --git a/dtwc/solver/SparseSimplex.cpp b/dtwc/solver/SparseSimplex.cpp deleted file mode 100644 index 23a8623..0000000 --- a/dtwc/solver/SparseSimplex.cpp +++ /dev/null @@ -1,255 +0,0 @@ -/* - * SparseSimplex.cpp - * - * LP solution - - * Created on: 26 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#include "SparseSimplex.hpp" -#include "SimplexTable.hpp" - -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "../Problem.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - - -namespace dtwc::solver { - -void SparseSimplex::warmStartPhaseOne() -{ - if (Nb == 0 || Nc == 0) return; // it is not defined by clustering problem so we dont know what to do. - dtwc::Clock clk{}; - std::cout << "Warmstart starting!" << std::endl; - // Just activate first Nc variables as clusters and fill all to the first cluster. - for (int i = 0; i < Nb; i++) { - const int p = (i < Nc) ? (i * (Nb + 1)) : i; - const int q = table.findMinStep(p); - if (q != -1) - table.pivoting(p, q); - } - - // const int slack_begin = (Nb * Nb) + (Nb - 1); - // const int slack_end = (Nb * Nb) + (Nb - 1) * Nb; - - // for (int p = slack_begin; p < slack_end; p++) { // Make slack variables one. - // const int q = table.findMinStep(p); - // if (q != -1) - // table.pivoting(p, q); - - // if (p % 100 == 0) - // std::cout << "Pivoting variable " << p << " of " << slack_end << std::endl; - // } - - std::cout << "Warmstart ended in " << clk << std::endl; -}; - - -std::tuple SparseSimplex::simplex() -{ - // bool unbounded{}, optimal{}; - const int m = eq.A.rows(), n = eq.A.cols(); - - // Ensure all elements of b are non-negative - for (int i = 0; i < m; ++i) - if (eq.b[i] < 0) { - eq.b[i] = -eq.b[i]; - for (auto it = eq.A.row_begin(i), it_end = eq.A.row_end(i); it != it_end; ++it) - it->second = -it->second; - } - std::cout << "Creating Phase-I table." << std::endl; - table.createPhaseOneTableau(eq); - // warmStartPhaseOne(); - - std::cout << "Running algorithm with Phase-I table." << std::endl; - auto [optimal, unbounded] = table.simplexAlgorithmTableau(); - - if (unbounded) return { true, false }; - if (table.getObjective() > epsilon) return { false, true }; // Infeasible problem - - std::vector basicRows(m + n); - std::set basicIndices; - std::set tobeCleaned; - - for (int k = 0; k < m + n; ++k) { - basicRows[k] = table.getRow(k); - if (basicRows[k] != -1) { - basicIndices.insert(k); - if (k >= n) tobeCleaned.insert(k); // If k>= n are basic indices then clean them. - } - } - - if (!tobeCleaned.empty()) { - std::cout << "There are things to be cleaned!" << std::endl; - throw "There are things to be cleaned!\n"; - } - - // while (!tobeCleaned.empty()) { - // int auxiliaryColumn = *tobeCleaned.begin(); - // tobeCleaned.erase(tobeCleaned.begin()); - - // int rowpivotIndex = basicRows[auxiliaryColumn]; - // VectorXd rowpivot = phaseOneTableau.row(rowpivotIndex); - - // std::vector originalNonbasic; - // for (int i = 0; i < n; ++i) - // if (basicIndices.find(i) == basicIndices.end()) - // originalNonbasic.push_back(i); - - // int colpivot = -1; - // double maxVal = std::numeric_limits::epsilon(); - - // for (int col : originalNonbasic) - // if (std::abs(rowpivot(col)) > maxVal) { - // maxVal = std::abs(rowpivot(col)); - // colpivot = col; - // } - - // if (colpivot != -1) { - // pivoting(phaseOneTableau, colpivot, rowpivotIndex); - // basicRows[colpivot] = rowpivotIndex; - // basicRows[auxiliaryColumn] = -1; - // } else { - // phaseOneTableau.conservativeResize(phaseOneTableau.rows() - 1, phaseOneTableau.cols()); - // for (int k = 0; k < m + n; ++k) - // basicRows[k] = phaseOneTableau.getRow(k); - // } - // } - - std::cout << "Creating Phase-II table." << std::endl; - - table.removeColumns(n, n + m); - - basicRows.resize(n); - for (int k = 0; k < n; ++k) - basicRows[k] = table.getRow(k); - - // Calculate last row - std::vector basicIndicesList, nonbasicIndicesList; - for (int i = 0; i < n; ++i) - if (basicRows[i] != -1) - basicIndicesList.push_back(i); - else - nonbasicIndicesList.push_back(i); - - for (int k : nonbasicIndicesList) { - double sumVal = 0.0; - for (int j : basicIndicesList) { - if (j >= c.size()) break; // Because after this is zero. - sumVal += c[j] * table.inner(basicRows[j], k); - } - const auto c_k = (k < c.size()) ? c[k] : 0.0; - table.setReducedCost(k) = c_k - sumVal; - } - - double lastRowSum = 0.0; - for (int j : basicIndicesList) { - if (j >= c.size()) break; // Because after this is zero. - lastRowSum += c[j] * table.getRHS(basicRows[j]); - } - - table.setNegativeObjective(-lastRowSum); - - // Phase II - std::cout << "Running algorithm with Phase-II table." << std::endl; - auto [optimal2, unbounded2] = table.simplexAlgorithmTableau(); // it was startPhaseTwo - return { unbounded2, false }; -} - -std::pair, double> SparseSimplex::getResults() const -{ - int n = table.cols() - 1; // finalTableau - - std::vector solution(n); - for (int k = 0; k < n; ++k) - solution[k] = table.getValue(k); - - return { solution, table.getObjective() }; // Solution and optimal cost -} - - -void SparseSimplex::gomory() -{ - fmt::println("=================================="); - fmt::println("Problem with {} variables and {} constraints", eq.A.cols(), eq.A.rows()); - - bool unbounded{}, infeasible{}; - std::tie(unbounded, infeasible) = simplex(); // Assuming simplex is defined - - if (unbounded) { - fmt::println("Unbounded problem"); - nGomory = 0; - return; - } - - if (infeasible) { - fmt::println("Infeasible problem"); - nGomory = 0; - return; - } - - auto [solution, copt] = getResults(); // Assuming getResults is defined - if constexpr (settings::debug_Simplex) - fmt::println("Solution: {} and Copt = [{}]\n", solution, copt); - - - int m_now = eq.A.rows(); - int n_now = eq.A.cols(); - - nGomory = 0; - for (int i = 0; i < (table.rows() - 1); i++) { - double rhs_now = table.getRHS(i); - if (isFractional(rhs_now)) // Scan last column - { - eq.A(m_now + nGomory, n_now + nGomory) = 1.0; // Slack variable for new row. - - for (auto j = 0; j < eq.A.cols(); j++) { - auto val = std::floor(table.inner(i, j)); // gamma - - if (!isAround(val)) - eq.A(m_now + nGomory, j) = val; - } - - eq.b.emplace_back(std::floor(rhs_now)); - nGomory++; // One more gomory cut. - } - } - - eq.A.expand(m_now + nGomory, n_now + nGomory); - - c.resize(n_now + nGomory); // add zeros - fmt::println("Number of Gomory cuts: {}", nGomory); -} - -SparseSimplex::SparseSimplex(Problem &prob) -{ - Nb = prob.data.size(); - Nc = prob.cluster_size(); - - eq = defaultConstraints(Nb, Nc); - - c.resize(Nb * Nb); - - for (size_t j{ 0 }; j < Nb; j++) - for (size_t i{ 0 }; i < Nb; i++) - c[i + j * Nb] = prob.distByInd_scaled(i, j); -} - - -} // namespace dtwc::solver diff --git a/dtwc/solver/SparseSimplex.hpp b/dtwc/solver/SparseSimplex.hpp deleted file mode 100644 index 3d36d96..0000000 --- a/dtwc/solver/SparseSimplex.hpp +++ /dev/null @@ -1,68 +0,0 @@ -/* - * SparseSimplex.hpp - * - * LP solution - - * Created on: 26 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - - -#pragma once - -#include "../settings.hpp" -#include "../utility.hpp" -#include "solver_util.hpp" -#include "SimplexTable.hpp" -#include "SimplexRowTable.hpp" -#include "SimplexFlatRowTable.hpp" -#include "SimplexFlatTable.hpp" - -#include -#include -#include -#include -#include - - -#include -#include -#include -#include -#include -#include - - -namespace dtwc { -class Problem; -} - -namespace dtwc::solver { - -class SparseSimplex -{ - using VectorType = std::vector; - // -1 means None. - int nGomory{ -1 }; - int Nb{}, Nc{}; - SimplexFlatTable table; - EqualityConstraints eq; - VectorType c; // c*x = cost. - -public: - SparseSimplex() = default; - SparseSimplex(Problem &prob); - SparseSimplex(EqualityConstraints &eq_, const VectorType &c_) : eq(eq_), c(c_) {} - - void gomory(); - void gomoryAlgorithm() - { - while (nGomory != 0) gomory(); - } - - std::pair, double> getResults() const; - std::tuple simplex(); - void warmStartPhaseOne(); // Warm-starts Phase-I as we know a feasible solution. -}; - -} // namespace dtwc::solver diff --git a/dtwc/solver/test.hpp b/dtwc/solver/test.hpp deleted file mode 100644 index 61b58cd..0000000 --- a/dtwc/solver/test.hpp +++ /dev/null @@ -1,36 +0,0 @@ -/* - * test.hpp - * - * Test problems - - * Created on: 29 Oct 2023 - * Author(s): Volkan Kumtepeli, Becky Perriment - */ - -#pragma once - -#include "EqualityConstraints.hpp" -#include -#include - -namespace dtwc::solver { -auto get_prob_small() -{ - EqualityConstraints eq(2, 4); - std::vector c{ 1, -2, 0, 0 }; - - eq.A(0, 0) = -4; - eq.A(0, 1) = 6; - eq.A(0, 2) = 1; - - eq.A(1, 0) = 1; - eq.A(1, 1) = 1; - eq.A(1, 3) = 1; - - eq.b[0] = 5; - eq.b[1] = 5; - - return std::pair(eq, c); -} - -} // namespace dtwc::solver \ No newline at end of file diff --git a/dtwc/solver/sparse_util.hpp b/dtwc/types/element_types.hpp similarity index 96% rename from dtwc/solver/sparse_util.hpp rename to dtwc/types/element_types.hpp index 4d85078..1975d87 100644 --- a/dtwc/solver/sparse_util.hpp +++ b/dtwc/types/element_types.hpp @@ -1,5 +1,5 @@ /* - * sparse_util.hpp + * element_types.hpp * * Helper class for sparse matrix utilities. @@ -10,7 +10,7 @@ #pragma once -#include "solver_util.hpp" +#include "types_util.hpp" namespace dtwc::solver { @@ -23,6 +23,20 @@ struct Element Element(int index_, double value_) : index(index_), value(value_) {} }; +struct Coordinate +{ + int row{}, col{}; // Row and column of the value +}; + +struct Triplet +{ + int row{}, col{}; // Row and column of the value + double val{}; + + Triplet() = default; + Triplet(int row_, int col_, double val_) : row(row_), col(col_), val(val_) {} +}; + struct CompElementIndices { bool operator()(const Element &c1, const Element &c2) const @@ -43,21 +57,6 @@ struct CompElementValuesAndIndices } }; - -struct Coordinate -{ - int row{}, col{}; // Row and column of the value -}; - -struct Triplet -{ - int row{}, col{}; // Row and column of the value - double val{}; - - Triplet() = default; - Triplet(int row_, int col_, double val_) : row(row_), col(col_), val(val_) {} -}; - struct RowMajor { bool operator()(const auto &c1, const auto &c2) const @@ -74,5 +73,4 @@ struct ColumnMajor } }; - } // namespace dtwc::solver \ No newline at end of file diff --git a/dtwc/types/types.hpp b/dtwc/types/types.hpp index 896fe20..02abc51 100644 --- a/dtwc/types/types.hpp +++ b/dtwc/types/types.hpp @@ -4,3 +4,4 @@ #include "BandMatrix.hpp" #include "SkewedBandMatrix.hpp" #include "Range.hpp" +#include "element_types.hpp" diff --git a/dtwc/solver/solver_util.hpp b/dtwc/types/types_util.hpp similarity index 70% rename from dtwc/solver/solver_util.hpp rename to dtwc/types/types_util.hpp index 3c07bf1..0764aa5 100644 --- a/dtwc/solver/solver_util.hpp +++ b/dtwc/types/types_util.hpp @@ -1,7 +1,7 @@ /* - * solver_util.hpp + * types_util.hpp * - * Utility functions for solver + * Utility functions for types * Created on: 19 Dec 2022 * Author(s): Volkan Kumtepeli, Becky Perriment @@ -17,7 +17,6 @@ namespace dtwc::solver { - constexpr data_t int_threshold = 0.01; constexpr double epsilon = 1e-8; @@ -27,13 +26,6 @@ bool inline isAround(double x, double y = 0.0, double tolerance = epsilon) } bool inline isFractional(double x) { return std::abs(x - std::round(x)) > epsilon; } -enum class ConvergenceFlag { - error_sizeNotSet = -1, // inline bool is_one(T x) { return x > (1 - int_threshold); } diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 72667cc..b02ffde 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -1,2 +1,3 @@ # Add executables and tests with coverage flags applied -add_executable_with_coverage_and_test(unit_test_warping unit_test_warping.cpp) +add_executable_with_coverage_and_test(unit_test_warping) +add_executable_with_coverage_and_test(unit_test_Problem) diff --git a/tests/unit/unit_test_Problem.cpp b/tests/unit/unit_test_Problem.cpp new file mode 100644 index 0000000..c9d7cf0 --- /dev/null +++ b/tests/unit/unit_test_Problem.cpp @@ -0,0 +1,88 @@ +/* + * unit_test_warping.cpp + * + * Unit test file for time warping functions + * Created on: 03 Dec 2023 + * Author(s): Volkan Kumtepeli, Becky Perriment + */ + +#include + +#include +#include +#include + +using Catch::Matchers::WithinAbs; + +using namespace dtwc; + +TEST_CASE("dtwFull_test", "[dtwFull]") +{ + dtwc::Clock clk; // Create a clock object + std::string probName = "DTW_kMeans_results"; + + auto Nc = 3; // Number of clusters + + constexpr int N_repetition = 5; + constexpr int maxIter = 100; + constexpr int Ndata_max = 10; + + dtwc::DataLoader dl{ settings::dataPath / "dummy", Ndata_max }; + dl.startColumn(1).startRow(1); // Since dummy files are in Pandas format skip first row/column. + + dtwc::Problem prob{ probName, dl }; // Create a problem. + + prob.set_numberOfClusters(Nc); // Nc = number of clusters. + + + REQUIRE(prob.cluster_size() == Nc); + REQUIRE(prob.name == probName); + + // prob.cluster_by_kMedoidsPAM_repetetive(N_repetition, maxIter); +} + +TEST_CASE("dtwFull_L_test", "[dtwFull_L]") +{ + using data_t = double; + std::vector x{ 1, 2, 3 }, y{ 3, 4, 5, 6, 7 }, z{ 1, 2, 3 }, empty{}; + constexpr double ground_truth = 13; + + // Zero distance between same vectors: + REQUIRE_THAT(dtwFull_L(x, x), WithinAbs(0, 1e-15)); + REQUIRE_THAT(dtwFull_L(x, z), WithinAbs(0, 1e-15)); + REQUIRE_THAT(dtwFull_L(z, x), WithinAbs(0, 1e-15)); + + // Some distance between others: 13 + REQUIRE_THAT(dtwFull_L(x, y), WithinAbs(ground_truth, 1e-15)); + REQUIRE_THAT(dtwFull_L(y, x), WithinAbs(ground_truth, 1e-15)); + + // Empty vector should give infinite cost. + REQUIRE(dtwFull_L(x, empty) > 1e10); + REQUIRE(dtwFull_L(empty, x) > 1e10); +} + +TEST_CASE("dtwBanded_test", "[dtwBanded]") +{ + using data_t = double; + std::vector x{ 1, 2, 3 }, y{ 3, 4, 5, 6, 7 }, z{ 1, 2, 3 }, empty{}; + constexpr double ground_truth = 13; + + // Zero distance between same vectors: + REQUIRE_THAT(dtwBanded(x, x), WithinAbs(0, 1e-15)); + REQUIRE_THAT(dtwBanded(x, z), WithinAbs(0, 1e-15)); + REQUIRE_THAT(dtwBanded(z, x), WithinAbs(0, 1e-15)); + + // Some distance between others with too large band, should be same as unbanded. + int band = 100; + REQUIRE_THAT(dtwBanded(x, y, band), WithinAbs(ground_truth, 1e-15)); + REQUIRE_THAT(dtwBanded(y, x, band), WithinAbs(ground_truth, 1e-15)); + + // Banded distance: + band = 2; + REQUIRE_THAT(dtwBanded(x, y, band), WithinAbs(ground_truth, 1e-15)); + REQUIRE_THAT(dtwBanded(y, x, band), WithinAbs(ground_truth, 1e-15)); + + // Empty vector should give infinite cost. + REQUIRE(dtwBanded(x, empty) > 1e10); + REQUIRE(dtwBanded(empty, x) > 1e10); +} \ No newline at end of file