first part of high level operations work (only implemented in context of op4utest)

src/include/dynamic_operation.h

@@ -127,11 +127,15 @@ namespace MultiArrayTools
 	static constexpr bool CONT = true;
 	
 	DynamicO() = default;
-	DynamicO(const DynamicO& in) : mOp(in.mOp->deepCopy()) {}
+	DynamicO(const DynamicO& in) : mOp(in.mOp ? in.mOp->deepCopy() : nullptr) {}
-	DynamicO(DynamicO&& in) : mOp(in.mOp->deepCopy()) {}
+	DynamicO(DynamicO&& in) : mOp(in.mOp ? in.mOp->deepCopy() : nullptr) {}
-	DynamicO& operator=(const DynamicO& in) { mOp = in.mOp->deepCopy(); }
+	DynamicO& operator=(const DynamicO& in)
-	DynamicO& operator=(DynamicO&& in) { mOp = in.mOp->deepCopy(); }
+	{ mOp = in.mOp ? in.mOp->deepCopy() : nullptr; return *this; }
+	DynamicO& operator=(DynamicO&& in)
+	{ mOp = in.mOp ? in.mOp->deepCopy() : nullptr; return *this; }
 
+	bool init() const { return mOp != nullptr; }
+	
 	template <class Op>
 	DynamicO(const Op& op) : mOp(std::make_shared<DynamicOperation<T,Op>>(op)) {}
 

src/include/helper_tools.cc.h

@@ -374,5 +374,11 @@ namespace MultiArrayTools
         return MultiArray<T,Ranges...>(rs..., val);
     }
 
+    template <typename T, class... Ranges>
+    auto mkArrayPtr(const std::shared_ptr<Ranges>&... rs)
+        -> std::shared_ptr<MultiArray<T,Ranges...>>
+    {
+	return std::make_shared<MultiArray<T,Ranges...>>(rs...);
+    }
 
 }

src/include/helper_tools.h

@@ -155,6 +155,10 @@ namespace MultiArrayTools
     auto mkArray(const std::shared_ptr<Ranges>&... rs)
         -> MultiArray<T,Ranges...>;
 
+    template <typename T, class... Ranges>
+    auto mkArrayPtr(const std::shared_ptr<Ranges>&... rs)
+        -> std::shared_ptr<MultiArray<T,Ranges...>>;
+
     template <typename T, class... Ranges>
     auto mkArray(const std::shared_ptr<Ranges>&... rs, const T& val)
         -> MultiArray<T,Ranges...>;

src/include/multi_array_operation.cc.h

@@ -211,6 +211,39 @@ namespace MultiArrayTools
         return nullptr; //???!!!
     }
 
+    template <typename T, class... Ops>
+    MOp<T,Ops...>::MOp(const Ops&... exprs) : mOps(exprs...)
+    {
+	static_assert(SIZE == sizeof...(Ops), "size missmatch");
+    }
+
+    template <typename T, class... Ops>
+    inline size_t MOp<T,Ops...>::get(ExtType last) const
+    {
+	return RootSumN<sizeof...(Ops)-1>::get(last,mOps);
+    }
+    
+    template <typename T, class... Ops>
+    inline MOp<T,Ops...>& MOp<T,Ops...>::set(ExtType last)
+    {
+	RootSumN<sizeof...(Ops)-1>::set(last,mOps);
+	return *this;
+    }
+    
+    template <typename T, class... Ops>
+    template <class Expr>
+    auto MOp<T,Ops...>::loop(Expr exp) const
+	-> decltype(PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp))
+    {
+	return PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp);
+    }
+
+    template <typename T, class... Ops>
+    auto MOp<T,Ops...>::rootSteps(std::intptr_t iPtrNum) const -> ExtType
+    {
+	return RootSumN<sizeof...(Ops)-1>::rootSteps(mOps,iPtrNum);
+    }
+
     template <class OpClass, class NextExpr>
     GetExpr<OpClass,NextExpr>::GetExpr(const OpClass& sec, const NextExpr& nexpr) :
         mSec(sec), mNExpr(nexpr) {}
