sim 0.1.0

Sim provides a discrete event simulation engine, to facilitate Rust- and npm-based simulation products and projects
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
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

use std::f64::INFINITY;

use serde::{Deserialize, Serialize};
use wasm_bindgen::prelude::*;

pub mod t_scores;
use super::utils;

#[wasm_bindgen]
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ConfidenceInterval {
    lower: f64,
    upper: f64,
}

#[wasm_bindgen]
impl ConfidenceInterval {
    pub fn lower(&self) -> f64 {
        self.lower
    }

    pub fn upper(&self) -> f64 {
        self.upper
    }

    pub fn half_width(&self) -> f64 {
        (self.upper - self.lower) / 2.0
    }
}

#[wasm_bindgen]
#[derive(Debug, Default, Clone, Serialize, Deserialize)]
pub struct IndependentSample {
    // Independent observations only - use time series for auto-correlated data
    points: Vec<f64>,
    mean: f64,
    variance: f64,
}

#[wasm_bindgen]
impl IndependentSample {
    pub fn post(points: Vec<f64>) -> IndependentSample {
        let mean = utils::sample_mean(&points);
        let variance = utils::sample_variance(&points, &mean);
        IndependentSample {
            points,
            mean,
            variance,
        }
    }

    pub fn confidence_interval_mean(&self, alpha: f64) -> ConfidenceInterval {
        if self.points.len() == 1 {
            return ConfidenceInterval {
                lower: self.mean,
                upper: self.mean,
            };
        }
        ConfidenceInterval {
            lower: self.mean
                - t_scores::t_score(alpha, self.points.len() - 1) * self.variance.sqrt()
                    / (self.points.len() as f64).sqrt(),
            upper: self.mean
                + t_scores::t_score(alpha, self.points.len() - 1) * self.variance.sqrt()
                    / (self.points.len() as f64).sqrt(),
        }
    }

    pub fn point_estimate_mean(&self) -> f64 {
        self.mean
    }

    pub fn variance(&self) -> f64 {
        self.variance
    }
}

#[wasm_bindgen]
#[derive(Debug, Default, Clone, Serialize, Deserialize)]
pub struct TerminatingSimulationOutput {
    time_series_replications: Vec<Vec<f64>>,
    replication_means: Vec<f64>,
    replications_mean: Option<f64>,
    replications_variance: Option<f64>,
}

#[wasm_bindgen]
impl TerminatingSimulationOutput {
    pub fn new() -> TerminatingSimulationOutput {
        TerminatingSimulationOutput {
            ..Default::default()
        }
    }

    pub fn post_time_series(&mut self, time_series: Vec<f64>) {
        self.time_series_replications.push(time_series);
    }
}

#[wasm_bindgen]
#[derive(Debug, Default, Clone, Serialize, Deserialize)]
pub struct SteadyStateSimulationOutput {
    time_series: Vec<f64>,
    // Delete points at the beginning of the sample for initialization bias reduction
    deletion_point: Option<usize>,
    // Divide remining points into independent-ish batches
    batch_size: Option<usize>,
    batch_count: Option<usize>,
    batch_means: Vec<f64>,
    batches_mean: Option<f64>,
    batches_variance: Option<f64>,
}

#[wasm_bindgen]
impl SteadyStateSimulationOutput {
    pub fn new() -> SteadyStateSimulationOutput {
        SteadyStateSimulationOutput {
            ..Default::default()
        }
    }

    pub fn post(time_series: Vec<f64>) -> SteadyStateSimulationOutput {
        SteadyStateSimulationOutput {
            time_series,
            ..Default::default()
        }
    }

    pub fn put(&mut self, time_series: Vec<f64>) {
        self.time_series = time_series;
    }

