im com

src/base_def.h

@@ -114,7 +114,7 @@ namespace MultiArrayTools
     class MultiArrayOperation;
 
     // multi_array_operation.h
-    template <typename T, class Operation, class Range, class... MAOps>
+    template <typename T, class ContractOperation, class Range, class... MAOps>
     class MultiArrayContraction;
     
     // slice.h

src/multi_array.h

@@ -174,6 +174,12 @@ namespace MultiArrayTools
 
 	virtual bool isConst() const override;
 
+	template <class... NameTypes>
+	ConstMultiArrayOperationRoot<T,Range> operator()(bool x, const NameTypes&... str) const
+	{
+	    return MAB::operator()(str...);
+	}
+	
 	template <class... NameTypes>
 	MultiArrayOperationRoot<T,Range> operator()(const NameTypes&... str);
 

src/multi_array_operation.cc

@@ -147,17 +147,21 @@ namespace MultiArrayTools
     MultiArrayOperationRoot<T,Range>&
     MultiArrayOperationRoot<T,Range>::operator=(const MultiArrayOperation<T,Operation,MAOps...>& in)
     {
-	//CHECK;
-	//if(mArrayRef.isSlice() and not mArrayRef.isInit()){
-	    // NO SLICE CREATION !!! (total array not initialized!!)
-	    // throw !
-	//    assert(0);
-	//}
    	performAssignment(in);
 	freeIndex();
 	return *this;
     }
-
+    
+    template <typename T, class Range>
+    template <class Operation, class Range2, class... MAOps>
+    MultiArrayOperationRoot<T,Range>&
+    MultiArrayOperationRoot<T,Range>::operator=(const MultiArrayContraction<T,Operation,Range2,MAOps...>& in)
+    {
+   	performAssignment(in);
+	freeIndex();
+	return *this;
+    }
+    
     template <typename T, class Range>
     template <class Operation, class... MAOps>
     MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>
@@ -293,6 +297,12 @@ namespace MultiArrayTools
 	return 1;
     }
 
+    template <typename T, class Range>
+    IndefinitIndexBase* MultiArrayOperationRoot<T,Range>::getLinked(const std::string& name) const
+    {
+	return mIndex.getLinked(name);
+    }
+    
     template <typename T, class Range>
     void MultiArrayOperationRoot<T,Range>::linkIndicesTo(IndefinitIndexBase* target) const
     {
@@ -534,6 +544,12 @@ namespace MultiArrayTools
 	return 1;
     }
 
+    template <typename T, class Range>
+    IndefinitIndexBase* ConstMultiArrayOperationRoot<T,Range>::getLinked(const std::string& name) const
+    {
+	return mIndex.getLinked(name);
+    }
+    
     template <typename T, class Range>
     void ConstMultiArrayOperationRoot<T,Range>::linkIndicesTo(IndefinitIndexBase* target) const
     {
@@ -630,6 +646,37 @@ namespace MultiArrayTools
 	}
     };
 
+    template <size_t N>
+    struct LinkedIndexGetter
+    {
+	template <class Tuple>
+	static IndefinitIndexBase* getLinked(const Tuple& optuple,
+					     const std::string& name,
+					     IndefinitIndexBase* current)
+	{
+	    if(current == nullptr){
+		current = std::get<N>(optuple).getLinked(name);
+		LinkedIndexGetter<N-1>::getLinked(optuple, name, current);
+	    }
+	    return current;
+	}
+    };
+
+    template <>
+    struct LinkedIndexGetter<0>
+    {
+	template <class Tuple>
+	static IndefinitIndexBase* getLinked(const Tuple& optuple,
+					     const std::string& name,
+					     IndefinitIndexBase* current)
+	{
+	    if(current == nullptr){
+		current = std::get<0>(optuple).getLinked(name);
+	    }
+	    return current;
+	}
+    };
+    
     template <typename T, class Operation, class... MAOps>
     MultiArrayOperation<T,Operation,MAOps...>::
     MultiArrayOperation(Operation& op, const MAOps&... args) :
@@ -662,37 +709,17 @@ namespace MultiArrayTools
 
     
     template <typename T, class Operation, class... MAOps>
-    template <class Range2, class ContractOperation>
+    template <class Range2, class ContractOperation, class... MAOps2>
-    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...> >
+    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...>,MAOps2...>
     MultiArrayOperation<T,Operation,MAOps...>::
