181 lines
7.4 KiB
Python
181 lines
7.4 KiB
Python
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() |