From 69027d31aee313263b87163161e199fe1489e400 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 11 May 2023 21:29:49 +0200 Subject: [PATCH 1/9] Extend / cleanup event tests Three new event test cases and some cleanup. To reduce changes in #1539 Co-authored-by: Paul Stapor --- python/tests/test_events.py | 339 ++++++++++++++++++++++++++++---- python/tests/test_heavisides.py | 89 ++++----- python/tests/util.py | 15 +- 3 files changed, 345 insertions(+), 98 deletions(-) diff --git a/python/tests/test_events.py b/python/tests/test_events.py index a6f5334b83..da5487effa 100644 --- a/python/tests/test_events.py +++ b/python/tests/test_events.py @@ -4,15 +4,20 @@ import numpy as np import pytest +from amici.testing import skip_on_valgrind from util import (check_trajectories_with_forward_sensitivities, check_trajectories_without_sensitivities, create_amici_model, create_sbml_model) -from amici.testing import skip_on_valgrind @pytest.fixture(params=[ pytest.param('events_plus_heavisides', marks=skip_on_valgrind), - 'nested_events', + pytest.param('piecewise_plus_event_simple_case', marks=skip_on_valgrind), + pytest.param('piecewise_plus_event_semi_complicated', + marks=skip_on_valgrind), + pytest.param('piecewise_plus_event_trigger_depends_on_state', + marks=skip_on_valgrind), + pytest.param('nested_events', marks=skip_on_valgrind), ]) def model(request): """Returns the requested AMICI model and analytical expressions.""" @@ -23,8 +28,8 @@ def model(request): species, events, timepoints, - x_pected, - sx_pected + x_expected, + sx_expected ) = get_model_definition(request.param) # SBML model @@ -45,22 +50,28 @@ def model(request): ) amici_model.setTimepoints(timepoints) - return amici_model, parameters, timepoints, x_pected, sx_pected + return amici_model, parameters, timepoints, x_expected, sx_expected def get_model_definition(model_name): + if model_name == 'piecewise_plus_event_simple_case': + return model_definition_piecewise_plus_event_simple_case() + if model_name == 'piecewise_plus_event_semi_complicated': + return model_definition_piecewise_plus_event_semi_complicated() + if model_name == 'piecewise_plus_event_trigger_depends_on_state': + return model_definition_piecewise_plus_event_trigger_depends_on_state() if model_name == 'events_plus_heavisides': - return model_definition_events_plus_heavisides() + return model_definition_events_plus_heavisides() if model_name == 'nested_events': return model_definition_nested_events() - else: - raise NotImplementedError( - f'Model with name {model_name} is not implemented.' - ) + + raise NotImplementedError( + f'Model with name {model_name} is not implemented.' + ) def model_definition_events_plus_heavisides(): - """Test model for state- and parameter-dependent heavisides. + """Test model for state- and parameter-dependent Heavisides. ODEs ---- @@ -124,7 +135,7 @@ def model_definition_events_plus_heavisides(): timepoints = np.linspace(0, 8, 400) # Analytical solution - def x_pected(t, k1, k2, k3, alpha, beta, gamma, delta, eta, zeta): + def x_expected(t, k1, k2, k3, alpha, beta, gamma, delta, eta, zeta): # The system reads dx/dt = Ax + b # x0 = (k1, k2, k3) x0 = np.array([[k1], [k2], [k3]]) @@ -180,25 +191,27 @@ def get_early_x(t): x3 += np.array([[0], [0], [zeta / 3]]) hom_x = np.matmul(expm((t - event_3_time) * A), x3) - inhom_x = [[0], [0], - [-np.exp(-eta * (t - event_3_time)) / (eta) - + 1 / (eta)]] + inhom_x = [ + [0], + [0], + [-np.exp(-eta * (t - event_3_time)) / eta + 1 / eta], + ] x = (hom_x + inhom_x).flatten() return np.array(x) - def sx_pected(t, parameters): - # get sx, w.r.t. parameters, via finite differences + def sx_expected(t, parameters): + """get sx, w.r.t. parameters, via finite differences""" sx = [] + eps = 1e-6 for ip in parameters: - eps = 1e-6 perturbed_params = deepcopy(parameters) perturbed_params[ip] += eps - sx_p = x_pected(t, **perturbed_params) + sx_p = x_expected(t, **perturbed_params) perturbed_params[ip] -= 2*eps - sx_m = x_pected(t, **perturbed_params) + sx_m = x_expected(t, **perturbed_params) sx.append((sx_p - sx_m) / (2 * eps)) return np.array(sx) @@ -210,8 +223,8 @@ def sx_pected(t, parameters): species, events, timepoints, - x_pected, - sx_pected + x_expected, + sx_expected ) @@ -252,7 +265,7 @@ def model_definition_nested_events(): 'inflow_1': 4, 'decay_1': 2, 'decay_2': 5, - 'bolus': 0, # for bolus != 0, nested event sensitivities are off! + 'bolus': 0, # for bolus != 0, nested event sensitivities are off! } events = { 'event_1': { @@ -269,13 +282,13 @@ def model_definition_nested_events(): timepoints = np.linspace(0, 1, 101) # Analytical solution - def x_pected(t, k1, k2, inflow_1, decay_1, decay_2, bolus): + def x_expected(t, k1, k2, inflow_1, decay_1, decay_2, bolus): # gather temporary variables # event_time = x_1 > inflow_1 / decay_2 equil = inflow_1 / decay_1 tmp1 = inflow_1 / decay_2 - inflow_1 / decay_1 tmp2 = k1 - inflow_1 / decay_1 - event_time = (- 1 / decay_1) * np.log( tmp1 / tmp2) + event_time = (- 1 / decay_1) * np.log(tmp1 / tmp2) def get_early_x(t): # compute dynamics before event @@ -301,17 +314,17 @@ def get_early_x(t): return x.flatten() - def sx_pected(t, parameters): - # get sx, w.r.t. parameters, via finite differences + def sx_expected(t, parameters): + """get sx, w.r.t. parameters, via finite differences""" sx = [] + eps = 1e-6 for ip in parameters: - eps = 1e-6 perturbed_params = deepcopy(parameters) perturbed_params[ip] += eps - sx_p = x_pected(t, **perturbed_params) + sx_p = x_expected(t, **perturbed_params) perturbed_params[ip] -= 2*eps - sx_m = x_pected(t, **perturbed_params) + sx_m = x_expected(t, **perturbed_params) sx.append((sx_p - sx_m) / (2 * eps)) return np.array(sx) @@ -323,20 +336,275 @@ def sx_pected(t, parameters): species, events, timepoints, - x_pected, - sx_pected + x_expected, + sx_expected + ) + + +def model_definition_piecewise_plus_event_simple_case(): + """Test model for boolean operations in a piecewise condition. + + ODEs + ---- + d/dt x_1: + - { 1, (alpha <= t and t < beta) + - { 0, otherwise + """ + # Model components + species = ['x_1'] + initial_assignments = {'x_1': 'x_1_0'} + rate_rules = { + 'x_1': 'piecewise(1, (alpha < time && time < beta), 0)' + } + parameters = { + 'alpha': 2, + 'beta': 3, + 'gamma': 4.5, + 'x_1_0': 1, + } + timepoints = np.linspace(0., 5., 100) # np.array((0.0, 4.0,)) + events = { + 'event_1': { + 'trigger': 'time > alpha', + 'target': 'x_1', + 'assignment': 'gamma' + }, + 'event_2': { + 'trigger': 'time > beta', + 'target': 'x_1', + 'assignment': 'x_1 + 2.5' + } + } + + # Analytical solution + def x_expected(t, x_1_0, alpha, beta, gamma): + t_event_1 = alpha + t_event_2 = beta + + if t < t_event_1: + x = x_1_0 + elif t < t_event_2: + x = gamma + t - t_event_1 + else: + x = gamma + t_event_2 - t_event_1 + 2.5 + + return np.array((x,)) + + def sx_expected(t, parameters): + """get sx, w.r.t. parameters, via finite differences""" + sx = [] + eps = 1e-6 + + for ip in parameters: + perturbed_params = deepcopy(parameters) + perturbed_params[ip] += eps + sx_p = np.array(x_expected(t, **perturbed_params)) + perturbed_params[ip] -= 2*eps + sx_m = np.array(x_expected(t, **perturbed_params)) + sx.append((sx_p - sx_m) / (2 * eps)) + + return np.array(sx) + + return ( + initial_assignments, + parameters, + rate_rules, + species, + events, + timepoints, + x_expected, + sx_expected + ) + + +def model_definition_piecewise_plus_event_semi_complicated(): + """Test model for boolean operations in a piecewise condition, discrete + events and a non-vanishing quadrature for the adjoint state. + """ + # Model components + species = ['x_1', 'x_2'] + initial_assignments = {'x_1': 'x_1_0', + 'x_2': 'x_2_0'} + rate_rules = { + 'x_1': + 'piecewise(delta * x_1, (alpha < time && time < beta), - x_1)', + 'x_2': '- eta * x_2', + } + parameters = { + 'alpha': 2, + 'beta': 3, + 'gamma': 4.5, + 'x_1_0': 1, + 'x_2_0': 5, + 'delta': 2.5, + 'eta': 1.4 + } + timepoints = np.linspace(0., 5., 100) + events = { + 'event_1': { + 'trigger': 'time > alpha / 2', + 'target': 'x_1', + 'assignment': 'gamma' + }, + 'event_2': { + 'trigger': 'time > beta', + 'target': 'x_1', + 'assignment': 'x_1 + x_2' + } + } + + # Analytical solution + def x_expected(t, x_1_0, x_2_0, alpha, beta, gamma, delta, eta): + t_event_1 = alpha / 2 + t_event_2 = beta + heaviside_1 = alpha + + x_2 = x_2_0 * np.exp(- eta * t) + + if t < t_event_1: + x_1 = x_1_0 * np.exp(- t) + elif t < heaviside_1: + x_1 = gamma * np.exp(- (t - t_event_1)) + elif t < t_event_2: + x_1_heaviside_1 = gamma * np.exp(- (heaviside_1 - t_event_1)) + x_1 = x_1_heaviside_1 * np.exp(delta * (t - heaviside_1)) + else: + x_1_heaviside_1 = gamma * np.exp(- (heaviside_1 - t_event_1)) + x_1_at_event_2 = ( + x_1_heaviside_1 * np.exp(delta * (t_event_2 - heaviside_1)) + ) + x_2_at_event_2 = x_2_0 * np.exp(- eta * t_event_2) + x1_after_event_2 = x_1_at_event_2 + x_2_at_event_2 + x_1 = x1_after_event_2 * np.exp(- (t - t_event_2)) + + return np.array((x_1, x_2)) + + def sx_expected(t, parameters): + """get sx, w.r.t. parameters, via finite differences""" + sx = [] + eps = 1e-6 + + for ip in parameters: + perturbed_params = deepcopy(parameters) + perturbed_params[ip] += eps + sx_p = np.array(x_expected(t, **perturbed_params)) + perturbed_params[ip] -= 2*eps + sx_m = np.array(x_expected(t, **perturbed_params)) + sx.append((sx_p - sx_m) / (2 * eps)) + + return np.array(sx) + + return ( + initial_assignments, + parameters, + rate_rules, + species, + events, + timepoints, + x_expected, + sx_expected + ) + + +def model_definition_piecewise_plus_event_trigger_depends_on_state(): + """Test model for boolean operations in a piecewise condition. + + ODEs + ---- + d/dt x_1: + - { 1, (alpha <= t and t < beta) + - { 0, otherwise + """ + # Model components + species = ['x_1', 'x_2'] + initial_assignments = {'x_1': 'x_1_0', + 'x_2': 'x_2_0'} + rate_rules = { + 'x_1': 'piecewise(1, (alpha < time && time < beta), 0)', + 'x_2': '- x_2', + } + parameters = { + 'alpha': 2, + 'beta': 3, + 'gamma': 4.5, + 'x_1_0': 1, + 'x_2_0': 5, + } + timepoints = np.linspace(0., 5., 100) + events = { + 'event_1': { + 'trigger': 'x_1 > 1.4', + 'target': 'x_1', + 'assignment': 'x_1 + gamma' + }, + 'event_2': { + 'trigger': 'time > beta', + 'target': 'x_1', + 'assignment': 'x_1 + x_2' + } + } + + # Analytical solution + def x_expected(t, x_1_0, x_2_0, alpha, beta, gamma): + + heaviside_1 = alpha + t_event_1 = alpha + 1.4 - x_1_0 + t_event_2 = beta + # This should hold in order that the analytical solution is correct + assert heaviside_1 < t_event_1 + + # x_2 never gets perturbed + x_2 = x_2_0 * np.exp(- t) + + if t < heaviside_1: + x_1 = x_1_0 + elif t < t_event_1: + x_1 = (t - heaviside_1) + x_1_0 + elif t < t_event_2: + x_1 = gamma + (t - heaviside_1) + x_1_0 + else: + x_2_at_event_2 = x_2_0 * np.exp(- t_event_2) + x_1_at_event_2 = gamma + (t_event_2 - heaviside_1) + x_1_0 + x_1 = x_1_at_event_2 + x_2_at_event_2 + + return np.array((x_1, x_2)) + + def sx_expected(t, parameters): + """get sx, w.r.t. parameters, via finite differences""" + sx = [] + eps = 1e-6 + + for ip in parameters: + perturbed_params = deepcopy(parameters) + perturbed_params[ip] += eps + sx_p = np.array(x_expected(t, **perturbed_params)) + perturbed_params[ip] -= 2*eps + sx_m = np.array(x_expected(t, **perturbed_params)) + sx.append((sx_p - sx_m) / (2 * eps)) + + return np.array(sx) + + return ( + initial_assignments, + parameters, + rate_rules, + species, + events, + timepoints, + x_expected, + sx_expected ) def test_models(model): - amici_model, parameters, timepoints, x_pected, sx_pected = model + amici_model, parameters, timepoints, x_expected, sx_expected = model result_expected_x = np.array([ - x_pected(t, **parameters) + x_expected(t, **parameters) for t in timepoints ]) result_expected_sx = np.array([ - sx_pected(t, parameters) + sx_expected(t, parameters) for t in timepoints ]) @@ -347,7 +615,6 @@ def test_models(model): result_expected_x, result_expected_sx) - def expm(x): """``expm`` wrapper diff --git a/python/tests/test_heavisides.py b/python/tests/test_heavisides.py index 3d1a3564b4..8f1040fe67 100644 --- a/python/tests/test_heavisides.py +++ b/python/tests/test_heavisides.py @@ -2,13 +2,10 @@ import numpy as np import pytest +from util import (check_trajectories_with_forward_sensitivities, + check_trajectories_without_sensitivities, create_amici_model, + create_sbml_model) -from util import ( - create_sbml_model, - create_amici_model, - check_trajectories_without_sensitivities, - check_trajectories_with_forward_sensitivities, -) @pytest.fixture(params=[ 'state_and_param_dep_heavisides', @@ -24,8 +21,8 @@ def model(request): species, events, timepoints, - x_pected, - sx_pected + x_expected, + sx_expected ) = get_model_definition(request.param) # SBML model @@ -46,18 +43,18 @@ def model(request): ) amici_model.setTimepoints(timepoints) - return amici_model, parameters, timepoints, x_pected, sx_pected + return amici_model, parameters, timepoints, x_expected, sx_expected def test_models(model): - amici_model, parameters, timepoints, x_pected, sx_pected = model + amici_model, parameters, timepoints, x_expected, sx_expected = model result_expected_x = np.array([ - x_pected(t, **parameters) + x_expected(t, **parameters) for t in timepoints ]) result_expected_sx = np.array([ - sx_pected(t, **parameters) + sx_expected(t, **parameters) for t in timepoints ]) @@ -115,7 +112,7 @@ def model_definition_state_and_parameter_dependent_heavisides(): events = {} # Analytical solution - def x_pected(t, alpha, beta, gamma, delta, eta, zeta): + def x_expected(t, alpha, beta, gamma, delta, eta, zeta): # get x_1 tau_1 = (np.exp(gamma * delta) - delta * eta) / (1 - eta) if t < tau_1: @@ -130,9 +127,9 @@ def x_pected(t, alpha, beta, gamma, delta, eta, zeta): else: x_2 = np.exp(gamma*delta) + eta*(t-delta) - return (x_1, x_2) + return x_1, x_2 - def sx_pected(t, alpha, beta, gamma, delta, eta, zeta): + def sx_expected(t, alpha, beta, gamma, delta, eta, zeta): # get sx_1, w.r.t. parameters tau_1 = (np.exp(gamma * delta) - delta * eta) / (1 - eta) if t < tau_1: @@ -172,20 +169,17 @@ def sx_pected(t, alpha, beta, gamma, delta, eta, zeta): # get sx_2, w.r.t. parameters tau_2 = delta + sx_2_alpha = 0 + sx_2_beta = 0 + sx_2_zeta = 0 if t < tau_2: - sx_2_alpha = 0 - sx_2_beta = 0 sx_2_gamma = t * np.exp(gamma*t) sx_2_delta = 0 sx_2_eta = 0 - sx_2_zeta = 0 else: - sx_2_alpha = 0 - sx_2_beta = 0 sx_2_gamma = delta * np.exp(gamma*delta) sx_2_delta = gamma*np.exp(gamma*delta) - eta sx_2_eta = t - delta - sx_2_zeta = 0 sx_1 = (sx_1_alpha, sx_1_beta, sx_1_gamma, sx_1_delta, sx_1_eta, sx_1_zeta) @@ -201,8 +195,8 @@ def sx_pected(t, alpha, beta, gamma, delta, eta, zeta): species, events, timepoints, - x_pected, - sx_pected + x_expected, + sx_expected ) @@ -239,37 +233,26 @@ def model_definition_piecewise_with_boolean_operations(): events = {} # Analytical solution - def x_pected(t, x_1_0, alpha, beta, gamma, delta): + def x_expected(t, x_1_0, alpha, beta, gamma, delta): if t < alpha: - return (x_1_0,) + return x_1_0, elif alpha <= t < beta: - return (x_1_0 + (t - alpha),) + return x_1_0 + (t - alpha), elif beta <= t < gamma: - return (x_1_0 + (beta - alpha),) + return x_1_0 + (beta - alpha), elif gamma <= t < delta: - return (x_1_0 + (beta - alpha) + (t - gamma),) + return x_1_0 + (beta - alpha) + (t - gamma), else: - return (x_1_0 + (beta - alpha) + (delta - gamma), ) + return x_1_0 + (beta - alpha) + (delta - gamma), - def sx_pected(t, x_1_0, alpha, beta, gamma, delta): + def sx_expected(t, x_1_0, alpha, beta, gamma, delta): # x0 is very simple... sx_x0 = 1 - sx_alpha = 0 - sx_beta = 0 - sx_gamma = 0 - sx_delta = 0 - - if t >= alpha: - sx_alpha = -1 - if t >= beta: - sx_beta = 1 - if t >= gamma: - sx_gamma = -1 - if t >= delta: - sx_delta = 1 - + sx_alpha = -1 if t >= alpha else 0 + sx_beta = 1 if t >= beta else 0 + sx_gamma = -1 if t >= gamma else 0 + sx_delta = 1 if t >= delta else 0 sx = (sx_alpha, sx_beta, sx_gamma, sx_delta, sx_x0) - return np.array((sx,)).transpose() return ( @@ -279,8 +262,8 @@ def sx_pected(t, x_1_0, alpha, beta, gamma, delta): species, events, timepoints, - x_pected, - sx_pected + x_expected, + sx_expected ) @@ -316,13 +299,13 @@ def model_definition_piecewise_many_conditions(): events = {} # Analytical solution - def x_pected(t, x_1_0): + def x_expected(t, x_1_0): if np.floor(t) % 2 == 1: - return (x_1_0 + (np.floor(t)-1)/2 + (t-np.floor(t)), ) + return x_1_0 + (np.floor(t) - 1) / 2 + (t - np.floor(t)), else: - return (x_1_0 + np.floor(t)/2, ) + return x_1_0 + np.floor(t) / 2, - def sx_pected(t, x_1_0): + def sx_expected(t, x_1_0): return np.array([[1, ], ]) return ( @@ -332,6 +315,6 @@ def sx_pected(t, x_1_0): species, events, timepoints, - x_pected, - sx_pected + x_expected, + sx_expected ) diff --git a/python/tests/util.py b/python/tests/util.py index c51070a8b0..65c9fb77b3 100644 --- a/python/tests/util.py +++ b/python/tests/util.py @@ -1,7 +1,8 @@ """Tests for SBML events, including piecewise expressions.""" +from pathlib import Path + import libsbml import numpy as np -from pathlib import Path from amici import ( AmiciModel, @@ -26,13 +27,12 @@ def create_amici_model(sbml_model, model_name, **kwargs) -> AmiciModel: output_dir = sbml_test_models_output_dir / model_name sbml_importer.sbml2amici( model_name=model_name, - output_dir=str(output_dir), + output_dir=output_dir, **kwargs ) - model_module = import_model_module(model_name, str(output_dir.resolve())) - model = model_module.getModel() - return model + model_module = import_model_module(model_name, output_dir) + return model_module.getModel() def create_sbml_model( @@ -113,10 +113,7 @@ def creat_event_assignment(target, assignment): event_def['assignment']) if to_file: - libsbml.writeSBMLToFile( - document, - str(to_file), - ) + libsbml.writeSBMLToFile(document, to_file) # Need to return document, else AMICI throws an error. # (possibly due to garbage collection?) From f08b3f1ea10d3a9b1b65c5e9d008f7f0e54e4edb Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 12 May 2023 07:35:33 +0200 Subject: [PATCH 2/9] typo --- python/tests/util.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/tests/util.py b/python/tests/util.py index 65c9fb77b3..2731333e03 100644 --- a/python/tests/util.py +++ b/python/tests/util.py @@ -97,7 +97,7 @@ def create_sbml_model( trigger.setPersistent(True) trigger.setInitialValue(True) - def creat_event_assignment(target, assignment): + def create_event_assignment(target, assignment): ea = event.createEventAssignment() ea.setVariable(target) ea.setMath(libsbml.parseL3Formula(assignment)) @@ -106,11 +106,11 @@ def creat_event_assignment(target, assignment): for event_target, event_assignment in zip( event_def['target'], event_def['assignment'] ): - creat_event_assignment(event_target, event_assignment) + create_event_assignment(event_target, event_assignment) else: - creat_event_assignment(event_def['target'], - event_def['assignment']) + create_event_assignment(event_def['target'], + event_def['assignment']) if to_file: libsbml.writeSBMLToFile(document, to_file) From 05892374d766ea940313368838c95384fbd6d3f7 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 12 May 2023 13:51:40 +0200 Subject: [PATCH 3/9] windows --- python/tests/util.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/python/tests/util.py b/python/tests/util.py index 566bf50c41..95e93a8472 100644 --- a/python/tests/util.py +++ b/python/tests/util.py @@ -1,4 +1,6 @@ """Tests for SBML events, including piecewise expressions.""" +import sys +import tempfile from pathlib import Path import libsbml @@ -19,12 +21,17 @@ def create_amici_model(sbml_model, model_name, **kwargs) -> AmiciModel: """ Import an sbml file and create an AMICI model from it """ - sbml_test_models = Path("sbml_test_models") - sbml_test_models_output_dir = sbml_test_models / "amici_models" + sbml_test_models_output_dir = Path("amici_models") sbml_test_models_output_dir.mkdir(parents=True, exist_ok=True) sbml_importer = SbmlImporter(sbml_model) - output_dir = sbml_test_models_output_dir / model_name + # try not to exceed the stupid maximum path length on windows 💩 + output_dir = ( + sbml_test_models_output_dir / model_name + if sys.platform != "win32" + else tempfile.mkdtemp() + ) + sbml_importer.sbml2amici(model_name=model_name, output_dir=output_dir, **kwargs) model_module = import_model_module(model_name, output_dir) From 4381ac3ca9d969071f0389c575c24371e5eef317 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 15 May 2023 16:54:32 +0200 Subject: [PATCH 4/9] fix? --- python/sdist/amici/de_export.py | 221 +++++++++++++++++++++++--------- 1 file changed, 163 insertions(+), 58 deletions(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index 30f2cc99ba..bc924c6720 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -70,7 +70,9 @@ # Template for model/swig/CMakeLists.txt SWIG_CMAKE_TEMPLATE_FILE = os.path.join(amiciSwigPath, "CMakeLists_model.cmake") # Template for model/CMakeLists.txt -MODEL_CMAKE_TEMPLATE_FILE = os.path.join(amiciSrcPath, "CMakeLists.template.cmake") +MODEL_CMAKE_TEMPLATE_FILE = os.path.join( + amiciSrcPath, "CMakeLists.template.cmake" +) IDENTIFIER_PATTERN = re.compile(r"^[a-zA-Z_]\w*$") DERIVATIVE_PATTERN = re.compile(r"^d(x_rdata|xdot|\w+?)d(\w+?)(?:_explicit)?$") @@ -283,7 +285,8 @@ def arguments(self, ode: bool = True) -> str: " const realtype *k, const int ip", ), "sigmaz": _FunctionInfo( - "realtype *sigmaz, const realtype t, const realtype *p, " "const realtype *k", + "realtype *sigmaz, const realtype t, const realtype *p, " + "const realtype *k", ), "sroot": _FunctionInfo( "realtype *stau, const realtype t, const realtype *x, " @@ -324,7 +327,8 @@ def arguments(self, ode: bool = True) -> str: assume_pow_positivity=True, ), "x0": _FunctionInfo( - "realtype *x0, const realtype t, const realtype *p, " "const realtype *k" + "realtype *x0, const realtype t, const realtype *p, " + "const realtype *k" ), "x0_fixedParameters": _FunctionInfo( "realtype *x0_fixedParameters, const realtype t, " @@ -1016,10 +1020,14 @@ def transform_dxdt_to_concentration(species_id, dxdt): # we may end up with a time derivative of the compartment # volume due to parameter rate rules comp_rate_vars = [ - p for p in v.free_symbols if p in si.symbols[SymbolId.SPECIES] + p + for p in v.free_symbols + if p in si.symbols[SymbolId.SPECIES] ] for var in comp_rate_vars: - dv_dt += v.diff(var) * si.symbols[SymbolId.SPECIES][var]["dt"] + dv_dt += ( + v.diff(var) * si.symbols[SymbolId.SPECIES][var]["dt"] + ) dv_dx = v.diff(species_id) xdot = (dxdt - dv_dt * species_id) / (dv_dx * species_id + v) return xdot @@ -1038,7 +1046,9 @@ def transform_dxdt_to_concentration(species_id, dxdt): return dxdt / v # create dynamics without respecting conservation laws first - dxdt = smart_multiply(si.stoichiometric_matrix, MutableDenseMatrix(fluxes)) + dxdt = smart_multiply( + si.stoichiometric_matrix, MutableDenseMatrix(fluxes) + ) for ix, ((species_id, species), formula) in enumerate( zip(symbols[SymbolId.SPECIES].items(), dxdt) ): @@ -1048,7 +1058,9 @@ def transform_dxdt_to_concentration(species_id, dxdt): if species["amount"]: species["dt"] = formula else: - species["dt"] = transform_dxdt_to_concentration(species_id, formula) + species["dt"] = transform_dxdt_to_concentration( + species_id, formula + ) # create all basic components of the DE model and add them. for symbol_name in symbols: @@ -1096,7 +1108,8 @@ def transform_dxdt_to_concentration(species_id, dxdt): # fill in 'self._sym' based on prototypes and components in ode_model self.generate_basic_variables() self._has_quadratic_nllh = all( - llh["dist"] in ["normal", "lin-normal", "log-normal", "log10-normal"] + llh["dist"] + in ["normal", "lin-normal", "log-normal", "log10-normal"] for llh in si.symbols[SymbolId.LLHY].values() ) @@ -1178,14 +1191,21 @@ def add_conservation_law( )[0] except StopIteration: raise ValueError( - f"Specified state {state} was not found in the " f"model states." + f"Specified state {state} was not found in the " + f"model states." ) state_id = self._differential_states[ix].get_id() # \sum_{i≠j}(a_i * x_i)/a_j target_expression = ( - sp.Add(*(c_i * x_i for x_i, c_i in coefficients.items() if x_i != state)) + sp.Add( + *( + c_i * x_i + for x_i, c_i in coefficients.items() + if x_i != state + ) + ) / coefficients[state] ) @@ -1382,7 +1402,9 @@ def sparseeq(self, name) -> sp.Matrix: self._generate_sparse_symbol(name) return self._sparseeqs[name] - def colptrs(self, name: str) -> Union[List[sp.Number], List[List[sp.Number]]]: + def colptrs( + self, name: str + ) -> Union[List[sp.Number], List[List[sp.Number]]]: """ Returns (and constructs if necessary) the column pointers for a sparsified symbolic variable. @@ -1399,7 +1421,9 @@ def colptrs(self, name: str) -> Union[List[sp.Number], List[List[sp.Number]]]: self._generate_sparse_symbol(name) return self._colptrs[name] - def rowvals(self, name: str) -> Union[List[sp.Number], List[List[sp.Number]]]: + def rowvals( + self, name: str + ) -> Union[List[sp.Number], List[List[sp.Number]]]: """ Returns (and constructs if necessary) the row values for a sparsified symbolic variable. @@ -1651,7 +1675,9 @@ def get_appearance_counts(self, idxs: List[int]) -> List[int]: return [ free_symbols_dt.count(str(self._differential_states[idx].get_id())) - + free_symbols_expr.count(str(self._differential_states[idx].get_id())) + + free_symbols_expr.count( + str(self._differential_states[idx].get_id()) + ) for idx in idxs ] @@ -1736,7 +1762,9 @@ def _compute_equation(self, name: str) -> None: time_symbol = sp.Matrix([symbol_with_assumptions("t")]) if name in self._equation_prototype: - self._equation_from_components(name, self._equation_prototype[name]()) + self._equation_from_components( + name, self._equation_prototype[name]() + ) elif name in self._total_derivative_prototypes: args = self._total_derivative_prototypes[name] @@ -1826,7 +1854,9 @@ def _compute_equation(self, name: str) -> None: if any(sym in eq.free_symbols for sym in k) ] eq = self.eq("x0") - self._eqs[name] = sp.Matrix([eq[ix] for ix in self._x0_fixedParameters_idx]) + self._eqs[name] = sp.Matrix( + [eq[ix] for ix in self._x0_fixedParameters_idx] + ) elif name == "dtotal_cldx_rdata": x_rdata = self.sym("x_rdata") @@ -1839,7 +1869,9 @@ def _compute_equation(self, name: str) -> None: elif name == "dtcldx": # this is always zero - self._eqs[name] = sp.zeros(self.num_cons_law(), self.num_states_solver()) + self._eqs[name] = sp.zeros( + self.num_cons_law(), self.num_states_solver() + ) elif name == "dtcldp": # force symbols @@ -1864,15 +1896,21 @@ def _compute_equation(self, name: str) -> None: elif name == "dx_rdatadp": if self.num_cons_law(): - self._eqs[name] = smart_jacobian(self.eq("x_rdata"), self.sym("p")) + self._eqs[name] = smart_jacobian( + self.eq("x_rdata"), self.sym("p") + ) else: # so far, dx_rdatadp is only required for sx_rdata # in case of no conservation laws, C++ code will directly use # sx, we don't need this - self._eqs[name] = sp.zeros(self.num_states_rdata(), self.num_par()) + self._eqs[name] = sp.zeros( + self.num_states_rdata(), self.num_par() + ) elif name == "dx_rdatadtcl": - self._eqs[name] = smart_jacobian(self.eq("x_rdata"), self.sym("tcl")) + self._eqs[name] = smart_jacobian( + self.eq("x_rdata"), self.sym("tcl") + ) elif name == "dxdotdx_explicit": # force symbols @@ -1935,7 +1973,9 @@ def _compute_equation(self, name: str) -> None: self._eqs[name] = event_eqs elif name == "z": - event_observables = [sp.zeros(self.num_eventobs(), 1) for _ in self._events] + event_observables = [ + sp.zeros(self.num_eventobs(), 1) for _ in self._events + ] event_ids = [e.get_id() for e in self._events] # TODO: get rid of this stupid 1-based indexing as soon as we can # the matlab interface @@ -1943,7 +1983,9 @@ def _compute_equation(self, name: str) -> None: event_ids.index(event_obs.get_event()) + 1 for event_obs in self._event_observables ] - for (iz, ie), event_obs in zip(enumerate(z2event), self._event_observables): + for (iz, ie), event_obs in zip( + enumerate(z2event), self._event_observables + ): event_observables[ie - 1][iz] = event_obs.get_val() self._eqs[name] = event_observables @@ -1961,7 +2003,10 @@ def _compute_equation(self, name: str) -> None: ] if name == "dzdx": for ie in range(self.num_events()): - dtaudx = -self.eq("drootdx")[ie, :] / self.eq("drootdt_total")[ie] + dtaudx = ( + -self.eq("drootdx")[ie, :] + / self.eq("drootdt_total")[ie] + ) for iz in range(self.num_eventobs()): if ie != self._z2event[iz] - 1: continue @@ -2024,7 +2069,9 @@ def _compute_equation(self, name: str) -> None: # This part only works if we use self.eq('xdot') # instead of self.sym('xdot'). Not immediately clear # why that is. - tmp_dxdp += smart_multiply(self.eq("xdot"), self.sym("stau").T) + tmp_dxdp += smart_multiply( + self.eq("xdot_old"), self.sym("stau").T + ) # finish chain rule for the state variables tmp_eq += smart_multiply(self.eq("ddeltaxdx")[ie], tmp_dxdp) @@ -2046,7 +2093,9 @@ def _compute_equation(self, name: str) -> None: # that we need to reverse the order here for cl in reversed(self._conservation_laws) ] - ).col_join(smart_jacobian(self.eq("w")[self.num_cons_law() :, :], x)) + ).col_join( + smart_jacobian(self.eq("w")[self.num_cons_law() :, :], x) + ) elif match_deriv: self._derivative(match_deriv[1], match_deriv[2], name) @@ -2120,7 +2169,10 @@ def _derivative(self, eq: str, var: str, name: str = None) -> None: and cv not in self._lock_total_derivative and var != cv and min(self.sym(cv).shape) - and ((eq, var) not in ignore_chainrule or ignore_chainrule[(eq, var)] != cv) + and ( + (eq, var) not in ignore_chainrule + or ignore_chainrule[(eq, var)] != cv + ) ] if len(chainvars): self._lock_total_derivative += chainvars @@ -2223,9 +2275,14 @@ def _total_derivative( dxdz = self.sym_or_eq(name, dxdz_name) # Save time for large models if one multiplicand is zero, # which is not checked for by sympy - if not smart_is_zero_matrix(dydx) and not smart_is_zero_matrix(dxdz): + if not smart_is_zero_matrix(dydx) and not smart_is_zero_matrix( + dxdz + ): dydx_times_dxdz = smart_multiply(dydx, dxdz) - if dxdz.shape[1] == 1 and self._eqs[name].shape[1] != dxdz.shape[1]: + if ( + dxdz.shape[1] == 1 + and self._eqs[name].shape[1] != dxdz.shape[1] + ): for iz in range(self._eqs[name].shape[1]): self._eqs[name][:, iz] += dydx_times_dxdz else: @@ -2251,7 +2308,9 @@ def sym_or_eq(self, name: str, varname: str) -> sp.Matrix: # within a column may differ from the initialization of symbols here, # so those are not safe to use. Not removing them from signature as # this would break backwards compatibility. - if var_in_function_signature(name, varname, self.is_ode()) and varname not in [ + if var_in_function_signature( + name, varname, self.is_ode() + ) and varname not in [ "dwdx", "dwdp", ]: @@ -2375,7 +2434,8 @@ def state_has_fixed_parameter_initial_condition(self, ix: int) -> bool: if not isinstance(ic, sp.Basic): return False return any( - fp in (c.get_id() for c in self._constants) for fp in ic.free_symbols + fp in (c.get_id() for c in self._constants) + for fp in ic.free_symbols ) def state_has_conservation_law(self, ix: int) -> bool: @@ -2692,7 +2752,9 @@ def __init__( # include/amici/model.h for details) self.model: DEModel = de_model self.model._code_printer.known_functions.update( - splines.spline_user_functions(self.model.splines, self._get_index("p")) + splines.spline_user_functions( + self.model.splines, self._get_index("p") + ) ) # To only generate a subset of functions, apply subselection here @@ -2708,7 +2770,9 @@ def generate_model_code(self) -> None: Generates the native C++ code for the loaded model and a Matlab script that can be run to compile a mex file from the C++ code """ - with _monkeypatched(sp.Pow, "_eval_derivative", _custom_pow_eval_derivative): + with _monkeypatched( + sp.Pow, "_eval_derivative", _custom_pow_eval_derivative + ): self._prepare_model_folder() self._generate_c_code() self._generate_m_code() @@ -2759,7 +2823,9 @@ def _generate_c_code(self) -> None: name in self.functions and not self.functions[name].body and name not in nobody_functions - ) or (name not in self.functions and len(self.model.sym(name)) == 0): + ) or ( + name not in self.functions and len(self.model.sym(name)) == 0 + ): continue self._write_index_files(name) @@ -2770,7 +2836,9 @@ def _generate_c_code(self) -> None: self._write_swig_files() self._write_module_setup() - shutil.copy(CXX_MAIN_TEMPLATE_FILE, os.path.join(self.model_path, "main.cpp")) + shutil.copy( + CXX_MAIN_TEMPLATE_FILE, os.path.join(self.model_path, "main.cpp") + ) def _compile_c_code( self, @@ -2848,8 +2916,10 @@ def _generate_m_code(self) -> None: lines = [ "% This compile script was automatically created from" " Python SBML import.", - "% If mex compiler is set up within MATLAB, it can be run" " from MATLAB ", - "% in order to compile a mex-file from the Python" " generated C++ files.", + "% If mex compiler is set up within MATLAB, it can be run" + " from MATLAB ", + "% in order to compile a mex-file from the Python" + " generated C++ files.", "", f"modelName = '{self.model_name}';", "amimodel.compileAndLinkModel(modelName, '', [], [], [], []);", @@ -2881,7 +2951,10 @@ def _get_index(self, name: str) -> Dict[sp.Symbol, int]: else: raise ValueError(f"Unknown symbolic array: {name}") - return {strip_pysb(symbol).name: index for index, symbol in enumerate(symbols)} + return { + strip_pysb(symbol).name: index + for index, symbol in enumerate(symbols) + } def _write_index_files(self, name: str) -> None: """ @@ -2934,7 +3007,8 @@ def _write_function_file(self, function: str) -> None: if function in sparse_functions: equations = self.model.sparseeq(function) elif ( - not self.allow_reinit_fixpar_initcond and function == "sx0_fixedParameters" + not self.allow_reinit_fixpar_initcond + and function == "sx0_fixedParameters" ): # Not required. Will create empty function body. equations = sp.Matrix() @@ -2983,7 +3057,10 @@ def _write_function_file(self, function: str) -> None: lines.append(f'#include "{sym}.h"') # include return symbols - if function in self.model.sym_names() and function not in non_unique_id_symbols: + if ( + function in self.model.sym_names() + and function not in non_unique_id_symbols + ): lines.append(f'#include "{function}.h"') lines.extend( @@ -3010,7 +3087,9 @@ def _write_function_file(self, function: str) -> None: body = [ # execute this twice to catch cases where the ending '(' would # be the starting (^|\W) for the following match - pow_rx.sub(r"\1amici::pos_pow(", pow_rx.sub(r"\1amici::pos_pow(", line)) + pow_rx.sub( + r"\1amici::pos_pow(", pow_rx.sub(r"\1amici::pos_pow(", line) + ) for line in body ] @@ -3140,7 +3219,9 @@ def _write_function_index(self, function: str, indextype: str) -> None: with open(filename, "w") as fileout: fileout.write("\n".join(lines)) - def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: + def _get_function_body( + self, function: str, equations: sp.Matrix + ) -> List[str]: """ Generate C++ code for body of function ``function``. @@ -3179,7 +3260,9 @@ def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: + str(len(self.model._x0_fixedParameters_idx)) + "> _x0_fixedParameters_idxs = {", " " - + ", ".join(str(x) for x in self.model._x0_fixedParameters_idx), + + ", ".join( + str(x) for x in self.model._x0_fixedParameters_idx + ), " };", "", # Set all parameters that are to be reset to 0, so that the @@ -3216,7 +3299,9 @@ def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: lines.extend(get_switch_statement("ip", cases, 1)) elif function == "x0_fixedParameters": - for index, formula in zip(self.model._x0_fixedParameters_idx, equations): + for index, formula in zip( + self.model._x0_fixedParameters_idx, equations + ): lines.append( f" if(std::find(reinitialization_state_idxs.cbegin(), " f"reinitialization_state_idxs.cend(), {index}) != " @@ -3250,7 +3335,10 @@ def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: outer_cases[ie] = copy.copy(inner_lines) lines.extend(get_switch_statement("ie", outer_cases, 1)) - elif function in sensi_functions and equations.shape[1] == self.model.num_par(): + elif ( + function in sensi_functions + and equations.shape[1] == self.model.num_par() + ): cases = { ipar: self.model._code_printer._get_sym_lines_array( equations[:, ipar], function, 0 @@ -3283,7 +3371,8 @@ def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: lines.extend(get_switch_statement(iterator, cases, 1)) elif ( - function in self.model.sym_names() and function not in non_unique_id_symbols + function in self.model.sym_names() + and function not in non_unique_id_symbols ): if function in sparse_functions: symbols = list(map(sp.Symbol, self.model.sparsesym(function))) @@ -3315,14 +3404,14 @@ def _get_create_splines_body(self): nodes = f"{ind8}{{{', '.join(map(str, spline.nodes))}}}, " # vector with the node values - values = f"{ind8}{{{', '.join(map(str, spline.values_at_nodes))}}}, " + values = ( + f"{ind8}{{{', '.join(map(str, spline.values_at_nodes))}}}, " + ) # vector with the slopes if spline.derivatives_by_fd: slopes = f"{ind8}{{}}," else: - slopes = ( - f"{ind8}{{{', '.join(map(str, spline.derivatives_at_nodes))}}}," - ) + slopes = f"{ind8}{{{', '.join(map(str, spline.derivatives_at_nodes))}}}," body.extend( [ @@ -3345,7 +3434,8 @@ def _get_create_splines_body(self): body.append(ind8 + bc_to_cpp[bc]) except KeyError: raise ValueError( - f"Unknown boundary condition '{bc}' " "found in spline object" + f"Unknown boundary condition '{bc}' " + "found in spline object" ) extrapolate_to_cpp = { None: "SplineExtrapolation::noExtrapolation, ", @@ -3359,12 +3449,15 @@ def _get_create_splines_body(self): body.append(ind8 + extrapolate_to_cpp[extr]) except KeyError: raise ValueError( - f"Unknown extrapolation '{extr}' " "found in spline object" + f"Unknown extrapolation '{extr}' " + "found in spline object" ) line = ind8 line += "true, " if spline.derivatives_by_fd else "false, " line += ( - "true, " if isinstance(spline.nodes, splines.UniformGrid) else "false, " + "true, " + if isinstance(spline.nodes, splines.UniformGrid) + else "false, " ) line += "true" if spline.logarithmic_parametrization else "false" body.append(line) @@ -3443,10 +3536,12 @@ def _write_model_header_cpp(self) -> None: "NK": self.model.num_const(), "O2MODE": "amici::SecondOrderMode::none", # using code printer ensures proper handling of nan/inf - "PARAMETERS": self.model._code_printer.doprint(self.model.val("p"))[1:-1], - "FIXED_PARAMETERS": self.model._code_printer.doprint(self.model.val("k"))[ + "PARAMETERS": self.model._code_printer.doprint(self.model.val("p"))[ 1:-1 ], + "FIXED_PARAMETERS": self.model._code_printer.doprint( + self.model.val("k") + )[1:-1], "PARAMETER_NAMES_INITIALIZER_LIST": self._get_symbol_name_initializer_list( "p" ), @@ -3461,12 +3556,16 @@ def _write_model_header_cpp(self) -> None: ), "OBSERVABLE_TRAFO_INITIALIZER_LIST": "\n".join( f"ObservableScaling::{trafo.value}, // y[{idx}]" - for idx, trafo in enumerate(self.model.get_observable_transformations()) + for idx, trafo in enumerate( + self.model.get_observable_transformations() + ) ), "EXPRESSION_NAMES_INITIALIZER_LIST": self._get_symbol_name_initializer_list( "w" ), - "PARAMETER_IDS_INITIALIZER_LIST": self._get_symbol_id_initializer_list("p"), + "PARAMETER_IDS_INITIALIZER_LIST": self._get_symbol_id_initializer_list( + "p" + ), "STATE_IDS_INITIALIZER_LIST": self._get_symbol_id_initializer_list( "x_rdata" ), @@ -3547,16 +3646,22 @@ def _write_model_header_cpp(self) -> None: indexfield, nobody=True, ) - tpl_data[f"{func_name.upper()}_{indexfield.upper()}_DEF"] = "" + tpl_data[ + f"{func_name.upper()}_{indexfield.upper()}_DEF" + ] = "" tpl_data[ f"{func_name.upper()}_{indexfield.upper()}_IMPL" ] = impl continue - tpl_data[f"{func_name.upper()}_DEF"] = get_function_extern_declaration( + tpl_data[ + f"{func_name.upper()}_DEF" + ] = get_function_extern_declaration( func_name, self.model_name, self.model.is_ode() ) - tpl_data[f"{func_name.upper()}_IMPL"] = get_model_override_implementation( + tpl_data[ + f"{func_name.upper()}_IMPL" + ] = get_model_override_implementation( func_name, self.model_name, self.model.is_ode() ) if func_name in sparse_functions: From d9d41bdc6dd98800880603e7f095c6f53134412e Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 15 May 2023 16:58:51 +0200 Subject: [PATCH 5/9] Revert "fix?" This reverts commit 4381ac3ca9d969071f0389c575c24371e5eef317. --- python/sdist/amici/de_export.py | 221 +++++++++----------------------- 1 file changed, 58 insertions(+), 163 deletions(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index bc924c6720..30f2cc99ba 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -70,9 +70,7 @@ # Template for model/swig/CMakeLists.txt SWIG_CMAKE_TEMPLATE_FILE = os.path.join(amiciSwigPath, "CMakeLists_model.cmake") # Template for model/CMakeLists.txt -MODEL_CMAKE_TEMPLATE_FILE = os.path.join( - amiciSrcPath, "CMakeLists.template.cmake" -) +MODEL_CMAKE_TEMPLATE_FILE = os.path.join(amiciSrcPath, "CMakeLists.template.cmake") IDENTIFIER_PATTERN = re.compile(r"^[a-zA-Z_]\w*$") DERIVATIVE_PATTERN = re.compile(r"^d(x_rdata|xdot|\w+?)d(\w+?)(?:_explicit)?$") @@ -285,8 +283,7 @@ def arguments(self, ode: bool = True) -> str: " const realtype *k, const int ip", ), "sigmaz": _FunctionInfo( - "realtype *sigmaz, const realtype t, const realtype *p, " - "const realtype *k", + "realtype *sigmaz, const realtype t, const realtype *p, " "const realtype *k", ), "sroot": _FunctionInfo( "realtype *stau, const realtype t, const realtype *x, " @@ -327,8 +324,7 @@ def arguments(self, ode: bool = True) -> str: assume_pow_positivity=True, ), "x0": _FunctionInfo( - "realtype *x0, const realtype t, const realtype *p, " - "const realtype *k" + "realtype *x0, const realtype t, const realtype *p, " "const realtype *k" ), "x0_fixedParameters": _FunctionInfo( "realtype *x0_fixedParameters, const realtype t, " @@ -1020,14 +1016,10 @@ def transform_dxdt_to_concentration(species_id, dxdt): # we may end up with a time derivative of the compartment # volume due to parameter rate rules comp_rate_vars = [ - p - for p in v.free_symbols - if p in si.symbols[SymbolId.SPECIES] + p for p in v.free_symbols if p in si.symbols[SymbolId.SPECIES] ] for var in comp_rate_vars: - dv_dt += ( - v.diff(var) * si.symbols[SymbolId.SPECIES][var]["dt"] - ) + dv_dt += v.diff(var) * si.symbols[SymbolId.SPECIES][var]["dt"] dv_dx = v.diff(species_id) xdot = (dxdt - dv_dt * species_id) / (dv_dx * species_id + v) return xdot @@ -1046,9 +1038,7 @@ def transform_dxdt_to_concentration(species_id, dxdt): return dxdt / v # create dynamics without respecting conservation laws first - dxdt = smart_multiply( - si.stoichiometric_matrix, MutableDenseMatrix(fluxes) - ) + dxdt = smart_multiply(si.stoichiometric_matrix, MutableDenseMatrix(fluxes)) for ix, ((species_id, species), formula) in enumerate( zip(symbols[SymbolId.SPECIES].items(), dxdt) ): @@ -1058,9 +1048,7 @@ def transform_dxdt_to_concentration(species_id, dxdt): if species["amount"]: species["dt"] = formula else: - species["dt"] = transform_dxdt_to_concentration( - species_id, formula - ) + species["dt"] = transform_dxdt_to_concentration(species_id, formula) # create all basic components of the DE model and add them. for symbol_name in symbols: @@ -1108,8 +1096,7 @@ def transform_dxdt_to_concentration(species_id, dxdt): # fill in 'self._sym' based on prototypes and components in ode_model self.generate_basic_variables() self._has_quadratic_nllh = all( - llh["dist"] - in ["normal", "lin-normal", "log-normal", "log10-normal"] + llh["dist"] in ["normal", "lin-normal", "log-normal", "log10-normal"] for llh in si.symbols[SymbolId.LLHY].values() ) @@ -1191,21 +1178,14 @@ def add_conservation_law( )[0] except StopIteration: raise ValueError( - f"Specified state {state} was not found in the " - f"model states." + f"Specified state {state} was not found in the " f"model states." ) state_id = self._differential_states[ix].get_id() # \sum_{i≠j}(a_i * x_i)/a_j target_expression = ( - sp.Add( - *( - c_i * x_i - for x_i, c_i in coefficients.items() - if x_i != state - ) - ) + sp.Add(*(c_i * x_i for x_i, c_i in coefficients.items() if x_i != state)) / coefficients[state] ) @@ -1402,9 +1382,7 @@ def sparseeq(self, name) -> sp.Matrix: self._generate_sparse_symbol(name) return self._sparseeqs[name] - def colptrs( - self, name: str - ) -> Union[List[sp.Number], List[List[sp.Number]]]: + def colptrs(self, name: str) -> Union[List[sp.Number], List[List[sp.Number]]]: """ Returns (and constructs if necessary) the column pointers for a sparsified symbolic variable. @@ -1421,9 +1399,7 @@ def colptrs( self._generate_sparse_symbol(name) return self._colptrs[name] - def rowvals( - self, name: str - ) -> Union[List[sp.Number], List[List[sp.Number]]]: + def rowvals(self, name: str) -> Union[List[sp.Number], List[List[sp.Number]]]: """ Returns (and constructs if necessary) the row values for a sparsified symbolic variable. @@ -1675,9 +1651,7 @@ def get_appearance_counts(self, idxs: List[int]) -> List[int]: return [ free_symbols_dt.count(str(self._differential_states[idx].get_id())) - + free_symbols_expr.count( - str(self._differential_states[idx].get_id()) - ) + + free_symbols_expr.count(str(self._differential_states[idx].get_id())) for idx in idxs ] @@ -1762,9 +1736,7 @@ def _compute_equation(self, name: str) -> None: time_symbol = sp.Matrix([symbol_with_assumptions("t")]) if name in self._equation_prototype: - self._equation_from_components( - name, self._equation_prototype[name]() - ) + self._equation_from_components(name, self._equation_prototype[name]()) elif name in self._total_derivative_prototypes: args = self._total_derivative_prototypes[name] @@ -1854,9 +1826,7 @@ def _compute_equation(self, name: str) -> None: if any(sym in eq.free_symbols for sym in k) ] eq = self.eq("x0") - self._eqs[name] = sp.Matrix( - [eq[ix] for ix in self._x0_fixedParameters_idx] - ) + self._eqs[name] = sp.Matrix([eq[ix] for ix in self._x0_fixedParameters_idx]) elif name == "dtotal_cldx_rdata": x_rdata = self.sym("x_rdata") @@ -1869,9 +1839,7 @@ def _compute_equation(self, name: str) -> None: elif name == "dtcldx": # this is always zero - self._eqs[name] = sp.zeros( - self.num_cons_law(), self.num_states_solver() - ) + self._eqs[name] = sp.zeros(self.num_cons_law(), self.num_states_solver()) elif name == "dtcldp": # force symbols @@ -1896,21 +1864,15 @@ def _compute_equation(self, name: str) -> None: elif name == "dx_rdatadp": if self.num_cons_law(): - self._eqs[name] = smart_jacobian( - self.eq("x_rdata"), self.sym("p") - ) + self._eqs[name] = smart_jacobian(self.eq("x_rdata"), self.sym("p")) else: # so far, dx_rdatadp is only required for sx_rdata # in case of no conservation laws, C++ code will directly use # sx, we don't need this - self._eqs[name] = sp.zeros( - self.num_states_rdata(), self.num_par() - ) + self._eqs[name] = sp.zeros(self.num_states_rdata(), self.num_par()) elif name == "dx_rdatadtcl": - self._eqs[name] = smart_jacobian( - self.eq("x_rdata"), self.sym("tcl") - ) + self._eqs[name] = smart_jacobian(self.eq("x_rdata"), self.sym("tcl")) elif name == "dxdotdx_explicit": # force symbols @@ -1973,9 +1935,7 @@ def _compute_equation(self, name: str) -> None: self._eqs[name] = event_eqs elif name == "z": - event_observables = [ - sp.zeros(self.num_eventobs(), 1) for _ in self._events - ] + event_observables = [sp.zeros(self.num_eventobs(), 1) for _ in self._events] event_ids = [e.get_id() for e in self._events] # TODO: get rid of this stupid 1-based indexing as soon as we can # the matlab interface @@ -1983,9 +1943,7 @@ def _compute_equation(self, name: str) -> None: event_ids.index(event_obs.get_event()) + 1 for event_obs in self._event_observables ] - for (iz, ie), event_obs in zip( - enumerate(z2event), self._event_observables - ): + for (iz, ie), event_obs in zip(enumerate(z2event), self._event_observables): event_observables[ie - 1][iz] = event_obs.get_val() self._eqs[name] = event_observables @@ -2003,10 +1961,7 @@ def _compute_equation(self, name: str) -> None: ] if name == "dzdx": for ie in range(self.num_events()): - dtaudx = ( - -self.eq("drootdx")[ie, :] - / self.eq("drootdt_total")[ie] - ) + dtaudx = -self.eq("drootdx")[ie, :] / self.eq("drootdt_total")[ie] for iz in range(self.num_eventobs()): if ie != self._z2event[iz] - 1: continue @@ -2069,9 +2024,7 @@ def _compute_equation(self, name: str) -> None: # This part only works if we use self.eq('xdot') # instead of self.sym('xdot'). Not immediately clear # why that is. - tmp_dxdp += smart_multiply( - self.eq("xdot_old"), self.sym("stau").T - ) + tmp_dxdp += smart_multiply(self.eq("xdot"), self.sym("stau").T) # finish chain rule for the state variables tmp_eq += smart_multiply(self.eq("ddeltaxdx")[ie], tmp_dxdp) @@ -2093,9 +2046,7 @@ def _compute_equation(self, name: str) -> None: # that we need to reverse the order here for cl in reversed(self._conservation_laws) ] - ).col_join( - smart_jacobian(self.eq("w")[self.num_cons_law() :, :], x) - ) + ).col_join(smart_jacobian(self.eq("w")[self.num_cons_law() :, :], x)) elif match_deriv: self._derivative(match_deriv[1], match_deriv[2], name) @@ -2169,10 +2120,7 @@ def _derivative(self, eq: str, var: str, name: str = None) -> None: and cv not in self._lock_total_derivative and var != cv and min(self.sym(cv).shape) - and ( - (eq, var) not in ignore_chainrule - or ignore_chainrule[(eq, var)] != cv - ) + and ((eq, var) not in ignore_chainrule or ignore_chainrule[(eq, var)] != cv) ] if len(chainvars): self._lock_total_derivative += chainvars @@ -2275,14 +2223,9 @@ def _total_derivative( dxdz = self.sym_or_eq(name, dxdz_name) # Save time for large models if one multiplicand is zero, # which is not checked for by sympy - if not smart_is_zero_matrix(dydx) and not smart_is_zero_matrix( - dxdz - ): + if not smart_is_zero_matrix(dydx) and not smart_is_zero_matrix(dxdz): dydx_times_dxdz = smart_multiply(dydx, dxdz) - if ( - dxdz.shape[1] == 1 - and self._eqs[name].shape[1] != dxdz.shape[1] - ): + if dxdz.shape[1] == 1 and self._eqs[name].shape[1] != dxdz.shape[1]: for iz in range(self._eqs[name].shape[1]): self._eqs[name][:, iz] += dydx_times_dxdz else: @@ -2308,9 +2251,7 @@ def sym_or_eq(self, name: str, varname: str) -> sp.Matrix: # within a column may differ from the initialization of symbols here, # so those are not safe to use. Not removing them from signature as # this would break backwards compatibility. - if var_in_function_signature( - name, varname, self.is_ode() - ) and varname not in [ + if var_in_function_signature(name, varname, self.is_ode()) and varname not in [ "dwdx", "dwdp", ]: @@ -2434,8 +2375,7 @@ def state_has_fixed_parameter_initial_condition(self, ix: int) -> bool: if not isinstance(ic, sp.Basic): return False return any( - fp in (c.get_id() for c in self._constants) - for fp in ic.free_symbols + fp in (c.get_id() for c in self._constants) for fp in ic.free_symbols ) def state_has_conservation_law(self, ix: int) -> bool: @@ -2752,9 +2692,7 @@ def __init__( # include/amici/model.h for details) self.model: DEModel = de_model self.model._code_printer.known_functions.update( - splines.spline_user_functions( - self.model.splines, self._get_index("p") - ) + splines.spline_user_functions(self.model.splines, self._get_index("p")) ) # To only generate a subset of functions, apply subselection here @@ -2770,9 +2708,7 @@ def generate_model_code(self) -> None: Generates the native C++ code for the loaded model and a Matlab script that can be run to compile a mex file from the C++ code """ - with _monkeypatched( - sp.Pow, "_eval_derivative", _custom_pow_eval_derivative - ): + with _monkeypatched(sp.Pow, "_eval_derivative", _custom_pow_eval_derivative): self._prepare_model_folder() self._generate_c_code() self._generate_m_code() @@ -2823,9 +2759,7 @@ def _generate_c_code(self) -> None: name in self.functions and not self.functions[name].body and name not in nobody_functions - ) or ( - name not in self.functions and len(self.model.sym(name)) == 0 - ): + ) or (name not in self.functions and len(self.model.sym(name)) == 0): continue self._write_index_files(name) @@ -2836,9 +2770,7 @@ def _generate_c_code(self) -> None: self._write_swig_files() self._write_module_setup() - shutil.copy( - CXX_MAIN_TEMPLATE_FILE, os.path.join(self.model_path, "main.cpp") - ) + shutil.copy(CXX_MAIN_TEMPLATE_FILE, os.path.join(self.model_path, "main.cpp")) def _compile_c_code( self, @@ -2916,10 +2848,8 @@ def _generate_m_code(self) -> None: lines = [ "% This compile script was automatically created from" " Python SBML import.", - "% If mex compiler is set up within MATLAB, it can be run" - " from MATLAB ", - "% in order to compile a mex-file from the Python" - " generated C++ files.", + "% If mex compiler is set up within MATLAB, it can be run" " from MATLAB ", + "% in order to compile a mex-file from the Python" " generated C++ files.", "", f"modelName = '{self.model_name}';", "amimodel.compileAndLinkModel(modelName, '', [], [], [], []);", @@ -2951,10 +2881,7 @@ def _get_index(self, name: str) -> Dict[sp.Symbol, int]: else: raise ValueError(f"Unknown symbolic array: {name}") - return { - strip_pysb(symbol).name: index - for index, symbol in enumerate(symbols) - } + return {strip_pysb(symbol).name: index for index, symbol in enumerate(symbols)} def _write_index_files(self, name: str) -> None: """ @@ -3007,8 +2934,7 @@ def _write_function_file(self, function: str) -> None: if function in sparse_functions: equations = self.model.sparseeq(function) elif ( - not self.allow_reinit_fixpar_initcond - and function == "sx0_fixedParameters" + not self.allow_reinit_fixpar_initcond and function == "sx0_fixedParameters" ): # Not required. Will create empty function body. equations = sp.Matrix() @@ -3057,10 +2983,7 @@ def _write_function_file(self, function: str) -> None: lines.append(f'#include "{sym}.h"') # include return symbols - if ( - function in self.model.sym_names() - and function not in non_unique_id_symbols - ): + if function in self.model.sym_names() and function not in non_unique_id_symbols: lines.append(f'#include "{function}.h"') lines.extend( @@ -3087,9 +3010,7 @@ def _write_function_file(self, function: str) -> None: body = [ # execute this twice to catch cases where the ending '(' would # be the starting (^|\W) for the following match - pow_rx.sub( - r"\1amici::pos_pow(", pow_rx.sub(r"\1amici::pos_pow(", line) - ) + pow_rx.sub(r"\1amici::pos_pow(", pow_rx.sub(r"\1amici::pos_pow(", line)) for line in body ] @@ -3219,9 +3140,7 @@ def _write_function_index(self, function: str, indextype: str) -> None: with open(filename, "w") as fileout: fileout.write("\n".join(lines)) - def _get_function_body( - self, function: str, equations: sp.Matrix - ) -> List[str]: + def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: """ Generate C++ code for body of function ``function``. @@ -3260,9 +3179,7 @@ def _get_function_body( + str(len(self.model._x0_fixedParameters_idx)) + "> _x0_fixedParameters_idxs = {", " " - + ", ".join( - str(x) for x in self.model._x0_fixedParameters_idx - ), + + ", ".join(str(x) for x in self.model._x0_fixedParameters_idx), " };", "", # Set all parameters that are to be reset to 0, so that the @@ -3299,9 +3216,7 @@ def _get_function_body( lines.extend(get_switch_statement("ip", cases, 1)) elif function == "x0_fixedParameters": - for index, formula in zip( - self.model._x0_fixedParameters_idx, equations - ): + for index, formula in zip(self.model._x0_fixedParameters_idx, equations): lines.append( f" if(std::find(reinitialization_state_idxs.cbegin(), " f"reinitialization_state_idxs.cend(), {index}) != " @@ -3335,10 +3250,7 @@ def _get_function_body( outer_cases[ie] = copy.copy(inner_lines) lines.extend(get_switch_statement("ie", outer_cases, 1)) - elif ( - function in sensi_functions - and equations.shape[1] == self.model.num_par() - ): + elif function in sensi_functions and equations.shape[1] == self.model.num_par(): cases = { ipar: self.model._code_printer._get_sym_lines_array( equations[:, ipar], function, 0 @@ -3371,8 +3283,7 @@ def _get_function_body( lines.extend(get_switch_statement(iterator, cases, 1)) elif ( - function in self.model.sym_names() - and function not in non_unique_id_symbols + function in self.model.sym_names() and function not in non_unique_id_symbols ): if function in sparse_functions: symbols = list(map(sp.Symbol, self.model.sparsesym(function))) @@ -3404,14 +3315,14 @@ def _get_create_splines_body(self): nodes = f"{ind8}{{{', '.join(map(str, spline.nodes))}}}, " # vector with the node values - values = ( - f"{ind8}{{{', '.join(map(str, spline.values_at_nodes))}}}, " - ) + values = f"{ind8}{{{', '.join(map(str, spline.values_at_nodes))}}}, " # vector with the slopes if spline.derivatives_by_fd: slopes = f"{ind8}{{}}," else: - slopes = f"{ind8}{{{', '.join(map(str, spline.derivatives_at_nodes))}}}," + slopes = ( + f"{ind8}{{{', '.join(map(str, spline.derivatives_at_nodes))}}}," + ) body.extend( [ @@ -3434,8 +3345,7 @@ def _get_create_splines_body(self): body.append(ind8 + bc_to_cpp[bc]) except KeyError: raise ValueError( - f"Unknown boundary condition '{bc}' " - "found in spline object" + f"Unknown boundary condition '{bc}' " "found in spline object" ) extrapolate_to_cpp = { None: "SplineExtrapolation::noExtrapolation, ", @@ -3449,15 +3359,12 @@ def _get_create_splines_body(self): body.append(ind8 + extrapolate_to_cpp[extr]) except KeyError: raise ValueError( - f"Unknown extrapolation '{extr}' " - "found in spline object" + f"Unknown extrapolation '{extr}' " "found in spline object" ) line = ind8 line += "true, " if spline.derivatives_by_fd else "false, " line += ( - "true, " - if isinstance(spline.nodes, splines.UniformGrid) - else "false, " + "true, " if isinstance(spline.nodes, splines.UniformGrid) else "false, " ) line += "true" if spline.logarithmic_parametrization else "false" body.append(line) @@ -3536,12 +3443,10 @@ def _write_model_header_cpp(self) -> None: "NK": self.model.num_const(), "O2MODE": "amici::SecondOrderMode::none", # using code printer ensures proper handling of nan/inf - "PARAMETERS": self.model._code_printer.doprint(self.model.val("p"))[ + "PARAMETERS": self.model._code_printer.doprint(self.model.val("p"))[1:-1], + "FIXED_PARAMETERS": self.model._code_printer.doprint(self.model.val("k"))[ 1:-1 ], - "FIXED_PARAMETERS": self.model._code_printer.doprint( - self.model.val("k") - )[1:-1], "PARAMETER_NAMES_INITIALIZER_LIST": self._get_symbol_name_initializer_list( "p" ), @@ -3556,16 +3461,12 @@ def _write_model_header_cpp(self) -> None: ), "OBSERVABLE_TRAFO_INITIALIZER_LIST": "\n".join( f"ObservableScaling::{trafo.value}, // y[{idx}]" - for idx, trafo in enumerate( - self.model.get_observable_transformations() - ) + for idx, trafo in enumerate(self.model.get_observable_transformations()) ), "EXPRESSION_NAMES_INITIALIZER_LIST": self._get_symbol_name_initializer_list( "w" ), - "PARAMETER_IDS_INITIALIZER_LIST": self._get_symbol_id_initializer_list( - "p" - ), + "PARAMETER_IDS_INITIALIZER_LIST": self._get_symbol_id_initializer_list("p"), "STATE_IDS_INITIALIZER_LIST": self._get_symbol_id_initializer_list( "x_rdata" ), @@ -3646,22 +3547,16 @@ def _write_model_header_cpp(self) -> None: indexfield, nobody=True, ) - tpl_data[ - f"{func_name.upper()}_{indexfield.upper()}_DEF" - ] = "" + tpl_data[f"{func_name.upper()}_{indexfield.upper()}_DEF"] = "" tpl_data[ f"{func_name.upper()}_{indexfield.upper()}_IMPL" ] = impl continue - tpl_data[ - f"{func_name.upper()}_DEF" - ] = get_function_extern_declaration( + tpl_data[f"{func_name.upper()}_DEF"] = get_function_extern_declaration( func_name, self.model_name, self.model.is_ode() ) - tpl_data[ - f"{func_name.upper()}_IMPL" - ] = get_model_override_implementation( + tpl_data[f"{func_name.upper()}_IMPL"] = get_model_override_implementation( func_name, self.model_name, self.model.is_ode() ) if func_name in sparse_functions: From 06b5d7bca54cb18d89f53a97b3ab45acf91eb6f9 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 15 May 2023 16:59:24 +0200 Subject: [PATCH 6/9] fix? --- python/sdist/amici/de_export.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index 30f2cc99ba..151315ac1c 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -2024,7 +2024,7 @@ def _compute_equation(self, name: str) -> None: # This part only works if we use self.eq('xdot') # instead of self.sym('xdot'). Not immediately clear # why that is. - tmp_dxdp += smart_multiply(self.eq("xdot"), self.sym("stau").T) + tmp_dxdp += smart_multiply(self.eq("xdot_old"), self.sym("stau").T) # finish chain rule for the state variables tmp_eq += smart_multiply(self.eq("ddeltaxdx")[ie], tmp_dxdp) From 97baddc700fdd169e608a884d670b7d48441f9c2 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 15 May 2023 18:17:17 +0200 Subject: [PATCH 7/9] sym vs eq? --- python/sdist/amici/de_export.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index 151315ac1c..87af013964 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -2021,10 +2021,7 @@ def _compute_equation(self, name: str) -> None: ) # additional part of chain rule state variables - # This part only works if we use self.eq('xdot') - # instead of self.sym('xdot'). Not immediately clear - # why that is. - tmp_dxdp += smart_multiply(self.eq("xdot_old"), self.sym("stau").T) + tmp_dxdp += smart_multiply(self.sym("xdot_old"), self.sym("stau").T) # finish chain rule for the state variables tmp_eq += smart_multiply(self.eq("ddeltaxdx")[ie], tmp_dxdp) From 4ab881330615bc1412d853ff0cbe7f5e03a27fb7 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 16 May 2023 07:12:10 +0200 Subject: [PATCH 8/9] Why test twice? --- python/tests/util.py | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/python/tests/util.py b/python/tests/util.py index aa29cce88d..d5907da29c 100644 --- a/python/tests/util.py +++ b/python/tests/util.py @@ -129,14 +129,6 @@ def check_trajectories_without_sensitivities( Check whether the AMICI simulation matches a known solution (ideally an analytically calculated one). """ - - # Does the AMICI simulation match the analytical solution? - solver = amici_model.getSolver() - solver.setAbsoluteTolerance(1e-15) - rdata = runAmiciSimulation(amici_model, solver=solver) - _check_close(rdata["x"], result_expected_x, field="x", rtol=5e-5, atol=1e-13) - - # Show that we can do arbitrary precision here (test 8 digits) solver = amici_model.getSolver() solver.setAbsoluteTolerance(1e-15) solver.setRelativeTolerance(1e-12) @@ -153,17 +145,6 @@ def check_trajectories_with_forward_sensitivities( Check whether the forward sensitivities of the AMICI simulation match a known solution (ideally an analytically calculated one). """ - - # Show that we can do arbitrary precision here (test 8 digits) - solver = amici_model.getSolver() - solver.setAbsoluteTolerance(1e-15) - solver.setSensitivityOrder(SensitivityOrder.first) - solver.setSensitivityMethod(SensitivityMethod.forward) - rdata = runAmiciSimulation(amici_model, solver=solver) - _check_close(rdata["x"], result_expected_x, field="x", rtol=1e-5, atol=1e-13) - _check_close(rdata["sx"], result_expected_sx, field="sx", rtol=1e-5, atol=1e-7) - - # Show that we can do arbitrary precision here (test 8 digits) solver = amici_model.getSolver() solver.setSensitivityOrder(SensitivityOrder.first) solver.setSensitivityMethod(SensitivityMethod.forward) From f9c197b3a7b7853036c27a68fb748c0c10cefea1 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 16 May 2023 12:05:03 +0200 Subject: [PATCH 9/9] increase tolerance --- python/tests/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/util.py b/python/tests/util.py index d5907da29c..14f514c997 100644 --- a/python/tests/util.py +++ b/python/tests/util.py @@ -154,4 +154,4 @@ def check_trajectories_with_forward_sensitivities( solver.setRelativeToleranceFSA(1e-13) rdata = runAmiciSimulation(amici_model, solver=solver) _check_close(rdata["x"], result_expected_x, field="x", rtol=1e-10, atol=1e-12) - _check_close(rdata["sx"], result_expected_sx, field="sx", rtol=1e-10, atol=1e-9) + _check_close(rdata["sx"], result_expected_sx, field="sx", rtol=1e-7, atol=1e-9)