MRS  1.0
A C++ Class Library for Statistical Set Processing
Implementation of Applied Interval Analysis procedures under C-XSC

Regular subpavings and set computation

[AIA2001] show how subpavings can approximate compact sets in a guaranteed way and how computation on subpavings allows approximate computation on these compact sets. For any full compact set X we can find a subpaving 'lower bound' of X (which can be thought of an an 'inner paving' of X) and a subpaving 'upper bound' of X (or 'outer paving' of X) as 'close' to X as we desire (distance and ordering of sets and hence what is meant by 'close' has to be defined of course). X is enclosed by the inner and outer paving. [AIA2001, pp. 48-49]

To get 'close' to the set we want to enclose, we can take a large box (the root node of the binary tree representing the subpaving) which we know encloses the set and then progressively subdivide the box and the subboxes, etc, checking to see if we can discard any subbox (prune that branch) to get a subpaving/tree representation of the subpaving which meets our criteria for 'closeness' to our target compact set.


Set inversion using interval analysis

Set inversion is the computation of a reciprocal image X = f-1(Y) where Y is regular subpaving of Rm. [AIA2001] develop an algorithm SIVIA (Set Inverter Via Interval Analysis) for finding an outer subpaving of X to a specified level of precision. Starting with a large search box [x](0) to which the outer subpaving of X is guaranteed to belong and given a inclusion function [f]([x]) for f, we form a subpaving by progressively bisecting and then testing the resulting subboxes to see if they should be included in our solution subpaving.

The test for any box [x] is a comparison of [f]([x]) to Y. This test has three possible outcomes

  1. If [f]([x]) has an empty intersection with Y then [x] does not belong to X and can be cut off from the solution tree
  2. If [f]([x]) is entirely in Y then [x] belongs to the solution and is kept as part of the solution tree
  3. If [f]([x]) has a non-empty intersection with Y but is not entirely in Y then [x] is undetermined. If [x] has width greater than the specified precision parameter it is bisected and each subbox (child in the tree structure) is then tested. If [x] has width less than or equal to the specified precision parameter it included in the solution, ie the suppaving outer bound of the set X [AIA2001, pp. 55-56].

Thus the solution, the outer subpaving of the set X, includes the 'uncertainty' layer of boxes that cross the boundaries of the set but whose width is small enough to satisfy our precision criteria. We can make the uncertainty layer thinner by having a smaller precision and subdividing and testing the boxes in that layer further [AIA2001, p. 56].

An implementation of SIVIA using C++ is described in [AIA2001, pp. 339-342]. This includes an extension of the usual booleans TRUE and FALSE to a set of interval booleans which includes notion of indeterminate.

SIVIA can also be used to evaluate an outer subpaving of the image of set by a function provided that the function is invertible in the usual sense. This involves specifying the inclusion function [f-1] for f-1 the inverse of f and taking as an initial search box some box guaranteed to contain the image of the set.


Image evaluation

Computation of the direct image of a subpaving Y = f(X) where X is a regular subpaving of Rn (image evaluation) is more complex in the case where the function is not invertible in the usual sense. AIA2001 develop an algorithm ImageSP for finding an outer subpaving of Y to a specified level of precision.

Given a subpaving X for which we seek an outer approximation of the image and an inclusion function [f]([x]) for f, ImageSP proceeds as follows:

  1. Mince X (recursively bisect each box) until each box in X has width less than the specified level of precision. (X will no longer be a minimal subpaving.)
  2. For each [x] in the minced X, add [f]([x]) to a list of image boxes U and compute the interval hull of the union of all these image boxes.
  3. If [f]([x]) Merge these boxes into a single, minimal, subpaving the root of which corresponds to the hull and which only contains boxes with width lower than the specified level of precision.

An implementation of ImageSP using C++ is described in [AIA2001, pp. 342-347].

We have reimplemented the algorithms for set inversion and image evaluation using the C-XSC interval library with as little change to the structure and code used in [Applied Interval Analysis, Springer, 2001] as possible, other than that necessitated by the use of C-XSC. Our class is called AIASPnode (a pointer to an AIASPnode is aliased as AIASubPaving)


Examples from Applied Interval Analysis using AIASPnodes

