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

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



%% 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 secondo la dinamica ritenuta corretta

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

percTeor=0.98;
%trova gli indici dei sensori con il minimo errore quadratico medio
%rispetto alla proiezione dei sensori sulle componenti principali
[Ybas,n,percR,Tsd]=KLT1(T,percTeor);
YbasTot(1,:)=Ybas(1,:);
[mT,Tc]=centra(T);
[I,eqmP]=trovaMinEQ(Tc,Ybas);
Itot(1)=I(1);
eqmTot(1)=eqmP(1);
Ttmp=T;

plot(Ybas')
figure
plot(Tc')
figure
plot(Tsd')

% for i=2:n
%     Ttmp=[Ttmp(1:I(1)-1,:);Ttmp(I(1)+1:size(Ttmp,1),:)];    
%     [Ybas,ntmp,percR]=KLT1(Ttmp,percTeor);
%     YbasTot(i,:)=Ybas(1,:);
%     [I,eqmP]=trovaMinEQ(Ttmp,Ybas);
%     It(i)=I(1);
%     %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    
%     eqmTot(i)=eqmP(1);
% end
% 
% for i=2:length(It)
%     [r,c,k]=find(It(1:i-1)<=It(i));
%     Itot(i)=It(i)+sum(k);
% 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!!!!
% for i=1:n
%     figure(i)
%     plot(YbasTot(i,:))
%     hold on
%     plot(T(Itot(i),:),'--r')
%     grid
%     hold off
% end

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



