Routines from first three notebooks

Reordering Indices

NamedTranspose[mat_,orig_,final_]:=Module[{ord=Flatten[Position[final,#]&/@orig]},Transpose[mat,ord]]

NamedTranspose::usage="NamedTranspose[mat_,orig_,final_] takes a tensor `mat' whose indices are named `orig' and shuffles the order, so they match `final'.  For example NamedTranspose[SparseArray[{1,2,3}->5,{3,3,3}],{'a','b','c'},{'b','c','a'}]";

ReOrder[tensordescription_List,finalnames_]:=Module[{
tensor="tensor"/.tensordescription, (* Extract Tensor *)
names="names"/.tensordescription, (* Extract names of indices *)
newtensor (* Declare local variables *)
},
(* Calculate new Tensor *)
newtensor=NamedTranspose[tensor,names,finalnames];
(* Assemble into rules *)
{"tensor"->newtensor,"names"->finalnames}
]

ReOrder::usage="ReOrder[tensordescription_List,finalnames_] takes a NamedTensor object of the form {'tensor'->tensor,'names'->names} and returns a NamedTensor, with the indices permuted to match finalnames"

Combining Indices

CombineIndices[namedtensor_List,indicestocombine_List,newindex_]:=
Module[{oldnames="names"/.namedtensor,names,targpos,newnamedtensor,dims,tensor,newdims,newtensor,finalnamedtensor,ind,indexkeys,refine,newnames,targnames},
(* Shift the indices we want to the right *)
targpos=First@Flatten@Position[oldnames,Alternatives@@indicestocombine];
names=Join[Complement[oldnames,indicestocombine],indicestocombine];
newnamedtensor=ReOrder[namedtensor,names];
(* Next extract the tensor *)
tensor="tensor"/.newnamedtensor;
(* Calculate the lengths of each of the indices *)
dims=Dimensions[tensor];
newdims=Append[dims[[;;-Length[indicestocombine]-1]],Times@@dims[[-Length[indicestocombine];;]]];
(* Do the reshaping *)
newtensor=ArrayReshape[tensor,newdims];
newnames=Append[names[[;;-Length[indicestocombine]-1]],newindex];
(* Do a final reordering *)
targnames=Join[newnames[[;;targpos-1]],{newnames[[-1]]},newnames[[targpos;;-2]]];
ReOrder[{"tensor"->newtensor,"names"->newnames},targnames]]

rightCombineIndices[namedtensor_List,indicestocombine_List,newindex_]:=
Module[{oldnames="names"/.namedtensor,names,newnamedtensor,dims,tensor,newdims,newtensor,finalnamedtensor,ind,indexkeys,refine,newnames,targnames},
(* Shift the indices we want to the right *)
names=Join[SortBy[Complement[oldnames,indicestocombine],Position[oldnames,#]&],indicestocombine];
newnamedtensor=ReOrder[namedtensor,names];
(* Next extract the tensor *)
tensor="tensor"/.newnamedtensor;
(* Calculate the lengths of each of the indices *)
dims=Dimensions[tensor];
newdims=Append[dims[[;;-Length[indicestocombine]-1]],Times@@dims[[-Length[indicestocombine];;]]];
(* Do the reshaping *)
newtensor=ArrayReshape[tensor,newdims];
newnames=Append[names[[;;-Length[indicestocombine]-1]],newindex];
{"tensor"->newtensor,"names"->newnames}]

leftCombineIndices[namedtensor_List,indicestocombine_List,newindex_]:=
Module[{oldnames="names"/.namedtensor,names,newnamedtensor,dims,tensor,newdims,newtensor,finalnamedtensor,ind,indexkeys,refine,newnames,targnames},
(* Shift the indices we want to the right *)
names=Join[indicestocombine,SortBy[Complement[oldnames,indicestocombine],Position[oldnames,#]&]];
newnamedtensor=ReOrder[namedtensor,names];
(* Next extract the tensor *)
tensor="tensor"/.newnamedtensor;
(* Calculate the lengths of each of the indices *)
dims=Dimensions[tensor];
newdims=Prepend[dims[[Length[indicestocombine]+1;;]],Times@@dims[[;;Length[indicestocombine]]]];
(* Do the reshaping *)
newtensor=ArrayReshape[tensor,newdims];
newnames=Prepend[names[[Length[indicestocombine]+1;;]],newindex];
{"tensor"->newtensor,"names"->newnames}]

makename[indicestocombine_List]:=StringJoin@@Riffle[ToString/@indicestocombine,Table["_",{Length[indicestocombine]-1}]]

breakname[compoundname_String]:=StringSplit[compoundname,"_"]

CombineIndices[namedtensor_List,indicestocombine_List]:=
Module[{newindex=makename[indicestocombine]},
CombineIndices[namedtensor,indicestocombine,newindex]]

leftCombineIndices[namedtensor_List,indicestocombine_List]:=
Module[{newindex=makename[indicestocombine]},
leftCombineIndices[namedtensor,indicestocombine,newindex]]

rightCombineIndices[namedtensor_List,indicestocombine_List]:=
Module[{newindex=makename[indicestocombine]},
rightCombineIndices[namedtensor,indicestocombine,newindex]]

BreakIndices[namedtensor_List,indextobreak_,newnames_List,newdimensions_List]:=
Module[{tensor="tensor"/.namedtensor,
names="names"/.namedtensor,
olddims,
pos,newtensor,newdims,finalnames,dimtobereplaced
},
olddims=Dimensions[tensor];
pos=First@Flatten@Position[names,indextobreak];
(* Do a quick consistency check *)
dimtobereplaced=olddims[[pos]];
If[dimtobereplaced≠Times@@newdimensions,
Return[Failure["Inconsistent Dimensions",
Association["MessageTemplate"->"Trying to break an index of length `l1` into indices of length `dims`",
"MessageParameters"->Association[
"l1"->dimtobereplaced,
"dims"->newdimensions
]]]];
];
newdims=Join[olddims[[;;pos-1]],newdimensions,olddims[[pos+1;;]]];
newtensor=ArrayReshape[tensor,newdims];
finalnames=Join[names[[;;pos-1]],newnames,names[[pos+1;;]]];
{"tensor"->newtensor,"names"->finalnames}
]

Outer Products

disambiguate[names_List]:=Module[{seive,count},
count[_]=0;
seive[a_]:=Module[{astr=ToString[a],root,numticks=0},
root=astr;
While[StringTake[root,-1]=="'",
numticks+=1;
root=StringTake[root,{1,-2}];
];
count[root]+=1;
root<>StringRepeat["'",count[root]-1]];
seive/@names]

Combine[namedtensors__List]:={"tensor"->TensorProduct@@("tensor"/.{namedtensors}),"names"->disambiguate[Join@@("names"/.{namedtensors})]}

Contractions (Inner Products)

To Contract two tensors, the most natural way is to use “Dot” -- one puts the indices that you want contracted at the beginning and end of the two tensor.

Contract[namedtensor1_List,namedtensor2_List,indextocontract_]:=
Module[{t1, 
t2,
names1="names"/.namedtensor1,
names2="names"/.namedtensor2,
newnames1,newnames2,newnamedtensor1,newnamedtensor2,newt,newnames},
(*Get new ordering of names*)
newnames1=RotateLeft[names1,Flatten@Position[names1,indextocontract]];
newnames2=RotateLeft[names2,Flatten@Position[names2,indextocontract]-1];
(*Reorder tensor indices*)
newnamedtensor1=ReOrder[namedtensor1,newnames1];
newnamedtensor2=ReOrder[namedtensor2,newnames2];
(*Extract Tensors*)
t1="tensor"/.newnamedtensor1;
t2="tensor"/.newnamedtensor2;
(*Do Dot Product *)
newt=Dot[t1,t2];
newnames=Join[newnames1[[;;-2]],newnames2[[2;;]]];
{"tensor"->newt,"names"->newnames}
]

Contract[namedtensor_List,indicestocontract_List]:=
Module[{tensor="tensor"/.namedtensor,
names="names"/.namedtensor,
positions,newtensor,keep,newnames},
(* Find out which indices we need to contract *)
positions=Flatten[Position[names,#]&/@indicestocontract];
(* Do contraction *)
newtensor=TensorContract[tensor,positions];
(* Find out names of indices we are keeping *)
keep=Complement[Range[Length[names]],positions];
newnames=names[[keep]];
(* Reassemble *)
{"tensor"->newtensor,"names"->newnames}
]

Fancier contractions

autoContract[namedtensor_List]:=
Module[{tensor="tensor"/.namedtensor,newtensor,
names="names"/.namedtensor,signature,repeats,newindices},
signature=Last@Reap[MapIndexed[Sow[#2,#1]&,names],_,Rule];
repeats=Cases[signature[[All,2]],(locs:{_,_}):>Flatten[locs]];
newtensor=TensorContract[tensor,repeats];
newindices=Delete[names,Transpose[{Flatten[repeats]}]];
{"tensor"-> newtensor,"names"->newindices}
]

autoContract[namedtensor1_List,namedtensor2_List]:=
Module[{names1="names"/.namedtensor1,
names2="names"/.namedtensor2,commonnames,left,right,leftnames,rightnames,
ltens,rtens,newtensor,newnames,pos1,pos2},
commonnames=Intersection[names1,names2];
left=rightCombineIndices[namedtensor1,commonnames];
right=leftCombineIndices[namedtensor2,commonnames];
ltens="tensor"/.left;
rtens="tensor"/.right;
newtensor=Dot[ltens,rtens];
leftnames="names"/.left;
rightnames="names"/.right;
newnames=Join[leftnames[[;;-2]],rightnames[[2;;]]];
{"tensor"->newtensor,"names"->newnames}]

autoContract[namedtensor1_List,namedtensor2_List,moretensors__List]:=
autoContract[autoContract[namedtensor1,namedtensor2],moretensors]

Adding Tensors

addTensor[namedtensor1_List,namedtensor2_List]:=
Module[{tens1="tensor"/.namedtensor1,
tens2="tensor"/.namedtensor2,
names1="names"/.namedtensor1,
names2="names"/.namedtensor2,ordnamedtens2},
If[names1==names2,
Return[{"tensor"->tens1+tens2,"names"->names1}]];
If[Sort[names1]==Sort[names2],
ordnamedtens2=ReOrder[namedtensor2,names1];
Return[{"tensor"->tens1+("tensor"/.ordnamedtens2),"names"->names1}]];
Failure["Shape Mismatch",Association["MessageTemplate"->
"Trying to add a tensor with indices `ind1` to one with indices `ind2`",
"MessageParameters"->Association["ind1"->names1,"ind2"->names2]]]
]

Applying operations to tensors

ApplytoTensor[namedtensor_List,op_]:=namedtensor /.("tensor"->t_):>("tensor"->op[t])

Multiplying a tensor by a number

multiply[namedtensor_List,number_]:=namedtensor /.("tensor"->t_):>("tensor"->number t)

Useful matrices

pauliZ=SparseArray@{{1,0},{0,-1}};
pauliX=SparseArray@{{0,1},{1,0}};
pauliY=SparseArray@{{0,-I},{I,0}};
pauliI = SparseArray@{{1,0},{0,1}};

id[leftname_,rightname_,length_Integer]:={"tensor"->SparseArray[{i_,i_}->1,{length,length}],"names"->{leftname,rightname}}

SVD

svd[namedtensor_,newbonds_:{"b1","b2"},trunc___]:=
Module[{tensor="tensor"/.namedtensor,leftbond,rightbond,csvd,maxkeep,tol,u,d,dvals,v,numsv,b1=newbonds[[1]],b2=newbonds[[-1]]},
{leftbond,rightbond}="names"/.namedtensor;
maxkeep="trunc"/.{trunc} //."trunc"->Min[Dimensions[tensor]];
tol="tolerance"/.{trunc} //."tolerance"->0;
csvd=SingularValueDecomposition[tensor,UpTo[maxkeep],Tolerance->tol];
(* I find the SVD doesn't always truncate properly, so lets do that *)
dvals=Take[Select[Tr[csvd[[2]],List],(#>tol&)],UpTo[maxkeep]];
numsv=Length[dvals];
u={"tensor"->csvd[[1,All,;;numsv]],"names"->{leftbond, b1}};
d={"tensor"->SparseArray[{Band[{1,1}]->dvals}],"names"->{b1,b2}};
v={"tensor"->csvd[[3,All,;;numsv]],"names"->{ rightbond,b2}};
{u,d,v}]

Eigensolver

eig[nTensor_List,λ_]:=
Module[{mat="tensor"/.nTensor,names="names"/.nTensor,eigsys,ord},
If[Length[mat]>15,
eigsys=Eigensystem[mat,1,Method->{"Arnoldi","Shift"->λ}],
(* if the matrix is small enough, just do a dense routine *)
eigsys=Eigensystem[Normal[mat]];
ord=Ordering[eigsys[[1]]];
eigsys=eigsys[[All,{ord[[1]]}]];
];
{"tensor"->eigsys[[2,1]],"names"->{names[[-1]]},"eigenvalue"->eigsys[[1,1]]}]

eig[nTensor_List,λ_,guess_List]:=
Module[{mat="tensor"/.nTensor,
vec="tensor"/.guess
,names="names"/.nTensor,eigsys,ord},
If[Length[mat]>2,
eigsys=Eigensystem[mat,1,Method->{"Arnoldi","Shift"->λ,"StartingVector"->vec}],
(* if the matrix is small enough, just use dense matrix techniques *)
eigsys=Eigensystem[Normal[mat]];
ord=Ordering[eigsys[[1]]];
eigsys=eigsys[[All,{ord[[1]]}]];
];
{"tensor"->eigsys[[2,1]],"names"->{names[[-1]]},"eigenvalue"->eigsys[[1,1]]}]