function [x,n,speed,popedge,core_avg,edge_avg] = run_simulation(lambda,b,a,v,mu,tau,ngens,sflag,kflag,beta,ICflag)
% USAGE: [x,n,speed,popedge,core_avg,edge_avg] = run_simulation(lambda,b,a,v,mu,tau,ngens,sflag,kflag,beta,ICflag)
%
% model written by Allison Shaw (contact for assistance: ashaw@umn.edu)
%   last updated: 3 January 2023
%
% Simulates a population with heterogeneity in dispersal propensity. All
%   dispersers have a fixed mortality cost. Simulations can be run with or
%   without spatial sorting.
%
% INPUTS:
%   lambda = density-independent growth rate
%   b = density-dependence parameter
%   a = Allee threshold
%   v = dispersal variance
%   mu = vector of length tau, mortality of each disperser type
%   tau = number of types to track in the population
%   ngens = number of generations
%   sflag = whether to run without (0) or with (1) spatial sorting
%   kflag = what dispersal kernel to use: laplace (0), gaussian (1), or
%      gaussian with distance-dependent dispersal mortality (2)
%   beta = dispersal mortality constant, to use only if kflag is 2
%   ICflag = what initial conditions to use: uniform (0), skewed to first type (1), skewed to last type (2)
%
% OUTPUTS:
%   x = location of population points
%   n = density of each type, over time and space
%   speed = instantaneous spread rate of host population
%   popedge = x location of the right-most edge of the total population for
%        each gen
%   core_avg = average dispersal phenotype at core
%   edge_avg = average dispersal phenotype at edge

%-----MODEL SET UP--------------------------------------------------------%
    eps1 = 1e-15;    % threshold for a 'zero' population density
    eps2 = 1e-1;     % threshold for good kernel integration
    ncrit = 0.001;   % threshold for the edge of the population

    nodes = (2^15)+1;      % number of nodes/bins to have in the domain (2^m + 1)
    diameter = 700;        % length of domain/space to simulate across
    radius = diameter/2;

    % set up 1-D spatial domain
    x = linspace(-radius,radius,nodes);          % vector of node locations
    x2 = linspace(-diameter,diameter,2*nodes-1); % vector of node locations, extended domain
    dx = diameter/(nodes-1);                     % distance between nodes

    pvec = linspace(0.1,1,tau); % propensity to disperse
    
    if kflag==0
        % dispersal kernel (laplace)
        k = exp(-sqrt((2*(x2).^2)/v))./sqrt(2*v);
    elseif kflag==1
        % dispersal kernel (gaussian)
        k = exp(-0.5*(x2.^2/v))./sqrt(2*pi*v);
    elseif kflag==2
        % dispersal kernel (gaussian)
        k = exp(-0.5*(x2.^2/v))./sqrt(2*pi*v);
    else
        error('unrecognized kflag')
    end
    if abs(1-trapz(x2,k)) > eps2; error('bad dispersal kernel integration: increase nodes'); end
    
    % add in dispersal mortality of using gaussian with mortality
    if kflag==2; k = exp(-beta)*k; end
    
    % matrix to store population densities across time (dim 1), type (dim 2),
    % and space (dim 3)
    n = zeros(ngens+1,tau,length(x));
    
    % variables to store information
    popedge = NaN(ngens+1,1);
    speed = NaN(ngens,1);
    core_avg = NaN;
    edge_avg = NaN;
%-----MODEL SET UP--------------------------------------------------------%

%set up initial conditions
IC_rad = 0.5; % radius of initial population
temp = find(abs(x) <= IC_rad);
if ICflag==0 % uniform
    IC_dens = 1;  % initial density
    n(1,:,temp) = IC_dens*ones(size(n(1,:,temp))); % set pop with initial conditions
elseif ICflag==1 % first
    IC_dens = [5 (5/(tau-1))*ones(1,tau-1)];  % initial density - skewed to first value of p
    n(1,:,temp) = repmat(IC_dens,1,1,length(temp)).*ones(size(n(1,:,temp))); % set pop with initial conditions
elseif ICflag==2 % last
    IC_dens = [(5/(tau-1))*ones(1,tau-1) 5];  % initial density - skewed to last value of p
    n(1,:,temp) = repmat(IC_dens,1,1,length(temp)).*ones(size(n(1,:,temp))); % set pop with initial conditions
else
    error('unrecognized ICflag')
end
clear temp

Nvia = sum(squeeze(sum(n(1,:,:)))>ncrit); % number viable pop locations
i = 0;

while i<ngens && Nvia > 0
    i = i+1;
    
    % INITIAL CONDITIONS
    n0 = squeeze(n(i,:,:));  % tau x nodes

    % GROWTH
    ntot = sum(n0,1);
    g = b./(b+ntot); % strength of density dependence
    g(ntot<a) = 0;     % account for Allee effect
    
    if sflag==0 % no spatial sorting
        n1 = repmat(lambda*ntot.*g./tau,tau,1); % offspring are evenly distributed
    elseif sflag==1 % with spatial sorting
        n1 = lambda*n0.*repmat(g,tau,1); % offspring inherit parent's type
    else
        error('unrecognized sflag')
    end
    
    % DISPERSAL
    % calculate convolutions
    % only pvec(j) of each type disperse
    % only 1-mu of dispersers survive
    for j = 1:tau
        nx(j,:) = fft_conv(k,(1-mu(j))*pvec(j).*n1(j,:));
    end
    % convert to population density
    n2 = dx*nx(:,nodes:length(x2));
    % fix ends
    n2(:,1) = n2(:,1)/2;
    n2(:,nodes) = n2(:,nodes)/2;
    % add in non-dispersing individuals
    n2 = n2 + (1-repmat(pvec',1,nodes)).*n1;   
    % set any densities that are too low to zero
    temp = find(n2 < eps1);
    n2(temp) = zeros(size(n2(temp)));
    clear j temp
    
    n(i+1,:,:) = n2;
    
    % reset temporary density vectors
    n0(:) = 0;
    n1(:) = 0;
    n2(:) = 0;
    nx(:) = 0;
        
    % overall population size
    N = squeeze(sum(n(i+1,:,:),2))';
    % find the edge of the total population
    jj = find(N >= ncrit,1,'last');
    % if there is a population edge
    if ~isempty(jj)
        % interpolate to get the location
        popedge(i+1) = interp1(N(jj:jj+1),x(jj:jj+1),ncrit);
    end
    speed(i) = popedge(i+1)-popedge(i);
    
    Nvia = sum(squeeze(sum(n(i,:,:)))>ncrit); % number viable pop locations

end

% get average dispersal phenotype at core and edge
% only save if there are some individuals left
if Nvia > 0
    core_avg = sum(n(i+1,:,ceil(nodes/2)).*pvec)./sum(n(i+1,:,ceil(nodes/2)));
    edge_avg = sum(n(i+1,:,jj).*pvec)./sum(n(i+1,:,jj));
end