
% pem con ingressi selezionati da sseh. l'ultima parte è da far girare per
% la validazione.


clc
clear


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


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



%% 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 = true
% % SSGLAR = true
% % GLAR = true
PEM = true


% 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


%% PEM TEMPERATURA 

if (TEM && PEM)
    
    %%%%%%%% parametri variabili:
    n_vet=[500];
    n_val_fin = 1200; % fine dati per val
    %%%%%%%%%%%%
    
    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('stato.txt', 'w' );
        fprintf(fid, strt);
        fclose(fid);
           
        n = n_vet(num_n);
        
        % da Udati_id a Udati_id_sel:
        for jj=1:size(Dati(1,num_sens).t.sseh_yy_norm,2) 
            mm = Dati(1,num_sens).t.sseh_yy_norm(2,jj); % al variare di n, p=25
            mod_id = mm.Mnpf;
            coef_non_zer = find(mod_id.nb);         
            Udati_id_sel = [];
            for kk = 1:length(coef_non_zer)
                Udati_id_sel = [Udati_id_sel Dati(num_sens).t.Udati_id(:,coef_non_zer(kk))];  
            end
            U = Udati_id_sel(1:n,:); % <-ora contiene solo le col selez

            y = Dati(num_sens).t.y_id(1:n,:);
            
            %dati per val
            Uf = Udati_id_sel(n+1:n_val_fin,:);
            yf = Dati(num_sens).t.y_id(n+1:n_val_fin,:);
            dat_val = iddata(yf,Uf);


            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(1,jj) = struct('Mpem_BIC',Mpem_BIC,'Mpem_AICC',Mpem_AICC,'dat_val',dat_val,'n',n,'n_sel_sseh',mm.n_,'in_selez',coef_non_zer);

        end
        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

% 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);

%%
%%%%%%%%%%%%%%%%%%%% AICC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid = fopen('ValidazioneChiuso/pem_norm_Mpem_selez_AICC1500.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_AICC'];
    fprintf(fid, [str_nome_id, '\n']);
    
    disp(['Sensore: Dati(1,' num2str(num_sens) ').t.pem_norm(p,n).Mpem_AICC']);
    
    for k=1:30 % passo di predizione
        

        str_dati = [str_dati, num2str(k), '\t'];
        
        for num_n=3:3%1:len_n % 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_AICC;
                     
                    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);
