Skip to content

RealInterval Example Algos

Serge Stinckwich edited this page Jun 26, 2018 · 1 revision

There are a few algorithms implemented in this project, that I'll show you here. But first let's have a short look how these algos are generally working.

FunctionRefinement

suppose you have the function f(x) = x - x² defined over a domain x ∈ [0,1] and you want to know the range of this function.

f:=[:x|x - x squared ].
domain:=0 hull:1.0.

For the theoretically inclined reader, an interval extension F of a function f is simply an interval function, where if the number x is in an interval X then f(x) is always in the interval F(X).

Of course we can use f itself as an interval extension of the original function f(x) = x - x² and apply domain.

f value: domain."--> [-1.0,1.0]"

The result certainly encompasses the range, but it is not a given that every number in the result actually is in the range.

Since x - x² = x (1 - x ), we can also use x (1 - x ) as an interval extension of the original function f(x):

domain * (1 - domain)."--> [0.0,1.0]"

Now this result is tighter, hence better. And one can see that functions that are equal, when defined on real numbers, are unequal when defined on intervals. Let's try another possibility, x - x² = ¼ - ( ½ - x )²:

(1/4)-((1/2)-domain)squared."--> [0.0,(1/4)]"

This result is now rather compact and actually the real range of f(x). One sees that it is better to rearrange interval equations in such a way that each variable appears just once. Of course this not always possible or can be complicated to do by hand. What one wants is an algo that can produce a similar result with the original block. Let's split the domain in two parts, and calc the union of the results:

(f value: (0.5 hull:1))union: (f value: (0 hull:0.5))."--> [-0.5,0.75]"

This result, although not as good as the last one, is nevertheless better than the first one. Of course if one would split the domain in more parts, the result would be better. The simple implementation FunctionRefinement0 does exactly that, splitting the domain in thousand intervals of equal width and calculating the result. Have a look at the naive implementation.

ff:=FunctionRefinement0 function: f interval: domain .
ff evaluate ."--> [-0.0010000000000000009,0.251]" 

This result in a way nicely encompasses the best solution [0.0,(1/4)]. Practically all interval algos are based on this trick. They usually split an interval only in two subintervals; if the result of a subinterval is good enough they keep it, if its not good enough then they split it again, and so on. And in many algos the result of a subinterval is bad enough or empty so that this interval can be thrown away. That's the whole trick. In FunctionRefinement1 a much faster and better algo is used:

ff:=FunctionRefinement1 function: f box: domain .
ff evaluate."--> [-1.1102230246251565e-16,0.2500211522783466]"

There is also a difference in the interface in so far as FunctionRefinement1 not only accepts RealIntervals but also IntervalBoxes, which are a special kind of Arrays of RealIintervals or IntervalUnions, in other words a multi-dimensional vector of intervals. A normal Array of RealIntervals would not work here as a normal Array does not understand split. The use of IntervalBoxes is shown in the next part anyway.

An important Remark

While we will use the normal function f(x) defined for Numbers as interval extension in the following examples, you should be aware of the fact that this is not always possible. Especially if f(x) uses constructions similar to x<y ifTrue: block1 ifFalse: block2 you have to construct your interval extension by hand: If x or y are intervals, x<y will not produce a Boolean but an error, because the result could be partly true. This is the reason why < and similar methods are not implemented for RealIntervals. Even constructions like x=y ifTrue: do not necessarily produce the intended result, if eg x includes:y, but x~=y; in other words if x=y should be partly true (for only some numbers in x), but returns false.

SkelboeMoore

This algo calculates the minimum of a function in a given interval essentially by repeatedly splitting this interval, calculating f of a value in this split and then throwing away all splits, whose infimum is higher than this value. There are two implementations, that have different advantages. We will check the minimum of De Jongs first function, which is just the square of a vector and has a minimum if every element of the vector is zero.

     f(x) = ∑i=1,…,n xi 2

f:=[:x| x squared sum].