    pub fn set_to_fixed_budget(&mut self) {
        let mut s = 0.0;
        let mut q = 0.0;
        let mut d = self.time_series.len() - 2;
        let mut mser = vec![0.0; self.time_series.len() - 1];
        loop {
            s += self.time_series[d + 1];
            q += self.time_series[d + 1].powi(2);
            mser[d] = (q - s.powi(2) / ((self.time_series.len() - d) as f64))
                / ((self.time_series.len() - d).pow(2) as f64);
            if d == 0 {
                // Find the minimum MSER in the first half of the time series
                let min_mser = (0..(self.time_series.len() - 1) / 2)
                    .fold(INFINITY, |min_mser, mser_index| {
                        f64::min(min_mser, mser[mser_index])
                    });
                // Use that point for deletion determination
                self.deletion_point = mser
                    .iter()
                    .position(|mser_value| utils::equivalent_f64(*mser_value, min_mser));
                break;
            } else {
                d -= 1;
            }
        }
        // Schmeiser [1982] found that, for a fixed total sample size, there is little benefit from dividing it into
        // more than k = 30 batches, even if we could do so and still retain independence between the batch means.
        self.batch_count = Some(f64::min(
            ((self.time_series.len() - self.deletion_point.unwrap()) as f64).sqrt(),
            30.0,
        ) as usize);
        self.batch_size = Some(
            (self.time_series.len() - self.deletion_point.unwrap()) / self.batch_count.unwrap(),
        );
        // if data are left over, eliminate from the beginning
        self.deletion_point =
            Some(self.time_series.len() - self.batch_count.unwrap() * self.batch_size.unwrap());
    }

    pub fn calculate_batch_statistics(&mut self) {
        if self.batch_count.is_none() {
            self.set_to_fixed_budget();
        }
        self.batch_means = (0..self.batch_count.unwrap())
            .map(|batch_index| {
                let batch_start_index =
                    self.deletion_point.unwrap() + self.batch_size.unwrap() * batch_index;
                let batch_end_index =
                    self.deletion_point.unwrap() + self.batch_size.unwrap() * (batch_index + 1);
                let points: Vec<f64> = (batch_start_index..batch_end_index)
                    .map(|index| self.time_series[index])
                    .collect();
                utils::sample_mean(&points)
            })
            .collect();
        self.batches_mean = Some(utils::sample_mean(&self.batch_means));
        self.batches_variance = Some(utils::sample_variance(
            &self.batch_means,
            &self.batches_mean.unwrap(),
        ));
    }

    pub fn confidence_interval_mean(&mut self, alpha: f64) -> ConfidenceInterval {
        if self.batches_mean.is_none() {
            self.calculate_batch_statistics();
        }
        if self.batch_count.unwrap() == 1 {
            return ConfidenceInterval {
                lower: self.batches_mean.unwrap(),
                upper: self.batches_mean.unwrap(),
            };
        }
        ConfidenceInterval {
            lower: self.batches_mean.unwrap()
                - t_scores::t_score(alpha, self.batch_count.unwrap() - 1)
                    * self.batches_variance.unwrap().sqrt()
                    / (self.batch_count.unwrap() as f64).sqrt(),
            upper: self.batches_mean.unwrap()
                + t_scores::t_score(alpha, self.batch_count.unwrap() - 1)
                    * self.batches_variance.unwrap().sqrt()
                    / (self.batch_count.unwrap() as f64).sqrt(),
        }
    }

    pub fn point_estimate_mean(&mut self) -> f64 {
        if self.batches_mean.is_none() {
            self.calculate_batch_statistics();
        }
        self.batches_mean.unwrap()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn epsilon() -> f64 {
        1.0e-12
    }

    #[test]
    fn confidence_interval_mean() {
        let sample = IndependentSample::post(vec![
            1.02, 0.73, 3.20, 0.23, 1.76, 0.47, 1.89, 1.45, 0.44, 0.23,
        ]);
        let confidence_interval = sample.confidence_interval_mean(0.1);
        assert!((confidence_interval.lower - 0.7492630635369267).abs() < epsilon());
        assert!((confidence_interval.upper - 1.534736936463073).abs() < epsilon());
    }
}