From 6bb1d5294b777c310acc9ac1d664ef6188c38ffb Mon Sep 17 00:00:00 2001 From: Jan Grewe Date: Mon, 14 Sep 2020 17:37:37 +0200 Subject: [PATCH] add baseline and no-other responses --- punit_responses.py | 173 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 133 insertions(+), 40 deletions(-) diff --git a/punit_responses.py b/punit_responses.py index 77202a7..4220ac6 100644 --- a/punit_responses.py +++ b/punit_responses.py @@ -3,7 +3,7 @@ import nixio as nix import os from numpy.core.fromnumeric import repeat from traitlets.traitlets import Instance -from chirp_ams import get_signals, other_freq, other_signal, self_freq +from chirp_ams import get_signals from model import simulate, load_models from IPython import embed import matplotlib.pyplot as plt @@ -39,83 +39,176 @@ def save(filename, name, stimulus_settings, model_settings, self_signal, other_s b.metadata = mdata # save stimulus + stim_da = b.create_data_array("complete_stimulus", "nix.timeseries.sampled", dtype=nix.DataType.Float, data=complete_stimulus) stim_da.label = "voltage" - stim_da.label = "mv/cm" + stim_da.label = "mV/cm" dim = stim_da.append_sampled_dimension(model_settings["deltat"]) dim.label = "time" dim.unit = "s" + self_freq_da = None + if self_freq is not None: + self_freq_da = b.create_data_array("self frequency", "nix.timeseries.sampled", dtype=nix.DataType.Float, + data=self_freq) + self_freq_da.label = "frequency" + self_freq_da.label = "Hz" + dim = self_freq_da.append_sampled_dimension(model_settings["deltat"]) + dim.label = "time" + dim.unit = "s" + + other_freq_da = None + if other_freq is not None: + other_freq_da = b.create_data_array("other frequency", "nix.timeseries.sampled", dtype=nix.DataType.Float, + data=self_freq) + self_freq_da.label = "frequency" + self_freq_da.label = "Hz" + dim = other_freq_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]) - + dim = da.append_alias_range_dimension() + # bind it all together + tag = b.create_tag("chirp stimulation", "nix.stimulus", [0.0]) + tag.extent = [stim_da.shape[0] * stim_da.dimensions[0].sampling_interval] + for da in b.data_arrays: + if "response" in da.name.lower(): + tag.references.append(da) + tag.create_feature(stim_da, nix.LinkType.Untagged) + if self_freq_da is not None: + tag.create_feature(self_freq_da, nix.LinkType.Untagged) + if other_freq_da is not None: + tag.create_feature(other_freq_da, nix.LinkType.Untagged) + nf.close() + + +def save_baseline_response(filename, block_name, spikes, model_settings, overwrite=False): + 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 block_name in nf.blocks: + print("Data with this name is already stored! ", block_name) + nf.close() + return + mdata = nf.create_section(block_name, "nix.simulation") + append_settings(mdata, "model parameter", "nix.model.settings", model_settings) + b = nf.create_block(block_name, "nix.simulation") + b.metadata = mdata + da = b.create_data_array("baseline_response", "nix.timeseries.events.spike_times", + dtype=nix.DataType.Float, data=spikes) + da.label = "time" + da.unit = "s" + dim = da.append_alias_range_dimension() nf.close() +def get_baseline_response(model_params, duration=10): + eodf = model_params["EODf"] + dt = model_params["deltat"] + cell_params = model_params.copy() + del cell_params["cell"] + del cell_params["EODf"] + time, pre_stim = get_pre_stimulus(eodf) + baseline_time = np.arange(0.0, duration, dt) + eod = np.sin(baseline_time * eodf * 2 * np.pi) + stim = np.hstack((pre_stim, eod)) + v_0 = np.random.rand(1)[0] + cell_params["v_zero"] = v_0 + spikes = simulate(stim, **cell_params) + spikes = spikes[spikes > time[-1]] - time[-1] + return spikes + + +def get_pre_stimulus(eodf, duration=2, df=20, contrast=20, dt=1./20000): + time = np.arange(0.0, duration, dt) + eod = np.sin(time * eodf * 2 * np.pi) + eod_2 = np.sin(time * (eodf + df) * 2 * np.pi) * contrast/100 + eod_2 *= np.hanning(len(eod_2)) + stim = eod + eod_2 + return time, stim + + def simulate_responses(stimulus_params, model_params, repeats=10, deltaf=20): cell_params = model_params.copy() cell = cell_params["cell"] + filename = "cell_%s.nix" % cell del cell_params["cell"] del cell_params["EODf"] - for c in stimulus_params["contrasts"]: - print("\t\tcreating stimuli ... ") + conditions = ["other", "self"] + + pre_time, pre_stim = get_pre_stimulus(stimulus_params["eodfs"]["self"], dt=model_params["deltat"]) + for contrast in stimulus_params["contrasts"]: params = stimulus_params.copy() + params["contrast"] = contrast 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("\t\tcreating p-unit responses ...") - - spikes = [] - no_other_spikes = [] - for r in range(repeats): - spikes.append(simulate(signal, **cell_params)) - no_other_spikes.append(simulate, self_signal)) - print("\t\tsaving ...") - save("cell_%s.nix" % cell, "contrast_%.3f_condition_%s_deltaf_%i" %(c, stimulus_params["condition"], deltaf), params, - cell_params, self_signal, other_signal, self_freq, other_freq, signal, spikes) + + for condition in conditions: + print("\tcontrast: %s, condition: %s" %(contrast, condition), " "*10, end="\r") + block_name = "contrast_%.3f_condition_%s_deltaf_%i" %(contrast, condition, deltaf) + params["condition"] = condition + time, self_signal, self_freq, other_signal, other_freq = get_signals(**params) + full_signal = (self_signal + other_signal) + + spikes = [] + no_other_spikes = [] + for _ in range(repeats): + v_0 = np.random.rand(1)[0] + cell_params["v_zero"] = v_0 + sp = simulate(np.hstack((pre_stim, full_signal)), **cell_params) + spikes.append(sp[sp > pre_time[-1]] - pre_time[-1]) + if condition == "self": + v_0 = np.random.rand(1)[0] + cell_params["v_zero"] = v_0 + no_other_spikes.append(simulate(np.hstack((pre_stim, self_signal)), **cell_params)) + if condition == "self": + name = "contrast_%.3f_condition_no_other_deltaf_%i" %(contrast, deltaf) + save(filename, name, params, cell_params, self_signal, None, self_freq, None, self_signal, no_other_spikes) + save(filename, block_name, params, cell_params, self_signal, other_signal, self_freq, other_freq, full_signal, spikes) + print("\n") def main(): models = load_models("models.csv") deltafs = [-200, -100, -20, 20, 100, 200] # Hz, difference frequency between self and other + stimulus_params = { "eodfs": {"self": 0.0, "other": 0.0}, + "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": 1, # s, stepsize of the simulation + } for cell_id in range(len(models)): model_params = models[cell_id] + + + baseline_spikes = get_baseline_response(model_params, duration=30) + save_baseline_response( "cell_%s.nix" % model_params["cell"], "baseline response", baseline_spikes, model_params) + print("Cell: %s" % model_params["cell"]) for deltaf in deltafs: + stimulus_params["eodfs"] = {"self": model_params["EODf"], "other": model_params["EODf"] + deltaf} + stimulus_params["dt"] = model_params["deltat"] + print("\t Deltaf: %i" % deltaf) - 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: - print("\t\tcondition: %s" % c) - stimulus_params["condition"] = c - simulate_responses(stimulus_params, model_params, repeats=25, - deltaf=deltaf) - + stimulus_params["chirp_times"] = chirp_times + simulate_responses(stimulus_params, model_params, repeats=25, deltaf=deltaf) + exit() if __name__ == "__main__": main() \ No newline at end of file