AIA2001 show how computation on subpavings allows approximate computation on compact sets. Their SUBPAVINGS class and various supporting procedures are used to implement their two main algorithms:

  • SIVIA (Set Inversion Via Interval Analysis). Set inversion is the computation of a reciprocal image X = f-1(Y) where Y is regular subpaving of Rm
  • ImageSp (Image evaluation). Computation of the direct image of a subpaving Y = f(X) where X is a regular subpaving of Rn

The following examples show how we have replicated the SUBPAVING class and the implementations of SIVIA and ImageSp created by AIA2001 and available on their website http://www.lss.supelec.fr/books/intervals/)

Our implementation uses the C-XSC library but otherwise makes as little alteration to the structure and code used in Applied Interval Analysis as possible, other than that necessitated by the use of C-XSC. Our class is called AIASPnode (a pointer to an AIASPnode is aliased as AIASubPaving).

The AIASPnode class declarations and inline definitions are in the header file AIAsubpaving.hpp. Other definitions are in the file AIAsubpaving.cpp

Exercise 11.33

See AIA2001, p. 342.

Exercise 11.33 explores the use of SIVIA and the AIASPnode class for set inversion. In particular, this example demonstrates that it is possible to use to evaluate the direct image of a set by a function, provided that the function is invertible in the usual sense.

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

Exr_11_33.hpp shows the includes for this example

#include "AIAsubpaving.hpp"
#include <time.h>
#include <fstream>

The file AIAsubpaving.hpp is included in order to be able to use the AIASPnode class it declares. <time.h> is included in order to use the clock method and get information on the time take to run the example. <fstream> is included in order to be able to output data to a file.

Turning now to Exm_11_33.cpp,

First we include the header file

#include "Exr_11_33.hpp"

Then we declare some AIASubPavings as global (AIASubPaving is declared as an alias for a pointer to an AIASPnode in the AIAsubpavings.hpp header file). Declaring them as global means that they are available to all the functions in the file (they have file scope) and so they do not need to be passed as parameters to the functions that use them.

// These AIASubPavings are declared as global
AIASubPaving Sc, Sc1, Sc2;

Remarks:
Use of globals is normally discouraged in good programming,

Then we specify the interval boolean tests we are going to use in the example. 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 Rm (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).

The first test IBTAnnular tests a 2-d box in IR2 for inclusion in the set corresponding to the area between circles centred at the origin and with radii 1 and 2. The test is performed on the interval image [f]([x]) of the box [x] supplied as the function argument. The test itself specifies the inclusion function f [f] (x[1]2 + x[2]2). In this case the test also specifies the subpaving to invert - in this case this is a very simple subpaving being just a single interval [1 2]. The test can return BI_TRUE, BI_FALSE, or BI_INDET.

// specifications of example interval boolean tests
// The boolean interval test can return BI_TRUE, BI_FALSE, or BI_INDET

AIA_BOOL_INTERVAL IBTAnnular(const ivector& x)
{
  // here we test a 2-d box for inclusion in the area between circles centred 
  // on the origin with radii 1 and 2
  interval ToInvert(1.0,2.0),Temp;
  interval Img = sqr(x[1]) + sqr(x[2]);

  if (!Intersection(Temp,Img,ToInvert)) return BI_FALSE;
  if ( Img<=ToInvert ) return BI_TRUE;

  return BI_INDET;
}

The second test IBTFdirect is also an interval boolean test but on first sight looks rather different to IBTAnnular. Taking it apart, we can how it does the same thing. The test takes a ivector argument and calculates the interval image of this using an inclusion function. The inclusion function is again specified within the test. In this case it is in fact the inverse of an invertible function: this test will be used to run SIVIA 'backwards'. Note that the subpaving we are inverting is one of the global AIASubPaving variables, Sc. Because the subpaving is not a simple interval as in the previous test, the test uses the AIASPnode::operator<=(const ivector&, AIASubPaving), which returns an AIA_BOOL_INTERVAL type just like the test above.

AIA_BOOL_INTERVAL IBTFdirect(const ivector& x)
{
  // A boolean interval test to illustrate SIVIA being used to evaluate the 
  // direct image of a set by a function provided the function is invertible
  // ie SIVIA will invert inverse_f

  interval Temp;
  ivector Img(2);

  // taking the function f : (x1,x2) -> (2x1-x2,-x1+2x2)
  // f is invertible and the inverse is 
  // inverse_f : (x1,x2) -> (2x1 + x2,x1+2x2)/3
  // this is an inclusion function for inverse_f

  Img[1] = (2.0*x[1]+x[2]) / 3.0;
  Img[2] = (x[1]+2.0*x[2]) / 3.0;

  return (Img<=Sc); // Sc is an AIASubPaving set up by the first example
}

