diffsol 0.17.0

A library for solving ordinary differential equations (ODEs) in Rust.
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

use crate::{
    jacobian::{find_constant_sens_non_zeros, JacobianColoring},
    matrix::sparsity::MatrixSparsity,
    ConstantOp, ConstantOpSens, Matrix, Op, Vector,
};

use super::{BuilderOp, ParameterisedOp};

pub struct ConstantClosureWithSens<M, I, J>
where
    M: Matrix,
    I: Fn(&[M::T], M::T, &mut [M::T]),
    J: Fn(&[M::T], M::T, &[M::T], &mut [M::T]),
{
    func: I,
    func_sens: J,
    nout: usize,
    nparams: usize,
    sens_sparsity: Option<M::Sparsity>,
    sens_coloring: Option<JacobianColoring<M>>,
    ctx: M::C,
}

impl<M, I, J> ConstantClosureWithSens<M, I, J>
where
    M: Matrix,
    I: Fn(&[M::T], M::T, &mut [M::T]),
    J: Fn(&[M::T], M::T, &[M::T], &mut [M::T]),
{
    pub fn new(func: I, func_sens: J, nout: usize, nparams: usize, ctx: M::C) -> Self {
        Self {
            func,
            func_sens,
            nout,
            nparams,
            sens_sparsity: None,
            sens_coloring: None,
            ctx,
        }
    }

    /// Find the sparsity of the parameter Jacobian dy0/dp, and a coloring to compute it with.
    pub fn calculate_sens_sparsity(&mut self, t0: M::T, p: &M::V) {
        let op = ParameterisedOp { op: self, p };
        let non_zeros = find_constant_sens_non_zeros(&op, t0);
        self.sens_sparsity = Some(
            MatrixSparsity::try_from_indices(self.nout(), p.len(), non_zeros.clone())
                .expect("invalid sparsity pattern"),
        );
        self.sens_coloring = Some(JacobianColoring::new(
            self.sens_sparsity.as_ref().unwrap(),
            &non_zeros,
            self.ctx.clone(),
        ));
    }
}

impl<M, I, J> Op for ConstantClosureWithSens<M, I, J>
where
    M: Matrix,
    I: Fn(&[M::T], M::T, &mut [M::T]),
    J: Fn(&[M::T], M::T, &[M::T], &mut [M::T]),
{
    type V = M::V;
    type T = M::T;
    type M = M;
    type C = M::C;
    fn nstates(&self) -> usize {
        0
    }
    fn nout(&self) -> usize {
        self.nout
    }
    fn nparams(&self) -> usize {
        self.nparams
    }
    fn context(&self) -> &Self::C {
        &self.ctx
    }
}

impl<M, I, J> BuilderOp for ConstantClosureWithSens<M, I, J>
where
    M: Matrix,
    I: Fn(&[M::T], M::T, &mut [M::T]),
    J: Fn(&[M::T], M::T, &[M::T], &mut [M::T]),
{
    fn calculate_sparsity(&mut self, _y0: &Self::V, _t0: Self::T, _p: &Self::V) {
        // do nothing
    }
    fn calculate_augmented_sparsity(&mut self, _y0: &Self::V, t0: Self::T, p: &Self::V) {
        ConstantClosureWithSens::calculate_sens_sparsity(self, t0, p);
    }
    fn set_nstates(&mut self, _nstates: usize) {
        // do nothing
    }
    fn set_nout(&mut self, nout: usize) {
        self.nout = nout;
    }
    fn set_nparams(&mut self, nparams: usize) {
        self.nparams = nparams;
    }
}

impl<M, I, J> ConstantOp for ParameterisedOp<'_, ConstantClosureWithSens<M, I, J>>
where
    M: Matrix,
    I: Fn(&[M::T], M::T, &mut [M::T]),
    J: Fn(&[M::T], M::T, &[M::T], &mut [M::T]),
{
    fn call_inplace(&self, t: Self::T, y: &mut Self::V) {
        y.for_each_batch_host([self.p], |y, [p], _| (self.op.func)(p, t, y));
    }
}

impl<M, I, J> ConstantOpSens for ParameterisedOp<'_, ConstantClosureWithSens<M, I, J>>
where
    M: Matrix,
    I: Fn(&[M::T], M::T, &mut [M::T]),
    J: Fn(&[M::T], M::T, &[M::T], &mut [M::T]),
{
    fn sens_mul_inplace(&self, t: Self::T, v: &Self::V, y: &mut Self::V) {
        y.for_each_batch_host([self.p, v], |y, [p, v], _| (self.op.func_sens)(p, t, v, y));
    }

    fn sens_inplace(&self, t: Self::T, y: &mut Self::M) {
        if let Some(coloring) = self.op.sens_coloring.as_ref() {
            coloring.constant_sens_inplace(self, t, y);
        } else {
            self._default_sens_inplace(t, y);
        }
    }

    fn sens_sparsity(&self) -> Option<<Self::M as Matrix>::Sparsity> {
        self.op.sens_sparsity.clone()
    }
}