diff --git a/CMakeLists.txt b/CMakeLists.txt
index 13d46a1..b32352a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -6,6 +6,8 @@ include(GNUInstallDirs)
 option(TESTING OFF)
 option(STATIC OFF)
 
+set(CMAKE_BUILD_TYPE Debug)
+
 cmake_policy(SET CMP0048 NEW) # project() command manages VERSION variables
 set(CMAKE_CXX_STANDARD 17)
 set(UDBM_PACKAGE_STRING "${PACKAGE_NAME} ${PACKAGE_VERSION}")
@@ -14,10 +16,13 @@ set(ENABLE_STORE_MINGRAPH 1)
 CONFIGURE_FILE("src/config.h.cmake" "include/dbm/config.h")
 
 set(CMAKE_PREFIX_PATH "${CMAKE_PREFIX_PATH};${CMAKE_CURRENT_SOURCE_DIR}/libs")
+find_package(doctest 2.4.11 REQUIRED)
 find_package(xxHash 0.8.0 CONFIG REQUIRED)
 find_package(UUtils 1.1.1 REQUIRED COMPONENTS base hash debug)
 
-add_library(UDBM src/DBMAllocator.cpp src/dbm.c src/fed_dbm.cpp src/mingraph.c src/mingraph_read.c src/partition.cpp src/print.cpp src/gen.c src/mingraph_cache.cpp src/mingraph_relation.c src/pfed.cpp src/fed.cpp src/infimum.cpp src/mingraph_equal.c src/mingraph_write.c src/priced.cpp)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/libs)
+
+add_library(UDBM src/DBMAllocator.cpp src/dbm.c src/fed_dbm.cpp src/mingraph.c src/mingraph_read.c src/partition.cpp src/print.cpp src/gen.c src/mingraph_cache.cpp src/mingraph_relation.c src/pfed.cpp src/fed.cpp src/infimum.cpp src/mingraph_equal.c src/mingraph_write.c src/priced.cpp src/cost_type.cpp)
 target_link_libraries(UDBM UUtils::base UUtils::udebug UUtils::hash)
 
 if(STATIC)
@@ -25,6 +30,11 @@ if(STATIC)
 	set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Wl,-Bstatic,--whole-archive -lwinpthread -Wl,--no-whole-archive")
 endif(STATIC)
 
+
+set(CMAKE_BUILD_TYPE Release)
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3")
+
+
 target_include_directories(UDBM
     PRIVATE
         # where the library itself will look for its internal headers
diff --git a/getlibs.sh b/getlibs.sh
index 85b0684..dbc014b 100755
--- a/getlibs.sh
+++ b/getlibs.sh
@@ -35,3 +35,20 @@ cd build
 cmake $CMAKE_ARGS -DCMAKE_INSTALL_PREFIX="$SOURCE_DIR/libs/UUtils" ..
 cmake --build . --config Release
 cmake --install . --config Release
+
+cd $SOURCE_DIR/libs/sources;
+wget -nc https://boostorg.jfrog.io/artifactory/main/release/1.81.0/source/boost_1_81_0.tar.gz
+tar -xvf boost_1_81_0.tar.gz
+mkdir -p "$SOURCE_DIR/libs/sources/boost_1_81_0/build"
+cd "$SOURCE_DIR/libs/sources/boost_1_81_0/build"
+cp -r "$SOURCE_DIR/libs/sources/boost_1_81_0/boost" "$SOURCE_DIR/libs/"
+
+cd $SOURCE_DIR/libs/sources;
+wget https://github.com/doctest/doctest/archive/refs/tags/v2.4.11.tar.gz
+tar -xvf v2.4.11.tar.gz
+mkdir -p "$SOURCE_DIR/libs/sources/doctest-2.4.11/build"
+cd "$SOURCE_DIR/libs/sources/doctest-2.4.11"
+cd build
+cmake $CMAKE_ARGS -DOCTEST_WITH_TESTS=OFF -DDOCTEST_WITH_MAIN_IN_STATIC_LIB=ON -DDOCTEST_USE_STD_HEADERS=ON -DCMAKE_INSTALL_PREFIX="$SOURCE_DIR/libs" ..
+cmake --build . --config Release
+cmake --install . --config Release
diff --git a/include/dbm/cost_type.h b/include/dbm/cost_type.h
new file mode 100644
index 0000000..233ee9c
--- /dev/null
+++ b/include/dbm/cost_type.h
@@ -0,0 +1,18 @@
+//
+// Created by ragusa on 3/17/23.
+//
+
+#ifndef UDBM_COST_TYPE_H
+#define UDBM_COST_TYPE_H
+
+#include <boost/rational.hpp>
+
+/**
+ * The type used to represent cost in the pdbm
+ */
+typedef boost::rational<int> CostType;
+const CostType INFINITE_COST = (CostType)INT_MAX;
+
+uint32_t hash_cost_type(const CostType& val);
+
+#endif  // UDBM_COST_TYPE_H
diff --git a/include/dbm/fed.h b/include/dbm/fed.h
index ac2d7aa..9f9313c 100644
--- a/include/dbm/fed.h
+++ b/include/dbm/fed.h
@@ -955,6 +955,9 @@ namespace dbm
         fed_t& unionWith(fed_t& arg);
         fed_t& unionWithC(fed_t arg);  // dummy wrapper
 
+        // Performs unionWith operation and returns whether the fed changed
+        bool unionWithChanged(const fed_t& arg);
+
         /// Simply add (list concatenation) DBMs to this federation.
         /// @pre same dimension.
 
diff --git a/include/dbm/inline_fed.h b/include/dbm/inline_fed.h
index d164c71..5004393 100644
--- a/include/dbm/inline_fed.h
+++ b/include/dbm/inline_fed.h
@@ -512,6 +512,12 @@ namespace dbm
             return *this;
         }
 
+        bool unionWithChanged(dbmlist_t& arg) {
+            removeIncluded(arg);
+            append(arg);
+            return arg.size() > 0;
+        }
+
         /// Simple reduction by inclusion check of DBMs.
         void reduce(cindex_t dim);
 
diff --git a/include/dbm/partition.h b/include/dbm/partition.h
index 774f996..ae9e813 100644
--- a/include/dbm/partition.h
+++ b/include/dbm/partition.h
@@ -15,6 +15,7 @@
 #define INCLUDE_DBM_PARTITION_H
 
 #include "dbm/fed.h"
+#include "dbm/pfed.h"
 #include "base/intutils.h"
 
 /**
diff --git a/include/dbm/pfed.h b/include/dbm/pfed.h
index 930c8d2..97dee89 100644
--- a/include/dbm/pfed.h
+++ b/include/dbm/pfed.h
@@ -11,14 +11,65 @@
 #define DBM_PFED_H
 
 #include "Valuation.h"
+#include "fed.h"
 
 #include "dbm/priced.h"
+#include "dbm/cost_type.h"
 
 #include <list>
+
 #include <stdexcept>
 
+
+
 namespace dbm
 {
+
+    struct CostPlaneOperation
+    {
+        enum Type {
+            Delay, // delay from a facet of clock i with rate value
+            DelayKeep, // delay but keeping the rates of previous zone, value is delay that was overwritten
+            DiscreteOffset, // increase offset by value
+            ContinuousOffset, // increase offset by value * rate(i)
+            RelativeReset, // reset clock i on a facet relative to clock j
+        };
+        Type type{};
+        cindex_t clock_i;
+        union {
+            cindex_t clock_j;
+            CostType value;
+        };
+        CostPlaneOperation(Type type, cindex_t clock_i, cindex_t clock_j)
+                : type(type), clock_i(clock_i), clock_j(clock_j) {}
+
+        CostPlaneOperation(Type type, cindex_t clock_i, CostType value)
+                : type(type), clock_i(clock_i), value(value) {}
+
+        friend std::ostream& operator<<(std::ostream& os, const CostPlaneOperation& op) {
+            switch (op.type) {
+                case Delay:
+                    os << "↑" << op.clock_i << "," << op.value << "";
+                    break;
+                case DelayKeep:
+                    os << "↑" << op.clock_i << ",?";
+                    break;
+                case DiscreteOffset:
+                    os << "Δ" << op.value << "";
+                    break;
+                case ContinuousOffset:
+                    os << "Δ"  << op.value << "⋅r(" << op.clock_i << ")";
+                    break;
+                case RelativeReset:
+                    os << "↓" << op.clock_i << "," << op.clock_j << "";
+                    break;
+            }
+            return os;
+        }
+    };
+
+    class pdbm_t;
+
     /**
      * Small C++ wrapper for the PDBM priced DBM type.
      *
@@ -41,8 +92,10 @@ namespace dbm
     protected:
         PDBM pdbm;
         cindex_t dim;
-
     public:
+
+        std::vector<CostPlaneOperation> cost_plane_operations;
+
         /**
          * The default constructor constructs an empty priced DBM of
          * dimension zero.  This must not be assigned to another
@@ -92,6 +145,27 @@ namespace dbm
          */
         operator PDBM() const { return pdbm; }
 
+        /**
+         * Returns the cost of the offset point in the pdbm.
+         */
+        [[nodiscard]] CostType getOffsetCost() const { return pdbm_getCostAtOffset(pdbm, dim); };
+
+        [[nodiscard]] CostType getInfimum() const { return pdbm_getInfimum(pdbm, dim); };
+
+        [[nodiscard]] CostType getSupremum() const { return pdbm_getSupremum(pdbm, dim); };
+
+        [[nodiscard]] const CostType* getRates() const { return pdbm_getRates(pdbm, dim); };
+
+        void setRates(CostType* rates) {
+            for (cindex_t x = 1; x < dim; x++) {
+                pdbm_setRate(pdbm, dim, x, rates[x]);
+            }
+        };
+
+        [[nodiscard]] bool canDelay() const;
+
+        bool constrain(cindex_t i, cindex_t j, raw_t c);
+
         /**
          * Returns a hash value for the priced DBM.
          */
@@ -114,6 +188,11 @@ namespace dbm
         /** @see dbm_t::getDBM() */
         raw_t* getDBM();
 
+        /** Converts the pdbm_t to a dbm_t, discarding the cost information */
+        dbm_t to_dbm() const { return dbm_t(const_dbm(), dim); };
+
+        static pdbm_t from_dbm(const dbm_t& dbm) { return pdbm_t(pdbm_from_dbm(dbm.const_dbm(), dbm.getDimension()), dbm.getDimension()); };
+
         /** @see dbm_t::analyseForMinDBM() */
         size_t analyzeForMinDBM(uint32_t* bitMatrix) const;
 
@@ -124,6 +203,31 @@ namespace dbm
         /** @see pdbm_relationWithMinDBM() */
         relation_t relation(const int32_t*, raw_t*) const;
 
+
+        /*
+         * Same as relation
+         */
+        [[nodiscard]] relation_t relation(const pdbm_t& other) const;
+
+        /*
+         * Strict relation, requires cost to be stricly cheaper. Warning return relation may not be as expected.
+         * Returns
+         * (Z,c) ⊃ (Z',c') iff (Z ⊇ Z' ∧ c < c')                        // (Z,c) strictly dominates (Z',c')
+         * (Z,c) ⊂ (Z',c') iff (Z ⊆ Z' ∧ c > c') ∨ (Z ⊂ Z' ∧ c ≥ c')    // (Z',c') weakly dominates (Z,c)
+         * (Z,c) = (Z',c') iff Z = Z' ∧ c = c'                          // (Z,c) = (Z',c')
+         */
+        [[nodiscard]] relation_t strict_relation(const pdbm_t& other) const;
+
+
+        /*
+         * Compares the cost of this pdbm with the cost of the other pdbm, where their zones are identical.
+         * Returns a pair of the relation between the two costs and a strictness.
+         */
+        [[nodiscard]] std::pair<relation_t, strictness_t> compare_to_identical(const pdbm_t& other) const {
+            auto [rel, is_strict] = pdbm_compare_cost_identical_pdbms(pdbm, other.pdbm, dim);
+            return std::make_pair(rel, is_strict ? strictness_t::dbm_STRICT : strictness_t::dbm_WEAK);
+        };
+
         /** @see dbm_t::freeClock() */
         pdbm_t& freeClock(cindex_t clock);
 
@@ -132,6 +236,30 @@ namespace dbm
 
         /** @see dbm_t::readFromMinDBM() */
         static pdbm_t readFromMinDBM(cindex_t dimension, const int32_t*);
+
+        /** Sets the offset to <cost> and all rates to 0 */
+        void setUniformCost(CostType cost);
+
+
+        pdbm_t& operator&=(const dbm::pdbm_t& other);
+
+
+        /// @return string representation of the
+        /// constraints of this PDBM. A clock
+        /// is always positive, so "true" simply means
+        /// all clocks positive.
+        std::string toString(const ClockAccessor&, bool full = false) const;
+
+        /// @return true if it is empty.
+        bool isEmpty() const;
+
+        /// @return dimension
+        cindex_t pdim() const;
+
+        void simplify_rational_cost();
+
+        bool contains(const int32_t* point, cindex_t dim) const;
+        bool contains(const double* point, cindex_t dim) const;
     };
 
     inline pdbm_t::pdbm_t(): pdbm(nullptr), dim(0) {}
@@ -140,9 +268,9 @@ namespace dbm
 
     inline pdbm_t::pdbm_t(PDBM pdbm, cindex_t dim): pdbm(pdbm), dim(dim) { pdbm_incRef(pdbm); }
 
-    inline pdbm_t::pdbm_t(const pdbm_t& other): pdbm(other.pdbm), dim(other.dim) { pdbm_incRef(pdbm); }
+    inline pdbm_t::pdbm_t(const pdbm_t& other): pdbm(other.pdbm), dim(other.dim), cost_plane_operations(other.cost_plane_operations) { pdbm_incRef(pdbm); }
 
-    inline pdbm_t::~pdbm_t() { pdbm_decRef(pdbm); }
+        inline pdbm_t::~pdbm_t() { pdbm_decRef(pdbm); }
 
     inline pdbm_t& pdbm_t::operator=(const pdbm_t& other)
     {
@@ -202,6 +330,7 @@ namespace dbm
         return *this;
     }
 
+
     /**
      * Implementation of priced federations.
      *
@@ -215,9 +344,115 @@ namespace dbm
     class pfed_t
     {
     public:
-        typedef std::list<pdbm_t>::const_iterator const_iterator;
+//
+//        typedef std::list<dbm::pdbm_t>::const_iterator const_list_iterator;
+//        typedef std::list<dbm::pdbm_t>::iterator list_iterator;
+
+        /// Mutable iterator -> iterate though dbm_t
+        class iterator
+        {
+        public:
+            iterator();
+            /// Initialize the iterator of a federation.
+            /// @param _pfed: federation.
+            explicit iterator(pfed_t* _pfed);
+
+            iterator& end();
+
+            /// Dereference to dbm_t, @pre !null()
+            pdbm_t& operator*() const;
+
+            /// Dereference to dbm_t*, @pre !null()
+            pdbm_t* operator->() const;
+
+            /// Mutable access to the matrix as for fed_t, @pre !null()
+            raw_t* operator()() const;
+            raw_t operator()(cindex_t i, cindex_t j) const;
+
+            /// Increment iterator, @pre !null()
+            iterator& operator++();
+
+            // Post increment iterator
+            const iterator operator++(int);
+
+            /// Test if there are DBMs left on the list.
+            bool null() const;
 
-        typedef std::list<pdbm_t>::iterator iterator;
+            /// @return true if there is another DBM after, @pre !null()
+            bool hasNext() const;
+
+            /// Equality test of the internal fdbm_t*
+            bool operator==(const iterator& arg) const;
+            bool operator!=(const iterator& arg) const;
+
+            /// Remove (and deallocate) current dbm_t.
+            void remove();
+
+            /// Remove (and deallocate) current empty dbm_t.
+            void removeEmpty();
+
+            /// Extract the current DBM from the list.
+            /// The result->getNext() points to the rest of the list.
+            pfed_t::iterator extract();
+
+            /// Insert a DBM in the list at the current position.
+            void insert(pdbm_t& pdbm);
+
+            const std::list<pdbm_t>::iterator& getInternalIterator() { return it; };
+
+        private:
+            std::list<pdbm_t>* zones;  /// list of DBMs
+            std::list<pdbm_t>::iterator it;
+            pfed_t* pfed;   /// to update the size
+        };
+
+        /// Const iterator -> iterate though dbm_t
+        class const_iterator
+        {
+        public:
+            const_iterator();
+
+            /// Constructor: @param fed: federation.
+            explicit const_iterator(const pfed_t* pfed);
+
+            // Move iterator to end
+            const_iterator& end();
+
+            /// Dereference to dbm_t
+            const pdbm_t& operator*() const;
+
+            /// Dereference to dbm_t*, @pre !null()
+            const pdbm_t* operator->() const;
+
+            /// Access to the matrix as for fed_t
+            const raw_t* operator()() const;
+            raw_t operator()(cindex_t i, cindex_t j) const;
+
+            // New iterator at end
+
+            /// Increment iterator, @pre !null()
+            const_iterator& operator++();
+
+            /// Post increment iterator
+            const const_iterator operator++(int);
+
+            /// Test if there are DBMs left on the list.
+            bool null() const;
+
+            /// @return true if there is another DBM after, @pre !null()
+            bool hasNext() const;
+
+            /// Equality test of the internal fdbm_t*
+            bool operator==(const const_iterator& arg) const;
+            bool operator!=(const const_iterator& arg) const;
+
+            const std::list<pdbm_t>::const_iterator& getInternalIterator() { return it; };
+
+        private:
+            const std::list<pdbm_t>* zones;  /// list of DBMs
+            std::list<pdbm_t>::const_iterator it;
+            const pfed_t* pfed;   /// to update the size
+        };
 
     protected:
         struct pfed_s
@@ -244,18 +479,19 @@ namespace dbm
         /** Copy-on-write: Creates an unshared copy of the federation. */
         void cow();
 
