hc4 contractor

[AdBedouhene]

Hi, I've got an issue while trying to make an approximation of a circle using bisections and the hc4 contractor. At some points, it seems that the circle doesn't belong to the set of boxes I generated due to an inaccurate contraction. My constraint is : x² + y² = 1 (I added the constraint x>=0 and y >=0 to generate less boxes and to isolate the problem to the up right part of the circle) The program behaves as follow:

1. After contraction :([0, 1] ; [0, 1])
2. ([0, 1] ; [0, 1]) is not accepted -> bisection:
3. after bisection : ([0, 0.5] ; [0, 1])([0.5, 1] ; [0, 1])
4. Before contraction :([0.5, 1] ; [0, 1])
5. After contraction :([0.5, 1] ; [0, 0.866025])
6. ([0.5, 1] ; [0, 0.866025]) is not accepted -> bisection:
7. after bisection : ([0.5, 1] ; [0, 0.433013])([0.5, 1] ; [0.433013, 0.866025])
8. Before contraction :([0.5, 1] ; [0.433013, 0.866025])
9. After contraction :([0.5, 0.901388] ; [0.433013, 0.866025])
10. ([0.5, 0.901388] ; [0.433013, 0.866025]) is accepted and added to the list
11. Before contraction :([0.5, 1] ; [0, 0.433013])
12. After contraction :([0.901388, 1] ; [0, 0.433013])
13. ([0.901388, 1] ; [0, 0.433013]) is accepted and added to the list
14. Before contraction :([0, 0.5] ; [0, 1])
15. After contraction :([0, 0.5] ; [0.866025, 1])
16. ([0, 0.5] ; [0.866025, 1]) is not accepted -> bisection:
17. after bisection : ([0, 0.25] ; [0.866025, 1])([0.25, 0.5] ; [0.866025, 1])
18. Before contraction :([0.25, 0.5] ; [0.866025, 1])
19. After contraction :([0.25, 0.5] ; [0.866025, 0.968246])
20. ([0.25, 0.5] ; [0.866025, 0.968246]) is accepted and added to the list
21. Before contraction :([0, 0.25] ; [0.866025, 1])
22. After contraction :([0, 0.25] ; [0.968246, 1])
23. ([0, 0.25] ; [0.968246, 1]) is accepted and added to the list

The resulting boxes are:

• r([0.50000000000000011102, 0.90138781886599739224] ; [0.43301270189221935381, 0.86602540378443870761])
• g([0.90138781886599717019, 1.000000000000000222] ; [0, 0.43301270189221935381])
• b([0.25000000000000005551, 0.50000000000000011102] ; [0.86602540378443837454, 0.96824583655185436637])
• y([0, 0.25000000000000005551] ; [0.9682458365518540333, 1.000000000000000222])

The boxes b and y perfectly meet at the point (0.25000000000000005551, 0.9682458365518540333) but the distance of this point to the centre (0.0) is [0.99999999999999955591, 0.99999999999999977796] < 1. It means that the circle goes outside the y box then goes inside the b box.

As we can see it in the line 18 and 19 of the result, the box b results from a non accurate contraction from [0.25, 0.5] ; [0.866025, 1] to [0.866025, 0.968246].

So I don't know if I'm doing the bisections and contractions wrong, or if it is due to a bug.

There is the code:

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

using namespace std;
using namespace ibex;
using namespace vibes;

void hc4contract(IntervalVector & bounds){
    Variable vx, vy;
    SystemFactory fac;
    fac.add_var(vx);
    fac.add_var(vy);
    Interval cercle = Interval(1.).inflate(0.);
    fac.add_ctr(pow(vx,2) + pow(vy,2) = cercle);
    fac.add_ctr(vx>= 0.);
    fac.add_ctr(vy>=0.);

    System sys(fac);
    ibex::CtcHC4 hc4(sys);
    cout << " Before contraction  :" << bounds << endl;

    hc4.contract(bounds);
    cout << " After contraction  :" << bounds << endl;
}

void border(list<IntervalVector> &sol, IntervalVector box, double epsilon){

    list <IntervalVector> s;
    s.emplace_back(box);

    while(!s.empty()){
        IntervalVector bis = s.front();
        s.pop_front();
        hc4contract(bis);

        bool emptiness;
        emptiness=bis.is_empty();
        if(emptiness) {
            cout << "contraction is empty !"<<endl;
            continue;
        }

        if(bis.max_diam()<=epsilon){
            cout << bis <<" is accepted and added to the list " << endl;
            sol.emplace_back(bis);
        }
        else{
            cout << bis <<" is not accepted -> bisection: "<<endl;
            int i = bis.extr_diam_index(false);
            pair<IntervalVector, IntervalVector> b=bis.bisect(i);
            cout << "after bisection : "<<b.first << b.second << endl;
            s.emplace_front(b.first);
            s.emplace_front(b.second);

        }

    }
}

int main(){
    vibes::beginDrawing();
    newFigure ("fig_boxes");
    const int n=2;
    Interval domain(0,1);

    IntervalVector boxes(n);

    list <IntervalVector> solution;
    double epsilon(0.5);
    border(solution,boxes,epsilon);

    list <IntervalVector>::iterator i;

//Color parameter to show the resulting boxes in diffrent colors
//this parameter should be removed for a epsilon != 0.5
    char col[4];
    col[0] = 'r';
    col[1] = 'g';
    col[2] = 'b';
    col[3] = 'y';
    int ii(0);
    for(i=solution.begin(); i!=solution.end(); i++){

        IntervalVector a=*i;
        string sii;
        sii= col[ii];
        cout << col[ii] << " : ";

        cout<< a << endl;

        drawBox(a[0].lb(), a[0].ub(),a[1].lb(),a[1].ub(),sii);
        ii++;
    }
    drawCircle(0.,0.,1);
    vibes::endDrawing();

    return 0;
}

[amarendet]

Hi,

Just taking your example for boxes b and y, they do not meet perfectly at the same point, they slightly overlap : b & y = (<0.25, 0.2500000000000001> ; [0.968245836551854, 0.9682458365518544]). And in fact, with f(x, y) = x^2+y^2-1, f(b & y) = [-4.440892098500627e-16, 4.440892098500627e-16].

You can verify that with this code:

#include "ibex.h"
using namespace ibex;

int main(int argc, char** argv) {
    IntervalVector a{{0.25000000000000005551, 0.50000000000000011102},{0.86602540378443837454, 0.9682458365518543663}};
    IntervalVector b{{0, 0.25000000000000005551} ,{ 0.9682458365518540333, 1.000000000000000222}};
    Function f("x", "y", "x^2+y^2-1");
    std::cout << (a&b) << std::endl;
    Interval res = f.eval(a & b);
    std::cout << res << std::endl;
    return 0;
}

Regards

[AdBedouhene]

Hi, Thanks for answering so quickly. You're right, the circle goes across the overlapping part. The vibes drawing shows the circle going outside the boxes, I shouldn't have trusted that. Regards