
% usa la funzione ssspline2 in cui si possono vincolare degli ingressi
% attraverso il vettore index (1 i-esima componente -> i-esimo ingresso vincolato)
% id della casa di chiuso


clc
clear


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


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


iniz % qui si effettua la normalizzazione ecc.



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

if (TEM && SSEH)
	cont_tmp = 0;
    
    
    sens_canc = [7 9 10 14 15 17 18 19 21 22 24];     % selezionare i sensori da eliminare
    sens_selez = [1 2 3 4 5 6 8 11 12 13 16 20 23];
    sens_selez_num = [1 2 3 4 5 6 8 3 4 5 8 4 7]; % per asse x grafico excel

        y_id = Dati(1).t.y_id; 
        U_id = Dati(1).t.Udati_id;
        U_id = [y_id U_id];
        
    for k=1:len    % SCELTA SENSORI
    
        % qui chiamo solo la funzione di andrea stat
        U = cvcfu(U_id, k, sens_canc);
        y = Dati(k).t.y_id;
        
        % se k fa parte della selezione identificala usando un ny
        if (length(find(sens_canc==k)) == 0) % se non lo trova restituisce empty matrix
            ny_yy = 'yy'; 
        else
            ny_yy = 'ny'; 
        end


        i_start= 10;  
        i_stop = 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)
        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 e index %%%%%%%%%%%%%%%%%%%%%%%
        p=[25];
        n=[500];
        
        ind_term = find (Dati(k).t.Colonne_U == 5);
        ind_est = find (Dati(k).t.Colonne_U == 8);
        
        index = zeros(1,23);  % per vincolare delle variabili
% commentare per non effettuare vincoli       
%         if ind_term(1) < 8 
%             index(1,ind_term(1)) = 1;
%         end
%         if ind_est(1) < 8
%             index(1,ind_est(1)) = 1;
%         end
        
        index = [0 index]; % aggiunta posiz y
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        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('stato.txt', 'w' );
                fprintf(fid, strt);
                fclose(fid);
                
         
                
                %%%% stima modello ny_yy
                tic 
                try 
                    index=zeros(1,13);
                    [Mnpf,iperf,~,Mnpe,ipere,~] = stimanonparametrica_mod(yp,Up,p(i),n(j),ny_yy,index); % chiamata principale.stima sseh.
                catch exception
                Mnpf = NaN;
                Mnpf_toc = NaN;
                iperf = NaN;
                Mnpe = NaN;
                ipere = NaN;
                nve=NaN;
                nvf=NaN;
                disp(['ite: sens ' num2str(k) ' n: ' num2str(j) 'p:' num2str(i) '. ' exception.message])                   
                end
                                 
                
                Dati_id(1,k).id_t_sseh(i,j).Mnpf = Mnpf; 
                Dati_id(1,k).id_t_sseh(i,j).Mnpe = Mnpe;
                Dati_id(1,k).id_t_sseh(i,j).ny_yy = ny_yy;
                Dati_id(1,k).id_t_sseh(i,j).sens_selez = sens_selez;
                Dati_id(1,k).id_t_sseh(i,j).p_ = p(i);
                Dati_id(1,k).id_t_sseh(i,j).n_ = n(j);
                Dati_id(1,k).id_t_sseh(i,j).dat_val = dat_val;
                Dati_id(1,k).id_t_sseh(i,j).index=index;
                Dati_id(1,k).id_t_sseh(i,j).iperf = iperf;
                Dati_id(1,k).id_t_sseh(i,j).ipere = ipere;

                toc
               

            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'

        for i=1:len
            Dati(i).t.sseh_norm = Dati_id(i).id_t_sseh;
        end

    save(str,'Dati')
    

end





%% sseh_norm validazione  Mnpf

% attenzione in fid ai / o \


%%%%%%%%%%%%%%%% intervallo dati di validazione
val_start=i_stop+1;
val_stop=length(Dati(1,2).t.sseh_norm.dat_val.y);

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


fid = fopen('ValidazioneChiuso\sseh_norm_Mnpf_selez.txt', 'w+');

str_spar = 'Sparsit�:\t';
str_num_np = '';
str_nome_dati = '';
str_dati = '';
str_dati_sim = 'Spar_sim\t';

