//----------------------------------------------------------------------
// 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]);
}
}
}