zone.hh

#pragma once
 
#include <cmath>
#include <limits>
#include <algorithm>
#include <unordered_map>
#include <boost/unordered_map.hpp>
#include <boost/range/adaptor/reversed.hpp>
 
#include "common_types.hh"
#include "bounds.hh"
#include "constraint.hh"
 
#include <Eigen/Core>
#include <utility>
 
namespace learnta {
 
  struct Zone {
    Eigen::Matrix <Bounds, Eigen::Dynamic, Eigen::Dynamic> value;
    Bounds M;
    std::vector<double> maxConstraints;
 
    Zone() = default;
 
    explicit Zone(const Eigen::Matrix <Bounds, Eigen::Dynamic, Eigen::Dynamic> &value) :
            Zone(Eigen::Matrix<Bounds, Eigen::Dynamic, Eigen::Dynamic>(value)) {}
 
    explicit Zone(Eigen::Matrix <Bounds, Eigen::Dynamic, Eigen::Dynamic> &&value) : value(std::move(value)) {}
 
    Zone(Eigen::Matrix <Bounds, Eigen::Dynamic, Eigen::Dynamic> value, Bounds m) :
            value(std::move(value)), M(std::move(m)) {
      maxConstraints.resize(this->value.cols() - 1);
      std::fill(maxConstraints.begin(), maxConstraints.end(), M.first);
    }
 
    explicit Zone(const std::vector<double> &valuation, Bounds M) : M(std::move(M)) {
      value.resize(valuation.size() + 1, valuation.size() + 1);
      value.fill(Bounds(std::numeric_limits<double>::max(), false));
      for (std::size_t i = 0; i < valuation.size(); ++i) {
        this->tighten(i, -1, Bounds{valuation.at(i), true});
        this->tighten(-1, i, Bounds{-valuation.at(i), true});
      }
      maxConstraints.resize(this->value.cols() - 1);
      std::fill(maxConstraints.begin(), maxConstraints.end(), this->M.first);
    }
 
    [[nodiscard]] inline std::size_t getNumOfVar() const noexcept {
      return value.cols() - 1;
    }
 
    static Zone zero(int size) {
      static Zone zeroZone;
      if (zeroZone.value.cols() == size) {
        return zeroZone;
      }
      zeroZone.value.resize(size, size);
      zeroZone.value.fill(Bounds(0, true));
      return zeroZone;
    }
 
    static Zone top(std::size_t size) {
      static Zone topZone;
      if (static_cast<std::size_t>(topZone.value.cols()) == size) {
        return topZone;
      }
      topZone.value.resize(size, size);
      topZone.value.fill(Bounds(std::numeric_limits<double>::max(), false));
      return topZone;
    }
 
    void tighten(ClockVariables x, ClockVariables y, Bounds c) {
      x++;
      y++;
      value(x, y) = std::min(value(x, y), c);
      close1(x);
      close1(y);
    }
 
    void tighten(const Constraint &constraint) {
      switch (constraint.odr) {
        case learnta::Constraint::Order::ge:
          // lower bound
          this->tighten(-1, constraint.x, Bounds{-constraint.c, true});
          break;
        case learnta::Constraint::Order::gt:
          // lower bound
          this->tighten(-1, constraint.x, Bounds{-constraint.c, false});
          break;
        case learnta::Constraint::Order::le:
          // upper bound
          this->tighten(constraint.x, -1, Bounds{constraint.c, true});
          break;
        case learnta::Constraint::Order::lt:
          // upper bound
          this->tighten(constraint.x, -1, Bounds{constraint.c, false});
          break;
      }
    }
 
    void tighten(const std::vector<Constraint> &constraints) {
      for (const auto &constraint: constraints) {
        this->tighten(constraint);
      }
    }
 
    void applyResets(const std::vector<std::pair<ClockVariables, std::variant<double, ClockVariables>>> &resets) {
      // We apply renaming first
      for (const auto &[resetVariable, updatedVariable]: resets) {
        if (updatedVariable.index() == 1 && resetVariable != std::get<ClockVariables>(updatedVariable)) {
          this->unconstrain(resetVariable);
          this->value(resetVariable + 1, std::get<ClockVariables>(updatedVariable) + 1) = Bounds{0.0, true};
          this->value(std::get<ClockVariables>(updatedVariable) + 1, resetVariable + 1) = Bounds{0.0, true};
          canonize();
        }
      }
      // Then, assign a value
      for (const auto &[resetVariable, updatedVariable]: resets) {
        if (updatedVariable.index() == 0) {
          this->unconstrain(resetVariable);
          this->value(0, resetVariable + 1) = Bounds(-std::get<double>(updatedVariable), true);
          this->value(resetVariable + 1, 0) = Bounds(std::get<double>(updatedVariable), true);
          canonize();
        }
      }
    }
 
