(* Airplane-blip example from the BLOG papers of Brian Milch et al. This example was mentioned in Stuart Russell's talk at the NIPS2008 Probabilistic Programming workshop. See also example 1.3 of http://people.csail.mit.edu/milch/papers/blog-chapter.pdf This example is meant to show-case open worlds and BLOG's ability to reason with open worlds. The blip example is HMM with a twist: the number of airplanes is not fixed; rather, it is described by a geometric distribution. Any number of airplanes is possible (although with progressively decreasing probabilities). We write our model explicitly as a HMM with a variable number of objects, evolving stochastically. We assert things like blip() as the evidence. So, we can answer AILog queries like observe blip(3,3,0) & blip(3,4,next(0)) as computing the probability of evidence on these two blips. But we can answer more queries: the distribution on the number of airplanes that are consistent with the observation; prediction on the trajectories of selected planes; probability that two blips are for the same plane, etc. In this file, we do not use laziness. This works for only one, max two planes. We also interpret the evidence about blips as partial: if we are told of the observed blip at position (3,3) we assume nothing is known about other positions on the radar screen. There may be blips at the other positions, too -- but we are not told about them. *) open Ptypes open ProbM type dir = North | East | South | West;; (* --------------- Priors *) (* We use the parameters from the AILog2 implementation of the example. See also blip_ailog.ml. *) let number_aircrafts () = geometric 0.85;; let npos = 10;; (* all coordinates are within 0..npos-1 *) (* Initial (x,y) positions of the plane i *) let xpos0 i = uniform npos;; let ypos0 i = uniform npos;; (* Initial direction of the plane i *) let dir0 i = uniformly [|North; East; South; West|];; (* x- or y-Derivatives for the plane i at time t travelling in the direction dir *) let xderiv i t = function | North -> dist [(0.1, -1); (0.8, 0); (0.1, 1)] | East -> dist [(0.2, 0); (0.7, 1); (0.1, 2)] | South -> dist [(0.1, -1); (0.8, 0); (0.1, 1)] | West -> dist [(0.2, 0); (0.7, -1); (0.1, 2)] ;; let yderiv i t = function | North -> dist [(0.2, 0); (0.7, 1); (0.1, 2)] | East -> dist [(0.1, -1); (0.8, 0); (0.1, 1)] | South -> dist [(0.2, 0); (0.7, -1); (0.1, 2)] | West -> dist [(0.1, -1); (0.8, 0); (0.1, 1)] ;; (* The new direction for the plane i travelling in the direction dir at the time t *) let new_dir i t = function | North -> dist [(0.2, West); (0.6, North); (0.2, East)] | East -> dist [(0.2, North); (0.6, East); (0.2, South)] | South -> dist [(0.2, East); (0.6, South); (0.2, West)] | West -> dist [(0.2, South); (0.6, West); (0.2, North)] ;; type pos = int * int;; (* --------------- The equations of motion and radar detection. *) (* The state of each plane at a given time moment *) type plane_state = {plane_idx : int; plane_pos : pos; plane_dir : dir};; (* obtain a sample initial plane state, for a plane number i *) let plane_state_0 i = {plane_idx = i; plane_pos = (xpos0 i, ypos0 i); plane_dir = dir0 i};; (* Evolve the state of one plane to the next time moment *) let plane_fly t pstate : plane_state = let i = pstate.plane_idx in let x' = xderiv i t pstate.plane_dir in let y' = yderiv i t pstate.plane_dir in let (x,y) = pstate.plane_pos in let newx = x + x' in let () = if newx < 0 || newx > 9 then fail () in let newy = y + y' in let () = if newy < 0 || newy > 9 then fail () in {plane_idx = i; plane_pos = (newx,newy); plane_dir = new_dir i t pstate.plane_dir} ;; type plane_states = plane_state list;; let planes_fly t pstates : plane_states = List.map (plane_fly t) pstates;; (* Asserting the evidence *) (* First, the simplest case (which was implemented in AILog) We are given the set of blips that have been observed at given position. Nothing is said about the blips at other positions. That is, the given set of blips is not the exhaustive set; there may have been other blips. *) let blips_at_least (blips : pos list) (pstates : plane_states) = let check_blip blip = flip 0.02 || (* randomly occurring blip *) List.fold_left (fun acc pstate -> acc || (pstate.plane_pos = blip && flip 0.9)) false pstates in List.iter (fun blip -> if check_blip blip then () else fail ()) blips ;; (* The ideal case with no noise *) let blips_at_least_ideal (blips : pos list) (pstates : plane_states) = let check_blip blip = List.fold_left (fun acc pstate -> acc || pstate.plane_pos = blip) false pstates in List.iter (fun blip -> if check_blip blip then () else fail ()) blips ;; (* Predict the observation at the current time moment: the blip given the states of the planes. *) (* Parameters of the model: nplanes -- Some n: the maximal expected number of planes None : the possible number of planes is unlimited ntimes -- the number of time steps to run the model (the last time step) observations: a function that checks the evidence. It receives the time moment and plane_states outcomes: the function that computes the final outcome *) let blip_HMM_model nplanes ntimes evidence foutcome () = let np = match nplanes with | Some n -> geometric_bounded n 0.85 | None -> geometric 0.85 in let rec state0 = function 0 -> [] | i -> plane_state_0 i :: state0 (pred i) in let rec iter t = let new_ps = if t = 0 then state0 np else let ps = variable_elim iter (pred t) in planes_fly t ps in let () = evidence new_ps t in (* check the evidence *) new_ps in foutcome (iter ntimes) ;; (* Check the prob of evidence of two blips, at pos (3,3) and (3,4) at successive time moments *) let two_blips ps = function 0 -> blips_at_least [(3,3)] ps | 1 -> blips_at_least [(3,4)] ps | _ -> () ;; (* a few tests first *) let [(0.0100000000000000019, V ())] = exact_reify (fun () -> if (plane_state_0 1).plane_pos = (3,3) then () else fail ());; let [(0.130434782608695676, V 1); (0.869565217391304324, V 0)] = exact_reify (fun () -> geometric_bounded 1 0.85);; (* explanation: 0.02 +. 0.98 *. 0.00130434782608695684 *. 0.9;; *) let [(0.0211504347826093099, V ())] = exact_reify (blip_HMM_model (Some 1) 0 (fun ps -> function 0 -> blips_at_least [(3,3)] ps) (fun _ -> ())) ;; (* the result is the product of prob of one plane times prob of that plane be at pos (3,3) *) let [(0.00130434782608695684, V ())] = exact_reify (blip_HMM_model (Some 1) 0 (fun ps -> function 0 -> blips_at_least_ideal [(3,3)] ps) (fun _ -> ())) ;; (* Two planes at the same time *) let [(0.0191897654584221797, V 2); (0.127931769722814503, V 1); (0.852878464818763282, V 0)] = exact_reify (fun () -> geometric_bounded 2 0.85);; let [(0.000199999999999999901, V ())] = exact_reify (let evidence = blips_at_least_ideal [(3,3); (3,5)] in fun () -> evidence [plane_state_0 1; plane_state_0 2]);; (* Indeed: plane 1 may be at (3,3) and plane 2 at (3,5) -- or the other way around. *) let [] = exact_reify (blip_HMM_model (Some 1) 0 (fun ps -> function 0 -> blips_at_least_ideal [(3,3); (3,5)] ps) (fun _ -> ())) ;; let [(3.83795309168443282e-06, V 2)] = exact_reify (blip_HMM_model (Some 2) 0 (let ev = blips_at_least_ideal [(3,3); (3,5)] in fun ps -> function 0 -> ev ps) (fun ps -> List.length ps)) ;; (* The exact result: prob of two planes is 0.0191897654584221797 The computed result is correct *) (* Alas, already the following takes too long... exact_reify (blip_HMM_model (Some 3) 0 (let ev = blips_at_least_ideal [(3,3); (3,5); (3,7)] in fun ps -> function 0 -> ev ps) (fun ps -> List.length ps)) ;; *) (* Planes in motion *) let [(0.00059257469565211177, V ())] = exact_reify (blip_HMM_model (Some 1) 1 two_blips (fun _ -> ())) ;; let [(0.000626031812173919, V ())] = sample_importance (random_selector 17) 1000 (blip_HMM_model (Some 1) 1 two_blips (fun _ -> ())) ;; (* But already with two planes the model becomes too complex to solve exactly. We remove the variable elimination. It seems sampling 5000 worlds gives a good enough approximation, and it finishes within a second. *) let blip_HMM_model nplanes ntimes evidence foutcome () = let np = match nplanes with | Some n -> geometric_bounded n 0.85 | None -> geometric 0.85 in let rec state0 = function 0 -> [] | i -> plane_state_0 i :: state0 (pred i) in let rec iter t = let new_ps = if t = 0 then state0 np else let ps = iter (pred t) in planes_fly t ps in let () = evidence new_ps t in (* check the evidence *) new_ps in foutcome (iter ntimes) ;; let [(0.000537223652173907085, V ())] = sample_importance (random_selector 17) 5000 (blip_HMM_model (Some 1) 1 two_blips (fun _ -> ())) ;; let [(0.000858385704719978691, V ())] = sample_importance (random_selector 17) 5000 (blip_HMM_model None 1 two_blips (fun _ -> ())) ;; (* This is the result we obtained in blip_ailog. But now we can show the distribution in the number of planes. *) let [(2.53840000000000075e-08, V 5); (7.1733848e-07, V 4); (2.89022758400000202e-05, V 3); (0.000126835698400000349, V 2); (0.000357885007999987398, V 1); (0.000344019999999996154, V 0)] = sample_importance (random_selector 17) 5000 (blip_HMM_model None 1 two_blips (fun ps -> List.length ps)) ;; (* The exact result for the query below is 3.83795309168443282e-06 As we increase the number of samples, the estimates converge to the exacty result. *) let [(5.88486140724946778e-06, V 2)] = sample_importance (random_selector 17) 5000 (blip_HMM_model (Some 2) 0 (let ev = blips_at_least_ideal [(3,3); (3,5)] in fun ps -> function 0 -> ev ps) (fun ps -> List.length ps)) ;; (* Alas, even with 10,000 samples, we don't get any result for three planes... We really should use laziness. *) let [] = sample_importance (random_selector 17) 5000 (blip_HMM_model (Some 3) 0 (let ev = blips_at_least_ideal [(3,3); (3,5); (3,7)] in fun ps -> function 0 -> ev ps) (fun ps -> List.length ps)) ;; (* However, if noise is present, we get an estimate ... *) let [(6.67493667269055322e-07, V 3); (9.12817797810127662e-07, V 2); (2.37663148208755515e-06, V 1); (4.19746996917189332e-06, V 0)] = sample_importance (random_selector 17) 5000 (blip_HMM_model (Some 3) 0 (let ev = blips_at_least [(3,3); (3,5); (3,7)] in fun ps -> function 0 -> ev ps) (fun ps -> List.length ps)) ;; (* The exact result is: 1e-6 *. (* join prob of three planes being at three cells of the grid *) 6.0 *. (* 3! *) 0.00287020304028914744;;(* prob that there are three planes, geom_bounded 3*) = 1.72212182417348842e-08 With 25000 samples, the estimated prob of three planes is 4.04041263569676812e-07. So we see a slow convergence to the exact result. Without noise, 25000 samples produce nothing. *) (* Let's try the example from Stuart Russell's talk: at each time moment, three marks are observed. *) let blip333 ps = function | 0 -> blips_at_least [(3,5); (3,7); (3,9)] ps | 1 -> blips_at_least [(4,5); (4,7); (4,9)] ps | 2 -> blips_at_least [(5,5); (5,7); (5,9)] ps | _ -> () ;; let (7.95996296929161396e-14, [(1.40953355226456531e-07, V 5); (0.0158591404344528564, V 4); (0.901408459486833902, V 3); (0.0250208630427655628, V 2); (0.0519967805529113342, V 1); (0.00571461552968104743, V 0)]) = Inference.normalize (sample_importance (random_selector 17) 5000 (blip_HMM_model None 2 blip333 (fun ps -> List.length ps))) ;; (* Three planes seem to be most likely *) (* Future work: modify HMM and evidence checking to return the index of a plane which produced the observed blip. The model should accumulate these facts and return the distribution: essentially the probability of different identifications of the plane. *)