Note

This tutorial is for the Ibex core library.
If you just want to solve **equations** or an **optimization problem**,
jump to the documentation of the *IbexOpt plugin*.

To write a program with Ibex, use the following canvas:

```
#include "ibex.h"
using namespace std;
using namespace ibex;
int main(int argc, char** argv) {
// write your own code here
}
```

You can execute by yourself all the code snippets of this tutorial, using this canvas.

To compile a program, the easiest way is to copy-paste the `makefile` of the `examples/` subfolder of Ibex.
See also *Compiling a Test Program*.

Here are examples of intervals

```
Interval x(1,2); // create the interval [1,2]
Interval y; // create the interval (-oo,oo)
Interval z=Interval::ALL_REALS; // create the interval (-oo,oo)
Interval w=Interval::EMPTY_SET; // create the empty interval
```

C++ operator overloading allows you to calculate the sum of two intervals by using directly the “+” symbol:

```
// - create the inteval x=[1,2] and y=[3,4]
// - calculate the interval sum x+y
Interval x(1,2);
Interval y(3,4);
cout << "x+y=" << x+y << endl; // display [4,6]
```

You can use the other operators similarly (`-`, `*`, `/`).

All the elementary functions can be applied to intervals, and composed in an arbitrarly way:

```
Interval x(0,1);
Interval y=exp(x+1); //y is [1,7.389...]
```

You can create an interval vector by using an intermediate array of n*2 `double`, representing the lower and
uppoer bounds of each components.
The first argument of the constructor of `IntervalVector in this case is the dimension (here, 3),
the second the array of ``double`.

```
double _x[3][2]={{0,1},{2,3},{4,5}};
IntervalVector x(3,_x); //create ([0,1],[2,3],[4,5])
```

You can also create an interval vector by duplicating a given interval or simply create the empty interval vector.

```
IntervalVector x(3,Interval(1,2)); //create ([1,2],[1,2],[1,2])
IntervalVector y=IntervalVector::empty(3); //create a vector of 3 empty intervals
```

Interval matrices can be created in a similar way. However, since we cannot build 3-dimensional arrays in C++,
all the bounds must be set in a single n*2 array representing the matrix row by row (and n is the total number of entries
of the matrix). The two first arguments of the constructor are the number of rows and columns respectively. The last one
is the array of `double`.
Here is an example of a 3x3 matrix:

```
double _M[9][2]={{0,1},{0,1},{0,1},
{0,2},{0,2},{0,2},
{0,3},{0,3},{0,3}};
IntervalMatrix M(3,3,_M);
//create (([0,1] [0,1] [0,1]) ; ([0,2] [0,2] [0,2]) ; ([0,3] [0,3] [0,3]))
```

You can use the usual operations of linear algebra between matrices and vectors (*sum of vectors,
transpose of vectors, sum of matrices, left multiplication of a matrix by a scalar, etc.*).

```
// ------------------------------------------------
// Vector/matrix interval arithmetic
// - create an interval vector x
// - create an interval matrix M
// - calculate M*x
// - calculate M'*x, where M' is the transpose of M
// ------------------------------------------------
double _x[3][2]={{0,1},{2,3},{4,5}};
IntervalVector x(3,_x);
double _M[9][2]={{0,1},{0,1},{0,1}, // 3*3 matrix of intervals
{0,2},{0,2},{0,2},
{0,1},{0,1},{0,1}};
IntervalMatrix M(3,3,_M);
IntervalVector y=M*x; // matrix-vector multiplication
IntervalMatrix N=M.transpose(); // N is M^T
```

These usual properties can be obtained for intervals. They are also all extended to interval vectors or matrices componentwise. For instance, the radius of an interval matrix is the (real) matrix of the radii.

As a consequence, Ibex also has classes to handle real (versus interval) vectors and matrices. Mathematical Operations (like the sum) can also be applied to these objects but, of course, using this times floating-point arithmetic (not interval).

