%% MAKE a second script for radial basis functions predictions, online prediction and machine learning
clear 
close all 
clc
addpath('/Users/simenbjorvand/Documents/OneDrive - NTNU/project copy/matlab')



%% load the data here



time=[];
z=[];
wGlab = [];
wLlab = [];
Pin = [];
Prt = [];
Prb = [];

t=true;
NLs = [];

index =1:63;

for i = index
    filePath = ['/Users/simenbjorvand/Documents/OneDrive - NTNU/project copy/matlab/filterData/timeSeries',num2str(i)];
    load(filePath)
    if ((data(2,1)-data(1,1))<0.8)
        NLs=[NLs;length(data(1:10:end,1))];
        if t
            time = [time;data(1:10:end,1)];
            t=false;
        else
            time = [time;(time(end)+data(1:10:end,1))];
        end
        z = [z;data(1:10:end,4)];
        wGlab = [wGlab;data(1:10:end,2)];
        wLlab = [wLlab;data(1:10:end,3)];
        Pin = [Pin;data(1:10:end,5)]; %use pressure in buffertank as inlet pressure
        Prt = [Prt;data(1:10:end,8)];
        Prb = [Prb;data(1:10:end,7)];
    else
        NLs=[NLs;length(data(:,1))];
        if t
            time = [time;data(:,1)];
            t=false;
        else
            time = [time;(time(end)+data(:,1))];
        end
        z = [z;data(:,4)];
        wGlab = [wGlab;data(:,2)];
        wLlab = [wLlab;data(:,3)];
        Pin = [Pin;data(:,5)]; %use pressure in buffertank as inlet pressure
        Prt = [Prt;data(:,8)];
        Prb = [Prb;data(:,7)];
    end
    
end


%%
NL = length(z);
Nd=80;


%% do some initial simulations using simplyfied model
wGsim = NaN(Nd,1);
wLsim = NaN(Nd,1);

u = [z,Pin,Prt,Prb];
Nu = length(u(1,:));

simulate=true;
if(simulate)
    h = waitbar(0,'Simulating');
    for i = 1:Nd
        y = simplifiedModel(u(i,:));
        wGsim(i) = y(1);
        wLsim(i) = y(2);
        waitbar(i/Nd);
    end
    close(h)
end

initLabData = [wGlab(1:Nd),wLlab(1:Nd)];
initSimData = [wGsim,wLsim];




%% model further down from here

errorNN = radialBasisNetwork2(initLabData,initSimData,u(1:Nd,:),Nd,5);
predErrG=[];
predErrL=[];
yhat = nan(NL-Nd,2);

h = waitbar(0,'Simulating');
for i = (Nd+1):NL
    ysim = simplifiedModel(u(i,:));
    predictedError = errorNN.predictError(u(i,:));
    yhat((i-Nd),:)=ysim+predictedError;
    errorNN = errorNN.uppdateWithNewData(u(i,:),[wGlab(i),wLlab(i)],ysim);
    
    predErrG=[predErrG;predictedError(1)];
    predErrL=[predErrL;predictedError(2)];
    wGsim = [wGsim;ysim(1)];
    wLsim = [wLsim;ysim(2)];  
    
    waitbar(i/NL)
end
close(h)


errorG = wGlab((Nd+1):end,:)-wGsim((Nd+1):end,:);
errorL = wLlab((Nd+1):end,:)-wLsim((Nd+1):end,:);

