%% Test per la scelta del ritardo migliore 
clc
clear all
close all

n_VL = 5;
maxOrdine = 10;
% ogni ordine scelto 
% si calcolano le matrici dei coefficienti. 

%% Carica i dati e imposta path
     tempoCampionamento = 2;
     frazione=3/4;

     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
    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)); %frazione delle uscite di indentificazione 54 x 5668
    [my,Yid] = centra(Yid);
    Yid = Yid(1:10,:); %prendo 3 uscite 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Uid = U(:,1:floor(size(U,2)*frazione)); % frazione degli ingressi di identificazione 14 x 5668
    [my,Uid] = centra(Uid);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Yval = T(:,ceil(size(T,2)*frazione):end); %frazione delle uscite di validazione 54 x 1890
    [my,Yval]= centra(Yval);
    Yval = Yval(1:10,:);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Uval = U(:,ceil(size(U,2)*frazione):end); 
    [mu,Uval]= centra(Uval);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     model=[];
 for n=1:maxOrdine
     [CPLS,erM(1,:,n)] = OrARX(n,n_VL);
     A = CPLS(1:n*10,:);  %54 uscite
     B = CPLS((n*10)+1:end,:);
     AA(:,:,1)= eye(10);
     for k = 1:n
         AA(:,:,k+1) = A((((k-1)*10))+1:k*10,:)';
         BB(:,:,k+1) = B((((k-1)*14))+1:k*14,:)';   % +1:k*14,: 14 ingressi
     end
    x0 = Yid(:,1:n);
    Ysim = simulazARX(AA,BB,Uid,x0);
   model.type{n} = arx([Ysim' Uid'],[[ones(10).*n],[ones(10,14).*n],[ones(10,14)]]);
end
fp = fpe(model.type{1},model.type{2},model.type{3},model.type{4},model.type{5},model.type{6},model.type{7},model.type{8},model.type{9},model.type{10});
figure;
plot(fp,'r','linewidth',1);
grid on;
h2 = title('FPE di 10 modelli ARX con 14 ingressi e 10 uscite');
set(h2,'interpreter','latex','fontsize',12);
h2 = xlabel('ordine del modello');
set(h2,'interpreter','latex','fontsize',12);
h2 = ylabel('FPE');
set(h2,'interpreter','latex','fontsize',12);
% si graficano gli MSE di ogni sensore per ordine di scelta del modello.
for s=1:5
    for n=1:maxOrdine
        sen(n) = erM(1,s,n);
    end
    figure(s);
    stem(sen);
    grid on
end
