clear all
close all

%% Données à modifier

patient = 5 ;

type = 1 ;
complement = '_'; %modèle elliptique

% type = 2 ;
% complement = '_realiste_'; %modèle réaliste

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

% Définir les bornes de l'espace parcouru
% E1 = E_RA ; E2 = E_LM

E1_min = 0.1 ;
E1_max = 5 ;
E2_min = 0.1 ;
E2_max = 5 ;
E3 = 1 ;

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

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

% On range dans X_mes et Y_mes les coordonnées à gonflement max, dans
% l'ordre où les résultats de simu arrivent (= par ordre croissant des
% numéros de noeuds)

if patient == 5 && type==1
    % 7 = noeud 13 ; 6 = noeud 5 ; 9 = noeud 20
    X_mes(1) = coord.X_max3(6) ; % point 6 correspond à noeud 5 donc 1er dans résultat simu
    X_mes(2) = coord.X_max3(7) ; % point 7 correspond à noeud 13 donc 2e dans résultat simu 
    X_mes(3) = coord.X_max3(9) ; % peu importe, pas besoin du point 9 qui est confondu avec 6, on choisit le 14

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

% elseif patient == 5 && type ==2
%     X_mes(1) = coord.X_max(7) ; % point 7 correspond à noeud 2 donc 1er dans résultat simu
%     X_mes(2) = coord.X_max(6) ; % point 6 correspond à noeud 7 donc 2e dans résultat simu (si on exclut le noeud 11)
%     X_mes(3) = coord.X_max(9) ; % point 9 correspond à noeud 10 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
    % 7 = noeud 12 ; 6 = noeud 17 ; 9 = noeud 7
    X_mes(1) = coord.X_max3(9) ; % point 9 correspond à noeud 7 donc 1er dans résultat simu
    X_mes(2) = coord.X_max3(7) ; % point 7 correspond à noeud 12 donc 2e dans résultat simu
    X_mes(3) = coord.X_max3(6) ; % point 6 correspond à noeud 17 donc 3e dans résultat simu

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

elseif patient == 20 && type == 1
    % 7 = noeud 9 ; 6 = noeud 20 ; 9 = noeud 14
    X_mes(1) = coord.X_max3(7) ; % point 7 correspond à noeud 9 donc 1er dans résultat simu
    X_mes(2) = coord.X_max3(9) ; % point 9 correspond à noeud 14 donc 2e dans résultat simu 
    X_mes(3) = coord.X_max3(6) ; % point 6 correspond à noeud 20 donc 3e dans résultat simu

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

% elseif patient == 20 && type == 2
%     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 == 18
    % 7 = noeud 6 ; 6 = noeud 14 ; 9 = noeud 18
    X_mes(1) = coord.X_max3(7) ; % point 7 correspond à noeud 6 donc 1er dans résultat simu
    X_mes(2) = coord.X_max3(6) ; % point 6 correspond à noeud 14 donc 2e dans résultat simu 
    X_mes(3) = coord.X_max3(9) ; % point 9 correspond à noeud 18 donc 3e dans résultat simu

    Y_mes(1) = coord.Y_max3(7) ;
    Y_mes(2) = coord.Y_max3(6) ;
    Y_mes(3) = coord.Y_max3(9) ;
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',complement,num2str(patient),'_PS_new.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

%% Création de la base de données

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

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

filename_0 = [num2str(patient),complement,'0'];
nom_fichier_inp = append(filename_0,'_PS_new.inp');

