
%% Script principale modificato
% seleziona ingressi e uscite, divide e scala le misure, setta intervalli di validazione/identificazione,
% e elimina le misure corrotte (o ritenute tali!!)
%clear all; 
close all; 
clc;


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




%% 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

% 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));


%% Imposta il tempo di campionamento e preleva ingressi ed uscite

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

%---------------------------------------------------------------
%% INGRESSI 
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

% USCITE
T = double(data.T(pos_y,1:camp:end));
%---------------------------------------------------------------
%---------------------------------------------------------------
%% prova di PCA sulle uscite (T)
figure
plot(T')
title(['Temperature interne (uscite)']);
xlabel('campioni');
ylabel('°C');
axis tight;
grid


percTeor=0.998;
%trova gli indici dei sensori con il minimo errore quadratico medio
%rispetto alla proiezione dei sensori sulle componenti principali
Ttmp=T;
[Ybas,n,percR,BaseY]=KLT(Ttmp,percTeor);

figure
plot(Ybas')
title(['PC delle uscite con ' num2str(percR) '% di informazione spiegata']);
xlabel('campioni');
axis tight;
grid

%
Coeff=T*pinv(Ybas);
figure
stem(Coeff)
title(['Loadings delle prime ' num2str(n) ' componenti principali' ]);
xlabel('numero sensore');
axis tight;
grid

%% PROVA:selezionare i sensori con i coefficienti più alti
%sensori n°:2,6,22,30,38,47
for i=1:n
    [Val,Nsens(i)]=max(abs(Coeff(:,i)));
    Coeff=[Coeff(1:Nsens(i)-1,:);Coeff(Nsens(i)+1:size(Coeff,1),:)];
end

Nsens1=[2,6,22,30,38,47];

Tsel=T(Nsens,:);
figure
plot(Tsel')
title(['Uscite selezionate date le ' num2str(n) ' componenti principali']);
xlabel('campioni');
ylabel('°C');
axis tight;
grid

% E' CMQ UN APPROCCIO GREEDY MA SEMBRA VADA BENE CMQ PER DESCRIVERE IL
% COMPORTAMENTO DEL SISTEMA COMPLESSIVO

%% STRUTTURA 3D DELLA DISPOSIZIONE DEI SENSORI
nomi = data.name(pos_y,:);
nomiSel=nomi(Nsens,:);
stampamotes(nomiSel,data)



%% KALMAN DOPO LA SELEZIONE DEI SENSORI

% x0 = pinv(C)*ones(length(C),1)*mean(Yval(1,Nsens)); % se nn conosco tutte le misure, uso la media
%x0 = pinv(C)*y(:,1); % se conosco le misure approssimo meglio che posso x0 (least squares)
YvalKF = fK( A,B,C,K,VarErr,Base'*UvalC,YvalC,Nsens,eye(size(A)),x0 );

%escludo i primi 50 campioni che sono relativi al tempo di assestamento del filtro

figure 
stem(sum((YvalC(:,50:size(YvalC,2))'-YvalKF(:,50:size(YvalC,2))').^2)./size(YvalC,1))
title('MSE dei sensori rispetto alle uscite di KF')
xlabel('numero sensore')
grid


figure
subplot(2,1,1)
plot(YvalC')
grid
title('Uscite di validazione')
xlabel('campioni');
ylabel('°C');
axis tight
subplot(2,1,2)
plot(YvalKF')
grid
title('Uscite con KF')
xlabel('campioni');
ylabel('°C');
axis tight


%% PROVE VARIE

% YbasTot(1,:)=Ybas(1,:);
% [mT,Ttmp]=centra(T);
% [I,eqmP]=trovaMinEQ(Ttmp,Ybas);
% 
% 
% Itot(1)=I(1);
% eqmTot(1)=eqmP(1);
% 
% 
% 
% for i=2:n
%     Ttmp=[Ttmp(1:I(1)-1,:);Ttmp(I(1)+1:size(Ttmp,1),:)];
%     [Ybas,ntmp,percR]=KLT(Ttmp,percTeor);    
%     YbasTot(i,:)=Ybas(1,:);
%     [I,eqmP]=trovaMinEQ(Ttmp,Ybas);       
%     %N.B.:se prendo un sensore di indice > di Itot(i-1) su Ttmp
%     %     =>devo incrementare di 1 il num di sensore scelto su T
%     [r,c,k]=find(Itot<=I(1));
%     Itot(i)=I(1)+sum(k);
%     eqmTot(i)=eqmP(1);
%     
% %     figure(i)
% %     plot(T(Itot(i),:))
% %     hold on
% %     plot(Ttmp(I(1),:),'--r')
% %     grid
% %     hold off
%     
%     %c'è qualche errore, perchè il Ttmp selezionato e T non corrispondono!!    
%     
% end
% 
% 
% %% ANDAMENTI TEMPORALI di Ybas e dei relativi sensori scelti
% % N.B.B.B:non sono reali andamenti temporali, perche' scalati secondo la
% %         dev std di ogniuno e a media nulla!!!!
% [mT,Ttmp]=centra(T);
% for i=1:n
%     figure(i)
%     plot(YbasTot(i,:))
%     hold on
%     plot(Ttmp(Itot(i),:),'--r')
%     grid
%     hold off
% end
% 
% %STRUTTURA 3D DELLA DISPOSIZIONE DEI SENSORI
% nomi = data.name(pos_y,:);
% nomiSel=nomi(Itot,:);
% stampamotes(nomiSel,data)