@@ -230,7 +263,13 @@ namespace MultiArrayTools
         mSec.get(last);
         mNExpr(start,last.next());
     }
-    
+
+    template <class OpClass, class NextExpr>
+    inline void GetExpr<OpClass,NextExpr>::get(ExtType last)
+    {
+	(*this)(0,last);
+    }
+
     template <class OpClass, class NextExpr>
     typename GetExpr<OpClass,NextExpr>::ExtType
     GetExpr<OpClass,NextExpr>::rootSteps(std::intptr_t iPtrNum) const

src/include/multi_array_operation.h

@@ -98,6 +98,98 @@ namespace MultiArrayTools
 	friend OperationClass;
     };
 
+    template <class Op>
+    size_t sumRootNum()
+    {
+	return typename Op::rootNum();
+    }
+    
+    template <class Op1, class Op2, class... Ops>
+    size_t sumRootNum()
+    {
+	return typename Op1::rootNum() + sumRootNum<Op2,Ops...>();
+    }
+
+    template <size_t N>
+    struct RootSumN
+    {
+	template <class Op1, class... Ops>
+	struct rs
+	{
+	    static constexpr size_t SIZE = Op1::SIZE + RootSumN<N-1>::template rs<Ops...>::SIZE;
+	};
+
+	template <class... Exprs>
+	static inline auto rootSteps(const std::tuple<Exprs...>& etp, std::intptr_t i)
+	{
+	    return RootSumN<N-1>::rootSteps(etp,i).extend( std::get<N>(etp).rootSteps(i) );
+	}
+
+	template <class ExtType, class... Exprs>
+	static inline void exec( size_t start, ExtType last, std::tuple<Exprs...>& etp)
+	{
+	    std::get<sizeof...(Exprs)-N-1>(etp)(start,last);
+	    RootSumN<N-1>::exec(start,last.next(),etp);
+	}
+
+	template <class ExtType, class... Exprs>
+	static inline size_t get( ExtType last, const std::tuple<Exprs...>& etp)
+	{
+	    std::get<sizeof...(Exprs)-N-1>(etp).get(last);
+	    return RootSumN<N-1>::get(last.next(),etp);
+	}
+
+	template <class ExtType, class... Exprs>
+	static inline void set( ExtType last, std::tuple<Exprs...>& etp)
+	{
+	    std::get<sizeof...(Exprs)-N-1>(etp).set(last);
+	    RootSumN<N-1>::set(last.next(),etp);
+	}
+    };
+
+    template <>
+    struct RootSumN<0>
+    {
+	template <class Op1>
+	struct rs
+	{
+	    static constexpr size_t SIZE = Op1::SIZE;
+	};
+
+	template <class... Exprs>
+	static inline auto rootSteps(const std::tuple<Exprs...>& etp, std::intptr_t i)
+	{
+	    return std::get<0>(etp).rootSteps(i);
+	}
+
+	template <class ExtType, class... Exprs>
+	static inline void exec( size_t start, ExtType last, std::tuple<Exprs...>& etp)
+	{
+	    std::get<sizeof...(Exprs)-1>(etp)(start,last);
+	}
+
+	template <class ExtType, class... Exprs>
+	static inline size_t get( ExtType last, const std::tuple<Exprs...>& etp)
+	{
+	    std::get<sizeof...(Exprs)-1>(etp).get(last);
+	    return 0;
+	}
+
+	template <class ExtType, class... Exprs>
+	static inline void set( ExtType last, std::tuple<Exprs...>& etp)
+	{
+	    std::get<sizeof...(Exprs)-1>(etp).set(last);
+	}
+    };
+
+    
+    template <class... Ops>
+    struct RootSum
+    {
+	static constexpr size_t SIZE = RootSumN<sizeof...(Ops)-1>::template rs<Ops...>::SIZE;
+    };
+
+    
     template <typename T>
     struct SelfIdentity
     {
@@ -149,6 +241,38 @@ namespace MultiArrayTools
         inline DExt dExtension() const override final;
     };
 
