# Do it Yourself!¶

The examples in this page are presented under the form of “labs” so that you can use them for practicing.

The complete source codes are available under the `examples/`

folder.

In these labs, we use Vibes to plot boxes but you can easily adapt the code to use your own graphical tool.

Fast instructions for installing and using Vibes are given here.

## Set image¶

The complete code can be found here: `examples/lab/lab1.cpp`

.

**Introduction**

The goal of this first lab is to calculate the image of a box by a function using interval arithmetic:

If we denote by f the function and [x] the initial box, then the set S to calculate is:

Applying directly an interval evaluation of the function f to [x] will give a single box that only represents an enclosure of S.

To fight with the wrapping effect, we will split [x] into smaller boxes and evaluate the function with every little box as argument. This will result in a better description of the set that will eventually converge to S as the size of the boxes tend to zero.

This will consist in three tasks:

- creating the function f
- creating the initial box [x] and splitting it into small boxes
- evaluating the function of every boxes
- plotting the results

**Question 1**

Create in the `main`

the function

**Question 2**

Create the box ([x],[y])=([0,6],[0,6]) and split each dimension into n slices, where n is a constant.

show/hide solution**Question 3**

Evaluate the function on each box and plot the result with Vibes.

show/hide solution**Question 4**

Compare the result with n=15, n=80 and n=500.

You should obtain the following pictures:

## Set inversion (basic)¶

The complete code can be found here: `examples/lab/lab2.cpp`

.

**Introduction**

The goal of this lab is to program **Sivia** (*set inversion with interval analysis*) [Jaulin & Walter 1993] [Jaulin 2001], an algorithm that draws a paving
representing a set E defined implicitely as the preimage of an interval [z] by a non-linear function \(f:\mathbb{R}^n\to\mathbb{R}\) (here n=2).

The Sivia algorithm performs a recursive exploration of the initial box by applying the following steps:

**inner test**: if the image of ([x],[y]) by f is a subset of [z], the box is painted in green;**outer test**: if the image does not intersect [z], the box is painted in red;- if none of these test succeeds and if ([x],[y]) has a maximal diameter greater than \(\varepsilon\), the box is split and the procedure is recursively called on the two subboxes.

**Question 1 (Initialisation)**

Create the `Function`

object that represents

and the initial bounding box ([-10,10],[-10,10]).

show/hide solution**Question 2 (Initialisation)**

We shall use a stack for implementing the recursivity. This stack is a container that will be used to store boxes.

Create a C++ stack and set the precision of bisection to 0.1.

Push the initial box in the stack. Define the image interval [z] and initialize it to [0,2].

show/hide solution**Question 3**

Create the loop that pop boxes from the stack until it is empty.
Define a local variable `box`

to be the current box (the one on top of the stack).

*Hint: use the top() and pop() functions of the* `stack`

*class*.

**Question 4**

Implement the inner test (see above).

*Hint*: use is_subset.

**Question 5**

Implement the outer test (see above).

*Hint*: use intersects.

**Question 6**

If none of these test succeeds, split the box. We will split the box on the axis of its largest size. Finally, the two subboxes are pushed on the stack.

*Hint*: use extr_diam_index and bisect.

## Set inversion (with contractors)¶

The complete code can be found here: `examples/lab/lab3.cpp`

.

**Introduction**

We will improve the **Sivia** algorithm by replacing in the loop the inner and outer tests by contractions.
This leads to a more compact paving and a smaller number of boxes (see figure below).

The first part of the code is unchanged:

```
int main() {
vibes::beginDrawing ();
vibes::newFigure("lab3");
// Create the function we want to apply SIVIA on.
Variable x,y;
Function f(x,y,sin(x+y)-0.1*x*y);
// Build the initial box
IntervalVector box(2);
box[0]=Interval(-10,10);
box[1]=Interval(-10,10);
// Create a stack (for depth-first search)
stack<IntervalVector> s;
// Precision (boxes of size less than eps are not processed)
double eps=0.1;
// Push the initial box in the stack
s.push(box);
...
```

The idea is to contract the current box either with respect to the constraint

in which case the contracted part will be painted in red, or

in which case the contracted part will be painted in green.

Given a contractor c, the contracted part is also called the *trace* of the contraction and is defined as \([x]\backslash c([x])\).

**Question 1**

Build forward-backward contractors for the four constraints (see the tutorial):

**Question 2**

Thanks to the composition, build a contractor w.r.t. \(f(x)\in[0,2]\).

Similarly, thanks to the union, build a contractor w.r.t. \(f(x)\not\in[0,2]\).

show/hide solution**Question 3**

Create the function `contract_and_draw`

with the following signature:

```
void contract_and_draw(Ctc& c, IntervalVector& box, const char* color);
```

This function must contract the box `box`

in argument with the contractor `c`

and plot the trace of the contraction (see above) with Vibes,
with the specified color `color`

.

*Hints: use the* diff *function of* `IntervalVector`

*to calculate the set difference between two boxes.*

**Question 4**

Replace in the loop the inner/outer tests by contractions.

show/hide solution## Set Inversion (using “Sets”)¶

The complete code can be found here: `examples/lab/lab4.cpp`

.

**Introduction**

The purpose of this exercice is just to get familiar with the structure proposed in Ibex for representing sets (or pavings).

The set inversion is naturally one of the main features proposed in this part of the library.
We will solve the same problem as before but this time with the `Set`

class directly.
This will take only a few lines of code.

Give first a look at the documentation on sets.

**Question 1**

Create the function \((x,y)\mapsto \sin(x+y)-0.1\times x\times y\) and the forward-backward separator associated to the constraint