+    public:
+        /// Initialize a pfed_t to empty federation of a given dimension.
+        /// @param dim: dimension of the federation.
+        /// @post isEmpty()
+        explicit pfed_t(cindex_t dim = 1);
+
         /**
          * Adds \a pdbm to the federation. The reference count on pdbm
          * is incremented by one.
          */
         void add(const PDBM pdbm, cindex_t dim);
 
-    public:
-        /** Allocate empty priced federation of dimension 0. */
-        pfed_t();
-
-        /** Allocate empty priced federation of dimension \a dim. */
-        explicit pfed_t(cindex_t dim);
+        void add(pdbm_t pdbm);
 
         /**
          * Allocate a priced federation of dimension \a dim initialised to
@@ -272,6 +508,18 @@ namespace dbm
          */
         ~pfed_t();
 
+        static pfed_t from_fed(const fed_t& fed) {
+            pfed_t pfed(fed.getDimension());
+            for (const dbm::dbm_t& dbm : fed) {
+                pfed.add(pdbm_t::from_dbm(dbm));
+            }
+            return pfed;
+        }
+
+        const pdbm_t& first() {
+            return ptr->zones.front();
+        }
+
         /**
          * Constrain x(i) to value.
          *
@@ -314,7 +562,16 @@ namespace dbm
         bool constrain(const constraint_t& c);
 
         /** Returns the infimum of the federation. */
-        int32_t getInfimum() const;
+        CostType getInfimum() const;
+
+        /** Returns the supremum of the federation */
+        CostType getSupremum() const;
+
+
+        bool canDelay() const;
+
+        /** Sets the offset to <cost> and all rates to 0 */
+        void setUniformCost(CostType cost);
 
         /**
          * Check if the federation satisfies a given constraint.
@@ -337,7 +594,7 @@ namespace dbm
         bool isUnbounded() const;
 
         /** Contains a hash of the federation. */
-        uint32_t hash(uint32_t seed) const;
+        uint32_t hash(uint32_t seed = 0) const;
 
         /** Returns true iff the federation contains \v. */
         bool contains(const DoubleValuation& v) const;
@@ -345,6 +602,15 @@ namespace dbm
         /** Returns true iff the federation contains \v. */
         bool contains(const IntValuation& v) const;
 
+        /// Not implemented
+        bool contains(const int32_t* point, cindex_t dim) const;
+        bool contains(const double* point, cindex_t dim) const;
+
+        pfed_t& predt(const pfed_t& bad, const raw_t* restrict = NULL);
+
+
+        [[nodiscard]] fed_t to_fed() const;
+
         /**
          * Returns true iff the federation contains \v, ignoring
          * strictness of constraints.
@@ -352,31 +618,36 @@ namespace dbm
         bool containsWeakly(const IntValuation& v) const;
 
         /** Delay with the current delay rate. */
-        void up();
+        pfed_t& up();
 
         /** Delay with rate \a rate. */
-        void up(int32_t rate);
+        pfed_t& up(CostType rate);
 
         /** Set x(clock) to \a value. */
-        void updateValue(cindex_t clock, uint32_t value);
+        pfed_t& updateValue(cindex_t clock, uint32_t value);
 
         /**
          * Set x(clock) to \a value. x(zero) must be on a zero-cycle with
          * x(clock).
          */
-        void updateValueZero(cindex_t clock, int32_t value, cindex_t zero);
+        pfed_t& updateValueZero(cindex_t clock, int32_t value, cindex_t zero);
 
         void extrapolateMaxBounds(int32_t* max);
         void diagonalExtrapolateMaxBounds(int32_t* max);
         void diagonalExtrapolateLUBounds(int32_t* lower, int32_t* upper);
 
-        void incrementCost(int32_t value);
+        void incrementCost(CostType value);
+
+        CostType getCostOfValuation(const IntValuation& valuation) const;
+
 
-        int32_t getCostOfValuation(const IntValuation& valuation) const;
         void relax();
 
         void freeClock(cindex_t clock);
 
+        // Simplifies the rational cost to avoid large numbers, no real guarantees of what the cost is afterwards
+        void simplify_rational_cost();
+
         /**
          * Resets the federation to the federation containing only the
          * origin, with a cost of 0.
@@ -419,7 +690,7 @@ namespace dbm
          * Erases a DBM from the federation and returns an iterator to
          * the successor element.
          */
-        iterator erase(iterator);
+        iterator erase(iterator& iter);
 
         /** Assignment operator. */
         pfed_t& operator=(const pfed_t&);
@@ -442,9 +713,15 @@ namespace dbm
         /** Not implemented. */
         pfed_t& operator&=(const pfed_t&);
 
+        /** Not implemented. */
+        pfed_t& operator&=(const pdbm_t&);
+
         /** Not implemented. */
         bool intersects(const pfed_t&) const;
 
+        /** Equals Operator */
+        bool operator==(const pfed_t&) const;
+
         /** Not implemented. */
         pfed_t operator&(const pfed_t& b) const;
 
@@ -452,13 +729,24 @@ namespace dbm
         pfed_t operator-(const pfed_t& b) const;
 
         /** Not implemented. */
-        void down();
+        pfed_t& down();
+
+        /** Resize all the DBMs of this federation, @see dbm_t.
+         * @see dbm_shrinkExpand in dbm.h.
+         * @param bitSrc,bitDst,bitSize: bit strings of (int) size bitSize
+         * that mark the source and destination active clocks.
+         * @param table: redirection table to write.
+         * @pre bitSrc & bitDst are uint32_t[bitSize] and
+         * table is a uint32_t[32*bitSize]
+         * @post the indirection table is written.
+         */
+        void resize(const uint32_t* bitSrc, const uint32_t* bitDst, size_t bitSize, cindex_t* table);
 
         /** Not implemented. */
         int32_t getUpperMinimumCost(cindex_t) const;
 
         /** Not implemented. */
-        void relaxUp();
+        pfed_t& relaxUp();
 
         /** Not implemented. */
         void getValuation(double* cval, size_t dim, bool* freeC = nullptr) const;
@@ -481,6 +769,96 @@ namespace dbm
         /** Not implemented. */
         pfed_t& convexHull();
 
+        /** Not implemented
+         * @return the max upper bound (raw) of a clock.
+         */
+        raw_t getMaxUpper(cindex_t) const;
+
+        /** Not implemented
+         * @return the max lower bound (raw) of a clock.
+         */
+        raw_t getMaxLower(cindex_t) const;
+
+        /*
+         * Prints string
+         */
+        std::string toString(const ClockAccessor&, bool full = false) const;
+
+        /// Not implemented
+        pfed_t& upStop(const uint32_t* stopped);
+
+        /// Not implemented
+        void updateClock(cindex_t x, cindex_t y);
+
+        /// Not implemented
+        pfed_t& append(pfed_t& arg);
+
+        /// Not implemented
+        pfed_t& steal(pfed_t& arg);
+
+        /// Not implemented
+        pfed_t& toLowerBounds() const;
+
+        /// Not implemented
+        pfed_t& toUpperBounds() const;
+
+        /// Not implemented
+        bool isConstrainedBy(cindex_t, cindex_t, raw_t) const;
+
+        /// Not implemented
+        bool getDelay(const double* point, cindex_t dim, double* min, double* max, double* minVal = NULL,
+                      bool* minStrict = NULL, double* maxVal = NULL, bool* maxStrict = NULL,
+                      const uint32_t* stopped = NULL) const;
+
+        /// Not implemented
+        pfed_t& mergeReduce(size_t skip = 0, int level = 0);
+
+        /// Not implemented
+        bool le(const pfed_t& arg) const;
+
+        /// Not implemented
+        bool lt(const pfed_t& arg) const;
+
+        /// Not implemented
+        bool gt(const pfed_t& arg) const;
+
+        /// Not implemented
+        bool ge(const pfed_t& arg) const;
+
+        /// Not implemented
+        bool eq(const pfed_t& arg) const;
+
+        void clear_cost_plane_operations() {
+            cow();
+            std::for_each(ptr->zones.begin(), ptr->zones.end(), [](pdbm_t& z) { z.cost_plane_operations.clear(); });
+        }
+
+
+        /*
+         * Non-exact union. Keeps all zones that are not subseteq of another zone. Returns true if the zones changed.
+         */
+        bool unionWith(const pfed_t& arg);
+
+        // Performs unionWith operation and returns whether the fed changed
+        bool unionWithChanged(const pfed_t& arg) { return unionWith(pfed_t(arg)); };
+
+        /// Not implemented
+        int32_t maxOnZero(cindex_t clock);
+
+        /// Not implemented
+        pfed_t& downStop(const uint32_t* stopped);
+
+        /// Not implemented
+        bool hasZero() const;
+
+        /// Not implemented
+        void swap(pfed_t&);
+
+        /// Not implemented
+        void intern();
+
+        void nil();
+
         /**
          * Relation between two priced federations: SUBSET is returned
          * if all zones of this federation are contained in some zone
@@ -505,7 +883,7 @@ namespace dbm
          * Returns the infimum of the federation given a partial
          * valuation.
          */
-        int32_t getInfimumValuation(IntValuation& valuation, const bool* free = nullptr) const;
+        CostType getInfimumValuation(IntValuation& valuation, const bool* free = nullptr) const;
 
         /** Returns the number of zones in the federation. */
         size_t size() const;
@@ -582,21 +960,23 @@ namespace dbm
         }
     }
 
+
+
     inline pfed_t::iterator pfed_t::beginMutable()
     {
         prepare();
-        return ptr->zones.begin();
+        return iterator(this);
     }
 
     inline pfed_t::iterator pfed_t::endMutable()
     {
         prepare();
-        return ptr->zones.end();
+        return iterator(this).end();
     }
 
-    inline pfed_t::const_iterator pfed_t::begin() const { return ptr->zones.begin(); }
+    inline pfed_t::const_iterator pfed_t::begin() const { return const_iterator(this); }
 
-    inline pfed_t::const_iterator pfed_t::end() const { return ptr->zones.end(); }
+    inline pfed_t::const_iterator pfed_t::end() const { return const_iterator(this).end(); }
 
     inline cindex_t pfed_t::getDimension() const { return ptr->dim; }
 
@@ -634,6 +1014,190 @@ namespace dbm
     inline bool pfed_t::satisfies(const constraint_t& c) const { return satisfies(c.i, c.j, c.value); }
 
     std::ostream& operator<<(std::ostream&, const pfed_t&);
+
+    pfed_t operator|(const pfed_t& a, const pfed_t& b);
+
+    /***********************************************
+     *  Inlined implementations of pfed_t::iterator *
+     ***********************************************/
+
+    inline pfed_t::iterator::iterator(pfed_t* _pfed): zones(&_pfed->ptr->zones), it(_pfed->ptr->zones.begin()), pfed(_pfed) {
+        assert(_pfed);
+    }
+
+    inline pfed_t::iterator::iterator(): zones(nullptr), pfed(nullptr) {}
+
+    inline pdbm_t& pfed_t::iterator::operator*() const
+    {
+        assert(pfed && zones);
+        return *it;
+    }
+
+    inline pdbm_t* pfed_t::iterator::operator->() const
+    {
+        assert(pfed && zones);
+        return &*it;
+    }
+
+    inline raw_t* pfed_t::iterator::operator()() const
+    {
+        assert(pfed && zones);
+        return it->getDBM();
+    }
+
+    inline raw_t pfed_t::iterator::operator()(cindex_t i, cindex_t j) const
+    {
+        assert(pfed && zones);
+        return (*it)(i, j);
+    }
+
+    inline pfed_t::iterator& pfed_t::iterator::operator++()
+    {
+        assert(pfed && zones);
+        ++it;
+        return *this;
+    }
+
+    inline const pfed_t::iterator pfed_t::iterator::operator++(int)
+    {
+        assert(pfed && zones);
+        auto before = iterator(*this);
+        ++it;
+        return before;
+    }
+
+    inline bool pfed_t::iterator::null() const
+    {
+        return pfed == nullptr || zones == nullptr || it == zones->end();
+    }
+
+    inline bool pfed_t::iterator::hasNext() const
+    {
+        assert(pfed && zones);
+        return it != zones->end();
+    }
+
+    inline bool pfed_t::iterator::operator==(const iterator& arg) const
+    {
+        assert(pfed && zones);
+        return it == arg.it;
+    }
+
+    inline bool pfed_t::iterator::operator!=(const iterator& arg) const { return !(*this == arg); }
+
+    inline void pfed_t::iterator::remove()
+    {
+        assert(pfed && zones);
+        it = zones->erase(it);
+    }
+
+    inline void pfed_t::iterator::removeEmpty()
+    {
+        //        assert(fdbm && *fdbm);
+        //        *fdbm = (*fdbm)->removeEmptyAndNext();
+        //        ifed->decSize();
+        throw std::logic_error("pfed_t::iterator::removeEmpty not implemented");
+    }
+
+    inline pfed_t::iterator pfed_t::iterator::extract()
+    {
+        //        assert(fdbm && *fdbm);
+        //        fpdbm_t* current = *fdbm;
+        //        *fdbm = current->getNext();
+        //        ifed->decSize();
+        //        return current;
+        throw std::logic_error("pfed_t::iterator::extract not implemented");
+    }
+
+    inline void pfed_t::iterator::insert(pdbm_t& pdbm)
+    {
+        assert(pfed->getDimension() == pdbm.getDimension());
+        pfed->prepare();
+        it = zones->insert(it, pdbm);
+    }
+
+    inline pfed_t::iterator& pfed_t::iterator::end() {
+        assert(pfed && zones);
+        it = zones->end();
+        return *this;
+    }
+
+    inline pfed_t::iterator pfed_t::erase(iterator& iter)
+    {
+        iter.remove();
+        return iter;
+    }
+
+    /*****************************************************
+     *  Inlined implementations of pfed_t::const_iterator *
+     *****************************************************/
+
+    inline pfed_t::const_iterator::const_iterator(const pfed_t* _pfed) : zones(&_pfed->ptr->zones), it(_pfed->ptr->zones.begin()), pfed(_pfed) {
+        assert(_pfed && zones);
+    }
+
+    inline pfed_t::const_iterator::const_iterator(): zones(nullptr), pfed(nullptr) {}
+
+    inline const pdbm_t& pfed_t::const_iterator::operator*() const
+    {
+        assert(pfed && zones);
+        return *it;
+    }
+
+    inline const pdbm_t* pfed_t::const_iterator::operator->() const
+    {
+        assert(pfed && zones);
+        return &*it;
+    }
+
+    inline const raw_t* pfed_t::const_iterator::operator()() const
+    {
+        assert(pfed && zones);
+        return it->const_dbm();
+    }
+
+    inline raw_t pfed_t::const_iterator::operator()(cindex_t i, cindex_t j) const
+    {
+        assert(pfed && zones);
+        return (*it)(i, j);
+    }
+
+    inline pfed_t::const_iterator& pfed_t::const_iterator::operator++()
+    {
+        assert(pfed && zones);
+        ++it;
+        return *this;
+    }
+
+    inline const pfed_t::const_iterator pfed_t::const_iterator::operator++(int)
+    {
+        assert(pfed && zones);
+        auto before = const_iterator(*this);
+        ++it;
+        return before;
+    }
+
+    inline bool pfed_t::const_iterator::null() const
+    {
+        return pfed == nullptr || zones == nullptr || it == zones->end();
+    }
+
+    inline bool pfed_t::const_iterator::hasNext() const
+    {
+        assert(pfed && zones);
+        return it != zones->end();
+    }
+
+    inline bool pfed_t::const_iterator::operator==(const const_iterator& arg) const { return it == arg.it; }
+
+    inline bool pfed_t::const_iterator::operator!=(const const_iterator& arg) const { return !(*this == arg); }
+
+    inline pfed_t::const_iterator& pfed_t::const_iterator::end()
+    {
+        assert(pfed && zones);
+        it = zones->end();
+        return *this;
+    }
 }  // namespace dbm
 
 #endif /* DBM_PFED_H */