+    template <typename T, class... Ops>
+    class MOp
+    {
+    private:
+	MOp() = default;
+	std::tuple<Ops...> mOps;
+
+    public:
+        static constexpr size_t LAYER = 0;
+        static constexpr size_t SIZE = RootSum<Ops...>::SIZE;
+        typedef decltype(RootSumN<sizeof...(Ops)-1>::rootSteps(mOps,0) ) ExtType;
+
+	MOp(const Ops&... exprs);
+	
+	MOp(const MOp& in) = default;
+	MOp(MOp&& in) = default;
+	MOp& operator=(const MOp& in) = default;
+	MOp& operator=(MOp&& in) = default;
+
+        inline size_t get(ExtType last) const;
+	inline MOp& set(ExtType last);
+        auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
+
+	template <class Expr>
+	auto loop(Expr exp) const
+	    -> decltype(PackNum<sizeof...(Ops)-1>::mkLoop( mOps, exp));
+
+	T* data() const { assert(0); return nullptr; }
+
+    };
+
+    
     template <class OpClass, class NextExpr>
     class GetExpr : public ExpressionBase
     {
@@ -176,6 +300,7 @@ namespace MultiArrayTools
 	}
 
         inline void operator()(size_t start = 0);
+        inline void get(ExtType last);
         inline void operator()(size_t start, ExtType last);
         auto rootSteps(std::intptr_t iPtrNum = 0) const -> ExtType;
 
@@ -190,7 +315,13 @@ namespace MultiArrayTools
     {
         return GetExpr<OpClass,NextExpr>(op, nexpr);
     }
-    
+
+    template <typename T, class... Ops>
+    auto mkMOp(const Ops&... exprs)
+    {
+        return MOp<T,Ops...>(exprs...);
+    }
+
     //template <typename T, class OpClass>
     template <typename T, class Target, class OpClass, OpIndexAff OIA=OpIndexAff::EXTERN>
     class AddExpr : public ExpressionBase
@@ -513,44 +644,6 @@ namespace MultiArrayTools
 	T mVal;
     };
     
-    template <class Op>
-    size_t sumRootNum()
-    {
-	return typename Op::rootNum();
-    }
-    
-    template <class Op1, class Op2, class... Ops>
-    size_t sumRootNum()
-    {
-	return typename Op1::rootNum() + sumRootNum<Op2,Ops...>();
-    }
-
-    template <size_t N>
-    struct RootSumN
-    {
-	template <class Op1, class... Ops>
-	struct rs
-	{
-	    static constexpr size_t SIZE = Op1::SIZE + RootSumN<N-1>::template rs<Ops...>::SIZE;
-	};
-    };
-
-    template <>
-    struct RootSumN<0>
-    {
-	template <class Op1>
-	struct rs
-	{
-	    static constexpr size_t SIZE = Op1::SIZE;
-	};
-    };
-
-    
-    template <class... Ops>
-    struct RootSum
-    {
-	static constexpr size_t SIZE = RootSumN<sizeof...(Ops)-1>::template rs<Ops...>::SIZE;
-    };
 
     
     template <typename T, class OpFunction, class... Ops>

src/tests/op4_unit_test.cc

@@ -299,24 +299,26 @@ namespace
         resx2(i1,di4) = mkDynOp(ma1(i1,di1) * exp(ma2(i1,di2)));
 	resx3(i1,di4) = mkDynOp(mkDynOp(ma1(i1,di1)) * mkDynOp(exp(mkDynOp(ma2(i1,di2)))));
 	
-	auto xx = std::make_shared<decltype(resx4)>(resx4);
+	//auto xx = std::make_shared<decltype(resx4)>(resx4);
+	auto xx = mkArrayPtr<double>(nullr());
 	auto mi = mkMIndex(i1,di4a);
-	
+
-	auto op1 = ma2(i1,di2);
+	auto hop1 = mkHLO(ma2(i1,di2));
-	auto op3 = ma1(i1,di1);
+	auto hop3 = mkHLO(ma1(i1,di1));
+	auto hop2 = exp(hop1);
+	auto hop4 = hop3 * hop2;
+
 	auto opr = resx4(i1,di4);
 	auto loop = mkPILoop
