(* This provides a demonstration for the programs used for the project "Canonical tensor model through data analysis". The program below represents more or less the actual programs used in the project. This treats the spcial case of a fuzzy circle. *)

(* Construction of real sym 3-way tensor p for fuzzy circle *)

L=7; (* N=2L+1 *)

Print["Starting construction of a symmetric 3-way tensor for a fuzzy circle with N=",2L+1];

al=1;

pmom=Table[Exp[-al/L^2(i^2+j^2+k^2)]KroneckerDelta[i+j+k,0],{i,-L,L},{j,-L,L},{k,-L,L}];

maxL=2 L+1;

T=Table[0,{maxL},{maxL}];

Do[
T[[i,i]]=1/Sqrt[2];
T[[i,maxL-i+1]]=1/Sqrt[2];
,{i,1,L}];

T[[L+1,L+1]]=1;

Do[
T[[i,i]]=1/(Sqrt[2] I);
T[[i,maxL-i+1]]=-1/(Sqrt[2] I);
,{i,L+2,maxL}];

p=Chop[1.*T.Transpose[T.pmom.Transpose[T],{2,1,3}]];


(* 
Tensor-rank decomposition: Mathematica version. Paper used much faster C++version.

Input: p ; three-index real symmetric tensor p .
Output: vlist; list of N-dim vectors v^i (i=1,2,...,r=Rank).
  *)

Clear[minv,v,u,w,vert,f,r];

r= Floor[maxL]; (* The rank you assume*)
repnum=100; (* The number of searching cycles. The larger, the better result but slower *)

Print["Starting tensor-rank decomposition with Rank=",r];

minv=Compile[{{a,_Real,1},{b,_Real,3}},Module[{tor,countmax,count,prevw,w},

tor=.0000001; 
countmax=100;

count=0;

w=a;

While[
True,
count ++;
prevw=w;
w=(b.w).w/(w.w);
If[(prevw-w).(prevw-w)/w.w<tor || count>countmax,Break[] ]] ;

(w.w )^(-1/3) w]];


time=AbsoluteTime[];

workp=p;

w=Table[x[i],{i,1,maxL}];

plist={p};
vlist={};

Do[

v=minv[Table[RandomReal[{-1,1}],{maxL}],workp];

vlist=Append[vlist,v];

workp=workp-Outer[Times,v,v,v];

,{r}];

plist=Append[plist,workp];

countermax=Floor[repnum/100];
counter=0;

Do[

condlist=Table[vert=vlist[[i]];If[vert.(vert.p.vert)<0,0,1],{i,1,Length[vlist]}];

vlist=condlist*vlist;

workp=p-Sum[Outer[Times,vlist[[i]],vlist[[i]],vlist[[i]]],{i,1,r}];

Do[
workp=workp+Outer[Times,vlist[[vnum]],vlist[[vnum]],vlist[[vnum]]];

v=minv[vlist[[vnum]]+0.1*Table[RandomReal[{-1,1}],{maxL}],workp];

vlist[[vnum]]=v;

workp=workp-Outer[Times,v,v,v];

,{vnum,1,r}];

If[counter==countermax,
plist=Append[plist,workp];counter=0,
counter++];

,{repnum}];

Print["Ended with time=",AbsoluteTime[]-time, " sec"];

conmat=vlist.Transpose[vlist];

condlist=Table[vert.(vert.p.vert)<0,{i,1,Length[vlist]}];

vert=Table[TensorContract[plist[[i]].plist[[i]],{{1,3},{2,4}}],{i,1,Length[plist],Max[1,Floor[Length[plist]/10]]}];

Print["\n Development of error ΔP^2:"];
ListLogPlot[vert,PlotRange->All]

Print["Final error ratio Δp/p =",Sqrt[vert[[-1]]/vert[[1]]]]



(*Showing inner products among points *)

conmat=vlist.Transpose[vlist];

vert=Flatten[conmat];

Print["\n Histogram of v^i・v^j :"];
Histogram[vert,PlotRange->All]


(* 3D Image of the object made by connecing points satsifying v^i・v^j > cut *)

cut=.4; (* The parameter defining neighborhoods *)

cmat={};

r=Length[conmat];

Do[
If[conmat[[i,j]]>cut,cmat=Append[cmat,{i,j}]];,
{i,1,r-1},{j,i+1,r}];

adjmat=SparseArray[Table[cmat[[i]]->1,{i,1,Length[cmat]}],{r,r}];

adjmat=adjmat+Transpose[adjmat];

vert=ConnectedComponents[AdjacencyGraph[adjmat]];
llist=Length /@ vert;
lmax=Max[llist];
pos=Position[llist,lmax];

concomp=vert[[pos[[1,1]]]];

adjmat0=Table[adjmat[[concomp[[i]],concomp[[j]]]],{i,1,Length[concomp]},{j,1,Length[concomp]}];

size=Length[adjmat0];

vert=0;

pos=Table[z[++vert],{i,1,size},{j,1,3}];

totz=vert;

al=1;

potweight[i_,j_]:=Abs[conmat[[i,j]]];

pot=Sum[
Module[{r2},
r2=pos[[i]]-pos[[j]];
r2=Total[r2^2];
al/Sqrt[r2]+If[adjmat0[[i,j]]==1, potweight[i,j]r2,0]],{i,1,size-1},{j,i+1,size}]+
.01*Sum[pos[[i,3]]^2,{i,1,size}];

minz=FindMinimum[pot,Table[{z[i],RandomReal[{-10,10}]},{i,1,totz}],MaxIterations->100];

pos0=pos /. minz[[2]];

pos0=pos0-Table[Mean[pos0],{Length[pos0]}];

linelist={};

Do[
If[adjmat0[[i,j]]≠0,
AppendTo[linelist,{pos0[[i]],pos0[[j]]}]];
,{i,1,size-1},{j,i+1,size}];

Print["\n 3D image with cut=",cut," :"];

Graphics3D[{{Thick,Blue,Line[linelist]},{PointSize[Large],Red,Point[pos0]}},Boxed->False,ViewPoint->Above]