n_ary_operation.hpp

Go to the documentation of this file.

/* This file is part of OpenKalman, a header-only C++ library for
 * Kalman filters and other recursive filters.
 *
 * Copyright (c) 2022 Christopher Lee Ogden <ogden@gatech.edu>
 *
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at https://mozilla.org/MPL/2.0/.
 */

#ifndef OPENKALMAN_N_ARY_OPERATION_HPP
#define OPENKALMAN_N_ARY_OPERATION_HPP

namespace OpenKalman
{
  // ----------------- //
  //  n_ary_operation  //
  // ----------------- //

  namespace detail
  {

    // Check that dimensions of arguments Args are compatible with \ref coordinates::pattern Ds.
    template<std::size_t...ixs, typename DTup, typename...Args>
    constexpr void check_n_ary_dims(std::index_sequence<ixs...>, const DTup& d_tup, const Args&...args)
    {
      return ([](const DTup& d_tup, const Args&...args){
        constexpr auto ix = ixs;
        return ([](const DTup& d_tup, const auto& arg){
          using Arg = decltype(arg);
          if constexpr (dynamic_dimension<Arg, ix> or dynamic_pattern<std::tuple_element_t<ix, DTup>>)
          {
            auto arg_d = get_vector_space_descriptor<ix>(arg);
            auto tup_d = std::get<ix>(d_tup);
            auto dim_arg_d = get_dimension(arg_d);
            auto dim_tup_d = get_dimension(tup_d);
            if (not (arg_d == tup_d) and not internal::is_uniform_component_of(arg_d, tup_d))
              throw std::logic_error {"In an argument to n_ary_operation, the dimension of index " +
                std::to_string(ix) + " is " + std::to_string(dim_arg_d) + ", but should be 1 " +
                (dim_tup_d == 1 ? "" : "or " + std::to_string(dim_tup_d)) +
                "(the dimension of Ds number " + std::to_string(ix)};
          }
          else
          {
            using D_Arg = vector_space_descriptor_of_t<Arg, ix>;
            using D = std::tuple_element_t<ix, DTup>;
            static_assert(compares_with<D_Arg, D> or equivalent_to_uniform_static_vector_space_descriptor_component_of<D_Arg, D>or
              (ix >= index_count_v<Arg> and uniform_static_vector_space_descriptor<D>),
              "In argument to n_ary_operation, the dimension of each index must be either 1 or that of Ds.");
          }
        }(d_tup, args),...);
      }(d_tup, args...),...);
    }


    template<typename Op, typename...Args, std::size_t...I>
    constexpr bool is_invocable_with_indices(std::index_sequence<I...>)
    {
      return std::is_invocable_v<Op, Args..., decltype(I)...>;
    }


    template<typename Op, std::size_t...I, typename...Args>
    constexpr decltype(auto) n_ary_invoke_op(const Op& op, std::index_sequence<I...> seq, Args&&...args)
    {
      if constexpr (is_invocable_with_indices<const Op&, Args&&...>(seq))
        return op(std::forward<Args>(args)..., static_cast<decltype(I)>(0)...);
      else
        return op(std::forward<Args>(args)...);
    }


#ifdef __cpp_concepts
    template<typename Operation, std::size_t indices, typename...Args>
#else
    template<typename Operation, std::size_t indices, typename = void, typename...Args>
#endif
    struct n_ary_operator_traits_impl {};


#ifdef __cpp_concepts
    template<typename Op, std::size_t indices, typename...Args>
    requires (is_invocable_with_indices<Op, Args...>(std::make_index_sequence<indices> {})) or
      std::is_invocable_v<Op, Args...>
    struct n_ary_operator_traits_impl<Op, indices, Args...>
#else
    template<typename Op, std::size_t indices, typename...Args>
    struct n_ary_operator_traits_impl<Op, indices, std::enable_if_t<
      is_invocable_with_indices<Op, Args...>(std::make_index_sequence<indices> {}) or
      std::is_invocable_v<Op, Args...>>, Args...>
#endif
    {
      using type = decltype(n_ary_invoke_op(std::declval<Op>(), std::make_index_sequence<indices> {}, std::declval<Args>()...));
    };


