Force root state to be equal to the root of the subtree as modeled by _zz

Author: melkebir

src/ilppmhctisolver.cpp

@@ -110,7 +110,58 @@ void IlpPmhCtiSolver::processSolution()
   const int nrAnatomicalSites = _anatomicalSiteToIndex.size();
   const int nrSamples = _F.getNrSamples();
   const int nrNodes = _indexToNode.size();
+//  const int nrArcs = _indexToArc.size();
 
+//  for (int ij = 0; ij < nrArcs; ++ij)
+//  {
+//    for (int s = 0; s < nrAnatomicalSites; ++s)
+//    {
+//      int size_L_s = _L[s].size();
+//      for (int c = 0; c < size_L_s; ++c)
+//      {
+//        for (int t = 0; t < nrAnatomicalSites; ++t)
+//        {
+//          int size_L_t = _L[t].size();
+//          for (int d = 0; d < size_L_t; ++d)
+//          {
+//            if (_xx[ij][s][c][t][d].get(GRB_DoubleAttr_X) >= 0.4)
+//            {
+//              std::cout << _xx[ij][s][c][t][d].get(GRB_StringAttr_VarName)
+//              << " = " << _xx[ij][s][c][t][d].get(GRB_DoubleAttr_X) << std::endl;
+//            }
+//          }
+//        }
+//      }
+//    }
+//  }
+//  
+//  std::cout << std::endl;
+//  
+//  for (int i = 0; i < nrNodes; ++i)
+//  {
+//    for (int s = 0; s < nrAnatomicalSites; ++s)
+//    {
+//      int size_L_s = _L[s].size();
+//      for (int c = 0; c < size_L_s; ++c)
+//      {
+//        for (int t = 0; t < nrAnatomicalSites; ++t)
+//        {
+//          int size_L_t = _L[t].size();
+//          for (int d = 0; d < size_L_t; ++d)
+//          {
+//            if (_zz[i][s][c][t][d].get(GRB_DoubleAttr_X) >= 0.4)
+//            {
+//              std::cout << _zz[i][s][c][t][d].get(GRB_StringAttr_VarName)
+//              << " = " << _zz[i][s][c][t][d].get(GRB_DoubleAttr_X) << std::endl;
+//            }
+//          }
+//        }
+//      }
+//    }
+//  }
+//
+//  std::cout << std::endl;
+//  
 //  for (int i = 0; i < nrNodes; ++i)
 //  {
 //    for (int s = 0; s < nrAnatomicalSites; ++s)
@@ -572,7 +623,9 @@ void IlpPmhCtiSolver::initConstraints()
 //  {
 //    std::cout << i << " : " << getLabel(_indexToNode[i]) << std::endl;
 //  }
-//  _model.addConstr(_x[28][nrAnatomicalSites][0] == 1);
+//  _model.addConstr(_x[21][0][0] == 1);
+//  _model.addConstr(_x[21][2][0] == 1);
+//  _model.addConstr(_x[21][3][0] == 1);
 }
 
 void IlpPmhCtiSolver::initSingleSourceSeedingConstraints()

src/ilppmhprsolver.cpp

@@ -552,6 +552,56 @@ void IlpPmhPrSolver::initConstraintsG()
     sum2.clear();
   }
 
+  for (NodeIt v_i(getTree()); v_i != lemon::INVALID; ++v_i)
+  {
+    if (isLeaf(v_i)) continue;
+    const int i = (*_pNodeToIndex)[v_i];
+    
+    for (int t = 0; t < nrAnatomicalSites; ++t)
+    {
+      const int size_L_t = _L[t].size();
+      for (int d = 0; d < size_L_t; ++d)
+      {
+        for (int s = 0; s < nrAnatomicalSites; ++s)
+        {
+          const int size_L_s = _L[s].size();
+          for (int c = 0; c < size_L_s; ++c)
+          {
+            sum += _zz[i][s][c][t][d];
+          }
+        }
+        _model.addConstr(sum == _x[i][t][d] - _r[i][t][d]);
+        sum.clear();
+      }
+    }
+  }
+  
+  // one arc in G
+  for (ArcIt a_ij(getTree()); a_ij != lemon::INVALID; ++a_ij)
+  {
+    Node v_j = getTree().target(a_ij);
+    if (isLeaf(v_j)) continue;
+    const int ij = (*_pArcToIndex)[a_ij];
+    for (int t = 0; t < nrAnatomicalSites; ++t)
+    {
+      const int size_L_t = _L[t].size();
+      for (int d = 0; d < size_L_t; ++d)
+      {
+        for (int s = 0; s < nrAnatomicalSites; ++s)
+        {
+          const int size_L_s = _L[s].size();
+          for (int c = 0; c < size_L_s; ++c)
+          {
+            sum += _xx[ij][s][c][t][d];
+          }
+        }
+        
+//        _model.addConstr(sum <= 1);
+        sum.clear();
+      }
+    }
+  }
+  
   // relate _zz to _z
   for (NodeIt v_i(getTree()); v_i != lemon::INVALID; ++v_i)
   {

src/ilppmhsolver.cpp

@@ -487,8 +487,8 @@ bool IlpPmhSolver::solve(int nrThreads, int timeLimit)
     else if (status == GRB_INFEASIBLE)
     {
       //std::cerr << "Model is infeasible" << std::endl;
-//      _model.computeIIS();
-//      _model.write("/tmp/model_IIS.ilp");
+      _model.computeIIS();
+      _model.write("/tmp/model_IIS.ilp");
       return false;
     }
     else if (status == GRB_UNBOUNDED)