affine_subspace.hpp

Go to the documentation of this file.

/*
 * Copyright © 2025-2026 Akreon LLC
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
 
#include "convex_set.hpp"
 
namespace convex {
 
template <linalg::RnSpace Space>
using AffineDom = ConvexDomain<Space, VarType::FREE>;
 
template <linalg::RnSpace Space>
using AffineCodom = ConvexDomain<Space, VarType::FREE>;
 
template <linalg::RnSpace Dom, linalg::RnSpace Codom>
using AffineMeta =
    ConvexAtomMeta<AffineDom<Dom>, AffineCodom<Codom>, CurvatureType::AFFINE>;
 
template <typename Meta, typename Child>
struct AffineSubspace : ConvexSet<Meta, AffineSubspace<Meta, Child>, Child> {
  using ConvexBase = ConvexSet<Meta, AffineSubspace<Meta, Child>, Child>;
  using ConvexBase::M;
  using ConvexBase::N;
  using typename ConvexBase::PointT;
  using typename ConvexBase::T;
 
  static_assert(M == N, "Affine subspace's linear map must be square");
 
  using Mat = linalg::Matrix<T, M, N>;
 
  using AmbientVec = linalg::Vector<T, M>;
 
  constexpr explicit AffineSubspace(const Mat& A_in, const AmbientVec& b_in,
                                    const Child& child_in)
      : ConvexBase(child_in), A(A_in), b(b_in) {}
 
  constexpr bool contains_impl(const PointT& x, T tol) const {
    return linalg::L2<AmbientVec>()(A * x - b) < tol;
  }
 
  constexpr bool is_well_formed() const {
    const auto rank = linalg::QR(A).rank();
    return rank == N;
  }
 
  template <typename Other>
  constexpr bool is_equal_impl(const Other& other) const {
    return A == other.self().A && b == other.self().b &&
           other.coord_map == this->coord_map;
  }
 
  template <linalg::SquareMatrix TMat>
    requires(linalg::mat_cols_v<TMat> == N &&
             std::same_as<linalg::matrix_scalar_t<TMat>, T>)
  constexpr auto transform_impl(const affine::AffineMap<TMat>& map) const {
    const auto [U_inv, event] =
        linalg::MatrixInverse<TMat>()(map.linear_part());
    static_cast<void>(event);
    const auto A_tilde = A * U_inv;
    const auto b_tilde = b + A_tilde * map.translation();
    return AffineSubspace<Meta, Child>(A_tilde, b_tilde,
                                       std::get<0>(this->children));
  }
 
  template <AtomType Type, conic::BuilderLike Builder>
  constexpr auto canonicalize_set(const Builder& builder) const {
    static_cast<void>(Type);
    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 eq_map = affine::AffineMap<Mat>(A, b).compose(child_coord);
    const auto emission = conic::LocalEqualityEmission(eq_map);
    return child_builder.realize(emission);
  }
 
  Mat A;        //< basis that spans the space of R^n
  AmbientVec b; //< offset vector which translates the basis to a new origin
};
 
template <typename Mat, typename Vec, typename Child>
AffineSubspace(const Mat&, const Vec&, const Child&)
    -> AffineSubspace<AffineMeta<dom_space_t<meta_codom_t<atom_meta_t<Child>>>,
                                 dom_space_t<meta_codom_t<atom_meta_t<Child>>>>,
                      atom_decay_t<Child>>;
 
} // namespace convex