PropagationGuidedNeighborhood computes the wrong value of logSum

[GustavBjordal]

It looks like someone forgot to delete lines 106 to 108 in PropagationGuidedNeighborhood.java.

In previous versions, line 105 was logSum = 0 but was at some point this updated to be logSum = Arrays.stream(variables).mapToDouble(v -> Math.log(v.getDomainSize())).sum();, which directly computes the logSum and makes lines 106 to 108 obsolete.
However, because lines 106 to 108 remain and increase the value of logSum, the resulting value is actually two times the intended value.

[cprudhom]

Thank you for the bug report.
That will be integrated in the next release.

Kind regards,
CP

[GustavBjordal]

Great!

One more thing, should line 150 be .sorted(Comparator.comparingInt(i -> -all[i])) as sorted sorts in increasing order and we want the first element to be the ones with the most propagation rather than the least, which is what the current implementation would do?

[cprudhom]

The article states that "the list keeps the most affected variables".
Since we put in all[i] the difference between dsize[i] (domain size of the _i_th variable at root node) and ds[i] (current domain size of the _i_th), the smaller bigger the value of all[i], the more the variable i is affected (ie, reduction is bigger). So, you are right the current code is not correct.
(edited)

[cprudhom]

It seems that ReversePropagationGuidedNeighborhood is not correct too.

[mr-tl]

Great!

One more thing, should line 150 be .sorted(Comparator.comparingInt(i -> -all[i])) as sorted sorts in increasing order and we want the first element to be the ones with the most propagation rather than the least, which is what the current implementation would do?

A more clear way to change this sorting may be to use reversed, as follows:

                // 4. update variable list
                candidates = IntStream.range(0, n)
                    .filter(i -> fragment.get(i) && all[i] > 0)
                    .boxed()
                    .sorted(Comparator.comparingInt((Integer i) -> all[i]).reversed())
                    .limit(listSize)
                    .collect(Collectors.toList());

/Magnus

[cprudhom]

@mr-tl
I suppose you mean that:

candidates = IntStream.range(0, n)
        .filter(i -> fragment.get(i) && all[i] > 0)
        .boxed()
        .sorted(Comparator.comparingInt(i -> all[(int) i]).reversed())
        .limit(listSize)
        .collect(Collectors.toList());

What about :

candidates = IntStream.range(0, n)
        .filter(i -> fragment.get(i) && all[i] > 0)
        .boxed()
        .sorted(Comparator.comparingInt(i -> all[i]))
        .sorted(Comparator.reverseOrder())
        .limit(listSize)
        .collect(Collectors.toList());

?

[mr-tl]

@mr-tl
I suppose you mean that:

candidates = IntStream.range(0, n)
        .filter(i -> fragment.get(i) && all[i] > 0)
        .boxed()
        .sorted(Comparator.comparingInt(i -> all[(int) i]).reversed())
        .limit(listSize)
        .collect(Collectors.toList());

Yes, the only difference being that I stated the type of i in place of the cast.

What about :

candidates = IntStream.range(0, n)
        .filter(i -> fragment.get(i) && all[i] > 0)
        .boxed()
        .sorted(Comparator.comparingInt(i -> all[i]))
        .sorted(Comparator.reverseOrder())
        .limit(listSize)
        .collect(Collectors.toList());

?

Hmm, but that would double the amount of work, no? Two sorted passes instead of one.

And your second sorted is on the indices, and not on the all entries. So I believe this is not correct.

/Magnus

[cprudhom]

Hmm, but that would double the amount of work, no?

I agree, I'm not a huge fan of this, even if that works.

Yes, the only difference being that I stated the type of i in place of the cast.

Yes, and that the code you gave does not compile :)

I will try JMH to check what is the best option, and let us know.

[cprudhom]

The JMH file I executed is here.

With n = 10000:

MyBenchmark.met1  thrpt    2  1145,850          ops/s
MyBenchmark.met2  thrpt    2  1120,376          ops/s
MyBenchmark.met3  thrpt    2  1116,059          ops/s
MyBenchmark.met4  thrpt    2   619,583          ops/s

With n = 1000 (random change density of fragment):

MyBenchmark.met1  thrpt    2  19271,951          ops/s
MyBenchmark.met2  thrpt    2  17883,254          ops/s
MyBenchmark.met3  thrpt    2  17599,031          ops/s
MyBenchmark.met4  thrpt    2   8888,857          ops/s

Results came in. I suppose we should ignore the last option.
Since performances are also expected here, I vote for the first version, even if it is not the clearest

CP

[mr-tl]

I see that your commit for the reversed order uses both negating the sorting value:

.sorted(Comparator.comparingInt(i -> -all[(int) i]))

and then the Comparator.reverseOrder().

However, as I said in my previous comment, Comparator.reverseOrder() will in your case only sort the entries according to the (reversed) indices. So your result will have nothing to do with the values of all. This cannot be what you want.

See the documentation for Comparator.reverseOrder():

     * Returns a comparator that imposes the reverse of the <em>natural
     * ordering</em>.
     *
     * <p>The returned comparator is serializable and throws {@link
     * NullPointerException} when comparing {@code null}.
     *
     * @param  <T> the {@link Comparable} type of element to be compared
     * @return a comparator that imposes the reverse of the <i>natural
     *         ordering</i> on {@code Comparable} objects.
     * @see Comparable
     * @since 1.8
     */

