src/simplex/HDual.h

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                       */
/*    This file is part of the HiGHS linear optimization suite           */
/*                                                                       */
/*    Written and engineered 2008-2021 at the University of Edinburgh    */
/*                                                                       */
/*    Available as open-source under the MIT License                     */
/*                                                                       */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/**@file simplex/HDual.h
 * @brief Dual simplex solver for HiGHS
 * @author Julian Hall, Ivet Galabova, Qi Huangfu and Michael Feldmeier
 */
#ifndef SIMPLEX_HDUAL_H_
#define SIMPLEX_HDUAL_H_

#include <set>
#include <string>
#include <vector>

#include "HConfig.h"
#include "lp_data/HighsModelObject.h"
#include "simplex/HCrash.h"
#include "simplex/HDualRHS.h"
#include "simplex/HDualRow.h"
#include "simplex/HSimplex.h"
#include "simplex/HVector.h"

class HFactor;

/**
 * Limit on the number of column slices for parallel calculations. SIP uses
 * num_threads-2 slices; PAMI uses num_threads-1 slices
 */
const int HIGHS_SLICED_LIMIT =
    HIGHS_THREAD_LIMIT;  // Was 100, but can't see why this should be higher
                         // than HIGHS_THREAD_LIMIT;

/**
 * Parameters controlling number of Devex iterations.
 *
 * There is a new Devex framework if either
 *
 * 1) The weight inaccuracy ratio exceeds maxAllowedDevexWeightRatio
 *
 * 2) There have been max(minAbsNumberDevexIterations,
 * numRow/minRlvNumberDevexIterations) Devex iterations
 */
const int minAbsNumberDevexIterations = 25;
const double minRlvNumberDevexIterations = 1e-2;
const double maxAllowedDevexWeightRatio = 3.0;

/**
 * Multiplier used in running average calculations
 */
const double runningAverageMu = 0.05;

/**
 * Candidate persistence cut-off in PAMI
 */
const double pami_cutoff = 0.95;

/**
 * @brief Dual simplex solver for HiGHS
 */
class HDual {
 public:
  HDual(HighsModelObject& model_object)
      : workHMO(model_object), dualRow(model_object), dualRHS(model_object) {
    dualRow.setup();
    for (int i = 0; i < HIGHS_SLICED_LIMIT; i++)
      slice_dualRow.push_back(HDualRow(model_object));
    dualRHS.setup();
  }

  /**
   * @brief Solve a model instance
   */
  HighsStatus solve();

  const SimplexAlgorithm algorithm = SimplexAlgorithm::DUAL;

 public:
  /**
   * @brief Set solver options from simplex options
   */
  void options();
  /**
   * @brief Initialise a dual simplex instance
   *
   * Copy dimensions and pointers to matrix, factor and solver-related
   * model data, plus tolerances. Sets up local std::vectors (columnDSE,
   * columnBFRT, column, row_ep and row_ap), scalars for their average
   * density and buffers for dualRow and dualRHS.
   */
  void init();

  /**
   * @brief Initialise parallel aspects of a dual simplex instance
   *
   * Sets up data structures for SIP or PAMI
   */
  void initParallel();

  /**
   * @brief Initialise matrix slices and slices of row_ap or dualRow for SIP or
   * PAMI
   *
   * TODO generalise call slice_matrix[i].setup_lgBs so slice can be
   * used with non-logical initial basis
   */
  void initSlice(
      const int init_sliced_num  //!< Ideal number of slices - true number
                                 //!< is modified in light of limits
  );

  /**
   * @brief Perform Phase 1 dual simplex iterations
   */
  void solvePhase1();

  /**
   * @brief Perform Phase 2 dual simplex iterations
   */
  void solvePhase2();

  /**
   * @brief Reinvert if INVERT not fresh, then recompute dual and primal values
   *
   * Also collects primal infeasibilities and computes the dual objective value
   */

  void rebuild();

