clear();

//========================================================================================================================
//Ce programme modélise une statistique d'épidémie
//========================================================================================================================

//========================================================================================================================
//Saisie des paramètres
//      - N: Nombre d'individus
//      - T: Durée de la simulation
//      - n: Nombre de subdivisions en temps
//      - sigma: Amplitude des déplacements des individus
//      - eps: Paramètre de distanciation sociale
//      - TR: Temps de guérison (en unité de temps: T/n)
//========================================================================================================================

N = 500;
T = 20;
n = 50;
sigma = 0.01;
eps = 0.05;
TR = 10;

//========================================================================================================================
//Introduction d'autres paramètres
//      - tau: Pas de subdivision en temps
//      - XX, YY: Abscisses et ordonnées des points de tous les individus à un instant donné (vecteurs ligne)
//      - X,Y: Matrices donnant les coordonnées de tous les points à tous les instants
//========================================================================================================================

tau = T/n; //Pas de subdivision en temps

XX = grand(1,N,'unf',0,1);
YY = grand(1,N,'unf',0,1);

XX = XX - floor(XX); //Pour se ramener au segment [0,1[
YY = YY - floor(YY); //Pour se ramener au segment [0,1[

X = XX;
Y = YY;

//========================================================================================================================
//Création de matrices
//      - EX, EY: Matrices donnant les écarts entre abscisses et ordonnées des individus deux par deux à tout instant
//      - Dist: Matrice donnant la distance entre les individus deux par deux (notons que c'est une matrice symétrique)
//      - Etatd: Vecteur colonne donnant l'état de chaque individu à un instant donné: 0 (malade) ou 1 (non malade)
//      - Etat: Matrice donnant l'état de chaque individu à tous les instants
//      - Ret: Matrice donnant l'immunité de chaque individu à l'infection: 0 (immunisé) ou 1 (non immunisé). Noter qu'
//après infection, un individu devient immunisé. Tout individu est considéré comme contagieux durant l'infection.
//      - I: Matrice donnant la proportion d'infectés à tous les instants.
//========================================================================================================================

EX = X(1,:)'*ones(1,N) - ones(N,1)*X(1,:);
EY = Y(1,:)'*ones(1,N) - ones(N,1)*Y(1,:);

Dist = sqrt(EX.^2 + EY.^2);

Etat = [1;zeros(N-1,1)];
Etatd = Etat;

Ret = ones(N,n+TR);

I = []; //Proportion d'infectés

//========================================================================================================================
//Simulation de l'épidémie
//========================================================================================================================

for j=1:n
    XX = XX+sigma*grand(1,N,'nor',0,tau);
    YY = YY+sigma*grand(1,N,'nor',0,tau);
    XX = XX - floor(XX);
    YY = YY - floor(YY);
    X = [X;XX];
    Y = [Y;YY];
    EX = X(j,:)'*ones(1,N) - ones(N,1)*X(j,:); //Matrice des écarts en abscisse
    EY = Y(j,:)'*ones(1,N) - ones(N,1)*Y(j,:); //Matrice des écarts en ordonnée
    Dist = [Dist,sqrt(EX.^2 + EY.^2)]; //Matrice des distances entre les individus
    Ret(:,j+TR) = floor(0.5*sign(Ret(:,j) - Etatd)+0.5);
    Etatd = sign(round(sign(Dist(:,(j-1)*N+1:j*N)-clean(Dist(:,(j-1)*N+1:j*N),eps)))*(Ret(:,j).*Etatd) + Ret(:,j+1).*Etatd);
    Etat = [Etat,Ret(:,j+1).*Etatd];
    //disp(sign(Dist(:,(j-1)*N+1:j*N)-clean(Dist(:,(j-1)*N+1:j*N),eps)));
    //disp(Etat);
    //disp(Ret);
    I = [I,sum(Etat(:,j))/N];
    disp('Calculs: '+string(0.1*round(1000*j/n))+'%');
end

Ret = floor(Ret);
Etat = floor(Etat);

//========================================================================================================================
//Introduction d'un nouveau paramètre
//      - Rep: Vecteur colonne donnant l'indice de reproduction de l' infection (appelé R effectif) à tous les instants
//========================================================================================================================