```
// ------------------------------------------------
// Mixing real/interval vector/matrices
// - calculate the magnitude of an interval matrix (a real matrix)
// - calculate the midvector of an interval vector (a real vector)
// - multiply the latters (floating point arithmetic)
// ------------------------------------------------
double _x[][2]={{0,1},{0,1},{0,1}};
IntervalVector x(3,_x);
double _M[9][2]={{0,1},{0,1},{0,1},
{0,2},{0,2},{0,2},
{0,1},{0,1},{0,1}};
IntervalMatrix M(3,3,_M);
Matrix M2=M.mag(); // the matrix of magnitudes
Vector x2=x.mid(); // the vector of midpoints
Vector y=M2*x2; // a matrix-vector product (subject to roundoff errors)
```

The easiest way to create a function is with a string directly:

```
Function f("x","y","sin(x+y)"); // create the function (x,y)->sin(x+y)
```

However, this has some limitations (see *Advanced examples*).
Another (more flexible) way to create function is using C++ operator overloading.
The only difference is that you must have to first build *variables*:

```
Variable x("x");
Variable y("y");
Function f(x,y,sin(x+y)); // create the function (x,y)->sin(x+y)
```

You can insert interval constants in the expresion of a function, even in C++-style. For instance, if you want to create the function \(x\mapsto\sin(2x)\), just write:

```
Variable x;
Function f(x,sin(2*x)); // create the function (x,y)->sin(2*x)
```

Assume now that the function to be created is \(x\mapsto\sin(\pi x)\). It is still possible to
use a `double` representing approximately \(\pi\); but to keep
numerical reliability, it is required in this case to use an interval constant enclosing
\(\pi\). Next function must be seen as a “thick” function that rigorously encloses \(\sin(\pi x)\):

```
Interval pi(3.1415,3.1416);
Variable x;
Function f(x,sin(pi*x)); // create the function (x,y)->sin(π*x)
```

Or with strings directly:

```
Function f("x","sin([3.1415,3.1416]*x)"); // create the function (x,y)->sin(π*x)
```

Arguments of a function are not necessarily scalar variables. They can also be vectors or matrices. In the following example, we build the distance function: \(dist:(a,b)\mapsto\| a-b \|\) where a and b are 2-dimensional vectors.

```
Variable a(2);
Variable b(2);
Function dist(a,b,sqrt(sqr(a[0]-b[0])+sqr(a[1]-b[1])));
```

We can also create the same function with string directly (note that the syntax quite *differs*).

```
Function dist("a[2]","b[2]","sqrt((a(1)-b(1)^2)+(a(2)-b(2))^2)");
```

*Note*: *Evaluation* of a thick function will necessarily result in an interval with non-null diameter, even if the argument is reduced to a point.

You can compose functions to build new functions. We build here the function that maps a vector x to its distance with a constant point (1,2). To this end, we first define a generic distance function dist(a,b) as above.

```
/* create the distance function with 2 arguments */
Variable a(2);
Variable b(2);
Function dist(a,b,sqrt(sqr(a[0]-b[0])+sqr(a[1]-b[1])));
/* create the constant vector pt=(1,2) */
Vector pt(2);
pt[0]=1;
pt[1]=2;
/* create the function x->dist(x,pt). */
Variable x(2);
Function f(x,dist(x,pt));
```

The display is as folllows. Note that constant values like 0 are automatically replaced by degenerated intervals (like [0,0]):

```
f:(x)->(dist(x,(<0, 0> ; <0, 0>));dist(x,(<1, 1> ; <1, 1>)))
```

Let us start with a basic example: the function \(x\mapsto(x-1,x+1)\).

With strings:

```
Function f("x","(x-1,x+1)");
```

With operator overloading:

```
Variable x;
Function f(x,Return(x-1,x+1));
```

*Note:* The `Return` keyword is only necessary when the output of a function is a vector (or a matrix).

Now, in line with the previous sections, let us define a more complicated example:
the function that associates to a vector x its distance with two fixed points `pt1` and `pt2`
initialized in our program to (0,0) and (1,1):

