modppl 0.3.1

a experimental library for probabilistic programming in Rust.
Documentation 
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

use modppl::real::consts::PI;
use modppl::{inference::ParticleSystem, normal, u01, Distribution, DynTrie, Real};
use nalgebra::dvector;
use rand::rngs::ThreadRng;
use std::fs::{create_dir_all, write};
use std::sync::Arc;

pub mod pointed_model;
use pointed_model::types_2d::{Bounds, Point};

pub mod dyngenfns;
use dyngenfns::spiral_model;

fn simulate_loop(rng: &mut ThreadRng, bounds: &Bounds, timesteps: i64) -> Vec<DynTrie> {
    let init_angle = u01(rng) * 2. * PI;

    let xrange = (bounds.xmax - bounds.xmin) as Real;
    let yrange = (bounds.ymax - bounds.ymin) as Real;
    let center = dvector![xrange / 2. + bounds.xmin, yrange / 2. + bounds.ymin];
    let radius = Real::max(bounds.xmax - bounds.xmin, bounds.ymax - bounds.ymin) / 5.;

    let mut observations = vec![];
    let perturb_means = (0..timesteps)
        .filter(|_| u01(rng) < 0.3)
        .collect::<Vec<_>>();
    for t in 0..timesteps {
        let mut deformation = 0.;
        for perturb_t in &perturb_means {
            deformation += normal
                .logpdf(&(t as Real), (perturb_t.clone() as Real, 1.))
                .exp()
        }
        let r = radius + deformation;
        let t = 2. * PI * (t as Real) / timesteps as Real;
        let obs = dvector![
            center[0] + r * (t + init_angle).cos(),
            center[1] + r * (t + init_angle).sin()
        ];
        let mut constraints = DynTrie::new();
        constraints.observe("obs", Arc::new(obs));
        observations.push(constraints);
    }
    observations
}

#[test]
fn test_smc() -> std::io::Result<()> {
    create_dir_all("../data")?;

    let mut rng = ThreadRng::default();
    const NUM_TIMESTEPS: i64 = 20;
    const NUM_PARTICLES: usize = 500;

    let bounds = Bounds {
        xmin: -1.,
        xmax: 1.,
        ymin: -1.,
        ymax: 1.,
    };
    let data = simulate_loop(&mut rng, &bounds, NUM_TIMESTEPS);
    let obs = data
        .iter()
        .map(|t| t.read::<Point>("obs"))
        .collect::<Vec<Point>>();
    let obs_strs = obs
        .iter()
        .map(|obs| format!("[{},{}]", obs[0], obs[1]))
        .collect::<Vec<String>>();
    let json = "[".to_owned() + &obs_strs.join(", ") + "]";
    write(format!("../data/smc_obs.json"), json)?;

    let mut filter = ParticleSystem::new(spiral_model, NUM_PARTICLES, rng);
    let mut data_it = data.into_iter();
    filter.init_step(dvector![0., 0.], vec![data_it.next().unwrap()]);

    let states = filter
        .traces
        .iter()
        .map(|vtr| vtr.retv.as_ref().unwrap().last().unwrap().clone())
        .collect::<Vec<_>>();
    let state_strs = states
        .iter()
        .map(|latent| format!("[{},{}]", latent[0], latent[1]))
        .collect::<Vec<String>>();
    let json = "[".to_owned() + &state_strs.join(", ") + "]";
    write(format!("../data/smc_traces_before_resample_0.json"), json)?;
    filter.resample();

    let states = filter
        .traces
        .iter()
        .map(|vtr| vtr.retv.as_ref().unwrap().last().unwrap().clone())
        .collect::<Vec<_>>();
    let state_strs = states
        .iter()
        .map(|latent| format!("[{},{}]", latent[0], latent[1]))
        .collect::<Vec<String>>();
    let json = "[".to_owned() + &state_strs.join(", ") + "]";
    write(format!("../data/smc_traces_{}.json", 0), json)?;

    for (t, constraints) in data_it.enumerate() {
        filter = filter.step(vec![constraints]);
        let states = filter
            .traces
            .iter()
            .map(|vtr| vtr.retv.as_ref().unwrap().last().unwrap().clone())
            .collect::<Vec<_>>();
        let state_strs = states
            .iter()
            .map(|latent| format!("[{},{}]", latent[0], latent[1]))
            .collect::<Vec<String>>();
        let json = "[".to_owned() + &state_strs.join(", ") + "]";
        write(
            format!("../data/smc_traces_before_resample_{}.json", t + 1),
            json,
        )?;
        filter.resample();
        let states = filter
            .traces
            .iter()
            .map(|vtr| vtr.retv.as_ref().unwrap().last().unwrap().clone())
            .collect::<Vec<_>>();
        let state_strs = states
            .iter()
            .map(|latent| format!("[{},{}]", latent[0], latent[1]))
            .collect::<Vec<String>>();
        let json = "[".to_owned() + &state_strs.join(", ") + "]";
        write(format!("../data/smc_traces_{}.json", t + 1), json)?;
    }

    Ok(())
}