diff --git a/FiCurve.py b/FiCurve.py index d32e1a7..74a1c6c 100644 --- a/FiCurve.py +++ b/FiCurve.py @@ -2,6 +2,7 @@ from CellData import CellData import numpy as np from scipy.optimize import curve_fit +from scipy.stats import linregress import matplotlib.pyplot as plt from warnings import warn import functions as fu @@ -33,14 +34,24 @@ class FICurve: def all_calculate_frequency_points(self): mean_frequencies = self.cell_data.get_mean_isi_frequencies() + time_axes = self.cell_data.get_time_axes_mean_frequencies() if len(mean_frequencies) == 0: warn("FICurve:all_calculate_frequency_points(): mean_frequencies is empty.\n" "Was all_calculate_mean_isi_frequencies already called?") - for freq in mean_frequencies: - self.f_zeros.append(self.__calculate_f_zero__(freq)) - self.f_baselines.append(self.__calculate_f_baseline__(freq)) - self.f_infinities.append(self.__calculate_f_infinity__(freq)) + for i in range(len(mean_frequencies)): + + if time_axes[i][0] > self.cell_data.get_stimulus_start(): + # TODO + warn("TODO: Deal with to strongly cut frequency traces in cell data! ") + self.f_zeros.append(-1) + self.f_baselines.append(-1) + self.f_infinities.append(-1) + continue + + self.f_zeros.append(self.__calculate_f_zero__(time_axes[i], mean_frequencies[i])) + self.f_baselines.append(self.__calculate_f_baseline__(time_axes[i], mean_frequencies[i])) + self.f_infinities.append(self.__calculate_f_infinity__(time_axes[i], mean_frequencies[i])) def fit_line(self): popt, pcov = curve_fit(fu.clipped_line, self.stimulus_value, self.f_infinities) @@ -48,7 +59,7 @@ class FICurve: def fit_boltzmann(self): max_f0 = float(max(self.f_zeros)) - min_f0 = float(min(self.f_zeros)) + min_f0 = 0.1 # float(min(self.f_zeros)) mean_int = float(np.mean(self.stimulus_value)) total_increase = max_f0 - min_f0 @@ -57,78 +68,93 @@ class FICurve: popt, pcov = curve_fit(fu.full_boltzmann, self.stimulus_value, self.f_zeros, p0=(max_f0, min_f0, start_k, mean_int), - maxfev=10000, bounds=([0, 0, -np.inf, -np.inf], [3000, 3000, np.inf, np.inf])) + maxfev=10000, bounds=([0, 0, -np.inf, -np.inf], [5000, 1, np.inf, np.inf])) self.boltzmann_fit_vars = popt - def plot_fi_curve(self, savepath: str = None): - min_x = min(self.stimulus_value) - max_x = max(self.stimulus_value) - step = (max_x - min_x) / 5000 - x_values = np.arange(min_x, max_x, step) + def __calculate_f_baseline__(self, time, frequency, buffer=0.025): - plt.plot(self.stimulus_value, self.f_baselines, color='blue', label='f_base') + stim_start = self.cell_data.get_stimulus_start() - time[0] + sampling_interval = self.cell_data.get_sampling_interval() + if stim_start < 0.1: + warn("FICurve:__calculate_f_baseline__(): Quite short delay at the start.") - plt.plot(self.stimulus_value, self.f_infinities, 'o', color='lime', label='f_inf') - plt.plot(x_values, [fu.clipped_line(x, self.f_infinity_fit[0], self.f_infinity_fit[1]) for x in x_values], - color='darkgreen', label='f_inf_fit') + start_idx = 0 + end_idx = int((stim_start-buffer)/sampling_interval) + f_baseline = np.mean(frequency[start_idx:end_idx]) - plt.plot(self.stimulus_value, self.f_zeros, 'o', color='orange', label='f_zero') - popt = self.boltzmann_fit_vars - plt.plot(x_values, [fu.full_boltzmann(x, popt[0], popt[1], popt[2], popt[3]) for x in x_values], - color='red', label='f_0_fit') + plt.plot((start_idx, end_idx), (f_baseline, f_baseline), label="f_baseline") + return f_baseline - plt.legend() - plt.ylabel("Frequency [Hz]") - if self.using_contrast: - plt.xlabel("Stimulus contrast") - else: - plt.xlabel("Stimulus intensity [mv]") - if savepath is None: - plt.show() - else: - plt.savefig(savepath + "fi_curve.png") - plt.close() + def __calculate_f_zero__(self, time, frequency, peak_buffer_percent=0.05, buffer=0.025): - def __calculate_f_baseline__(self, frequency, buffer=0.025): - delay = self.cell_data.get_delay() + stimulus_start = self.cell_data.get_stimulus_start() - time[0] # time start is generally != 0 and != delay sampling_interval = self.cell_data.get_sampling_interval() - if delay < 0.1: - warn("FICurve:__calculate_f_baseline__(): Quite short delay at the start.") - idx_start = int(buffer/sampling_interval) - idx_end = int((delay-buffer)/sampling_interval) - return np.mean(frequency[idx_start:idx_end]) + freq_before = frequency[0:int((stimulus_start - buffer) / sampling_interval)] + min_before = min(freq_before) + max_before = max(freq_before) + mean_before = np.mean(freq_before) - def __calculate_f_zero__(self, frequency, length_of_mean=0.1, buffer=0.025): - stimulus_start = self.cell_data.get_delay() + self.cell_data.get_stimulus_start() - sampling_interval = self.cell_data.get_sampling_interval() + # time where the f-zero is searched in + start_idx = int((stimulus_start-0.1*buffer) / sampling_interval) + end_idx = int((stimulus_start + buffer) / sampling_interval) - start_idx = int((stimulus_start - buffer) / sampling_interval) - end_idx = int((stimulus_start + buffer*2) / sampling_interval) + min_during_start_of_stim = min(frequency[start_idx:end_idx]) + max_during_start_of_stim = max(frequency[start_idx:end_idx]) - freq_before = frequency[start_idx-(int(length_of_mean/sampling_interval)):start_idx] - fb_mean = np.mean(freq_before) - fb_std = np.std(freq_before) + if abs(mean_before-min_during_start_of_stim) > abs(max_during_start_of_stim-mean_before): + f_zero = min_during_start_of_stim + else: + f_zero = max_during_start_of_stim + + peak_buffer = (max_before - min_before) * peak_buffer_percent + if min_before - peak_buffer <= f_zero <= max_before + peak_buffer: + end_idx = start_idx + int((end_idx-start_idx)/2) + f_zero = np.mean(frequency[start_idx:end_idx]) + + plt.plot(frequency) + plt.plot((start_idx, end_idx), (f_zero, f_zero), label="f_zero, {:.2f}".format(peak_buffer)) + + return f_zero + + # start_idx = int(stimulus_start / sampling_interval) + # end_idx = int((stimulus_start + buffer*2) / sampling_interval) + # + # freq_before = frequency[start_idx-(int(length_of_mean/sampling_interval)):start_idx] + # fb_mean = np.mean(freq_before) + # fb_std = np.std(freq_before) + # + # peak_frequency = fb_mean + # count = 0 + # for i in range(start_idx + 1, end_idx): + # if fb_mean-3*fb_std <= frequency[i] <= fb_mean+3*fb_std: + # continue + # + # if abs(frequency[i] - fb_mean) > abs(peak_frequency - fb_mean): + # peak_frequency = frequency[i] + # count += 1 + + # return peak_frequency + + def __calculate_f_infinity__(self, time, frequency, length=0.1, buffer=0.025): + stimulus_end_time = self.cell_data.get_stimulus_start() + self.cell_data.get_stimulus_duration() - time[0] - peak_frequency = fb_mean - count = 0 - for i in range(start_idx + 1, end_idx): - if fb_mean-3*fb_std <= frequency[i] <= fb_mean+3*fb_std: - continue + start_idx = int((stimulus_end_time - length - buffer) / self.cell_data.get_sampling_interval()) + end_idx = int((stimulus_end_time - buffer) / self.cell_data.get_sampling_interval()) - if abs(frequency[i] - fb_mean) > abs(peak_frequency - fb_mean): - peak_frequency = frequency[i] - count += 1 + x = np.arange(start_idx, end_idx, 1) # time[start_idx:end_idx] + slope, intercept, r_value, p_value, std_err = linregress(x, frequency[start_idx:end_idx]) - return peak_frequency + if p_value < 0.0001: + plt.title("significant slope: {:.2f}, p: {:.5f}, r: {:.5f}".format(slope, p_value, r_value)) + plt.plot(x, [i*slope + intercept for i in x], color="black") - def __calculate_f_infinity__(self, frequency, length=0.2, buffer=0.025): - stimulus_end_time = \ - self.cell_data.get_delay() + self.cell_data.get_stimulus_start() + self.cell_data.get_stimulus_duration() - start_idx = int((stimulus_end_time - length - buffer) / self.cell_data.get_sampling_interval()) - end_idx = int((stimulus_end_time - buffer) / self.cell_data.get_sampling_interval()) + plt.plot((start_idx, end_idx), (np.mean(frequency[start_idx:end_idx]), np.mean(frequency[start_idx:end_idx])), label="f_inf") + plt.legend() + plt.show() + plt.close() return np.mean(frequency[start_idx:end_idx]) @@ -150,3 +176,48 @@ class FICurve: def get_fi_curve_slope_of_straight(self): fit_vars = self.boltzmann_fit_vars return fu.full_boltzmann_straight_slope(fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3]) + + def plot_fi_curve(self, savepath: str = None): + min_x = min(self.stimulus_value) + max_x = max(self.stimulus_value) + step = (max_x - min_x) / 5000 + x_values = np.arange(min_x, max_x, step) + + plt.plot(self.stimulus_value, self.f_baselines, color='blue', label='f_base') + + plt.plot(self.stimulus_value, self.f_infinities, 'o', color='lime', label='f_inf') + plt.plot(x_values, [fu.clipped_line(x, self.f_infinity_fit[0], self.f_infinity_fit[1]) for x in x_values], + color='darkgreen', label='f_inf_fit') + + plt.plot(self.stimulus_value, self.f_zeros, 'o', color='orange', label='f_zero') + popt = self.boltzmann_fit_vars + plt.plot(x_values, [fu.full_boltzmann(x, popt[0], popt[1], popt[2], popt[3]) for x in x_values], + color='red', label='f_0_fit') + + plt.legend() + plt.ylabel("Frequency [Hz]") + if self.using_contrast: + plt.xlabel("Stimulus contrast") + else: + plt.xlabel("Stimulus intensity [mv]") + if savepath is None: + plt.show() + else: + plt.savefig(savepath + "fi_curve.png") + plt.close() + + def plot_f_point_detections(self): + mean_frequencies = np.array(self.cell_data.get_mean_isi_frequencies()) + time_axes = self.cell_data.get_time_axes_mean_frequencies() + + + for i in range(len(mean_frequencies)): + fig, axes = plt.subplots(1, 1, sharex="all") + axes.plot(time_axes[i], mean_frequencies[i], label="voltage") + axes.plot((time_axes[i][0],time_axes[i][-1]), (self.f_zeros[i], self.f_zeros[i]), label="f_zero") + axes.plot((time_axes[i][0],time_axes[i][-1]), (self.f_infinities[i], self.f_infinities[i]), '--', label="f_inf") + axes.plot((time_axes[i][0],time_axes[i][-1]), (self.f_baselines[i], self.f_baselines[i]), label="f_base") + axes.set_title(str(self.stimulus_value[i])) + plt.legend() + + plt.show()