# Solvers¶

This documentation is for using and building solvers in C++.

You should first know how solvers work in Ibex. Read for this the user guide.

## Calling IbexSolve from C++¶

You can call IbexSolve (the default solver) and get the solutions from C++.

Two objects must be built: the first represents the problem (namely, a system), the second the solver itself. Then, we just run the solver. Here is a simple example:

```
/* Build a system of equations from the file */
System system(IBEX_BENCHS_DIR "/others/kolev36.bch");
/* Build a default solver for the system and with a precision set to 1e-07 */
DefaultSolver solver(system,1e-07);
solver.solve(system.box); // Run the solver
/* Display the solutions. */
output << solver.get_data() << endl;
```

The output is:

```
solution n°1 = ([0.1173165676349099, 0.1173165676349106] ; [0.4999999999999989, 0.5000000000000014] ; [0.8826834323650891, 0.8826834323650912] ; [-0.2071067811865488, -0.2071067811865466] ; [1.207106781186544, 1.207106781186553] ; [-2.000000000000004, -1.999999999999998])
```

## The generic solver¶

The generic solver is the main C++ class behind the implementation of IbexSolve. It is a classical branch and prune algorithm that interleaves contraction and branching (bisection) until boxes get sufficiently small. However, it performs a more general task that just finding solution points of square systems of equations: it also knows how to deal with under-constrained systems and handle manifolds.

Compared to IbexSolve, the generic solver allows the following additional operators as inputs:

a

**contractor**Operator that contracts boxes by removing non-solution points. The contraction operator must be compatible with the system given (equations/inequalities). The solver performs no check (it is the user responsability). See Contractors.

a

**bisector**Operator that splits a box. Note that some bisectors have a

