data_objects.hpp

#ifndef LPINTERFACE_DATA_OBJECTS_H
#define LPINTERFACE_DATA_OBJECTS_H

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <set>
#include <type_traits>
#include <vector>

#include "common.hpp"
#include "errors.hpp"

namespace lpint {

// matrix entry is templated over T, with T restricted to
// arithmetic types i.e. numbers
template <typename T>
class MatrixEntry {
  static_assert(std::is_arithmetic<T>::value,
                "MatrixEntry<T> requires T to be arithmetic");

 private:
  using iterator = typename std::vector<T>::iterator;
  using const_iterator = typename std::vector<T>::const_iterator;

 public:
  using Index = int;
  using SizeType = std::size_t;

  using iterator_category = std::forward_iterator_tag;
  using value_type = T;
  using difference_type = int;
  using pointer = value_type*;
  using reference = value_type&;

  MatrixEntry() = default;
  MatrixEntry(MatrixEntry<T>&&) = default;
  MatrixEntry(const MatrixEntry<T>&) = delete;
  MatrixEntry<T>& operator=(const MatrixEntry<T>&) = delete;
  MatrixEntry<T>& operator=(MatrixEntry<T>&&) = default;

  explicit MatrixEntry(const std::size_t size)
      : values_(size), nonzero_indices_(size) {}

  MatrixEntry(const std::vector<T>& values, const std::vector<Index>& indices)
      : values_(values), nonzero_indices_(indices) {
    if (values_.size() != nonzero_indices_.size()) {
      throw MismatchedDimensionsException();
    }
    check_entry_valid();
  }
  virtual ~MatrixEntry() = default;

  T operator[](const SizeType index) const {
    auto index_in_data =
        std::find(nonzero_indices_.begin(), nonzero_indices_.end(), index);
    if (index_in_data == nonzero_indices_.end()) {
      return T();
    } else {
#if NDEBUG
      return values_[static_cast<SizeType>(index_in_data -
                                           nonzero_indices_.begin())];
#else
      return values_.at(
          static_cast<SizeType>(index_in_data - nonzero_indices_.begin()));
#endif
    }
  }

  T lower_bound() const {
    return std::min_element(values_.begin(), values_.end());
  }

  T upper_bound() const {
    return std::max_element(values_.begin(), values_.end());
  }

  typename std::vector<T>::size_type num_nonzero() const {
    return values_.size();
  }

  const std::vector<T>& values() const { return values_; }

  std::vector<T>& values() { return values_; }

  const std::vector<Index>& nonzero_indices() const { return nonzero_indices_; }

  std::vector<Index>& nonzero_indices() { return nonzero_indices_; }

  iterator begin() { return values_.begin(); }

  const_iterator begin() const { return values_.begin(); }

  iterator end() { return values_.end(); }

  const_iterator end() const { return values_.end(); }

 private:
  std::vector<T> values_;
  std::vector<Index> nonzero_indices_;  // indices of nonzero entries

  void check_entry_valid() {
    // matrix entries are invalid if there are
    // duplicate nonzero indices present
    std::set<Index> index_set(nonzero_indices_.begin(), nonzero_indices_.end());
    if (nonzero_indices_.size() != index_set.size()) {
      throw InvalidMatrixEntryException();
    }
  }
};

// TODO: fix this for the case that left and right are non-equal permutations
// of each other, e.g. left == Entry {[1, 2] [0, 1]}, right == Entry {[2, 1],
// [0, 1]}. This example currently incorrectly evaluates as equal.
template <class T>
bool operator==(const MatrixEntry<T>& left, const MatrixEntry<T>& right) {
  return std::is_permutation(left.nonzero_indices().begin(),
                             left.nonzero_indices().end(),
                             right.nonzero_indices().begin()) &&
         std::is_permutation(left.values().begin(), left.values().end(),
                             right.values().begin());
}

template <typename T>
class Column : public MatrixEntry<T> {
 public:
  using Index = typename MatrixEntry<T>::Index;
  using SizeType = typename MatrixEntry<T>::SizeType;