  /**
   * @brief Remove perturbation and recompute the dual solution
   *
   * Also collects primal infeasibilities and computes the dual objective value
   */
  void cleanup();

  /**
   * @brief Perform a single serial dual simplex iteration
   *
   * All the methods it calls have as their first line "if (invertHint)
   * return;", where invertHint is, for example, set to 1 when CHUZR
   * finds no candidate. This causes a break from the inner loop of
   * solve_phase% and, hence, a call to rebuild().
   */
  void iterate();

  /**
   * @brief Perform a single SIP dual simplex iteration
   */
  void iterateTasks();

  /**
   * @brief Perform a single PAMI dual simplex iteration - source code in
   * HDualMulti.cpp
   */
  void iterateMulti();  // in HDualMulti.cpp

  /**
   * @brief Pass the data for the serial iteration analysis, report and rebuild
   * report
   */
  void iterationAnalysisData();

  /**
   * @brief Perform the serial iteration analysis
   */
  void iterationAnalysis();

  /**
   * @brief Pass the data for the PAMI iteration analysis for a minor iteration,
   * report and rebuild report
   */
  void iterationAnalysisMinorData();

  /**
   * @brief Perform the PAMI iteration analysis for a minor iteration
   */
  void iterationAnalysisMinor();

  /**
   * @brief Pass the data for the PAMI iteration analysis for a major iteration
   */
  void iterationAnalysisMajorData();

  /**
   * @brief Perform the PAMI iteration analysis for a major iteration
   */
  void iterationAnalysisMajor();

  /**
   * @brief Single line report after rebuild
   */
  void reportRebuild(const int rebuild_invert_hint = -1);

  /**
   * @brief Choose the index of a good row to leave the basis (CHUZR)
   */
  void chooseRow();

  /**
   * @brief Determine whether the updated_edge_weight is accurate enough to
   * be accepted, and update the analysis of weight errors
   */
  bool acceptDualSteepestEdgeWeight(const double updated_edge_weight);

  /**
   * @brief Determine whether the updated_edge_weight error should trigger a new
   * Devex framework
   */
  bool newDevexFramework(const double updated_edge_weight);

  /**
   * @brief Compute pivot row (PRICE) and choose the index of a good column to
   * enter the basis (CHUZC)
   */
  void chooseColumn(HVector* row_ep);

  /**
   * @brief Choose the index of a good column to enter the basis (CHUZC) by
   * exploiting slices of the pivotal row - for SIP and PAMI
   */
  void chooseColumnSlice(HVector* row_ep);

  /**
   * @brief Compute the pivotal column (FTRAN)
   */
  void updateFtran();

  /**
   * @brief Compute the RHS changes corresponding to the BFRT
   * (FTRAN-BFRT)
   */
  void updateFtranBFRT();

  /**
   * @brief Compute the std::vector required to update DSE weights - being
   * FTRAN applied to the pivotal column (FTRAN-DSE)
   */
  void updateFtranDSE(HVector* DSE_Vector  //!< Pivotal column as RHS for FTRAN
  );
  /**
   * @brief Compare the pivot value computed row-wise and column-wise
   * and determine whether reinversion is advisable
   */
  void updateVerify();

  /**
   * @brief Update the dual values
   */
  void updateDual();

  /**
   * @brief Update the primal values and any edge weights
   */
  void updatePrimal(HVector* DSE_Vector  //!< FTRANned pivotal column
  );

  /**
   * @brief Update the basic and nonbasic variables, iteration count,
   * invertible representation of the basis matrix and row-wise
   * representation of the nonbasic columns, delete the Freelist entry
   * for the entering column, update the primal value for the row
   * where the basis change has occurred, and set the corresponding
   * primal infeasibility value in dualRHS.work_infeasibility, and
   * then determine whether to reinvert according to the synthetic
   * clock
   */
  void updatePivots();

