[BUG] Arithmetic bug with Ibex and power operator

[schmittjoaopedro]

I'm having a bug with the following program:

Model model = new Model();
RealVar x = model.realVar("x", -10000, 10000, 0.001);
RealVar y = model.realVar("y", -10000, 10000, 0.001);

y.eq(x.pow(2.0)).ibex(0.001).post();
model.getSolver().propagate();
System.out.println(x); // Domain should be -> [-100 .. 100], but returns -> [0 .. 100]

Expected behavior

The domain of X is not being filtered correctly. It should return the negative part (-100.0 .. 0) as well.

Context

I'm using this for product configuration, it's very important because influences in the dimensional calculations of engineering products. I found this error during a calculation of physical properties between two different mechanical componentes.

Environment

• Choco-solver version: 4.10.3
• Java version: 1.8 x32 bits

[cprudhom]

Hi,

I suppose, this is related to Ibex:

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();

Let's check with them.

[cprudhom]

Follow Ibex issue

[schmittjoaopedro]

As commented by @gchabert in the issue IBEX-LIB-471, the problem is related about how choco-solver sends the power operations to ibex, as ibex can't handle complex numbers (it only work with positive numbers in the domain of power functions), then the expressions should be changed from:

model.realIbexGenericConstraint("{0}={1}^{2}", y, x, k).post();

to:

model.realIbexGenericConstraint("{0}={1}^" + k, y, x).post();

By this way, the power functions with integer exponents could consider the domain's negative side.