**Question 2**

Build the initial set [-10,10]x[-10,10] and contract it with the separator with a precision of 0.1.

show/hide solution**Question 3**

Plot the set with Vibes using a `SetVisitor`

.

Solution: copy-paste the code given here.

## Parameter Estimation¶

The complete code can be found here: `examples/lab/lab5.cpp`

.

**Introduction**

This exercice is inspired by this video.

The problem is to find the values of two parameters \((p_1,p_2)\) of a physical process that are consistent with some measurements. Measurements are subject to error and we want a garanteed enclosure of all the feasible parameters values.

The physical process is modeled by a function \(f_{p_1,p_2}:t\mapsto y\) and a measurement is a couple of input-ouput \((t_i,y_i)\). We assume the input has no error. The error on the output is represented by an interval.

The model is:

We have the following series of measurements:

t | y |
---|---|

1 | [4.5,7.5] |

2 | [0.67,4.6] |

3 | [-1,2.8] |

4 | [-1.7,1.7] |

5 | [-1.9,0.93] |

6 | [-1.8,0.5] |

7 | [-1.6,0.24] |

8 | [-1.4,0.09] |

9 | [-1.2,0.0089] |

10 | [-1,-0.031] |

**Question 1**

Build the function f as a mapping of 3 variables, p1, p2 and t.

show/hide solution**Question 2**

Build two interval vectors `t`

and `y`

of size 10 that contain the measurements data (even if the input has no error, we will enter times as
degenerated intervals).

*Hint:* build interval vectors from array of double.

**Question 3**

Build a system using a system factory. The system must contain the 10 constraints that represent each measurements and the additional bound constraints on the parameters:

**Question 4**

Calculate the parameter values using set inversion (see lab n°4). You should obtain the following picture:

show/hide solution## Parameter Estimation (advanced)¶

The complete code can be found here: `examples/lab/lab6.cpp`

.

**Introduction**

This lab is a follow-up of the previous one.

We now introduce uncertainty on the input variable t. However, we will try somehow to make our parameter
estimation *robust* with respect to this uncertainty. This means that the values of our parameters should
be consistent with our output whatever is the actual value of the input. Mathematically, we require
p1 and p2 to respect the following constraints.

As a contractor-oriented library, Ibex does not provide quantifiers at the modeling stage. This means that you cannot write directly a constraint like this one. You have to build contractors and apply “quantifiers” on contractors. Read the documentation about contractors and quantifiers.

**Question 1**

Like in the previous lab, create the function and the vector of measurements.

show/hide solutionThen, define a constant `delta_t`

(the uncertainty on time) and inflate the
vector of representing input times by this constant:

**Question 2**

We will use the generic constructor of `CtcForAll`

.
Create a bitset that will indicate among the arguments of the function f which ones will be treated as variables
and which ones will be treated as quantified parameters.

**Question 3**

Create the inner and outer contractor for the “robust” parameter estimation problem.
Set the splitting precision \(\varepsilon\) of the parameter to `tdelta`

/5.

*Hints: follow the same idea as in the contractor variant of set inversion* (see lab n°3).
*Also, we need here to deal with n inner/outer contractors. To perform the composition/union of several
contractors, see* here.

**Question 4**

- Create a separator from the two contractors;
- Create a set from the box [0,1]x[0,1];
- Contract the set with the separator.

The pictures below show the results obtained for increasing values of `deltat`

(0,0.1,0.2 and 0.3).
As expected, the larger the error on input, the smaller the set of feasible
parameters.

deltat=0 | deltat=0.1 | deltat=0.2 | deltat=0.3 |

## Stability¶

The complete code can be found here: `examples/lab/lab7.cpp`

.

**Introduction**

The goal of this lab is to cast a classical problem in control theory into a set inversion problem.

We have a dynamical system y(t) governed by the following linear differential equation:

where a and b are two unknown parameters.

Our goal is to find the set of couples (a,b) that makes the origin y=0 stable. It is depicted in the figure:

*Hint: apply the Routh-Hurwitz criterion to the caracteristic polynomial of the system.*

## Unstructured Mapping¶

The complete code can be found here: `examples/lab/lab8.cpp`

.

**Introduction**

A robot is moving in a rectangular area [-L,L]x[-L,L] (with L=2) and tries to build a map while avoiding obstacles. The map is precisely the shape of obstacles inside the area.

The only information we have is when its euclidian distance to an obstacle
get smaller than 0.9. It receives an alert, a vector of measurements which contain
its own position (named `x_rob`

and `y_rob`

in the code)
and the position of the detected obstacle point
(named `x_obs`

and `y_obs`

in the code). We also know that
all the points that are less distant than 0.1 from the detected point belong to the obstacle.

The robot has a series of n=10 measurements. The goal is to build an approximation of the map using set intervals.

You can first copy-paste the data:

```
// positions of the robot
double x_rob[n]={2, 1.61803, 0.618034, -0.618034, -1.61803, -2, -1.61803, -0.618034, 0.618034, 1.61803};
double y_rob[n]={0, 1.17557, 1.90211, 1.90211, 1.17557, 0, -1.17557, -1.90211, -1.90211, -1.17557};
// positions of the obstacle point
double x_obs[n]={0.9, 0.728115, 0.278115, -0.278115, -0.728115, -0.9, -0.728115, -0.278115, 0.278115, 0.728115};
double y_obs[n]={0, 0.529007, 0.855951, 0.855951, 0.529007, 0, -0.529007, -0.855951, -0.855951, -0.529007};
```