t=true;
for i = 1:length(index)
    %% calculate the total mse between predicted error and model error
    itStart = NL-sum(NLs(i:end))+1;
    if t
        itStart = Nd+1;
    end
    
    itEnd = sum(NLs(1:i));
    
    
    %% calculate some mean and covarinace matrices here
    up = u(itStart:itEnd,:);
    yp = [wGlab(itStart:itEnd),wLlab(itStart:itEnd)];
    ysp = [wGsim(itStart:itEnd),wLsim(itStart:itEnd)];
    epsilon = yp-ysp;
    yhatp = yhat((itStart-Nd:itEnd-Nd),:);
    epsilonhat = [predErrG(itStart-Nd:itEnd-Nd),predErrL(itStart-Nd:itEnd-Nd)];

    %calculate mean of these 
    mean_up = mean(up);
    mean_yp = mean(yp);
    mean_ysp = mean(ysp);
    mean_yhat = mean(yhatp);
    mean_epsilon = mean(epsilon);
    mean_epsilonhat = mean(epsilonhat);

    %calculate covarinace matrix
    cov_up = cov(up);
    cov_yp = cov(yp);
    cov_ysp = cov(ysp);
    cov_yhat = cov(yhatp);
    cov_epsilon = cov(epsilon);
    cov_epsilonhat = cov(epsilonhat);

    % extract the variance from the covarinace
    var_up = diag(cov_up)';
    var_yp = diag(cov_yp)';
    var_ysp = diag(cov_ysp)';
    var_yhat = diag(cov_yhat)';
    var_epsilon = diag(cov_epsilon)';
    var_epsilonhat = diag(cov_epsilonhat)';

    std_up = sqrt(var_up);
    std_yp = sqrt(var_yp);
    std_ysp = sqrt(var_ysp);
    std_yhat = sqrt(var_yhat);
    std_epsilon = sqrt(var_epsilon);
    std_epsilonhat = sqrt(var_epsilonhat);
    
    if t
        mu_yp = mean_yp;
        mu_yhat = mean_yhat;
        mu_ysp = mean_ysp;
        mu_epsilon = mean_epsilon;
        mu_epsilonhat = mean_epsilonhat;
        
        s_yp = std_yp;
        s_yhat = std_yhat;
        s_ysp = std_ysp;
        s_epsilon = std_epsilon;
        s_epsilonhat = std_epsilonhat;
        t = false;
    else
        mu_yp =[mu_yp;mean_yp];
        mu_yhat = [mu_yhat;mean_yhat];
        mu_ysp = [mu_ysp;mean_ysp];
        mu_epsilon = [mu_epsilon;mean_epsilon];
        mu_epsilonhat = [mu_epsilonhat;mean_epsilonhat];
        
        s_yp = [s_yp;std_yp];
        s_yhat = [s_yhat;std_yhat];
        s_ysp = [s_ysp;std_ysp];
        s_epsilon = [s_epsilon;std_epsilon];
        s_epsilonhat = [s_epsilonhat;std_epsilonhat];
    end
        

end


figure(1)
subplot(211)
plot(2:7,2:7)
hold on
% errorbar(mu_yp(:,1),mu_yhat(:,1),-2*s_yhat(:,1),2*s_yhat(:,1),-2*s_yp(:,1),2*s_yp(:,1),'x')
scatter(mu_yp(:,1),mu_yhat(:,1),'x','MarkerEdgeColor','r')
scatter(mu_yp(:,1),mu_ysp(:,1),'o','MarkerEdgeColor','b')
%xlim([4.485,4.515])
xlim([3.9,4.6])
ylim([2,7])
title('Gas flowrate')
xlabel('Target')
ylabel('Prediction')
legend({'','$\overline{\hat{y}}$','$\overline{y^{*}}$'},'Interpreter','latex')
grid on
hold off

subplot(212)
plot(1:6,1:6)
hold on
% errorbar(mu_yp(:,2),mu_yhat(:,2),-2*s_yhat(:,2),2*s_yhat(:,2),-2*s_yp(:,2),2*s_yp(:,2),'x')
scatter(mu_yp(:,2),mu_yhat(:,2),'x','MarkerEdgeColor','r')
scatter(mu_yp(:,2),mu_ysp(:,2),'o','MarkerEdgeColor','b')
xlim([1,6])
ylim([1,6])
title('Liquid flowrate')
xlabel('Target')
ylabel('Prediction')
legend({'','$\overline{\hat{y}}$','$\overline{y^{*}}$'},'Interpreter','latex')
grid on
hold off

