Interface.m

classdef Interface < handle
    properties
        %% données d'entrée
        T %temps (ou numéro de frame)
        Y %observation : volume du ventricule gauche. Une observation (frame) par ligne de Y
        N
        
        %% résultats en sortie
        meilleure_sequence
        meilleure_sequence_models
        loglikelihood
        
        %% modèle statistique
        interpolant_filenames
        
        %% morceaux courants (cad non finis) A, B1, B2 et C
        dynamics_models
        first_frame_of_models
        
        %% paramètres de la méthode
        initial_alphas
        
        sigma_estimation
        sigma_test
        
        threshold_dynamics
        
        display_screen
        
        %% Viterbi
        transition
        
        proba_morceaux          % probabilité courante pour chaque morceau, avec son meilleur path Viterbi
        morceaux_finis          % pour chaque morceau courant : liste des morceaux finis du path. Cette liste contient une copie du morceau dans l'état où il était quand il a été terminé
        proba_morceaux_finis    % pour chaque morceau courant : liste des probabilités associées aux morceaux finis. Elle combine la proba de chaque observation du morceau et la proba de transition du morceau au suivant
        
        model_indexes           % pour chaque morceau courant : liste des indices des morceaux du path Viterbi pour chaque frame
        scores_finis            % pour chaque morceau courant : liste des proba d'observation pour chaque frame (uniquement les frames des morceaux finis, car le morceau courant n'est pas encore stable et les x peuvent encore changer)
        
    end
    methods
        function obj = Interface(interpolant_filenames, display_screen)
            
            obj.interpolant_filenames = interpolant_filenames;
            
            obj.transition = [  0.5 0.5 0 0;   % départ de A
                                0 1/3 1/3 1/3; % départ de B1
                                0 0 0.5 0.5;   % départ de B2
                                0 0 0 1        % départ de C
                             ];
            
            obj.display_screen = display_screen;
        end
        function process_new_sequence(obj, filename, initial_alphas, sigma_estimation, sigma_test, threshold_dynamics)
            obj.threshold_dynamics = threshold_dynamics;
            
            obj.dynamics_models = {NaN, NaN, NaN, NaN}; % modèles courants A B1 B2 et C
            obj.first_frame_of_models = [1 NaN NaN NaN]; % frame où chaque modèle commence
            
            if obj.display_screen == -1
                obj.dynamics_models{1} = Dynamics_Model(obj.interpolant_filenames{1}, obj.display_screen);
            else
                obj.dynamics_models{1} = Dynamics_Model(obj.interpolant_filenames{1}, obj.display_screen+1);
            end
            
            obj.proba_morceaux = [1 0 0 0]; % proba des modèles A, B1, B2, et C
            
            obj.T = [];
            obj.Y = [];
            obj.N = [];
            
            obj.scores_finis = {[], [], [], []};
            obj.model_indexes = {[], [], [], []};
            
            obj.initial_alphas = initial_alphas;
            obj.sigma_estimation = sigma_estimation;
            obj.sigma_test = sigma_test;
            
            obj.morceaux_finis = {[], [], [], []};
            obj.proba_morceaux_finis = {[], [], [], []};
            
            obj.read_data(filename);
            obj.N = size(obj.T,1);
            
            obj.dynamics_models{1}.initialise(0, initial_alphas(1), sigma_estimation, sigma_test);
            %ligne 14
            for frame=1:obj.N
                obj.process_frame(frame);
                
                if obj.display_screen ~= -1
                    for model=1:4
                        if obj.proba_morceaux(model) == 0 || isnan(obj.proba_morceaux(model))
                            continue
                        end
                        obj.dynamics_models{model}.plot_estimated_X();
                    end
                end
            end
            
            
            %% choix de la meilleure sequence (peut être pour plus tard: toujours choisir la séquence qui se termine sur le modèle C ?)
            [max_proba, indice] = max(obj.proba_morceaux);
            obj.meilleure_sequence = obj.model_indexes{indice}
            obj.meilleure_sequence_models = [obj.morceaux_finis{indice} obj.dynamics_models{indice}];
            
            %% calcul les scores pour les frames du dernier morceau
            % mise à jour des x pour bien les étaler
            obj.meilleure_sequence_models(end).spread_X();
            new_scores = obj.meilleure_sequence_models(end).compute_loglikelihood_dynamics([1:obj.N-obj.first_frame_of_models(indice)+1]);
            obj.loglikelihood = [obj.scores_finis{indice} new_scores];
            
            %% à faire : afficher la séquence des modèles et la séquence de scores comme dans plot
            close all
            for i=1:length(obj.meilleure_sequence_models)
                obj.meilleure_sequence_models(i).plot_estimated_X();
            end