We'll use a 4-dimensional search range (that is n = 4), from -5 to 5 for each of the variables xi :

a:= -5 hull:5.
searchrange := IntervalBox with: a with: a with: a with: a."--> 
an IntervalBox([-5,5] [-5,5] [-5,5] [-5,5])"
f value: searchrange."-->  [0,100]"

The minimum is of course 0 for #(0 0 0 0).

minimizer := SkelboeMoore2 function: f box: searchrange.
r := minimizer evaluate."--> an Array(an IntervalBox([-0.001220703125,0.001220703125] [-0.00244140625,0.00244140625] [-0.00244140625,0.00244140625] [-0.00244140625,0.00244140625]) [0,1.9371509552001953e-5])" 

The result is an Array of a 4-dimensional interval-vector of the argument that results in a minimum and of a interval in which that minimum lies. We'll test the midpoint and an extreme point of the intervals:

r first midPoint."--> #(0.0 0.0 0.0 0.0)"
t:=r first inf."--> #(-0.001220703125 -0.00244140625 -0.00244140625 -0.00244140625)"
f value:t."--> 1.9371509552001953e-5"
(f value:t) = r second sup."--> true"

The minimum interval [0,1.9371509552001953e-5] is relatively wide since the defaultPrecision for SkelboeMoore2 is 2e-5 and we can get tighter results this way:

minimizer desiredPrecision: 1e-10.
r:=minimizer evaluate. "--> an Array(an IntervalBox([-4.76837158203125e-6,4.76837158203125e-6] [-4.76837158203125e-6,4.76837158203125e-6] [-4.76837158203125e-6,4.76837158203125e-6] [-4.76837158203125e-6,4.76837158203125e-6]) [0,9.094947017729282e-11])" 

If the desired precision is too tiny, then the maximumIterations will be to low and there will be regions that are not explored by the algo:

minimizer desiredPrecision: 1e-11.
minimizer evaluate."--> an Array(an IntervalBox([-1.1920928955078125e-6,5] [-1.1920928955078125e-6,5] [-5,1.1920928955078125e-6] [-5,7.62939453125e-5]) [0,100])"

You see the minimum lies in [0,100]. Not very helpful. You can also check whether the desired precision is reached, this way:

minimizer desiredPrecision."--> false"
"and this shows the number of not yet evaluated split-results"
minimizer inExactResults ."--> 2"
"but if you want to look at these results, you have to set this before evaluation:"
minimizer useSloppySimplification: false.
r:=minimizer evaluate.
minimizer inExactResults ."--> 2" 
"finding those inexact results is only in this case easy, they are here the first two:"
r first first."--> an IntervalBox([0,5] [0,5] [-5,0] [-5,0])"
r first second." an IntervalBox([0,3.814697265625e-5] [0,3.814697265625e-5] [-3.814697265625e-5,0] [0,7.62939453125e-5])"
"using the later explained BoxVisualizer can sometimes help understanding the problem here"
" but anyway, usually raising the maximumIterations helps"
minimizer maximumIterations: 2000.
minimizer useSloppySimplification: true. "back to default"
minimizer evaluate."--> an Array(an IntervalBox([-1.1920928955078125e-6,1.1920928955078125e-6] [-1.1920928955078125e-6,1.1920928955078125e-6] [-1.1920928955078125e-6,1.1920928955078125e-6] [-2.384185791015625e-6,2.384185791015625e-6]) [0,9.947598300641403e-12])"
minimizer desiredPrecision."--> true"

One word about the useSloppySimplification: true default: Without this you will usually get a rather messy output. The first part of the resulting Array will be a usually huge collection of IntervalBoxes that result in the desired minimum and the sloppySimplification lays a kind of hull, occasionally something a little smaller, around these results. Of course this is somehow wrong, but then one can shut it off.

