%-------------------------------------------------------------------------%
% 1SN - TP Optimisation                                                   %
% INP Toulouse - ENSEEIHT                                                 %
% Novembre 2020                                                           %
%                                                                         %
% Ce fichier contient le programme principal permettant l'estimation      %
% des parametres de la fonction de desintegration radioactive du          %
% carbone 14 par une approche des moindres carres.                        %                                      %
% Modele : A(t)= A0*exp(-lambda*t)                                        %
% Les algorithmes utilises pour la minimisation sont :                    %
% - l'algorithme de Gauss-Newton                                          %
% - l'algorithme de Newton.                                               %
%-------------------------------------------------------------------------%

clear
close all
clc
format shortG
taille_ecran = get(0,'ScreenSize');
L = taille_ecran(3);
H = taille_ecran(4);

% Initialisation

% Donnees
if exist('C14_results.txt','file')
    delete('C14_results.txt');
end
diary C14_results.txt
Ti = [ 500; 1000; 2000; 3000; 4000; 5000; 6300];
Ai = [14.5; 13.5; 12.0; 10.8;  9.9;  8.9;  8.0];
Donnees = [Ti, Ai];

% Estimation a priori des parametres du modele : beta0 = [A0, lambda]
beta0 = [10; 0.0001];   % Newton, Gauus-Newton, fminunc et leastsq converge
% beta0 = [15; 0.001];    % Newton, Gauss-Newton, fminunc et leastsq divergent
% beta0 = [15; 0.0005];   % Newton diverge, Gauss-Newton, fminunc et leastsq convergent
% beta0 = [10; 0.0005];   % Gauss-Newton converge

%% Calcul et affichage du modele initial ----------------------------------
xmin = 9; xmax = 20;
x = linspace(xmin, xmax, 100);

ymin = -0.0001; ymax = 0.0005;
y = linspace(ymin, ymax, 100);

[A0_plot, lambda_plot] = meshgrid(x, y);
[m , n] = size(A0_plot);
f_plot = zeros(m, n);

for i = 1:m
   for j = 1:n
      res_plot = residu_C14([A0_plot(i,j); lambda_plot(i,j)], Donnees);
      f_plot(i,j) = res_plot.'*res_plot /2;
   end
end

figure('Position',[0.1*L,0.1*H,0.8*L,0.8*H]); 
subplot(1,3,1) % Fonction f en 3D
    mesh(A0_plot, lambda_plot, f_plot)
    axis('square')
    title('Representation de la fonction f des moindres carres')
    xlabel('A_0');
    ylabel('\lambda')
    zlabel('f(A_0,\lambda)')

T = linspace(0,6500,100);
A = beta0(1)*exp(-beta0(2)*T);
txt_legend{1} = 'donnees';
txt_legend{2} = 'depart';

subplot(2, 3, 2) % Donnees + modele initial pour Gauss-Newton
    title('Desintegration radioactive du carbone 14 (Gauss-Newton)')
    axis([0 max(T) 0 18])
    xlabel('dur�e T')
    ylabel('radioactivit� A')
    hold on
    grid on
    plot(Ti, Ai, 'ok')
    plot(T, A)
    legend(txt_legend{1:2}, 'Location', 'SouthWest')

subplot(2, 3, 5) % Courbes de niveaux de f pour Gauss-Newton
    title('Recherche des parametres (Gauss-Newton)')
    xlabel('A_0')
    ylabel('\lambda')
    hold on
    contour(A0_plot, lambda_plot, f_plot, 100);
    plot(beta0(1),beta0(2),'ok')
    text(beta0(1),beta0(2),'  depart \beta^{(0)}')
    
subplot(2, 3, 3) % Donnees + modele initial pour Newton
    title('Desintegration radioactive du carbone 14 (Newton)')
    axis([0 max(T) 0 18])
    xlabel('dur�e T')
    ylabel('radioactivit� A')
    hold on
    grid on
    plot(Ti, Ai, 'ok')
    plot(T, A)
    legend(txt_legend{1:2}, 'Location', 'SouthWest')