The third test IBTFinverse is similar again. The function f is the same function as that for which we used the inverse in the test IBTFdirect above. The subpaving we test for inclusion in is another of the globals, Sc1, and again we use the AIASPnode::operator<=(const ivector&, AIASubPaving) to returns an AIA_BOOL_INTERVAL.

AIA_BOOL_INTERVAL IBTFinverse(const ivector& x)
{
  interval Temp;
  ivector Img(2);

  // the function f is f : (x1,x2) -> (2x1-x2,-x1+2x2)

  Img[1] = 2.0*x[1]-x[2];
  Img[2] = -x[1]+2.0*x[2];

  return (Img<=Sc1);
}

// end specification of example boolean interval tests

Now we say that we are using the std and cxsc namespaces to avoid having to type cxsc::ivector or std::cout etc.

using namespace cxsc;
using namespace std;

In the main process we start by declare some of the variables we will be using

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

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

  AIASubPaving A;
  A = new AIASPnode(x);

A 2-dimensional interval vector x is set up which is the 2-d box [-5 5]2. This is used to initialise a newed AIASPnode in dynamic memory, with A as an AIASubPaving or pointer to an AIASPnode. Creating the node object in dynamic memory allows us to use pointers to it outside the scope of the function which created the object

Now we use the SIVIA algorithm (see Set inversion using interval analysis) to find a 'subpaving characterisation' Sc containing the area between circles centred on the origin with radii 1 and 2 in 2 dimensional space. We say 'subpaving characterisation' because the subpaving is an outer subpaving of the actual area between circles centred on the origin with radii 1 and 2.

  // Using SIVIA for set inversion
  //find an AIASubPaving characterisation Sc containing the area between 
  // circles centred on the origin
  // with radii 1 and 2 (in 2 dimensional space)
  cout << "Characterization of the set Sc={(x1,x2) | 1 <= x1^2+x2^2 <= 2 }" 
       << endl;
  cout << "Enter a precision (between 1 and 0.001): ";

A description of what we are doing is sent to standard output, The program asks the user to input a precision between 1 and 0.001.

The clock is started and Sivia(AIA_PIBT BoolTest, AIASubPaving A, double eps) is called. The AIA_PIBT (pointer to an interval boolean test) is IBTAnnular, as described above. A is provided as an initial search box, and prec is given as the argument for the eps parameter. The clock is stopped when Sivia has returned a value for Sc.

  start = clock();

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

  Sc = Sivia(IBTAnnular,A,prec);
  end = clock();

Sc now points to the root node of a tree representing the subpaving constructed with Sivia.

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 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.

We stop the clock after calling Sivia and report the computing time and the volume and number of leaf boxes of the subpaving.

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

Then we create an output file of giving the subpaving as a list of interval vectors. The file name is specified in the example as AIAannular.txt and the first three lines of the file give, respectively, the dimensions we are working in, the initial search box, and the precision used. The subpaving itself is then output. The format used for outputting the subpaving is specified in operator<<(std::ostream &, AIASubPaving)

  // To realize a file output of the AIASubPaving Sc
                    // Filename
  ofstream os("AIAannular.txt");
  os << 2 << endl;  // Dimension of the AIASubPaving
                    // Root box
  os << interval(-5.0,5.0) << " "
    << interval(-5.0,5.0) << " " << endl;
                    // Precision used
  os << "Precision is " << prec << endl;
  os << Sc << endl; // AIASubPaving itself
  cout << "The output AIASubPaving has been written to AIAannular.txt" 
       << endl << endl;

  // end of testing to find reciprocal image

  // the AIASubPaving Sc that we have created is the regular AIASubPaving 
  // that covers the set
  // X = {(x1,x2) in R2 | sqr(x1) + sqr(x2) is in [1,2]} 
  // (remember that it contains this area rather than being this area, 
  // and that eps has determined how small we go in the AIASubPaving 

This example now moves on to use this subpaving in further demonstrations of the SIVIA algorithm. First, we delete the AIASPnode that A currently points to and replace it with a new one, again providing an initial search box [-5.0 5.0]2

  // make a new AIASubPaving to provide an initial source box for the test next