-    contract(const ContractOperation& cop, const std::string& indexName) const
+    contract(const ContractOperation& cop,
+	     const std::string& indexName,
+	     const MAOps2&... mao) const
     {
-#error "HERE"
+	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( getLinked(indexName) );
-	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
+	return MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...>,
-	//typename Range2::IndexType ind = Range2().begin();
+				     MAOps2...>(cop, *ind, *this, mao...);
-	//ind.name(indexName);
-	return MultiArrayContraction<T,ContractOperation,Range2,
-				     MultiArrayOperation<T,Operation,MAOps...> >(cop, *this, *ind);
-	
     }
-
-    /*    
-    template <typename T, class Operation, class... MAOps>
-    template <class Range2, class ContractOperation>
-    MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...> >
-    MultiArrayOperation<T,Operation,MAOps...>::
-    contract(const ContractOperation& cop, const std::string& indexName,
-	     const typename Range2::IndexType& begin,
-	     const typename Range2::IndexType& end) const
-    {
-	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( mIndex.getLinked(indexName) );
-	//typename Range2::IndexType ind = Range2().begin();
-	//ind.name(indexName);
-	return MultiArrayContraction<T,ContractOperation,Range2,
-				     MultiArrayOperation<T,Operation,MAOps...> >(cop, *this, *ind, begin, end);
-
-    }
-    */
     
     template <typename T, class Operation, class... MAOps>
     template <class MAOp2>
@@ -732,6 +759,12 @@ namespace MultiArrayTools
 	return sizeof...(MAOps) + 1;
     }
 