    template<typename Op, std::size_t indices, typename...Args>
    struct n_ary_operator_traits
#ifdef __cpp_concepts
      : n_ary_operator_traits_impl<Op, indices, Args...> {};
#else
      : n_ary_operator_traits_impl<Op, indices, void, Args...> {};
#endif


#ifndef __cpp_concepts
    template<typename Op, std::size_t Indices, typename = void, typename...Args>
    struct n_ary_operator_impl : std::false_type {};

    template<typename Op, std::size_t Indices, typename...Args>
    struct n_ary_operator_impl<Op, Indices, std::enable_if_t<
      std::is_invocable<Op, typename std::add_lvalue_reference<typename scalar_type_of<Args>::type>::type...>::value or
      is_invocable_with_indices<Op, typename std::add_lvalue_reference<typename scalar_type_of<Args>::type>::type...>(
        std::make_index_sequence<Indices> {})>, Args...>
    : std::true_type {};
#endif


    template<typename Op, std::size_t Indices, typename...Args>
#ifdef __cpp_concepts
    concept n_ary_operator = std::is_invocable_v<Op, std::add_lvalue_reference_t<scalar_type_of_t<Args>>...> or
      is_invocable_with_indices<Op, std::add_lvalue_reference_t<scalar_type_of_t<Args>>...>(std::make_index_sequence<Indices> {});
#else
    constexpr bool n_ary_operator = n_ary_operator_impl<Op, Indices, void, Args...>::value;
#endif


    template<typename Arg, std::size_t...I, typename...J>
    inline auto n_ary_operation_get_component_impl(Arg&& arg, std::index_sequence<I...>, J...j)
    {
      if constexpr (sizeof...(I) == sizeof...(J))
        return get_component(std::forward<Arg>(arg), (j < get_index_dimension_of<I>(arg) ? j : 0)...);
      else
        return get_component(std::forward<Arg>(arg), [](auto dim, const auto& j_tup){
          auto j = std::get<I>(j_tup);
          if (j < dim) return j;
          else return 0_uz;
        }(get_index_dimension_of<I>(arg), std::tuple {j...})...);
    }


    template<typename Op, typename ArgsTup, std::size_t...ArgI, typename...J>
    inline auto n_ary_operation_get_component(const Op& op, ArgsTup&& args_tup, std::index_sequence<ArgI...>, J...j)
    {
      if constexpr (std::is_invocable_v<const Op&, scalar_type_of_t<std::tuple_element_t<ArgI, std::decay_t<ArgsTup>>>..., J...>)
        return op(n_ary_operation_get_component_impl(
          std::get<ArgI>(std::forward<ArgsTup>(args_tup)),
          std::make_index_sequence<index_count_v<std::tuple_element_t<ArgI, std::decay_t<ArgsTup>>>> {},
          j...)..., j...);
      else
        return op(n_ary_operation_get_component_impl(
          std::get<ArgI>(std::forward<ArgsTup>(args_tup)),
          std::make_index_sequence<index_count_v<std::tuple_element_t<ArgI, std::decay_t<ArgsTup>>>> {},
          j...)...);
    }


    template<typename M, typename Op, typename ArgsTup, typename...J>
    inline void n_ary_operation_iterate(M& m, const Op& op, ArgsTup&& args_tup, std::index_sequence<>, J...j)
    {
      std::make_index_sequence<std::tuple_size_v<ArgsTup>> seq;
      set_component(m, n_ary_operation_get_component(op, std::forward<ArgsTup>(args_tup), seq, j...), j...);
    }


