polyhedron.hpp

Go to the documentation of this file.

/*
 * Copyright © 2025-2026 Akreon LLC
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
 
#include "halfspace.hpp"
#include "hyperplane.hpp"
 
namespace convex {
 
template <linalg::RnSpace Space>
using PolyhedronDom = ConvexDomain<Space, VarType::FREE>;
 
template <linalg::RnSpace Space>
using PolyhedronCodom = ConvexDomain<Space, VarType::FREE>;
 
template <linalg::RnSpace Dom, linalg::RnSpace Codom>
using PolyhedronMeta =
    ConvexAtomMeta<PolyhedronDom<Dom>, PolyhedronCodom<Codom>,
                   CurvatureType::CONVEX>;
 
template <typename Meta, typename Child>
struct Polyhedron : ConvexSet<Meta, Polyhedron<Meta, Child>, Child> {
  using ConvexBase = ConvexSet<Meta, Polyhedron<Meta, Child>, Child>;
  using ConvexBase::M;
  using ConvexBase::N;
  using typename ConvexBase::PointT;
  using typename ConvexBase::T;
 
  using Mat = linalg::Matrix<T, M, N>;
  using Bias = linalg::Vector<T, M>;
  using Space = linalg::Matrix<T, 1, N>;
 
  constexpr explicit Polyhedron(const Mat& A_in, const Bias& b_in,
                                const Child& child_in)
      : ConvexBase(child_in), A(A_in), b(b_in) {}
 
  constexpr bool is_well_formed() const {
    for (std::size_t i = 0; i < M; ++i)
      if (!Halfspace<HalfspaceMeta<Space>, Child, HalfspaceType::LH>(
               A.row(i), b(i), std::get<0>(this->children))
               .is_well_formed())
        return false;
    return true;
  }
 
  constexpr bool contains_impl(const PointT& x, T tol) const {
    for (std::size_t i = 0; i < M; ++i) {
      const auto half =
          Halfspace<HalfspaceMeta<Space>, Child, HalfspaceType::LH>(
              A.row(i), b(i), std::get<0>(this->children));
      if (!half.contains(x, tol)) return false;
    }
    return true;
  }
 
  template <linalg::SquareMatrix Mat>
    requires(std::same_as<T, linalg::matrix_scalar_t<Mat>> &&
             linalg::mat_cols_v<Mat> == N)
  constexpr auto transform_impl(const affine::AffineMap<Mat>& map) const {
    const auto& U = map.linear_part();
    const auto& trans = map.translation();
 
    const auto [U_inv, event] = linalg::MatrixInverse<Mat>()(U);
    static_cast<void>(event);
    const auto A_new = A * U_inv;
    const auto b_new = b + A_new * trans;
 
    return Polyhedron<Meta, Child>(A_new, b_new, std::get<0>(this->children));
  }
 
  template <typename Other>
  constexpr bool is_equal_impl(const Other& other) const {
    return other.self().A == A && other.self().b == b;
  }
 
  template <AtomType Type, conic::BuilderLike Builder>
  constexpr auto canonicalize_set(const Builder& builder) const {
    const auto child_builder =
        std::get<0>(this->children)
            .template canonicalize<AtomType::SET>(builder);
    const auto child_coord = child_builder.coord_map;
 
    static constexpr std::size_t ChildCodom =
        linalg::mat_cols_v<decltype(child_coord.linear_part())>;
    static_assert(ChildCodom == N,
                  "Convex set domain does not equal the child coordinate "
                  "linear map's codomain in size");
 
    const auto map = affine::AffineMap<decltype(A)>(T(-1) * A, b);
    const auto rplus = conic::NonnegOrthantCone<linalg::Vector<T, M>>();
    const auto cones = std::make_tuple(rplus);
    const auto emission =
        conic::LocalEmission(map, cones, child_builder.eq_map);
 
    return child_builder.template realize<Type>(emission);
  }
 
  Mat A;  //< set of finite number of halfspaces
  Bias b; //< codomain vector describing halfspace inequalities
};
 
template <typename Mat, typename Bias, typename Child>
Polyhedron(const Mat&, const Bias&, const Child&) -> Polyhedron<
    PolyhedronMeta<dom_space_t<meta_codom_t<atom_meta_t<Child>>>,
                   linalg::Matrix<dom_type_t<meta_codom_t<atom_meta_t<Child>>>,
                                  1, linalg::mat_rows_v<Mat>>>,
    atom_decay_t<Child>>;
 
} // namespace convex