rand.cpp

Go to the documentation of this file.

//----------------------------------------------------------------------
//      File:                   rand.cpp
//      Programmer:             Sunil Arya and David Mount
//      Description:    Routines for random point generation
//      Last modified:  08/04/06 (Version 1.1.1)
//----------------------------------------------------------------------
// Copyright (c) 1997-2005 University of Maryland and Sunil Arya and
// David Mount.  All Rights Reserved.
// 
// This software and related documentation is part of the Approximate
// Nearest Neighbor Library (ANN).  This software is provided under
// the provisions of the Lesser GNU Public License (LGPL).  See the
// file ../ReadMe.txt for further information.
// 
// The University of Maryland (U.M.) and the authors make no
// representations about the suitability or fitness of this software for
// any purpose.  It is provided "as is" without express or implied
// warranty.
//----------------------------------------------------------------------
// History:
//      Revision 0.1  03/04/98
//              Initial release
//      Revision 0.2  03/26/98
//              Changed random/srandom declarations for SGI's.
//      Revision 1.0  04/01/05
//              annClusGauss centers distributed over [-1,1] rather than [0,1]
//              Added annClusOrthFlats distribution
//              Changed procedure names to avoid namespace conflicts
//              Added annClusFlats distribution
//              Added rand/srand option and fixed annRan0() initialization.
//      Revision 1.1.1  08/04/06
//              Added planted distribution
//----------------------------------------------------------------------

#include "rand.h"                                               // random generator declarations

using namespace std;                                    // make std:: accessible

//----------------------------------------------------------------------
//      Globals
//----------------------------------------------------------------------
int             annIdum = 0;                                    // used for random number generation

//------------------------------------------------------------------------
//      annRan0 - (safer) uniform random number generator
//
//      The code given here is taken from "Numerical Recipes in C" by
//      William Press, Brian Flannery, Saul Teukolsky, and William
//      Vetterling. The task of the code is to do an additional randomizing
//      shuffle on the system-supplied random number generator to make it
//      safer to use. 
//
//      Returns a uniform deviate between 0.0 and 1.0 using the
//      system-supplied routine random() or rand(). Set the global
//      annIdum to any negative value to initialise or reinitialise
//      the sequence.
//------------------------------------------------------------------------

double annRan0()
{
        const int TAB_SIZE = 97;                        // table size: any large number
        int j;

        static double y, v[TAB_SIZE];
        static int iff = 0;
        const double RAN_DIVISOR = double(ANN_RAND_MAX + 1UL);
        if (RAN_DIVISOR < 0) {
                cout << "RAN_DIVISOR " << RAN_DIVISOR << endl;
                exit(0);
        }

        //--------------------------------------------------------------------
        // As a precaution against misuse, we will always initialize on the
        // first call, even if "annIdum" is not set negative. Determine
        // "maxran", the next integer after the largest representable value
        // of type int. We assume this is a factor of 2 smaller than the
        // corresponding value of type unsigned int. 
        //--------------------------------------------------------------------

        if (annIdum < 0 || iff == 0) {          // initialize
                iff = 1;
                ANN_SRAND(annIdum);                             // (re)seed the generator
                annIdum = 1;

                for (j = 0; j < TAB_SIZE; j++)  // exercise the system routine
                        ANN_RAND();                                     // (values intentionally ignored)

                for (j = 0; j < TAB_SIZE; j++)  // then save TAB_SIZE-1 values
                        v[j] = ANN_RAND();
                y = ANN_RAND();                                 // generate starting value
         }

        //--------------------------------------------------------------------
        // This is where we start if not initializing. Use the previously
        // saved random number y to get an index j between 1 and TAB_SIZE-1.
        // Then use the corresponding v[j] for both the next j and as the
        // output number.
        //--------------------------------------------------------------------

        j = int(TAB_SIZE * (y / RAN_DIVISOR));
        y = v[j];
        v[j] = ANN_RAND();                                      // refill the table entry
        return y / RAN_DIVISOR;
}

//------------------------------------------------------------------------
//      annRanInt - generate a random integer from {0,1,...,n-1}
//
//              If n == 0, then -1 is returned.
//------------------------------------------------------------------------

static int annRanInt(
        int                                     n)
{
        int r = (int) (annRan0()*n);
        if (r == n) r--;                                        // (in case annRan0() == 1 or n == 0)
        return r;
}

//------------------------------------------------------------------------
//      annRanUnif - generate a random uniform in [lo,hi]
//------------------------------------------------------------------------

