associativity.hh

#ifndef __DUNE_ACFEM_TENSORS_OPERATIONS_ASSOCIATIVITY_HH__
#define __DUNE_ACFEM_TENSORS_OPERATIONS_ASSOCIATIVITY_HH__
 
#include "../tensorbase.hh"
#include "../expressionoperations.hh"
#include "einsum.hh"
#include "product.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    namespace Tensor {
 
      namespace ProductOperations {
 
        template<class F, class T0, class T1, class SFINAE = void>
        struct IsRightAssociative
          : FalseType
        {};
 
        template<class F, class T0, class T1>
        struct IsRightAssociative<F, T0, T1, std::enable_if_t<!IsDecay<F>::value> >
          : IsRightAssociative<std::decay_t<F>, T0, T1>
        {};
 
        template<class Seq0, class Seq1, class Dims, class T0, class T1>
        struct IsRightAssociative<OperationTraits<EinsumOperation<Seq0, Seq1, Dims> >, T0, T1>
          : FunctorHas<IsEinsumOperation, Functor<T0> >
        {};
 
        template<class Seq0, class Seq1, class Dims, class T0, class T1>
        struct IsRightAssociative<
          OperationTraits<TensorProductOperation<Seq0, Seq1, Dims> >, T0, T1,
          std::enable_if_t<(FunctorHas<IsTensorProductOperation, Functor<T0> >::value
                            && isSimple(Seq0{})
                            && std::is_same<Seq0, Seq1>::value
                            && std::is_same<typename TensorProductTraits<Operation<T0> >::LeftIndexPositions,
                                            typename TensorProductTraits<Operation<T0> >::RightIndexPositions>::value
                            && std::is_same<Seq0,
                                            typename TensorProductTraits<Operation<T0> >::RightIndexPositions>::value
            )> >
          : TrueType
        {};
 
        template<class F, class T0, class T1, class SFINAE = void>
        struct IsLeftAssociative
          : FalseType
        {};
 
        template<class F, class T0, class T1>
        struct IsLeftAssociative<F, T0, T1, std::enable_if_t<!IsDecay<F>::value> >
          : IsLeftAssociative<std::decay_t<F>, T0, T1>
        {};
 
        template<class Seq0, class Seq1, class Dims, class T0, class T1>
        struct IsLeftAssociative<OperationTraits<EinsumOperation<Seq0, Seq1, Dims> >, T0, T1>
          : FunctorHas<IsEinsumOperation, Functor<T1> >
        {};
 
        template<class Seq0, class Seq1, class Dims, class T0, class T1>
        struct IsLeftAssociative<
          OperationTraits<TensorProductOperation<Seq0, Seq1, Dims> >, T0, T1,
          std::enable_if_t<(FunctorHas<IsTensorProductOperation, Functor<T1> >::value
                            && isSimple(Seq0{})
                            && std::is_same<Seq0, Seq1>::value
                            && std::is_same<typename TensorProductTraits<Operation<T1> >::LeftIndexPositions,
                                            typename TensorProductTraits<Operation<T1> >::RightIndexPositions>::value
                            && std::is_same<Seq0,
                                            typename TensorProductTraits<Operation<T1> >::RightIndexPositions>::value
            )> >
          : TrueType
        {};
 
        template<class InnerLeftIndexPositions, class InnerRightIndexPositions, class InnerDims,
                 class LeftIndexPositions, class RightIndexPositions, class OuterDims,
                 std::size_t Rank1, std::size_t Rank2, std::size_t Rank3,
                 std::enable_if_t<Rank1 * Rank2 * Rank3 != 0, int> = 0>
        constexpr auto associateRightOperations(
          OperationTraits<EinsumOperation<InnerLeftIndexPositions, InnerRightIndexPositions, InnerDims> >,
          OperationTraits<EinsumOperation<LeftIndexPositions, RightIndexPositions, OuterDims> >,
          IndexSequence<Rank1, Rank2, Rank3>)
        {
          // this ain't pretty ...
 
          constexpr std::size_t innerDefectRank = InnerDims::size();
          constexpr std::ptrdiff_t rightRank = Rank3;
 
          constexpr std::ptrdiff_t innerLeftRank = Rank1;
          constexpr std::ptrdiff_t innerLeftDefectRank = innerLeftRank - innerDefectRank;
          constexpr std::ptrdiff_t innerRightRank = Rank2;
          constexpr std::ptrdiff_t innerRightDefectRank = innerRightRank - innerDefectRank;
 
          // Indices of the outer contraction positions referring to type One and Two
          using InnerLeftPosIndices = TransformSequence<LeftIndexPositions, Seq<>, IndexFunctor, AcceptInputInRangeFunctor<std::size_t, 0, innerLeftDefectRank> >;
          using InnerRightPosIndices = TransformSequence<LeftIndexPositions, Seq<>, IndexFunctor, AcceptInputInRangeFunctor<std::size_t, innerLeftDefectRank, innerLeftDefectRank+innerRightDefectRank> >;
 
          // Parts of the outer left contraction positions referring to type One and Two
          using InnerLeftPos = SequenceSlice<LeftIndexPositions, InnerLeftPosIndices>;
          using InnerRightPos = OffsetSequence<-innerLeftDefectRank, SequenceSlice<LeftIndexPositions, InnerRightPosIndices> >;
 
          // Parts of the outer right contraction positions referring to type Two and Three
          using OuterLeftPos = SequenceSlice<RightIndexPositions, InnerLeftPosIndices>;
          using OuterRightPos = SequenceSlice<RightIndexPositions, InnerRightPosIndices>;
 
          // Lookup table to convert to original indices for left (inner) tensors
          using InnerLeftLookup = SequenceSliceComplement<MakeIndexSequence<innerLeftRank>, InnerLeftIndexPositions>;
          using InnerRightLookup = SequenceSliceComplement<MakeIndexSequence<innerRightRank>, InnerRightIndexPositions>;
 
          // Map the outer contraction indices to the original indices
          using InnerLeftMappedPos = TransformedSequence<MapSequenceFunctor<InnerLeftLookup>, InnerLeftPos>;
          using InnerRightMappedPos = TransformedSequence<MapSequenceFunctor<InnerRightLookup>, InnerRightPos>;
 
          // Contraction indices for One. Concat the inner contraction
          // indices with the mapped outer contraction indices
          // referring to one
          using OnePos = SequenceCat<InnerLeftIndexPositions, InnerLeftMappedPos>;
          using OneDims = SequenceCat<InnerDims, SequenceSlice<OuterDims, InnerLeftPosIndices> >;
 
          // Contraction indices for Two. Map the indices referring to
          // Two to the original indices.
          using TwoPos = InnerRightMappedPos;
 
          // contraction position for Three which refer to B
          using ThreePos = OuterRightPos;
          using ThreeDims = SequenceSlice<OuterDims, InnerRightPosIndices>;
 
          // Contraction indices for the contraction of One with (Two
          // * Three). This refers to the contraction of One with
          // Two. We need the indices referring to Two in the
          // contraction of (Two * Three).
          using TwoThreeTwoPosLookup = SequenceSliceComplement<MakeIndexSequence<innerRightRank>, TwoPos>;
          using TwoThreeTwoMappedPos = TransformedSequence<InverseMapSequenceFunctor<TwoThreeTwoPosLookup>, InnerRightIndexPositions>;
 
          using TwoThreeThreePosLookup = SequenceSliceComplement<MakeIndexSequence<rightRank>, ThreePos>;
          using TwoThreeThreeMappedPos = TransformedSequence<InverseMapSequenceFunctor<TwoThreeThreePosLookup>, OuterLeftPos>;
 
          using TwoThreePos = SequenceCat<TwoThreeTwoMappedPos, OffsetSequence<innerRightRank-TwoPos::size() ,TwoThreeThreeMappedPos> >;
 
#if 0
          std::clog << "One: " << OnePos{} << std::endl;
          std::clog << "TwoThree: " << TwoThreePos{} << std::endl;
          std::clog << "Two: " << TwoPos{} << std::endl;
          std::clog << "Three: " << ThreePos{} << std::endl;
#endif
 
          return std::make_pair(OperationTraits<EinsumOperation<OnePos, TwoThreePos, OneDims> >{},
                                OperationTraits<EinsumOperation<TwoPos, ThreePos, ThreeDims> >{});
        }
 
        template<class InnerLeftIndexPositions, class InnerRightIndexPositions, class InnerDims,
                 class LeftIndexPositions, class RightIndexPositions, class OuterDims,
                 std::size_t Rank1, std::size_t Rank2, std::size_t Rank3,
                 std::enable_if_t<Rank1 * Rank2 * Rank3 == 0, int> = 0>
        constexpr auto associateRightOperations(
          OperationTraits<EinsumOperation<InnerLeftIndexPositions, InnerRightIndexPositions, InnerDims> > f1,
          OperationTraits<EinsumOperation<LeftIndexPositions, RightIndexPositions, OuterDims> > f2,
          IndexSequence<Rank1, Rank2, Rank3>)
        {
          // (A *1 B) *2 C
          // A scalar: A *1 (B *2 C)
          // B scalar: A *2 (B *1 C)
          // C scalar: A *1 (B *2 C)
          if constexpr (Rank1 == 0) {
            return std::make_pair(f1, f2);
          } else if constexpr (Rank2 == 0) {
            return std::make_pair(f2, f1);
          } else {
            return std::make_pair(f1, f2);
          }
        }
 
        template<class LeftIndexPositions, class RightIndexPositions, class OuterDims,
                 class InnerLeftIndexPositions, class InnerRightIndexPositions, class InnerDims,
                 std::size_t Rank1, std::size_t Rank2, std::size_t Rank3,
                 std::enable_if_t<Rank1 * Rank2 * Rank3 != 0, int> = 0>
        constexpr auto associateLeftOperations(
          OperationTraits<EinsumOperation<LeftIndexPositions, RightIndexPositions, OuterDims> >,
          OperationTraits<EinsumOperation<InnerLeftIndexPositions, InnerRightIndexPositions, InnerDims> >,
          IndexSequence<Rank1, Rank2, Rank3>)
        {
          // this ain't pretty ...
          constexpr std::size_t innerDefectRank = InnerDims::size();
          constexpr std::ptrdiff_t leftRank = Rank1;
 
          constexpr std::ptrdiff_t innerLeftRank = Rank2;
          constexpr std::ptrdiff_t innerLeftDefectRank = innerLeftRank - innerDefectRank;
          constexpr std::ptrdiff_t innerRightRank = Rank3;
          constexpr std::ptrdiff_t innerRightDefectRank = innerRightRank - innerDefectRank;
 
          // Indices of the outer contraction positions referring to type two and three
          using InnerLeftPosIndices = TransformSequence<RightIndexPositions, Seq<>, IndexFunctor, AcceptInputInRangeFunctor<std::size_t, 0, innerLeftDefectRank> >;
          using InnerRightPosIndices = TransformSequence<RightIndexPositions, Seq<>, IndexFunctor, AcceptInputInRangeFunctor<std::size_t, innerLeftDefectRank, innerLeftDefectRank+innerRightDefectRank> >;
 
          // Parts of the outer right contraction positions referring to type Two and Three
          using InnerLeftPos = SequenceSlice<RightIndexPositions, InnerLeftPosIndices>;
          using InnerRightPos = OffsetSequence<-innerLeftDefectRank, SequenceSlice<RightIndexPositions, InnerRightPosIndices> >;
 
          // Parts of the outer left contraction positions referring to type Two and Three
          using OuterLeftPos = SequenceSlice<LeftIndexPositions, InnerLeftPosIndices>;
          using OuterRightPos = SequenceSlice<LeftIndexPositions, InnerRightPosIndices>;
 
          // Lookup table to convert to original indices for left (inner) tensors
          using InnerLeftLookup = SequenceSliceComplement<MakeIndexSequence<innerLeftRank>, InnerLeftIndexPositions>;
          using InnerRightLookup = SequenceSliceComplement<MakeIndexSequence<innerRightRank>, InnerRightIndexPositions>;
 
          // Map the outer contraction indices to the original indices
          using InnerLeftMappedPos = TransformedSequence<MapSequenceFunctor<InnerLeftLookup>, InnerLeftPos>;
          using InnerRightMappedPos = TransformedSequence<MapSequenceFunctor<InnerRightLookup>, InnerRightPos>;
 
          // Part of left-most contraction referring to two
          using OnePos = OuterLeftPos;
          using OneDims = SequenceSlice<OuterDims, InnerLeftPosIndices>;
 
          // Contraction indices for Two.
          using TwoPos = InnerLeftMappedPos;
 
          // Contraction indices for Three. Concat the inner contraction
          // indices with the mapped outer contraction indices
          // referring to two and sort the result.
          using ThreePos = SequenceCat<InnerRightMappedPos, InnerRightIndexPositions>;
          using ThreeDims = SequenceCat<SequenceSlice<OuterDims, InnerRightPosIndices>, InnerDims>;
 
          // Contraction indices for the contraction of (One*Two) with
          // Three. This refers to the contraction of One with Three
          // and Two with Three. We need the indices referring to Two
          // in the contraction of (One * Two).
 
          // indices for A
          using OneTwoOnePosLookup = SequenceSliceComplement<MakeIndexSequence<leftRank>, OnePos>;
          using OneTwoOneMappedPos = TransformedSequence<InverseMapSequenceFunctor<OneTwoOnePosLookup>, OuterRightPos>;
 
          // indices for B
          using OneTwoTwoPosLookup = SequenceSliceComplement<MakeIndexSequence<innerLeftRank>, TwoPos>;
          using OneTwoTwoMappedPos = TransformedSequence<InverseMapSequenceFunctor<OneTwoTwoPosLookup>, InnerLeftIndexPositions>;
 
          using OneTwoPos = SequenceCat<OneTwoOneMappedPos, OffsetSequence<leftRank-OnePos::size(), OneTwoTwoMappedPos> >;
 
          return std::make_pair(
            OperationTraits<EinsumOperation<OnePos, TwoPos, OneDims> >{},
            OperationTraits<EinsumOperation<OneTwoPos, ThreePos, ThreeDims> >{});
        }
 
        template<class LeftIndexPositions, class RightIndexPositions, class OuterDims,
                 class InnerLeftIndexPositions, class InnerRightIndexPositions, class InnerDims,
                 std::size_t Rank1, std::size_t Rank2, std::size_t Rank3,
                 std::enable_if_t<Rank1 * Rank2 * Rank3 == 0, int> = 0>
        constexpr auto associateLeftOperations(
          OperationTraits<EinsumOperation<LeftIndexPositions, RightIndexPositions, OuterDims> > f1,
          OperationTraits<EinsumOperation<InnerLeftIndexPositions, InnerRightIndexPositions, InnerDims> > f2,
          IndexSequence<Rank1, Rank2, Rank3>)
        {
          // A *1 (B *2 C)
          // A scalar: (A *1 B) *2 C
          // B scalar: (A *2 B) *1 C
          // C scalar: (A *1 B) *2 C
          if constexpr (Rank1 == 0) {
            return std::make_pair(f1, f2);
          } else if constexpr (Rank2 == 0) {
            return std::make_pair(f2, f1);
          } else {
            return std::make_pair(f1, f2);
          }
        }
 
        template<class Defects, class Dims,
                 std::size_t Rank1, std::size_t Rank2, std::size_t Rank3>
        constexpr auto associateRightOperations(
          OperationTraits<TensorProductOperation<Defects, Defects, Dims> > f1,
          OperationTraits<TensorProductOperation<Defects, Defects, Dims> > f2,
          IndexSequence<Rank1, Rank2, Rank3>)
        {
          return std::make_pair(f1, f2);
        }
 
        template<class Defects, class Dims,
                 std::size_t Rank1, std::size_t Rank2, std::size_t Rank3>
        constexpr auto associateLeftOperations(
          OperationTraits<TensorProductOperation<Defects, Defects, Dims> > f1,
          OperationTraits<TensorProductOperation<Defects, Defects, Dims> > f2,
          IndexSequence<Rank1, Rank2, Rank3>)
        {
          return std::make_pair(f1, f2);
        }
 
        template<class F, class T0, class T1,
                 class OptimizeInner = Expressions::OptimizeTop,
                 std::enable_if_t<IsRightAssociative<F, T0, T1>::value, int> = 0>
        constexpr auto associateRightExpression(F&& f, T0&& t0, T1&& t1, OptimizeInner = OptimizeInner{})
        {
          auto operations = associateRightOperations(
            std::forward<T0>(t0).operation(), std::forward<F>(f),
            IndexSequence<TensorTraits<Operand<0, T0> >::rank, TensorTraits<Operand<1, T0> >::rank, TensorTraits<T1>::rank>{});
 
          return Expressions::storage(
            std::move(operations.first),
            std::forward<T0>(t0).operand(0_c),
            operate(
              OptimizeInner{},
              std::move(operations.second),
              std::forward<T0>(t0).operand(1_c),
              std::forward<T1>(t1)
              )
            );
        }
 
        template<class F, class T0, class T1,
                 std::enable_if_t<(!IsRightAssociative<F, T0, T1>::value
                                   && FunctorHas<IsProductOperation, F>::value
                                   && IsProductExpression<T1>::value
                   ), int> = 0>
        constexpr auto associateRightExpression(F&& f, T0&& t0, T1&& t1, DontOptimize = DontOptimize{})
        {
          return Expressions::storage(std::forward<F>(f), std::forward<T0>(t0), std::forward<T1>(t1));
        }
 
#if 1
        template<class F, class T0, class T1,
                 class OptimizeOuter = Expressions::OptimizeTop,
                 class OptimizeInner = Expressions::OptimizeTop,
                 std::enable_if_t<AreProperTensors<T0, T1>::value, int> = 0>
        constexpr decltype(auto) associateRight(F&& f, T0&& t0, T1&& t1,
                                                OptimizeOuter = OptimizeOuter{}, OptimizeInner = OptimizeInner{})
        {
          auto expr = associateRightExpression(std::forward<F>(f), std::forward<T0>(t0), std::forward<T1>(t1), OptimizeInner{});
          return operate(
            OptimizeOuter{},
            std::move(expr).operation(),
            std::move(expr).operand(0_c),
            std::move(expr).operand(1_c));
        }
#else
        template<class F, class T0, class T1,
                 class OptimizeOuter = Expressions::OptimizeTop,
                 class OptimizeInner = Expressions::OptimizeTop,
                 std::enable_if_t<(// && AreProperTensors<T0, T1>::value
                                   IsRightAssociative<F, T0, T1>::value), int> = 0>
        constexpr decltype(auto) associateRight(F&& f, T0&& t0, T1&& t1,
                                      OptimizeOuter = OptimizeOuter{}, OptimizeInner = OptimizeInner{})
        {
          auto operations = associateRightOperations(
            std::forward<T0>(t0).operation(), std::forward<F>(f),
            IndexSequence<TensorTraits<Operand<0, T0> >::rank, TensorTraits<Operand<1, T0> >::rank, TensorTraits<T1>::rank>{});
 
          // Do we need forwardReturnValue here?
          return operate(
            OptimizeOuter{},
            std::move(operations.first),
            std::forward<T0>(t0).operand(0_c),
            operate(
              OptimizeInner{},
              std::move(operations.second),
              std::forward<T0>(t0).operand(1_c),
              std::forward<T1>(t1)
              )
            );
        }
 
        template<class F, class T0, class T1,
                 class OptimizeOuter = OptimizeTop,
                 std::enable_if_t<(!IsRightAssociative<F, T0, T1>::value
                                   && FunctorHas<IsProductOperation, F>::value
                                   && IsProductExpression<T1>::value
                   ), int> = 0>
        constexpr decltype(auto) associateRight(F&& f, T0&& t0, T1&& t1,
                                                OptimizeOuter = OptimizeOuter{}, OptimizeTop = OptimizeTop{})
        {
          return operate(OptimizeOuter{}, std::forward<F>(f), std::forward<T0>(t0), std::forward<T1>(t1));
        }
#endif
 
        template<class F, class T0, class T1,
                 class OptimizeInner = OptimizeTop,
                 std::enable_if_t<IsLeftAssociative<F, T0, T1>::value, int> = 0>
        constexpr auto associateLeftExpression(F&& f, T0&& t0, T1&& t1,
                                               OptimizeInner = OptimizeInner{})
        {
          // A*(B*C) -> (A*B)*C
          auto operations = associateLeftOperations(
            std::forward<F>(f), std::forward<T1>(t1).operation(),
            IndexSequence<TensorTraits<T0>::rank, TensorTraits<Operand<0, T1> >::rank, TensorTraits<Operand<1, T1> >::rank>{});
 
          return Expressions::storage(
            std::move(operations.second), // outer
            operate(
              OptimizeInner{},
              std::move(operations.first), // inner
              std::forward<T0>(t0),
              std::forward<T1>(t1).operand(0_c)
              ),
            std::forward<T1>(t1).operand(1_c)
            );
        }
 
        template<class F, class T0, class T1,
                 std::enable_if_t<(!IsLeftAssociative<F, T0, T1>::value
                                   && IsProductExpression<T0>::value
                                   && FunctorHas<IsProductOperation, F>::value
          ), int> = 0>
        constexpr auto associateLeftExpression(F&& f, T0&& t0, T1&& t1,
                                               DontOptimize = DontOptimize{})
        {
          return Expressions::storage(std::forward<F>(f), std::forward<T0>(t0), std::forward<T1>(t1));
        }
 
#if 1
        template<class F, class T0, class T1,
                 class OptimizeOuter = OptimizeTop, class OptimizeInner = OptimizeTop,
                 std::enable_if_t<AreProperTensors<T0, T1>::value, int> = 0>
        constexpr decltype(auto) associateLeft(F&& f, T0&& t0, T1&& t1,
                                               OptimizeOuter = OptimizeOuter{}, OptimizeInner = OptimizeInner{})
        {
          auto expr = associateLeftExpression(std::forward<F>(f), std::forward<T0>(t0), std::forward<T1>(t1), OptimizeInner{});
 
          return operate(
            OptimizeOuter{},
            std::move(expr).operation(),
            std::move(expr).operand(0_c),
            std::move(expr).operand(1_c));
        }
#else
        template<class F, class T0, class T1,
                 class OptimizeOuter = OptimizeTop, class OptimizeInner = OptimizeTop,
                 std::enable_if_t<(// && AreProperTensors<T0, T1>::value
                                   IsLeftAssociative<F, T0, T1>::value), int> = 0>
        constexpr decltype(auto) associateLeft(F&& f, T0&& t0, T1&& t1,
                                               OptimizeOuter = OptimizeOuter{}, OptimizeInner = OptimizeInner{})
        {
          // A*(B*C) -> (A*B)*C
          auto operations = associateLeftOperations(
            std::forward<F>(f), std::forward<T1>(t1).operation(),
            IndexSequence<TensorTraits<T0>::rank, TensorTraits<Operand<0, T1> >::rank, TensorTraits<Operand<1, T1> >::rank>{});
 
          return operate(
            OptimizeOuter{},
            std::move(operations.second), // outer
            operate(
              OptimizeInner{},
              std::move(operations.first), // inner
              std::forward<T0>(t0),
              std::forward<T1>(t1).operand(0_c)
              ),
            std::forward<T1>(t1).operand(1_c)
            );
        }
 
        template<class F, class T0, class T1,
                 class OptimizeOuter = OptimizeTop,
                 std::enable_if_t<(// && AreProperTensors<T0, T1>::value
                                   !IsLeftAssociative<F, T0, T1>::value
                                   && IsProductExpression<T0>::value
                                   && FunctorHas<IsProductOperation, F>::value
          ), int> = 0>
        constexpr decltype(auto) associateLeft(F&& f, T0&& t0, T1&& t1,
                                     OptimizeOuter = OptimizeOuter{}, DontOptimize = DontOptimize{})
        {
          return operate(OptimizeOuter{}, std::forward<F>(f), std::forward<T0>(t0), std::forward<T1>(t1));
        }
#endif
 
        template<std::size_t N, class F, class T0, class T1>
        using AssociateRightOperation =
          TupleElement<N,
                       decltype(
                         associateRightOperations(
                           std::declval<Functor<T0> >(), std::declval<F>(),
                           IndexSequence<TensorTraits<Operand<0, T0> >::rank,
                                         TensorTraits<Operand<1, T0> >::rank,
                                         TensorTraits<T1>::rank>{})
                         )>;
 
        template<std::size_t N, class F, class T0, class T1>
        using AssociateLeftOperation =
          TupleElement<N,
                       decltype(
                         associateLeftOperations(
                           std::declval<F>(), std::declval<Functor<T1> >(),
                           IndexSequence<TensorTraits<T0>::rank,
                                         TensorTraits<Operand<0, T1> >::rank,
                                         TensorTraits<Operand<1, T1> >::rank>{})
                         )>;
 
        template<class T>
        constexpr decltype(auto) rightMostFactor(T&& t)
        {
          if constexpr (!IsProductExpression<T>::value) {
            return std::forward<T>(t);
          } else if constexpr (!IsLeftAssociative<Functor<T>, Operand<0, T>, Operand<1, T> >::value) {
            return std::forward<T>(t).operand(1_c);
          } else {
            return rightMostFactor(std::forward<T>(t).operand(1_c));
          }
        }
 
        template<class T>
        constexpr decltype(auto) leftMostFactor(T&& t)
        {
          if constexpr (!IsProductExpression<T>::value) {
            return std::forward<T>(t);
          } else if constexpr (!IsRightAssociative<Functor<T>, Operand<0, T>, Operand<1, T> >::value) {
            return std::forward<T>(t).operand(0_c);
          } else {
            return leftMostFactor(std::forward<T>(t).operand(0_c));
          }
        }
 
        template<class F, class T0, class T1, class Optimize = OptimizeTop,
                 std::enable_if_t<FunctorHas<IsProductOperation, F>::value, int> = 0>
        constexpr auto factorOutRight(const F& f, T0&& t0, T1&& t1, Optimize = Optimize{})
        {
          if constexpr (!IsLeftAssociative<F, T0, T1>::value) {
            return Expressions::Storage<F, T0, T1>(f, std::forward<T0>(t0), std::forward<T1>(t1));
          } else {
            // A*(B*C) -> (A*B)*C
            auto operations = associateLeftOperations(
              f, std::forward<T1>(t1).operation(),
              IndexSequence<TensorTraits<T0>::rank, TensorTraits<Operand<0, T1> >::rank, TensorTraits<Operand<1, T1> >::rank>{});
 
            return factorOutRight(std::move(operations.second),
                                  operate(
                                    Optimize{},
                                    std::move(operations.first), // inner
                                    std::forward<T0>(t0),
                                    std::forward<T1>(t1).operand(0_c)
                                    ),
                                  std::forward<T1>(t1).operand(1_c));
          }
        }
 
        template<class F, class T0, class T1, class Optimize = OptimizeTop,
                 std::enable_if_t<FunctorHas<IsProductOperation, F>::value, int> = 0>
        constexpr auto factorOutLeft(const F& f, T0&& t0, T1&& t1, Optimize = Optimize{})
        {
          if constexpr (!IsRightAssociative<F, T0, T1>::value) {
            return Expressions::Storage<F, T0, T1>(f, std::forward<T0>(t0), std::forward<T1>(t1));
          } else {
            // (A*B)*C -> A*(B*C)
            auto operations = associateRightOperations(
              std::forward<T0>(t0).operation(), f,
              IndexSequence<TensorTraits<Operand<0, T0> >::rank, TensorTraits<Operand<1, T0> >::rank, TensorTraits<T1>::rank>{});
 
            return factorOutLeft(std::move(operations.first),
                                 std::forward<T0>(t0).operand(0_c),
                                 operate(
                                   Optimize{},
                                   std::move(operations.second), // inner
                                   std::forward<T0>(t0).operand(1_c),
                                   std::forward<T1>(t1)));
          }
        }
 
        template<class F, class T0, class T1>
        using FactorOutRight = std::decay_t<decltype(factorOutRight(std::declval<F>(), std::declval<T0>(), std::declval<T1>()))>;
 
        template<class F, class T0, class T1>
        using FactorOutLeft = std::decay_t<decltype(factorOutLeft(std::declval<F>(), std::declval<T0>(), std::declval<T1>()))>;
 
#if 1
        template<class T>
        using RightMostFactor = decltype(rightMostFactor(std::declval<T>()));
 
        template<class T>
        using LeftMostFactor = decltype(leftMostFactor(std::declval<T>()));
#else
        namespace {
          template<class T, class SFINAE = void>
          struct RightMostFactorHelper
          {
            using Type = T;
          };
 
          template<class T>
          struct RightMostFactorHelper<T, std::enable_if_t<IsProductExpression<T>::value> >
          {
            using Type = ConditionalType<IsLeftAssociative<Functor<T>, Operand<0, T>, Operand<1, T> >::value,
                                         typename RightMostFactorHelper<Operand<1, T> >::Type,
                                         Operand<1, T> >;
          };
 
          template<class T, class SFINAE = void>
          struct LeftMostFactorHelper
          {
            using Type = T;
          };
 
          template<class T>
          struct LeftMostFactorHelper<T, std::enable_if_t<IsProductExpression<T>::value> >
          {
            using Type = ConditionalType<IsRightAssociative<Functor<T>, Operand<0, T>, Operand<1, T> >::value,
                                         typename LeftMostFactorHelper<Operand<0, T> >::Type,
                                         Operand<0, T> >;
          };
        }
 
        template<class T>
        using RightMostFactor = typename RightMostFactorHelper<T>::Type;
 
        template<class T>
        using LeftMostFactor = typename LeftMostFactorHelper<T>::Type;
#endif
 
      } // NS ProductOperations
 
      using ProductOperations::IsRightAssociative;
      using ProductOperations::IsLeftAssociative;
      using ProductOperations::AssociateRightOperation;
      using ProductOperations::associateRightExpression;
      using ProductOperations::associateRight;
      using ProductOperations::associateRight;
      using ProductOperations::AssociateLeftOperation;
      using ProductOperations::associateLeftOperations;
      using ProductOperations::associateLeftExpression ;
      using ProductOperations::associateLeft;
      using ProductOperations::factorOutRight;
      using ProductOperations::rightMostFactor;
      using ProductOperations::RightMostFactor;
      using ProductOperations::FactorOutRight;
      using ProductOperations::factorOutLeft;
      using ProductOperations::leftMostFactor;
      using ProductOperations::LeftMostFactor;
      using ProductOperations::FactorOutLeft;
 
    } // NS Tensor
 
  } // NS ACFem
 
} // NS Dune
 
#endif // __DUNE_ACFEM_TENSORS_OPERATIONS_ASSOCIATIVITY_HH__