%% Hour by hour analysis script

% Script loads all mice data
% - calculates hour by hour analysis for specific mouse
% - time consuming, not optimised

clear; 
close all;

%% Load data

% clear;
load data\travelled\dlc_travelled_jc20181211_1_05.mat
dlc_travelled_jc20181211_1_05 = DLC_travelled;
travelled_ourDLC_5fps = medfilt1(dlc_travelled_jc20181211_1_05.travelled_ourDLC, 10*5)*5/100;
dlc_travelled_jc20181211_1_05.travelled_ourDLC = travelled_ourDLC_5fps;

load data\travelled\dlc_travelled_jc20181218_1_05.mat
dlc_travelled_jc20181218_1_05 = DLC_travelled;
travelled_ourDLC_5fps = medfilt1(dlc_travelled_jc20181218_1_05.travelled_ourDLC, 10*5)*5/100;
dlc_travelled_jc20181218_1_05.travelled_ourDLC = travelled_ourDLC_5fps;

load data\travelled\dlc_travelled_jc20181219_2_05.mat
dlc_travelled_jc20181219_2_05 = DLC_travelled;
travelled_ourDLC_5fps = medfilt1(dlc_travelled_jc20181219_2_05.travelled_ourDLC, 10*5)*5/100;
dlc_travelled_jc20181219_2_05.travelled_ourDLC = travelled_ourDLC_5fps;

load data\travelled\pix_travelled_jc20181211_1_05.mat
pix_travelled_jc20181211_1_05 = pixels_travelled;
travelled_pixels = medfilt1(pix_travelled_jc20181211_1_05.travelled_pixels, 10*5);
pix_travelled_jc20181211_1_05.travelled_pixels = travelled_pixels;

load data\travelled\pix_travelled_jc20181218_1_05.mat
pix_travelled_jc20181218_1_05 = pixels_travelled;
travelled_pixels = medfilt1(pix_travelled_jc20181218_1_05.travelled_pixels, 10*5);
pix_travelled_jc20181218_1_05.travelled_pixels = travelled_pixels;

load data\travelled\pix_travelled_jc20181219_2_05.mat
pix_travelled_jc20181219_2_05 = pixels_travelled;   
travelled_pixels = medfilt1(pix_travelled_jc20181219_2_05.travelled_pixels, 10*5);
pix_travelled_jc20181219_2_05.travelled_pixels = travelled_pixels;

% Load seizure timestamps
load data\jc20181211_1\jc20181211_1-seizureOnsetsDatenum.mat
jc20181211_1_seizures = seizureOnsetsDatenum;
load data\jc20181218_1\jc20181218_1-seizureOnsetsDatenum.mat
jc20181218_1_seizures = seizureOnsetsDatenum;
load data\jc20181219_2\jc20181219_2-seizureOnsetsDatenum.mat
jc20181219_2_seizures = seizureOnsetsDatenum;

% subtract day of birth 'jc20181211_1'	
date_of_birth = '20181030';
Y = str2double(date_of_birth(1:4));
M = str2double(date_of_birth(5:6));
D = str2double(date_of_birth(7:8));
jc20181211_1_dob = datenum(datetime(Y,M,D));

jc20181211_cut = (dlc_travelled_jc20181211_1_05.t_dlc - jc20181211_1_dob) <= 70;
jc20181211_cut_pix = (pix_travelled_jc20181211_1_05.t_pix - jc20181211_1_dob) <= 70;

dlc_travelled_jc20181211_1_05.travelled_ourDLC = dlc_travelled_jc20181211_1_05.travelled_ourDLC(jc20181211_cut);
dlc_travelled_jc20181211_1_05.t_dlc = dlc_travelled_jc20181211_1_05.t_dlc(jc20181211_cut);

pix_travelled_jc20181211_1_05.travelled_pixels = pix_travelled_jc20181211_1_05.travelled_pixels(jc20181211_cut_pix);
pix_travelled_jc20181211_1_05.t_pix = pix_travelled_jc20181211_1_05.t_pix(jc20181211_cut_pix);
%% hour in datenum

Y = 1990;
M = 1;
D = 1;
H = 0;
MI = 0;
S = 0;
date_start = datetime(Y,M,D,H,MI,S);
H = 1;
MI = 0;
S = 0;
date_end = datetime(Y,M,D,H,MI,S);
hour_length = datenum(date_end) - datenum(date_start);

%%
% windows = [1, 3, 6, 12, 24]*hour_length;
wlen = 24*hour_length;

% pozn: 72 hodin uz je moc dlouhy okno -> vyradi to data z nonseizure

curr_mouse = pix_travelled_jc20181219_2_05;
seizures = jc20181219_2_seizures;