    template<typename M, typename Op, typename ArgsTup, std::size_t I, std::size_t...Is, typename...J>
    inline void n_ary_operation_iterate(M& m, const Op& op, ArgsTup&& args_tup, std::index_sequence<I, Is...>, J...j)
    {
      for (std::size_t i = 0; i < get_index_dimension_of<I>(m); i++)
        n_ary_operation_iterate(m, op, std::forward<ArgsTup>(args_tup), std::index_sequence<Is...> {}, j..., i);
    }


    template<typename...Ds, typename Arg, std::size_t...Ix_Ds>
    static decltype(auto)
    replicate_arg(const std::tuple<Ds...>& d_tup, Arg&& arg, std::index_sequence<Ix_Ds...>)
    {
      OpenKalman::broadcast(std::forward<Arg>(arg),
        values::operation {
          std::divides<scalar_type_of_t<Arg>>{},
          get_dimension(std::get<Ix_Ds>(d_tup)),
          get_index_dimension_of<Ix_Ds>(arg)}...);
    }


    template<typename PatternMatrix, typename...Ds, typename Op, typename...Args>
    static constexpr auto
    n_ary_operation_impl(const std::tuple<Ds...>& d_tup, Op&& op, Args&&...args)
    {
      constexpr std::index_sequence_for<Ds...> seq;

      // constant_matrix:
      if constexpr (sizeof...(Args) > 0 and (constant_matrix<Args> and ...) and not is_invocable_with_indices<Op, scalar_type_of_t<Args>...>(seq))
      {
        values::operation c {op, constant_coefficient {std::forward<Args>(args)}...};
        return std::apply(
          [](auto&&...as){ return make_constant<PatternMatrix>(std::forward<decltype(as)>(as)...); },
          std::tuple_cat(std::tuple{std::move(c)}, d_tup));
      }
      // Library handles n-ary operation.
      else if constexpr (is_invocable_with_indices<Op, std::add_lvalue_reference_t<scalar_type_of_t<Args>>...>(seq) and
        interface::n_ary_operation_defined_for<PatternMatrix, const std::tuple<Ds...>&, Op&&, Args&&...>)
      {
        using Interface = interface::library_interface<std::decay_t<PatternMatrix>>;
        auto ret {Interface::n_ary_operation(d_tup, std::forward<Op>(op), replicate_arg(d_tup, std::forward<Args>(args), seq)...)};
        return std::apply([](auto&& a, auto&&...vs){
          return make_vector_space_adapter(std::forward<decltype(a)>(a), std::forward<decltype(vs)>(vs)...);
          }, std::tuple_cat(std::forward_as_tuple(std::move(ret)), d_tup));
      }
      else // Catch-all: library does not provide for this n-ary operation.
      {
        using Scalar = std::decay_t<typename n_ary_operator_traits<Op, sizeof...(Ds),
          std::add_lvalue_reference_t<scalar_type_of_t<Args>>...>::type>;

        if constexpr (((coordinates::dimension_of_v<Ds> == 1) and ...))
        {
          // one-by-one matrix
          auto e = op(get_component(std::forward<Args>(args))...);
          return make_dense_object_from<PatternMatrix, Layout::none, Scalar>(d_tup, e);
        }
        else
        {
          auto m = std::apply([](auto&&...ds){
            return make_dense_object<PatternMatrix, Layout::none, Scalar>(std::forward<decltype(ds)>(ds)...);
          }, d_tup);
          n_ary_operation_iterate(m, op, std::forward_as_tuple(std::forward<Args>(args)...), seq);
          return m;
        }
      }
    }

  } // namespace detail


#ifdef __cpp_concepts
  template<coordinates::pattern...Ds, typename Operation, indexible...Args> requires (sizeof...(Args) > 0) and
    detail::n_ary_operator<Operation, sizeof...(Ds), Args...> and (... and (coordinates::dimension_of_v<Ds> != 0))
  constexpr compatible_with_vector_space_descriptor_collection<std::tuple<Ds...>> auto
#else
  template<typename...Ds, typename Operation, typename...Args, std::enable_if_t<
    (coordinates::pattern<Ds> and ...) and (indexible<Args> and ...) and (sizeof...(Args) > 0) and
    detail::n_ary_operator<Operation, sizeof...(Ds), Args...> and (... and (coordinates::dimension_of_v<Ds> != 0)), int> = 0>
  constexpr auto
#endif
  n_ary_operation(const std::tuple<Ds...>& d_tup, Operation&& operation, Args&&...args)
  {
    detail::check_n_ary_dims(std::index_sequence_for<Ds...> {}, d_tup, args...);
    using Arg0 = std::decay_t<std::tuple_element_t<0, std::tuple<Args...>>>; // \todo Pick the first appropriate pattern matrix, even if not the first one.
    return detail::n_ary_operation_impl<Arg0>(d_tup, std::forward<Operation>(operation), std::forward<Args>(args)...);
  }


