Source file bisec_tree.ml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
open Printf
module A = BatArray
module L = List
module type Point = sig
type t
val dist: t -> t -> float
end
type vp_heuristic = One_band | Two_bands
type direction = Left | Right
type step = L of float
| R of float
let string_of_addr addr =
let char_of_step = function
| L _ -> '0'
| R _ -> '1' in
let buff = Buffer.create 80 in
L.iter (fun a ->
Buffer.add_char buff (char_of_step a)
) addr;
Buffer.contents buff
let string_of_path path =
let char_of_dir = function
| Left -> '0'
| Right -> '1' in
let buff = Buffer.create 80 in
L.iter (fun d ->
Buffer.add_char buff (char_of_dir d)
) path;
Buffer.contents buff
module Make = functor (P: Point) -> struct
type bucket = { vp: P.t;
sup: float;
points: P.t array }
type node =
{
l_vp: P.t;
l_sup: float;
r_vp: P.t;
r_sup: float;
left: t;
right: t }
and t = Empty
| Node of node
| Bucket of bucket
let rng = Random.State.make_self_init ()
let rand_int n =
Random.State.int rng n
let fcmp (x: float) (y: float): int =
if x < y then -1
else if x > y then 1
else 0
let fmax (x: float) (y: float): float =
if x > y then x else y
type point1 = { p: P.t;
d1: float }
let point1_cmp (x: point1) (y: point1): int =
fcmp x.d1 y.d1
let enr (vp: P.t) (p: P.t): point1 =
{ p; d1 = P.dist vp p }
type point2 = { p: P.t;
d1: float;
d2: float }
let enr2 (vp: P.t) (p: point1): point2 =
{ p = p.p; d1 = p.d1; d2 = P.dist vp p.p }
let strip2 (points: point2 array): P.t array =
A.map (fun x -> x.p) points
let max1 (points: point2 array): float =
let maxi = ref 0.0 in
A.iter (fun x ->
maxi := fmax !maxi x.d1
) points;
!maxi
let max2 (points: point2 array): float =
let maxi = ref 0.0 in
A.iter (fun x ->
maxi := fmax !maxi x.d2
) points;
!maxi
type pre_bucket = { vp: P.t;
points: point2 array }
type pre_node = { l_vp: P.t;
points: point2 array;
r_vp: P.t }
type pre = Pre_bucket of pre_bucket
| Pre_node of pre_node
| Pre_empty
let pre_bucket_length (b: pre_bucket): int =
1 + A.length b.points
let bucket_length (b: bucket): int =
1 + A.length b.points
let rand_vp (points: P.t array): point1 array =
let n = A.length points in
assert(n > 0);
if n = 1 then [|{ p = points.(0); d1 = 0.0 }|]
else
let i = rand_int n in
let vp = points.(i) in
let enr_points = A.map (enr vp) points in
A.sort point1_cmp enr_points;
enr_points
let one_band (k: int) (points: P.t array) =
let n = A.length points in
if n = 0 then Pre_empty
else if n = 1 then Pre_bucket { vp = points.(0); points = [||] }
else
let enr_points = rand_vp points in
let vp1 = enr_points.(0).p in
let vp2 = enr_points.(n - 1).p in
let lr_gap = enr_points.(n - 1).d1 in
if n = 2 || n <= k || lr_gap = 0.0 then
let enr_rem = A.sub enr_points 0 (n - 1) in
let rem = A.map (enr2 vp2) enr_rem in
Pre_bucket { vp = vp2; points = rem }
else
let enr_rem = A.sub enr_points 1 (n - 2) in
let rem = A.map (enr2 vp2) enr_rem in
Pre_node { l_vp = vp1; points = rem; r_vp = vp2 }
let two_bands (k: int) (points: P.t array) =
let n = A.length points in
if n = 0 then Pre_empty
else if n = 1 then Pre_bucket { vp = points.(0); points = [||] }
else
let enr_points = rand_vp points in
let vp = enr_points.(n - 1).p in
let enr_points1 = A.map (enr vp) points in
A.sort point1_cmp enr_points1;
let vp1 = enr_points1.(0).p in
let vp2 = enr_points1.(n - 1).p in
let lr_gap = enr_points1.(n - 1).d1 in
if n = 2 || n <= k || lr_gap = 0.0 then
let enr_rem = A.sub enr_points1 0 (n - 1) in
let rem = A.map (enr2 vp2) enr_rem in
Pre_bucket { vp = vp2; points = rem }
else
let enr_rem = A.sub enr_points1 1 (n - 2) in
let rem = A.map (enr2 vp2) enr_rem in
Pre_node { l_vp = vp1; points = rem; r_vp = vp2 }
let sample_distances (sample_size: int) (points: P.t array): float array =
let n = A.length points in
assert(n > 0);
let distances =
A.init sample_size (fun _ ->
P.dist points.(rand_int n) points.(rand_int n)
) in
A.sort fcmp distances;
distances
let create ?(progress_callback = fun _x _y -> ())
(k: int) (h: vp_heuristic) (points': P.t array): t =
let nb_points = A.length points' in
let indexed = ref 0 in
let heuristic = match h with
| One_band -> one_band
| Two_bands -> two_bands in
let rec loop points = match heuristic k points with
| Pre_empty -> Empty
| Pre_bucket b ->
begin
indexed := !indexed + (pre_bucket_length b);
progress_callback !indexed nb_points;
Bucket { vp = b.vp; sup = max2 b.points; points = strip2 b.points }
end
| Pre_node pn ->
let lpoints, rpoints = A.partition (fun p -> p.d1 < p.d2) pn.points in
indexed := !indexed + 2;
progress_callback !indexed nb_points;
Node { l_vp = pn.l_vp;
l_sup = max1 lpoints;
r_vp = pn.r_vp;
r_sup = max2 rpoints;
left = loop (strip2 lpoints);
right = loop (strip2 rpoints) } in
loop points'
let rec to_list_loop acc = function
| Empty -> acc
| Node n ->
let acc' = to_list_loop acc n.right in
to_list_loop (n.l_vp :: n.r_vp :: acc') n.left
| Bucket b ->
A.fold_left (fun acc' x ->
x :: acc'
) (b.vp :: acc) b.points
let to_list t =
to_list_loop [] t
let length t =
let rec loop acc = function
| Empty -> acc
| Node n ->
let acc' = loop acc n.right in
loop (acc' + 2) n.left
| Bucket b ->
1 + acc + (A.length b.points) in
loop 0 t
let dump max_depth t =
let rec loop acc path curr_depth = function
| Empty -> acc
| Bucket b ->
let points = to_list (Bucket b) in
(L.rev path, points) :: acc
| Node n ->
if curr_depth = max_depth then
let l_points = n.l_vp :: to_list n.left in
let l_path = Left :: path in
let r_points = n.r_vp :: to_list n.right in
let r_path = Right :: path in
(L.rev l_path, l_points) ::
(L.rev r_path, r_points) :: acc
else
let depth' = curr_depth + 1 in
let l_path = Left :: path in
let r_path = Right :: path in
let acc' = (L.rev l_path, [n.l_vp]) :: acc in
let acc'' = loop acc' l_path depth' n.left in
let acc''' = (L.rev r_path, [n.r_vp]) :: acc'' in
loop acc''' r_path depth' n.right in
loop [] [] 1 t
let is_empty = function
| Empty -> true
| _ -> false
let root = function
| Empty -> raise Not_found
| Node n -> n.l_vp
| Bucket b -> b.vp
let nearest_neighbor query tree =
let rec loop ((_x, d) as acc) = function
| Empty -> acc
| Bucket b ->
let b_d = P.dist query b.vp in
let x', d' = if b_d < d then (b.vp, b_d) else acc in
if b_d -. b.sup >= d' then (x', d')
else
A.fold_left (fun (nearest_p, nearest_d) y ->
let y_d = P.dist query y in
if y_d < nearest_d then (y, y_d) else (nearest_p, nearest_d)
) (x', d') b.points
| Node n ->
let l_d = P.dist query n.l_vp in
let x', d' = if l_d < d then (n.l_vp, l_d) else acc in
let x'', d'' =
if l_d -. n.l_sup >= d' then (x', d')
else loop (x', d') n.left in
let r_d = P.dist query n.r_vp in
let x''', d''' = if r_d < d'' then (n.r_vp, r_d) else (x'', d'') in
if r_d -. n.r_sup >= d''' then (x''', d''')
else loop (x''', d''') n.right in
match tree with
| Empty -> raise Not_found
| not_empty ->
let x = root not_empty in
loop (x, P.dist query x) not_empty
let neighbors query tol tree =
let rec loop acc = function
| Empty -> acc
| Bucket b ->
let b_d = P.dist query b.vp in
if b_d -. b.sup > tol then acc
else if b_d +. b.sup <= tol then
to_list_loop acc (Bucket b)
else
let acc' = if b_d <= tol then b.vp :: acc else acc in
A.fold_left (fun acc'' y ->
let y_d = P.dist query y in
if y_d <= tol then y :: acc'' else acc''
) acc' b.points
| Node n ->
let l_d = P.dist query n.l_vp in
let acc'' =
if l_d -. n.l_sup > tol then acc
else if l_d +. n.l_sup <= tol then
to_list_loop (n.l_vp :: acc) n.left
else
let acc' = if l_d <= tol then n.l_vp :: acc else acc in
loop acc' n.left in
let r_d = P.dist query n.r_vp in
if r_d -. n.r_sup > tol then acc''
else if r_d +. n.r_sup <= tol then
to_list_loop (n.r_vp :: acc'') n.right
else
let acc''' = if r_d <= tol then (n.r_vp :: acc'') else acc'' in
loop acc''' n.right in
loop [] tree
let partition query tol tree =
let rec loop included excluded = function
| Empty -> (included, excluded)
| Bucket b ->
let b_d = P.dist query b.vp in
if b_d -. b.sup > tol then
(included, to_list_loop excluded (Bucket b))
else if b_d +. b.sup <= tol then
(to_list_loop included (Bucket b), excluded)
else
let included', excluded' =
if b_d <= tol then
(b.vp :: included, excluded)
else
(included, b.vp :: excluded) in
A.fold_left (fun (included'', excluded'') y ->
let y_d = P.dist query y in
if y_d <= tol then
(y :: included'', excluded'')
else
(included'', y :: excluded'')
) (included', excluded') b.points
| Node n ->
let l_d = P.dist query n.l_vp in
let included', excluded' =
if l_d -. n.l_sup > tol then
(included, to_list_loop (n.l_vp :: excluded) n.left)
else if l_d +. n.l_sup <= tol then
(to_list_loop (n.l_vp :: included) n.left, excluded)
else
let included'', excluded'' =
if l_d <= tol then
(n.l_vp :: included, excluded)
else
(included, n.l_vp :: excluded) in
loop included'' excluded'' n.left in
let r_d = P.dist query n.r_vp in
if r_d -. n.r_sup > tol then
(included', to_list_loop (n.r_vp :: excluded') n.right)
else if r_d +. n.r_sup <= tol then
(to_list_loop (n.r_vp :: included') n.right, excluded')
else
let included'', excluded'' =
if r_d <= tol then
(n.r_vp :: included', excluded')
else
(included', n.r_vp :: excluded') in
loop included'' excluded'' n.right in
loop [] [] tree
let rec check = function
| Empty -> true
| Bucket b ->
A.for_all (fun x ->
let d = P.dist b.vp x in
d <= b.sup
) b.points
| Node n ->
L.for_all (fun x ->
let l_d = P.dist n.l_vp x in
let r_d = P.dist n.r_vp x in
(l_d <= n.l_sup && l_d < r_d)
) (to_list n.left) &&
L.for_all (fun x ->
let r_d = P.dist n.r_vp x in
let l_d = P.dist n.l_vp x in
(r_d <= n.r_sup && r_d <= l_d)
) (to_list n.right) &&
check n.left && check n.right
let inspect tree =
let rec loop acc = function
| Empty -> acc
| Bucket b -> b.vp :: acc
| Node n ->
let acc' = n.l_vp :: n.r_vp :: acc in
let acc'' = loop acc' n.left in
loop acc'' n.right in
loop [] tree
let find query tree =
let nearest_p, nearest_d = nearest_neighbor query tree in
if nearest_d = 0.0 then nearest_p else raise Not_found
let mem query tree =
try let _ = find query tree in true
with Not_found -> false
let get_addr query tree =
let rec loop acc = function
| Empty | Bucket _ -> L.rev acc
| Node n ->
let l_d = P.dist query n.l_vp in
let r_d = P.dist query n.r_vp in
if l_d < r_d then
loop (L l_d :: acc) n.left
else
loop (R r_d :: acc) n.right in
loop [] tree
let add query address tree =
let rec loop addr = function
| Empty ->
begin match addr with
| [] -> Bucket { vp = query; sup = 0.0; points = [||] }
| _ -> assert(false)
end
| Bucket b ->
begin match addr with
| [] ->
let d = P.dist query b.vp in
let points = A.append b.points [|query|] in
Bucket { vp = b.vp; sup = max b.sup d; points }
| _ -> assert(false)
end
| Node n ->
begin match addr with
| [] -> assert(false)
| L l_d :: rest ->
Node { l_vp = n.l_vp;
l_sup = max n.l_sup l_d;
r_vp = n.r_vp;
r_sup = n.r_sup;
left = loop rest n.left;
right = n.right }
| R r_d :: rest ->
Node { l_vp = n.l_vp;
l_sup = n.l_sup;
r_vp = n.r_vp;
r_sup = max n.r_sup r_d;
left = n.left;
right = loop rest n.right }
end
in
loop address tree
let to_string t =
let rec loop path acc = function
| Empty ->
let str = sprintf "%s 0" (string_of_path (L.rev path)) in
str :: acc
| Bucket b ->
let str =
sprintf "%s %d"
(string_of_path (L.rev path))
(bucket_length b) in
str :: acc
| Node n ->
let acc' = loop (Left :: path) acc n.left in
loop (Right :: path) acc' n.right
in
let unsorted = loop [] [] t in
let sorted = List.sort compare unsorted in
String.concat "\n" sorted
end