static double annRanUnif(
        double                          lo,
        double                          hi)
{
        return annRan0()*(hi-lo) + lo;
}

//------------------------------------------------------------------------
//      annRanGauss - Gaussian random number generator
//              Returns a normally distributed deviate with zero mean and unit
//              variance, using annRan0() as the source of uniform deviates.
//------------------------------------------------------------------------

static double annRanGauss()
{
        static int iset=0;
        static double gset;

        if (iset == 0) {                                        // we don't have a deviate handy
                double v1, v2;
                double r = 2.0;
                while (r >= 1.0) {
                        //------------------------------------------------------------
                        // Pick two uniform numbers in the square extending from -1 to
                        // +1 in each direction, see if they are in the circle of radius
                        // 1.  If not, try again 
                        //------------------------------------------------------------
                        v1 = annRanUnif(-1, 1);
                        v2 = annRanUnif(-1, 1);
                        r = v1 * v1 + v2 * v2;
                }
                double fac = sqrt(-2.0 * log(r) / r);
                //-----------------------------------------------------------------
                // Now make the Box-Muller transformation to get two normal
                // deviates.  Return one and save the other for next time.
                //-----------------------------------------------------------------
                gset = v1 * fac;
                iset = 1;                                               // set flag
                return v2 * fac;
        }
        else {                                                          // we have an extra deviate handy
                iset = 0;                                               // so unset the flag
                return gset;                                    // and return it
        }
}

//------------------------------------------------------------------------
//      annRanLaplace - Laplacian random number generator
//              Returns a Laplacian distributed deviate with zero mean and
//              unit variance, using annRan0() as the source of uniform deviates. 
//
//                              prob(x) = b/2 * exp(-b * |x|).
//
//              b is chosen to be sqrt(2.0) so that the variance of the Laplacian
//              distribution [2/(b^2)] becomes 1. 
//------------------------------------------------------------------------

static double annRanLaplace() 
{
        const double b = 1.4142136;

        double laprand = -log(annRan0()) / b;
        double sign = annRan0();
        if (sign < 0.5) laprand = -laprand;
        return(laprand);
}

//----------------------------------------------------------------------
//      annUniformPts - Generate uniformly distributed points
//              A uniform distribution over [-1,1].
//----------------------------------------------------------------------

void annUniformPts(                             // uniform distribution
        ANNpointArray           pa,                             // point array (modified)
        int                                     n,                              // number of points
        int                                     dim)                    // dimension
{
        for (int i = 0; i < n; i++) {
                for (int d = 0; d < dim; d++) {
                        pa[i][d] = (ANNcoord) (annRanUnif(-1,1));
                }
        }
}

//----------------------------------------------------------------------
//      annGaussPts - Generate Gaussian distributed points
//              A Gaussian distribution with zero mean and the given standard
//              deviation.
//----------------------------------------------------------------------

void annGaussPts(                                               // Gaussian distribution
        ANNpointArray           pa,                             // point array (modified)
        int                                     n,                              // number of points
        int                                     dim,                    // dimension
        double                          std_dev)                // standard deviation
{
        for (int i = 0; i < n; i++) {
                for (int d = 0; d < dim; d++) {
                        pa[i][d] = (ANNcoord) (annRanGauss() * std_dev);
                }
        }
}

//----------------------------------------------------------------------
//      annLaplacePts - Generate Laplacian distributed points
//              Generates a Laplacian distribution (zero mean and unit variance).
//----------------------------------------------------------------------

void annLaplacePts(                             // Laplacian distribution
        ANNpointArray           pa,                             // point array (modified)
        int                                     n,                              // number of points
        int                                     dim)                    // dimension
{
        for (int i = 0; i < n; i++) {
                for (int d = 0; d < dim; d++) {
                        pa[i][d] = (ANNcoord) annRanLaplace();
                }
        }
}

//----------------------------------------------------------------------
//      annCoGaussPts - Generate correlated Gaussian distributed points
//              Generates a Gauss-Markov distribution of zero mean and unit
//              variance.
//----------------------------------------------------------------------

void annCoGaussPts(                             // correlated-Gaussian distribution
        ANNpointArray           pa,                             // point array (modified)
        int                                     n,                              // number of points
        int                                     dim,                    // dimension
        double                          correlation)    // correlation
{
        double std_dev_w = sqrt(1.0 - correlation * correlation);
        for (int i = 0; i < n; i++) {
                double previous = annRanGauss();
                pa[i][0] = (ANNcoord) previous;
                for (int d = 1; d < dim; d++) {
                        previous = correlation*previous + std_dev_w*annRanGauss();
                        pa[i][d] = (ANNcoord) previous;
                } 
        }
}

