The purpose of this chapter is to show how to create and manipulate
objects corresponding to the mathematical concept of *function*.

Let us rule out a potential ambiguity.

Since we are in the C++ programming language, the term *variable* and
*function* already refers to something precise. For instance, the following
piece of code introduces a *function* sum and a *variable* x:

```
int sum(int x, int y) {
return x+y;
}
int x=2;
```

The variable x may represent, say, the balance of a bank account.
The account number is what we call the *semantic* of x, that is, what x is supposed
to represent in the user’s mind. So, on one side, we have *what we write*, that is, a program with variables and functions, and on the other side, *what we represent*, that is, concepts
like a bank account.

With IBEX, we write programs to represent mathematical concepts
that are also called *variables* and *functions*.
The mapping \((x,y)\mapsto\sin(x+y)\) is an example of function that
we want to represent. It shall not be confused with the function sum
above.

To avoid ambiguity, we shall talk about *mathematical*
variables (resp. functions) versus *program* variables (resp. functions).
We will also use italic symbol like *x* to denote a mathematical variable
and postscript symbols like x for program variables.
In most of our discussions, variables and functions will refer
to the mathematical objects so that the mathematical meaning will be the implicit one.

Mathematical functions are represented by objects of the class `Function`.

A (mathematical) variable does not necessarily represent a single real value. It can also be a vector, a matrix or an array-of-matrices. One can, e.g., build the following function

\[\begin{split}\begin{array}{cccc}
f: & \mathbb{R}^2\times\mathbb{R}^3 & \to & \mathbb{R}\\
& (x,y) & \mapsto & x_1\times y_1+x_2\times y_2 - x_3
\end{array}.\end{split}\]

In this case, *x* and *y* are vector variables with 2 and 3 components respectively.

We see, at this point, that the term *variable* becomes ambiguous.
For instance, if I say that the function *f* takes 2 variables, we don’t really know if
it means that the function takes two arguments (that might be vectors or matrices) or if
the total input size is a vector of \(\mathbb{R}^2\).

For this reason, from now on, we will call **argument** the formal parameters
or input symbols the function has been defined with and **variable** a component of the
latters.

Hence, the function *f* in the previous paragraphs has two arguments, *x* and *y*
and 5 variables \(x_1, x_2, y_1, y_2\) and \(y_3\).

Note that, as a consequence, variables are always real-valued.

Before telling you which class represents the arguments of a function, let us say first that
this class does not play a big role.
Indeed, the only purpose of declaring an argument *x* in IBEX is
for building a function right after, like \(x\mapsto x+1\).
Functions play, in contrast, a big role.

In other words, *x* is nothing but a syntaxic leaf in the expression
of the function. In particular, an argument is not a slot for
representing domain.
E.g, if you want to calculate the range of *f* for \(x\in[0,1]\),
you just call a (program) function `eval` with a plain box in argument.
It’s just as if *f* was the function that takes one argument and
increment it, whatever the name of this argument is.

Once *f* has been built, we can almost say that *x* is no longer useful.
Arguments must be seen only as temporary objects, in the process of function construction.

Before going on, let us slightly moderate this point.
We have assumed here that, as a user of IBEX the operations you are interested in are: *evaluate* *f* on a box,
calculate *f’* on a box, solve *f(x)=0* and so on. All these operations can be qualified as numerical: they take
intervals and return intervals. You don’t need to deal again with the expression of the function, once built.
But if you need to handle, for any reason, the symbolic form of the function then you have to inspect the syntax
and arguments appear again.

We have said in the previous paragraph that an argument *x* can actually represent
n variables \(x_1,\ldots,x_n\). So each argument has some associated information about
its dimension(s).

Let us consider again this function:

\[\begin{split}\begin{array}{cccc}
f: & \mathbb{R}^2\times\mathbb{R}^3 & \to & \mathbb{R}\\
& (x,y) & \mapsto & x_1\times y_1+x_2\times y_2 - x_3
\end{array}.\end{split}\]

From the user standpoint, the function *f* (once built) is “flattened” (or “serialized”) to a mapping from \(\mathbb{R}^5\) to \(\mathbb{R}\).
Each C++ function (eval, etc.) expects a 5-dimensional box as parameter.

The way intervals are mapped to the variables components follows a straightforward ordering:
everytime we call a (program) function of *f* with the box \([b]=([b]_1,\ldots,[b]_5)\) in argument, we simply enforce

\[x\in[b]_1\times[b]_2 \quad \mbox{and} \quad y\in[b]_3\times[b]_4\times[b]_5.\]

If you don’t want to create functions in C++, you can move now to *function operations*.

