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 IntervalVectorinthiscaseisthedimension(here,3),thesecondthearrayof``double.

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:

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).

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:

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)\):

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.

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.

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

With strings:

Functionf("x","(x-1,x+1)");

With operator overloading:

Variablex;Functionf(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):

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:

functionf(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:

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:

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

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]).

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.

Variablex,y,z;NumConstraintc(x,y,z,x+y=z);CtcFwdBwdctc(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).

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.

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.

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:

Variablex;NumConstraintc1(x,x<=-1);NumConstraintc2(x,x>=1);CtcFwdBwdctc1(c1);CtcFwdBwdctc2(c2);IntervalVectorbox(1,Interval::POS_REALS);// the box [0,oo)CtcUnionctc3(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):

Variablex;NumConstraintc1(x,x>=-1);NumConstraintc2(x,x<=1);CtcFwdBwdctc1(c1);CtcFwdBwdctc2(c2);IntervalVectorbox(1,Interval::ALL_REALS);// the box (-oo,oo)CtcCompoctc3(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:

Variablex,y;Functionf(x,y,Return(sqrt(sqr(x)+sqr(y))-1,sqrt(sqr(x-1.0)+sqr(y-1.0))-1));doubleinit_box[][2]={{0.999,1.001},{-0.001,0.001}};IntervalVectorbox(2,init_box);// Build an interval Newton iteration// for solving f(x)=0 where f is// a vector-valued function representing// the system.CtcNewtonnewton(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 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:

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.

constintN=6;/* The measurements (coordinates of the points and distances) */doublebx[N]={5.09392,4.51835,0.76443,7.6879,0.823486,1.70958};doubleby[N]={0.640775,7.25862,0.417032,8.74453,3.48106,4.42533};doublebd[N]={5.0111,2.5197,7.5308,3.52119,5.85707,4.73568};

define the measurement intervals (with uncertainty taken into account)

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

/* The initial box: [0,10]x[0,10] */IntervalVectorinitbox(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 */CtcQInterq(array,5);// 2 is the number of variables, 5 the number of correct measurement/* Perform a first contraction */IntervalVectorbox=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 */CtcFixPointfix(q);/* Perform a stronger contraction with the fixpoint */fix.contract(box);cout<<"after fix+q-inter ="<<box<<endl;

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.