  /**
   * @brief Initialise a Devex framework: reference set is all basic
   * variables
   */
  void initialiseDevexFramework(const bool parallel = false);

  /**
   * @brief Interpret the dual edge weight strategy as setting of a mode and
   * other actions
   */
  void interpretDualEdgeWeightStrategy(
      const int simplex_dual_edge_weight_strategy);

  bool reachedExactDualObjectiveValueUpperBound();
  double computeExactDualObjectiveValue();

  /**
   * @brief PAMI: Choose the indices of a good set of rows to leave the
   * basis (CHUZR)
   */
  void majorChooseRow();

  /**
   * @brief PAMI: Perform multiple BTRAN
   */
  void majorChooseRowBtran();

  /**
   * @brief PAMI: Choose the index (from the set of indices) of a good
   * row to leave the basis (CHUZR-MI)
   */
  void minorChooseRow();

  /**
   * @brief PAMI: Update the data during minor iterations
   */
  void minorUpdate();

  /**
   * @brief PAMI: Update the dual values during minor iterations
   */
  void minorUpdateDual();

  /**
   * @brief PAMI: Update the primal values during minor iterations
   */
  void minorUpdatePrimal();

  /**
   * @brief PAMI: Perform a basis change during minor iterations
   */
  void minorUpdatePivots();

  /**
   * @brief PAMI: Update the tableau rows during minor iterations
   */
  void minorUpdateRows();

  /**
   * @brief PAMI: Initialise a new Devex framework during minor iterations
   */
  void minorInitialiseDevexFramework();

  /**
   * @brief PAMI: Perform updates after a set of minor iterations
   */
  void majorUpdate();

  /**
   * @brief PAMI: Prepare for the FTRANs after a set of minor iterations
   */
  void majorUpdateFtranPrepare();

  /**
   * @brief PAMI: Perform the parallel part of multiple FTRANs after a
   * set of minor iterations
   */
  void majorUpdateFtranParallel();

  /**
   * @brief PAMI: Perform the final part of multiple FTRANs after a set
   * of minor iterations
   */
  void majorUpdateFtranFinal();

  /**
   * @brief PAMI: Update the primal values after a set of minor
   * iterations
   */
  void majorUpdatePrimal();

  /**
   * @brief PAMI: Update the invertible representation of the basis
   * matrix after a set of minor iterations
   */
  void majorUpdateFactor();

  /**
   * @brief PAMI: Roll back some iterations if numerical trouble
   * detected when updating the invertible representation of the basis
   * matrix after a set of minor iterations
   */
  void majorRollback();

  // private:
  HighsStatus returnFromSolve(const HighsStatus return_status);
  void saveDualRay();
  bool getNonsingularInverse();
  bool getBacktrackingBasis(vector<double>& scattered_edge_weights);
  void putBacktrackingBasis();
  void putBacktrackingBasis(const vector<int>& basicIndex_before_compute_factor,
                            const vector<double>& scattered_edge_weights);

  void assessPhase1Optimality();
  void exitPhase1ResetDuals();
  void reportOnPossibleLpDualInfeasibility();

  bool checkNonUnitWeightError(std::string message);
  bool dualInfoOk(const HighsLp& lp);
  bool bailoutReturn();
  bool bailoutOnTimeIterations();
  bool bailoutOnDualObjective();

  bool solve_bailout;  //!< Set true if control is to be returned immediately to
                       //!< calling function

  // Devex scalars
  int num_devex_iterations =
      0;  //!< Number of Devex iterations with the current framework
  bool new_devex_framework = false;  //!< Set a new Devex framework
  bool minor_new_devex_framework =
      false;  //!< Set a new Devex framework in PAMI minor iterations

  // Model
  HighsModelObject& workHMO;
  int solver_num_row;
  int solver_num_col;
  int solver_num_tot;

  const HMatrix* matrix;
  const HFactor* factor;
  HighsSimplexAnalysis* analysis;

