more work on SAT

makefiles/Makefile.gen.mk

@@ -1108,6 +1108,7 @@ SAT_DEPS = \
  $(SRC_DIR)/ortools/sat/drat_proof_handler.h \
  $(SRC_DIR)/ortools/sat/drat_writer.h \
  $(SRC_DIR)/ortools/sat/encoding.h \
+ $(SRC_DIR)/ortools/sat/feasibility_pump.h \
  $(SRC_DIR)/ortools/sat/implied_bounds.h \
  $(SRC_DIR)/ortools/sat/integer_expr.h \
  $(SRC_DIR)/ortools/sat/integer.h \
@@ -1173,6 +1174,7 @@ SAT_LIB_OBJS = \
  $(OBJ_DIR)/sat/drat_proof_handler.$O \
  $(OBJ_DIR)/sat/drat_writer.$O \
  $(OBJ_DIR)/sat/encoding.$O \
+ $(OBJ_DIR)/sat/feasibility_pump.$O \
  $(OBJ_DIR)/sat/implied_bounds.$O \
  $(OBJ_DIR)/sat/integer.$O \
  $(OBJ_DIR)/sat/integer_expr.$O \
@@ -1716,6 +1718,23 @@ objs/sat/encoding.$O: ortools/sat/encoding.cc ortools/sat/encoding.h \
  ortools/util/integer_pq.h | $(OBJ_DIR)/sat
 	$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Ssat$Sencoding.cc $(OBJ_OUT)$(OBJ_DIR)$Ssat$Sencoding.$O
 
+objs/sat/feasibility_pump.$O: ortools/sat/feasibility_pump.cc ortools/sat/feasibility_pump.h \
+ ortools/base/int_type.h ortools/base/macros.h \
+ ortools/base/integral_types.h ortools/base/logging.h \
+ ortools/gen/ortools/sat/boolean_problem.pb.h ortools/sat/pb_constraint.h \
+ ortools/base/int_type_indexed_vector.h ortools/sat/model.h \
+ ortools/base/map_util.h ortools/base/typeid.h ortools/sat/sat_base.h \
+ ortools/util/bitset.h ortools/gen/ortools/sat/sat_parameters.pb.h \
+ ortools/util/stats.h ortools/base/timer.h ortools/base/basictypes.h \
+ ortools/sat/sat_solver.h ortools/base/hash.h ortools/sat/clause.h \
+ ortools/sat/drat_proof_handler.h ortools/sat/drat_checker.h \
+ ortools/sat/drat_writer.h ortools/base/file.h \
+ ortools/util/random_engine.h ortools/util/time_limit.h \
+ ortools/base/commandlineflags.h ortools/util/running_stat.h \
+ ortools/sat/restart.h ortools/sat/sat_decision.h \
+ ortools/util/integer_pq.h | $(OBJ_DIR)/sat
+	$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sortools$Ssat$Sfeasibility_pump.cc $(OBJ_OUT)$(OBJ_DIR)$Ssat$Sfeasibility_pump.$O
+
 objs/sat/implied_bounds.$O: ortools/sat/implied_bounds.cc \
  ortools/sat/implied_bounds.h ortools/base/int_type.h \
  ortools/base/macros.h ortools/base/int_type_indexed_vector.h \
@@ -4100,4 +4119,3 @@ $(GEN_DIR)/ortools/constraint_solver/solver_parameters.pb.h: \
 $(OBJ_DIR)/constraint_solver/solver_parameters.pb.$O: \
  $(GEN_DIR)/ortools/constraint_solver/solver_parameters.pb.cc | $(OBJ_DIR)/constraint_solver
 	$(CCC) $(CFLAGS) -c $(GEN_PATH)$Sortools$Sconstraint_solver$Ssolver_parameters.pb.cc $(OBJ_OUT)$(OBJ_DIR)$Sconstraint_solver$Ssolver_parameters.pb.$O
-

ortools/sat/BUILD

@@ -170,6 +170,7 @@ cc_library(
         ":cp_model_utils",
         ":drat_checker",
         ":drat_proof_handler",
+        ":feasibility_pump",
         ":integer",
         ":integer_expr",
         ":integer_search",
@@ -1132,6 +1133,27 @@ cc_library(
     ],
 )
 
