1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the HiGHS linear optimization suite */
4 /* */
5 /* Written and engineered 2008-2021 at the University of Edinburgh */
6 /* */
7 /* Available as open-source under the MIT License */
8 /* */
9 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10 /**@file lp_data/Highs.cpp
11 * @brief
12 * @author Julian Hall, Ivet Galabova, Qi Huangfu and Michael Feldmeier
13 */
14 #include "Highs.h"
15
16 #include <algorithm>
17 #include <iostream>
18 #include <memory>
19 #include <sstream>
20
21 #include "HConfig.h"
22 #include "io/Filereader.h"
23 #include "io/HighsIO.h"
24 #include "io/LoadOptions.h"
25 #include "lp_data/HighsLpUtils.h"
26 #include "lp_data/HighsModelUtils.h"
27 #include "lp_data/HighsSolution.h"
28 #include "lp_data/HighsSolve.h"
29 #include "simplex/HSimplexDebug.h"
30 #include "simplex/HighsSimplexInterface.h"
31 #include "util/HighsMatrixPic.h"
32
33 #ifdef OPENMP
34 #include "omp.h"
35 #endif
36
Highs()37 Highs::Highs() {
38 hmos_.clear();
39 hmos_.push_back(HighsModelObject(lp_, options_, timer_));
40 }
41
setHighsOptionValue(const std::string & option,const bool value)42 HighsStatus Highs::setHighsOptionValue(const std::string& option,
43 const bool value) {
44 if (setOptionValue(options_.logfile, option, options_.records, value) ==
45 OptionStatus::OK)
46 return HighsStatus::OK;
47 return HighsStatus::Error;
48 }
49
setHighsOptionValue(const std::string & option,const int value)50 HighsStatus Highs::setHighsOptionValue(const std::string& option,
51 const int value) {
52 if (setOptionValue(options_.logfile, option, options_.records, value) ==
53 OptionStatus::OK)
54 return HighsStatus::OK;
55 return HighsStatus::Error;
56 }
57
setHighsOptionValue(const std::string & option,const double value)58 HighsStatus Highs::setHighsOptionValue(const std::string& option,
59 const double value) {
60 if (setOptionValue(options_.logfile, option, options_.records, value) ==
61 OptionStatus::OK)
62 return HighsStatus::OK;
63 return HighsStatus::Error;
64 }
65
setHighsOptionValue(const std::string & option,const std::string value)66 HighsStatus Highs::setHighsOptionValue(const std::string& option,
67 const std::string value) {
68 if (setOptionValue(options_.logfile, option, options_.records, value) ==
69 OptionStatus::OK)
70 return HighsStatus::OK;
71 return HighsStatus::Error;
72 }
73
setHighsOptionValue(const std::string & option,const char * value)74 HighsStatus Highs::setHighsOptionValue(const std::string& option,
75 const char* value) {
76 if (setOptionValue(options_.logfile, option, options_.records, value) ==
77 OptionStatus::OK)
78 return HighsStatus::OK;
79 return HighsStatus::Error;
80 }
81
setHighsLogfile(FILE * logfile)82 HighsStatus Highs::setHighsLogfile(FILE* logfile) {
83 options_.logfile = logfile;
84 return HighsStatus::OK;
85 }
86
setHighsOutput(FILE * output)87 HighsStatus Highs::setHighsOutput(FILE* output) {
88 options_.output = output;
89 return HighsStatus::OK;
90 }
91
readHighsOptions(const std::string filename)92 HighsStatus Highs::readHighsOptions(const std::string filename) {
93 if (filename.size() <= 0) {
94 HighsLogMessage(options_.logfile, HighsMessageType::WARNING,
95 "Empty file name so not reading options");
96 return HighsStatus::Warning;
97 }
98 options_.options_file = filename;
99 if (!loadOptionsFromFile(options_)) return HighsStatus::Error;
100 return HighsStatus::OK;
101 }
102
passHighsOptions(const HighsOptions & options)103 HighsStatus Highs::passHighsOptions(const HighsOptions& options) {
104 if (passOptions(options_.logfile, options, options_) == OptionStatus::OK)
105 return HighsStatus::OK;
106 return HighsStatus::Error;
107 }
108
getHighsOptions()109 const HighsOptions& Highs::getHighsOptions() { return options_; }
110
getHighsOptionValue(const std::string & option,bool & value)111 HighsStatus Highs::getHighsOptionValue(const std::string& option, bool& value) {
112 if (getOptionValue(options_.logfile, option, options_.records, value) ==
113 OptionStatus::OK)
114 return HighsStatus::OK;
115 return HighsStatus::Error;
116 }
117
getHighsOptionValue(const std::string & option,int & value)118 HighsStatus Highs::getHighsOptionValue(const std::string& option, int& value) {
119 if (getOptionValue(options_.logfile, option, options_.records, value) ==
120 OptionStatus::OK)
121 return HighsStatus::OK;
122 return HighsStatus::Error;
123 }
124
getHighsOptionValue(const std::string & option,double & value)125 HighsStatus Highs::getHighsOptionValue(const std::string& option,
126 double& value) {
127 if (getOptionValue(options_.logfile, option, options_.records, value) ==
128 OptionStatus::OK)
129 return HighsStatus::OK;
130 return HighsStatus::Error;
131 }
132
getHighsOptionValue(const std::string & option,std::string & value)133 HighsStatus Highs::getHighsOptionValue(const std::string& option,
134 std::string& value) {
135 if (getOptionValue(options_.logfile, option, options_.records, value) ==
136 OptionStatus::OK)
137 return HighsStatus::OK;
138 return HighsStatus::Error;
139 }
140
getHighsOptionType(const std::string & option,HighsOptionType & type)141 HighsStatus Highs::getHighsOptionType(const std::string& option,
142 HighsOptionType& type) {
143 if (getOptionType(options_.logfile, option, options_.records, type) ==
144 OptionStatus::OK)
145 return HighsStatus::OK;
146 return HighsStatus::Error;
147 }
148
resetHighsOptions()149 HighsStatus Highs::resetHighsOptions() {
150 resetOptions(options_.records);
151 return HighsStatus::OK;
152 }
153
writeHighsOptions(const std::string filename,const bool report_only_non_default_values)154 HighsStatus Highs::writeHighsOptions(
155 const std::string filename, const bool report_only_non_default_values) {
156 HighsStatus return_status = HighsStatus::OK;
157 HighsLp lp = this->lp_;
158 FILE* file;
159 bool html;
160 return_status = interpretCallStatus(
161 openWriteFile(filename, "writeHighsOptions", file, html), return_status,
162 "openWriteFile");
163 if (return_status == HighsStatus::Error) return return_status;
164
165 return_status = interpretCallStatus(
166 writeOptionsToFile(file, options_.records, report_only_non_default_values,
167 html),
168 return_status, "writeOptionsToFile");
169 return return_status;
170 }
171
getHighsOptions() const172 const HighsOptions& Highs::getHighsOptions() const { return options_; }
173
getHighsInfo() const174 const HighsInfo& Highs::getHighsInfo() const { return info_; }
175
getHighsInfoValue(const std::string & info,int & value)176 HighsStatus Highs::getHighsInfoValue(const std::string& info, int& value) {
177 if (getInfoValue(options_, info, info_.records, value) == InfoStatus::OK)
178 return HighsStatus::OK;
179 return HighsStatus::Error;
180 }
181
getHighsInfoValue(const std::string & info,double & value) const182 HighsStatus Highs::getHighsInfoValue(const std::string& info,
183 double& value) const {
184 if (getInfoValue(options_, info, info_.records, value) == InfoStatus::OK)
185 return HighsStatus::OK;
186 return HighsStatus::Error;
187 }
188
writeHighsInfo(const std::string filename)189 HighsStatus Highs::writeHighsInfo(const std::string filename) {
190 HighsStatus return_status = HighsStatus::OK;
191 HighsLp lp = this->lp_;
192 FILE* file;
193 bool html;
194 return_status =
195 interpretCallStatus(openWriteFile(filename, "writeHighsInfo", file, html),
196 return_status, "openWriteFile");
197 if (return_status == HighsStatus::Error) return return_status;
198
199 return_status =
200 interpretCallStatus(writeInfoToFile(file, info_.records, html),
201 return_status, "writeInfoToFile");
202 return return_status;
203 }
204
205 // Methods below change the incumbent model or solver infomation
206 // associated with it. Hence returnFromHighs is called at the end of
207 // each
reset()208 HighsStatus Highs::reset() {
209 HighsStatus return_status = HighsStatus::OK;
210 // Clear the status, solution, basis and info associated with any previous
211 // model
212 return_status =
213 interpretCallStatus(clearSolver(), return_status, "clearSolver");
214 if (return_status == HighsStatus::Error) return return_status;
215 // Clear any HiGHS model object
216 hmos_.clear();
217 // Create a HiGHS model object for this LP
218 hmos_.push_back(HighsModelObject(lp_, options_, timer_));
219
220 presolve_.clear();
221
222 return returnFromHighs(return_status);
223 }
224
passModel(HighsLp lp)225 HighsStatus Highs::passModel(HighsLp lp) {
226 HighsStatus return_status = HighsStatus::OK;
227 // move the copy of the LP to the internal LP
228 lp_ = std::move(lp);
229 // Check validity of the LP, normalising its values
230 return_status =
231 interpretCallStatus(assessLp(lp_, options_), return_status, "assessLp");
232 if (return_status == HighsStatus::Error) return return_status;
233 // Clear solver status, solution, basis and info associated with any
234 // previous model; clear any HiGHS model object; create a HiGHS
235 // model object for this LP
236 return_status = interpretCallStatus(reset(), return_status, "reset");
237 return returnFromHighs(return_status);
238 }
239
passModel(const int num_col,const int num_row,const int num_nz,const double * costs,const double * col_lower,const double * col_upper,const double * row_lower,const double * row_upper,const int * astart,const int * aindex,const double * avalue)240 HighsStatus Highs::passModel(const int num_col, const int num_row,
241 const int num_nz, const double* costs,
242 const double* col_lower, const double* col_upper,
243 const double* row_lower, const double* row_upper,
244 const int* astart, const int* aindex,
245 const double* avalue) {
246 HighsLp lp;
247 lp.numCol_ = num_col;
248 lp.numRow_ = num_row;
249 if (num_col > 0) {
250 assert(costs != NULL);
251 assert(col_lower != NULL);
252 assert(col_upper != NULL);
253 lp.colCost_.assign(costs, costs + num_col);
254 lp.colLower_.assign(col_lower, col_lower + num_col);
255 lp.colUpper_.assign(col_upper, col_upper + num_col);
256 }
257 if (num_row > 0) {
258 assert(row_lower != NULL);
259 assert(row_upper != NULL);
260 lp.rowLower_.assign(row_lower, row_lower + num_row);
261 lp.rowUpper_.assign(row_upper, row_upper + num_row);
262 }
263 if (num_nz > 0) {
264 assert(num_col > 0);
265 assert(num_row > 0);
266 assert(astart != NULL);
267 assert(aindex != NULL);
268 assert(avalue != NULL);
269 lp.Astart_.assign(astart, astart + num_col);
270 lp.Aindex_.assign(aindex, aindex + num_nz);
271 lp.Avalue_.assign(avalue, avalue + num_nz);
272 }
273 lp.Astart_.resize(num_col + 1);
274 lp.Astart_[num_col] = num_nz;
275 return this->passModel(std::move(lp));
276 }
277
readModel(const std::string filename)278 HighsStatus Highs::readModel(const std::string filename) {
279 HighsStatus return_status = HighsStatus::OK;
280 Filereader* reader = Filereader::getFilereader(filename);
281 if (reader == NULL) {
282 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
283 "Model file %s not supported", filename.c_str());
284 return HighsStatus::Error;
285 }
286
287 HighsLp model;
288 this->options_.model_file = filename;
289
290 FilereaderRetcode call_code =
291 reader->readModelFromFile(this->options_, model);
292 delete reader;
293 if (call_code != FilereaderRetcode::OK) {
294 interpretFilereaderRetcode(this->options_.logfile, filename.c_str(),
295 call_code);
296 return_status = interpretCallStatus(HighsStatus::Error, return_status,
297 "readModelFromFile");
298 if (return_status == HighsStatus::Error) return return_status;
299 }
300 model.model_name_ = extractModelName(filename);
301 return_status =
302 interpretCallStatus(this->passModel(model), return_status, "passModel");
303 return returnFromHighs(return_status);
304 }
305
clearModel()306 HighsStatus Highs::clearModel() {
307 HighsStatus return_status = HighsStatus::OK;
308 // Remove all HighsModelObject entries
309 hmos_.clear();
310 // Set up with an empty LP so that addrows/cols can be used to build
311 // HighsLp empty_lp; lp_ = empty_lp;
312 lp_.clear();
313 hmos_.push_back(HighsModelObject(lp_, options_, timer_));
314 return_status =
315 interpretCallStatus(this->clearSolver(), return_status, "clearSolver");
316 if (return_status == HighsStatus::Error) return return_status;
317 return returnFromHighs(return_status);
318 }
319
readBasis(const std::string filename)320 HighsStatus Highs::readBasis(const std::string filename) {
321 HighsStatus return_status = HighsStatus::OK;
322 // Try to read basis file into read_basis
323 HighsBasis read_basis = this->basis_;
324 return_status =
325 interpretCallStatus(readBasisFile(options_, read_basis, filename),
326 return_status, "readBasis");
327 if (return_status != HighsStatus::OK) return return_status;
328 // Basis read OK: check whether it's consistent with the LP
329 if (!isBasisConsistent(lp_, read_basis)) {
330 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
331 "readBasis: invalid basis");
332 return HighsStatus::Error;
333 }
334 // Update the HiGHS basis and invalidate any simplex basis for the model
335 this->basis_ = read_basis;
336 this->basis_.valid_ = true;
337 if (hmos_.size() > 0) {
338 HighsSimplexInterface interface(hmos_[0]);
339 interface.clearBasis();
340 }
341 // Can't use returnFromHighs since...
342 return HighsStatus::OK;
343 }
344
writeModel(const std::string filename)345 HighsStatus Highs::writeModel(const std::string filename) {
346 HighsStatus return_status = HighsStatus::OK;
347 HighsLp model = this->lp_;
348
349 if (filename == "") {
350 // Empty file name: report model on stdout
351 reportLp(options_, model, 2);
352 return_status = HighsStatus::OK;
353 } else {
354 Filereader* writer = Filereader::getFilereader(filename);
355 if (writer == NULL) {
356 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
357 "Model file %s not supported", filename.c_str());
358 return HighsStatus::Error;
359 }
360 return_status =
361 interpretCallStatus(writer->writeModelToFile(options_, filename, model),
362 return_status, "writeModelToFile");
363 delete writer;
364 }
365 return returnFromHighs(return_status);
366 }
367
writeBasis(const std::string filename)368 HighsStatus Highs::writeBasis(const std::string filename) {
369 HighsStatus return_status = HighsStatus::OK;
370 return_status =
371 interpretCallStatus(writeBasisFile(options_, this->basis_, filename),
372 return_status, "writeBasis");
373 return returnFromHighs(return_status);
374 }
375
376 // Checks the options calls presolve and postsolve if needed. Solvers are called
377 // with runLpSolver(..)
run()378 HighsStatus Highs::run() {
379 #ifdef HiGHSDEV
380 const int min_highs_debug_level = HIGHS_DEBUG_LEVEL_MIN;
381 // HIGHS_DEBUG_LEVEL_CHEAP;
382 // HIGHS_DEBUG_LEVEL_COSTLY;
383 // HIGHS_DEBUG_LEVEL_MAX;
384 if (options_.highs_debug_level < min_highs_debug_level) {
385 printf(
386 "Highs::run() HiGHSDEV define so switching options_.highs_debug_level "
387 "from %d to %d\n",
388 options_.highs_debug_level, min_highs_debug_level);
389 options_.highs_debug_level = min_highs_debug_level;
390 }
391 writeModel("HighsRunModel.mps");
392 // if (lp_.numRow_>0 && lp_.numCol_>0) writeLpMatrixPicToFile(options_,
393 // "LpMatrix", lp_);
394 #endif
395
396 #ifdef OPENMP
397 omp_max_threads = omp_get_max_threads();
398 assert(omp_max_threads > 0);
399 #ifdef HiGHSDEV
400 if (omp_max_threads <= 0)
401 printf("WARNING: omp_get_max_threads() returns %d\n", omp_max_threads);
402 printf("Running with %d OMP thread(s)\n", omp_max_threads);
403 #endif
404 #endif
405 HighsStatus return_status = HighsStatus::OK;
406 HighsStatus call_status;
407 // Zero the HiGHS iteration counts
408 zeroHighsIterationCounts(info_);
409 /*
410 if (options_.message_level >= 0) {
411 printf("\n!! Actually solving an LP with %d cols, %d rows", lp_.numCol_,
412 lp_.numRow_); if (lp_.numCol_) printf(" and %d nonzeros",
413 lp_.Astart_[lp_.numCol_]); printf(":basis.valid_ = %d: basis_.valid_ = %d:
414 simplex_lp_status_.has_basis = %d!!\n\n", basis_.valid_, hmos_[0].basis_.valid_,
415 hmos_[0].simplex_lp_status_.has_basis);
416 if (basis_.valid_ != hmos_[0].basis_.valid_) {
417 printf("NB %d = basis_.valid_ != hmos_[0].basis_.valid_ = %d\n",
418 basis_.valid_, hmos_[0].basis_.valid_);
419 }
420 }
421 */
422 // Determine whether a model has been loaded.
423 assert((int)hmos_.size() <= 1);
424 if (hmos_.size() == 0) {
425 // No Highs model object, so load model according to value of
426 // model_file
427 if (options_.model_file.compare(FILENAME_DEFAULT) == 0) {
428 // model_file is still default value, so return with error
429 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
430 "No model can be loaded in run()");
431 return_status = HighsStatus::Error;
432 return returnFromRun(return_status);
433 } else {
434 std::string model_file = options_.model_file;
435 call_status = readModel(model_file);
436 return_status =
437 interpretCallStatus(call_status, return_status, "readModel");
438 if (return_status == HighsStatus::Error)
439 return returnFromRun(return_status);
440 }
441 }
442 // Ensure that there is exactly one Highs model object
443 assert((int)hmos_.size() == 1);
444
445 // Initialise the HiGHS model status values
446 hmos_[0].scaled_model_status_ = HighsModelStatus::NOTSET;
447 hmos_[0].unscaled_model_status_ = HighsModelStatus::NOTSET;
448 model_status_ = hmos_[0].scaled_model_status_;
449 scaled_model_status_ = hmos_[0].unscaled_model_status_;
450
451 #ifdef HIGHSDEV
452 // Shouldn't have to check validity of the LP since this is done when it is
453 // loaded or modified
454 call_status = assessLp(lp_, options_);
455 // If any errors have been found or normalisation carried out,
456 // call_status will be ERROR or WARNING, so only valid return is OK.
457 assert(call_status == HighsStatus::OK);
458 return_status = interpretCallStatus(call_status, return_status, "assessLp");
459 if (return_status == HighsStatus::Error) return returnFromRun(return_status);
460 #endif
461
462 // Return immediately if the LP has no columns
463 if (!lp_.numCol_) {
464 model_status_ = HighsModelStatus::MODEL_EMPTY;
465 scaled_model_status_ = model_status_;
466 hmos_[0].unscaled_model_status_ = model_status_;
467 hmos_[0].scaled_model_status_ = model_status_;
468 return_status = highsStatusFromHighsModelStatus(model_status_);
469 return returnFromRun(return_status);
470 }
471
472 HighsSetIO(options_);
473 #ifdef HiGHSDEV
474 if (checkOptions(options_.logfile, options_.records) != OptionStatus::OK) {
475 return_status = HighsStatus::Error;
476 return returnFromRun(return_status);
477 }
478 #endif
479 HighsPrintMessage(options_.output, options_.message_level, ML_VERBOSE,
480 "Solving %s\n", lp_.model_name_.c_str());
481
482 double this_presolve_time = -1;
483 double this_solve_presolved_lp_time = -1;
484 double this_postsolve_time = -1;
485 double this_solve_original_lp_time = -1;
486
487 // Running as LP solver: start the HiGHS clock unless it's already running
488 bool run_highs_clock_already_running = timer_.runningRunHighsClock();
489 if (!run_highs_clock_already_running) timer_.startRunHighsClock();
490 // Record the initial time and set the postsolve iteration count to
491 // -1 to identify whether it's not required
492 double initial_time = timer_.readRunHighsClock();
493 int postsolve_iteration_count = -1;
494 // Define identifiers to refer to the HMO of the original LP (0) and
495 // the HMO created when using presolve. The index of this HMO is 1
496 // when solving a one-off LP, but greater than one if presolve has
497 // been called multiple times. It's equal to the size of HMO
498 const int original_hmo = 0;
499 const int presolve_hmo = hmos_.size();
500 // Keep track of the hmo that is the most recently solved. By default it's the
501 // original LP
502 int solved_hmo = original_hmo;
503
504 if (!basis_.valid_ && options_.presolve != off_string) {
505 // No HiGHS basis so consider presolve
506 //
507 // If using IPX to solve the reduced LP, crossover must be run
508 // since a basic solution is required by postsolve
509 if (options_.solver == ipm_string && !options_.run_crossover) {
510 HighsLogMessage(options_.logfile, HighsMessageType::WARNING,
511 "Forcing IPX to use crossover after presolve");
512 options_.run_crossover = true;
513 }
514
515 hmos_[original_hmo].scaled_model_status_ = HighsModelStatus::NOTSET;
516 // Presolve. runPresolve handles the level of presolving (0 = don't
517 // presolve).
518
519 // printf("Writing before_presolve.mps\n");
520 // writeModel("before_presolve.mps");
521
522 // Run and time presolve.
523 const double from_presolve_time = timer_.read(timer_.presolve_clock);
524 this_presolve_time = -from_presolve_time;
525 timer_.start(timer_.presolve_clock);
526
527 HighsPresolveStatus presolve_status = runPresolve();
528 timer_.stop(timer_.presolve_clock);
529 const double to_presolve_time = timer_.read(timer_.presolve_clock);
530 this_presolve_time += to_presolve_time;
531 presolve_.info_.presolve_time = this_presolve_time;
532
533 // Set an illegal local pivot threshold value that's updated after
534 // solving the presolved LP - if simplex is used
535 double factor_pivot_threshold = -1;
536
537 // Run solver.
538 switch (presolve_status) {
539 case HighsPresolveStatus::NotPresolved: {
540 hmos_[solved_hmo].lp_.lp_name_ = "Original LP";
541 this_solve_original_lp_time = -timer_.read(timer_.solve_clock);
542 timer_.start(timer_.solve_clock);
543 call_status = runLpSolver(solved_hmo, "Not presolved: solving the LP");
544 timer_.stop(timer_.solve_clock);
545 this_solve_original_lp_time += timer_.read(timer_.solve_clock);
546 return_status =
547 interpretCallStatus(call_status, return_status, "runLpSolver");
548 if (return_status == HighsStatus::Error)
549 return returnFromRun(return_status);
550 break;
551 }
552 case HighsPresolveStatus::NotReduced: {
553 hmos_[solved_hmo].lp_.lp_name_ = "Unreduced LP";
554 // Log the presolve reductions
555 reportPresolveReductions(hmos_[original_hmo].options_,
556 hmos_[original_hmo].lp_, false);
557 this_solve_original_lp_time = -timer_.read(timer_.solve_clock);
558 timer_.start(timer_.solve_clock);
559 call_status = runLpSolver(
560 solved_hmo, "Problem not reduced by presolve: solving the LP");
561 timer_.stop(timer_.solve_clock);
562 this_solve_original_lp_time += timer_.read(timer_.solve_clock);
563 return_status =
564 interpretCallStatus(call_status, return_status, "runLpSolver");
565 if (return_status == HighsStatus::Error)
566 return returnFromRun(return_status);
567 break;
568 }
569 case HighsPresolveStatus::Reduced: {
570 HighsLp& reduced_lp = presolve_.getReducedProblem();
571 // Validate the reduced LP
572 assert(assessLp(reduced_lp, options_) == HighsStatus::OK);
573 call_status = cleanBounds(options_, reduced_lp);
574 // Ignore any warning from clean bounds since the original LP
575 // is still solved after presolve
576 if (interpretCallStatus(call_status, return_status, "cleanBounds") ==
577 HighsStatus::Error)
578 return HighsStatus::Error;
579 // Add reduced lp object to vector of HighsModelObject,
580 // so the last one in lp_ is the presolved one.
581
582 hmos_.push_back(HighsModelObject(reduced_lp, options_, timer_));
583 // Log the presolve reductions
584 reportPresolveReductions(hmos_[original_hmo].options_,
585 hmos_[original_hmo].lp_,
586 hmos_[presolve_hmo].lp_);
587 // Record the HMO to be solved
588 solved_hmo = presolve_hmo;
589 hmos_[solved_hmo].lp_.lp_name_ = "Presolved LP";
590 // Don't try dual cut-off when solving the presolved LP, as the
591 // objective values aren't correct
592 // HighsOptions& options = hmos_[solved_hmo].options_;
593 // HighsOptions save_options = options;
594 const double save_dual_objective_value_upper_bound =
595 options_.dual_objective_value_upper_bound;
596 options_.dual_objective_value_upper_bound = HIGHS_CONST_INF;
597 this_solve_presolved_lp_time = -timer_.read(timer_.solve_clock);
598 timer_.start(timer_.solve_clock);
599 call_status = runLpSolver(solved_hmo, "Solving the presolved LP");
600 timer_.stop(timer_.solve_clock);
601 this_solve_presolved_lp_time += timer_.read(timer_.solve_clock);
602 if (hmos_[solved_hmo].simplex_lp_status_.valid) {
603 // Record the pivot threshold resulting from solving the presolved LP
604 // with simplex
605 factor_pivot_threshold =
606 hmos_[solved_hmo].simplex_info_.factor_pivot_threshold;
607 }
608 // Restore the dual objective cut-off
609 options_.dual_objective_value_upper_bound =
610 save_dual_objective_value_upper_bound;
611 return_status =
612 interpretCallStatus(call_status, return_status, "runLpSolver");
613 if (return_status == HighsStatus::Error)
614 return returnFromRun(return_status);
615
616 break;
617 }
618 case HighsPresolveStatus::ReducedToEmpty: {
619 reportPresolveReductions(hmos_[original_hmo].options_,
620 hmos_[original_hmo].lp_, true);
621 hmos_[original_hmo].unscaled_model_status_ = HighsModelStatus::OPTIMAL;
622 hmos_[original_hmo].scaled_model_status_ =
623 hmos_[original_hmo].unscaled_model_status_;
624 // Proceed to postsolve.
625 break;
626 }
627 // printf("\nHighs::run() 3: presolve status = %d\n",
628 //(int)presolve_status);fflush(stdout);
629 case HighsPresolveStatus::Infeasible:
630 case HighsPresolveStatus::Unbounded: {
631 if (presolve_status == HighsPresolveStatus::Infeasible) {
632 model_status_ = HighsModelStatus::PRIMAL_INFEASIBLE;
633 } else {
634 model_status_ = HighsModelStatus::PRIMAL_UNBOUNDED;
635 }
636 HighsLogMessage(options_.logfile, HighsMessageType::INFO,
637 "Problem status detected on presolve: %s",
638 highsModelStatusToString(model_status_).c_str());
639
640 // Report this way for the moment. May modify after merge with
641 // OSIinterface branch which has new way of setting up a
642 // HighsModelObject and can support multiple calls to run(). Stop and
643 // read the HiGHS clock, then work out time for this call
644 if (!run_highs_clock_already_running) timer_.stopRunHighsClock();
645
646 // Transfer the model status to the scaled model status and orriginal
647 // HMO statuses;
648 scaled_model_status_ = model_status_;
649 hmos_[original_hmo].unscaled_model_status_ = model_status_;
650 hmos_[original_hmo].scaled_model_status_ = model_status_;
651 return_status = HighsStatus::OK;
652 return returnFromRun(return_status);
653 }
654 case HighsPresolveStatus::Timeout: {
655 model_status_ = HighsModelStatus::PRESOLVE_ERROR;
656 HighsPrintMessage(options_.output, options_.message_level, ML_ALWAYS,
657 "Presolve reached timeout\n");
658 if (run_highs_clock_already_running) timer_.stopRunHighsClock();
659 return HighsStatus::Warning;
660 }
661 case HighsPresolveStatus::OptionsError: {
662 model_status_ = HighsModelStatus::PRESOLVE_ERROR;
663 HighsPrintMessage(options_.output, options_.message_level, ML_ALWAYS,
664 "Presolve options error.\n");
665 if (run_highs_clock_already_running) timer_.stopRunHighsClock();
666 return HighsStatus::Warning;
667 }
668 default: {
669 // case HighsPresolveStatus::Error
670 model_status_ = HighsModelStatus::PRESOLVE_ERROR;
671 HighsPrintMessage(options_.output, options_.message_level, ML_ALWAYS,
672 "Presolve failed.\n");
673 if (run_highs_clock_already_running) timer_.stopRunHighsClock();
674 // Transfer the model status to the scaled model status and orriginal
675 // HMO statuses;
676 scaled_model_status_ = model_status_;
677 hmos_[original_hmo].unscaled_model_status_ = model_status_;
678 hmos_[original_hmo].scaled_model_status_ = model_status_;
679 return_status = HighsStatus::Error;
680 return returnFromRun(return_status);
681 }
682 }
683 // Postsolve. Does nothing if there were no reductions during presolve.
684 if (hmos_[solved_hmo].scaled_model_status_ == HighsModelStatus::OPTIMAL) {
685 if (presolve_status == HighsPresolveStatus::Reduced ||
686 presolve_status == HighsPresolveStatus::ReducedToEmpty) {
687 // If presolve is nontrivial, extract the optimal solution
688 // and basis for the presolved problem in order to generate
689 // the solution and basis for postsolve to use to generate a
690 // solution(?) and basis that is, hopefully, optimal. This is
691 // confirmed or corrected by hot-starting the simplex solver
692 if (presolve_status == HighsPresolveStatus::ReducedToEmpty) {
693 clearSolutionUtil(hmos_[solved_hmo].solution_);
694 clearBasisUtil(hmos_[solved_hmo].basis_);
695 }
696
697 presolve_.data_.reduced_solution_ = hmos_[solved_hmo].solution_;
698 presolve_.data_.reduced_basis_.col_status =
699 hmos_[solved_hmo].basis_.col_status;
700 presolve_.data_.reduced_basis_.row_status =
701 hmos_[solved_hmo].basis_.row_status;
702
703 this_postsolve_time = -timer_.read(timer_.postsolve_clock);
704 timer_.start(timer_.postsolve_clock);
705 HighsPostsolveStatus postsolve_status = runPostsolve();
706 timer_.stop(timer_.postsolve_clock);
707 this_postsolve_time += -timer_.read(timer_.postsolve_clock);
708 presolve_.info_.postsolve_time = this_postsolve_time;
709
710 if (postsolve_status == HighsPostsolveStatus::SolutionRecovered) {
711 HighsPrintMessage(options_.output, options_.message_level, ML_VERBOSE,
712 "Postsolve finished\n");
713 //
714 // Now hot-start the simplex solver for the original_hmo:
715 //
716 // The original model hasn't been solved, so set up its solution
717 // parameters
718 resetModelStatusAndSolutionParams(hmos_[original_hmo]);
719 // Set solution and its status
720 hmos_[original_hmo].solution_ = presolve_.data_.recovered_solution_;
721
722 // Set basis and its status
723 hmos_[original_hmo].basis_.valid_ = true;
724 hmos_[original_hmo].basis_.col_status =
725 presolve_.data_.recovered_basis_.col_status;
726 hmos_[original_hmo].basis_.row_status =
727 presolve_.data_.recovered_basis_.row_status;
728
729 // Possibly force debug to perform KKT check on what's
730 // returned from postsolve
731 const bool force_debug = false;
732 int save_highs_debug_level = options_.highs_debug_level;
733 if (force_debug)
734 options_.highs_debug_level = HIGHS_DEBUG_LEVEL_COSTLY;
735 debugHighsBasicSolution("After returning from postsolve", options_,
736 lp_, hmos_[original_hmo].basis_,
737 hmos_[original_hmo].solution_);
738 options_.highs_debug_level = save_highs_debug_level;
739
740 // Now hot-start the simplex solver for the original_hmo
741 solved_hmo = original_hmo;
742 // Save the options to allow the best simplex strategy to
743 // be used
744 HighsOptions& options = hmos_[solved_hmo].options_;
745 HighsOptions save_options = options;
746 const bool full_logging = false;
747 if (full_logging) options.message_level = ML_ALWAYS;
748 // Force the use of simplex to clean up if IPM has been used
749 // to solve the presolved problem
750 if (options.solver == ipm_string) options.solver = simplex_string;
751 options.simplex_strategy = SIMPLEX_STRATEGY_CHOOSE;
752 // Ensure that the parallel solver isn't used
753 options.highs_min_threads = 1;
754 options.highs_max_threads = 1;
755 // Use any pivot threshold resulting from solving the presolved LP
756 if (factor_pivot_threshold > 0)
757 options.factor_pivot_threshold = factor_pivot_threshold;
758
759 hmos_[solved_hmo].lp_.lp_name_ = "Postsolve LP";
760 int iteration_count0 = info_.simplex_iteration_count;
761 this_solve_original_lp_time = -timer_.read(timer_.solve_clock);
762 timer_.start(timer_.solve_clock);
763 call_status = runLpSolver(
764 solved_hmo,
765 "Solving the original LP from the solution after postsolve");
766 timer_.stop(timer_.solve_clock);
767 this_solve_original_lp_time += timer_.read(timer_.solve_clock);
768 return_status =
769 interpretCallStatus(call_status, return_status, "runLpSolver");
770 // Recover the options
771 options = save_options;
772 if (return_status == HighsStatus::Error)
773 return returnFromRun(return_status);
774 postsolve_iteration_count =
775 info_.simplex_iteration_count - iteration_count0;
776 } else {
777 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
778 "Postsolve return status is %d\n",
779 (int)postsolve_status);
780 model_status_ = HighsModelStatus::POSTSOLVE_ERROR;
781 scaled_model_status_ = model_status_;
782 hmos_[0].unscaled_model_status_ = model_status_;
783 hmos_[0].scaled_model_status_ = model_status_;
784 return_status = HighsStatus::Error;
785 return returnFromRun(return_status);
786 }
787 }
788 } else {
789 // Optimal solution of presolved problem has not been found
790 // The original model inherits the solved model's status
791 hmos_[original_hmo].unscaled_model_status_ =
792 hmos_[solved_hmo].unscaled_model_status_;
793 hmos_[original_hmo].scaled_model_status_ =
794 hmos_[solved_hmo].scaled_model_status_;
795 }
796 } else {
797 // There is a valid basis for the problem or presolve is off
798 solved_hmo = original_hmo;
799 hmos_[solved_hmo].lp_.lp_name_ = "LP without presolve or with basis";
800 // There is a valid HiGHS basis, so use it to initialise the basis
801 // in the HMO to be solved
802 if (basis_.valid_) hmos_[solved_hmo].basis_ = basis_;
803 this_solve_original_lp_time = -timer_.read(timer_.solve_clock);
804 timer_.start(timer_.solve_clock);
805 call_status =
806 runLpSolver(solved_hmo, "Solving LP without presolve or with basis");
807 timer_.stop(timer_.solve_clock);
808 this_solve_original_lp_time += timer_.read(timer_.solve_clock);
809 return_status =
810 interpretCallStatus(call_status, return_status, "runLpSolver");
811 if (return_status == HighsStatus::Error)
812 return returnFromRun(return_status);
813 }
814 // else if (reduced problem failed to solve) {
815 // todo: handle case when presolved problem failed to solve. Try to solve
816 // again with no presolve.
817 // }
818
819 // assert(solved_hmo == original_hmo);
820 // solved_hmo will be original_hmo unless the presolved LP is found to be
821 // infeasible or unbounded
822
823 if (!getHighsModelStatusAndInfo(solved_hmo)) {
824 return_status = HighsStatus::Error;
825 return returnFromRun(return_status);
826 }
827
828 // Copy HMO solution/basis to HiGHS solution/basis: this resizes solution_ and
829 // basis_ The HiGHS solution and basis have to come from the original_hmo for
830 // them to have the right dimension.
831 solution_ = hmos_[original_hmo].solution_;
832 basis_ = hmos_[original_hmo].basis_;
833 // Stop and read the HiGHS clock, then work out time for this call
834 if (!run_highs_clock_already_running) timer_.stopRunHighsClock();
835
836 double lp_solve_final_time = timer_.readRunHighsClock();
837 double this_solve_time = lp_solve_final_time - initial_time;
838 if (postsolve_iteration_count < 0) {
839 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
840 "Postsolve : \n");
841 } else {
842 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
843 "Postsolve : %d\n", postsolve_iteration_count);
844 }
845 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
846 "Time : %8.2f\n", this_solve_time);
847 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
848 "Time Pre : %8.2f\n", this_presolve_time);
849 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
850 "Time PreLP : %8.2f\n", this_solve_presolved_lp_time);
851 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
852 "Time PostLP: %8.2f\n", this_solve_original_lp_time);
853 if (this_solve_time > 0) {
854 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
855 "For LP %16s",
856 hmos_[original_hmo].lp_.model_name_.c_str());
857 double sum_time = 0;
858 if (this_presolve_time > 0) {
859 sum_time += this_presolve_time;
860 int pct = (100 * this_presolve_time) / this_solve_time;
861 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
862 ": Presolve %8.2f (%3d%%)", this_presolve_time, pct);
863 }
864 if (this_solve_presolved_lp_time > 0) {
865 sum_time += this_solve_presolved_lp_time;
866 int pct = (100 * this_solve_presolved_lp_time) / this_solve_time;
867 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
868 ": Solve presolved LP %8.2f (%3d%%)",
869 this_solve_presolved_lp_time, pct);
870 }
871 if (this_postsolve_time > 0) {
872 sum_time += this_postsolve_time;
873 int pct = (100 * this_postsolve_time) / this_solve_time;
874 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
875 ": Postsolve %8.2f (%3d%%)", this_postsolve_time, pct);
876 }
877 if (this_solve_original_lp_time > 0) {
878 sum_time += this_solve_original_lp_time;
879 int pct = (100 * this_solve_original_lp_time) / this_solve_time;
880 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
881 ": Solve original LP %8.2f (%3d%%)",
882 this_solve_original_lp_time, pct);
883 }
884 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
885 "\n");
886 double rlv_time_difference =
887 fabs(sum_time - this_solve_time) / this_solve_time;
888 if (rlv_time_difference > 0.1)
889 HighsPrintMessage(options_.output, options_.message_level, ML_MINIMAL,
890 "Strange: Solve time = %g; Sum times = %g: relative "
891 "difference = %g\n",
892 this_solve_time, sum_time, rlv_time_difference);
893 }
894 // Assess success according to the scaled model status, unless
895 // something worse has happened earlier
896 call_status = highsStatusFromHighsModelStatus(scaled_model_status_);
897 return_status = interpretCallStatus(call_status, return_status);
898 return returnFromRun(return_status);
899 }
900
getLp() const901 const HighsLp& Highs::getLp() const { return lp_; }
902
getSolution() const903 const HighsSolution& Highs::getSolution() const { return solution_; }
904
getBasis() const905 const HighsBasis& Highs::getBasis() const { return basis_; }
906
getModelStatus(const bool scaled_model) const907 const HighsModelStatus& Highs::getModelStatus(const bool scaled_model) const {
908 if (scaled_model) {
909 return scaled_model_status_;
910 } else {
911 return model_status_;
912 }
913 }
914
getDualRay(bool & has_dual_ray,double * dual_ray_value)915 HighsStatus Highs::getDualRay(bool& has_dual_ray, double* dual_ray_value) {
916 if (!haveHmo("getDualRay")) return HighsStatus::Error;
917 HighsSimplexInterface simplex_interface(hmos_[0]);
918 return simplex_interface.getDualRay(has_dual_ray, dual_ray_value);
919 }
920
getPrimalRay(bool & has_primal_ray,double * primal_ray_value)921 HighsStatus Highs::getPrimalRay(bool& has_primal_ray,
922 double* primal_ray_value) {
923 underDevelopmentLogMessage("getPrimalRay");
924 if (!haveHmo("getPrimalRay")) return HighsStatus::Error;
925 HighsSimplexInterface simplex_interface(hmos_[0]);
926 return simplex_interface.getPrimalRay(has_primal_ray, primal_ray_value);
927 }
928
getRanging(HighsRanging & ranging)929 HighsStatus Highs::getRanging(HighsRanging& ranging) {
930 underDevelopmentLogMessage("getRanging");
931 if (!haveHmo("getRanging")) return HighsStatus::Error;
932 return getHighsRanging(ranging, hmos_[0]);
933 }
934
getBasicVariables(int * basic_variables)935 HighsStatus Highs::getBasicVariables(int* basic_variables) {
936 if (!haveHmo("getBasicVariables")) return HighsStatus::Error;
937 if (basic_variables == NULL) {
938 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
939 "getBasicVariables: basic_variables is NULL");
940 return HighsStatus::Error;
941 }
942 HighsSimplexInterface simplex_interface(hmos_[0]);
943 return simplex_interface.getBasicVariables(basic_variables);
944 }
945
getBasisInverseRow(const int row,double * row_vector,int * row_num_nz,int * row_indices)946 HighsStatus Highs::getBasisInverseRow(const int row, double* row_vector,
947 int* row_num_nz, int* row_indices) {
948 if (!haveHmo("getBasisInverseRow")) return HighsStatus::Error;
949 if (row_vector == NULL) {
950 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
951 "getBasisInverseRow: row_vector is NULL");
952 return HighsStatus::Error;
953 }
954 // row_indices can be NULL - it's the trigger that determines
955 // whether they are identified or not
956 int numRow = hmos_[0].lp_.numRow_;
957 if (row < 0 || row >= numRow) {
958 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
959 "Row index %d out of range [0, %d] in getBasisInverseRow",
960 row, numRow - 1);
961 return HighsStatus::Error;
962 }
963 if (!hmos_[0].simplex_lp_status_.has_invert) {
964 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
965 "No invertible representation for getBasisInverseRow");
966 return HighsStatus::Error;
967 }
968 // Compute a row i of the inverse of the basis matrix by solving B^Tx=e_i
969 vector<double> rhs;
970 rhs.assign(numRow, 0);
971 rhs[row] = 1;
972 HighsSimplexInterface simplex_interface(hmos_[0]);
973 simplex_interface.basisSolve(rhs, row_vector, row_num_nz, row_indices, true);
974 return HighsStatus::OK;
975 }
976
getBasisInverseCol(const int col,double * col_vector,int * col_num_nz,int * col_indices)977 HighsStatus Highs::getBasisInverseCol(const int col, double* col_vector,
978 int* col_num_nz, int* col_indices) {
979 if (!haveHmo("getBasisInverseCol")) return HighsStatus::Error;
980 if (col_vector == NULL) {
981 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
982 "getBasisInverseCol: col_vector is NULL");
983 return HighsStatus::Error;
984 }
985 // col_indices can be NULL - it's the trigger that determines
986 // whether they are identified or not
987 int numRow = hmos_[0].lp_.numRow_;
988 if (col < 0 || col >= numRow) {
989 HighsLogMessage(
990 options_.logfile, HighsMessageType::ERROR,
991 "Column index %d out of range [0, %d] in getBasisInverseCol", col,
992 numRow - 1);
993 return HighsStatus::Error;
994 }
995 if (!hmos_[0].simplex_lp_status_.has_invert) {
996 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
997 "No invertible representation for getBasisInverseCol");
998 return HighsStatus::Error;
999 }
1000 // Compute a col i of the inverse of the basis matrix by solving Bx=e_i
1001 vector<double> rhs;
1002 rhs.assign(numRow, 0);
1003 rhs[col] = 1;
1004 HighsSimplexInterface simplex_interface(hmos_[0]);
1005 simplex_interface.basisSolve(rhs, col_vector, col_num_nz, col_indices, false);
1006 return HighsStatus::OK;
1007 }
1008
getBasisSolve(const double * Xrhs,double * solution_vector,int * solution_num_nz,int * solution_indices)1009 HighsStatus Highs::getBasisSolve(const double* Xrhs, double* solution_vector,
1010 int* solution_num_nz, int* solution_indices) {
1011 if (!haveHmo("getBasisSolve")) return HighsStatus::Error;
1012 if (Xrhs == NULL) {
1013 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1014 "getBasisSolve: Xrhs is NULL");
1015 return HighsStatus::Error;
1016 }
1017 if (solution_vector == NULL) {
1018 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1019 "getBasisSolve: solution_vector is NULL");
1020 return HighsStatus::Error;
1021 }
1022 // solution_indices can be NULL - it's the trigger that determines
1023 // whether they are identified or not
1024 if (!hmos_[0].simplex_lp_status_.has_invert) {
1025 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1026 "No invertible representation for getBasisSolve");
1027 return HighsStatus::Error;
1028 }
1029 int numRow = hmos_[0].lp_.numRow_;
1030 vector<double> rhs;
1031 rhs.assign(numRow, 0);
1032 for (int row = 0; row < numRow; row++) rhs[row] = Xrhs[row];
1033 HighsSimplexInterface simplex_interface(hmos_[0]);
1034 simplex_interface.basisSolve(rhs, solution_vector, solution_num_nz,
1035 solution_indices, false);
1036 return HighsStatus::OK;
1037 }
1038
getBasisTransposeSolve(const double * Xrhs,double * solution_vector,int * solution_num_nz,int * solution_indices)1039 HighsStatus Highs::getBasisTransposeSolve(const double* Xrhs,
1040 double* solution_vector,
1041 int* solution_num_nz,
1042 int* solution_indices) {
1043 if (!haveHmo("getBasisTransposeSolve")) return HighsStatus::Error;
1044 if (Xrhs == NULL) {
1045 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1046 "getBasisTransposeSolve: Xrhs is NULL");
1047 return HighsStatus::Error;
1048 }
1049 if (solution_vector == NULL) {
1050 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1051 "getBasisTransposeSolve: solution_vector is NULL");
1052 return HighsStatus::Error;
1053 }
1054 // solution_indices can be NULL - it's the trigger that determines
1055 // whether they are identified or not
1056 if (!hmos_[0].simplex_lp_status_.has_invert) {
1057 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1058 "No invertible representation for getBasisTransposeSolve");
1059 return HighsStatus::Error;
1060 }
1061 int numRow = hmos_[0].lp_.numRow_;
1062 vector<double> rhs;
1063 rhs.assign(numRow, 0);
1064 for (int row = 0; row < numRow; row++) rhs[row] = Xrhs[row];
1065 HighsSimplexInterface simplex_interface(hmos_[0]);
1066 simplex_interface.basisSolve(rhs, solution_vector, solution_num_nz,
1067 solution_indices, true);
1068 return HighsStatus::OK;
1069 }
1070
getReducedRow(const int row,double * row_vector,int * row_num_nz,int * row_indices,const double * pass_basis_inverse_row_vector)1071 HighsStatus Highs::getReducedRow(const int row, double* row_vector,
1072 int* row_num_nz, int* row_indices,
1073 const double* pass_basis_inverse_row_vector) {
1074 if (!haveHmo("getReducedRow")) return HighsStatus::Error;
1075 if (row_vector == NULL) {
1076 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1077 "getReducedRow: row_vector is NULL");
1078 return HighsStatus::Error;
1079 }
1080 // row_indices can be NULL - it's the trigger that determines
1081 // whether they are identified or not pass_basis_inverse_row_vector
1082 // NULL - it's the trigger to determine whether it's computed or not
1083 if (row < 0 || row >= hmos_[0].lp_.numRow_) {
1084 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1085 "Row index %d out of range [0, %d] in getReducedRow", row,
1086 hmos_[0].lp_.numRow_ - 1);
1087 return HighsStatus::Error;
1088 }
1089 if (!hmos_[0].simplex_lp_status_.has_invert) {
1090 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1091 "No invertible representation for getReducedRow");
1092 return HighsStatus::Error;
1093 }
1094 HighsLp& lp = hmos_[0].lp_;
1095 int numRow = lp.numRow_;
1096 vector<double> basis_inverse_row;
1097 double* basis_inverse_row_vector = (double*)pass_basis_inverse_row_vector;
1098 if (basis_inverse_row_vector == NULL) {
1099 vector<double> rhs;
1100 vector<int> col_indices;
1101 rhs.assign(numRow, 0);
1102 rhs[row] = 1;
1103 basis_inverse_row.resize(numRow, 0);
1104 HighsSimplexInterface simplex_interface(hmos_[0]);
1105 // Form B^{-T}e_{row}
1106 simplex_interface.basisSolve(rhs, &basis_inverse_row[0], NULL, NULL, true);
1107 basis_inverse_row_vector = &basis_inverse_row[0];
1108 }
1109 bool return_indices = row_num_nz != NULL;
1110 if (return_indices) *row_num_nz = 0;
1111 for (int col = 0; col < lp.numCol_; col++) {
1112 double value = 0;
1113 for (int el = lp.Astart_[col]; el < lp.Astart_[col + 1]; el++) {
1114 int row = lp.Aindex_[el];
1115 value += lp.Avalue_[el] * basis_inverse_row_vector[row];
1116 }
1117 row_vector[col] = 0;
1118 if (fabs(value) > HIGHS_CONST_TINY) {
1119 if (return_indices) row_indices[(*row_num_nz)++] = col;
1120 row_vector[col] = value;
1121 }
1122 }
1123 return HighsStatus::OK;
1124 }
1125
getReducedColumn(const int col,double * col_vector,int * col_num_nz,int * col_indices)1126 HighsStatus Highs::getReducedColumn(const int col, double* col_vector,
1127 int* col_num_nz, int* col_indices) {
1128 if (!haveHmo("getReducedColumn")) return HighsStatus::Error;
1129 if (col_vector == NULL) {
1130 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1131 "getReducedColumn: col_vector is NULL");
1132 return HighsStatus::Error;
1133 }
1134 // col_indices can be NULL - it's the trigger that determines
1135 // whether they are identified or not
1136 if (col < 0 || col >= hmos_[0].lp_.numCol_) {
1137 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1138 "Column index %d out of range [0, %d] in getReducedColumn",
1139 col, hmos_[0].lp_.numCol_ - 1);
1140 return HighsStatus::Error;
1141 }
1142 if (!hmos_[0].simplex_lp_status_.has_invert) {
1143 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1144 "No invertible representation for getReducedColumn");
1145 return HighsStatus::Error;
1146 }
1147 HighsLp& lp = hmos_[0].lp_;
1148 int numRow = lp.numRow_;
1149 vector<double> rhs;
1150 rhs.assign(numRow, 0);
1151 for (int el = lp.Astart_[col]; el < lp.Astart_[col + 1]; el++)
1152 rhs[lp.Aindex_[el]] = lp.Avalue_[el];
1153 HighsSimplexInterface simplex_interface(hmos_[0]);
1154 simplex_interface.basisSolve(rhs, col_vector, col_num_nz, col_indices, false);
1155 return HighsStatus::OK;
1156 }
1157
setSolution(const HighsSolution & solution)1158 HighsStatus Highs::setSolution(const HighsSolution& solution) {
1159 HighsStatus return_status = HighsStatus::OK;
1160 // Check if solution is valid.
1161 assert((int)solution_.col_value.size() != 0 ||
1162 (int)solution_.col_value.size() != lp_.numCol_);
1163 assert((int)solution.col_dual.size() == 0 ||
1164 (int)solution.col_dual.size() == lp_.numCol_);
1165 assert((int)solution.row_dual.size() == 0 ||
1166 (int)solution.row_dual.size() == lp_.numRow_);
1167
1168 if (solution.col_value.size()) solution_.col_value = solution.col_value;
1169 if (solution.col_dual.size()) solution_.col_dual = solution.col_dual;
1170 if (solution.row_dual.size()) solution_.row_dual = solution.row_dual;
1171
1172 if (solution.col_value.size() > 0) {
1173 return_status = interpretCallStatus(calculateRowValues(lp_, solution_),
1174 return_status, "calculateRowValues");
1175 if (return_status == HighsStatus::Error) return return_status;
1176 }
1177 if (solution.row_dual.size() > 0) {
1178 return_status = interpretCallStatus(calculateColDuals(lp_, solution_),
1179 return_status, "calculateColDuals");
1180 if (return_status == HighsStatus::Error) return return_status;
1181 }
1182 return returnFromHighs(return_status);
1183 }
1184
setBasis(const HighsBasis & basis)1185 HighsStatus Highs::setBasis(const HighsBasis& basis) {
1186 // Check the user-supplied basis
1187 if (!isBasisConsistent(lp_, basis)) {
1188 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1189 "setBasis: invalid basis");
1190 return HighsStatus::Error;
1191 }
1192 // Update the HiGHS basis
1193 this->basis_ = basis;
1194 this->basis_.valid_ = true;
1195 // Follow implications of a new HiGHS basis
1196 newHighsBasis();
1197 // Can't use returnFromHighs since...
1198 return HighsStatus::OK;
1199 }
1200
setBasis()1201 HighsStatus Highs::setBasis() {
1202 // Invalidate the basis for HiGHS Don't set to logical basis since
1203 // that causes presolve to be skipped
1204 this->basis_.valid_ = false;
1205 // Follow implications of a new HiGHS basis
1206 newHighsBasis();
1207 // Can't use returnFromHighs since...
1208 return HighsStatus::OK;
1209 }
1210
addRow(const double lower_bound,const double upper_bound,const int num_new_nz,const int * indices,const double * values)1211 bool Highs::addRow(const double lower_bound, const double upper_bound,
1212 const int num_new_nz, const int* indices,
1213 const double* values) {
1214 int starts = 0;
1215 return addRows(1, &lower_bound, &upper_bound, num_new_nz, &starts, indices,
1216 values);
1217 }
1218
addRows(const int num_new_row,const double * lower_bounds,const double * upper_bounds,const int num_new_nz,const int * starts,const int * indices,const double * values)1219 bool Highs::addRows(const int num_new_row, const double* lower_bounds,
1220 const double* upper_bounds, const int num_new_nz,
1221 const int* starts, const int* indices,
1222 const double* values) {
1223 HighsStatus return_status = HighsStatus::OK;
1224 // Check that there is a HighsModelObject
1225 if (!haveHmo("addRows")) return false;
1226 HighsSimplexInterface interface(hmos_[0]);
1227 return_status = interpretCallStatus(
1228 interface.addRows(num_new_row, lower_bounds, upper_bounds, num_new_nz,
1229 starts, indices, values),
1230 return_status, "addRows");
1231 if (return_status == HighsStatus::Error) return false;
1232 return returnFromHighs(return_status) != HighsStatus::Error;
1233 }
1234
addCol(const double cost,const double lower_bound,const double upper_bound,const int num_new_nz,const int * indices,const double * values)1235 bool Highs::addCol(const double cost, const double lower_bound,
1236 const double upper_bound, const int num_new_nz,
1237 const int* indices, const double* values) {
1238 int starts = 0;
1239 return addCols(1, &cost, &lower_bound, &upper_bound, num_new_nz, &starts,
1240 indices, values);
1241 }
1242
addCols(const int num_new_col,const double * costs,const double * lower_bounds,const double * upper_bounds,const int num_new_nz,const int * starts,const int * indices,const double * values)1243 bool Highs::addCols(const int num_new_col, const double* costs,
1244 const double* lower_bounds, const double* upper_bounds,
1245 const int num_new_nz, const int* starts, const int* indices,
1246 const double* values) {
1247 HighsStatus return_status = HighsStatus::OK;
1248 if (!haveHmo("addCols")) return false;
1249 HighsSimplexInterface interface(hmos_[0]);
1250 return_status = interpretCallStatus(
1251 interface.addCols(num_new_col, costs, lower_bounds, upper_bounds,
1252 num_new_nz, starts, indices, values),
1253 return_status, "addCols");
1254 if (return_status == HighsStatus::Error) return false;
1255 return returnFromHighs(return_status) != HighsStatus::Error;
1256 }
1257
changeObjectiveSense(const ObjSense sense)1258 bool Highs::changeObjectiveSense(const ObjSense sense) {
1259 HighsStatus return_status = HighsStatus::OK;
1260 if (!haveHmo("changeObjectiveSense")) return false;
1261 HighsSimplexInterface interface(hmos_[0]);
1262 return_status = interpretCallStatus(interface.changeObjectiveSense(sense),
1263 return_status, "changeObjectiveSense");
1264 if (return_status == HighsStatus::Error) return false;
1265 return returnFromHighs(return_status) != HighsStatus::Error;
1266 }
1267
changeColCost(const int col,const double cost)1268 bool Highs::changeColCost(const int col, const double cost) {
1269 return changeColsCost(1, &col, &cost);
1270 }
1271
changeColsCost(const int num_set_entries,const int * set,const double * cost)1272 bool Highs::changeColsCost(const int num_set_entries, const int* set,
1273 const double* cost) {
1274 if (num_set_entries <= 0) return true;
1275 HighsStatus return_status = HighsStatus::OK;
1276 HighsStatus call_status;
1277 // Create a local set that is not const since index_collection.set_
1278 // cannot be const as it may change if the set is not ordered
1279 vector<int> local_set{set, set + num_set_entries};
1280 HighsIndexCollection index_collection;
1281 index_collection.dimension_ = lp_.numCol_;
1282 index_collection.is_set_ = true;
1283 index_collection.set_ = &local_set[0];
1284 index_collection.set_num_entries_ = num_set_entries;
1285 if (!haveHmo("changeColsCost")) return false;
1286 HighsSimplexInterface interface(hmos_[0]);
1287 call_status = interface.changeCosts(index_collection, cost);
1288 return_status =
1289 interpretCallStatus(call_status, return_status, "changeCosts");
1290 if (return_status == HighsStatus::Error) return false;
1291 return returnFromHighs(return_status) != HighsStatus::Error;
1292 }
1293
changeColsCost(const int * mask,const double * cost)1294 bool Highs::changeColsCost(const int* mask, const double* cost) {
1295 HighsStatus return_status = HighsStatus::OK;
1296 HighsStatus call_status;
1297 // Create a local mask that is not const since
1298 // index_collection.mask_ cannot be const as it changes when
1299 // deleting rows/columns
1300 vector<int> local_mask{mask, mask + lp_.numCol_};
1301 HighsIndexCollection index_collection;
1302 index_collection.dimension_ = lp_.numCol_;
1303 index_collection.is_mask_ = true;
1304 index_collection.mask_ = &local_mask[0];
1305 if (!haveHmo("changeColsCost")) return false;
1306 HighsSimplexInterface interface(hmos_[0]);
1307 call_status = interface.changeCosts(index_collection, cost);
1308 return_status =
1309 interpretCallStatus(call_status, return_status, "changeCosts");
1310 if (return_status == HighsStatus::Error) return false;
1311 return returnFromHighs(return_status) != HighsStatus::Error;
1312 }
1313
changeColBounds(const int col,const double lower,const double upper)1314 bool Highs::changeColBounds(const int col, const double lower,
1315 const double upper) {
1316 return changeColsBounds(1, &col, &lower, &upper);
1317 }
1318
changeColsBounds(const int from_col,const int to_col,const double * lower,const double * upper)1319 bool Highs::changeColsBounds(const int from_col, const int to_col,
1320 const double* lower, const double* upper) {
1321 HighsStatus return_status = HighsStatus::OK;
1322 HighsStatus call_status;
1323 HighsIndexCollection index_collection;
1324 index_collection.dimension_ = lp_.numCol_;
1325 index_collection.is_interval_ = true;
1326 index_collection.from_ = from_col;
1327 index_collection.to_ = to_col;
1328 if (!haveHmo("changeColsBounds")) return false;
1329 HighsSimplexInterface interface(hmos_[0]);
1330 call_status = interface.changeColBounds(index_collection, lower, upper);
1331 return_status =
1332 interpretCallStatus(call_status, return_status, "changeColBounds");
1333 if (return_status == HighsStatus::Error) return false;
1334 return returnFromHighs(return_status) != HighsStatus::Error;
1335 }
1336
changeColsBounds(const int num_set_entries,const int * set,const double * lower,const double * upper)1337 bool Highs::changeColsBounds(const int num_set_entries, const int* set,
1338 const double* lower, const double* upper) {
1339 if (num_set_entries <= 0) return true;
1340 HighsStatus return_status = HighsStatus::OK;
1341 HighsStatus call_status;
1342 // Create a local set that is not const since index_collection.set_
1343 // cannot be const as it may change if the set is not ordered
1344 vector<int> local_set{set, set + num_set_entries};
1345 HighsIndexCollection index_collection;
1346 index_collection.dimension_ = lp_.numCol_;
1347 index_collection.is_set_ = true;
1348 index_collection.set_ = &local_set[0];
1349 index_collection.set_num_entries_ = num_set_entries;
1350 if (!haveHmo("changeColsBounds")) return false;
1351 HighsSimplexInterface interface(hmos_[0]);
1352 call_status = interface.changeColBounds(index_collection, lower, upper);
1353 return_status =
1354 interpretCallStatus(call_status, return_status, "changeColBounds");
1355 if (return_status == HighsStatus::Error) return false;
1356 return returnFromHighs(return_status) != HighsStatus::Error;
1357 }
1358
changeColsBounds(const int * mask,const double * lower,const double * upper)1359 bool Highs::changeColsBounds(const int* mask, const double* lower,
1360 const double* upper) {
1361 HighsStatus return_status = HighsStatus::OK;
1362 HighsStatus call_status;
1363 // Create a local mask that is not const since
1364 // index_collection.mask_ cannot be const as it changes when
1365 // deleting rows/columns
1366 vector<int> local_mask{mask, mask + lp_.numCol_};
1367 HighsIndexCollection index_collection;
1368 index_collection.dimension_ = lp_.numCol_;
1369 index_collection.is_mask_ = true;
1370 index_collection.mask_ = &local_mask[0];
1371 if (!haveHmo("changeColsBounds")) return false;
1372 HighsSimplexInterface interface(hmos_[0]);
1373 call_status = interface.changeColBounds(index_collection, lower, upper);
1374 return_status =
1375 interpretCallStatus(call_status, return_status, "changeColBounds");
1376 if (return_status == HighsStatus::Error) return false;
1377 return returnFromHighs(return_status) != HighsStatus::Error;
1378 }
1379
changeRowBounds(const int row,const double lower,const double upper)1380 bool Highs::changeRowBounds(const int row, const double lower,
1381 const double upper) {
1382 return changeRowsBounds(1, &row, &lower, &upper);
1383 }
1384
changeRowsBounds(const int num_set_entries,const int * set,const double * lower,const double * upper)1385 bool Highs::changeRowsBounds(const int num_set_entries, const int* set,
1386 const double* lower, const double* upper) {
1387 if (num_set_entries <= 0) return true;
1388 HighsStatus return_status = HighsStatus::OK;
1389 HighsStatus call_status;
1390 // Create a local set that is not const since index_collection.set_
1391 // cannot be const as it may change if the set is not ordered
1392 vector<int> local_set{set, set + num_set_entries};
1393 HighsIndexCollection index_collection;
1394 index_collection.dimension_ = lp_.numRow_;
1395 index_collection.is_set_ = true;
1396 index_collection.set_ = &local_set[0];
1397 index_collection.set_num_entries_ = num_set_entries;
1398 if (!haveHmo("changeRowsBounds")) return false;
1399 HighsSimplexInterface interface(hmos_[0]);
1400 call_status = interface.changeRowBounds(index_collection, lower, upper);
1401 return_status =
1402 interpretCallStatus(call_status, return_status, "changeRowBounds");
1403 if (return_status == HighsStatus::Error) return false;
1404 return returnFromHighs(return_status) != HighsStatus::Error;
1405 }
1406
changeRowsBounds(const int * mask,const double * lower,const double * upper)1407 bool Highs::changeRowsBounds(const int* mask, const double* lower,
1408 const double* upper) {
1409 HighsStatus return_status = HighsStatus::OK;
1410 HighsStatus call_status;
1411 // Create a local mask that is not const since
1412 // index_collection.mask_ cannot be const as it changes when
1413 // deleting rows/columns
1414 vector<int> local_mask{mask, mask + lp_.numRow_};
1415 HighsIndexCollection index_collection;
1416 index_collection.dimension_ = lp_.numRow_;
1417 index_collection.is_mask_ = true;
1418 index_collection.mask_ = &local_mask[0];
1419 if (!haveHmo("changeRowsBounds")) return false;
1420 HighsSimplexInterface interface(hmos_[0]);
1421 call_status = interface.changeRowBounds(index_collection, lower, upper);
1422 return_status =
1423 interpretCallStatus(call_status, return_status, "changeRowBounds");
1424 if (return_status == HighsStatus::Error) return false;
1425 return returnFromHighs(return_status) != HighsStatus::Error;
1426 }
1427
changeCoeff(const int row,const int col,const double value)1428 bool Highs::changeCoeff(const int row, const int col, const double value) {
1429 HighsStatus return_status = HighsStatus::OK;
1430 HighsStatus call_status;
1431 if (!haveHmo("changeCoeff")) return false;
1432 HighsSimplexInterface interface(hmos_[0]);
1433 call_status = interface.changeCoefficient(row, col, value);
1434 return_status =
1435 interpretCallStatus(call_status, return_status, "changeCoefficient");
1436 if (return_status == HighsStatus::Error) return false;
1437 return returnFromHighs(return_status) != HighsStatus::Error;
1438 }
1439
getObjectiveSense(ObjSense & sense)1440 bool Highs::getObjectiveSense(ObjSense& sense) {
1441 if (!haveHmo("getObjectiveSense")) return false;
1442 sense = hmos_[0].lp_.sense_;
1443 return true;
1444 }
1445
getCols(const int from_col,const int to_col,int & num_col,double * costs,double * lower,double * upper,int & num_nz,int * start,int * index,double * value)1446 bool Highs::getCols(const int from_col, const int to_col, int& num_col,
1447 double* costs, double* lower, double* upper, int& num_nz,
1448 int* start, int* index, double* value) {
1449 HighsStatus return_status = HighsStatus::OK;
1450 HighsStatus call_status;
1451 HighsIndexCollection index_collection;
1452 index_collection.dimension_ = lp_.numCol_;
1453 index_collection.is_interval_ = true;
1454 index_collection.from_ = from_col;
1455 index_collection.to_ = to_col;
1456 if (!haveHmo("getCols")) return false;
1457 HighsSimplexInterface interface(hmos_[0]);
1458 call_status = interface.getCols(index_collection, num_col, costs, lower,
1459 upper, num_nz, start, index, value);
1460 return_status = interpretCallStatus(call_status, return_status, "getCols");
1461 if (return_status == HighsStatus::Error) return false;
1462 return returnFromHighs(return_status) != HighsStatus::Error;
1463 }
1464
getCols(const int num_set_entries,const int * set,int & num_col,double * costs,double * lower,double * upper,int & num_nz,int * start,int * index,double * value)1465 bool Highs::getCols(const int num_set_entries, const int* set, int& num_col,
1466 double* costs, double* lower, double* upper, int& num_nz,
1467 int* start, int* index, double* value) {
1468 if (num_set_entries <= 0) return true;
1469 HighsStatus return_status = HighsStatus::OK;
1470 HighsStatus call_status;
1471 // Create a local set that is not const since index_collection.set_
1472 // cannot be const as it may change if the set is not ordered
1473 vector<int> local_set{set, set + num_set_entries};
1474 HighsIndexCollection index_collection;
1475 index_collection.dimension_ = lp_.numCol_;
1476 index_collection.is_set_ = true;
1477 index_collection.set_ = &local_set[0];
1478 index_collection.set_num_entries_ = num_set_entries;
1479 if (!haveHmo("getCols")) return false;
1480 HighsSimplexInterface interface(hmos_[0]);
1481 call_status = interface.getCols(index_collection, num_col, costs, lower,
1482 upper, num_nz, start, index, value);
1483 return_status = interpretCallStatus(call_status, return_status, "getCols");
1484 if (return_status == HighsStatus::Error) return false;
1485 return returnFromHighs(return_status) != HighsStatus::Error;
1486 }
1487
getCols(const int * mask,int & num_col,double * costs,double * lower,double * upper,int & num_nz,int * start,int * index,double * value)1488 bool Highs::getCols(const int* mask, int& num_col, double* costs, double* lower,
1489 double* upper, int& num_nz, int* start, int* index,
1490 double* value) {
1491 HighsStatus return_status = HighsStatus::OK;
1492 HighsStatus call_status;
1493 // Create a local mask that is not const since
1494 // index_collection.mask_ cannot be const as it changes when
1495 // deleting rows/columns
1496 vector<int> local_mask{mask, mask + lp_.numCol_};
1497 HighsIndexCollection index_collection;
1498 index_collection.dimension_ = lp_.numCol_;
1499 index_collection.is_mask_ = true;
1500 index_collection.mask_ = &local_mask[0];
1501 if (!haveHmo("getCols")) return false;
1502 HighsSimplexInterface interface(hmos_[0]);
1503 call_status = interface.getCols(index_collection, num_col, costs, lower,
1504 upper, num_nz, start, index, value);
1505 return_status = interpretCallStatus(call_status, return_status, "getCols");
1506 if (return_status == HighsStatus::Error) return false;
1507 return returnFromHighs(return_status) != HighsStatus::Error;
1508 }
1509
getRows(const int from_row,const int to_row,int & num_row,double * lower,double * upper,int & num_nz,int * start,int * index,double * value)1510 bool Highs::getRows(const int from_row, const int to_row, int& num_row,
1511 double* lower, double* upper, int& num_nz, int* start,
1512 int* index, double* value) {
1513 HighsStatus return_status = HighsStatus::OK;
1514 HighsStatus call_status;
1515 HighsIndexCollection index_collection;
1516 index_collection.dimension_ = lp_.numRow_;
1517 index_collection.is_interval_ = true;
1518 index_collection.from_ = from_row;
1519 index_collection.to_ = to_row;
1520 if (!haveHmo("getRows")) return false;
1521 HighsSimplexInterface interface(hmos_[0]);
1522 call_status = interface.getRows(index_collection, num_row, lower, upper,
1523 num_nz, start, index, value);
1524 return_status = interpretCallStatus(call_status, return_status, "getRows");
1525 if (return_status == HighsStatus::Error) return false;
1526 return returnFromHighs(return_status) != HighsStatus::Error;
1527 }
1528
getRows(const int num_set_entries,const int * set,int & num_row,double * lower,double * upper,int & num_nz,int * start,int * index,double * value)1529 bool Highs::getRows(const int num_set_entries, const int* set, int& num_row,
1530 double* lower, double* upper, int& num_nz, int* start,
1531 int* index, double* value) {
1532 if (num_set_entries <= 0) return true;
1533 HighsStatus return_status = HighsStatus::OK;
1534 HighsStatus call_status;
1535 // Create a local set that is not const since index_collection.set_
1536 // cannot be const as it may change if the set is not ordered
1537 vector<int> local_set{set, set + num_set_entries};
1538 HighsIndexCollection index_collection;
1539 index_collection.dimension_ = lp_.numRow_;
1540 index_collection.is_set_ = true;
1541 index_collection.set_ = &local_set[0];
1542 index_collection.set_num_entries_ = num_set_entries;
1543 if (!haveHmo("getRows")) return false;
1544 HighsSimplexInterface interface(hmos_[0]);
1545 call_status = interface.getRows(index_collection, num_row, lower, upper,
1546 num_nz, start, index, value);
1547 return_status = interpretCallStatus(call_status, return_status, "getRows");
1548 if (return_status == HighsStatus::Error) return false;
1549 return returnFromHighs(return_status) != HighsStatus::Error;
1550 }
1551
getRows(const int * mask,int & num_row,double * lower,double * upper,int & num_nz,int * start,int * index,double * value)1552 bool Highs::getRows(const int* mask, int& num_row, double* lower, double* upper,
1553 int& num_nz, int* start, int* index, double* value) {
1554 HighsStatus return_status = HighsStatus::OK;
1555 HighsStatus call_status;
1556 // Create a local mask that is not const since
1557 // index_collection.mask_ cannot be const as it changes when
1558 // deleting rows/columns
1559 vector<int> local_mask{mask, mask + lp_.numRow_};
1560 HighsIndexCollection index_collection;
1561 index_collection.dimension_ = lp_.numRow_;
1562 index_collection.is_mask_ = true;
1563 index_collection.mask_ = &local_mask[0];
1564 if (!haveHmo("getRows")) return false;
1565 HighsSimplexInterface interface(hmos_[0]);
1566 call_status = interface.getRows(index_collection, num_row, lower, upper,
1567 num_nz, start, index, value);
1568 return_status = interpretCallStatus(call_status, return_status, "getRows");
1569 if (return_status == HighsStatus::Error) return false;
1570 return returnFromHighs(return_status) != HighsStatus::Error;
1571 }
1572
getCoeff(const int row,const int col,double & value)1573 bool Highs::getCoeff(const int row, const int col, double& value) {
1574 HighsStatus return_status = HighsStatus::OK;
1575 HighsStatus call_status;
1576 if (!haveHmo("getCoeff")) return false;
1577 HighsSimplexInterface interface(hmos_[0]);
1578
1579 call_status = interface.getCoefficient(row, col, value);
1580 return_status =
1581 interpretCallStatus(call_status, return_status, "getCoefficient");
1582 if (return_status == HighsStatus::Error) return false;
1583 return returnFromHighs(return_status) != HighsStatus::Error;
1584 }
1585
deleteCols(const int from_col,const int to_col)1586 bool Highs::deleteCols(const int from_col, const int to_col) {
1587 HighsStatus return_status = HighsStatus::OK;
1588 HighsStatus call_status;
1589 HighsIndexCollection index_collection;
1590 index_collection.dimension_ = lp_.numCol_;
1591 index_collection.is_interval_ = true;
1592 index_collection.from_ = from_col;
1593 index_collection.to_ = to_col;
1594 if (!haveHmo("deleteCols")) return false;
1595 HighsSimplexInterface interface(hmos_[0]);
1596 call_status = interface.deleteCols(index_collection);
1597 return_status = interpretCallStatus(call_status, return_status, "deleteCols");
1598 if (return_status == HighsStatus::Error) return false;
1599 return returnFromHighs(return_status) != HighsStatus::Error;
1600 }
1601
deleteCols(const int num_set_entries,const int * set)1602 bool Highs::deleteCols(const int num_set_entries, const int* set) {
1603 if (num_set_entries <= 0) return true;
1604 HighsStatus return_status = HighsStatus::OK;
1605 HighsStatus call_status;
1606 // Create a local set that is not const since index_collection.set_
1607 // cannot be const as it may change if the set is not ordered
1608 vector<int> local_set{set, set + num_set_entries};
1609 HighsIndexCollection index_collection;
1610 index_collection.dimension_ = lp_.numCol_;
1611 index_collection.is_set_ = true;
1612 index_collection.set_ = &local_set[0];
1613 index_collection.set_num_entries_ = num_set_entries;
1614 if (!haveHmo("deleteCols")) return false;
1615 HighsSimplexInterface interface(hmos_[0]);
1616 call_status = interface.deleteCols(index_collection);
1617 return_status = interpretCallStatus(call_status, return_status, "deleteCols");
1618 if (return_status == HighsStatus::Error) return false;
1619 return returnFromHighs(return_status) != HighsStatus::Error;
1620 }
1621
deleteCols(int * mask)1622 bool Highs::deleteCols(int* mask) {
1623 HighsStatus return_status = HighsStatus::OK;
1624 HighsStatus call_status;
1625 HighsIndexCollection index_collection;
1626 index_collection.dimension_ = lp_.numCol_;
1627 index_collection.is_mask_ = true;
1628 index_collection.mask_ = &mask[0];
1629 if (!haveHmo("deleteCols")) return false;
1630 HighsSimplexInterface interface(hmos_[0]);
1631 call_status = interface.deleteCols(index_collection);
1632 return_status = interpretCallStatus(call_status, return_status, "deleteCols");
1633 if (return_status == HighsStatus::Error) return false;
1634 return returnFromHighs(return_status) != HighsStatus::Error;
1635 }
1636
deleteRows(const int from_row,const int to_row)1637 bool Highs::deleteRows(const int from_row, const int to_row) {
1638 HighsStatus return_status = HighsStatus::OK;
1639 HighsStatus call_status;
1640 HighsIndexCollection index_collection;
1641 index_collection.dimension_ = lp_.numRow_;
1642 index_collection.is_interval_ = true;
1643 index_collection.from_ = from_row;
1644 index_collection.to_ = to_row;
1645 if (!haveHmo("deleteRows")) return false;
1646 HighsSimplexInterface interface(hmos_[0]);
1647 call_status = interface.deleteRows(index_collection);
1648 return_status = interpretCallStatus(call_status, return_status, "deleteRows");
1649 if (return_status == HighsStatus::Error) return false;
1650 return returnFromHighs(return_status) != HighsStatus::Error;
1651 }
1652
deleteRows(const int num_set_entries,const int * set)1653 bool Highs::deleteRows(const int num_set_entries, const int* set) {
1654 if (num_set_entries <= 0) return true;
1655 HighsStatus return_status = HighsStatus::OK;
1656 HighsStatus call_status;
1657 // Create a local set that is not const since index_collection.set_
1658 // cannot be const as it may change if the set is not ordered
1659 vector<int> local_set{set, set + num_set_entries};
1660 HighsIndexCollection index_collection;
1661 index_collection.dimension_ = lp_.numRow_;
1662 index_collection.is_set_ = true;
1663 index_collection.set_ = &local_set[0];
1664 index_collection.set_num_entries_ = num_set_entries;
1665 if (!haveHmo("deleteRows")) return false;
1666 HighsSimplexInterface interface(hmos_[0]);
1667 call_status = interface.deleteRows(index_collection);
1668 return_status = interpretCallStatus(call_status, return_status, "deleteRows");
1669 if (return_status == HighsStatus::Error) return false;
1670 return returnFromHighs(return_status) != HighsStatus::Error;
1671 }
1672
deleteRows(int * mask)1673 bool Highs::deleteRows(int* mask) {
1674 HighsStatus return_status = HighsStatus::OK;
1675 HighsStatus call_status;
1676 HighsIndexCollection index_collection;
1677 index_collection.dimension_ = lp_.numRow_;
1678 index_collection.is_mask_ = true;
1679 index_collection.mask_ = &mask[0];
1680 if (!haveHmo("deleteRows")) return false;
1681 HighsSimplexInterface interface(hmos_[0]);
1682 call_status = interface.deleteRows(index_collection);
1683 return_status = interpretCallStatus(call_status, return_status, "deleteRows");
1684 if (return_status == HighsStatus::Error) return false;
1685 return returnFromHighs(return_status) != HighsStatus::Error;
1686 }
1687
scaleCol(const int col,const double scaleval)1688 bool Highs::scaleCol(const int col, const double scaleval) {
1689 HighsStatus return_status = HighsStatus::OK;
1690 HighsStatus call_status;
1691 if (!haveHmo("scaleCol")) return false;
1692 HighsSimplexInterface interface(hmos_[0]);
1693 call_status = interface.scaleCol(col, scaleval);
1694 return_status = interpretCallStatus(call_status, return_status, "scaleCol");
1695 if (return_status == HighsStatus::Error) return false;
1696 return returnFromHighs(return_status) != HighsStatus::Error;
1697 }
1698
scaleRow(const int row,const double scaleval)1699 bool Highs::scaleRow(const int row, const double scaleval) {
1700 HighsStatus return_status = HighsStatus::OK;
1701 HighsStatus call_status;
1702 if (!haveHmo("scaleRow")) return false;
1703 HighsSimplexInterface interface(hmos_[0]);
1704 call_status = interface.scaleRow(row, scaleval);
1705 return_status = interpretCallStatus(call_status, return_status, "scaleRow");
1706 if (return_status == HighsStatus::Error) return false;
1707 return returnFromHighs(return_status) != HighsStatus::Error;
1708 }
1709
getHighsInfinity()1710 double Highs::getHighsInfinity() { return HIGHS_CONST_INF; }
1711
getHighsRunTime()1712 double Highs::getHighsRunTime() { return timer_.readRunHighsClock(); }
1713
clearSolver()1714 HighsStatus Highs::clearSolver() {
1715 clearModelStatus();
1716 clearSolution();
1717 clearBasis();
1718 clearInfo();
1719 return HighsStatus::OK;
1720 }
1721
1722 #ifdef HiGHSDEV
reportModelStatusSolutionBasis(const std::string message,const int hmo_ix)1723 void Highs::reportModelStatusSolutionBasis(const std::string message,
1724 const int hmo_ix) {
1725 HighsModelStatus& model_status = model_status_;
1726 HighsModelStatus& scaled_model_status = scaled_model_status_;
1727 HighsSolution& solution = solution_;
1728 HighsBasis& basis = basis_;
1729 int unscaled_primal_status = info_.primal_status;
1730 int scaled_primal_status = unscaled_primal_status;
1731 int unscaled_dual_status = info_.dual_status;
1732 int scaled_dual_status = unscaled_dual_status;
1733 HighsLp& lp = lp_;
1734 if (hmo_ix >= 0) {
1735 assert(hmo_ix < (int)hmos_.size());
1736 model_status = hmos_[hmo_ix].unscaled_model_status_;
1737 scaled_model_status = hmos_[hmo_ix].scaled_model_status_;
1738 solution = hmos_[hmo_ix].solution_;
1739 basis = hmos_[hmo_ix].basis_;
1740 unscaled_primal_status =
1741 hmos_[hmo_ix].unscaled_solution_params_.primal_status;
1742 scaled_primal_status = hmos_[hmo_ix].scaled_solution_params_.primal_status;
1743 unscaled_dual_status = hmos_[hmo_ix].unscaled_solution_params_.dual_status;
1744 scaled_dual_status = hmos_[hmo_ix].scaled_solution_params_.dual_status;
1745 lp = hmos_[hmo_ix].lp_;
1746 }
1747 printf(
1748 "\n%s\nModel status = %s; Scaled model status = %s; LP(%d, %d); solution "
1749 "([%d:%d] %d, %d; [%d:%d] %d, %d); basis %d "
1750 "(%d, %d)\n\n",
1751 message.c_str(), utilHighsModelStatusToString(model_status).c_str(),
1752 utilHighsModelStatusToString(scaled_model_status).c_str(), lp.numCol_,
1753 lp.numRow_, unscaled_primal_status, scaled_primal_status,
1754 (int)solution.col_value.size(), (int)solution.row_value.size(),
1755 unscaled_dual_status, scaled_dual_status, (int)solution.col_dual.size(),
1756 (int)solution.row_dual.size(), basis.valid_, (int)basis.col_status.size(),
1757 (int)basis.row_status.size());
1758 }
1759 #endif
1760
highsModelStatusToString(const HighsModelStatus model_status) const1761 std::string Highs::highsModelStatusToString(
1762 const HighsModelStatus model_status) const {
1763 return utilHighsModelStatusToString(model_status);
1764 }
1765
primalDualStatusToString(const int primal_dual_status)1766 std::string Highs::primalDualStatusToString(const int primal_dual_status) {
1767 return utilPrimalDualStatusToString(primal_dual_status);
1768 }
1769
1770 // Private methods
runPresolve()1771 HighsPresolveStatus Highs::runPresolve() {
1772 // Exit if the problem is empty or if presolve is set to off.
1773 if (options_.presolve == off_string) return HighsPresolveStatus::NotPresolved;
1774 if (lp_.numCol_ == 0 && lp_.numRow_ == 0)
1775 return HighsPresolveStatus::NullError;
1776
1777 // Clear info from previous runs if lp_ has been modified.
1778 if (presolve_.has_run_) presolve_.clear();
1779 double start_presolve = timer_.readRunHighsClock();
1780
1781 // Set time limit.
1782 if (options_.time_limit > 0 && options_.time_limit < HIGHS_CONST_INF) {
1783 double left = options_.time_limit - start_presolve;
1784 if (left <= 0) {
1785 HighsPrintMessage(options_.output, options_.message_level, ML_VERBOSE,
1786 "Time limit reached while reading in matrix\n");
1787 return HighsPresolveStatus::Timeout;
1788 }
1789
1790 HighsPrintMessage(options_.output, options_.message_level, ML_VERBOSE,
1791 "Time limit set: reading matrix took %.2g, presolve "
1792 "time left: %.2g\n",
1793 start_presolve, left);
1794 presolve_.options_.time_limit = left;
1795 }
1796
1797 // Presolve.
1798 presolve_.init(lp_, timer_);
1799 if (options_.time_limit > 0 && options_.time_limit < HIGHS_CONST_INF) {
1800 double current = timer_.readRunHighsClock();
1801 double time_init = current - start_presolve;
1802 double left = presolve_.options_.time_limit - time_init;
1803 if (left <= 0) {
1804 HighsPrintMessage(
1805 options_.output, options_.message_level, ML_VERBOSE,
1806 "Time limit reached while copying matrix into presolve.\n");
1807 return HighsPresolveStatus::Timeout;
1808 }
1809
1810 HighsPrintMessage(options_.output, options_.message_level, ML_VERBOSE,
1811 "Time limit set: copying matrix took %.2g, presolve "
1812 "time left: %.2g\n",
1813 time_init, left);
1814 presolve_.options_.time_limit = options_.time_limit;
1815 }
1816
1817 presolve_.data_.presolve_[0].message_level = options_.message_level;
1818 presolve_.data_.presolve_[0].output = options_.output;
1819
1820 HighsPresolveStatus presolve_return_status = presolve_.run();
1821
1822 // Handle max case.
1823 if (presolve_return_status == HighsPresolveStatus::Reduced &&
1824 lp_.sense_ == ObjSense::MAXIMIZE) {
1825 presolve_.negateReducedLpCost();
1826 presolve_.data_.reduced_lp_.sense_ = ObjSense::MAXIMIZE;
1827 }
1828
1829 // Update reduction counts.
1830 switch (presolve_.presolve_status_) {
1831 case HighsPresolveStatus::Reduced: {
1832 HighsLp& reduced_lp = presolve_.getReducedProblem();
1833 presolve_.info_.n_cols_removed = lp_.numCol_ - reduced_lp.numCol_;
1834 presolve_.info_.n_rows_removed = lp_.numRow_ - reduced_lp.numRow_;
1835 presolve_.info_.n_nnz_removed =
1836 (int)lp_.Avalue_.size() - (int)reduced_lp.Avalue_.size();
1837 break;
1838 }
1839 case HighsPresolveStatus::ReducedToEmpty: {
1840 presolve_.info_.n_cols_removed = lp_.numCol_;
1841 presolve_.info_.n_rows_removed = lp_.numRow_;
1842 presolve_.info_.n_nnz_removed = (int)lp_.Avalue_.size();
1843 break;
1844 }
1845 default:
1846 break;
1847 }
1848 return presolve_return_status;
1849 }
1850
runPostsolve()1851 HighsPostsolveStatus Highs::runPostsolve() {
1852 assert(presolve_.has_run_);
1853 bool solution_ok = isSolutionRightSize(presolve_.getReducedProblem(),
1854 presolve_.data_.reduced_solution_);
1855 if (!solution_ok) return HighsPostsolveStatus::ReducedSolutionDimenionsError;
1856
1857 if (presolve_.presolve_status_ != HighsPresolveStatus::Reduced &&
1858 presolve_.presolve_status_ != HighsPresolveStatus::ReducedToEmpty)
1859 return HighsPostsolveStatus::NoPostsolve;
1860
1861 // Handle max case.
1862 if (lp_.sense_ == ObjSense::MAXIMIZE) presolve_.negateReducedLpColDuals(true);
1863
1864 // Run postsolve
1865 HighsPostsolveStatus postsolve_status =
1866 presolve_.data_.presolve_[0].postsolve(
1867 presolve_.data_.reduced_solution_, presolve_.data_.reduced_basis_,
1868 presolve_.data_.recovered_solution_,
1869 presolve_.data_.recovered_basis_);
1870
1871 if (postsolve_status != HighsPostsolveStatus::SolutionRecovered)
1872 return postsolve_status;
1873
1874 if (lp_.sense_ == ObjSense::MAXIMIZE)
1875 presolve_.negateReducedLpColDuals(false);
1876
1877 return HighsPostsolveStatus::SolutionRecovered;
1878 }
1879
1880 // The method below runs calls solveLp to solve the LP associated with
1881 // a particular model, integrating the iteration counts into the
1882 // overall values in HighsInfo
runLpSolver(const int model_index,const string message)1883 HighsStatus Highs::runLpSolver(const int model_index, const string message) {
1884 HighsStatus return_status = HighsStatus::OK;
1885 HighsStatus call_status;
1886
1887 // Check that the model index is OK
1888 bool model_index_ok = model_index >= 0 && model_index < (int)hmos_.size();
1889 assert(model_index_ok);
1890 if (!model_index_ok) return HighsStatus::Error;
1891
1892 HighsModelObject& model = hmos_[model_index];
1893
1894 // Transfer the LP solver iteration counts to this model
1895 HighsIterationCounts& iteration_counts = hmos_[model_index].iteration_counts_;
1896 copyHighsIterationCounts(info_, iteration_counts);
1897
1898 // Solve the LP
1899 call_status = solveLp(model, message);
1900 return_status = interpretCallStatus(call_status, return_status, "solveLp");
1901 if (return_status == HighsStatus::Error) return return_status;
1902
1903 // Transfer this model's LP solver iteration counts to HiGHS
1904 copyHighsIterationCounts(iteration_counts, info_);
1905
1906 return returnFromHighs(return_status);
1907 }
1908
writeSolution(const std::string filename,const bool pretty) const1909 HighsStatus Highs::writeSolution(const std::string filename,
1910 const bool pretty) const {
1911 HighsStatus return_status = HighsStatus::OK;
1912 HighsStatus call_status;
1913 HighsLp lp = this->lp_;
1914 HighsBasis basis = this->basis_;
1915 HighsSolution solution = this->solution_;
1916 FILE* file;
1917 bool html;
1918 call_status = openWriteFile(filename, "writeSolution", file, html);
1919 return_status =
1920 interpretCallStatus(call_status, return_status, "openWriteFile");
1921 if (return_status == HighsStatus::Error) return return_status;
1922
1923 // std::cout << "warning: Feature under development" << std::endl;
1924
1925 writeSolutionToFile(file, lp, basis, solution, pretty);
1926
1927 // return HighsStatus::Warning;
1928 return HighsStatus::OK;
1929 }
1930
1931 // Actions to take if there is a new Highs basis
newHighsBasis()1932 void Highs::newHighsBasis() {
1933 if (hmos_.size() > 0) {
1934 // Copy this basis to the HMO basis and clear any simplex basis
1935 hmos_[0].basis_ = this->basis_;
1936 HighsSimplexInterface interface(hmos_[0]);
1937 interface.clearBasis();
1938 }
1939 }
1940
1941 // Ensure that the HiGHS solution and basis have the same size as the
1942 // model, and that the HiGHS basis is kept up-to-date with any solved
1943 // basis
forceHighsSolutionBasisSize()1944 void Highs::forceHighsSolutionBasisSize() {
1945 // Ensure that the HiGHS solution vectors are the right size
1946 solution_.col_value.resize(lp_.numCol_);
1947 solution_.row_value.resize(lp_.numRow_);
1948 solution_.col_dual.resize(lp_.numCol_);
1949 solution_.row_dual.resize(lp_.numRow_);
1950 // Ensure that the HiGHS basis vectors are the right size,
1951 // invalidating the basis if they aren't
1952 if ((int)basis_.col_status.size() != lp_.numCol_) {
1953 basis_.col_status.resize(lp_.numCol_);
1954 basis_.valid_ = false;
1955 }
1956 if ((int)basis_.row_status.size() != lp_.numRow_) {
1957 basis_.row_status.resize(lp_.numRow_);
1958 basis_.valid_ = false;
1959 }
1960 }
1961
getHighsModelStatusAndInfo(const int solved_hmo)1962 bool Highs::getHighsModelStatusAndInfo(const int solved_hmo) {
1963 if (!haveHmo("getHighsModelStatusAndInfo")) return false;
1964
1965 model_status_ = hmos_[solved_hmo].unscaled_model_status_;
1966 scaled_model_status_ = hmos_[solved_hmo].scaled_model_status_;
1967
1968 HighsSolutionParams& solution_params =
1969 hmos_[solved_hmo].unscaled_solution_params_;
1970
1971 info_.primal_status = solution_params.primal_status;
1972 info_.dual_status = solution_params.dual_status;
1973 info_.objective_function_value = solution_params.objective_function_value;
1974 info_.num_primal_infeasibilities = solution_params.num_primal_infeasibilities;
1975 info_.max_primal_infeasibility = solution_params.max_primal_infeasibility;
1976 info_.sum_primal_infeasibilities = solution_params.sum_primal_infeasibilities;
1977 info_.num_dual_infeasibilities = solution_params.num_dual_infeasibilities;
1978 info_.max_dual_infeasibility = solution_params.max_dual_infeasibility;
1979 info_.sum_dual_infeasibilities = solution_params.sum_dual_infeasibilities;
1980 return true;
1981 }
1982
openWriteFile(const string filename,const string method_name,FILE * & file,bool & html) const1983 HighsStatus Highs::openWriteFile(const string filename,
1984 const string method_name, FILE*& file,
1985 bool& html) const {
1986 html = false;
1987 if (filename == "") {
1988 // Empty file name: use stdout
1989 file = stdout;
1990 } else {
1991 file = fopen(filename.c_str(), "w");
1992 if (file == 0) {
1993 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
1994 "Cannot open writeable file \"%s\" in %s",
1995 filename.c_str(), method_name.c_str());
1996 return HighsStatus::Error;
1997 }
1998 const char* dot = strrchr(filename.c_str(), '.');
1999 if (dot && dot != filename) html = strcmp(dot + 1, "html") == 0;
2000 }
2001 return HighsStatus::OK;
2002 }
2003
getUseModelStatus(HighsModelStatus & use_model_status,const double unscaled_primal_feasibility_tolerance,const double unscaled_dual_feasibility_tolerance,const bool rerun_from_logical_basis)2004 HighsStatus Highs::getUseModelStatus(
2005 HighsModelStatus& use_model_status,
2006 const double unscaled_primal_feasibility_tolerance,
2007 const double unscaled_dual_feasibility_tolerance,
2008 const bool rerun_from_logical_basis) {
2009 if (model_status_ != HighsModelStatus::NOTSET) {
2010 use_model_status = model_status_;
2011 } else {
2012 // Handle the case where the status of the unscaled model is not set
2013 HighsStatus return_status = HighsStatus::OK;
2014 HighsStatus call_status;
2015 const double report = false; // true;//
2016 if (unscaledOptimal(unscaled_primal_feasibility_tolerance,
2017 unscaled_dual_feasibility_tolerance, report)) {
2018 use_model_status = HighsModelStatus::OPTIMAL;
2019 } else if (rerun_from_logical_basis) {
2020 std::string save_presolve = options_.presolve;
2021 basis_.valid_ = false;
2022 options_.presolve = on_string;
2023 call_status = run();
2024 return_status = interpretCallStatus(call_status, return_status, "run()");
2025 options_.presolve = save_presolve;
2026 if (return_status == HighsStatus::Error) return return_status;
2027
2028 if (report)
2029 printf(
2030 "Unscaled model status was NOTSET: after running from logical "
2031 "basis it is %s\n",
2032 highsModelStatusToString(model_status_).c_str());
2033
2034 if (model_status_ != HighsModelStatus::NOTSET) {
2035 use_model_status = model_status_;
2036 } else if (unscaledOptimal(unscaled_primal_feasibility_tolerance,
2037 unscaled_dual_feasibility_tolerance, report)) {
2038 use_model_status = HighsModelStatus::OPTIMAL;
2039 }
2040 } else {
2041 // Nothing to be done: use original unscaled model status
2042 use_model_status = model_status_;
2043 }
2044 }
2045 return HighsStatus::OK;
2046 }
2047
unscaledOptimal(const double unscaled_primal_feasibility_tolerance,const double unscaled_dual_feasibility_tolerance,const bool report)2048 bool Highs::unscaledOptimal(const double unscaled_primal_feasibility_tolerance,
2049 const double unscaled_dual_feasibility_tolerance,
2050 const bool report) {
2051 if (scaled_model_status_ == HighsModelStatus::OPTIMAL) {
2052 const double max_primal_infeasibility = info_.max_primal_infeasibility;
2053 const double max_dual_infeasibility = info_.max_dual_infeasibility;
2054 if (report)
2055 printf(
2056 "Scaled model status is OPTIMAL: max unscaled (primal / dual) "
2057 "infeasibilities are (%g / %g)\n",
2058 max_primal_infeasibility, max_dual_infeasibility);
2059 if ((max_primal_infeasibility > unscaled_primal_feasibility_tolerance) ||
2060 (max_dual_infeasibility > unscaled_dual_feasibility_tolerance)) {
2061 printf(
2062 "Use model status of NOTSET since max unscaled (primal / dual) "
2063 "infeasibilities are (%g / %g)\n",
2064 max_primal_infeasibility, max_dual_infeasibility);
2065 } else {
2066 if (report)
2067 printf(
2068 "Set unscaled model status to OPTIMAL since unscaled "
2069 "infeasibilities are tolerable\n");
2070 return true;
2071 }
2072 }
2073 return false;
2074 }
2075
haveHmo(const string method_name) const2076 bool Highs::haveHmo(const string method_name) const {
2077 bool have_hmo = hmos_.size() > 0;
2078 #ifdef HiGHSDEV
2079 if (!have_hmo)
2080 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
2081 "Method %s called without any HighsModelObject",
2082 method_name.c_str());
2083 #endif
2084 assert(have_hmo);
2085 return have_hmo;
2086 }
2087
clearModelStatus()2088 void Highs::clearModelStatus() {
2089 model_status_ = HighsModelStatus::NOTSET;
2090 scaled_model_status_ = HighsModelStatus::NOTSET;
2091 }
2092
clearSolution()2093 void Highs::clearSolution() {
2094 info_.primal_status = (int)PrimalDualStatus::STATUS_NOTSET;
2095 info_.dual_status = (int)PrimalDualStatus::STATUS_NOTSET;
2096 clearSolutionUtil(solution_);
2097 }
2098
clearBasis()2099 void Highs::clearBasis() { clearBasisUtil(basis_); }
2100
clearInfo()2101 void Highs::clearInfo() { info_.clear(); }
2102
2103 // Applies checks before returning from run()
returnFromRun(const HighsStatus run_return_status)2104 HighsStatus Highs::returnFromRun(const HighsStatus run_return_status) {
2105 bool have_solution = false;
2106 HighsStatus return_status = run_return_status;
2107 if (hmos_.size() == 0) {
2108 // No model has been loaded: ensure that the status, solution,
2109 // basis and info associated with any previous model are cleared
2110 clearSolver();
2111 return returnFromHighs(return_status);
2112 } else {
2113 // A model has been loaded: remove any additional HMO created when solving
2114 if (hmos_.size() > 1) hmos_.pop_back();
2115 // There should be only one entry in hmos_
2116 assert((int)hmos_.size() == 1);
2117 // Make sure that the unscaled status, solution, basis and info
2118 // are consistent with the scaled status
2119 #ifdef HiGHSDEV
2120 reportModelStatusSolutionBasis("returnFromRun(HiGHS)");
2121 reportModelStatusSolutionBasis("returnFromRun(HMO_0)", 0);
2122 #endif
2123 switch (scaled_model_status_) {
2124 // First consider the error returns
2125 case HighsModelStatus::NOTSET:
2126 case HighsModelStatus::LOAD_ERROR:
2127 case HighsModelStatus::MODEL_ERROR:
2128 case HighsModelStatus::PRESOLVE_ERROR:
2129 case HighsModelStatus::SOLVE_ERROR:
2130 case HighsModelStatus::POSTSOLVE_ERROR:
2131 clearSolver();
2132 assert(return_status == HighsStatus::Error);
2133 break;
2134
2135 // Then consider the OK returns
2136 case HighsModelStatus::MODEL_EMPTY:
2137 clearSolution();
2138 clearBasis();
2139 clearInfo();
2140 assert(model_status_ == scaled_model_status_);
2141 assert(return_status == HighsStatus::OK);
2142 break;
2143
2144 case HighsModelStatus::PRIMAL_INFEASIBLE:
2145 clearSolution();
2146 // May have a basis, according to whether infeasibility was
2147 // detected in presolve or solve
2148 assert(model_status_ == scaled_model_status_);
2149 assert(return_status == HighsStatus::OK);
2150 break;
2151
2152 case HighsModelStatus::PRIMAL_UNBOUNDED:
2153 clearSolution();
2154 // May have a basis, according to whether infeasibility was
2155 // detected in presolve or solve
2156 clearInfo();
2157 assert(model_status_ == scaled_model_status_);
2158 assert(return_status == HighsStatus::OK);
2159 break;
2160
2161 case HighsModelStatus::OPTIMAL:
2162 have_solution = true;
2163 // The following is an aspiration
2164 // assert(info_.primal_status ==
2165 // (int)PrimalDualStatus::STATUS_FEASIBLE_POINT);
2166 // assert(info_.dual_status ==
2167 // (int)PrimalDualStatus::STATUS_FEASIBLE_POINT);
2168 assert(model_status_ == HighsModelStatus::NOTSET ||
2169 model_status_ == HighsModelStatus::OPTIMAL);
2170 assert(return_status == HighsStatus::OK);
2171 break;
2172
2173 case HighsModelStatus::PRIMAL_DUAL_INFEASIBLE:
2174 clearSolution();
2175 // May have a basis, according to whether infeasibility was
2176 // detected in presolve or solve
2177 clearInfo();
2178 assert(model_status_ == scaled_model_status_);
2179 assert(return_status == HighsStatus::OK);
2180 break;
2181
2182 case HighsModelStatus::REACHED_DUAL_OBJECTIVE_VALUE_UPPER_BOUND:
2183 clearSolution();
2184 clearBasis();
2185 clearInfo();
2186 assert(model_status_ == scaled_model_status_);
2187 assert(return_status == HighsStatus::OK);
2188 break;
2189
2190 // Finally consider the warning returns
2191 case HighsModelStatus::REACHED_TIME_LIMIT:
2192 case HighsModelStatus::REACHED_ITERATION_LIMIT:
2193 clearSolution();
2194 clearBasis();
2195 clearInfo();
2196 assert(model_status_ == scaled_model_status_);
2197 assert(return_status == HighsStatus::Warning);
2198 break;
2199 case HighsModelStatus::DUAL_INFEASIBLE:
2200 clearSolution();
2201 // May have a basis, according to whether infeasibility was
2202 // detected in presolve or solve
2203 clearInfo();
2204 assert(model_status_ == scaled_model_status_);
2205 assert(return_status == HighsStatus::Warning);
2206 break;
2207 }
2208 }
2209 if (have_solution) debugSolutionRightSize(options_, lp_, solution_);
2210 bool have_basis = basis_.valid_;
2211 if (have_basis) {
2212 if (debugBasisRightSize(options_, lp_, basis_) ==
2213 HighsDebugStatus::LOGICAL_ERROR)
2214 return_status = HighsStatus::Error;
2215 }
2216
2217 if (have_solution && have_basis) {
2218 if (debugHighsBasicSolution("Return from run()", options_, lp_, basis_,
2219 solution_, info_, model_status_) ==
2220 HighsDebugStatus::LOGICAL_ERROR)
2221 return_status = HighsStatus::Error;
2222 }
2223 return returnFromHighs(return_status);
2224 }
2225
2226 // Applies checks before returning from HiGHS
returnFromHighs(HighsStatus highs_return_status)2227 HighsStatus Highs::returnFromHighs(HighsStatus highs_return_status) {
2228 HighsStatus return_status = highs_return_status;
2229
2230 forceHighsSolutionBasisSize();
2231
2232 const bool consistent = debugBasisConsistent(options_, lp_, basis_) !=
2233 HighsDebugStatus::LOGICAL_ERROR;
2234 if (!consistent) {
2235 HighsLogMessage(
2236 options_.logfile, HighsMessageType::ERROR,
2237 "returnFromHighs: Supposed to be a HiGHS basis, but not consistent");
2238 assert(consistent);
2239 return_status = HighsStatus::Error;
2240 }
2241
2242 if (hmos_.size()) {
2243 const bool simplex_lp_ok =
2244 debugSimplexLp(hmos_[0]) != HighsDebugStatus::LOGICAL_ERROR;
2245 if (!simplex_lp_ok) {
2246 HighsLogMessage(options_.logfile, HighsMessageType::ERROR,
2247 "returnFromHighs: Simplex LP not OK");
2248 assert(simplex_lp_ok);
2249 return_status = HighsStatus::Error;
2250 }
2251 }
2252
2253 return return_status;
2254 }
underDevelopmentLogMessage(const string method_name)2255 void Highs::underDevelopmentLogMessage(const string method_name) {
2256 HighsLogMessage(
2257 options_.logfile, HighsMessageType::WARNING,
2258 "Method %s is still under development and behaviour may be unpredictable",
2259 method_name.c_str());
2260 }
2261
getPresolveReductionCounts(int & rows,int & cols,int & nnz) const2262 void Highs::getPresolveReductionCounts(int& rows, int& cols, int& nnz) const {
2263 rows = presolve_.info_.n_rows_removed;
2264 cols = presolve_.info_.n_cols_removed;
2265 nnz = presolve_.info_.n_nnz_removed;
2266 }
2267