% ANALISI PEM

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

str_num_np = '';
str_nome_dati = '';
str_dati = '';
str_dati_sim = 'Spar_sim\t';

for num_sens=1:length(Dati) % scelta sensore 
    
    
    [len_p len_n] = size(Dati(1,num_sens).t.pem_norm.mod);
    
    str_nome_id = ['\tDati(1,' num2str(num_sens) ').t.pem_norm(p,n).Mpem_BIC'];
    fprintf(fid, [str_nome_id, '\n']);
    
    disp(['Sensore: Dati(1,' num2str(num_sens) ').t.pem_norm(p,n).Mpem_BIC']);
    
    for k=1:30 % passo di predizione
        

        str_dati = [str_dati, num2str(k), '\t'];
        
        for num_n=3:3 % solo con sseh che ha n=500
            for num_p=1:len_p
                
                try
                if k==1
                    
                    str_num_np = [str_num_np 'n_sel_sseh=' num2str(Dati(1,num_sens).t.pem_norm.mod(num_p,num_n).n_sel_sseh)  '\t\t\t'];
                    str_nome_dati = [str_nome_dati 'fit\trms\tcod\t'];
                    
                    dat_val = Dati(1,num_sens).t.pem_norm.mod(num_p,num_n).dat_val;
                    yf = dat_val.y;
                    uf = dat_val.u;
                    
                    Model = Dati(1,num_sens).t.pem_norm.mod(num_p,num_n).Mpem_BIC;
                     
                    spars = sparsi(Model);
                    
                    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_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, ['\t' str_nome_dati, '\n']);
    fprintf(fid, [str_dati]);
    fprintf(fid, [str_dati_sim, '\n \n \n' ]);
    
    str_num_np = '';
    str_nome_dati = '';
    str_dati = '';
    str_dati_sim = 'Spar_sim\t';
    
    
    
    
end

fclose(fid);