% ciclo main. Si impostano le costanti TEM, HUM, NOR per scegliere se
% identificare la temperatura, l'umidita' e se utilizzare i dati con
% varianza unitaria o no. la media è nulla.
% il file passa in rassegna i vari metodi di identificazione e salva i
% risultati in una struttura dati contenuta nel file Dati_Matrice_id.mat



clc
clear


load Dati_Matrice_chiuso.mat; %  QUI SI CARICANO I DATI PER L'IDENTIFICAZIONE


str=['Dati_Matrice_iperfnonsel_k2ssmod'];%  Nome del file di salvataggio



% NORMALIZZAZIONE DATI PER L'IDENTIFICAZIONE VAR1 med0
% prendo i dati di Dati_Matrice e creo Dati_Matrice_norm che contiene i
% dati normalizzati per poter poi usare il parfor nelle simulazioni

for f=1:1
% UMIDITA' MEDIA 0 VAR 1
for k=1:length(Dati)
  
    y = Dati(k).h.y_org; 
    U = Dati(k).h.Udati_org;   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%      
    varU = var(U);
    meanU = mean(U);

    for i=1:size(U,2)
       u(:,i) = (U(:,i)-meanU(i))/sqrt(varU(i));    
    end
    Dati(k).h.Udati_norm = u(:,:); % <- aggiunta del campo

    varY = var(y);
    meanY = mean(y);
    
    Dati(k).h.y_norm=(y(:,1)-meanY(1))/sqrt(varY(1));  % <- aggiunta del campo
end

save(str,'Dati')

% TEMPERATURA MEDIA 0 VAR 1
for k=1:length(Dati)
  
    y = Dati(k).t.y_org; 
    U = Dati(k).t.Udati_org;   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%   
    varU = var(U);
    meanU = mean(U);

    for i=1:size(U,2)
       u(:,i) = (U(:,i)-meanU(i))/sqrt(varU(i));    
    end
    Dati(k).t.Udati_norm = u(:,:); % <- aggiunta del campo

    varY = var(y);
    meanY = mean(y);
    
    Dati(k).t.y_norm=(y(:,1)-meanY(1))/sqrt(varY(1));  % <- aggiunta del campo    
end

% ora salvo la Dati_Matrice a cui ho aggiunto i due campi.
save(str,'Dati')

% TEMPERATURA MEDIA 0 VAR !=1
% si aggiunge l'ident. sseh sulla temperatura con modello 'yy' e 'ny'
for k=1:length(Dati)
  
    y = Dati(k).t.y_org; 
    U = Dati(k).t.Udati_org;   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    meanU = mean(U);

    for i=1:size(U,2)
       u(:,i) = (U(:,i)-meanU(i));    
    end
    Dati(k).t.Udati= u(:,:); % <- aggiunta del campo

    meanY = mean(y);
    
    Dati(k).t.y = (y(:,1)-meanY(1));  % <- aggiunta del campo    
end

% ora salvo la Dati_Matrice a cui ho aggiunto i due campi.
save(str,'Dati')


% UMIDITA' MEDIA 0 VAR !=1
% si aggiunge l'ident. sseh sulla temperatura con modello 'yy' e 'ny'
for k=1:length(Dati)
  
    y = Dati(k).h.y_org; 
    U = Dati(k).h.Udati_org;   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    meanU = mean(U);

    for i=1:size(U,2)
       u(:,i) = (U(:,i)-meanU(i));    
    end
    Dati(k).h.Udati= u(:,:); % <- aggiunta del campo

    meanY = mean(y);
    
    Dati(k).h.y = (y(:,1)-meanY(1));  % <- aggiunta del campo    
    %%%%%%%%%%%%%%%%%%%%%%%%%%% 
end


save(str,'Dati')

end % for inutile. solo per chiudere questa parte di codice
% fine normalizzazione dei dati.
 

%% Inizio identificazione dei dati:


% contiene una struttura Dati che viene via via aggiornata con i risultati
% dei vari metodi di indentificazione. Poi però viene salvata in un fil
% chiamato Dati_Matrice_id

len = length(Dati);

i_end = 2406;    % fine dei dati utilizzati

Ts=60*5; % tempo di campionamento

TEM = true % mettere true se si vuole effettuare una id sulla temper 
HUM = false % mettere true se si vuole effettuare una id sulla umidità

%Per scegliere algoritmo di identificazione
SSEH = true
N4SID = false
SSGLAR = false
GLAR = false
PEM = false


% NOR 1: uso dati con media 0 e var 1 (con nome_norm)
% NOR 2: uso dati con media 0 e var!=1 (con nome )

NOR=1;
switch NOR
    case 1
        for k=1:len 
            Dati(k).t.y_id = Dati(k).t.y_norm;
            Dati(k).t.Udati_id = Dati(k).t.Udati_norm;   

            Dati(k).h.y_id = Dati(k).h.y_norm;
            Dati(k).h.Udati_id = Dati(k).h.Udati_norm;   
        end
    case 2
        for k=1:len 
            Dati(k).t.y_id = Dati(k).t.y;
            Dati(k).t.Udati_id = Dati(k).t.Udati;   

            Dati(k).h.y_id = Dati(k).h.y;
            Dati(k).h.Udati_id = Dati(k).h.Udati;   
        end