As we have just said, arguments are just symbols in expression. For this reason,
they are represented by a class named `ExprSymbol`.
In fact, there is also another class we introduced for convenience, called `Variable`.
It is, of course, a very confusing name from the programer’s viewpoint since a `Variable`
does actually not represent *a variable but an argument*. However, from the user’s viewpoint,
this distinction is not visible and “variable” is more meaningful than “argument”.
Anyway, the programer never has to deal with a “Variable” object. Without going further into details, the `Variable`
class must be seen as a kind of “macro” that generates `ExprSymbol` objects.
This macro is only useful if you *build arguments in C++*.

Once built, an argument is always typed `ExprSymbol`.

If `x` is an `ExprSymbol` object, you can obtain the information about its dimensions via `x`.dim.
The `dim` field is of type `Dim`, a class that simply contains 3 integers (one for each dimension, see
the API for further details).

Finally, an argument also has a name, that is only useful for displaying. It is a
regular C string (`char*`) stored in the field `name`.

Various interval computations can be performed with a function. We detail below the main ones.

Take a look first at the *tutorial* for introductory examples.

Since function overloading does not work for return types in C++, you have to either
call `eval`, `eval_vector` or `eval_matrix` depending if your function
respectively returns a scalar, a vector or a matrix.

All `eval_XXX` functions expects a single box in argument that represents all the arguments (scalars, vectors, matrices) stored in a single flat array (see *Dimensions and ordering*).

To build this vector, the best is to use *backward projection functions*.

Here is an example with f(A,B,C)=A+B-C where A, B and C are matrices.

```
const int nb_rows=2;
const int nb_cols=2;
Variable a(nb_rows,nb_cols),b(nb_rows,nb_cols),c(nb_rows,nb_cols);
Function pA(a,b,c,a);
Function pB(a,b,c,b);
Function pC(a,b,c,c);
Function f(a,b,c, (a+b-c));
double _A[nb_rows*nb_cols][2]={{2,2},{-1,-1},{-1,-1},{2,2}};
IntervalMatrix MA(2,2,_A);
double _B[nb_rows*nb_cols][2]={{1,1},{-1,-1},{1,1},{1,1}};
IntervalMatrix MB(2,2,_B);
double _C[nb_rows*nb_cols][2]={{1,1},{0,0},{0,0},{1,1}};
IntervalMatrix MC(2,2,_C);
IntervalVector box(3*nb_rows*nb_cols);
// the backward call on pA will force the sub-vector of
// "box" that represents the domain of "a" to contain the
// interval matrix "MA".
pA.backward(MA,box);
// idem
pB.backward(MB,box);
pC.backward(MC,box);
IntervalMatrix M = f.eval_matrix(box);
output << "A+B-C=" << M << endl;
```

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].
The range [y] can be any constant (real value, interval, inteval vector, etc.).
Rigorously, given two intervals [x] and [y], the contraction gives a new interval [z] such that

\[\forall x\in[x], \quad f(x)\in[y] \Longrightarrow x \in[z] \subseteq [x].\]

