Saving enumerated mutation trees in a single file

Author: melkebir

src/enumeratecanonicalclonetrees.cpp

@@ -37,7 +37,7 @@ void EnumerateCanonicalCloneTrees::enumerate(const std::string& outputDirectory,
 
   gm::RootedCladisticNoisySparseEnumeration enumerate(G, limit, -1,
                                                       nrThreads,
-                                                      _F.getNrCharacters(),
+                                                      1, //_F.getNrCharacters(),
                                                       true,
                                                       false,
                                                       IntSet());

src/enumeratemutationtrees.cpp

@@ -17,6 +17,7 @@ EnumerateMutationTrees::EnumerateMutationTrees(const FrequencyMatrix& F)
 void EnumerateMutationTrees::enumerate(const std::string& outputDirectory,
                                        const int nrThreads,
                                        const int limit,
+                                       const int timeLimit,
                                        TreeVector& mutationTrees) const
 {
   // 1. Construct frequency tensors F_lb and F_ub
@@ -35,15 +36,17 @@ void EnumerateMutationTrees::enumerate(const std::string& outputDirectory,
   G.init();
   G.setLabels(F_lb);
 
-  gm::RootedCladisticNoisyEnumeration enumerate(G, limit, -1,
+  gm::RootedCladisticNoisyEnumeration enumerate(G, limit, timeLimit,
                                                 nrThreads,
-                                                _F.getNrCharacters(),
+                                                1, //_F.getNrCharacters(),
                                                 true,
                                                 false,
                                                 IntSet());
 
   enumerate.run();
 
+  int size = enumerate.objectiveValue();
+  
   gm::SolutionSet solutionSet;
   enumerate.populateSolutionSet(solutionSet);
 
@@ -112,7 +115,7 @@ void EnumerateMutationTrees::enumerate(const std::string& outputDirectory,
     }
   }
 
-  std::cerr << "Found " << mutationTrees.size() << " mutation trees." << std::endl;
+  std::cerr << "Found " << mutationTrees.size() << " mutation trees with " << size << " out of " << _F.getNrCharacters() << " mutations" << std::endl;
 }
 
 void EnumerateMutationTrees::getFrequencyTensor(RealTensor& F_lb,
@@ -142,3 +145,63 @@ void EnumerateMutationTrees::getFrequencyTensor(RealTensor& F_lb,
     }
   }
 }
+
+std::ostream& operator<<(std::ostream& out,
+                         const EnumerateMutationTrees::TreeVector& trees)
+{
+  out << trees.size() << " #trees" << std::endl;
+  int idx = 0;
+  for (const CloneTree& T : trees)
+  {
+    out << lemon::countArcs(T.tree()) << " # edges, tree " << idx + 1 << std::endl;
+    T.write(out);
+    ++idx;
+  }
+  return out;
+}
+
+std::istream& operator>>(std::istream& in,
+                         EnumerateMutationTrees::TreeVector& trees)
+{
+  int nrTrees = -1;
+  
+  std::string line;
+  getline(in, line);
+  
+  std::stringstream ss(line);
+  ss >> nrTrees;
+  
+  if (nrTrees < 0)
+  {
+    throw std::runtime_error("Error: number of trees should be nonnegative");
+  }
+  
+  for (int i = 0; i < nrTrees; ++i)
+  {
+    int nrEdges = -1;
+    
+    getline(in, line);
+    ss.clear();
+    ss.str(line);
+    ss >> nrEdges;
+
+    if (nrEdges < 0)
+    {
+      throw std::runtime_error("Error: number of edges should be nonnegative");
+    }
+    
+    ss.clear();
+    for (int j = 0; j < nrEdges; ++j)
+    {
+      getline(in, line);
+      ss << line;
+    }
+    
+    CloneTree T;
+    T.read(ss);
+    
+    trees.push_back(T);
+  }
+  
+  return in;
+}

src/enumeratemutationtrees.h

@@ -28,10 +28,12 @@ class EnumerateMutationTrees
   /// @param outputDirectory Output directory, may be empty
   /// @param nrThreads Number of threads to use in the enumeration
   /// @param limit Maximum number of trees to enumerate
+  /// @param timeLimit Time limit for the enumeration
   /// @param mutationTrees Vector to store enumerated mutation trees
   void enumerate(const std::string& outputDirectory,
                  const int nrThreads,
                  const int limit,
+                 const int timeLimit,
                  TreeVector& mutationTrees) const;
 
 private:
@@ -47,4 +49,7 @@ class EnumerateMutationTrees
   const FrequencyMatrix& _F;
 };
 