 public:
  explicit Column(const std::size_t size) : MatrixEntry<T>(size) {}
  Column() = default;
  Column(const std::vector<T>& values, const std::vector<Index>& indices)
      : MatrixEntry<T>(values, indices) {}
  explicit Column(MatrixEntry<T>&& m) : MatrixEntry<T>(std::move(m)) {}
};

template <typename T>
class Row : public MatrixEntry<T> {
 public:
  using Index = typename MatrixEntry<T>::Index;
  using SizeType = typename MatrixEntry<T>::SizeType;

 public:
  explicit Row(const std::size_t size) : MatrixEntry<T>(size) {}
  Row() = default;
  Row(const std::vector<T>& values, const std::vector<Index>& indices)
      : MatrixEntry<T>(values, indices) {}
  explicit Row(MatrixEntry<T>&& m) : MatrixEntry<T>(std::move(m)) {}
};

template <typename T>
struct Constraint {
  static_assert(std::is_arithmetic<T>::value,
                "T must be arithmetic in order to be ordered");
  Constraint() : row(), lower_bound(), upper_bound() {}
  Constraint(Row<T>&& r, T lb, T ub)
      : row(std::move(r)), lower_bound(lb), upper_bound(ub) {}

  Row<T> row;
  T lower_bound;
  T upper_bound;
};

template <class T>
bool operator==(const Constraint<T>& left, const Constraint<T>& right) {
  return fabs(left.upper_bound - right.upper_bound) < DOUBLE_TOLERANCE &&
         fabs(left.lower_bound - right.lower_bound) < DOUBLE_TOLERANCE &&
         left.row == right.row;
}

template <typename T>
struct Objective {
  Objective() = default;
  // We can leave variable_types empty, since
  // Gurobi assumes continuous variables by default.
  // This might be different in other solvers, so
  // we'll have to carefully check when adding support.
  explicit Objective(std::vector<T>&& vals) : values(std::move(vals)) {}
  std::vector<T> values;
};

template <class T>
bool operator==(const Objective<T>& left, const Objective<T>& right) {
  return left.values == right.values;
}

class Variable {
 public:
  Variable() = default;
  Variable(double lb, double ub) : lower_bound_(lb), upper_bound_(ub) 
  {
    if (lb > ub) {
      throw InvalidVariableBoundsException();
    }
  }

  double upper() const { return upper_bound_; }
  double lower() const { return lower_bound_; }

 private:
  double lower_bound_ = 0.0;
  double upper_bound_ = LPINT_INFINITY;
};

inline bool operator==(const Variable& left, const Variable& right) {
  return left.lower() == right.lower() && left.upper() == right.upper();
}

inline std::ostream& operator<<(std::ostream& os, const Variable& var) {
  os << "Variable { lower = " << var.lower() << ", upper = " << var.upper()
     << " }";
  return os;
}

template <typename T>
struct Solution {
  std::vector<T> primal;
  std::vector<T> dual;
  T objective_value;
};

// LCOV_EXCL_START
template <typename T>
inline std::ostream& operator<<(std::ostream& os,
                                const Constraint<T>& constraint) {
  os << constraint.lower_bound << " <= " << constraint.row
     << " <= " << constraint.upper_bound;
  return os;
}

template <typename T>
inline std::ostream& operator<<(std::ostream& os, const Objective<T>& obj) {
  if (obj.values.size() == 0) {
    os << "Objective {}";
    return os;
  }
  os << "Objective {";
  std::size_t n = obj.values.size();
  for (std::size_t i = 0; i < n; i++) {
    os << obj.values[i] << ", ";
  }
  os << "\b\b}";
  return os;
}

template <typename T>
inline std::ostream& operator<<(std::ostream& os, const MatrixEntry<T>& row) {
  if (row.num_nonzero() == 0) {
    os << "Entry {[] []}";
    return os;
  }
  os << "Entry {[";
  for (const auto& val : row.values()) {
    os << val << " ";
  }
  os << "\b] [";
  for (const auto& val : row.nonzero_indices()) {
    os << val << " ";
  }
  os << "\b]}";
  return os;
}
// LCOV_EXCL_STOP

}  // namespace lpint

#endif  // LPINTERFACE_DATA_OBJECTS_H