end
% da qui si usano solo i due nuovi campi che contengono i dati con media  e
% varianza a scelta



%% SSEH DI TEMPERATURA %%%%%%%%%%%%%%%%%%%%%%%%%%%

if (TEM && SSEH)
	cont_tmp = 0;
    for k=1:len 

        y = Dati(k).t.y_id; 
        U = Dati(k).t.Udati_id;   


        i_start= 10;  % 800 per ident, 600 per validazione. si scartano gli ultimi.
        i_stop = 1500; % il  max è 1500

        
        
        Up=U(i_start:i_stop,:);     % p per passato
        yp=y(i_start:i_stop,:);

        Uf=U(i_stop+1:i_end,:);       % f per futuro (per validazione) % qua tagliare la parte rumorosa finale.TO DO.
        yf=y(i_stop+1:i_end,:);

        dat_id = iddata(yp,Up,Ts);
        dat_val = iddata(yf,Uf,Ts);

    %%%%%%%%%%%% QUI SI IMPOSTANO p E n  %%%%%%%%%%%%%%%%%%%%%%%
        p=[25];
        n=[500];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for j=1:length(n)
            for i=1:length(p)

                cont_tmp = cont_tmp+1;
                strt = [' sseh temperatura ' num2str(cont_tmp) ' di ' num2str(len*length(n)*length(p))];
                disp(strt)
                fid = fopen('stato1.txt', 'w' );
                fprintf(fid, strt);
                fclose(fid);
                
                
                %%%%%%%%%%% yy
                tic 
                try 
                    
%                    [Mnpf_yy,iperf,Mnpe_yy,ipere] = stimanonparametrica_ric(yp,Up,p(i),n(j),'yy'); % chiamata principale.stima sseh.
                 index = zeros(1,24); % utilizzati per il punto 2
                   [Mnpf_yy,iperf,Mnpf_toc,Mnpe_yy,ipere,Mnpe_toc] = stimanonparametrica_mod(yp,Up,p(i),n(j),'yy',index); % chiamata principale.stima sseh.
                    
                catch exception
                Mnpf_yy = NaN; 
                iperf = NaN;
                Mnpe_yy = NaN;
                ipere = NaN;
                disp(['ite: sens ' num2str(k) ' n: ' num2str(j) 'p:' num2str(i) '. ' exception.message])                   
                end
                toc
                
                Dati_id(1,k).id_t_sseh_yy(i,j).Mnpf = Mnpf_yy;   
                Dati_id(1,k).id_t_sseh_yy(i,j).Mnpe = Mnpe_yy;
                Dati_id(1,k).id_t_sseh_yy(i,j).iperf = iperf;
                Dati_id(1,k).id_t_sseh_yy(i,j).ipere = ipere;
                Dati_id(1,k).id_t_sseh_yy(i,j).p_ = p(i);
                Dati_id(1,k).id_t_sseh_yy(i,j).n_ = n(j);
                Dati_id(1,k).id_t_sseh_yy(i,j).dati_val = dat_val;
                

                
               %%%%%%%%%%% ny
                tic
                try                   
                    [Mnpf_ny,iperf,Mnpe_ny,ipere] = stimanonparametrica_ric(yp,Up,p(i),n(j),'ny'); % chiamata principale.stima sseh.
                
                catch exception
                Mnpf_ny = NaN; 
                iperf = NaN;
                Mnpe_ny = NaN;
                ipere = NaN;
                disp(['ite: sens ' num2str(k) ' n: ' num2str(j) 'p:' num2str(i) '. ' exception.message])                   
                end
                toc
                
                    
                Dati_id(1,k).id_t_sseh_ny(i,j).Mnpf = Mnpf_ny; 
                Dati_id(1,k).id_t_sseh_ny(i,j).Mnpe = Mnpe_ny; 
                Dati_id(1,k).id_t_sseh_ny(i,j).iperf = iperf;
                Dati_id(1,k).id_t_sseh_ny(i,j).ipere = ipere;
                Dati_id(1,k).id_t_sseh_ny(i,j).p_ = p(i);
                Dati_id(1,k).id_t_sseh_ny(i,j).n_ = n(j);
                Dati_id(1,k).id_t_sseh_ny(i,j).dati_val = dat_val;
                
                
                

            end
        end

        Dati_id(1,k).aName = ['sensore numero ' num2str(k)];

        display(' ')
    end
    

    %%%%%%%%%%%%%%%%% AGGIUNTA SSEH TEMPERATURA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % aggiunge i campi dell' identificazione eseguita alla struttura Dati
    % si aggiunge l'ident. sseh sulla temperatura con modello 'yy' e 'ny'
    if NOR == 1 
        for i=2:2%1:len
            Dati(i).t.sseh_yy_norm = Dati_id(i).id_t_sseh_yy;
            %Dati(i).t.sseh_ny_norm = Dati_id(i).id_t_sseh_ny;
        end
    elseif NOR == 2
        for i=1:len
            Dati(i).t.sseh_yy = Dati_id(i).id_t_sseh_yy;
            Dati(i).t.sseh_ny = Dati_id(i).id_t_sseh_ny;
        end
    end

    save(str,'Dati')
    
  
