algorithm.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_SOLVERS_COMMON_ALGORITHM_HH
#define DUNE_SOLVERS_COMMON_ALGORITHM_HH

#include <dune/common/indices.hh>

#include <dune/common/typeutilities.hh>
#include <dune/common/typetraits.hh>


#include <dune/istl/multitypeblockvector.hh>

namespace Dune {
namespace Solvers {
namespace Hybrid {

namespace Imp {

  // Try if tuple_size is implemented for class
  template<class T, int i>
  constexpr auto size(const Dune::FieldVector<T, i>*, const PriorityTag<5>&)
    -> decltype(std::integral_constant<std::size_t,i>())
  {
    return {};
  }

  // Try if we have an instance of std::integer_sequence
  template<class T, T... t, class Index>
  constexpr auto size(std::integer_sequence<T, t...>, PriorityTag<4>)
  {
    using sizeAsType = std::tuple_size<decltype(std::make_tuple(t...))>;
    return std::integral_constant<std::size_t, sizeAsType::value>();
  }

  // Try if tuple_size is implemented for class
  template<class T>
  constexpr auto size(const T*, const PriorityTag<3>&)
    -> decltype(std::integral_constant<std::size_t,std::tuple_size<T>::value>())
  {
    return {};
  }

  // Try if there's a static constexpr size()
  template<class T>
  constexpr auto size(const T*, const PriorityTag<1>&)
    -> decltype(std::integral_constant<std::size_t,T::size()>())
  {
    return {};
  }

  // As a last resort try if there's a static constexpr size()
  template<class T>
  constexpr auto size(const T* t, const PriorityTag<0>&)
  {
    return t->size();
  }

} // namespace Imp



/**
 * \brief Size query
 *
 * \tparam T Type of container whose size is queried
 *
 * \param t Container whose size is queried
 *
 * \return Size of t
 *
 * If the size of t is known at compile type the size is
 * returned as std::integral_constant<std::size_t, size>.
 * Otherwise the result of t.size() is returned.
 *
 * Supported types for deriving the size at compile time are:
 * * instances of std::integer_sequence
 * * all types std::tuple_size is implemented for
 * * all typed that have a static method ::size()
 * * instances of Dune::FieldVector
 */
template<class T>
constexpr auto size(const T& t)
{
  return Imp::size(&t, PriorityTag<42>());
}



namespace Imp {

  template<class Container, class Index,
    std::enable_if_t<IsTuple<std::decay_t<Container>>::value, int> = 0>
  constexpr decltype(auto) elementAt(Container&& c, Index&&, PriorityTag<2>)
  {
    return std::get<Index::value>(c);
  }

  template<class T, T... t, class Index>
  constexpr decltype(auto) elementAt(std::integer_sequence<T, t...> c, Index&&, PriorityTag<1>)
  {
    return std::get<Index::value>(std::make_tuple(std::integral_constant<T, t>()...));
  }

  template<class Container, class Index>
  constexpr decltype(auto) elementAt(Container&& c, Index&& i, PriorityTag<0>)
  {
    return c[i];
  }

} // namespace Imp



/**
 * \brief Get element at given position from container
 *
 * \tparam Container Type of given container
 * \tparam Index Type of index
 *
 * \param c Given container
 * \param i Index of element to obtain
 *
 * \return The element at position i, i.e. c[i]
 *
 * If this returns the i-th entry of c. It supports the following
 * containers
 * * Containers providing dynamic access via operator[]
 * * Heterogenous containers providing access via operator[](integral_constant<...>)
 * * std::tuple<...>
 * * std::integer_sequence
 */
template<class Container, class Index>
constexpr decltype(auto) elementAt(Container&& c, Index&& i)
{
  return Imp::elementAt(std::forward<Container>(c), std::forward<Index>(i), PriorityTag<42>());
}



namespace Imp {

  template<class Begin, class End>
  class StaticIntegralRange
  {
  public:

    template<std::size_t i>
    constexpr auto operator[](Dune::index_constant<i>) const
    {
      return Dune::index_constant<Begin::value+i>();
    }

