OpenKalman
eigen-library-interface.hpp
Go to the documentation of this file.
1 /* This file is part of OpenKalman, a header-only C++ library for
2  * Kalman filters and other recursive filters.
3  *
4  * Copyright (c) 2019-2024 Christopher Lee Ogden <ogden@gatech.edu>
5  *
6  * This Source Code Form is subject to the terms of the Mozilla Public
7  * License, v. 2.0. If a copy of the MPL was not distributed with this
8  * file, You can obtain one at https://mozilla.org/MPL/2.0/.
9  */
10 
16 #ifndef OPENKALMAN_EIGEN_LIBRARY_INTERFACE_HPP
17 #define OPENKALMAN_EIGEN_LIBRARY_INTERFACE_HPP
18 
19 #include <type_traits>
20 #include <tuple>
21 #include <random>
23 
24 namespace OpenKalman::interface
25 {
26 #ifdef __cpp_concepts
27  template<Eigen3::eigen_general<true> T>
28  struct library_interface<T>
29 #else
30  template<typename T>
31  struct library_interface<T, std::enable_if_t<Eigen3::eigen_general<T, true>>>
32 #endif
33  {
34  template<typename Derived>
35  using LibraryBase = Eigen3::EigenAdapterBase<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>>>>;
38 
39  private:
40 
41  template<typename Arg, std::size_t...I, typename...Ind>
42  static constexpr void
43  check_index_bounds(const Arg& arg, std::index_sequence<I...>, Ind...ind)
44  {
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),...);
49  }
50 
51 
52  template<typename Arg>
53  static constexpr decltype(auto)
54  get_coeff(Arg&& arg, Eigen::Index i, Eigen::Index j)
55  {
57  //check_index_bounds(arg, std::make_index_sequence<2> {}, i, j);
58  using Scalar = scalar_type_of_t<Arg>;
59 
60  if constexpr (Eigen3::eigen_DiagonalMatrix<Arg> or Eigen3::eigen_DiagonalWrapper<Arg>)
61  {
62  if (i == j)
63  return static_cast<Scalar>(get_coeff(nested_object(std::forward<Arg>(arg)), i, 0));
64  else
65  return static_cast<Scalar>(0);
66  }
67  else if constexpr (Eigen3::eigen_TriangularView<Arg>)
68  {
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);
72  else
73  return static_cast<Scalar>(get_coeff(nested_object(std::forward<Arg>(arg)), i, j));
74  }
75  else if constexpr (Eigen3::eigen_SelfAdjointView<Arg>)
76  {
77  return get_coeff(nested_object(std::forward<Arg>(arg)), i, j);
78  }
79  else
80  {
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);
86  }
87  }
88 
89 
90  template<typename Indices>
91  static constexpr std::tuple<Eigen::Index, Eigen::Index>
92  extract_indices(const Indices& indices)
93  {
94 #ifdef __cpp_lib_ranges
95  namespace ranges = std::ranges;
96 #endif
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");
105  }
106 
107  public:
108 
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)
114 #else
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)
118 #endif
119  get_component(Arg&& arg, const Indices& indices)
120  {
121  using Scalar = scalar_type_of_t<Arg>;
122  auto [i, j] = extract_indices(indices);
123  if constexpr (Eigen3::eigen_SelfAdjointView<Arg>)
124  {
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>)
128  {
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));
131  }
132  else
133  {
134  if (transp) return get_coeff(std::forward<Arg>(arg), j, i);
135  else return get_coeff(std::forward<Arg>(arg), i, j);
136  }
137  }
138  else
139  {
140  return get_coeff(std::forward<Arg>(arg), i, j);
141  }
142  }
143 
144 
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>)
152 #else
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>
158 #endif
159  static void
160  set_component(Arg& arg, const scalar_type_of_t<Arg>& s, const Indices& indices)
161  {
162  using Scalar = scalar_type_of_t<Arg>;
163  auto [i, j] = extract_indices(indices);
164  if constexpr (Eigen3::eigen_SelfAdjointView<Arg>)
165  {
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>)
169  {
170  if (transp) get_coeff(arg, j, i) = values::conj(s);
171  else get_coeff(arg, i, j) = s;
172  }
173  else
174  {
175  if (transp) get_coeff(arg, j, i) = s;
176  else get_coeff(arg, i, j) = s;
177  }
178  }
179  else
180  {
181  get_coeff(arg, i, j) = s;
182  }
183  }
184 
185  private:
186 
187  template<typename Arg>
188  static decltype(auto)
189  wrap_if_nests_by_reference(Arg&& arg)
190  {
191  if constexpr (Eigen3::eigen_general<Arg>)
192  {
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);
196  else
197  return Eigen3::make_eigen_wrapper(std::forward<Arg>(arg));
198  }
199  else
200  {
201  return Eigen3::make_eigen_wrapper(std::forward<Arg>(arg));
202  }
203  }
204 
205  public:
206 
207  template<typename Arg>
208  static decltype(auto)
209  to_native_matrix(Arg&& arg)
210  {
211  if constexpr (Eigen3::eigen_wrapper<Arg>)
212  {
213  return std::forward<Arg>(arg);
214  }
215  else if constexpr (internal::library_wrapper<Arg>)
216  {
217  return to_native_matrix(OpenKalman::nested_object(std::forward<Arg>(arg)));
218  }
219  else if constexpr (Eigen3::eigen_ArrayWrapper<Arg>)
220  {
221  return to_native_matrix(nested_object(std::forward<Arg>(arg)));
222  }
223  else if constexpr (Eigen3::eigen_array_general<Arg>)
224  {
225  return wrap_if_nests_by_reference(std::forward<Arg>(arg)).matrix();
226  }
227  else if constexpr (Eigen3::eigen_matrix_general<Arg>)
228  {
229  return wrap_if_nests_by_reference(std::forward<Arg>(arg));
230  }
231  else if constexpr (not Eigen3::eigen_general<Arg> and directly_accessible<Arg> and std::is_lvalue_reference_v<Arg>)
232  {
233  using Scalar = scalar_type_of_t<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>;
236 
237  if constexpr (layout_of_v<Arg> == Layout::stride)
238  {
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;
245  }();
246  constexpr int es1 = []() -> int {
247  if constexpr (values::fixed<S1>) return static_cast<std::ptrdiff_t>(S1{});
248  else return Eigen::Dynamic;
249  }();
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));
253 
254  if constexpr (values::fixed<S0> and
255  (es0 == 1 or (values::fixed<S1> and es0 < es1)))
256  {
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};
260  }
261  else if constexpr (values::fixed<S1> and
262  (es1 == 1 or (values::fixed<S0> and es1 < es0)))
263  {
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};
267  }
268  else
269  {
270  if (is1 < is0)
271  {
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};
275  }
276  else
277  {
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};
281  }
282  }
283  }
284  else
285  {
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};
289  }
290  }
291  else
292  {
293  return Eigen3::make_eigen_wrapper(std::forward<Arg>(arg));
294  }
295  }
296 
297  private:
298 
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>>;
302 
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>;
307 
308  public:
309 
310 #ifdef __cpp_concepts
311  template<typename To, Eigen3::eigen_general From> requires (std::assignable_from<To&, From&&>)
312 #else
313  template<typename To, typename From, std::enable_if_t<Eigen3::eigen_general<From> and std::is_assignable_v<To&, From&&>, int> = 0>
314 #endif
315  static void assign(To& a, From&& b)
316  {
317  if constexpr (Eigen3::eigen_DiagonalWrapper<From>)
318  {
319  if constexpr (vector<nested_object_of_t<From>>)
320  a = std::forward<From>(b);
321  else
322  a = diagonal_of(std::forward<From>(b)).asDiagonal();
323  }
324  else
325  {
326  a = std::forward<From>(b);
327  }
328  }
329 
330  private:
331 
332  template<typename Descriptors>
333  static constexpr decltype(auto)
334  extract_descriptors(Descriptors&& descriptors)
335  {
336  if constexpr (pattern_tuple<Descriptors>)
337  {
338  constexpr auto dim = std::tuple_size_v<std::decay_t<Descriptors>>;
339  static_assert(dim <= 2);
340  if constexpr (dim == 0) return std::tuple {coordinates::Axis{}, coordinates::Axis{}};
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);
343  }
344  else
345  {
346 #ifdef __cpp_lib_ranges
347  namespace ranges = std::ranges;
348 #endif
349  auto it = ranges::begin(descriptors);
350  auto e = ranges::end(descriptors);
351  if (it == e) return std::tuple {coordinates::Axis{}, coordinates::Axis{}};
352  auto i = *it;
353  if (++it == e) return std::tuple {i, coordinates::Axis{}};
354  auto j = *it;
355  if (++it == e) return std::tuple {i, j};
356  throw std::logic_error("Wrong number of vector space descriptors");
357  }
358  }
359 
360  public:
361 
362 #ifdef __cpp_concepts
363  template<Layout layout, values::number Scalar, coordinates::euclidean_pattern_collection Ds>
364 #else
365  template<Layout layout, typename Scalar, typename Ds, std::enable_if_t<values::number<Scalar> and
366  coordinates::euclidean_pattern_collection<Ds>, int> = 0>
367 #endif
368  static auto make_default(Ds&& ds)
369  {
370  using IndexType = typename Eigen::Index;
371 
372  if constexpr (pattern_tuple<Ds> and collections::size_of_v<Ds> > 2)
373  {
374  constexpr auto options = layout == Layout::right ? Eigen::RowMajor : Eigen::ColMajor;
375 
376  if constexpr (fixed_pattern_tuple<Ds>)
377  {
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));
381 
382  return Eigen::TensorFixedSize<Scalar, decltype(sizes), options, IndexType> {};
383  }
384  else
385  {
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));
390  }
391  }
392  else
393  {
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)>;
397 
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.");
402 
403  constexpr auto options =
404  layout == Layout::right or (layout == Layout::none and dim0 == 1 and dim1 != 1) ?
405  Eigen::RowMajor : Eigen::ColMajor;
406 
407  using M = writable_type<Scalar, dim0, dim1, options>;
408 
409  if constexpr (dim0 == dynamic_size or dim1 == dynamic_size)
410  return M {static_cast<IndexType>(get_dimension(d0)), static_cast<IndexType>(get_dimension(d1))};
411  else
412  return M {};
413  }
414 
415 
416  }
417 
418 
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>
421  requires (layout == Layout::right) or (layout == Layout::left)
422 #else
423  template<Layout layout, typename Arg, typename S, typename...Ss, std::enable_if_t<writable<Arg> and
424  (layout == Layout::right or layout == Layout::left) 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>
427 #endif
428  static void fill_components(Arg& arg, const S s, const Ss ... ss)
429  {
430  if constexpr (layout == Layout::left) ((arg.transpose() << s), ... , ss);
431  else ((arg << s), ... , ss);
432  }
433 
434 
435 #ifdef __cpp_lib_ranges
436  template<values::dynamic C, coordinates::euclidean_pattern_collection Ds> requires
437  (collections::size_of_v<Ds> != dynamic_size) and (collections::size_of_v<Ds> <= 2)
438  static constexpr constant_matrix auto
439 #else
440  template<typename C, typename Ds, std::enable_if_t<values::dynamic<C> and
441  coordinates::euclidean_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
444 #endif
445  make_constant(C&& c, Ds&& ds)
446  {
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)>;
450 
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>;
455 
456  using IndexType = typename Eigen::Index;
457 
458  if constexpr (fixed_pattern_collection<Ds>)
459  return M::Constant(value);
460  else
461  return M::Constant(static_cast<IndexType>(dim0), static_cast<IndexType>(get_dimension(d1)), value);
462  // Note: We don't address orders greater than 2 because Eigen::Tensor's constant has substantial limitations.
463  }
464 
465 
466 #ifdef __cpp_concepts
467  template<typename Scalar, coordinates::euclidean_pattern_collection Ds> requires
468  (collections::size_of_v<Ds> != dynamic_size) and (collections::size_of_v<Ds> <= 2)
469  static constexpr identity_matrix auto
470 #else
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
474 #endif
475  make_identity_matrix(Ds&& ds)
476  {
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)>;
480 
481  constexpr auto options = (dim0 == 1 and dim1 != 1) ? Eigen::RowMajor : Eigen::ColMajor;
482  using M = writable_type<Scalar, dim0, dim1, options>;
483 
484  using IndexType = typename Eigen::Index;
485 
486  if constexpr (fixed_pattern_collection<Ds>)
487  return M::Identity();
488  else
489  return M::Identity(static_cast<IndexType>(get_dimension(d0)), static_cast<IndexType>(get_dimension(d1)));
490  }
491 
492 
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>)
496 #else
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>
499 #endif
500  static constexpr auto
501  make_triangular_matrix(Arg&& arg)
502  {
503  constexpr auto Mode = t == TriangleType::upper ? Eigen::Upper : Eigen::Lower;
504  return arg.template triangularView<Mode>();
505  }
506 
507 
508 #ifdef __cpp_concepts
509  template<HermitianAdapterType t, Eigen3::eigen_dense_general Arg> requires std::is_lvalue_reference_v<Arg>
510 #else
511  template<HermitianAdapterType t, typename Arg, std::enable_if_t<Eigen3::eigen_dense_general<Arg> and std::is_lvalue_reference_v<Arg>, int> = 0>
512 #endif
513  static constexpr auto
514  make_hermitian_adapter(Arg&& arg)
515  {
516  constexpr auto Mode = t == HermitianAdapterType::upper ? Eigen::Upper : Eigen::Lower;
517  return arg.template selfadjointView<Mode>();
518  }
519 
520 
521  // to_euclidean not defined--rely on default
522 
523  // from_euclidean not defined--rely on default
524 
525  // wrap_angles not defined--rely on default
526 
527 
528 #ifdef __cpp_concepts
529  template<typename Arg, typename...Begin, typename...Size> requires (sizeof...(Begin) <= 2)
530 #else
531  template<typename Arg, typename...Begin, typename...Size, std::enable_if_t<(sizeof...(Begin) <= 2), int> = 0>
532 #endif
533  static auto
534  get_slice(Arg&& arg, const std::tuple<Begin...>& begin, const std::tuple<Size...>& size)
535  {
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));
540  }(begin);
541 
542  auto b1 = [](const auto& begin){
543  if constexpr (sizeof...(Begin) < 2) return std::integral_constant<Eigen::Index, 0>{};
544  else
545  {
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));
549  }
550  }(begin);
551 
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));
556  }(size);
557 
558  auto s1 = [](const auto& size){
559  if constexpr (sizeof...(Size) < 2) return std::integral_constant<Eigen::Index, 1>{};
560  else
561  {
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));
565  }
566  }(size);
567 
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);
570 
571  decltype(auto) m = to_native_matrix(std::forward<Arg>(arg));
572  using M = decltype(m);
573 
574  constexpr auto Flags = Eigen::internal::traits<std::remove_reference_t<M>>::Flags;
575 
576  if constexpr (directly_accessible<M> and not (std::is_lvalue_reference_v<M> and static_cast<bool>(Flags & Eigen::NestByRefBit)))
577  {
578  // workaround for Eigen::Block's special treatment of directly accessible nested types:
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)};
583  else
584  return B {std::move(rep), std::move(b0), std::move(b1), std::move(s0), std::move(s1)};
585  }
586  else
587  {
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)};
593  else
594  return B {std::forward<M>(m), std::move(b0), std::move(b1), std::move(s0), std::move(s1)};
595  }
596  }
597 
598 
599 #ifdef __cpp_concepts
600  template<Eigen3::eigen_dense_general Arg, Eigen3::eigen_general Block, typename...Begin> requires (sizeof...(Begin) <= 2)
601 #else
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>
604 #endif
605  static void
606  set_slice(Arg& arg, Block&& block, const Begin&...begin)
607  {
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);
610 
611  if constexpr (Eigen3::eigen_Block<Block>)
612  {
613  if (std::addressof(arg) == std::addressof(block.nestedExpression()) and b0 == block.startRow() and b1 == block.startCol())
614  return;
615  }
616 
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>;
620 
621  if constexpr (not has_dynamic_dimensions<Block>)
622  {
623  Bk {arg, b0, b1} = std::forward<Block>(block);
624  }
625  else
626  {
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);
630  }
631  }
632 
633 
634 #ifdef __cpp_concepts
635  template<TriangleType t, Eigen3::eigen_SelfAdjointView 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&&>
637 #else
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>
641 #endif
642  static void
643  set_triangle(A&& a, B&& b)
644  {
645  // A SelfAdjointView won't always be a hermitian_adapter
648  internal::set_triangle<t>(OpenKalman::nested_object(std::forward<A>(a)), OpenKalman::adjoint(std::forward<B>(b)));
649  else
650  internal::set_triangle<t>(OpenKalman::nested_object(std::forward<A>(a)), std::forward<B>(b));
651  }
652 
653 
654 #ifdef __cpp_concepts
655  template<TriangleType t, typename A, Eigen3::eigen_general B> requires
656  (diagonal_adapter<A> or t == TriangleType::diagonal) and
657  std::assignable_from<decltype(OpenKalman::diagonal_of(std::declval<A&&>())), decltype(OpenKalman::diagonal_of(std::declval<B&&>()))>
658 #else
659  template<TriangleType t, typename A, typename B, std::enable_if_t<
660  Eigen3::eigen_general<B> and (diagonal_adapter<A> or t == TriangleType::diagonal) and
661  std::is_assignable_v<decltype(OpenKalman::diagonal_of(std::declval<A&&>())), decltype(OpenKalman::diagonal_of(std::declval<B&&>()))>, int> = 0>
662 #endif
663  static void
664  set_triangle(A&& a, B&& b)
665  {
666  assign(OpenKalman::diagonal_of(std::forward<A>(a)), OpenKalman::diagonal_of(std::forward<B>(b)));
667  }
668 
669 
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&&>
674 #else
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>
678 #endif
679  static void
680  set_triangle(A&& a, B&& b)
681  {
682  internal::set_triangle<t>(OpenKalman::nested_object(std::forward<A>(a)), std::forward<B>(b));
683  }
684 
685 
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
689  writable<A&&> and (t != TriangleType::diagonal)
690 #else
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
694  writable<A&&> and (t != TriangleType::diagonal), int> = 0>
695 #endif
696  static void
697  set_triangle(A&& a, B&& b)
698  {
699  a.template triangularView<t == TriangleType::upper ? Eigen::Upper : Eigen::Lower>() = std::forward<B>(b);
700  }
701 
702 
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>)
706 #else
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>
710 #endif
711  static constexpr auto
712  to_diagonal(Arg&& arg)
713  {
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))};
716 
717  if constexpr (square_shaped<Arg> or dimension_size_of_index_is<Arg, 0, 1>)
718  {
719  return internal::make_fixed_size_adapter(std::forward<Arg>(arg)); // Make one-dimensional matrix
720  }
721  else if constexpr (Eigen3::eigen_array_general<Arg>)
722  {
723  return arg.matrix().asDiagonal();
724  }
725  else
726  {
727  return arg.asDiagonal();
728  }
729  }
730 
731 
732 #ifdef __cpp_concepts
733  template<Eigen3::eigen_SelfAdjointView Arg> requires (hermitian_adapter_type_of_v<Arg> == HermitianAdapterType::lower)
734 #else
735  template<typename Arg, std::enable_if_t<Eigen3::eigen_SelfAdjointView<Arg>, int> = 0>
736 #endif
737  static constexpr auto
738  to_diagonal(Arg&& arg)
739  {
740  // If it is a column vector, a "lower SelfAdjointView wrapper is irrelevant
741  return OpenKalman::to_diagonal(nested_object(std::forward<Arg>(arg)));
742  }
743 
744 
745 #ifdef __cpp_concepts
746  template<Eigen3::eigen_TriangularView Arg> requires (triangle_type_of_v<Arg> == TriangleType::lower)
747 #else
748  template<typename Arg, std::enable_if_t<Eigen3::eigen_TriangularView<Arg> and (triangle_type_of_v<Arg> == TriangleType::lower), int> = 0>
749 #endif
750  static constexpr auto
751  to_diagonal(Arg&& arg)
752  {
753  // If it is a column vector, a "lower" TriangularView wrapper is irrelevant
754  return OpenKalman::to_diagonal(nested_object(std::forward<Arg>(arg)));
755  }
756 
757 
758  template<typename Arg>
759 #ifdef __cpp_concepts
760  static constexpr vector decltype(auto)
761 #else
762  static constexpr decltype(auto)
763 #endif
764  diagonal_of(Arg&& arg)
765  {
766  using Scalar = scalar_type_of_t<Arg>;
767 
768  if constexpr (Eigen3::eigen_DiagonalWrapper<Arg>)
769  {
770  using Diag = decltype(nested_object(std::forward<Arg>(arg))); //< must be nested_object(...) rather than .diagonal() because of const_cast
771  using EigenTraits = Eigen::internal::traits<std::decay_t<Diag>>;
772  constexpr auto rows = EigenTraits::RowsAtCompileTime;
773  constexpr auto cols = EigenTraits::ColsAtCompileTime;
774 
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)
778  {
779  auto ret {nested_object(std::forward<Arg>(arg))};
780  return ret;
781  }
782  else if constexpr (rows == 1 or rows == 0)
783  {
784  auto ret {OpenKalman::transpose(nested_object(std::forward<Arg>(arg)))};
785  return ret;
786  }
787  else if constexpr (rows == Eigen::Dynamic or cols == Eigen::Dynamic)
788  {
789  decltype(auto) diag = nested_object(std::forward<Arg>(arg));
790  using M = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
791  return M {M::Map(to_dense_object(std::forward<Diag>(diag)).data(),
792  get_index_dimension_of<0>(diag) * get_index_dimension_of<1>(diag))};
793  }
794  else // rows > 1 and cols > 1
795  {
796  using M = Eigen::Matrix<Scalar, rows * cols, 1>;
797  return M {M::Map(to_dense_object(nested_object(std::forward<Arg>(arg))).data())};
798  }
799  }
800  else if constexpr (Eigen3::eigen_SelfAdjointView<Arg> or Eigen3::eigen_TriangularView<Arg>)
801  {
802  return OpenKalman::diagonal_of(nested_object(std::forward<Arg>(arg)));
803  }
804  else if constexpr (Eigen3::eigen_Identity<Arg>)
805  {
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);
809  }
810  else if constexpr (Eigen3::eigen_dense_general<Arg>)
811  {
812  decltype(auto) m = to_native_matrix(std::forward<Arg>(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)};
817  }
818  else
819  {
820  return OpenKalman::diagonal_of(to_native_matrix(std::forward<Arg>(arg)));
821  }
822  }
823 
824 
825  template<typename Arg, typename Factor0 = std::integral_constant<std::size_t, 1>, typename Factor1 = std::integral_constant<std::size_t, 1>>
826  static auto
827  broadcast(Arg&& arg, const Factor0& factor0 = Factor0{}, const Factor1& factor1 = Factor1{})
828  {
829  constexpr int F0 = []{
830  if constexpr (values::fixed<Factor0>) return static_cast<std::size_t>(Factor0{});
831  else return Eigen::Dynamic;
832  }();
833  constexpr int F1 = []{
834  if constexpr (values::fixed<Factor1>) return static_cast<std::size_t>(Factor1{});
835  else return Eigen::Dynamic;
836  }();
837 
838  using IndexType = typename std::decay_t<Arg>::Index;
839 
840  decltype(auto) m = to_native_matrix(std::forward<Arg>(arg));
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)};
845  else
846  return R {std::forward<M>(m), static_cast<IndexType>(factor0), static_cast<IndexType>(factor1)};
847  }
848 
849  private:
850 
851  // Only to be used in a non-evaluated context
852  template<typename Op, typename...S>
853  static constexpr auto dummy_op(Op op, S...s)
854  {
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});
858  }
859 
860 
861  template<typename Op, typename...Args>
862  struct EigenNaryOp;
863 
864  template<typename Op, typename Arg>
865  struct EigenNaryOp<Op, Arg> { using type = Eigen::CwiseUnaryOp<Op, Arg>; };
866 
867  template<typename Op, typename Arg1, typename Arg2>
868  struct EigenNaryOp<Op, Arg1, Arg2> { using type = Eigen::CwiseBinaryOp<Op, Arg1, Arg2>; };
869 
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>; };
872 
873  public:
874 
875 #ifdef __cpp_concepts
876  template<coordinates::pattern...Ds, typename Operation, indexible...Args> requires
877  (sizeof...(Ds) <= 2) and (sizeof...(Args) <= 3) and
878  (values::number<std::invoke_result_t<Operation, scalar_type_of_t<Args>...>> or
879  (sizeof...(Args) == 0 and
880  (values::number<std::invoke_result_t<Operation, std::conditional_t<true, std::size_t, Ds>...>> or
881  values::number<std::invoke_result_t<Operation, std::size_t>>)))
882 #else
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
885  (values::number<typename std::invoke_result<Operation, typename scalar_type_of<Args>::type...>::type> or
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>
889 #endif
890  static auto
891  n_ary_operation(const std::tuple<Ds...>& tup, Operation&& operation, Args&&...args)
892  {
893  decltype(auto) op = Eigen3::native_operation(std::forward<Operation>(operation));
894  using Op = decltype(op);
895  using Scalar = decltype(dummy_op(operation, std::declval<scalar_type_of_t<Args>>()...));
896 
897  if constexpr (sizeof...(Args) == 0)
898  {
899  using P = dense_writable_matrix_t<T, Layout::none, Scalar, std::tuple<Ds...>>;
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)};
904  }
905  else
906  {
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)};
910  }
911  }
912 
913 
914 #ifdef __cpp_concepts
915  template<std::size_t...indices, typename BinaryFunction, typename Arg> requires ((indices <= 1) and ...)
916 #else
917  template<std::size_t...indices, typename BinaryFunction, typename Arg, std::enable_if_t<((indices <= 1) and ...), int> = 0>
918 #endif
919  static auto
920  reduce(BinaryFunction&& b, Arg&& arg)
921  {
922  if constexpr (Eigen3::eigen_dense_general<Arg>)
923  {
924  auto&& op = Eigen3::native_operation(std::forward<BinaryFunction>(b));
925  using Op = decltype(op);
926 
927  if constexpr (((indices == 0) or ...) and ((indices == 1) or ...))
928  {
929  return std::forward<Arg>(arg).redux(std::forward<Op>(op)); // Complete reduction, which will be a scalar.
930  }
931  else
932  {
933  constexpr auto dir = ((indices == 0) and ...) ? Eigen::Vertical : Eigen::Horizontal;
934  using ROp = Eigen::internal::member_redux<std::decay_t<Op>, scalar_type_of_t<Arg>>;
935  decltype(auto) m = to_native_matrix(std::forward<Arg>(arg));
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)}};
940  }
941  }
942  else
943  {
944  return reduce<indices...>(std::forward<BinaryFunction>(b), to_native_matrix(std::forward<Arg>(arg)));
945  }
946  }
947 
948 
949  template<typename Arg>
950  static constexpr decltype(auto)
951  conjugate(Arg&& arg)
952  {
953  if constexpr (Eigen3::eigen_dense_general<Arg>)
954  {
955  using UnaryOp = Eigen::internal::scalar_conjugate_op<scalar_type_of_t<Arg>>;
956  decltype(auto) m = to_native_matrix(std::forward<Arg>(arg));
957  using M = decltype(m);
958  return Eigen::CwiseUnaryOp<UnaryOp, std::remove_reference_t<M>> {std::forward<M>(m)};
959  }
960  else if constexpr (Eigen3::eigen_TriangularView<Arg> or Eigen3::eigen_SelfAdjointView<Arg>)
961  {
962  return std::forward<Arg>(arg).conjugate();
963  }
964  else if constexpr (triangular_matrix<Arg>)
965  {
966  return OpenKalman::conjugate(TriangularAdapter {std::forward<Arg>(arg)});
967  }
968  else
969  {
970  return conjugate(to_native_matrix(std::forward<Arg>(arg)));
971  }
972  // Note: the global conjugate function already handles Eigen::DiagonalMatrix and Eigen::DiagonalWrapper
973  }
974 
975 
976  template<typename Arg>
977  static constexpr decltype(auto)
978  transpose(Arg&& arg)
979  {
980  if constexpr (Eigen3::eigen_matrix_general<Arg>)
981  {
982  decltype(auto) m = to_native_matrix(std::forward<Arg>(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)};
987  }
988  else if constexpr (Eigen3::eigen_TriangularView<Arg> or Eigen3::eigen_SelfAdjointView<Arg>)
989  {
990  return std::forward<Arg>(arg).transpose();
991  }
992  else if constexpr (triangular_matrix<Arg>)
993  {
994  static_assert(triangle_type_of_v<Arg> == TriangleType::upper or triangle_type_of_v<Arg> == TriangleType::lower);
995  constexpr auto t = triangular_matrix<Arg, TriangleType::upper> ? TriangleType::lower : TriangleType::upper;
996  return OpenKalman::make_triangular_matrix<t>(transpose(to_native_matrix(std::forward<Arg>(arg))));
997  }
998  else
999  {
1000  return transpose(to_native_matrix(std::forward<Arg>(arg)));
1001  }
1002  // Note: the global transpose function already handles zero, constant, diagonal, constant-diagonal, and symmetric cases.
1003  }
1004 
1005 
1006 #ifdef __cpp_concepts
1007  template<typename Arg> requires (not Eigen3::eigen_dense_general<Arg>)
1008 #else
1009  template<typename Arg, std::enable_if_t<not Eigen3::eigen_dense_general<Arg>, int> = 0>
1010 #endif
1011  static constexpr decltype(auto)
1012  adjoint(Arg&& arg)
1013  {
1014  if constexpr (Eigen3::eigen_TriangularView<Arg> or Eigen3::eigen_SelfAdjointView<Arg>)
1015  {
1016  return std::forward<Arg>(arg).adjoint();
1017  }
1018  else if constexpr (triangular_matrix<Arg>)
1019  {
1020  static_assert(triangle_type_of_v<Arg> == TriangleType::upper or triangle_type_of_v<Arg> == TriangleType::lower);
1021  constexpr auto t = triangular_matrix<Arg, TriangleType::upper> ? TriangleType::lower : TriangleType::upper;
1022  return OpenKalman::make_triangular_matrix<t>(OpenKalman::adjoint(to_native_matrix(std::forward<Arg>(arg))));
1023  }
1024  else
1025  {
1026  return OpenKalman::adjoint(to_native_matrix(std::forward<Arg>(arg)));
1027  }
1028  // Note: the global adjoint function already handles zero, constant, diagonal, non-complex, and hermitian cases.
1029  }
1030 
1031 
1032  template<typename Arg>
1033  static constexpr auto
1034  determinant(Arg&& arg)
1035  {
1036  if constexpr (Eigen3::eigen_matrix_general<Arg, true>)
1037  return std::forward<Arg>(arg).determinant();
1038  else if constexpr (Eigen3::eigen_array_general<Arg, true>)
1039  return std::forward<Arg>(arg).matrix().determinant();
1040  else
1041  return to_native_matrix(std::forward<Arg>(arg)).determinant();
1042  // Note: the global determinant function already handles Eigen::TriangularView, Eigen::DiagonalMatrix, and Eigen::DiagonalWrapper
1043  }
1044 
1045 
1046 #ifdef __cpp_concepts
1047  template<typename A, Eigen3::eigen_general B>
1048 #else
1049  template<typename A, typename B, std::enable_if_t<Eigen3::eigen_general<B>, int> = 0>
1050 #endif
1051  static constexpr auto
1052  sum(A&& a, B&& b)
1053  {
1054  constexpr HermitianAdapterType h = hermitian_adapter_type_of_v<A, B> == HermitianAdapterType::any ?
1055  HermitianAdapterType::lower : hermitian_adapter_type_of_v<A, B>;
1056 
1057  auto f = [](auto&& x) -> decltype(auto) {
1058  constexpr bool herm = hermitian_matrix<A> and hermitian_matrix<B>;
1059  if constexpr ((triangle_type_of_v<A, B> != TriangleType::any and triangular_adapter<decltype(x)>) or (herm and hermitian_adapter<decltype(x), h>))
1060  return nested_object(std::forward<decltype(x)>(x));
1061  else if constexpr (herm and hermitian_adapter<decltype(x)>)
1062  return OpenKalman::transpose(nested_object(std::forward<decltype(x)>(x)));
1063  else
1064  return std::forward<decltype(x)>(x);
1065  };
1066 
1067  auto op = Eigen::internal::scalar_sum_op<scalar_type_of_t<A>, scalar_type_of_t<B>>{};
1068 
1069  decltype(auto) a_wrap = to_native_matrix(f(std::forward<A>(a)));
1070  using AWrap = decltype(a_wrap);
1071  decltype(auto) b_wrap = to_native_matrix(f(std::forward<B>(b)));
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)};
1075 
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));
1078  else return s;
1079  }
1080 
1081  private:
1082 
1083  template<typename Arg, typename S, typename Op>
1084  static constexpr auto
1085  scalar_op_impl(Arg&& arg, S&& s, Op&& op)
1086  {
1087  using Scalar = scalar_type_of_t<Arg>;
1088  using ConstOp = Eigen::internal::scalar_constant_op<Scalar>;
1089  decltype(auto) m {to_native_matrix(std::forward<Arg>(arg))};
1090  using M = decltype(m);
1091  using PlainObjectType = dense_writable_matrix_t<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)),
1096  ConstOp {static_cast<Scalar>(values::to_number(s))}}};
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)};
1099  }
1100 
1101  public:
1102 
1103 #ifdef __cpp_concepts
1104  template<indexible Arg, values::scalar S>
1105  static constexpr vector_space_descriptors_may_match_with<Arg> auto
1106 #else
1107  template<typename Arg, typename S>
1108  static constexpr auto
1109 #endif
1110  scalar_product(Arg&& arg, S&& s)
1111  {
1112  using Scalar = scalar_type_of_t<Arg>;
1113  using Op = Eigen::internal::scalar_product_op<Scalar, Scalar>;
1114  return scalar_op_impl(std::forward<Arg>(arg), std::forward<S>(s), Op{});
1115  }
1116 
1117 
1118 #ifdef __cpp_concepts
1119  template<indexible Arg, values::scalar S>
1120  static constexpr vector_space_descriptors_may_match_with<Arg> auto
1121 #else
1122  template<typename Arg, typename S>
1123  static constexpr auto
1124 #endif
1125  scalar_quotient(Arg&& arg, S&& s)
1126  {
1127  using Scalar = scalar_type_of_t<Arg>;
1128  using Op = Eigen::internal::scalar_quotient_op<Scalar, Scalar>;
1129  return scalar_op_impl(std::forward<Arg>(arg), std::forward<S>(s), Op{});
1130  }
1131 
1132 
1133 #ifdef __cpp_concepts
1134  template<Eigen3::eigen_general A, Eigen3::eigen_general B>
1135 #else
1136  template<typename A, typename B, std::enable_if_t<(Eigen3::eigen_general<A> and Eigen3::eigen_general<B>), int> = 0>
1137 #endif
1138  static constexpr auto
1139  contract(A&& a, B&& b)
1140  {
1141  if constexpr (diagonal_matrix<A> and not Eigen3::eigen_DiagonalWrapper<A> and not Eigen3::eigen_DiagonalMatrix<A>)
1142  {
1143  return contract(to_native_matrix(diagonal_of(std::forward<A>(a))).asDiagonal(), std::forward<B>(b));
1144  }
1145  else if constexpr (diagonal_matrix<B> and not Eigen3::eigen_DiagonalWrapper<B> and not Eigen3::eigen_DiagonalMatrix<B>)
1146  {
1147  return contract(std::forward<A>(a), to_native_matrix(diagonal_of(std::forward<B>(b))).asDiagonal());
1148  }
1149  else if constexpr (diagonal_matrix<A> or diagonal_matrix<B>)
1150  {
1151  decltype(auto) a_wrap = to_native_matrix(std::forward<A>(a));
1152  using AWrap = decltype(a_wrap);
1153  decltype(auto) b_wrap = to_native_matrix(std::forward<B>(b));
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)};
1157  }
1158  else if constexpr (not Eigen3::eigen_matrix_general<B>)
1159  {
1160  return contract(std::forward<A>(a), to_native_matrix(std::forward<B>(b)));
1161  }
1162  else if constexpr (Eigen3::eigen_matrix_general<A> or Eigen3::eigen_TriangularView<A> or Eigen3::eigen_SelfAdjointView<A>)
1163  {
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)});
1166  }
1167  else
1168  {
1169  return contract(to_native_matrix(std::forward<A>(a)), std::forward<B>(b));
1170  }
1171  }
1172 
1173 
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>
1177 #else
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>
1180 #endif
1181  static A&
1182  contract_in_place(A& a, B&& b)
1183  {
1184  if constexpr (Eigen3::eigen_DiagonalMatrix<A> or Eigen3::eigen_DiagonalWrapper<A>)
1185  {
1186  static_assert(diagonal_matrix<B>);
1187  nested_object(a) = diagonal_of(a).array() * diagonal_of(std::forward<B>(b)).array();
1188  }
1189  else if constexpr (Eigen3::eigen_TriangularView<A>)
1190  {
1191  static_assert(triangular_matrix<B, triangle_type_of_v<A>>);
1192  nested_object(a) = contract(a, std::forward<B>(b));
1193  }
1194  else
1195  {
1196  auto&& ma = [](A& a) -> decltype(auto) {
1197  if constexpr (Eigen3::eigen_array_general<A>) return a.matrix();
1198  else return (a);
1199  }(a);
1200 
1201  if constexpr (on_the_right)
1202  return ma.applyOnTheRight(OpenKalman::to_native_matrix<T>(std::forward<B>(b)));
1203  else
1204  return ma.applyOnTheLeft(OpenKalman::to_native_matrix<T>(std::forward<B>(b)));
1205  return a;
1206  }
1207  }
1208 
1209 
1210 #ifdef __cpp_concepts
1211  template<TriangleType triangle_type, Eigen3::eigen_SelfAdjointView A>
1212 #else
1213  template<TriangleType triangle_type, typename A, std::enable_if_t<Eigen3::eigen_SelfAdjointView<A>, int> = 0>
1214 #endif
1215  static constexpr auto
1216  cholesky_factor(A&& a)
1217  {
1218  using NestedMatrix = std::decay_t<nested_object_of_t<A>>;
1219  using Scalar = scalar_type_of_t<A>;
1220  auto dim = *is_square_shaped(a);
1221 
1222  if constexpr (constant_matrix<NestedMatrix>)
1223  {
1224  // If nested matrix is a positive constant matrix, construct the Cholesky factor using a shortcut.
1225 
1226  auto s = constant_coefficient {a};
1227 
1228  if (values::to_number(s) < Scalar(0))
1229  {
1230  // Cholesky factor elements are complex, so throw an exception.
1231  throw (std::runtime_error("cholesky_factor of constant HermitianAdapter: result is indefinite"));
1232  }
1233 
1234  if constexpr(triangle_type == TriangleType::diagonal)
1235  {
1236  static_assert(diagonal_matrix<A>);
1237  return OpenKalman::to_diagonal(make_constant<A>(values::sqrt(s), dim, Dimensions<1>{}));
1238  }
1239  else if constexpr(triangle_type == TriangleType::lower)
1240  {
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>{});
1244  return make_vector_space_adapter(OpenKalman::make_triangular_matrix<triangle_type>(concatenate_horizontal(col0, othercols)), dim, dim);
1245  }
1246  else
1247  {
1248  static_assert(triangle_type == TriangleType::upper);
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);
1252  return make_vector_space_adapter(OpenKalman::make_triangular_matrix<triangle_type>(concatenate_vertical(row0, otherrows)), dim, dim);
1253  }
1254  }
1255  else
1256  {
1257  // For the general case, perform an LLT Cholesky decomposition.
1258  using M = dense_writable_matrix_t<A>;
1259  M b;
1260  auto LL_x = a.llt();
1261  if (LL_x.info() == Eigen::Success)
1262  {
1263  if constexpr((triangle_type == TriangleType::upper and hermitian_adapter_type_of_v<A> == HermitianAdapterType::upper) or
1264  triangle_type == TriangleType::lower and hermitian_adapter_type_of_v<A> == HermitianAdapterType::lower)
1265  {
1266  b = std::move(LL_x.matrixLLT());
1267  }
1268  else
1269  {
1270  constexpr unsigned int uplo = triangle_type == TriangleType::upper ? Eigen::Upper : Eigen::Lower;
1271  b.template triangularView<uplo>() = LL_x.matrixLLT().adjoint();
1272  }
1273  }
1274  else [[unlikely]]
1275  {
1276  // If covariance is not positive definite, use the more robust LDLT decomposition.
1277  auto LDL_x = nested_object(std::forward<A>(a)).ldlt();
1278  if ((not LDL_x.isPositive() and not LDL_x.isNegative()) or LDL_x.info() != Eigen::Success) [[unlikely]]
1279  {
1280  if (LDL_x.isPositive() and LDL_x.isNegative()) // Covariance is zero, even though decomposition failed.
1281  {
1282  if constexpr(triangle_type == TriangleType::lower)
1283  b.template triangularView<Eigen::Lower>() = make_zero(nested_object(a));
1284  else
1285  b.template triangularView<Eigen::Upper>() = make_zero(nested_object(a));
1286  }
1287  else // Covariance is indefinite, so throw an exception.
1288  {
1289  throw (std::runtime_error("cholesky_factor of HermitianAdapter: covariance is indefinite"));
1290  }
1291  }
1292  else if constexpr(triangle_type == TriangleType::lower)
1293  {
1294  b.template triangularView<Eigen::Lower>() =
1295  LDL_x.matrixL().toDenseMatrix() * LDL_x.vectorD().cwiseSqrt().asDiagonal();
1296  }
1297  else
1298  {
1299  b.template triangularView<Eigen::Upper>() =
1300  LDL_x.vectorD().cwiseSqrt().asDiagonal() * LDL_x.matrixU().toDenseMatrix();
1301  }
1302  }
1303  return TriangularAdapter<M, triangle_type> {std::move(b)};
1304  }
1305  }
1306 
1307 
1308  template<HermitianAdapterType significant_triangle, typename A, typename U, typename Alpha>
1309  static decltype(auto)
1310  rank_update_hermitian(A&& a, U&& u, const Alpha alpha)
1311  {
1312  if constexpr (Eigen3::eigen_matrix_general<A>)
1313  {
1314  static_assert(writable<A>);
1315  constexpr auto s = significant_triangle == HermitianAdapterType::lower ? Eigen::Lower : Eigen::Upper;
1316  a.template selfadjointView<s>().template rankUpdate(std::forward<U>(u), alpha);
1317  return std::forward<A>(a);
1318  }
1319  else
1320  {
1321  return rank_update_hermitian<significant_triangle>(to_native_matrix(std::forward<A>(a)), std::forward<U>(u), alpha);
1322  }
1323  }
1324 
1325 
1326  template<TriangleType triangle, typename A, typename U, typename Alpha>
1327  static decltype(auto)
1328  rank_update_triangular(A&& a, U&& u, const Alpha alpha)
1329  {
1330  if constexpr (Eigen3::eigen_matrix_general<A>)
1331  {
1332  static_assert(writable<A>);
1333  constexpr auto t = triangle == TriangleType::lower ? Eigen::Lower : Eigen::Upper;
1334  using Scalar = scalar_type_of_t<A>;
1335  for (std::size_t i = 0; i < get_index_dimension_of<1>(u); i++)
1336  {
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"));
1339  }
1340  return std::forward<A>(a);
1341  }
1342  else
1343  {
1344  return rank_update_triangular<triangle>(to_native_matrix(std::forward<A>(a)), std::forward<U>(u), alpha);
1345  }
1346  }
1347 
1348 
1349 #ifdef __cpp_concepts
1350  template<bool must_be_unique, bool must_be_exact, typename A, typename B> requires Eigen3::eigen_matrix_general<B>
1351 #else
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>
1353 #endif
1354  static constexpr auto
1355  solve(A&& a, B&& b)
1356  {
1357  using Scalar = scalar_type_of_t<A>;
1358 
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>;
1362 
1363  if constexpr (Eigen3::eigen_TriangularView<A>)
1364  {
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;
1367  else return to_dense_object(std::move(ret));
1368  }
1369  else if constexpr (Eigen3::eigen_SelfAdjointView<A>)
1370  {
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>()};
1373  auto llt {v.llt()};
1374 
1376  if (llt.info() == Eigen::Success)
1377  {
1378  ret = Eigen::Solve {llt, std::forward<B>(b)};
1379  }
1380  else [[unlikely]]
1381  {
1382  // A is semidefinite. Use LDLT decomposition instead.
1383  auto ldlt {v.ldlt()};
1384  if ((not ldlt.isPositive() and not ldlt.isNegative()) or ldlt.info() != Eigen::Success)
1385  {
1386  throw (std::runtime_error("Eigen solve (hermitian case): A is indefinite"));
1387  }
1388  ret = Eigen::Solve {ldlt, std::forward<B>(b)};
1389  }
1390  return ret;
1391  }
1392  else if constexpr (Eigen3::eigen_matrix_general<A>)
1393  {
1395  if constexpr (must_be_exact or must_be_unique)
1396  {
1397  auto a_cols_rt = get_index_dimension_of<1>(a);
1398  auto qr {Eigen::ColPivHouseholderQR<Mat> {std::forward<A>(a)}};
1399 
1400  if constexpr (must_be_unique)
1401  {
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"};
1404  }
1405 
1406  auto res {to_dense_object(Eigen::Solve {std::move(qr), std::forward<B>(b)})};
1407 
1408  if constexpr (must_be_exact)
1409  {
1410  bool a_solution_exists = (a * res).isApprox(b, a_cols_rt * std::numeric_limits<scalar_type_of_t<A>>::epsilon());
1411 
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"};
1415  }
1416  else
1417  {
1418  return res;
1419  }
1420  }
1421  else
1422  {
1423  return to_dense_object(Eigen::Solve {Eigen::HouseholderQR<Mat> {std::forward<A>(a)}, std::forward<B>(b)});
1424  }
1425  }
1426  else
1427  {
1428  return solve<must_be_unique, must_be_exact>(to_native_matrix(std::forward<A>(a)), std::forward<B>(b));
1429  }
1430  }
1431 
1432  private:
1433 
1434  template<typename A>
1435  static constexpr auto
1436  QR_decomp_impl(A&& a)
1437  {
1438  using Scalar = scalar_type_of_t<A>;
1439  constexpr auto rows = index_dimension_of_v<A, 0>;
1440  constexpr auto cols = index_dimension_of_v<A, 1>;
1441  using MatrixType = Eigen3::eigen_matrix_t<Scalar, rows, cols>;
1442  using ResultType = Eigen3::eigen_matrix_t<Scalar, cols, cols>;
1443 
1444  Eigen::HouseholderQR<MatrixType> QR {std::forward<A>(a)};
1445 
1446  if constexpr (dynamic_dimension<A, 1>)
1447  {
1448  auto rt_cols = get_index_dimension_of<1>(a);
1449 
1450  ResultType ret {rt_cols, rt_cols};
1451 
1452  if constexpr (dynamic_dimension<A, 0>)
1453  {
1454  auto rt_rows = get_index_dimension_of<0>(a);
1455 
1456  if (rt_rows < rt_cols)
1457  ret << QR.matrixQR().topRows(rt_rows),
1459  else
1460  ret = QR.matrixQR().topRows(rt_cols);
1461  }
1462  else
1463  {
1464  if (rows < rt_cols)
1465  ret << QR.matrixQR().template topRows<rows>(),
1467  else
1468  ret = QR.matrixQR().topRows(rt_cols);
1469  }
1470 
1471  return ret;
1472  }
1473  else
1474  {
1475  ResultType ret;
1476 
1477  if constexpr (dynamic_dimension<A, 0>)
1478  {
1479  auto rt_rows = get_index_dimension_of<0>(a);
1480 
1481  if (rt_rows < cols)
1482  ret << QR.matrixQR().topRows(rt_rows),
1484  else
1485  ret = QR.matrixQR().template topRows<cols>();
1486  }
1487  else
1488  {
1489  if constexpr (rows < cols)
1490  ret << QR.matrixQR().template topRows<rows>(), Eigen3::eigen_matrix_t<Scalar, cols - rows, cols>::Zero();
1491  else
1492  ret = QR.matrixQR().template topRows<cols>();
1493  }
1494 
1495  return ret;
1496  }
1497  }
1498 
1499  public:
1500 
1501  template<typename A>
1502  static constexpr auto
1503  LQ_decomposition(A&& a)
1504  {
1505  return OpenKalman::make_triangular_matrix<TriangleType::lower>(OpenKalman::adjoint(QR_decomp_impl(OpenKalman::adjoint(std::forward<A>(a)))));
1506  }
1507 
1508 
1509  template<typename A>
1510  static constexpr auto
1511  QR_decomposition(A&& a)
1512  {
1513  return OpenKalman::make_triangular_matrix<TriangleType::upper>(QR_decomp_impl(std::forward<A>(a)));
1514  }
1515 
1516  };
1517 
1518 
1519 } // namespace OpenKalman::interface
1520 
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&#39;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&#39;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