The derived class (based on the SPnode class) used to reimplement the set computations described in AIA2001 is SPMinimalnode. A pointer to an SPMinimalnode is aliased as MinimalSubPaving. This class is specifically for minimal trees of nodes (representing minimal subpavings). As discussed above, a tree (or the subpaving it represents) is minimal if it has no sibling leaves; that is, no node has more than one leaf-child.

Equivalence to the AIASPnode class

In order to demonstrate that this implementation is equivalent to that of AIA2001, we reimplement the set inversion and image evaluation functions and using our new SPMinimalnode class.

For examples using the new SPMinimalnode class for set computation and demonstrating the equivalence of the results to those in Examples from Applied Interval Analysis using AIASPnodes see Example of set computation with SPMinimalnodes

Example of set computation with SPMinimalnodes

SPMinimalnodes are re-implemention in C++ of a subpaving-as-a-binary-tree object based on [AIA2001]. SPMinimalnodes behave exactly as AIASPnodes do but we highlight here some of the minor changes that we have made the way that this class is coded and used.

Example 3.4

In this example show how example 3.4 from AIA2001, pp. 61-63, can be performed using the SPMinimalnode class (for the same example using AIASPnodes, code closely based on AIA2001, see Example 3.4 using AIAPSnodes). This example uses our implementations of both Sivia (see Set inversion using interval analysis) for set inversion and ImageSp (see Image evaluation) for image evaluation. The functions we use here are Sivia (PIBT BoolTest, const SPMinimalnode * const toInvert, SPMinimalnode * const search, const double eps) and ImageSp(PIVF f, SPMinimalnode* spn, double eps).

Our implementation of this example is in the file NewExm_3_4.cpp, which has header file NewExm_3_4.hpp. The header file is

#ifndef ___SPTEST_EXM3_4_HPP__
#define ___SPTEST_EXM3_4_HPP__

// include all headers for using spnodes
#include "spnodeall.hpp"

#include <time.h>
#include <fstream>

using namespace cxsc;
using namespace std;
using namespace subpavings;


#endif

We include headers and libraries we want to be able to use: <time.h> for a clock, <fstream> for file output, and suppavings.hpp for our class declarations. We also specify the namespaces we want to be able to use without qualification. Note that subpavings is now a namespace. If we did not specify that we are using this namespace we would have to qualify all references to entities declared within it, for example subpavings::MinimalSubPaving (MinimalSubPaving is an alias for a pointer to an SPMinimalnode).

In NewExm_3_4.cpp we include the header file

// implementation of AIA Example 3.4 using class SPMinimalnode
#include "NewExm_3_4.hpp"

Then we define some functions used in our main program.

The first is a boolean interval test, IBT_ex3_4.

