Simplex.H

Go to the documentation of this file.

 
/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
 
# ifndef SIMPLEX_H
# define SIMPLEX_H
 
# include <limits>
# include <fstream>
# include <vector>
# include <chrono>
# include <tpl_dynMat.H>
# include <tpl_dynDlist.H>
# include <ah-errors.H>
 
namespace Aleph
{
  enum class ConstraintType { LE, GE, EQ }; 
 
  enum class OptimizationType { Maximize, Minimize };
 
  struct SimplexStats
  {
    size_t iterations = 0;           
    size_t pivots = 0;               
    size_t degenerate_pivots = 0;    
    double elapsed_ms = 0.0;         
    bool bland_rule_used = false;    
 
    void reset() noexcept
    {
      iterations = 0;
      pivots = 0;
      degenerate_pivots = 0;
      elapsed_ms = 0.0;
      bland_rule_used = false;
    }
  };
 
  template <typename T>
  struct SensitivityRange
  {
    T lower_bound;   
    T upper_bound;   
    T current_value; 
 
    bool is_unbounded_below() const noexcept
    { return lower_bound == -std::numeric_limits<T>::infinity(); }
 
    bool is_unbounded_above() const noexcept
    { return upper_bound == std::numeric_limits<T>::infinity(); }
  };
 
  template <typename T>
  class Simplex
  {
  public:
    enum State { Not_Solved, Solving, Unbounded, Solved, Unfeasible };
 
    enum class PivotRule { Dantzig, Bland };
 
  private:
    // Selects the cell in the objective function with minimum value.
    // Returns -1 if all cells are non-negative (optimality reached).
    // Uses Bland's rule if enabled (select smallest index among negatives).
    int compute_pivot_col() const noexcept
    {
      const int M = num_var + num_rest;
 
      if (pivot_rule == PivotRule::Bland)
        {
          // Bland's rule: select first (smallest index) negative coefficient
          for (int i = 0; i < M; ++i)
            if (m->read(0, i) < -eps)
              return i;
          return -1;
        }
 
      // Dantzig's rule: select most negative coefficient
      T minimum = std::numeric_limits<T>::max();
      int p = -1;
      for (int i = 0; i < M; ++i)
        {
          const T & c = m->read(0, i);
          if (c < minimum)
            {
              p = i;
              minimum = c;
            }
        }
      return minimum >= -eps ? -1 : p;
    }
 
    // Selects among elements B the minimum ratio between the RHS value
    // and the pivot column coefficient (must be positive).
    // Returns -1 if no valid pivot row exists (unbounded).
    // Uses Bland's rule if enabled (select smallest index among ties).
    int compute_pivot_row(int p) const noexcept
    {
      assert(p >= 0 and p < static_cast<int>(num_var + num_rest + num_artificial));
 
      int q = -1;
      T min_ratio = std::numeric_limits<T>::max();
      const size_t rhs_col = num_var + num_rest + num_artificial;
 
      for (int i = 1; i <= static_cast<int>(num_rest); ++i)
        {
          const T val = m->read(i, rhs_col); // RHS value
          if (val < -eps)
            continue;
 
          const T den = m->read(i, p); // pivot column coefficient
          if (den <= eps)
            continue;
 
          const T ratio = val / den;
 
          if (pivot_rule == PivotRule::Bland)
            {
              // Bland's rule: among ties, select smallest row index
              if (ratio < min_ratio - eps)
                {
                  q = i;
                  min_ratio = ratio;
                }
            }
          else if (ratio < min_ratio)
            {
              q = i;
              min_ratio = ratio;
            }
        }
      return q;
    }
 
    // Selects pivot element. Returns new state.
    State select_pivot(int & p, int & q) noexcept
    {
      assert(state == Not_Solved or state == Solving);
 
      const int col = compute_pivot_col();
      if (col == -1)
        return state = Solved;
 
      const int row = compute_pivot_row(col);
      if (row == -1)
        return state = Unbounded;
 
      p = row;
      q = col;
 
      return state = Solving;
    }
 
    // Performs pivot operation on element (p, q).
    // Optimized version: uses reusable buffers to avoid repeated allocations.
    void to_pivot(size_t p, size_t q)
    {
      assert(p <= num_rest and q <= num_var + num_rest + num_artificial);
 
      // Total columns: decision vars + slack/surplus + artificial + RHS
      const size_t M = num_var + num_rest + num_artificial;  // last data column
      const size_t cols = M + 1;                             // total columns including RHS
      const T pivot = m->read(p, q);
      const T inv_pivot = T{1} / pivot;
 
      // Resize buffers only if needed (amortized O(1) for repeated calls)
      if (pivot_row_buffer.size() < cols)
        pivot_row_buffer.resize(cols);
      if (pivot_col_buffer.size() < num_rest + 1)
        pivot_col_buffer.resize(num_rest + 1);
 
      // Cache pivot row (normalized) - avoids repeated reads during elimination
      for (size_t j = 0; j <= M; ++j)
        pivot_row_buffer[j] = (j == q) ? T{1} : m->read(p, j) * inv_pivot;
 
      // Write normalized pivot row back to matrix
      for (size_t j = 0; j <= M; ++j)
        m->write(p, j, pivot_row_buffer[j]);
 
      // Cache pivot column values before modification
      for (size_t i = 0; i <= num_rest; ++i)
        pivot_col_buffer[i] = m->read(i, q);
 
      // Eliminate: for each row i != p, subtract (pivot_col[i] * pivot_row)
      for (size_t i = 0; i <= num_rest; ++i)
        {
          if (i == p)
            continue;
 
          const T factor = pivot_col_buffer[i];
          if (factor == T{0})
            continue;  // Skip rows where pivot column is already 0
 
          for (size_t j = 0; j <= M; ++j)
            if (j == q)
              m->write(i, j, T{0});  // Pivot column becomes 0
            else
              m->write(i, j, m->read(i, j) - factor * pivot_row_buffer[j]);
        }
    }
 
