From 6fcf549f9fd238dc333a7404203b60b0b097b16d Mon Sep 17 00:00:00 2001
From: Laurie Boffelli <laurie.boffelli@orange.fr>
Date: Mon, 24 Apr 2023 10:09:03 +0200
Subject: [PATCH] Correction Dynamics_Model (modification: fmincon)

---
 parmorceauxMax-git/Dynamics_Model.m | 59 ++++++++++++++++++++++++-----
 1 file changed, 49 insertions(+), 10 deletions(-)

diff --git a/parmorceauxMax-git/Dynamics_Model.m b/parmorceauxMax-git/Dynamics_Model.m
index a5d35bb..e2a7ca6 100644
--- a/parmorceauxMax-git/Dynamics_Model.m
+++ b/parmorceauxMax-git/Dynamics_Model.m
@@ -161,7 +161,13 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
             X_ini = [X_ini; new_X_ini];
             X_ini=obj.normalise_X(X_ini)
 
-            [new_X_variable,loglikelihood] = fminunc(f,X_ini);
+            lb = repmat(0,1,length(X_ini));
+            ub = repmat(1,1,length(X_ini));
+            A = [];
+            b = [];
+            Aeq = [];
+            beq = [];
+            [new_X_variable,loglikelihood] = fmincon(f,X_ini,A,b,Aeq,beq,lb,ub);
             obj.confidence = [obj.confidence, -loglikelihood / size(new_X_variable,1)];
             
             new_X_variable = obj.normalise_X(new_X_variable);
@@ -178,7 +184,7 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
                     obj.nb_times_fixed = [obj.nb_times_fixed; 0];
                 end
                 
-                if abs(obj.X(i) - new_X(i)) < 1e-3
+                if abs(obj.X(i) - new_X(i)) < 1e-3%1e-2
                     obj.nb_times_fixed(i) = obj.nb_times_fixed(i)+1;
                 else
                     obj.nb_times_fixed(i) = 0;
@@ -202,21 +208,21 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
             if new_N_fixed > N-3
                 new_N_fixed = N-3;
             end
-
+% 
             if new_N_fixed < 0
                 new_N_fixed = 0;
             end
-
-            if new_N_fixed < N-13
-                new_N_fixed = N-13;
+% 
+            if new_N_fixed < N-13 %N-8
+                new_N_fixed = N-13; %N-8;
             end
             
-            if size(new_X, 1) < 12
+            if size(new_X, 1) < 4 %12
                 new_N_fixed = size(new_X, 1) - 2;
                 if new_N_fixed < 0
                     new_N_fixed = 0;
                 end
-            elseif size(new_X, 1) == 12
+            elseif size(new_X, 1) == 4 %12
                 new_N_fixed = 0;
             end
 
@@ -262,8 +268,13 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
 
             
             X_ini = X_tmp(2:end-1)';
-            
-            [new_X,loglikelihood] = fminunc(f,X_ini);
+            lb = repmat(0,1,length(X_ini));
+            ub = repmat(1,1,length(X_ini));
+            A = [];
+            b = [];
+            Aeq = [];
+            beq = [];
+            [new_X,loglikelihood] = fmincon(f,X_ini,A,b,Aeq,beq,lb,ub);
             
             new_X = [0; obj.normalise_X(new_X); 1];
             obj.X = new_X;
@@ -288,6 +299,8 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
         function loglikelihood = logP_x_knowing_y_partialUpdate(obj, last_fixed_X, variable_X, alpha, sigma)
             N_variable = size(variable_X,1);
 
+ %           variable_X = obj.normalise_X(variable_X);
+
             % if N == 1, we only maximise P_cond_y(Y(1),X(1))
             % if N == 2, we maximise P_cond_y(Y(2),X(2)) * P_x_Markov_simple(X(2),X(1)) * P_cond_y(Y(1),X(1)) with P_x_Markov_simple(X(2),X(1)) = 0 if delta_X is negative, and uniform proba otherwise
 
@@ -415,11 +428,37 @@ classdef Dynamics_Model < matlab.mixin.Copyable %handle
                     ylabel(['Manifold coordinate Y_' i])
                     hold on
                     plot(obj.normalise_X(obj.X(start_frame:N)),obj.Y(start_frame:N,i),'kx')
+                    
+                end
+                %pour faire des gifs
+                %saveas(obj.display_screen, strcat('estimatedX_frame', int2str(N), '.png'))
+            end
+
+        end
+
+        function [im]=final_plot_estimated_X(obj)
+            if obj.display_screen ~= -1
+                N = size(obj.X,1);
+                start_frame = max(N - 50, 1);
+                
+                figure(obj.display_screen)
+                
+                
+                for i=1:obj.num_dim_Y
+                    subplot(obj.num_dim_Y,1,i)
+                    hold off
+                    im=imagesc([0 1],[obj.limits_Yi(i,1) obj.limits_Yi(i,2)],obj.marginal_X_Yi(:,:,i),'CDataMapping','scaled')
+                    %set(gca,'YDir','normal')
+                    hold on
+                    plot(obj.normalise_X(obj.X(start_frame:N)),obj.Y(start_frame:N,i),'kx')
+                    axis off
+                    
                 end
                 %pour faire des gifs
                 %saveas(obj.display_screen, strcat('estimatedX_frame', int2str(N), '.png'))
             end
         end
+
         function X = normalise_X(obj, X)
             X(X<0) = 0;
             X(X>1) = 1;
-- 
GitLab