Algorithm to generate stars-and-bars alike solutions from a known one?

Question

In my pursue of doing a Simulated Annealing heuristic to solve a problem, I am trying to find the best way to generate neighbours of my current proposed solution.

Solutions come in form of a vector of integer positions (p1,...,pN), which I understand as a binary chain

0 0 0 0 1 0 ... 0 0 1 0 ... 0 1 0 0
       p1          pj        pN

With some restrictions (pj - p(j-1) > D for all j, and p1 > D/2, length - pN > D/2).

Now, my idea is to use something similar to the Levenshtein distance to create new solutions, so if I have [0,1,0,0,1,0] (D=3) and I want a new state within a distance lesser or equal than 1, then I can get [1,0,0,0,1,0], for example, but not [1,0,0,1,0,0].

What I do (in R) is the following:

GenNewSeq <- function(seq, D, dist){
for(i in 1:dist){
    Diffs <- c((ceiling(D/2)+seq[1]),diff(seq),(ceiling(D/2)+seq[length(seq)]))
    position <- sample((2:length(seq))[Diffs > D], size=1)
    move <- sample(c(-1*as.integer(Diffs[position-1]>D),0,1*as.integer(Diffs[position]>D)), size = 1)
    seq[position-1] <- seq[position-1]+move
    }
seq
}

Maybe it is a bit obscure, if you want I can explain better what it does. The thing is that this is 1) slow (I don't know how can I avoid the for), 2) weirdly not working as intended. It tends too much to move only the last positions and/or stabilizing moving forward and backward the same element all the time, so I get biased results on my Simulated Annealing.

I have thought of removing the restriction of distances and put it in the fitness function (something like exp(D-(pj-p(j-1)))), so I can simply move them with normals, or make them move altogether and then oscillate... and I am starting to think that it would be the easiest way. However, I would appreciate very much a reference to how can I do an efficient and reliable algorithm that does what I ask for, I don't mind if I have to do it in C. I have checked this but I wasn't able to solve my doubts.

Thank you very much for your help.

Answer 1

The bug in your program is this. When you select position at random, you are choosing a segment at random from the set of segments of length at least D. The element you are going to end up moving is the right-hand endpoint of this segment.

And, although it seems as though you are choosing the direction of the move at random, in fact the move is more likely to be in the downward direction than upward. This is because Diffs[position-1] is guaranteed to be greater than D (due to the way position was selected), but Diffs[position] is not. Which means that in some cases move is going to be chosen at random from c(-1,0,1) and in other cases it is going to be chosen at random from c(-1,0,0). So, over time, downwards moves will occur more than upwards moves.

Your algorithm can be fixed by selecting at random among all points for which either adjacent segment has length at least D, that where there won't be any bias in move direction:

GenNewSeq2 <- function(seq, D, dist){
  for(i in 1:dist){
    Diffs <- c((ceiling(D/2)+seq[1]),diff(seq))
    bigGaps <- Diffs>D
    moveable <- bigGaps[-1] | head(bigGaps,-1)
    position <- sample(which(moveable),1)
    move <- sample(c(-1*(Diffs[position]>D),1*(Diffs[position+1]>D)), size = 1)
    seq[position] <- seq[position]+move
  }
  seq
}

It is also possible to generate a random new sequence without a for loop. Here is one possible implementation:

newseq<-function(seq,D,dist){
   diffs <- c((ceiling(D/2)+seq[1]),diff(seq))
   bigGaps<-diffs>D
   selected<-sample(which(bigGaps),min(length(bigGaps),dist))
   directions<-sample(c(-1,1),length(selected),T)
   down<-directions<0
   up<-directions>0
   selected[up]<-selected[up]-1
   move<-rep(0,length(seq))
   move[selected[up]]<-1
   move[selected[down]]<-move[selected[down]]-1

   move[length(seq)]<-0  ## the last element of seq stays fixed always
   seq+move
}

This implementation is more efficient, and it doesn't slow down nearly as much when dist grows.

> set.seed(123)
> seq<-sort(sample(1000,20))
> microbenchmark(newseq(seq,20,3),GenNewSeq2(seq,20,3))
Unit: microseconds
                   expr     min       lq  median      uq     max neval
     newseq(seq, 20, 3)  53.503  55.0965  56.026  56.761  68.804   100
 GenNewSeq2(seq, 20, 3) 183.091 188.0490 189.492 191.249 367.094   100
> microbenchmark(newseq(seq,20,6),GenNewSeq2(seq,20,6))
Unit: microseconds
                   expr     min       lq   median       uq     max neval
     newseq(seq, 20, 6)  54.027  56.4960  57.3865  58.2955  70.258   100
 GenNewSeq2(seq, 20, 6) 368.306 373.7745 377.5225 381.4565 559.037   100
> 

We can also verify that GenNewSeq2 and newseq don't drift towards zero by running the following code for each of the three functions, and then plotting the mean of seq over time:

set.seed(12345)
seq<-sort(sample(1000,20))
x<-rep(0,20000)
for(i in 1:20000){
  x[i]<-mean(seq)
  seq<-GenNewSeq(seq,20,3)
}
plot(x,type='l')