1 /***********************************************************************
2  * Software License Agreement (BSD License)
3  *
4  * Copyright 2008-2009  Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
5  * Copyright 2008-2009  David G. Lowe (lowe@cs.ubc.ca). All rights reserved.
6  *
7  * THE BSD LICENSE
8  *
9  * Redistribution and use in source and binary forms, with or without
10  * modification, are permitted provided that the following conditions
11  * are met:
12  *
13  * 1. Redistributions of source code must retain the above copyright
14  *    notice, this list of conditions and the following disclaimer.
15  * 2. Redistributions in binary form must reproduce the above copyright
16  *    notice, this list of conditions and the following disclaimer in the
17  *    documentation and/or other materials provided with the distribution.
18  *
19  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
20  * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
21  * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
22  * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
23  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
24  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
25  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
26  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
28  * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29  *************************************************************************/
30 #ifndef OPENCV_FLANN_AUTOTUNED_INDEX_H_
31 #define OPENCV_FLANN_AUTOTUNED_INDEX_H_
32 
33 #include <sstream>
34 
35 #include "general.h"
36 #include "nn_index.h"
37 #include "ground_truth.h"
38 #include "index_testing.h"
39 #include "sampling.h"
40 #include "kdtree_index.h"
41 #include "kdtree_single_index.h"
42 #include "kmeans_index.h"
43 #include "composite_index.h"
44 #include "linear_index.h"
45 #include "logger.h"
46 
47 namespace cvflann
48 {
49 
50 template<typename Distance>
51 NNIndex<Distance>* create_index_by_type(const Matrix<typename Distance::ElementType>& dataset, const IndexParams& params, const Distance& distance);
52 
53 
54 struct AutotunedIndexParams : public IndexParams
55 {
56     AutotunedIndexParams(float target_precision = 0.8, float build_weight = 0.01, float memory_weight = 0, float sample_fraction = 0.1)
57     {
58         (*this)["algorithm"] = FLANN_INDEX_AUTOTUNED;
59         // precision desired (used for autotuning, -1 otherwise)
60         (*this)["target_precision"] = target_precision;
61         // build tree time weighting factor
62         (*this)["build_weight"] = build_weight;
63         // index memory weighting factor
64         (*this)["memory_weight"] = memory_weight;
65         // what fraction of the dataset to use for autotuning
66         (*this)["sample_fraction"] = sample_fraction;
67     }
68 };
69 
70 
71 template <typename Distance>
72 class AutotunedIndex : public NNIndex<Distance>
73 {
74 public:
75     typedef typename Distance::ElementType ElementType;
76     typedef typename Distance::ResultType DistanceType;
77 
78     AutotunedIndex(const Matrix<ElementType>& inputData, const IndexParams& params = AutotunedIndexParams(), Distance d = Distance()) :
dataset_(inputData)79         dataset_(inputData), distance_(d)
80     {
81         target_precision_ = get_param(params, "target_precision",0.8f);
82         build_weight_ =  get_param(params,"build_weight", 0.01f);
83         memory_weight_ = get_param(params, "memory_weight", 0.0f);
84         sample_fraction_ = get_param(params,"sample_fraction", 0.1f);
85         bestIndex_ = NULL;
86         speedup_ = 0;
87     }
88 
89     AutotunedIndex(const AutotunedIndex&);
90     AutotunedIndex& operator=(const AutotunedIndex&);
91 
~AutotunedIndex()92     virtual ~AutotunedIndex()
93     {
94         if (bestIndex_ != NULL) {
95             delete bestIndex_;
96             bestIndex_ = NULL;
97         }
98     }
99 
100     /**
101      *          Method responsible with building the index.
102      */
buildIndex()103     virtual void buildIndex() CV_OVERRIDE
104     {
105         std::ostringstream stream;
106         bestParams_ = estimateBuildParams();
107         print_params(bestParams_, stream);
108         Logger::info("----------------------------------------------------\n");
109         Logger::info("Autotuned parameters:\n");
110         Logger::info("%s", stream.str().c_str());
111         Logger::info("----------------------------------------------------\n");
112 
113         bestIndex_ = create_index_by_type(dataset_, bestParams_, distance_);
114         bestIndex_->buildIndex();
115         speedup_ = estimateSearchParams(bestSearchParams_);
116         stream.str(std::string());
117         print_params(bestSearchParams_, stream);
118         Logger::info("----------------------------------------------------\n");
119         Logger::info("Search parameters:\n");
120         Logger::info("%s", stream.str().c_str());
121         Logger::info("----------------------------------------------------\n");
122     }
123 
124     /**
125      *  Saves the index to a stream
126      */
saveIndex(FILE * stream)127     virtual void saveIndex(FILE* stream) CV_OVERRIDE
128     {
129         save_value(stream, (int)bestIndex_->getType());
130         bestIndex_->saveIndex(stream);
131         save_value(stream, get_param<int>(bestSearchParams_, "checks"));
132     }
133 
134     /**
135      *  Loads the index from a stream
136      */
loadIndex(FILE * stream)137     virtual void loadIndex(FILE* stream) CV_OVERRIDE
138     {
139         int index_type;
140 
141         load_value(stream, index_type);
142         IndexParams params;
143         params["algorithm"] = (flann_algorithm_t)index_type;
144         bestIndex_ = create_index_by_type<Distance>(dataset_, params, distance_);
145         bestIndex_->loadIndex(stream);
146         int checks;
147         load_value(stream, checks);
148         bestSearchParams_["checks"] = checks;
149     }
150 
151     /**
152      *      Method that searches for nearest-neighbors
153      */
findNeighbors(ResultSet<DistanceType> & result,const ElementType * vec,const SearchParams & searchParams)154     virtual void findNeighbors(ResultSet<DistanceType>& result, const ElementType* vec, const SearchParams& searchParams) CV_OVERRIDE
155     {
156         int checks = get_param<int>(searchParams,"checks",FLANN_CHECKS_AUTOTUNED);
157         if (checks == FLANN_CHECKS_AUTOTUNED) {
158             bestIndex_->findNeighbors(result, vec, bestSearchParams_);
159         }
160         else {
161             bestIndex_->findNeighbors(result, vec, searchParams);
162         }
163     }
164 
165 
getParameters()166     IndexParams getParameters() const CV_OVERRIDE
167     {
168         return bestIndex_->getParameters();
169     }
170 
getSearchParameters()171     SearchParams getSearchParameters() const
172     {
173         return bestSearchParams_;
174     }
175 
getSpeedup()176     float getSpeedup() const
177     {
178         return speedup_;
179     }
180 
181 
182     /**
183      *      Number of features in this index.
184      */
size()185     virtual size_t size() const CV_OVERRIDE
186     {
187         return bestIndex_->size();
188     }
189 
190     /**
191      *  The length of each vector in this index.
192      */
veclen()193     virtual size_t veclen() const CV_OVERRIDE
194     {
195         return bestIndex_->veclen();
196     }
197 
198     /**
199      * The amount of memory (in bytes) this index uses.
200      */
usedMemory()201     virtual int usedMemory() const CV_OVERRIDE
202     {
203         return bestIndex_->usedMemory();
204     }
205 
206     /**
207      * Algorithm name
208      */
getType()209     virtual flann_algorithm_t getType() const CV_OVERRIDE
210     {
211         return FLANN_INDEX_AUTOTUNED;
212     }
213 
214 private:
215 
216     struct CostData
217     {
218         float searchTimeCost;
219         float buildTimeCost;
220         float memoryCost;
221         float totalCost;
222         IndexParams params;
223     };
224 
evaluate_kmeans(CostData & cost)225     void evaluate_kmeans(CostData& cost)
226     {
227         StartStopTimer t;
228         int checks;
229         const int nn = 1;
230 
231         Logger::info("KMeansTree using params: max_iterations=%d, branching=%d\n",
232                      get_param<int>(cost.params,"iterations"),
233                      get_param<int>(cost.params,"branching"));
234         KMeansIndex<Distance> kmeans(sampledDataset_, cost.params, distance_);
235         // measure index build time
236         t.start();
237         kmeans.buildIndex();
238         t.stop();
239         float buildTime = (float)t.value;
240 
241         // measure search time
242         float searchTime = test_index_precision(kmeans, sampledDataset_, testDataset_, gt_matches_, target_precision_, checks, distance_, nn);
243 
244         float datasetMemory = float(sampledDataset_.rows * sampledDataset_.cols * sizeof(float));
245         cost.memoryCost = (kmeans.usedMemory() + datasetMemory) / datasetMemory;
246         cost.searchTimeCost = searchTime;
247         cost.buildTimeCost = buildTime;
248         Logger::info("KMeansTree buildTime=%g, searchTime=%g, build_weight=%g\n", buildTime, searchTime, build_weight_);
249     }
250 
251 
evaluate_kdtree(CostData & cost)252     void evaluate_kdtree(CostData& cost)
253     {
254         StartStopTimer t;
255         int checks;
256         const int nn = 1;
257 
258         Logger::info("KDTree using params: trees=%d\n", get_param<int>(cost.params,"trees"));
259         KDTreeIndex<Distance> kdtree(sampledDataset_, cost.params, distance_);
260 
261         t.start();
262         kdtree.buildIndex();
263         t.stop();
264         float buildTime = (float)t.value;
265 
266         //measure search time
267         float searchTime = test_index_precision(kdtree, sampledDataset_, testDataset_, gt_matches_, target_precision_, checks, distance_, nn);
268 
269         float datasetMemory = float(sampledDataset_.rows * sampledDataset_.cols * sizeof(float));
270         cost.memoryCost = (kdtree.usedMemory() + datasetMemory) / datasetMemory;
271         cost.searchTimeCost = searchTime;
272         cost.buildTimeCost = buildTime;
273         Logger::info("KDTree buildTime=%g, searchTime=%g\n", buildTime, searchTime);
274     }
275 
276 
277     //    struct KMeansSimpleDownhillFunctor {
278     //
279     //        Autotune& autotuner;
280     //        KMeansSimpleDownhillFunctor(Autotune& autotuner_) : autotuner(autotuner_) {}
281     //
282     //        float operator()(int* params) {
283     //
284     //            float maxFloat = numeric_limits<float>::max();
285     //
286     //            if (params[0]<2) return maxFloat;
287     //            if (params[1]<0) return maxFloat;
288     //
289     //            CostData c;
290     //            c.params["algorithm"] = KMEANS;
291     //            c.params["centers-init"] = CENTERS_RANDOM;
292     //            c.params["branching"] = params[0];
293     //            c.params["max-iterations"] = params[1];
294     //
295     //            autotuner.evaluate_kmeans(c);
296     //
297     //            return c.timeCost;
298     //
299     //        }
300     //    };
301     //
302     //    struct KDTreeSimpleDownhillFunctor {
303     //
304     //        Autotune& autotuner;
305     //        KDTreeSimpleDownhillFunctor(Autotune& autotuner_) : autotuner(autotuner_) {}
306     //
307     //        float operator()(int* params) {
308     //            float maxFloat = numeric_limits<float>::max();
309     //
310     //            if (params[0]<1) return maxFloat;
311     //
312     //            CostData c;
313     //            c.params["algorithm"] = KDTREE;
314     //            c.params["trees"] = params[0];
315     //
316     //            autotuner.evaluate_kdtree(c);
317     //
318     //            return c.timeCost;
319     //
320     //        }
321     //    };
322 
323 
324 
optimizeKMeans(std::vector<CostData> & costs)325     void optimizeKMeans(std::vector<CostData>& costs)
326     {
327         Logger::info("KMEANS, Step 1: Exploring parameter space\n");
328 
329         // explore kmeans parameters space using combinations of the parameters below
330         int maxIterations[] = { 1, 5, 10, 15 };
331         int branchingFactors[] = { 16, 32, 64, 128, 256 };
332 
333         int kmeansParamSpaceSize = FLANN_ARRAY_LEN(maxIterations) * FLANN_ARRAY_LEN(branchingFactors);
334         costs.reserve(costs.size() + kmeansParamSpaceSize);
335 
336         // evaluate kmeans for all parameter combinations
337         for (size_t i = 0; i < FLANN_ARRAY_LEN(maxIterations); ++i) {
338             for (size_t j = 0; j < FLANN_ARRAY_LEN(branchingFactors); ++j) {
339                 CostData cost;
340                 cost.params["algorithm"] = FLANN_INDEX_KMEANS;
341                 cost.params["centers_init"] = FLANN_CENTERS_RANDOM;
342                 cost.params["iterations"] = maxIterations[i];
343                 cost.params["branching"] = branchingFactors[j];
344 
345                 evaluate_kmeans(cost);
346                 costs.push_back(cost);
347             }
348         }
349 
350         //         Logger::info("KMEANS, Step 2: simplex-downhill optimization\n");
351         //
352         //         const int n = 2;
353         //         // choose initial simplex points as the best parameters so far
354         //         int kmeansNMPoints[n*(n+1)];
355         //         float kmeansVals[n+1];
356         //         for (int i=0;i<n+1;++i) {
357         //             kmeansNMPoints[i*n] = (int)kmeansCosts[i].params["branching"];
358         //             kmeansNMPoints[i*n+1] = (int)kmeansCosts[i].params["max-iterations"];
359         //             kmeansVals[i] = kmeansCosts[i].timeCost;
360         //         }
361         //         KMeansSimpleDownhillFunctor kmeans_cost_func(*this);
362         //         // run optimization
363         //         optimizeSimplexDownhill(kmeansNMPoints,n,kmeans_cost_func,kmeansVals);
364         //         // store results
365         //         for (int i=0;i<n+1;++i) {
366         //             kmeansCosts[i].params["branching"] = kmeansNMPoints[i*2];
367         //             kmeansCosts[i].params["max-iterations"] = kmeansNMPoints[i*2+1];
368         //             kmeansCosts[i].timeCost = kmeansVals[i];
369         //         }
370     }
371 
372 
optimizeKDTree(std::vector<CostData> & costs)373     void optimizeKDTree(std::vector<CostData>& costs)
374     {
375         Logger::info("KD-TREE, Step 1: Exploring parameter space\n");
376 
377         // explore kd-tree parameters space using the parameters below
378         int testTrees[] = { 1, 4, 8, 16, 32 };
379 
380         // evaluate kdtree for all parameter combinations
381         for (size_t i = 0; i < FLANN_ARRAY_LEN(testTrees); ++i) {
382             CostData cost;
383             cost.params["algorithm"] = FLANN_INDEX_KDTREE;
384             cost.params["trees"] = testTrees[i];
385 
386             evaluate_kdtree(cost);
387             costs.push_back(cost);
388         }
389 
390         //         Logger::info("KD-TREE, Step 2: simplex-downhill optimization\n");
391         //
392         //         const int n = 1;
393         //         // choose initial simplex points as the best parameters so far
394         //         int kdtreeNMPoints[n*(n+1)];
395         //         float kdtreeVals[n+1];
396         //         for (int i=0;i<n+1;++i) {
397         //             kdtreeNMPoints[i] = (int)kdtreeCosts[i].params["trees"];
398         //             kdtreeVals[i] = kdtreeCosts[i].timeCost;
399         //         }
400         //         KDTreeSimpleDownhillFunctor kdtree_cost_func(*this);
401         //         // run optimization
402         //         optimizeSimplexDownhill(kdtreeNMPoints,n,kdtree_cost_func,kdtreeVals);
403         //         // store results
404         //         for (int i=0;i<n+1;++i) {
405         //             kdtreeCosts[i].params["trees"] = kdtreeNMPoints[i];
406         //             kdtreeCosts[i].timeCost = kdtreeVals[i];
407         //         }
408     }
409 
410     /**
411      *  Chooses the best nearest-neighbor algorithm and estimates the optimal
412      *  parameters to use when building the index (for a given precision).
413      *  Returns a dictionary with the optimal parameters.
414      */
estimateBuildParams()415     IndexParams estimateBuildParams()
416     {
417         std::vector<CostData> costs;
418 
419         int sampleSize = int(sample_fraction_ * dataset_.rows);
420         int testSampleSize = std::min(sampleSize / 10, 1000);
421 
422         Logger::info("Entering autotuning, dataset size: %d, sampleSize: %d, testSampleSize: %d, target precision: %g\n", dataset_.rows, sampleSize, testSampleSize, target_precision_);
423 
424         // For a very small dataset, it makes no sense to build any fancy index, just
425         // use linear search
426         if (testSampleSize < 10) {
427             Logger::info("Choosing linear, dataset too small\n");
428             return LinearIndexParams();
429         }
430 
431         // We use a fraction of the original dataset to speedup the autotune algorithm
432         sampledDataset_ = random_sample(dataset_, sampleSize);
433         // We use a cross-validation approach, first we sample a testset from the dataset
434         testDataset_ = random_sample(sampledDataset_, testSampleSize, true);
435 
436         // We compute the ground truth using linear search
437         Logger::info("Computing ground truth... \n");
438         gt_matches_ = Matrix<int>(new int[testDataset_.rows], testDataset_.rows, 1);
439         StartStopTimer t;
440         t.start();
441         compute_ground_truth<Distance>(sampledDataset_, testDataset_, gt_matches_, 0, distance_);
442         t.stop();
443 
444         CostData linear_cost;
445         linear_cost.searchTimeCost = (float)t.value;
446         linear_cost.buildTimeCost = 0;
447         linear_cost.memoryCost = 0;
448         linear_cost.params["algorithm"] = FLANN_INDEX_LINEAR;
449 
450         costs.push_back(linear_cost);
451 
452         // Start parameter autotune process
453         Logger::info("Autotuning parameters...\n");
454 
455         optimizeKMeans(costs);
456         optimizeKDTree(costs);
457 
458         float bestTimeCost = costs[0].searchTimeCost;
459         for (size_t i = 0; i < costs.size(); ++i) {
460             float timeCost = costs[i].buildTimeCost * build_weight_ + costs[i].searchTimeCost;
461             if (timeCost < bestTimeCost) {
462                 bestTimeCost = timeCost;
463             }
464         }
465 
466         float bestCost = costs[0].searchTimeCost / bestTimeCost;
467         IndexParams bestParams = costs[0].params;
468         if (bestTimeCost > 0) {
469             for (size_t i = 0; i < costs.size(); ++i) {
470                 float crtCost = (costs[i].buildTimeCost * build_weight_ + costs[i].searchTimeCost) / bestTimeCost +
471                                 memory_weight_ * costs[i].memoryCost;
472                 if (crtCost < bestCost) {
473                     bestCost = crtCost;
474                     bestParams = costs[i].params;
475                 }
476             }
477         }
478 
479         delete[] gt_matches_.data;
480         delete[] testDataset_.data;
481         delete[] sampledDataset_.data;
482 
483         return bestParams;
484     }
485 
486 
487 
488     /**
489      *  Estimates the search time parameters needed to get the desired precision.
490      *  Precondition: the index is built
491      *  Postcondition: the searchParams will have the optimum params set, also the speedup obtained over linear search.
492      */
estimateSearchParams(SearchParams & searchParams)493     float estimateSearchParams(SearchParams& searchParams)
494     {
495         const int nn = 1;
496         const size_t SAMPLE_COUNT = 1000;
497 
498         assert(bestIndex_ != NULL); // must have a valid index
499 
500         float speedup = 0;
501 
502         int samples = (int)std::min(dataset_.rows / 10, SAMPLE_COUNT);
503         if (samples > 0) {
504             Matrix<ElementType> testDataset = random_sample(dataset_, samples);
505 
506             Logger::info("Computing ground truth\n");
507 
508             // we need to compute the ground truth first
509             Matrix<int> gt_matches(new int[testDataset.rows], testDataset.rows, 1);
510             StartStopTimer t;
511             t.start();
512             compute_ground_truth<Distance>(dataset_, testDataset, gt_matches, 1, distance_);
513             t.stop();
514             float linear = (float)t.value;
515 
516             int checks;
517             Logger::info("Estimating number of checks\n");
518 
519             float searchTime;
520             float cb_index;
521             if (bestIndex_->getType() == FLANN_INDEX_KMEANS) {
522                 Logger::info("KMeans algorithm, estimating cluster border factor\n");
523                 KMeansIndex<Distance>* kmeans = (KMeansIndex<Distance>*)bestIndex_;
524                 float bestSearchTime = -1;
525                 float best_cb_index = -1;
526                 int best_checks = -1;
527                 for (cb_index = 0; cb_index < 1.1f; cb_index += 0.2f) {
528                     kmeans->set_cb_index(cb_index);
529                     searchTime = test_index_precision(*kmeans, dataset_, testDataset, gt_matches, target_precision_, checks, distance_, nn, 1);
530                     if ((searchTime < bestSearchTime) || (bestSearchTime == -1)) {
531                         bestSearchTime = searchTime;
532                         best_cb_index = cb_index;
533                         best_checks = checks;
534                     }
535                 }
536                 searchTime = bestSearchTime;
537                 cb_index = best_cb_index;
538                 checks = best_checks;
539 
540                 kmeans->set_cb_index(best_cb_index);
541                 Logger::info("Optimum cb_index: %g\n", cb_index);
542                 bestParams_["cb_index"] = cb_index;
543             }
544             else {
545                 searchTime = test_index_precision(*bestIndex_, dataset_, testDataset, gt_matches, target_precision_, checks, distance_, nn, 1);
546             }
547 
548             Logger::info("Required number of checks: %d \n", checks);
549             searchParams["checks"] = checks;
550 
551             speedup = linear / searchTime;
552 
553             delete[] gt_matches.data;
554             delete[] testDataset.data;
555         }
556 
557         return speedup;
558     }
559 
560 private:
561     NNIndex<Distance>* bestIndex_;
562 
563     IndexParams bestParams_;
564     SearchParams bestSearchParams_;
565 
566     Matrix<ElementType> sampledDataset_;
567     Matrix<ElementType> testDataset_;
568     Matrix<int> gt_matches_;
569 
570     float speedup_;
571 
572     /**
573      * The dataset used by this index
574      */
575     const Matrix<ElementType> dataset_;
576 
577     /**
578      * Index parameters
579      */
580     float target_precision_;
581     float build_weight_;
582     float memory_weight_;
583     float sample_fraction_;
584 
585     Distance distance_;
586 
587 
588 };
589 }
590 
591 #endif /* OPENCV_FLANN_AUTOTUNED_INDEX_H_ */
592