spin range contraction test + some simplifying tools

src/container_range.h

@@ -188,7 +188,7 @@ namespace MultiArrayTools
 
 	typedef ContainerRange<Ranges...> oType;
 	
-	ContainerRangeFactory() = delete;
+	ContainerRangeFactory();
 	ContainerRangeFactory(const std::shared_ptr<Ranges>&... rs);
 	ContainerRangeFactory(const typename ContainerRange<Ranges...>::SpaceType& space);
 	
@@ -239,6 +239,7 @@ namespace MultiArrayTools
 
 	friend ContainerRangeFactory<Ranges...>;
 
+	static const bool defaultable = false;
     };
    
 } // end namespace MultiArrayTools

src/helper_tools.h

@@ -11,6 +11,19 @@ namespace MultiArrayTools
     auto getIndex(std::shared_ptr<RangeType> range)
 	-> std::shared_ptr<typename RangeType::IndexType>;
 
+    // only if 'RangeType' is defaultable and unique (Singleton)
+    template <class RangeType>
+    auto getIndex() -> std::shared_ptr<typename RangeType::IndexType>;
+
+    template <class... RangeTypes>
+    auto mkMulti(std::shared_ptr<RangeTypes>... ranges)
+	-> std::shared_ptr<MultiRange<RangeTypes...> >;
+
+    template <class... IndexTypes>
+    auto mkMIndex(std::shared_ptr<IndexTypes>... indices)
+	-> decltype( getIndex( mkMulti( indices.range()... ) ) );
+
+    
 }
 
 /* ========================= *
@@ -26,6 +39,34 @@ namespace MultiArrayTools
 	return std::dynamic_pointer_cast<IndexWrapper<typename RangeType::IndexType> >
 	    ( range->index() )->get();
     }
+
+    template <class RangeType>
+    auto getIndex() -> std::shared_ptr<typename RangeType::IndexType>
+    {
+	static_assert( RangeType::defaultable,
+		       /*typeid(typename RangeType).name() + */" is not defaultable" );
+	static auto f = RangeType::factory();
+	static auto r = std::dynamic_pointer_cast<RangeType>( f.create() );
+	return std::dynamic_pointer_cast<IndexWrapper<typename RangeType::IndexType> >
+	    ( r->index() )->get();
+    }
+    
+    template <class... RangeTypes>
+    auto mkMulti(std::shared_ptr<RangeTypes>... ranges)
+	-> std::shared_ptr<MultiRange<RangeTypes...> >
+    {
+	MultiRangeFactory<RangeTypes...> mrf( ranges... );
+	return std::dynamic_pointer_cast<MultiRange<RangeTypes...> >( mrf.create() );
+    }
+
+    template <class... IndexTypes>
+    auto mkMIndex(std::shared_ptr<IndexTypes>... indices)
+	-> decltype( getIndex( mkMulti( indices->range()... ) ) )
+    {
+	auto mi = getIndex( mkMulti( indices->range()... ) );	    
+	(*mi)( indices... );
+	return mi;
+    }
 }
 
 #endif

src/multi_array.h

@@ -202,6 +202,10 @@ namespace MultiArrayTools
 	MultiArray(const std::shared_ptr<SRanges>&... ranges, const std::vector<T>& vec);
 	MultiArray(const std::shared_ptr<SRanges>&... ranges, std::vector<T>&& vec);
 
+	// Only if ALL ranges have default extensions:
+	//MultiArray(const std::vector<T>& vec);
+	//MultiArray(std::vector<T>&& vec);
+	
 	// template <class Range2, class Range3>
 	// MultiArray(const MultiArray<MultiArray<T,Range2>,Range3> in);
 
@@ -716,6 +720,7 @@ namespace MultiArrayTools
 	    mCont.erase(mCont.begin() + MAB::mRange->size(), mCont.end());
 	}
     }
