Faiss

- Main Page
- Namespaces
- Classes
- Files

- File List

All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends
/data/users/matthijs/github_faiss/faiss/hamming.cpp
1
2 /**
3 * Copyright (c) 2015-present, Facebook, Inc.
4 * All rights reserved.
5 *
6 * This source code is licensed under the CC-by-NC license found in the
7 * LICENSE file in the root directory of this source tree.
8 */
9
10 /* Copyright 2004-present Facebook. All Rights Reserved.
11 *
12 * Implementation of Hamming related functions (distances, smallest distance
13 * selection with regular heap|radix and probabilistic heap|radix.
14 *
15 * IMPLEMENTATION NOTES
16 * Bitvectors are generally assumed to be multiples of 64 bits.
17 *
18 * hamdis_t is used for distances because at this time
19 * it is not clear how we will need to balance
20 * - flexibility in vector size (unclear more than 2^16 or even 2^8 bitvectors)
21 * - memory usage
22 * - cache-misses when dealing with large volumes of data (lower bits is better)
23 *
24 * The hamdis_t should optimally be compatibe with one of the Torch Storage
25 * (Byte,Short,Long) and therefore should be signed for 2-bytes and 4-bytes
26 */
27
28 #include "hamming.h"
29
30 #include <stdlib.h>
31 #include <stdio.h>
32 #include <math.h>
33 #include <assert.h>
34 #include <limits.h>
35
36 #include "Heap.h"
37 #include "FaissAssert.h"
38
39 static const size_t BLOCKSIZE_QUERY = 8192;
40
41
42 namespace faiss {
43
44 static const uint8_t hamdis_tab_ham_bytes[256] = {
45 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4,
46 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
47 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
48 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
49 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
50 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
51 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
52 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
53 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 5,
54 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
55 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
56 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
57 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 6,
58 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
59 3, 4, 4, 5, 4, 5, 5, 6, 4, 5, 5, 6, 5, 6, 6, 7,
60 4, 5, 5, 6, 5, 6, 6, 7, 5, 6, 6, 7, 6, 7, 7, 8
61 };
62
63
64 /* Elementary Hamming distance computation: unoptimized */
65 template <size_t nbits, typename T>
66 T hamming (const uint8_t *bs1,
67 const uint8_t *bs2)
68 {
69 const size_t nbytes = nbits / 8;
70 size_t i;
71 T h = 0;
72 for (i = 0; i < nbytes; i++)
73 h += (T) hamdis_tab_ham_bytes[bs1[i]^bs2[i]];
74 return h;
75 }
76
77
78 /* Hamming distances for multiples of 64 bits */
79 template <size_t nbits>
80 hamdis_t hamming (const uint64_t * bs1, const uint64_t * bs2)
81 {
82 const size_t nwords = nbits / 64;
83 size_t i;
84 hamdis_t h = 0;
85 for (i = 0; i < nwords; i++)
86 h += popcount64 (bs1[i] ^ bs2[i]);
87 return h;
88 }
89
90
91
92 /* specialized (optimized) functions */
93 template <>
94 hamdis_t hamming<64> (const uint64_t * pa, const uint64_t * pb)
95 {
96 return popcount64 (pa[0] ^ pb[0]);
97 }
98
99
100 template <>
101 hamdis_t hamming<128> (const uint64_t *pa, const uint64_t *pb)
102 {
103 return popcount64 (pa[0] ^ pb[0]) + popcount64(pa[1] ^ pb[1]);
104 }
105
106
107 template <>
108 hamdis_t hamming<256> (const uint64_t * pa, const uint64_t * pb)
109 {
110 return popcount64 (pa[0] ^ pb[0])
111 + popcount64 (pa[1] ^ pb[1])
112 + popcount64 (pa[2] ^ pb[2])
113 + popcount64 (pa[3] ^ pb[3]);
114 }
115
116
117 /* Hamming distances for multiple of 64 bits */
118 hamdis_t hamming (
119 const uint64_t * bs1,
120 const uint64_t * bs2,
121 size_t nwords)
122 {
123 size_t i;
124 hamdis_t h = 0;
125 for (i = 0; i < nwords; i++)
126 h += popcount64 (bs1[i] ^ bs2[i]);
127 return h;
128 }
129
130
131
132 template <size_t nbits>
133 void hammings (
134 const uint64_t * bs1,
135 const uint64_t * bs2,
136 size_t n1, size_t n2,
137 hamdis_t * dis)
138
139 {
140 size_t i, j;
141 const size_t nwords = nbits / 64;
142 for (i = 0; i < n1; i++) {
143 const uint64_t * __restrict bs1_ = bs1 + i * nwords;
144 hamdis_t * __restrict dis_ = dis + i * n2;
145 for (j = 0; j < n2; j++)
146 dis_[j] = hamming<nbits>(bs1_, bs2 + j * nwords);
147 }
148 }
149
150
151
152 void hammings (
153 const uint64_t * bs1,
154 const uint64_t * bs2,
155 size_t n1,
156 size_t n2,
157 size_t nwords,
158 hamdis_t * __restrict dis)
159 {
160 size_t i, j;
161 n1 *= nwords;
162 n2 *= nwords;
163 for (i = 0; i < n1; i+=nwords) {
164 const uint64_t * bs1_ = bs1+i;
165 for (j = 0; j < n2; j+=nwords)
166 dis[j] = hamming (bs1_, bs2+j, nwords);
167 }
168 }
169
170
171
172
173 /* Count number of matches given a max threshold */
174 template <size_t nbits>
175 void hamming_count_thres (
176 const uint64_t * bs1,
177 const uint64_t * bs2,
178 size_t n1,
179 size_t n2,
180 hamdis_t ht,
181 size_t * nptr)
182 {
183 const size_t nwords = nbits / 64;
184 size_t i, j, posm = 0;
185 const uint64_t * bs2_ = bs2;
186
187 for (i = 0; i < n1; i++) {
188 bs2 = bs2_;
189 for (j = 0; j < n2; j++) {
190 /* collect the match only if this satisfies the threshold */
191 if (hamming <nbits> (bs1, bs2) <= ht)
192 posm++;
193 bs2 += nwords;
194 }
195 bs1 += nwords; /* next signature */
196 }
197 *nptr = posm;
198 }
199
200
201 template <size_t nbits>
202 void crosshamming_count_thres (
203 const uint64_t * dbs,
204 size_t n,
205 int ht,
206 size_t * nptr)
207 {
208 const size_t nwords = nbits / 64;
209 size_t i, j, posm = 0;
210 const uint64_t * bs1 = dbs;
211 for (i = 0; i < n; i++) {
212 const uint64_t * bs2 = bs1 + 2;
213 for (j = i + 1; j < n; j++) {
214 /* collect the match only if this satisfies the threshold */
215 if (hamming <nbits> (bs1, bs2) <= ht)
216 posm++;
217 bs2 += nwords;
218 }
219 bs1 += nwords;
220 }
221 *nptr = posm;
222 }
223
224
225 template <size_t nbits>
226 size_t match_hamming_thres (
227 const uint64_t * bs1,
228 const uint64_t * bs2,
229 size_t n1,
230 size_t n2,
231 int ht,
232 long * idx,
233 hamdis_t * hams)
234 {
235 const size_t nwords = nbits / 64;
236 size_t i, j, posm = 0;
237 hamdis_t h;
238 const uint64_t * bs2_ = bs2;
239 for (i = 0; i < n1; i++) {
240 bs2 = bs2_;
241 for (j = 0; j < n2; j++) {
242 /* Here perform the real work of computing the distance */
243 h = hamming <nbits> (bs1, bs2);
244
245 /* collect the match only if this satisfies the threshold */
246 if (h <= ht) {
247 /* Enough space to store another match ? */
248 *idx = i; idx++;
249 *idx = j; idx++;
250 *hams = h;
251 hams++;
252 posm++;
253 }
254 bs2+=nwords; /* next signature */
255 }
256 bs1+=nwords;
257 }
258 return posm;
259 }
260
261
262 /* Return closest neighbors w.r.t Hamming distance */
263 template <class HammingComputer>
264 static
265 void hammings_knn_hc (
266 int bytes_per_code,
267 int_maxheap_array_t * ha,
268 const uint8_t * bs1,
269 const uint8_t * bs2,
270 size_t n2,
271 bool order = true,
272 bool init_heap = true)
273 {
274 size_t k = ha->k;
275
276
277 if (init_heap) ha->heapify ();
278
279 /* The computation here does not involved any blockization,
280 which is suboptimal for many queries in parallel. */
281 #pragma omp parallel for
282 for (size_t i = 0; i < ha->nh; i++) {
283 HammingComputer hc (bs1 + i * bytes_per_code, bytes_per_code);
284
285 const uint8_t * bs2_ = bs2;
286 hamdis_t dis;
287 hamdis_t * __restrict bh_val_ = ha->val + i * k;
288 long * __restrict bh_ids_ = ha->ids + i * k;
289 size_t j;
290 for (j = 0; j < n2; j++, bs2_+= bytes_per_code) {
291 dis = hc.hamming (bs2_);
292 if (dis < bh_val_[0]) {
293 faiss::maxheap_pop<hamdis_t> (k, bh_val_, bh_ids_);
294 faiss::maxheap_push<hamdis_t> (k, bh_val_, bh_ids_, dis, j);
295 }
296 }
297 }
298 if (order) ha->reorder ();
299 }
300
301
302
303 // works faster than the template version
304 static
305 void hammings_knn_1 (
306 int_maxheap_array_t * ha,
307 const uint64_t * bs1,
308 const uint64_t * bs2,
309 size_t n2,
310 bool order = true,
311 bool init_heap = true)
312 {
313 const size_t nwords = 1;
314 size_t k = ha->k;
315
316
317 if (init_heap) {
318 ha->heapify ();
319 }
320
321 #pragma omp parallel for
322 for (size_t i = 0; i < ha->nh; i++) {
323 const uint64_t bs1_ = bs1 [i];
324 const uint64_t * bs2_ = bs2;
325 hamdis_t dis;
326 hamdis_t * bh_val_ = ha->val + i * k;
327 hamdis_t bh_val_0 = bh_val_[0];
328 long * bh_ids_ = ha->ids + i * k;
329 size_t j;
330 for (j = 0; j < n2; j++, bs2_+= nwords) {
331 dis = popcount64 (bs1_ ^ *bs2_);
332 if (dis < bh_val_0) {
333 faiss::maxheap_pop<hamdis_t> (k, bh_val_, bh_ids_);
334 faiss::maxheap_push<hamdis_t> (k, bh_val_, bh_ids_, dis, j);
335 bh_val_0 = bh_val_[0];
336 }
337 }
338 }
339 if (order) {
340 ha->reorder ();
341 }
342 }
343
344
345
346
347 /* Functions to maps vectors to bits. Assume proper allocation done beforehand,
348 meaning that b should be be able to receive as many bits as x may produce. */
349
350 /*
351 * dimension 0 corresponds to the least significant bit of b[0], or
352 * equivalently to the lsb of the first byte that is stored.
353 */
354 void fvec2bitvec (const float * x, uint8_t * b, size_t d)
355 {
356 for (int i = 0; i < d; i += 8) {
357 uint8_t w = 0;
358 uint8_t mask = 1;
359 int nj = i + 8 <= d ? 8 : d - i;
360 for (int j = 0; j < nj; j++) {
361 if (x[i + j] >= 0)
362 w |= mask;
363 mask <<= 1;
364 }
365 *b = w;
366 b++;
367 }
368 }
369
370
371
372 /* Same but for n vectors.
373 Ensure that the ouptut b is byte-aligned (pad with 0s). */
374 void fvecs2bitvecs (const float * x, uint8_t * b, size_t d, size_t n)
375 {
376 const long ncodes = ((d + 7) / 8);
377 #pragma omp parallel for
378 for (size_t i = 0; i < n; i++)
379 fvec2bitvec (x + i * d, b + i * ncodes, d);
380 }
381
382
383 /* Reverse bit (NOT a optimized function, only used for print purpose) */
384 static uint64_t uint64_reverse_bits (uint64_t b)
385 {
386 int i;
387 uint64_t revb = 0;
388 for (i = 0; i < 64; i++) {
389 revb <<= 1;
390 revb |= b & 1;
391 b >>= 1;
392 }
393 return revb;
394 }
395
396
397 /* print the bit vector */
398 void bitvec_print (const uint8_t * b, size_t d)
399 {
400 size_t i, j;
401 for (i = 0; i < d; ) {
402 uint64_t brev = uint64_reverse_bits (* (uint64_t *) b);
403 for (j = 0; j < 64 && i < d; j++, i++) {
404 printf ("%d", (int) (brev & 1));
405 brev >>= 1;
406 }
407 b += 8;
408 printf (" ");
409 }
410 }
411
412
413
414
415
416 /*----------------------------------------*/
417 /* Hamming distance computation and k-nn */
418
419
420 #define C64(x) ((uint64_t *)x)
421
422
423 /* Compute a set of Hamming distances */
424 void hammings (
425 const uint8_t * a,
426 const uint8_t * b,
427 size_t na, size_t nb,
428 size_t ncodes,
429 hamdis_t * __restrict dis)
430 {
431 FAISS_ASSERT (ncodes % 8 == 0);
432 switch (ncodes) {
433 case 8:
434 faiss::hammings <64> (C64(a), C64(b), na, nb, dis); return;
435 case 16:
436 faiss::hammings <128> (C64(a), C64(b), na, nb, dis); return;
437 case 32:
438 faiss::hammings <256> (C64(a), C64(b), na, nb, dis); return;
439 case 64:
440 faiss::hammings <512> (C64(a), C64(b), na, nb, dis); return;
441 default:
442 faiss::hammings (C64(a), C64(b), na, nb, ncodes * 8, dis); return;
443 }
444 }
445
446
447 void hammings_knn_core (
448 int_maxheap_array_t * ha,
449 const uint8_t * a,
450 const uint8_t * b,
451 size_t nb,
452 size_t ncodes)
453 {
454 FAISS_ASSERT (ncodes % 8 == 0);
455
456 switch (ncodes) {
457 hammings_knn_1 (ha, C64(a), C64(b), nb, false, true);
458 // hammings_knn_hc<faiss::HammingComputer8> (8, ha, a, b, nb, false, true);
459 break;
460 case 16:
461 hammings_knn_hc<faiss::HammingComputer16> (16, ha, a, b, nb, false, true);
462 break;
463 case 32:
464 hammings_knn_hc<faiss::HammingComputer32> (32, ha, a, b, nb, false, true);
465 break;
466 default:
467 hammings_knn_hc<faiss::HammingComputerM8> (ncodes, ha, a, b, nb, false, true);
468 }
469 }
470
471 void hammings_knn (
472 int_maxheap_array_t * ha,
473 const uint8_t * a,
474 const uint8_t * b,
475 size_t nb,
476 size_t ncodes,
477 int order)
478 {
479 switch (ncodes) {
480 case 4:
481 hammings_knn_hc<faiss::HammingComputer4> (4, ha, a, b, nb, order, true);
482 break;
483 case 8:
484 hammings_knn_1 (ha, C64(a), C64(b), nb, order, true);
485 // hammings_knn_hc<faiss::HammingComputer8> (8, ha, a, b, nb, order, true);
486 break;
487 case 16:
488 hammings_knn_hc<faiss::HammingComputer16> (16, ha, a, b, nb, order, true);
489 break;
490 case 32:
491 hammings_knn_hc<faiss::HammingComputer32> (32, ha, a, b, nb, order, true);
492 break;
493 default:
494 FAISS_ASSERT (ncodes % 8 == 0);
495 hammings_knn_hc<faiss::HammingComputerM8> (ncodes, ha, a, b, nb, order, true);
496
497 }
498 }
499
500
501
502
503 /* Count number of matches given a max threshold */
504 void hamming_count_thres (
505 const uint8_t * bs1,
506 const uint8_t * bs2,
507 size_t n1,
508 size_t n2,
509 hamdis_t ht,
510 size_t ncodes,
511 size_t * nptr)
512 {
513 switch (ncodes) {
514 case 8:
515 faiss::hamming_count_thres <64> (C64(bs1), C64(bs2),
516 n1, n2, ht, nptr);
517 return;
518 case 16:
519 faiss::hamming_count_thres <128> (C64(bs1), C64(bs2),
520 n1, n2, ht, nptr);
521 return;
522 case 32:
523 faiss::hamming_count_thres <256> (C64(bs1), C64(bs2),
524 n1, n2, ht, nptr);
525 return;
526 case 64:
527 faiss::hamming_count_thres <512> (C64(bs1), C64(bs2),
528 n1, n2, ht, nptr);
529 return;
530 default:
531 FAISS_ASSERT (!"not-implemented for this number of bits");
532 }
533 }
534
535
536 /* Count number of cross-matches given a threshold */
537 void crosshamming_count_thres (
538 const uint8_t * dbs,
539 size_t n,
540 hamdis_t ht,
541 size_t ncodes,
542 size_t * nptr)
543 {
544 switch (ncodes) {
545 case 8:
546 faiss::crosshamming_count_thres <64> (C64(dbs), n, ht, nptr);
547 return;
548 case 16:
549 faiss::crosshamming_count_thres <128> (C64(dbs), n, ht, nptr);
550 return;
551 case 32:
552 faiss::crosshamming_count_thres <256> (C64(dbs), n, ht, nptr);
553 return;
554 case 64:
555 faiss::crosshamming_count_thres <512> (C64(dbs), n, ht, nptr);
556 return;
557 default:
558 FAISS_ASSERT (!"not-implemented for this number of bits");
559 }
560 }
561
562
563 /* Returns all matches given a threshold */
564 size_t match_hamming_thres (
565 const uint8_t * bs1,
566 const uint8_t * bs2,
567 size_t n1,
568 size_t n2,
569 hamdis_t ht,
570 size_t ncodes,
571 long * idx,
572 hamdis_t * dis)
573 {
574 switch (ncodes) {
575 case 8:
576 return faiss::match_hamming_thres <64> (C64(bs1), C64(bs2),
577 n1, n2, ht, idx, dis);
578 case 16:
579 return faiss::match_hamming_thres <128> (C64(bs1), C64(bs2),
580 n1, n2, ht, idx, dis);
581 case 32:
582 return faiss::match_hamming_thres <256> (C64(bs1), C64(bs2),
583 n1, n2, ht, idx, dis);
584 case 64:
585 return faiss::match_hamming_thres <512> (C64(bs1), C64(bs2),
586 n1, n2, ht, idx, dis);
587 default:
588 FAISS_ASSERT (!"not-implemented for this number of bits");
589 }
590 }
591
592
593 #undef C64
594
595
596
597 /*************************************
598 * generalized Hamming distances
599 ************************************/
600
601
602
603 template <class HammingComputer>
604 static void hamming_dis_inner_loop (
605 const uint8_t *ca,
606 const uint8_t *cb,
607 size_t nb,
608 size_t code_size,
609 int k,
610 hamdis_t * bh_val_,
611 long * bh_ids_)
612 {
613
614 HammingComputer hc (ca, code_size);
615
616 for (size_t j = 0; j < nb; j++) {
617 int ndiff = hc.hamming (cb);
618 cb += code_size;
619 if (ndiff < bh_val_[0]) {
620 maxheap_pop<hamdis_t> (k, bh_val_, bh_ids_);
621 maxheap_push<hamdis_t> (k, bh_val_, bh_ids_, ndiff, j);
622 }
623 }
624 }
625
626 void generalized_hammings_knn (
627 int_maxheap_array_t * ha,
628 const uint8_t * a,
629 const uint8_t * b,
630 size_t nb,
631 size_t code_size,
632 int ordered)
633 {
634 int na = ha->nh;
635 int k = ha->k;
636
637 if (ordered)
638 ha->heapify ();
639
640 #pragma omp parallel for
641 for (int i = 0; i < na; i++) {
642 const uint8_t *ca = a + i * code_size;
643 const uint8_t *cb = b;
644
645 hamdis_t * bh_val_ = ha->val + i * k;
646 long * bh_ids_ = ha->ids + i * k;
647
648 switch (code_size) {
649 case 8:
650 hamming_dis_inner_loop<GenHammingComputer8>
651 (ca, cb, nb, 8, k, bh_val_, bh_ids_);
652 break;
653 case 16:
654 hamming_dis_inner_loop<GenHammingComputer16>
655 (ca, cb, nb, 16, k, bh_val_, bh_ids_);
656 break;
657 case 32:
658 hamming_dis_inner_loop<GenHammingComputer32>
659 (ca, cb, nb, 32, k, bh_val_, bh_ids_);
660 break;
661 default:
662 hamming_dis_inner_loop<GenHammingComputerM8>
663 (ca, cb, nb, code_size, k, bh_val_, bh_ids_);
664 break;
665 }
666 }
667
668 if (ordered)
669 ha->reorder ();
670
671 }
672
673
674
675
676
677 } // namespace faiss
faiss::HeapArray::k
size_t k
allocated size per heap
Definition: Heap.h:356
faiss::generalized_hammings_knn
void generalized_hammings_knn(int_maxheap_array_t *ha, const uint8_t *a, const uint8_t *b, size_t nb, size_t code_size, int ordered)
Definition: hamming.cpp:626
faiss::HeapArray::reorder
void reorder()
reorder all the heaps
Definition: Heap.cpp:35
faiss::HeapArray
Definition: Heap.h:351
faiss::hammings_knn
void hammings_knn(int_maxheap_array_t *ha, const uint8_t *a, const uint8_t *b, size_t nb, size_t ncodes, int order)
Definition: hamming.cpp:471
faiss::HeapArray::ids
TI * ids
identifiers (size nh * k)
Definition: Heap.h:357
faiss::HeapArray::heapify
void heapify()
prepare all the heaps before adding
Definition: Heap.cpp:27
faiss::hammings
void hammings(const uint8_t *a, const uint8_t *b, size_t na, size_t nb, size_t nbytespercode, hamdis_t *dis)
faiss::HeapArray::val
T * val
values (distances or similarities), size nh * k
Definition: Heap.h:358
faiss::HeapArray::nh
size_t nh
number of heaps
Definition: Heap.h:355