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 | } |
---|