diff --git a/LICENSE.md b/LICENSE.md
new file mode 100644
index 0000000000000000000000000000000000000000..4ee22f49fb4d7f9a461564fa2bcb292d60212fd2
--- /dev/null
+++ b/LICENSE.md
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (©) 2012-2022 INRIA
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/matlab/circ_vmpdf.m b/matlab/circ_vmpdf.m
new file mode 100644
index 0000000000000000000000000000000000000000..52bc381112461ac5937426e0df5c69f71fa6e32b
--- /dev/null
+++ b/matlab/circ_vmpdf.m
@@ -0,0 +1,60 @@
+function [p alpha] = circ_vmpdf(alpha, thetahat, kappa)
+
+% [p alpha] = circ_vmpdf(alpha, w, p)
+% Computes the circular von Mises pdf with preferred direction thetahat
+% and concentration kappa at each of the angles in alpha
+%
+% The vmpdf is given by f(phi) =
+% (1/(2pi*I0(kappa))*exp(kappa*cos(phi-thetahat)
+%
+% Input:
+% alpha angles to evaluate pdf at, if empty alphas are chosen to
+% 100 uniformly spaced points around the circle
+% [thetahat preferred direction, default is 0]
+% [kappa concentration parameter, default is 1]
+%
+% Output:
+% p von Mises pdf evaluated at alpha
+% alpha angles at which pdf was evaluated
+%
+%
+% References:
+% Statistical analysis of circular data, Fisher
+%
+% Circular Statistics Toolbox for Matlab
+
+% By Philipp Berens and Marc J. Velasco, 2009
+% velasco@ccs.fau.edu
+
+% if no angles are supplied, 100 evenly spaced points around the circle are
+% chosen
+if nargin < 1 || isempty(alpha)
+ alpha = linspace(0, 2*pi, 101)';
+ alpha = alpha(1:end-1);
+end
+if nargin < 3
+ kappa = 1;
+end
+if nargin < 2
+ thetahat = 0;
+end
+
+alpha = alpha(:);
+
+% evaluate pdf
+C = 1/(2*pi*besseli(0,kappa));
+p = C * exp(kappa*cos(alpha-thetahat));
+
+
+% Toolbox for circular statistics with Matlab.
+%
+% Authors: Philipp Berens
+% Email: philipp@bethgelab.org
+% Homepage: http://philippberens.wordpress.com/code/circstats/
+%
+% Contributors:
+% Marc Velasco, Tal Krasovsky
+%
+% Reference:
+% P. Berens, CircStat: A Matlab Toolbox for Circular Statistics, Journal of Statistical Software, Volume 31, Issue 10, 2009
+% http://www.jstatsoft.org/v31/i10
diff --git a/matlab/circ_vmrnd.m b/matlab/circ_vmrnd.m
new file mode 100644
index 0000000000000000000000000000000000000000..9a9d596fab706666b9922b4ab25c9d4e1d11bb76
--- /dev/null
+++ b/matlab/circ_vmrnd.m
@@ -0,0 +1,99 @@
+function alpha = circ_vmrnd(theta, kappa, n)
+
+% alpha = circ_vmrnd(theta, kappa, n)
+% Simulates n random angles from a von Mises distribution, with preferred
+% direction thetahat and concentration parameter kappa.
+%
+% Input:
+% [theta preferred direction, default is 0]
+% [kappa width, default is 1]
+% [n number of samples, default is 10]
+%
+% If n is a vector with two entries (e.g. [2 10]), the function creates
+% a matrix output with the respective dimensionality.
+%
+% Output:
+% alpha samples from von Mises distribution
+%
+%
+% References:
+% Statistical analysis of circular data, Fisher, sec. 3.3.6, p. 49
+%
+% Circular Statistics Toolbox for Matlab
+
+% By Philipp Berens and Marc J. Velasco, 2009
+% velasco@ccs.fau.edu
+
+
+% default parameter
+if nargin < 3
+ n = 10;
+end
+
+if nargin < 2
+ kappa = 1;
+end
+
+if nargin < 1
+ theta = 0;
+end
+
+if numel(n) > 2
+ error('n must be a scalar or two-entry vector!')
+elseif numel(n) == 2
+ m = n;
+ n = n(1) * n(2);
+end
+
+% if kappa is small, treat as uniform distribution
+if kappa < 1e-6
+ alpha = 2*pi*rand(n,1);
+ return
+end
+
+% other cases
+a = 1 + sqrt((1+4*kappa.^2));
+b = (a - sqrt(2*a))/(2*kappa);
+r = (1 + b^2)/(2*b);
+
+alpha = zeros(n,1);
+for j = 1:n
+ while true
+ u = rand(3,1);
+
+ z = cos(pi*u(1));
+ f = (1+r*z)/(r+z);
+ c = kappa*(r-f);
+
+ if u(2) < c * (2-c) || ~(log(c)-log(u(2)) + 1 -c < 0)
+ break
+ end
+
+
+ end
+
+ alpha(j) = theta + sign(u(3) - 0.5) * acos(f);
+ alpha(j) = angle(exp(i*alpha(j)));
+end
+
+if exist('m','var')
+ alpha = reshape(alpha,m(1),m(2));
+end
+
+% Toolbox for circular statistics with Matlab.
+%
+% Authors: Philipp Berens
+% Email: philipp@bethgelab.org
+% Homepage: http://philippberens.wordpress.com/code/circstats/
+%
+% Contributors:
+% Marc Velasco, Tal Krasovsky
+%
+% Reference:
+% P. Berens, CircStat: A Matlab Toolbox for Circular Statistics, Journal of Statistical Software, Volume 31, Issue 10, 2009
+% http://www.jstatsoft.org/v31/i10
+
+
+
+
+
diff --git a/matlab/main_clonal.m b/matlab/main_clonal.m
new file mode 100644
index 0000000000000000000000000000000000000000..7c3d14823369d29a0b7e6d6879f5518b5cc0581e
--- /dev/null
+++ b/matlab/main_clonal.m
@@ -0,0 +1,965 @@
+function main_clonal(run_mode_choice)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% main routine:
+% run_mode = 0 : silent mod / all result in a new directory
+% run_mode = 1 : verbose mode without graphics (no new directory created)
+% run_mode = 2 : verbose mode with graphics (no new directory created)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% content : simulation of a clonal plant growth on a given resource
+% map resources(x,y) this map is not affected by the IBM
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Main routine :
+% main_clonal
+%
+% Authors :
+% Fabien Campillo and Nicolas Champagnat (INRIA)
+% Fabien.Campillo@inria.fr
+% http://www-sop.inria.fr/members/Fabien.Campillo/
+%
+% Reference :
+% F. Campillo and N. Champagnat. "Simulation and analysis of an
+% individual-based model for graph-structured plant dynamics",
+% Ecological Modelling, Volume 234, 10 June 2012, Pages 93-105
+% http://www.sciencedirect.com/science/article/pii/S030438001200124X
+%
+% see
+% http://www-sop.inria.fr/members/Fabien.Campillo/software/ibm-clonal/
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% The software was developed within the ANR Syscomm (SYSt�mes COmplexes et
+% Mod�lisation Math�matique) project MODECOL (MOD�lisation ECOLogique
+% de prairies virtuelles) [ANR-08-SYSC-012].
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% the main routine main_clonal uses two routines :
+% circ_vmrnd : Simulates n random angles from a von Mises
+% distribution, with preferred direction thetahat and
+% concentration parameter kappa
+% circ_vmpdf : Computes the circular von Mises pdf with preferred
+% direction thetahat and concentration kappa at each of
+% the angles in alpha
+% from the Toolbox for circular statistics with Matlab:
+%
+% Authors: Philipp Berens
+% Email: philipp@bethgelab.org
+% Homepage: http://philippberens.wordpress.com/code/circstats/
+% Contributors: Marc Velasco, Tal Krasovsky
+% Reference: P. Berens, CircStat: A Matlab Toolbox for Circular
+% Statistics, Journal of Statistical Software, Volume 31, Issue 10,
+% 2009. http://www.jstatsoft.org/v31/i10
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% First version May 31, 2011
+% Last version August 20, 2012
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+path('./src/',path)
+
+%% --- data structures -----------------------------------------------------
+
+% >>>> INPUTS
+
+% run_mode : 0 or '' silent mode .... (do run(999))
+% muLN s2LN : log-normal lengtn dist. of the shoots
+% muVM kaVM : Von Mises angle dist. of the shoots (wrt the gradient)
+% instants : instants of event
+% LL : side size of the field
+% LLzoom : zoomed area for output
+% iseed : seed of the generator
+% n_events : maximal nb of events
+% plot_event : number of event between 2 outpus
+% lambda0 alpha : parameters of the birth rate
+% mu0 beta : parameters of the death rate
+% initial_ramet : position of the initial ramet
+% N_ini : # initial ramets
+% diffusion : 0/1 if diffusion or not
+% sigma_dif : diffusion cefficient
+
+% cons_std_dev : consomption parameters std deviation
+% cons_max : maximum consomption
+% cons_th : minimim threshold
+
+% >>>> GRAPHICS
+
+% node_size : size of the node
+% node_color : color of the node
+% link_color : color of the link
+% type_graph : 1 node, 2 links, 2 nodes+links
+% contourlevels : number of contour level [useless]
+
+% >>>> OUTPUTS
+
+% population : structure of the population at a given time
+% result : instants and nodes at each times of output
+% resources : array of available resources
+% x_grid y_grid : grid of the available resources
+% pop_sizes : evolution of the size of pop
+% pop_fitness : evolution of the fiteness (res(x_i)/N)
+% total_res : evolution of the available resources
+% total_event : nb of realized events (neq n_events in case of extinct)
+
+% fid : id of the output file
+
+%% --- global settings ----------------------------------------------------
+
+global run_mode population resources resources_max result x_grid y_grid
+global muLN s2LN muVM kaVM lambda0 alpha mu0 beta N_ini
+global LL LLzoom gridsize contourlevels iseed n_events plot_event
+global instants pop_sizes pop_fitness total_res total_event initial_ramet
+global work_dir fid
+global cons_std_dev cons_max cons_th
+global node_size node_color link_color type_graph
+global diffusion sigma_dif
+
+%% --- structures ---------------------------------------------------------
+
+% --- population structure at time t
+
+population = struct(... % the plant container:
+ 'pos',{},... % ramets' position x,y
+ 'resou',{},... % resource at x,y
+ 'links',{}); % connected ramets indices
+
+% --- result structure
+
+result = struct(... % the result container:
+ 'time',{},... % time
+ 'pos',{},... % nodes' positions at the given time
+ 'res',{}); % resources at the given time
+
+%% --- parameters ---------------------------------------------------------
+
+LL = 10; % field size, hyp: square field [0,LL]^2
+LLzoom = [0 LL 0 LL];%[2.5 6 6 9.5];%[3 6 6 9];
+iseed = 12;%1111;
+
+
+n_events = 100000;%100000; % number of events
+plot_event = 1; %round(n_events/200);%1000; % events between plots
+
+
+% --- coef for the individual birth rates :
+% lambda0+alpha*r(x) [r(x) resources in x]
+lambda0 = 2; %1
+alpha = 300; %10
+
+% --- coef for the individual death rates :
+% mu0+beta*[rmax-r(x)] [r(x) resources in x]
+mu0 = 2;
+beta = 100;
+
+% --- consomption parameters
+cons_std_dev = 0.03;
+cons_max = 0.5;
+cons_th = 0.0001;
+
+% --- parameters of the dispersion law
+% ----- lognormal for the length of the link
+% muLN parameter mu of the lognormal
+% s2LN parameter sigma2 of the lognormal
+% ----- von mises distribution for the direction ofthe link
+% muVM parameter mu of the von mises distribution
+% kaVM parameter kappa of the von mises distribution (inverse variance)
+% ----- example of parameters for PHALANX vs GUERILLA behavior
+% muLN kaVM
+% PHALANX -2.5 0.0000000001
+% GUERILLA -1.0 100
+%
+
+muLN = -2.5;
+s2LN = 0.1 ;
+muVM = 0;
+kaVM = 0.0000000001;
+
+% --- diffusion
+diffusion = 1;
+sigma_dif = 0.3;
+
+% --- initial ramet
+initial_ramet = [2 4];
+N_ini=10;
+
+% --- contour plot parameters for the resources
+gridsize = 200; % grid size of the field
+%contourlevels = 20 ; % nb of contour levels (useless)
+[x_grid, y_grid] = meshgrid(0:LL/gridsize:LL);
+
+% --- graphics
+node_size = 1;
+node_color = 'r';
+link_color = 'g';
+type_graph = 3;
+
+%% --- runs ---------------------------------------------------------------
+
+if nargin == 0, run_mode_choice = 0; end
+run_mode = run_mode_choice;
+
+switch run_mode
+
+ case 0 % --- silent mode, simulation, no graphics
+
+ work_dir = [pwd '/essai_' datestr(now,30)];
+ unix(['mkdir ' work_dir]);
+ unix(['cp ' pwd '/main_clonal.m ' work_dir]);
+ fid = fopen([work_dir '/out.txt'],'w');
+ simu();
+ fclose(fid);
+ save([work_dir '/data.mat'])
+
+ case 1 % --- verbose mode, simulation, no graphics
+
+ disp('-------------------------------------------------------------')
+ disp('--- simulation + no graphics --------------------------------')
+ disp(' ')
+ fid = 1 ; % std output
+ simu();
+
+ case 2 % --- verbose mode, simulation, graphics (on-line animation)
+
+ disp('-------------------------------------------------------------')
+ disp('--- simulation + on-line animation --------------------------')
+ disp(' ')
+ fid = 1 ; % std output
+ simu();
+
+ case 3 % --- verbose, no simulation, interface for movies creation
+
+ disp('-------------------------------------------------------------')
+ disp('--- no simulation + movie creation --------------------------')
+ disp(' ')
+ liste_essais = dir('essai*');
+ nb_essais = size(liste_essais,1);
+ disp(' ')
+ disp(' >> list of tests :')
+ for i=1:nb_essais
+ disp([' ' int2str(i) ' ' liste_essais(i).name])
+ end
+ answer=str2double(input(' choice: ','s'));
+ work_dir = [pwd '/' liste_essais(answer).name];
+ fprintf(' loading the data... ')
+ load([work_dir '/data.mat'])
+ fprintf('done \n')
+ work_dir = [pwd '/' liste_essais(answer).name]; % modified by load
+ disp(' creating the video: ')
+ create_video()
+
+ case 4 % --- verbose, no simulation, pre_graphics
+
+ disp('-------------------------------------------------------------')
+ disp('--- no simulation + pre-graphics (no sim needed) ------------')
+ disp(' ')
+
+ pre_graphics()
+
+ case 5 % --- verbose, no simulation, post_graphics
+
+ disp('-------------------------------------------------------------')
+ disp('--- no simulation + post-graphics (sim needed) --------------')
+ disp(' + some snaptshots ')
+ disp(' ')
+ liste_essais = dir('essai*');
+ nb_essais = size(liste_essais,1);
+ disp(' ')
+ disp(' >> list of tests :')
+ for i=1:nb_essais
+ disp([' ' int2str(i) ' ' liste_essais(i).name])
+ end
+ answer=str2double(input(' choice: ','s'));
+ work_dir = [pwd '/' liste_essais(answer).name];
+ fprintf(' loading the data... ')
+ load([work_dir '/data.mat'])
+ fprintf('done \n')
+ work_dir = [pwd '/' liste_essais(answer).name]; % modified by load
+ some_snapshots()
+ full_graphics()
+
+ otherwise
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ disp(' ')
+ disp(' run : ')
+ disp(' sim basic ani movie verbose')
+ disp(' graphics mation creation ')
+ disp(' ------------------------------------------------------------------')
+ disp(' 0 x x - - - (for batch: see matbg)')
+ disp(' 1 x x - - x')
+ disp(' 2 x x x - x')
+ disp(' 3 - - x x x (with interface)')
+ disp(' 4 - - - - x (pre graphics)')
+ disp(' 5 - x - - x (post graphics)')
+ disp(' otherwise this message')
+ disp(' ')
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+end
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% sub functions
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+%%
+function simu()
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% simulation of a clonal plant on a given resource map resources(x,y)
+% this map is not affected by the IBM
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% can be improved: memory preallocation for structures
+% (aviod structure growth phenomenon)
+% TODO for
+% population : population structure at a given time
+% result : structure along time (nodes states at each times)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+global run_mode population resources resources_max result x_grid y_grid
+global lambda0 alpha mu0 beta N_ini
+global LL gridsize iseed n_events total_event plot_event
+global instants pop_sizes pop_fitness total_res
+global fid diffusion sigma_dif
+
+% --- setting the random stream and seeding it (only for R1020)
+
+%mtstream = RandStream('mt19937ar','Seed',iseed);
+%RandStream.setDefaultStream(mtstream);
+
+% --- for older versions:
+
+randn('state',iseed)
+rand('state',iseed)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% --- resources map -------------------------------------------------------
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+resources = 3+abs(peaks(7*(x_grid-0.5*LL)/LL,7*(y_grid-0.5*LL)/LL));
+resources = min(resources,5);
+% resources = 9*0.5*(1+cos(5*x_grid))...
+% .*min(x_grid,1).*min(LL-x_grid,1)...
+% .*min(y_grid,1).*min(LL-y_grid,1); % alternae resource map
+resources_max = max(max(resources));
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% --- dynamics ------------------------------------------------------------
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+tic;
+ini_population();
+
+if run_mode==2
+ plot_landscape()
+ set(gcf,'DoubleBuffer','on')
+end
+
+t = 0;
+i_save = 1;
+step = LL/gridsize;
+
+result(i_save).time = t;
+result(i_save).pos = reshape([population.pos],2,N_ini);
+result(i_save).res = resources;
+
+instants = zeros(n_events+1,1);
+pop_sizes = zeros(n_events+1,1);
+pop_fitness = zeros(n_events+1,1);
+total_res = zeros(n_events+1,1);
+instants(1) = t;
+pop_sizes(1) = size(population,2);
+pop_fitness(1) = sum([population.resou])/size(population,2);
+total_res(1) = sum(sum(resources))*LL*LL/gridsize/gridsize;
+
+
+fprintf(fid,['\n\n'...
+ ' -------------------------- \n'...
+ ' clonal test / ' date() '\n']);
+
+nb_birth = 0;
+nb_death = 0;
+
+
+
+for i_event=1:n_events
+
+ % --- birth/death rates
+
+ lambda = lambda0+alpha*[population.resou]; % birth rates
+ %lambda = lambda0+alpha*exp([population.resou]); %FUN % birth rates
+ mu = mu0+beta*(resources_max-[population.resou]); % death rates
+ llambda = sum(lambda);
+ mmu = sum(mu);
+
+ delta = exprnd(1/(llambda+mmu));
+ t = t+delta; % time increase
+ instants(i_event+1) = t;
+
+ if diffusion ==1
+ % finite difference scheme
+ criteria = sigma_dif*delta/(step*step);
+ if criteria > 0.1
+ disp([' >>> critera ' num2str(sigma_dif*delta/(step*step))]);
+ end
+ Laplace = del2(resources,step);
+ resources = resources + sigma_dif*delta*Laplace;
+ end
+
+ if (rand()<llambda/(llambda+mmu))
+
+ % --- birth -------------------------------------------------------
+
+ nb_birth = nb_birth+1;
+ chosen = f_finiternd(lambda,1); % will give birth
+ add_individual(chosen) % add a daughter node to the node "chosen"
+
+ else
+
+ % --- death -------------------------------------------------------
+
+ nb_death = nb_death+1;
+ chosen = f_finiternd(mu,1); % choose an individual
+ remove_individual(chosen) % removed from population
+
+ end
+
+ pop_sizes(i_event+1) = size(population,2);
+ pop_fitness(i_event+1) = sum([population.resou])/size(population,2);
+ total_res(i_event+1) = sum(sum(resources))*LL*LL/gridsize/gridsize;
+
+ if isempty(population)
+ fprintf('population extinction \n');
+ break
+ end
+
+
+ % --- plot
+
+ if (mod(i_event,plot_event)==0)
+ i_save = i_save+1;
+ cc = (LL/gridsize)*(LL/gridsize);
+ fprintf(fid,' #ite %i of %i / #pop %i / t %g / res %g \n',...
+ i_event,n_events,pop_sizes(i_event),t,sum(sum(resources))*cc);
+
+ result(i_save).time = t;
+ result(i_save).pos = reshape([population.pos],2,[]);
+ result(i_save).res = resources;
+
+ if run_mode==2
+ plot_landscape()
+ drawnow
+ end
+ end
+
+end
+
+total_event = i_event+1;
+
+fprintf(fid,'\n');
+fprintf(fid,[' CPU time ' num2str(toc) ' sec\n']);
+fprintf(fid,' final time %g \n',t);
+fprintf(fid,' birth %i / death %i \n',nb_birth,nb_death);
+fprintf(fid,'\n');
+
+
+
+
+
+function ini_population()
+%--------------------------------------------------------------------------
+% initialise the structure population (1 center node + 3 daughter nodes)
+%--------------------------------------------------------------------------
+global population N_ini LL muLN s2LN initial_ramet
+
+population(1).pos = initial_ramet;
+population(1).resou = f_resources(population(1).pos);
+
+N_ini=4;
+
+for i=2:4
+ angle = 2*pi*rand();
+ length = random('logn',muLN,s2LN);
+ population(i).pos = population(1).pos+length*[cos(angle) sin(angle)];
+ population(i).resou = f_resources(population(i).pos);
+ population(i).links = 1;
+ population(1).links = [population(1).links i];
+end
+
+% for i=1:N_ini
+% population(i).pos = LL*rand(1,2);
+% population(i).resou = f_resources(population(i).pos);
+% population(i).links = [];
+% end
+
+function ini_population2()
+%--------------------------------------------------------------------------
+% initialise the structure population (3 nodes, 2 links)
+%--------------------------------------------------------------------------
+global population muLN s2LN muVM kaVM
+
+population(1).pos = [2 2];
+population(1).resou = f_resources(population(1).pos);
+population(1).links = 2;
+
+angle = circ_vmrnd(muVM,kaVM,1);
+length = random('logn',muLN,s2LN);
+population(2).pos = population(1).pos+length*[cos(angle) sin(angle)];
+population(2).resou = f_resources(population(2).pos);
+population(2).links = [1 3];
+
+angle = circ_vmrnd(muVM,kaVM,1)+pi/2;
+length = random('logn',muLN,s2LN);
+population(3).pos = population(2).pos+length*[cos(angle) sin(angle)];
+population(3).resou = f_resources(population(3).pos);
+population(3).links = 2;
+population(3).links = 2;
+
+
+function remove_individual(chosen)
+%--------------------------------------------------------------------------
+% remove and individual from population
+%--------------------------------------------------------------------------
+
+global population
+
+pop_size = size(population,2);
+
+indices = 1:pop_size;
+trans = [indices(indices<=chosen) indices(indices>chosen)-1];
+indices = indices(indices~=chosen);
+population = population(indices); % suppression of the individual
+
+for i=1:pop_size-1
+ llinks = population(i).links ;
+ % suppress indice chosen and shift higher indice from -1
+ population(i).links = trans(llinks(llinks~=chosen));
+end
+
+
+function add_individual(chosen)
+%--------------------------------------------------------------------------
+% add a daughter node the the node "chosen" of population
+%--------------------------------------------------------------------------
+
+global population resources x_grid y_grid ...
+ muVM kaVM muLN s2LN LL ...
+ cons_std_dev cons_max cons_th
+
+pop_size = size(population,2);
+
+grad = approx_gradient_angle(chosen); % in this direction
+angle = grad+circ_vmrnd(muVM,kaVM,1);
+length = random('logn',muLN,s2LN);
+
+new_pos = population(chosen).pos+length*[cos(angle) sin(angle)];
+population(pop_size+1).pos = min(max(new_pos,0),LL); % forced in the field
+population(pop_size+1).resou = f_resources(population(pop_size+1).pos);
+population(pop_size+1).links = chosen;
+population(chosen).links = [population(chosen).links pop_size+1];
+
+% --- resources update
+
+coef = 1/(2*cons_std_dev*cons_std_dev);
+cons = cons_max*exp(-coef*((x_grid-new_pos(1)).^2 ...
+ + (y_grid-new_pos(2)).^2));
+cons = (cons>cons_th).*cons;
+resources = max(resources-cons,0);
+
+
+function gradient=approx_gradient(indice)
+%--------------------------------------------------------------------------
+% compute an approximation of the gradient of the resources in the
+% individual "indice" (if no neighbor: random shoot)
+%--------------------------------------------------------------------------
+global population
+neighbor = population(indice).links;
+nb_neighbor = length(neighbor);
+if (nb_neighbor>0)
+ vec = reshape([population(neighbor).pos],2,[])...
+ -repmat(population(indice).pos',1,nb_neighbor);
+ norm_vec2 = sum(vec.^2);
+ vec_resou = [population(neighbor).resou]...
+ -repmat(population(indice).resou,1,nb_neighbor);
+ vec_resou = vec_resou./norm_vec2;
+ gradient = sum(repmat(vec_resou,2,1).*vec,2)/nb_neighbor;
+else
+ % no point: random shoot
+ angle = 2*pi*rand();
+ gradient =[cos(angle);sin(angle)];
+end
+
+
+function angle=approx_gradient_angle(indice)
+%--------------------------------------------------------------------------
+% compute an approximation of the angle of the gradient of the resources
+% in the individual "indice" (if no neighbor: random shoot)
+%--------------------------------------------------------------------------
+global population
+neighbor = population(indice).links;
+nb_neighbor = length(neighbor);
+if (nb_neighbor>0)
+ vec = reshape([population(neighbor).pos],2,[])...
+ -repmat(population(indice).pos',1,nb_neighbor);
+ norm_vec2 = sum(vec.^2);
+ vec_resou = [population(neighbor).resou]...
+ -repmat(population(indice).resou,1,nb_neighbor);
+ vec_resou = vec_resou./norm_vec2;
+ gradient = sum(repmat(vec_resou,2,1).*vec,2)/nb_neighbor;
+ angle = atan2(gradient(2),gradient(1));
+else
+ % no point: random shoot
+ angle =2*pi*rand();
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% --- technical -----------------------------------------------------------
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function value = f_resources(pos)
+%--------------------------------------------------------------------------
+% given resources/x_grid/y_grid give the interpolated value in pos
+%--------------------------------------------------------------------------
+global resources x_grid y_grid
+value = interp2(x_grid,y_grid,resources,pos(1),pos(2));
+
+
+function x = f_finiternd(finite_distribution,n_sample)
+%--------------------------------------------------------------------------
+% sampling from a finite distribution
+% finite_distribution : finite distribution (not need to be normalised)
+% n_sample : sampling size
+%--------------------------------------------------------------------------
+if nargin == 1
+ n_sample = 1;
+end
+n_outcomes = length(finite_distribution); % number of outcomes
+finite_distribution = reshape(finite_distribution,n_outcomes,1);
+x = sum(repmat(rand(1,n_sample),n_outcomes,1)> ...
+ repmat(cumsum(finite_distribution)/sum(finite_distribution), ...
+ 1,n_sample),1)+1;
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% --- plotting routines ---------------------------------------------------
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+function create_video()
+%--------------------------------------------------------------------------
+% create a avi from the data.mat
+%--------------------------------------------------------------------------
+global result LL LLzoom work_dir resources_max
+
+num_frame = size(result,2);
+
+fprintf(' %g frames available \n',num_frame)
+pdf_frames=str2double(input(' how many frames to pdf ? ','s'));
+frames_indices=floor(1:(num_frame-1)/pdf_frames:num_frame);
+
+
+fprintf(' makeing the avi... \n')
+
+ffig = figure('visible','off'); %turns visibility of figure off
+aviobj = avifile([work_dir '/animation.avi']);
+
+
+for i_frame=frames_indices
+ disp([' frame ' num2str(i_frame) ' / ' num2str(num_frame)])
+ imagesc([0 LL], [0 LL],result(i_frame).res,[0 resources_max])
+ colormap(flipud(gray));
+ axis xy;
+ axis(LLzoom);
+ axis square;
+ hold on
+ plot_population_nodes(i_frame)
+ drawnow
+ aviobj=addframe(aviobj,getframe(gca)); %adds frames to the AVI file
+end
+
+
+aviobj=close(aviobj); %closes the AVI file
+close(ffig); %closes the handle to invisible figure
+fprintf('done \n')
+
+
+
+function plot_population_nodes(i_frame)
+%--------------------------------------------------------------------------
+% ploting the nodes in the population with results
+%--------------------------------------------------------------------------
+global result node_size node_color
+
+positions = reshape([result(i_frame).pos],2,[]);
+plot(positions(1,:),positions(2,:),...
+ 'go','LineStyle','none',...
+ 'MarkerFaceColor',node_color,'MarkerEdgeColor',...
+ node_color,'MarkerSize',node_size)
+
+
+
+function plot_landscape()
+%--------------------------------------------------------------------------
+% plot the landscape: the resource map and the plant
+%--------------------------------------------------------------------------
+
+global population resources resources_max LL LLzoom ...
+ node_size node_color link_color type_graph
+
+clf
+imagesc([0 LL], [0 LL],resources,[0 resources_max])
+colormap(flipud(gray));
+hold on
+if ~(isempty(population))
+ pop_size = size(population,2);
+ if type_graph~=1
+ % --- plot the links
+ for i=1:pop_size
+ for j = population(i).links
+ if (i<=j)
+ line([population(i).pos(1) population(j).pos(1)],...
+ [population(i).pos(2) population(j).pos(2)],...
+ 'Color',link_color)
+ end
+ end
+ end
+ end
+ if type_graph~=2
+ % --- plot the nodes
+ positions = reshape([population.pos],2,[]);
+ plot(positions(1,:),positions(2,:),...
+ 'gs','LineStyle','none',...
+ 'MarkerFaceColor',node_color,'MarkerEdgeColor',...
+ node_color,'MarkerSize',node_size)
+ end
+end
+
+axis xy;
+axis(LLzoom);
+axis square;
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% --- fixe graphics
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function full_graphics()
+
+global work_dir node_size node_color link_color type_graph
+
+hfig=figure('Visible','off');
+clf
+node_size = 2;
+node_color = 'r';
+link_color = 'g';
+type_graph = 3;
+plot_landscape()
+saveas(gcf,[work_dir '/graph_last_landscape.pdf'])
+saveas(gcf,[work_dir '/graph_last_landscape.fig'])
+clf
+plot_links_stat()
+saveas(gcf,[work_dir '/graph_last_links_stat.pdf'])
+saveas(gcf,[work_dir '/graph_last_links_stat.fig'])
+clf
+plot_shoot_profile()
+saveas(gcf,[work_dir '/graph_shoot_profile.pdf'])
+saveas(gcf,[work_dir '/graph_shoot_profile.fig'])
+clf
+plot_pop_size_res()
+saveas(gcf,[work_dir '/graph_pop_size_res.pdf'])
+saveas(gcf,[work_dir '/graph_pop_size_res.fig'])
+clf
+plot_pop_fitness()
+saveas(gcf,[work_dir '/graph_pop_fitness.pdf'])
+saveas(gcf,[work_dir '/graph_pop_fitness.fig'])
+clf
+plot_link_length_pdf()
+saveas(gcf,[work_dir '/graph_link_length_pdf.pdf'])
+saveas(gcf,[work_dir '/graph_link_length_pdf.fig'])
+clf
+plot_link_angle_pdf()
+saveas(gcf,[work_dir '/graph_link_angle_pdf.pdf'])
+saveas(gcf,[work_dir '/graph_link_angle_pdf.fig'])
+close(hfig)
+
+function some_snapshots()
+%--------------------------------------------------------------------------
+% some landcsape at a given number of time
+%--------------------------------------------------------------------------
+global result LL LLzoom work_dir resources_max
+
+
+num_frame = size(result,2);
+fprintf(' %g frames available \n',num_frame)
+pdf_frames=str2double(input(' how many frames to pdf ? ','s'));
+frames_indices=floor(1:(num_frame-1)/pdf_frames:num_frame);
+
+
+fprintf(' makeing the avi... \n')
+
+ffig = figure('visible','off'); %turns visibility of figure off
+
+for i_frame=frames_indices
+ disp([' frame ' num2str(i_frame) ' / ' num2str(num_frame)])
+ imagesc([0 LL], [0 LL],result(i_frame).res,[0 resources_max])
+ colormap(flipud(gray));
+ axis xy;
+ axis(LLzoom);
+ axis square;
+ hold on
+ plot_population_nodes(i_frame)
+ saveas(gcf,[work_dir '/graph_landscape' num2str(i_frame,'%04i\n') '.pdf'])
+end
+
+
+close(ffig); %closes the handle to invisible figure
+fprintf('done \n')
+
+
+function pre_graphics()
+
+clf
+subplot(3,1,1)
+plot_shoot_profile()
+subplot(3,1,2)
+plot_link_length_pdf()
+subplot(3,1,3)
+plot_link_angle_pdf()
+
+
+function plot_links_stat()
+%--------------------------------------------------------------------------
+% hist of # links in the last landsapce
+%--------------------------------------------------------------------------
+global population
+pp = size(population,2);
+if pp>10
+ nb_links = zeros(pp,1);
+ for i=1:pp
+ nb_links(i) = size(population(i).links,2);
+ end
+ lmin = min(nb_links);
+ lmax = max(nb_links);
+ hist(nb_links,lmin:lmax);
+ h = findobj(gca,'Type','patch');
+ set(h,'FaceColor','b','EdgeColor','w')
+ box on
+ xlabel('number of links')
+ ylabel('frequency')
+end
+
+function plot_shoot_profile()
+%--------------------------------------------------------------------------
+% plot some shoots
+%--------------------------------------------------------------------------
+global muVM kaVM muLN s2LN
+for i=1:20
+ angle = circ_vmrnd(muVM,kaVM,1);
+ length = random('logn',muLN,s2LN);
+ line([0 length*cos(angle)],...
+ [0 length*sin(angle)],'Color','r')
+end
+box on
+axis square
+daspect([1 1 1])
+title('shoot profile')
+
+
+
+function plot_pop_size_res()
+%--------------------------------------------------------------------------
+% evolution of the size and fitness of pop
+%--------------------------------------------------------------------------
+global instants pop_sizes total_res total_event
+
+% [AX,H1,H2] = plotyy(instants(1:total_event),pop_sizes(1:total_event),...
+% instants(1:total_event),total_res(1:total_event));
+% set(get(AX(1),'Ylabel'),'String','population size')
+% set(get(AX(2),'Ylabel'),'String','resource evolution')
+% set(AX(2),'XTick',instants)
+% xlabel('time')
+% xlim(AX(1),[0 instants(total_event)])
+% xlim(AX(2),[0 instants(total_event)])
+%axis([AX(1) AX(2)],'tight')
+
+hl1 = line(instants(1:total_event),pop_sizes(1:total_event),'Color','r');
+ax1 = gca;
+set(ax1,'XColor','k','YColor','r')
+ax2 = axes('Position',get(ax1,'Position'),...
+ 'XAxisLocation','top',...
+ 'YAxisLocation','right',...
+ 'Color','none',...
+ 'XColor','k','YColor','b','TickDir','in','XTickLabel',{});
+% ax2 = axes('Position',get(ax1,'Position'),...
+% 'XAxisLocation','top',...
+% 'YAxisLocation','right',...
+% 'Color','none',...
+% 'XColor','k','YColor','b','TickDir','in','XTickLabel',{},...
+% 'XTick',instants(1:total_event));
+hl2 = line(instants(1:total_event),total_res(1:total_event),'Color','b','Parent',ax2);
+set(get(ax1,'Ylabel'),'String','population size')
+set(get(ax2,'Ylabel'),'String','resource evolution')
+xlim(ax1,[0 instants(total_event)])
+xlim(ax2,[0 instants(total_event)])
+
+
+
+
+function plot_pop_fitness()
+%--------------------------------------------------------------------------
+% evolution of the ~fitness of pop
+%--------------------------------------------------------------------------
+global instants pop_fitness
+plot(instants,pop_fitness);
+ylabel('fitness')
+xlabel('time')
+axis tight
+
+
+function plot_link_length_pdf()
+%--------------------------------------------------------------------------
+% plot link length pdf
+%--------------------------------------------------------------------------
+global muLN s2LN
+xx = 0:(0.5/100):0.5 ;
+plot(xx,lognpdf(xx,muLN,s2LN))
+title('link lenght pdf')
+axis tight
+
+
+function plot_link_angle_pdf()
+%--------------------------------------------------------------------------
+% plot_link_angle_pdf
+%--------------------------------------------------------------------------
+global muVM kaVM
+xx = -pi:(2*pi/100):pi ;
+plot(xx,circ_vmpdf(xx,muVM,kaVM))
+title('link angle pdf')
+axis tight
+
+
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Attribution-NonCommercial 3.0 Unported (CC BY-NC 3.0)
+% Paternit� - Pas d'Utilisation Commerciale CC BY-NC
+% You are free
+% to Share -- to copy, distribute and transmit the work
+% to Remix -- to adapt the work
+% Under the following conditions:
+% Attribution -- You must attribute the work in the manner specified by
+% the author or licensor (but not in any way that suggests
+% that they endorse you or your use of the work).
+% Noncommercial -- You may not use this work for commercial purposes.
+%
+% see details in http://creativecommons.org/licenses/by-nc/3.0/
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+