Tutorial¶
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 IbexSolve or IbexOpt plugins.
Basic Interval computations¶
Start a program¶
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.
Creating intervals¶
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
Operation between intervals¶
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 (-, *, /).
Applying a function to an interval¶
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...]
Interval vectors¶
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¶
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]))
Operations between matrices and vectors¶
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
Midpoint, radius, magnitude, etc.¶
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)
Functions¶
Creating functions¶
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)
Constants inside functions¶
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)
Functions with vector arguments¶
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.
Composing functions¶
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>)))
Vector-valued functions¶
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):
// ------------------------------------------------
// 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.
Matrix-valued functions¶
You can also create functions that return matrices. Here is an example of a function from \(R\) to \(R^{2\times 2}\) where:
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,ExprVector::ROW),Return(-x,3*x,ExprVector::ROW)));
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.
Using the Minibex syntax¶
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");
Minibex syntax with intermediate variables¶
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
Evaluation over floating-point numbers¶
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.
Interval evaluation¶
The interval evaluation of f is the image of the given input interval vector [x] by f, this range being noted by f([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
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);
Interval gradient¶
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);
Interval Jacobian matrix¶
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);
// ------------------------------------------------
Backward (or contraction)¶
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]).
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);
Constraints¶
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);
Contractors¶
What is a contractor programming?¶
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”:
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.
Forward-Backward¶
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.
Fixpoint¶
The fixpoint operator applies a contractor C iteratively:
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;
Intersection, union & composition¶
Given two of more contractors, we can apply the two logical operators union and intersection:
However, the latter operation is barely used and usually replaced by the composition:
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;
Interval Newton¶
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;
Propagation¶
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:
(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.
Q-Intersection (robustness w.r.t. outliers)¶
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:
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.
Build your own contractor¶
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], ...