(************ Basic Pauli Spin Matrices ********************)
σ[0] = {{1,0}, {0,1}};
σ[1] = {{0,1}, {1,0}};
σ[2] = {{0, -I}, {I,0}};
σ[3] = {{1,0}, {0,-1}};

(******** This defines the number of spins (n), and the number of operators,
we are allowed to insert, (k) ***************)
n = 5;
k = 1;

(********* Counts the number of constraints we expect to get **********)
constraintcount[n_, k_] := Sum[Binomial[n,j] 3^j, {j, 0, k}]

constraintcount[4,1]
2^4

(****** Now, define the spin matrices acting at each site, given by s[i,a] ******)
For[i = 1, i ≤ n, i++, For[a = 1, a ≤ 3, a++, list = Table[σ[0], {n}]; list[[i]] = σ[a]; s[i,a] = KroneckerProduct@@list ]]
idn = KroneckerProduct@@Table[σ[0], {n}];
s[x_, 0] := idn;

(********* Now construct all possible products of up to K insertions. We use
a slightly inefficient algorithm, but this is done only once*)

Clear[a]

makeproduct[ curcount_,  pos_, curop_] := Module[{a}, If[(curcount == 0) && (pos == 1), (*Initialize*)
opcount = 0;
listops = {}; ];
If[(curcount ==  k)  || (pos > n), 
Print[opcount]; (*Comment this line to stop the printouts *)
opcount++; listops = Append[listops, curop],(*Else continue*)makeproduct[ curcount, pos+1, curop];
For[a=1, a ≤  3, a++, 
(*Not the identity. Multiply the s, and pass it ahead *)
makeproduct[ curcount+1, pos+1, s[pos,a].curop]
]]]

Print["The counter below must reach ", constraintcount[n,k]]
Print["*****************"]
makeproduct[0,1,idn]
Print["*****************"]

(********************** Now Generate a Random Vector ************)

psi = Array[RandomComplex[{-1 - I, 1 + I}]&, 2^n];
psi = psi/Sqrt[psi.Conjugate[psi]];
MatrixForm[psi]
(*Check that psi is correctly normalized *)
Chop[psi.Conjugate[psi]];

(****** Now generate the vectors of Hpsi ***********)
v = Range[opcount];
For[i = 1, i≤ opcount, i++, v[[i]] = listops[[i]].psi];
(***** Now, generate the metric on Hpsi *************)
g = Array[0&, {opcount, opcount}];
For[i = 1, i ≤ opcount, i++, For[j = 1, j ≤ opcount, j++, g[[i,j]] = Conjugate[v[[i]]].v[[j]]]]
(****** Now, find the inverse metric *******)
ginv = Inverse[g];

(********* Now, to find the tilde of a given element. Finding S, lets us find all these at the same time, but we have to go one higher in k. So, here, since increasing k is expensive, we find the tildes one by one *******)
(************************ TO FIND THE TILDES FOR ANOTHER OPERATOR, START HERE. THE STEPS ABOVE THIS LINE DO NOT NEED TO BE REPEATED *************)
tsi = 1;
tsa = 1;
(****** Now generate the right hand side of the linear equations ***********)
vans = Range[opcount];
For[i = 1, i≤ opcount, i++, vans[[i]] = -listops[[i]].s[tsi, tsa].psi];
(******** Now, calculate ts[tsi, tsa] *****)
vdual = ginv.Conjugate[v];
ts[tsi, tsa] = Sum[KroneckerProduct[vans[[i]], vdual[[i]]], {i, 1, opcount}];

MatrixForm[s[tsi,tsa]]

(************ Now, we perform a few verification checks *********)
MatrixForm[ts[tsi,tsa]]

(********** Check the Correlation Function **********)
Chop[Conjugate[psi].ts[tsi,tsa].s[1, 1].psi]

(****** Check that it commutes with ordinary Pauli matrices, while acting on |\Psi> ***********)
arbop = s[5,2];
Chop[(arbop.ts[tsi, tsa] - ts[tsi, tsa].arbop).psi]

(********* Note, the commutator of these guys may be quite complicated *****)


comm = arbop.ts[tsi, tsa] - ts[tsi, tsa].arbop;
MatrixForm[comm]

(*remove the semi-colon to see the complicated commutator*)
Chop[comm.psi]

(***** But the commutator annihilates the state *******)