lutest + corresponding fixes

src/include/array/array_base.cc.h

@@ -30,8 +30,8 @@ namespace CNORXZ
     template <typename I, typename M>
     const T& CArrayBase<T>::operator[](const IndexInterface<I,M>& i) const
     {
-	if(formatIsTrivial()){
+	if(formatIsTrivial() or not i.formatIsTrivial()){
-	    return data()[i.lex()];
+	    return data()[i.pos()];
 	}
 	else {
 	    auto ai = itLex(i);
@@ -130,19 +130,12 @@ namespace CNORXZ
 	else {
 	    if(i->formatIsTrivial()){
 		// try to apply container format.
-		// if the reformat changes the index type in any manner
-		// the format is not applicable:
-		if constexpr(std::is_same<decltype(i->reformat( Vector<SizeT>(), Vector<SizeT>() )),Index>::value){
 		auto beg = begin();
 		auto aformat = beg.deepFormat();
 		auto amax = beg.deepMax();
-		    auto fi = i->reformat( aformat, amax );
+		//auto j = std::make_shared<Index>(*i); // use copy, otherwise would change i
-		    return coproot(*this, moveToPtr( fi ) );
+		i->reformat( aformat, amax );
-		}
+		return coproot(*this, i );
-		else {
-		    this->checkFormatCompatibility(*i);
-		    return coproot(*this, i);
-		}
 	    }
 	    else {
 		// check if format is compatible
@@ -192,7 +185,7 @@ namespace CNORXZ
 	auto j = begin();
 	CXZ_ASSERT(acc.lmax().val() == j.lmax().val(),
 		   "got index of iteration space size = " << acc.lmax().val()
-		   << ", expected size = " << acc.lmax().val());
+		   << ", expected size = " << j.lmax().val());
 	Vector<SizeT> f1 = toVec(acc.deepFormat());
 	Vector<SizeT> f2 = j.deepFormat();
 	std::sort(f1.begin(),f1.end());
@@ -234,9 +227,8 @@ namespace CNORXZ
     template <typename I, typename M>
     T& ArrayBase<T>::operator[](const IndexInterface<I,M>& i)
     {
-	// TODO: if format of i is non-trivial, use that format without check (checks only in at()!)
+	if(this->formatIsTrivial() or not i.formatIsTrivial()){
-	if(this->formatIsTrivial()){
+	    return data()[i.pos()];
-	    return data()[i.lex()];
 	}
 	else {
 	    auto ai = itLex(i);
@@ -323,19 +315,12 @@ namespace CNORXZ
 	else {
 	    if(i->formatIsTrivial()){
 		// try to apply container format.
-		// if the reformat changes the index type in any manner
-		// the format is not applicable:
-		if constexpr(std::is_same<decltype(i->reformat( Vector<SizeT>(), Vector<SizeT>() )),Index>::value){
 		auto beg = begin();
 		auto aformat = beg.deepFormat();
 		auto amax = beg.deepMax();
-		    auto fi = i->reformat( aformat, amax );
+		//auto j = std::make_shared<Index>(*i); // use copy, otherwise would change i
-		    return oproot(*this, moveToPtr( fi ) );
+		i->reformat( aformat, amax );
-		}
+		return oproot(*this, i );
-		else {
-		    this->checkFormatCompatibility(*i);
-		    return oproot(*this, i);
-		}
 	    }
 	    else {
 		// check if format is compatible

src/include/array/array_base.h

@@ -45,7 +45,9 @@ namespace CNORXZ
 	/** default destructor */
 	virtual ~CArrayBase() = default;
 
-	/** const data element access
+	/** const data element access.
+	    If the array's format is trivial or the index is already non-trivially formatted,
+	    this is equivalent to data()[i.pos()]. Otherwise, the array's format is applied.
 	    @tparam I index type
 	    @tparam M meta data type
 	    @param i index
@@ -100,7 +102,8 @@ namespace CNORXZ
 	Sptr<CArrayBase<T>> sl(const IndexInterface<I,M>& begin,
 			       const IndexInterface<I,M>& end) const;
 
-	/** create operation on this container
+	/** create operation on this container.
+	    Caution: might modify the index format.
 	    @tparam Index type of operation index
 	    @param i operation index
 	 */

src/include/base/utils.h

@@ -36,78 +36,113 @@ namespace CNORXZ
 	return Vector<T> { a };
     }
 
+    template <typename T, typename U>
+    struct Concat
+    {
+	static constexpr Arr<T,2> cat(const T& a1, const U& a2)
+	{
+	    static_assert( std::is_same<T,U>::value, "types have to be vector, array or equal" );
+	    return Arr<T,2> { a1, a2 };
+	}
+    };
+
     template <typename T, SizeT N1, SizeT N2>
-    constexpr Arr<T,N1+N2> cat2(const Arr<T,N1>& a1, const Arr<T,N2>& a2)
+    struct Concat<Arr<T,N1>,Arr<T,N2>>
+    {
+	static constexpr Arr<T,N1+N2> cat(const Arr<T,N1>& a1, const Arr<T,N2>& a2)
 	{
 	    return iter<0,N1+N2>
 		( [&](auto i) { if constexpr(i < N1) { return std::get<i>(a1); } else { return std::get<i-N1>(a2); } },
 		  [](const auto&... e) { return Arr<T,N1+N2> { e... }; } );
 	}
+    };
 
     template <typename T, SizeT N1>
-    constexpr Arr<T,N1+1> cat2(const Arr<T,N1>& a1, const T& a2)
+    struct Concat<Arr<T,N1>,T>
+    {
+	static constexpr Arr<T,N1+1> cat(const Arr<T,N1>& a1, const T& a2)
 	{
 	    return iter<0,N1>
 		( [&](auto i) { return std::get<i>(a1); },
 		  [&](const auto&... e) { return Arr<T,N1+1> { e..., a2 }; } );
 	}
+    };
 
     template <typename T, SizeT N1>
-    constexpr Arr<T,N1+1> cat2(const T& a1, const Arr<T,N1>& a2)
+    struct Concat<T,Arr<T,N1>>
+    {
+	static constexpr Arr<T,N1+1> cat(const T& a1, const Arr<T,N1>& a2)
 	{
 	    return iter<0,N1>
 		( [&](auto i) { return std::get<i>(a2); },
 		  [&](const auto&... e) { return Arr<T,N1+1> { a1, e... }; } );
 	}
+    };
 
-    template <typename T>
-    constexpr Arr<T,2> cat2(const T& a1, const T& a2)
-    {
-	return Arr<T,2> { a1, a2 };
-    }
 
     template <typename T, SizeT N2>
-    Vector<T> cat2(const Vector<T>& a1, const Arr<T,N2>& a2)
+    struct Concat<Vector<T>,Arr<T,N2>>
+    {
+	static Vector<T> cat(const Vector<T>& a1, const Arr<T,N2>& a2)
 	{
 	    Vector<T> o(a1.size()+N2);
 	    std::copy(a1.begin(), a1.end(), o.begin());
 	    std::copy(a2.begin(), a2.end(), o.begin()+a1.size());
 	    return o;
 	}
+    };
 
     template <typename T, SizeT N1>
-    Vector<T> cat2(const Arr<T,N1>& a1, const Vector<T>& a2)
+    struct Concat<Arr<T,N1>,Vector<T>>
+    {
+	static Vector<T> cat(const Arr<T,N1>& a1, const Vector<T>& a2)
 	{
 	    Vector<T> o(N1+a2.size());
 	    std::copy(a1.begin(), a1.end(), o.begin());
 	    std::copy(a2.begin(), a2.end(), o.begin()+N1);
 	    return o;
 	}
+    };
 
     template <typename T>
-    Vector<T> cat2(const Vector<T>& a1, const Vector<T>& a2)
+    struct Concat<Vector<T>,Vector<T>>
+    {
+	static Vector<T> cat(const Vector<T>& a1, const Vector<T>& a2)
 	{
 	    Vector<T> o(a1.size()+a2.size());
 	    std::copy(a1.begin(), a1.end(), o.begin());
 	    std::copy(a2.begin(), a2.end(), o.begin()+a1.size());
 	    return o;
 	}
+    };
 
     template <typename T>
-    Vector<T> cat2(const Vector<T>& a1, const T& a2)
+    struct Concat<Vector<T>,T>
+    {
+	static Vector<T> cat(const Vector<T>& a1, const T& a2)
 	{
 	    Vector<T> o(a1);
 	    o.push_back(a2);
 	    return o;
 	}
+    };
 
-    template <typename T, SizeT N1>
+    template <typename T>
-    Vector<T> cat2(const T& a1, const Vector<T>& a2)
+    struct Concat<T,Vector<T>>
+    {
+	static Vector<T> cat(const T& a1, const Vector<T>& a2)
 	{
 	    Vector<T> o { a1 };
 	    o.insert(o.end(), a2.begin(), a2.end());
 	    return o;
 	}
+    };
+
+    template <typename T1, typename T2>
+    decltype(auto) cat2(const T1& a1, const T2& a2)
+    {
+	return Concat<T1,T2>::cat(a1, a2);
+    }
 
     template <typename T1, typename T2, typename... Ts>
     decltype(auto) concat(const T1& a1, const T2& a2, const Ts&... as)

src/include/operation/op_types.cc.h

@@ -509,7 +509,41 @@ namespace CNORXZ
 	return i->xpr(i);
     }
 
+    /*============+
+     |    MapOp   |
+     +============*/
 
+    template <class Index, class F>
+    constexpr MapOp<Index,F>::MapOp(const Sptr<Index>& i, F&& f) :
+	mI(i),
+	mF(std::forward<F>(f))
+    {}
+
+    template <class Index, class F>
+    template <class PosT>
+    constexpr decltype(auto) MapOp<Index,F>::operator()(const PosT& pos) const
+    {
+	return mF(pos.val());
+    }
+    
+    template <class Index, class F>
+    constexpr decltype(auto) MapOp<Index,F>::operator()() const
+    {
+	return mF(SPos<0>());
+    }
+
+    template <class Index, class F>
+    template <SizeT I>
+    constexpr decltype(auto) MapOp<Index,F>::rootSteps(const IndexId<I>& id) const
+    {
+	return mI->stepSize(id);
+    }
+
+    template <class Index, class F>
+    constexpr decltype(auto) mapop(const Sptr<Index>& i, F&& f)
+    {
+	return MapOp<Index,F>(i, std::forward<F>(f));
+    }
 }
 
 #endif

src/include/operation/op_types.h

@@ -290,9 +290,35 @@ namespace CNORXZ
     template <class F, class Op, class IndexT>
     constexpr decltype(auto) contraction(F&& f, Op&& op, const Sptr<IndexT>& i);
 
+    // redundant? see xpr()
     template <class IndexT>
     constexpr decltype(auto) indexOp(const Sptr<IndexT>& i);
 
+    template <class Index, class F>
+    class MapOp : public OpInterface<MapOp<Index,F>>
+    {
+    public:
+	typedef OpInterface<MapOp<Index,F>> OI;
+
+	constexpr MapOp() = default;
+	
+	constexpr MapOp(const Sptr<Index>& i, F&& f);
+
+	template <class PosT>
+	constexpr decltype(auto) operator()(const PosT& pos) const;
+
+	constexpr decltype(auto) operator()() const;
+
+	template <SizeT I>
+	constexpr decltype(auto) rootSteps(const IndexId<I>& id) const;
+
+    private:
+	Sptr<Index> mI;
+	F mF;
+    };
+
+    template <class Index, class F>
+    constexpr decltype(auto) mapop(const Sptr<Index>& i, F&& f);
 }
 
 #endif

src/include/ranges/mrange.cc.h

@@ -79,7 +79,7 @@ namespace CNORXZ
     {
 	return gmformat
 	    ( iter<Is,NI>
-	      ( [&](auto i) { return ipack[i]->pmax(); },
+	      ( [&](auto i) { return ipack[i]->lmax(); },
 		[](const auto&... as) { return (as * ...); } )...,
 	      SPos<1>() );
     }
@@ -637,6 +637,12 @@ namespace CNORXZ
 	return std::make_shared<GMIndex<FormatT,Indices...>>(bs, pack);
     }    
     
+    template <class... Indices>
+    constexpr decltype(auto) gmindexPtr(const SPack<Indices...>& pack)
+    {
+	return std::make_shared<GMIndex<MFormat<sizeof...(Indices)>,Indices...>>(pack);
+    }
+
     template <class I1, class FormatT, class... Indices>
     decltype(auto) operator*(const Sptr<GMIndex<FormatT,Indices...>>& a, const Sptr<I1>& b)
     {

src/include/ranges/mrange.h

@@ -299,6 +299,12 @@ namespace CNORXZ
     template <class FormatT, class... Indices>
     constexpr decltype(auto) gmindexPtr(const FormatT& bs, const SPack<Indices...>& pack);
 
+    /** Create pointer to GMIndex from index pack (default format).
+	@param pack Pack of input indices.
+     */
+    template <class... Indices>
+    constexpr decltype(auto) gmindexPtr(const SPack<Indices...>& pack);
+
     /** Specialization for index multiplication with GMIndex on the l.h.s.
 	@param a First index of type GMIndex.
 	@param b Second index of arbitrary type.

src/tests/l_unit_test.cc

@@ -38,7 +38,7 @@ namespace
 	    mLocTR = std::dynamic_pointer_cast<CRange> ( CRangeFactory(8).create() );
 	    mRaSR = std::dynamic_pointer_cast<CRange> ( CRangeFactory(2).create() );
 	    mRaTR = std::dynamic_pointer_cast<CRange> ( CRangeFactory(3).create() );
-	    mData = MArray<Double>(yrange({mRaTR,mRaSR,mRaSR,mRaSR,mLocTR,mLocSR,mLocSR,mLocSR,mSpinR,mColorR,mSpinR,mColorR}));
+	    mData = MArray<Arr<Int,4>>(yrange({mRaTR,mRaSR,mRaSR,mRaSR,mLocTR,mLocSR,mLocSR,mLocSR,mSpinR,mColorR,mSpinR,mColorR}));
 	    mDataFormat = mData.begin().format();
 	    auto vf = mDataFormat.all();
 	    vf[0] = mDataFormat[0];
@@ -50,7 +50,13 @@ namespace
 	    vf[6] = mDataFormat[3];
 	    vf[7] = mDataFormat[7];
 	    mViewFormat = YFormat(vf);
-	    mView = Slice<Double>(yrange({mRaTR,mLocTR,mRaSR,mLocSR,mRaSR,mLocSR,mRaSR,mLocSR,mSpinR,mColorR,mSpinR,mColorR}), &mData, mViewFormat, 0);
+	    mView = Slice<Arr<Int,4>>(yrange({mRaTR,mLocTR,mRaSR,mLocSR,mRaSR,mLocSR,mRaSR,mLocSR,mSpinR,mColorR,mSpinR,mColorR}), &mData, mViewFormat, 0);
+	    mLocSI = CIndex(mLocSR);
+	    mLocTI = CIndex(mLocTR);
+	    mRaSI = CIndex(mRaSR);
+	    mRaTI = CIndex(mRaTR);
+	    mSpinI = SIndex<SizeT,4>(mSpinR);
+	    mColorI = SIndex<SizeT,3>(mColorR);
 	}
 	
 	Sptr<SRange<SizeT,4>> mSpinR;
@@ -62,16 +68,69 @@ namespace
 	Sptr<CRange> mRaSR;
 	Sptr<CRange> mRaTR;
 
-	MArray<Double> mData;
+	MArray<Arr<Int,4>> mData;
-	Slice<Double> mView;
+	Slice<Arr<Int,4>> mView;
 	YFormat mDataFormat;
 	YFormat mViewFormat;
+
+	CIndex mLocSI;
+	CIndex mLocTI;
+	CIndex mRaSI;
+	CIndex mRaTI;
+	SIndex<SizeT,4> mSpinI;
+	SIndex<SizeT,3> mColorI;
+	MIndex<CIndex,CIndex> mSpatialI;
+	MIndex<CIndex,CIndex> mTemporalI;
+	//GMIndex<Arr<SizeT,2>,CIndex,CIndex> mSpatialDI;
+	//GMIndex<Arr<SizeT,2>,CIndex,CIndex> mTemporalDI;
     };
 
     TEST_F(L_Test, Basic)
     {
 	EXPECT_EQ(mData.size(), mView.size());
-	VCHECK(mData.size());
     }
 
+    TEST_F(L_Test, Assign)
+    {
+	auto x0l = std::make_shared<CIndex>(mLocTI);
+	auto x0r = std::make_shared<CIndex>(mRaTI);
+	auto x0 = gmindexPtr(x0r*x0l);
+	auto x0a = gmindexPtr(x0r*x0l);
+	auto x1l = std::make_shared<CIndex>(mLocSI);
+	auto x1r = std::make_shared<CIndex>(mRaSI);
+	auto x1 = gmindexPtr(x1r*x1l);
+	auto x1a = gmindexPtr(x1r*x1l);
+	auto x2l = std::make_shared<CIndex>(mLocSI);
+	auto x2r = std::make_shared<CIndex>(mRaSI);
+	auto x2 = gmindexPtr(x2r*x2l);
+	auto x2a = gmindexPtr(x2r*x2l);
+	auto x3l = std::make_shared<CIndex>(mLocSI);
+	auto x3r = std::make_shared<CIndex>(mRaSI);
+	auto x3 = gmindexPtr(x3r*x3l);
+	auto x3a = gmindexPtr(x3r*x3l);
+	auto x = gmindexPtr(x0*x1*x2*x3);
+	auto xx = gmindexPtr(x0a*x1a*x2a*x3a);
+	auto al = std::make_shared<SIndex<SizeT,4>>(mSpinI);
+	auto be = std::make_shared<SIndex<SizeT,4>>(mSpinI);
+	auto a = std::make_shared<SIndex<SizeT,3>>(mColorI);
+	auto b = std::make_shared<SIndex<SizeT,3>>(mColorI);
+
+	EXPECT_TRUE(x0->formatIsTrivial());
+	EXPECT_TRUE(x1->formatIsTrivial());
+	EXPECT_TRUE(x2->formatIsTrivial());
+	EXPECT_TRUE(x3->formatIsTrivial());
+	EXPECT_TRUE(x->formatIsTrivial());
+
+	auto fs = [&](SizeT i) { return i >= mSpatialR->size()/2 ? static_cast<Int>(i - mSpatialR->size()) : static_cast<Int>(i); };
+	auto ft = [&](SizeT i) { return i >= mTemporalR->size()/2 ? static_cast<Int>(i - mTemporalR->size()) : static_cast<Int>(i); };
+        EXPECT_TRUE(mindexPtr(xx*al*a*be*b)->formatIsTrivial());
+
+	mView(gmindexPtr(xx*al*a*be*b)) = operation( [](Int a0, Int a1, Int a2, Int a3) { Arr<Int,4> x{a0,a1,a2,a3}; return x; },
+						     mapop(x0, ft), mapop(x1, fs), mapop(x2, fs), mapop(x3, fs) );
+	EXPECT_TRUE(x->formatIsTrivial());
+	for(*x = 0; x->lex() != x->lmax().val(); ++(*x)){
+	    Arr<Int,4> m{ ft(x0->lex()), fs(x1->lex()), fs(x2->lex()), fs(x3->lex()) };
+	    EXPECT_EQ( mView[x*al*a*be*b], m );
+	}
+    }
 }