\input texinfo   @c -*-texinfo-*-
@c %**start of header
@setfilename ranlip.info
@settitle ranlip Manual 1.0
@c %**end of header

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

@copying
		    GNU GENERAL PUBLIC LICENSE
		       Version 2, June 1991

 Copyright @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.
@end copying


@c =============================================================================
@c the title page witch only appears if info page is passed to latex to produce
@c printed document.

@titlepage
@title Library ranlip for Multivariate Non-uniform Random Variate Generation
@vskip 0pt plus 1filll
@insertcopying
@end titlepage

@c==============================================================================
@c Output the table of the contents at the beginning.

@contents

@c =============================================================================
@c declaring the top node of the info page.

@ifnottex
@node Top
@top 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 @strong{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 @strong{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.

@insertcopying
@end ifnottex

@c =============================================================================
@c contents of the top node note that every entry in the menu must have a node
@c defined for it.

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

@menu
* Introduction::     				Overview of software.
* Description of Library:: 			Library interface methods. 
* Examples::					Library in action.  
* Performance::					Computational complexity and performance of the algorithms.
* Index::            				Complete index.
@end menu


@c ==============================================================================
@c first node definition must introduce chapter and index tags if you wish the
@c node to appear in the contents and index. Index tag can be placed anywhere
@c in the document and will serve as a link for quick referencing from the info
@c page index to the location of the tag.

@c @node Introduction, Features of the Interpolant, Top, Top
@node Introduction
@chapter Introduction

@cindex chapter, Introduction

This manual describes the programming library @strong{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 @i{rho} is known, or can be approximated, and that computation of the values of @i{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 @i{rho}. Lipschitz properties of @i{rho}allow one to overestimate @i{rho} at all other points, and thus to overestimate the absolute maxima of @i{rho} on the elements of the partition.

The class library @strong{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 @i{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 @i{rho} to be given analytically, to be differentiable, or to be normalized. 


@c===============================================================================
@c third node

@c @node Description of Library, Examples, Features of the Interpolant, Top
@node Description of Library
@chapter Description of Library
@cindex chapter, ranlip Description

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

@menu
* Class CRanLip::		interface of class CRanLip, inheriting from CRanLip
* Closer Look at Interface::	more detailed view of the interface.
@end menu


@c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@c chapter three section

@c @node Class STCInterpolant, Description of Library, Closer Look at Interface, Description of Library
@node Class CRanLip
@section Class CRanLip
@cindex section, interface of 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.

@example
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();
... 
@}
@end example


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

@example
// 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;
@end example


@c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@c chapter three section

@c @node Closer Look at Interface, Class STCInterpolant, Closer Look at Interface, Description of Library
@node Closer Look at Interface
@section Closer Look at Interface
@cindex section, RanLip library Interface

@subsection Init(int dim, double* left, double* right)

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

@subsection SetUniformGenerator(UFunction gen)

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

@subsection Seed(int seed)

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

@subsection PrepareHatFunction(int num, int numfine, double Lip)

Builds the hat function, using the Lipschitz constant supplied in @i{Lip}. Parameters @i{num} and @i{numfine} determine the rough and fine partitions. @i{num} is the number of subdivisions in each variable to partition the domain @i{D} into hyper rectangles @i{D_k}. On each @i{D_k}, the hat function will have a constant value @i{h_k}. numfine(>1) is the number of  subdivisions in the finer partition in each variable. Each set @i{D_k} is subdivided into @i{(numfine-1)^dim} smaller hyper rectangles, in order to improve the quality of the overestimate @i{h_k}. There will be in total @i{(num*numfine)^dim} evaluations of @i{rho} (calls to @i{Distribution()}). @i{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). @i{numfine} can be 2, in which case the fine partition is not used.

@subsection PrepareHatFunctionAuto(int num, int numfine, double minLip=0)

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

@subsection RandomVec(double* p)

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

@subsection FreeMemory()

Frees the memory occupied by the data structures, which can be very large. @i{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 @i{CRanLip} still exists. 


@c ==============================================================================
@c third node

@node Examples
@chapter Exemples

@cindex chapter, exemples

There are several examples of the usage of @strong{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.

@menu
* Common Use Example::	Common usage.
* Another Example::	Common usage.
* Procedural Example::	Uses C syntax.
@end menu


@c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@c third node subnode
@node Common Use Example
@section Common Use Example
@cindex Examples, Common Use Example

@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);
   @}
@}
@end example

@c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@c third node subnode

@c @node TNT Example, Examples, Common Use Example, Examples
@node Another Example
@section Another Example
@cindex Examples, Another Example

@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
   @}
@}
@end example

@c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@c third node subnode

@c @node TNT Example, Examples, Common Use Example, Examples
@node Procedural Example
@section Procedural Example
@cindex Examples, Procedural Example

@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);
@}
@end example

@c===============================================================================
@c second node
@node Performance
@chapter Performance

@cindex chapter, Computational complexity and performance

It is important to estimate the computing time and memory requirements when using @strong{ranlip}, especially for the case of several variables. The quality of
the hat function directly depends on the number of values of @i{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 @i{PrepareHatFunction(num, numfine, LipConst)} uses the first two parameters to establish the rough partition of the domain, sets @i{D_k}, on which  the hat function will have a constant value @i{h_k}, and the fine partition, used to compute this value. In total, @i{(num*numfine)^dim} evaluations of @i{rho} will be performed. The memory required by the algorihm is @i{numfine^dim + 2*num^dim} values of type @i{double} (8 bytes each).

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



@c ==============================================================================
@c here the index is declared as a node
@node Index
@unnumbered Index

@printindex cp

@bye