diff --git a/include/dbm/priced.h b/include/dbm/priced.h
index 40d0b5d..9867fa6 100644
--- a/include/dbm/priced.h
+++ b/include/dbm/priced.h
@@ -15,6 +15,7 @@
 #include "dbm/dbm.h"
 #include "dbm/mingraph.h"
 #include "base/inttypes.h"
+#include "dbm/cost_type.h"
 
 #include <iosfwd>
 #include <utility>
@@ -245,6 +246,20 @@ bool pdbm_constrainToFacet(PDBM& pdbm, cindex_t dim, cindex_t i, cindex_t j);
  */
 relation_t pdbm_relation(const PDBM pdbm1, const PDBM pdbm2, cindex_t dim);
 
+relation_t pdbm_relation_strict(const PDBM pdbm1, const PDBM pdbm2, cindex_t dim);
+
+/*
+ * Compares the costs of two identical DBMs.
+ * If the zones are not identical, it returns {base_DIFFERENT, false}
+ *
+ * @see    relation_t
+ * @param  pdbm1       is a closed priced DBM of dimension \a dim.
+ * @param  pdbm2       is a closed priced DBM of dimension \a dim.
+ * @param  dim         is the dimension of \a pdbm1 and \a pdbm2.
+ * @return The cost relation between pdbm1 and pdbm2 and if the relation is strict.
+ */
+std::pair<relation_t, bool> pdbm_compare_cost_identical_pdbms(const PDBM pdbm1, const PDBM 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
@@ -267,7 +282,17 @@ relation_t pdbm_relationWithMinDBM(const PDBM pdbm, cindex_t dim, const mingraph
  * @param  dim  is the dimension of \a pdbm.
  * @return The infimum cost of \a pdbm.
  */
-int32_t pdbm_getInfimum(const PDBM pdbm, cindex_t dim);
+CostType pdbm_getInfimum(const PDBM pdbm, cindex_t dim);
+
+/**
+ * Computes the supremum cost of the priced DBM.
+ * This does not cache the result and it is recomputed every time!
+ *
+ * @param  pdbm is a closed priced DBM of dimension \a dim
+ * @param  dim  is the dimension of \a pdbm.
+ * @return The supremum cost of \a pdbm.
+ */
+CostType pdbm_getSupremum(const PDBM pdbm, cindex_t dim);
 
 /**
  * Generates a valuation which has the infimum cost of the priced DBM.
@@ -291,7 +316,7 @@ int32_t pdbm_getInfimum(const PDBM pdbm, cindex_t dim);
  * @throw out_of_range if no valuation with the given constraints can
  * be found.
  */
-int32_t pdbm_getInfimumValuation(const PDBM pdbm, cindex_t dim, int32_t* valuation, const bool* free);
+CostType pdbm_getInfimumValuation(const PDBM pdbm, cindex_t dim, int32_t* valuation, const bool* free);
 
 /**
  * Check if a priced DBM satisfies a given constraint.
@@ -322,6 +347,16 @@ bool pdbm_isEmpty(const PDBM pdbm, cindex_t dim);
  */
 bool pdbm_isUnbounded(const PDBM pdbm, cindex_t dim);
 
+
+/**
+ * Sets the cost of a priced zone to be uniform
+ *
+ * @param pdbm is a priced DBM of \a dim
+ * @param dim is the dimension of \a pdbm
+ * @param is the cost to set
+ */
+void pdbm_setUniformCost(PDBM& pdbm, cindex_t dim, CostType cost);
+
 /**
  * Compute a hash value for a priced DBM.
  *
@@ -387,7 +422,7 @@ void pdbm_up(PDBM& pdbm, cindex_t dim);
  *              reference clock.
  * @post  The priced DBM is closed.
  */
-void pdbm_upZero(PDBM& pdbm, cindex_t dim, int32_t rate, cindex_t zero);
+void pdbm_upZero(PDBM& pdbm, cindex_t dim, CostType rate, cindex_t zero);
 
 /**
  * Updates \a clock to \a value. This is only legitimate if the
@@ -458,7 +493,7 @@ void pdbm_diagonalExtrapolateLUBounds(PDBM& pdbm, cindex_t dim, int32_t* lower,
  * @post  The priced DBM is closed.
  * @pre   value >= 0
  */
-void pdbm_incrementCost(PDBM& pdbm, cindex_t dim, int32_t value);
+void pdbm_incrementCost(PDBM& pdbm, cindex_t dim, CostType value);
 
 /**
  * Compute the closure of a priced DBM. This function is only relevant
@@ -571,7 +606,7 @@ bool pdbm_findNextZeroCycle(const PDBM pdbm, cindex_t dim, cindex_t x, cindex_t*
  * @param pdbm  is a closed priced DBM of dimension \a dim.
  * @param dim   is the dimension of \a pdbm.
  */
-int32_t pdbm_getSlopeOfDelayTrajectory(const PDBM pdbm, cindex_t dim);
+CostType pdbm_getSlopeOfDelayTrajectory(const PDBM pdbm, cindex_t dim);
 
 /**
  * Returns the rate (coefficient of the hyperplane) of \a clock.
@@ -581,9 +616,9 @@ int32_t pdbm_getSlopeOfDelayTrajectory(const PDBM pdbm, cindex_t dim);
  * @param  clock is the clock for which to return the coefficient.
  * @return the rate of \a clock.
  */
-int32_t pdbm_getRate(const PDBM pdbm, cindex_t dim, cindex_t clock);
+CostType pdbm_getRate(const PDBM pdbm, cindex_t dim, cindex_t clock);
 
-const int32_t* pdbm_getRates(const PDBM pdbm, cindex_t dim);
+const CostType* pdbm_getRates(const PDBM pdbm, cindex_t dim);
 
 /**
  * Returns the cost of the offset point.
@@ -591,7 +626,7 @@ const int32_t* pdbm_getRates(const PDBM pdbm, cindex_t dim);
  * @param pdbm  is a closed priced DBM of dimension \a dim.
  * @param dim   is the dimension of \a pdbm.
  */
-uint32_t pdbm_getCostAtOffset(const PDBM pdbm, cindex_t dim);
+CostType pdbm_getCostAtOffset(const PDBM pdbm, cindex_t dim);
 
 /**
  * Sets the cost at the offset point.
@@ -605,7 +640,7 @@ uint32_t pdbm_getCostAtOffset(const PDBM pdbm, cindex_t dim);
  * @param dim   is the dimension of \a pdbm.
  * @param value is the new cost of the offset point.
  */
-void pdbm_setCostAtOffset(PDBM& pdbm, cindex_t dim, uint32_t value);
+void pdbm_setCostAtOffset(PDBM& pdbm, cindex_t dim, CostType value);
 
 /**
  * Returns true if the DBM is valid. Useful for debugging.
@@ -674,7 +709,7 @@ uint32_t pdbm_getUpperFacets(PDBM& pdbm, cindex_t dim, cindex_t* facets);
  * @pre    pdbm_containsInt(pdbm, dim, valuation)
  * @return The cost of \a valuation in \a pdbm.
  */
-int32_t pdbm_getCostOfValuation(const PDBM pdbm, cindex_t dim, const int32_t* valuation);
+CostType pdbm_getCostOfValuation(const PDBM pdbm, cindex_t dim, const int32_t* valuation);
 
 /**
  * Makes all strong constraints of a priced DBM weak.
@@ -703,7 +738,7 @@ void pdbm_getOffset(const PDBM pdbm, cindex_t dim, int32_t* valuation);
  * @param clock is the index of a clock for which to set the coefficient.
  * @param rate  is the coefficient.
  */
-void pdbm_setRate(PDBM& pdbm, cindex_t dim, cindex_t clock, int32_t rate);
+void pdbm_setRate(PDBM& pdbm, cindex_t dim, cindex_t clock, CostType rate);
 
 /**
  * Returns the inner matrix of a priced DBM. The matrix can be
@@ -788,6 +823,18 @@ bool pdbm_hasNormalForm(PDBM pdbm, cindex_t dim);
  */
 void pdbm_normalise(PDBM pdbm, cindex_t dim);
 
+PDBM pdbm_from_dbm(const int32_t* dbm, cindex_t dim);
+
+/**
+ * Calculates the intersection of dst and src and saves it in dst using the cost plane of dst and ignoring the cost plane of src.
+ */
+void pdbm_intersect(PDBM& dst, const PDBM& src, cindex_t dim);
+
+/**
+ * Simplifies the rational cost values in the pdbm by approximating them as integers
+ */
+void pdbm_simplify_rational_cost(PDBM& pdbm, cindex_t dim);
+
 ///////////////////////////////////////////////////////////////////////////
 
 #endif /* INCLUDE_DBM_PRICED_H */
diff --git a/src/cost_type.cpp b/src/cost_type.cpp
new file mode 100644
index 0000000..52d1b68
--- /dev/null
+++ b/src/cost_type.cpp
@@ -0,0 +1,9 @@
+//
+// Created by ragusa on 4/5/23.
+//
+
+#include "dbm/cost_type.h"
+
+uint32_t hash_cost_type(const CostType& val) {
+    return ((uint32_t)val.numerator() << 16) ^ (uint32_t)val.denominator();
+}
diff --git a/src/fed.cpp b/src/fed.cpp
index 271ce2a..7239644 100644
--- a/src/fed.cpp
+++ b/src/fed.cpp
@@ -1374,6 +1374,22 @@ namespace dbm
         return *this;
     }
 
+    bool fed_t::unionWithChanged(const fed_t& arg) {
+        assert(getDimension() == arg.getDimension());
+        assert(isOK());
+        assert(arg.isOK());
+
+        if (!arg.isEmpty()) {
+            setMutable();
+            auto copy = fed_t(arg);
+            copy.setMutable();
+            bool changed = ifed()->unionWithChanged(*copy.ifed());
+            copy.ifed()->reset();
+            return changed;
+        }
+        return false;
+    }
+
     fed_t& fed_t::append(fed_t& arg)
     {
         assert(getDimension() == arg.getDimension());
diff --git a/src/infimum.cpp b/src/infimum.cpp
index 6e73506..0eeeb26 100644
--- a/src/infimum.cpp
+++ b/src/infimum.cpp
@@ -8,6 +8,7 @@
 
 #include <stdexcept>
 #include <cstdlib>
+#include <limits>
 
 /*
  * constraintValue(i,j)
@@ -40,7 +41,7 @@ struct Node
     uint32_t depth;    /**< Depth from root (Always node0). */
     Node* thread;      /**< Next according to the preorder. */
     bool inbound;      /**< Is the arc pointing to me or away. */
-    int32_t flow;      /**< Flow for the current solution. */
+    CostType flow;      /**< Flow for the current solution. */
     int32_t potential; /**< Node potential for spanning tree solution. */
 };
 
@@ -61,14 +62,14 @@ static void printSimpleNodeInfo(Node* nodes, uint32_t i)
               << ",  inbound: " << inbound(i) << ", flow: " << flow(i) << ", potential: " << potential(i);
 }
 
-static void printNodeInfo(Node* nodes, const int32_t* rates, uint32_t i)
+static void printNodeInfo(Node* nodes, const CostType* rates, uint32_t i)
 {
     printSimpleNodeInfo(nodes, i);
     std::cerr << ", s/d: " << b(i) << " ";
     std::cerr << std::endl;
 }
 
