Coverage for physioblocks / computing / assembling.py: 100%
74 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-09 16:40 +0100
1# SPDX-FileCopyrightText: Copyright INRIA
2#
3# SPDX-License-Identifier: LGPL-3.0-only
4#
5# Copyright INRIA
6#
7# This file is part of PhysioBlocks, a library mostly developed by the
8# [Ananke project-team](https://team.inria.fr/ananke) at INRIA.
9#
10# Authors:
11# - Colin Drieu
12# - Dominique Chapelle
13# - François Kimmig
14# - Philippe Moireau
15#
16# PhysioBlocks is free software: you can redistribute it and/or modify it under the
17# terms of the GNU Lesser General Public License as published by the Free Software
18# Foundation, version 3 of the License.
19#
20# PhysioBlocks is distributed in the hope that it will be useful, but WITHOUT ANY
21# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
22# PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
23#
24# You should have received a copy of the GNU Lesser General Public License along with
25# PhysioBlocks. If not, see <https://www.gnu.org/licenses/>.
27"""
28Allows to build systems that computes residual and gradient from collections
29of functions.
30"""
32from __future__ import annotations
34from dataclasses import dataclass
35from pprint import pformat
36from typing import Any
38import numpy as np
39from numpy.typing import NDArray
41from physioblocks.computing.models import ModelComponent, SystemFunction
44@dataclass
45class _EqSystemPart:
46 """
47 Part of the global equation system.
49 It computes a term that can be summed in the global system at the provided index.
50 """
52 line_index: int
53 """Index where to sum the part in the global system"""
54 res_part_size: int
55 """Number of lines in the result of the residual part function"""
56 model_component: ModelComponent
57 """The model that holds the data use to compute the residual and gradient part"""
58 res_func: SystemFunction
59 """The function to compute the residual part"""
60 column_index: int
61 """Start column index of the gradient part"""
62 grad_part_size: int
63 """Number of colums in the result of the gradient part function"""
64 grad_funcs: dict[int, SystemFunction]
65 """A collection of fonctions to compute the gradient part for each variables,
66 associated with their column index"""
68 def compute_residual_part(self) -> np.float64 | NDArray[np.float64]:
69 """
70 Compute the residual part numerical value.
72 Can be a scalar or vector of size _res_part_size.
74 :return: the residual part value
75 :rtype: np.float64 | NDArray[np.float64]
76 """
77 return self.res_func(self.model_component)
79 def compute_gradient_part(self) -> NDArray[np.float64]:
80 """
81 Compute the gradient part numerical value.
83 It a matrix of size (_res_part_size, _grad_part_size) to sum in the
84 global gradient at index (_index, _min_grad_index).
86 The numerical values are evaluated where the gradient functions are provided,
87 the remaing values are set to 0.
89 :return: the gradient part numerical value
90 :rtype: np.float64 | NDArray[np.float64]
91 """
92 grad_part = np.zeros(
93 (self.res_part_size, self.grad_part_size),
94 )
96 for grad_index in self.grad_funcs:
97 gradient_func = self.grad_funcs[grad_index]
98 grad_colum_index = grad_index - self.column_index
99 grad_part[:, grad_colum_index] += gradient_func(self.model_component)
101 return grad_part
104def _create_eq_system_part(
105 residual_index: int,
106 residual_size: int,
107 residual_func: SystemFunction,
108 gradients_func: dict[int, SystemFunction],
109 model_component: ModelComponent,
110) -> _EqSystemPart:
111 """
112 Create an equation system part.
114 :param index: the line index where to sum the expression in the global
115 system
116 :type index: int
118 :param residual_size: the size of the result of the residual_func
119 :type residual_size: int
121 :param residual_func: the function to compute the residual part
122 :type residual_func: SystemFunction
124 :param gradients_func: the functions to compute the gradients parts associated with
125 the index of the variable in the state vector
127 :param model_component: the model holding the data needed to compute the expression
128 :type model_component: ModelComponent
130 :return: an equation system part
131 :rtype: _EqSystemPart
132 """
133 grad_part_size = 0
134 min_grad_index = 0
135 if len(gradients_func) > 0:
136 min_grad_index = min(gradients_func.keys())
137 grad_part_size = max(gradients_func) - min_grad_index + 1
139 return _EqSystemPart(
140 residual_index,
141 residual_size,
142 model_component,
143 residual_func,
144 min_grad_index,
145 grad_part_size,
146 gradients_func,
147 )
150class EqSystem:
151 """
152 Global Equation system
154 Allows to compute the numerical value of the residual and the gradient
155 for a :class:`~physioblocks.description.nets.Net`.
157 The system is built from a collection of **System Parts**.
158 The parts are summed in a global system to compute the numerical value of
159 the residual and the gradient.
160 """
162 _system_size: int
163 """The system size, fixed at system creation"""
164 _system_parts: list[_EqSystemPart]
165 """The collection of system parts"""
167 def __init__(self, size: int):
168 """
169 EqSystem constructor
171 :param size: system fixed size
172 :type size: int
173 """
174 self._system_size = size
175 self._system_parts = []
177 def __str__(self) -> str:
178 """
179 Get the equation system string representation.
181 :return: the equation system string representation
182 :rtype: str
183 """
185 eq_system_dict = {
186 index: [
187 (eq_part.res_func.__qualname__, "size " + str(eq_part.res_part_size))
188 for eq_part in self._system_parts
189 if eq_part.line_index == index
190 ]
191 for index in range(0, self.system_size)
192 }
194 return pformat(eq_system_dict, indent=2, compact=False)
196 @property
197 def system_size(self) -> int:
198 """
199 Get the system size
201 :return: the system size
202 :rtype: int
203 """
204 return self._system_size
206 def add_system_part(
207 self,
208 residual_index: int,
209 residual_size: int,
210 residual_func: SystemFunction,
211 gradients_func: dict[int, SystemFunction],
212 parameters: Any,
213 ) -> None:
214 """add_system_part(self, residual_index: int, residual_size: int, residual_func: SystemFunction, gradients_func: dict[int, SystemFunction], parameters: Any)
215 Add an equation system part to the global system.
217 :param residual_index:
218 the line index where to sum the function in the global system
219 :type residual_index: int
221 :param residual_size: the size of the result of the residual_func
222 :type residual_size: int
224 :param residual_func: the function to compute the residual part
225 :type residual_func: SystemFunction
227 :param gradients_func:
228 the functions to compute the gradients parts associated with
229 the index of the variable in the state vector
231 :param parameters: the parameters needed to compute the SystemFunctions
232 (both residual and gradient)
233 (should be the same type as the argument of the residual_func and the
234 gradients_func input type)
235 :type parameters: Any
237 :raise ValueError: Raises a ValueError when either the size of the residual
238 or the gradient part will exceed the global system size.
239 """ # noqa: E501
241 part = _create_eq_system_part(
242 residual_index, residual_size, residual_func, gradients_func, parameters
243 )
245 if part.line_index + part.res_part_size > self.system_size:
246 raise ValueError("The residual part exceed system size")
248 if part.column_index + part.grad_part_size > self.system_size:
249 raise ValueError("The gradient part exceed system size")
251 self._system_parts.append(part)
253 def compute_residual(self) -> NDArray[np.float64]:
254 """
255 Compute the numerical value of the residual.
257 :return: the residual
258 :rtype: NDArray[np.float64]
259 """
261 residual = np.zeros(
262 shape=self.system_size,
263 )
265 for part in self._system_parts:
266 residual_part = part.compute_residual_part()
267 residual_column_index = part.line_index + part.res_part_size
268 residual[part.line_index : residual_column_index] += residual_part
270 return residual
272 def compute_gradient(self) -> NDArray[np.float64]:
273 """
274 Compute the numerical value of the gradient.
276 :return: the gradient
277 :rtype: NDArray[np.float64]
278 """
279 grad = np.zeros(
280 shape=(self.system_size, self.system_size),
281 )
283 for part in self._system_parts:
284 if part.grad_part_size > 0:
285 gradient_part = part.compute_gradient_part()
286 grad[
287 part.line_index : part.line_index + part.res_part_size,
288 part.column_index : part.column_index + part.grad_part_size,
289 ] += gradient_part
291 return grad