+cc_library(
+    name = "feasibility_pump",
+    srcs = ["feasibility_pump.cc"],
+    hdrs = ["feasibility_pump.h"],
+    visibility = ["//visibility:public"],
+    deps = [
+        ":cp_model_loader",
+        ":cp_model_cc_proto",
+        ":integer",
+        ":linear_constraint",
+        ":synchronization",
+        ":util",
+        "//ortools/base",
+        "//ortools/glop:revised_simplex",
+        "//ortools/lp_data",
+        "//ortools/lp_data:base",
+        "//ortools/lp_data:lp_data_utils",
+        "//ortools/util:saturated_arithmetic",
+    ],
+)
+
 cc_library(
     name = "rins",
     srcs = ["rins.cc"],

ortools/sat/cp_model_lns.cc

@@ -13,6 +13,7 @@
 
 #include "ortools/sat/cp_model_lns.h"
 
+#include <limits>
 #include <numeric>
 #include <random>
 #include <vector>
@@ -586,6 +587,10 @@ Neighborhood SchedulingTimeWindowNeighborhoodGenerator::Generate(
 }
 
 bool RelaxationInducedNeighborhoodGenerator::ReadyToGenerate() const {
+  if (incomplete_solutions_ != nullptr) {
+    return incomplete_solutions_->HasNewSolution();
+  }
+
   if (response_manager_ != nullptr) {
     if (response_manager_->SolutionsRepository().NumSolutions() == 0) {
       return false;
@@ -618,6 +623,10 @@ Neighborhood RelaxationInducedNeighborhoodGenerator::Generate(
       (relaxation_solutions_ != nullptr &&
        relaxation_solutions_->NumSolutions() > 0);
 
+  const bool incomplete_solution_available =
+      (incomplete_solutions_ != nullptr &&
+       incomplete_solutions_->HasNewSolution());
+
   if (!lp_solution_available && !relaxation_solution_available) {
     return neighborhood;
   }
@@ -632,14 +641,19 @@ Neighborhood RelaxationInducedNeighborhoodGenerator::Generate(
           ? random_bool(*random)
           : lp_solution_available;
   if (use_lp_relaxation) {
-    rins_neighborhood = GetRINSNeighborhood(response_manager_,
-                                            /*relaxation_solutions=*/nullptr,
-                                            lp_solutions_, random);
+    rins_neighborhood =
+        GetRINSNeighborhood(response_manager_,
+                            /*relaxation_solutions=*/nullptr, lp_solutions_,
+                            incomplete_solutions_, random);
+    neighborhood.source_info =
+        incomplete_solution_available ? "incomplete" : "lp";
   } else {
     CHECK(relaxation_solution_available);
-    rins_neighborhood =
-        GetRINSNeighborhood(response_manager_, relaxation_solutions_,
-                            /*lp_solutions=*/nullptr, random);
+    rins_neighborhood = GetRINSNeighborhood(
+        response_manager_, relaxation_solutions_,
+        /*lp_solutions=*/nullptr, incomplete_solutions_, random);
+    neighborhood.source_info =
+        incomplete_solution_available ? "incomplete" : "relaxation";
   }
 
   if (rins_neighborhood.fixed_vars.empty() &&

ortools/sat/cp_model_lns.h

@@ -49,6 +49,10 @@ struct Neighborhood {
   // TODO(user): Make sure that the id is unique for each generated
   // neighborhood for each generator.
   int64 id = 0;
+
+  // Used for identifying the source of the neighborhood if it is generated
+  // using solution repositories.
+  std::string source_info = "";
 };
 
 // Contains pre-computed information about a given CpModelProto that is meant
@@ -404,27 +408,36 @@ class SchedulingTimeWindowNeighborhoodGenerator : public NeighborhoodGenerator {
                         double difficulty, random_engine_t* random) final;
 };
 
-// Generates a neighborhood by fixing the variables who have same solution value
-// as their linear relaxation. This was published in "Exploring relaxation
-// induced neighborhoods to improve MIP solutions" 2004 by E. Danna et.
+// Generates a neighborhood by fixing the variables to solutions reported in
+// various repositories. This is inspired from RINS published in "Exploring
+// relaxation induced neighborhoods to improve MIP solutions" 2004 by E. Danna
+// et.
 //
-// If use_only_relaxation_values is true, this generates a neighborhood using
-// only the linear/general relaxation values. The domain of the variables are
-// reduced to the integer values around their lp solution/relaxation solution
-// values. This was published in "RENS – The Relaxation Enforced Neighborhood"
-// 2009 by Timo Berthold.
+// If incomplete_solutions is provided, this generates a neighborhood by fixing
+// the variable values to a solution in the SharedIncompleteSolutionManager and
+// ignores the other repositories.
+//
+// Otherwise, if response_manager is not provided, this generates a neighborhood
+// using only the linear/general relaxation values. The domain of the variables
+// are reduced to the integer values around their lp solution/relaxation
+// solution values. This was published in "RENS – The Relaxation Enforced
+// Neighborhood" 2009 by Timo Berthold.
 class RelaxationInducedNeighborhoodGenerator : public NeighborhoodGenerator {
  public:
   explicit RelaxationInducedNeighborhoodGenerator(
       NeighborhoodGeneratorHelper const* helper,
       const SharedResponseManager* response_manager,
       const SharedRelaxationSolutionRepository* relaxation_solutions,
-      const SharedLPSolutionRepository* lp_solutions, const std::string& name)
+      const SharedLPSolutionRepository* lp_solutions,
+      SharedIncompleteSolutionManager* incomplete_solutions,
+      const std::string& name)
       : NeighborhoodGenerator(name, helper),
         response_manager_(response_manager),
         relaxation_solutions_(relaxation_solutions),
-        lp_solutions_(lp_solutions) {
-    CHECK(lp_solutions_ != nullptr || relaxation_solutions_ != nullptr);
+        lp_solutions_(lp_solutions),
+        incomplete_solutions_(incomplete_solutions) {
+    CHECK(lp_solutions_ != nullptr || relaxation_solutions_ != nullptr ||
+          incomplete_solutions != nullptr);
   }
 
   // Both initial solution and difficulty values are ignored.
@@ -434,9 +447,11 @@ class RelaxationInducedNeighborhoodGenerator : public NeighborhoodGenerator {
   // Returns true if the required solutions are available.
   bool ReadyToGenerate() const override;
 
+ private:
   const SharedResponseManager* response_manager_;
   const SharedRelaxationSolutionRepository* relaxation_solutions_;
   const SharedLPSolutionRepository* lp_solutions_;
+  SharedIncompleteSolutionManager* incomplete_solutions_;
 };
 
 // Generates a relaxation of the original model by removing a consecutive span

ortools/sat/cp_model_solver.cc

@@ -25,6 +25,7 @@
 #include <vector>
 
 #include "ortools/base/cleanup.h"
+#include "ortools/sat/feasibility_pump.h"
 
 #if !defined(__PORTABLE_PLATFORM__)
 #include "absl/synchronization/notification.h"
@@ -551,9 +552,8 @@ void TryToAddCutGenerators(const CpModelProto& model_proto,
 
 }  // namespace
 
-// Adds one LinearProgrammingConstraint per connected component of the model.
-IntegerVariable AddLPConstraints(const CpModelProto& model_proto,
-                                 int linearization_level, Model* m) {
+LinearRelaxation ComputeLinearRelaxation(const CpModelProto& model_proto,
+                                         int linearization_level, Model* m) {
   LinearRelaxation relaxation;
 
   // Linearize the constraints.
@@ -670,6 +670,14 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto,
   VLOG(3) << relaxation.linear_constraints.size()
           << " constraints in the LP relaxation.";
   VLOG(3) << relaxation.cut_generators.size() << " cuts generators.";
+  return relaxation;
+}
+
+// Adds one LinearProgrammingConstraint per connected component of the model.
+IntegerVariable AddLPConstraints(const CpModelProto& model_proto,
+                                 int linearization_level, Model* m) {
+  const LinearRelaxation relaxation =
+      ComputeLinearRelaxation(model_proto, linearization_level, m);
 
   // The bipartite graph of LP constraints might be disconnected:
   // make a partition of the variables into connected components.
@@ -719,6 +727,7 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto,
   // as much as possible the objective bound by using any bounds the LP give
   // us on one of its components. This is critical on the zephyrus problems for
   // instance.
+  auto* mapping = m->GetOrCreate<CpModelMapping>();
   for (int i = 0; i < model_proto.objective().coeffs_size(); ++i) {
     const IntegerVariable var =
         mapping->Integer(model_proto.objective().vars(i));
@@ -1044,7 +1053,7 @@ void RegisterObjectiveBestBoundExport(
 
 // Registers a callback to import new objective bounds. It will be called each
 // time the search main loop is back to level zero. Note that it the presence of
-// assumptions, this will not happend until the set of assumptions is changed.
+// assumptions, this will not happen until the set of assumptions is changed.
 void RegisterObjectiveBoundsImport(
     SharedResponseManager* shared_response_manager, Model* model) {
   auto* solver = model->GetOrCreate<SatSolver>();
@@ -1098,6 +1107,55 @@ void RegisterObjectiveBoundsImport(
       import_objective_bounds);
 }
 
+void LoadFeasibilityPump(const CpModelProto& model_proto,
+                         SharedResponseManager* shared_response_manager,
+                         Model* model) {
+  CHECK(shared_response_manager != nullptr);
+
+  auto* sat_solver = model->GetOrCreate<SatSolver>();
+
+  // Simple function for the few places where we do "return unsat()".
+  const auto unsat = [shared_response_manager, sat_solver, model] {
+    sat_solver->NotifyThatModelIsUnsat();
+    shared_response_manager->NotifyThatImprovingProblemIsInfeasible(
+        absl::StrCat(model->GetOrCreate<WorkerInfo>()->worker_name,
+                     " [loading]"));
+  };
+
+  auto* mapping = model->GetOrCreate<CpModelMapping>();
+  const SatParameters& parameters = *(model->GetOrCreate<SatParameters>());
+  const bool view_all_booleans_as_integers =
+      (parameters.linearization_level() >= 2) ||
+      (parameters.search_branching() == SatParameters::FIXED_SEARCH &&
+       model_proto.search_strategy().empty());
+  mapping->CreateVariables(model_proto, view_all_booleans_as_integers, model);
+
+  // Check the model is still feasible before continuing.
+  if (sat_solver->IsModelUnsat()) return unsat();
+
+  // We don't load any constraint here.
+
+  if (parameters.linearization_level() == 0) return;
+
+  // Add linear constraints to Feasibility Pump.
+  const LinearRelaxation relaxation = ComputeLinearRelaxation(
+      model_proto, parameters.linearization_level(), model);
+  const int num_lp_constraints = relaxation.linear_constraints.size();
+  auto* feasibility_pump = model->GetOrCreate<FeasibilityPump>();
+  for (int i = 0; i < num_lp_constraints; i++) {
+    feasibility_pump->AddLinearConstraint(relaxation.linear_constraints[i]);
+  }
+
+  if (model_proto.has_objective()) {
+    for (int i = 0; i < model_proto.objective().coeffs_size(); ++i) {
+      const IntegerVariable var =
+          mapping->Integer(model_proto.objective().vars(i));
+      const int64 coeff = model_proto.objective().coeffs(i);
+      feasibility_pump->SetObjectiveCoefficient(var, IntegerValue(coeff));
+    }
+  }
+}
+
 // Loads a CpModelProto inside the given model.
 // This should only be called once on a given 'Model' class.
 //
@@ -1836,6 +1894,7 @@ struct SharedClasses {
   SharedResponseManager* response;
   SharedRelaxationSolutionRepository* relaxation_solutions;
   SharedLPSolutionRepository* lp_solutions;
+  SharedIncompleteSolutionManager* incomplete_solutions;
 
   bool SearchIsDone() {
     if (response->ProblemIsSolved()) return true;
@@ -1877,6 +1936,11 @@ class FullProblemSolver : public SubSolver {
       local_model_->Register<SharedLPSolutionRepository>(shared->lp_solutions);
     }
 
+    if (shared->incomplete_solutions != nullptr) {
+      local_model_->Register<SharedIncompleteSolutionManager>(
+          shared->incomplete_solutions);
+    }
+
     // Level zero variable bounds sharing.
     if (shared_->bounds != nullptr) {
       RegisterVariableBoundsLevelZeroExport(
@@ -1982,6 +2046,118 @@ class FullProblemSolver : public SubSolver {
   bool previous_task_is_completed_ ABSL_GUARDED_BY(mutex_) = true;
 };
 
+class FeasibilityPumpSolver : public SubSolver {
+ public:
+  FeasibilityPumpSolver(int id, const SatParameters& local_parameters,
+                        SharedClasses* shared)
+      : SubSolver(id, "feasibility_pump"),
+        shared_(shared),
+        local_model_(absl::make_unique<Model>()) {
+    // Setup the local model parameters and time limit.
+    local_model_->Add(NewSatParameters(local_parameters));
+    shared_->time_limit->UpdateLocalLimit(
+        local_model_->GetOrCreate<TimeLimit>());
+
+    // Stores info that will be used for logs in the local model.
+    WorkerInfo* worker_info = local_model_->GetOrCreate<WorkerInfo>();
+    worker_info->worker_name = "feasibility_pump";
+    worker_info->worker_id = id;
+
+    if (shared->response != nullptr) {
+      local_model_->Register<SharedResponseManager>(shared->response);
+    }
+
+    if (shared->relaxation_solutions != nullptr) {
+      local_model_->Register<SharedRelaxationSolutionRepository>(
+          shared->relaxation_solutions);
+    }
+
+    if (shared->lp_solutions != nullptr) {
+      local_model_->Register<SharedLPSolutionRepository>(shared->lp_solutions);
+    }
+
+    if (shared->incomplete_solutions != nullptr) {
+      local_model_->Register<SharedIncompleteSolutionManager>(
+          shared->incomplete_solutions);
+    }
+
+    // Level zero variable bounds sharing.
+    if (shared_->bounds != nullptr) {
+      RegisterVariableBoundsLevelZeroImport(
+          *shared_->model_proto, shared_->bounds, local_model_.get());
+    }
+  }
+
+  bool TaskIsAvailable() override {
+    if (shared_->SearchIsDone()) return false;
+    absl::MutexLock mutex_lock(&mutex_);
+    return previous_task_is_completed_;
+  }
+
+  std::function<void()> GenerateTask(int64 task_id) override {
+    return [this]() {
+      {
+        absl::MutexLock mutex_lock(&mutex_);
+        if (!previous_task_is_completed_) return;
+        previous_task_is_completed_ = false;
+      }
+      {
+        absl::MutexLock mutex_lock(&mutex_);
+        if (solving_first_chunk_) {
+          LoadFeasibilityPump(*shared_->model_proto, shared_->response,
+                              local_model_.get());
+          solving_first_chunk_ = false;
+          // Abort first chunk and allow to schedule the next.
+          previous_task_is_completed_ = true;
+          return;
+        }
+      }
+
+      auto* time_limit = local_model_->GetOrCreate<TimeLimit>();
+      const double saved_dtime = time_limit->GetElapsedDeterministicTime();
+      auto* feasibility_pump = local_model_->Mutable<FeasibilityPump>();
+      feasibility_pump->Solve();
+
+      {
+        absl::MutexLock mutex_lock(&mutex_);
+        deterministic_time_since_last_synchronize_ +=
+            time_limit->GetElapsedDeterministicTime() - saved_dtime;
+      }
+
+      // Abort if the problem is solved.
+      if (shared_->SearchIsDone()) {
+        shared_->time_limit->Stop();
+        return;
+      }
+
+      absl::MutexLock mutex_lock(&mutex_);
+      previous_task_is_completed_ = true;
+    };
+  }
+
+  void Synchronize() override {
+    absl::MutexLock mutex_lock(&mutex_);
+    deterministic_time_ += deterministic_time_since_last_synchronize_;
+    shared_->time_limit->AdvanceDeterministicTime(
+        deterministic_time_since_last_synchronize_);
+    deterministic_time_since_last_synchronize_ = 0.0;
+  }
+
+ private:
+  SharedClasses* shared_;
+  std::unique_ptr<Model> local_model_;
+
+  absl::Mutex mutex_;
+
+  // The first chunk is special. It is the one in which we load the linear
+  // constraints.
+  bool solving_first_chunk_ ABSL_GUARDED_BY(mutex_) = true;
+
+  double deterministic_time_since_last_synchronize_ ABSL_GUARDED_BY(mutex_) =
+      0.0;
+  bool previous_task_is_completed_ ABSL_GUARDED_BY(mutex_) = true;
+};
+
 // A Subsolver that generate LNS solve from a given neighborhood.
 class LnsSolver : public SubSolver {
  public:
@@ -2058,9 +2234,13 @@ class LnsSolver : public SubSolver {
       const double fully_solved_proportion =
           static_cast<double>(generator_->num_fully_solved_calls()) /
           std::max(int64{1}, generator_->num_calls());
+      std::string source_info = name();
+      if (!neighborhood.source_info.empty()) {
+        absl::StrAppend(&source_info, "_", neighborhood.source_info);
+      }
       const std::string solution_info = absl::StrFormat(
-          "%s(d=%0.2f s=%i t=%0.2f p=%0.2f)", name(), data.difficulty, task_id,
-          data.deterministic_limit, fully_solved_proportion);
+          "%s(d=%0.2f s=%i t=%0.2f p=%0.2f)", source_info, data.difficulty,
+          task_id, data.deterministic_limit, fully_solved_proportion);
 
       SatParameters local_params(parameters_);
       local_params.set_max_deterministic_time(data.deterministic_limit);
@@ -2283,6 +2463,14 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
       /*num_solutions_to_keep=*/10);
   global_model->Register<SharedLPSolutionRepository>(shared_lp_solutions.get());
 
+  std::unique_ptr<SharedIncompleteSolutionManager> shared_incomplete_solutions;
+  if (global_model->GetOrCreate<SatParameters>()->use_feasibility_pump()) {
+    shared_incomplete_solutions =
+        absl::make_unique<SharedIncompleteSolutionManager>();
+    global_model->Register<SharedIncompleteSolutionManager>(
+        shared_incomplete_solutions.get());
+  }
+
   SharedClasses shared;
   shared.model_proto = &model_proto;
   shared.wall_timer = wall_timer;
@@ -2291,6 +2479,7 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
   shared.response = shared_response_manager;
   shared.relaxation_solutions = shared_relaxation_solutions.get();
   shared.lp_solutions = shared_lp_solutions.get();
+  shared.incomplete_solutions = shared_incomplete_solutions.get();
 
   // The list of all the SubSolver that will be used in this parallel search.
   std::vector<std::unique_ptr<SubSolver>> subsolvers;
@@ -2353,6 +2542,12 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
     }
   }
 
+  // Add FeasibilityPumpSolver if enabled.
+  if (parameters.use_feasibility_pump() && !parameters.interleave_search()) {
+    subsolvers.push_back(absl::make_unique<FeasibilityPumpSolver>(
+        /*id=*/subsolvers.size(), parameters, &shared));
+  }
+
   // Add LNS SubSolver(s).
 
   // Add the NeighborhoodGeneratorHelper as a special subsolver so that its
@@ -2418,7 +2613,8 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
           /*id=*/subsolvers.size(),
           absl::make_unique<RelaxationInducedNeighborhoodGenerator>(
               helper, shared.response, shared.relaxation_solutions,
-              shared.lp_solutions, absl::StrCat("rins_lns_", strategy_name)),
+              shared.lp_solutions, /*incomplete_solutions=*/nullptr,
+              absl::StrCat("rins_lns_", strategy_name)),
           local_params, helper, &shared));
 
       // RENS.
@@ -2426,9 +2622,11 @@ void SolveCpModelParallel(const CpModelProto& model_proto,
           /*id=*/subsolvers.size(),
           absl::make_unique<RelaxationInducedNeighborhoodGenerator>(
               helper, /*respons_manager=*/nullptr, shared.relaxation_solutions,
-              shared.lp_solutions, absl::StrCat("rens_lns_", strategy_name)),
+              shared.lp_solutions, shared.incomplete_solutions,
+              absl::StrCat("rens_lns_", strategy_name)),
           local_params, helper, &shared));
     }
+
     if (parameters.use_relaxation_lns()) {
       subsolvers.push_back(absl::make_unique<LnsSolver>(
           /*id=*/subsolvers.size(),

ortools/sat/feasibility_pump.cc

@@ -0,0 +1,441 @@
+// Copyright 2010-2018 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#include "ortools/sat/feasibility_pump.h"
+
+#include <limits>
+
+#include "ortools/base/integral_types.h"
+#include "ortools/lp_data/lp_types.h"
+#include "ortools/sat/cp_model.pb.h"
+#include "ortools/util/saturated_arithmetic.h"
+
+namespace operations_research {
+namespace sat {
+
+using glop::ColIndex;
+using glop::Fractional;
+using glop::RowIndex;
+
+const double FeasibilityPump::kCpEpsilon = 1e-4;
+
+FeasibilityPump::FeasibilityPump(Model* model)
+    : sat_parameters_(*(model->GetOrCreate<SatParameters>())),
+      time_limit_(model->GetOrCreate<TimeLimit>()),
+      integer_trail_(model->GetOrCreate<IntegerTrail>()),
+      trail_(model->GetOrCreate<Trail>()),
+      integer_encoder_(model->GetOrCreate<IntegerEncoder>()),
+      incomplete_solutions_(model->Mutable<SharedIncompleteSolutionManager>()),
+      mapping_(model->Get<CpModelMapping>()) {
+  // Tweak the default parameters to make the solve incremental.
+  glop::GlopParameters parameters;
+  // TODO(user): Determine which simplex is better for this, dual or primal.
+  parameters.set_use_dual_simplex(false);
+  parameters.set_max_number_of_iterations(2000);
+  simplex_.SetParameters(parameters);
+  lp_data_.Clear();
+  integer_lp_.clear();
+}
+
+FeasibilityPump::~FeasibilityPump() {
+  VLOG(1) << "Feasibility Pump Total number of simplex iterations: "
+          << total_num_simplex_iterations_;
+}
+
+void FeasibilityPump::AddLinearConstraint(const LinearConstraint& ct) {
+  // We still create the mirror variable right away though.
+  for (const IntegerVariable var : ct.vars) {
+    GetOrCreateMirrorVariable(PositiveVariable(var));
+  }
+
+  integer_lp_.push_back(LinearConstraintInternal());
+  LinearConstraintInternal& new_ct = integer_lp_.back();
+  new_ct.lb = ct.lb;
+  new_ct.ub = ct.ub;
+  const int size = ct.vars.size();
+  CHECK_LE(ct.lb, ct.ub);
+  for (int i = 0; i < size; ++i) {
+    // We only use positive variable inside this class.
+    IntegerVariable var = ct.vars[i];
+    IntegerValue coeff = ct.coeffs[i];
+    if (!VariableIsPositive(var)) {
+      var = NegationOf(var);
+      coeff = -coeff;
+    }
+    new_ct.terms.push_back({GetOrCreateMirrorVariable(var), coeff});
+  }
+  // Important to keep lp_data_ "clean".
+  std::sort(new_ct.terms.begin(), new_ct.terms.end());
+}
+
+void FeasibilityPump::SetObjectiveCoefficient(IntegerVariable ivar,
+                                              IntegerValue coeff) {
+  objective_is_defined_ = true;
+  const IntegerVariable pos_var =
+      VariableIsPositive(ivar) ? ivar : NegationOf(ivar);
+  if (ivar != pos_var) coeff = -coeff;
+
+  const auto it = mirror_lp_variable_.find(pos_var);
+  if (it == mirror_lp_variable_.end()) return;
+  const glop::ColIndex col = it->second;
+  integer_objective_.push_back({col, coeff});
+  objective_infinity_norm_ =
+      std::max(objective_infinity_norm_, IntTypeAbs(coeff));
+}
+
+glop::ColIndex FeasibilityPump::GetOrCreateMirrorVariable(
+    IntegerVariable positive_variable) {
+  DCHECK(VariableIsPositive(positive_variable));
+
+  const auto it = mirror_lp_variable_.find(positive_variable);
+  if (it == mirror_lp_variable_.end()) {
+    const int model_var =
+        mapping_->GetProtoVariableFromIntegerVariable(positive_variable);
+    model_vars_size_ = std::max(model_vars_size_, model_var + 1);
+
+    const glop::ColIndex col(integer_variables_.size());
+    mirror_lp_variable_[positive_variable] = col;
+    integer_variables_.push_back(positive_variable);
+    lp_solution_.push_back(std::numeric_limits<double>::infinity());
+    integer_solution_.push_back(kMaxIntegerValue.value());
+
+    return col;
+  }
+  return it->second;
+}
+
+void FeasibilityPump::PrintStats() {
+  if (lp_solution_is_set_) {
+    VLOG(2) << "Fractionality: " << lp_solution_fractionality_;
+  } else {
+    VLOG(2) << "Fractionality: NA";
+    VLOG(2) << "simplex status: " << simplex_.GetProblemStatus();
+  }
+
+  if (integer_solution_is_set_) {
+    VLOG(2) << "#Infeasible const: " << num_infeasible_constraints_;
+    VLOG(2) << "Infeasibility: " << integer_solution_infeasibility_;
+  } else {
+    VLOG(2) << "Infeasibility: NA";
+  }
+}
+
+void FeasibilityPump::Solve() {
+  if (lp_data_.num_variables() == 0) {
+    InitializeWorkingLP();
+  }
+  UpdateBoundsOfLpVariables();
+  lp_solution_is_set_ = false;
+  integer_solution_is_set_ = false;
+
+  // Restore the original objective
+  for (ColIndex col(0); col < lp_data_.num_variables(); ++col) {
+    lp_data_.SetObjectiveCoefficient(col, 0.0);
+  }
+  for (const auto& term : integer_objective_) {
+    lp_data_.SetObjectiveCoefficient(term.first, ToDouble(term.second));
+  }
+
+  // TODO(user): Add cycle detection.
+  mixing_factor_ = 1.0;
+  for (int i = 0; i < max_fp_iterations_; ++i) {
+    L1DistanceMinimize();
+    if (!SolveLp()) break;
+    if (lp_solution_is_integer_) break;
+    SimpleRounding();
+    // We don't end this loop if the integer solutions is feasible in hope to
+    // get better solution.
+    if (integer_solution_is_feasible_) MaybePushToRepo();
+  }
+
+  PrintStats();
+  MaybePushToRepo();
+}
+
+void FeasibilityPump::MaybePushToRepo() {
+  if (incomplete_solutions_ == nullptr) return;
+
+  std::vector<double> lp_solution(model_vars_size_,
+                                  std::numeric_limits<double>::infinity());
+  // TODO(user): Consider adding solutions that have low fractionality.
+  if (lp_solution_is_integer_) {
+    // Fill the solution using LP solution values.
+    for (const IntegerVariable positive_var : integer_variables_) {
+      const int model_var =
+          mapping_->GetProtoVariableFromIntegerVariable(positive_var);
+      if (model_var >= 0 && model_var < model_vars_size_) {
+        lp_solution[model_var] = GetLPSolutionValue(positive_var);
+      }
+    }
+    incomplete_solutions_->AddNewSolution(lp_solution);
+  }
+
+  if (integer_solution_is_feasible_) {
+    // Fill the solution using Integer solution values.
+    for (const IntegerVariable positive_var : integer_variables_) {
+      const int model_var =
+          mapping_->GetProtoVariableFromIntegerVariable(positive_var);
+      if (model_var >= 0 && model_var < model_vars_size_) {
+        lp_solution[model_var] = GetIntegerSolutionValue(positive_var);
+      }
+    }
+    incomplete_solutions_->AddNewSolution(lp_solution);
+  }
+}
+
+// ----------------------------------------------------------------
+// -------------------LPSolving------------------------------------
+// ----------------------------------------------------------------
+
+void FeasibilityPump::InitializeWorkingLP() {
+  lp_data_.Clear();
+  // Create variables.
+  for (int i = 0; i < integer_variables_.size(); ++i) {
+    CHECK_EQ(glop::ColIndex(i), lp_data_.CreateNewVariable());
+    lp_data_.SetVariableType(glop::ColIndex(i),
+                             glop::LinearProgram::VariableType::INTEGER);
+  }
+
+  // Add constraints.
+  for (const LinearConstraintInternal& ct : integer_lp_) {
+    const ConstraintIndex row = lp_data_.CreateNewConstraint();
+    lp_data_.SetConstraintBounds(row, ToDouble(ct.lb), ToDouble(ct.ub));
+    for (const auto& term : ct.terms) {
+      lp_data_.SetCoefficient(row, term.first, ToDouble(term.second));
+    }
+  }
+
+  // Add objective.
+  for (const auto& term : integer_objective_) {
+    lp_data_.SetObjectiveCoefficient(term.first, ToDouble(term.second));
+  }
+
+  const int num_vars = integer_variables_.size();
+  for (int i = 0; i < num_vars; i++) {
+    const IntegerVariable cp_var = integer_variables_[i];
+    const double lb = ToDouble(integer_trail_->LevelZeroLowerBound(cp_var));
+    const double ub = ToDouble(integer_trail_->LevelZeroUpperBound(cp_var));
+    lp_data_.SetVariableBounds(glop::ColIndex(i), lb, ub);
+  }
+
+  objective_normalization_factor_ = 0.0;
+  glop::ColIndexVector integer_variables;
+  const ColIndex num_cols = lp_data_.num_variables();
+  for (ColIndex col : lp_data_.IntegerVariablesList()) {
+    if (!lp_data_.IsVariableBinary(col)) {
+      integer_variables.push_back(col);
+    }
+
+    // The aim of this normalization value is to compute a coefficient of the
+    // d_i variables that should be minimized.
+    objective_normalization_factor_ +=
+        std::abs(lp_data_.GetObjectiveCoefficientForMinimizationVersion(col));
+  }
+  CHECK_GT(lp_data_.IntegerVariablesList().size(), 0);
+  objective_normalization_factor_ =
+      objective_normalization_factor_ / lp_data_.IntegerVariablesList().size();
+
+  if (!integer_variables.empty()) {
+    // Update the LpProblem with norm variables and constraints.
+    norm_variables_.assign(num_cols, ColIndex(-1));
+    norm_lhs_constraints_.assign(num_cols, RowIndex(-1));
+    norm_rhs_constraints_.assign(num_cols, RowIndex(-1));
+    // For each integer non-binary variable x_i we introduce one new variable
+    // d_i subject to two new constraints:
+    //   d_i - x_i >= -round(x'_i)
+    //   d_i + x_i >= +round(x'_i)
+    // That's round(x'_i) - d_i <= x_i <= round(x'_i) + d_i, where d_i is an
+    // unbounded non-negative, and x'_i is the value of variable i from the
+    // previous solution obtained during the projection step. Consequently
+    // coefficients of the constraints are set here, but bounds of the
+    // constraints are updated at each iteration of the feasibility pump. Also
+    // coefficients of the objective are set here: x_i's are not present in the
+    // objective (i.e., coefficients set to 0.0), and d_i's are present in the
+    // objective with coefficients set to 1.0.
+    // Note that the treatment of integer non-binary variables is different
+    // from the treatment of binary variables. Binary variables do not impose
+    // any extra variables, nor extra constraints, but their objective
+    // coefficients are changed in the linear projection steps.
+    for (const ColIndex col : integer_variables) {
+      const ColIndex norm_variable = lp_data_.CreateNewVariable();
+      norm_variables_[col] = norm_variable;
+      lp_data_.SetVariableBounds(norm_variable, 0.0, glop::kInfinity);
+      const RowIndex row_a = lp_data_.CreateNewConstraint();
+      norm_lhs_constraints_[col] = row_a;
+      lp_data_.SetCoefficient(row_a, norm_variable, 1.0);
+      lp_data_.SetCoefficient(row_a, col, -1.0);
+      const RowIndex row_b = lp_data_.CreateNewConstraint();
+      norm_rhs_constraints_[col] = row_b;
+      lp_data_.SetCoefficient(row_b, norm_variable, 1.0);
+      lp_data_.SetCoefficient(row_b, col, 1.0);
+    }
+  }
+
+  // TODO(user): Try to add scaling.
+  lp_data_.AddSlackVariablesWhereNecessary(
+      /*detect_integer_constraints=*/false);
+}
+
+void FeasibilityPump::L1DistanceMinimize() {
+  std::vector<double> new_obj_coeffs(lp_data_.num_variables().value(), 0.0);
+
+  // Set the original subobjective. The coefficients are scaled by mixing factor
+  // and the offset remains at 0 (because it does not affect the solution).
+  const ColIndex num_cols(lp_data_.objective_coefficients().size());
+  for (ColIndex col(0); col < num_cols; ++col) {
+    new_obj_coeffs[col.value()] =
+        mixing_factor_ * lp_data_.objective_coefficients()[col];
+  }
+
+  // Set the norm subobjective. The coefficients are scaled by 1 - mixing factor
+  // and the offset remains at 0 (because it does not affect the solution).
+  for (const ColIndex col : lp_data_.IntegerVariablesList()) {
+    if (lp_data_.IsVariableBinary(col)) {
+      const Fractional objective_coefficient =
+          mixing_factor_ * lp_data_.objective_coefficients()[col] +
+          (1 - mixing_factor_) * objective_normalization_factor_ *
+              (1 - 2 * integer_solution_[col.value()]);
+      new_obj_coeffs[col.value()] = objective_coefficient;
+    } else {  // The variable is integer.
+      // Update the bounds of the constraints added in
+      // InitializeIntegerVariables() (see there for more details):
+      //   d_i - x_i >= -round(x'_i)
+      //   d_i + x_i >= +round(x'_i)
+
+      // TODO(user): We change both the objective and the bounds, thus
+      // breaking the incrementality. Handle integer variables differently,
+      // e.g., intensify rounding, or use soft fixing from: Fischetti, Lodi,
+      // "Local Branching", Math Program Ser B 98:23-47 (2003).
+      const Fractional objective_coefficient =
+          (1 - mixing_factor_) * objective_normalization_factor_;
+      new_obj_coeffs[norm_variables_[col].value()] = objective_coefficient;
+      // At this point, constraint bounds have already been transformed into
+      // bounds of slack variables. Instead of updating the constraints, we need
+      // to update the slack variables corresponding to them.
+      const ColIndex norm_a_slack_variable =
+          lp_data_.GetSlackVariable(norm_lhs_constraints_[col]);
+      lp_data_.SetVariableBounds(norm_a_slack_variable, -glop::kInfinity,
+                                 integer_solution_[col.value()]);
+      const ColIndex norm_b_slack_variable =
+          lp_data_.GetSlackVariable(norm_rhs_constraints_[col]);
+      lp_data_.SetVariableBounds(norm_b_slack_variable, -glop::kInfinity,
+                                 -integer_solution_[col.value()]);
+    }
+  }
+  for (ColIndex col(0); col < lp_data_.num_variables(); ++col) {
+    lp_data_.SetObjectiveCoefficient(col, new_obj_coeffs[col.value()]);
+  }
+  // TODO(user): Tune this or expose as parameter.
+  mixing_factor_ *= 0.8;
+}
+
+bool FeasibilityPump::SolveLp() {
+  const int num_vars = integer_variables_.size();
+  VLOG(3) << "LP relaxation: " << lp_data_.GetDimensionString() << ".";
+
+  const auto status = simplex_.Solve(lp_data_, time_limit_);
+  total_num_simplex_iterations_ += simplex_.GetNumberOfIterations();
+  if (!status.ok()) {
+    VLOG(1) << "The LP solver encountered an error: " << status.error_message();
+    simplex_.ClearStateForNextSolve();
+    return false;
+  }
+
+  VLOG(3) << "simplex status: " << simplex_.GetProblemStatus();
+  lp_solution_fractionality_ = 0.0;
+  if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL ||
+      simplex_.GetProblemStatus() == glop::ProblemStatus::DUAL_FEASIBLE ||
+      simplex_.GetProblemStatus() == glop::ProblemStatus::PRIMAL_FEASIBLE ||
+      simplex_.GetProblemStatus() == glop::ProblemStatus::IMPRECISE) {
+    lp_solution_is_set_ = true;
+    for (int i = 0; i < num_vars; i++) {
+      const double value = simplex_.GetVariableValue(glop::ColIndex(i));
+      lp_solution_[i] = value;
+      lp_solution_fractionality_ = std::max(
+          lp_solution_fractionality_, std::abs(value - std::round(value)));
+    }
+
+    // Compute the objective value.
+    lp_objective_ = 0;
+    for (const auto& term : integer_objective_) {
+      lp_objective_ += lp_solution_[term.first.value()] * term.second.value();
+    }
+    lp_solution_is_integer_ = lp_solution_fractionality_ < kCpEpsilon;
+  }
+  return true;
+}
+
+void FeasibilityPump::UpdateBoundsOfLpVariables() {
+  const int num_vars = integer_variables_.size();
+  for (int i = 0; i < num_vars; i++) {
+    const IntegerVariable cp_var = integer_variables_[i];
+    const double lb = ToDouble(integer_trail_->LowerBound(cp_var));
+    const double ub = ToDouble(integer_trail_->UpperBound(cp_var));
+    lp_data_.SetVariableBounds(glop::ColIndex(i), lb, ub);
+  }
+}
+
+double FeasibilityPump::GetLPSolutionValue(IntegerVariable variable) const {
+  return lp_solution_[gtl::FindOrDie(mirror_lp_variable_, variable).value()];
+}
+
+// ----------------------------------------------------------------
+// -------------------Rounding-------------------------------------
+// ----------------------------------------------------------------
+
+int64 FeasibilityPump::GetIntegerSolutionValue(IntegerVariable variable) const {
+  return integer_solution_[gtl::FindOrDie(mirror_lp_variable_, variable)
+                               .value()];
+}
+
+void FeasibilityPump::SimpleRounding() {
+  if (!lp_solution_is_set_) return;
+  integer_solution_is_set_ = true;
+  for (int i = 0; i < lp_solution_.size(); ++i) {
+    integer_solution_[i] = static_cast<int64>(std::round(lp_solution_[i]));
+  }
+  // Compute the objective value.
+  integer_solution_objective_ = 0;
+  for (const auto& term : integer_objective_) {
+    integer_solution_objective_ +=
+        integer_solution_[term.first.value()] * term.second.value();
+  }
+
+  integer_solution_is_feasible_ = true;
+  num_infeasible_constraints_ = 0;
+  integer_solution_infeasibility_ = 0;
+  for (glop::RowIndex i(0); i < integer_lp_.size(); ++i) {
+    int64 activity = 0;
+    for (const auto& term : integer_lp_[i].terms) {
+      const int64 prod =
+          CapProd(integer_solution_[term.first.value()], term.second.value());
+      if (prod <= kint64min || prod >= kint64max) {
+        activity = prod;
+        break;
+      }
+      activity = CapAdd(activity, prod);
+      if (activity <= kint64min || activity >= kint64max) break;
+    }
+    if (activity > integer_lp_[i].ub || activity < integer_lp_[i].lb) {
+      integer_solution_is_feasible_ = false;
+      num_infeasible_constraints_++;
+      integer_solution_infeasibility_ =
+          std::max(integer_solution_infeasibility_,
+                   std::max(activity - integer_lp_[i].ub.value(),
+                            integer_lp_[i].lb.value() - activity));
+    }
+  }
+}
+
+}  // namespace sat
+}  // namespace operations_research

ortools/sat/feasibility_pump.h

@@ -0,0 +1,185 @@
+// Copyright 2010-2018 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef OR_TOOLS_SAT_FEASIBILITY_PUMP_H_
+#define OR_TOOLS_SAT_FEASIBILITY_PUMP_H_
+
+#include "ortools/glop/revised_simplex.h"
+#include "ortools/lp_data/lp_data.h"
+#include "ortools/lp_data/lp_data_utils.h"
+#include "ortools/lp_data/lp_types.h"
+#include "ortools/sat/cp_model_loader.h"
+#include "ortools/sat/integer.h"
+#include "ortools/sat/linear_constraint.h"
+#include "ortools/sat/synchronization.h"
+#include "ortools/sat/util.h"
+
+namespace operations_research {
+namespace sat {
+
+struct FeasibilityPumpLpSolution
+    : public gtl::ITIVector<IntegerVariable, double> {
+  FeasibilityPumpLpSolution() {}
+};
+
+class FeasibilityPump {
+ public:
+  explicit FeasibilityPump(Model* model);
+  ~FeasibilityPump();
+
+  typedef glop::RowIndex ConstraintIndex;
+
+  void SetMaxFPIterations(int max_iter) {
+    max_fp_iterations_ = std::max(1, max_iter);
+  }
+
+  // Add a new linear constraint to this LP.
+  void AddLinearConstraint(const LinearConstraint& ct);
+
+  // Set the coefficient of the variable in the objective. Calling it twice will
+  // overwrite the previous value. Note that this doesn't set the objective
+  // coefficient if the variable doesn't appear in any constraints. So this has
+  // to be called after all the constraints are added.
+  void SetObjectiveCoefficient(IntegerVariable ivar, IntegerValue coeff);
+
+  // Returns the LP value of a variable in the current
+  // solution. These functions should only be called when HasSolution() is true.
+  bool HasLPSolution() const { return lp_solution_is_set_; }
+  double LPSolutionObjectiveValue() const { return lp_objective_; }
+  double GetLPSolutionValue(IntegerVariable variable) const;
+  bool LPSolutionIsInteger() const { return lp_solution_is_integer_; }
+  double LPSolutionFractionality() const { return lp_solution_fractionality_; }
+
+  // Returns the Integer solution value of a variable in the current rounded
+  // solution. These functions should only be called when HasIntegerSolution()
+  // is true.
+  bool HasIntegerSolution() const { return integer_solution_is_set_; }
+  int64 IntegerSolutionObjectiveValue() const {
+    return integer_solution_objective_;
+  }
+  bool IntegerSolutionIsFeasible() const {
+    return integer_solution_is_feasible_;
+  }
+  int64 GetIntegerSolutionValue(IntegerVariable variable) const;
+
+  void Solve();
+
+ private:
+  // Solve the LP, returns false if something went wrong in the LP solver.
+  bool SolveLp();
+
+  // Round the fractional LP solution values to nearest integer values.
+  void SimpleRounding();
+
+  // Loads the lp_data_.
+  void InitializeWorkingLP();
+
+  // Changes the LP objective and bounds of the norm constraints so the new
+  // objective also tries to minimize the distance to the rounded solution.
+  void L1DistanceMinimize();
+
+  // Stores the solutions in the shared repository. Stores LP solution if it is
+  // integer and stores the integer solution if it is feasible.
+  void MaybePushToRepo();
+
+  void PrintStats();
+
+  // Shortcut for an integer linear expression type.
+  using LinearExpression = std::vector<std::pair<glop::ColIndex, IntegerValue>>;
+
+  // Gets or creates an LP variable that mirrors a model variable.
+  // The variable should be a positive reference.
+  glop::ColIndex GetOrCreateMirrorVariable(IntegerVariable positive_variable);
+
+  // Updates the bounds of the LP variables from the CP bounds.
+  void UpdateBoundsOfLpVariables();
+
+  // This epsilon is related to the precision of the value returned by the LP
+  // once they have been scaled back into the CP domain. So for large domain or
+  // cost coefficient, we may have some issues.
+  static const double kCpEpsilon;
+
+  // Initial problem in integer form.
+  // We always sort the inner vectors by increasing glop::ColIndex.
+  struct LinearConstraintInternal {
+    IntegerValue lb;
+    IntegerValue ub;
+    LinearExpression terms;
+  };
+  LinearExpression integer_objective_;
+  IntegerValue objective_infinity_norm_ = IntegerValue(0);
+  double objective_normalization_factor_ = 0.0;
+  double mixing_factor_ = 1.0;
+
+  gtl::ITIVector<glop::RowIndex, LinearConstraintInternal> integer_lp_;
+  int model_vars_size_ = 0;
+
+  // Underlying LP solver API.
+  glop::LinearProgram lp_data_;
+  glop::RevisedSimplex simplex_;
+
+  glop::ColMapping norm_variables_;
+  glop::ColToRowMapping norm_lhs_constraints_;
+  glop::ColToRowMapping norm_rhs_constraints_;
+
+  // Structures used for mirroring IntegerVariables inside the underlying LP
+  // solver: an integer variable var is mirrored by mirror_lp_variable_[var].
+  // Note that these indices are dense in [0, mirror_lp_variable_.size()] so
+  // they can be used as vector indices.
+  std::vector<IntegerVariable> integer_variables_;
+  absl::flat_hash_map<IntegerVariable, glop::ColIndex> mirror_lp_variable_;
+
+  // We need to remember what to optimize if an objective is given, because
+  // then we will switch the objective between feasibility and optimization.
+  bool objective_is_defined_ = false;
+
+  // Singletons from Model.
+  const SatParameters& sat_parameters_;
+  TimeLimit* time_limit_;
+  IntegerTrail* integer_trail_;
+  Trail* trail_;
+  IntegerEncoder* integer_encoder_;
+  SharedIncompleteSolutionManager* incomplete_solutions_;
+  const CpModelMapping* mapping_;
+
+  // Last OPTIMAL/Feasible solution found by a call to the underlying LP solver.
+  bool lp_solution_is_set_ = false;
+  bool lp_solution_is_integer_ = false;
+  double lp_objective_;
+  std::vector<double> lp_solution_;
+  std::vector<double> best_lp_solution_;
+  // We use max fractionality of all variables.
+  double lp_solution_fractionality_;
+
+  // Rounded Integer solution. This might not be feasible.
+  bool integer_solution_is_set_ = false;
+  bool integer_solution_is_feasible_ = false;
+  int64 integer_solution_objective_;
+  std::vector<int64> integer_solution_;
+  std::vector<int64> best_integer_solution_;
+  int num_infeasible_constraints_;
+  // We use max infeasibility of all constraints.
+  int64 integer_solution_infeasibility_;
+
+  // Sum of all simplex iterations performed by this class. This is useful to
+  // test the incrementality and compare to other solvers.
+  int64 total_num_simplex_iterations_ = 0;
+
+  // TODO(user): Tune default value. Expose as parameter.
+  int max_fp_iterations_ = 20;
+};
+
+}  // namespace sat
+}  // namespace operations_research
+
+#endif  // OR_TOOLS_SAT_FEASIBILITY_PUMP_H_

ortools/sat/rins.cc

@@ -85,33 +85,50 @@ std::vector<double> GetGeneralRelaxationValues(
   }
   return relaxation_values;
 }
+
+std::vector<double> GetIncompleteSolutionValues(
+    SharedIncompleteSolutionManager* incomplete_solutions) {
+  std::vector<double> empty_solution_values;
+
+  if (incomplete_solutions == nullptr ||
+      !incomplete_solutions->HasNewSolution()) {
+    return empty_solution_values;
+  }
+
+  return incomplete_solutions->GetNewSolution();
+}
 }  // namespace
 
 RINSNeighborhood GetRINSNeighborhood(
     const SharedResponseManager* response_manager,
     const SharedRelaxationSolutionRepository* relaxation_solutions,
-    const SharedLPSolutionRepository* lp_solutions, random_engine_t* random) {
+    const SharedLPSolutionRepository* lp_solutions,
+    SharedIncompleteSolutionManager* incomplete_solutions,
+    random_engine_t* random) {
   RINSNeighborhood rins_neighborhood;
 
   const bool use_only_relaxation_values =
       (response_manager == nullptr ||
        response_manager->SolutionsRepository().NumSolutions() == 0);
 
-  if (use_only_relaxation_values && lp_solutions == nullptr) {
+  if (use_only_relaxation_values && lp_solutions == nullptr &&
+      incomplete_solutions == nullptr) {
     // As of now RENS doesn't generate good neighborhoods from integer
     // relaxation solutions.
     return rins_neighborhood;
   }
 
   std::vector<double> relaxation_values;
-  if (lp_solutions == nullptr) {
+  if (incomplete_solutions != nullptr) {
+    relaxation_values = GetIncompleteSolutionValues(incomplete_solutions);
+  } else if (lp_solutions != nullptr) {
+    relaxation_values = GetLPRelaxationValues(lp_solutions, random);
+  } else {
     CHECK(relaxation_solutions != nullptr)
         << "No relaxation solutions repository or lp solutions repository "
            "provided.";
     relaxation_values =
         GetGeneralRelaxationValues(relaxation_solutions, random);
-  } else {
-    relaxation_values = GetLPRelaxationValues(lp_solutions, random);
   }
   if (relaxation_values.empty()) return rins_neighborhood;
 

ortools/sat/rins.h

@@ -60,17 +60,24 @@ struct RINSNeighborhood {
 };
 
 // Helper method to create a RINS neighborhood by fixing variables with same
-// values in relaxation solution and last solution. Uses lp values if
-// 'use_lp_relaxation' is true, otherwise uses general relaxation value stored
-// in the repository.
+// values in relaxation solution and the current best solution in the
+// response_manager. Prioritizes repositories in following order to get a
+// relaxation solution.
+//  1. incomplete_solutions
+//  2. lp_solutions
+//  3. relaxation_solutions
 //
-// If use_only_relaxation_values is true, this Generates a RENS neighborhood by
+// If response_manager is not provided, this generates a RENS neighborhood by
 // ignoring the solutions and using the relaxation values. The domain of the
-// variables are reduced to integer values around relaxation values.
+// variables are reduced to integer values around relaxation values. If the
+// relaxation value is integer, then we fix the domain of the variable to that
+// value.
 RINSNeighborhood GetRINSNeighborhood(
     const SharedResponseManager* response_manager,
     const SharedRelaxationSolutionRepository* relaxation_solutions,
-    const SharedLPSolutionRepository* lp_solutions, random_engine_t* random);
+    const SharedLPSolutionRepository* lp_solutions,
+    SharedIncompleteSolutionManager* incomplete_solutions,
+    random_engine_t* random);
 
 // Adds the current LP solution to the pool.
 void RecordLPRelaxationValues(Model* model);

ortools/sat/samples/AssignmentSat.java

@@ -0,0 +1,118 @@
+// Copyright 2010-2018 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// CP-SAT example that solves an assignment problem.
+// [START program]
+package com.google.ortools.sat.samples;
+// [START import]
+import com.google.ortools.sat.CpModel;
+import com.google.ortools.sat.CpSolver;
+import com.google.ortools.sat.CpSolverStatus;
+import com.google.ortools.sat.IntVar;
+import com.google.ortools.sat.LinearExpr;
+// [END import]
+
+/** Assignment problem. */
+public class AssignmentSat {
+  static {
+    System.loadLibrary("jniortools");
+  }
+
+  public static void main(String[] args) {
+    // Data
+    // [START data_model]
+    int[][] costs = {
+        {90, 80, 75, 70},
+        {35, 85, 55, 65},
+        {125, 95, 90, 95},
+        {45, 110, 95, 115},
+        {50, 100, 90, 100},
+    };
+    final int numWorkers = costs.length;
+    final int numTasks = costs[0].length;
+    // [END data_model]
+
+    // Model
+    // [START model]
+    CpModel model = new CpModel();
+    // [END model]
+
+    // Variables
+    // [START variables]
+    IntVar[][] x = new IntVar[numWorkers][numTasks];
+    // Variables in a 1-dim array.
+    IntVar[] xFlat = new IntVar[numWorkers * numTasks];
+    int[] costsFlat = new int[numWorkers * numTasks];
+    for (int i = 0; i < numWorkers; ++i) {
+      for (int j = 0; j < numTasks; ++j) {
+        x[i][j] = model.newIntVar(0, 1, "");
+        int k = i * numTasks + j;
+        xFlat[k] = x[i][j];
+        costsFlat[k] = costs[i][j];
+      }
+    }
+    // [END variables]
+
+    // Constraints
+    // [START constraints]
+    // Each worker is assigned to at most one task.
+    for (int i = 0; i < numWorkers; ++i) {
+      IntVar[] vars = new IntVar[numTasks];
+      for (int j = 0; j < numTasks; ++j) {
+        vars[j] = x[i][j];
+      }
+      model.addLessOrEqual(LinearExpr.sum(vars), 1);
+    }
+    // Each task is assigned to exactly one worker.
+    for (int j = 0; j < numTasks; ++j) {
+      // LinearExpr taskSum;
+      IntVar[] vars = new IntVar[numWorkers];
+      for (int i = 0; i < numWorkers; ++i) {
+        vars[i] = x[i][j];
+      }
+      model.addEquality(LinearExpr.sum(vars), 1);
+    }
+    // [END constraints]
+
+    // Objective
+    // [START objective]
+    model.minimize(LinearExpr.scalProd(xFlat, costsFlat));
+    // [END objective]
+
+    // Solve
+    // [START solve]
+    CpSolver solver = new CpSolver();
+    CpSolverStatus status = solver.solve(model);
+    // [END solve]
+
+    // Print solution.
+    // [START print_solution]
+    // Check that the problem has a feasible solution.
+    if (status == CpSolverStatus.OPTIMAL || status == CpSolverStatus.FEASIBLE) {
+      System.out.println("Total cost: " + solver.objectiveValue() + "\n");
+      for (int i = 0; i < numWorkers; ++i) {
+        for (int j = 0; j < numTasks; ++j) {
+          if (solver.value(x[i][j]) == 1) {
+            System.out.println(
+                "Worker " + i + " assigned to task " + j + ".  Cost: " + costs[i][j]);
+          }
+        }
+      }
+    } else {
+      System.err.println("No solution found.");
+    }
+    // [END print_solution]
+  }
+
+  private AssignmentSat() {}
+}

ortools/sat/samples/assignment_sat.cc

@@ -0,0 +1,110 @@
+// Copyright 2010-2018 Google LLC
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// [START program]
+// [START import]
+#include "ortools/sat/cp_model.h"
+// [END import]
+namespace operations_research {
+namespace sat {
+
+void IntegerProgrammingExample() {
+  // Data
+  // [START data_model]
+  const std::vector<std::vector<double>> costs{
+      {90, 80, 75, 70},   {35, 85, 55, 65},   {125, 95, 90, 95},
+      {45, 110, 95, 115}, {50, 100, 90, 100},
+  };
+  const int num_workers = costs.size();
+  const int num_tasks = costs[0].size();
+  // [END data_model]
+
+  // Model
+  // [START model]
+  CpModelBuilder cp_model;
+  // [END model]
+
+  // Variables
+  // [START variables]
+  // x[i][j] is an array of Boolean variables. x[i][j] is true
+  // if worker i is assigned to task j.
+  std::vector<std::vector<BoolVar>> x(num_workers,
+                                      std::vector<BoolVar>(num_tasks));
+  for (int i = 0; i < num_workers; ++i) {
+    for (int j = 0; j < num_tasks; ++j) {
+      x[i][j] = cp_model.NewBoolVar();
+    }
+  }
+  // [END variables]
+
+  // Constraints
+  // [START constraints]
+  // Each worker is assigned to at most one task.
+  for (int i = 0; i < num_workers; ++i) {
+    LinearExpr worker_sum;
+    for (int j = 0; j < num_tasks; ++j) {
+      worker_sum.AddTerm(x[i][j], 1);
+    }
+    cp_model.AddLessOrEqual(worker_sum, 1);
+  }
+  // Each task is assigned to exactly one worker.
+  for (int j = 0; j < num_tasks; ++j) {
+    LinearExpr task_sum;
+    for (int i = 0; i < num_workers; ++i) {
+      task_sum.AddTerm(x[i][j], 1);
+    }
+    cp_model.AddEquality(task_sum, 1);
+  }
+  // [END constraints]
+
+  // Objective
+  // [START objective]
+  LinearExpr total_cost;
+  for (int i = 0; i < num_workers; ++i) {
+    for (int j = 0; j < num_tasks; ++j) {
+      total_cost.AddTerm(x[i][j], costs[i][j]);
+    }
+  }
+  cp_model.Minimize(total_cost);
+  // [END objective]
+
+  // Solve
+  // [START solve]
+  const CpSolverResponse response = Solve(cp_model.Build());
+  // [END solve]
+
+  // Print solution.
+  // [START print_solution]
+  if (response.status() == CpSolverStatus::INFEASIBLE) {
+    LOG(FATAL) << "No solution found.";
+  }
+
+  LOG(INFO) << "Total cost: " << response.objective_value();
+  LOG(INFO);
+  for (int i = 0; i < num_workers; ++i) {
+    for (int j = 0; j < num_tasks; ++j) {
+      if (SolutionBooleanValue(response, x[i][j])) {
+        LOG(INFO) << "Task " << i << " assigned to worker " << j
+                  << ".  Cost: " << costs[i][j];
+      }
+    }
+  }
+  // [END print_solution]
+}
+}  // namespace sat
+}  // namespace operations_research
+
+int main(int argc, char** argv) {
+  operations_research::sat::IntegerProgrammingExample();
+  return EXIT_SUCCESS;
+}

ortools/sat/samples/assignment_sat.py

@@ -0,0 +1,92 @@
+# Copyright 2010-2018 Google LLC
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Solve a simple assignment problem."""
+# [START program]
+# [START import]
+from ortools.sat.python import cp_model
+# [END import]
+
+
+def main():
+    # Data
+    # [START data_model]
+    costs = [
+        [90, 80, 75, 70],
+        [35, 85, 55, 65],
+        [125, 95, 90, 95],
+        [45, 110, 95, 115],
+        [50, 100, 90, 100],
+    ]
+    num_workers = len(costs)
+    num_tasks = len(costs[0])
+    # [END data_model]
+
+    # Model
+    # [START model]
+    model = cp_model.CpModel()
+    # [END model]
+
+    # Variables
+    # [START variables]
+    x = []
+    for i in range(num_workers):
+        t = []
+        for j in range(num_tasks):
+            t.append(model.NewBoolVar('x[%i,%i]' % (i, j)))
+        x.append(t)
+    # [END variables]
+
+    # Constraints
+    # [START constraints]
+    # Each worker is assigned to at most one task.
+    for i in range(num_workers):
+        model.Add(sum(x[i][j] for j in range(num_tasks)) <= 1)
+
+    # Each task is assigned to exactly one worker.
+    for j in range(num_tasks):
+        model.Add(sum(x[i][j] for i in range(num_workers)) == 1)
+    # [END constraints]
+
+    # Objective
+    # [START objective]
+    objective_terms = []
+    for i in range(num_workers):
+        for j in range(num_tasks):
+            objective_terms.append(costs[i][j] * x[i][j])
+    model.Minimize(sum(objective_terms))
+    # [END objective]
+
+    # Solve
+    # [START solve]
+    solver = cp_model.CpSolver()
+    status = solver.Solve(model)
+    # [END solve]
+
+    # Print solution.
+    # [START print_solution]
+    if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
+        print('Total cost = %i' % solver.ObjectiveValue())
+        print()
+        for i in range(num_workers):
+            for j in range(num_tasks):
+                if solver.BooleanValue(x[i][j]):
+                    print('Worker ', i, ' assigned to task ', j, '  Cost = ',
+                          costs[i][j])
+    else:
+        print('No solution found.')
+    # [END print_solution]
+
+
+if __name__ == '__main__':
+    main()
+# [END program]

ortools/sat/sat_parameters.proto

@@ -21,7 +21,7 @@ option java_multiple_files = true;
 // Contains the definitions for all the sat algorithm parameters and their
 // default values.
 //
-// NEXT TAG: 164
+// NEXT TAG: 165
 message SatParameters {
   // ==========================================================================
   // Branching and polarity
@@ -808,6 +808,9 @@ message SatParameters {
   // Turns on relaxation induced neighborhood generator.
   optional bool use_rins_lns = 129 [default = true];
 
+  // Adds a feasibility pump subsolver along with lns subsolvers.
+  optional bool use_feasibility_pump = 164 [default = false];
+
   // Turns on a lns worker which solves relaxed version of the original problem
   // by removing constraints from the problem in order to get better bounds.
   optional bool use_relaxation_lns = 150 [default = false];

ortools/sat/synchronization.cc

@@ -63,7 +63,7 @@ void SharedRelaxationSolutionRepository::NewRelaxationSolution(
 }
 
 void SharedLPSolutionRepository::NewLPSolution(
-    std::vector<double> lp_solution) {
+    const std::vector<double>& lp_solution) {
   // Note that the Add() method already applies mutex lock. So we don't need it
   // here.
 
@@ -80,6 +80,27 @@ void SharedLPSolutionRepository::NewLPSolution(
   AddInternal(solution);
 }
 
+bool SharedIncompleteSolutionManager::HasNewSolution() const {
+  absl::MutexLock mutex_lock(&mutex_);
+  return !solutions_.empty();
+}
+
+std::vector<double> SharedIncompleteSolutionManager::GetNewSolution() {
+  absl::MutexLock mutex_lock(&mutex_);
+  std::vector<double> solution;
+  if (solutions_.empty()) return solution;
+
+  solution = std::move(solutions_.back());
+  solutions_.pop_back();
+  return solution;
+}
+
+void SharedIncompleteSolutionManager::AddNewSolution(
+    const std::vector<double>& lp_solution) {
+  absl::MutexLock mutex_lock(&mutex_);
+  solutions_.push_back(lp_solution);
+}
+
 // TODO(user): Experiments and play with the num_solutions_to_keep parameter.
 SharedResponseManager::SharedResponseManager(bool log_updates,
                                              bool enumerate_all_solutions,

ortools/sat/synchronization.h

@@ -136,7 +136,28 @@ class SharedLPSolutionRepository : public SharedSolutionRepository<double> {
   explicit SharedLPSolutionRepository(int num_solutions_to_keep)
       : SharedSolutionRepository<double>(num_solutions_to_keep) {}
 
-  void NewLPSolution(std::vector<double> lp_solution);
+  void NewLPSolution(const std::vector<double>& lp_solution);
+};
+
+// Set of partly filled solutions. They are meant to be finished by some lns
+// worker.
+//
+// The solutions are stored as a vector of doubles. The value at index i
+// represents the solution value of model variable indexed i. Note that some
+// values can be infinity which should be interpreted as 'unknown' solution
+// value for that variable. These solutions can not necessarily be completed to
+// complete feasible solutions.
+class SharedIncompleteSolutionManager {
+ public:
+  bool HasNewSolution() const;
+  std::vector<double> GetNewSolution();
+
+  void AddNewSolution(const std::vector<double>& lp_solution);
+
+ private:
+  // New solutions are added and removed from the back.
+  std::vector<std::vector<double>> solutions_;
+  mutable absl::Mutex mutex_;
 };
 
 // Manages the global best response kept by the solver.