//----------------------------------------------------------------------
//      annCoLaplacePts - Generate correlated Laplacian distributed points
//              Generates a Laplacian-Markov distribution of zero mean and unit
//              variance.
//----------------------------------------------------------------------

void annCoLaplacePts(                   // correlated-Laplacian distribution
        ANNpointArray           pa,                             // point array (modified)
        int                                     n,                              // number of points
        int                                     dim,                    // dimension
        double                          correlation)    // correlation
{
        double wn;
        double corr_sq = correlation * correlation;

        for (int i = 0; i < n; i++) {
                double previous = annRanLaplace();
                pa[i][0] = (ANNcoord) previous;
                for (int d = 1; d < dim; d++) {
                        double temp = annRan0();
                        if (temp < corr_sq)
                                wn = 0.0;
                        else
                                wn = annRanLaplace();
                        previous = correlation * previous + wn;
                        pa[i][d] = (ANNcoord) previous;
                } 
        }
}

//----------------------------------------------------------------------
//      annClusGaussPts - Generate clusters of Gaussian distributed points
//              Cluster centers are uniformly distributed over [-1,1], and the
//              standard deviation within each cluster is fixed.
//
//              Note: Once cluster centers have been set, they are not changed,
//              unless new_clust = true.  This is so that subsequent calls generate
//              points from the same distribution.  It follows, of course, that any
//              attempt to change the dimension or number of clusters without
//              generating new clusters is asking for trouble.
//
//              Note: Cluster centers are not generated by a call to uniformPts().
//              Although this could be done, it has been omitted for
//              compatibility with annClusGaussPts() in the colored version,
//              rand_c.cc.
//----------------------------------------------------------------------

void annClusGaussPts(                   // clustered-Gaussian distribution
        ANNpointArray           pa,                             // point array (modified)
        int                                     n,                              // number of points
        int                                     dim,                    // dimension
        int                                     n_clus,                 // number of colors
        ANNbool                         new_clust,              // generate new clusters.
        double                          std_dev)                // standard deviation within clusters
{
        static ANNpointArray clusters = NULL;// cluster storage

        if (clusters == NULL || new_clust) {// need new cluster centers
                if (clusters != NULL)                   // clusters already exist
                        annDeallocPts(clusters);        // get rid of them
                clusters = annAllocPts(n_clus, dim);
                                                                                // generate cluster center coords
                for (int i = 0; i < n_clus; i++) {
                        for (int d = 0; d < dim; d++) {
                                clusters[i][d] = (ANNcoord) annRanUnif(-1,1);
                        }
                }
        }

        for (int i = 0; i < n; i++) {
                int c = annRanInt(n_clus);                              // generate cluster index
                for (int d = 0; d < dim; d++) {
                  pa[i][d] = (ANNcoord) (std_dev*annRanGauss() + clusters[c][d]);
                }
        }
}

//----------------------------------------------------------------------
//      annClusOrthFlats - points clustered along orthogonal flats
//
//              This distribution consists of a collection points clustered
//              among a collection of axis-aligned low dimensional flats in
//              the hypercube [-1,1]^d.  A set of n_clus orthogonal flats are
//              generated, each whose dimension is a random number between 1
//              and max_dim.  The points are evenly distributed among the clusters.
//              For each cluster, we generate points uniformly distributed along
//              the flat within the hypercube.
//
//              This is done as follows.  Each cluster is defined by a d-element
//              control vector whose components are either:
//              
//                              CO_FLAG indicating that this component is to be generated
//                                              uniformly in [-1,1],
//                              x               a value other than CO_FLAG in the range [-1,1],
//                                              which indicates that this coordinate is to be
//                                              generated as x plus a Gaussian random deviation
//                                              with the given standard deviation.
//                                              
//              The number of zero components is the dimension of the flat, which
//              is a random integer in the range from 1 to max_dim.  The points
//              are disributed between clusters in nearly equal sized groups.
//
//              Note: Once cluster centers have been set, they are not changed,
//              unless new_clust = true.  This is so that subsequent calls generate
//              points from the same distribution.  It follows, of course, that any
//              attempt to change the dimension or number of clusters without
//              generating new clusters is asking for trouble.
//
//              To make this a bad scenario at query time, query points should be
//              selected from a different distribution, e.g. uniform or Gaussian.
//
//              We use a little programming trick to generate groups of roughly
//              equal size.  If n is the total number of points, and n_clus is
//              the number of clusters, then the c-th cluster (0 <= c < n_clus)
//              is given floor((n+c)/n_clus) points.  It can be shown that this
//              will exactly consume all n points.
//
//              This procedure makes use of the utility procedure, genOrthFlat
//              which generates points in one orthogonal flat, according to
//              the given control vector.
//
//----------------------------------------------------------------------
const double CO_FLAG = 999;                                             // special flag value