end


%% SSEH DI UMIDITÀ

if (HUM && SSEH)
    cont_tmp = 0;
    for k=1:len 

        y = Dati(k).h.y_id; 
        U = Dati(k).h.Udati_id;   

        i_start= 50;  % 800 per ident, 600 per validazione. si scartano gli ultimi.
        i_stop = 850;
        
        Up=U(i_start:i_stop,:);     % p per passato
        yp=y(i_start:i_stop,:);

        Uf=U(i_stop+1:i_end,:);       % f per futuro (per validazione) % qua tagliare la parte rumorosa finale.TO DO.
        yf=y(i_stop+1:i_end,:);

        dat_id = iddata(yp,Up,Ts);
        dat_val = iddata(yf,Uf,Ts);

    %%%%%%%%%%%% QUI SI IMPOSTANO p E n  %%%%%%%%%%%%%%%%%%%%%%%
%         p=[ 5 10 20 40];
%         n=[ 10 50 100];
        p=[20];
        n=[200 300];
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


        for j=1:length(n)
            for i=1:length(p)

                cont_tmp = cont_tmp+1;
                strt = [' sseh umidita ' num2str(cont_tmp) ' di ' num2str(len*length(n)*length(p))];
                disp(strt)
                fid = fopen('stato1.txt', 'w' );
                fprintf(fid, strt);
                fclose(fid);

                [Mnpf_yy,iperf,Mnpe,ipere] = stimanonparametrica_ric(yp,Up,p(i),n(j),'yy'); % chiamata principale.stima sseh.
                Dati_id(1,k).id_h_sseh_yy(i,j).Mnpf = Mnpf_yy;   
                Dati_id(1,k).id_h_sseh_yy(i,j).p_ = p(i);
                Dati_id(1,k).id_h_sseh_yy(i,j).n_ = n(j);
                Dati_id(1,k).id_h_sseh_yy(i,j).dati_val = dat_val;

                [Mnpf_ny,iperf,Mnpe,ipere] = stimanonparametrica_ric(yp,Up,p(i),n(j),'ny'); % chiamata principale.stima sseh.
                Dati_id(1,k).id_h_sseh_ny(i,j).Mnpf = Mnpf_ny;   
                Dati_id(1,k).id_h_sseh_ny(i,j).p_ = p(i);
                Dati_id(1,k).id_h_sseh_ny(i,j).n_ = n(j);
                Dati_id(1,k).id_h_sseh_ny(i,j).dati_val = dat_val;

                %cont = cont+1;
            end
        end
    %     Dati_id(1,k).id_t_sseh_yy = id_t_sseh_yy;
    %     Dati_id(1,k).id_t_sseh_ny = id_t_sseh_ny;
        Dati_id(1,k).aName = ['sensore numero ' num2str(k)];

        display(' ')
    end
    

    %%%%%%%%%% AGGIUNTA SSEH UMIDITA' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % aggiunge i campi dei risultati delle nuove identificazioni alla struttura Dati
    % alla fine non modifica Dati_Matrice.mat ma salva l'unione in una nuova
    % struttura chiamata Dati_Matrice_id.mat

    % si aggiunge l'ident. sseh sulla temperatura con modello 'yy' e 'ny'
    if NOR == 1
        for i=1:len
            Dati(i).h.sseh_yy_norm = Dati_id(i).id_h_sseh_yy;
            Dati(i).h.sseh_ny_norm = Dati_id(i).id_h_sseh_ny;
        end
    else
        for i=1:len
            Dati(i).h.sseh_yy = Dati_id(i).id_h_sseh_yy;
            Dati(i).h.sseh_ny = Dati_id(i).id_h_sseh_ny;
        end
    end

    save(str,'Dati')

end


%% N4SID TEMPERATURA

if (TEM && N4SID)

    fit_m =[]; % salvo tutti i fit accetabili
    cont_tmp = 0;
    for k=1:len % per ogni sensore

        i_start = 10;  % si eliminano i primi dati

        y = Dati(k).t.y_id; 
        U = Dati(k).t.Udati_id;   

        %%%%%%%%%%%% ora si cicla al variare della finestra di dati usati per l'ident.