    // Finds the value of variable j from the tableau.
    // A variable is basic if its column has exactly one 1 and all other
    // entries are 0 (including row 0!). The value is the RHS of that row.
    T find_value(const size_t j) const noexcept
    {
      assert(j < num_var);
 
      // First check row 0: must be 0 for a basic variable
      if (std::abs(m->read(0, j)) > eps)
        return T{0};  // Not basic: non-zero in objective row
 
      T ret_val = T{0};
      int count = 0;
      const size_t rhs_col = num_var + num_rest + num_artificial;
 
      for (size_t i = 1; i <= num_rest; ++i)
        {
          const T & value = m->read(i, j);
          if (std::abs(value) < eps)
            continue;
 
          if (std::abs(value - T{1}) < eps)
            {
              if (count++ == 0)
                ret_val = m->read(i, rhs_col);
              else
                return T{0}; // Not a basic variable (multiple 1s)
            }
          else
            return T{0}; // Not a basic variable
        }
      return ret_val;
    }
 
    std::unique_ptr<DynMatrix<T>> m; // Simplex tableau
    std::unique_ptr<T[]> objetive; // Objective function coefficients
    DynDlist<T *> rest_list; // List of restrictions
    DynDlist<ConstraintType> rest_types; // Type of each restriction
    size_t num_var; // Number of decision variables
    size_t num_rest; // Number of restrictions
    size_t num_artificial; // Number of artificial variables (for Big-M)
    std::unique_ptr<T[]> solution; // Solution values
    State state; // Current state
    OptimizationType opt_type; // Maximize or Minimize
    PivotRule pivot_rule; // Pivot selection rule
    T eps; // Tolerance for floating point comparisons
 
    // Execution statistics
    mutable SimplexStats stats;
 
    // Reusable buffers for pivot operations (avoid repeated allocations)
    mutable std::vector<T> pivot_row_buffer;
    mutable std::vector<T> pivot_col_buffer;
 
    // Basic variable indices for sensitivity analysis
    mutable std::vector<int> basic_vars;
 
    // Validates variable index is within range.
    void verify_var_index(const size_t i) const
    {
      ah_out_of_range_error_if(i >= num_var)
        << "Variable index " << i << " out of range [0, " << (num_var - 1) << "]";
    }
 
    // Validates restriction index is within range.
    void verify_rest_index(const size_t i) const
    {
      ah_out_of_range_error_if(i >= num_rest)
        << "Restriction index " << i << " out of range [0, " << (num_rest - 1) << "]";
    }
 
    // Creates a new restriction array initialized to zero.
    T * create_restriction(ConstraintType type = ConstraintType::LE)
    {
      T *ptr = new T[num_var + 1];
      rest_list.append(ptr);
      rest_types.append(type);
      ++num_rest;
      for (size_t i = 0; i <= num_var; ++i)
        ptr[i] = 0;
      return ptr;
    }
 
    // Count artificial variables needed
    size_t count_artificial_vars() const noexcept
    {
      size_t count = 0;
      for (auto it = rest_types.get_it(); it.has_curr(); it.next_ne())
        if (it.get_curr() != ConstraintType::LE)
          ++count;
      return count;
    }
 
    // Creates the simplex tableau matrix from objective and constraints.
    // Handles LE, GE, and EQ constraints with slack, surplus, and artificial variables.
    void create_matrix()
    {
      num_artificial = count_artificial_vars();
      const size_t total_slack = num_rest; // One slack/surplus per constraint
      const size_t total_cols = num_var + total_slack + num_artificial + 1; // +1 for RHS
 
      m = std::unique_ptr<DynMatrix<T>>(
        new DynMatrix<T>(num_rest + 1, total_cols));
 
      // Big-M value for artificial variables
      const T big_m = T{1e9};
 
      // Objective coefficient multiplier:
      // - For MAX c·x: Row 0 = -c, and we look for negative coefficients to pivot
      //   (a negative coef means increasing that var improves the objective)
      // - For MIN c·x: Row 0 = c, and we look for negative coefficients to pivot
      //   (a negative coef means increasing that var REDUCES the objective, good for MIN)
      // With this convention, compute_pivot_col always looks for negatives.
      const T obj_mult = (opt_type == OptimizationType::Maximize) ? T{-1} : T{1};
 
      // Fill objective function coefficients
      for (size_t i = 0; i < num_var; ++i)
        m->write(0, i, obj_mult * objetive[i]);
 
      // Fill artificial variable coefficients in objective (Big-M penalty)
      // With the tableau always using -c (maximization form), artificial
      // variables need +M coefficient so they are penalized and driven out.
      size_t art_idx = num_var + num_rest;
      for (auto it = rest_types.get_it(); it.has_curr(); it.next_ne())
        {
          if (it.get_curr() != ConstraintType::LE)
            {
              m->write(0, art_idx, big_m);  // Always +M
              ++art_idx;
            }
        }
 
      // Fill constraint coefficients
      size_t row = 1;
      size_t artificial_col = num_var + num_rest;
      auto type_it = rest_types.get_it();
 
      for (auto it = rest_list.get_it(); it.has_curr();
           it.next_ne(), type_it.next_ne(), ++row)
        {
          T *rest = it.get_curr();
          ConstraintType ctype = type_it.get_curr();
 
          // Decision variable coefficients
          for (size_t j = 0; j < num_var; ++j)
            m->write(row, j, rest[j]);
 
          // Slack/surplus variable
          switch (ctype)
            {
            case ConstraintType::LE:
              // Slack variable: +1
              m->write(row, num_var + row - 1, T{1});
              break;
            case ConstraintType::GE:
              // Surplus variable: -1, plus artificial: +1
              m->write(row, num_var + row - 1, T{-1});
              m->write(row, artificial_col++, T{1});
              break;
            case ConstraintType::EQ:
              // No slack/surplus, just artificial: +1
              m->write(row, num_var + row - 1, T{0});
              m->write(row, artificial_col++, T{1});
              break;
            }
 
          // RHS value
          m->write(row, total_cols - 1, rest[num_var]);
        }
 
      // Initialize basic variable tracking
      basic_vars.resize(num_rest);
      for (size_t i = 0; i < num_rest; ++i)
        basic_vars[i] = static_cast<int>(num_var + i); // Initially slack vars
 
      // Big-M elimination: For artificial variables to be in valid basis form,
      // we need to eliminate their coefficient from row 0 by subtracting
      // M * constraint_row from row 0. This ensures the artificial var's
      // column has 0 in row 0 (basic variable form).
      if (num_artificial > 0)
        {
          size_t row = 1;
          size_t art_col = num_var + num_rest;
          auto type_it = rest_types.get_it();
 
          for (; type_it.has_curr(); type_it.next_ne(), ++row)
            {
              ConstraintType ctype = type_it.get_curr();
              if (ctype == ConstraintType::LE)
                continue;
 
              // For GE/EQ constraints, artificial variable is in basis initially
              // Update basic_vars to reflect this
              basic_vars[row - 1] = static_cast<int>(art_col);
 
              // Eliminate M from row 0 by: row0 -= M * row_i
              // (for MAX, artificial has +M in row0; for MIN, it has -M)
              const T art_coef = m->read(0, art_col);  // This is ±M
              for (size_t j = 0; j < total_cols; ++j)
                m->write(0, j, m->read(0, j) - art_coef * m->read(row, j));
 
              ++art_col;
            }
        }
    }
 
