clear all
close all

%% Données à modifier

patient = 18 ;

% Définir le nombre de jeux de paramètres
n_samples = 100 ;

% Définir les bornes de l'espace parcouru

E1_min = 0.001 ;
E1_max = 0.1 ;
E2_min = 0.001 ;
E2_max = 0.1 ;
E3_min = 0.001 ;
E3_max = 0.1 ;

% Définir le nom du dossier
folder_name = ['patient',num2str(patient),'_',num2str(n_samples),'_3params'];

% Créer le dossier si celui-ci n'existe pas déjà
if ~exist(folder_name, 'dir')
    mkdir(folder_name);
end

ell = table(E1_min,E1_max,E2_min,E2_max,E3_min,E3_max) ;
file_name = fullfile(folder_name,'bornes.txt');
writetable(ell,file_name);

% coordonnées / angles
filename = ['coord_',num2str(patient),'.txt']; % Remplace par le chemin de ton fichier
coord = readtable(filename, 'Delimiter', ','); % Lire le fichier en spécifiant que le séparateur est une virgule

if patient == 5
    X_mes(1) = coord.X_max(7) ; % point 7 correspond à noeud 6 donc 1er dans résultat simu
    X_mes(2) = coord.X_max(6) ; % point 6 correspond à noeud 14 donc 2e dans résultat simu (si on exclut le noeud 11)
    X_mes(3) = coord.X_max(9) ; % point 9 correspond à noeud 18 donc 3e dans résultat simu

    Y_mes(1) = coord.Y_max(7) ; 
    Y_mes(2) = coord.Y_max(6) ; 
    Y_mes(3) = coord.Y_max(9) ;
elseif patient == 13
    X_mes(1) = coord.X_max(7) ; % point 7 correspond à noeud 6 donc 1er dans résultat simu
    X_mes(2) = coord.X_max(6) ; % point 6 correspond à noeud 14 donc 2e dans résultat simu (si on exclut le noeud 11)
    X_mes(3) = coord.X_max(9) ; % point 9 correspond à noeud 18 donc 3e dans résultat simu

    Y_mes(1) = coord.Y_max(7) ; 
    Y_mes(2) = coord.Y_max(6) ; 
    Y_mes(3) = coord.Y_max(9) ;

elseif patient == 20
    X_mes(1) = coord.X_max(7) ; % point 7 correspond à noeud 9 donc 1er dans résultat simu
    X_mes(2) = coord.X_max(9) ; % point 9 correspond à noeud 14 donc 2e dans résultat simu (si on exclut le noeud 11)
    X_mes(3) = coord.X_max(6) ; % point 6 correspond à noeud 20 donc 3e dans résultat simu

    Y_mes(1) = coord.Y_max(7) ;
    Y_mes(2) = coord.Y_max(9) ;
    Y_mes(3) = coord.Y_max(6) ;

elseif patient == 18
    X_mes(1) = coord.X_max(7) ; % point 7 correspond à noeud 9 donc 1er dans résultat simu
    X_mes(2) = coord.X_max(9) ; % point 9 correspond à noeud 14 donc 2e dans résultat simu (si on exclut le noeud 11)
    X_mes(3) = coord.X_max(6) ; % point 6 correspond à noeud 20 donc 3e dans résultat simu

    Y_mes(1) = coord.Y_max(7) ;
    Y_mes(2) = coord.Y_max(9) ;
    Y_mes(3) = coord.Y_max(6) ;
end

% grand axe / petit axe de l'ellipse
filename = ['ell_',num2str(patient),'.txt']; % Remplace par le chemin de ton fichier
ellipse = readtable(filename, 'Delimiter', ','); % Lire le fichier en spécifiant que le séparateur est une virgule

% module d'Young et épaisseurs
filename = ['thickness_',num2str(patient),'.txt']; % Remplace par le chemin de ton fichier
thick_young = readtable(filename, 'Delimiter', ','); % Lire le fichier en spécifiant que le séparateur est une virgule

% coordonnées des noeuds d'intérêt
filename = ['coord_simu_',num2str(patient),'.txt']; % Remplace par le chemin de ton fichier
coord_simu = readtable(filename, 'Delimiter', ','); % Lire le fichier en spécifiant que le séparateur est une virgule

