============================================================================== = this is the copyright label which can be duplicated later on using the insert copyright command

==============================================================================

Class Library ranlip for Multivariate Non-uniform Random Variate Generation
***************************************************************************

This manual describes generation of non-uniform random variates from Lipschitz-continuous densities using acceptance/ rejection, and the class library ranlip which implements this method. It is assumed that the required distribution has Lipschitz-continuous density, which is either given analytically or as a black box. The algorithm builds a piecewise constant upper approximation to the density (the hat function), using a large number of its values and subdivision of the domain into hyper rectangles.

The class library ranlip provides very competitive preprocessing and generation times, and yields  small rejection constant, which is a measure of efficiency of the generation step. It exhibits good performance for up to five variables, and provides the user with a black box non-uniform random variate generator for a large class of distributions, in particular, multimodal distributions.

GNU GENERAL PUBLIC LICENSE Version 2, June 1991

Copyright © 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 sections which give a brief background, overview of the library and illustrative examples on how to use it.

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

This manual describes the programming library ranlip, which implements the method of acceptance/ rejection in the multivariate case, for Lipschitz continuous  densities. It assumes that the Lipschitz constant of the density rho is known, or can be approximated, and that computation of the values of rho at distinct points is not expensive. The method builds a piecewise constant hat function, by subdividing the domain into hyper rectangles, and by using a large number of values of rho. Lipschitz properties of rhoallow one to overestimate rho at all other points, and thus to overestimate the absolute maxima of rho on the elements of the partition.

The class library ranlip implements computation of the hat function and generation of random variates, and makes this process transparent to the user. The user needs to provide a method of evaluation of rho at a given point, and the number of elements in the subdivision of the domain, which is the  parameter characterizing the quality of the hat function and the number of computations at the preprocessing step.

The class of Lipschitz-continuous densities is very broad, and includes many multimodal densities, which are hard to deal with. No other properties beyond Lipschitz continuity are required, and the Lipschitz constant, if not provided, can be estimated automatically. The algorithm does not require rho to be given analytically, to be differentiable, or to be normalized.

===============================================================================

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

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

Class CRanLip
=============

The main class which provides the interface to the preprocessing and random variate generation is called CRanLip. This is an abstract class, from which the user must derive his own class which overrides the method for rho(x), and declare an instance of that class. The interface of Class CRandLip is illustrated bellow followed by an example of how to inherit from CRanLip.

     class CRanLip  {
     public:
     // initialises the tables and arrays for the generator
     // should be the first method to call
         void Init(int dim, double* left, double* right);

     // sets the pointer to the uniform random number generator.
     // The default is ranlux generator
         void SetUniformGenerator(UFunction gen);

     // sets the seed for the uniform random number generator
         void Seed(int seed);

     // computes the hat function given the Lipschitz constant
     // num and numfine determine the rough and fine partitions
         void PrepareHatFunction(int num, int numfine, double Lip);

     // computes the hat function, and automatically computes
     // the Lipschitz constant
         void PrepareHatFunctionAuto(int num, int numfine, double minLip);

     // generates a random variate with the required density
         void RandomVec(double* p);

     // frees the memory occupied by the partition and various tables
          void FreeMem();
     ...
     }

The example bellow illustrates how the user needs to inherit from CRanLip.

     // declare a derived class MyRnumGen, with one method to override
     class MyRnumGen:public CRanLip {
     	 public: virtual double	Distribution(double* p) ;
     };

     double	MyRnumGen::Distribution(double* p)
     { // example: multivriate normal distribution
        double r;
        for(int j=0;j<Dimension;j++) {
           r+=p[j]*p[j];
        }
        return exp(-r);;
     }
     // declare an instance of this class
     MyRnumGen MyGen;

Closer Look at Interface
========================

Init(int dim, double* left, double* right)
------------------------------------------

Initialises the internal variables of this class. dim is the dimension, left and right are arrays of size dim which determine the domain of rho: left_i <= x_i <= right_i. Init must be called only once before any other method.

SetUniformGenerator(UFunction gen)
----------------------------------

Sets a pointer to the uniform random number generator on (0,1). The default is M. Luescher's ranlux generator, but the user can override it and use his own preferred generator.

Seed(int seed)
--------------

Sets the seed of the default uniform random number generator ranlux. If the user has supplied his own generator in SetUniformGenerator, that generator's seed() function should be called instead.

PrepareHatFunction(int num, int numfine, double Lip)
----------------------------------------------------

Builds the hat function, using the Lipschitz constant supplied in Lip. Parameters num and numfine determine the rough and fine partitions. num is the number of subdivisions in each variable to partition the domain D into hyper rectangles D_k. On each D_k, the hat function will have a constant value h_k. numfine(>1) is the number of  subdivisions in the finer partition in each variable. Each set D_k is subdivided into (numfine-1)^dim smaller hyper rectangles, in order to improve the quality of the overestimate h_k. There will be in total (num*numfine)^dim evaluations of rho (calls to Distribution()). numfine should be a power of 2 for numerical efficiency reasons (if not, it will be automatically changed to a power of 2 larger than the supplied value). numfine can be 2, in which case the fine partition is not used.

PrepareHatFunctionAuto(int num, int numfine, double minLip=0)
-------------------------------------------------------------

Builds the hat function and automatically computes an estimate to the Lipschitz constant.  Parameters num and numfine determine the rough and fine partitions, and are described in PrepareHatFunction() method. minLip denotes the lower bound on the value of the computed Lipschitz constant, the default value is 0.

RandomVec(double* p)
--------------------

Generates a random variate with the density rho. Should be called after PrepareHatFunction() or PrepareHatFunctionAuto(). The parameter p is an array of size dim, it will contain the components of the computed random vector.

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