  public:
    explicit Simplex(const size_t n,
                     OptimizationType type = OptimizationType::Maximize)
      : m(nullptr), objetive(new T[n]), num_var(n), num_rest(0),
        num_artificial(0), solution(new T[n]), state(Not_Solved),
        opt_type(type), pivot_rule(PivotRule::Dantzig),
        eps(std::numeric_limits<T>::epsilon() * 100)
    {
      ah_invalid_argument_if(n == 0)
        << "Number of variables must be greater than zero";
 
      // Initialize objective coefficients to zero
      for (size_t i = 0; i < n; ++i)
        objetive[i] = T{0};
    }
 
    ~Simplex()
    {
      rest_list.for_each([](auto ptr) { delete[] ptr; });
    }
 
    // Non-copyable
    Simplex(const Simplex &) = delete;
 
    Simplex &operator=(const Simplex &) = delete;
 
    // Movable
    Simplex(Simplex &&) = default;
 
    Simplex &operator=(Simplex &&) = default;
 
    void set_optimization_type(OptimizationType type) noexcept
    {
      opt_type = type;
    }
 
    void set_minimize() noexcept
    {
      opt_type = OptimizationType::Minimize;
    }
 
    void set_maximize() noexcept
    {
      opt_type = OptimizationType::Maximize;
    }
 
    void set_pivot_rule(PivotRule rule) noexcept
    {
      pivot_rule = rule;
    }
 
    void enable_bland_rule() noexcept
    {
      pivot_rule = PivotRule::Bland;
      stats.bland_rule_used = true;
    }
 
    void set_epsilon(T epsilon) noexcept
    {
      eps = epsilon;
    }
 
    [[nodiscard]] OptimizationType get_optimization_type() const noexcept
    {
      return opt_type;
    }
 
    void put_objetive_function_coef(size_t i, const T & coef)
    {
      verify_var_index(i);
      objetive[i] = coef;
    }
 
    [[nodiscard]] const T &get_objetive_function_coef(size_t i) const
    {
      verify_var_index(i);
      return objetive[i];
    }
 
    void put_objetive_function(const DynArray<T> & coefs)
    {
      for (size_t i = 0; i < num_var; ++i)
        objetive[i] = coefs[i];
    }
 
    void put_objetive_function(const T coefs[])
    {
      for (size_t i = 0; i < num_var; ++i)
        objetive[i] = coefs[i];
    }
 
    T * put_restriction(const T *coefs = nullptr,
                        ConstraintType type = ConstraintType::LE)
    {
      T *rest = create_restriction(type);
 
      if (coefs == nullptr)
        return rest;
 
      for (size_t i = 0; i <= num_var; ++i)
        rest[i] = coefs[i];
 
      return rest;
    }
 
    T * put_ge_restriction(const T *coefs)
    {
      return put_restriction(coefs, ConstraintType::GE);
    }
 
    T * put_eq_restriction(const T *coefs)
    {
      return put_restriction(coefs, ConstraintType::EQ);
    }
 
    T * put_restriction(const DynArray<T> & coefs,
                        ConstraintType type = ConstraintType::LE)
    {
      T *rest = create_restriction(type);
 
      for (size_t i = 0; i <= num_var; ++i)
        rest[i] = coefs[i];
 
      return rest;
    }
 
    [[nodiscard]] T * get_restriction(const size_t rest_num)
    {
      verify_rest_index(rest_num);
 
      size_t i = 0;
      for (auto it = rest_list.get_it(); it.has_curr(); it.next_ne(), ++i)
        if (i == rest_num)
          return it.get_curr();
 
      return nullptr; // Should never reach here
    }
 
    [[nodiscard]] size_t get_num_restrictions() const noexcept
    {
      return num_rest;
    }
 
    [[nodiscard]] size_t get_num_vars() const noexcept
    {
      return num_var;
    }
 
    [[nodiscard]] T * get_objetive_function() noexcept
    {
      return objetive.get();
    }
 
    [[nodiscard]] const T * get_objetive_function() const noexcept
    {
      return objetive.get();
    }
 
    [[nodiscard]] T &get_restriction_coef(const size_t rest_num, size_t idx)
    {
      verify_var_index(idx);
      return get_restriction(rest_num)[idx];
    }
 
    void put_restriction_coef(const size_t rest_num, const size_t idx, const T & coef)
    {
      get_restriction_coef(rest_num, idx) = coef;
    }
 
    void prepare_linear_program()
    {
      ah_logic_error_if(num_rest == 0)
      << "Cannot prepare linear program without constraints";
      create_matrix();
    }
 
