0
$\begingroup$

With the following code, I am trying to get Mathematica to return a 3×3 magic square where the constraint is not that the rows, columns, and diagonals sum to the same value, but rather that they satisfy a given function, memoF[x, y, z].

If memoF[x, y, z] = x + y + z, it works very quickly and without issues. However, as soon as the function becomes more complex —such as in the case I am showing, where functions like Max or KroneckerDelta are used— the computation times become unmanageable. I assume this is because FindInstance starts considering all possible combinations.

I suppose that maybe I should approach the problem in a different way, or perhaps code it in a style more suited to a language where I can have more control over each step... I'm stuck Any ideas?

 Clear[memoF];

   memoF[x_, y_, z_] := 
  memoF[x, y, z] = 
   Max[Max[x, y] + 1 + KroneckerDelta[x, y], z] + 1 + 
    KroneckerDelta[Max[x, y] + 1 + KroneckerDelta[x, y], z];

magicSquareILP[k_, target_] := 
 Module[{vars, constraints}, 
  vars = Flatten@Table[a[i, j], {i, 3}, {j, 3}];
  constraints = Join[
    
    Thread[1 <= vars <= k],
    
    Table[memoF @@ Table[a[i, j], {j, 3}] == target, {i, 3}],
    Table[memoF @@ Table[a[j, i], {j, 3}] == target, {i, 3}],
    {memoF @@ Diagonal@Array[a, {3, 3}] == target, 
     memoF @@ Diagonal[Reverse /@ Array[a, {3, 3}]] == target}
    ];
  
  FindInstance[constraints, vars, Integers, 1]]

With[{k = 4, target = 5}, magicSquareILP[k, target]]

Thank you!

Improve this question
$\endgroup$
1
  • 1
    $\begingroup$ Perhaps related: 73131. $\endgroup$ Commented Feb 16 at 12:10

2 Answers 2

Reset to default
1
$\begingroup$

Here is a way using the SelectPermutations function by Sander Huisman. I assume you want each row/column/diagonal to be duplicate free (otherwise you can use SelectTuples in place of SelectPermutations). The whole matrix itself cannot be duplicate free since the elements can only take on k values, and k < n^2.

(*SelectPermutations generates tuples meeting a condition,\
without generating all possible Permutations*)
selectPerms = ResourceFunction["SelectPermutations"];

(*function to get rows,columns,diagonal,\
antidiagonal of candidate square*)
matParts[mat_] := 
 Join @@ Comap[{Transpose, List@*Diagonal, 
    List@Reverse[Diagonal[Reverse /@ #]] &}, mat]

(*parameter defintions;
n=dim of magic square;
s=number of magic squares you want;
k=maximum value of integer elements of square;
c=function value of each row/column/diag;
f=function;*)

findMagicSquares[n_, s_, k_, c_, f_] := Module[{rowFunc, possElems},
  rowFunc := Fold[f, #] &;
  possElems = selectPerms[Range[k], {n}, rowFunc[#] == c &];
  selectPerms[possElems, {n}, Union[rowFunc /@ matParts[#]] == {c} &, 
   s]
  ]

It selects all possible tuples of length n of the integers 1,2,...,k that equal c under f. Then it finds all n-tuples of those tuples such that the columns and diagonals also equal c under f.

I redefine your memoF in terms of a pairwise function that gets Folded along each row/column/diagonal:

f[x_, y_] := 1 + Max[x, y] + Boole[x == y]

So with n = 3, s =2, k=4, c=5 for example we generate 2 magic squares that meet the criteria

f[x_, y_] := 1 + Max[x, y] + Boole[x == y]
magicList = findMagicSquares[3, 2, 4, 5, f];
MatrixForm /@ magicList

enter image description here

And confirming that these two matrices indeed satisfy your original memoF function, and the matrix rows/columns/diagonals are all equal to your target of 5:

memoF[x_, y_, z_] := 
  memoF[x, y, z] = 
   Max[Max[x, y] + 1 + KroneckerDelta[x, y], z] + 1 + 
    KroneckerDelta[Max[x, y] + 1 + KroneckerDelta[x, y], z];

memoF @@@ matParts[#] & /@ magicList

(*{{5, 5, 5, 5, 5}, {5, 5, 5, 5, 5}}*)

Ways to improve

You could iteravely add one row at a time to a matrix, and check that the partial columns/diagonals are the beginning of one of the possElems (i.e. if you have two rows put together, the partial columns/diagonals of length 2 should all be the first two elements of at least one element of possElems).

Improve this answer
$\endgroup$
1
  • $\begingroup$ Amazing. Thank you very much! $\endgroup$ Commented Feb 19 at 7:05
0
$\begingroup$

Assume all your constraints are declared in a list: "constraints" and the variables in a list: "vars". Then we define a Do loop and iterate over all the variables , e.g. 1 ... 4.

In the do loop, we check the conditions and also if there are no variables with the same value. If both are True, we store the variable values in a list using "Sow":

iter = Thread[{vars, 1, 4}] ;
sol = Reap[
  Do[
   If[(And @@ constraints) || ( ! DuplicateFreeQ[vars]), Continue[]];
   Sow[vars];
   , Evaluate[Sequence @@ iter]]
  ][[2]]

{}

As you see, the solution is an empty list. Therefore, there is no magic square with these conditions.

Improve this answer
$\endgroup$
6
  • $\begingroup$ I´m sorry, Daniel...but..I can't understand how to use the code you provided... $\endgroup$ Commented Feb 16 at 16:05
  • $\begingroup$ I mean.. I would first have to run magicSquareILP, so I would still have the same performance issue... $\endgroup$ Commented Feb 16 at 16:15
  • $\begingroup$ Remove the local variables from: magicSquareILP[k_, target_] := Module[{} ...... And remove also "FindInstance[...]". Then run magicSquareILP[4,5t]. Now "vars" and "constraints" are defined. $\endgroup$ Commented Feb 16 at 21:04
  • $\begingroup$ Thanks, Daniel, I'll try it. $\endgroup$ Commented Feb 17 at 14:29
  • $\begingroup$ Yes, it really works very well... but when I change the function to the sum x+y+z with 9 numbers and set the target to 15, which runs like a charm in the original code, this code exceeds the computation time... I think I'm too stupid for this kind of language, I'm going back to Python, hahaha :-D" $\endgroup$ Commented Feb 17 at 14:47

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.