package prbnmcn-dagger

  1. Overview
  2. Docs

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 (* assert s1 = s2 *)
          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 (* assert s1 = s2 *)
          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) }

  (* This would work for any mappable container. *)
  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 _ ->
                                  (* A distribution can't dynamically switch from stateless to kernel
                                     (TODO: lift this) *)
                                  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
                                  (* Format.printf
                                   *   "reevaluating variable, fresh sample: %d@."
                                   *   uid ; *)
                                  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 _ ->
                                  (* A distribution can't dynamically switch from kernel to stateless
                                     (TODO: lift this) *)
                                  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
  (* Format.printf "evaluating@." ; *)
  let { Traced.value = first_value; trace = first_trace } = run v in
  (* Format.printf "initial trace: %a@." Trace.pp first_trace.samples ; *)
  let mcmc_move prev_value prev_trace _fwd_ll _bwd_ll undo =
    (* Format.printf "evaluating@." ; *)
    let { Traced.value = new_value; trace = new_trace } =
      Incremental_monad.run v
    in
    (* Format.printf
     *   "prev_trace: %a new trace: %a@."
     *   Trace.pp
     *   prev_trace.trace.samples
     *   Trace.pp
     *   new_trace.samples ; *)
    let new_trace = process new_trace in
    let intersection =
      Trace.intersect_samples prev_trace.trace.samples new_trace.trace.samples
    in
    (* Format.printf "old/new intersection: %a@." Trace.pp intersection ; *)
    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)
OCaml

Innovation. Community. Security.