% %% Etude paramétrique à la main
% X = [Ef_RA Ef_LM Ef_LA ; 10*Ef_RA Ef_LM Ef_LA ; Ef_RA/10 Ef_LM Ef_LA ; Ef_RA 10*Ef_LM Ef_LA ; Ef_RA Ef_LM/10 Ef_LA ; Ef_RA Ef_LM 10*Ef_LA ; Ef_RA Ef_LM Ef_LA/10 ];
% filename_0 = [num2str(patient),'_0'];
% nom_fichier_inp = append(filename_0,'.inp');
% 
% for i=1:7
%     filename = append(filename_0,num2str(i));
%     nom_fichier_out = append(filename,'_out');
%     nom_fichier_inp_out = append(nom_fichier_out,'.inp');
%     nom_fichier_dat_out = append(nom_fichier_out,'.dat');
%     nom_fichier_dat = append(filename,'.dat');
%     
%     modify_INP(nom_fichier_inp, nom_fichier_inp_out, X(i,1),X(i,2),X(i,3));
%     RunINP(nom_fichier_inp_out);  
%     [U0,Ux0,Uy0] = readAbaqusOutput(nom_fichier_dat_out);
% 
% %     Y_sim(1) = coord_simu.Var3(1)+Uy0(1) ;
% %     Y_sim(2) = coord_simu.Var3(2)+Uy0(3) ;
% %     Y_sim(3) = coord_simu.Var3(3)+Uy0(4) ;
% 
%     Y(i,:) = [coord_simu.Var3(1)+Uy0(1) coord_simu.Var3(2)+Uy0(3) coord_simu.Var3(3)+Uy0(4)] ;
% end

%% Création de la base de données
% % Définir le nombre de jeux de paramètres
% n_samples = 1;
% 
% % Définir les bornes de l'espace parcouru
% E1_min = 0.001 ;
% E1_max = 0.1 ;
% E2_min = 0.0001 ;
% E2_max = 0.1 ;
% E3_min = 0.005 ;
% E3_max = 0.1 ;

% Générer aléatoirement des points dans les plages définies
E1_values = E1_min + (E1_max - E1_min) * rand(n_samples, 1);
E2_values = E2_min + (E2_max - E2_min) * rand(n_samples, 1);
E3_values = E3_min + (E3_max - E3_min) * rand(n_samples, 1);

% Combinaison des paramètres générés
X = [E1_values, E2_values, E3_values];
%X = [E1_values, E2_values];

filename_0 = [num2str(patient),'_0'];
nom_fichier_inp = append(filename_0,'.inp');

for i=1:n_samples
    filename = append(filename_0,num2str(i));
    nom_fichier_out = append(filename,'_out');
    nom_fichier_inp_out = append(nom_fichier_out,'.inp');
    nom_fichier_dat_out = append(nom_fichier_out,'.dat');
    nom_fichier_dat = append(filename,'.dat');
    
    modify_INP(nom_fichier_inp, nom_fichier_inp_out, X(i,1),X(i,2),X(i,3));
    RunINP(nom_fichier_inp_out);  
    [U0,Ux0,Uy0] = readAbaqusOutput(nom_fichier_dat_out);

%     X_sim_0(1) = coord_simu.Var2(1)+Ux0(1) ;
%     X_sim_0(2) = coord_simu.Var2(2)+Ux0(3) ;
%     X_sim_0(3) = coord_simu.Var2(3)+Ux0(4) ;
%     Y_sim(1) = coord_simu.Var3(1)+Uy0(1) ;
%     Y_sim(2) = coord_simu.Var3(2)+Uy0(3) ;
%     Y_sim(3) = coord_simu.Var3(3)+Uy0(4) ;

    if patient == 5
        Y(i,:) = [coord_simu.Var2(1)+Ux0(1) coord_simu.Var2(2)+Ux0(2) coord_simu.Var2(3)+Ux0(3) coord_simu.Var3(1)+Uy0(1) coord_simu.Var3(2)+Uy0(2) coord_simu.Var3(3)+Uy0(3)] ;
        erreur(i) = ((Y(i,4)-Y_mes(1))^2+(Y(i,5)-Y_mes(2))^2)^0.5; % uniquement les points 6 et 7 e Y
    elseif patient == 13
        Y(i,:) = [coord_simu.Var2(1)+Ux0(1) coord_simu.Var2(2)+Ux0(2) coord_simu.Var2(3)+Ux0(3) coord_simu.Var3(1)+Uy0(1) coord_simu.Var3(2)+Uy0(2) coord_simu.Var3(3)+Uy0(3)] ;
        erreur(i) = ((Y(i,4)-Y_mes(1))^2+(Y(i,5)-Y_mes(2))^2)^0.5; % uniquement les points 6 et 7 e Y
    elseif patient == 20
        Y(i,:) = [coord_simu.Var2(1)+Ux0(1) coord_simu.Var2(2)+Ux0(2) coord_simu.Var2(3)+Ux0(3) coord_simu.Var3(1)+Uy0(1) coord_simu.Var3(2)+Uy0(2) coord_simu.Var3(3)+Uy0(3)] ;
        erreur(i) = ((Y(i,2)-X_mes(2))^2+(Y(i,4)-Y_mes(1))^2+(Y(i,5)-Y_mes(2))^2+(Y(i,6)-Y_mes(3))^2)^0.5; % Ux pour le point 8 et Uy pour 6, 7 et 8
    elseif patient == 18
        Y(i,:) = [coord_simu.Var2(1)+Ux0(1) coord_simu.Var2(2)+Ux0(2) coord_simu.Var2(3)+Ux0(3) coord_simu.Var3(1)+Uy0(1) coord_simu.Var3(2)+Uy0(2) coord_simu.Var3(3)+Uy0(3)] ;
        erreur(i) = ((Y(i,4)-Y_mes(1))^2+(Y(i,5)-Y_mes(2))^2)^0.5; % uniquement les points 6 et 7 e Y
    end
    
