Arithmetic bug with Ibex and power operator (wrt Choco)

[cprudhom]

Dear developers,

It seems that something is going wrong with Ibex+Choco, when dealing with pow operators.

Running this with Ibex-2.8.7 (without choco):

Ibex ibex = new Ibex(new double[]{0.001, 0.001, 0.001});
        
ibex.add_ctr("{0}={1}^{2}");
ibex.build();
double domains[] = {-16.0, 16.0, -4.0, 4.0, 2.0, 2.0};
System.out.println("Before contract:");
System.out.println("([" + domains[0] + "," + domains[1] + "] ; [" + domains[2] + "," + domains[3] + "])");
        
int result = ibex.contract(0, domains);
        
if (result == Ibex.FAIL) {
    System.out.println("Failed!");
} else if (result == Ibex.CONTRACT) {
    System.out.println("After contract:");
    System.out.println("([" + domains[0] + "," + domains[1] + "] ; [" + domains[2] + "," + domains[3] + "])");
} else {
    System.out.println("Nothing.");
}
ibex.release();

prints this:

Before contract:
([-16.0,16.0] ; [-4.0,4.0])
After contract:
([0.0,16.0] ; [0.0,4.0])

(cf. #702)

Any ideas?

Best,
CP

[schmittjoaopedro]

I've done some experiments in C++ using ibex 2.8.7. Could I suppose that it's a bug on the integration of choco and ibex?
Because it seems that ibex is working correctly (Ref: contractors), and choco is calling ibex correctly as well.

Following, the program I tested (foo.cpp):

#include "ibex.h"
#include <iostream>

using namespace std;
using namespace ibex;

int main(int argc, char** argv) {
	// create a constraint on (x,y)
	Variable x, y;
	NumConstraint c(x,y,y=sqr(x));

	// create domains for x and y
	IntervalVector box_x(1, Interval(-10000,10000));
	IntervalVector box_y(1, Interval(-10000,10000));

	// set the precision that controls how much y will be bisected
	double epsilon=0.01;

	// create a contractor on x by transforming y into an
	// existentially-quantified parameter.
	CtcExist exist_y(c,y,box_y,epsilon);

	cout << "Before: X interval:" << box_x << endl;
	cout << "Before: Y interval:" << box_y << endl;
	
	exist_y.contract(box_x);
	
	cout << "After: X interval:" << box_x << endl;
	cout << "After: Y interval:" << box_y << endl;
}

The program result is:

jschmitt@BRJGSD333007 /c/MinGW/msys/1.0/home/jschmitt/Ibex/ibex-2.8.7/examples
$ make foo
g++  -frounding-math -ffloat-store -msse3 -Ic:/MinGW/msys/1.0/home/jschmitt/Ibex/ibex-2.8.7/include -Ic:/MinGW/msys/1.0/home/jschmitt/Ibex/ibex-2.8.7/include/ibex -Ic:/MinGW/msys/1.0/home/jschmitt/Ibex/ibex-2.8.7/include\\ibex\\3rd    -O3 -DNDEBUG  -std=c++11 -DIBEX_BENCHS_DIR=\"../benchs/solver\" -U__STRICT_ANSI__ -o foo foo.cpp -Lc:/MinGW/msys/1.0/home/jschmitt/Ibex/ibex-2.8.7/lib -libex

jschmitt@BRJGSD333007 /c/MinGW/msys/1.0/home/jschmitt/Ibex/ibex-2.8.7/examples
$ ./foo
Before: X interval:([-10000, 10000])
Before: Y interval:([-10000, 10000])
After: X interval:([-100, 100])
After: Y interval:([-10000, 10000])

I'm really depending of this fix, could you sugest a workaround until the released of the fixed version?

[gchabert]

Thanks Charles for reporting and João Pedro for investigating
I don't have a simple fix for that.

The problem is that negative values fox x are not allowed in the expression x^y because Ibex is working with continuous domains. Just take as example the value x=-1. As long as y is a positive integer, you can perfectly write (-1)^y, the result is a real number, either 1 or -1. But if y is real, (-1)^y can be a complex number. E.g., (-1)^0.5 = sqrt(-1) = i or -i !
Of course, Ibex does not work with complex numbers so negative values are considered as outside the definition domain of the power function, unless the power is a constant positive integer. This is why your example works João Pedro.

I cannot change that because x^y is rewritten as exp(y*log(x)) to be processed, a correct transformation if x>0, and actually, a definition for power with real exponents.

You may suggest to check dynamically the domain of y and apply a different procedure in the case of an integral singleton; but that would immediately leads to inconsistent non-monotonic behavior (y = [2,2.000001] would be more contracting than y =[2,2]...)

I am afraid there is no workaround here because we cannot change maths!
The only solution I see is to only allow constant exponents in choco/ibex.

ibex.add_ctr("{0}={1}^2");

Would that be so restricting?

[schmittjoaopedro]

I see, for me that's fine. I can handle the constant exponent by at the Java side.

So, I did that to execute some tests, and I got stuck in another problem. If I execute the following code:

Ibex ibex = new Ibex(new double[]{0.001, 0.001});

ibex.add_ctr("{0}={1}^3");
ibex.build();
double domains[] = {9.5, 1.0E12, 0, 1.0E5}; // Contracted: ([9.5,1.0E12] , [2.117911792127444,10000.000000000126])
//double domains[] = {9.5, 1.0E12, 0, 1.0E4}; // Didn't contract: ([9.5,1.0E12] , [0.0,10000.0])
System.out.println("Before contract:");
System.out.println("([" + domains[0] + "," + domains[1] + "] , [" + domains[2] + "," + domains[3] + "])");

int result = ibex.contract(0, domains);

if (result == Ibex.FAIL) {
	System.out.println("Failed!");
} else if (result == Ibex.CONTRACT) {
	System.out.println("After contract:");
	System.out.println("([" + domains[0] + "," + domains[1] + "] , [" + domains[2] + "," + domains[3] + "])");
} else {
	System.out.println("Nothing.");
}
ibex.release();

The domain is not contracted for 1.0E4, I think that it's related to the contraction ratio from IBEX (already seen in this ISSUE-435. Is that correct? I'm kind really need this parameter, is it possible to create a patch for this version be integrated in choco-solver @gchabert ? What do you think about it @cprudhom?

[gchabert]

@schmittjoaopedro It's not yet documented but the Java plugin is now located in a separate github project. You can already download and install the current code under the ibex/plugins folder, to benefit from the contraction ratio parameter.
A new release (2.8.8) is going to be issued in a couple of days with the new installation instructions.

[gchabert]

By the way, I'm closing and moving this issue to the other repository.
In the future, it's preferable to post issues there. We can separate this way problems related to the core lib and those related to the java interface.