add model and p-unit respones

model.py

@@ -0,0 +1,83 @@
+from IPython.terminal.embed import embed
+from numba.core.types import float64
+import numpy as np
+
+try:
+    from numba import jit
+except ImportError:
+    def jit(nopython):
+        def decorator_jit(func):
+            return func
+        return decorator_jit
+
+
+def load_models(file):
+    """ Load model parameter from csv file.
+
+    Parameters
+    ----------
+    file: string
+        Name of file with model parameters.
+
+    Returns
+    -------
+    parameters: list of dict
+        For each cell a dictionary with model parameters.
+    """
+    parameters = []
+    with open(file, 'r') as file:
+        header_line = file.readline()
+        header_parts = header_line.strip().split(",")
+        keys = header_parts
+        for line in file:
+            line_parts = line.strip().split(",")
+            parameter = {}
+            for i in range(len(keys)):
+                parameter[keys[i]] = float(line_parts[i]) if i > 0 else line_parts[i]
+            parameters.append(parameter)
+    return parameters
+
+
+@jit(nopython=True)
+def simulate(stimulus, deltat=0.00005, v_zero=0.0, a_zero=2.0, threshold=1.0, v_base=0.0,
+             delta_a=0.08, tau_a=0.1, v_offset=-10.0, mem_tau=0.015, noise_strength=0.05,
+             input_scaling=60.0, dend_tau=0.001, ref_period=0.001, **kwargs):
+    """ Simulate a P-unit.
+
+    Returns
+    -------
+    spike_times: 1-D array
+        Simulated spike times in seconds.
+    """    
+    # initial conditions:
+    v_dend = stimulus[0]
+    v_mem = v_zero
+    adapt = a_zero
+
+    # prepare noise:    
+    noise = np.random.randn(len(stimulus))
+    noise *= noise_strength / np.sqrt(deltat)
+
+    # rectify stimulus array:
+    stimulus = stimulus.copy()
+    stimulus[stimulus < 0.0] = 0.0
+
+    # integrate:
+    spike_times = []
+    for i in range(len(stimulus)):
+        v_dend += (-v_dend + stimulus[i]) / dend_tau * deltat
+        v_mem += (v_base - v_mem + v_offset + (
+                    v_dend * input_scaling) - adapt + noise[i]) / mem_tau * deltat
+        adapt += -adapt / tau_a * deltat
+
+        # refractory period:
+        if len(spike_times) > 0 and (deltat * i) - spike_times[-1] < ref_period + deltat/2:
+            v_mem = v_base
+
+        # threshold crossing:
+        if v_mem > threshold:
+            v_mem = v_base
+            spike_times.append(i * deltat)
+            adapt += delta_a / tau_a
+
+    return np.array(spike_times)

punit_responses.py

@@ -0,0 +1,119 @@
+import numpy as np
+import nixio as nix
+import os
+from numpy.core.fromnumeric import repeat
+from traitlets.traitlets import Instance
+from chirp_ams import get_signals
+from model import simulate, load_models
+from IPython import embed
+import matplotlib.pyplot as plt
+
+def append_settings(section, sec_name, sec_type, settings):
+    section = section.create_section(sec_name, sec_type)
+    for k in settings.keys():
+        if isinstance(settings[k], dict):
+            append_settings(section, k, "settings", settings[k])
+        else:
+            if isinstance(settings[k], np.ndarray):
+                if len(settings[k].shape) == 1:
+                    section[k] = list(settings[k])
+            else:
+                section[k] = settings[k]
+
+    
+def save(filename, name, stimulus_settings, model_settings, stimulus, responses, overwrite=False):
+    print("saving! ", filename, name)
+    if os.path.exists(filename) and not overwrite:
+        nf = nix.File.open(filename, nix.FileMode.ReadWrite)
+    else:
+        nf = nix.File.open(filename, mode=nix.FileMode.Overwrite,
+                               compression=nix.Compression.DeflateNormal)
+    
+    if name in nf.blocks:
+        print("Data with this name is already stored! ", name)
+        nf.close()
+        return
+    mdata = nf.create_section(name, "nix.simulation")
+    append_settings(mdata, "model parameter", "nix.model.settings", model_settings)
+    append_settings(mdata, "stimulus parameter", "nix.stimulus.settings", stimulus_settings)
+    b = nf.create_block(name, "nix.simulation")
+    b.metadata = mdata
+    
+    # save stimulus
+    stim_da = b.create_data_array("stimulus", "nix.timeseries.sampled", dtype=nix.DataType.Float,
+                                  data=stimulus)
+    stim_da.label = "voltage"
+    stim_da.label = "mv/cm"
+    dim = stim_da.append_sampled_dimension(model_settings["deltat"])
+    dim.label = "time"
+    dim.unit = "s"
+    
+    # save responses
+    for i in range(len(responses)):
+        da = b.create_data_array("response_%i" %i, "nix.timeseries.events.spike_times", 
+                                 dtype=nix.DataType.Float, data=responses[i])
+        da.label = "time"
+        da.unit = "s"
+        dim = da.append_range_dimension()
+        dim.link_data_array(da, [-1])
+    nf.close()
+    pass
+
+
+def plot_responses():
+    pass
+
+
+def simulate_responses(stimulus_params, model_params, repeats=10):
+    cell_params = model_params.copy()
+    cell = cell_params["cell"]
+    del cell_params["cell"]
+    del cell_params["EODf"]
+    for c in stimulus_params["contrasts"]:
+        print("creating stimuli\n\tcontrast: ", str(c), "\t condition: ", 
+              stimulus_params["condition"])
+        params = stimulus_params.copy()
+        del params["contrasts"]
+        del params["chirp_frequency"]
+        params["contrast"] = c
+        time, self_signal, self_freq, other_signal, other_freq = get_signals(**params)
+        signal = (self_signal + other_signal)
+        signal /= np.max(signal)
+        print("create p-unit responses for cell: ", cell)
+        
+        spikes = []
+        for r in range(repeats):
+            spikes.append(simulate(signal, **cell_params))
+        save("test.nix", "contrast_%.3f_condition_%s" %(c, stimulus_params["condition"]), params,
+             cell_params, signal, spikes)
+
+
+def main():
+    cell_id = 20
+    models = load_models("models.csv")
+    deltaf = 20.  # Hz, difference frequency between self and other
+    model_params = models[cell_id]
+    stimulus_params = {"eodfs": {"self": model_params["EODf"],
+                                           "other": model_params["EODf"] + deltaf},
+                       "contrasts": [20, 10, 5, 2.5, 1.25, 0.625, 0.3125],
+                       "chirp_size": 100,  # Hz, frequency excursion
+                       "chirp_duration": 0.015,  # s, chirp duration
+                       "chirp_amplitude_dip": 0.05,  # %, amplitude drop during chirp
+                       "chirp_frequency": 5,  # Hz, how often does the fish chirp
+                       "duration": 5.,  # s, total duration of simulation
+                       "dt": model_params["deltat"],  # s, stepsize of the simulation
+    }
+    chirp_times = np.arange(stimulus_params["chirp_duration"],
+                            stimulus_params["duration"] - stimulus_params["chirp_duration"], 
+                            1./stimulus_params["chirp_frequency"])
+    stimulus_params["chirp_times"] = chirp_times
+    conditions = ["other", "self"]
+   
+    for c in conditions:
+        stimulus_params["condition"] = c
+        simulate_responses(stimulus_params, model_params, repeats=25)
+    pass
+
+
+if __name__ == "__main__":
+    main()