SkelboeMoore1 can be used in essentialy the same way. Its results are only more exact because its default precision is set higher. Those two implementations have different advantages: An important difference is that SM1 returns the most exact result it can get within the given number of maxIterations - this is the reason why its default precision is set higher than in SM2 -, while SM2 - as we have just seen - returns nonsense, if it can't get the desired precision within the given number of maxIterations (one might want to set its default maxIterations higher than in SM1). SM2 is generally faster though. A simple explanation of SM1 can be found in http://arxiv.org/abs/physics/0302034 page 4, SM2 can be found eg in http://www.theinfinitycomputer.com/jily6.ps.gz page 2.

A comparison between the two SkelboeMoore algos and AnotherGeneticOptimizer you can find in the comparison Notes

Constraints

The following algos use Constraints, hence i first explain them: They consist of a Block and the admissible results of that Block. The Block may have only one argument, which of course can be a SequenceableCollection (as a container of several arguments). The admissible results, set with admissibleImage: possibleValues must in the general case be something that understands includes: and intersection: and the result of intersection: has to understand isEmpty (iow any collection would be sufficient). In our cases admissibleImage will always be a RealInterval or IntervalUnion. A few examples, how a Constraint can be constructed:

c:=Constraint block: [ :x | x squared ]."--> a Constraint([ :x | x squared ] <= 0)"
c min: 2."--> [2,Float infinity]" 
c."--> a Constraint(2 <= [ :x | x squared ])"
"or:"
c max: 3."--> [Float infinity negated,3]"
c."--> a Constraint([ :x | x squared ] <= 3)"
"or:"
c equalToNumber:  3."--> [3,3]"
c."--> a Constraint([ :x | x squared ] = 3)"
"or:"
c admissibleImage: (0 hull:3)."--> [0,3]"
c."--> a Constraint(0 <= [ :x | x squared ] <= 3)"

With a Constraint one can test wether a value (in our case a Number, RealInterval, IntervalUnion or IntervalBox) fullfills that constraint completely, partly or not at all; the Constraint c in its last form tests   0 ≤ arg² ≤ 3 .

c isCompletelyAdmissibleValue: (0 hull:1.5)."--> true"
c isCompletelyAdmissibleValue: (0 hull:2)."--> false"
c isPerhapsAdmissibleValue: (0 hull:1.5)."--> true"
c isPerhapsAdmissibleValue: (0 hull:2)."--> true"
c isPerhapsAdmissibleValue: 5."--> false"
c isNotAdmissibleValue: 5."--> true"
"or if you want to test for all three possibilities, it is faster to use #admissibleNumberForValue: aValue,
which returns 0 if aValue is not admissible, 1 if it is partly admissible and 2 if it is completely admissible"
c admissibleNumberForValue: 5."--> 0"
c admissibleNumberForValue: (1 hull:5)."--> 1"
c admissibleNumberForValue: 1."--> 2"

Now back to the algos:

BoxSplicing & ConstraintRefinement

BoxSplicing is an algo that repeatedly tries to shave off unnecessary parts of an interval at the outer borders as long as possible. An description of this algo can eg be found in http://arxiv.org/pdf/physics/0302034v2.pdf page 7.

Lets say we want to solve this equation: x³ = 5x + 1. We can make a simple constraint out of that:

c:=Constraint block: [ :x ||x1|x1:=x first.(x1 raisedToInteger: 3) - (5 * x1) ] .
"i did the complication via x1, since the constraint related algos at the moment only accept IntervalBoxes,
even in the 1-dimensional case! i could change that but it would slow down computation, or i was just
too lazy."
c equalToNumber: 1. 
c."--> a Constraint([ :x |  | x1 | x1 := x first. (x1 raisedToInteger: 3) - (5 * x1) ] = 1)"
b:=BoxSplicing constraints:c.
b doRepeatedSplicing."-->Error:Box is not initialized'"
"well, we have first to give the algo a region, where it looks for a solution, lets try this one" 
b box: (IntervalBox  with: (RealInterval inf: 0 sup:1)) . "use an IntervalBox here!"

Now Boxsplicing will shave off the unnecessary parts and make that interval smaller:

b doRepeatedSplicing."--> an IntervalBox([empty])" 
"obviously no solution in there. one would like to use an infinite RealInterval, but unfortunately this algo cant really deal with them. 

well ok, lets try something else:"
b box: (IntervalBox  with: (RealInterval inf: 0 negated sup: 1000 )) .
b doRepeatedSplicing.
"--> an IntervalBox([2.3300587395679817,2.3300587395679826])"
" lets see whether this is the best result:"
b result first inf successor successor."--> 2.3300587395679826"
"iow here it does not get the smallest interval. but without further information a better result is not possible; and the result is ok:"
(c block value:  (IntervalBox  with: (RealInterval fromNumber: 2.3300587395679817 )))asNumber."---> 0.9999999999999964"
(c block value:  (IntervalBox  with: (RealInterval fromNumber: 2.3300587395679826 )))asNumber."---> 1.000000000000007"

An unimportant hint (more a reminder for myself) why a better result is not possible, is here.

Well, let's try another box:

b box: (IntervalBox  with: (RealInterval inf: -100  sup: 100 )) .
b doRepeatedSplicing.
"--> an IntervalBox([-2.1284190638445786,2.330058739567983])"
"there is obviously some interval splitting necessary, which this algo cant do. lets do it by hand:"
b box: (IntervalBox  with: (RealInterval inf: -10  sup: 0 )) .
b doRepeatedSplicing." an IntervalBox([-2.1284190638445777,-0.2016396757234043])"

"do the splitting once more:"
b box: (IntervalBox  with: (RealInterval inf: -10  sup: -2 )) .
b doRepeatedSplicing."---> an IntervalBox([-2.1284190638445777,-2.128419063844577])"
"here it again does not get the smallest intervall:"
b result first inf successor successor. "-2.128419063844577 " 
x := IntervalBox  with: (RealInterval fromNumber: -2.1284190638445777 ).
(c block value: x)."---> [0.9999999999999964,0.9999999999999964]"
x := IntervalBox  with: (RealInterval fromNumber: -2.128419063844577 ).
(c block value: x)asNumber." 1.0000000000000053"

"the other side of the split:"
b box: (IntervalBox  with: (RealInterval inf: -2 sup: 0 )) .
b doRepeatedSplicing."-->an IntervalBox([-0.2016396757234047,-0.20163967572340466])"
"a best possible result:"
b result first isAtomic ."--> true" "or:"
b result first inf successor =b result first sup."--> true"
x := IntervalBox  with: (RealInterval fromNumber: b result first inf ).
(c block value: x)asNumber." 1.0000000000000002"
x := IntervalBox  with: (RealInterval fromNumber: b result first sup  ).
(c block value: x)asNumber." 0.9999999999999999"
"iow the result is not an exact Float!"

On the other hand all this splitting can be automated with ConstraintRefinement:

b :=ConstraintRefinement constraints: c box: (IntervalBox  with: (RealInterval inf: -10  sup: 10 )) .
" -->a ConstraintRefinement(box: an IntervalBox([-10,10]) constraints: a ConstraintsCollection(a Constraint([ :x |  | x1 | x1 := x first. (x1 raisedToInteger: 3) - (5 * x1) ] = 1)))"
r:=b evaluate. "an OrderedCollection(
an IntervalBox([2.3300587395679817,2.3300587395679826]) 
an IntervalBox([-0.2016396757234047,-0.20163967572340466]) 
an IntervalBox([-2.1284190638445777,-2.128419063844577]))"
"the same 3 solutions for the equation as above." 
"or for an approximate Number result:"
r collect:[:i| i first midPoint ].
 "an OrderedCollection(2.330058739567982 -0.2016396757234047 -2.1284190638445772)"

Lets have a look at a more complicated 2-dimensional example. At http://icwww.epfl.ch/~sam/Coconut-benchs/Transcendental/descartesfolium.mod i found a problem that obviously can't be solved by every program.
"Domains
var X >= -100, <= 100;
var Y >= -100, <= 100;