-	    ( [&op1,&op3,&opr,&xx,&ic_1,&ic_2,this](){
+	    ( [&opr,&hop4,&xx,&ic_1,&ic_2,this](){
-		auto dop1 = mkDynOutOp(exp(op1), ic_1, ic_2);
+		auto hop4x = hop4;
-		auto op2 = *dop1.data()->mOp;
+		auto dop2 = hop4x.create(ic_1,ic_2);
-		auto dop2 = mkDynOutOp(op3 * op2, ic_1, ic_2);
+		auto gexp = mkDynOp1<size_t>(mkMOp<size_t>(dop2.outer,dop2.op));
-		return mkGetExpr
+		auto xloop = mkILoop(std::make_tuple(*dop2.op.data()->mOp), std::make_tuple(ic_1, ic_2),
-		    (dop1,mkGetExpr
+				     std::make_tuple(xx),
-		     (dop2,mkILoop(std::make_tuple(*dop2.data()->mOp), std::make_tuple(ic_1, ic_2),
+				     std::make_tuple(opr.assign( *dop2.op.data()->mOp, mkMIndex(ic_1, ic_2) )),
-				   std::make_tuple(xx),
+				     std::array<size_t,1>({1}), std::array<size_t,1>({0}));
-				   std::make_tuple(opr.assign( *dop2.data()->mOp, mkMIndex(ic_1, ic_2) )),
+		return mkGetExpr(gexp, xloop); });
-				   std::array<size_t,1>({1}), std::array<size_t,1>({0})))); } );
-	
         mi->pifor(1,loop)();
         
 	auto i2_1 = imap.at("i2_1");

src/tests/op4_unit_test.h

@@ -7,6 +7,13 @@ typedef OperationRoot<double,DR> OpD;
 typedef CR::IndexType CI;
 typedef std::shared_ptr<CI> CIP;
 
+template <typename T, class Op>
+DynamicO<T> mkDynOp1(const Op& op)
+{
+    return DynamicO<T>(op);
+}
+
+
 class OpHolderBase
 {
 public:
@@ -26,7 +33,7 @@ struct OpName
     {
         return 0;
     }
-}
+};
 
 template <class Op>
 class OpHolder : public OpHolderBase
@@ -53,39 +60,117 @@ public:
 
 class HighLevelOpBase
 {
-protected:
-    std::shared_ptr<HighLevelOpBase> mPrev
-    
 public:
     
     typedef OperationRoot<double,CR,CR> OpCC;
     typedef OperationRoot<double,CR,DR> OpCD;
     typedef OperationRoot<double,DR> OpD;
 
+    struct RetT
+    {
+	DynamicO<OpH<OpCC>> op;
+	DynamicO<size_t> outer;
+
+	template <class... Ops>
+	void appendOuterM(const Ops&... ops)
+	{
+	    // does not check anything regarding input !!!
+	    if(outer.init()){
+		outer = mkDynOp1<size_t>(mkMOp<size_t>(outer,ops...));
+	    }
+	    else {
+		outer = mkDynOp1<size_t>(mkMOp<size_t>(ops...));
+	    }
+	}
+	
+	void appendOuter(const DynamicO<size_t>& in)
+	{
+	    if(in.init()){
+		if(outer.init()){
+		    outer = mkDynOp1<size_t>(mkMOp<size_t>(outer,in));
+		}
+		else {
+		    outer = in;
+		}
+	    }
+	}
+
+	void appendOuter(const RetT& in)
+	{
+	    appendOuter(in.outer);
+	}
+
+    };
+    
     virtual bool root() const = 0;
     
-    virtual DynamicO<OpH<OpRR>>
+    virtual RetT create(const std::shared_ptr<CI> ind1,
-    create(const std::shared_ptr<CI> ind1,
+			const std::shared_ptr<CI> ind2) = 0;
-           const std::shared_ptr<CI> ind2) = 0;
 
-    virtual OperationRoot<double,DR> get1() = 0;
+    virtual const OperationRoot<double,DR>* get1() const = 0;
-    virtual OperationRoot<double,CR,DR> get2() = 0;
+    virtual const OperationRoot<double,CR,DR>* get2() const = 0;
 
-    void appendPrev(const std::shared_ptr<HighLevelOpBase>& in)
+};
+
+template <bool SAME>
+struct Fwd
+{
+    template <class O1, class O2>
+    static inline const O1* fwd(const O2* in)
     {
-        if(mPrev){
+	assert(0);
-            mPrev->appendPrev(in)
+	return nullptr;
-        }
-        else {
-            mPrev = in;
-        }
     }
 };
 