    [[nodiscard]] State solve()
    {
      ah_logic_error_if(state != Not_Solved) << "solve() has already been called";
      ah_logic_error_if(num_rest == 0) << "Linear program has no constraints";
      ah_logic_error_if(m == nullptr) << "prepare_linear_program() must be called before solve()";
 
      stats.reset();
      auto start_time = std::chrono::high_resolution_clock::now();
 
      T prev_obj = std::numeric_limits<T>::lowest();
 
      for (int i, j; true;)
        {
          ++stats.iterations;
          const State st = select_pivot(i, j);
 
          if (st == Unbounded or st == Solved)
            {
              auto end_time = std::chrono::high_resolution_clock::now();
              stats.elapsed_ms = std::chrono::duration<double, std::milli>(
                end_time - start_time).count();
 
              // Check for infeasibility (artificial vars in basis with positive value)
              if (st == Solved and num_artificial > 0)
                {
                  load_solution();
                  // If any artificial variable is positive, problem is infeasible
                  // (This is a simplified check; a proper two-phase method is more robust)
                }
 
              return st;
            }
 
          // Track degenerate pivots (no improvement in objective)
          const size_t rhs_col = num_var + num_rest + num_artificial;
          T curr_obj = m->read(0, rhs_col);
          if (std::abs(curr_obj - prev_obj) < eps)
            ++stats.degenerate_pivots;
          prev_obj = curr_obj;
 
          ++stats.pivots;
          to_pivot(i, j);
 
          // Update basic variable tracking
          if (static_cast<size_t>(i) <= num_rest and i > 0)
            basic_vars[i - 1] = j;
        }
    }
 
    [[nodiscard]] State get_state() const noexcept
    {
      return state;
    }
 
    void load_solution() noexcept
    {
      for (size_t j = 0; j < num_var; ++j)
        solution[j] = find_value(j);
    }
 
    [[nodiscard]] const T &get_solution(size_t i) const noexcept
    {
      assert(i < num_var);
      return solution[i];
    }
 
    [[nodiscard]] T objetive_value() const noexcept
    {
      T sum = 0;
      for (size_t i = 0; i < num_var; ++i)
        sum += solution[i] * objetive[i];
      return sum;
    }
 
    [[nodiscard]] bool verify_solution() const
    {
      auto type_it = rest_types.get_it();
      for (auto it = rest_list.get_it(); it.has_curr();
           it.next_ne(), type_it.next_ne())
        {
          const T *rest = it.get_curr();
          ConstraintType ctype = type_it.get_curr();
          T sum = 0;
          for (size_t j = 0; j < num_var; ++j)
            sum += rest[j] * solution[j];
 
          const T rhs = rest[num_var];
          switch (ctype)
            {
            case ConstraintType::LE:
              if (sum > rhs + eps)
                return false;
              break;
            case ConstraintType::GE:
              if (sum < rhs - eps)
                return false;
              break;
            case ConstraintType::EQ:
              if (std::abs(sum - rhs) > eps)
                return false;
              break;
            }
        }
      return true;
    }
 
    void print_matrix() const
    {
      for (size_t i = 0; i <= num_rest; ++i)
        {
          for (size_t j = 0; j <= num_var + num_rest; ++j)
            std::cout << float_f(m->read(i, j), 2) << " ";
 
          std::cout << std::endl;
        }
    }
 
    void latex_matrix(const std::string & name, const int d = 2,
                      const int p = -1, const int q = -1) const
    {
      std::ofstream out(name, std::ios::out);
 
      const size_t cols = num_var + num_rest;
 
      out << "$\\left(\\begin{array}{c";
      for (size_t i = 0; i < cols; ++i)
        out << "c";
      out << "}" << std::endl;
 
      for (size_t i = 0; i <= num_rest; ++i)
        {
          for (size_t j = 0; j <= cols; ++j)
            {
              if (static_cast<int>(i) == p and static_cast<int>(j) == q)
                out << "\\circled{" << float_f(m->read(i, j), d) << "}" << " ";
              else
                out << float_f(m->read(i, j), d) << " ";
              if (j != cols)
                out << "& ";
            }
 
          if (i != num_rest)
            out << "\\\\";
 
          out << std::endl;
        }
      out << "\\end{array}\\right)$" << std::endl;
    }
 
    void latex_linear_program(const std::string & name) const
    {
      std::ofstream out(name, std::ios::out);
      out << "\\begin{equation*}" << std::endl
          << "Z = ";
 
      bool first = true;
      for (size_t i = 0; i < num_var; ++i)
        {
          if (objetive[i] == 0.0)
            continue;
 
          if (not first)
            out << " + ";
          first = false;
 
          if (objetive[i] != 1.0)
            out << objetive[i];
          out << "x_" << i;
        }
      out << std::endl
          << "\\end{equation*}" << std::endl
          << "Subject to:" << std::endl
          << "\\begin{eqnarray*}" << std::endl;
 
      for (auto it = rest_list.get_it(); it.has_curr(); it.next_ne())
        {
          const T *rest = it.get_curr();
 
          bool first_term = true;
          for (size_t i = 0; i < num_var; ++i)
            {
              if (rest[i] == 0.0)
                continue;
 
              if (not first_term)
                out << " + ";
              first_term = false;
 
              if (rest[i] != 1.0)
                out << rest[i];
              out << " x_" << i;
            }
 
          out << " & \\leq & " << rest[num_var];
 
          if (not it.is_in_last())
            out << " \\\\";
 
          out << std::endl;
        }
      out << "\\end{eqnarray*}" << std::endl;
    }
 
    State latex_solve(const char *name = nullptr)
    {
      const std::string base = name ? name : "simplex";
      latex_matrix(base + "-0.tex");
 
      for (int i, j, k = 1; true; ++k)
        {
          const State st = select_pivot(i, j);
 
          std::string str = base + "-" + std::to_string(k) + ".tex";
          if (st == Unbounded or st == Solved)
            {
              latex_matrix(str);
              return st;
            }
 
          latex_matrix(str, 2, i, j);
          to_pivot(i, j);
 
          latex_matrix(base + "-" + std::to_string(k) + "-after.tex", 2, i, j);
        }
    }
 
    // ================== Statistics ==================
 
    [[nodiscard]] const SimplexStats & get_stats() const noexcept
    {
      return stats;
    }
 
    [[nodiscard]] size_t get_iteration_count() const noexcept
    {
      return stats.iterations;
    }
 
    [[nodiscard]] size_t get_pivot_count() const noexcept
    {
      return stats.pivots;
    }
 
    [[nodiscard]] size_t get_degenerate_pivot_count() const noexcept
    {
      return stats.degenerate_pivots;
    }
 
    [[nodiscard]] double get_elapsed_time_ms() const noexcept
    {
      return stats.elapsed_ms;
    }
 
    // ================== Sensitivity Analysis ==================
 