The interval boolean test is the key to set inversion with interval analysis. Recall that set inversion is the computation of a reciprocal image X = f^-1(Y) where Y is regular subpaving of R^m (see Set inversion using interval analysis). Y is the subpaving we want to invert. The interval boolean test takes a box in 'x-space' and tests whether the image of the box in 'y-space' ie the image under the inclusion function [f], is in the subpaving Y. The interval boolean test can return one of the special interval boolean types BI_TRUE (the image of the box is inside Y), BI_FALSE (the image of the box is outside Y), or BI_INDET (indeterminate: the image of the box [f]([x] is partly in and partly out of the subpaving Y, ie overlaps the boundary of Y).

This test IBT_ex3_4 tests the image of a box x for inclusion in a subpaving represented by the SPMinimalnode pointer spn. Values for x and spn are given by the caller as the function arguments. f is specified in the test as f: R² -> R such that f(x₁, x₂) = x₁⁴ - x₁² + 4x₂².

// specifications of example interval boolean tests
// The boolean interval test can return BI_TRUE, BI_FALSE, or BI_INDET
BOOL_INTERVAL IBT_ex3_4(const ivector& x, const SPMinimalnode * const spn)
{
    ivector Img(1);
  // the inclusion function image of the given interval x
    Img[1] = sqr(sqr(x[1])) - sqr(x[1]) + 4*sqr(x[2]);

    return Img<=spn;

}

Note:: Interval boolean tests in the subpavings namespace take a pointer to subpaving we want to invert as a parameter. This avoids the use of global subpaving variables in these tests as in example 3.4 with AIASPnodes and makes the functions somewhat clearer to interpret.

Then we have another boolean interval test, IBTFinverse_ex3_4. Again this specifies a function f, this time

f: R² -> R², f₁(x₁, x₂) = (x₁ - 1)² -1 + x₂, f₂(x₁, x₂) = -x₁² + (x₂ - 1)². Values for x, the box whose image we test, and spn, a pointer to the subpaving to invert, are given by the caller as the function arguments.

// A boolean interval test to illustrate SIVIA being used
// to evaluate the domain corresponding to an
// image subpaving spn and a function f
BOOL_INTERVAL IBTFinverse_ex3_4(const ivector& x, const SPMinimalnode * const spn)
{

    //f1(x1, x2) = (x1-1)^2 - 1+ x2
    //f2(x1, x2) = -(x1^2) + (x2-1)^2

    ivector Img(2);

    Img[1] = sqr(x[1]) - 2*x[1] + x[2];
    Img[2] = -sqr(x[1]) + sqr(x[2]) - 2*x[2] + 1;

    return (Img<=spn);
}

// end specification of example boolean interval tests

The next function is an interval vector function, IVF_ex3_4. An interval vector function returns the interval vector image of an interval vector x under f for f as specified in the function and x supplied as the function argument. In this case f is exactly the same as in the boolean interval test IBTFinverse_ex3_4 above.

// specification of example interval vector function
// (ie, vector inclusion function)

ivector IVF_ex3_4(const ivector& x)
{
    // example in 2-d space R2
    // for f: R2 -> R2
    //f1(x1, x2) = (x1-1)^2 - 1+ x2
    //f2(x1, x2) = -(x1^2) + (x2-1)^2

    ivector Img(2);

    Img[1] = sqr(x[1]) - 2*x[1] + x[2];
    Img[2] = -sqr(x[1]) + sqr(x[2]) - 2*x[2] + 1;

    return (Img);
}

// end specification of example interval vector functions

Then we start the main program, declare some variables, and set up an initial search box [-3.0 3.0]². This is used to initialise a newed SPMinimalnode in dynamic memory, with A as a MinimalSubPaving or pointer to a SPMinimalnode. Creating the node object in dynamic memory allows us to use pointers to it outside the scope of the function which created the object

int main()
{
    double prec;
    clock_t start, end;

    ivector x(2);
    x[1] = interval(-3.0,3.0);
    x[2] = interval(-3.0,3.0);

    MinimalSubPaving A;
    A = new SPMinimalnode(x);

We then set up the subpaving we want to invert (again, in dynamic memory).

    ivector t(1);
    // a restriction on the range of the function - we want to find the x
    t[1] = interval(-0.1,0.1);
                    // for which f(x) is in this interval
                    // ie acts like a simple MinimalSubPaving for SIVIA

And then prepare to use SIVIA to do the set inversion. We make some standard output to say what we are doing and ask the user to enter a precision.

    // Using SIVIA for set inversion

    //find an MinimalSubPaving characterisation Sc5 as in example 3.4

    cout << "Characterization of the set Sc5={(x1,x2) | -0.1 "
       << "<= x1^4-x1^2+4x2^2 <= 0.1 }" << endl;
    cout << "Enter a precision (between 1 and 0.001): ";

Now we use the SIVIA algorithm (see Set inversion using interval analysis) to find a 'subpaving characterisation' Sc5 of the set we want with Sivia(PIBT BoolTest, const SPMinimalnode * const toInvert, SPMinimalnode * const search, const double eps). A PIBT is a pointer to a boolean interval test, in this case IBT_ex3_4. The set we want is defined by the combination of the SubPaving ToInvert and the f specifed in the boolean interval test IBT_ex3_4. We are finding a characterisation for the set (x₁, x₂) such that f(x₁, x₂) = x₁⁴ - x₁² + 4x₂² is in ToInvert, which is the interval [-0.1 0.1].

    start = clock();

    // when we start we give A a box big enough to guarantee to contain
    // the characterisation of Sc

    MinimalSubPaving Sc5 = Sivia(IBT_ex3_4, ToInvert, A, prec);
    end = clock();

The SIVIA algorithm works by taking an initial search box, testing it, rejecting those returning BI_FALSE, including those returning BI_TRUE in the subpaving it is building, and bisecting again if the interval boolean test returns BI_INDET, and recursively sending each subbox to SIVIA to be tested similarly. Each subbox continues to be bisected until either the test returns a clear BI_TRUE or BI_FALSE or the test is BI_INDET but the width of the box is the value given for eps, which means that it is thin enough to be included in the subpaving to be returned. The eps parameter (precision) provides a 'stopping rule' prevents SIVIA recursing endlessly. A larger value for eps will mean a thicker uncertainty layer in the outer subpaving of the area we seek to characterise.

Note:: This version of Sivia takes a parameter for the SubPaving to invert because it must pass this on to the boolean interval test.

Sc5 now points to the subpaving characterisation we have created. We say 'subpaving characterisation' because the subpaving is an outer subpaving of the actual set we want.

We print the computing time, volume, and number of leaves of the subpaving we create to standard output.

    cout << "Computing time : "
       << ((static_cast<double>(end - start)) / CLOCKS_PER_SEC) << " s."<< endl;
    cout << "Volume: " << spVolume(Sc5) << endl;
    cout << "The number of leaves is " << spLeaves(Sc5) << endl;

And send the output to a txt file.

    // To realize a file output of the MinimalSubPaving Sc5
    ofstream os5("New3_4a.txt");         // Filename
    os5 << 2 << endl;                    // Dimension of the MinimalSubPaving
    os5 << interval(-3.0,3.0) << " "         // Root box
        << interval(-3.0,3.0) << " " << endl;
    os5 << "Precision is " << prec << endl; // Precision used
    os5 << Sc5 << endl;                   // MinimalSubPaving itself
    cout << "The output MinimalSubPaving has been written to New3_4a.txt"

Now we want to find the image of this set using the function specified in IVF_ex3_4. f: R² -> R², f₁(x₁, x₂) = (x₁ - 1)² -1 + x₂, f₂(x₁, x₂) = -x₁² + (x₂ - 1)².

    // Using ImageSp to find the image of Sc5 using
    //  f1(x) = (x1-1)^2 - 1+ x2
    // "    f2(x) = -(x1^2) + (x2-1)^2


    cout << "Characterization of the set Sc6=f(Sc5)" << endl
        << " with Sc5 from our previous example and "<< endl
        << "with f1(x) = (x1-1)^2 - 1 +x2," << endl
        << "     f2(x) = -(x1^2) + (x2-1)^2" << endl;
    cout << "by realizing the image of f by ImageSp" << endl;
    cout << "Enter a precision (between 1 and 0.01): ";

This function is not invertible in the usual sense and we have to use ImageSp to find the image. The version of ImageSp used is ImageSp(PIVF f, SPMinimalnode* spn, double eps). We supply the interval vector function IVF_ex3_4, the subpaving Sc5 and the precision input by the user as arguments in the call to ImageSp.

    start = clock();

    // use Image SP to find a characterisation of the image of Sc5
  // using the function in IVF
    MinimalSubPaving Sc6 = ImageSp(IVF_ex3_4, Sc5, prec);

    end = clock();

Sc6 now points to the root node of a tree representing the subpaving constructed with ImageSp.

The ImageSp algorithm works by taking an initial subpaving (represented by Sc5 in this case), mincing it up into a fine subpaving where every box has width less than the precision specified, and finding the image of each of these boxes with the specified interval vector test. It then forms a minimal regular subpaving which covers the union of all these image boxes, again with precision as specified.

Sc6 now points to a subpaving representation in 'y-space', the image of Sc5 under the function f in IVF_ex3_4.

We report computing time, volume and number of leaves and send the suppaving output to a txt file.

    cout << "Computing time : "
       << ((static_cast<double>(end - start)) / CLOCKS_PER_SEC) << " s."<< endl;

    cout << "The volume is " << spVolume(Sc6) << endl;
    cout << "The number of leaves is " << spLeaves(Sc6) << endl;

    // To realize a file output of the MinimalSubPaving Sc
    ofstream os6("New3_4b.txt");            // Filename
    os6 << 2 << endl;                    // Dimension of the MinimalSubPaving
    os6 << interval(-3.0,3.0) << " "         // Domain MinimalSubPaving
        << interval(-3.0,3.0) << " " << endl;
    os6 << "Precision is " << prec << endl; // Precision used
    os6 << Sc6 << endl;                   // Image MinimalSubPaving itself

    cout << "The output MinimalSubPaving has been written to New3_4b.txt"

Subpaving A, which will have been minced and mangled in the process of forming Sc6, is now deleted and remade to provide another initial search box.

    //remake A
    delete A;

    x[1] = interval(-5.0,5.0);
    x[2] = interval(-5.0,5.0);

    A = new SPMinimalnode(x);

What happens if we now use Sivia on our ImageSp-created image Sc6 in 'y-space', inverting our image to get back to 'x-space'?

    // this uses the MinimalSubPaving Sc6 created by the above example
    // create an MinimalSubPaving Sc7 which contains inverse_f(Sc6)

    cout << "Characterization of the set Sc7=f-1(Sc6)" << endl
        << "with f as above  and Sc6 as above " << endl;
    cout << "by realizing the inversion of f by Sivia" << endl;
    cout << "Enter a precision (between 1 and 0.01): ";

The boolean interval test IBTFinverse_ex3_4 specifies the function. We pass a pointer to this function to Sivia by supplying it as the value for the PIBT parameter. We give Sc6 as representing the MinimalSubPaving to invert and A as reprsenting the initial search box, and we also give the precision supplied by user.

    start = clock();
    MinimalSubPaving Sc7 = Sivia(IBTFinverse_ex3_4, Sc6, A, prec);
    end = clock();

Sc7 is now points to 'x-space' subpaving characterisation of the reciprocal image of Sc6, which was in turn a subpaving characterisation of Sc5.

We report computing time, volume and number of leaves and output the subpaving to a txt file.

    cout << "Computing time : "
       << ((static_cast<double>(end - start)) / CLOCKS_PER_SEC) << " s."<< endl;
    cout << "Volume: " << spVolume(Sc7) << endl;
    cout << "The number of leaves is " << spLeaves(Sc7) << endl;

    // To realize a file output of the MinimalSubPaving Sc
    ofstream os7("New3_4c.txt");            // Filename
    os7 << 2 << endl;                    // Dimension of the MinimalSubPaving
    os7 << interval(-5.0,5.0) << " "         // Root box
        << interval(-5.0,5.0) << " " << endl;
    os7 << "Precision is " << prec << endl; // Precision used
    os7 << Sc7 << endl;                   // MinimalSubPaving itself

    cout << "The output MinimalSubPaving has been written to New3_4c.txt"

We delete our subpavings created in dynamic memory before we end the program

    delete A;
    delete Sc5;
    delete Sc6;
    delete Sc7;

    return 0;
}

The subpavings produced by this program run using precision 0.05 in each case are shown below.

As well as showing the initial subpaving, the image subpaving, and the reciprocal image of the image, we capture an intermediate step in the process of creating Sc6 from Sc5 with ImageSp. This is the evaluation step, where we have a large set of (possibly overlapping) image boxes formed from all the minced up subboxes of the initial box [-3.0 3.0]² chopped up so that each one is less than 0.05 wide. This set can be compared to the Sc6, the regular minimal subpaving characterisation of the image.

The most interesting comparison though is between the initial subpaving and the subpaving for a reciprocal image of the image subpaving. The initial set (characterised by Sc5) is in there but the result is fatter due to error accumulation in the process of going to 'y-space' and back, and the final subpaving has additional parts which appear because f(.) is only invertible in the set-theoretic sense (AIA2001, pp. 63)

Results for Example 3.4 using precision 0.05

We can rerun the program to explore the effect of the precision. The image below was created using a precision of 0.001 for the initial subpaving and then precisions of 0.01 to create the image of this in 'y-space' and and 0.01 to go back again to the reciprocal image of the image in 'x-space'. The evaluated images of the initial subpaving after it has been minced up very finely (precision 0.001) are also shown. The reciprocal image of the image is a much better comparison to the initial subpaving as the errors accumulated are smaller, but the additional parts due to are again apparent - they are caused by the process, not by the precision to which we carry it out.

Results for Example 3.4 using precision 0.001 for initial paving, then 0.01