[cprudhom]

and then the Comparator.reverseOrder().

Fixed, thank you for your vigilance.

See the documentation

I should have checked that in the first place, you're right.

[GustavBjordal]

The implementation looks correct now, thank for taking care of it!

However, I have been going through the original PGLNS paper and the implementation here, and there something that's been bugging me: why do we have the candidates list in the first place?
The candidates list is rebuilt from scratch after each variable is fixed, meaning an array is allocated, populated, sorted, and trimmed. This seems unnecessarily expensive considering how it is used later on in the code, as I will try to explain here.

But first, I would like to point out that the original paper is actually not very clear on how this list is to be used. The paper essentially states that there is a list, it is of size 10, it is sorted by how much propagation that has occurred, and it is used to select which variable to freeze. However, it is not described how it is used to select a variable, except for in the experiments section which states that they "always select the most affected variable", i.e., the first element of the list. But the paper never explains why the list is then of length 10 and not of length 1, nor does it explain when any but the first element of the list is selected.

The current implementation in Choco can essentially be boiled down into the following code w.r.t the candidates list, sorry for the crappy pseudo-code:

function update() {
  while(not done) {
    v = selectVariable();
    if(v can be frozen to the previous value) {
      throw away the old candidates list and compute a new one from scratch;
    } else {
      ignore v;
    }
  }
}

function selectVariable(){
  if(candidate list empty){
    return a random variable to freeze;
  }else{
    pop and return the head of the candidates list;
  }
}

Here I would argue that the else ignore v case happens very rarely, but can happen due to the bounds of the objective variable being different compared to when the previous value of v was feasible. Therefore, I suspect that the code in practice almost never looks at even the second element of the candidate list. Note also that the return a random variable to freeze only happens when the first variable is selected and when the else ignore v has happened 10 times in a row.

It therefore seems much more reasonable for me to simply compute the head of the candidates list on demand in the selectVariable function by just computing the most affected variable there, which would take O(n) time and O(1) space, as oppose to (n logn) time and O(n) space. One can also just add a counter in the else ignore v and change the line if(candidate list empty){ to if(first variable or counter is equal to 10){.
Not only does it seem like this would be more efficient but it would, in my opinion, also be clearer which variable we select, specifically using some sort of maxBy(), as oppose to sorting a list (which caused the bug in the first place) and selecting from that list somewhere else.

Sorry for the long post, normally I would make a pull request but this was faster than setting up a development environment, doing tests, etc. Still, I hope the point of my suggestion comes across.

[GustavBjordal]

One more thing, related to Reverse PGLNS (RPGLNS):

When going through the PGLNS paper yet again I noticed that RPGLNS is ment to be used together with PGLNS and not as a stand-alone thing as in the Choco implementation.
It seems that, in order to run RPGLNS, one first runs PGLNS for an unspecified amount of time to collect the closeness data, and then run RPGLNS, possibly repeating this loop (this is not clear). Note the line in section 4.3: "We decided to define this closeness relation as the average volume of propagation [...] other [sic; should be 'over'] the previous runs of direct [emphasis added] PGLNS".
Actually, in the experiments they apparently also interleave this with random LNS in order to actually get the good performance.

Furthermore, it seems like the epsilon parameter of RPGLNS should also be used in PGLNS and this is in fact claimed to be crucial: "Experimental results show that this epsilon is a crucial improvement to the PGLNS and the reverse PGLNS as only by using this method were we able to beat our hand-written neighbourhoods."

[cprudhom]

Well, there are many things to discuss.
First the easy one:

Epsilon should be used in PGLNS

I confess that I removed that in the last commit, based on diagonal reading of the article.
Somehow, dynamicity was maintained by relaxing fgmtSize as long as no solution is found.
But, re-introducing epsilon is more consistent whit both the article and history.

Then the more tricky ones.

the else ignore v case happens very rarely

It is not clear to me how often this case happens. If you have any clues, I'll take it.

the code in practice almost never looks at even the second element of the candidate list.

Once again, I don't know.
In the article (§ 4.1), it is said that "Then for each variable that does exhibit a change in its domain size, we update the list with these variables". How to update the list is not defined, it supposes that a function that merges propagation effects exists...
In our implementation though, candidates list is rebuild each time a variable is moved to the fragment.
Without merging operator, I suppose that keeping the most affected variable is enough.

Note also that the return a random variable to freeze only happens (...) when the else ignore v has happened 10 times in a row.

In our implementation, I would say no. Indeed, candidates list is build from scratch at each step of the while-loop. In practice, the random call is made once for all at the beginning, and very rarely after.

to simply compute the head of the candidates list on demand in the selectVariable function

I agree, under those conditions, it is a better idea, unless we are able to propose a new operator.

One can also just add a counter

candidates list is build from scratch at each step of the while-loop, so this is not possible.

RPGLNS is ment to be used together with PGLNS (...) possibly repeating this loop (this is not clear).

You are right, we need to add a map that maintains the closeness relation per pair of variables.