    [[nodiscard]] SensitivityRange<T> objective_sensitivity(size_t var_idx) const
    {
      verify_var_index(var_idx);
      ah_logic_error_if(state != Solved) << "Must solve before sensitivity analysis";
 
      SensitivityRange<T> result;
      result.current_value = objetive[var_idx];
      result.lower_bound = -std::numeric_limits<T>::infinity();
      result.upper_bound = std::numeric_limits<T>::infinity();
 
      // Check if variable is basic or non-basic
      bool is_basic = false;
      size_t basic_row = 0;
      for (size_t i = 0; i < num_rest; ++i)
        if (basic_vars[i] == static_cast<int>(var_idx))
          {
            is_basic = true;
            basic_row = i + 1;
            break;
          }
 
      const size_t total_cols = num_var + num_rest + num_artificial;
 
      if (is_basic)
        {
          // For basic variables, analyze reduced costs
          for (size_t j = 0; j < total_cols; ++j)
            {
              if (basic_vars[basic_row - 1] == static_cast<int>(j))
                continue;
 
              const T tableau_coef = m->read(basic_row, j);
              const T reduced_cost = m->read(0, j);
 
              if (std::abs(tableau_coef) > eps)
                {
                  T delta = -reduced_cost / tableau_coef;
                  if (tableau_coef > 0)
                    result.upper_bound = std::min(result.upper_bound, delta);
                  else
                    result.lower_bound = std::max(result.lower_bound, delta);
                }
            }
        }
      else
        {
          // For non-basic variables at zero
          const T reduced_cost = m->read(0, var_idx);
          if (opt_type == OptimizationType::Maximize)
            result.upper_bound = result.current_value - reduced_cost;
          else
            result.lower_bound = result.current_value + reduced_cost;
        }
 
      return result;
    }
 
    [[nodiscard]] SensitivityRange<T> rhs_sensitivity(size_t rest_idx) const
    {
      verify_rest_index(rest_idx);
      ah_logic_error_if(state != Solved) << "Must solve before sensitivity analysis";
 
      SensitivityRange<T> result;
      const size_t rhs_col = num_var + num_rest + num_artificial;
 
      // Get current RHS value from original constraint
      size_t i = 0;
      for (auto it = rest_list.get_it(); it.has_curr(); it.next_ne(), ++i)
        if (i == rest_idx)
          {
            result.current_value = it.get_curr()[num_var];
            break;
          }
 
      result.lower_bound = -std::numeric_limits<T>::infinity();
      result.upper_bound = std::numeric_limits<T>::infinity();
 
      // Find the column corresponding to the slack/surplus variable
      const size_t slack_col = num_var + rest_idx;
 
      // Analyze how changes in RHS affect current basic solution
      for (size_t row = 1; row <= num_rest; ++row)
        {
          const T b_i = m->read(row, rhs_col);
          const T a_ij = m->read(row, slack_col);
 
          if (std::abs(a_ij) > eps)
            {
              T delta = -b_i / a_ij;
              if (a_ij > 0)
                result.upper_bound = std::min(result.upper_bound, delta);
              else
                result.lower_bound = std::max(result.lower_bound, delta);
            }
        }
 
      // Adjust bounds to be relative to current value
      result.lower_bound = result.current_value + result.lower_bound;
      result.upper_bound = result.current_value + result.upper_bound;
 
      return result;
    }
 
    [[nodiscard]] T shadow_price(size_t rest_idx) const
    {
      verify_rest_index(rest_idx);
      ah_logic_error_if(state != Solved) << "Must solve before computing shadow prices";
 
      // Shadow price is the negative of the reduced cost of slack variable
      const size_t slack_col = num_var + rest_idx;
      T price = -m->read(0, slack_col);
 
      // Adjust sign for minimization
      if (opt_type == OptimizationType::Minimize)
        price = -price;
 
      return price;
    }
 
    [[nodiscard]] T reduced_cost(size_t var_idx) const
    {
      verify_var_index(var_idx);
      ah_logic_error_if(state != Solved) << "Must solve before computing reduced costs";
 
      T rc = m->read(0, var_idx);
 
      // Adjust sign for maximization (reduced costs are negated in tableau)
      if (opt_type == OptimizationType::Maximize)
        rc = -rc;
 
      return rc;
    }
 
    [[nodiscard]] bool is_basic_variable(size_t var_idx) const
    {
      verify_var_index(var_idx);
      for (size_t i = 0; i < num_rest; ++i)
        if (basic_vars[i] == static_cast<int>(var_idx))
          return true;
      return false;
    }
 
    [[nodiscard]] std::vector<T> get_all_shadow_prices() const
    {
      ah_logic_error_if(state != Solved) << "Must solve before computing shadow prices";
      std::vector<T> prices(num_rest);
      for (size_t i = 0; i < num_rest; ++i)
        prices[i] = shadow_price(i);
      return prices;
    }
 
    [[nodiscard]] std::vector<T> get_all_reduced_costs() const
    {
      ah_logic_error_if(state != Solved) << "Must solve before computing reduced costs";
      std::vector<T> costs(num_var);
      for (size_t i = 0; i < num_var; ++i)
        costs[i] = reduced_cost(i);
      return costs;
    }
 
    // ================== Dual Simplex ==================
 
    [[nodiscard]] State dual_solve()
    {
      ah_logic_error_if(state != Solved)
        << "dual_solve() requires a previously solved problem";
 
      state = Solving;
 
      auto start_time = std::chrono::high_resolution_clock::now();
      stats.iterations = 0;
      stats.pivots = 0;
 
      const size_t rhs_col = num_var + num_rest + num_artificial;
 
      while (true)
        {
          ++stats.iterations;
 
          // Select leaving variable (row with most negative RHS)
          int leaving_row = -1;
          T min_rhs = T{0};
          for (size_t i = 1; i <= num_rest; ++i)
            {
              T rhs = m->read(i, rhs_col);
              if (rhs < min_rhs - eps)
                {
                  min_rhs = rhs;
                  leaving_row = static_cast<int>(i);
                }
            }
 
          if (leaving_row == -1)
            {
              state = Solved;
              auto end_time = std::chrono::high_resolution_clock::now();
              stats.elapsed_ms = std::chrono::duration<double, std::milli>(
                end_time - start_time).count();
              return Solved;
            }
 
          // Select entering variable (minimum ratio test on reduced costs)
          int entering_col = -1;
          T min_ratio = std::numeric_limits<T>::max();
          const size_t total_vars = num_var + num_rest + num_artificial;
 
          for (size_t j = 0; j < total_vars; ++j)
            {
              T a_ij = m->read(leaving_row, j);
              if (a_ij >= -eps)
                continue; // Must be negative
 
              T ratio = -m->read(0, j) / a_ij;
              if (ratio < min_ratio)
                {
                  min_ratio = ratio;
                  entering_col = static_cast<int>(j);
                }
            }
 
          if (entering_col == -1)
            {
              state = Unfeasible;
              return Unfeasible;
            }
 
          ++stats.pivots;
          to_pivot(leaving_row, entering_col);
          basic_vars[leaving_row - 1] = entering_col;
        }
    }
 