    static constexpr auto size()
    {
      return std::integral_constant<typename Begin::value_type, End::value - Begin::value>();
    }
  };

  template<class T>
  class DynamicIntegralRange
  {
  public:
    constexpr DynamicIntegralRange(const T& begin, const T& end):
      begin_(begin),
      end_(end)
    {}

    const T& begin() const
    { return begin_; }

    const T& end() const
    { return end_; }

    constexpr auto size() const
    {
      return end() - begin();
    }

    constexpr T operator[](const T&i) const
    { return begin()+i; }

  private:
    T begin_;
    T end_;
  };

  template<class Begin, class End,
    std::enable_if_t<IsIntegralConstant<Begin>::value and IsIntegralConstant<End>::value, int> = 0>
  constexpr auto integralRange(const Begin& begin, const End& end, const PriorityTag<1>&)
  {
    static_assert(Begin::value <= End::value, "You cannot create an integralRange where end<begin");
    return Imp::StaticIntegralRange<Begin,End>();
  }

  template<class Begin, class End>
  constexpr auto integralRange(const Begin& begin, const End& end, const PriorityTag<0>&)
  {
    assert(begin <= end);
    return Imp::DynamicIntegralRange<End>(begin, end);
  }

} // namespace Imp



/**
 * \brief Create an integral range
 *
 * \tparam Begin Type of begin entry of the range
 * \tparam End Type of end entry of the range
 *
 * \param begin First entry of the range
 * \param end One past the last entry of the range
 *
 * \returns An object encoding the given range
 *
 * If Begin and End are both instances of type
 * std::integral_constant, the returnes range
 * encodes begin and end statically.
 */
template<class Begin, class End>
constexpr auto integralRange(const Begin& begin, const End& end)
{
  return Imp::integralRange(begin, end, PriorityTag<42>());
}

/**
 * \brief Create an integral range starting from 0
 *
 * \tparam End Type of end entry of the range
 *
 * \param end One past the last entry of the range
 *
 * \returns An object encoding the given range
 *
 * This is a short cut for integralRange(_0, end).
 */
template<class End>
constexpr auto integralRange(const End& end)
{
  return Imp::integralRange(Dune::Indices::_0, end, PriorityTag<42>());
}



namespace Imp {

  template<class Range, class F, class Index, Index... i>
  constexpr void forEachIndex(Range&& range, F&& f, std::integer_sequence<Index, i...>)
  {
    std::initializer_list<int>{(f(Hybrid::elementAt(range, std::integral_constant<Index,i>())), 0)...};
  }

  template<class Range, class F,
    std::enable_if_t<IsIntegralConstant<decltype(Hybrid::size(std::declval<Range>()))>::value, int> = 0>
  constexpr void forEach(Range&& range, F&& f, PriorityTag<1>)
  {
    auto size = Hybrid::size(range);
    auto indices = std::make_index_sequence<size>();
    forEachIndex(std::forward<Range>(range), std::forward<F>(f), indices);
  }

  template<class Range, class F>
  constexpr void forEach(Range&& range, F&& f, PriorityTag<0>)
  {
    for(std::size_t i=0; i<range.size(); ++i)
      f(range[i]);
  // \ToDo Why does the following not compile?
  //  for(auto e : range)
  //    f(e);
  }

} // namespace Imp



/**
 * \brief Range based for loop
 *
 * \tparam Range Type of given range
 * \tparam F Type of given predicate
 *
 * \param range The range to loop over
 * \param f A predicate that will be called with each entry of the range
 *
 * This supports looping over the following ranges
 * * ranges obtained from integralRange()
 * * all ranges that provide Hybrid::size() and Hybrid::elementAt()
 *
 * This especially included instances of std::integer_sequence,
 * std::tuple, Dune::TupleVector, and Dune::MultiTypeBlockVector.
 */
template<class Range, class F>
constexpr void forEach(Range&& range, F&& f)
{
  Imp::forEach(std::forward<Range>(range), std::forward<F>(f), PriorityTag<42>());
}



namespace Imp {

