#include "vibes.h"
#include "interval.h"
#include "box.h"
#include <list>

//   more more explanations, see : https://www.ensta-bretagne.fr/jaulin/easibex.html  Example 2


using namespace std;
int m=10;  // number of marks
int q=10;  // number of inliers
vector<double> Ax(m),Ay(m),D(m);  //marks and distances

void Sep(box& X, bool outer)
{   vector<box> Y(m);
    for (int i=0; i<m;i++)
    {  Y[i]=X;
       SinRing(Y[i][1],Y[i][2],Ax[i],Ay[i],D[i]+interval(-0.5,0.5),outer);
    }
    S_q_in(X,q,Y,outer);
}

void sivia()
{   box X(interval(-10,10),interval(-10,10));
    list<box> L;
    int k=0;
    L.push_back (X);
    while ( !L.empty() )
    {   k++; X=L.front();   L.pop_front();
        vibes::drawBox(X[1].inf,X[1].sup,X[2].inf,X[2].sup,"blue");
        Sep(X,1);
        vibes::drawBox(X[1].inf,X[1].sup,X[2].inf,X[2].sup,"magenta");
        Sep(X,0);
        if (k<200)
        {   box X1(2);  box X2(2);
            Bisect (X,X1,X2);
            L.push_back(X1);  L.push_back(X2);
        }
        else vibes::drawBox(X[1].inf,X[1].sup,X[2].inf,X[2].sup,"yellow[yellow]");
    }
}

int main()
{   vibes::beginDrawing();
    vibes::newFigure("SIVIA");
    vibes::clearFigure();
    vector<double> Cx,Cy,Rx,Ry,W;
    for (int k=0;k<m;k++)
    {   Cx.push_back(rand()%20-10);        Cy.push_back(rand()%20-10);
        Rx.push_back(rand()%10+1);         Ry.push_back(rand()%10+1);
        W.push_back(rand()%4);
    }
    for (double t=0;t<1;t=t+0.01)//jusqu a 1
    {   vibes::clearFigure();
        double Px=3*sin(t+1);    double Py=1+4*cos(2*t);
        for (int k=0; k<m;k++)
        {  Ax[k]=Cx[k]+Rx[k]*cos(W[k]*t);
           Ay[k]=Cy[k]+Ry[k]*sin(W[k]*t);
           D[k]=round(sqrt(pow(Px-Ax[k],2)+pow(Py-Ay[k],2)));
           vibes::drawCircle(Ax[k],Ay[k],0.1,"black[black]");
        }
        for (int k1=0;k1<m-q;k1++)
           D[rand()%m]=10; //outlier;
        vibes::drawCircle(Px,Py,0.2,"red[green]");
        sivia();
    }
    vibes::endDrawing();
}
