diff --git a/MODS.md b/MODS.md index aac270e0c6..8e617a50bc 100644 --- a/MODS.md +++ b/MODS.md @@ -1,22 +1,20 @@ # Modifications for v1.3.0 since v1.2.2 * Added (partial) Python wrapper `highspy` - * Highs::setSolution can now be used to give a solution to the simplex solver (#775) - * Highs::addVar; Highs::addVars; Highs::deleteVars(Interval/set/mask) introduced for more natural modelling - * logHeader now written as first output, even when using libraries (#784) - * Highs::presolve and Highs::postsolve completed - * Highs::resetGlobalScheduler added to reset the global scheduler - -# Planned modifications for v1.3.0 - +* write_solution_style option modified: + * "Old" HiGHS raw solution style now indicated by value `kSolutionStyleOldRaw = -1` + * Raw and pretty solution formats for Glpsol now indicated by values `kSolutionStyleGlpsolRaw = 2` and `kSolutionStyleGlpsolPretty = 3` +* Many minor fixes handling marginal LP fle behaviour * Highs::crossover completed (#815) - * scaled_model_status_ removed from Highs (#814) +* Major revisions of CMake build system + +# Planned modifications for v1.3.0 # Planned modifications beyond v1.3.0 diff --git a/src/io/FilereaderMps.cpp b/src/io/FilereaderMps.cpp index 5ae2218c68..88fe681e01 100644 --- a/src/io/FilereaderMps.cpp +++ b/src/io/FilereaderMps.cpp @@ -40,6 +40,7 @@ FilereaderRetcode FilereaderMps::readModelFromFile(const HighsOptions& options, switch (result) { case FreeFormatParserReturnCode::kSuccess: lp.ensureColwise(); + assert(model.lp_.objective_name_ != ""); return FilereaderRetcode::kOk; case FreeFormatParserReturnCode::kParserError: return FilereaderRetcode::kParserError; @@ -63,13 +64,15 @@ FilereaderRetcode FilereaderMps::readModelFromFile(const HighsOptions& options, readMps(options.log_options, filename, -1, -1, lp.num_row_, lp.num_col_, lp.sense_, lp.offset_, lp.a_matrix_.start_, lp.a_matrix_.index_, lp.a_matrix_.value_, lp.col_cost_, lp.col_lower_, lp.col_upper_, - lp.row_lower_, lp.row_upper_, lp.integrality_, lp.col_names_, - lp.row_names_, hessian.dim_, hessian.start_, hessian.index_, - hessian.value_, options.keep_n_rows); + lp.row_lower_, lp.row_upper_, lp.integrality_, lp.objective_name_, + lp.col_names_, lp.row_names_, hessian.dim_, hessian.start_, + hessian.index_, hessian.value_, lp.cost_row_location_, + options.keep_n_rows); if (return_code == FilereaderRetcode::kOk) lp.ensureColwise(); // Comment on existence of names with spaces hasNamesWithSpaces(options.log_options, lp.num_col_, lp.col_names_); hasNamesWithSpaces(options.log_options, lp.num_row_, lp.row_names_); + assert(model.lp_.objective_name_ != ""); return return_code; } diff --git a/src/io/HMPSIO.cpp b/src/io/HMPSIO.cpp index fc22f32b2e..ff4d6259a4 100644 --- a/src/io/HMPSIO.cpp +++ b/src/io/HMPSIO.cpp @@ -34,24 +34,25 @@ using std::map; // // Read file called filename. Returns 0 if OK and 1 if file can't be opened // -FilereaderRetcode readMps(const HighsLogOptions& log_options, - const std::string filename, HighsInt mxNumRow, - HighsInt mxNumCol, HighsInt& numRow, HighsInt& numCol, - ObjSense& objSense, double& objOffset, - vector& Astart, vector& Aindex, - vector& Avalue, vector& colCost, - vector& colLower, vector& colUpper, - vector& rowLower, vector& rowUpper, - vector& integerColumn, - vector& col_names, vector& row_names, - HighsInt& Qdim, vector& Qstart, - vector& Qindex, vector& Qvalue, - const HighsInt keep_n_rows) { +FilereaderRetcode readMps( + const HighsLogOptions& log_options, const std::string filename, + HighsInt mxNumRow, HighsInt mxNumCol, HighsInt& numRow, HighsInt& numCol, + ObjSense& objSense, double& objOffset, vector& Astart, + vector& Aindex, vector& Avalue, vector& colCost, + vector& colLower, vector& colUpper, + vector& rowLower, vector& rowUpper, + vector& integerColumn, std::string& objective_name, + vector& col_names, vector& row_names, + HighsInt& Qdim, vector& Qstart, vector& Qindex, + vector& Qvalue, HighsInt& cost_row_location, + const HighsInt keep_n_rows) { // MPS file buffer numRow = 0; numCol = 0; + cost_row_location = -1; objOffset = 0; objSense = ObjSense::kMinimize; + objective_name = ""; // Astart.clear() added since setting Astart.push_back(0) in // setup_clearModel() messes up the MPS read @@ -117,7 +118,12 @@ FilereaderRetcode readMps(const HighsLogOptions& log_options, if (flag[0] == 'N' && (objName == 0 || keep_n_rows == kKeepNRowsDeleteRows)) { // N-row: take the first as the objective and possibly ignore any others - if (objName == 0) objName = data[1]; + if (objName == 0) { + objName = data[1]; + std::string name(&line[4], &line[4] + 8); + objective_name = trim(name); + cost_row_location = numRow; + } } else { if (mxNumRow > 0 && numRow >= mxNumRow) return FilereaderRetcode::kParserError; @@ -540,7 +546,6 @@ HighsStatus writeModelAsMps(const HighsOptions& options, std::vector local_row_names; local_col_names.resize(lp.num_col_); local_row_names.resize(lp.num_row_); - // // Initialise the local names to any existing names if (have_col_names) local_col_names = lp.col_names_; if (have_row_names) local_row_names = lp.row_names_; @@ -576,17 +581,20 @@ HighsStatus writeModelAsMps(const HighsOptions& options, warning_found = true; } } + // Set a local objective name, creating one if necessary + const std::string local_objective_name = + findModelObjectiveName(&lp, &hessian); // If there is Hessian data to write out, writeMps assumes that hessian is // triangular if (hessian.dim_) assert(hessian.format_ == HessianFormat::kTriangular); - HighsStatus write_status = - writeMps(options.log_options, filename, lp.model_name_, lp.num_row_, - lp.num_col_, hessian.dim_, lp.sense_, lp.offset_, lp.col_cost_, - lp.col_lower_, lp.col_upper_, lp.row_lower_, lp.row_upper_, - lp.a_matrix_.start_, lp.a_matrix_.index_, lp.a_matrix_.value_, - hessian.start_, hessian.index_, hessian.value_, lp.integrality_, - local_col_names, local_row_names, use_free_format); + HighsStatus write_status = writeMps( + options.log_options, filename, lp.model_name_, lp.num_row_, lp.num_col_, + hessian.dim_, lp.sense_, lp.offset_, lp.col_cost_, lp.col_lower_, + lp.col_upper_, lp.row_lower_, lp.row_upper_, lp.a_matrix_.start_, + lp.a_matrix_.index_, lp.a_matrix_.value_, hessian.start_, hessian.index_, + hessian.value_, lp.integrality_, local_objective_name, local_col_names, + local_row_names, use_free_format); if (write_status == HighsStatus::kOk && warning_found) return HighsStatus::kWarning; return write_status; @@ -602,7 +610,7 @@ HighsStatus writeMps( const vector& a_start, const vector& a_index, const vector& a_value, const vector& q_start, const vector& q_index, const vector& q_value, - const vector& integrality, + const vector& integrality, const std::string objective_name, const vector& col_names, const vector& row_names, const bool use_free_format) { const bool write_zero_no_cost_columns = true; @@ -629,6 +637,7 @@ HighsStatus writeMps( max_name_length); return HighsStatus::kError; } + assert(objective_name != ""); vector r_ty; vector rhs, ranges; bool have_rhs = false; @@ -748,7 +757,7 @@ HighsStatus writeMps( fprintf(file, "NAME %s\n", model_name.c_str()); fprintf(file, "ROWS\n"); - fprintf(file, " N COST\n"); + fprintf(file, " N %-8s\n", objective_name.c_str()); for (HighsInt r_n = 0; r_n < num_row; r_n++) { if (r_ty[r_n] == MPS_ROW_TY_E) { fprintf(file, " E %-8s\n", row_names[r_n].c_str()); @@ -770,7 +779,8 @@ HighsStatus writeMps( if (write_zero_no_cost_columns) { // Give the column a presence by writing out a zero cost double v = 0; - fprintf(file, " %-8s COST %.15g\n", col_names[c_n].c_str(), v); + fprintf(file, " %-8s %-8s %.15g\n", col_names[c_n].c_str(), + objective_name.c_str(), v); } continue; } @@ -793,7 +803,8 @@ HighsStatus writeMps( } if (col_cost[c_n] != 0) { double v = (HighsInt)sense * col_cost[c_n]; - fprintf(file, " %-8s COST %.15g\n", col_names[c_n].c_str(), v); + fprintf(file, " %-8s %-8s %.15g\n", col_names[c_n].c_str(), + objective_name.c_str(), v); } for (HighsInt el_n = a_start[c_n]; el_n < a_start[c_n + 1]; el_n++) { double v = a_value[el_n]; @@ -814,7 +825,7 @@ HighsStatus writeMps( if (offset) { // Handle the objective offset as a RHS entry for the cost row double v = -(HighsInt)sense * offset; - fprintf(file, " RHS_V COST %.15g\n", v); + fprintf(file, " RHS_V %-8s %.15g\n", objective_name.c_str(), v); } for (HighsInt r_n = 0; r_n < num_row; r_n++) { double v = rhs[r_n]; diff --git a/src/io/HMPSIO.h b/src/io/HMPSIO.h index 701d5a8112..c0b79d6cac 100644 --- a/src/io/HMPSIO.h +++ b/src/io/HMPSIO.h @@ -53,9 +53,10 @@ FilereaderRetcode readMps( vector& Aindex, vector& Avalue, vector& colCost, vector& colLower, vector& colUpper, vector& rowLower, vector& rowUpper, - vector& integerColumn, vector& col_names, - vector& row_names, HighsInt& Qdim, vector& Qstart, - vector& Qindex, vector& Qvalue, + vector& integerColumn, std::string& objective_name, + vector& col_names, vector& row_names, + HighsInt& Qdim, vector& Qstart, vector& Qindex, + vector& Qvalue, HighsInt& cost_row_location, const HighsInt keep_n_rows = 0); HighsStatus writeMps( @@ -68,7 +69,7 @@ HighsStatus writeMps( const vector& a_start, const vector& a_index, const vector& a_value, const vector& q_start, const vector& q_index, const vector& q_value, - const vector& integrality, + const vector& integrality, std::string objective_name, const vector& col_names, const vector& row_names, const bool use_free_format = true); diff --git a/src/io/HMpsFF.cpp b/src/io/HMpsFF.cpp index 97f2f2445f..1375594d72 100644 --- a/src/io/HMpsFF.cpp +++ b/src/io/HMpsFF.cpp @@ -13,6 +13,8 @@ #include "io/HMpsFF.h" +#include "lp_data/HighsModelUtils.h" + #ifdef ZLIB_FOUND #include "zstr.hpp" #endif @@ -98,6 +100,7 @@ FreeFormatParserReturnCode HMpsFF::loadProblem( lp.row_lower_ = std::move(row_lower); lp.row_upper_ = std::move(row_upper); + lp.objective_name_ = objective_name; lp.row_names_ = std::move(row_names); lp.col_names_ = std::move(col_names); @@ -120,6 +123,10 @@ FreeFormatParserReturnCode HMpsFF::loadProblem( // 0 if (hessian.start_.size() == 0) hessian.clear(); + // Set the objective name, creating one if necessary + lp.objective_name_ = findModelObjectiveName(&lp, &hessian); + lp.cost_row_location_ = cost_row_location; + return FreeFormatParserReturnCode::kSuccess; } @@ -228,6 +235,7 @@ FreeFormatParserReturnCode HMpsFF::parse(const HighsLogOptions& log_options, num_row = 0; num_col = 0; num_nz = 0; + cost_row_location = -1; // Indicate that no duplicate rows or columns have been found has_duplicate_row_name_ = false; has_duplicate_col_name_ = false; @@ -521,7 +529,8 @@ HMpsFF::Parsekey HMpsFF::parseRows(const HighsLogOptions& log_options, std::istream& file) { std::string strline, word; bool hasobj = false; - objective_name = ""; + // Assign a default objective name + objective_name = "Objective"; assert(num_row == 0); assert(row_lower.size() == 0); @@ -565,9 +574,10 @@ HMpsFF::Parsekey HMpsFF::parseRows(const HighsLogOptions& log_options, row_upper.push_back(0.0); row_type.push_back(Boundtype::kLe); } else if (strline[start] == 'N') { - if (objective_name == "") { + if (!hasobj) { isobj = true; hasobj = true; + cost_row_location = num_row; } else { isFreeRow = true; } diff --git a/src/io/HMpsFF.h b/src/io/HMpsFF.h index 654a7ba7d7..0b855bf0b5 100644 --- a/src/io/HMpsFF.h +++ b/src/io/HMpsFF.h @@ -101,6 +101,9 @@ class HMpsFF { // any LI or UI flags in the BOUNDS section std::vector col_binary; + // Record where the cost row is encountered + HighsInt cost_row_location; + // Record whether there are duplicate row or column names, and the // name and indices of the first duplicates bool has_duplicate_row_name_; diff --git a/src/lp_data/HConst.h b/src/lp_data/HConst.h index d07937c02a..268e6edb3c 100644 --- a/src/lp_data/HConst.h +++ b/src/lp_data/HConst.h @@ -122,11 +122,20 @@ enum BasisValidity { }; enum SolutionStyle { + kSolutionStyleOldRaw = -1, kSolutionStyleRaw = 0, - kSolutionStylePretty, // 1; - kSolutionStyleOldRaw, // 2; - kSolutionStyleMin = kSolutionStyleRaw, - kSolutionStyleMax = kSolutionStyleOldRaw + kSolutionStylePretty, // 1; + kSolutionStyleGlpsolRaw, // 2; + kSolutionStyleGlpsolPretty, // 3; + kSolutionStyleMin = kSolutionStyleOldRaw, + kSolutionStyleMax = kSolutionStyleGlpsolPretty +}; + +enum GlpsolCostRowLocation { + kGlpsolCostRowLocationLast = -2, + kGlpsolCostRowLocationNone, // -1 + kGlpsolCostRowLocationNoneIfEmpty, // 0 + kGlpsolCostRowLocationMin = kGlpsolCostRowLocationLast }; const std::string kHighsFilenameDefault = ""; @@ -187,8 +196,13 @@ const HighsInt kMaxAllowedMatrixPow2Scale = 30; // Illegal values of num/max/sum infeasibility - used to indicate that true // values aren't known -const HighsInt kHighsIllegalInfeasibilityCount = -1; const double kHighsIllegalInfeasibilityMeasure = kHighsInf; +const HighsInt kHighsIllegalInfeasibilityCount = -1; + +// Illegal values for HighsError - used to indicate that true +// values aren't known +const double kHighsIllegalErrorValue = kHighsInf; +const HighsInt kHighsIllegalErrorIndex = -1; // Maximum upper bound on semi-variables const double kMaxSemiVariableUpper = 1e5; diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp index ff50ce2c88..1a05ce8701 100644 --- a/src/lp_data/Highs.cpp +++ b/src/lp_data/Highs.cpp @@ -1201,17 +1201,22 @@ HighsStatus Highs::run() { basis_.col_status = presolve_.data_.recovered_basis_.col_status; basis_.row_status = presolve_.data_.recovered_basis_.row_status; basis_.debug_origin_name += ": after postsolve"; - // Possibly force debug to perform KKT check on what's - // returned from postsolve - const bool force_debug = false; - HighsInt save_highs_debug_level = options_.highs_debug_level; - if (force_debug) - options_.highs_debug_level = kHighsDebugLevelCostly; - if (debugHighsSolution("After returning from postsolve", options_, - model_, solution_, - basis_) == HighsDebugStatus::kLogicalError) - return returnFromRun(HighsStatus::kError); - options_.highs_debug_level = save_highs_debug_level; + // Basic primal activities are wrong after postsolve, so + // possibly skip KKT check + const bool perform_kkt_check = true; + if (perform_kkt_check) { + // Possibly force debug to perform KKT check on what's + // returned from postsolve + const bool force_debug = false; + HighsInt save_highs_debug_level = options_.highs_debug_level; + if (force_debug) + options_.highs_debug_level = kHighsDebugLevelCostly; + if (debugHighsSolution("After returning from postsolve", options_, + model_, solution_, + basis_) == HighsDebugStatus::kLogicalError) + return returnFromRun(HighsStatus::kError); + options_.highs_debug_level = save_highs_debug_level; + } // Save the options to allow the best simplex strategy to // be used HighsOptions save_options = options_; @@ -2390,8 +2395,8 @@ HighsStatus Highs::writeSolution(const std::string filename, return_status = interpretCallStatus(options_.log_options, call_status, return_status, "openWriteFile"); if (return_status == HighsStatus::kError) return return_status; - writeSolutionFile(file, model_.lp_, basis_, solution_, info_, model_status_, - style); + writeSolutionFile(file, options_, model_, basis_, solution_, info_, + model_status_, style); if (style == kSolutionStyleRaw) { fprintf(file, "\n# Basis\n"); writeBasisFile(file, basis_); @@ -2560,6 +2565,8 @@ HighsPostsolveStatus Highs::runPostsolve() { presolve_.data_.postSolveStack.undo(options_, presolve_.data_.recovered_solution_, presolve_.data_.recovered_basis_); + // Compute the row activities + calculateRowValuesQuad(model_.lp_, presolve_.data_.recovered_solution_); if (have_dual_solution && model_.lp_.sense_ == ObjSense::kMaximize) presolve_.negateReducedLpColDuals(true); diff --git a/src/lp_data/HighsDeprecated.cpp b/src/lp_data/HighsDeprecated.cpp index f833a2f479..f9ff0a4d3d 100644 --- a/src/lp_data/HighsDeprecated.cpp +++ b/src/lp_data/HighsDeprecated.cpp @@ -162,7 +162,8 @@ HighsStatus Highs::writeSolution(const std::string filename, } else { style = kSolutionStyleRaw; } - writeSolutionFile(file, model_.lp_, basis_, solution_, info_, model_status_, + writeSolutionFile(file, options_, + model_, basis_, solution_, info_, model_status_, style); if (file != stdout) fclose(file); return HighsStatus::kOk; diff --git a/src/lp_data/HighsLp.cpp b/src/lp_data/HighsLp.cpp index 69ad86ac61..2a57383422 100644 --- a/src/lp_data/HighsLp.cpp +++ b/src/lp_data/HighsLp.cpp @@ -43,6 +43,7 @@ bool HighsLp::hasSemiVariables() const { bool HighsLp::operator==(const HighsLp& lp) { bool equal = equalButForNames(lp); + equal = this->objective_name_ == lp.objective_name_ && equal; equal = this->row_names_ == lp.row_names_ && equal; equal = this->col_names_ == lp.col_names_ && equal; return equal; @@ -127,6 +128,7 @@ void HighsLp::clear() { this->offset_ = 0; this->model_name_ = ""; + this->objective_name_ = ""; this->col_names_.clear(); this->row_names_.clear(); @@ -136,7 +138,7 @@ void HighsLp::clear() { this->clearScale(); this->is_scaled_ = false; this->is_moved_ = false; - + this->cost_row_location_ = -1; this->mods_.clear(); } diff --git a/src/lp_data/HighsLp.h b/src/lp_data/HighsLp.h index caf52e281e..65183f56eb 100644 --- a/src/lp_data/HighsLp.h +++ b/src/lp_data/HighsLp.h @@ -41,6 +41,7 @@ class HighsLp { double offset_; std::string model_name_; + std::string objective_name_; std::vector col_names_; std::vector row_names_; @@ -50,6 +51,7 @@ class HighsLp { HighsScale scale_; bool is_scaled_; bool is_moved_; + HighsInt cost_row_location_; HighsLpMods mods_; bool operator==(const HighsLp& lp); diff --git a/src/lp_data/HighsLpUtils.cpp b/src/lp_data/HighsLpUtils.cpp index 3d36c4e981..aa714c22d6 100644 --- a/src/lp_data/HighsLpUtils.cpp +++ b/src/lp_data/HighsLpUtils.cpp @@ -1952,39 +1952,6 @@ void analyseLp(const HighsLogOptions& log_options, const HighsLp& lp) { lp.row_upper_); } -void writeSolutionFile(FILE* file, const HighsLp& lp, const HighsBasis& basis, - const HighsSolution& solution, const HighsInfo& info, - const HighsModelStatus model_status, - const HighsInt style) { - const bool have_primal = solution.value_valid; - const bool have_dual = solution.dual_valid; - const bool have_basis = basis.valid; - const double double_tolerance = 1e-13; - if (style == kSolutionStylePretty) { - const HighsVarType* integrality_ptr = - lp.integrality_.size() > 0 ? &lp.integrality_[0] : NULL; - writeModelBoundSolution(file, true, lp.num_col_, lp.col_lower_, - lp.col_upper_, lp.col_names_, have_primal, - solution.col_value, have_dual, solution.col_dual, - have_basis, basis.col_status, integrality_ptr); - writeModelBoundSolution(file, false, lp.num_row_, lp.row_lower_, - lp.row_upper_, lp.row_names_, have_primal, - solution.row_value, have_dual, solution.row_dual, - have_basis, basis.row_status); - fprintf(file, "\nModel status: %s\n", - utilModelStatusToString(model_status).c_str()); - std::array objStr = highsDoubleToString( - (double)info.objective_function_value, double_tolerance); - fprintf(file, "\nObjective value: %s\n", objStr.data()); - } else if (style == kSolutionStyleOldRaw) { - writeOldRawSolution(file, lp, basis, solution); - } else { - fprintf(file, "Model status\n"); - fprintf(file, "%s\n", utilModelStatusToString(model_status).c_str()); - writeModelSolution(file, lp, solution, info); - } -} - HighsStatus readSolutionFile(const std::string filename, const HighsOptions& options, const HighsLp& lp, HighsBasis& basis, HighsSolution& solution, @@ -2742,68 +2709,3 @@ void removeRowsOfCountOne(const HighsLogOptions& log_options, HighsLp& lp) { highsLogUser(log_options, HighsLogType::kWarning, "Removed %d rows of count 1\n", (int)num_row_count_1); } - -void writeOldRawSolution(FILE* file, const HighsLp& lp, const HighsBasis& basis, - const HighsSolution& solution) { - const bool have_value = solution.value_valid; - const bool have_dual = solution.dual_valid; - const bool have_basis = basis.valid; - vector use_col_value; - vector use_row_value; - vector use_col_dual; - vector use_row_dual; - vector use_col_status; - vector use_row_status; - if (have_value) { - use_col_value = solution.col_value; - use_row_value = solution.row_value; - } - if (have_dual) { - use_col_dual = solution.col_dual; - use_row_dual = solution.row_dual; - } - if (have_basis) { - use_col_status = basis.col_status; - use_row_status = basis.row_status; - } - if (!have_value && !have_dual && !have_basis) return; - fprintf(file, - "%" HIGHSINT_FORMAT " %" HIGHSINT_FORMAT - " : Number of columns and rows for primal or dual solution " - "or basis\n", - lp.num_col_, lp.num_row_); - if (have_value) { - fprintf(file, "T"); - } else { - fprintf(file, "F"); - } - fprintf(file, " Primal solution\n"); - if (have_dual) { - fprintf(file, "T"); - } else { - fprintf(file, "F"); - } - fprintf(file, " Dual solution\n"); - if (have_basis) { - fprintf(file, "T"); - } else { - fprintf(file, "F"); - } - fprintf(file, " Basis\n"); - fprintf(file, "Columns\n"); - for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { - if (have_value) fprintf(file, "%.15g ", use_col_value[iCol]); - if (have_dual) fprintf(file, "%.15g ", use_col_dual[iCol]); - if (have_basis) - fprintf(file, "%" HIGHSINT_FORMAT "", (HighsInt)use_col_status[iCol]); - fprintf(file, "\n"); - } - fprintf(file, "Rows\n"); - for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { - if (have_value) fprintf(file, "%.15g ", use_row_value[iRow]); - if (have_dual) fprintf(file, "%.15g ", use_row_dual[iRow]); - if (have_basis) - fprintf(file, "%" HIGHSINT_FORMAT "", (HighsInt)use_row_status[iRow]); - fprintf(file, "\n"); - } -} diff --git a/src/lp_data/HighsLpUtils.h b/src/lp_data/HighsLpUtils.h index 0040bcc5d7..0977ad8b89 100644 --- a/src/lp_data/HighsLpUtils.h +++ b/src/lp_data/HighsLpUtils.h @@ -196,11 +196,6 @@ void getLpMatrixCoefficient(const HighsLp& lp, const HighsInt row, // Analyse the data in an LP problem void analyseLp(const HighsLogOptions& log_options, const HighsLp& lp); -void writeSolutionFile(FILE* file, const HighsLp& lp, const HighsBasis& basis, - const HighsSolution& solution, const HighsInfo& info, - const HighsModelStatus model_status, - const HighsInt style); - HighsStatus readSolutionFile(const std::string filename, const HighsOptions& options, const HighsLp& lp, HighsBasis& basis, HighsSolution& solution, @@ -238,6 +233,4 @@ HighsLp withoutSemiVariables(const HighsLp& lp); void removeRowsOfCountOne(const HighsLogOptions& log_options, HighsLp& lp); -void writeOldRawSolution(FILE* file, const HighsLp& lp, const HighsBasis& basis, - const HighsSolution& solution); #endif // LP_DATA_HIGHSLPUTILS_H_ diff --git a/src/lp_data/HighsModelUtils.cpp b/src/lp_data/HighsModelUtils.cpp index a7afb580a3..8cdc8583ca 100644 --- a/src/lp_data/HighsModelUtils.cpp +++ b/src/lp_data/HighsModelUtils.cpp @@ -17,13 +17,19 @@ #include "lp_data/HighsModelUtils.h" #include +#include #include #include -#include "HConfig.h" -#include "io/HighsIO.h" -#include "lp_data/HConst.h" -#include "util/HighsUtils.h" +//#include "model/HighsModel.h" +//#include "HConfig.h" +//#include "io/HighsIO.h" +//#include "lp_data/HConst.h" +#include "lp_data/HighsSolution.h" +#include "util/stringutil.h" + +const double kHighsDoubleTolerance = 1e-13; +const double kGlpsolDoubleTolerance = 1e-12; void analyseModelBounds(const HighsLogOptions& log_options, const char* message, HighsInt numBd, const std::vector& lower, @@ -210,7 +216,6 @@ void writeModelSolution(FILE* file, const HighsLp& lp, assert((int)solution.row_dual.size() >= lp.num_row_); assert(info.dual_solution_status != kSolutionStatusNone); } - const double double_tolerance = 1e-13; fprintf(file, "\n# Primal solution values\n"); if (!have_primal || info.primal_solution_status == kSolutionStatusNone) { fprintf(file, "None\n"); @@ -224,13 +229,13 @@ void writeModelSolution(FILE* file, const HighsLp& lp, HighsCDouble objective_function_value = lp.offset_; for (HighsInt i = 0; i < lp.num_col_; ++i) objective_function_value += lp.col_cost_[i] * solution.col_value[i]; - std::array objStr = - highsDoubleToString((double)objective_function_value, double_tolerance); + std::array objStr = highsDoubleToString( + (double)objective_function_value, kHighsDoubleTolerance); fprintf(file, "Objective %s\n", objStr.data()); fprintf(file, "# Columns %" HIGHSINT_FORMAT "\n", lp.num_col_); for (HighsInt ix = 0; ix < lp.num_col_; ix++) { std::array valStr = - highsDoubleToString(solution.col_value[ix], double_tolerance); + highsDoubleToString(solution.col_value[ix], kHighsDoubleTolerance); // Create a column name ss.str(std::string()); ss << "C" << ix; @@ -240,7 +245,7 @@ void writeModelSolution(FILE* file, const HighsLp& lp, fprintf(file, "# Rows %" HIGHSINT_FORMAT "\n", lp.num_row_); for (HighsInt ix = 0; ix < lp.num_row_; ix++) { std::array valStr = - highsDoubleToString(solution.row_value[ix], double_tolerance); + highsDoubleToString(solution.row_value[ix], kHighsDoubleTolerance); // Create a row name ss.str(std::string()); ss << "R" << ix; @@ -261,7 +266,7 @@ void writeModelSolution(FILE* file, const HighsLp& lp, fprintf(file, "# Columns %" HIGHSINT_FORMAT "\n", lp.num_col_); for (HighsInt ix = 0; ix < lp.num_col_; ix++) { std::array valStr = - highsDoubleToString(solution.col_dual[ix], double_tolerance); + highsDoubleToString(solution.col_dual[ix], kHighsDoubleTolerance); ss.str(std::string()); ss << "C" << ix; const std::string name = have_col_names ? lp.col_names_[ix] : ss.str(); @@ -270,7 +275,7 @@ void writeModelSolution(FILE* file, const HighsLp& lp, fprintf(file, "# Rows %" HIGHSINT_FORMAT "\n", lp.num_row_); for (HighsInt ix = 0; ix < lp.num_row_; ix++) { std::array valStr = - highsDoubleToString(solution.row_dual[ix], double_tolerance); + highsDoubleToString(solution.row_dual[ix], kHighsDoubleTolerance); ss.str(std::string()); ss << "R" << ix; const std::string name = have_row_names ? lp.row_names_[ix] : ss.str(); @@ -351,6 +356,717 @@ HighsStatus normaliseNames(const HighsLogOptions& log_options, return HighsStatus::kOk; } +void writeSolutionFile(FILE* file, const HighsOptions& options, + const HighsModel& model, const HighsBasis& basis, + const HighsSolution& solution, const HighsInfo& info, + const HighsModelStatus model_status, + const HighsInt style) { + const bool have_primal = solution.value_valid; + const bool have_dual = solution.dual_valid; + const bool have_basis = basis.valid; + const HighsLp& lp = model.lp_; + if (style == kSolutionStyleOldRaw) { + writeOldRawSolution(file, lp, basis, solution); + } else if (style == kSolutionStylePretty) { + const HighsVarType* integrality_ptr = + lp.integrality_.size() > 0 ? &lp.integrality_[0] : NULL; + writeModelBoundSolution(file, true, lp.num_col_, lp.col_lower_, + lp.col_upper_, lp.col_names_, have_primal, + solution.col_value, have_dual, solution.col_dual, + have_basis, basis.col_status, integrality_ptr); + writeModelBoundSolution(file, false, lp.num_row_, lp.row_lower_, + lp.row_upper_, lp.row_names_, have_primal, + solution.row_value, have_dual, solution.row_dual, + have_basis, basis.row_status); + fprintf(file, "\nModel status: %s\n", + utilModelStatusToString(model_status).c_str()); + std::array objStr = highsDoubleToString( + (double)info.objective_function_value, kHighsDoubleTolerance); + fprintf(file, "\nObjective value: %s\n", objStr.data()); + } else if (style == kSolutionStyleGlpsolRaw || + style == kSolutionStyleGlpsolPretty) { + const bool raw = style == kSolutionStyleGlpsolRaw; + writeGlpsolSolution(file, options, model, basis, solution, model_status, + info, raw); + } else { + fprintf(file, "Model status\n"); + fprintf(file, "%s\n", utilModelStatusToString(model_status).c_str()); + writeModelSolution(file, lp, solution, info); + } +} + +void writeGlpsolCostRow(FILE* file, const bool raw, const bool is_mip, + const HighsInt row_id, const std::string objective_name, + const double objective_function_value) { + if (raw) { + double double_value = objective_function_value; + std::array double_string = + highsDoubleToString(double_value, kGlpsolDoubleTolerance); + // Last term of 0 for dual should (also) be blank when not MIP + fprintf(file, "i %d %s%s%s\n", (int)row_id, is_mip ? "" : "b ", + double_string.data(), is_mip ? "" : " 0"); + } else { + fprintf(file, "%6d ", (int)row_id); + if (objective_name.length() <= 12) { + fprintf(file, "%-12s ", objective_name.c_str()); + } else { + fprintf(file, "%s\n%20s", objective_name.c_str(), ""); + } + if (is_mip) { + fprintf(file, " "); + } else { + fprintf(file, "B "); + } + fprintf(file, "%13.6g %13s %13s \n", objective_function_value, "", ""); + } +} + +void writeGlpsolSolution(FILE* file, const HighsOptions& options, + const HighsModel& model, const HighsBasis& basis, + const HighsSolution& solution, + const HighsModelStatus model_status, + const HighsInfo& info, const bool raw) { + const bool have_value = solution.value_valid; + const bool have_dual = solution.dual_valid; + const bool have_basis = basis.valid; + if (!have_value) return; + const double kGlpsolHighQuality = 1e-9; + const double kGlpsolMediumQuality = 1e-6; + const double kGlpsolLowQuality = 1e-3; + const double kGlpsolPrintAsZero = 1e-9; + const HighsLp& lp = model.lp_; + const bool have_col_names = lp.col_names_.size(); + const bool have_row_names = lp.row_names_.size(); + // Determine number of nonzeros including the objective function + // and, hence, determine whether there is an objective function + HighsInt num_nz = lp.a_matrix_.numNz(); + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) + if (lp.col_cost_[iCol]) num_nz++; + bool empty_cost_row = num_nz == lp.a_matrix_.numNz(); + bool has_objective = !empty_cost_row; + if (model.hessian_.dim_) has_objective = true; + // When writing out the row information (and hence the number of + // rows and nonzeros), the case of the cost row is tricky + // (particularly if it's empty) if HiGHS is to be able to reproduce + // the (inconsistent) behaviour of Glpsol. + // + // If Glpsol is run from a .mod file then the cost row is reported + // unless it is empty. However, its position depends on where the + // objecive appears in the .mod file. If Glpsol is run from a .mod + // file, and reads a .sol file, it must be in the right format. + // + // HiGHS can't read ..mod files, so works from an MPS or LP file + // generated by glpsol. + // + // An MPS file generated by glpsol will have the cost row in the + // same position as it was in the .mod file + // + // An LP file generated by glpsol will have the objective defined + // first, so the desired position of the cost row in the .sol file + // is unavailable. The only option with this route is to define the + // cost row location "by hand" using glpsol_cost_row_location + // + // If Glpsol is run from an LP or MPS file then the cost row is not + // reported. This behaviour is defined by setting + // glpsol_cost_row_location = -1; + // + // This inconsistent behaviour means that it must be possible to + // tell HiGHS to suppress the cost row + // + HighsInt cost_row_option = options.glpsol_cost_row_location; + // Define cost_row_location + // + // It is indexed from 1 so that it matches the index printed on that + // row... + // + // ... hence a location of zero means that the cost row isn't + // reported + HighsInt cost_row_location = 0; + + if (cost_row_option <= kGlpsolCostRowLocationLast || + cost_row_option > lp.num_row_) { + // Place the cost row last + cost_row_location = lp.num_row_ + 1; + } else if (cost_row_option == kGlpsolCostRowLocationNone) { + // Don't report the cost row + assert(cost_row_location == 0); + } else if (cost_row_option == kGlpsolCostRowLocationNoneIfEmpty) { + // This option allows the cost row to be omitted if it's empty. + if (empty_cost_row) { + // The cost row is empty, so don't report it + assert(cost_row_location == 0); + } else { + // Place the cost row according to lp.cost_row_location_ + if (lp.cost_row_location_ >= 0) { + // The cost row location is known from the MPS file. NB To + // index from zero whenever possible, lp.cost_row_location_ = + // 0 if the cost row came first + assert(lp.cost_row_location_ <= lp.num_row_); + cost_row_location = lp.cost_row_location_ + 1; + } else { + // The location isn't known from an MPS file, so place it + // last, giving a warning + cost_row_location = lp.num_row_ + 1; + highsLogUser( + options.log_options, HighsLogType::kWarning, + "The cost row for the Glpsol solution file is reported last since " + "there is no indication of where it should be\n"); + } + } + } else { + // Place the cost row according to the option value + cost_row_location = cost_row_option; + } + assert(0 <= cost_row_location && cost_row_location <= lp.num_row_ + 1); + // Despite being written in C, GLPSOL indexes rows (columns) from + // 1..m (1..n) with - bizarrely! - m being one more than the number + // of constraints if the cost vector is reported. + const HighsInt num_row = lp.num_row_; + const HighsInt num_col = lp.num_col_; + // There's one more row and more nonzeros if the cost row is + // reported + const HighsInt delta_num_row = cost_row_location > 0; + const HighsInt glpsol_num_row = num_row + delta_num_row; + // If the cost row isn't reported, then the number of nonzeros is + // just the number in the constraint matrix + if (cost_row_location <= 0) num_nz = lp.a_matrix_.numNz(); + // Record the discrete nature of the model + HighsInt num_integer = 0; + HighsInt num_binary = 0; + bool is_mip = false; + if (lp.integrality_.size() == lp.num_col_) { + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { + if (lp.integrality_[iCol] != HighsVarType::kContinuous) { + is_mip = true; + num_integer++; + if (lp.col_lower_[iCol] == 0 && lp.col_upper_[iCol] == 1) num_binary++; + } + } + } + // Raw and pretty are the initially the same, but for the "c " + // prefix to raw lines + std::string line_prefix = ""; + if (raw) line_prefix = "c "; + fprintf(file, "%s%-12s%s\n", line_prefix.c_str(), + "Problem:", lp.model_name_.c_str()); + fprintf(file, "%s%-12s%d\n", line_prefix.c_str(), + "Rows:", (int)glpsol_num_row); + fprintf(file, "%s%-12s%d", line_prefix.c_str(), "Columns:", (int)num_col); + if (!raw && is_mip) + fprintf(file, " (%d integer, %d binary)", (int)num_integer, + (int)num_binary); + fprintf(file, "\n"); + fprintf(file, "%s%-12s%d\n", line_prefix.c_str(), "Non-zeros:", (int)num_nz); + std::string model_status_text = ""; + std::string model_status_char = ""; // Just for raw MIP solution + switch (model_status) { + case HighsModelStatus::kOptimal: + if (is_mip) model_status_text = "INTEGER "; + model_status_text += "OPTIMAL"; + model_status_char = "o"; + break; + case HighsModelStatus::kInfeasible: + if (is_mip) { + model_status_text = "INTEGER INFEASIBLE"; + } else { + model_status_text = "INFEASIBLE (FINAL)"; + } + model_status_char = "n"; + break; + case HighsModelStatus::kUnboundedOrInfeasible: + model_status_text += "INFEASIBLE (INTERMEDIATE)"; + model_status_char = "?"; // Glpk has no equivalent + break; + case HighsModelStatus::kUnbounded: + if (is_mip) model_status_text = "INTEGER "; + model_status_text += "UNBOUNDED"; + model_status_char = "u"; // Glpk has no equivalent + break; + default: + model_status_text = "????"; + model_status_char = "?"; + break; + } + assert(model_status_text != ""); + assert(model_status_char != ""); + fprintf(file, "%s%-12s%s\n", line_prefix.c_str(), + "Status:", model_status_text.c_str()); + // If info is not valid, then cannot write more + if (!info.valid) return; + // Now write out the numerical information + // + // Determine the objective name to write out + std::string objective_name = lp.objective_name_; + // Make sure that no objective name is written out if there are rows + // and no row names + if (lp.num_row_ && !have_row_names) objective_name = ""; + // if there are row names to be written out, there must be a + // non-trivial objective name + if (have_row_names) assert(lp.objective_name_ != ""); + const bool has_objective_name = lp.objective_name_ != ""; + fprintf(file, "%s%-12s%s%.10g (%s)\n", line_prefix.c_str(), "Objective:", + !(has_objective && has_objective_name) + ? "" + : (objective_name + " = ").c_str(), + has_objective ? info.objective_function_value : 0, + lp.sense_ == ObjSense::kMinimize ? "MINimum" : "MAXimum"); + // No space after "c" on blank line! + if (raw) line_prefix = "c"; + fprintf(file, "%s\n", line_prefix.c_str()); + // Detailed lines are rather different + if (raw) { + fprintf(file, "s %s %d %d ", is_mip ? "mip" : "bas", (int)glpsol_num_row, + (int)num_col); + if (is_mip) { + fprintf(file, "%s", model_status_char.c_str()); + } else { + if (info.primal_solution_status == kSolutionStatusNone) { + fprintf(file, "u"); + } else if (info.primal_solution_status == kSolutionStatusInfeasible) { + fprintf(file, "i"); + } else if (info.primal_solution_status == kSolutionStatusFeasible) { + fprintf(file, "f"); + } else { + fprintf(file, "?"); + } + fprintf(file, " "); + if (info.dual_solution_status == kSolutionStatusNone) { + fprintf(file, "u"); + } else if (info.dual_solution_status == kSolutionStatusInfeasible) { + fprintf(file, "i"); + } else if (info.dual_solution_status == kSolutionStatusFeasible) { + fprintf(file, "f"); + } else { + fprintf(file, "?"); + } + } + double double_value = has_objective ? info.objective_function_value : 0; + std::array double_string = + highsDoubleToString(double_value, kHighsDoubleTolerance); + fprintf(file, " %s\n", double_string.data()); + } + if (!raw) { + fprintf(file, + " No. Row name %s Activity Lower bound " + " Upper bound", + have_basis ? "St" : " "); + if (have_dual) fprintf(file, " Marginal"); + fprintf(file, "\n"); + + fprintf(file, + "------ ------------ %s ------------- ------------- " + "-------------", + have_basis ? "--" : " "); + if (have_dual) fprintf(file, " -------------"); + fprintf(file, "\n"); + } + + HighsInt row_id = 0; + for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { + row_id++; + if (row_id == cost_row_location) { + writeGlpsolCostRow(file, raw, is_mip, row_id, objective_name, + info.objective_function_value); + row_id++; + } + if (raw) { + fprintf(file, "i %d ", (int)row_id); + if (is_mip) { + // Complete the line if for a MIP + double double_value = solution.row_value[iRow]; + std::array double_string = + highsDoubleToString(double_value, kHighsDoubleTolerance); + fprintf(file, "%s\n", double_string.data()); + continue; + } + } else { + fprintf(file, "%6d ", (int)row_id); + std::string row_name = ""; + if (have_row_names) row_name = lp.row_names_[iRow]; + if (row_name.length() <= 12) { + fprintf(file, "%-12s ", row_name.c_str()); + } else { + fprintf(file, "%s\n%20s", row_name.c_str(), ""); + } + } + const double lower = lp.row_lower_[iRow]; + const double upper = lp.row_upper_[iRow]; + const double value = solution.row_value[iRow]; + const double dual = have_dual ? solution.row_dual[iRow] : 0; + std::string status_text = " "; + std::string status_char = ""; + if (have_basis) { + const HighsBasisStatus status = basis.row_status[iRow]; + switch (basis.row_status[iRow]) { + case HighsBasisStatus::kBasic: + status_text = "B "; + status_char = "b"; + break; + case HighsBasisStatus::kLower: + status_text = lower == upper ? "NS" : "NL"; + status_char = lower == upper ? "s" : "l"; + break; + case HighsBasisStatus::kUpper: + status_text = lower == upper ? "NS" : "NU"; + status_char = lower == upper ? "s" : "u"; + break; + case HighsBasisStatus::kZero: + status_text = "NF"; + status_char = "f"; + break; + default: + status_text = "??"; + status_char = "?"; + break; + } + } + if (raw) { + fprintf(file, "%s ", status_char.c_str()); + double double_value = solution.row_value[iRow]; + std::array double_string = + highsDoubleToString(double_value, kHighsDoubleTolerance); + fprintf(file, "%s ", double_string.data()); + } else { + fprintf(file, "%s ", status_text.c_str()); + fprintf(file, "%13.6g ", fabs(value) <= kGlpsolPrintAsZero ? 0.0 : value); + if (lower > -kHighsInf) + fprintf(file, "%13.6g ", lower); + else + fprintf(file, "%13s ", ""); + if (lower != upper && upper < kHighsInf) + fprintf(file, "%13.6g ", upper); + else + fprintf(file, "%13s ", lower == upper ? "=" : ""); + } + if (have_dual) { + if (raw) { + double double_value = solution.row_dual[iRow]; + std::array double_string = + highsDoubleToString(double_value, kHighsDoubleTolerance); + fprintf(file, "%s", double_string.data()); + } else { + // If the row is known to be basic, don't print the dual + // value. If there's no basis, row cannot be known to be basic + bool not_basic = have_basis; + if (have_basis) + not_basic = basis.row_status[iRow] != HighsBasisStatus::kBasic; + if (not_basic) { + if (fabs(dual) <= kGlpsolPrintAsZero) + fprintf(file, "%13s", "< eps"); + else + fprintf(file, "%13.6g ", dual); + } + } + } + fprintf(file, "\n"); + } + + if (cost_row_location == lp.num_row_ + 1) { + row_id++; + writeGlpsolCostRow(file, raw, is_mip, row_id, objective_name, + info.objective_function_value); + } + if (!raw) fprintf(file, "\n"); + + if (!raw) { + fprintf(file, + " No. Column name %s Activity Lower bound " + " Upper bound", + have_basis ? "St" : " "); + if (have_dual) fprintf(file, " Marginal"); + fprintf(file, "\n"); + fprintf(file, + "------ ------------ %s ------------- ------------- " + "-------------", + have_basis ? "--" : " "); + if (have_dual) fprintf(file, " -------------"); + fprintf(file, "\n"); + } + + if (raw) line_prefix = "j "; + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { + if (raw) { + fprintf(file, "%s%d ", line_prefix.c_str(), (int)(iCol + 1)); + if (is_mip) { + double double_value = solution.col_value[iCol]; + std::array double_string = + highsDoubleToString(double_value, kHighsDoubleTolerance); + fprintf(file, "%s\n", double_string.data()); + continue; + } + } else { + fprintf(file, "%6d ", (int)(iCol + 1)); + std::string col_name = ""; + if (have_col_names) col_name = lp.col_names_[iCol]; + if (!have_col_names || col_name.length() <= 12) { + fprintf(file, "%-12s ", !have_col_names ? "" : col_name.c_str()); + } else { + fprintf(file, "%s\n%20s", col_name.c_str(), ""); + } + } + const double lower = lp.col_lower_[iCol]; + const double upper = lp.col_upper_[iCol]; + const double value = solution.col_value[iCol]; + const double dual = have_dual ? solution.col_dual[iCol] : 0; + std::string status_text = " "; + std::string status_char = ""; + if (have_basis) { + const HighsBasisStatus status = basis.col_status[iCol]; + switch (basis.col_status[iCol]) { + case HighsBasisStatus::kBasic: + status_text = "B "; + status_char = "b"; + break; + case HighsBasisStatus::kLower: + status_text = lower == upper ? "NS" : "NL"; + status_char = lower == upper ? "s" : "l"; + break; + case HighsBasisStatus::kUpper: + status_text = lower == upper ? "NS" : "NU"; + status_char = lower == upper ? "s" : "u"; + break; + case HighsBasisStatus::kZero: + status_text = "NF"; + status_char = "f"; + break; + default: + status_text = "??"; + status_char = "?"; + break; + } + } else if (is_mip) { + if (lp.integrality_[iCol] != HighsVarType::kContinuous) + status_text = "* "; + } + if (raw) { + fprintf(file, "%s ", status_char.c_str()); + double double_value = solution.col_value[iCol]; + std::array double_string = + highsDoubleToString(double_value, kHighsDoubleTolerance); + fprintf(file, "%s ", double_string.data()); + } else { + fprintf(file, "%s ", status_text.c_str()); + fprintf(file, "%13.6g ", fabs(value) <= kGlpsolPrintAsZero ? 0.0 : value); + if (lower > -kHighsInf) + fprintf(file, "%13.6g ", lower); + else + fprintf(file, "%13s ", ""); + if (lower != upper && upper < kHighsInf) + fprintf(file, "%13.6g ", upper); + else + fprintf(file, "%13s ", lower == upper ? "=" : ""); + } + if (have_dual) { + if (raw) { + double double_value = solution.col_dual[iCol]; + std::array double_string = + highsDoubleToString(double_value, kHighsDoubleTolerance); + fprintf(file, "%s", double_string.data()); + } else { + // If the column is known to be basic, don't print the dual + // value. If there's no basis, column cannot be known to be + // basic + bool not_basic = have_basis; + if (have_basis) + not_basic = basis.col_status[iCol] != HighsBasisStatus::kBasic; + if (not_basic) { + if (fabs(dual) <= kGlpsolPrintAsZero) + fprintf(file, "%13s", "< eps"); + else + fprintf(file, "%13.6g ", dual); + } + } + } + fprintf(file, "\n"); + } + if (raw) { + fprintf(file, "e o f\n"); + return; + } + HighsPrimalDualErrors errors; + HighsInfo local_info; + HighsInt absolute_error_index; + double absolute_error_value; + HighsInt relative_error_index; + double relative_error_value; + getKktFailures(options, model, solution, basis, local_info, errors, true); + fprintf(file, "\n"); + if (is_mip) { + fprintf(file, "Integer feasibility conditions:\n"); + } else { + fprintf(file, "Karush-Kuhn-Tucker optimality conditions:\n"); + } + fprintf(file, "\n"); + // Primal residual + absolute_error_value = errors.max_primal_residual.absolute_value; + absolute_error_index = errors.max_primal_residual.absolute_index + 1; + relative_error_value = errors.max_primal_residual.relative_value; + relative_error_index = errors.max_primal_residual.relative_index + 1; + if (!absolute_error_value) absolute_error_index = 0; + if (!relative_error_value) relative_error_index = 0; + fprintf(file, "KKT.PE: max.abs.err = %.2e on row %d\n", absolute_error_value, + absolute_error_index == 0 ? 0 : (int)absolute_error_index); + fprintf(file, " max.rel.err = %.2e on row %d\n", relative_error_value, + absolute_error_index == 0 ? 0 : (int)relative_error_index); + fprintf(file, "%8s%s\n", "", + relative_error_value <= kGlpsolHighQuality + ? "High quality" + : relative_error_value <= kGlpsolMediumQuality + ? "Medium quality" + : relative_error_value <= kGlpsolLowQuality + ? "Low quality" + : "PRIMAL SOLUTION IS WRONG"); + fprintf(file, "\n"); + + // Primal infeasibility + absolute_error_value = errors.max_primal_infeasibility.absolute_value; + absolute_error_index = errors.max_primal_infeasibility.absolute_index + 1; + relative_error_value = errors.max_primal_infeasibility.relative_value; + relative_error_index = errors.max_primal_infeasibility.relative_index + 1; + if (!absolute_error_value) absolute_error_index = 0; + if (!relative_error_value) relative_error_index = 0; + bool on_col = absolute_error_index > 0 && absolute_error_index <= lp.num_col_; + fprintf(file, "KKT.PB: max.abs.err = %.2e on %s %d\n", absolute_error_value, + on_col ? "column" : "row", + absolute_error_index <= lp.num_col_ + ? (int)absolute_error_index + : (int)(absolute_error_index - lp.num_col_)); + on_col = relative_error_index > 0 && relative_error_index <= lp.num_col_; + fprintf(file, " max.rel.err = %.2e on %s %d\n", relative_error_value, + on_col ? "column" : "row", + relative_error_index <= lp.num_col_ + ? (int)relative_error_index + : (int)(relative_error_index - lp.num_col_)); + fprintf(file, "%8s%s\n", "", + relative_error_value <= kGlpsolHighQuality + ? "High quality" + : relative_error_value <= kGlpsolMediumQuality + ? "Medium quality" + : relative_error_value <= kGlpsolLowQuality + ? "Low quality" + : "PRIMAL SOLUTION IS INFEASIBLE"); + fprintf(file, "\n"); + + if (have_dual) { + // Dual residual + absolute_error_value = errors.max_dual_residual.absolute_value; + absolute_error_index = errors.max_dual_residual.absolute_index + 1; + relative_error_value = errors.max_dual_residual.relative_value; + relative_error_index = errors.max_dual_residual.relative_index + 1; + if (!absolute_error_value) absolute_error_index = 0; + if (!relative_error_value) relative_error_index = 0; + fprintf(file, "KKT.DE: max.abs.err = %.2e on column %d\n", + absolute_error_value, (int)absolute_error_index); + fprintf(file, " max.rel.err = %.2e on column %d\n", + relative_error_value, (int)relative_error_index); + fprintf(file, "%8s%s\n", "", + relative_error_value <= kGlpsolHighQuality + ? "High quality" + : relative_error_value <= kGlpsolMediumQuality + ? "Medium quality" + : relative_error_value <= kGlpsolLowQuality + ? "Low quality" + : "DUAL SOLUTION IS WRONG"); + fprintf(file, "\n"); + + // Dual infeasibility + absolute_error_value = errors.max_dual_infeasibility.absolute_value; + absolute_error_index = errors.max_dual_infeasibility.absolute_index + 1; + relative_error_value = errors.max_dual_infeasibility.relative_value; + relative_error_index = errors.max_dual_infeasibility.relative_index + 1; + if (!absolute_error_value) absolute_error_index = 0; + if (!relative_error_value) relative_error_index = 0; + bool on_col = + absolute_error_index > 0 && absolute_error_index <= lp.num_col_; + fprintf(file, "KKT.DB: max.abs.err = %.2e on %s %d\n", absolute_error_value, + on_col ? "column" : "row", + absolute_error_index <= lp.num_col_ + ? (int)absolute_error_index + : (int)(absolute_error_index - lp.num_col_)); + on_col = relative_error_index > 0 && relative_error_index <= lp.num_col_; + fprintf(file, " max.rel.err = %.2e on %s %d\n", relative_error_value, + on_col ? "column" : "row", + relative_error_index <= lp.num_col_ + ? (int)relative_error_index + : (int)(relative_error_index - lp.num_col_)); + fprintf(file, "%8s%s\n", "", + relative_error_value <= kGlpsolHighQuality + ? "High quality" + : relative_error_value <= kGlpsolMediumQuality + ? "Medium quality" + : relative_error_value <= kGlpsolLowQuality + ? "Low quality" + : "DUAL SOLUTION IS INFEASIBLE"); + fprintf(file, "\n"); + } + fprintf(file, "End of output\n"); +} + +void writeOldRawSolution(FILE* file, const HighsLp& lp, const HighsBasis& basis, + const HighsSolution& solution) { + const bool have_value = solution.value_valid; + const bool have_dual = solution.dual_valid; + const bool have_basis = basis.valid; + vector use_col_value; + vector use_row_value; + vector use_col_dual; + vector use_row_dual; + vector use_col_status; + vector use_row_status; + if (have_value) { + use_col_value = solution.col_value; + use_row_value = solution.row_value; + } + if (have_dual) { + use_col_dual = solution.col_dual; + use_row_dual = solution.row_dual; + } + if (have_basis) { + use_col_status = basis.col_status; + use_row_status = basis.row_status; + } + if (!have_value && !have_dual && !have_basis) return; + fprintf(file, + "%" HIGHSINT_FORMAT " %" HIGHSINT_FORMAT + " : Number of columns and rows for primal or dual solution " + "or basis\n", + lp.num_col_, lp.num_row_); + if (have_value) { + fprintf(file, "T"); + } else { + fprintf(file, "F"); + } + fprintf(file, " Primal solution\n"); + if (have_dual) { + fprintf(file, "T"); + } else { + fprintf(file, "F"); + } + fprintf(file, " Dual solution\n"); + if (have_basis) { + fprintf(file, "T"); + } else { + fprintf(file, "F"); + } + fprintf(file, " Basis\n"); + fprintf(file, "Columns\n"); + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { + if (have_value) fprintf(file, "%.15g ", use_col_value[iCol]); + if (have_dual) fprintf(file, "%.15g ", use_col_dual[iCol]); + if (have_basis) + fprintf(file, "%" HIGHSINT_FORMAT "", (HighsInt)use_col_status[iCol]); + fprintf(file, "\n"); + } + fprintf(file, "Rows\n"); + for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { + if (have_value) fprintf(file, "%.15g ", use_row_value[iRow]); + if (have_dual) fprintf(file, "%.15g ", use_row_dual[iRow]); + if (have_basis) + fprintf(file, "%" HIGHSINT_FORMAT "", (HighsInt)use_row_status[iRow]); + fprintf(file, "\n"); + } +} + HighsBasisStatus checkedVarHighsNonbasicStatus( const HighsBasisStatus ideal_status, const double lower, const double upper) { @@ -535,3 +1251,52 @@ HighsStatus highsStatusFromHighsModelStatus(HighsModelStatus model_status) { return HighsStatus::kError; } } + +std::string findModelObjectiveName(const HighsLp* lp, + const HighsHessian* hessian) { + // Return any non-trivial current objective name + if (lp->objective_name_ != "") return lp->objective_name_; + + std::string objective_name = ""; + // Determine whether there is a nonzero cost vector + bool has_objective = false; + for (HighsInt iCol = 0; iCol < lp->num_col_; iCol++) { + if (lp->col_cost_[iCol]) { + has_objective = true; + break; + } + } + if (!has_objective && hessian) { + // Zero cost vector, so only chance of an objective comes from any + // Hessian + has_objective = hessian->dim_; + } + HighsInt pass = 0; + for (;;) { + // Loop until a valid name is found. Vanishingly unlikely to have + // to pass more than once, since check for objective name + // duplicating a row name is very unlikely to fail + // + // So set up an appropriate name (stem) + if (has_objective) { + objective_name = "Obj"; + } else { + objective_name = "NoObj"; + } + if (pass) objective_name += pass; + // Ensure that the objective name doesn't clash with any row names + bool ok_objective_name = true; + for (HighsInt iRow = 0; iRow < lp->num_row_; iRow++) { + std::string trimmed_name = lp->row_names_[iRow]; + trimmed_name = trim(trimmed_name); + if (objective_name == trimmed_name) { + ok_objective_name = false; + break; + } + } + if (ok_objective_name) break; + pass++; + } + assert(objective_name != ""); + return objective_name; +} diff --git a/src/lp_data/HighsModelUtils.h b/src/lp_data/HighsModelUtils.h index 54dfa90053..7d5bff59c9 100644 --- a/src/lp_data/HighsModelUtils.h +++ b/src/lp_data/HighsModelUtils.h @@ -18,10 +18,12 @@ //#include "Highs.h" //#include "lp_data/HighsStatus.h" -#include "lp_data/HStruct.h" #include "lp_data/HighsInfo.h" -#include "lp_data/HighsLp.h" -#include "lp_data/HighsOptions.h" +#include "model/HighsModel.h" +//#include "lp_data/HStruct.h" +//#include "lp_data/HighsInfo.h" +//#include "lp_data/HighsLp.h" +//#include "lp_data/HighsOptions.h" // Analyse lower and upper bounds of a model void analyseModelBounds(const HighsLogOptions& log_options, const char* message, @@ -48,6 +50,25 @@ HighsStatus normaliseNames(const HighsLogOptions& log_options, std::vector& names, HighsInt& max_name_length); +void writeSolutionFile(FILE* file, const HighsOptions& options, + const HighsModel& model, const HighsBasis& basis, + const HighsSolution& solution, const HighsInfo& info, + const HighsModelStatus model_status, + const HighsInt style); + +void writeGlpsolCostRow(FILE* file, const bool raw, const bool is_mip, + const HighsInt row_id, const std::string objective_name, + const double objective_function_value); + +void writeGlpsolSolution(FILE* file, const HighsOptions& options, + const HighsModel& model, const HighsBasis& basis, + const HighsSolution& solution, + const HighsModelStatus model_status, + const HighsInfo& info, const bool raw); + +void writeOldRawSolution(FILE* file, const HighsLp& lp, const HighsBasis& basis, + const HighsSolution& solution); + HighsBasisStatus checkedVarHighsNonbasicStatus( const HighsBasisStatus ideal_status, const double lower, const double upper); @@ -65,4 +86,8 @@ HighsStatus highsStatusFromHighsModelStatus(HighsModelStatus model_status); std::string statusToString(const HighsBasisStatus status, const double lower, const double upper); std::string typeToString(const HighsVarType type); + +std::string findModelObjectiveName(const HighsLp* lp, + const HighsHessian* hessian = nullptr); + #endif diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index b769504e85..7d8b588709 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -302,6 +302,7 @@ struct HighsOptionsStruct { bool write_model_to_file; bool write_solution_to_file; HighsInt write_solution_style; + HighsInt glpsol_cost_row_location; // Control of HiGHS log bool output_flag; bool log_to_console; @@ -636,12 +637,24 @@ class HighsOptions : public HighsOptionsStruct { record_int = new OptionRecordInt("write_solution_style", - "Write the solution in style: 0=>Raw " - "(computer-readable); 1=>Pretty (human-readable) ", + "Style of solution file Raw (computer-readable); " + "Pretty (human-readable): " + "0 => HiGHS raw; 1 => HiGHS pretty; 2 => Glpsol " + "raw; 3 => Glpsol pretty; ", advanced, &write_solution_style, kSolutionStyleMin, kSolutionStyleRaw, kSolutionStyleMax); records.push_back(record_int); + record_int = new OptionRecordInt( + "glpsol_cost_row_location", + "Location of cost row for Glpsol file: " + "-2 => Last; -1 => None; 0 => None if empty, otherwise data file " + "location; 1 <= n <= num_row => Location n; n > " + "num_row => Last", + advanced, &glpsol_cost_row_location, kGlpsolCostRowLocationMin, 0, + kHighsIInf); + records.push_back(record_int); + record_bool = new OptionRecordBool("icrash", "Run iCrash", advanced, &icrash, false); records.push_back(record_bool); diff --git a/src/lp_data/HighsSolution.cpp b/src/lp_data/HighsSolution.cpp index 60f626a80b..c5dc4770db 100644 --- a/src/lp_data/HighsSolution.cpp +++ b/src/lp_data/HighsSolution.cpp @@ -71,21 +71,38 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, double primal_feasibility_tolerance = options.primal_feasibility_tolerance; double dual_feasibility_tolerance = options.dual_feasibility_tolerance; // highs_info are the values computed in this method. + HighsInt& num_primal_infeasibility = highs_info.num_primal_infeasibilities; - double& max_primal_infeasibility = highs_info.max_primal_infeasibility; + + HighsInt& max_absolute_primal_infeasibility_index = + primal_dual_errors.max_primal_infeasibility.absolute_index; + double& max_absolute_primal_infeasibility_value = + highs_info.max_primal_infeasibility; + HighsInt& max_relative_primal_infeasibility_index = + primal_dual_errors.max_primal_infeasibility.relative_index; + double& max_relative_primal_infeasibility_value = + primal_dual_errors.max_primal_infeasibility.relative_value; + double& sum_primal_infeasibility = highs_info.sum_primal_infeasibilities; + HighsInt& num_dual_infeasibility = highs_info.num_dual_infeasibilities; - double& max_dual_infeasibility = highs_info.max_dual_infeasibility; + + HighsInt& max_dual_infeasibility_index = + primal_dual_errors.max_dual_infeasibility.absolute_index; + double& max_dual_infeasibility_value = highs_info.max_dual_infeasibility; + double& sum_dual_infeasibility = highs_info.sum_dual_infeasibilities; num_primal_infeasibility = kHighsIllegalInfeasibilityCount; - max_primal_infeasibility = kHighsIllegalInfeasibilityMeasure; + max_absolute_primal_infeasibility_value = kHighsIllegalInfeasibilityMeasure; sum_primal_infeasibility = kHighsIllegalInfeasibilityMeasure; + primal_dual_errors.max_primal_infeasibility.invalidate(); highs_info.primal_solution_status = kSolutionStatusNone; num_dual_infeasibility = kHighsIllegalInfeasibilityCount; - max_dual_infeasibility = kHighsIllegalInfeasibilityMeasure; + max_dual_infeasibility_value = kHighsIllegalInfeasibilityMeasure; sum_dual_infeasibility = kHighsIllegalInfeasibilityMeasure; + primal_dual_errors.max_dual_infeasibility.invalidate(); highs_info.dual_solution_status = kSolutionStatusNone; const bool& have_primal_solution = solution.value_valid; @@ -103,25 +120,45 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, assert((int)solution.col_value.size() >= lp.num_col_); assert((int)solution.row_value.size() >= lp.num_row_); num_primal_infeasibility = 0; - max_primal_infeasibility = 0; + max_absolute_primal_infeasibility_value = 0; sum_primal_infeasibility = 0; + primal_dual_errors.max_primal_infeasibility.reset(); if (have_dual_solution) { // There's a dual solution, so check its size and initialise the // infeasiblilty counts assert((int)solution.col_dual.size() >= lp.num_col_); assert((int)solution.row_dual.size() >= lp.num_row_); num_dual_infeasibility = 0; - max_dual_infeasibility = 0; + max_dual_infeasibility_value = 0; sum_dual_infeasibility = 0; + primal_dual_errors.max_dual_infeasibility.reset(); } } HighsInt& num_primal_residual = primal_dual_errors.num_primal_residual; - double& max_primal_residual = primal_dual_errors.max_primal_residual; + + HighsInt& max_absolute_primal_residual_index = + primal_dual_errors.max_primal_residual.absolute_index; + double& max_absolute_primal_residual_value = + primal_dual_errors.max_primal_residual.absolute_value; + HighsInt& max_relative_primal_residual_index = + primal_dual_errors.max_primal_residual.relative_index; + double& max_relative_primal_residual_value = + primal_dual_errors.max_primal_residual.relative_value; + double& sum_primal_residual = primal_dual_errors.sum_primal_residual; HighsInt& num_dual_residual = primal_dual_errors.num_dual_residual; - double& max_dual_residual = primal_dual_errors.max_dual_residual; + + HighsInt& max_absolute_dual_residual_index = + primal_dual_errors.max_dual_residual.absolute_index; + double& max_absolute_dual_residual_value = + primal_dual_errors.max_dual_residual.absolute_value; + HighsInt& max_relative_dual_residual_index = + primal_dual_errors.max_dual_residual.relative_index; + double& max_relative_dual_residual_value = + primal_dual_errors.max_dual_residual.relative_value; + double& sum_dual_residual = primal_dual_errors.sum_dual_residual; HighsInt& num_nonzero_basic_duals = @@ -139,21 +176,25 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, if (have_primal_solution && get_residuals) { num_primal_residual = 0; - max_primal_residual = 0; + max_absolute_primal_residual_value = 0; sum_primal_residual = 0; + primal_dual_errors.max_primal_residual.reset(); } else { num_primal_residual = kHighsIllegalInfeasibilityCount; - max_primal_residual = kHighsIllegalInfeasibilityMeasure; + max_absolute_primal_residual_value = kHighsIllegalInfeasibilityMeasure; sum_primal_residual = kHighsIllegalInfeasibilityMeasure; + primal_dual_errors.max_primal_residual.invalidate(); } if (have_dual_solution && get_residuals) { num_dual_residual = 0; - max_dual_residual = 0; + max_absolute_dual_residual_value = 0; sum_dual_residual = 0; + primal_dual_errors.max_dual_residual.reset(); } else { num_dual_residual = kHighsIllegalInfeasibilityCount; - max_dual_residual = kHighsIllegalInfeasibilityMeasure; + max_absolute_dual_residual_value = kHighsIllegalInfeasibilityMeasure; sum_dual_residual = kHighsIllegalInfeasibilityMeasure; + primal_dual_errors.max_dual_residual.invalidate(); } if (have_basis) { num_nonzero_basic_duals = 0; @@ -174,13 +215,20 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, max_off_bound_nonbasic = kHighsIllegalInfeasibilityMeasure; sum_off_bound_nonbasic = kHighsIllegalInfeasibilityMeasure; } + // Without a primal solution, nothing can be done! if (!have_primal_solution) return; - std::vector primal_activities; - std::vector dual_activities; + std::vector primal_positive_sum; + std::vector primal_negative_sum; + std::vector dual_positive_sum; + std::vector dual_negative_sum; if (get_residuals) { - primal_activities.assign(lp.num_row_, 0); - if (have_dual_solution) dual_activities.resize(lp.num_col_); + primal_positive_sum.assign(lp.num_row_, 0); + primal_negative_sum.assign(lp.num_row_, 0); + if (have_dual_solution) { + dual_positive_sum.resize(lp.num_col_); + dual_negative_sum.resize(lp.num_col_); + } } HighsInt num_basic_var = 0; HighsInt num_non_basic_var = 0; @@ -192,7 +240,8 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, // entry of basis since this is const. const HighsBasisStatus* status_pointer = have_basis ? &status : NULL; - double primal_infeasibility; + double absolute_primal_infeasibility; + double relative_primal_infeasibility; double dual_infeasibility; double value_residual; double lower; @@ -221,23 +270,33 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, } // Flip dual according to lp.sense_ dual *= (HighsInt)lp.sense_; - getVariableKktFailures( // const HighsLogOptions& log_options, + getVariableKktFailures( primal_feasibility_tolerance, dual_feasibility_tolerance, lower, upper, - value, dual, status_pointer, integrality, primal_infeasibility, - dual_infeasibility, value_residual); + value, dual, status_pointer, integrality, absolute_primal_infeasibility, + relative_primal_infeasibility, dual_infeasibility, value_residual); // Accumulate primal infeasiblilties - if (primal_infeasibility > primal_feasibility_tolerance) + if (absolute_primal_infeasibility > primal_feasibility_tolerance) num_primal_infeasibility++; - max_primal_infeasibility = - std::max(primal_infeasibility, max_primal_infeasibility); - sum_primal_infeasibility += primal_infeasibility; + if (max_absolute_primal_infeasibility_value < + absolute_primal_infeasibility) { + max_absolute_primal_infeasibility_index = iVar; + max_absolute_primal_infeasibility_value = absolute_primal_infeasibility; + } + if (max_relative_primal_infeasibility_value < + relative_primal_infeasibility) { + max_relative_primal_infeasibility_index = iVar; + max_relative_primal_infeasibility_value = relative_primal_infeasibility; + } + sum_primal_infeasibility += absolute_primal_infeasibility; if (have_dual_solution) { // Accumulate dual infeasiblilties if (dual_infeasibility > dual_feasibility_tolerance) num_dual_infeasibility++; - max_dual_infeasibility = - std::max(dual_infeasibility, max_dual_infeasibility); + if (max_dual_infeasibility_value < dual_infeasibility) { + max_dual_infeasibility_value = dual_infeasibility; + max_dual_infeasibility_index = iVar; + } sum_dual_infeasibility += dual_infeasibility; } if (have_basis) { @@ -263,35 +322,89 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, } if (iVar < lp.num_col_ && get_residuals) { HighsInt iCol = iVar; - if (have_dual_solution) dual_activities[iCol] = gradient[iCol]; + if (have_dual_solution) { + if (gradient[iCol] > 0) { + dual_positive_sum[iCol] = gradient[iCol]; + } else { + dual_negative_sum[iCol] = -gradient[iCol]; + } + } for (HighsInt el = lp.a_matrix_.start_[iCol]; el < lp.a_matrix_.start_[iCol + 1]; el++) { HighsInt iRow = lp.a_matrix_.index_[el]; double Avalue = lp.a_matrix_.value_[el]; - primal_activities[iRow] += value * Avalue; + double term = value * Avalue; + if (term > 0) { + primal_positive_sum[iRow] += term; + } else { + primal_negative_sum[iRow] -= term; + } // @FlipRowDual += became -= - if (have_dual_solution) - dual_activities[iCol] -= solution.row_dual[iRow] * Avalue; + if (have_dual_solution) { + double term = -solution.row_dual[iRow] * Avalue; + if (term > 0) { + dual_positive_sum[iCol] += term; + } else { + dual_negative_sum[iCol] -= term; + } + } } } } if (get_residuals) { const double large_residual_error = 1e-12; for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { - double primal_residual_error = - std::fabs(primal_activities[iRow] - solution.row_value[iRow]); - if (primal_residual_error > large_residual_error) num_primal_residual++; - max_primal_residual = - std::max(primal_residual_error, max_primal_residual); - sum_primal_residual += primal_residual_error; + double term = -solution.row_value[iRow]; + if (term > 0) { + primal_positive_sum[iRow] += term; + } else { + primal_negative_sum[iRow] -= term; + } + assert(primal_positive_sum[iRow] >= 0); + assert(primal_negative_sum[iRow] >= 0); + double absolute_primal_residual = + std::fabs(primal_positive_sum[iRow] - primal_negative_sum[iRow]); + double relative_primal_residual = + absolute_primal_residual / + (1 + primal_positive_sum[iRow] + primal_negative_sum[iRow]); + + if (absolute_primal_residual > large_residual_error) + num_primal_residual++; + if (max_absolute_primal_residual_value < absolute_primal_residual) { + max_absolute_primal_residual_value = absolute_primal_residual; + max_absolute_primal_residual_index = iRow; + } + if (max_relative_primal_residual_value < relative_primal_residual) { + max_relative_primal_residual_value = relative_primal_residual; + max_relative_primal_residual_index = iRow; + } + sum_primal_residual += absolute_primal_residual; } if (have_dual_solution) { for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { - double dual_residual_error = - std::fabs(dual_activities[iCol] - solution.col_dual[iCol]); - if (dual_residual_error > large_residual_error) num_dual_residual++; - max_dual_residual = std::max(dual_residual_error, max_dual_residual); - sum_dual_residual += dual_residual_error; + double term = -solution.col_dual[iCol]; + if (term > 0) { + dual_positive_sum[iCol] += term; + } else { + dual_negative_sum[iCol] -= term; + } + assert(dual_positive_sum[iCol] >= 0); + assert(dual_negative_sum[iCol] >= 0); + double absolute_dual_residual = + std::fabs(dual_positive_sum[iCol] - dual_negative_sum[iCol]); + double relative_dual_residual = + absolute_dual_residual / + (1 + dual_positive_sum[iCol] + dual_negative_sum[iCol]); + if (absolute_dual_residual > large_residual_error) num_dual_residual++; + if (max_absolute_dual_residual_value < absolute_dual_residual) { + max_absolute_dual_residual_value = absolute_dual_residual; + max_absolute_dual_residual_index = iCol; + } + if (max_relative_dual_residual_value < relative_dual_residual) { + max_relative_dual_residual_value = relative_dual_residual; + max_relative_dual_residual_index = iCol; + } + sum_dual_residual += absolute_dual_residual; } } } @@ -309,13 +422,26 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, highs_info.dual_solution_status = kSolutionStatusFeasible; } } + // Assign the two entries in primal_dual_errors that are accumulated + // in highs_info + primal_dual_errors.max_primal_infeasibility.absolute_value = + highs_info.max_primal_infeasibility; + ; + primal_dual_errors.max_dual_infeasibility.absolute_value = + highs_info.max_dual_infeasibility; + // Relative dual_infeasibility data is the same as absolute + primal_dual_errors.max_dual_infeasibility.relative_value = + primal_dual_errors.max_dual_infeasibility.absolute_value; + primal_dual_errors.max_dual_infeasibility.relative_index = + primal_dual_errors.max_dual_infeasibility.absolute_index; } + // Gets the KKT failures for a variable. The lack of a basis status is // indicated by status_pointer being null. // // Value and dual are used compute the primal and dual infeasibility - -// according to the basis status (if valid) or primal value.It's up to -// the calling method to ignore these if the value or dual are not +// according to the basis status (if valid) or primal value. It's up +// to the calling method to ignore these if the value or dual are not // valid. // // If the basis status is valid, then the numbers of basic and @@ -327,25 +453,33 @@ void getVariableKktFailures(const double primal_feasibility_tolerance, const double value, const double dual, const HighsBasisStatus* status_pointer, const HighsVarType integrality, - double& primal_infeasibility, + double& absolute_primal_infeasibility, + double& relative_primal_infeasibility, double& dual_infeasibility, double& value_residual) { const double middle = (lower + upper) * 0.5; // @primal_infeasibility calculation - primal_infeasibility = 0; + absolute_primal_infeasibility = 0; + relative_primal_infeasibility = 0; if (value < lower - primal_feasibility_tolerance) { // Below lower - primal_infeasibility = lower - value; + absolute_primal_infeasibility = lower - value; + relative_primal_infeasibility = + absolute_primal_infeasibility / (1 + std::fabs(lower)); } else if (value > upper + primal_feasibility_tolerance) { // Above upper - primal_infeasibility = value - upper; + absolute_primal_infeasibility = value - upper; + relative_primal_infeasibility = + absolute_primal_infeasibility / (1 + std::fabs(upper)); } // Account for semi-variables - if (primal_infeasibility > 0 && + if (absolute_primal_infeasibility > 0 && (integrality == HighsVarType::kSemiContinuous || integrality == HighsVarType::kSemiInteger) && - std::fabs(value) < primal_feasibility_tolerance) - primal_infeasibility = 0; + std::fabs(value) < primal_feasibility_tolerance) { + absolute_primal_infeasibility = 0; + relative_primal_infeasibility = 0; + } value_residual = std::min(std::fabs(lower - value), std::fabs(value - upper)); // Determine whether the variable is at a bound using the basis // status (if valid) or value residual. @@ -377,6 +511,28 @@ void getVariableKktFailures(const double primal_feasibility_tolerance, } } +void HighsError::print(std::string message) { + printf( + "\n%s\nAbsolute value = %11.4g; index = %9d\nRelative value = %11.4g; " + "index = %9d\n", + message.c_str(), this->absolute_value, (int)this->absolute_index, + this->relative_value, (int)this->relative_index); +} + +void HighsError::reset() { + this->absolute_value = 0; + this->absolute_index = 0; + this->relative_value = 0; + this->relative_index = 0; +} + +void HighsError::invalidate() { + this->absolute_value = kHighsIllegalErrorValue; + this->absolute_index = kHighsIllegalErrorIndex; + this->relative_value = kHighsIllegalErrorValue; + this->relative_index = kHighsIllegalErrorIndex; +} + double computeObjectiveValue(const HighsLp& lp, const HighsSolution& solution) { double objective_value = 0; for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) diff --git a/src/lp_data/HighsSolution.h b/src/lp_data/HighsSolution.h index 4f8c8c6c0c..b2137449cf 100644 --- a/src/lp_data/HighsSolution.h +++ b/src/lp_data/HighsSolution.h @@ -32,6 +32,16 @@ class HighsOptions; using std::string; +struct HighsError { + double absolute_value; + HighsInt absolute_index; + double relative_value; + HighsInt relative_index; + void print(std::string message); + void reset(); + void invalidate(); +}; + struct HighsPrimalDualErrors { HighsInt num_nonzero_basic_duals; HighsInt num_large_nonzero_basic_duals; @@ -41,11 +51,13 @@ struct HighsPrimalDualErrors { double max_off_bound_nonbasic; double sum_off_bound_nonbasic; HighsInt num_primal_residual; - double max_primal_residual; double sum_primal_residual; HighsInt num_dual_residual; - double max_dual_residual; double sum_dual_residual; + HighsError max_primal_residual; + HighsError max_primal_infeasibility; + HighsError max_dual_residual; + HighsError max_dual_infeasibility; }; void getKktFailures(const HighsOptions& options, const HighsModel& model, @@ -81,7 +93,8 @@ void getVariableKktFailures(const double primal_feasibility_tolerance, const double value, const double dual, const HighsBasisStatus* status_pointer, const HighsVarType integrality, - double& primal_infeasibility, + double& absolute_primal_infeasibility, + double& relative_primal_infeasibility, double& dual_infeasibility, double& value_residual); double computeObjectiveValue(const HighsLp& lp, const HighsSolution& solution); diff --git a/src/lp_data/HighsSolutionDebug.cpp b/src/lp_data/HighsSolutionDebug.cpp index 703ff65ecc..71eeb13b56 100644 --- a/src/lp_data/HighsSolutionDebug.cpp +++ b/src/lp_data/HighsSolutionDebug.cpp @@ -292,7 +292,7 @@ HighsDebugStatus debugAnalysePrimalDualErrors( if (force_report) report_level = HighsLogType::kInfo; highsLogDev( options.log_options, report_level, - "PrDuErrors : %-9s Nonzero basic duals: num = %2" HIGHSINT_FORMAT + "PrDuErrors : %-9s Nonzero basic duals: num = %7" HIGHSINT_FORMAT "; " "max = %9.4g; sum = %9.4g\n", value_adjective.c_str(), primal_dual_errors.num_nonzero_basic_duals, @@ -312,7 +312,7 @@ HighsDebugStatus debugAnalysePrimalDualErrors( if (force_report) report_level = HighsLogType::kInfo; highsLogDev( options.log_options, report_level, - "PrDuErrors : %-9s Off-bound nonbasic values: num = %2" HIGHSINT_FORMAT + "PrDuErrors : %-9s Off-bound nonbasic values: num = %7" HIGHSINT_FORMAT "; " "max = %9.4g; sum = %9.4g\n", value_adjective.c_str(), primal_dual_errors.num_off_bound_nonbasic, @@ -320,11 +320,13 @@ HighsDebugStatus debugAnalysePrimalDualErrors( primal_dual_errors.sum_off_bound_nonbasic); } if (primal_dual_errors.num_primal_residual >= 0) { - if (primal_dual_errors.max_primal_residual > excessive_residual_error) { + if (primal_dual_errors.max_primal_residual.absolute_value > + excessive_residual_error) { value_adjective = "Excessive"; report_level = HighsLogType::kError; return_status = HighsDebugStatus::kError; - } else if (primal_dual_errors.max_primal_residual > large_residual_error) { + } else if (primal_dual_errors.max_primal_residual.absolute_value > + large_residual_error) { value_adjective = "Large"; report_level = HighsLogType::kDetailed; return_status = HighsDebugStatus::kWarning; @@ -336,19 +338,21 @@ HighsDebugStatus debugAnalysePrimalDualErrors( if (force_report) report_level = HighsLogType::kInfo; highsLogDev( options.log_options, report_level, - "PrDuErrors : %-9s Primal residual: num = %2" HIGHSINT_FORMAT + "PrDuErrors : %-9s Primal residual: num = %7" HIGHSINT_FORMAT "; " "max = %9.4g; sum = %9.4g\n", value_adjective.c_str(), primal_dual_errors.num_primal_residual, - primal_dual_errors.max_primal_residual, + primal_dual_errors.max_primal_residual.absolute_value, primal_dual_errors.sum_primal_residual); } if (primal_dual_errors.num_dual_residual >= 0) { - if (primal_dual_errors.max_dual_residual > excessive_residual_error) { + if (primal_dual_errors.max_dual_residual.absolute_value > + excessive_residual_error) { value_adjective = "Excessive"; report_level = HighsLogType::kError; return_status = HighsDebugStatus::kError; - } else if (primal_dual_errors.max_dual_residual > large_residual_error) { + } else if (primal_dual_errors.max_dual_residual.absolute_value > + large_residual_error) { value_adjective = "Large"; report_level = HighsLogType::kDetailed; return_status = HighsDebugStatus::kWarning; @@ -360,11 +364,11 @@ HighsDebugStatus debugAnalysePrimalDualErrors( if (force_report) report_level = HighsLogType::kInfo; highsLogDev( options.log_options, report_level, - "PrDuErrors : %-9s Dual residual: num = %2" HIGHSINT_FORMAT + "PrDuErrors : %-9s Dual residual: num = %7" HIGHSINT_FORMAT "; " "max = %9.4g; sum = %9.4g\n", value_adjective.c_str(), primal_dual_errors.num_dual_residual, - primal_dual_errors.max_dual_residual, + primal_dual_errors.max_dual_residual.absolute_value, primal_dual_errors.sum_dual_residual); } return return_status; diff --git a/src/model/HighsModel.cpp b/src/model/HighsModel.cpp index 1a115ae414..d717954dcb 100644 --- a/src/model/HighsModel.cpp +++ b/src/model/HighsModel.cpp @@ -34,6 +34,6 @@ void HighsModel::objectiveGradient(const std::vector& solution, } else { gradient.assign(this->lp_.num_col_, 0); } - for (HighsInt iCol = 0; iCol < this->hessian_.dim_; iCol++) + for (HighsInt iCol = 0; iCol < this->lp_.num_col_; iCol++) gradient[iCol] += this->lp_.col_cost_[iCol]; }