for num_sens=1:length(Dati) % scelta sensore 
    

    dat_val = Dati(1,num_sens).t.sseh_norm.dat_val;
    
    
    [len_p len_n] = size(Dati(1,num_sens).t.sseh_norm);
    
    str_nome_id = ['\tDati(1,' num2str(num_sens) ').t.sseh_norm(p,n).Mnpf'];
    fprintf(fid, [str_nome_id, '\n']);
    
    disp(['Sensore: Dati(1,' num2str(num_sens) ').t.sseh_norm(p,n).Mnpf']);
    
    for k=1:30 % passo di predizione
        

        str_dati = [str_dati, num2str(k), '\t'];
        
        for num_n=1:len_n
            for num_p=1:len_p
                
                try
                if k==1
                    
                    Model = Dati(1,num_sens).t.sseh_norm(num_p,num_n).Mnpf;
                    
                    spars = sparsi(Model);
                    str_spar = [str_spar num2str(spars) '%%\t\t\t'];
                    str_num_np = [str_num_np 'n=' num2str(Dati(1,num_sens).t.sseh_norm(num_p,num_n).n_) ' p=' num2str(Dati(1,num_sens).t.sseh_norm(num_p,num_n).p_) '\t\t\t'];
                    str_nome_dati = [str_nome_dati 'fit\trms\tcod\t'];
                    
                    y_sim = sim(Model,dat_val.u);
                    y_sim = iddata(y_sim);
                    y_sim = y_sim.y;

                    %Calcolo indici
                    fit = (1 - norm(yf-y_sim)/norm(yf-mean(yf)));
                    rms = sqrt(mean((yf-y_sim).^2));
                    cod = 1-(rms^2/(mean((yf-mean(yf)).^2)));

                    str_dati_sim = [str_dati_sim, num2str(fit) '\t', num2str(rms) '\t', num2str(cod) '\t'];
                    
                end
                
                y_pred = predict(Model,dat_val,k);
                y_pred = iddata(y_pred);
                y_pred = y_pred.y;
                
                %Calcolo indici
                fit = (1 - norm(yf-y_pred)/norm(yf-mean(yf)));
                rms = sqrt(mean((yf-y_pred).^2));
                cod = 1-(rms^2/(mean((yf-mean(yf)).^2)));
                
                
                catch exception
                
                fit = NaN;
                rms = NaN;
                cod = NaN;
                
                if k==1
                    str_spar = [str_spar num2str(NaN) '%%\t\t\t'];
                    str_num_np = [str_num_np 'n=' num2str(Dati(1,num_sens).t.sseh_norm(num_p,num_n).n_) ' p=' num2str(Dati(1,num_sens).t.sseh_norm(num_p,num_n).p_) '\t\t\t'];
                    str_nome_dati = [str_nome_dati 'fit\trms\tcod\t'];
                    str_dati_sim = [str_dati_sim, num2str(fit) '\t', num2str(rms) '\t', num2str(cod) '\t'];
                end
                                   
                end
                
                str_dati = [str_dati, num2str(fit) '\t', num2str(rms) '\t', num2str(cod) '\t'];
                
            end
        end
        str_dati = [str_dati, '\n'];
    end
    
    fprintf(fid, ['\t' str_num_np, '\n']);    
    fprintf(fid, [str_spar '\n']);
    fprintf(fid, ['\t' str_nome_dati, '\n']);
    fprintf(fid, [str_dati]);
    fprintf(fid, [str_dati_sim, '\n \n \n' ]);
    
    str_spar = 'Sparsit�:\t';
    str_num_np = '';
    str_nome_dati = '';
    str_dati = '';
    str_dati_sim = 'Spar_sim\t';
    
end%% sseh_norm

fclose(fid);


%% sseh_norm validazione  Mnpe

% attenzione in fid ai / o \


%%%%%%%%%%%%%%%% intervallo dati di validazione
val_start=i_stop+1;
val_stop=length(Dati(1,2).t.sseh_norm.dat_val.y);

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


fid = fopen('ValidazioneChiuso\sseh_norm_Mnpe_selez.txt', 'w+');

str_spar = 'Sparsit�:\t';
str_num_np = '';
str_nome_dati = '';
str_dati = '';
str_dati_sim = 'Spar_sim\t';

