From 33314b18699f2e360aa6ab796c56feb67d508e5c Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Wed, 5 Jul 2017 16:27:00 -0700 Subject: [PATCH] cleanup sat code after change to allow linear objective --- cmake/utils.cmake | 1 + ortools/sat/boolean_problem.cc | 25 ++--- ortools/sat/cp_model_presolve.cc | 86 ++++++++++------- ortools/sat/cp_model_solver.cc | 74 +++++++++----- ortools/sat/optimization.cc | 159 ++++--------------------------- ortools/sat/optimization.h | 22 +---- 6 files changed, 138 insertions(+), 229 deletions(-) diff --git a/cmake/utils.cmake b/cmake/utils.cmake index f3b0a142d1..68a08d9ce8 100644 --- a/cmake/utils.cmake +++ b/cmake/utils.cmake @@ -33,6 +33,7 @@ FUNCTION(GET_VERSION_FROM_GIT OUTPUT_VAR VERSION_MAJOR_ VERSION_MINOR_ VERSION_P ERROR_QUIET) STRING(STRIP _VERSION_PATCH ${_VERSION_PATCH}) STRING(REGEX REPLACE "\n$" "" _VERSION_PATCH ${_VERSION_PATCH}) + STRING(REGEX REPLACE " " "" _VERSION_PATCH ${_VERSION_PATCH}) STRING(REGEX REPLACE "^v([0-9]+)\\..*" "\\1" _VERSION_MAJOR "${VERSION_FULL}") STRING(REGEX REPLACE "^v[0-9]+\\.([0-9]+).*" "\\1" _VERSION_MINOR "${VERSION_FULL}") diff --git a/ortools/sat/boolean_problem.cc b/ortools/sat/boolean_problem.cc index 6d22131cc4..449e470cef 100644 --- a/ortools/sat/boolean_problem.cc +++ b/ortools/sat/boolean_problem.cc @@ -232,20 +232,21 @@ bool LoadAndConsumeBooleanProblem(LinearBooleanProblem* problem, void UseObjectiveForSatAssignmentPreference(const LinearBooleanProblem& problem, SatSolver* solver) { const LinearObjective& objective = problem.objective(); - double max_weight = 0; - for (int i = 0; i < objective.literals_size(); ++i) { - max_weight = std::max(max_weight, - fabs(static_cast(objective.coefficients(i)))); + CHECK_EQ(objective.literals_size(), objective.coefficients_size()); + int64 max_abs_weight = 0; + for (const int64 coefficient : objective.coefficients()) { + max_abs_weight = std::max(max_abs_weight, std::abs(coefficient)); } + const double max_abs_weight_double = max_abs_weight; for (int i = 0; i < objective.literals_size(); ++i) { - const double weight = - fabs(static_cast(objective.coefficients(i))) / max_weight; - if (objective.coefficients(i) > 0) { - solver->SetAssignmentPreference(Literal(objective.literals(i)).Negated(), - weight); - } else { - solver->SetAssignmentPreference(Literal(objective.literals(i)), weight); - } + const Literal literal(objective.literals(i)); + const int64 coefficient = objective.coefficients(i); + const double abs_weight = std::abs(coefficient) / max_abs_weight_double; + // Because this is a minimization problem, we prefer to assign a Boolean + // variable to its "low" objective value. So if a literal has a positive + // weight when true, we want to set it to false. + solver->SetAssignmentPreference( + coefficient > 0 ? literal.Negated() : literal, abs_weight); } } diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 2a4249bfef..73fb8ab79f 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -1474,15 +1474,14 @@ void PresolveCpModel(const CpModelProto& initial_model, // TODO(user): Insert in main loop. if (context.working_model->has_objective() && context.working_model->objective().vars_size() == 1) { - CpObjectiveProto* const mutable_objective = - context.working_model->mutable_objective(); - const int initial_obj_var = mutable_objective->vars(0); - const int64 initial_coeff = mutable_objective->coeffs(0); - const double initial_offset = mutable_objective->offset(); + const int initial_obj_ref = context.working_model->objective().vars(0); + // TODO(user): Expand the linear equation recursively in order to have // as much term as possible? This would also enable expanding an objective // with multiple terms. int expanded_linear_index = -1; + int64 objective_coeff_in_expanded_constraint; + int64 size_of_expanded_constraint = 0; for (int ct_index = 0; ct_index < context.working_model->constraints_size(); ++ct_index) { const ConstraintProto& ct = context.working_model->constraints(ct_index); @@ -1494,59 +1493,72 @@ void PresolveCpModel(const CpModelProto& initial_model, if (ct.linear().domain().size() != 2) continue; if (ct.linear().domain(0) != ct.linear().domain(1)) continue; - // Find out if initial_obj_var appear in this constraint. + // Find out if initial_obj_ref appear in this constraint. bool present = false; int64 objective_coeff; const int num_terms = ct.linear().vars_size(); for (int i = 0; i < num_terms; ++i) { const int ref = ct.linear().vars(i); const int64 coeff = ct.linear().coeffs(i); - if (PositiveRef(ref) == PositiveRef(initial_obj_var)) { + if (PositiveRef(ref) == PositiveRef(initial_obj_ref)) { CHECK(!present) << "Duplicate variables not supported"; present = true; - objective_coeff = ref == initial_obj_var ? coeff : -coeff; + objective_coeff = (ref == initial_obj_ref) ? coeff : -coeff; } } // We use the longest equality we can find. // TODO(user): Deal with objective_coeff with a magnitude greater than 1? - // Accept when initial_coeff divides objective_coeff. if (present && std::abs(objective_coeff) == 1 && - num_terms > mutable_objective->vars_size() + 1) { + num_terms > size_of_expanded_constraint) { expanded_linear_index = ct_index; - mutable_objective->clear_coeffs(); - mutable_objective->clear_vars(); - const int64 rhs = ct.linear().domain(0); - if (rhs != 0) { - mutable_objective->set_offset(rhs * initial_coeff * objective_coeff + - initial_offset); - } - for (int i = 0; i < num_terms; ++i) { - const int ref = ct.linear().vars(i); - if (PositiveRef(ref) != PositiveRef(initial_obj_var)) { - mutable_objective->add_vars(ref); - mutable_objective->add_coeffs( - -1 * initial_coeff * ct.linear().coeffs(i) * objective_coeff); - } - } + size_of_expanded_constraint = num_terms; + objective_coeff_in_expanded_constraint = objective_coeff; } } if (expanded_linear_index != -1) { context.UpdateRuleStats("objective: expanded single objective"); - ConstraintProto* const ct = - context.working_model->mutable_constraints(expanded_linear_index); - // Remove the objective variable special case and make sure the new - // objective variables cannot be removed: - for (int ref : ct->linear().vars()) { - context.var_to_constraints[PositiveRef(ref)].insert(-1); - } - context.var_to_constraints[PositiveRef(initial_obj_var)].erase(-1); - // This function will detect that the old objective is not used - // elsewhere and remove it from the equation. - PresolveLinear(ct, &context); - context.UpdateConstraintVariableUsage(expanded_linear_index); + // Rewrite the objective. + const int64 multiplier = context.working_model->objective().coeffs(0) / + objective_coeff_in_expanded_constraint; + const ConstraintProto& ct = + context.working_model->constraints(expanded_linear_index); + CpObjectiveProto* const mutable_objective = + context.working_model->mutable_objective(); + mutable_objective->set_offset(mutable_objective->offset() + + ct.linear().domain(0) * multiplier); + mutable_objective->clear_coeffs(); + mutable_objective->clear_vars(); + const int num_terms = ct.linear().vars_size(); + for (int i = 0; i < num_terms; ++i) { + const int ref = ct.linear().vars(i); + if (PositiveRef(ref) != PositiveRef(initial_obj_ref)) { + mutable_objective->add_vars(ref); + mutable_objective->add_coeffs(-ct.linear().coeffs(i) * multiplier); + } + } + + // TODO(user): for now this seems to lower our perf on the minzinc + // benchmark. The likely explanation is a different search heuristic, not + // taking into account the objective. In any case, we need to investigate + // more. + if (/*DISABLES CODE*/ false) { + ConstraintProto* const ct = + context.working_model->mutable_constraints(expanded_linear_index); + // Remove the objective variable special case and make sure the new + // objective variables cannot be removed: + for (int ref : ct->linear().vars()) { + context.var_to_constraints[PositiveRef(ref)].insert(-1); + } + context.var_to_constraints[PositiveRef(initial_obj_ref)].erase(-1); + + // This function will detect that the old objective is not used + // elsewhere and remove it from the equation. + PresolveLinear(ct, &context); + context.UpdateConstraintVariableUsage(expanded_linear_index); + } } } diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 2c23def0c3..e7c79a307f 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -81,6 +81,14 @@ VariableUsage ComputeVariableUsage(const CpModelProto& model_proto) { } } + // Make sure a Booleans is created for each [0,1] Boolean variables. + for (int i = 0; i < model_proto.variables_size(); ++i) { + if (model_proto.variables(i).domain_size() != 2) continue; + if (model_proto.variables(i).domain(0) != 0) continue; + if (model_proto.variables(i).domain(1) != 1) continue; + references.literals.insert(i); + } + std::vector used_integers; for (const int var : references.variables) { used_integers.push_back(PositiveRef(var)); @@ -1498,20 +1506,11 @@ void FillSolutionInResponse(const CpModelProto& model_proto, } namespace { -IntegerVariable GetOrCreateVariableEqualToSumOf( - Model* model, const std::vector>& terms) { - if (terms.empty()) return model->Add(ConstantIntegerVariable(0)); - if (terms.size() == 1 && terms.front().second == 1) { - return terms.front().first; - } - if (terms.size() == 1 && terms.front().second == -1) { - return NegationOf(terms.front().first); - } - // Create a new variable equal to the sum, with a tight domain. +IntegerVariable CreateVariableWithTightBound( + Model* model, const std::vector>& terms) { int64 sum_min = 0; int64 sum_max = 0; - for (const std::pair var_coeff : terms) { const int64 min_domain = model->Get(LowerBound(var_coeff.first)); const int64 max_domain = model->Get(UpperBound(var_coeff.first)); @@ -1521,9 +1520,21 @@ IntegerVariable GetOrCreateVariableEqualToSumOf( sum_min += std::min(prod1, prod2); sum_max += std::max(prod1, prod2); } - IntegerVariable new_var = model->Add(NewIntegerVariable(sum_min, sum_max)); + return model->Add(NewIntegerVariable(sum_min, sum_max)); +} - // Link new variables with the linear terms. +IntegerVariable GetOrCreateVariableGreaterOrEqualToSumOf( + Model* model, const std::vector>& terms) { + if (terms.empty()) return model->Add(ConstantIntegerVariable(0)); + if (terms.size() == 1 && terms.front().second == 1) { + return terms.front().first; + } + if (terms.size() == 1 && terms.front().second == -1) { + return NegationOf(terms.front().first); + } + + // Create a new variable and link it with the linear terms. + const IntegerVariable new_var = CreateVariableWithTightBound(model, terms); std::vector vars; std::vector coeffs; for (const auto& term : terms) { @@ -1532,7 +1543,7 @@ IntegerVariable GetOrCreateVariableEqualToSumOf( } vars.push_back(new_var); coeffs.push_back(-1); - model->Add(FixedWeightedSum(vars, coeffs, 0)); + model->Add(WeightedSumLowerOrEqual(vars, coeffs, 0)); return new_var; } } // namespace @@ -1636,15 +1647,21 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto, FindOrDie(representative_to_lp_constraint, id); const std::vector>& terms = it.second; const IntegerVariable sub_obj_var = - GetOrCreateVariableEqualToSumOf(m->model(), terms); + GetOrCreateVariableGreaterOrEqualToSumOf(m->model(), terms); top_level_cp_terms.push_back(std::make_pair(sub_obj_var, 1)); lp->SetMainObjectiveVariable(sub_obj_var); num_components_containing_objective++; } } - const IntegerVariable main_objective_var = - GetOrCreateVariableEqualToSumOf(m->model(), top_level_cp_terms); + IntegerVariable main_objective_var; + if (m->GetOrCreate()->parameters().optimize_with_core()) { + main_objective_var = + CreateVariableWithTightBound(m->model(), top_level_cp_terms); + } else { + main_objective_var = GetOrCreateVariableGreaterOrEqualToSumOf( + m->model(), top_level_cp_terms); + } // Register LP constraints. Note that this needs to be done after all the // constraints have been added. @@ -1652,11 +1669,11 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto, lp_constraint->RegisterWith(m->GetOrCreate()); } + VLOG(1) << top_level_cp_terms.size() + << " terms in the main objective linear equation (" + << num_components_containing_objective << " from LP constraints)."; VLOG_IF(1, !lp_constraints.empty()) << "Added " << lp_constraints.size() << " LP constraints."; - VLOG_IF(1, num_components_containing_objective > 1) - << "Objective split into " << num_components_containing_objective - << " components"; return main_objective_var; } @@ -1847,7 +1864,6 @@ CpSolverResponse SolveCpModelInternal(const CpModelProto& model_proto, // Create an objective variable and its associated linear constraint if // needed. IntegerVariable objective_var = kNoIntegerVariable; - if (parameters.use_global_lp_constraint()) { // Linearize some part of the problem and register LP constraint(s). objective_var = AddLPConstraints(model_proto, &m); @@ -1858,7 +1874,12 @@ CpSolverResponse SolveCpModelInternal(const CpModelProto& model_proto, for (int i = 0; i < obj.vars_size(); ++i) { terms.push_back(std::make_pair(m.Integer(obj.vars(i)), obj.coeffs(i))); } - objective_var = GetOrCreateVariableEqualToSumOf(m.model(), terms); + if (parameters.optimize_with_core()) { + objective_var = CreateVariableWithTightBound(m.model(), terms); + } else { + objective_var = + GetOrCreateVariableGreaterOrEqualToSumOf(m.model(), terms); + } } model->GetOrCreate() @@ -1914,13 +1935,20 @@ CpSolverResponse SolveCpModelInternal(const CpModelProto& model_proto, } else { // Optimization problem. const CpObjectiveProto& obj = model_proto.objective(); + VLOG(1) << obj.vars_size() << " terms in the proto objective."; const auto solution_observer = [&model_proto, &response, &num_solutions, &obj, &m, objective_var](const Model& sat_model) { num_solutions++; FillSolutionInResponse(model_proto, m, &response); + int64 objective_value = 0; + for (int i = 0; i < model_proto.objective().vars_size(); ++i) { + objective_value += model_proto.objective().coeffs(i) * + sat_model.Get(Value( + m.Integer(model_proto.objective().vars(i)))); + } response.set_objective_value( - ScaleObjectiveValue(obj, sat_model.Get(Value(objective_var)))); + ScaleObjectiveValue(obj, objective_value)); VLOG(1) << "Solution #" << num_solutions << " obj:" << response.objective_value() << " num_bool:" << sat_model.Get()->NumVariables(); diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc index bd404cce06..df91f2b592 100644 --- a/ortools/sat/optimization.cc +++ b/ortools/sat/optimization.cc @@ -1175,7 +1175,14 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( int num_solutions = 0; IntegerValue best_objective = integer_trail->UpperBound(objective_var); const auto process_solution = [&]() { - const IntegerValue objective(model->Get(Value(objective_var))); + // We don't assume that objective_var is linked with its linear term, so + // we recompute the objective here. + IntegerValue objective(0); + for (int i = 0; i < variables.size(); ++i) { + objective += + coefficients[i] * IntegerValue(model->Get(Value(variables[i]))); + } + if (objective >= best_objective && num_solutions > 0) return true; ++num_solutions; @@ -1292,6 +1299,12 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( stratified_threshold.value(), max_depth); } + // Abort if we have a solution and a gap of zero. + if (implied_objective_lb == best_objective && num_solutions > 0) { + result = SatSolver::MODEL_SAT; + break; + } + // Solve under the assumptions. std::vector> cores; result = FindCores(assumptions, next_decision, model, &cores); @@ -1451,7 +1464,13 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( int num_solutions = 0; IntegerValue best_objective = integer_trail->UpperBound(objective_var); const auto process_solution = [&]() { - const IntegerValue objective(model->Get(Value(objective_var))); + // We don't assume that objective_var is linked with its linear term, so + // we recompute the objective here. + IntegerValue objective(0); + for (int i = 0; i < variables.size(); ++i) { + objective += + coefficients[i] * IntegerValue(model->Get(Value(variables[i]))); + } if (objective >= best_objective) return; ++num_solutions; best_objective = objective; @@ -1646,141 +1665,5 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( } #endif // defined(USE_CBC) || defined(USE_SCIP) -SatSolver::Status MinimizeWeightedLiteralSumWithCoreAndLazyEncoding( - bool log_info, const std::vector& literals, - const std::vector& int64_coeffs, - const std::function& next_decision, - const std::function& feasible_solution_observer, - Model* model) { - // Timing. - WallTimer wall_timer; - UserTimer user_timer; - wall_timer.Start(); - user_timer.Start(); - - // TODO(user): Rename Coefficient into IntegerValue everywhere. - std::vector coeffs(int64_coeffs.begin(), int64_coeffs.end()); - - // Create one initial nodes per variables with cost. - Coefficient lower_bound(0); - Coefficient upper_bound(kint64max); - Coefficient offset(0); - std::deque repository; - std::vector nodes = - CreateInitialEncodingNodes(literals, coeffs, &offset, &repository); - - // Print the number of variables with a non-zero cost. - SatSolver* sat_solver = model->GetOrCreate(); - if (log_info) { - LOG(INFO) << StringPrintf("c #weights:%zu #vars:%d", nodes.size(), - sat_solver->NumVariables()); - } - - // This is used by the "stratified" approach. - Coefficient stratified_lower_bound(0); - if (sat_solver->parameters().max_sat_stratification() == - SatParameters::STRATIFICATION_DESCENT) { - // In this case, we initialize it to the maximum assumption weights. - for (EncodingNode* n : nodes) { - stratified_lower_bound = std::max(stratified_lower_bound, n->weight()); - } - } - - // Start the algorithm. - int max_depth = 0; - std::string previous_core_info = ""; - SatSolver::Status result = SatSolver::MODEL_SAT; // Not important. - for (int iter = 0;; ++iter) { - const std::vector assumptions = ReduceNodesAndExtractAssumptions( - upper_bound, stratified_lower_bound, &lower_bound, &nodes, sat_solver); - - // No assumptions with the current stratified_lower_bound, lower it if - // possible. - if (assumptions.empty()) { - stratified_lower_bound = - MaxNodeWeightSmallerThan(nodes, stratified_lower_bound); - if (stratified_lower_bound > 0) { - --iter; - continue; - } - } - - // Display the progress. - const std::string gap_string = - (upper_bound == kCoefficientMax) - ? "" - : StringPrintf(" gap:%lld", (upper_bound - lower_bound).value()); - if (log_info) { - LOG(INFO) << StringPrintf( - "c iter:%d [%s] lb:%lld%s assumptions:%zu depth:%d", iter, - previous_core_info.c_str(), lower_bound.value() - offset.value(), - gap_string.c_str(), nodes.size(), max_depth); - } - - // No assumptions means that there is no solution with cost < upper_bound. - if (assumptions.empty()) { - if (log_info) LOG(INFO) << "c no assumptions."; - result = (lower_bound == upper_bound) ? SatSolver::MODEL_SAT - : SatSolver::MODEL_UNSAT; - break; - } - - // Solve under the assumptions. - result = - SolveIntegerProblemWithLazyEncoding(assumptions, next_decision, model); - if (result == SatSolver::MODEL_SAT) { - // Extract the new solution and save it if it is the best found so far. - Coefficient obj(0); - for (int i = 0; i < literals.size(); ++i) { - if (sat_solver->Assignment().LiteralIsTrue(literals[i])) { - obj += coeffs[i]; - } - } - if (obj + offset < upper_bound) { - if (feasible_solution_observer != nullptr) { - feasible_solution_observer(*model); - } - upper_bound = obj + offset; - if (log_info) LOG(INFO) << "c ub:" << upper_bound; - } - - // If not all assumptions where taken, continue with a lower stratified - // bound. Otherwise we have an optimal solution. - stratified_lower_bound = - MaxNodeWeightSmallerThan(nodes, stratified_lower_bound); - if (stratified_lower_bound == 0) break; - } - if (result != SatSolver::ASSUMPTIONS_UNSAT) break; - - // We have a new core. - std::vector core = sat_solver->GetLastIncompatibleDecisions(); - if (sat_solver->parameters().minimize_core()) { - MinimizeCore(sat_solver, &core); - } - - // Compute the min weight of all the nodes in the core. - // The lower bound will be increased by that much. - const Coefficient min_weight = ComputeCoreMinWeight(nodes, core); - previous_core_info = - StringPrintf("core:%zu mw:%lld", core.size(), min_weight.value()); - - // Increase stratified_lower_bound according to the parameters. - if (stratified_lower_bound < min_weight && - sat_solver->parameters().max_sat_stratification() == - SatParameters::STRATIFICATION_ASCENT) { - stratified_lower_bound = min_weight; - } - - ProcessCore(core, min_weight, &repository, &nodes, sat_solver); - max_depth = std::max(max_depth, nodes.back()->depth()); - } - - if (log_info) { - LogSolveInfo(result, *sat_solver, wall_timer, user_timer, - upper_bound.value(), lower_bound.value() - offset.value()); - } - return result; -} - } // namespace sat } // namespace operations_research diff --git a/ortools/sat/optimization.h b/ortools/sat/optimization.h index 27ead2ad82..8fcc8775f1 100644 --- a/ortools/sat/optimization.h +++ b/ortools/sat/optimization.h @@ -132,12 +132,9 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( Model* model); // Same as MinimizeIntegerVariableWithLinearScanAndLazyEncoding() but use -// a core-based approach instead. The given objective_var must be equal to the -// sum of the given variables using the given coefficients. -// -// TODO(user): It is not needed to have objective_var and the linear objective -// constraint encoded in the model. Remove this precondition in order to -// improve the solving time. +// a core-based approach instead. Note that the given objective_var is just used +// for reporting the lower-bound and do not need to be linked with its linear +// representation. SatSolver::Status MinimizeWithCoreAndLazyEncoding( bool log_info, IntegerVariable objective_var, const std::vector& variables, @@ -170,19 +167,6 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( Model* model); #endif // defined(USE_CBC) || defined(USE_SCIP) -// Similar to MinimizeIntegerVariableWithLinearScanAndLazyEncoding() but use -// a core based approach. Note that this require the objective to be given as -// a weighted sum of literals -// -// TODO(user): The function above is more general, remove this one after -// checking that the performances are similar. -SatSolver::Status MinimizeWeightedLiteralSumWithCoreAndLazyEncoding( - bool log_info, const std::vector& literals, - const std::vector& coeffs, - const std::function& next_decision, - const std::function& feasible_solution_observer, - Model* model); - } // namespace sat } // namespace operations_research