From 3f93e916246d8c03d916596ed5653f5ad2037316 Mon Sep 17 00:00:00 2001 From: jean-pierre both Date: Sat, 18 Apr 2020 11:23:32 +0200 Subject: [PATCH] Reinitialize .git due to large packed file of 76Mb --- .gitignore | 5 + Cargo.toml | 81 ++ LICENSE-APACHE | 13 + LICENSE-MIT | 25 + README.md | 95 ++ examples/ann-glove25-angular.rs | 147 +++ examples/ann-mnist-784-euclidean.rs | 107 +++ examples/random.rs | 66 ++ src/annhdf5.rs | 217 +++++ src/api.rs | 126 +++ src/dist.rs | 869 ++++++++++++++++++ src/hnsw.rs | 1300 +++++++++++++++++++++++++++ src/hnswio.rs | 689 ++++++++++++++ src/lib.rs | 43 + src/libext.rs | 782 ++++++++++++++++ src/prelude.rs | 5 + src/test.rs | 316 +++++++ 17 files changed, 4886 insertions(+) create mode 100644 .gitignore create mode 100644 Cargo.toml create mode 100644 LICENSE-APACHE create mode 100644 LICENSE-MIT create mode 100644 README.md create mode 100644 examples/ann-glove25-angular.rs create mode 100644 examples/ann-mnist-784-euclidean.rs create mode 100644 examples/random.rs create mode 100644 src/annhdf5.rs create mode 100644 src/api.rs create mode 100644 src/dist.rs create mode 100644 src/hnsw.rs create mode 100644 src/hnswio.rs create mode 100644 src/lib.rs create mode 100644 src/libext.rs create mode 100644 src/prelude.rs create mode 100644 src/test.rs diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..34bfb7b --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +target/** +Runs +Cargo.lock +rls* +dumpreloadtest* diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..50d6cf3 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,81 @@ +[package] +name = "hnsw_rs" +version = "0.1.3" +authors = ["jeanpierre.both@gmail.com"] +description = "Ann based on Hierarchical Navigable Small World Graphs from Yu.A. Malkov and D.A Yashunin" +license = "MIT/Apache-2.0" +readme = "README.md" +keywords = ["algorithms", "ann", "hnsw"] +repository = "https://github.com/jean-pierreBoth/hnswlib-rs" +edition= "2018" + +[features] + +# declare a feature with no dependancy to get some modulated debug print +# to be run with cargo build --features verbose_1 +#verbose_1 = [ ] + +[profile.release] +lto = true +opt-level = 3 + +[lib] +# cargo rustc --lib -- --crate-type dylib [or staticlib] or rlib (default) +# if we want to avoid specifying in advance crate-type +lib = "hnswlib" +path = "src/lib.rs" +#crate-type = ["dylib"] + + +[[examples]] +name = "random" +path = "examples/random.rs" + +[[examples]] +name = "ann-glove" +path = "examples/ann-glove25-angular.rs" + + + +[[examples]] +name = "ann-mnist" +path = "examples/ann-mnist-784-euclidean.rs" + +#[[example]] + +[dependencies] +# default is version spec is ^ meaning can update up to max non null version number +# cargo doc --no-deps avoid dependencies doc generation +# + +serde= {version = "1.0", doc = false} +serde_derive={ version = "1.0", doc = false} + + +# for // +crossbeam-utils = "0.7" +crossbeam-channel = "0.4" +parking_lot = "0.9" +rayon = {version = "0.9.0", doc = false} +num_cpus = {version = "1.8.0", doc = false} +simdeez = {version = "0.6", doc = false} + +cpu-time = {version = "0.1", doc = false} +time = {version = "0.1.39", doc = false} + +ndarray = {version = "0.12", doc = false} + +clap = {version = "2.29", doc = false} +# for hashing . hashbrown still needed beccause of get_key_value(&key) +hashbrown = {version = "0.3", doc = false} +skiplist= {version = "0.2.10", doc=false} + +rand = {version = "0.7", doc = false} +lazy_static = { version = "1.3", doc = false} +typename = {version = "0.1", doc = false} +# for benchmark reading +hdf5 = {version="0.5", doc = false} +# decreasing order of log for debug build : (max_level_)trace debug info warn error off +# decreasing order of log for release build (release_max_level_) .. idem +log = { version = "0.4", features = ["max_level_debug", "release_max_level_info"] } +simple_logger = { version = "1.0"} diff --git a/LICENSE-APACHE b/LICENSE-APACHE new file mode 100644 index 0000000..d8afa4c --- /dev/null +++ b/LICENSE-APACHE @@ -0,0 +1,13 @@ +Copyright 2020 jean-pierre.both + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. diff --git a/LICENSE-MIT b/LICENSE-MIT new file mode 100644 index 0000000..531ac46 --- /dev/null +++ b/LICENSE-MIT @@ -0,0 +1,25 @@ +Copyright (c) 2020 jean-pierre.both + +Permission is hereby granted, free of charge, to any +person obtaining a copy of this software and associated +documentation files (the "Software"), to deal in the +Software without restriction, including without +limitation the rights to use, copy, modify, merge, +publish, distribute, sublicense, and/or sell copies of +the Software, and to permit persons to whom the Software +is furnished to do so, subject to the following +conditions: + +The above copyright notice and this permission notice +shall be included in all copies or substantial portions +of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF +ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED +TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A +PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT +SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION +OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR +IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +DEALINGS IN THE SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..3419151 --- /dev/null +++ b/README.md @@ -0,0 +1,95 @@ +# hnsw-rs + +This crate provides a Rust implementation of the paper by Yu.A. Malkov and D.A Yashunin: + +"Efficient and Robust approximate nearest neighbours using Hierarchical Navigable Small World Graphs" (2016,2018) +[https://arxiv.org/abs/1603.09320] + +## License + +Licensed under either of + +* Apache License, Version 2.0, [LICENSE-APACHE](LICENSE-APACHE) or +* MIT license [LICENSE-MIT](LICENSE-MIT) or + +at your option. + +This software was written on my own while working at [CEA](http://www.cea.fr/), [CEA-LIST](http://www-list.cea.fr/en/) + +## Functionalities + +The crate provides: + +* usual distances as L1, L2, Cosine, Jaccard, Hamming for vectors of standard numeric types. + +* Hellinger and Jeffreys distances between probability distributions (f32 and f64). It must be noted that the Jeffreys distance +(a symetrized Kullback-Leibler divergence) do not satisfy the triangle inequality. (Neither Cosine distance !). + +* A structure to enable the user to implement its own distances. It takes as data, vectors of types T:Copy+Clone+Send+Sync. + +* An interface towards C and more specifically to the [Julia](https://julialang.org/) language. +See the companion Julia package [HnswAnn.jl](https://gitlab.com/jpboth/HnswAnn.jl) and the building paragraph for some help for Julia users. + +* Dump and reload functions (Cf module hnswio) to store the graph once it is built. As the time necessary to compute the graph can be important it can be useful to store it for future use. + +## Implementation + +The graph construction and searches are multithreaded with the **parking_lot** crate (See **parallel_insert_data** and **parallel_search_neighbours** functions and also examples files). +For the heavily used case (f32) we provide simd avx2 implementation in distance computations +currently based on the **simdeez** crate. + +## Building + +By default the crate is a standalone project and builds a static libray and executable. +To be used with the companion Julia package it is necessary to build a dynamic library. +This can be done by just uncommenting (i.e get rid of the #) in file Cargo.toml the line: + +*#crate-type = ["dylib"]* + +and rerun the command: cargo build --release. + +This will generate a .so file in the target/release directory. + +## Algorithm and Input Parameters + +The algorithm stores points in layers (at most 16), and a graph is constructed to enable a search from less densely populated levels to most densely populated levels by constructing links from less dense layers to the most dense layer (level 0). + +Roughly the algorithm goes along runs as follows: + +Upon insertion, the level ***l*** of a new point is sampled with an exponential law, limiting the number of levels to 16, +so that level 0 is the most densely populated layer, upper layers being exponentially less populated as level increases. +The nearest neighbour of the point is searched in lookup tables from the upper level to the level just above its layer (***l***), so we should arrive near the new point at its level at a relatively low cost. Then the ***max_nb_connection*** nearest neighbours are searched in neighbours of neighbours table (with a reverse updating of tables) recursively from its layer ***l*** down to the most populated level 0. + +The scale parameter of the exponential law depends on the maximum number of connection possible for a point (parameter ***max_nb_connection***) to others. +Explicitly the scale parameter is chosen as : `scale=1/ln(max_nb_connection)`. + +The main parameters occuring in constructing the graph or in searching are: + +* max_nb_connection (in hnsw initialization) + The maximum number of links from one point to others. Values ranging from 16 to 64 are standard initialising values, the higher the more time consuming. + +* ef_construction (in hnsw initialization) + This parameter controls the width of the search for neighbours during insertion. Values from 200 to 800 are standard initialising values, the higher the more time consuming. + +* max_layer (in hnsw initialization) + The maximum number of layers in graph. Must be less or equal than 16. + +* ef_arg (in search methods) + This parameter controls the width of the search in the lowest level, it must be greater than number of neighbours asked but can be less than ***ef_construction***. + As a rule of thumb could be between the number of neighbours we will ask for (knbn arg in search method) and max_nb_connection. + +* keep_pruned and extend_candidates. + These parameters are described in the paper by Malkov and Yashunin can be used to + modify the search strategy. The interested user should check the paper to see the impact. By default + the values are as recommended in the paper. + +## Examples and Benchmarks + +Some examples are taken from the [ann-benchmarks site](https://github.com/erikbern/ann-benchmarks) +and recall rates and request/s are given in comments in the examples files for some input parameters. +The annhdf5 module implements reading the standardized data files +of the [ann-benchmarks site](https://github.com/erikbern/ann-benchmarks), +just download the necessary benchmark data files and modify path in sources accordingly. +Then run: cargo build --examples --release. +It is possible in these examples to change from parallel searches to serial searches to check for speeds +or modify parameters to see the impact on performance. diff --git a/examples/ann-glove25-angular.rs b/examples/ann-glove25-angular.rs new file mode 100644 index 0000000..2934d13 --- /dev/null +++ b/examples/ann-glove25-angular.rs @@ -0,0 +1,147 @@ +use std::time::{Duration, SystemTime}; +use cpu_time::ProcessTime; + + +use typename::TypeName; + + +// glove 25 // 2.7 Ghz 4 cores 8Mb L3 k = 10 +// +// max_nb_conn ef_cons ef_search scale_factor extend keep pruned recall req/s last ratio +// 24 800 64 1. 1 0 0.928 4090 1.003 +// 24 800 64 1. 1 1 0.927 4594 1.003 +// 24 400, 48 1. 1 0 0.919 6349 1.0044 +// 24 800 48 1 1 1 0.918 5785 1.005 +// 24 400 32 1. 0 0 0.898 8662 +// 24 400 64 1. 1 0 0.930 4711 1.0027 +// 24 400 64 1. 1 1 0.921 4550 1.0039 +// 24 400 64 0.5 0 0 0.916 4896 1.0046 +// 24 1600 48 1 1 0 0.924 5380 1.0034 + +// 32 400 48 1 1 0 0.93 4706 1.0026 +// 32 800 64 1 1 0 0.94 3780. 1.0015 +// 32 1600 48 1 1 0 0.934 4455 1.0023 +// 48 1600 48 1 1 0 0.945 3253 1.00098 + +// 24 400 48 1 1 0 0.92 6036. 1.0038 +// 48 800 48 1 1 0 0.935 4018 1.002 +// 48 800 64 1 1 0 0.942 3091 1.0014 +// 48 800 64 0.5 1 0 0.9395 3234 1.00167 +// 48 800 64 1 1 1 0.9435 2640 1.00126 + + +// k = 100 + +// 24 800 48 1 1 0 0.96 2432 1.004 +// 48 800 128 1 1 0 0.979 1626 1.001 + +extern crate hnsw_rs; + +use hnsw_rs::prelude::*; + + +pub fn main() { + + let _res = simple_logger::init(); + let parallel = true; + // + let fname = String::from("/home.2/Data/ANN/glove-25-angular.hdf5"); + println!("\n\n test_load_hdf5 {:?}", fname); + // now recall that data are stored in row order. + let mut anndata = AnnBenchmarkData::new(fname).unwrap(); + // pre normalisation to use Dot computations instead of Cosine + anndata.do_l2_normalization(); + // run bench + let nb_elem = anndata.train_data.len(); + let max_nb_connection = 48; + let ef_c = 800; + println!(" max_nb_conn : {:?}, ef_construction : {:?} ", max_nb_connection, ef_c); + let nb_layer = 16.min((nb_elem as f32).ln().trunc() as usize); + println!(" number of elements to insert {:?} , setting max nb layer to {:?} ef_construction {:?}", nb_elem, nb_layer, ef_c); + let nb_search = anndata.test_data.len(); + println!(" number of search {:?}", nb_search); + // Hnsw allocation + let mut hnsw = Hnsw::::new(max_nb_connection, nb_elem, nb_layer, ef_c, DistDot{}); + hnsw.set_extend_candidates(true); + // parallel insertion + let start = ProcessTime::now(); + let now = SystemTime::now(); + let data_for_par_insertion = anndata.train_data.iter().map( |x| (&x.0, x.1)).collect(); + if parallel { + println!(" \n parallel insertion"); + hnsw.parallel_insert(&data_for_par_insertion); + } + else { + println!(" \n serial insertion"); + for d in data_for_par_insertion { + hnsw.insert(d); + } + } + let cpu_time: Duration = start.elapsed(); + // + println!("\n hnsw data insertion cpu time {:?} system time {:?} ", cpu_time, now.elapsed()); + hnsw.dump_layer_info(); + println!(" hnsw data nb point inserted {:?}", hnsw.get_nb_point()); + // + // Now the bench with 10 neighbours + // + let knbn = 10; + let ef_search = 48; + search(&mut hnsw, knbn, ef_search, &anndata); + + let knbn = 100; + let ef_search = 128; + search(&mut hnsw, knbn, ef_search, &anndata); +} + + +pub fn search(hnsw: &mut Hnsw, knbn : usize, ef_search: usize, anndata : & AnnBenchmarkData) + where Dist : Distance + Send + Sync + TypeName { + + println!("\n\n ef_search : {:?} knbn : {:?} ", ef_search, knbn); + let parallel = true; + // + let nb_elem = anndata.train_data.len(); + let nb_search = anndata.test_data.len(); + // + let mut recalls = Vec::::with_capacity(nb_elem); + let mut nb_returned = Vec::::with_capacity(nb_elem); + let mut last_distances_ratio = Vec::::with_capacity(nb_elem); + let mut knn_neighbours_for_tests = Vec::>::with_capacity(nb_elem); + hnsw.set_searching_mode(true); + println!("searching with ef : {:?}", ef_search); + let start = ProcessTime::now(); + let now = SystemTime::now(); + // search + if parallel { + println!(" \n parallel search"); + knn_neighbours_for_tests = hnsw.parallel_search(&anndata.test_data, knbn, ef_search); + } else { + println!(" \n serial search"); + for i in 0..anndata.test_data.len() { + let knn_neighbours : Vec = hnsw.search(&anndata.test_data[i], knbn, ef_search); + knn_neighbours_for_tests.push(knn_neighbours); + } + } + let cpu_time = start.elapsed(); + let search_cpu_time = cpu_time.as_micros() as f32; + let search_sys_time = now.elapsed().unwrap().as_micros() as f32; + println!("total cpu time for search requests {:?} , system time {:?} ", search_cpu_time, now.elapsed()); + // now compute recall rate + for i in 0..anndata.test_data.len() { + let max_dist = anndata.test_distances.row(i)[knbn-1]; + let knn_neighbours_d : Vec = knn_neighbours_for_tests[i].iter().map(|p| p.distance).collect(); + nb_returned.push(knn_neighbours_d.len()); + let recall = knn_neighbours_d.iter().filter(|d| *d <= &max_dist).count(); + recalls.push(recall); + let mut ratio = 0.; + if knn_neighbours_d.len() >= 1 { + ratio = knn_neighbours_d[knn_neighbours_d.len()-1]/max_dist; + } + last_distances_ratio.push(ratio); + } + let mean_recall = (recalls.iter().sum::() as f32)/((knbn * recalls.len()) as f32); + println!("\n mean fraction nb returned by search {:?} ", (nb_returned.iter().sum::() as f32)/ ((nb_returned.len() * knbn) as f32)); + println!("\n last distances ratio {:?} ", last_distances_ratio.iter().sum::() / last_distances_ratio.len() as f32); + println!("\n recall rate for {:?} is {:?} , nb req /s {:?}", anndata.fname, mean_recall, (nb_search as f32) * 1.0e+6_f32/search_sys_time); +} \ No newline at end of file diff --git a/examples/ann-mnist-784-euclidean.rs b/examples/ann-mnist-784-euclidean.rs new file mode 100644 index 0000000..914f2f4 --- /dev/null +++ b/examples/ann-mnist-784-euclidean.rs @@ -0,0 +1,107 @@ +use std::time::{Duration, SystemTime}; +use cpu_time::ProcessTime; +// search in serial mode +// max_nb_conn ef_cons ef_search scale_factor extend keep pruned recall req/s last ratio +// +// 12 400 12 1 0 0 0.917 6486 1.005 +// 24 400 24 1 1 0 0.9779 3456 1.001 + +// parallel mode +// max_nb_conn ef_cons ef_search scale_factor extend keep pruned recall req/s last ratio +// 12 400 12 1 0 0 0.977 12566 1.005 + +extern crate hnsw_rs; + +use hnsw_rs::prelude::*; + + +pub fn main() { + + let _res = simple_logger::init(); + let mut parallel = true; + // + let fname = String::from("/home.2/Data/ANN/fashion-mnist-784-euclidean.hdf5"); + println!("\n\n test_load_hdf5 {:?}", fname); + // now recall that data are stored in row order. + let anndata = AnnBenchmarkData::new(fname).unwrap(); + // run bench + let nb_elem = anndata.train_data.len(); + let max_nb_connection = 24; + let nb_layer = 16.min((nb_elem as f32).ln().trunc() as usize); + let ef_c = 400; + println!(" number of elements to insert {:?} , setting max nb layer to {:?} ef_construction {:?}", nb_elem, nb_layer, ef_c); + println!(" ====================================================================================="); + let nb_search = anndata.test_data.len(); + println!(" number of search {:?}", nb_search); + + let mut hnsw = Hnsw::::new(max_nb_connection, nb_elem, nb_layer, ef_c, DistL2{}); + hnsw.set_extend_candidates(false); + // parallel insertion + let mut start = ProcessTime::now(); + let mut now = SystemTime::now(); + let data_for_par_insertion = anndata.train_data.iter().map( |x| (&x.0, x.1)).collect(); + if parallel { + println!(" \n parallel insertion"); + hnsw.parallel_insert(&data_for_par_insertion); + } + else { + println!(" \n serial insertion"); + for d in data_for_par_insertion { + hnsw.insert(d); + } + } + let mut cpu_time: Duration = start.elapsed(); + // + println!("\n hnsw data insertion cpu time {:?} system time {:?} ", cpu_time, now.elapsed()); + hnsw.dump_layer_info(); + println!(" hnsw data nb point inserted {:?}", hnsw.get_nb_point()); + // + // Now the bench with 10 neighbours + // + let mut recalls = Vec::::with_capacity(nb_elem); + let mut nb_returned = Vec::::with_capacity(nb_elem); + let mut last_distances_ratio = Vec::::with_capacity(nb_elem); + let mut knn_neighbours_for_tests = Vec::>::with_capacity(nb_elem); + hnsw.set_searching_mode(false); + let knbn = 10; + let ef_c = max_nb_connection; + println!("\n searching with ef : {:?}", ef_c); + start = ProcessTime::now(); + now = SystemTime::now(); + // search + parallel = true; + if parallel { + println!(" \n parallel search"); + knn_neighbours_for_tests = hnsw.parallel_search(&anndata.test_data, knbn, ef_c); + } else { + println!(" \n serial search"); + for i in 0..anndata.test_data.len() { + let knn_neighbours : Vec = hnsw.search(&anndata.test_data[i], knbn, ef_c); + knn_neighbours_for_tests.push(knn_neighbours); + } + } + cpu_time = start.elapsed(); + let search_sys_time = now.elapsed().unwrap().as_micros() as f32; + let search_cpu_time = cpu_time.as_micros() as f32; + println!("total cpu time for search requests {:?} , system time {:?} ", search_cpu_time, search_sys_time); + // now compute recall rate + for i in 0..anndata.test_data.len() { + let true_distances = anndata.test_distances.row(i); + let max_dist = true_distances[knbn-1]; + let mut _knn_neighbours_id : Vec = knn_neighbours_for_tests[i].iter().map(|p| p.d_id).collect(); + let knn_neighbours_dist : Vec = knn_neighbours_for_tests[i].iter().map(|p| p.distance).collect(); + nb_returned.push(knn_neighbours_dist.len()); + // count how many distances of knn_neighbours_dist are less than + let recall = knn_neighbours_dist.iter().filter(|x| *x <= &max_dist).count(); + recalls.push(recall); + let mut ratio = 0.; + if knn_neighbours_dist.len() >= 1 { + ratio = knn_neighbours_dist[knn_neighbours_dist.len()-1]/max_dist; + } + last_distances_ratio.push(ratio); + } + let mean_recall = (recalls.iter().sum::() as f32)/((knbn * recalls.len()) as f32); + println!("\n mean fraction nb returned by search {:?} ", (nb_returned.iter().sum::() as f32)/ ((nb_returned.len() * knbn) as f32)); + println!("\n last distances ratio {:?} ", last_distances_ratio.iter().sum::() / last_distances_ratio.len() as f32); + println!("\n recall rate for {:?} is {:?} , nb req /s {:?}", anndata.fname, mean_recall, (nb_search as f32)*1.0e+6_f32/search_sys_time); +} diff --git a/examples/random.rs b/examples/random.rs new file mode 100644 index 0000000..e4d53b3 --- /dev/null +++ b/examples/random.rs @@ -0,0 +1,66 @@ +extern crate hnsw_rs; + + +use std::time::{Duration, SystemTime}; +use cpu_time::ProcessTime; +use rand::distributions::{Uniform}; +use rand::prelude::*; + +use hnsw_rs::prelude::*; + +//use crate::point::*; + +use hnsw_rs::test::*; + + + +fn main() { + let _res = simple_logger::init(); + // + let nb_elem = 500000; + let dim = 25; + let data = gen_random_matrix_f32(dim, nb_elem); + let data_with_id = data.iter().zip(0..data.len()).collect(); + + let ef_c = 200; + let max_nb_connection = 15; + let nb_layer = 16.min((nb_elem as f32).ln().trunc() as usize); + let hns = Hnsw::::new(max_nb_connection, nb_elem, nb_layer, ef_c, DistL2{}); + let mut start = ProcessTime::now(); + let mut begin_t = SystemTime::now(); + hns.parallel_insert(&data_with_id); + let mut cpu_time: Duration = start.elapsed(); + println!(" hnsw data insertion cpu time {:?}", cpu_time); + println!(" hnsw data insertion parallel, system time {:?} \n", begin_t.elapsed().unwrap()); + hns.dump_layer_info(); + println!(" parallel hnsw data nb point inserted {:?}", hns.get_nb_point()); + // + // serial insertion + // + let hns = Hnsw::::new(max_nb_connection, nb_elem, nb_layer, ef_c, DistL2{}); + start = ProcessTime::now(); + begin_t = SystemTime::now(); + for _i in 0..data_with_id.len() { + hns.insert(data_with_id[_i]); + } + cpu_time = start.elapsed(); + println!("\n\n serial hnsw data insertion {:?}", cpu_time); + println!(" hnsw data insertion serial, system time {:?}", begin_t.elapsed().unwrap()); + hns.dump_layer_info(); + println!(" serial hnsw data nb point inserted {:?}", hns.get_nb_point()); + + let ef_search = max_nb_connection * 2; + let knbn = 10; + // + for _iter in 0..100 { + let mut r_vec = Vec::::with_capacity(dim); + let mut rng = thread_rng(); + let unif = Uniform::::new(0.,1.); + for _ in 0..dim { + r_vec.push(rng.sample(unif)); + } + // + let _neighbours = hns.search(&r_vec, knbn, ef_search); + } + +} \ No newline at end of file diff --git a/src/annhdf5.rs b/src/annhdf5.rs new file mode 100644 index 0000000..efd604d --- /dev/null +++ b/src/annhdf5.rs @@ -0,0 +1,217 @@ +// +//! This file provides hdf5 utilities to load ann-benchmarks hdf5 data files + +use ndarray::{Array2}; +use ::hdf5::*; + + + +use crate::hnsw; + +#[allow(unused_imports)] // necessary for rls +use crate::dist; + +pub use self::hnsw::*; + +// datasets +// . distances (nbojects, dim) f32 matrix for tests objects +// . neighbors (nbobjects, nbnearest) int32 matrix giving the num of nearest neighbors in train data +// . test (nbobjects, dim) f32 matrix test data +// . train (nbobjects, dim) f32 matrix train data + +/// a structure to load hdf5 data file benchmarks from https://github.com/erikbern/ann-benchmarks +pub struct AnnBenchmarkData { + pub fname: String, + /// distances from each test object to its nearest neighbours. + pub test_distances : Array2, + // for each test data , id of its nearest neighbours + pub test_neighbours: Array2, + /// list of vectors for which we will search ann. + pub test_data: Vec>, + /// list of data vectors and id + pub train_data: Vec<(Vec, usize)>, + /// searched results. first neighbours for each test data. + pub searched_neighbours : Vec>, + /// distances of neighbours obtained of each test + pub searched_distances: Vec>, +} + +impl AnnBenchmarkData { + pub fn new(fname: String) -> Result { + let res = hdf5::File::open(fname.clone(), "r"); + if res.is_err() { + println!("you must download file {:?}", fname); + panic!("download benchmark file some where and modify examples source file accordingly"); + } + let file = res.ok().unwrap(); + // + // get test distances + // + let res_distances = file.dataset("distances"); + if res_distances.is_err() { + // let reader = hdf5::Reader::::new(&test_distance); + panic!("error getting distances dataset"); + } + let distances = res_distances.unwrap(); + let shape = distances.shape(); + assert_eq!(shape.len(), 2); + let dataf32 = distances.dtype().unwrap().is::(); + if !dataf32 { + // error + panic!("error getting type distances dataset"); + } + // read really data + let res = distances.read_2d:: (); + if res.is_err() { + // some error + panic!("error reading distances dataset"); + } + let test_distances = res.unwrap(); + // a check for row order + log::debug!(" first 2 distances for first test {:?} {:?} ", test_distances.get((0,0)).unwrap(), test_distances.get((0,1)).unwrap()); + // + // read neighbours + // + let res_neighbours = file.dataset("neighbors"); + if res_neighbours.is_err() { + // let reader = hdf5::Reader::::new(&test_distance); + panic!("error getting neighbours"); + } + let neighbours = res_neighbours.unwrap(); + let shape = neighbours.shape(); + assert_eq!(shape.len(), 2); + println!("neighbours shape : {:?}", shape); + let datai32 = neighbours.dtype().unwrap().is::(); + if !datai32 { + // error + panic!("error getting type neighbours"); + } + // read really data + let res = neighbours.read_2d:: (); + if res.is_err() { + // some error + panic!("error reading neighbours dataset"); + } + let test_neighbours = res.unwrap(); + log::debug!(" first 2 neighbours for first test {:?} {:?} ", test_neighbours.get((0,0)).unwrap(), test_neighbours.get((0,1)).unwrap()); + println!("\n 10 first neighbours for first vector : "); + for i in 0..10 { + print!(" {:?} ", test_neighbours.get((0,i)).unwrap()); + } + println!("\n 10 first neighbours for second vector : "); + for i in 0..10 { + print!(" {:?} ", test_neighbours.get((1,i)).unwrap()); + } + // + // read test data + // =============== + // + let res_testdata = file.dataset("test"); + if res_testdata.is_err() { + panic!("error getting test de notataset"); + } + let test_data = res_testdata.unwrap(); + let shape = test_data.shape(); // nota shape returns a slice, dim returns a t-uple + assert_eq!(shape.len(), 2); + let dataf32 = test_data.dtype().unwrap().is::(); + if !dataf32 { + panic!("error getting type de notistances dataset"); + } + // read really datae not + let res = test_data.read_2d:: (); + if res.is_err() { + // some error + panic!("error reading distances dataset"); + } + let test_data_2d = res.unwrap(); + let mut test_data = Vec::>::with_capacity(shape[1]); + let (nbrow, nbcolumn) = test_data_2d.dim(); + println!(" test data, nb element {:?}, dim : {:?}", nbrow, nbcolumn); + for i in 0..nbrow { + let mut vec = Vec::with_capacity(nbcolumn); + for j in 0..nbcolumn { + vec.push(*test_data_2d.get((i,j)).unwrap()); + } + test_data.push(vec); + } + // + // loaf train data + // + let res_traindata = file.dataset("train"); + if res_traindata.is_err() { + panic!("error getting distances dataset"); + } + let train_data = res_traindata.unwrap(); + let train_shape = train_data.shape(); + assert_eq!(shape.len(), 2); + if test_data_2d.dim().1 != train_shape[1] { + println!("test and train have not the same dimension"); + panic!(); + } + println!("\n train data shape : {:?}, nbvector {:?} ", train_shape, train_shape[0]); + let dataf32 = train_data.dtype().unwrap().is::(); + if !dataf32 { + // error + panic!("error getting type distances dataset"); + } + // read really data + let res = train_data.read_2d:: (); + if res.is_err() { + // some error + panic!("error reading distances dataset"); + } + let train_data_2d = res.unwrap(); + let mut train_data = Vec::<(Vec,usize)>::with_capacity(shape[1]); + let (nbrow, nbcolumn) = train_data_2d.dim(); + for i in 0..nbrow { + let mut vec = Vec::with_capacity(nbcolumn); + for j in 0..nbcolumn { + vec.push(*train_data_2d.get((i,j)).unwrap()); + } + train_data.push((vec,i)); + } + // + // now allocate array's for result + // + println!(" allocating vector for search neighbours answer : {:?}", test_data.len()); + let searched_neighbours = Vec::>::with_capacity(test_data.len()); + let searched_distances = Vec::>::with_capacity(test_data.len()); + // searched_distances + Ok(AnnBenchmarkData{fname:fname.clone(), test_distances, test_neighbours, test_data, train_data, searched_neighbours, searched_distances}) + } // end new + + /// do l2 normalisation of test and train vector to use DistDot metrinc instead DistCosine to spare cpu + pub fn do_l2_normalization(&mut self) { + for i in 0..self.test_data.len() { + dist::l2_normalize(&mut self.test_data[i]); + } + for i in 0..self.train_data.len() { + dist::l2_normalize(&mut self.train_data[i].0); + } + } +} // end of impl block + + + + +#[cfg(test)] + + +mod tests { + +use super::*; + + +#[test] + +fn test_load_hdf5() { + let _res = simple_logger::init(); + // + let fname = String::from("/home.2/Data/ANN/glove-25-angular.hdf5"); + println!("\n\n test_load_hdf5 {:?}", fname); + // now recall that data are stored in row order. + let _anndata = AnnBenchmarkData::new(fname).unwrap(); + // +} // end of test_load_hdf5 + +} // end of module test \ No newline at end of file diff --git a/src/api.rs b/src/api.rs new file mode 100644 index 0000000..3f8b43e --- /dev/null +++ b/src/api.rs @@ -0,0 +1,126 @@ +//! Api for external language. +//! This file provides a trait to be used as an opaque pointer for C or Julia calls used in file libext.rs + + +use std::io::prelude::*; +use std::fs::OpenOptions; +use std::io::{BufWriter}; +use std::path::{PathBuf}; + +use typename::TypeName; + +use crate::hnsw::*; +use crate::hnswio::*; + + +pub trait AnnT { + /// type of data vectors + type Val; + /// + fn insert_data(&mut self, data: &Vec, id: usize); + /// + fn search_neighbours(&self, data :&Vec , knbn : usize, ef_s: usize) -> Vec; + /// + fn parallel_insert_data(&mut self, data: &Vec<(&Vec, usize)> ); + /// + fn parallel_search_neighbours(&self, data: &Vec >, knbn : usize, ef_s: usize) -> Vec> ; + /// + /// dumps a data and graph in 2 files. + /// Datas are dumped in file filename.hnsw.data and graph in filename.hnsw.graph + fn file_dump(&self, filename: &String) -> Result; + } + + +impl AnnT for Hnsw where T:Copy+Clone+Send+Sync+TypeName , D: Distance+TypeName+Send+Sync { + type Val= T; + /// + fn insert_data(&mut self, data: &Vec, id: usize) { + self.insert((data, id)); + } + /// + fn search_neighbours(&self, data : &Vec, knbn : usize, ef_s: usize) -> Vec { + self.search(data, knbn, ef_s) + } + fn parallel_insert_data(&mut self, data: &Vec<(&Vec, usize)> ) { + self.parallel_insert(data); + } + + fn parallel_search_neighbours(&self, data: &Vec>, knbn : usize, ef_s: usize) -> Vec> { + self.parallel_search(data, knbn, ef_s) + } + /// + fn file_dump(&self, filename: &String) -> Result { + let mut graphname = filename.clone(); + graphname.push_str(".hnsw.graph"); + let graphpath = PathBuf::from(graphname); + let fileres = OpenOptions::new().write(true).create(true).truncate(true).open(&graphpath); + if fileres.is_err() { + println!("could not open file {:?}", graphpath.as_os_str()); + return Err("could not open file".to_string()); + } + let graphfile = fileres.unwrap(); + // + let mut dataname = filename.clone(); + dataname.push_str(".hnsw.data"); + let datapath = PathBuf::from(dataname); + let fileres = OpenOptions::new().write(true).create(true).truncate(true).open(&datapath); + if fileres.is_err() { + println!("could not open file {:?}", datapath.as_os_str()); + return Err("could not open file".to_string()); + } + let datafile = fileres.unwrap(); + let mut graphbufw = BufWriter::with_capacity(10000000 , graphfile); + let mut databufw = BufWriter::with_capacity(10000000 , datafile); + let res = self.dump(DumpMode::Full, &mut graphbufw, &mut databufw); + graphbufw.flush().unwrap(); + databufw.flush().unwrap(); + log::debug!("\n end of dump"); + return res; + } // end of dump +} // end of impl block AnnT for Hnsw + + +// macro export makes the macro export t the root of the crate +#[macro_export] +macro_rules! mapdist_t( + ("DistL1") => (crate::dist::DistL1); + ("DistL2") => (crate::dist::DistL2); + ("DistL2") => (crate::dist::DistL2); + ("DistDot") => (crate::dist::DistDot); + ("DistHamming") => (crate::dist::DistHamming); + ("DistJaccard") => (crate::dist::DistJaccard); + ("DistPtr") => (crate::dist::DistPtr); +); + + +#[macro_export] +macro_rules! mapdist_c( + ("DistL1") => (crate::dist::DistL1{}); +); + + +#[macro_export] +macro_rules! genhnsw( + ($val:ty, $dist:tt, $max_nb_conn:tt, $ef_const:tt) => ( + Hnsw::<$val, mapdist_t!($dist) >::new($max_nb_conn, 10000, 16, $ef_const, mapdist_c!($dist)) + ) +); + + + + + +////////////////////////////////////////////////////////////////////////////////////////// + +#[cfg(test)] +mod tests { + +use super::*; + + +#[test] + fn test_genhnsw() { + let h = genhnsw!(f32, "DistL1", 24, 48); + println!("test constructed Hnsw with distance : {:?} ", h.get_distance_name()); + } +} // end of module test diff --git a/src/dist.rs b/src/dist.rs new file mode 100644 index 0000000..973e052 --- /dev/null +++ b/src/dist.rs @@ -0,0 +1,869 @@ +//! Some standard distances as L1, L2, Cosine, Jaccard, Hamming +//! and a structure to enable the user to implement its own distances. +//! For the heavily used case (f32) we provide simd avx2 implementation. + +#![allow(dead_code)] +#![allow(unused_macros)] + + +extern crate simdeez; +use simdeez::avx2::*; +use simdeez::sse2::*; +use simdeez::*; + + +/// The trait describing distance. +/// For example for the L1 distance +/// +/// pub struct DistL1; +/// +/// implement Distance for DistL1 { +/// } +/// +/// +/// The L1 and Cosine distance are implemented for u16, i32, i64, f32, f64 +/// +/// + + +use typename::TypeName; + +use std::os::raw::*; + +enum DistKind { + DistL1(String), + DistL2(String), + /// This is the same as Cosine dist but all data L2-normalized to 1. + DistDot(String), + DistCosine(String), + DistHamming(String), + DistJaccard(String), + DistHellinger(String), + DistJeffreys(String), + /// To store a distance defined by a C pointer function + DistPtr(String), +} + + +/// This is the basic Trait describing a distance. The structure Hnsw can be instantiated by anything +/// satisfying this Trait. The crate provides for L1, L2 , Cosine, Jaccard, Hamming. +/// For other distance the structure `struct DistFn` +/// can use any closure satisfying the trait Fn(&[T], &[T]) -> f32> +pub trait Distance { + fn eval(&self, va:&[T], vb: &[T]) -> f32; +} + +#[derive(TypeName, Default)] +pub struct DistL1; + + +macro_rules! implementL1Distance ( + ($ty:ty) => ( + + impl Distance<$ty> for DistL1 { + fn eval(&self, va:&[$ty], vb: &[$ty]) -> f32 { + // RUSTFLAGS = "-C opt-level=3 -C target-cpu=native" + va.iter().zip(vb.iter()).map(|t| (*t.0 as f32- *t.1 as f32).abs()).sum() + } // end of compute + } // end of impl block + ) // end of pattern matching +); + + +implementL1Distance!(i32); +implementL1Distance!(f64); +implementL1Distance!(i64); +implementL1Distance!(u32); +implementL1Distance!(u16); +implementL1Distance!(u8); + + +unsafe fn distance_l1_f32 (va:&[f32], vb: &[f32]) -> f32 { + let mut dist = 0_f32; + // + let nb_simd = va.len() / S::VF32_WIDTH; + let simd_length = nb_simd * S::VF32_WIDTH; + let mut i = 0; + while i < simd_length { + let a = S::loadu_ps(&va[i]); + let b = S::loadu_ps(&vb[i]); + let delta = S::abs_ps(a-b); + dist += S::horizontal_add_ps(delta); + // + i += S::VF32_WIDTH; + } + for i in simd_length..va.len() { +// log::debug!("distance_l1_f32, i {:?} len {:?} nb_simd {:?} VF32_WIDTH {:?}", i, va.len(), nb_simd, S::VF32_WIDTH); + dist += (va[i]- vb[i]).abs(); + } + assert!( dist >= 0.); + dist +} // end of distance_l1_f32 + + + #[target_feature(enable = "avx2")] +unsafe fn distance_l1_f32_avx2(va:&[f32], vb: &[f32]) -> f32 { + distance_l1_f32::(va,vb) +} + + + + +impl Distance for DistL1 { + fn eval(&self, va:&[f32], vb: &[f32]) -> f32 { + // + // assert_eq!(va.len(), vb.len()); + // + if is_x86_feature_detected!("avx2") { + unsafe {distance_l1_f32_avx2(va,vb)} + } + else { + va.iter().zip(vb.iter()).map(|t| (*t.0 as f32- *t.1 as f32).abs()).sum() + } + } + +} +//======================================================================== + +#[derive(TypeName, Default)] +pub struct DistL2; + + +macro_rules! implementL2Distance ( + ($ty:ty) => ( + + impl Distance<$ty> for DistL2 { + fn eval(&self, va:&[$ty], vb: &[$ty]) -> f32 { + let norm : f32 = va.iter().zip(vb.iter()).map(|t| (*t.0 as f32- *t.1 as f32) * (*t.0 as f32- *t.1 as f32)).sum(); + norm.sqrt() + } // end of compute + } // end of impl block + ) // end of pattern matching +); + + +//implementL2Distance!(f32); +implementL2Distance!(i32); +implementL2Distance!(f64); +implementL2Distance!(i64); +implementL2Distance!(u32); +implementL2Distance!(u16); +implementL2Distance!(u8); + + + + +unsafe fn distance_l2_f32 (va:&[f32], vb: &[f32]) -> f32 { + let mut dist = 0_f32; + // + let nb_simd = va.len() / S::VF32_WIDTH; + let simd_length = nb_simd * S::VF32_WIDTH; + let mut i = 0; + while i < simd_length { + let a = S::loadu_ps(&va[i]); + let b = S::loadu_ps(&vb[i]); + let mut delta = a-b; + delta *= delta; + dist += S::horizontal_add_ps(delta); + // + i += S::VF32_WIDTH; + } + for i in simd_length..va.len() { + dist += (va[i]- vb[i]) * (va[i]- vb[i]); + } + assert!( dist >= 0.); + dist.sqrt() +} // end of distance_l2_f32 + + + #[target_feature(enable = "avx2")] +unsafe fn distance_l2_f32_avx2(va:&[f32], vb: &[f32]) -> f32 { + distance_l2_f32::(va,vb) +} + + + +impl Distance for DistL2 { + fn eval(&self, va:&[f32], vb: &[f32]) -> f32 { + // + if is_x86_feature_detected!("avx2") { + unsafe {distance_l2_f32_avx2(va,vb)} + } + else { + let norm : f32 = va.iter().zip(vb.iter()).map(|t| (*t.0 as f32- *t.1 as f32) * (*t.0 as f32- *t.1 as f32)).sum(); + assert!(norm >= 0.); + norm.sqrt() + } + } + +} + +//========================================================================= + +#[derive(TypeName, Default)] +pub struct DistCosine; + + +macro_rules! implementCosDistance( + ($ty:ty) => ( + impl Distance<$ty> for DistCosine { + fn eval(&self, va:&[$ty], vb: &[$ty]) -> f32 { + let dist:f32; + let zero:f32 = 0.; + // to // by rayon + let res = va.iter().zip(vb.iter()).map(|t| ((*t.0 * *t.1) as f32, (*t.0 * *t.0) as f32, (*t.1 * *t.1) as f32)). + fold((0., 0., 0.), |acc , t| (acc.0 + t.0, acc.1 + t.1, acc.2 + t.2)); + // + if res.1 > zero && res.2 > zero { + dist = 1. - res.0 / (res.1 * res.2).sqrt(); + } + else { + dist = 0.; + } + // + return dist; + } // end of function + } // end of impl block + ) // end of matching +); + + +implementCosDistance!(f32); +implementCosDistance!(f64); +implementCosDistance!(i64); +implementCosDistance!(i32); +implementCosDistance!(u16); + + + +//========================================================================= + +/// This is essentially the Cosine distance but we suppose +/// all vectors (graph construction and request vectors have been l2 normalized to unity +/// BEFORE INSERTING in HNSW!. +/// No control is made, so it is the user responsability to send normalized vectors +/// everywhere in inserting and searching. +/// +/// In large dimensions (hundreds) this pre-normalization spare cpu time. +/// At low dimensions (a few ten's there is not a significant gain). +/// This distance makes sense only for f16, f32 or f64 +/// We provide for avx2 implementations for f32 that provides consequent gains +/// in large dimensions + +#[derive(TypeName, Default)] +pub struct DistDot; + + + +macro_rules! implementDotDistance( + ($ty:ty) => ( + impl Distance<$ty> for DistDot { + fn eval(&self, va:&[$ty], vb: &[$ty]) -> f32 { + let zero:f32 = 0f32; + // to // by rayon + let dot = va.iter().zip(vb.iter()).map(|t| (*t.0 * *t.1) as f32).fold(0., |acc , t| (acc + t)); + // + assert(dot <= 1.); + return 1. - dot; + } // end of function + } // end of impl block + ) // end of matching +); + + + +unsafe fn distance_dot_f32 (va:&[f32], vb: &[f32]) -> f32 { + let mut dot = 0_f32; + // + let mut i = 0; + let nb_simd = va.len() / S::VF32_WIDTH; + let simd_length = nb_simd * S::VF32_WIDTH; + while i < simd_length { + let a = S::loadu_ps(&va[i]); + let b = S::loadu_ps(&vb[i]); + let delta = a*b; + dot += S::horizontal_add_ps(delta); + // + i += S::VF32_WIDTH; + } + for i in simd_length..va.len() { + dot += va[i]*vb[i]; + } + assert!(dot <= 1.000002); + (1. - dot).max(0.) +} // end of distance_dot_f32 + + + + + #[target_feature(enable = "avx2")] +unsafe fn distance_dot_f32_avx2(va:&[f32], vb: &[f32]) -> f32 { + distance_dot_f32::(va,vb) +} + + + #[target_feature(enable = "sse2")] +unsafe fn distance_dot_f32_sse2(va:&[f32], vb: &[f32]) -> f32 { + distance_dot_f32::(va,vb) +} + + +impl Distance for DistDot { + fn eval(&self, va:&[f32], vb: &[f32]) -> f32 { + // + if is_x86_feature_detected!("avx2") { + unsafe { distance_dot_f32_avx2(va,vb) } + } + else if is_x86_feature_detected!("sse2") { + unsafe { distance_dot_f32_sse2(va,vb) } + } + else { + let dot = 1. - va.iter().zip(vb.iter()).map(|t| (*t.0 * *t.1) as f32).fold(0., |acc , t| (acc + t)); + assert!( dot >= 0.); + dot + } + } +} + +pub fn l2_normalize(va:& mut [f32]) { + let l2norm = va.iter().map(|t| (*t * *t) as f32).sum::().sqrt(); + if l2norm > 0. { + for i in 0..va.len() { + va[i] = va[i]/l2norm; + } + } +} + + +//======================================================================================= + +/// +/// A structure to compute Hellinger distance between probalilities. +/// Vector must be >= 0 and normalized to 1. +/// +/// The distance computation does not check that +/// and in fact simplifies the expression of distance assuming vectors are positive and L1 normalised to 1. +/// The user must enforce these conditions before inserting otherwise results will be meaningless +/// at best or code will panic! +/// +/// For f32 a simd implementation is provided if avx2 is detected. +#[derive(TypeName, Default)] +pub struct DistHellinger; + + +// default implementation +macro_rules! implementHellingerDistance ( + ($ty:ty) => ( + + impl Distance<$ty> for DistHellinger { + fn eval(&self, va:&[$ty], vb: &[$ty]) -> f32 { + // RUSTFLAGS = "-C opt-level=3 -C target-cpu=native" + // to // by rayon + let mut dist = va.iter().zip(vb.iter()).map(|t| ((*t.0).sqrt() * (*t.1).sqrt()) as f32).fold(0., |acc , t| (acc + t*t)); + dist = (1. - dist).sqrt(); + dist + } // end of compute + } // end of impl block + ) // end of pattern matching +); + + +implementHellingerDistance!(f64); + + +unsafe fn distance_hellinger_f32 (va:&[f32], vb: &[f32]) -> f32 { + let mut dist = 0_f32; + // + let mut i = 0; + let nb_simd = va.len() / S::VF32_WIDTH; + let simd_length = nb_simd * S::VF32_WIDTH; + while i < simd_length { + let a = S::loadu_ps(&va[i]); + let b = S::loadu_ps(&vb[i]); + let prod = a * b; + let prod_s = S::sqrt_ps(prod); + dist += S::horizontal_add_ps(prod_s); + // + i += S::VF32_WIDTH; + } + for i in simd_length..va.len() { + dist += va[i].sqrt() * vb[i].sqrt(); + } + assert!(1. - dist >= -0.000001); + dist = (1. - dist).max(0.).sqrt(); + dist +} // end of distance_hellinger_f32 + + + +#[target_feature(enable = "avx2")] +unsafe fn distance_hellinger_f32_avx2(va:&[f32], vb: &[f32]) -> f32 { + distance_hellinger_f32::(va,vb) +} + + + +impl Distance for DistHellinger { + fn eval(&self, va:&[f32], vb: &[f32]) -> f32 { + // + if is_x86_feature_detected!("avx2") { + unsafe { distance_hellinger_f32_avx2(va,vb) } + } + else { + let mut dist = va.iter().zip(vb.iter()).map(|t| ((*t.0).sqrt() * (*t.1).sqrt()) as f32).fold(0., |acc , t| (acc + t*t)); + // if too far away from >= panic else reset! + assert!(1. - dist >= -0.000001); + dist = (1. - dist).max(0.).sqrt(); + dist + } + } +} + + +//======================================================================================= + + +/// +/// A structure to compute Jeffreys distance between probalilities. +/// If p and q 2 probability distributions +/// the distance is computed as: +/// sum (p[i] - q[i]) * ln(p[i]/q[i]) +/// +/// To take care of null probabilities in the formula we use max(x[i],1.E-30) +/// for x = p and q in the log compuations +/// +/// Vector must be >= 0 and normalized to 1! +/// The distance computation does not check that, +/// The user must enforce these conditions before inserting in the hnws structure, +/// otherwise results will be meaningless at best or code will panic! +/// +/// For f32 a simd implementation is provided if avx2 is detected. +#[derive(TypeName, Default)] +pub struct DistJeffreys; + + +const M_MIN:f32 = 1.0e-30; + + +// default implementation +macro_rules! implementJeffreysDistance ( + ($ty:ty) => ( + + impl Distance<$ty> for DistJeffreys { + fn eval(&self, va:&[$ty], vb: &[$ty]) -> f32 { + // RUSTFLAGS = "-C opt-level=3 -C target-cpu=native" + let dist = va.iter().zip(vb.iter()).map(|t| (*t.0 - *t.1) * ((*t.0).max(M_MIN as f64)/ (*t.1).max(M_MIN as f64)).ln() as f64).fold(0., |acc , t| (acc + t*t)); + dist as f32 + } // end of compute + } // end of impl block + ) // end of pattern matching +); + + +implementJeffreysDistance!(f64); + + + +unsafe fn distance_jeffreys_f32 (va:&[f32], vb: &[f32]) -> f32 { + let mut dist = 0_f32; + // + let mut i = 0; + let mut logslice = Vec::::with_capacity(S::VF32_WIDTH as usize); + let nb_simd = va.len() / S::VF32_WIDTH; + let simd_length = nb_simd * S::VF32_WIDTH; + while i < simd_length { + let a = S::loadu_ps(&va[i]); + let b = S::loadu_ps(&vb[i]); + let delta = a - b; + for j in 0..S::VF32_WIDTH { + // take care of zeros! + logslice.push((va[i+j].max(M_MIN)/vb[i+j].max(M_MIN)).ln()); + } + let prod_s = delta * S::loadu_ps(&logslice.as_slice()[0]); + dist += S::horizontal_add_ps(prod_s); + // + i += S::VF32_WIDTH; + } + for i in simd_length..va.len() { + if vb[i] > 0. { + dist += (va[i] - vb[i]) * (va[i].max(M_MIN) / vb[i].max(M_MIN)).ln(); + } + } + dist +} // end of distance_hellinger_f32 + + + + +#[target_feature(enable = "avx2")] +unsafe fn distance_jeffreys_f32_avx2(va:&[f32], vb: &[f32]) -> f32 { + distance_jeffreys_f32::(va,vb) +} + + +impl Distance for DistJeffreys { + fn eval(&self, va:&[f32], vb: &[f32]) -> f32 { + // + if is_x86_feature_detected!("avx2") { + unsafe { distance_jeffreys_f32_avx2(va,vb) } + } + else { + let dist = va.iter().zip(vb.iter()).map(|t| (*t.0 - *t.1) * ((*t.0).max(M_MIN)/ (*t.1).max(M_MIN)).ln() as f32).fold(0., |acc , t| (acc + t*t)); + dist + } + } +} + + +//======================================================================================= + + +#[derive(TypeName, Default)] +pub struct DistHamming; + + + +macro_rules! implementHammingDistance ( + ($ty:ty) => ( + + impl Distance<$ty> for DistHamming { + fn eval(&self, va:&[$ty], vb: &[$ty]) -> f32 { + // RUSTFLAGS = "-C opt-level=3 -C target-cpu=native" + let norm : f32 = va.iter().zip(vb.iter()).filter(|t| t.0 != t.1).count() as f32; + norm + } // end of compute + } // end of impl block + ) // end of pattern matching +); + + +implementHammingDistance!(u8); +implementHammingDistance!(u16); +implementHammingDistance!(u32); + + +//==================================================================================== +// Jaccard Distance + + +#[derive(TypeName, Default)] +pub struct DistJaccard; + + +// contruct a 2-uple accumulator that has sum of max in first component , and sum of min in 2 component +// stay in integer as long as possible + +macro_rules! implementJaccardDistance ( + ($ty:ty) => ( + + impl Distance<$ty> for DistJaccard { + fn eval(&self, va:&[$ty], vb: &[$ty]) -> f32 { + let (max,min) : (u32, u32) = va.iter().zip(vb.iter()).fold((0u32,0u32), |acc, t| if t.0 > t.1 { + (acc.0 + *t.0 as u32, acc.1 + *t.1 as u32) } + else { + (acc.0 + *t.1 as u32 , acc.1 + *t.0 as u32) + } + ); + if max > 0 { + 1. - (min as f32)/ (max as f32) + } + else { + 0. + } + } // end of compute + } // end of impl block + ) // end of pattern matching +); + + +implementJaccardDistance!(u8); +implementJaccardDistance!(u16); +implementJaccardDistance!(u32); + + +//======================================================================================= +// Case of function pointers (cover Trait Fn , FnOnce ...) +// The book (Function item types): " There is a coercion from function items to function pointers with the same signature " +// The book (Call trait and coercions): "Non capturing closures can be coerced to function pointers with the same signature" + + +/// This type is for function with a C-API +/// Distances can be computed by such a function. It +/// takes as arguments the two (C, rust, julia) pointers to primitive type vectos and length +/// passed as a unsignedlonlong (64 bits) (which is called c_ulong in Rust!) and Culonglong in Julia +/// +type DistCFnPtr = extern "C" fn(*const T, *const T, len : c_ulong) -> f32; + + +/// A structure to implement Distance Api for type DistCFnPtr +/// It must be noted that this can be used in Julia via the macro @cfunction +/// to define interactiveley a distance function , compile it on the fly and sent it +/// to Rust via the init_hnsw_{f32, i32, u16, u32, u8} function +/// defined in libext +/// +#[derive(TypeName)] +pub struct DistCFFI { + dist_function : DistCFnPtr, +} + +impl DistCFFI { + pub fn new(f:DistCFnPtr) -> Self { + DistCFFI{ dist_function:f} + } +} + +impl Distance for DistCFFI { + fn eval(&self, va:&[T], vb: &[T]) -> f32 { + // get pointers + let len = va.len(); + let ptr_a = va.as_ptr(); + let ptr_b = vb.as_ptr(); + let dist = (self.dist_function)(ptr_a, ptr_b, len as c_ulong); + log::trace!("DistCFFI dist_function_ptr {:?} returning {:?} ", self.dist_function, dist); + dist + } // end of compute + } // end of impl block + + +//======================================================================================================== + + +/// This structure is to let user define their own distance with closures. +/// An example is given in tests. +#[derive(TypeName)] +pub struct DistFn { + dist_function : Box f32>, +} + +impl DistFn { + /// construction of a DistFn + pub fn new(f : Box f32>) -> Self { + DistFn{ dist_function : f } + } + +} + +impl Distance for DistFn { + fn eval(&self, va:&[T], vb: &[T]) -> f32 { + (self.dist_function)(va,vb) + } +} + +//======================================================================================= + +#[cfg(test)] + +mod tests { +use super::*; + + +#[test] +fn test_access_to_dist_l1() { + let distl1 = DistL1; + // + let v1: Vec = vec![1, 2, 3]; + let v2: Vec = vec![2, 2, 3]; + + let d1 = Distance::eval(&distl1, &v1,&v2); + assert_eq!(d1, 1 as f32); + + let v3: Vec = vec![1., 2., 3.]; + let v4: Vec = vec![2., 2., 3.]; + let d2 = distl1.eval(&v3,&v4); + assert_eq!(d2, 1 as f32); + + +} + + +#[test] + +fn have_avx2() { + if is_x86_feature_detected!("avx2") { + println!("I have avx2"); + } + else { + println!(" ************ I DO NOT have avx2 ***************"); + } +} + +#[test] + +fn have_sse2() { + if is_x86_feature_detected!("sse2") { + println!("I have sse2"); + } + else { + println!(" ************ I DO NOT have SSE2 ***************"); + } +} + + +#[test] +fn test_access_to_dist_cos() { + let distcos = DistCosine; + // + let v1: Vec = vec![1,-1, 1]; + let v2: Vec = vec![2, 1 , -1]; + + let d1 = Distance::eval(&distcos, &v1, &v2); + assert_eq!(d1, 1. as f32); + // + let v1: Vec = vec![1.234,-1.678, 1.367]; + let v2: Vec = vec![4.234,-6.678, 10.367]; + let d1 = Distance::eval(&distcos, &v1,&v2); + + let mut normv1 = 0.; + let mut normv2 = 0.; + let mut prod = 0.; + for i in 0..v1.len() { + prod += v1[i] * v2[i]; + normv1 += v1[i]*v1[i]; + normv2 += v2[i] * v2[i]; + } + let dcos = 1. - prod/(normv1*normv2).sqrt(); + println!("dist cos avec macro = {:?} , avec for {:?}", d1 , dcos); +} + +#[test] +fn test_dot_distances () { + let mut v1: Vec = vec![1.234,-1.678, 1.367]; + let mut v2: Vec = vec![4.234,-6.678, 10.367]; + + let mut normv1 = 0.; + let mut normv2 = 0.; + let mut prod = 0.; + for i in 0..v1.len() { + prod += v1[i] * v2[i]; + normv1 += v1[i]*v1[i]; + normv2 += v2[i] * v2[i]; + } + let dcos = 1. - prod/(normv1*normv2).sqrt(); + // + l2_normalize(&mut v1); + l2_normalize(&mut v2); + + println!( " after normalisation v1 = {:?}" , v1); + + let dot = DistDot.eval(&v1,&v2); + + println!("dot cos avec prenormalisation = {:?} , avec for {:?}", dot , dcos); +} + +#[test] +fn test_jaccard_u16() { + let v1: Vec = vec![1, 2 , 1, 4, 3]; + let v2: Vec = vec![2, 2 , 1, 5, 6]; + + let dist = DistJaccard.eval(&v1, &v2); + println!("dist jaccard = {:?}", dist); + assert_eq!(dist, 1. - 11./16.); +} // end of test_jaccard + + + +extern "C" fn dist_func_float(va : *const f32, vb : *const f32, len : c_ulonglong) -> f32 { + let mut dist : f32 = 0.; + let sa = unsafe {std::slice::from_raw_parts(va, len as usize) }; + let sb = unsafe { std::slice::from_raw_parts(vb, len as usize) }; + + for i in 0..len { + dist += (sa[i as usize] - sb[i as usize]).abs().sqrt(); + } + dist +} + + +#[test] + +fn test_dist_ext_float() { + let va : Vec:: = vec! [1. , 2., 3.]; + let vb : Vec:: = vec! [1. , 2., 3.]; + println!("in test_dist_ext_float"); + let dist1 = dist_func_float(va.as_ptr(), vb.as_ptr(), va.len() as c_ulong); + println!("test_dist_ext_float computed : {:?}", dist1); + + let mydist = DistCFFI::::new(dist_func_float); + + let dist2 = mydist.eval(&va, &vb); + assert_eq!(dist1, dist2); +} // end test_dist_ext_float + + +#[test] + +fn test_my_closure() { + let weight = vec![0.1, 0.8, 0.1]; + let my_fn = move | va : &[f32] , vb: &[f32] | -> f32 { + // should check that we work with same size for va, vb, and weight... + let mut dist : f32 = 0.; + for i in 0..va.len() { + dist += weight[i] * (va[i] - vb[i]).abs(); + } + dist + }; + let my_boxed_f = Box::new(my_fn); + let my_boxed_dist = DistFn::::new(my_boxed_f); + let va : Vec:: = vec! [1. , 2., 3.]; + let vb : Vec:: = vec! [2. , 2., 4.]; + let dist = my_boxed_dist.eval(&va, &vb); + println!("test_my_closure computed : {:?}", dist); + assert_eq!(dist, 0.2); + +} + +#[test] + +fn test_hellinger() { + let length = 9; + let mut p_data = Vec::with_capacity(length); + let mut q_data = Vec::with_capacity(length); + for _ in 0..length { + p_data.push(1./length as f32); + q_data.push(1./length as f32); + } + p_data[0] -= 1./(2*length) as f32; + p_data[1] += 1./(2*length) as f32; + // + let dist = DistHellinger.eval(&p_data, &q_data); + + let dist_exact_fn = | n : usize | -> f32 { let d1 = (4. - (6 as f32).sqrt() - (2 as f32).sqrt())/n as f32 ; + d1.sqrt()/(2 as f32).sqrt() + }; + let dist_exact = dist_exact_fn(length); + // + log::info!("dist computed {:?} dist exact{:?} ", dist, dist_exact); + println!("dist computed {:?} , dist exact {:?} ", dist, dist_exact); + // + assert!((dist-dist_exact).abs() < 1.0e-5 ); + +} + + +#[test] + +fn test_jeffreys() { + + let length = 9; + let mut p_data : Vec = Vec::with_capacity(length); + let mut q_data : Vec = Vec::with_capacity(length); + for _ in 0..length { + p_data.push(1./length as f32); + q_data.push(1./length as f32); + } + p_data[0] -= 1./(2*length) as f32; + p_data[1] += 1./(2*length) as f32; + // + let dist_eval = DistJeffreys.eval(&p_data, &q_data); + let mut dist_test = 0.; + for i in 0..length { + dist_test += (p_data[i] - q_data[i]) * (p_data[i].max(M_MIN)/q_data[i].max(M_MIN)).ln(); + } + // + log::info!("dist eval {:?} dist test{:?} ", dist_eval, dist_test); + println!("dist eval {:?} , dist test {:?} ", dist_eval, dist_test); + assert!(dist_test >= 0.); + assert!((dist_eval-dist_test).abs() < 1.0e-5 ); +} + +} // end of module tests diff --git a/src/hnsw.rs b/src/hnsw.rs new file mode 100644 index 0000000..05bc52e --- /dev/null +++ b/src/hnsw.rs @@ -0,0 +1,1300 @@ +#![allow(dead_code)] + + +//! A rust implementation of Approximate NN search from: +//! Efficient and robust approximate nearest neighbour search using Hierarchical Navigable +//! small World graphs. +//! Yu. A. Malkov, D.A Yashunin 2016, 2018 + + +use std::cmp::Ordering; + +use parking_lot::{RwLock, Mutex}; +use std::sync::Arc; +use rayon::prelude::*; +use std::sync::mpsc::channel; + +use typename::TypeName; + +use hashbrown::HashMap; +use std::collections::binary_heap::BinaryHeap; + +use log::debug; +use log::trace; + +pub use crate::dist::Distance; + + +// TODO +// Profiling. + + +/// basic data in nodes + + + +/// maximum number of layers +pub(crate) const NB_LAYER_MAX:u8 = 16; // so max layer is 15!! + +#[derive(Debug, Default, Clone,Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] +/// The 2-uple represent layer as u8 and rank in layer as a i32 as stored in our structure +pub struct PointId(pub u8, pub i32); + + +/// this type is for an identificateur of each data vector, given by client. +/// can be the rank of data in an array, a hash value or anything that permits +/// retrieving the data +pub type DataId = usize; + + +pub type PointDistance = Box< dyn Distance >; + + +/// A structure containing internal pointId with distance to this pointId. +/// The order is given by ordering the distance to the point it refers to. +/// So points ordering has a meaning only has points refers to the same point +#[derive(Debug, Clone, Copy)] +pub struct PointIdWithOrder { + /// the identificateur of the point for which we store a distance + pub point_id: PointId, + /// The distance to a reference point (not represented in the structure) + pub dist_to_ref : f32, +} + + + +impl PartialEq for PointIdWithOrder { + fn eq(&self, other: &PointIdWithOrder) -> bool { + return self.dist_to_ref == other.dist_to_ref + } // end eq +} + + + +// order points by distance to self. +impl PartialOrd for PointIdWithOrder { + fn partial_cmp(&self, other: &PointIdWithOrder) -> Option { + self.dist_to_ref.partial_cmp(& other.dist_to_ref) + } // end cmp +} // end impl PartialOrd + + + +impl From< &PointWithOrder > for PointIdWithOrder { + fn from(point : &PointWithOrder ) -> PointIdWithOrder { + PointIdWithOrder::new(point.point_ref.p_id, point.dist_to_ref) + } +} + +impl PointIdWithOrder { + pub fn new(point_id: PointId, dist_to_ref:f32) -> Self { + PointIdWithOrder{point_id, + dist_to_ref, + } + } + +} // end of impl block + + +//======================================================================================= +/// The struct giving an answer point to a search request. +/// This structure is exported to other language API. +/// First field is origin id of the request point, second field is distance to request point +#[repr(C)] +#[derive(Debug, Copy, Clone, PartialEq, Default)] +pub struct Neighbour { + /// identification of data vector as given in initializing hnsw + pub d_id : DataId, + /// distance of neighbours + pub distance :f32, + /// point identification inside layers + pub p_id :PointId, +} + +impl Neighbour { + pub fn new(d_id : DataId, distance: f32, p_id: PointId) -> Neighbour { + Neighbour{d_id, distance, p_id} + } +} + +//======================================================================================= + +/// The basestructure representing a data point. +/// Its constains data as coming from the client, its client id, +/// and position in layer representation and neighbours. +/// +// neighbours table : one vector by layer so neighbours is allocated to NB_LAYER_MAX +// +#[derive(Debug, Clone)] +pub struct Point { + /// The data of this point, coming from hnsw client and associated to origin_id, + v: Vec, + /// an id coming from client using hnsw + origin_id : usize, + /// a point id identifying point as stored in our structure + p_id: PointId, + /// neighbours info + pub(crate) neighbours:Arc> >> >>, +} + +impl Point { + pub fn new(v: &Vec, origin_id: usize, p_id:PointId) -> Self { + let mut neighbours = Vec::with_capacity(NB_LAYER_MAX as usize); + // CAVEAT, perhaps pass nb layer as arg ? + for _ in 0..NB_LAYER_MAX { + neighbours.push(Vec::> >::new()); + } + Point{v:v.clone(), origin_id, p_id, neighbours: Arc::new(RwLock::new(neighbours))} + } + + + /// get a reference to vector data + pub fn get_v(&self) -> &[T] { + & self.v.as_slice() + } + + /// return coordinates in indexation + pub fn get_point_id(&self) -> PointId { + self.p_id + } + + /// returns external (or client id) id of point + pub fn get_origin_id(&self) -> usize { + self.origin_id + } + + // returns for each layer, a vector Neighbour of a point, one vector by layer + // + pub fn get_neighborhood_id(&self) -> Vec> { + let ref_neighbours = self.neighbours.read(); + let nb_layer = ref_neighbours.len(); + let mut neighborhood = Vec::>::with_capacity(nb_layer); + for i in 0..nb_layer { + let mut neighbours = Vec::::new(); + let nb_ngbh = ref_neighbours[i].len(); + if nb_ngbh > 0usize { + neighbours.reserve(nb_ngbh); + for pointwo in &ref_neighbours[i] { + neighbours.push(Neighbour::new(pointwo.point_ref.get_origin_id(),pointwo.dist_to_ref, + pointwo.point_ref.get_point_id())); + } + } + neighborhood.push(neighbours); + } + neighborhood + } + + + pub fn debug_dump(&self) { + println!(" \n dump of point id : {:?}", self.p_id); + println!("\n origin id : {:?} ", self.origin_id); + println!(" neighbours : ..."); + let ref_neighbours = self.neighbours.read(); + for i in 0..ref_neighbours.len() { + if ref_neighbours[i].len() > 0usize { + println!("neighbours at layer {:?}", i); + for n in &ref_neighbours[i] { + println!(" {:?}",n.point_ref.p_id); + } + } + } + println!(" neighbours dump : end"); + } + +} // end of block + +//=========================================================================================== + +/// A structure to store neighbours for of a point. +#[derive(Debug, Clone)] +pub(crate) struct PointWithOrder { + /// the identificateur of the point for which we store a distance to a point for which + /// we made a request. + point_ref: Arc>, + /// The distance to a point_ref to the request point (not represented in the structure) + dist_to_ref : f32, +} + + + +impl PartialEq for PointWithOrder { + fn eq(&self, other: &PointWithOrder) -> bool { + return self.dist_to_ref == other.dist_to_ref + } // end eq +} + + + +impl Eq for PointWithOrder {} + + +// order points by distance to self. +impl PartialOrd for PointWithOrder { + fn partial_cmp(&self, other: &PointWithOrder) -> Option { + self.dist_to_ref.partial_cmp(& other.dist_to_ref) + } // end cmp +} // end impl PartialOrd + + + +impl Ord for PointWithOrder { + fn cmp(&self, other: &PointWithOrder) -> Ordering { + if !self.dist_to_ref.is_nan() && !other.dist_to_ref.is_nan() { + self.dist_to_ref.partial_cmp(& other.dist_to_ref).unwrap() + } + else { + panic!("got a NaN in a distance"); + } + } // end cmp +} + + +impl PointWithOrder { + pub fn new(point_ref: &Arc>, dist_to_ref:f32) -> Self { + PointWithOrder{point_ref: Arc::clone(point_ref), + dist_to_ref, + } + } +} // end of impl block + + + + +//============================================================================================ + + + +// LayerGenerator +use rand::prelude::*; +use rand::distributions::{Uniform}; + +/// a struct to randomly generate a level for an item according to an exponential law +/// of parameter given by scale. +/// The distribution is constrained to be in [0..maxlevel[ +pub struct LayerGenerator { + rng: Arc>, + unif: Uniform, + scale: f64, + maxlevel: usize, +} + +impl LayerGenerator { + pub fn new(max_nb_connection: usize, maxlevel:usize) -> Self { + let scale = 1./(max_nb_connection as f64).ln(); + LayerGenerator{rng: Arc::new(Mutex::new(StdRng::from_entropy())), + unif: Uniform::::new(0.,1.), scale:scale, maxlevel:maxlevel} + } + // + // l=0 most densely packed layer + // if S is scale we sample so that P(l=n) = exp(-n/S) - exp(- (n+1)/S) + // with S = 1./ln(max_nb_connection) P(l >= maxlevel) = exp(-maxlevel * ln(max_nb_connection)) + // for nb_conn = 10, even with maxlevel = 10, we get P(l >= maxlevel) = 1.E-13 + // In Malkov(2016) S = 1./log(max_nb_connection) + // + /// generate a layer with given maxlevel. upper layers (higher index) are of decreasing probabilities. + /// thread safe method. + fn generate(&self) -> usize { + let mut protected_rng = self.rng.lock(); + let xsi = protected_rng.sample(self.unif); + let level = -xsi.ln() * self.scale; + let mut ulevel = level.floor() as usize; + // we redispatch possibly sampled level >= maxlevel to required range + if ulevel >= self.maxlevel { + // This occurs with very low probability. Cf commentary above. + ulevel = protected_rng.sample(Uniform::::new(0,self.maxlevel)); + } + ulevel + } + + /// just to try some variations on exponential level sampling. Unused. + pub fn set_scale_modification(&mut self, scale_modification : f64) { + self.scale = 1./( (1./self.scale) + scale_modification.ln()); + } + +} // end impl for LayerGenerator + + +// ==================================================================== + +/// a structure for indexation of points in layer +pub struct PointIndexation { + /// max number of connection for a point at a layer + pub(crate) max_nb_connection: usize, + /// + pub(crate) max_layer: usize, + /// needs at least one representation of points. points_by_layers[i] gives the points in layer i + pub(crate) points_by_layer: Arc> >> >>, + /// utility to generate a level + pub(crate) layer_g: LayerGenerator, + /// number of points in indexed structure + pub(crate) nb_point: Arc>, + /// curent enter_point: an Arc RwLock on a possible Arc Point + pub(crate) entry_point : Arc> >> >, +} + + +impl PointIndexation { + pub fn new(max_nb_connection: usize, max_layer:usize, max_elements:usize) -> Self { + let mut points_by_layer = Vec::with_capacity(max_layer); + for i in 0..max_layer { // recall that range are right extremeity excluded + // compute fraction of points going into layer i and do expected memory reservation + let frac = (-(i as f64)/(max_layer as f64)).exp() - (-((i+1) as f64)/(max_layer as f64)); + let expected_size = ((frac * max_elements as f64).round()) as usize; + points_by_layer.push(Vec::with_capacity(expected_size)); + } + let layer_g = LayerGenerator::new(max_nb_connection, max_layer); + PointIndexation{ + max_nb_connection, + max_layer, + points_by_layer: Arc::new(RwLock::new(points_by_layer)), + layer_g, + nb_point:Arc::new(RwLock::new(0)), + entry_point: Arc::new(RwLock::new(None)), + } + } // end of new + + /// + fn get_max_level_observed(&self) -> u8 { + let opt = self.entry_point.read(); + match opt.as_ref() { + Some(arc_point) => arc_point.p_id.0, + None => 0, + } + } + /// + fn debug_dump(&self) { + println!(" debug dump of PointIndexation"); + let max_level_observed = self.get_max_level_observed(); + // CAVEAT a lock once + for l in 0..=max_level_observed as usize { + println!(" layer {} : length : {} ", l, self.points_by_layer.read()[l].len()); + } + } + + + /// real insertion of point in point indexation + // generate a new Point/ArcPoint (with neigbourhood info empty) and store it in global table + // The function is called by Hnsw insert method + fn generate_new_point(& self, data : &Vec, origin_id : usize) -> (Arc> , usize) { + // get a write lock at the beginning of the function + let level = self.layer_g.generate(); + let new_point; + { // open a write lock on points_by_layer + let mut points_by_layer_ref = self.points_by_layer.write(); + let mut p_id = PointId{0:level as u8, 1:-1}; + p_id.1 = points_by_layer_ref[p_id.0 as usize].len() as i32; + // make a Point and then an Arc + let point = Point::new(data, origin_id, p_id); + new_point = Arc::new(point); + log::trace!("definitive pushing of point {:?}", p_id); + points_by_layer_ref[p_id.0 as usize].push(Arc::clone(&new_point)); + } // close write lock on points_by_layer + // + let nb_point; + { + let mut lock_nb_point = self.nb_point.write(); + *lock_nb_point += 1; + nb_point = *lock_nb_point; + if nb_point % 50000 == 0 { + println!(" setting number of points {:?} ", nb_point); + } + } + log::trace!(" setting number of points {:?} ", *self.nb_point); + // Now possibly this is a point on a new layer that will have no neighbours in its layer + return (Arc::clone(&new_point), nb_point); + } // end of insert + + + /// check if entry_point is modified + fn check_entry_point(&self, new_point : &Arc>) { + // + // take directly a write lock so that we are sure nobody can change anything between read and write + // of entry_point_id + log::trace!("trying to get a lock on entry point"); + let mut entry_point_ref = self.entry_point.write(); + match entry_point_ref.as_ref() { + Some(arc_point) => { + if new_point.p_id.0 > arc_point.p_id.0 { + log::debug!("Hnsw , inserting entry point {:?} ", new_point.p_id); + log::debug!("\n PointIndexation insert setting max level from {:?} to {:?}", arc_point.p_id.0 , new_point.p_id.0); + *entry_point_ref = Some(Arc::clone(new_point)); + } + } + None => { + log::trace!("initializing entry point"); + log::debug!("Hnsw , inserting entry point {:?} ", new_point.p_id); + *entry_point_ref = Some(Arc::clone(new_point)); + } + } + } // end of check_entry_point + + /// returns the number of points in layered structure + pub fn get_nb_point(&self) ->usize { + *self.nb_point.read() + } + + /// returns the size of data vector in graph if any, else return 0 + pub fn get_data_dimension(&self) -> usize { + let ep = self.entry_point.read(); + let dim = match ep.as_ref() { + Some(point) => point.get_v().len(), + None => 0, + }; + dim + } +} // end of impl PointIndexation + + + +//============================================================================================ + +/// an iterator on points stored. +/// The iteration begins at level 0 (most populated level) and goes upward in levels. +/// Must not be used during parallel insertion. +pub struct IterPoint<'a,T:Clone+Copy+Send+Sync> { + point_indexation : &'a PointIndexation, + layer:i64, + slot_in_layer:i64, +} + +impl <'a, T:Clone+Copy+Send+Sync> IterPoint<'a,T>{ + pub fn new(point_indexation : &'a PointIndexation) -> Self { + IterPoint{ point_indexation, layer:-1, slot_in_layer : -1 } + } +} // end of block impl IterPoint + + +/// iterator for layer 0 to upper layer. +impl <'a,T:Clone+Copy+Send+Sync> Iterator for IterPoint<'a,T> { + type Item = Arc>; + // + fn next(&mut self) -> Option{ + if self.layer == -1 { + self.layer = 0; + self.slot_in_layer = 0; + } + if (self.slot_in_layer as usize) < self.point_indexation.points_by_layer.read()[self.layer as usize].len() { + let slot = self.slot_in_layer as usize; + self.slot_in_layer += 1; + return Some(self.point_indexation.points_by_layer.read()[self.layer as usize][slot].clone()); + } + else { + self.slot_in_layer = 0; + self.layer += 1; + // must reach a non empty layer if possible + let entry_point_ref = self.point_indexation.entry_point.read(); + let points_by_layer = self.point_indexation.points_by_layer.read(); + let entry_point_level = entry_point_ref.as_ref().unwrap().p_id.0; + while (self.layer as u8) <= entry_point_level && + points_by_layer[self.layer as usize ].len() == 0 { + self.layer += 1; + } + // now here either (self.layer as u8) > self.point_indexation.max_level_observed + // or self.point_indexation.points_by_layer[self.layer as usize ].len() > 0 + if (self.layer as u8) <= entry_point_level { + let slot = self.slot_in_layer as usize; + self.slot_in_layer += 1; + return Some(points_by_layer[self.layer as usize][slot].clone()); + } + else { + return None; + } + } + } // end of next + +} // end of impl Iterator + + +impl<'a,T:Clone+Copy+Send+Sync> IntoIterator for &'a PointIndexation { + type Item = Arc>; + type IntoIter = IterPoint<'a,T>; + // + fn into_iter(self) -> Self::IntoIter { + IterPoint::new(self) + } + +} + + +// ============================================================================================ + +// The fields are made pub(crate) to be able to initialize struct from hnswio +/// The Base structure for hnsw implementation. +/// The main useful functions are : new, insert, insert_parallel, search and parallel_search. +/// Other functions are mainly for others crate to get access to some fields. +#[allow(dead_code)] +pub struct Hnsw+TypeName > { + /// asked number of candidates in search + pub(crate) ef_construction : usize, + /// maximum number of connection by layer for a point + pub(crate) max_nb_connection : usize, + /// flag to enforce that we have ef candidates as pruning strategy can discard some points + /// Can be set to true with method :set_extend_candidates + /// When set to true used only in base layer. + pub(crate) extend_candidates: bool, + /// defuault to false + pub(crate) keep_pruned: bool, + /// max layer , recall rust is in 0..maxlevel right bound excluded + pub(crate) max_layer: usize, + /// The global table containing points + pub(crate) layer_indexed_points: PointIndexation, + /// dimension data stored in points + pub(crate) data_dimension : usize, + /// distance between points. initialized at first insertion + pub(crate) dist_f : D, + /// insertion mode or searching mode. This flag prevents a internal thread to do a write when searching with other threads. + pub(crate) searching : bool, +} // end of Hnsw + + +impl +Send+Sync+TypeName > Hnsw { + /// allocation function + /// . max_nb_connection : number of neighbours stored in tables. + /// . ef_construction : controls numbers of neighbours explored during construction. See README or paper. + /// . max_elements : hint to speed up allocation tables. number of elements expected. + /// . f : the distance function + pub fn new(max_nb_connection : usize , max_elements:usize, max_layer:usize, ef_construction: usize, f:D) -> Self { + let adjusted_max_layer = (NB_LAYER_MAX as usize).min(max_layer); + let layer_indexed_points = PointIndexation::::new(max_nb_connection, adjusted_max_layer, max_elements); + let extend_candidates = false; + let keep_pruned = false; + // + log::info!("Hnsw max_nb_connection {:?}", max_nb_connection); + log::info!("Hnsw nb elements {:?}", max_elements); + log::info!("Hnsw ef_construction {:?}", ef_construction); + log::info!("Hnsw distance {:?}", f.type_name_of()); + log::info!("Hnsw extend candidates {:?}", extend_candidates); + // + Hnsw{max_nb_connection, + ef_construction:ef_construction, + extend_candidates:extend_candidates, + keep_pruned: keep_pruned, + max_layer: adjusted_max_layer, + layer_indexed_points: layer_indexed_points, + data_dimension : 0, + dist_f: f, + searching : false, + } + } // end of new + + /// get ef_construction used in graph creation + pub fn get_ef_construction(&self) -> usize { + self.ef_construction + } + /// + pub fn get_max_level(&self) -> usize { + self.max_layer + } + pub fn get_max_level_observed(&self) -> u8 { + self.layer_indexed_points.get_max_level_observed() as u8 + } + /// + pub fn get_max_nb_connection(&self) -> u8 { + self.max_nb_connection as u8 + } + /// returns number of points stored in hnsw structure + pub fn get_nb_point(&self) -> usize { + self.layer_indexed_points.get_nb_point() + } + /// set searching mode. + /// It is not possible to do parallel insertion and parallel searching simultaneously in different threads + /// so to enable searching after parallel insertion the flag must be set to true. + /// To resume parallel insertion reset the flag to false and so on. + pub fn set_searching_mode(&mut self, flag: bool) { + // must use an atomic! + self.searching = flag; + } + /// get name if distance + pub fn get_distance_name(&self) -> String { + D::type_name().to_string() + } + /// set the flag asking to keep pruned vectors by Navarro's heuristic (see Paper). + pub fn set_keeping_pruned(&mut self, flag: bool) { + self.keep_pruned = flag; + } + + /// set extend_candidates to given flag. By default it is false. + /// Only used in the level 0 layer during insertion (see the paper) + /// flag to enforce that we have ef candidates neighbours examined as pruning strategy + /// can discard some points + pub fn set_extend_candidates(&mut self, flag:bool) { + self.extend_candidates = flag; + } + + // multiplicative factor applied to default scale. Must between 0.5 and 1. + // more than 1. gives more occupied layers. This is just to experiment + // parameters variations on the algorithm but not general use. + fn set_scale_modification(&mut self, scale_modification : f64) { + println!("\n scale modification factor {:?}, scale value : {:?} (factor must be between 0.5 and 2.)", + scale_modification, self.layer_indexed_points.layer_g.scale); + // + if scale_modification >= 0.5 && scale_modification <= 2. { + self.layer_indexed_points.layer_g.set_scale_modification(scale_modification); + println!(" new scale value {:?}", self.layer_indexed_points.layer_g.scale); + } + else { + println!("\n scale modificationarg {:?} not valid , factor must be between 0.5 and 2.)", scale_modification); + } + } // end of set_scale_modification + + // here we could pass a point_id_with_order instead of entry_point_id: PointId + // The efficacity depends on greedy part depends on how near entry point is from point. + // ef is the number of points to return + // The method returns a BinaryHeap with positive distances. The caller must transforms it according its need + /// + /// Greedy algorithm n° 2 in Malkov paper. + /// search in a layer (layer) for the ef points nearest a point to be inserted in hnsw. + fn search_layer(& self, point: &Vec, entry_point: Arc> , ef:usize, layer: u8) -> BinaryHeap> > { + // + trace!("entering search_layer with entry_point_id {:?} layer : {:?} ef {:?} ", entry_point.p_id, layer, ef); + // + // here we allocate a skiplist on values not on reference beccause we want to return + // log2(skiplist_size) must be greater than 1. + let skiplist_size = ef.max(2); + // we will store positive distances in this one + let mut return_points = BinaryHeap::> >::with_capacity(skiplist_size); + // + if self.layer_indexed_points.points_by_layer.read()[layer as usize].len() == 0 { + // at the beginning we can have nothing in layer + trace!("search layer {:?}, empty layer", layer); + return return_points; + } + if entry_point.p_id.1 < 0 { + trace!("search layer negative point id : {:?}", entry_point.p_id); + return return_points; + } + // initialize visited points + let dist_to_entry_point = self.dist_f.eval(point, & entry_point.v); + log::trace!(" distance to entry point: {:?} ", dist_to_entry_point); + // keep a list of id visited + let mut visited_point_id = HashMap::>>::new(); + visited_point_id.insert(entry_point.p_id, Arc::clone(&entry_point)); + // + let mut candidate_points = BinaryHeap::> >::with_capacity(skiplist_size); + + candidate_points.push(Arc::new(PointWithOrder::new(&entry_point, -dist_to_entry_point))); + return_points.push(Arc::new(PointWithOrder::new(&entry_point, dist_to_entry_point))); + // at the beginning candidate_points contains point passed as arg in layer entry_point_id.0 + while candidate_points.len() > 0 { + // get nearest point in candidate_points + let c = candidate_points.pop().unwrap(); + // f farthest point to + let f = return_points.peek().unwrap(); + assert!(f.dist_to_ref >= 0.); + assert!(c.dist_to_ref <= 0.); + log::trace!("comparaing c : {:?} f : {:?}", -(c.dist_to_ref), f.dist_to_ref); + if -(c.dist_to_ref) > f.dist_to_ref { + // this comparison requires that we are sure that distances compared are distances to the same point : + // This is the case we compare distance to point passed as arg. + log::trace!("fast return from search_layer, nb points : {:?} \n \t c {:?} \n \t f {:?} dists: {:?} {:?}", + return_points.len(), c.point_ref.p_id, f.point_ref.p_id, -(c.dist_to_ref), f.dist_to_ref); + return return_points; + } + // now we scan neighborhood of c in layer and increment visited_point, candidate_points + // and optimize candidate_points so that it contains points with lowest distances to point arg + // + let neighbours_c_l = &c.point_ref.neighbours.read()[layer as usize]; + let c_pid = c.point_ref.p_id; + log::trace!(" search_layer, {:?} has nb neighbours : {:?} ", c_pid, neighbours_c_l.len()); + for e in neighbours_c_l { + // HERE WE sEE THAT neighbours should be stored as PointIdWithOrder !! + // CAVEAT what if several point_id with same distance to ref point? + if visited_point_id.contains_key(&e.point_ref.p_id) != true { + visited_point_id.insert(e.point_ref.p_id, Arc::clone(&e.point_ref)); + log::trace!(" visited insertion {:?}", e.point_ref.p_id); + let f_opt = return_points.peek(); + if f_opt.is_none() { + // do some debug info, dumped distance is from e to c! as e is in c neighbours + debug!("return points empty when inserting {:?}", e.point_ref.p_id); + return return_points; + } + let f = f_opt.unwrap(); + let e_dist_to_p = self.dist_f.eval(point, & e.point_ref.v); + let f_dist_to_p = f.dist_to_ref; + if e_dist_to_p < f_dist_to_p || return_points.len() < ef { + let e_prime = Arc::new(PointWithOrder::new(&e.point_ref, e_dist_to_p)); + // a neighbour of neighbour is better, we insert it into candidate with the distance to point + log::trace!(" inserting new candidate {:?}", e_prime.point_ref.p_id); + candidate_points.push(Arc::new(PointWithOrder::new(&e.point_ref, -e_dist_to_p))); + return_points.push(Arc::clone(&e_prime)); + if return_points.len() > ef { + return_points.pop(); + } + } // end if e.dist_to_ref < f.dist_to_ref + } + } // end of for on neighbours_c + } // end of while in candidates + // + trace!("return from search_layer, nb points : {:?}", return_points.len()); + return_points + } // end of search_layer + + + // Hnsw insert. + /// insert a data vector with its external id as given by the client. + /// The insertion method will give the point an internal id. + pub fn insert(&self, data_with_id: (&Vec,usize)) { + // + let (data , origin_id) = data_with_id; + let keep_pruned = self.keep_pruned; + // insert in indexation and get point_id adn generate a new entry_point if necessary + let (new_point, point_rank) = self.layer_indexed_points.generate_new_point(data, origin_id); + log::trace!("\n\n Hnsw insert generated new point {:?} ", new_point.p_id); + // now real work begins + // allocate a binary heap + let level = new_point.p_id.0; + let mut enter_point_copy = None; + let mut max_level_observed = 0; + // entry point has been set in + { // I open a read lock on an option + if let Some(arc_point) = self.layer_indexed_points.entry_point.read().as_ref() { + enter_point_copy = Some(Arc::clone(&arc_point)); + if point_rank == 1 { + log::debug!("Hnsw stored first point , direct return {:?} ", new_point.p_id); + return; + } + max_level_observed = enter_point_copy.as_ref().unwrap().p_id.0; + } + } + if enter_point_copy.is_none() { + self.layer_indexed_points.check_entry_point(&new_point); + return; + } + let mut dist_to_entry = self.dist_f.eval(data , & enter_point_copy.as_ref().unwrap().v); + // we go from self.max_level_observed to level+1 included + for l in ((level+1)..(max_level_observed+1)).rev() { + // CAVEAT could bypass when layer empty, avoid allocation.. + let mut sorted_points = self.search_layer(&data, Arc::clone(enter_point_copy.as_ref().unwrap()), 1, l); + log::trace!("in insert :search_layer layer {:?}, returned {:?} points ", l, sorted_points.len()); + if sorted_points.len() > 1 { + panic!("in insert : search_layer layer {:?}, returned {:?} points ", l, sorted_points.len()); + } + // the heap conversion is useless beccause of the preceding test. + // sorted_points = from_positive_binaryheap_to_negative_binary_heap(&mut sorted_points); + // + if let Some(ep) = sorted_points.pop() { + // get the lowest distance point + let tmp_dist = self.dist_f.eval(data , & ep.point_ref.v); + if tmp_dist < dist_to_entry { + enter_point_copy = Some(Arc::clone(& ep.point_ref)); + dist_to_entry = tmp_dist; + } + } + else { + // this layer is not yet filled + log::trace!("layer still empty {} : got null list", l); + } + } + // now enter_point_id_copy contains id of nearest + // now loop down to 0 + for l in (0..level+1).rev() { + let ef = self.ef_construction; + // when l == level, we cannot get new_point in sorted_points as it is seen only from declared neighbours + let mut sorted_points = self.search_layer(&data, Arc::clone(enter_point_copy.as_ref().unwrap()), ef, l); + log::trace!("in insert :search_layer layer {:?}, returned {:?} points ", l, sorted_points.len()); + sorted_points = from_positive_binaryheap_to_negative_binary_heap(&mut sorted_points); + if sorted_points.len() > 0 { + let nb_conn; + let extend_c; + if l==0 { + nb_conn = 2*self.max_nb_connection; + extend_c = self.extend_candidates; + } + else { + nb_conn = self.max_nb_connection; + extend_c = false; + } + let mut neighbours = Vec::> >::with_capacity(nb_conn); + self.select_neighbours(&data, &mut sorted_points, nb_conn, extend_c, l, keep_pruned, &mut neighbours); + // sort neighbours + neighbours.sort_unstable(); + // we must add bidirecti*onal from data i.e new_point_id to neighbours + new_point.neighbours.write()[l as usize] = neighbours.clone(); + // this reverse neighbour update could be done here but we put it at end to gather all code + // requiring a mutex guard for multi threading. + // update ep for loop iteration. As we sorted neighbours the nearest + if neighbours.len() > 0 { + enter_point_copy = Some(Arc::clone(&neighbours[0].point_ref)); + } + } + } // for l + // + // new_point has been inserted at the beginning in table + // so that we can call reverse_update_neighborhoodwe consitently + // now reverse update of neighbours. + self.reverse_update_neighborhood_simple (Arc::clone(&new_point)); + // + self.layer_indexed_points.check_entry_point(&new_point); + // + log::trace!("Hnsw exiting insert new point {:?} ", new_point.p_id); + } // end of insert + + + /// insert in parallel a slice of Vec. Uses Rayon. + /// The number of insertions asked for must be large enough to be efficient. + /// Typically 1000 * the number of threads. + pub fn parallel_insert(&self, datas: &Vec<(&Vec, usize)> ) { + datas.par_iter().for_each( |&item| self.insert(item)); + } // end of parallel_insert + + + + + /// insert new_point in neighbourhood info of point + fn reverse_update_neighborhood_simple(&self, new_point: Arc>) { + // println!("reverse update neighbourhood for new point {:?} ", new_point.p_id); + log::trace!("reverse update neighbourhood for new point {:?} ", new_point.p_id); + let level = new_point.p_id.0; + for l in (0..level+1).rev() { + for q in &new_point.neighbours.read()[l as usize] { + if new_point.p_id != q.point_ref.p_id { + // as new point is in global table, do not loop and deadlock!! + let q_point = &q.point_ref; + let mut q_point_neighbours = q_point.neighbours.write(); + let n_to_add = PointWithOrder::::new(&Arc::clone(&new_point), q.dist_to_ref); + // must be sure that we add a point at the correct level + let l_n = n_to_add.point_ref.p_id.0 as usize; + let already = q_point_neighbours[l_n as usize].iter().position(|old| old.point_ref.p_id == new_point.p_id); + if already.is_some() { + // log::debug!(" new_point.p_id {:?} already in neighbourhood of q_point {:?} at index {:?}", new_point.p_id, q_point.p_id, already.unwrap()); + // q_point.debug_dump(); cannot be called as its neighbours are locked write by this method. + // new_point.debug_dump(); + // panic!(); + continue; + } + q_point_neighbours[l_n].push(Arc::new(n_to_add)); + let nbn_at_l = q_point_neighbours[l_n].len(); + + // + // if l < level, update upward chaining, insert does a sort! t_q has a neighbour not yet in global table of points! + let threshold_shrinking; // TO DO optimize threshold + if l_n > 0 { + threshold_shrinking = self.max_nb_connection; + } + else { + threshold_shrinking = 2 * self.max_nb_connection; + } + let shrink = if nbn_at_l > threshold_shrinking { + true + } else { + false + }; + { // sort and shring if necessary + q_point_neighbours[l_n as usize].sort_unstable(); + if shrink { + q_point_neighbours[l_n as usize].pop(); + } + } + } // end protection against point identity + } + } + // println!(" exitingreverse update neighbourhood for new point {:?} ", new_point.p_id); + } // end of reverse_update_neighborhood_simple + + + + + + pub fn get_point_indexation(&self) -> &PointIndexation { + &self.layer_indexed_points + } + // This is best explained in : Navarro. Searching in metric spaces by spatial approximation. + /// simplest searh neighbours + // The binary heaps here is with negative distance sorted. + fn select_neighbours(&self, data: &[T], candidates: &mut BinaryHeap> >, + nb_neighbours_asked: usize, + extend_candidates_asked: bool, + layer:u8, + keep_pruned: bool, + neighbours_vec : &mut Vec> >) + { + // + log::trace!("entering select_neighbours : nb candidates: {}", candidates.len()); + // + neighbours_vec.clear(); + // we will extend if we do not have enough candidates and it is explicitly asked in arg + let mut extend_candidates = false; + if candidates.len() <= nb_neighbours_asked { + if !extend_candidates_asked { + // just transfer taking care of signs + while !candidates.is_empty() { + let p = candidates.pop().unwrap(); + assert!(-p.dist_to_ref >= 0.); + neighbours_vec.push(Arc::new(PointWithOrder::new(&p.point_ref, -p.dist_to_ref))); + } + return; + } + else { + extend_candidates = true; + } + } + // + // + //extend_candidates = true; + // + if extend_candidates { + let mut candidates_set = HashMap::> >::new(); + for c in candidates.iter() { + candidates_set.insert(c.point_ref.p_id, Arc::clone(&c.point_ref)); + } + let mut new_candidates_set = HashMap::> >::new(); + // get a list of all neighbours of candidates + for (_p_id, p_point) in candidates_set.iter() { + let n_p_layer = &p_point.neighbours.read()[layer as usize]; + for q in n_p_layer { + if !candidates_set.contains_key(& q.point_ref.p_id) && !new_candidates_set.contains_key(& q.point_ref.p_id) { + new_candidates_set.insert(q.point_ref.p_id, Arc::clone(&q.point_ref)); + } + } + } // end of for p + log::trace!("select neighbours extend candidates from : {:?} adding : {:?}", candidates.len(), new_candidates_set.len()); + for (_p_id, p_point) in new_candidates_set.iter() { + let dist_topoint = self.dist_f.eval(data, &p_point.v); + candidates.push(Arc::new(PointWithOrder::new(p_point, -dist_topoint))); + } + } // end if extend_candidates + // + let mut discarded_points = BinaryHeap:: > >::new(); + while candidates.len() > 0 && neighbours_vec.len() < nb_neighbours_asked { + // compare distances of e to data. we do not need to recompute dists! + if let Some(e_p) = candidates.pop() { + let mut e_to_insert = true; + let e_point_v = &e_p.point_ref.v; + assert!(e_p.dist_to_ref <= 0.); + // is e_p the nearest to reference? data than to previous neighbours + if neighbours_vec.len() > 0 { + e_to_insert = !neighbours_vec.iter().any(|d| + self.dist_f.eval(e_point_v, &(d.point_ref.v)) + <= -e_p.dist_to_ref); + } + if e_to_insert { + log::trace!("inserting neighbours : {:?} ", e_p.point_ref.p_id); + neighbours_vec.push(Arc::new(PointWithOrder::new(&e_p.point_ref, -e_p.dist_to_ref))); + } else { + log::trace!("discarded neighbours : {:?} ", e_p.point_ref.p_id); + // ep is taken from a binary heap, so it has a negative sign, we keep its sign + // to store it in another binary heap will possibly need to retain the best ones from the discarde binaryHeap + if keep_pruned { + discarded_points.push(Arc::new(PointWithOrder::new(&e_p.point_ref, e_p.dist_to_ref))); + } + } + } + } + // now this part of neighbours is the most interesting and is distance sorted. + + // not pruned are at the end of neighbours_vec which is not re-sorted , but discarded are sorted. + if keep_pruned { + while discarded_points.len() > 0 && neighbours_vec.len() < nb_neighbours_asked { + let best_point = discarded_points.pop().unwrap(); + // do not forget to reverse sign + assert!(best_point.dist_to_ref <= 0.); + neighbours_vec.push(Arc::new(PointWithOrder::new(&best_point.point_ref, -best_point.dist_to_ref))); + } + }; + // + if log::log_enabled!(log::Level::Trace) { + log::trace!("exiting select_neighbours : nb candidates: {}", neighbours_vec.len()); + for n in neighbours_vec { + println!(" neighbours {:?} ", n.point_ref.p_id); + } + } + // + } // end of select_neighbours + + + pub fn dump_layer_info(&self) { + self.layer_indexed_points.debug_dump(); + } + + + // search the first knbn nearest neigbours of a data, but can modify ef for layer > 1 + // This function return Vec >> + // The parameter ef controls the width of the search in the lowest level, it must be greater + // than number of neighbours asked. A rule of thumb could be between knbn and max_nb_connection. + fn search_general(&self, data :&Vec , knbn:usize, ef_arg:usize) -> Vec { + // + let mut entry_point; + { // a lock on an option an a Arc + let entry_point_opt_ref = self.layer_indexed_points.entry_point.read(); + if entry_point_opt_ref.is_none() { + return Vec::::new(); + } + else { + entry_point = Arc::clone((*entry_point_opt_ref).as_ref().unwrap()); + } + } + // + let mut dist_to_entry = self.dist_f.eval(data , & entry_point.as_ref().v); + for layer in (1..=entry_point.p_id.0).rev() { + let mut neighbours = self.search_layer(data, Arc::clone(&entry_point), 1, layer); + neighbours = from_positive_binaryheap_to_negative_binary_heap(&mut neighbours); + if let Some(entry_point_tmp) = neighbours.pop() { + // get the lowest distance point. + let tmp_dist = self.dist_f.eval(data , & entry_point_tmp.point_ref.v); + if tmp_dist < dist_to_entry { + entry_point = Arc::clone(&entry_point_tmp.point_ref); + dist_to_entry = tmp_dist; + } + } + } + // ef must be greater than knbn. Possibly it should be between knbn and self.max_nb_connection + let ef = ef_arg.max(knbn); + // now search with asked ef in layer 0 + let neighbours_heap = self.search_layer(data, entry_point, ef, 0); + // go from heap of points with negative dist to a sorted vec of increasing points with > 0 distances. + let neighbours = neighbours_heap.into_sorted_vec(); + // get the min of K and ef points into a vector. + // + let last = knbn.min(ef).min(neighbours.len()); + let knn_neighbours : Vec = + neighbours[0..last].iter().map(|p| Neighbour::new(p.as_ref().point_ref.origin_id, p.as_ref().dist_to_ref, p.as_ref().point_ref.p_id)).collect(); + + knn_neighbours + } // end of knn_search + + + /// search the first knbn nearest neigbours of a data and returns a Vector of Neighbour. + /// The parameter ef controls the width of the search in the lowest level, it must be greater + /// than number of neighbours asked. + /// A rule of thumb could be between knbn and max_nb_connection. + pub fn search(&self, data :&Vec , knbn:usize, ef_arg:usize) -> Vec { + // + let entry_point; + { // a lock on an option an a Arc + let entry_point_opt_ref = self.layer_indexed_points.entry_point.read(); + if entry_point_opt_ref.is_none() { + return Vec::::new(); + } + else { + entry_point = Arc::clone((*entry_point_opt_ref).as_ref().unwrap()); + } + } + // + let mut dist_to_entry = self.dist_f.eval(data , & entry_point.as_ref().v); + let mut pivot = Arc::clone(&entry_point); + let mut new_pivot = None; + // + for layer in (1..=entry_point.p_id.0).rev() { + let mut has_changed = false; + // search in stored neighbours + { + let neighbours = &pivot.neighbours.read()[layer as usize]; + for n in neighbours { + // get the lowest distance point. + let tmp_dist = self.dist_f.eval(data , & n.point_ref.v); + if tmp_dist < dist_to_entry { + new_pivot = Some(Arc::clone(&n.point_ref)); + has_changed = true; + dist_to_entry = tmp_dist; + } + } // end of for on neighbours + } + if has_changed { + pivot = Arc::clone(new_pivot.as_ref().unwrap()); + } + } // end on for on layers + // ef must be greater than knbn. Possibly it should be between knbn and self.max_nb_connection + let ef = ef_arg.max(knbn); + // now search with asked ef in layer 0 + let neighbours_heap = self.search_layer(data, pivot, ef, 0); + // go from heap of points with negative dist to a sorted vec of increasing points with > 0 distances. + let neighbours = neighbours_heap.into_sorted_vec(); + // get the min of K and ef points into a vector. + // + let last = knbn.min(ef).min(neighbours.len()); + let knn_neighbours : Vec = + neighbours[0..last].iter().map(|p| Neighbour::new(p.as_ref().point_ref.origin_id, p.as_ref().dist_to_ref, p.as_ref().point_ref.p_id)).collect(); + + knn_neighbours + } // end of knn_search + + + + fn search_with_id(&self, request : (usize, &Vec) , knbn:usize, ef:usize) -> (usize, Vec) { + (request.0, self.search(request.1, knbn, ef)) + } + + /// knbn is the number of nearest neigbours asked for. Returns for each data vector + /// a Vector of Neighbour + pub fn parallel_search(&self, datas: &Vec>, knbn:usize, ef:usize) -> Vec> { + let (sender, receiver) = channel(); + // make up requests + let nb_request = datas.len(); + let requests : Vec<(usize, &Vec)> = (0..nb_request).into_iter().zip(datas.iter()).collect(); + // + requests.par_iter().for_each_with(sender, |s, item| s.send(self.search_with_id(*item, knbn, ef)).unwrap()); + let req_res : Vec::<(usize, Vec) > = receiver.iter().collect(); + // now sort to respect the key order of input + let mut answers = Vec::>::with_capacity(datas.len()); + // get a map from request id to rank + let mut req_hash = HashMap::::new(); + for i in 0..req_res.len() { + // the response of request req_res[i].0 is at rank i + req_hash.insert(req_res[i].0, i); + } + for i in 0..datas.len() { + let answer_i = req_hash.get_key_value(&i).unwrap().1; + answers.push((req_res[*answer_i].1).clone()); + } + answers + } // end of insert_parallel + +} // end of Hnsw + + + +// This function takes a binary heap with points declared with a negative distance +// and returns a vector of points with their correct positive distance to some reference distance +// The vector is sorted by construction +fn from_negative_binaryheap_to_sorted_vector(heap_points : &mut BinaryHeap >> ) -> Vec > > { + let nb_points = heap_points.len(); + let mut vec_points = Vec:: >>::with_capacity(nb_points); + // + for p in heap_points.iter() { + assert!(p.dist_to_ref <= 0.); + let reverse_p = Arc::new(PointWithOrder::new(&p.point_ref, -p.dist_to_ref)); + vec_points.push(reverse_p); + } + log::trace!("from_negative_binaryheap_to_sorted_vector nb points in out {:?} {:?} ", nb_points, vec_points.len()); + vec_points +} + + + +// This function takes a binary heap with points declared with a positive distance +// and returns a binary_heap of points with their correct negative distance to some reference distance +// +fn from_positive_binaryheap_to_negative_binary_heap(positive_heap : &mut BinaryHeap >> ) -> BinaryHeap > > { + let nb_points = positive_heap.len(); + let mut negative_heap = BinaryHeap:: >>::with_capacity(nb_points); + // + for p in positive_heap.iter() { + assert!(p.dist_to_ref >= 0.); + let reverse_p = Arc::new(PointWithOrder::new(&p.point_ref, -p.dist_to_ref)); + negative_heap.push(reverse_p); + } + log::trace!("from_positive_binaryheap_to_negative_binary_heap nb points in out {:?} {:?} ", nb_points, negative_heap.len()); + negative_heap +} + + + +// essentialy to check dump/reload conssistency + +pub(crate) fn check_equality(hnsw1:&Hnsw, hnsw2: &Hnsw) + where T:Copy+Clone+Send+Sync+TypeName, D:Distance+TypeName+Default+Send+Sync { + // + log::debug!("\n in check_equality"); + // + assert_eq!(hnsw1.get_nb_point(), hnsw2.get_nb_point()); + // check for entry point + assert!(hnsw1.layer_indexed_points.entry_point.read().is_some() || + hnsw1.layer_indexed_points.entry_point.read().is_some(), "one entry point is None"); + let ep1_read = hnsw1.layer_indexed_points.entry_point.read(); + let ep2_read = hnsw2.layer_indexed_points.entry_point.read(); + let ep1 = ep1_read.as_ref().unwrap(); + let ep2 = ep2_read.as_ref().unwrap(); + assert_eq!(ep1.origin_id, ep2.origin_id, "different entry points {:?} {:?}", ep1.origin_id, ep2.origin_id); + assert_eq!(ep1.p_id, ep2.p_id, "origin id {:?} ", ep1.origin_id); + assert_eq!(ep1.v.len(), ep2.v.len(), "different length {:?} {:?}", ep1.v.len(), ep2.v.len()); + // check layers + let layers_1 = hnsw1.layer_indexed_points.points_by_layer.read(); + let layers_2 = hnsw2.layer_indexed_points.points_by_layer.read(); + let mut nb_point_checked = 0; + let mut nb_neighbours_checked = 0; + for i in 0..NB_LAYER_MAX as usize { + log::debug!("\n checking layer {:?}", i); + assert_eq!(layers_1[i].len() , layers_2[i].len()); + for j in 0..layers_1[i].len() { + let p1 = &layers_1[i][j]; + let p2 = &layers_2[i][j]; + assert_eq!(p1.origin_id, p2.origin_id); + assert_eq!(p1.p_id, p2.p_id, "\n checking origin_id point {:?} ", p1.origin_id); + nb_point_checked += 1; + // check neighborhood + let nbgh1 = p1.neighbours.read(); + let nbgh2 = p2.neighbours.read(); + assert_eq!(nbgh1.len(), nbgh2.len()); + for k in 0..nbgh1.len() { + assert_eq!(nbgh1[k].len(), nbgh2[k].len()); + for l in 0..nbgh1[k].len() { + assert_eq!(nbgh1[k][l].point_ref.origin_id, nbgh2[k][l].point_ref.origin_id); + assert_eq!(nbgh1[k][l].point_ref.p_id, nbgh2[k][l].point_ref.p_id); + // CAVEAT for precision with f32 + assert_eq!(nbgh1[k][l].dist_to_ref, nbgh2[k][l].dist_to_ref); + nb_neighbours_checked += 1; + } + } + } // end of for j + } // end of for i + assert_eq!(nb_point_checked, hnsw1.get_nb_point()); + log::debug!("nb neighbours checked {:?}", nb_neighbours_checked); + log::debug!("exiting check_equality"); + +} // end of check_reload + + + + +#[cfg(test)] + + +mod tests { + +use super::*; +use crate::dist; + +use rand::distributions::{Uniform}; + +use cpu_time::ProcessTime; + +#[test] + +fn test_iter_point() { + let _res = simple_logger::init(); + // + println!("\n\n test_iter_point"); + // + let mut rng = rand::thread_rng(); + let unif = Uniform::::new(0.,1.); + let nbcolumn = 5000; + let nbrow = 10; + let mut xsi; + let mut data = Vec::with_capacity(nbcolumn); + for j in 0..nbcolumn { + data.push(Vec::with_capacity(nbrow)); + for _ in 0..nbrow { + xsi = rng.sample(unif); + data[j].push(xsi); + } + } + // + // check insertion + let ef_construct= 25; + let nb_connection = 10; + let start = ProcessTime::now(); + let hns = Hnsw::::new(nb_connection, nbcolumn, 16, ef_construct, dist::DistL1{}); + for i in 0..data.len() { + hns.insert((&data[i], i)); + } + let cpu_time = start.elapsed(); + println!(" test_insert_iter_point time inserting {:?}", cpu_time); + + hns.dump_layer_info(); + // now check iteration + let mut ptiter = hns.get_point_indexation().into_iter(); + let mut nb_dumped = 0; + loop { + if let Some(_point) = ptiter.next() { + // println!("point : {:?}", _point.p_id); + nb_dumped += 1; + } + else { + break; + } + } // end while + // + assert_eq!(nb_dumped, nbcolumn); +} // end of test_insert_iter_point + +} // end of module test \ No newline at end of file diff --git a/src/hnswio.rs b/src/hnswio.rs new file mode 100644 index 0000000..ce9ed8e --- /dev/null +++ b/src/hnswio.rs @@ -0,0 +1,689 @@ + +// +//! This file provides io dump/ reload of computed graph. +/// +/// The main structure to dump (or reload) is PointIndexation. +/// +/// a dump is constituted of 2 files. One dumps jut the graph (or topology) +/// with id of points and another file dumps the ids and vector in point. +/// The graph file is suffixed hnsw.graph the other is suffixed hnsw.data +/// +/// +// datafile +// MAGICDATAP : u32 +// dimension : u32 +// The for each point the triplet: (MAGICDATAP, origin_id , array of values.) ( u32, u64, ....) +// +// A point is dumped in graph file as given by its external id (type DataId i.e : a usize, possibly a hash value) +// and layer (u8) and rank_in_layer:i32. +// In the data file the point dump consist in the triplet: (MAGICDATAP, origin_id , array of values.) +// + + +use std::io; +use std::mem; + +use parking_lot::{RwLock}; +use std::sync::Arc; +use std::collections::HashMap; +use std::path::{PathBuf}; + +use typename::TypeName; + +use std::io::prelude::*; +use crate::hnsw; +use self::hnsw::*; +use crate::dist::Distance; + +// magic before each graph point data for each point +const MAGICPOINT : u32 = 0x000a678f; +// magic at beginning of description +const MAGICDESCR : u32 = 0x000a677f; +// magic at beginning of a layer dump +const MAGICLAYER : u32 = 0x000a676f; +// magic head of data file and before each data vector +const MAGICDATAP : u32 = 0xa67f0000; + +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum DumpMode { + Light, + Full, +} + + +/// The main interface for dumping struct Hnsw. +pub trait HnswIO { + fn dump(&self, mode : DumpMode, outgraph : &mut io::BufWriter, outdata: &mut io::BufWriter) -> Result; +} + + + + +/// structure describing main parameters for hnsnw data and written at the beginning of a dump file. +/// +/// Name of distance and type of data must be encoded in the dump file for a coherent reload. +#[repr(C)] +pub struct Description { + /// value is 1 for Full 0 for Light + pub dumpmode : u8, + /// max number of connections in layers != 0 + pub max_nb_connection : u8, + /// number of observed layers + pub nb_layer : u8, + /// search parameter + pub ef: usize, + /// total number of points + pub nb_point: usize, + /// data dimension + pub dimension : usize, + /// name of distance + pub distname : String, + /// T typename + pub t_name : String, +} + +impl Description { + /// The dump of Description consists in : + /// . The value MAGICDESCR as a u32 (4 u8) + /// . The type of dump as u8 + /// . max_nb_connection as u8 + /// . ef (search parameter used in construction) as usize + /// . nb_point (the number points dumped) as a usize + /// . the name of distance used. (nb byes as a usize then list of bytes) + /// + fn dump(&self, argmode : DumpMode, out : &mut io::BufWriter) -> Result { + log::info!("in dump of description"); + out.write(unsafe { &mem::transmute::()]>(MAGICDESCR) } ).unwrap(); + let mode : u8 = match argmode { + DumpMode::Full => 1, + _ => 0, + }; + // CAVEAT should check mode == self.mode + out.write(unsafe { &mem::transmute::(mode) } ).unwrap(); + out.write(unsafe { &mem::transmute::(self.max_nb_connection) } ).unwrap(); + out.write(unsafe { &mem::transmute::(self.nb_layer) } ).unwrap(); + if self.nb_layer != NB_LAYER_MAX { + println!("dump of Description, nb_layer != NB_MAX_LAYER"); + return Err(String::from("dump of Description, nb_layer != NB_MAX_LAYER")); + } + out.write(unsafe { &mem::transmute::()]>(self.ef) } ).unwrap(); + log::info!("dumping nb point {:?}", self.nb_point); + // + out.write(unsafe { &mem::transmute::()]>(self.nb_point) } ).unwrap(); + log::info!("dumping dimension of data {:?}", self.dimension); + out.write(unsafe { &mem::transmute::()]>(self.dimension) } ).unwrap(); + // dump of distance name + let namelen : usize = self.distname.len(); + log::info!("distance name {:?} ", self.distname); + out.write(unsafe { &mem::transmute::()]>(namelen) } ).unwrap(); + out.write(self.distname.as_bytes()).unwrap(); + // dump of T value typename + let namelen : usize = self.t_name.len(); + log::info!("T name {:?} ", self.t_name); + out.write(unsafe { &mem::transmute::()]>(namelen) } ).unwrap(); + out.write(self.t_name.as_bytes()).unwrap(); + // + return Ok(1); + } // end fo dump + +} // end of HnswIO impl for Descr + + +/// This method is a preliminary to do a full reload from a dump. +/// The method load_hnsw needs to know the typename , distance used, and construction parameters. +/// So the reload is made in two steps. +pub fn load_description(io_in: &mut dyn Read) -> io::Result { + // + let mut descr = Description{ dumpmode: 0, max_nb_connection: 0, nb_layer: 0, + ef: 0, nb_point: 0, dimension : 0, distname: String::from(""), t_name : String::from("")}; + let magic : u32 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&magic as *const u32) as *mut u8, ::std::mem::size_of::() )}; + io_in.read_exact(it_slice)?; + log::debug!(" magic {:X} ", magic); + if magic != MAGICDESCR { + return Err(io::Error::new(io::ErrorKind::Other, "bad magic at descr beginning")); + } + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&descr.dumpmode as *const u8) as *mut u8, ::std::mem::size_of::() )}; + io_in.read_exact(it_slice)?; + log::info!(" dumpmode {:?} ", descr.dumpmode); + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&descr.max_nb_connection as *const u8) as *mut u8, ::std::mem::size_of::() )}; + io_in.read_exact(it_slice)?; + log::info!(" max_nb_connection {:?} ", descr.max_nb_connection); + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&descr.nb_layer as *const u8) as *mut u8, ::std::mem::size_of::() )}; + io_in.read_exact(it_slice)?; + log::info!("nb_layer {:?} ", descr.nb_layer); + // ef + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&descr.ef as *const usize) as *mut u8, ::std::mem::size_of::() )}; + io_in.read_exact(it_slice)?; + log::info!("ef {:?} ", descr.ef); + // nb_point + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&descr.nb_point as *const usize) as *mut u8, ::std::mem::size_of::() )}; + io_in.read_exact(it_slice)?; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&descr.dimension as *const usize) as *mut u8, ::std::mem::size_of::() )}; + io_in.read_exact(it_slice)?; + log::info!("nb_point {:?} dimension {:?}", descr.nb_point, descr.dimension); + // distance name + let len : usize = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&len as *const usize) as *mut u8, ::std::mem::size_of::() )}; + io_in.read_exact(it_slice)?; + log::debug!("length of distance name {:?} ", len); + if len > 256 { + println!(" length of distance name should not exceed 256"); + return Err(io::Error::new(io::ErrorKind::Other, "bad lenght for distance name")); + } + let mut distv = Vec::::new(); + distv.resize(len , 0); + io_in.read_exact(distv.as_mut_slice())?; + let distname = String::from_utf8(distv).unwrap(); + log::debug!("distance name {:?} ", distname); + descr.distname = distname; + // reload of type name + let len : usize = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&len as *const usize) as *mut u8, ::std::mem::size_of::() )}; + io_in.read_exact(it_slice)?; + log::debug!("length of T name {:?} ", len); + if len > 256 { + println!(" length of T name should not exceed 256"); + return Err(io::Error::new(io::ErrorKind::Other, "bad lenght for T name")); + } + let mut tnamev = Vec::::new(); + tnamev.resize(len , 0); + io_in.read_exact(tnamev.as_mut_slice())?; + let t_name = String::from_utf8(tnamev).unwrap(); + log::debug!("T type name {:?} ", t_name); + descr.t_name = t_name; + log::debug!(" end of description load \n"); + // + Ok(descr) +} +// +// dump and load of Point +// ========================== +// + + /// Graph part of point dump + /// dump of a point consist in + /// 1. The value MAGICPOINT + /// 2. its identity ( a usize rank in original data , hash value or else , and PointId) + /// 3. for each layer dump of the number of neighbours followed by : + /// for each neighbour dump of its identity (: usize) and then distance (): u32) to point dumped. + /// + /// identity of a point is in full mode the triplet origin_id (: usize), layer (: u8) rank_in_layer (: u32) + /// light mode only origin_id (: usize) + /// For data dump + /// 1. The value MAGICDATAP (u32) + /// 2. origin_id as a u64 + /// 3. The vector of data (the length is known from Description) + +fn dump_point(point : &Point , mode : DumpMode, + graphout : &mut io::BufWriter, dataout : &mut io::BufWriter) -> Result { + // + graphout.write(unsafe { &mem::transmute::(MAGICPOINT) } ).unwrap(); + // dump ext_id: usize , layer : u8 , rank in layer : i32 + graphout.write(unsafe { &mem::transmute::(point.get_origin_id()) } ).unwrap(); + let p_id = point.get_point_id(); + if mode == DumpMode::Full { + graphout.write(unsafe { &mem::transmute::(p_id.0) } ).unwrap(); + graphout.write(unsafe { &mem::transmute::(p_id.1) } ).unwrap(); + } +// log::debug!(" point dump {:?} {:?} ", p_id, self.get_origin_id()); + // then dump neighborhood info : nb neighbours : u32 , then list of origin_id, layer, rank_in_layer + let neighborhood = point.get_neighborhood_id(); + // in any case nb_layers are dumped with possibly 0 neighbours at a layer, but this does not occur by construction + for l in 0..neighborhood.len() { + let neighbours_at_l = &neighborhood[l]; + graphout.write(unsafe { &mem::transmute::(neighbours_at_l.len() as u8) } ).unwrap(); + for n in neighbours_at_l { // dump d_id : uszie , distance : f32, layer : u8, rank in layer : i32 + graphout.write(unsafe { &mem::transmute::(n.d_id) } ).unwrap(); + if mode == DumpMode::Full { + graphout.write(unsafe { &mem::transmute::(n.p_id.0) } ).unwrap(); + graphout.write(unsafe { &mem::transmute::()]>(n.p_id.1) } ).unwrap(); + } + graphout.write(unsafe { &mem::transmute::()]>(n.distance) } ).unwrap(); +// log::debug!(" voisins {:?} {:?} {:?}", n.p_id, n.d_id , n.distance); + } + } + // now we dump data vector! + dataout.write(unsafe { &mem::transmute::(MAGICDATAP) } ).unwrap(); + dataout.write(unsafe { &mem::transmute::(point.get_origin_id() as u64) } ).unwrap(); + let ref_v = point.get_v(); + let v_as_u8 = unsafe { std::slice::from_raw_parts(ref_v.as_ptr() as *const u8, mem::size_of::() * ref_v.len() )}; + dataout.write_all(v_as_u8).unwrap(); + // + return Ok(1); +} // end of dump for Point + + + +// +// Reload a point from a dump. +// +// The graph part is loaded from graph_in file +// the data vector itself is loaded from data_in +// +fn load_point(graph_in: &mut dyn Read, descr: &Description, + data_in: &mut dyn Read) -> io::Result<(Arc>, Vec >) > { + // + // read and check magic + let magic : u32 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&magic as *const u32) as *mut u8, ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + if magic != MAGICPOINT { + log::debug!("got instead of MAGICPOINT {:x}", magic); + return Err(io::Error::new(io::ErrorKind::Other, "bad magic at point beginning")); + } + let origin_id : DataId = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&origin_id as *const DataId) as *mut u8, + ::std::mem::size_of::())}; + graph_in.read_exact(it_slice)?; + // + let layer : u8 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&layer as *const u8) as *mut u8, + ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + // + let rank_in_l : i32 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&rank_in_l as *const i32) as *mut u8, + ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + let p_id = PointId{0: layer, 1:rank_in_l}; +// log::debug!(" point load {:?} {:?} ", p_id, origin_id); + // Now for each layer , read neighbours + let nb_layer = descr.nb_layer; + let nb_neighbours : u8 = 0; + let mut neighborhood = Vec:: >::with_capacity(NB_LAYER_MAX as usize); + for _l in 0..nb_layer { + let neighbour : Neighbour = Default::default(); + // read nb_neighbour : u8, then nb_neighbours times identity(depends on Full or Light) distance : f32 + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&nb_neighbours as *const u8) as *mut u8, + ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + let mut neighborhood_l : Vec = Vec::with_capacity(nb_neighbours as usize); + for _j in 0..nb_neighbours { + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&neighbour.d_id as *const DataId) as *mut u8, + ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + if descr.dumpmode == 1 { + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&neighbour.p_id.0 as *const u8) as *mut u8, + ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&neighbour.p_id.1 as *const i32) as *mut u8, + ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + } + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&neighbour.distance as *const f32) as *mut u8, + ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + // log::debug!(" voisins load {:?} {:?} {:?} ", neighbour.p_id, neighbour.d_id , neighbour.distance); + // now we have a new neighbour, we must really fill neighbourhood info, so it means going from Neighbour to PointWithOrder + neighborhood_l.push(neighbour); + } + neighborhood.push(neighborhood_l); + } + for _l in nb_layer..NB_LAYER_MAX { + neighborhood.push(Vec::::new()); + } + // + // construct a point from data_in + // + let magic : u32 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut( (&magic as *const u32) as *mut u8, + ::std::mem::size_of::() )}; + data_in.read_exact(it_slice)?; + assert_eq!(magic, MAGICDATAP, "magic not equal to MAGICDATAP in load_point, point_id : {:?} ", origin_id); + // read origin id + let origin_id_data : usize = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut( (&origin_id_data as *const usize) as *mut u8, ::std::mem::size_of::() )}; + data_in.read_exact(it_slice)?; + assert_eq!(origin_id, origin_id_data, "origin_id incoherent between graph and data"); + // now read data + let mut v : Vec = Vec::with_capacity(descr.dimension); + data_in.read_exact(unsafe { ::std::slice::from_raw_parts_mut(v.as_mut_ptr() as *mut u8, ::std::mem::size_of::() * descr.dimension)})?; + unsafe { v.set_len(descr.dimension);}; + let point = Point::::new(&v, origin_id as usize, p_id); + log::trace!("load_point origin {:?} allocated size {:?}, dim {:?}", origin_id, point.get_v().len(), descr.dimension); + // + return Ok((Arc::new(point), neighborhood)); +} // end of load_point + + + +// a dump of id and data vector of point in file suffixed by hnsw.data +#[allow(unused)] +fn dump_point_data(point : &Point, out : &mut io::BufWriter) { + out.write(unsafe { &mem::transmute::(MAGICDATAP) } ).unwrap(); + out.write(unsafe { &mem::transmute::(point.get_origin_id() as u64) } ).unwrap(); + let ref_v = point.get_v(); + let v_as_u8 = unsafe { std::slice::from_raw_parts(ref_v.as_ptr() as *const u8, mem::size_of::() * ref_v.len()) }; + out.write_all(v_as_u8).unwrap(); +} // end of dump_point_data + + +// +// dump and load of PointIndexation +// =================================== +// +// +// nb_layer : 8 +// a magick at each Layer : u32 +// . number of points in layer (usize), +// . list of point of layer +// dump entry point +// +impl HnswIO for PointIndexation { + fn dump(&self, mode : DumpMode, graphout : &mut io::BufWriter, dataout : &mut io::BufWriter) -> Result { + // dump max_layer + let layers = self.points_by_layer.read(); + let nb_layer = layers.len() as u8; + graphout.write(unsafe { &mem::transmute::(nb_layer) } ).unwrap(); + // dump layers from lower (most populatated to higher level) + for i in 0..layers.len() { + let nb_point = layers[i].len(); + log::debug!("dumping layer {:?}, nb_point {:?}", i, nb_point); + graphout.write(unsafe { &mem::transmute::(MAGICLAYER) } ).unwrap(); + graphout.write(unsafe { &mem::transmute::(nb_point) } ).unwrap(); + for j in 0..layers[i].len() { + assert_eq!(layers[i][j].get_point_id() , PointId{0: i as u8,1:j as i32 }); + dump_point(&layers[i][j], mode, graphout, dataout)?; + } + } + // dump id of entry point + let ep_read = self.entry_point.read(); + assert!(ep_read.is_some()); + let ep = ep_read.as_ref().unwrap(); + graphout.write(unsafe { &mem::transmute::()] >(ep.get_origin_id()) } ).unwrap(); + let p_id = ep.get_point_id(); + if mode == DumpMode::Full { + graphout.write(unsafe { &mem::transmute::(p_id.0) } ).unwrap(); + graphout.write(unsafe { &mem::transmute::(p_id.1) } ).unwrap(); + } + log::info!("dumped entry_point origin_d {:?}, p_id {:?} ", ep.get_origin_id(), p_id); + // + Ok(1) + } // end of dump for PointIndexation +} // end of impl HnswIO + + +fn load_point_indexation(graph_in: &mut dyn Read, + descr : &Description, + data_in: &mut dyn Read) -> io::Result > { + // + log::debug!(" in load_point_indexation"); + let mut points_by_layer : Vec> > >= Vec::with_capacity(NB_LAYER_MAX as usize); + let mut neighbourhood_map : HashMap> > = HashMap::new(); + // load max layer + let nb_layer : u8 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&nb_layer as *const u8) as *mut u8, ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + log::debug!("nb layer {:?}", nb_layer); + if nb_layer > NB_LAYER_MAX { + return Err(io::Error::new(io::ErrorKind::Other, "inconsistent number of layErrers")); + } + // + let mut nb_points_loaded : usize = 0; + // + for l in 0..nb_layer as usize { + // read and check magic + log::debug!("loading layer {:?}", l); + let magic : u32 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&magic as *const u32) as *mut u8, ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + if magic != MAGICLAYER { + return Err(io::Error::new(io::ErrorKind::Other, "bad magic at layer beginning")); + } + let nbpoints : usize = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&nbpoints as *const usize) as *mut u8, ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + log::debug!(" layer {:?} , nb points {:?}", l , nbpoints); + let mut vlayer : Vec>> = Vec::with_capacity(nbpoints); + for r in 0..nbpoints { + // load graph and data part of point. Points are dumped in the same order. + let load_point_res = load_point(graph_in, descr, data_in)?; + let point = load_point_res.0; + let p_id = point.get_point_id(); + // some checks + assert_eq!(l, p_id.0 as usize); + if r != p_id.1 as usize { + log::debug!("\n\n origin= {:?}, p_id = {:?}", point.get_origin_id(), p_id); + log::debug!("storing at l {:?}, r {:?}", l, r); + } + assert_eq!(r , p_id.1 as usize); + // store neoghbour info of this point + neighbourhood_map.insert(p_id, load_point_res.1); + vlayer.push(Arc::clone(&point)); + } + points_by_layer.push(vlayer); + nb_points_loaded += nbpoints; + } + // at this step all points are loaded , but without their neighbours fileds are not yet initialized + for (p_id , neighbours) in &neighbourhood_map { + let point = &points_by_layer[p_id.0 as usize][p_id.1 as usize]; + for l in 0..neighbours.len() { + for n in &neighbours[l] { + let n_point = &points_by_layer[n.p_id.0 as usize][n.p_id.1 as usize]; + // now n_point is the Arc corresponding to neighbour n of point, + // construct a corresponding PointWithOrder + let n_pwo = PointWithOrder::::new(&Arc::clone(&n_point), n.distance); + point.neighbours.write()[l].push(Arc::new(n_pwo)); + } // end of for n + // must sort + point.neighbours.write()[l].sort_unstable(); + } // end of for l + } // end loop in neighbourhood_map + // + // get id of entry_point + // load entry point + log::info!("\n end of layer loading, allocating PointIndexation, nb points loaded {:?}", nb_points_loaded); + // + let origin_id : usize = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&origin_id as *const DataId) as *mut u8, ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + // + let layer : u8 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&layer as *const u8) as *mut u8, ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + // + let rank_in_l : i32 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut((&rank_in_l as *const i32) as *mut u8, ::std::mem::size_of::() )}; + graph_in.read_exact(it_slice)?; + // + log::info!("found entry point, origin_id {:?} , layer {:?}, rank in layer {:?} ", origin_id, layer, rank_in_l); + let entry_point = Arc::clone(&points_by_layer[layer as usize][rank_in_l as usize]); + log::info!(" loaded entry point, origin_id {:} p_id {:?}", entry_point.get_origin_id(),entry_point.get_point_id()); + + // + let point_indexation = PointIndexation { + max_nb_connection : descr.max_nb_connection as usize, + max_layer : NB_LAYER_MAX as usize, + points_by_layer : Arc::new(RwLock::new(points_by_layer)), + layer_g : LayerGenerator::new(descr.max_nb_connection as usize , NB_LAYER_MAX as usize), + nb_point : Arc::new(RwLock::new(nb_points_loaded)), // CAVEAT , we should increase , the whole thing is to be able to increment graph ? + entry_point : Arc::new(RwLock::new(Some(entry_point))), + }; + // + log::debug!("\n exiting load_pointIndexation"); + Ok(point_indexation) +} // end of load_pointIndexation + + + +// +// dump and load of Hnsw +// ========================= +// +// + +impl +TypeName+Send+Sync> HnswIO for Hnsw { + fn dump(&self, mode : DumpMode, graphout : &mut io::BufWriter, dataout : &mut io::BufWriter) -> Result { + // dump description , then PointIndexation + let dumpmode : u8 = match mode { + DumpMode::Full => 1, + _ => 0, + }; + let datadim : usize = self.layer_indexed_points.get_data_dimension(); + + let description = Description { + /// value is 1 for Full 0 for Light + dumpmode : dumpmode, + max_nb_connection : self.get_max_nb_connection(), + nb_layer : self.get_max_level() as u8, + ef: self.get_ef_construction(), + nb_point: self.get_nb_point(), + dimension : datadim, + distname : self.get_distance_name(), + t_name: T::type_name(), + }; + log::debug!("dump obtained typename {:?}", T::type_name()); + description.dump(mode, graphout)?; + // We must dump a header for dataout. + dataout.write(unsafe { &mem::transmute::(MAGICDATAP) } ).unwrap(); + dataout.write(unsafe { &mem::transmute::()]>(datadim) } ).unwrap(); + // + self.layer_indexed_points.dump(mode, graphout, dataout)?; + Ok(1) + } +} // end impl block for Hnsw + + + +/// The reload is made in two steps. +/// First a call to load_description must be used to get basic information +/// about structure to reload (Typename, distance type, construction parameters). +/// Cf ```fn load_description(io_in: &mut dyn Read) -> io::Result``` +/// +pub fn load_hnsw+TypeName+Default+Send+Sync>(graph_in: &mut dyn Read, + description: &Description, + data_in : &mut dyn Read) -> io::Result > { + // In datafile , we must read MAGICDATAP and dimension and check + let magic : u32 = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut( (&magic as *const u32) as *mut u8, + ::std::mem::size_of::() )}; + data_in.read_exact(it_slice)?; + assert_eq!(magic, MAGICDATAP, "magic not equal to MAGICDATAP in load_point"); + let dimension : usize = 0; + let it_slice = unsafe {::std::slice::from_raw_parts_mut( (&dimension as *const usize) as *mut u8, + ::std::mem::size_of::() )}; + data_in.read_exact(it_slice)?; + assert_eq!(dimension, description.dimension, "data dimension incoherent {:?} {:?} ", + dimension, description.dimension); + // + let _mode = description.dumpmode; + let distname = description.distname.clone(); + // We must ensure that the distance stored matches the one asked for in loading hnsw + // for that we check for short names equality stripping + log::debug!("distance = {:?}", distname); + let d_type_name = D::type_name(); + let v: Vec<&str> = d_type_name.rsplit_terminator("::").collect(); + for s in v { + log::debug!(" distname part {:?}", s); + } + if d_type_name != distname { + let mut errmsg = String::from("error in distances : dumped distance is : "); + errmsg.push_str(&distname); + errmsg.push_str(" asked distance in loading is : "); + errmsg.push_str(&d_type_name); + return Err(io::Error::new(io::ErrorKind::Other, errmsg)); + } + let t_type = description.t_name.clone(); + log::debug!("T type name = {:?}", t_type); + let layer_point_indexation = load_point_indexation(graph_in, &description, data_in)?; + let data_dim = layer_point_indexation.get_data_dimension(); + // + let hnsw : Hnsw:: = Hnsw{ max_nb_connection : description.max_nb_connection as usize, + ef_construction : description.ef, + extend_candidates : true, + keep_pruned: false, + max_layer: description.nb_layer as usize, + layer_indexed_points: layer_point_indexation, + data_dimension : data_dim, + dist_f: D::default(), + searching : false, + } ; + // + Ok(hnsw) +} // end of load_hnsw + + + + +//=============================================================================================================== + +#[cfg(test)] + + +mod tests { + +use super::*; +use crate::dist; + + +use std::fs::OpenOptions; +use std::io::{BufReader}; +pub use crate::dist::*; +pub use crate::api::AnnT; + +use rand::distributions::{Distribution, Uniform}; +#[test] + +fn test_dump_reload() { + println!("\n\n test_dump_reload"); + let _res = simple_logger::init(); + // generate a random test + let mut rng = rand::thread_rng(); + let unif = Uniform::::new(0.,1.); + let nbcolumn = 1000; + let nbrow = 10; + let mut xsi; + let mut data = Vec::with_capacity(nbcolumn); + for j in 0..nbcolumn { + data.push(Vec::with_capacity(nbrow)); + for _ in 0..nbrow { + xsi = unif.sample(&mut rng); + data[j].push(xsi); + } + } + // define hnsw + let ef_construct= 25; + let nb_connection = 10; + let hnsw = Hnsw::::new(nb_connection, nbcolumn, 16, ef_construct, dist::DistL1{}); + for i in 0..data.len() { + hnsw.insert((&data[i], i)); + } + hnsw.dump_layer_info(); + // dump in a file + let fname = String::from("dumpreloadtest"); + let _res = hnsw.file_dump(&fname); + // This will dump in 2 files named dumpreloadtest.hnsw.graph and dumpreloadtest.hnsw.data + // + // reload + log::debug!("\n\n hnsw reload"); + // we will need a procedural macro to get from distance name to its instanciation. + // from now on we test with DistL1 + let graphfname = String::from("dumpreloadtest.hnsw.graph"); + let graphpath = PathBuf::from(graphfname); + let graphfileres = OpenOptions::new().read(true).open(&graphpath); + if graphfileres.is_err() { + println!("could not open file {:?}", graphpath.as_os_str()); + panic!("could not open file".to_string()); + } + let graphfile = graphfileres.unwrap(); + // + let datafname = String::from("dumpreloadtest.hnsw.data"); + let datapath = PathBuf::from(datafname); + let datafileres = OpenOptions::new().read(true).open(&datapath); + if datafileres.is_err() { + println!("could not open file {:?}", datapath.as_os_str()); + panic!("could not open file".to_string()); + } + let datafile = datafileres.unwrap(); + // + let mut graph_in = BufReader::new(graphfile); + let mut data_in = BufReader::new(datafile); + // we need to call load_description first to get distance name + let hnsw_description = load_description(&mut graph_in).unwrap(); + let hnsw_loaded : Hnsw= load_hnsw(&mut graph_in, &hnsw_description, &mut data_in).unwrap(); + // test equality + check_equality(&hnsw_loaded, &hnsw); +} // end of test_dump_reload + +} // end module tests \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..4f13a49 --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,43 @@ +extern crate rand; + +// for logging (debug mostly, switched at compile time in cargo.toml) +extern crate log; +extern crate simple_logger; + +#[macro_use] +extern crate lazy_static; + + +pub mod hnsw; +pub mod dist; +pub mod annhdf5; +pub mod hnswio; +pub mod test; +pub mod prelude; +pub mod api; +pub mod libext; + +lazy_static! { + #[allow(dead_code)] + static ref LOG: u64 = { + let res = init_log(); + res + }; +} + +// install a logger facility +fn init_log() -> u64 { + simple_logger::init().unwrap(); + println!("\n ************** initializing logger *****************\n"); + return 1; +} + +#[cfg(test)] +mod tests { + use super::*; + #[test] + // initialize once log system for tests. + fn init_log() { + let _res = simple_logger::init(); + } +} // end of tests diff --git a/src/libext.rs b/src/libext.rs new file mode 100644 index 0000000..35d0a6a --- /dev/null +++ b/src/libext.rs @@ -0,0 +1,782 @@ +//! This file contains lib to call hnsw from julia (or any language providing a C api) +//! The AnnT trait is implemented with macros for u32, u16, u8, f32, f64 and i32. +//! an declare_myapi_type! and procudes struct HnswApif32 and so on. +//! + +use std::ptr; +use std::fs::OpenOptions; +use std::io::{BufReader}; +use std::path::{PathBuf}; + +use crate::api::*; +use crate::dist::*; +use crate::hnsw::*; +use crate::hnswio::*; + + + +// the export macro makes the macro global in crate and accecssible via crate::declare_myapi_type! +#[macro_export] +macro_rules! declare_myapi_type( + ($name:ident, $ty:ty) => ( + pub struct $name { +#[allow(dead_code)] + pub(crate) opaque: Box>, + } // end struct + impl $name { + pub fn new(arg: Box>) -> Self { + $name{ opaque:arg} + } // end new + } // end impl + ) +); + + +// declare_myapi_type!(HnswApif64, f64); +// declare_myapi_type!(HnswApif32, f32); + +/// to be able to return a vector from rust in a julia struct before converting to a julia Vector +#[repr(C)] +pub struct Vec_api { + len : i64, + ptr : *const T, +} // end struct + + + +#[repr(C)] +/// The basic Neighbour info returned by api +pub struct Neighbour_api { + /// id of neighbour + pub id : usize, + /// distance of data sent in request to this neighbour + pub d : f32, +} + + +impl From<&Neighbour> for Neighbour_api { + fn from(neighbour : &Neighbour) -> Self { + Neighbour_api{ id : neighbour.d_id, d: neighbour.distance} + } +} + + +#[repr(C)] +/// The response to a neighbour search requests +pub struct Neighbourhood_api { + pub nbgh : i64, + pub neighbours : *const Neighbour_api, +} + + + +#[repr(C)] +pub struct Neighbour_api_parsearch_answer { + /// The number of answers (o request), i.e size of both vectors nbgh and neighbours + pub nb_answer : usize, + /// for each request, we get a Neighbourhood_api + pub neighbourhoods : *const Neighbourhood_api, +} + +//===================================== f32 type ===================================== + + +// macros have been exported to the root of the crate so we do not refer to them via api:: +super::declare_myapi_type!(HnswApif32, f32); +super::declare_myapi_type!(HnswApif64, f64); + + + +//=================================================================================================== +// These are the macros to generate trait implementation for useful numeric types +#[allow(unused_macros)] + +macro_rules! generate_insert( +($function_name:ident, $api_name:ty, $type_val:ty) => ( + #[no_mangle] + pub extern "C" fn $function_name(hnsw_api : *mut $api_name, len:usize, data : *const $type_val, id : usize) { + log::trace!("entering insert, type {:?} vec len is {:?}, id : {:?} ", stringify!($type_val), len, id); + // construct vector: Rust clones and take ownership. + let data_v : Vec<$type_val>; + unsafe { + let slice = std::slice::from_raw_parts(data, len); + data_v = Vec::from(slice); + log::trace!("calling insert data"); + (*hnsw_api).opaque.insert_data(&data_v, id); + } + log::trace!("exiting insert for type {:?}", stringify!($type_val)); + } // end of insert + ) +); + + + +macro_rules! generate_parallel_insert( +($function_name:ident, $api_name:ty, $type_val:ty) => ( + #[no_mangle] + pub extern "C" fn $function_name(hnsw_api : *mut $api_name, nb_vec: usize, vec_len : usize, + datas : *mut *const $type_val, ids : *const usize) { + // + log::trace!("entering parallel_insert type {:?} , vec len is {:?}, nb_vec : {:?}", stringify!($type_val), vec_len, nb_vec); + let data_ids : Vec; + let data_ptrs : Vec<*const $type_val>; + unsafe { + let slice = std::slice::from_raw_parts(ids, nb_vec); + data_ids = Vec::from(slice); + } + // log::debug!("got ids"); + unsafe { + let slice = std::slice::from_raw_parts(datas, nb_vec); + data_ptrs = Vec::from(slice); + } + // log::debug!("got data ptrs"); + let mut data_v = Vec::>::with_capacity(nb_vec); + for i in 0..nb_vec { + unsafe { + let slice = std::slice::from_raw_parts(data_ptrs[i], vec_len); + let v = Vec::from(slice); + data_v.push(v); + } + } + // log::debug!("sending request"); + let mut request : Vec<(&Vec<$type_val>, usize)> = Vec::with_capacity(nb_vec); + for i in 0..nb_vec { + request.push((&data_v[i], data_ids[i])); + } + // + unsafe { (*hnsw_api).opaque.parallel_insert_data(&request); }; + log::trace!("exiting parallel_insert"); + } // end of parallel_insert + ) +); + + +macro_rules! generate_search_neighbours( +($function_name:ident, $api_name:ty, $type_val:ty) => ( + #[no_mangle] + pub extern "C" fn $function_name(hnsw_api : *const $api_name, len:usize, data : *const $type_val, + knbn : usize, ef_search : usize) -> *const Neighbourhood_api { + // + log::trace!("entering search_neighbours type {:?}, vec len is {:?}, id : {:?} ef_search {:?}", stringify!($type_val), len, knbn, ef_search); + let data_v : Vec<$type_val>; + let neighbours : Vec; + unsafe { + let slice = std::slice::from_raw_parts(data, len); + data_v = Vec::from(slice); + log::trace!("calling search neighbours {:?}", data_v); + neighbours = (*hnsw_api).opaque.search_neighbours(&data_v, knbn, ef_search); + } + let neighbours_api : Vec = neighbours.iter().map(|n| Neighbour_api::from(n)).collect(); + log::trace!(" got nb neighbours {:?}", neighbours_api.len()); + // for i in 0..neighbours_api.len() { + // println!(" id {:?} dist : {:?} ", neighbours_api[i].id, neighbours_api[i].d); + // } + let nbgh_i = neighbours.len() as i64; + let neighbours_ptr = neighbours_api.as_ptr(); + std::mem::forget(neighbours_api); + let answer = Neighbourhood_api { + nbgh : nbgh_i, + neighbours : neighbours_ptr, + }; + log::trace!("search_neighbours returning nb neighbours {:?} id ptr {:?} ", nbgh_i, neighbours_ptr); + Box::into_raw(Box::new(answer)) + } + ) +); + + + + +macro_rules! generate_parallel_search_neighbours( +($function_name:ident, $api_name:ty, $type_val:ty) => ( + #[no_mangle] + /// search nb_vec of size vzc_len. The the searches will be done in // as far as possible. + pub extern "C" fn $function_name(hnsw_api : *const $api_name, nb_vec : usize, vec_len :i64, + data : *mut *const $type_val, knbn : usize, ef_search : usize) -> *const Vec_api { + // + // must build a Vec to build request + log::trace!("recieving // search request for type: {:?} with {:?} vectors", stringify!($type_val), nb_vec); + let neighbours : Vec >; + let mut data_v = Vec::>::with_capacity(nb_vec); + unsafe { + let slice = std::slice::from_raw_parts(data, nb_vec); + let ptr_list : Vec<*const $type_val> = Vec::from(slice); + for i in 0..nb_vec { + let slice_i = std::slice::from_raw_parts(ptr_list[i], vec_len as usize); + let v = Vec::from(slice_i); + data_v.push(v); + } + // log::debug!(" reconstructed input vectors"); + neighbours = (*hnsw_api).opaque.parallel_search_neighbours(&data_v, knbn, ef_search); + } + // construct a vector of Neighbourhood_api + // reverse work, construct 2 arrays, one vector of Neighbours, and one vectors of number of returned neigbours by input a vector. + let mut neighbour_lists = Vec::::with_capacity(nb_vec); + for v in neighbours { + let neighbours_api : Vec = v.iter().map(|n| Neighbour_api::from(n)).collect(); + let nbgh = neighbours_api.len(); + let neighbours_api_ptr = neighbours_api.as_ptr(); + std::mem::forget(neighbours_api); + let v_answer = Neighbourhood_api { + nbgh : nbgh as i64, + neighbours: neighbours_api_ptr, + }; + neighbour_lists.push(v_answer); + } + log::trace!(" reconstructed output pointers to vectors"); + let neighbour_lists_ptr = neighbour_lists.as_ptr(); + std::mem::forget(neighbour_lists); + let answer = Vec_api:: { + len : nb_vec as i64, + ptr : neighbour_lists_ptr, + }; + Box::into_raw(Box::new(answer)) + } // end of parallel_search_neighbours_f32 for HnswApif32 + ) +); + +#[allow(unused_macros)] +macro_rules! generate_file_dump( + ($function_name:ident, $api_name:ty, $type_val:ty) => ( + #[no_mangle] + pub extern "C" fn $function_name(hnsw_api : *const $api_name, namelen : usize, filename :*const u8) -> i64 { + log::info!("receiving request for file dump"); + let slice = unsafe { std::slice::from_raw_parts(filename, namelen) } ; + let fstring = String::from_utf8_lossy(slice).into_owned(); + let res = unsafe { (*hnsw_api).opaque.file_dump(&fstring) } ; + if res.is_ok() { + return 1; + } + else { return -1; } + } // end of function_name + ) +); + + +#[allow(unused_macros)] +macro_rules! generate_loadhnsw( + ($function_name:ident, $api_name:ty, $type_val:ty) => ( + #[no_mangle] + pub extern "C" fn $function_name(flen : usize, name : *const u8) -> *const $api_name { + let _res = simple_logger::init(); + let slice = unsafe { std::slice::from_raw_parts(name, flen)} ; + let filename = String::from_utf8_lossy(slice).into_owned(); + // + let buffers = make_readers(&filename); + let mut graph_in = buffers.0; + let mut data_in = buffers.1; + // we need to call load_description first to get distance name + let hnsw_description = load_description(&mut graph_in).unwrap(); + // here we need a macro to dispatch + let hnsw_loaded_res = load_hnsw::<$type_val, crate::mapdist_t!("DistL1")>(&mut graph_in, &hnsw_description, &mut data_in); + if let Ok(hnsw_loaded) = hnsw_loaded_res { + let api = <$api_name>::new(Box::new(hnsw_loaded)); + return Box::into_raw(Box::new(api)); + } + else { + log::warn!("an error occured, could not reload data from {:?}", filename); + } + return ptr::null(); + } // end of load_hnswdump_ + ) +); + + +generate_loadhnsw!(load_hnswdump_f32, HnswApif32, f32); +generate_loadhnsw!(load_hnswdump_i32, HnswApii32, i32); +generate_loadhnsw!(load_hnswdump_u32, HnswApiu32, u32); +generate_loadhnsw!(load_hnswdump_u16, HnswApiu16, u16); +generate_loadhnsw!(load_hnswdump_u8, HnswApiu8, u8); + +//=============== implementation for i32 + + + +#[no_mangle] +pub extern "C" fn init_hnsw_f32(max_nb_conn : usize, ef_const:usize, namelen: usize, cdistname : *const u8) -> *const HnswApif32 { + let _res = simple_logger::init(); + log::info!("entering init_hnsw_f32"); + let slice = unsafe { std::slice::from_raw_parts(cdistname, namelen)} ; + let dname = String::from_utf8_lossy(slice).into_owned(); + // map distname to sthg. This whole block will go to a macro + match dname.as_str() { + "DistL1" => { + log::info!(" received DistL1"); + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL1{}); + let api = HnswApif32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + "DistL2" => { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL2{}); + let api = HnswApif32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + "DistDot" => { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistDot{}); + let api = HnswApif32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + "DistHellinger" => { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistHellinger{}); + let api = HnswApif32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + "DistJeffreys" => { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistJeffreys{}); + let api = HnswApif32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + _ => { + log::warn!("init_hnsw_f32 received unknow distance {:?} ", dname); + let p = ptr::null::< HnswApif32 >(); + return p; + } + } // znd match +} // end of init_hnsw_f32 + + + +#[no_mangle] +pub extern "C" fn init_hnsw_ptrdist_f32(max_nb_conn : usize, ef_const:usize, c_func : extern "C" fn(*const f32, *const f32, u64) -> f32 ) -> *const HnswApif32 { + let _res = simple_logger::init(); + log::info!("init_ hnsw_ptrdist: ptr {:?}", c_func); + let c_dist = DistCFFI::::new(c_func); + let h = Hnsw:: >::new(max_nb_conn, 10000, 16, ef_const, c_dist); + let api = HnswApif32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); +} + +#[no_mangle] +pub extern "C" fn insert_f32(hnsw_api : *mut HnswApif32, len:usize, data : *const f32, id : usize) { + log::trace!("entering insert_f32, vec len is {:?}, id : {:?} ", len, id); + // construct vector: Rust clones and take ownership. + let data_v : Vec; + unsafe { + let slice = std::slice::from_raw_parts(data, len); + data_v = Vec::from(slice); + // log::debug!("calling insert data"); + (*hnsw_api).opaque.insert_data(&data_v, id); + } + // log::trace!("exiting insert_f32"); +} // end of insert_f32 + + +#[no_mangle] +pub extern "C" fn parallel_insert_f32(hnsw_api : *mut HnswApif32, nb_vec: usize, vec_len : usize, + datas : *mut *const f32, ids : *const usize) { + // + // log::debug!("entering parallel_insert_f32 , vec len is {:?}, nb_vec : {:?}", vec_len, nb_vec); + let data_ids : Vec; + let data_ptrs : Vec<*const f32>; + unsafe { + let slice = std::slice::from_raw_parts(ids, nb_vec); + data_ids = Vec::from(slice); + } + unsafe { + let slice = std::slice::from_raw_parts(datas, nb_vec); + data_ptrs = Vec::from(slice); + } + // log::debug!("got data ptrs"); + let mut data_v = Vec::>::with_capacity(nb_vec); + for i in 0..nb_vec { + unsafe { + let slice = std::slice::from_raw_parts(data_ptrs[i], vec_len); + let v = Vec::from(slice); + data_v.push(v); + } + } + // log::debug!("sending request"); + let mut request : Vec<(&Vec, usize)> = Vec::with_capacity(nb_vec); + for i in 0..nb_vec { + request.push((&data_v[i], data_ids[i])); + } + // + unsafe { (*hnsw_api).opaque.parallel_insert_data(&request); }; + log::trace!("exiting parallel_insert"); +} // end of parallel_insert_f32 + + + + +#[no_mangle] +pub extern "C" fn search_neighbours_f32(hnsw_api : *const HnswApif32, len:usize, data : *const f32, + knbn : usize, ef_search : usize) -> *const Neighbourhood_api { + // + log::trace!("entering search_neighbours , vec len is {:?}, id : {:?} ef_search {:?}", len, knbn, ef_search); + let data_v : Vec; + let neighbours : Vec; + unsafe { + let slice = std::slice::from_raw_parts(data, len); + data_v = Vec::from(slice); + log::trace!("calling search neighbours {:?}", data_v); + neighbours = (*hnsw_api).opaque.search_neighbours(&data_v, knbn, ef_search); + } + let neighbours_api : Vec = neighbours.iter().map(|n| Neighbour_api::from(n)).collect(); + log::trace!(" got nb neighbours {:?}", neighbours_api.len()); + // for i in 0..neighbours_api.len() { + // println!(" id {:?} dist : {:?} ", neighbours_api[i].id, neighbours_api[i].d); + // } + let nbgh_i = neighbours.len() as i64; + let neighbours_ptr = neighbours_api.as_ptr(); + std::mem::forget(neighbours_api); + let answer = Neighbourhood_api { + nbgh : nbgh_i, + neighbours : neighbours_ptr, + }; + log::trace!("search_neighbours returning nb neighbours {:?} id ptr {:?} ", nbgh_i, neighbours_ptr); + Box::into_raw(Box::new(answer)) +} + // end of search_neighbours for HnswApif32 + + + + +generate_parallel_search_neighbours!(parallel_search_neighbours_f32, HnswApif32, f32); +generate_file_dump!(file_dump_f32, HnswApif32, f32); + + + +//=============== implementation for i32 + +#[no_mangle] +pub extern "C" fn init_hnsw_i32(max_nb_conn : usize, ef_const:usize, namelen: usize, cdistname : *const u8) -> *const HnswApii32 { + let _res = simple_logger::init(); + log::info!("entering init_hnsw_i32"); + let slice = unsafe { std::slice::from_raw_parts(cdistname, namelen)} ; + let dname = String::from_utf8_lossy(slice); + // map distname to sthing. This whole block will go to a macro + if dname == "DistL1" { + log::info!(" received DistL1"); + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL1{}); + let api = HnswApii32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + else if dname == "DistL2" { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL2{}); + let api = HnswApii32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + let p = ptr::null::< HnswApii32 >(); + p +} // end of init_hnsw_i32 + + + + +#[no_mangle] +pub extern "C" fn init_hnsw_ptrdist_i32(max_nb_conn : usize, ef_const:usize, + c_func : extern "C" fn(*const i32, *const i32, u64) -> f32 ) -> *const HnswApii32 { + let _res = simple_logger::init(); + log::debug!("init_ hnsw_ptrdist: ptr {:?}", c_func); + let c_dist = DistCFFI::::new(c_func); + let h = Hnsw:: >::new(max_nb_conn, 10000, 16, ef_const, c_dist); + let api = HnswApii32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); +} + + +//==generation of function for i32 + +super::declare_myapi_type!(HnswApii32, i32); + +generate_insert!(insert_i32, HnswApii32, i32); +generate_parallel_insert!(parallel_insert_i32, HnswApii32, i32); +generate_search_neighbours!(search_neighbours_i32, HnswApii32, i32); +generate_parallel_search_neighbours!(parallel_search_neighbours_i32, HnswApii32, i32); +generate_file_dump!(file_dump_i32, HnswApii32, i32); + +//========== generation for u32 + + +#[no_mangle] +pub extern "C" fn init_hnsw_u32(max_nb_conn : usize, ef_const:usize, namelen: usize, cdistname : *const u8) -> *const HnswApiu32 { + let _res = simple_logger::init(); + log::debug!("entering init_hnsw_u32"); + let slice = unsafe { std::slice::from_raw_parts(cdistname, namelen)} ; + let dname = String::from_utf8_lossy(slice); + // map distname to sthg. This whole block will go to a macro + if dname == "DistL1" { + log::debug!(" received DistL1"); + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL1{}); + let api = HnswApiu32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + else if dname == "DistL2" { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL2{}); + let api = HnswApiu32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + else if dname == "DistJaccard" { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistJaccard{}); + let api = HnswApiu32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + let p = ptr::null::< HnswApiu32 >(); + p +} // end of init_hnsw_i32 + + + +#[no_mangle] +pub extern "C" fn init_hnsw_ptrdist_u32(max_nb_conn : usize, ef_const:usize, + c_func : extern "C" fn(*const u32, *const u32, u64) -> f32 ) -> *const HnswApiu32 { + let _res = simple_logger::init(); + log::info!("init_ hnsw_ptrdist: ptr {:?}", c_func); + let c_dist = DistCFFI::::new(c_func); + let h = Hnsw:: >::new(max_nb_conn, 10000, 16, ef_const, c_dist); + let api = HnswApiu32{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); +} + + +super::declare_myapi_type!(HnswApiu32, u32); + +generate_insert!(insert_u32, HnswApiu32, u32); +generate_parallel_insert!(parallel_insert_u32, HnswApiu32, u32); +generate_search_neighbours!(search_neighbours_u32, HnswApiu32, u32); +generate_parallel_search_neighbours!(parallel_search_neighbours_u32, HnswApiu32, u32); +generate_file_dump!(file_dump_u32, HnswApiu32, u32); + + + +//============== generation of function for u16 ===================== + +super::declare_myapi_type!(HnswApiu16, u16); + + +#[no_mangle] +pub extern "C" fn init_hnsw_u16(max_nb_conn : usize, ef_const:usize, namelen: usize, cdistname : *const u8) -> *const HnswApiu16 { + let _res = simple_logger::init(); + log::info!("entering init_hnsw_u16"); + let slice = unsafe { std::slice::from_raw_parts(cdistname, namelen)} ; + let dname = String::from_utf8_lossy(slice); + // map distname to sthg. This whole block will go to a macro + if dname == "DistL1" { + log::info!(" received DistL1"); + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL1{}); + let api = HnswApiu16{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + else if dname == "DistL2" { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL2{}); + let api = HnswApiu16{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + else if dname == "DistHamming" { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistHamming{}); + let api = HnswApiu16{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + else if dname == "DistJaccard" { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistJaccard{}); + let api = HnswApiu16{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + let p = ptr::null::< HnswApiu16 >(); + p +} // end of init_hnsw_u16 + + + +#[no_mangle] +pub extern "C" fn init_hnsw_ptrdist_u16(max_nb_conn : usize, ef_const:usize, + c_func : extern "C" fn(*const u16, *const u16, u64) -> f32 ) -> *const HnswApiu16 { + let _res = simple_logger::init(); + log::info!("init_ hnsw_ptrdist: ptr {:?}", c_func); + let c_dist = DistCFFI::::new(c_func); + let h = Hnsw:: >::new(max_nb_conn, 10000, 16, ef_const, c_dist); + let api = HnswApiu16{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); +} + + +generate_insert!(insert_u16, HnswApiu16, u16); +generate_parallel_insert!(parallel_insert_u16, HnswApiu16, u16); +generate_search_neighbours!(search_neighbours_u16, HnswApiu16, u16); +generate_parallel_search_neighbours!(parallel_search_neighbours_u16, HnswApiu16, u16); +generate_file_dump!(file_dump_u16, HnswApiu16, u16); + + + +//============== generation of function for u8 ===================== + +super::declare_myapi_type!(HnswApiu8, u8); + + +#[no_mangle] +pub extern "C" fn init_hnsw_u8(max_nb_conn : usize, ef_const:usize, namelen: usize, cdistname : *const u8) -> *const HnswApiu8 { + let _res = simple_logger::init(); + log::debug!("entering init_hnsw_u8"); + let slice = unsafe { std::slice::from_raw_parts(cdistname, namelen)} ; + let dname = String::from_utf8_lossy(slice); + // map distname to sthg. This whole block will go to a macro + if dname == "DistL1" { + log::info!(" received DistL1"); + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL1{}); + let api = HnswApiu8{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + else if dname == "DistL2" { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistL2{}); + let api = HnswApiu8{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + else if dname == "DistHamming" { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistHamming{}); + let api = HnswApiu8{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + else if dname == "DistJaccard" { + let h = Hnsw::::new(max_nb_conn, 10000, 16, ef_const, DistJaccard{}); + let api = HnswApiu8{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); + } + let p = ptr::null::< HnswApiu8 >(); + p +} // end of init_hnsw_u16 + + + +#[no_mangle] +pub extern "C" fn init_hnsw_ptrdist_u8(max_nb_conn : usize, ef_const:usize, + c_func : extern "C" fn(*const u8, *const u8, u64) -> f32 ) -> *const HnswApiu8 { + let _res = simple_logger::init(); + log::info!("init_ hnsw_ptrdist: ptr {:?}", c_func); + let c_dist = DistCFFI::::new(c_func); + let h = Hnsw:: >::new(max_nb_conn, 10000, 16, ef_const, c_dist); + let api = HnswApiu8{opaque: Box::new(h)}; + return Box::into_raw(Box::new(api)); +} + +generate_insert!(insert_u8, HnswApiu8, u8); +generate_parallel_insert!(parallel_insert_u8, HnswApiu8, u8); +generate_search_neighbours!(search_neighbours_u8, HnswApiu8, u8); +generate_parallel_search_neighbours!(parallel_search_neighbours_u8, HnswApiu8, u8); +generate_file_dump!(file_dump_u8, HnswApiu8, u8); + + + +//=========================== dump restore functions + +/// This structure provides a light description of the graph to be passed to C compatible languages. +#[repr(C)] +pub struct DescriptionFFI { + /// value is 1 for Full 0 for Light + pub dumpmode : u8, + /// max number of connections in layers != 0 + pub max_nb_connection : u8, + /// number of observed layers + pub nb_layer : u8, + /// search parameter + pub ef: usize, + /// total number of points + pub nb_point: usize, + /// dimension of data vector + pub data_dimension : usize, + /// length and pointer on dist name + pub distname_len: usize, + pub distname : *const u8, + /// T typename + pub t_name_len : usize, + pub t_name : *const u8, +} + +impl DescriptionFFI { + pub fn new() -> Self { + DescriptionFFI { + dumpmode : 0, + max_nb_connection : 0, + nb_layer : 0, + ef : 0, + nb_point : 0, + data_dimension : 0, + distname_len : 0, + distname : ptr::null(), + t_name_len : 0, + t_name : ptr::null() + } + } // end of new +} + + +/// returns a const pointer to a DescriptionFFI from a dump file, given filename length and pointer (*const u8) +#[no_mangle] +pub extern "C" fn load_hnsw_description(flen : usize, name : *const u8) -> *const DescriptionFFI { + // opens file + let slice = unsafe { std::slice::from_raw_parts(name, flen)} ; + let filename = String::from_utf8_lossy(slice).into_owned(); + let fpath = PathBuf::from(filename); + let fileres = OpenOptions::new().read(true).open(&fpath); + // + let mut ffi_description = DescriptionFFI::new(); + match fileres { + Ok(file) => { + // + let mut bufr = BufReader::with_capacity(10000000 , file); + let res = load_description(&mut bufr); + if let Ok(description) = res { + + let distname = String::clone(&description.distname); + let distname_ptr = distname.as_ptr(); + let distname_len = distname.len(); + std::mem::forget(distname); + + let t_name = String::clone(&description.t_name); + let t_name_ptr = t_name.as_ptr(); + let t_name_len = t_name.len(); + std::mem::forget(t_name); + + ffi_description.dumpmode = 1; // CAVEAT + ffi_description.max_nb_connection = description.max_nb_connection; + ffi_description.nb_layer = description.nb_layer; + ffi_description.ef = description.ef; + ffi_description.data_dimension = description.dimension; + ffi_description.distname_len = distname_len; + ffi_description.distname = distname_ptr; + ffi_description.t_name_len = t_name_len; + ffi_description.t_name = t_name_ptr; + return Box::into_raw(Box::new(ffi_description)); + } + else { + println!("could not get descrption of hnsw from file {:?} ", fpath.as_os_str()); + return ptr::null(); + } + }, + Err(_e) => { + println!("could not open file {:?}", fpath.as_os_str()); + return ptr::null(); + } + } +} // end of load_hnsw_description + + +// This function provides the 2 buffers needed to read the dump of a graph +// in file basename generated by +fn make_readers(basename: &String) -> (BufReader, BufReader) { + let mut graphfname = String::clone(basename); + graphfname.push_str(".hnsw.graph"); + let graphpath = PathBuf::from(graphfname); + let graphfileres = OpenOptions::new().read(true).open(&graphpath); + if graphfileres.is_err() { + println!("could not open file {:?}", graphpath.as_os_str()); + panic!("could not open file".to_string()); + } + let graphfile = graphfileres.unwrap(); + // + let mut datafname = String::clone(basename); + datafname.push_str(".hnsw.data"); + let datapath = PathBuf::from(datafname); + let datafileres = OpenOptions::new().read(true).open(&datapath); + if datafileres.is_err() { + println!("could not open file {:?}", datapath.as_os_str()); + panic!("could not open file".to_string()); + } + let datafile = datafileres.unwrap(); + // + let graph_in = BufReader::new(graphfile); + let data_in = BufReader::new(datafile); + (graph_in, data_in) +} + + + diff --git a/src/prelude.rs b/src/prelude.rs new file mode 100644 index 0000000..b69d134 --- /dev/null +++ b/src/prelude.rs @@ -0,0 +1,5 @@ +// gathers modules to include + +pub use crate::hnsw::*; +pub use crate::dist::*; +pub use crate::annhdf5::*; \ No newline at end of file diff --git a/src/test.rs b/src/test.rs new file mode 100644 index 0000000..de26d2a --- /dev/null +++ b/src/test.rs @@ -0,0 +1,316 @@ +//! some testing utilities +//! run with to get output statistics : cargo test --release -- --nocapture --test test_parallel +//! serial test corresponds to random-10nn-euclidean(k=10) +//! par test corresponds to random data in 25 dimensions k = 10, dist Cosine + +#![allow(dead_code)] +use rand::distributions::{Uniform}; +use rand::prelude::*; + +use ndarray::*; +use skiplist::OrderedSkipList; + +use crate::hnsw; + +#[allow(unused_imports)] // necessary for rls +use crate::dist; + +pub use self::hnsw::*; + + +pub fn gen_random_vector_f32(nbrow :usize) -> Vec { + let mut rng = thread_rng(); + let unif = Uniform::::new(0.,1.); + let mut vec = Vec::::with_capacity(nbrow); + for _ in 0..nbrow { + let xsi = rng.sample(unif); + vec.push(xsi); + } + vec +} + + + +/// return nbcolumn vectors of dimension nbrow +pub fn gen_random_matrix_f32(nbrow:usize, nbcolumn:usize) -> Vec> { + let mut rng = thread_rng(); + let unif = Uniform::::new(0.,1.); + + let mut xsi; + let mut data = Vec::with_capacity(nbcolumn); + for j in 0..nbcolumn { + data.push(Vec::with_capacity(nbrow)); + for _ in 0..nbrow { + xsi = rng.sample(unif); + data[j].push(xsi); + } + } + return data; +} + +/// return brute force distance matrix +pub fn brute_force_1(data:&Vec>, dist: &PointDistance ) -> Array2 { + let nbelem = data.len(); + + let mut dist_matrix = Array2::::zeros((nbelem, nbelem)); + + for i in 0..nbelem { + for j in 0..i { + dist_matrix[[i,j] ] = dist.eval(&data[i], &data[j]); + } + } + dist_matrix +} + + +fn brute_force_neighbours(nb_neighbours: usize, refdata: &PointIndexation, distance: PointDistance, data : &Vec) -> OrderedSkipList { + + let mut neighbours = OrderedSkipList::::with_capacity(refdata.get_nb_point()); + + let mut ptiter = refdata.into_iter(); + let mut more = true; + while more { + if let Some(point) = ptiter.next() { + let dist_p = distance.eval(data, point.get_v()); + let ordered_point = PointIdWithOrder::new(point.get_point_id(), dist_p); +// log::debug!(" brute force inserting {:?}", ordered_point); + if neighbours.len() < nb_neighbours { + neighbours.insert(ordered_point); + } + else { + neighbours.insert(ordered_point); + neighbours.pop_back(); + } + } + else { + more = false; + } + } // end while + neighbours +} // end of brute_force_2 + +//================================================================================================ + + +#[cfg(test)] +mod tests { + +use super::*; +use std::time::Duration; +use cpu_time::ProcessTime; + +use dist::l2_normalize; + +#[test] +fn test_serial() { + // + let _res = simple_logger::init(); + // + let nb_elem = 1000; + let dim = 10; + let knbn = 10; + let ef = 20; + let parallel = true; + // + println!("\n\n test_serial nb_elem {:?}", nb_elem); + // + let data = gen_random_matrix_f32(dim, nb_elem); + let data_with_id = data.iter().zip(0..data.len()).collect(); + + let ef_c = 400; + let max_nb_connection = 32; + let nb_layer = 16.min((nb_elem as f32).ln().trunc() as usize); + let mut hns = Hnsw::::new(max_nb_connection, nb_elem, nb_layer, ef_c, dist::DistL1{}); + hns.set_extend_candidates(true); + hns.set_keeping_pruned(true); + let mut start = ProcessTime::now(); + if parallel { + println!("parallel insertion"); + hns.parallel_insert(&data_with_id); + } + else { + println!("serial insertion"); + for i in 0..data.len() { + hns.insert((&data[i],i)); + } + } + let mut cpu_time: Duration = start.elapsed(); + println!(" hnsw serial data insertion {:?}", cpu_time); + hns.dump_layer_info(); + println!(" hnsw data nb point inserted {:?}", hns.get_nb_point()); + // + + let nbtest=300; + let mut recalls = Vec::::with_capacity(nbtest); + let mut nb_returned = Vec::::with_capacity(nb_elem); + let mut search_times = Vec::::with_capacity(nbtest); + for _itest in 0..nbtest { + // + let mut r_vec = Vec::::with_capacity(dim); + let mut rng = thread_rng(); + let unif = Uniform::::new(0.,1.); + for _ in 0..dim { + r_vec.push(rng.sample(unif)); + } + start = ProcessTime::now(); + let brute_neighbours = brute_force_neighbours(knbn, hns.get_point_indexation() , Box::new(dist::DistL1{}) , &r_vec); + cpu_time = start.elapsed(); + if nbtest <= 100 { + println!("\n\n **************** test {:?}", _itest); + println!("\n brute force neighbours :" ); + println!("======================"); + println!(" brute force computing {:?} \n ", cpu_time); + for i in 0..brute_neighbours.len() { + let p = brute_neighbours[i].point_id; + println!(" {:?} {:?} ", p, brute_neighbours[i].dist_to_ref); + } + } + // + hns.set_searching_mode(true); + start = ProcessTime::now(); + let knn_neighbours = hns.search(&r_vec, knbn, ef); + cpu_time = start.elapsed(); + search_times.push(cpu_time.as_micros() as f32); + if nbtest <= 100 { + println!("\n\n hnsw searching {:?} \n", cpu_time); + println!("\n knn neighbours"); + println!("======================"); + for n in &knn_neighbours { + println!(" {:?} {:?} {:?} ", n.d_id, n.p_id, n.distance); + } + } + // compute recall + let knn_neighbours_dist : Vec = knn_neighbours.iter().map(|p| p.distance).collect(); + let max_dist = brute_neighbours[knbn-1].dist_to_ref; + let recall = knn_neighbours_dist.iter().filter(|d| *d <= &max_dist).count(); + if nbtest <= 100 { + println!("recall {:?}", (recall as f32)/ (knbn as f32)); + } + recalls.push(recall); + nb_returned.push(knn_neighbours.len()); + } // end on nbtest + // + // compute recall + // + + let mean_recall = (recalls.iter().sum::() as f32)/((knbn * recalls.len()) as f32); + let mean_search_time = (search_times.iter().sum::() as f32)/(search_times.len() as f32); + println!("\n mean fraction (of knbn) returned by search {:?} ", (nb_returned.iter().sum::() as f32)/ ((nb_returned.len() * knbn) as f32)); + println!("\n nb element {:?} nb search : {:?} recall rate is {:?} search time inverse {:?} ", nb_elem, nbtest, mean_recall, 1.0e+6_f32/mean_search_time); + +} // end test1 + +#[test] +fn test_parallel() { + // + let _res = simple_logger::init(); + // + let nb_elem = 1000; + let dim = 25; + let knbn = 10; + let ef_c = 800; + let max_nb_connection = 48; + let ef = 20; + // + // + let mut data = gen_random_matrix_f32(dim, nb_elem); + for i in 0..data.len() { + l2_normalize(&mut data[i]); + } + let data_with_id = data.iter().zip(0..data.len()).collect(); + let nb_layer = 16.min((nb_elem as f32).ln().trunc() as usize); + let mut hns = Hnsw::::new(max_nb_connection, nb_elem, nb_layer, ef_c, dist::DistDot{}); + // ! + // hns.set_extend_candidates(true); + let mut start = ProcessTime::now(); + let now = std::time::SystemTime::now(); + // parallel insertion + hns.parallel_insert(&data_with_id); + let mut cpu_time: Duration = start.elapsed(); + println!("\n hnsw data parallel insertion cpu time {:?} , system time {:?}", cpu_time, now.elapsed()); + // one serial more to check + let mut v = gen_random_vector_f32(dim); + l2_normalize(&mut v); + hns.insert((&v, hns.get_nb_point() +1)); + // + hns.dump_layer_info(); + println!(" hnsw data nb point inserted {:?}", hns.get_nb_point()); + // + println!("\n hnsw testing requests ..."); + let nbtest=100; + let mut recalls = Vec::::with_capacity(nbtest); + let mut recalls_id = Vec::::with_capacity(nbtest); + + let mut search_times = Vec::::with_capacity(nbtest); + for _itest in 0..nbtest { + let mut r_vec = Vec::::with_capacity(dim); + let mut rng = thread_rng(); + let unif = Uniform::::new(0.,1.); + for _ in 0..dim { + r_vec.push(rng.sample(unif)); + } + l2_normalize(&mut r_vec); + + start = ProcessTime::now(); + let brute_neighbours = brute_force_neighbours(knbn, hns.get_point_indexation() , Box::new(dist::DistDot) , &r_vec); + cpu_time = start.elapsed(); + if nbtest <= 100 { + println!("\n\n test_par nb_elem {:?}", nb_elem); + println!("\n brute force neighbours :" ); + println!("======================"); + println!(" brute force computing {:?} \n", cpu_time); + for i in 0..brute_neighbours.len() { + println!(" {:?} {:?} ", brute_neighbours[i].point_id, brute_neighbours[i].dist_to_ref); + } + } + // + let knbn = 10; + hns.set_searching_mode(true); + start = ProcessTime::now(); + let knn_neighbours = hns.search(&r_vec, knbn, ef); + cpu_time = start.elapsed(); + search_times.push(cpu_time.as_micros() as f32); + if nbtest <= 100 { + println!("\n knn neighbours"); + println!("======================"); + println!(" hnsw searching {:?} \n", cpu_time); + for n in &knn_neighbours { + println!(" {:?} \t {:?} \t {:?}", n.d_id, n.p_id, n.distance); + } + } + // compute recall with balls + let knn_neighbours_dist : Vec = knn_neighbours.iter().map(|p| p.distance).collect(); + let max_dist = brute_neighbours[knbn-1].dist_to_ref; + let recall = knn_neighbours_dist.iter().filter(|d| *d <= &max_dist).count(); + if nbtest <= 100 { + println!("recall {:?}", (recall as f32)/ (knbn as f32)); + } + recalls.push(recall); + // compute recall with id + let mut recall_id = 0; + let mut knn_neighbours_id : Vec = knn_neighbours.iter().map(|p| p.p_id).collect(); + knn_neighbours_id.sort_unstable(); + let snbn = knbn.min(brute_neighbours.len()); + for j in 0..snbn { + let to_search = brute_neighbours[j].point_id; + if knn_neighbours_id.binary_search(&to_search).is_ok() { + recall_id +=1; + } + } + recalls_id.push(recall_id); + } // end on nbtest + // + // compute recall + // + + let mean_recall = (recalls.iter().sum::() as f32)/((knbn * recalls.len()) as f32); + let mean_search_time = (search_times.iter().sum::() as f32)/(search_times.len() as f32); + println!("\n nb search {:?} recall rate is {:?} search time inverse {:?} ", nbtest, mean_recall, 1.0e+6_f32/mean_search_time); + let mean_recall_id = (recalls.iter().sum::() as f32)/((knbn * recalls.len()) as f32); + println!("mean recall rate with point ids {:?}", mean_recall_id); + // + // assert!(1==0); +} // end test_par + + +} // end mod tests \ No newline at end of file