% FKPP - version individu-centre 1D remplacement aleatoire
% Des individus sont distribues sur un intervalle
% Certains mutants sont porteurs d'un gene favorable qui donne
% un avantage de croissance r
% Tous les individus de deplacent selon une marche aleatoire avec
% des sauts en sqrt(2*D)
% Les wild-types (non-porteurs) on un taux de croissance 0
% A chaque reproduction d'un mutant, l'individu le plus
% proche est remplace par la progeniture
% Le nombre total d'individu est constant = N0
% la croissance est limitee par le fait la probabilite de remplacer un
% wild-type est de 1 - umut, ou umut est
% la densite de population de mutant
% L'avantage de ce code est qu'il ne depend pas de la densite umut
% umut n'est calculee que pour l'affichage
%
% prob(mutant en x se reproduise) = dt*r
% prob(remplacer un wild-type) ~ (1-umut(x))

%% parametre des equations
r = 0.5;        % taux de croissance
D = 0.1;        % coefficient de diffusion

%% parametres de simulation, espace
% x in [x0,x1]
h = 0.1;        % intervalle de discretisation spatiale
S = 20.0;       % longueur de l'intervalle
x0 = 0;         % bord gauche de l'intervalle
x1 = S;         % bord droit de l'intervalle
xi = x0:h:x1;   % discretisation
J = length(xi); % nombre de points de discretisation

%% Individu-centre
N0 = 10000;   % population totale

% chaque individu possede une position x et trait g
% g = 1 : gene favorable (mutants)
% g = 0 : pas de gene favorable
x = x0 + S*rand([1,N0]); % individus distribues sur [x0,x1]
g = false(1,N0);
g( x < S/10 ) = true; % individus sur la gauche de l'intervalle sont porteurs
w = zeros(1,N0);     % probabilite pour l'individu de se reproduire

%% parametres de simulation, temps
t0 = 0;
tfinal = 30.0;
t = t0;
dt = 0.1

%% Difference finie  pour comparaison
% methode implicite Crank-Nicolson
ufd = zeros(J,1); % stocke seulement l'etat au temps t
newufd = zeros(J,1);
% Le Laplacien discretise est une matrice L de taille JxJ
% En 1D elle est symetrique et tridiagonale
L0 = sparse(1:J,1:J,-2); % matrice creuse, compacte en memoire
L0 = spdiags(ones(J,2),[-1 1],L0); % forme la matrice tridiagonale
L0(1,:) = 0; % u(x0) est donne par les conditions au bord
L0(J,:) = 0; % u(x1) est donne par les conditions au bord
% condition initiales
ufd = zeros(J,1);
ufd( xi < S/10 ) = 1.0;
% Schema implicite Crank-Nicolson
A = (speye(J) - dt/h^2*D/2*L0);

%% Initialisation
figure(1); clf;
umut = ksdensity(x(g),xi,'support',[x0,x1]); % kernel density approximation des mutants
plot(xi,sum(g)*umut*S/N0,xi,ufd)

%% BOUCLE PRINCIPALE
tic
while t < tfinal
    drawnow;
    % on calcule pour chaque mutant une probabilite de se reproduire
    % (les individus non porteurs ne se reproduisent pas)
    w  = dt*r*g;                                   % probabilite de se reproduire
    irep = find(rand(1,N0) < w);                 % realisation
    % un individu se reproduisant remplace l'individu le plus proche
    % (autre que lui-meme)
    % la probabilite de choisir un wild-type est environ 1-umut
    for i = 1:length(irep)
        [~,rempl] = min(abs(x - x(irep(i))) + 10*(x == x(irep(i))));
        g(rempl) = true;
    end

    % on deplace les individu en marche aleatoire
    x = x + sqrt(dt)*sqrt(2*D)*randn(size(x));
    x = abs(x); % condition en x0 reflechissant

    umut = ksdensity(x(g),xi,'support',[x0,x1]); % kernel density approximation des mutants

    % solution differences finies implicite pour comparaison
    newufd =A\(ufd + dt*r*ufd.*(1-ufd) + dt/h^2*D/2*L0*ufd);
    newufd(1) = newufd(2);
    newufd(J) = newufd(J-1);
    ufd = newufd;

    % affichage
    plot(xi,sum(g)*umut*S/N0,xi,ufd);
    axis([x0 x1 0 1.1]);
    t = t + dt;
    fprintf("t = %.5f, Nmut = %d\n",t, sum(g));
end
toc