%         v = [500 600 700 800]; % <- qui si impostano le lunghezze dei dati per l'ident.
        v = [100 200 300 500 600 800];
        %%%%%%%%%%%%%%%
        
        
        for i=1:length(v) %[10 50 100 200 300 400]
            
            cont_tmp = cont_tmp+1;
            strt = [' n4sid temperatura ' num2str(cont_tmp) ' di ' num2str(len*length(v))];
            disp(strt)
            fid = fopen('stato1.txt', 'w' );
            fprintf(fid, strt);
            fclose(fid);
            
            i_stop = i_start+v(i);

            Up=U(i_start:i_stop,:);     % p per passato
            yp=y(i_start:i_stop,:);

            Uf=U(i_stop+1:i_end,:);       % f per futuro (per validazione) % qua tagliare la parte rumorosa finale.TO DO.
            yf=y(i_stop+1:i_end,:);

            dat_id = iddata(yp,Up,Ts);
            dat_val = iddata(yf,Uf,Ts);


            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%% N4SID
            % MODEL = N4SID(DATA) or MODEL = N4SID(DATA,ORDER)
            
            try            
                mod_n4sid = n4sid(dat_id,'best');
          
            catch exception
             mod_n4sid = NaN;                    
            disp(['ite: sens ' num2str(k) ' n: ' num2str(v(i)) '. ' exception.message])                   
            end
            
            
   

            %%%%%%%%% SALVATAGGIO DATI IN STRUTTURA
            lung_fin = i_stop-i_start;
           
            id_n4sid(k).struc_n4sid(i).dati_id = dat_id;
            id_n4sid(k).struc_n4sid(i).dati_val = dat_val;
            id_n4sid(k).struc_n4sid(i).iden_n4sid = mod_n4sid;
%             id_n4sid(k).struc_n4sid(i).fit = fit;
            id_n4sid(k).struc_n4sid(i).n = v(i);
           
            % trova il modello con il fit maggiore (accettabile) e salva in struttura.
%             vv = [];
%             for j=1:length(id_n4sid(k).struc_n4sid)
%                 vv = [vv id_n4sid(k).struc_n4sid(1,j).fit];
%             end

%             [fit_max i_max] = max(vv);

%             mod_fit_max = id_n4sid(k).struc_n4sid(1,i_max);     % modello con il fit maggiore 
            %display(['la finestra di ' num2str(v(i_max)) ' ha il fit max di ' num2str(fit_max)])

            % SALVO IN STRUTTURA il modello migliore (in base al fit) per quel sensore
%             id_n4sid(k).best_mod = mod_fit_max;

        end
    end
  
    %%%%%%%%%% AGGIUNTA N4SID TEMPERATURA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % aggiunge i campi dei risultati delle nuove identificazioni alla struttura Dati
    % alla fine non modifica Dati_Matrice.mat ma salva l'unione in una nuova
    % struttura chiamata Dati_Matrice_id.mat

    if NOR == 1
        for i=1:len
            Dati(i).t.n4sid_norm = id_n4sid(i);
        end
    else
        for i=1:len
            Dati(i).t.n4sid = id_n4sid(i);
        end
    end

    save(str,'Dati')
  
end


%% N4SID UMIDITA'
if (HUM && N4SID)

    fit_m =[]; % salvo tutti i fit accetabili
    
    cont_tmp = 0;
    for k=1:len % per ogni sensore

        i_start= 50;  % si eliminano i primi dati
        %i_stop = 850; % ora ha il significato di finestra temporale1

        
        y = Dati(k).h.y_id; 
        U = Dati(k).h.Udati_id;   

        %%%%%%%%%%%% ora si cicla al variare della finestra di dati usati per l'ident.
