Le Monde puzzle [#882]

A terrific Le Monde mathematical puzzle:

All integers between 1 and n² are written in an (n,n)  matrix under the constraint that two consecutive integers are adjacent (i.e. 15 and 13 are two of the four neighbours of 14). What is the maximal value for the sum of the diagonal of this matrix?

Indeed, when considering a simulation resolution (for small values of m), it constitutes an example of self-avoiding random walk: when inserting the integers one by one at random, one produces a random walk over the (n,n) grid.

While the solution is trying to stick as much as possible to the diagonal vicinity for the descending sequence n²,n²-1, &tc., leaving space away from the diagonal for the terminal values, as in this example for n=5,

25 22 21 14 13
24 23 20 15 12
01 02 19 16 11
04 03 18 17 10
05 06 07 08 09

simulating such a random walk is a bit challenging as the brute force solution does not work extremely well:

n=5
n2=n^2

init=function(){
  board=invoard=matrix(0,n,n)
  set=rep(0,n2)

  start=val=n2 #sample(1:n2,1)
  neigh=board[start]=val
  invoard[val]=start
  set[val]=1

  return(list(board=board,invoard=invoard,
  set=set,neigh=neigh))
  }

voisi=function(i){
 a=arrayInd(i,c(n,n))
 b=a[2];a=a[1]

 if ((a>1)&(a<n)){
  voizin=(a+c(-1,1))+(b-1)*n}else{
  if (a==1) voizin=a+1+(b-1)*n
  if (a==n) voizin=a-1+(b-1)*n}

 if ((b>1)&(b<=n)){
  voizin=c(voizin,a+(b+c(-2,0))*n)}else{
  if (b==1) voizin=c(voizin,a+b*n)
  if (b>n) voizin=c(voizin,a+(b-2)*n)}

 voizin=voizin[(voizin>0)&(voizin<=n2)]
 return(voizin)
}

a=init()
board=a$board
invoard=a$invoard
set=a$set
neigh=a$neigh

while (sum(set)<n2){

  if (length(neigh)==1) neighb=neigh+c(-1,1)
  if (length(neigh)==2) neighb=c(min(neigh)-1,
       max(neigh)+1)
  neighb=neighb[(neighb>0)&(neighb<=n2)]
  neighb=neighb[set[neighb]==0]

  for (i in 1:length(neighb)){

    j=1
    if (i==2) j=1+(length(neigh)>1)
    loc=voisi(invoard[neigh[j]])
    loc=loc[board[loc]==0]

    if (length(loc)==0) break() #no solution
    if (length(loc)==1){
         board[loc]=neighb[i]
         invoard[neighb[i]]=loc
         set[neighb[i]]=1}
    if (length(loc)>1){ #2 or more solutions
      val=sample(loc,1)
      board[val]=neighb[i]
      invoard[neighb[i]]=val
      set[neighb[i]]=1}
    }

  if (min(set[neighb])==0){#start afresco
    a=init()
    board=a$board
    invoard=a$invoard
    set=a$set
    neigh=a$neigh
  }else{
   if (length(neighb)==1) neigh=neighb
    if (length(neighb)>1){
      neigh=sort(neighb)
      neigh=neigh[(neigh>0)&(neigh<=n2)]
      }}
}

the reason being that the chain often has to restart a fresco, the more the larger n is… For n=5, I still recovered the optimal solution:

> while (sum(diag(board))<93) 
+ source("lemonde882.R")
     [,1] [,2] [,3] [,4] [,5]
[1,]    9    8    7    6    1
[2,]   10   17   18    5    2
[3,]   11   16   19    4    3
[4,]   12   15   20   23   24
[5,]   13   14   21   22   25

but running the R code for n=7 towards finding the maximum (259?) takes quite a while and 50 proposals for n=8 took the whole night… If I run a simple log-log regression on the values obtained for n=2,…,7, the prediction for n=10 is 768. A non-stochastic prediction is 870.

As a last ditch attempt to recover the sequence n²-2k along the diagonal for k=0,1,…, I modified my code to favour simulations close to the diagonal at the start, as

if (length(loc)>1){#two or more solutions
   val=sample(loc,1,
   proba=exp(-abs(loc%%n-1-loc%/%n)/sum(set)))

which produces higher diagonal values but also more rejections.