  namespace detail
  {
    template<std::size_t ix, typename Arg, typename...Args>
    constexpr auto find_max_dim(const Arg& arg, const Args&...args)
    {
      if constexpr (sizeof...(Args) == 0)
      {
        auto ret = get_vector_space_descriptor<ix>(arg);
        if constexpr (dynamic_pattern<decltype(ret)>)
        {
          if (get_dimension(ret) == 0) throw std::invalid_argument {"A dimension of an arguments "
            "to n_ary_operation is zero for at least index " + std::to_string(ix) + "."};
        }
        else static_assert(index_dimension_of_v<Arg, ix> != 0, "Arguments to n_ary_operation cannot have zero dimensions");
        return ret;
      }
      else
      {
        auto max_d = find_max_dim<ix>(args...);
        using Arg_D = vector_space_descriptor_of_t<Arg, ix>;
        using Max_D = decltype(max_d);

        if constexpr (fixed_pattern<Arg_D> and fixed_pattern<Max_D>)
        {
          constexpr auto dim_arg_d = coordinates::dimension_of_v<Arg_D>;
          if constexpr (compares_with<Arg_D, Max_D>or (dim_arg_d == 1 and equivalent_to_uniform_static_vector_space_descriptor_component_of<Arg_D, Max_D>))
          {
            return max_d;
          }
          else
          {
            constexpr auto dim_max_d = coordinates::dimension_of_v<Max_D>;
            static_assert(dim_max_d != 1 or not equivalent_to_uniform_static_vector_space_descriptor_component_of<Max_D, Arg_D>,
              "The dimension of arguments to n_ary_operation are not compatible with each other for at least one index.");
            return get_vector_space_descriptor<ix>(arg);
          }
        }
        else if constexpr (coordinates::euclidean_pattern<Arg_D> and coordinates::euclidean_pattern<Max_D>)
        {
          if constexpr (fixed_pattern<Arg_D>)
          {
            constexpr std::size_t a = coordinates::dimension_of_v<Arg_D>;
            std::size_t m = get_dimension(max_d);
            if (a != m and a != 1 and m != 1) throw std::invalid_argument {"The dimension of arguments to n_ary_operation "
                "are not compatible with each other for at least index " + std::to_string(ix) + "."};

            if constexpr (a == 1) return max_d;
            else return get_vector_space_descriptor<ix>(arg);
          }
          else if constexpr (fixed_pattern<Max_D>)
          {
            auto arg_d = get_vector_space_descriptor<ix>(arg);
            std::size_t a = get_dimension(arg_d);
            constexpr std::size_t m = coordinates::dimension_of_v<Max_D>;
            if (a != m and a != 1 and m != 1) throw std::invalid_argument {"The dimension of arguments to n_ary_operation "
                "are not compatible with each other for at least index " + std::to_string(ix) + "."};

            if constexpr (m == 1) return arg_d;
            else return max_d;
          }
          else
          {
            std::size_t a = get_index_dimension_of<ix>(arg);
            std::size_t m = get_dimension(max_d);
            if (a == m or a == 1) return m;
            else if (m == 1 and m <= a) return a;
            else throw std::invalid_argument {"The dimension of arguments to n_ary_operation are not compatible with "
              "each other for at least index " + std::to_string(ix) + "."};
          }
        }
        else
        {
          auto arg_d = get_vector_space_descriptor<ix>(arg);
          //using Scalar = scalar_type_of_t<Arg>;
          if (internal::is_uniform_component_of(arg_d, max_d))
          {
            //if constexpr (internal::is_DynamicDescriptor<Max_D>::value) return DynamicDescriptor {max_d};
            //else return DynamicDescriptor<Scalar> {max_d};
            return max_d;
          }
          else if (internal::is_uniform_component_of(max_d, arg_d))
          {
            //if constexpr (internal::is_DynamicDescriptor<Max_D>::value) return DynamicDescriptor {arg_d};
            //else return DynamicDescriptor<Scalar> {arg_d};
            return arg_d;
          }
          else throw std::invalid_argument {"The dimension of arguments to n_ary_operation are not compatible with "
            "each other for at least index " + std::to_string(ix) + "."};
        }
      }
    }