subplot(2, 3, 6) % Courbes de niveaux de f pour Newton
    title('Recherche des parametres (Newton)')
    xlabel('A_0')
    ylabel('\lambda')
    hold on
    contour(A0_plot, lambda_plot, f_plot, 100);
    plot(beta0(1),beta0(2),'ok')
    text(beta0(1),beta0(2),'  depart \beta^{(0)}')

pause(0.5)

%% Algorithmes
% Choix du nombre d'iteration maximal

nb_iterations_max = 8;
txt_legend = cell(1,nb_iterations_max+2);
txt_legend{1} = 'donnees';
txt_legend{2} = 'depart';


%% Gauss-Newton
% -------------
% Initialisation de l'affichage
disp('Algorithme de Gauss-Newton')
disp('--------------------------------------------------------------------------------------------')
disp('     nb_iter         A0        lambda    ||f''(beta)||     f(beta)    ||delta||     exitflag ')
disp('--------------------------------------------------------------------------------------------')

% Calcul et affichage des valeurs initiales
res_beta = residu_C14(beta0, Donnees);
f_beta = 0.5*(res_beta.')*res_beta;
J_res_beta = J_residu_C14(beta0, Donnees);
norm_grad_f_beta = norm((J_res_beta.')*res_beta);
disp([0 beta0(1) beta0(2) norm_grad_f_beta f_beta]);

options = [sqrt(eps) 1.e-12 0];

for i = 1:nb_iterations_max
    txt_legend{i+2} = ['iteration ' num2str(i)];
	options(3) = i;
	[beta, norm_grad_f_beta, f_beta, norm_delta, k, exitflag] = ...
	Algo_Gauss_Newton(@(beta) residu_C14(beta, Donnees), ...
                      @(beta) J_residu_C14(beta, Donnees), ...
                      beta0, options);
	disp([k beta(1) beta(2) norm_grad_f_beta f_beta norm_delta exitflag])
    % Visualisation
    A = beta(1)*exp(-beta(2)*T);
	subplot(2, 3, 2)
	plot(T,A)
    legend(txt_legend{1:i+2})
    %  eval(['print -depsc fig_GN_courbe' int2str(i) '_C14'])
	subplot(2, 3, 5)
    plot(beta(1),beta(2),'ok')
	text(beta(1),beta(2),['  \beta^{(' num2str(i) ')}'])
  
	pause(0.5)
end

% Affichage des itérés de beta et sauvegarde des graphique
disp('--------------------------------------------------------------------------------------------')

%% Newton
% -------
% Initialisation de l'affichage
disp('Algorithme de Newton')
disp('--------------------------------------------------------------------------------------------')
disp('     nb_iter         A0        lambda    ||f''(beta)||     f(beta)    ||delta||     exitflag ')
disp('--------------------------------------------------------------------------------------------')

% Calcul et affichage des valeurs initiales
res_beta = residu_C14(beta0, Donnees);
f_beta = 0.5*(res_beta.')*res_beta;
J_res_beta = J_residu_C14(beta0, Donnees);
norm_grad_f_beta = norm((J_res_beta.')*res_beta);
disp([0 beta0(1) beta0(2) norm_grad_f_beta f_beta]);

options = [sqrt(eps) 1.e-12 0];

for i = 1:nb_iterations_max
    options(3) = i;
    % Algorithme de Newton
    [beta, norm_grad_f_beta, f_beta, norm_delta, k, exitflag] = ...
    Algo_Newton(@(beta) Hess_f_C14(beta, Donnees, ...
                                   @(beta) residu_C14(beta, Donnees), ...
                                   @(beta) J_residu_C14(beta, Donnees)), ...
                beta0, options);
    % Affichage des valeurs
    disp([k beta(1) beta(2) norm_grad_f_beta f_beta norm_delta exitflag])
    % Visualisation
    A = beta(1)*exp(-beta(2)*T);
    subplot(2, 3, 3)
    plot(T,A)
    legend(txt_legend{1:i+2})
    subplot(2, 3, 6)
    plot(beta(1),beta(2),'ok')
	text(beta(1),beta(2),['  \beta^{(' num2str(i) ')}'])
  
	pause(0.5)
end

% Affichage des iteres de beta et sauvegarde des courbes
disp('--------------------------------------------------------------------------------------------')
print('C14_figures','-dpng')
diary