tensor: working on add_permuted_contractions

t.cc

@@ -3135,20 +3135,25 @@ generate antisymmetrization operator for Brandow diagrams in the Kucharski conve
                 #correspond to already existing antisymmetry
                 #in principle it would be more logical to printout permutation of places rather than indices
                 #but this is more convenient for further processing
+P "i/j,k l" input is 
+2
+3 1 2 2
+1 3
 */
 
 NRVec<NRVec_from1<int> > groups;
 cin >> groups;
 int ntot=0;
 for(int i=0; i<groups.size(); ++i) ntot+=groups[i].size();
-NRVec_from1<string> indices(ntot);
+NRVec_from1<INDEXNAME> indices(ntot);
-cin >>indices;
+for(int i=1; i<=ntot; ++i) sprintf(indices[i].name,"%c",'i'+i-1);
 
 PermutationAlgebra<int,int> a=general_antisymmetrizer(groups,-2,true);
+cout <<"PA = "<<a<<endl;
 for(int i=0; i<a.size(); ++i)
 	{
 	cout << (a[i].weight==1?'+':'-');
-	NRVec_from1<string> pindices=applypermutation(a[i].perm,indices,false);
+	NRVec_from1<INDEXNAME> pindices=applypermutation(a[i].perm,indices,false);
 	for(int j=1; j<=pindices.size(); ++j) cout <<pindices[j];
 	cout <<endl;
 	}
@@ -4164,7 +4169,7 @@ cout<<"resulting index order = "<<y.names<<endl;
 }
 
 
-if(1)
+if(0)
 {
 int nn;
 cin>>nn;
@@ -4208,6 +4213,51 @@ for(int k=0; k<nn; ++k) for(int l=0; l<nn; ++l) for(int m=0; m<nn; ++m) for(int
 	}
 
 cout <<"Error = "<<(z-zz).norm()<<endl;
+}
+
+if(1)
+{
+int nn;
+cin>>nn;
+int r=4;
+NRVec<INDEXGROUP> s(r);
+for(int i=0; i<r; ++i)
+        {
+        s[i].number=1;
+        s[i].symmetry=0;
+        s[i].range=nn;
+        s[i].offset=0;
+        }
+Tensor<double> x(s); x.randomize(1.);
+INDEXNAME xlist[] = {"i","j","k","l"};
+x.names=NRVec<INDEXNAME>(xlist);
+Tensor<double> y(s); y.randomize(1.);
+INDEXNAME ylist[] = {"n","m","j","i"};
+y.names=NRVec<INDEXNAME>(ylist);
+
+Tensor<double>z(s); z.clear();
+INDEXNAME zlist[] = {"k","l","m","n"};
+z.names=NRVec<INDEXNAME>(zlist);
+Tensor<double>zz(z);
+
+
+z.add_permuted_contractions("k/l,m",x,y,1,0,false,false);
+
+zz.copyonwrite();
+for(int k=0; k<nn; ++k) for(int l=0; l<nn; ++l) for(int m=0; m<nn; ++m) for(int n=0; n<nn; ++n)
+	{
+	double s=0;
+	for(int i=0; i<nn; ++i) for(int j=0; j<nn; ++j) s += x(i,j,k,l)*y(n,m,j,i);
+	zz.lhs(k,l,m,n) = s;
+	}
+
+Tensor<double> zzz(s);
+for(int k=0; k<nn; ++k) for(int l=0; l<nn; ++l) for(int m=0; m<nn; ++m) for(int n=0; n<nn; ++n)
+	{
+	zzz.lhs(k,l,m,n) = 1/3.*(zz(k,l,m,n)-zz(l,k,m,n)+zz(m,k,l,n));
+	}
+
+cout <<"Error = "<<(z-zzz).norm()<<endl;
 
 }
 

tensor.cc

@@ -23,6 +23,7 @@
 #include "miscfunc.h"
 #include <complex>
 #include "bitvector.h"
