From 0d8a6aeb39c69ac8b3d10900c2703bfb7a89b67b Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Thu, 18 Apr 2019 16:34:36 +0200 Subject: [PATCH] improve robustness of GLOP; continue rewriting of CP-SAT search internals; fix table extraction --- ortools/lp_data/matrix_utils.cc | 25 ++++---- ortools/sat/cp_model_presolve.cc | 5 +- ortools/sat/optimization.cc | 102 +++++++++---------------------- ortools/sat/optimization.h | 8 ++- ortools/sat/synchronization.cc | 10 +++ ortools/sat/table.cc | 13 +++- 6 files changed, 74 insertions(+), 89 deletions(-) diff --git a/ortools/lp_data/matrix_utils.cc b/ortools/lp_data/matrix_utils.cc index 898cf9e0a3..01a2bf539b 100644 --- a/ortools/lp_data/matrix_utils.cc +++ b/ortools/lp_data/matrix_utils.cc @@ -39,13 +39,13 @@ bool AreColumnsProportional(const SparseColumn& a, const SparseColumn& b, const Fractional coeff_a = a.EntryCoefficient(i); const Fractional coeff_b = b.EntryCoefficient(i); if (multiple == 0.0) { - a_is_larger = fabs(coeff_a) > fabs(coeff_b); + a_is_larger = std::abs(coeff_a) > std::abs(coeff_b); multiple = a_is_larger ? coeff_a / coeff_b : coeff_b / coeff_a; } else { if (a_is_larger) { - if (fabs(coeff_a / coeff_b - multiple) > tolerance) return false; + if (std::abs(coeff_a / coeff_b - multiple) > tolerance) return false; } else { - if (fabs(coeff_b / coeff_a - multiple) > tolerance) return false; + if (std::abs(coeff_b / coeff_a - multiple) > tolerance) return false; } } } @@ -54,10 +54,10 @@ bool AreColumnsProportional(const SparseColumn& a, const SparseColumn& b, // A column index together with its fingerprint. See ComputeFingerprint(). struct ColumnFingerprint { - ColumnFingerprint(ColIndex _col, int32 _hash, double _value) + ColumnFingerprint(ColIndex _col, int64 _hash, double _value) : col(_col), hash(_hash), value(_value) {} ColIndex col; - int32 hash; + int64 hash; double value; // This order has the property that if AreProportionalCandidates() is true for @@ -78,7 +78,7 @@ struct ColumnFingerprint { bool AreProportionalCandidates(ColumnFingerprint a, ColumnFingerprint b, Fractional tolerance) { if (a.hash != b.hash) return false; - return fabs(a.value - b.value) < tolerance; + return std::abs(a.value - b.value) < tolerance; } // The fingerprint of a column has two parts: @@ -86,16 +86,16 @@ bool AreProportionalCandidates(ColumnFingerprint a, ColumnFingerprint b, // - A double value which should be the same for two proportional columns // modulo numerical errors. ColumnFingerprint ComputeFingerprint(ColIndex col, const SparseColumn& column) { - int32 non_zero_pattern_hash = 0; + int64 non_zero_pattern_hash = 0; Fractional min_abs = std::numeric_limits::max(); Fractional max_abs = 0.0; Fractional sum = 0.0; for (const SparseColumn::Entry e : column) { non_zero_pattern_hash = - Hash32NumWithSeed(e.row().value(), non_zero_pattern_hash); + util_hash::Hash(e.row().value(), non_zero_pattern_hash); sum += e.coefficient(); - min_abs = std::min(min_abs, fabs(e.coefficient())); - max_abs = std::max(max_abs, fabs(e.coefficient())); + min_abs = std::min(min_abs, std::abs(e.coefficient())); + max_abs = std::max(max_abs, std::abs(e.coefficient())); } // The two scaled values are in [0, 1]. @@ -104,7 +104,8 @@ ColumnFingerprint ComputeFingerprint(ColIndex col, const SparseColumn& column) { DCHECK_NE(0.0, max_abs); const double inverse_dynamic_range = min_abs / max_abs; const double scaled_average = - fabs(sum) / (static_cast(column.num_entries().value()) * max_abs); + std::abs(sum) / + (static_cast(column.num_entries().value()) * max_abs); return ColumnFingerprint(col, non_zero_pattern_hash, inverse_dynamic_range + scaled_average); } @@ -172,8 +173,10 @@ ColMapping FindProportionalColumnsUsingSimpleAlgorithm( const ColIndex num_cols = matrix.num_cols(); ColMapping mapping(num_cols, kInvalidCol); for (ColIndex col_a(0); col_a < num_cols; ++col_a) { + if (matrix.column(col_a).IsEmpty()) continue; if (mapping[col_a] != kInvalidCol) continue; for (ColIndex col_b(col_a + 1); col_b < num_cols; ++col_b) { + if (matrix.column(col_b).IsEmpty()) continue; if (mapping[col_b] != kInvalidCol) continue; if (AreColumnsProportional(matrix.column(col_a), matrix.column(col_b), tolerance)) { diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index b6a7853227..fb1ebce0b2 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -3947,8 +3947,9 @@ void RemoveUnusedEquivalentVariables(PresolveContext* context) { if (domain_modified) { LOG(WARNING) << "Domain of " << r.representative << " was not fully propagated using the affine relation " - << "(representative =" << r.representative - << ", coeff = " << r.coeff << ", offset = " << r.offset + << "(var = " << var + << " representative = " << r.representative + << " coeff = " << r.coeff << " offset = " << r.offset << ")"; } proto = context->mapping_model; diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc index 12455f69bd..766bbe3c7c 100644 --- a/ortools/sat/optimization.cc +++ b/ortools/sat/optimization.cc @@ -1085,29 +1085,6 @@ SatSolver::Status SolveWithCardinalityEncodingAndCore( } } -namespace { -void LogSolveInfo(SatSolver::Status result, const SatSolver& sat_solver, - const WallTimer& wall_timer, const UserTimer& user_timer, - int64 objective, int64 best_bound) { - printf("status: %s\n", result == SatSolver::FEASIBLE - ? "OPTIMAL" - : SatStatusString(result).c_str()); - if (objective < kint64max) { - absl::PrintF("objective: %d\n", objective); - } else { - printf("objective: NA\n"); - } - absl::PrintF("best_bound: %d\n", best_bound); - printf("booleans: %d\n", sat_solver.NumVariables()); - absl::PrintF("conflicts: %d\n", sat_solver.num_failures()); - absl::PrintF("branches: %d\n", sat_solver.num_branches()); - absl::PrintF("propagations: %d\n", sat_solver.num_propagations()); - printf("walltime: %f\n", wall_timer.Get()); - printf("usertime: %f\n", user_timer.Get()); - printf("deterministic_time: %f\n", sat_solver.deterministic_time()); -} -} // namespace - SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( IntegerVariable objective_var, const std::function& feasible_solution_observer, Model* model) { @@ -1119,13 +1096,12 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( // Simple linear scan algorithm to find the optimal. SatSolver::Status result; bool model_is_feasible = false; - IntegerValue objective(kint64max); while (true) { result = ResetAndSolveIntegerProblem(/*assumptions=*/{}, model); if (result != SatSolver::FEASIBLE) break; // The objective is the current lower bound of the objective_var. - objective = integer_trail->LowerBound(objective_var); + const IntegerValue objective = integer_trail->LowerBound(objective_var); // We have a solution! model_is_feasible = true; @@ -1308,8 +1284,6 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( // This will be called each time a feasible solution is found. Returns false // if a conflict was detected while trying to constrain the objective to a // smaller value. - int num_solutions = 0; - IntegerValue best_objective = integer_trail->UpperBound(objective_var); const auto process_solution = [&]() { // We don't assume that objective_var is linked with its linear term, so // we recompute the objective here. @@ -1318,14 +1292,8 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( objective += coefficients[i] * IntegerValue(model->Get(LowerBound(variables[i]))); } + if (objective > integer_trail->UpperBound(objective_var)) return true; - // The situation is tricky when best_objective is an externally provided - // bound, but we don't have any solution lower or equal to this bound yet. - if (objective > best_objective) return true; - if (objective >= best_objective && num_solutions > 0) return true; - - ++num_solutions; - best_objective = objective; if (feasible_solution_observer != nullptr) { feasible_solution_observer(); } @@ -1427,8 +1395,7 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( if (!integer_trail->Enqueue( IntegerLiteral::GreaterOrEqual(objective_var, implied_objective_lb), {}, {})) { - result = SatSolver::INFEASIBLE; - break; + return SatSolver::INFEASIBLE; } // No assumptions with the current stratified_threshold? use the new one. @@ -1453,22 +1420,20 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( model->Add(WeightedSumLowerOrEqual(constraint_vars, constraint_coeffs, -objective_offset.value())); - result = MinimizeIntegerVariableWithLinearScanAndLazyEncoding( + return MinimizeIntegerVariableWithLinearScanAndLazyEncoding( objective_var, feasible_solution_observer, model); - break; } // Display the progress. - const IntegerValue objective_lb = integer_trail->LowerBound(objective_var); if (VLOG_IS_ON(1)) { - const int64 lb = objective_lb.value(); - const int64 ub = best_objective.value(); + const int64 lb = integer_trail->LowerBound(objective_var).value(); + const int64 ub = integer_trail->UpperBound(objective_var).value(); const int gap = lb == ub ? 0 : static_cast(std::ceil( 100.0 * (ub - lb) / std::max(std::abs(ub), std::abs(lb)))); - VLOG(1) << absl::StrCat("unscaled_objective:[", lb, ",", ub, + VLOG(1) << absl::StrCat("unscaled_next_obj_range:[", lb, ",", ub, "]" " gap:", gap, "%", " assumptions:", term_indices.size(), @@ -1476,12 +1441,6 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( " depth:", max_depth); } - // Abort if we have a solution and a gap of zero. - if (objective_lb == best_objective && num_solutions > 0) { - result = SatSolver::FEASIBLE; - break; - } - // Bulk cover optimization. bool some_cover_opt = false; if (parameters.cover_optimization()) { @@ -1528,10 +1487,7 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( VLOG(1) << "cover_opt var:" << var << " domain:[" << integer_trail->LevelZeroLowerBound(var) << "," << best << "]"; - if (!process_solution()) { - result = SatSolver::INFEASIBLE; - break; - } + if (!process_solution()) return SatSolver::INFEASIBLE; if (parameters.stop_after_first_solution()) { return SatSolver::LIMIT_REACHED; } @@ -1544,8 +1500,7 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( // order to abort early if the optimal is proved. if (!integer_trail->Enqueue(IntegerLiteral::GreaterOrEqual(var, best), {}, {})) { - result = SatSolver::INFEASIBLE; - break; + return SatSolver::INFEASIBLE; } } else if (result != SatSolver::FEASIBLE) { break; @@ -1573,7 +1528,10 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( std::vector> cores; result = FindCores(assumptions, model, &cores); if (result == SatSolver::FEASIBLE) { - process_solution(); + if (!process_solution()) return SatSolver::INFEASIBLE; + + // TODO(user): We actually do not know if this solution was not out of + // the requested objective bound. if (parameters.stop_after_first_solution()) { return SatSolver::LIMIT_REACHED; } @@ -1670,10 +1628,7 @@ SatSolver::Status MinimizeWithCoreAndLazyEncoding( } } - // Returns FEASIBLE if we found the optimal. - return num_solutions > 0 && result == SatSolver::INFEASIBLE - ? SatSolver::FEASIBLE - : result; + return result; } // TODO(user): take the MPModelRequest or MPModelProto directly, so that we can @@ -1690,8 +1645,6 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( IntegerEncoder* integer_encoder = model->GetOrCreate(); // This will be called each time a feasible solution is found. - int num_solutions = 0; - IntegerValue best_objective = integer_trail->UpperBound(objective_var); const auto process_solution = [&]() { // We don't assume that objective_var is linked with its linear term, so // we recompute the objective here. @@ -1700,12 +1653,22 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( objective += coefficients[i] * IntegerValue(model->Get(Value(variables[i]))); } - if (objective >= best_objective) return; - ++num_solutions; - best_objective = objective; + if (objective > integer_trail->UpperBound(objective_var)) return true; + if (feasible_solution_observer != nullptr) { feasible_solution_observer(); } + + // Constrain objective_var. This has a better result when objective_var is + // used in an LP relaxation for instance. + sat_solver->Backtrack(0); + sat_solver->SetAssumptionLevel(0); + if (!integer_trail->Enqueue( + IntegerLiteral::LowerOrEqual(objective_var, objective - 1), {}, + {})) { + return false; + } + return true; }; // This is the "generalized" hitting set problem we will solve. Each time @@ -1812,9 +1775,7 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( --iter; // "false" iteration, the lower bound does not increase. continue; } else { - result = - num_solutions > 0 ? SatSolver::FEASIBLE : SatSolver::INFEASIBLE; - break; + return SatSolver::INFEASIBLE; } } @@ -1823,7 +1784,7 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( std::vector> cores; result = FindCores(assumptions, model, &cores); if (result == SatSolver::FEASIBLE) { - process_solution(); + if (!process_solution()) return SatSolver::INFEASIBLE; if (parameters.stop_after_first_solution()) { return SatSolver::LIMIT_REACHED; } @@ -1888,10 +1849,7 @@ SatSolver::Status MinimizeWithHittingSetAndLazyEncoding( } } - // Returns FEASIBLE if we found the optimal. - return num_solutions > 0 && result == SatSolver::INFEASIBLE - ? SatSolver::FEASIBLE - : result; + return result; #else // !__PORTABLE_PLATFORM__ && (USE_CBC || USE_SCIP) LOG(FATAL) << "Not supported."; #endif // __PORTABLE_PLATFORM || (!USE_CBC && !USE_SCIP) diff --git a/ortools/sat/optimization.h b/ortools/sat/optimization.h index 03059dd816..aa46886f58 100644 --- a/ortools/sat/optimization.h +++ b/ortools/sat/optimization.h @@ -143,8 +143,12 @@ void RestrictObjectiveDomainWithBinarySearch( // Same as MinimizeIntegerVariableWithLinearScanAndLazyEncoding() but use // 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. +// for reporting the lower-bound/upper-bound and do not need to be linked with +// its linear representation. +// +// Unlike MinimizeIntegerVariableWithLinearScanAndLazyEncoding() this function +// just return the last solver status. In particular if it is INFEASIBLE but +// feasible_solution_observer() was called, it means we are at OPTIMAL. SatSolver::Status MinimizeWithCoreAndLazyEncoding( IntegerVariable objective_var, const std::vector& variables, diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index 0281786901..96cd7c031b 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -47,6 +47,16 @@ void SharedResponseManager::UpdateInnerObjectiveBounds( change = true; inner_objective_upper_bound_ = ub.value(); } + if (inner_objective_lower_bound_ > inner_objective_upper_bound_) { + if (best_response_.status() == CpSolverStatus::FEASIBLE || + best_response_.status() == CpSolverStatus::OPTIMAL) { + best_response_.set_status(CpSolverStatus::OPTIMAL); + } else { + best_response_.set_status(CpSolverStatus::INFEASIBLE); + } + if (log_updates_) LogNewSatSolution("Done", wall_timer_.Get(), worker_info); + return; + } if (log_updates_ && change) { const CpObjectiveProto& obj = model_proto_.objective(); double new_lb = ScaleObjectiveValue(obj, inner_objective_lower_bound_); diff --git a/ortools/sat/table.cc b/ortools/sat/table.cc index 83286541ef..13ca5fa444 100644 --- a/ortools/sat/table.cc +++ b/ortools/sat/table.cc @@ -219,6 +219,10 @@ void AddTableConstraint(const std::vector& vars, } // The variable domains have been computed. Fully encode variables. + // Note that in some corner cases (like duplicate vars), as we call + // UpdateInitialDomain(), the domain of other variable could become more + // restricted that values_per_var. For now, we do not try to reach a fixed + // point here. std::vector> encodings(n); for (int i = 0; i < n; ++i) { const std::vector reached_values(values_per_var[i].begin(), @@ -293,6 +297,7 @@ void AddTableConstraint(const std::vector& vars, std::vector clause; for (int j = 0; j < tuples.size(); ++j) { clause.clear(); + bool tuple_is_valid = true; for (int i = 0; i + 1 < n; ++i) { // Ignore fixed variables. if (values_per_var[i].size() == 1) continue; @@ -302,13 +307,17 @@ void AddTableConstraint(const std::vector& vars, if (v == any_value) continue; const IntegerValue value(v); - CHECK(encodings[i].contains(value)); + if (!encodings[i].contains(value)) { + tuple_is_valid = false; + break; + } clause.push_back(gtl::FindOrDie(encodings[i], value).Negated()); } + if (!tuple_is_valid) continue; // Add the target of the implication. const IntegerValue target_value = IntegerValue(tuples[j][n - 1]); - CHECK(encodings[n - 1].contains(target_value)); + if (!encodings[n - 1].contains(target_value)) continue; const Literal target_literal = gtl::FindOrDie(encodings[n - 1], target_value); clause.push_back(target_literal);