Normalize rate matrix in DecisionDiagramModel for DTMC, so that it re…

…presents

the transition probability matrix.

cudd/cudd/cudd.h

@@ -981,7 +981,7 @@ extern int Cudd_NextCube (DdGen *gen, int **cube, CUDD_VALUE_TYPE *value);
 extern DdGen * Cudd_FirstPrime(DdManager *dd, DdNode *l, DdNode *u, int **cube);
 extern int Cudd_NextPrime(DdGen *gen, int **cube);
 extern DdNode * Cudd_bddComputeCube (DdManager *dd, DdNode * const *vars, int *phase, int n);
-extern DdNode * Cudd_addComputeCube (DdManager *dd, DdNode **vars, int *phase, int n);
+extern DdNode * Cudd_addComputeCube (DdManager *dd, DdNode * const *vars, int *phase, int n);
 extern DdNode * Cudd_CubeArrayToBdd (DdManager *dd, int *array);
 extern int Cudd_BddToCubeArray (DdManager *dd, DdNode *cube, int *array);
 extern DdGen * Cudd_FirstNode (DdManager *dd, DdNode *f, DdNode **node);

cudd/cudd/cuddUtil.c

@@ -2186,7 +2186,7 @@ Cudd_bddComputeCube(
 DdNode *
 Cudd_addComputeCube(
   DdManager * dd,
-  DdNode ** vars,
+  DdNode * const * vars,
   int * phase,
   int  n)
 {

models.cc

@@ -611,7 +611,12 @@ DecisionDiagramModel DecisionDiagramModel::Make(
     // TODO(hlsyounes): implement.
     LOG(FATAL) << "not implemented";
   }
-  // TODO(hlsyounes): normalize rate matrix for DTMC.
+  if (model.type() == CompiledModelType::DTMC) {
+    // Normalize rate matrix.
+    ADD col_cube = manager->GetCube(
+        manager->GetAddVariableArray(1, 2, manager->GetVariableCount()));
+    rate_matrix = rate_matrix / rate_matrix.ExistAbstract(col_cube);
+  }
 
   std::cout << manager->GetVariableCount() << " variables." << std::endl;
 

src/ddutil.cc

@@ -153,17 +153,25 @@ DecisionDiagram& DecisionDiagram::operator=(const DecisionDiagram& dd) {
 
 DecisionDiagram::~DecisionDiagram() { Deref(); }
 
+void DecisionDiagram::Ref(DdNode* node) {
+  if (node) {
+    Cudd_Ref(node);
+  }
+}
+
 void DecisionDiagram::Ref() {
-  if (node_) {
-    Cudd_Ref(node_);
+  Ref(node_);
+}
+
+void DecisionDiagram::Deref(DdManager* manager, DdNode* node) {
+  if (node) {
+    Cudd_RecursiveDeref(manager, node);
   }
 }
 
 void DecisionDiagram::Deref() {
-  if (node_) {
-    Cudd_RecursiveDeref(manager_, node_);
-    node_ = nullptr;
-  }
+  Deref(manager_, node_);
+  node_ = nullptr;
 }
 
 bool DecisionDiagram::IsConstant() const { return Cudd_IsConstant(node_); }
@@ -270,6 +278,10 @@ BDD ADD::StrictThreshold(double threshold) const {
              Cudd_addBddStrictThreshold(manager(), node(), threshold));
 }
 
+ADD ADD::ExistAbstract(const ADD& cube) const {
+  return ADD(manager(), Cudd_addExistAbstract(manager(), node(), cube.node()));
+}
+
 ADD ADD::operator-() const { return MonadicApply(AddNegate, *this); }
 
 ADD ADD::operator+(const ADD& dd) const {
@@ -403,11 +415,26 @@ VariableArray<BDD> DecisionDiagramManager::GetBddVariableArray(int start,
   return VariableArray<BDD>(variables);
 }
 
+VariableArray<ADD> DecisionDiagramManager::GetAddVariableArray(int start,
+                                                               int incr,
+                                                               int end) const {
+  std::vector<ADD> variables;
+  for (int i = start; i < end; i += incr) {
+    variables.push_back(GetAddVariable(i));
+  }
+  return VariableArray<ADD>(variables);
+}
+
 BDD DecisionDiagramManager::GetCube(const VariableArray<BDD>& variables) const {
   return BDD(manager_, Cudd_bddComputeCube(manager_, variables.data(), nullptr,
                                            variables.size()));
 }
 
+ADD DecisionDiagramManager::GetCube(const VariableArray<ADD>& variables) const {
+  return ADD(manager_, Cudd_addComputeCube(manager_, variables.data(), nullptr,
+                                           variables.size()));
+}
+
 ODD::ODD(const OddNode* root, int node_count)
     : root_(root), node_count_(node_count) {}
 

src/ddutil.h

@@ -76,7 +76,9 @@ class DecisionDiagram {
   DdNode* node() const { return node_; }
 
  private:
+  static void Ref(DdNode* node);
   void Ref();
+  static void Deref(DdManager* manager, DdNode* node);
   void Deref();
 
   // The decision-diagram manager.
@@ -179,6 +181,10 @@ class ADD : public DecisionDiagram {
   // Returns the BDD representing *this > threshold.
   BDD StrictThreshold(double threshold) const;
 
+  // Returns the ADD that abstracts through summation all the variables in cube
+  // from this ADD.
+  ADD ExistAbstract(const ADD& cube) const;
+
   // Arithmetic operators for ADDs.
   ADD operator-() const;
   ADD operator+(const ADD& dd) const;
@@ -247,13 +253,16 @@ class VariableArray {
  public:
   VariableArray();
 
+  ~VariableArray();
+
   std::vector<DdNode*>::size_type size() const { return variables_.size(); }
 
  private:
   VariableArray(std::vector<DD> variables);
 
   DdNode* const* data() const { return variables_.data(); }
 
+  DdManager* manager_;
   std::vector<DdNode*> variables_;
 
   friend class DecisionDiagramManager;
@@ -304,9 +313,16 @@ class DecisionDiagramManager {
   // end (exclusive).
   VariableArray<BDD> GetBddVariableArray(int start, int incr, int end) const;
 
+  // Returns an array of every incr ADD variables between start (inclusive) and
+  // end (exclusive).
+  VariableArray<ADD> GetAddVariableArray(int start, int incr, int end) const;
+
   // Returns the positive unate cube of the given variables.
   BDD GetCube(const VariableArray<BDD>& variables) const;
 
+  // Returns the positive unate cube of the given variables.
+  ADD GetCube(const VariableArray<ADD>& variables) const;
+
   // TODO(hlsyounes): remove once all code is using wrapper classes.
   DdManager* manager() const { return manager_; }
 
@@ -350,13 +366,26 @@ class ODD {
 int Log2(int n);
 
 template <typename DD>
-VariableArray<DD>::VariableArray() = default;
+VariableArray<DD>::VariableArray()
+    : manager_(nullptr) {}
 
 template <typename DD>
-VariableArray<DD>::VariableArray(std::vector<DD> variables) {
+VariableArray<DD>::VariableArray(std::vector<DD> variables)
+    : manager_(nullptr) {
   variables_.reserve(variables.size());
   for (const DD& dd : variables) {
+    if (manager_ == nullptr) {
+      manager_ = dd.manager();
+    }
     variables_.push_back(dd.node());
+    DecisionDiagram::Ref(dd.node());
+  }
+}
+
+template <typename DD>
+VariableArray<DD>::~VariableArray() {
+  for (DdNode* node : variables_) {
+    DecisionDiagram::Deref(manager_, node);
   }
 }
 

src/ddutil_test.cc

@@ -672,6 +672,26 @@ TEST(DecisionDiagramTest, ComputesCubes) {
   EXPECT_FALSE(dd4.ValueInState({true, false}));
   EXPECT_FALSE(dd4.ValueInState({false, true}));
   EXPECT_FALSE(dd4.ValueInState({false, false}));
+  const ADD dd5 = manager.GetCube(VariableArray<ADD>());
+  EXPECT_EQ(1, dd5.ValueInState({true, true}));
+  EXPECT_EQ(1, dd5.ValueInState({true, false}));
+  EXPECT_EQ(1, dd5.ValueInState({false, true}));
+  EXPECT_EQ(1, dd5.ValueInState({false, false}));
+  const ADD dd6 = manager.GetCube(manager.GetAddVariableArray(0, 1, 1));
+  EXPECT_EQ(1, dd6.ValueInState({true, true}));
+  EXPECT_EQ(1, dd6.ValueInState({true, false}));
+  EXPECT_EQ(0, dd6.ValueInState({false, true}));
+  EXPECT_EQ(0, dd6.ValueInState({false, false}));
+  const ADD dd7 = manager.GetCube(manager.GetAddVariableArray(1, 1, 2));
+  EXPECT_EQ(1, dd7.ValueInState({true, true}));
+  EXPECT_EQ(0, dd7.ValueInState({true, false}));
+  EXPECT_EQ(1, dd7.ValueInState({false, true}));
+  EXPECT_EQ(0, dd7.ValueInState({false, false}));
+  const ADD dd8 = manager.GetCube(manager.GetAddVariableArray(0, 1, 2));
+  EXPECT_EQ(1, dd8.ValueInState({true, true}));
+  EXPECT_EQ(0, dd8.ValueInState({true, false}));
+  EXPECT_EQ(0, dd8.ValueInState({false, true}));
+  EXPECT_EQ(0, dd8.ValueInState({false, false}));
 }
 
 TEST(DecisionDiagramTest, BddValueInState) {