contraction compiles (untested)

src/block.cc

@@ -22,6 +22,20 @@ namespace MultiArrayHelper
 	}
 	return res;
     }
+
+    template <typename T, class OpFunc, class BlockClass>
+    BlockBinaryOpSelf<T,OpFunc,BlockClass>::BlockBinaryOpSelf(BlockResult<T>& res) : mRes(res) {}
+
+    template <typename T, class OpFunc, class BlockClass>
+    void BlockBinaryOpSelf<T,OpFunc,BlockClass>::operator()(const BlockClass& arg)
+    {
+	static OpFunc f;
+	assert(mRes.size() == arg.size());
+	for(size_t i = 0; i != arg.size(); ++i){
+	    mRes[i] = f(mRes[i], arg[i]);
+	}
+    }
+
     
     /*****************
      *   BlockBase   *
@@ -35,7 +49,7 @@ namespace MultiArrayHelper
     {
 	return mSize;
     }
-
+    /*
     template <typename T>
     template <class OpFunction>
     BlockResult<T> BlockBase<T>::operate(const BlockBase<T>& in)
@@ -43,6 +57,7 @@ namespace MultiArrayHelper
 	assert(mSize == in.size());
 	OpFunction f;
 	BlockResult<T> res(mSize);
+	CHECK;
 	for(size_t i = 0; i != mSize; ++i){
 	    res[i] = f((*this)[i], in[i]);
 	}
@@ -73,7 +88,7 @@ namespace MultiArrayHelper
     {
 	return operate<std::divides<T> >(in);
     }
-
+    */
     /************************
      *   MutableBlockBase   *
      ************************/
@@ -190,6 +205,7 @@ namespace MultiArrayHelper
     template <class BlockClass>
     BlockResult<T>& BlockResult<T>::operator=(const BlockClass& in)
     {
+	//CHECK;
 	for(size_t i = 0; i != BlockBase<T>::mSize; ++i){
 	    (*this)[i] = in[i];
 	}
@@ -228,25 +244,33 @@ namespace MultiArrayHelper
     }
 
     template <typename T>
-    BlockResult<T>& BlockResult<T>::operator+=(const BlockBase& in)
+    BlockResult<T>& BlockResult<T>::assing(const T& val)
+    {
+	mRes.assing(BB::mSize, val);
+	return *this;
+    }
+    
+    /*
+    template <typename T>
+    BlockResult<T>& BlockResult<T>::operator+=(const BlockBase<T>& in)
     {
 	return operateSelf<std::plus<T> >(in);
     }
 
     template <typename T>
-    BlockResult<T>& BlockResult<T>::operator-=(const BlockBase& in)
+    BlockResult<T>& BlockResult<T>::operator-=(const BlockBase<T>& in)
     {
 	return operateSelf<std::minus<T> >(in);
     }
 
     template <typename T>
-    BlockResult<T>& BlockResult<T>::operator*=(const BlockBase& in)
+    BlockResult<T>& BlockResult<T>::operator*=(const BlockBase<T>& in)
     {
 	return operateSelf<std::multiplies<T> >(in);
     }
 
     template <typename T>
-    BlockResult<T>& BlockResult<T>::operator/=(const BlockBase& in)
+    BlockResult<T>& BlockResult<T>::operator/=(const BlockBase<T>& in)
     {
 	return operateSelf<std::divides<T> >(in);
     }
@@ -255,14 +279,14 @@ namespace MultiArrayHelper
     template <class OpFunction>
     BlockResult<T>& BlockResult<T>::operateSelf(const BlockBase<T>& in)
     {
-	assert(mSize == in.size());
+	assert(BB::mSize == in.size());
 	OpFunction f;
 	//BlockResult<T> res(mSize);
-	for(size_t i = 0; i != mSize; ++i){
+	for(size_t i = 0; i != BB::mSize; ++i){
 	    (*this)[i] = f((*this)[i], in[i]);
 	}
 	return *this;
     }
-
+    */
     
 } // end namespace MultiArrayHelper

src/block.h

@@ -29,6 +29,17 @@ namespace MultiArrayHelper
 	BlockResult<T> operator()(const BlockClass1& arg1, const BlockClass2& arg2);
     };
 
+    template <typename T, class OpFunc, class BlockClass>
+    class BlockBinaryOpSelf
+    {
+    public:
+	BlockBinaryOpSelf(BlockResult<T>& res);
+	void operator()(const BlockClass& arg);
+    private:
+	BlockResult<T>& mRes;
+    };
+
+    
     // EVERYTHING IN HERE MUST  N O T  BE VITUAL !!
     
     template <typename T>
@@ -39,7 +50,7 @@ namespace MultiArrayHelper
 	BlockBase(size_t size);
 
 	size_t size() const;
-	
+	/*
 	template <class OpFunction>
 	BlockResult<T> operate(const BlockBase& in);
 	
@@ -47,7 +58,7 @@ namespace MultiArrayHelper
 	BlockResult<T> operator-(const BlockBase& in);
 	BlockResult<T> operator*(const BlockBase& in);
 	BlockResult<T> operator/(const BlockBase& in);
-	
+	*/
     protected:
 	size_t mSize;
     };
