add scripts to create plots and read info from data

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

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

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

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