reHC*, added option for searching a solution with at least a given nu…

…mber of recombinations.

docs/README.md

@@ -224,6 +224,14 @@ Here `XX` is a number between `0.0` and `1.0` that represents the
 maximum number of recombinations *r* as a fraction of the total number
 of possible recombination loci, while `YY` is (directly) the maximum
 number of recombinations *r*.
+Moreover, if option `--global-recomb` is enabled and
+`--global-recomb-number` is used, it is also possible to search for a
+haplotype configuration with a given minimum number of recombinations by
+specifying the option `--global-recomb-min-number=ZZ`, where `ZZ` is the
+sought lower bound.
+This option should only be used to specify a lower bound that has been
+already proved since the resulting haplotype configuration could induce
+unnecessary recombination in order to satisfy the given lower bound.
 
 Similarly, to enable genotyping errors in the computed haplotype
 configuration, the program options `--global-error` and

include/ped2cnf-constraints.hpp

@@ -162,6 +162,24 @@ class at_most_global_constraints_abs_t
 
 };
 
+class interval_global_constraints_abs_t
+  :public constraint_handler_t
+{
+private:
+  const unsigned int _vmin;
+  const unsigned int _vmax;
+
+  virtual void
+  _handle_constraints(pedcnf_t& cnf, const individuals_variables_t& variables) const;
+
+public:
+  interval_global_constraints_abs_t(const unsigned int vmin, const unsigned int vmax,
+												const std::string& description="true variables")
+		:constraint_handler_t(description), _vmin(vmin), _vmax(vmax)
+  {};
+
+};
+
 
 class at_most_windowed_constraints_t
   :public constraint_handler_t

include/pedcnf.hpp

@@ -364,6 +364,11 @@ add_card_constraint_less_or_equal_than(pedcnf_t& cnf,
 													const std::vector<var_t>& in_vars,
 													const size_t k);
 