+std::ostream& operator<<(std::ostream& out, const EnumerateMutationTrees::TreeVector& trees);
+std::istream& operator>>(std::istream& in, EnumerateMutationTrees::TreeVector& trees);
+
 #endif // ENUMERATEMUTATIONTREES_H

src/generatemutationtreesmain.cpp

@@ -20,10 +20,12 @@ int main(int argc, char** argv)
   bool canonicalCloneTrees = false;
   int nrThreads = 1;
   int limit = -1;
+  int timeLimit = -1;
 
   lemon::ArgParser ap(argc, argv);
   ap.refOption("o", "Output directory" , outputDirectory);
-  ap.refOption("l", "Maxim number of mutation trees to enumerate (default: -1, unlimited)" , limit);
+  ap.refOption("l", "Maximum number of mutation trees to enumerate (default: -1, unlimited)" , limit);
+  ap.refOption("tl", "Time limit in seconds (default: -1, unlimited)" , timeLimit);
   ap.refOption("C", "Output canonical clone trees", canonicalCloneTrees);
   ap.refOption("t", "Number of threads (default: 1)", nrThreads);
   ap.other("frequencies", "Frequencies");
@@ -73,7 +75,16 @@ int main(int argc, char** argv)
   {
     EnumerateMutationTrees::TreeVector trees;
     EnumerateMutationTrees enumerate(F);
-    enumerate.enumerate(outputDirectory, nrThreads, limit, trees);
+    enumerate.enumerate(outputDirectory,
+                        nrThreads,
+                        limit,
+                        timeLimit,
+                        trees);
+    
+    if (outputDirectory.empty())
+    {
+      std::cout << trees;
+    }
   }
 
   return 0;

src/ilpbinarizationsolver.cpp

