@c This is part of the GNU Octave Interval Package Manual.
@c Copyright 2015 Oliver Heimlich.
@c Copyright 2017 Joel Dahne
@c See the file manual.texinfo for copying conditions.
@documentencoding UTF-8
@include macros.texinfo

@chapter Examples

This chapter presents some more or less exotic use cases for the interval package.

@section Arithmetic with System-independent Accuracy

According to IEEE Std 754 only the most basic floating-point operations must be provided with high accuracy.  This is also true for the arithmetic functions in Octave.  It is no surprise that many arithmetic functions fail to provide perfect results and their output may be system dependent.

We compute the cosecant for 100 different values.

@example
@group
x = vec (1 ./ magic (10));
sum (subset (csc (x), csc (infsupdec (x))))
  @result{} ans =  98
@end group
@end example

Due to the general containment rule of interval arithmetic @code{x ∈ X ⇒ f (x) ∈ f (X)} one would expect the @code{csc (x)} to always be contained in the interval version of the cosecant for the same input.  However, the classic cosecant is not very accurate whereas the interval version is.  In 2 out 100 cases the built-in cosecant is less accurate than 1 ULP.

@section Prove the Existence of a Fixed Point

A weaker formulation of Brower's fixed-point theorem goes: If @var{x} is a bounded interval and function @var{f} is continuous and @var{f} (@var{x}) ⊂ @var{x}, then there exists a point @var{x}₀ ∈ @var{x} such that @var{f} (@var{x}₀) = @var{x}₀.

These properties can be tested automatically.  Decorated intervals can even prove that the function is continuous.

@example
@group
x = infsupdec ("[-1, +1]");
f = @@cos;
subset (f (x), x)
  @result{} ans =  1
iscommoninterval (x)
  @result{} ans =  1
continuous = strcmp (decorationpart (f (x)), "com")
  @result{} continuous =  1
@end group
@end example

Furthermore it is sometimes possible to approximate the fixed-point by repetitive evaluation of the function, although there are better methods to do so in general.

@example
@group
for i = 1 : 20
    x = f (x);
endfor
display (x)
  @result{} x ⊂ [0.73893, 0.73919]_com
@end group
@end example

@menu
Further Examples
* Floating-point Numbers:: Analyze properties of binary64 numbers with intervals
* Root Finding:: Find guaranteed enclosures for roots of a function
* Parameter Estimation:: Examples of set inversion via interval analysis
* Path Planning:: Find a feasible path between two points
@end menu


@node Floating-point Numbers
@section Floating-point Numbers

Floating-point numbers are most commonly used in binary64 format, a.k.a. double precision.  Internally they are stored in the form @code{± @var{m} * 2 ^ @var{e}} with some integral mantissa @var{m} and exponent @var{e}.  Most decimal fractions can only be stored approximately in this format.

The @funref{intervaltotext} function can be used to output the approximate value up to the last decimal digit.

@example
@group
intervaltotext (infsup (0.1), "[.55g]")
  @result{} ans = [0.1000000000000000055511151231257827021181583404541015625]
@end group
@end example

It can be seen that 0.1 is converted into the most accurate floating-point number.  In this case that value is greater than 0.1.  The next lower value can be seen after producing an interval enclosure around 0.1 with the nearest floating-point numbers in each direction.

@example
@group
intervaltotext (infsup ("0.1"), "[.55g]")
  @result{} ans =
    [0.09999999999999999167332731531132594682276248931884765625,
     0.1000000000000000055511151231257827021181583404541015625]
@end group
@end example

The error of this approximation can be examined with the @funref{@@infsup/wid} function.

@example
@group
wid (infsup ("0.1"))
  @result{} ans =    1.3878e-17
@end group
@end example

With the @funref{@@infsup/nextout} function an interval can be enlarged in each direction up to the next floating-point number.  Around zero the distance towards the next floating point number is very small, but gets bigger for numbers of higher magnitude.

@example
@group
wid (nextout (infsup ([0, 1, 1e10, 1e100])))
  @result{} ans =

       9.8813e-324    3.3307e-16    3.8147e-06    3.8853e+84