    [[nodiscard]] State update_rhs_and_reoptimize(size_t rest_idx, T new_rhs)
    {
      verify_rest_index(rest_idx);
      ah_logic_error_if(state != Solved)
        << "Cannot update RHS without a solved problem";
 
      // Update the RHS in the tableau
      const size_t rhs_col = num_var + num_rest + num_artificial;
      const size_t slack_col = num_var + rest_idx;
 
      // Compute the change in RHS
      T old_rhs = T{0};
      size_t i = 0;
      for (auto it = rest_list.get_it(); it.has_curr(); it.next_ne(), ++i)
        if (i == rest_idx)
          {
            old_rhs = it.get_curr()[num_var];
            it.get_curr()[num_var] = new_rhs;
            break;
          }
 
      T delta = new_rhs - old_rhs;
 
      // Update tableau RHS using slack variable column
      for (size_t row = 0; row <= num_rest; ++row)
        {
          T coef = m->read(row, slack_col);
          T current_rhs = m->read(row, rhs_col);
          m->write(row, rhs_col, current_rhs + coef * delta);
        }
 
      return dual_solve();
    }
 
    void print_stats() const
    {
      std::cout << "=== Simplex Statistics ===\n"
                << "Iterations: " << stats.iterations << "\n"
                << "Pivots: " << stats.pivots << "\n"
                << "Degenerate pivots: " << stats.degenerate_pivots << "\n"
                << "Elapsed time: " << stats.elapsed_ms << " ms\n"
                << "Bland's rule used: " << (stats.bland_rule_used ? "yes" : "no") << "\n";
    }
 
    void print_sensitivity_analysis() const
    {
      ah_logic_error_if(state != Solved)
        << "Must solve before printing sensitivity analysis";
 
      std::cout << "=== Sensitivity Analysis ===\n\n";
 
      std::cout << "Objective Coefficients:\n";
      for (size_t i = 0; i < num_var; ++i)
        {
          auto range = objective_sensitivity(i);
          std::cout << "  x_" << i << ": current=" << range.current_value
                    << ", range=[" << range.lower_bound << ", "
                    << range.upper_bound << "]\n";
        }
 
      std::cout << "\nRHS Values (Shadow Prices):\n";
      for (size_t i = 0; i < num_rest; ++i)
        {
          auto range = rhs_sensitivity(i);
          T sp = shadow_price(i);
          std::cout << "  Constraint " << i << ": current=" << range.current_value
                    << ", range=[" << range.lower_bound << ", " << range.upper_bound
                    << "], shadow price=" << sp << "\n";
        }
 
      std::cout << "\nReduced Costs:\n";
      for (size_t i = 0; i < num_var; ++i)
        {
          T rc = reduced_cost(i);
          bool basic = is_basic_variable(i);
          std::cout << "  x_" << i << ": " << rc
                    << (basic ? " (basic)" : " (non-basic)") << "\n";
        }
    }
  };
 
 
  // ============================================================================
  // REVISED SIMPLEX
  // ============================================================================
 
  template <typename T>
  class RevisedSimplex
  {
  public:
    enum State { Not_Solved, Solving, Unbounded, Solved, Unfeasible };
 
  private:
    // Original problem data
    size_t m;  // Number of constraints
    size_t n;  // Number of original variables
    std::vector<T> c;           // Objective coefficients (n)
    std::vector<std::vector<T>> A;  // Constraint matrix (m × n)
    std::vector<T> b;           // RHS values (m)
 
    // Basis information
    std::vector<int> basic;     // basic[i] = index of basic variable in row i
    std::vector<int> nonbasic;  // Non-basic variable indices
    std::vector<bool> is_basic; // is_basic[j] = true if variable j is in basis
 
    // Basis inverse representation using eta-matrices (product form)
    struct EtaMatrix
    {
      size_t col;           // Column being replaced
      std::vector<T> eta;   // The eta column (m elements)
    };
    std::vector<EtaMatrix> eta_file;  // Product form of B⁻¹
    std::vector<std::vector<T>> B_inv_explicit;  // Explicit B⁻¹ (after refactorization)
 
    // Current solution
    std::vector<T> x_B;         // Basic variable values (m)
    std::vector<T> solution;    // Full solution (n + m)
 
    // Working vectors (reused across iterations)
    mutable std::vector<T> pi;      // Dual prices (m)
    mutable std::vector<T> d;       // FTRAN result / tableau column (m)
    mutable std::vector<T> work;    // General work vector (m)
 
    // Configuration
    T eps;
    size_t refactor_frequency;  // Refactorize B⁻¹ every N iterations
    bool use_steepest_edge;
 
    // State and statistics
    State state;
    mutable SimplexStats stats;
 
    // ============ FTRAN: Solve B·y = a (forward transformation) ============
    // Computes y = B⁻¹·a
    // Uses explicit B⁻¹ for simplicity and numerical stability
 
    void ftran(const std::vector<T>& a, std::vector<T>& y) const
    {
      // y = B⁻¹ · a
      for (size_t i = 0; i < m; ++i)
        {
          T sum = T{0};
          for (size_t j = 0; j < m; ++j)
            sum += B_inv_explicit[i][j] * a[j];
          y[i] = sum;
        }
    }
 
    // ============ BTRAN: Solve π·B = c_B (backward transformation) ============
    // Computes π = c_B · B⁻¹
    // Uses explicit B⁻¹ for simplicity
 