Rep = I(2:$)./(I(1:$-1)+1E-10); //Indice de reproduction du virus
Rep = [Rep,Rep($)];
Rep = 0.01*round(100*Rep);

//========================================================================================================================
//Affichage de la simulation
//========================================================================================================================

f0 = scf(1);
xclick();   
delete();
aa=get("current_axes")
aa.data_bounds=[min(X),min(Y);max(X),max(Y)];
j=1;
drawlater();
    b=get("current_entity");
    delete(b.children);
    subplot(2,2,1);
    p=scatter(X(j,:),Y(j,:),100,[Etat(:,j),0.7*(ones(N,1)-Ret(:,j)).*(ones(N,1)-Etat(:,j))+Ret(:,j).*(ones(N,1)-Etat(:,j)),Ret(:,j).*(ones(N,1)-Etat(:,j))],"fill");
    p.mark_foreground = color("black");
    title('Evolution de la population à t = '+string(round(100*(j-1)*tau)/100)+', Sains: '+string(100-0.1*round(1000*sum(Etat(:,j))/N))+'%'+', Infectés: '+string(0.1*round(1000*sum(Etat(:,j))/N))+'%'+', Rétablis: '+string(0)+'%');
    
    subplot(2,2,2);
    set(gca(),"grid",[1 1]);
    aa=get("current_axes")
    aa.data_bounds=[0,0;T,100];
    xlabel('Temps');
    ylabel('Infectés (pourcentage)');
    title('Evolution du pourcentage d infectés à t= '+string(round(100*(j-1)*tau)/100));
    plot(tau,100*I(1),'color','red','thickness',3);
    
    subplot(2,2,3);
    set(gca(),"grid",[1 1]);
    aa=get("current_axes")
    aa.data_bounds=[0,0;T,max(Rep)];
    xlabel('Temps');
    ylabel('R effectif');
    title('Evolution de l indice de reproduction à t= '+string(round(100*(j-1)*tau)/100));
    plot(tau,1,'color','orange','thickness',3);
    
drawnow();
sleep(50);
for j=2:n+1
    drawlater();
        b=get("current_figure");
        delete(b.children);
        subplot(2,2,1);
        aa=get("current_axes")
        aa.data_bounds=[min(X),min(Y);max(X),max(Y)];
        p=scatter(X(j,:),Y(j,:),100,[Etat(:,j),0.7*(ones(N,1)-Ret(:,j)).*(ones(N,1)-Etat(:,j))+Ret(:,j).*(ones(N,1)-Etat(:,j)),Ret(:,j).*(ones(N,1)-Etat(:,j))],"fill");
        p.mark_foreground = color("black");
        title('Evolution de la population à t = '+string(round(100*(j-1)*tau)/100)+', Sains: '+string(0.1*round(1000*sum(Ret(:,j))/N)-0.1*round(1000*sum(Etat(:,j))/N))+'%'+', Infectés: '+string(0.1*round(1000*sum(Etat(:,j))/N))+'%'+', Rétablis: '+string(100-0.1*round(1000*sum(Ret(:,j))/N))+'%');
    
        subplot(2,2,2);
        set(gca(),"grid",[1 1]);
        aa=get("current_axes")
        aa.data_bounds=[0,0;T,100];
        xlabel('Temps');
        ylabel('Infectés (pourcentage)');
        title('Evolution du pourcentage d infectés à t= '+string(round(100*(j-1)*tau)/100));
        plot((1:j-1)*tau,100*I(1:j-1),'color','red','thickness',3);
        
        subplot(2,2,3);
        set(gca(),"grid",[1 1]);
        aa=get("current_axes")
        aa.data_bounds=[0,0;T,max(Rep)];
        xlabel('Temps');
        ylabel('R effectif');
        title('Evolution de l indice de reproduction à t= '+string(round(100*(j-1)*tau)/100));
        plot((1:j-1)*tau,Rep(1:j-1),'color','orange','thickness',3);
        legend('Taux de reproduction: R='+string(round(100*Rep(j-1))/100));
    drawnow();
    //xclick();
    sleep(50);
    
end
