From 29bc8e8e100435646907cbb02883c0b0a0d34e7e Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Wed, 19 Aug 2026 10:05:47 +0200 Subject: [PATCH 01/18] small modifs to run it on my computer --- CMakeLists.txt | 2 +- cmake/boost.cmake | 3 +++ cmake/doctest.cmake | 4 ++-- test/CMakeLists.txt | 1 + 4 files changed, 7 insertions(+), 3 deletions(-) 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/test/CMakeLists.txt b/test/CMakeLists.txt index 7dc05ec..4f7db19 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -14,6 +14,7 @@ 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 From 4aac36fd6a2c9abfb71e49e638462543aa461171 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Thu, 20 Aug 2026 11:11:35 +0200 Subject: [PATCH 02/18] PPDBM data structure first draft --- include/dbm/ParamPricedDBM.h | 763 +++++++++++++++++++++++++++++++++++ src/ParamPricedDBM.cpp | 0 2 files changed, 763 insertions(+) create mode 100644 include/dbm/ParamPricedDBM.h create mode 100644 src/ParamPricedDBM.cpp diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h new file mode 100644 index 0000000..046724b --- /dev/null +++ b/include/dbm/ParamPricedDBM.h @@ -0,0 +1,763 @@ +// -*- mode: C++; c-file-style: "stroustrup"; c-basic-offset: 4; indent-tabs-mode: nil; -*- +//////////////////////////////////////////////////////////////////// +// +// This file is a part of the UPPAAL toolkit. +// Copyright (c) 1995 - 2005, Uppsala University and Aalborg University. +// All right reserved. +// +// $Id: priced.h,v 1.15 2005/05/31 20:54:59 behrmann Exp $ +// +/////////////////////////////////////////////////////////////////// + +#ifndef INCLUDE_DBM_PARAMPRICEDDBM_H +#define INCLUDE_DBM_PARAMPRICEDDBM_H + +#include "dbm/dbm.h" +#include "dbm/mingraph.h" +#include "base/inttypes.h" + +#include +#include +#include +#include + +/** + * @file + * + * Functions for handling priced DBMs. + * + * A priced DBM is a data structure representing a zone with an affine + * hyperplane over the zone assigning a cost to each point in the + * zone. + * + * A few facts about priced DBMs: + * + * - Priced DBMs have a positive dimension. + * + * - Priced DBMs are stored continously in memory. + * + * - The clock with index zero is the reference clock. + * + * - The affine hyperplane is given by the cost in the offset point + * and the coefficients of the hyperplane. + * + * - A closed DBM is never empty. + * + * Allocation can be handled by the library by using \c + * pdbm_allocate() and \c pdbm_deallocate(). This can optionally be + * combined with reference counting using \c pdbm_incRef() and \c + * pdbm_decRef(), in which case \c pdbm_deallocate() is called + * automatically as soon as the reference count reaches zero. + * Reference counting is used to implement copy on write + * semantics. Thus whenever you modify a priced DBM with a reference + * count larger than 1, then that DBM is automatically copied. For + * this reason, all functions modifying priced DBMs take a reference + * parameter to the priced DBM to modify. + * + * A NULL pointer is equivalent to an empty DBM with a reference count + * of 1. Functions that can cause the priced DBM to become empty can + * choose to deallocate the DBM rather than to copy it. All functions + * that do not require the input DBM to be closed will allocate a new + * DBM when given a NULL pointer as input. + * + * Alternatively, memory can be allocated outside the library. In that + * case \c pdbm_size() bytes need to be allocated and preinitialized + * for the use by the library by calling \c pdbm_reserve(). Priced + * DBMs allocated in this manner must not be reference counted. + */ + +/** + * Data type for priced dbm. + */ +struct PDBM_s; +using PDBMPtr = std::shared_ptr; +using PDBMCPtr = std::shared_ptr; + +/** + * Computes the size in number of bytes of a priced dbm of the given + * dimension. + * + * @param dim is a dimension. + * @pre dim > 0 + * @return Amount of memory in bytes. + */ +size_t pdbm_size(cindex_t dim); + +/** + * Reserves a memory area for use as a priced DBM. The priced DBM will + * be unitialised, hence further initialised with e.g. \c pdbm_init() + * or \c pdbm_zero() is required. + * + * @param dim is the dimension of the priced DBM to reserve. + * @param p is a memory area of size \c pdbm_size(dim). + * + * @return An unitialised priced DBM of dimension \a dim. + */ +PDBMPtr pdbm_reserve(cindex_t dim, void* p); + +/** + * Allocates a new priced DBM. The reference count is initialised to + * 0. No other initialisation is performed. + * + * @param dim is the dimension. + * @return A newly allocated priced DBM of dimension \a dim. + * @pre dim is larger than 0. + */ +PDBMPtr pdbm_allocate(cindex_t dim); + +/** + * Deallocates a priced DBM. + * + * @pre + * - \a pdbm was allocated with \c pdbm_allocate(). + * - The reference count of \a pdbm is zero. + */ +void pdbm_deallocate(PDBM_s* pdbm); + +/** + * Copy a priced DBM. + * + * Priced DBMs normally use a copy-on-write scheme, thus is no need to + * call this function. The only use of it (besides internally in the + * library as part of the copy-on-write implementation) is for priced + * DBMs, for which you choose not to use reference counting. Hence, a + * call to this function is only permissible if the reference count of + * \a dst is zero. If \a dst is NULL, a new DBM is allocate with \c + * pdbm_allocate(). + * + * @param dst The destination. + * @param src The source. + * @param dim The dimension of \a dst and \a src. + * @pre The reference count of \a dst is zero or dst is NULL. + * @post The reference count of the return value is zero. + * @return The destination. + */ +PDBMPtr pdbm_copy(PDBMPtr dst, const PDBMCPtr& src, cindex_t dim); + +/** + * Initialises a priced DBM to the DBM containing all valuations. The + * cost of all valuations will be zero. If \a pdbm is NULL, a new DBM + * is allocated with \c pdbm_allocate(). + * + * @param pdbm is the priced DBM to initialize. + * @param dim is the dimension of \a pdbm. + */ +void pdbm_init(PDBMPtr& pdbm, cindex_t dim); + +/** + * Initialize a priced DBM to only contain the origin with a cost of + * 0. If \a pdbm is NULL, a new DBM is allocated with \c + * pdbm_allocate(). + * + * @param pdbm is the priced DBM to initialise. + * @param dim is the dimension of \a pdbm. + */ +void pdbm_zero(PDBMPtr& pdbm, cindex_t dim); + +/** + * Constrain a priced DBM. + * + * After the call, the DBM will be constrained such that the + * difference between clock \a i and clock \a j is smaller than \a + * constraint. + * + * @see dbm_constrain1 + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param i is a clock index. + * @param j is a clock index. + * @param constraint is a raw_t bound. + * @pre i != j + * @post The DBM is empty or closed. + * @return True if and only if the result is not empty. + */ +bool pdbm_constrain1(PDBMPtr& pdbm, cindex_t dim, cindex_t i, cindex_t j, raw_t constraint); + +/** + * Constrain a priced DBM with multiple constraints. + * + * @see dbm_constraint1 + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param constraints is an array of \a n constraints. + * @param n is the length of \a constraints. + * @pre the constraints are valid + * @post The DBM is empty or closed. + * @return + * - true if the DBM is non empty + + * closed non empty DBM + updated consistent rates + * - or false if the DBM is empty + + * empty DBM + inconsistent rates + */ +bool pdbm_constrainN(PDBMPtr& pdbm, cindex_t dim, const constraint_t* constraints, size_t n); + +/** + * Constrain a priced DBM to a facet (\a i, \a j). + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param i is a clock index + * @param j is a clock index + * @post The DBM is empty or closed. + * @return + * true if and only if the result is non empty. + */ +bool pdbm_constrainToFacet(PDBMPtr& pdbm, cindex_t dim, cindex_t i, cindex_t j); + +/** + * Relation between two priced DBMs. + * + * @see relation_t + * @param pdbm1 is a closed priced DBMs of dimension \a dim. + * @param pdbm2 is a closed priced DBMs of dimension \a dim. + * @param dim is the dimension of \a pdbm1 and \a pdbm2. + * @return The relation between pdbm1 and pdbm2 + */ +relation_t pdbm_relation(const PDBMCPtr& pdbm1, const PDBMCPtr& pdbm2, cindex_t dim); + +/** + * Relation between 2 priced dbms where one is in compressed. Notice + * that in constract to dbm_relationWithMinDBM, \a buffer may not be + * NULL. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param minDBM is the compressed priced dbm. + * @param buffer is a buffer of size dim*dim. + * @return The relation between pdbm1 and pdbm2 + * @see dbm_relationWithMinDBM + * @see relation_t + */ +relation_t pdbm_relationWithMinDBM(const PDBMCPtr& pdbm, cindex_t dim, const mingraph_t minDBM, raw_t* buffer); + +/** + * Computes the infimum cost of the priced DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim + * @param dim is the dimension of \a pdbm. + * @return The infimum cost of \a pdbm. + */ +int32_t pdbm_getInfimum(const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Generates a valuation which has the infimum cost of the priced DBM. + * + * There is no guarantee that the valuation will be contained in the + * priced DBM, but if it is not it is arbitrarily close to a valuation + * that is contained in the priced DBM. + * + * A \a false entry in \a free indicates that the value of this clock + * has already been set in \a valuation and that this clock must not + * be modified. The function will only modify free clocks. + * + * @param pdbm is a closed priced DBM of dimension \a dim + * @param dim is the dimension of \a pdbm. + * @param valuation is an array of at least \a dim elements to which the + * the valuation will be written. + * @param free is an array of at least \a dim elements. + * + * @return The cost of \a valuation. + * + * @throw out_of_range if no valuation with the given constraints can + * be found. + */ +int32_t pdbm_getInfimumValuation(const PDBMCPtr& pdbm, cindex_t dim, int32_t* valuation, const bool* free); + +/** + * Check if a priced DBM satisfies a given constraint. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param i,j are theindices of clocks for the clock constraint. + * @param constraint is the raw_t bound. + * @return true if the DBM satisfies the constraint. + */ +bool pdbm_satisfies(const PDBMCPtr& pdbm, cindex_t dim, cindex_t i, cindex_t j, raw_t constraint); + +/** + * Returns true if the priced DBM is empty. + * + * @param pdbm is a closed or empty priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @return true if and only if the priced dbm is empty. + */ +bool pdbm_isEmpty(const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Check if at least one point can delay infinitely. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @return true if unbounded, false otherwise. + */ +bool pdbm_isUnbounded(const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Compute a hash value for a priced DBM. + * + * ISSUE: canonical form needed. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param seed is a seed for the hash function. + * @return hash value. + */ +uint32_t pdbm_hash(const PDBMCPtr& pdbm, cindex_t dim, uint32_t seed); + +/** + * Test if a point is included in the priced DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param pt is a clock valuation. + * @return true if \a pt satisfies the constraints of dbm. + */ +bool pdbm_containsInt(const PDBMCPtr& pdbm, cindex_t dim, const int32_t* pt); + +/** + * Test if a point is included in the priced DBM when strictness of + * constraints are ignored. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param pt is a clock valuation. + * @return true if \a pt satisfies the constraints of dbm + * (ignoring strictness) + */ +bool pdbm_containsIntWeakly(const PDBMCPtr& pdbm, cindex_t dim, const int32_t* pt); + +/** + * Test if a point is included in the priced DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param pt is a clock valuation. + * @return true if \a pt satisfies the constraints of dbm. + */ +bool pdbm_containsDouble(const PDBMCPtr& pdbm, cindex_t dim, const double* pt); + +/** + * Delay with the current delay rate. + * + * @param pdbm is a closed priced dbm of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @see pdbm_delayRate() + * @post The priced DBM is closed. + */ +void pdbm_up(PDBMPtr& pdbm, cindex_t dim); + +/** + * Delay with delay rate \a rate. There must be at least one clock + * forming a zero cycle with the reference clock. + * + * @param pdbm is a closed priced dbm of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param rate is the cost of delaying. + * @param zero is the index of a clock forming a zero cycle with the + * reference clock. + * @post The priced DBM is closed. + */ +void pdbm_upZero(PDBMPtr& pdbm, cindex_t dim, int32_t rate, cindex_t zero); + +/** + * Updates \a clock to \a value. This is only legitimate if the + * current rate of \a clock is zero. + * + * @param pdbm is a closed priced dbm of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param clock is the index of the clock to reset. + * @param value is the value to which to set \a clock. + * @pre pdbm_getRate(pdbm, dim, clock) == 0 + * @post The priced DBM is closed. + */ +void pdbm_updateValue(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, uint32_t value); + +/** + * Updates \a clock to \a value. This is only legitimate if the clock + * forms a zero cycle with another clock (possibly the reference + * clock). + * + * @param pdbm is a closed priced dbm of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param clock is the index of the clock to reset. + * @param value is the value to which to set \a clock. + * @param zero is a clock (possibly zero) forming a zero cycle + * with \a clock. + * @post The priced DBM is closed. + * @post pdbm_getRate(pdbm, dim, clock) == 0 + */ +void pdbm_updateValueZero(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, uint32_t value, cindex_t zero); + +/** + * Unfinished extrapolation function. + * + * @see dbm_extrapolateMaxBounds + */ +void pdbm_extrapolateMaxBounds(PDBMPtr& pdbm, cindex_t dim, int32_t* max); + +/** + * Unfinished extrapolation function. + * + * @see dbm_diagonalExtrapolateMaxBounds + */ +void pdbm_diagonalExtrapolateMaxBounds(PDBMPtr& pdbm, cindex_t dim, int32_t* max); + +/** + * Extrapolate a priced zone. The extrapolation is based on a lower + * and an upper bound for each clock. The output zone simulates the + * input zone. + * + * The implementation is not finished. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param lower is an array of lower bounds for each clock. + * @param upper is an array of upper bounds for each clock. + * @pre lower[0] = upper[0] = 0 + * @post The priced DBM is closed. + * @see dbm_diagonalExtrapolateLUBounds + */ +void pdbm_diagonalExtrapolateLUBounds(PDBMPtr& pdbm, cindex_t dim, int32_t* lower, int32_t* upper); + +/** + * Increments the cost of each point in a priced DBM by \a value. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param value is the amount with which to increase the cost. + * @post The priced DBM is closed. + * @pre value >= 0 + */ +void pdbm_incrementCost(PDBMPtr& pdbm, cindex_t dim, int32_t value); + +/** + * Compute the closure of a priced DBM. This function is only relevant + * if the DBM has been modified with \c pdbm_setBound(). + * + * @param pdbm is a priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @post The priced DBM is closed or empty. + * @see pdbm_setBound + */ +void pdbm_close(PDBMPtr& pdbm, cindex_t dim); + +/** + * Analyze a priced DBM for its minimal graph representation. Computes + * the smallest number of constraints needed to represent the same + * zone. + * + * The result is returned as a bit matrix, where each bit indicates + * whether that entry of the DBM is relevant. I.e. if the bit \f$ i + * \cdot dim + j\f$ is set, then the constraint \f$(i,j)\f$ of \a dbm + * is needed. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param bitMatrix is bit matrix of size dim*dim + * @return The number of bits marked one in \a bitMatrix. + */ +size_t pdbm_analyzeForMinDBM(const PDBMCPtr& pdbm, cindex_t dim, uint32_t* bitMatrix); + +/** + * Convert the DBM to a more compact representation. + * + * The API supports allocation of larger data structures than needed + * for the actual zone representation. When the \a offset argument is + * bigger than zero, \a offset extra integers are allocated and the + * zone is written with the given offset. Thus when \c + * int32_t[data_size] is needed to represent the reduced zone, an \c + * int32_t array of size \c offset+data_size is allocated. The first + * \a offset elements can be used by the caller. It is important to + * notice that the other functions typically expect a pointer to the + * actual zone data and not to the beginning of the allocated + * block. Thus in the following piece of code, most functions expect + * \c mg and not \c memory: + * + * \code + * int32_t *memory = dbm_writeToMinDBMWithOffset(...); + * mingraph_t mg = &memory[offset]; + * \endcode + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param minimizeGraph when true, try to use minimal constraint form. + * @param tryConstraints16 when true, try to save constraints on 16 bits. + * @param c_alloc is a C allocator wrapper. + * @param offset is the offset for allocation. + * @return The converted priced DBM. The first \a offset integers are unused. + * @pre allocFunction allocates memory in integer units + */ +int32_t* pdbm_writeToMinDBMWithOffset(const PDBMCPtr& pdbm, cindex_t dim, bool minimizeGraph, bool tryConstraints16, + allocator_t c_alloc, uint32_t offset); + +/** + * Uncompresses a compressed priced DBM. The compressed priced DBM \a + * src is written to \a dst. The destination does not need to + * initialised beforehand. The dimension \a dim must match the + * dimension of the compressed priced DBM. + * + * @param dst is a memory area of size pdbm_size(dim). + * @param dim is the dimension of the priced DBM. + * @param src is the compressed priced DBM. + * @post dst is a closed priced DBM. + */ +void pdbm_readFromMinDBM(PDBMPtr& dst, cindex_t dim, mingraph_t src); + +/** + * Finds a clock that is on a zero cycle with \a clock. Returns true + * if and only if such a clock is found. If multiple clocks are on a + * zero cycle with \a clock, then the clock with the smallest index is + * returned. + * + * This function is equivalent to calling \c pdbm_findNextZeroCyle + * with \a output initialised to zero. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param clock is the index of a clock. + * @param out is where the clock found is written. + * @return true if and only if a clock is found. + */ +bool pdbm_findZeroCycle(const PDBMCPtr& pdbm, cindex_t dim, cindex_t clock, cindex_t* out); + +/** + * Finds a clock that is on a zero cycle with \a clock. Returns true + * if and only if such a clock is found. Only clocks with a value + * equal to or greater than the existing value of \a output are + * considered. If multiple clocks are on a zero cycle with \a clock, + * then the clock with the smallest index is returned. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param clock is the index of a clock. + * @param out is where the clock found is written. + * @return true if and only if a clock is found. + */ +bool pdbm_findNextZeroCycle(const PDBMCPtr& pdbm, cindex_t dim, cindex_t x, cindex_t* out); + +/** + * Returns the slope of the cost plane along the delay trajectory. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + */ +int32_t pdbm_getSlopeOfDelayTrajectory(const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Returns the rate (coefficient of the hyperplane) of \a clock. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param clock is the clock for which to return the coefficient. + * @return the rate of \a clock. + */ +int32_t pdbm_getRate(const PDBMCPtr& pdbm, cindex_t dim, cindex_t clock); + +const int32_t* pdbm_getRates(const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Returns the cost of the offset point. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + */ +uint32_t pdbm_getCostAtOffset(const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Sets the cost at the offset point. + * + * Notice that the value must be large enough such that the infimum + * cost of the priced DBM is not negative. Temporary inconsistencies + * are allowed as long as no operations are performed on the priced + * DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param value is the new cost of the offset point. + */ +void pdbm_setCostAtOffset(PDBMPtr& pdbm, cindex_t dim, uint32_t value); + +/** + * Returns true if the DBM is valid. Useful for debugging. + * + * @param pdbm is a priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + */ +bool pdbm_isValid(const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Computes the lower facets of a priced DBM relative to \a clock. As + * a side-effect, all lower facets relative to \a clock are converted + * to weak facets. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param clock is the clock for which to return the lower facets + * @param facets is an array of at least \a dim elements to which + * the lower facets will be written. + * @return The number of facets written to \a facets. + */ +uint32_t pdbm_getLowerRelativeFacets(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, cindex_t* facets); + +/** + * Computes the upper facets of a priced DBM relative to \a clock. As + * a side-effect, all upper facets relative to \a clock are converted + * to weak facets. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param clock is the clock for which to return the upper facets + * @param facets is an array of at least \a dim elements to which + * the upper facets will be written. + * @return The number of facets written to \a facets. + */ +uint32_t pdbm_getUpperRelativeFacets(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, cindex_t* facets); + +/** + * Computes the lower facets of a priced DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param facets is an array of at least \a dim elements to which the + * lower facets will be written. + * @return The number of facets written to \a facets. + */ +uint32_t pdbm_getLowerFacets(PDBMPtr& pdbm, cindex_t dim, cindex_t* facets); + +/** + * Computes the upper facets of a priced DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param facets is an array of at least \a dim elements to which the + * upper facets will be written. + * @return The number of facets written to \a facets. + */ +uint32_t pdbm_getUpperFacets(PDBMPtr& pdbm, cindex_t dim, cindex_t* facets); + +/** + * Computes the cost of a valuation in a priced DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param valuation is a valuation in \a pdbm. + * @pre pdbm_containsInt(pdbm, dim, valuation) + * @return The cost of \a valuation in \a pdbm. + */ +int32_t pdbm_getCostOfValuation(const PDBMCPtr& pdbm, cindex_t dim, const int32_t* valuation); + +/** + * Makes all strong constraints of a priced DBM weak. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @post The priced DBM is closed. + */ +void pdbm_relax(PDBMPtr& pdbm, cindex_t dim); + +/** + * Computes the offset point of a priced DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param valuation is an array of at least \a dim elements to which the + * offset point is written. + */ +void pdbm_getOffset(const PDBMCPtr& pdbm, cindex_t dim, int32_t* valuation); + +/** + * Sets a coefficient of the hyperplane of a priced DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param clock is the index of a clock for which to set the coefficient. + * @param rate is the coefficient. + */ +void pdbm_setRate(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, int32_t rate); + +/** + * Returns the inner matrix of a priced DBM. The matrix can be + * modified as long as \c pdbm_close() is called before any other + * operations are performed on the priced DBM. + * + * @param pdbm is a priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + */ +raw_t* pdbm_getMutableMatrix(PDBMPtr& pdbm, cindex_t dim); + +/** + * Returns the inner matrix of a priced DBM. The matrix is read-only. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + */ +const raw_t* pdbm_getMatrix(const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Frees a clock of a priced DBM. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param clock is the index of the clock to free. + */ +void pdbm_freeClock(PDBMPtr& pdbm, cindex_t dim, cindex_t clock); + +/** + * Prints a priced DBM to a stream. + * + * @param f is the stream to print to. + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @see dbm_print + */ +void pdbm_print(FILE* f, const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Prints a priced DBM to a stream. + * + * @param o is the stream to print to. + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @see dbm_print + */ +std::ostream& pdbm_print(std::ostream& os, const PDBMCPtr& pdbm, cindex_t dim); + +/** + * Implementation of the free up operation for priced DBMs. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param index is an index of a clock. + * @see dbm_freeUp + */ +void pdbm_freeUp(PDBMPtr& pdbm, cindex_t dim, cindex_t index); + +/** + * Implementation of the free down operation for priced DBMs. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + * @param index is an index of a clock. * + * @see dbm_freeDown + */ +void pdbm_freeDown(PDBMPtr& pdbm, cindex_t dim, cindex_t index); + +/** + * Checks whether a priced DBM is in normal form. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + */ +bool pdbm_hasNormalForm(const PDBMPtr& pdbm, cindex_t dim); + +/** + * Brings a priced DBM into normal form. + * + * @param pdbm is a closed priced DBM of dimension \a dim. + * @param dim is the dimension of \a pdbm. + */ +void pdbm_normalise(const PDBMPtr& pdbm, cindex_t dim); + +/////////////////////////////////////////////////////////////////////////// + +#endif /* INCLUDE_DBM_PRICED_H */ diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp new file mode 100644 index 0000000..e69de29 From 8a51c4e1a7aef547b74e942e3dc90acb1f8e97d3 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Thu, 20 Aug 2026 11:15:21 +0200 Subject: [PATCH 03/18] PPDBM data structure first draft --- include/dbm/ParamPricedDBM.h | 783 +++-------------------------------- src/ParamPricedDBM.cpp | 92 ++++ 2 files changed, 144 insertions(+), 731 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index 046724b..854f3b7 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -1,763 +1,84 @@ -// -*- mode: C++; c-file-style: "stroustrup"; c-basic-offset: 4; indent-tabs-mode: nil; -*- -//////////////////////////////////////////////////////////////////// -// -// This file is a part of the UPPAAL toolkit. -// Copyright (c) 1995 - 2005, Uppsala University and Aalborg University. -// All right reserved. -// -// $Id: priced.h,v 1.15 2005/05/31 20:54:59 behrmann Exp $ -// -/////////////////////////////////////////////////////////////////// - #ifndef INCLUDE_DBM_PARAMPRICEDDBM_H #define INCLUDE_DBM_PARAMPRICEDDBM_H -#include "dbm/dbm.h" -#include "dbm/mingraph.h" -#include "base/inttypes.h" +#include "dbm.h" +#include "mingraph.h" +#include "constraints.h" #include #include #include #include +#include -/** - * @file - * - * Functions for handling priced DBMs. - * - * A priced DBM is a data structure representing a zone with an affine - * hyperplane over the zone assigning a cost to each point in the - * zone. - * - * A few facts about priced DBMs: - * - * - Priced DBMs have a positive dimension. - * - * - Priced DBMs are stored continously in memory. - * - * - The clock with index zero is the reference clock. - * - * - The affine hyperplane is given by the cost in the offset point - * and the coefficients of the hyperplane. - * - * - A closed DBM is never empty. - * - * Allocation can be handled by the library by using \c - * pdbm_allocate() and \c pdbm_deallocate(). This can optionally be - * combined with reference counting using \c pdbm_incRef() and \c - * pdbm_decRef(), in which case \c pdbm_deallocate() is called - * automatically as soon as the reference count reaches zero. - * Reference counting is used to implement copy on write - * semantics. Thus whenever you modify a priced DBM with a reference - * count larger than 1, then that DBM is automatically copied. For - * this reason, all functions modifying priced DBMs take a reference - * parameter to the priced DBM to modify. - * - * A NULL pointer is equivalent to an empty DBM with a reference count - * of 1. Functions that can cause the priced DBM to become empty can - * choose to deallocate the DBM rather than to copy it. All functions - * that do not require the input DBM to be closed will allocate a new - * DBM when given a NULL pointer as input. - * - * Alternatively, memory can be allocated outside the library. In that - * case \c pdbm_size() bytes need to be allocated and preinitialized - * for the use by the library by calling \c pdbm_reserve(). Priced - * DBMs allocated in this manner must not be reference counted. - */ - -/** - * Data type for priced dbm. - */ -struct PDBM_s; -using PDBMPtr = std::shared_ptr; -using PDBMCPtr = std::shared_ptr; - -/** - * Computes the size in number of bytes of a priced dbm of the given - * dimension. - * - * @param dim is a dimension. - * @pre dim > 0 - * @return Amount of memory in bytes. - */ -size_t pdbm_size(cindex_t dim); - -/** - * Reserves a memory area for use as a priced DBM. The priced DBM will - * be unitialised, hence further initialised with e.g. \c pdbm_init() - * or \c pdbm_zero() is required. - * - * @param dim is the dimension of the priced DBM to reserve. - * @param p is a memory area of size \c pdbm_size(dim). - * - * @return An unitialised priced DBM of dimension \a dim. - */ -PDBMPtr pdbm_reserve(cindex_t dim, void* p); - -/** - * Allocates a new priced DBM. The reference count is initialised to - * 0. No other initialisation is performed. - * - * @param dim is the dimension. - * @return A newly allocated priced DBM of dimension \a dim. - * @pre dim is larger than 0. - */ -PDBMPtr pdbm_allocate(cindex_t dim); -/** - * Deallocates a priced DBM. - * - * @pre - * - \a pdbm was allocated with \c pdbm_allocate(). - * - The reference count of \a pdbm is zero. - */ -void pdbm_deallocate(PDBM_s* pdbm); - -/** - * Copy a priced DBM. - * - * Priced DBMs normally use a copy-on-write scheme, thus is no need to - * call this function. The only use of it (besides internally in the - * library as part of the copy-on-write implementation) is for priced - * DBMs, for which you choose not to use reference counting. Hence, a - * call to this function is only permissible if the reference count of - * \a dst is zero. If \a dst is NULL, a new DBM is allocate with \c - * pdbm_allocate(). - * - * @param dst The destination. - * @param src The source. - * @param dim The dimension of \a dst and \a src. - * @pre The reference count of \a dst is zero or dst is NULL. - * @post The reference count of the return value is zero. - * @return The destination. - */ -PDBMPtr pdbm_copy(PDBMPtr dst, const PDBMCPtr& src, cindex_t dim); -/** - * Initialises a priced DBM to the DBM containing all valuations. The - * cost of all valuations will be zero. If \a pdbm is NULL, a new DBM - * is allocated with \c pdbm_allocate(). - * - * @param pdbm is the priced DBM to initialize. - * @param dim is the dimension of \a pdbm. - */ -void pdbm_init(PDBMPtr& pdbm, cindex_t dim); -/** - * Initialize a priced DBM to only contain the origin with a cost of - * 0. If \a pdbm is NULL, a new DBM is allocated with \c - * pdbm_allocate(). - * - * @param pdbm is the priced DBM to initialise. - * @param dim is the dimension of \a pdbm. - */ -void pdbm_zero(PDBMPtr& pdbm, cindex_t dim); +constexpr int INF = INT_MAX >> 1; -/** - * Constrain a priced DBM. - * - * After the call, the DBM will be constrained such that the - * difference between clock \a i and clock \a j is smaller than \a - * constraint. - * - * @see dbm_constrain1 - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param i is a clock index. - * @param j is a clock index. - * @param constraint is a raw_t bound. - * @pre i != j - * @post The DBM is empty or closed. - * @return True if and only if the result is not empty. - */ -bool pdbm_constrain1(PDBMPtr& pdbm, cindex_t dim, cindex_t i, cindex_t j, raw_t constraint); /** - * Constrain a priced DBM with multiple constraints. - * - * @see dbm_constraint1 - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param constraints is an array of \a n constraints. - * @param n is the length of \a constraints. - * @pre the constraints are valid - * @post The DBM is empty or closed. - * @return - * - true if the DBM is non empty + - * closed non empty DBM + updated consistent rates - * - or false if the DBM is empty + - * empty DBM + inconsistent rates + * A parameter constraint is of the form (a_0 + sum^n_{i=1} a_i p_i) <= 0 */ -bool pdbm_constrainN(PDBMPtr& pdbm, cindex_t dim, const constraint_t* constraints, size_t n); +class ParameterConstraint { public: std::vector coeffs; }; -/** - * Constrain a priced DBM to a facet (\a i, \a j). - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param i is a clock index - * @param j is a clock index - * @post The DBM is empty or closed. - * @return - * true if and only if the result is non empty. - */ -bool pdbm_constrainToFacet(PDBMPtr& pdbm, cindex_t dim, cindex_t i, cindex_t j); +class Polyhedron { public: std::vector constraints; }; /** - * Relation between two priced DBMs. - * - * @see relation_t - * @param pdbm1 is a closed priced DBMs of dimension \a dim. - * @param pdbm2 is a closed priced DBMs of dimension \a dim. - * @param dim is the dimension of \a pdbm1 and \a pdbm2. - * @return The relation between pdbm1 and pdbm2 + * a constraint stored in a DBM is of the form x_i - x_j ~ value with ~ being < or <= */ -relation_t pdbm_relation(const PDBMCPtr& pdbm1, const PDBMCPtr& pdbm2, cindex_t dim); +class DbmBound +{ +public: + int value; // constant coeff of the constraint + bool strict; // is the constraint strict? -/** - * Relation between 2 priced dbms where one is in compressed. Notice - * that in constract to dbm_relationWithMinDBM, \a buffer may not be - * NULL. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param minDBM is the compressed priced dbm. - * @param buffer is a buffer of size dim*dim. - * @return The relation between pdbm1 and pdbm2 - * @see dbm_relationWithMinDBM - * @see relation_t - */ -relation_t pdbm_relationWithMinDBM(const PDBMCPtr& pdbm, cindex_t dim, const mingraph_t minDBM, raw_t* buffer); + constexpr DbmBound(int value, bool strict): value(value), strict(strict){} + bool isInfinite() const { return value == INF; } // is this constraint unrestraining? + bool isLooserThan(const DbmBound& other) const { + return value > other.value || (value == other.value && strict && !other.strict); + } + DbmBound negate() const { assert(!isInfinite()); return DbmBound(-value, !strict); } +}; +constexpr DbmBound INF_BOUND(INF,true); +constexpr DbmBound DIAG_BOUND(0,false); -/** - * Computes the infimum cost of the priced DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim - * @param dim is the dimension of \a pdbm. - * @return The infimum cost of \a pdbm. - */ -int32_t pdbm_getInfimum(const PDBMCPtr& pdbm, cindex_t dim); +class DBMatrix +{ +public: + DBMatrix(int dim): dim(dim), data(dim * dim, INF_BOUND) + { + for (int i = 0; i < dim; ++i) { + data[i * dim + i] = DIAG_BOUND; + }; + }; -/** - * Generates a valuation which has the infimum cost of the priced DBM. - * - * There is no guarantee that the valuation will be contained in the - * priced DBM, but if it is not it is arbitrarily close to a valuation - * that is contained in the priced DBM. - * - * A \a false entry in \a free indicates that the value of this clock - * has already been set in \a valuation and that this clock must not - * be modified. The function will only modify free clocks. - * - * @param pdbm is a closed priced DBM of dimension \a dim - * @param dim is the dimension of \a pdbm. - * @param valuation is an array of at least \a dim elements to which the - * the valuation will be written. - * @param free is an array of at least \a dim elements. - * - * @return The cost of \a valuation. - * - * @throw out_of_range if no valuation with the given constraints can - * be found. - */ -int32_t pdbm_getInfimumValuation(const PDBMCPtr& pdbm, cindex_t dim, int32_t* valuation, const bool* free); + DbmBound& operator()(int i, int j){ return data[i * dim + j]; }; + const DbmBound& operator()(int i, int j) const { return data[i * dim + j]; }; -/** - * Check if a priced DBM satisfies a given constraint. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param i,j are theindices of clocks for the clock constraint. - * @param constraint is the raw_t bound. - * @return true if the DBM satisfies the constraint. - */ -bool pdbm_satisfies(const PDBMCPtr& pdbm, cindex_t dim, cindex_t i, cindex_t j, raw_t constraint); + int dim; // number of clocks + 1 (for the ref clock) + std::vector data; +}; -/** - * Returns true if the priced DBM is empty. - * - * @param pdbm is a closed or empty priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @return true if and only if the priced dbm is empty. - */ -bool pdbm_isEmpty(const PDBMCPtr& pdbm, cindex_t dim); /** - * Check if at least one point can delay infinitely. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @return true if unbounded, false otherwise. + * Parametric priced timed difference bound matrix */ -bool pdbm_isUnbounded(const PDBMCPtr& pdbm, cindex_t dim); +class Ppdbm +{ + int dim, param = 0; // number of clocks and parameters + DBMatrix DBM; // the classic zone + bool emptyZone = false; // is the zone represented by the dbm empty? -/** - * Compute a hash value for a priced DBM. - * - * ISSUE: canonical form needed. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param seed is a seed for the hash function. - * @return hash value. - */ -uint32_t pdbm_hash(const PDBMCPtr& pdbm, cindex_t dim, uint32_t seed); - -/** - * Test if a point is included in the priced DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param pt is a clock valuation. - * @return true if \a pt satisfies the constraints of dbm. - */ -bool pdbm_containsInt(const PDBMCPtr& pdbm, cindex_t dim, const int32_t* pt); - -/** - * Test if a point is included in the priced DBM when strictness of - * constraints are ignored. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param pt is a clock valuation. - * @return true if \a pt satisfies the constraints of dbm - * (ignoring strictness) - */ -bool pdbm_containsIntWeakly(const PDBMCPtr& pdbm, cindex_t dim, const int32_t* pt); - -/** - * Test if a point is included in the priced DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param pt is a clock valuation. - * @return true if \a pt satisfies the constraints of dbm. - */ -bool pdbm_containsDouble(const PDBMCPtr& pdbm, cindex_t dim, const double* pt); - -/** - * Delay with the current delay rate. - * - * @param pdbm is a closed priced dbm of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @see pdbm_delayRate() - * @post The priced DBM is closed. - */ -void pdbm_up(PDBMPtr& pdbm, cindex_t dim); - -/** - * Delay with delay rate \a rate. There must be at least one clock - * forming a zero cycle with the reference clock. - * - * @param pdbm is a closed priced dbm of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param rate is the cost of delaying. - * @param zero is the index of a clock forming a zero cycle with the - * reference clock. - * @post The priced DBM is closed. - */ -void pdbm_upZero(PDBMPtr& pdbm, cindex_t dim, int32_t rate, cindex_t zero); - -/** - * Updates \a clock to \a value. This is only legitimate if the - * current rate of \a clock is zero. - * - * @param pdbm is a closed priced dbm of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param clock is the index of the clock to reset. - * @param value is the value to which to set \a clock. - * @pre pdbm_getRate(pdbm, dim, clock) == 0 - * @post The priced DBM is closed. - */ -void pdbm_updateValue(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, uint32_t value); - -/** - * Updates \a clock to \a value. This is only legitimate if the clock - * forms a zero cycle with another clock (possibly the reference - * clock). - * - * @param pdbm is a closed priced dbm of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param clock is the index of the clock to reset. - * @param value is the value to which to set \a clock. - * @param zero is a clock (possibly zero) forming a zero cycle - * with \a clock. - * @post The priced DBM is closed. - * @post pdbm_getRate(pdbm, dim, clock) == 0 - */ -void pdbm_updateValueZero(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, uint32_t value, cindex_t zero); - -/** - * Unfinished extrapolation function. - * - * @see dbm_extrapolateMaxBounds - */ -void pdbm_extrapolateMaxBounds(PDBMPtr& pdbm, cindex_t dim, int32_t* max); - -/** - * Unfinished extrapolation function. - * - * @see dbm_diagonalExtrapolateMaxBounds - */ -void pdbm_diagonalExtrapolateMaxBounds(PDBMPtr& pdbm, cindex_t dim, int32_t* max); - -/** - * Extrapolate a priced zone. The extrapolation is based on a lower - * and an upper bound for each clock. The output zone simulates the - * input zone. - * - * The implementation is not finished. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param lower is an array of lower bounds for each clock. - * @param upper is an array of upper bounds for each clock. - * @pre lower[0] = upper[0] = 0 - * @post The priced DBM is closed. - * @see dbm_diagonalExtrapolateLUBounds - */ -void pdbm_diagonalExtrapolateLUBounds(PDBMPtr& pdbm, cindex_t dim, int32_t* lower, int32_t* upper); - -/** - * Increments the cost of each point in a priced DBM by \a value. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param value is the amount with which to increase the cost. - * @post The priced DBM is closed. - * @pre value >= 0 - */ -void pdbm_incrementCost(PDBMPtr& pdbm, cindex_t dim, int32_t value); - -/** - * Compute the closure of a priced DBM. This function is only relevant - * if the DBM has been modified with \c pdbm_setBound(). - * - * @param pdbm is a priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @post The priced DBM is closed or empty. - * @see pdbm_setBound - */ -void pdbm_close(PDBMPtr& pdbm, cindex_t dim); - -/** - * Analyze a priced DBM for its minimal graph representation. Computes - * the smallest number of constraints needed to represent the same - * zone. - * - * The result is returned as a bit matrix, where each bit indicates - * whether that entry of the DBM is relevant. I.e. if the bit \f$ i - * \cdot dim + j\f$ is set, then the constraint \f$(i,j)\f$ of \a dbm - * is needed. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param bitMatrix is bit matrix of size dim*dim - * @return The number of bits marked one in \a bitMatrix. - */ -size_t pdbm_analyzeForMinDBM(const PDBMCPtr& pdbm, cindex_t dim, uint32_t* bitMatrix); + std::vector offsetCost; // the cost of the offset (affine function of parameters with int coeffs) + std::vector> rates; // the cost rates of clocks (affine functions of parameters with int coeffs) + Polyhedron PC; // parametric constraints set -/** - * Convert the DBM to a more compact representation. - * - * The API supports allocation of larger data structures than needed - * for the actual zone representation. When the \a offset argument is - * bigger than zero, \a offset extra integers are allocated and the - * zone is written with the given offset. Thus when \c - * int32_t[data_size] is needed to represent the reduced zone, an \c - * int32_t array of size \c offset+data_size is allocated. The first - * \a offset elements can be used by the caller. It is important to - * notice that the other functions typically expect a pointer to the - * actual zone data and not to the beginning of the allocated - * block. Thus in the following piece of code, most functions expect - * \c mg and not \c memory: - * - * \code - * int32_t *memory = dbm_writeToMinDBMWithOffset(...); - * mingraph_t mg = &memory[offset]; - * \endcode - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param minimizeGraph when true, try to use minimal constraint form. - * @param tryConstraints16 when true, try to save constraints on 16 bits. - * @param c_alloc is a C allocator wrapper. - * @param offset is the offset for allocation. - * @return The converted priced DBM. The first \a offset integers are unused. - * @pre allocFunction allocates memory in integer units - */ -int32_t* pdbm_writeToMinDBMWithOffset(const PDBMCPtr& pdbm, cindex_t dim, bool minimizeGraph, bool tryConstraints16, - allocator_t c_alloc, uint32_t offset); - -/** - * Uncompresses a compressed priced DBM. The compressed priced DBM \a - * src is written to \a dst. The destination does not need to - * initialised beforehand. The dimension \a dim must match the - * dimension of the compressed priced DBM. - * - * @param dst is a memory area of size pdbm_size(dim). - * @param dim is the dimension of the priced DBM. - * @param src is the compressed priced DBM. - * @post dst is a closed priced DBM. - */ -void pdbm_readFromMinDBM(PDBMPtr& dst, cindex_t dim, mingraph_t src); - -/** - * Finds a clock that is on a zero cycle with \a clock. Returns true - * if and only if such a clock is found. If multiple clocks are on a - * zero cycle with \a clock, then the clock with the smallest index is - * returned. - * - * This function is equivalent to calling \c pdbm_findNextZeroCyle - * with \a output initialised to zero. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param clock is the index of a clock. - * @param out is where the clock found is written. - * @return true if and only if a clock is found. - */ -bool pdbm_findZeroCycle(const PDBMCPtr& pdbm, cindex_t dim, cindex_t clock, cindex_t* out); - -/** - * Finds a clock that is on a zero cycle with \a clock. Returns true - * if and only if such a clock is found. Only clocks with a value - * equal to or greater than the existing value of \a output are - * considered. If multiple clocks are on a zero cycle with \a clock, - * then the clock with the smallest index is returned. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param clock is the index of a clock. - * @param out is where the clock found is written. - * @return true if and only if a clock is found. - */ -bool pdbm_findNextZeroCycle(const PDBMCPtr& pdbm, cindex_t dim, cindex_t x, cindex_t* out); - -/** - * Returns the slope of the cost plane along the delay trajectory. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - */ -int32_t pdbm_getSlopeOfDelayTrajectory(const PDBMCPtr& pdbm, cindex_t dim); - -/** - * Returns the rate (coefficient of the hyperplane) of \a clock. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param clock is the clock for which to return the coefficient. - * @return the rate of \a clock. - */ -int32_t pdbm_getRate(const PDBMCPtr& pdbm, cindex_t dim, cindex_t clock); - -const int32_t* pdbm_getRates(const PDBMCPtr& pdbm, cindex_t dim); - -/** - * Returns the cost of the offset point. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - */ -uint32_t pdbm_getCostAtOffset(const PDBMCPtr& pdbm, cindex_t dim); - -/** - * Sets the cost at the offset point. - * - * Notice that the value must be large enough such that the infimum - * cost of the priced DBM is not negative. Temporary inconsistencies - * are allowed as long as no operations are performed on the priced - * DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param value is the new cost of the offset point. - */ -void pdbm_setCostAtOffset(PDBMPtr& pdbm, cindex_t dim, uint32_t value); - -/** - * Returns true if the DBM is valid. Useful for debugging. - * - * @param pdbm is a priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - */ -bool pdbm_isValid(const PDBMCPtr& pdbm, cindex_t dim); - -/** - * Computes the lower facets of a priced DBM relative to \a clock. As - * a side-effect, all lower facets relative to \a clock are converted - * to weak facets. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param clock is the clock for which to return the lower facets - * @param facets is an array of at least \a dim elements to which - * the lower facets will be written. - * @return The number of facets written to \a facets. - */ -uint32_t pdbm_getLowerRelativeFacets(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, cindex_t* facets); - -/** - * Computes the upper facets of a priced DBM relative to \a clock. As - * a side-effect, all upper facets relative to \a clock are converted - * to weak facets. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param clock is the clock for which to return the upper facets - * @param facets is an array of at least \a dim elements to which - * the upper facets will be written. - * @return The number of facets written to \a facets. - */ -uint32_t pdbm_getUpperRelativeFacets(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, cindex_t* facets); - -/** - * Computes the lower facets of a priced DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param facets is an array of at least \a dim elements to which the - * lower facets will be written. - * @return The number of facets written to \a facets. - */ -uint32_t pdbm_getLowerFacets(PDBMPtr& pdbm, cindex_t dim, cindex_t* facets); - -/** - * Computes the upper facets of a priced DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param facets is an array of at least \a dim elements to which the - * upper facets will be written. - * @return The number of facets written to \a facets. - */ -uint32_t pdbm_getUpperFacets(PDBMPtr& pdbm, cindex_t dim, cindex_t* facets); - -/** - * Computes the cost of a valuation in a priced DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param valuation is a valuation in \a pdbm. - * @pre pdbm_containsInt(pdbm, dim, valuation) - * @return The cost of \a valuation in \a pdbm. - */ -int32_t pdbm_getCostOfValuation(const PDBMCPtr& pdbm, cindex_t dim, const int32_t* valuation); - -/** - * Makes all strong constraints of a priced DBM weak. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @post The priced DBM is closed. - */ -void pdbm_relax(PDBMPtr& pdbm, cindex_t dim); - -/** - * Computes the offset point of a priced DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param valuation is an array of at least \a dim elements to which the - * offset point is written. - */ -void pdbm_getOffset(const PDBMCPtr& pdbm, cindex_t dim, int32_t* valuation); - -/** - * Sets a coefficient of the hyperplane of a priced DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param clock is the index of a clock for which to set the coefficient. - * @param rate is the coefficient. - */ -void pdbm_setRate(PDBMPtr& pdbm, cindex_t dim, cindex_t clock, int32_t rate); - -/** - * Returns the inner matrix of a priced DBM. The matrix can be - * modified as long as \c pdbm_close() is called before any other - * operations are performed on the priced DBM. - * - * @param pdbm is a priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - */ -raw_t* pdbm_getMutableMatrix(PDBMPtr& pdbm, cindex_t dim); - -/** - * Returns the inner matrix of a priced DBM. The matrix is read-only. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - */ -const raw_t* pdbm_getMatrix(const PDBMCPtr& pdbm, cindex_t dim); - -/** - * Frees a clock of a priced DBM. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param clock is the index of the clock to free. - */ -void pdbm_freeClock(PDBMPtr& pdbm, cindex_t dim, cindex_t clock); - -/** - * Prints a priced DBM to a stream. - * - * @param f is the stream to print to. - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @see dbm_print - */ -void pdbm_print(FILE* f, const PDBMCPtr& pdbm, cindex_t dim); - -/** - * Prints a priced DBM to a stream. - * - * @param o is the stream to print to. - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @see dbm_print - */ -std::ostream& pdbm_print(std::ostream& os, const PDBMCPtr& pdbm, cindex_t dim); - -/** - * Implementation of the free up operation for priced DBMs. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param index is an index of a clock. - * @see dbm_freeUp - */ -void pdbm_freeUp(PDBMPtr& pdbm, cindex_t dim, cindex_t index); - -/** - * Implementation of the free down operation for priced DBMs. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - * @param index is an index of a clock. * - * @see dbm_freeDown - */ -void pdbm_freeDown(PDBMPtr& pdbm, cindex_t dim, cindex_t index); - -/** - * Checks whether a priced DBM is in normal form. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - */ -bool pdbm_hasNormalForm(const PDBMPtr& pdbm, cindex_t dim); - -/** - * Brings a priced DBM into normal form. - * - * @param pdbm is a closed priced DBM of dimension \a dim. - * @param dim is the dimension of \a pdbm. - */ -void pdbm_normalise(const PDBMPtr& pdbm, cindex_t dim); +public: + Ppdbm(int dim, int param); + void reset(); + bool constrain_DBM(int i, int j, DbmBound constraint); +}; -/////////////////////////////////////////////////////////////////////////// -#endif /* INCLUDE_DBM_PRICED_H */ +#endif diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index e69de29..3a9ffe2 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -0,0 +1,92 @@ +/* -*- mode: C++; c-file-style: "stroustrup"; c-basic-offset: 4; indent-tabs-mode: nil; -*- */ + +#include "../include/dbm/ParamPricedDBM.h" + +#include + +#include // transform +#include +#include +#include +#include +#include + +/** + * constructs an unrestrained matrix + * @param dim + * @param param + */ +Ppdbm::Ppdbm(int dim, int param) + : dim(dim), + param(param), + DBM(dim), + offsetCost(param + 1), + rates(dim, std::vector(param + 1)) +{} + +/** + * resets a DBM to its unrestrained state. + */ +void Ppdbm::reset() +{ + std::fill(DBM.data.begin(), DBM.data.end(), INF_BOUND); + for (int i = 0; i < dim; ++i) { + DBM(i,i)= DIAG_BOUND; + } + + std::fill(offsetCost.begin(), offsetCost.end(), 0); + + for (auto& rate : rates) { + std::fill(rate.begin(), rate.end(), 0); + } + + PC.constraints.clear(); +} + +/** + * + * @param i + * @param j + * @param constraint + * @return + */ +bool Ppdbm::constrain_DBM(int i, int j, DbmBound constraint) +{ + assert(i >= 0 && i < dim); + assert(j >= 0 && j < dim); + + // The constraint is useless? + if (constraint.isLooserThan(DBM(i,j))) { + return true; + } + + // The constraint empties the zone? + if (constraint.negate().isLooserThan(DBM(j, i))) { + DBM(i, j) = constraint; + emptyZone = true; + return false; + } + + // remove the cost at the previous offset + for (int h = 1; h < dim; h++) { // for every clock, + for (int k = 0; k <= param; k++) { // for every parameter, + offsetCost[k] -= rates[h][k] * (-DBM(0, h).value); // x_0 - x_k <= value so x_k >= -value + } + } + + // Add the constraint and close the DBM + DBM(i, j) = constraint; + DBM.close(i,j); + + // Add the cost at the new offset + for (int h = 1; h < dim; h++) { // for every clock, + for (int k = 0; k <= param; k++) { // for every parameter, + offsetCost[k] += rates[h][k] * (-DBM(0, h).value); // x_0 - x_k <= value so x_k >= -value + } + } + + return true; +} + + + From 1a0a21ec42f96141cf1bb3145f1b33e3a5e83e4d Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Thu, 20 Aug 2026 15:33:05 +0200 Subject: [PATCH 04/18] before first review --- src/CMakeLists.txt | 2 +- src/ParamPricedDBM.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8d0be52..d727c27 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' diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index 3a9ffe2..d23e449 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -1,6 +1,6 @@ /* -*- mode: C++; c-file-style: "stroustrup"; c-basic-offset: 4; indent-tabs-mode: nil; -*- */ -#include "../include/dbm/ParamPricedDBM.h" +#include "dbm/ParamPricedDBM.h" #include From 05625c46b3ed198b61f002f1b395e2f7b4a34f4f Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Thu, 20 Aug 2026 19:53:29 +0200 Subject: [PATCH 05/18] review taken into account, and some logic fixed in constrain. --- include/dbm/ParamPricedDBM.h | 113 +++++++++++++++++++++++------------ src/ParamPricedDBM.cpp | 63 +++++++------------ 2 files changed, 94 insertions(+), 82 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index 854f3b7..4b37d4a 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -9,75 +9,110 @@ #include #include #include +#include #include +#include +#include +constexpr auto INF = std::numeric_limits::max() >> 1u; -constexpr int INF = INT_MAX >> 1; - +template +struct Matrix +{ + constexpr Matrix(uint32_t rows, uint32_t cols, const T x): _rows{rows}, _cols{cols}, _data(rows * cols, x){} + constexpr auto getRows() const { return _rows; } + constexpr auto getCols() const { return _cols; } + constexpr auto& getData() const { return _data; } + constexpr auto& getData() { return _data; } + constexpr T& operator()(int i, int j) { return _data[i * _cols + j]; } + constexpr const T& operator()(int i, int j) const { return _data[i * _cols + j]; } + constexpr 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 */ -class ParameterConstraint { public: std::vector coeffs; }; +struct ParameterConstraint { std::vector coeffs; }; -class Polyhedron { public: std::vector constraints; }; +struct Polyhedron { std::vector constraints; }; /** * a constraint stored in a DBM is of the form x_i - x_j ~ value with ~ being < or <= */ -class DbmBound +struct DbmBound { -public: - int value; // constant coeff of the constraint - bool strict; // is the constraint strict? - - constexpr DbmBound(int value, bool strict): value(value), strict(strict){} - bool isInfinite() const { return value == INF; } // is this constraint unrestraining? - bool isLooserThan(const DbmBound& other) const { - return value > other.value || (value == other.value && strict && !other.strict); - } - DbmBound negate() const { assert(!isInfinite()); return DbmBound(-value, !strict); } + constexpr auto getValue() const { return _value; } + constexpr auto getStrict() const { return _strict; } + + constexpr DbmBound( int32_t value, bool strict): _value{value}, _strict{strict}{} + + constexpr bool isInfinite() const { return _value == INF; } ///< is this constraint unrestraining? + constexpr bool dominates(const DbmBound& other) const { + return _value < other._value || (_value == other._value && (_strict || !other._strict)); } ///< is this constraint useless? + 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())); + } ///< do the two constraints force the zone to be empty? + +private: + int32_t _value : 31; + bool _strict : 1; // constant coeff of the constraint, is the constraint strict? }; -constexpr DbmBound INF_BOUND(INF,true); -constexpr DbmBound DIAG_BOUND(0,false); -class DBMatrix +constexpr auto INF_BOUND = DbmBound{INF,true}; +constexpr auto DIAG_BOUND = DbmBound{0,false}; + +/** + * + */ +struct DBMatrix { -public: - DBMatrix(int dim): dim(dim), data(dim * dim, INF_BOUND) + constexpr DBMatrix(uint32_t dim) : _matrix{dim, dim, INF_BOUND} { - for (int i = 0; i < dim; ++i) { - data[i * dim + i] = DIAG_BOUND; - }; - }; + for (uint32_t i = 0; i < dim; ++i) { _matrix(i, i) = DIAG_BOUND; } + } - DbmBound& operator()(int i, int j){ return data[i * dim + j]; }; - const DbmBound& operator()(int i, int j) const { return data[i * dim + j]; }; + constexpr auto getDim() const { return _matrix.getRows(); } + constexpr DbmBound& operator()(int i, int j) { return _matrix(i, j); } + constexpr const DbmBound& operator()(int i, int j) const { return _matrix(i, j); } + void reset() + { + _matrix.fill(INF_BOUND); + for (uint32_t i = 0; i < _matrix.getRows(); ++i) { _matrix(i,i)= DIAG_BOUND; } + } - int dim; // number of clocks + 1 (for the ref clock) - std::vector data; +private: + Matrix _matrix; }; /** * Parametric priced timed difference bound matrix */ -class Ppdbm +struct Ppdbm { - int dim, param = 0; // number of clocks and parameters - DBMatrix DBM; // the classic zone - bool emptyZone = false; // is the zone represented by the dbm empty? + Ppdbm(uint32_t dim, uint32_t param): _dim{dim}, _param{param}, _dbm{dim}, + _offsetCost(param + 1), + _rates(dim, param + 1, 0) {} - std::vector offsetCost; // the cost of the offset (affine function of parameters with int coeffs) - std::vector> rates; // the cost rates of clocks (affine functions of parameters with int coeffs) - Polyhedron PC; // parametric constraints set - -public: - Ppdbm(int dim, int param); void reset(); - bool constrain_DBM(int i, int j, DbmBound constraint); + bool constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint); + +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 }; diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index d23e449..e31ff17 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -4,43 +4,23 @@ #include -#include // transform +#include #include #include #include #include #include -/** - * constructs an unrestrained matrix - * @param dim - * @param param - */ -Ppdbm::Ppdbm(int dim, int param) - : dim(dim), - param(param), - DBM(dim), - offsetCost(param + 1), - rates(dim, std::vector(param + 1)) -{} - /** * resets a DBM to its unrestrained state. */ void Ppdbm::reset() { - std::fill(DBM.data.begin(), DBM.data.end(), INF_BOUND); - for (int i = 0; i < dim; ++i) { - DBM(i,i)= DIAG_BOUND; - } - - std::fill(offsetCost.begin(), offsetCost.end(), 0); - - for (auto& rate : rates) { - std::fill(rate.begin(), rate.end(), 0); - } - - PC.constraints.clear(); + _dbm.reset(); + std::fill(_offsetCost.begin(), _offsetCost.end(), 0); + _rates.fill(0); + _PC.constraints.clear(); + _emptyZone = false; } /** @@ -50,38 +30,35 @@ void Ppdbm::reset() * @param constraint * @return */ -bool Ppdbm::constrain_DBM(int i, int j, DbmBound constraint) +bool Ppdbm::constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint) { - assert(i >= 0 && i < dim); - assert(j >= 0 && j < dim); + assert(i < _dim && j < _dim); // The constraint is useless? - if (constraint.isLooserThan(DBM(i,j))) { - return true; - } + if (_dbm(i,j).dominates(constraint)) { return true; } // The constraint empties the zone? - if (constraint.negate().isLooserThan(DBM(j, i))) { - DBM(i, j) = constraint; - emptyZone = true; + if (DbmBound::isIncompatible(constraint,_dbm(j,i))) { + _dbm(i, j) = constraint; + _emptyZone = true; return false; } // remove the cost at the previous offset - for (int h = 1; h < dim; h++) { // for every clock, - for (int k = 0; k <= param; k++) { // for every parameter, - offsetCost[k] -= rates[h][k] * (-DBM(0, h).value); // x_0 - x_k <= value so x_k >= -value + for (uint32_t h = 1; h < _dim; h++) { // for every clock, + for (uint32_t k = 0; k <= _param; k++) { // for every parameter, + _offsetCost[k] -= _rates(h,k) * (-_dbm(0, h).getValue()); // x_0 - x_k <= value so x_k >= -value } } // Add the constraint and close the DBM - DBM(i, j) = constraint; - DBM.close(i,j); + _dbm(i, j) = constraint; + //_dbm.close(i,j); // Add the cost at the new offset - for (int h = 1; h < dim; h++) { // for every clock, - for (int k = 0; k <= param; k++) { // for every parameter, - offsetCost[k] += rates[h][k] * (-DBM(0, h).value); // x_0 - x_k <= value so x_k >= -value + for (uint32_t h = 1; h < _dim; h++) { // for every clock, + for (uint32_t k = 0; k <= _param; k++) { // for every parameter, + _offsetCost[k] += _rates(h,k) * (-_dbm(0, h).getValue()); // x_0 - x_k <= value so x_k >= -value } } From e78a2804c93533a1f03e9232cccd5980d699fa61 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Fri, 21 Aug 2026 12:37:17 +0200 Subject: [PATCH 06/18] constraint and constraintN applied to PPDBM working. --- include/dbm/ParamPricedDBM.h | 212 ++++++++++++++++++++--------------- src/ParamPricedDBM.cpp | 193 ++++++++++++++++++++++++------- 2 files changed, 276 insertions(+), 129 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index 4b37d4a..e4a4b2c 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -16,104 +16,140 @@ -constexpr auto INF = std::numeric_limits::max() >> 1u; - -template -struct Matrix -{ - constexpr Matrix(uint32_t rows, uint32_t cols, const T x): _rows{rows}, _cols{cols}, _data(rows * cols, x){} - constexpr auto getRows() const { return _rows; } - constexpr auto getCols() const { return _cols; } - constexpr auto& getData() const { return _data; } - constexpr auto& getData() { return _data; } - constexpr T& operator()(int i, int j) { return _data[i * _cols + j]; } - constexpr const T& operator()(int i, int j) const { return _data[i * _cols + j]; } - constexpr 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 ParameterConstraint { std::vector coeffs; }; - -struct Polyhedron { 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}{} - - constexpr bool isInfinite() const { return _value == INF; } ///< is this constraint unrestraining? - constexpr bool dominates(const DbmBound& other) const { - return _value < other._value || (_value == other._value && (_strict || !other._strict)); } ///< is this constraint useless? - 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())); - } ///< do the two constraints force the zone to be empty? - -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 DIAG_BOUND = DbmBound{0,false}; - -/** - * - */ -struct DBMatrix +namespace PPDBM { - constexpr DBMatrix(uint32_t dim) : _matrix{dim, dim, INF_BOUND} + constexpr auto INF = std::numeric_limits::max() >> 1u; + + template + struct Matrix { - for (uint32_t i = 0; i < dim; ++i) { _matrix(i, i) = DIAG_BOUND; } - } + constexpr Matrix(uint32_t rows, uint32_t cols, const T x): _rows{rows}, _cols{cols}, _data(rows * cols, x){} + constexpr auto getRows() const { return _rows; } + constexpr auto getCols() const { return _cols; } + constexpr auto& getData() const { return _data; } + constexpr auto& getData() { return _data; } + constexpr T& operator()(int i, int j) { return _data[i * _cols + j]; } + constexpr const T& operator()(int i, int j) const { return _data[i * _cols + j]; } + constexpr 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 ParameterConstraint { std::vector coeffs; }; + + struct Polyhedron { std::vector constraints; }; + + - constexpr auto getDim() const { return _matrix.getRows(); } - constexpr DbmBound& operator()(int i, int j) { return _matrix(i, j); } - constexpr const DbmBound& operator()(int i, int j) const { return _matrix(i, j); } - void reset() + + /** + * a constraint stored in a DBM is of the form x_i - x_j ~ value with ~ being < or <= + */ + struct DbmBound { - _matrix.fill(INF_BOUND); - for (uint32_t i = 0; i < _matrix.getRows(); ++i) { _matrix(i,i)= DIAG_BOUND; } - } + constexpr auto getValue() const { return _value; } + constexpr auto getStrict() const { return _strict; } -private: - Matrix _matrix; -}; + constexpr DbmBound( int32_t value, bool strict): _value{value}, _strict{strict}{} + constexpr bool isInfinite() const { return _value == INF; } ///< is this constraint unrestraining? -/** - * Parametric priced timed difference bound matrix - */ -struct Ppdbm -{ - Ppdbm(uint32_t dim, uint32_t param): _dim{dim}, _param{param}, _dbm{dim}, - _offsetCost(param + 1), - _rates(dim, param + 1, 0) {} + constexpr bool dominates(const DbmBound& other) const { + return _value < other._value || (_value == other._value && (_strict || !other._strict)); } ///< is the new constraint useless? + + 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())); + } ///< do the two constraints force the zone to be empty? - void reset(); - bool constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint); + private: + int32_t _value : 31; + bool _strict : 1; ///< constant coeff of the constraint, is the constraint strict? + }; -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? + constexpr auto INF_BOUND = DbmBound{INF,true}; + constexpr auto DIAG_BOUND = DbmBound{0,false}; + + 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() }; + } ///< x_i - x_j <= a and x_j - x_k < b gives x_i - x_k < a + b - 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 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) = DIAG_BOUND; } + } + + constexpr auto getDim() const { return _matrix.getRows(); } + constexpr DbmBound& operator()(int i, int j) { return _matrix(i, j); } + constexpr const DbmBound& operator()(int i, int j) const { return _matrix(i, j); } + void reset() + { + _matrix.fill(INF_BOUND); + for (uint32_t i = 0; i < _matrix.getRows(); ++i) { _matrix(i,i)= DIAG_BOUND; } + } + void close(uint32_t i, uint32_t j); + + private: + Matrix _matrix; + }; + + + + /** + * Parametric priced timed difference bound matrix + */ + struct Ppdbm + { + Ppdbm(uint32_t dim, uint32_t param): _dim{dim}, _param{param}, _dbm{dim}, + _offsetCost(param + 1), + _rates(dim, param + 1, 0) {} + + 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; + + 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 + }; +} // namespace PPDBM #endif diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index e31ff17..5384f83 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -10,60 +10,171 @@ #include #include #include +#include -/** - * resets a DBM to its unrestrained state. - */ -void Ppdbm::reset() +namespace PPDBM { - _dbm.reset(); - std::fill(_offsetCost.begin(), _offsetCost.end(), 0); - _rates.fill(0); - _PC.constraints.clear(); - _emptyZone = false; -} - -/** - * - * @param i - * @param j - * @param constraint - * @return - */ -bool Ppdbm::constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint) -{ - assert(i < _dim && j < _dim); + /** + * + * @param i + * @param j + */ + void DBMatrix::close(uint32_t i, uint32_t j) + { + const auto dim = getDim(); + const auto dij = _matrix(i, j); - // The constraint is useless? - if (_dbm(i,j).dominates(constraint)) { return true; } + 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; + } + } + } + } - // The constraint empties the zone? - if (DbmBound::isIncompatible(constraint,_dbm(j,i))) { - _dbm(i, j) = constraint; - _emptyZone = true; - return false; + + /** + * resets a DBM to its unrestrained state. + */ + void Ppdbm::reset() + { + _dbm.reset(); + std::fill(_offsetCost.begin(), _offsetCost.end(), 0); + _rates.fill(0); + _PC.constraints.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()); + } + } } - // remove the cost at the previous offset - for (uint32_t h = 1; h < _dim; h++) { // for every clock, - for (uint32_t k = 0; k <= _param; k++) { // for every parameter, - _offsetCost[k] -= _rates(h,k) * (-_dbm(0, h).getValue()); // x_0 - x_k <= value so x_k >= -value + /** + * + */ + 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()); + } } } - // Add the constraint and close the DBM - _dbm(i, j) = constraint; - //_dbm.close(i,j); + /** + * + * @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); - // Add the cost at the new offset - for (uint32_t h = 1; h < _dim; h++) { // for every clock, - for (uint32_t k = 0; k <= _param; k++) { // for every parameter, - _offsetCost[k] += _rates(h,k) * (-_dbm(0, h).getValue()); // x_0 - x_k <= value so x_k >= -value + // 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; } - return true; -} + /** + * + * @param constraints + * @return + */ + bool Ppdbm::constrain_DBM_N(const std::vector& constraints) + { + assert(!constraints.empty()); + + // 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; + } +} // namespace PPBDM From f1e55add64a9462dc20bfb295d7315860adf507b Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Fri, 21 Aug 2026 15:33:31 +0200 Subject: [PATCH 07/18] first draft of relation function --- include/dbm/ParamPricedDBM.h | 58 +++++++++---- src/ParamPricedDBM.cpp | 153 +++++++++++++++++++++++++++++++++-- 2 files changed, 186 insertions(+), 25 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index e4a4b2c..7eee8cd 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -61,29 +61,33 @@ namespace PPDBM constexpr DbmBound( int32_t value, bool strict): _value{value}, _strict{strict}{} - constexpr bool isInfinite() const { return _value == INF; } ///< is this constraint unrestraining? + /// 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)); } ///< is the new constraint useless? + 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())); - } ///< do the two constraints force the zone to be empty? + } private: int32_t _value : 31; - bool _strict : 1; ///< constant coeff of the constraint, is the constraint strict? + bool _strict : 1; /// constant coeff of the constraint, is the constraint strict? }; constexpr auto INF_BOUND = DbmBound{INF,true}; constexpr auto DIAG_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() }; - } ///< x_i - x_j <= a and x_j - x_k < b gives x_i - x_k < a + b + } /** @@ -111,12 +115,15 @@ namespace PPDBM constexpr auto getDim() const { return _matrix.getRows(); } constexpr DbmBound& operator()(int i, int j) { return _matrix(i, j); } constexpr const DbmBound& operator()(int i, int j) const { return _matrix(i, j); } - void reset() + constexpr void reset() { _matrix.fill(INF_BOUND); for (uint32_t i = 0; i < _matrix.getRows(); ++i) { _matrix(i,i)= DIAG_BOUND; } } - void close(uint32_t i, uint32_t j); + constexpr void close(uint32_t i, uint32_t j); + + /// Relation géométrique pure entre deux zones (ignorant le prix). + relation_t relation(const DBMatrix& other) const; private: Matrix _matrix; @@ -124,32 +131,49 @@ namespace PPDBM + using Vertex = std::vector; // taille _dim, valuation d'horloges + using CostVector = std::vector; // taille _param + 1, coeffs sur les paramètres + /** * Parametric priced timed difference bound matrix */ struct Ppdbm { - Ppdbm(uint32_t dim, uint32_t param): _dim{dim}, _param{param}, _dbm{dim}, + constexpr Ppdbm(uint32_t dim, uint32_t param): _dim{dim}, _param{param}, _dbm{dim}, _offsetCost(param + 1), _rates(dim, param + 1, 0) {} - 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; + constexpr bool isEmpty() const { return _emptyZone; } + constexpr void reset(); + constexpr bool constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint); + constexpr void removeCostAtOffset(); + constexpr void addCostAtOffset(); + constexpr bool isDiagonalNegative(uint32_t i) const; + constexpr bool constrain_DBM_N(const std::vector& constraints); + constexpr std::vector costAtOtherOffset(const DBMatrix& otherDbm) const; + relation_t relation(const Ppdbm& other) const; 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 + + /// Sommets de la zone representee par `dbm`. TODO: a ecrire ensemble. + std::vector sommetsDBM(const DBMatrix& dbm) const; + + /// Evalue le vecteur de cout (coeffs sur les parametres) de *this* au sommet `v`. + CostVector costAtVertex(const Vertex& v) const; + + /// a domine b (a moins cher ou egal) coefficient par coefficient. + static bool dominatesVector(const CostVector& a, const CostVector& b); + static bool equalVector(const CostVector& a, const CostVector& b); + + /// Compare deux listes de vecteurs de cout alignees sommet par sommet. + relation_t compareCostLists(const std::vector& c1, const std::vector& c2) const; }; } // namespace PPDBM -#endif +#endif \ No newline at end of file diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index 5384f83..39d8abc 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -19,7 +19,7 @@ namespace PPDBM * @param i * @param j */ - void DBMatrix::close(uint32_t i, uint32_t j) + constexpr void DBMatrix::close(uint32_t i, uint32_t j) { const auto dim = getDim(); const auto dij = _matrix(i, j); @@ -38,7 +38,7 @@ namespace PPDBM /** * resets a DBM to its unrestrained state. */ - void Ppdbm::reset() + constexpr void Ppdbm::reset() { _dbm.reset(); std::fill(_offsetCost.begin(), _offsetCost.end(), 0); @@ -50,7 +50,7 @@ namespace PPDBM /** * */ - void Ppdbm::removeCostAtOffset() + constexpr void Ppdbm::removeCostAtOffset() { for (uint32_t h = 1; h < _dim; h++) { for (uint32_t k = 0; k <= _param; k++) { @@ -62,7 +62,7 @@ namespace PPDBM /** * */ - void Ppdbm::addCostAtOffset() + constexpr void Ppdbm::addCostAtOffset() { for (uint32_t h = 1; h < _dim; h++) { for (uint32_t k = 0; k <= _param; k++) { @@ -76,7 +76,7 @@ namespace PPDBM * @param i * @return */ - bool Ppdbm::isDiagonalNegative(uint32_t i) const + constexpr bool Ppdbm::isDiagonalNegative(uint32_t i) const { const auto& d = _dbm(i, i); return d.getValue() < 0 || (d.getValue() == 0 && d.getStrict()); @@ -89,7 +89,7 @@ namespace PPDBM * @param constraint * @return */ - bool Ppdbm::constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint) + constexpr bool Ppdbm::constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint) { assert(i < _dim && j < _dim); @@ -122,7 +122,7 @@ namespace PPDBM * @param constraints * @return */ - bool Ppdbm::constrain_DBM_N(const std::vector& constraints) + constexpr bool Ppdbm::constrain_DBM_N(const std::vector& constraints) { assert(!constraints.empty()); @@ -158,7 +158,7 @@ namespace PPDBM * @param otherDbm * @return */ - std::vector Ppdbm::costAtOtherOffset(const DBMatrix& otherDbm) const + constexpr std::vector Ppdbm::costAtOtherOffset(const DBMatrix& otherDbm) const { std::vector cost = _offsetCost; @@ -174,6 +174,143 @@ namespace PPDBM return cost; } + + + + + + + + + + + + + + ///////////////////////////////////////////////////attempt at relation///////////////////////////////////////////// + + relation_t DBMatrix::relation(const DBMatrix& other) const + { + assert(getDim() == other.getDim()); + const auto dim = getDim(); + + bool oneLeqTwo = true; // this <= other partout -> this subset of other + bool twoLeqOne = true; // other <= this partout -> 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))) { oneLeqTwo = false; } + if (!other(i, j).dominates((*this)(i, j))) { twoLeqOne = false; } + } + } + + if (oneLeqTwo && twoLeqOne) { return base_EQUAL; } + if (oneLeqTwo) { return base_SUBSET; } + if (twoLeqOne) { return base_SUPERSET; } + return base_DIFFERENT; + } + + CostVector Ppdbm::costAtVertex(const Vertex& v) const + { + assert(v.size() == _dim); + + CostVector cost = _offsetCost; // copie: coeffs a l'offset + for (uint32_t h = 1; h < _dim; ++h) { + const int32_t offsetH = -_dbm(0, h).getValue(); // valeur min de l'horloge h + const int32_t delta = v[h] - offsetH; + for (uint32_t k = 0; k <= _param; ++k) { + cost[k] += _rates(h, k) * delta; + } + } + return cost; + } + + bool Ppdbm::dominatesVector(const CostVector& a, const CostVector& b) + { + assert(a.size() == b.size()); + for (size_t k = 0; k < a.size(); ++k) { + if (a[k] > b[k]) { return false; } + } + return true; + } + + bool Ppdbm::equalVector(const CostVector& a, const CostVector& b) + { + return a == b; + } + + relation_t Ppdbm::compareCostLists(const std::vector& c1, + const std::vector& c2) const + { + assert(c1.size() == c2.size()); + + bool oneLeqTwo = true; // cout de *this* <= cout de other, sur tous les sommets + bool twoLeqOne = true; + + for (size_t v = 0; v < c1.size(); ++v) { + if (!dominatesVector(c1[v], c2[v])) { oneLeqTwo = false; } + if (!dominatesVector(c2[v], c1[v])) { twoLeqOne = false; } + if (!oneLeqTwo && !twoLeqOne) { break; } // court-circuit: plus rien a decider + } + + if (oneLeqTwo && twoLeqOne) { return base_EQUAL; } + if (oneLeqTwo) { return base_SUPERSET; } // *this* moins cher partout -> domine + if (twoLeqOne) { return base_SUBSET; } + return base_DIFFERENT; // pas de domination (cf. ta preuve, on ne coupe pas la branche) + } + + 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 = sommetsDBM(_dbm); // meme zone des deux cotes + 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* est la zone la plus grande, other la plus petite. + auto vertices = sommetsDBM(other._dbm); + 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 est la zone la plus grande, *this* la plus petite. + auto vertices = sommetsDBM(_dbm); + 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; + } + } } // namespace PPBDM From 2f140131b8c271e2ec951bcd2a8ba970a8f6a239 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Fri, 21 Aug 2026 17:34:19 +0200 Subject: [PATCH 08/18] relation function finished! 1st version, needs tests. --- include/dbm/ParamPricedDBM.h | 4 +- src/ParamPricedDBM.cpp | 92 +++++++++++++++++++++++++----------- 2 files changed, 66 insertions(+), 30 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index 7eee8cd..269c799 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -131,8 +131,8 @@ namespace PPDBM - using Vertex = std::vector; // taille _dim, valuation d'horloges - using CostVector = std::vector; // taille _param + 1, coeffs sur les paramètres + using Vertex = std::vector; // size _dim, clock valuations + using CostVector = std::vector; // size _param + 1, coeffs on parameters /** * Parametric priced timed difference bound matrix diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index 39d8abc..1d71eba 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -189,76 +189,112 @@ namespace PPDBM ///////////////////////////////////////////////////attempt at 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 oneLeqTwo = true; // this <= other partout -> this subset of other - bool twoLeqOne = true; // other <= this partout -> other subset of this + 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))) { oneLeqTwo = false; } - if (!other(i, j).dominates((*this)(i, j))) { twoLeqOne = false; } + if (!(*this)(i, j).dominates(other(i, j))) { thisLeqOther = false; } + if (!other(i, j).dominates((*this)(i, j))) { otherLeqThis = false; } } } - if (oneLeqTwo && twoLeqOne) { return base_EQUAL; } - if (oneLeqTwo) { return base_SUBSET; } - if (twoLeqOne) { return base_SUPERSET; } + 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; // copie: coeffs a l'offset - for (uint32_t h = 1; h < _dim; ++h) { - const int32_t offsetH = -_dbm(0, h).getValue(); // valeur min de l'horloge h - const int32_t delta = v[h] - offsetH; + 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(h, k) * delta; + 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 (size_t k = 0; k < a.size(); ++k) { + 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; } - relation_t Ppdbm::compareCostLists(const std::vector& c1, - const std::vector& c2) const + /** + * 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) const { assert(c1.size() == c2.size()); - bool oneLeqTwo = true; // cout de *this* <= cout de other, sur tous les sommets - bool twoLeqOne = true; + 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 (size_t v = 0; v < c1.size(); ++v) { - if (!dominatesVector(c1[v], c2[v])) { oneLeqTwo = false; } - if (!dominatesVector(c2[v], c1[v])) { twoLeqOne = false; } - if (!oneLeqTwo && !twoLeqOne) { break; } // court-circuit: plus rien a decider + 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 (oneLeqTwo && twoLeqOne) { return base_EQUAL; } - if (oneLeqTwo) { return base_SUPERSET; } // *this* moins cher partout -> domine - if (twoLeqOne) { return base_SUBSET; } - return base_DIFFERENT; // pas de domination (cf. ta preuve, on ne coupe pas la branche) + if (thisLeqOther && otherLeqThis) { return base_EQUAL; } + if (thisLeqOther) { return base_SUPERSET; } // *this* cheaper everywhere -> dominates + if (otherLeqThis) { return base_SUBSET; } + return base_DIFFERENT; // no domination --> no branch cutting } + /** + * most important function, 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); @@ -268,7 +304,7 @@ namespace PPDBM return base_DIFFERENT; case base_EQUAL: { - auto vertices = sommetsDBM(_dbm); // meme zone des deux cotes + auto vertices = sommetsDBM(_dbm); // same zone on both sides std::vector c1, c2; c1.reserve(vertices.size()); c2.reserve(vertices.size()); @@ -280,7 +316,7 @@ namespace PPDBM } case base_SUPERSET: { - // *this* est la zone la plus grande, other la plus petite. + // *this* is the bigger zone, other is the smaller auto vertices = sommetsDBM(other._dbm); std::vector c1, c2; c1.reserve(vertices.size()); @@ -294,7 +330,7 @@ namespace PPDBM } case base_SUBSET: { - // other est la zone la plus grande, *this* la plus petite. + // other is the bigger zone, *this* is the smaller: auto vertices = sommetsDBM(_dbm); std::vector c1, c2; c1.reserve(vertices.size()); From 0e534a0ddccc4dddeabdae5fc7d86ea9f2bbac8a Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Mon, 24 Aug 2026 14:36:53 +0200 Subject: [PATCH 09/18] verticesDBM finished as well --- include/dbm/ParamPricedDBM.h | 45 ++++---- src/ParamPricedDBM.cpp | 206 +++++++++++++++++++++++++++++++---- 2 files changed, 208 insertions(+), 43 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index 269c799..04450c4 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -2,16 +2,11 @@ #define INCLUDE_DBM_PARAMPRICEDDBM_H #include "dbm.h" -#include "mingraph.h" -#include "constraints.h" -#include #include #include -#include #include #include -#include #include @@ -102,6 +97,8 @@ namespace PPDBM + using Vertex = std::vector; + /** * */ @@ -122,8 +119,10 @@ namespace PPDBM } constexpr void close(uint32_t i, uint32_t j); - /// Relation géométrique pure entre deux zones (ignorant le prix). - relation_t relation(const DBMatrix& other) const; + /// pure zone comparison (no cost involved). + constexpr relation_t relation(const DBMatrix& other) const; + constexpr std::vector verticesDBM() const; + private: Matrix _matrix; @@ -132,7 +131,7 @@ namespace PPDBM using Vertex = std::vector; // size _dim, clock valuations - using CostVector = std::vector; // size _param + 1, coeffs on parameters + using CostVector = std::vector; // size _param + 1, coeffs on parameters /** * Parametric priced timed difference bound matrix @@ -143,6 +142,10 @@ namespace PPDBM _offsetCost(param + 1), _rates(dim, param + 1, 0) {} + constexpr uint32_t getDim() const { return _dim; } + constexpr DbmBound& operator()(int i, int j) { return _dbm(i, j); } + constexpr const DbmBound& operator()(int i, int j) const { return _dbm(i, j); } + constexpr bool isEmpty() const { return _emptyZone; } constexpr void reset(); constexpr bool constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint); @@ -151,7 +154,18 @@ namespace PPDBM constexpr bool isDiagonalNegative(uint32_t i) const; constexpr bool constrain_DBM_N(const std::vector& constraints); constexpr std::vector costAtOtherOffset(const DBMatrix& otherDbm) const; - relation_t relation(const Ppdbm& other) const; + + /// Is there a domination relation between the two zones? + constexpr relation_t relation(const Ppdbm& other) const; + /// Vertices of the zone. + constexpr std::vector verticesDBM(const DBMatrix& dbm) const{ return _dbm.verticesDBM(); }; + /// Evaluates the cost vector(coeffs on parameters) of *this* at vertex `v`. + constexpr CostVector costAtVertex(const Vertex& v) const; + /// a dominates b coefficient by coefficient? + constexpr static bool dominatesVector(const CostVector& a, const CostVector& b); + constexpr static bool equalVector(const CostVector& a, const CostVector& b); + /// Compare two lists of cost vectors aligned vertex by vertex. + constexpr static relation_t compareCostLists(const std::vector& c1, const std::vector& c2); private: uint32_t _dim, _param = 0; // number of clocks and parameters @@ -160,19 +174,6 @@ namespace PPDBM 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 - - /// Sommets de la zone representee par `dbm`. TODO: a ecrire ensemble. - std::vector sommetsDBM(const DBMatrix& dbm) const; - - /// Evalue le vecteur de cout (coeffs sur les parametres) de *this* au sommet `v`. - CostVector costAtVertex(const Vertex& v) const; - - /// a domine b (a moins cher ou egal) coefficient par coefficient. - static bool dominatesVector(const CostVector& a, const CostVector& b); - static bool equalVector(const CostVector& a, const CostVector& b); - - /// Compare deux listes de vecteurs de cout alignees sommet par sommet. - relation_t compareCostLists(const std::vector& c1, const std::vector& c2) const; }; } // namespace PPDBM diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index 1d71eba..d3717a6 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -6,8 +6,6 @@ #include #include -#include -#include #include #include #include @@ -178,16 +176,7 @@ namespace PPDBM - - - - - - - - - - ///////////////////////////////////////////////////attempt at relation///////////////////////////////////////////// + /////////////////////////////////////////////////// relation ///////////////////////////////////////////////////// /** @@ -195,7 +184,7 @@ namespace PPDBM * @param other * @return */ - relation_t DBMatrix::relation(const DBMatrix& other) const + constexpr relation_t DBMatrix::relation(const DBMatrix& other) const { assert(getDim() == other.getDim()); const auto dim = getDim(); @@ -221,7 +210,7 @@ namespace PPDBM * @param v * @return */ - CostVector Ppdbm::costAtVertex(const Vertex& v) const + constexpr CostVector Ppdbm::costAtVertex(const Vertex& v) const { assert(v.size() == _dim); @@ -242,7 +231,7 @@ namespace PPDBM * @param b * @return */ - bool Ppdbm::dominatesVector(const CostVector& a, const CostVector& b) + constexpr bool Ppdbm::dominatesVector(const CostVector& a, const CostVector& b) { assert(a.size() == b.size()); for (uint32_t k = 0; k < a.size(); ++k) { @@ -257,7 +246,7 @@ namespace PPDBM * @param b * @return */ - bool Ppdbm::equalVector(const CostVector& a, const CostVector& b) + constexpr bool Ppdbm::equalVector(const CostVector& a, const CostVector& b) { return a == b; } @@ -269,7 +258,7 @@ namespace PPDBM * @param c2 * @return */ - relation_t Ppdbm::compareCostLists(const std::vector& c1, const std::vector& c2) const + constexpr relation_t Ppdbm::compareCostLists(const std::vector& c1, const std::vector& c2) { assert(c1.size() == c2.size()); @@ -295,7 +284,7 @@ namespace PPDBM * @param other * @return */ - relation_t Ppdbm::relation(const Ppdbm& other) const + constexpr relation_t Ppdbm::relation(const Ppdbm& other) const { assert(_dim == other._dim && _param == other._param); @@ -304,7 +293,7 @@ namespace PPDBM return base_DIFFERENT; case base_EQUAL: { - auto vertices = sommetsDBM(_dbm); // same zone on both sides + auto vertices = _dbm.verticesDBM(); // same zone on both sides std::vector c1, c2; c1.reserve(vertices.size()); c2.reserve(vertices.size()); @@ -317,7 +306,7 @@ namespace PPDBM case base_SUPERSET: { // *this* is the bigger zone, other is the smaller - auto vertices = sommetsDBM(other._dbm); + auto vertices = other._dbm.verticesDBM(); std::vector c1, c2; c1.reserve(vertices.size()); c2.reserve(vertices.size()); @@ -331,7 +320,7 @@ namespace PPDBM case base_SUBSET: { // other is the bigger zone, *this* is the smaller: - auto vertices = sommetsDBM(_dbm); + auto vertices = _dbm.verticesDBM(); std::vector c1, c2; c1.reserve(vertices.size()); c2.reserve(vertices.size()); @@ -347,6 +336,181 @@ namespace PPDBM return base_DIFFERENT; } } + + + + + + + + + + + + + +//////////////////////////////////////////////// sommetsDBM() with claude's help ///////////////////////////////////// + + + /** + * finds the vertices of a dbm! very important function! + * @return + */ + constexpr 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 : énumère tous les arbres couvrants sur '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; + } } // namespace PPBDM From 22b2b6fbc28baf05d809f08b916d5264299193fa Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Mon, 24 Aug 2026 20:21:32 +0200 Subject: [PATCH 10/18] test added --- include/dbm/ParamPricedDBM.h | 1 + src/ParamPricedDBM.cpp | 3 - test/CMakeLists.txt | 4 +- test/test_ParamPricedDBM.cpp | 330 +++++++++++++++++++++++++++++++++++ 4 files changed, 334 insertions(+), 4 deletions(-) create mode 100644 test/test_ParamPricedDBM.cpp diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index 04450c4..2e27af5 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -145,6 +145,7 @@ namespace PPDBM constexpr uint32_t getDim() const { return _dim; } 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 setRates(const Matrix& rates) { _rates = rates; } constexpr bool isEmpty() const { return _emptyZone; } constexpr void reset(); diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index d3717a6..9adc296 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -346,9 +346,6 @@ namespace PPDBM - - - //////////////////////////////////////////////// sommetsDBM() with claude's help ///////////////////////////////////// diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4f7db19..6cdd62e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -9,7 +9,7 @@ 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}) @@ -37,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..ce68fd0 --- /dev/null +++ b/test/test_ParamPricedDBM.cpp @@ -0,0 +1,330 @@ +#include "dbm/ParamPricedDBM.h" + +#include + +#include +#include + +using namespace PPDBM; + +TEST_CASE("DbmBound dominates") +{ + // valeur plus petite => contrainte plus forte => domine + const auto a = DbmBound{3, false}; // <= 3 + const auto b = DbmBound{5, false}; // <= 5 + CHECK(a.dominates(b)); + CHECK(!b.dominates(a)); + + // meme valeur : strict domine non-strict (< 3 est plus fort que <= 3) + const auto strict3 = DbmBound{3, true}; + const auto loose3 = DbmBound{3, false}; + CHECK(strict3.dominates(loose3)); + CHECK(!loose3.dominates(strict3)); + + // une borne se domine toujours elle-meme + CHECK(a.dominates(a)); + + // INF ne domine jamais une borne finie, l'inverse est toujours vrai + CHECK(!INF_BOUND.dominates(a)); + CHECK(a.dominates(INF_BOUND)); +} + +TEST_CASE("DbmBound isIncompatible") +{ + // x_i - x_j <= 3 et x_j - x_i <= -4 => x_i - x_j >= 4, incompatible avec <=3 + const auto a = DbmBound{3, false}; + const auto b = DbmBound{-4, false}; + CHECK(DbmBound::isIncompatible(a, b)); + + // x_i - x_j <= 3 et x_j - x_i <= -3 => point unique {x_i - x_j = 3}, compatible + const auto c = DbmBound{-3, false}; + CHECK(!DbmBound::isIncompatible(a, c)); + + // si l'une des deux bornes est stricte, la somme nulle devient incompatible + const auto aStrict = DbmBound{3, true}; + CHECK(DbmBound::isIncompatible(aStrict, c)); + + // deux bornes larges, somme > 0 => compatible + const auto d = DbmBound{1, false}; + const auto e = DbmBound{1, false}; + CHECK(!DbmBound::isIncompatible(d, e)); +} + +TEST_CASE("DbmBound operator+") +{ + const auto a = DbmBound{3, false}; + const auto b = DbmBound{2, false}; + const auto sum = a + b; + CHECK(sum.getValue() == 5); + CHECK(!sum.getStrict()); + + // si l'une des deux est stricte, le resultat est strict + const auto bStrict = DbmBound{2, true}; + const auto sum2 = a + bStrict; + CHECK(sum2.getValue() == 5); + CHECK(sum2.getStrict()); + + // INF + n'importe quoi = INF + const auto sum3 = a + INF_BOUND; + CHECK(sum3.isInfinite()); +} + +TEST_CASE("DBMatrix construction and reset") +{ + auto m = DBMatrix{3}; + REQUIRE(m.getDim() == 3); + + SUBCASE("etat initial : diagonale a 0, reste a INF") + { + for (uint32_t i = 0; i < 3; ++i) { + CHECK(m(i, i).getValue() == 0); + CHECK(!m(i, i).getStrict()); + for (uint32_t j = 0; j < 3; ++j) { + if (i != j) CHECK(m(i, j).isInfinite()); + } + } + } + + SUBCASE("reset revient a l'etat initial") + { + m(0, 1) = DbmBound{7, false}; + m.reset(); + CHECK(m(0, 1).isInfinite()); + CHECK(m(0, 0).getValue() == 0); + CHECK(m(1, 1).getValue() == 0); + } +} + +TEST_CASE("DBMatrix close") +{ + SUBCASE("propagation simple : x1-x0<=3, x2-x1<=2 => x2-x0<=5") + { + auto m = DBMatrix{3}; + m(1, 0) = DbmBound{3, false}; + m.close(1, 0); + m(2, 1) = DbmBound{2, false}; + m.close(2, 1); + + CHECK(m(2, 0).getValue() == 5); + CHECK(!m(2, 0).getStrict()); + } + + SUBCASE("propagation d'une borne stricte") + { + // x1-x0<3 (stricte), x2-x1<=2 (large) => x2-x0<5 (stricte) + auto m = DBMatrix{3}; + m(1, 0) = DbmBound{3, true}; + m.close(1, 0); + m(2, 1) = DbmBound{2, false}; + m.close(2, 1); + + CHECK(m(2, 0).getValue() == 5); + CHECK(m(2, 0).getStrict()); + } +} + +TEST_CASE("DBMatrix relation") +{ + // m1: x1<=5, m2: x1<=3 => m2 subset de m1 + auto m1 = DBMatrix{2}; + m1(1, 0) = DbmBound{5, false}; + m1.close(1, 0); + + auto m2 = DBMatrix{2}; + m2(1, 0) = DbmBound{3, false}; + m2.close(1, 0); + + CHECK(m1.relation(m2) == base_SUPERSET); + CHECK(m2.relation(m1) == base_SUBSET); + + auto m3 = DBMatrix{2}; + m3(1, 0) = DbmBound{5, false}; + m3.close(1, 0); + CHECK(m1.relation(m3) == base_EQUAL); +} + +TEST_CASE("DBMatrix verticesDBM") +{ + // 2 horloges : x0 = reference (toujours 0), x1 dans [2,5] + auto m = DBMatrix{2}; + m(1, 0) = DbmBound{5, false}; // x1 - x0 <= 5 + m(0, 1) = DbmBound{-2, false}; // x0 - x1 <= -2, donc x1 >= 2 + + const auto vertices = m.verticesDBM(); + const auto v25 = Vertex{0, 5}; + const auto v22 = Vertex{0, 2}; + + REQUIRE(vertices.size() == 2); + CHECK(std::find(vertices.begin(), vertices.end(), v25) != vertices.end()); + CHECK(std::find(vertices.begin(), vertices.end(), v22) != vertices.end()); +} + +TEST_CASE("Ppdbm reset") +{ + auto p = Ppdbm{2, 1}; + p.constrain_DBM(1, 0, DbmBound{4, false}); + p.reset(); + CHECK(!p.isEmpty()); + CHECK(p(1, 0).isInfinite()); +} + +TEST_CASE("Ppdbm constrain_DBM") +{ + SUBCASE("une contrainte plus faible que l'existante est ignoree") + { + auto p = Ppdbm{2, 0}; + CHECK(p.constrain_DBM(1, 0, DbmBound{5, false})); + CHECK(p.constrain_DBM(1, 0, DbmBound{9, false})); + CHECK(p(1, 0).getValue() == 5); + } + + SUBCASE("une contrainte plus forte resserre la zone") + { + auto p = Ppdbm{2, 0}; + CHECK(p.constrain_DBM(1, 0, DbmBound{5, false})); + CHECK(p.constrain_DBM(1, 0, DbmBound{3, false})); + CHECK(p(1, 0).getValue() == 3); + CHECK(!p.isEmpty()); + } + + SUBCASE("une contrainte incompatible vide la zone") + { + // x1-x0<=3 puis x0-x1<=-4 (x1>=4) => incompatible avec x1<=3 + auto p = Ppdbm{2, 0}; + CHECK(p.constrain_DBM(1, 0, DbmBound{3, false})); + CHECK(!p.constrain_DBM(0, 1, DbmBound{-4, false})); + CHECK(p.isEmpty()); + } +} + +TEST_CASE("Ppdbm constrain_DBM_N") +{ + SUBCASE("plusieurs contraintes compatibles appliquees d'un coup") + { + // x1 dans [2,5], x2 dans [1,4] + auto p = Ppdbm{3, 0}; + const auto constraints = std::vector{ + {1, 0, DbmBound{5, false}}, + {0, 1, DbmBound{-2, false}}, + {2, 0, DbmBound{4, false}}, + {0, 2, DbmBound{-1, false}}, + }; + CHECK(p.constrain_DBM_N(constraints)); + CHECK(!p.isEmpty()); + CHECK(p(1, 0).getValue() == 5); + CHECK(p(0, 1).getValue() == -2); + CHECK(p(2, 0).getValue() == 4); + CHECK(p(0, 2).getValue() == -1); + } + + SUBCASE("une contrainte incompatible dans la liste vide la zone") + { + auto p = Ppdbm{2, 0}; + const auto constraints = std::vector{ + {1, 0, DbmBound{3, false}}, + {0, 1, DbmBound{-4, false}}, + }; + CHECK(!p.constrain_DBM_N(constraints)); + CHECK(p.isEmpty()); + } +} + +TEST_CASE("Ppdbm isDiagonalNegative") +{ + auto p = Ppdbm{2, 0}; + CHECK(!p.isDiagonalNegative(0)); + CHECK(!p.isDiagonalNegative(1)); + + p.constrain_DBM(1, 0, DbmBound{3, false}); + p.constrain_DBM(0, 1, DbmBound{-4, false}); + CHECK(p.isDiagonalNegative(1)); +} + +TEST_CASE("Ppdbm costAtVertex") +{ + // 1 parametre. rate(1,0)=2 (constante), rate(1,1)=1 (coeff parametre) + auto p = Ppdbm{2, 1}; + auto rates = Matrix{2, 2, 0}; + rates(1, 0) = 2; + rates(1, 1) = 1; + p.setRates(rates); + p.constrain_DBM(1, 0, DbmBound{5, false}); + p.constrain_DBM(0, 1, DbmBound{-2, false}); // x1 dans [2,5] + + const auto v = Vertex{0, 5}; // delta = 5 - 2 = 3 + const auto cost = p.costAtVertex(v); + REQUIRE(cost.size() == 2); + CHECK(cost[0] == 6); // 0 + 2*3 + CHECK(cost[1] == 3); // 0 + 1*3 +} + +TEST_CASE("Ppdbm dominatesVector and equalVector") +{ + const auto a = CostVector{1, 2, 3}; + const auto b = CostVector{2, 2, 4}; + const auto c = CostVector{1, 2, 3}; + const auto d = CostVector{1, 2, 4}; + + CHECK(Ppdbm::dominatesVector(a, b)); + CHECK(!Ppdbm::dominatesVector(b, a)); + CHECK(Ppdbm::dominatesVector(a, c)); // egalite => domine (<=) + + CHECK(Ppdbm::equalVector(a, c)); + CHECK(!Ppdbm::equalVector(a, d)); +} + +TEST_CASE("Ppdbm compareCostLists") +{ + const auto c1 = std::vector{{1, 1}, {2, 2}}; + const auto c2 = std::vector{{1, 1}, {2, 2}}; + CHECK(Ppdbm::compareCostLists(c1, c2) == base_EQUAL); + + const auto cheaper = std::vector{{1, 1}, {1, 1}}; + const auto costly = std::vector{{2, 2}, {3, 3}}; + // cheaper <= costly partout => cheaper "domine" => base_SUPERSET + CHECK(Ppdbm::compareCostLists(cheaper, costly) == base_SUPERSET); + CHECK(Ppdbm::compareCostLists(costly, cheaper) == base_SUBSET); + + const auto mixed1 = std::vector{{1, 5}, {5, 1}}; + const auto mixed2 = std::vector{{5, 1}, {1, 5}}; + CHECK(Ppdbm::compareCostLists(mixed1, mixed2) == base_DIFFERENT); +} + +TEST_CASE("Ppdbm relation") +{ + SUBCASE("meme zone, couts differents => domination par le cout") + { + auto p1 = Ppdbm{2, 1}; + auto p2 = Ppdbm{2, 1}; + p1.constrain_DBM(1, 0, DbmBound{5, false}); + p1.constrain_DBM(0, 1, DbmBound{-2, false}); + p2.constrain_DBM(1, 0, DbmBound{5, false}); + p2.constrain_DBM(0, 1, DbmBound{-2, false}); + + auto rates1 = Matrix{2, 2, 0}; + rates1(1, 1) = 1; // cout = x1 + p1.setRates(rates1); + + auto rates2 = Matrix{2, 2, 0}; + rates2(1, 1) = 2; // cout = 2*x1, toujours plus cher (ou egal) + p2.setRates(rates2); + + // meme zone, p1 moins cher partout => p1 domine (SUPERSET), p2 est domine (SUBSET) + CHECK(p1.relation(p2) == base_SUPERSET); + CHECK(p2.relation(p1) == base_SUBSET); + } + + SUBCASE("zones differentes, le resultat reste stable") + { + auto p1 = Ppdbm{2, 0}; + auto p2 = Ppdbm{2, 0}; + p1.constrain_DBM(1, 0, DbmBound{5, false}); + p2.constrain_DBM(1, 0, DbmBound{3, false}); + p2.constrain_DBM(0, 1, DbmBound{-1, false}); + + const auto r1 = p1.relation(p2); + const auto r2 = p2.relation(p1); + CHECK((r1 == base_SUPERSET || r1 == base_DIFFERENT)); + CHECK((r2 == base_SUBSET || r2 == base_DIFFERENT)); + } +} \ No newline at end of file From 93c1f33929d9f75be928bc4b43ada3a339ad73d1 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Tue, 25 Aug 2026 09:12:10 +0200 Subject: [PATCH 11/18] every function and their tests up to now working --- include/dbm/ParamPricedDBM.h | 38 +-- src/ParamPricedDBM.cpp | 30 +-- test/test_ParamPricedDBM.cpp | 500 +++++++++++++++++------------------ 3 files changed, 275 insertions(+), 293 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index 2e27af5..7ee9a05 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -106,7 +106,7 @@ namespace PPDBM { constexpr DBMatrix(uint32_t dim) : _matrix{dim, dim, INF_BOUND} { - for (uint32_t i = 0; i < dim; ++i) { _matrix(i, i) = DIAG_BOUND; } + for (uint32_t i = 0; i < dim; ++i) { _matrix(i, i) = DIAG_BOUND; _matrix(0,i) = DIAG_BOUND; } } constexpr auto getDim() const { return _matrix.getRows(); } @@ -115,13 +115,15 @@ namespace PPDBM constexpr void reset() { _matrix.fill(INF_BOUND); - for (uint32_t i = 0; i < _matrix.getRows(); ++i) { _matrix(i,i)= DIAG_BOUND; } + for (uint32_t i = 0; i < _matrix.getRows(); ++i) { + _matrix(i,i)= DIAG_BOUND; + _matrix(0,i) = DIAG_BOUND; + } } - constexpr void close(uint32_t i, uint32_t j); + void close(uint32_t i, uint32_t j); - /// pure zone comparison (no cost involved). - constexpr relation_t relation(const DBMatrix& other) const; - constexpr std::vector verticesDBM() const; + relation_t relation(const DBMatrix& other) const; ///< pure zone comparison (no cost involved). + std::vector verticesDBM() const; private: @@ -148,25 +150,25 @@ namespace PPDBM void setRates(const Matrix& rates) { _rates = rates; } constexpr bool isEmpty() const { return _emptyZone; } - constexpr void reset(); - constexpr bool constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint); - constexpr void removeCostAtOffset(); - constexpr void addCostAtOffset(); - constexpr bool isDiagonalNegative(uint32_t i) const; - constexpr bool constrain_DBM_N(const std::vector& constraints); - constexpr std::vector costAtOtherOffset(const DBMatrix& otherDbm) const; + 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; /// Is there a domination relation between the two zones? - constexpr relation_t relation(const Ppdbm& other) const; + relation_t relation(const Ppdbm& other) const; /// Vertices of the zone. constexpr std::vector verticesDBM(const DBMatrix& dbm) const{ return _dbm.verticesDBM(); }; /// Evaluates the cost vector(coeffs on parameters) of *this* at vertex `v`. - constexpr CostVector costAtVertex(const Vertex& v) const; + CostVector costAtVertex(const Vertex& v) const; /// a dominates b coefficient by coefficient? - constexpr static bool dominatesVector(const CostVector& a, const CostVector& b); - constexpr static bool equalVector(const CostVector& a, const CostVector& b); + 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. - constexpr static relation_t compareCostLists(const std::vector& c1, const std::vector& c2); + static relation_t compareCostLists(const std::vector& c1, const std::vector& c2); private: uint32_t _dim, _param = 0; // number of clocks and parameters diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index 9adc296..0590dfe 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -17,7 +17,7 @@ namespace PPDBM * @param i * @param j */ - constexpr void DBMatrix::close(uint32_t i, uint32_t j) + void DBMatrix::close(uint32_t i, uint32_t j) { const auto dim = getDim(); const auto dij = _matrix(i, j); @@ -36,7 +36,7 @@ namespace PPDBM /** * resets a DBM to its unrestrained state. */ - constexpr void Ppdbm::reset() + void Ppdbm::reset() { _dbm.reset(); std::fill(_offsetCost.begin(), _offsetCost.end(), 0); @@ -48,7 +48,7 @@ namespace PPDBM /** * */ - constexpr void Ppdbm::removeCostAtOffset() + void Ppdbm::removeCostAtOffset() { for (uint32_t h = 1; h < _dim; h++) { for (uint32_t k = 0; k <= _param; k++) { @@ -60,7 +60,7 @@ namespace PPDBM /** * */ - constexpr void Ppdbm::addCostAtOffset() + void Ppdbm::addCostAtOffset() { for (uint32_t h = 1; h < _dim; h++) { for (uint32_t k = 0; k <= _param; k++) { @@ -74,7 +74,7 @@ namespace PPDBM * @param i * @return */ - constexpr bool Ppdbm::isDiagonalNegative(uint32_t i) const + bool Ppdbm::isDiagonalNegative(uint32_t i) const { const auto& d = _dbm(i, i); return d.getValue() < 0 || (d.getValue() == 0 && d.getStrict()); @@ -87,7 +87,7 @@ namespace PPDBM * @param constraint * @return */ - constexpr bool Ppdbm::constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint) + bool Ppdbm::constrain_DBM(uint32_t i, uint32_t j, DbmBound constraint) { assert(i < _dim && j < _dim); @@ -120,7 +120,7 @@ namespace PPDBM * @param constraints * @return */ - constexpr bool Ppdbm::constrain_DBM_N(const std::vector& constraints) + bool Ppdbm::constrain_DBM_N(const std::vector& constraints) { assert(!constraints.empty()); @@ -156,7 +156,7 @@ namespace PPDBM * @param otherDbm * @return */ - constexpr std::vector Ppdbm::costAtOtherOffset(const DBMatrix& otherDbm) const + std::vector Ppdbm::costAtOtherOffset(const DBMatrix& otherDbm) const { std::vector cost = _offsetCost; @@ -184,7 +184,7 @@ namespace PPDBM * @param other * @return */ - constexpr relation_t DBMatrix::relation(const DBMatrix& other) const + relation_t DBMatrix::relation(const DBMatrix& other) const { assert(getDim() == other.getDim()); const auto dim = getDim(); @@ -210,7 +210,7 @@ namespace PPDBM * @param v * @return */ - constexpr CostVector Ppdbm::costAtVertex(const Vertex& v) const + CostVector Ppdbm::costAtVertex(const Vertex& v) const { assert(v.size() == _dim); @@ -231,7 +231,7 @@ namespace PPDBM * @param b * @return */ - constexpr bool Ppdbm::dominatesVector(const CostVector& a, const CostVector& b) + bool Ppdbm::dominatesVector(const CostVector& a, const CostVector& b) { assert(a.size() == b.size()); for (uint32_t k = 0; k < a.size(); ++k) { @@ -246,7 +246,7 @@ namespace PPDBM * @param b * @return */ - constexpr bool Ppdbm::equalVector(const CostVector& a, const CostVector& b) + bool Ppdbm::equalVector(const CostVector& a, const CostVector& b) { return a == b; } @@ -258,7 +258,7 @@ namespace PPDBM * @param c2 * @return */ - constexpr relation_t Ppdbm::compareCostLists(const std::vector& c1, const std::vector& c2) + relation_t Ppdbm::compareCostLists(const std::vector& c1, const std::vector& c2) { assert(c1.size() == c2.size()); @@ -284,7 +284,7 @@ namespace PPDBM * @param other * @return */ - constexpr relation_t Ppdbm::relation(const Ppdbm& other) const + relation_t Ppdbm::relation(const Ppdbm& other) const { assert(_dim == other._dim && _param == other._param); @@ -353,7 +353,7 @@ namespace PPDBM * finds the vertices of a dbm! very important function! * @return */ - constexpr std::vector DBMatrix::verticesDBM() const + std::vector DBMatrix::verticesDBM() const { const uint32_t dim = getDim(); std::vector vertices; diff --git a/test/test_ParamPricedDBM.cpp b/test/test_ParamPricedDBM.cpp index ce68fd0..4115b20 100644 --- a/test/test_ParamPricedDBM.cpp +++ b/test/test_ParamPricedDBM.cpp @@ -1,3 +1,5 @@ +/* -*- mode: C++; c-file-style: "stroustrup"; c-basic-offset: 4; indent-tabs-mode: nil; -*- */ + #include "dbm/ParamPricedDBM.h" #include @@ -7,324 +9,302 @@ using namespace PPDBM; -TEST_CASE("DbmBound dominates") +TEST_CASE("Matrix basic operations") { - // valeur plus petite => contrainte plus forte => domine - const auto a = DbmBound{3, false}; // <= 3 - const auto b = DbmBound{5, false}; // <= 5 - CHECK(a.dominates(b)); - CHECK(!b.dominates(a)); - - // meme valeur : strict domine non-strict (< 3 est plus fort que <= 3) - const auto strict3 = DbmBound{3, true}; - const auto loose3 = DbmBound{3, false}; - CHECK(strict3.dominates(loose3)); - CHECK(!loose3.dominates(strict3)); - - // une borne se domine toujours elle-meme - CHECK(a.dominates(a)); - - // INF ne domine jamais une borne finie, l'inverse est toujours vrai - CHECK(!INF_BOUND.dominates(a)); - CHECK(a.dominates(INF_BOUND)); + Matrix m{2, 3, 7}; + CHECK(m.getRows() == 2); + CHECK(m.getCols() == 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); // untouched 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 isIncompatible") +TEST_CASE("DbmBound dominates and isIncompatible") { - // x_i - x_j <= 3 et x_j - x_i <= -4 => x_i - x_j >= 4, incompatible avec <=3 - const auto a = DbmBound{3, false}; - const auto b = DbmBound{-4, false}; - CHECK(DbmBound::isIncompatible(a, b)); - - // x_i - x_j <= 3 et x_j - x_i <= -3 => point unique {x_i - x_j = 3}, compatible - const auto c = DbmBound{-3, false}; - CHECK(!DbmBound::isIncompatible(a, c)); - - // si l'une des deux bornes est stricte, la somme nulle devient incompatible - const auto aStrict = DbmBound{3, true}; - CHECK(DbmBound::isIncompatible(aStrict, c)); - - // deux bornes larges, somme > 0 => compatible - const auto d = DbmBound{1, false}; - const auto e = DbmBound{1, false}; - CHECK(!DbmBound::isIncompatible(d, e)); + 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 operator+") +TEST_CASE("DbmBound addition (transitive composition)") { - const auto a = DbmBound{3, false}; - const auto b = DbmBound{2, false}; + const DbmBound a{3, false}; + const DbmBound b{2, false}; + const auto sum = a + b; CHECK(sum.getValue() == 5); - CHECK(!sum.getStrict()); + CHECK_FALSE(sum.getStrict()); - // si l'une des deux est stricte, le resultat est strict - const auto bStrict = DbmBound{2, true}; - const auto sum2 = a + bStrict; - CHECK(sum2.getValue() == 5); - CHECK(sum2.getStrict()); + const DbmBound strictB{2, true}; + const auto sumStrict = a + strictB; + CHECK(sumStrict.getValue() == 5); + CHECK(sumStrict.getStrict()); // strictness propagates - // INF + n'importe quoi = INF - const auto sum3 = a + INF_BOUND; - CHECK(sum3.isInfinite()); + const auto sumInf = a + INF_BOUND; + CHECK(sumInf.isInfinite()); } -TEST_CASE("DBMatrix construction and reset") +TEST_CASE("DBMatrix initial state") { - auto m = DBMatrix{3}; - REQUIRE(m.getDim() == 3); - - SUBCASE("etat initial : diagonale a 0, reste a INF") - { - for (uint32_t i = 0; i < 3; ++i) { - CHECK(m(i, i).getValue() == 0); - CHECK(!m(i, i).getStrict()); - for (uint32_t j = 0; j < 3; ++j) { - if (i != j) CHECK(m(i, j).isInfinite()); - } - } - } + DBMatrix m{3}; + CHECK(m.getDim() == 3); - SUBCASE("reset revient a l'etat initial") - { - m(0, 1) = DbmBound{7, false}; - m.reset(); - CHECK(m(0, 1).isInfinite()); - CHECK(m(0, 0).getValue() == 0); - CHECK(m(1, 1).getValue() == 0); + 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") +TEST_CASE("DBMatrix close propagates transitive bounds") { - SUBCASE("propagation simple : x1-x0<=3, x2-x1<=2 => x2-x0<=5") - { - auto m = DBMatrix{3}; - m(1, 0) = DbmBound{3, false}; - m.close(1, 0); - m(2, 1) = DbmBound{2, false}; - m.close(2, 1); - - CHECK(m(2, 0).getValue() == 5); - CHECK(!m(2, 0).getStrict()); - } + 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 - SUBCASE("propagation d'une borne stricte") - { - // x1-x0<3 (stricte), x2-x1<=2 (large) => x2-x0<5 (stricte) - auto m = DBMatrix{3}; - m(1, 0) = DbmBound{3, true}; - m.close(1, 0); - m(2, 1) = DbmBound{2, false}; - m.close(2, 1); - - CHECK(m(2, 0).getValue() == 5); - CHECK(m(2, 0).getStrict()); - } -} + CHECK(m(2, 0).isInfinite()); // not propagated yet -TEST_CASE("DBMatrix relation") -{ - // m1: x1<=5, m2: x1<=3 => m2 subset de m1 - auto m1 = DBMatrix{2}; - m1(1, 0) = DbmBound{5, false}; - m1.close(1, 0); - - auto m2 = DBMatrix{2}; - m2(1, 0) = DbmBound{3, false}; - m2.close(1, 0); - - CHECK(m1.relation(m2) == base_SUPERSET); - CHECK(m2.relation(m1) == base_SUBSET); - - auto m3 = DBMatrix{2}; - m3(1, 0) = DbmBound{5, false}; - m3.close(1, 0); - CHECK(m1.relation(m3) == base_EQUAL); + m.close(2, 1); + + CHECK(m(2, 0).getValue() == 5); // x2 <= x1 + 2 <= 5 + CHECK_FALSE(m(2, 0).getStrict()); } -TEST_CASE("DBMatrix verticesDBM") +TEST_CASE("DBMatrix relation between boxes") { - // 2 horloges : x0 = reference (toujours 0), x1 dans [2,5] - auto m = DBMatrix{2}; - m(1, 0) = DbmBound{5, false}; // x1 - x0 <= 5 - m(0, 1) = DbmBound{-2, false}; // x0 - x1 <= -2, donc x1 >= 2 + DBMatrix big{2}; // x1 in [0,5] + big(1, 0) = DbmBound{5, false}; - const auto vertices = m.verticesDBM(); - const auto v25 = Vertex{0, 5}; - const auto v22 = Vertex{0, 2}; + 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}; - REQUIRE(vertices.size() == 2); - CHECK(std::find(vertices.begin(), vertices.end(), v25) != vertices.end()); - CHECK(std::find(vertices.begin(), vertices.end(), v22) != vertices.end()); + CHECK(overlapA.relation(overlapB) == base_DIFFERENT); + CHECK(overlapB.relation(overlapA) == base_DIFFERENT); } -TEST_CASE("Ppdbm reset") +TEST_CASE("DBMatrix verticesDBM on a trivial (single-clock) dimension") { - auto p = Ppdbm{2, 1}; - p.constrain_DBM(1, 0, DbmBound{4, false}); - p.reset(); - CHECK(!p.isEmpty()); - CHECK(p(1, 0).isInfinite()); + DBMatrix m{1}; + const auto vertices = m.verticesDBM(); + REQUIRE(vertices.size() == 1); + CHECK(vertices[0] == Vertex{0}); } -TEST_CASE("Ppdbm constrain_DBM") +TEST_CASE("DBMatrix verticesDBM computes the corners of a box") { - SUBCASE("une contrainte plus faible que l'existante est ignoree") - { - auto p = Ppdbm{2, 0}; - CHECK(p.constrain_DBM(1, 0, DbmBound{5, false})); - CHECK(p.constrain_DBM(1, 0, DbmBound{9, false})); - CHECK(p(1, 0).getValue() == 5); - } + 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 - SUBCASE("une contrainte plus forte resserre la zone") - { - auto p = Ppdbm{2, 0}; - CHECK(p.constrain_DBM(1, 0, DbmBound{5, false})); - CHECK(p.constrain_DBM(1, 0, DbmBound{3, false})); - CHECK(p(1, 0).getValue() == 3); - CHECK(!p.isEmpty()); - } + auto vertices = m.verticesDBM(); - SUBCASE("une contrainte incompatible vide la zone") - { - // x1-x0<=3 puis x0-x1<=-4 (x1>=4) => incompatible avec x1<=3 - auto p = Ppdbm{2, 0}; - CHECK(p.constrain_DBM(1, 0, DbmBound{3, false})); - CHECK(!p.constrain_DBM(0, 1, DbmBound{-4, false})); - CHECK(p.isEmpty()); - } + 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_N") +TEST_CASE("Ppdbm constrain_DBM tightens bounds and detects emptiness") { - SUBCASE("plusieurs contraintes compatibles appliquees d'un coup") - { - // x1 dans [2,5], x2 dans [1,4] - auto p = Ppdbm{3, 0}; - const auto constraints = std::vector{ - {1, 0, DbmBound{5, false}}, - {0, 1, DbmBound{-2, false}}, - {2, 0, DbmBound{4, false}}, - {0, 2, DbmBound{-1, false}}, - }; - CHECK(p.constrain_DBM_N(constraints)); - CHECK(!p.isEmpty()); - CHECK(p(1, 0).getValue() == 5); - CHECK(p(0, 1).getValue() == -2); - CHECK(p(2, 0).getValue() == 4); - CHECK(p(0, 2).getValue() == -1); - } + Ppdbm p{2, 0}; // 1 clock, no parameters - SUBCASE("une contrainte incompatible dans la liste vide la zone") - { - auto p = Ppdbm{2, 0}; - const auto constraints = std::vector{ - {1, 0, DbmBound{3, false}}, - {0, 1, DbmBound{-4, false}}, - }; - CHECK(!p.constrain_DBM_N(constraints)); - CHECK(p.isEmpty()); - } + 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 isDiagonalNegative") +TEST_CASE("Ppdbm constrain_DBM_N applies a batch of constraints atomically") { - auto p = Ppdbm{2, 0}; - CHECK(!p.isDiagonalNegative(0)); - CHECK(!p.isDiagonalNegative(1)); + 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") +{ + Ppdbm p{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(); - p.constrain_DBM(1, 0, DbmBound{3, false}); - p.constrain_DBM(0, 1, DbmBound{-4, false}); - CHECK(p.isDiagonalNegative(1)); + CHECK_FALSE(p.isEmpty()); + CHECK(p(1, 0).isInfinite()); + CHECK(p(0, 1).getValue() == 0); } -TEST_CASE("Ppdbm costAtVertex") +TEST_CASE("Ppdbm costAtVertex reflects rates and the zone's lower bound") { - // 1 parametre. rate(1,0)=2 (constante), rate(1,1)=1 (coeff parametre) - auto p = Ppdbm{2, 1}; - auto rates = Matrix{2, 2, 0}; - rates(1, 0) = 2; + Ppdbm p{2, 1}; // 1 clock, 1 parameter + + Matrix rates{2, 2, 0}; + rates(1, 0) = 1; rates(1, 1) = 1; p.setRates(rates); - p.constrain_DBM(1, 0, DbmBound{5, false}); - p.constrain_DBM(0, 1, DbmBound{-2, false}); // x1 dans [2,5] - const auto v = Vertex{0, 5}; // delta = 5 - 2 = 3 - const auto cost = p.costAtVertex(v); - REQUIRE(cost.size() == 2); - CHECK(cost[0] == 6); // 0 + 2*3 - CHECK(cost[1] == 3); // 0 + 1*3 + 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 dominatesVector and equalVector") +TEST_CASE("Ppdbm costAtOtherOffset translates the cost to another zone's offset") { - const auto a = CostVector{1, 2, 3}; - const auto b = CostVector{2, 2, 4}; - const auto c = CostVector{1, 2, 3}; - const auto d = CostVector{1, 2, 4}; + Ppdbm p{2, 1}; + Matrix rates{2, 2, 0}; + rates(1, 0) = 1; + rates(1, 1) = 1; + p.setRates(rates); + REQUIRE(p.constrain_DBM(1, 0, DbmBound{5, false})); + REQUIRE(p.constrain_DBM(0, 1, DbmBound{-2, false})); // x1 in [2,5] - CHECK(Ppdbm::dominatesVector(a, b)); - CHECK(!Ppdbm::dominatesVector(b, a)); - CHECK(Ppdbm::dominatesVector(a, c)); // egalite => domine (<=) + DBMatrix otherDbm{2}; // x1 in [0,5] : offset is 0 instead of 2 + otherDbm(1, 0) = DbmBound{5, false}; - CHECK(Ppdbm::equalVector(a, c)); - CHECK(!Ppdbm::equalVector(a, d)); + const auto cost = p.costAtOtherOffset(otherDbm); + CHECK(cost == CostVector{0, 0}); } -TEST_CASE("Ppdbm compareCostLists") +TEST_CASE("Ppdbm::dominatesVector and equalVector") { - const auto c1 = std::vector{{1, 1}, {2, 2}}; - const auto c2 = std::vector{{1, 1}, {2, 2}}; - CHECK(Ppdbm::compareCostLists(c1, c2) == base_EQUAL); - - const auto cheaper = std::vector{{1, 1}, {1, 1}}; - const auto costly = std::vector{{2, 2}, {3, 3}}; - // cheaper <= costly partout => cheaper "domine" => base_SUPERSET - CHECK(Ppdbm::compareCostLists(cheaper, costly) == base_SUPERSET); - CHECK(Ppdbm::compareCostLists(costly, cheaper) == base_SUBSET); - - const auto mixed1 = std::vector{{1, 5}, {5, 1}}; - const auto mixed2 = std::vector{{5, 1}, {1, 5}}; - CHECK(Ppdbm::compareCostLists(mixed1, mixed2) == base_DIFFERENT); + 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 relation") +TEST_CASE("Ppdbm::compareCostLists aggregates per-vertex domination") { - SUBCASE("meme zone, couts differents => domination par le cout") - { - auto p1 = Ppdbm{2, 1}; - auto p2 = Ppdbm{2, 1}; - p1.constrain_DBM(1, 0, DbmBound{5, false}); - p1.constrain_DBM(0, 1, DbmBound{-2, false}); - p2.constrain_DBM(1, 0, DbmBound{5, false}); - p2.constrain_DBM(0, 1, DbmBound{-2, false}); - - auto rates1 = Matrix{2, 2, 0}; - rates1(1, 1) = 1; // cout = x1 - p1.setRates(rates1); - - auto rates2 = Matrix{2, 2, 0}; - rates2(1, 1) = 2; // cout = 2*x1, toujours plus cher (ou egal) - p2.setRates(rates2); - - // meme zone, p1 moins cher partout => p1 domine (SUPERSET), p2 est domine (SUBSET) - CHECK(p1.relation(p2) == base_SUPERSET); - CHECK(p2.relation(p1) == base_SUBSET); - } + 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); +} - SUBCASE("zones differentes, le resultat reste stable") - { - auto p1 = Ppdbm{2, 0}; - auto p2 = Ppdbm{2, 0}; - p1.constrain_DBM(1, 0, DbmBound{5, false}); - p2.constrain_DBM(1, 0, DbmBound{3, false}); - p2.constrain_DBM(0, 1, DbmBound{-1, false}); - - const auto r1 = p1.relation(p2); - const auto r2 = p2.relation(p1); - CHECK((r1 == base_SUPERSET || r1 == base_DIFFERENT)); - CHECK((r2 == base_SUBSET || r2 == 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.setRates(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.setRates(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.setRates(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); } \ No newline at end of file From 93114ff5eae4a6e0037d0e2fd9789c281ee39261 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Wed, 26 Aug 2026 15:06:03 +0200 Subject: [PATCH 12/18] post operators written. --- include/dbm/ParamPricedDBM.h | 81 ++++++- src/ParamPricedDBM.cpp | 414 ++++++++++++++++++++++++++++++++++- 2 files changed, 477 insertions(+), 18 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index 7ee9a05..3b9dabc 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -13,8 +13,12 @@ namespace PPDBM { + ///to facilitate notations + using Vertex = std::vector; // size _dim, clock valuations + using CostVector = std::vector; // size _param + 1, coeffs on parameters constexpr auto INF = std::numeric_limits::max() >> 1u; + /// generic matrix type (used everywhere) template struct Matrix { @@ -39,9 +43,26 @@ namespace PPDBM /** * A parameter constraint is of the form (a_0 + sum^n_{i=1} a_i p_i) <= 0 */ - struct ParameterConstraint { std::vector coeffs; }; + struct ParametricConstraint + { + constexpr ParametricConstraint(std::vector coeffs, bool strict): + _coeffs{std::move(coeffs)}, + _strict{strict} {} + constexpr const bool& getStrict() const { return _strict; } + private: + std::vector _coeffs; + bool _strict; + }; - struct Polyhedron { std::vector constraints; }; + 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); } + private: + std::vector _constraints; + }; @@ -97,7 +118,6 @@ namespace PPDBM - using Vertex = std::vector; /** * @@ -121,6 +141,13 @@ namespace PPDBM } } 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; @@ -132,9 +159,15 @@ namespace PPDBM - using Vertex = std::vector; // size _dim, clock valuations - using CostVector = std::vector; // size _param + 1, coeffs on parameters + + struct Facet; + + + + + + //////////////////////////// The PPDBM class that represents Parametric priced zones ////////////////////////////// /** * Parametric priced timed difference bound matrix */ @@ -145,9 +178,15 @@ namespace PPDBM _rates(dim, param + 1, 0) {} 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 setRates(const Matrix& rates) { _rates = rates; } constexpr bool isEmpty() const { return _emptyZone; } void reset(); @@ -156,7 +195,26 @@ namespace PPDBM void addCostAtOffset(); bool isDiagonalNegative(uint32_t i) const; bool constrain_DBM_N(const std::vector& constraints); - std::vector costAtOtherOffset(const DBMatrix& otherDbm) const; + std::vector costAtOtherOffset(const DBMatrix& otherDbm) const; + CostVector getSlope(); + void addConstantCost(const CostVector& q); + + /// Delay operation on a facet put in the form Z \land (clock - otherClock = value). + void delay(uint32_t clock, 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 succesors of the ppdbm by delay. + std::vector post_delta(CostVector p, const std::vector& J); + /// operator post, computes the succesors of the ppdbm by transition e. + std::vector post_e(const std::vector& g, std::optional resetClock, const CostVector& q) const; + /// Is there a domination relation between the two zones? relation_t relation(const Ppdbm& other) const; @@ -174,10 +232,17 @@ namespace PPDBM 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) + 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; + }; } // namespace PPDBM #endif \ No newline at end of file diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index 0590dfe..08119d4 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -9,6 +9,7 @@ #include #include #include +#include namespace PPDBM { @@ -41,7 +42,7 @@ namespace PPDBM _dbm.reset(); std::fill(_offsetCost.begin(), _offsetCost.end(), 0); _rates.fill(0); - _PC.constraints.clear(); + _PC.getConstraints().clear(); _emptyZone = false; } @@ -122,7 +123,7 @@ namespace PPDBM */ bool Ppdbm::constrain_DBM_N(const std::vector& constraints) { - assert(!constraints.empty()); + if (constraints.empty()) return true; // nothing to do --> done // remove the cost at the previous offset removeCostAtOffset(); @@ -156,9 +157,9 @@ namespace PPDBM * @param otherDbm * @return */ - std::vector Ppdbm::costAtOtherOffset(const DBMatrix& otherDbm) const + std::vector Ppdbm::costAtOtherOffset(const DBMatrix& otherDbm) const { - std::vector cost = _offsetCost; + std::vector cost = _offsetCost; for (uint32_t x = 1; x < _dim; ++x) { // difference in x coordinate of the offset @@ -173,10 +174,406 @@ namespace PPDBM return cost; } + /** + * + * @param q + */ + void Ppdbm::addConstantCost(const CostVector& q) + { + for (uint32_t k = 0; k <= _param; ++k) _offsetCost[k] += q[k]; + } - /////////////////////////////////////////////////// relation ///////////////////////////////////////////////////// + + /////////////////////////////////////////////// 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() + { + 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, 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) = DIAG_BOUND; + _dbm(0, clockToReset) = DIAG_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) + { + 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 + 1; ++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, std::optional resetClock, const CostVector& 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) { + // 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); + + { + // cas 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)); + } + } + { + // cas 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; + } + + + + ////////////////////////////////////////////// WQO relation //////////////////////////////////////////////////////// /** @@ -343,10 +740,7 @@ namespace PPDBM - - - -//////////////////////////////////////////////// sommetsDBM() with claude's help ///////////////////////////////////// + //////////////////////////////////////////////// sommetsDBM() with claude's help ///////////////////////////////////// /** @@ -486,7 +880,7 @@ namespace PPDBM }; - // Backtracking : énumère tous les arbres couvrants sur 'candidates' + // 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; From a974bcb51f493ee3be19b0beb99dd425661a814e Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Wed, 26 Aug 2026 15:09:05 +0200 Subject: [PATCH 13/18] previous tests working --- test/test_ParamPricedDBM.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_ParamPricedDBM.cpp b/test/test_ParamPricedDBM.cpp index 4115b20..4f73949 100644 --- a/test/test_ParamPricedDBM.cpp +++ b/test/test_ParamPricedDBM.cpp @@ -218,7 +218,7 @@ TEST_CASE("Ppdbm costAtVertex reflects rates and the zone's lower bound") Matrix rates{2, 2, 0}; rates(1, 0) = 1; rates(1, 1) = 1; - p.setRates(rates); + 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 @@ -236,7 +236,7 @@ TEST_CASE("Ppdbm costAtOtherOffset translates the cost to another zone's offset" Matrix rates{2, 2, 0}; rates(1, 0) = 1; rates(1, 1) = 1; - p.setRates(rates); + 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] @@ -281,7 +281,7 @@ TEST_CASE("Ppdbm::relation compares zones and, when equal, their costs") Matrix ratesA{2, 2, 0}; ratesA(1, 0) = 1; ratesA(1, 1) = 1; - a.setRates(ratesA); + 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] @@ -289,7 +289,7 @@ TEST_CASE("Ppdbm::relation compares zones and, when equal, their costs") Matrix ratesB{2, 2, 0}; ratesB(1, 0) = 2; // strictly more expensive on the first coefficient ratesB(1, 1) = 1; - b.setRates(ratesB); + 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 @@ -302,7 +302,7 @@ TEST_CASE("Ppdbm::relation compares zones and, when equal, their costs") Matrix ratesC{2, 2, 0}; ratesC(1, 0) = 1; ratesC(1, 1) = 1; - c.setRates(ratesC); + 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 From 3dcc570ee9cf9f5ca6db4ea1de23057e9ff7f062 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Thu, 27 Aug 2026 12:33:13 +0200 Subject: [PATCH 14/18] working tests for every function including post operators. --- include/dbm/ParamPricedDBM.h | 3 +- src/ParamPricedDBM.cpp | 2 +- test/test_ParamPricedDBM.cpp | 316 ++++++++++++++++++++++++++++++++++- 3 files changed, 318 insertions(+), 3 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index 3b9dabc..de1b552 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -8,6 +8,7 @@ #include #include #include +#include @@ -196,7 +197,7 @@ namespace PPDBM bool isDiagonalNegative(uint32_t i) const; bool constrain_DBM_N(const std::vector& constraints); std::vector costAtOtherOffset(const DBMatrix& otherDbm) const; - CostVector getSlope(); + CostVector getSlope() const; void addConstantCost(const CostVector& q); /// Delay operation on a facet put in the form Z \land (clock - otherClock = value). diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index 08119d4..272e58c 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -398,7 +398,7 @@ namespace PPDBM * returns sum_{x \in X} r(x) * @return */ - CostVector Ppdbm::getSlope() + CostVector Ppdbm::getSlope() const { CostVector slope(_param + 1,0); for (uint32_t i = 1; i < _dim; ++i) { diff --git a/test/test_ParamPricedDBM.cpp b/test/test_ParamPricedDBM.cpp index 4f73949..ba8d535 100644 --- a/test/test_ParamPricedDBM.cpp +++ b/test/test_ParamPricedDBM.cpp @@ -307,4 +307,318 @@ TEST_CASE("Ppdbm::relation compares zones and, when equal, their costs") // overlapping but neither zone contains the other -> DIFFERENT, regardless of cost CHECK(a.relation(c) == base_DIFFERENT); -} \ No newline at end of file +} + +//////////////////////////////////// Nouvelles fonctions : relax / findLinkedClocks ////////////////////////////////// + +TEST_CASE("DBMatrix relaxDown and relaxUp remove strictness only") +{ + DBMatrix m{3}; + m(0, 1) = DbmBound{-2, true}; // x1 > 2 (borne basse stricte) + m(0, 2) = DbmBound{0, false}; // x2 >= 0 (deja faible) + m(1, 0) = DbmBound{5, true}; // x1 < 5 (borne haute stricte) + m(2, 0) = DbmBound{3, false}; // x2 <= 3 (deja faible) + + m.relaxDown(); + CHECK_FALSE(m(0, 1).getStrict()); + CHECK(m(0, 1).getValue() == -2); // la valeur ne change pas, seule la stricture saute + CHECK_FALSE(m(0, 2).getStrict()); // deja faible : inchange + CHECK(m(1, 0).getStrict()); // relaxDown ne touche pas les bornes hautes + + m.relaxUp(); + CHECK_FALSE(m(1, 0).getStrict()); + CHECK(m(1, 0).getValue() == 5); + CHECK_FALSE(m(2, 0).getStrict()); + + // une borne infinie reste infinie et n'est jamais "relachee" + 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); // relache (i,1) pour tout i != 1 + CHECK_FALSE(m(0, 1).getStrict()); + CHECK_FALSE(m(2, 1).getStrict()); + CHECK(m(1, 0).getStrict()); // pas touche : ce n'est pas une entree (*, 1) + CHECK(m(1, 2).getStrict()); // idem + + m.relaxUpClock(1); // relache (1,i) pour tout 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 exactement : x1-x2<=3 et 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 est lie a 2 + CHECK(next[0] == 0); // la reference n'est liee a personne ici + CHECK(next[2] == 0); // 2 est un representant (aucun j>2 a verifier) + + // si la difference n'est pas exactement nulle, pas de lien + DBMatrix m2{3}; + m2(1, 2) = DbmBound{3, false}; + m2(2, 1) = DbmBound{-2, false}; // x1-x2 dans [2,3], pas fixe + 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 par defaut, aucun lien entre x1 et x2 + + const auto lower = p.lowerFacets(); + REQUIRE(lower.size() == 2); // une facette par horloge, aucune redondance + for (const auto& f : lower) { + CHECK(f._ppdbm(f._constraint.i, 0).getValue() == 0); // chaque facette fixe l'horloge a 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 est implique par x2<=5 et 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); // seule x2 <= 5 est une vraie facette + 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, pas de lien x1/x2 + + const auto lower = p.lowerFacetsRelativeTo(1); + REQUIRE(lower.size() == 1); // seule la reference (0) donne une facette non redondante + CHECK(lower[0]._constraint.i == 1); + CHECK(lower[0]._constraint.j == 0); + CHECK(lower[0]._ppdbm(1, 0).getValue() == 2); // x1 fixe a sa borne basse + + 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 fixe a sa borne haute +} + +////////////////////////////////////////////////////// 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()); + // le taux de clock1 absorbe la difference : la pente totale devient exactement 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); // la diagonale de reference n'est pas touchee + CHECK_FALSE(p(0, 0).getStrict()); + CHECK(p(0, 1).getValue() == 0); // x1 >= 0 toujours vrai +} + +//////////////////////////////////////////////////// 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]; // copie : x1 est deja fixe a 2 sur cette facette + CHECK(f._ppdbm(1, 0).getValue() == 2); + + f._ppdbm.clockReset(1, f._constraint); + + // la zone : x1 est reinitialise a 0, et sa relation a x2 recopie celle de la 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()); + + // le taux de x1 est transfere (ici vers la reference, donc "perdu" pour le calcul du cout futur) + CHECK(f._ppdbm.getSlope() == CostVector{3}); // 5 - 2 : seul rates(x2) reste observable +} + +/////////////////////////////////////////////////////// post_delta ////////////////////////////////////////////////// + +TEST_CASE("Ppdbm post_delta produces the three symbolic successors of the formula") +{ + Ppdbm p{2, 0}; // une seule horloge, pas de parametre (constantes) + 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}, {}); // J vide : teste le chemin constrain_DBM_N(vide) + REQUIRE(result.size() == 3); + + // successeur issu de LF(Z) (cas p <= slope) : x1 devient non borne, la pente devient 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 + + // successeur { (l, Z) } du cas p > slope : zone et pente inchangees + 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 + + // successeur issu de UF(Z) (cas p > slope) : x1 partait de 5, devient non borne, meme pente 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}); // pas de garde, reset de x1, cout d'arc 10 + + REQUIRE(result.size() == 2); // une branche par facette relative a x1 (basse puis haute) + + for (const auto& r : result) { + // dans les deux cas x1 est reinitialise a 0, et la zone ne depend plus que de 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}); // le taux de x1 (2) est retire de la pente observable + CHECK(r.getPC().getConstraints().size() == 1); + } + + // branche issue de lowerFacetsRelativeTo (x1 valait 2 avant reset) : cout = 4 (deja accumule) + 10 + CHECK(result[0].getOffsetCost() == CostVector{14}); + CHECK_FALSE(result[0].getPC().getConstraints()[0].getStrict()); + + // branche issue de upperFacetsRelativeTo (x1 valait 5 avant reset) : cout = 10 (deja accumule) + 10 + CHECK(result[1].getOffsetCost() == CostVector{20}); + CHECK(result[1].getPC().getConstraints()[0].getStrict()); +} + From 5c6d84eb3179ccd6de1b28444344c164234bc4c0 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Fri, 28 Aug 2026 15:28:48 +0200 Subject: [PATCH 15/18] all comments in english now --- test/test_ParamPricedDBM.cpp | 103 +++++++++++++++++------------------ 1 file changed, 51 insertions(+), 52 deletions(-) diff --git a/test/test_ParamPricedDBM.cpp b/test/test_ParamPricedDBM.cpp index ba8d535..3f36de8 100644 --- a/test/test_ParamPricedDBM.cpp +++ b/test/test_ParamPricedDBM.cpp @@ -20,7 +20,7 @@ TEST_CASE("Matrix basic operations") m(1, 2) = 42; CHECK(m(1, 2) == 42); - CHECK(m(0, 0) == 7); // untouched cell + CHECK(m(0, 0) == 7); // unchanged cell m.fill(3); for (uint32_t i = 0; i < 2; ++i) @@ -54,7 +54,7 @@ TEST_CASE("DbmBound dominates and isIncompatible") 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 + // 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})); } @@ -168,7 +168,7 @@ TEST_CASE("Ppdbm constrain_DBM tightens bounds and detects emptiness") 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 + // x1 >= 6 is incompatible with x1 <= 4: the zone becomes empty CHECK_FALSE(p.constrain_DBM(0, 1, DbmBound{-6, false})); CHECK(p.isEmpty()); } @@ -223,7 +223,7 @@ TEST_CASE("Ppdbm costAtVertex reflects rates and the zone's lower bound") 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 + // 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) @@ -309,28 +309,28 @@ TEST_CASE("Ppdbm::relation compares zones and, when equal, their costs") CHECK(a.relation(c) == base_DIFFERENT); } -//////////////////////////////////// Nouvelles fonctions : relax / findLinkedClocks ////////////////////////////////// +//////////////////////////////////// New functions: relax / findLinkedClocks ////////////////////////////////// TEST_CASE("DBMatrix relaxDown and relaxUp remove strictness only") { DBMatrix m{3}; - m(0, 1) = DbmBound{-2, true}; // x1 > 2 (borne basse stricte) - m(0, 2) = DbmBound{0, false}; // x2 >= 0 (deja faible) - m(1, 0) = DbmBound{5, true}; // x1 < 5 (borne haute stricte) - m(2, 0) = DbmBound{3, false}; // x2 <= 3 (deja faible) + 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); // la valeur ne change pas, seule la stricture saute - CHECK_FALSE(m(0, 2).getStrict()); // deja faible : inchange - CHECK(m(1, 0).getStrict()); // relaxDown ne touche pas les bornes hautes + 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()); - // une borne infinie reste infinie et n'est jamais "relachee" + // an infinite bound remains infinite and is never "relaxed" DBMatrix m2{2}; CHECK(m2(1, 0).isInfinite()); m2.relaxUp(); @@ -345,13 +345,13 @@ TEST_CASE("DBMatrix relaxDownClock and relaxUpClock only touch the given clock's m(1, 0) = DbmBound{5, true}; // x1 < 5 m(1, 2) = DbmBound{1, true}; // x1 - x2 < 1 - m.relaxDownClock(1); // relache (i,1) pour tout i != 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()); // pas touche : ce n'est pas une entree (*, 1) - CHECK(m(1, 2).getStrict()); // idem + CHECK(m(1, 0).getStrict()); // untouched: this is not an entry (*, 1) + CHECK(m(1, 2).getStrict()); // same - m.relaxUpClock(1); // relache (1,i) pour tout i != 1 + m.relaxUpClock(1); // relax (1,i) for every i != 1 CHECK_FALSE(m(1, 0).getStrict()); CHECK_FALSE(m(1, 2).getStrict()); } @@ -359,20 +359,20 @@ TEST_CASE("DBMatrix relaxDownClock and relaxUpClock only touch the given clock's TEST_CASE("DBMatrix findLinkedClocks detects zero cycles") { DBMatrix m{3}; - // x1 - x2 = 3 exactement : x1-x2<=3 et x2-x1<=-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 est lie a 2 - CHECK(next[0] == 0); // la reference n'est liee a personne ici - CHECK(next[2] == 0); // 2 est un representant (aucun j>2 a verifier) + 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) - // si la difference n'est pas exactement nulle, pas de lien + // 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 dans [2,3], pas fixe + m2(2, 1) = DbmBound{-2, false}; // x1-x2 in [2,3], not fixed CHECK(m2.findLinkedClocks() == std::vector{0, 0, 0}); } @@ -383,12 +383,12 @@ 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 par defaut, aucun lien entre x1 et x2 + // x1, x2 >= 0 by default, no link between x1 and x2 const auto lower = p.lowerFacets(); - REQUIRE(lower.size() == 2); // une facette par horloge, aucune redondance + REQUIRE(lower.size() == 2); // one facet per clock, no redundancy for (const auto& f : lower) { - CHECK(f._ppdbm(f._constraint.i, 0).getValue() == 0); // chaque facette fixe l'horloge a 0 + CHECK(f._ppdbm(f._constraint.i, 0).getValue() == 0); // each facet fixes the clock to 0 } const auto upper = p.upperFacets(); @@ -404,14 +404,14 @@ TEST_CASE("Ppdbm lowerFacets and upperFacets on an independent box") TEST_CASE("Ppdbm upperFacets discards a redundant facet") { - // x1 <= 5, x2 <= 5, x1 - x2 <= 0 => x1 <= 5 est implique par x2<=5 et x1<=x2 + // 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); // seule x2 <= 5 est une vraie facette + 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); } @@ -423,19 +423,19 @@ TEST_CASE("Ppdbm lowerFacetsRelativeTo and upperFacetsRelativeTo isolate the res 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, pas de lien x1/x2 + 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); // seule la reference (0) donne une facette non redondante + 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 fixe a sa borne basse + 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 fixe a sa borne haute + CHECK(upper[0]._ppdbm(1, 0).getValue() == 5); // x1 fixed to its upper bound } ////////////////////////////////////////////////////// getSlope ///////////////////////////////////////////////////// @@ -469,7 +469,7 @@ TEST_CASE("Ppdbm delay(clock, c) removes upper bounds and rebalances the rate so CHECK(p(1, 0).isInfinite()); CHECK(p(2, 0).isInfinite()); - // le taux de clock1 absorbe la difference : la pente totale devient exactement c + // the rate of clock1 absorbs the difference: the total slope becomes exactly c CHECK(p.getSlope() == CostVector{5, 7}); } @@ -481,9 +481,9 @@ TEST_CASE("Ppdbm delay() only removes upper bounds") p.delay(); CHECK(p(1, 0).isInfinite()); - CHECK(p(0, 0).getValue() == 0); // la diagonale de reference n'est pas touchee + 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 toujours vrai + CHECK(p(0, 1).getValue() == 0); // x1 >= 0 is still true } //////////////////////////////////////////////////// clockReset ///////////////////////////////////////////////////// @@ -503,12 +503,12 @@ TEST_CASE("Ppdbm clockReset applied on a facet from lowerFacetsRelativeTo") const auto facets = p.lowerFacetsRelativeTo(1); REQUIRE(facets.size() == 1); - Facet f = facets[0]; // copie : x1 est deja fixe a 2 sur cette facette + 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); - // la zone : x1 est reinitialise a 0, et sa relation a x2 recopie celle de la reference + // 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); @@ -517,37 +517,37 @@ TEST_CASE("Ppdbm clockReset applied on a facet from lowerFacetsRelativeTo") CHECK(f._ppdbm(2, 1).getValue() == f._ppdbm(2, 0).getValue()); CHECK(f._ppdbm(2, 1).getStrict() == f._ppdbm(2, 0).getStrict()); - // le taux de x1 est transfere (ici vers la reference, donc "perdu" pour le calcul du cout futur) - CHECK(f._ppdbm.getSlope() == CostVector{3}); // 5 - 2 : seul rates(x2) reste observable + // 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}; // une seule horloge, pas de parametre (constantes) + 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}, {}); // J vide : teste le chemin constrain_DBM_N(vide) + const auto result = p.post_delta(CostVector{3}, {}); // empty J: tests the constrain_DBM_N(empty) path REQUIRE(result.size() == 3); - // successeur issu de LF(Z) (cas p <= slope) : x1 devient non borne, la pente devient p + // 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 - // successeur { (l, Z) } du cas p > slope : zone et pente inchangees + // 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 - // successeur issu de UF(Z) (cas p > slope) : x1 partait de 5, devient non borne, meme pente p + // 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); @@ -577,7 +577,7 @@ TEST_CASE("Ppdbm post_e returns nothing when the guard empties the zone") 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 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()); @@ -596,12 +596,12 @@ TEST_CASE("Ppdbm post_e with a reset generates the lower- and upper-facet branch CHECK(p.getOffsetCost() == CostVector{4}); // rates(x1)=2 * lowerbound(x1)=2 - const auto result = p.post_e({}, 1u, CostVector{10}); // pas de garde, reset de x1, cout d'arc 10 + const auto result = p.post_e({}, 1u, CostVector{10}); // no guard, reset x1, edge cost 10 - REQUIRE(result.size() == 2); // une branche par facette relative a x1 (basse puis haute) + REQUIRE(result.size() == 2); // one branch per facet relative to x1 (lower then upper) for (const auto& r : result) { - // dans les deux cas x1 est reinitialise a 0, et la zone ne depend plus que de x2 + // 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); @@ -609,16 +609,15 @@ TEST_CASE("Ppdbm post_e with a reset generates the lower- and upper-facet branch 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}); // le taux de x1 (2) est retire de la pente observable + CHECK(r.getSlope() == CostVector{3}); // the rate of x1 (2) is removed from the observable slope CHECK(r.getPC().getConstraints().size() == 1); } - // branche issue de lowerFacetsRelativeTo (x1 valait 2 avant reset) : cout = 4 (deja accumule) + 10 + // 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()); - // branche issue de upperFacetsRelativeTo (x1 valait 5 avant reset) : cout = 10 (deja accumule) + 10 + // 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()); } - From 7101c530dc358c319212697222100fce4ac746a6 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Sat, 29 Aug 2026 22:57:13 +0200 Subject: [PATCH 16/18] trying to make my own model in order to have some early results for this algorithm --- src/ParamPricedDBM.cpp | 240 +++++++++++++++++++++++--- src/ParametricOptimalReachability.cpp | 74 ++++++++ 2 files changed, 294 insertions(+), 20 deletions(-) create mode 100644 src/ParametricOptimalReachability.cpp diff --git a/src/ParamPricedDBM.cpp b/src/ParamPricedDBM.cpp index 272e58c..9d719af 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -3,16 +3,170 @@ #include "dbm/ParamPricedDBM.h" #include - +#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; + } + + + + + /** * * @param i @@ -178,7 +332,7 @@ namespace PPDBM * * @param q */ - void Ppdbm::addConstantCost(const CostVector& q) + void Ppdbm::addConstantCost(const std::vector& q) { for (uint32_t k = 0; k <= _param; ++k) _offsetCost[k] += q[k]; } @@ -252,7 +406,7 @@ namespace PPDBM // 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; @@ -298,7 +452,8 @@ namespace PPDBM 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 + 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 @@ -403,7 +558,11 @@ namespace PPDBM 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); + if (_rates(i,j) == INF) { + slope[j] = INF; + } else { + slope[j] += _rates(i,j); + } } } return slope; @@ -414,7 +573,7 @@ namespace PPDBM * @param clock * @param c */ - void Ppdbm::delay(uint32_t clock, CostVector c) + void Ppdbm::delay(uint32_t clock, const CostVector& c) { // modify the zone for (uint32_t i = 1; i < getDim(); ++i) @@ -452,8 +611,8 @@ namespace PPDBM assert(clockToReset < _dim); // modify the zone - _dbm(clockToReset, 0) = DIAG_BOUND; - _dbm(0, clockToReset) = DIAG_BOUND; + _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); @@ -472,14 +631,14 @@ namespace PPDBM * @param J * @return */ - std::vector Ppdbm::post_delta(CostVector p, const std::vector& J) + 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 + 1; ++i) { + for (uint32_t i = 0; i <= _param; ++i) { temp[i] = p[i] - temp[i]; } ParametricConstraint pconstraint{temp,false}; @@ -521,7 +680,7 @@ namespace PPDBM * 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, std::optional resetClock, const CostVector& q) const + std::vector Ppdbm::post_e(const std::vector& g, uint32_t resetClock, const std::vector& q) const { std::vector result; @@ -530,19 +689,19 @@ namespace PPDBM if (!g.empty() && !guarded.constrain_DBM_N(g)) return result; // g empties the zone if (guarded.isEmpty()) return result; - if (!resetClock) { + 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; + const uint32_t x = resetClock; CostVector rx(_param + 1); for (uint32_t k = 0; k <= _param; ++k) rx[k] = guarded._rates(x,k); { - // cas r(x) >= 0 <=> -r(x) <= 0 + // 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}; @@ -556,7 +715,7 @@ namespace PPDBM } } { - // cas r(x) < 0 + // case r(x) < 0 ParametricConstraint pconstraint{rx, true}; for (Facet f : guarded.upperFacetsRelativeTo(x)) { @@ -572,6 +731,31 @@ namespace PPDBM } + 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 //////////////////////////////////////////////////////// @@ -670,12 +854,12 @@ namespace PPDBM if (thisLeqOther && otherLeqThis) { return base_EQUAL; } if (thisLeqOther) { return base_SUPERSET; } // *this* cheaper everywhere -> dominates - if (otherLeqThis) { return base_SUBSET; } + if (otherLeqThis) { return base_SUBSET; } // *this* more expensive everywhere -> dominated return base_DIFFERENT; // no domination --> no branch cutting } /** - * most important function, decides when the loops stop / when we cut branches of the exploration graph. + * 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 @@ -734,13 +918,29 @@ namespace PPDBM } } + /** + * + * @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() with claude's help ///////////////////////////////////// + //////////////////////////////////////////////////// sommetsDBM() //////////////////////////////////////////////// /** @@ -902,7 +1102,7 @@ namespace PPDBM vertices.erase(std::unique(vertices.begin(), vertices.end()), vertices.end()); return vertices; } -} // namespace PPBDM +}// namespace PPBDM diff --git a/src/ParametricOptimalReachability.cpp b/src/ParametricOptimalReachability.cpp new file mode 100644 index 0000000..0f77576 --- /dev/null +++ b/src/ParametricOptimalReachability.cpp @@ -0,0 +1,74 @@ +#include "../include/dbm/ParamPricedDBM.h" + +#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{1,0,PPDBM::DbmBound{3,false}}; + auto guards = std::vector{}; + guards.emplace_back(guard); + + // resets (none for the example) + uint32_t clockToReset = 0; + + // edges + auto edge01 = PPDBM::Edge{0,1,clockToReset,std::vector{}}; + auto edge02 = PPDBM::Edge{0,2,clockToReset,guards}; + auto edge12 = PPDBM::Edge{1,2,clockToReset,guards}; + 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(0,PPDBM::Ppdbm{nbClocks + 1,nbParam,true})}; + WAITING[0].second = WAITING[0].second.post_delta(automaton.getCostFunction(0), automaton.getInvariants(0))[0]; + + 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); + } + } + } + + for (auto res : RES) { + res.display(); + } + + return 0; +} + From b487859eb11de692227a8a28b1c10f4d6bb4d33c Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Sat, 29 Aug 2026 22:57:40 +0200 Subject: [PATCH 17/18] trying to make my own model in order to have some early results for this algorithm - 2 --- include/dbm/ParamPricedDBM.h | 126 +++++++++++++++++++++++++++-------- test/test_ParamPricedDBM.cpp | 6 +- 2 files changed, 103 insertions(+), 29 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index de1b552..dfe88ab 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -16,21 +16,21 @@ namespace PPDBM { ///to facilitate notations using Vertex = std::vector; // size _dim, clock valuations - using CostVector = std::vector; // size _param + 1, coeffs on parameters - constexpr auto INF = std::numeric_limits::max() >> 1u; + using CostVector = std::vector; // size _param + 1, coeffs on parameters + constexpr auto INF = std::numeric_limits::max() >> 1u; /// generic matrix type (used everywhere) template struct Matrix { - constexpr Matrix(uint32_t rows, uint32_t cols, const T x): _rows{rows}, _cols{cols}, _data(rows * cols, x){} - constexpr auto getRows() const { return _rows; } - constexpr auto getCols() const { return _cols; } - constexpr auto& getData() const { return _data; } - constexpr auto& getData() { return _data; } - constexpr T& operator()(int i, int j) { return _data[i * _cols + j]; } - constexpr const T& operator()(int i, int j) const { return _data[i * _cols + j]; } - constexpr auto fill(T x){ std::fill(_data.begin(), _data.end(), x); } + 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; @@ -50,6 +50,7 @@ namespace PPDBM _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; @@ -61,6 +62,8 @@ namespace PPDBM constexpr const std::vector& getConstraints() const { return _constraints; } constexpr void addParametricConstraint(const ParametricConstraint& constraint) { _constraints.push_back(constraint); } + + bool isEmpty(); private: std::vector _constraints; }; @@ -97,7 +100,7 @@ namespace PPDBM }; constexpr auto INF_BOUND = DbmBound{INF,true}; - constexpr auto DIAG_BOUND = DbmBound{0,false}; + 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) @@ -127,18 +130,18 @@ namespace PPDBM { constexpr DBMatrix(uint32_t dim) : _matrix{dim, dim, INF_BOUND} { - for (uint32_t i = 0; i < dim; ++i) { _matrix(i, i) = DIAG_BOUND; _matrix(0,i) = DIAG_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.getRows(); } + 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.getRows(); ++i) { - _matrix(i,i)= DIAG_BOUND; - _matrix(0,i) = DIAG_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); @@ -163,7 +166,8 @@ namespace PPDBM struct Facet; - + struct Edge; + struct ParamPricedTimedAutomata; @@ -174,9 +178,16 @@ namespace PPDBM */ struct Ppdbm { - constexpr Ppdbm(uint32_t dim, uint32_t param): _dim{dim}, _param{param}, _dbm{dim}, + Ppdbm(uint32_t dim, uint32_t param, bool init = false): + _dim{dim}, + _param{param}, + _dbm{dim}, _offsetCost(param + 1), - _rates(dim, param + 1, 0) {} + _rates(dim, param + 1, INF) + { + 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; } @@ -188,6 +199,7 @@ namespace PPDBM 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(); @@ -198,10 +210,12 @@ namespace PPDBM bool constrain_DBM_N(const std::vector& constraints); std::vector costAtOtherOffset(const DBMatrix& otherDbm) const; CostVector getSlope() const; - void addConstantCost(const CostVector& q); + 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, CostVector c); + 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); @@ -211,16 +225,17 @@ namespace PPDBM std::vector lowerFacetsRelativeTo(uint32_t clock) const; std::vector upperFacetsRelativeTo(uint32_t clock) const; - /// operator post, computes the succesors of the ppdbm by delay. - std::vector post_delta(CostVector p, const std::vector& J); + /// 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, std::optional resetClock, const CostVector& q) const; - + 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 DBMatrix& dbm) const{ return _dbm.verticesDBM(); }; + 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? @@ -244,6 +259,65 @@ namespace PPDBM Ppdbm _ppdbm; DbmConstraint _constraint; }; + + + //////////////////////////////////////////////// automata model /////////////////////////////////////////////////// + + struct Edge + { + Edge(uint32_t from, uint32_t to, uint32_t clockToReset, const std::vector& guards, + std::vector costFunction = std::vector{}): + _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/test/test_ParamPricedDBM.cpp b/test/test_ParamPricedDBM.cpp index 3f36de8..7d07833 100644 --- a/test/test_ParamPricedDBM.cpp +++ b/test/test_ParamPricedDBM.cpp @@ -12,8 +12,8 @@ using namespace PPDBM; TEST_CASE("Matrix basic operations") { Matrix m{2, 3, 7}; - CHECK(m.getRows() == 2); - CHECK(m.getCols() == 3); + 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); @@ -199,7 +199,7 @@ TEST_CASE("Ppdbm constrain_DBM_N applies a batch of constraints atomically") TEST_CASE("Ppdbm reset returns the zone to its unconstrained state") { - Ppdbm p{2, 0}; + 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()); From 1797815c19f52075fa9255008eec7b129f262489 Mon Sep 17 00:00:00 2001 From: Mathieu Doucen Date: Mon, 31 Aug 2026 11:40:35 +0200 Subject: [PATCH 18/18] small example finally workinggit add .! --- include/dbm/ParamPricedDBM.h | 40 ++++- src/CMakeLists.txt | 3 + src/ParamPricedDBM.cpp | 249 +++++++++++++++++++++++++- src/ParametricOptimalReachability.cpp | 44 ++++- 4 files changed, 319 insertions(+), 17 deletions(-) diff --git a/include/dbm/ParamPricedDBM.h b/include/dbm/ParamPricedDBM.h index dfe88ab..56fea82 100644 --- a/include/dbm/ParamPricedDBM.h +++ b/include/dbm/ParamPricedDBM.h @@ -17,8 +17,20 @@ 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 @@ -64,6 +76,9 @@ namespace PPDBM { _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; }; @@ -151,8 +166,6 @@ namespace PPDBM 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; @@ -183,7 +196,7 @@ namespace PPDBM _param{param}, _dbm{dim}, _offsetCost(param + 1), - _rates(dim, param + 1, INF) + _rates(dim, param + 1, 0) { if (init) for (uint32_t i = 0; i < _dim; ++i) {_dbm(i,0) = ZERO_BOUND;} @@ -244,6 +257,9 @@ namespace PPDBM /// 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 @@ -260,13 +276,27 @@ namespace 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 = std::vector{}): + std::vector costFunction): _from{from}, _to{to}, _clockToReset{clockToReset}, @@ -318,6 +348,8 @@ namespace PPDBM 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 d727c27..430e5e6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 index 9d719af..4ed05eb 100644 --- a/src/ParamPricedDBM.cpp +++ b/src/ParamPricedDBM.cpp @@ -14,7 +14,6 @@ namespace PPDBM { - /** * */ @@ -163,6 +162,49 @@ namespace PPDBM 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); + } + @@ -558,11 +600,7 @@ namespace PPDBM CostVector slope(_param + 1,0); for (uint32_t i = 1; i < _dim; ++i) { for (uint32_t j = 0; j < _param + 1; ++j) { - if (_rates(i,j) == INF) { - slope[j] = INF; - } else { - slope[j] += _rates(i,j); - } + slope[j] += _rates(i,j); } } return slope; @@ -1102,6 +1140,205 @@ namespace PPDBM 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 index 0f77576..9675624 100644 --- a/src/ParametricOptimalReachability.cpp +++ b/src/ParametricOptimalReachability.cpp @@ -2,13 +2,14 @@ #include #include +#include int main() { // the goal of what follows is to modelize the following automaton: --> (0, +p1) ------> (1,+p2) // \ / - // x<=3 \ / x<=3 + // x>=3 \ / x>=3 // \ / // \ / // \/ \/ @@ -19,7 +20,7 @@ int main() uint32_t nbParam = 2; // x - x0 <= 3 - auto guard = PPDBM::DbmConstraint{1,0,PPDBM::DbmBound{3,false}}; + auto guard = PPDBM::DbmConstraint{0,1,PPDBM::DbmBound{-3,false}}; auto guards = std::vector{}; guards.emplace_back(guard); @@ -27,9 +28,10 @@ int main() uint32_t clockToReset = 0; // edges - auto edge01 = PPDBM::Edge{0,1,clockToReset,std::vector{}}; - auto edge02 = PPDBM::Edge{0,2,clockToReset,guards}; - auto edge12 = PPDBM::Edge{1,2,clockToReset,guards}; + 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) @@ -46,8 +48,11 @@ int main() // algo: auto RES = std::vector{}; auto PASSED = std::vector>{}; - auto WAITING = std::vector>{std::make_pair(0,PPDBM::Ppdbm{nbClocks + 1,nbParam,true})}; - WAITING[0].second = WAITING[0].second.post_delta(automaton.getCostFunction(0), automaton.getInvariants(0))[0]; + 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(); @@ -65,9 +70,34 @@ int main() } } + + // 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; }