In my tests it worked really well for huge contractions of the domain (let's say greater than 1%). I'm quoting this because there's still one requirement for a fix in ibex that relies on the contraction ratio, already seen in these issues IBEX-LIB-435 and #653 . There's a solution but it's still not available to use, it would be really nice if we could create a patch to solve this issue.

For while, to overcome this problem I've implemented my own power propagator. @cprudhom could you take a look if it seems ok?

package net.weg.maestro.configurator.solver.propagators;

import org.chocosolver.solver.Model;
import org.chocosolver.solver.constraints.Constraint;
import org.chocosolver.solver.constraints.Propagator;
import org.chocosolver.solver.constraints.PropagatorPriority;
import org.chocosolver.solver.exception.ContradictionException;
import org.chocosolver.solver.expression.continuous.arithmetic.CArExpression;
import org.chocosolver.solver.expression.continuous.arithmetic.RealIntervalConstant;
import org.chocosolver.solver.variables.RealVar;
import org.chocosolver.solver.variables.events.RealEventType;
import org.chocosolver.util.ESat;
import org.chocosolver.util.objects.RealInterval;
import org.chocosolver.util.tools.RealUtils;

/**
 * y = x^k
 * <p>
 * Just in case, this is an external library for interval arithmetic.
 * Reference: https://jitpack.io/p/jinterval/jinterval#readme
 */
public class RealPowerProp extends Propagator<RealVar> {

    private RealVar yVar;

    private RealVar xVar;

    private int k;

    /**
     * The integration between choco-solver and ibex have a problem with the
     * contraction ratio. The ibex contractor only updates the domains if the
     * BOUND is increased (or reduced) more than 1%. For example:
     * x = [0, 100]
     * IF x > 1.0 -> ibex updates x to [1, 100]
     * However, IF x > 0.5 -> ibex doesn't updates x to [0, 100]
     * <p>
     * When ibex-lib and choco-solver release the following issues, the ratio parameter
     * will be available and we will be able to start using ibex again.
     * See the following:
     * https://github.com/chocoteam/choco-solver/issues/702
     * https://github.com/chocoteam/choco-solver/issues/653
     * https://github.com/ibex-team/ibex-lib/issues/435
     * https://github.com/ibex-team/ibex-lib/issues/471
     */
    public static boolean isIbex = false;

    public RealPowerProp(RealVar yVar, RealVar xVar, int k) {
        super(new RealVar[]{yVar, xVar}, PropagatorPriority.BINARY, false);
        this.yVar = vars[0];
        this.xVar = vars[1];
        this.k = k;
    }

    @Override
    public int getPropagationConditions(int vIdx) {
        return RealEventType.BOUND.getMask();
    }

    @Override
    public ESat isEntailed() {
        ESat result = ESat.UNDEFINED;
        if (isCompletelyInstantiated()) {
            double lb = IntervalUtils.offsetLB(yVar) - Math.pow(IntervalUtils.offsetLB(xVar), k);
            double ub = IntervalUtils.offsetUB(yVar) - Math.pow(IntervalUtils.offsetUB(xVar), k);
            result = ESat.eval(lb <= 0 && 0 <= ub);
        }
        return result;
    }

    @Override
    public final void propagate(int evtmask) throws ContradictionException {
        boolean hasChanged = true;
        // Projection + propagation
        while (hasChanged) {
            hasChanged = kPower();
            hasChanged |= kRoot();
        }
    }

    private boolean kPower() throws ContradictionException {
        RealInterval x = new RealIntervalConstant(xVar.getLB(), xVar.getUB());
        RealInterval y = new RealIntervalConstant(yVar.getLB(), yVar.getUB());
        RealInterval power = RealUtils.iPower(x, k); // y = x^k
        // For even power, y will range [0, +Inf]
        // For odd power, y will range [-Inf, +Inf]
        // So it doesn't require special handling to update the bounds
        y = IntervalUtils.updateBounds(y, power);
        return yVar.updateBounds(y.getLB(), y.getUB(), this);
    }

    private boolean kRoot() throws ContradictionException {
        boolean isIntervalPositive = IntervalUtils.offsetLB(xVar) >= 0;
        boolean isIntervalNegative = IntervalUtils.offsetUB(xVar) <= 0;
        RealInterval x = new RealIntervalConstant(xVar.getLB(), xVar.getUB());
        RealInterval y = new RealIntervalConstant(yVar.getLB(), yVar.getUB());
        RealInterval root = RealUtils.iRoot(y, k); // x = root(y, k)
        boolean isEven = k % 2 == 0;
        // In all cases, bounds ar calculate as intersection between x and temp
        if (isEven) {
            if (isIntervalPositive) {
                x = IntervalUtils.intersectLowerBound(x, root.getLB());
                x = IntervalUtils.intersectUpperBound(x, root.getUB());
            } else if (isIntervalNegative) {
                x = IntervalUtils.intersectLowerBound(x, -root.getUB());
                x = IntervalUtils.intersectUpperBound(x, -root.getLB());
            } else { // It's crossing zero
                x = IntervalUtils.intersectLowerBound(x, -root.getUB());
                x = IntervalUtils.intersectUpperBound(x, root.getUB());
            }
        } else {
            // For odd roots isn't necessary specific handling as roots can be negative naturally
            x = IntervalUtils.updateBounds(x, root);
        }
        return xVar.updateLowerBound(x.getLB(), this) |
                xVar.updateUpperBound(x.getUB(), this);
    }

    public static CArExpression create(CArExpression exp, int k, double precision) {
        RealVar yVar;
        RealVar xVar = getXVar(exp, precision);
        Model model = xVar.getModel();
        if (isIbex) {
            String expression;
            if (k >= 0) {
                expression = "{0}={1}^" + k;
            } else {
                // y = 1/x^k
                k = -k;
                expression = "{0}=1/{1}^" + k;
            }
            RealInterval powerResult = RealUtils.iPower(xVar, k);
            yVar = model.realVar(powerResult.getLB(), powerResult.getUB(), precision);
            model.realIbexGenericConstraint(expression, yVar, xVar).post();
        } else {
            if (k >= 0) {
                // y = x^k
                RealInterval powerResult = RealUtils.iPower(xVar, k);
                yVar = model.realVar(powerResult.getLB(), powerResult.getUB(), precision);
                new Constraint(model.generateName(), new RealPowerProp(yVar, xVar, k)).post();
            } else {
                // y = 1/x^k
                k = -k;
                RealInterval powerResult = RealUtils.iPower(xVar, k);
                RealVar xPowerK = model.realVar(powerResult.getLB(), powerResult.getUB(), precision);
                new Constraint(model.generateName(), new RealPowerProp(xPowerK, xVar, k)).post();
                RealVar one = model.realVar(1.0, 1.0, precision);
                yVar = model.realVar(powerResult.getLB(), powerResult.getUB(), precision);
                yVar.eq(one.div(xPowerK)).post();
            }
        }
        return yVar;
    }

    private static RealVar getXVar(CArExpression exp, double precision) {
        RealVar xVar;
        if (exp instanceof RealVar) {
            xVar = (RealVar) exp;
        } else {
            xVar = exp.realVar(precision);
        }
        return xVar;
    }

}

[schmittjoaopedro]

@cprudhom I've created the pull request #710 to handle this problem, could you take a look?