    template<std::size_t...ixs, typename...Args>
    constexpr auto find_max_dims(std::index_sequence<ixs...>, const Args&...args)
    {
      return std::tuple {find_max_dim<ixs>(args...)...};
    }

  } // namespace detail


#ifdef __cpp_concepts
  template<typename Operation, indexible...Args> requires (sizeof...(Args) > 0) and
    detail::n_ary_operator<Operation, std::max({index_count_v<Args>...}), Args...>
#else
  template<typename Operation, typename...Args, std::enable_if_t<(indexible<Args> and ...) and
    (sizeof...(Args) > 0) and detail::n_ary_operator<Operation, std::max({index_count<Args>::value...}), Args...>, int> = 0>
#endif
  constexpr auto
  n_ary_operation(Operation&& operation, Args&&...args)
  {
    auto d_tup = detail::find_max_dims(std::make_index_sequence<std::max({index_count_v<Args>...})> {}, args...);
    using Arg0 = std::decay_t<std::tuple_element_t<0, std::tuple<Args...>>>;
    return detail::n_ary_operation_impl<Arg0>(std::move(d_tup), std::forward<Operation>(operation), std::forward<Args>(args)...);
  }


  namespace detail
  {
    template<typename M, typename Operation, typename Vs_tuple, typename Index_seq, typename K_seq, typename...Is>
    void nullary_set_components(M& m, const Operation& op, const Vs_tuple&, Index_seq, K_seq, Is...is)
    {
      constexpr auto seq = std::make_index_sequence<sizeof...(Is)> {};
      if constexpr (detail::is_invocable_with_indices<Operation>(seq))
        set_component(m, op(is...), is...); //< Operation takes a full set of indices.
      else
        set_component(m, op(), is...); //< Operation takes no indices.
    }


    template<std::size_t DsIndex, std::size_t...DsIndices, typename M, typename Operation,
      typename Vs_tuple, std::size_t...indices, std::size_t...Ks, typename...Is>
    void nullary_set_components(M& m, const Operation& op, const Vs_tuple& ds_tup,
      std::index_sequence<indices...> index_seq, std::index_sequence<Ks...> k_seq, Is...is)
    {
      if constexpr (((DsIndex == indices) or ...))
      {
        constexpr std::integral_constant<size_t, ((DsIndex == indices ? Ks : 0) + ...)> i;
        nullary_set_components<DsIndices...>(m, op, ds_tup, index_seq, k_seq, is..., i);
      }
      else
      {
        // Iterate through the dimensions of the current DsIndex and add set elements for each dimension iteratively.
        for (std::size_t i = 0; i < get_dimension(std::get<DsIndex>(ds_tup)); ++i)
        {
          nullary_set_components<DsIndices...>(m, op, ds_tup, index_seq, k_seq, is..., i);
        }
      }
    }


