jacquelinegoebl/line_tracking_of_fish_movement
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Upload files to ''

Browse Source
This commit is contained in:
jacquelinegoebl 2021-07-16 12:08:26 +02:00
parent b810e4bd77
commit e564736eef
5 changed files with 695 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

134
params.py Normal file
View File

@ -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

132
plot_activity_freq.py Normal file
View File

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

181
plot_activity_time.py Normal file
View File

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

99
plot_aifl.py Normal file
View File

@ -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])

149
plot_direction.py Normal file
View File

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