1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /*                                                                       */
3 /*    This file is part of the HiGHS linear optimization suite           */
4 /*                                                                       */
5 /*    Written and engineered 2008-2020 at the University of Edinburgh    */
6 /*                                                                       */
7 /*    Available as open-source under the MIT License                     */
8 /*                                                                       */
9 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
10 /**@file simplex/HighsSimplexAnalysis.cpp
11  * @brief
12  * @author Julian Hall, Ivet Galabova, Qi Huangfu and Michael Feldmeier
13  */
14 #include <cmath>
15 //#include <cstdio>
16 #include "HConfig.h"
17 #include "simplex/FactorTimer.h"
18 #include "simplex/HFactor.h"
19 #include "simplex/HighsSimplexAnalysis.h"
20 #include "simplex/SimplexTimer.h"
21 
setup(const HighsLp & lp,const HighsOptions & options,const int simplex_iteration_count_)22 void HighsSimplexAnalysis::setup(const HighsLp& lp, const HighsOptions& options,
23                                  const int simplex_iteration_count_) {
24   // Copy Problem size
25   numRow = lp.numRow_;
26   numCol = lp.numCol_;
27   numTot = numRow + numCol;
28   // Copy tolerances from options
29   allow_dual_steepest_edge_to_devex_switch =
30       options.simplex_dual_edge_weight_strategy ==
31       SIMPLEX_DUAL_EDGE_WEIGHT_STRATEGY_CHOOSE;
32   dual_steepest_edge_weight_log_error_threshold =
33       options.dual_steepest_edge_weight_log_error_threshold;
34   //
35   AnIterIt0 = simplex_iteration_count_;
36   AnIterCostlyDseFq = 0;
37   AnIterNumCostlyDseIt = 0;
38   // Copy messaging parameter from options
39   messaging(options.logfile, options.output, options.message_level);
40   // Initialise the densities
41   col_aq_density = 0;
42   row_ep_density = 0;
43   row_ap_density = 0;
44   row_DSE_density = 0;
45   col_BFRT_density = 0;
46   primal_col_density = 0;
47   // Set the row_dual_density to 1 since it's assumed all costs are at
48   // least perturbed from zero, if not initially nonzero
49   dual_col_density = 1;
50   // Set up the data structures for scatter data
51   tran_stage.resize(NUM_TRAN_STAGE_TYPE);
52   tran_stage[TRAN_STAGE_FTRAN_LOWER].name_ = "FTRAN lower";
53   tran_stage[TRAN_STAGE_FTRAN_UPPER_FT].name_ = "FTRAN upper FT";
54   tran_stage[TRAN_STAGE_FTRAN_UPPER].name_ = "FTRAN upper";
55   tran_stage[TRAN_STAGE_BTRAN_UPPER].name_ = "BTRAN upper";
56   tran_stage[TRAN_STAGE_BTRAN_UPPER_FT].name_ = "BTRAN upper FT";
57   tran_stage[TRAN_STAGE_BTRAN_LOWER].name_ = "BTRAN lower";
58   for (int tran_stage_type = 0; tran_stage_type < NUM_TRAN_STAGE_TYPE;
59        tran_stage_type++) {
60     TranStageAnalysis& stage = tran_stage[tran_stage_type];
61     initialiseScatterData(20, stage.rhs_density_);
62     stage.num_decision_ = 0;
63     stage.num_wrong_original_sparse_decision_ = 0;
64     stage.num_wrong_original_hyper_decision_ = 0;
65     stage.num_wrong_new_sparse_decision_ = 0;
66     stage.num_wrong_new_hyper_decision_ = 0;
67   }
68   original_start_density_tolerance.resize(NUM_TRAN_STAGE_TYPE);
69   new_start_density_tolerance.resize(NUM_TRAN_STAGE_TYPE);
70   historical_density_tolerance.resize(NUM_TRAN_STAGE_TYPE);
71   predicted_density_tolerance.resize(NUM_TRAN_STAGE_TYPE);
72 
73   for (int tran_stage_type = 0; tran_stage_type < NUM_TRAN_STAGE_TYPE;
74        tran_stage_type++) {
75     original_start_density_tolerance[tran_stage_type] = 0.05;
76     new_start_density_tolerance[tran_stage_type] = 0.05;
77   }
78   historical_density_tolerance[TRAN_STAGE_FTRAN_LOWER] = 0.15;
79   historical_density_tolerance[TRAN_STAGE_FTRAN_UPPER] = 0.10;
80   historical_density_tolerance[TRAN_STAGE_BTRAN_UPPER] = 0.10;
81   historical_density_tolerance[TRAN_STAGE_BTRAN_LOWER] = 0.15;
82   predicted_density_tolerance[TRAN_STAGE_FTRAN_LOWER] = 0.10;
83   predicted_density_tolerance[TRAN_STAGE_FTRAN_UPPER] = 0.10;
84   predicted_density_tolerance[TRAN_STAGE_BTRAN_UPPER] = 0.10;
85   predicted_density_tolerance[TRAN_STAGE_BTRAN_LOWER] = 0.10;
86 
87   // Initialise the measures used to analyse accuracy of steepest edge weights
88   //
89   const int dual_edge_weight_strategy =
90       options.simplex_dual_edge_weight_strategy;
91   if (dual_edge_weight_strategy == SIMPLEX_DUAL_EDGE_WEIGHT_STRATEGY_CHOOSE ||
92       dual_edge_weight_strategy ==
93           SIMPLEX_DUAL_EDGE_WEIGHT_STRATEGY_STEEPEST_EDGE ||
94       dual_edge_weight_strategy ==
95           SIMPLEX_DUAL_EDGE_WEIGHT_STRATEGY_STEEPEST_EDGE_UNIT_INITIAL) {
96     // Initialise the measures used to analyse accuracy of steepest edge weights
97     num_dual_steepest_edge_weight_check = 0;
98     num_dual_steepest_edge_weight_reject = 0;
99     num_wrong_low_dual_steepest_edge_weight = 0;
100     num_wrong_high_dual_steepest_edge_weight = 0;
101     average_frequency_low_dual_steepest_edge_weight = 0;
102     average_frequency_high_dual_steepest_edge_weight = 0;
103     average_log_low_dual_steepest_edge_weight_error = 0;
104     average_log_high_dual_steepest_edge_weight_error = 0;
105     max_average_frequency_low_dual_steepest_edge_weight = 0;
106     max_average_frequency_high_dual_steepest_edge_weight = 0;
107     max_sum_average_frequency_extreme_dual_steepest_edge_weight = 0;
108     max_average_log_low_dual_steepest_edge_weight_error = 0;
109     max_average_log_high_dual_steepest_edge_weight_error = 0;
110     max_sum_average_log_extreme_dual_steepest_edge_weight_error = 0;
111   }
112   num_devex_framework = 0;
113 
114   num_iteration_report_since_last_header = -1;
115   num_invert_report_since_last_header = -1;
116 
117   // Set following averages to illegal values so that first average is
118   // set equal to first value
119   average_num_threads = -1;
120   average_fraction_of_possible_minor_iterations_performed = -1;
121   sum_multi_chosen = 0;
122   sum_multi_finished = 0;
123 
124 #ifdef HiGHSDEV
125   AnIterPrevIt = simplex_iteration_count_;
126 
127   SimplexTimer simplex_timer;
128   for (HighsTimerClock& clock : thread_simplex_clocks)
129     simplex_timer.initialiseSimplexClocks(clock);
130 
131   FactorTimer factor_timer;
132   for (HighsTimerClock& clock : thread_factor_clocks)
133     factor_timer.initialiseFactorClocks(clock);
134 
135   AnIterOpRec* AnIter;
136   AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_BTRAN_EP];
137   AnIter->AnIterOpName = "BTRAN e_p";
138   AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_PRICE_AP];
139   AnIter->AnIterOpName = "PRICE a_p";
140   AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_BTRAN_FULL];
141   AnIter->AnIterOpName = "BTRAN Full";
142   AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_PRICE_FULL];
143   AnIter->AnIterOpName = "PRICE Full";
144   AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_FTRAN];
145   AnIter->AnIterOpName = "FTRAN";
146   AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_FTRAN_BFRT];
147   AnIter->AnIterOpName = "FTRAN BFRT";
148   AnIter = &AnIterOp[ANALYSIS_OPERATION_TYPE_FTRAN_DSE];
149   AnIter->AnIterOpName = "FTRAN DSE";
150   for (int k = 0; k < NUM_ANALYSIS_OPERATION_TYPE; k++) {
151     AnIter = &AnIterOp[k];
152     if ((k == ANALYSIS_OPERATION_TYPE_PRICE_AP) ||
153         (k == ANALYSIS_OPERATION_TYPE_PRICE_FULL)) {
154       AnIter->AnIterOpHyperCANCEL = 1.0;
155       AnIter->AnIterOpHyperTRAN = 1.0;
156       AnIter->AnIterOpRsDim = numCol;
157     } else {
158       if ((k == ANALYSIS_OPERATION_TYPE_BTRAN_EP) ||
159           (k == ANALYSIS_OPERATION_TYPE_BTRAN_FULL)) {
160         AnIter->AnIterOpHyperCANCEL = hyperCANCEL;
161         AnIter->AnIterOpHyperTRAN = hyperBTRANU;
162       } else {
163         AnIter->AnIterOpHyperCANCEL = hyperCANCEL;
164         AnIter->AnIterOpHyperTRAN = hyperFTRANL;
165       }
166       AnIter->AnIterOpRsDim = numRow;
167     }
168     AnIter->AnIterOpNumCa = 0;
169     AnIter->AnIterOpNumHyperOp = 0;
170     AnIter->AnIterOpNumHyperRs = 0;
171     AnIter->AnIterOpSumLog10RsDensity = 0;
172     initialiseValueDistribution("", "density ", 1e-8, 1.0, 10.0,
173                                 AnIter->AnIterOp_density);
174   }
175   int last_invert_hint = INVERT_HINT_Count - 1;
176   for (int k = 1; k <= last_invert_hint; k++) AnIterNumInvert[k] = 0;
177   num_col_price = 0;
178   num_row_price = 0;
179   num_row_price_with_switch = 0;
180   int last_dual_edge_weight_mode = (int)DualEdgeWeightMode::STEEPEST_EDGE;
181   for (int k = 0; k <= last_dual_edge_weight_mode; k++) AnIterNumEdWtIt[k] = 0;
182   AnIterTraceNumRec = 0;
183   AnIterTraceIterDl = 1;
184   AnIterTraceRec* lcAnIter = &AnIterTrace[0];
185   lcAnIter->AnIterTraceIter = AnIterIt0;
186   lcAnIter->AnIterTraceTime = timer_->getWallTime();
187   initialiseValueDistribution("Primal step summary", "", 1e-16, 1e16, 10.0,
188                               primal_step_distribution);
189   initialiseValueDistribution("Dual step summary", "", 1e-16, 1e16, 10.0,
190                               dual_step_distribution);
191   initialiseValueDistribution("Simplex pivot summary", "", 1e-8, 1e16, 10.0,
192                               simplex_pivot_distribution);
193   initialiseValueDistribution("Factor pivot threshold summary", "",
194                               min_pivot_threshold, max_pivot_threshold,
195                               pivot_threshold_change_factor,
196                               factor_pivot_threshold_distribution);
197   initialiseValueDistribution("Numerical trouble summary", "", 1e-16, 1.0, 10.0,
198                               numerical_trouble_distribution);
199   initialiseValueDistribution("", "1 ", 1e-16, 1e16, 10.0,
200                               cost_perturbation1_distribution);
201   initialiseValueDistribution("", "2 ", 1e-16, 1e16, 10.0,
202                               cost_perturbation2_distribution);
203   initialiseValueDistribution("FTRAN upper sparse summary - before", "", 1e-8,
204                               1.0, 10.0, before_ftran_upper_sparse_density);
205   initialiseValueDistribution("FTRAN upper sparse summary - after", "", 1e-8,
206                               1.0, 10.0, before_ftran_upper_hyper_density);
207   initialiseValueDistribution("FTRAN upper hyper-sparse summary - before", "",
208                               1e-8, 1.0, 10.0, ftran_upper_sparse_density);
209   initialiseValueDistribution("FTRAN upper hyper-sparse summary - after", "",
210                               1e-8, 1.0, 10.0, ftran_upper_hyper_density);
211   initialiseValueDistribution("Cleanup dual change summary", "", 1e-16, 1e16,
212                               10.0, cleanup_dual_change_distribution);
213   initialiseValueDistribution("Cleanup primal change summary", "", 1e-16, 1e16,
214                               10.0, cleanup_primal_step_distribution);
215   initialiseValueDistribution("Cleanup primal step summary", "", 1e-16, 1e16,
216                               10.0, cleanup_primal_change_distribution);
217   initialiseValueDistribution("Cleanup dual step summary", "", 1e-16, 1e16,
218                               10.0, cleanup_dual_step_distribution);
219 #endif
220 }
221 
messaging(FILE * logfile_,FILE * output_,const int message_level_)222 void HighsSimplexAnalysis::messaging(FILE* logfile_, FILE* output_,
223                                      const int message_level_) {
224   logfile = logfile_;
225   output = output_;
226   message_level = message_level_;
227 }
228 
updateOperationResultDensity(const double local_density,double & density)229 void HighsSimplexAnalysis::updateOperationResultDensity(
230     const double local_density, double& density) {
231   density = (1 - running_average_multiplier) * density +
232             running_average_multiplier * local_density;
233 }
234 
iterationReport()235 void HighsSimplexAnalysis::iterationReport() {
236   if (!(iteration_report_message_level & message_level)) return;
237   const bool header = (num_iteration_report_since_last_header < 0) ||
238                       (num_iteration_report_since_last_header > 49);
239   if (header) {
240     iterationReport(header);
241     num_iteration_report_since_last_header = 0;
242   }
243   iterationReport(false);
244 }
245 
invertReport()246 void HighsSimplexAnalysis::invertReport() {
247   if (!(invert_report_message_level & message_level)) return;
248   const bool header = (num_invert_report_since_last_header < 0) ||
249                       (num_invert_report_since_last_header > 49) ||
250                       (num_iteration_report_since_last_header >= 0);
251   if (header) {
252     invertReport(header);
253     num_invert_report_since_last_header = 0;
254   }
255   invertReport(false);
256   // Force an iteration report header if this is an INVERT report without an
257   // invert_hint
258   if (!invert_hint) num_iteration_report_since_last_header = -1;
259 }
260 
invertReport(const bool header)261 void HighsSimplexAnalysis::invertReport(const bool header) {
262   if (!(invert_report_message_level & message_level)) return;
263   reportAlgorithmPhaseIterationObjective(header, invert_report_message_level);
264 #ifdef HiGHSDEV
265   if (simplex_strategy == SIMPLEX_STRATEGY_DUAL_MULTI) {
266     // Report on threads and PAMI
267     reportThreads(header, invert_report_message_level);
268     reportMulti(header, invert_report_message_level);
269   }
270   reportDensity(header, invert_report_message_level);
271   reportInvert(header, invert_report_message_level);
272   //  reportCondition(header, invert_report_message_level);
273 #endif
274   reportInfeasibility(header, invert_report_message_level);
275   HighsPrintMessage(output, message_level, invert_report_message_level, "\n");
276   if (!header) num_invert_report_since_last_header++;
277 }
278 
dualSteepestEdgeWeightError(const double computed_edge_weight,const double updated_edge_weight)279 void HighsSimplexAnalysis::dualSteepestEdgeWeightError(
280     const double computed_edge_weight, const double updated_edge_weight) {
281   const bool accept_weight =
282       updated_edge_weight >= accept_weight_threshold * computed_edge_weight;
283   int low_weight_error = 0;
284   int high_weight_error = 0;
285   double weight_error;
286 #ifdef HiGHSDEV
287   string error_type = "  OK";
288 #endif
289   num_dual_steepest_edge_weight_check++;
290   if (!accept_weight) num_dual_steepest_edge_weight_reject++;
291   if (updated_edge_weight < computed_edge_weight) {
292     // Updated weight is low
293     weight_error = computed_edge_weight / updated_edge_weight;
294     if (weight_error > weight_error_threshold) {
295       low_weight_error = 1;
296 #ifdef HiGHSDEV
297       error_type = " Low";
298 #endif
299     }
300     average_log_low_dual_steepest_edge_weight_error =
301         0.99 * average_log_low_dual_steepest_edge_weight_error +
302         0.01 * log(weight_error);
303   } else {
304     // Updated weight is correct or high
305     weight_error = updated_edge_weight / computed_edge_weight;
306     if (weight_error > weight_error_threshold) {
307       high_weight_error = 1;
308 #ifdef HiGHSDEV
309       error_type = "High";
310 #endif
311     }
312     average_log_high_dual_steepest_edge_weight_error =
313         0.99 * average_log_high_dual_steepest_edge_weight_error +
314         0.01 * log(weight_error);
315   }
316   average_frequency_low_dual_steepest_edge_weight =
317       0.99 * average_frequency_low_dual_steepest_edge_weight +
318       0.01 * low_weight_error;
319   average_frequency_high_dual_steepest_edge_weight =
320       0.99 * average_frequency_high_dual_steepest_edge_weight +
321       0.01 * high_weight_error;
322   max_average_frequency_low_dual_steepest_edge_weight =
323       max(max_average_frequency_low_dual_steepest_edge_weight,
324           average_frequency_low_dual_steepest_edge_weight);
325   max_average_frequency_high_dual_steepest_edge_weight =
326       max(max_average_frequency_high_dual_steepest_edge_weight,
327           average_frequency_high_dual_steepest_edge_weight);
328   max_sum_average_frequency_extreme_dual_steepest_edge_weight =
329       max(max_sum_average_frequency_extreme_dual_steepest_edge_weight,
330           average_frequency_low_dual_steepest_edge_weight +
331               average_frequency_high_dual_steepest_edge_weight);
332   max_average_log_low_dual_steepest_edge_weight_error =
333       max(max_average_log_low_dual_steepest_edge_weight_error,
334           average_log_low_dual_steepest_edge_weight_error);
335   max_average_log_high_dual_steepest_edge_weight_error =
336       max(max_average_log_high_dual_steepest_edge_weight_error,
337           average_log_high_dual_steepest_edge_weight_error);
338   max_sum_average_log_extreme_dual_steepest_edge_weight_error =
339       max(max_sum_average_log_extreme_dual_steepest_edge_weight_error,
340           average_log_low_dual_steepest_edge_weight_error +
341               average_log_high_dual_steepest_edge_weight_error);
342 #ifdef HiGHSDEV
343   const bool report_weight_error = false;
344   if (report_weight_error && weight_error > 0.5 * weight_error_threshold) {
345     printf(
346         "DSE Wt Ck |%8d| OK = %1d (%4d / %6d) (c %10.4g, u %10.4g, er %10.4g - "
347         "%s): Low (Fq %10.4g, Er %10.4g); High (Fq%10.4g, Er%10.4g) | %10.4g "
348         "%10.4g %10.4g %10.4g %10.4g %10.4g\n",
349         simplex_iteration_count, accept_weight,
350         num_dual_steepest_edge_weight_check,
351         num_dual_steepest_edge_weight_reject, computed_edge_weight,
352         updated_edge_weight, weight_error, error_type.c_str(),
353         average_frequency_low_dual_steepest_edge_weight,
354         average_log_low_dual_steepest_edge_weight_error,
355         average_frequency_high_dual_steepest_edge_weight,
356         average_log_high_dual_steepest_edge_weight_error,
357         max_average_frequency_low_dual_steepest_edge_weight,
358         max_average_frequency_high_dual_steepest_edge_weight,
359         max_sum_average_frequency_extreme_dual_steepest_edge_weight,
360         max_average_log_low_dual_steepest_edge_weight_error,
361         max_average_log_high_dual_steepest_edge_weight_error,
362         max_sum_average_log_extreme_dual_steepest_edge_weight_error);
363   }
364 #endif
365 }
366 
switchToDevex()367 bool HighsSimplexAnalysis::switchToDevex() {
368   bool switch_to_devex = false;
369   // Firstly consider switching on the basis of NLA cost
370   double AnIterCostlyDseMeasureDen;
371   AnIterCostlyDseMeasureDen =
372       max(max(row_ep_density, col_aq_density), row_ap_density);
373   if (AnIterCostlyDseMeasureDen > 0) {
374     AnIterCostlyDseMeasure = row_DSE_density / AnIterCostlyDseMeasureDen;
375     AnIterCostlyDseMeasure = AnIterCostlyDseMeasure * AnIterCostlyDseMeasure;
376   } else {
377     AnIterCostlyDseMeasure = 0;
378   }
379   bool CostlyDseIt = AnIterCostlyDseMeasure > AnIterCostlyDseMeasureLimit &&
380                      row_DSE_density > AnIterCostlyDseMnDensity;
381   AnIterCostlyDseFq = (1 - running_average_multiplier) * AnIterCostlyDseFq;
382   if (CostlyDseIt) {
383     AnIterNumCostlyDseIt++;
384     AnIterCostlyDseFq += running_average_multiplier * 1.0;
385     int lcNumIter = simplex_iteration_count - AnIterIt0;
386     // Switch to Devex if at least 5% of the (at least) 0.1NumTot iterations
387     // have been costly
388     switch_to_devex =
389         allow_dual_steepest_edge_to_devex_switch &&
390         (AnIterNumCostlyDseIt > lcNumIter * AnIterFracNumCostlyDseItbfSw) &&
391         (lcNumIter > AnIterFracNumTot_ItBfSw * numTot);
392 #ifdef HiGHSDEV
393     if (switch_to_devex) {
394       HighsLogMessage(
395           logfile, HighsMessageType::INFO,
396           "Switch from DSE to Devex after %d costly DSE iterations of %d: "
397           "C_Aq_Density = %11.4g; R_Ep_Density = %11.4g; R_Ap_Density = "
398           "%11.4g; DSE_Density = %11.4g",
399           AnIterNumCostlyDseIt, lcNumIter, col_aq_density, row_ep_density,
400           row_ap_density, row_DSE_density);
401     }
402 #endif
403   }
404   if (!switch_to_devex) {
405     // Secondly consider switching on the basis of weight accuracy
406     double dse_weight_error_measure =
407         average_log_low_dual_steepest_edge_weight_error +
408         average_log_high_dual_steepest_edge_weight_error;
409     double dse_weight_error_threshold =
410         dual_steepest_edge_weight_log_error_threshold;
411     switch_to_devex = allow_dual_steepest_edge_to_devex_switch &&
412                       dse_weight_error_measure > dse_weight_error_threshold;
413 #ifdef HiGHSDEV
414     if (switch_to_devex) {
415       HighsLogMessage(logfile, HighsMessageType::INFO,
416                       "Switch from DSE to Devex with log error measure of %g > "
417                       "%g = threshold",
418                       dse_weight_error_measure, dse_weight_error_threshold);
419     }
420 #endif
421   }
422   return switch_to_devex;
423 }
424 
predictEndDensity(const int tran_stage_type,const double start_density,double & end_density)425 bool HighsSimplexAnalysis::predictEndDensity(const int tran_stage_type,
426                                              const double start_density,
427                                              double& end_density) {
428   return predictFromScatterData(tran_stage[tran_stage_type].rhs_density_,
429                                 start_density, end_density);
430 }
431 
afterTranStage(const int tran_stage_type,const double start_density,const double end_density,const double historical_density,const double predicted_end_density,const bool use_solve_sparse_original_HFactor_logic,const bool use_solve_sparse_new_HFactor_logic)432 void HighsSimplexAnalysis::afterTranStage(
433     const int tran_stage_type, const double start_density,
434     const double end_density, const double historical_density,
435     const double predicted_end_density,
436     const bool use_solve_sparse_original_HFactor_logic,
437     const bool use_solve_sparse_new_HFactor_logic) {
438   TranStageAnalysis& stage = tran_stage[tran_stage_type];
439   const double rp = false;
440   if (predicted_end_density > 0) {
441     stage.num_decision_++;
442     if (end_density <= max_hyper_density) {
443       // Should have done hyper-sparse TRAN
444       if (use_solve_sparse_original_HFactor_logic) {
445         // Original logic makes wrong decision to use sparse TRAN
446         if (rp) {
447           printf("Original: Wrong sparse: ");
448           const double start_density_tolerance =
449               original_start_density_tolerance[tran_stage_type];
450           const double this_historical_density_tolerance =
451               historical_density_tolerance[tran_stage_type];
452           if (start_density > start_density_tolerance) {
453             printf("(start = %10.4g >  %4.2f)  or ", start_density,
454                    start_density_tolerance);
455           } else {
456             printf(" start = %10.4g              ", start_density);
457           }
458           if (historical_density > this_historical_density_tolerance) {
459             printf("(historical = %10.4g  > %4.2f); ", historical_density,
460                    this_historical_density_tolerance);
461           } else {
462             printf(" historical = %10.4g           ", historical_density);
463           }
464           printf("end = %10.4g", end_density);
465           if (end_density < 0.1 * historical_density) printf(" !! OG");
466           printf("\n");
467         }
468         stage.num_wrong_original_sparse_decision_++;
469       }
470       if (use_solve_sparse_new_HFactor_logic) {
471         // New logic makes wrong decision to use sparse TRAN
472         if (rp) {
473           printf("New     : Wrong sparse: ");
474           const double start_density_tolerance =
475               original_start_density_tolerance[tran_stage_type];
476           const double end_density_tolerance =
477               predicted_density_tolerance[tran_stage_type];
478           if (start_density > start_density_tolerance) {
479             printf("(start = %10.4g >  %4.2f)  or ", start_density,
480                    start_density_tolerance);
481           } else {
482             printf(" start = %10.4g                       ", start_density);
483           }
484           if (predicted_end_density > end_density_tolerance) {
485             printf("( predicted = %10.4g  > %4.2f); ", predicted_end_density,
486                    end_density_tolerance);
487           } else {
488             printf("  predicted = %10.4g           ", predicted_end_density);
489           }
490           printf("end = %10.4g", end_density);
491           if (end_density < 0.1 * predicted_end_density) printf(" !! NW");
492           printf("\n");
493         }
494         stage.num_wrong_new_sparse_decision_++;
495       }
496     } else {
497       // Should have done sparse TRAN
498       if (!use_solve_sparse_original_HFactor_logic) {
499         // Original logic makes wrong decision to use hyper TRAN
500         if (rp) {
501           printf(
502               "Original: Wrong  hyper: (start = %10.4g <= %4.2f) and "
503               "(historical = %10.4g <= %4.2f); end = %10.4g",
504               start_density, original_start_density_tolerance[tran_stage_type],
505               historical_density, historical_density_tolerance[tran_stage_type],
506               end_density);
507           if (end_density > 10.0 * historical_density) printf(" !! OG");
508           printf("\n");
509         }
510         stage.num_wrong_original_hyper_decision_++;
511       }
512       if (!use_solve_sparse_new_HFactor_logic) {
513         // New logic makes wrong decision to use hyper TRAN
514         if (rp) {
515           printf(
516               "New     : Wrong  hyper: (start = %10.4g <= %4.2f) and ( "
517               "predicted = %10.4g <= %4.2f); end = %10.4g",
518               start_density, new_start_density_tolerance[tran_stage_type],
519               predicted_end_density,
520               predicted_density_tolerance[tran_stage_type], end_density);
521           if (end_density > 10.0 * predicted_end_density) printf(" !! NW");
522           printf("\n");
523         }
524         stage.num_wrong_new_hyper_decision_++;
525       }
526     }
527   }
528   updateScatterData(start_density, end_density, stage.rhs_density_);
529   regressScatterData(stage.rhs_density_);
530 }
531 
summaryReportFactor()532 void HighsSimplexAnalysis::summaryReportFactor() {
533   for (int tran_stage_type = 0; tran_stage_type < NUM_TRAN_STAGE_TYPE;
534        tran_stage_type++) {
535     TranStageAnalysis& stage = tran_stage[tran_stage_type];
536     //    printScatterData(stage.name_, stage.rhs_density_);
537     printScatterDataRegressionComparison(stage.name_, stage.rhs_density_);
538     if (!stage.num_decision_) return;
539     printf("Of %10d Sps/Hyper decisions made using regression:\n",
540            stage.num_decision_);
541     printf(
542         "   %10d wrong sparseTRAN; %10d wrong hyperTRAN: using original "
543         "logic\n",
544         stage.num_wrong_original_sparse_decision_,
545         stage.num_wrong_original_hyper_decision_);
546     printf(
547         "   %10d wrong sparseTRAN; %10d wrong hyperTRAN: using new      "
548         "logic\n",
549         stage.num_wrong_new_sparse_decision_,
550         stage.num_wrong_new_hyper_decision_);
551   }
552 }
553 
simplexTimerStart(const int simplex_clock,const int thread_id)554 void HighsSimplexAnalysis::simplexTimerStart(const int simplex_clock,
555                                              const int thread_id) {
556 #ifdef HiGHSDEV
557   thread_simplex_clocks[thread_id].timer_.start(
558       thread_simplex_clocks[thread_id].clock_[simplex_clock]);
559 #endif
560 }
561 
simplexTimerStop(const int simplex_clock,const int thread_id)562 void HighsSimplexAnalysis::simplexTimerStop(const int simplex_clock,
563                                             const int thread_id) {
564 #ifdef HiGHSDEV
565   thread_simplex_clocks[thread_id].timer_.stop(
566       thread_simplex_clocks[thread_id].clock_[simplex_clock]);
567 #endif
568 }
569 
simplexTimerRunning(const int simplex_clock,const int thread_id)570 bool HighsSimplexAnalysis::simplexTimerRunning(const int simplex_clock,
571                                                const int thread_id) {
572   bool argument = false;
573 #ifdef HiGHSDEV
574   argument =
575       thread_simplex_clocks[thread_id]
576           .timer_
577           .clock_start[thread_simplex_clocks[thread_id].clock_[simplex_clock]] <
578       0;
579 #endif
580   return argument;
581 }
582 
simplexTimerNumCall(const int simplex_clock,const int thread_id)583 int HighsSimplexAnalysis::simplexTimerNumCall(const int simplex_clock,
584                                               const int thread_id) {
585   int argument = 0;
586 #ifdef HiGHSDEV
587   argument = thread_simplex_clocks[thread_id].timer_.clock_num_call
588                  [thread_simplex_clocks[thread_id].clock_[simplex_clock]];
589 #endif
590   return argument;
591 }
592 
simplexTimerRead(const int simplex_clock,const int thread_id)593 double HighsSimplexAnalysis::simplexTimerRead(const int simplex_clock,
594                                               const int thread_id) {
595   double argument = 0;
596 #ifdef HiGHSDEV
597   argument = thread_simplex_clocks[thread_id].timer_.read(
598       thread_simplex_clocks[thread_id].clock_[simplex_clock]);
599 #endif
600   return argument;
601 }
602 
getThreadFactorTimerClockPointer()603 HighsTimerClock* HighsSimplexAnalysis::getThreadFactorTimerClockPointer() {
604   HighsTimerClock* factor_timer_clock_pointer = NULL;
605 #ifdef HiGHSDEV
606   int thread_id = 0;
607 #ifdef OPENMP
608   thread_id = omp_get_thread_num();
609 #endif
610   factor_timer_clock_pointer = &thread_factor_clocks[thread_id];
611 #endif
612   return factor_timer_clock_pointer;
613 }
614 
615 #ifdef HiGHSDEV
reportFactorTimer()616 void HighsSimplexAnalysis::reportFactorTimer() {
617   FactorTimer factor_timer;
618   int omp_max_threads = 1;
619 #ifdef OPENMP
620   omp_max_threads = omp_get_max_threads();
621 #endif
622   for (int i = 0; i < omp_max_threads; i++) {
623     //  for (HighsTimerClock clock : thread_factor_clocks) {
624     printf("reportFactorTimer: HFactor clocks for OMP thread %d / %d\n", i,
625            omp_max_threads - 1);
626     factor_timer.reportFactorClock(thread_factor_clocks[i]);
627   }
628   if (omp_max_threads > 1) {
629     HighsTimer& timer = thread_factor_clocks[0].timer_;
630     HighsTimerClock all_factor_clocks(timer);
631     vector<int>& clock = all_factor_clocks.clock_;
632     factor_timer.initialiseFactorClocks(all_factor_clocks);
633     for (int i = 0; i < omp_max_threads; i++) {
634       vector<int>& thread_clock = thread_factor_clocks[i].clock_;
635       for (int clock_id = 0; clock_id < FactorNumClock; clock_id++) {
636         int all_factor_iClock = clock[clock_id];
637         int thread_factor_iClock = thread_clock[clock_id];
638         timer.clock_num_call[all_factor_iClock] +=
639             timer.clock_num_call[thread_factor_iClock];
640         timer.clock_time[all_factor_iClock] +=
641             timer.clock_time[thread_factor_iClock];
642       }
643     }
644     printf("reportFactorTimer: HFactor clocks for all %d threads\n",
645            omp_max_threads);
646     factor_timer.reportFactorClock(all_factor_clocks);
647   }
648 }
649 
iterationRecord()650 void HighsSimplexAnalysis::iterationRecord() {
651   int AnIterCuIt = simplex_iteration_count;
652   if (invert_hint > 0) AnIterNumInvert[invert_hint]++;
653   if (AnIterCuIt > AnIterPrevIt)
654     AnIterNumEdWtIt[(int)edge_weight_mode] += (AnIterCuIt - AnIterPrevIt);
655 
656   AnIterTraceRec& lcAnIter = AnIterTrace[AnIterTraceNumRec];
657   //  if (simplex_iteration_count ==
658   //  AnIterTraceIterRec[AnIterTraceNumRec]+AnIterTraceIterDl) {
659   if (simplex_iteration_count == lcAnIter.AnIterTraceIter + AnIterTraceIterDl) {
660     if (AnIterTraceNumRec == AN_ITER_TRACE_MX_NUM_REC) {
661       for (int rec = 1; rec <= AN_ITER_TRACE_MX_NUM_REC / 2; rec++)
662         AnIterTrace[rec] = AnIterTrace[2 * rec];
663       AnIterTraceNumRec = AnIterTraceNumRec / 2;
664       AnIterTraceIterDl = AnIterTraceIterDl * 2;
665     } else {
666       AnIterTraceNumRec++;
667       AnIterTraceRec& lcAnIter = AnIterTrace[AnIterTraceNumRec];
668       lcAnIter.AnIterTraceIter = simplex_iteration_count;
669       lcAnIter.AnIterTraceTime = timer_->getWallTime();
670       if (average_fraction_of_possible_minor_iterations_performed > 0) {
671         lcAnIter.AnIterTraceMulti =
672             average_fraction_of_possible_minor_iterations_performed;
673       } else {
674         lcAnIter.AnIterTraceMulti = 0;
675       }
676       lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN] =
677           col_aq_density;
678       lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_BTRAN_EP] =
679           row_ep_density;
680       lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_PRICE_AP] =
681           row_ap_density;
682       lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_BFRT] =
683           col_aq_density;
684       if (edge_weight_mode == DualEdgeWeightMode::STEEPEST_EDGE) {
685         lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE] =
686             row_DSE_density;
687         lcAnIter.AnIterTraceCostlyDse = AnIterCostlyDseMeasure;
688       } else {
689         lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE] = 0;
690         lcAnIter.AnIterTraceCostlyDse = 0;
691       }
692       lcAnIter.AnIterTrace_dual_edge_weight_mode = (int)edge_weight_mode;
693     }
694   }
695   AnIterPrevIt = AnIterCuIt;
696   updateValueDistribution(primal_step, cleanup_primal_step_distribution);
697   updateValueDistribution(dual_step, cleanup_dual_step_distribution);
698   updateValueDistribution(primal_step, primal_step_distribution);
699   updateValueDistribution(dual_step, dual_step_distribution);
700   updateValueDistribution(pivot_value_from_column, simplex_pivot_distribution);
701   updateValueDistribution(factor_pivot_threshold,
702                           factor_pivot_threshold_distribution);
703   // Only update the distribution of legal values for
704   // numerical_trouble. Illegal values are set in PAMI since it's not
705   // known in minor iterations
706   if (numerical_trouble >= 0)
707     updateValueDistribution(numerical_trouble, numerical_trouble_distribution);
708 }
709 
iterationRecordMajor()710 void HighsSimplexAnalysis::iterationRecordMajor() {
711   sum_multi_chosen += multi_chosen;
712   sum_multi_finished += multi_finished;
713   assert(multi_chosen > 0);
714   const double fraction_of_possible_minor_iterations_performed =
715       1.0 * multi_finished / multi_chosen;
716   if (average_fraction_of_possible_minor_iterations_performed < 0) {
717     average_fraction_of_possible_minor_iterations_performed =
718         fraction_of_possible_minor_iterations_performed;
719   } else {
720     average_fraction_of_possible_minor_iterations_performed =
721         running_average_multiplier *
722             fraction_of_possible_minor_iterations_performed +
723         (1 - running_average_multiplier) *
724             average_fraction_of_possible_minor_iterations_performed;
725   }
726   if (average_num_threads < 0) {
727     average_num_threads = num_threads;
728   } else {
729     average_num_threads =
730         running_average_multiplier * num_threads +
731         (1 - running_average_multiplier) * average_num_threads;
732   }
733 }
734 
operationRecordBefore(const int operation_type,const HVector & vector,const double historical_density)735 void HighsSimplexAnalysis::operationRecordBefore(
736     const int operation_type, const HVector& vector,
737     const double historical_density) {
738   operationRecordBefore(operation_type, vector.count, historical_density);
739 }
740 
operationRecordBefore(const int operation_type,const int current_count,const double historical_density)741 void HighsSimplexAnalysis::operationRecordBefore(
742     const int operation_type, const int current_count,
743     const double historical_density) {
744   double current_density = 1.0 * current_count / numRow;
745   AnIterOpRec& AnIter = AnIterOp[operation_type];
746   AnIter.AnIterOpNumCa++;
747   if (current_density <= AnIter.AnIterOpHyperCANCEL &&
748       historical_density <= AnIter.AnIterOpHyperTRAN)
749     AnIter.AnIterOpNumHyperOp++;
750 }
751 
operationRecordAfter(const int operation_type,const HVector & vector)752 void HighsSimplexAnalysis::operationRecordAfter(const int operation_type,
753                                                 const HVector& vector) {
754   operationRecordAfter(operation_type, vector.count);
755 }
756 
operationRecordAfter(const int operation_type,const int result_count)757 void HighsSimplexAnalysis::operationRecordAfter(const int operation_type,
758                                                 const int result_count) {
759   AnIterOpRec& AnIter = AnIterOp[operation_type];
760   const double result_density = 1.0 * result_count / AnIter.AnIterOpRsDim;
761   if (result_density <= hyperRESULT) AnIter.AnIterOpNumHyperRs++;
762   if (result_density > 0) {
763     AnIter.AnIterOpSumLog10RsDensity += log(result_density) / log(10.0);
764   } else {
765     /*
766     // TODO Investigate these zero norms
767     double vectorNorm = 0;
768 
769     for (int index = 0; index < AnIter.AnIterOpRsDim; index++) {
770       double vectorValue = vector.array[index];
771       vectorNorm += vectorValue * vectorValue;
772     }
773     vectorNorm = sqrt(vectorNorm);
774     printf("Strange: operation %s has result density = %g: ||vector|| = %g\n",
775     AnIter.AnIterOpName.c_str(), result_density, vectorNorm);
776     */
777   }
778   updateValueDistribution(result_density, AnIter.AnIterOp_density);
779 }
780 
summaryReport()781 void HighsSimplexAnalysis::summaryReport() {
782   int AnIterNumIter = simplex_iteration_count - AnIterIt0;
783   if (AnIterNumIter <= 0) return;
784   printf("\nAnalysis of %d iterations (%d to %d)\n", AnIterNumIter,
785          AnIterIt0 + 1, simplex_iteration_count);
786   if (AnIterNumIter <= 0) return;
787   int lc_EdWtNumIter;
788   lc_EdWtNumIter = AnIterNumEdWtIt[(int)DualEdgeWeightMode::STEEPEST_EDGE];
789   if (lc_EdWtNumIter > 0)
790     printf("DSE for %12d (%3d%%) iterations\n", lc_EdWtNumIter,
791            (100 * lc_EdWtNumIter) / AnIterNumIter);
792   lc_EdWtNumIter = AnIterNumEdWtIt[(int)DualEdgeWeightMode::DEVEX];
793   if (lc_EdWtNumIter > 0)
794     printf("Dvx for %12d (%3d%%) iterations\n", lc_EdWtNumIter,
795            (100 * lc_EdWtNumIter) / AnIterNumIter);
796   lc_EdWtNumIter = AnIterNumEdWtIt[(int)DualEdgeWeightMode::DANTZIG];
797   if (lc_EdWtNumIter > 0)
798     printf("Dan for %12d (%3d%%) iterations\n", lc_EdWtNumIter,
799            (100 * lc_EdWtNumIter) / AnIterNumIter);
800   for (int k = 0; k < NUM_ANALYSIS_OPERATION_TYPE; k++) {
801     AnIterOpRec& AnIter = AnIterOp[k];
802     int lcNumCa = AnIter.AnIterOpNumCa;
803     printf("\n%-10s performed %d times\n", AnIter.AnIterOpName.c_str(),
804            AnIter.AnIterOpNumCa);
805     if (lcNumCa > 0) {
806       int lcHyperOp = AnIter.AnIterOpNumHyperOp;
807       int lcHyperRs = AnIter.AnIterOpNumHyperRs;
808       int pctHyperOp = (100 * lcHyperOp) / lcNumCa;
809       int pctHyperRs = (100 * lcHyperRs) / lcNumCa;
810       double lcRsDensity =
811           pow(10.0, AnIter.AnIterOpSumLog10RsDensity / lcNumCa);
812       int lcAnIterOpRsDim = AnIter.AnIterOpRsDim;
813       int lcNumNNz = lcRsDensity * lcAnIterOpRsDim;
814       printf("%12d hyper-sparse operations (%3d%%)\n", lcHyperOp, pctHyperOp);
815       printf("%12d hyper-sparse results    (%3d%%)\n", lcHyperRs, pctHyperRs);
816       printf("%12g density of result (%d / %d nonzeros)\n", lcRsDensity,
817              lcNumNNz, lcAnIterOpRsDim);
818       printValueDistribution(AnIter.AnIterOp_density, AnIter.AnIterOpRsDim);
819     }
820   }
821   int NumInvert = 0;
822 
823   int last_invert_hint = INVERT_HINT_Count - 1;
824   for (int k = 1; k <= last_invert_hint; k++) NumInvert += AnIterNumInvert[k];
825   if (NumInvert > 0) {
826     int lcNumInvert = 0;
827     printf("\nInvert    performed %d times: average frequency = %d\n",
828            NumInvert, AnIterNumIter / NumInvert);
829     lcNumInvert = AnIterNumInvert[INVERT_HINT_UPDATE_LIMIT_REACHED];
830     if (lcNumInvert > 0)
831       printf("%12d (%3d%%) Invert operations due to update limit reached\n",
832              lcNumInvert, (100 * lcNumInvert) / NumInvert);
833     lcNumInvert = AnIterNumInvert[INVERT_HINT_SYNTHETIC_CLOCK_SAYS_INVERT];
834     if (lcNumInvert > 0)
835       printf("%12d (%3d%%) Invert operations due to pseudo-clock\n",
836              lcNumInvert, (100 * lcNumInvert) / NumInvert);
837     lcNumInvert = AnIterNumInvert[INVERT_HINT_POSSIBLY_OPTIMAL];
838     if (lcNumInvert > 0)
839       printf("%12d (%3d%%) Invert operations due to possibly optimal\n",
840              lcNumInvert, (100 * lcNumInvert) / NumInvert);
841     lcNumInvert = AnIterNumInvert[INVERT_HINT_POSSIBLY_PRIMAL_UNBOUNDED];
842     if (lcNumInvert > 0)
843       printf(
844           "%12d (%3d%%) Invert operations due to possibly primal unbounded\n",
845           lcNumInvert, (100 * lcNumInvert) / NumInvert);
846     lcNumInvert = AnIterNumInvert[INVERT_HINT_POSSIBLY_DUAL_UNBOUNDED];
847     if (lcNumInvert > 0)
848       printf("%12d (%3d%%) Invert operations due to possibly dual unbounded\n",
849              lcNumInvert, (100 * lcNumInvert) / NumInvert);
850     lcNumInvert = AnIterNumInvert[INVERT_HINT_POSSIBLY_SINGULAR_BASIS];
851     if (lcNumInvert > 0)
852       printf("%12d (%3d%%) Invert operations due to possibly singular basis\n",
853              lcNumInvert, (100 * lcNumInvert) / NumInvert);
854     lcNumInvert =
855         AnIterNumInvert[INVERT_HINT_PRIMAL_INFEASIBLE_IN_PRIMAL_SIMPLEX];
856     if (lcNumInvert > 0)
857       printf(
858           "%12d (%3d%%) Invert operations due to primal infeasible in primal "
859           "simplex\n",
860           lcNumInvert, (100 * lcNumInvert) / NumInvert);
861   }
862   int suPrice = num_col_price + num_row_price + num_row_price_with_switch;
863   if (suPrice > 0) {
864     printf("\n%12d Price operations:\n", suPrice);
865     printf("%12d Col Price      (%3d%%)\n", num_col_price,
866            (100 * num_col_price) / suPrice);
867     printf("%12d Row Price      (%3d%%)\n", num_row_price,
868            (100 * num_row_price) / suPrice);
869     printf("%12d Row PriceWSw   (%3d%%)\n", num_row_price_with_switch,
870            (100 * num_row_price_with_switch / suPrice));
871   }
872   printf("\n%12d (%3d%%) costly DSE        iterations\n", AnIterNumCostlyDseIt,
873          (100 * AnIterNumCostlyDseIt) / AnIterNumIter);
874 
875   // Look for any Devex data to summarise
876   if (num_devex_framework) {
877     printf("\nDevex summary\n");
878     printf("%12d Devex frameworks\n", num_devex_framework);
879     printf(
880         "%12d average number of iterations\n",
881         AnIterNumEdWtIt[(int)DualEdgeWeightMode::DEVEX] / num_devex_framework);
882   }
883   // Look for any PAMI data to summarise
884   if (sum_multi_chosen > 0) {
885     const int pct_minor_iterations_performed =
886         (100 * sum_multi_finished) / sum_multi_chosen;
887     printf("\nPAMI summary: for average of %0.1g threads \n",
888            average_num_threads);
889     printf("%12d Major iterations\n", multi_iteration_count);
890     printf("%12d Minor iterations\n", sum_multi_finished);
891     printf(
892         "%12d Total rows chosen: performed %3d%% of possible minor "
893         "iterations\n\n",
894         sum_multi_chosen, pct_minor_iterations_performed);
895   }
896 
897   printf("\nCost perturbation summary\n");
898   printValueDistribution(cost_perturbation1_distribution);
899   printValueDistribution(cost_perturbation2_distribution);
900 
901   printValueDistribution(before_ftran_upper_sparse_density, numRow);
902   printValueDistribution(ftran_upper_sparse_density, numRow);
903   printValueDistribution(before_ftran_upper_hyper_density, numRow);
904   printValueDistribution(ftran_upper_hyper_density, numRow);
905   printValueDistribution(primal_step_distribution);
906   printValueDistribution(dual_step_distribution);
907   printValueDistribution(simplex_pivot_distribution);
908   printValueDistribution(factor_pivot_threshold_distribution);
909   printValueDistribution(numerical_trouble_distribution);
910   printValueDistribution(cleanup_dual_change_distribution);
911   printValueDistribution(cleanup_primal_step_distribution);
912   printValueDistribution(cleanup_dual_step_distribution);
913   printValueDistribution(cleanup_primal_change_distribution);
914 
915   if (AnIterTraceIterDl >= 100) {
916     // Possibly (usually) add a temporary record for the final
917     // iterations: may end up with one more than
918     // AN_ITER_TRACE_MX_NUM_REC records, so ensure that there is
919     // enough space in the arrays
920     //
921     const bool add_extra_record =
922         simplex_iteration_count >
923         AnIterTrace[AnIterTraceNumRec].AnIterTraceIter;
924     if (add_extra_record) {
925       AnIterTraceNumRec++;
926       AnIterTraceRec& lcAnIter = AnIterTrace[AnIterTraceNumRec];
927       lcAnIter.AnIterTraceIter = simplex_iteration_count;
928       lcAnIter.AnIterTraceTime = timer_->getWallTime();
929       if (average_fraction_of_possible_minor_iterations_performed > 0) {
930         lcAnIter.AnIterTraceMulti =
931             average_fraction_of_possible_minor_iterations_performed;
932       } else {
933         lcAnIter.AnIterTraceMulti = 0;
934       }
935       lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN] =
936           col_aq_density;
937       lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_BTRAN_EP] =
938           row_ep_density;
939       lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_PRICE_AP] =
940           row_ap_density;
941       lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_BFRT] =
942           col_aq_density;
943       if (edge_weight_mode == DualEdgeWeightMode::STEEPEST_EDGE) {
944         lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE] =
945             row_DSE_density;
946         lcAnIter.AnIterTraceCostlyDse = AnIterCostlyDseMeasure;
947       } else {
948         lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE] = 0;
949         lcAnIter.AnIterTraceCostlyDse = 0;
950       }
951       lcAnIter.AnIterTrace_dual_edge_weight_mode = (int)edge_weight_mode;
952     }
953     // Determine whether the Multi and steepest edge columns should be reported
954     double su_multi_values = 0;
955     double su_dse_values = 0;
956     for (int rec = 1; rec <= AnIterTraceNumRec; rec++) {
957       AnIterTraceRec& lcAnIter = AnIterTrace[rec];
958       su_multi_values += fabs(lcAnIter.AnIterTraceMulti);
959       su_dse_values +=
960           fabs(lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE]);
961     }
962     const bool report_multi = su_multi_values > 0;
963     const bool rp_dual_steepest_edge = su_dse_values > 0;
964     printf("\n Iteration speed analysis\n");
965     AnIterTraceRec& lcAnIter = AnIterTrace[0];
966     int fmIter = lcAnIter.AnIterTraceIter;
967     double fmTime = lcAnIter.AnIterTraceTime;
968     printf("        Iter (      FmIter:      ToIter)      Time      Iter/sec ");
969     if (report_multi) printf("| PAMI ");
970     printf("| C_Aq R_Ep R_Ap ");
971     if (rp_dual_steepest_edge) printf(" DSE ");
972     printf("| EdWt ");
973     if (rp_dual_steepest_edge) {
974       printf("| CostlyDse\n");
975     } else {
976       printf("\n");
977     }
978 
979     for (int rec = 1; rec <= AnIterTraceNumRec; rec++) {
980       AnIterTraceRec& lcAnIter = AnIterTrace[rec];
981       int toIter = lcAnIter.AnIterTraceIter;
982       double toTime = lcAnIter.AnIterTraceTime;
983       int dlIter = toIter - fmIter;
984       if (rec < AnIterTraceNumRec && dlIter != AnIterTraceIterDl)
985         printf("STRANGE: %d = dlIter != AnIterTraceIterDl = %d\n", dlIter,
986                AnIterTraceIterDl);
987       double dlTime = toTime - fmTime;
988       int iterSpeed = 0;
989       if (dlTime > 0) iterSpeed = dlIter / dlTime;
990       int lc_dual_edge_weight_mode = lcAnIter.AnIterTrace_dual_edge_weight_mode;
991       std::string str_dual_edge_weight_mode;
992       if (lc_dual_edge_weight_mode == (int)DualEdgeWeightMode::STEEPEST_EDGE)
993         str_dual_edge_weight_mode = "DSE";
994       else if (lc_dual_edge_weight_mode == (int)DualEdgeWeightMode::DEVEX)
995         str_dual_edge_weight_mode = "Dvx";
996       else if (lc_dual_edge_weight_mode == (int)DualEdgeWeightMode::DANTZIG)
997         str_dual_edge_weight_mode = "Dan";
998       else
999         str_dual_edge_weight_mode = "XXX";
1000       printf("%12d (%12d:%12d) %9.4f  %12d ", dlIter, fmIter, toIter, dlTime,
1001              iterSpeed);
1002       if (report_multi) {
1003         const int pct = (100 * lcAnIter.AnIterTraceMulti);
1004         printf("|  %3d ", pct);
1005       }
1006       printf("|");
1007       reportOneDensity(
1008           lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN]);
1009       reportOneDensity(
1010           lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_BTRAN_EP]);
1011       reportOneDensity(
1012           lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_PRICE_AP]);
1013       double use_row_DSE_density;
1014       if (rp_dual_steepest_edge) {
1015         if (lc_dual_edge_weight_mode ==
1016             (int)DualEdgeWeightMode::STEEPEST_EDGE) {
1017           use_row_DSE_density =
1018               lcAnIter.AnIterTraceDensity[ANALYSIS_OPERATION_TYPE_FTRAN_DSE];
1019         } else {
1020           use_row_DSE_density = 0;
1021         }
1022         reportOneDensity(use_row_DSE_density);
1023       }
1024       printf(" |  %3s ", str_dual_edge_weight_mode.c_str());
1025       if (rp_dual_steepest_edge) {
1026         double use_costly_dse;
1027         printf("|     ");
1028         if (lc_dual_edge_weight_mode ==
1029             (int)DualEdgeWeightMode::STEEPEST_EDGE) {
1030           use_costly_dse = lcAnIter.AnIterTraceCostlyDse;
1031         } else {
1032           use_costly_dse = 0;
1033         }
1034         reportOneDensity(use_costly_dse);
1035       }
1036       printf("\n");
1037       fmIter = toIter;
1038       fmTime = toTime;
1039     }
1040     printf("\n");
1041     // Remove any temporary record added for the final iterations
1042     if (add_extra_record) AnIterTraceNumRec--;
1043   }
1044 }
1045 #endif
1046 
iterationReport(const bool header)1047 void HighsSimplexAnalysis::iterationReport(const bool header) {
1048   if (!(iteration_report_message_level & message_level)) return;
1049   if (!header && (pivotal_row_index < 0 || entering_variable < 0)) return;
1050   reportAlgorithmPhaseIterationObjective(header,
1051                                          iteration_report_message_level);
1052 #ifdef HiGHSDEV
1053   reportDensity(header, iteration_report_message_level);
1054   reportIterationData(header, iteration_report_message_level);
1055 #endif
1056   HighsPrintMessage(output, message_level, iteration_report_message_level,
1057                     "\n");
1058   if (!header) num_iteration_report_since_last_header++;
1059 }
1060 
reportAlgorithmPhaseIterationObjective(const bool header,const int this_message_level)1061 void HighsSimplexAnalysis::reportAlgorithmPhaseIterationObjective(
1062     const bool header, const int this_message_level) {
1063   if (header) {
1064     HighsPrintMessage(output, message_level, this_message_level,
1065                       "       Iteration        Objective    ");
1066   } else {
1067     std::string algorithm;
1068     if (dualAlgorithm()) {
1069       algorithm = "Du";
1070     } else {
1071       algorithm = "Pr";
1072     }
1073     HighsPrintMessage(output, message_level, this_message_level,
1074                       "%2sPh%1d %10d %20.10e", algorithm.c_str(), solve_phase,
1075                       simplex_iteration_count, objective_value);
1076   }
1077 }
1078 
reportInfeasibility(const bool header,const int this_message_level)1079 void HighsSimplexAnalysis::reportInfeasibility(const bool header,
1080                                                const int this_message_level) {
1081   if (header) {
1082     HighsPrintMessage(output, message_level, this_message_level,
1083                       " Infeasibilities num(sum)");
1084   } else {
1085     if (solve_phase == 1) {
1086       HighsPrintMessage(output, message_level, this_message_level,
1087                         " Ph1: %d(%g)", num_primal_infeasibilities,
1088                         sum_primal_infeasibilities);
1089     } else {
1090       HighsPrintMessage(output, message_level, this_message_level,
1091                         " Pr: %d(%g)", num_primal_infeasibilities,
1092                         sum_primal_infeasibilities);
1093     }
1094     if (sum_dual_infeasibilities > 0) {
1095       HighsPrintMessage(output, message_level, this_message_level,
1096                         "; Du: %d(%g)", num_dual_infeasibilities,
1097                         sum_dual_infeasibilities);
1098     }
1099   }
1100 }
1101 
1102 #ifdef HiGHSDEV
reportThreads(const bool header,const int this_message_level)1103 void HighsSimplexAnalysis::reportThreads(const bool header,
1104                                          const int this_message_level) {
1105   if (header) {
1106     HighsPrintMessage(output, message_level, this_message_level, "  Threads");
1107   } else if (num_threads > 0) {
1108     HighsPrintMessage(output, message_level, this_message_level, " %2d|%2d|%2d",
1109                       min_threads, num_threads, max_threads);
1110   } else {
1111     HighsPrintMessage(output, message_level, this_message_level, "   |  |  ");
1112   }
1113 }
1114 
reportMulti(const bool header,const int this_message_level)1115 void HighsSimplexAnalysis::reportMulti(const bool header,
1116                                        const int this_message_level) {
1117   if (header) {
1118     HighsPrintMessage(output, message_level, this_message_level, "  Multi");
1119   } else if (average_fraction_of_possible_minor_iterations_performed >= 0) {
1120     HighsPrintMessage(
1121         output, message_level, this_message_level, "   %3d%%",
1122         (int)(100 * average_fraction_of_possible_minor_iterations_performed));
1123   } else {
1124     HighsPrintMessage(output, message_level, this_message_level, "       ");
1125   }
1126 }
1127 
reportOneDensity(const int this_message_level,const double density)1128 void HighsSimplexAnalysis::reportOneDensity(const int this_message_level,
1129                                             const double density) {
1130   const int log_10_density = intLog10(density);
1131   if (log_10_density > -99) {
1132     HighsPrintMessage(output, message_level, this_message_level, " %4d",
1133                       log_10_density);
1134   } else {
1135     HighsPrintMessage(output, message_level, this_message_level, "     ");
1136   }
1137 }
1138 
reportOneDensity(const double density)1139 void HighsSimplexAnalysis::reportOneDensity(const double density) {
1140   const int log_10_density = intLog10(density);
1141   if (log_10_density > -99) {
1142     printf(" %4d", log_10_density);
1143   } else {
1144     printf("     ");
1145   }
1146 }
1147 
reportDensity(const bool header,const int this_message_level)1148 void HighsSimplexAnalysis::reportDensity(const bool header,
1149                                          const int this_message_level) {
1150   const bool rp_dual_steepest_edge =
1151       edge_weight_mode == DualEdgeWeightMode::STEEPEST_EDGE;
1152   if (header) {
1153     HighsPrintMessage(output, message_level, this_message_level,
1154                       " C_Aq R_Ep R_Ap");
1155     if (rp_dual_steepest_edge) {
1156       HighsPrintMessage(output, message_level, this_message_level, "  DSE");
1157     } else {
1158       HighsPrintMessage(output, message_level, this_message_level, "     ");
1159     }
1160   } else {
1161     reportOneDensity(this_message_level, col_aq_density);
1162     reportOneDensity(this_message_level, row_ep_density);
1163     reportOneDensity(this_message_level, row_ap_density);
1164     double use_row_DSE_density;
1165     if (rp_dual_steepest_edge) {
1166       use_row_DSE_density = row_DSE_density;
1167     } else {
1168       use_row_DSE_density = 0;
1169     }
1170     reportOneDensity(this_message_level, use_row_DSE_density);
1171   }
1172 }
1173 
reportInvert(const bool header,const int this_message_level)1174 void HighsSimplexAnalysis::reportInvert(const bool header,
1175                                         const int this_message_level) {
1176   if (header) {
1177     HighsPrintMessage(output, message_level, this_message_level, " Inv");
1178   } else {
1179     HighsPrintMessage(output, message_level, this_message_level, "  %2d",
1180                       invert_hint);
1181   }
1182 }
1183 
1184 #ifdef HiGHSDEV
reportCondition(const bool header,const int this_message_level)1185 void HighsSimplexAnalysis::reportCondition(const bool header,
1186                                            const int this_message_level) {
1187   if (header) {
1188     HighsPrintMessage(output, message_level, this_message_level, "       k(B)");
1189   } else {
1190     HighsPrintMessage(output, message_level, this_message_level, " %10.4g",
1191                       basis_condition);
1192   }
1193 }
1194 #endif
1195 
reportIterationData(const bool header,const int this_message_level)1196 void HighsSimplexAnalysis::reportIterationData(const bool header,
1197                                                const int this_message_level) {
1198   if (header) {
1199     HighsPrintMessage(output, message_level, this_message_level,
1200                       "       NumCk     LvR     LvC     EnC        DlPr    "
1201                       "    ThDu        ThPr          Aa");
1202   } else {
1203     HighsPrintMessage(output, message_level, this_message_level,
1204                       " %11.4g %7d %7d %7d %11.4g %11.4g %11.4g %11.4g",
1205                       numerical_trouble, pivotal_row_index, leaving_variable,
1206                       entering_variable, primal_delta, dual_step, primal_step,
1207                       pivot_value_from_column);
1208   }
1209 }
1210 
intLog10(const double v)1211 int HighsSimplexAnalysis::intLog10(const double v) {
1212   int intLog10V = -99;
1213   if (v > 0) intLog10V = log(v) / log(10.0);
1214   return intLog10V;
1215 }
1216 
1217 #endif
1218 
dualAlgorithm()1219 bool HighsSimplexAnalysis::dualAlgorithm() {
1220   return (simplex_strategy == SIMPLEX_STRATEGY_DUAL ||
1221           simplex_strategy == SIMPLEX_STRATEGY_DUAL_TASKS ||
1222           simplex_strategy == SIMPLEX_STRATEGY_DUAL_MULTI);
1223 }
1224