Skip to content

Commit 6761859

Browse files
committed
Merge branch 'custom_cost_func'
2 parents 18ebf5b + 447eae5 commit 6761859

File tree

2 files changed

+84
-31
lines changed

2 files changed

+84
-31
lines changed

README.md

Lines changed: 30 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,8 @@
22

33
![Python version badge](https://img.shields.io/badge/python-3.6-blue.svg)
44
[![license](https://img.shields.io/github/license/LoLab-VU/swarm_it.svg)](LICENSE)
5-
![version](https://img.shields.io/badge/version-0.1.1-orange.svg)
6-
[![release](https://img.shields.io/github/release-pre/LoLab-VU/swarm_it.svg)](https://github.com/LoLab-VU/swarm_it/releases/tag/v0.1.1)
5+
![version](https://img.shields.io/badge/version-0.2.0-orange.svg)
6+
[![release](https://img.shields.io/github/release-pre/LoLab-VU/swarm_it.svg)](https://github.com/LoLab-VU/swarm_it/releases/tag/v0.2.0)
77

88
**swarm_it** is a python utility designed to abstract away some of the effort in setting up Particle Swarm Optimization (PSO)-based calibrations of biological models in the [PySB](http://pysb.org/) format using [simplePSO](https://github.com/LoLab-VU/ParticleSwarmOptimization).
99

@@ -77,6 +77,7 @@ data = np.load('my_data.npy')
7777
data_sd = np.load('my_data_sd.npy')
7878
observable_data = dict()
7979
time_idxs = list(range(len(timespan)))
80+
# my_observable is the name of an Observable in the model.
8081
observable_data['my_observable'] = (data, data_sd, time_idxs)
8182
# Initialize the SwarmIt instance with the model details.
8283
swarmit = SwarmIt(my_model, observable_data, timespan)
@@ -91,31 +92,44 @@ pso.run(num_particles, num_iterations)
9192
print("best_theta: ",pso.best)
9293
```
9394

94-
SwarmIt constructs the PSO object to sample all of a model's kinetic rate parameters. It assumes that the priors are uniform with size 4 orders of magnitude and centered on the values defined in the model.
95+
#### SwarmParam
96+
By default SwarmIt constructs the PSO object to sample all of a model's kinetic rate parameters; it assumes that the priors are uniform with size 4 orders of magnitude and centered on the values defined in the model. However, the swarm_it module has a built-in helper class, SwarmParam, which can be used in conjunction of with SwarmIt class to set which model parameters are sampled and the corresponding sampling ranges. SwarmParam can be used at the level of PySB model definition to log which parameters to include in a Particle Swarm Optimization run, or it can be used after model definition at as long as it is defined before defining a SwarmIt object. It can be imported from swarm_it:
97+
```python
98+
from swarm_it import SwarmParam
99+
```
100+
It is passed at instantiation of the SwarmIt class via the `swarm_param` input parameter,
101+
```python
102+
swarmit = SwarmIt(model, observable_data, tspan, swarm_param=swarm_param)
103+
```
104+
which uses it to build the set of parameters to sample, their staring PSO position and bounds, as well as the parameter mask for use with built-in cost functions.
95105

96-
In addition, SwarmIt crrently has three pre-defined loglikelihood functions with different estimators, specified with the keyword parameter cost_type:
106+
Note that if you flag a parameter for sampling without setting sampling bounds, SwarmParam will by default assign the parameter bounds centered on the parameter's value (as defined in the model) with a width of 4 orders of magnitude.
107+
108+
For a more in depth example usage of SwarmParam, see the [simplepso_dimerization_model_SwarmIt_SwarmParam example](./examples/pysb_dimerization_model/simplepso_dimerization_model_SwarmIt_SwarmParam.py) and corresponding model [dimerization_model_swarmit](./examples/pysb_dimerization_model/dimerization_model_swarmit.py).
109+
110+
##### Built-in cost functions
111+
SwarmIt currently has three pre-defined cost functions with different estimators, specified with the keyword parameter cost_type:
97112
```python
98113
# Now build the PSO object.
99114
pso = swarmit(cost_type='mse')
100115
```
101116
The options are
102-
* 'norm_logpdf'=> (default) Compute the cost a the negative sum of normal distribution logpdf's over each data point
103-
* 'mse'=>Compute the cost using the mean squared error estimator
104-
* 'sse'=>Compute the cost using the sum of squared errors estimator.
117+
* 'norm_logpdf'=> (default) Compute the cost as the negative sum of normal distribution logpdf's over each timecourse data points and each given model observable.
118+
* 'mse'=>Compute the cost using the mean squared error estimator over the timecourse data and each given model observable.
119+
* 'sse'=>Compute the cost using the sum of squared errors estimator over the timecourse data and each given observable.
105120

106121
Each of these functions computes the cost estimate using the timecourse output of a model simulation for each observable defined in the `observable_data` dictionary.
107-
If you want to use a different or more complicated likelihood function with SwarmIt then you'll need to subclass it and either override one of the existing cost functions or define a new one.
108122

109-
#### SwarmParam
110-
The swarm_it module has a built-in helper class, SwarmParam, which can be used in conjunction of with SwarmIt class. SwarmParam can be used at the level of PySB model definition to log which parameters to include in
111-
a Particle Swarm Optimization run, or it can be used after model definition at as long as it is defined before defining a SwarmIt object. It can be imported from swarm_it:
123+
##### Custom cost function
124+
If you want to use a different or more complicated cost function than the built-ins, the `SwarmIt` class object can also accept a custom cost function defined by the user with the parameters `cost_type='custom'` and `custom_cost`:
112125
```python
113-
from swarm_it import SwarmParam
126+
# Now build the PSO object with a custom cost function.
127+
def user_cost(model, simulation_output):
128+
# Do something and compute the cost for this model and its evaluation.
129+
return cost
130+
pso = swarmit(cost_type='custom', custom_cost=user_cost)
114131
```
115-
It is passed at instantiation to the SwarmIt class, which uses it
116-
to build the set of parameters to use, their staring PSO position and bounds, as well as the parameter mask for the cost function.
117-
118-
Note that if you flag a parameter for sampling without setting sampling bounds, SwarmParam will by default assign the parameter bounds centered on the parameter's value (as defined in the model) with a width of 4 orders of magnitude.
132+
Note that the custom cost function defined by the user, `user_cost`, should accept two parameters: the model itself and the simulation output (or trajectory). The user can then define how to compute the cost using these input objects or using other data/objects defined outside `user_cost`.
119133

120134
### Examples
121135
Additional example scripts that show how to setup and launch PSO runs using **swarm_it** can be found under [examples](./examples).

swarm_it.py

Lines changed: 54 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@
1111
"""
1212

1313
import importlib
14+
import warnings
1415
import os.path
1516
try:
1617
import pysb
@@ -153,9 +154,15 @@ def upper(self):
153154
def add_all_kinetic_params(self, pysb_model):
154155
for rule in pysb_model.rules:
155156
if rule.rate_forward:
156-
self.__call__(rule.rate_forward)
157+
try:
158+
self.__call__(rule.rate_forward)
159+
except:
160+
pass
157161
if rule.rate_reverse:
158-
self.__call__(rule.rate_reverse)
162+
try:
163+
self.__call__(rule.rate_reverse)
164+
except:
165+
pass
159166
return
160167

161168
def add_all_nonkinetic_params(self, pysb_model):
@@ -228,14 +235,17 @@ def __init__(self, model, observable_data, timespan,
228235
self._like_data = dict()
229236
self._data = dict()
230237
self._data_mask = dict()
231-
for observable_key in observable_data.keys():
232-
self._like_data[observable_key] = norm(loc=observable_data[observable_key][0],
233-
scale=observable_data[observable_key][1])
234-
self._data[observable_key] = observable_data[observable_key][0]
235-
self._data_mask[observable_key] = observable_data[observable_key][2]
236-
# print(observable_data[observable_key][2])
237-
if observable_data[observable_key][2] is None:
238-
self._data_mask[observable_key] = range(len(self.timespan))
238+
if observable_data is not None:
239+
for observable_key in observable_data.keys():
240+
self._like_data[observable_key] = norm(loc=observable_data[observable_key][0],
241+
scale=observable_data[observable_key][1])
242+
self._data[observable_key] = observable_data[observable_key][0]
243+
self._data_mask[observable_key] = observable_data[observable_key][2]
244+
# print(observable_data[observable_key][2])
245+
if observable_data[observable_key][2] is None:
246+
self._data_mask[observable_key] = range(len(self.timespan))
247+
else:
248+
warnings.warn('The observable_data input was set to None, which is only compatible with the use custom cost functions.')
239249
self._model_solver = solver(self.model, tspan=self.timespan, **solver_kwargs)
240250
if swarm_param is not None:
241251
parm_mask = swarm_param.mask(model.parameters)
@@ -246,17 +256,18 @@ def __init__(self, model, observable_data, timespan,
246256
else:
247257
swarm_param = SwarmParam()
248258
for rule in model.rules:
249-
if rule.rate_forward:
250-
swarm_param(rule.rate_forward)
251-
if rule.rate_reverse:
252-
swarm_param(rule.rate_reverse)
259+
if rule.rate_forward and (rule.rate_forward in model.parameters):
260+
swarm_param(rule.rate_forward)
261+
if rule.rate_reverse and (rule.rate_reverse in model.parameters):
262+
swarm_param(rule.rate_reverse)
253263
parm_mask = swarm_param.mask(model.parameters)
254264
# self._sampled_parameters = [SampledParameter(parm.name, *swarm_param[parm.name]) for i,parm in enumerate(model.parameters) if parm_mask[i]]
255265
self._rate_mask = parm_mask
256266
self._starting_position = swarm_param.centers()
257267
self._lower = swarm_param.lower()
258268
self._upper = swarm_param.upper()
259269
self._param_values = np.array([param.value for param in model.parameters])
270+
self._custom_cost = None
260271
return
261272

262273
def norm_logpdf_cost(self, position):
@@ -327,8 +338,33 @@ def sse_cost(self, position):
327338
return np.inf
328339
return logl,
329340

341+
def custom_cost(self, position):
342+
"""Compute the cost using the negative sum of squared errors estimator.
343+
344+
Args:
345+
position (numpy.array): The parameter vector the compute cost
346+
of.
347+
348+
Returns:
349+
float: The natural logarithm of the likelihood estimate.
350+
351+
"""
352+
Y = np.copy(position)
353+
params = self._param_values.copy()
354+
params[self._rate_mask] = 10.**Y
355+
sim = self._model_solver.run(param_values=[params]).all
356+
logl = self._custom_cost(self.model, sim)
357+
if np.isnan(logl):
358+
return np.inf
359+
return logl,
360+
361+
def set_custom_cost(self, cost_func):
362+
363+
self._custom_cost = cost_func
364+
return
365+
330366
def __call__(self, pso_kwargs=None,
331-
cost_type='norm_logpdf'):
367+
cost_type='norm_logpdf', custom_cost=None):
332368
"""Call the SwarmIt instance to construct to instance of the NestedSampling object.
333369
334370
Args:
@@ -355,6 +391,9 @@ def __call__(self, pso_kwargs=None,
355391
cost = self.mse_cost
356392
elif cost_type == 'sse':
357393
cost = self.sse_cost
394+
elif (cost_type == 'custom') and (custom_cost is not None):
395+
self.set_custom_cost(custom_cost)
396+
cost = self.custom_cost
358397
else:
359398
cost = self.norm_logpdf_cost
360399

0 commit comments

Comments
 (0)