[37] | 1 | //---------------------------------------------------------------------- |
---|
| 2 | // File: ann_test.cpp |
---|
| 3 | // Programmer: Sunil Arya and David Mount |
---|
| 4 | // Description: test program for ANN (approximate nearest neighbors) |
---|
| 5 | // Last modified: 08/04/06 (Version 1.1.1) |
---|
| 6 | //---------------------------------------------------------------------- |
---|
| 7 | // Copyright (c) 1997-2005 University of Maryland and Sunil Arya and |
---|
| 8 | // David Mount. All Rights Reserved. |
---|
| 9 | // |
---|
| 10 | // This software and related documentation is part of the Approximate |
---|
| 11 | // Nearest Neighbor Library (ANN). This software is provided under |
---|
| 12 | // the provisions of the Lesser GNU Public License (LGPL). See the |
---|
| 13 | // file ../ReadMe.txt for further information. |
---|
| 14 | // |
---|
| 15 | // The University of Maryland (U.M.) and the authors make no |
---|
| 16 | // representations about the suitability or fitness of this software for |
---|
| 17 | // any purpose. It is provided "as is" without express or implied |
---|
| 18 | // warranty. |
---|
| 19 | //---------------------------------------------------------------------- |
---|
| 20 | // History: |
---|
| 21 | // Revision 0.1 03/04/98 |
---|
| 22 | // Initial release |
---|
| 23 | // Revision 0.2 06/26/98 |
---|
| 24 | // Added CLOCKS_PER_SEC definition if needed |
---|
| 25 | // Revision 1.0 04/01/05 |
---|
| 26 | // Added comments (from "#" to eol) |
---|
| 27 | // Added clus_orth_flats and clus_ellipsoids distributions |
---|
| 28 | // Fixed order of fair and midpt in split_table |
---|
| 29 | // Added dump/load operations |
---|
| 30 | // Cleaned up C++ for modern compilers |
---|
| 31 | // Revision 1.1 05/03/05 |
---|
| 32 | // Added fixed radius kNN search |
---|
| 33 | // Revision 1.1.1 08/04/06 |
---|
| 34 | // Added planted distribution |
---|
| 35 | //---------------------------------------------------------------------- |
---|
| 36 | |
---|
| 37 | #include <ctime> // clock |
---|
| 38 | #include <cmath> // math routines |
---|
| 39 | #include <string> // C string ops |
---|
| 40 | #include <fstream> // file I/O |
---|
| 41 | |
---|
| 42 | #include <ANN/ANN.h> // ANN declarations |
---|
| 43 | #include <ANN/ANNx.h> // more ANN declarations |
---|
| 44 | #include <ANN/ANNperf.h> // performance evaluation |
---|
| 45 | |
---|
| 46 | #include "rand.h" // random point generation |
---|
| 47 | |
---|
| 48 | #ifndef CLOCKS_PER_SEC // define clocks-per-second if needed |
---|
| 49 | #define CLOCKS_PER_SEC 1000000 |
---|
| 50 | #endif |
---|
| 51 | |
---|
| 52 | using namespace std; // make std:: available |
---|
| 53 | |
---|
| 54 | //---------------------------------------------------------------------- |
---|
| 55 | // ann_test |
---|
| 56 | // |
---|
| 57 | // This program is a driver for testing and evaluating the ANN library |
---|
| 58 | // for computing approximate nearest neighbors. It allows the user to |
---|
| 59 | // generate data and query sets of various sizes, dimensions, and |
---|
| 60 | // distributions, to build kd- and bbd-trees of various types, and then |
---|
| 61 | // run queries and outputting various performance statistics. |
---|
| 62 | // |
---|
| 63 | // Overview: |
---|
| 64 | // --------- |
---|
| 65 | // The test program is run as follows: |
---|
| 66 | // |
---|
| 67 | // ann_test < test_input > test_output |
---|
| 68 | // |
---|
| 69 | // where the test_input file contains a list of directives as described |
---|
| 70 | // below. Directives consist of a directive name, followed by list of |
---|
| 71 | // arguments (depending on the directive). Arguments and directives are |
---|
| 72 | // separated by white space (blank, tab, and newline). String arguments |
---|
| 73 | // are not quoted, and consist of a string of nonwhite chacters. A |
---|
| 74 | // character "#" denotes a comment. The following characters up to |
---|
| 75 | // the end of line are ignored. Comments may only be inserted between |
---|
| 76 | // directives (not within the argument list of a directive). |
---|
| 77 | // |
---|
| 78 | // Basic operations: |
---|
| 79 | // ----------------- |
---|
| 80 | // The test program can perform the following operations. How these |
---|
| 81 | // operations are performed depends on the options which are described |
---|
| 82 | // later. |
---|
| 83 | // |
---|
| 84 | // Data Generation: |
---|
| 85 | // ---------------- |
---|
| 86 | // read_data_pts <file> Create a set of data points whose |
---|
| 87 | // coordinates are input from file <file>. |
---|
| 88 | // gen_data_pts Create a set of data points whose |
---|
| 89 | // coordinates are generated from the |
---|
| 90 | // current point distribution. |
---|
| 91 | // |
---|
| 92 | // Building the tree: |
---|
| 93 | // ------------------ |
---|
| 94 | // build_ann Generate an approximate nearest neighbor |
---|
| 95 | // structure for the current data set, using |
---|
| 96 | // the selected splitting rules. Any existing |
---|
| 97 | // tree will be destroyed. |
---|
| 98 | // |
---|
| 99 | // Query Generation/Searching: |
---|
| 100 | // --------------------------- |
---|
| 101 | // read_query_pts <file> Create a set of query points whose |
---|
| 102 | // coordinates are input from file <file>. |
---|
| 103 | // gen_query_pts Create a set of query points whose |
---|
| 104 | // coordinates are generated from the |
---|
| 105 | // current point distribution. |
---|
| 106 | // run_queries <string> Apply nearest neighbor searching to the |
---|
| 107 | // query points using the approximate nearest |
---|
| 108 | // neighbor structure and the given search |
---|
| 109 | // strategy. Possible strategies are: |
---|
| 110 | // standard = standard kd-tree search |
---|
| 111 | // priority = priority search |
---|
| 112 | // |
---|
| 113 | // Miscellaneous: |
---|
| 114 | // -------------- |
---|
| 115 | // output_label Output a label to the output file. |
---|
| 116 | // dump <file> Dump the current structure to given file. |
---|
| 117 | // (The dump format is explained further in |
---|
| 118 | // the source file kd_tree.cc.) |
---|
| 119 | // load <file> Load a tree from a data file which was |
---|
| 120 | // created by the dump operation. Any |
---|
| 121 | // existing tree will be destroyed. |
---|
| 122 | // |
---|
| 123 | // Options: |
---|
| 124 | // -------- |
---|
| 125 | // How these operations are performed depends on a set of options. |
---|
| 126 | // If an option is not specified, a default value is used. An option |
---|
| 127 | // retains its value until it is set again. String inputs are not |
---|
| 128 | // enclosed in quotes, and must contain no embedded white space (sorry, |
---|
| 129 | // this is C++'s convention). |
---|
| 130 | // |
---|
| 131 | // Options affecting search tree structure: |
---|
| 132 | // ---------------------------------------- |
---|
| 133 | // split_rule <type> Type of splitting rule to use in building |
---|
| 134 | // the search tree. Choices are: |
---|
| 135 | // kd = optimized kd-tree |
---|
| 136 | // midpt = midpoint split |
---|
| 137 | // fair = fair split |
---|
| 138 | // sl_midpt = sliding midpt split |
---|
| 139 | // sl_fair = sliding fair split |
---|
| 140 | // suggest = authors' choice for best |
---|
| 141 | // The default is "suggest". See the file |
---|
| 142 | // kd_split.cc for more detailed information. |
---|
| 143 | // |
---|
| 144 | // shrink_rule <type> Type of shrinking rule to use in building |
---|
| 145 | // a bd-tree data structure. If "none" is |
---|
| 146 | // given, then no shrinking is performed and |
---|
| 147 | // the result is a kd-tree. Choices are: |
---|
| 148 | // none = perform no shrinking |
---|
| 149 | // simple = simple shrinking |
---|
| 150 | // centroid = centroid shrinking |
---|
| 151 | // suggest = authors' choice for best |
---|
| 152 | // The default is "none". See the file |
---|
| 153 | // bd_tree.cc for more information. |
---|
| 154 | // bucket_size <int> Bucket size, that is, the maximum number of |
---|
| 155 | // points stored in each leaf node. |
---|
| 156 | // |
---|
| 157 | // Options affecting data and query point generation: |
---|
| 158 | // -------------------------------------------------- |
---|
| 159 | // dim <int> Dimension of space. |
---|
| 160 | // seed <int> Seed for random number generation. |
---|
| 161 | // data_size <int> Number of data points. When reading data |
---|
| 162 | // points from a file, this indicates the |
---|
| 163 | // maximum number of points for storage |
---|
| 164 | // allocation. Default = 100. |
---|
| 165 | // query_size <int> Same as data_size for query points. |
---|
| 166 | // std_dev <float> Standard deviation (used in gauss, |
---|
| 167 | // planted, and clustered distributions). |
---|
| 168 | // This is the "small" distribution for |
---|
| 169 | // clus_ellipsoids. Default = 1. |
---|
| 170 | // std_dev_lo <float> Low and high standard deviations (used in |
---|
| 171 | // std_dev_hi <float> clus_ellipsoids). Default = 1. |
---|
| 172 | // corr_coef <float> Correlation coefficient (used in co-gauss |
---|
| 173 | // and co_lapace distributions). Default = 0.05. |
---|
| 174 | // colors <int> Number of color classes (clusters) (used |
---|
| 175 | // in the clustered distributions). Default = 5. |
---|
| 176 | // new_clust Once generated, cluster centers are not |
---|
| 177 | // normally regenerated. This is so that both |
---|
| 178 | // query points and data points can be generated |
---|
| 179 | // using the same set of clusters. This option |
---|
| 180 | // forces new cluster centers to be generated |
---|
| 181 | // with the next generation of either data or |
---|
| 182 | // query points. |
---|
| 183 | // max_clus_dim <int> Maximum dimension of clusters (used in |
---|
| 184 | // clus_orth_flats and clus_ellipsoids). |
---|
| 185 | // Default = 1. |
---|
| 186 | // distribution <string> Type of input distribution |
---|
| 187 | // uniform = uniform over cube [-1,1]^d. |
---|
| 188 | // gauss = Gaussian with mean 0 |
---|
| 189 | // laplace = Laplacian, mean 0 and var 1 |
---|
| 190 | // co_gauss = correlated Gaussian |
---|
| 191 | // co_laplace = correlated Laplacian |
---|
| 192 | // clus_gauss = clustered Gaussian |
---|
| 193 | // clus_orth_flats = clusters of orth flats |
---|
| 194 | // clus_ellipsoids = clusters of ellipsoids |
---|
| 195 | // planted = planted distribution |
---|
| 196 | // See the file rand.cpp for further information. |
---|
| 197 | // |
---|
| 198 | // Options affecting nearest neighbor search: |
---|
| 199 | // ------------------------------------------ |
---|
| 200 | // epsilon <float> Error bound for approx. near neigh. search. |
---|
| 201 | // near_neigh <int> Number of nearest neighbors to compute. |
---|
| 202 | // max_pts_visit <int> Maximum number of points to visit before |
---|
| 203 | // terminating. (Used in applications where |
---|
| 204 | // real-time performance is important.) |
---|
| 205 | // (Default = 0, which means no limit.) |
---|
| 206 | // radius_bound <float> Sets an upper bound on the nearest |
---|
| 207 | // neighbor search radius. If the bound is |
---|
| 208 | // positive, then fixed-radius nearest |
---|
| 209 | // neighbor searching is performed, and the |
---|
| 210 | // count of the number of points in the |
---|
| 211 | // range is returned. If the bound is |
---|
| 212 | // zero, then standard search is used. |
---|
| 213 | // This can only be used with standard, not |
---|
| 214 | // priority, search. (Default = 0, which |
---|
| 215 | // means standard search.) |
---|
| 216 | // |
---|
| 217 | // Options affection general program behavior: |
---|
| 218 | // ------------------------------------------- |
---|
| 219 | // stats <string> Level of statistics output |
---|
| 220 | // silent = no output, |
---|
| 221 | // exec_time += execution time only |
---|
| 222 | // prep_stats += preprocessing statistics |
---|
| 223 | // query_stats += query performance stats |
---|
| 224 | // query_res += results of queries |
---|
| 225 | // show_pts += show the data points |
---|
| 226 | // show_struct += print search structure |
---|
| 227 | // validate <string> Validate experiment and compute average |
---|
| 228 | // error. Since validation causes exact |
---|
| 229 | // nearest neighbors to be computed by the |
---|
| 230 | // brute force method, this can take a long |
---|
| 231 | // time. Valid arguments are: |
---|
| 232 | // on = turn validation on |
---|
| 233 | // off = turn validation off |
---|
| 234 | // true_near_neigh <int> Number of true nearest neighbors to compute. |
---|
| 235 | // When validating, we compute the difference |
---|
| 236 | // in rank between each reported nearest neighbor |
---|
| 237 | // and the true nearest neighbor of the same |
---|
| 238 | // rank. Thus it is necessary to compute a |
---|
| 239 | // few more true nearest neighbors. By default |
---|
| 240 | // we compute 10 more than near_neigh. With |
---|
| 241 | // this option the exact number can be set. |
---|
| 242 | // (Used only when validating.) |
---|
| 243 | // |
---|
| 244 | // Example: |
---|
| 245 | // -------- |
---|
| 246 | // output_label test_run_0 # output label for this run |
---|
| 247 | // validate off # do not perform validation |
---|
| 248 | // dim 16 # points in dimension 16 |
---|
| 249 | // stats query_stats # output performance statistics for queries |
---|
| 250 | // seed 121212 # random number seed |
---|
| 251 | // data_size 1000 |
---|
| 252 | // distribution uniform |
---|
| 253 | // gen_data_pts # 1000 uniform data points in dim 16 |
---|
| 254 | // query_size 100 |
---|
| 255 | // std_dev 0.05 |
---|
| 256 | // distribution clus_gauss |
---|
| 257 | // gen_query_pts # 100 points in 10 clusters with std_dev 0.05 |
---|
| 258 | // bucket_size 2 |
---|
| 259 | // split_rule kd |
---|
| 260 | // shrink_rule none |
---|
| 261 | // build_ann # kd-tree, bucket size 2 |
---|
| 262 | // epsilon 0.1 |
---|
| 263 | // near_neigh 5 |
---|
| 264 | // max_pts_visit 100 # stop search if more than 100 points seen |
---|
| 265 | // run_queries standard # run queries; 5 nearest neighbors, 10% error |
---|
| 266 | // data_size 500 |
---|
| 267 | // read_data_pts data.in # read up to 500 points from file data.in |
---|
| 268 | // split_rule sl_midpt |
---|
| 269 | // shrink_rule simple |
---|
| 270 | // build_ann # bd-tree; simple shrink, sliding midpoint split |
---|
| 271 | // epsilon 0 |
---|
| 272 | // run_queries priority # run same queries; 0 allowable error |
---|
| 273 | // |
---|
| 274 | //------------------------------------------------------------------------ |
---|
| 275 | |
---|
| 276 | //------------------------------------------------------------------------ |
---|
| 277 | // Constants |
---|
| 278 | //------------------------------------------------------------------------ |
---|
| 279 | |
---|
| 280 | const int STRING_LEN = 500; // max string length |
---|
| 281 | const double ERR = 0.00001; // epsilon (for float compares) |
---|
| 282 | |
---|
| 283 | //------------------------------------------------------------------------ |
---|
| 284 | // Enumerated values and conversions |
---|
| 285 | //------------------------------------------------------------------------ |
---|
| 286 | |
---|
| 287 | typedef enum {DATA, QUERY} PtType; // point types |
---|
| 288 | |
---|
| 289 | //------------------------------------------------------------------------ |
---|
| 290 | // Statistics output levels |
---|
| 291 | //------------------------------------------------------------------------ |
---|
| 292 | |
---|
| 293 | typedef enum { // stat levels |
---|
| 294 | SILENT, // no output |
---|
| 295 | EXEC_TIME, // just execution time |
---|
| 296 | PREP_STATS, // preprocessing info |
---|
| 297 | QUERY_STATS, // query performance |
---|
| 298 | QUERY_RES, // query results |
---|
| 299 | SHOW_PTS, // show data points |
---|
| 300 | SHOW_STRUCT, // show tree structure |
---|
| 301 | N_STAT_LEVELS} // number of levels |
---|
| 302 | StatLev; |
---|
| 303 | |
---|
| 304 | const char stat_table[N_STAT_LEVELS][STRING_LEN] = { |
---|
| 305 | "silent", // SILENT |
---|
| 306 | "exec_time", // EXEC_TIME |
---|
| 307 | "prep_stats", // PREP_STATS |
---|
| 308 | "query_stats", // QUERY_STATS |
---|
| 309 | "query_res", // QUERY_RES |
---|
| 310 | "show_pts", // SHOW_PTS |
---|
| 311 | "show_struct"}; // SHOW_STRUCT |
---|
| 312 | |
---|
| 313 | //------------------------------------------------------------------------ |
---|
| 314 | // Distributions |
---|
| 315 | //------------------------------------------------------------------------ |
---|
| 316 | |
---|
| 317 | typedef enum { // distributions |
---|
| 318 | UNIFORM, // uniform over cube [-1,1]^d. |
---|
| 319 | GAUSS, // Gaussian with mean 0 |
---|
| 320 | LAPLACE, // Laplacian, mean 0 and var 1 |
---|
| 321 | CO_GAUSS, // correlated Gaussian |
---|
| 322 | CO_LAPLACE, // correlated Laplacian |
---|
| 323 | CLUS_GAUSS, // clustered Gaussian |
---|
| 324 | CLUS_ORTH_FLATS, // clustered on orthog flats |
---|
| 325 | CLUS_ELLIPSOIDS, // clustered on ellipsoids |
---|
| 326 | PLANTED, // planted distribution |
---|
| 327 | N_DISTRIBS} |
---|
| 328 | Distrib; |
---|
| 329 | |
---|
| 330 | const char distr_table[N_DISTRIBS][STRING_LEN] = { |
---|
| 331 | "uniform", // UNIFORM |
---|
| 332 | "gauss", // GAUSS |
---|
| 333 | "laplace", // LAPLACE |
---|
| 334 | "co_gauss", // CO_GAUSS |
---|
| 335 | "co_laplace", // CO_LAPLACE |
---|
| 336 | "clus_gauss", // CLUS_GAUSS |
---|
| 337 | "clus_orth_flats", // CLUS_ORTH_FLATS |
---|
| 338 | "clus_ellipsoids", // CLUS_ELLIPSOIS |
---|
| 339 | "planted"}; // PLANTED |
---|
| 340 | |
---|
| 341 | //------------------------------------------------------------------------ |
---|
| 342 | // Splitting rules for kd-trees (see ANN.h for types) |
---|
| 343 | //------------------------------------------------------------------------ |
---|
| 344 | |
---|
| 345 | const int N_SPLIT_RULES = 6; |
---|
| 346 | const char split_table[N_SPLIT_RULES][STRING_LEN] = { |
---|
| 347 | "standard", // standard optimized kd-tree |
---|
| 348 | "midpt", // midpoint split |
---|
| 349 | "fair", // fair split |
---|
| 350 | "sl_midpt", // sliding midpt split |
---|
| 351 | "sl_fair", // sliding fair split |
---|
| 352 | "suggest"}; // authors' choice for best |
---|
| 353 | |
---|
| 354 | //------------------------------------------------------------------------ |
---|
| 355 | // Shrinking rules for bd-trees (see ANN.h for types) |
---|
| 356 | //------------------------------------------------------------------------ |
---|
| 357 | |
---|
| 358 | const int N_SHRINK_RULES = 4; |
---|
| 359 | const char shrink_table[N_SHRINK_RULES][STRING_LEN] = { |
---|
| 360 | "none", // perform no shrinking (kd-tree) |
---|
| 361 | "simple", // simple shrinking |
---|
| 362 | "centroid", // centroid shrinking |
---|
| 363 | "suggest"}; // authors' choice for best |
---|
| 364 | |
---|
| 365 | //---------------------------------------------------------------------- |
---|
| 366 | // Short utility functions |
---|
| 367 | // Error - general error routine |
---|
| 368 | // printPoint - print a point to standard output |
---|
| 369 | // lookUp - look up a name in table and return index |
---|
| 370 | //---------------------------------------------------------------------- |
---|
| 371 | |
---|
| 372 | void Error( // error routine |
---|
| 373 | char *msg, // error message |
---|
| 374 | ANNerr level) // abort afterwards |
---|
| 375 | { |
---|
| 376 | if (level == ANNabort) { |
---|
| 377 | cerr << "ann_test: ERROR------->" << msg << "<-------------ERROR\n"; |
---|
| 378 | exit(1); |
---|
| 379 | } |
---|
| 380 | else { |
---|
| 381 | cerr << "ann_test: WARNING----->" << msg << "<-------------WARNING\n"; |
---|
| 382 | } |
---|
| 383 | } |
---|
| 384 | |
---|
| 385 | void printPoint( // print point |
---|
| 386 | ANNpoint p, // the point |
---|
| 387 | int dim) // the dimension |
---|
| 388 | { |
---|
| 389 | cout << "["; |
---|
| 390 | for (int i = 0; i < dim; i++) { |
---|
| 391 | cout << p[i]; |
---|
| 392 | if (i < dim-1) cout << ","; |
---|
| 393 | } |
---|
| 394 | cout << "]"; |
---|
| 395 | } |
---|
| 396 | |
---|
| 397 | int lookUp( // look up name in table |
---|
| 398 | const char *arg, // name to look up |
---|
| 399 | const char (*table)[STRING_LEN], // name table |
---|
| 400 | int size) // table size |
---|
| 401 | { |
---|
| 402 | int i; |
---|
| 403 | for (i = 0; i < size; i++) { |
---|
| 404 | if (!strcmp(arg, table[i])) return i; |
---|
| 405 | } |
---|
| 406 | return i; |
---|
| 407 | } |
---|
| 408 | |
---|
| 409 | //------------------------------------------------------------------------ |
---|
| 410 | // Function declarations |
---|
| 411 | //------------------------------------------------------------------------ |
---|
| 412 | |
---|
| 413 | void generatePts( // generate data/query points |
---|
| 414 | ANNpointArray &pa, // point array (returned) |
---|
| 415 | int n, // number of points |
---|
| 416 | PtType type, // point type |
---|
| 417 | ANNbool new_clust, // new cluster centers desired? |
---|
| 418 | ANNpointArray src = NULL, // source array (for PLANTED) |
---|
| 419 | int n_src = 0); // source size (for PLANTED) |
---|
| 420 | |
---|
| 421 | void readPts( // read data/query points from file |
---|
| 422 | ANNpointArray &pa, // point array (returned) |
---|
| 423 | int &n, // number of points |
---|
| 424 | char *file_nm, // file name |
---|
| 425 | PtType type); // point type (DATA, QUERY) |
---|
| 426 | |
---|
| 427 | void doValidation(); // perform validation |
---|
| 428 | void getTrueNN(); // compute true nearest neighbors |
---|
| 429 | |
---|
| 430 | void treeStats( // print statistics on kd- or bd-tree |
---|
| 431 | ostream &out, // output stream |
---|
| 432 | ANNbool verbose); // print stats |
---|
| 433 | |
---|
| 434 | //------------------------------------------------------------------------ |
---|
| 435 | // Default execution parameters |
---|
| 436 | //------------------------------------------------------------------------ |
---|
| 437 | const int extra_nn = 10; // how many extra true nn's? |
---|
| 438 | |
---|
| 439 | const int def_dim = 2; // def dimension |
---|
| 440 | const int def_data_size = 100; // def data size |
---|
| 441 | const int def_query_size = 100; // def number of queries |
---|
| 442 | const int def_n_color = 5; // def number of colors |
---|
| 443 | const ANNbool def_new_clust = ANNfalse; // def new clusters flag |
---|
| 444 | const int def_max_dim = 1; // def max flat dimension |
---|
| 445 | const Distrib def_distr = UNIFORM; // def distribution |
---|
| 446 | const double def_std_dev = 1.00; // def standard deviation |
---|
| 447 | const double def_corr_coef = 0.05; // def correlation coef |
---|
| 448 | const int def_bucket_size = 1; // def bucket size |
---|
| 449 | const double def_epsilon = 0.0; // def error bound |
---|
| 450 | const int def_near_neigh = 1; // def number of near neighbors |
---|
| 451 | const int def_max_visit = 0; // def number of points visited |
---|
| 452 | const int def_rad_bound = 0; // def radius bound |
---|
| 453 | // def number of true nn's |
---|
| 454 | const int def_true_nn = def_near_neigh + extra_nn; |
---|
| 455 | const int def_seed = 0; // def seed for random numbers |
---|
| 456 | const ANNbool def_validate = ANNfalse; // def validation flag |
---|
| 457 | // def statistics output level |
---|
| 458 | const StatLev def_stats = QUERY_STATS; |
---|
| 459 | const ANNsplitRule // def splitting rule |
---|
| 460 | def_split = ANN_KD_SUGGEST; |
---|
| 461 | const ANNshrinkRule // def shrinking rule |
---|
| 462 | def_shrink = ANN_BD_NONE; |
---|
| 463 | |
---|
| 464 | //------------------------------------------------------------------------ |
---|
| 465 | // Global variables - Execution options |
---|
| 466 | //------------------------------------------------------------------------ |
---|
| 467 | |
---|
| 468 | int dim; // dimension |
---|
| 469 | int data_size; // data size |
---|
| 470 | int query_size; // number of queries |
---|
| 471 | int n_color; // number of colors |
---|
| 472 | ANNbool new_clust; // generate new clusters? |
---|
| 473 | int max_dim; // maximum flat dimension |
---|
| 474 | Distrib distr; // distribution |
---|
| 475 | double corr_coef; // correlation coef |
---|
| 476 | double std_dev; // standard deviation |
---|
| 477 | double std_dev_lo; // low standard deviation |
---|
| 478 | double std_dev_hi; // high standard deviation |
---|
| 479 | int bucket_size; // bucket size |
---|
| 480 | double epsilon; // error bound |
---|
| 481 | int near_neigh; // number of near neighbors |
---|
| 482 | int max_pts_visit; // max number of points to visit |
---|
| 483 | double radius_bound; // maximum radius search bound |
---|
| 484 | int true_nn; // number of true nn's |
---|
| 485 | ANNbool validate; // validation flag |
---|
| 486 | StatLev stats; // statistics output level |
---|
| 487 | ANNsplitRule split; // splitting rule |
---|
| 488 | ANNshrinkRule shrink; // shrinking rule |
---|
| 489 | |
---|
| 490 | //------------------------------------------------------------------------ |
---|
| 491 | // More globals - pointers to dynamically allocated arrays and structures |
---|
| 492 | // |
---|
| 493 | // It is assumed that all these values are set to NULL when nothing |
---|
| 494 | // is allocated. |
---|
| 495 | // |
---|
| 496 | // data_pts, query_pts The data and query points |
---|
| 497 | // the_tree Points to the kd- or bd-tree for |
---|
| 498 | // nearest neighbor searching. |
---|
| 499 | // apx_nn_idx, apx_dists Record approximate near neighbor |
---|
| 500 | // indices and distances |
---|
| 501 | // apx_pts_in_range Counts of the number of points in |
---|
| 502 | // the in approx range, for fixed- |
---|
| 503 | // radius NN searching. |
---|
| 504 | // true_nn_idx, true_dists Record true near neighbor |
---|
| 505 | // indices and distances |
---|
| 506 | // min_pts_in_range, max_... Min and max counts of the number |
---|
| 507 | // of points in the in approximate |
---|
| 508 | // range. |
---|
| 509 | // valid_dirty To avoid repeated validation, |
---|
| 510 | // we only validate query results |
---|
| 511 | // once. This validation becomes |
---|
| 512 | // invalid, if a new tree, new data |
---|
| 513 | // points or new query points have |
---|
| 514 | // been generated. |
---|
| 515 | // tree_data_size The number of points in the |
---|
| 516 | // current tree. (This will be the |
---|
| 517 | // same a data_size unless points have |
---|
| 518 | // been added since the tree was |
---|
| 519 | // built.) |
---|
| 520 | // |
---|
| 521 | // The approximate and true nearest neighbor results are stored |
---|
| 522 | // in: apx_nn_idx, apx_dists, and true_nn_idx, true_dists. |
---|
| 523 | // They are really flattened 2-dimensional arrays. Each of these |
---|
| 524 | // arrays consists of query_size blocks, each of which contains |
---|
| 525 | // near_neigh (or true_nn) entries, one for each of the nearest |
---|
| 526 | // neighbors for a given query point. |
---|
| 527 | //------------------------------------------------------------------------ |
---|
| 528 | |
---|
| 529 | ANNpointArray data_pts; // data points |
---|
| 530 | ANNpointArray query_pts; // query points |
---|
| 531 | ANNbd_tree* the_tree; // kd- or bd-tree search structure |
---|
| 532 | ANNidxArray apx_nn_idx; // storage for near neighbor indices |
---|
| 533 | ANNdistArray apx_dists; // storage for near neighbor distances |
---|
| 534 | int* apx_pts_in_range; // storage for no. of points in range |
---|
| 535 | ANNidxArray true_nn_idx; // true near neighbor indices |
---|
| 536 | ANNdistArray true_dists; // true near neighbor distances |
---|
| 537 | int* min_pts_in_range; // min points in approx range |
---|
| 538 | int* max_pts_in_range; // max points in approx range |
---|
| 539 | |
---|
| 540 | ANNbool valid_dirty; // validation is no longer valid |
---|
| 541 | |
---|
| 542 | //------------------------------------------------------------------------ |
---|
| 543 | // Initialize global parameters |
---|
| 544 | //------------------------------------------------------------------------ |
---|
| 545 | |
---|
| 546 | void initGlobals() |
---|
| 547 | { |
---|
| 548 | dim = def_dim; // init execution parameters |
---|
| 549 | data_size = def_data_size; |
---|
| 550 | query_size = def_query_size; |
---|
| 551 | distr = def_distr; |
---|
| 552 | corr_coef = def_corr_coef; |
---|
| 553 | std_dev = def_std_dev; |
---|
| 554 | std_dev_lo = def_std_dev; |
---|
| 555 | std_dev_hi = def_std_dev; |
---|
| 556 | new_clust = def_new_clust; |
---|
| 557 | max_dim = def_max_dim; |
---|
| 558 | n_color = def_n_color; |
---|
| 559 | bucket_size = def_bucket_size; |
---|
| 560 | epsilon = def_epsilon; |
---|
| 561 | near_neigh = def_near_neigh; |
---|
| 562 | max_pts_visit = def_max_visit; |
---|
| 563 | radius_bound = def_rad_bound; |
---|
| 564 | true_nn = def_true_nn; |
---|
| 565 | validate = def_validate; |
---|
| 566 | stats = def_stats; |
---|
| 567 | split = def_split; |
---|
| 568 | shrink = def_shrink; |
---|
| 569 | annIdum = -def_seed; // init. global seed for ran0() |
---|
| 570 | |
---|
| 571 | data_pts = NULL; // initialize storage pointers |
---|
| 572 | query_pts = NULL; |
---|
| 573 | the_tree = NULL; |
---|
| 574 | apx_nn_idx = NULL; |
---|
| 575 | apx_dists = NULL; |
---|
| 576 | apx_pts_in_range = NULL; |
---|
| 577 | true_nn_idx = NULL; |
---|
| 578 | true_dists = NULL; |
---|
| 579 | min_pts_in_range = NULL; |
---|
| 580 | max_pts_in_range = NULL; |
---|
| 581 | |
---|
| 582 | valid_dirty = ANNtrue; // (validation must be done) |
---|
| 583 | } |
---|
| 584 | |
---|
| 585 | //------------------------------------------------------------------------ |
---|
| 586 | // getDirective - skip comments and read next directive |
---|
| 587 | // Returns ANNtrue if directive read, and ANNfalse if eof seen. |
---|
| 588 | //------------------------------------------------------------------------ |
---|
| 589 | |
---|
| 590 | ANNbool skipComment( // skip any comments |
---|
| 591 | istream &in) // input stream |
---|
| 592 | { |
---|
| 593 | char ch = 0; |
---|
| 594 | // skip whitespace |
---|
| 595 | do { in.get(ch); } while (isspace(ch) && !in.eof()); |
---|
| 596 | while (ch == '#' && !in.eof()) { // comment? |
---|
| 597 | // skip to end of line |
---|
| 598 | do { in.get(ch); } while(ch != '\n' && !in.eof()); |
---|
| 599 | // skip whitespace |
---|
| 600 | do { in.get(ch); } while(isspace(ch) && !in.eof()); |
---|
| 601 | } |
---|
| 602 | if (in.eof()) return ANNfalse; // end of file |
---|
| 603 | in.putback(ch); // put character back |
---|
| 604 | return ANNtrue; |
---|
| 605 | } |
---|
| 606 | |
---|
| 607 | ANNbool getDirective( |
---|
| 608 | istream &in, // input stream |
---|
| 609 | char *directive) // directive storage |
---|
| 610 | { |
---|
| 611 | if (!skipComment(in)) // skip comments |
---|
| 612 | return ANNfalse; // found eof along the way? |
---|
| 613 | in >> directive; // read directive |
---|
| 614 | return ANNtrue; |
---|
| 615 | } |
---|
| 616 | |
---|
| 617 | |
---|
| 618 | //------------------------------------------------------------------------ |
---|
| 619 | // main program - driver |
---|
| 620 | // The main program reads input options, invokes the necessary |
---|
| 621 | // routines to process them. |
---|
| 622 | //------------------------------------------------------------------------ |
---|
| 623 | |
---|
| 624 | int main(int argc, char** argv) |
---|
| 625 | { |
---|
| 626 | long clock0; // clock time |
---|
| 627 | char directive[STRING_LEN]; // input directive |
---|
| 628 | char arg[STRING_LEN]; // all-purpose argument |
---|
| 629 | |
---|
| 630 | cout << "------------------------------------------------------------\n" |
---|
| 631 | << "ann_test: Version " << ANNversion << " " << ANNversionCmt << "\n" |
---|
| 632 | << " Copyright: " << ANNcopyright << ".\n" |
---|
| 633 | << " Latest Revision: " << ANNlatestRev << ".\n" |
---|
| 634 | << "------------------------------------------------------------\n\n"; |
---|
| 635 | |
---|
| 636 | initGlobals(); // initialize global values |
---|
| 637 | |
---|
| 638 | //-------------------------------------------------------------------- |
---|
| 639 | // Main input loop |
---|
| 640 | //-------------------------------------------------------------------- |
---|
| 641 | // read input directive |
---|
| 642 | while (getDirective(cin, directive)) { |
---|
| 643 | //---------------------------------------------------------------- |
---|
| 644 | // Read options |
---|
| 645 | //---------------------------------------------------------------- |
---|
| 646 | if (!strcmp(directive,"dim")) { |
---|
| 647 | cin >> dim; |
---|
| 648 | } |
---|
| 649 | else if (!strcmp(directive,"colors")) { |
---|
| 650 | cin >> n_color; |
---|
| 651 | } |
---|
| 652 | else if (!strcmp(directive,"new_clust")) { |
---|
| 653 | new_clust = ANNtrue; |
---|
| 654 | } |
---|
| 655 | else if (!strcmp(directive,"max_clus_dim")) { |
---|
| 656 | cin >> max_dim; |
---|
| 657 | } |
---|
| 658 | else if (!strcmp(directive,"std_dev")) { |
---|
| 659 | cin >> std_dev; |
---|
| 660 | } |
---|
| 661 | else if (!strcmp(directive,"std_dev_lo")) { |
---|
| 662 | cin >> std_dev_lo; |
---|
| 663 | } |
---|
| 664 | else if (!strcmp(directive,"std_dev_hi")) { |
---|
| 665 | cin >> std_dev_hi; |
---|
| 666 | } |
---|
| 667 | else if (!strcmp(directive,"corr_coef")) { |
---|
| 668 | cin >> corr_coef; |
---|
| 669 | } |
---|
| 670 | else if (!strcmp(directive, "data_size")) { |
---|
| 671 | cin >> data_size; |
---|
| 672 | } |
---|
| 673 | else if (!strcmp(directive,"query_size")) { |
---|
| 674 | cin >> query_size; |
---|
| 675 | } |
---|
| 676 | else if (!strcmp(directive,"bucket_size")) { |
---|
| 677 | cin >> bucket_size; |
---|
| 678 | } |
---|
| 679 | else if (!strcmp(directive,"epsilon")) { |
---|
| 680 | cin >> epsilon; |
---|
| 681 | } |
---|
| 682 | else if (!strcmp(directive,"max_pts_visit")) { |
---|
| 683 | cin >> max_pts_visit; |
---|
| 684 | valid_dirty = ANNtrue; // validation must be redone |
---|
| 685 | } |
---|
| 686 | else if (!strcmp(directive,"radius_bound")) { |
---|
| 687 | cin >> radius_bound; |
---|
| 688 | valid_dirty = ANNtrue; // validation must be redone |
---|
| 689 | } |
---|
| 690 | else if (!strcmp(directive,"near_neigh")) { |
---|
| 691 | cin >> near_neigh; |
---|
| 692 | true_nn = near_neigh + extra_nn; // also reset true near neighs |
---|
| 693 | valid_dirty = ANNtrue; // validation must be redone |
---|
| 694 | } |
---|
| 695 | else if (!strcmp(directive,"true_near_neigh")) { |
---|
| 696 | cin >> true_nn; |
---|
| 697 | valid_dirty = ANNtrue; // validation must be redone |
---|
| 698 | } |
---|
| 699 | //---------------------------------------------------------------- |
---|
| 700 | // seed option |
---|
| 701 | // The seed is reset by setting the global annIdum to the |
---|
| 702 | // negation of the seed value. See rand.cpp. |
---|
| 703 | //---------------------------------------------------------------- |
---|
| 704 | else if (!strcmp(directive,"seed")) { |
---|
| 705 | cin >> annIdum; |
---|
| 706 | annIdum = -annIdum; |
---|
| 707 | } |
---|
| 708 | //---------------------------------------------------------------- |
---|
| 709 | // validate option |
---|
| 710 | //---------------------------------------------------------------- |
---|
| 711 | else if (!strcmp(directive,"validate")) { |
---|
| 712 | cin >> arg; // input argument |
---|
| 713 | if (!strcmp(arg, "on")) { |
---|
| 714 | validate = ANNtrue; |
---|
| 715 | cout << "validate = on " |
---|
| 716 | << "(Warning: this may slow execution time.)\n"; |
---|
| 717 | } |
---|
| 718 | else if (!strcmp(arg, "off")) { |
---|
| 719 | validate = ANNfalse; |
---|
| 720 | } |
---|
| 721 | else { |
---|
| 722 | cerr << "Argument: " << arg << "\n"; |
---|
| 723 | Error("validate argument must be \"on\" or \"off\"", ANNabort); |
---|
| 724 | } |
---|
| 725 | } |
---|
| 726 | //---------------------------------------------------------------- |
---|
| 727 | // distribution option |
---|
| 728 | //---------------------------------------------------------------- |
---|
| 729 | else if (!strcmp(directive,"distribution")) { |
---|
| 730 | cin >> arg; // input name and translate |
---|
| 731 | distr = (Distrib) lookUp(arg, distr_table, N_DISTRIBS); |
---|
| 732 | if (distr >= N_DISTRIBS) { // not something we recognize |
---|
| 733 | cerr << "Distribution: " << arg << "\n"; |
---|
| 734 | Error("Unknown distribution", ANNabort); |
---|
| 735 | } |
---|
| 736 | } |
---|
| 737 | //---------------------------------------------------------------- |
---|
| 738 | // stats option |
---|
| 739 | //---------------------------------------------------------------- |
---|
| 740 | else if (!strcmp(directive,"stats")) { |
---|
| 741 | cin >> arg; // input name and translate |
---|
| 742 | stats = (StatLev) lookUp(arg, stat_table, N_STAT_LEVELS); |
---|
| 743 | if (stats >= N_STAT_LEVELS) { // not something we recognize |
---|
| 744 | cerr << "Stats level: " << arg << "\n"; |
---|
| 745 | Error("Unknown statistics level", ANNabort); |
---|
| 746 | } |
---|
| 747 | if (stats > SILENT) |
---|
| 748 | cout << "stats = " << arg << "\n"; |
---|
| 749 | } |
---|
| 750 | //---------------------------------------------------------------- |
---|
| 751 | // split_rule option |
---|
| 752 | //---------------------------------------------------------------- |
---|
| 753 | else if (!strcmp(directive,"split_rule")) { |
---|
| 754 | cin >> arg; // input split_rule name |
---|
| 755 | split = (ANNsplitRule) lookUp(arg, split_table, N_SPLIT_RULES); |
---|
| 756 | if (split >= N_SPLIT_RULES) { // not something we recognize |
---|
| 757 | cerr << "Splitting rule: " << arg << "\n"; |
---|
| 758 | Error("Unknown splitting rule", ANNabort); |
---|
| 759 | } |
---|
| 760 | } |
---|
| 761 | //---------------------------------------------------------------- |
---|
| 762 | // shrink_rule option |
---|
| 763 | //---------------------------------------------------------------- |
---|
| 764 | else if (!strcmp(directive,"shrink_rule")) { |
---|
| 765 | cin >> arg; // input split_rule name |
---|
| 766 | shrink = (ANNshrinkRule) lookUp(arg, shrink_table, N_SHRINK_RULES); |
---|
| 767 | if (shrink >= N_SHRINK_RULES) { // not something we recognize |
---|
| 768 | cerr << "Shrinking rule: " << arg << "\n"; |
---|
| 769 | Error("Unknown shrinking rule", ANNabort); |
---|
| 770 | } |
---|
| 771 | } |
---|
| 772 | //---------------------------------------------------------------- |
---|
| 773 | // label operation |
---|
| 774 | //---------------------------------------------------------------- |
---|
| 775 | else if (!strcmp(directive,"output_label")) { |
---|
| 776 | cin >> arg; |
---|
| 777 | if (stats > SILENT) |
---|
| 778 | cout << "<" << arg << ">\n"; |
---|
| 779 | } |
---|
| 780 | //---------------------------------------------------------------- |
---|
| 781 | // gen_data_pts operation |
---|
| 782 | //---------------------------------------------------------------- |
---|
| 783 | else if (!strcmp(directive,"gen_data_pts")) { |
---|
| 784 | if (distr == PLANTED) { // planted distribution |
---|
| 785 | Error("Cannot use planted distribution for data points", ANNabort); |
---|
| 786 | } |
---|
| 787 | generatePts( // generate data points |
---|
| 788 | data_pts, // data points |
---|
| 789 | data_size, // data size |
---|
| 790 | DATA, // data points |
---|
| 791 | new_clust); // new clusters flag |
---|
| 792 | valid_dirty = ANNtrue; // validation must be redone |
---|
| 793 | new_clust = ANNfalse; // reset flag |
---|
| 794 | } |
---|
| 795 | //---------------------------------------------------------------- |
---|
| 796 | // gen_query_pts operation |
---|
| 797 | // If the distribution is PLANTED, then the query points |
---|
| 798 | // are planted near the data points (which must already be |
---|
| 799 | // generated). |
---|
| 800 | //---------------------------------------------------------------- |
---|
| 801 | else if (!strcmp(directive,"gen_query_pts")) { |
---|
| 802 | if (distr == PLANTED) { // planted distribution |
---|
| 803 | if (data_pts == NULL) { |
---|
| 804 | Error("Must generate data points before query points for planted distribution", ANNabort); |
---|
| 805 | } |
---|
| 806 | generatePts( // generate query points |
---|
| 807 | query_pts, // point array |
---|
| 808 | query_size, // number of query points |
---|
| 809 | QUERY, // query points |
---|
| 810 | new_clust, // new clusters flag |
---|
| 811 | data_pts, // plant around data pts |
---|
| 812 | data_size); |
---|
| 813 | } |
---|
| 814 | else { // all other distributions |
---|
| 815 | generatePts( // generate query points |
---|
| 816 | query_pts, // point array |
---|
| 817 | query_size, // number of query points |
---|
| 818 | QUERY, // query points |
---|
| 819 | new_clust); // new clusters flag |
---|
| 820 | } |
---|
| 821 | valid_dirty = ANNtrue; // validation must be redone |
---|
| 822 | new_clust = ANNfalse; // reset flag |
---|
| 823 | } |
---|
| 824 | //---------------------------------------------------------------- |
---|
| 825 | // read_data_pts operation |
---|
| 826 | //---------------------------------------------------------------- |
---|
| 827 | else if (!strcmp(directive,"read_data_pts")) { |
---|
| 828 | cin >> arg; // input file name |
---|
| 829 | readPts( |
---|
| 830 | data_pts, // point array |
---|
| 831 | data_size, // number of points |
---|
| 832 | arg, // file name |
---|
| 833 | DATA); // data points |
---|
| 834 | valid_dirty = ANNtrue; // validation must be redone |
---|
| 835 | } |
---|
| 836 | //---------------------------------------------------------------- |
---|
| 837 | // read_query_pts operation |
---|
| 838 | //---------------------------------------------------------------- |
---|
| 839 | else if (!strcmp(directive,"read_query_pts")) { |
---|
| 840 | cin >> arg; // input file name |
---|
| 841 | readPts( |
---|
| 842 | query_pts, // point array |
---|
| 843 | query_size, // number of points |
---|
| 844 | arg, // file name |
---|
| 845 | QUERY); // query points |
---|
| 846 | valid_dirty = ANNtrue; // validation must be redone |
---|
| 847 | } |
---|
| 848 | //---------------------------------------------------------------- |
---|
| 849 | // build_ann operation |
---|
| 850 | // We always invoke the constructor for bd-trees. Note |
---|
| 851 | // that when the shrinking rule is NONE (which is true by |
---|
| 852 | // default), then this constructs a kd-tree. |
---|
| 853 | //---------------------------------------------------------------- |
---|
| 854 | else if (!strcmp(directive,"build_ann")) { |
---|
| 855 | //------------------------------------------------------------ |
---|
| 856 | // Build the tree |
---|
| 857 | //------------------------------------------------------------ |
---|
| 858 | if (the_tree != NULL) { // tree exists already |
---|
| 859 | delete the_tree; // get rid of it |
---|
| 860 | } |
---|
| 861 | clock0 = clock(); // start time |
---|
| 862 | |
---|
| 863 | the_tree = new ANNbd_tree( // build it |
---|
| 864 | data_pts, // the data points |
---|
| 865 | data_size, // number of points |
---|
| 866 | dim, // dimension of space |
---|
| 867 | bucket_size, // maximum bucket size |
---|
| 868 | split, // splitting rule |
---|
| 869 | shrink); // shrinking rule |
---|
| 870 | |
---|
| 871 | //------------------------------------------------------------ |
---|
| 872 | // Print summary |
---|
| 873 | //------------------------------------------------------------ |
---|
| 874 | long prep_time = clock() - clock0; // end of prep time |
---|
| 875 | |
---|
| 876 | if (stats > SILENT) { |
---|
| 877 | cout << "[Build ann-structure:\n"; |
---|
| 878 | cout << " split_rule = " << split_table[split] << "\n"; |
---|
| 879 | cout << " shrink_rule = " << shrink_table[shrink] << "\n"; |
---|
| 880 | cout << " data_size = " << data_size << "\n"; |
---|
| 881 | cout << " dim = " << dim << "\n"; |
---|
| 882 | cout << " bucket_size = " << bucket_size << "\n"; |
---|
| 883 | |
---|
| 884 | if (stats >= EXEC_TIME) { // output processing time |
---|
| 885 | cout << " process_time = " |
---|
| 886 | << double(prep_time)/CLOCKS_PER_SEC << " sec\n"; |
---|
| 887 | } |
---|
| 888 | |
---|
| 889 | if (stats >= PREP_STATS) // output or check tree stats |
---|
| 890 | treeStats(cout, ANNtrue); // print tree stats |
---|
| 891 | else |
---|
| 892 | treeStats(cout, ANNfalse); // check stats |
---|
| 893 | |
---|
| 894 | if (stats >= SHOW_STRUCT) { // print the whole tree |
---|
| 895 | cout << " (Structure Contents:\n"; |
---|
| 896 | the_tree->Print(ANNfalse, cout); |
---|
| 897 | cout << " )\n"; |
---|
| 898 | } |
---|
| 899 | cout << "]\n"; |
---|
| 900 | } |
---|
| 901 | } |
---|
| 902 | //---------------------------------------------------------------- |
---|
| 903 | // dump operation |
---|
| 904 | //---------------------------------------------------------------- |
---|
| 905 | else if (!strcmp(directive,"dump")) { |
---|
| 906 | cin >> arg; // input file name |
---|
| 907 | if (the_tree == NULL) { // no tree |
---|
| 908 | Error("Cannot dump. No tree has been built yet", ANNwarn); |
---|
| 909 | } |
---|
| 910 | else { // there is a tree |
---|
| 911 | // try to open file |
---|
| 912 | ofstream out_dump_file(arg); |
---|
| 913 | if (!out_dump_file) { |
---|
| 914 | cerr << "File name: " << arg << "\n"; |
---|
| 915 | Error("Cannot open dump file", ANNabort); |
---|
| 916 | } |
---|
| 917 | // dump the tree and points |
---|
| 918 | the_tree->Dump(ANNtrue, out_dump_file); |
---|
| 919 | if (stats > SILENT) { |
---|
| 920 | cout << "(Tree has been dumped to file " << arg << ")\n"; |
---|
| 921 | } |
---|
| 922 | } |
---|
| 923 | } |
---|
| 924 | //---------------------------------------------------------------- |
---|
| 925 | // load operation |
---|
| 926 | // Since this not only loads a tree, but loads a new set |
---|
| 927 | // of data points. |
---|
| 928 | //---------------------------------------------------------------- |
---|
| 929 | else if (!strcmp(directive,"load")) { |
---|
| 930 | cin >> arg; // input file name |
---|
| 931 | if (the_tree != NULL) { // tree exists already |
---|
| 932 | delete the_tree; // get rid of it |
---|
| 933 | } |
---|
| 934 | if (data_pts != NULL) { // data points exist already |
---|
| 935 | delete data_pts; // get rid of them |
---|
| 936 | } |
---|
| 937 | |
---|
| 938 | ifstream in_dump_file(arg); // try to open file |
---|
| 939 | if (!in_dump_file) { |
---|
| 940 | cerr << "File name: " << arg << "\n"; |
---|
| 941 | Error("Cannot open file for loading", ANNabort); |
---|
| 942 | } |
---|
| 943 | // build tree by loading |
---|
| 944 | the_tree = new ANNbd_tree(in_dump_file); |
---|
| 945 | |
---|
| 946 | dim = the_tree->theDim(); // new dimension |
---|
| 947 | data_size = the_tree->nPoints(); // number of points |
---|
| 948 | data_pts = the_tree->thePoints(); // new points |
---|
| 949 | |
---|
| 950 | valid_dirty = ANNtrue; // validation must be redone |
---|
| 951 | |
---|
| 952 | if (stats > SILENT) { |
---|
| 953 | cout << "(Tree has been loaded from file " << arg << ")\n"; |
---|
| 954 | } |
---|
| 955 | if (stats >= SHOW_STRUCT) { // print the tree |
---|
| 956 | cout << " (Structure Contents:\n"; |
---|
| 957 | the_tree->Print(ANNfalse, cout); |
---|
| 958 | cout << " )\n"; |
---|
| 959 | } |
---|
| 960 | } |
---|
| 961 | //---------------------------------------------------------------- |
---|
| 962 | // run_queries operation |
---|
| 963 | // This section does all the query processing. It consists |
---|
| 964 | // of the following subsections: |
---|
| 965 | // |
---|
| 966 | // ** input the argument (standard or priority) and output |
---|
| 967 | // the header describing the essential information. |
---|
| 968 | // ** allocate space for the results to be stored. |
---|
| 969 | // ** run the queries by invoking the appropriate search |
---|
| 970 | // procedure on the query points. Print nearest neighbor |
---|
| 971 | // if requested. |
---|
| 972 | // ** print final summaries |
---|
| 973 | // |
---|
| 974 | // The approach for processing multiple nearest neighbors is |
---|
| 975 | // pretty crude. We allocate an array whose size is the |
---|
| 976 | // product of the total number of queries times the number of |
---|
| 977 | // nearest neighbors (k), and then use each k consecutive |
---|
| 978 | // entries to store the results of each query. |
---|
| 979 | //---------------------------------------------------------------- |
---|
| 980 | else if (!strcmp(directive,"run_queries")) { |
---|
| 981 | |
---|
| 982 | //------------------------------------------------------------ |
---|
| 983 | // Input arguments and print summary |
---|
| 984 | //------------------------------------------------------------ |
---|
| 985 | enum {STANDARD, PRIORITY} method; |
---|
| 986 | |
---|
| 987 | cin >> arg; // input argument |
---|
| 988 | if (!strcmp(arg, "standard")) { |
---|
| 989 | method = STANDARD; |
---|
| 990 | } |
---|
| 991 | else if (!strcmp(arg, "priority")) { |
---|
| 992 | method = PRIORITY; |
---|
| 993 | } |
---|
| 994 | else { |
---|
| 995 | cerr << "Search type: " << arg << "\n"; |
---|
| 996 | Error("Search type must be \"standard\" or \"priority\"", |
---|
| 997 | ANNabort); |
---|
| 998 | } |
---|
| 999 | if (data_pts == NULL || query_pts == NULL) { |
---|
| 1000 | Error("Either data set and query set not constructed", ANNabort); |
---|
| 1001 | } |
---|
| 1002 | if (the_tree == NULL) { |
---|
| 1003 | Error("No search tree built.", ANNabort); |
---|
| 1004 | } |
---|
| 1005 | |
---|
| 1006 | //------------------------------------------------------------ |
---|
| 1007 | // Set up everything |
---|
| 1008 | //------------------------------------------------------------ |
---|
| 1009 | |
---|
| 1010 | #ifdef ANN_PERF // performance only |
---|
| 1011 | annResetStats(data_size); // reset statistics |
---|
| 1012 | #endif |
---|
| 1013 | |
---|
| 1014 | clock0 = clock(); // start time |
---|
| 1015 | // deallocate existing storage |
---|
| 1016 | if (apx_nn_idx != NULL) delete [] apx_nn_idx; |
---|
| 1017 | if (apx_dists != NULL) delete [] apx_dists; |
---|
| 1018 | if (apx_pts_in_range != NULL) delete [] apx_pts_in_range; |
---|
| 1019 | // allocate apx answer storage |
---|
| 1020 | apx_nn_idx = new ANNidx[near_neigh*query_size]; |
---|
| 1021 | apx_dists = new ANNdist[near_neigh*query_size]; |
---|
| 1022 | apx_pts_in_range = new int[query_size]; |
---|
| 1023 | |
---|
| 1024 | annMaxPtsVisit(max_pts_visit); // set max points to visit |
---|
| 1025 | |
---|
| 1026 | //------------------------------------------------------------ |
---|
| 1027 | // Run the queries |
---|
| 1028 | //------------------------------------------------------------ |
---|
| 1029 | // pointers for current query |
---|
| 1030 | ANNidxArray curr_nn_idx = apx_nn_idx; |
---|
| 1031 | ANNdistArray curr_dists = apx_dists; |
---|
| 1032 | |
---|
| 1033 | for (int i = 0; i < query_size; i++) { |
---|
| 1034 | #ifdef ANN_PERF |
---|
| 1035 | annResetCounts(); // reset counters |
---|
| 1036 | #endif |
---|
| 1037 | apx_pts_in_range[i] = 0; |
---|
| 1038 | |
---|
| 1039 | if (radius_bound == 0) { // no radius bound |
---|
| 1040 | if (method == STANDARD) { |
---|
| 1041 | the_tree->annkSearch( |
---|
| 1042 | query_pts[i], // query point |
---|
| 1043 | near_neigh, // number of near neighbors |
---|
| 1044 | curr_nn_idx, // nearest neighbors (returned) |
---|
| 1045 | curr_dists, // distance (returned) |
---|
| 1046 | epsilon); // error bound |
---|
| 1047 | } |
---|
| 1048 | else if (method == PRIORITY) { |
---|
| 1049 | the_tree->annkPriSearch( |
---|
| 1050 | query_pts[i], // query point |
---|
| 1051 | near_neigh, // number of near neighbors |
---|
| 1052 | curr_nn_idx, // nearest neighbors (returned) |
---|
| 1053 | curr_dists, // distance (returned) |
---|
| 1054 | epsilon); // error bound |
---|
| 1055 | } |
---|
| 1056 | else { |
---|
| 1057 | Error("Internal error - invalid method", ANNabort); |
---|
| 1058 | } |
---|
| 1059 | } |
---|
| 1060 | else { // use radius bound |
---|
| 1061 | if (method != STANDARD) { |
---|
| 1062 | Error("A nonzero radius bound assumes standard search", |
---|
| 1063 | ANNwarn); |
---|
| 1064 | } |
---|
| 1065 | apx_pts_in_range[i] = the_tree->annkFRSearch( |
---|
| 1066 | query_pts[i], // query point |
---|
| 1067 | ANN_POW(radius_bound), // squared radius search bound |
---|
| 1068 | near_neigh, // number of near neighbors |
---|
| 1069 | curr_nn_idx, // nearest neighbors (returned) |
---|
| 1070 | curr_dists, // distance (returned) |
---|
| 1071 | epsilon); // error bound |
---|
| 1072 | } |
---|
| 1073 | curr_nn_idx += near_neigh; // increment current pointers |
---|
| 1074 | curr_dists += near_neigh; |
---|
| 1075 | |
---|
| 1076 | #ifdef ANN_PERF |
---|
| 1077 | annUpdateStats(); // update stats |
---|
| 1078 | #endif |
---|
| 1079 | } |
---|
| 1080 | |
---|
| 1081 | long query_time = clock() - clock0; // end of query time |
---|
| 1082 | |
---|
| 1083 | if (validate) { // validation requested |
---|
| 1084 | if (valid_dirty) getTrueNN(); // get true near neighbors |
---|
| 1085 | doValidation(); // validate |
---|
| 1086 | } |
---|
| 1087 | |
---|
| 1088 | //------------------------------------------------------------ |
---|
| 1089 | // Print summaries |
---|
| 1090 | //------------------------------------------------------------ |
---|
| 1091 | |
---|
| 1092 | if (stats > SILENT) { |
---|
| 1093 | cout << "[Run Queries:\n"; |
---|
| 1094 | cout << " query_size = " << query_size << "\n"; |
---|
| 1095 | cout << " dim = " << dim << "\n"; |
---|
| 1096 | cout << " search_method = " << arg << "\n"; |
---|
| 1097 | cout << " epsilon = " << epsilon << "\n"; |
---|
| 1098 | cout << " near_neigh = " << near_neigh << "\n"; |
---|
| 1099 | if (max_pts_visit != 0) |
---|
| 1100 | cout << " max_pts_visit = " << max_pts_visit << "\n"; |
---|
| 1101 | if (radius_bound != 0) |
---|
| 1102 | cout << " radius_bound = " << radius_bound << "\n"; |
---|
| 1103 | if (validate) |
---|
| 1104 | cout << " true_nn = " << true_nn << "\n"; |
---|
| 1105 | |
---|
| 1106 | if (stats >= EXEC_TIME) { // print exec time summary |
---|
| 1107 | cout << " query_time = " << |
---|
| 1108 | double(query_time)/(query_size*CLOCKS_PER_SEC) |
---|
| 1109 | << " sec/query"; |
---|
| 1110 | #ifdef ANN_PERF |
---|
| 1111 | cout << " (biased by perf measurements)"; |
---|
| 1112 | #endif |
---|
| 1113 | cout << "\n"; |
---|
| 1114 | } |
---|
| 1115 | |
---|
| 1116 | if (stats >= QUERY_STATS) { // output performance stats |
---|
| 1117 | #ifdef ANN_PERF |
---|
| 1118 | cout.flush(); |
---|
| 1119 | annPrintStats(validate); |
---|
| 1120 | #else |
---|
| 1121 | cout << " (Performance statistics unavailable.)\n"; |
---|
| 1122 | #endif |
---|
| 1123 | } |
---|
| 1124 | |
---|
| 1125 | if (stats >= QUERY_RES) { // output results |
---|
| 1126 | cout << " (Query Results:\n"; |
---|
| 1127 | cout << " Pt\tANN\tDist\n"; |
---|
| 1128 | curr_nn_idx = apx_nn_idx; // subarray pointers |
---|
| 1129 | curr_dists = apx_dists; |
---|
| 1130 | // output nearest neighbors |
---|
| 1131 | for (int i = 0; i < query_size; i++) { |
---|
| 1132 | cout << " " << setw(4) << i; |
---|
| 1133 | for (int j = 0; j < near_neigh; j++) { |
---|
| 1134 | // exit if no more neighbors |
---|
| 1135 | if (curr_nn_idx[j] == ANN_NULL_IDX) { |
---|
| 1136 | cout << "\t[no other pts in radius bound]\n"; |
---|
| 1137 | break; |
---|
| 1138 | } |
---|
| 1139 | else { // output point info |
---|
| 1140 | cout << "\t" << curr_nn_idx[j] |
---|
| 1141 | << "\t" << ANN_ROOT(curr_dists[j]) |
---|
| 1142 | << "\n"; |
---|
| 1143 | } |
---|
| 1144 | } |
---|
| 1145 | // output range count |
---|
| 1146 | if (radius_bound != 0) { |
---|
| 1147 | cout << " pts_in_radius_bound = " |
---|
| 1148 | << apx_pts_in_range[i] << "\n"; |
---|
| 1149 | } |
---|
| 1150 | // increment subarray pointers |
---|
| 1151 | curr_nn_idx += near_neigh; |
---|
| 1152 | curr_dists += near_neigh; |
---|
| 1153 | } |
---|
| 1154 | cout << " )\n"; |
---|
| 1155 | } |
---|
| 1156 | cout << "]\n"; |
---|
| 1157 | } |
---|
| 1158 | } |
---|
| 1159 | //---------------------------------------------------------------- |
---|
| 1160 | // Unknown directive |
---|
| 1161 | //---------------------------------------------------------------- |
---|
| 1162 | else { |
---|
| 1163 | cerr << "Directive: " << directive << "\n"; |
---|
| 1164 | Error("Unknown directive", ANNabort); |
---|
| 1165 | } |
---|
| 1166 | } |
---|
| 1167 | //-------------------------------------------------------------------- |
---|
| 1168 | // End of input loop (deallocate stuff that was allocated) |
---|
| 1169 | //-------------------------------------------------------------------- |
---|
| 1170 | if (the_tree != NULL) delete the_tree; |
---|
| 1171 | if (data_pts != NULL) annDeallocPts(data_pts); |
---|
| 1172 | if (query_pts != NULL) annDeallocPts(query_pts); |
---|
| 1173 | if (apx_nn_idx != NULL) delete [] apx_nn_idx; |
---|
| 1174 | if (apx_dists != NULL) delete [] apx_dists; |
---|
| 1175 | if (apx_pts_in_range != NULL) delete [] apx_pts_in_range; |
---|
| 1176 | |
---|
| 1177 | annClose(); // close ANN |
---|
| 1178 | |
---|
| 1179 | return EXIT_SUCCESS; |
---|
| 1180 | } |
---|
| 1181 | |
---|
| 1182 | //------------------------------------------------------------------------ |
---|
| 1183 | // generatePts - call appropriate routine to generate points of a |
---|
| 1184 | // given distribution. |
---|
| 1185 | //------------------------------------------------------------------------ |
---|
| 1186 | |
---|
| 1187 | void generatePts( |
---|
| 1188 | ANNpointArray &pa, // point array (returned) |
---|
| 1189 | int n, // number of points to generate |
---|
| 1190 | PtType type, // point type |
---|
| 1191 | ANNbool new_clust, // new cluster centers desired? |
---|
| 1192 | ANNpointArray src, // source array (if distr=PLANTED) |
---|
| 1193 | int n_src) // source size (if distr=PLANTED) |
---|
| 1194 | { |
---|
| 1195 | if (pa != NULL) annDeallocPts(pa); // get rid of any old points |
---|
| 1196 | pa = annAllocPts(n, dim); // allocate point storage |
---|
| 1197 | |
---|
| 1198 | switch (distr) { |
---|
| 1199 | case UNIFORM: // uniform over cube [-1,1]^d. |
---|
| 1200 | annUniformPts(pa, n, dim); |
---|
| 1201 | break; |
---|
| 1202 | case GAUSS: // Gaussian with mean 0 |
---|
| 1203 | annGaussPts(pa, n, dim, std_dev); |
---|
| 1204 | break; |
---|
| 1205 | case LAPLACE: // Laplacian, mean 0 and var 1 |
---|
| 1206 | annLaplacePts(pa, n, dim); |
---|
| 1207 | break; |
---|
| 1208 | case CO_GAUSS: // correlated Gaussian |
---|
| 1209 | annCoGaussPts(pa, n, dim, corr_coef); |
---|
| 1210 | break; |
---|
| 1211 | case CO_LAPLACE: // correlated Laplacian |
---|
| 1212 | annCoLaplacePts(pa, n, dim, corr_coef); |
---|
| 1213 | break; |
---|
| 1214 | case CLUS_GAUSS: // clustered Gaussian |
---|
| 1215 | annClusGaussPts(pa, n, dim, n_color, new_clust, std_dev); |
---|
| 1216 | break; |
---|
| 1217 | case CLUS_ORTH_FLATS: // clustered on orthog flats |
---|
| 1218 | annClusOrthFlats(pa, n, dim, n_color, new_clust, std_dev, max_dim); |
---|
| 1219 | break; |
---|
| 1220 | case CLUS_ELLIPSOIDS: // clustered ellipsoids |
---|
| 1221 | annClusEllipsoids(pa, n, dim, n_color, new_clust, std_dev, |
---|
| 1222 | std_dev_lo, std_dev_hi, max_dim); |
---|
| 1223 | break; |
---|
| 1224 | case PLANTED: // planted distribution |
---|
| 1225 | annPlanted(pa, n, dim, src, n_src, std_dev); |
---|
| 1226 | break; |
---|
| 1227 | default: |
---|
| 1228 | Error("INTERNAL ERROR: Unknown distribution", ANNabort); |
---|
| 1229 | break; |
---|
| 1230 | } |
---|
| 1231 | |
---|
| 1232 | if (stats > SILENT) { |
---|
| 1233 | if(type == DATA) cout << "[Generating Data Points:\n"; |
---|
| 1234 | else cout << "[Generating Query Points:\n"; |
---|
| 1235 | cout << " number = " << n << "\n"; |
---|
| 1236 | cout << " dim = " << dim << "\n"; |
---|
| 1237 | cout << " distribution = " << distr_table[distr] << "\n"; |
---|
| 1238 | if (annIdum < 0) |
---|
| 1239 | cout << " seed = " << annIdum << "\n"; |
---|
| 1240 | if (distr == GAUSS || distr == CLUS_GAUSS |
---|
| 1241 | || distr == CLUS_ORTH_FLATS) |
---|
| 1242 | cout << " std_dev = " << std_dev << "\n"; |
---|
| 1243 | if (distr == CLUS_ELLIPSOIDS) { |
---|
| 1244 | cout << " std_dev = " << std_dev << " (small) \n"; |
---|
| 1245 | cout << " std_dev_lo = " << std_dev_lo << "\n"; |
---|
| 1246 | cout << " std_dev_hi = " << std_dev_hi << "\n"; |
---|
| 1247 | } |
---|
| 1248 | if (distr == CO_GAUSS || distr == CO_LAPLACE) |
---|
| 1249 | cout << " corr_coef = " << corr_coef << "\n"; |
---|
| 1250 | if (distr == CLUS_GAUSS || distr == CLUS_ORTH_FLATS |
---|
| 1251 | || distr == CLUS_ELLIPSOIDS) { |
---|
| 1252 | cout << " colors = " << n_color << "\n"; |
---|
| 1253 | if (new_clust) |
---|
| 1254 | cout << " (cluster centers regenerated)\n"; |
---|
| 1255 | } |
---|
| 1256 | if (distr == CLUS_ORTH_FLATS || distr == CLUS_ELLIPSOIDS) { |
---|
| 1257 | cout << " max_dim = " << max_dim << "\n"; |
---|
| 1258 | } |
---|
| 1259 | } |
---|
| 1260 | // want to see points? |
---|
| 1261 | if ((type == DATA && stats >= SHOW_PTS) || |
---|
| 1262 | (type == QUERY && stats >= QUERY_RES)) { |
---|
| 1263 | if(type == DATA) cout << "(Data Points:\n"; |
---|
| 1264 | else cout << "(Query Points:\n"; |
---|
| 1265 | for (int i = 0; i < n; i++) { |
---|
| 1266 | cout << " " << setw(4) << i << "\t"; |
---|
| 1267 | printPoint(pa[i], dim); |
---|
| 1268 | cout << "\n"; |
---|
| 1269 | } |
---|
| 1270 | cout << " )\n"; |
---|
| 1271 | } |
---|
| 1272 | cout << "]\n"; |
---|
| 1273 | } |
---|
| 1274 | |
---|
| 1275 | //------------------------------------------------------------------------ |
---|
| 1276 | // readPts - read a collection of data or query points. |
---|
| 1277 | //------------------------------------------------------------------------ |
---|
| 1278 | |
---|
| 1279 | void readPts( |
---|
| 1280 | ANNpointArray &pa, // point array (returned) |
---|
| 1281 | int &n, // number of points |
---|
| 1282 | char *file_nm, // file name |
---|
| 1283 | PtType type) // point type (DATA, QUERY) |
---|
| 1284 | { |
---|
| 1285 | int i; |
---|
| 1286 | //-------------------------------------------------------------------- |
---|
| 1287 | // Open input file and read points |
---|
| 1288 | //-------------------------------------------------------------------- |
---|
| 1289 | ifstream in_file(file_nm); // try to open data file |
---|
| 1290 | if (!in_file) { |
---|
| 1291 | cerr << "File name: " << file_nm << "\n"; |
---|
| 1292 | Error("Cannot open input data/query file", ANNabort); |
---|
| 1293 | } |
---|
| 1294 | // allocate storage for points |
---|
| 1295 | if (pa != NULL) annDeallocPts(pa); // get rid of old points |
---|
| 1296 | pa = annAllocPts(n, dim); |
---|
| 1297 | |
---|
| 1298 | for (i = 0; i < n; i++) { // read the data |
---|
| 1299 | if (!(in_file >> pa[i][0])) break; |
---|
| 1300 | for (int d = 1; d < dim; d++) { |
---|
| 1301 | in_file >> pa[i][d]; |
---|
| 1302 | } |
---|
| 1303 | } |
---|
| 1304 | |
---|
| 1305 | char ignore_me; // character for EOF test |
---|
| 1306 | in_file >> ignore_me; // try to get one more character |
---|
| 1307 | if (!in_file.eof()) { // exhausted space before eof |
---|
| 1308 | if (type == DATA) |
---|
| 1309 | Error("`data_size' too small. Input file truncated.", ANNwarn); |
---|
| 1310 | else |
---|
| 1311 | Error("`query_size' too small. Input file truncated.", ANNwarn); |
---|
| 1312 | } |
---|
| 1313 | n = i; // number of points read |
---|
| 1314 | |
---|
| 1315 | //-------------------------------------------------------------------- |
---|
| 1316 | // Print summary |
---|
| 1317 | //-------------------------------------------------------------------- |
---|
| 1318 | if (stats > SILENT) { |
---|
| 1319 | if (type == DATA) { |
---|
| 1320 | cout << "[Read Data Points:\n"; |
---|
| 1321 | cout << " data_size = " << n << "\n"; |
---|
| 1322 | } |
---|
| 1323 | else { |
---|
| 1324 | cout << "[Read Query Points:\n"; |
---|
| 1325 | cout << " query_size = " << n << "\n"; |
---|
| 1326 | } |
---|
| 1327 | cout << " file_name = " << file_nm << "\n"; |
---|
| 1328 | cout << " dim = " << dim << "\n"; |
---|
| 1329 | // print if results requested |
---|
| 1330 | if ((type == DATA && stats >= SHOW_PTS) || |
---|
| 1331 | (type == QUERY && stats >= QUERY_RES)) { |
---|
| 1332 | cout << " (Points:\n"; |
---|
| 1333 | for (i = 0; i < n; i++) { |
---|
| 1334 | cout << " " << i << "\t"; |
---|
| 1335 | printPoint(pa[i], dim); |
---|
| 1336 | cout << "\n"; |
---|
| 1337 | } |
---|
| 1338 | cout << " )\n"; |
---|
| 1339 | } |
---|
| 1340 | cout << "]\n"; |
---|
| 1341 | } |
---|
| 1342 | } |
---|
| 1343 | |
---|
| 1344 | //------------------------------------------------------------------------ |
---|
| 1345 | // getTrueNN |
---|
| 1346 | // Computes the true nearest neighbors. For purposes of validation, |
---|
| 1347 | // this intentionally done in a rather dumb (but safe way), by |
---|
| 1348 | // invoking the brute-force search. |
---|
| 1349 | // |
---|
| 1350 | // The number of true nearest neighbors is somewhat larger than |
---|
| 1351 | // the number of nearest neighbors. This is so that the validation |
---|
| 1352 | // can determine the expected difference in element ranks. |
---|
| 1353 | // |
---|
| 1354 | // This procedure is invoked just prior to running queries. Since |
---|
| 1355 | // the operation takes a long time, it is performed only if needed. |
---|
| 1356 | // In particular, once generated, it will be regenerated only if |
---|
| 1357 | // new query or data points are generated, or if the requested number |
---|
| 1358 | // of true near neighbors or approximate near neighbors has changed. |
---|
| 1359 | // |
---|
| 1360 | // To validate fixed-radius searching, we compute two counts, one |
---|
| 1361 | // with the original query radius (trueSqRadius) and the other with |
---|
| 1362 | // a radius shrunken by the error factor (minSqradius). We then |
---|
| 1363 | // check that the count of points inside the approximate range is |
---|
| 1364 | // between these two bounds. Because fixed-radius search is |
---|
| 1365 | // allowed to ignore points within the shrunken radius, we only |
---|
| 1366 | // compute exact neighbors within this smaller distance (for we |
---|
| 1367 | // cannot guarantee that we will even visit the other points). |
---|
| 1368 | //------------------------------------------------------------------------ |
---|
| 1369 | |
---|
| 1370 | void getTrueNN() // compute true nearest neighbors |
---|
| 1371 | { |
---|
| 1372 | if (stats > SILENT) { |
---|
| 1373 | cout << "(Computing true nearest neighbors for validation. This may take time.)\n"; |
---|
| 1374 | } |
---|
| 1375 | // deallocate existing storage |
---|
| 1376 | if (true_nn_idx != NULL) delete [] true_nn_idx; |
---|
| 1377 | if (true_dists != NULL) delete [] true_dists; |
---|
| 1378 | if (min_pts_in_range != NULL) delete [] min_pts_in_range; |
---|
| 1379 | if (max_pts_in_range != NULL) delete [] max_pts_in_range; |
---|
| 1380 | |
---|
| 1381 | if (true_nn > data_size) { // can't get more nn than points |
---|
| 1382 | true_nn = data_size; |
---|
| 1383 | } |
---|
| 1384 | |
---|
| 1385 | // allocate true answer storage |
---|
| 1386 | true_nn_idx = new ANNidx[true_nn*query_size]; |
---|
| 1387 | true_dists = new ANNdist[true_nn*query_size]; |
---|
| 1388 | min_pts_in_range = new int[query_size]; |
---|
| 1389 | max_pts_in_range = new int[query_size]; |
---|
| 1390 | |
---|
| 1391 | ANNidxArray curr_nn_idx = true_nn_idx; // current locations in arrays |
---|
| 1392 | ANNdistArray curr_dists = true_dists; |
---|
| 1393 | |
---|
| 1394 | // allocate search structure |
---|
| 1395 | ANNbruteForce *the_brute = new ANNbruteForce(data_pts, data_size, dim); |
---|
| 1396 | // compute nearest neighbors |
---|
| 1397 | for (int i = 0; i < query_size; i++) { |
---|
| 1398 | if (radius_bound == 0) { // standard kNN search |
---|
| 1399 | the_brute->annkSearch( // compute true near neighbors |
---|
| 1400 | query_pts[i], // query point |
---|
| 1401 | true_nn, // number of nearest neighbors |
---|
| 1402 | curr_nn_idx, // where to put indices |
---|
| 1403 | curr_dists); // where to put distances |
---|
| 1404 | } |
---|
| 1405 | else { // fixed radius kNN search |
---|
| 1406 | // search radii limits |
---|
| 1407 | ANNdist trueSqRadius = ANN_POW(radius_bound); |
---|
| 1408 | ANNdist minSqRadius = ANN_POW(radius_bound / (1+epsilon)); |
---|
| 1409 | min_pts_in_range[i] = the_brute->annkFRSearch( |
---|
| 1410 | query_pts[i], // query point |
---|
| 1411 | minSqRadius, // shrunken search radius |
---|
| 1412 | true_nn, // number of near neighbors |
---|
| 1413 | curr_nn_idx, // nearest neighbors (returned) |
---|
| 1414 | curr_dists); // distance (returned) |
---|
| 1415 | max_pts_in_range[i] = the_brute->annkFRSearch( |
---|
| 1416 | query_pts[i], // query point |
---|
| 1417 | trueSqRadius, // true search radius |
---|
| 1418 | 0, NULL, NULL); // (ignore kNN info) |
---|
| 1419 | } |
---|
| 1420 | curr_nn_idx += true_nn; // increment nn index pointer |
---|
| 1421 | curr_dists += true_nn; // increment nn dist pointer |
---|
| 1422 | } |
---|
| 1423 | delete the_brute; // delete brute-force struct |
---|
| 1424 | valid_dirty = ANNfalse; // validation good for now |
---|
| 1425 | } |
---|
| 1426 | |
---|
| 1427 | //------------------------------------------------------------------------ |
---|
| 1428 | // doValidation |
---|
| 1429 | // Compares the approximate answers to the k-nearest neighbors |
---|
| 1430 | // against the true nearest neighbors (computed earlier). It is |
---|
| 1431 | // assumed that the true nearest neighbors and indices have been |
---|
| 1432 | // computed earlier. |
---|
| 1433 | // |
---|
| 1434 | // First, we check that all the results are within their allowed |
---|
| 1435 | // limits, and generate an internal error, if not. For the sake of |
---|
| 1436 | // performance evaluation, we also compute the following two |
---|
| 1437 | // quantities for nearest neighbors: |
---|
| 1438 | // |
---|
| 1439 | // Average Error |
---|
| 1440 | // ------------- |
---|
| 1441 | // The relative error between the distance to a reported nearest |
---|
| 1442 | // neighbor and the true nearest neighbor (of the same rank), |
---|
| 1443 | // |
---|
| 1444 | // Rank Error |
---|
| 1445 | // ---------- |
---|
| 1446 | // The difference in rank between the reported nearest neighbor and |
---|
| 1447 | // its position (if any) among the true nearest neighbors. If we |
---|
| 1448 | // cannot find this point among the true nearest neighbors, then |
---|
| 1449 | // it assumed that the rank of the true nearest neighbor is true_nn+1. |
---|
| 1450 | // |
---|
| 1451 | // Because of the possibility of duplicate distances, this is computed |
---|
| 1452 | // as follows. For the j-th reported nearest neighbor, we count the |
---|
| 1453 | // number of true nearest neighbors that are at least this close. Let |
---|
| 1454 | // this be rnk. Then the rank error is max(0, j-rnk). (In the code |
---|
| 1455 | // below, j is an array index and so the first item is 0, not 1. Thus |
---|
| 1456 | // we take max(0, j+1-rnk) instead.) |
---|
| 1457 | // |
---|
| 1458 | // For the results of fixed-radious range count, we verify that the |
---|
| 1459 | // reported number of points in the range lies between the actual |
---|
| 1460 | // number of points in the shrunken and the true search radius. |
---|
| 1461 | //------------------------------------------------------------------------ |
---|
| 1462 | |
---|
| 1463 | void doValidation() // perform validation |
---|
| 1464 | { |
---|
| 1465 | int* curr_apx_idx = apx_nn_idx; // approx index pointer |
---|
| 1466 | ANNdistArray curr_apx_dst = apx_dists; // approx distance pointer |
---|
| 1467 | int* curr_tru_idx = true_nn_idx; // true index pointer |
---|
| 1468 | ANNdistArray curr_tru_dst = true_dists; // true distance pointer |
---|
| 1469 | int i, j; |
---|
| 1470 | |
---|
| 1471 | if (true_nn < near_neigh) { |
---|
| 1472 | Error("Cannot validate with fewer true near neighbors than actual", ANNabort); |
---|
| 1473 | } |
---|
| 1474 | |
---|
| 1475 | for (i = 0; i < query_size; i++) { // validate each query |
---|
| 1476 | //---------------------------------------------------------------- |
---|
| 1477 | // Compute result errors |
---|
| 1478 | // In fixed radius search it is possible that not all k |
---|
| 1479 | // nearest neighbors were computed. Because the true |
---|
| 1480 | // results are computed over the shrunken radius, we should |
---|
| 1481 | // have at least as many true nearest neighbors as |
---|
| 1482 | // approximate nearest neighbors. (If not, an infinite |
---|
| 1483 | // error will be generated, and so an internal error will |
---|
| 1484 | // will be generated. |
---|
| 1485 | // |
---|
| 1486 | // Because nearest neighbors are sorted in increasing order |
---|
| 1487 | // of distance, as soon as we see a null index, we can |
---|
| 1488 | // terminate the distance checking. The error in the |
---|
| 1489 | // result should not exceed epsilon. However, if |
---|
| 1490 | // max_pts_visit is nonzero (meaning that the search is |
---|
| 1491 | // terminated early) this might happen. |
---|
| 1492 | //---------------------------------------------------------------- |
---|
| 1493 | for (j = 0; j < near_neigh; j++) { |
---|
| 1494 | if (curr_tru_idx[j] == ANN_NULL_IDX)// no more true neighbors? |
---|
| 1495 | break; |
---|
| 1496 | // true i-th smallest distance |
---|
| 1497 | double true_dist = ANN_ROOT(curr_tru_dst[j]); |
---|
| 1498 | // reported i-th smallest |
---|
| 1499 | double rept_dist = ANN_ROOT(curr_apx_dst[j]); |
---|
| 1500 | // better than optimum? |
---|
| 1501 | if (rept_dist < true_dist*(1-ERR)) { |
---|
| 1502 | Error("INTERNAL ERROR: True nearest neighbor incorrect", |
---|
| 1503 | ANNabort); |
---|
| 1504 | } |
---|
| 1505 | |
---|
| 1506 | double resultErr; // result error |
---|
| 1507 | if (true_dist == 0.0) { // let's not divide by zero |
---|
| 1508 | if (rept_dist != 0.0) resultErr = ANN_DBL_MAX; |
---|
| 1509 | else resultErr = 0.0; |
---|
| 1510 | } |
---|
| 1511 | else { |
---|
| 1512 | resultErr = (rept_dist - true_dist) / ((double) true_dist); |
---|
| 1513 | } |
---|
| 1514 | |
---|
| 1515 | if (resultErr > epsilon && max_pts_visit == 0) { |
---|
| 1516 | Error("INTERNAL ERROR: Actual error exceeds epsilon", |
---|
| 1517 | ANNabort); |
---|
| 1518 | } |
---|
| 1519 | #ifdef ANN_PERF |
---|
| 1520 | ann_average_err += resultErr; // update statistics error |
---|
| 1521 | #endif |
---|
| 1522 | } |
---|
| 1523 | //-------------------------------------------------------------------- |
---|
| 1524 | // Compute rank errors (only needed for perf measurements) |
---|
| 1525 | //-------------------------------------------------------------------- |
---|
| 1526 | #ifdef ANN_PERF |
---|
| 1527 | for (j = 0; j < near_neigh; j++) { |
---|
| 1528 | if (curr_tru_idx[i] == ANN_NULL_IDX) // no more true neighbors? |
---|
| 1529 | break; |
---|
| 1530 | |
---|
| 1531 | double rnkErr = 0.0; // rank error |
---|
| 1532 | // reported j-th distance |
---|
| 1533 | ANNdist rept_dist = curr_apx_dst[j]; |
---|
| 1534 | int rnk = 0; // compute rank of this item |
---|
| 1535 | while (rnk < true_nn && curr_tru_dst[rnk] <= rept_dist) |
---|
| 1536 | rnk++; |
---|
| 1537 | if (j+1-rnk > 0) rnkErr = (double) (j+1-rnk); |
---|
| 1538 | ann_rank_err += rnkErr; // update average rank error |
---|
| 1539 | } |
---|
| 1540 | #endif |
---|
| 1541 | //---------------------------------------------------------------- |
---|
| 1542 | // Check range counts from fixed-radius query |
---|
| 1543 | //---------------------------------------------------------------- |
---|
| 1544 | if (radius_bound != 0) { // fixed-radius search |
---|
| 1545 | if (apx_pts_in_range[i] < min_pts_in_range[i] || |
---|
| 1546 | apx_pts_in_range[i] > max_pts_in_range[i]) |
---|
| 1547 | Error("INTERNAL ERROR: Invalid fixed-radius range count", |
---|
| 1548 | ANNabort); |
---|
| 1549 | } |
---|
| 1550 | |
---|
| 1551 | curr_apx_idx += near_neigh; |
---|
| 1552 | curr_apx_dst += near_neigh; |
---|
| 1553 | curr_tru_idx += true_nn; // increment current pointers |
---|
| 1554 | curr_tru_dst += true_nn; |
---|
| 1555 | } |
---|
| 1556 | } |
---|
| 1557 | |
---|
| 1558 | //---------------------------------------------------------------------- |
---|
| 1559 | // treeStats |
---|
| 1560 | // Computes a number of statistics related to kd_trees and |
---|
| 1561 | // bd_trees. These statistics are printed if in verbose mode, |
---|
| 1562 | // and otherwise they are only printed if they are deemed to |
---|
| 1563 | // be outside of reasonable operating bounds. |
---|
| 1564 | //---------------------------------------------------------------------- |
---|
| 1565 | |
---|
| 1566 | #define log2(x) (log(x)/log(2.0)) // log base 2 |
---|
| 1567 | |
---|
| 1568 | void treeStats( |
---|
| 1569 | ostream &out, // output stream |
---|
| 1570 | ANNbool verbose) // print stats |
---|
| 1571 | { |
---|
| 1572 | const int MIN_PTS = 20; // min no. pts for checking |
---|
| 1573 | const float MAX_FRAC_TL = 0.50; // max frac of triv leaves |
---|
| 1574 | const float MAX_AVG_AR = 20; // max average aspect ratio |
---|
| 1575 | |
---|
| 1576 | ANNkdStats st; // statistics structure |
---|
| 1577 | |
---|
| 1578 | the_tree->getStats(st); // get statistics |
---|
| 1579 | // total number of nodes |
---|
| 1580 | int n_nodes = st.n_lf + st.n_spl + st.n_shr; |
---|
| 1581 | // should be O(n/bs) |
---|
| 1582 | int opt_n_nodes = (int) (2*(float(st.n_pts)/st.bkt_size)); |
---|
| 1583 | int too_many_nodes = 10*opt_n_nodes; |
---|
| 1584 | if (st.n_pts >= MIN_PTS && n_nodes > too_many_nodes) { |
---|
| 1585 | out << "-----------------------------------------------------------\n"; |
---|
| 1586 | out << "Warning: The tree has more than 10x as many nodes as points.\n"; |
---|
| 1587 | out << "You may want to consider a different split or shrink method.\n"; |
---|
| 1588 | out << "-----------------------------------------------------------\n"; |
---|
| 1589 | verbose = ANNtrue; |
---|
| 1590 | } |
---|
| 1591 | // fraction of trivial leaves |
---|
| 1592 | float frac_tl = (st.n_lf == 0 ? 0 : ((float) st.n_tl)/ st.n_lf); |
---|
| 1593 | if (st.n_pts >= MIN_PTS && frac_tl > MAX_FRAC_TL) { |
---|
| 1594 | out << "-----------------------------------------------------------\n"; |
---|
| 1595 | out << "Warning: A significant fraction of leaves contain no points.\n"; |
---|
| 1596 | out << "You may want to consider a different split or shrink method.\n"; |
---|
| 1597 | out << "-----------------------------------------------------------\n"; |
---|
| 1598 | verbose = ANNtrue; |
---|
| 1599 | } |
---|
| 1600 | // depth should be O(dim*log n) |
---|
| 1601 | int too_many_levels = (int) (2.0 * st.dim * log2((double) st.n_pts)); |
---|
| 1602 | int opt_levels = (int) log2(double(st.n_pts)/st.bkt_size); |
---|
| 1603 | if (st.n_pts >= MIN_PTS && st.depth > too_many_levels) { |
---|
| 1604 | out << "-----------------------------------------------------------\n"; |
---|
| 1605 | out << "Warning: The tree is more than 2x as deep as (dim*log n).\n"; |
---|
| 1606 | out << "You may want to consider a different split or shrink method.\n"; |
---|
| 1607 | out << "-----------------------------------------------------------\n"; |
---|
| 1608 | verbose = ANNtrue; |
---|
| 1609 | } |
---|
| 1610 | // average leaf aspect ratio |
---|
| 1611 | if (st.n_pts >= MIN_PTS && st.avg_ar > MAX_AVG_AR) { |
---|
| 1612 | out << "-----------------------------------------------------------\n"; |
---|
| 1613 | out << "Warning: Average aspect ratio of cells is quite large.\n"; |
---|
| 1614 | out << "This may slow queries depending on the point distribution.\n"; |
---|
| 1615 | out << "-----------------------------------------------------------\n"; |
---|
| 1616 | verbose = ANNtrue; |
---|
| 1617 | } |
---|
| 1618 | |
---|
| 1619 | //------------------------------------------------------------------ |
---|
| 1620 | // Print summaries if requested |
---|
| 1621 | //------------------------------------------------------------------ |
---|
| 1622 | if (verbose) { // output statistics |
---|
| 1623 | out << " (Structure Statistics:\n"; |
---|
| 1624 | out << " n_nodes = " << n_nodes |
---|
| 1625 | << " (opt = " << opt_n_nodes |
---|
| 1626 | << ", best if < " << too_many_nodes << ")\n" |
---|
| 1627 | << " n_leaves = " << st.n_lf |
---|
| 1628 | << " (" << st.n_tl << " contain no points)\n" |
---|
| 1629 | << " n_splits = " << st.n_spl << "\n" |
---|
| 1630 | << " n_shrinks = " << st.n_shr << "\n"; |
---|
| 1631 | out << " empty_leaves = " << frac_tl*100 |
---|
| 1632 | << " percent (best if < " << MAX_FRAC_TL*100 << " percent)\n"; |
---|
| 1633 | out << " depth = " << st.depth |
---|
| 1634 | << " (opt = " << opt_levels |
---|
| 1635 | << ", best if < " << too_many_levels << ")\n"; |
---|
| 1636 | out << " avg_aspect_ratio = " << st.avg_ar |
---|
| 1637 | << " (best if < " << MAX_AVG_AR << ")\n"; |
---|
| 1638 | out << " )\n"; |
---|
| 1639 | } |
---|
| 1640 | } |
---|