subject to cons1 : (X^2)/Y + (Y^2)/X = 2;
cons2 : Y = exp(-X);"

box:=IntervalBox  with: (RealInterval inf: -100 sup: 100) with: (RealInterval inf: -100 sup: 100)."-->
an IntervalBox([-100,100] [-100,100])"
f1:=[:x|(x at:1)squared / (x at:2) +((x at:2) squared /(x at:1))].
(cons1:=Constraint block: f1) equalToNumber:  2.
cons1."--> a Constraint([ :x | (x at: 1) squared / (x at: 2) + ((x at: 2) squared / (x at: 1)) ] = 2)"
f2:=[:x|(x at:1)negated exp - (x at:2)].
(cons2:=Constraint block: f2) equalToNumber:  0.
cons2."--> a Constraint([ :x | (x at: 1) negated exp - (x at: 2) ] = 0)" 
cr:=ConstraintRefinement constraints: {cons1. cons2} box: box."--> a ConstraintRefinement(box: an IntervalBox([-100,100] [-100,100]) constraints: a ConstraintsCollection(a Constraint([ :x | (x at: 1) squared / (x at: 2) + ((x at: 2) squared / (x at: 1)) ] = 2) a Constraint([ :x | (x at: 1) negated exp - (x at: 2) ] = 0)))"
"attention, on my netbook this took almost a minute:"
cr evaluate. "an OrderedCollection(an IntervalBox(
[0.29456271159628045,0.2945627115962816] 
[0.7448572336907798,0.7448572336907809]) an IntervalBox(
[0.8684182784952086,0.8684182784952094] 
[0.41961473827445583,0.4196147382744562]))"
"iow there seem to exist two solutions
x~0.29456271159628, y~0.74485723369078
and
x~0.86841827849521, y~0.419614738274456"

But with constraints you won't always get those simple solutions. Lets take this equation:

     2 = x² + y²

We can make a Constraint like this:

f1:=[ :x | (x at: 1) squared  + (x at: 2) squared ].
c1:=Constraint block: f1.
c1 equalToNumber:2.
box:=IntervalBox  with: (RealInterval inf: -3 sup: 3) with: (RealInterval inf: -3 sup: 3) .
cr:=ConstraintRefinement constraints: c1 box: box.

Here we would get so many solutions that we have to raise the number of possible iterations and lower the desired precision. fast:true makes sure that BoxSplicing does not run too often.

cr maximumIterations: 2000.
cr desiredPrecision: 0.02.
cr fast:true.
[cr evaluate]timeToRun ."--> 0:00:00:53.761"
"that is still acceptable:"
cr desiredPrecision."--> true" "this works like in SkelboeMoore2"
cr iterations."--> 1948"
cr result size."--> 832"

The result size is too high to look at it 'by hand'. Therefore we have the BoxVisualizer:

v:=(BoxVisualizer result: cr result).
v show .

bv1
Of course '2 = x² + y²' just describes a circle; if you zoom in, you can see that the solution is made up of a lot of little boxes.
bv2
The BoxVisualizer also comes in handy for debugging purposes, for example for estimating what number of maximumIterations to use and such things.

Or let's take a look at a problem where you have several constraints:

     x² + y² ≤ 2
     ( 1.2 abs(x) - 0.62 )² + (y - 0.48)² ≥ 0.13
     (x/10)² + ( 2.5 + 3 y - x² )² ≥ 0.01
and
     -3 ≤ (x,y) ≤ 3

Let's make the constraints:

f1:=[ :x | (x at: 1) squared + (x at: 2) squared ].
f2 :=[ :x | ((x at: 1)abs * 1.2 - 0.62) squared  + ((x at: 2)-0.48) squared ].
f3:=[ :x | ((x at: 1)/10) squared  + ((x at: 2)*3+2.5-(x at: 1)squared) squared ].
c1:=Constraint block: f1.
c1 max:2.
c2:=Constraint block: f2.
c2 min:0.13. 
c3:=Constraint block: f3.
c3 min: 0.01.
box:=IntervalBox  with: (-3 hull:3) with: (-3 hull:3) .