Frees the memory occupied by the data structures, which can be very large. It destroys the hat function, and RandomVec() method cannot be called after FreeMemory(). Automatically called from the destructor. This method is useful to deallocate memory while the object CRanLip still exists.

Exemples
********

There are several examples of the usage of ranlip provided in the distribution. There are three basic steps: to supply the required distribution, to build the hat function, and to generate random variates. As illustrated in following examples.

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

     #include "ranlip.h"
     #define dim 4             // the dimension
     // declare a derived class MyRnumGen, with one method to override
     class MyRnumGen:public CRanLip {
     	 public: virtual double	Distribution(double* p) ;
     };
     double	MyRnumGen::Distribution(double* p)
     { // example: multivriate normal distribution
        double r=0.0;
        for(int j=0;j<Dimension;j++)  r+=p[j]*p[j];
        return exp(-r);
     }
     void main(int argc, char *argv[]){
        double LipConst = 4.0;
        MyRnumGen MyGen;
        double left[dim], right[dim], p[dim];
        int i;
     // set the domain to be the unit hypercube
        for(i=0;i<dim;i++) {left[i]=0; right[i]=1;}

        MyGen.Init(dim,left,right);
        MyGen.PrepareHatFunction(10,8,LipConst);
        MyGen.Seed(10);
        for(i=0;i<1000;i++) {
           MyGen.RandomVec(p);
        // do something with p
        }
        MyGen.FreeMemory();
     // now MyGen can be reused
        MyGen.Init(dim,left,right);
        MyGen.PrepareHatFunctionAuto(10,32);
        for(i=0;i<1000;i++) {
           MyGen.RandomVec(p);
        }
     }

Another Example
===============

     // Example 2
     #include "ranlip.h"
     #define dim 3             // the dimension
     // declare a derived class MyRnumGen, with one method to override
     class MyRnumGen:public CRanLip {
     	 public: virtual double	Distribution(double* p) ;
     };
     double	MyRnumGen::Distribution(double* p)
     { // example: multivriate normal distribution
        double r=0.0;
        for(int j=0;j<Dimension;j++)  r+=p[j]*p[j];
        return exp(-r);
     }
     void main(int argc, char *argv[]){
        double LipConst;
        MyRnumGen MyGen;
        double left[dim], right[dim], p[dim];
        int i;
     // set the domain to be a hypercube
        for(i=0;i<dim;i++) {left[i]=-2; right[i]=2;}
        MyGen.Init(dim,left,right);
        MyGen.PrepareHatFunctionAuto(10,32,0.01);
        MyGen.Seed(10);

        for(i=0;i<1000;i++) {
           MyGen.RandomVec(p);
        // do something with p
        }
        cout<<"acceptance ratio is "<<1000.0/ MyGen.count_total<<endl;
        LipConst=MyGen.Lipschitz; // computed Lipschitz constant
        if(MyGen.count_error>0) { // Lipschitz constant was too low
           LipConst*=2;
           MyGen.PrepareHatFunction(10,32,LipConst);
        }
        for(i=0;i<1000;i++) {
           MyGen.RandomVec(p);
        // do something with p
        }
     }

Procedural Example
==================

     // Example 3  using procedural interface
     #include "ranlipproc.h"
     #define Dim 3             // the dimension
     // implement calculation of the density in MyDist function
     double MyDist(double* p, int dim)
     { // example: multivariate normal distribution
        double r=0.0;
        for(int j=0;j<dim;j++)  r+=p[j]*p[j];
        return exp(-r);
     }
     double MyRand() // use my own random number generator
     { return  (double)rand()/(RAND_MAX+0.001); } //random numbers in [0,1)

     void main(int argc, char *argv[]){
        double LipConst; int i;
        double left[Dim], right[Dim], p[Dim];
     // set the domain to be a hypercube
        for(i=0;i<Dim;i++) {left[i]=-2; right[i]=2;}

        InitRanLip(Dim,a,b);
     // pass the address of the density function (required)
        SetDistFunctionRanLip(&MyDist);
     // pass the address of my generator (optional)
        SetUniformGeneratorRanLip(&MyRand);
        srand(10); // use srand, not SeedRanLip

        PrepareHatFunctionAutoRanLip(10,8);
        for(i=0;i<1000;i++) {
           RandomVecRanLip(p);
        // do something with p
        }
        cout<<"acceptance ratio is "<<1000.0/ Count_totalRanLip() <<endl;
        LipConst=LipschitzRanLip(); // computed Lipschitz constant
        if(Count_errorRanLip()>0) { // Lipschitz constant was too low
           LipConst*=2;
           PrepareHatFunctionRanLip(10,32,LipConst);
        }
        for(i=0;i<1000;i++)
           RandomVecRanLip(p);
     }

===============================================================================

Performance
***********

It is important to estimate the computing time and memory requirements when using ranlip, especially for the case of several variables. The quality of the hat function directly depends on the number of values of rho(x) used for its computation. The higher this number is, the longer is preprocessing step (building the hat function), but the more efficient is the generation step (less rejected variates). A good estimate of the Lipschitz constant is also important, as it improves  the quality of the hat function.

The method PrepareHatFunction(num, numfine, LipConst) uses the first two parameters to establish the rough partition of the domain, sets D_k, on which  the hat function will have a constant value h_k, and the fine partition, used to compute this value. In total, (num*numfine)^dim evaluations of rho will be performed. The memory required by the algorihm is numfine^dim + 2*num^dim values of type double (8 bytes each).

For numerical efficiency, the second parameter numfine should be a power of 2. This facilitates computation of the neighbours  in an n-dimensional mesh by using binary arithmetic.