Hot Start

Author: melkebir

src/ilppmhprsolver.cpp

@@ -734,28 +734,28 @@ GRBLinExpr IlpPmhPrSolver::initObjective(const IntTriple& bounds)
   obj += migrationNumber
       + (1. / (nrAnatomicalSites * nrAnatomicalSites)) * comigrationNumber;
 
-  double factor = (1. / (nrAnatomicalSites * nrAnatomicalSites))
-                * (1. / (nrNodes * nrNodes));
-  
-  for (NodeIt v_i(getTree()); v_i != lemon::INVALID; ++v_i)
-  {
-    const int i = (*_pNodeToIndex)[v_i];
-    for (int s = 0; s < nrAnatomicalSites; ++s)
-    {
-      const int size_L_s = _L[s].size();
-      for (int c = 0; c < size_L_s; ++c)
-      {
-        for (int t = 0; t < nrAnatomicalSites; ++t)
-        {
-          const int size_L_t = _L[t].size();
-          for (int d = 0; d < size_L_t; ++d)
-          {
-            obj += factor * _x[i][s][c];
-          }
-        }
-      }
-    }
-  }
+//  double factor = (1. / (nrAnatomicalSites * nrAnatomicalSites))
+//                * (1. / (nrNodes * nrNodes));
+//
+//  for (NodeIt v_i(getTree()); v_i != lemon::INVALID; ++v_i)
+//  {
+//    const int i = (*_pNodeToIndex)[v_i];
+//    for (int s = 0; s < nrAnatomicalSites; ++s)
+//    {
+//      const int size_L_s = _L[s].size();
+//      for (int c = 0; c < size_L_s; ++c)
+//      {
+//        for (int t = 0; t < nrAnatomicalSites; ++t)
+//        {
+//          const int size_L_t = _L[t].size();
+//          for (int d = 0; d < size_L_t; ++d)
+//          {
+//            obj += factor * _x[i][s][c];
+//          }
+//        }
+//      }
+//    }
+//  }
 
   GRBLinExpr seedingSiteNumber;
   for (int s = 0; s < nrAnatomicalSites; ++s)
@@ -768,9 +768,9 @@ GRBLinExpr IlpPmhPrSolver::initObjective(const IntTriple& bounds)
     _model.addConstr(seedingSiteNumber <= bounds.second.second);
   }
 
-  factor = (1. / (nrAnatomicalSites * nrAnatomicalSites))
-         * (1. / (nrNodes * nrNodes))
-         * (1. / (nrAnatomicalSites + 1));
+  double factor = (1. / (nrAnatomicalSites * nrAnatomicalSites))
+//         * (1. / (nrNodes * nrNodes))
+                * (1. / (nrAnatomicalSites + 1));
 
   obj += factor * seedingSiteNumber;
 

src/ilppmhsolver.cpp

@@ -64,6 +64,33 @@ void IlpPmhSolver::init(const IntTriple& bounds)
   initLeafConstraints();
   initObjective(bounds);
   initCallbacks();
+  initWarmStart();
+}
+
+void IlpPmhSolver::initWarmStart()
+{
+  const int nrAnatomicalSites = _anatomicalSiteToIndex.size();
+  for (NodeIt v_i(getTree()); v_i != lemon::INVALID; ++v_i)
+  {
+    if (isLeaf(v_i)) continue;
+    
+    const int i = (*_pNodeToIndex)[v_i];
+    for (int s = 0; s < nrAnatomicalSites; ++s)
+    {
+      const int size_L_s = _L[s].size();
+      for (int c = 0; c < size_L_s; ++c)
+      {
+        if (s == _primaryIndex && c == 0)
+        {
+          _x[i][s][c].set(GRB_DoubleAttr_Start, 1);
+        }
+        else
+        {
+          _x[i][s][c].set(GRB_DoubleAttr_Start, 0);
+        }
+      }
+    }
+  }
 }
 
 void IlpPmhSolver::initIndices()

src/ilppmhsolver.h

@@ -155,6 +155,9 @@ class IlpPmhSolver
   /// Initialize ILP S constraints
   virtual void initSingleSourceSeedingConstraints();
 
+  /// Initialize initial solution
+  virtual void initWarmStart();
+  
   /// Initialize ILP objective function
   ///
   /// @param bounds Upper bounds on mu, gamma and sigma