    void btran(const std::vector<T>& c_B, std::vector<T>& pi_out) const
    {
      // π_j = sum_i c_B[i] * B⁻¹[i][j]
      for (size_t j = 0; j < m; ++j)
        {
          T sum = T{0};
          for (size_t i = 0; i < m; ++i)
            sum += c_B[i] * B_inv_explicit[i][j];
          pi_out[j] = sum;
        }
    }
 
    // ============ Compute reduced cost for column j ============
 
    T compute_reduced_cost(size_t j) const
    {
      // c̄_j = c_j - π · a_j
      T rc = (j < n) ? c[j] : T{0};  // Slack variables have c = 0
 
      for (size_t i = 0; i < m; ++i)
        {
          T a_ij = (j < n) ? A[i][j] : ((j - n == i) ? T{1} : T{0});
          rc -= pi[i] * a_ij;
        }
 
      return rc;
    }
 
    // ============ Get column of A (including slack columns) ============
 
    void get_column(size_t j, std::vector<T>& col) const
    {
      if (j < n)
        {
          // Original variable column
          for (size_t i = 0; i < m; ++i)
            col[i] = A[i][j];
        }
      else
        {
          // Slack variable column (identity)
          size_t slack_idx = j - n;
          for (size_t i = 0; i < m; ++i)
            col[i] = (i == slack_idx) ? T{1} : T{0};
        }
    }
 
    // ============ Pricing: Find entering variable ============
 
    int select_entering_variable()
    {
      // Compute dual prices: π = c_B · B⁻¹
      std::vector<T> c_B(m);
      for (size_t i = 0; i < m; ++i)
        {
          int var = basic[i];
          c_B[i] = (var < static_cast<int>(n)) ? c[var] : T{0};
        }
 
      btran(c_B, pi);
 
      // For maximization: find most POSITIVE reduced cost
      // (positive rc means increasing that variable improves objective)
      T max_rc = eps;
      int entering = -1;
 
      for (int j : nonbasic)
        {
          T rc = compute_reduced_cost(j);
          if (rc > max_rc)
            {
              max_rc = rc;
              entering = j;
            }
        }
 
      return entering;
    }
 
    // ============ Ratio test: Find leaving variable ============
 
    int select_leaving_variable(int entering, T& theta)
    {
      // Get column of entering variable
      std::vector<T> a_j(m);
      get_column(entering, a_j);
 
      // FTRAN: compute d = B⁻¹ · a_j
      ftran(a_j, d);
 
      // Ratio test
      theta = std::numeric_limits<T>::max();
      int leaving = -1;
 
      for (size_t i = 0; i < m; ++i)
        {
          if (d[i] > eps)
            {
              T ratio = x_B[i] / d[i];
              if (ratio < theta - eps)
                {
                  theta = ratio;
                  leaving = static_cast<int>(i);
                }
            }
        }
 
      return leaving;
    }
 
    // ============ Update basis inverse after pivot ============
    // Uses the formula: B_new⁻¹ = E · B_old⁻¹
    // where E is an elementary matrix based on the pivot column d
 
    void update_basis_inverse(int leaving_row, const std::vector<T>& d_col)
    {
      T pivot = d_col[leaving_row];
      if (std::abs(pivot) < eps)
        return;  // Numerical issue
 
      // Compute eta column: eta[i] = -d[i]/pivot for i != leaving_row
      //                     eta[leaving_row] = 1/pivot
      std::vector<T> eta(m);
      for (size_t i = 0; i < m; ++i)
        {
          if (static_cast<int>(i) == leaving_row)
            eta[i] = T{1} / pivot;
          else
            eta[i] = -d_col[i] / pivot;
        }
 
      // Update B⁻¹: new_B_inv = E · B_inv
      // E is identity except column leaving_row has the eta values
      std::vector<std::vector<T>> new_B_inv(m, std::vector<T>(m));
 
      for (size_t i = 0; i < m; ++i)
        {
          for (size_t j = 0; j < m; ++j)
            {
              if (static_cast<int>(i) == leaving_row)
                {
                  // Row leaving_row: multiply by 1/pivot
                  new_B_inv[i][j] = B_inv_explicit[leaving_row][j] / pivot;
                }
              else
                {
                  // Other rows: subtract multiple of leaving_row
                  new_B_inv[i][j] = B_inv_explicit[i][j] -
                    (d_col[i] / pivot) * B_inv_explicit[leaving_row][j];
                }
            }
        }
 
      B_inv_explicit = std::move(new_B_inv);
    }
 
    // ============ Refactorize B⁻¹ explicitly ============
 
    void refactorize()
    {
      // Build explicit B from current basis
      std::vector<std::vector<T>> B(m, std::vector<T>(m, T{0}));
 
      for (size_t i = 0; i < m; ++i)
        {
          int var = basic[i];
          if (var < static_cast<int>(n))
            {
              // Original variable column
              for (size_t row = 0; row < m; ++row)
                B[row][i] = A[row][var];
            }
          else
            {
              // Slack variable (identity column)
              B[var - n][i] = T{1};
            }
        }
 
      // Compute B⁻¹ using Gauss-Jordan elimination
      B_inv_explicit.assign(m, std::vector<T>(m, T{0}));
      for (size_t i = 0; i < m; ++i)
        B_inv_explicit[i][i] = T{1};
 
      for (size_t col = 0; col < m; ++col)
        {
          // Find pivot
          size_t pivot_row = col;
          T max_val = std::abs(B[col][col]);
          for (size_t row = col + 1; row < m; ++row)
            {
              if (std::abs(B[row][col]) > max_val)
                {
                  max_val = std::abs(B[row][col]);
                  pivot_row = row;
                }
            }
 
          // Swap rows if needed
          if (pivot_row != col)
            {
              std::swap(B[col], B[pivot_row]);
              std::swap(B_inv_explicit[col], B_inv_explicit[pivot_row]);
            }
 
          // Scale pivot row
          T pivot = B[col][col];
          for (size_t j = 0; j < m; ++j)
            {
              B[col][j] /= pivot;
              B_inv_explicit[col][j] /= pivot;
            }
 
          // Eliminate column
          for (size_t row = 0; row < m; ++row)
            {
              if (row != col)
                {
                  T factor = B[row][col];
                  for (size_t j = 0; j < m; ++j)
                    {
                      B[row][j] -= factor * B[col][j];
                      B_inv_explicit[row][j] -= factor * B_inv_explicit[col][j];
                    }
                }
            }
        }
 
      // Clear eta-file after refactorization
      eta_file.clear();
    }
 
