From c64add286fba71682ca60809422534e036a53616 Mon Sep 17 00:00:00 2001 From: "a.ott" Date: Tue, 14 Jul 2020 17:25:42 +0200 Subject: [PATCH] add scripts to create plots and read info from data --- Figures_Baseline.py | 115 ++++++++++++++++++++++++++++++++++++++++++ Figures_Model.py | 56 +++++++++++++++++++++ Figures_Stimuli.py | 67 +++++++++++++++++++++++++ read_data_infos.py | 118 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 356 insertions(+) create mode 100644 Figures_Baseline.py create mode 100644 Figures_Model.py create mode 100644 Figures_Stimuli.py create mode 100644 read_data_infos.py diff --git a/Figures_Baseline.py b/Figures_Baseline.py new file mode 100644 index 0000000..df4ab0a --- /dev/null +++ b/Figures_Baseline.py @@ -0,0 +1,115 @@ + + +import matplotlib.pyplot as plt +import numpy as np +import os + +import helperFunctions as hF + +from CellData import CellData +from Baseline import BaselineCellData +from FiCurve import FICurveCellData + +MODEL_COLOR = "blue" +DATA_COLOR = "orange" + +DATA_SAVE_PATH = "data/figure_data/" + +def main(): + # data_isi_histogram() + # data_mean_freq_step_stimulus_examples() + # data_mean_freq_step_stimulus_with_detections() + data_fi_curve() + pass + + +def data_fi_curve(): + cell = "data/invivo/2013-04-17-ac-invivo-1/" + cell_data = CellData(cell) + fi = FICurveCellData(cell_data, cell_data.get_fi_contrasts()) + fi.plot_fi_curve() + +def data_mean_freq_step_stimulus_with_detections(): + cell = "data/invivo/2013-04-17-ac-invivo-1/" + cell_data = CellData(cell) + fi = FICurveCellData(cell_data, cell_data.get_fi_contrasts()) + + time_traces, freq_traces = fi.get_time_and_freq_traces() + mean_times, mean_freqs = fi.get_mean_time_and_freq_traces() + idx = -1 + time = np.array(mean_times[idx]) + freq = np.array(mean_freqs[idx]) + f_inf = fi.f_inf_frequencies[idx] + f_zero = fi.f_zero_frequencies[idx] + + plt.plot(time, freq, color=DATA_COLOR) + plt.plot(time[freq == f_zero][0], f_zero, "o", color="black") + f_inf_time = time[(0.2 < time) & (time < 0.4)] + plt.plot(f_inf_time, [f_inf for _ in f_inf_time], color="black") + plt.xlim((-0.1, 0.6)) + plt.show() + + +def data_mean_freq_step_stimulus_examples(): + # todo smooth! add f_0, f_inf, f_base to it? + cell = "data/invivo/2013-04-17-ac-invivo-1/" + cell_data = CellData(cell) + fi = FICurveCellData(cell_data, cell_data.get_fi_contrasts()) + + time_traces, freq_traces = fi.get_time_and_freq_traces() + mean_times, mean_freqs = fi.get_mean_time_and_freq_traces() + + used_idicies = (0, 7, -1) + fig, axes = plt.subplots(len(used_idicies), figsize=(8, 12), sharex=True, sharey=True) + for ax_idx, idx in enumerate(used_idicies): + sv = fi.stimulus_values[idx] + + # for j in range(len(time_traces[i])): + # axes[i].plot(time_traces[i][j], freq_traces[i][j], color="gray", alpha=0.5) + + axes[ax_idx].plot(mean_times[idx], mean_freqs[idx], color=DATA_COLOR) + # plt.plot(mean_times[i], mean_freqs[i], color="black") + + axes[ax_idx].set_ylabel("Frequency [Hz]") + axes[ax_idx].set_xlim((-0.2, 0.6)) + axes[ax_idx].set_title("Contrast {:.2f} ({:} trials)".format(sv, len(time_traces[idx]))) + axes[ax_idx].set_xlabel("Time [s]") + plt.show() + + +def data_isi_histogram(recalculate=True): + # if isis loadable - load + name = "isi_cell_data.npy" + path = os.path.join(DATA_SAVE_PATH, name) + if os.path.exists(path) and not recalculate: + isis = np.load(path) + print("loaded") + else: + # if not get them from the cell + cell = "data/invivo/2013-04-17-ac-invivo-1/" # not bursty + # cell = "data/invivo/2014-12-03-ad-invivo-1/" # half bursty + # cell = "data/invivo/2015-01-20-ad-invivo-1/" # does triple peaks... + # cell = "data/invivo/2018-05-08-ae-invivo-1/" # a bit bursty + # cell = "data/invivo/2013-04-10-af-invivo-1/" # a bit bursty + cell_data = CellData(cell) + base = BaselineCellData(cell_data) + isis = np.array(base.get_interspike_intervals()) + # base.plot_baseline(position=0,time_length=10) + + # save isis + np.save(path, isis) + isis = isis * 1000 + # plot histogram + bins = np.arange(0, 30.1, 0.1) + plt.hist(isis, bins=bins, color=DATA_COLOR) + plt.xlabel("Inter spike intervals [ms]") + plt.ylabel("Count") + plt.tight_layout() + + plt.show() + + + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/Figures_Model.py b/Figures_Model.py new file mode 100644 index 0000000..c33de81 --- /dev/null +++ b/Figures_Model.py @@ -0,0 +1,56 @@ + +from models.LIFACnoise import LifacNoiseModel +from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus + +import matplotlib.pyplot as plt +import numpy as np + + +def main(): + plot_model_example() + pass + + +def plot_model_example(): + + parameter = {} + model = LifacNoiseModel(parameter) + # frequency, contrast, start_time=0, duration=np.inf, amplitude=1) + frequency = 350 + contrast = 0 + start_time = 5 + duration = 0 + stimulus = SinusoidalStepStimulus(frequency, contrast, start_time, duration) + + time_start = 0 + time_duration = 0.5 + time_step = model.get_sampling_interval() + v1, spikes = model.simulate_fast(stimulus, total_time_s=time_duration, time_start=time_start) + adaption = model.get_adaption_trace() + time = np.arange(time_start, time_start+time_duration, time_step) + fig, axes = plt.subplots(2, sharex=True) + + # axes[0].plot(time, stimulus.as_array(time_start, time_duration, time_step)) + + start = 0.26 + end = 0.29 + start_idx = int(start / time_step) + end_idx = int(end / time_step) + time_part = np.arange(0, end_idx-start_idx, 1) * time_step *1000 + # axes[0].plot(time[start_idx:end_idx], v1[start_idx:end_idx]) + axes[0].plot(time_part, v1[start_idx:end_idx]) + axes[0].set_ylabel("Membrane voltage [mV]") + # axes[1].plot(time[start_idx:end_idx], adaption[start_idx:end_idx]) + axes[1].plot(time_part, adaption[start_idx:end_idx]) + axes[1].set_ylabel("Adaption current [mV]") + axes[1].set_xlabel("Time [ms]") + axes[1].set_xlim((0, 30)) + plt.show() + plt.close() + + + + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/Figures_Stimuli.py b/Figures_Stimuli.py new file mode 100644 index 0000000..aeef33c --- /dev/null +++ b/Figures_Stimuli.py @@ -0,0 +1,67 @@ + +from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus +from stimuli.SinusAmplitudeModulation import SinusAmplitudeModulationStimulus +import numpy as np +import matplotlib.pyplot as plt + + +def main(): + plot_step_stimulus() + plot_sam_stimulus() + pass + + +def plot_step_stimulus(): + start = 0 + end = 1 + + time_start = -0.2 + time_end = 1.2 + step_size = 0.00005 + + frequency = 20 + contrast = 0.5 + # frequency, contrast, start_time=0, duration=np.inf, amplitude=1 + step_stim= SinusoidalStepStimulus(frequency, contrast, start, end-start) + + values = step_stim.as_array(time_start, time_end-time_start, step_size) + time = np.arange(time_start, time_end, step_size) + + plt.plot(time, values) + plt.xlabel("Time [s]") + plt.ylabel("Voltage [mV]") + plt.savefig("thesis/figures/sin_step_stim_example.pdf") + plt.close() + + +def plot_sam_stimulus(): + start = 0 + end = 1 + + time_start = -0.2 + time_end = 1.2 + step_size = 0.00005 + + contrast = 0.5 + mod_freq = 10 + carrier_freq = 53 + # carrier_frequency, contrast, modulation_frequency, start_time=0, duration=np.inf, amplitude=1 + step_stim = SinusAmplitudeModulationStimulus(carrier_freq, contrast, mod_freq, start, end - start) + + values = step_stim.as_array(time_start, time_end - time_start, step_size) + time = np.arange(time_start, time_end, step_size) + plt.plot(time, values) + + beat_time = np.arange(start, end, step_size) + beat_values = np.sin(beat_time*2*np.pi*mod_freq) * contrast + 1 + plt.plot(beat_time, beat_values) + + plt.xlabel("Time [s]") + plt.ylabel("Voltage [mV]") + # plt.show() + plt.savefig("thesis/figures/sam_stim_example.pdf") + plt.close() + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/read_data_infos.py b/read_data_infos.py new file mode 100644 index 0000000..9a1c80b --- /dev/null +++ b/read_data_infos.py @@ -0,0 +1,118 @@ +from CellData import icelldata_of_dir, CellData +import numpy as np +import os +import pyrelacs.DataLoader as Dl + + +def main(): + fish_info() + # eod_info() + # fi_recording_times() + + +def eod_info(): + cells = [] + for item in sorted(os.listdir("data/invivo/")): + cells.append(os.path.join("data/invivo/", item)) + + for item in sorted(os.listdir("data/invivo_bursty/")): + cells.append(os.path.join("data/invivo_bursty/", item)) + + eod_freq = [] + + for cell in cells: + data = CellData(cell) + eod_f = data.get_eod_frequency() + if not np.isnan(eod_f): + eod_freq.append(eod_f) + else: + print(cell) + + print("eod Freq: min {}, max {}, mean: {:.2f}, std: {:.2f}".format(min(eod_freq), max(eod_freq), np.mean(eod_freq), + np.std(eod_freq))) + +def fish_info(): + cells = [] + for item in sorted(os.listdir("data/invivo/")): + cells.append(os.path.join("data/invivo/", item)) + + for item in sorted(os.listdir("data/invivo_bursty/")): + cells.append(os.path.join("data/invivo_bursty/", item)) + cell_type = [] + weight = [] + size = [] + eod_freq = [] + preparation = [] + for cell in cells: + info_file = os.path.join(cell, "info.dat") + for metadata in Dl.load(info_file): + if "CellType" not in metadata.keys(): + cell_type.append(metadata["Cell"]["CellType"]) + if cell_type[-1] != "P-unit": + print("not P-unit?", cell) + if "Weight" in metadata["Subject"].keys(): + weight.append(float(metadata["Subject"]["Weight"][:-1])) + size.append(float(metadata["Subject"]["Size"][:-2])) + if "CellProperties" in metadata.keys(): + eod_freq.append(float(metadata["CellProperties"]["EOD Frequency"][:-2])) + elif "Cell properties" in metadata.keys(): + eod_freq.append(float(metadata["Cell properties"]["EOD Frequency"][:-2])) + preparation.append(metadata["Preparation"]) + # print(metadata) + else: + cell_type.append(metadata["CellType"]) + if cell_type[-1] != "P-unit": + print("not P-unit?", cell) + + weight.append(float(metadata["Weight"][:-1])) + size.append(float(metadata["Size"][:-2])) + # 'LocalAnaesthesia': 'true', 'AnaestheticDose': '120mg/l', 'Anaesthetic': 'MS 222', 'LocalAnaesthetic': 'Lidocaine', 'Anaesthesia': 'true', 'Type': 'in vivo', 'Immobilization': 'Tubocurarin' + eod_freq.append(float(metadata["EOD Frequency"][:-2])) + prep_dict = {} + for key in ('LocalAnaesthesia', 'AnaestheticDose', 'Anaesthetic', 'LocalAnaesthetic', 'Anaesthesia', 'Immobilization'): + prep_dict[key] = metadata[key] + preparation.append(prep_dict) + # print(metadata) + + print("Size: min {}, max {}".format(min(size), max(size))) + print("weight: min {}, max {}".format(min(weight), max(weight))) + print("eod Freq: min {}, max {}, mean: {:.2f}, std: {:.2f}".format(min(eod_freq), max(eod_freq), np.mean(eod_freq), np.std(eod_freq))) + print("anaesthetics:", np.unique([x['Anaesthetic'] for x in preparation])) + print("anaesthetic dosages:", np.unique([x['AnaestheticDose'] for x in preparation])) + print("local anaesthetic:", np.unique([x['LocalAnaesthesia'] for x in preparation])) + print("Immobilization:", np.unique([x['Immobilization'] for x in preparation])) + + +def fi_recording_times(): + + recording_times = [] + for cell_data in icelldata_of_dir("data/invivo/", test_for_v1_trace=False): + # time_start, stimulus_start, stimulus_duration, after_stimulus_duration + recording_times.append(cell_data.get_recording_times()) + for cell_data in icelldata_of_dir("data/invivo_bursty/", test_for_v1_trace=False): + # time_start, stimulus_start, stimulus_duration, after_stimulus_duration + recording_times.append(cell_data.get_recording_times()) + + recording_times = np.array(recording_times) + + time_starts = recording_times[:, 0] + stimulus_starts = recording_times[:, 1] + stimulus_durations = recording_times[:, 2] + after_durations = recording_times[:, 3] + + print("Fi-curve stimulus recording times:") + print("time_starts:", np.unique(time_starts)) + print("stimulus_starts:", np.unique(stimulus_starts)) + unique_durations = np.unique(stimulus_durations) + print("stimulus_durations:", unique_durations) + + for d in unique_durations: + print("cells with stimulus duration {}: {}".format(d, np.sum(stimulus_durations == d))) + + print("after_durations:", np.unique(after_durations)) + + + + +if __name__ == '__main__': + main() \ No newline at end of file