%%%         save plot
     
        end
        function read_data(obj, filename)
            obj.Y =load(filename);
            %obj.T = [1:75]';
            %obj.Y = obj.Y / max(obj.Y);
            %obj.Y = repmat(obj.Y, 3, 1);
            obj.T = [1:25]';
            obj.Y = obj.Y / max(obj.Y);
        end
        function process_frame(obj, frame_number)
            
            frame_number
            
            Yt = obj.Y(frame_number,:);
            Tt = obj.T(frame_number);
            
            %% tableaux temporaires dans lesquels on prépare la prochaine version des tableaux de la classe sans les modifier avant d'avoir terminé de traiter la frame pour tous les morceaux
            new_proba_morceaux = [NaN NaN NaN NaN];
            new_morceaux_finis = {[], [], [], []};
            new_proba_morceaux_finis = {[], [], [], []};
            
            new_models = {NaN, NaN, NaN, NaN};
            new_first_frame_of_models = [NaN NaN NaN NaN];
            
            new_scores_finis = {[], [], [], []};
            new_model_indexes = {[], [], [], []};
            
            
            %% boucle sur les 4 modèles possibles
            for morceau = 1:4
                
                % on va regarder toutes les branches qui arrivent sur le modèle
                
                %% tableaux temporaires pour le morceau courant et toutes les branches possibles
                modele_branche = {NaN, NaN, NaN, NaN};               % copie des modèles pour leur montrer la nouvelle frame, sans modifier ceux de obj.dynamics_models
                first_frame_of_model_branche = [NaN NaN NaN NaN];    % pour les copies des modèles
                proba_obs_branche = [NaN NaN NaN NaN];               % proba d'observation donnée par les copies des modèles ayant vu la frame courante
                nouveau_morceau_fini_branche = {NaN, NaN, NaN, NaN}; % si une branche implique un nouveau morceau fini, on le sauvegarde ici en attendant Viterbi et le choix de la branche (cad si on le garde ou pas)
                proba_branche = [NaN NaN NaN NaN];                   % proba totale de chaque branche
                proba_morceau_prec = [NaN NaN NaN NaN];              % proba totale des frames vues par le morceau qui était courant à la frame précédente (donc d'indice=branche). Ces proba sont combinées avec les probas de transition au sein du morceau
                new_scores_branche = {[] [] [] []};                  % score des frames vues par le morceau qui était courant à la frame précédente
                
                if frame_number == 1
                    if obj.proba_morceaux(morceau) == 0
                        continue
                    end
                    
                    new_models{morceau} = copy(obj.dynamics_models{morceau});
                    new_first_frame_of_models(morceau) = obj.first_frame_of_models(morceau);
                    
                    new_models{morceau}.process_new_frame(Yt, Tt);
                    proba_obs = new_models{morceau}.compute_loglikelihood_dynamics(1);
                    
                    new_proba_morceaux(morceau) = log(obj.proba_morceaux(morceau)) + proba_obs;
                    new_morceaux_finis{morceau} = [];
                    
                    new_model_indexes{morceau} = [morceau];
                else
                    %% on calcule les probas de toutes les branches qui arrivent sur le morceau courant
                    for branche = 1:4
                        trans = obj.transition(branche, morceau);
                        
                        % on commence par ignorer les branches impossibles
                        if trans == 0
                            continue
                        end
                        if obj.proba_morceaux(branche) == 0 || isnan(obj.proba_morceaux(branche))
                            continue
                        end
                        
                        if branche == morceau   % cas sans changement de modèle, donc le modèle est déjà initialisé
                            modele_branche{branche} = copy(obj.dynamics_models{branche});
                            first_frame_of_model_branche(branche) = obj.first_frame_of_models(morceau);
                        else            % cas avec changement de modèle, donc il faut créer et initialiser le modèle
                            %% on crée le modèle et on l'initialise
                            if obj.display_screen == -1
                                modele_branche{branche} = Dynamics_Model(obj.interpolant_filenames{morceau}, obj.display_screen);
                            else
                                modele_branche{branche} = Dynamics_Model(obj.interpolant_filenames{morceau}, obj.display_screen+morceau);
                            end
                            modele_branche{branche}.initialise(0, obj.initial_alphas(morceau), obj.sigma_estimation, obj.sigma_test);
                            %modele_branche{branche}.initialise(0, obj.dynamics_models{branche}.current_alpha, obj.sigma_estimation, obj.sigma_test);
                            first_frame_of_model_branche(branche) = frame_number;
                        end
                            
                        %% on donne la frame au modèle, et on récupère la proba d'observation (dans la fenêtre)
                        modele_branche{branche}.process_new_frame(Yt, Tt);
                        %print branche et morceau
                        BRANCHE=branche
                        MORCEAU=morceau
                        proba_obs_branche(branche) = modele_branche{branche}.compute_loglikelihood_dynamics(frame_number-first_frame_of_model_branche(branche)+1);
                        
                        %% on garde en mémoire les morceaux qui sont finis, pour l'affichage final des X
%%%%                    %%%%%%%%%%%%%%%% a faire : mettre à jour les x pour bien les étaler
                        if branche == morceau
                            nouveau_morceau_fini_branche{branche} = 0;
                        else
                            nouveau_morceau_fini_branche{branche} = copy(obj.dynamics_models{branche}); % ce modèle n'a pas encore vu la frame courante, on peut donc le sauvegader dans son état à la frame précédente
%%%%                        %%%%%% mettre à jour la copie ici
                            nouveau_morceau_fini_branche{branche}.spread_X();
                        end
                        
                        %% on récupère les scores des frames vues par le modèle
                        if branche == morceau
                            %% on récupère les scores pour les frames précédentes, qui ont pu bouger avec le argmax
                            proba_morceau_prec(branche) = sum(modele_branche{branche}.compute_loglikelihood_dynamics([1:frame_number-first_frame_of_model_branche(branche)]));
                            proba_morceau_prec(branche) = proba_morceau_prec(branche) + log(trans) * (frame_number-first_frame_of_model_branche(branche)-1);
                        else
                            %% on récupère les scores pour les frames précédentes, qui ont pu bouger avec le argmax
                            new_scores_branche{branche} = nouveau_morceau_fini_branche{branche}.compute_loglikelihood_dynamics([1:frame_number-obj.first_frame_of_models(branche)]);
                            proba_morceau_prec(branche) = sum(new_scores_branche{branche});
                            trans_morceau_prec = obj.transition(branche, branche);
                            proba_morceau_prec(branche) = proba_morceau_prec(branche) + log(trans_morceau_prec) * (frame_number-obj.first_frame_of_models(branche)-1);
                        end
          
                        proba_branche(branche) = sum(obj.proba_morceaux_finis{branche}) + proba_morceau_prec(branche) + log(trans) + proba_obs_branche(branche);
                    end
                    
                    %% on choisi la meilleure branche
                    [max_val, indice] = max(proba_branche);
                    
                    if isnan(max_val)
                        continue
                    end
                    
                    new_proba_morceaux(morceau) = max_val;
                    
                    new_models{morceau} = modele_branche{indice};
                    new_first_frame_of_models(morceau) = first_frame_of_model_branche(indice);
                    
                    if nouveau_morceau_fini_branche{indice} ~= 0
                        new_morceaux_finis{morceau} = [obj.morceaux_finis{indice} nouveau_morceau_fini_branche{indice}];
                        trans = obj.transition(indice, morceau);
                        new_proba_morceaux_finis{morceau} = [obj.proba_morceaux_finis{indice} proba_morceau_prec(indice) + log(trans)];
                        new_scores_finis{morceau} = [obj.scores_finis{indice} new_scores_branche{indice}];
                    else
                        new_morceaux_finis{morceau} = obj.morceaux_finis{indice};
                        new_proba_morceaux_finis{morceau} = obj.proba_morceaux_finis{indice};
                        new_scores_finis{morceau} = obj.scores_finis{indice};
                    end
                    
                    new_model_indexes{morceau} = [obj.model_indexes{indice} morceau];
                end
            end
            
            obj.proba_morceaux = new_proba_morceaux;
            obj.dynamics_models = new_models;
            obj.first_frame_of_models = new_first_frame_of_models;
            
            obj.scores_finis = new_scores_finis;
            
            obj.morceaux_finis = new_morceaux_finis;
            obj.proba_morceaux_finis = new_proba_morceaux_finis;
            
            obj.model_indexes = new_model_indexes;
         
            
            %% réflexions pour plus tard
            %% gérer : alpha différents pour les différents modèles (alpha_ini + marges)
            %%         fin des modèles : duppliquer les dernières colonnes pour que x ait le temps de se fixer ? (avec pénalité lors de la comparaison des modèles)
            %%         ou alors : on n'attend pas que les x soient stables
            %%         initialiser le modèle concurrent à chaque fois (tant qu'on n'a pas basculé dessus)    
        end
        function plot(obj) %% cette fonction est à modifier car maintenant nous avons plusieurs modèles / morceaux
            obj.dynamics_model.plot_estimated_X();
            
            figure(obj.display_screen);
            n = length(obj.loglikelihood_dynamics_window);
            subplot(4,1,1)
            hold off
            plot(obj.T(1:n),obj.Y(1:n,1))
            title('Pose representation')
            xlabel('frame number')
            ylabel('1st dimension')
            XL = xlim;
            
            subplot(4,1,2)
            hold off
            plot(obj.T(1:n),obj.loglikelihood_pose)
            hold on
            indexes_p = (obj.loglikelihood_pose < obj.threshold_pose);
            plot(obj.T(indexes_p),obj.loglikelihood_pose(indexes_p),'.r')
            title('Likelihood of the pose to be normal')
            xlabel('Frame number')
            ylabel('Log likelihood')
            xlim(XL)
            YL = ylim;
            YL1 = max(YL(1),-5);
            ylim([YL1 YL(2)])
            hold on
            
            subplot(4,1,3)
            hold off
            plot(obj.T(1:n),obj.loglikelihood_dynamics_window)
            hold on
            indexes_w = (obj.loglikelihood_dynamics_window < obj.threshold_dynamics);
            plot(obj.T(indexes_w),obj.loglikelihood_dynamics_window(indexes_w),'.r')
            title('Likelihood of the dynamics to be normal, computed in the dynamic window')
            xlabel('Frame number')
            ylabel('Log likelihood')
            xlim(XL)
            YL = ylim;
            YL1 = max(YL(1),-20);
            ylim([YL1 YL(2)])
            hold on
            
            subplot(4,1,4)
            hold off
            n2 = length(obj.loglikelihood_dynamics);
            plot(obj.T(1:n2),obj.loglikelihood_dynamics)
            hold on
            indexes_d = (obj.loglikelihood_dynamics < obj.threshold_dynamics);
            plot(obj.T(indexes_d),obj.loglikelihood_dynamics(indexes_d),'.r')
            title('Likelihood of the dynamics to be normal, computed for each converged frame')
            xlabel('Frame number')
            ylabel('Log likelihood')
            xlim(XL)
            YL = ylim;
            YL1 = max(YL(1),-20);
            YL2 = min(YL(2),10);
            ylim([YL1 YL2])
            hold on
            
            subplot(4,1,1)
            hold on
            plot(obj.T(indexes_p),obj.Y(indexes_p,1),'.r')
            %plot(obj.T(indexes_w),obj.Y(indexes_w,1),'.r')
            plot(obj.T(indexes_d),obj.Y(indexes_d,1),'.r')
            %indexes = indexes_p == 0 & indexes_w == 0 & [indexes_d ones(1,n-n2)] == 0;
            indexes = indexes_p == 0 & [indexes_d ones(1,n-n2)] == 0;
            plot(obj.T(indexes),obj.Y(indexes,1),'.g')
        end
    end
end