One way to do this is by using the famous *forward-backward* (alias `HC4Revise`).
It is quick since it runs in linear time w.r.t. the size
of the constraint syntax and optimal when arguments have all one occurrence
in this syntax.
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);
```

One can indeed check that the resulting box is a consistent narrowing of the initial one.

Consider \(f:(x,y)\mapsto x\times y\). The first and most simple way of calculating the gradient is:

```
double init_xy[][2] = { {1,2}, {3,4} };
IntervalVector box(2,init_xy);
cout << "gradient=" << f.gradient(box) << endl;
```

Since \(\frac{\partial{f}}{\partial{x}}=y\) and \(\frac{\partial{f}}{\partial{y}}=x\) we get:

```
gradient=([3,4] ; [1,2])
```

In this first variant, the returned vector is a new object created each time the function is called.
When we have to compute many times different values of the gradient for the same function, we can
also build a vector once for all and ask the `gradient` to store the result in this slot:

```
IntervalVector g(4);
f.gradient(box,g);
cout << "gradient=" << g << endl;
```

The interval Jacobian matrix of a function *f* on a box *[x]* is

\[\begin{split}J=\left(\begin{array}{ccc}
\frac{\partial{f_1}}{\partial{x_1}}([x]) & \ldots & \frac{\partial{f_1}}{\partial{x_n}}([x])\\
\vdots \\
\frac{\partial{f_m}}{\partial{x_1}}([x]) & \ldots & \frac{\partial{f_m}}{\partial{x_n}}([x])\\
\end{array}\right)\end{split}\]

The interval Jacobian matrix is obtained exactly as for the gradient. Just write:

```
f.jacobian(box)
```

to get an `IntervalMatrix` containing an enclosure of the Jacobian matrix of *f* on the box in argument.

There is also a variant where the matrix is passed as parameter (as for the gradient) in order to avoid allocating memory for the calculated matrix:

```
f.jacobian(box,J)
```

You can also compute with IBEX the “Hansen matrix”. This matrix
is another *slope* matrix, thiner than the interval Jacobian (but slower to be calculated).
It is, for example, used inside the interval Newton operator.
The Hansen matrix corresponds to the following matrix, where \((x_1,\ldots,x_n)\) denotes
the midvector of \([x]\).

\[\begin{split}\left(\begin{array}{cccc}
\frac{\partial{f_1}}{\partial{x_1}}([x]_1,x_2,\ldots,x_n) &
\frac{\partial{f_1}}{\partial{x_2}}([x]_1,[x]_2,\ldots,x_n) &
\ldots &
\frac{\partial{f_1}}{\partial{x_n}}([x]_1,[x]_2,\ldots,[x]_n) \\
\vdots & \\
\frac{\partial{f_m}}{\partial{x_1}}([x]_1,x_2,\ldots,x_n) &
\frac{\partial{f_n}}{\partial{x_2}}([x]_1,[x]_2,\ldots,x_n) &
\ldots &
\frac{\partial{f_m}}{\partial{x_n}}([x]_1,[x]_2,\ldots,[x]_n) \\
\end{array}\right)\end{split}\]

Here is an example:

```
Variable x,y;
Function f(x,y,Return(sqr(x)*y,sqr(y)*x));
IntervalMatrix H(2,2);
IntervalVector box(2,Interval(1,2));
f.hansen_matrix(box,H);
output << "Hansen matrix:\n" << H << endl;
```

The display is:

```
Hansen matrix:
(([3, 6] ; [1, 4])
([2.25, 2.25] ; [2, 8]))
```

`Function` objects are very easy to build.

This section explains how to build them
using C++ operator overloading but using
the *Minibex* syntax is even simpler.

The following piece of code creates an argument `x` and prints it:

```
Variable x;
cout << x << endl;
```

The first instruction creates a (program) variable `x`. It is initialized by default, since
nothing is given here to the constructor.
By default, the argument is real (or *scalar)*, meaning it is not a vector nor a matrix.
Furthermore, the argument has a name that is automatically
generated. Of course, the name of the argument does not necessarily correspond to the name of the
program variable.
For instance, `x` is the name of a C++ variable but the corresponding argument is named *_x_0*.
The second instruction prints the name of the argument on the standard output:

```
_x_0
```

It is possible to rename arguments, see below.

To create a n-dimensional vector argument, just give the number n as a parameter to the constructor:

```
Variable y(3); // creates a 3-dimensional vector
```

To create a mxn matrix, give m (number of rows) and n (number of columns) as parameters:

```
Variable z(2,3); // creates a 2*3-dimensional matrix
```

We can go like this up to 3 dimensional arrays:

```
Variable t(2,3,4); // creates a 2*3*4-dimensional array
```

Usually, you don’t really care about the names of arguments since you handle program variables in your code. However, if you want a more user-friendly display, you can specify the name of the argument as a last parameter to the constructor.

In the following example, we create a scalar, a vector and a matrix argument each time with a chosen name.:

```
Variable x("x"); // creates a real argument named "x"
Variable y(3,"y"); // creates a vector argument named "y"
Variable z(2,3,"z"); // creates a matrix argument named "z"
cout << x << " " << y << " " << z << endl;
```

Now, the display is:

```
x y z
```

The following piece of code creates the function \((x,y)\mapsto \sin(x+y)\):

```
Variable x("x");
Variable y("y");
Function f(x,y,sin(x+y)); // create the function (x,y)->sin(x+y)
output << f << endl;
```

The display is:

```
_f_7:(x,y)->sin((x+y))
```

You can directly give up to 20 variables in argument of the `Function` constructor:

```
Variable a,b,c,d,e,f;
Function _f(a,b,c,d,e,f,a+b+c+d+e+f);
```

If more than 20 variables are needed, you need to build an intermediate array for collecting the arguments.
More precisely, this intermediate object is an `Array<const ExprSymbol>`. The usage is summarized below. In this
example, we have 7 variables. But instead of creating the function

\[x\mapsto x_1+\ldots+x_7\]

with one argument (a vector with 7 components), we decide to create the function

\[(x_1,\ldots,x_7)\mapsto x_1+\ldots+x_7.\]

with 7 arguments (7 scalar variables):

```
Variable x[7]; // not to be confused with x(7)
Array<const ExprSymbol> vars(7);
for (int i=0; i<7; i++)
vars.set_ref(i,x[i]);
Function f(vars, x[0]+x[1]+x[2]+x[3]+x[4]+x[5]+x[6]);
output << f << endl;
```

The display is:

```
_f_9:(_x_11,_x_12,_x_13,_x_14,_x_15,_x_16,_x_17)->((((((_x_11+_x_12)+_x_13)+_x_14)+_x_15)+_x_16)+_x_17)
```

**Note:** Because of a potential conflict with `std::min` (or `std::max`), you might be forced to prefix the min (max) function with `ibex::`:

```
Variable x,y;
Function f(x,y,ibex::min(x,y));
output << f << endl;
```

If arguments are vectors, you can refer to the component of an argument using square brackets. Indices start by 0, following the convention of the C language.

We rewrite here the previous distance function using 2-dimensional
arguments `a` and `b` instead:

```
Variable a(2);
Variable b(2);
Function dist(a,b,sqrt(sqr(a[0]-b[0])+sqr(a[1]-b[1])),"dist");
```

To define a vector-valued function, the `Return` keword allows
you to list the function’s components.

See the example in the *tutorial*.

C++ operator overloading allows you to create a DAG instead of an expression tree. This will result in a gain in performance.
For that, you need to handle references of shared subexpressions with variables types `const ExprNode&`.

In the following example we create the function :

\[f:x\mapsto ((\cos(x)+1)^2, (\cos(x)+1)^3)\]

and we want the subexpression cos(x)+1 to be shared:

```
Variable x;
const ExprNode& e=cos(x)+1;
Function f(x,Return(pow(e,2),pow(e,3)));
```

Let us build a function that returns the sum of the square of `N` variables, where `N` is some constant.

The only difficulty is that we cannot assign references in C++, so we need to use pointers to `(const ExprNode&`) instead:

```
int N=10;
Variable x(N,"x");
const ExprNode* e=&(sqr(x[0]));
for (int i=1; i<N; i++)
e = & (*e + sqr(x[i]));
Function f(x,*e,"f");
output << f << endl;
```

The display is:

```
f:(x)->(((((((((x[0]^2+x[1]^2)+x[2]^2)+x[3]^2)+x[4]^2)+x[5]^2)+x[6]^2)+x[7]^2)+x[8]^2)+x[9]^2)
```

By default, function names are also generated. But you can also set your own function name, as the last parameter of the constructor:

```
Function f(x,y,sin(x+y),"f");
```

The following symbols are allowed in expressions:

sign, min, max, sqr, sqrt, exp, log, pow, cos, sin, tan, acos, asin, atan, cosh, sinh, tanh, acosh, asinh, atanh atan2

Power symbols `^` are not allowed. You must
either use `pow(x,y)`, or simply `sqr(x)` for the square function.

Here is an example of the distance function between `(xa,ya)` and `(xb,yb)`:

```
Variable xa,xb,ya,yb;
Function dist(xa,xb,ya,yb, sqrt(sqr(xa-xb)+sqr(ya-yb)));
```

(to do)

You can compose functions. Each argument of the called function can be substitued by an argument of the calling function, a subexpression or a constant value.

See the example in the *tutorial*.

We have explained how to create a function with *an arbitrary number of arguments*.
We explain now how to call (perform composition) with such function.

It is as simple as storing all the actual arguments in an array structure, namely, a structure
of expression nodes (typed `Array<const ExprNode>`).

However, when an actual argument is not a formal expression but a numerical constant (data), it is necessary to
explicitly encapsulate this constant in a expression node. This is what the `ExprConstant` class stands for.

Here is an example. We create the function \(f:(x,y)\mapsto x+y\) and apply it to the hybrid couple (z,1) where z is another variable. We do it in the generic way, using arrays:

```
Variable x,y;
// formal arguments
Array<const ExprSymbol> vars(2);
vars.set_ref(0,x);
vars.set_ref(1,y);
Function f(vars, x+y,"f");
// actual arguments
const ExprSymbol& z=ExprSymbol::new_("z"); // <=> "Variable z" (but more "safe")
const ExprConstant& c=ExprConstant::new_scalar(1.0);
// =============================================
// before release 2.1.6:
// const ExprNode* args[2];
// args[0]=&z;
// args[1]=&c;
// before release 2.1.6:
// from release 2.1.6 and subsequents:
Array<const ExprNode> args(2);
args.set_ref(0,z);
args.set_ref(1,c);
// =============================================
Function g(z,f(args),"g");
output << g << endl;
```

The display is:

```
g:(z)->f(z,1)
```

Differentiation of a function is another function.
So symbolic differentiation is obtained via a copy constructor where the copy “mode” is set to the special value `Function::DIFF`:

```
Function f("x","y","z","x*y*z");
Function df(f,Function::DIFF);
output << "df=" << df << endl;
```

The output is

```
df=d_f_13:(x,y,z)->((z*y),(z*x),(x*y))
```