DelpiLpSolver.h

#pragma once
 
#include <string>
#include <utility>
#include <vector>
 
#include "delpi/libs/eigen.h"
#include "delpi/libs/gmp.h"
#include "delpi/solver/LpSolver.h"
#include "delpi/solver/internal/Basis.h"
#include "delpi/symbolic/Expression.h"
#include "delpi/symbolic/Variable.h"
#include "internal/LpProblem.h"
 
namespace delpi {

class DelpiLpSolver final : public LpSolver {
 public:
  explicit DelpiLpSolver(Config config = {}, const std::string& class_name = "DelpiLpSolver");
 
  [[nodiscard]] int num_columns() const override;
  [[nodiscard]] int num_rows() const override;
  [[nodiscard]] const internal::LpProblem& problem() const { return problem_; }
 
  [[nodiscard]] Column column(ColumnIndex column_idx) const override;
  [[nodiscard]] Row row(RowIndex row_idx) const override;
  void ReserveColumns(int num_columns) override;
  void ReserveRows(int num_rows) override;
  ColumnIndex AddColumn(const Variable& var, const mpq_class& obj, const mpq_class& lb, const mpq_class& ub) override;
  RowIndex AddRow(const std::vector<Expression::Addend>& addends, const mpq_class& lb, const mpq_class& ub) override;
  RowIndex AddRow(const Expression::Addends& lhs, FormulaKind sense, const mpq_class& rhs) override;
  void SetBound(Variable var, const mpq_class& lb, const mpq_class& ub) override;
  void SetCoefficient(RowIndex row, ColumnIndex column, const mpq_class& value) override;
  void SetObjective(int column, const mpq_class& value) override;
 
#ifndef NDEBUG
  void Dump() override;
#endif
 
 private:
  LpResult SolveCore() override;
  template <IsAnyOf<double, mpq_class> T>
  LpResult InternalSolve(const Matrix<T>& A, const Vector<T>& b, const Vector<T>& c, const T& tolerance,
                         internal::Basis<T>& basis, Vector<T>* x = nullptr, T* obj = nullptr);
  LpResult FeasibilitySolve(Matrix<mpq_class>& slack_A, Vector<mpq_class>& slack_b,
                            internal::Basis<mpq_class>& slack_basis);
  LpResult OptimalitySolve(const Matrix<mpq_class>& slack_A, const Vector<mpq_class>& slack_b,
                           const Vector<mpq_class>& slack_c, internal::Basis<mpq_class>& slack_basis);
  LpResult FeasibilityCheck(const Matrix<mpq_class>& aux_A, const Vector<mpq_class>& slack_b,
                            const Vector<mpq_class>& aux_c, const internal::Basis<mpq_class>& feas_basis) const;
  LpResult OptimalityCheck(const Matrix<mpq_class>& slack_A, const Vector<mpq_class>& slack_b,
                           const Vector<mpq_class>& slack_c, const internal::Basis<mpq_class>& basis);
  LpResult UnboundednessCheck(const Matrix<mpq_class>& slack_A, const Vector<mpq_class>& slack_b,
                              const Vector<mpq_class>& slack_c, const internal::Basis<mpq_class>& basis) const;
  void RemoveAuxiliaryColumns(const internal::Basis<mpq_class>& feas_basis, const std::vector<Index>& aux_columns,
                              Matrix<mpq_class>& slack_A, Vector<mpq_class>& slack_b,
                              internal::Basis<mpq_class>& slack_basis) const;

  void SlackForm(Matrix<mpq_class>& slack_A, Vector<mpq_class>& slack_c) const;

  internal::Basis<mpq_class> AuxForm(const Matrix<mpq_class>& slack_A, const Vector<mpq_class>& slack_b,
                                     Matrix<mpq_class>& aux_A, Vector<mpq_class>& aux_c,
                                     std::vector<Index>& aux_columns) const;

  void ComputeSlackAndAuxVariables();
 
 protected:
  void EnsureSenseCore() override;
 
 private:
#if 0
  void UpdateInfeasible();
#endif
 
  Vector<mpq_class> x_;          
  internal::LpProblem problem_;  
  mpq_class delta_;              
};
 
std::ostream& operator<<(std::ostream& os, const DelpiLpSolver& solver);
 
}  // namespace delpi