Now we prepare to use SIVIA to find the direct image of an invertible function, which can be thought of as using SIVIA in reverse: normally SIVIA characterises the reciprocal image X = f-1(Y) where Y, the subpaving we want to invert, is regular subpaving of Rm. However, if we have X but want to find Y (or a characterisation for Y) and f is invertible so that we can specify an inclusion function [f-1] for f-1 then we can use this inclusion function to characterise Y. SIVIA can only be used to find a direct image under f when f is invertible in the normal sense because we need to be able to specify [f-1] in our interval boolean test.

A description of what we are doing is printed to standard output and the user again asked to enter a precision, and the clock is started

  // testing using SIVIA to find the direct image of an invertible function
  // remember that we are only finding some upper enclosure of the direct 
  // image really
  // Note that this example will use the AIASubPaving Sc we created 
  // above - see IBTFdirect

  // ie create an AIASubPaving Sc1 containing f(Sc), where Sc was found above

  cout << "Characterization of the set Sc1=f(Sc)" << endl
    << "with f1(x) = 2*x1-x2," << endl
    << "      f2(x) = -x1+2*x2," << endl;
  cout << "by realizing the inversion of f-1 by Sivia" << endl;
  cout << "Enter a precision (between 1 and 0.01): ";

Now we run Sivia as above, but the argument supplied for the interval boolean test parameter is IBTFdirect as given at the top of the file. This specifies [f-1] for f1(x) = 2x1 - x2, f2(x) = -x1 + 2x2. Recall that IBTFdirect specifies that the subpaving to invert is Sc, ie the subpaving we found using Sivia above.

  start = clock();
                    // Sc1 will be used by the following example
  Sc1 = Sivia(IBTFdirect,A,prec);
  end = clock();

Sc1 now points to the root node of a tree representing the subpaving constructed with Sivia.

The computing time, volume and number of leaves are reported and an output file produced.

  cout << "Computing time : " 
       << ((static_cast<double>(end - start)) / CLOCKS_PER_SEC) << " s."<< endl;
  cout << "Volume: " << Volume(Sc1) << endl;
  cout << "Number of leaves: " << NbLeaves(Sc1) << endl;

  // To realize a file output of the AIASubPaving Sc1
                    // Filename
  ofstream os1("AIAdirect.txt");
  os1 << 2 << endl; // Dimension of the AIASubPaving
                    // Root box
  os1 << interval(-5.0,5.0) << " "
    << interval(-5.0,5.0) << " " << endl;
                    // Precision used
  os1 << "Precision is " << prec << endl;
                    // AIASubPaving itself
  os1 << Sc1 << endl;

  cout << "The output AIASubPaving has been written to AIAdirect.txt" 
       << endl << endl;

  // end of example for finding direct image

We have now found a subpaving Sc in 'x-space' and used Sivia on an function inverse to get a characterisation Sc1 of the image of Sc in 'y-space'. The final step is to use Sivia again to go back to 'x-space' and find the reciprocal image Sc2 of Sc1!

Start by resetting A again, describe what we are doing, and get the precision

  // get a new AIASubPaving A to provide initial source box for next example
  delete A;

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

  A = new AIASPnode(x);

  // Image evaluation using set inversion
  // this uses the AIASubPaving Sc1 created by the above example
  // create an AIASubPaving Sc2 which contains inverse_f(Sc1)

  cout << "Characterization of the set Sc2=f-1(Sc1)" << endl
    << "with f^-1_1(x) = (2*x1+x2)/3," << endl
    << "     f^-1_2(x) = (x1+2*x2)/3," << endl;
  cout << "by realizing the inversion of f by Sivia" << endl;
  cout << "Enter a precision (between 1 and 0.01): ";

Now the argument supplied for the interval boolean test parameter is IBTFinverse from the top of the file. This specifies [f] for f1(x) = 2x1 - x2,f2(x) = -x1 + 2x2. IBTFinverse also specifies that the subpaving to invert is Sc1, ie the subpaving in 'y-space' we have just found using Sivia above.

  start = clock();
  Sc2 = Sivia(IBTFinverse,A,prec);
  end = clock();

Sc2 now points to the root node of a tree representing the subpaving constructed with Sivia.

