From 12cf5b76a5bdd6b59b90acd77a29ac6b164ed2ef Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Wed, 22 Oct 2025 16:50:02 +0200 Subject: [PATCH] bugfix unwind_index; Tucker works --- t.cc | 83 +++++++++++++++++++++++++-- tensor.cc | 167 +++++++++++++++++++++++++++++++++++++----------------- tensor.h | 17 +++++- 3 files changed, 206 insertions(+), 61 deletions(-) diff --git a/t.cc b/t.cc index 3ecd74c..fa8e767 100644 --- a/t.cc +++ b/t.cc @@ -3275,13 +3275,20 @@ if(0) int n=5; INDEXGROUP g; g.number=4; -g.symmetry= -1; +g.symmetry= 0; +//g.symmetry= 1; +//g.symmetry= -1; g.offset=0; g.range=n; Tensor e(g); e.randomize(1.); Tensor eu = e.unwind_index(0,1); +Tensor eu2 = e.unwind_index(0,2); +Tensor eu3 = e.unwind_index(0,3); + +cout < g(4); +for(int i=0; i<4; ++i) + { + g[i].number=1; + g[i].symmetry= 0; + g[i].offset=0; + g[i].range=n+i; + } + +Tensor e(g); +e.randomize(1.); + +cout<< "E shape = "< eu = e.unwind_index(1,0); +cout<< "Eu shape = "< eu2 = e.unwind_index(2,0); +cout<< "Eu2 shape = "< eu3 = e.unwind_index(3,0); +cout<< "Eu3 shape = "< a; @@ -3576,5 +3619,33 @@ cout << "Error "<<(u*sdiag*vt-abak).norm()<>r>>n; +NRVec shape(r); +for(int i=0; i x(shape); +x.randomize(1.); +cout< x0(x); +x0.copyonwrite(); +bool inv=true; +NRVec > dec=x.Tucker(1e-12,inv); +cout<<"Tucker\n"< y = x.inverseTucker(dec,inv); +cout <<"invTucker\n"<number<=0) laerror("empty index group"); //we do not support scalar as a trivial case r+=sh->number; } myrank=r; @@ -46,6 +46,7 @@ return r; template LA_largeindex Tensor::calcsize() { +if(shape.size()==0) laerror("tensor must have rank at least 1"); groupsizes.resize(shape.size()); cumsizes.resize(shape.size()); LA_largeindex s=1; @@ -330,11 +331,11 @@ calcsize(); } template -Tensor::Tensor(const NRMat &x) +Tensor::Tensor(const NRMat &x, bool flat) : data(&x(0,0),x.nrows()*x.ncols()) { myrank=2; -if(x.nrows()==x.ncols()) +if(x.nrows()==x.ncols() && !flat) { shape.resize(1); shape[0].number=2; @@ -542,7 +543,9 @@ help_t->data[target] = *v; template Tensor Tensor::permute_index_groups(const NRPerm &p) const { -NRVec newshape=shape.permuted(p); +//std::cout <<"permute_index_groups permutation = "< newshape=shape.permuted(p,true); +//std::cout <<"permute_index_groups newshape = "< r(newshape); //prepare statics for the callback @@ -568,13 +571,8 @@ for(int i=0; i &shape) -{ -int ii=0; -for(int g=0; g &shape) +{return flatposition(i.group,i.index,shape);} int flatposition(int group, int index, const NRVec &shape) { @@ -584,12 +582,26 @@ ii += index; return ii; } +INDEX indexposition(int flatindex, const NRVec &shape) +{ +INDEX I={0,0}; +if(flatindex<0) laerror("illegal index in indexposition"); +for(int g=0; g static void unwind_callback(const SUPERINDEX &I, T *v) { FLATINDEX J = superindex2flat(I); -FLATINDEX JP = J.permuted(*help_p,true); +FLATINDEX JP = J.permuted(*help_p,false); //std::cout <<"TEST unwind_callback: from "<)(JP); //rhs operator() generates the redundant elements for the unwinded lhs tensor } @@ -637,6 +649,8 @@ for(int i=0; i u = rhs1.unwind_index(group,index); -if(conjugate1) u.conjugateme(); -Tensor rhsu = rhs.unwind_index(rhsgroup,rhsindex); +const Tensor u = conjugate1? (rhs1.unwind_index(group,index)).conjugate() : rhs1.unwind_index(group,index); +const Tensor rhsu = rhs.unwind_index(rhsgroup,rhsindex); NRVec newshape(u.shape.size()+rhsu.shape.size()-2); @@ -1077,51 +1092,99 @@ return r; template -NRVec > Tensor::Tucker(typename LA_traits::normtype thr) +NRVec > Tensor::Tucker(typename LA_traits::normtype thr, bool inverseorder) { int r=rank(); NRVec > ret(r); -if(r<2) return ret; +if(r<1) laerror("illegal rank in Tucker"); +copyonwrite(); -int rr=0; -for(int i=0; i(data,data.size(),1); + shape[0].range=1; + data.resize(calcsize()); + calcrank(); + data[0]=1; + return ret; + } + +//loop over all indices; relies on the fact tha unwinding does not change order of remaining indices +for(int i=0; i um; + NRVec ushape; + { + Tensor u=unwind_index(I); + ushape=u.shape; ushape.copyonwrite(); + um=u.matrix(); + } + int mini=um.nrows(); if(um.ncols() u(um.nrows(),mini),vt(mini,um.ncols()); + NRVec::normtype> w(mini); + singular_decomposition(um,&u,w,&vt,0); + um.resize(0,0); //deallocate + int preserve=mini; + for(int k=0; k umnew; + if(preserve um; - NRVec ushape; - { - Tensor u=unwind_index(i,j); - ushape=u.shape; - um=u.matrix(); - } - int mini=um.nrows(); if(um.ncols() u(um.nrows(),mini),vt(mini,um.ncols()); - NRVec::normtype> w(mini); - singular_decomposition(um,&u,w,&vt,0); - um.resize(0,0); //deallocate - int preserve=mini; - for(int k=0; k umnew; - if(preserve newdata(umnew); - *this = Tensor(ushape,newdata); + vt=vt.submatrix(0,preserve-1,0,um.ncols()-1); + w=w.subvector(0,preserve-1); + umnew=u.submatrix(0,um.nrows()-1,0,preserve-1); } + else umnew=u; + ret[(inverseorder? r-i-1 : i)]=vt.transpose(true); + umnew.diagmultr(w); + //rebuild tensor of the preserved shape from matrix + ushape[0].range=preserve; + { + NRVec newdata(umnew); + umnew.resize(0,0);//deallocate + *this = Tensor(ushape,newdata); + } + } +if(!is_flat()) laerror("this should not happen"); +if(!inverseorder) + { + NRPerm p(r); + for(int i=1; i<=r; ++i) p[r-i+1]=i; + *this = permute_index_groups(p); + } return ret; } +template +Tensor Tensor::inverseTucker(const NRVec > &x, bool inverseorder) const +{ +if(rank()!=x.size()) laerror("input of inverseTucker does not match rank"); +Tensor tmp(*this); +Tensor r; +if(!is_flat()) laerror("inverseTucker only for flat tensors as produced by Tucker"); +for(int i=0; i mat(x[i],true); + r= tmp.contraction(i,0,mat,0,0,(T)1,false,false); + if(i p(r.rank()); + for(int i=1; i<=r.rank(); ++i) p[r.rank()-i+1]=i; + return r.permute_index_groups(p); + } +else + return r; +} + diff --git a/tensor.h b/tensor.h index c69bc43..87c8dd0 100644 --- a/tensor.h +++ b/tensor.h @@ -35,6 +35,7 @@ #include "vec.h" #include "mat.h" #include "smat.h" +#include "fourindex.h" #include "miscfunc.h" //TODO: @@ -44,6 +45,10 @@ //@@@ will need to store vector of INDEX to the original tensor for the result's flatindex //@@@ will not be particularly efficient // +//@@@conversions to/from fourindex +// +//@@@!!!!!!!!!!!const loopover and grouploopover +// //@@@?general permutation of individual indices - check the indices in sym groups remain adjacent, calculate result's shape, loopover the result and permute using unwind_callback // // @@ -122,7 +127,9 @@ int index; }; typedef NRVec INDEXLIST; //collection of several indices +int flatposition(int group, int index, const NRVec &shape); int flatposition(const INDEX &i, const NRVec &shape); //position of that index in FLATINDEX +INDEX indexposition(int flatindex, const NRVec &shape); //inverse to flatposition FLATINDEX superindex2flat(const SUPERINDEX &I); @@ -150,7 +157,7 @@ public: Tensor(const Tensor &rhs): myrank(rhs.myrank), shape(rhs.shape), groupsizes(rhs.groupsizes), cumsizes(rhs.cumsizes), data(rhs.data) {}; Tensor(int xrank, const NRVec &xshape, const NRVec &xgroupsizes, const NRVec xcumsizes, const NRVec &xdata) : myrank(xrank), shape(xshape), groupsizes(xgroupsizes), cumsizes(xcumsizes), data(xdata) {}; explicit Tensor(const NRVec &x); - explicit Tensor(const NRMat &x); + explicit Tensor(const NRMat &x, bool flat=false); explicit Tensor(const NRSMat &x); NRMat matrix() const {return NRMat(data,data.size()/groupsizes[0],groupsizes[0],0);}; //reinterpret as matrix with column index being the tensor's leftmost index group (typically the unwound single index) @@ -162,6 +169,7 @@ public: LA_largeindex size() const {return data.size();}; void copyonwrite() {shape.copyonwrite(); groupsizes.copyonwrite(); cumsizes.copyonwrite(); data.copyonwrite();}; void resize(const NRVec &s) {shape=s; data.resize(calcsize()); calcrank();}; + void deallocate() {data.resize(0); shape.resize(0); groupsizes.resize(0); cumsizes.resize(0);}; inline Signedpointer lhs(const SUPERINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer(&data[i],sign);}; inline T operator()(const SUPERINDEX &I) const {int sign; LA_largeindex i=index(&sign,I); if(sign==0) return 0; return sign>0 ?data[i] : -data[i];}; inline Signedpointer lhs(const FLATINDEX &I) {int sign; LA_largeindex i=index(&sign,I); return Signedpointer(&data[i],sign);}; @@ -176,6 +184,8 @@ public: inline Tensor& operator/=(const T &a) {data/=a; return *this;}; inline Tensor operator/(const T &a) const {Tensor r(*this); r /=a; return r;}; + typename LA_traits::normtype norm() const {return data.norm();}; + inline Tensor operator*(const Tensor &rhs) const {return Tensor(rhs.shape.concat(shape),data.otimes2vec(rhs.data));} //outer product, rhs indices will be the less significant Tensor& conjugateme() {data.conjugateme(); return *this;}; @@ -215,6 +225,7 @@ public: Tensor permute_index_groups(const NRPerm &p) const; //rearrange the tensor storage permuting index groups as a whole Tensor unwind_index(int group, int index) const; //separate an index from a group and expand it to full range as the least significant one (the leftmost one) + Tensor unwind_index(const INDEX &I) const {return unwind_index(I.group,I.index);}; Tensor unwind_indices(const INDEXLIST &il) const; //the same for a list of indices void addcontraction(const Tensor &rhs1, int group, int index, const Tensor &rhs2, int rhsgroup, int rhsindex, T alpha=1, T beta=1, bool doresize=false, bool conjugate1=false, bool conjugate=false); //rhs1 will have more significant non-contracted indices in the result than rhs2 inline Tensor contraction(int group, int index, const Tensor &rhs, int rhsgroup, int rhsindex, T alpha=1, bool conjugate1=false, bool conjugate=false) const {Tensor r; r.addcontraction(*this,group,index,rhs,rhsgroup,rhsindex,alpha,0,true, conjugate1, conjugate); return r; }; @@ -232,8 +243,8 @@ public: void split_index_group(int group); //formal split of a non-symmetric index group WITHOUT the need for data reorganization void merge_adjacent_index_groups(int groupfrom, int groupto); //formal merge of non-symmetric index groups WITHOUT the need for data reorganization Tensor merge_index_groups(const NRVec &groups) const; - NRVec > Tucker(typename LA_traits::normtype thr=1e-12); //HOSVD-Tucker decomposition, return core tensor in *this, flattened - + NRVec > Tucker(typename LA_traits::normtype thr=1e-12, bool inverseorder=true); //HOSVD-Tucker decomposition, return core tensor in *this, flattened + Tensor inverseTucker(const NRVec > &x, bool inverseorder=true) const; //rebuild the original tensor from Tucker };