@@ -116,12 +127,16 @@ namespace MultiArrayHelper
     class BlockResult : public MutableBlockBase<T>
     {
     public:
+	typedef BlockBase<T> BB;
+	
 	DEFAULT_MEMBERS(BlockResult);
 	
 	BlockResult(size_t size);
 	
 	template <class BlockClass>
 	BlockResult& operator=(const BlockClass& in);
+
+	BlockResult& assing(const T& val);
 	
 	BlockType type() const;
 	const T& operator[](size_t pos) const;
@@ -129,13 +144,13 @@ namespace MultiArrayHelper
 	BlockResult& set(size_t npos);
 	size_t stepSize() const;
 	
-	BlockResult<T>& operator+=(const BlockBase& in);
+	//BlockResult<T>& operator+=(const BlockBase<T>& in);
-	BlockResult<T>& operator-=(const BlockBase& in);
+	//BlockResult<T>& operator-=(const BlockBase<T>& in);
-	BlockResult<T>& operator*=(const BlockBase& in);
+	//BlockResult<T>& operator*=(const BlockBase<T>& in);
-	BlockResult<T>& operator/=(const BlockBase& in);
+	//BlockResult<T>& operator/=(const BlockBase<T>& in);
 
-	template <class OpFunction>
+	//template <class OpFunction>
-	BlockResult<T>& operateSelf(const BlockBase& in);
+	//BlockResult<T>& operateSelf(const BlockBase<T>& in);
 	
     protected:
 	std::vector<T> mRes;

src/multi_array_operation.cc

@@ -160,7 +160,7 @@ namespace MultiArrayTools
     {
     	return Operation<T,std::divides<T>,OperationClass,Second>(*mOc, in);
     }
-
+    /*
     template <typename T, class OperationClass>
     template <class IndexType>
     auto OperationTemplate<T,OperationClass>::c(std::shared_ptr<IndexType>& ind) const
@@ -168,7 +168,7 @@ namespace MultiArrayTools
     {
 	return Contraction<T,OperationClass,IndexType>(*mOc, ind);
     }
-
+    */
     
     /*************************
      *   OperationMaster     *
@@ -324,6 +324,7 @@ namespace MultiArrayTools
     const BlockResult<T>& Operation<T,OpFunction,Ops...>::get() const
     {
 	mRes = std::move( PackNum<sizeof...(Ops)-1>::template unpackArgs<T,OpFunction>(mOps) );
+	//CHECK;
 	return mRes;
     }
 
@@ -347,15 +348,18 @@ namespace MultiArrayTools
      *********************/
 
     template <typename T, class Op, class IndexType>
-    Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
+    Contraction<T,Op,IndexType>::Contraction(const Op& op, std::shared_ptr<IndexType> ind) :
 	OperationTemplate<T,Contraction<T,Op,IndexType> >(this),
 	mOp(op) {}
 
-    const BlockResult<T>& get() const
+    template <typename T, class Op, class IndexType>
+    const BlockResult<T>& Contraction<T,Op,IndexType>::get() const
     {
-	// set mRes = 0 !!!
+	BlockBinaryOpSelf<T,std::plus<T>,BlockResult<T> > f(mRes);
-	for(*mIndex = 0; mIndex->pos() != mIndex->max(); ++(*mIndex)){
+	mRes.assign( static_cast<T>(0) );
-	    mRes += mOp.get();
+	for(*mInd = 0; mInd->pos() != mInd->max(); ++(*mInd)){
+	    //mRes += mOp.get();
+	    f(mOp.get());
 	}
 	return mRes;
     }
@@ -364,12 +368,12 @@ namespace MultiArrayTools
     std::vector<BTSS> Contraction<T,Op,IndexType>::block(const std::shared_ptr<IndexBase> blockIndex) const
     {
 	std::vector<BTSS> btv;
-	PackNum<sizeof...(Ops)-1>::makeBlockTypeVec(btv, mOp, blockIndex);
+	PackNum<0>::makeBlockTypeVec(btv, std::make_tuple( mOp ), blockIndex);
 	return btv;
     }
 
     template <typename T, class Op, class IndexType>
-    const Contraction& block() const
+    const Contraction<T,Op,IndexType>& Contraction<T,Op,IndexType>::block() const
     {
 	return *this;
     }

src/multi_array_operation.h

@@ -108,11 +108,11 @@ namespace MultiArrayTools
 	template <class Second>
 	auto operator/(const Second& in) const
 	    -> Operation<T,std::divides<T>,OperationClass,Second>;
-
+	/*
 	template <class IndexType>
 	auto c(std::shared_ptr<IndexType>& ind) const
 	    -> Contraction<T,OperationClass,IndexType>;
-	
+	*/
     private:
 	OperationClass* mOc;
     };
@@ -226,9 +226,9 @@ namespace MultiArrayTools
 	std::tuple<Ops...> mOps;
 	mutable BlockResult<T> mRes;
     };
-
+    
     template <typename T, class Op, class IndexType>
-    class Contraction : public OperationTemplate<T,Contraction<T,Op> >
+    class Contraction : public OperationTemplate<T,Contraction<T,Op,IndexType> >
     {
     public:
 

src/op_unit_test.cc

@@ -198,6 +198,7 @@ namespace {
 	std::cout << "ratio: " << static_cast<double>( end - begin ) / static_cast<double>( end2 - begin2 ) << std::endl;
 
 	EXPECT_EQ( xround( res.at(mkt(7,9)) ), xround(res2[7*vs1 + 9]) );
+	//EXPECT_EQ( xround( res.at(mkt(700,900)) ), xround(res2[700*vs1 + 900]) );
 	
     }