    template<std::size_t CurrentOpIndex, std::size_t factor, typename M, typename Operations_tuple,
      typename Vs_tuple, typename UniqueIndicesSeq, std::size_t...AllDsIndices, typename K_seq>
    void nullary_iterate(M& m, const Operations_tuple& op_tup, const Vs_tuple& ds_tup,
      UniqueIndicesSeq unique_indices_seq, std::index_sequence<AllDsIndices...>, K_seq k_seq)
    {
      nullary_set_components<AllDsIndices...>(m, std::get<CurrentOpIndex>(op_tup), ds_tup, unique_indices_seq, k_seq);
    }


    template<std::size_t CurrentOpIndex, std::size_t factor, std::size_t index, std::size_t...indices,
      typename M, typename Operations_tuple, typename Vs_tuple, typename UniqueIndicesSeq, typename AllDsSeq,
      std::size_t...Ks, std::size_t...Js, typename...J_seqs>
    void nullary_iterate(M& m, const Operations_tuple& op_tup, const Vs_tuple& ds_tup,
      UniqueIndicesSeq unique_indices_seq, AllDsSeq all_ds_seq, std::index_sequence<Ks...>, std::index_sequence<Js...>,
      J_seqs...j_seqs)
    {
      constexpr std::size_t new_factor = factor / coordinates::dimension_of_v<std::tuple_element_t<index, Vs_tuple>>;

      (nullary_iterate<CurrentOpIndex + new_factor * Js, new_factor, indices...>(
        m, op_tup, ds_tup, unique_indices_seq, all_ds_seq, std::index_sequence<Ks..., Js>{}, j_seqs...),...);
    }

  } // namespace detail


  #ifdef __cpp_concepts
  template<indexible PatternMatrix, std::size_t...indices, coordinates::pattern...Ds, typename...Operations>
  requires ((fixed_pattern<std::tuple_element_t<indices, std::tuple<Ds...>>>) and ...) and
    (sizeof...(Operations) == (1 * ... * coordinates::dimension_of_v<std::tuple_element_t<indices, std::tuple<Ds...>>>)) and
    (detail::n_ary_operator<Operations, sizeof...(Ds)> and ...)
  #else
  template<typename PatternMatrix, std::size_t...indices, typename...Ds, typename...Operations, std::enable_if_t<
    indexible<PatternMatrix> and (coordinates::pattern<Ds> and ...) and
    ((fixed_pattern<std::tuple_element_t<indices, std::tuple<Ds...>>>) and ...) and
    (sizeof...(Operations) == (1 * ... * coordinates::dimension_of<std::tuple_element_t<indices, std::tuple<Ds...>>>::value)) and
    (detail::n_ary_operator<Operations, sizeof...(Ds)> and ...), int> = 0>
  #endif
  constexpr auto
  n_ary_operation(const std::tuple<Ds...>& d_tup, const Operations&...operations)
  {
    using Scalar = std::common_type_t<std::decay_t<typename detail::n_ary_operator_traits<Operations, sizeof...(Ds)>::type>...>;

    // One operation for all elements combined:
    if constexpr (sizeof...(Operations) == 1)
    {
      return detail::n_ary_operation_impl<std::decay_t<PatternMatrix>>(d_tup, operations...);
    }
    // One operation for each element, and the operations are not invocable with indices:
    else if constexpr (((not dynamic_pattern<Ds>) and ...) and
      sizeof...(operations) == (coordinates::dimension_of_v<Ds> * ...) and
      not (detail::is_invocable_with_indices<const Operations&>(std::make_index_sequence<sizeof...(Ds)> {}) or ...))
    {
      return make_dense_object_from<PatternMatrix, Layout::none, Scalar>(d_tup, operations()...);
    }
    // All other cases:
    else
    {
      auto m = make_dense_object<PatternMatrix, Layout::none, Scalar>(d_tup);
      detail::nullary_iterate<0, sizeof...(Operations), indices...>(
        m,
        std::forward_as_tuple(operations...),
        d_tup,
        std::index_sequence<indices...> {},
        std::index_sequence_for<Ds...> {},
        std::index_sequence<> {},
        std::make_index_sequence<coordinates::dimension_of_v<std::tuple_element_t<indices, std::tuple<Ds...>>>> {}...);
      return m;
    }
  }


#ifdef __cpp_concepts
  template<indexible PatternMatrix, std::size_t...indices, coordinates::pattern...Ds, typename...Operations>
  requires ((fixed_pattern<vector_space_descriptor_of_t<PatternMatrix, indices>>) and ...) and
    (sizeof...(Operations) == (1 * ... * index_dimension_of_v<PatternMatrix, indices>))
#else
  template<typename PatternMatrix, std::size_t...indices, typename...Ds, typename...Operations, std::enable_if_t<
    indexible<PatternMatrix> and (coordinates::pattern<Ds> and ...) and
    ((fixed_pattern<typename vector_space_descriptor_of<PatternMatrix, indices>::type>) and ...) and
    (sizeof...(Operations) == (1 * ... * index_dimension_of<PatternMatrix, indices>::value)), int> = 0>
#endif
  constexpr auto
  n_ary_operation(const Operations&...operations)
  {
    auto d_tup = all_vector_space_descriptors<PatternMatrix>();
    return n_ary_operation<PatternMatrix, indices...>(d_tup, operations...);
  }