end

figure
plot(erreur)
file_name = fullfile(folder_name,'erreur.png');
saveas(gcf, file_name);

%%
% Créer un maillage pour la surface
x = X(:,1);
y = X(:,2);
z = erreur'; % Utiliser l'erreur comme hauteur de la surface

% Pour tracer une surface, il faut transformer x et y en grilles
[Xgrid, Ygrid] = meshgrid(unique(x), unique(y));

% Interpolation pour créer les valeurs correspondantes de z
Zgrid = griddata(x, y, z, Xgrid, Ygrid);

% Tracer la surface
figure
surf(Xgrid, Ygrid, Zgrid,'EdgeColor', 'none');

% Ajouter des labels
xlabel('E_1');
ylabel('E_2');
zlabel('Error');
title('Surface Plot of Error');
file_name = fullfile(folder_name,'E1_E2.png');
saveas(gcf, file_name);

% Créer un maillage pour la surface
x = X(:,2);
y = X(:,3);
z = erreur'; % Utiliser l'erreur comme hauteur de la surface

% Pour tracer une surface, il faut transformer x et y en grilles
[Xgrid, Ygrid] = meshgrid(unique(x), unique(y));

% Interpolation pour créer les valeurs correspondantes de z
Zgrid = griddata(x, y, z, Xgrid, Ygrid);

% Tracer la surface
figure
surf(Xgrid, Ygrid, Zgrid,'EdgeColor', 'none');

% Ajouter des labels
xlabel('E_2');
ylabel('E_3');
zlabel('Error');
title('Surface Plot of Error');
file_name = fullfile(folder_name,'E2_E3.png');
saveas(gcf, file_name);

% Créer un maillage pour la surface
x = X(:,1);
y = X(:,3);
z = erreur'; % Utiliser l'erreur comme hauteur de la surface

% Pour tracer une surface, il faut transformer x et y en grilles
[Xgrid, Ygrid] = meshgrid(unique(x), unique(y));

% Interpolation pour créer les valeurs correspondantes de z
Zgrid = griddata(x, y, z, Xgrid, Ygrid);

% Tracer la surface
figure
surf(Xgrid, Ygrid, Zgrid,'EdgeColor', 'none');

% Ajouter des labels
xlabel('E_1');
ylabel('E_3');
zlabel('Error');
title('Surface Plot of Error');
file_name = fullfile(folder_name,'E1_E3.png');
saveas(gcf, file_name);

%% Réseau de neurones
net = feedforwardnet(10); % 10 neurones dans la couche cachée
net = train(net, X', Y');

%% Sauvegarde
file_name_txt = fullfile(folder_name, 'mon_reseau_neurones.mat');
save(file_name,'net');
file_name_txt = fullfile(folder_name, 'Y.mat');
save(file_name,'Y');
file_name_txt = fullfile(folder_name, 'X.mat');
save(file_name,'X');
file_name_txt = fullfile(folder_name, 'erreur.mat');
save(file_name,'erreur');


%% Optimisation
% Bornes inférieures et supérieures des paramètres
lb = [E1_min, E2_min, E3_min];
ub = [E1_max, E2_max, E3_max];

% Paramètres initiaux pour l'optimisation
init_params = [mean([E1_min, E1_max]), mean([E2_min, E2_max]),mean([E3_min, E3_max])];

% Options pour l'optimisation
options = optimoptions('fmincon', 'Algorithm', 'sqp');

