sbml-sim 0.1.2

Simulator for SBML models.
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

use crate::structs::bindings::Bindings;
use std::collections::HashMap;

#[allow(non_snake_case)]
pub fn runge_kutta_4(bindings: &Bindings, step_size: f64) -> Result<HashMap<String, f64>, String> {
    //
    let original_bindings = bindings.values();
    let mut k: HashMap<String, Vec<f64>> = HashMap::new();

    // K1
    let k1_assignments = HashMap::<String, f64>::new();
    for ODE in &bindings.ODEs {
        let dependent_variable = &ODE.id;
        k.insert(dependent_variable.clone(), Vec::new());
        let ode_value = ODE.evaluate(&k1_assignments, bindings)?;
        let k1 = step_size * ode_value;
        k.entry(dependent_variable.clone())
            .and_modify(|v| v.push(k1));
    }
    drop(k1_assignments);

    let mut k2_assignments: HashMap<String, f64> = HashMap::new();
    for (dependent_variable, k_current) in &k {
        let k1 = k_current[0];
        if let Some(original_value) = original_bindings.get(dependent_variable) {
            k2_assignments.insert(dependent_variable.to_owned(), original_value + k1 / 2.0);
        }
    }

    for (key, val) in bindings.emulate_assignment_rules(&k2_assignments) {
        k2_assignments.insert(key, val);
    }

    for ODE in &bindings.ODEs {
        let dependent_variable = &ODE.id;
        let ode_value = ODE.evaluate(&k2_assignments, bindings)?;
        let k2 = step_size * ode_value;
        k.entry(dependent_variable.clone())
            .and_modify(|v| v.push(k2));
    }
    drop(k2_assignments);

    let mut k3_assignments = HashMap::<String, f64>::new();
    for (dependent_variable, k_current) in &k {
        let k2 = k_current[1];
        if let Some(original_value) = original_bindings.get(dependent_variable) {
            k3_assignments.insert(dependent_variable.to_owned(), original_value + k2 / 2.0);
        }
    }

    for (key, val) in bindings.emulate_assignment_rules(&k3_assignments) {
        k3_assignments.insert(key, val);
    }

    // Calculate k3 values
    for ODE in &bindings.ODEs {
        let dependent_variable = &ODE.id;
        let ode_value = ODE.evaluate(&k3_assignments, bindings)?;
        let k3 = step_size * ode_value;
        k.entry(dependent_variable.clone())
            .and_modify(|v| v.push(k3));
    }
    drop(k3_assignments);

    let mut k4_assignments = HashMap::<String, f64>::new();
    for (dependent_variable, k_current) in &k {
        let k3 = k_current[2];
        if let Some(original_value) = original_bindings.get(dependent_variable) {
            k4_assignments.insert(dependent_variable.to_owned(), original_value + k3);
        }
    }

    for (key, val) in bindings.emulate_assignment_rules(&k4_assignments) {
        k4_assignments.insert(key, val);
    }

    // Calculate k4 values
    for ODE in &bindings.ODEs {
        let dependent_variable = &ODE.id;
        let ode_value = ODE.evaluate(&k4_assignments, bindings)?;
        let k4 = step_size * ode_value;
        k.entry(dependent_variable.clone())
            .and_modify(|v| v.push(k4));
    }
    drop(k4_assignments);

    // Calculate final changes to derivatives
    let mut deltas: HashMap<String, f64> = HashMap::new();
    for (dependent_variable, k_current) in &k {
        let k1 = k_current[0];
        let k2 = k_current[1];
        let k3 = k_current[2];
        let k4 = k_current[3];

        let delta = (k1 + k4) / 6.0 + (k2 + k3) / 3.0;
        deltas.insert(dependent_variable.to_owned(), delta);
    }

    Ok(deltas)
}