Source file lmh_incremental_inference.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
module Trace = struct
type sample =
| Sample :
{ uid : int;
dist : 'a Dist.dist;
var : 'a Cgraph.Var.t;
score : Log_space.t
}
-> sample
| Kernel_sample :
{ uid : int;
dist : 'a Dist.kernel;
var : ('a * 'a) Cgraph.Var.t;
score : Log_space.t
}
-> sample
type sample_trace = sample list
type score = Score of { uid : int; score : Log_space.t }
type score_trace = score list
type t = { samples : sample_trace; scores : score_trace }
let empty = { samples = []; scores = [] }
let uid = function Sample { uid; _ } | Kernel_sample { uid; _ } -> uid
[@@inline]
let score = function
| Sample { score; _ } | Kernel_sample { score; _ } -> score
[@@inline]
module Internal_for_tests = struct
let rec equal_trace (trace1 : sample_trace) (trace2 : sample_trace) =
match (trace1, trace2) with
| ([], []) -> true
| ([], _) | (_, []) -> false
| (s1 :: tl1, s2 :: tl2) -> uid s1 = uid s2 && equal_trace tl1 tl2
let pp fmtr trace =
let open Format in
pp_print_list
~pp_sep:(fun fmtr () -> fprintf fmtr ", ")
(fun fmtr s ->
let uid = uid s in
fprintf fmtr "{%d}" uid)
fmtr
trace
end
[@@ocaml.warning "-32"]
let total_sample trace =
let rec loop list acc =
match list with
| [] -> acc
| hd :: tl -> loop tl (Log_space.mul (score hd) acc)
in
loop trace Log_space.one
let total_score trace =
let rec loop list acc =
match list with
| [] -> acc
| Score { score; _ } :: tl -> loop tl (Log_space.mul score acc)
in
loop trace Log_space.one
let total trace =
let total_sampling_score = total_sample trace.samples in
( total_sampling_score,
Log_space.mul total_sampling_score (total_score trace.scores) )
let cardinal { samples; _ } = List.length samples
let rec add_sample s trace =
match trace with
| [] -> [s]
| (Kernel_sample { uid = uid'; _ } as hd) :: tl
| (Sample { uid = uid'; _ } as hd) :: tl ->
let uid = uid s in
if uid < uid' then hd :: add_sample s tl
else if uid > uid' then s :: trace
else trace
let rec add_score (Score { uid; _ } as s) trace =
match trace with
| [] -> [s]
| (Score { uid = uid'; _ } as hd) :: tl ->
if uid < uid' then hd :: add_score s tl
else if uid > uid' then s :: trace
else trace
let add_sample (s : sample) trace =
{ trace with samples = add_sample s trace.samples }
let add_score (s : score) trace =
{ trace with scores = add_score s trace.scores }
let rec intersect_samples trace1 trace2 =
match (trace1, trace2) with
| ([], _) | (_, []) -> []
| (s1 :: tl1, s2 :: tl2) ->
let uid1 = uid s1 in
let uid2 = uid s2 in
if uid1 < uid2 then intersect_samples trace1 tl2
else if uid1 > uid2 then intersect_samples tl1 trace2
else s1 :: intersect_samples tl1 tl2
[@@ocaml.warning "-32"]
let rec union_samples trace1 trace2 =
match (trace1, trace2) with
| ([], t) | (t, []) -> t
| (s1 :: tl1, s2 :: tl2) ->
let uid1 = uid s1 in
let uid2 = uid s2 in
if uid1 < uid2 then s2 :: union_samples trace1 tl2
else if uid1 > uid2 then s1 :: union_samples tl1 trace2
else
s1 :: union_samples tl1 tl2
let rec union_scores trace1 trace2 =
match (trace1, trace2) with
| ([], t) | (t, []) -> t
| ( (Score { uid = uid1; _ } as s1) :: tl1,
(Score { uid = uid2; _ } as s2) :: tl2 ) ->
if uid1 < uid2 then s2 :: union_scores trace1 tl2
else if uid1 > uid2 then s1 :: union_scores tl1 trace2
else
s1 :: union_scores tl1 tl2
[@@ocaml.warning "-32"]
let union t1 t2 =
{ samples = union_samples t1.samples t2.samples;
scores = union_scores t1.scores t2.scores
}
end
module Counter = struct
let x = ref 0
let gen () =
let v = !x in
incr x ;
v
end
module Traced = Traced_monad.Make (Incremental_monad) (Trace)
module Syntax = struct
include Cps_monad.Make (Traced)
type 'a shared = 'a Traced.t
let with_shared (m : 'a t) (f : 'a shared -> 'b t) : 'b t =
fun ~handler ->
let m = (m ~handler).cont Fun.id in
{ cont = (fun k -> (f m ~handler).cont k) }
let with_shared_list (ms : 'a t list) (f : 'a shared list -> 'b t) : 'b t =
fun ~handler ->
let ms = List.map (fun m -> (m ~handler).cont Fun.id) ms in
{ cont = (fun k -> (f ms ~handler).cont k) }
let with_shared_array (ms : 'a t array) (f : 'a shared array -> 'b t) : 'b t =
fun ~handler ->
let ms = Array.map (fun m -> (m ~handler).cont Fun.id) ms in
{ cont = (fun k -> (f ms ~handler).cont k) }
let use : 'a shared -> 'a t =
fun node ~handler:_ -> { cont = (fun k -> k node) }
module Make_shared (C : sig
type 'a t
val map : 'a t -> ('a -> 'b) -> 'b t
end) =
struct
let with_shared (ms : 'a t C.t) (f : 'a shared C.t -> 'b t) : 'b t =
fun ~handler ->
let ms = C.map ms (fun m -> (m ~handler).cont Fun.id) in
{ cont = (fun k -> (f ms ~handler).cont k) }
end
module Infix = struct
include Infix
let ( let*! ) = with_shared
let use = use
end
end
let handler : RNG.t -> Syntax.handler =
let open Syntax in
fun rng_state ->
{ handler =
(fun (type a) (dist : a payload) ->
{ cont =
(fun k ->
match dist with
| Dist dist -> (
match (Cgraph.get dist).value with
| Stateless { sample; ll = _ } ->
let pos = sample rng_state in
let var = Cgraph.Var.create pos in
let node = Cgraph.var var in
let node =
Cgraph.map2 node dist (fun sample dist ->
match dist.value with
| Kernel _ ->
failwith
"Lmh_incremental_inference.handler: \
distribution switched from Stateless to \
Kernel"
| Stateless ({ ll; _ } as d) ->
let score = ll sample in
let uid = Counter.gen () in
let trace =
Trace.add_sample
(Trace.Sample
{ uid; dist = d; var; score })
dist.trace
in
{ Traced.value = sample; trace })
in
k node
| Kernel ({ start; sample; ll = _ } as d) ->
let pos = sample start rng_state in
let var = Cgraph.Var.create (start, pos) in
let node = Cgraph.var var in
let node =
Cgraph.map2 node dist (fun (prev, current) dist ->
match dist.value with
| Stateless _ ->
failwith
"Lmh_incremental_inference.handler: \
distribution switched from Kernel to \
Stateless"
| Kernel { ll; _ } ->
let score = ll prev current in
let uid = Counter.gen () in
let trace =
Trace.add_sample
(Trace.Kernel_sample
{ uid; dist = d; var; score })
dist.trace
in
{ Traced.value = current; trace })
in
k node)
| Score (m, ll) ->
k
@@ Cgraph.map2
m
ll
(fun
{ Traced.value; trace }
{ Traced.value = score; trace = trace' }
->
let uid = Counter.gen () in
{ Traced.value;
trace =
Trace.union
trace'
(Trace.add_score
(Trace.Score { uid; score })
trace)
}))
})
}
type processed_trace =
{ trace : Trace.t;
card : int;
samples : Trace.sample array Lazy.t;
sampling_score : Log_space.t;
score : Log_space.t
}
let to_dot fname (model : 'a Syntax.t) =
let oc = open_out fname in
let handler = handler (RNG.make [| 0x1337; 0x533D |]) in
Cgraph.Internal.set_debug true ;
let model = (model ~handler).cont Fun.id in
let _ = Cgraph.get model in
Cgraph.Internal.set_debug false ;
Cgraph.Internal.(to_dot ~mode:Full (Cgraph.ex (Obj.magic model)) oc) ;
close_out oc
let process : Trace.t -> processed_trace =
fun trace ->
let samples = lazy (Array.of_list trace.samples) in
let card = Trace.cardinal trace in
let (sampling_score, score) = Trace.total trace in
{ trace; card; samples; sampling_score; score }
let stream_samples (type a) (v : a Syntax.t) rng_state : a Seq.t =
let handler = handler rng_state in
let v = (v ~handler).cont Fun.id in
let select_resampling ({ samples; card; _ } : processed_trace) rng_state =
if card = 0 then None
else
let samples = Lazy.force samples in
Some samples.(RNG.int rng_state card)
[@@inline]
in
let run (v : a Traced.t) = Cgraph.get v in
let { Traced.value = first_value; trace = first_trace } = run v in
let mcmc_move prev_value prev_trace _fwd_ll _bwd_ll undo =
let { Traced.value = new_value; trace = new_trace } =
Incremental_monad.run v
in
let new_trace = process new_trace in
let intersection =
Trace.intersect_samples prev_trace.trace.samples new_trace.trace.samples
in
let intersection_score = Trace.total_sample intersection in
let forward_sampling_score =
Log_space.mul
Log_space.one
(Log_space.div new_trace.sampling_score intersection_score)
in
let backward_sampling_score =
Log_space.mul
Log_space.one
(Log_space.div prev_trace.sampling_score intersection_score)
in
let forward_flow =
Log_space.(
mul
prev_trace.score
(mul
(of_float (1. /. float_of_int prev_trace.card))
forward_sampling_score))
in
let backward_flow =
Log_space.(
mul
new_trace.score
(mul
(of_float (1. /. float_of_int new_trace.card))
backward_sampling_score))
in
let ratio = Log_space.div backward_flow forward_flow in
let acceptance = Log_space.(to_float (min one ratio)) in
if RNG.float rng_state 1.0 < acceptance then (new_value, new_trace)
else
let () = Cgraph.undo undo in
(prev_value, prev_trace)
in
let sample_step (prev_value : a) (prev_trace : processed_trace) rng_state =
match select_resampling prev_trace rng_state with
| None -> (prev_value, prev_trace)
| Some (Trace.Kernel_sample { uid = _; dist; var; score = _ }) ->
let (_previous, current) = Cgraph.Var.peek var in
let sample = dist.sample current rng_state in
let fwd_ll = dist.ll current sample in
let undo = Cgraph.Var.set_with_undo var (current, sample) in
let bwd_ll = dist.ll sample current in
mcmc_move prev_value prev_trace fwd_ll bwd_ll undo
| Some (Trace.Sample { uid = _; dist; var; score = bwd_ll }) ->
let (undo, fwd_ll) =
let sample = dist.sample rng_state in
let ll = dist.ll sample in
let undo = Cgraph.Var.set_with_undo var sample in
(undo, ll)
in
mcmc_move prev_value prev_trace fwd_ll bwd_ll undo
in
Seq.unfold
(fun (prev_value, prev_trace) ->
let next = sample_step prev_value prev_trace rng_state in
Some (prev_value, next))
(first_value, process first_trace)
include Syntax
module List_ops = Foldable.Make_list (Syntax)
module Seq_ops = Foldable.Make_seq (Syntax)
module Array_ops = Foldable.Make_array (Syntax)