From 5964c5e2ddfd52b574ae039b3b8f08a3c780454e Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Fri, 8 Jan 2021 09:49:56 +0100 Subject: [PATCH] improve primal integral computation --- ortools/sat/cp_model.proto | 2 +- ortools/sat/cp_model_presolve.cc | 2 +- ortools/sat/cp_model_solver.cc | 24 +++++++++++++++++-- ortools/sat/samples/minimal_jobshop_sat.py | 15 ++++++------ ortools/sat/synchronization.cc | 27 ++++++++++++++++++---- ortools/sat/synchronization.h | 11 +++++++++ 6 files changed, 66 insertions(+), 15 deletions(-) diff --git a/ortools/sat/cp_model.proto b/ortools/sat/cp_model.proto index da746ef91a..03860919c8 100644 --- a/ortools/sat/cp_model.proto +++ b/ortools/sat/cp_model.proto @@ -558,7 +558,7 @@ enum CpSolverStatus { // detailed error by calling ValidateCpModel(model_proto). MODEL_INVALID = 1; - // A feasible solution as been found. But the search was stopped before we + // A feasible solution has been found. But the search was stopped before we // could prove optimality or before we enumerated all solutions of a // feasibility problem (if asked). FEASIBLE = 2; diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index c5f5f7196d..c06678f6bf 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -1706,7 +1706,7 @@ bool CpModelPresolver::PropagateDomainsInLinear(int c, ConstraintProto* ct) { // Can we perform some substitution? // - // TODO(user): there is no guarante we will not miss some since we might + // TODO(user): there is no guarantee we will not miss some since we might // not reprocess a constraint once other have been deleted. // Skip affine constraint. It is more efficient to substitute them lazily diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 5470e13ed7..7c33852add 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -2878,6 +2878,9 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { *model->GetOrCreate()); SharedTimeLimit shared_time_limit(model->GetOrCreate()); #if defined(_MSC_VER) + // On windows, The final_response is optimized out in the return part, and is + // swapped out before the cleanup callback is called. A workaround is to + // create a unique ptr that will forbid this optimization. std::unique_ptr final_response_ptr(new CpSolverResponse()); CpSolverResponse& final_response = *final_response_ptr.get(); #else @@ -3093,8 +3096,7 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { const auto& observers = model->GetOrCreate()->observers; if (!observers.empty()) { shared_response_manager.AddSolutionCallback( - [&model_proto, &observers, &wall_timer, &user_timer, - &postprocess_solution, &shared_time_limit]( + [&model_proto, &observers, &postprocess_solution]( const CpSolverResponse& response_of_presolved_problem) { CpSolverResponse response = response_of_presolved_problem; postprocess_solution(&response); @@ -3113,6 +3115,23 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { }); } + // If specified, we load the initial objective domain right away in the + // response manager. Note that the presolve will always fill it with the + // trivial min/max value if the user left it empty. This avoids to display + // [-infinity, infinity] for the initial objective search space. + if (new_cp_model_proto.has_objective()) { + const Domain domain = ReadDomainFromProto(new_cp_model_proto.objective()); + if (!domain.IsEmpty()) { + shared_response_manager.UpdateInnerObjectiveBounds( + "initial domain", IntegerValue(domain.Min()), + IntegerValue(domain.Max())); + } + } + + // Start counting the primal integral from the current determistic time and + // initial objective domain gap that we just filled. + shared_response_manager.UpdatePrimalIntegral(); + #if !defined(__PORTABLE_PLATFORM__) if (absl::GetFlag(FLAGS_cp_model_dump_models)) { const std::string presolved_file = absl::StrCat( @@ -3170,6 +3189,7 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { LOG(INFO) << absl::StrFormat("*** starting to load the model at %.2fs", wall_timer.Get()); } + shared_response_manager.SetUpdatePrimalIntegralOnEachChange(true); LoadCpModel(new_cp_model_proto, &shared_response_manager, model); shared_response_manager.LoadDebugSolution(model); if (log_search) { diff --git a/ortools/sat/samples/minimal_jobshop_sat.py b/ortools/sat/samples/minimal_jobshop_sat.py index 0992660d09..4189652bc0 100644 --- a/ortools/sat/samples/minimal_jobshop_sat.py +++ b/ortools/sat/samples/minimal_jobshop_sat.py @@ -60,8 +60,9 @@ def MinimalJobshopSat(): end_var = model.NewIntVar(0, horizon, 'end' + suffix) interval_var = model.NewIntervalVar(start_var, duration, end_var, 'interval' + suffix) - all_tasks[job_id, task_id] = task_type( - start=start_var, end=end_var, interval=interval_var) + all_tasks[job_id, task_id] = task_type(start=start_var, + end=end_var, + interval=interval_var) machine_to_intervals[machine].append(interval_var) # [END variables] @@ -101,11 +102,11 @@ def MinimalJobshopSat(): for task_id, task in enumerate(job): machine = task[0] assigned_jobs[machine].append( - assigned_task_type( - start=solver.Value(all_tasks[job_id, task_id].start), - job=job_id, - index=task_id, - duration=task[1])) + assigned_task_type(start=solver.Value( + all_tasks[job_id, task_id].start), + job=job_id, + index=task_id, + duration=task[1])) # Create per machine output lines. output = '' diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index 4870e9fa18..74dd862b49 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -132,13 +132,21 @@ void LogNewSatSolution(const std::string& event_or_solution_count, } // namespace +void SharedResponseManager::SetUpdatePrimalIntegralOnEachChange(bool set) { + absl::MutexLock mutex_lock(&mutex_); + update_integral_on_each_change_ = set; +} + void SharedResponseManager::UpdatePrimalIntegral() { absl::MutexLock mutex_lock(&mutex_); + UpdatePrimalIntegralInternal(); +} + +void SharedResponseManager::UpdatePrimalIntegralInternal() { if (!model_proto_.has_objective()) return; const double current_time = shared_time_limit_->GetElapsedDeterministicTime(); const double time_delta = current_time - last_primal_integral_time_stamp_; - last_primal_integral_time_stamp_ = current_time; // We use the log of the absolute objective gap. // @@ -148,10 +156,14 @@ void SharedResponseManager::UpdatePrimalIntegral() { const CpObjectiveProto& obj = model_proto_.objective(); const double factor = obj.scaling_factor() != 0.0 ? std::abs(obj.scaling_factor()) : 1.0; - const double bounds_delta = std::log( - 1 + factor * std::abs(static_cast(inner_objective_upper_bound_) - - static_cast(inner_objective_lower_bound_))); + const double bounds_delta = std::log(1 + factor * last_absolute_gap_); primal_integral_ += time_delta * bounds_delta; + + // Update with new value. + last_primal_integral_time_stamp_ = current_time; + last_absolute_gap_ = + std::max(0.0, static_cast(inner_objective_upper_bound_) - + static_cast(inner_objective_lower_bound_)); } void SharedResponseManager::SetGapLimitsFromParameters( @@ -163,6 +175,11 @@ void SharedResponseManager::SetGapLimitsFromParameters( } void SharedResponseManager::TestGapLimitsIfNeeded() { + // This is called on each internal limit change, so it is a good place to + // update the integral. Note that this is not called at the end of the search + // though. + if (update_integral_on_each_change_) UpdatePrimalIntegralInternal(); + if (absolute_gap_limit_ == 0 && relative_gap_limit_ == 0) return; if (best_solution_objective_value_ >= kMaxIntegerValue) return; if (inner_objective_lower_bound_ <= kMinIntegerValue) return; @@ -228,6 +245,7 @@ void SharedResponseManager::UpdateInnerObjectiveBounds( } else { best_response_.set_status(CpSolverStatus::INFEASIBLE); } + if (update_integral_on_each_change_) UpdatePrimalIntegralInternal(); if (log_updates_) LogNewSatSolution("Done", wall_timer_.Get(), worker_info); return; } @@ -264,6 +282,7 @@ void SharedResponseManager::NotifyThatImprovingProblemIsInfeasible( // We just proved that the best solution cannot be improved uppon, so we // have a new lower bound. inner_objective_lower_bound_ = best_solution_objective_value_; + if (update_integral_on_each_change_) UpdatePrimalIntegralInternal(); } else { CHECK_EQ(num_solutions_, 0); best_response_.set_status(CpSolverStatus::INFEASIBLE); diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index 40e7f017dc..b49851d0b8 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -215,6 +215,9 @@ class SharedResponseManager { // particular instance. Or to evaluate how efficient our LNS code is improving // solution. // + // Note: The integral will start counting on the first UpdatePrimalIntegral() + // call, since before the difference is assumed to be zero. + // // Important: To report a proper deterministic integral, we only update it // on UpdatePrimalIntegral() which should be called in the main subsolver // synchronization loop. @@ -225,6 +228,11 @@ class SharedResponseManager { double PrimalIntegral() const; void UpdatePrimalIntegral(); + // Sets this to true to have the "real" but non-deterministic primal integral. + // If this is true, then there is no need to manually call + // UpdatePrimalIntegral() but it is not an issue to do so. + void SetUpdatePrimalIntegralOnEachChange(bool set); + // Updates the inner objective bounds. void UpdateInnerObjectiveBounds(const std::string& worker_info, IntegerValue lb, IntegerValue ub); @@ -291,6 +299,7 @@ class SharedResponseManager { ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_); void SetStatsFromModelInternal(Model* model) ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_); + void UpdatePrimalIntegralInternal() ABSL_EXCLUSIVE_LOCKS_REQUIRED(mutex_); const bool log_updates_; const bool enumerate_all_solutions_; @@ -317,7 +326,9 @@ class SharedResponseManager { IntegerValue synchronized_inner_objective_upper_bound_ ABSL_GUARDED_BY(mutex_) = IntegerValue(kint64max); + bool update_integral_on_each_change_ ABSL_GUARDED_BY(mutex_) = false; double primal_integral_ ABSL_GUARDED_BY(mutex_) = 0.0; + double last_absolute_gap_ ABSL_GUARDED_BY(mutex_) = 0.0; double last_primal_integral_time_stamp_ ABSL_GUARDED_BY(mutex_) = 0.0; int next_callback_id_ ABSL_GUARDED_BY(mutex_) = 0;