% Optimisation des paramètres E1, E2, E3 pour minimiser l'écart avec les données expérimentales
[E_opt, fval] = fmincon(@(E) objective_function(E, net, X_mes(3), Y_mes(1),Y_mes(2)), init_params, [], [], [], [], lb, ub, [], options);

% Affichage des résultats
fprintf('Les paramètres optimisés sont : E_RA = %f, E_LM = %f, E_LA = %f\n', E_opt(1), E_opt(2), E_opt(3));


%% Vérification
filename = append(filename_0,'_verif3');
nom_fichier_out = append(filename,'_out');
nom_fichier_inp_out = append(nom_fichier_out,'.inp');
nom_fichier_dat_out = append(nom_fichier_out,'.dat');
nom_fichier_dat = append(filename,'.dat');

modify_INP(nom_fichier_inp, nom_fichier_inp_out, E_opt(1), E_opt(2), E_opt(3));
RunINP(nom_fichier_inp_out);
[U0,Ux0,Uy0] = readAbaqusOutput(nom_fichier_dat_out);

X_sim(1) = coord_simu.Var2(1)+Ux0(1) ;
X_sim(2) = coord_simu.Var2(2)+Ux0(2) ;
X_sim(3) = coord_simu.Var2(3)+Ux0(3) ;

Y_sim(1) = coord_simu.Var3(1)+Uy0(1) ;
Y_sim(2) = coord_simu.Var3(2)+Uy0(2) ;
Y_sim(3) = coord_simu.Var3(3)+Uy0(3) ;

figure
plot([0,max(X_mes)],[0,max(X_mes)])
hold on
plot(X_mes,X_sim,'gx')
hold on
plot(Y_mes,Y_sim,'rx')
hold on
test = net(E_opt') ;
plot(X_sim,test(1:3),'go');
hold on
plot(Y_sim,test(4:6),'ro');
legend('','Xmes/Xsim','Ymes/Ysim','Xsim/Xnet','Ysim/Ynet')
file_name = fullfile(folder_name, 'comparaison_exp_simu_net.png');
saveas(gcf,file_name);

figure
plot(coord_simu.Var2(:),coord_simu.Var3(:),'mo')
hold on
plot(X_mes,Y_mes,'gx')
hold on
plot(X_sim,Y_sim,'ro')
hold on
test = net(E_opt') ;
plot(test(1:3),test(4:6),'bx');
xlim([0 100]);
ylim([0 100]);
legend('coord init','coord def','coord simu','coord net')
file_name = fullfile(folder_name, 'comparaison_exp_simu_net_2.png');
saveas(gcf,file_name);

% %% Sensibilité
% % Valeurs de référence des paramètres d'entrée
% E_ref = E_opt';
% 
% % Initialiser les variations (par exemple, augmenter de 10% chaque paramètre)
% delta = 0.3 * E_ref;
% 
% % Pour chaque paramètre, on le fait varier d'une certaine quantité et on mesure les sorties
% E1_varied = E_ref + [delta(1), 0]';  % Variation de E1
% Y1_pred = net(E1_varied);  % Prédictions correspondantes pour Ua, Ub, Uc
% 
% E2_varied = E_ref + [0, delta(2)]';  % Variation de E2
% Y2_pred = net(E2_varied);
% 
% % Comparaison des résultats avec la sortie pour les paramètres de référence
% Y_ref = net(E_ref);  % Prédiction avec les paramètres de référence
% 
% % Affichage des résultats
% fprintf('Variation de E1 : Uya = %f, Uyb = %f, Uyc = %f\n', Y1_pred(4), Y1_pred(5), Y1_pred(6));
% fprintf('Variation de E2 : Uya = %f, Uyb = %f, Uyc = %f\n', Y2_pred(4), Y2_pred(5), Y2_pred(6));
% 
% % Comparaison avec la prédiction de référence
% fprintf('Valeurs de référence : Ua = %f, Ub = %f, Uc = %f\n', Y_ref(4), Y_ref(5), Y_ref(6));
% 
% % Affichage des résultats
% fprintf('Variation de E1 : Uxa = %f, Uxb = %f, Uxc = %f\n', Y1_pred(1), Y1_pred(2), Y1_pred(3));
% fprintf('Variation de E2 : Uxa = %f, Uxb = %f, Uxc = %f\n', Y2_pred(1), Y2_pred(2), Y2_pred(3));
% 
% % Comparaison avec la prédiction de référence
% fprintf('Valeurs de référence : Uxa = %f, Uxb = %f, Uxc = %f\n', Y_ref(1), Y_ref(2), Y_ref(3));
%