+
     /*
     template <typename T, class... SRanges>
     template <class Range2, class Range3>

src/multi_range.h

@@ -251,6 +251,7 @@ namespace MultiArrayTools
 
 	friend MultiRangeFactory<Ranges...>;
 
+	static const bool defaultable = false;
     };
     
 }

src/op_unit_test.cc

@@ -57,6 +57,12 @@ namespace {
 	return std::make_tuple(ts...);
     }
 
+    template <typename... Ts>
+    auto mkts(Ts&&... ts) -> decltype(std::make_tuple(ts...))
+    {
+	return std::make_tuple(static_cast<size_t>( ts )...);
+    }
+    
     class OpTest_1Dim : public ::testing::Test
     {
     protected:
@@ -177,46 +183,89 @@ namespace {
     {
     protected:
 
-	typedef SpinRangeFactory SRF;
+	typedef SpinRF SRF;
 	typedef SpinRange SR;
 	typedef MultiRangeFactory<SR,SR,SR,SR,SR,SR,SR,SR> SR8F;
 	typedef SR8F::oType SR8;
 
 	static const size_t s = 65536;
 	
-	OpTest_Spin() : data(s)
+	OpTest_Spin()
 	{
+	    data.resize(s);
 	    for(size_t i = 0; i != s; ++i){
 		double arg = static_cast<double>( i - s ) - 0.1; 
 		data[i] = sin(arg)/arg;
 	    }
-
+	    SRF f;
-	    swapFactory<SRF>(rfbptr);
+	    sr = std::dynamic_pointer_cast<SR>(f.create());
-	    srptr = std::dynamic_pointer_cast<SR>( rfbptr->create() );
-	    /*
-	    swapMFactory<SR8F>(rfbptr, srptr, srptr, srptr, srptr,
-			       srptr, srptr, srptr, srptr);
-	    mrptr = std::dynamic_pointer_cas<SR8>( rfbptr->create() );
-	    */
 	}
 
-	std::shared_ptr<RangeFactoryBase> rfbptr;
-	std::shared_ptr<SR> srptr;
-	//std::shared_ptr<SR8> mrptr;
 	std::vector<double> data;
+	std::shared_ptr<SR> sr;
     };
 
     TEST_F(OpTest_Spin, Contract)
     {
-	MultiArray<double,SR,SR,SR,SR,SR,SR,SR,SR> ma(srptr, srptr, srptr, srptr,
+	MultiArray<double,SR,SR,SR,SR,SR,SR,SR,SR> ma(sr, sr, sr, sr, sr, sr, sr, sr, data);
-						      srptr, srptr, srptr, srptr, data);
+	MultiArray<double,SR,SR> res1( sr, sr );
 	
-	auto alpha = MAT::getIndex(srptr);
+	auto alpha = MAT::getIndex<SR>();
-	auto beta = MAT::getIndex(srptr);
+	auto beta = MAT::getIndex<SR>();
-	auto gamma = MAT::getIndex(srptr);
+	auto gamma = MAT::getIndex<SR>();
-	auto delta = MAT::getIndex(srptr);
+	auto delta = MAT::getIndex<SR>();
+	auto deltap = MAT::getIndex<SR>();
 	
-	// !!!
+	auto mix = MAT::mkMIndex( alpha, beta, gamma );
+
+	std::clock_t begin = std::clock();
+	res1(delta, deltap) = ma(delta, alpha, alpha, beta, beta, gamma, gamma, deltap).c(mix);
+	std::clock_t end = std::clock();
+	std::cout << "MultiArray time: " << static_cast<double>( end - begin ) / CLOCKS_PER_SEC
+		  << std::endl;
+	
+	std::vector<double> vres(4*4);
+
+	std::clock_t begin2 = std::clock();
+	for(size_t d = 0; d != 4; ++d){
+	    for(size_t p = 0; p != 4; ++p){
+		const size_t tidx = d*4 + p;
+		vres[tidx] = 0.;
+		for(size_t a = 0; a != 4; ++a){
+		    for(size_t b = 0; b != 4; ++b){
+			for(size_t c = 0; c != 4; ++c){
+			    const size_t sidx = d*4*4*4*4*4*4*4 + a*5*4*4*4*4 + b*5*4*4*4 + c*5*4 + p;
+			    vres[tidx] += data[sidx];
+			}
+		    }
+		}
+	    }
+	}
+	std::clock_t end2 = std::clock();
+
+	EXPECT_EQ( xround(res1.at(mkts(0,0))), xround(vres[0]) );
+	EXPECT_EQ( xround(res1.at(mkts(0,1))), xround(vres[1]) );
+	EXPECT_EQ( xround(res1.at(mkts(0,2))), xround(vres[2]) );
+	EXPECT_EQ( xround(res1.at(mkts(0,3))), xround(vres[3]) );
+	
+	EXPECT_EQ( xround(res1.at(mkts(1,0))), xround(vres[4]) );
+	EXPECT_EQ( xround(res1.at(mkts(1,1))), xround(vres[5]) );
+	EXPECT_EQ( xround(res1.at(mkts(1,2))), xround(vres[6]) );
+	EXPECT_EQ( xround(res1.at(mkts(1,3))), xround(vres[7]) );
+
+	EXPECT_EQ( xround(res1.at(mkts(2,0))), xround(vres[8]) );
+	EXPECT_EQ( xround(res1.at(mkts(2,1))), xround(vres[9]) );
+	EXPECT_EQ( xround(res1.at(mkts(2,2))), xround(vres[10]) );
+	EXPECT_EQ( xround(res1.at(mkts(2,3))), xround(vres[11]) );
+
+	EXPECT_EQ( xround(res1.at(mkts(3,0))), xround(vres[12]) );
+	EXPECT_EQ( xround(res1.at(mkts(3,1))), xround(vres[13]) );
+	EXPECT_EQ( xround(res1.at(mkts(3,2))), xround(vres[14]) );
+	EXPECT_EQ( xround(res1.at(mkts(3,3))), xround(vres[15]) );
+
+	std::cout << "std::vector - for loop time: " << static_cast<double>( end2 - begin2 ) / CLOCKS_PER_SEC
+		  << std::endl;
+	std::cout << "ratio: " << static_cast<double>( end - begin ) / static_cast<double>( end2 - begin2 ) << std::endl;
     }    
     
     TEST_F(OpTest_Performance, PCheck)

