   =insert

Class Library lip for Multivariate Scattered Data Interpolation
***************************************************************

   This info manual describes a method of multivariate scattered data
interpolation and the programming library *lip* which implements this
method.

   We assume that the data set is generated by a Lipschitz continuous
function f, with the Lipschitz constant M.This method will interpolate
the data with a Lipschitz function g, which is the best approximation
to f in the worst case scenario. It provides best uniform error bounds,
and thus translates into reliable learning of f. The class library lip
is written in C++ language. It provides an interface to various methods
of calculating g.The choice of the method depends on the dimension of
the data, the size of the data set and on how many evaluations of the
interpolant are required. The usage of *lip* and its numerical
performance are illustrated on several examples.

   GNU GENERAL PUBLIC LICENCE

   Version 2, June 1991

   Copyright (C) 1989, 1991 Free Software Foundation, Inc. 59 Temple
Place, Suite 330, Boston, MA  02111-1307  USA Everyone is permitted to
copy and distribute verbatim copies of this license document, but
changing it is not allowed.

   The menu below lists the major nodes, which give a brief background,
overview of the library and illustrative examples on how to use it.

Introduction
************

   This document describes the programming library *lip*, which
implements a method of reliable multivariate interpolation of scattered
data.The underlying interpolation method assumes that the interpolated
function is Lipschitz continuous, and relies on the knowledge (or
estimate) of its Lipschitz constant. The Lipschitz constant is easily
interpreted as the upper bound on the rate of change of the function,
but the function needs not be differentiable. The computed interpolant
has a distinctive feature of providing the best approximation to the
unknown function that generated the data in the worst case scenario.

   The Lipschitz interpolant possesses a number of desirable features,
such as continuous dependence on the data, preservation of Lipschitz
properties and of the range of the data, uniform approximation and best
error bounds. On the practical side, construction and evaluation of the
interpolant is computationally stable. There are two alternative ways
to compute the interpolant: fast and explicit. The fast method involves
a preprocessing step after which the speed of evaluation is
proportional to the logarithm of the size of the data set. For higher
dimensions, the preprocessing step becomes too expensive, and limited
RAM may prevent efficient storage and use of the data structures. In
this case there is the second alternative to use the explicit
evaluation method, which does not require preprocessing. The
availability of two alternative evaluation methods is a distinctive
feature of the presented approach.

   The described interpolation method is highly competitive to the
alternative approaches in terms of efficiency and accuracy, and works
irrespectively of the dimension or distribution of data points.

Features of the Interpolant
***************************

   The interpolant g(x) implemented in this library provides the best
estimate of the unknown function f in the worst case scenario, based on
the provided data and its Lipschitz constant M. An estimate of M should
be provided by the user, but can also be computed from the data. In
many cases an educated guess is enough.

   In addition, this method has several useful features.

   The interpolant g provides a local approximation scheme (i.e.,
values of g depend only on the nearest data points).

  1. Preservation of the range of the data: min_k{y^k} <= g(x) <=
     max_k{y^k}.

  2. g approximates f uniformly. The upper bound on the error of
     approximation is M max_x min_k d(x,x^k), i.e., proportional to the
     distance between the most remote x and its nearest neighbour in
     the data set. This upper bound provides a guarantee on the quality
     of approximation regardless the distribution of data points or
     which particular function f element of Lip(M) generated these data.

  3. For simplicial distance d(x,x^k) the interpolant is a piecewise
     continuous linear function.

  4. The Lipschitz constant of g with respect to d is M.

  5. The interpolant g depends continuously on the data.

  6. The interpolant g provides a local approximation scheme (i.e.,
     values of g depend only on the nearest data points).

Description of Library
**********************

   The main class which provides the interface to the preprocessing and
