//----------------------------------------------------------------------
// File: ann_test.cpp
// Programmer: Sunil Arya and David Mount
// Description: test program for ANN (approximate nearest neighbors)
// 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 06/26/98
// Added CLOCKS_PER_SEC definition if needed
// Revision 1.0 04/01/05
// Added comments (from "#" to eol)
// Added clus_orth_flats and clus_ellipsoids distributions
// Fixed order of fair and midpt in split_table
// Added dump/load operations
// Cleaned up C++ for modern compilers
// Revision 1.1 05/03/05
// Added fixed radius kNN search
// Revision 1.1.1 08/04/06
// Added planted distribution
//----------------------------------------------------------------------
#include <ctime> // clock
#include <cmath> // math routines
#include <string> // C string ops
#include <fstream> // file I/O
#include <cstring>
#include <ANN/ANN.h> // ANN declarations
#include <ANN/ANNx.h> // more ANN declarations
#include <ANN/ANNperf.h> // performance evaluation
#include "rand.h" // random point generation
#ifndef CLOCKS_PER_SEC // define clocks-per-second if needed
#define CLOCKS_PER_SEC 1000000
#endif
using namespace std; // make std:: available
//----------------------------------------------------------------------
// ann_test
//
// This program is a driver for testing and evaluating the ANN library
// for computing approximate nearest neighbors. It allows the user to
// generate data and query sets of various sizes, dimensions, and
// distributions, to build kd- and bbd-trees of various types, and then
// run queries and outputting various performance statistics.
//
// Overview:
// ---------
// The test program is run as follows:
//
// ann_test < test_input > test_output
//
// where the test_input file contains a list of directives as described
// below. Directives consist of a directive name, followed by list of
// arguments (depending on the directive). Arguments and directives are
// separated by white space (blank, tab, and newline). String arguments
// are not quoted, and consist of a string of nonwhite chacters. A
// character "#" denotes a comment. The following characters up to
// the end of line are ignored. Comments may only be inserted between
// directives (not within the argument list of a directive).
//
// Basic operations:
// -----------------
// The test program can perform the following operations. How these
// operations are performed depends on the options which are described
// later.
//
// Data Generation:
// ----------------
// read_data_pts <file> Create a set of data points whose
// coordinates are input from file <file>.
// gen_data_pts Create a set of data points whose
// coordinates are generated from the
// current point distribution.
//
// Building the tree:
// ------------------
// build_ann Generate an approximate nearest neighbor
// structure for the current data set, using
// the selected splitting rules. Any existing
// tree will be destroyed.
//
// Query Generation/Searching:
// ---------------------------
// read_query_pts <file> Create a set of query points whose
// coordinates are input from file <file>.
// gen_query_pts Create a set of query points whose
// coordinates are generated from the
// current point distribution.
// run_queries <string> Apply nearest neighbor searching to the
// query points using the approximate nearest
// neighbor structure and the given search
// strategy. Possible strategies are:
// standard = standard kd-tree search
// priority = priority search
//
// Miscellaneous:
// --------------
// output_label Output a label to the output file.
// dump <file> Dump the current structure to given file.
// (The dump format is explained further in
// the source file kd_tree.cc.)
// load <file> Load a tree from a data file which was
// created by the dump operation. Any
// existing tree will be destroyed.
//
// Options:
// --------
// How these operations are performed depends on a set of options.
// If an option is not specified, a default value is used. An option
// retains its value until it is set again. String inputs are not
// enclosed in quotes, and must contain no embedded white space (sorry,
// this is C++'s convention).
//
// Options affecting search tree structure:
// ----------------------------------------
// split_rule <type> Type of splitting rule to use in building
// the search tree. Choices are:
// kd = optimized kd-tree
// midpt = midpoint split
// fair = fair split
// sl_midpt = sliding midpt split
// sl_fair = sliding fair split
// suggest = authors' choice for best
// The default is "suggest". See the file
// kd_split.cc for more detailed information.
//
// shrink_rule <type> Type of shrinking rule to use in building
// a bd-tree data structure. If "none" is
// given, then no shrinking is performed and
// the result is a kd-tree. Choices are:
// none = perform no shrinking
// simple = simple shrinking
// centroid = centroid shrinking
// suggest = authors' choice for best
// The default is "none". See the file
// bd_tree.cc for more information.
// bucket_size <int> Bucket size, that is, the maximum number of
// points stored in each leaf node.
//
// Options affecting data and query point generation:
// --------------------------------------------------
// dim <int> Dimension of space.
// seed <int> Seed for random number generation.
// data_size <int> Number of data points. When reading data
// points from a file, this indicates the
// maximum number of points for storage
// allocation. Default = 100.
// query_size <int> Same as data_size for query points.
// std_dev <float> Standard deviation (used in gauss,
// planted, and clustered distributions).
// This is the "small" distribution for
// clus_ellipsoids. Default = 1.
// std_dev_lo <float> Low and high standard deviations (used in
// std_dev_hi <float> clus_ellipsoids). Default = 1.
// corr_coef <float> Correlation coefficient (used in co-gauss
// and co_lapace distributions). Default = 0.05.
// colors <int> Number of color classes (clusters) (used
// in the clustered distributions). Default = 5.
// new_clust Once generated, cluster centers are not
// normally regenerated. This is so that both
// query points and data points can be generated
// using the same set of clusters. This option
// forces new cluster centers to be generated
// with the next generation of either data or
// query points.
// max_clus_dim <int> Maximum dimension of clusters (used in
// clus_orth_flats and clus_ellipsoids).
// Default = 1.
// distribution <string> Type of input distribution
// uniform = uniform over cube [-1,1]^d.
// gauss = Gaussian with mean 0
// laplace = Laplacian, mean 0 and var 1
// co_gauss = correlated Gaussian
// co_laplace = correlated Laplacian
// clus_gauss = clustered Gaussian
// clus_orth_flats = clusters of orth flats
// clus_ellipsoids = clusters of ellipsoids
// planted = planted distribution
// See the file rand.cpp for further information.
//
// Options affecting nearest neighbor search:
// ------------------------------------------
// epsilon <float> Error bound for approx. near neigh. search.
// near_neigh <int> Number of nearest neighbors to compute.
// max_pts_visit <int> Maximum number of points to visit before
// terminating. (Used in applications where
// real-time performance is important.)
// (Default = 0, which means no limit.)
// radius_bound <float> Sets an upper bound on the nearest
// neighbor search radius. If the bound is
// positive, then fixed-radius nearest
// neighbor searching is performed, and the
// count of the number of points in the
// range is returned. If the bound is
// zero, then standard search is used.
// This can only be used with standard, not
// priority, search. (Default = 0, which
// means standard search.)
//
// Options affection general program behavior:
// -------------------------------------------
// stats <string> Level of statistics output
// silent = no output,
// exec_time += execution time only
// prep_stats += preprocessing statistics
// query_stats += query performance stats
// query_res += results of queries
// show_pts += show the data points
// show_struct += print search structure
// validate <string> Validate experiment and compute average
// error. Since validation causes exact
// nearest neighbors to be computed by the
// brute force method, this can take a long
// time. Valid arguments are:
// on = turn validation on
// off = turn validation off
// true_near_neigh <int> Number of true nearest neighbors to compute.
// When validating, we compute the difference
// in rank between each reported nearest neighbor
// and the true nearest neighbor of the same
// rank. Thus it is necessary to compute a
// few more true nearest neighbors. By default
// we compute 10 more than near_neigh. With
// this option the exact number can be set.
// (Used only when validating.)
//
// Example:
// --------
// output_label test_run_0 # output label for this run
// validate off # do not perform validation
// dim 16 # points in dimension 16
// stats query_stats # output performance statistics for queries
// seed 121212 # random number seed
// data_size 1000
// distribution uniform
// gen_data_pts # 1000 uniform data points in dim 16
// query_size 100
// std_dev 0.05
// distribution clus_gauss
// gen_query_pts # 100 points in 10 clusters with std_dev 0.05
// bucket_size 2
// split_rule kd
// shrink_rule none
// build_ann # kd-tree, bucket size 2
// epsilon 0.1
// near_neigh 5
// max_pts_visit 100 # stop search if more than 100 points seen
// run_queries standard # run queries; 5 nearest neighbors, 10% error
// data_size 500
// read_data_pts data.in # read up to 500 points from file data.in
// split_rule sl_midpt
// shrink_rule simple
// build_ann # bd-tree; simple shrink, sliding midpoint split
// epsilon 0
// run_queries priority # run same queries; 0 allowable error
//
//------------------------------------------------------------------------
//------------------------------------------------------------------------
// Constants
//------------------------------------------------------------------------
const int STRING_LEN = 500; // max string length
const double ERR = 0.00001; // epsilon (for float compares)
//------------------------------------------------------------------------
// Enumerated values and conversions
//------------------------------------------------------------------------
typedef enum {DATA, QUERY} PtType; // point types
//------------------------------------------------------------------------
// Statistics output levels
//------------------------------------------------------------------------
typedef enum { // stat levels
SILENT, // no output
EXEC_TIME, // just execution time
PREP_STATS, // preprocessing info
QUERY_STATS, // query performance
QUERY_RES, // query results
SHOW_PTS, // show data points
SHOW_STRUCT, // show tree structure
N_STAT_LEVELS} // number of levels
StatLev;
const char stat_table[N_STAT_LEVELS][STRING_LEN] = {
"silent", // SILENT
"exec_time", // EXEC_TIME
"prep_stats", // PREP_STATS
"query_stats", // QUERY_STATS
"query_res", // QUERY_RES
"show_pts", // SHOW_PTS
"show_struct"}; // SHOW_STRUCT
//------------------------------------------------------------------------
// Distributions
//------------------------------------------------------------------------
typedef enum { // distributions
UNIFORM, // uniform over cube [-1,1]^d.
GAUSS, // Gaussian with mean 0
LAPLACE, // Laplacian, mean 0 and var 1
CO_GAUSS, // correlated Gaussian
CO_LAPLACE, // correlated Laplacian
CLUS_GAUSS, // clustered Gaussian
CLUS_ORTH_FLATS, // clustered on orthog flats
CLUS_ELLIPSOIDS, // clustered on ellipsoids
PLANTED, // planted distribution
N_DISTRIBS}
Distrib;
const char distr_table[N_DISTRIBS][STRING_LEN] = {
"uniform", // UNIFORM
"gauss", // GAUSS
"laplace", // LAPLACE
"co_gauss", // CO_GAUSS
"co_laplace", // CO_LAPLACE
"clus_gauss", // CLUS_GAUSS
"clus_orth_flats", // CLUS_ORTH_FLATS
"clus_ellipsoids", // CLUS_ELLIPSOIS
"planted"}; // PLANTED
//------------------------------------------------------------------------
// Splitting rules for kd-trees (see ANN.h for types)
//------------------------------------------------------------------------
const int N_SPLIT_RULES = 6;
const char split_table[N_SPLIT_RULES][STRING_LEN] = {
"standard", // standard optimized kd-tree
"midpt", // midpoint split
"fair", // fair split
"sl_midpt", // sliding midpt split
"sl_fair", // sliding fair split
"suggest"}; // authors' choice for best
//------------------------------------------------------------------------
// Shrinking rules for bd-trees (see ANN.h for types)
//------------------------------------------------------------------------
const int N_SHRINK_RULES = 4;
const char shrink_table[N_SHRINK_RULES][STRING_LEN] = {
"none", // perform no shrinking (kd-tree)
"simple", // simple shrinking
"centroid", // centroid shrinking
"suggest"}; // authors' choice for best
//----------------------------------------------------------------------
// Short utility functions
// Error - general error routine
// printPoint - print a point to standard output
// lookUp - look up a name in table and return index
//----------------------------------------------------------------------
void Error( // error routine
char *msg, // error message
ANNerr level) // abort afterwards
{
if (level == ANNabort) {
cerr << "ann_test: ERROR------->" << msg << "<-------------ERROR\n";
exit(1);
}
else {
cerr << "ann_test: WARNING----->" << msg << "<-------------WARNING\n";
}
}
void printPoint( // print point
ANNpoint p, // the point
int dim) // the dimension
{
cout << "[";
for (int i = 0; i < dim; i++) {
cout << p[i];
if (i < dim-1) cout << ",";
}
cout << "]";
}
int lookUp( // look up name in table
const char *arg, // name to look up
const char (*table)[STRING_LEN], // name table
int size) // table size
{
int i;
for (i = 0; i < size; i++) {
if (!strcmp(arg, table[i])) return i;
}
return i;
}
//------------------------------------------------------------------------
// Function declarations
//------------------------------------------------------------------------
void generatePts( // generate data/query points
ANNpointArray &pa, // point array (returned)
int n, // number of points
PtType type, // point type
ANNbool new_clust, // new cluster centers desired?
ANNpointArray src = NULL, // source array (for PLANTED)
int n_src = 0); // source size (for PLANTED)
void readPts( // read data/query points from file
ANNpointArray &pa, // point array (returned)
int &n, // number of points
char *file_nm, // file name
PtType type); // point type (DATA, QUERY)
void doValidation(); // perform validation
void getTrueNN(); // compute true nearest neighbors
void treeStats( // print statistics on kd- or bd-tree
ostream &out, // output stream
ANNbool verbose); // print stats
//------------------------------------------------------------------------
// Default execution parameters
//------------------------------------------------------------------------
const int extra_nn = 10; // how many extra true nn's?
const int def_dim = 2; // def dimension
const int def_data_size = 100; // def data size
const int def_query_size = 100; // def number of queries
const int def_n_color = 5; // def number of colors
const ANNbool def_new_clust = ANNfalse; // def new clusters flag
const int def_max_dim = 1; // def max flat dimension
const Distrib def_distr = UNIFORM; // def distribution
const double def_std_dev = 1.00; // def standard deviation
const double def_corr_coef = 0.05; // def correlation coef
const int def_bucket_size = 1; // def bucket size
const double def_epsilon = 0.0; // def error bound
const int def_near_neigh = 1; // def number of near neighbors
const int def_max_visit = 0; // def number of points visited
const int def_rad_bound = 0; // def radius bound
// def number of true nn's
const int def_true_nn = def_near_neigh + extra_nn;
const int def_seed = 0; // def seed for random numbers
const ANNbool def_validate = ANNfalse; // def validation flag
// def statistics output level
const StatLev def_stats = QUERY_STATS;
const ANNsplitRule // def splitting rule
def_split = ANN_KD_SUGGEST;
const ANNshrinkRule // def shrinking rule
def_shrink = ANN_BD_NONE;
//------------------------------------------------------------------------
// Global variables - Execution options
//------------------------------------------------------------------------
int dim; // dimension
int data_size; // data size
int query_size; // number of queries
int n_color; // number of colors
ANNbool new_clust; // generate new clusters?
int max_dim; // maximum flat dimension
Distrib distr; // distribution
double corr_coef; // correlation coef
double std_dev; // standard deviation
double std_dev_lo; // low standard deviation
double std_dev_hi; // high standard deviation
int bucket_size; // bucket size
double epsilon; // error bound
int near_neigh; // number of near neighbors
int max_pts_visit; // max number of points to visit
double radius_bound; // maximum radius search bound
int true_nn; // number of true nn's
ANNbool validate; // validation flag
StatLev stats; // statistics output level
ANNsplitRule split; // splitting rule
ANNshrinkRule shrink; // shrinking rule
//------------------------------------------------------------------------
// More globals - pointers to dynamically allocated arrays and structures
//
// It is assumed that all these values are set to NULL when nothing
// is allocated.
//
// data_pts, query_pts The data and query points
// the_tree Points to the kd- or bd-tree for
// nearest neighbor searching.
// apx_nn_idx, apx_dists Record approximate near neighbor
// indices and distances
// apx_pts_in_range Counts of the number of points in
// the in approx range, for fixed-
// radius NN searching.
// true_nn_idx, true_dists Record true near neighbor
// indices and distances
// min_pts_in_range, max_... Min and max counts of the number
// of points in the in approximate
// range.
// valid_dirty To avoid repeated validation,
// we only validate query results
// once. This validation becomes
// invalid, if a new tree, new data
// points or new query points have
// been generated.
// tree_data_size The number of points in the
// current tree. (This will be the
// same a data_size unless points have
// been added since the tree was
// built.)
//
// The approximate and true nearest neighbor results are stored
// in: apx_nn_idx, apx_dists, and true_nn_idx, true_dists.
// They are really flattened 2-dimensional arrays. Each of these
// arrays consists of query_size blocks, each of which contains
// near_neigh (or true_nn) entries, one for each of the nearest
// neighbors for a given query point.
//------------------------------------------------------------------------
ANNpointArray data_pts; // data points
ANNpointArray query_pts; // query points
ANNbd_tree* the_tree; // kd- or bd-tree search structure
ANNidxArray apx_nn_idx; // storage for near neighbor indices
ANNdistArray apx_dists; // storage for near neighbor distances
int* apx_pts_in_range; // storage for no. of points in range
ANNidxArray true_nn_idx; // true near neighbor indices
ANNdistArray true_dists; // true near neighbor distances
int* min_pts_in_range; // min points in approx range
int* max_pts_in_range; // max points in approx range
ANNbool valid_dirty; // validation is no longer valid
//------------------------------------------------------------------------
// Initialize global parameters
//------------------------------------------------------------------------
void initGlobals()
{
dim = def_dim; // init execution parameters
data_size = def_data_size;
query_size = def_query_size;
distr = def_distr;
corr_coef = def_corr_coef;
std_dev = def_std_dev;
std_dev_lo = def_std_dev;
std_dev_hi = def_std_dev;
new_clust = def_new_clust;
max_dim = def_max_dim;
n_color = def_n_color;
bucket_size = def_bucket_size;
epsilon = def_epsilon;
near_neigh = def_near_neigh;
max_pts_visit = def_max_visit;
radius_bound = def_rad_bound;
true_nn = def_true_nn;
validate = def_validate;
stats = def_stats;
split = def_split;
shrink = def_shrink;
annIdum = -def_seed; // init. global seed for ran0()
data_pts = NULL; // initialize storage pointers
query_pts = NULL;
the_tree = NULL;
apx_nn_idx = NULL;
apx_dists = NULL;
apx_pts_in_range = NULL;
true_nn_idx = NULL;
true_dists = NULL;
min_pts_in_range = NULL;
max_pts_in_range = NULL;
valid_dirty = ANNtrue; // (validation must be done)
}
//------------------------------------------------------------------------
// getDirective - skip comments and read next directive
// Returns ANNtrue if directive read, and ANNfalse if eof seen.
//------------------------------------------------------------------------
ANNbool skipComment( // skip any comments
istream &in) // input stream
{
char ch = 0;
// skip whitespace
do { in.get(ch); } while (isspace(ch) && !in.eof());
while (ch == '#' && !in.eof()) { // comment?
// skip to end of line
do { in.get(ch); } while(ch != '\n' && !in.eof());
// skip whitespace
do { in.get(ch); } while(isspace(ch) && !in.eof());
}
if (in.eof()) return ANNfalse; // end of file
in.putback(ch); // put character back
return ANNtrue;
}
ANNbool getDirective(
istream &in, // input stream
char *directive) // directive storage
{
if (!skipComment(in)) // skip comments
return ANNfalse; // found eof along the way?
in >> directive; // read directive
return ANNtrue;
}
//------------------------------------------------------------------------
// main program - driver
// The main program reads input options, invokes the necessary
// routines to process them.
//------------------------------------------------------------------------
int main(int argc, char** argv)
{
long clock0; // clock time
char directive[STRING_LEN]; // input directive
char arg[STRING_LEN]; // all-purpose argument
cout << "------------------------------------------------------------\n"
<< "ann_test: Version " << ANNversion << " " << ANNversionCmt << "\n"
<< " Copyright: " << ANNcopyright << ".\n"
<< " Latest Revision: " << ANNlatestRev << ".\n"
<< "------------------------------------------------------------\n\n";
initGlobals(); // initialize global values
//--------------------------------------------------------------------
// Main input loop
//--------------------------------------------------------------------
// read input directive
while (getDirective(cin, directive)) {
//----------------------------------------------------------------
// Read options
//----------------------------------------------------------------
if (!strcmp(directive,"dim")) {
cin >> dim;
}
else if (!strcmp(directive,"colors")) {
cin >> n_color;
}
else if (!strcmp(directive,"new_clust")) {
new_clust = ANNtrue;
}
else if (!strcmp(directive,"max_clus_dim")) {
cin >> max_dim;
}
else if (!strcmp(directive,"std_dev")) {
cin >> std_dev;
}
else if (!strcmp(directive,"std_dev_lo")) {
cin >> std_dev_lo;
}
else if (!strcmp(directive,"std_dev_hi")) {
cin >> std_dev_hi;
}
else if (!strcmp(directive,"corr_coef")) {
cin >> corr_coef;
}
else if (!strcmp(directive, "data_size")) {
cin >> data_size;
}
else if (!strcmp(directive,"query_size")) {
cin >> query_size;
}
else if (!strcmp(directive,"bucket_size")) {
cin >> bucket_size;
}
else if (!strcmp(directive,"epsilon")) {
cin >> epsilon;
}
else if (!strcmp(directive,"max_pts_visit")) {
cin >> max_pts_visit;
valid_dirty = ANNtrue; // validation must be redone
}
else if (!strcmp(directive,"radius_bound")) {
cin >> radius_bound;
valid_dirty = ANNtrue; // validation must be redone
}
else if (!strcmp(directive,"near_neigh")) {
cin >> near_neigh;
true_nn = near_neigh + extra_nn; // also reset true near neighs
valid_dirty = ANNtrue; // validation must be redone
}
else if (!strcmp(directive,"true_near_neigh")) {
cin >> true_nn;
valid_dirty = ANNtrue; // validation must be redone
}
//----------------------------------------------------------------
// seed option
// The seed is reset by setting the global annIdum to the
// negation of the seed value. See rand.cpp.
//----------------------------------------------------------------
else if (!strcmp(directive,"seed")) {
cin >> annIdum;
annIdum = -annIdum;
}
//----------------------------------------------------------------
// validate option
//----------------------------------------------------------------
else if (!strcmp(directive,"validate")) {
cin >> arg; // input argument
if (!strcmp(arg, "on")) {
validate = ANNtrue;
cout << "validate = on "
<< "(Warning: this may slow execution time.)\n";
}
else if (!strcmp(arg, "off")) {
validate = ANNfalse;
}
else {
cerr << "Argument: " << arg << "\n";
Error("validate argument must be \"on\" or \"off\"", ANNabort);
}
}
//----------------------------------------------------------------
// distribution option
//----------------------------------------------------------------
else if (!strcmp(directive,"distribution")) {
cin >> arg; // input name and translate
distr = (Distrib) lookUp(arg, distr_table, N_DISTRIBS);
if (distr >= N_DISTRIBS) { // not something we recognize
cerr << "Distribution: " << arg << "\n";
Error("Unknown distribution", ANNabort);
}
}
//----------------------------------------------------------------
// stats option
//----------------------------------------------------------------
else if (!strcmp(directive,"stats")) {
cin >> arg; // input name and translate
stats = (StatLev) lookUp(arg, stat_table, N_STAT_LEVELS);
if (stats >= N_STAT_LEVELS) { // not something we recognize
cerr << "Stats level: " << arg << "\n";
Error("Unknown statistics level", ANNabort);
}
if (stats > SILENT)
cout << "stats = " << arg << "\n";
}
//----------------------------------------------------------------
// split_rule option
//----------------------------------------------------------------
else if (!strcmp(directive,"split_rule")) {
cin >> arg; // input split_rule name
split = (ANNsplitRule) lookUp(arg, split_table, N_SPLIT_RULES);
if (split >= N_SPLIT_RULES) { // not something we recognize
cerr << "Splitting rule: " << arg << "\n";
Error("Unknown splitting rule", ANNabort);
}
}
//----------------------------------------------------------------
// shrink_rule option
//----------------------------------------------------------------
else if (!strcmp(directive,"shrink_rule")) {
cin >> arg; // input split_rule name
shrink = (ANNshrinkRule) lookUp(arg, shrink_table, N_SHRINK_RULES);
if (shrink >= N_SHRINK_RULES) { // not something we recognize
cerr << "Shrinking rule: " << arg << "\n";
Error("Unknown shrinking rule", ANNabort);
}
}
//----------------------------------------------------------------
// label operation
//----------------------------------------------------------------
else if (!strcmp(directive,"output_label")) {
cin >> arg;
if (stats > SILENT)
cout << "<" << arg << ">\n";
}
//----------------------------------------------------------------
// gen_data_pts operation
//----------------------------------------------------------------
else if (!strcmp(directive,"gen_data_pts")) {
if (distr == PLANTED) { // planted distribution
Error("Cannot use planted distribution for data points", ANNabort);
}
generatePts( // generate data points
data_pts, // data points
data_size, // data size
DATA, // data points
new_clust); // new clusters flag
valid_dirty = ANNtrue; // validation must be redone
new_clust = ANNfalse; // reset flag
}
//----------------------------------------------------------------
// gen_query_pts operation
// If the distribution is PLANTED, then the query points
// are planted near the data points (which must already be
// generated).
//----------------------------------------------------------------
else if (!strcmp(directive,"gen_query_pts")) {
if (distr == PLANTED) { // planted distribution
if (data_pts == NULL) {
Error("Must generate data points before query points for planted distribution", ANNabort);
}
generatePts( // generate query points
query_pts, // point array
query_size, // number of query points
QUERY, // query points
new_clust, // new clusters flag
data_pts, // plant around data pts
data_size);
}
else { // all other distributions
generatePts( // generate query points
query_pts, // point array
query_size, // number of query points
QUERY, // query points
new_clust); // new clusters flag
}
valid_dirty = ANNtrue; // validation must be redone
new_clust = ANNfalse; // reset flag
}
//----------------------------------------------------------------
// read_data_pts operation
//----------------------------------------------------------------
else if (!strcmp(directive,"read_data_pts")) {
cin >> arg; // input file name
readPts(
data_pts, // point array
data_size, // number of points
arg, // file name
DATA); // data points
valid_dirty = ANNtrue; // validation must be redone
}
//----------------------------------------------------------------
// read_query_pts operation
//----------------------------------------------------------------
else if (!strcmp(directive,"read_query_pts")) {
cin >> arg; // input file name
readPts(
query_pts, // point array
query_size, // number of points
arg, // file name
QUERY); // query points
valid_dirty = ANNtrue; // validation must be redone
}
//----------------------------------------------------------------
// build_ann operation
// We always invoke the constructor for bd-trees. Note
// that when the shrinking rule is NONE (which is true by
// default), then this constructs a kd-tree.
//----------------------------------------------------------------
else if (!strcmp(directive,"build_ann")) {
//------------------------------------------------------------
// Build the tree
//------------------------------------------------------------
if (the_tree != NULL) { // tree exists already
delete the_tree; // get rid of it
}
clock0 = clock(); // start time
the_tree = new ANNbd_tree( // build it
data_pts, // the data points
data_size, // number of points
dim, // dimension of space
bucket_size, // maximum bucket size
split, // splitting rule
shrink); // shrinking rule
//------------------------------------------------------------
// Print summary
//------------------------------------------------------------
long prep_time = clock() - clock0; // end of prep time
if (stats > SILENT) {
cout << "[Build ann-structure:\n";
cout << " split_rule = " << split_table[split] << "\n";
cout << " shrink_rule = " << shrink_table[shrink] << "\n";
cout << " data_size = " << data_size << "\n";
cout << " dim = " << dim << "\n";
cout << " bucket_size = " << bucket_size << "\n";
if (stats >= EXEC_TIME) { // output processing time
cout << " process_time = "
<< double(prep_time)/CLOCKS_PER_SEC << " sec\n";
}
if (stats >= PREP_STATS) // output or check tree stats
treeStats(cout, ANNtrue); // print tree stats
else
treeStats(cout, ANNfalse); // check stats
if (stats >= SHOW_STRUCT) { // print the whole tree
cout << " (Structure Contents:\n";
the_tree->Print(ANNfalse, cout);
cout << " )\n";
}
cout << "]\n";
}
}
//----------------------------------------------------------------
// dump operation
//----------------------------------------------------------------
else if (!strcmp(directive,"dump")) {
cin >> arg; // input file name
if (the_tree == NULL) { // no tree
Error("Cannot dump. No tree has been built yet", ANNwarn);
}
else { // there is a tree
// try to open file
ofstream out_dump_file(arg);
if (!out_dump_file) {
cerr << "File name: " << arg << "\n";
Error("Cannot open dump file", ANNabort);
}
// dump the tree and points
the_tree->Dump(ANNtrue, out_dump_file);
if (stats > SILENT) {
cout << "(Tree has been dumped to file " << arg << ")\n";
}
}
}
//----------------------------------------------------------------
// load operation
// Since this not only loads a tree, but loads a new set
// of data points.
//----------------------------------------------------------------
else if (!strcmp(directive,"load")) {
cin >> arg; // input file name
if (the_tree != NULL) { // tree exists already
delete the_tree; // get rid of it
}
if (data_pts != NULL) { // data points exist already
delete data_pts; // get rid of them
}
ifstream in_dump_file(arg); // try to open file
if (!in_dump_file) {
cerr << "File name: " << arg << "\n";
Error("Cannot open file for loading", ANNabort);
}
// build tree by loading
the_tree = new ANNbd_tree(in_dump_file);
dim = the_tree->theDim(); // new dimension
data_size = the_tree->nPoints(); // number of points
data_pts = the_tree->thePoints(); // new points
valid_dirty = ANNtrue; // validation must be redone
if (stats > SILENT) {
cout << "(Tree has been loaded from file " << arg << ")\n";
}
if (stats >= SHOW_STRUCT) { // print the tree
cout << " (Structure Contents:\n";
the_tree->Print(ANNfalse, cout);
cout << " )\n";
}
}
//----------------------------------------------------------------
// run_queries operation
// This section does all the query processing. It consists
// of the following subsections:
//
// ** input the argument (standard or priority) and output
// the header describing the essential information.
// ** allocate space for the results to be stored.
// ** run the queries by invoking the appropriate search
// procedure on the query points. Print nearest neighbor
// if requested.
// ** print final summaries
//
// The approach for processing multiple nearest neighbors is
// pretty crude. We allocate an array whose size is the
// product of the total number of queries times the number of
// nearest neighbors (k), and then use each k consecutive
// entries to store the results of each query.
//----------------------------------------------------------------
else if (!strcmp(directive,"run_queries")) {
//------------------------------------------------------------
// Input arguments and print summary
//------------------------------------------------------------
enum {STANDARD, PRIORITY} method;
cin >> arg; // input argument
if (!strcmp(arg, "standard")) {
method = STANDARD;
}
else if (!strcmp(arg, "priority")) {
method = PRIORITY;
}
else {
cerr << "Search type: " << arg << "\n";
Error("Search type must be \"standard\" or \"priority\"",
ANNabort);
}
if (data_pts == NULL || query_pts == NULL) {
Error("Either data set and query set not constructed", ANNabort);
}
if (the_tree == NULL) {
Error("No search tree built.", ANNabort);
}
//------------------------------------------------------------
// Set up everything
//------------------------------------------------------------
#ifdef ANN_PERF // performance only
annResetStats(data_size); // reset statistics
#endif
clock0 = clock(); // start time
// deallocate existing storage
if (apx_nn_idx != NULL) delete [] apx_nn_idx;
if (apx_dists != NULL) delete [] apx_dists;
if (apx_pts_in_range != NULL) delete [] apx_pts_in_range;
// allocate apx answer storage
apx_nn_idx = new ANNidx[near_neigh*query_size];
apx_dists = new ANNdist[near_neigh*query_size];
apx_pts_in_range = new int[query_size];
annMaxPtsVisit(max_pts_visit); // set max points to visit
//------------------------------------------------------------
// Run the queries
//------------------------------------------------------------
// pointers for current query
ANNidxArray curr_nn_idx = apx_nn_idx;
ANNdistArray curr_dists = apx_dists;
for (int i = 0; i < query_size; i++) {
#ifdef ANN_PERF
annResetCounts(); // reset counters
#endif
apx_pts_in_range[i] = 0;
if (radius_bound == 0) { // no radius bound
if (method == STANDARD) {
the_tree->annkSearch(
query_pts[i], // query point
near_neigh, // number of near neighbors
curr_nn_idx, // nearest neighbors (returned)
curr_dists, // distance (returned)
epsilon); // error bound
}
else if (method == PRIORITY) {
the_tree->annkPriSearch(
query_pts[i], // query point
near_neigh, // number of near neighbors
curr_nn_idx, // nearest neighbors (returned)
curr_dists, // distance (returned)
epsilon); // error bound
}
else {
Error("Internal error - invalid method", ANNabort);
}
}
else { // use radius bound
if (method != STANDARD) {
Error("A nonzero radius bound assumes standard search",
ANNwarn);
}
apx_pts_in_range[i] = the_tree->annkFRSearch(
query_pts[i], // query point
ANN_POW(radius_bound), // squared radius search bound
near_neigh, // number of near neighbors
curr_nn_idx, // nearest neighbors (returned)
curr_dists, // distance (returned)
epsilon); // error bound
}
curr_nn_idx += near_neigh; // increment current pointers
curr_dists += near_neigh;
#ifdef ANN_PERF
annUpdateStats(); // update stats
#endif
}
long query_time = clock() - clock0; // end of query time
if (validate) { // validation requested
if (valid_dirty) getTrueNN(); // get true near neighbors
doValidation(); // validate
}
//------------------------------------------------------------
// Print summaries
//------------------------------------------------------------
if (stats > SILENT) {
cout << "[Run Queries:\n";
cout << " query_size = " << query_size << "\n";
cout << " dim = " << dim << "\n";
cout << " search_method = " << arg << "\n";
cout << " epsilon = " << epsilon << "\n";
cout << " near_neigh = " << near_neigh << "\n";
if (max_pts_visit != 0)
cout << " max_pts_visit = " << max_pts_visit << "\n";
if (radius_bound != 0)
cout << " radius_bound = " << radius_bound << "\n";
if (validate)
cout << " true_nn = " << true_nn << "\n";
if (stats >= EXEC_TIME) { // print exec time summary
cout << " query_time = " <<
double(query_time)/(query_size*CLOCKS_PER_SEC)
<< " sec/query";
#ifdef ANN_PERF
cout << " (biased by perf measurements)";
#endif
cout << "\n";
}
if (stats >= QUERY_STATS) { // output performance stats
#ifdef ANN_PERF
cout.flush();
annPrintStats(validate);
#else
cout << " (Performance statistics unavailable.)\n";
#endif
}
if (stats >= QUERY_RES) { // output results
cout << " (Query Results:\n";
cout << " Pt\tANN\tDist\n";
curr_nn_idx = apx_nn_idx; // subarray pointers
curr_dists = apx_dists;
// output nearest neighbors
for (int i = 0; i < query_size; i++) {
cout << " " << setw(4) << i;
for (int j = 0; j < near_neigh; j++) {
// exit if no more neighbors
if (curr_nn_idx[j] == ANN_NULL_IDX) {
cout << "\t[no other pts in radius bound]\n";
break;
}
else { // output point info
cout << "\t" << curr_nn_idx[j]
<< "\t" << ANN_ROOT(curr_dists[j])
<< "\n";
}
}
// output range count
if (radius_bound != 0) {
cout << " pts_in_radius_bound = "
<< apx_pts_in_range[i] << "\n";
}
// increment subarray pointers
curr_nn_idx += near_neigh;
curr_dists += near_neigh;
}
cout << " )\n";
}
cout << "]\n";
}
}
//----------------------------------------------------------------
// Unknown directive
//----------------------------------------------------------------
else {
cerr << "Directive: " << directive << "\n";
Error("Unknown directive", ANNabort);
}
}
//--------------------------------------------------------------------
// End of input loop (deallocate stuff that was allocated)
//--------------------------------------------------------------------
if (the_tree != NULL) delete the_tree;
if (data_pts != NULL) annDeallocPts(data_pts);
if (query_pts != NULL) annDeallocPts(query_pts);
if (apx_nn_idx != NULL) delete [] apx_nn_idx;
if (apx_dists != NULL) delete [] apx_dists;
if (apx_pts_in_range != NULL) delete [] apx_pts_in_range;
annClose(); // close ANN
return EXIT_SUCCESS;
}
//------------------------------------------------------------------------
// generatePts - call appropriate routine to generate points of a
// given distribution.
//------------------------------------------------------------------------
void generatePts(
ANNpointArray &pa, // point array (returned)
int n, // number of points to generate
PtType type, // point type
ANNbool new_clust, // new cluster centers desired?
ANNpointArray src, // source array (if distr=PLANTED)
int n_src) // source size (if distr=PLANTED)
{
if (pa != NULL) annDeallocPts(pa); // get rid of any old points
pa = annAllocPts(n, dim); // allocate point storage
switch (distr) {
case UNIFORM: // uniform over cube [-1,1]^d.
annUniformPts(pa, n, dim);
break;
case GAUSS: // Gaussian with mean 0
annGaussPts(pa, n, dim, std_dev);
break;
case LAPLACE: // Laplacian, mean 0 and var 1
annLaplacePts(pa, n, dim);
break;
case CO_GAUSS: // correlated Gaussian
annCoGaussPts(pa, n, dim, corr_coef);
break;
case CO_LAPLACE: // correlated Laplacian
annCoLaplacePts(pa, n, dim, corr_coef);
break;
case CLUS_GAUSS: // clustered Gaussian
annClusGaussPts(pa, n, dim, n_color, new_clust, std_dev);
break;
case CLUS_ORTH_FLATS: // clustered on orthog flats
annClusOrthFlats(pa, n, dim, n_color, new_clust, std_dev, max_dim);
break;
case CLUS_ELLIPSOIDS: // clustered ellipsoids
annClusEllipsoids(pa, n, dim, n_color, new_clust, std_dev,
std_dev_lo, std_dev_hi, max_dim);
break;
case PLANTED: // planted distribution
annPlanted(pa, n, dim, src, n_src, std_dev);
break;
default:
Error("INTERNAL ERROR: Unknown distribution", ANNabort);
break;
}
if (stats > SILENT) {
if(type == DATA) cout << "[Generating Data Points:\n";
else cout << "[Generating Query Points:\n";
cout << " number = " << n << "\n";
cout << " dim = " << dim << "\n";
cout << " distribution = " << distr_table[distr] << "\n";
if (annIdum < 0)
cout << " seed = " << annIdum << "\n";
if (distr == GAUSS || distr == CLUS_GAUSS
|| distr == CLUS_ORTH_FLATS)
cout << " std_dev = " << std_dev << "\n";
if (distr == CLUS_ELLIPSOIDS) {
cout << " std_dev = " << std_dev << " (small) \n";
cout << " std_dev_lo = " << std_dev_lo << "\n";
cout << " std_dev_hi = " << std_dev_hi << "\n";
}
if (distr == CO_GAUSS || distr == CO_LAPLACE)
cout << " corr_coef = " << corr_coef << "\n";
if (distr == CLUS_GAUSS || distr == CLUS_ORTH_FLATS
|| distr == CLUS_ELLIPSOIDS) {
cout << " colors = " << n_color << "\n";
if (new_clust)
cout << " (cluster centers regenerated)\n";
}
if (distr == CLUS_ORTH_FLATS || distr == CLUS_ELLIPSOIDS) {
cout << " max_dim = " << max_dim << "\n";
}
}
// want to see points?
if ((type == DATA && stats >= SHOW_PTS) ||
(type == QUERY && stats >= QUERY_RES)) {
if(type == DATA) cout << "(Data Points:\n";
else cout << "(Query Points:\n";
for (int i = 0; i < n; i++) {
cout << " " << setw(4) << i << "\t";
printPoint(pa[i], dim);
cout << "\n";
}
cout << " )\n";
}
cout << "]\n";
}
//------------------------------------------------------------------------
// readPts - read a collection of data or query points.
//------------------------------------------------------------------------
void readPts(
ANNpointArray &pa, // point array (returned)
int &n, // number of points
char *file_nm, // file name
PtType type) // point type (DATA, QUERY)
{
int i;
//--------------------------------------------------------------------
// Open input file and read points
//--------------------------------------------------------------------
ifstream in_file(file_nm); // try to open data file
if (!in_file) {
cerr << "File name: " << file_nm << "\n";
Error("Cannot open input data/query file", ANNabort);
}
// allocate storage for points
if (pa != NULL) annDeallocPts(pa); // get rid of old points
pa = annAllocPts(n, dim);
for (i = 0; i < n; i++) { // read the data
if (!(in_file >> pa[i][0])) break;
for (int d = 1; d < dim; d++) {
in_file >> pa[i][d];
}
}
char ignore_me; // character for EOF test
in_file >> ignore_me; // try to get one more character
if (!in_file.eof()) { // exhausted space before eof
if (type == DATA)
Error("`data_size' too small. Input file truncated.", ANNwarn);
else
Error("`query_size' too small. Input file truncated.", ANNwarn);
}
n = i; // number of points read
//--------------------------------------------------------------------
// Print summary
//--------------------------------------------------------------------
if (stats > SILENT) {
if (type == DATA) {
cout << "[Read Data Points:\n";
cout << " data_size = " << n << "\n";
}
else {
cout << "[Read Query Points:\n";
cout << " query_size = " << n << "\n";
}
cout << " file_name = " << file_nm << "\n";
cout << " dim = " << dim << "\n";
// print if results requested
if ((type == DATA && stats >= SHOW_PTS) ||
(type == QUERY && stats >= QUERY_RES)) {
cout << " (Points:\n";
for (i = 0; i < n; i++) {
cout << " " << i << "\t";
printPoint(pa[i], dim);
cout << "\n";
}
cout << " )\n";
}
cout << "]\n";
}
}
//------------------------------------------------------------------------
// getTrueNN
// Computes the true nearest neighbors. For purposes of validation,
// this intentionally done in a rather dumb (but safe way), by
// invoking the brute-force search.
//
// The number of true nearest neighbors is somewhat larger than
// the number of nearest neighbors. This is so that the validation
// can determine the expected difference in element ranks.
//
// This procedure is invoked just prior to running queries. Since
// the operation takes a long time, it is performed only if needed.
// In particular, once generated, it will be regenerated only if
// new query or data points are generated, or if the requested number
// of true near neighbors or approximate near neighbors has changed.
//
// To validate fixed-radius searching, we compute two counts, one
// with the original query radius (trueSqRadius) and the other with
// a radius shrunken by the error factor (minSqradius). We then
// check that the count of points inside the approximate range is
// between these two bounds. Because fixed-radius search is
// allowed to ignore points within the shrunken radius, we only
// compute exact neighbors within this smaller distance (for we
// cannot guarantee that we will even visit the other points).
//------------------------------------------------------------------------
void getTrueNN() // compute true nearest neighbors
{
if (stats > SILENT) {
cout << "(Computing true nearest neighbors for validation. This may take time.)\n";
}
// deallocate existing storage
if (true_nn_idx != NULL) delete [] true_nn_idx;
if (true_dists != NULL) delete [] true_dists;
if (min_pts_in_range != NULL) delete [] min_pts_in_range;
if (max_pts_in_range != NULL) delete [] max_pts_in_range;
if (true_nn > data_size) { // can't get more nn than points
true_nn = data_size;
}
// allocate true answer storage
true_nn_idx = new ANNidx[true_nn*query_size];
true_dists = new ANNdist[true_nn*query_size];
min_pts_in_range = new int[query_size];
max_pts_in_range = new int[query_size];
ANNidxArray curr_nn_idx = true_nn_idx; // current locations in arrays
ANNdistArray curr_dists = true_dists;
// allocate search structure
ANNbruteForce *the_brute = new ANNbruteForce(data_pts, data_size, dim);
// compute nearest neighbors
for (int i = 0; i < query_size; i++) {
if (radius_bound == 0) { // standard kNN search
the_brute->annkSearch( // compute true near neighbors
query_pts[i], // query point
true_nn, // number of nearest neighbors
curr_nn_idx, // where to put indices
curr_dists); // where to put distances
}
else { // fixed radius kNN search
// search radii limits
ANNdist trueSqRadius = ANN_POW(radius_bound);
ANNdist minSqRadius = ANN_POW(radius_bound / (1+epsilon));
min_pts_in_range[i] = the_brute->annkFRSearch(
query_pts[i], // query point
minSqRadius, // shrunken search radius
true_nn, // number of near neighbors
curr_nn_idx, // nearest neighbors (returned)
curr_dists); // distance (returned)
max_pts_in_range[i] = the_brute->annkFRSearch(
query_pts[i], // query point
trueSqRadius, // true search radius
0, NULL, NULL); // (ignore kNN info)
}
curr_nn_idx += true_nn; // increment nn index pointer
curr_dists += true_nn; // increment nn dist pointer
}
delete the_brute; // delete brute-force struct
valid_dirty = ANNfalse; // validation good for now
}
//------------------------------------------------------------------------
// doValidation
// Compares the approximate answers to the k-nearest neighbors
// against the true nearest neighbors (computed earlier). It is
// assumed that the true nearest neighbors and indices have been
// computed earlier.
//
// First, we check that all the results are within their allowed
// limits, and generate an internal error, if not. For the sake of
// performance evaluation, we also compute the following two
// quantities for nearest neighbors:
//
// Average Error
// -------------
// The relative error between the distance to a reported nearest
// neighbor and the true nearest neighbor (of the same rank),
//
// Rank Error
// ----------
// The difference in rank between the reported nearest neighbor and
// its position (if any) among the true nearest neighbors. If we
// cannot find this point among the true nearest neighbors, then
// it assumed that the rank of the true nearest neighbor is true_nn+1.
//
// Because of the possibility of duplicate distances, this is computed
// as follows. For the j-th reported nearest neighbor, we count the
// number of true nearest neighbors that are at least this close. Let
// this be rnk. Then the rank error is max(0, j-rnk). (In the code
// below, j is an array index and so the first item is 0, not 1. Thus
// we take max(0, j+1-rnk) instead.)
//
// For the results of fixed-radious range count, we verify that the
// reported number of points in the range lies between the actual
// number of points in the shrunken and the true search radius.
//------------------------------------------------------------------------
void doValidation() // perform validation
{
int* curr_apx_idx = apx_nn_idx; // approx index pointer
ANNdistArray curr_apx_dst = apx_dists; // approx distance pointer
int* curr_tru_idx = true_nn_idx; // true index pointer
ANNdistArray curr_tru_dst = true_dists; // true distance pointer
int i, j;
if (true_nn < near_neigh) {
Error("Cannot validate with fewer true near neighbors than actual", ANNabort);
}
for (i = 0; i < query_size; i++) { // validate each query
//----------------------------------------------------------------
// Compute result errors
// In fixed radius search it is possible that not all k
// nearest neighbors were computed. Because the true
// results are computed over the shrunken radius, we should
// have at least as many true nearest neighbors as
// approximate nearest neighbors. (If not, an infinite
// error will be generated, and so an internal error will
// will be generated.
//
// Because nearest neighbors are sorted in increasing order
// of distance, as soon as we see a null index, we can
// terminate the distance checking. The error in the
// result should not exceed epsilon. However, if
// max_pts_visit is nonzero (meaning that the search is
// terminated early) this might happen.
//----------------------------------------------------------------
for (j = 0; j < near_neigh; j++) {
if (curr_tru_idx[j] == ANN_NULL_IDX)// no more true neighbors?
break;
// true i-th smallest distance
double true_dist = ANN_ROOT(curr_tru_dst[j]);
// reported i-th smallest
double rept_dist = ANN_ROOT(curr_apx_dst[j]);
// better than optimum?
if (rept_dist < true_dist*(1-ERR)) {
Error("INTERNAL ERROR: True nearest neighbor incorrect",
ANNabort);
}
double resultErr; // result error
if (true_dist == 0.0) { // let's not divide by zero
if (rept_dist != 0.0) resultErr = ANN_DBL_MAX;
else resultErr = 0.0;
}
else {
resultErr = (rept_dist - true_dist) / ((double) true_dist);
}
if (resultErr > epsilon && max_pts_visit == 0) {
Error("INTERNAL ERROR: Actual error exceeds epsilon",
ANNabort);
}
#ifdef ANN_PERF
ann_average_err += resultErr; // update statistics error
#endif
}
//--------------------------------------------------------------------
// Compute rank errors (only needed for perf measurements)
//--------------------------------------------------------------------
#ifdef ANN_PERF
for (j = 0; j < near_neigh; j++) {
if (curr_tru_idx[i] == ANN_NULL_IDX) // no more true neighbors?
break;
double rnkErr = 0.0; // rank error
// reported j-th distance
ANNdist rept_dist = curr_apx_dst[j];
int rnk = 0; // compute rank of this item
while (rnk < true_nn && curr_tru_dst[rnk] <= rept_dist)
rnk++;
if (j+1-rnk > 0) rnkErr = (double) (j+1-rnk);
ann_rank_err += rnkErr; // update average rank error
}
#endif
//----------------------------------------------------------------
// Check range counts from fixed-radius query
//----------------------------------------------------------------
if (radius_bound != 0) { // fixed-radius search
if (apx_pts_in_range[i] < min_pts_in_range[i] ||
apx_pts_in_range[i] > max_pts_in_range[i])
Error("INTERNAL ERROR: Invalid fixed-radius range count",
ANNabort);
}
curr_apx_idx += near_neigh;
curr_apx_dst += near_neigh;
curr_tru_idx += true_nn; // increment current pointers
curr_tru_dst += true_nn;
}
}
//----------------------------------------------------------------------
// treeStats
// Computes a number of statistics related to kd_trees and
// bd_trees. These statistics are printed if in verbose mode,
// and otherwise they are only printed if they are deemed to
// be outside of reasonable operating bounds.
//----------------------------------------------------------------------
#define log2(x) (log(x)/log(2.0)) // log base 2
void treeStats(
ostream &out, // output stream
ANNbool verbose) // print stats
{
const int MIN_PTS = 20; // min no. pts for checking
const float MAX_FRAC_TL = 0.50; // max frac of triv leaves
const float MAX_AVG_AR = 20; // max average aspect ratio
ANNkdStats st; // statistics structure
the_tree->getStats(st); // get statistics
// total number of nodes
int n_nodes = st.n_lf + st.n_spl + st.n_shr;
// should be O(n/bs)
int opt_n_nodes = (int) (2*(float(st.n_pts)/st.bkt_size));
int too_many_nodes = 10*opt_n_nodes;
if (st.n_pts >= MIN_PTS && n_nodes > too_many_nodes) {
out << "-----------------------------------------------------------\n";
out << "Warning: The tree has more than 10x as many nodes as points.\n";
out << "You may want to consider a different split or shrink method.\n";
out << "-----------------------------------------------------------\n";
verbose = ANNtrue;
}
// fraction of trivial leaves
float frac_tl = (st.n_lf == 0 ? 0 : ((float) st.n_tl)/ st.n_lf);
if (st.n_pts >= MIN_PTS && frac_tl > MAX_FRAC_TL) {
out << "-----------------------------------------------------------\n";
out << "Warning: A significant fraction of leaves contain no points.\n";
out << "You may want to consider a different split or shrink method.\n";
out << "-----------------------------------------------------------\n";
verbose = ANNtrue;
}
// depth should be O(dim*log n)
int too_many_levels = (int) (2.0 * st.dim * log2((double) st.n_pts));
int opt_levels = (int) log2(double(st.n_pts)/st.bkt_size);
if (st.n_pts >= MIN_PTS && st.depth > too_many_levels) {
out << "-----------------------------------------------------------\n";
out << "Warning: The tree is more than 2x as deep as (dim*log n).\n";
out << "You may want to consider a different split or shrink method.\n";
out << "-----------------------------------------------------------\n";
verbose = ANNtrue;
}
// average leaf aspect ratio
if (st.n_pts >= MIN_PTS && st.avg_ar > MAX_AVG_AR) {
out << "-----------------------------------------------------------\n";
out << "Warning: Average aspect ratio of cells is quite large.\n";
out << "This may slow queries depending on the point distribution.\n";
out << "-----------------------------------------------------------\n";
verbose = ANNtrue;
}
//------------------------------------------------------------------
// Print summaries if requested
//------------------------------------------------------------------
if (verbose) { // output statistics
out << " (Structure Statistics:\n";
out << " n_nodes = " << n_nodes
<< " (opt = " << opt_n_nodes
<< ", best if < " << too_many_nodes << ")\n"
<< " n_leaves = " << st.n_lf
<< " (" << st.n_tl << " contain no points)\n"
<< " n_splits = " << st.n_spl << "\n"
<< " n_shrinks = " << st.n_shr << "\n";
out << " empty_leaves = " << frac_tl*100
<< " percent (best if < " << MAX_FRAC_TL*100 << " percent)\n";
out << " depth = " << st.depth
<< " (opt = " << opt_levels
<< ", best if < " << too_many_levels << ")\n";
out << " avg_aspect_ratio = " << st.avg_ar
<< " (best if < " << MAX_AVG_AR << ")\n";
out << " )\n";
}
}