@end group
@end example


@node Root Finding
@section Root Finding
@subsection Interval Newton Method

In numerical analysis, @uref{https://en.wikipedia.org/wiki/Newton%27s_method,Newton's method} can find an approximation to a root of a function.  Starting at a location @var{x}₀ the algorithms executes the following step to produce a better approximation:

@display
@var{x}₁ = @var{x}₀ @minus{} @var{f} (@var{x}₀) / @var{f}' (@var{x}₀)
@end display

The step can be interpreted geometrically as an intersection of the graph's tangent with the x-axis.  Eventually, this may converge to a single root. In interval analysis, we start with an interval @var{x}₀ and utilize the following interval Newton step:

@display
@var{x}₁ = (mid (@var{x}₀) @minus{} @var{f} (mid (@var{x}₀)) / @var{f}' (@var{x}₀)) ∩ @var{x}₀
@end display

Here we use the pivot element @code{mid (@var{x}₀)} and produce an enclosure of all possible tangents with the x-axis.  In special cases the division with @code{@var{f}' (@var{x}₀)} yields two intervals and the algorithm bisects the search range.  Eventually this algorithm produces enclosures for all possible roots of the function @var{f} in the interval @var{x}₀.  The interval newton method is implemented by the function @funref{@@infsup/fzero}.

To produce the derivative of function @var{f}, the automatic differentiation from the symbolic package bears a helping hand.  However, be careful since this may introduce numeric errors with coefficients.

@example
@group
@c doctest: +SKIP_IF(strcmp ('Not installed', nthargout(2, 'pkg', 'describe', 'symbolic')@{1@}) || compare_versions (ver ("symbolic").Version, "2.5.0", "<"))
f = @@(x) sqrt (x) + (x + 1) .* cos (x);
@end group
@group
pkg load symbolic
df = function_handle (diff (formula (f (sym ("x")))))
  @result{} df =

   @@(x) -(x + 1) .* sin (x) + cos (x) + 1 ./ (2 * sqrt (x))
@end group
@group
@c doctest: +SKIP_IF(strcmp ('Not installed', nthargout(2, 'pkg', 'describe', 'symbolic')@{1@}) || compare_versions (ver ("symbolic").Version, "2.5.0", "<"))
fzero (f, infsup ("[0, 6]"), df)
  @result{} ans ⊂ 2×1 interval vector

        [2.059, 2.0591]
       [4.3107, 4.3108]

@end group
@end example

We could find two roots in the interval [0, 6].


@subsection Bisection

Consider the function @code{f (@var{x}, @var{y}) = -(5*@var{y} - 20*@var{y}^3 + 16*@var{y}^5)^6 + (-(5*@var{x} - 20*@var{x}^3 + 16*@var{x}^5)^3 + 5*@var{y} - 20*@var{y}^3 + 16*@var{y}^5)^2}, which has several roots in the area @var{x}, @var{y} ∈ [-1, 1].

@myimage{image/poly-example-surf.m,Surface plot of @code{f (@var{x}, @var{y})} which shows a lot of roots for the function}

The function is particular difficult to compute with intervals, because its variables appear several times in the expression, which benefits overestimation from the dependency problem.  Computing root enclosures with the @funref{@@infsup/fsolve} function is unfeasible, because many bisections would be necessary until the algorithm terminates with a useful result.  It is possible to reduce the overestimation with the @funref{@@infsup/polyval} function to some degree, but since this function is quite costly to compute, it does not speed up the bisecting algorithm.

@include image/poly-example-roots-simple.m.texinfo
@myimage{image/poly-example-roots-simple.m,Enclosures of roots for the function @code{f (@var{x}, @var{y})}}

Now we use the same algorithm with the same number of iterations, but also utilize the @emph{mean value theorem} to produce better enclosures of the function value with first order approximation of the function. The function is evaluated at the interval's midpoint and a range evaluation of the derivative can be used to produce an enclosure of possible function values.

@include image/poly-example-roots-with-deriv.m.texinfo

By using the derivative, it is possible to reduce overestimation errors and achieve a much better convergence behavior.

@myimage{image/poly-example-roots-with-deriv.m,Enclosures of roots for the function @code{f (@var{x}, @var{y})}}


@node Parameter Estimation
@section Parameter Estimation
@subsection Small Search Space

Consider the model @code{y (@var{t}) = @var{p1} * exp (@var{p2} * t)}.  The parameters @var{p1} and @var{p2} are unknown, but it is known that the model fulfills the following constraints, which have been obtained using measurements with known error bounds.

@display
@verbatim
p1, p2  ∈   [-3, 3]
y (0.2) ∈  [1.5, 2]
y (1)   ∈  [0.7, 0.8]
y (2)   ∈  [0.1, 0.3]
y (4)   ∈ [-0.1, 0.03]
@end verbatim
@end display

A better enclosure of the parameters @var{p1} and @var{p2} can be estimated with the @funref{@@infsup/fsolve} function.

@example
@group
## Model
y = @@(p1, p2, t) p1 .* exp (p2 .* t);
## Observations / Constraints
t = [0.2; 1; 2; 4];
y_t = infsup ("[1.5, 2]; [0.7, 0.8]; [0.1, 0.3]; [-0.1, 0.03]");
## Estimate parameters
f = @@(p1, p2) y (p1, p2, t);
p = fsolve (f, infsup ("[-3, 3]; [-3, 3]"), y_t)
  @result{} p ⊂ 2×1 interval vector

         [1.9863, 2.6075]
       [-1.3243, -1.0429]

@end group
@end example

The resulting @code{p} guarantees to contain all parameters @code{[@var{p1}; @var{p2}]} which satisfy all constraints on @var{y}. It is no surprise that @code{f (p)} intersects the constraints for @var{y}.

@example
@group
f (p(1), p(2))
  @result{} ans ⊂ 4×1 interval vector

            [1.5241, 2.1166]
          [0.52838, 0.91888]
          [0.14055, 0.32382]
       [0.0099459, 0.040216]

@end group
@end example


@subsection Larger Search Space

Consider the function @code{f (x) = @var{p1} ^ x * (@var{p2} + @var{p3} * x + @var{p4} * x^2)}.  Let's say we have some known function values (measurements) and want to find matching parameters @var{p1} through @var{p4}.  The data sets (@var{x}, @var{y}) can be simulated.  The parameters shall be reconstructed from the observed values on the search range @var{p}.

Using plain @funref{@@infsup/fsolve} would take considerably longer, because the search range has 4 dimensions.  Bisecting intervals requires an exponential number of steps and can easily become inefficient.  Thus we use a contractor for function @var{f}, which in addition to the function value can produce a refinement for its parameter constraints.  Contractors can easily be build using interval reverse operations like @funref{@@infsup/mulrev}, @funref{@@infsup/sqrrev}, @funref{@@infsup/powrev1}, etc.

@example
@c doctest: -TEXINFO_SKIP_BLOCKS_WO_OUTPUT
@group
## Simulate some data sets and add uncertainty
x = -6 : 3 : 18;
f = @@(p1, p2, p3, p4) ...
    p1 .^ x .* (p2 + p3 .* x + p4 .* x .^ 2);
y = f (1.5, 1, -3, 0.5) .* infsup ("[0.999, 1.001]");
@end group
@group
function [fval, p1, p2, p3, p4] = ...
    contractor (y, p1, p2, p3, p4)
    x = -6 : 3 : 18;
    ## Forward evaluation
    a = p1 .^ x;
    b = p3 .* x;
    c = p2 + b;
    d = p4 .* x .^ 2;
    e = c + d;
    fval = a .* e;
    ## Reverse evaluation and
    ## undo broadcasting of x
    y = intersect (y, fval);
    a = mulrev (e, y, a);
    e = mulrev (a, y, e);
    p1 = powrev1 (x, a, p1);
    p1 = intersect (p1, [], 2);
    c = intersect (c, e - d);
    d = intersect (d, e - c);
    p2 = intersect (p2, c - b);
    p2 = intersect (p2, [], 2);
    b = intersect (b, c - p2);
    p3 = mulrev (x, b, p3);
    p3 = intersect (p3, [], 2);
    p4 = mulrev (x .^ 2, d, p4);
    p4 = intersect (p4, [], 2);
endfunction
@end group
@end example

Now, search for solutions in the range of @code{p} and try to restore the function parameters.

@example
@group
p = infsup ("[1.1, 2] [1, 5] [-5, -1] [0.1, 5]");
p = fsolve (@@contractor, ...
            p, y, ...
            struct ("Contract", true))'
  @result{} p ⊂ 4×1 interval vector

         [1.4991, 1.5009]
              [1, 1.0011]
       [-3.0117, -2.9915]
       [0.49772, 0.50578]

@end group
@end example

The function parameters 1.5, 1, @minus{}3, and 0.5 from above could be restored.  The contractor function could significantly improve the convergence speed of the algorithm.


@subsection Combination of Functions

Sometimes it is hard to express the search range in terms of a single function and its constraints, when the preimage of the function consists of a union or intersection of different parts.  Several contractor functions can be combined using @funref{ctc_union} or @funref{ctc_intersect} to make a contractor function for more complicated sets.  The combined contractor function allows one to solve for more complicated sets in a single step.

@include image/contractor-rings-union.m.texinfo
@myimage{image/contractor-rings-union.m,Set inversion for two rings}

Intersections of contractor functions are especially useful to apply several constraints at once.  For example, when it is known that a particular location has a distance of @var{a} ∈ [3, 4] from object A, located at coordinates (1, 3), and a distance of @var{b} ∈ [5, 6] from object B, located at coordinates (2, -1), the intersection of both rings yields all possible locations in the search range.  The combined contractor function enables fast convergence of the search algorithm.

@include image/contractor-rings-intersect.m.texinfo
@myimage{image/contractor-rings-intersect.m,Set inversion for intersection of two rings}


@node Path Planning
@section Path Planning

@float Figure,cameleon-problem
@caption{Cameleon Problem: The polygon has to be moved from the left to the right without touching any obstacles along the path.}
@shortcaption{Cameleon Problem Description}
@myimage{image/cameleon-start-end.svg,Cameleon Problem: Start and End Position}
@end float

The problem presented here is a simplified version from the paper L. Jaulin (2001). @uref{https://www.ensta-bretagne.fr/jaulin/cameleon.html,Path planning using intervals and graphs.} Reliable Computing, issue 1, volume 7, 1–15.

There is an object, a simple polygon in this case, which shall be moved from a starting position to a specified target position.  Along the way there are obstacles which may not be touched by the polygon.  The polygon can be moved in one direction (left to right or right to left) and may be rotated around its lower left corner.

This makes a two dimensional parameter space and any feasible positions can be determined using interval arithmetic like in the examples above.

Then we use a simple path planning algorithm: We move along the centers of adjacent and feasible boxes in the parameter space until we have a closed path from the start position to the end position. The path is guaranteed to be feasible, that is, there will be no collisions if we follow the path.

@include image/cameleon.m.texinfo

The script visualizes the solution in the parameter space.  Unfeasible parameters are white, and uncertain combinations of parameters are red.  The algorithm's accuracy is just good enough to find a closed path, which is drawn in green color.  The uncertain red area is quite big because we have used a very simple check for verification whether the polygon overlaps the obstacles.  This could be improved.

@myimage{image/cameleon.m,Computed feasible path in parameter space}

The solution is not optimal, please refer to Luc Jaulin's paper for more sophisticated approaches.  However, we could find a valid solution that moves the polygon as desired without touching any obstacles.

@float Figure,cameleon-solution
@caption{Cameleon Problem: A possible solution which moves the polygon from the left to the right without touching obstacles.}
@shortcaption{Cameleon Problem Solution}
@html
<object data="image/cameleon-animation.svg" type="image/svg+xml">
<param name="src" value="image/cameleon-animation.svg" />
@end html
@myimage{image/cameleon-transition.svg,Cameleon Problem: Transition from Start to End Position}
@html
</object>
@end html
@end float