diff --git a/params.py b/params.py new file mode 100644 index 0000000..531dcfb --- /dev/null +++ b/params.py @@ -0,0 +1,134 @@ +import matplotlib.colors as mcolors +import matplotlib as mpl +import matplotlib.dates as mdates +import matplotlib.pyplot as plt +import datetime +import numpy as np +from IPython import embed + +############################################################################################# +# default fontsize +plt.rcParams['font.size'] = 11 +plt.rcParams['axes.linewidth'] = 2 + +############################################################################################# +# for boxplots +props = dict(linewidth=2, color='black') +props_a = dict(linewidth=2, color='#884ea0') +props_e = dict(linewidth=2, color='#656565') +flierprops = dict(marker='o', markeredgecolor='#656565', markersize=6, linestyle='none') + +############################################################################################# +# time boundaries of day, dusk, night and dawn +time_bound = ['1900_01_01_6_30', '1900_01_01_16_30', '1900_01_01_18_30', + '1900_01_02_4_30', '1900_01_02_6_30', '1900_01_02_16_30'] +time_boxes = [] +for i in range(len(time_bound)): + time_boxes.append(datetime.datetime.strptime(time_bound[i], '%Y_%m_%d_%H_%M')) +datetime_box = mdates.date2num(time_boxes) + +# bin length for activity calculations in second +bin_len = 5 +############################################################################################# +# colors +color_efm = ['#2e4053', '#ab1717', '#004d8d', '#000000'] +color_diffdays = ['#1b2631', '#2e4053', '#5d6d7e', '#aeb6bf'] + +cd = ['#566573', '#2c3e50', '#abb2b9', '#2471a3', '#1a5276', '#a9cce3'] +cm20c = plt.cm.get_cmap('tab20c') +color2 = ['#f39c12', '#d35400', '#fad7a0', '#e59866', '#c0392b', '#d98880', '#f1c40f'] +color_dadunida = ['#DFD13A', '#c0392b', '#62285B', '#CE7227'] + +farben = list(mcolors.CSS4_COLORS) +color_vec = [] +for i in range(len(farben)): + color_vec.append(mcolors.CSS4_COLORS[farben[i]]) +color_vec = sorted(color_vec) + +############################################################################################# +# labels +labels = ['$\it{Eigenmannia\,sp.}$', '$\it{A.\,macrostomus}\,\u2640$', '$\it{A.\,macrostomus}\,\u2642$', 'unknown'] + +############################################################################################# +# plot params +fs = 11 # font_size +fs_ticks = 11 # fs_ticks +lw = 2 # linewidth +a = 0.3 # alpha +ms = 7 +h_bins = 30 # hist_bins +save_path = '../../thesis/Figures/Results/' +save_path_pres = '../../pres_ergebnisse/figures/' + +inch = 2.45 + +############################################################################################# +# for nice plots +def nice_ax(ax): + ''' + Makes nice x and y axis: + top and right axis invisible and axis width = 2 and axis ticks labelsize = 11 + :param ax: + :return: + ''' + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.tick_params(width=2) + ax.tick_params(axis='both', which='major', labelsize=11) + + +mpl.axes.Axes.make_nice_ax = nice_ax + + +def beau_time_channel_axis(ax): + ''' + Makes nice axis for trajectory plots: + top and right axis invisible and axis width = 2 and axis ticks labelsize = 11 + white line for the gap: y = 7.5 + y ticks and yticklabels: ticks at each electrode, labels at every second + x and y labels + :param ax: + :return: + ''' + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.tick_params(width=2) + ax.tick_params(axis='both', which='major', labelsize=11) + # ax.set_ylim([0, 15]) + ax.axhline(7.5, xmin=0, xmax=15, color='white', lw=4) + # ax.set_yticklabels([1, 3, 5, 7, 14, 16, 18, 20]) + ax.set_yticks([0, 1, 2, 3, 4, 5, 6, 7, 7.5, 8, 9, 10, 11, 12, 13, 14, 15]) + ax.set_yticklabels(['1', '', '3', '', '5', '', '7', '', 'gap', '', '10', '', '12', '', '14', '', '16']) + # ax.set_yticklabels([0, 0, 2, 4, 6, 13, 15, 17, 19]) + ax.invert_yaxis() + ax.set_xlabel('Time', fontsize=11) + ax.set_ylabel('Electrode', fontsize=11) + + +mpl.axes.Axes.beautimechannelaxis = beau_time_channel_axis + + +def timeaxis(ax): + ''' + Makes nice time axis for datetime x values: + top and right axis invisible and axis width = 2 and axis ticks labelsize = 11 + converts the x ticks in the date format into a real time axis of the form H:M + :param ax: + :return: + ''' + + ax.spines["top"].set_visible(False) + ax.spines["right"].set_visible(False) + ax.tick_params(width=2) + ax.tick_params(axis='both', which='major', labelsize=11) + ax.xaxis_date() + date_format = mdates.DateFormatter('%H:%M') + # date_format = mdates.DateFormatter('%d.%m.%Y %H:%M') + ax.xaxis.set_major_formatter(date_format) + + +mpl.axes.Axes.timeaxis = timeaxis + + + + diff --git a/plot_activity_freq.py b/plot_activity_freq.py new file mode 100644 index 0000000..b5efb6e --- /dev/null +++ b/plot_activity_freq.py @@ -0,0 +1,132 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +import matplotlib.colors as mcolors +import matplotlib.gridspec as gridspec +from IPython import embed +from scipy import stats +import os +from IPython import embed +from params import * +import datetime +import itertools +from statisitic_functions import * + + +if __name__ == '__main__': + ################################################################################################################### + # parameter + # save_path = '../../thesis/Figures/Results/' + # inch = 2.45 + cat_v1 = [0, 0, 750, 0] + cat_v2 = [750, 750, 1000, 1000] + cat_v = [0, 580, 750, 1000] + cat_n = ['Eigenmannia', 'Apteronotus', 'Apteronotus'] + pos = [1, 2, 3, 4] + + ################################################################################################################### + # load data + freq = np.load('../data/f.npy', allow_pickle=True) + c = np.load('../data/all_changes.npy', allow_pickle=True) + names = np.load('../data/n.npy', allow_pickle=True) + stl = np.load('../data/stl.npy', allow_pickle=True) + + ################################################################################################################### + # figure + fig = plt.figure(figsize=[16 / inch, 8 / inch]) + spec = gridspec.GridSpec(ncols=3, nrows=1, figure=fig, hspace=0.3, wspace=0.2, + width_ratios=[4,1,2], left=0.08, bottom=0.15, right=0.95, top=0.90) + ax = fig.add_subplot(spec[0, 0]) + ax1 = fig.add_subplot(spec[0, 2]) + ax2 = fig.add_subplot(spec[0, 1]) + + ############################################################################################################### + # analysis: separates the activity and frequencies by the species and sex + activitys = [] + freqs = [] + for p in range(3): + f = freq[(names == cat_n[p]) & (freq >= cat_v1[p]) & (freq < cat_v2[p])] + a = c[(names == cat_n[p]) & (freq >= cat_v1[p]) & (freq < cat_v2[p])] + + # f = freq[(names == cat_n[p]) & (freq >= cat_v1[p]) & (freq < cat_v2[p]) & (stl==False)] + # a = c[(names == cat_n[p]) & (freq >= cat_v1[p]) & (freq < cat_v2[p]) & (stl==False)] + # + # f_m = freq[(names == cat_n[p]) & (freq >= cat_v1[p]) & (freq < cat_v2[p]) & stl] + # a_m = c[(names == cat_n[p]) & (freq >= cat_v1[p]) & (freq < cat_v2[p]) & stl] + + activitys.append(a) + freqs.append(f) + + ############################################################################################################### + # point plot of freq vs activity + ax.plot(f, a, 'o', alpha=0.7, color=color_efm[p], label=labels[p]) + # ax.plot(f_m, a_m, 'o', alpha=0.7, color=color_efm[p], markeredgecolor='k') + ax.set_xlim(400, 1000) + ax.set_xticks([400, 600, 800, 1000]) + ax.set_xticklabels([400, 600, 800, 1000]) + ax.set_xlabel('EODf [Hz]', fontsize=fs) + ax.set_ylabel('Activity [Changes/min]', fontsize=fs) + ax.set_ylim(-0.2, 6) + + ############################################################################################################### + # boxplot + bp = ax1.boxplot(activitys, widths=0.6, showfliers=False) + for element in ['boxes', 'medians']: + for idx in range(3): + bp[element][idx].set(linewidth=2, color=color_efm[idx]) + + for element in ['whiskers', 'caps']: + for idx, num in enumerate([0,2,4]): + for i in [0,1]: + bp[element][num+i].set(linewidth=2, color=color_efm[idx]) + if element == 'whiskers' and i == 1: + ax1.text(idx+1, bp['whiskers'][num+i].get_ydata()[1]+0.3, str(len(activitys[idx])), + horizontalalignment='center') + + ############################################################################################################### + # histogram: activity distribution + for j in [1,2,0]: + n, bins = np.histogram(activitys[j], bins=np.arange(0,8.2,0.2)) + ax2.hist(activitys[j], bins=np.arange(0,8.2,0.25), orientation='horizontal', histtype='step', + color=color_efm[j], linewidth=2) + + ################################################################################################################ + # statisitics: cohens d' and mann whitney u + + for subset in itertools.combinations([0,1,2], 2): + print(labels[subset[0]], labels[subset[1]]) + U, p = stats.mannwhitneyu(activitys[subset[0]], activitys[subset[1]]) + print('U = ', U, 'p = ', p*3) + r = 1-(2*U)/(len(activitys[subset[0]])*len(activitys[subset[1]])) + print('r = ',r) + print('d = ',cohen_d(activitys[subset[0]], activitys[subset[1]])) + significance_bar(ax1, p*3, None, subset[0]+1, subset[1]+1, 3.8) + + ################################################################################################################# + # nice axis + tagx = [-0.15, -0.4, -0.17] + tagy = [1.05,1.05,1.05] + for ax_idx, axis in enumerate([ax,ax2,ax1]): + axis.text(tagx[ax_idx], tagy[ax_idx], chr(ord('A') + ax_idx), transform=axis.transAxes, fontsize='large') + axis.make_nice_ax() + + ax1.set_xlim(0.5,3.5) + ax1.set_xticks([1, 2, 3]) + ax1.set_xticklabels(['$\it{E}$', '$\it{A}\u2640$', '$\it{A}\u2642$']) + ax1.set_xlabel('Species', fontsize=fs) + + ax1.set_ylim(ax.get_ylim()) + ax1.spines["left"].set_visible(False) + ax1.get_yaxis().set_visible(False) + + ax2.set_ylim(ax.get_ylim()) + ax2.set_yticklabels('') + ax2.set_xlabel('n', fontsize=fs) + + fig.align_xlabels([ax,ax1,ax2]) + fig.legend(loc='upper right', bbox_to_anchor=(0.9, 0.85, 0.1, 0.1)) + + # fig.savefig(save_path_pres+'activ_freq.pdf') + # fig.savefig(save_path+'activ_freq.pdf') + plt.show() + embed() diff --git a/plot_activity_time.py b/plot_activity_time.py new file mode 100644 index 0000000..9014740 --- /dev/null +++ b/plot_activity_time.py @@ -0,0 +1,181 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +import matplotlib.colors as mcolors +import matplotlib.gridspec as gridspec +from IPython import embed +from scipy import stats +import math +import os +from IPython import embed +import helper_functions as hf +from params import * +from statisitic_functions import * +import itertools +import random + + +def get_recording_number_in_time_bins(time_bins): + """ + Calculates the number of the recordings in the time bins + + :param time_bins: numpy array with borders of the time bins + :return: time_bins_recording: numpy array with the number of recordings to that specific time bin + """ + # variables + time_bins_recordings = np.zeros(len(time_bins) - 1) + + # load data + for index, filename_idx in enumerate([0, 1, 2, 3]): + filename = sorted(os.listdir('../data/'))[filename_idx] + time_points = np.load('../data/' + filename + '/all_hms.npy', allow_pickle=True) + + # in which bins is this recording, fill time_bins_recordings + unique_time_points = np.unique(np.hstack(time_points)) + for idx, tb in enumerate(time_bins[:-1]): + if np.any((unique_time_points >= tb) & (unique_time_points <= time_bins[idx + 1])): + time_bins_recordings[idx] += 1 + + return time_bins_recordings + + +if __name__ == '__main__': + ################################################################################################################### + # parameter and variables + # plot params + c = 0 + cc = 1 + + all_changes = [] + all_first = [] + all_last = [] + + # time_bins 5 min + time_factor = 60 * 60 + time_bins = np.arange(0, 24 * time_factor + 1, bin_len) + # time_bins_recordings = np.zeros(len(time_bins) - 1) + + # percent roaming + pr = [] + + # kernel + kernel_size = 120 + kernel = np.ones(kernel_size) / kernel_size + + ################################################################################################################### + # load data + ################################################################################################################### + # load all the data of one day + cbf = np.load('../data/cbf5.npy', allow_pickle=True) + # embed() + # quit() + N_rec_time_bins = get_recording_number_in_time_bins(time_bins) + + ################################################################################################################### + # time zones + time_edges = np.array([4.5, 6.5, 16.5, 18.5]) * time_factor + day = (time_bins[:-1] >= time_edges[1]) & (time_bins[:-1] <= time_edges[2]) + dusk = (time_bins[:-1] >= time_edges[2]) & (time_bins[:-1] <= time_edges[3]) + night = (time_bins[:-1] <= time_edges[0]) | (time_bins[:-1] >= time_edges[3]) + dawn = (time_bins[:-1] >= time_edges[0]) & (time_bins[:-1] <= time_edges[1]) + + ################################################################################################################### + # activity in time bins + activity = np.full(len(time_bins) - 1, np.nan) + + for i in range(len(activity)): + activity[i] = np.nansum(cbf[:2, i]) # * (60/bin_len) + + # activity for boxplot in time zones + activ_boxes = [[], [], [], []] + activ_boxes_corr = [[], [], [], []] + activC = activity / N_rec_time_bins + conv_activ = np.convolve(activity, kernel, 'same') + conv_N_rec_tb = np.convolve(N_rec_time_bins, kernel, 'same') + + activ_boxes[2].extend(activC[night]) # 18:30 - 04:30 + activ_boxes[1].extend(activC[dawn]) # 04:30 - 06:30 + activ_boxes[0].extend(activC[day]) # 06:30 - 16:30 + activ_boxes[3].extend(activC[dusk]) # 16:30 - 18:30 + + # correction if nans in array + for j in range(4): + b = np.array(activ_boxes[j]) + activ_boxes_corr[j] = b[~np.isnan(b)] + t_b = np.array(time_bins[:-1][day]) + t_b = t_b[~np.isnan(np.array(activ_boxes[0]))] + + # rolled time axis for nicer plot midnight in the middle start noon + rolled_time = np.roll(time_bins[:-1] / time_factor, int(len(time_bins[:-1]) / 2)) + rolled_activ = np.roll(activC, int(len(time_bins[:-1]) / 2)) + rolled_conv_activ = np.roll(conv_activ / conv_N_rec_tb, int(len(time_bins[:-1]) / 2) + 1) + + ################################################################################################################### + # activity over time + fig5 = plt.figure(constrained_layout=True, figsize=[15 / inch, 8 / inch]) + gs = gridspec.GridSpec(ncols=2, nrows=1, figure=fig5, hspace=0.05, wspace=0.0, + width_ratios=[6, 3], left=0.1, bottom=0.15, right=0.95, top=0.98) + + ax5 = fig5.add_subplot(gs[0, 0]) + ax6 = fig5.add_subplot(gs[0, 1]) + + ax5.plot(rolled_activ, color=color2[5]) + ax5.plot(rolled_conv_activ, color=color2[4]) + + ax5.plot([16.5*time_factor/5, 6.5*time_factor/bin_len], [9, 9], color=color_diffdays[0], lw=7) + ax5.plot([16.5*time_factor/5, 18.5*time_factor/bin_len], [9, 9], color=color_diffdays[3], lw=7) + ax5.plot([4.5*time_factor/5, 6.5*time_factor/bin_len], [9, 9], color=color_diffdays[3], lw=7) + + # ax5.set_xlim(0, 24) + ax5.set_xticks(np.arange(0, 25, 6) * time_factor / bin_len) + ax5.set_xticklabels(['12:00', '18:00', '00:00', '06:00', '12:00']) + # ax4.set_yticks([0.00, 0.01, 0.02, 0.03]) + ax5.set_xlabel('Time', fontsize=fs) + ax5.set_ylabel('Activity [Changes/ ' + str(bin_len) + ' s]', fontsize=fs) + + for idx_zone, plot_zone, color_zone, day_zone, pos_zone in \ + zip([0, 1, 2, 3], [day, dusk, night, dawn], [6, 1, 4, 0], ['day', 'dusk', 'night', 'dawn'], [1, 2, 3, 4]): + + # boxplot ax1 + props_e = dict(linewidth=2, color=color_dadunida[idx_zone]) + bp = ax6.boxplot(activ_boxes_corr[idx_zone], positions=[pos_zone], widths=0.7, showfliers=False, + boxprops=props_e, medianprops=props_e, capprops=props_e, whiskerprops=props_e) + # ax6.boxplot(activ_boxes_corr, positions=[1, 2, 3, 4], widths=0.7, showfliers=False, + # boxprops=props_e, medianprops=props_e, capprops=props_e, whiskerprops=props_e) + + ax6.set_xticks([1, 2, 3, 4]) + ax6.set_xticklabels(['day', 'dusk', 'night', 'dawn']) + + for ax_idx, axis in enumerate([ax5, ax6]): + axis.text(-0.12, 1, chr(ord('A') + ax_idx), transform=axis.transAxes, fontsize='large') + axis.make_nice_ax() + + # ax5.set_ylim([ax5.get_ylim()[0], 9.5]) + ax6.set_ylim(ax5.get_ylim()) + + dc = [6.3, 7.2, 8.1, 6.3, 9.0, 6.3] + c = 0 + for subset in itertools.combinations([0, 1, 2, 3], 2): + print(subset[0], subset[1]) + # print(stats.mannwhitneyu(activ_boxes_corr[subset[0]], activ_boxes_corr[subset[1]])) + U, p = stats.mannwhitneyu(activ_boxes_corr[subset[0]], activ_boxes_corr[subset[1]]) + print(U, p * 4) + r = 1 - (2 * U) / (len(activ_boxes_corr[subset[0]]) * len(activ_boxes_corr[subset[1]])) + # print(r) + d = cohen_d(activ_boxes_corr[subset[0]], activ_boxes_corr[subset[1]]) + significance_bar(ax6, None, d, subset[0] + 1, subset[1] + 1, dc[c]) + c += 1 + + # fig5.savefig(save_path + 'activ_time.png') + fig5.savefig(save_path + 'activ_time_box.pdf') + plt.show() + + ################################################################################################################### + # statistic + + print(stats.pearsonr(activ_boxes_corr[0], t_b)) + print(stats.pearsonr(activ_boxes_corr[1], time_bins[:-1][dawn])) + print(stats.pearsonr(activ_boxes_corr[2], time_bins[:-1][night])) + print(stats.pearsonr(activ_boxes_corr[3], time_bins[:-1][dusk])) + + embed() \ No newline at end of file diff --git a/plot_aifl.py b/plot_aifl.py new file mode 100644 index 0000000..72d383a --- /dev/null +++ b/plot_aifl.py @@ -0,0 +1,99 @@ +import numpy as np +import matplotlib.pyplot as plt +from IPython import embed +import helper_functions as hf +import do_check_for_overlap as cfo +from params import * +import matplotlib.colors as mcolors +import os + +if __name__ == '__main__': + + ################################################################################################################### + # load data + ################################################################################################################### + # load all the data of one day + filename = sorted(os.listdir('../data/'))[2] + + ident = np.load('../../../data_masterthesis/'+filename+'/all_ident_v.npy', allow_pickle=True) + freq = np.load('../../../data_masterthesis/'+filename+'/all_fund_v.npy', allow_pickle=True) + timeidx = np.load('../../../data_masterthesis/'+filename+'/all_idx_v.npy', allow_pickle=True) + t = np.load('../../../data_masterthesis/'+filename+'/all_times.npy', allow_pickle=True) + + aifl = np.load('../data/'+filename+'/aifl2.npy', allow_pickle=True) + faifl = np.load('../data/'+filename+'/faifl2.npy', allow_pickle=True) + oofl = np.load('../data/'+filename+'/oofl.npy', allow_pickle=True) + faifl = np.delete(faifl, [0], axis=0) + + fish_in_aifl = list(np.unique(np.where(~np.isnan(aifl[:, :, 0]))[1])) + fish_in_faifl = list(np.unique(faifl[:, [1, 3]])) + correct_aifl = sorted(list(set(fish_in_aifl) - set(fish_in_faifl))) + + dt = datetime.datetime.strptime(filename[-5:], '%H_%M') + embed() + quit() + + ################################################################################################################### + # plot traces in oofl + counter = 0 + oofl = np.array(oofl) + for i in range(len(oofl[:, 0])): + channel = int(oofl[i, 0]) + time_diff = timeidx[channel][ident[channel] == oofl[i][1]][-1] - timeidx[channel][ident[channel] == oofl[i][1]][0] + if time_diff >= 0: #4800 + plt.plot(timeidx[channel][ident[channel] == oofl[i][1]], + freq[channel][ident[channel] == oofl[i][1]], Linewidth=2) + plt.text(timeidx[channel][ident[channel] == oofl[i][1]][0] + np.random.rand(1) * 0.3, + freq[channel][ident[channel] == oofl[i][1]][0] + np.random.rand(1) * 0.3, + str(oofl[i][0]) + '_' + str(oofl[i][1]), color='blue') + counter = counter + 1 + plt.show() + + ################################################################################################################### + # plot overlapping traces + new_sorting = cfo.get_list_of_fishN_with_overlap(aifl, fish_in_aifl, timeidx, ident) + for fish_number in new_sorting: + hf.plot_together(timeidx, freq, ident, aifl, int(fish_number), color_vec[0]) + + ################################################################################################################### + # plot fish in faifl + for i in range(len(faifl)): + fishid1 = int(faifl[i, 1]) + fishid2 = int(faifl[i, 3]) + hf.plot_all_channels(timeidx, freq, ident, aifl, fishid1, fishN2=fishid2) + + ################################################################################################################### + # plot all traces + + fig, ax = plt.subplots(1, 1, figsize=(15 / inch, 8 / inch)) + fig.subplots_adjust(left=0.12, bottom=0.15, right=0.98, top=0.98) + for color_counter, fish_number in enumerate(fish_in_aifl): + for channel_idx in [13]: + fish1 = aifl[channel_idx, fish_number, ~np.isnan(aifl[channel_idx, fish_number])] + r1 = len(fish1) + print(fish1) + for len_idx1 in range(r1): + zeit = t[timeidx[channel_idx][ident[channel_idx] == fish1[len_idx1]]] + plot_zeit = [] + for i in range(len(zeit)): + plot_zeit.append(dt + datetime.timedelta(seconds=zeit[i])) + plt.plot(plot_zeit, + freq[channel_idx][ident[channel_idx] == fish1[len_idx1]], + Linewidth=1, label=fish_number, color=color_vec[color_counter+40]) + + ax.set_ylim([450, 1000]) + ax.set_xlabel('Time', fontsize=fs) + ax.set_ylabel('EOD frequency [Hz]', fontsize=fs) + ax.make_nice_ax() + ax.timeaxis() + fig.savefig(save_path_pres+'EOD_sorter.pdf') + fig.savefig('../../thesis/Figures/Methods/EOD_sorter.pdf') + plt.show() + + ################################################################################################################### + # plot + u = np.unique(faifl[:, [1, 3]]) + for fish_number in range(len(u)): + hf.plot_together(timeidx, freq, ident, aifl, int(u[fish_number]), color_vec[fish_number]) + + diff --git a/plot_direction.py b/plot_direction.py new file mode 100644 index 0000000..a050181 --- /dev/null +++ b/plot_direction.py @@ -0,0 +1,149 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +import matplotlib.colors as mcolors +import matplotlib.gridspec as gridspec + +from IPython import embed +from scipy import stats +import os +from IPython import embed +import helper_functions as hf +from params import * +import itertools +from statisitic_functions import * + + +def get_recording_number_in_time_bins(time_bins): + """ + Calculates the number of the recordings in the time bins + + :param time_bins: numpy array with borders of the time bins + :return: time_bins_recording: numpy array with the number of recordings to that specific time bin + """ + # variables + time_bins_recordings = np.zeros(len(time_bins) - 1) + + # load data + for index, filename_idx in enumerate([0, 1, 2, 3]): + filename = sorted(os.listdir('../data/'))[filename_idx] + time_points = np.load('../data/' + filename + '/all_hms.npy', allow_pickle=True) + + # in which bins is this recording, fill time_bins_recordings + unique_time_points = np.unique(np.hstack(time_points)) + for idx, tb in enumerate(time_bins[:-1]): + if np.any((unique_time_points >= tb) & (unique_time_points <= time_bins[idx + 1])): + time_bins_recordings[idx] += 1 + + return time_bins_recordings + + +if __name__ == '__main__': + ################################################################################################################### + # parameter and variables + # plot params + inch = 2.45 + save_path = '../../thesis/Figures/Results/' + kernel_size = 120 + c = 0 + cc = 1 + all_changes = [] + + all_first = [] + all_last = [] + + # time_bins 5 min + time_factor = 60 * 60 + time_bins = np.arange(0, 24 * time_factor + 1, bin_len) + + # changes per bin of 5min for all fish + cbf = np.full([3, len(time_bins) - 1, 300], np.nan) + cbf_counter = 0 + + # kernel + kernel = np.ones(kernel_size) / kernel_size + + ################################################################################################################### + # load data + ################################################################################################################### + # load all the data of one day + cbf = np.load('../data/cbf5.npy', allow_pickle=True) + # embed() + # quit() + N_rec_time_bins = get_recording_number_in_time_bins(time_bins) + + ################################################################################################################### + # direction and act plot + no = np.full(len(time_bins) - 1, np.nan) + down = np.full(len(time_bins) - 1, np.nan) + up = np.full(len(time_bins) - 1, np.nan) + + for i in range(len(up)): + up[i] = np.nansum(cbf[1, i]) # * (60/bin_len) + down[i] = np.nansum(cbf[0, i]) # * (60/bin_len) + no[i] = np.nansum(cbf[2, i]) # * (60/bin_len) + conv_direct = np.convolve(up - down, kernel, 'same') + conv_up = np.convolve(up, kernel, 'same') + conv_down = np.convolve(down, kernel, 'same') + conv_N_rec_tb = np.convolve(N_rec_time_bins, kernel, 'same') + + # rolled time axis for nicer plot midnight in the middle start noon + rolled_time = np.roll(time_bins[:-1] / time_factor, int(len(time_bins[:-1]) / 2)) + rolled_up = np.roll(up / N_rec_time_bins, int(len(time_bins[:-1]) / 2)) + rolled_down = np.roll(down * -1 / N_rec_time_bins, int(len(time_bins[:-1]) / 2)) + rolled_conv_dir = np.roll(conv_direct / conv_N_rec_tb, int(len(time_bins[:-1]) / 2)) + rolled_conv_up = np.roll(conv_up / conv_N_rec_tb, int(len(time_bins[:-1]) / 2)) + rolled_conv_down = np.roll(conv_down * -1 / conv_N_rec_tb, int(len(time_bins[:-1]) / 2)) + ################################################################################################################### + # time zones + time_edges = np.array([4.5, 6.5, 16.5, 18.5]) * time_factor + day = (time_bins[:-1] >= time_edges[1]) & (time_bins[:-1] <= time_edges[2]) + dusk = (time_bins[:-1] >= time_edges[2]) & (time_bins[:-1] <= time_edges[3]) + night = (time_bins[:-1] <= time_edges[0]) | (time_bins[:-1] >= time_edges[3]) + dawn = (time_bins[:-1] >= time_edges[0]) & (time_bins[:-1] <= time_edges[1]) + + ################################################################################################################### + # activity up down over time + fig2, ax2 = plt.subplots(1, 1, figsize=(15 / inch, 8 / inch)) + fig2.subplots_adjust(left=0.1, bottom=0.15, right=0.95, top=0.98) + + ax2.plot(rolled_up, color=color2[2]) + ax2.plot(rolled_down, color=color2[3]) + ax2.plot(rolled_conv_dir, 'k', lw=3) + ax2.plot(rolled_conv_up, color=color2[0], lw=3, label='upstream') + ax2.plot(rolled_conv_down, color=color2[1], lw=3, label='downstream') + + ax2.plot([16.5*time_factor/5, 6.5*time_factor/bin_len], [5.5, 5.5], color=color_diffdays[0], lw=7) + ax2.plot([16.5*time_factor/5, 18.5*time_factor/bin_len], [5.5, 5.5], color=color_diffdays[3], lw=7) + ax2.plot([4.5*time_factor/5, 6.5*time_factor/bin_len], [5.5, 5.5], color=color_diffdays[3], lw=7) + + ax2.set_xlim([0, 24 * time_factor / bin_len]) + ax2.set_xticks(np.arange(0, 25, 3) * time_factor / bin_len) + ax2.set_xticklabels(['12:00', '15:00', '18:00', '21:00', '00:00', '03:00', '06:00', '09:00', '12:00']) + ax2.make_nice_ax() + + # ax2.set_yticks([-7., -5., -2.5, 0, 2.5, 5., 7.]) + ax2.set_xlabel('Time', fontsize=fs) + ax2.set_ylabel('Activity [Changes/ ' + str(bin_len) + ' s]', fontsize=fs) + plt.legend() + + # fig2.savefig(save_path+'activ_updown_time.png') + fig2.savefig(save_path + 'activ_updown_time.pdf') + fig2.savefig(save_path_pres + 'activ_updown_time.pdf') + plt.show() + + ################################################################################################################### + # statistics + + # print('night down: ', np.nanmedian(conv_down[night]/conv_N_rec_tb[night])) + # print('day down: ', np.nanmedian(conv_down[day]/conv_N_rec_tb[day])) + + print('night', stats.pearsonr(conv_down[night]/conv_N_rec_tb[night], conv_up[night]/conv_N_rec_tb[night])) + print('day', stats.pearsonr(conv_down[day]/conv_N_rec_tb[day], conv_up[day]/conv_N_rec_tb[day])) + print('dusk', stats.pearsonr(conv_down[dusk]/conv_N_rec_tb[dusk], conv_up[dusk]/conv_N_rec_tb[dusk])) + print('dawn', stats.pearsonr(conv_down[dawn]/conv_N_rec_tb[dawn], conv_up[dawn]/conv_N_rec_tb[dawn])) + # print('night up: ', np.nanmedian(conv_up[night]/conv_N_rec_tb[night])) + # print('day up: ', np.nanmedian(conv_up[day]/conv_N_rec_tb[day])) + + embed() + quit()