dlc = 0; pix = 1;

if dlc == 1

    len = length(curr_mouse.travelled_ourDLC);
    
    % num_win (x) 24 hours (x) seizure/noseizure/wasseizure (x) approx length
    data = NaN(24, 3, ceil(len / 23));
    
    curr_seizure = 1;
    curr_nonseizure = 1;
    curr_wasseizure = 1;
    
    tic
    for i = 1:len
        
        time = curr_mouse.t_dlc(i);
            
        x = datetime(time, 'ConvertFrom', 'datenum');
        % muzu priradit pomoci x.Hour +1
        
        value = curr_mouse.travelled_ourDLC(i);    
            
        % look if there is seizure occuring
        is_seizure = sum(seizures > time & seizures < time + wlen, 'all');
        was_seizure = sum(seizures < time & seizures > time - wlen, 'all');
        
        if is_seizure % or was seizure + rozdelit na dva pri zapisu
            % write seizure preceding value
            data(x.Hour+1, 2, curr_seizure) = value;
            curr_seizure = curr_seizure + 1;
        else
            % write seizure non preceding value
            data(x.Hour+1, 1, curr_nonseizure) = value;
            curr_nonseizure = curr_nonseizure + 1;
        end
        % klidovy (no seizure) a ty zbyly dva se musi vylucovat
        if was_seizure
            data(x.Hour+1, 3, curr_nonseizure) = value;
            curr_wasseizure = curr_wasseizure + 1;
        end
    end
    toc
end

if pix == 1

    len = length(curr_mouse.travelled_pixels);
    
    % num_win (x) 24 hours (x) seizure/noseizure/wasseizure (x) approx length
    data = NaN(24, 3, ceil(len / 23));
    
    curr_seizure = 1;
    curr_nonseizure = 1;
    curr_wasseizure = 1;
    
    tic
    for i = 1:len
        
        time = curr_mouse.t_pix(i);
            
        x = datetime(time, 'ConvertFrom', 'datenum');
        % muzu priradit pomoci x.Hour +1
        
        value = curr_mouse.travelled_pixels(i);    
            
        % look if there is seizure occuring
        is_seizure = sum(seizures > time & seizures < time + wlen, 'all');
        was_seizure = sum(seizures < time & seizures > time - wlen, 'all');
        
        if is_seizure % or was seizure + rozdelit na dva pri zapisu
            % write seizure preceding value
            data(x.Hour+1, 2, curr_seizure) = value;
            curr_seizure = curr_seizure + 1;
        else
            % write seizure non preceding value
            data(x.Hour+1, 1, curr_nonseizure) = value;
            curr_nonseizure = curr_nonseizure + 1;
        end
        % klidovy (no seizure) a ty zbyly dva se musi vylucovat
        if was_seizure
            data(x.Hour+1, 3, curr_nonseizure) = value;
            curr_wasseizure = curr_wasseizure + 1;
        end
    end
    toc

end

%%

xx = mean(data, 3, 'omitnan');
figure(70);
imagesc(xx);
colorbar
ax = gca;
ax.XTick = [1,2,3];
ax.XTickLabel = {'No Seizure', 'Pre Seizure', 'After Seizure'};
ax.YTick = 1:24;
ax.YTickLabel = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23};
title(['Activity by day hour - ', num2str(wlen/hour_length), ' wlen.']);
fontsize(gcf,scale=2.5);

cmap = colormap;
cmap(1,:) = [0,0,0];

colormap(cmap);

%%

figure(71);
bar(xx);
legend('No seizure', 'Pre Seizure', 'After Seizure','Location', 'best');

seiz_larger = sum(xx(:,2) > xx (:,1));
title(['Seizure larger in ', num2str(seiz_larger), ' / ', num2str(length(xx)), ' cases.']);

%%

nans = sum(isnan(data), 'all');
nonans = sum(~isnan(data), 'all');
disp(['nans: ', num2str(nans), ' x nonans: ', num2str(nonans)])

%%
xx = mean(data_20181211, 3, 'omitnan');
yy = mean(data_20181218, 3, 'omitnan');
zz = mean(data_20181219, 3, 'omitnan');

xx(isnan(xx))=0;
yy(isnan(yy))=0;
zz(isnan(zz))=0;
data_mean = (xx + yy + zz) / 3;

figure(70);
imagesc(data_mean);
colorbar
ax = gca;
ax.XTick = [1,2,3];
ax.XTickLabel = {'No Seizure', 'Pre Seizure', 'After Seizure'};
ax.YTick = 1:24;
ax.YTickLabel = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23};
title(['Activity by day hour - ', num2str(wlen/hour_length), ' wlen.']);
fontsize(gcf,scale=2.5);