The computing time, volume and number of leaves are reported and an output file produced.

  cout << "Computing time : " 
       << ((static_cast<double>(end - start)) / CLOCKS_PER_SEC) << " s."<< endl;
  cout << "Volume: " << Volume(Sc2) << endl;
  cout << "Number of leaves: " << NbLeaves(Sc2) << endl;

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

  cout << "The output AIASubPaving has been written to AIAinverse.txt" 
       << endl << endl;

  // end of testing SIVIA for set inversion
  // we should compare Sc2 to Sc in terms of volume and look at the effects of 
  // the precision variable

Finally we need to delete all the AIASPnode trees we created in dynamic memory, then the program can return.

  delete A;         // delete subpavings newed in dyamic memory
  delete Sc;
  delete Sc1;
  delete Sc2;

  return 0;
}

The example can be run with different values supplied for the precision in each case, to examine the effect that this has on the volume and number of leaves of the subpavings.

One of the interesting points about this example is that we can see the effect of pessimism in inclusion functions (AIA2001, pp. 15-17, p. 342). We start by making a subpaving Sc in 'x-space' and then find Sc1, a subpaving characterisation for the image of Sc under a function f (ie, a subpaving in 'y-space'), and then go back to 'x-space' again and with Sc2 as a subpaving characterisation for the reciprocal image of Sc1 under f.

Under the output format currently specified, the the first part of the output file AIAdirect.txt looks like this:

2
[ -5.000000,  5.000000] [ -5.000000,  5.000000]
Precision is 0.05
[  -1.953125 ,  -1.875000 ] , [  -0.078125 ,   0.000000 ]
[  -1.289062 ,  -1.250000 ] , [  -0.820312 ,  -0.781250 ]
[  -1.406250 ,  -1.328125 ] , [  -0.703125 ,  -0.625000 ]
[  -1.328125 ,  -1.289062 ] , [  -0.742188 ,  -0.703125 ]

All the output from this example rendered graphically looks like this.

AIAexample11_33.png
Results for Exercise 11.33 using precision 0.05

Exercise 11.35

See examples/AIA/Exr_11_35

Example 3.3

See AIA2001, p. 61.

Example 3.3 explores the use of ImageSp and the AIASPnode class for image evaluation. In Example 3.3 used Sivia to find the image of a subpaving under a function f, but this is only possible when f is invertible in the usual sense.

Computation of the direct image of a subpaving Y = f(X) where X is a regular subpaving of Rn (image evaluation) is more complex in the case where the function is not invertible in the usual sense.

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

Exm_3_3.hpp shows is similar to the header for Example 11.33 so we will move straight on the describing the file Exm_3_3.cpp.

First we include the header file and declare more global subpavings.

#include "Exm_3_3.hpp"

// These AIASubPavings are declared as global
AIASubPaving Sc, Sc4;

Now we look at the functions used in our main program.

The first is a boolean interval test IBTAnnular is exactly the same as that used for Exercise 11.33.

// specifications of example interval boolean tests
// The boolean interval test can return BI_TRUE, BI_FALSE, or BI_INDET

AIA_BOOL_INTERVAL IBTAnnular(const ivector& x)
{
  // here we test a 2-d box for inclusion in the area between circles centred 
  // on the origin with radii 1 and 2
  interval ToInvert(1.0,2.0),Temp;
  interval Img = sqr(x[1]) + sqr(x[2]);

  if (!Intersection(Temp,Img,ToInvert)) return BI_FALSE;
  if ( Img<=ToInvert ) return BI_TRUE;

  return BI_INDET;
}

// end specification of example boolean interval tests

The next function is an interval vector function, IVFex3_3. 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: R2 -> R2, f1(x1, x2) = x1x2, f2(x1, x2) = x1 + x2

// specification of example interval vector function,
// ie, vector inclusion function
// as in Example 3.3 of AIA
ivector IVF_ex3_3(const ivector& x)
{
  // example in 2-d space R2
  // for f: R2 -> R2
  //f1(x1, x2) = x1*x2
  //f2(x1, x2) = x1 + x2

  ivector Img(2);

  Img[1] = x[1] * x[2];
  Img[2] = x[1] + x[2];

  return (Img);
}

// end specification of example interval vector funciions

Then we say what namespaces we are using, start the main program, and set up an initial search box [-3.0 3.0]2 and a pointer A to an AIASPnode based on this search box.