computation is called `STCInterpolant'. It is illustrated in the
following sections, together with amore extencive description of the
interface.  Can be viewed via the following menu:

Class STCInterpolant
====================

     class STCInterpolant { //saw-tooth cover interpolant
     public:
      // supplies the data set to this class.
      // All other methods are called after SetData.
           void SetData(int dim, int K, double* x, double* y, int test=0);
      // Determines an estimate of the Lipschitz constant from the data
           double DetermineLipschitz();
      // Constructs the interpolant for either fast or explicit evaluation
           void Construct();
           void ConstructExplicit();
      // Computes the value of the interpolant
           double Value(int dim, double* x);
           double ValueExplicit(int dim, double* x);
     // Sets the value of the Lipschitz constant (or its estimate)
           void SetConstants(double newconst);
     ...
     };

     class SLipInt { //simple Lipschitz interpolant
     public:
     //  Computes the interpolant value by explicit method in l-2,l-Inf norm
           double  ValueInf(int dim,int K,double* x,double* X,double* Y,double LC);
      double  Value2(int dim,int K, double* x, double* X,double* Y,double LC);
     // Estimates the Lipschitz constant in l-2 or l-Inf norm
           void    ComputeLipschitzInf(int dim, int K, double* X, double* Y);
           void    ComputeLipschitz2(int dim, int K, double* X, double* Y);
     ...
     };

STCInterpolant member methods
=============================

SetData(int dim, int K, double* x, double* y, int test=0)
---------------------------------------------------------

   This method supplies the data set to be interpolated. The data set
is supplied in the variables x and y. x is a two-dimensional array K x
dim that contains values x^k_i in its rows, and y is an array of size K
which contains the values y^k. K is the size of the data set and $dim$
is the dimension. The last parameter test may be omitted, but if it is
set to 1, then the data will be screened for repeated x (which slows
down the process marginally). The construction algorithm assumes that
all data points are unique. If the user is unsure, test must be set to
1. The original data in x, y is not needed after the call to SetData.
SetData must be called only once before any other method.

double DetermineLipschitz()
---------------------------

   This method determines an (under)estimate of the Lipschitz constant
based on the data. Bear in mind that DetermineLipschitz is a rather
slow procedure O(mK^2) and should be avoided if there are other means
to estimate the Lipschitz constant. Returns the computed value of the
Lipschitz constant.

SetConstants(double newconst) and SetConstants()
------------------------------------------------

   Supplies the Lipschitz constant to the algorithm, either provided by
the user, or calculated from the data. It can be called with no
parameters after DetermineLipschitz only.

Construct()
-----------

   Performs the preprocessing step, necessary for the subsequent
evaluation by the fast method. This method may require significant
computational effort. *Should not be called after ConstructExplicit().*

ConstructExplicit()
-------------------

   Prepares the data for the subsequent evaluation of the interpolant
using explicit evaluation. No significant preprocessing. *Should not be
called after Construct().*

double Value(int dim, double* x )
---------------------------------

   Computes and returns the value of the interpolant g(x) using the
fast method. Should be called after preprocessing by Construct. If
called after ConstructExplicit, the explicit evaluation routine
ValueExplicit will be executed instead. Parameter x is an array of size
dim. dim should be specified as the dimension of the data m (the
algorithm will automatically compute the required slack variable). The
user may precompute the slack variable himself (see examples in the
next section), and specify dim as m+1, which is marginally faster.

double ValueExplicit(int dim, double* x )
-----------------------------------------

   Computes the value of the interpolant at x using explicit method.
Can be called after Construct or ConstructExplicit. As in Value, dim
should be specified as the dimension of the data m (the algorithm will
automatically compute the required slack variable). The user may
precompute the slack variable himself, and specify dim as m+1, which is
marginally faster. Parameter x is an array of size dim. Returns g(x).

FreeMemory()
------------

   Frees the memory occupied by the data structures computed in
Construct(), which can be very large. It destroys the interpolant, and
Value() methods cannot be called after FreeMemory(). Automatically
called from the destructor. This method is useful to deallocate memory
while the object STCInterpolant still exists.

int LastError()
---------------

   Returns the error code of the last operation. Useful to check
whether the construction or evaluation of the interpolant were
successful. 0 indicates a successful operation. Other possible values
defined in the interpol.h file are: 1 - lower interpolation failed, 2-
upper interpolation failed, 3 - both failed, 10 - Lipschitz constant
too small, or repeated data.

   *NOTE:* There are variations of Value and ValueExplicit methods,
which use parameter x of type dVec. This type is defined as
TNT::Vector<double> in common.h, and it uses TNT library (included).
TNT library available from (<http://math.nist.gov/tnt/>) provides a
convenient way of coding vectors and matrices, and is used internally
by *lip*. Since the user may prefer other ways to code vectors, generic
type double*  is used in this document. See interpol.h and examples for
more detailed information.

SLipInt member methods
======================

double  ValueInf(int dim,int K,double* x,double* X,double* Y,double LC)
-----------------------------------------------------------------------

   Computes and returns the value of the interpolant g(x) with
d(x,z)=max(|x_i -z_i|). Does not require any preprocessing. dim is the
dimension, K is the number of data, x is the vector of size dim, X is
the vector of data of size K * dim which contains values x^k_i$ in its
rows, y is the vector of size K of values to be interpolated, LC is the
Lipschitz constant.

double  ValueInf(int dim,int K,double* x,double* X,double* Y,double*LC)
-----------------------------------------------------------------------

   Variation of the above, but LC is a vector of size dim of Lipschitz
constants with respect to the i-th variable.

double  ValueInf(int dim,int K,double* x,double* X,double* Y)
-------------------------------------------------------------

   Variation of the above, uses Lipschitz constants automaticlly
identified from the data. *Should only be called after
ComputeLipschitzInf()*.

double  Value2(int dim,int K, double* x, double* X,double* Y,double LC)
-----------------------------------------------------------------------

   Computes and returns the value of the interpolant g(x) with d(x,z)
being Euclidean distance. The parameters are the same as in ValueInf().

double  Value2(int dim,int K, double* x, double* X,double* Y)
-------------------------------------------------------------

   Variation of the above, uses Lipschitz constants automaticlly
identified from the data. *Should only be called after
ComputeLipschitz2()*.

void  ComputeLipschitzInf(int dim,int K,double* X,double* Y)
------------------------------------------------------------

   Computes an (under)estimate of the Lipschitz constant in
$l-\infty$-norm from the data. Bear in mind that ComputeLipschitzInf is
a rather slow procedure O(nK^2) and should be avoided if there are
other means to estimate the Lipschitz constant. The value(s) are kept
in LipConst member variable.

void  ComputeLipschitz2(int dim,int K,double* X,double* Y)
----------------------------------------------------------

   Same as above, but uses Euclidean norm.

   There are versions of the above methods with Col suffix in their
names, these are used when the calling routine is written in another
language, like Fortran, which stores the data columnwise (not row wise,
like C). The only difference is that the matrix X is transposed.

Examples
********

   There are several examples of the usage of *lip* provided in the
distribution. The best way to use *lip* is to declare an instance of
the class STCInterpolant and call its members directly. There are four
basic steps: to supply the data, to supply the Lipschitz constant, to
construct the interpolant, and to evaluate it at the desired point.

Common Use Example
==================

     #include "interpol.h"
     int dim=4;             // the dimension and size of the data set
     int	npts=1000;

     void main(int argc, char *argv[]){
        double LipConst;
        double *x, *XData, *YData;
     // arrays to store the data
        x=(double*)malloc(dim*sizeof(double));
        XData=(double*)malloc(dim*npts*sizeof(double));
        YData=(double*)malloc(npts*sizeof(double));

        STCInterpolant LipInt;
        for(i=0;i<npts;i++) {
          for(j=0;j<dim;j++) // generate random data in [0,3]^m
            XData[i*dim + j]=x[j]=random(3.0,0);
          YData[i]=fun(x); // some function values
         }
     // suppy the data and eliminate repeated values
        LipInt.SetData(dim,npts, XData,YData,1);
        LipConst=LipInt.DetermineLipschitz();
        LipInt.SetConstants();  // supply Lipschitz constant
        LipInt.Construct();     // construct the interpolant
        free(XData); free(YData); // may now destroy the data

        double w,s,x1[10];  // reserve space for at least dim+1 components
        for(j=0;j<dim;j++) x1[j]=random(3.0,0); // some random x
        w=LipInt.Value(dim,x1); // calculate the value

     // alternatively, pre-compute the slack variable
        for(s=0,j=0; j<dim; j++) s+=x1[j];
        x1[dim]= 1.0-s;
        w=LipInt.Value(dim+1,x1);         // calculate the value
        w=LipInt.ValueExplicit(dim+1,x1); // same using explicit method
        LipInt.FreeMemory();  // destroys the interpolant
     }

TNT Example
===========

     // another example of usage with TNT library
     #include "interpol.h"
     #include "common.h"    // various types, such as dVec, dMat
     int dim=4;             // the dimension and the data set
     int	npts=1000;
     double LipConst=10;

     void main(int argc, char *argv[]){
     // arrays to store the data
        dVec  x(dim), YData(npts); // these types defined in common.h
        dMat  XData(npts,dim);
        STCInterpolant LipInt;
        for(i=0;i<npts;i++) {
          for(j=0;j<dim;j++) {     // generate random data in [0,3]^dim
            XData[i][j]=x[j]=random(3.0,0);
          }
          YData[i]=fun(x.begin()); // some function values
        }
     // suppy the data
        LipInt.SetData(dim,npts,&(XData[0][0]),YData.begin());
        LipInt.SetConstants(LipConst);  // supply Lipschitz constant
        LipInt.Construct();             // construct the interpolant

        double s,w;
        dVec x1(dim+1);
        for(j=0;j<dim;j++) x[j]=random(3.0,0);  // some random x
        w=LipInt.ValueSlack(x);         // calculate the value of g(x)
        w=LipInt.Value(dim,x.begin());  // the same calculation

     // alternatively, pre-compute the slack variable
        for(s=0,j=0; j<dim; j++) { x1[j]=random(3.0,0); s+=x1[j];}
        x1[dim]= 1.0-s;
        w=LipInt.Value(x1);         // calculate the value
        w=LipInt.ValueExplicit(x1); // same using explicit method
     }

Example Using SLipInt
=====================

     // example of usage of SLipInt class
     #include "slipint.h"
     int dim=4;        // the dimension and size of the data set
     int npts=1000;

     void main(int argc, char *argv[]){
        double LipConst=4;
        double *x, *XData, *YData;
     // arrays to store the data
        x=(double*)malloc(dim*sizeof(double));
        XData=(double*)malloc(dim*npts*sizeof(double));
        YData=(double*)malloc(npts*sizeof(double));

        SlipInt LipInt;
        for(i=0;i<npts;i++) {
          for(j=0;j<dim;j++) // generate random data in [0,3]^m
            XData[i*dim + j]=x[j]=random(3.0,0);
          YData[i]=fun(x); // some function values
         }

        for(j=0;j<dim;j++) x[j]=random(3.0,0); // some random x

     // calculate the value
        w=LipInt.ValueInf(dim,npts,x,XData, YData,LipConst);
     // estimate Lipschitz constant
        LipInt.ComputeLipschitzInf(dim,npts,XData, YData);
     // uses the computed Lipschitz constant
        w=LipInt.ValueInf(dim,npts,x,XData, YData);

     // the same in Euclidean norm
        LipInt.ComputeLipschitz2(dim,npts,XData, YData);
     // calculate the value
        w=LipInt.ValueInf(dim,npts,x,XData, YData);

        // Overwrite Lipschitz constant and calculate value
        w=LipInt.ValueInf(dim,npts,x,XData, YData,LipConst);

         free(XData); free(YData); free(x);
     }