Various changes while revising/rerunning experiments (#52)

* progress on all_optical fig

* improve action spectrum extrapolation

* start Newman 15 with alternate opsins

* update Overview doc

* finish newman15 validation w/ alternate opsins

* Update paper plot style

* style_plots_for_paper docstring

* change Jupyter to Python for linguist

* change OGB-1 name

* add warning to scope when sensor not injected

* add showcase to docs

* implement LIOP reset in _base_reset

* add showcase doc to index

* add _base_reset to sim.reset

* fresh run through of all-optical-fig notebook

* add README to notebooks, env info to all_optical_fig.ipynb

* bump to v0.15.0

* remove images [skip ci]

.gitattributes

@@ -1 +1,2 @@
 *.ipynb merge=nbdev-merge
+*.ipynb linguist-language=Python

.gitignore

@@ -6,6 +6,7 @@ tmp/*
 # figures from notebooks not saved by default
 notebooks/img/*
 notebooks/results/*
+docs/tutorials/*.svg
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

cleo/base.py

@@ -245,13 +245,14 @@ def preprocess_ctrl_signals(
         """
         return {}
 
-    def get_intersample_ctrl_signal(self, query_time_ms: float) -> dict:
-        """Get per-stimulator control signal between samples. I.e., for implementing
-        a time-varying waveform based on parameters from the last sample.
-        Such parameters will need to be stored in the :class:`~cleo.IOProcessor`."""
-        return {}
+    def _base_reset(self):
+        """Gets called on :meth:`~cleo.CLSimulator.reset`. The user should not implement this,
+        but rather :meth:`~reset`."""
+        pass
 
     def reset(self, **kwargs) -> None:
+        """Gets called on :meth:`~cleo.CLSimulator.reset`. The user should implement this if
+        their processor has a state that needs to be reset."""
         pass
 
 
@@ -504,6 +505,7 @@ def reset(self, **kwargs):
         for device in self.devices:
             device.reset(**kwargs)
         if self.io_processor is not None:
+            self.io_processor._base_reset()
             self.io_processor.reset(**kwargs)
 
     def to_neo(self) -> neo.core.Block:

cleo/imaging/__init__.py

@@ -1,27 +1,26 @@
 """Contains Scope and sensors for two-photon microscopy"""
 from cleo.imaging.scope import Scope, target_neurons_in_plane
-
 from cleo.imaging.sensors import (
-    Sensor,
     GECI,
-    geci,
     CalBindingActivationModel,
-    NullBindingActivation,
+    CalciumModel,
     DoubExpCalBindingActivation,
+    DynamicCalcium,
     ExcitationModel,
     LightExcitation,
+    NullBindingActivation,
     NullExcitation,
-    CalciumModel,
-    DynamicCalcium,
     PreexistingCalcium,
+    Sensor,
     gcamp3,
+    gcamp6_rs06,
+    gcamp6_rs09,
     gcamp6f,
     gcamp6s,
-    gcamp6rs09,
-    gcamp6rs06,
-    ogb1,
-    jgcamp7f,
-    jgcamp7s,
+    geci,
     jgcamp7b,
     jgcamp7c,
+    jgcamp7f,
+    jgcamp7s,
+    ogb_1,
 )

cleo/imaging/scope.py

@@ -1,18 +1,18 @@
 from __future__ import annotations
-from typing import Any, Callable
+
 import warnings
+from typing import Any, Callable
 
-from attrs import define, field
-from brian2 import NeuronGroup, Unit, mm, np, Quantity, um, meter, ms
 import matplotlib as mpl
+from attrs import define, field
+from brian2 import NeuronGroup, Quantity, Unit, meter, mm, ms, np, um
 from matplotlib.artist import Artist
 from mpl_toolkits.mplot3d import Axes3D
 from nptyping import NDArray
 
 from cleo.base import Recorder
-
-from cleo.imaging.sensors import Sensor
 from cleo.coords import coords_from_ng
+from cleo.imaging.sensors import Sensor
 from cleo.utilities import normalize_coords, rng
 
 
@@ -235,6 +235,11 @@ def get_state(self) -> NDArray[(Any,), float]:
         # separately for each injection, so we'll recover that here
         n_prev_targets_for_ng = {}
         for ng, i_targets in zip(self.neuron_groups, self.i_targets_per_injct):
+            if ng.name not in signal_per_ng:
+                raise RuntimeError(
+                    f"Sensor {self.sensor.name} has no signal for neuron group {ng.name}."
+                    " Did you forget to call inject_sensor_for_targets() after scope injections??"
+                )
             subset_start = n_prev_targets_for_ng.get(ng, 0)
             subset_for_injct = slice(subset_start, subset_start + len(i_targets))
             signal.append(signal_per_ng[ng.name][subset_for_injct])

cleo/imaging/sensors.py

@@ -1,10 +1,11 @@
 from __future__ import annotations
 
-from attrs import define, field, fields_dict, asdict
-from brian2 import Synapses, np, Quantity, NeuronGroup, second, umolar, nmolar
+from attrs import asdict, define, field, fields_dict
+from brian2 import NeuronGroup, Quantity, Synapses, nmolar, np, second, umolar
 
 from cleo.base import SynapseDevice
 from cleo.light import LightDependent
+from cleo.utilities import brian_safe_name
 
 
 @define(eq=False)
@@ -389,45 +390,45 @@ def geci_fn(
 
     geci_fn.__doc__ += "\n\n" + (" " * 8) + extra_doc
 
-    globals()[name] = geci_fn
+    globals()[brian_safe_name(name.lower())] = geci_fn
 
 
 # from NAOMi:
 #                      K_d, n_H, dFF_max, sigma_noise_rel, dFF_1AP_rel, ca_amp, t_on, t_off
-_create_geci_fn("gcamp6f", 290, 2.7, 25.2, 1.24, 0.735, 76.1251, 0.8535, 98.6173)
-_create_geci_fn("gcamp6s", 147, 2.45, 27.2, 1, 1, 54.6943, 0.4526, 68.5461)
+_create_geci_fn("GCaMP6f", 290, 2.7, 25.2, 1.24, 0.735, 76.1251, 0.8535, 98.6173)
+_create_geci_fn("GCaMP6s", 147, 2.45, 27.2, 1, 1, 54.6943, 0.4526, 68.5461)
 # noise, dFF for GCaMP3 from Dana 2019, since not in Zhang 2023
 _create_geci_fn(
-    "gcamp3", 287, 2.52, 12, (3.9 / 2.1) / (13.3 / 4.4), 3.9 / 13.3, 0.05, 1, 1
+    "GCaMP3", 287, 2.52, 12, (3.9 / 2.1) / (13.3 / 4.4), 3.9 / 13.3, 0.05, 1, 1
 )
 
 # from NAOMi, but
 # don't have double exponential convolution parameters for these:
 _create_geci_fn(
-    "ogb1",
+    "OGB-1",
     250,
     1,
     14,
     1,
     extra_doc="*Don't know sigma_noise. Default is the GCaMP6s value.*",
 )
 _create_geci_fn(
-    "gcamp6rs09",
+    "GCaMP6-RS09",
     520,
     3.2,
     25,
     1,
     extra_doc="*Don't know sigma_noise. Default is the GCaMP6s value.*",
 )
 _create_geci_fn(
-    "gcamp6rs06",
+    "GCaMP6-RS06",
     320,
     3,
     15,
     1,
     extra_doc="*Don't know sigma_noise. Default is the GCaMP6s value.*",
 )
-_create_geci_fn("jgcamp7f", 174, 2.3, 30.2, 0.72, 1.71)
-_create_geci_fn("jgcamp7s", 68, 2.49, 40.4, 0.33, 4.96)
-_create_geci_fn("jgcamp7b", 82, 3.06, 22.1, 0.25, 4.64)
-_create_geci_fn("jgcamp7c", 298, 2.44, 145.6, 0.39, 1.85)
+_create_geci_fn("jGCaMP7f", 174, 2.3, 30.2, 0.72, 1.71)
+_create_geci_fn("jGCaMP7s", 68, 2.49, 40.4, 0.33, 4.96)
+_create_geci_fn("jGCaMP7b", 82, 3.06, 22.1, 0.25, 4.64)
+_create_geci_fn("jGCaMP7c", 298, 2.44, 145.6, 0.39, 1.85)

cleo/ioproc/base.py

@@ -232,6 +232,11 @@ def process(self, state_dict: dict, sample_time_ms: float) -> Tuple[dict, float]
         """
         pass
 
+    def _base_reset(self):
+        self.t_samp_ms = []
+        self.out_buffer = deque()
+        self._needs_off_schedule_sample = False
+
 
 class RecordOnlyProcessor(LatencyIOProcessor):
     """Take samples without performing any control.

cleo/light/__init__.py

@@ -10,6 +10,11 @@
     cubic_interpolator,
     equal_photon_flux_spectrum,
     linear_interpolator,
+    log_linear_interpolator,
+    log_makima_interpolator,
+    log_pchip_interpolator,
+    makima_interpolator,
+    pchip_interpolator,
     plot_spectra,
 )
 from cleo.light.two_photon import (

cleo/light/light_dependence.py

@@ -3,22 +3,58 @@
 import warnings
 from typing import Callable, Tuple
 
+import matplotlib.pyplot as plt
 from attrs import define, field
 from brian2 import NeuronGroup, mm, np
-from scipy.interpolate import CubicSpline
+from scipy.interpolate import (
+    Akima1DInterpolator,
+    CubicSpline,
+    PchipInterpolator,
+    interp1d,
+)
 
 from cleo.coords import assign_xyz
 from cleo.utilities import brian_safe_name, wavelength_to_rgb
 
 
 def linear_interpolator(lambdas_nm, epsilons, lambda_new_nm):
-    return np.interp(lambda_new_nm, lambdas_nm, epsilons)
+    # return np.interp(lambda_new_nm, lambdas_nm, epsilons)
+    return interp1d(lambdas_nm, epsilons, fill_value="extrapolate")(lambda_new_nm)
 
 
 def cubic_interpolator(lambdas_nm, epsilons, lambda_new_nm):
     return CubicSpline(lambdas_nm, epsilons)(lambda_new_nm)
 
 
+def pchip_interpolator(lambdas_nm, epsilons, lambda_new_nm):
+    return PchipInterpolator(lambdas_nm, epsilons)(lambda_new_nm)
+
+
+def makima_interpolator(lambdas_nm, epsilons, lambda_new_nm):
+    intrp = Akima1DInterpolator(lambdas_nm, np.log(epsilons))
+    intrp.extrapolate = True
+    intrp.method = "makima"
+    return intrp(lambda_new_nm)
+
+
+def _log_(fn, lambdas_nm, epsilons, lambda_new_nm):
+    assert isinstance(epsilons, np.ndarray)
+    epsilons[epsilons == 0] = 1e-10
+    return np.exp(fn(lambdas_nm, np.log(epsilons), lambda_new_nm))
+
+
+def log_linear_interpolator(lambdas_nm, epsilons, lambda_new_nm):
+    return _log_(linear_interpolator, lambdas_nm, epsilons, lambda_new_nm)
+
+
+def log_pchip_interpolator(lambdas_nm, epsilons, lambda_new_nm):
+    return _log_(pchip_interpolator, lambdas_nm, epsilons, lambda_new_nm)
+
+
+def log_makima_interpolator(lambdas_nm, epsilons, lambda_new_nm):
+    return _log_(makima_interpolator, lambdas_nm, epsilons, lambda_new_nm)
+
+
 # hacky MRO stuff...multiple inheritance only works because slots=False,
 # and must be placed *before* SynapseDevice to work right
 @define(eq=False, slots=False)
@@ -48,10 +84,13 @@ def _default_spectrum(self):
         )
         return [(-1e10, 1), (1e10, 1)]
 
-    spectrum_interpolator: Callable = field(default=cubic_interpolator, repr=False)
+    spectrum_interpolator: Callable = field(default=log_pchip_interpolator, repr=False)
     """Function of signature (lambdas_nm, epsilons, lambda_new_nm) that interpolates
     the action spectrum data and returns :math:`\\varepsilon \\in [0,1]` for the new
     wavelength."""
+    extrapolate: bool = False
+    """Whether or not to attempt to extrapolate (using :attr:`spectrum_interpolator`)
+    outside of the provided excitation/action spectrum."""
 
     @property
     def light_agg_ngs(self):
@@ -81,19 +120,28 @@ def _get_source_for_synapse(
     def epsilon(self, lambda_new) -> float:
         """Returns the :math:`\\varepsilon` value for a given lambda (in nm)
         representing the relative sensitivity of the opsin to that wavelength."""
-        action_spectrum = np.array(self.spectrum)
-        lambdas = action_spectrum[:, 0]
-        epsilons = action_spectrum[:, 1]
-        if lambda_new < min(lambdas) or lambda_new > max(lambdas):
-            warnings.warn(
-                f"λ = {lambda_new} nm is outside the range of the action spectrum data"
-                f" for {self.name}. Assuming ε = 0."
-            )
-            return 0
+        lambdas, epsilons = np.array(self.spectrum).T
         eps_new = self.spectrum_interpolator(lambdas, epsilons, lambda_new)
+
+        # out of data range
+        if lambda_new < min(lambdas) or lambda_new > max(lambdas):
+            if not self.extrapolate:
+                warnings.warn(
+                    f"λ = {lambda_new} nm is outside the range of the action spectrum data"
+                    f" for {self.name} and extrapolate=False. Assuming ε = 0."
+                )
+                return 0
+            elif self.extrapolate:
+                warnings.warn(
+                    f"λ = {lambda_new} nm is outside the range of the action spectrum data"
+                    f" for {self.name}. Extrapolating: ε = {eps_new:.3f}."
+                )
         if eps_new < 0:
             warnings.warn(f"ε = {eps_new} < 0 for {self.name}. Setting ε = 0.")
             eps_new = 0
+        if eps_new > 1:
+            warnings.warn(f"ε = {eps_new} > 1 for {self.name}. Setting ε = 1.")
+            eps_new = 1
         return eps_new
 
 
@@ -109,19 +157,27 @@ def equal_photon_flux_spectrum(
     return list(zip(lambdas, eps_Irr))
 
 
-def plot_spectra(*ldds: LightDependent) -> tuple[plt.Figure, plt.Axes]:
+def plot_spectra(
+    *ldds: LightDependent, extrapolate=False
+) -> tuple[plt.Figure, plt.Axes]:
     """Plots the action/excitation spectra for multiple light-dependent devices."""
     import matplotlib.pyplot as plt
 
+    if extrapolate:
+        all_lambdas = [np.array(ldd.spectrum)[:, 0] for ldd in ldds]
+        lambda_min = min([lambdas.min() for lambdas in all_lambdas])
+        lambda_max = max([lambdas.max() for lambdas in all_lambdas])
+
     fig, ax = plt.subplots()
     for ldd in ldds:
-        spectrum = np.array(ldd.spectrum)
-        lambdas = spectrum[:, 0]
-        epsilons = spectrum[:, 1]
-        lambdas_new = np.linspace(min(lambdas), max(lambdas), 100)
+        lambdas, epsilons = np.array(ldd.spectrum).T
+        if not extrapolate:
+            lambda_min = np.min(lambdas)
+            lambda_max = np.max(lambdas)
+        lambdas_new = np.linspace(lambda_min, lambda_max, 100)
         epsilons_new = ldd.spectrum_interpolator(lambdas, epsilons, lambdas_new)
         c_points = [wavelength_to_rgb(l) for l in lambdas]
-        c_line = wavelength_to_rgb(lambdas_new[np.argmax(epsilons_new)])
+        c_line = wavelength_to_rgb(lambdas[np.argmax(epsilons)])
         ax.plot(lambdas_new, epsilons_new, c=c_line, label=ldd.name)
         ax.scatter(lambdas, epsilons, marker="o", s=50, color=c_points)
     title = (

cleo/opto/opsin_library.py

@@ -324,5 +324,6 @@ def enphr3_3s():
         chrimson_4s(),
         gtacr2_4s(),
         enphr3_3s(),
+        extrapolate=True,
     )
     plt.show()

cleo/utilities.py

@@ -295,13 +295,22 @@ def style_plots_for_docs(dark=True):
     plt.rc("font", **{"sans-serif": "Open Sans"})
 
 
-def style_plots_for_paper():
+def style_plots_for_paper(fontscale=5 / 6):
+    """
+    fontscale=5/6 goes from default paper font size of 9.6 down to 8
+    """
     # some hacky workaround for params not being updated until after first plot
     f = plt.figure()
     plt.plot()
     plt.close(f)
 
-    plt.style.use("seaborn-v0_8-paper")
+    try:
+        import seaborn as sns
+
+        sns.set_context("paper", font_scale=fontscale)
+    except ImportError:
+        plt.style.use("seaborn-v0_8-paper")
+        warnings.warn("Seaborn not found, using matplotlib style, ignoring fontscale")
     plt.rc("savefig", transparent=True, bbox="tight", dpi=300)
     plt.rc("svg", fonttype="none")
     plt.rc("axes.spines", top=False, right=False)

docs/index.rst

@@ -15,6 +15,7 @@ Documentation contents
 
    overview
    tutorials/index
+   showcase
    reference/index