You can enter those constraints simply with any collection:

cr:=ConstraintRefinement constraints: {c1.c2.c3} box: box.
cr maximumIterations: 3000.
cr desiredPrecision: 0.05. "less restrictive precision, so that this does not take too long"
cr fast:true.
[cr evaluate]timeToRun ." 0:00:02:18.481" "phhh"
cr iterations." 1375"
cr desiredPrecision." true"
cr result size." 1102"

And you can look at the result and change the visualization this way:

v:=(BoxVisualizer result: cr result).
v backgroundColor: Color black .
v color:Color yellow.
v withBorder: false.
v show.

bv3

If you have a more than 2-dimensional output, BoxVisualizer by default shows the projection on the first 2 dimensions but you can choose the displayed dimensions with dimensions: anArrayOfDimensions, for example you can exchange the output dimensions of that picture with:

v dimensions: #(2 1).
"the visual effect will be applied immediately"

bv4

Optimization with Constraints

For this we have ConstraintsMinimizer1. This is a SubObject of SkelboeMoore1, in other words it works exactly like that algo, simply with Constraints added. Hence some of the tricks I'll show here also apply to the Skelboe-Moore algos.

Let's try a simple problem, say we search for the lowest point within the closed unit disk centered at the origin:

minimize
     f(x,y) = y
with constraint
     x2 + y2 ≤ 1

The minimum obviously is at 0@(-1).

f1:=[ :x | (x at: 1) squared  + (x at: 2) squared ].
c1:=Constraint block: f1.
c1 max:1.
box:=IntervalBox  with: (RealInterval inf: -3 sup: 3) with: (RealInterval inf: -3 sup: 3) .
c1. "-->a Constraint([ :x | (x at: 1) squared + (x at: 2) squared ] <= 1)"
a:=ConstraintsMinimizer1 function: [:x|x second] box: box constraints: {c1}.
[r:=a evaluate]timeToRun."-->0:00:00:01.392"
r. "-->an Array(an IntervalBox(an IntervalUnion([-0.00633800165035087,0.006355185716093608]u[0.007807221955565399,0.009085458267685547]) [-1.000001083542924,-0.9999587780455154]) [-1.000001083542924,-0.9999587780455154])"
r first midPoint. "-->#(0.0013737283086673385 -0.9999799307942197)"

ConstraintsMinimizer1 returns, exactly like SkelboeMoore1, an Array with the argument of the function that minimizes that function, and the return-value of that function. Hence it says with
x = an IntervalUnion([-0.00633,0.00635]u[0.0078,0.0090])
y = [-1.000001083542924,-0.9999587780455154]
you get the minimum of y. The true minimum at 0@(-1) is obviously encompassed, i just wonder why the interval [0.007807221955565399,0.009085458267685547] for x was not thrown away. It looks a bit ugly. Obviously the precision of the calculation is not good enough. Since the default precision in SkelboeMoore1 is set relatively high it is is enough to set maximumIterations higher to get a more precise result.

a maximumIterations:40000.
[r:=a evaluate]timeToRun ."-->0:00:00:19.902"
r. "-->an Array(an IntervalBox(an IntervalUnion([-0.00038281297364662655,0.00038281297364662655]) [-1.000000091987346,-0.9999999267281785]) [-1.000000091987346,-0.9999999267281785])"
r first midPoint ."-->#(0.0 -1.0000000093577621)"

This is better:
     x = [-0.00038281297364662655,0.00038281297364662655]
     y = [-1.000000091987346,-0.9999999267281785]
and the midPoint is also nearer to the correct value, although it does not exactly fulfill the constraint. But you can't expect that, you can only expect that the real minimum lies somewhere in the intervals, in other words the midPoint is only an approximation.

