diff --git a/CMakeLists.txt b/CMakeLists.txt index 03513d4..43369e7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.30) +cmake_minimum_required(VERSION 3.28) project(UDBM VERSION 2.0.15 LANGUAGES CXX C) include(CMakePackageConfigHelpers) include(GNUInstallDirs) diff --git a/cmake/boost.cmake b/cmake/boost.cmake index ce7d4df..f6dc1b4 100644 --- a/cmake/boost.cmake +++ b/cmake/boost.cmake @@ -9,6 +9,9 @@ set(Boost_USE_STATIC_RUNTIME ON) # Mac insists on ON for boost_program_options #set(Boost_DEBUG ON) set(Boost_VERSION 1.86.0) +if (POLICY CMP0167) + cmake_policy (SET CMP0167 NEW) +endif () if (BOOST_INCLUDE_LIBRARIES) find_package(Boost ${Boost_VERSION} COMPONENTS ${BOOST_INCLUDE_LIBRARIES} QUIET) else(BOOST_INCLUDE_LIBRARIES) diff --git a/cmake/doctest.cmake b/cmake/doctest.cmake index acbd69a..39196ed 100644 --- a/cmake/doctest.cmake +++ b/cmake/doctest.cmake @@ -1,4 +1,4 @@ -find_package(doctest 2.4.11 QUIET) +find_package(doctest 2.5.2 QUIET) if (doctest_FOUND) if (TARGET doctest::doctest_with_main) @@ -25,7 +25,7 @@ else(doctest_FOUND) FetchContent_Declare( doctest GIT_REPOSITORY https://github.com/doctest/doctest - GIT_TAG v2.4.11 + GIT_TAG v2.5.2 GIT_SHALLOW TRUE # get only the last commit version GIT_PROGRESS TRUE # show progress of download # FIND_PACKAGE_ARGS NAMES doctest diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h new file mode 100644 index 0000000..56fea82 --- /dev/null +++ b/include/dbm/ParamPricedDBM.h @@ -0,0 +1,355 @@ +#ifndef INCLUDE_DBM_PARAMPRICEDDBM_H +#define INCLUDE_DBM_PARAMPRICEDDBM_H + +#include "dbm.h" + +#include +#include +#include +#include +#include +#include + + + +namespace PPDBM +{ + ///to facilitate notations + using Vertex = std::vector; // size _dim, clock valuations + using CostVector = std::vector; // size _param + 1, coeffs on parameters + ///copied from relation.h because include problems otherwise + using relation_t = enum { /* EXACT relation | NON EXACT relation */ + /*--------------------|--------------------*/ + base_DIFFERENT = 0, /**< incomparable | not (set1 <= set2) */ + base_SUPERSET = 1, /**< set1 > set2 | not used */ + base_GREATER = 1, /**< same as superset | */ + base_SUBSET = 2, /**< set1 < set2 | set1 <= set2 */ + base_LESS = 2, /**< same as subset | */ + base_EQUAL = 3 /**< set1 == set2 | not used */ + }; + ///to represent infinity + constexpr auto INF = std::numeric_limits::max() >> 1u; + + + /// generic matrix type (used everywhere) + template + struct Matrix + { + Matrix(uint32_t rows, uint32_t cols, const T x): _rows{rows}, _cols{cols}, _data(rows * cols, x){} + auto rows() const { return _rows; } + auto cols() const { return _cols; } + const std::vector& getData() const { return _data; } + std::vector& getData() { return _data; } + T& operator()(int i, int j) { return _data[i * _cols + j]; } + const T& operator()(int i, int j) const { return _data[i * _cols + j]; } + auto fill(T x){ std::fill(_data.begin(), _data.end(), x); } + + private: + uint32_t _rows; + uint32_t _cols; + std::vector _data; + }; + + + + + /** + * A parameter constraint is of the form (a_0 + sum^n_{i=1} a_i p_i) <= 0 + */ + struct ParametricConstraint + { + constexpr ParametricConstraint(std::vector coeffs, bool strict): + _coeffs{std::move(coeffs)}, + _strict{strict} {} + constexpr const bool& getStrict() const { return _strict; } + constexpr const std::vector& getCoeffs() const { return _coeffs; } + private: + std::vector _coeffs; + bool _strict; + }; + + struct Polyhedron + { + constexpr std::vector& getConstraints(){ return _constraints; } + constexpr const std::vector& getConstraints() const { return _constraints; } + constexpr void addParametricConstraint(const ParametricConstraint& constraint) + { _constraints.push_back(constraint); } + + bool isEmpty(); + /// removes every constraint implied by the others (does not change the represented set) + void minimize(); + static void printPolyhedron(const std::vector& constraints); + private: + std::vector _constraints; + }; + + + + + /** + * a constraint stored in a DBM is of the form x_i - x_j ~ value with ~ being < or <= + */ + struct DbmBound + { + constexpr auto getValue() const { return _value; } + constexpr auto getStrict() const { return _strict; } + + constexpr DbmBound( int32_t value, bool strict): _value{value}, _strict{strict}{} + + /// is this constraint unrestraining? + constexpr bool isInfinite() const { return _value == INF; } + + /// is the new constraint useless? + constexpr bool dominates(const DbmBound& other) const { + return _value < other._value || (_value == other._value && (_strict || !other._strict)); } + + /// do the two constraints force the zone to be empty? + static constexpr bool isIncompatible(const DbmBound& a, const DbmBound& b) { + return a.getValue() + b.getValue() < 0 || + (a.getValue() + b.getValue() == 0 && (a.getStrict() || b.getStrict())); + } + + private: + int32_t _value : 31; + bool _strict : 1; /// constant coeff of the constraint, is the constraint strict? + }; + + constexpr auto INF_BOUND = DbmBound{INF,true}; + constexpr auto ZERO_BOUND = DbmBound{0,false}; + + /// x_i - x_j <= a and x_j - x_k < b gives x_i - x_k < a + b + constexpr DbmBound operator+(const DbmBound& a, const DbmBound& b) + { + if (a.isInfinite() || b.isInfinite()) return INF_BOUND; + return DbmBound{ a.getValue() + b.getValue(), a.getStrict() || b.getStrict() }; + } + + + /** + * + */ + struct DbmConstraint + { + uint32_t i; + uint32_t j; + DbmBound bound; + }; + + + + + /** + * + */ + struct DBMatrix + { + constexpr DBMatrix(uint32_t dim) : _matrix{dim, dim, INF_BOUND} + { + for (uint32_t i = 0; i < dim; ++i) { _matrix(i, i) = ZERO_BOUND; _matrix(0,i) = ZERO_BOUND; } + } + + constexpr auto getDim() const { return _matrix.rows(); } + constexpr DbmBound& operator()(int i, int j) { return _matrix(i, j); } + constexpr const DbmBound& operator()(int i, int j) const { return _matrix(i, j); } + constexpr void reset() + { + _matrix.fill(INF_BOUND); + for (uint32_t i = 0; i < _matrix.rows(); ++i) { + _matrix(i,i)= ZERO_BOUND; + _matrix(0,i) = ZERO_BOUND; + } + } + void close(uint32_t i, uint32_t j); + void relaxDown(); + void relaxUp(); + void relaxDownClock(uint32_t clock); + void relaxUpClock(uint32_t clock); + std::vector findLinkedClocks() const; + + relation_t relation(const DBMatrix& other) const; ///< pure zone comparison (no cost involved). + std::vector verticesDBM() const; + + + private: + Matrix _matrix; + }; + + + + + + struct Facet; + struct Edge; + struct ParamPricedTimedAutomata; + + + + + //////////////////////////// The PPDBM class that represents Parametric priced zones ////////////////////////////// + /** + * Parametric priced timed difference bound matrix + */ + struct Ppdbm + { + Ppdbm(uint32_t dim, uint32_t param, bool init = false): + _dim{dim}, + _param{param}, + _dbm{dim}, + _offsetCost(param + 1), + _rates(dim, param + 1, 0) + { + if (init) + for (uint32_t i = 0; i < _dim; ++i) {_dbm(i,0) = ZERO_BOUND;} + } + + constexpr uint32_t getDim() const { return _dim; } + constexpr uint32_t getParam() const { return _param; } + constexpr Polyhedron& getPC() {return _PC;} + constexpr const Polyhedron& getPC() const {return _PC;} + constexpr Matrix& getRates() { return _rates; } + constexpr const Matrix& getRates() const { return _rates; } + constexpr std::vector& getOffsetCost() { return _offsetCost; } + constexpr const std::vector& getOffsetCost() const { return _offsetCost; } + constexpr DbmBound& operator()(int i, int j) { return _dbm(i, j); } + constexpr const DbmBound& operator()(int i, int j) const { return _dbm(i, j); } + void display() const; + + constexpr bool isEmpty() const { return _emptyZone; } + void reset(); + bool constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint); + void removeCostAtOffset(); + void addCostAtOffset(); + bool isDiagonalNegative(uint32_t i) const; + bool constrain_DBM_N(const std::vector& constraints); + std::vector costAtOtherOffset(const DBMatrix& otherDbm) const; + CostVector getSlope() const; + void addConstantCost(const std::vector& q); + static bool alreadyCoveredBy(std::pair symbolicState, + const std::vector>& symbolicStateListe); + + /// Delay operation on a facet put in the form Z \land (clock - otherClock = value). + void delay(uint32_t clock, const CostVector& c); + void delay(); + /// clock reset operation on a facet put in the form Z \land (clockToReset - otherClock = value). + void clockReset(uint32_t clockToReset, DbmConstraint facetConstraint); + /// lower & upper facets functions for the operators + std::vector lowerFacets() const; + std::vector upperFacets() const; + std::vector lowerFacetsRelativeTo(uint32_t clock) const; + std::vector upperFacetsRelativeTo(uint32_t clock) const; + + /// operator post, computes the successors of the ppdbm by delay. + std::vector post_delta(CostVector p, const std::vector& J) const; + /// operator post, computes the succesors of the ppdbm by transition e. + std::vector post_e(const std::vector& g, uint32_t resetClock, const std::vector& q) const; + /// + static std::vector> post(const std::pair& symbolicState, ParamPricedTimedAutomata automaton); + + /// Is there a domination relation between the two zones? + relation_t relation(const Ppdbm& other) const; + /// Vertices of the zone. + constexpr std::vector verticesDBM() const{ return _dbm.verticesDBM(); }; + /// Evaluates the cost vector(coeffs on parameters) of *this* at vertex `v`. + CostVector costAtVertex(const Vertex& v) const; + /// a dominates b coefficient by coefficient? + static bool dominatesVector(const CostVector& a, const CostVector& b); + static bool equalVector(const CostVector& a, const CostVector& b); + /// Compare two lists of cost vectors aligned vertex by vertex. + static relation_t compareCostLists(const std::vector& c1, const std::vector& c2); + + /// removes every result whose (cost, region) is subsumed by another result with the SAME cost function + static std::vector pruneRedundant(std::vector results); + + private: + uint32_t _dim, _param = 0; // number of clocks and parameters + DBMatrix _dbm; // the classic zone + bool _emptyZone = false; // is the zone represented by the dbm empty? + std::vector _offsetCost; // the cost of the offset (affine function of parameters with int coeffs) + Matrix _rates; // the cost rates of clocks (affine functions of parameters with int coeffs) + Polyhedron _PC; // parametric constraints set + }; + + + struct Facet + { + Ppdbm _ppdbm; + DbmConstraint _constraint; + }; + + struct OptimalRegion + { + static std::vector pruneRedundantRegions(std::vector regions); + + CostVector cost; + std::vector constraints; + }; + + /// splits a set of (cost, region) results into the regions where each cost is + /// actually the minimum among all of them (the "lower envelope"). + std::vector resolveOptimalPartition(const std::vector& results); + + + + + //////////////////////////////////////////////// automata model /////////////////////////////////////////////////// + + struct Edge + { + Edge(uint32_t from, uint32_t to, uint32_t clockToReset, const std::vector& guards, + std::vector costFunction): + _from{from}, + _to{to}, + _clockToReset{clockToReset}, + _guards{guards}, + _costFunction{std::move(costFunction)} {} + + uint32_t _from; + uint32_t _to; + uint32_t _clockToReset; + std::vector _guards; + std::vector _costFunction; + }; + + struct ParamPricedTimedAutomata + { + ParamPricedTimedAutomata(uint32_t nbLocations, uint32_t initLocation, uint32_t nbClocks, uint32_t nbParam, + std::vector edges, PPDBM::Matrix inv, PPDBM::Matrix costs): + _nbLocations{nbLocations}, + _initialLocation{initLocation}, + _nbClocks{nbClocks}, + _nbParam{nbParam}, + _edges{std::move(edges)}, + _invariants{std::move(inv)}, + _locationCosts{std::move(costs)} {} + + std::vector getCostFunction(uint32_t location) + { + std::vector costVector{}; + for (uint32_t j = 0; j <= _nbParam; ++j) { + costVector.push_back(_locationCosts(location,j)); + } + return costVector; + }; + + std::vector getInvariants(uint32_t location) + { + std::vector invariantsVector{}; + for (uint32_t j = 0; j <= _nbClocks; ++j) { + if (_invariants(location,j).getValue() != PPDBM::INF) invariantsVector.push_back(PPDBM::DbmConstraint{j,0, _invariants(location,j)}); + } + return invariantsVector; + } + + uint32_t _nbLocations; + uint32_t _initialLocation; + uint32_t _nbClocks; + uint32_t _nbParam; + std::vector _edges; + PPDBM::Matrix _invariants; + PPDBM::Matrix _locationCosts; + }; + + +} // namespace PPDBM + +#endif \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8d0be52..430e5e6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,6 +1,6 @@ add_library(UDBM STATIC DBMAllocator.cpp dbm.c fed_dbm.cpp mingraph.c mingraph_read.c partition.cpp print.cpp gen.c mingraph_cache.cpp mingraph_relation.c pfed.cpp fed.cpp infimum.cpp mingraph_equal.c mingraph_write.c - priced.cpp valuation.cpp) + priced.cpp valuation.cpp ParamPricedDBM.cpp) set_property(TARGET UDBM PROPERTY C_VISIBILITY_PRESET hidden) set_property(TARGET UDBM PROPERTY VISIBILITY_INLINES_HIDDEN ON) if (NOT CMAKE_SYSTEM_NAME STREQUAL Windows) # unknown argument: '-fno-keep-inline-dllexport' @@ -18,3 +18,6 @@ target_include_directories(UDBM # where external projects will look for the library's public headers $ ) + +add_executable(param_reachability ParametricOptimalReachability.cpp) +target_link_libraries(param_reachability PRIVATE UDBM glpk) \ No newline at end of file diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp new file mode 100644 index 0000000..4ed05eb --- /dev/null +++ b/src/ParamPricedDBM.cpp @@ -0,0 +1,1345 @@ +/* -*- mode: C++; c-file-style: "stroustrup"; c-basic-offset: 4; indent-tabs-mode: nil; -*- */ + +#include "dbm/ParamPricedDBM.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace PPDBM +{ + /** + * + */ + void Ppdbm::display() const + { + std::cout << "==================== Ppdbm ====================\n"; + std::cout << "dim = " << _dim << " | param = " << _param + << " | empty = " << (_emptyZone ? "true" : "false") << "\n\n"; + + // ---- DBM ---- + std::cout << "-- DBM (x_i - x_j ~ value) --\n"; + std::cout << std::setw(6) << " "; + for (uint32_t j = 0; j < _dim; ++j) + std::cout << std::setw(10) << ("j=" + std::to_string(j)); + std::cout << "\n"; + + for (uint32_t i = 0; i < _dim; ++i) { + std::cout << std::setw(6) << ("i=" + std::to_string(i)); + for (uint32_t j = 0; j < _dim; ++j) { + const auto& b = _dbm(i, j); + std::ostringstream cell; + if (b.isInfinite()) { + cell << "inf"; + } else { + cell << (b.getStrict() ? "<" : "<=") << b.getValue(); + } + std::cout << std::setw(10) << cell.str(); + } + std::cout << "\n"; + } + + // ---- Rates ---- + std::cout << "\n-- Rates (cost rate per clock, coeffs on parameters) --\n"; + for (uint32_t i = 0; i < _dim; ++i) { + std::cout << "clock " << i << ": [ "; + for (uint32_t k = 0; k <= _param; ++k) { + std::cout << _rates(i, k); + if (k != _param) std::cout << ", "; + } + std::cout << " ]\n"; + } + + // ---- Offset cost ---- + std::cout << "\n-- Offset cost (coeffs on parameters) --\n[ "; + for (uint32_t k = 0; k < _offsetCost.size(); ++k) { + std::cout << _offsetCost[k]; + if (k + 1 != _offsetCost.size()) std::cout << ", "; + } + std::cout << " ]\n"; + + // ---- Parametric constraints (PC) ---- + std::cout << "\n-- Polyhedron PC (" << _PC.getConstraints().size() << " constraints) --\n"; + for (const auto& c : _PC.getConstraints()) { + const auto& coeffs = c.getCoeffs(); + std::cout << " "; + bool first = true; + for (size_t k = 0; k < coeffs.size(); ++k) { + if (coeffs[k] == 0) continue; + if (!first) std::cout << " + "; + if (k == 0) + std::cout << coeffs[k]; + else + std::cout << coeffs[k] << "*p" << k; + first = false; + } + if (first) std::cout << "0"; // toutes les coeffs sont nulles + std::cout << (c.getStrict() ? " < 0" : " <= 0") << "\n"; + } + std::cout << "================================================\n"; + } + + /** + * determines if a polyhedron is empty using glpk, generated by claude. Dont know how to do it otherwise :/ + * @return + */ + bool Polyhedron::isEmpty() + { + if (_constraints.empty()) + return false; // aucune contrainte => tout N^param, non vide + + // coeffs = [a_0, a_1, ..., a_n] => nbParam = taille - 1 + const int nbParam = static_cast(_constraints[0].getCoeffs().size()) - 1; + + glp_prob* lp = glp_create_prob(); + glp_set_obj_dir(lp, GLP_MIN); // pas d'objectif réel, on ne teste que la faisabilité + + // colonnes = p_1..p_n, entiers, >= 0 + glp_add_cols(lp, nbParam); + for (int j = 1; j <= nbParam; ++j) + { + glp_set_col_kind(lp, j, GLP_IV); + glp_set_col_bnds(lp, j, GLP_LO, 0.0, 0.0); // p_j >= 0 + glp_set_obj_coef(lp, j, 0.0); + } + + const int nbCons = static_cast(_constraints.size()); + glp_add_rows(lp, nbCons); + + std::vector ia{0}, ja{0}; + std::vector ar{0.0}; // GLPK : indices 1-based, l'entrée 0 est ignorée + + for (int i = 0; i < nbCons; ++i) + { + const auto& coeffs = _constraints[i].getCoeffs(); // [a_0, a_1, ..., a_n] + + // a_0 + sum a_i p_i <= 0 (non stricte) + // a_0 + sum a_i p_i <= -1 (stricte, exact car tout est entier) + const double rhs = -coeffs[0] - (_constraints[i].getStrict() ? 1 : 0); + glp_set_row_bnds(lp, i + 1, GLP_UP, 0.0, rhs); + + for (int j = 1; j <= nbParam; ++j) + { + ia.push_back(i + 1); + ja.push_back(j); + ar.push_back(static_cast(coeffs[j])); + } + } + + glp_load_matrix(lp, static_cast(ia.size()) - 1, ia.data(), ja.data(), ar.data()); + + // 1) relaxation LP (rationnelle) -- requise par GLPK avant le MIP + glp_smcp lpParams; + glp_init_smcp(&lpParams); + lpParams.msg_lev = GLP_MSG_OFF; + glp_simplex(lp, &lpParams); + + bool empty; + if (glp_get_status(lp) != GLP_OPT) + { + empty = true; // infaisable même sur Q => infaisable sur N + } + else + { + // 2) résolution entière (branch & bound) + glp_iocp mipParams; + glp_init_iocp(&mipParams); + mipParams.msg_lev = GLP_MSG_OFF; + mipParams.presolve = GLP_ON; + glp_intopt(lp, &mipParams); + + const int status = glp_mip_status(lp); + empty = !(status == GLP_OPT || status == GLP_FEAS); + } + + glp_delete_prob(lp); + return empty; + } + + /** + * negation of a constraint: a0 + sum a_i p_i <= 0 (resp. < 0) + * becomes a0 + sum a_i p_i > 0 (resp. >= 0) + * i.e. -a0 - sum a_i p_i < 0 (resp. <= 0) + */ + static ParametricConstraint negate(const ParametricConstraint& c) + { + std::vector negCoeffs(c.getCoeffs().size()); + for (size_t k = 0; k < negCoeffs.size(); ++k) + negCoeffs[k] = -c.getCoeffs()[k]; + return ParametricConstraint{negCoeffs, !c.getStrict()}; + } + + /** + * removes redundant constraints: a constraint is redundant if the rest of the + * polyhedron already forbids violating it (i.e. rest ∧ ¬c is empty). + * checked against the progressively-reduced set so mutually-redundant pairs + * don't both get removed. + */ + void Polyhedron::minimize() + { + std::vector kept; + kept.reserve(_constraints.size()); + + for (size_t i = 0; i < _constraints.size(); ++i) { + Polyhedron test; + for (const auto& c : kept) test.addParametricConstraint(c); + for (size_t j = i + 1; j < _constraints.size(); ++j) + test.addParametricConstraint(_constraints[j]); + test.addParametricConstraint(negate(_constraints[i])); + + if (!test.isEmpty()) { + // there's a point satisfying everything else but violating c_i + // => c_i is NOT implied, we must keep it + kept.push_back(_constraints[i]); + } + // else: test is empty => impossible to violate c_i without violating + // something already kept/remaining => c_i is redundant, drop it + } + + _constraints = std::move(kept); + } + + + + + + /** + * + * @param i + * @param j + */ + void DBMatrix::close(uint32_t i, uint32_t j) + { + const auto dim = getDim(); + const auto dij = _matrix(i, j); + + for (uint32_t a = 0; a < dim; ++a) { + for (uint32_t b = 0; b < dim; ++b) { + auto candidate = _matrix(a, i) + dij + _matrix(j, b); + if (candidate.dominates(_matrix(a, b))) { + _matrix(a, b) = candidate; + } + } + } + } + + + /** + * resets a DBM to its unrestrained state. + */ + void Ppdbm::reset() + { + _dbm.reset(); + std::fill(_offsetCost.begin(), _offsetCost.end(), 0); + _rates.fill(0); + _PC.getConstraints().clear(); + _emptyZone = false; + } + + /** + * + */ + void Ppdbm::removeCostAtOffset() + { + for (uint32_t h = 1; h < _dim; h++) { + for (uint32_t k = 0; k <= _param; k++) { + _offsetCost[k] -= _rates(h,k) * (-_dbm(0, h).getValue()); + } + } + } + + /** + * + */ + void Ppdbm::addCostAtOffset() + { + for (uint32_t h = 1; h < _dim; h++) { + for (uint32_t k = 0; k <= _param; k++) { + _offsetCost[k] += _rates(h,k) * (-_dbm(0, h).getValue()); + } + } + } + + /** + * + * @param i + * @return + */ + bool Ppdbm::isDiagonalNegative(uint32_t i) const + { + const auto& d = _dbm(i, i); + return d.getValue() < 0 || (d.getValue() == 0 && d.getStrict()); + } + + /** + * + * @param i + * @param j + * @param constraint + * @return + */ + bool Ppdbm::constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint) + { + assert(i < _dim && j < _dim); + + // The constraint is useless? + if (_dbm(i,j).dominates(constraint)) { return true; } + + // The constraint empties the zone? + if (DbmBound::isIncompatible(constraint, _dbm(j,i))) { + _dbm(i, j) = constraint; + _emptyZone = true; + return false; + } + + // remove the cost at the previous offset + removeCostAtOffset(); + + // Add the constraint and close the DBM + _dbm(i, j) = constraint; + _dbm.close(i, j); + if (isDiagonalNegative(i)) { _emptyZone = true; return false; } + + // Add the cost at the new offset + addCostAtOffset(); + + return true; + } + + /** + * + * @param constraints + * @return + */ + bool Ppdbm::constrain_DBM_N(const std::vector& constraints) + { + if (constraints.empty()) return true; // nothing to do --> done + + // remove the cost at the previous offset + removeCostAtOffset(); + + bool empty = false; + for (auto it = constraints.begin(); it != constraints.end() && !empty; ++it) { + const auto& [i, j, bound] = *it; + + if (_dbm(i,j).dominates(bound)) { continue; } // useless constraint + + if (DbmBound::isIncompatible(bound, _dbm(j,i))) { + _dbm(i, j) = bound; + empty = true; + break; + } // constraint empties the zone + + _dbm(i, j) = bound; + _dbm.close(i, j); + empty = isDiagonalNegative(i); + } + + // add the cost of the new offset + addCostAtOffset(); + + if (empty) { _emptyZone = true; return false; } + return true; + } + + /** + * cost function of another DBM's offset calculated with actual rates. + * @param otherDbm + * @return + */ + std::vector Ppdbm::costAtOtherOffset(const DBMatrix& otherDbm) const + { + std::vector cost = _offsetCost; + + for (uint32_t x = 1; x < _dim; ++x) { + // difference in x coordinate of the offset + const int32_t offsetDifference = -otherDbm(0, x).getValue() - (-_dbm(0, x).getValue()); + + // multiply the difference by the associated rate and add it to the cost + for (uint32_t k = 0; k <= _param; ++k) { + cost[k] += _rates(x, k) * offsetDifference; + } + } + + return cost; + } + + /** + * + * @param q + */ + void Ppdbm::addConstantCost(const std::vector& q) + { + for (uint32_t k = 0; k <= _param; ++k) _offsetCost[k] += q[k]; + } + + + + + /////////////////////////////////////////////// Post operators //////////////////////////////////////////////////// + /** + * removes the strictness of lower bounds + */ + void DBMatrix::relaxDown() + { + for (uint32_t i = 1; i < getDim(); ++i) { + auto& b = _matrix(0, i); + if (!b.isInfinite() && b.getStrict()) b = DbmBound{b.getValue(), false}; + } + } + + /** + * removes the strictness of upper bounds + */ + void DBMatrix::relaxUp() + { + for (uint32_t i = 1; i < getDim(); ++i) { + auto& b = _matrix(i, 0); + if (!b.isInfinite() && b.getStrict()) b = DbmBound{b.getValue(), false}; + } + } + + /** + * + * @param clock + */ + void DBMatrix::relaxDownClock(uint32_t clock) + { + // relax (i, clock) for each i : "x_i - x_clock <= b" + for (uint32_t i = 0; i < getDim(); ++i) { + if (i == clock) continue; + auto& b = _matrix(i, clock); + if (!b.isInfinite() && b.getStrict()) b = DbmBound{b.getValue(), false}; + } + } + + /** + * + * @param clock + */ + void DBMatrix::relaxUpClock(uint32_t clock) + { + // relax (clock, i) for each i : "x_clock - x_i <= b" + for (uint32_t i = 0; i < getDim(); ++i) { + if (i == clock) continue; + auto& b = _matrix(clock, i); + if (!b.isInfinite() && b.getStrict()) b = DbmBound{b.getValue(), false}; + } + } + + /** + * finds the clocks that are linked together. ( x_i = x_j + constant) + * @return + */ + std::vector DBMatrix::findLinkedClocks() const + { + const auto dim = getDim(); + std::vector next(dim, 0); + for (uint32_t i = 0; i < dim; ++i) { + for (uint32_t j = i + 1; j < dim; ++j) { + if (!(*this)(i, j).isInfinite() && !(*this)(j, i).isInfinite() && + (*this)(i, j).getValue() + (*this)(j, i).getValue() == 0) { + // if DBM(i,j) + DBM(j,i) = 0 then clocks i and j are linked together by a constant. + next[i] = j; + break; + } + } + } + return next; + } + + /** + * if DBM(i,j) >= DBM(i,k) + DBM(k,j) then the facet is already covered by another one --> we dont consider it + * @param dbm + * @param i + * @param j + * @param next + * @return + */ + static bool isRedundant(const DBMatrix& dbm, uint32_t i, uint32_t j, const std::vector& next) + { + if (i == j) return true; + const auto bij = dbm(i, j); + if (bij.isInfinite()) return true; + + for (uint32_t k = 0; k < dbm.getDim(); ++k) { + if (k == i || k == j || next[k] != 0) continue; + const auto bik = dbm(i, k); + const auto bkj = dbm(k, j); + if (!bik.isInfinite() && !bkj.isInfinite() && + bij.getValue() >= bik.getValue() + bkj.getValue()) { + return true; + } + } + return false; + } + + /** + * returns the list of lower facets (defined by Z \land (x_0 - x_i <= -L)) + * @return + */ + std::vector Ppdbm::lowerFacets() const + { + std::vector result; + + Ppdbm relaxed = *this; + relaxed._dbm.relaxDown(); // closure + const auto next = relaxed._dbm.findLinkedClocks(); + + for (uint32_t i = 1; i < _dim; ++i) { + if (next[i] != 0) continue; + if (isRedundant(relaxed._dbm, 0, i, next)) + continue; // Facet covered by others + + Ppdbm facetZone = relaxed; + const auto lower = relaxed(0, i); // x0 - xi <= -L_i already set + const DbmBound pin{-lower.getValue(), false}; + const bool ok = facetZone.constrain_DBM(i, 0, pin); // xi - x0 <= L_i as well, so now xi = L_i + assert(ok); // should never empty the zone + + result.push_back(Facet{std::move(facetZone), DbmConstraint{i, 0, pin}}); + } + return result; + } + + /** + * returns the list of upper facets (defined by Z \land (x_i - x_0 <= U)) + * @return + */ + std::vector Ppdbm::upperFacets() const + { + std::vector result; + + Ppdbm relaxed = *this; + relaxed._dbm.relaxUp(); + const auto next = relaxed._dbm.findLinkedClocks(); + + for (uint32_t i = 1; i < _dim; ++i) { + if (next[i] != 0) continue; + if (isRedundant(relaxed._dbm, i, 0, next)) continue; + + Ppdbm facetZone = relaxed; + const auto upper = relaxed(i, 0); // xi - x0 <= U_i already set + const DbmBound pin{-upper.getValue(), false}; + const bool ok = facetZone.constrain_DBM(0, i, pin); // x0 - xi <= -U_i as well, so now xi = U_i + assert(ok); + + result.push_back(Facet{std::move(facetZone), DbmConstraint{i, 0, upper}}); + } + return result; + } + + /** + * + * @param clock + * @return + */ + std::vector Ppdbm::lowerFacetsRelativeTo(uint32_t clock) const + { + std::vector result; + + Ppdbm relaxed = *this; + relaxed._dbm.relaxDownClock(clock); + const auto next = relaxed._dbm.findLinkedClocks(); + + for (uint32_t i = 0; i < _dim; ++i) { + if (next[i] != 0) continue; + if (isRedundant(relaxed._dbm, i, clock, next)) continue; + + Ppdbm facetZone = relaxed; + const auto bound = relaxed(i, clock); // x_i - x_clock <= bound + const DbmBound pin{-bound.getValue(), false}; // x_clock - x_i <= -bound + const bool ok = facetZone.constrain_DBM(clock, i, pin); + assert(ok); + + result.push_back(Facet{std::move(facetZone), DbmConstraint{clock, i, pin}}); + } + return result; + } + + /** + * + * @param clock + * @return + */ + std::vector Ppdbm::upperFacetsRelativeTo(uint32_t clock) const + { + std::vector result; + + Ppdbm relaxed = *this; + relaxed._dbm.relaxUpClock(clock); + const auto next = relaxed._dbm.findLinkedClocks(); + + for (uint32_t i = 0; i < _dim; ++i) { + if (next[i] != 0) continue; + if (isRedundant(relaxed._dbm, clock, i, next)) continue; + + Ppdbm facetZone = relaxed; + const auto bound = relaxed(clock, i); // x_clock - x_i <= bound + const DbmBound pin{-bound.getValue(), false}; // x_i - x_clock <= -bound + const bool ok = facetZone.constrain_DBM(i, clock, pin); + assert(ok); + + result.push_back(Facet{std::move(facetZone), DbmConstraint{clock, i, bound}}); + } + return result; + } + + /** + * returns sum_{x \in X} r(x) + * @return + */ + CostVector Ppdbm::getSlope() const + { + CostVector slope(_param + 1,0); + for (uint32_t i = 1; i < _dim; ++i) { + for (uint32_t j = 0; j < _param + 1; ++j) { + slope[j] += _rates(i,j); + } + } + return slope; + } + + /** + * suppress the clock constraints of the form x_i <= w. Delay operation on a facet. + * @param clock + * @param c + */ + void Ppdbm::delay(uint32_t clock, const CostVector& c) + { + // modify the zone + for (uint32_t i = 1; i < getDim(); ++i) + _dbm(i, 0) = INF_BOUND; + + // adapt the price + for (uint32_t j = 0; j < _param + 1; ++j) { + // for every parameter, add the coeff of c and substract the coeffs of every other rate + _rates(clock, j) = c[j]; + for (uint32_t k = 1; k < _dim; ++k) + if (k != clock) + _rates(clock,j) -= _rates(k,j); + } + } + + /** + * suppress the clock constraints of the form x_i <= w. Delay operation on a facet. + */ + void Ppdbm::delay() + { + // modify the zone + for (uint32_t i = 1; i < getDim(); ++i) + _dbm(i, 0) = INF_BOUND; + } + + /** + * clock reset operation on a facet. + * assumes that clockToReset is the clock to reset and facetConstraint is the active constraint defining the facet + * on which we apply the reset, put in the form (clockToReset - anotherClock = value). + * @param clockToReset + * @param facetConstraint + */ + void Ppdbm::clockReset(uint32_t clockToReset, DbmConstraint facetConstraint) + { + assert(clockToReset < _dim); + + // modify the zone + _dbm(clockToReset, 0) = ZERO_BOUND; + _dbm(0, clockToReset) = ZERO_BOUND; + for (uint32_t i = 1; i < getDim(); ++i) { + _dbm(clockToReset, i) = _dbm(0,i); + _dbm(i,clockToReset) = _dbm(i,0); + } + + // adapt the price + for (uint32_t k = 0; k < _param + 1; ++k) { + _rates(facetConstraint.j,k) += _rates(facetConstraint.i,k); + _rates(facetConstraint.i,k) = 0; + } + } + + /** + * returns every successor by delay + * @param p + * @param J + * @return + */ + std::vector Ppdbm::post_delta(CostVector p, const std::vector& J) const + { + std::vector result; + + { + // case 1 : { (l, F↑p ∧ J) | F ∈ LF(Z) } and the constraint is p <= ∑r(x) thus p - ∑r(x) <= 0 + CostVector temp = getSlope(); + for (uint32_t i = 0; i <= _param; ++i) { + temp[i] = p[i] - temp[i]; + } + ParametricConstraint pconstraint{temp,false}; + + std::vector LF = lowerFacets(); + for (Facet f : LF ) { + auto [ppdbm, constraint] = f; + ppdbm.delay(constraint.i, p); + ppdbm.constrain_DBM_N(J); + ppdbm.getPC().addParametricConstraint(pconstraint); + result.push_back(std::move(ppdbm)); + } + } + { + // case 2 : { (l, Z) } ∪ { (l, F↑p ∧ J) | F ∈ UF(Z) } and the constraint is p > ∑r(x) thus ∑r(x) - p < 0 + CostVector temp = getSlope(); + for (uint32_t i = 0; i < _param + 1; ++i) { + temp[i] -= p[i]; + } + ParametricConstraint pconstraint{temp,true}; + + Ppdbm copy = *this; + copy.getPC().addParametricConstraint(pconstraint); + result.push_back(std::move(copy)); + + std::vector UF = upperFacets(); + for (Facet f : UF ) { + auto [ppdbm, constraint] = f; + ppdbm.delay(constraint.i, p); + ppdbm.constrain_DBM_N(J); + ppdbm.getPC().addParametricConstraint(pconstraint); + result.push_back(std::move(ppdbm)); + } + } + return result; + } + + /** + * returns every successor by transition e = (l, g, {x}, l'), with edge cost q. + * resetClock = std::nullopt si e ne reset aucune horloge. + */ + std::vector Ppdbm::post_e(const std::vector& g, uint32_t resetClock, const std::vector& q) const + { + std::vector result; + + // Z ∧ g + Ppdbm guarded = *this; + if (!g.empty() && !guarded.constrain_DBM_N(g)) return result; // g empties the zone + if (guarded.isEmpty()) return result; + + if (resetClock == 0) { + // if no reset : { (l', Z ∧ g + q) } + guarded.addConstantCost(q); + result.push_back(std::move(guarded)); + return result; + } + + const uint32_t x = resetClock; + CostVector rx(_param + 1); + for (uint32_t k = 0; k <= _param; ++k) rx[k] = guarded._rates(x,k); + + { + // case r(x) >= 0 <=> -r(x) <= 0 + CostVector temp(_param + 1); + for (uint32_t k = 0; k <= _param; ++k) temp[k] = -rx[k]; + ParametricConstraint pconstraint{temp, false}; + + for (Facet f : guarded.lowerFacetsRelativeTo(x)) { + auto [ppdbm, constraint] = f; + ppdbm.clockReset(x, constraint); + ppdbm.addConstantCost(q); + ppdbm.getPC().addParametricConstraint(pconstraint); + result.push_back(std::move(ppdbm)); + } + } + { + // case r(x) < 0 + ParametricConstraint pconstraint{rx, true}; + + for (Facet f : guarded.upperFacetsRelativeTo(x)) { + auto [ppdbm, constraint] = f; + ppdbm.clockReset(x, constraint); + ppdbm.addConstantCost(q); + ppdbm.getPC().addParametricConstraint(pconstraint); + result.push_back(std::move(ppdbm)); + } + } + + return result; + } + + + std::vector> Ppdbm::post(const std::pair& symbolicState, ParamPricedTimedAutomata automaton) + { + // will contain the successors by transition + delay + auto result = std::vector>{}; + // temporary for the successors by transition + auto temp = std::vector>{}; + for (auto edge : automaton._edges) { + auto [from, to, clockToReset, guards, costFunction] = edge; + if (from == symbolicState.first) { + auto successors = symbolicState.second.post_e(guards,clockToReset,costFunction); + for (auto successor : successors) { + temp.emplace_back(to,successor); + } + } + } + for (auto state : temp) { + auto successors = state.second.post_delta(automaton.getCostFunction(state.first),automaton.getInvariants(state.first)); + for (auto successor : successors) { + result.emplace_back(state.first,successor); + } + } + return result; + } + + + + ////////////////////////////////////////////// WQO relation //////////////////////////////////////////////////////// + + + /** + * compares each constraints, and if every constraint of one dominates the other, then subset. + * @param other + * @return + */ + relation_t DBMatrix::relation(const DBMatrix& other) const + { + assert(getDim() == other.getDim()); + const auto dim = getDim(); + + bool thisLeqOther = true; // this <= other everywhere -> this subset of other + bool otherLeqThis = true; // other <= this everywhere -> other subset of this + + for (uint32_t i = 0; i < dim; ++i) { + for (uint32_t j = 0; j < dim; ++j) { + if (!(*this)(i, j).dominates(other(i, j))) { thisLeqOther = false; } + if (!other(i, j).dominates((*this)(i, j))) { otherLeqThis = false; } + } + } + + if (thisLeqOther && otherLeqThis) { return base_EQUAL; } + if (thisLeqOther) { return base_SUBSET; } + if (otherLeqThis) { return base_SUPERSET; } + return base_DIFFERENT; + } + + /** + * computes the coefficients of the parameters in the cost function a the vertex v + * @param v + * @return + */ + CostVector Ppdbm::costAtVertex(const Vertex& v) const + { + assert(v.size() == _dim); + + CostVector cost = _offsetCost; // copy: offset coeffs + for (uint32_t x = 1; x < _dim; ++x) { + const int32_t offsetX = -_dbm(0, x).getValue(); // min value of clock h + const int32_t delta = v[x] - offsetX; + for (uint32_t k = 0; k <= _param; ++k) { + cost[k] += _rates(x, k) * delta; + } + } + return cost; + } + + /** + * is every coefficient of a lesser or equal than every coefficient ob b? + * @param a + * @param b + * @return + */ + bool Ppdbm::dominatesVector(const CostVector& a, const CostVector& b) + { + assert(a.size() == b.size()); + for (uint32_t k = 0; k < a.size(); ++k) { + if (a[k] > b[k]) { return false; } + } + return true; + } + + /** + * vector equality test + * @param a + * @param b + * @return + */ + bool Ppdbm::equalVector(const CostVector& a, const CostVector& b) + { + return a == b; + } + + /** + * Compares both versions of every cost function of a summit. If One DBM has lower coefficients everywhere, + * then better cost functions and domination. + * @param c1 + * @param c2 + * @return + */ + relation_t Ppdbm::compareCostLists(const std::vector& c1, const std::vector& c2) + { + assert(c1.size() == c2.size()); + + bool thisLeqOther = true; // cost of *this* <= cost of other, on every summit + bool otherLeqThis = true; // cost of other <= cost of *this*, on every summit + + for (uint32_t v = 0; v < c1.size(); ++v) { + if (!dominatesVector(c1[v], c2[v])) { thisLeqOther = false; } + if (!dominatesVector(c2[v], c1[v])) { otherLeqThis = false; } + if (!thisLeqOther && !otherLeqThis) { break; } // shortcut + } + + if (thisLeqOther && otherLeqThis) { return base_EQUAL; } + if (thisLeqOther) { return base_SUPERSET; } // *this* cheaper everywhere -> dominates + if (otherLeqThis) { return base_SUBSET; } // *this* more expensive everywhere -> dominated + return base_DIFFERENT; // no domination --> no branch cutting + } + + /** + * one of the most important functions, decides when the loops stop / when we cut branches of the exploration graph. + * is one or the other DBM not interesting to consider? meaning is there a zone who is a subset of the other and + * with less interesting cost functions ? (a zone already covered basically?) + * @param other + * @return + */ + relation_t Ppdbm::relation(const Ppdbm& other) const + { + assert(_dim == other._dim && _param == other._param); + + switch (_dbm.relation(other._dbm)) { + case base_DIFFERENT: + return base_DIFFERENT; + + case base_EQUAL: { + auto vertices = _dbm.verticesDBM(); // same zone on both sides + std::vector c1, c2; + c1.reserve(vertices.size()); + c2.reserve(vertices.size()); + for (const auto& v : vertices) { + c1.push_back(costAtVertex(v)); + c2.push_back(other.costAtVertex(v)); + } + return compareCostLists(c1, c2); + } + + case base_SUPERSET: { + // *this* is the bigger zone, other is the smaller + auto vertices = other._dbm.verticesDBM(); + std::vector c1, c2; + c1.reserve(vertices.size()); + c2.reserve(vertices.size()); + for (const auto& v : vertices) { + c1.push_back(costAtVertex(v)); + c2.push_back(other.costAtVertex(v)); + } + auto costRel = compareCostLists(c1, c2); + return (costRel == base_SUPERSET) ? base_SUPERSET : base_DIFFERENT; + } + + case base_SUBSET: { + // other is the bigger zone, *this* is the smaller: + auto vertices = _dbm.verticesDBM(); + std::vector c1, c2; + c1.reserve(vertices.size()); + c2.reserve(vertices.size()); + for (const auto& v : vertices) { + c1.push_back(costAtVertex(v)); + c2.push_back(other.costAtVertex(v)); + } + auto costRel = compareCostLists(c1, c2); + return (costRel == base_SUBSET) ? base_SUBSET : base_DIFFERENT; + } + + default: + return base_DIFFERENT; + } + } + + /** + * + * @param symbolicState + * @param symbolicStateListe + * @return + */ + bool Ppdbm::alreadyCoveredBy(std::pair symbolicState, + const std::vector>& symbolicStateListe) + { + auto [location, ppdbm] = symbolicState; + for (auto state : symbolicStateListe) { + auto [previousLocation, previousPpdbm] = state; + if (location == previousLocation && ppdbm.relation(previousPpdbm) == base_SUBSET) { return true; } + } + return false; + } + + + + + + + //////////////////////////////////////////////////// sommetsDBM() //////////////////////////////////////////////// + + + /** + * finds the vertices of a dbm! very important function! + * @return + */ + std::vector DBMatrix::verticesDBM() const + { + const uint32_t dim = getDim(); + std::vector vertices; + if (dim <= 1) { + vertices.push_back(Vertex{0}); return vertices; + } + + // Matrix of bound values + Matrix boundValues{dim,dim,0}; + for (uint32_t i = 0; i < dim; ++i) + for (uint32_t j = 0; j < dim; ++j) + if (i != j) boundValues(i,j) = _matrix(i,j).getValue(); + + // candidate edges {i,j} (i candidates; + for (uint32_t i = 0; i < dim; ++i) + for (uint32_t j = i + 1; j < dim; ++j) + if (boundValues(i,j) != INF || boundValues(j,i) != INF) + candidates.push_back({i, j}); + + // Union-find to build the spanning trees + std::vector parents(dim); + for (uint32_t i = 0; i < dim; ++i) {parents[i] = i;}; + + std::vector chosen; + + // + auto find = [&](uint32_t x) { + while (parents[x] != x) + x = parents[x]; + return x; + }; + + // evaluates a spanning tree : tries every possible direction + // on every edge, propagates from 0, checks global feasibility. + auto evaluateTree = [&](const std::vector& tree) { + const uint32_t nbedge = tree.size(); + + // we register the options of active constraint for each couple (x_i, x_j) with i> options(nbedge,2, {0, 0, INF}); + for (uint32_t k = 0; k < nbedge; ++k) { + auto [i, j] = tree[k]; + if (boundValues(i,j) != INF) options(k,0) = {i, j, boundValues(i,j)}; // x_i - x_j = w + if (boundValues(j,i) != INF) options(k,1) = {j, i, boundValues(j,i)}; // x_j - x_i = w + } + + // initial combination: first valid option for everyone + std::vector combination(nbedge); + for (uint32_t k = 0; k < nbedge; ++k) { + // option 0 is valid + if (std::get<2>(options(k, 0)) != INF) { combination[k] = 0; } + // otherwise use option 1 + else { combination[k] = 1; } + } + + // function to calculate the next combination of valid options for the edges + auto nextCombination = [&]() { + uint32_t k = 0; + while (k < nbedge) { + ++combination[k]; + // We reached the end of the options for this edge? + if (combination[k] >= 2) { combination[k] = 0; ++k; continue; } + // Skip this option if it does not exist. + auto [i, j, w] = options(k, combination[k]); + if (w != INF) return true; + } + return false; + }; + + // we try every possible combination of active constraints for couples (x_i, x_j) + for (;;) { + // directed adjacency matrix (to go from x to y, costs -3, and from y to x, costs +3...) + Matrix directedAdj(dim, dim, INF); + + // fill the adjacency matrix + for (uint32_t k = 0; k < nbedge; ++k) { + // example: x_0 - x_1 = 3. + auto [from, to, w] = options(k,combination[k]); + directedAdj(from,to) = -w; // x_1 = x_0 - 3. remove 3 to go from x0 to x1. + directedAdj(to, from) = w; // x_0 = x_1 + 3. add 3 to go from x1 to x0. + } + + // propagate the clock values through the spanning tree + std::vector val(dim, INF); + val[0] = 0; // the ref clock is always 0 + std::vector stack{0}; + + while (!stack.empty()) { + uint32_t u = stack.back(); + stack.pop_back(); + + for (uint32_t v = 0; v < dim; ++v) { + const int delta = directedAdj(u, v); + // if there is an active constraint linking u and v and v was not already visited + if (delta != INF && val[v] == INF) { + // visit v and add it to the stack + val[v] = val[u] + delta; + stack.push_back(v); + } + } + } + + // is every constraint verified with this configuration? + bool ok = true; + for (uint32_t i = 0; ok && i < dim; ++i) + for (uint32_t j = 0; ok && j < dim; ++j) + if (i != j && boundValues(i, j) != INF && val[i] - val[j] > boundValues(i, j)) + ok = false; + + // is the spanning tree covering every clock? + for (uint32_t i = 0; i < dim; ++i) { + if (val[i] == INF) { + ok = false; + break; + } + } + + // if everything fits, we found a vertex. + if (ok) { + Vertex v(dim); + for (uint32_t k = 0; k < dim; ++k) v[k] = val[k]; + vertices.push_back(std::move(v)); + } + + // then we check next combination + if (!nextCombination()) + break; + } + }; + + + // Backtracking : find all spanning trees on 'candidates' + std::function backtrack = [&](size_t start) { + if (chosen.size() == static_cast(dim - 1)) { evaluateTree(chosen); return; } + if (candidates.size() - start < (dim - 1 - chosen.size())) return; + + for (size_t k = start; k < candidates.size(); ++k) { + auto [i, j] = candidates[k]; + cindex_t ri = find(i), rj = find(j); + if (ri == rj) continue; // cycle + parents[ri] = rj; + chosen.push_back({i, j}); + backtrack(k + 1); + chosen.pop_back(); + parents[ri] = ri; + } + }; + backtrack(0); + + std::sort(vertices.begin(), vertices.end()); + vertices.erase(std::unique(vertices.begin(), vertices.end()), vertices.end()); + return vertices; + } + + + + /////////////////////////////////////////////////// result manipulation /////////////////////////////////////////// + + + /** + * is the feasible set of A included in that of B? + * true iff, for every constraint c of B, no point of A can violate c. + */ + static bool polyhedronIncludedIn(const std::vector& A, + const std::vector& B) + { + for (const auto& cb : B) { + Polyhedron test; + for (const auto& ca : A) test.addParametricConstraint(ca); + test.addParametricConstraint(negate(cb)); + if (!test.isEmpty()) return false; // found a point of A violating cb + } + return true; + } + + /** + * do two Ppdbm have the exact same cost function (offset + rates)? + */ + static bool sameCost(const Ppdbm& a, const Ppdbm& b) + { + return a.getOffsetCost() == b.getOffsetCost() && + a.getRates().getData() == b.getRates().getData(); + } + + /** + * removes results that are redundant: same cost function AND region included + * in another kept result's region. Keeps the more general one. + */ + std::vector Ppdbm::pruneRedundant(std::vector results) + { + std::vector removed(results.size(), false); + + for (size_t i = 0; i < results.size(); ++i) { + if (removed[i]) continue; + for (size_t j = i + 1; j < results.size(); ++j) { + if (removed[j]) continue; + if (!sameCost(results[i], results[j])) continue; + + const auto& pcI = results[i].getPC().getConstraints(); + const auto& pcJ = results[j].getPC().getConstraints(); + + const bool iInJ = polyhedronIncludedIn(pcI, pcJ); + const bool jInI = polyhedronIncludedIn(pcJ, pcI); + + if (iInJ && jInI) { + removed[j] = true; // identical regions, drop the duplicate + } else if (iInJ) { + removed[i] = true; // i is a special case of j, same cost -> drop i + break; // i is gone, stop comparing it further + } else if (jInI) { + removed[j] = true; // j is a special case of i, same cost -> drop j + } + // else: neither included in the other -> keep both + } + } + + std::vector kept; + for (size_t i = 0; i < results.size(); ++i) + if (!removed[i]) kept.push_back(std::move(results[i])); + return kept; + } + + using PCList = std::vector; + + static bool pcEmpty(const PCList& pc) + { + Polyhedron p; + for (const auto& c : pc) p.addParametricConstraint(c); + return p.isEmpty(); + } + + static PCList pcAnd(const PCList& a, const PCList& b) + { + PCList out = a; + out.insert(out.end(), b.begin(), b.end()); + return out; + } + + static CostVector costDiff(const CostVector& a, const CostVector& b) + { + CostVector r(a.size()); + for (size_t k = 0; k < a.size(); ++k) r[k] = a[k] - b[k]; + return r; + } + + /// A \ B : classic complement-of-convex decomposition into disjoint convex fragments. + static std::vector subtractRegion(const PCList& A, const PCList& B) + { + std::vector result; + PCList satisfiedSoFar; // constraints of B already forced true in previous branches + for (const auto& bc : B) { + PCList fragment = pcAnd(A, satisfiedSoFar); + fragment.push_back(negate(bc)); // this branch violates bc + if (!pcEmpty(fragment)) result.push_back(std::move(fragment)); + satisfiedSoFar.push_back(bc); + } + return result; + } + + std::vector resolveOptimalPartition(const std::vector& results) + { + struct Group { CostVector cost; std::vector fragments; }; + std::vector groups; + + for (const auto& r : results) { + const auto& cost = r.getOffsetCost(); + auto it = std::find_if(groups.begin(), groups.end(), + [&](const Group& g) { return g.cost == cost; }); + if (it == groups.end()) groups.push_back({cost, {r.getPC().getConstraints()}}); + else it->fragments.push_back(r.getPC().getConstraints()); + } + + // snapshot of each group's ORIGINAL region, used as the fixed conflict source + std::vector> original; + for (const auto& g : groups) original.push_back(g.fragments); + + auto refine = [&](std::vector& fragments, const std::vector& otherOriginal, + const ParametricConstraint& winCondition) { + for (const auto& otherFragment : otherOriginal) { + std::vector nextGen; + for (const auto& piece : fragments) { + PCList overlap = pcAnd(piece, otherFragment); + if (pcEmpty(overlap)) { nextGen.push_back(piece); continue; } + + for (auto& outside : subtractRegion(piece, otherFragment)) + nextGen.push_back(std::move(outside)); + + PCList winningPart = overlap; + winningPart.push_back(winCondition); + if (!pcEmpty(winningPart)) nextGen.push_back(std::move(winningPart)); + } + fragments = std::move(nextGen); + } + }; + + for (size_t i = 0; i < groups.size(); ++i) { + for (size_t j = i + 1; j < groups.size(); ++j) { + const CostVector diff = costDiff(groups[i].cost, groups[j].cost); + const ParametricConstraint iWinsOrTie{diff, false}; // cost_i <= cost_j + const ParametricConstraint jWinsStrict = negate(iWinsOrTie); // cost_i > cost_j + + refine(groups[i].fragments, original[j], iWinsOrTie); + refine(groups[j].fragments, original[i], jWinsStrict); + } + } + + std::vector out; + for (const auto& g : groups) + for (const auto& f : g.fragments) + if (!pcEmpty(f)) out.push_back({g.cost, f}); + return out; + } + + + void Polyhedron::printPolyhedron(const std::vector& constraints) + { + for (const auto& c : constraints) { + const auto& coeffs = c.getCoeffs(); + std::cout << " "; + bool first = true; + for (size_t k = 0; k < coeffs.size(); ++k) { + if (coeffs[k] == 0) continue; + if (!first) std::cout << " + "; + if (k == 0) std::cout << coeffs[k]; + else std::cout << coeffs[k] << "*p" << k; + first = false; + } + if (first) std::cout << "0"; + std::cout << (c.getStrict() ? " < 0" : " <= 0") << "\n"; + } + } + + std::vector OptimalRegion::pruneRedundantRegions(std::vector regions) + { + std::vector removed(regions.size(), false); + for (size_t i = 0; i < regions.size(); ++i) { + if (removed[i]) continue; + for (size_t j = i + 1; j < regions.size(); ++j) { + if (removed[j] || regions[i].cost != regions[j].cost) continue; + if (polyhedronIncludedIn(regions[i].constraints, regions[j].constraints)) { + removed[i] = true; break; + } + if (polyhedronIncludedIn(regions[j].constraints, regions[i].constraints)) { + removed[j] = true; + } + } + } + std::vector kept; + for (size_t i = 0; i < regions.size(); ++i) + if (!removed[i]) kept.push_back(std::move(regions[i])); + return kept; + } +}// namespace PPBDM + + + diff --git a/src/ParametricOptimalReachability.cpp b/src/ParametricOptimalReachability.cpp new file mode 100644 index 0000000..9675624 --- /dev/null +++ b/src/ParametricOptimalReachability.cpp @@ -0,0 +1,104 @@ +#include "../include/dbm/ParamPricedDBM.h" + +#include +#include +#include + + +int main() +{ + // the goal of what follows is to modelize the following automaton: --> (0, +p1) ------> (1,+p2) + // \ / + // x>=3 \ / x>=3 + // \ / + // \ / + // \/ \/ + // (2,+0) + uint32_t nbLocations = 3; + uint32_t initLocation = 0; + uint32_t nbClocks = 1; + uint32_t nbParam = 2; + + // x - x0 <= 3 + auto guard = PPDBM::DbmConstraint{0,1,PPDBM::DbmBound{-3,false}}; + auto guards = std::vector{}; + guards.emplace_back(guard); + + // resets (none for the example) + uint32_t clockToReset = 0; + + // edges + auto defaultCost = std::vector(nbParam + 1,0); + auto edge01 = PPDBM::Edge{0,1,clockToReset,std::vector{},defaultCost}; + auto edge02 = PPDBM::Edge{0,2,clockToReset,guards,defaultCost}; + auto edge12 = PPDBM::Edge{1,2,clockToReset,guards,defaultCost}; + auto edges = std::vector{edge01,edge02, edge12}; + + // invariants (none for the example) + auto inv = PPDBM::Matrix(nbLocations,nbClocks + 1,PPDBM::INF_BOUND); + + // costs + auto costs = PPDBM::Matrix(nbLocations,nbParam+1,0); + costs(0,1) = 1; + costs(1,2) = 1; + + // automaton + auto automaton = PPDBM::ParamPricedTimedAutomata{nbLocations,initLocation,nbClocks,nbParam,edges,inv,costs}; + + // algo: + auto RES = std::vector{}; + auto PASSED = std::vector>{}; + auto WAITING = std::vector>{std::make_pair(initLocation,PPDBM::Ppdbm{nbClocks + 1,nbParam,true})}; + auto [location, initPpdbm] = WAITING[0]; + for (auto ppdbm : initPpdbm.post_delta(automaton.getCostFunction(location), automaton.getInvariants(location))) { + WAITING.emplace_back(location,ppdbm); + } + + while (WAITING.size()>0) { + auto symbolicState = WAITING.back(); + WAITING.pop_back(); + + if (!symbolicState.second.getPC().isEmpty() && !PPDBM::Ppdbm::alreadyCoveredBy(symbolicState,PASSED)) { + PASSED.emplace_back(symbolicState); + auto successors = PPDBM::Ppdbm::post(symbolicState, automaton); + for (auto successor : successors) { + WAITING.emplace_back(successor); + } + if (symbolicState.first == 2) { + RES.emplace_back(symbolicState.second); + } + } + } + + + // results displaying + /* + // RES = PPDBM::Ppdbm::pruneRedundant(std::move(RES)); + for (auto res : RES) { + // res.getPC().minimize(); + res.display(); + } + */ + + RES = PPDBM::Ppdbm::pruneRedundant(std::move(RES)); + for (auto& res : RES) res.getPC().minimize(); + + auto optimal = PPDBM::resolveOptimalPartition(RES); + optimal = PPDBM::OptimalRegion::pruneRedundantRegions(std::move(optimal)); + + std::cout << "\n\n========== OPTIMAL PARTITION ==========\n"; + for (auto& region : optimal) { + PPDBM::Polyhedron p; + for (auto& c : region.constraints) p.addParametricConstraint(c); + p.minimize(); + + std::cout << "cost = [ "; + for (auto v : region.cost) std::cout << v << " "; + std::cout << "]\nregion:\n"; + PPDBM::Polyhedron::printPolyhedron(p.getConstraints()); + std::cout << "---------------------------------------\n"; + } + + return 0; +} + diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7dc05ec..6cdd62e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -9,11 +9,12 @@ foreach(source ${test_c_sources}) target_link_libraries(${test_target} PRIVATE ${libs}) endforeach() -file(GLOB test_cpp_sources test_fed.cpp test_fed_dbm.cpp test_fp_intersection.cpp test_valuation.cpp test_constraint.cpp) +file(GLOB test_cpp_sources test_fed.cpp test_fed_dbm.cpp test_fp_intersection.cpp test_valuation.cpp test_constraint.cpp test_ParamPricedDBM.cpp) foreach(source ${test_cpp_sources}) get_filename_component(test_target ${source} NAME_WE) add_executable(${test_target} ${source}) target_link_libraries(${test_target} ${libs} doctest_with_main) + target_compile_definitions(${test_target} PUBLIC DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN) endforeach() # comments contain expected time of Linux debug build, YMMV, whereas timeouts @@ -36,3 +37,5 @@ add_test(NAME test_allocation COMMAND test_allocation) add_test(NAME test_constraint COMMAND test_constraint) set_tests_properties(test_dbm_1_10 test_fed PROPERTIES TIMEOUT 1200) + +add_test(NAME test_ParamPricedDBM COMMAND test_ParamPricedDBM) diff --git a/test/test_ParamPricedDBM.cpp b/test/test_ParamPricedDBM.cpp new file mode 100644 index 0000000..7d07833 --- /dev/null +++ b/test/test_ParamPricedDBM.cpp @@ -0,0 +1,623 @@ +/* -*- mode: C++; c-file-style: "stroustrup"; c-basic-offset: 4; indent-tabs-mode: nil; -*- */ + +#include "dbm/ParamPricedDBM.h" + +#include + +#include +#include + +using namespace PPDBM; + +TEST_CASE("Matrix basic operations") +{ + Matrix m{2, 3, 7}; + CHECK(m.rows() == 2); + CHECK(m.cols() == 3); + for (uint32_t i = 0; i < 2; ++i) + for (uint32_t j = 0; j < 3; ++j) + CHECK(m(i, j) == 7); + + m(1, 2) = 42; + CHECK(m(1, 2) == 42); + CHECK(m(0, 0) == 7); // unchanged cell + + m.fill(3); + for (uint32_t i = 0; i < 2; ++i) + for (uint32_t j = 0; j < 3; ++j) + CHECK(m(i, j) == 3); +} + +TEST_CASE("DbmBound dominates and isIncompatible") +{ + const DbmBound b5{5, false}; // x - y <= 5 + const DbmBound b5strict{5, true}; // x - y < 5 + const DbmBound b3{3, false}; // x - y <= 3 + + // a tighter bound dominates a looser one + CHECK(b3.dominates(b5)); + CHECK_FALSE(b5.dominates(b3)); + + // same value: strict dominates non-strict + CHECK(b5strict.dominates(b5)); + CHECK_FALSE(b5.dominates(b5strict)); + + // a bound dominates itself + CHECK(b5.dominates(b5)); + + // infinite bound is dominated by anything finite + CHECK(INF_BOUND.isInfinite()); + CHECK_FALSE(b5.isInfinite()); + CHECK(b3.dominates(INF_BOUND)); + + // x - y <= 3 and y - x <= -5 => x <= y - 5 and x >= y + 3 : impossible + CHECK(DbmBound::isIncompatible(DbmBound{-5, false}, b3)); + CHECK_FALSE(DbmBound::isIncompatible(b3, b5)); + + // boundary case a+b == 0: incompatible only if one side is strict + CHECK_FALSE(DbmBound::isIncompatible(DbmBound{3, false}, DbmBound{-3, false})); + CHECK(DbmBound::isIncompatible(DbmBound{3, true}, DbmBound{-3, false})); +} + +TEST_CASE("DbmBound addition (transitive composition)") +{ + const DbmBound a{3, false}; + const DbmBound b{2, false}; + + const auto sum = a + b; + CHECK(sum.getValue() == 5); + CHECK_FALSE(sum.getStrict()); + + const DbmBound strictB{2, true}; + const auto sumStrict = a + strictB; + CHECK(sumStrict.getValue() == 5); + CHECK(sumStrict.getStrict()); // strictness propagates + + const auto sumInf = a + INF_BOUND; + CHECK(sumInf.isInfinite()); +} + +TEST_CASE("DBMatrix initial state") +{ + DBMatrix m{3}; + CHECK(m.getDim() == 3); + + for (uint32_t i = 0; i < 3; ++i) { + CHECK(m(i, i).getValue() == 0); + CHECK_FALSE(m(i, i).getStrict()); + CHECK(m(0, i).getValue() == 0); // every clock >= 0 by construction + } + // no upper bound has been set yet + CHECK(m(1, 0).isInfinite()); + CHECK(m(2, 0).isInfinite()); +} + +TEST_CASE("DBMatrix close propagates transitive bounds") +{ + DBMatrix m{3}; // clocks: 0 (reference), 1, 2 + m(1, 0) = DbmBound{3, false}; // x1 <= 3 + m(2, 1) = DbmBound{2, false}; // x2 - x1 <= 2 + + CHECK(m(2, 0).isInfinite()); // not propagated yet + + m.close(2, 1); + + CHECK(m(2, 0).getValue() == 5); // x2 <= x1 + 2 <= 5 + CHECK_FALSE(m(2, 0).getStrict()); +} + +TEST_CASE("DBMatrix relation between boxes") +{ + DBMatrix big{2}; // x1 in [0,5] + big(1, 0) = DbmBound{5, false}; + + DBMatrix small{2}; // x1 in [0,3] + small(1, 0) = DbmBound{3, false}; + + CHECK(big.relation(small) == base_SUPERSET); + CHECK(small.relation(big) == base_SUBSET); + CHECK(big.relation(big) == base_EQUAL); + + DBMatrix overlapA{2}; // x1 in [0,3] + overlapA(1, 0) = DbmBound{3, false}; + + DBMatrix overlapB{2}; // x1 in [2,5] + overlapB(1, 0) = DbmBound{5, false}; + overlapB(0, 1) = DbmBound{-2, false}; + + CHECK(overlapA.relation(overlapB) == base_DIFFERENT); + CHECK(overlapB.relation(overlapA) == base_DIFFERENT); +} + +TEST_CASE("DBMatrix verticesDBM on a trivial (single-clock) dimension") +{ + DBMatrix m{1}; + const auto vertices = m.verticesDBM(); + REQUIRE(vertices.size() == 1); + CHECK(vertices[0] == Vertex{0}); +} + +TEST_CASE("DBMatrix verticesDBM computes the corners of a box") +{ + DBMatrix m{3}; // reference clock 0, x1, x2 + m(1, 0) = DbmBound{5, false}; // x1 <= 5 + m(2, 0) = DbmBound{3, false}; // x2 <= 3 + // x1 >= 0 and x2 >= 0 come from the initial state + + auto vertices = m.verticesDBM(); + + std::vector expected{{0, 0, 0}, {0, 0, 3}, {0, 5, 0}, {0, 5, 3}}; + std::sort(expected.begin(), expected.end()); + + REQUIRE(vertices.size() == expected.size()); + CHECK(vertices == expected); +} + +TEST_CASE("Ppdbm constrain_DBM tightens bounds and detects emptiness") +{ + Ppdbm p{2, 0}; // 1 clock, no parameters + + CHECK(p.constrain_DBM(1, 0, DbmBound{5, false})); // x1 <= 5 + CHECK_FALSE(p.isEmpty()); + + // a weaker constraint is useless: nothing changes + CHECK(p.constrain_DBM(1, 0, DbmBound{10, false})); + CHECK(p(1, 0).getValue() == 5); + + // a tighter constraint is applied + CHECK(p.constrain_DBM(1, 0, DbmBound{4, false})); + CHECK(p(1, 0).getValue() == 4); + + // x1 >= 6 is incompatible with x1 <= 4: the zone becomes empty + CHECK_FALSE(p.constrain_DBM(0, 1, DbmBound{-6, false})); + CHECK(p.isEmpty()); +} + +TEST_CASE("Ppdbm constrain_DBM_N applies a batch of constraints atomically") +{ + Ppdbm p{3, 0}; + + const std::vector batch{ + {1, 0, DbmBound{4, false}}, // x1 <= 4 + {2, 0, DbmBound{6, false}}, // x2 <= 6 + {0, 1, DbmBound{-1, false}}, // x1 >= 1 + }; + + CHECK(p.constrain_DBM_N(batch)); + CHECK_FALSE(p.isEmpty()); + CHECK(p(1, 0).getValue() == 4); + CHECK(p(2, 0).getValue() == 6); + CHECK(p(0, 1).getValue() == -1); + + // x1 >= 10 is incompatible with the x1 <= 4 constraint set above + const std::vector contradiction{ + {0, 1, DbmBound{-10, false}}, + }; + CHECK_FALSE(p.constrain_DBM_N(contradiction)); + CHECK(p.isEmpty()); +} + +TEST_CASE("Ppdbm reset returns the zone to its unconstrained state") +{ + auto p = Ppdbm{2, 0}; + p.constrain_DBM(1, 0, DbmBound{5, false}); + p.constrain_DBM(0, 1, DbmBound{-100, false}); // makes the zone empty + REQUIRE(p.isEmpty()); + + p.reset(); + + CHECK_FALSE(p.isEmpty()); + CHECK(p(1, 0).isInfinite()); + CHECK(p(0, 1).getValue() == 0); +} + +TEST_CASE("Ppdbm costAtVertex reflects rates and the zone's lower bound") +{ + Ppdbm p{2, 1}; // 1 clock, 1 parameter + + Matrix rates{2, 2, 0}; + rates(1, 0) = 1; + rates(1, 1) = 1; + p.getRates() = rates; + + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); // x1 <= 5 + REQUIRE(p.constrain_DBM(0, 1, DbmBound{-2, false})); // x1 >= 2 + + // at the lower corner of the zone, delta == 0: cost == offset cost + CHECK(p.costAtVertex(Vertex{0, 2}) == CostVector{2, 2}); + + // at the upper corner, cost == offset + rate * (5 - 2) + CHECK(p.costAtVertex(Vertex{0, 5}) == CostVector{5, 5}); +} + +TEST_CASE("Ppdbm costAtOtherOffset translates the cost to another zone's offset") +{ + Ppdbm p{2, 1}; + Matrix rates{2, 2, 0}; + rates(1, 0) = 1; + rates(1, 1) = 1; + p.getRates() = rates; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); + REQUIRE(p.constrain_DBM(0, 1, DbmBound{-2, false})); // x1 in [2,5] + + DBMatrix otherDbm{2}; // x1 in [0,5] : offset is 0 instead of 2 + otherDbm(1, 0) = DbmBound{5, false}; + + const auto cost = p.costAtOtherOffset(otherDbm); + CHECK(cost == CostVector{0, 0}); +} + +TEST_CASE("Ppdbm::dominatesVector and equalVector") +{ + const CostVector a{1, 2, 3}; + const CostVector b{1, 2, 4}; + const CostVector c{2, 1, 4}; + + CHECK(Ppdbm::dominatesVector(a, a)); + CHECK(Ppdbm::dominatesVector(a, b)); // a <= b on every coefficient + CHECK_FALSE(Ppdbm::dominatesVector(b, a)); + CHECK_FALSE(Ppdbm::dominatesVector(a, c)); // incomparable vectors + CHECK_FALSE(Ppdbm::dominatesVector(c, a)); + + CHECK(Ppdbm::equalVector(a, a)); + CHECK_FALSE(Ppdbm::equalVector(a, b)); +} + +TEST_CASE("Ppdbm::compareCostLists aggregates per-vertex domination") +{ + const std::vector cheaper{{1, 1}, {2, 2}}; + const std::vector pricier{{2, 2}, {3, 3}}; + const std::vector mixed{{0, 5}, {5, 0}}; + + CHECK(Ppdbm::compareCostLists(cheaper, cheaper) == base_EQUAL); + CHECK(Ppdbm::compareCostLists(cheaper, pricier) == base_SUPERSET); // cheaper dominates + CHECK(Ppdbm::compareCostLists(pricier, cheaper) == base_SUBSET); + CHECK(Ppdbm::compareCostLists(cheaper, mixed) == base_DIFFERENT); +} + +TEST_CASE("Ppdbm::relation compares zones and, when equal, their costs") +{ + Ppdbm a{2, 1}; + Matrix ratesA{2, 2, 0}; + ratesA(1, 0) = 1; + ratesA(1, 1) = 1; + a.getRates() = ratesA; + REQUIRE(a.constrain_DBM(1, 0, DbmBound{5, false})); + REQUIRE(a.constrain_DBM(0, 1, DbmBound{-2, false})); // x1 in [2,5] + + Ppdbm b{2, 1}; + Matrix ratesB{2, 2, 0}; + ratesB(1, 0) = 2; // strictly more expensive on the first coefficient + ratesB(1, 1) = 1; + b.getRates() = ratesB; + REQUIRE(b.constrain_DBM(1, 0, DbmBound{5, false})); + REQUIRE(b.constrain_DBM(0, 1, DbmBound{-2, false})); // same zone as a + + // same zone, a is cheaper everywhere -> a dominates (SUPERSET), b is dominated (SUBSET) + CHECK(a.relation(b) == base_SUPERSET); + CHECK(b.relation(a) == base_SUBSET); + CHECK(a.relation(a) == base_EQUAL); + + Ppdbm c{2, 1}; + Matrix ratesC{2, 2, 0}; + ratesC(1, 0) = 1; + ratesC(1, 1) = 1; + c.getRates() = ratesC; + REQUIRE(c.constrain_DBM(1, 0, DbmBound{3, false})); // x1 in [0,3] + + // overlapping but neither zone contains the other -> DIFFERENT, regardless of cost + CHECK(a.relation(c) == base_DIFFERENT); +} + +//////////////////////////////////// New functions: relax / findLinkedClocks ////////////////////////////////// + +TEST_CASE("DBMatrix relaxDown and relaxUp remove strictness only") +{ + DBMatrix m{3}; + m(0, 1) = DbmBound{-2, true}; // x1 > 2 (strict lower bound) + m(0, 2) = DbmBound{0, false}; // x2 >= 0 (already weak) + m(1, 0) = DbmBound{5, true}; // x1 < 5 (strict upper bound) + m(2, 0) = DbmBound{3, false}; // x2 <= 3 (already weak) + + m.relaxDown(); + CHECK_FALSE(m(0, 1).getStrict()); + CHECK(m(0, 1).getValue() == -2); // the value does not change, only the strictness is removed + CHECK_FALSE(m(0, 2).getStrict()); // already weak: unchanged + CHECK(m(1, 0).getStrict()); // relaxDown does not touch upper bounds + + m.relaxUp(); + CHECK_FALSE(m(1, 0).getStrict()); + CHECK(m(1, 0).getValue() == 5); + CHECK_FALSE(m(2, 0).getStrict()); + + // an infinite bound remains infinite and is never "relaxed" + DBMatrix m2{2}; + CHECK(m2(1, 0).isInfinite()); + m2.relaxUp(); + CHECK(m2(1, 0).isInfinite()); +} + +TEST_CASE("DBMatrix relaxDownClock and relaxUpClock only touch the given clock's column/row") +{ + DBMatrix m{3}; + m(0, 1) = DbmBound{-2, true}; // x1 > 2 + m(2, 1) = DbmBound{4, true}; // x2 - x1 < 4 + m(1, 0) = DbmBound{5, true}; // x1 < 5 + m(1, 2) = DbmBound{1, true}; // x1 - x2 < 1 + + m.relaxDownClock(1); // relax (i,1) for every i != 1 + CHECK_FALSE(m(0, 1).getStrict()); + CHECK_FALSE(m(2, 1).getStrict()); + CHECK(m(1, 0).getStrict()); // untouched: this is not an entry (*, 1) + CHECK(m(1, 2).getStrict()); // same + + m.relaxUpClock(1); // relax (1,i) for every i != 1 + CHECK_FALSE(m(1, 0).getStrict()); + CHECK_FALSE(m(1, 2).getStrict()); +} + +TEST_CASE("DBMatrix findLinkedClocks detects zero cycles") +{ + DBMatrix m{3}; + // x1 - x2 = 3 exactly: x1-x2<=3 and x2-x1<=-3 + m(1, 2) = DbmBound{3, false}; + m(2, 1) = DbmBound{-3, false}; + + const auto next = m.findLinkedClocks(); + REQUIRE(next.size() == 3); + CHECK(next[1] == 2); // 1 is linked to 2 + CHECK(next[0] == 0); // the reference is not linked to anything here + CHECK(next[2] == 0); // 2 is a representative (no j>2 to check) + + // if the difference is not exactly zero, there is no link + DBMatrix m2{3}; + m2(1, 2) = DbmBound{3, false}; + m2(2, 1) = DbmBound{-2, false}; // x1-x2 in [2,3], not fixed + CHECK(m2.findLinkedClocks() == std::vector{0, 0, 0}); +} + +//////////////////////////////////////////// lowerFacets / upperFacets ////////////////////////////////////////////// + +TEST_CASE("Ppdbm lowerFacets and upperFacets on an independent box") +{ + Ppdbm p{3, 0}; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); // x1 <= 5 + REQUIRE(p.constrain_DBM(2, 0, DbmBound{3, false})); // x2 <= 3 + // x1, x2 >= 0 by default, no link between x1 and x2 + + const auto lower = p.lowerFacets(); + REQUIRE(lower.size() == 2); // one facet per clock, no redundancy + for (const auto& f : lower) { + CHECK(f._ppdbm(f._constraint.i, 0).getValue() == 0); // each facet fixes the clock to 0 + } + + const auto upper = p.upperFacets(); + REQUIRE(upper.size() == 2); + bool sawX1 = false, sawX2 = false; + for (const auto& f : upper) { + if (f._constraint.i == 1) { CHECK(f._ppdbm(1, 0).getValue() == 5); sawX1 = true; } + if (f._constraint.i == 2) { CHECK(f._ppdbm(2, 0).getValue() == 3); sawX2 = true; } + } + CHECK(sawX1); + CHECK(sawX2); +} + +TEST_CASE("Ppdbm upperFacets discards a redundant facet") +{ + // x1 <= 5, x2 <= 5, x1 - x2 <= 0 => x1 <= 5 is implied by x2<=5 and x1<=x2 + Ppdbm p{3, 0}; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); + REQUIRE(p.constrain_DBM(2, 0, DbmBound{5, false})); + REQUIRE(p.constrain_DBM(1, 2, DbmBound{0, false})); + + const auto upper = p.upperFacets(); + REQUIRE(upper.size() == 1); // only x2 <= 5 is a true facet + CHECK(upper[0]._constraint.i == 2); + CHECK(upper[0]._ppdbm(2, 0).getValue() == 5); +} + +///////////////////////////////////// lowerFacetsRelativeTo / upperFacetsRelativeTo /////////////////////////////////// + +TEST_CASE("Ppdbm lowerFacetsRelativeTo and upperFacetsRelativeTo isolate the reset clock's own bounds") +{ + Ppdbm p{3, 0}; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); // x1 <= 5 + REQUIRE(p.constrain_DBM(0, 1, DbmBound{-2, false})); // x1 >= 2 + REQUIRE(p.constrain_DBM(2, 0, DbmBound{3, false})); // x2 <= 3, x2 >= 0, no link between x1 and x2 + + const auto lower = p.lowerFacetsRelativeTo(1); + REQUIRE(lower.size() == 1); // only the reference (0) gives a non-redundant facet + CHECK(lower[0]._constraint.i == 1); + CHECK(lower[0]._constraint.j == 0); + CHECK(lower[0]._ppdbm(1, 0).getValue() == 2); // x1 fixed to its lower bound + + const auto upper = p.upperFacetsRelativeTo(1); + REQUIRE(upper.size() == 1); + CHECK(upper[0]._constraint.i == 1); + CHECK(upper[0]._constraint.j == 0); + CHECK(upper[0]._ppdbm(1, 0).getValue() == 5); // x1 fixed to its upper bound +} + +////////////////////////////////////////////////////// getSlope ///////////////////////////////////////////////////// + +TEST_CASE("Ppdbm getSlope sums the rates of every clock") +{ + Ppdbm p{3, 1}; + Matrix rates{3, 2, 0}; + rates(1, 0) = 2; + rates(1, 1) = 3; + rates(2, 0) = 1; + rates(2, 1) = -1; + p.getRates() = rates; + + CHECK(p.getSlope() == CostVector{3, 2}); +} + +////////////////////////////////////////////////////////// delay //////////////////////////////////////////////////// + +TEST_CASE("Ppdbm delay(clock, c) removes upper bounds and rebalances the rate so the slope equals c") +{ + Ppdbm p{3, 1}; + Matrix rates{3, 2, 0}; + rates(1, 0) = 1; rates(1, 1) = 0; + rates(2, 0) = 2; rates(2, 1) = 1; + p.getRates() = rates; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); + REQUIRE(p.constrain_DBM(2, 0, DbmBound{3, false})); + + p.delay(1, CostVector{5, 7}); + + CHECK(p(1, 0).isInfinite()); + CHECK(p(2, 0).isInfinite()); + // the rate of clock1 absorbs the difference: the total slope becomes exactly c + CHECK(p.getSlope() == CostVector{5, 7}); +} + +TEST_CASE("Ppdbm delay() only removes upper bounds") +{ + Ppdbm p{2, 0}; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); + + p.delay(); + + CHECK(p(1, 0).isInfinite()); + CHECK(p(0, 0).getValue() == 0); // the reference diagonal is not touched + CHECK_FALSE(p(0, 0).getStrict()); + CHECK(p(0, 1).getValue() == 0); // x1 >= 0 is still true +} + +//////////////////////////////////////////////////// clockReset ///////////////////////////////////////////////////// + +TEST_CASE("Ppdbm clockReset applied on a facet from lowerFacetsRelativeTo") +{ + Ppdbm p{3, 0}; + Matrix rates{3, 1, 0}; + rates(1, 0) = 2; + rates(2, 0) = 3; + p.getRates() = rates; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); // x1 <= 5 + REQUIRE(p.constrain_DBM(0, 1, DbmBound{-2, false})); // x1 >= 2 + REQUIRE(p.constrain_DBM(2, 0, DbmBound{3, false})); // x2 <= 3 + + CHECK(p.getSlope() == CostVector{5}); // 2 + 3 + + const auto facets = p.lowerFacetsRelativeTo(1); + REQUIRE(facets.size() == 1); + Facet f = facets[0]; // copy: x1 is already fixed to 2 on this facet + CHECK(f._ppdbm(1, 0).getValue() == 2); + + f._ppdbm.clockReset(1, f._constraint); + + // the zone: x1 is reset to 0, and its relation to x2 is copied from the reference + CHECK(f._ppdbm(1, 0).getValue() == 0); + CHECK_FALSE(f._ppdbm(1, 0).getStrict()); + CHECK(f._ppdbm(0, 1).getValue() == 0); + CHECK(f._ppdbm(1, 2).getValue() == f._ppdbm(0, 2).getValue()); + CHECK(f._ppdbm(1, 2).getStrict() == f._ppdbm(0, 2).getStrict()); + CHECK(f._ppdbm(2, 1).getValue() == f._ppdbm(2, 0).getValue()); + CHECK(f._ppdbm(2, 1).getStrict() == f._ppdbm(2, 0).getStrict()); + + // the rate of x1 is transferred (here to the reference, so it is "lost" for future cost calculations) + CHECK(f._ppdbm.getSlope() == CostVector{3}); // 5 - 2: only rates(x2) remains observable +} + +/////////////////////////////////////////////////////// post_delta ////////////////////////////////////////////////// + +TEST_CASE("Ppdbm post_delta produces the three symbolic successors of the formula") +{ + Ppdbm p{2, 0}; // one clock, no parameter (constants) + Matrix rates{2, 1, 0}; + rates(1, 0) = 2; + p.getRates() = rates; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); // x1 in [0,5] + + const auto result = p.post_delta(CostVector{3}, {}); // empty J: tests the constrain_DBM_N(empty) path + REQUIRE(result.size() == 3); + + // successor from LF(Z) (case p <= slope): x1 becomes unbounded, the slope becomes p + const auto& fromLower = result[0]; + CHECK(fromLower(1, 0).isInfinite()); + CHECK(fromLower(0, 1).getValue() == 0); + CHECK(fromLower.getPC().getConstraints().size() == 1); + CHECK_FALSE(fromLower.getPC().getConstraints()[0].getStrict()); // p - slope <= 0 + + // successor { (l, Z) } for the case p > slope: zone and slope unchanged + const auto& untouched = result[1]; + CHECK(untouched(1, 0).getValue() == 5); + CHECK(untouched.getPC().getConstraints().size() == 1); + CHECK(untouched.getPC().getConstraints()[0].getStrict()); // slope - p < 0 + + // successor from UF(Z) (case p > slope): x1 started at 5, becomes unbounded, same slope p + const auto& fromUpper = result[2]; + CHECK(fromUpper(1, 0).isInfinite()); + CHECK(fromUpper(0, 1).getValue() == -5); + CHECK(fromUpper.getPC().getConstraints().size() == 1); + CHECK(fromUpper.getPC().getConstraints()[0].getStrict()); +} + +///////////////////////////////////////////////////////// post_e //////////////////////////////////////////////////// + +TEST_CASE("Ppdbm post_e without reset applies the guard and adds the constant cost") +{ + Ppdbm p{2, 0}; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); // x1 <= 5 + + const std::vector g{{1, 0, DbmBound{3, false}}}; // x1 <= 3 + const auto result = p.post_e(g, std::nullopt, CostVector{7}); + + REQUIRE(result.size() == 1); + CHECK_FALSE(result[0].isEmpty()); + CHECK(result[0](1, 0).getValue() == 3); + CHECK(result[0].getOffsetCost() == CostVector{7}); +} + +TEST_CASE("Ppdbm post_e returns nothing when the guard empties the zone") +{ + Ppdbm p{2, 0}; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{3, false})); // x1 <= 3 + REQUIRE(p.constrain_DBM(0, 1, DbmBound{-1, false})); // x1 >= 1 + + const std::vector g{{0, 1, DbmBound{-10, false}}}; // x1 >= 10: incompatible + const auto result = p.post_e(g, std::nullopt, CostVector{0}); + + CHECK(result.empty()); +} + +TEST_CASE("Ppdbm post_e with a reset generates the lower- and upper-facet branches") +{ + Ppdbm p{3, 0}; + Matrix rates{3, 1, 0}; + rates(1, 0) = 2; + rates(2, 0) = 3; + p.getRates() = rates; + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); // x1 <= 5 + REQUIRE(p.constrain_DBM(0, 1, DbmBound{-2, false})); // x1 >= 2 + REQUIRE(p.constrain_DBM(2, 0, DbmBound{3, false})); // x2 <= 3 + + CHECK(p.getOffsetCost() == CostVector{4}); // rates(x1)=2 * lowerbound(x1)=2 + + const auto result = p.post_e({}, 1u, CostVector{10}); // no guard, reset x1, edge cost 10 + + REQUIRE(result.size() == 2); // one branch per facet relative to x1 (lower then upper) + + for (const auto& r : result) { + // in both cases x1 is reset to 0, and the zone only depends on x2 + CHECK(r(1, 0).getValue() == 0); + CHECK_FALSE(r(1, 0).getStrict()); + CHECK(r(0, 1).getValue() == 0); + CHECK(r(1, 2).getValue() == r(0, 2).getValue()); + CHECK(r(1, 2).getStrict() == r(0, 2).getStrict()); + CHECK(r(2, 1).getValue() == r(2, 0).getValue()); + CHECK(r(2, 1).getStrict() == r(2, 0).getStrict()); + CHECK(r.getSlope() == CostVector{3}); // the rate of x1 (2) is removed from the observable slope + CHECK(r.getPC().getConstraints().size() == 1); + } + + // branch from lowerFacetsRelativeTo (x1 was 2 before reset): cost = 4 (already accumulated) + 10 + CHECK(result[0].getOffsetCost() == CostVector{14}); + CHECK_FALSE(result[0].getPC().getConstraints()[0].getStrict()); + + // branch from upperFacetsRelativeTo (x1 was 5 before reset): cost = 10 (already accumulated) + 10 + CHECK(result[1].getOffsetCost() == CostVector{20}); + CHECK(result[1].getPC().getConstraints()[0].getStrict()); +}