nonlinearoperator.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
#include <functional>
 
#include <dune/common/hybridutilities.hh>
 
#include <ikarus/utils/traits.hh>
 
namespace Ikarus {
 
namespace Impl {
  template <class F, class Tuple, std::size_t... I>
  constexpr decltype(auto) applyAndRemoveRefererenceWrapper(F&& f, Tuple&& t, std::index_sequence<I...>) {
    return std::invoke(std::forward<F>(f),
                       std::get<I>(std::forward<Tuple>(t)).get()...); //.get gets the impl type of std::referenceWrapper
  }
 
  template <class F, class Tuple>
  constexpr decltype(auto) applyAndRemoveReferenceWrapper(F&& f, Tuple&& t) {
    return applyAndRemoveRefererenceWrapper(
        std::forward<F>(f), std::forward<Tuple>(t),
        std::make_index_sequence<std::tuple_size_v<std::remove_reference_t<Tuple>>>{});
  }
 
  template <typename T>
  auto forwardasReferenceWrapperIfIsReference(T&& t) {
    if constexpr (std::is_lvalue_reference_v<decltype(t)>)
      return std::ref(t);
    else
      return t;
  }
 
  template <class Pars, class Tuple, std::size_t... I>
  constexpr decltype(auto) makeTupleOfValuesAndReferences(Tuple&& t, Pars&& p, std::index_sequence<I...>) {
    return std::make_tuple(
        forwardasReferenceWrapperIfIsReference(applyAndRemoveReferenceWrapper(std::get<I>(t), p))...);
  }
 
  template <typename... Args>
  struct Functions
  {
    std::tuple<Args...> args;
  };
 
  template <typename... Args>
  struct Parameter
  {
    std::tuple<Args...> args;
  };
 
} // namespace Impl
 
template <typename... Args>
auto parameter(Args&&... args) {
  return Impl::Parameter<Args&&...>{std::forward_as_tuple(std::forward<Args>(args)...)};
}
 
template <typename... Args>
auto functions(Args&&... args) {
  return Impl::Functions<Args&&...>{std::forward_as_tuple(std::forward<Args>(args)...)};
}
 
template <typename... DerivativeArgs, typename... ParameterArgs>
auto initResults(const std::tuple<DerivativeArgs...>& derivativesFunctions,
                 const std::tuple<ParameterArgs...>& parameter) {
  return Impl::makeTupleOfValuesAndReferences(
      derivativesFunctions, parameter,
      std::make_index_sequence<std::tuple_size_v<std::remove_reference_t<decltype(derivativesFunctions)>>>{});
}
 
template <typename TypeListOne, typename TypeListTwo>
class NonLinearOperator
{
public:
  NonLinearOperator([[maybe_unused]] const TypeListOne& derivativesFunctions,
                    [[maybe_unused]] const TypeListTwo& args) {
    static_assert(!sizeof(TypeListOne),
                  "This type should not be instantiated. check that your arguments satisfies the template below");
  }
};
 
template <typename... DerivativeArgs, typename... ParameterArgs>
class NonLinearOperator<Impl::Functions<DerivativeArgs...>, Impl::Parameter<ParameterArgs...>>
{
public:
  using FunctionReturnValues =
      std::tuple<Ikarus::traits::ReturnType<DerivativeArgs, ParameterArgs&...>...>; 
  using ParameterValues = std::tuple<ParameterArgs...>;                             
 
  static constexpr int numberOfFunctions  = sizeof...(DerivativeArgs); 
  static constexpr int numberOfParameters = sizeof...(ParameterArgs);  
 
  template <int n>
  using FunctionReturnType = std::tuple_element_t<n, FunctionReturnValues>;
 
  template <int n>
  using ParameterValue = std::remove_cvref_t<std::tuple_element_t<n, ParameterValues>>;
 
  using ValueType =
      std::remove_cvref_t<std::tuple_element_t<0, FunctionReturnValues>>; 
  using DerivativeType =
      std::remove_cvref_t<std::tuple_element_t<1, FunctionReturnValues>>; 
 
  template <typename U = void>
  requires(not std::is_rvalue_reference_v<DerivativeArgs> and ...)
  explicit NonLinearOperator(const Impl::Functions<DerivativeArgs...>& derivativesFunctions,
                             const Impl::Parameter<ParameterArgs...>& parameterI)
      : derivatives_{derivativesFunctions.args},
        args_{parameterI.args},
        derivativesEvaluated_(initResults(derivatives_, args_)) {}
 