for num_sens=1:length(Dati) % scelta sensore 
    

    dat_val = Dati(1,num_sens).t.sseh_norm.dat_val;
    
    
    [len_p len_n] = size(Dati(1,num_sens).t.sseh_norm);
    
    str_nome_id = ['\tDati(1,' num2str(num_sens) ').t.sseh_norm(p,n).Mnpe'];
    fprintf(fid, [str_nome_id, '\n']);
    
    disp(['Sensore: Dati(1,' num2str(num_sens) ').t.sseh_norm(p,n).Mnpe']);
    
    for k=1:30 % passo di predizione
        

        str_dati = [str_dati, num2str(k), '\t'];
        
        for num_n=1:len_n
            for num_p=1:len_p
                
                try
                if k==1
                    
                    Model = Dati(1,num_sens).t.sseh_norm(num_p,num_n).Mnpe;
                    
                    spars = sparsi(Model);
                    str_spar = [str_spar num2str(spars) '%%\t\t\t'];
                    str_num_np = [str_num_np 'n=' num2str(Dati(1,num_sens).t.sseh_norm(num_p,num_n).n_) ' p=' num2str(Dati(1,num_sens).t.sseh_norm(num_p,num_n).p_) '\t\t\t'];
                    str_nome_dati = [str_nome_dati 'fit\trms\tcod\t'];
                    
                    y_sim = sim(Model,dat_val.u);
                    y_sim = iddata(y_sim);
                    y_sim = y_sim.y;

                    %Calcolo indici
                    fit = (1 - norm(yf-y_sim)/norm(yf-mean(yf)));
                    rms = sqrt(mean((yf-y_sim).^2));
                    cod = 1-(rms^2/(mean((yf-mean(yf)).^2)));

                    str_dati_sim = [str_dati_sim, num2str(fit) '\t', num2str(rms) '\t', num2str(cod) '\t'];
                    
                end
                
                y_pred = predict(Model,dat_val,k);
                y_pred = iddata(y_pred);
                y_pred = y_pred.y;
                
                %Calcolo indici
                fit = (1 - norm(yf-y_pred)/norm(yf-mean(yf)));
                rms = sqrt(mean((yf-y_pred).^2));
                cod = 1-(rms^2/(mean((yf-mean(yf)).^2)));
                
                
                catch exception
                
                fit = NaN;
                rms = NaN;
                cod = NaN;
                
                if k==1
                    str_spar = [str_spar num2str(NaN) '%%\t\t\t'];
                    str_num_np = [str_num_np 'n=' num2str(Dati(1,num_sens).t.sseh_norm(num_p,num_n).n_) ' p=' num2str(Dati(1,num_sens).t.sseh_norm(num_p,num_n).p_) '\t\t\t'];
                    str_nome_dati = [str_nome_dati 'fit\trms\tcod\t'];
                    str_dati_sim = [str_dati_sim, num2str(fit) '\t', num2str(rms) '\t', num2str(cod) '\t'];
                end
                                   
                end
                
                str_dati = [str_dati, num2str(fit) '\t', num2str(rms) '\t', num2str(cod) '\t'];
                
            end
        end
        str_dati = [str_dati, '\n'];
    end
    
    fprintf(fid, ['\t' str_num_np, '\n']);    
    fprintf(fid, [str_spar '\n']);
    fprintf(fid, ['\t' str_nome_dati, '\n']);
    fprintf(fid, [str_dati]);
    fprintf(fid, [str_dati_sim, '\n \n \n' ]);
    
    str_spar = 'Sparsit�:\t';
    str_num_np = '';
    str_nome_dati = '';
    str_dati = '';
    str_dati_sim = 'Spar_sim\t';
    
end%% sseh_norm

fclose(fid);





  %% norma degli ingressi selezionati






%%%%%%%%%%% conf per Mnpf

% stampiamo per ogni sensore (la norma di) gli ingressi usati (per gli Mnpf)
for k=1:3%1:length(Dati)    
    in_sel_tmp = [];
    tt = length(Dati(1,k).t.sseh_norm.Mnpf.nb);
    for i=1:tt
        norm_i = norm(Dati(1,k).t.sseh_norm.Mnpf(1,i).b);
        in_sel_tmp = [in_sel_tmp norm_i];
    end
    
    Conf_ingressi.in_sseh(1,k).in_selez = in_sel_tmp;
end

% stampa in txt

fid = fopen('ValidazioneChiuso\conf_norme_ingressi_sseh_selez_Mnpf.txt', 'w+');

caso = ['Mnpf sseh' '\t'];
fprintf(fid, [ caso, '\n']);    
for k=2:2%1:length(Dati)
    v = ['Sens \t'];
    tt = length(Dati(1,k).t.sseh_norm.Mnpf.nb);
    for i=1:tt
       v = [v num2str(Dati(k).t.sseh_norm.sens_selez(i)) '\t' ];  
    end
    fprintf(fid, [v  '\n']);
    v=[];
    for i=1:tt
       v = [v num2str(Conf_ingressi.in_sseh(1,k).in_selez(i),3) '\t' ]; 
    end
    fprintf(fid, [num2str(k) '\t' v  '\n']);
end 


fclose(fid);


disp('fine')