  template<class IfFunc, class ElseFunc>
  constexpr void ifElse(std::true_type, IfFunc&& ifFunc, ElseFunc&& elseFunc)
  {
    ifFunc([](auto&& x) { return std::forward<decltype(x)>(x);});
  }

  template<class IfFunc, class ElseFunc>
  constexpr void ifElse(std::false_type, IfFunc&& ifFunc, ElseFunc&& elseFunc)
  {
    elseFunc([](auto&& x) { return std::forward<decltype(x)>(x);});
  }

  template<class IfFunc, class ElseFunc>
  constexpr void ifElse(const bool& condition, IfFunc&& ifFunc, ElseFunc&& elseFunc)
  {
    if (condition)
      ifFunc([](auto&& x) { return std::forward<decltype(x)>(x);});
    else
      elseFunc([](auto&& x) { return std::forward<decltype(x)>(x);});
  }

} // namespace Imp



/**
 * \brief A conditional expression
 *
 * This will call either ifFunc or elseFunc depending
 * on the condition. In any case a single argument
 * will be passed to the called function. This will always
 * be the indentity function. Passing an expression through
 * this function will lead to lazy evaluation. This way both
 * 'branches' can contain expressions that are only valid
 * within this branch if the condition is a std::integral_constant<bool,*>.
 *
 * In order to do this, the passed functors must have a single
 * argument of type auto.
 *
 * Due to the lazy evaluation mechanism and support for
 * std::integral_constant<bool,*> this allows to emulate
 * a static if statement.
 */
template<class Condition, class IfFunc, class ElseFunc>
constexpr void ifElse(const Condition& condition, IfFunc&& ifFunc, ElseFunc&& elseFunc)
{
  Imp::ifElse(condition, std::forward<IfFunc>(ifFunc), std::forward<ElseFunc>(elseFunc));
}

/**
 * \brief A conditional expression
 *
 * This provides an ifElse conditional with empty else clause.
 */
template<class Condition, class IfFunc>
constexpr void ifElse(const Condition& condition, IfFunc&& ifFunc)
{
  ifElse(condition, std::forward<IfFunc>(ifFunc), [](auto&& i) {});
}



namespace Imp {

  template<class T1, class T2>
  constexpr auto equals(const T1& t1, const T2& t2, PriorityTag<1>) -> decltype(T1::value, T2::value, std::integral_constant<bool,T1::value == T2::value>())
  { return {}; }