Now we will look at a problem where the minimum has infinitely many solutions (from the book Introduction to Interval analysis)

minimize
     f(x) = (| x12 + x22 - 1 | + 0.001 ) * | x12 + x22 - 0.25 |
with constraint
     max{ ( 1 - max{ | x1 | / 0.6 , | x2 | / 0.25 } ) , ( 1 - max{ | x1 | / 0.25 , | x2 - 0.4 | / 0.3 } ) } ≤ 0
and
     -1.2 ≤ ( x1, x2 ) ≤ 1.2

f1:=[ :x | (((x at: 1) squared  + (x at: 2) squared -1)abs +0.0001)*((x at: 1) squared  + (x at: 2) squared -0.25)abs].
f2:=[ :x | (1 - ((x at: 1) abs  / 0.6 max: (x at: 2) abs  / 0.25))max:(1 - ((x at: 1) abs  / 0.25 max: ((x at: 2)-0.4) abs  / 0.3))].
c1:=Constraint block: f2.
box:=IntervalBox  with: (RealInterval inf: -1.2 sup: 1.2) with: (RealInterval inf: -1.2 sup: 1.2) .
a:=ConstraintsMinimizer1 function: f1 box: box constraints: {c1}.
a iterations.
a evaluate."-->an Array(an IntervalBox(an IntervalUnion([-0.433554242093114,0.433554242093114]) 
an IntervalUnion([-0.5018130538924789,-0.24955148704901042]u[0.24955148704901042,0.433554242093114])) 
[0.0,0.0038966431729056243])"
"if you have difficulties visualizing a result like this, you can
always use BoxVisualizer:"
v:=BoxVisualizer result: r first .
v show.

bv3

This looks like sloppySimplification, which lays hulls around the result to make the output easier to digest, would have been a bit too sloppy.

a useSloppySimplification: false.
[r:=a evaluate]timeToRun . "-->0:00:00:02.256" 
r first size. "-->887" 
v:=(BoxVisualizer result: r first ).
v color:Color red.
v withBorder:false.
v backgroundColor: Color gray.
v show.

bv3

Well, this looks like a smiley face, as advertised in the mentioned book. Ok, I'll leave it at that. Well, perhaps one remark why my implementation of this algo is essentially wrong. A few further test-functions are here

Final Remarks

These algos can of course be made more efficient and faster. I did not do it, since I think that by keeping them in their original stupid simple version, they make up nice little bricks of this experimental kit. They are easily understandable and comparable to the explanations in introductory books and papers for interval arithmetic like eg (first two links are written by R. A. Moore, one of the pioneers of interval arithmetic) :
http://interval.louisiana.edu/Moores_early_papers/Moore_in_Rall_V1.pdf
http://www-sbras.nsc.ru/interval/Library/InteBooks/IntroIntervAn.pdf
http://arxiv.org/pdf/physics/0302034v2.pdf
http://www.mat.univie.ac.at/~neum/ms/glopt03.pdf
http://pages.cpsc.ucalgary.ca/~rokne/global_book.pdf
http://eprints.ma.man.ac.uk/1204/01/narep416.pdf
And of course a whole lot of other things on this subject can be found in the internet, especially at:
http://www.cs.utep.edu/interval-comp/

I should perhaps mention that this project generally works nicely together with ArbitraryPrecisionFloat. You simply lay a hull around two ArbitraryPrecisionFloats and do your calculations with the resulting RealInterval. Two things I noticed are, that binary operations on ArbitraryPrecisionFloats with RealIntervals as arguments of course don't work, in that case you could use Numbers like Fractions or so. Another problem can be that the RealInterval implementation of some functions (eg trigonometric functions) uses Float constants which necessarily lowers the precision of calculation with more precise ArbitraryPrecisionFloats.

RealInterval works also nicely together with AutomaticDifferentiation, I could detect no incompatibilities. A small example is here.

Generally, if one has to deal with a special problem, one should not forget that adjusting the interval extension by hand can sometimes help and simplify the problem.