-static void printAllNodeInfo(Node* nodes, const int32_t* rates, uint32_t dim)
+static void printAllNodeInfo(Node* nodes, const CostType* rates, uint32_t dim)
 {
     std::cerr << "Nodes:" << std::endl;
     for (uint32_t i = 0; i < dim; i++) {
@@ -87,11 +88,11 @@ static void printSimpleAllNodeInfo(Node *nodes, uint32_t dim)
 }
 */
 
-static bool checkTreeIntegrity(const raw_t* dbm, const int32_t* rates, Node* nodes, uint32_t dim)
+static bool checkTreeIntegrity(const raw_t* dbm, const CostType* rates, Node* nodes, uint32_t dim)
 {
     bool treeOK = true;
     // Check that all nodes get their flow
-    int32_t sum[dim], total = 0;
+    CostType sum[dim], total = 0;
     uint32_t i;
     sum[0] = 0;
 
@@ -203,7 +204,7 @@ static bool checkTreeIntegrity(const raw_t* dbm, const int32_t* rates, Node* nod
 
 static inline bool nodeHasNoChildren(Node* node) { return (node->thread->depth > node->depth); }
 
-static void printSolution(int32_t* valuation, const int32_t* rates, uint32_t dim)
+static void printSolution(int32_t* valuation, const CostType* rates, uint32_t dim)
 {
     for (uint32_t i = 0; i < dim; i++) {
         std::cerr << "Val " << i << ": " << valuation[i] << " (rate: " << rates[i] << ")" << std::endl;
@@ -249,7 +250,7 @@ static void printAllRates(const int32_t *rates, uint32_t dim)
  * The initial solution includes computing the intial node potentials.
  *
  */
-static void findInitialSpanningTreeSolution(const raw_t* dbm, uint32_t dim, const int32_t* rates, Node* nodes)
+static void findInitialSpanningTreeSolution(const raw_t* dbm, uint32_t dim, const CostType* rates, Node* nodes)
 {
     Node* last = nodes + dim;
     Node* node = nodes;
@@ -362,8 +363,7 @@ static inline bool isPredecessorOf(Node* n, Node* m) { return n == findNthPredec
  * I.e., pred, depth, and thread.
  *
  */
-static void updateNonRootSubtree(Node* rootNode, Node* nonRootNode, Node* leave, bool sourceInRootSubtree,
-                                 uint32_t flow)
+static void updateNonRootSubtree(Node* rootNode, Node* nonRootNode, Node* leave, bool sourceInRootSubtree, CostType flow)
 {
     /*
      * Update thread information NOT COMPLETELY SURE ABOUT CORRECTNESS
@@ -427,7 +427,7 @@ static void updateNonRootSubtree(Node* rootNode, Node* nonRootNode, Node* leave,
      */
 
     Node *tmppred1, *tmppred2, *newi;
-    int32_t tmpflow1, tmpflow2;
+    CostType tmpflow1, tmpflow2;
     bool tmpinbound1, tmpinbound2;
     tmppred1 = rootNode;
     tmpflow1 = flow;
@@ -475,7 +475,7 @@ static void updateNonRootSubtree(Node* rootNode, Node* nonRootNode, Node* leave,
  * and subtract flow from arcs in the opposite direction
  * of (k,l).
  */
-static void updateFlowInCycle(Node* k, Node* l, Node* root, int32_t flowToAugment)
+static void updateFlowInCycle(Node* k, Node* l, Node* root, CostType flowToAugment)
 {
     if (flowToAugment > 0) {
         while (k != root) {
@@ -498,7 +498,7 @@ static void updateFlowInCycle(Node* k, Node* l, Node* root, int32_t flowToAugmen
 static void updateSpanningTree(Node* k, Node* l, Node* leave, Node* root, int32_t costEnter)
 {
     int32_t reducedCostEnter = costEnter - k->potential + l->potential;
-    int32_t flowToAugment = leave->flow;
+    CostType flowToAugment = leave->flow;
     updateFlowInCycle(k, l, root, flowToAugment);
     if (!isPredecessorOf(leave, k)) {
         updatePotentials(leave, -reducedCostEnter);
@@ -598,7 +598,7 @@ static Node* discoverCycleRoot(Node* k, Node* l)
  */
 static Node* findLeavingArc(Node* k, Node* l, Node* root)
 {
-    int32_t smallestFlow = INT_MAX;
+    CostType smallestFlow = INFINITE_COST;
     /*
      * Node with the lowest depth of the leaving arc.
      */
@@ -631,7 +631,7 @@ static Node* findLeavingArc(Node* k, Node* l, Node* root)
         }
         l = l->pred;
     }
-    ASSERT(smallestFlow != INT_MAX, std::cerr << "No oppositely directed arcs" << std::endl);
+    ASSERT(smallestFlow != INFINITE_COST, std::cerr << "No oppositely directed arcs" << std::endl);
 
     return smallestFlowNode;
 }
@@ -645,14 +645,14 @@ static Node* findLeavingArc(Node* k, Node* l, Node* root)
  * If it is not a transshipment node, we need to follow the
  * thread to update the potentials. We do it by computing
  */
-static void testAndRemoveArtificialArcs(const raw_t* dbm, const int32_t* rates, Node* nodes, uint32_t dim)
+static void testAndRemoveArtificialArcs(const raw_t* dbm, const CostType* rates, Node* nodes, uint32_t dim)
 {
     for (uint32_t i = 1; i < dim; i++) {
         if (potential(i) == dbm_INFINITY && pred(i) == 0 && flow(i) == 0) {
             inbound(i) = true;
 
             Node* tmp = nodes[i].thread;
-            int32_t rateSum = b(i);
+            CostType rateSum = b(i);
 
             int32_t minPotential = dbm_INFINITY + constraintValue(0, i);
 
@@ -688,7 +688,7 @@ static void testAndRemoveArtificialArcs(const raw_t* dbm, const int32_t* rates,
  * Returns true if and only if all elements in [first, last) are
  * non-negative.
  */
-static bool allPositive(const int32_t* first, const int32_t* last)
+static bool allPositive(const CostType* first, const CostType* last)
 {
     while (first != last) {
         if (*first < 0) {
@@ -704,7 +704,7 @@ static bool allPositive(const int32_t* first, const int32_t* last)
  * network flow problem. I.e. the cost of the offset and the costless
  * moves to the offset is not taken into account.
  */
-static void infimumNetSimplex(const raw_t* dbm, uint32_t dim, const int32_t* rates, Node* nodes)
+static void infimumNetSimplex(const raw_t* dbm, uint32_t dim, const CostType* rates, Node* nodes)
 {
     /* Find and store the minimal set of arcs. The netsimplex
      * algorithm is polynomial in the number of arcs.
@@ -766,7 +766,7 @@ static void infimumNetSimplex(const raw_t* dbm, uint32_t dim, const int32_t* rat
  * the infimum-achieving point in the zone.
  *
  */
-int32_t pdbm_infimum(const raw_t* dbm, uint32_t dim, uint32_t offsetCost, const int32_t* rates)
+CostType pdbm_infimum(const raw_t* dbm, uint32_t dim, CostType offsetCost, const CostType* rates)
 {
     if (allPositive(rates, rates + dim)) {
         return offsetCost;
@@ -775,7 +775,7 @@ int32_t pdbm_infimum(const raw_t* dbm, uint32_t dim, uint32_t offsetCost, const
     Node nodes[dim];
     infimumNetSimplex(dbm, dim, rates, nodes);
 
-    int32_t solution = offsetCost;
+    CostType solution = offsetCost;
     for (uint32_t i = 0; i < dim; i++) {
         ASSERT(potential(i) >= 0, std::cerr << "Node: " << i << std::endl; printAllNodeInfo(nodes, rates, dim);
                printClockLowerBounds(dbm, dim));
@@ -784,7 +784,7 @@ int32_t pdbm_infimum(const raw_t* dbm, uint32_t dim, uint32_t offsetCost, const
          * if so, return minus infinity as the solution is unbounded.
          */
         if (potential(i) == dbm_INFINITY && pred(i) == 0 && flow(i) > 0) {
-            return -dbm_INFINITY;
+            return -INFINITE_COST;
         }
         solution += rates[i] * (potential(i) + constraintValue(0, i));
     }
@@ -792,7 +792,7 @@ int32_t pdbm_infimum(const raw_t* dbm, uint32_t dim, uint32_t offsetCost, const
     return solution;
 }
 
-void pdbm_infimum(const raw_t* dbm, uint32_t dim, uint32_t offsetCost, const int32_t* rates, int32_t* valuation)
+void pdbm_infimum(const raw_t* dbm, uint32_t dim, CostType offsetCost, const CostType* rates, int32_t* valuation)
 {
     if (allPositive(rates, rates + dim)) {
         valuation[0] = 0;
diff --git a/src/infimum.h b/src/infimum.h
index 126e720..876b9c2 100644
--- a/src/infimum.h
+++ b/src/infimum.h
@@ -12,7 +12,11 @@
 
 #include "dbm/constraints.h"
 
-int32_t pdbm_infimum(const raw_t* dbm, uint32_t dim, uint32_t offsetCost, const int32_t* rates);
-void pdbm_infimum(const raw_t* dbm, uint32_t dim, uint32_t offsetCost, const int32_t* rates, int32_t* valuation);
+#include <cmath>
+#include "dbm/cost_type.h"
+
+
+CostType pdbm_infimum(const raw_t* dbm, uint32_t dim, CostType offsetCost, const CostType* rates);
+void pdbm_infimum(const raw_t* dbm, uint32_t dim, CostType offsetCost, const CostType* rates, int32_t* valuation);
 
 #endif
diff --git a/src/pfed.cpp b/src/pfed.cpp
index 6862c4c..472bd5a 100644
--- a/src/pfed.cpp
+++ b/src/pfed.cpp
@@ -13,6 +13,7 @@
 #include <functional>  // std::bind
 #include <numeric>
 #include <stdexcept>
+#include <float.h>
 
 using std::bind;
 using std::for_each;
@@ -45,14 +46,6 @@ namespace dbm
         ptr->count = 1;
     }
 
-    pfed_t::pfed_t()
-    {
-        ptr = alloc.allocate(1);
-        alloc.construct(ptr, pfed_s());
-        ptr->dim = 0;
-        ptr->count = 1;
-    }
-
     pfed_t::pfed_t(cindex_t dim)
     {
         ptr = alloc.allocate(1);
@@ -70,12 +63,6 @@ namespace dbm
         add(pdbm, dim);
     }
 
-    pfed_t::iterator pfed_t::erase(iterator i)
-    {
-        assert(ptr->count <= 1);
-        return ptr->zones.erase(i);
-    }
-
 #define for_each__(Z)     for (iterator Z = beginMutable(), Z_ = endMutable(); Z != Z_; zone++)
 #define for_each_(Z)      for (iterator Z = beginMutable(), Z_ = endMutable(); Z != Z_;)
 #define for_each_const(Z) for (const_iterator Z = begin(), Z_ = end(); Z != Z_; zone++)
@@ -89,37 +76,106 @@ namespace dbm
 
     bool pfed_t::constrain(cindex_t i, cindex_t j, raw_t constraint)
     {
-        // lambda x . pdbm_constrain1(x, ptr->dim, i, j, constraint)
-        erase_if_not(bind(pdbm_constrain1, _1, ptr->dim, i, j, constraint));
+//        // lambda x . pdbm_constrain1(x, ptr->dim, i, j, constraint)
+//        erase_if_not(bind(pdbm_constrain1, _1, ptr->dim, i, j, constraint));
+
+        auto it = beginMutable();
+        auto end = endMutable();
+        for (; it != end;) {
+            int prev_offset[ptr->dim];
+            pdbm_getOffset(*it, ptr->dim, prev_offset);
+
+            if (!pdbm_constrain1(*it, ptr->dim, i, j, constraint)) {
+                it = erase(it);
+                continue;
+            }
+
+            int new_offset[ptr->dim];
+            pdbm_getOffset(*it, ptr->dim, new_offset);
+            for (int x = 1; x < ptr->dim; x++) {
+                if (prev_offset[x] != new_offset[x]) {
+                    it->cost_plane_operations.emplace_back(CostPlaneOperation::Type::ContinuousOffset, x, CostType(new_offset[x] - prev_offset[x]));
+                }
+            }
+            it++;
+        }
         return !isEmpty();
     }
 
     bool pfed_t::constrain(const constraint_t* constraints, size_t n)
     {
-        erase_if_not(bind(pdbm_constrainN, _1, ptr->dim, constraints, n));
+//        erase_if_not(bind(pdbm_constrainN, _1, ptr->dim, constraints, n));
+
+        auto it = beginMutable();
+        auto end = endMutable();
+        for (; it != end;) {
+            int prev_offset[ptr->dim];
+            pdbm_getOffset(*it, ptr->dim, prev_offset);
+
+            if (!pdbm_constrainN(*it, ptr->dim, constraints, n)) {
+                it = erase(it);
+                continue;
+            }
+
+            int new_offset[ptr->dim];
+            pdbm_getOffset(*it, ptr->dim, new_offset);
+            for (int x = 1; x < ptr->dim; x++) {
+                if (prev_offset[x] != new_offset[x]) {
+                    it->cost_plane_operations.emplace_back(CostPlaneOperation::Type::ContinuousOffset, x, CostType(new_offset[x] - prev_offset[x]));
+                }
+            }
+            it++;
+        }
         return !isEmpty();
     }
 
-    static int32_t min(int32_t a, int32_t b) { return a < b ? a : b; }
+    static CostType min(CostType a, CostType b) { return a < b ? a : b; }
+    static CostType max(CostType a, CostType b) { return a > b ? a : b; }
 
-    int32_t pfed_t::getInfimum() const
+    CostType pfed_t::getInfimum() const
     {
-        /* The binary operator used in the following accumulation is
-         * lambda inf zone.min(inf, pdbm_getInfimum(zone, ptr->dim));
-         */
-        return accumulate(begin(), end(), INT_MAX, bind(min, _1, bind(pdbm_getInfimum, _2, ptr->dim)));
+        CostType inf = INFINITE_COST;
+        for (const auto& zone : *this) {
+            inf = min(inf, pdbm_getInfimum(zone, ptr->dim));
+        }
+        return inf;
     }
 
-    int32_t pfed_t::getInfimumValuation(IntValuation& valuation, const bool* free) const
+    CostType pfed_t::getSupremum() const
+    {
+        CostType sup = -INFINITE_COST;
+        for (const auto& zone : *this) {
+            sup = max(sup, pdbm_getSupremum(zone, ptr->dim));
+        }
+        return sup;
+    }
+
+
+    bool pfed_t::canDelay() const {
+        for (const auto& zone : *this) {
+            if (zone.canDelay()) {
+                return true;
+            }
+        }
+        return false;
+    }
+
+    void pfed_t::setUniformCost(CostType cost) {
+        for (auto it = beginMutable(); it != endMutable(); ++it) {
+            it->setUniformCost(cost);
+        }
+    }
+
+    CostType pfed_t::getInfimumValuation(IntValuation& valuation, const bool* free) const
     {
         uint32_t dim = ptr->dim;
-        int32_t infimum = INT_MAX;
+        CostType infimum = INFINITE_COST;
         int32_t copy[dim];
 
         for_each_const(zone)
         {
             std::copy(valuation.begin(), valuation.end(), copy);
-            int32_t inf = pdbm_getInfimumValuation(*zone, dim, copy, free);
+            CostType inf = pdbm_getInfimumValuation(*zone, dim, copy, free);
             if (inf < infimum) {
                 infimum = inf;
                 std::copy(copy, copy + dim, valuation.begin());
@@ -130,10 +186,10 @@ namespace dbm
 
     bool pfed_t::satisfies(cindex_t i, cindex_t j, raw_t constraint) const
     {
-        return find_if(begin(), end(), bind(pdbm_satisfies, _1, ptr->dim, i, j, constraint)) != end();
+        return find_if(ptr->zones.begin(), ptr->zones.end(), bind(pdbm_satisfies, _1, ptr->dim, i, j, constraint)) != ptr->zones.end();
     }
 
-    bool pfed_t::isUnbounded() const { return find_if(begin(), end(), bind(pdbm_isUnbounded, _1, ptr->dim)) != end(); }
+    bool pfed_t::isUnbounded() const { return find_if(ptr->zones.begin(), ptr->zones.end(), bind(pdbm_isUnbounded, _1, ptr->dim)) != ptr->zones.end(); }
 
     uint32_t pfed_t::hash(uint32_t seed) const
     {
@@ -145,17 +201,17 @@ namespace dbm
 
     bool pfed_t::contains(const IntValuation& valuation) const
     {
-        return find_if(begin(), end(), bind(pdbm_containsInt, _1, ptr->dim, valuation())) != end();
+        return find_if(ptr->zones.begin(), ptr->zones.end(), bind(pdbm_containsInt, _1, ptr->dim, valuation())) != ptr->zones.end();
     }
 
     bool pfed_t::contains(const DoubleValuation& valuation) const
     {
-        return find_if(begin(), end(), bind(pdbm_containsDouble, _1, ptr->dim, valuation())) != end();
+        return find_if(ptr->zones.begin(), ptr->zones.end(), bind(pdbm_containsDouble, _1, ptr->dim, valuation())) != ptr->zones.end();
     }
 
     bool pfed_t::containsWeakly(const IntValuation& valuation) const
     {
-        return find_if(begin(), end(), bind(pdbm_containsIntWeakly, _1, ptr->dim, valuation())) != end();
+        return find_if(ptr->zones.begin(), ptr->zones.end(), bind(pdbm_containsIntWeakly, _1, ptr->dim, valuation())) != ptr->zones.end();
     }
 
     relation_t pfed_t::relation(const pfed_t& b) const
@@ -219,12 +275,13 @@ namespace dbm
         return relation(b);
     }
 
-    void pfed_t::up() { for_each(beginMutable(), endMutable(), bind(pdbm_up, _1, ptr->dim)); }
+    pfed_t& pfed_t::up() {
+        for_each(beginMutable(), endMutable(), bind(pdbm_up, _1, ptr->dim));
+        return *this;
+    }
 
-    void pfed_t::up(int32_t rate)
+    pfed_t& pfed_t::up(CostType rate)
     {
-        assert(rate >= 0);
-
         uint32_t dim = ptr->dim;
 
         prepare();
@@ -232,17 +289,36 @@ namespace dbm
         for_each__(zone)
         {
             cindex_t x;
-            int32_t oldrate = pdbm_getSlopeOfDelayTrajectory(*zone, dim);
-
-            if (rate == oldrate) {
-                /* This is a simple case which does not require splits.
-                 */
-                pdbm_up(*zone, dim);
-            } else if (pdbm_findZeroCycle(*zone, dim, 0, &x)) {
+            CostType oldrate = pdbm_getSlopeOfDelayTrajectory(*zone, dim);
+            int prev_offset[ptr->dim];
+            pdbm_getOffset(*zone, ptr->dim, prev_offset);
+
+//            pdbm_print(std::cout, zone->operator PDBM(), dim);
+
+            auto constrain_and_up = [&](pdbm_t& z, cindex_t i, cindex_t j, cindex_t zero) {
+                pdbm_constrainToFacet(z, dim, i, j);
+                int new_offset[ptr->dim];
+                pdbm_getOffset(z, ptr->dim, new_offset);
+                for (int x = 1; x < ptr->dim; x++) {
+                    if (prev_offset[x] != new_offset[x]) {
+                        z.cost_plane_operations.emplace_back(CostPlaneOperation::Type::ContinuousOffset, x, CostType(new_offset[x] - prev_offset[x]));
+                    }
+                }
+                pdbm_upZero(z, dim, rate, zero);
+                z.cost_plane_operations.emplace_back(CostPlaneOperation::Delay, zero, rate);
+            };
+
+//            if (rate == oldrate) {
+//                /* This is a simple case which does not require splits.
+//                 */
+//                pdbm_up(*zone, dim);
+//            } else
+            if (pdbm_findZeroCycle(*zone, dim, 0, &x)) {
                 /* This is a simple case which does not require splits.
                  */
                 pdbm_upZero(*zone, dim, rate, x);
-            } else if (rate < oldrate) {
+                zone->cost_plane_operations.emplace_back(CostPlaneOperation::Delay, x, rate);
+            } else if (rate <= oldrate) {
                 /* New rate is smaller than old rate, so we delay from the
                  * lower facets
                  */
@@ -251,11 +327,9 @@ namespace dbm
                 assert(count > 0 && count <= dim);
                 for (uint32_t j = 0; j < count - 1; j++) {
                     ptr->zones.push_front(*zone);
-                    pdbm_constrainToFacet(ptr->zones.front(), dim, 0, facets[j]);
-                    pdbm_upZero(ptr->zones.front(), dim, rate, facets[j]);
+                    constrain_and_up(ptr->zones.front(), 0, facets[j], facets[j]);
                 }
-                pdbm_constrainToFacet(*zone, dim, 0, facets[count - 1]);
-                pdbm_upZero(*zone, dim, rate, facets[count - 1]);
+                constrain_and_up(*zone, 0, facets[count - 1], facets[count - 1]);
             } else {
                 assert(rate > oldrate);
 
@@ -267,21 +341,23 @@ namespace dbm
                 assert(count > 0 && count <= dim);
                 for (uint32_t j = 0; j < count; j++) {
                     ptr->zones.push_front(*zone);
-                    pdbm_constrainToFacet(ptr->zones.front(), dim, facets[j], 0);
-                    pdbm_upZero(ptr->zones.front(), dim, rate, facets[j]);
+                    constrain_and_up(ptr->zones.front(), facets[j], 0, facets[j]);
                 }
-                pdbm_constrainToFacet(*zone, dim, facets[count - 1], 0);
-                pdbm_upZero(*zone, dim, rate, facets[count - 1]);
+                zone->cost_plane_operations.emplace_back(CostPlaneOperation::DelayKeep, 0, rate);
+//                pdbm_constrainToFacet(*zone, dim, facets[count - 1], 0);
+//                pdbm_upZero(*zone, dim, rate, facets[count - 1]);
             }
         }
+        return *this;
     }
 
-    void pfed_t::updateValueZero(cindex_t clock, int32_t value, cindex_t zero)
+    pfed_t& pfed_t::updateValueZero(cindex_t clock, int32_t value, cindex_t zero)
     {
         for_each(beginMutable(), endMutable(), bind(pdbm_updateValueZero, _1, ptr->dim, clock, value, zero));
+        return *this;
     }
 
-    void pfed_t::updateValue(cindex_t clock, uint32_t value)
+    pfed_t& pfed_t::updateValue(cindex_t clock, uint32_t value)
     {
         uint32_t dim = ptr->dim;
 
@@ -290,26 +366,48 @@ namespace dbm
         for_each__(zone)
         {
             cindex_t k;
-            int32_t rate = pdbm_getRate(*zone, dim, clock);
+            CostType rate = pdbm_getRate(*zone, dim, clock);
+
+            uint32_t facets[dim];
+            int prev_offset[ptr->dim];
+            pdbm_getOffset(*zone, ptr->dim, prev_offset);
+
+            // Constrain to the facet i - j < m where zero is the clock in the constraint that is not clock; reset clock
+            auto constrain_and_reset = [&](pdbm_t& z, uint32_t i, uint32_t j, uint32_t zero) {
+                pdbm_constrainToFacet(z, dim, i, j);
+                int new_offset[ptr->dim];
+                pdbm_getOffset(z, ptr->dim, new_offset);
+                for (int x = 1; x < ptr->dim; x++) {
+                    if (prev_offset[x] != new_offset[x]) {
+                      z.cost_plane_operations.emplace_back(CostPlaneOperation::Type::ContinuousOffset, x, CostType(new_offset[x] - prev_offset[x]));
+                    }
+                }
+                pdbm_updateValueZero(z, dim, clock, value, zero);
+                if (0 != zero) {
+                    z.cost_plane_operations.emplace_back(CostPlaneOperation::Type::RelativeReset, clock, zero);
+                }
+            };
 
-            if (rate == 0) {
-                pdbm_updateValue(*zone, dim, clock, value);
-            } else if (pdbm_findZeroCycle(*zone, dim, clock, &k)) {
+            // Remove this because we cannot know which facet was reset
+//            if (rate == 0) {
+//                pdbm_updateValue(*zone, dim, clock, value);
+//            } else
+            if (pdbm_findZeroCycle(*zone, dim, clock, &k)) {
                 pdbm_updateValueZero(*zone, dim, clock, value, k);
+                if (0 != k) {
+                    zone->cost_plane_operations.emplace_back(CostPlaneOperation::RelativeReset, clock, k);
+                }
             } else if (rate > 0) {
                 /* Find and reset lower facets when rate is positive.
                  */
-                uint32_t facets[dim];
                 int32_t count = pdbm_getLowerRelativeFacets(*zone, dim, clock, facets);
 
                 assert(count >= 1);
                 for (int32_t j = 0; j < count - 1; j++) {
                     ptr->zones.push_front(*zone);
-                    pdbm_constrainToFacet(ptr->zones.front(), dim, facets[j], clock);
-                    pdbm_updateValueZero(ptr->zones.front(), dim, clock, value, facets[j]);
+                    constrain_and_reset(ptr->zones.front(), facets[j], clock, facets[j]);
                 }
-                pdbm_constrainToFacet(*zone, dim, facets[count - 1], clock);
-                pdbm_updateValueZero(*zone, dim, clock, value, facets[count - 1]);
+                constrain_and_reset(*zone, facets[count - 1], clock, facets[count - 1]);
             } else {
                 /* Find and reset upper facets when rate is negative.
                  */
@@ -319,13 +417,12 @@ namespace dbm
                 assert(count >= 1);
                 for (int32_t j = 0; j < count - 1; j++) {
                     ptr->zones.push_front(*zone);
-                    pdbm_constrainToFacet(ptr->zones.front(), dim, clock, facets[j]);
-                    pdbm_updateValueZero(ptr->zones.front(), dim, clock, value, facets[j]);
+                    constrain_and_reset(ptr->zones.front(), clock, facets[j], facets[j]);
                 }
-                pdbm_constrainToFacet(*zone, dim, clock, facets[count - 1]);
-                pdbm_updateValueZero(*zone, dim, clock, value, facets[count - 1]);
+                constrain_and_reset(*zone, clock, facets[count - 1], facets[count - 1]);
             }
         }
+        return *this;
     }
 
     void pfed_t::extrapolateMaxBounds(int32_t* max)
@@ -343,18 +440,17 @@ namespace dbm
         for_each(beginMutable(), endMutable(), bind(pdbm_diagonalExtrapolateLUBounds, _1, ptr->dim, lower, upper));
     }
 
-    void pfed_t::incrementCost(int32_t value)
+    void pfed_t::incrementCost(CostType value)
     {
+        for (auto& pdbm : ptr->zones) {
+            pdbm.cost_plane_operations.emplace_back(CostPlaneOperation::DiscreteOffset, 0, value);
+        }
         for_each(beginMutable(), endMutable(), bind(pdbm_incrementCost, _1, ptr->dim, value));
     }
 
-    int32_t pfed_t::getCostOfValuation(const IntValuation& valuation) const
+    CostType pfed_t::getCostOfValuation(const IntValuation& valuation) const
     {
-        /* The binary operator used in the following accumulation is
-         * lambda x y . min(x, pdbm_getCostOfValuation(y, dim, val))
-         */
-        return accumulate(begin(), end(), INT_MAX,
-                          bind(min, _1, bind(pdbm_getCostOfValuation, _2, ptr->dim, valuation())));
+        throw std::logic_error("pfed_t::getCostOfValuation not implemented");
     }
 
     void pfed_t::relax() { for_each(beginMutable(), endMutable(), bind(pdbm_relax, _1, ptr->dim)); }
@@ -371,11 +467,18 @@ namespace dbm
         ptr->zones.push_front(pdbm_t(pdbm, dim));
     }
 
+    void pfed_t::add(pdbm_t pdbm){
+        assert(pdbm.getDimension() == ptr->dim);
+        prepare();
+        ptr->zones.push_front(pdbm);
+    }
+
     void pfed_t::setZero()
     {
         PDBM pdbm = nullptr;
         pdbm_zero(pdbm, ptr->dim);
         *this = pfed_t(pdbm, ptr->dim);
+        pdbm_decRef(pdbm); // We do not reference this, the pdbm_t does now.
     }
 
     void pfed_t::setInit()
@@ -383,6 +486,12 @@ namespace dbm
         PDBM pdbm = nullptr;
         pdbm_init(pdbm, ptr->dim);
         *this = pfed_t(pdbm, ptr->dim);
+        pdbm_decRef(pdbm); // We do not reference this, the pdbm_t does now.
+    }
+
+    void pfed_t::nil()
+    {
+        setDimension(1);
     }
 
     void pfed_t::setEmpty() { *this = pfed_t(ptr->dim); }
@@ -401,12 +510,20 @@ namespace dbm
         return *this;
     }
 
-    pfed_t& pfed_t::operator|=(const pfed_t& fed)
+    fed_t pfed_t::to_fed() const {
+        fed_t fed(getDimension());
+        for (const pdbm_t& pdbm : *this) {
+            fed.add(pdbm.to_dbm());
+        }
+        return fed;
+    }
+
+    pfed_t& pfed_t::operator|=(const pfed_t& pfed)
     {
         // REVISIT: Eliminate included zones.
-        assert(fed.ptr->dim == ptr->dim);
+        assert(pfed.ptr->dim == ptr->dim);
         prepare();
-        ptr->zones.insert(ptr->zones.begin(), fed.begin(), fed.end());
+        ptr->zones.insert(ptr->zones.begin(), pfed.ptr->zones.begin(), pfed.ptr->zones.end());
         return *this;
     }
 
@@ -417,10 +534,14 @@ namespace dbm
         return *this;
     }
 
+    pfed_t operator|(const pfed_t& a, const pfed_t& b) { return pfed_t(a) |= b; }
+
     bool pfed_t::intersects(const pfed_t&) const { throw std::logic_error("pfed_t::intersects not implemented"); }
 
     pfed_t& pfed_t::operator&=(const pfed_t&) { throw std::logic_error("pfed_t::operator &= not implemented"); }
 
+    pfed_t& pfed_t::operator&=(const pdbm_t&) { throw std::logic_error("pfed_t::operator &= not implemented"); }
+
     pfed_t pfed_t::operator&(const pfed_t&) const { throw std::logic_error("pfed_t::operator & not implemented"); }
 
     pfed_t pfed_t::operator-(const pfed_t&) const { throw std::logic_error("pfed_t::operator - not implemented"); }
@@ -429,14 +550,18 @@ namespace dbm
 
     pfed_t& pfed_t::operator+=(const pfed_t&) { throw std::logic_error("pfed_t::operator += not implemented"); }
 
-    void pfed_t::down() { throw std::logic_error("pfed_t::down not implemented"); }
+    pfed_t&  pfed_t::down() { throw std::logic_error("pfed_t::down not implemented"); }
+
+    void pfed_t::resize(const uint32_t* bitSrc, const uint32_t* bitDst, size_t bitSize, cindex_t* table) { throw std::logic_error("pfed_t::resize not implemented"); }
+
+    bool pfed_t::operator==(const pfed_t& arg) const { return relation(arg) == base_EQUAL; }
 
     int32_t pfed_t::getUpperMinimumCost(cindex_t) const
     {
         throw std::logic_error("pfed_t::getUpperMinimumCost not implemented");
     }
 
-    void pfed_t::relaxUp() { throw std::logic_error("pfed_t::relaxUp not implemented"); }
+    pfed_t& pfed_t::relaxUp() { throw std::logic_error("pfed_t::relaxUp not implemented"); }
 
     void pfed_t::getValuation(double*, size_t, bool* freeC) const
     {
@@ -464,8 +589,270 @@ namespace dbm
         pfed_t::const_iterator last = f.end();
         while (first != last) {
             pdbm_print(o, *first, f.getDimension());
-            first++;
+            ++first;
         }
         return o;
     }
+
+    raw_t pfed_t::getMaxUpper(cindex_t clock) const
+    {
+        throw std::logic_error("pfed_t::getMaxUpper not implemented");
+    }
+
+    raw_t pfed_t::getMaxLower(cindex_t clock) const
+    {
+        throw std::logic_error("pfed_t::getMaxLower not implemented");
+    }
+
+    std::string pfed_t::toString(const ClockAccessor& access, bool full) const{
+        bool dbmprinted = false;
+        bool brackets = size() > 1;
+        uint32_t size = getDimension();
+        std::stringstream out;
+        for (auto z = begin(); z != end(); ++z) {
+            bool oneprinted = false;
+            if (dbmprinted) {
+                out << " or ";
+            }
+            if (brackets) {
+                out << "(";
+            }
+            for (uint32_t i = 0; i < size; i++) {
+                for (uint32_t j = 0; j < size; j++) {
+                    raw_t raw = z(i, j);
+                    bool strict = dbm_rawIsStrict(raw);
+                    int32_t bound = dbm_raw2bound(raw);
+                    if (raw < dbm_LS_INFINITY) {
+                        if (i != 0u) {
+                            if (i != j) {
+                                out << (oneprinted ? ", " : "") << access.getClockName(i);
+                                if (j > 0) {
+                                    out << "-" << access.getClockName(j);
+                                }
+                                out << (strict ? "<" : "<=") << bound;
+                                oneprinted = true;
+                            }
+                        } else if ((j != 0u) && ((bound != 0) || strict)) {
+                            // print for i==0,j!=0,bound!=0
+                            out << (oneprinted ? ", " : "") << access.getClockName(j) << (strict ? ">" : ">=")
+                                << -bound;
+                            oneprinted = true;
+                        }
+                    }
+                }
+            }
+            oneprinted = false;
+            out << " :";
+            for (uint32_t i = 1; i < size; ++i) {
+                out << (oneprinted? ",": "") << " r(" << access.getClockName(i) << ") = " << pdbm_getRate(z->operator PDBM(), size, i);
+                oneprinted = true;
+            }
+            if (brackets) {
+                out << ")";
+            }
+            dbmprinted = true;
+        }
+        return out.str();
+    }
+
+    void pfed_t::updateClock(cindex_t x, cindex_t y){
+        throw std::logic_error("pfed_t::updateClock not implemented");
+    }
+
+    pfed_t& pfed_t::upStop(const uint32_t* stopped){
+        throw std::logic_error("pfed_t::upStop not implemented");
+    }
+
+    pfed_t& pfed_t::append(pfed_t& arg){
+        throw std::logic_error("pfed_t::append not implemented");
+    }
+
+    pfed_t& pfed_t::steal(pfed_t& arg){
+        throw std::logic_error("pfed_t::steal not implemented");
+    }
+
+    pfed_t& pfed_t::toLowerBounds() const {
+        throw std::logic_error("pfed_t::toLowerBounds not implemented");
+    }
+
+    pfed_t& pfed_t::toUpperBounds() const {
+        throw std::logic_error("pfed_t::toUpperBounds not implemented");
+    }
+
+    bool pfed_t::isConstrainedBy(cindex_t, cindex_t, raw_t) const {
+        throw std::logic_error("pfed_t::isConstrainedBy not implemented");
+    }
+    bool pfed_t::getDelay(const double* point, cindex_t dim, double* min, double* max, double* minVal, bool* minStrict,
+                          double* maxVal, bool* maxStrict, const uint32_t* stopped) const
+    {
+        throw std::logic_error("pfed_t::getDelay not implemented");
+    }
+
+    bool pdbm_t::constrain(cindex_t i, cindex_t j, raw_t c) {
+        if (!isEmpty()) {
+            return pdbm_constrain1(pdbm, dim, i, j, c);
+        }
+        return false;
+    }
+
+    void pdbm_t::setUniformCost(CostType cost) {
+        if (!isEmpty()) {
+            pdbm_setUniformCost(pdbm, dim, cost);
+        }
+    }
+
+    bool pdbm_t::isEmpty() const {
+        return pdbm_isEmpty(pdbm, dim);
+    }
+
+    bool pfed_t::le(const pfed_t& arg) const {
+        throw std::logic_error("pfed_t::le not implemented");
+    }
+
+    bool pfed_t::lt(const pfed_t& arg) const {
+        throw std::logic_error("pfed_t::lt not implemented");
+    }
+
+    bool pfed_t::gt(const pfed_t& arg) const {
+        throw std::logic_error("pfed_t::gt not implemented");
+    }
+
+    bool pfed_t::ge(const pfed_t& arg) const {
+        throw std::logic_error("pfed_t::ge not implemented");
+    }
+
+    bool pfed_t::eq(const pfed_t& arg) const {
+        throw std::logic_error("pfed_t::eq not implemented");
+    }
+
+    bool pfed_t::unionWith(const pfed_t& other) {
+        assert(ptr->dim == other.ptr->dim);
+        bool changed = false;
+
+        for (auto ot = other.begin(); ot != other.end(); ++ot) {
+            bool maybe_insert = true;
+            bool definitely_insert = false;
+            for (auto it = beginMutable(); it != endMutable();) {
+                switch (pdbm_relation(*ot, *it, ptr->dim)) {
+                case base_SUPERSET:
+                    it.remove();
+                    definitely_insert = true;
+                    break;
+                case base_SUBSET:
+                case base_EQUAL:
+                    maybe_insert = false;
+                case base_DIFFERENT:
+                    ++it;
+                    break;
+                }
+            }
+            if (definitely_insert || maybe_insert) {
+                add(*ot);
+                changed = true;
+            }
+        }
+        return changed;
+    }
+
+
+    int32_t pfed_t::maxOnZero(cindex_t x){
+        throw std::logic_error("pfed_t::maxOnZero not implemented");
+    }
+
+    pfed_t& pfed_t::downStop(const uint32_t* stopped){
+        throw std::logic_error("pfed_t::downStop not implemented");
+    }
+
+    bool pdbm_t::contains(const int32_t* point, cindex_t dim) const{
+        throw std::logic_error("pfed_t::contains not implemented");
+    }
+
+    bool pdbm_t::contains(const double* point, cindex_t dim) const{
+        throw std::logic_error("pfed_t::contains not implemented");
+    }
+
+    std::string pdbm_t::toString(const ClockAccessor&, bool full) const
+    {
+        std::stringstream ss;
+        pdbm_print(ss, *this, dim);
+        return ss.str();
+    }
+
+    bool pfed_t::hasZero() const{
+        throw std::logic_error("pfed_t::hasZero not implemented");
+    }
+
+    void pfed_t::swap(pfed_t& arg){
+        throw std::logic_error("pfed_t::swap not implemented");
+    }
+
+    bool pfed_t::contains(const int32_t* point, cindex_t dim) const{
+        throw std::logic_error("pfed_t::contains not implemented");
+    }
+
+    bool pfed_t::contains(const double * point, cindex_t dim) const{
+        throw std::logic_error("pfed_t::contains not implemented");
+    }
+
+    pfed_t& pfed_t::predt(const pfed_t& bad, const raw_t* restrict){
+        throw std::logic_error("pfed_t::predt not implemented");
+    }
+
+    void pfed_t::intern(){
+        throw std::logic_error("pfed_t::intern not implemented");
+    }
+
+    pfed_t& pfed_t::mergeReduce(size_t skip, int expensiveTry){
+        throw std::logic_error("pfed_t::mergeReduce not implemented");
+    }
+    void pfed_t::simplify_rational_cost() {
+        for (auto it = beginMutable(); it != endMutable(); ++it) {
+            it->simplify_rational_cost();
+        }
+    }
+
+    cindex_t pdbm_t::pdim() const { return dim; }
+
+
+#define DBM_IJ(DBM, I, J) DBM[(I)*dim + (J)]
+#define DBM(I, J)         DBM_IJ(dbm, I, J)
+
+    bool pdbm_t::canDelay() const {
+        if (isEmpty())
+            return false;
+        const raw_t* dbm = const_dbm();
+        cindex_t dim = pdim();
+        for (cindex_t i = 1; i < dim; ++i) {
+            if (dbm_rawIsWeak(DBM(i, 0)) && dbm_rawIsWeak(DBM(0, i)) && DBM(i, 0) == dbm_weakNegRaw(DBM(0, i))) {
+                return false;
+            }
+        }
+        return true;
+    }
+    relation_t pdbm_t::relation(const pdbm_t& other) const {
+        assert(dim == other.dim);
+        return pdbm_relation(pdbm, other.pdbm, dim);
+    }
+
+    relation_t pdbm_t::strict_relation(const pdbm_t& other) const {
+        assert(dim == other.dim);
+        return pdbm_relation_strict(pdbm, other.pdbm, dim);
+    }
+
+    pdbm_t& pdbm_t::operator&=(const dbm::pdbm_t& arg) {
+        assert(getDimension() == arg.getDimension());
+        if (arg.isEmpty()) {
+            pdbm_decRef(pdbm);
+            pdbm = nullptr;
+        } else if (!isEmpty()) {
+            pdbm_intersect(pdbm, arg.pdbm, dim);
+        }
+        return *this;
+    }
+    void pdbm_t::simplify_rational_cost() {
+        if (!isEmpty()) {
+            pdbm_simplify_rational_cost(pdbm, dim);
+        }
+    }
 }  // namespace dbm
+
diff --git a/src/priced.cpp b/src/priced.cpp
index 0fabd9d..443f4b2 100644
--- a/src/priced.cpp
+++ b/src/priced.cpp
@@ -21,12 +21,14 @@
 #include <stdexcept>
 #include <cassert>
 #include <cstdlib>
+#include <limits>
 
 struct PDBM_s
 {
     uint32_t count;
-    uint32_t cost;
-    uint32_t infimum;
+    CostType cost;
+    CostType infimum;
+    CostType* rates;
     int32_t data[];
 };
 
@@ -37,10 +39,10 @@ struct PDBM_s
 #define pdbm_cost(pdbm) ((pdbm)->cost)
 
 /** Returns the vectors of coefficients. */
-#define pdbm_rates(pdbm) ((pdbm)->data)
+#define pdbm_rates(pdbm) ((pdbm)->rates)
 
 /** Returns the matrix. */
-#define pdbm_matrix(pdbm) ((pdbm)->data + dim)
+#define pdbm_matrix(pdbm) ((pdbm)->data)
 
 /** Returns the cache infimum. */
 #define pdbm_cache(pdbm) ((pdbm)->infimum)
@@ -49,7 +51,7 @@ struct PDBM_s
 #define pdbm_count(pdbm) ((pdbm)->count)
 
 /** Constant to mark the cached infimum as void. */
-#define INVALID INT_MAX
+#define INVALID INFINITE_COST
 
 #ifndef NDEBUG
 /**
@@ -132,7 +134,8 @@ static void pdbm_blank(PDBM& pdbm, cindex_t dim)
 size_t pdbm_size(cindex_t dim)
 {
     assert(dim);
-    return sizeof(struct PDBM_s) + (dim * dim + dim) * sizeof(int32_t);
+    return sizeof(struct PDBM_s) + (dim * dim) * sizeof(int32_t);
+    //return sizeof(struct PDBM_s) + (dim * dim) * sizeof(int32_t) + dim * sizeof(CostType);
 }
 
 PDBM pdbm_reserve(cindex_t dim, void* p)
@@ -142,7 +145,10 @@ PDBM pdbm_reserve(cindex_t dim, void* p)
     PDBM pdbm = (PDBM)p;
     pdbm_count(pdbm) = 0;
     pdbm_cache(pdbm) = INVALID;
-    pdbm_rates(pdbm)[0] = 0;
+    pdbm_rates(pdbm) = (CostType*)malloc(sizeof(CostType) * dim);
+    pdbm_cost(pdbm) = 0;
+    std::fill(pdbm_rates(pdbm), pdbm_rates(pdbm) + dim, 0);
+
     return pdbm;
 }
 
@@ -156,6 +162,7 @@ void pdbm_deallocate(PDBM& pdbm)
 {
     assert(pdbm == nullptr || pdbm_count(pdbm) == 0);
 
+    free(pdbm->rates);
     free(pdbm);
 
     /* Setting the pointer to NULL protects against accidental use of
@@ -172,7 +179,11 @@ PDBM pdbm_copy(PDBM dst, const PDBM src, cindex_t dim)
         dst = pdbm_allocate(dim);
     }
 
+    free(dst->rates); // We overwrite the rates pointer in the memcpy below, so make sure to free it first.
+
     memcpy(dst, src, pdbm_size(dim));
+    dst->rates = (CostType*)malloc(sizeof(CostType) * dim);
+    std::copy(src->rates, src->rates + dim, dst->rates);
     pdbm_count(dst) = 0;
 
     return dst;
@@ -186,7 +197,7 @@ void pdbm_init(PDBM& pdbm, cindex_t dim)
     dbm_init(pdbm_matrix(pdbm), dim);
     pdbm_cost(pdbm) = 0;
     pdbm_cache(pdbm) = 0;
-    auto* r = pdbm_rates(pdbm);
+    CostType* r = pdbm_rates(pdbm);
     std::fill(r, r + dim, 0);
 
     assertx(pdbm_isValid(pdbm, dim));
@@ -200,7 +211,7 @@ void pdbm_zero(PDBM& pdbm, cindex_t dim)
     dbm_zero(pdbm_matrix(pdbm), dim);
     pdbm_cost(pdbm) = 0;
     pdbm_cache(pdbm) = 0;
-    auto* r = pdbm_rates(pdbm);
+    CostType* r = pdbm_rates(pdbm);
     std::fill(r, r + dim, 0);
 
     assertx(pdbm_isValid(pdbm, dim));
@@ -238,8 +249,8 @@ bool pdbm_constrain1(PDBM& pdbm, cindex_t dim, cindex_t i, cindex_t j, raw_t con
 
     /* Compute the cost at the origin.
      */
-    uint32_t cost = pdbm_cost(pdbm);
-    int32_t* rates = pdbm_rates(pdbm);
+    CostType cost = pdbm_cost(pdbm);
+    CostType* rates = pdbm_rates(pdbm);
     for (uint32_t k = 1; k < dim; k++) {
         cost += rates[k] * dbm_raw2bound(DBM(0, k));
     }
@@ -273,8 +284,8 @@ bool pdbm_constrainN(PDBM& pdbm, cindex_t dim, const constraint_t* constraints,
             pdbm_prepare(pdbm, dim);
 
             dbm = pdbm_matrix(pdbm);
-            uint32_t cost = pdbm_cost(pdbm);
-            int32_t* rates = pdbm_rates(pdbm);
+            CostType cost = pdbm_cost(pdbm);
+            CostType* rates = pdbm_rates(pdbm);
 
             for (uint32_t k = 1; k < dim; k++) {
                 cost += rates[k] * dbm_raw2bound(DBM(0, k));
@@ -306,7 +317,7 @@ bool pdbm_constrainN(PDBM& pdbm, cindex_t dim, const constraint_t* constraints,
  *               we want to find the cost.
  * @param dim    is the dimension of \a dbm1 and \a dbm2.
  */
-static int32_t costAtOtherOffset(const raw_t* dbm1, const int32_t* rates1, uint32_t cost1, const raw_t* dbm2,
+static CostType costAtOtherOffset(const raw_t* dbm1, const CostType* rates1, CostType cost1, const raw_t* dbm2,
                                  cindex_t dim)
 {
     assert(dbm1 && dbm2 && rates1 && dim);
@@ -328,11 +339,56 @@ static int32_t costAtOtherOffset(const raw_t* dbm1, const int32_t* rates1, uint3
  * @param last  is a pointer to the last coefficient of the first plane.
  * @param rate  is a pointer to the first coefficient of the second plane.
  */
-inline static bool leq(const int32_t* first, const int32_t* last, const int32_t* rate)
+inline static bool leq(const CostType* first, const CostType* last, const CostType* rate)
+{
+    assert(first && last && rate);
+
+    return std::equal(first, last, rate, std::less_equal<>());
+}
+
+/**
+ * Given two planes, returns true if the slope of the first is strictly less
+ * than to the slope of the other.
+ *
+ * @param first is a pointer to the first coefficient of the first plane.
+ * @param last  is a pointer to the last coefficient of the first plane.
+ * @param rate  is a pointer to the first coefficient of the second plane.
+ */
+inline static bool le(const CostType* first, const CostType* last, const CostType* rate)
 {
     assert(first && last && rate);
 
-    return std::equal(first, last, rate, std::less_equal<int32_t>());
+    return std::equal(first, last, rate, std::less<>());
+}
+
+/**
+ * Given two planes, returns true if the slope of the first is equal to
+ * the slope of the other.
+ *
+ * @param first is a pointer to the first coefficient of the first plane.
+ * @param last  is a pointer to the last coefficient of the first plane.
+ * @param rate  is a pointer to the first coefficient of the second plane.
+ */
+inline static bool eq(const CostType* first, const CostType* last, const CostType* rate)
+{
+    assert(first && last && rate);
+
+    return std::equal(first, last, rate, std::equal_to<>());
+}
+
+/**
+ * Given two planes, returns true if the slope of the first is less
+ * than to the slope of the other.
+ *
+ * @param first is a pointer to the first coefficient of the first plane.
+ * @param last  is a pointer to the last coefficient of the first plane.
+ * @param rate  is a pointer to the first coefficient of the second plane.
+ */
+inline static bool less(const CostType* first, const CostType* last, const CostType* rate)
+{
+    assert(first && last && rate);
+
+    return std::equal(first, last, rate, std::less<>());
 }
 
 /**
@@ -346,30 +402,78 @@ inline static bool leq(const int32_t* first, const int32_t* last, const int32_t*
  * @param cost2 The cost at the offset for the second plane.
  * @param rate2 The rates of the second plane.
  */
-static int32_t infOfDiff(const raw_t* dbm, uint32_t dim, int32_t cost1, const int32_t* rate1, int32_t cost2,
-                         const int32_t* rate2)
+static CostType infOfDiff(const raw_t* dbm, uint32_t dim, CostType cost1, const CostType* rate1, CostType cost2,
+                        const CostType* rate2)
 {
     assert(dbm && dim && rate1 && rate2);
 
-    int32_t cost = cost1 - cost2;
-    int32_t rates[dim];
-    std::transform(rate1, rate1 + dim, rate2, rates, std::minus<int32_t>());
+    CostType cost = cost1 - cost2;
+    CostType rates[dim];
+    std::transform(rate1, rate1 + dim, rate2, rates, std::minus<>());
 
     return pdbm_infimum(dbm, dim, cost, rates);
 }
 
+std::pair<relation_t, bool> pdbm_compare_cost_identical_pdbms(const PDBM pdbm1, const PDBM pdbm2, cindex_t dim)
+{
+    assert(pdbm1 && pdbm2 && dim);
+
+    raw_t* dbm1 = pdbm_matrix(pdbm1);
+    raw_t* dbm2 = pdbm_matrix(pdbm2);
+    CostType* rates1 = pdbm_rates(pdbm1);
+    CostType* rates2 = pdbm_rates(pdbm2);
+    CostType cost1 = pdbm_cost(pdbm1);
+    CostType cost2 = pdbm_cost(pdbm2);
+
+    /* Function should only be called for identical dbms.
+     * In case they are not, return costs are different. */
+    if (dbm_relation(dbm1, dbm2, dim) != base_EQUAL) {
+        return {base_DIFFERENT, false};
+    }
+
+    /* Both have the same size. We need to compare the planes to
+         * see which one is cheaper.
+     */
+    CostType c = infOfDiff(dbm1, dim, cost2, rates2, cost1, rates1);
+    if (c > 0) {
+        return {base_SUPERSET, true};
+    }
+
+    CostType d = infOfDiff(dbm1, dim, cost1, rates1, cost2, rates2);
+    if (c == 0 && d == 0) {
+        return {base_EQUAL, false};
+    }
+    if (d > 0) {
+        return {base_SUBSET, true};
+    }
+    if (c >= 0) {
+        return {base_SUPERSET, false};
+    }
+    if (d >= 0) {
+        return {base_SUBSET, false};
+    }
+    return {base_DIFFERENT, false};
+}
+
+/*
+ * Returns
+ * (Z,c) ⊃ (Z',c') iff (Z ⊃ Z' ∧ c ≤ c') ∨ (Z = Z' ∧ c ≤ c')
+ * (Z,c) ⊂ (Z',c') iff (Z ⊂ Z' ∧ c ≥ c') ∨ (Z = Z' ∧ c ≥ c')
+ * (Z,c) = (Z',c') iff Z = Z' ∧ c = c'
+ */
 relation_t pdbm_relation(const PDBM pdbm1, const PDBM pdbm2, cindex_t dim)
 {
     assert(pdbm1 && pdbm2 && dim);
 
     raw_t* dbm1 = pdbm_matrix(pdbm1);
     raw_t* dbm2 = pdbm_matrix(pdbm2);
-    int32_t* rates1 = pdbm_rates(pdbm1);
-    int32_t* rates2 = pdbm_rates(pdbm2);
-    uint32_t cost1 = pdbm_cost(pdbm1);
-    uint32_t cost2 = pdbm_cost(pdbm2);
+    CostType* rates1 = pdbm_rates(pdbm1);
+    CostType* rates2 = pdbm_rates(pdbm2);
+    CostType cost1 = pdbm_cost(pdbm1);
+    CostType cost2 = pdbm_cost(pdbm2);
 
-    int32_t c, d;
+    bool a, b;
+    CostType c, d;
 
     switch (dbm_relation(dbm1, dbm2, dim)) {
     case base_SUPERSET:
@@ -407,14 +511,14 @@ relation_t pdbm_relation(const PDBM pdbm1, const PDBM pdbm2, cindex_t dim)
 
         /* Do sound and cheap comparison first.
          */
-        c = cost1 <= cost2 && leq(rates1, rates1 + dim, rates2);
-        d = cost2 <= cost1 && leq(rates2, rates2 + dim, rates1);
+        a = cost1 <= cost2 && leq(rates1, rates1 + dim, rates2);
+        b = cost2 <= cost1 && leq(rates2, rates2 + dim, rates1);
 
-        if (c & d) {
+        if (a && b) {
             return base_EQUAL;
-        } else if (c) {
+        } else if (a) {
             return base_SUPERSET;
-        } else if (d) {
+        } else if (b) {
             return base_SUBSET;
         }
 
@@ -451,57 +555,39 @@ relation_t pdbm_relation(const PDBM pdbm1, const PDBM pdbm2, cindex_t dim)
     }
 }
 
-relation_t pdbm_relationWithMinDBM(const PDBM pdbm1, cindex_t dim, const mingraph_t pdbm2, raw_t* dbm2)
+/*
+ * Returns
+ * (Z,c) ⊃ (Z',c') iff (Z ⊇ Z' ∧ c < c')                        // (Z,c) strictly dominates (Z',c')
+ * (Z,c) ⊂ (Z',c') iff (Z ⊆ Z' ∧ c > c') ∨ (Z ⊂ Z' ∧ c ≥ c')    // (Z',c') weakly dominates (Z,c)
+ * (Z,c) = (Z',c') iff Z = Z' ∧ c = c'                          // (Z,c) = (Z',c')
+ */
+relation_t pdbm_relation_strict(const PDBM pdbm1, const PDBM pdbm2, cindex_t dim)
 {
-    assert(pdbm1 && pdbm2 && dim && dbm2);
+    assert(pdbm1 && pdbm2 && dim);
 
     raw_t* dbm1 = pdbm_matrix(pdbm1);
-    uint32_t cost1 = pdbm_cost(pdbm1);
-    int32_t* rates1 = pdbm_rates(pdbm1);
-
-    int32_t c, d;
+    raw_t* dbm2 = pdbm_matrix(pdbm2);
+    CostType* rates1 = pdbm_rates(pdbm1);
+    CostType* rates2 = pdbm_rates(pdbm2);
+    CostType cost1 = pdbm_cost(pdbm1);
+    CostType cost2 = pdbm_cost(pdbm2);
 
-    /* We know how the cost and the rates are encoded in a mingraph:
-     * see writeToMinDBMWithOffset and readFromMinDBM.
-     */
-    uint32_t cost2 = pdbm2[0];
-    const int32_t* rates2 = pdbm2 + 2;
-
-    /* dbm_relationWithMinDBM will in some cases unpack pdbm2 into
-     * dbm2. In order to detect whether this has happened, we set the
-     * first entry of dbm2 to -1. If dbm_relationWithMinDBM did indeed
-     * unpack pdbm2, this entry will have a different value
-     * afterwards.
-     */
-    dbm2[0] = -1;
+    bool a, b;
+    CostType c, d;
 
-    switch (dbm_relationWithMinDBM(dbm1, dim, pdbm2 + dim + 2, dbm2)) {
+    switch (dbm_relation(dbm1, dbm2, dim)) {
     case base_SUPERSET:
-        /* pdbm2 is smaller. Check whether it is also the most
-         * expensive: This is the case when the difference between
-         * pdbm2 and pdbm1 is always non-negative (the infimum is not
-         * smaller than 0).
-         */
-        if (dbm2[0] == -1) {
-            dbm_readFromMinDBM(dbm2, pdbm2 + dim + 2);
-        }
+        // pdbm2 is smaller. Check whether it is also strictly more expensive:
         cost1 = costAtOtherOffset(dbm1, rates1, cost1, dbm2, dim);
-        if (cost1 <= cost2 &&
-            (leq(rates1, rates1 + dim, rates2) || infOfDiff(dbm2, dim, cost2, rates2, cost1, rates1) >= 0)) {
+        if (cost1 < cost2 &&
+            (leq(rates1, rates1 + dim, rates2) || infOfDiff(dbm2, dim, cost2, rates2, cost1, rates1) > 0)) {
             return base_SUPERSET;
         } else {
             return base_DIFFERENT;
         }
 
     case base_SUBSET:
-        /* pdbm2 is bigger. Check whether it is also the cheapest:
-         * This is the case when the difference between pdbm1 and
-         * pdbm2 is always non-negative (the infimum is not smaller
-         * than 0).
-         */
-        if (dbm2[0] == -1) {
-            dbm_readFromMinDBM(dbm2, pdbm2 + dim + 2);
-        }
+        // pdbm2 is bigger. Check whether it is also weakly less expensive:
         cost2 = costAtOtherOffset(dbm2, rates2, cost2, dbm1, dim);
         if (cost2 <= cost1 &&
             (leq(rates2, rates2 + dim, rates1) || infOfDiff(dbm1, dim, cost1, rates1, cost2, rates2) >= 0)) {
@@ -517,14 +603,11 @@ relation_t pdbm_relationWithMinDBM(const PDBM pdbm1, cindex_t dim, const mingrap
 
         /* Do sound and cheap comparison first.
          */
-        c = cost1 <= cost2 && leq(rates1, rates1 + dim, rates2);
-        d = cost2 <= cost1 && leq(rates2, rates2 + dim, rates1);
-
-        if (c & d) {
+        if (cost1 == cost2 && eq(rates1, rates1 + dim, rates2)) {
             return base_EQUAL;
-        } else if (c) {
+        } else if (cost1 < cost2 && leq(rates1, rates1 + dim, rates2)) {
             return base_SUPERSET;
-        } else if (d) {
+        } else if (cost2 <= cost1 && leq(rates2, rates2 + dim, rates1)) {
             return base_SUBSET;
         }
 
@@ -535,7 +618,7 @@ relation_t pdbm_relationWithMinDBM(const PDBM pdbm1, cindex_t dim, const mingrap
          * situation where both infima are zero.
          *
          * Notice that dbm1 equals dbm2, hence we do not need to
-         * unpacked dbm2.
+         * unpack dbm2.
          */
         c = infOfDiff(dbm1, dim, cost2, rates2, cost1, rates1);
         if (c > 0) {
@@ -549,9 +632,6 @@ relation_t pdbm_relationWithMinDBM(const PDBM pdbm1, cindex_t dim, const mingrap
         if (c == 0 && d == 0) {
             return base_EQUAL;
         }
-        if (c >= 0) {
-            return base_SUPERSET;
-        }
         if (d >= 0) {
             return base_SUBSET;
         }
@@ -561,29 +641,162 @@ relation_t pdbm_relationWithMinDBM(const PDBM pdbm1, cindex_t dim, const mingrap
     }
 }
 
-int32_t pdbm_getInfimum(const PDBM pdbm, cindex_t dim)
+relation_t pdbm_relationWithMinDBM(const PDBM pdbm1, cindex_t dim, const mingraph_t pdbm2, raw_t* dbm2)
+{
+    throw std::logic_error("pdbm_relationWithMinDBM not implemented!");
+
+//    assert(pdbm1 && pdbm2 && dim && dbm2);
+//
+//    raw_t* dbm1 = pdbm_matrix(pdbm1);
+//    CostType cost1 = pdbm_cost(pdbm1);
+//    CostType* rates1 = pdbm_rates(pdbm1);
+//
+//    CostType c, d;
+//
+//    /* We know how the cost and the rates are encoded in a mingraph:
+//     * see writeToMinDBMWithOffset and readFromMinDBM.
+//     */
+//    CostType cost2 = pdbm2[0];
+//    const CostType* rates2 = pdbm2 + 2;
+//
+//    /* dbm_relationWithMinDBM will in some cases unpack pdbm2 into
+//     * dbm2. In order to detect whether this has happened, we set the
+//     * first entry of dbm2 to -1. If dbm_relationWithMinDBM did indeed
+//     * unpack pdbm2, this entry will have a different value
+//     * afterwards.
+//     */
+//    dbm2[0] = -1;
+//
+//    switch (dbm_relationWithMinDBM(dbm1, dim, pdbm2 + dim + 2, dbm2)) {
+//    case base_SUPERSET:
+//        /* pdbm2 is smaller. Check whether it is also the most
+//         * expensive: This is the case when the difference between
+//         * pdbm2 and pdbm1 is always non-negative (the infimum is not
+//         * smaller than 0).
+//         */
+//        if (dbm2[0] == -1) {
+//            dbm_readFromMinDBM(dbm2, pdbm2 + dim + 2);
+//        }
+//        cost1 = costAtOtherOffset(dbm1, rates1, cost1, dbm2, dim);
+//        if (cost1 <= cost2 &&
+//            (leq(rates1, rates1 + dim, rates2) || infOfDiff(dbm2, dim, cost2, rates2, cost1, rates1) >= 0)) {
+//            return base_SUPERSET;
+//        } else {
+//            return base_DIFFERENT;
+//        }
+//
+//    case base_SUBSET:
+//        /* pdbm2 is bigger. Check whether it is also the cheapest:
+//         * This is the case when the difference between pdbm1 and
+//         * pdbm2 is always non-negative (the infimum is not smaller
+//         * than 0).
+//         */
+//        if (dbm2[0] == -1) {
+//            dbm_readFromMinDBM(dbm2, pdbm2 + dim + 2);
+//        }
+//        cost2 = costAtOtherOffset(dbm2, rates2, cost2, dbm1, dim);
+//        if (cost2 <= cost1 &&
+//            (leq(rates2, rates2 + dim, rates1) || infOfDiff(dbm1, dim, cost1, rates1, cost2, rates2) >= 0)) {
+//            return base_SUBSET;
+//        } else {
+//            return base_DIFFERENT;
+//        }
+//
+//    case base_EQUAL:
+//        /* Both have the same size. We need to compare the planes to
+//         * see which one is cheaper.
+//         */
+//
+//        /* Do sound and cheap comparison first.
+//         */
+//        c = cost1 <= cost2 && leq(rates1, rates1 + dim, rates2);
+//        d = cost2 <= cost1 && leq(rates2, rates2 + dim, rates1);
+//
+//        if (c & d) {
+//            return base_EQUAL;
+//        } else if (c) {
+//            return base_SUPERSET;
+//        } else if (d) {
+//            return base_SUBSET;
+//        }
+//
+//        /* The planes are incomparable, so we need to do the
+//         * subtraction. Notice that priced zones are not canonical,
+//         * so the two zones can in fact be equal even though the rates
+//         * are different.  Therefore we must also check for the
+//         * situation where both infima are zero.
+//         *
+//         * Notice that dbm1 equals dbm2, hence we do not need to
+//         * unpacked dbm2.
+//         */
+//        c = infOfDiff(dbm1, dim, cost2, rates2, cost1, rates1);
+//        if (c > 0) {
+//            /* Early return to avoid unnecessary computation of the
+//             * second subtraction.
+//             */
+//            return base_SUPERSET;
+//        }
+//
+//        d = infOfDiff(dbm1, dim, cost1, rates1, cost2, rates2);
+//        if (c == 0 && d == 0) {
+//            return base_EQUAL;
+//        }
+//        if (c >= 0) {
+//            return base_SUPERSET;
+//        }
+//        if (d >= 0) {
+//            return base_SUBSET;
+//        }
+//        return base_DIFFERENT;
+//
+//    default: return base_DIFFERENT;
+//    }
+}
+
+CostType pdbm_getInfimum(const PDBM pdbm, cindex_t dim)
 {
     assert(pdbm && dim);
     assert(dbm_isValid(pdbm_matrix(pdbm), dim));
-    uint32_t cache = pdbm_cache(pdbm);
+    CostType cache = pdbm_cache(pdbm);
     if (cache == INVALID) {
         pdbm_cache((PDBM)pdbm) = cache = pdbm_infimum(pdbm_matrix(pdbm), dim, pdbm_cost(pdbm), pdbm_rates(pdbm));
     }
     return cache;
 }
 
-int32_t pdbm_getInfimumValuation(const PDBM pdbm, cindex_t dim, int32_t* valuation, const bool* free)
+CostType pdbm_getSupremum(const PDBM pdbm, cindex_t dim)
+{
+    assert(pdbm && dim);
+    assert(dbm_isValid(pdbm_matrix(pdbm), dim));
+
+    // Supremum is infimum of zone where rates are inversed
+    CostType inv_rates[dim];
+    std::transform(pdbm_rates(pdbm), pdbm_rates(pdbm) + dim, inv_rates, [](CostType c){ return -c; });
+    return -pdbm_infimum(pdbm_matrix(pdbm), dim, -pdbm_cost(pdbm), inv_rates);
+}
+
+void pdbm_setUniformCost(PDBM& pdbm, cindex_t dim, CostType cost) {
+    pdbm_prepare(pdbm, dim);
+    pdbm_setCostAtOffset(pdbm, dim, cost);
+    for (cindex_t x = 1; x < dim; ++x) {
+        pdbm_setRate(pdbm, dim, x, 0);
+    }
+    pdbm_cache(pdbm) = cost;
+}
+
+
+CostType pdbm_getInfimumValuation(const PDBM pdbm, cindex_t dim, int32_t* valuation, const bool* free)
 {
     assert(pdbm && dim && valuation);
     assert(pdbm_isValid(pdbm, dim));
 
     raw_t copy[dim * dim];
     raw_t* dbm = pdbm_matrix(pdbm);
-    int32_t* rates = pdbm_rates(pdbm);
+    CostType* rates = pdbm_rates(pdbm);
 
     /* Compute the cost of the origin.
      */
-    int32_t cost = pdbm_cost(pdbm);
+    CostType cost = pdbm_cost(pdbm);
     for (uint32_t i = 1; i < dim; i++) {
         cost -= rates[i] * -dbm_raw2bound(DBM(0, i));
     }
@@ -623,7 +836,7 @@ bool pdbm_satisfies(const PDBM pdbm, cindex_t dim, cindex_t i, cindex_t j, raw_t
 
 bool pdbm_isEmpty(const PDBM pdbm, cindex_t dim)
 {
-    assert(pdbm && dim);
+    assert(dim);
     return pdbm == nullptr || dbm_isEmpty(pdbm_matrix(pdbm), dim);
 }
 
@@ -636,7 +849,12 @@ bool pdbm_isUnbounded(const PDBM pdbm, cindex_t dim)
 uint32_t pdbm_hash(const PDBM pdbm, cindex_t dim, uint32_t seed)
 {
     assert(pdbm && dim && !(pdbm_size(dim) & 3));
-    return hash_computeI32((int32_t*)pdbm, pdbm_size(dim) >> 2, seed);
+    uint32_t hash = seed ^ dbm_hash(pdbm_matrix(pdbm), dim);
+    hash ^= hash_cost_type(pdbm_cost(pdbm));
+    for (cindex_t i = 1; i < dim; ++i) {
+        hash ^= (hash ^ hash_cost_type(pdbm_rates(pdbm)[i])) << 4;
+    }
+    return hash;
 }
 
 static bool isPointIncludedWeakly(const int32_t* pt, const raw_t* dbm, cindex_t dim)
@@ -682,11 +900,14 @@ void pdbm_up(PDBM& pdbm, cindex_t dim)
 
     pdbm_prepare(pdbm, dim);
     dbm_up(pdbm_matrix(pdbm), dim);
+    if (pdbm_getSlopeOfDelayTrajectory(pdbm, dim) < 0) {
+        pdbm->infimum = INVALID;
+    }
 
     assertx(pdbm_isValid(pdbm, dim));
 }
 
-void pdbm_upZero(PDBM& pdbm, cindex_t dim, int32_t rate, cindex_t zero)
+void pdbm_upZero(PDBM& pdbm, cindex_t dim, CostType rate, cindex_t zero)
 {
     assert(pdbm && dim && zero > 0 && zero < dim);
     assert(pdbm_areOnZeroCycle(pdbm, dim, 0, zero));
@@ -694,11 +915,12 @@ void pdbm_upZero(PDBM& pdbm, cindex_t dim, int32_t rate, cindex_t zero)
     pdbm_prepare(pdbm, dim);
 
     raw_t* dbm = pdbm_matrix(pdbm);
-    int32_t* rates = pdbm_rates(pdbm);
+    CostType* rates = pdbm_rates(pdbm);
 
     dbm_up(dbm, dim);
     rates[zero] = 0;
     rates[zero] = rate - pdbm_getSlopeOfDelayTrajectory(pdbm, dim);
+    pdbm_cache(pdbm) = INVALID;
 
     assertx(pdbm_isValid(pdbm, dim));
 }
@@ -722,12 +944,12 @@ void pdbm_updateValueZero(PDBM& pdbm, cindex_t dim, cindex_t clock, uint32_t val
 
     pdbm_prepare(pdbm, dim);
 
-    int32_t* rates = pdbm_rates(pdbm);
+    CostType* rates = pdbm_rates(pdbm);
 
     if (zero) {
         rates[zero] += rates[clock];
     }
-    rates[clock] = 0;
+//    rates[clock] = 0;
     dbm_updateValue(pdbm_matrix(pdbm), dim, clock, value);
 
     assertx(pdbm_isValid(pdbm, dim));
@@ -781,7 +1003,7 @@ void pdbm_extrapolateMaxBounds(PDBM& pdbm, cindex_t dim, int32_t* max)
 {
     assert(pdbm && dim);
 
-    int32_t* rates = pdbm_rates(pdbm);
+    CostType* rates = pdbm_rates(pdbm);
 
     for (uint32_t i = 1; i < dim; i++) {
         if (rates[i] != 0) {
@@ -804,7 +1026,7 @@ void pdbm_diagonalExtrapolateMaxBounds(PDBM& pdbm, cindex_t dim, int32_t* max)
 {
     assert(pdbm && dim);
 
-    int32_t* rates = pdbm_rates(pdbm);
+    CostType* rates = pdbm_rates(pdbm);
 
     /* If possible, transfer the cost rate of inactive clocks to other
      * clocks (clocks which are on a zero cycle with the inactive
@@ -834,11 +1056,11 @@ void pdbm_diagonalExtrapolateMaxBounds(PDBM& pdbm, cindex_t dim, int32_t* max)
     }
 
     /* It is only safe to extrapolate clocks with a zero cost rate. */
-    for (cindex_t i = 1; i < dim; i++) {
-        if (rates[i] != 0) {
-            max[i] = dbm_INFINITY;
-        }
-    }
+//    for (cindex_t i = 1; i < dim; i++) {
+//        if (rates[i] != 0) {
+//            max[i] = dbm_INFINITY;
+//        }
+//    }
 
     pdbm_prepare(pdbm, dim);
     dbm_diagonalExtrapolateMaxBounds(pdbm_matrix(pdbm), dim, max);
@@ -848,7 +1070,7 @@ void pdbm_diagonalExtrapolateLUBounds(PDBM& pdbm, cindex_t dim, int32_t* lower,
 {
     assert(pdbm && dim);
 
-    int32_t* rates = pdbm_rates(pdbm);
+    CostType* rates = pdbm_rates(pdbm);
 
     for (uint32_t i = 1; i < dim; i++) {
         //          if (rates[i] < 0)
@@ -869,13 +1091,15 @@ void pdbm_diagonalExtrapolateLUBounds(PDBM& pdbm, cindex_t dim, int32_t* lower,
     dbm_diagonalExtrapolateLUBounds(pdbm_matrix(pdbm), dim, lower, upper);
 }
 
-void pdbm_incrementCost(PDBM& pdbm, cindex_t dim, int32_t value)
+void pdbm_incrementCost(PDBM& pdbm, cindex_t dim, CostType value)
 {
-    assert(pdbm && dim && value >= 0);
+    assert(pdbm && dim);
 
     pdbm_prepare(pdbm, dim);
     pdbm_cost(pdbm) += value;
-    pdbm_cache(pdbm) += value;
+    if (pdbm_cache(pdbm) != INVALID && pdbm_cache(pdbm) != -INFINITE_COST) {
+        pdbm_cache(pdbm) += value;
+    }
 
     assertx(pdbm_isValid(pdbm, dim));
 }
@@ -902,28 +1126,30 @@ size_t pdbm_analyzeForMinDBM(const PDBM pdbm, cindex_t dim, uint32_t* bitMatrix)
 int32_t* pdbm_writeToMinDBMWithOffset(const PDBM pdbm, cindex_t dim, bool minimizeGraph, bool tryConstraints16,
                                       allocator_t allocator, uint32_t offset)
 {
-    assert(pdbm && dim);
-
-    int32_t* graph = dbm_writeToMinDBMWithOffset(pdbm_matrix(pdbm), dim, minimizeGraph, tryConstraints16, allocator,
-                                                 offset + dim + 2);
-    graph[offset] = pdbm_cost(pdbm);
-    graph[offset + 1] = pdbm_cache(pdbm);
-    auto* r = pdbm_rates(pdbm);
-    std::copy(r, r + dim, graph + offset + 2);
-    return graph;
+    throw std::logic_error("pdbm_writeToMinDBMWithOffset not implemented");
+//    assert(pdbm && dim);
+//
+//    int32_t* graph = dbm_writeToMinDBMWithOffset(pdbm_matrix(pdbm), dim, minimizeGraph, tryConstraints16, allocator,
+//                                                 offset + dim + 2);
+//    graph[offset] = pdbm_cost(pdbm);
+//    graph[offset + 1] = pdbm_cache(pdbm);
+//    auto* r = pdbm_rates(pdbm);
+//    std::copy(r, r + dim, graph + offset + 2);
+//    return graph;
 }
 
 void pdbm_readFromMinDBM(PDBM& dst, cindex_t dim, mingraph_t src)
 {
-    assert(dst && dim);
-
-    pdbm_blank(dst, dim);
-    pdbm_cost(dst) = src[0];
-    pdbm_cache(dst) = src[1];
-    std::copy(src + 2, src + 2 + dim, pdbm_rates(dst));
-    dbm_readFromMinDBM(pdbm_matrix(dst), src + dim + 2);
-
-    assertx(pdbm_isValid(dst, dim));
+    throw std::logic_error("pdbm_readFromMinDBM not implemented");
+//    assert(dst && dim);
+//
+//    pdbm_blank(dst, dim);
+//    pdbm_cost(dst) = src[0];
+//    pdbm_cache(dst) = src[1];
+//    std::copy(src + 2, src + 2 + dim, pdbm_rates(dst));
+//    dbm_readFromMinDBM(pdbm_matrix(dst), src + dim + 2);
+//
+//    assertx(pdbm_isValid(dst, dim));
 }
 
 bool pdbm_findNextZeroCycle(const PDBM pdbm, cindex_t dim, cindex_t x, cindex_t* out)
@@ -951,33 +1177,33 @@ bool pdbm_findZeroCycle(const PDBM pdbm, cindex_t dim, cindex_t x, cindex_t* out
     return pdbm_findNextZeroCycle(pdbm, dim, x, out);
 }
 
-int32_t pdbm_getSlopeOfDelayTrajectory(const PDBM pdbm, cindex_t dim)
+CostType pdbm_getSlopeOfDelayTrajectory(const PDBM pdbm, cindex_t dim)
 {
     assert(pdbm && dim);
 
-    int32_t* rates = pdbm_rates(pdbm);
-    int32_t sum = 0;
+    CostType* rates = pdbm_rates(pdbm);
+    CostType sum = 0;
     for (uint32_t i = 1; i < dim; i++) {
         sum += rates[i];
     }
     return sum;
 }
 
-int32_t pdbm_getRate(const PDBM pdbm, cindex_t dim, cindex_t clock)
+CostType pdbm_getRate(const PDBM pdbm, cindex_t dim, cindex_t clock)
 {
     assert(pdbm && dim && clock > 0 && clock < dim);
 
     return pdbm_rates(pdbm)[clock];
 }
 
-uint32_t pdbm_getCostAtOffset(const PDBM pdbm, cindex_t dim)
+CostType pdbm_getCostAtOffset(const PDBM pdbm, cindex_t dim)
 {
     assert(pdbm && dim);
 
     return pdbm_cost(pdbm);
 }
 
-void pdbm_setCostAtOffset(PDBM& pdbm, cindex_t dim, uint32_t value)
+void pdbm_setCostAtOffset(PDBM& pdbm, cindex_t dim, CostType value)
 {
     assert(pdbm && dim);
 
@@ -1140,13 +1366,13 @@ uint32_t pdbm_getUpperFacets(PDBM& pdbm, cindex_t dim, cindex_t* facets)
     return cnt;
 }
 
-int32_t pdbm_getCostOfValuation(const PDBM pdbm, cindex_t dim, const int32_t* valuation)
+CostType pdbm_getCostOfValuation(const PDBM pdbm, cindex_t dim, const int32_t* valuation)
 {
     assert(pdbm && dim && valuation);
     assert(pdbm_containsInt(pdbm, dim, valuation));
 
     raw_t* dbm = pdbm_matrix(pdbm);
-    int32_t cost = pdbm_cost(pdbm);
+    CostType cost = pdbm_cost(pdbm);
     for (uint32_t i = 1; i < dim; i++) {
         cost += (valuation[i] - -dbm_raw2bound(DBM(0, i))) * pdbm_rates(pdbm)[i];
     }
@@ -1179,13 +1405,19 @@ bool pdbm_isValid(const PDBM pdbm, cindex_t dim)
     }
 
     raw_t* dbm = pdbm_matrix(pdbm);
-    int32_t* rates = pdbm_rates(pdbm);
-    uint32_t cost = pdbm_cost(pdbm);
-    int32_t cache = pdbm_cache(pdbm);
-    int32_t inf = pdbm_infimum(dbm, dim, cost, rates);
+    CostType* rates = pdbm_rates(pdbm);
+    CostType cost = pdbm_cost(pdbm);
+    CostType cache = pdbm_cache(pdbm);
+    CostType inf = pdbm_infimum(dbm, dim, cost, rates);
+
+//    std::cout << "inf: " << inf << std::endl;
+//    std::cout << "(cache == INVALID || cache == inf): " << (cache == INVALID || cache == inf) << std::endl;
+//    std::cout << "dbm_isValid(dbm, dim): " << dbm_isValid(dbm, dim) << std::endl;
+//    std::cout << " rates[0] == 0: " << (rates[0] == 0) << std::endl;
+//    std::cout << "(!(pdbm_isUnbounded(pdbm, dim) && pdbm_getSlopeOfDelayTrajectory(pdbm, dim) < 0) || inf == -INFINITE_COST): " << (!(pdbm_isUnbounded(pdbm, dim) && pdbm_getSlopeOfDelayTrajectory(pdbm, dim) < 0) || inf == -INFINITE_COST) << std::endl;
 
     return (cache == INVALID || cache == inf) && dbm_isValid(dbm, dim) && rates[0] == 0 &&
-           (!pdbm_isUnbounded(pdbm, dim) || pdbm_getSlopeOfDelayTrajectory(pdbm, dim) >= 0);
+           (!(pdbm_isUnbounded(pdbm, dim) && pdbm_getSlopeOfDelayTrajectory(pdbm, dim) < 0) || inf == -INFINITE_COST);
 }
 
 void pdbm_freeClock(PDBM& pdbm, cindex_t dim, cindex_t clock)
@@ -1209,7 +1441,7 @@ void pdbm_getOffset(const PDBM pdbm, cindex_t dim, int32_t* valuation)
     }
 }
 
-void pdbm_setRate(PDBM& pdbm, cindex_t dim, cindex_t clock, int32_t rate)
+void pdbm_setRate(PDBM& pdbm, cindex_t dim, cindex_t clock, CostType rate)
 {
     assert(pdbm && dim && clock > 0 && clock < dim);
 
@@ -1239,7 +1471,7 @@ const raw_t* pdbm_getMatrix(const PDBM pdbm, cindex_t dim)
     return pdbm_matrix(pdbm);
 }
 
-const int32_t* pdbm_getRates(const PDBM pdbm, cindex_t dim)
+const CostType* pdbm_getRates(const PDBM pdbm, cindex_t dim)
 {
     assert(pdbm && dim);
     return pdbm_rates(pdbm);
@@ -1254,19 +1486,18 @@ bool pdbm_constrainToFacet(PDBM& pdbm, cindex_t dim, cindex_t i, cindex_t j)
 
 void pdbm_print(FILE* f, const PDBM pdbm, cindex_t dim)
 {
-    int32_t infimum = pdbm_getInfimum(pdbm, dim);
+    CostType infimum = pdbm_getInfimum(pdbm, dim);
     dbm_print(f, pdbm_matrix(pdbm), dim);
     fprintf(f, "Rates:");
     for (cindex_t i = 1; i < dim; i++) {
-        fprintf(f, " %d", pdbm_getRate(pdbm, dim, i));
+//        fprintf(f, " %f", pdbm_getRate(pdbm, dim, i));
     }
     fprintf(f, "\n");
-    fprintf(f, "Offset: %d Infimum: %d\n", pdbm_cost(pdbm), infimum);
 }
 
 void pdbm_print(std::ostream& o, const PDBM pdbm, cindex_t dim)
 {
-    int32_t infimum = pdbm_getInfimum(pdbm, dim);
+    CostType infimum = pdbm_getInfimum(pdbm, dim);
     dbm_cppPrint(o, pdbm_matrix(pdbm), dim);
     o << "Rates:";
     for (cindex_t i = 1; i < dim; i++) {
@@ -1293,7 +1524,7 @@ void pdbm_freeDown(PDBM& pdbm, cindex_t dim, cindex_t index)
 
     /* Move offset point.
      */
-    int32_t rate = pdbm_rates(pdbm)[index];
+    CostType rate = pdbm_rates(pdbm)[index];
     int32_t bound = -dbm_raw2bound(DBM(0, index));
     pdbm_cost(pdbm) -= bound * rate;
 
@@ -1304,7 +1535,7 @@ void pdbm_freeDown(PDBM& pdbm, cindex_t dim, cindex_t index)
 
 void pdbm_normalise(PDBM pdbm, cindex_t dim)
 {
-    int32_t* rates = pdbm_rates(pdbm);
+    CostType* rates = pdbm_rates(pdbm);
     cindex_t next[dim];
     cindex_t i;
 
@@ -1330,7 +1561,7 @@ void pdbm_normalise(PDBM pdbm, cindex_t dim)
 
 bool pdbm_hasNormalForm(PDBM pdbm, cindex_t dim)
 {
-    int32_t* rates = pdbm_rates(pdbm);
+    CostType* rates = pdbm_rates(pdbm);
     cindex_t next[dim];
     cindex_t i;
 
@@ -1358,3 +1589,51 @@ bool pdbm_hasNormalForm(PDBM pdbm, cindex_t dim)
 
     return true;
 }
+
+
+PDBM pdbm_from_dbm(const int32_t* dbm, cindex_t dim) {
+    PDBM pdbm = pdbm_allocate(dim);
+    memcpy(pdbm->data, dbm, dim * dim * sizeof(int32_t));
+    return pdbm;
+}
+
+void pdbm_intersect(PDBM& dst, const PDBM& src, cindex_t dim) {
+    assert(dst && src && dim > 0);
+
+    /* Compute the cost at the origin.
+     */
+    CostType cost = pdbm_cost(dst);
+    CostType* rates = pdbm_rates(dst);
+    raw_t* dbm = pdbm_matrix(dst);
+    for (uint32_t k = 1; k < dim; k++) {
+        cost += rates[k] * dbm_raw2bound(DBM(0, k));
+    }
+
+    if (!dbm_intersection(dbm, pdbm_matrix(src), dim)) {
+        // If there was no intersection, dst becomes empty.
+        pdbm_decRef(dst);
+        dst = nullptr;
+    } else {
+        pdbm_prepare(dst, dim);
+        // Compute the cost at the new offset point and invalidate the cache.
+        for (uint32_t k = 1; k < dim; k++) {
+            cost -= rates[k] * dbm_raw2bound(DBM(0, k));
+        }
+        pdbm_cost(dst) = cost;
+        pdbm_cache(dst) = INVALID;
+
+        assertx(pdbm_isValid(dst, dim));
+
+        pdbm_cache(dst) = INVALID;
+    }
+
+}
+
+void pdbm_simplify_rational_cost(PDBM& pdbm, cindex_t dim) {
+    pdbm_prepare(pdbm, dim);
+    pdbm_cost(pdbm) = pdbm_cost(pdbm).numerator() / pdbm_cost(pdbm).denominator();
+    for (cindex_t i = 1; i < dim; i++) {
+        pdbm_rates(pdbm)[i] = pdbm_rates(pdbm)[i].numerator() / pdbm_rates(pdbm)[i].denominator();
+    }
+    pdbm_cache(pdbm) = INVALID;
+}
\ No newline at end of file
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index c356200..2ff3041 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -3,13 +3,12 @@
 set(libs UDBM UUtils::udebug UUtils::base UUtils::hash m)
 
 if (TESTING)
-  file(GLOB test_sources "test_*.c" "test_*.cpp")
+  add_executable(test_dbm test_dbm.c)
+  target_link_libraries(test_dbm ${libs})
 
-  foreach(source ${test_sources})
-    get_filename_component(test_target ${source} NAME_WE)
-    add_executable(${test_target} ${source})
-    target_link_libraries(${test_target} ${libs})
-  endforeach()
+  add_executable(test_pfed test_pfed.cpp)
+  target_compile_definitions(test_pfed PRIVATE DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN)
+  target_link_libraries(test_pfed PRIVATE ${libs} doctest::doctest)
 
   # comments contain expected time of Linux debug build, YMMV, whereas timeouts are worst case
   add_test(NAME dbm_dbm_1_10 COMMAND test_dbm 1 10)
diff --git a/test/test_pfed.cpp b/test/test_pfed.cpp
new file mode 100644
index 0000000..643afab
--- /dev/null
+++ b/test/test_pfed.cpp
@@ -0,0 +1,79 @@
+#include "dbm/cost_type.h"
+#include "dbm/pfed.h"
+
+#include <doctest/doctest.h>
+
+// Asserts that the valuation is contained in some pdbm in the pfed and that the cost of that valuation is the expected cost
+void any_contains_with_cost(const dbm::pfed_t& pfed, int* valuation, size_t dim, CostType expected) {
+    bool any_contains = false;
+    for (const auto& pdbm : pfed) {
+        if (pdbm_containsInt(pdbm, dim, valuation)) {
+            CHECK_MESSAGE((pdbm_getCostOfValuation(pdbm, dim, valuation) == expected), "Expected cost " << expected << " but got " << pdbm_getCostOfValuation(pdbm, dim, valuation));
+            any_contains = true;
+        }
+    }
+    CHECK(any_contains);
+}
+
+TEST_SUITE("pfed")
+{
+    TEST_CASE("delay_constraint")
+    {
+        dbm::pfed_t pfed(2);
+        pfed.setZero();
+        pfed.up(4);
+        pfed.constrain(1, 0, 2, false);
+        pfed.constrain(0, 1, -1, false);
+
+        REQUIRE((pfed.getInfimum() == 4));
+    }
+
+    TEST_CASE("reset")
+    {
+        dbm::pfed_t pfed(3);
+        pfed.setZero();
+        pfed.up(1);
+        pfed.constrain(1, 0, 3, false);
+        pfed.updateValue(1, 0);
+        pfed.up(2);
+        pfed.constrain(0, 2, -1, false);
+        pfed.updateValue(2, 0);
+        pfed.up(3);
+
+        REQUIRE((pfed.size() == 2));
+        REQUIRE((pfed.getInfimum() == 1));
+
+        // Assert that some pdbm in pfed contains the valuation (1,1) and has cost 4
+        {
+            int valuation[3] = {0, 1, 1};
+            any_contains_with_cost(pfed, valuation, 3, 4);
+        }
+
+        {
+            int valuation[3] = {0, 2, 2};
+            any_contains_with_cost(pfed, valuation, 3, 7);
+        }
+
+        {
+            int valuation[3] = {0, 2, 0};
+            any_contains_with_cost(pfed, valuation, 3, 4);
+        }
+    }
+
+    TEST_CASE("set_uniform_cost") {
+        dbm::pfed_t pfed(3);
+        pfed.setZero();
+        pfed.up(1);
+        pfed.setUniformCost(3);
+        // All valuations in the zone should now have cost 3
+        {
+            int valuation[3] = {0, 0, 0};
+            any_contains_with_cost(pfed, valuation, 3, 3);
+        }
+
+        {
+            int valuation[3] = {0, 5, 5};
+            any_contains_with_cost(pfed, valuation, 3, 3);
+        }
+    }
+}
\ No newline at end of file