Value of a function over permutations of two lists

Question

I want to do the sum (in Ising lattice gauge theory)

n = 3;
energy[sx_,sy_] := Sum[sx[[i, j]] sy[[i, j]] sx[[i + 1, j]] sy[[i, j + 1]], {i, 1, n}, {j, 1, n}]
sx = Tuples[{-1,1}, {n+1,n}];
sy = Tuples[{-1,1}, {n,n+1}];

I want the values of this function over all possible permutations of sx and sy but I am unable to find a way to do so. Basically this is a variant of the usual Ising model but this time the spins placed on the bonds instead of the nodes. Spins still take values $\pm 1$ but the Hamiltonian is now the product of spins on a 'square'. I am trying to calculate the partition function of this model and hence this sum arose. I divided the spins into ones along the $x$ axis and others along the $y$ axis and hence am using two lists. Any help would be greatly appreciated. Thanks.

Edit: as pointed out in the comment the definitions of sx and sy are wrong. sx is the name of the matrix of horizontal bonds and sy is the matrix of vertical bonds.

Answer 1

You can do the multiplications more efficiently without needing Sum.

Remove["Global`*"];
energy[sx_, sy_] := 
 Total[Most[sx]*sy[[All, ;; -2]]*Rest[sx]*sy[[All, 2 ;;]], 2]

SeedRandom[123];
sx = RandomChoice[{-1, 1}, {4, 3}];
sy = RandomChoice[{-1, 1}, {3, 4}];
energy[sx, sy]
(* 5 *)

Summing over all permutations is easy, but I wouldn't make n much bigger due to memory constraints:

Remove["Global`*"];

energy[sx_, sy_] := 
 Total[Most[sx]*sy[[All, ;; -2]]*Rest[sx]*sy[[All, 2 ;;]], 2]
n = 3;
sxPerms = Tuples[{-1, 1}, n*(n + 1)];
syPerms = Tuples[{-1, 1}, n*(n + 1)];
result = Table[
   energy[ArrayReshape[sx, {n + 1, n}], ArrayReshape[sy, {n, n + 1}]]
   , {sx, sxPerms}, {sy, syPerms}];
MatrixPlot[result]

This image looks cool, but it's not very meaningful. The row/col indices correspond to permutations.

Instead we can compute the grid products and preserve them, then sum later to add all possible grids together. The result is a zero matrix:

energy[sx_, sy_] := 
 Most[sx]*sy[[All, ;; -2]]*Rest[sx]*sy[[All, 2 ;;]]
n = 3;
sxPerms = Tuples[{-1, 1}, n*(n + 1)];
syPerms = Tuples[{-1, 1}, n*(n + 1)];
result = 
  Sum[energy[ArrayReshape[sx, {n + 1, n}], 
    ArrayReshape[sy, {n, n + 1}]], {sx, sxPerms}, {sy, syPerms}];

(* {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}} *)

---

Also, you can avoid needing to store the tuples by generating them based off an index. This is slower, but uses less memory:

Remove["Global`*"];

energy[sx_, sy_] := 
 Total[Most[sx]*sy[[All, ;; -2]]*Rest[sx]*sy[[All, 2 ;;]], 2]
generateMatrix[i_, dim_] := With[{d = Times @@ dim},
  ArrayReshape[
   2 IntegerDigits[i, 2, IntegerLength[2^d - 1, 2]] - 1
   , dim]]

n = 3;
result = 
  ParallelTable[
   energy[generateMatrix[ix, {n + 1, n}], 
    generateMatrix[iy, {n, n + 1}]]
   , {ix, 0, 2^(n*(n + 1)) - 1}, {iy, 0, 2^(n*(n + 1)) - 1}];