16 #ifndef OPENKALMAN_EIGEN_LIBRARY_INTERFACE_HPP 17 #define OPENKALMAN_EIGEN_LIBRARY_INTERFACE_HPP 19 #include <type_traits> 27 template<Eigen3::eigen_general<true> T>
28 struct library_interface<T>
34 template<
typename Derived>
36 std::conditional_t<Eigen3::eigen_array_general<T>, Eigen::ArrayBase<Derived>,
37 std::conditional_t<Eigen3::eigen_matrix_general<T>, Eigen::MatrixBase<Derived>, Eigen::EigenBase<Derived>>>>;
41 template<
typename Arg, std::size_t...I,
typename...Ind>
43 check_index_bounds(
const Arg& arg, std::index_sequence<I...>, Ind...ind)
45 ([](
auto max,
auto ind){
if (ind < 0 or ind >= max)
throw std::out_of_range {
46 (
"Index " + std::to_string(I) +
" is out of bounds: it is " + std::to_string(ind) +
47 " but should be in range [0..." + std::to_string(max - 1) +
"].")};
48 }(get_index_dimension_of<I>(arg), ind),...);
52 template<
typename Arg>
53 static constexpr decltype(
auto)
54 get_coeff(Arg&& arg, Eigen::Index i, Eigen::Index j)
60 if constexpr (Eigen3::eigen_DiagonalMatrix<Arg> or Eigen3::eigen_DiagonalWrapper<Arg>)
63 return static_cast<Scalar
>(get_coeff(
nested_object(std::forward<Arg>(arg)), i, 0));
65 return static_cast<Scalar
>(0);
67 else if constexpr (Eigen3::eigen_TriangularView<Arg>)
69 constexpr
int Mode {std::decay_t<Arg>::Mode};
70 if ((i > j and (Mode &
int{Eigen::Upper}) != 0x0) or (i < j and (Mode &
int{Eigen::Lower}) != 0x0))
71 return static_cast<Scalar
>(0);
73 return static_cast<Scalar
>(get_coeff(
nested_object(std::forward<Arg>(arg)), i, j));
75 else if constexpr (Eigen3::eigen_SelfAdjointView<Arg>)
77 return get_coeff(
nested_object(std::forward<Arg>(arg)), i, j);
81 auto evaluator = Eigen::internal::evaluator<std::decay_t<Arg>>(arg);
82 constexpr
bool use_coeffRef = (Eigen::internal::traits<std::decay_t<Arg>>::Flags & Eigen::LvalueBit) != 0x0 and
83 std::is_lvalue_reference_v<Arg> and not std::is_const_v<std::remove_reference_t<Arg>>;
84 if constexpr (use_coeffRef)
return evaluator.coeffRef(i, j);
85 else return evaluator.coeff(i, j);
90 template<
typename Indices>
91 static constexpr std::tuple<Eigen::Index, Eigen::Index>
92 extract_indices(
const Indices& indices)
94 #ifdef __cpp_lib_ranges 95 namespace ranges = std::ranges;
97 auto it = ranges::begin(indices);
98 auto e = ranges::end(indices);
99 if (it == e)
return {0, 0};
100 auto i =
static_cast<std::size_t
>(*it);
101 if (++it == e)
return {i, 0};
102 auto j =
static_cast<std::size_t
>(*it);
103 if (++it == e)
return {i, j};
104 throw std::logic_error(
"Wrong number of indices on component access");
109 #ifdef __cpp_lib_ranges 110 template<
typename Arg, std::ranges::input_range Indices> requires
111 std::convertible_to<std::ranges::range_value_t<Indices>,
const typename std::decay_t<Arg>::Index> and
112 (collections::size_of_v<Indices> ==
dynamic_size or collections::size_of_v<Indices> <= 2)
113 static constexpr values::scalar decltype(
auto)
115 template<
typename Arg,
typename Indices, std::enable_if_t<
116 (collections::size_of_v<Indices> ==
dynamic_size or collections::size_of_v<Indices> <= 2),
int> = 0>
117 static constexpr decltype(
auto)
122 auto [i, j] = extract_indices(indices);
123 if constexpr (Eigen3::eigen_SelfAdjointView<Arg>)
125 constexpr
int Mode {std::decay_t<Arg>::Mode};
126 bool transp = (i > j and (Mode &
int{Eigen::Upper}) != 0) or (i < j and (Mode &
int{Eigen::Lower}) != 0);
127 if constexpr (values::complex<Scalar>)
129 if (transp)
return static_cast<Scalar
>(
values::conj(get_coeff(std::as_const(arg), j, i)));
130 else return static_cast<Scalar
>(get_coeff(std::as_const(arg), i, j));
134 if (transp)
return get_coeff(std::forward<Arg>(arg), j, i);
135 else return get_coeff(std::forward<Arg>(arg), i, j);
140 return get_coeff(std::forward<Arg>(arg), i, j);
145 #ifdef __cpp_lib_ranges 146 template<
typename Arg, std::ranges::input_range Indices> requires (not std::is_const_v<Arg>) and
147 std::convertible_to<std::ranges::range_value_t<Indices>,
const typename Arg::Index> and
148 (collections::size_of_v<Indices> ==
dynamic_size or collections::size_of_v<Indices> <= 2) and
149 std::assignable_from<decltype(get_coeff(std::declval<Arg&>(), 0, 0)),
const scalar_type_of_t<Arg>&> and
150 ((Eigen::internal::traits<std::decay_t<Arg>>::Flags & Eigen::LvalueBit) != 0) and
151 (not Eigen3::eigen_DiagonalWrapper<Arg>) and (not Eigen3::eigen_TriangularView<Arg>)
153 template<
typename Arg,
typename Indices, std::enable_if_t<(not std::is_const_v<Arg>) and
154 (collections::size_of_v<Indices> == dynamic_size or collections::size_of_v<Indices> <= 2) and
155 std::is_assignable<decltype(get_coeff(std::declval<Arg&>(), 0, 0)), const scalar_type_of_t<Arg>&>::value and
156 (Eigen::
internal::traits<std::decay_t<Arg>>::Flags & Eigen::LvalueBit) != 0 and
157 (not Eigen3::eigen_DiagonalWrapper<Arg>) and (not Eigen3::eigen_TriangularView<Arg>),
int> = 0>
160 set_component(Arg& arg, const scalar_type_of_t<Arg>& s, const Indices& indices)
162 using Scalar = scalar_type_of_t<Arg>;
163 auto [i, j] = extract_indices(indices);
164 if constexpr (Eigen3::eigen_SelfAdjo
intView<Arg>)
166 constexpr
int Mode {std::decay_t<Arg>::Mode};
167 bool transp = (i > j and (Mode &
int{Eigen::Upper}) != 0) or (i < j and (Mode &
int{Eigen::Lower}) != 0);
168 if constexpr (values::complex<Scalar>)
171 else get_coeff(arg, i, j) = s;
175 if (transp) get_coeff(arg, j, i) = s;
176 else get_coeff(arg, i, j) = s;
181 get_coeff(arg, i, j) = s;
187 template<
typename Arg>
188 static decltype(
auto)
189 wrap_if_nests_by_reference(Arg&& arg)
191 if constexpr (Eigen3::eigen_general<Arg>)
193 constexpr
auto Flags = Eigen::internal::traits<std::remove_reference_t<Arg>>::Flags;
194 if constexpr (std::is_lvalue_reference_v<Arg> and static_cast<bool>(Flags & Eigen::NestByRefBit))
195 return std::forward<Arg>(arg);
197 return Eigen3::make_eigen_wrapper(std::forward<Arg>(arg));
201 return Eigen3::make_eigen_wrapper(std::forward<Arg>(arg));
207 template<
typename Arg>
208 static decltype(
auto)
211 if constexpr (Eigen3::eigen_wrapper<Arg>)
213 return std::forward<Arg>(arg);
215 else if constexpr (internal::library_wrapper<Arg>)
219 else if constexpr (Eigen3::eigen_ArrayWrapper<Arg>)
223 else if constexpr (Eigen3::eigen_array_general<Arg>)
225 return wrap_if_nests_by_reference(std::forward<Arg>(arg)).matrix();
227 else if constexpr (Eigen3::eigen_matrix_general<Arg>)
229 return wrap_if_nests_by_reference(std::forward<Arg>(arg));
231 else if constexpr (not Eigen3::eigen_general<Arg> and directly_accessible<Arg> and std::is_lvalue_reference_v<Arg>)
234 constexpr
int rows = dynamic_dimension<Arg, 0> ? Eigen::Dynamic : index_dimension_of_v<Arg, 0>;
235 constexpr
int cols = dynamic_dimension<Arg, 1> ? Eigen::Dynamic : index_dimension_of_v<Arg, 1>;
239 auto [s0, s1] = internal::strides(arg);
240 using S0 = decltype(s0);
241 using S1 = decltype(s1);
242 constexpr
int es0 = []() ->
int {
243 if constexpr (values::fixed<S0>)
return static_cast<std::ptrdiff_t
>(S0{});
244 else return Eigen::Dynamic;
246 constexpr
int es1 = []() ->
int {
247 if constexpr (values::fixed<S1>)
return static_cast<std::ptrdiff_t
>(S1{});
248 else return Eigen::Dynamic;
250 using IndexType =
typename std::decay_t<Arg>::Index;
251 auto is0 =
static_cast<IndexType
>(
static_cast<std::ptrdiff_t
>(s0));
252 auto is1 =
static_cast<IndexType
>(
static_cast<std::ptrdiff_t
>(s1));
254 if constexpr (values::fixed<S0> and
255 (es0 == 1 or (values::fixed<S1> and es0 < es1)))
257 using M = Eigen::Matrix<Scalar, rows, cols, Eigen::ColMajor>;
258 Eigen::Stride<es1, es0> strides {is1, is0};
259 return Eigen::Map<const M, Eigen::Unaligned, decltype(strides)> {internal::raw_data(arg), rows, cols, strides};
261 else if constexpr (values::fixed<S1> and
262 (es1 == 1 or (values::fixed<S0> and es1 < es0)))
264 using M = Eigen::Matrix<Scalar, rows, cols, Eigen::RowMajor>;
265 Eigen::Stride<es0, es1> strides {is0, is1};
266 return Eigen::Map<const M, Eigen::Unaligned, decltype(strides)> {internal::raw_data(arg), rows, cols, strides};
272 using M = Eigen::Matrix<Scalar, rows, cols, Eigen::RowMajor>;
273 Eigen::Stride<es0, es1> strides {is0, is1};
274 return Eigen::Map<const M, Eigen::Unaligned, decltype(strides)> {internal::raw_data(arg), rows, cols, strides};
278 using M = Eigen::Matrix<Scalar, rows, cols, Eigen::ColMajor>;
279 Eigen::Stride<es1, es0> strides {is1, is0};
280 return Eigen::Map<const M, Eigen::Unaligned, decltype(strides)> {internal::raw_data(arg), rows, cols, strides};
286 constexpr
auto l = layout_of_v<Arg> ==
Layout::right ? Eigen::RowMajor : Eigen::ColMajor;
287 using M = Eigen::Matrix<Scalar, rows, cols, l>;
288 return Eigen::Map<const M> {internal::raw_data(arg), rows, cols};
293 return Eigen3::make_eigen_wrapper(std::forward<Arg>(arg));
299 template<
typename Scalar,
int rows,
int cols, auto options>
300 using dense_type = std::conditional_t<Eigen3::eigen_array_general<T, true>,
301 Eigen::Array<Scalar, rows, cols, options>, Eigen::Matrix<Scalar, rows, cols, options>>;
303 template<
typename Scalar, std::
size_t rows, std::
size_t cols, auto options>
304 using writable_type = dense_type<Scalar,
305 (rows ==
dynamic_size ? Eigen::Dynamic :
static_cast<int>(rows)),
306 (cols ==
dynamic_size ? Eigen::Dynamic :
static_cast<int>(cols)), options>;
310 #ifdef __cpp_concepts 311 template<
typename To, Eigen3::eigen_general From> requires (std::assignable_from<To&, From&&>)
313 template<
typename To,
typename From, std::enable_if_t<Eigen3::eigen_general<From> and std::is_assignable_v<To&, From&&>,
int> = 0>
315 static void assign(To& a, From&& b)
317 if constexpr (Eigen3::eigen_DiagonalWrapper<From>)
320 a = std::forward<From>(b);
322 a =
diagonal_of(std::forward<From>(b)).asDiagonal();
326 a = std::forward<From>(b);
332 template<
typename Descriptors>
333 static constexpr decltype(
auto)
334 extract_descriptors(Descriptors&& descriptors)
336 if constexpr (pattern_tuple<Descriptors>)
338 constexpr
auto dim = std::tuple_size_v<std::decay_t<Descriptors>>;
339 static_assert(dim <= 2);
341 else if constexpr (dim == 1)
return std::tuple {std::get<0>(std::forward<Descriptors>(descriptors)),
coordinates::Axis{}};
342 else return std::forward<Descriptors>(descriptors);
346 #ifdef __cpp_lib_ranges 347 namespace ranges = std::ranges;
349 auto it = ranges::begin(descriptors);
350 auto e = ranges::end(descriptors);
355 if (++it == e)
return std::tuple {i, j};
356 throw std::logic_error(
"Wrong number of vector space descriptors");
362 #ifdef __cpp_concepts 363 template<Layout layout, values::number Scalar, coordinates::eucl
idean_pattern_collection Ds>
365 template<Layout layout,
typename Scalar,
typename Ds, std::enable_if_t<values::number<Scalar> and
366 coordinates::eucl
idean_pattern_collection<Ds>,
int> = 0>
368 static auto make_default(Ds&& ds)
370 using IndexType =
typename Eigen::Index;
372 if constexpr (pattern_tuple<Ds> and collections::size_of_v<Ds> > 2)
374 constexpr
auto options = layout ==
Layout::right ? Eigen::RowMajor : Eigen::ColMajor;
376 if constexpr (fixed_pattern_tuple<Ds>)
378 auto sizes = std::apply([](
auto&&...d){
379 return Eigen::Sizes<static_cast<std::ptrdiff_t>(coordinates::dimension_of_v<decltype(d)>)...> {};
380 }, std::forward<Ds>(ds));
382 return Eigen::TensorFixedSize<Scalar, decltype(sizes), options, IndexType> {};
386 return std::apply([](
auto&&...d){
387 using Ten = Eigen::Tensor<Scalar, collections::size_of_v<Ds>, options, IndexType>;
388 return Ten {
static_cast<IndexType
>(get_dimension(d))...};
389 }, std::forward<Ds>(ds));
394 auto [d0, d1] = extract_descriptors(std::forward<Ds>(ds));
395 constexpr
auto dim0 = coordinates::dimension_of_v<decltype(d0)>;
396 constexpr
auto dim1 = coordinates::dimension_of_v<decltype(d1)>;
398 static_assert(layout !=
Layout::right or dim0 == 1 or dim1 != 1,
399 "Eigen does not allow creation of a row-major column vector.");
400 static_assert(layout !=
Layout::left or dim0 != 1 or dim1 == 1,
401 "Eigen does not allow creation of a column-major row vector.");
403 constexpr
auto options =
405 Eigen::RowMajor : Eigen::ColMajor;
407 using M = writable_type<Scalar, dim0, dim1, options>;
410 return M {
static_cast<IndexType
>(get_dimension(d0)), static_cast<IndexType>(get_dimension(d1))};
419 #ifdef __cpp_concepts 420 template<Layout layout, writable Arg, std::convertible_to<scalar_type_of_t<Arg>> S, std::convertible_to<scalar_type_of_t<Arg>>...Ss>
423 template<
Layout layout,
typename Arg,
typename S,
typename...Ss, std::enable_if_t<writable<Arg> and
425 std::conjunction<std::is_convertible<S, typename scalar_type_of<Arg>::type>,
426 std::is_convertible<Ss, typename scalar_type_of<Arg>::type>...>
::value,
int> = 0>
430 if constexpr (layout ==
Layout::left) ((arg.transpose() << s), ... , ss);
431 else ((arg << s), ... , ss);
435 #ifdef __cpp_lib_ranges 436 template<values::dynamic C, coordinates::eucl
idean_pattern_collection Ds> requires
437 (collections::size_of_v<Ds> !=
dynamic_size) and (collections::size_of_v<Ds> <= 2)
440 template<
typename C,
typename Ds, std::enable_if_t<values::dynamic<C> and
441 coordinates::eucl
idean_pattern_collection<Ds> and
442 (collections::size_of_v<Ds> != dynamic_size) and (collections::size_of_v<Ds> <= 2)),
int> = 0>
443 static constexpr auto
445 make_constant(C&& c, Ds&& ds)
447 auto [d0, d1] = extract_descriptors(std::forward<Ds>(ds));
448 constexpr auto dim0 = coordinates::dimension_of_v<decltype(d0)>;
449 constexpr auto dim1 = coordinates::dimension_of_v<decltype(d1)>;
451 auto value = values::to_number(std::forward<C>(c));
452 using Scalar = std::decay_t<decltype(value)>;
453 constexpr auto options = (dim0 == 1 and dim1 != 1) ? Eigen::RowMajor : Eigen::ColMajor;
454 using M = writable_type<Scalar, dim0, dim1, options>;
456 using IndexType =
typename Eigen::Index;
458 if constexpr (fixed_pattern_collection<Ds>)
459 return M::Constant(value);
461 return M::Constant(static_cast<IndexType>(dim0), static_cast<IndexType>(get_dimension(d1)), value);
466 #ifdef __cpp_concepts
467 template<
typename Scalar, coordinates::eucl
idean_pattern_collection Ds> requires
468 (collections::size_of_v<Ds> !=
dynamic_size) and (collections::size_of_v<Ds> <= 2)
471 template<
typename Scalar,
typename...Ds, std::enable_if_t<coordinates::euclidean_pattern_collection<Ds> and
472 (collections::size_of_v<Ds> !=
dynamic_size) and (collections::size_of_v<Ds> <= 2),
int> = 0>
473 static constexpr
auto 475 make_identity_matrix(Ds&& ds)
477 auto [d0, d1] = extract_descriptors(std::forward<Ds>(ds));
478 constexpr
auto dim0 = coordinates::dimension_of_v<decltype(d0)>;
479 constexpr
auto dim1 = coordinates::dimension_of_v<decltype(d1)>;
481 constexpr
auto options = (dim0 == 1 and dim1 != 1) ? Eigen::RowMajor : Eigen::ColMajor;
482 using M = writable_type<Scalar, dim0, dim1, options>;
484 using IndexType =
typename Eigen::Index;
486 if constexpr (fixed_pattern_collection<Ds>)
487 return M::Identity();
489 return M::Identity(static_cast<IndexType>(get_dimension(d0)), static_cast<IndexType>(get_dimension(d1)));
493 #ifdef __cpp_concepts 494 template<TriangleType t, Eigen3::eigen_dense_general Arg> requires std::is_lvalue_reference_v<Arg> or
495 (not std::is_lvalue_reference_v<
typename Eigen::internal::ref_selector<std::decay_t<Arg>>::non_const_type>)
497 template<TriangleType t,
typename Arg, std::enable_if_t<Eigen3::eigen_dense_general<Arg> and (std::is_lvalue_reference_v<Arg> or
498 not std::is_lvalue_reference_v<
typename Eigen::
internal::ref_selector<std::remove_reference_t<Arg>>::non_const_type>),
int> = 0>
500 static constexpr
auto 504 return arg.template triangularView<Mode>();
508 #ifdef __cpp_concepts 509 template<HermitianAdapterType t, Eigen3::eigen_dense_general Arg> requires std::is_lvalue_reference_v<Arg>
511 template<HermitianAdapterType t,
typename Arg, std::enable_if_t<Eigen3::eigen_dense_general<Arg> and std::is_lvalue_reference_v<Arg>,
int> = 0>
513 static constexpr
auto 514 make_hermitian_adapter(Arg&& arg)
517 return arg.template selfadjointView<Mode>();
528 #ifdef __cpp_concepts 529 template<
typename Arg,
typename...Begin,
typename...Size> requires (
sizeof...(Begin) <= 2)
531 template<
typename Arg,
typename...Begin,
typename...Size, std::enable_if_t<(
sizeof...(Begin) <= 2),
int> = 0>
534 get_slice(Arg&& arg,
const std::tuple<Begin...>& begin,
const std::tuple<Size...>&
size)
536 auto b0 = [](
const auto& begin){
537 using Begin0 = std::decay_t<decltype(std::get<0>(begin))>;
538 if constexpr (values::fixed<Begin0>)
return std::integral_constant<Eigen::Index, Begin0::value>{};
539 else return static_cast<Eigen::Index
>(std::get<0>(begin));
542 auto b1 = [](
const auto& begin){
543 if constexpr (
sizeof...(Begin) < 2)
return std::integral_constant<Eigen::Index, 0>{};
546 using Begin1 = std::decay_t<decltype(std::get<1>(begin))>;
547 if constexpr (values::fixed<Begin1>)
return std::integral_constant<Eigen::Index, Begin1::value>{};
548 else return static_cast<Eigen::Index
>(std::get<1>(begin));
552 auto s0 = [](
const auto&
size){
553 using Size0 = std::decay_t<decltype(std::get<0>(
size))>;
554 if constexpr (values::fixed<Size0>)
return std::integral_constant<Eigen::Index, Size0::value>{};
555 else return static_cast<Eigen::Index
>(std::get<0>(
size));
558 auto s1 = [](
const auto&
size){
559 if constexpr (
sizeof...(Size) < 2)
return std::integral_constant<Eigen::Index, 1>{};
562 using Size1 = std::decay_t<decltype(std::get<1>(
size))>;
563 if constexpr (values::fixed<Size1>)
return std::integral_constant<Eigen::Index, Size1::value>{};
564 else return static_cast<Eigen::Index
>(std::get<1>(
size));
568 constexpr
int S0 =
static_cast<int>(values::fixed<decltype(s0)> ?
static_cast<Eigen::Index
>(s0) : Eigen::Dynamic);
569 constexpr
int S1 =
static_cast<int>(values::fixed<decltype(s1)> ?
static_cast<Eigen::Index
>(s1) : Eigen::Dynamic);
572 using M = decltype(m);
574 constexpr
auto Flags = Eigen::internal::traits<std::remove_reference_t<M>>::Flags;
576 if constexpr (directly_accessible<M> and not (std::is_lvalue_reference_v<M> and static_cast<bool>(Flags & Eigen::NestByRefBit)))
579 auto rep {std::forward<M>(m).
template replicate<1,1>()};
580 using B = Eigen::Block<const decltype(rep), S0, S1>;
581 if constexpr ((values::fixed<Size> and ...))
582 return B {std::move(rep), std::move(b0), std::move(b1)};
584 return B {std::move(rep), std::move(b0), std::move(b1), std::move(s0), std::move(s1)};
588 using M_noref = std::remove_reference_t<M>;
589 using XprType = std::conditional_t<std::is_lvalue_reference_v<M>, M_noref,
const M_noref>;
590 using B = Eigen::Block<XprType, S0, S1>;
591 if constexpr ((values::fixed<Size> and ...))
592 return B {std::forward<M>(m), std::move(b0), std::move(b1)};
594 return B {std::forward<M>(m), std::move(b0), std::move(b1), std::move(s0), std::move(s1)};
599 #ifdef __cpp_concepts 600 template<Eigen3::eigen_dense_general Arg, Eigen3::eigen_general Block,
typename...Begin> requires (
sizeof...(Begin) <= 2)
602 template<
typename Arg,
typename Block,
typename...Begin, std::enable_if_t<
603 Eigen3::eigen_dense_general<Arg> and Eigen3::eigen_general<Block> and (
sizeof...(Begin) <= 2),
int> = 0>
606 set_slice(Arg& arg, Block&& block,
const Begin&...begin)
608 auto [b0, b1] = [](Eigen::Index bs0, Eigen::Index bs1,
const auto&...){
return std::tuple {bs0, bs1}; }
609 (
static_cast<std::size_t
>(begin)..., 0_uz, 0_uz);
611 if constexpr (Eigen3::eigen_Block<Block>)
613 if (std::addressof(arg) == std::addressof(block.nestedExpression()) and b0 == block.startRow() and b1 == block.startCol())
617 constexpr
auto Bx0 =
static_cast<int>(index_dimension_of_v<Block, 0>);
618 constexpr
auto Bx1 =
static_cast<int>(index_dimension_of_v<Block, 1>);
619 using Bk = Eigen::Block<std::remove_reference_t<Arg>, Bx0, Bx1>;
621 if constexpr (not has_dynamic_dimensions<Block>)
623 Bk {arg, b0, b1} = std::forward<Block>(block);
627 auto s0 =
static_cast<Eigen::Index
>(get_index_dimension_of<0>(block));
628 auto s1 =
static_cast<Eigen::Index
>(get_index_dimension_of<1>(block));
629 Bk {arg, b0, b1, s0, s1} = std::forward<Block>(block);
634 #ifdef __cpp_concepts 635 template<TriangleType t, Eigen3::eigen_SelfAdjo
intView A, Eigen3::eigen_general B> requires (not hermitian_matrix<A>) and
636 set_triangle_defined_for<T, t, decltype(OpenKalman::nested_object(std::declval<A>())), B&&>
638 template<
TriangleType t,
typename A,
typename B, std::enable_if_t<
639 Eigen3::eigen_general<B> and Eigen3::eigen_SelfAdjointView<A> and (not hermitian_matrix<A>) and
640 set_triangle_defined_for<T, t, decltype(OpenKalman::nested_object(std::declval<A>())), B&&>,
int> = 0>
643 set_triangle(
A&& a, B&& b)
654 #ifdef __cpp_concepts 655 template<TriangleType t,
typename A, Eigen3::eigen_general B> requires
657 std::assignable_from<decltype(OpenKalman::diagonal_of(std::declval<A&&>())), decltype(
OpenKalman::diagonal_of(std::declval<B&&>()))>
659 template<
TriangleType t,
typename A,
typename B, std::enable_if_t<
661 std::is_assignable_v<decltype(OpenKalman::diagonal_of(std::declval<A&&>())), decltype(
OpenKalman::diagonal_of(std::declval<B&&>()))>,
int> = 0>
664 set_triangle(
A&& a, B&& b)
670 #ifdef __cpp_concepts 671 template<TriangleType t,
typename A, Eigen3::eigen_general B> requires
672 (Eigen3::eigen_MatrixWrapper<A> or Eigen3::eigen_ArrayWrapper<A>) and
673 set_triangle_defined_for<T, t, decltype(OpenKalman::nested_object(std::declval<A>())), B&&>
675 template<
TriangleType t,
typename A,
typename B, std::enable_if_t<Eigen3::eigen_general<B> and
676 (Eigen3::eigen_MatrixWrapper<A> or Eigen3::eigen_ArrayWrapper<A>) and
677 set_triangle_defined_for<T, t, decltype(OpenKalman::nested_object(std::declval<A>())), B&&>,
int> = 0>
680 set_triangle(
A&& a, B&& b)
686 #ifdef __cpp_concepts 687 template<TriangleType t, Eigen3::eigen_dense_general A, Eigen3::eigen_general B> requires
688 (not Eigen3::eigen_MatrixWrapper<A>) and (not Eigen3::eigen_ArrayWrapper<A>) and
691 template<
TriangleType t,
typename A,
typename B, std::enable_if_t<
692 Eigen3::eigen_dense_general<A> and Eigen3::eigen_general<B> and
693 (not Eigen3::eigen_MatrixWrapper<A>) and (not Eigen3::eigen_ArrayWrapper<A>) and
697 set_triangle(
A&& a, B&& b)
699 a.template triangularView<t == TriangleType::upper ? Eigen::Upper : Eigen::Lower>() = std::forward<B>(b);
703 #ifdef __cpp_concepts 704 template<Eigen3::eigen_dense_general Arg> requires square_shaped<Arg> or dimension_size_of_index_is<Arg, 0, 1> or
705 std::is_lvalue_reference_v<Arg> or (not std::is_lvalue_reference_v<
typename std::decay_t<Arg>::Nested>)
707 template<
typename Arg, std::enable_if_t<Eigen3::eigen_dense_general<Arg> and
708 (square_shaped<Arg> or dimension_size_of_index_is<Arg, 0, 1> or std::is_lvalue_reference_v<Arg> or
709 not std::is_lvalue_reference_v<
typename std::decay_t<Arg>::Nested>),
int> = 0>
711 static constexpr
auto 714 if constexpr (not vector<Arg>)
if (not
is_vector(arg))
throw std::invalid_argument {
715 "Argument of to_diagonal must have 1 column; instead it has " + std::to_string(get_index_dimension_of<1>(arg))};
717 if constexpr (square_shaped<Arg> or dimension_size_of_index_is<Arg, 0, 1>)
719 return internal::make_fixed_size_adapter(std::forward<Arg>(arg));
721 else if constexpr (Eigen3::eigen_array_general<Arg>)
723 return arg.matrix().asDiagonal();
727 return arg.asDiagonal();
732 #ifdef __cpp_concepts 735 template<
typename Arg, std::enable_if_t<Eigen3::eigen_SelfAdjo
intView<Arg>,
int> = 0>
737 static constexpr
auto 745 #ifdef __cpp_concepts 746 template<Eigen3::eigen_TriangularView Arg> requires (triangle_type_of_v<Arg> ==
TriangleType::lower)
748 template<
typename Arg, std::enable_if_t<Eigen3::eigen_TriangularView<Arg> and (triangle_type_of_v<Arg> == TriangleType::lower),
int> = 0>
750 static constexpr
auto 758 template<
typename Arg>
759 #ifdef __cpp_concepts 760 static constexpr
vector decltype(
auto)
762 static constexpr decltype(
auto)
768 if constexpr (Eigen3::eigen_DiagonalWrapper<Arg>)
770 using Diag = decltype(
nested_object(std::forward<Arg>(arg)));
771 using EigenTraits = Eigen::internal::traits<std::decay_t<Diag>>;
772 constexpr
auto rows = EigenTraits::RowsAtCompileTime;
773 constexpr
auto cols = EigenTraits::ColsAtCompileTime;
775 static_assert(cols != 1,
"For Eigen::DiagonalWrapper<T> interface, T should never be a column vector " 776 "because diagonal_of function handles this case.");
777 if constexpr (cols == 0)
782 else if constexpr (rows == 1 or rows == 0)
787 else if constexpr (rows == Eigen::Dynamic or cols == Eigen::Dynamic)
790 using M = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
792 get_index_dimension_of<0>(diag) * get_index_dimension_of<1>(diag))};
796 using M = Eigen::Matrix<Scalar, rows * cols, 1>;
800 else if constexpr (Eigen3::eigen_SelfAdjointView<Arg> or Eigen3::eigen_TriangularView<Arg>)
804 else if constexpr (Eigen3::eigen_Identity<Arg>)
806 auto f = [](
const auto& a,
const auto& b) {
return std::min(a, b); };
807 auto dim =
values::operation{f, get_index_dimension_of<0>(arg), get_index_dimension_of<1>(arg)};
808 return make_constant<Arg, Scalar, 1>(dim);
810 else if constexpr (Eigen3::eigen_dense_general<Arg>)
813 using M = decltype(m);
814 using M_noref = std::remove_reference_t<M>;
815 using MatrixType = std::conditional_t<std::is_lvalue_reference_v<M>, M_noref,
const M_noref>;
816 return Eigen::Diagonal<MatrixType, 0> {std::forward<M>(m)};
825 template<
typename Arg,
typename Factor0 = std::
integral_constant<std::
size_t, 1>,
typename Factor1 = std::
integral_constant<std::
size_t, 1>>
827 broadcast(Arg&& arg,
const Factor0& factor0 = Factor0{},
const Factor1& factor1 = Factor1{})
829 constexpr
int F0 = []{
830 if constexpr (values::fixed<Factor0>)
return static_cast<std::size_t
>(Factor0{});
831 else return Eigen::Dynamic;
833 constexpr
int F1 = []{
834 if constexpr (values::fixed<Factor1>)
return static_cast<std::size_t
>(Factor1{});
835 else return Eigen::Dynamic;
838 using IndexType =
typename std::decay_t<Arg>::Index;
841 using M = decltype(m);
842 using R = Eigen::Replicate<std::remove_reference_t<M>, F0, F1>;
843 if constexpr (values::fixed<Factor0> and values::fixed<Factor1>)
844 return R {std::forward<M>(m)};
846 return R {std::forward<M>(m), static_cast<IndexType>(factor0),
static_cast<IndexType
>(factor1)};
852 template<
typename Op,
typename...S>
853 static constexpr
auto dummy_op(Op op, S...s)
855 if constexpr (std::is_invocable_v<Op, S...>)
return op(s...);
856 else if constexpr (std::is_invocable_v<Op, std::size_t, std::size_t>)
return op(std::size_t{0}, std::size_t{0});
857 else return op(std::size_t{0});
861 template<
typename Op,
typename...Args>
864 template<
typename Op,
typename Arg>
865 struct EigenNaryOp<Op, Arg> {
using type = Eigen::CwiseUnaryOp<Op, Arg>; };
867 template<
typename Op,
typename Arg1,
typename Arg2>
868 struct EigenNaryOp<Op, Arg1, Arg2> {
using type = Eigen::CwiseBinaryOp<Op, Arg1, Arg2>; };
870 template<
typename Op,
typename Arg1,
typename Arg2,
typename Arg3>
871 struct EigenNaryOp<Op, Arg1, Arg2, Arg3> {
using type = Eigen::CwiseTernaryOp<Op, Arg1, Arg2, Arg3>; };
875 #ifdef __cpp_concepts 876 template<coordinates::pattern...Ds,
typename Operation,
indexible...Args> requires
877 (
sizeof...(Ds) <= 2) and (
sizeof...(Args) <= 3) and
879 (
sizeof...(Args) == 0 and
880 (
values::number<std::invoke_result_t<Operation, std::conditional_t<true, std::size_t, Ds>...>> or
883 template<
typename...Ds,
typename Operation,
typename...Args, std::enable_if_t<
sizeof...(Ds) <= 2 and
sizeof...(Args) <= 3 and
884 (coordinates::pattern<Ds> and ...) and (indexible<Args> and ...) and
886 (
sizeof...(Args) == 0 and
887 (
values::number<
typename std::invoke_result<Operation, std::conditional_t<true, std::size_t, Ds>...>::type> or
888 values::number<
typename std::invoke_result<Operation, std::size_t>::type>))),
int> = 0>
893 decltype(
auto) op = Eigen3::native_operation(std::forward<Operation>(
operation));
894 using Op = decltype(op);
897 if constexpr (
sizeof...(Args) == 0)
900 return Eigen::CwiseNullaryOp<std::remove_reference_t<Op>, P> {
901 static_cast<typename P::Index
>(get_dimension(std::get<0>(tup))),
902 static_cast<typename P::Index>(get_dimension(std::get<1>(tup))),
903 std::forward<Op>(op)};
907 auto seq = std::index_sequence_for<Ds...>{};
908 using CW =
typename EigenNaryOp<std::decay_t<Op>, std::remove_reference_t<Args>...>::type;
909 return CW {std::forward<Args>(args)..., std::forward<Op>(op)};
914 #ifdef __cpp_concepts 915 template<std::size_t...indices,
typename BinaryFunction,
typename Arg> requires ((indices <= 1) and ...)
917 template<std::size_t...indices,
typename BinaryFunction,
typename Arg, std::enable_if_t<((indices <= 1) and ...),
int> = 0>
920 reduce(BinaryFunction&& b, Arg&& arg)
922 if constexpr (Eigen3::eigen_dense_general<Arg>)
924 auto&& op = Eigen3::native_operation(std::forward<BinaryFunction>(b));
925 using Op = decltype(op);
927 if constexpr (((indices == 0) or ...) and ((indices == 1) or ...))
929 return std::forward<Arg>(arg).redux(std::forward<Op>(op));
933 constexpr
auto dir = ((indices == 0) and ...) ? Eigen::Vertical : Eigen::Horizontal;
936 using M = decltype(m);
937 using M_noref = std::remove_reference_t<M>;
938 using MatrixType = std::conditional_t<std::is_lvalue_reference_v<M>, M_noref,
const M_noref>;
939 return Eigen::PartialReduxExpr<MatrixType, ROp, dir> {std::forward<M>(m), ROp{std::forward<Op>(op)}};
949 template<
typename Arg>
950 static constexpr decltype(
auto)
953 if constexpr (Eigen3::eigen_dense_general<Arg>)
955 using UnaryOp = Eigen::internal::scalar_conjugate_op<scalar_type_of_t<Arg>>;
957 using M = decltype(m);
958 return Eigen::CwiseUnaryOp<UnaryOp, std::remove_reference_t<M>> {std::forward<M>(m)};
960 else if constexpr (Eigen3::eigen_TriangularView<Arg> or Eigen3::eigen_SelfAdjointView<Arg>)
962 return std::forward<Arg>(arg).
conjugate();
964 else if constexpr (triangular_matrix<Arg>)
976 template<
typename Arg>
977 static constexpr decltype(
auto)
980 if constexpr (Eigen3::eigen_matrix_general<Arg>)
983 using M = decltype(m);
984 using M_noref = std::remove_reference_t<M>;
985 using MatrixType = std::conditional_t<std::is_lvalue_reference_v<M>, M_noref,
const M_noref>;
986 return Eigen::Transpose<MatrixType> {std::forward<M>(m)};
988 else if constexpr (Eigen3::eigen_TriangularView<Arg> or Eigen3::eigen_SelfAdjointView<Arg>)
990 return std::forward<Arg>(arg).
transpose();
992 else if constexpr (triangular_matrix<Arg>)
1006 #ifdef __cpp_concepts 1007 template<
typename Arg> requires (not Eigen3::eigen_dense_general<Arg>)
1009 template<
typename Arg, std::enable_if_t<not Eigen3::eigen_dense_general<Arg>,
int> = 0>
1011 static constexpr decltype(
auto)
1014 if constexpr (Eigen3::eigen_TriangularView<Arg> or Eigen3::eigen_SelfAdjointView<Arg>)
1016 return std::forward<Arg>(arg).
adjoint();
1018 else if constexpr (triangular_matrix<Arg>)
1032 template<
typename Arg>
1033 static constexpr
auto 1036 if constexpr (Eigen3::eigen_matrix_general<Arg, true>)
1038 else if constexpr (Eigen3::eigen_array_general<Arg, true>)
1039 return std::forward<Arg>(arg).matrix().determinant();
1046 #ifdef __cpp_concepts 1047 template<
typename A, Eigen3::eigen_general B>
1049 template<
typename A,
typename B, std::enable_if_t<Eigen3::eigen_general<B>,
int> = 0>
1051 static constexpr
auto 1057 auto f = [](
auto&& x) -> decltype(
auto) {
1058 constexpr
bool herm = hermitian_matrix<A> and hermitian_matrix<B>;
1064 return std::forward<decltype(x)>(x);
1070 using AWrap = decltype(a_wrap);
1072 using BWrap = decltype(b_wrap);
1073 using CW = Eigen::CwiseBinaryOp<decltype(op), std::remove_reference_t<AWrap>, std::remove_reference_t<BWrap>>;
1074 CW s {std::forward<AWrap>(a_wrap), std::forward<BWrap>(b_wrap), std::move(op)};
1076 if constexpr (triangle_type_of_v<A, B> !=
TriangleType::any)
return OpenKalman::make_triangular_matrix<triangle_type_of_v<A, B>>(std::move(s));
1077 else if constexpr (hermitian_matrix<A> and hermitian_matrix<B>)
return make_hermitian_matrix<h>(std::move(s));
1083 template<
typename Arg,
typename S,
typename Op>
1084 static constexpr
auto 1085 scalar_op_impl(Arg&& arg, S&& s, Op&& op)
1088 using ConstOp = Eigen::internal::scalar_constant_op<Scalar>;
1090 using M = decltype(m);
1092 using Index =
typename PlainObjectType::Index;
1093 auto c {Eigen::CwiseNullaryOp<ConstOp, PlainObjectType> {
1094 static_cast<Index
>(get_index_dimension_of<0>(m)),
1095 static_cast<Index
>(get_index_dimension_of<1>(m)),
1097 using CW = Eigen::CwiseBinaryOp<Op, std::remove_reference_t<M>, decltype(c)>;
1098 return CW {std::forward<M>(m), c, std::forward<Op>(op)};
1103 #ifdef __cpp_concepts 1104 template<indexible Arg, values::scalar S>
1105 static constexpr vector_space_descriptors_may_match_with<Arg>
auto 1107 template<
typename Arg,
typename S>
1108 static constexpr
auto 1110 scalar_product(Arg&& arg, S&& s)
1113 using Op = Eigen::internal::scalar_product_op<Scalar, Scalar>;
1114 return scalar_op_impl(std::forward<Arg>(arg), std::forward<S>(s), Op{});
1118 #ifdef __cpp_concepts 1119 template<indexible Arg, values::scalar S>
1120 static constexpr vector_space_descriptors_may_match_with<Arg>
auto 1122 template<
typename Arg,
typename S>
1123 static constexpr
auto 1125 scalar_quotient(Arg&& arg, S&& s)
1128 using Op = Eigen::internal::scalar_quotient_op<Scalar, Scalar>;
1129 return scalar_op_impl(std::forward<Arg>(arg), std::forward<S>(s), Op{});
1133 #ifdef __cpp_concepts 1134 template<Eigen3::eigen_general A, Eigen3::eigen_general B>
1136 template<
typename A,
typename B, std::enable_if_t<(Eigen3::eigen_general<A> and Eigen3::eigen_general<B>),
int> = 0>
1138 static constexpr
auto 1141 if constexpr (diagonal_matrix<A> and not Eigen3::eigen_DiagonalWrapper<A> and not Eigen3::eigen_DiagonalMatrix<A>)
1145 else if constexpr (diagonal_matrix<B> and not Eigen3::eigen_DiagonalWrapper<B> and not Eigen3::eigen_DiagonalMatrix<B>)
1149 else if constexpr (diagonal_matrix<A> or diagonal_matrix<B>)
1152 using AWrap = decltype(a_wrap);
1154 using BWrap = decltype(b_wrap);
1155 using Prod = Eigen::Product<std::remove_reference_t<AWrap>, std::remove_reference_t<BWrap>, Eigen::LazyProduct>;
1156 return Prod {std::forward<AWrap>(a_wrap), std::forward<BWrap>(b_wrap)};
1158 else if constexpr (not Eigen3::eigen_matrix_general<B>)
1162 else if constexpr (Eigen3::eigen_matrix_general<A> or Eigen3::eigen_TriangularView<A> or Eigen3::eigen_SelfAdjointView<A>)
1164 using Prod = Eigen::Product<std::remove_reference_t<A>, std::remove_reference_t<B>, Eigen::DefaultProduct>;
1165 return to_dense_object(Prod {std::forward<A>(a), std::forward<B>(b)});
1174 #ifdef __cpp_concepts 1175 template<
bool on_the_right, writable A, indexible B> requires Eigen3::eigen_dense_general<A> or
1176 Eigen3::eigen_DiagonalMatrix<A> or Eigen3::eigen_DiagonalWrapper<A>
1178 template<
bool on_the_right,
typename A,
typename B, std::enable_if_t<writable<A> and
1179 (Eigen3::eigen_dense_general<A> or Eigen3::eigen_DiagonalMatrix<A> or (Eigen3::eigen_DiagonalWrapper<A> and diagonal_adapter<A, 0>)),
int> = 0>
1184 if constexpr (Eigen3::eigen_DiagonalMatrix<A> or Eigen3::eigen_DiagonalWrapper<A>)
1186 static_assert(diagonal_matrix<B>);
1189 else if constexpr (Eigen3::eigen_TriangularView<A>)
1196 auto&& ma = [](
A& a) -> decltype(
auto) {
1197 if constexpr (Eigen3::eigen_array_general<A>)
return a.matrix();
1201 if constexpr (on_the_right)
1202 return ma.applyOnTheRight(OpenKalman::to_native_matrix<T>(std::forward<B>(b)));
1204 return ma.applyOnTheLeft(OpenKalman::to_native_matrix<T>(std::forward<B>(b)));
1210 #ifdef __cpp_concepts 1211 template<TriangleType triangle_type, Eigen3::eigen_SelfAdjo
intView A>
1213 template<TriangleType triangle_type,
typename A, std::enable_if_t<Eigen3::eigen_SelfAdjo
intView<A>,
int> = 0>
1215 static constexpr
auto 1218 using NestedMatrix = std::decay_t<nested_object_of_t<A>>;
1222 if constexpr (constant_matrix<NestedMatrix>)
1231 throw (std::runtime_error(
"cholesky_factor of constant HermitianAdapter: result is indefinite"));
1236 static_assert(diagonal_matrix<A>);
1241 auto euc_dim = get_dimension(dim);
1242 auto col0 = make_constant<A>(
values::sqrt(s), euc_dim, Dimensions<1>{});
1243 auto othercols = make_zero<A>(euc_dim, euc_dim - Dimensions<1>{});
1249 auto euc_dim = get_dimension(dim);
1250 auto row0 = make_constant<A>(
values::sqrt(s), Dimensions<1>{}, dim);
1251 auto otherrows = make_zero<A>(euc_dim - Dimensions<1>{}, euc_dim);
1260 auto LL_x = a.llt();
1261 if (LL_x.info() == Eigen::Success)
1266 b = std::move(LL_x.matrixLLT());
1270 constexpr
unsigned int uplo = triangle_type ==
TriangleType::upper ? Eigen::Upper : Eigen::Lower;
1271 b.template triangularView<uplo>() = LL_x.matrixLLT().adjoint();
1278 if ((not LDL_x.isPositive() and not LDL_x.isNegative()) or LDL_x.info() != Eigen::Success) [[unlikely]]
1280 if (LDL_x.isPositive() and LDL_x.isNegative())
1289 throw (std::runtime_error(
"cholesky_factor of HermitianAdapter: covariance is indefinite"));
1294 b.template triangularView<Eigen::Lower>() =
1295 LDL_x.matrixL().toDenseMatrix() * LDL_x.vectorD().cwiseSqrt().asDiagonal();
1299 b.template triangularView<Eigen::Upper>() =
1300 LDL_x.vectorD().cwiseSqrt().asDiagonal() * LDL_x.matrixU().toDenseMatrix();
1308 template<HermitianAdapterType significant_triangle,
typename A,
typename U,
typename Alpha>
1309 static decltype(
auto)
1312 if constexpr (Eigen3::eigen_matrix_general<A>)
1314 static_assert(writable<A>);
1316 a.template selfadjointView<s>().
template rankUpdate(std::forward<U>(u), alpha);
1317 return std::forward<A>(a);
1321 return rank_update_hermitian<significant_triangle>(
to_native_matrix(std::forward<A>(a)), std::forward<U>(u), alpha);
1326 template<TriangleType triangle,
typename A,
typename U,
typename Alpha>
1327 static decltype(
auto)
1330 if constexpr (Eigen3::eigen_matrix_general<A>)
1332 static_assert(writable<A>);
1335 for (std::size_t i = 0; i < get_index_dimension_of<1>(u); i++)
1337 if (Eigen::internal::llt_inplace<Scalar, t>::rankUpdate(a, get_chip<1>(u, i), alpha) >= 0)
1338 throw (std::runtime_error(
"rank_update_triangular: product is not positive definite"));
1340 return std::forward<A>(a);
1344 return rank_update_triangular<triangle>(
to_native_matrix(std::forward<A>(a)), std::forward<U>(u), alpha);
1349 #ifdef __cpp_concepts 1350 template<
bool must_be_unique,
bool must_be_exact,
typename A,
typename B> requires Eigen3::eigen_matrix_general<B>
1352 template<
bool must_be_unique,
bool must_be_exact,
typename A,
typename B, std::enable_if_t<Eigen3::eigen_matrix_general<B>,
int> = 0>
1354 static constexpr
auto 1359 constexpr std::size_t a_rows = dynamic_dimension<A, 0> ? index_dimension_of_v<B, 0> : index_dimension_of_v<A, 0>;
1360 constexpr std::size_t a_cols = index_dimension_of_v<A, 1>;
1361 constexpr std::size_t b_cols = index_dimension_of_v<B, 1>;
1363 if constexpr (Eigen3::eigen_TriangularView<A>)
1365 auto ret {Eigen::Solve {std::forward<A>(a), std::forward<B>(b)}};
1366 if constexpr (std::is_lvalue_reference_v<A> and std::is_lvalue_reference_v<B>)
return ret;
1369 else if constexpr (Eigen3::eigen_SelfAdjointView<A>)
1371 constexpr
auto uplo = hermitian_adapter_type_of_v<A> ==
TriangleType::upper ? Eigen::Upper : Eigen::Lower;
1372 auto v {std::forward<A>(a).
template selfadjointView<uplo>()};
1376 if (llt.info() == Eigen::Success)
1378 ret = Eigen::Solve {llt, std::forward<B>(b)};
1383 auto ldlt {v.ldlt()};
1384 if ((not ldlt.isPositive() and not ldlt.isNegative()) or ldlt.info() != Eigen::Success)
1386 throw (std::runtime_error(
"Eigen solve (hermitian case): A is indefinite"));
1388 ret = Eigen::Solve {ldlt, std::forward<B>(b)};
1392 else if constexpr (Eigen3::eigen_matrix_general<A>)
1395 if constexpr (must_be_exact or must_be_unique)
1397 auto a_cols_rt = get_index_dimension_of<1>(a);
1398 auto qr {Eigen::ColPivHouseholderQR<Mat> {std::forward<A>(a)}};
1400 if constexpr (must_be_unique)
1402 if (qr.rank() < a_cols_rt)
throw std::runtime_error {
"solve function requests a " 1403 "unique solution, but A is rank-deficient, so result X is not unique"};
1406 auto res {
to_dense_object(Eigen::Solve {std::move(qr), std::forward<B>(b)})};
1408 if constexpr (must_be_exact)
1410 bool a_solution_exists = (a * res).isApprox(b, a_cols_rt * std::numeric_limits<
scalar_type_of_t<A>>::epsilon());
1412 if (a_solution_exists)
return res;
1413 else throw std::runtime_error {
"solve function requests an exact solution, " 1414 "but the solution is only an approximation"};
1423 return to_dense_object(Eigen::Solve {Eigen::HouseholderQR<Mat> {std::forward<A>(a)}, std::forward<B>(b)});
1428 return solve<must_be_unique, must_be_exact>(
to_native_matrix(std::forward<A>(a)), std::forward<B>(b));
1434 template<
typename A>
1435 static constexpr
auto 1436 QR_decomp_impl(
A&& a)
1439 constexpr
auto rows = index_dimension_of_v<A, 0>;
1440 constexpr
auto cols = index_dimension_of_v<A, 1>;
1444 Eigen::HouseholderQR<MatrixType> QR {std::forward<A>(a)};
1446 if constexpr (dynamic_dimension<A, 1>)
1448 auto rt_cols = get_index_dimension_of<1>(a);
1450 ResultType ret {rt_cols, rt_cols};
1452 if constexpr (dynamic_dimension<A, 0>)
1454 auto rt_rows = get_index_dimension_of<0>(a);
1456 if (rt_rows < rt_cols)
1457 ret << QR.matrixQR().topRows(rt_rows),
1460 ret = QR.matrixQR().topRows(rt_cols);
1465 ret << QR.matrixQR().template topRows<rows>(),
1468 ret = QR.matrixQR().topRows(rt_cols);
1477 if constexpr (dynamic_dimension<A, 0>)
1479 auto rt_rows = get_index_dimension_of<0>(a);
1482 ret << QR.matrixQR().topRows(rt_rows),
1485 ret = QR.matrixQR().template topRows<cols>();
1489 if constexpr (rows < cols)
1492 ret = QR.matrixQR().template topRows<cols>();
1501 template<
typename A>
1502 static constexpr
auto 1509 template<
typename A>
1510 static constexpr
auto 1513 return OpenKalman::make_triangular_matrix<TriangleType::upper>(QR_decomp_impl(std::forward<A>(a)));
1521 #endif //OPENKALMAN_EIGEN_LIBRARY_INTERFACE_HPP decltype(auto) constexpr concatenate_vertical(V &&v, Vs &&... vs)
Concatenate one or more typed matrices objects vertically.
Definition: typed-matrix-overloads.hpp:270
typename nested_object_of< T >::type nested_object_of_t
Helper type for nested_object_of.
Definition: nested_object_of.hpp:66
constexpr auto n_ary_operation(const std::tuple< Ds... > &d_tup, Operation &&operation, Args &&...args)
Perform a component-wise n-ary operation, using broadcasting to match the size of a pattern matrix...
Definition: n_ary_operation.hpp:319
decltype(auto) constexpr contract(A &&a, B &&b)
Matrix multiplication of A * B.
Definition: contract.hpp:54
TriangleType
The type of a triangular matrix.
Definition: global-definitions.hpp:60
Definition: indexible_object_traits.hpp:36
Definition: basics.hpp:41
constexpr bool hermitian_adapter
Specifies that a type is a hermitian matrix adapter of a particular type.
Definition: hermitian_adapter.hpp:54
decltype(auto) rank_update_hermitian(A &&a, U &&u, scalar_type_of_t< A > alpha=1)
Do a rank update on a hermitian matrix.
Definition: rank_update_hermitian.hpp:45
Row-major storage (C or C++ style): contiguous storage in which the right-most index has a stride of ...
Definition: eigen-forward-declarations.hpp:90
constexpr bool constant_matrix
Specifies that all components of an object are the same constant value.
Definition: constant_matrix.hpp:31
An operation involving some number of values.
Definition: operation.hpp:69
No storage layout (e.g., if the elements are calculated rather than stored).
decltype(auto) constexpr conjugate(Arg &&arg)
Take the conjugate of a matrix.
Definition: conjugate.hpp:33
auto make_vector_space_adapter(Arg &&arg, Descriptors &&descriptors)
If necessary, wrap an object in a wrapper that adds vector space descriptors for each index...
Definition: make_vector_space_adapter.hpp:37
decltype(auto) constexpr get_slice(Arg &&arg, const std::tuple< Offset... > &offsets, const std::tuple< Extent... > &extents)
Extract a slice from a matrix or tensor.
Definition: get_slice.hpp:101
Forward declarations for OpenKalman's Eigen interface.
typename scalar_type_of< T >::type scalar_type_of_t
helper template for scalar_type_of.
Definition: scalar_type_of.hpp:54
A triangular_adapter, where components above or below the diagonal (or both) are zero.
Definition: forward-class-declarations.hpp:257
Lower, upper, or diagonal matrix.
constexpr bool indexible
T is a generalized tensor type.
Definition: indexible.hpp:32
decltype(auto) to_native_matrix(Arg &&arg)
If it isn't already, convert Arg to a native object in the library associated with LibraryObject...
Definition: to_native_matrix.hpp:35
decltype(auto) constexpr QR_decomposition(A &&a)
Perform a QR decomposition of matrix A=Q[U,0], U is a upper-triangular matrix, and Q is orthogonal...
Definition: QR_decomposition.hpp:33
constexpr auto is_square_shaped(const T &t)
Determine whether an object is square_shaped at runtime.
Definition: is_square_shaped.hpp:63
constexpr bool number
T is a numerical type.
Definition: number.hpp:33
Definition: tuple_reverse.hpp:103
HermitianAdapterType
The type of a hermitian adapter, indicating which triangle of the nested matrix is used...
Definition: global-definitions.hpp:78
constexpr bool triangular_adapter
Specifies that a type is a triangular adapter of triangle type triangle_type.
Definition: triangular_adapter.hpp:45
constexpr bool value
T is numerical value or is reducible to a numerical value.
Definition: value.hpp:31
decltype(auto) constexpr to_diagonal(Arg &&arg)
Convert an indexible object into a diagonal matrix.
Definition: to_diagonal.hpp:32
decltype(auto) constexpr concatenate_horizontal(V &&v, Vs &&... vs)
Concatenate one or more matrix objects vertically.
Definition: typed-matrix-overloads.hpp:308
decltype(auto) constexpr to_dense_object(Arg &&arg)
Convert the argument to a dense, writable matrix of a particular scalar type.
Definition: to_dense_object.hpp:37
decltype(auto) constexpr reduce(BinaryFunction &&b, Arg &&arg)
Perform a partial reduction based on an associative binary function, across one or more indices...
Definition: reduce.hpp:143
An upper-right triangular matrix.
std::conditional_t< sizeof...(dims)==1, Eigen::Matrix< Scalar, detail::eigen_index_convert< dims >..., detail::eigen_index_convert< 1 > >, Eigen::Matrix< Scalar, detail::eigen_index_convert< dims >... > > eigen_matrix_t
An alias for a self-contained, writable, native Eigen matrix.
Definition: eigen-forward-declarations.hpp:491
decltype(auto) constexpr broadcast(Arg &&arg, const Factors &...factors)
Broadcast an object by replicating it by factors specified for each index.
Definition: broadcast.hpp:49
constexpr auto to_number(Arg arg)
Convert any values::value to a values::number.
Definition: to_number.hpp:34
Type scalar type (e.g., std::float, std::double, std::complex<double>) of a tensor.
Definition: scalar_type_of.hpp:32
constexpr bool triangular_matrix
Specifies that a type is a triangular matrix (upper, lower, or diagonal).
Definition: triangular_matrix.hpp:37
The constant associated with T, assuming T is a constant_matrix.
Definition: constant_coefficient.hpp:36
decltype(auto) constexpr transpose(Arg &&arg)
Take the transpose of a matrix.
Definition: transpose.hpp:58
An interface to various routines from the linear algebra library associated with indexible object T...
Definition: library_interface.hpp:37
Column-major storage (Fortran, Matlab, or Eigen style): contiguous storage in which the left-most ext...
constexpr bool identity_matrix
Specifies that a type is an identity matrix.
Definition: identity_matrix.hpp:45
constexpr A && contract_in_place(A &&a, B &&b)
In-place matrix multiplication of A * B, storing the result in A.
Definition: contract_in_place.hpp:38
constexpr auto solve(A &&a, B &&b)
Solve the equation AX = B for X, which may or may not be a unique solution.
Definition: solve.hpp:87
constexpr auto conj(const Arg &arg)
A constexpr function for the complex conjugate of a (complex) number.
Definition: conj.hpp:39
constexpr auto sqrt(const Arg &arg)
A constexpr alternative to std::sqrt.
Definition: sqrt.hpp:46
constexpr auto determinant(Arg &&arg)
Take the determinant of a matrix.
Definition: determinant.hpp:44
decltype(auto) constexpr diagonal_of(Arg &&arg)
Extract a column vector (or column slice for rank>2 tensors) comprising the diagonal elements...
Definition: diagonal_of.hpp:33
A generalization of the above: a custom stride is specified for each index.
constexpr bool size
T is either an index representing a size, or void which represents that there is no size...
Definition: size.hpp:32
decltype(auto) constexpr adjoint(Arg &&arg)
Take the adjoint of a matrix.
Definition: adjoint.hpp:33
Arg && fill_components(Arg &&arg, S...s)
This is an overloaded member function, provided for convenience. It differs from the above function o...
Definition: fill_components.hpp:44
Layout
The layout format of a multidimensional array.
Definition: global-definitions.hpp:47
decltype(auto) constexpr sum(Ts &&...ts)
Element-by-element sum of one or more objects.
Definition: sum.hpp:112
decltype(auto) constexpr LQ_decomposition(A &&a)
Perform an LQ decomposition of matrix A=[L,0]Q, L is a lower-triangular matrix, and Q is orthogonal...
Definition: LQ_decomposition.hpp:33
decltype(auto) rank_update_triangular(A &&a, U &&u, scalar_type_of_t< A > alpha=1)
Do a rank update on triangular matrix.
Definition: rank_update_triangular.hpp:48
std::decay_t< decltype(make_dense_object< T, layout, S >(std::declval< D >()))> dense_writable_matrix_t
An alias for a dense, writable matrix, patterned on parameter T.
Definition: dense_writable_matrix_t.hpp:38
A structure representing the dimensions associated with of a particular index.
Definition: Dimensions.hpp:42
A diagonal matrix (both a lower-left and an upper-right triangular matrix).
constexpr bool is_vector(const T &t)
Return true if T is a vector at runtime.
Definition: is_vector.hpp:44
constexpr auto make_zero(Descriptors &&descriptors)
Make a zero associated with a particular library.
Definition: make_zero.hpp:36
constexpr To && assign(To &&a, From &&b)
Assign a writable object from an indexible object.
Definition: assign.hpp:51
decltype(auto) constexpr get_component(Arg &&arg, const Indices &indices)
Get a component of an object at a particular set of indices.
Definition: get_component.hpp:54
operation(const Operation &, const Args &...) -> operation< Operation, Args... >
Deduction guide.
constexpr std::size_t dynamic_size
A constant indicating that a size or index is dynamic.
Definition: global-definitions.hpp:33
decltype(auto) constexpr nested_object(Arg &&arg)
Retrieve a nested object of Arg, if it exists.
Definition: nested_object.hpp:34
A lower-left triangular matrix.
constexpr bool vector
T is a vector (e.g., column or row vector).
Definition: vector.hpp:67
A matrix with typed rows and columns.
Definition: forward-class-declarations.hpp:448
decltype(auto) constexpr make_triangular_matrix(Arg &&arg)
Create a triangular_matrix from a general matrix.
Definition: make_triangular_matrix.hpp:35
decltype(auto) constexpr cholesky_factor(A &&a)
Take the Cholesky factor of a matrix.
Definition: cholesky_factor.hpp:38
constexpr Arg && set_slice(Arg &&arg, Block &&block, const Begin &...begin)
Assign an object to a particular slice of a matrix or tensor.
Definition: set_slice.hpp:56