v0.5.2 (#159)

* bump version

* add back fast Nanduri model cascade

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
 

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).

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)

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,36 +35,40 @@ 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
 
         # 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

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',

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