    void revertResets(const std::vector<std::pair<ClockVariables, std::variant<double, ClockVariables>>> &resets) {
      std::vector<ClockVariables> resetVariables;
      resetVariables.reserve(resets.size());
      std::unordered_map<ClockVariables, ClockVariables> reverseAssignments;
      for (const auto &[resetVariable, updatedVariable]: resets) {
        resetVariables.push_back(resetVariable);
        if (updatedVariable.index() == 1 && resetVariable != std::get<ClockVariables>(updatedVariable)) {
          reverseAssignments[std::get<ClockVariables>(updatedVariable)] = resetVariable;
        }
      }
      for (const auto &resetVariable: boost::adaptors::reverse(resetVariables)) {
        this->unconstrain(resetVariable);
        auto it = reverseAssignments.find(resetVariable);
        if (it != reverseAssignments.end()) {
          this->value(it->first + 1, it->second + 1) = Bounds{0.0, true};
          this->value(it->second + 1, it->first + 1) = Bounds{0.0, true};
        }
        canonize();
      }
    }
 
    Zone operator&&(const Zone &another) const {
      assert(this->value.cols() == another.value.cols());
      assert(this->value.rows() == another.value.rows());
      auto result = Zone{this->value.array().min(another.value.array())};
      result.M = this->M;
      result.canonize();
 
      return result;
    }
 
    Zone operator&=(const Zone &another) {
      assert(this->value.cols() == another.value.cols());
      assert(this->value.rows() == another.value.rows());
      this->value = this->value.cwiseMin(another.value);
      canonize();
 
      return *this;
    }
 
    Zone operator^(const Zone &another) const {
      Zone result;
      const size_t N = this->getNumOfVar();
      const size_t O = another.getNumOfVar();
      result = Zone::top(N + O);
      // Copy \f$\mathbb{T}'\f$
      result.value.block(N, N, O, O) = another.value.block(1, 1, O, O);
      result.value.block(0, N + 1, 1, O - 1) = another.value.block(0, 1, 1, O - 1);
      result.value.block(N + 1, 0, O - 1, 1) = another.value.block(1, 0, O - 1, 1);
      // Copy \f$\mathbb{T}\f$
      result.value.block(0, 0, N + 1, N + 1) = this->value;
      result.canonize();
 
      return result;
    }
 
    [[nodiscard]] std::vector<double> sample() {
      this->canonize();
      assert(this->isSatisfiableNoCanonize());
      std::vector<double> valuation;
      std::size_t N = this->getNumOfVar();
      valuation.resize(N);
      for (std::size_t i = 0; i < N; i++) {
        Bounds lowerBound = this->value(0, i + 1);
        Bounds upperBound = this->value(i + 1, 0);
        if (isPoint(upperBound, lowerBound)) {
          valuation[i] = upperBound.first;
        } else {
          double lower = std::max(0.0, -lowerBound.first);
          double upper = upperBound.first;
          for (std::size_t j = 0; j < i; j++) {
            Bounds tmpLowerBound = this->value(j + 1, i + 1);
            Bounds tmpUpperBound = this->value(i + 1, j + 1);
            lower = std::max(lower, -tmpLowerBound.first + valuation[j]);
            upper = std::min(upper, tmpUpperBound.first + valuation[j]);
          }
          assert(lower <= upper);
          if (lowerBound.second) {
            valuation[i] = lower;
          } else if (upper - lower > 0.5) {
            valuation[i] = lower + 0.5;
          } else {
            valuation[i] = (lower + upper) * 0.5;
          }
        }
      }
 
      return valuation;
    }
 
    void close1(ClockVariables x) {
      for (int i = 0; i < value.rows(); i++) {
        value.row(i) = value.row(i).array().min(value.row(x).array() + value(i, x));
        //      for (int j = 0; j < value.cols(); j++) {
        //        value(i, j) = std::min(value(i, j), value(i, x) + value(x, j));
        //      }
      }
    }
 