static void genOrthFlat(                // generate points on an orthog flat
        ANNpointArray           pa,                             // point array
        int                                     n,                              // number of points
        int                                     dim,                    // dimension
        double                          *control,               // control vector
        double                          std_dev)                // standard deviation
{
        for (int i = 0; i < n; i++) {                           // generate each point
                for (int d = 0; d < dim; d++) {                 // generate each coord
                        if (control[d] == CO_FLAG)                      // dimension on flat
                                pa[i][d] = (ANNcoord) annRanUnif(-1,1);
                        else                                                            // dimension off flat
                                pa[i][d] = (ANNcoord) (std_dev*annRanGauss() + control[d]);
                }
        }
}

void annClusOrthFlats(                  // clustered along orthogonal flats
        ANNpointArray           pa,                             // point array (modified)
        int                                     n,                              // number of points
        int                                     dim,                    // dimension
        int                                     n_clus,                 // number of colors
        ANNbool                         new_clust,              // generate new clusters.
        double                          std_dev,                // standard deviation within clusters
        int                                     max_dim)                // maximum dimension of the flats
{
        static ANNpointArray control = NULL;            // control vectors

        if (control == NULL || new_clust) {                     // need new cluster centers
                if (control != NULL) {                                  // clusters already exist
                        annDeallocPts(control);                         // get rid of them
                }
                control = annAllocPts(n_clus, dim);

                for (int c = 0; c < n_clus; c++) {              // generate clusters
                        int n_dim = 1 + annRanInt(max_dim); // number of dimensions in flat
                        for (int d = 0; d < dim; d++) {         // generate side locations
                                                                                                // prob. of picking next dim
                                double Prob = ((double) n_dim)/((double) (dim-d));
                                if (annRan0() < Prob) {                 // add this one to flat
                                        control[c][d] = CO_FLAG;        // flag this entry
                                        n_dim--;                                        // one fewer dim to fill
                                }
                                else {                                                  // don't take this one
                                        control[c][d] = annRanUnif(-1,1);// random value in [-1,1]
                                }
                        }
                }
        }
        int offset = 0;                                                         // offset in pa array
        for (int c = 0; c < n_clus; c++) {                      // generate clusters
                int pick = (n+c)/n_clus;                                // number of points to pick
                                                                                                // generate the points
                genOrthFlat(pa+offset, pick, dim, control[c], std_dev);
                offset += pick;                                                 // increment offset
        }
}

//----------------------------------------------------------------------
//      annClusEllipsoids - points clustered around axis-aligned ellipsoids
//
//              This distribution consists of a collection points clustered
//              among a collection of low dimensional ellipsoids whose axes
//              are alligned with the coordinate axes in the hypercube [-1,1]^d.
//              The objective is to model distributions in which the points are
//              distributed in lower dimensional subspaces, and within this
//              lower dimensional space the points are distributed with a
//              Gaussian distribution (with no correlation between the
//              dimensions).
//
//              The distribution is given the number of clusters or "colors"
//              (n_clus), maximum number of dimensions (max_dim) of the lower
//              dimensional subspace, a "small" standard deviation
//              (std_dev_small), and a "large" standard deviation range
//              (std_dev_lo, std_dev_hi).
//
//              The algorithm generates n_clus cluster centers uniformly from
//              the hypercube [-1,1]^d.  For each cluster, it selects the
//              dimension of the subspace as a random number r between 1 and
//              max_dim.  These are the dimensions of the ellipsoid.  Then it
//              generates a d-element std dev vector whose entries are the
//              standard deviation for the coordinates of each cluster in the
//              distribution.  Among the d-element control vector, r randomly
//              chosen values are chosen uniformly from the range [std_dev_lo,
//              std_dev_hi].  The remaining values are set to std_dev_small.
//
//              Note that annClusGaussPts is a special case of this in which
//              max_dim = 0, and std_dev = std_dev_small.
//
//              If the flag new_clust is set, then new cluster centers are
//              generated.
//
//              This procedure makes use of the utility procedure genGauss
//              which generates points distributed according to a Gaussian
//              distribution.
//
//----------------------------------------------------------------------

