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

delete("saves/exercice3.gif");

% Paramètres :
N = 0; % Nombre de disques d'une configuration
R = 10; % Rayon des disques
nb_points_affichage_disque = 30;
increment_angulaire = 2 * pi / nb_points_affichage_disque;
theta = 0:increment_angulaire:2 * pi;
rose = [253 108 158] / 255;
q_max = 1000;
nb_affichages = 20;
pas_entre_affichages = floor(q_max / nb_affichages);
temps_pause = 0.001;

alpha = 0.99;
beta = 1;
lambda = 100;
T = 0.1;
gamma = 5;
S = 130;

a1 = 10;
a2 = 7;
e1 = 0.8;

% Lecture et affichage de l'image :
I = imread('colonie.png');
I = rgb2gray(I);
I = double(I);
I = I(1:400, 100:450);
[nb_lignes, nb_colonnes] = size(I);
poog = figure('Name', ['Detection de ' num2str(N) ' flamants roses'], 'Position', [0.25 * L, 0, 0.4 * L, 0.5 * H]);

pause

% Tirage aléatoire d'une configuration initiale et calcul des niveaux de gris moyens :
c = zeros(N, 2);
a = zeros(N, 1);
e = zeros(N, 1);
teta = zeros(N, 1);
I_moyen = zeros(N, 1);

for i = 1:N
    c(i, :) = [nb_colonnes * rand nb_lignes * rand];
    a(i) = a1 + a2 * rand();
    e(i) = e1 * rand();
    teta(i) = 2 * pi * rand();
    I_moyen(i) = calcul_I_moyen_elli(I, c(i), a(i), e(i), teta(i));
end

liste_q = 0;
liste_U = 0;
liste_N = N;
liste_I_moyen_config = 0;

% Recherche de la configuration optimale :
q = 1;
qq = 1;

while 1
    c_ancien = c;
    a_ancien = a;
    e_ancien = e;
    teta_ancien = teta;

    % 1 Naissances
    N_til = poissrnd(lambda);
    N = length(c);
    c = [c; zeros(N_til, 2)];
    a = [a; zeros(N_til, 1)];
    e = [e; zeros(N_til, 1)];
    teta = [teta; zeros(N_til, 1)];
    I_moyen = [I_moyen; zeros(N_til, 1)];

    for i = N + 1:N + N_til
        c(i, :) = [nb_colonnes * rand nb_lignes * rand];
        a(i) = a1 + a2 * rand();
        e(i) = e1 * rand();
        teta(i) = 2 * pi * rand();
        I_moyen(i) = calcul_I_moyen_elli(I, c(i, :), a(i), e(i), teta(i));
    end

    % 2 Tri des disques
    N = length(c);
    U = 1 - 2 ./ (1 + exp(-gamma * (I_moyen / S - 1)));
    [U, ordre] = sort(U, 1, "descend");
    I_moyen = I_moyen(ordre);
    c = c(ordre, :);
    a = a(ordre);
    e = e(ordre);
    teta = teta(ordre);

    % 3 Morts
    i = 1;
    calcul_U = 1;

    while i < length(c)
        N = length(c);

        if calcul_U
            diffX = c(:, 1) - c(:, 1)';
            diffY = c(:, 2) - c(:, 2)';
            dist = sqrt(diffX.^2 + diffY.^2);
            delta = sum(dist < sqrt(2) * (a1 + a2 / 2), "all") - N;
            U_tot = sum(U) + beta * delta;
        end

        diffX = c([1:i - 1 i + 1:N], 1) - c([1:i - 1 i + 1:N], 1)';
        diffY = c([1:i - 1 i + 1:N], 2) - c([1:i - 1 i + 1:N], 2)';
        dist = sqrt(diffX.^2 + diffY.^2);
        delta = sum(dist < sqrt(2) * (a1 + a2 / 2), "all") - (N - 1);
        U_c = sum(U([1:i - 1 i + 1:N])) + beta * delta;

        proba = lambda / (lambda + exp((U_c - U_tot) / T));

        if proba > rand()
            c(i, :) = [];
            a(i) = [];
            e(i) = [];
            teta(i) = [];
            I_moyen(i) = [];
            U(i) = [];
            calcul_U = 1;
        else
            i = i + 1;
            calcul_U = 0;
        end

    end

    N = length(c);

    % 4 Test de convergence
    if q > q_max
        break
    end

    T = alpha * T;
    lambda = alpha * lambda;

    % Si le disque proposé est "meilleur", mises à jour :
    hold off;
    subplot(2, 2, 1);
    imagesc(I);
    axis image;
    axis off;
    colormap gray;
    hold on;

    for j = 1:N
        AZE = a(j) * [cos(theta); sqrt(1 - e(j)^2) * sin(theta)];
        AZE = c(j, :)' + [cos(teta(j)) (-sin(teta(j))); sin(teta(j)) cos(teta(j))] * AZE;
        x_affich = AZE(1, :);
        y_affich = AZE(2, :);
        indices = find(x_affich > 0 & x_affich < nb_colonnes & y_affich > 0 & y_affich < nb_lignes);
        subplot(2, 2, 1);
        plot(x_affich(indices), y_affich(indices), 'Color', rose, 'LineWidth', 3);
    end

    pause(temps_pause);

    % Courbe d'évolution du niveau d'énergie :
    diffX = c(:, 1) - c(:, 1)';
    diffY = c(:, 2) - c(:, 2)';
    dist = sqrt(diffX.^2 + diffY.^2);
    delta = sum(dist < sqrt(2) * (a1 + a2 / 2), "all") - N;
    U_tot = sum(U) + beta * delta;

    if rem(q, pas_entre_affichages) == 0 || q == 1
        liste_q = [liste_q q];
        liste_U = [liste_U U_tot];
        liste_N = [liste_N N];

        I_moyen_config = mean(I_moyen);
        liste_I_moyen_config = [liste_I_moyen_config I_moyen_config];
        subplot(2, 2, 2);
        plot(liste_q, liste_I_moyen_config, '.-', 'Color', rose, 'LineWidth', 3);
        axis([0 q_max 0 255]);
        xlabel('Nombre d''iterations');
        ylabel('Niveau de gris moyen');

        subplot(2, 2, 3);
        plot(liste_q, liste_U, '.-', 'Color', rose, 'LineWidth', 3);
        xlabel('Nombre d''iterations');
        ylabel('Energie');
        xlim([0, q_max]);

        subplot(2, 2, 4);
        plot(liste_q, liste_N, '.-', 'Color', rose, 'LineWidth', 3);
        xlim([0, q_max]);
        xlabel('Nombre de''iteration');
        ylabel('Nombre de point');

        export_fig(poog, "saves/exercice3.gif", '-png', '-painters', '-m2', '-append');
    end

    q = q + 1;
end