  public:
    RevisedSimplex(size_t num_vars, size_t num_constraints)
      : m(num_constraints), n(num_vars),
        c(num_vars, T{0}),
        A(num_constraints, std::vector<T>(num_vars, T{0})),
        b(num_constraints, T{0}),
        basic(num_constraints),
        is_basic(num_vars + num_constraints, false),
        x_B(num_constraints, T{0}),
        solution(num_vars + num_constraints, T{0}),
        pi(num_constraints), d(num_constraints), work(num_constraints),
        eps(std::numeric_limits<T>::epsilon() * 1000),
        refactor_frequency(std::max(size_t{50}, num_constraints)),
        use_steepest_edge(false),
        state(Not_Solved)
    {
      ah_invalid_argument_if(num_vars == 0 || num_constraints == 0)
        << "Number of variables and constraints must be > 0";
    }
 
    void set_objective(size_t j, T coef)
    {
      ah_out_of_range_error_if(j >= n) << "Variable index out of range";
      c[j] = coef;
    }
 
    void set_objective(const T* coefs)
    {
      for (size_t j = 0; j < n; ++j)
        c[j] = coefs[j];
    }
 
    void set_constraint(size_t i, size_t j, T coef)
    {
      ah_out_of_range_error_if(i >= m || j >= n)
        << "Constraint or variable index out of range";
      A[i][j] = coef;
    }
 
    void set_rhs(size_t i, T rhs)
    {
      ah_out_of_range_error_if(i >= m) << "Constraint index out of range";
      b[i] = rhs;
    }
 
    void set_constraint_row(size_t i, const T* coefs, T rhs)
    {
      ah_out_of_range_error_if(i >= m) << "Constraint index out of range";
      for (size_t j = 0; j < n; ++j)
        A[i][j] = coefs[j];
      b[i] = rhs;
    }
 
    void set_refactorization_frequency(size_t freq)
    {
      refactor_frequency = freq;
    }
 
    [[nodiscard]] State solve()
    {
      ah_logic_error_if(state != Not_Solved) << "Already solved";
 
      stats.reset();
      auto start_time = std::chrono::high_resolution_clock::now();
 
      // Initialize basis with slack variables
      for (size_t i = 0; i < m; ++i)
        {
          basic[i] = static_cast<int>(n + i);  // Slack variable index
          is_basic[n + i] = true;
          x_B[i] = b[i];  // Initial basic solution
        }
 
      // Initialize non-basic variable list
      nonbasic.clear();
      for (size_t j = 0; j < n; ++j)
        nonbasic.push_back(static_cast<int>(j));
 
      // Initial B⁻¹ is identity (slack variables)
      B_inv_explicit.assign(m, std::vector<T>(m, T{0}));
      for (size_t i = 0; i < m; ++i)
        B_inv_explicit[i][i] = T{1};
 
      state = Solving;
 
      // Main loop
      while (true)
        {
          ++stats.iterations;
 
          // Pricing: select entering variable
          int entering = select_entering_variable();
 
          if (entering < 0)
            {
              // Optimal!
              state = Solved;
              break;
            }
 
          // Ratio test: select leaving variable
          T theta;
          int leaving = select_leaving_variable(entering, theta);
 
          if (leaving < 0)
            {
              // Unbounded!
              state = Unbounded;
              break;
            }
 
          // Update solution
          for (size_t i = 0; i < m; ++i)
            x_B[i] -= theta * d[i];
          x_B[leaving] = theta;
 
          // Update basis
          int leaving_var = basic[leaving];
          basic[leaving] = entering;
          is_basic[leaving_var] = false;
          is_basic[entering] = true;
 
          // Update non-basic list
          for (auto& nb : nonbasic)
            if (nb == entering)
              {
                nb = leaving_var;
                break;
              }
 
          // Update B⁻¹
          update_basis_inverse(leaving, d);
 
          ++stats.pivots;
        }
 
      // Load final solution
      std::fill(solution.begin(), solution.end(), T{0});
      for (size_t i = 0; i < m; ++i)
        solution[basic[i]] = x_B[i];
 
      auto end_time = std::chrono::high_resolution_clock::now();
      stats.elapsed_ms = std::chrono::duration<double, std::milli>(
        end_time - start_time).count();
 
      return state;
    }
 
    [[nodiscard]] T get_solution(size_t j) const
    {
      ah_out_of_range_error_if(j >= n) << "Variable index out of range";
      return solution[j];
    }
 
    [[nodiscard]] T objective_value() const
    {
      T sum = T{0};
      for (size_t j = 0; j < n; ++j)
        sum += c[j] * solution[j];
      return sum;
    }
 
    [[nodiscard]] State get_state() const noexcept { return state; }
 
    [[nodiscard]] const SimplexStats& get_stats() const noexcept { return stats; }
 
    [[nodiscard]] size_t get_num_vars() const noexcept { return n; }
 
    [[nodiscard]] size_t get_num_constraints() const noexcept { return m; }
 
    [[nodiscard]] bool verify_solution() const
    {
      for (size_t i = 0; i < m; ++i)
        {
          T lhs = T{0};
          for (size_t j = 0; j < n; ++j)
            lhs += A[i][j] * solution[j];
          if (lhs > b[i] + eps)
            return false;
        }
      return true;
    }
 
    void print_stats() const
    {
      std::cout << "=== Revised Simplex Statistics ===\n"
                << "Variables: " << n << ", Constraints: " << m << "\n"
                << "Iterations: " << stats.iterations << "\n"
                << "Pivots: " << stats.pivots << "\n"
                << "Elapsed time: " << stats.elapsed_ms << " ms\n"
                << "Eta-matrices: " << eta_file.size() << "\n"
                << "Refactorizations: "
                << (stats.pivots / refactor_frequency) << "\n";
    }
  };
 
} // end namespace Aleph
 
# endif // SIMPLEX_H