for i=1:n_samples
    filename = append(filename_0,num2str(i));
    nom_fichier_out = append(filename,'_PS_new_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),E3);
    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 en 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 en 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 9 et Uy pour 6, 7 et 9
        erreur(i) = ((Y(i,4)-Y_mes(1))^2+(Y(i,5)-Y_mes(2))^2)^0.5; % uniquement les points 6 et 7 en Y
    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 en 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_{RA} (MPa)');
ylabel('E_{LM} (MPa)');
zlabel('Error (mm)');
title('Surface Plot of Error');
file_name = fullfile(folder_name,'E1_E2.fig');
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 = fullfile(folder_name, 'mon_reseau_neurones.mat');
save(file_name,'net');
file_name = fullfile(folder_name, 'Y.mat');
save(file_name,'Y');
file_name = fullfile(folder_name, 'X.mat');
save(file_name,'X');
file_name = fullfile(folder_name, 'erreur.mat');
save(file_name,'erreur');


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

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

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

% Optimisation des paramètres E1, E2 pour minimiser l'écart avec les données expérimentales
% :!\ bien choisir la fonction d'optimisation cohérente avec la mesure
% d'erreur faite précédemment
if patient == 5
    [E_opt, fval] = fmincon(@(E) objective_function_Y1Y2(E, net,Y_mes(1),Y_mes(2)), init_params, [], [], [], [], lb, ub, [], options);
elseif patient ==13
    [E_opt, fval] = fmincon(@(E) objective_function_Y1Y2(E, net,Y_mes(1),Y_mes(2)), init_params, [], [], [], [], lb, ub, [], options);
elseif patient == 18
    [E_opt, fval] = fmincon(@(E) objective_function_Y1Y2(E, net,Y_mes(1),Y_mes(2)), init_params, [], [], [], [], lb, ub, [], options);
elseif patient == 20
    [E_opt, fval] = fmincon(@(E) objective_function_X2Y1Y2Y3(E, net,X_mes(2),Y_mes(1),Y_mes(2),Y_mes(3)), init_params, [], [], [], [], lb, ub, [], options);
end

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


%% Vérification
filename = append(filename_0,'_PS_verif');
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), E3);
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) ;

erreur = ((Y_sim(1)-Y_mes(1))^2+(Y_sim(2)-Y_mes(2))^2)^0.5; % uniquement les points 6 et 7 en Y
fprintf('L erreur associée est e = %f \n', erreur);

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.X(6:10),coord.Y(6:10),'k^')
% hold on
% plot(coord_simu.Var2(:),coord_simu.Var3(:),'m^')
% hold on
% plot(X_mes,Y_mes,'kx')
% hold on
% plot(X_sim,Y_sim,'mx')
% xlim([0 100]);
% ylim([0 100]);
% %tit = ['Patient',num2str(patient),'   E2 = ',num2str(E_opt),' MPa'];
% tit = sprintf('Patient %d   E2 = %.2e MPa',patient,E_opt);
% title(tit);
% legend('coord exp init','coord simu init','coord exp def','coord simu def')
% file_name = fullfile(folder_name, 'comparaison_exp_simu_net_2.fig');
% saveas(gcf,file_name);

if patient == 20
    testx = [coord.X3(7);coord.X3(9);coord.X3(6)];
    testy = [coord.Y3(7);coord.Y3(9);coord.Y3(6)];
elseif patient == 5
    testx = [coord.X3(6);coord.X3(7);coord.X3(9)];
    testy = [coord.Y3(6);coord.Y3(7);coord.Y3(9)];
elseif patient == 13
    testx = [coord.X3(9);coord.X3(7);coord.X3(6)];
    testy = [coord.Y3(9);coord.Y3(7);coord.Y3(6)];
elseif patient == 18
    testx = [coord.X3(7);coord.X3(6);coord.X3(9)];
    testy = [coord.Y3(7);coord.Y3(6);coord.Y3(9)];
end

figure
if patient == 5
    %plot(coord.X3(6:8),coord.Y3(6:8),'r^')
    plot(testx(1:2),testy(1:2),'r^')
    hold on
    plot(coord_simu.Var2(1:2),coord_simu.Var3(1:2),'k^')
    hold on
    plot(X_mes(1:2),Y_mes(1:2),'b^')
    hold on
    plot(X_sim(1:2),Y_sim(1:2),'m^')
else
    plot(testx,testy,'r^')
    hold on
    plot(coord_simu.Var2(:),coord_simu.Var3(:),'k^')
    hold on
    plot(X_mes,Y_mes,'b^')
    hold on
    plot(X_sim,Y_sim,'m^')
end

tit = sprintf('Participant %d  E_{RA} = %.2e MPa E_{LM} = %.2e MPa error = %.2e',patient,E_opt(1), E_opt(2), erreur);
xlim([0 100]);
ylim([0 100]);
title(tit);
xlabel('X (mm)') 
ylabel('Y (mm)') 
legend('init center','init model','def center','def model')
file_name = fullfile(folder_name, 'comparaison_exp_simu_net_2.fig');
saveas(gcf,file_name);

% init exp, init simu, def exp, def simu
trace_final = [testx,testy,coord_simu.Var2(:),coord_simu.Var3(:),X_mes',Y_mes',X_sim',Y_sim']
file_name = fullfile(folder_name, 'trace_final.mat');
save(file_name,'trace_final');

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