246 lines
10 KiB
Python
246 lines
10 KiB
Python
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from IPython import embed
|
|
from model import simulate, load_models
|
|
import matplotlib.gridspec as gridspec
|
|
from plot_eod_chirp import power_parameters
|
|
from scipy.ndimage import gaussian_filter
|
|
import matplotlib.mlab as ml
|
|
import scipy.integrate as si
|
|
import pandas as pd
|
|
"""
|
|
Dependencies:
|
|
numpy
|
|
matplotlib
|
|
numba (optional, speeds simulation up: pre-compiles functions to machine code)
|
|
"""
|
|
|
|
|
|
def main():
|
|
# tiny example program:
|
|
|
|
example_cell_idx = 11
|
|
|
|
# load model parameter:
|
|
parameters = load_models("models.csv")
|
|
|
|
counter = 0
|
|
beat_results = [[]] * (3 * len(parameters))
|
|
beat_results = []
|
|
for d in range(len(parameters)):
|
|
model_params = parameters[d]
|
|
cell = model_params.pop('cell')
|
|
eod_fr = model_params.pop('EODf')
|
|
print("Example with cell:", cell)
|
|
step = 20
|
|
eod_fe = np.arange(0, eod_fr * 5, step)
|
|
# generate EOD-like stimulus with an amplitude step:
|
|
deltat = model_params["deltat"]
|
|
stimulus_length = 11.0 # in seconds
|
|
#time = np.arange(0, stimulus_length, deltat)
|
|
time = np.arange(0, stimulus_length, deltat)
|
|
# baseline EOD with amplitude 1:
|
|
a_fr = 1 # amplitude fish reciever
|
|
a_fe = 0.2 # amplitude fish emitter
|
|
# results = [[]]*
|
|
p_array = [[]]*len(eod_fe)
|
|
for e in range(len(eod_fe)):
|
|
f_corr = create_beat_corr(np.array([eod_fe[e]]), np.array([eod_fr]))
|
|
time_fish_r = time * 2 * np.pi * eod_fr
|
|
eod_fish_r = a_fr * np.sin(time_fish_r)
|
|
time_fish_e = time * 2 * np.pi * eod_fe[e]
|
|
eod_fish_e = a_fe * np.sin(time_fish_e)
|
|
stimulus = eod_fish_e+eod_fish_r
|
|
|
|
# integrate the model:
|
|
spikes = simulate(stimulus, **model_params)
|
|
spikes_new = spikes[spikes > 1]
|
|
sampling_rate = 1/deltat
|
|
spikes_mat = np.zeros(int(spikes_new[-1] * sampling_rate) + 2)
|
|
spikes_idx = np.round((spikes_new) * sampling_rate)
|
|
for spike in spikes_idx:
|
|
spikes_mat[int(spike)] = 1
|
|
spikes_mat = spikes_mat*sampling_rate
|
|
#window005 = 0.00005 * sampling_rate
|
|
window05 = 0.0005 * sampling_rate
|
|
window2 = 0.002 * sampling_rate
|
|
#smoothened_spikes_mat005 = gaussian_filter(spikes_mat, sigma=window005)
|
|
smoothened_spikes_mat05 = gaussian_filter(spikes_mat, sigma=window05)
|
|
smoothened_spikes_mat2 = gaussian_filter(spikes_mat, sigma=window2)
|
|
nfft = 4096*2
|
|
freq_step = sampling_rate / nfft
|
|
|
|
name = 'threshold'
|
|
p_array, f = ml.psd(spikes_mat - np.mean(spikes_mat), Fs=sampling_rate, NFFT=nfft,
|
|
noverlap=nfft / 2) #
|
|
p05_array, f = ml.psd(smoothened_spikes_mat05 - np.mean(smoothened_spikes_mat05), Fs=sampling_rate, NFFT=nfft,
|
|
noverlap=nfft / 2) #
|
|
p2_array, f = ml.psd(smoothened_spikes_mat2 - np.mean(smoothened_spikes_mat2), Fs=sampling_rate, NFFT=nfft,
|
|
noverlap=nfft / 2) #
|
|
|
|
beat_results.append({})
|
|
#embed()
|
|
#p_array = np.nanmean(p_array, axis=1)
|
|
beat_results[-1]['result_frequency_original'] = f[np.argmax(p_array[f < 0.5 * eod_fr])] # f ist immer gleich
|
|
beat_results[-1]['result_amplitude_original'] = np.sqrt(si.trapz(p_array, f, freq_step))
|
|
beat_results[-1]['result_amplitude_max_original'] = np.sqrt(
|
|
np.max(p_array[f < 0.5 * eod_fr]) * freq_step)
|
|
beat_results[-1]['result_toblerone_max_original'] = np.sqrt(
|
|
(p_array[np.argmin(np.abs(f - f_corr[0]))]) * freq_step)
|
|
|
|
beat_results[-1]['result_frequency_whole'] = f[np.argmax(p_array)] # f ist immer gleich
|
|
beat_results[-1]['result_amplitude_whole'] = np.sqrt(si.trapz(p_array, f, freq_step))
|
|
beat_results[-1]['result_amplitude_max_whole'] = np.sqrt(np.max(p_array) * freq_step)
|
|
beat_results[-1]['result_toblerone_max_whole'] = np.sqrt(
|
|
(p_array[np.argmin(np.abs(f - f_corr[0]))]) * freq_step)
|
|
#if b != 'no':
|
|
# beat_results[-1]['result_baseline_max_whole'] = np.sqrt(
|
|
# (p_array[np.argmin(np.abs(f - b))]) * freq_step)
|
|
# beat_results[-1]['result_baseline_max_original'] = np.sqrt(
|
|
# (p_array[np.argmin(np.abs(f - b))]) * freq_step)
|
|
# arrays = [p05_array, p005_array, p2_array, p_inst_array]
|
|
# titles = ['05', '005', '2', 'inst']
|
|
#array = [smoothened_spikes_mat05,smoothened_spikes_mat2]
|
|
arrays = [p05_array, p2_array]
|
|
#arrays = [[]]*len(array)
|
|
titles = ['05','2']
|
|
for a in range(len(arrays)):
|
|
|
|
beat_results[-1]['result_frequency_' + titles[a]] = f[np.argmax(arrays[a])]
|
|
beat_results[-1]['result_frequency_half' + titles[a]] = f[np.argmax(arrays[a][f < 0.5 * eod_fr])]
|
|
beat_results[-1]['result_amplitude_' + titles[a]] = np.sqrt(si.trapz(arrays[a], f, freq_step))
|
|
beat_results[-1]['result_amplitude_max_half' + titles[a]] = np.sqrt(
|
|
np.max(arrays[a][f < 0.5 * eod_fr]) * freq_step)
|
|
beat_results[-1]['result_amplitude_max_' + titles[a]] = np.sqrt(np.max(arrays[a]) * freq_step)
|
|
beat_results[-1]['result_toblerone_max_' + titles[a]] = np.sqrt(
|
|
arrays[a][np.argmin(np.abs(f - f_corr[0]))] * freq_step)
|
|
#if b != 'no':
|
|
# beat_results[-1]['result_baseline_max_' + titles[a]] = np.sqrt(
|
|
# arrays[a][np.argmin(np.abs(f - b))] * freq_step)
|
|
|
|
# try:
|
|
# beat_results[-1],p,f = power_parameters(beat_results[-1], array[i], 1/deltat, nfft, name, eod_fr,title = titles[i],max_all = 'result_amplitude_max_',max = 'result_amplitude_max_half',integral = 'result_amplitude_',plot = False,freq_half = 'result_frequency_half' ,freq_whole = 'result_frequency_' )
|
|
# except:
|
|
# embed()
|
|
beat_results[-1]['dataset'] = cell
|
|
beat_results[-1]['eodf'] = eod_fr
|
|
#beat_results[counter]['baseline_fr'] = fr
|
|
beat_results[-1]['beat_corr'] = f_corr
|
|
beat_results[-1]['delta_f'] = eod_fe[e] -eod_fr
|
|
#if key[0] == 50:
|
|
# embed()
|
|
beat_results[-1]['contrast'] = a_fe
|
|
beat_results[-1]['type'] = name
|
|
#beat_results[-1]['f_original'] = max_peak
|
|
beat_results[-1]['f'] = f
|
|
#counter += 1
|
|
#embed()
|
|
df_beat = pd.DataFrame(beat_results)
|
|
df_beat = df_beat.dropna()
|
|
df_beat.to_pickle('modell_all_cell.pkl')
|
|
df_beat.to_csv('modell_all_cell.csv')
|
|
np.save('modell_all_cell.npy', df_beat)
|
|
embed()
|
|
ax = {}
|
|
for i in range(len(results)):
|
|
ax[i] = plt.subplot(2,3,i+1)
|
|
plt.plot((eod_fe -eod_fr)/(eod_fr)+1,results[i]['f'],color = 'red')
|
|
#plt.scatter((eod_fe - eod_fr) / (eod_fr) + 1, results[i]['f'],color = 'red')
|
|
ax[0].set_ylabel('[Hz]')
|
|
ax[i].set_ylim([0,eod_fr/2])
|
|
for i in range(len(results)):
|
|
ax[i+len(results)] = plt.subplot(2,3,i+len(results)+1)
|
|
plt.plot((eod_fe -eod_fr)/(eod_fr)+1,results[i]['max'],color = 'steelblue')
|
|
#plt.scatter((eod_fe - eod_fr) / (eod_fr) + 1, results[i]['max'],color = 'blue')
|
|
ax[len(results)].set_ylabel('Modulation')
|
|
ax[len(results)+i].set_xlabel('EOD multiples')
|
|
plt.subplots_adjust(wspace = 0.3)
|
|
plt.savefig('modell_single_cell.pdf')
|
|
plt.savefig('../highbeats_pdf/cell_simulations/modell_single_cell'+cell+'.pdf')
|
|
plt.show()
|
|
embed()
|
|
# some analysis an dplotting:
|
|
#embed()
|
|
grid = gridspec.GridSpec(int(np.sqrt(len(parameters))), int(np.ceil(np.sqrt(len(parameters)))))
|
|
parameters = load_models("models.csv")
|
|
for i in range(4):
|
|
#for i in range(len(parameters)):
|
|
model_params = parameters[i]
|
|
print(cell)
|
|
cell = model_params.pop('cell')
|
|
EODf = model_params.pop('EODf')
|
|
|
|
|
|
# generate EOD-like stimulus
|
|
deltat = model_params["deltat"]
|
|
stimulus_length = 11.0 # in seconds
|
|
time = np.arange(0, stimulus_length, deltat)
|
|
# baseline EOD with amplitude 1:
|
|
stimulus = np.sin(2 * np.pi * EODf * time)
|
|
# das lasse ich eine sekunde integrieren dann weitere 10 sekunden integrieren und das nehmen
|
|
spikes = simulate(stimulus, **model_params)
|
|
|
|
# cut off first second of response
|
|
new_spikes = spikes[spikes >1]
|
|
|
|
freq,isis = calculate_isi_frequency(new_spikes, deltat)
|
|
|
|
#embed()
|
|
plt.subplot(grid[i])
|
|
plt.title('B:'+np.mean(freq))
|
|
plt.hist(isis, bins = 100, density = True)
|
|
plt.savefig('isi_model.pdf')
|
|
plt.savefig('../highbeats_pdf/isi_model.pdf')
|
|
plt.show()
|
|
|
|
freq_time = np.arange(spikes[0], spikes[-1], deltat)
|
|
|
|
fig, axs = plt.subplots(2, 1, sharex="col")
|
|
|
|
axs[0].plot(time, stimulus)
|
|
axs[0].set_title("Stimulus")
|
|
axs[0].set_ylabel("Amplitude in mV")
|
|
|
|
axs[1].plot(freq_time, freq)
|
|
axs[1].set_title("Model Frequency")
|
|
axs[1].set_ylabel("Frequency in Hz")
|
|
axs[1].set_xlabel("Time in s")
|
|
plt.show()
|
|
plt.close()
|
|
|
|
def create_beat_corr(hz_range, eod_fr):
|
|
beat_corr = hz_range%eod_fr
|
|
beat_corr[beat_corr>eod_fr/2] = eod_fr[beat_corr>eod_fr/2] - beat_corr[beat_corr>eod_fr/2]
|
|
return beat_corr
|
|
|
|
def calculate_isi_frequency(spikes, deltat):
|
|
"""
|
|
calculates inter-spike interval frequency
|
|
(wasn't tested a lot may give different length than time = np.arange(spikes[0], spikes[-1], deltat),
|
|
or raise an index error for some inputs)
|
|
|
|
:param spikes: spike time points
|
|
:param deltat: integration time step of the model
|
|
|
|
:return: the frequency trace:
|
|
starts at the time of first spike
|
|
ends at the time of the last spike.
|
|
"""
|
|
|
|
isis = np.diff(spikes)
|
|
freq_points = 1 / isis
|
|
freq = np.zeros(int((spikes[-1] - spikes[0]) / deltat))
|
|
|
|
current_idx = 0
|
|
for i, isi in enumerate(isis):
|
|
end_idx = int(current_idx + np.rint(isi / deltat))
|
|
freq[current_idx:end_idx] = freq_points[i]
|
|
current_idx = end_idx
|
|
|
|
return freq,isis
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main()
|