\[f:x\mapsto(\| x-(1,1)\|, \| x-(0,0)\|).\]

```
// ------------------------------------------------
// Vector-valued functions
// ------------------------------------------------
/* create the distance function with 2 arguments */
Variable x(2,"x");
Variable pt(2,"p");
Function dist(x,pt,sqrt(sqr(x[0]-pt[0])+sqr(x[1]-pt[1])),"dist");
/* create the two constant vectors */
Vector pt1=Vector::zeros(2);
Vector pt2=Vector::ones(2);
/* create the function x->(dist(x,pt1),dist(x,pt2)). */
Function f(x,Return(dist(x,pt1),dist(x,pt2)),"f");
cout << f << endl;
```

The last construction is much more cumbersome with strings.

You can also create functions that return matrices. Here is an example of a function from \(R\) to \(R^{2\times 2}\) where:

\[f: x \mapsto ( (2x, -x) ; (-x,3x) ).\]

With strings:

```
Function f("x","((2*x,x);(-x,3*x))");
```

With C++ operator overloading:

```
Variable x("x");
Function f(x,Return(Return(2*x,x,true),Return(-x,3*x,true)));
```

The boolean value `true` given here to the two embedded `Return`
means that, each time, the two components must be put in rows, and not in column as it is by default.
In contrast, the enclosing `Return` keeps the default behaviour since the two rows are
put in column in order to form a 2x2 matrix.

To create sophisticated functions we advice you to use an intermediate “minibex” input file as follows instead of embedding the function directly in your C++ program. The previous example can be written in a plain text file:

```
function f(x)
return ((2*x,-x);(-x,3*x));
end
```

Save this file under the name “myfunction.txt”. Now, you can load this function in your C++ program:

```
Function f("myfunction.txt");
```

When several occurrences of the same subexpression occur in a function, it is a good idea for readibility (and, actually, efficiency) to put this subexpression into intermediate variables.

The following example is the function thar returns the rotation matrix from the three Euler angles. In this function an experssion like cos(phi) occurs several times.:

```
/* Computes the rotation matrix from the Euler angles:
roll(phi), the pitch (theta) and the yaw (psi) */
function euler(phi,theta,psi)
cphi = cos(phi);
sphi = sin(phi);
ctheta = cos(theta);
stheta = sin(theta);
cpsi = cos(psi);
spsi = sin(psi);
return
( (ctheta*cpsi, -cphi*spsi+stheta*cpsi*sphi, spsi*sphi+stheta*cpsi*cphi) ;
(ctheta*spsi, cpsi*cphi+stheta*spsi*sphi, -cpsi*sphi+stheta*cphi*spsi) ;
(-stheta, ctheta*sphi, ctheta*cphi) );
end
```

Given input `double` values x, you can obtain a rigorous inclusion of f(x) either using
`eval`, `eval_vector` or `eval_matrix`.
These functions return interval enclosures of the true result.

These functions are presented below in a more general setting where the inputs are intervals as well.

So, to get the image by f of fixed floating-point values, simply create degenerated intervals in the next examples.

The interval evaluation of f is the image of the given input interval vector [x] by f, this range being noted by f([x]):

\[f([x]) := \{ f(x), x\in[x] \}.\]

Let us start with a real-valued function f with scalar arguments:

```
Variable x;
Variable y;
Function f(x,y,sin(x+y));
double _x[2][2]={{1,2},{3,4}};
IntervalVector xy(2,_x); // build xy=([1,2],[3,4])
Interval z=f.eval(xy); // z=f(xy)=sin([4,6])=[-1, -0.27941]
```

The sine function is not monotonic on [4,6] and actually reaches its minimum at \(3\pi/2\).

Note that the `eval` takes an `IntervalVector`
as argument, even if there is only one variable. So, in the latter case, you have to build a vector
reduced to a single component.

We consider now a vector-valued function.
Since the return type of an evaluation is not anymore an `Interval` but
an `IntervalVector`, we have to use a method with a different
signature, namely, `eval_vector`:

```
Variable a;
Function f(a,Return(sqr(a),-a));
IntervalVector x(1,Interval(1,2)); // build x=([1,2])
/* calculate y=f(x)=([1, 4] ; [-2, -1]) */
IntervalVector y=f.eval_vector(x);
```

Finally, for a matrix-valued function, the evaluation is obtained via `eval_matrix`.
We assume again that the following matrix-valued function

\[f: x \mapsto ( (2x, -x) ; (-x,3x) )\]

has been written in a “minibex” input file (see above).

```
Function f("myfunction.txt");
IntervalVector x(1,Interval(0,1));
// calculate M=f(x)=(([0, 2] , [-1, -0]) ; ([-1, -0] , [0, 3]))
IntervalMatrix M=f.eval_matrix(x);
```

For a scalar-valued function, you can get an interval enclosure of the gradient:

```
Variable x,y,z;
Function f(x,y,z,x*y+z*y);
double _xyz[3][2]={{0,1},{0,2},{0,3}};
IntervalVector xyz(3,_xyz);
/* calculate g=grad_f(x)=(y,x+z,y)=[0, 2] ; [0, 4] ; [0, 2]) */
IntervalVector g=f.gradient(xyz);
```

For a vector-valued function, you can get an interval enclosure of the Jacobian matrix:

```
// ------------------------------------------------
// Vector-valued functions, Jacobian matrix
//
// > create the function dist:(x,pt)->||x-pt||
// > create the function f:x->(dist(x,pt1),dist(x,pt2)
// > calculate the Jacobian matrix of f over the box
// ------------------------------------------------
Variable x(2,"x");
Variable pt(2,"p");
Function dist(x,pt,sqrt(sqr(x[0]-pt[0])+sqr(x[1]-pt[1])),"dist");
Vector pt1=Vector::zeros(2);
Vector pt2=Vector::ones(2);
Function f(x,Return(dist(x,pt1),dist(x,pt2)));
double init_box[][2] = { {-10,10},{-10,10} };
IntervalVector box(2,init_box);
/* calculate J as a m*n interval enclosure of the Jacobian matrix */
IntervalMatrix J=f.jacobian(box);
// ------------------------------------------------
```

One of the main feature of Ibex is the ability to *contract* a box representing the domain of a variable
x with respect to the constraint that f(x) belongs to a restricted input range [y].
Rigorously, given two intervals [x] and [y], the contraction gives a new interval [z] such that

\[\forall x\in[x], \quad f(x)\in[y] \Longrightarrow x \in[z] \subseteq [x]\]

One way to do this is by using the famous *backward* algorithm.
This algorithm does not return a new interval [z] but contract the input interval [x] which is therefore
an input-output argument.

In the following snippet we require the function sin(x+y) to take the value -1 (a degenerated interval). With an initial box (x,y)=([1,2],[3,4]), we obtain the result that (x,y) must lie in the subdomain ([1, 1.7123] ; [3, 3.7124]).

```
Variable x;
Variable y;
Function f(x,y,sin(x+y));
double _box[2][2]={{1,2},{3,4}};
IntervalVector box(2,_box);
/* the backward sets box to ([1, 1.7123] ; [3, 3.7124]) */
f.backward(-1.0,box);
```

To create a constraint, you can also either use strings or C++ objects:

With strings:

```
NumConstraint c("x","y","z","x+y<=z");
```

With C++ objects:

```
Variable x,y,z;
NumConstraint c(x,y,z,x+y<=z);
```

You can also refer to a previously defined function f to create, e.g., f(x)<=0:

```
Variable x,y,z;
Function f(x,y,z,x+y-z);
NumConstraint c(f,LEQ);
```

The key idea behind *contractor programming* *[Chabert & Jaulin, 2009]* is to abstract the algorithm
from the underlying constraint and to view it a function “C”:

\[C: \mathbb{IR}^n \to \mathbb{IR}^n \ \mbox{such that} \ C([x])\subseteq[x],\]

where \(\mathbb{IR}\) denotes the set of real intervals.