+    template <typename T, class Operation, class... MAOps>
+    IndefinitIndexBase* MultiArrayOperation<T,Operation,MAOps...>::getLinked(const std::string& name) const
+    {
+	return LinkedIndexGetter<sizeof...(MAOps)-1>::getLinked(mArgs, name, nullptr);
+    }
+    
     template <typename T, class Operation, class... MAOps>    
     void MultiArrayOperation<T,Operation,MAOps...>::linkIndicesTo(IndefinitIndexBase* target) const
     {
@@ -744,19 +777,131 @@ namespace MultiArrayTools
     {
 	mVal = OperationCall<sizeof...(MAOps)-1>::
 	    template callOperation(mOp, mArgs);
+	std::cout << mVal << std::endl;
 	return mVal;
     }
 
+       
     /*******************************
      *   MultiArrayContraction     *
      *******************************/
 
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    MultiArrayContraction<T,ContractOperation,Range,MAOps...>::
+    MultiArrayContraction(const ContractOperation& op,
+			  const typename Range::IndexType& runIndex,
+			  const MAOps&... args) :
+	mOp(op),
+	mArgs(std::make_tuple(args...)) {}
+
+    // !!!!!
+    
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    template <class Operation2, class... MAOps2>
+    MultiArrayOperation<T,Operation2,MultiArrayContraction<T,ContractOperation,Range,MAOps...>,MAOps2...>
+    MultiArrayContraction<T,ContractOperation,Range,MAOps...>::operator()(Operation2& op, const MAOps2&... secs) const
+    {
+	return MultiArrayOperation<T,Operation2,MultiArrayContraction<T,ContractOperation,Range,MAOps...>,
+				     MAOps2...>(op, *this, secs...);
+    }
+
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    template <class Operation2, class... MAOps2>
+    MultiArrayOperation<T,Operation2,MultiArrayContraction<T,ContractOperation,Range,MAOps...>,MAOps2...>
+    MultiArrayContraction<T,ContractOperation,Range,MAOps...>::operator()(const Operation2& op, const MAOps2&... secs) const
+    {
+	return MultiArrayOperation<T,Operation2,MultiArrayContraction<T,ContractOperation,Range,MAOps...>,
+				     MAOps2...>(op, *this, secs...);
+    }
+
+    
+    template <typename T, class ContractOperation2, class Range, class... MAOps>
+    template <class Range2, class ContractOperation, class... MAOps2>
+    MultiArrayContraction<T,ContractOperation2,Range2,
+			  MultiArrayContraction<T,ContractOperation,Range,MAOps...>,MAOps2...>
+    MultiArrayContraction<T,ContractOperation,Range,MAOps...>::
+    contract(const ContractOperation2& cop,
+	     const std::string& indexName,
+	     const MAOps2&... mao) const
+    {
+	typename Range2::IndexType* ind = dynamic_cast<typename Range2::IndexType*>( getLinked(indexName) );
+	return MultiArrayContraction<T,ContractOperation2,Range2,MultiArrayContraction<T,ContractOperation,Range,MAOps...>,
+				     MAOps2...>(cop, *ind, *this, mao...);
+    }
+    
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    template <class MAOp2>
+    auto MultiArrayContraction<T,ContractOperation,Range,MAOps...>::operator+(const MAOp2& sec)
+	-> decltype(operator()(std::plus<T>(), sec))
+    {
+	return operator()(std::plus<T>(), sec);
+    }
+
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    template <class MAOp2>
+    auto MultiArrayContraction<T,ContractOperation,Range,MAOps...>::operator-(const MAOp2& sec)
+	-> decltype(operator()(std::minus<T>(), sec))
+    {
+	return operator()(std::minus<T>(), sec);
+    }
+
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    template <class MAOp2>
+    auto MultiArrayContraction<T,ContractOperation,Range,MAOps...>::operator*(const MAOp2& sec)
+	-> decltype(operator()(std::multiplies<T>(), sec))
+    {
+	return operator()(std::multiplies<T>(), sec);
+    }
+
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    template <class MAOp2>
+    auto MultiArrayContraction<T,ContractOperation,Range,MAOps...>::operator/(const MAOp2& sec)
+	-> decltype(operator()(std::divides<T>(), sec))
+    {
+	return operator()(std::divides<T>(), sec);
+    }
+    
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    size_t MultiArrayContraction<T,ContractOperation,Range,MAOps...>::argNum() const
+    {
+	return sizeof...(MAOps) + 1;
+    }
+
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    IndefinitIndexBase* MultiArrayContraction<T,ContractOperation,Range,MAOps...>::getLinked(const std::string& name) const
+    {
+	return LinkedIndexGetter<sizeof...(MAOps)-1>::getLinked(mArgs, name, nullptr);
+    }
+    
+    template <typename T, class ContractOperation, class Range, class... MAOps>    
+    void MultiArrayContraction<T,ContractOperation,Range,MAOps...>::linkIndicesTo(IndefinitIndexBase* target) const
+    {
+	TupleIndicesLinker<sizeof...(MAOps)-1>::linkTupleIndicesTo(mArgs, target);
+    }
+
+
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    const T& MultiArrayContraction<T,ContractOperation,Range,MAOps...>::get() const
+    {
+	mOp.reset();
+	for(mRunIndex.copyPos( mBeginIndex ); mRunIndex.pos() != mEndIndex.pos(); ++mRunIndex){
+	    OperationCall<sizeof...(MAOps)-1>::
+		template callOperation(mOp, mArgs);
+	    //MAO::mOp(std::get<0>(MAO::mArgs).get() );
+	}
+	mOp.endOp();
+	std::cout << MAO::mOp() << std::endl;
+	return mOp();
+    }
+
+    
+    /*
     template <typename T, class ContractOperation, class Range, class... MAOps>
     MultiArrayContraction<T,ContractOperation,Range,MAOps...>::
     MultiArrayContraction(const ContractOperation& cop,
 			  const typename Range::IndexType& runIndex,
 			  const MAOps&... mao) :
-	MultiArrayOperation<T,ContractOperation,MAOps...>(cop, mao...),
+	MultiArrayContraction<T,ContractOperation,MAOps...>(cop, mao...),
 	mBeginIndex(runIndex), mEndIndex(runIndex),
 	mRunIndex(runIndex)
     {
@@ -772,7 +917,7 @@ namespace MultiArrayTools
 			  size_t begin,
 			  size_t end,
 			  const MAOps&... mao) :
-	MultiArrayOperation<T,ContractOperation,MAOps...>(cop, mao...),
+	MultiArrayContraction<T,ContractOperation,MAOps...>(cop, mao...),
 	mBeginIndex(runIndex), mEndIndex(runIndex),
 	mRunIndex(runIndex)
     {
@@ -793,6 +938,7 @@ namespace MultiArrayTools
 	    //MAO::mOp(std::get<0>(MAO::mArgs).get() );
 	}
 	MAO::mOp.endOp();
+	std::cout << MAO::mOp() << std::endl;
 	return MAO::mOp();
-    }
+	}*/
 }

src/multi_array_operation.h

@@ -22,6 +22,7 @@ namespace MultiArrayTools
 
 	virtual size_t argNum() const = 0;
 	const IndefinitIndexBase& index() const;
+	virtual IndefinitIndexBase* getLinked(const std::string& name) const = 0;
 	virtual void linkIndicesTo(IndefinitIndexBase* target) const = 0;
 	
 	virtual const T& get() const = 0;
