#ifndef SIMPLEX_HAPP_H_
#define SIMPLEX_HAPP_H_
#include <cmath>
#include "lp_data/HighsLpSolverObject.h"
#include "lp_data/HighsLpUtils.h"
#include "lp_data/HighsSolution.h"
#include "lp_data/HighsSolve.h"
#include "simplex/HEkk.h"
#include "simplex/HSimplex.h"
inline HighsStatus returnFromSolveLpSimplex(HighsLpSolverObject& solver_object,
HighsStatus return_status) {
HighsOptions& options = solver_object.options_;
HEkk& ekk_instance = solver_object.ekk_instance_;
HighsLp& incumbent_lp = solver_object.lp_;
solver_object.highs_info_.simplex_iteration_count =
ekk_instance.iteration_count_;
HighsInt sub_solver_ix = -1;
if (std::signbit(solver_object.sub_solver_call_time_
.run_time[kSubSolverDuSimplexBasis]))
sub_solver_ix = kSubSolverDuSimplexBasis;
if (std::signbit(solver_object.sub_solver_call_time_
.run_time[kSubSolverDuSimplexNoBasis]))
sub_solver_ix = kSubSolverDuSimplexNoBasis;
if (std::signbit(solver_object.sub_solver_call_time_
.run_time[kSubSolverPrSimplexBasis]))
sub_solver_ix = kSubSolverPrSimplexBasis;
if (std::signbit(solver_object.sub_solver_call_time_
.run_time[kSubSolverPrSimplexNoBasis]))
sub_solver_ix = kSubSolverPrSimplexNoBasis;
assert(sub_solver_ix >= 0);
if (sub_solver_ix != kSubSolverDuSimplexBasis)
assert(!std::signbit(solver_object.sub_solver_call_time_
.run_time[kSubSolverDuSimplexBasis]));
if (sub_solver_ix != kSubSolverDuSimplexNoBasis)
assert(!std::signbit(solver_object.sub_solver_call_time_
.run_time[kSubSolverDuSimplexNoBasis]));
if (sub_solver_ix != kSubSolverPrSimplexBasis)
assert(!std::signbit(solver_object.sub_solver_call_time_
.run_time[kSubSolverPrSimplexBasis]));
if (sub_solver_ix != kSubSolverPrSimplexNoBasis)
assert(!std::signbit(solver_object.sub_solver_call_time_
.run_time[kSubSolverPrSimplexNoBasis]));
solver_object.sub_solver_call_time_.num_call[sub_solver_ix]++;
solver_object.sub_solver_call_time_.run_time[sub_solver_ix] +=
solver_object.timer_.read();
assert(!incumbent_lp.is_moved_);
assert(!incumbent_lp.is_scaled_);
if (return_status == HighsStatus::kError) {
ekk_instance.clear();
return return_status;
}
assert(ekk_instance.status_.has_invert);
assert(ekk_instance.status_.has_nla);
ekk_instance.setNlaPointersForLpAndScale(incumbent_lp);
assert(ekk_instance.debugNlaScalingOk(incumbent_lp));
HighsInt alt_debug_level = -1;
if (ekk_instance.debugNlaCheckInvert("HApp: returnFromSolveLpSimplex",
alt_debug_level) ==
HighsDebugStatus::kError) {
highsLogUser(options.log_options, HighsLogType::kError,
"Error in basis matrix inverse after solving the LP\n");
return_status = HighsStatus::kError;
}
if (solver_object.model_status_ == HighsModelStatus::kOptimal) {
solver_object.highs_info_.num_complementarity_violations = 0;
solver_object.highs_info_.max_complementarity_violation = 0;
}
return return_status;
}
inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
HighsStatus return_status = HighsStatus::kOk;
HighsStatus call_status;
HighsOptions& options = solver_object.options_;
HighsLp& incumbent_lp = solver_object.lp_;
HighsSolution& solution = solver_object.solution_;
HighsModelStatus& model_status = solver_object.model_status_;
HighsModelStatus scaled_model_status = HighsModelStatus::kUnknown;
HighsInfo& highs_info = solver_object.highs_info_;
HighsBasis& basis = solver_object.basis_;
HEkk& ekk_instance = solver_object.ekk_instance_;
HighsLp& ekk_lp = ekk_instance.lp_;
HighsSimplexInfo& ekk_info = ekk_instance.info_;
SimplexBasis& ekk_basis = ekk_instance.basis_;
HighsSimplexStatus& status = ekk_instance.status_;
bool retained_ekk_data_ok = ekk_instance.debugRetainedDataOk(incumbent_lp) !=
HighsDebugStatus::kLogicalError;
if (!retained_ekk_data_ok) {
highsLogUser(options.log_options, HighsLogType::kError,
"solveLpSimplex: Retained Ekk data not OK on entry\n");
assert(retained_ekk_data_ok);
return_status = HighsStatus::kError;
}
HighsInt sub_solver_ix = -1;
if (options.simplex_strategy == kSimplexStrategyPrimal) {
if (basis.valid) {
sub_solver_ix = kSubSolverPrSimplexBasis;
} else {
sub_solver_ix = kSubSolverPrSimplexNoBasis;
}
} else {
if (basis.valid) {
sub_solver_ix = kSubSolverDuSimplexBasis;
} else {
sub_solver_ix = kSubSolverDuSimplexNoBasis;
}
}
assert(sub_solver_ix >= 0);
assert(solver_object.sub_solver_call_time_.run_time.size() > 0);
solver_object.sub_solver_call_time_.run_time[sub_solver_ix] =
-solver_object.timer_.read();
ekk_instance.iteration_count_ = highs_info.simplex_iteration_count;
resetModelStatusAndHighsInfo(solver_object);
ekk_instance.initialiseSimplexStats();
bool positive_num_row = solver_object.lp_.num_row_ > 0;
assert(positive_num_row);
if (!positive_num_row) {
highsLogUser(
options.log_options, HighsLogType::kError,
"solveLpSimplex called for LP with non-positive (%" HIGHSINT_FORMAT
") "
"number of constraints\n",
incumbent_lp.num_row_);
return returnFromSolveLpSimplex(solver_object, HighsStatus::kError);
}
assert(!incumbent_lp.is_scaled_);
assert(!incumbent_lp.is_moved_);
considerScaling(options, incumbent_lp);
const bool was_scaled = incumbent_lp.is_scaled_;
if (!status.has_basis && !basis.valid && basis.useful) {
assert(basis.col_status.size() ==
static_cast<size_t>(incumbent_lp.num_col_));
assert(basis.row_status.size() ==
static_cast<size_t>(incumbent_lp.num_row_));
HighsStatus return_status = formSimplexLpBasisAndFactor(solver_object);
if (return_status != HighsStatus::kOk)
return returnFromSolveLpSimplex(solver_object, HighsStatus::kError);
refineBasis(incumbent_lp, solution, basis);
basis.valid = true;
}
assert(incumbent_lp.is_scaled_ == was_scaled);
ekk_instance.moveLp(solver_object);
if (!status.has_basis) {
if (basis.valid) {
call_status = ekk_instance.setBasis(basis);
if (call_status == HighsStatus::kError) {
incumbent_lp.moveBackLpAndUnapplyScaling(ekk_lp);
return returnFromSolveLpSimplex(solver_object, call_status);
}
} else {
if (options.simplex_dualize_strategy == kHighsOptionChoose ||
options.simplex_dualize_strategy == kHighsOptionOn) {
bool dualize_lp = true;
if (options.simplex_dualize_strategy == kHighsOptionChoose) {
if (incumbent_lp.num_row_ < 10 * incumbent_lp.num_col_)
dualize_lp = false;
}
if (dualize_lp) ekk_instance.dualize();
}
if (options.simplex_permute_strategy == kHighsOptionChoose ||
options.simplex_permute_strategy == kHighsOptionOn) {
ekk_instance.permute();
}
}
}
HighsInt num_unscaled_primal_infeasibilities =
kHighsIllegalInfeasibilityCount;
HighsInt num_unscaled_dual_infeasibilities = kHighsIllegalInfeasibilityCount;
bool solve_unscaled_lp = false;
bool solved_unscaled_lp = false;
if (!incumbent_lp.scale_.has_scaling) {
solve_unscaled_lp = true;
return_status = ekk_instance.solve();
solved_unscaled_lp = true;
ekk_instance.unpermute();
ekk_instance.undualize();
assert(!ekk_instance.status_.is_permuted &&
!ekk_instance.status_.is_dualized);
if (options.cost_scale_factor) {
double cost_scale_factor = pow(2.0, -options.cost_scale_factor);
highsLogDev(options.log_options, HighsLogType::kInfo,
"Objective = %11.4g\n",
cost_scale_factor * ekk_instance.info_.dual_objective_value);
ekk_instance.model_status_ = HighsModelStatus::kNotset;
return_status = HighsStatus::kError;
}
} else {
bool refine_solution = false;
if (options.simplex_unscaled_solution_strategy ==
kSimplexUnscaledSolutionStrategyNone ||
options.simplex_unscaled_solution_strategy ==
kSimplexUnscaledSolutionStrategyRefine) {
assert(incumbent_lp.scale_.has_scaling);
assert(incumbent_lp.is_scaled_);
return_status = ekk_instance.solve();
ekk_instance.unpermute();
ekk_instance.undualize();
assert(!ekk_instance.status_.is_permuted &&
!ekk_instance.status_.is_dualized);
if (options.cost_scale_factor) {
double cost_scale_factor = pow(2.0, -options.cost_scale_factor);
highsLogDev(
options.log_options, HighsLogType::kInfo, "Objective = %11.4g\n",
cost_scale_factor * ekk_instance.info_.dual_objective_value);
ekk_instance.model_status_ = HighsModelStatus::kNotset;
return_status = HighsStatus::kError;
}
if (return_status == HighsStatus::kError) {
incumbent_lp.moveBackLpAndUnapplyScaling(ekk_lp);
return returnFromSolveLpSimplex(solver_object, return_status);
}
scaled_model_status = ekk_instance.model_status_;
highs_info.objective_function_value = ekk_info.primal_objective_value;
highs_info.simplex_iteration_count = ekk_instance.iteration_count_;
solution = ekk_instance.getSolution();
basis = ekk_instance.getHighsBasis(ekk_lp);
assert(basis.valid);
highs_info.basis_validity = kBasisValidityValid;
incumbent_lp.moveBackLpAndUnapplyScaling(ekk_lp);
ekk_instance.setNlaPointersForLpAndScale(incumbent_lp);
unscaleSolution(solution, incumbent_lp.scale_);
getUnscaledInfeasibilities(options, incumbent_lp.scale_, ekk_basis,
ekk_info, highs_info);
num_unscaled_primal_infeasibilities =
highs_info.num_primal_infeasibilities;
num_unscaled_dual_infeasibilities = highs_info.num_dual_infeasibilities;
const bool scaled_optimality_but_unscaled_infeasibilities =
scaled_model_status == HighsModelStatus::kOptimal &&
(num_unscaled_primal_infeasibilities ||
num_unscaled_dual_infeasibilities);
const bool scaled_objective_target_but_unscaled_primal_infeasibilities =
scaled_model_status == HighsModelStatus::kObjectiveTarget &&
highs_info.num_primal_infeasibilities > 0;
const bool scaled_objective_bound_but_unscaled_dual_infeasibilities =
scaled_model_status == HighsModelStatus::kObjectiveBound &&
highs_info.num_dual_infeasibilities > 0;
if (scaled_optimality_but_unscaled_infeasibilities ||
scaled_objective_target_but_unscaled_primal_infeasibilities ||
scaled_objective_bound_but_unscaled_dual_infeasibilities)
highsLogDev(options.log_options, HighsLogType::kInfo,
"After unscaling with status %s, have num/max/sum primal "
"(%" HIGHSINT_FORMAT "/%g/%g) and dual (%" HIGHSINT_FORMAT
"/%g/%g) "
"unscaled infeasibilities\n",
utilModelStatusToString(scaled_model_status).c_str(),
highs_info.num_primal_infeasibilities,
highs_info.max_primal_infeasibility,
highs_info.sum_primal_infeasibilities,
highs_info.num_dual_infeasibilities,
highs_info.max_dual_infeasibility,
highs_info.sum_dual_infeasibilities);
refine_solution =
options.simplex_unscaled_solution_strategy ==
kSimplexUnscaledSolutionStrategyRefine &&
(scaled_optimality_but_unscaled_infeasibilities ||
scaled_model_status == HighsModelStatus::kInfeasible ||
scaled_model_status == HighsModelStatus::kUnboundedOrInfeasible ||
scaled_model_status == HighsModelStatus::kUnbounded ||
scaled_objective_bound_but_unscaled_dual_infeasibilities ||
scaled_objective_target_but_unscaled_primal_infeasibilities ||
scaled_model_status == HighsModelStatus::kUnknown);
if (!refine_solution) {
model_status = scaled_model_status;
return_status = highsStatusFromHighsModelStatus(model_status);
return returnFromSolveLpSimplex(solver_object, return_status);
}
} else {
assert(options.simplex_unscaled_solution_strategy ==
kSimplexUnscaledSolutionStrategyDirect);
incumbent_lp.moveBackLpAndUnapplyScaling(ekk_lp);
}
assert(options.simplex_unscaled_solution_strategy ==
kSimplexUnscaledSolutionStrategyDirect ||
refine_solution);
assert(!incumbent_lp.is_moved_);
assert(!incumbent_lp.is_scaled_);
ekk_instance.moveLp(solver_object);
solve_unscaled_lp = true;
if (scaled_model_status == HighsModelStatus::kInfeasible &&
ekk_instance.exit_algorithm_ == SimplexAlgorithm::kDual)
assert(ekk_instance.dual_ray_record_.index != kNoRayIndex);
if (scaled_model_status == HighsModelStatus::kInfeasible &&
ekk_instance.dual_ray_record_.index != kNoRayIndex) {
ekk_instance.setNlaPointersForLpAndScale(ekk_lp);
if (ekk_instance.proofOfPrimalInfeasibility()) solve_unscaled_lp = false;
}
if (solve_unscaled_lp) {
HighsInt simplex_strategy = options.simplex_strategy;
double dual_simplex_cost_perturbation_multiplier =
options.dual_simplex_cost_perturbation_multiplier;
HighsInt simplex_dual_edge_weight_strategy =
ekk_info.dual_edge_weight_strategy;
if (num_unscaled_primal_infeasibilities == 0 ||
scaled_model_status == HighsModelStatus::kObjectiveTarget) {
options.simplex_strategy = kSimplexStrategyPrimal;
if (scaled_model_status == HighsModelStatus::kObjectiveTarget) {
highsLogDev(options.log_options, HighsLogType::kInfo,
"solveLpSimplex: Calling primal simplex after "
"scaled_model_status == "
"HighsModelStatus::kObjectiveTarget: solve "
"= %d; tick = %d; iter = %d\n",
(int)ekk_instance.debug_solve_call_num_,
(int)ekk_instance.debug_initial_build_synthetic_tick_,
(int)ekk_instance.iteration_count_);
}
} else {
assert(num_unscaled_primal_infeasibilities > 0);
if (options.simplex_strategy != kSimplexStrategyDual)
highsLogDev(
options.log_options, HighsLogType::kInfo,
"Forcing change from %s to %s\n",
ekk_instance.simplexStrategyToString(options.simplex_strategy)
.c_str(),
ekk_instance.simplexStrategyToString(kSimplexStrategyDual)
.c_str());
options.simplex_strategy = kSimplexStrategyDual;
if ((status.has_basis || basis.valid) &&
!status.has_dual_steepest_edge_weights) {
ekk_info.dual_edge_weight_strategy = kSimplexEdgeWeightStrategyDevex;
}
}
const bool force_phase1 =
(options.simplex_unscaled_solution_strategy ==
kSimplexUnscaledSolutionStrategyDirect) ||
(scaled_model_status == HighsModelStatus::kObjectiveTarget);
const bool force_phase2 =
(options.simplex_unscaled_solution_strategy !=
kSimplexUnscaledSolutionStrategyDirect) &&
(scaled_model_status != HighsModelStatus::kObjectiveTarget);
assert(force_phase2 == !force_phase1);
return_status = ekk_instance.solve(force_phase2);
solved_unscaled_lp = true;
if (scaled_model_status != HighsModelStatus::kObjectiveBound &&
ekk_instance.model_status_ == HighsModelStatus::kObjectiveBound) {
const bool objective_bound_refinement =
ekk_info.num_dual_infeasibilities > 0;
if (objective_bound_refinement) {
options.simplex_strategy = kSimplexStrategyPrimal;
return_status = ekk_instance.solve(force_phase2);
}
}
options.simplex_strategy = simplex_strategy;
options.dual_simplex_cost_perturbation_multiplier =
dual_simplex_cost_perturbation_multiplier;
ekk_info.dual_edge_weight_strategy = simplex_dual_edge_weight_strategy;
}
}
if (solved_unscaled_lp) {
scaled_model_status = ekk_instance.model_status_;
highs_info.objective_function_value = ekk_info.primal_objective_value;
highs_info.simplex_iteration_count = ekk_instance.iteration_count_;
solution = ekk_instance.getSolution();
basis = ekk_instance.getHighsBasis(ekk_lp);
assert(basis.valid);
highs_info.basis_validity = kBasisValidityValid;
}
incumbent_lp = std::move(ekk_lp);
incumbent_lp.is_moved_ = false;
ekk_instance.setNlaPointersForLpAndScale(incumbent_lp);
if (return_status == HighsStatus::kError) {
model_status = scaled_model_status;
return_status = highsStatusFromHighsModelStatus(model_status);
assert(return_status == HighsStatus::kError);
return returnFromSolveLpSimplex(solver_object, HighsStatus::kError);
}
if (solved_unscaled_lp) {
assert(solve_unscaled_lp);
highs_info.num_primal_infeasibilities = ekk_info.num_primal_infeasibilities;
highs_info.max_primal_infeasibility = ekk_info.max_primal_infeasibility;
highs_info.sum_primal_infeasibilities = ekk_info.sum_primal_infeasibilities;
highs_info.num_dual_infeasibilities = ekk_info.num_dual_infeasibilities;
highs_info.max_dual_infeasibility = ekk_info.max_dual_infeasibility;
highs_info.sum_dual_infeasibilities = ekk_info.sum_dual_infeasibilities;
} else {
assert(!solve_unscaled_lp);
assert(scaled_model_status == HighsModelStatus::kInfeasible);
}
setSolutionStatus(highs_info);
model_status = scaled_model_status;
return_status = highsStatusFromHighsModelStatus(model_status);
return returnFromSolveLpSimplex(solver_object, return_status);
}
#endif