%%
figure(3)
subplot(211)
plot(2:7,2:7)
hold on
% errorbar(mu_yp(:,1),mu_yhat(:,1),-2*s_yhat(:,1),2*s_yhat(:,1),-2*s_yp(:,1),2*s_yp(:,1),'x','MarkerEdgeColor','r','Color','r')
% errorbar(mu_yp(:,1),mu_ysp(:,1),-2*s_ysp(:,1),2*s_ysp(:,1),-2*s_yp(:,1),2*s_yp(:,1),'o','MarkerEdgeColor','b','Color','b')
errorbar(mu_yp(:,1),mu_yhat(:,1),-2*s_yhat(:,1),2*s_yhat(:,1),'x','MarkerEdgeColor','r','Color','r')
errorbar(mu_yp(:,1),mu_ysp(:,1),-2*s_ysp(:,1),2*s_ysp(:,1),'o','MarkerEdgeColor','b','Color','b')
errorbar(mu_yp(:,1),mu_yp(:,1),-2*s_yp(:,1),2*s_yp(:,1),'o','MarkerEdgeColor','g','Color','g')

xlim([4.485,4.515])
ylim([0,10])
xlabel('labData')
ylabel('Prediction')
legend({'','$\overline{\hat{y}}$','$\overline{y^{*}}$','$\overline{y}$'},'Interpreter','latex')
grid on
hold off

subplot(212)
plot(1:6,1:6)
hold on
% errorbar(mu_yp(:,2),mu_yhat(:,2),-2*s_yhat(:,2),2*s_yhat(:,2),-2*s_yp(:,2),2*s_yp(:,2),'x','MarkerEdgeColor','r','Color','r')
% errorbar(mu_yp(:,2),mu_ysp(:,2),-2*s_ysp(:,2),2*s_ysp(:,2),-2*s_yp(:,2),2*s_yp(:,2),'o','MarkerEdgeColor','b','Color','b')
errorbar(mu_yp(:,2),mu_yhat(:,2),-2*s_yhat(:,2),2*s_yhat(:,2),'x','MarkerEdgeColor','r','Color','r')
errorbar(mu_yp(:,2),mu_ysp(:,2),-2*s_ysp(:,2),2*s_ysp(:,2),'o','MarkerEdgeColor','b','Color','b')
errorbar(mu_yp(:,2),mu_yp(:,2),-2*s_yp(:,2),2*s_yp(:,2),'o','MarkerEdgeColor','g','Color','g')

xlim([1,6])
ylim([1,6])
xlabel('labData')
ylabel('Prediction')
legend({'','$\overline{\hat{y}}$','$\overline{y^{*}}$','$\overline{y}$'},'Interpreter','latex')
grid on
hold off


%%





%% make some performance plots here: 
figure(10000)
subplot(221)
scatter(wGlab(Nd+1:end),yhat(:,1))
hold on
errorbar(mu_yp(:,1),mu_yhat(:,1),-2*s_yhat(:,1),2*s_yhat(:,1),-2*s_yp(:,1),2*s_yp(:,1),'x')
xlabel('labData')
ylabel('Prediction')
grid on

subplot(222)
scatter(errorG,predErrG)
hold on
errorbar(mu_epsilon(:,1),mu_epsilonhat(:,1),-2*s_epsilonhat(:,1),2*s_epsilonhat(:,1),-2*s_epsilon(:,1),2*s_epsilon(:,1),'x')
xlabel('model error')
ylabel('Prediction error')
grid on

subplot(223)
scatter(wLlab(Nd+1:end),yhat(:,2))
hold on
errorbar(mu_yp(:,2),mu_yhat(:,2),-2*s_yhat(:,2),2*s_yhat(:,2),-2*s_yp(:,2),2*s_yp(:,2),'x')
xlabel('labData')
ylabel('Prediction')
grid on

subplot(224)
scatter(errorL,predErrL)
hold on
errorbar(mu_epsilon(:,2),mu_epsilonhat(:,2),-2*s_epsilonhat(:,2),2*s_epsilonhat(:,2),-2*s_epsilon(:,2),2*s_epsilon(:,2),'x')
xlabel('model error')
ylabel('Prediction error')
grid on