@@ -65,7 +66,11 @@ namespace MultiArrayTools
 	
 	template <class Operation, class... MAOps>
 	MultiArrayOperationRoot& operator=(const MultiArrayOperation<T,Operation,MAOps...>& in);
-	
+
+	template <class Operation, class Range2, class... MAOps>
+	MultiArrayOperationRoot<T,Range>&
+	operator=(const MultiArrayContraction<T,Operation,Range2,MAOps...>& in);
+
 	//template <class Operation, class... MAOps>
 	//MultiArrayOperation<T,Operation,MultiArrayOperationRoot<T,Range>, MAOps...>
 	//operator()(Operation& op, const MAOps&... secs) const;
@@ -122,7 +127,8 @@ namespace MultiArrayTools
 	// set index -> implement !!!!!
 	MultiArrayOperationRoot<T,Range>& operator[](const IndexType& ind);
 	const MultiArrayOperationRoot<T,Range>& operator[](const IndexType& ind) const;
-	
+
+	virtual IndefinitIndexBase* getLinked(const std::string& name) const override;
 	virtual void linkIndicesTo(IndefinitIndexBase* target) const override;
 
 	virtual T& get() override;
@@ -207,7 +213,8 @@ namespace MultiArrayTools
 
 	// set index -> implement !!!!!
 	const ConstMultiArrayOperationRoot<T,Range>& operator[](const IndexType& ind) const;
-	
+
+	virtual IndefinitIndexBase* getLinked(const std::string& name) const override;
 	virtual void linkIndicesTo(IndefinitIndexBase* target) const override;
 
 	virtual const T& get() const override;
@@ -259,17 +266,11 @@ namespace MultiArrayTools
 	operator()(const Operation2& op, const MAOps2&... secs) const;
 
 	
-	template <class Range2, class ContractOperation>
+	template <class Range2, class ContractOperation, class... MAOps2>
-	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...> >
+	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...>,MAOps2...>
-	contract(const ContractOperation& cop, const std::string& indexName) const;
-
-	/*
-	template <class Range2, class ContractOperation>
-	MultiArrayContraction<T,ContractOperation,Range2,MultiArrayOperation<T,Operation,MAOps...> >
 	contract(const ContractOperation& cop, const std::string& indexName,
-		 const typename Range2::IndexType& begin,
+		 const MAOps2&... mao) const;
-		 const typename Range2::IndexType& end) const;
+
-	*/
 	
 	template <class MAOp2>
 	auto operator+(const MAOp2& sec) -> decltype(operator()(std::plus<T>(), sec));
@@ -284,7 +285,8 @@ namespace MultiArrayTools
 	auto operator/(const MAOp2& sec) -> decltype(operator()(std::divides<T>(), sec));
 	
 	virtual size_t argNum() const override;
-	
+
+	virtual IndefinitIndexBase* getLinked(const std::string& name) const override;
 	virtual void linkIndicesTo(IndefinitIndexBase* target) const override;
 
 	virtual const T& get() const override;
@@ -296,6 +298,64 @@ namespace MultiArrayTools
 	OBT mArgs; // include first arg also here !!!
     };
 
