% MR 21.05.2024

clc; clear; close all;

%% Data loading

load lkdata

% lkdata"X""state": matrix of double  [number of subjects x scale k] containg
% the natural logaritmic values of the fractal curve L(k)

% if "X" == 1 : subjects from 0 to 1 y.o.
% if "X" == 2: subjects from 1 to 17 y.o.
% if "state" == sleep : computed from EEG during sleep
% if "state" == wake : computed from EEG during wakefulness

%% Wakefulness

x_axis = log(1:40);

lkdata1wake_avg = mean(lkdata1wake); % average L curve in subjects aged 0-1 y.o.
lkdata2wake_avg = zeros(40,4); % to contain average L curve in subjects aged [2-4], [5-8], [9-12], [13-17] y.o.

index = 1;
for k = 1:40:121
    lkdata2wake_avg(:,index) = mean(lkdata2wake(k:k+39,:));
    index = index+1;
end

% computation of the first derivative of the L curves during wakefulness
d1_lkdata1wake_avg = gradient(lkdata1wake_avg',x_axis); % 0-1 y.o.
d1_lkdata2wake_avg = gradient(lkdata2wake_avg',x_axis,1); % [2-4], [5-8], [9-12], [13-17] y.o.

% computation of the second derivative of the L curves during wakefulness
d2_lkdata1wake_avg = gradient(d1_lkdata1wake_avg',x_axis); % 0-1 y.o.
d2_lkdata2wake_avg = gradient(d1_lkdata2wake_avg,x_axis,1); % [2-4], [5-8], [9-12], [13-17] y.o.

%% Sleep

lkdata1sleep_avg = mean(lkdata1sleep);
lkdata2sleep_avg = [mean(lkdata2sleep(1:20,:))' mean(lkdata2sleep(21:36,:))' mean(lkdata2sleep(37:42,:))' mean(lkdata2sleep(43:45,:))'];

d1_lkdata1sleep_avg = gradient(lkdata1sleep_avg',x_axis);
d1_lkdata2sleep_avg = gradient(lkdata2sleep_avg',x_axis,1);

d2_lkdata1sleep_avg = gradient(d1_lkdata1sleep_avg',x_axis);
d2_lkdata2sleep_avg = gradient(d1_lkdata2sleep_avg,x_axis,1);

%% Choice of k lin from average L curve

% computation of the tuning parameter klin - wakefulness

d1_all_wake = d1_lkdata1wake_avg; % average L curve for all subjects (0-1 y) during wakefulness
avg_mob_wake = cumsum(d1_all_wake)./(1:40); % cumulative average that is an appoximation of the HFD
k_lin_wake = zeros(size(d1_lkdata2wake_avg,1)+1,1);
k_lin_wake (1,1) = find(floor(100* abs((d1_all_wake./avg_mob_wake)-1) ) > 5,1,'first'); % first value of the first derivative k that overcomes the cumulative average of 5%

for w = 1:size(d1_lkdata2wake_avg,1)
    d1_all_wake = d1_lkdata2wake_avg(w,:); % average L curve for all subjects (>1y) during wakefulness
    avg_mob_wake = cumsum(d1_all_wake)./(1:40); % cumulative average that is an appoximation of the HFD
    k_lin_wake (w+1,1) = find(floor(100* abs((d1_all_wake./avg_mob_wake)-1) ) > 5,1,'first'); % first value of the first derivative k that overcomes the cumulative average of 5%
end

k_lin_wake_final = floor(mean(k_lin_wake));

% computation of the tuning parameter klin - sleep

d1_all_sleep = d1_lkdata1sleep_avg; 
avg_mob_sleep = cumsum(d1_all_sleep)./(1:40); 
k_lin_sleep = zeros(size(d1_lkdata2sleep_avg,1)+1,1);
k_lin_sleep (1,1) = find(floor(100* abs((d1_all_sleep./avg_mob_sleep)-1) ) > 5,1,'first'); 

for w = 1:size(d1_lkdata2sleep_avg,1)
    d1_all_sleep = d1_lkdata2sleep_avg(w,:); 
    avg_mob_sleep = cumsum(d1_all_sleep)./(1:40); 
    k_lin_sleep (w+1,1) = find(floor(100* abs((d1_all_sleep./avg_mob_sleep)-1) ) > 5,1,'first'); 
end

k_lin_sleep_final = floor(mean(k_lin_sleep));

%% Plot
figure()
% plot of L curve for each age group during sleep 
subplot(321)
plot(x_axis,lkdata1sleep_avg,'.-');
hold on
plot(x_axis,lkdata2sleep_avg,'.-');
xline(x_axis(k_lin_sleep_final), '--', 'Color', [0 0 0 0.2], 'LineWidth', 0.8); % k=6 line
ylabel('log(L(k))');
xlabel('log(k)');
title('Sleep');
axis([0 x_axis(end) 0 7])
grid on
legend('0-1 y','2-4 y','5-8 y','9-12 y','13-17 y')

subplot(323)
plot(x_axis, d1_lkdata1sleep_avg, '.-');
hold on
plot(x_axis, d1_lkdata2sleep_avg', '.-');
xline(x_axis(k_lin_sleep_final), '--', 'Color', [0 0 0 0.2], 'LineWidth', 0.8);
xlabel('log(k)');
ylabel('\delta(log(L(k)))/\delta(log(k))');
title('1^{st} derivative');
axis([0 x_axis(end) -3.5 3.5])
grid on

subplot(325)
plot(x_axis, d2_lkdata1sleep_avg, '.-');
hold on
plot(x_axis, d2_lkdata2sleep_avg', '.-');
xline(x_axis(k_lin_sleep_final), '--', 'Color', [0 0 0 0.2], 'LineWidth', 0.8);
xlabel('log(k))');
ylabel('\delta^2(log(L(k)))/\delta^2(log(k))');
title('2^{nd} derivative');
grid on
ylim([-6 8])
xlim([0 x_axis(end)])

% plot of L curve for each age group during wakefulness 
subplot(322)
plot(x_axis,lkdata1wake_avg,'.-');
hold on
plot(x_axis,lkdata2wake_avg,'.-');
xline(x_axis(k_lin_wake_final), '--', 'Color', [0 0 0 0.2], 'LineWidth', 0.8); % k=6 line
ylabel('log(L(k))');
xlabel('log(k)');
title('Wakefulness');
axis([0 x_axis(end) 0 7])
legend('0-1 y','2-4 y','5-8 y','9-12 y','13-17 y')
grid on

subplot(324)
% plot of the first derivative during wakefulness
plot(x_axis, d1_lkdata1wake_avg, '.-');
hold on
plot(x_axis, d1_lkdata2wake_avg', '.-');
xline(x_axis(k_lin_wake_final), '--', 'Color', [0 0 0 0.2], 'LineWidth', 0.8);
xlabel('log(k)');
ylabel('\delta(log(L(k)))/\delta(log(k))');
title('1^{st} derivative');
axis([0 x_axis(end) -3.5 3.5])
grid on

subplot(326)
% plot of the second derivative during wakefulness
plot(x_axis, d2_lkdata1wake_avg, '.-');
hold on
plot(x_axis, d2_lkdata2wake_avg', '.-');
xline(x_axis(k_lin_wake_final), '--', 'Color', [0 0 0 0.2], 'LineWidth', 0.8);
xlabel('log(k))');
ylabel('\delta^2(log(L(k)))/\delta^2(log(k))');
title('2^{nd} derivative');
grid on
ylim([-6 8])
xlim([0 x_axis(end)])


set(gcf, 'Units', 'Inches', 'Position', [0, 0, 8, 7]);

% Assicurati che la modalità PaperPosition rispetti le dimensioni su schermo
set(gcf, 'PaperPositionMode', 'auto');

% Esporta TUTTA la figura in TIFF 600 dpi
exportgraphics(gcf, 'mioPlot.tiff', 'Resolution', 600);