Coverage for physioblocks / utils / gradient_test_utils.py: 91%

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/>. 

26 

27""" 

28Define functions to perform gradient tests to check the coherence between an 

29function and its gradient. 

30 

31The gradient test principle is to compare the result of a **computed** and estimated 

32gradient from a :class:`~physioblocks.computing.assembling.EqSystem` object. 

33""" 

34 

35from __future__ import annotations 

36 

37import logging 

38from os import linesep 

39from typing import Any 

40 

41import numpy as np 

42from numpy.typing import NDArray 

43from scipy.optimize import approx_fprime 

44 

45import physioblocks.simulation.setup as setup 

46from physioblocks.computing.assembling import EqSystem 

47from physioblocks.computing.models import Block, Expression, ModelComponent 

48from physioblocks.configuration.aliases import unwrap_aliases 

49from physioblocks.configuration.functions import load 

50from physioblocks.io.configuration import read_json 

51from physioblocks.simulation.runtime import AbstractSimulation 

52from physioblocks.simulation.state import State 

53 

54_logger = logging.getLogger(__name__) 

55 

56_ZERO = 1e-16 

57_REL_TOL = 1e-6 

58_NEW_STATE_SHIFT_FACTOR = 1e-3 

59_SHIFT_FACTOR = 1e-4 

60_CONVERGENCE_FACTOR = 0.1 

61_CONVERGENCE_TOL = 0.11 

62 

63 

64def gradient_test_from_file(config_file_path: str) -> bool: 

65 """ 

66 Read the given configuration file and perform a gradient test on the full Net 

67 with the provided parameters and variables. 

68 

69 .. warning:: 

70 

71 Every type (especially **Blocks**, **ModelComponents** and **Functions**) and 

72 alias used in the configuration file have to be loaded. 

73 Follow this 

74 :ref:`user guide section <user_guide_level_3_block_test_dynamic_import>` 

75 to import libraries and aliases dynamically. 

76 

77 :param config_file_path: the file path to the simulation configuration file. 

78 :type config_file_path: str 

79 

80 :return: True if the gradient test is successfull, false otherwise. 

81 """ 

82 configuration = read_json(config_file_path) 

83 configuration = unwrap_aliases(configuration) 

84 sim: AbstractSimulation = load(configuration) 

85 sim._initialize() # noqa SLF001 

86 return gradient_test( 

87 sim.eq_system, 

88 sim.state, 

89 sim.solver._get_state_magnitude(sim.state, sim.magnitudes), # noqa: SLF001 

90 ) 

91 

92 

93def gradient_test_from_model( 

94 model: ModelComponent, state: State, state_magnitude: NDArray[np.float64] 

95) -> bool: 

96 """ 

97 Create an equation system for the given block only and perform a gradient test. 

98 

99 .. note:: 

100 

101 It does not test the submodels of the model. 

102 

103 :param model: the model to test. 

104 :type model: str 

105 

106 :param state: the state used to determine the variables in the model. 

107 :type state: str 

108 

109 :param state_magnitude: the state variables magnitudes 

110 :type state_magnitude: str 

111 

112 :return: True if the gradient test is successfull, false otherwise. 

113 """ 

114 line_index = 0 

115 expressions = setup.SystemExpressions() 

116 for internal_expr_def in type(model).internal_expressions: 

117 expressions.append((line_index, internal_expr_def.expression, model)) 

118 line_index += internal_expr_def.expression.size 

119 

120 if isinstance(model, Block): 

121 for fluxes_expr_def in type(model).fluxes_expressions.values(): 

122 expressions.append((line_index, fluxes_expr_def.expression, model)) 

123 line_index += fluxes_expr_def.expression.size 

124 

125 eq_system = setup.build_eq_system(expressions, state) 

126 model.initialize() 

127 return gradient_test(eq_system, state, state_magnitude) 

128 

129 