+#include <string.h>
 
 
 namespace LA {
@@ -1908,6 +1909,109 @@ s>>x.group>>x.index;
 return s;
 }
 
+static NRPerm<int> find_indexperm(const NRVec<INDEXNAME> &what, const NRVec<INDEXNAME> &where)
+{
+if(what.size()!=where.size()) laerror("sizes of index name vectors do not match in find_indexperm");
+NRPerm<int> basicperm(what.size());
+basicperm.clear();
+for(int i=0; i<what.size(); ++i) 
+	{
+	for(int j=0; j<what.size(); ++j)
+             if(what[i] == where[j])
+		{
+		basicperm[i+1]=j+1;
+		break;
+		}
+	}
+if(!basicperm.is_valid()) laerror("indices mismatch between lhs and rhs in add_permuted_contractions");
+return basicperm;
+}
+
+
+static void parse_antisymmetrizer(const char *txt0, NRVec<NRVec_from1<int> > &groups0, NRVec<INDEXNAME> &names0)
+{
+int igroup= -1;
+int iname= -1;
+int iclass=0;
+
+char *txt=strdup(txt0);
+char *p=txt;
+std::list<std::list<int> > groups;
+std::list<INDEXNAME> names;
+std::list<int> currentgroup;
+
+foundspace:
+	while(isspace(*p)) ++p;
+	if(*p==0) goto finished;
+	if(!isalpha(*p)) laerror("letter expected in parse_antisymmetrizer");
+	++iname;
+	++igroup;
+	++iclass;
+foundname:
+	{
+	INDEXNAME n;
+	int i=0;
+	while(isalnum(*p)) 
+		{
+		n.name[i]=*p;
+		++i;
+		if(i>N_INDEXNAME) laerror("too long index name in parse_antisymmetrizer");
+		++p;
+		}
+	if(i>=N_INDEXNAME) laerror("too long index name in parse_antisymmetrizer");
+	n.name[i]=0;
+
+	//store
+	names.push_back(n);	
+	currentgroup.push_back(iclass);
+	
+	//now decide based on next character
+	if(*p==',')
+		{
+		++iname;
+		++p;
+		if(!isalpha(*p)) laerror("letter expected in parse_antisymmetrizer");
+		goto foundname;
+		}
+	if(*p=='/')
+		{
+		++iname;
+		++p;
+		++iclass;
+		if(!isalpha(*p)) laerror("letter expected in parse_antisymmetrizer");
+                goto foundname;
+		}
+	if(isspace(*p))
+		{
+		groups.push_back(currentgroup);
+		currentgroup.clear();
+		++p;
+		goto foundspace;
+		}
+	if(*p==0) 
+		{
+		groups.push_back(currentgroup);
+		goto finished;
+		}
+	laerror("unexpected character in parse_antisymmetrizer");
+	}
+	
+finished:
+
+names0= NRVec<INDEXNAME> (names);
+std::cout <<"names parsed "<<names0;
+
+groups0.resize(groups.size());
+int i=0;
+for(typename std::list<std::list<int> >::const_iterator ii=groups.begin(); ii!=groups.end(); ++ii)
+	{
+	groups0[i++] = NRVec_from1<int>(*ii);
+	}
+std::cout<<"groups parsed "<<groups0;
+
+free(txt);
+}
+
 
 template<typename T>
 void Tensor<T>::add_permuted_contractions(const char *antisymmetrizer, const Tensor &rhs1, const Tensor &rhs2, T alpha, T beta, bool conjugate1, bool conjugate2)
@@ -1945,33 +2049,76 @@ Tensor<T> tmp=rhs1.contractions(il1,rhs2,il2,alpha,conjugate1,conjugate2);
 if(rank()!=tmp.rank()) laerror("rank mismatch in add_permuted_contractions");
 
 //generate the antisymmetrizer, adding also indices not involved as a constant subpermutation
