Skip to content

Commit 676e545

Browse files
mtaomtao
committed
delegate simplex kernels to quiver, remove old balsa transforms, fix includes
- volume.hpp, circumcenter.hpp, point_in_circumcircle.hpp: rewritten as thin wrappers delegating to quiver::simplex. Eigen overloads convert via as_zipper() then delegate. - Delete balsa::geometry::quiver::simplex_transforms.hpp (superseded by quiver/simplex/transforms.hpp) - batch ops (volumes.hpp, circumcenters.hpp): replace balsa::zipper:: type aliases with ::zipper:: types directly - tests: update imports/namespaces, add explicit includes for types that were previously provided transitively
1 parent 89fe5f2 commit 676e545

9 files changed

Lines changed: 75 additions & 446 deletions

File tree

‎geometry/include/balsa/geometry/quiver/simplex_transforms.hpp‎

Lines changed: 0 additions & 73 deletions
This file was deleted.

‎geometry/include/balsa/geometry/simplex/circumcenter.hpp‎

Lines changed: 14 additions & 198 deletions
Original file line numberDiff line numberDiff line change
@@ -1,226 +1,42 @@
11
#if !defined(BALSA_GEOMETRY_SIMPLEX_CIRCUMCENTER_HPP)
22
#define BALSA_GEOMETRY_SIMPLEX_CIRCUMCENTER_HPP
3-
#include <zipper/utils/format.hpp>
4-
#include <spdlog/spdlog.h>
5-
#include <stdexcept>
6-
#include <zipper/utils/extents/extent_arithmetic.hpp>
7-
#include <zipper/expression/nullary/Constant.hpp>
8-
#include <zipper/utils/decomposition/qr.hpp>
9-
#include <tuple>
10-
#include "balsa/zipper/types.hpp"
11-
#include "zipper/concepts/Matrix.hpp"
12-
#include "zipper/concepts/Vector.hpp"
13-
#include "zipper/concepts/Zipper.hpp"
14-
15-
namespace balsa::geometry::simplex {
16-
17-
18-
template<::zipper::concepts::Matrix MatType>
19-
requires(MatType::extents_traits::is_dynamic || MatType::extents_type::static_extent(0) + 1 == MatType::extents_type::static_extent(1))
20-
auto circumcenter_spd(const MatType &V);
21-
template<::zipper::concepts::Matrix MatType>
22-
auto circumcenter_spsd(const MatType &V);
23-
24-
25-
template<::zipper::concepts::Matrix MatType>
26-
auto circumcenter(const MatType &V);
27-
28-
template<::zipper::concepts::Matrix SimplexVertices>
29-
auto circumcenter_with_squared_radius(const SimplexVertices &S);
30-
31-
template<::zipper::concepts::Matrix SimplexVertices>
32-
auto circumcenter_with_radius(const SimplexVertices &S);
33-
34-
35-
namespace detail {
36-
37-
// Back-substitution for upper triangular system R*x = b.
38-
// Near-zero pivots are treated as zero (for rank-deficient systems).
39-
template<::zipper::concepts::Matrix RMat, ::zipper::concepts::Zipper BVec>
40-
auto upper_triangular_solve(const RMat &R, const BVec &b) {
41-
using value_type = typename RMat::value_type;
42-
constexpr auto N = RMat::static_extent(1);
43-
::zipper::Vector<value_type, N> x(R.extent(1));
44-
const value_type pivot_tol = std::numeric_limits<value_type>::epsilon() * 100;
45-
for (::zipper::index_type ii = 0; ii < R.extent(1); ++ii) {
46-
::zipper::index_type i = R.extent(1) - 1 - ii;
47-
value_type s = b(i);
48-
for (::zipper::index_type j = i + 1; j < R.extent(1); ++j) {
49-
s -= R(i, j) * x(j);
50-
}
51-
x(i) = (std::abs(R(i, i)) > pivot_tol) ? s / R(i, i) : value_type(0);
52-
}
53-
return x;
54-
}
55-
56-
// Solve A*x = b via QR decomposition with back-substitution.
57-
template<::zipper::concepts::Matrix AMat, ::zipper::concepts::Zipper BVec>
58-
auto qr_solve(const AMat &A, const BVec &b) {
59-
auto [Q, R] = ::zipper::utils::decomposition::qr(A);
60-
auto Qtb = (Q.transpose() * b).eval();
61-
return upper_triangular_solve(R, Qtb);
62-
}
63-
64-
// Check if the QR R factor indicates a rank-deficient matrix
65-
// (any diagonal element near zero).
66-
template<::zipper::concepts::Matrix RMat>
67-
bool qr_is_degenerate(const RMat &R) {
68-
using value_type = typename RMat::value_type;
69-
const value_type tol = std::numeric_limits<value_type>::epsilon() * 100;
70-
::zipper::index_type n = std::min(R.extent(0), R.extent(1));
71-
for (::zipper::index_type i = 0; i < n; ++i) {
72-
if (std::abs(R(i, i)) < tol) {
73-
return true;
74-
}
75-
}
76-
return false;
77-
}
78-
79-
template<::zipper::concepts::Matrix MatType>
80-
auto circumcenter_spsd(const MatType &V) {
81-
using value_type = typename MatType::value_type;
82-
constexpr static ::zipper::index_type compile_cols = MatType::static_extent(1);
83-
constexpr static ::zipper::index_type Msize = ::zipper::utils::extents::plus(compile_cols, 1);
84-
using MType = ::zipper::Matrix<value_type, Msize, Msize>;
853

4+
// Circumcenter computation for simplices.
5+
//
6+
// The zipper-concept implementations live in quiver::simplex and are
7+
// re-exported here for backward compatibility.
868

87-
constexpr static bool static_cols = !MatType::extents_traits::is_dynamic_extent(1);
88-
zipper::index_type simplex_count = V.cols();
89-
auto A = MType(simplex_count + 1, simplex_count + 1);
9+
#include <quiver/simplex/circumcenter.hpp>
9010

11+
namespace balsa::geometry::simplex {
9112

92-
auto ones = ::zipper::expression::nullary::Constant<value_type>(1);
93-
A.row(simplex_count).head(A.extent(0) - 1) = ones;
94-
A.col(simplex_count).head(A.extent(1) - 1) = ones;
95-
A(simplex_count, simplex_count) = 0;
96-
if constexpr (static_cols) {
97-
auto m = A.template slice<
98-
::zipper::static_slice_t<0, compile_cols, 1>,//
99-
::zipper::static_slice_t<0, compile_cols, 1>//
100-
>();
101-
m = 2 * V.transpose() * V;
102-
103-
104-
} else {
105-
auto m = A.slice(
106-
::zipper::slice({}, V.extent(1)),
107-
::zipper::slice({}, V.extent(1)));
108-
m = 2 * V.transpose() * V;
109-
}
110-
auto b = V.colwise().norm_powered().homogeneous().eval();
111-
112-
auto x = qr_solve(A, b);
113-
114-
return (V * x.head(V.cols())).eval();
115-
}
116-
117-
// should be a N,N+1 shape matrix
118-
// if RowCol static then get a shape guarantee
119-
template<::zipper::concepts::Matrix MatType>
120-
requires(MatType::extents_traits::is_dynamic || (MatType::extents_type::static_extent(0) + 1 == MatType::extents_type::static_extent(1)))
121-
auto circumcenter_spd(const MatType &V) {
122-
// 2 V.dot(C) = sum(V.colwise().squaredNorm()).transpose()
123-
124-
assert(V.extent(0) + 1 == V.extent(1));
125-
// probably dont really need this temporary
126-
using ET = MatType::extents_type;
127-
using ETraits = MatType::extents_traits;
128-
constexpr static ::zipper::index_type static_cols = !ETraits::is_dynamic_extent(0) ? ET::static_extent(0) + 1 : ET::static_extent(1);
129-
constexpr static bool has_static_cols = static_cols != std::dynamic_extent;
130-
131-
// todo fix this construction in zipper
132-
auto m = [](const auto &VV) {
133-
if constexpr (has_static_cols) {
134-
using slice_t = ::zipper::static_slice_t<1, static_cols - 1, 1>;
135-
return VV.template slice<::zipper::full_extent_t, slice_t>();
136-
} else {
137-
return VV.template slice<::zipper::full_extent_t>(
138-
::zipper::full_extent_t{},
139-
::zipper::slice(
140-
::zipper::static_index_t<1>{},
141-
VV.extent(1) - 1));
142-
}
143-
}(V)
144-
.eval();
145-
146-
147-
// todo fix this construction in zipper
148-
for (::zipper::index_type j = 0; j < m.extent(1); ++j) {
149-
auto mv = m.col(j);
150-
mv = mv - V.col(0);
151-
}
152-
auto b = m.colwise().norm_powered().eval();
153-
auto A = (2 * m.transpose() * m).eval();
154-
155-
auto [Q, R] = ::zipper::utils::decomposition::qr(A);
156-
if (qr_is_degenerate(R)) {
157-
spdlog::debug("circumcenter_spd: Degenerate simplex detected, using circumcenter_spsd");
158-
return circumcenter_spsd(V);
159-
}
160-
161-
auto Qtb = (Q.transpose() * b).eval();
162-
auto sol = upper_triangular_solve(R, Qtb);
163-
164-
return (V.col(0) + m * sol).eval();
165-
}
166-
167-
}// namespace detail
168-
13+
// ── Zipper overloads (delegate to quiver) ──────────────────────────────
16914

17015
template<::zipper::concepts::Matrix MatType>
17116
requires(MatType::extents_traits::is_dynamic || MatType::extents_type::static_extent(0) + 1 == MatType::extents_type::static_extent(1))
17217
auto circumcenter_spd(const MatType &V) {
173-
174-
return detail::circumcenter_spd(V);
18+
return ::quiver::simplex::circumcenter_spd(V);
17519
}
20+
17621
template<::zipper::concepts::Matrix MatType>
17722
auto circumcenter_spsd(const MatType &V) {
178-
return detail::circumcenter_spsd(V);
23+
return ::quiver::simplex::circumcenter_spsd(V);
17924
}
18025

181-
18226
template<::zipper::concepts::Matrix MatType>
18327
auto circumcenter(const MatType &V) {
184-
// auto v = eigen::as_eigen(V);
185-
// auto c = circumcenter(v);
186-
// return eigen::as_zipper(c).eval();
187-
using extents_type = typename MatType::extents_type;
188-
using extents_traits = typename MatType::extents_traits;
189-
constexpr static ::zipper::index_type static_rows = extents_type::static_extent(0);
190-
constexpr static ::zipper::index_type static_cols = extents_type::static_extent(1);
191-
192-
if constexpr (extents_traits::is_static) {
193-
if constexpr (static_cols == 2) {
194-
return ((V.col(0) + V.col(1)) / typename MatType::value_type(2)).eval();
195-
} else if constexpr (static_cols == static_rows + 1) {
196-
return detail::circumcenter_spd(V);
197-
} else {
198-
return detail::circumcenter_spsd(V);
199-
}
200-
} else {
201-
using RetType = ::zipper::Vector<typename MatType::value_type, extents_type::static_extent(0)>;
202-
if (V.extent(1) == 2) {
203-
return RetType((V.col(0) + V.col(1)) / typename MatType::value_type(2));
204-
} else if (V.extent(0) + 1 == V.extent(1)) {
205-
return RetType(detail::circumcenter_spd(V));
206-
} else {
207-
return RetType(detail::circumcenter_spsd(V));
208-
}
209-
}
28+
return ::quiver::simplex::circumcenter(V);
21029
}
21130

21231
template<::zipper::concepts::Matrix SimplexVertices>
21332
auto circumcenter_with_squared_radius(const SimplexVertices &S) {
214-
auto C = circumcenter(S);
215-
return std::make_tuple(C, (C - S.col(0)).template norm_powered<2>());
33+
return ::quiver::simplex::circumcenter_with_squared_radius(S);
21634
}
21735

21836
template<::zipper::concepts::Matrix SimplexVertices>
21937
auto circumcenter_with_radius(const SimplexVertices &S) {
220-
auto C = circumcenter(S);
221-
return std::make_tuple(C, (C - S.col(0)).norm());
38+
return ::quiver::simplex::circumcenter_with_radius(S);
22239
}
22340

224-
22541
}// namespace balsa::geometry::simplex
226-
#endif// CIRCUMCENTER_H
42+
#endif// BALSA_GEOMETRY_SIMPLEX_CIRCUMCENTER_HPP

‎geometry/include/balsa/geometry/simplex/point_in_circumcircle.hpp‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,17 @@
11
#if !defined(BALSA_GEOMETRY_SIMPLEX_POINT_IN_CIRCUMCIRCLE_HPP)
22
#define BALSA_GEOMETRY_SIMPLEX_POINT_IN_CIRCUMCIRCLE_HPP
33

4-
#include <zipper/concepts/Matrix.hpp>
5-
#include <zipper/concepts/Vector.hpp>
6-
#include "circumcenter.hpp"
7-
namespace balsa::geometry::simplex {
4+
// Point-in-circumcircle test.
5+
//
6+
// Delegates to quiver::simplex::point_in_circumcircle.
7+
8+
#include <quiver/simplex/point_in_circumcircle.hpp>
89

10+
namespace balsa::geometry::simplex {
911

1012
template<::zipper::concepts::Matrix SimplexVertices, ::zipper::concepts::Vector PointType>
1113
bool point_in_circumcircle(const SimplexVertices &S, const PointType &P) {
12-
13-
auto [C, r2] = circumcenter_with_squared_radius(S);
14-
return (C - P).template norm_powered<2>() < r2;
14+
return ::quiver::simplex::point_in_circumcircle(S, P);
1515
}
1616

1717
}// namespace balsa::geometry::simplex

0 commit comments

Comments
 (0)