@@ -258,7 +258,7 @@ void IlpBinarizationSolver::run(const CloneTree& T,
                                 bool outputILP,
                                 bool outputSearchGraph,
                                 int timeLimit,
-                                double UB,
+                                const IntTriple& bounds,
                                 const StringPairList& forcedComigrations)
 {
   char buf[1024];
@@ -289,7 +289,7 @@ void IlpBinarizationSolver::run(const CloneTree& T,
                                pattern,
                                filenameGurobiLog,
                                forcedComigrations);
-  solver.init(UB);
+  solver.init(bounds);
 
   if (outputSearchGraph)
   {

src/ilpbinarizationsolver.h

@@ -77,7 +77,7 @@ class IlpBinarizationSolver : public IlpSolver
   /// @param outputILP Output ILP model
   /// @param outputSearchGraph Output search graph
   /// @param timeLimit Time limit in seconds
-  /// @param UB Upper bound on objective value
+  /// @param bounds Upper bounds on mu, gamma and sigma
   /// @param forcedComigrations List of ordered pairs of anatomical sites
   /// that must be present
   static void run(const CloneTree& T,
@@ -89,7 +89,7 @@ class IlpBinarizationSolver : public IlpSolver
                   bool outputILP,
                   bool outputSearchGraph,
                   int timeLimit,
-                  double UB,
+                  const IntTriple& bounds,
                   const StringPairList& forcedComigrations);
 
 protected:

src/ilpsolver.cpp

@@ -36,14 +36,14 @@ IlpSolver::IlpSolver(const CloneTree& T,
 {
 }
 
-void IlpSolver::init(double upperBound)
+void IlpSolver::init(const IntTriple& bounds)
 {
   initIndices();
   initVariables();
   initConstraints();
   initLeafConstraints();
   initForcedComigrations();
-  initObjective(upperBound);
+  initObjective(bounds);
   _model.update();
 }
 
@@ -454,12 +454,12 @@ void IlpSolver::initParallelSingleSourceSeedingConstraints()
   }
 }
 
-void IlpSolver::initObjective(double upperBound)
+void IlpSolver::initObjective(const IntTriple& bounds)
 {
   const int nrArcs = _indexToArc.size();
   const int nrAnatomicalSites = _indexToAnatomicalSite.size();
 
-  GRBLinExpr obj, sum_migrations, sum_comigrations;
+  GRBLinExpr obj, sum_migrations, sum_comigrations, sum_seeding_sites;
   // migrations
   for (int ij = 0; ij < nrArcs; ++ij)
   {
@@ -468,6 +468,11 @@ void IlpSolver::initObjective(double upperBound)
   }
   _model.addConstr(sum_migrations >= nrAnatomicalSites - 1);
 
+  if (bounds.first != -1)
+  {
+    _model.addConstr(sum_migrations <= bounds.first);
+  }
+  
   if (_forcedComigrations.empty()
       || _pattern == MigrationGraph::M
       || _pattern == MigrationGraph::R)
@@ -487,19 +492,25 @@ void IlpSolver::initObjective(double upperBound)
     _model.addConstr(sum_comigrations >= nrAnatomicalSites - 1);
     _model.addConstr(sum_migrations >= sum_comigrations);
 
+    if (bounds.second.first != -1)
+    {
+      _model.addConstr(sum_comigrations <= bounds.second.first);
+    }
+    
     // seeding sites
     {
       const double g = (1. / (nrArcs + 1)) * (1. / (nrAnatomicalSites + 1));
       for (int s = 0; s < nrAnatomicalSites; ++s)
       {
         obj += g * _d[s];
+        sum_seeding_sites += _d[s];
       }
     }
-  }
-  
-  if (upperBound >= 0)
-  {
-    _model.addConstr(obj <= upperBound);
+    
+    if (bounds.second.second != -1)
+    {
+      _model.addConstr(sum_seeding_sites <= bounds.second.second);
+    }
   }
 
   _model.setObjective(obj, GRB_MINIMIZE);
@@ -564,7 +575,7 @@ bool IlpSolver::solve(int nrThreads, int timeLimit)
     _model.getEnv().set(GRB_IntParam_LogToConsole, 0);
     _model.optimize();
     int status = _model.get(GRB_IntAttr_Status);
-    if (status == GRB_OPTIMAL || status == GRB_TIME_LIMIT)
+    if (status == GRB_OPTIMAL || status == GRB_SUBOPTIMAL)
     {
       _LB = _model.get(GRB_DoubleAttr_ObjBound);
       _UB = _model.get(GRB_DoubleAttr_ObjVal);
@@ -588,6 +599,17 @@ bool IlpSolver::solve(int nrThreads, int timeLimit)
       std::cerr << "Model is unbounded" << std::endl;
       return false;
     }
+    else if (status == GRB_TIME_LIMIT)
+    {
+      _LB = _model.get(GRB_DoubleAttr_ObjBound);
+      _UB = _model.get(GRB_DoubleAttr_ObjVal);
+      if (_UB <= lemon::countNodes(_T.tree()) * _T.getNrAnatomicalSites())
+      {
+        processSolution();
+        return true;
+      }
+      return false;
+    }
   }
   catch (const GRBException& e)
   {
@@ -613,7 +635,7 @@ void IlpSolver::run(const CloneTree& T,
                     bool outputILP,
                     bool outputSearchGraph,
                     int timeLimit,
-                    double UB,
+                    const IntTriple& bounds,
                     const StringPairList& forcedComigrations)
 {
   char buf[1024];
@@ -642,7 +664,7 @@ void IlpSolver::run(const CloneTree& T,
                    pattern,
                    filenameGurobiLog,
                    forcedComigrations);
-  solver.init(UB);
+  solver.init(bounds);
 
   if (!outputDirectory.empty() && outputILP)
   {

src/ilpsolver.h

@@ -59,7 +59,9 @@ class IlpSolver
   }
 
   /// Initialize ILP
-  virtual void init(double upperBound);
+  ///
+  /// @param bounds Upper bounds on mu, gamma and sigma
+  virtual void init(const IntTriple& bounds);
 
   /// Return lower bound
   double LB() const
@@ -84,8 +86,7 @@ class IlpSolver
   /// @param outputILP Output ILP model
   /// @param outputSearchGraph Output search graph
   /// @param timeLimit Time limit in seconds
-  ///
-  /// @param UB Upper bound on objective value
+  /// @param bounds Upper bounds on mu, gamma and sigma
   /// @param forcedComigrations List of ordered pairs of anatomical sites
   /// that must be present
   static void run(const CloneTree& T,
@@ -97,7 +98,7 @@ class IlpSolver
                   bool outputILP,
                   bool outputSearchGraph,
                   int timeLimit,
-                  double UB,
+                  const IntTriple& bounds,
                   const StringPairList& forcedComigrations);
 
 protected:
@@ -124,8 +125,8 @@ class IlpSolver
 
   /// Initialize ILP objective
   ///
-  /// @param upperBound Upper bound on objective function
-  virtual void initObjective(double upperBound);
+  /// @param bounds Upper bounds on mu, gamma and sigma
+  virtual void initObjective(const IntTriple& bounds);
 
   /// Initialize ILP constraints regarding primary tumor anatomical site
   /// (only applicable for PS and S)