%% prova PCR su ingressi
% idea:fare pcr su ogni gruppo di ingressi
% gruppo 1: temperatura esterna
% gruppo 2: radiazione attiva
% gruppo 3: radiazione totale
% N.B: forse i gruppi 2 e 3 vanno uniti in un unico gruppo!!
clear all 
close all
clc

%% Carica i dati e imposta path
load data;

lum = ['051';'007';'137';'033';'132'];   % insieme sensori esterni luminositÃ 
tempE = ['131';'147';'157';'010'];       % insieme sensori esterni temperatura

% tempo di campionamento [in MINUTI]
tempoCampionamento = 10;

% spezza i dati in IDENTIFICAZIONE (frazione*totale) / VALIDAZIONE (1-frazione)*totale
frazione = 0.8; % tra zero e uno

% trovo le posizioni nel database
pos_uT = name2pos(tempE,data); 
pos_uL = name2pos(lum,data);

% e le tolgo dall'insieme delle misure
pos_y = setdiff([1:length(data.name)],[name2pos(tempE,data) name2pos(lum,data)]);

% tolgo il 102 e il 154 causa misure corrotte
pos_y = setdiff(pos_y,name2pos(['102';'154'],data));

%% INGRESSI 

camp = round(tempoCampionamento*2); % tempo campionamento dati: 10 minuti
% SOTTOCAMPIONA 

Text = double(data.T(pos_uT,1:camp:end));
ARext = double(data.aR(pos_uL,1:camp:end))./500000; % scala la luminosita' attiva
TRext = double(data.tR(pos_uL,1:camp:end))./50000;  % scala la luminosita' totale
U=[Text;ARext;TRext];%vettore con tutti gli ingressi


% USCITE
T = double(data.T(pos_y,1:camp:end));

%% PARTIZIONAMENTO DEI DATI

% identificazione
Yid = T(:,1:floor(size(T,2)*frazione)); 
Uid = U(:,1:floor(size(U,2)*frazione));
% validazione
Yval = T(:,ceil(size(T,2)*frazione):end);
Uval = U(:,ceil(size(U,2)*frazione):end);
figure
plot(Yid')


figure
subplot(2,1,1)
plot(Uid(1:4,:)')
title('Temperature esterne','fontsize',11);
%xlabel('campioni');
ylabel('°C');
axis tight;
grid
subplot(2,1,2)
plot(Uid(5:end,:)')
title('Luminositá attive e solari totali','fontsize',11);
xlabel('campioni');
ylabel('lux');
axis tight;
grid

%% SELEZIONE COMPONENTI PRINCIPALI DI U(t)
percT=0.99;

[mYid,YidC]=centra(Yid);
[mUid,UidC]=centra(Uid);
[mYval,YvalC]=centra(Yval);
[mUval,UvalC]=centra(Uval);


[Ucp,n,percR,Base]=KLT(Uid,percT);

figure
plot(Ucp')
title(['PC degli ingressi con ' num2str(percR) '% di informazione spiegata']);
xlabel('campioni');
axis tight;
grid

%% IDENTIFICAZIONE
p=7;
f=4;

[A,B,C,D,K,VarErr] = PBSID_opt_ARX_weighted(YidC,Ucp,[ p f ]);
%selezionando un sistema di ordine 30/35 vengono molto bene gli errori medi

%stima lo stato iniziale
m_idss = idss(A,B,C,D,K);
d1 = iddata(YvalC',(Base'*UvalC)');
[e,x0] = pe(m_idss,d1,'e');


S = ss(A,B,C,D,1);
Yca=lsim(S,Base'*UvalC,0:length(Uval)-1,x0);
%Yca=Yca+(mYval'*ones(size(Yca,2),size(Yca,1)))';

Srsd=idss(A,B,C,D,K);

%% PRESENTAZIONE DEI RISULTATI
% figure
% subplot(2,1,1)
% plot(Yca)
% ylabel('Y catena aperta')
% grid
% subplot(2,1,2)
% plot(YvalC')
% ylabel('Y validazione')
% grid
% 
% figure
% plot(Yca-YvalC')
% grid
% 
% figure
% stem(sum((Yca-YvalC').^2)./size(Yca,1))
% grid


% z = iddata(YidC',Ucp');
% S=n4sid(z);
% Yca1=lsim(S,Base(:,1:n)'*Uval,0:length(Uval)-1);
% % Yca=Yca+mYval'*ones(size(Yca,1),size(Yca,2));
% figure
% subplot(2,1,1)
% xlabel('Y catena aperta 1')
% plot(Yca1)
% grid
% subplot(2,1,2)
% xlabel('Y validazione')
% plot(YvalC')
% grid

