Skip to content

fix rounding errors in Hill functions #183

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
murrayrm wants to merge 1 commit into
base: master
Choose a base branch
from
Open

fix rounding errors in Hill functions #183

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions bioscrape/types.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,10 @@ from bioscrape.sbmlutil import add_species, add_parameter, add_reaction, add_rul

from libc.math cimport log, sqrt, cos, round, exp, fabs

# Define static epsilon for handling rounding errors
from libcpp.limits cimport numeric_limits
cdef float float_epsilon = numeric_limits[float].epsilon()

################################################## ####################################################
###################################### PROPENSITY TYPES ##############################
################################################# ################################################
Expand Down Expand Up @@ -222,13 +226,17 @@ cdef class PositiveHillPropensity(Propensity):
cdef double K = params[self.K_index]
cdef double n = params[self.n_index]
cdef double rate = params[self.rate_index]
if X < 0 and abs(X) < float_epsilon:
X = 0. # Fix rounding errors
return rate * (X / K) ** n / (1 + (X/K)**n)

cdef double get_volume_propensity(self, double *state, double *params, double volume, double time):
cdef double X = state[self.s1_index] / volume
cdef double K = params[self.K_index]
cdef double n = params[self.n_index]
cdef double rate = params[self.rate_index]
if X < 0 and abs(X) < float_epsilon:
X = 0. # Fix rounding errors
return rate * (X / K) ** n / (1 + (X/K)**n)

def initialize(self, dict param_dictionary, dict species_indices, dict parameter_indices):
Expand Down Expand Up @@ -266,6 +274,8 @@ cdef class PositiveProportionalHillPropensity(Propensity):
cdef double n = params[self.n_index]
cdef double rate = params[self.rate_index]
cdef double d = state[self.d_index]
if X < 0 and abs(X) < float_epsilon:
X = 0. # Fix rounding errors
return rate * d * (X / K) ** n / (1 + (X/K)**n)

cdef double get_volume_propensity(self, double *state, double *params, double volume, double time):
Expand All @@ -274,6 +284,8 @@ cdef class PositiveProportionalHillPropensity(Propensity):
cdef double n = params[self.n_index]
cdef double d = state[self.d_index]
cdef double rate = params[self.rate_index]
if X < 0 and abs(X) < float_epsilon:
X = 0. # Fix rounding errors
return d * rate * (X / K) ** n / (1 + (X/K)**n)


Expand Down Expand Up @@ -357,6 +369,8 @@ cdef class NegativeProportionalHillPropensity(Propensity):
cdef double n = params[self.n_index]
cdef double rate = params[self.rate_index]
cdef double d = state[self.d_index]
if X < 0 and abs(X) < float_epsilon:
X = 0. # Fix rounding errors
return rate * d * 1/ (1 + (X/K)**n)

cdef double get_volume_propensity(self, double *state, double *params, double volume, double time):
Expand All @@ -365,6 +379,8 @@ cdef class NegativeProportionalHillPropensity(Propensity):
cdef double n = params[self.n_index]
cdef double d = state[self.d_index]
cdef double rate = params[self.rate_index]
if X < 0 and abs(X) < float_epsilon:
X = 0. # Fix rounding errors
return d * rate * 1 / (1 + (X/K)**n)


Expand Down