130def gradient_test_from_expression( 

131 expr: Expression, 

132 expr_params: Any, 

133 state: State, 

134 state_magnitude: NDArray[np.float64], 

135) -> bool: 

136 """ 

137 Create an equation system for the given expression only and perform a gradient test. 

138 

139 :param expr: the expression to test. 

140 :type expr: Expression 

141 

142 :param expr_params: the parameters to pass to the expression 

143 :type expr_params: Any 

144 

145 :param state: the state used to determine the variables in the expression. 

146 :type state: str 

147 

148 :param state_magnitude: the state variables magnitudes 

149 :type state_magnitude: str 

150 

151 :return: True if the gradient test is successfull, false otherwise. 

152 :rtype: bool 

153 """ 

154 

155 eq_system = setup.build_eq_system([(0, expr, expr_params)], state) 

156 return gradient_test(eq_system, state, state_magnitude) 

157 

158 

159def gradient_test( 

160 eq_system: EqSystem, state: State, state_magnitude: NDArray[np.float64] 

161) -> bool: 

162 """ 

163 Test the computed gradient for the equation system by comparing it to 

164 a gradient estimated with finite differences. 

165 

166 :param eq_system: the equation system providing method to compute a residual and a 

167 gradient 

168 :type eq_system: EqSystem 

169 

170 :param state: system state 

171 :type state: State 

172 

173 :param state_magnitude: the state variables magnitudes 

174 :type state_magnitude: str 

175 

176 :return: True if the estimated and computed gradient meet tolerance, 

177 False otherwise. 

178 """ 

179 

180 _logger.info(str.format("State:{0}{1}", linesep, state)) 

181 _logger.info(str.format("System:{0}{1}", linesep, eq_system)) 

182 

183 new_state = state.state_vector + _NEW_STATE_SHIFT_FACTOR * state_magnitude 

184 

185 shift_1 = _SHIFT_FACTOR * state_magnitude 

186 state.update_state_vector(new_state) 

187 

188 res = eq_system.compute_residual() 

189 grad = eq_system.compute_gradient() 

190 

191 grad_estimated_1 = _estimate_gradient(eq_system, shift_1, state) 

192 

193 shift_2 = _CONVERGENCE_FACTOR * shift_1 

194 state.update_state_vector(new_state) 

195 grad_estimated_2 = _estimate_gradient(eq_system, shift_2, state) 

196 

197 state.update_state_vector(new_state) 

198 

199 with np.printoptions(precision=9, linewidth=10000, suppress=True): 

200 _logger.debug(str.format("Residual:{0}{1}", linesep, res)) 

201 _logger.debug(str.format("Gradient{0}{1}", linesep, grad)) 

202 _logger.debug( 

203 str.format("Gradient Estimated:{0}{1}", linesep, grad_estimated_1) 

204 ) 

205 

206 grad_error_pos = get_errors_gradient(grad, grad_estimated_1, grad_estimated_2) 

207 

208 error_message = str.format("State:{0}{1}", linesep, state) 

209 error_message += str.format("System:{0}{1}{0}{0}", linesep, eq_system) 

210 

211 for line_error, col_error in grad_error_pos: 

212 error_msg_line = str.format( 

213 "Error in gradient at ({0},{1}).{2}", line_error, col_error, linesep 

214 ) 

215 _logger.debug(error_msg_line) 

216 error_message += error_msg_line 

217 

218 error_msg_line = str.format( 

219 "Variable Id: {0}.{1}", state.get_variable_id(col_error), linesep 

220 ) 

221 _logger.debug(error_msg_line) 

222 

223 error_msg_line = str.format("Residual Index: {0}.{1}{1}", line_error, linesep) 

224 _logger.debug(error_msg_line) 

225 

226 if len(grad_error_pos) != 0: 

227 raise GradientError(error_message) 

228 

229 return len(grad_error_pos) == 0 

230 

231 

