From e9c1879e22a3349db92b82f5446dae3af3795eaa Mon Sep 17 00:00:00 2001 From: Michael Beyeler Date: Tue, 25 Feb 2020 12:29:05 -0800 Subject: [PATCH] v0.5.2 (#159) * bump version * add back fast Nanduri model cascade --- README.md | 2 +- doc/index.rst | 4 ++-- pulse2percept/api.py | 5 +++- pulse2percept/fast_retina.pyx | 44 +++++++++++++++++++---------------- pulse2percept/retina.py | 11 ++------- pulse2percept/version.py | 2 +- 6 files changed, 34 insertions(+), 34 deletions(-) diff --git a/README.md b/README.md index 64af7a33..c8477d3f 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,7 @@ [![Documentation Status](https://readthedocs.org/projects/pulse2percept/badge/?version=latest)](https://pulse2percept.readthedocs.io/en/latest/?badge=latest) [![license](https://img.shields.io/badge/License-BSD%203--Clause-blue.svg)](https://github.com/uwescience/pulse2percept/blob/master/LICENSE) -# pulse2percept: A Python-based simulation framework for bionic vision +# pulse2percept 0.5.2: A Python-based simulation framework for bionic vision ## Summary diff --git a/doc/index.rst b/doc/index.rst index 035c02e7..c41f2270 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -15,8 +15,8 @@ developers/contributing developers/releases -pulse2percept: Models for sight restoration -=========================================== +pulse2percept |version| documentation +===================================== By 2020 roughly 200 million people worldwide will suffer from retinal diseases such as retinitis pigmentosa (RP) and age-related macular degeneration (AMD). diff --git a/pulse2percept/api.py b/pulse2percept/api.py index 53511d87..af75c7f0 100644 --- a/pulse2percept/api.py +++ b/pulse2percept/api.py @@ -271,6 +271,9 @@ def set_ganglion_cell_layer(self, model, **kwargs): shift : float, optional, default: 16.0 Shift of the logistic function in the stationary nonlinearity stage. + use_cython : bool, optional, default: True + Flag whether to use the fast Cython implementation (True) + or the slower NumPy implementation (False) - 'Horsager2009': A model of temporal sensitivity as described in [2]_. @@ -486,7 +489,7 @@ def pulse2percept(self, stim, t_percept=None, tol=0.05, engine=self.engine, scheduler=self.scheduler, func_args=[pt_data, layers, self.use_jit]) bm = np.zeros(self.ofl.gridx.shape - + (sr_list[0].data.shape[-1], )) + + (sr_list[0].data.shape[-1], )) idxer = tuple(np.array(idx_list)[:, i] for i in range(2)) bm[idxer] = [sr.data for sr in sr_list] percept = utils.TimeSeries(sr_list[0].tsample, bm) diff --git a/pulse2percept/fast_retina.pyx b/pulse2percept/fast_retina.pyx index 7cdac503..f9dc4801 100644 --- a/pulse2percept/fast_retina.pyx +++ b/pulse2percept/fast_retina.pyx @@ -2,11 +2,11 @@ import numpy as np cimport numpy as np import cython -cdef extern from "math.h": - cpdef float expf(float x) +from libc.math cimport(exp as c_exp) + cdef inline float expit(float x): - return 1.0 / (1.0 + expf(-x)) + return 1.0 / (1.0 + c_exp(-x)) cdef inline float float_max(float a, float b): return a if a >= b else b @@ -35,25 +35,20 @@ def nanduri2012_calc_layer_current(double[:] in_arr, double[:, ::1] pt_arr): def nanduri2012_model_cascade(double[:] stimdata, float dt, float tau1, float tau2, float tau3, float asymptote, float shift, - float slope, float eps, float max_R3): + float slope, float eps): """Cython implementation of the Nanduri 2012 model cascade""" - cdef float tmp_ca = 0 - cdef float tmp_R1 = 0 - cdef float tmp_R2 = 0 - cdef float tmp_R3 = 0 - cdef float tmp_R4a = 0 - cdef float tmp_R4b = 0 - cdef float tmp_R4c = 0 - # Stationary nonlinearity: used to depend on future values of the - # intermediary response, now has to be passed through `max_R3` - # (because we can't look into the future): - cdef float scale = asymptote * expit((max_R3 - shift) / slope) + # Because the stationary nonlinearity depends on `max_R3`, which is the + # largest value of R3 over all time points, we have to process the stimulus + # in two steps: - # Output array + # Step 1: Calculate `tmp_R3` for all time points and extract `max_R3`: + cdef float tmp_ca = 0.0 + cdef float tmp_R1 = 0.0 + cdef float tmp_R2 = 0.0 + cdef float max_R3 = 1e-12 # avoid division by zero cdef np.intp_t arr_size = len(stimdata) - cdef double[:] out_R4 = np.empty(arr_size, dtype=float) - + cdef double[:] tmp_R3 = np.empty(arr_size, dtype=float) for i in range(arr_size): # Fast ganglion cell response: tmp_R1 += dt * (-stimdata[i] - tmp_R1) / tau1 @@ -61,10 +56,19 @@ def nanduri2012_model_cascade(double[:] stimdata, # Leaky integrated charge accumulation: tmp_ca += dt * float_max(stimdata[i], 0) tmp_R2 += dt * (tmp_ca - tmp_R2) / tau2 - tmp_R3 = float_max(tmp_R1 - eps * tmp_R2, 0) / max_R3 * scale + tmp_R3[i] = float_max(tmp_R1 - eps * tmp_R2, 0) + if tmp_R3[i] > max_R3: + max_R3 = tmp_R3[i] + # Step 2: Calculate `out_R4` from `tmp_R3` + cdef float tmp_R4a = 0.0 + cdef float tmp_R4b = 0.0 + cdef float tmp_R4c = 0.0 + cdef float scale = asymptote * expit((max_R3 - shift) / slope) / max_R3 + cdef double[:] out_R4 = np.empty(arr_size, dtype=float) + for i in range(arr_size): # Slow response: 3-stage leaky integrator - tmp_R4a += dt * (tmp_R3 - tmp_R4a) / tau3 + tmp_R4a += dt * (tmp_R3[i] * scale - tmp_R4a) / tau3 tmp_R4b += dt * (tmp_R4a - tmp_R4b) / tau3 tmp_R4c += dt * (tmp_R4b - tmp_R4c) / tau3 out_R4[i] = tmp_R4c diff --git a/pulse2percept/retina.py b/pulse2percept/retina.py index bcf1ce76..1b513aa6 100644 --- a/pulse2percept/retina.py +++ b/pulse2percept/retina.py @@ -575,13 +575,7 @@ def __init__(self, **kwargs): self.asymptote = 14.0 self.slope = 3.0 self.shift = 16.0 - self.use_cython = False - - # Nanduri (2012) has a term in the stationary nonlinearity step that - # depends on future values of R3: max_t(R3). Because the finite - # difference model cannot look into the future, we need to set a - # scaling factor here: - self.maxR3 = 100.0 + self.use_cython = True # perform one-time setup calculations # gamma1 is used for the fast response @@ -651,8 +645,7 @@ def model_cascade(self, in_arr, pt_list, layers, use_jit): self.tau1, self.tau2, self.tau3, self.asymptote, self.shift, - self.slope, self.eps, - self.maxR3) + self.slope, self.eps) else: # Fast response b2 = self.tsample * utils.conv(-b1, self.gamma1, mode='full', diff --git a/pulse2percept/version.py b/pulse2percept/version.py index 50de1329..4f5181d4 100644 --- a/pulse2percept/version.py +++ b/pulse2percept/version.py @@ -4,7 +4,7 @@ # Format expected by setup.py and doc/source/conf.py: string of form "X.Y.Z" _version_major = 0 _version_minor = 5 -_version_micro = 1 # use '' for first of series, number for 1 and above +_version_micro = 2 # use '' for first of series, number for 1 and above _version_extra = '' # _version_extra = '' # Uncomment this for full releases