IbexSolve

Getting started

IbexSolve is a end-user program that solves a system of nonlinear equations rigorously (that is, it does not lose any solution and return each solution under the form of a small box enclosing the true value). It resorts to a unique black-box strategy (whatever the input problem is) and with a very limited number of parameters. Needless to say, this strategy is a kind of compromise and not the best one for a given problem.

Note that this program is based on the generic solver, a C++ class that allows to build a more customizable solver.

You can directly apply this solver on one of the benchmark problems distributed with Ibex. The benchmarks are all written in the Minibex syntax and stored in an arborescence under plugins/solver/benchs/.

Open a terminal, move to the bin subfolder and run IbexSolve with, for example, the problem named kolev36 located at the specified path:

~/ibex-2.8.0$ cd bin
~/ibex-2.8.0/bin$ ./ibexsolve ../plugins/solver/benchs/others/kolev36.bch

The following result should be displayed:

***************************** setup *****************************
  file loaded:         ../plugins/solver/benchs/others/kolev36.bch
  output file:         ../plugins/solver/benchs/others/kolev36.mnf
*****************************************************************

running............

solving successful!

number of inner boxes:         1
number of boundary boxes:      0
number of unknown boxes:       0
number of pending boxes:       0
cpu time used:                 0.0640001s
number of cells:               25

results written in ../plugins/solver/benchs/others/kolev36.mnf

The number of “inner boxes” correspond to the number of solutions found (there is just one here). To see the solution, use the option -s.

In the report, the “number of cells” correspond to the number of hypothesis (bisections) that was required to solve the problem.

Options

-e<float>, –eps-min=<float>

Minimal width of output boxes. This is a criterion to stop bisection: a non-validated box will not be larger than ‘eps-min’. Default value is 1e-3.

-E<float>, –eps-max=<float>

Maximal width of output boxes. This is a criterion to force bisection: a validated box will not be larger than ‘eps-max’ (unless there is no equality and it is fully inside inequalities). Default value is +oo (none)

-t<float>, –timeout=<float>

Timeout (time in seconds). Default value is +oo (none).

-i<filename>, –input=<filename>

COV input file. The file contains a (intermediate) description of the manifold with boxes in the COV (binary) format.

-o<filename>, –output=<filename>

COV output file. The file will contain the description of the manifold with boxes in the COV (binary) format.

–simpl=…

Expression simplification level. Possible values are:

  • 0: no simplification at all (fast).

  • 1: basic simplifications (fairly fast).

    E.g. x+1+1 –> x+2

  • 2: more advanced simplifications without developing (can be slow).

    E.g. x*x + x^2 –> 2x^2. Note that the DAG structure can be lost.

  • 3: simplifications with full polynomial developing (can blow up!).

    E.g. x*(x-1) + x –> x^2. Note that the DAG structure can be lost.

Default value is : 1.

-s, –sols

Display the “solutions” (output boxes) on the standard output.

–bfs

Perform breadth-first search (instead of depth-first search, by default)

–trace

Activate trace. “Solutions” (output boxes) are displayed as and when they are found.

–boundary=…

Boundary test strength. Possible values are:

  • true: always satisfied. Set by default for under constrained problems (0<m<n).

  • full-rank: the gradients of all constraints (equalities and potentially activated inequalities) must be linearly independent.

  • half-ball: (not implemented yet) the intersection of the box and the solution set is homeomorphic to a half-ball of R^n

  • false: never satisfied. Set by default if m=0 or m=n (inequalities only/square systems)

–random-seed=<float>

Random seed (useful for reproducibility). Default value is 1.

-q, –quiet

Print no report on the standard output.

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.

Note

A more detailed documentation about under-constrained systems will be available soon.

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

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

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

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

  1. HC4

  2. ACID

  3. Interval Newton (only if it is a square system of equations)

  4. 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;