+template <>
+struct Fwd<true>
+{
+    template <class O1, class O2>
+    static inline const O1* fwd(const O2* in)
+    {
+	return in;
+    }
+};
+
+template <class OR>
 class HighLevelOpRoot : public HighLevelOpBase
 {
+private:
+    typedef OperationRoot<double,DR> OType1;
+    typedef OperationRoot<double,CR,DR> OType2;
+    typedef HighLevelOpBase B;
     
-}
+    OR mOp;
+public:
+
+    HighLevelOpRoot(const OR& op) : mOp(op) {}
+    
+    virtual bool root() const override final
+    {
+	return true;
+    }
+    
+    virtual B::RetT create(const std::shared_ptr<CI> ind1,
+			   const std::shared_ptr<CI> ind2) override final
+    {
+	assert(0);
+	return B::RetT();
+    }
+
+    virtual const OType1* get1() const override final
+    {
+	return Fwd<std::is_same<OType1,OR>::value>::template fwd<OType1>(&mOp);
+    }
+    
+    virtual const OType2* get2() const override final
+    {
+	return Fwd<std::is_same<OType2,OR>::value>::template fwd<OType2>(&mOp);
+    }
+
+   
+};
+
 
 template <class OpF>
 class HighLevelOp2 : public HighLevelOpBase
@@ -93,37 +178,136 @@ class HighLevelOp2 : public HighLevelOpBase
 private:
     std::shared_ptr<HighLevelOpBase> mIn1;
     std::shared_ptr<HighLevelOpBase> mIn2;
-
+    
 public:
     typedef HighLevelOpBase B;
     typedef Operation<double,OpF,B::OpCC,B::OpCC> OxCCxCC;
     typedef Operation<double,OpF,B::OpCD,B::OpCD> OxCDxCD;
+    typedef Operation<double,OpF,B::OpCC,B::OpCD> OxCCxCD;
+    typedef Operation<double,OpF,B::OpCD,B::OpCC> OxCDxCC;
 
+    HighLevelOp2(const std::shared_ptr<HighLevelOpBase>& in1,
+		 const std::shared_ptr<HighLevelOpBase>& in2) : mIn1(in1), mIn2(in2) {}
+    
     virtual bool root() const override final
     {
         return false;
     }
+    virtual const OperationRoot<double,DR>* get1() const override final
+    { assert(0); return nullptr; }
+    virtual const OperationRoot<double,CR,DR>* get2() const override final
+    { assert(0); return nullptr; }
     