src/pack_num.h

@@ -260,6 +260,13 @@ namespace MultiArrayHelper
 	    std::get<N>(ip)->print(offset);
 	    PackNum<N-1>::printIndex(ip, offset);
 	}
+
+	template <class Range, class... Ranges>
+	static void checkDefaultable()
+	{
+	    static_assert( Range::defaultable, "not defaultable" );
+	    PackNum<N-1>::template checkDefaultable<Ranges...>();
+	}
     };
     
     template<>
@@ -454,6 +461,12 @@ namespace MultiArrayHelper
 	    std::get<0>(ip)->print(offset);
 	}
 
+	template <class Range>
+	static void checkDefaultable()
+	{
+	    static_assert( Range::defaultable, "not defaultable" );
+	}
+	
     };
 
     template <typename... T>

src/range_types/header.h

@@ -10,12 +10,14 @@
 
 #ifdef __incl_this__
 
-#ifndef __ranges_header__
 #define __ranges_header__
+//#ifndef __ranges_header__
+//#define __ranges_header__
 
 #include "spin_range.h"
 
-#endif
+#undef __ranges_header__
+//#endif
 #endif
 
 #undef __incl_this__

src/range_types/spin_range.h

@@ -7,8 +7,8 @@ include_range_type(SPIN,2)
 #ifdef __ranges_header__
 // assert, that this is only used within range_types/header.h
 
