%% Create average graph for each muscle for each case

clear all;
close all;

% Load files created by "GetStepsFromEMG.mlx". If you have a different file structure you have to adjust the paths.

kind = 'Gastro';
speed = "PWS"; %SlowWalking

% Slow Walking
rl_sub1 = load("DataAnalysis\sub1\"+speed+"\"+kind + "\aiControl.mat","corrected_graphs_final_ai").corrected_graphs_final_ai;
t_new_rl_sub1=(0:0.0005:0.0005*(length(rl_sub1)-1));
emp_sub1 = load("DataAnalysis\sub1\"+speed+"\"+kind + "\empDev.mat","corrected_graphs_final_self").corrected_graphs_final_self;
t_new_emp_sub1=(0:0.0005:0.0005*(length(emp_sub1)-1));
nc_sub1 = load("DataAnalysis\sub1\"+speed+"\"+kind + "\withoutControl.mat","corrected_graphs_final_withoutControl").corrected_graphs_final_withoutControl;
t_new_nc_sub1=(0:0.0005:0.0005*(length(nc_sub1)-1));
ne_sub1 = load("DataAnalysis\sub1\"+speed+"\"+kind + "\withoutExo.mat","corrected_graphs_final_ref").corrected_graphs_final_ref;
t_new_ne_sub1=(0:0.0005:0.0005*(length(ne_sub1)-1));

rl_sub2 = load("DataAnalysis\sub2\"+speed+"\"+kind + "\rlControl.mat","corrected_graphs_final_ai").corrected_graphs_final_ai;
t_new_rl_sub2=(0:0.0005:0.0005*(length(rl_sub2)-1));
emp_sub2 = load("DataAnalysis\sub2\"+speed+"\"+kind + "\empDev.mat","corrected_graphs_final_self").corrected_graphs_final_self;
t_new_emp_sub2=(0:0.0005:0.0005*(length(emp_sub2)-1));
nc_sub2 = load("DataAnalysis\sub2\"+speed+"\"+kind + "\withoutControl.mat","corrected_graphs_final_withoutControl").corrected_graphs_final_withoutControl;
t_new_nc_sub2=(0:0.0005:0.0005*(length(nc_sub2)-1));
ne_sub2 = load("DataAnalysis\sub2\"+speed+"\"+kind + "\withoutExo.mat","corrected_graphs_final_ref").corrected_graphs_final_ref;
t_new_ne_sub2=(0:0.0005:0.0005*(length(ne_sub2)-1));

% NOT USED IN PWS
% rl_sub3 = load("DataAnalysis\sub3\"+speed+"\"+kind + "\rlControl.mat","corrected_graphs_final_ai").corrected_graphs_final_ai;
% t_new_rl_sub3=(0:0.0005:0.0005*(length(rl_sub3)-1));
% emp_sub3 = load("DataAnalysis\sub3\"+speed+"\"+kind + "\empDev.mat","corrected_graphs_final_self").corrected_graphs_final_self;
% t_new_emp_sub3=(0:0.0005:0.0005*(length(emp_sub3)-1));
% nc_sub3 = load("DataAnalysis\sub3\"+speed+"\"+kind + "\withoutControl.mat","corrected_graphs_final_withoutControl").corrected_graphs_final_withoutControl;
% t_new_nc_sub3=(0:0.0005:0.0005*(length(nc_sub3)-1));
% ne_sub3 = load("DataAnalysis\sub3\"+speed+"\"+kind + "\withoutExo.mat","corrected_graphs_final_ref").corrected_graphs_final_ref;
% t_new_ne_sub3=(0:0.0005:0.0005*(length(ne_sub3)-1));

rl_sub4 = load("DataAnalysis\sub4\"+speed+"\"+kind + "\rlControl.mat","corrected_graphs_final_ai").corrected_graphs_final_ai;
t_new_rl_sub4=(0:0.0005:0.0005*(length(rl_sub4)-1));
emp_sub4 = load("DataAnalysis\sub4\"+speed+"\"+kind + "\empDev.mat","corrected_graphs_final_self").corrected_graphs_final_self;
t_new_emp_sub4=(0:0.0005:0.0005*(length(emp_sub4)-1));
nc_sub4 = load("DataAnalysis\sub4\"+speed+"\"+kind + "\withoutControl.mat","corrected_graphs_final_withoutControl").corrected_graphs_final_withoutControl;
t_new_nc_sub4=(0:0.0005:0.0005*(length(nc_sub4)-1));
ne_sub4 = load("DataAnalysis\sub4\"+speed+"\"+kind + "\withoutExo.mat","corrected_graphs_final_ref").corrected_graphs_final_ref;
t_new_ne_sub4=(0:0.0005:0.0005*(length(ne_sub4)-1));

% Create the average step of all subjects by using the averages of each subject (excluding some outliers). When adjusting which data of the different subjects should be used or not, make sure that you also change the necessary parts for "time_corrected_subx" and "steps_corrected_xx", as is demonstrated below.

%RL

% not considered for RF in PWS and Slow Walking     :   sub2

mean_duration = round(mean([t_new_rl_sub1(end), t_new_rl_sub2(end),t_new_rl_sub3(end),t_new_rl_sub4(end)]),3);
factor = mean_duration./[t_new_rl_sub1(end), t_new_rl_sub2(end),t_new_rl_sub3(end),t_new_rl_sub4(end)];
t_new_rl=(0:0.0005:mean_duration);

time_corrected_sub1 = t_new_rl_sub1.*factor(1);
time_corrected_sub2 = t_new_rl_sub2.*factor(2);
time_corrected_sub3 = t_new_rl_sub3.*factor(3);
time_corrected_sub4 = t_new_rl_sub4.*factor(4);

steps_corrected_rl = [];
steps_corrected_rl(1,:) = interp1(time_corrected_sub1,rl_sub1,t_new_rl, 'linear','extrap');
steps_corrected_rl(2,:) = interp1(time_corrected_sub2,rl_sub2,t_new_rl, 'linear','extrap');
steps_corrected_rl(3,:) = interp1(time_corrected_sub3, rl_sub3, t_new_rl, 'linear','extrap');
steps_corrected_rl(4,:) = interp1(time_corrected_sub4, rl_sub4, t_new_rl, 'linear','extrap');

average_rl = mean(steps_corrected_rl);

%Emp

% not considered for RF and VAS in PWS  :   sub 2
% not considered for BF in PWS          :   sub 4

% not considered for GAS, RF, and VAS in Slow Walking    :   sub2

mean_duration = round(mean([t_new_emp_sub1(end),t_new_emp_sub2(end),t_new_emp_sub3(end),t_new_emp_sub4(end)]),4);
factor = mean_duration./[t_new_emp_sub1(end),t_new_emp_sub2(end),t_new_emp_sub3(end),t_new_emp_sub4(end)];
t_new_emp=(0:0.0005:mean_duration);

time_corrected_sub1 = t_new_emp_sub1.*factor(1);
time_corrected_sub2 = t_new_emp_sub2.*factor(2);
time_corrected_sub3 = t_new_emp_sub3.*factor(3);
time_corrected_sub4 = t_new_emp_sub4.*factor(4);

steps_corrected_emp = [];
steps_corrected_emp(1,:) = interp1(time_corrected_sub1,emp_sub1,t_new_emp, 'linear','extrap');
steps_corrected_emp(2,:) = interp1(time_corrected_sub2,emp_sub2,t_new_emp, 'linear','extrap');
steps_corrected_emp(3,:) = interp1(time_corrected_sub3, emp_sub3, t_new_emp, 'linear','extrap');
steps_corrected_emp(4,:) = interp1(time_corrected_sub4, emp_sub4, t_new_emp, 'linear','extrap');

average_emp = mean(steps_corrected_emp);

% NC

% not considered for GAS, VAS and RF in Slow Walking: sub2
% not considered for RF in PWS      :   sub2
% not considered for BF in PWS      :   sub4

mean_duration = round(mean([t_new_nc_sub1(end),t_new_nc_sub3(end),t_new_nc_sub4(end)]),3);
factor = mean_duration./[t_new_nc_sub1(end),t_new_nc_sub3(end),t_new_nc_sub4(end)];
t_new_nc=(0:0.0005:mean_duration);

time_corrected_sub1 = t_new_nc_sub1.*factor(1);
% time_corrected_sub2 = t_new_nc_sub2.*factor(2);
time_corrected_sub3 = t_new_nc_sub3.*factor(2);
time_corrected_sub4 = t_new_nc_sub4.*factor(3);

steps_corrected_nc = [];
steps_corrected_nc(1,:) = interp1(time_corrected_sub1,nc_sub1,t_new_nc, 'linear','extrap');
% steps_corrected_nc(2,:) = interp1(time_corrected_sub2,nc_sub2,t_new_nc, 'linear','extrap');
steps_corrected_nc(2,:) = interp1(time_corrected_sub3, nc_sub3, t_new_nc, 'linear','extrap');
steps_corrected_nc(3,:) = interp1(time_corrected_sub4, nc_sub4, t_new_nc, 'linear','extrap');

average_nc = mean(steps_corrected_nc);

% NE

% not considered for RF in PWS     :    sub1

mean_duration = round(mean([t_new_ne_sub1(end), t_new_ne_sub2(end),t_new_ne_sub3(end), t_new_ne_sub4(end)]),4);
factor = mean_duration./[t_new_ne_sub1(end),t_new_ne_sub2(end), t_new_ne_sub3(end),t_new_ne_sub4(end)];
t_new_ne=(0:0.0005:mean_duration);

time_corrected_sub1 = t_new_ne_sub1.*factor(1);
time_corrected_sub2 = t_new_ne_sub2.*factor(2);
time_corrected_sub3 = t_new_ne_sub3.*factor(3);
time_corrected_sub4 = t_new_ne_sub4.*factor(4);

steps_corrected_ne = [];
steps_corrected_ne(1,:) = interp1(time_corrected_sub1,ne_sub1,t_new_ne, 'linear','extrap');
steps_corrected_ne(2,:) = interp1(time_corrected_sub2,ne_sub2,t_new_ne, 'linear','extrap');
steps_corrected_ne(3,:) = interp1(time_corrected_sub3, ne_sub3, t_new_ne, 'linear','extrap');
steps_corrected_ne(4,:) = interp1(time_corrected_sub4, ne_sub4, t_new_ne, 'linear','extrap');

average_ne = mean(steps_corrected_ne);

%remove bias
average_ne = average_ne - min(average_ne);
average_nc = average_nc - min(average_nc);
average_rl = average_rl - min(average_rl);
average_emp = average_emp - min(average_emp);

% save the average steps

% kind = 'GAS';
% emp = "DataAnalysis\General\"+speed+"\"+kind+"\emp";
% rl = "DataAnalysis\General\"+speed+"\"+kind + "\rl";
% withoutExo = "DataAnalysis\General\"+speed+"\"+kind + "\ne";
% withoutControl = "DataAnalysis\General\"+speed+"\"+kind + "\nc";
% 
% save(emp, "average_emp");
% save(rl, "average_rl");
% save(withoutExo, "average_ne");
% save(withoutControl, "average_nc");

% create a graph for demonstration

set(groot,'defaultAxesTickLabelInterpreter','latex'); set(groot,'defaulttextinterpreter','latex'); set(groot,'defaultLegendInterpreter','latex');
withoutExo_color = '#77A7BB';
withoutControl_color = '#DE9536';
emp_color = '#AFC86A';
rl_color = '#7B6093';
fontsize = 18.5;

plot(t_new_ne/(t_new_ne(end))*100, average_ne, Color=withoutExo_color, LineWidth=2.5);
hold on
plot(t_new_emp/(t_new_emp(end))*100, average_emp, Color=emp_color, LineWidth=2.5);
plot(t_new_rl/(t_new_rl(end))*100, average_rl, Color=rl_color, LineWidth=2.5);
plot(t_new_nc/(t_new_nc(end))*100, average_nc, Color = withoutControl_color, LineWidth=2.5);
legende = {'NE','Emp','RL', 'NC'};
ylim([-0,1.6])
xline(65, LineStyle="--",Color=withoutControl_color,LineWidth=2.5);
xline(64.25, LineStyle="--",Color=withoutExo_color,LineWidth=2.5);
xline(64.75, LineStyle="--",Color=emp_color,LineWidth=2.5);
xline(64, LineStyle="--",Color=rl_color,LineWidth=2.5);
xlabel("Gait Cycle [\%]");
ylabel('GAS','Interpreter','latex','FontSize', fontsize)
set(gcf,'Units','centimeters','PaperUnits','centimeters'); Pos=get(gcf,'Position'); 
set(gcf,'PaperPosition', [0 0 Pos(3) Pos(4)],'PaperSize', [Pos(3) Pos(4)]); 
set(gca,'FontSize',fontsize);
% legend(legende, 'Location','northeast');
hold off

% saveas(gcf,['GAS_Act','.pdf']);

%
% Create Plot

gas_nc = load("DataAnalysis\General\PWS\GAS\nc.mat","average_nc").average_nc;
gas_ne = load("DataAnalysis\General\PWS\GAS\ne.mat","average_ne").average_ne;
gas_rl = load("DataAnalysis\General\PWS\GAS\rl.mat","average_rl").average_rl;
gas_emp = load("DataAnalysis\General\PWS\GAS\emp.mat","average_emp").average_emp;

vas_nc = load("DataAnalysis\General\PWS\VAS\nc.mat","average_nc").average_nc;
vas_ne = load("DataAnalysis\General\PWS\VAS\ne.mat","average_ne").average_ne;
vas_rl = load("DataAnalysis\General\PWS\VAS\rl.mat","average_rl").average_rl;
vas_emp = load("DataAnalysis\General\PWS\VAS\emp.mat","average_emp").average_emp;

sol_nc = load("DataAnalysis\General\PWS\SOL\nc.mat","average_nc").average_nc;
sol_ne = load("DataAnalysis\General\PWS\SOL\ne.mat","average_ne").average_ne;
sol_rl = load("DataAnalysis\General\PWS\SOL\rl.mat","average_rl").average_rl;
sol_emp = load("DataAnalysis\General\PWS\SOL\emp.mat","average_emp").average_emp;

rf_nc = load("DataAnalysis\General\PWS\RF\nc.mat","average_nc").average_nc;
rf_ne = load("DataAnalysis\General\PWS\RF\ne.mat","average_ne").average_ne;
rf_rl = load("DataAnalysis\General\PWS\RF\rl.mat","average_rl").average_rl;
rf_emp = load("DataAnalysis\General\PWS\RF\emp.mat","average_emp").average_emp;

bf_nc = load("DataAnalysis\General\PWS\BF\nc.mat","average_nc").average_nc;
bf_ne = load("DataAnalysis\General\PWS\BF\ne.mat","average_ne").average_ne;
bf_rl = load("DataAnalysis\General\PWS\BF\rl.mat","average_rl").average_rl;
bf_emp = load("DataAnalysis\General\PWS\BF\emp.mat","average_emp").average_emp;

max_gas_nc = max(gas_nc);
max_gas_ne = max(gas_ne);
max_gas_rl = max(gas_rl);
max_gas_emp = max(gas_emp);

max_vas_nc = max(vas_nc);
max_vas_ne = max(vas_ne);
max_vas_rl = max(vas_rl);
max_vas_emp = max(vas_emp);

max_sol_nc = max(sol_nc);
max_sol_ne = max(sol_ne);
max_sol_rl = max(sol_rl);
max_sol_emp = max(sol_emp);

max_rf_nc = max(rf_nc);
max_rf_ne = max(rf_ne);
max_rf_rl = max(rf_rl);
max_rf_emp = max(rf_emp);

max_bf_nc = max(bf_nc);
max_bf_ne = max(bf_ne);
max_bf_rl = max(bf_rl);
max_bf_emp = max(bf_emp);


X=categorical({' HAM ',' RF ',' VAS ', ' GAS ', ' SOL '});
X=reordercats(X,{' HAM ',' RF ',' VAS ', ' GAS ', ' SOL  '});

y = [max_bf_nc-1 max_bf_emp-1 max_bf_rl-1;
     max_rf_nc-1 max_rf_emp-1 max_rf_rl-1;
     max_vas_nc-1 max_vas_emp-1 max_vas_rl-1;
     max_gas_nc-1 max_gas_emp-1 max_gas_rl-1;
     max_sol_nc-1 max_sol_emp-1 max_sol_rl-1 ]*100;

legende = {'WithoutControl',  'Experimental', 'AI-developed' };
h_bar = bar(X,y, 'FaceColor','flat');
hold on
text(h_bar(2).XEndPoints(4) ,0,'*','HorizontalAlignment','center',...
    'VerticalAlignment','bottom','FontSize',fontsize)
% text(h_bar(3).XEndPoints(1) ,0,'*','HorizontalAlignment','center',...
        % 'VerticalAlignment','bottom', 'FontSize',fontsize)
ylabel('Peak Value EMG [\%]','Interpreter','latex','FontSize', fontsize)
% set(h_bar(1),'FaceColor',withoutExo_color);
set(h_bar(1),'FaceColor',withoutControl_color);
set(h_bar(2),'FaceColor',emp_color);
set(h_bar(3),'FaceColor',rl_color);
set(gcf,'Units','centimeters','PaperUnits','centimeters'); Pos=get(gcf,'Position'); 
set(gcf,'PaperPosition', [0 0 Pos(3) Pos(4)],'PaperSize', [Pos(3) Pos(4)]); 
set(gca,'FontSize',fontsize);
% ylim([-inf,0.6])
axis("padded")
% legend(legende, 'Location','northeast');
hold off

saveas(gcf,['Max_values_PWS','.pdf']);

%
%
% MAV

mav_gas_nc = mean(abs(gas_nc));
mav_gas_ne = mean(abs(gas_ne));
mav_gas_rl = mean(abs(gas_rl));
mav_gas_emp = mean(abs(gas_emp));

mav_vas_nc = mean(abs(vas_nc));
mav_vas_ne = mean(abs(vas_ne));
mav_vas_rl = mean(abs(vas_rl));
mav_vas_emp = mean(abs(vas_emp));

mav_sol_nc = mean(abs(sol_nc));
mav_sol_ne = mean(abs(sol_ne));
mav_sol_rl = mean(abs(sol_rl));
mav_sol_emp = mean(abs(sol_emp));

mav_rf_nc = mean(abs(rf_nc));
mav_rf_ne = mean(abs(rf_ne));
mav_rf_rl = mean(abs(rf_rl));
mav_rf_emp = mean(abs(rf_emp));

mav_bf_nc = mean(abs(bf_nc));
mav_bf_ne = mean(abs(bf_ne));
mav_bf_rl = mean(abs(bf_rl));
mav_bf_emp = mean(abs(bf_emp));

X=categorical({' HAM ',' RF ',' VAS ', ' GAS ', ' SOL '});
X=reordercats(X,{' HAM ',' RF ',' VAS ', ' GAS ', ' SOL  '});

y = [(mav_bf_nc/mav_bf_ne)-1 (mav_bf_emp/mav_bf_ne)-1 (mav_bf_rl/mav_bf_ne)-1;
     (mav_rf_nc/mav_rf_ne)-1 (mav_rf_emp/mav_rf_ne)-1 (mav_rf_rl/mav_rf_ne)-1;
     (mav_vas_nc/mav_vas_ne)-1 (mav_vas_emp/mav_vas_ne)-1 (mav_vas_rl/mav_vas_ne)-1;
     (mav_gas_nc/mav_gas_ne)-1 (mav_gas_emp/mav_gas_ne)-1 (mav_gas_rl/mav_gas_ne)-1;
     (mav_sol_nc/mav_sol_ne)-1 (mav_sol_emp/mav_sol_ne)-1 (mav_sol_rl/mav_sol_ne)-1 ]*100;

legende = {'NC',  'Emp', 'RL' };
h_bar = bar(X,y, 'FaceColor','flat');
hold on
% text(h_bar(2).XEndPoints(1) ,0,'*','HorizontalAlignment','center',...
    % 'VerticalAlignment','bottom','FontSize',fontsize)
% text(h_bar(3).XEndPoints(3) ,75,'*','HorizontalAlignment','center',...
        % 'VerticalAlignment','bottom', 'FontSize',fontsize)
ylabel('MAV','Interpreter','latex','FontSize', fontsize)
% set(h_bar(1),'FaceColor',withoutExo_color);
set(h_bar(1),'FaceColor',withoutControl_color);
set(h_bar(2),'FaceColor',emp_color);
set(h_bar(3),'FaceColor',rl_color);
set(gcf,'Units','centimeters','PaperUnits','centimeters'); Pos=get(gcf,'Position'); 
set(gcf,'PaperPosition', [0 0 Pos(3) Pos(4)],'PaperSize', [Pos(3) Pos(4)]); 
set(gca,'FontSize',fontsize);
% ylim([-90,5])
axis("padded")
legend(legende, 'Location','northwest');
hold off

% saveas(gcf,['MAV','.pdf']);