#include "sivia.h"


void Sivia::contract_and_draw(Ctc& c, IntervalVector& box, const QColor & pencolor, const QColor & brushcolor) {
    IntervalVector initbox=box;       // get a copy
    try {
        c.contract(box);
        if (box==initbox) return;     // nothing contracted.
        IntervalVector* rest;
        int n=initbox.diff(box,rest); // calculate the set difference
        for (int i=0; i<n; i++) {     // display the boxes
            frame.DrawBox(rest[i][0],rest[i][1],QPen(pencolor),QBrush(brushcolor));
        }
        delete[] rest;
    } catch(EmptyBoxException&) {
        frame.DrawBox(initbox[0],initbox[1],QPen(pencolor),QBrush(brushcolor));
    }
}


Sivia::Sivia(Frame& frame, double epsilon) : frame(frame) {

    // Create the function we want to apply SIVIA on.
    Variable x,y;
    Function f(x,y,sin(x+y)-0.1*x*y);

    NumConstraint c1(x,y,f(x,y)<=2);
    NumConstraint c2(x,y,f(x,y)>=0);
    NumConstraint c3(x,y,f(x,y)>2);
    NumConstraint c4(x,y,f(x,y)<0);

    // Create contractors with respect to each
    // of the previous constraints.
    CtcFwdBwd out1(c1);
    CtcFwdBwd out2(c2);
    CtcFwdBwd in1(c3);
    CtcFwdBwd in2(c4);

    // Create a contractor that removes all the points
    // that do not satisfy either f(x,y)<=2 or f(x,y)>=0.
    // These points are "outside" of the solution set.
    CtcCompo outside(out1,out2);

    // Create a contractor that removes all the points
    // that do not satisfy both f(x,y)>2 or f(x,y)<0.
    // These points are "inside" the solution set.
    CtcUnion inside(in1,in2);

    // Build the initial box.
    IntervalVector box(2);
    box[0]=Interval(-10,10);
    box[1]=Interval(-10,10);

    // Build the way boxes will be bisected.
    // "LargestFirst" means that the dimension bisected
    // is always the largest one.
    LargestFirst lf;

    stack<IntervalVector> s;
    s.push(box);

    while (!s.empty()) {
        // Get a copy of the current box (on top of the stack)
        IntervalVector box=s.top();

        // Remove the box from the stack
        s.pop();
        try {
            // Remove the part that is outside
            contract_and_draw(inside,box,Qt::magenta,Qt::red);
            if (box.is_empty()) { continue; }

            // Remove the part that is inside
            contract_and_draw(outside,box,Qt::darkBlue,Qt::cyan);
            if (box.is_empty()) { continue; }

            // Check if the box is small enough
            if (box.max_diam()<epsilon) {
                frame.DrawBox(box[0],box[1],QPen(Qt::yellow),QBrush(Qt::NoBrush));
            } else {
                // otherwise, bisect it and
                // push the two subboxes on the stack.
                pair<IntervalVector,IntervalVector> boxes=lf.bisect(box);
                s.push(boxes.first);
                s.push(boxes.second);
            }
        } catch(EmptyBoxException&) { }
    }


    frame.Save("paving");
}