232def _estimate_gradient( 

233 eq_system: EqSystem, shift: NDArray[np.float64], state: State 

234) -> NDArray[np.float64]: 

235 """ 

236 Gradient finite difference estimation 

237 

238 :param eq_system: the equation system providing method to compute a residual 

239 and a gradient 

240 :type eq_system: EqSystem 

241 

242 :param shift_factor: the shift factor to apply on the state (can be 0). 

243 It is multiplied with the variables magnitude then added to the state. 

244 :type shift_factor: NDArray 

245 

246 :param state: system state 

247 :type state: State 

248 """ 

249 

250 def residual(x: NDArray[np.float64]) -> NDArray[np.float64]: 

251 state.update_state_vector(x) 

252 res = eq_system.compute_residual() 

253 

254 return res 

255 

256 x0 = state.state_vector 

257 grad_est = np.atleast_2d(approx_fprime(x0, residual, shift)) 

258 

259 return grad_est 

260 

261 

262def get_errors_gradient( 

263 computed: NDArray[np.float64], 

264 estimated_1: NDArray[np.float64], 

265 estimated_2: NDArray[np.float64], 

266) -> list[tuple[int, int]]: 

267 """ 

268 Get the ``(line, column)`` positions of the errors in the gradient matrix. 

269 

270 :param computed: computed gradient 

271 :type computed: NDArray[np.float64] 

272 

273 :param estimated_1: first gradient estimate 

274 :type estimated_1: NDArray[np.float64] 

275 

276 :param estimated_2: second gradient estimate 

277 :type estimated_2: NDArray[np.float64] 

278 

279 :return: the position of the errors in the gradient matrix 

280 :rtype: list[tuple[int, int]] 

281 """ 

282 abs_computed = np.abs(computed) 

283 

284 error_1 = np.empty(computed.shape) 

285 error_1[:] = np.inf 

286 error_2 = error_1.copy() 

287 

288 for grad_line in range(0, computed.shape[0]): 

289 for grad_col in range(0, computed.shape[1]): 

290 if abs_computed[grad_line, grad_col] < _ZERO: 

291 error_1[grad_line, grad_col] = ( 

292 np.abs( 

293 computed[grad_line, grad_col] - estimated_1[grad_line, grad_col] 

294 ) 

295 / _ZERO 

296 ) 

297 error_2[grad_line, grad_col] = ( 

298 np.abs( 

299 computed[grad_line, grad_col] - estimated_2[grad_line, grad_col] 

300 ) 

301 / _ZERO 

302 ) 

303 else: 

304 error_1[grad_line, grad_col] = ( 

305 np.abs( 

306 computed[grad_line, grad_col] - estimated_1[grad_line, grad_col] 

307 ) 

308 / abs_computed[grad_line, grad_col] 

309 ) 

310 error_2[grad_line, grad_col] = ( 

311 np.abs( 

312 computed[grad_line, grad_col] - estimated_2[grad_line, grad_col] 

313 ) 

314 / abs_computed[grad_line, grad_col] 

315 ) 

316 

317 with np.printoptions(precision=9, linewidth=10000, suppress=False): 

318 _logger.debug(str.format("Error:{0}{1}", error_1, linesep)) 

319 

320 error_positions = [] 

321 

322 for grad_line in range(0, computed.shape[0]): 

323 for grad_col in range(0, computed.shape[1]): 

324 if error_1[grad_line, grad_col] > _REL_TOL: 

325 err_div = error_2[grad_line, grad_col] / error_1[grad_line, grad_col] 

326 if err_div > _CONVERGENCE_TOL: 

327 error_positions.append((grad_line, grad_col)) 

328 elif np.isnan(error_1[grad_line, grad_col]): 

329 error_positions.append((grad_line, grad_col)) 

330 

331 return error_positions 

332 

333 

334class GradientError(Exception): 

335 """Error raised when the estimated and computed gradients do not match.""" 

336 

337 pass