%         v = [500 600 700 800]; % <- qui si impostano le lunghezze dei dati per l'ident.
        v = [600 800 1000];
        %%%%%%%%%%%%%%%
        for i=1:length(v) %[10 50 100 200 300 400]
                  
            cont_tmp = cont_tmp+1;
            strt = [' n4sid umidita ' num2str(cont_tmp) ' di ' num2str(len*length(v))];
            disp(strt)
            fid = fopen('stato1.txt', 'w' );
            fprintf(fid, strt);
            fclose(fid);
            
            
            i_stop = i_start+v(i);

            Up=U(i_start:i_stop,:);     % p per passato
            yp=y(i_start:i_stop,:);

            Uf=U(i_stop+1:i_end,:);       % f per futuro (per validazione) % qua tagliare la parte rumorosa finale.TO DO.
            yf=y(i_stop+1:i_end,:);

            dat_id = iddata(yp,Up,Ts);
            dat_val = iddata(yf,Uf,Ts);


            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%% N4SID

            mod_n4sid = n4sid(dat_id,'best');
            
            %mod_n4sid = PBSID_opt_ARX_weighted()
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
            %%% calcolo del fit
            y_stim=sim(mod_n4sid,Uf);         % generale. sim il mod id con gli u futuri
            err_qm=mean((yf-y_stim).^2);      % errore quadratico medio 
            fit = 100*(1 - norm(yf-y_stim)/norm(yf-mean(yf)));
            p_tmp = length(yp);
            dat_val = iddata(yf,Uf,Ts);
            %predict (mod_n4sid,dat_val,1)

            %%%%%%%%% SALVATAGGIO DATI IN STRUTTURA
            lung_fin = i_stop-i_start;

            id_n4sid(k).struc_n4sid(i).dati_id = dat_id;
            id_n4sid(k).struc_n4sid(i).dati_val = dat_val;
            id_n4sid(k).struc_n4sid(i).iden_n4sid = mod_n4sid;
            id_n4sid(k).struc_n4sid(i).n = v(i);
            id_n4sid(k).struc_n4sid(i).fit = fit;
            %id_n4sid(k).num_sens = ['numero sensore' num2str(k)];
           
            % trova il modello con il fit maggiore (accettabile) e salva in struttura.
            vv = [];
            for j=1:length(id_n4sid(k).struc_n4sid)
                vv = [vv id_n4sid(k).struc_n4sid(1,j).fit];
            end

            [fit_max i_max] = max(vv);

            mod_fit_max = id_n4sid(k).struc_n4sid(1,i_max);     % modello con il fit maggiore 
            %display(['la finestra di ' num2str(v(i_max)) ' ha il fit max di ' num2str(fit_max)])

            % SALVO IN STRUTTURA il modello migliore (in base al fit) per quel sensore
            id_n4sid(k).best_mod = mod_fit_max;

        end
    end
  
    %%%%%%%%%% AGGIUNTA N4SID UMIDITA' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % aggiunge i campi dei risultati delle nuove identificazioni alla struttura Dati
    % alla fine non modifica Dati_Matrice.mat ma salva l'unione in una nuova
    % struttura chiamata Dati_Matrice_id.mat

    if NOR == 1
        for i=1:len
            Dati(i).h.n4sid_norm = id_n4sid(i);
        end
    else
        for i=1:len
            Dati(i).h.n4sid = id_n4sid(i);
        end
    end

    save(str,'Dati')

end


%% SSGLAR TEMPERATURA

if (TEM && SSGLAR)
    
    %%%%%%%%% PARAMETRI DA CAMBIARE:  
    pRKHS_vet=[5 25]; % length of eigenfunctions
    n_vet=[100 300 500];%input('# of measurements = ');%numero dati effettivi per la stima (50,300)
    %%%%%%%%%%
  
    frac=1/3;% fraction of data used for validation/selection
    cont_tmp = 0;
    for num_sens=1:len; %disp(['num_sens ' num2str(num_sens)])

        U = Dati(num_sens).t.Udati_id;
        y = Dati(num_sens).t.y_id;

        for num_n=1:length(n_vet)
            n=n_vet(num_n);     

            ntval=floor(frac*(n));%dimensione validation set

            for num_pRKHS=1:length(pRKHS_vet)
                
                cont_tmp = cont_tmp+1;
                strt = [' ssglar temperatura ' num2str(cont_tmp) ' di ' num2str(len*length(n_vet)*length(pRKHS_vet))];
                disp(strt)
                fid = fopen('stato1.txt', 'w' );
                fprintf(fid, strt);
                fclose(fid);
                
                pRKHS=pRKHS_vet(num_pRKHS);  

                ya=y(1:(n-ntval+pRKHS),:);
                Ua=U(1:(n-ntval+pRKHS),:);
                yval=y((n-ntval+1+pRKHS):n,:);
                Uval=U((n-ntval+1+pRKHS):n,:);

                clear AIK 
                tic
                try

                    %%%%%%%%%%%% STIMA NONPARAMETRICA 2%%%%%%%%%%%%%%%%%%%%%%%
                    [Mnpg,beta,sigma,NT,Mnpc,beta2,sigma2,NT2] = RKHS_Pred_LARS_Validation_all_k_v2(ya,Ua,1,30,yval,Uval,pRKHS);
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                catch exception                   
                    Mnpg = NaN;
                    beta = NaN;
                    sigma = NaN;
                    NT = NaN;
                    Mnpc = NaN;
                    beta2 = NaN;
                    sigma2 = NaN;
                    NT2 = NaN;
                    disp(['ite: sens ' num2str(num_sens) ' n: ' num2str(num_n) ' p ' num2str(pRKHS) '. ' exception.message])
                    
                end
                toc
                id_ssglar(1,num_sens).mod(num_pRKHS,num_n) = struct('Mnpg',Mnpg,'beta',beta,'sigma',sigma,'NT',NT,'Mnpc',Mnpc,'beta2',beta2,'sigma2',sigma2,'NT2',NT2,'n',n,'pRKHS',pRKHS);

            end       
        end

    % save('mod_RKHS_stand','mod_RKHS');
    end

    %%%%%%%%%% AGGIUNTA SSGLAR TEMPERATURA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if NOR == 1
        for i=1:len
            Dati(i).t.ssglar_norm = id_ssglar(i);
        end
    else
        for i=1:len
            Dati(i).t.ssglar = id_ssglar(i);
        end
    end

    save(str,'Dati')

end


