-
Notifications
You must be signed in to change notification settings - Fork 43
/
amici_calculator.py
219 lines (193 loc) · 6.92 KB
/
amici_calculator.py
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
import numpy as np
from typing import Dict, List, Sequence, Union
from .constants import (
MODE_FUN, MODE_RES, FVAL, GRAD, HESS, RES, SRES, RDATAS, CHI2
)
from .amici_util import (
add_sim_grad_to_opt_grad, add_sim_hess_to_opt_hess,
sim_sres_to_opt_sres, log_simulation, get_error_output)
try:
import amici
import amici.petab_objective
import amici.parameter_mapping
from amici.parameter_mapping import ParameterMapping
except ImportError:
pass
AmiciModel = Union['amici.Model', 'amici.ModelPtr']
AmiciSolver = Union['amici.Solver', 'amici.SolverPtr']
class AmiciCalculator:
"""
Class to perform the actual call to AMICI and obtain requested objective
function values.
"""
def __init__(self):
self._known_least_squares_safe = False
def initialize(self):
"""Initialize the calculator. Default: Do nothing."""
def __call__(self,
x_dct: Dict,
sensi_order: int,
mode: str,
amici_model: AmiciModel,
amici_solver: AmiciSolver,
edatas: List['amici.ExpData'],
n_threads: int,
x_ids: Sequence[str],
parameter_mapping: 'ParameterMapping',
fim_for_hess: bool):
"""Perform the actual AMICI call.
Called within the :func:`AmiciObjective.__call__` method.
Parameters
----------
x_dct:
Parameters for which to compute function value and derivatives.
sensi_order:
Maximum sensitivity order.
mode:
Call mode (function value or residual based).
amici_model:
The AMICI model.
amici_solver:
The AMICI solver.
edatas:
The experimental data.
n_threads:
Number of threads for AMICI call.
x_ids:
Ids of optimization parameters.
parameter_mapping:
Mapping of optimization to simulation parameters.
fim_for_hess:
Whether to use the FIM (if available) instead of the Hessian (if
requested).
"""
# set order in solver
if sensi_order == 2 and fim_for_hess:
# we use the FIM
amici_solver.setSensitivityOrder(sensi_order-1)
else:
amici_solver.setSensitivityOrder(sensi_order)
# fill in parameters
# TODO (#226) use plist to compute only required derivatives
amici.parameter_mapping.fill_in_parameters(
edatas=edatas,
problem_parameters=x_dct,
scaled_parameters=True,
parameter_mapping=parameter_mapping,
amici_model=amici_model
)
# run amici simulation
rdatas = amici.runAmiciSimulations(
amici_model,
amici_solver,
edatas,
num_threads=min(n_threads, len(edatas)),
)
if not self._known_least_squares_safe and mode == MODE_RES and \
sensi_order > 0:
if any(
((r['ssigmay'] is not None and np.any(r['ssigmay']))
or
(r['ssigmaz'] is not None and np.any(r['ssigmaz'])))
for r in rdatas
):
raise RuntimeError('Cannot use least squares solver with'
'parameter dependent sigma!')
self._known_least_squares_safe = True # don't check this again
return calculate_function_values(
rdatas=rdatas, sensi_order=sensi_order, mode=mode,
amici_model=amici_model, amici_solver=amici_solver, edatas=edatas,
x_ids=x_ids, parameter_mapping=parameter_mapping,
fim_for_hess=fim_for_hess)
def calculate_function_values(rdatas,
sensi_order: int,
mode: str,
amici_model: AmiciModel,
amici_solver: AmiciSolver,
edatas: List['amici.ExpData'],
x_ids: Sequence[str],
parameter_mapping: 'ParameterMapping',
fim_for_hess: bool):
# full optimization problem dimension (including fixed parameters)
dim = len(x_ids)
# check if the simulation failed
if any(rdata['status'] < 0.0 for rdata in rdatas):
return get_error_output(amici_model, edatas, rdatas, dim)
# prepare outputs
nllh = 0.0
snllh = None
s2nllh = None
if mode == MODE_FUN and sensi_order > 0:
snllh = np.zeros(dim)
s2nllh = np.zeros([dim, dim])
chi2 = None
res = None
sres = None
if mode == MODE_RES:
chi2 = 0.0
res = np.zeros([0])
if sensi_order > 0:
sres = np.zeros([0, dim])
par_sim_ids = list(amici_model.getParameterIds())
sensi_method = amici_solver.getSensitivityMethod()
# iterate over return data
for data_ix, rdata in enumerate(rdatas):
log_simulation(data_ix, rdata)
condition_map_sim_var = \
parameter_mapping[data_ix].map_sim_var
# add objective value
nllh -= rdata['llh']
if mode == MODE_FUN:
if sensi_order > 0:
# add gradient
add_sim_grad_to_opt_grad(
x_ids,
par_sim_ids,
condition_map_sim_var,
rdata['sllh'],
snllh,
coefficient=-1.0
)
# Hessian
if sensi_order > 1:
if sensi_method == amici.SensitivityMethod_forward \
and fim_for_hess:
# add FIM for Hessian
add_sim_hess_to_opt_hess(
x_ids,
par_sim_ids,
condition_map_sim_var,
rdata['FIM'],
s2nllh,
coefficient=+1.0
)
else:
raise ValueError("AMICI cannot compute Hessians yet.")
elif mode == MODE_RES:
chi2 += rdata['chi2']
res = np.hstack([res, rdata['res']]) \
if res.size else rdata['res']
if sensi_order > 0:
opt_sres = sim_sres_to_opt_sres(
x_ids,
par_sim_ids,
condition_map_sim_var,
rdata['sres'],
coefficient=1.0
)
sres = np.vstack([sres, opt_sres]) \
if sres.size else opt_sres
ret = {
FVAL: nllh,
CHI2: chi2,
GRAD: snllh,
HESS: s2nllh,
RES: res,
SRES: sres,
RDATAS: rdatas
}
return {
key: val
for key, val in ret.items()
if val is not None
}