+    template <typename T, class ContractOperation, class Range, class... MAOps>
+    class MultiArrayContraction : public MultiArrayOperationBase<T>
+    {
+    public:
+
+	typedef MultiArrayOperationBase<T> MAOB;
+	typedef std::tuple<MAOps...> OBT;
+	
+	MultiArrayContraction(ContractOperation& op, const MAOps&... secs);
+	MultiArrayContraction(const ContractOperation& op, const MAOps&... secs);
+
+	template <class Operation2, class... MAOps2>
+	MultiArrayOperation<T,Operation2,MultiArrayContraction<T,ContractOperation,Range,MAOps...>,MAOps2...>
+	operator()(Operation2& op, const MAOps2&... secs) const;
+	
+	template <class Operation2, class... MAOps2>
+	MultiArrayOperation<T,Operation2,MultiArrayContraction<T,ContractOperation,Range,MAOps...>,MAOps2...>
+	operator()(const Operation2& op, const MAOps2&... secs) const;
+
+	
+	template <class Range2, class ContractOperation2, class... MAOps2>
+	MultiArrayContraction<T,ContractOperation2,Range2,
+			      MultiArrayContraction<T,ContractOperation,Range,MAOps...>,MAOps2...>
+	contract(const ContractOperation2& cop, const std::string& indexName,
+		 const MAOps2&... mao) const;
+
+	
+	template <class MAOp2>
+	auto operator+(const MAOp2& sec) -> decltype(operator()(std::plus<T>(), sec));
+
+	template <class MAOp2>
+	auto operator-(const MAOp2& sec) -> decltype(operator()(std::minus<T>(), sec));
+
+	template <class MAOp2>
+	auto operator*(const MAOp2& sec) -> decltype(operator()(std::multiplies<T>(), sec));
+
+	template <class MAOp2>
+	auto operator/(const MAOp2& sec) -> decltype(operator()(std::divides<T>(), sec));
+	
+	virtual size_t argNum() const override;
+
+	virtual IndefinitIndexBase* getLinked(const std::string& name) const override;
+	virtual void linkIndicesTo(IndefinitIndexBase* target) const override;
+
+	virtual const T& get() const override;
+
+    protected:
+
+ 	mutable T mVal;
+	ContractOperation mOp;
+	OBT mArgs; // include first arg also here !!!
+	typename Range::IndexType mBeginIndex;
+	typename Range::IndexType mEndIndex;
+        mutable typename Range::IndexType mRunIndex;
+
+    };
+
+    /*
     template <typename T, class ContractOperation, class Range, class... MAOps>
     class MultiArrayContraction : public MultiArrayOperation<T,ContractOperation,MAOps...>
     {
@@ -312,7 +372,7 @@ namespace MultiArrayTools
 			      size_t begin,
 			      size_t end,
 			      const MAOps&... mao);
-	
+
 	virtual const T& get() const override;
 
     protected:
@@ -320,7 +380,7 @@ namespace MultiArrayTools
 	typename Range::IndexType mEndIndex;
         mutable typename Range::IndexType mRunIndex;
     };
-    
+    */
 }
 
 #include "multi_array_operation.cc"

src/unit_test.cc

@@ -178,7 +178,8 @@ namespace {
 			    rb(r1,r2),
 			    r3d(r1,r2,r3),
 			    ma(r3d, {-5,6,2,1,9,54,27,-7,-13,32,90,-67,
-					-10,16,-2,101,39,-64,81,-22,14,34,95,-62}) {}
+					-10,16,-2,101,39,-64,81,-22,14,34,95,-62}),
+			    max(ra, {-5,6,2,1,9,54}){}
 	
 	Range1dAny r1;
 	Range1dAny r2;
@@ -187,8 +188,10 @@ namespace {
 	Range2dAny rb;
 	Range3dAny r3d;
 	MultiArray3dAny ma;
+	MultiArray2dAny max;
     };
-    
+
+    /*
     TEST_F(OneDimTest, CorrectExtensions)
     {
 	EXPECT_EQ(ma.size(), 5);
@@ -507,7 +510,7 @@ namespace {
 	EXPECT_EQ(sl[j(j1 = 2, j2 = 0)], 14);
 	EXPECT_EQ(sl[j(j1 = 2, j2 = 1)], 34);
     }
-
+    */
     TEST_F(ContractionTest, ContractionWorks)
     {
 	MultiArray2dAny ma2(ra);
@@ -528,6 +531,27 @@ namespace {
 	EXPECT_EQ(ma2[i(i1 = 2, i2 = 1)], -114);
     }
     
+    TEST_F(ContractionTest, ContractionWorks_2)
+    {
+	MultiArray2dAny ma2(ra);
+
+	ma2("alpha","gamma") = ma("alpha","beta","gamma").contract<Range1dAny>(sum<int>(),"beta")
+	    * max("alpha","gamma");
+
+	auto i = ma2.beginIndex();
+	auto i1 = i.template getIndex<0>();
+	auto i2 = i.template getIndex<1>();
+
+	EXPECT_EQ(ma2[i(i1 = 0, i2 = 0)], -275);
+	EXPECT_EQ(ma2[i(i1 = 0, i2 = 1)], 324);
+
+	EXPECT_EQ(ma2[i(i1 = 1, i2 = 0)], 130);
+	EXPECT_EQ(ma2[i(i1 = 1, i2 = 1)], 82);
+
+	EXPECT_EQ(ma2[i(i1 = 2, i2 = 0)], 2061);
+	EXPECT_EQ(ma2[i(i1 = 2, i2 = 1)], -6156);
+    }
+    
 } // end namespace 
 
 int main(int argc, char** argv)