diff --git a/makefiles/Makefile.cpp.mk b/makefiles/Makefile.cpp.mk index 6b66ecb071..9e9b4a2185 100644 --- a/makefiles/Makefile.cpp.mk +++ b/makefiles/Makefile.cpp.mk @@ -474,7 +474,9 @@ LINEAR_SOLVER_LIB_OBJS = \ $(OBJ_DIR)/gurobi_interface.$O \ $(OBJ_DIR)/linear_solver.$O \ $(OBJ_DIR)/linear_solver.pb.$O \ + $(OBJ_DIR)/linear_solver2.pb.$O \ $(OBJ_DIR)/model_exporter.$O \ + $(OBJ_DIR)/proto_tools.$O \ $(OBJ_DIR)/scip_interface.$O \ $(OBJ_DIR)/sulum_interface.$O @@ -491,7 +493,7 @@ $(OBJ_DIR)/glpk_interface.$O:$(SRC_DIR)/linear_solver/glpk_interface.cc $(OBJ_DIR)/gurobi_interface.$O:$(SRC_DIR)/linear_solver/gurobi_interface.cc $(CCC) $(CFLAGS) -c $(SRC_DIR)/linear_solver/gurobi_interface.cc $(OBJ_OUT)$(OBJ_DIR)$Sgurobi_interface.$O -$(OBJ_DIR)/linear_solver.$O:$(SRC_DIR)/linear_solver/linear_solver.cc $(GEN_DIR)/linear_solver/linear_solver.pb.h +$(OBJ_DIR)/linear_solver.$O:$(SRC_DIR)/linear_solver/linear_solver.cc $(GEN_DIR)/linear_solver/linear_solver.pb.h $(GEN_DIR)/linear_solver/linear_solver2.pb.h $(CCC) $(CFLAGS) -c $(SRC_DIR)/linear_solver/linear_solver.cc $(OBJ_OUT)$(OBJ_DIR)$Slinear_solver.$O $(OBJ_DIR)/linear_solver.pb.$O:$(GEN_DIR)/linear_solver/linear_solver.pb.cc @@ -502,9 +504,20 @@ $(GEN_DIR)/linear_solver/linear_solver.pb.cc:$(SRC_DIR)/linear_solver/linear_sol $(GEN_DIR)/linear_solver/linear_solver.pb.h:$(GEN_DIR)/linear_solver/linear_solver.pb.cc +$(OBJ_DIR)/linear_solver2.pb.$O:$(GEN_DIR)/linear_solver/linear_solver2.pb.cc + $(CCC) $(CFLAGS) -c $(GEN_DIR)/linear_solver/linear_solver2.pb.cc $(OBJ_OUT)$(OBJ_DIR)$Slinear_solver2.pb.$O + +$(GEN_DIR)/linear_solver/linear_solver2.pb.cc:$(SRC_DIR)/linear_solver/linear_solver2.proto + $(PROTOBUF_DIR)/bin/protoc --proto_path=$(INC_DIR) --cpp_out=$(GEN_DIR) $(SRC_DIR)/linear_solver/linear_solver2.proto + +$(GEN_DIR)/linear_solver/linear_solver2.pb.h:$(GEN_DIR)/linear_solver/linear_solver2.pb.cc + $(OBJ_DIR)/model_exporter.$O:$(SRC_DIR)/linear_solver/model_exporter.cc $(GEN_DIR)/linear_solver/linear_solver.pb.h $(CCC) $(CFLAGS) -c $(SRC_DIR)/linear_solver/model_exporter.cc $(OBJ_OUT)$(OBJ_DIR)$Smodel_exporter.$O +$(OBJ_DIR)/proto_tools.$O:$(SRC_DIR)/linear_solver/proto_tools.cc $(GEN_DIR)/linear_solver/linear_solver.pb.h $(GEN_DIR)/linear_solver/linear_solver2.pb.h + $(CCC) $(CFLAGS) -c $(SRC_DIR)/linear_solver/proto_tools.cc $(OBJ_OUT)$(OBJ_DIR)$Sproto_tools.$O + $(OBJ_DIR)/scip_interface.$O:$(SRC_DIR)/linear_solver/scip_interface.cc $(CCC) $(CFLAGS) -c $(SRC_DIR)/linear_solver/scip_interface.cc $(OBJ_OUT)$(OBJ_DIR)$Sscip_interface.$O diff --git a/src/linear_solver/glpk_interface.cc b/src/linear_solver/glpk_interface.cc index 08f3d1c563..a07c07d5c5 100644 --- a/src/linear_solver/glpk_interface.cc +++ b/src/linear_solver/glpk_interface.cc @@ -35,18 +35,17 @@ #if defined(USE_GLPK) extern "C" { - #include "glpk.h" +#include "glpk.h" } DECLARE_double(solver_timeout_in_seconds); DECLARE_string(solver_write_model); namespace operations_research { - // Class to store information gathered in the callback class GLPKInformation { public: - explicit GLPKInformation(bool maximize): num_all_nodes_(0) { + explicit GLPKInformation(bool maximize) : num_all_nodes_(0) { ResetBestObjectiveBound(maximize); } void Reset(bool maximize) { @@ -64,8 +63,6 @@ class GLPKInformation { double best_objective_bound_; }; - - // Function to be called in the GLPK callback void GLPKGatherInformationCallback(glp_tree* tree, void* info) { CHECK_NOTNULL(tree); @@ -122,8 +119,7 @@ class GLPKInterface : public MPSolverInterface { // Change a coefficient in a constraint. virtual void SetCoefficient(MPConstraint* const constraint, const MPVariable* const variable, - double new_value, - double old_value); + double new_value, double old_value); // Clear a constraint from all its terms. virtual void ClearConstraint(MPConstraint* const constraint); // Change a coefficient in the linear objective @@ -166,9 +162,7 @@ class GLPKInterface : public MPSolverInterface { return StringPrintf("GLPK %s", glp_version()); } - virtual void* underlying_solver() { - return reinterpret_cast(lp_); - } + virtual void* underlying_solver() { return reinterpret_cast(lp_); } virtual double ComputeExactConditionNumber() const; @@ -187,8 +181,7 @@ class GLPKInterface : public MPSolverInterface { virtual void SetLpAlgorithm(int value); void ExtractOldConstraints(); - void ExtractOneConstraint(MPConstraint* const constraint, - int* const indices, + void ExtractOneConstraint(MPConstraint* const constraint, int* const indices, double* const coefs); // Transforms basis status from GLPK integer code to MPSolver::BasisStatus. MPSolver::BasisStatus TransformGLPKBasisStatus(int glpk_basis_status) const; @@ -196,16 +189,16 @@ class GLPKInterface : public MPSolverInterface { // Computes the L1-norm of the current scaled basis. // The L1-norm |A| is defined as max_j sum_i |a_ij| // This method is available only for continuous problems. - double ComputeScaledBasisL1Norm( - int num_rows, int num_cols, - double* row_scaling_factor, double* column_scaling_factor) const; + double ComputeScaledBasisL1Norm(int num_rows, int num_cols, + double* row_scaling_factor, + double* column_scaling_factor) const; // Computes the L1-norm of the inverse of the current scaled // basis. // This method is available only for continuous problems. - double ComputeInverseScaledBasisL1Norm( - int num_rows, int num_cols, - double* row_scaling_factor, double* column_scaling_factor) const; + double ComputeInverseScaledBasisL1Norm(int num_rows, int num_cols, + double* row_scaling_factor, + double* column_scaling_factor) const; glp_prob* lp_; bool mip_; @@ -281,7 +274,7 @@ void GLPKInterface::SetVariableInteger(int var_index, bool integer) { if (mip_) { if (var_index != kNoIndex) { // Not cached if the variable has been extracted. - glp_set_col_kind(lp_, var_index, integer ? GLP_IV : GLP_CV); + glp_set_col_kind(lp_, var_index, integer ? GLP_IV : GLP_CV); } else { sync_status_ = MUST_RELOAD; } @@ -316,8 +309,7 @@ void GLPKInterface::SetConstraintBounds(int index, double lb, double ub) { void GLPKInterface::SetCoefficient(MPConstraint* const constraint, const MPVariable* const variable, - double new_value, - double old_value) { + double new_value, double old_value) { InvalidateSolutionSynchronization(); // GLPK does not allow to modify one coefficient at a time, so we // extract the whole constraint again, if it has been extracted @@ -405,15 +397,15 @@ void GLPKInterface::ExtractNewVariables() { // Extract again existing constraints if they contain new variables. void GLPKInterface::ExtractOldConstraints() { - int max_constraint_size = solver_->ComputeMaxConstraintSize( - 0, last_constraint_index_); + int max_constraint_size = + solver_->ComputeMaxConstraintSize(0, last_constraint_index_); // The first entry in the following arrays is dummy, to be // consistent with glpk API. std::unique_ptr indices(new int[max_constraint_size + 1]); std::unique_ptr coefs(new double[max_constraint_size + 1]); for (int i = 0; i < last_constraint_index_; ++i) { - MPConstraint* const ct = solver_->constraints_[i]; + MPConstraint* const ct = solver_->constraints_[i]; DCHECK_NE(kNoIndex, ct->index()); const int size = ct->coefficients_.size(); if (size == 0) { @@ -441,7 +433,7 @@ void GLPKInterface::ExtractOneConstraint(MPConstraint* const constraint, coefs[k] = entry.second; ++k; } - glp_set_mat_row(lp_, constraint->index(), k-1, indices, coefs); + glp_set_mat_row(lp_, constraint->index(), k - 1, indices, coefs); } // Define new constraints on old and new variables. @@ -595,15 +587,13 @@ MPSolver::ResultStatus GLPKInterface::Solve(const MPSolverParameters& param) { if (mip_) { const double row_activity = glp_mip_row_val(lp_, ct->index()); ct->set_activity(row_activity); - VLOG(4) << "row " << ct->index() - << ": activity = " << row_activity; + VLOG(4) << "row " << ct->index() << ": activity = " << row_activity; } else { const double row_activity = glp_get_row_prim(lp_, ct->index()); ct->set_activity(row_activity); const double dual_value = glp_get_row_dual(lp_, ct->index()); ct->set_dual_value(dual_value); - VLOG(4) << "row " << ct->index() - << ": activity = " << row_activity + VLOG(4) << "row " << ct->index() << ": activity = " << row_activity << ": dual value = " << dual_value; } } @@ -633,8 +623,7 @@ MPSolver::ResultStatus GLPKInterface::Solve(const MPSolverParameters& param) { result_status_ = MPSolver::OPTIMAL; } else if (tmp_status == GLP_FEAS) { result_status_ = MPSolver::FEASIBLE; - } else if (tmp_status == GLP_NOFEAS || - tmp_status == GLP_INFEAS) { + } else if (tmp_status == GLP_NOFEAS || tmp_status == GLP_INFEAS) { // For infeasible problems, GLPK actually seems to return // GLP_UNDEF. So this is never (?) reached. Return infeasible // in case GLPK returns a correct status in future versions. @@ -654,8 +643,8 @@ MPSolver::ResultStatus GLPKInterface::Solve(const MPSolverParameters& param) { return result_status_; } -MPSolver::BasisStatus -GLPKInterface::TransformGLPKBasisStatus(int glpk_basis_status) const { +MPSolver::BasisStatus GLPKInterface::TransformGLPKBasisStatus( + int glpk_basis_status) const { switch (glpk_basis_status) { case GLP_BS: return MPSolver::BASIC; @@ -676,17 +665,12 @@ GLPKInterface::TransformGLPKBasisStatus(int glpk_basis_status) const { // ------ Query statistics on the solution and the solve ------ int64 GLPKInterface::iterations() const { - if (mip_) { - LOG(WARNING) << "Total number of iterations is not available"; - return kUnknownNumberOfIterations; - } else { - if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfIterations; -#if GLP_MINOR_VERSION >= 49 - return kUnknownNumberOfIterations; -#else +#if GLP_MINOR_VERSION < 49 + if (!mip_ && CheckSolutionIsSynchronized()) return lpx_get_int_parm(lp_, LPX_K_ITCNT); #endif - } + LOG(WARNING) << "Total number of iterations is not available"; + return kUnknownNumberOfIterations; } int64 GLPKInterface::nodes() const { @@ -732,7 +716,6 @@ MPSolver::BasisStatus GLPKInterface::column_status(int variable_index) const { return TransformGLPKBasisStatus(glpk_basis_status); } - bool GLPKInterface::CheckSolutionExists() const { if (result_status_ == MPSolver::ABNORMAL) { LOG(WARNING) << "Ignoring ABNORMAL status from GLPK: This status may or may" @@ -759,9 +742,8 @@ bool GLPKInterface::CheckBestObjectiveBoundExists() const { double GLPKInterface::ComputeExactConditionNumber() const { if (!IsContinuous()) { // TODO(user): support MIP. - LOG(DFATAL) - << "ComputeExactConditionNumber not implemented for" - << " GLPK_MIXED_INTEGER_PROGRAMMING"; + LOG(DFATAL) << "ComputeExactConditionNumber not implemented for" + << " GLPK_MIXED_INTEGER_PROGRAMMING"; return 0.0; } if (!CheckSolutionIsSynchronized()) return 0.0; @@ -779,18 +761,16 @@ double GLPKInterface::ComputeExactConditionNumber() const { for (int col = 1; col <= num_cols; ++col) { column_scaling_factor[col] = glp_get_sjj(lp_, col); } - return - ComputeInverseScaledBasisL1Norm( - num_rows, num_cols, - row_scaling_factor.get(), column_scaling_factor.get()) * - ComputeScaledBasisL1Norm( - num_rows, num_cols, - row_scaling_factor.get(), column_scaling_factor.get()); + return ComputeInverseScaledBasisL1Norm(num_rows, num_cols, + row_scaling_factor.get(), + column_scaling_factor.get()) * + ComputeScaledBasisL1Norm(num_rows, num_cols, row_scaling_factor.get(), + column_scaling_factor.get()); } double GLPKInterface::ComputeScaledBasisL1Norm( - int num_rows, int num_cols, - double* row_scaling_factor, double* column_scaling_factor) const { + int num_rows, int num_cols, double* row_scaling_factor, + double* column_scaling_factor) const { double norm = 0.0; std::unique_ptr values(new double[num_rows + 1]); std::unique_ptr indices(new int[num_rows + 1]); @@ -826,8 +806,8 @@ double GLPKInterface::ComputeScaledBasisL1Norm( } double GLPKInterface::ComputeInverseScaledBasisL1Norm( - int num_rows, int num_cols, - double* row_scaling_factor, double* column_scaling_factor) const { + int num_rows, int num_cols, double* row_scaling_factor, + double* column_scaling_factor) const { // Compute the LU factorization if it doesn't exist yet. if (!glp_bf_exists(lp_)) { const int factorize_status = glp_factorize(lp_); @@ -929,7 +909,7 @@ void GLPKInterface::ConfigureGLPKParameters(const MPSolverParameters& param) { glp_scale_prob(lp_, GLP_SF_AUTO); // Use advanced initial basis (options: standard / advanced / Bixby's). - glp_adv_basis(lp_, nullptr); + glp_adv_basis(lp_, 0); // Set parameters specified by the user. SetParameters(param); @@ -955,9 +935,7 @@ void GLPKInterface::SetPrimalTolerance(double value) { lp_param_.tol_bnd = value; } -void GLPKInterface::SetDualTolerance(double value) { - lp_param_.tol_dj = value; -} +void GLPKInterface::SetDualTolerance(double value) { lp_param_.tol_dj = value; } void GLPKInterface::SetPresolveMode(int value) { switch (value) { @@ -1004,6 +982,5 @@ MPSolverInterface* BuildGLPKInterface(bool mip, MPSolver* const solver) { return new GLPKInterface(solver, mip); } - } // namespace operations_research #endif // #if defined(USE_GLPK) diff --git a/src/linear_solver/linear_solver.cc b/src/linear_solver/linear_solver.cc index 9d0ac31550..5b6a98c15c 100644 --- a/src/linear_solver/linear_solver.cc +++ b/src/linear_solver/linear_solver.cc @@ -35,7 +35,9 @@ #include "base/hash.h" #include "base/accurate_sum.h" #include "linear_solver/linear_solver.pb.h" +#include "linear_solver/linear_solver2.pb.h" #include "linear_solver/model_exporter.h" +#include "linear_solver/proto_tools.h" #include "util/fp_utils.h" @@ -553,6 +555,47 @@ MPSolver::LoadStatus MPSolver::LoadModel(const MPModelProto& input_model) { return MPSolver::NO_ERROR; } +MPSolver::LoadStatus MPSolver::LoadModelFromProto( + const new_proto::MPModelProto& input_model) { + for (int i = 0; i < input_model.variable_size(); ++i) { + const new_proto::MPVariableProto& var_proto = input_model.variable(i); + + // TODO(user): The variable name var_proto.name() is lost at this point. + // This is because MPSolver has no notion of name, just of ids, and these + // must be unique which is not necessarily the case of names in the new + // proto. MakeNumVar() will just assign an automated variable id when the + // the given name argument is empty. + // + // Note(user): The auto assigned id limits the number of variables to 10^9. + MPVariable* variable = MakeNumVar( + var_proto.lower_bound(), var_proto.upper_bound(), /*name=*/ ""); + variable->SetInteger(var_proto.is_integer()); + SetObjectiveCoefficient(variable, var_proto.objective_coefficient()); + } + + for (int i = 0; i < input_model.constraint_size(); ++i) { + const new_proto::MPConstraintProto& ct_proto = input_model.constraint(i); + MPConstraint* const ct = MakeRowConstraint(ct_proto.lower_bound(), + ct_proto.upper_bound(), + ct_proto.name()); + for (int j = 0; j < ct_proto.linear_term_size(); ++j) { + const new_proto::MPConstraintProto::UnaryTerm& term_proto = + ct_proto.linear_term(j); + if (term_proto.var_index() >= variables_.size()) { + LOG(ERROR) << "Variable index out of bound in constraint named " + << ct_proto.name() << "."; + return MPSolver::UNKNOWN_VARIABLE_ID; + } + ct->SetCoefficient(variables_[term_proto.var_index()], + term_proto.coefficient()); + } + } + SetOptimizationDirection(input_model.maximize()); + if (input_model.has_objective_offset()) { + MutableObjective()->SetOffset(input_model.objective_offset()); + } + return MPSolver::NO_ERROR; +} void MPSolver::FillSolutionResponse(MPSolutionResponse* response) const { CHECK_NOTNULL(response); @@ -612,6 +655,40 @@ void MPSolver::FillSolutionResponse(MPSolutionResponse* response) const { } } +namespace { +new_proto::MPSolutionResponse::Status ResultStatusToMPSolutionResponse( + MPSolver::ResultStatus status) { + switch (status) { + case MPSolver::OPTIMAL : + return new_proto::MPSolutionResponse::OPTIMAL; + case MPSolver::FEASIBLE : + return new_proto::MPSolutionResponse::FEASIBLE; + case MPSolver::INFEASIBLE : + return new_proto::MPSolutionResponse::INFEASIBLE; + case MPSolver::UNBOUNDED : + return new_proto::MPSolutionResponse::UNBOUNDED; + case MPSolver::ABNORMAL : + return new_proto::MPSolutionResponse::ABNORMAL; + default: + return new_proto::MPSolutionResponse::UNKNOWN; + } +} +} // namespace + +void MPSolver::FillSolutionResponseProto( + new_proto::MPSolutionResponse* response) const { + CHECK_NOTNULL(response); + response->Clear(); + response->set_status( + ResultStatusToMPSolutionResponse(interface_->result_status_)); + if (interface_->result_status_ == MPSolver::OPTIMAL || + interface_->result_status_ == MPSolver::FEASIBLE) { + response->set_objective_value(objective_value()); + for (int i = 0; i < variables_.size(); ++i) { + response->add_variable_value(variables_[i]->solution_value()); + } + } +} // static void MPSolver::SolveWithProtocolBuffers(const MPModelRequest& model_request, @@ -624,6 +701,8 @@ void MPSolver::SolveWithProtocolBuffers(const MPModelRequest& model_request, const MPSolver::LoadStatus loadStatus = solver.LoadModel(model); if (loadStatus != MPSolver::NO_ERROR) { LOG(WARNING) << "Loading model from protocol buffer failed, " + << "load status = " + << Error::Code_Name(static_cast(loadStatus)) << " (" << loadStatus << ")"; response->set_result_status(MPSolutionResponse::ABNORMAL); return; @@ -636,6 +715,33 @@ void MPSolver::SolveWithProtocolBuffers(const MPModelRequest& model_request, solver.FillSolutionResponse(response); } +// static +void MPSolver::SolveWithProto( + const new_proto::MPModelRequest& model_request, + new_proto::MPSolutionResponse* response) { + CHECK_NOTNULL(response); + const new_proto::MPModelProto& model = model_request.model(); + MPSolver solver(model.name(), + static_cast( + model_request.solver_type())); + const MPSolver::LoadStatus loadStatus = solver.LoadModelFromProto(model); + if (loadStatus != MPSolver::NO_ERROR) { + LOG(WARNING) << "Loading model from protocol buffer failed, " + << "load status = " + << Error::Code_Name(static_cast(loadStatus)) + << " (" << loadStatus << ")"; + response->set_status(new_proto::MPSolutionResponse::ABNORMAL); + return; + } + if (model_request.has_solver_time_limit_seconds()) { + // static_cast avoids a warning with -Wreal-conversion. This + // helps catching bugs with unwanted conversions from double to ints. + solver.set_time_limit( + static_cast(model_request.solver_time_limit_seconds()) * 1000); + } + solver.Solve(); + solver.FillSolutionResponseProto(response); +} namespace { // Outputs the terms in var_coeff_map to output_term, by sorting the @@ -673,32 +779,32 @@ void OutputTermsToProto( void MPSolver::ExportModel(MPModelProto* output_model) const { DCHECK(output_model != NULL); - if (output_model->variables_size() > 0 || - output_model->has_maximize() || - output_model->objective_terms_size() > 0 || - output_model->constraints_size() > 0 || - output_model->has_name() || - output_model->has_objective_offset()) { - LOG(WARNING) << "The model protocol buffer is not empty, " - << "it will be overwritten."; - output_model->clear_variables(); - output_model->clear_maximize(); - output_model->clear_objective_terms(); - output_model->clear_constraints(); - output_model->clear_name(); - } + new_proto::MPModelProto new_proto_model; + ExportModelToNewProto(&new_proto_model); + ConvertNewMPModelProtoToOld(new_proto_model, output_model); +} +void MPSolver::ExportModelToNewProto( + new_proto::MPModelProto* output_model) const { + DCHECK(output_model != NULL); + output_model->Clear(); // Name output_model->set_name(Name()); // Variables for (int j = 0; j < variables_.size(); ++j) { const MPVariable* const var = variables_[j]; - MPVariableProto* const variable_proto = output_model->add_variables(); - DCHECK(!var->name().empty()); - variable_proto->set_id(var->name()); - variable_proto->set_lb(var->lb()); - variable_proto->set_ub(var->ub()); - variable_proto->set_integer(var->integer()); + new_proto::MPVariableProto* const variable_proto + = output_model->add_variable(); + // TODO(user): Add option to avoid filling the var name to avoid overly + // large protocol buffers. + variable_proto->set_name(var->name()); + variable_proto->set_lower_bound(var->lb()); + variable_proto->set_upper_bound(var->ub()); + variable_proto->set_is_integer(var->integer()); + if (objective_->GetCoefficient(var) != 0.0) { + variable_proto->set_objective_coefficient( + objective_->GetCoefficient(var)); + } } // Map the variables to their indices. This is needed to output the @@ -707,37 +813,56 @@ void MPSolver::ExportModel(MPModelProto* output_model) const { // This step is needed as long as the variable indices are given by the // underlying solver at the time of model extraction. // TODO(user): remove this step. - hash_map var_name_to_index; + hash_map var_to_index; for (int j = 0; j < variables_.size(); ++j) { - var_name_to_index[variables_[j]] = j; + var_to_index[variables_[j]] = j; } // Constraints for (int i = 0; i < constraints_.size(); ++i) { MPConstraint* const constraint = constraints_[i]; - MPConstraintProto* const constraint_proto = output_model->add_constraints(); - // Constraint names need to be non-empty. - DCHECK(!constraint->name().empty()); - constraint_proto->set_id(constraint->name()); - constraint_proto->set_lb(constraint->lb()); - constraint_proto->set_ub(constraint->ub()); - OutputTermsToProto(variables_, var_name_to_index, constraint->coefficients_, - constraint_proto->mutable_terms()); + new_proto::MPConstraintProto* const constraint_proto + = output_model->add_constraint(); + constraint_proto->set_name(constraint->name()); + constraint_proto->set_lower_bound(constraint->lb()); + constraint_proto->set_upper_bound(constraint->ub()); + constraint_proto->set_is_lazy(constraint->is_lazy()); + // Vector linear_term will contain pairs (variable index, coeff), that will + // be sorted by variable index. + std::vector > linear_term; + for (CoeffEntry entry : constraint->coefficients_) { + const MPVariable* const var = entry.first; + const int var_index = FindWithDefault(var_to_index, var, -1); + DCHECK_NE(-1, var_index); + const double coeff = entry.second; + linear_term.push_back(std::pair(var_index, coeff)); + } + // The cost of std::sort is expected to be low as constraints usually have very + // few terms. + std::sort(linear_term.begin(), linear_term.end()); + // Now use linear term. + for (int k = 0; k < linear_term.size(); ++k) { + const std::pair& p = linear_term[k]; + const int var_index = p.first; + const double coeff = p.second; + new_proto::MPConstraintProto_UnaryTerm* unary_term + = constraint_proto->add_linear_term(); + unary_term->set_var_index(var_index); + unary_term->set_coefficient(coeff); + constraint_proto->add_var_index(var_index); + constraint_proto->add_coefficient(coeff); + } } - - // Objective. - OutputTermsToProto(variables_, var_name_to_index, objective_->coefficients_, - output_model->mutable_objective_terms()); output_model->set_maximize(Objective().maximization()); output_model->set_objective_offset(Objective().offset()); } - -bool MPSolver::LoadSolutionFromProto(const MPSolutionResponse& response) { +bool MPSolver::LoadSolutionFromNewProto( + const new_proto::MPSolutionResponse& response) { interface_->result_status_ = - static_cast(response.result_status()); - if (response.result_status() != MPSolutionResponse::OPTIMAL && - response.result_status() != MPSolutionResponse::FEASIBLE) { + static_cast(response.status()); + if (response.status() != new_proto::MPSolutionResponse::OPTIMAL && + response.status() != new_proto::MPSolutionResponse::FEASIBLE) { LOG(ERROR) << "Cannot load a solution unless its status is OPTIMAL or FEASIBLE."; return false; @@ -745,55 +870,37 @@ bool MPSolver::LoadSolutionFromProto(const MPSolutionResponse& response) { // Before touching the variables, verify that the solution looks legit: // each variable of the MPSolver must have its value listed exactly once, and // each listed solution should correspond to a known variable. - if (response.solution_values_size() != variables_.size()) { + if (response.variable_value_size() != variables_.size()) { LOG(ERROR) << "Trying to load a solution whose number of variables does not" << " correspond to the Solver."; return false; } - std::vector variable_has_solution_value(variables_.size(), false); + double largest_error = 0; + interface_->ExtractModel(); int num_vars_out_of_bounds = 0; - for (int i = 0; i < response.solution_values_size(); ++i) { - const string& var_name = response.solution_values(i).variable_id(); - const int var_index = - FindWithDefault(variable_name_to_index_, var_name, -1); - if (var_index < 0) { - LOG(ERROR) - << "Trying to load a solution with unknown var '" << var_name << "'."; - return false; - } - if (variable_has_solution_value[var_index]) { - LOG(ERROR) - << "Trying to load a solution value where variable '" << var_name - << "' has its solution value listed twice."; - return false; - } - variable_has_solution_value[var_index] = true; + const double tolerance = MPSolverParameters::kDefaultPrimalTolerance; + for (int i = 0; i< response.variable_value_size(); ++i) { // Look further: verify the bounds. Since linear solvers yield (small) // numerical errors, though, we just log a warning if the variables look // like they are out of their bounds. The user should inspect the values. - const double var_value = response.solution_values(i).value(); - const MPVariable* var = variables_[var_index]; - DCHECK_EQ(var->name(), var_name); - if (var_value < var->lb() || var_value > var->ub()) { + const double var_value = response.variable_value(i); + MPVariable* var = variables_[i]; + // TODO(user): Use parameter when they become available in this class. + const double lb_error = var->lb() - var_value; + const double ub_error = var_value - var->ub(); + if (lb_error > tolerance || ub_error > tolerance) { ++num_vars_out_of_bounds; + largest_error = std::max(largest_error, std::max(lb_error, ub_error)); } + var->set_solution_value(var_value); } if (num_vars_out_of_bounds > 0) { LOG(WARNING) << "Loaded a solution whose variables matched the solver's, but " << num_vars_out_of_bounds << " out of " << variables_.size() - << " had values outside their bounds."; - } - // At this point, we know for sure that the solution can be safely loaded. - // We extract the model onto the underlying solver, so that the accessors work - // properly later (as if we had called Solve()). - interface_->ExtractModel(); - // Do a second pass, where we actually set the variable values. - for (int i = 0; i < response.solution_values_size(); ++i) { - variables_[FindOrDie(variable_name_to_index_, - response.solution_values(i).variable_id())] - ->set_solution_value(response.solution_values(i).value()); + << " exceed one of their bounds by more tahn the primal tolerance: " + << tolerance; } // Set the objective value, if is known. if (response.has_objective_value()) { @@ -1579,5 +1686,81 @@ int MPSolverParameters::GetIntegerParam( } +bool MPSolver::LoadSolutionFromProto(const MPSolutionResponse& response) { + // CHANGES TO THIS METHOD ARE DISCOURAGED. + // This method will be deprecated in 2014. + // Please try using LoadSolutionFromNewProto() instead. + // TODO(user): Deprecate then remove this code. + interface_->result_status_ = + static_cast(response.result_status()); + if (response.result_status() != MPSolutionResponse::OPTIMAL && + response.result_status() != MPSolutionResponse::FEASIBLE) { + LOG(ERROR) + << "Cannot load a solution unless its status is OPTIMAL or FEASIBLE."; + return false; + } + // Before touching the variables, verify that the solution looks legit: + // each variable of the MPSolver must have its value listed exactly once, and + // each listed solution should correspond to a known variable. + if (response.solution_values_size() != variables_.size()) { + LOG(ERROR) + << "Trying to load a solution whose number of variables does not" + << " correspond to the Solver."; + return false; + } + std::vector variable_has_solution_value(variables_.size(), false); + int num_vars_out_of_bounds = 0; + for (int i = 0; i < response.solution_values_size(); ++i) { + const string& var_name = response.solution_values(i).variable_id(); + const int var_index = + FindWithDefault(variable_name_to_index_, var_name, -1); + if (var_index < 0) { + LOG(ERROR) + << "Trying to load a solution with unknown var '" << var_name << "'."; + return false; + } + if (variable_has_solution_value[var_index]) { + LOG(ERROR) + << "Trying to load a solution value where variable '" << var_name + << "' has its solution value listed twice."; + return false; + } + variable_has_solution_value[var_index] = true; + // Look further: verify the bounds. Since linear solvers yield (small) + // numerical errors, though, we just log a warning if the variables look + // like they are out of their bounds. The user should inspect the values. + const double var_value = response.solution_values(i).value(); + const MPVariable* var = variables_[var_index]; + DCHECK_EQ(var->name(), var_name); + if (var_value < var->lb() || var_value > var->ub()) { + ++num_vars_out_of_bounds; + } + } + if (num_vars_out_of_bounds > 0) { + LOG(WARNING) + << "Loaded a solution whose variables matched the solver's, but " + << num_vars_out_of_bounds << " out of " << variables_.size() + << " had values outside their bounds."; + } + // At this point, we know for sure that the solution can be safely loaded. + // We extract the model onto the underlying solver, so that the accessors work + // properly later (as if we had called Solve()). + interface_->ExtractModel(); + // Do a second pass, where we actually set the variable values. + for (int i = 0; i < response.solution_values_size(); ++i) { + variables_[FindOrDie(variable_name_to_index_, + response.solution_values(i).variable_id())] + ->set_solution_value(response.solution_values(i).value()); + } + // Set the objective value, if is known. + if (response.has_objective_value()) { + interface_->objective_value_ = response.objective_value(); + } + // Mark the status as SOLUTION_SYNCHRONIZED, so that users may inspect the + // solution normally. + interface_->sync_status_ = MPSolverInterface::SOLUTION_SYNCHRONIZED; + return true; +} + } // namespace operations_research diff --git a/src/linear_solver/linear_solver.h b/src/linear_solver/linear_solver.h index 0d55b592dd..2c6ff947c5 100644 --- a/src/linear_solver/linear_solver.h +++ b/src/linear_solver/linear_solver.h @@ -158,6 +158,13 @@ class MPSolverInterface; class MPSolverParameters; class MPVariable; +// The new MP protocol buffer format. New clients that want to work with +// protocol buffers should use this. +namespace new_proto { + class MPModelProto; + class MPModelRequest; + class MPSolutionResponse; +} // namespace new_proto // This mathematical programming (MP) solver class is the main class // though which users build and solve problems. @@ -351,12 +358,12 @@ class MPSolver { }; // Loads model from protocol buffer. - LoadStatus LoadModel(const MPModelProto& input_model); + LoadStatus LoadModelFromProto(const new_proto::MPModelProto& input_model); // Encodes the current solution in a solution response protocol buffer. // Only nonzero variable values are stored in order to reduce the // size of the MPSolutionResponse protocol buffer. - void FillSolutionResponse(MPSolutionResponse* response) const; + void FillSolutionResponseProto(new_proto::MPSolutionResponse* response) const; // Solves the model encoded by a MPModelRequest protocol buffer and // fills the solution encoded as a MPSolutionResponse. @@ -364,27 +371,28 @@ class MPSolver { // end. If you want to keep the MPSolver alive (for debugging, or for // incremental solving), you should write another version of this function // that creates the MPSolver object on the heap and returns it. - static void SolveWithProtocolBuffers(const MPModelRequest& model_request, - MPSolutionResponse* response); - + static void SolveWithProto( + const new_proto::MPModelRequest& model_request, + new_proto::MPSolutionResponse* response); // Exports model to protocol buffer. - void ExportModel(MPModelProto* output_model) const; + // TODO(user): rename to ExportModelToProto when possible. + void ExportModelToNewProto(new_proto::MPModelProto* output_model) const; - - // Load a solution encoded in a protocol buffer onto this solver. + // Load a solution encoded in a protocol buffer onto this solver for easy + // access via the MPSolver interface. // // IMPORTANT: This may only be used in conjunction with ExportModel(), // following this example: // MPSolver my_solver; // ... add variables and constraints ... - // MPModelProto model_proto; - // my_solver.ExportModel(&model_proto); - // MPSolutionResponse solver_response; + // new_proto::MPModelProto model_proto; + // my_solver.ExportModelToNewProto(&model_proto); + // new_proto::MPSolutionResponse solver_response; // // This can be replaced by a stubby call to the linear solver server. - // MPSolver::SolveWithProtocolBuffers(model_proto, &solver_response); + // MPSolver::SolveWithProto(model_proto, &solver_response); // if (solver_response.result_status() == MPSolutionResponse::OPTIMAL) { - // CHECK(my_solver.LoadSolutionFromProto(solver_response)); + // CHECK(my_solver.LoadSolutionFromNewProto(solver_response)); // ... inspect the solution using the usual API: solution_value(), etc... // } // @@ -400,12 +408,14 @@ class MPSolver { // - loading a solution with a status other than OPTIMAL / FEASIBLE. // Note: the variable and objective values aren't checked. You can use // VerifySolution() for that. - bool LoadSolutionFromProto(const MPSolutionResponse& response); + // TODO(user): rename to LoadSolutionFromProto when possible. + bool LoadSolutionFromNewProto( + const new_proto::MPSolutionResponse& response); // ----- Export model to files or strings ----- // Shortcuts to the homonymous MPModelProtoExporter methods, via - // exporting to a MPModelProto with ExportModel() (see above). + // exporting to a MPModelProto with ExportModelToNewProto() (see above). bool ExportModelAsLpFormat(bool obfuscated, string* model_str); bool ExportModelAsMpsFormat( bool fixed_format, bool obfuscated, string* model_str); @@ -548,6 +558,16 @@ class MPSolver { bool Maximization() const; bool Minimization() const; + // DEPRECATED methods. Use the explicitly listed replacement. + LoadStatus LoadModel(const MPModelProto& input_model); // LoadModelFromProto + void FillSolutionResponse( // FillSolutionResponseProto + MPSolutionResponse* response) const; + static void SolveWithProtocolBuffers( // SolveWithProto + const MPModelRequest& model_request, MPSolutionResponse* response); + void ExportModel(MPModelProto* output_model) const; // ExportModelToNewProto + bool LoadSolutionFromProto( // LoadSolutionFromNewProto + const MPSolutionResponse& response); + private: // Computes the size of the constraint with the largest number of // coefficients with index in [min_constraint_index, diff --git a/src/linear_solver/linear_solver.proto b/src/linear_solver/linear_solver.proto index 1e76f0b587..244c3fdd18 100644 --- a/src/linear_solver/linear_solver.proto +++ b/src/linear_solver/linear_solver.proto @@ -87,3 +87,14 @@ message MPSolutionResponse { repeated MPSolutionValue solution_values = 3; } +message Error { + enum Code { + NO_ERROR = 0; + INVALID_PROBLEM_TYPE = 1; + DUPLICATE_VARIABLE_ID = 2; + UNKNOWN_VARIABLE_ID = 3; + REQUEST_IS_QOD = 4; + RPC_DEADLINE_TOO_SMALL = 5; + } +} + diff --git a/src/linear_solver/linear_solver.swig b/src/linear_solver/linear_solver.swig index 5db179c089..9e82d8158e 100644 --- a/src/linear_solver/linear_solver.swig +++ b/src/linear_solver/linear_solver.swig @@ -20,6 +20,7 @@ %{ #include "linear_solver/linear_solver.h" #include "linear_solver/linear_solver.pb.h" +#include "linear_solver/linear_solver2.pb.h" #include "linear_solver/linear_solver_ext.h" %} @@ -526,6 +527,9 @@ namespace operations_research { // The following 3 methods that use protocol buffers as output // arguments are replaced by methods that return a protocol buffer, // see code below. +%ignore MPSolver::ExportModelToNewProto; +%ignore MPSolver::FillSolutionResponseProto; +%ignore MPSolver::SolveWithProto; %ignore MPSolver::ExportModel; %ignore MPSolver::FillSolutionResponse; %ignore MPSolver::SolveWithProtocolBuffers; @@ -635,6 +639,7 @@ namespace operations_research { // arguments are replaced by methods that return a protocol buffer, // see code below. %ignore MPSolver::ExportModel; +%ignore MPSolver::ExportModelToNewProto; %ignore MPSolver::FillSolutionResponse; %ignore MPSolver::SolveWithProtocolBuffers;