-#ifndef __spin_range_h__
+//#ifndef __spin_range_h__
-#define __spin_range_h__
+//#define __spin_range_h__
 
 namespace MultiArrayTools
 {
@@ -19,10 +19,9 @@ namespace MultiArrayTools
     {
     public:
 	
-	typedef SingleRange<U,TYPE> oType;
+	typedef SingleRange<size_t,RangeType::SPIN> oType;
 
-	SingleRangeFactory() = delete;
+	SingleRangeFactory();
-	SingleRangeFactory(size_t spinNum); // = 4 :)
 	std::shared_ptr<RangeBase> create();
 	
     };
@@ -32,7 +31,7 @@ namespace MultiArrayTools
     {
     public:
 	typedef RangeBase RB;
-	typedef typename RangeInterface<SingleIndex<U,TYPE> >::IndexType IndexType;
+	typedef typename RangeInterface<SingleIndex<size_t,RangeType::SPIN> >::IndexType IndexType;
 	
 	virtual size_t size() const override;
 	virtual size_t dim() const override;
@@ -42,18 +41,22 @@ namespace MultiArrayTools
 	
 	virtual IndexType begin() const override;
 	virtual IndexType end() const override;
-	virtual std::shared_ptr<IndexBase> index() const override;
+	virtual std::shared_ptr<VIWB> index() const override;
 	
 	friend SingleRangeFactory<size_t,RangeType::SPIN>;
 
+	static const bool defaultable = true;
+	static const size_t mSpinNum = 4;
+	
+	static SingleRangeFactory<size_t, RangeType::SPIN> factory()
+	{ return SingleRangeFactory<size_t, RangeType::SPIN>(); }
+	
     protected:
 
 	SingleRange() = default;
 	SingleRange(const SingleRange& in) = delete;
 
-	SingleRange(size_t spinNum);
+	//SingleRange(size_t spinNum);
-
-	size_t mSpinNum = 4;
     };
 
     typedef SingleRange<size_t,RangeType::SPIN> SpinRange;
@@ -70,14 +73,17 @@ namespace MultiArrayTools
      *   SingleRange    *
      ********************/
     
-    SingleRangeFactory<size_t,RangeType::SPIN>::SingleRangeFactory(const std::vector<U>& space)
+    SingleRangeFactory<size_t,RangeType::SPIN>::SingleRangeFactory()
     {
-	mProd = std::shared_ptr<oType>( new SingleRange<size_t,RangeType::SPIN>( space ) );
+	// Quasi Singleton
+	if(not mProd){
+	    mProd = std::shared_ptr<oType>( new SingleRange<size_t,RangeType::SPIN>() );
+	    setSelf();
+	}
     }
 
     std::shared_ptr<RangeBase> SingleRangeFactory<size_t,RangeType::SPIN>::create()
     {
-	setSelf();
 	return mProd;
     }
     
@@ -85,12 +91,6 @@ namespace MultiArrayTools
      *   SingleRange    *
      ********************/
     
-    SingleRange<size_t,RangeType::SPIN>::SingleRange(size_t spinNum) :
-	RangeInterface<SingleIndex<size_t,RangeType::SPIN> >()
-    {
-	mSpinNum = spinNum;
-    }
-        
     size_t SingleRange<size_t,RangeType::SPIN>::get(size_t pos) const
     {
 	return pos;
@@ -128,15 +128,18 @@ namespace MultiArrayTools
     }
 
     // put this in the interface class !!!
-    std::shared_ptr<IndexBase> SingleRange<size_t,RangeType::SPIN>::index() const
+    std::shared_ptr<VIWB> SingleRange<size_t,RangeType::SPIN>::index() const
     {
-	return std::make_shared<SingleIndex<size_t,RangeType::SPIN> >
+	typedef IndexWrapper<IndexType> IW;
-	    ( std::dynamic_pointer_cast<SingleRange<size_t,RangeType::SPIN> >
+	return std::make_shared<IW>
-	      ( std::shared_ptr<RangeBase>( RB::mThis ) ) );
+            ( std::make_shared<IndexType>
+	      ( std::dynamic_pointer_cast<SingleRange<size_t,RangeType::SPIN> >
+		( std::shared_ptr<RangeBase>( RB::mThis ) ) ) );
     }
+
 }
 
-#endif // #ifndef __spin_range_h__
+//#endif // #ifndef __spin_range_h__
 
 #endif // #ifdef __ranges_header__
 

src/single_range.h

@@ -162,6 +162,8 @@ namespace MultiArrayTools
 	
 	friend SingleRangeFactory<U,TYPE>;
 
+	static const bool defaultable = false;
+	
     protected:
 
 	SingleRange() = delete;