diff --git a/plot_eigen_trajectories.py b/plot_eigen_trajectories.py new file mode 100644 index 0000000..c770803 --- /dev/null +++ b/plot_eigen_trajectories.py @@ -0,0 +1,90 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +import matplotlib.gridspec as gridspec + +from IPython import embed +import helper_functions as hf +from params import * +import os + +if __name__ == '__main__': + ################################################################################################################### + # parameter and variables + kernel_size = 100 + kernel = np.ones(kernel_size) / kernel_size + + fig1 = plt.figure(constrained_layout=True, figsize=[15 / inch, 15 / inch]) + gs = gridspec.GridSpec(ncols=2, nrows=3, figure=fig1, hspace=0.05, wspace=0.05, + width_ratios=[1, 1], height_ratios=[1, 1, 1], left=0.1, bottom=0.15, right=0.95, top=0.98) + c = 0 + ax_counter = 0 + ################################################################################################################### + # load all the data of one day + + for filename_idx in [0, 1, 2, 3]: + filename = sorted(os.listdir('../data/'))[filename_idx] + + all_max_ch_means = np.load('../data/' + filename + '/all_max_ch.npy', allow_pickle=True) + all_xticks = np.load('../data/' + filename + '/all_xtickses.npy', allow_pickle=True) + all_ipp = np.load('../data/' + filename + '/all_ipp.npy', allow_pickle=True) + power_means = np.load('../data/' + filename + '/power_means.npy', allow_pickle=True) + freq = np.load('../data/' + filename + '/fish_freq_q10.npy', allow_pickle=True) + names = np.load('../data/' + filename + '/fish_species.npy', allow_pickle=True) + + ############################################################################################################### + # get fish + print(filename) + + for fish_number in range(len(power_means)): + if names[fish_number] == 'Eigenmannia' and power_means[fish_number] >= -90.0: + ipp = all_ipp[fish_number] + x_tickses = all_xticks[fish_number] + max_ch_mean = all_max_ch_means[fish_number] + + # smoothing of max channel mean + smooth_mcm = np.convolve(max_ch_mean, kernel, 'valid') + + try: + smooth_x = x_tickses[int(np.ceil(kernel_size / 2)):-int(np.floor(kernel_size / 2))] + except: + embed() + quit() + + ##################################################################################################### + # plot traces + ax1 = fig1.add_subplot(gs[c]) + ax1.imshow(ipp[::20].T[::-1], vmin=-100, vmax=-50, aspect='auto', interpolation='gaussian', + extent=[x_tickses[0], x_tickses[-1], -0.5, 15.5]) + + ax1.plot(smooth_x[::20], smooth_mcm[::20], '.', color=color2[4]) + + ax1.make_nice_ax() + ax1.axhline(7.5, xmin=0, xmax=15, color='white', lw=4) + ax1.set_yticks([0, 1, 2, 3, 4, 5, 6, 7, 7.5, 8, 9, 10, 11, 12, 13, 14, 15]) + ax1.set_yticklabels(['1', '', '3', '', '5', '', '7', '', 'g', '', '10', '', '12', '', '14', '', '16'], + fontsize=9) + + print(ax_counter) + ax1.text(-0.17, 1, chr(ord('A') + ax_counter), transform=ax1.transAxes, fontsize='large') + ax_counter += 1 + ax1.invert_yaxis() + ax1.xaxis_date() + date_format = mdates.DateFormatter('%H:%M') + ax1.xaxis.set_major_formatter(date_format) + + + if c in [3, 4]: + ax1.set_xlabel('Time', fontsize=11) + if c in [1, 2]: + ax1.set_xticks(ax1.get_xticks()[::2]) + + if c in [0, 2, 4]: + ax1.set_ylabel('Electrode', fontsize=11) + + c += 1 + + ################################################################################################################### + # save plot + fig1.savefig(save_path + 'eigen_trajectories.pdf') + plt.show() diff --git a/plot_hist_fish_freq.py b/plot_hist_fish_freq.py new file mode 100644 index 0000000..a234e92 --- /dev/null +++ b/plot_hist_fish_freq.py @@ -0,0 +1,79 @@ +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 +import pandas as pd + +if __name__ == '__main__': + ################################################################################################################### + # parameter + save_path = '../../thesis/Figures/Results/' + inch=2.54 + # color_emf = ['#a2a2a2', '#729ece', '#ed665d'] + + + # lists + eigen = [] + males = [] + females = [] + ################################################################################################################### + # load data + ################################################################################################################### + # load all the data of one day + for index, filename_idx in enumerate([0,1,2,3]): + filename = sorted(os.listdir('../data/'))[filename_idx] + all_q10 = np.load('../data/'+filename+'/fish_freq_q10.npy', allow_pickle=True) + power_means = np.load('../data/' + filename + '/power_means.npy', allow_pickle=True) + names = np.load('../data/' + filename + '/fish_species.npy', allow_pickle=True) + + for fish_number in range(len(power_means)): + if power_means[fish_number] >= -90.0: + + # if all_q10[fish_number,2] <560: + # eigen.append(all_q10[fish_number,2]) + # elif all_q10[fish_number,2] >=760: + # males.append(all_q10[fish_number, 2]) + # elif (all_q10[fish_number,2]>=560) and (all_q10[fish_number,2]<760): + # females.append(all_q10[fish_number, 2]) + if names[fish_number] == 'Eigenmannia': + eigen.append(all_q10[fish_number,2]) + elif names[fish_number] == 'Apteronotus' and all_q10[fish_number,2] >=760: + males.append(all_q10[fish_number, 2]) and all_q10[fish_number,2]>=400 and all_q10[fish_number,2]<760 + elif names[fish_number] == 'Apteronotus': + females.append(all_q10[fish_number, 2]) + + + # eigen.extend(all_q10[:,2][all_q10[:,2]<500]) + # males.extend(all_q10[:,2][all_q10[:,2]>=760]) + # females.extend(all_q10[:,2][(all_q10[:,2]>=500)&(all_q10[:,2]<760)]) + + # plot + fig, ax = plt.subplots(1, 1, figsize=(16 /inch, 8 /inch)) + fig.subplots_adjust(left=0.1, bottom=0.15, right=0.95, top=0.95) + + # ax.hist(eigen, bins=np.arange(400, 1000, 20), color=color_efm[0], label=labels[0]) + # ax.hist(females, bins=np.arange(400, 1000, 20), color=color_efm[1], label=labels[1]) + # ax.hist(males, bins=np.arange(400, 1000, 20), color=color_efm[2], label=labels[2]) + ax.hist([eigen,females,males], bins=np.arange(400, 1000, 20), histtype='bar', stacked=True, label=labels[0:3], color=color_efm[0:3]) + # ax.hist(females, bins=np.arange(400, 1000, 20), color=color_efm[1], label=labels[1]) + # ax.hist(males, bins=np.arange(400, 1000, 20), color=color_efm[2], label=labels[2]) + + ax.make_nice_ax() + ax.set_xlabel('EODf [Hz]', fontsize=fs) + ax.set_ylabel('n', fontsize=fs) + ax.set_yticks([0,5,10,15,20,25]) + ax.set_xlim([400,1000]) + fig.legend(loc='upper right') + # fig.savefig(save_path+'EODf_hist.png') + fig.savefig(save_path+'EODf_hist.pdf') + plt.show() + + embed() diff --git a/plot_jans_pdf.py b/plot_jans_pdf.py new file mode 100644 index 0000000..d7dfe99 --- /dev/null +++ b/plot_jans_pdf.py @@ -0,0 +1,396 @@ +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 mpl_toolkits.axes_grid1.inset_locator import inset_axes +from IPython import embed +from scipy import stats, optimize +import pandas as pd +import math +import os +from IPython import embed + +from eventdetection import threshold_crossings, merge_events +import helper_functions as hf +from params import * +from statisitic_functions import significance_bar, cohen_d +import itertools + + +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 + + +def func(x, a, tau, c): + return a * np.exp(-x / tau) + c + + +def calc_movement(cbf, i): + movement = cbf[0, :, i] + cbf[1, :, i] + movement[np.isnan(movement)] = 0 + re_mov = cbf[0, :, i] - cbf[1, :, i] + re_mov[np.isnan(re_mov)] = 0 + + return movement, re_mov + + +if __name__ == '__main__': + ################################################################################################################### + # parameter and variables + # plot params + inch = 2.54 + + c = 0 + cat_v1 = [0, 0, 750, 0] + cat_v2 = [750, 750, 1000, 1000] + cat_n = ['Eigenmannia', 'Apteronotus', 'Apteronotus'] + + # time + # time_bins 5 min + time_factor = 60 * 60 + # tb2 = np.arange(0, 24 * time_factor + 1, 2) + tb5 = np.arange(0, 24 * time_factor + 1, 5) + # tb10 = np.arange(0, 24 * time_factor + 1, 10) + tb15 = np.arange(0, 24 * time_factor + 1, 15) + # tb30 = np.arange(0, 24 * time_factor + 1, 30) + tb60 = np.arange(0, 24 * time_factor + 1, 60) + tb150 = np.arange(0, 24 * time_factor + 1, 150) + # tb180 = np.arange(0, 24 * time_factor + 1, 180) + tb300 = np.arange(0, 24 * time_factor + 1, 300) + + + # time_edges = np.array([4.5, 6.5, 16.5, 18.5]) * time_factor + # day = time_bins[:-1][(time_bins[:-1] >= time_edges[1]) & (time_bins[:-1] <= time_edges[2])] + # dusk = time_bins[:-1][(time_bins[:-1] >= time_edges[2]) & (time_bins[:-1] <= time_edges[3])] + # night = time_bins[:-1][(time_bins[:-1] <= time_edges[0]) | (time_bins[:-1] >= time_edges[3])] + # dawn = time_bins[:-1][(time_bins[:-1] >= time_edges[0]) & (time_bins[:-1] <= time_edges[1])] + + ################################################################################################################### + # load data + ################################################################################################################### + # load all the data of one day + # cbf2 = np.load('../data/cbf2.npy', allow_pickle=True) + cbf5 = np.load('../data/cbf5.npy', allow_pickle=True) + # cbf10 = np.load('../data/cbf10.npy', allow_pickle=True) + cbf15 = np.load('../data/cbf15.npy', allow_pickle=True) + # cbf30 = np.load('../data/cbf30.npy', allow_pickle=True) + cbf60 = np.load('../data/cbf60.npy', allow_pickle=True) + cbf150 = np.load('../data/cbf150.npy', allow_pickle=True) + # cbf180 = np.load('../data/cbf180.npy', allow_pickle=True) + cbf300 = np.load('../data/cbf300.npy', allow_pickle=True) + + stl = np.load('../data/stl.npy', allow_pickle=True) + names = np.load('../data/n.npy', allow_pickle=True) + freq = np.load('../data/f.npy', allow_pickle=True) + + trajectories = np.load('../data/trajectories.npy', allow_pickle=True) + trajec_x = np.load('../data/trajec_x.npy', allow_pickle=True) + + ############################################################################################################### + # variables + for index, filename_idx in enumerate([0]): + filename = sorted(os.listdir('../data/'))[filename_idx] + all_Ctime_v = np.load('../data/' + filename + '/all_Ctime_v.npy', allow_pickle=True) + sampling_rate = 1 / np.diff(all_Ctime_v[0])[0] # in sec + + cbf_counter = 0 + ################################################################################################################### + # analysis + for i in range(len(trajectories)): + if names[i] == 'unknown': + continue + + + # mov2, re_mov2 = calc_movement(cbf2, cbf_counter) + mov5, re_mov5 = calc_movement(cbf5, cbf_counter) + # mov10, re_mov10 = calc_movement(cbf10, cbf_counter) + mov15, re_mov15 = calc_movement(cbf15, cbf_counter) + # mov30, re_mov30 = calc_movement(cbf30, cbf_counter) + mov60, re_mov60 = calc_movement(cbf60, cbf_counter) + mov150, re_mov150 = calc_movement(cbf150, cbf_counter) + # mov180, re_mov180 = calc_movement(cbf180, cbf_counter) + mov300, re_mov300 = calc_movement(cbf300, cbf_counter) + cbf_counter += 1 + + trajec = trajectories[i] + t_x = trajec_x[i] + + fig = plt.figure(constrained_layout=True, figsize=[20 / inch, 26 / inch]) + gs = gridspec.GridSpec(ncols=2, nrows=6, figure=fig, hspace=0.01, wspace=0.01, + height_ratios=[1, 1, 1, 1, 1, 1], width_ratios=[4,1],left=0.1, bottom=0.15, right=0.95, + top=0.95) + + ax0 = fig.add_subplot(gs[0, 0]) + ax1 = fig.add_subplot(gs[1, 0], sharex=ax0) + ax2 = fig.add_subplot(gs[2, 0], sharex=ax0) + ax3 = fig.add_subplot(gs[3, 0], sharex=ax0) + ax4 = fig.add_subplot(gs[4, 0], sharex=ax0) + ax5 = fig.add_subplot(gs[5, 0], sharex=ax0) + # ax6 = fig.add_subplot(gs[6, 0], sharex=ax0) + ax11 = fig.add_subplot(gs[1, 1]) + ax21 = fig.add_subplot(gs[2, 1]) + ax31 = fig.add_subplot(gs[3, 1]) + ax41 = fig.add_subplot(gs[4, 1]) + ax51 = fig.add_subplot(gs[5, 1]) + # ax61 = fig.add_subplot(gs[6, 1]) + + ax0.plot(t_x/60/60, trajec) + # ax1.plot(tb2[:-1]/60/60, mov2) + ax1.plot(tb5[:-1]/60/60, mov5) + # ax3.plot(tb10[:-1]/60/60, mov10) + ax2.plot(tb15[:-1]/60/60, mov15) + # ax3.plot(tb30[:-1]/60/60, mov30) + ax3.plot(tb60[:-1]/60/60, mov60) + ax4.plot(tb150[:-1]/60/60, mov150) + ax5.plot(tb300[:-1]/60/60, mov300) + + # ax11.hist(mov2, bins=np.linspace(1,np.max(mov2),int(np.max(mov2)))) + ax11.hist(mov5, bins=np.linspace(1,np.max(mov5)+1,int(np.max(mov5)+1))) + # ax31.hist(mov10, bins=np.linspace(1,np.max(mov10),int(np.max(mov10)))) + ax21.hist(mov15, bins=np.linspace(1,np.max(mov15)+1,int(np.max(mov15)+1))) + # ax31.hist(mov30, bins=np.linspace(1,np.max(mov30),int(np.max(mov30)))) + ax31.hist(mov60, bins=np.linspace(1,np.max(mov60)+1,int(np.max(mov60)+1))) + ax41.hist(mov150, bins=np.linspace(1,np.max(mov150)+1,int(np.max(mov150)+1))) + ax51.hist(mov300, bins=np.linspace(1,np.max(mov300)+1,int(np.max(mov300)+1))) + # ax7.hist(mov2) + + tag = ['trajectory', '5', '15', '60', '150', '300'] + for idx, ax in enumerate([ax0, ax1, ax2, ax3, ax4, ax5]): + xl_min=np.min(t_x)/60/60 + xl_max=np.max(t_x)/60/60 + ax.set_xlim([xl_min ,xl_max]) + ax.text(0.01, 0.7, tag[idx], transform=ax.transAxes, fontsize='small') + if ax != ax0: + ax.set_ylabel('n') + + ax0.set_ylim([0,15]) + ax0.invert_yaxis() + ax0.set_ylabel('electrode') + ax5.set_xlabel('Time [h]') + fig.suptitle('EODf '+str(np.round(freq[i],2))+' '+names[i], fontsize=12) + # embed() + # quit() + fig.savefig('../../../jan_plots/trajec'+str(i)+'.pdf') + + plt.close() + + + + + # ############################################################################################################### + # # roll time axis + # start = [] + # stop = [] + # for j in range(len(roaming_events)): + # start.extend(roaming_events[j][0]) + # stop.extend(roaming_events[j][1]) + # + # N_rec_time_bins = get_recording_number_in_time_bins(time_bins[::int((60 / bin_len) * 60)]) + # + # # rolled time axis for nicer plot midnight in the middle start noon + # N_start, bin_edges = np.histogram(np.array(start) * 5, bins=time_bins[::int((60 / bin_len) * 60)]) + # N_stop, bin_edges2 = np.histogram(np.array(stop) * 5, bins=time_bins[::int((60 / bin_len) * 60)]) + # rolled_start = np.roll(N_start / N_rec_time_bins, int(len(N_start) / 2)) + # rolled_stop = np.roll(N_stop / N_rec_time_bins, int(len(N_stop) / 2)) + # rolled_bins = (bin_edges[:-1] / time_factor) + 0.5 + # + # ############################################################################################################### + # # figure 1: max_channel_changes per time zone and per duration of the roaming event + # fig = plt.figure(constrained_layout=True, figsize=[15 / inch, 14 / inch]) + # gs = gridspec.GridSpec(ncols=6, nrows=3, figure=fig, hspace=0.01, wspace=0.01, + # height_ratios=[1, 1, 1], width_ratios=[1, 1, 1, 1, 1, 1], left=0.1, bottom=0.15, right=0.95, + # top=0.95) + # + # ax0 = fig.add_subplot(gs[0, :]) + # ax1 = fig.add_subplot(gs[1, :3]) + # ax2 = fig.add_subplot(gs[1, 3:], sharex=ax1) + # ax3 = fig.add_subplot(gs[2, :2], sharey=ax2) + # ax4 = fig.add_subplot(gs[2, 2:4], sharey=ax2) + # ax5 = fig.add_subplot(gs[2, 4:]) + # + # # axins = inset_axes(ax1, width='30%', height='60%') + # + # # bar plot + # ax0.bar(rolled_bins, rolled_start, color=color2[4]) + # print('bar plot') + # print('day: mean ', np.round(np.mean([rolled_start[:6], rolled_start[18:]]), 2), + # ' std: ', np.round(np.std([rolled_start[:6], rolled_start[18:]]), 2)) + # + # print('night: mean ', np.round(np.mean(rolled_start[6:18]), 2), + # ' std: ', np.round(np.std(rolled_start[6:18]), 2)) + # + # ax0.plot([16.5, 6.5], [20, 20], color=color_diffdays[0], lw=7) + # ax0.plot([16.5, 18.5], [20, 20], color=color_diffdays[3], lw=7) + # ax0.plot([4.5, 6.5], [20, 20], color=color_diffdays[3], lw=7) + # + # ############################################################################################################### + # # curve_fit: tau, std, n + # curvefit_stat = [] + # + # xdata = np.linspace(0.0, 10., 500) + # y_speeds = [] + # for plot_zone, color_zone, day_zone, pos_zone in \ + # zip([day, dusk, night, dawn], [6, 1, 4, 0], ['day', 'dusk', 'night', 'dawn'], [1, 2, 3, 4]): + # + # # boxplot ax1 + # props_e = dict(linewidth=2, color=color2[color_zone]) + # bp = ax1.boxplot(dauer[np.in1d(wann * 5, plot_zone)], positions=[pos_zone], widths=0.7, + # showfliers=False, vert=False, + # boxprops=props_e, medianprops=props_e, capprops=props_e, whiskerprops=props_e) + # + # x_n = [item.get_xdata() for item in bp['whiskers']][1][1] + # n = len(dauer[np.in1d(wann * 5, plot_zone)]) + # ax1.text(x_n + 2, pos_zone, str(n), ha='left', va='center') + # print('dauer: ', day_zone, np.median(dauer[np.in1d(wann * 5, plot_zone)]), + # ' 25, 75: ', np.percentile(dauer[np.in1d(wann * 5, plot_zone)], [25, 75])) + # + # # curve fit + # x_dauer = dauer[dauer <= 10][np.in1d(wann[dauer <= 10] * 5, plot_zone)] + # y_speed = speeds[dauer <= 10][np.in1d(wann[dauer <= 10] * 5, plot_zone)] + # y_speeds.append(y_speed) + # + # popt, pcov = optimize.curve_fit(func, x_dauer, y_speed) + # perr = np.sqrt(np.diag(pcov)) + # print(day_zone, popt, 'perr', perr[1]) + # curvefit_stat.append(np.array([popt[1], perr[1], n])) + # + # # plot dauer vs speed + # ax2.plot(x_dauer, y_speed, 'o', alpha=0.3, color=color2[color_zone]) + # + # ax3.plot(x_dauer, y_speed, 'o', alpha=0.3, color=color2[color_zone]) + # + # # plot curve fit + # ax4.plot(xdata, func(xdata, *popt), '-', color=color2[color_zone], label=day_zone) + # ax4.set_ylim(ax2.get_ylim()) + # + # curvefit_stat = np.array(curvefit_stat) + # # plot std of tau + # ax5.bar([0, 1, 2, 3], curvefit_stat[:, 0], yerr=curvefit_stat[:, 1], color=color2[4]) + # + # ############################################################################################################### + # # statistic + # day_group = [day, dusk, night, dawn] + # for subset in itertools.combinations([0, 1, 2, 3], 2): + # mean1, std1, n1 = curvefit_stat[subset[0]] + # mean2, std2, n2 = curvefit_stat[subset[1]] + # t, p = stats.ttest_ind_from_stats(mean1, std1, n1, mean2, std2, n2) + # d = cohen_d(y_speeds[subset[0]], y_speeds[subset[1]]) + # print(['day', 'dusk', 'night', 'dawn'][subset[0]], ['day', 'dusk', 'night', 'dawn'][subset[1]], 't: ', + # np.round(t, 2), 'p: ', np.round(p, 4), 'd: ', d) + # + # print(stats.mannwhitneyu(dauer[dauer <= 100][np.in1d(wann[dauer <= 100] * 5, day_group[subset[0]])], + # dauer[dauer <= 100][np.in1d(wann[dauer <= 100] * 5, day_group[subset[1]])])) + # if subset[0] == 0 and subset[1] == 2: + # significance_bar(ax5, p, None, subset[0], subset[1], 4.) + # + # ############################################################################################################### + # # labels + # ax0.set_ylabel('# Roaming Events', fontsize=fs) + # ax0.set_xticks([0, 6, 12, 18, 24]) + # ax0.set_xticklabels(['12:00', '18:00', '00:00', '06:00', '12:00']) + # ax0.set_xlabel('Time', fontsize=fs) + # + # ax1.set_yticks([1, 2, 3, 4]) + # ax1.set_yticklabels(['day', 'dusk', 'night', 'dawn']) + # ax1.set_xlabel('Duration [min]', fontsize=fs) + # ax1.invert_yaxis() + # + # ax2.set_xlabel('Duration [min]', fontsize=fs) + # ax2.set_ylabel('Speed [m/min]', fontsize=fs) + # ax2.set_ylim([0, 27]) + # + # ax3.set_ylabel('Speed [m/min]', fontsize=fs) + # ax3.set_xlabel('Duration [min]', fontsize=fs) + # ax3.set_xlim([0, 10]) + # + # ax4.set_xlabel('Duration [min]', fontsize=fs) + # ax4.set_xlim([0, 10]) + # + # ax5.set_xticks([0, 1, 2, 3]) + # ax5.set_xticklabels(['day', 'dusk', 'night', 'dawn'], rotation=45) + # ax5.set_ylabel(r'$\tau$') + # + # tagx = [-0.05, -0.07, -0.07, -0.17, -0.17, -0.17] + # for idx, ax in enumerate([ax0, ax1, ax2, ax3, ax4, ax5]): + # ax.make_nice_ax() + # ax.text(tagx[idx], 1.05, chr(ord('A') + idx), transform=ax.transAxes, fontsize='large') + # + # # fig.align_ylabels() + # # fig.savefig(save_path + 'roaming_events.pdf') + # # fig.savefig(save_path_pres + 'roaming_events.pdf') + # + # ############################################################################################################### + # # figure 2: + # linregress_stat = [] + # fig2 = plt.figure(constrained_layout=True, figsize=[15 / inch, 10 / inch]) + # gs = gridspec.GridSpec(ncols=1, nrows=2, figure=fig2, hspace=0.05, wspace=0.0, + # height_ratios=[1, 2], left=0.1, bottom=0.15, right=0.95, top=0.95) + # + # ax21 = fig2.add_subplot(gs[0, 0]) + # ax23 = fig2.add_subplot(gs[1, 0]) + # + # for plot_zone, color_zone, day_zone, bar_pos, pos_zone in \ + # zip([day, dusk, night, dawn], [6, 1, 4, 0], ['day', 'dusk', 'night', 'dawn'], [-0.3, -0.1, 0.1, 0.3], + # [0, 1, 2, 3]): + # # pdf + # N_roam, bin_roam = np.histogram(roam_dist[np.in1d(wann * 5, plot_zone)], bins=np.linspace(0, 15, 16)) + # N_roam = N_roam / np.sum(N_roam) / (bin_roam[1] - bin_roam[0]) + # ax21.plot(bin_roam[:-1], N_roam, color=color2[color_zone], label=day_zone) + # ax21.set_xlabel('Distance [m]') + # ax21.set_ylabel('PDF') + # ax21.set_xlim([1, 15]) + # + # # duration vs distance + # ax23.plot(dauer[np.in1d(wann * 5, plot_zone)], roam_dist[np.in1d(wann * 5, plot_zone)], 'o', + # color=color2[color_zone], alpha=0.3) + # res = stats.linregress(dauer[np.in1d(wann * 5, plot_zone)], roam_dist[np.in1d(wann * 5, plot_zone)]) + # print(day_zone, res.slope) + # linregress_stat.append(np.array([res.slope, res.stderr, len(dauer[np.in1d(wann * 5, plot_zone)])])) + # ax23.set_xlabel('Duration [min]') + # ax23.set_ylabel('Distance [m]') + # ax23.set_xlim([0, 100]) + # + # print('linregress') + # for subset in itertools.combinations([0, 1, 2, 3], 2): + # mean1, std1, n1 = linregress_stat[subset[0]] + # mean2, std2, n2 = linregress_stat[subset[1]] + # t, p = stats.ttest_ind_from_stats(mean1, std1, n1, mean2, std2, n2) + # d = cohen_d(y_speeds[subset[0]], y_speeds[subset[1]]) + # # print(['day', 'dusk', 'night', 'dawn'][subset[0]], ['day', 'dusk', 'night', 'dawn'][subset[1]], 't: ', + # # np.round(t, 2), 'p: ', np.round(p, 4), 'd: ', d) + # # print(np.round(0.05 / 6, 4)) + # + # for axis in [ax21, ax23]: + # axis.make_nice_ax() + # + # ax21.legend(loc='best', bbox_to_anchor=(0.5, 0.7, 0.5, 0.5), ncol=2) + # + # fig2.savefig(save_path_pres + 'roaming_distance.pdf') + # fig2.savefig(save_path + 'roaming_distance.pdf') + # + # plt.show() + # + # # df = pd.DataFrame({'duration': dauer, 'speed': speeds, 'distance': roam_dist}) + # embed() + # quit() diff --git a/plot_location.py b/plot_location.py new file mode 100644 index 0000000..a983d0d --- /dev/null +++ b/plot_location.py @@ -0,0 +1,152 @@ +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 * + + +if __name__ == '__main__': + + ################################################################################################################### + # parameter and variables + # plot params + inch = 2.45 + save_path = '../../thesis/Figures/Results/' + + # lists + afx = [] # all_flat_x + + ################################################################################################################### + # load data + for filename_idx in [0, 1, 2, 3]: + filename = sorted(os.listdir('../data/'))[filename_idx] + all_xticks = np.load('../data/' + filename + '/all_xtickses.npy', allow_pickle=True) + + # lists + afx.extend(np.unique(np.hstack(all_xticks))) + afx = sorted(np.unique(afx)) + c = 0 + ################################################################################################################### + # make alm: all location matrix + alm = np.full([len(afx), 16], 0) + for filename_idx in [0, 1, 2, 3]: + filename = sorted(os.listdir('../data/'))[filename_idx] + all_max_ch_means = np.load('../data/' + filename + '/all_max_ch.npy', allow_pickle=True) + all_xticks = np.load('../data/' + filename + '/all_xtickses.npy', allow_pickle=True) + power_means = np.load('../data/' + filename + '/power_means.npy', allow_pickle=True) + names = np.load('../data/' + filename + '/fish_species.npy', allow_pickle=True) + + for fish_number in range(len(power_means)): + if power_means[fish_number] >= -90 and names[fish_number] != 'unknown': + c += 1 + x_tickses = all_xticks[fish_number] + max_ch_mean = all_max_ch_means[fish_number] + + alm[np.flatnonzero(np.isin(afx, x_tickses)), max_ch_mean] += 1 + + hist = np.sum(alm, axis=0)/len(alm) + # hist = np.mean(alm, axis=0) #/len(alm) + # hist_err = np.std(alm, axis=0) #/len(alm) + + ################################################################################################################### + # figure + # fig1, [ax1, ax11] = plt.subplots(1, 2, figsize=(16 / inch, 10 / inch)) + # fig1.subplots_adjust(left=0.12, bottom=0.15, right=0.99, top=0.99, wspace=0.0) + + fig1 = plt.figure(figsize=[13 / inch, 10 / inch]) + spec = gridspec.GridSpec(ncols=2, nrows=1, figure=fig1, hspace=0.5, wspace=0.05, + width_ratios=[7,1], left=0.12, bottom=0.15, right=0.99, top=0.95) + ax1 = fig1.add_subplot(spec[0, 0]) + ax11 = fig1.add_subplot(spec[0, 1]) + ################################################################################################################### + # plot + ax1.imshow(alm[::20].T[::-1], vmin=0.0, vmax=7.0, aspect='auto', interpolation='gaussian', + extent=[afx[0], afx[-1], -0.5, 15.5]) + ax1.beautimechannelaxis() + ax1.timeaxis() + # fig1.autofmt_xdate() + + ax11.barh(np.arange(0,16), hist, color=color2[4]) + # ax11.barh(np.arange(0,16), hist, color=color2[4], xerr=hist_err) + ax11.make_nice_ax() + + ax11.set_ylim(-0.5,15.5) + ax11.set_yticks([]) + ax11.set_xticks([]) + ax11.invert_yaxis() + + fig1.savefig(save_path+'all_trajectories.pdf') + + x = np.array(afx) + x1 = np.where((x>=datetime_box[0])&(x<=datetime_box[2]))[0] + x2 = np.where((x>=datetime_box[4])&(x<=datetime_box[5]))[0] + y = np.concatenate((alm[x1],alm[x2])) + y_a=np.sum(y, axis=0)/len(y) + + x3 = np.where((x>=datetime_box[2])&(x<=datetime_box[4]))[0] + y_b = np.sum(alm[x3], axis=0)/len(alm[x3]) + + fig2 = plt.figure(figsize=[13 / inch, 10 / inch]) + spec = gridspec.GridSpec(ncols=2, nrows=1, figure=fig2, hspace=0.5, wspace=0.05, + width_ratios=[1, 1], left=0.12, bottom=0.15, right=0.99, top=0.95) + ax2 = fig2.add_subplot(spec[0, 0]) + ax22 = fig2.add_subplot(spec[0, 1]) + ################################################################################################################### + # plot + ax2.barh(np.arange(1, 17), y_a, color=color2[4]) + ax22.barh(np.arange(1, 17), y_b, color=color2[4]) + + for ax in [ax2, ax22]: + ax.make_nice_ax() + ax.set_ylim(-0.5, 15.5) + ax.set_ylim(0, 6.0) + # ax.set_xticks([]) + 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', '', 'g', '', '10', '', '12', '', '14', '', '16'], fontsize=9) + ax.invert_yaxis() + ax22.set_yticklabels([]) + + plt.show() + + embed() + quit() + + ################################################################################################################### + # figure 2: each day on its own + # for filename_idx in [1, 4, 6]: + # filename = sorted(os.listdir('../../../data/'))[filename_idx] + # all_max_ch_means = np.load('../data/' + filename + '/all_max_ch.npy', allow_pickle=True) + # all_xticks = np.load('../data/' + filename + '/all_xtickses.npy', allow_pickle=True) + # power_means = np.load('../data/' + filename + '/power_means.npy', allow_pickle=True) + # + # # lists + # flat_x = np.unique(np.hstack(all_xticks)) + # afx.extend(flat_x) + # location_matrix = np.full([len(flat_x), 16], 0) + # + # for fish_number in range(len(power_means)): + # # if power_means[fish_number] >= -65: + # x_tickses = all_xticks[fish_number] + # max_ch_mean = all_max_ch_means[fish_number] + # + # location_matrix[np.flatnonzero(np.isin(flat_x, x_tickses)), max_ch_mean] += 1 + # + # fig2 = plt.figure(figsize=[16 / inch, 12 / inch]) + # spec = gridspec.GridSpec(ncols=1, nrows=1, figure=fig2, hspace=0.5, wspace=0.5) + # ax2 = fig2.add_subplot(spec[0, 0]) + # + # ax2.imshow(location_matrix[::20].T[::-1], aspect='auto', interpolation='gaussian', cmap='jet', + # extent=[flat_x[0], flat_x[-1], 0, 15]) + # ax2.beautimechannelaxis() + # ax2.timeaxis() + # fig2.autofmt_xdate() + # + # plt.show() + # embed() \ No newline at end of file diff --git a/plot_material_connect_ID.py b/plot_material_connect_ID.py new file mode 100644 index 0000000..3cd1b7f --- /dev/null +++ b/plot_material_connect_ID.py @@ -0,0 +1,136 @@ +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +import matplotlib.colors as mcolors +import matplotlib.gridspec as gridspec +import math +from IPython import embed +from scipy import stats +import os +from IPython import embed +from params import * +import datetime +import itertools +import pandas as pd +import helper_functions as hf + +inch = 2.45 +save_path = '../../thesis/Figures/Methods/' + +# color2 = ['#f39c12', '#d35400', '#f8c471', '#dc7633', '#c0392b', '#f1c40f'] + + +# fig, ax = plt.subplots(1, 1, figsize=(16 / inch, 5 / inch)) +# fig.subplots_adjust(left=0.1, bottom=0.30, right=0.95, top=0.95) +# +# ax.plot([0,0.5,1,2], [1,1,1,1], '*', color=color2[0], label='Fish 1', markersize=7) +# ax.plot([0,1,1.5,2], [2,2,2,2], '*', color=color2[1], label='Fish 2', markersize=7) +# +# ax.set_ylim(0.5,2.5) +# ax.set_xlim(-0.2,2.7) +# ax.set_xticks([0,1,2]) +# ax.set_yticks([1,2]) +# ax.set_yticklabels([700,720]) +# ax.set_xlabel('Time', fontsize=fs) +# ax.set_ylabel('Frequency', fontsize=fs) +# ax.make_nice_ax() +# ax.legend(loc='upper right') +# +# # fig.savefig(save_path+'till_data_structure.png') +# fig.savefig(save_path + 'till_data_structure.pdf') +# +# plt.show() + +x1 = [2.5, 5] +y1 = [1, 1] + +x2 = [2.5, 4] +y2 = [1.2, 1.2] + +x3 = [1,3] +y3 = [2.5,2.5] + +x4 = [1,1.5] +y4 = [1.7,1.7] + +# +fig1, [ax1,ax2,ax3] = plt.subplots(3, 1, figsize=(16 / inch, 10 / inch)) +fig1.subplots_adjust(left=0.1, bottom=0.1, right=0.85, top=0.95) + +ax2.plot(x1,y1, '-', color=color2[0], label='Fish 1', linewidth=2, markersize=7) +ax2.plot(x2,y2, '-', color=color2[1], label='ID 2', linewidth=2, markersize=7) +ax2.plot(x3,y3, '-', color=color2[4], label='ID 3', linewidth=2, markersize=7) +ax2.plot(x4,y4, '-', color='gray', label='ID 4', linewidth=2, markersize=7) + +ax1.plot(x1,y1, '-', color=color2[0], linewidth=2, markersize=7) +ax1.plot(x2,y2, '-', color=color2[1], linewidth=2, markersize=7) +ax1.plot(x3,y3, '-', color=color2[4], linewidth=2, markersize=7) +ax1.plot(x4,y4, '-', color='gray', linewidth=2, markersize=7) + +ax3.plot(x1,y1, '-', color=color2[0], linewidth=2, markersize=7) +ax3.plot(x2,y2, '-', color=color2[1], linewidth=2, markersize=7) +ax3.plot(x3,y3, '-', color='gray', linewidth=2, markersize=7) +ax3.plot([2.5, 4], y1, 'o', color=color2[0], linewidth=2, markersize=7) +ax3.plot(x2,y2, 'o', color=color2[1], linewidth=2, markersize=7) +ax3.plot([2.5, 3], y3, 'o', color='gray', linewidth=2, markersize=7) + +ax1.vlines(2.2, ymin=0, ymax=5.5, color='k', linestyle='dashed') +ax1.vlines(5.3, ymin=0, ymax=5.5, color='k', linestyle='dashed') +ax1.plot([2.2,2.5], [0.8,0.8],'k', linestyle='dotted') +ax1.plot([5.0,5.3], [0.8,0.8],'k', linestyle='dotted') + +ax2.axvspan(2.2, 2.8, color=color2[0], alpha=0.1) +ax2.axvspan(2.2, 2.8, color=color2[1], alpha=0.1) +ax2.axvspan(2.2, 2.8, color=color2[4], alpha=0.1) +ax2.axvspan(3.7, 4.3, color=color2[1], alpha=0.1) +ax2.axvspan(2.7, 3.3, color=color2[4], alpha=0.1) + +ax3.hlines(0.6, xmin=0, xmax=5.5, color='k', linestyle='dashed') +ax3.hlines(1.4, xmin=0, xmax=5.5, color='k', linestyle='dashed') +# ax3.vlines(3.7, ymin=0.8, ymax=1.2, color='k') +# ax3.vlines(4.3, ymin=0.8, ymax=1.2, color='k') +# ax3.hlines(0.8, ymin=2.3, ymax=2.7, color='k') +# ax3.hlines(3.3, ymin=2.3, ymax=2.7, color='k') +# ax3.plot([2.2,2.8], [1.2,1.2],'k') +# ax3.plot([3.7,4.3], [1.2,1.2],'k') +# ax3.plot([2.7,3.3], [2.7,2.7],'k') +# ax3.plot([2.2,2.8], [0.8,0.8],'k') +# ax3.plot([3.7,4.3], [0.8,0.8],'k') +# ax3.plot([2.7,3.3], [2.3,2.3],'k') + +for idx, ax in enumerate([ax1,ax2,ax3]): + ax.text(-0.09, 1, chr(ord('A') + idx), transform=ax.transAxes, fontsize='large') + ax.set_ylim(0,3.5) + ax.set_xlim(0.5,5.5) + ax.set_yticks([]) + ax.set_xticks([]) + # ax.set_xticklabels([]) + # ax.set_yticklabels([]) + ax.set_xlabel('Time', fontsize=fs) + ax.set_ylabel('EODf', fontsize=fs) + ax.make_nice_ax() + +ax = plt.gca() +ax2.legend(loc='center right', bbox_to_anchor=(1.2, 0.5)) + +# fig1.savefig(save_path + 'connect_ID.png') +fig1.savefig(save_path + 'connect_ID.pdf') + +plt.show() + +exit() + + +# temperatur +t = [] + +for index, filename_idx in enumerate([1, 4, 6]): + filename = sorted(os.listdir('../../../data/'))[filename_idx] + + temp = pd.read_csv('../../../data/' + filename + '/temperatures.csv', sep=';') + t.append(np.array(temp.values.tolist())[:,1]) + + +print(np.mean(np.hstack(t))) +print(np.min(np.hstack(t))) +print(np.max(np.hstack(t)))