%% SSGLAR UMIDITA'
if (HUM && SSGLAR)
    
    %%%%%%%%%%%%%% PARAMETRI DA CAMBIARE:   
    pRKHS_vet=[30]; % length of eigenfunctions
    n_vet=[500];%input('# of measurements = ');%numero dati effettivi per la stima (50,300)
    %%%%%%%%%%%%%%%%%%
  
    frac=1/3;% fraction of data used for validation/selection
    cont_tmp = 0;
    for num_sens=1:len; %disp(['num_sens ' num2str(num_sens)])

        U = Dati(num_sens).h.Udati_id;
        y = Dati(num_sens).h.y_id;

        for num_n=1:length(n_vet)
            n=n_vet(num_n);

            ntval=floor(frac*(n));%dimensione validation set

            for num_pRKHS=1:length(pRKHS_vet)
                
                cont_tmp = cont_tmp+1;
                strt = [' ssglar umidita ' num2str(cont_tmp) ' di ' num2str(len*length(n_vet)*length(pRKHS_vet))];
                disp(strt)
                fid = fopen('stato1.txt', 'w' );
                fprintf(fid, strt);
                fclose(fid);
                
                pRKHS=pRKHS_vet(num_pRKHS);  
                
                ya=y(1:(n-ntval+pRKHS),:);
                Ua=U(1:(n-ntval+pRKHS),:);
                yval=y((n-ntval+1+pRKHS):n,:);
                Uval=U((n-ntval+1+pRKHS):n,:);

                clear AIK  
                try

                    %%%%%%%%%%%% STIMA NONPARAMETRICA 2%%%%%%%%%%%%%%%%%%%%%%%
                    [Mnpg,beta,sigma,NT,Mnpc,beta2,sigma2,NT2] = RKHS_Pred_LARS_Validation_all_k_v2(ya,Ua,1,30,yval,Uval,pRKHS);
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                catch exception
                    
                    Mnpg = NaN;
                    beta = NaN;
                    sigma = NaN;
                    NT = NaN;
                    Mnpc = NaN;
                    beta2 = NaN;
                    sigma2 = NaN;
                    NT2 = NaN; 
                    disp(['ite: sens ' num2str(num_sens) ' n: ' num2str(num_n) ' p ' num2str(pRKHS) '. ' exception.message])
                    
                end
                
                id_ssglar(1,num_sens).mod(num_pRKHS,num_n) = struct('Mnpg',Mnpg,'beta',beta,'sigma',sigma,'NT',NT,'Mnpc',Mnpc,'beta2',beta2,'sigma2',sigma2,'NT2',NT2,'n',n,'pRKHS',pRKHS);

            end       
        end

    % save('mod_RKHS_stand','mod_RKHS');
    end


    %%%%%%%%%% AGGIUNTA SSGLAR UMIDITA' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if NOR == 1
        for i=1:len
            Dati(i).h.ssglar_norm = id_ssglar(i);
        end
    else
        for i=1:len
            Dati(i).h.ssglar = id_ssglar(i);
        end
    end

    save(str,'Dati')

end


%% GLAR TEMPERATURA

if (TEM && GLAR)
    %%%%%% parametri variabili
    n_vet=[100 300 500 ];%input('# of measurements = ');%numero dati effettivi per la stima (50,300)
    maxp_GLARS = 30; % maximum order for GLARS;
    frac=1/3;% fraction of data used for validation/selection
    %%%%%%%%%%%%

    %creazione matrice ingressi U, vettore uscite y
    cont_tmp = 0;
    for num_sens=1:len

        for num_n=1:length(n_vet)
            
                 
            cont_tmp = cont_tmp+1;
            strt = [' glar temperatura ' num2str(cont_tmp) ' di ' num2str(len*length(n_vet))];
            disp(strt)
            fid = fopen('stato1.txt', 'w' );
            fprintf(fid, strt);
            fclose(fid);
            
            n=n_vet(num_n);

            ntval=floor(frac*(n));%dimensione validation set

            U = Dati(num_sens).t.Udati_id;
            y = Dati(num_sens).t.y_id;

            yaLARS=y(1:(n-ntval),:);
            UaLARS=U(1:(n-ntval),:);
            yvalLARS=y((n-ntval+1):n,:);
            UvalLARS=U((n-ntval+1):n,:);

            clear AIK    
            tic
            try
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%% STIMA GLARS VALIDATION %%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                dat_val=iddata(yvalLARS,UvalLARS);
                t=0;
                for pl=1:1:maxp_GLARS;
                    t=t+1;
                    ll(t,:)=pl;
                    [Mgl,Mgls] = Dyn_GLARS_all_k(yaLARS,[],UaLARS,pl,1,0,yvalLARS,[],UvalLARS);
                    predd=predict(Mgl,dat_val);predd=iddata(predd);predd=predd.y;
                    RM1(t)=sum((predd-yvalLARS).^2);
                    predd=predict(Mgls,dat_val);predd=iddata(predd);predd=predd.y;
                    RM2(t)=sum((predd-yvalLARS).^2);
                end
                [q,w1]=min(RM1);ll(w1,:)
                [Mgl,Mgls] = Dyn_GLARS_all_k(yaLARS,[],UaLARS,ll(w1,:),1,0,yvalLARS,[],UvalLARS);

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            
            catch exception
                 w1=1;
                 Mgl = NaN;
                 Mgls = NaN;
                 ll(w1,:) = NaN;
                 disp(['ite: sens ' num2str(num_sens) ' n: ' num2str(num_n) '. ' exception.message])
                    
            end
            toc
            id_glar(1,num_sens).mod(1,num_n) = struct('Mgl',Mgl,'Mgls',Mgls,'deg',ll(w1,:),'n',n);

        end
    end
    
        %%%%%%%%%% AGGIUNTA GLAR TEMPERATURA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if NOR == 1
        for i=1:len
            Dati(i).t.glar_norm = id_glar(i);
        end
    else
        for i=1:len
            Dati(i).t.glar = id_glar(i);
        end
    end

    save(str,'Dati')
    
    
