sim/
sim.rs

1use aphreco::prelude::*;
2
3fn main() {
4  let model = Model::new();
5  // let stepper = Stepper::Rk4(StepOptions::Default);
6  let step_options = StepOptions::Dopri45 {
7    h0: 1e-3,
8    abstol: 1e-6,
9    reltol: 1e-6,
10    hmin: 1e-6,
11    hmax: 1e-3,
12  };
13  let stepper = Stepper::Dopri45(step_options);
14
15  let simulator = Simulator::new(model, stepper);
16
17  let sampling_time = smptime();
18  clock!(let simres = simulator.run(&sampling_time));
19  simres.save("./");
20}
21
22const LEN_Y: usize = 4;
23const LEN_P: usize = 11;
24const LEN_B: usize = 2;
25
26#[allow(dead_code)]
27pub struct Model {
28  p: [f64; LEN_P],
29}
30
31#[allow(dead_code)]
32impl SimModelTrait<LEN_Y, LEN_P, LEN_B> for Model {
33  fn new() -> Self {
34    let p = [
35      0.2,      // p[0] k12
36      0.05,     // p[1] k21
37      0.0,      // p[2] ini_b1
38      100000.0, // p[3] end_b1
39      0.1,      // p[4] tau_b1
40      1.0,      // p[5] ini_b2
41      100000.0, // p[6] end_b2
42      0.2,      // p[7] tau_b2
43      2.0,      // p[8] R_cre
44      10.0,     // p[9] X_dose_A
45      500.0,    // p[10] MW_A
46    ];
47
48    Self { p }
49  }
50
51  fn init(&self) -> (f64, [f64; LEN_Y]) {
52    let t0 = 0.0;
53    let y0 = [
54      100.0, // y[0] A
55      0.0,   // y[1] B
56      10.0,  // y[2] C
57      0.0,   // y[3] D
58    ];
59    let _b0 = ();
60    let _a0 = ();
61    (t0, y0)
62  }
63
64  #[allow(unused_variables)]
65  fn ode(&self, t: &f64, y: &[f64; LEN_Y], deriv_y: &mut [f64; LEN_Y]) {
66    deriv_y[0] = -self.p[0] * y[0] + self.p[1] * y[1];
67    deriv_y[1] = self.p[0] * y[0] - self.p[1] * y[1];
68  }
69
70  #[allow(unused_variables)]
71  fn rec(&self, t: &f64, y: &[f64; LEN_Y], delta_y: &mut [f64; LEN_Y], act: &[bool; LEN_B]) {
72    if act[0] {
73      delta_y[2] += self.p[8];
74    }
75    if act[1] {
76      delta_y[2] -= self.p[8];
77    }
78  }
79
80  #[allow(unused_variables)]
81  fn cond(
82    &self,
83    dec_t: &Decimal,
84    act: &mut [bool; LEN_B],
85    next_t: &[Decimal; LEN_B],
86    y: &[f64; LEN_Y],
87  ) {
88    act[0] = if *dec_t == next_t[0] { true } else { false };
89    act[1] = if *dec_t == next_t[1] { true } else { false };
90  }
91
92  #[allow(unused_variables)]
93  fn beat(&self, t: &f64, y: &[f64; LEN_Y]) -> [[Decimal; 3]; LEN_B] {
94    [
95      beat![self.p[2], self.p[3], self.p[4]],
96      beat![self.p[5], self.p[6], self.p[7]],
97    ]
98  }
99
100  #[allow(unused_variables)]
101  fn cre(&self, t: &f64, y: &mut [f64; LEN_Y]) {
102    y[3] = self.p[8] * t;
103  }
104}
105
106fn smptime() -> Vec<f64> {
107  let mut vec_smptime = Vec::new();
108  for i in 0..=5000 {
109    vec_smptime.push(i as f64 / 100.0);
110  }
111  vec_smptime
112}
sim/sim.rs

sim/
sim.rs