    void reset(ClockVariables x) {
      // 0 is the special variable here
      x++;
      value(0, x) = Bounds(-0.0, true);
      value(x, 0) = Bounds(0.0, true);
      value.col(x).tail(value.rows() - 1) = value.col(0).tail(value.rows() - 1);
      value.row(x).tail(value.cols() - 1) = value.row(0).tail(value.cols() - 1);
    }
 
    void unconstrain(ClockVariables x) {
      // 0 is the special variable here
      x++;
      value.col(x).fill(Bounds(std::numeric_limits<double>::max(), false));
      value.row(x).fill(Bounds(std::numeric_limits<double>::max(), false));
    }
 
    void elapse() {
      static constexpr Bounds infinity = Bounds(std::numeric_limits<double>::max(), false);
      value.col(0).fill(Bounds(infinity));
    }
 
    void reverseElapse() {
      value.row(0).fill(Bounds(0, true));
    }
 
    void canonize() {
      for (int k = 0; k < value.cols(); k++) {
        close1(k);
      }
    }
 
    bool isSatisfiable() {
      canonize();
      return this->isSatisfiableNoCanonize();
    }
 
    [[nodiscard]] bool isSatisfiableNoCanonize() const {
      return (value + value.transpose()).minCoeff() >= Bounds(0.0, true);
    }
 
    explicit operator bool() {
      return isSatisfiable();
    }
 
    void abstractize() {
      static constexpr Bounds infinity = Bounds(std::numeric_limits<double>::max(), false);
      for (auto it = value.data(); it < value.data() + value.size(); it++) {
        if (*it > Bounds{M.first, true}) {
          *it = Bounds(infinity);
        } else if (*it < Bounds{-M.first, false}) {
          *it = Bounds(-M.first, false);
        }
      }
    }
 
    void extrapolate() {
      static constexpr Bounds infinity = Bounds(std::numeric_limits<double>::max(), false);
      for (std::size_t i = 0; i < this->maxConstraints.size(); ++i) {
        if (value(i + 1, 0).first > this->maxConstraints.at(i)) {
          value(i + 1, 0) = infinity;
        }
        if (-value(0, i + 1).first > this->maxConstraints.at(i)) {
          value(0, i + 1) = Bounds{-this->maxConstraints.at(i), false};
        }
        for (std::size_t j = 0; j < this->maxConstraints.size(); ++j) {
          if (value(i + 1, j + 1).first > this->maxConstraints.at(i)) {
            value(i + 1, j + 1) = infinity;
          } else if (-value(0, i + 1).first > this->maxConstraints.at(i)) {
            value(i + 1, j + 1) = infinity;
          } else if (-value(0, j + 1).first > this->maxConstraints.at(j)) {
            value(i + 1, j + 1) = infinity;
          }
        }
      }
    }
 
    void makeUnsat() {
      value(0, 0) = Bounds(-std::numeric_limits<double>::infinity(), false);
    }
 
    [[nodiscard]] bool includes(const Zone &zone) const {
      return this->value.cwiseMax(zone.value) == this->value;
    };
 
    bool operator==(const Zone &z) const {
      return value.cols() == z.value.cols() && value == z.value;
    }
 
    [[nodiscard]] bool equalIgnoreZero(Zone z) const {
      z.value(0, 0) = value(0, 0);
      return value == z.value;
    }
 
    [[nodiscard]] bool strictEqual(Zone z) const {
      z.value.diagonal() = value.diagonal();
      return value == z.value;
    }
  };
 
  static inline std::ostream &print(std::ostream &os, const learnta::Zone &zone) {
    for (int i = 0; i < zone.value.cols(); ++i) {
      for (int j = 0; j < zone.value.rows(); ++j) {
        print(os, zone.value(i, j));
        os << " ";
      }
      os << "\n";
    }
 
    return os;
  }
 
  static inline std::ostream &operator<<(std::ostream &os, const learnta::Zone &zone) {
    return learnta::print(os, zone);
  }
 
  inline std::size_t hash_value(learnta::Zone const &zone) {
    std::size_t seed = zone.value.array().size();
    const auto asVector = zone.value.array();
 
    union DI {
      double asD;
      uint64_t asI;
    };
 
    for (auto it = asVector.data(); it != asVector.data() + asVector.size(); it++) {
      DI value{};
      value.asD = it->first;
      seed ^= it->second + value.asI + 0x9e3779b9 + (seed << 6) + (seed >> 2);
    }
    return seed;
  }
}