end


%% GLAR UMIDITA'

if (HUM && GLAR)
    %%%%%%% parametri variabili
    n_vet=[ 200 500 ];%input('# of measurements = ');%numero dati effettivi per la stima (50,300)
    maxp_GLARS = 30; % maximum order for GLARS;
    frac=1/3;% fraction of data used for validation/selection
    %%%%%%%%%%%%%%%%

    %creazione matrice ingressi U, vettore uscite y
    cont_tmp = 0;
    for num_sens=1:len

        for num_n=1:length(n_vet)
            
            cont_tmp = cont_tmp+1;
            strt = [' glar umidita ' num2str(cont_tmp) ' di ' num2str(len*length(n_vet))];
            disp(strt)
            fid = fopen('stato1.txt', 'w' );
            fprintf(fid, strt);
            fclose(fid);            
            
            n=n_vet(num_n);

            ntval=floor(frac*(n));%dimensione validation set

            U = Dati(num_sens).h.Udati_id;
            y = Dati(num_sens).h.y_id;

            yaLARS=y(1:(n-ntval),:);
            UaLARS=U(1:(n-ntval),:);
            yvalLARS=y((n-ntval+1):n,:);
            UvalLARS=U((n-ntval+1):n,:);

            clear AIK 
            try
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%% STIMA GLARS VALIDATION %%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                dat_val=iddata(yvalLARS,UvalLARS);
                t=0;
                for pl=1:1:maxp_GLARS;
                    t=t+1;
                    ll(t,:)=pl;
                    [Mgl,Mgls] = Dyn_GLARS_all_k(yaLARS,[],UaLARS,pl,1,0,yvalLARS,[],UvalLARS);
                    predd=predict(Mgl,dat_val);predd=iddata(predd);predd=predd.y;
                    RM1(t)=sum((predd-yvalLARS).^2);
                    predd=predict(Mgls,dat_val);predd=iddata(predd);predd=predd.y;
                    RM2(t)=sum((predd-yvalLARS).^2);
                end
                [q,w1]=min(RM1);ll(w1,:)
                [Mgl,Mgls] = Dyn_GLARS_all_k(yaLARS,[],UaLARS,ll(w1,:),1,0,yvalLARS,[],UvalLARS);

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            catch exception
                 Mgl = NaN;
                 Mgls = NaN;
                 ll(w1,:) = NaN;
                 disp(['ite: sens ' num2str(num_sens) ' n: ' num2str(num_n) '. ' exception.message])                   
            end

            id_glar(1,num_sens).mod(1,num_n) = struct('Mgl',Mgl,'Mgls',Mgls,'deg',ll(w1,:),'n',n); 
        end
    end
    
        %%%%%%%%%% AGGIUNTA GLAR UMIDITA' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if NOR == 1
        for i=1:len
            Dati(i).h.glar_norm = id_glar(i);
        end
    else
        for i=1:len
            Dati(i).h.glar = id_glar(i);
        end
    end

    save(str,'Dati')


   
end


%% PEM TEMPERATURA 

