From ec468f2db2b16e80e6511c9ef2135a5be027f3e4 Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Mon, 8 Apr 2024 16:57:09 +0200
Subject: [PATCH] tensor: inverse_index

---
 t.cc      | 19 ++++++++++++++++---
 tensor.cc | 54 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 tensor.h  | 12 +++++++-----
 3 files changed, 77 insertions(+), 8 deletions(-)

diff --git a/t.cc b/t.cc
index bcdf2ef..a4c9d55 100644
--- a/t.cc
+++ b/t.cc
@@ -3176,7 +3176,7 @@ for(int l=1; l<k; ++l)
 cout <<count<<" "<<binom(n,4)<<endl;
 }
 
-if(1)
+if(0)
 {
 int d,n;
 cin>>d>>n;
@@ -3200,17 +3200,30 @@ cout <<d;
 }
 
 
-if(0)
+if(1)
 {
 INDEXGROUP g;
 g.number=3;
 g.symmetry= -1;
 g.offset=1;
-g.range=3;
+g.range=10;
 
 Tensor<double> epsilon(g);
 cout <<epsilon.size()<<endl;
 
+for(LA_largeindex s=0; s<epsilon.size(); ++s)
+	{
+	SUPERINDEX I = epsilon.inverse_index(s);
+	int sign;
+	LA_largeindex ss=epsilon.index(&sign,I);
+	if(ss!=s || sign!=1) 
+		{
+		cout <<"offset = "<<s<<endl;
+		cout <<"index = "<<I<<endl;
+		laerror("Internal error in tensor index calculation");	
+		}
+	}
+
 NRVec<LA_index> I({1,2,3});
 NRVec<NRVec<LA_index> > II(I,1);
 epsilon.lhs(II)=1;
diff --git a/tensor.cc b/tensor.cc
index 2e2d375..baeb778 100644
--- a/tensor.cc
+++ b/tensor.cc
@@ -109,6 +109,60 @@ return -1;
 }
 
 
+//inverse map of group superindex to canonically ordered index list
+NRVec<LA_index>  inverse_subindex(const INDEXGROUP &g, LA_largeindex s)
+{
+NRVec<LA_index> I(g.number);
+if(g.number==1) {I[0]=s+g.offset; return I;}
+switch(g.symmetry)
+	{
+	case 0:
+		for(int i=0; i<g.number; ++i)
+			{
+			I[i] = s%g.range;
+			s /= g.range;
+			}
+		break;
+	case 1:
+		for(int i=g.number; i>0; --i)
+			{
+			I[i-1] = inverse_simplicial(i,s);
+			s -= simplicial(i,I[i-1]);
+			}
+		break;
+	case -1:
+		for(int i=g.number-1; i>=0; --i)
+                        {
+                        I[i] = i + inverse_simplicial(i+1,s);
+                        s -= simplicial(i+1,I[i]-i);
+                        }
+		break;
+	default: laerror("illegal index symmetry");
+	}
+
+if(g.offset!=0) I += g.offset;
+return I;
+}
+
+
+template<typename T>
+SUPERINDEX Tensor<T>::inverse_index(LA_largeindex s) const
+{
+SUPERINDEX I(shape.size());
+for(int g=shape.size()-1; g>=0; --g)
+	{
+	LA_largeindex groupindex;
+	if(g>0) 
+		{
+		groupindex = s/cumsizes[g];
+		s %= cumsizes[g];
+		}
+	else groupindex=s;
+	I[g] = inverse_subindex(shape[g],groupindex);
+	}
+return I;
+}
+
 
 
 template<typename T>
diff --git a/tensor.h b/tensor.h
index a099ac0..4dc5504 100644
--- a/tensor.h
+++ b/tensor.h
@@ -93,13 +93,12 @@ class Tensor {
 	NRVec<LA_largeindex> cumsizes; //cumulative sizes of symmetry index groups (a function of shape but precomputed for efficiency)
 	NRVec<T> data;	
 
-private:
+public:
 	LA_largeindex index(int *sign, const SUPERINDEX &I) const; //map the tensor indices to the position in data
 	LA_largeindex index(int *sign, const  FLATINDEX &I) const; //map the tensor indices to the position in data
 	LA_largeindex vindex(int *sign, LA_index i1, va_list args) const; //map list of indices  to the position in data
-//@@@reversed index
+	SUPERINDEX inverse_index(LA_largeindex s) const; //inefficient, but possible if needed
 
-public:
 	//constructors
 	Tensor() : myrank(0) {};
 	Tensor(const NRVec<indexgroup> &s) : shape(s), data((int)calcsize()), myrank(calcrank()) {}; //general tensor
@@ -130,13 +129,16 @@ public:
 	void get(int fd);
 
 	//@@@TODO - unwinding to full size in a specified index
+	//@@@contraction by a whole index group
 	//@@@TODO - contractions - basic and efficient? first contraction in a single index; between a given group+index in group at each tensor
 	//@@@ dvojite rekurzivni loopover s callbackem - nebo iterator s funkci next???
-	//@@@nebo inverse index function?
 	//@@@ stream i/o na zaklade tohoto
 	//@@@permuteindexgroups
-	//@@@symmetreize a group, antisymmetrize a group, expand a (anti)symmetric grtoup
+	//@@@symmetrize a group, antisymmetrize a group, expand a (anti)symmetric grtoup - obecne symmetry change krome +1 na -1 vse mozne
+	//@@@outer product
+	//@@@explicit constructors from vec mat smat and dense fourindex
 	//@@@@@@+= -= + - on same shape
+	//@@@@@@ randomize
 };