operationtraits.hh

#ifndef __DUNE_ACFEM_TENSORS_OPERATIONTRAITS_HH__
#define __DUNE_ACFEM_TENSORS_OPERATIONTRAITS_HH__
 
#include "../mpl/compare.hh"
#include "../expressions/storage.hh"
#include "../expressions/expressionoperations.hh"
#include "expressionoperations.hh"
#include "operations/einsum.hh"
#include "operations/product.hh"
#include "operations/restriction.hh"
#include "operations/transpose.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<std::size_t... Dims, std::size_t... PivotIndices>
    struct OperationTraits<RestrictionOperation<IndexSequence<Dims...>, IndexSequence<PivotIndices...> > >
    {
      static_assert(sizeof...(Dims) == sizeof...(PivotIndices),
                    "Size of pivot indices must match size of co-dimension.");
 
      using PivotSequence = IndexSequence<PivotIndices...>;
      using InvertibleOperation = void; // restrictions are not invertible
      using OperationType = RestrictionOperation<IndexSequence<Dims...>, PivotSequence>;
 
      template<class T, class SFINAE = void>
      struct HasOperation
        : FalseType
      {
        // assume elementary type if no restriction operation is defined.
        using ResultType = T;
      };
 
      template<class T>
      struct HasOperation<T, std::enable_if_t<(sizeof(decltype(restriction<Dims...>(std::declval<T>(), PivotSequence{}))) >= 0)> >
        : TrueType
      {
        using ResultType = std::decay_t<decltype(restriction<Dims...>(std::declval<T>(), PivotSequence{}))>;
      };
 
      template<class T>
      using ResultType = typename HasOperation<T>::ResultType;
 
      // Assuming only real numbers.
      template<class Sign>
      static constexpr auto signPropagation(const Sign&)
      {
        return Sign{};
      }
 
      static std::string name()
      {
        return "rst<"+toString(IndexSequence<Dims...>{})+"@"+toString(PivotSequence{})+">";
      }
 
      template<class Order>
      static constexpr auto polynomialOrder(Order oldOrder)
      {
        return oldOrder;
      }
 
      template<class T, class... Rest, std::enable_if_t<HasOperation<T>::value, int> = 0>
      auto operator()(T&& t, Rest&&...) const
      {
        return restriction<Dims...>(std::forward<T>(t), PivotSequence{});
      }
 
      template<class T, class... Rest, std::enable_if_t<!HasOperation<T>::value, int> = 0>
      constexpr decltype(auto) operator()(T&& t, Rest&&...) const
      {
        return std::forward<T>(t);
      }
 
      using PivotType = PivotSequence;
      PivotType pivotIndices_;
    };
 
    template<std::size_t... Dims>
    struct OperationTraits<RestrictionOperation<IndexSequence<Dims...>, IndexSequence<> > >
    {
     private:
      using ThisType = OperationTraits;
     public:
      using InvertibleOperation = void; // restrictions are not invertible
      using OperationType = RestrictionOperation<IndexSequence<Dims...>, IndexSequence<> >;
 
      OperationTraits()
        : pivotIndices_({{}})
      {}
 
      OperationTraits(const OperationTraits& other)
        : pivotIndices_(other.pivotIndices_)
      {}
 
      OperationTraits(OperationTraits&& other)
        : pivotIndices_(std::move(other.pivotIndices_))
      {}
 
      template<class List, std::size_t... I, std::enable_if_t<IsTupleLike<List>::value, int> = 0>
      OperationTraits(const List& l, IndexSequence<I...>)
        : pivotIndices_({{ (IndexType)get<I>(l)... }})
      {}
 
      template<class List, std::enable_if_t<IsTupleLike<List>::value, int> = 0>
      OperationTraits(const List& l)
        : OperationTraits(l, MakeSequenceFor<List>{})
      {}
 
      template<class T, std::size_t N, std::size_t... I>
      OperationTraits(const T (&l)[N], IndexSequence<I...>)
        : pivotIndices_({{ (IndexType)get<I>(l)... }})
      {}
 
      template<class I, std::size_t N>
      OperationTraits(const I (&l)[N])
        : OperationTraits(l, MakeIndexSequence<N>{})
      {}
 
      template<class T, class SFINAE = void>
      struct HasOperation
        : FalseType
      {
        // assume elementary type if no restriction operation is defined.
        using ResultType = T;
      };
 
      template<class T>
      struct HasOperation<T, std::enable_if_t<sizeof(decltype(restriction<Dims...>(std::declval<T>()))) >= 0> >
        : TrueType
      {
        using ResultType = std::decay_t<decltype(restriction<Dims...>(std::declval<T>()))>;
      };
 
      template<class T>
      using ResultType = typename HasOperation<T>::ResultType;
 
      // Assuming only real numbers.
      template<class Sign>
      static constexpr auto signPropagation(const Sign&)
      {
        return Sign{};
      }
 
      std::string name() const
      {
        return "rst<"+toString(IndexSequence<Dims...>{})+"@["+toString(pivotIndices_)+"]>";
      }
 
      template<class Order>
      static constexpr auto polynomialOrder(Order oldOrder)
      {
        return oldOrder;
      }
 
     public:
      template<class T, class... Rest, std::enable_if_t<HasOperation<T>::value, int> = 0>
      constexpr auto operator()(T&& t, Rest&&...) const
      {
        return restriction<Dims...>(std::forward<T>(t), pivotIndices_);
      }
 
      template<class T, class... Rest, std::enable_if_t<!HasOperation<T>::value, int> = 0>
      constexpr decltype(auto) operator()(T&& t, Rest&&...) const
      {
        return std::forward<T>(t);
      }
 
      using IndexType = unsigned;
      using PivotType = std::array<IndexType, sizeof...(Dims)>;
      PivotType pivotIndices_;
    };
 
 
    template<class Perm>
    struct OperationTraits<TransposeOperation<Perm> >
    {
      using InvertibleOperation = TransposeOperation<InversePermutation<Perm> >;
      using OperationType = TransposeOperation<Perm>;
 
      template<class T, class SFINAE = void>
      struct HasOperation
        : FalseType
      {
        // assume elementary type if no restriction operation is defined.
        using ResultType = T;
      };
 
      template<class T>
      struct HasOperation<T, VoidType<decltype(transpose(std::declval<T>(), Perm{}))> >
        : TrueType
      {
        using ResultType = std::decay_t<decltype(transpose(std::declval<T>(), Perm{}))>;
      };
 
      // Assuming only real numbers.
      template<class Sign>
      static constexpr auto signPropagation(const Sign&)
      {
        return Sign{};
      }
 
      static std::string name()
      {
        return "tr<"+toString(Perm{})+">";
      }
 
      template<class Order>
      static constexpr auto polynomialOrder(Order oldOrder)
      {
        return oldOrder;
      }
 
      template<class T>
      auto operator()(T&& t) const
      {
        return transpose(std::forward<T>(t), Perm{});
      }
 
      template<class T, std::enable_if_t<!HasOperation<T>::value, int> = 0>
      constexpr decltype(auto) operator()(T&& t) const
      {
#if 0
        if (!IsScalar<T>::value) {
          std::clog << "no transpose " << toString(Perm{}) << " for " << t.name() << std::endl;
        }
#endif
        assert(IsScalar<T>::value);
        return std::forward<T>(t);
      }
    };
 
 
    template<class Sig>
    struct OperationTraits<ReshapeOperation<Sig> >
    {
      using InvertibleOperation = void;
      using OperationType = ReshapeOperation<Sig>;
 
      template<class T, class SFINAE = void>
      struct HasOperation
        : FalseType
      {
        // assume elementary type if no restriction operation is defined.
        using ResultType = T;
      };
 
      template<class T>
      struct HasOperation<T, VoidType<decltype(reshape(std::declval<T>(), Sig{}))> >
        : TrueType
      {
        using ResultType = std::decay_t<decltype(reshape(std::declval<T>(), Sig{}))>;
      };
 
      // Assuming only real numbers.
      template<class Sign>
      static constexpr auto signPropagation(const Sign&)
      {
        return Sign{};
      }
 
      static std::string name()
      {
        return "reshape<"+toString(Sig{})+">";
      }
 
      template<class Order>
      static constexpr auto polynomialOrder(Order oldOrder)
      {
        return oldOrder;
      }
 
      template<class T>
      auto operator()(T&& t) const
      {
        return reshape(std::forward<T>(t), Sig{});
      }
 
      template<class T, std::enable_if_t<!HasOperation<T>::value, int> = 0>
      constexpr decltype(auto) operator()(T&& t) const
      {
        assert(IsScalar<T>::value);
        return std::forward<T>(t);
      }
    };
 
 
    template<class Pos1, class Pos2, class Dims>
    struct OperationTraits<EinsumOperation<Pos1, Pos2, Dims> >
    {
      static_assert(Pos1::size() == Pos2::size() && Pos2::size() == Dims::size(),
                    "Number of contraction indices must coincide for Einstein-summation.");
 
      using InvertibleOperation = void;
      using OperationType = EinsumOperation<Pos1, Pos2, Dims>;
 
      using LeftIndexPositions = Pos1;
      using RightIndexPositions = Pos2;
      using Dimensions = Dims;
 
      template<class T0, class T1>
      using Signature = SequenceCat<SequenceSliceComplement<typename TensorTraits<T0>::Signature, LeftIndexPositions>,
                                    SequenceSliceComplement<typename TensorTraits<T1>::Signature, RightIndexPositions> >;
 
      static constexpr Dimensions dimensions()
      {
        return Dimensions{};
      }
 
      static constexpr std::size_t contractionDimension()
      {
        return prod(dimensions());
      }
 
      template<class T1, class T2, class SFINAE = void>
      struct HasOperation
        : FalseType
      {
        // assume elementary type if no restriction operation is defined.
        using ResultType = std::decay_t<decltype(std::declval<T1>() * std::declval<T2>())>;
      };
 
      template<class T1, class T2>
      struct HasOperation<T1, T2, VoidType<decltype(einsum<Pos1, Pos2>(std::declval<T1>(), std::declval<T2>()))> >
        : TrueType
      {
        using ResultType = std::decay_t<decltype(einsum<Pos1, Pos2>(std::declval<T1>(), std::declval<T2>()))>;
      };
 
      template<class T1, class T2>
      using ResultType = typename HasOperation<T1, T2>::ResultType;
 
      template<class Order1, class Order2>
      static constexpr auto polynomialOrder(Order1 order1, Order2 order2)
      {
        return order1 + order2;
      }
 
      static constexpr std::size_t dim = Prod<Dims>::value;
 
      // in arbitraty dimension this is a scalar product, hence the
      // nonZero property need not be maintained.
      template<class S1, class S2>
      static constexpr auto signPropagation(const S1&, const S2&)
      {
        return ExpressionSign<// non-singular
          (dim == 1
           &&
           S1::isNonSingular && S2::isNonSingular
           &&
           (((S1::reference >= 0) && S1::isSemiPositive)
            ||
            ((S1::reference <= 0) && S1::isSemiNegative))
           &&
           (((S2::reference >= 0) && S2::isSemiPositive)
            ||
            ((S2::reference <= 0) && S2::isSemiNegative))),
          // >= dim * S1::reference * S2::reference
          ((S1::reference*S2::reference >= 0)
           &&
           ((S1::isSemiPositive && S2::isSemiPositive)
            ||
            (S1::isSemiNegative && S2::isSemiNegative))),
          // <= dim * S1::reference * S2::reference
          ((S1::reference*S2::reference <= 0)
           &&
           ((S1::isSemiPositive && S2::isSemiNegative)
            ||
            (S1::isSemiNegative && S2::isSemiPositive))),
          // new reference
          dim * S1::reference * S2::reference>{};
      }
 
      static std::string name()
      {
        return "*"+(Pos1::size() > 0 ? "["+toString(Pos1{})+"]["+toString(Pos2{})+"]" : "");
      }
 
      template<class T1, class T2, std::enable_if_t<HasOperation<T1, T2>::value, int> = 0>
      constexpr auto operator()(T1&& t1, T2&& t2) const
      {
        return einsum<Pos1, Pos2>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      // degenerate to ordinary product for elementary operands
      template<class T1, class T2, std::enable_if_t<!HasOperation<T1, T2>::value, int> = 0>
      constexpr auto operator()(T1&& t1, T2&& t2) const
      {
        return std::forward<T1>(t1) * std::forward<T2>(t2);
      }
    };
 
 
    template<class Pos1, class Pos2, class Dims>
    struct OperationTraits<TensorProductOperation<Pos1, Pos2, Dims> >
    {
      static_assert(Pos1::size() == Pos2::size() && Pos2::size() == Dims::size(),
                    "Number of contraction indices must coincide for product operation.");
 
      using InvertibleOperation = void;
      using OperationType = TensorProductOperation<Pos1, Pos2, Dims>;
 
      using LeftIndexPositions = Pos1;
      using RightIndexPositions = Pos2;
      using Dimensions = Dims;
 
      template<class T1, class T2, class SFINAE = void>
      struct HasOperation
        : FalseType
      {
        // assume elementary type if no restriction operation is defined.
        using ResultType = std::decay_t<decltype(std::declval<T1>() * std::declval<T2>())>;
      };
 
      template<class T1, class T2>
      struct HasOperation<T1, T2, std::enable_if_t<(sizeof(decltype(multiply<Pos1, Pos2>(std::declval<T1>(), std::declval<T2>()))) >= 0)> >
        : TrueType
      {
        using ResultType = std::decay_t<decltype(multiply<Pos1, Pos2>(std::declval<T1>(), std::declval<T2>()))>;
      };
 
      template<class T1, class T2>
      using ResultType = typename HasOperation<T1, T2>::ResultType;
 
      template<class Order1, class Order2>
      static constexpr auto polynomialOrder(Order1 order1, Order2 order2)
      {
        return order1 + order2;
      }
 
      static constexpr std::size_t dim = 1; // pointwise, no summation
 
      // in arbitraty dimension this is a scalar product, hence the
      // nonZero property need not be maintained.
      template<class S1, class S2>
      static constexpr auto signPropagation(const S1&, const S2&)
      {
        return ExpressionSign<// non-singular
          (dim == 1
           &&
           S1::isNonSingular && S2::isNonSingular
           &&
           (((S1::reference >= 0) && S1::isSemiPositive)
            ||
            ((S1::reference <= 0) && S1::isSemiNegative))
           &&
           (((S2::reference >= 0) && S2::isSemiPositive)
            ||
            ((S2::reference <= 0) && S2::isSemiNegative))),
          // >= dim * S1::reference * S2::reference
          ((S1::reference*S2::reference >= 0)
           &&
           ((S1::isSemiPositive && S2::isSemiPositive)
            ||
            (S1::isSemiNegative && S2::isSemiNegative))),
          // <= dim * S1::reference * S2::reference
          ((S1::reference*S2::reference <= 0)
           &&
           ((S1::isSemiPositive && S2::isSemiNegative)
            ||
            (S1::isSemiNegative && S2::isSemiPositive))),
          // new reference
          dim * S1::reference * S2::reference>{};
      }
 
      static std::string name()
      {
        return "."+(Pos1::size() > 0 ? "["+toString(Pos1{})+"]["+toString(Pos2{})+"]" : "");
      }
 
      template<class T1, class T2>
      auto operator()(T1&& t1, T2&& t2) const
      {
        return multiply<Pos1, Pos2>(std::forward<T1>(t1), std::forward<T2>(t2));
      }
 
      // degenerate to ordinary product for elementary operands
      template<class T1, class T2, std::enable_if_t<!HasOperation<T1, T2>::value, int> = 0>
      constexpr auto operator()(T1&& t1, T2&& t2) const
      {
        return std::forward<T1>(t1) * std::forward<T2>(t2);
      }
    };
 
    template<class T>
    struct SquareTraits<T, std::enable_if_t<(HasTensorTraitsV<T> && TensorTraits<T>::rank > 0)> >
    {
      using ExplodeFunctor = OperationTraits<
        TensorProductOperation<
          MakeIndexSequence<TensorTraits<T>::rank>,
          MakeIndexSequence<TensorTraits<T>::rank>,
          typename TensorTraits<T>::Signature>
        >;
    };
 
    template<class T>
    struct SquareTraits<T, std::enable_if_t<(HasTensorTraitsV<T> && TensorTraits<T>::rank == 0)> >
    {
      using ExplodeFunctor = Tensor::ScalarEinsumFunctor;
    };
 
 
 
 
  } // NS ACFem
 
} // NS Dune
 
#endif // __DUNE_ACFEM_TENSORS_OPERATIONTRAITS_HH__