
% Pem su ingressi selezionati uguali per tutti, quindi U ridotta

clc
clear


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


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


iniz % qui si effettua la normalizzazione ecc.



%% PEM DI TEMPERATURA %%%%%%%%%%%%%%%%%%%%%%%%%%%

if (TEM && PEM)
	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 num_sens=1:len    % SCELTA SENSORI
    
        
        % riduzione della U
        U_sel = cvcfu(U_id, num_sens, sens_canc);
        y_sel = Dati(num_sens).t.y_id;
        
       
        % se num_sens fa parte della selezione identificala usando un ny
        if (length(find(sens_canc==num_sens)) == 0) % se non lo trova restituisce empty matrix
            ny_yy = 'yy'; 
        else
            ny_yy = 'ny'; 
        end

        % intervallo dati per l'identificazione:
        i_start= 10;  
        i_stop = 1520;
        
        
        Up=U_sel(i_start:i_stop,:);     % p per passato
        yp=y_sel(i_start:i_stop,:);

        Uf=U_sel(i_stop+1:i_end,:);       % f per futuro (per validazione)
        yf=y_sel(i_stop+1:i_end,:);

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

    %%%%%%%%%%%% QUI SI IMPOSTANO  n  %%%%%%%%%%%%%%%%%%%%%%%
        n_vet=[1500];
        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for num_n=1:length(n_vet)
                 
       %%%% stima modello pem

        cont_tmp = cont_tmp+1;
        strt = [' pem temperatura ' num2str(cont_tmp) ' di ' num2str(len*length(n_vet))];
        disp(strt)
        fid = fopen('stato.txt', 'w' );
        fprintf(fid, strt);
        fclose(fid);
           
        n = n_vet(num_n);
        
            % dati per l'iden.
                     
            U = dat_id.u(1:n,:); % <-ora contiene solo le col selez
            y = dat_id.y(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_sel(1,num_sens).mod = struct('Mpem_BIC',Mpem_BIC,'Mpem_AICC',Mpem_AICC,'dat_val',dat_val,'n',n,'sens_selez',sens_selez);
                %qui sopra non si salva in funz di num_n.
        end
    


    end

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

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

    save(str,'Dati')

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


end




%% ANALISI PEM


%%%%%%%%%%%%%%%%%%%% BIC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid = fopen('ValidazioneChiuso/pem_selez_BIC1500_modUrid.txt', 'w+');

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

for num_sens=1:length(Dati) % scelta sensore 
    
    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'];
      
                try
                if k==1

                    str_nome_dati = [str_nome_dati 'fit\trms\tcod\t'];
                    
                    dat_val = Dati(1,num_sens).t.pem_norm.mod.dat_val;
                    yf = dat_val.y;
                    uf = dat_val.u;
                    
                    Model = Dati(1,num_sens).t.pem_norm.mod.Mpem_BIC;
                    
                    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
                disp(['sens' num2str(num_sens) ': ' exception.message])
                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'];
                
            
        
        str_dati = [str_dati, '\n'];
    end
      
    fprintf(fid, ['\t' str_nome_dati, '\n']);
    fprintf(fid, [str_dati]);
    fprintf(fid, [str_dati_sim, '\n \n \n' ]);
    
    str_nome_dati = '';
    str_dati = '';
    str_dati_sim = 'Spar_sim\t';
    
    
   
    
end

fclose(fid);
 disp('fine')

