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