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/HSimplexDebug.cpp
11 * @brief
12 * @author Julian Hall, Ivet Galabova, Qi Huangfu and Michael Feldmeier
13 */
14
15 #include "simplex/HSimplexDebug.h"
16
17 #include <string>
18
19 #include "lp_data/HighsDebug.h"
20 #include "lp_data/HighsLpUtils.h"
21 #include "lp_data/HighsModelUtils.h"
22 #include "lp_data/HighsSolutionDebug.h"
23 #include "simplex/HDualRow.h"
24 #include "simplex/HFactorDebug.h"
25 #include "simplex/HSimplex.h"
26 #include "simplex/SimplexTimer.h"
27
28 const double excessive_absolute_primal_norm = 1e12;
29 const double excessive_relative_primal_norm = 1e6;
30 const double large_absolute_primal_norm = sqrt(excessive_absolute_primal_norm);
31 const double large_relative_primal_norm = sqrt(excessive_relative_primal_norm);
32
33 const double excessive_absolute_nonbasic_dual_norm = 1e12;
34 const double excessive_relative_nonbasic_dual_norm = 1e6;
35 const double large_absolute_nonbasic_dual_norm =
36 sqrt(excessive_absolute_nonbasic_dual_norm);
37 const double large_relative_nonbasic_dual_norm =
38 sqrt(excessive_relative_nonbasic_dual_norm);
39
40 const double large_absolute_basic_dual_norm = 1e-12;
41 const double large_relative_basic_dual_norm = 1e-14;
42 const double excessive_absolute_basic_dual_norm =
43 sqrt(large_absolute_basic_dual_norm);
44 const double excessive_relative_basic_dual_norm =
45 sqrt(large_relative_basic_dual_norm);
46
47 const double computed_primal_excessive_absolute_norm =
48 excessive_absolute_primal_norm;
49 const double computed_primal_excessive_relative_norm =
50 excessive_relative_primal_norm;
51 const double computed_primal_large_absolute_norm = large_absolute_primal_norm;
52 const double computed_primal_large_relative_norm = large_relative_primal_norm;
53
54 const double computed_dual_excessive_absolute_nonbasic_dual_norm =
55 excessive_absolute_nonbasic_dual_norm;
56 const double computed_dual_excessive_relative_nonbasic_dual_norm =
57 excessive_relative_nonbasic_dual_norm;
58 const double computed_dual_large_absolute_nonbasic_dual_norm =
59 large_absolute_nonbasic_dual_norm;
60 const double computed_dual_large_relative_nonbasic_dual_norm =
61 large_relative_nonbasic_dual_norm;
62
63 const double computed_dual_excessive_absolute_basic_dual_norm =
64 excessive_absolute_basic_dual_norm;
65 const double computed_dual_excessive_relative_basic_dual_norm =
66 excessive_relative_basic_dual_norm;
67 const double computed_dual_large_absolute_basic_dual_norm =
68 large_absolute_basic_dual_norm;
69 const double computed_dual_large_relative_basic_dual_norm =
70 large_relative_basic_dual_norm;
71
72 const double computed_dual_small_relative_nonbasic_dual_change_norm = 1e-12;
73 const double computed_dual_large_relative_nonbasic_dual_change_norm =
74 sqrt(computed_dual_small_relative_nonbasic_dual_change_norm);
75 const double computed_dual_small_absolute_nonbasic_dual_change_norm = 1e-6;
76 const double computed_dual_large_absolute_nonbasic_dual_change_norm =
77 sqrt(computed_dual_small_absolute_nonbasic_dual_change_norm);
78
79 const double updated_objective_small_relative_error = 1e-12;
80 const double updated_objective_large_relative_error =
81 sqrt(updated_objective_small_relative_error);
82 const double updated_objective_small_absolute_error = 1e-6;
83 const double updated_objective_large_absolute_error =
84 sqrt(updated_objective_small_absolute_error);
85
86 const double excessive_basis_condition = 1e16;
87 const double large_basis_condition = sqrt(excessive_basis_condition);
88 const double fair_basis_condition = sqrt(large_basis_condition);
89
90 const double cleanup_large_absolute_nonbasic_dual_change_norm = 1e-12;
91 const double cleanup_large_relative_nonbasic_dual_change_norm = 1e-6;
92 const double cleanup_excessive_absolute_nonbasic_dual_change_norm =
93 sqrt(cleanup_large_absolute_nonbasic_dual_change_norm);
94 const double cleanup_excessive_relative_nonbasic_dual_change_norm =
95 sqrt(cleanup_large_relative_nonbasic_dual_change_norm);
96
97 const double freelist_excessive_pct_num_entries = 25.0;
98 const double freelist_large_pct_num_entries = 10.0;
99 const double freelist_fair_pct_num_entries = 1.0;
100
debugSimplexLp(const HighsModelObject & highs_model_object)101 HighsDebugStatus debugSimplexLp(const HighsModelObject& highs_model_object) {
102 // Non-trivially expensive check that the .simplex_lp, if valid is .lp scaled
103 // according to .scale
104 const HighsSimplexLpStatus& simplex_lp_status =
105 highs_model_object.simplex_lp_status_;
106 if (!simplex_lp_status.valid ||
107 highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
108 return HighsDebugStatus::NOT_CHECKED;
109 HighsDebugStatus return_status = HighsDebugStatus::OK;
110 const HighsOptions& options = highs_model_object.options_;
111 const HighsLp& lp = highs_model_object.lp_;
112 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
113 const HighsScale& scale = highs_model_object.scale_;
114 const HFactor& factor = highs_model_object.factor_;
115
116 bool right_size = true;
117 right_size = (int)scale.col_.size() == lp.numCol_ && right_size;
118 right_size = (int)scale.row_.size() == lp.numRow_ && right_size;
119 if (!right_size) {
120 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
121 "scale size error");
122 assert(right_size);
123 return_status = HighsDebugStatus::LOGICAL_ERROR;
124 }
125 // Take a copy of the original LP
126 HighsLp check_lp = lp;
127 if (applyScalingToLp(options, check_lp, scale) != HighsStatus::OK) {
128 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
129 "debugSimplexLp: Error scaling check LP");
130 return HighsDebugStatus::LOGICAL_ERROR;
131 }
132 const bool simplex_lp_data_ok = check_lp == simplex_lp;
133 if (!simplex_lp_data_ok) {
134 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
135 "debugSimplexLp: Check LP and simplex LP not equal");
136 assert(simplex_lp_data_ok);
137 return_status = HighsDebugStatus::LOGICAL_ERROR;
138 }
139
140 if (simplex_lp_status.has_basis) {
141 const bool simplex_basis_correct =
142 debugDebugToHighsStatus(debugSimplexBasisCorrect(highs_model_object)) !=
143 HighsStatus::Error;
144 if (!simplex_basis_correct) {
145 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
146 "Supposed to be a Simplex basis, but incorrect");
147 assert(simplex_basis_correct);
148 return_status = HighsDebugStatus::LOGICAL_ERROR;
149 }
150 }
151
152 if (simplex_lp_status.has_invert) {
153 const bool invert_ok = debugDebugToHighsStatus(debugCheckInvert(
154 options, factor)) != HighsStatus::Error;
155 if (!invert_ok) {
156 HighsLogMessage(
157 options.logfile, HighsMessageType::ERROR,
158 "Supposed to be a Simplex basis inverse, but too inaccurate");
159 assert(invert_ok);
160 return_status = HighsDebugStatus::LOGICAL_ERROR;
161 }
162 }
163 return return_status;
164 }
165
debugBasisRightSize(const HighsOptions & options,const HighsLp & simplex_lp,const SimplexBasis & simplex_basis)166 HighsDebugStatus debugBasisRightSize(const HighsOptions& options,
167 const HighsLp& simplex_lp,
168 const SimplexBasis& simplex_basis) {
169 // Cheap analysis of a Simplex basis, checking vector sizes
170 if (options.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
171 return HighsDebugStatus::NOT_CHECKED;
172 HighsDebugStatus return_status = HighsDebugStatus::OK;
173 bool right_size = isBasisRightSize(simplex_lp, simplex_basis);
174 if (!right_size) {
175 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
176 "Simplex basis size error");
177 assert(right_size);
178 return_status = HighsDebugStatus::LOGICAL_ERROR;
179 }
180 return return_status;
181 }
182
debugSimplexInfoBasisRightSize(const HighsModelObject & highs_model_object)183 HighsDebugStatus debugSimplexInfoBasisRightSize(
184 const HighsModelObject& highs_model_object) {
185 // Trivially cheap check of dimensions and sizes
186 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
187 return HighsDebugStatus::NOT_CHECKED;
188
189 const HighsOptions& options = highs_model_object.options_;
190 const HighsLp& lp = highs_model_object.lp_;
191 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
192 const HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
193 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
194
195 int numCol = lp.numCol_;
196 int numRow = lp.numRow_;
197 int numTot = numCol + numRow;
198 HighsDebugStatus return_status = HighsDebugStatus::OK;
199 bool dimension_ok =
200 numCol == simplex_lp.numCol_ && numRow == simplex_lp.numRow_;
201 assert(dimension_ok);
202 if (!dimension_ok) {
203 HighsLogMessage(
204 options.logfile, HighsMessageType::ERROR,
205 "LP-SimplexLP dimension incompatibility (%d, %d) != (%d, %d)", numCol,
206 simplex_lp.numCol_, numRow, simplex_lp.numRow_);
207 return_status = HighsDebugStatus::LOGICAL_ERROR;
208 }
209 bool right_size = true;
210 right_size = (int)simplex_info.workCost_.size() == numTot && right_size;
211 right_size = (int)simplex_info.workDual_.size() == numTot && right_size;
212 right_size = (int)simplex_info.workShift_.size() == numTot && right_size;
213 right_size = (int)simplex_info.workLower_.size() == numTot && right_size;
214 right_size = (int)simplex_info.workUpper_.size() == numTot && right_size;
215 right_size = (int)simplex_info.workRange_.size() == numTot && right_size;
216 right_size = (int)simplex_info.workValue_.size() == numTot && right_size;
217 if (!right_size) {
218 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
219 "simplex_info work vector size error");
220 assert(right_size);
221 return_status = HighsDebugStatus::LOGICAL_ERROR;
222 }
223 if (debugBasisRightSize(options, simplex_lp, simplex_basis) !=
224 HighsDebugStatus::OK)
225 return_status = HighsDebugStatus::LOGICAL_ERROR;
226
227 return return_status;
228 }
229
debugComputePrimal(const HighsModelObject & highs_model_object,const std::vector<double> & primal_rhs)230 HighsDebugStatus debugComputePrimal(const HighsModelObject& highs_model_object,
231 const std::vector<double>& primal_rhs) {
232 // Non-trivially expensive analysis of computed primal values.
233 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
234 return HighsDebugStatus::NOT_CHECKED;
235 HighsDebugStatus return_status = HighsDebugStatus::NOT_CHECKED;
236 const std::vector<double>& primal_value =
237 highs_model_object.simplex_info_.baseValue_;
238
239 int num_row = highs_model_object.simplex_lp_.numRow_;
240
241 // Use the size of the RHS to determine whether to use it
242 const bool have_primal_rhs = (int)primal_rhs.size() == num_row;
243
244 double primal_rhs_norm = 0;
245 if (have_primal_rhs) {
246 for (int iRow = 0; iRow < num_row; iRow++)
247 primal_rhs_norm += fabs(primal_rhs[iRow]);
248 }
249 double computed_absolute_primal_norm = 0;
250 for (int iRow = 0; iRow < num_row; iRow++)
251 computed_absolute_primal_norm += fabs(primal_value[iRow]);
252
253 std::string value_adjective;
254 int report_level;
255 return_status = HighsDebugStatus::OK;
256 double computed_relative_primal_norm;
257 if (primal_rhs_norm) {
258 computed_relative_primal_norm =
259 computed_absolute_primal_norm / primal_rhs_norm;
260 } else {
261 computed_relative_primal_norm = -1;
262 }
263 if (computed_relative_primal_norm > computed_primal_excessive_relative_norm ||
264 computed_absolute_primal_norm > computed_primal_excessive_absolute_norm) {
265 value_adjective = "Excessive";
266 report_level = ML_ALWAYS;
267 return_status = HighsDebugStatus::ERROR;
268 } else if (computed_relative_primal_norm >
269 computed_primal_large_relative_norm ||
270 computed_absolute_primal_norm >
271 computed_primal_large_absolute_norm) {
272 value_adjective = "Large";
273 report_level = ML_DETAILED;
274 return_status = HighsDebugStatus::WARNING;
275 } else {
276 value_adjective = "SMALL";
277 report_level = ML_VERBOSE;
278 }
279 HighsPrintMessage(
280 highs_model_object.options_.output,
281 highs_model_object.options_.message_level, report_level,
282 "ComputePrimal: %-9s absolute (%9.4g) or relative (%9.4g) norm of "
283 "primal values\n",
284 value_adjective.c_str(), computed_absolute_primal_norm,
285 computed_relative_primal_norm);
286 return return_status;
287 }
debugComputeDual(const HighsModelObject & highs_model_object,const std::vector<double> & previous_dual,const std::vector<double> & basic_costs,const std::vector<double> & row_dual)288 HighsDebugStatus debugComputeDual(const HighsModelObject& highs_model_object,
289 const std::vector<double>& previous_dual,
290 const std::vector<double>& basic_costs,
291 const std::vector<double>& row_dual) {
292 // Non-trivially expensive analysis of computed dual values.
293 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
294 return HighsDebugStatus::NOT_CHECKED;
295 HighsDebugStatus return_status = HighsDebugStatus::NOT_CHECKED;
296 const std::vector<double>& new_dual =
297 highs_model_object.simplex_info_.workDual_;
298
299 int num_row = highs_model_object.simplex_lp_.numRow_;
300 int num_col = highs_model_object.simplex_lp_.numCol_;
301
302 const bool have_basic_costs = (int)basic_costs.size() == num_row;
303 const bool have_row_dual = (int)row_dual.size() == num_row;
304 const bool have_previous_dual =
305 (int)previous_dual.size() == num_col + num_row;
306
307 double basic_costs_norm = 0;
308 if (have_basic_costs) {
309 for (int iRow = 0; iRow < num_row; iRow++)
310 basic_costs_norm += fabs(basic_costs[iRow]);
311 }
312 double row_dual_norm = 0;
313 if (have_row_dual) {
314 for (int iRow = 0; iRow < num_row; iRow++)
315 row_dual_norm += fabs(row_dual[iRow]);
316 }
317 double computed_dual_absolute_basic_dual_norm = 0;
318 double computed_dual_absolute_nonbasic_dual_norm = 0;
319 for (int iVar = 0; iVar < num_row + num_col; iVar++) {
320 if (!highs_model_object.simplex_basis_.nonbasicFlag_[iVar]) {
321 computed_dual_absolute_basic_dual_norm += fabs(new_dual[iVar]);
322 continue;
323 }
324 computed_dual_absolute_nonbasic_dual_norm += fabs(new_dual[iVar]);
325 }
326 std::string value_adjective;
327 int report_level;
328 return_status = HighsDebugStatus::OK;
329 // Comment on the norm of the basic costs being zero
330 if (have_basic_costs && !basic_costs_norm) {
331 HighsLogMessage(
332 highs_model_object.options_.logfile, HighsMessageType::WARNING,
333 "ComputeDual: basic cost norm is = %9.4g", basic_costs_norm);
334 return_status = HighsDebugStatus::WARNING;
335 }
336 // Comment on the norm of the nonbasic duals being zero
337 if (!computed_dual_absolute_nonbasic_dual_norm) {
338 HighsLogMessage(highs_model_object.options_.logfile,
339 HighsMessageType::WARNING,
340 "ComputeDual: nonbasic dual norm is = %9.4g",
341 computed_dual_absolute_nonbasic_dual_norm);
342 return_status = HighsDebugStatus::WARNING;
343 }
344
345 // Comment on the norm of basic duals (relative to the norm of the
346 // basic costs) which, as c_B-BB^{-1}c_B, should be zero
347 double computed_dual_relative_basic_dual_norm;
348 if (basic_costs_norm) {
349 computed_dual_relative_basic_dual_norm =
350 computed_dual_absolute_basic_dual_norm / basic_costs_norm;
351 } else {
352 computed_dual_relative_basic_dual_norm = -1;
353 }
354 if (computed_dual_relative_basic_dual_norm >
355 computed_dual_excessive_relative_basic_dual_norm ||
356 computed_dual_absolute_basic_dual_norm >
357 computed_dual_excessive_absolute_basic_dual_norm) {
358 value_adjective = "Excessive";
359 report_level = ML_ALWAYS;
360 return_status = HighsDebugStatus::ERROR;
361 } else if (computed_dual_relative_basic_dual_norm >
362 computed_dual_large_relative_basic_dual_norm ||
363 computed_dual_absolute_basic_dual_norm >
364 computed_dual_large_absolute_basic_dual_norm) {
365 value_adjective = "Large";
366 report_level = ML_DETAILED;
367 return_status = HighsDebugStatus::WARNING;
368 } else {
369 value_adjective = "OK";
370 report_level = ML_VERBOSE;
371 }
372 HighsPrintMessage(
373 highs_model_object.options_.output,
374 highs_model_object.options_.message_level, report_level,
375 "ComputeDual: %-9s absolute (%9.4g) or relative (%9.4g) norm of "
376 " basic dual values\n",
377 value_adjective.c_str(), computed_dual_absolute_basic_dual_norm,
378 computed_dual_relative_basic_dual_norm);
379 // Comment on the norm of nonbasic duals relative to the norm of the
380 // basic costs
381 double computed_dual_relative_nonbasic_dual_norm;
382 if (basic_costs_norm) {
383 computed_dual_relative_nonbasic_dual_norm =
384 computed_dual_absolute_nonbasic_dual_norm / basic_costs_norm;
385 } else {
386 computed_dual_relative_nonbasic_dual_norm = -1;
387 }
388 if (computed_dual_relative_nonbasic_dual_norm >
389 computed_dual_excessive_relative_nonbasic_dual_norm ||
390 computed_dual_absolute_nonbasic_dual_norm >
391 computed_dual_excessive_absolute_nonbasic_dual_norm) {
392 value_adjective = "Excessive";
393 report_level = ML_ALWAYS;
394 return_status = HighsDebugStatus::ERROR;
395 } else if (computed_dual_relative_nonbasic_dual_norm >
396 computed_dual_large_relative_nonbasic_dual_norm ||
397 computed_dual_absolute_nonbasic_dual_norm >
398 computed_dual_large_absolute_nonbasic_dual_norm) {
399 value_adjective = "Large";
400 report_level = ML_DETAILED;
401 return_status = HighsDebugStatus::WARNING;
402 } else {
403 value_adjective = "OK";
404 report_level = ML_VERBOSE;
405 }
406 HighsPrintMessage(
407 highs_model_object.options_.output,
408 highs_model_object.options_.message_level, report_level,
409 "ComputeDual: %-9s absolute (%9.4g) or relative (%9.4g) norm of "
410 "nonbasic dual values\n",
411 value_adjective.c_str(), computed_dual_absolute_nonbasic_dual_norm,
412 computed_dual_relative_nonbasic_dual_norm);
413 double report_basic_costs_norm = -1;
414 if (basic_costs_norm) report_basic_costs_norm = basic_costs_norm;
415 double report_row_dual_norm = -1;
416 if (row_dual_norm) report_row_dual_norm = row_dual_norm;
417 HighsPrintMessage(highs_model_object.options_.output,
418 highs_model_object.options_.message_level, report_level,
419 "ComputeDual: B.pi=c_B has |c_B| = %9.4g; |pi| = %9.4g; "
420 "|pi^TA-c| = [basic %9.4g; nonbasic %9.4g]\n",
421 report_basic_costs_norm, report_row_dual_norm,
422 computed_dual_absolute_basic_dual_norm,
423 computed_dual_absolute_nonbasic_dual_norm);
424 if (have_previous_dual) {
425 // Comment on the change in the dual values
426 std::string change_adjective;
427 double computed_dual_absolute_nonbasic_dual_change_norm = 0;
428 for (int iVar = 0; iVar < num_row + num_col; iVar++) {
429 if (!highs_model_object.simplex_basis_.nonbasicFlag_[iVar]) continue;
430 computed_dual_absolute_nonbasic_dual_change_norm +=
431 fabs(new_dual[iVar] - previous_dual[iVar]);
432 }
433 double computed_dual_relative_nonbasic_dual_change_norm;
434 if (computed_dual_absolute_nonbasic_dual_norm) {
435 computed_dual_relative_nonbasic_dual_change_norm =
436 computed_dual_absolute_nonbasic_dual_change_norm /
437 computed_dual_absolute_nonbasic_dual_norm;
438 } else {
439 computed_dual_relative_nonbasic_dual_change_norm = -1;
440 }
441 if (computed_dual_relative_nonbasic_dual_change_norm >
442 computed_dual_large_relative_nonbasic_dual_change_norm ||
443 computed_dual_absolute_nonbasic_dual_change_norm >
444 computed_dual_large_absolute_nonbasic_dual_change_norm) {
445 change_adjective = "Large";
446 report_level = ML_ALWAYS;
447 return_status = HighsDebugStatus::WARNING;
448 } else if (computed_dual_relative_nonbasic_dual_change_norm >
449 computed_dual_small_relative_nonbasic_dual_change_norm ||
450 computed_dual_absolute_nonbasic_dual_change_norm >
451 computed_dual_small_absolute_nonbasic_dual_change_norm) {
452 change_adjective = "Small";
453 report_level = ML_DETAILED;
454 return_status = HighsDebugStatus::WARNING;
455 } else {
456 change_adjective = "OK";
457 report_level = ML_VERBOSE;
458 }
459 HighsPrintMessage(highs_model_object.options_.output,
460 highs_model_object.options_.message_level, report_level,
461 "ComputeDual: %-9s absolute (%9.4g) or relative "
462 "(%9.4g) nonbasic dual change\n",
463 change_adjective.c_str(),
464 computed_dual_absolute_nonbasic_dual_change_norm,
465 computed_dual_relative_nonbasic_dual_change_norm);
466 }
467 return return_status;
468 }
469
debugSimplexDualFeasibility(const HighsModelObject & highs_model_object,const std::string message,const bool force)470 HighsDebugStatus debugSimplexDualFeasibility(
471 const HighsModelObject& highs_model_object, const std::string message,
472 const bool force) {
473 // Non-trivially expensive check of dual feasibility.
474 if (highs_model_object.options_.highs_debug_level <
475 HIGHS_DEBUG_LEVEL_COSTLY &&
476 !force)
477 return HighsDebugStatus::NOT_CHECKED;
478 if (force)
479 HighsPrintMessage(highs_model_object.options_.output, 1, 1,
480 "SmplxDuFeas: Forcing debug\n");
481
482 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
483 const HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
484 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
485 double scaled_dual_feasibility_tolerance =
486 highs_model_object.scaled_solution_params_.dual_feasibility_tolerance;
487
488 int num_dual_infeasibilities = 0;
489 double sum_dual_infeasibilities = 0;
490 double max_dual_infeasibility = 0;
491 for (int iVar = 0; iVar < simplex_lp.numCol_ + simplex_lp.numRow_; iVar++) {
492 if (!simplex_basis.nonbasicFlag_[iVar]) continue;
493 // Nonbasic column
494 const double dual = simplex_info.workDual_[iVar];
495 const double lower = simplex_info.workLower_[iVar];
496 const double upper = simplex_info.workUpper_[iVar];
497 double dual_infeasibility = 0;
498 if (highs_isInfinity(-lower) && highs_isInfinity(upper)) {
499 // Free: any nonzero dual value is infeasible
500 dual_infeasibility = fabs(dual);
501 } else {
502 // Not free: any dual infeasibility is given by the dual value
503 // signed by nonbasicMove
504 dual_infeasibility = -simplex_basis.nonbasicMove_[iVar] * dual;
505 }
506 if (dual_infeasibility > 0) {
507 if (dual_infeasibility >= scaled_dual_feasibility_tolerance)
508 num_dual_infeasibilities++;
509 max_dual_infeasibility =
510 std::max(dual_infeasibility, max_dual_infeasibility);
511 sum_dual_infeasibilities += dual_infeasibility;
512 }
513 }
514 if (num_dual_infeasibilities) {
515 HighsPrintMessage(highs_model_object.options_.output,
516 highs_model_object.options_.message_level, ML_ALWAYS,
517 "SmplxDuFeas: num/max/sum simplex dual infeasibilities "
518 "= %d / %g / %g - %s\n",
519 num_dual_infeasibilities, max_dual_infeasibility,
520 sum_dual_infeasibilities, message.c_str());
521 return HighsDebugStatus::LOGICAL_ERROR;
522 }
523 return HighsDebugStatus::OK;
524 }
525
debugUpdatedObjectiveValue(HighsModelObject & highs_model_object,const SimplexAlgorithm algorithm,const int phase,const std::string message,const bool force)526 HighsDebugStatus debugUpdatedObjectiveValue(
527 HighsModelObject& highs_model_object, const SimplexAlgorithm algorithm,
528 const int phase, const std::string message, const bool force) {
529 // Non-trivially expensive check of updated objective value. Computes the
530 // exact objective value
531 if (highs_model_object.options_.highs_debug_level <
532 HIGHS_DEBUG_LEVEL_COSTLY &&
533 !force)
534 return HighsDebugStatus::NOT_CHECKED;
535 HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
536
537 static bool have_previous_exact_primal_objective_value;
538 static double previous_exact_primal_objective_value;
539 static double previous_updated_primal_objective_value;
540 static double updated_primal_objective_correction;
541
542 static bool have_previous_exact_dual_objective_value;
543 static double previous_exact_dual_objective_value;
544 static double previous_updated_dual_objective_value;
545 static double updated_dual_objective_correction;
546 if (phase < 0) {
547 if (algorithm == SimplexAlgorithm::PRIMAL) {
548 have_previous_exact_primal_objective_value = false;
549 } else {
550 have_previous_exact_dual_objective_value = false;
551 }
552 return HighsDebugStatus::OK;
553 }
554 double exact_objective_value;
555 double updated_objective_value;
556 bool have_previous_exact_objective_value;
557 // Assign values to prevent compiler warning
558 double previous_exact_objective_value = 0;
559 double previous_updated_objective_value = 0;
560 double updated_objective_correction = 0;
561 std::string algorithm_name;
562 if (algorithm == SimplexAlgorithm::PRIMAL) {
563 algorithm_name = "primal";
564 have_previous_exact_objective_value =
565 have_previous_exact_primal_objective_value;
566 if (have_previous_exact_objective_value) {
567 previous_exact_objective_value = previous_exact_primal_objective_value;
568 previous_updated_objective_value =
569 previous_updated_primal_objective_value;
570 updated_objective_correction = updated_primal_objective_correction;
571 }
572 updated_objective_value = simplex_info.updated_primal_objective_value;
573 // Save the current objective value so that it can be recovered
574 // after calling computePrimalObjectiveValue
575 double save_objective_value = simplex_info.primal_objective_value;
576 computePrimalObjectiveValue(highs_model_object);
577 exact_objective_value = simplex_info.primal_objective_value;
578 simplex_info.primal_objective_value = save_objective_value;
579 } else {
580 algorithm_name = "dual";
581 have_previous_exact_objective_value =
582 have_previous_exact_dual_objective_value;
583 if (have_previous_exact_objective_value) {
584 previous_exact_objective_value = previous_exact_dual_objective_value;
585 previous_updated_objective_value = previous_updated_dual_objective_value;
586 updated_objective_correction = updated_dual_objective_correction;
587 }
588 updated_objective_value = simplex_info.updated_dual_objective_value;
589 // Save the current objective value so that it can be recovered
590 // after calling computeDualObjectiveValue
591 double save_objective_value = simplex_info.dual_objective_value;
592 computeDualObjectiveValue(highs_model_object, phase);
593 exact_objective_value = simplex_info.dual_objective_value;
594 simplex_info.dual_objective_value = save_objective_value;
595 }
596 double change_in_objective_value = 0;
597 double change_in_updated_objective_value = 0;
598 if (have_previous_exact_objective_value) {
599 change_in_objective_value =
600 exact_objective_value - previous_exact_objective_value;
601 change_in_updated_objective_value =
602 updated_objective_value - previous_updated_objective_value;
603 updated_objective_value += updated_objective_correction;
604 } else {
605 updated_objective_correction = 0;
606 }
607 const double updated_objective_error =
608 exact_objective_value - updated_objective_value;
609 const double updated_objective_absolute_error = fabs(updated_objective_error);
610 const double updated_objective_relative_error =
611 updated_objective_absolute_error / max(1.0, fabs(exact_objective_value));
612 updated_objective_correction += updated_objective_error;
613
614 // Now update the records of previous objective value
615 if (algorithm == SimplexAlgorithm::PRIMAL) {
616 have_previous_exact_primal_objective_value = true;
617 previous_exact_primal_objective_value = exact_objective_value;
618 previous_updated_primal_objective_value = updated_objective_value;
619 updated_primal_objective_correction = updated_objective_correction;
620 } else {
621 have_previous_exact_dual_objective_value = true;
622 previous_exact_dual_objective_value = exact_objective_value;
623 previous_updated_dual_objective_value = updated_objective_value;
624 updated_dual_objective_correction = updated_objective_correction;
625 }
626
627 // Now analyse the error
628 HighsDebugStatus return_status = HighsDebugStatus::OK;
629 std::string error_adjective;
630 int report_level;
631 bool at_least_small_error =
632 updated_objective_relative_error >
633 updated_objective_small_relative_error ||
634 updated_objective_absolute_error > updated_objective_small_absolute_error;
635 if (!at_least_small_error) return return_status;
636 if (updated_objective_relative_error >
637 updated_objective_large_relative_error ||
638 updated_objective_absolute_error >
639 updated_objective_large_absolute_error) {
640 error_adjective = "Large";
641 report_level = ML_ALWAYS;
642 return_status = HighsDebugStatus::LARGE_ERROR;
643 } else if (updated_objective_relative_error >
644 updated_objective_small_relative_error ||
645 updated_objective_absolute_error >
646 updated_objective_small_absolute_error) {
647 error_adjective = "Small";
648 report_level = ML_DETAILED;
649 return_status = HighsDebugStatus::SMALL_ERROR;
650 } else {
651 error_adjective = "OK";
652 report_level = ML_VERBOSE;
653 return_status = HighsDebugStatus::OK;
654 }
655 HighsPrintMessage(
656 highs_model_object.options_.output,
657 highs_model_object.options_.message_level, report_level,
658 "UpdateObjVal: %-9s absolute (%9.4g) or relative (%9.4g) error in "
659 "updated %s objective value"
660 " - objective change - exact (%9.4g) updated (%9.4g) | %s\n",
661 error_adjective.c_str(), updated_objective_error,
662 updated_objective_relative_error, algorithm_name.c_str(),
663 change_in_objective_value, change_in_updated_objective_value,
664 message.c_str());
665 return return_status;
666 }
667
debugUpdatedObjectiveValue(const HighsModelObject & highs_model_object,const SimplexAlgorithm algorithm)668 HighsDebugStatus debugUpdatedObjectiveValue(
669 const HighsModelObject& highs_model_object,
670 const SimplexAlgorithm algorithm) {
671 // Cheap check of updated objective value - assumes that the
672 // objective value computed directly is correct, so only call after
673 // this has been done
674 if (highs_model_object.options_.highs_debug_level == HIGHS_DEBUG_LEVEL_NONE)
675 return HighsDebugStatus::NOT_CHECKED;
676 const HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
677 std::string algorithm_name = "dual";
678 if (algorithm == SimplexAlgorithm::PRIMAL) algorithm_name = "primal";
679 double exact_objective_value;
680 double updated_objective_value;
681 if (algorithm == SimplexAlgorithm::PRIMAL) {
682 assert(highs_model_object.simplex_lp_status_.has_primal_objective_value);
683 exact_objective_value = simplex_info.primal_objective_value;
684 updated_objective_value = simplex_info.updated_primal_objective_value;
685 } else {
686 assert(highs_model_object.simplex_lp_status_.has_dual_objective_value);
687 exact_objective_value = simplex_info.dual_objective_value;
688 updated_objective_value = simplex_info.updated_dual_objective_value;
689 }
690 const double updated_objective_error =
691 exact_objective_value - updated_objective_value;
692 const double updated_objective_absolute_error = fabs(updated_objective_error);
693 const double updated_objective_relative_error =
694 updated_objective_absolute_error / max(1.0, fabs(exact_objective_value));
695
696 // Now analyse the error
697 HighsDebugStatus return_status = HighsDebugStatus::OK;
698 std::string error_adjective;
699 int report_level;
700 if (updated_objective_relative_error >
701 updated_objective_large_relative_error ||
702 updated_objective_absolute_error >
703 updated_objective_large_absolute_error) {
704 error_adjective = "Large";
705 report_level = ML_ALWAYS;
706 return_status = HighsDebugStatus::LARGE_ERROR;
707 } else if (updated_objective_relative_error >
708 updated_objective_small_relative_error ||
709 updated_objective_absolute_error >
710 updated_objective_small_absolute_error) {
711 error_adjective = "Small";
712 report_level = ML_DETAILED;
713 return_status = HighsDebugStatus::SMALL_ERROR;
714 } else {
715 error_adjective = "OK";
716 report_level = ML_VERBOSE;
717 return_status = HighsDebugStatus::OK;
718 }
719 HighsPrintMessage(highs_model_object.options_.output,
720 highs_model_object.options_.message_level, report_level,
721 "UpdateObjVal: %-9s large absolute (%9.4g) or relative "
722 "(%9.4g) error in updated %s objective value\n",
723 error_adjective.c_str(), updated_objective_error,
724 updated_objective_relative_error, algorithm_name.c_str());
725 return return_status;
726 }
727
debugFixedNonbasicMove(const HighsModelObject & highs_model_object)728 HighsDebugStatus debugFixedNonbasicMove(
729 const HighsModelObject& highs_model_object) {
730 // Non-trivially expensive check of nonbasicMove for fixed variables
731 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
732 return HighsDebugStatus::NOT_CHECKED;
733 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
734 const HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
735 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
736 int num_fixed_variable_move_errors = 0;
737 for (int iVar = 0; iVar < simplex_lp.numCol_ + simplex_lp.numRow_; iVar++) {
738 if (!simplex_basis.nonbasicFlag_[iVar]) continue;
739 // Nonbasic column
740 if (simplex_info.workLower_[iVar] == simplex_info.workUpper_[iVar] &&
741 simplex_basis.nonbasicMove_[iVar])
742 num_fixed_variable_move_errors++;
743 }
744 assert(num_fixed_variable_move_errors == 0);
745 if (num_fixed_variable_move_errors) {
746 HighsPrintMessage(highs_model_object.options_.output,
747 highs_model_object.options_.message_level, ML_ALWAYS,
748 "There are %d fixed nonbasicMove errors",
749 num_fixed_variable_move_errors);
750 return HighsDebugStatus::LOGICAL_ERROR;
751 }
752 return HighsDebugStatus::OK;
753 }
754
debugNonbasicMove(const HighsModelObject & highs_model_object)755 HighsDebugStatus debugNonbasicMove(const HighsModelObject& highs_model_object) {
756 // Non-trivially expensive check of NonbasicMove
757 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
758 return HighsDebugStatus::NOT_CHECKED;
759 HighsDebugStatus return_status = HighsDebugStatus::OK;
760 const HighsOptions& options = highs_model_object.options_;
761 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
762 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
763 int num_free_variable_move_errors = 0;
764 int num_lower_bounded_variable_move_errors = 0;
765 int num_upper_bounded_variable_move_errors = 0;
766 int num_boxed_variable_move_errors = 0;
767 int num_fixed_variable_move_errors = 0;
768 const int numTot = simplex_lp.numCol_ + simplex_lp.numRow_;
769 bool right_size = (int)simplex_basis.nonbasicMove_.size() == numTot;
770 // Check consistency of nonbasicMove
771 if (!right_size) {
772 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
773 "nonbasicMove size error");
774 assert(right_size);
775 return_status = HighsDebugStatus::LOGICAL_ERROR;
776 }
777 double lower;
778 double upper;
779
780 for (int iVar = 0; iVar < numTot; iVar++) {
781 if (!simplex_basis.nonbasicFlag_[iVar]) continue;
782 // Nonbasic variable
783 if (iVar < simplex_lp.numCol_) {
784 lower = simplex_lp.colLower_[iVar];
785 upper = simplex_lp.colUpper_[iVar];
786 } else {
787 int iRow = iVar - simplex_lp.numCol_;
788 lower = -simplex_lp.rowUpper_[iRow];
789 upper = -simplex_lp.rowLower_[iRow];
790 }
791
792 if (highs_isInfinity(upper)) {
793 if (highs_isInfinity(-lower)) {
794 // Free
795 if (simplex_basis.nonbasicMove_[iVar]) {
796 num_free_variable_move_errors++;
797 }
798 } else {
799 // Only lower bounded
800 if (simplex_basis.nonbasicMove_[iVar] != NONBASIC_MOVE_UP) {
801 num_lower_bounded_variable_move_errors++;
802 }
803 }
804 } else {
805 if (highs_isInfinity(-lower)) {
806 // Only upper bounded
807 if (simplex_basis.nonbasicMove_[iVar] != NONBASIC_MOVE_DN) {
808 num_upper_bounded_variable_move_errors++;
809 }
810 } else {
811 // Boxed or fixed
812 if (lower != upper) {
813 // Boxed
814 if (!simplex_basis.nonbasicMove_[iVar]) {
815 num_boxed_variable_move_errors++;
816 }
817 } else {
818 // Fixed
819 if (simplex_basis.nonbasicMove_[iVar]) {
820 num_fixed_variable_move_errors++;
821 }
822 }
823 }
824 }
825 }
826 int num_errors =
827 num_free_variable_move_errors + num_lower_bounded_variable_move_errors +
828 num_upper_bounded_variable_move_errors + num_boxed_variable_move_errors +
829 num_fixed_variable_move_errors;
830
831 if (num_errors) {
832 HighsLogMessage(
833 options.logfile, HighsMessageType::ERROR,
834 "There are %d nonbasicMove errors: %d free; %d lower; %d upper; %d "
835 "boxed; %d fixed",
836 num_errors, num_free_variable_move_errors,
837 num_lower_bounded_variable_move_errors,
838 num_upper_bounded_variable_move_errors, num_boxed_variable_move_errors,
839 num_fixed_variable_move_errors);
840 assert(num_errors == 0);
841 return_status = HighsDebugStatus::LOGICAL_ERROR;
842 }
843 return return_status;
844 }
845
debugBasisCondition(const HighsModelObject & highs_model_object,const std::string message)846 HighsDebugStatus debugBasisCondition(const HighsModelObject& highs_model_object,
847 const std::string message) {
848 // Non-trivially expensive assessment of basis condition
849 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
850 return HighsDebugStatus::NOT_CHECKED;
851 double basis_condition = computeBasisCondition(highs_model_object);
852 std::string value_adjective;
853 int report_level;
854 HighsDebugStatus return_status = HighsDebugStatus::OK;
855 if (basis_condition > excessive_basis_condition) {
856 value_adjective = "Excessive";
857 report_level = ML_ALWAYS;
858 return_status = HighsDebugStatus::ERROR;
859 } else if (basis_condition > large_basis_condition) {
860 value_adjective = "Large";
861 report_level = ML_DETAILED;
862 return_status = HighsDebugStatus::WARNING;
863 } else if (basis_condition > fair_basis_condition) {
864 value_adjective = "Fair";
865 report_level = ML_VERBOSE;
866 return_status = HighsDebugStatus::OK;
867 } else {
868 value_adjective = "OK";
869 report_level = ML_VERBOSE;
870 return_status = HighsDebugStatus::OK;
871 }
872 HighsPrintMessage(
873 highs_model_object.options_.output,
874 highs_model_object.options_.message_level, report_level,
875 "BasisCond: %-9s basis condition estimate (%9.4g) - %s\n",
876 value_adjective.c_str(), basis_condition, message.c_str());
877 return return_status;
878 }
879
debugCleanup(HighsModelObject & highs_model_object,const std::vector<double> & original_dual)880 HighsDebugStatus debugCleanup(HighsModelObject& highs_model_object,
881 const std::vector<double>& original_dual) {
882 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
883 return HighsDebugStatus::NOT_CHECKED;
884 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
885 const HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
886 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
887 #ifdef HiGHSDEV
888 HighsSimplexAnalysis& analysis = highs_model_object.simplex_analysis_;
889 #endif
890 // Make sure that the original_dual has been set up
891 assert((int)original_dual.size() == simplex_lp.numCol_ + simplex_lp.numRow_);
892 const std::vector<double>& new_dual = simplex_info.workDual_;
893
894 const double dual_feasibility_tolerance =
895 highs_model_object.scaled_solution_params_.dual_feasibility_tolerance;
896 int num_dual_sign_change = 0;
897 double cleanup_absolute_nonbasic_dual_change_norm = 0;
898 double cleanup_absolute_nonbasic_dual_norm = 0;
899 for (int iVar = 0; iVar < simplex_lp.numCol_ + simplex_lp.numRow_; iVar++) {
900 if (!simplex_basis.nonbasicFlag_[iVar]) continue;
901 cleanup_absolute_nonbasic_dual_norm += std::fabs(new_dual[iVar]);
902 #ifdef HiGHSDEV
903 const double nonbasic_dual_change =
904 std::fabs(new_dual[iVar] - original_dual[iVar]);
905 updateValueDistribution(nonbasic_dual_change,
906 analysis.cleanup_dual_change_distribution);
907 cleanup_absolute_nonbasic_dual_change_norm += nonbasic_dual_change;
908 #endif
909 const double max_dual =
910 std::max(std::fabs(new_dual[iVar]), std::fabs(original_dual[iVar]));
911 if (max_dual > dual_feasibility_tolerance &&
912 new_dual[iVar] * original_dual[iVar] < 0)
913 num_dual_sign_change++;
914 }
915 // Comment on the norm of the nonbasic duals being zero
916 HighsDebugStatus return_status = HighsDebugStatus::OK;
917 if (!cleanup_absolute_nonbasic_dual_norm) {
918 HighsLogMessage(highs_model_object.options_.logfile,
919 HighsMessageType::WARNING,
920 "DualCleanup: dual norm is = %9.4g",
921 cleanup_absolute_nonbasic_dual_norm);
922 return_status = HighsDebugStatus::WARNING;
923 }
924 // Comment on the norm of the change being zero
925 if (!cleanup_absolute_nonbasic_dual_change_norm) {
926 HighsLogMessage(highs_model_object.options_.logfile,
927 HighsMessageType::WARNING,
928 "DualCleanup: dual norm is = %9.4g",
929 cleanup_absolute_nonbasic_dual_change_norm);
930 return_status = HighsDebugStatus::WARNING;
931 }
932 double cleanup_relative_nonbasic_dual_change_norm;
933 if (cleanup_absolute_nonbasic_dual_norm) {
934 cleanup_relative_nonbasic_dual_change_norm =
935 cleanup_absolute_nonbasic_dual_change_norm /
936 cleanup_absolute_nonbasic_dual_norm;
937 } else {
938 cleanup_relative_nonbasic_dual_change_norm = -1;
939 }
940 std::string value_adjective;
941 int report_level;
942 if (cleanup_absolute_nonbasic_dual_change_norm >
943 cleanup_excessive_absolute_nonbasic_dual_change_norm ||
944 cleanup_relative_nonbasic_dual_change_norm >
945 cleanup_excessive_relative_nonbasic_dual_change_norm) {
946 value_adjective = "Excessive";
947 report_level = ML_ALWAYS;
948 return_status = HighsDebugStatus::ERROR;
949 } else if (cleanup_absolute_nonbasic_dual_change_norm >
950 cleanup_large_absolute_nonbasic_dual_change_norm ||
951 cleanup_relative_nonbasic_dual_change_norm >
952 cleanup_large_relative_nonbasic_dual_change_norm) {
953 value_adjective = "Large";
954 report_level = ML_DETAILED;
955 return_status = HighsDebugStatus::WARNING;
956 } else {
957 value_adjective = "OK";
958 report_level = ML_VERBOSE;
959 return_status = HighsDebugStatus::OK;
960 }
961 HighsPrintMessage(
962 highs_model_object.options_.output,
963 highs_model_object.options_.message_level, report_level,
964 "DualCleanup: %-9s absolute (%9.4g) or relative (%9.4g) dual change, "
965 "with %d meaningful sign change(s)\n",
966 value_adjective.c_str(), cleanup_absolute_nonbasic_dual_change_norm,
967 cleanup_relative_nonbasic_dual_change_norm, num_dual_sign_change);
968 return return_status;
969 }
970
debugFreeListNumEntries(const HighsModelObject & highs_model_object,const std::set<int> & freeList)971 HighsDebugStatus debugFreeListNumEntries(
972 const HighsModelObject& highs_model_object, const std::set<int>& freeList) {
973 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
974 return HighsDebugStatus::NOT_CHECKED;
975
976 int freelist_num_entries = 0;
977 if (freeList.size() > 0) {
978 std::set<int>::iterator sit;
979 for (sit = freeList.begin(); sit != freeList.end(); sit++)
980 freelist_num_entries++;
981 }
982
983 const int numTot = highs_model_object.simplex_lp_.numCol_ +
984 highs_model_object.simplex_lp_.numRow_;
985 double pct_freelist_num_entries = (100.0 * freelist_num_entries) / numTot;
986
987 std::string value_adjective;
988 int report_level;
989 HighsDebugStatus return_status = HighsDebugStatus::NOT_CHECKED;
990
991 if (pct_freelist_num_entries > freelist_excessive_pct_num_entries) {
992 value_adjective = "Excessive";
993 report_level = ML_ALWAYS;
994 } else if (pct_freelist_num_entries > freelist_large_pct_num_entries) {
995 value_adjective = "Large";
996 report_level = ML_DETAILED;
997 } else if (pct_freelist_num_entries > freelist_fair_pct_num_entries) {
998 value_adjective = "Fair";
999 report_level = ML_VERBOSE;
1000 } else {
1001 value_adjective = "OK";
1002 if (freelist_num_entries) {
1003 report_level = ML_ALWAYS;
1004 } else {
1005 report_level = ML_VERBOSE;
1006 }
1007 return_status = HighsDebugStatus::OK;
1008 }
1009
1010 HighsPrintMessage(
1011 highs_model_object.options_.output,
1012 highs_model_object.options_.message_level, report_level,
1013 "FreeList : %-9s percentage (%6.2g) of %d variables on free list\n",
1014 value_adjective.c_str(), pct_freelist_num_entries, numTot);
1015
1016 return return_status;
1017 }
1018
debugDualChuzcFail(const HighsOptions & options,const int workCount,const std::vector<std::pair<int,double>> & workData,const double * workDual,const double selectTheta,const double remainTheta)1019 HighsDebugStatus debugDualChuzcFail(
1020 const HighsOptions& options, const int workCount,
1021 const std::vector<std::pair<int, double>>& workData, const double* workDual,
1022 const double selectTheta, const double remainTheta) {
1023 // Non-trivially expensive assessment of basis condition
1024 if (options.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
1025 return HighsDebugStatus::NOT_CHECKED;
1026
1027 HighsPrintMessage(options.output, options.message_level, ML_ALWAYS,
1028 "DualChuzC: No change in loop 2 so return error\n");
1029 double workDataNorm = 0;
1030 double dualNorm = 0;
1031 for (int i = 0; i < workCount; i++) {
1032 int iCol = workData[i].first;
1033 double value = workData[i].second;
1034 workDataNorm += value * value;
1035 value = workDual[iCol];
1036 dualNorm += value * value;
1037 }
1038 workDataNorm += sqrt(workDataNorm);
1039 dualNorm += sqrt(dualNorm);
1040 HighsPrintMessage(
1041 options.output, options.message_level, ML_ALWAYS,
1042 "DualChuzC: workCount = %d; selectTheta=%g; remainTheta=%g\n",
1043 workCount, selectTheta, remainTheta);
1044 HighsPrintMessage(options.output, options.message_level, ML_ALWAYS,
1045 "DualChuzC: workDataNorm = %g; dualNorm = %g\n",
1046 workDataNorm, dualNorm);
1047 return HighsDebugStatus::OK;
1048 }
1049
debugDualChuzcWorkDataAndGroupReport(const HighsModelObject & highs_model_object,const double workDelta,const double workTheta,const std::string message,const int report_workCount,const std::vector<std::pair<int,double>> & report_workData,const std::vector<int> & report_workGroup)1050 void debugDualChuzcWorkDataAndGroupReport(
1051 const HighsModelObject& highs_model_object, const double workDelta,
1052 const double workTheta, const std::string message,
1053 const int report_workCount,
1054 const std::vector<std::pair<int, double>>& report_workData,
1055 const std::vector<int>& report_workGroup) {
1056 const HighsOptions& options = highs_model_object.options_;
1057 const std::vector<int>& workMove =
1058 highs_model_object.simplex_basis_.nonbasicMove_;
1059 const std::vector<double>& workDual =
1060 highs_model_object.simplex_info_.workDual_;
1061 const std::vector<double>& workRange =
1062 highs_model_object.simplex_info_.workRange_;
1063 const double Td =
1064 highs_model_object.scaled_solution_params_.dual_feasibility_tolerance;
1065 double totalChange = initial_total_change;
1066 const double totalDelta = fabs(workDelta);
1067 HighsPrintMessage(
1068 options.output, options.message_level, ML_ALWAYS,
1069 "\n%s: totalDelta = %10.4g\nworkData\n En iCol Dual Value "
1070 " Ratio Change\n",
1071 message.c_str(), totalDelta);
1072 for (int i = 0; i < report_workCount; i++) {
1073 int iCol = report_workData[i].first;
1074 double value = report_workData[i].second;
1075 double dual = workMove[iCol] * workDual[iCol];
1076 totalChange += value * (workRange[iCol]);
1077 HighsPrintMessage(options.output, options.message_level, ML_ALWAYS,
1078 "%4d %4d %10.4g %10.4g %10.4g %10.4g\n", i, iCol, dual,
1079 value, dual / value, totalChange);
1080 }
1081 double selectTheta = workTheta;
1082 HighsPrintMessage(options.output, options.message_level, ML_ALWAYS,
1083 "workGroup\n Ix: selectTheta Entries\n");
1084 for (int group = 0; group < (int)report_workGroup.size() - 1; group++) {
1085 HighsPrintMessage(options.output, options.message_level, ML_ALWAYS,
1086 "%4d: selectTheta = %10.4g ", group, selectTheta);
1087 for (int en = report_workGroup[group]; en < report_workGroup[group + 1];
1088 en++) {
1089 HighsPrintMessage(options.output, options.message_level, ML_ALWAYS,
1090 "%4d ", en);
1091 }
1092 HighsPrintMessage(options.output, options.message_level, ML_ALWAYS, "\n");
1093 int en = report_workGroup[group + 1];
1094 int iCol = report_workData[en].first;
1095 double value = report_workData[en].second;
1096 double dual = workMove[iCol] * workDual[iCol];
1097 selectTheta = (dual + Td) / value;
1098 }
1099 }
1100
debugDualChuzcWorkDataAndGroup(const HighsModelObject & highs_model_object,const double workDelta,const double workTheta,const int workCount,const int alt_workCount,const int breakIndex,const int alt_breakIndex,const std::vector<std::pair<int,double>> & workData,const std::vector<std::pair<int,double>> & sorted_workData,const std::vector<int> & workGroup,const std::vector<int> & alt_workGroup)1101 HighsDebugStatus debugDualChuzcWorkDataAndGroup(
1102 const HighsModelObject& highs_model_object, const double workDelta,
1103 const double workTheta, const int workCount, const int alt_workCount,
1104 const int breakIndex, const int alt_breakIndex,
1105 const std::vector<std::pair<int, double>>& workData,
1106 const std::vector<std::pair<int, double>>& sorted_workData,
1107 const std::vector<int>& workGroup, const std::vector<int>& alt_workGroup) {
1108 // Cheap comparison and possible non-trivially expensive reporting
1109 // of the two sorting methods for BFRT nodes in dual CHUZC
1110 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
1111 return HighsDebugStatus::NOT_CHECKED;
1112 const HighsOptions& options = highs_model_object.options_;
1113 HighsDebugStatus return_status = HighsDebugStatus::OK;
1114 int workPivot = workData[breakIndex].first;
1115 int alt_workPivot = sorted_workData[alt_breakIndex].first;
1116 if (alt_workPivot != workPivot) {
1117 HighsPrintMessage(options.output, options.message_level, ML_ALWAYS,
1118 "Quad workPivot = %d; Heap workPivot = %d\n", workPivot,
1119 alt_workPivot);
1120 return_status = HighsDebugStatus::WARNING;
1121 if (highs_model_object.options_.highs_debug_level <
1122 HIGHS_DEBUG_LEVEL_COSTLY)
1123 return return_status;
1124 debugDualChuzcWorkDataAndGroupReport(highs_model_object, workDelta,
1125 workTheta, "Original", workCount,
1126 workData, workGroup);
1127 debugDualChuzcWorkDataAndGroupReport(
1128 highs_model_object, workDelta, workTheta, "Heap-derived", alt_workCount,
1129 sorted_workData, alt_workGroup);
1130 }
1131 return return_status;
1132 }
1133
debugSimplexBasicSolution(const string message,const HighsModelObject & highs_model_object)1134 HighsDebugStatus debugSimplexBasicSolution(
1135 const string message, const HighsModelObject& highs_model_object) {
1136 // Non-trivially expensive analysis of a simplex basic solution, starting from
1137 // solution_params
1138 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
1139 return HighsDebugStatus::NOT_CHECKED;
1140
1141 if (highsStatusFromHighsModelStatus(
1142 highs_model_object.scaled_model_status_) != HighsStatus::OK)
1143 return HighsDebugStatus::OK;
1144 HighsDebugStatus return_status = HighsDebugStatus::NOT_CHECKED;
1145
1146 const HighsLp& lp = highs_model_object.lp_;
1147 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
1148 const HighsScale& scale = highs_model_object.scale_;
1149 const HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
1150 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
1151
1152 return_status = debugSimplexInfoBasisRightSize(highs_model_object);
1153 if (return_status == HighsDebugStatus::LOGICAL_ERROR) return return_status;
1154
1155 // Determine a HiGHS basis from the simplex basis. Only basic/nonbasic is
1156 // needed
1157 HighsBasis basis;
1158 basis.col_status.resize(lp.numCol_);
1159 basis.row_status.resize(lp.numRow_);
1160 // Now scatter the indices of basic variables
1161 for (int iVar = 0; iVar < lp.numCol_ + lp.numRow_; iVar++) {
1162 if (iVar < lp.numCol_) {
1163 int iCol = iVar;
1164 if (simplex_basis.nonbasicFlag_[iVar] == NONBASIC_FLAG_TRUE) {
1165 basis.col_status[iCol] = HighsBasisStatus::NONBASIC;
1166 } else {
1167 basis.col_status[iCol] = HighsBasisStatus::BASIC;
1168 }
1169 } else {
1170 int iRow = iVar - lp.numCol_;
1171 if (simplex_basis.nonbasicFlag_[iVar] == NONBASIC_FLAG_TRUE) {
1172 basis.row_status[iRow] = HighsBasisStatus::NONBASIC;
1173 } else {
1174 basis.row_status[iRow] = HighsBasisStatus::BASIC;
1175 }
1176 }
1177 }
1178 basis.valid_ = true;
1179 // Possibly scaled model
1180 // Determine a HiGHS solution simplex solution
1181 HighsSolution solution;
1182 solution.col_value.resize(lp.numCol_);
1183 solution.col_dual.resize(lp.numCol_);
1184 solution.row_value.resize(lp.numRow_);
1185 solution.row_dual.resize(lp.numRow_);
1186
1187 for (int iVar = 0; iVar < lp.numCol_ + lp.numRow_; iVar++) {
1188 if (iVar < lp.numCol_) {
1189 int iCol = iVar;
1190 solution.col_value[iCol] = simplex_info.workValue_[iVar];
1191 solution.col_dual[iCol] =
1192 (int)simplex_lp.sense_ * simplex_info.workDual_[iVar];
1193 } else {
1194 int iRow = iVar - lp.numCol_;
1195 solution.row_value[iRow] = -simplex_info.workValue_[iVar];
1196 solution.row_dual[iRow] =
1197 (int)simplex_lp.sense_ * simplex_info.workDual_[iVar];
1198 }
1199 }
1200 // Now insert the basic values
1201 for (int ix = 0; ix < lp.numRow_; ix++) {
1202 int iVar = simplex_basis.basicIndex_[ix];
1203 if (iVar < lp.numCol_) {
1204 solution.col_value[iVar] = simplex_info.baseValue_[ix];
1205 solution.col_dual[iVar] = 0;
1206 } else {
1207 int iRow = iVar - lp.numCol_;
1208 solution.row_value[iRow] = -simplex_info.baseValue_[ix];
1209 solution.row_dual[iRow] = 0;
1210 }
1211 }
1212
1213 const std::string message_scaled = message + " - scaled";
1214 return_status = debugWorseStatus(
1215 debugHighsBasicSolution(message_scaled, highs_model_object.options_,
1216 simplex_lp, basis, solution,
1217 highs_model_object.scaled_solution_params_,
1218 highs_model_object.scaled_model_status_),
1219 return_status);
1220
1221 if (!highs_model_object.scale_.is_scaled_) return return_status;
1222
1223 // Doesn't work if simplex LP has permuted columns
1224 assert(!highs_model_object.simplex_lp_status_.is_permuted);
1225 for (int iCol = 0; iCol < lp.numCol_; iCol++) {
1226 solution.col_value[iCol] *= scale.col_[iCol];
1227 solution.col_dual[iCol] /= (scale.col_[iCol] / scale.cost_);
1228 }
1229 for (int iRow = 0; iRow < simplex_lp.numRow_; iRow++) {
1230 solution.row_value[iRow] /= scale.row_[iRow];
1231 solution.row_dual[iRow] *= (scale.row_[iRow] * scale.cost_);
1232 }
1233 // Cannot assume unscaled solution params or unscaled model status are known
1234 const std::string message_unscaled = message + " - unscaled";
1235 return_status = debugWorseStatus(
1236 debugHighsBasicSolution(message_unscaled, highs_model_object.options_, lp,
1237 basis, solution),
1238 return_status);
1239
1240 // Scaled model
1241 return return_status;
1242 }
1243
debugSimplexHighsSolutionDifferences(const HighsModelObject & highs_model_object)1244 HighsDebugStatus debugSimplexHighsSolutionDifferences(
1245 const HighsModelObject& highs_model_object) {
1246 // Nontrivially expensive check of dimensions and sizes
1247 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
1248 return HighsDebugStatus::NOT_CHECKED;
1249
1250 const HighsOptions& options = highs_model_object.options_;
1251 const HighsSolution& solution = highs_model_object.solution_;
1252 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
1253 const HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
1254 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
1255 const HighsScale& scale = highs_model_object.scale_;
1256
1257 HighsDebugStatus return_status = HighsDebugStatus::NOT_CHECKED;
1258
1259 // Go through the columns, finding the differences in nonbasic column values
1260 // and duals
1261 double max_nonbasic_col_value_difference = 0;
1262 double max_nonbasic_col_dual_difference = 0;
1263 for (int iCol = 0; iCol < simplex_lp.numCol_; iCol++) {
1264 int iVar = iCol;
1265 if (simplex_basis.nonbasicFlag_[iVar] == NONBASIC_FLAG_TRUE) {
1266 // Consider this nonbasic column
1267 double local_col_value = simplex_info.workValue_[iVar] * scale.col_[iCol];
1268 double local_col_dual = (int)simplex_lp.sense_ *
1269 simplex_info.workDual_[iVar] /
1270 (scale.col_[iCol] / scale.cost_);
1271 double value_difference =
1272 fabs(local_col_value - solution.col_value[iCol]);
1273 double dual_difference = fabs(local_col_dual - solution.col_dual[iCol]);
1274 max_nonbasic_col_value_difference =
1275 std::max(value_difference, max_nonbasic_col_value_difference);
1276 max_nonbasic_col_dual_difference =
1277 std::max(dual_difference, max_nonbasic_col_dual_difference);
1278 }
1279 }
1280 // Go through the rows, finding the differences in nonbasic and
1281 // basic row values and duals, as well as differences in basic
1282 // column values and duals
1283 double max_nonbasic_row_value_difference = 0;
1284 double max_nonbasic_row_dual_difference = 0;
1285 double max_basic_col_value_difference = 0;
1286 double max_basic_col_dual_difference = 0;
1287 double max_basic_row_value_difference = 0;
1288 double max_basic_row_dual_difference = 0;
1289
1290 for (int ix = 0; ix < simplex_lp.numRow_; ix++) {
1291 int iRow = ix;
1292 int iVar = simplex_lp.numCol_ + iRow;
1293 if (simplex_basis.nonbasicFlag_[iVar] == NONBASIC_FLAG_TRUE) {
1294 // Consider this nonbasic row
1295 double local_row_value =
1296 -simplex_info.workValue_[iVar] / scale.row_[iRow];
1297 double local_row_dual = (int)simplex_lp.sense_ *
1298 simplex_info.workDual_[iVar] *
1299 (scale.row_[iRow] * scale.cost_);
1300 double value_difference =
1301 fabs(local_row_value - solution.row_value[iRow]);
1302 double dual_difference = fabs(local_row_dual - solution.row_dual[iRow]);
1303 max_nonbasic_row_value_difference =
1304 std::max(value_difference, max_nonbasic_row_value_difference);
1305 max_nonbasic_row_dual_difference =
1306 std::max(dual_difference, max_nonbasic_row_dual_difference);
1307 }
1308 // Consider the basic variable associated with this row index
1309 iVar = simplex_basis.basicIndex_[ix];
1310 if (iVar < simplex_lp.numCol_) {
1311 // Consider this basic column
1312 int iCol = iVar;
1313 double local_col_value = simplex_info.baseValue_[ix] * scale.col_[iCol];
1314 double local_col_dual = 0;
1315 double value_difference =
1316 fabs(local_col_value - solution.col_value[iCol]);
1317 double dual_difference = fabs(local_col_dual - solution.col_dual[iCol]);
1318 max_basic_col_value_difference =
1319 std::max(value_difference, max_basic_col_value_difference);
1320 max_basic_col_dual_difference =
1321 std::max(dual_difference, max_basic_col_dual_difference);
1322 } else {
1323 // Consider this basic row
1324 iRow = iVar - simplex_lp.numCol_;
1325 double local_row_value = -simplex_info.baseValue_[ix] / scale.row_[iRow];
1326 double local_row_dual = 0;
1327 double value_difference =
1328 fabs(local_row_value - solution.row_value[iRow]);
1329 double dual_difference = fabs(local_row_dual - solution.row_dual[iRow]);
1330 max_basic_row_value_difference =
1331 std::max(value_difference, max_basic_row_value_difference);
1332 max_basic_row_dual_difference =
1333 std::max(dual_difference, max_basic_row_dual_difference);
1334 }
1335 }
1336
1337 HighsPrintMessage(options.output, options.message_level, ML_ALWAYS,
1338 "\nHiGHS-simplex solution differences\n");
1339 std::string value_adjective;
1340 int report_level;
1341 return_status = HighsDebugStatus::OK;
1342 if (max_nonbasic_col_value_difference > 0) {
1343 value_adjective = "Excessive";
1344 report_level = ML_ALWAYS;
1345 return_status = debugWorseStatus(HighsDebugStatus::ERROR, return_status);
1346 HighsPrintMessage(
1347 options.output, options.message_level, report_level,
1348 "HighsSimplexD: %-9s Nonbasic column value difference: %9.4g\n",
1349 value_adjective.c_str(), max_nonbasic_col_value_difference);
1350 }
1351 if (max_nonbasic_row_value_difference > 0) {
1352 value_adjective = "Excessive";
1353 report_level = ML_ALWAYS;
1354 return_status = debugWorseStatus(HighsDebugStatus::ERROR, return_status);
1355 HighsPrintMessage(
1356 options.output, options.message_level, report_level,
1357 "HighsSimplexD: %-9s Nonbasic row value difference: %9.4g\n",
1358 value_adjective.c_str(), max_nonbasic_row_value_difference);
1359 }
1360
1361 return_status = debugWorseStatus(
1362 debugAssessSolutionNormDifference(options, "Basic column value",
1363 max_basic_col_value_difference),
1364 return_status);
1365 return_status = debugWorseStatus(
1366 debugAssessSolutionNormDifference(options, "Basic row value",
1367 max_basic_row_value_difference),
1368 return_status);
1369 return_status = debugWorseStatus(
1370 debugAssessSolutionNormDifference(options, "Nonbasic column dual",
1371 max_nonbasic_col_dual_difference),
1372 return_status);
1373 return_status = debugWorseStatus(
1374 debugAssessSolutionNormDifference(options, "Nonbasic row dual",
1375 max_nonbasic_row_dual_difference),
1376 return_status);
1377
1378 if (max_basic_col_dual_difference > 0) {
1379 value_adjective = "Excessive";
1380 report_level = ML_ALWAYS;
1381 return_status = debugWorseStatus(HighsDebugStatus::ERROR, return_status);
1382 HighsPrintMessage(
1383 options.output, options.message_level, report_level,
1384 "HighsSimplexD: %-9s Basic column dual difference: %9.4g\n",
1385 value_adjective.c_str(), max_basic_col_dual_difference);
1386 }
1387 if (max_basic_row_dual_difference > 0) {
1388 value_adjective = "Excessive";
1389 report_level = ML_ALWAYS;
1390 return_status = debugWorseStatus(HighsDebugStatus::ERROR, return_status);
1391 HighsPrintMessage(
1392 options.output, options.message_level, report_level,
1393 "HighsSimplexD: %-9s Basic row dual difference: %9.4g\n",
1394 value_adjective.c_str(), max_basic_row_dual_difference);
1395 }
1396
1397 return return_status;
1398 }
1399
debugAssessSolutionNormDifference(const HighsOptions & options,const std::string type,const double difference)1400 HighsDebugStatus debugAssessSolutionNormDifference(const HighsOptions& options,
1401 const std::string type,
1402 const double difference) {
1403 const double small_difference = 1e-12;
1404 const double large_difference = 1e-8;
1405 const double excessive_difference = 1e-4;
1406 HighsDebugStatus return_status = HighsDebugStatus::OK;
1407 if (difference <= small_difference) return return_status;
1408 std::string value_adjective;
1409 int report_level;
1410
1411 if (difference > excessive_difference) {
1412 value_adjective = "Excessive";
1413 report_level = ML_ALWAYS;
1414 return_status = HighsDebugStatus::ERROR;
1415 } else if (difference > large_difference) {
1416 value_adjective = "Large";
1417 report_level = ML_DETAILED;
1418 return_status = HighsDebugStatus::WARNING;
1419 } else {
1420 value_adjective = "OK";
1421 report_level = ML_VERBOSE;
1422 }
1423 HighsPrintMessage(options.output, options.message_level, report_level,
1424 "HighsSimplexD: %-9s %s difference: %9.4g\n",
1425 value_adjective.c_str(), type.c_str(), difference);
1426 return return_status;
1427 }
1428
debugSimplexBasisCorrect(const HighsModelObject & highs_model_object)1429 HighsDebugStatus debugSimplexBasisCorrect(
1430 const HighsModelObject& highs_model_object) {
1431 // Nontrivially expensive analysis of a Simplex basis, checking
1432 // consistency, and then correctness of nonbasicMove
1433 const HighsOptions& options = highs_model_object.options_;
1434 if (options.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
1435 return HighsDebugStatus::NOT_CHECKED;
1436 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
1437 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
1438 HighsDebugStatus return_status = HighsDebugStatus::OK;
1439 const bool consistent =
1440 debugBasisConsistent(options, simplex_lp, simplex_basis) !=
1441 HighsDebugStatus::LOGICAL_ERROR;
1442 if (!consistent) {
1443 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1444 "Supposed to be a Simplex basis, but not consistent");
1445 assert(consistent);
1446 return_status = HighsDebugStatus::LOGICAL_ERROR;
1447 }
1448 if (options.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
1449 return return_status;
1450 const bool correct_nonbasicMove =
1451 debugNonbasicMove(highs_model_object) != HighsDebugStatus::LOGICAL_ERROR;
1452 if (!correct_nonbasicMove) {
1453 HighsLogMessage(
1454 options.logfile, HighsMessageType::ERROR,
1455 "Supposed to be a Simplex basis, but nonbasicMove is incorrect");
1456 assert(correct_nonbasicMove);
1457 return_status = HighsDebugStatus::LOGICAL_ERROR;
1458 }
1459 return return_status;
1460 }
1461
debugBasisConsistent(const HighsOptions & options,const HighsLp & simplex_lp,const SimplexBasis & simplex_basis)1462 HighsDebugStatus debugBasisConsistent(const HighsOptions& options,
1463 const HighsLp& simplex_lp,
1464 const SimplexBasis& simplex_basis) {
1465 // Cheap analysis of a Simplex basis, checking vector sizes, numbers
1466 // of basic/nonbasic variables and non-repetition of basic variables
1467 if (options.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
1468 return HighsDebugStatus::NOT_CHECKED;
1469 HighsDebugStatus return_status = HighsDebugStatus::OK;
1470 // Check consistency of nonbasicFlag
1471 if (debugNonbasicFlagConsistent(options, simplex_lp, simplex_basis) ==
1472 HighsDebugStatus::LOGICAL_ERROR) {
1473 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1474 "nonbasicFlag inconsistent");
1475 return_status = HighsDebugStatus::LOGICAL_ERROR;
1476 }
1477 const bool right_size =
1478 (int)simplex_basis.basicIndex_.size() == simplex_lp.numRow_;
1479 // Check consistency of basicIndex
1480 if (!right_size) {
1481 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1482 "basicIndex size error");
1483 assert(right_size);
1484 return_status = HighsDebugStatus::LOGICAL_ERROR;
1485 }
1486 // Use localNonbasicFlag so that duplicate entries in basicIndex can
1487 // be spotted
1488 vector<int> localNonbasicFlag = simplex_basis.nonbasicFlag_;
1489 for (int iRow = 0; iRow < simplex_lp.numRow_; iRow++) {
1490 int iCol = simplex_basis.basicIndex_[iRow];
1491 int flag = localNonbasicFlag[iCol];
1492 // Indicate that this column has been found in basicIndex
1493 localNonbasicFlag[iCol] = -1;
1494 if (flag) {
1495 // Nonzero value for localNonbasicFlag entry means that column is either
1496 if (flag == NONBASIC_FLAG_TRUE) {
1497 // Nonbasic...
1498 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1499 "Entry basicIndex_[%d] = %d is not basic", iRow, iCol);
1500 } else {
1501 // .. or is -1 since it has already been found in basicIndex
1502 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1503 "Entry basicIndex_[%d] = %d is already basic", iRow,
1504 iCol);
1505 assert(flag == -1);
1506 }
1507 assert(!flag);
1508 return_status = HighsDebugStatus::LOGICAL_ERROR;
1509 }
1510 }
1511 return return_status;
1512 }
1513
debugNonbasicFlagConsistent(const HighsOptions & options,const HighsLp & simplex_lp,const SimplexBasis & simplex_basis)1514 HighsDebugStatus debugNonbasicFlagConsistent(
1515 const HighsOptions& options, const HighsLp& simplex_lp,
1516 const SimplexBasis& simplex_basis) {
1517 if (options.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
1518 return HighsDebugStatus::NOT_CHECKED;
1519 HighsDebugStatus return_status = HighsDebugStatus::OK;
1520 int numTot = simplex_lp.numCol_ + simplex_lp.numRow_;
1521 const bool right_size = (int)simplex_basis.nonbasicFlag_.size() == numTot;
1522 if (!right_size) {
1523 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1524 "nonbasicFlag size error");
1525 assert(right_size);
1526 return_status = HighsDebugStatus::LOGICAL_ERROR;
1527 }
1528 int num_basic_variables = 0;
1529 for (int var = 0; var < numTot; var++) {
1530 if (simplex_basis.nonbasicFlag_[var] == NONBASIC_FLAG_FALSE) {
1531 num_basic_variables++;
1532 } else {
1533 assert(simplex_basis.nonbasicFlag_[var] == NONBASIC_FLAG_TRUE);
1534 }
1535 }
1536 bool right_num_basic_variables = num_basic_variables == simplex_lp.numRow_;
1537 if (!right_num_basic_variables) {
1538 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1539 "nonbasicFlag has %d, not %d basic variables",
1540 num_basic_variables, simplex_lp.numRow_);
1541 assert(right_num_basic_variables);
1542 return_status = HighsDebugStatus::LOGICAL_ERROR;
1543 }
1544 return return_status;
1545 }
1546
debugOkForSolve(const HighsModelObject & highs_model_object,const int phase)1547 HighsDebugStatus debugOkForSolve(const HighsModelObject& highs_model_object,
1548 const int phase) {
1549 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
1550 return HighsDebugStatus::NOT_CHECKED;
1551 const HighsDebugStatus return_status = HighsDebugStatus::OK;
1552 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
1553 const HighsSimplexLpStatus& simplex_lp_status =
1554 highs_model_object.simplex_lp_status_;
1555 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
1556 const HighsOptions& options = highs_model_object.options_;
1557 bool ok;
1558 // Minimal check - just look at flags. This means we trust them!
1559 ok = simplex_lp_status.has_basis && simplex_lp_status.has_matrix_col_wise &&
1560 simplex_lp_status.has_matrix_row_wise &&
1561 simplex_lp_status.has_factor_arrays &&
1562 simplex_lp_status.has_dual_steepest_edge_weights &&
1563 simplex_lp_status.has_invert;
1564 if (!ok) {
1565 if (!simplex_lp_status.has_basis)
1566 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1567 "Not OK to solve since simplex_lp_status.has_basis = %d",
1568 simplex_lp_status.has_basis);
1569 if (!simplex_lp_status.has_matrix_col_wise)
1570 HighsLogMessage(
1571 options.logfile, HighsMessageType::ERROR,
1572 "Not OK to solve since simplex_lp_status.has_matrix_col_wise "
1573 "= %d",
1574 simplex_lp_status.has_matrix_col_wise);
1575 if (!simplex_lp_status.has_matrix_row_wise)
1576 HighsLogMessage(
1577 options.logfile, HighsMessageType::ERROR,
1578 "Not OK to solve since simplex_lp_status.has_matrix_row_wise "
1579 "= %d",
1580 simplex_lp_status.has_matrix_row_wise);
1581 // if (!simplex_lp_status.has_factor_arrays)
1582 // HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1583 // "Not OK to solve since
1584 // simplex_lp_status.has_factor_arrays = %d",
1585 // simplex_lp_status.has_factor_arrays);
1586 if (!simplex_lp_status.has_dual_steepest_edge_weights)
1587 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1588 "Not OK to solve since "
1589 "simplex_lp_status.has_dual_steepest_edge_weights = %d",
1590 simplex_lp_status.has_dual_steepest_edge_weights);
1591 if (!simplex_lp_status.has_invert)
1592 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1593 "Not OK to solve since simplex_lp_status.has_invert = %d",
1594 simplex_lp_status.has_invert);
1595 }
1596 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_COSTLY)
1597 return return_status;
1598 // Basis and data check
1599 if (debugBasisConsistent(highs_model_object.options_, simplex_lp,
1600 highs_model_object.simplex_basis_) ==
1601 HighsDebugStatus::LOGICAL_ERROR)
1602 return HighsDebugStatus::LOGICAL_ERROR;
1603 if (!debugWorkArraysOk(highs_model_object, phase))
1604 return HighsDebugStatus::LOGICAL_ERROR;
1605 const int numTot = simplex_lp.numCol_ + simplex_lp.numRow_;
1606 for (int var = 0; var < numTot; ++var) {
1607 if (simplex_basis.nonbasicFlag_[var]) {
1608 // Nonbasic variable
1609 if (!debugOneNonbasicMoveVsWorkArraysOk(highs_model_object, var))
1610 return HighsDebugStatus::LOGICAL_ERROR;
1611 }
1612 }
1613 return return_status;
1614 }
1615
1616 // Methods below are not called externally
1617
debugWorkArraysOk(const HighsModelObject & highs_model_object,const int phase)1618 bool debugWorkArraysOk(const HighsModelObject& highs_model_object,
1619 const int phase) {
1620 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
1621 const HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
1622 const HighsOptions& options = highs_model_object.options_;
1623 bool ok = true;
1624 // Only check phase 2 bounds: others will have been set by solve() so can be
1625 // trusted
1626 if (phase == 2) {
1627 for (int col = 0; col < simplex_lp.numCol_; ++col) {
1628 int var = col;
1629 if (!highs_isInfinity(-simplex_info.workLower_[var])) {
1630 ok = simplex_info.workLower_[var] == simplex_lp.colLower_[col];
1631 if (!ok) {
1632 HighsLogMessage(
1633 options.logfile, HighsMessageType::ERROR,
1634 "For col %d, simplex_info.workLower_ should be %g but is %g", col,
1635 simplex_lp.colLower_[col], simplex_info.workLower_[var]);
1636 return ok;
1637 }
1638 }
1639 if (!highs_isInfinity(simplex_info.workUpper_[var])) {
1640 ok = simplex_info.workUpper_[var] == simplex_lp.colUpper_[col];
1641 if (!ok) {
1642 HighsLogMessage(
1643 options.logfile, HighsMessageType::ERROR,
1644 "For col %d, simplex_info.workUpper_ should be %g but is %g", col,
1645 simplex_lp.colUpper_[col], simplex_info.workUpper_[var]);
1646 return ok;
1647 }
1648 }
1649 }
1650 for (int row = 0; row < simplex_lp.numRow_; ++row) {
1651 int var = simplex_lp.numCol_ + row;
1652 if (!highs_isInfinity(-simplex_info.workLower_[var])) {
1653 ok = simplex_info.workLower_[var] == -simplex_lp.rowUpper_[row];
1654 if (!ok) {
1655 HighsLogMessage(
1656 options.logfile, HighsMessageType::ERROR,
1657 "For row %d, simplex_info.workLower_ should be %g but is %g", row,
1658 -simplex_lp.rowUpper_[row], simplex_info.workLower_[var]);
1659 return ok;
1660 }
1661 }
1662 if (!highs_isInfinity(simplex_info.workUpper_[var])) {
1663 ok = simplex_info.workUpper_[var] == -simplex_lp.rowLower_[row];
1664 if (!ok) {
1665 HighsLogMessage(
1666 options.logfile, HighsMessageType::ERROR,
1667 "For row %d, simplex_info.workUpper_ should be %g but is %g", row,
1668 -simplex_lp.rowLower_[row], simplex_info.workUpper_[var]);
1669 return ok;
1670 }
1671 }
1672 }
1673 }
1674 const int numTot = simplex_lp.numCol_ + simplex_lp.numRow_;
1675 for (int var = 0; var < numTot; ++var) {
1676 ok = simplex_info.workRange_[var] ==
1677 (simplex_info.workUpper_[var] - simplex_info.workLower_[var]);
1678 if (!ok) {
1679 HighsLogMessage(
1680 options.logfile, HighsMessageType::ERROR,
1681 "For variable %d, simplex_info.workRange_ should be %g = %g - %g "
1682 "but is %g",
1683 var, simplex_info.workUpper_[var] - simplex_info.workLower_[var],
1684 simplex_info.workUpper_[var], simplex_info.workLower_[var],
1685 simplex_info.workRange_[var]);
1686 return ok;
1687 }
1688 }
1689 // Don't check perturbed costs: these will have been set by solve() so can be
1690 // trusted
1691 if (!simplex_info.costs_perturbed) {
1692 for (int col = 0; col < simplex_lp.numCol_; ++col) {
1693 int var = col;
1694 ok = simplex_info.workCost_[var] ==
1695 (int)simplex_lp.sense_ * simplex_lp.colCost_[col];
1696 if (!ok) {
1697 HighsLogMessage(
1698 options.logfile, HighsMessageType::ERROR,
1699 "For col %d, simplex_info.workLower_ should be %g but is %g", col,
1700 simplex_lp.colLower_[col], simplex_info.workCost_[var]);
1701 return ok;
1702 }
1703 }
1704 for (int row = 0; row < simplex_lp.numRow_; ++row) {
1705 int var = simplex_lp.numCol_ + row;
1706 ok = simplex_info.workCost_[var] == 0.;
1707 if (!ok) {
1708 HighsLogMessage(
1709 options.logfile, HighsMessageType::ERROR,
1710 "For row %d, simplex_info.workCost_ should be zero but is %g", row,
1711 simplex_info.workCost_[var]);
1712 return ok;
1713 }
1714 }
1715 }
1716 // ok must be true if we reach here
1717 assert(ok);
1718 return ok;
1719 }
1720
debugOneNonbasicMoveVsWorkArraysOk(const HighsModelObject & highs_model_object,const int var)1721 bool debugOneNonbasicMoveVsWorkArraysOk(
1722 const HighsModelObject& highs_model_object, const int var) {
1723 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
1724 const HighsSimplexInfo& simplex_info = highs_model_object.simplex_info_;
1725 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
1726 const HighsOptions& options = highs_model_object.options_;
1727 assert(var >= 0);
1728 assert(var < simplex_lp.numCol_ + simplex_lp.numRow_);
1729 // Make sure we're not checking a basic variable
1730 if (!simplex_basis.nonbasicFlag_[var]) return true;
1731 bool ok;
1732 if (!highs_isInfinity(-simplex_info.workLower_[var])) {
1733 if (!highs_isInfinity(simplex_info.workUpper_[var])) {
1734 // Finite lower and upper bounds so nonbasic move depends on whether they
1735 // are equal
1736 if (simplex_info.workLower_[var] == simplex_info.workUpper_[var]) {
1737 // Fixed variable
1738 ok = simplex_basis.nonbasicMove_[var] == NONBASIC_MOVE_ZE;
1739 if (!ok) {
1740 HighsLogMessage(
1741 options.logfile, HighsMessageType::ERROR,
1742 "Fixed variable %d (simplex_lp.numCol_ = %d) [%11g, %11g, "
1743 "%11g] so nonbasic "
1744 "move should be zero but is %d",
1745 var, simplex_lp.numCol_, simplex_info.workLower_[var],
1746 simplex_info.workValue_[var], simplex_info.workUpper_[var],
1747 simplex_basis.nonbasicMove_[var]);
1748 return ok;
1749 }
1750 ok = simplex_info.workValue_[var] == simplex_info.workLower_[var];
1751 if (!ok) {
1752 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1753 "Fixed variable %d (simplex_lp.numCol_ = %d) so "
1754 "simplex_info.work value should be %g but "
1755 "is %g",
1756 var, simplex_lp.numCol_, simplex_info.workLower_[var],
1757 simplex_info.workValue_[var]);
1758 return ok;
1759 }
1760 } else {
1761 // Boxed variable
1762 ok = (simplex_basis.nonbasicMove_[var] == NONBASIC_MOVE_UP) ||
1763 (simplex_basis.nonbasicMove_[var] == NONBASIC_MOVE_DN);
1764 if (!ok) {
1765 HighsLogMessage(
1766 options.logfile, HighsMessageType::ERROR,
1767 "Boxed variable %d (simplex_lp.numCol_ = %d) [%11g, %11g, "
1768 "%11g] range %g so "
1769 "nonbasic move should be up/down but is %d",
1770 var, simplex_lp.numCol_, simplex_info.workLower_[var],
1771 simplex_info.workValue_[var], simplex_info.workUpper_[var],
1772 simplex_info.workUpper_[var] - simplex_info.workLower_[var],
1773 simplex_basis.nonbasicMove_[var]);
1774 return ok;
1775 }
1776 if (simplex_basis.nonbasicMove_[var] == NONBASIC_MOVE_UP) {
1777 ok = simplex_info.workValue_[var] == simplex_info.workLower_[var];
1778 if (!ok) {
1779 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1780 "Boxed variable %d (simplex_lp.numCol_ = %d) with "
1781 "NONBASIC_MOVE_UP so work "
1782 "value should be %g but is %g",
1783 var, simplex_lp.numCol_,
1784 simplex_info.workLower_[var],
1785 simplex_info.workValue_[var]);
1786 return ok;
1787 }
1788 } else {
1789 ok = simplex_info.workValue_[var] == simplex_info.workUpper_[var];
1790 if (!ok) {
1791 HighsLogMessage(options.logfile, HighsMessageType::ERROR,
1792 "Boxed variable %d (simplex_lp.numCol_ = %d) with "
1793 "NONBASIC_MOVE_DN so work "
1794 "value should be %g but is %g",
1795 var, simplex_lp.numCol_,
1796 simplex_info.workUpper_[var],
1797 simplex_info.workValue_[var]);
1798 return ok;
1799 }
1800 }
1801 }
1802 } else {
1803 // Infinite upper bound
1804 ok = simplex_basis.nonbasicMove_[var] == NONBASIC_MOVE_UP;
1805 if (!ok) {
1806 HighsLogMessage(
1807 options.logfile, HighsMessageType::ERROR,
1808 "Finite lower bound and infinite upper bound variable %d "
1809 "(simplex_lp.numCol_ = "
1810 "%d) [%11g, %11g, %11g] so nonbasic move should be up=%2d but is "
1811 "%d",
1812 var, simplex_lp.numCol_, simplex_info.workLower_[var],
1813 simplex_info.workValue_[var], simplex_info.workUpper_[var],
1814 NONBASIC_MOVE_UP, simplex_basis.nonbasicMove_[var]);
1815 return ok;
1816 }
1817 ok = simplex_info.workValue_[var] == simplex_info.workLower_[var];
1818 if (!ok) {
1819 HighsLogMessage(
1820 options.logfile, HighsMessageType::ERROR,
1821 "Finite lower bound and infinite upper bound variable %d "
1822 "(simplex_lp.numCol_ = "
1823 "%d) so work value should be %g but is %g",
1824 var, simplex_lp.numCol_, simplex_info.workLower_[var],
1825 simplex_info.workValue_[var]);
1826 return ok;
1827 }
1828 }
1829 } else {
1830 // Infinite lower bound
1831 if (!highs_isInfinity(simplex_info.workUpper_[var])) {
1832 ok = simplex_basis.nonbasicMove_[var] == NONBASIC_MOVE_DN;
1833 if (!ok) {
1834 HighsLogMessage(
1835 options.logfile, HighsMessageType::ERROR,
1836 "Finite upper bound and infinite lower bound variable %d "
1837 "(simplex_lp.numCol_ = "
1838 "%d) [%11g, %11g, %11g] so nonbasic move should be down but is "
1839 "%d",
1840 var, simplex_lp.numCol_, simplex_info.workLower_[var],
1841 simplex_info.workValue_[var], simplex_info.workUpper_[var],
1842 simplex_basis.nonbasicMove_[var]);
1843 return ok;
1844 }
1845 ok = simplex_info.workValue_[var] == simplex_info.workUpper_[var];
1846 if (!ok) {
1847 HighsLogMessage(
1848 options.logfile, HighsMessageType::ERROR,
1849 "Finite upper bound and infinite lower bound variable %d "
1850 "(simplex_lp.numCol_ = "
1851 "%d) so work value should be %g but is %g",
1852 var, simplex_lp.numCol_, simplex_info.workUpper_[var],
1853 simplex_info.workValue_[var]);
1854 return ok;
1855 }
1856 } else {
1857 // Infinite upper bound
1858 ok = simplex_basis.nonbasicMove_[var] == NONBASIC_MOVE_ZE;
1859 if (!ok) {
1860 HighsLogMessage(
1861 options.logfile, HighsMessageType::ERROR,
1862 "Free variable %d (simplex_lp.numCol_ = %d) [%11g, %11g, %11g] "
1863 "so nonbasic "
1864 "move should be zero but is %d",
1865 var, simplex_lp.numCol_, simplex_info.workLower_[var],
1866 simplex_info.workValue_[var], simplex_info.workUpper_[var],
1867 simplex_basis.nonbasicMove_[var]);
1868 return ok;
1869 }
1870 ok = simplex_info.workValue_[var] == 0.0;
1871 if (!ok) {
1872 HighsLogMessage(
1873 options.logfile, HighsMessageType::ERROR,
1874 "Free variable %d (simplex_lp.numCol_ = %d) so work value should "
1875 "be zero but "
1876 "is %g",
1877 var, simplex_lp.numCol_, simplex_info.workValue_[var]);
1878 return ok;
1879 }
1880 }
1881 }
1882 // ok must be true if we reach here
1883 assert(ok);
1884 return ok;
1885 }
1886
debugAllNonbasicMoveVsWorkArraysOk(const HighsModelObject & highs_model_object)1887 bool debugAllNonbasicMoveVsWorkArraysOk(
1888 const HighsModelObject& highs_model_object) {
1889 const HighsLp& simplex_lp = highs_model_object.simplex_lp_;
1890 // HighsSimplexInfo &simplex_info = highs_model_object.simplex_info_;
1891 const SimplexBasis& simplex_basis = highs_model_object.simplex_basis_;
1892 const HighsOptions& options = highs_model_object.options_;
1893 bool ok;
1894 const int numTot = simplex_lp.numCol_ + simplex_lp.numRow_;
1895 for (int var = 0; var < numTot; ++var) {
1896 HighsLogMessage(
1897 options.logfile, HighsMessageType::ERROR,
1898 "NonbasicMoveVsWorkArrays: var = %2d; simplex_basis.nonbasicFlag_[var] "
1899 "= %2d",
1900 var, simplex_basis.nonbasicFlag_[var]);
1901 if (!simplex_basis.nonbasicFlag_[var]) continue;
1902 ok = debugOneNonbasicMoveVsWorkArraysOk(highs_model_object, var);
1903 if (!ok) {
1904 HighsLogMessage(
1905 options.logfile, HighsMessageType::ERROR,
1906 "Error in NonbasicMoveVsWorkArrays for nonbasic variable %d", var);
1907 assert(ok);
1908 return ok;
1909 }
1910 }
1911 // ok must be true if we reach here
1912 assert(ok);
1913 return ok;
1914 }
1915
debugReportReinvertOnNumericalTrouble(const std::string method_name,const HighsModelObject & highs_model_object,const double numerical_trouble_measure,const double alpha_from_col,const double alpha_from_row,const double numerical_trouble_tolerance,const bool reinvert)1916 void debugReportReinvertOnNumericalTrouble(
1917 const std::string method_name, const HighsModelObject& highs_model_object,
1918 const double numerical_trouble_measure, const double alpha_from_col,
1919 const double alpha_from_row, const double numerical_trouble_tolerance,
1920 const bool reinvert) {
1921 if (highs_model_object.options_.highs_debug_level < HIGHS_DEBUG_LEVEL_CHEAP)
1922 return;
1923 const double abs_alpha_from_col = fabs(alpha_from_col);
1924 const double abs_alpha_from_row = fabs(alpha_from_row);
1925 const double abs_alpha_diff = fabs(abs_alpha_from_col - abs_alpha_from_row);
1926 const int iteration_count = highs_model_object.iteration_counts_.simplex;
1927 const int update_count = highs_model_object.simplex_info_.update_count;
1928 const std::string model_name = highs_model_object.simplex_lp_.model_name_;
1929
1930 const bool numerical_trouble =
1931 numerical_trouble_measure > numerical_trouble_tolerance;
1932 const bool near_numerical_trouble =
1933 10 * numerical_trouble_measure > numerical_trouble_tolerance;
1934
1935 const bool wrong_sign = alpha_from_col * alpha_from_row <= 0;
1936 if (!near_numerical_trouble && !wrong_sign) return;
1937 std::string adjective;
1938 if (numerical_trouble) {
1939 adjective = " exceeds";
1940 } else if (near_numerical_trouble) {
1941 adjective = "almost exceeds";
1942 } else {
1943 adjective = "clearly satisfies";
1944 }
1945 HighsLogMessage(highs_model_object.options_.logfile,
1946 HighsMessageType::WARNING,
1947 "%s (%s) [Iter %d; Update %d] Col: %11.4g; Row: %11.4g; Diff "
1948 "= %11.4g: Measure %11.4g %s %11.4g",
1949 method_name.c_str(), model_name.c_str(), iteration_count,
1950 update_count, abs_alpha_from_col, abs_alpha_from_row,
1951 abs_alpha_diff, numerical_trouble_measure, adjective.c_str(),
1952 numerical_trouble_tolerance);
1953 if (wrong_sign) {
1954 HighsLogMessage(highs_model_object.options_.logfile,
1955 HighsMessageType::WARNING,
1956 " Incompatible signs for Col: %11.4g and Row: %11.4g",
1957 alpha_from_col, alpha_from_row);
1958 }
1959 if ((numerical_trouble || wrong_sign) && !reinvert) {
1960 HighsLogMessage(highs_model_object.options_.logfile,
1961 HighsMessageType::WARNING,
1962 " Numerical trouble or wrong sign and not reinverting");
1963 }
1964 }
1965