-    virtual DynamicO<OpH<OperationRoot<double,CR,CR>>>
+    virtual B::RetT create(const std::shared_ptr<CI> ind1,
-    create(const std::shared_ptr<CI> ind1,
+			   const std::shared_ptr<CI> ind2) override final
-           const std::shared_ptr<CI> ind2) override final
     {
-        if(not mIn1->root()){
+	B::RetT res;
-            DynamicO<OpH<OperationRoot<double,CR,CR>>> dop1 = mIn1->create(ind1,ind2);
+	if(not mIn2->root()){
-            op1 = *dop1.data()->mOp;
+	    auto dop2 = mIn2->create(ind1,ind2);
-            if(not mIn2->root()){
+            auto op2 = *dop2.op.data()->mOp;
-                DynamicO<OpH<OperationRoot<double,CR,CR>>> dop2 = mIn2->create(ind1,ind2);
+	    res.outer = dop2.outer;
-                op2 = *dop2.data()->mOp;
+            if(not mIn1->root()){
-                return mkDynOutOp(OxRRRR(op1,op2), ind1, ind2);
+		auto dop1 = mIn1->create(ind1,ind2);
+                auto op1 = *dop1.op.data()->mOp;
+                res.op = mkDynOutOp(OxCCxCC(op1,op2), ind1, ind2);
+		res.appendOuter(dop1);
+		res.appendOuterM(dop1.op,dop2.op);
+		return res;
             }
             else {
-                OperationRoot<double,DR,CR>& op2 = mIn2->get2(ind1,ind2);
+                const OperationRoot<double,CR,DR>& op1 = *mIn1->get2();
-                return mkDynOutOp(OxRRRR(op1,op2), ind1, ind2);
+		res.op = mkDynOutOp(OxCDxCC(op1,op2), ind1, ind2);
+		res.appendOuterM(dop2.op);
+		return res;
             }
         }
-        return DynamicO<OpH<OperationRoot<double,CR,CR>>>();
+	assert(0);
+        return res;
     }
 
 };
 
+template <class OpF>
+class HighLevelOp1 : public HighLevelOpBase
+{
+private:
+    std::shared_ptr<HighLevelOpBase> mIn1;
+public:
+    typedef HighLevelOpBase B;
+    typedef Operation<double,OpF,B::OpCC> OxCC;
+    typedef Operation<double,OpF,B::OpCD> OxCD;
 
+    HighLevelOp1(const std::shared_ptr<HighLevelOpBase>& in1) : mIn1(in1) {}
+    
+    virtual bool root() const override final
+    {
+        return false;
+    }
+    virtual const OperationRoot<double,DR>* get1() const override final
+    { assert(0); return nullptr; }
+    virtual const OperationRoot<double,CR,DR>* get2() const override final
+    { assert(0); return nullptr; }
+    
+    virtual B::RetT create(const std::shared_ptr<CI> ind1,
+			   const std::shared_ptr<CI> ind2) override final
+    {
+	if(not mIn1->root()){
+	    auto dop1 = mIn1->create(ind1,ind2);
+	    auto op1 = *dop1.op.data()->mOp;
+	    dop1.appendOuterM(dop1.op);
+	    dop1.op = mkDynOutOp(OxCC(op1), ind1, ind2);
+	    return dop1;
+	}
+	else {
+	    const OperationRoot<double,CR,DR>& op1 = *mIn1->get2();
+	    B::RetT res;
+	    res.op = mkDynOutOp(OxCD(op1), ind1, ind2); 
+	    return res;
+	}
+    }
+
+
+};
+
+class HighLevelOpHolder
+{
+private:
+    std::shared_ptr<HighLevelOpBase> mOp;
+
+public:
+    HighLevelOpHolder() = default;
+    HighLevelOpHolder(const HighLevelOpHolder& in) = default; 
+    HighLevelOpHolder(HighLevelOpHolder&& in) = default;
+    HighLevelOpHolder& operator=(const HighLevelOpHolder& in) = default;
+    HighLevelOpHolder& operator=(HighLevelOpHolder&& in) = default;
+    
+    HighLevelOpHolder(const std::shared_ptr<HighLevelOpBase>& op) : mOp(op) {}
+
+    bool root() const { return mOp->root(); }
+    
+    auto create(const std::shared_ptr<CI> ind1,
+		const std::shared_ptr<CI> ind2) const { return mOp->create(ind1,ind2); }
+
+    auto get1() const { return mOp->get1(); }
+    auto get2() const { return mOp->get2(); }
+
+    std::shared_ptr<HighLevelOpBase> op() const { return mOp; }
+    
+    HighLevelOpHolder operator*(const HighLevelOpHolder in) const
+    {
+	return HighLevelOpHolder
+	    ( std::make_shared<HighLevelOp2<multipliesx<double,double>>>
+	      ( mOp, in.mOp ) );
+    }
+};
+
+template <class OR>
+HighLevelOpHolder mkHLO(const OR& op)
+{
+    return HighLevelOpHolder(std::make_shared<HighLevelOpRoot<OR>>( op ) );
+}
+
+HighLevelOpHolder exp(const HighLevelOpHolder& in)
+{
+    return HighLevelOpHolder( std::make_shared<HighLevelOp1<x_exp<double>>>( in.op() ) );
+}