  template <typename Funcs>
  explicit NonLinearOperator(const Funcs& derivativesFunctions, const Impl::Parameter<ParameterArgs...>& parameterI)
      : derivatives_{derivativesFunctions.args},
        args_{parameterI.args},
        derivativesEvaluated_(initResults(derivatives_, args_)) {}
 
  void updateAll() {
    Dune::Hybrid::forEach(
        Dune::Hybrid::integralRange(Dune::index_constant<sizeof...(DerivativeArgs)>()), [&](const auto i) {
          std::get<i>(derivativesEvaluated_) = Impl::applyAndRemoveReferenceWrapper(std::get<i>(derivatives_), args_);
        });
  }
 
  template <int n>
  void update() {
    std::get<n>(derivativesEvaluated_) = Impl::applyAndRemoveReferenceWrapper(std::get<n>(derivatives_), args_);
  }
 
  auto& value()
  requires(sizeof...(DerivativeArgs) > 0)
  {
    return nthDerivative<0>();
  }
 
  auto& derivative()
  requires(sizeof...(DerivativeArgs) > 1)
  {
    return nthDerivative<1>();
  }
 
  auto& secondDerivative()
  requires(sizeof...(DerivativeArgs) > 2)
  {
    return nthDerivative<2>();
  }
 
  template <int n>
  auto& nthDerivative()
  requires(sizeof...(DerivativeArgs) > n)
  {
    if constexpr (requires { std::get<n>(derivativesEvaluated_).get(); })
      return std::get<n>(derivativesEvaluated_).get();
    else
      return std::get<n>(derivativesEvaluated_);
  }
 
  auto& lastParameter() { return nthParameter<sizeof...(ParameterArgs) - 1>(); }
  auto& firstParameter()
  requires(sizeof...(ParameterArgs) > 0)
  {
    return nthParameter<0>();
  }
  auto& secondParameter()
  requires(sizeof...(ParameterArgs) > 1)
  {
    return nthParameter<1>();
  }
  template <int n>
  auto& nthParameter()
  requires(sizeof...(ParameterArgs) >= n)
  {
    return std::get<n>(args_).get();
  }
 
  template <int... Derivatives>
  auto subOperator() {
    auto derivatives = derivatives_;
    auto fs = functions([&derivatives]() -> decltype(auto) { return std::get<Derivatives>(derivatives); }()...);
    Ikarus::NonLinearOperator<Impl::Functions<std::tuple_element_t<Derivatives, decltype(derivatives_)>...>,
                              Impl::Parameter<ParameterArgs...>>
        subOp(std::move(fs), Impl::applyAndRemoveReferenceWrapper(parameter<ParameterArgs...>, args_));
 
    return subOp;
  }
 
private:
  using FunctionReturnValuesWrapper = std::tuple<std::conditional_t<
      std::is_reference_v<Ikarus::traits::ReturnType<DerivativeArgs, ParameterArgs&...>>,
      std::reference_wrapper<std::remove_reference_t<Ikarus::traits::ReturnType<DerivativeArgs, ParameterArgs&...>>>,
      std::remove_cvref_t<Ikarus::traits::ReturnType<DerivativeArgs, ParameterArgs&...>>>...>;
 
  std::tuple<std::conditional_t<std::is_lvalue_reference_v<DerivativeArgs>,
                                std::reference_wrapper<std::remove_reference_t<DerivativeArgs>>,
                                std::remove_reference_t<DerivativeArgs>>...>
      derivatives_;
 
  std::tuple<std::conditional_t<std::is_lvalue_reference_v<ParameterArgs>,
                                std::reference_wrapper<std::remove_reference_t<ParameterArgs>>,
                                std::remove_reference_t<ParameterArgs>>...>
      args_;
  FunctionReturnValuesWrapper derivativesEvaluated_{};
};
 
template <typename... DerivativeArgs, typename... ParameterArgs>
NonLinearOperator(const Impl::Functions<DerivativeArgs&&...>& a, const Impl::Parameter<ParameterArgs...>& b)
    -> NonLinearOperator<Impl::Functions<DerivativeArgs...>, Impl::Parameter<ParameterArgs...>>;
} // namespace Ikarus