  const int* jMove;
  const double* workRange;
  const double* baseLower;
  const double* baseUpper;
  double* baseValue;
  double* workDual;
  double* workValue;
  double* colLower;
  double* colUpper;
  double* rowLower;
  double* rowUpper;
  int* nonbasicFlag;

  // Options
  DualEdgeWeightMode dual_edge_weight_mode;
  bool initialise_dual_steepest_edge_weights;
  bool allow_dual_steepest_edge_to_devex_switch;

  const double min_dual_steepest_edge_weight = 1e-4;

  double Tp;  // Tolerance for primal
  double primal_feasibility_tolerance;

  double Td;  // Tolerance for dual
  double dual_feasibility_tolerance;
  double dual_objective_value_upper_bound;

  int solvePhase;
  int invertHint;

  HVector row_ep;
  HVector row_ap;
  HVector col_aq;
  HVector col_BFRT;
  HVector col_DSE;

  HDualRow dualRow;

  // Solving related buffers
  int dualInfeasCount;

  HDualRHS dualRHS;

  // Simplex pivotal information
  int rowOut;
  int columnOut;
  int sourceOut;  // -1 from small to lower, +1 to upper
  int columnIn;
  double deltaPrimal;
  double thetaDual;
  double thetaPrimal;
  double alpha;
  double alphaRow;
  double numericalTrouble;
  // (Local) value of computed weight
  double computed_edge_weight;

  bool check_invert_condition = false;

  // Partitioned coefficient matrix
  int slice_num;
  int slice_PRICE;
  int slice_start[HIGHS_SLICED_LIMIT + 1];
  HMatrix slice_matrix[HIGHS_SLICED_LIMIT];
  HVector slice_row_ap[HIGHS_SLICED_LIMIT];
  std::vector<HDualRow> slice_dualRow;

  /**
   * @brief Multiple CHUZR data
   */
  struct MChoice {
    int rowOut;
    double baseValue;
    double baseLower;
    double baseUpper;
    double infeasValue;
    double infeasEdWt;
    double infeasLimit;
    HVector row_ep;
    HVector col_aq;
    HVector col_BFRT;
  };

  /**
   * @brief Multiple minor iteration data
   */
  struct MFinish {
    int moveIn;
    double shiftOut;
    std::vector<int> flipList;

    int rowOut;
    int columnOut;
    int columnIn;
    double alphaRow;
    double thetaPrimal;
    double basicBound;
    double basicValue;
    double EdWt;
    HVector_ptr row_ep;
    HVector_ptr col_aq;
    HVector_ptr col_BFRT;
  };

  int multi_num;
  int multi_chosen;
  int multi_iChoice;
  int multi_nFinish;
  int multi_iteration;
  int multi_chooseAgain;
  MChoice multi_choice[HIGHS_THREAD_LIMIT];
  MFinish multi_finish[HIGHS_THREAD_LIMIT];

#ifdef HiGHSDEV
  const bool rp_reinvert_syntheticClock = false;  // true;//
  const bool rp_numericalTrouble = false;         // true;//
#endif
  const double original_multi_build_syntheticTick_mu = 1.5;
  const double multi_build_syntheticTick_mu = 1.0;
  // original_multi_build_syntheticTick_mu;//
  const double numerical_trouble_tolerance = 1e-7;
  const double original_multi_numerical_trouble_tolerance = 1e-8;
  const double multi_numerical_trouble_tolerance = 1e-7;
  // original_multi_numerical_trouble_tolerance;

  const int synthetic_tick_reinversion_min_update_count = 50;
  const int original_multi_synthetic_tick_reinversion_min_update_count = 201;
  const int multi_synthetic_tick_reinversion_min_update_count =
      synthetic_tick_reinversion_min_update_count;
  // original_multi_synthetic_tick_reinversion_min_update_count;

  double build_syntheticTick;
  double total_syntheticTick;
};

#endif /* SIMPLEX_HDUAL_H_ */