% funzione per il calcolo del numero di passi minimo in funzione dell'erore
% medio di uscita
function [C_pls,e] = OrARX(n,n_VL)
tempoCampionamento = 2;
m = 4000 ; %numero di campioni
mm = 1600;
ky = 100; % valore K-esimo INIZIALE dell'uscita del modello ARX 
frazione=3/4;

%% 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
    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);
    Uid = U(:,1:floor(size(U,2)*frazione)); % frazione degli ingressi di identificazione 14 x 5668
    [mu,Uid] = centra(Uid);
    % validazione
    Yval = T(:,ceil(size(T,2)*frazione):end); %frazione delle uscite di validazione 54 x 1890
    [my,Yval] = centra(Yval);
    Uval = U(:,ceil(size(U,2)*frazione):end); %frazione degli ingressi di validazione  14 x 1890
    [my,Uval] = centra(Uval);
    
     %prove su un set ridotto di ingressi
     %Uid = Uid(1:3,:); %prendo 3 ingressi  
     Yid = Yid(1:10,:); %prendo 3 uscite 
     %Uval = Uval(1:3,:);
     Yval = Yval(1:10,:);
    
    %% Suddivisione delle matrici di ingresso alla funzione PLSF DATI IDENTIFICAZIONE

    ingressiID=zeros(size(Uid,1),n); %inizializzo una matrice che verrà usata per salvare gli n_campioni precedenti all'istante Ky dell'ingressi
    usciteID = zeros(size(Yid,1),n); %inizializzo una matrice che verrà usata per salvare gli n_campioni precedenti all'istante Ky dell'uscita 
    for j=0:m-1                      %prendo m campioni 
        for i=1:n
            usciteID(:,i) = Yid(:,ky+j-i);  % prendo n_campioni precedenti al punto iniziale ky con n ordine del modello ARX
            usciteIDb(i,:)= usciteID(:,i)'; %li metto come vettore riga, alla fine avrò n righe  x sensori(colonne)
            ingressiID(:,i)= Uid(:,ky+j-i); %a partire da ky prendo n istanti precedenti con n ordine del modello ARX
            ingressiIDb(i,:)= ingressiID(:,i)'; 
        end
        X(:,:,j+1) = [usciteIDb,ingressiIDb]; % per ogni misura m_esima da Ky salvo in X gli n_campioni precedenti degli ingressi
        Y(:,:,j+1) = Yid(:,ky+j)';            % per ogni misura m_esima da Ky salvo in Y gli m_campioni successivi a partire da Ky
    end

    for j = 1:m            % con questi cicli si sistemano le matrici nel modo descritto da Qui, in Partial Least Squares Regression for recursive System Identification
        X1 = X(:,1:10,j);  % i cicli sono fatti su 'misura' per gli ingressi e le uscite usate ossia 14 ingressi e 54 uscite.                        
        X11 = zeros(1,size(X1,1)*size(X1,2));
        for i=1:size(X1,1) %ciclo per il numero di righe                                            
            X11(1,(i-1)*size(X1,2)+1:i*size(X1,2))= X1(i,:);%smanetto cn gli indici...ho controllato e funziona
        end    
        X_1(j,:) = X11(1,:); 
    end   

    for j = 1:m  % con questi cicli si sistemano le matrici nel modo descritto da Qui in Partial Least Squares Regression for recursive System Identification
        X2 = X(:,11:end,j); % i cicli sono fatti su misura per gli ingressi e le uscite usate ossia 14 ingressi e 54 uscite.        
        %alloco il vettore
        X22 = zeros(1,size(X2,1)*size(X2,2));
        for i=1:size(X2,1)%ciclo per il numero di righe
            %si lvora un pò con gli indici...
            X22(1,(i-1)*size(X2,2)+1:i*size(X2,2))= X2(i,:);
        end    
        X_2(j,:) = X22(1,:); 
    end   
    XplsID = [X_1,X_2]; %matrice di ingresso per l'algoritmo di PLS
    for j = 1:m
        YID(j,:)= Y(1,:,j);
    end
    %% sistemazione delle matrici per VERIFICA modello IDENTIFICATO VALIDAZIONE
    ingressiVAL=zeros(size(Uval,1),n);%inizializzo una matrice che verrà usata per salvare gli n_campioni precedenti all'istante Ky dell'ingressi
    usciteVAL = zeros(size(Yval,1),n);%inizializzo una matrice che verrà usata per salvare gli n_campioni precedenti all'istante Ky dell'uscita      
    for j=0:mm-1                      %prendo m campioni 
        for i=1:n
            usciteVAL(:,i) = Yval(:,ky+j-i); % prendo n_campioni precedenti al punto iniziale ky con n ordine del modello ARX
            usciteVALb(i,:)= usciteVAL(:,i)';      %li metto come vettore riga, alla fine avrò n righe  x sensori(colonne)
            ingressiVAL(:,i)= Uval(:,ky+j-i);   %a partire da ky prendo n istanti precedenti con n ordine del modello ARX
            ingressiVALb(i,:)= ingressiVAL(:,i)'; 
        end
        XX(:,:,j+1) = [usciteVALb,ingressiVALb]; % per ogni misura m_esima da Ky salvo in X gli n_campioni precedenti degli ingressi
        YY(:,:,j+1) = Yval(:,ky+j)';      % per ogni misura m_esima da Ky salvo in Y gli m_campioni successivi a partire da Ky
    end

    for j = 1:mm            % con questi cicli si sistemano le matrici nel modo descritto da Qui, in Partial Least Squares Regression for recursive System Identification
        X3 = XX(:,1:10,j);  % i cicli sono fatti su 'misura' per gli ingressi e le uscite usate ossia 14 ingressi e 54 uscite.                        
        X33 = zeros(1,size(X3,1)*size(X3,2));
        for i=1:size(X3,1) %ciclo per il numero di righe                                            
            X33(1,(i-1)*size(X3,2)+1:i*size(X3,2))= X3(i,:);%smanetto cn gli indici...ho controllato e funziona
        end    
        X_3(j,:) = X33(1,:); 
    end   

    for j = 1:mm  % con questi cicli si sistemano le matrici nel modo descritto da Qui in Partial Least Squares Regression for recursive System Identification
        X4 = XX(:,11:end,j); % i cicli sono fatti su misura per gli ingressi e le uscite usate ossia 14 ingressi e 54 uscite.        
        %alloco il vettore
        X44 = zeros(1,size(X4,1)*size(X4,2));
        for i=1:size(X4,1)%ciclo per il numero di righe
            %smanetto cn gli indici...ho controllato e funziona
            X44(1,(i-1)*size(X4,2)+1:i*size(X4,2))= X4(i,:);
        end    
        X_4(j,:) = X44(1,:); 
    end   
    XplsVAL = [X_3,X_4]; %matrice di ingresso per l'algoritmo di PLS
    for j = 1:mm
        YVAL(j,:)= YY(1,:,j);
    end
    %% Calcolo del modello ARX attraverso PLS

    [T1 P1 U1 Q1 W1 B1] = PLSF2(XplsID,YID,n_VL); % calcolo le matrici
    [nXplsID mXplsID] = size(XplsID);
   
    C_pls = pinv(P1*P1')*P1*B1*Q1'; % matrice dei coefficienti del modello ARX

    Y_ARX_ID = XplsID*C_pls;   %uscite del modello ARX ottenute con i dati di IDENTIFICAZIONE
    Y_ARX_VAL = XplsVAL*C_pls; %uscite del modello ARX ottenute con i dati di VALIDAZIONE
    for i=1: size(Y_ARX_VAL,2) %calcolo il MSE su dati di validazione.
        e(:,i)=((sum((Y_ARX_VAL(:,i)-YVAL(:,i)).^2))/size(Y_ARX_VAL,1));
    end
   
end