using namespace cxsc;
using namespace std;

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);

  AIASubPaving A;
  A = new AIASPnode(x);

The first subpaving we create is again Sc, the subpaving covering the set in R2 between circles centred at the origin and with radii 1 and 2.

  // Using SIVIA for set inversion
  // find an AIASubPaving characterisation Sc containing the area between 
  // circles centred on the origin
  // with radii 1 and 2 (in 2 dimensional space)
  cout << "Characterization of the set Sc={(x1,x2) | 1 <= x1^2+x2^2 <= 2 }" 
       << endl;
  cout << "Enter a precision (between 1 and 0.001): ";
  cin >> prec;

  start = clock();

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

  Sc = Sivia(IBTAnnular,A,prec);
  end = clock();

  cout << "Computing time : " 
       << ((static_cast<double>(end - start)) / CLOCKS_PER_SEC) << " s."<< endl;
  cout << "Volume: " << Volume(Sc) << endl;
  cout << "Number of leaves: " << NbLeaves(Sc) << endl;

  // To realize a file output of the AIASubPaving Sc
                    // Filename
  ofstream os("AIA3_3a.txt");
  os << 2 << endl;  // Dimension of the AIASubPaving
                    // Root box
  os << interval(-5.0,5.0) << " "
    << interval(-5.0,5.0) << " " << endl;
                    // Precision used
  os << "Precision is " << prec << endl;
  os << Sc << endl; // AIASubPaving itself
  cout << "The output AIASubPaving has been written to AIA3_3a.txt" 
       << endl << endl;

  // the AIASubPaving Sc that we have created is the regular 
  // AIASubPaving that covers the set
  // X = {(x1,x2) in R2 | sqr(x1) + sqr(x2) is in [1,2]} 
  // (remember that it contains this area rather than being this area, 
  // and that eps has determined how small we go in the AIASubPaving 

Now we move on to using ImageSp. We are going to find Sc4, a characterisation of the set f(Sc) using f as defined in IVF_ex3_3.

First we say what we are doing and ask for a precision.

  cout << "Characterization of the set Sc4=f(Sc)" << endl
    << " with Sc from our first example and "<< endl
    << "with f1(x) = x1*x2," << endl
    << "     f2(x) = x1 + x2," << endl;
  cout << "by realizing the image of f by ImageSp" << endl;
  cout << "Enter a precision (between 1 and 0.01): ";

Then start the clock and run ImageSp(AIA_PIVF, AIASubPaving A, double eps) giving as arguments our interval vector function, our subpaving node A, and our precision.

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

The ImageSp algorithm works by taking an initial subpaving (in this case, the box of A), mincing it up into a fine subpaving where every box has width less than precision, 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.

  start = clock();

  // use Image SP to find a characterisation of the image of Sc using the 
  // function in IVF
  Sc4 = ImageSp(IVF_ex3_3, Sc, prec);

  end = clock();

We report the running time, volume and number of leaves in the image subpaving and send the output to a txt file

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

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

  cout << "The output AIASubPaving has been written to AIA3_3d.txt" 

Finally, we delete the subpavings and end the program.

  delete A;         // Delete all subpavings newed in dynamic memory
  delete Sc;
  delete Sc4;

  return 0;
}

The subpavings produced by this program run using precision 0.05 for both subpavings is shown graphically below. As well as showing the initial subpaving (represented by Sc) and the image subpaving (represented by Sc4), we capture an intermediate step within 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]2 chopped up so that each one is less than 0.05 wide). This set can be compared to the final regular minimal subpaving characterisation of the image which is formed by the function Regularize.

AIAexample3_3.png
Results for Example 3.3 using precision 0.05

Example 3.4

See AIA2001, pp. 61-63.

Finally we show another example using both Sivia and ImageSp. This example shows why we need ImageSp when to evaluate the images of functions that are only invertible in a set-theoretic sense, and also illustrates the interesting effects that can occur in these circumstances.

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

Exm_3_4.hpp shows is similar to the other header files so we will move straight on the describing the file Exm_3_4.cpp.

First we include the header file and declare more global subpavings.

#include "Exm_3_4.hpp"

// These AIASubPavings are declared as global
AIASubPaving Sc5, Sc6, Sc7;

Now we look at the functions used in our main program.