-//@@@
+//
-PermutationAlgebra<int,T> pa(1);
+//first parse the antisymmetrizer
-pa[0].weight=1;
+NRVec<INDEXNAME> antinames;
-pa[0].perm.resize(rank());
+NRVec<NRVec_from1<int> > antigroups;
-pa[0].perm.identity();
+parse_antisymmetrizer(antisymmetrizer,antigroups,antinames);
+
+//check the names make sense and fill in the possibly missing ones as separate group
+if(antinames.size()>tmp.names.size()) laerror("too many indices in the antisymmetrizet");
+bitvector isexplicit(tmp.names.size());
+isexplicit.clear();
+for(int i=0; i<antinames.size(); ++i)
+	{
+	for(int j=0; j<i; ++j) if(antinames[i]==antinames[j]) laerror("repeated names in the antisymmetrizer");
+	for(int j=0; j<tmp.names.size(); ++j) 
+            if(antinames[i]==tmp.names[j]) 
+		{
+		isexplicit.set(j);
+		goto namefound;
+		}
+	laerror("index name from the antisymmetrized not found in tensor");
+	namefound: ;
+	}
+if(isexplicit.population()!=antinames.size()) laerror("internal error in add_permuted_contractions");
+
+//fill in additional names
+if(antinames.size()<tmp.names.size())
+	{
+	int lastgroup=antigroups.size()-1;
+	int lastclass;
+	if(antigroups.size()==0) lastclass=0;
+	else lastclass=antigroups[antigroups.size()-1][antigroups[antigroups.size()-1].size()];
+	int lastname=antinames.size()-1;
+	antinames.resize(tmp.names.size(),true);
+	antigroups.resize(antigroups.size()+tmp.names.size()-antinames.size(),true);
+	for(int j=0; j<names.size(); ++j)
+		if(!isexplicit[j])
+			{
+			++lastclass;
+			++lastname;
+			++lastgroup;
+			antigroups[lastgroup].resize(1);
+			antigroups[lastgroup][0]=lastclass;
+			antinames[lastname] = tmp.names[j];
+			}
+	}
+std::cout <<"final antisymmmetrizer names and groups"<<antinames<<antigroups;
+
+//prepare the antisymmetrizer
+PermutationAlgebra<int,int> pa = general_antisymmetrizer(antigroups,-2,true);
+
+//find permutation between antisym and TMP index order 
+NRPerm<int> antiperm=find_indexperm(antinames,tmp.names);
+
+//@@@conjugate the PA by antiperm or its inverse
+
+//@@@recast the PA
+
+//@@@permutationalgebra::is_identity() a pokd ano, neaplikovat ale vratit tmp tenzor resp. s nim udelat axpy
 
 //std::cout <<"LHS names = "<<names<<std::endl;
 //std::cout <<"TMP names = "<<tmp.names<<std::endl;
 
 //find permutation which will bring indices of tmp to the order as in *this: this->names[i] == tmp.names[p[i]]
-NRPerm<int> basicperm(rank());
+NRPerm<int> basicperm=find_indexperm(names,tmp.names);
-basicperm.clear();
-for(int i=0; i<rank(); ++i) 
-	{
-	for(int j=0; j<tmp.rank(); ++j)
-             if((* const_cast<const NRVec<INDEXNAME> *>(&names))[i] == tmp.names[j])
-		{
-		basicperm[i+1]=j+1;
-		break;
-		}
-	}
-if(!basicperm.is_valid()) laerror("indices mismatch between lhs and rhs in add_permuted_contractions");
 std::cout <<"Basic permutation = "<<basicperm<<std::endl;
 
+//@@@probably only onw of those will work
 PermutationAlgebra<int,T> pb = basicperm*pa;
 apply_permutation_algebra(tmp,pb,false,(T)1/(T)pb.size(),beta);
-//equivalently possible 
+//equivalently possible (for trivial pa)
 //PermutationAlgebra<int,T> pb = basicperm.inverse()*pa;
 //apply_permutation_algebra(tmp,pb,true,(T)1/(T)pb.size(),beta);
 }