  // -------------------------- //
  //  n_ary_operation_in_place  //
  // -------------------------- //

  namespace detail
  {
    template<typename Operation, typename Elem, typename...J>
    inline void do_elem_operation_in_place_impl(const Operation& operation, Elem& elem, J...j)
    {
      if constexpr (std::is_invocable_v<const Operation&, Elem&, J...>)
        operation(elem, j...);
      else
        operation(elem);
    }


    template<typename Operation, typename Arg, typename...J>
    inline void do_elem_operation_in_place(const Operation& operation, Arg& arg, J...j)
    {
      auto&& elem = get_component(arg, j...);
      if constexpr (std::is_assignable_v<decltype((elem)), std::decay_t<decltype(elem)>>)
      {
        do_elem_operation_in_place_impl(operation, elem, j...);
      }
      else
      {
        auto e {std::forward<decltype(elem)>(elem)}; // copy elem
        static_assert(std::is_assignable_v<decltype((e)), std::decay_t<decltype(elem)>>);
        do_elem_operation_in_place_impl(operation, e, j...);
        set_component(arg, std::move(e), j...);
      }
    }


    template<typename Operation, typename Arg, typename Count, typename...J>
    inline void unary_operation_in_place_impl(const Operation& operation, Arg& arg, const Count& count, J...j)
    {
      constexpr auto n = sizeof...(J);
      if constexpr (n < Count::value)
      {
        for (std::size_t i = 0; i < get_index_dimension_of<n>(arg); ++i)
          unary_operation_in_place_impl(operation, arg, count, j..., i);
      }
      else
      {
        do_elem_operation_in_place(operation, arg, j...);
      }
    }

  } // namespace detail


  #ifdef __cpp_concepts
  template<typename Operation, writable Arg> requires detail::n_ary_operator<Operation, index_count_v<Arg>, Arg>
  #else
  template<typename Operation, typename Arg, std::enable_if_t<writable<Arg> and
    detail::n_ary_operator<Operation, index_count_v<Arg>, Arg>, int> = 0>
  #endif
  constexpr decltype(auto)
  unary_operation_in_place(const Operation& operation, Arg&& arg)
  {
    // \todo If the native library has its own facilities for doing this, use them.

    detail::unary_operation_in_place_impl(operation, arg, count_indices(arg));
    return std::forward<Arg>(arg);
  }


} // namespace OpenKalman

#endif //OPENKALMAN_N_ARY_OPERATION_HPP