*precision*parameter: the box is bisected providing it is large enough. But this precision is not directly seen by the solver which has its own precision variables (see`-e`̀ and ``-E`

). If however the bisector does not split a box, this will generate an exception caught by the solver, which will not continue the search and backtrack. So fixing the bisector internal precision gives basically the same effect as fixing it with`--e`

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

Then, we build two contractors; a forward-bacwkard 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.

```
/* 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));
/* Create the system f(x,y)=0. */
SystemFactory factory;
factory.add_var(x);
factory.add_var(y);
factory.add_ctr(f(x,y)=0);
System system(factory);
/* Create the domain of variables */
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 0. */
RoundRobin bisector(0);
/* Create a "stack of boxes" (CellStack), which has the effect of
* performing a depth-first search. */
CellStack buff;
/* Vector precisions required on variables */
Vector prec(6, 1e-07);
/* Create a solver with the previous objects */
Solver s(system, compo, bisector, buff, prec, prec);
/* Run the solver */
s.solve(box);
/* Display the solutions */
output << s.get_data() << endl;
```

The output is:

```
solution n°1 = ([0.4999999999999996, 0.5000000000000003] ; [-0.8660254037844389, -0.8660254037844383])
solution n°2 = ([0.4999999999999998, 0.5000000000000005] ; [0.8660254037844383, 0.8660254037844389])
```

## Implementing IbexSolve (the default solver)¶

IbexSolve is an instance of the generic solver with (almost) all parameters set by default.

We already showed how to Calling IbexSolve from C++. To give a further insight into the generic solver and its possible settings, we explain now how to re-create the default solver by yourself.

The contractor of the default solver is obtained with the following receipe. This is a composition of

- HC4
- ACID
- Interval Newton (only if it is a square system of equations)
- A fixpoint of the Polytope Hull of two linear relaxations combined:
- the relaxation called X-Taylor;
- the relaxation generated by affine arithmetic. See Linearizations.

The bisector is based on the The Smear Function with maximal relative impact.

So the following program exactly reproduces the default solver.

```
System system(IBEX_BENCHS_DIR "/others/kolev36.bch");
/* ============================ building contractors ========================= */
CtcHC4 hc4(system,0.01);
CtcHC4 hc4_2(system,0.1,true);
CtcAcid acid(system, hc4_2);
CtcNewton newton(system.f_ctrs, 5e+08, 1e-07, 1e-04);
LinearizerXTaylor linear_relax(system);
CtcPolytopeHull polytope(linear_relax);
CtcCompo polytope_hc4(polytope, hc4);
CtcFixPoint fixpoint(polytope_hc4);
CtcCompo compo(hc4,acid,newton,fixpoint);
/* =========================================================================== */
/* Create a smear-function bisection heuristic. */
SmearSumRelative bisector(system, 1e-07);
/* Create a "stack of boxes" (CellStack) (depth-first search). */
CellStack buff;
/* Vector precisions required on variables */
Vector prec(6, 1e-07);
/* Create a solver with the previous objects */
Solver s(system, compo, bisector, buff, prec, prec);
/* Run the solver */
s.solve(system.box);
/* Display the solutions */
output << s.get_data() << endl;
/* Report performances */
output << "cpu time used=" << s.get_time() << "s."<< endl;
output << "number of cells=" << s.get_nb_cells() << endl;
```

## Parallelizing search¶

It is possible to parallelize the search by running (in parallel) solvers for different subboxes of the initial box.

Be aware however that Ibex has not been designed (so far) to be parallelized and the following lines only reports our preliminary experiments.

Here are the important observations:

The sub-library gaol is

**not**thread-safe. You must compile Ibex with**filib**which seems to be OK (see Configuration options).The linear solver Soplex (we have not tested yet with Cplex) seems to be thread-safe but sometimes generates error messages on the console like:

ISOLVE56 stop: 0, basis status: PRIMAL (2), solver status: RUNNING (-1)

So, calling Soplex several times simultaneously seems not to be allowed, but Soplex at least manages the case properly, that is, stops. As far as we have observed, we don’t lose solutions even when this kind of message appear.

Ibex objects are not thread-safe which means that the solvers run in parallel must share no information. In particular, each solver must have its

**own copy**of the system.

Here is an example:

```
// Get the system
System sys1(IBEX_BENCHS_DIR "/polynom/ponts-geo.bch");
// Create a copy for the second solver
System sys2(sys1,System::COPY);
// Precision of the solution boxes
double prec=1e-08;
// Create two solvers
DefaultSolver solver1(sys1,prec);
DefaultSolver solver2(sys2,prec);
// Create a partition of the initial box into two subboxes,
// by bisecting any variable (here, n°4)
pair<IntervalVector,IntervalVector> pair=sys1.box.bisect(4);
// =======================================================
// Run the solvers in parallel
// =======================================================
#pragma omp parallel sections
{
solver1.solve(pair.first);
#pragma omp section
solver2.solve(pair.second);
}
// =======================================================
output << "solver #1 found " << solver1.get_data().size() << endl;
output << "solver #2 found " << solver2.get_data().size() << endl;
```

If I remove the `#pragma`

the program displays:

```
solver #1 found 64
solver #2 found 64
real 0m5.121s // <-------- total time
user 0m5.088s
```

With the `#pragma`

, I obtain:

```
solver #1 found 64
solver #2 found 64
real 0m2.902s // <-------- total time
user 0m5.468s
```

**Note:** It is pure luck that by bisecting the 4th variable, we obtain exactly half of the solutions on each sub-box.
Also, looking for the 64 first solutions takes here around the same time than looking for the 64 subsequent ones, which
is particular to this example. So, contrary to what this example seems to prove, splitting the box in two subboxes does
not divide the running time by two in general. Of course :)

If you are afraid about the messages of the linear solver, you can replace the `DefaultSolver`

by your own dedicated solver
that does not resort to the simplex, ex:

```
Vector eps_min(sys1.nb_var,prec);
Vector eps_max(sys1.nb_var,POS_INFINITY);
Solver solver1(sys1,*new CtcCompo(*new CtcHC4(sys1),*new CtcNewton(sys1.f_ctrs)), *new RoundRobin(prec), *new CellStack(), eps_min, eps_max);
Solver solver2(sys2,*new CtcCompo(*new CtcHC4(sys2),*new CtcNewton(sys2.f_ctrs)), *new RoundRobin(prec), *new CellStack(), eps_min, eps_max);
```