problem.hpp

// Copyright (c) Sleipnir contributors
 
#pragma once
 
#include <algorithm>
#include <concepts>
#include <functional>
#include <iterator>
#include <optional>
#include <ranges>
#include <utility>
 
#include <Eigen/Core>
 
#include "sleipnir/autodiff/expression_type.hpp"
#include "sleipnir/autodiff/variable.hpp"
#include "sleipnir/autodiff/variable_matrix.hpp"
#include "sleipnir/optimization/solver/exit_status.hpp"
#include "sleipnir/optimization/solver/iteration_info.hpp"
#include "sleipnir/optimization/solver/options.hpp"
#include "sleipnir/util/small_vector.hpp"
#include "sleipnir/util/symbol_exports.hpp"
 
namespace slp {
 
class SLEIPNIR_DLLEXPORT Problem {
 public:
  Problem() noexcept = default;
 
  [[nodiscard]]
  Variable decision_variable() {
    m_decision_variables.emplace_back();
    return m_decision_variables.back();
  }
 
  [[nodiscard]]
  VariableMatrix decision_variable(int rows, int cols = 1) {
    m_decision_variables.reserve(m_decision_variables.size() + rows * cols);
 
    VariableMatrix vars{rows, cols};
 
    for (int row = 0; row < rows; ++row) {
      for (int col = 0; col < cols; ++col) {
        m_decision_variables.emplace_back();
        vars[row, col] = m_decision_variables.back();
      }
    }
 
    return vars;
  }
 
  [[nodiscard]]
  VariableMatrix symmetric_decision_variable(int rows) {
    // We only need to store the lower triangle of an n x n symmetric matrix;
    // the other elements are duplicates. The lower triangle has (n² + n)/2
    // elements.
    //
    //   n
    //   Σ k = (n² + n)/2
    //  k=1
    m_decision_variables.reserve(m_decision_variables.size() +
                                 (rows * rows + rows) / 2);
 
    VariableMatrix vars{rows, rows};
 
    for (int row = 0; row < rows; ++row) {
      for (int col = 0; col <= row; ++col) {
        m_decision_variables.emplace_back();
        vars[row, col] = m_decision_variables.back();
        vars[col, row] = m_decision_variables.back();
      }
    }
 
    return vars;
  }
 
  void minimize(const Variable& cost) { m_f = cost; }
 
  void minimize(Variable&& cost) { m_f = std::move(cost); }
 
  void maximize(const Variable& objective) {
    // Maximizing a cost function is the same as minimizing its negative
    m_f = -objective;
  }
 
  void maximize(Variable&& objective) {
    // Maximizing a cost function is the same as minimizing its negative
    m_f = -std::move(objective);
  }
 
  void subject_to(const EqualityConstraints& constraint) {
    m_equality_constraints.reserve(m_equality_constraints.size() +
                                   constraint.constraints.size());
    std::ranges::copy(constraint.constraints,
                      std::back_inserter(m_equality_constraints));
  }
 
  void subject_to(EqualityConstraints&& constraint) {
    m_equality_constraints.reserve(m_equality_constraints.size() +
                                   constraint.constraints.size());
    std::ranges::copy(constraint.constraints,
                      std::back_inserter(m_equality_constraints));
  }
 
  void subject_to(const InequalityConstraints& constraint) {
    m_inequality_constraints.reserve(m_inequality_constraints.size() +
                                     constraint.constraints.size());
    std::ranges::copy(constraint.constraints,
                      std::back_inserter(m_inequality_constraints));
  }
 
  void subject_to(InequalityConstraints&& constraint) {
    m_inequality_constraints.reserve(m_inequality_constraints.size() +
                                     constraint.constraints.size());
    std::ranges::copy(constraint.constraints,
                      std::back_inserter(m_inequality_constraints));
  }
 
  ExpressionType cost_function_type() const {
    if (m_f) {
      return m_f.value().type();
    } else {
      return ExpressionType::NONE;
    }
  }
 
  ExpressionType equality_constraint_type() const {
    if (!m_equality_constraints.empty()) {
      return std::ranges::max(m_equality_constraints, {}, &Variable::type)
          .type();
    } else {
      return ExpressionType::NONE;
    }
  }
 
  ExpressionType inequality_constraint_type() const {
    if (!m_inequality_constraints.empty()) {
      return std::ranges::max(m_inequality_constraints, {}, &Variable::type)
          .type();
    } else {
      return ExpressionType::NONE;
    }
  }
 
  ExitStatus solve(const Options& options = Options{},
                   [[maybe_unused]] bool spy = false);
 
  template <typename F>
    requires requires(F callback, const IterationInfo& info) {
      { callback(info) } -> std::same_as<void>;
    }
  void add_callback(F&& callback) {
    m_iteration_callbacks.emplace_back(
        [=, callback = std::forward<F>(callback)](const IterationInfo& info) {
          callback(info);
          return false;
        });
  }
 
  template <typename F>
    requires requires(F callback, const IterationInfo& info) {
      { callback(info) } -> std::same_as<bool>;
    }
  void add_callback(F&& callback) {
    m_iteration_callbacks.emplace_back(std::forward<F>(callback));
  }
 
  void clear_callbacks() { m_iteration_callbacks.clear(); }
 
 private:
  // The list of decision variables, which are the root of the problem's
  // expression tree
  small_vector<Variable> m_decision_variables;
 
  // The cost function: f(x)
  std::optional<Variable> m_f;
 
  // The list of equality constraints: cₑ(x) = 0
  small_vector<Variable> m_equality_constraints;
 
  // The list of inequality constraints: cᵢ(x) ≥ 0
  small_vector<Variable> m_inequality_constraints;
 
  // The iteration callbacks
  small_vector<std::function<bool(const IterationInfo& info)>>
      m_iteration_callbacks;
 
  void print_exit_conditions([[maybe_unused]] const Options& options);
  void print_problem_analysis();
};
 
}  // namespace slp