glpk_solver.hpp

#ifndef __GLPK_SOLVER_HPP_
#define __GLPK_SOLVER_HPP_

/*****************************************************************************\
* This file is part of DynGB.                                                 *
*                                                                             *
* DynGB is free software: you can redistribute it and/or modify               *
* it under the terms of the GNU General Public License as published by        *
* the Free Software Foundation, either version 2 of the License, or           *
* (at your option) any later version.                                         *
*                                                                             *
* DynGB is distributed in the hope that it will be useful,                    *
* but WITHOUT ANY WARRANTY; without even the implied warranty of              *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               *
* GNU General Public License for more details.                                *
*                                                                             *
* You should have received a copy of the GNU General Public License           *
* along with DynGB. If not, see <http://www.gnu.org/licenses/>.               *
\*****************************************************************************/

#include <glpk.h>

#include "system_constants.hpp"

#include "lp_solver.hpp"

namespace LP_Solvers {

/**
  @brief approximate skeleton of a polyhedral cone, using GLPK linear solver
  @author John Perry
  \version 1.0
  @date January 2017
  @copyright The University of Southern Mississippi
  @ingroup CLSSolvers
  @details This class serves as an interface to GLPK @cite glpk, which we can
      use to find an approximate skeleton to a polyhedral cone.
*/
class GLPK_Solver : public LP_Solver {
public:
  /** @name Construction */
  ///@{
  /** @brief initializes solver for @f$ n @f$ variables */
  explicit GLPK_Solver(NVAR_TYPE n);
  /** @brief copy constructor (deep copy) */
  explicit GLPK_Solver(const GLPK_Solver &);
  virtual bool copy(const LP_Solver *) override;
  ///@}
  /** @name Destruction */
  ///@{
  virtual ~GLPK_Solver();
  ///@}
  /** @name Basic properties */
  ///@{
  virtual NVAR_TYPE get_dimension() const override { return n; }
  virtual unsigned long get_number_of_rays() override {
    unsigned long result;
    if (dirty) result = get_rays().size();
    else result = rays.size();
    return result;
  }
  virtual const set<Ray> & get_rays() const override;
  virtual unsigned long get_number_of_constraints() const override { return m; }
  ///@}
  /** @name Modification */
  ///@{
  virtual bool solve(const Constraint &) override;
  virtual bool solve(const vector<Constraint> &) override;
  GLPK_Solver & operator=(const GLPK_Solver &);
  ///@}
  /** @name Computation */
  ///@{
  /** @brief tests for consistency of a constraint generated by two monomials. */
  virtual inline bool makes_consistent_constraint(
      const Monomial & t, const Monomial & u, bool show_data = false
  ) const override {
    bool inconsistent = true;
    for (
          auto riter = rays.begin();
          inconsistent and riter != rays.end();
          ++riter
    ) {
      int d = 0;
      for (int i = 0; i < riter->get_dimension(); ++i)
        d += (*riter)[i] * (t.degree(i) - u.degree(i));
      if (d > 0)
        inconsistent = false;
      if (show_data) {
        cout << d << ' ';
        if (!inconsistent) cout << *riter << endl;
      }
    }
    if (show_data) cout << endl;
    return not inconsistent;
  }
  ///@}
  /** @name I/O */
  ///@{
  friend ostream & operator<<(ostream & os, const GLPK_Solver & s);
  ///@}
protected:
  /**
    @brief overrides @c rays from @c LP_Solver ; necessary because for GLPK
      it is helpful to avoid updating the rays until they are needed
    @details You shouldn&rsquo;t care about the details, but here we go.
      When @c GLPK_Solver adds one or more constraints, it does not actually
      solve until requested. Until then, the rays are marked as
      &ldquo;dirty&rdquo;. Any attempt to access dirty rays prompts the
      solver. In some cases, the solver will have been marked as
      &ldquo;const&rdquo;; this is inherited from @c LP_Solver and I think
      it&rsquo;s A Good Idea<sup>TM</sup> for exact solvers.
      So I mark the rays here as @c mutable. Arguably I should do that in
      @c LP_Solver but I&rsquo;m not really in the mood to argue at the moment.
  */
  //mutable set<Ray> rays;
private:
  glp_prob * lp; /**< GLPK problem interface */
  glp_smcp smcp; /**< GLPK solver options */
  double * row_data; /**< an array to hold values when creating sparse matrix */
  int * row_indx; /**< an array to hold indices when creating sparse matrix */
  RAYENT_TYPE * ray_data; /**< to retrieve solution */
  int m; /**< number of constraints */
  int n; /**< number of variables */
  mutable bool dirty; /**< whether the rays are valid (here, @c false means valid) */
};

/**
  @brief prints out the constraints, then the rays, then the edges of @p s.
  @param os output stream to print to
  @param s skeleton to print
  @return the output stream
*/
ostream & operator<<(ostream & os, const GLPK_Solver & s);

}

#endif