In other word, we take as primary concept the *operational* definition of a constraint.

In this way, operators (like the intersection and the others below) can be extended to contractors.

Since contractors implicitly represent sets, the fundamental advantage of extending operations to contractors is that we actually extend these operations to sets.

All contractors in Ibex are algorithms represented by different classes.
See the strategy pattern for more information on this design choice.
Classes representing contractors are prefixed by `Ctc`.

The standard way to contract with respect to a constraint is by using the *forward-bacwkard* algorithm.
The corresponding class is `CtcFwdBwd`.

A constraint has to be built first using the `NumConstraint` class.
In the following piece of code, we build a forward-backward contractor with respect to x+y=z.

```
Variable x,y,z;
NumConstraint c(x,y,z,x+y=z);
CtcFwdBwd ctc(c);
```

Of course, the expression of a constraint can involve a previously defined function.
Furthermore, if the constraint is simply “f=0”, where f is a `Function` object, it is not
necessary in this case to build an intermediate `NumConstraint` object.
One can directly give the function f that has to be nullify to `CtcFwdBwd`.
In the next example, we consider the problem of finding the point which distance from
both (0,0) and (1,1) is sqrt(2)/2. The solution is (0.5,0.5).

```
Variable x,y;
double d=0.5*sqrt(2);
Function f(x,y,Return(sqrt(sqr(x)+sqr(y))-d, sqrt(sqr(x-1.0)+sqr(y-1.0))-d));
IntervalVector box(2,Interval(-10,10));
/* we give f directly (means that the constraint is f=0) */
CtcFwdBwd c(f);
c.contract(box);
/* display ([0.2929, 0.7072] ; [0.2929, 0.7072]) */
cout << box << endl;
```

Of course, the result is rather crude. Remember that the purpose of `CtcFwdBwd`
is to contract *quickly* with respect to *any* numerical constraint: it is widely applicable and takes
a time that is only proportional to the expression size. In the other hand, it is not accurate in general.

The fixpoint operator applies a contractor C iteratively:

\[fixpoint(C): [x] \mapsto C(\ldots C([x])\ldots),\]

while the “gain” is more than the given `ratio`. More precisely, the “gain”
is the relative Hausdorff distance between the input box [x] and the output box C([x]) but, often, you can ignore
the precise meaning of this gain and just consider that the procedure will loop until the contracted box
will roughly differ “by ratio” from the input one.

Let us now follow the previous example. As said, the solution is (0.5,0.5). We can see that simply embedding
the `CtcFwdBwd` contractor in a fixpoint loop (with a `ratio`
set to 0.1) gives a box with sharp bounds.

```
Variable x,y;
double d=0.5*sqrt(2);
Function f(x,y,Return(sqrt(sqr(x)+sqr(y))-d, sqrt(sqr(x-1.0)+sqr(y-1.0))-d));
IntervalVector box(2,Interval(-10,10));
CtcFwdBwd c(f);
CtcFixPoint fp(c,1e-03);
fp.contract(box);
/* display ([0.4990, 0.5001] ; [0.4990, 0.5001]) */
cout << "box after fixpoint=" << box << endl;
}
{
Variable x;
NumConstraint c1(x,x>=-1);
NumConstraint c2(x,x<=1);
CtcFwdBwd ctc1(c1);
CtcFwdBwd ctc2(c2);
IntervalVector box(1,Interval::ALL_REALS);
CtcCompo ctc3(ctc1,ctc2);
ctc3.contract(box);
cout << "compo: " << box << endl;
```

Given two of more contractors, we can apply the two logical operators *union* and *intersection*:

\[union(C_1,\ldots,C_n): [x] \mapsto C_1([x]) \cup\ldots\cup C_n([x]).\]\[inter(C_1,\ldots,C_n): [x] \mapsto C_1([x]) \cap\ldots\cap C_n([x]).\]

However, the latter operation is barely used and usually replaced by the *composition*:

\[compo(C_1,\ldots,C_n): [x] \mapsto C_n(\ldots(C_1([x])\ldots).\]

Indeed, one can see that the composition amounts to the same logical operation (the intersection of each contractor’s set), but in a more efficient way since we take advantage of the contraction performed by \(C_1,...,C_{i-1}\) when contracting with \(C_i\). In contrast, the intersection operator calls each contractor independently on the same initial box.

The corresponding classes are `CtcUnion` and `CtcCompo`.

As a rule of thumb, use `CtcUnion` for the union of two contractors and `CtcComp` for the intersection.
Here is an example with the union:

```
Variable x;
NumConstraint c1(x,x<=-1);
NumConstraint c2(x,x>=1);
CtcFwdBwd ctc1(c1);
CtcFwdBwd ctc2(c2);
IntervalVector box(1,Interval::POS_REALS); // the box [0,oo)
CtcUnion ctc3(ctc1,ctc2); // a contractor w.r.t. (x<=-1 or x>=1)
ctc3.contract(box); // box will be contracted to [1,oo)
cout << box << endl;
```

Here is an example with the intersection (composition):

```
Variable x;
NumConstraint c1(x,x>=-1);
NumConstraint c2(x,x<=1);
CtcFwdBwd ctc1(c1);
CtcFwdBwd ctc2(c2);
IntervalVector box(1,Interval::ALL_REALS);// the box (-oo,oo)
CtcCompo ctc3(ctc1,ctc2);// a contractor w.r.t. (x>=-1 and x<=1)
ctc3.contract(box); // box will be contracted to [-1,1]
cout << box << endl;
```

When a function is “square” (the dimension is the same as the codimension, i.e., \(f:\mathbb{R}^n\to\mathbb{R}^n\)), you can contract a box with respect to the constraint f(x)=0 using the interval Newton iteration.

You just have to build a `CtcNewton` object with the function and call `contract`.

This operator can give extremly accurate bounds proving that the input box is already “sufficiently” small (that is, “inside the convergence basin” of Newton’s iteration). In the following example, we give a box that encloses the solution (1,0) with a radius of 10^-3. Newton’s iteration contracts this box downto the maximal precision:

```
Variable x,y;
Function f(x,y,Return(sqrt(sqr(x)+sqr(y))-1, sqrt(sqr(x-1.0)+sqr(y-1.0))-1));
double init_box[][2]={{0.999,1.001},{-0.001,0.001}};
IntervalVector box(2,init_box);
// Build an interval Newton iteration
// for solving f(x)=0 where f is
// a vector-valued function representing
// the system.
CtcNewton newton(f);
/* Contract the box with Newton */
newton.contract(box);
/* display a very small box enclosing (1,0) */
cout << box << endl;
```

The *propagation* operator calculates the fixpoint of (the composition of) n contractors by using a more
sophisticated (“incremental”) strategy than a simple loop.
So, semantically, the propagation operator can be defined as follows:

\[propagation(C_1,\ldots,C_n):=fixpoint(compo(C_1,\ldots,C_n)).\]

(see above for the definition of the fixpoint and composition operators).

The key idea behind this operator is to avoid calling contractors that will certainly leave the box intact. Contractors that can potentially enforce a contraction are determined typically from the syntax of their underlying constraint. Consider for instance two contractors, C1 w.r.t. f(x,z)=0 and C2 w.r.t. g(x,y)=0. Assume that the fixpoint for C1 is reached with the current box ([x],[y],[z]). If a call to C2 only contracts the second interval (the one corresponding to y), it is then useless to call C1 again.

So, by using such principle, the propagation calculates the fixpoint by “awaking” contractors only when necessary. Of course, the more sparse the constraint system, the more valuable the propagation, when compared to a simple fixpoint.

The following example compares the number of contractors

The Q-intersection is typically used in a context where we have a set of contractors that result from measurements (each measurement enforces a constraint), some of which can be incorrect.

If we are sure that at least q measurements are correct (which amounts to say that the number of outliers is bounded by N-q) then we can contract the box in a robust way, by calculating the union of the boxes resulting from the contraction with all combinaisons of q contractors among N.

Mathematicaly, with \((i_1 , . . . , i_q)\) ranging over the set of all q distinct indices between 0 and N-1:

\[q-inter(C_1,\ldots,C_n,q):=union(\ldots,inter(Ci1,\ldots,Ciq),\ldots)\]

Here is a simple example inspired from parameter estimation.

We assume a point (x,y) has to be localized. We measure 4 distances “bD” from 6 (approximately known) points (bX,bY). Each position bX, bY and each distance bD has an uncertainty [-0.1,0.1]. We also know there may be at most one outlier.

The solution point is: x=6.32193 y=5.49908

First of all, let us enter the coordinates of the points (bX,bY) and the distances. This data will simulate our measurements.

```
const int N=6;
/* The measurements (coordinates of the points and distances) */
double bx[N]={5.09392,4.51835,0.76443,7.6879,0.823486,1.70958};
double by[N]={0.640775,7.25862,0.417032,8.74453,3.48106,4.42533};
double bd[N]={5.0111,2.5197,7.5308,3.52119,5.85707,4.73568};
```

define the measurement intervals (with uncertainty taken into account)

```
Interval bX[N];
Interval bY[N];
Interval bD[N];
/* add uncertainty on measurements */
for (int i=0; i<N; i++) {
bX[i]=bx[i]+Interval(-0.1,0.1);
bY[i]=by[i]+Interval(-0.1,0.1);
bD[i]=bd[i]+Interval(-0.1,0.1);
}
```

Now, we artificially introduce an outlier by shifting the interval for one measurement (here, x position n°5) by a large value:

```
bX[5]+=10;
```

Now, all our simulated data is set up. We just have to define the contractors. We first declare the distance function and then 6 contractors corresponding to the distance with each (bX,bY):

```
Variable x(2);
Variable px,py;
Function dist(x,px,py,sqrt(sqr(x[0]-px)+sqr(x[1]-py)));
Function f0(x,dist(x,bX[0],bY[0])-bD[0]);
Function f1(x,dist(x,bX[1],bY[1])-bD[1]);
Function f2(x,dist(x,bX[2],bY[2])-bD[2]);
Function f3(x,dist(x,bX[3],bY[3])-bD[3]);
Function f4(x,dist(x,bX[4],bY[4])-bD[4]);
Function f5(x,dist(x,bX[5],bY[5])-bD[5]);
CtcFwdBwd c0(f0);
CtcFwdBwd c1(f1);
CtcFwdBwd c2(f2);
CtcFwdBwd c3(f3);
CtcFwdBwd c4(f4);
CtcFwdBwd c5(f5);
```

We can contract now a box with the q-intersection of these contractors:

```
/* The initial box: [0,10]x[0,10] */
IntervalVector initbox(2,Interval(0,10));
/* Create the array of all the contractors */
Array<Ctc> array(c0,c1,c2,c3,c4,c5);
/* Create the q-intersection of the N contractors */
CtcQInter q(array,5); // 2 is the number of variables, 5 the number of correct measurement
/* Perform a first contraction */
IntervalVector box=initbox;
q.contract(box);
cout << "after q-inter =" << box << endl;
```

The displayed result is ([3.9667, 7.2381] ; [4.5389, 8.1479]). Of course, we can do better by calculating a fixpoint of the q-intersection:

```
/* Build a Fix-point of the q-intersection */
CtcFixPoint fix(q);
/* Perform a stronger contraction with the fixpoint */
fix.contract(box);
cout << "after fix+q-inter =" << box << endl;
```

The displayed result is ([5.9277, 6.8836] ; [5.0914, 5.7996]) which, indeed, better encloses the solution point x=6.32193 y=5.49908.

To create a contractor, you just have to
- declare a class that extends `Ctc`
- create inside a function `contract` that takes a reference to a box (`IntervalVector&`) and contracts it. The function returns `void`.

In the following example, we create a contractor that simply divides by two the radius of each component.

```
class MyContractor : public Ctc {
public:
MyContractor(int nb_var) : Ctc(nb_var) {}
void contract(IntervalVector& box) {
box=box.mid()+0.5*Interval(-1,1)*box.rad();
}
};
```

Then, if we create this contractor and applies it several time to the same box, we can observe the expected result:

```
/* build the contractor for 3-dimensional boxes. */
MyContractor c(3);
/* create the box [0,1]x[0,1]x[0,1] */
IntervalVector x(3,Interval(0,1));
c.contract(x);
cout << x << endl;// ([0.25, 0.75] ; [0.25, 0.75] ; [0.25, 0.75])
c.contract(x);
cout << x << endl;// ([0.375, 0.625] ; [0.375, 0.625] ; [0.375, 0.625])
c.contract(x);
cout << x << endl;// ([0.4375, 0.5625] ; [0.4375, 0.5625] ; [0.4375, 0.5625])
```

This contractor can now be combined with the ones built-in. For instance, we can decide to calculate the fixpoint. Then, the result is a small box enclosing (0.5,0.5,0.5):

```
CtcFixPoint fp(c,0.001);
fp.contract(x);
cout << x << endl;// ([0.4999999999999999, 0.5000000000000001], ...
```

The generic solver is a classical branch and prune algorithm that interleaves contraction and branching (bisection) until boxes get sufficiently small. It takes three operators in input:

a

**contractor**Operator that contracts boxes by removing non-solution points. It should be emphasized that the solver only sees the contraction operator, not the underlying system (of equations or inequalities). This is one of the reason why it is said to be “generic”. In particular, one can introduce contractors of other kind, that are not based on mathematical expressions.

a

**bisector**Operator that splits a box. The box is bisected providing it is large enough, so this operator takes also a

*precision*parameter. If the box is too small, the solver will not continue the search and backtrack. See*Bisectors*for more details.a

**cell buffer**Operator that manages the list of pending boxes (a

*cell*is a box with a little bit of extra information used by the search). See*Cell buffers*for more details.

Our next example creates a solver for the intersection of two circles of radius \(d\), one centered on \((0,0)\) and the other in \((1,0)\).

To this end we first create a vector-valued function:

\[\begin{split}(x,y) \mapsto \begin{pmatrix} x^2+y^2-d \\ (x-1)^2+y^2-d \end{pmatrix}\end{split}\]

Then, we build two contractors; a forward-backward contractor and (because the system is square), an interval Newton contractor.

We chose as bisection operator the *round-robin* operator, that splits each component in turn.
The precision of the solver is set to 1e-7.

Finally, the cell buffer is a stack, which leads to a depth-first search.

Warning

This solver is simple and not very efficient. The *IbexOpt* plugin comes with a
default solver that implements a more sophisticated strategy.

```
/* Create the function (x,y)->( ||(x,y)||-d, ||(x,y)-(1,0)||-d ) */
Variable x,y;
double d=1.0;
Function f(x,y,Return(sqrt(sqr(x)+sqr(y))-d,
sqrt(sqr(x-1.0)+sqr(y))-d));
double init_box[][2] = { {-10,10},{-10,10} };
IntervalVector box(2,init_box);
/* Create a first contractor w.r.t f(x,y)=0 (forward-backward) */
CtcFwdBwd fwdBwd(f);
/* Create a second contractor (interval Newton) */
CtcNewton newton(f);
/* Compose the two contractors */
CtcCompo compo(fwdBwd,newton);
/* Create a round-robin bisection heuristic and set the
* precision of boxes to 1e-07 */
RoundRobin bisector(1e-07);
/* Create a "stack of boxes" (CellStack), which has the effect of
* performing a depth-first search. */
CellStack buff;
/* Create a solver with the previous objects */
Solver s(compo,bisector,buff);
/* Run the solver */
vector<IntervalVector> sols=s.solve(box);
/* Display the solutions */
for (unsigned int i=0; i<sols.size(); i++)
cout << "solution n°" << i << "=\t" << sols[i] << endl;
```