if (TEM && PEM)
    
    %%%%%%%% parametri variabili:
    n_vet=[100 300 500 700 900 1500];
    %%%%%%%%%%%%
    
    cont_tmp = 0;

    for num_sens=1:len; 
        for num_n=1:length(n_vet)

        cont_tmp = cont_tmp+1;
        strt = [' pem temperatura ' num2str(cont_tmp) ' di ' num2str(len*length(n_vet))];
        disp(strt)
        fid = fopen('stato1.txt', 'w' );
        fprintf(fid, strt);
        fclose(fid);
           
        n = n_vet(num_n);
         
        y = Dati(num_sens).t.y_id(1:n,:);
        U = Dati(num_sens).t.Udati_id(1:n,:);

        m=size(U,2);
        mPEM=15;%20;%max ordine PEM
        
        tic
        try
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%% STIMA PEM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            clear L1 L2 ll1 Or0 Or1 Or2 RROrM RROr2M DOr SOr Orr0 Orr1 Orr2
            dat=iddata(y,U);
            t=0;
            for i1=1:mPEM;            
                    i2=i1*ones(1,m);
                    i3=i1;
                    t=t+1;
                    Mpem=pem(dat,[i1 i2 i3 zeros(1,m+1) ones(1,m)],'LimitError',0);
                    predd=predict(Mpem,dat);predd=iddata(predd);predd=predd.y;
                    RM=sum((predd-y).^2);npar=max(size(Mpem.ParameterVector));ndat=Mpem.EstimationInfo.DataLength;
                    L1(t)=ndat*log(RM/ndat)+log(ndat)*npar;%BIC
                    L2(t)=ndat*log(RM/ndat)+[ndat*(ndat+npar)]/[ndat-npar-2];%AICC
                    ll1(t,:)=[i1 i2 i3];

            end
            [q,w1]=min(L1);ll1(w1,:);
            Mpem_BIC=pem(dat,[ll1(w1,:) zeros(1,m+1) ones(1,m)],'LimitError',0);
            [q,w2]=min(L2);ll1(w2,:);
            Mpem_AICC=pem(dat,[ll1(w2,:) zeros(1,m+1) ones(1,m)],'LimitError',0);

        catch exception
            Mpem_BIC = NaN;
            Mpem_AICC = NaN;                    
            disp(['ite: sens ' num2str(num_sens) ' n: ' num2str(num_n) '. ' exception.message])                   
        end
        toc
            id_pem(1,num_sens).mod(1,num_n) = struct('Mpem_BIC',Mpem_BIC,'Mpem_AICC',Mpem_AICC,'n',n);

        end


    end

    %%%%%%%%%% AGGIUNTA PEM TEMPERATURA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if NOR == 1
        for i=1:len
            Dati(i).t.pem_norm = id_pem(i);
        end
    else
        for i=1:len
            Dati(i).t.pem = id_pem(i);
        end
    end

    save(str,'Dati')


    
end


%% PEM UMIDITA' 

if (HUM && PEM)
    
    %%%%%%% parametri varibili
    n_vet=[700 900 1300];
    %%%%%%%%%%%%%%%%%%%%
    
    cont_tmp = 0;
    for num_sens=1:len; 
        for num_n=1:length(n_vet)
            
        cont_tmp = cont_tmp+1;
        strt = [' pem umidita ' num2str(cont_tmp) ' di ' num2str(len*length(n_vet))];
        disp(strt)
        fid = fopen('stato1.txt', 'w' );
        fprintf(fid, strt);
        fclose(fid);
            
        n = n_vet(num_n);
         
        y = Dati(num_sens).h.y_id(1:n,:);
        U = Dati(num_sens).h.Udati_id(1:n,:);

        m=size(U,2);
        mPEM=15;%20;%max ordine PEM

        try
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%% STIMA PEM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            clear L1 L2 ll1 Or0 Or1 Or2 RROrM RROr2M DOr SOr Orr0 Orr1 Orr2
            dat=iddata(y,U);
            t=0;
            for i1=1:mPEM;                   
                    i2=i1*ones(1,m);
                    i3=i1;
                    t=t+1;
                    Mpem=pem(dat,[i1 i2 i3 zeros(1,m+1) ones(1,m)],'LimitError',0);
                    predd=predict(Mpem,dat);predd=iddata(predd);predd=predd.y;
                    RM=sum((predd-y).^2);npar=max(size(Mpem.ParameterVector));ndat=Mpem.EstimationInfo.DataLength;
                    L1(t)=ndat*log(RM/ndat)+log(ndat)*npar;%BIC
                    L2(t)=ndat*log(RM/ndat)+[ndat*(ndat+npar)]/[ndat-npar-2];%AICC
                    ll1(t,:)=[i1 i2 i3];

            end
            [q,w1]=min(L1);ll1(w1,:);
            Mpem_BIC=pem(dat,[ll1(w1,:) zeros(1,m+1) ones(1,m)],'LimitError',0);
            [q,w2]=min(L2);ll1(w2,:);
            Mpem_AICC=pem(dat,[ll1(w2,:) zeros(1,m+1) ones(1,m)],'LimitError',0);
            
        catch exception
            Mpem_BIC = NaN;
            Mpem_AICC = NaN;                    
            disp(['ite: sens ' num2str(num_sens) ' n: ' num2str(num_n) '. ' exception.message])                   
        end

            id_pem(1,num_sens).mod(1,num_n) = struct('Mpem_BIC',Mpem_BIC,'Mpem_AICC',Mpem_AICC,'n',n);

        end


    end

    %%%%%%%%%% AGGIUNTA PEM UMIDITA' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if NOR == 1
        for i=1:len
            Dati(i).h.pem_norm = id_pem(i);
        end
    else
        for i=1:len
            Dati(i).h.pem = id_pem(i);
        end
    end

    save(str,'Dati')

    
end