The first is a boolean interval test, IBT_ex3_4. The first specifies a function f: R2 -> R such that f(x1, x2) = x14 - x12 + 4x22 and tests the image of a box x under this function for inclusion in the subpaving to be inverted. As with IBTAnnular above, the subpaving to be inverted is also specifed in the test and here is just the interval [-0.1, 0.1].

// specifications of example interval boolean tests
// The boolean interval test can return BI_TRUE, BI_FALSE, or BI_INDET

AIA_BOOL_INTERVAL IBT_ex3_4(const ivector& x)
{
  // here we test a 2-d box for inclusion in the area 
  // such that x1^4 - x1^2 + 4x2^2 is in the interval [-0.1, 0.1]

  interval ToInvert(-0.1, 0.1),Temp;
  interval Img = power(x[1],4) - sqr(x[1]) + 4*sqr(x[2]);

  if (!Intersection(Temp,Img,ToInvert)) return BI_FALSE;
  if ( Img<=ToInvert ) return BI_TRUE;

  return BI_INDET;
}

Then we have another boolean interval test, IBTFinverse_ex3_4. Again this specifies a function f, this time f: R2 -> R2, f1(x1, x2) = (x1 - 1)2 -1 + x2, f2(x1, x2) = -x12 + (x2 - 1)2. The subpaving to invert is one of the globally defined subpaving Sc6

AIA_BOOL_INTERVAL IBTFinverse_ex3_4(const ivector& x)
{
  ivector Img(2);
  // 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

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

  return (Img<=Sc6);
}

// end specification of example boolean interval tests

The next function is an interval vector function, IVF_ex3_4. Recall that 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 say what namespaces we are using, start the main program, and set up an initial search box [-3.0 3.0]2 and a pointer A to an AIASPnode based on this search box.

using namespace cxsc;
using namespace std;

//int main(int argc, char* argv[])
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);

  AIASubPaving A;
  A = new AIASPnode(x);

First we use SIVIA to get a subpaving Sc5 in 'x-space' using the function f and subpaving to invert specified in IBT_ex3_4, ie Sc5 represents a subpaving characterisation of the set (x1, x2) such that x14 - x12 + 4x22 is in the interval [-0.1, 0.1].

  // Using SIVIA for set inversion

  // Find an AIASubPaving 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): ";
  cin >> prec;

  start = clock();

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

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

  cout << "Computing time : " 
       << ((static_cast<double>(end - start)) / CLOCKS_PER_SEC) << " s."<< endl;
  cout << "Volume: " << Volume(Sc5) << endl;
  cout << "Number of leaves: " << NbLeaves(Sc5) << endl;

And send the output to a txt file.

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

Now we want to find the image of this set using the function specified in IVF_ex3_4. f: R2 -> R2, f1(x1, x2) = (x1 - 1)2 -1 + x2, f2(x1, x2) = -x12 + (x2 - 1)2.

  // 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 certainly not invertible in the usual sense and we have to use ImageSp to find the image. 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
  Sc6 = ImageSp(IVF_ex3_4, Sc5, prec);

  end = clock();

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 " << Volume(Sc6) << endl;
  cout << "The number of leaves is " << NbLeaves(Sc6) << endl;

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

  cout << "The output AIASubPaving has been written to AIA3_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 AIASPnode(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'?

  // set inversion using SIVIA
  // this uses the AIASubPaving Sc6 created by the above example
  // create an AIASubPaving 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 and also specifies Sc6 as the subpaving to test for inclusion in. We give the new search box in the subpaving A, and the user-supplied precision.

  start = clock();
  Sc7 = Sivia(IBTFinverse_ex3_4,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: " << Volume(Sc7) << endl;
  cout << "Number of leaves: " << NbLeaves(Sc7) << endl;

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

  cout << "The output AIASubPaving has been written to AIA3_4c.txt" 

and then delete our subpavings and end the program

  delete A;         // delete all Subpavings newed in dynamic memory
  delete Sc5;
  delete Sc6;
  delete Sc7;

  return 0;
}

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

As well as showing the initial subpaving (represented by Sc5), the image subpaving (represented by Sc6), and the reciprocal image of the image (Sc7), 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]2 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 (Sc5) and the subpaving for a reciprocal image of the image subpaving (Sc7). 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)

AIAexample3_4.png
Results for Example 3.4 using precision 0.05
 All Classes Namespaces Functions Variables Typedefs Enumerations Friends