Skip to main content

gen_plots/
gen_plots.rs

1//! Generate Plotly HTML snippets for MkDocs documentation.
2//!
3//! Replaces the separate plot_ode / plot_interp examples and the gen_plots.py
4//! Python script. Run with:
5//!
6//! ```sh
7//! cargo run --example gen_plots --features all
8//! ```
9//!
10//! Writes interactive Plotly HTML snippets to `docs/includes/`.
11
12use numeris::control::{butterworth_lowpass, lag_compensator, lead_compensator, BiquadCascade};
13use numeris::interp::{CubicSpline, HermiteInterp, LagrangeInterp, LinearInterp};
14use numeris::ode::{AdaptiveSettings, RKAdaptive, Rosenbrock, RKTS54, RODAS4};
15use numeris::stats::{
16    Beta, Binomial, ContinuousDistribution, DiscreteDistribution, Gamma, Normal, Poisson,
17};
18use numeris::Vector;
19
20use std::f64::consts::PI;
21use std::fs;
22use std::path::Path;
23
24// ─── helpers ──────────────────────────────────────────────────────────────
25
26fn fmt_arr(v: &[f64]) -> String {
27    let inner: Vec<String> = v.iter().map(|x| format!("{x:.6}")).collect();
28    format!("[{}]", inner.join(","))
29}
30
31/// Shared axis and font styling applied to every plot layout.
32///
33/// Takes a partial layout JSON object (must start with `{` and end with `}`)
34/// and injects the common decoration keys. The caller provides title, axis
35/// titles, and any plot-specific overrides; this function adds:
36///
37/// - Serif font family at 15 px for body, 17 px for axis titles, 18 px bold
38///   for the plot title
39/// - Axis lines + mirror + ticks on all four sides
40/// - Light gridlines with dash pattern
41/// - Consistent legend and margin defaults (caller can override)
42fn decorate_layout(
43    title: &str,
44    x_title: &str,
45    y_title: &str,
46    extra: &str, // additional JSON key-value pairs (with leading comma)
47) -> String {
48    // We build the JSON manually to avoid pulling in serde.
49    format!(
50        concat!(
51            "{{",
52            "\"title\":{{\"text\":\"{title}\",\"font\":{{\"family\":\"Georgia, 'Times New Roman', serif\",\"size\":18}}}},",
53            "\"font\":{{\"family\":\"Georgia, 'Times New Roman', serif\",\"size\":15}},",
54            "\"xaxis\":{{",
55            "\"title\":{{\"text\":\"{x_title}\",\"font\":{{\"size\":16}}}},",
56            "\"showline\":true,\"linewidth\":1,\"linecolor\":\"black\",",
57            "\"mirror\":true,",
58            "\"ticks\":\"outside\",\"ticklen\":5,\"tickwidth\":1,\"tickcolor\":\"black\",",
59            "\"showgrid\":true,\"gridwidth\":1,\"gridcolor\":\"rgba(180,180,180,0.35)\",\"griddash\":\"dot\"",
60            "{x_extra}",
61            "}},",
62            "\"yaxis\":{{",
63            "\"title\":{{\"text\":\"{y_title}\",\"font\":{{\"size\":16}}}},",
64            "\"showline\":true,\"linewidth\":1,\"linecolor\":\"black\",",
65            "\"mirror\":true,",
66            "\"ticks\":\"outside\",\"ticklen\":5,\"tickwidth\":1,\"tickcolor\":\"black\",",
67            "\"showgrid\":true,\"gridwidth\":1,\"gridcolor\":\"rgba(180,180,180,0.35)\",\"griddash\":\"dot\"",
68            "{y_extra}",
69            "}},",
70            "\"legend\":{{\"orientation\":\"h\",\"y\":-0.22,\"font\":{{\"size\":14}}}},",
71            "\"margin\":{{\"t\":55,\"b\":70,\"l\":70,\"r\":30}},",
72            "\"plot_bgcolor\":\"white\",\"paper_bgcolor\":\"white\"",
73            "{extra}",
74            "}}",
75        ),
76        title = title,
77        x_title = x_title,
78        y_title = y_title,
79        x_extra = "",
80        y_extra = "",
81        extra = extra,
82    )
83}
84
85/// Like `decorate_layout` but with per-axis extra JSON fragments.
86fn decorate_layout_ex(
87    title: &str,
88    x_title: &str,
89    x_extra: &str,
90    y_title: &str,
91    y_extra: &str,
92    extra: &str,
93) -> String {
94    format!(
95        concat!(
96            "{{",
97            "\"title\":{{\"text\":\"{title}\",\"font\":{{\"family\":\"Georgia, 'Times New Roman', serif\",\"size\":18}}}},",
98            "\"font\":{{\"family\":\"Georgia, 'Times New Roman', serif\",\"size\":15}},",
99            "\"xaxis\":{{",
100            "\"title\":{{\"text\":\"{x_title}\",\"font\":{{\"size\":16}}}},",
101            "\"showline\":true,\"linewidth\":1,\"linecolor\":\"black\",",
102            "\"mirror\":true,",
103            "\"ticks\":\"outside\",\"ticklen\":5,\"tickwidth\":1,\"tickcolor\":\"black\",",
104            "\"showgrid\":true,\"gridwidth\":1,\"gridcolor\":\"rgba(180,180,180,0.35)\",\"griddash\":\"dot\"",
105            "{x_extra}",
106            "}},",
107            "\"yaxis\":{{",
108            "\"title\":{{\"text\":\"{y_title}\",\"font\":{{\"size\":16}}}},",
109            "\"showline\":true,\"linewidth\":1,\"linecolor\":\"black\",",
110            "\"mirror\":true,",
111            "\"ticks\":\"outside\",\"ticklen\":5,\"tickwidth\":1,\"tickcolor\":\"black\",",
112            "\"showgrid\":true,\"gridwidth\":1,\"gridcolor\":\"rgba(180,180,180,0.35)\",\"griddash\":\"dot\"",
113            "{y_extra}",
114            "}},",
115            "\"legend\":{{\"orientation\":\"h\",\"y\":-0.22,\"font\":{{\"size\":14}}}},",
116            "\"margin\":{{\"t\":55,\"b\":70,\"l\":70,\"r\":30}},",
117            "\"plot_bgcolor\":\"white\",\"paper_bgcolor\":\"white\"",
118            "{extra}",
119            "}}",
120        ),
121        title = title,
122        x_title = x_title,
123        x_extra = x_extra,
124        y_title = y_title,
125        y_extra = y_extra,
126        extra = extra,
127    )
128}
129
130fn plotly_snippet(div_id: &str, traces_json: &str, layout_json: &str, height: u32) -> String {
131    format!(
132        "<div id=\"{div_id}\" style=\"width:100%;height:{height}px;\"></div>\n\
133         <script>\n\
134         window.addEventListener('load',function(){{Plotly.newPlot('{div_id}',{traces_json},{layout_json},{{responsive:true}})}});\n\
135         </script>\n"
136    )
137}
138
139fn write_snippet(dir: &Path, name: &str, html: &str) {
140    let path = dir.join(format!("{name}.html"));
141    fs::write(&path, html).unwrap_or_else(|e| panic!("failed to write {}: {e}", path.display()));
142    println!("  → {}", path.display());
143}
144
145// ─── ODE plot ─────────────────────────────────────────────────────────────
146
147fn make_ode_plot() -> String {
148    let tau = 4.0 * PI;
149    let y0 = Vector::from_array([1.0_f64, 0.0]);
150    let settings = AdaptiveSettings {
151        dense_output: true,
152        ..AdaptiveSettings::default()
153    };
154    let sol = RKTS54::integrate(
155        0.0,
156        tau,
157        &y0,
158        |_t, y| Vector::from_array([y[1], -y[0]]),
159        &settings,
160    )
161    .expect("ODE integration failed");
162
163    const N: usize = 300;
164    let mut t = vec![0.0; N];
165    let mut x = vec![0.0; N];
166    let mut v = vec![0.0; N];
167    for i in 0..N {
168        let ti = tau * i as f64 / (N - 1) as f64;
169        let yi = RKTS54::interpolate(ti, &sol).unwrap();
170        t[i] = ti;
171        x[i] = yi[0];
172        v[i] = yi[1];
173    }
174
175    let traces = format!(
176        "[{{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"position x(t)\",\
177          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5}}}},\
178         {{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"velocity v(t)\",\
179          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5,\"dash\":\"dash\"}}}}]",
180        fmt_arr(&t),
181        fmt_arr(&x),
182        fmt_arr(&t),
183        fmt_arr(&v),
184    );
185
186    let layout = decorate_layout(
187        "Harmonic Oscillator — RKTS54 Dense Output",
188        "t",
189        "Value",
190        "",
191    );
192
193    plotly_snippet("plot-ode", &traces, &layout, 420)
194}
195
196// ─── Interpolation plot ───────────────────────────────────────────────────
197
198fn make_interp_plot() -> String {
199    let tau = 2.0 * PI;
200    let kx: [f64; 6] = core::array::from_fn(|i| tau * i as f64 / 5.0);
201    let ky: [f64; 6] = core::array::from_fn(|i| kx[i].sin());
202    let kd: [f64; 6] = core::array::from_fn(|i| kx[i].cos());
203
204    let linear = LinearInterp::new(kx, ky).unwrap();
205    let hermite = HermiteInterp::new(kx, ky, kd).unwrap();
206    let lagrange = LagrangeInterp::new(kx, ky).unwrap();
207    let spline = CubicSpline::new(kx, ky).unwrap();
208
209    const N: usize = 200;
210    let mut xv = vec![0.0; N];
211    let mut y_true = vec![0.0; N];
212    let mut y_lin = vec![0.0; N];
213    let mut y_her = vec![0.0; N];
214    let mut y_lag = vec![0.0; N];
215    let mut y_spl = vec![0.0; N];
216    for i in 0..N {
217        let xi = tau * i as f64 / (N - 1) as f64;
218        xv[i] = xi;
219        y_true[i] = xi.sin();
220        y_lin[i] = linear.eval(xi);
221        y_her[i] = hermite.eval(xi);
222        y_lag[i] = lagrange.eval(xi);
223        y_spl[i] = spline.eval(xi);
224    }
225
226    let traces = format!(
227        "[{{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"sin(x) exact\",\
228          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5,\"color\":\"rgba(120,120,120,0.6)\",\"dash\":\"dot\"}}}},\
229         {{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"Linear\",\
230          \"x\":{},\"y\":{},\"line\":{{\"width\":2}}}},\
231         {{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"Hermite\",\
232          \"x\":{},\"y\":{},\"line\":{{\"width\":2}}}},\
233         {{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"Lagrange\",\
234          \"x\":{},\"y\":{},\"line\":{{\"width\":2}}}},\
235         {{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"Cubic Spline\",\
236          \"x\":{},\"y\":{},\"line\":{{\"width\":2}}}},\
237         {{\"type\":\"scatter\",\"mode\":\"markers\",\"name\":\"knots\",\
238          \"x\":{},\"y\":{},\"marker\":{{\"size\":9,\"color\":\"black\",\"symbol\":\"diamond\"}}}}]",
239        fmt_arr(&xv), fmt_arr(&y_true),
240        fmt_arr(&xv), fmt_arr(&y_lin),
241        fmt_arr(&xv), fmt_arr(&y_her),
242        fmt_arr(&xv), fmt_arr(&y_lag),
243        fmt_arr(&xv), fmt_arr(&y_spl),
244        fmt_arr(&kx), fmt_arr(&ky),
245    );
246
247    let layout = decorate_layout(
248        "Interpolation Methods on sin(x) — 6 Knots",
249        "x",
250        "y",
251        "",
252    );
253
254    plotly_snippet("plot-interp", &traces, &layout, 440)
255}
256
257// ─── Control: Butterworth frequency response ──────────────────────────────
258
259fn biquad_cascade_freq_response<const N: usize>(
260    cascade: &BiquadCascade<f64, N>,
261    freq: f64,
262    fs: f64,
263) -> f64 {
264    let omega = 2.0 * PI * freq / fs;
265    let (sin_w, cos_w) = omega.sin_cos();
266    let cos_2w = 2.0 * cos_w * cos_w - 1.0;
267    let sin_2w = 2.0 * sin_w * cos_w;
268
269    let mut mag_sq = 1.0;
270    for section in &cascade.sections {
271        let (b, a) = section.coefficients();
272        let nr = b[0] + b[1] * cos_w + b[2] * cos_2w;
273        let ni = -b[1] * sin_w - b[2] * sin_2w;
274        let dr = a[0] + a[1] * cos_w + a[2] * cos_2w;
275        let di = -a[1] * sin_w - a[2] * sin_2w;
276        mag_sq *= (nr * nr + ni * ni) / (dr * dr + di * di);
277    }
278    mag_sq.sqrt()
279}
280
281fn make_control_plot() -> String {
282    let fs = 8000.0;
283    let fc = 1000.0;
284
285    let bw2: BiquadCascade<f64, 1> = butterworth_lowpass(2, fc, fs).unwrap();
286    let bw4: BiquadCascade<f64, 2> = butterworth_lowpass(4, fc, fs).unwrap();
287    let bw6: BiquadCascade<f64, 3> = butterworth_lowpass(6, fc, fs).unwrap();
288
289    const N: usize = 500;
290    let mut freqs = vec![0.0; N];
291    let mut db2 = vec![0.0; N];
292    let mut db4 = vec![0.0; N];
293    let mut db6 = vec![0.0; N];
294
295    let f_min: f64 = 10.0;
296    let f_max: f64 = 3900.0;
297    for i in 0..N {
298        let f = f_min * (f_max / f_min).powf(i as f64 / (N - 1) as f64);
299        freqs[i] = f;
300        db2[i] = 20.0 * biquad_cascade_freq_response(&bw2, f, fs).log10();
301        db4[i] = 20.0 * biquad_cascade_freq_response(&bw4, f, fs).log10();
302        db6[i] = 20.0 * biquad_cascade_freq_response(&bw6, f, fs).log10();
303    }
304
305    let traces = format!(
306        "[{{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"2nd order\",\
307          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5}}}},\
308         {{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"4th order\",\
309          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5}}}},\
310         {{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"6th order\",\
311          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5}}}}]",
312        fmt_arr(&freqs), fmt_arr(&db2),
313        fmt_arr(&freqs), fmt_arr(&db4),
314        fmt_arr(&freqs), fmt_arr(&db6),
315    );
316
317    let layout = decorate_layout_ex(
318        "Butterworth Lowpass — f<sub>c</sub> = 1 kHz, f<sub>s</sub> = 8 kHz",
319        "Frequency (Hz)",
320        ",\"type\":\"log\"",
321        "Magnitude (dB)",
322        ",\"range\":[-80,5]",
323        &format!(
324            ",\"shapes\":[{{\"type\":\"line\",\"x0\":{f_min},\"x1\":{f_max},\
325             \"y0\":-3,\"y1\":-3,\"line\":{{\"dash\":\"dot\",\"color\":\"rgba(160,80,80,0.5)\",\"width\":1.5}}}}]"
326        ),
327    );
328
329    plotly_snippet("plot-control", &traces, &layout, 420)
330}
331
332// ─── Control: Lead/Lag compensator Bode ───────────────────────────────────
333
334fn biquad_freq_response(b: &[f64; 3], a: &[f64; 3], freq: f64, fs: f64) -> (f64, f64) {
335    let omega = 2.0 * PI * freq / fs;
336    let (sin_w, cos_w) = omega.sin_cos();
337    let cos_2w = 2.0 * cos_w * cos_w - 1.0;
338    let sin_2w = 2.0 * sin_w * cos_w;
339    let nr = b[0] + b[1] * cos_w + b[2] * cos_2w;
340    let ni = -b[1] * sin_w - b[2] * sin_2w;
341    let dr = a[0] + a[1] * cos_w + a[2] * cos_2w;
342    let di = -a[1] * sin_w - a[2] * sin_2w;
343    let mag = ((nr * nr + ni * ni) / (dr * dr + di * di)).sqrt();
344    let phase = (ni.atan2(nr) - di.atan2(dr)).to_degrees();
345    (mag, phase)
346}
347
348fn make_lead_lag_plot() -> String {
349    let fs = 1000.0;
350    let lead = lead_compensator(std::f64::consts::FRAC_PI_4, 50.0, 1.0, fs).unwrap();
351    let lag = lag_compensator(10.0, 5.0, fs).unwrap();
352
353    let (b_lead, a_lead) = lead.coefficients();
354    let (b_lag, a_lag) = lag.coefficients();
355
356    const N: usize = 400;
357    let f_min: f64 = 0.1;
358    let f_max: f64 = 490.0;
359    let mut freqs = vec![0.0; N];
360    let mut lead_db = vec![0.0; N];
361    let mut lead_ph = vec![0.0; N];
362    let mut lag_db = vec![0.0; N];
363    let mut lag_ph = vec![0.0; N];
364
365    for i in 0..N {
366        let f = f_min * (f_max / f_min).powf(i as f64 / (N - 1) as f64);
367        freqs[i] = f;
368        let (m, p) = biquad_freq_response(&b_lead, &a_lead, f, fs);
369        lead_db[i] = 20.0 * m.log10();
370        lead_ph[i] = p;
371        let (m, p) = biquad_freq_response(&b_lag, &a_lag, f, fs);
372        lag_db[i] = 20.0 * m.log10();
373        lag_ph[i] = p;
374    }
375
376    // Magnitude plot
377    let traces_mag = format!(
378        "[{{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"Lead (45° @ 50 Hz)\",\
379          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5}}}},\
380         {{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"Lag (10× DC @ 5 Hz)\",\
381          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5}}}}]",
382        fmt_arr(&freqs), fmt_arr(&lead_db),
383        fmt_arr(&freqs), fmt_arr(&lag_db),
384    );
385    let layout_mag = decorate_layout_ex(
386        "Lead / Lag Compensators — Magnitude",
387        "Frequency (Hz)",
388        ",\"type\":\"log\"",
389        "Magnitude (dB)",
390        "",
391        "",
392    );
393
394    // Phase plot
395    let traces_ph = format!(
396        "[{{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"Lead phase\",\
397          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5}}}},\
398         {{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"Lag phase\",\
399          \"x\":{},\"y\":{},\"line\":{{\"width\":2.5}}}}]",
400        fmt_arr(&freqs), fmt_arr(&lead_ph),
401        fmt_arr(&freqs), fmt_arr(&lag_ph),
402    );
403    let layout_ph = decorate_layout_ex(
404        "Lead / Lag Compensators — Phase",
405        "Frequency (Hz)",
406        ",\"type\":\"log\"",
407        "Phase (°)",
408        "",
409        "",
410    );
411
412    let mut html = plotly_snippet("plot-lead-lag-mag", &traces_mag, &layout_mag, 380);
413    html.push_str(&plotly_snippet(
414        "plot-lead-lag-phase",
415        &traces_ph,
416        &layout_ph,
417        380,
418    ));
419    html
420}
421
422// ─── ODE: Van der Pol (stiff, RODAS4) ─────────────────────────────────────
423
424fn make_vanderpol_plot() -> String {
425    let mu = 20.0_f64;
426    let y0 = Vector::from_array([2.0, 0.0]);
427    let t_end = 120.0;
428
429    let settings = AdaptiveSettings {
430        abs_tol: 1e-8,
431        rel_tol: 1e-8,
432        max_steps: 100_000,
433        dense_output: true,
434        ..AdaptiveSettings::default()
435    };
436
437    let sol = RODAS4::integrate(
438        0.0,
439        t_end,
440        &y0,
441        |_t, y| {
442            Vector::from_array([y[1], mu * (1.0 - y[0] * y[0]) * y[1] - y[0]])
443        },
444        |_t, y| {
445            numeris::Matrix::new([
446                [0.0, 1.0],
447                [-2.0 * mu * y[0] * y[1] - 1.0, mu * (1.0 - y[0] * y[0])],
448            ])
449        },
450        &settings,
451    )
452    .expect("Van der Pol integration failed");
453
454    // Downsample accepted step points to a fixed grid for a manageable HTML size.
455    // The adaptive solver clusters points at sharp transitions; we keep enough
456    // resolution by picking the nearest stored point for each output sample.
457    let ds = sol.dense.as_ref().expect("no dense output");
458    let n_out = 2000usize;
459    let mut tv = Vec::with_capacity(n_out);
460    let mut xv = Vec::with_capacity(n_out);
461    let mut idx = 0usize;
462    for i in 0..n_out {
463        let t_want = t_end * i as f64 / (n_out - 1) as f64;
464        // advance index to nearest stored point
465        while idx + 1 < ds.t.len() && ds.t[idx + 1] <= t_want {
466            idx += 1;
467        }
468        tv.push(t_want);
469        xv.push(ds.y[idx][0]);
470    }
471
472    let traces = format!(
473        "[{{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"y₁(t)\",\
474          \"x\":{},\"y\":{},\"line\":{{\"width\":2}}}}]",
475        fmt_arr(&tv),
476        fmt_arr(&xv),
477    );
478
479    let layout = decorate_layout(
480        "Van der Pol Oscillator (μ = 20) — RODAS4",
481        "t",
482        "y₁",
483        "",
484    );
485
486    plotly_snippet("plot-vanderpol", &traces, &layout, 420)
487}
488
489// ─── Stats: continuous PDF plots ──────────────────────────────────────────
490
491/// Helper: build a single line trace JSON fragment.
492fn line_trace(name: &str, x: &[f64], y: &[f64]) -> String {
493    format!(
494        "{{\"type\":\"scatter\",\"mode\":\"lines\",\"name\":\"{name}\",\
495         \"x\":{},\"y\":{},\"line\":{{\"width\":2.5}}}}",
496        fmt_arr(x),
497        fmt_arr(y),
498    )
499}
500
501/// Helper: build a bar trace JSON fragment (for PMFs).
502fn bar_trace(name: &str, x: &[f64], y: &[f64]) -> String {
503    format!(
504        "{{\"type\":\"bar\",\"name\":\"{name}\",\
505         \"x\":{},\"y\":{},\"opacity\":0.7}}",
506        fmt_arr(x),
507        fmt_arr(y),
508    )
509}
510
511fn make_normal_pdf_plot() -> String {
512    let dists: Vec<(Normal<f64>, &str)> = vec![
513        (Normal::new(0.0, 1.0).unwrap(), "μ=0, σ=1"),
514        (Normal::new(0.0, 2.0).unwrap(), "μ=0, σ=2"),
515        (Normal::new(2.0, 0.7).unwrap(), "μ=2, σ=0.7"),
516    ];
517
518    const N: usize = 300;
519    let x_min = -6.0_f64;
520    let x_max = 6.0_f64;
521    let mut xv: Vec<f64> = (0..N)
522        .map(|i| x_min + (x_max - x_min) * i as f64 / (N - 1) as f64)
523        .collect();
524    // ensure exact 0
525    xv[N / 2] = 0.0;
526
527    let traces: Vec<String> = dists
528        .iter()
529        .map(|(d, name)| {
530            let yv: Vec<f64> = xv.iter().map(|&x| d.pdf(x)).collect();
531            line_trace(name, &xv, &yv)
532        })
533        .collect();
534
535    let layout = decorate_layout("Normal Distribution — PDF", "x", "f(x)", "");
536    plotly_snippet("plot-normal-pdf", &format!("[{}]", traces.join(",")), &layout, 400)
537}
538
539fn make_gamma_pdf_plot() -> String {
540    let dists: Vec<(Gamma<f64>, &str)> = vec![
541        (Gamma::new(1.0, 1.0).unwrap(), "α=1, β=1 (Exp)"),
542        (Gamma::new(2.0, 1.0).unwrap(), "α=2, β=1"),
543        (Gamma::new(5.0, 1.0).unwrap(), "α=5, β=1"),
544        (Gamma::new(2.0, 2.0).unwrap(), "α=2, β=2"),
545    ];
546
547    const N: usize = 300;
548    let x_min = 0.01_f64;
549    let x_max = 12.0_f64;
550    let xv: Vec<f64> = (0..N)
551        .map(|i| x_min + (x_max - x_min) * i as f64 / (N - 1) as f64)
552        .collect();
553
554    let traces: Vec<String> = dists
555        .iter()
556        .map(|(d, name)| {
557            let yv: Vec<f64> = xv.iter().map(|&x| d.pdf(x)).collect();
558            line_trace(name, &xv, &yv)
559        })
560        .collect();
561
562    let layout = decorate_layout("Gamma Distribution — PDF", "x", "f(x)", "");
563    plotly_snippet("plot-gamma-pdf", &format!("[{}]", traces.join(",")), &layout, 400)
564}
565
566fn make_beta_pdf_plot() -> String {
567    let dists: Vec<(Beta<f64>, &str)> = vec![
568        (Beta::new(0.5, 0.5).unwrap(), "α=0.5, β=0.5"),
569        (Beta::new(2.0, 2.0).unwrap(), "α=2, β=2"),
570        (Beta::new(2.0, 5.0).unwrap(), "α=2, β=5"),
571        (Beta::new(5.0, 2.0).unwrap(), "α=5, β=2"),
572    ];
573
574    const N: usize = 300;
575    let eps = 0.005;
576    let xv: Vec<f64> = (0..N)
577        .map(|i| eps + (1.0 - 2.0 * eps) * i as f64 / (N - 1) as f64)
578        .collect();
579
580    let traces: Vec<String> = dists
581        .iter()
582        .map(|(d, name)| {
583            let yv: Vec<f64> = xv.iter().map(|&x| d.pdf(x).min(8.0)).collect();
584            line_trace(name, &xv, &yv)
585        })
586        .collect();
587
588    let layout = decorate_layout_ex(
589        "Beta Distribution — PDF",
590        "x",
591        "",
592        "f(x)",
593        ",\"range\":[0,4.5]",
594        "",
595    );
596    plotly_snippet("plot-beta-pdf", &format!("[{}]", traces.join(",")), &layout, 400)
597}
598
599// ─── Stats: discrete PMF plots ───────────────────────────────────────────
600
601fn make_binomial_pmf_plot() -> String {
602    let dists: Vec<(Binomial<f64>, &str)> = vec![
603        (Binomial::new(10, 0.3).unwrap(), "n=10, p=0.3"),
604        (Binomial::new(10, 0.5).unwrap(), "n=10, p=0.5"),
605        (Binomial::new(20, 0.7).unwrap(), "n=20, p=0.7"),
606    ];
607
608    let k_max = 21u64;
609    let kv: Vec<f64> = (0..=k_max).map(|k| k as f64).collect();
610
611    let traces: Vec<String> = dists
612        .iter()
613        .map(|(d, name)| {
614            let yv: Vec<f64> = (0..=k_max).map(|k| d.pmf(k)).collect();
615            bar_trace(name, &kv, &yv)
616        })
617        .collect();
618
619    let layout = decorate_layout_ex(
620        "Binomial Distribution — PMF",
621        "k",
622        "",
623        "P(X = k)",
624        "",
625        ",\"barmode\":\"group\"",
626    );
627    plotly_snippet("plot-binomial-pmf", &format!("[{}]", traces.join(",")), &layout, 400)
628}
629
630fn make_poisson_pmf_plot() -> String {
631    let dists: Vec<(Poisson<f64>, &str)> = vec![
632        (Poisson::new(1.0).unwrap(), "λ = 1"),
633        (Poisson::new(4.0).unwrap(), "λ = 4"),
634        (Poisson::new(10.0).unwrap(), "λ = 10"),
635    ];
636
637    let k_max = 20u64;
638    let kv: Vec<f64> = (0..=k_max).map(|k| k as f64).collect();
639
640    let traces: Vec<String> = dists
641        .iter()
642        .map(|(d, name)| {
643            let yv: Vec<f64> = (0..=k_max).map(|k| d.pmf(k)).collect();
644            bar_trace(name, &kv, &yv)
645        })
646        .collect();
647
648    let layout = decorate_layout_ex(
649        "Poisson Distribution — PMF",
650        "k",
651        "",
652        "P(X = k)",
653        "",
654        ",\"barmode\":\"group\"",
655    );
656    plotly_snippet("plot-poisson-pmf", &format!("[{}]", traces.join(",")), &layout, 400)
657}
658
659fn make_continuous_cdf_plot() -> String {
660    let normal = Normal::new(0.0, 1.0).unwrap();
661    let gamma = Gamma::new(3.0, 1.0).unwrap();
662    let beta = Beta::new(2.0, 5.0).unwrap();
663
664    const N: usize = 300;
665
666    // Normal CDF on [-4, 4]
667    let xn: Vec<f64> = (0..N)
668        .map(|i| -4.0 + 8.0 * i as f64 / (N - 1) as f64)
669        .collect();
670    let yn: Vec<f64> = xn.iter().map(|&x| normal.cdf(x)).collect();
671
672    // Gamma CDF on [0, 10]
673    let xg: Vec<f64> = (0..N)
674        .map(|i| 0.01 + 10.0 * i as f64 / (N - 1) as f64)
675        .collect();
676    let yg: Vec<f64> = xg.iter().map(|&x| gamma.cdf(x)).collect();
677
678    // Beta CDF on [0, 1]
679    let xb: Vec<f64> = (0..N)
680        .map(|i| 0.005 + 0.99 * i as f64 / (N - 1) as f64)
681        .collect();
682    let yb: Vec<f64> = xb.iter().map(|&x| beta.cdf(x)).collect();
683
684    let traces = format!(
685        "[{},{},{}]",
686        line_trace("Normal(0, 1)", &xn, &yn),
687        line_trace("Gamma(3, 1)", &xg, &yg),
688        line_trace("Beta(2, 5)", &xb, &yb),
689    );
690
691    let layout = decorate_layout(
692        "Continuous Distributions — CDF",
693        "x",
694        "F(x)",
695        "",
696    );
697    plotly_snippet("plot-continuous-cdf", &traces, &layout, 400)
698}
699
700// ─── main ─────────────────────────────────────────────────────────────────
701
702fn main() {
703    let includes = Path::new(env!("CARGO_MANIFEST_DIR")).join("docs/includes");
704    fs::create_dir_all(&includes).expect("failed to create docs/includes/");
705
706    println!("Generating Plotly HTML snippets...");
707
708    write_snippet(&includes, "plot_ode", &make_ode_plot());
709    write_snippet(&includes, "plot_vanderpol", &make_vanderpol_plot());
710    write_snippet(&includes, "plot_interp", &make_interp_plot());
711    write_snippet(&includes, "plot_control", &make_control_plot());
712    write_snippet(&includes, "plot_lead_lag", &make_lead_lag_plot());
713    write_snippet(&includes, "plot_normal_pdf", &make_normal_pdf_plot());
714    write_snippet(&includes, "plot_gamma_pdf", &make_gamma_pdf_plot());
715    write_snippet(&includes, "plot_beta_pdf", &make_beta_pdf_plot());
716    write_snippet(&includes, "plot_binomial_pmf", &make_binomial_pmf_plot());
717    write_snippet(&includes, "plot_poisson_pmf", &make_poisson_pmf_plot());
718    write_snippet(&includes, "plot_continuous_cdf", &make_continuous_cdf_plot());
719
720    println!("Done.");
721}