static void genGauss(                   // generate points on a general Gaussian
        ANNpointArray           pa,                             // point array
        int                                     n,                              // number of points
        int                                     dim,                    // dimension
        double                          *center,                // center vector
        double                          *std_dev)               // standard deviation vector
{
        for (int i = 0; i < n; i++) {
                for (int d = 0; d < dim; d++) {
                        pa[i][d] = (ANNcoord) (std_dev[d]*annRanGauss() + center[d]);
                }
        }
}

void annClusEllipsoids(                 // clustered around ellipsoids
        ANNpointArray           pa,                             // point array (modified)
        int                                     n,                              // number of points
        int                                     dim,                    // dimension
        int                                     n_clus,                 // number of colors
        ANNbool                         new_clust,              // generate new clusters.
        double                          std_dev_small,  // small standard deviation
        double                          std_dev_lo,             // low standard deviation for ellipses
        double                          std_dev_hi,             // high standard deviation for ellipses
        int                                     max_dim)                // maximum dimension of the flats
{
        static ANNpointArray centers = NULL;            // cluster centers
        static ANNpointArray std_dev = NULL;            // standard deviations

        if (centers == NULL || new_clust) {                     // need new cluster centers
                if (centers != NULL)                                    // clusters already exist
                        annDeallocPts(centers);                         // get rid of them
                if (std_dev != NULL)                                    // std deviations already exist
                        annDeallocPts(std_dev);                         // get rid of them

                centers = annAllocPts(n_clus, dim);             // alloc new clusters and devs
                std_dev  = annAllocPts(n_clus, dim);

                for (int i = 0; i < n_clus; i++) {              // gen cluster center coords
                        for (int d = 0; d < dim; d++) {
                                centers[i][d] = (ANNcoord) annRanUnif(-1,1);
                        }
                }
                for (int c = 0; c < n_clus; c++) {              // generate cluster std dev
                        int n_dim = 1 + annRanInt(max_dim); // number of dimensions in flat
                        for (int d = 0; d < dim; d++) {         // generate std dev's
                                                                                                // prob. of picking next dim
                                double Prob = ((double) n_dim)/((double) (dim-d));
                                if (annRan0() < Prob) {                 // add this one to ellipse
                                                                                                // generate random std dev
                                        std_dev[c][d] = annRanUnif(std_dev_lo, std_dev_hi);
                                        n_dim--;                                        // one fewer dim to fill
                                }
                                else {                                                  // don't take this one
                                        std_dev[c][d] = std_dev_small;// use small std dev
                                }
                        }
                }
        }

        int offset = 0;                                                         // next slot to fill
        for (int c = 0; c < n_clus; c++) {                      // generate clusters
                int pick = (n+c)/n_clus;                                // number of points to pick
                                                                                                // generate the points
                genGauss(pa+offset, pick, dim, centers[c], std_dev[c]);
                offset += pick;                                                 // increment offset in array
        }
}

//----------------------------------------------------------------------
//      annPlanted - Generates points from a "planted" distribution
//              In high dimensional spaces, interpoint distances tend to be
//              highly clustered around the mean value.  Approximate nearest
//              neighbor searching makes little sense in this context, unless it
//              is the case that each query point is significantly closer to its
//              nearest neighbor than to other points.  Thus, the query points
//              should be planted close to the data points. Given a source data
//              set, this procedure generates a set of query points having this
//              property.
//
//              We are given a source data array and a standard deviation.  We
//              generate points as follows.  We select a random point from the
//              source data set, and we generate a Gaussian point centered about
//              this random point and perturbed by a normal distributed random
//              variable with mean zero and the given standard deviation along
//              each coordinate.
//
//              Note that this essentially the same a clustered Gaussian
//              distribution, but where the cluster centers are given by the
//              source data set.
//----------------------------------------------------------------------

void annPlanted(                                // planted nearest neighbors
        ANNpointArray   pa,                     // point array (modified)
        int                             n,                      // number of points
        int                             dim,            // dimension
        ANNpointArray   src,            // source point array
        int                             n_src,          // source size
        double                  std_dev)        // standard deviation about source
{
        for (int i = 0; i < n; i++) {
                int c = annRanInt(n_src);                       // generate source index
                for (int d = 0; d < dim; d++) {
                  pa[i][d] = (ANNcoord) (std_dev*annRanGauss() + src[c][d]);
                }
        }
}