+void
+add_card_constraint_between(pedcnf_t& cnf,
+									 const std::vector<var_t>& in_vars,
+									 const size_t k1, const size_t k2);
+
 void
 add_uniform_card_constraint_less_or_equal_than(pedcnf_t& cnf,
 															  const std::vector<var_t>& in_vars,

include/rehc_app.hpp

@@ -118,9 +118,17 @@ class rehcstar_t:
 // have been specified
 		if (vm.count("global-recomb-number") && !vm["global-recomb-number"].defaulted()) {
 		  const unsigned int recomb_number= vm["global-recomb-number"].as<unsigned int>();
-		  L_INFO("Enabling *GLOBAL* recombination handling ("
-					"recomb-number=" << recomb_number << ")");
-		  global_recomb_constraints.add(new at_most_global_constraints_abs_t(recomb_number, "recombinations"));
+		  if (vm.count("global-recomb-min-number") && !vm["global-recomb-min-number"].defaulted()) {
+			 const unsigned int recomb_min_number= vm["global-recomb-min-number"].as<unsigned int>();
+			 L_INFO("Enabling *GLOBAL* recombination handling ("
+					  "recomb-number=" << recomb_number << ", "
+					  "recomb-min-number=" << recomb_min_number << ")");
+			 global_recomb_constraints.add(new interval_global_constraints_abs_t(recomb_min_number, recomb_number, "recombinations"));
+		  } else {
+			 L_INFO("Enabling *GLOBAL* recombination handling ("
+					  "recomb-number=" << recomb_number << ")");
+			 global_recomb_constraints.add(new at_most_global_constraints_abs_t(recomb_number, "recombinations"));
+		  }
 		} else {
 		  const double recomb_rate= vm["global-recomb-rate"].as<double>();
 		  L_INFO("Enabling *GLOBAL* recombination handling ("

src/application/main.cpp

@@ -181,6 +181,13 @@ class rehcstar_application_t: public application_t {
 		("global-recomb-number", po::value< unsigned int >()->default_value(1),
 		 "Maximum number of recombinations in all the genotypes "
 		 "(used only if '--global-recomb' is specified, cannot be used with '--global-recomb-rate').")
+		("global-recomb-min-number", po::value< unsigned int >()->default_value(0),
+		 "Minimum number of recombinations in all the genotypes "
+		 "(used only if '--global-recomb' is specified, cannot be used with '--global-recomb-rate').\n"
+		 "NOTE: This lower bound is not strictly enforced, since the SAT solver could compute a "
+		 "solution with unnecessary recombinations. This option should be used to improve the "
+		 "efficiency of the search of the optimum and should be set to a value such that "
+		 "a solution with this number of recombinations does not exist.")
 		("individual-recomb", po::bool_switch()->default_value(false),
 		 "Enable INDIVIDUAL recombination handling (i.e., the recombination rate in each genotype is "
 		 "less than or equal to the specified recombination rate, computed over *ALL* loci).")
@@ -257,9 +264,19 @@ class rehcstar_application_t: public application_t {
 	 option_dependency(vm, "error-window-length", "uniform-error");
 	 option_dependency(vm, "global-recomb-rate", "global-recomb");
 	 option_dependency(vm, "global-recomb-number", "global-recomb");
+	 option_dependency(vm, "global-recomb-min-number", "global-recomb");
+	 conflicting_options(vm, "global-recomb-rate", "global-recomb-number");
+	 conflicting_options(vm, "global-recomb-rate", "global-recomb-min-number");
 	 option_dependency(vm, "individual-recomb-rate", "individual-recomb");
 	 option_dependency(vm, "max-recombs-in-window", "uniform-recomb");
 	 option_dependency(vm, "recomb-window-length", "uniform-recomb");
+	 if (vm["global-recomb"].as<bool>()) {
+		const unsigned int rmin= vm["global-recomb-min-number"].as<unsigned int>();
+		const unsigned int rmax= vm["global-recomb-number"].as<unsigned int>();
+		if (rmin > rmax) {
+		  throw std::logic_error(std::string("The minimum number of recombinations must be not greater than the maximum number of recombinations."));
+		}
+	 }
 	 if (vm["uniform-error"].as<bool>()) {
 		const unsigned int wlen= vm["error-window-length"].as<unsigned int>();
 		const unsigned int merr= vm["max-errors-in-window"].as<unsigned int>();

src/lib/ped2cnf-constraints.cpp

@@ -103,6 +103,18 @@ at_most_global_constraints_abs_t::_handle_constraints(pedcnf_t& cnf,
   add_card_constraint_less_or_equal_than(cnf, all_vars, _limit);
 };
 
+void
+interval_global_constraints_abs_t::_handle_constraints(pedcnf_t& cnf,
+																		 const individuals_variables_t& variables) const {
+  individual_variables_t all_vars;
+  BOOST_FOREACH(const individual_variables_t& ivar, variables) {
+	 all_vars.insert(all_vars.end(), ivar.begin(), ivar.end());
+  }
+  INFO("Globally, number of " << _description << " between " << _vmin <<
+		 " and " << _vmax << " (over " << all_vars.size() << ")");
+  add_card_constraint_between(cnf, all_vars, _vmin, _vmax);
+};
+
 void
 at_most_windowed_constraints_t::_handle_constraints(pedcnf_t& cnf,
 																	 const individuals_variables_t& variables) const {

src/lib/pedcnf_card.cpp

@@ -295,14 +295,14 @@ le_half_sorting_network(pedcnf_t& cnf,
 
 
 static void
-le_card_network(pedcnf_t& cnf,
-					 const std::vector<var_t>& in_vars,
-					 const size_t k) {
+le_card_network_int(pedcnf_t& cnf,
+						  const std::vector<var_t>& in_vars,
+						  const size_t k1, const size_t k2) {
   my_logger logger(get_my_logger("card_constraints"));
   DEBUG("Generating a cardinality network for encoding the cardinality constraint.");
-  const size_t p= pow2_of_ceiling_log2(k+1);
+  const size_t p= pow2_of_ceiling_log2(k2+1);
   DEBUG("Dividing " << in_vars.size() << " variables in blocks of size "
-		  << p << " since the upper bound is " << k << ".");
+		  << p << " since the upper bound is " << k2 << ".");
 // Ensure that the variables are multiples of p
   std::vector<var_t> all_vars(in_vars.begin(), in_vars.end());
   while (all_vars.size() % p != 0) {
@@ -321,7 +321,7 @@ le_card_network(pedcnf_t& cnf,
 	 input.insert(input.end(), next_in_var, next_in_var+p);
 	 std::vector<var_t> hsort_out;
 	 generate_hsort(cnf, input, hsort_out, logger);
-	 cnf.add_clause<1>((lit_t[]){ -hsort_out[k] });
+	 cnf.add_clause<1>((lit_t[]){ -hsort_out[k2] });
 	 // for (size_t limit= k; limit<p; ++limit) {
 	 // 	cnf.add_clause<1>((lit_t[]){ -hsort_out[limit] });
 	 // }
@@ -331,13 +331,25 @@ le_card_network(pedcnf_t& cnf,
 	 prev_ris.insert(prev_ris.end(), smerge_out.begin(), smerge_out.begin()+p);
 	 next_in_var += p;
 	 next_in_var_counter += p;
-	 cnf.add_clause<1>((lit_t[]){ -smerge_out[k] });
+	 cnf.add_clause<1>((lit_t[]){ -smerge_out[k2] });
 	 // for (size_t limit= k; limit<=p; ++limit) {
 	 // 	cnf.add_clause<1>((lit_t[]){ -smerge_out[limit] });
 	 // }
   }
+  if (k1 > 0) {
+	 DEBUG("Setting the lower bound >= " << k1);
+	 for (size_t i= 0; i<k1; ++i)
+		cnf.add_clause<1>((lit_t[]){ prev_ris[i] });
+  }
 };
 
+static void
+le_card_network(pedcnf_t& cnf,
+					 const std::vector<var_t>& in_vars,
+					 const size_t k) {
+  le_card_network_int(cnf, in_vars, 0, k);
+}
+
 static void
 le_card_network_tree_rec(pedcnf_t& cnf,
 								 std::vector<var_t>& out_vars,
@@ -469,6 +481,25 @@ add_card_constraint_less_or_equal_than(pedcnf_t& cnf,
   }
 };
 
+void
+add_card_constraint_between(pedcnf_t& cnf,
+									 const std::vector<var_t>& in_vars,
+									 const size_t k1, const size_t k2) {
+  my_logger logger(get_my_logger("card_constraints"));
+  if (k2==0) {
+	 eq_zero(cnf, in_vars);
+  } else if (k2 >= in_vars.size()) {
+	 INFO("The number of variables (" << in_vars.size() << ") "
+			"is not greater than the numeric upper bound (" << k2 << "). "
+			"Setting to " << in_vars.size() << ".");
+	 le_card_network_int(cnf, in_vars,
+								std::min(k1, in_vars.size()), in_vars.size());
+  } else {
+// Use Cardinality Networks
+	 le_card_network_int(cnf, in_vars, k1, k2);
+  }
+};
+
 
 static void
 add_uniform_card_constraint_less_or_equal_than_int(pedcnf_t& cnf,