  template<class T1, class T2>
  constexpr auto equals(const T1& t1, const T2& t2, PriorityTag<0>)
  {
    return t1==t2;
  }

} // namespace Imp



/**
 * \brief Equality comparison
 *
 * If both types have a static member value, the result of comparing
 * these is returned as std::integral_constant<bool, *>. Otherwise
 * the result of a runtime comparison of t1 and t2 is directly returned.
 */
template<class T1, class T2>
constexpr auto equals(T1&& t1,  T2&& t2)
{
  return Imp::equals(std::forward<T1>(t1), std::forward<T2>(t2), PriorityTag<1>());
}



} // namespace Hybrid







// Implementation of integralRangeFor
namespace Imp {

template<class ST, ST begin, ST end>
struct StaticForLoop
{
  template<class F, class...Args>
  static void apply(F&& f, Args&&... args)
  {
    f(std::integral_constant<ST, begin>(), std::forward<Args>(args)...);
    StaticForLoop<ST, begin+1, end>::apply(std::forward<F>(f), std::forward<Args>(args)...);
  }
};

template<class ST, ST end>
struct StaticForLoop<ST, end, end>
{
  template<class F, class...Args>
  static void apply(F&& f, Args&&...)
  {}
};

// Overload for static ranges
template<class Index, class Begin, class End, class F, class... Args,
    std::enable_if_t<IsIntegralConstant<Begin>::value and IsIntegralConstant<End>::value, int> = 0>
void integralRangeFor(Begin&& begin, End&& end, F&& f, Args&&... args)
{
  static const Index begin_t = begin;
  static const Index end_t = end;
  StaticForLoop<Index, begin_t, end_t>::apply(std::forward<F>(f), std::forward<Args>(args)...);
}

// Overload for dynamic ranges
template<class Index, class Begin, class End, class F, class... Args,
    std::enable_if_t<not(IsIntegralConstant<Begin>::value and IsIntegralConstant<End>::value), int> = 0>
void integralRangeFor(Begin&& begin, End&& end, F&& f, Args&&... args)
{
  for(Index i=begin; i != end; ++i)
    f(i, std::forward<Args>(args)...);
}

}

/**
 * \brief Hybrid for loop over integral range
 *
 * \tparam Index Raw type of used indices
 * \tparam Begin Type of begin index
 * \tparam End Type of end index
 * \tparam F Type of functor containing the loop body
 * \tparam Args Types of further arguments to the loop body
 *
 * \param begin Initial index
 * \param end One past last index
 * \param f Functor to call in each loop instance
 * \param args Additional arguments to be passed to the functor
 *
 * This is a hybrid for loop that can work on statically and dynamically
 * sized containers. The functor is called with index as first argument
 * and all additional arguments. If begin and end are both of type
 * std::integral_constant<*,*> than the loop is static with indices
 * of the form std::integral_constant<Index, *>, otherwise the loop
 * is dynamic with indices type Index.
 */
template<class Index, class Begin, class End, class F, class... Args>
void integralRangeFor(Begin&& begin, End&& end, F&& f, Args&&... args)
{
  Imp::integralRangeFor<Index>(std::forward<Begin>(begin), std::forward<End>(end), std::forward<F>(f), std::forward<Args>(args)...);
}



// Implementation of hybridEquals
namespace Imp {

// Compute t1==t2 either statically or dynamically
template<class T1, class T2>
constexpr auto hybridEquals(const T1& t1, const T2& t2, PriorityTag<1>) -> decltype(T1::value, T2::value, std::integral_constant<bool,T1::value == T2::value>())
{ return {}; }

template<class T1, class T2>
constexpr auto hybridEquals(const T1& t1, const T2& t2, PriorityTag<0>)
{
  return t1==t2;
}

} //end namespace Imp

/**
 * \brief Hybrid equality comparison
 *
 * If both types have a static member value, the result of comparing
 * these is returned as std::integral_constant<bool, *>. Otherwise
 * the result of a runtime comparison of t1 and t2 is directly returned.
 */
template<class T1, class T2>
constexpr auto hybridEquals(const T1& t1, const T2& t2)
{
  return Imp::hybridEquals(t1, t2, PriorityTag<1>());
}



// Implementation of hybridIf
namespace Imp {

template<class IfFunc, class ElseFunc>
constexpr void hybridIf(std::true_type, IfFunc&& ifFunc, ElseFunc&& elseFunc)
{
  ifFunc([](auto&& x) { return std::forward<decltype(x)>(x);});
}

template<class IfFunc, class ElseFunc>
constexpr void hybridIf(std::false_type, IfFunc&& ifFunc, ElseFunc&& elseFunc)
{
  elseFunc([](auto&& x) { return std::forward<decltype(x)>(x);});
}

template<class IfFunc, class ElseFunc>
constexpr void hybridIf(const bool& condition, IfFunc&& ifFunc, ElseFunc&& elseFunc)
{
  if (condition)
    ifFunc([](auto&& x) { return std::forward<decltype(x)>(x);});
  else
    elseFunc([](auto&& x) { return std::forward<decltype(x)>(x);});
}

} //end namespace Imp

/**
 * \brief Hybrid if
 *
 * This will call either ifFunc or elseFunc depending
 * on the condition. In any case a single argument
 * will be passed to the called function. This will always
 * be the indentity function. Passing an expression through
 * this function will lead to lazy evaluation. This way both
 * 'branches' can contain expressions that are only valid
 * within this branch if the condition is a std::integral_constant<bool,*>.
 *
 * In order to do this, the passed functors must have a single
 * argument of type auto.
 */
template<class Condition, class IfFunc, class ElseFunc>
constexpr void hybridIf(const Condition& condition, IfFunc&& ifFunc, ElseFunc&& elseFunc)
{
  Imp::hybridIf(condition, std::forward<IfFunc>(ifFunc), std::forward<ElseFunc>(elseFunc));
}

/**
 * \brief Hybrid if
 *
 * This provides a hybridIf with empty else clause.
 */
template<class Condition, class IfFunc>
constexpr void hybridIf(const Condition& condition, IfFunc&& ifFunc)
{
  hybridIf(condition, std::forward<IfFunc>(ifFunc), [](auto&& i) {});
}



// Everything in the next namespace block is just used to implement StaticSize, HasStaticSize, hybridSize
namespace Imp {

// As a last resort try if there's a static constexpr size()
template<class T>
constexpr auto staticSize(const T*, const PriorityTag<0>&)
  -> decltype(std::integral_constant<std::size_t,T::size()>())
{
  return {};
}

// Try if tuple_size is implemented for class
template<class T>
constexpr auto staticSize(const T*, const PriorityTag<2>&)
  -> decltype(std::integral_constant<std::size_t,std::tuple_size<T>::value>())
{
  return {};
}

// Try if tuple_size is implemented for class
template<class T, int i>
constexpr auto staticSize(const Dune::FieldVector<T, i>*, const PriorityTag<3>&)
  -> decltype(std::integral_constant<std::size_t,i>())
{
  return {};
}

template<class T>
constexpr std::false_type hasStaticSize(const T* t, const PriorityTag<0>& p)
{
  return {};
}

template<class T>
constexpr auto hasStaticSize(const T* t, const PriorityTag<1>& p)
  -> decltype(staticSize(t ,PriorityTag<42>()), std::true_type())
{
  return {};
}

}



/**
 * \brief Check if type is a statically sized container
 *
 * \ingroup Utility
 *
 * Derives from std::true_type or std::false_type
 */
template<class T>
struct HasStaticSize :
  public decltype(Imp::hasStaticSize((typename std::decay<T>::type*)(nullptr), PriorityTag<42>()))
{};



/**
 * \brief Obtain size of statically sized container
 *
 * \ingroup Utility
 *
 * Derives from std::integral_constant<std::size_t, size>
 */
template<class T>
struct StaticSize :
  public decltype(Imp::staticSize((typename std::decay<T>::type*)(nullptr), PriorityTag<42>()))
{};



/**
 * \brief Hybrid size query
 *
 * \tparam T Type of container whose size is queried
 *
 * \param t Container whose size is queried
 *
 * This function is hybrid in the sense that it returns a statically
 * encoded size, i.e., an integral_constant if possible and the
 * dynamic result of the t.size() method otherwise.
 *
 * This is the static-size overload which returns the size i
 * as std::integral_constant<std::size_t, i>.
 */
template<class T,
    std::enable_if_t<HasStaticSize<T>::value, int> = 0>
constexpr auto hybridSize(const T& t)
{
  return Imp::staticSize((T*)(nullptr), PriorityTag<42>());
}

/**
 * \brief Hybrid size query
 *
 * \tparam T Type of container whose size is queried
 *
 * \param t Container whose size is queried
 *
 * This function is hybrid in the sense that it returns a statically
 * encoded size, i.e., an integral_constant if possible and the
 * dynamic result of the *.size() method otherwise.
 *
 * This is the dynamic-size overload which returns the result
 * of t.size().
 */
template<class T,
    std::enable_if_t<not HasStaticSize<T>::value, int> = 0>
constexpr auto hybridSize(const T& t)
{
  return t.size();
}



/**
 * \brief Hybrid for loop over sparse range
 */
template<class... T, class F>
void sparseRangeFor(const Dune::MultiTypeBlockVector<T...>& range, F&& f)
{
  integralRangeFor<std::size_t>(Indices::_0, hybridSize(range), [&](auto&& i) {
      f(range[i], i);
  });
}



/**
 * \brief Hybrid for loop over sparse range
 */
template<class Range, class F>
void sparseRangeFor(Range&& range, F&& f)
{
  auto it = range.begin();
  auto end = range.end();
  for(; it!=end; ++it)
    f(*it, it.index());
}



} // namespace Solvers
} // namespace Dune


#endif// DUNE_SOLVERS_COMMON_FORLOOP_HH