function [m_res_strategy,ESSfound,I_res,N_res,m_res_realized,p_by_stage,p_by_distance,IavgNM,IavgM] = do_search_ESS(lambda,beta_min,gamma,alpha_s,mu,psi,phi_max,alpha_f,ns,pflag)
% written by Allison Shaw, Contact ashaw@umn.edu for assistance.
% started: 28 September 2021
% last updated: 18 October 2021
%
% This searches for the ESS migration strategy for a population where
%   infection and migration each incurs separate costs as well as 
%   potentially interactive costs.
%
% INPUTS:
%   lambda = time to half change in infection
%   beta_min = scalar of infection dynamics at x=1
%   gamma = midpoint of survival-v-intensity
%   alpha_s = slope of survival-v-intensity
%   mu = slope of survival-v-migration
%   psi = degree of migratory culling (interaction infection/migration)
%   phi_max = maximum per capita fecundity
%   alpha_f = slope of fecundity-v-intensity
%   ns = number of infection intensity stages
%   pflag = whether (1) or not (0) to plot time series when ESS not found
%
% OUTPUTS:
%   m_res_strategy = final resident migration strategy
%   ESSfound = whether (1) or not (0) the ESS was found
%   I_res =  final intensity of each intensity stage (column vector)
%   N_res = final number of individuals in each intensity stage (column vector)
%   m_res_realized = for residents, vector of the migration distances that each intensity stage
%      present in the population goes
%   p_by_stage = proportion of population in each infection intensity stage (column vector)
%   p_by_distance = proportion of population traveling each migration distance (column vector)
%   IavgNM = average infection intensity of Non-Migrants
%   IavgM = average infection intensity of Migrants

%%---MODEL SET UP---------------------------------------------------------%
    del_mig = 0.2; % spacing for different migration strategies
    eps2 = 1e-10; % error for invasion rate (must be within this to be same as fastest growing mutant)
    %
    nmut = 100;   % number of mutants to try each round of invasion
    count_tot_max = 100; % maximum number of times to attempt invasion
    count_res_max = 10;  % ESS threshold: strategy is an ESS if it successfully defends this many times
    %
    I_IC = 0; % starting intensity of each intensity stage (column vector)
    N_IC = 1; % starting number of individuals in each intensity stage (column vector)
    %
    m_res_strategy = zeros(ns,1); % initial resident migration strategy (column vector) for
    m_res_history = []; % to save the resident strategies for each invasion attempt
    p_stage_history = []; % to save the p-by-stage binning for each invasion attempt
    p_distance_history = []; % to save the p-by-stage binning for each invasion attempt
    %
    count_tot = 0; % number of total invasion attempts
    count_res = 0; % number of times resident strategy has resisted invasion
%%---MODEL SET UP---------------------------------------------------------%
    
% run invasion attempts until one of following happens:
%  * ESS is reached (count_res > count_res_max)
%  * run out of invasion attempts (count_tot > count_tot_max)
while count_tot < count_tot_max && count_res < count_res_max

    count_tot = count_tot+1; % count this invasion attempt

    % run ecological simulations for resident type
    eflag = 'eco'; % whether to run 'eco' or 'evo' simulations
    [I_res,N_res,m_res_realized,~,IavgNM,IavgM] = do_run_simulation(eflag,lambda,beta_min,gamma,alpha_s,mu,psi,phi_max,alpha_f,I_IC,N_IC,m_res_strategy,[]);
    
    % bin population by infection intensity stage and by migration distance
    % (column vectors)
    [p_by_stage,p_by_distance] = do_bin_population(I_res,N_res,m_res_realized,ns,del_mig);
    % save these in history
    p_stage_history = [p_stage_history p_by_stage];
    p_distance_history = [p_distance_history p_by_distance];

    if ~isempty(N_res)  % only invade if resident population is viable
        % create set of mutants with strategies all differ a bit from resident
        m_mut_mat = NaN(ns,nmut);
        m_mut_mat(:,1) = m_res_strategy; % set first mutant to be the same as resident
        del_mut = del_mig*randi([-1, 1],ns, nmut-1); % set mutation amount for others
        m_mut_mat(:,2:end) = min(1,max(0,m_res_strategy+del_mut)); % bound mutants between 0 and 1

        gro_mut = NaN(1,nmut); % empty matrix to hold growth rate for each mutant strategy

        % loop over each mutant
        for i = 1:nmut
            % get mutant strategy
            m_mut_strategy = m_mut_mat(:,i);
            % run simulations
            eflag = 'evo'; % whether to run 'eco' or 'evo' simulations
            [~,~,~,eigJ] = do_run_simulation(eflag,lambda,beta_min,gamma,alpha_s,mu,psi,phi_max,alpha_f,I_res,N_res,m_res_strategy,m_mut_strategy);
            % save growth rate of mutant
            gro_mut(i) = eigJ;
        end

        % find the index of the mutant that grew the most, within the error
        max_mut_ind = find(abs(gro_mut - max(gro_mut)) < eps2);
        % if more than one, take the first (which gives preference to the resident)
        max_mut_ind = max_mut_ind(1);

        % save best mutant as new resident
        m_res_strategy = m_mut_mat(:,max_mut_ind); 
        % save resident in history
        m_res_history = [m_res_history m_res_strategy];

        % tally whether resident resisted invasion or not
        if max_mut_ind == 1
            count_res = count_res + 1; % count as successful resistance to invasion
        else
            count_res = 0;             % reset resisting invasion counter
        end
    end

end

emptystages = find(p_by_stage<eps2); % find effectively empty
m_res_strategy(emptystages) = NaN;    % mask out these stages

% record whether simulations stopped because an ESS was found or we ran out
% of invasion attempts
% check for viable final population size
if sum(N_res)<1
    ESSfound = 0; % record no ESS
    error('population collapsed');
    m_res_strategy = NaN(ns,1);
elseif count_res == count_res_max
    ESSfound = 1; % record yes ESS
else
    ESSfound = 0; % record no ESS
    disp('ESS not found')
    
    % modify ESS migration strategy to be average
    m_res_strategy = mean(m_res_history(:,end-9:end),2);
    m_res_strategy(emptystages) = NaN;    % mask out these stages
    m_res_history(emptystages,:) = NaN;
    
    % modify binned p's to be average
    p_by_stage = mean(p_stage_history(:,end-9:end),2);
    p_by_distance = mean(p_distance_history(:,end-9:end),2);
    
    if pflag==1 % if plotting on
        % plot metric history if sims haven't converged to ESS
        figure(10); clf
        cmap = parula;
        cmap = [1 1 1; cmap; 0 0 0]; % add white for NaNs, black for inf
        colormap(cmap)
        
        subplot(2,2,1)
        imagesc(flipud(m_res_history),[-0.01 1.01])
        xlabel('time')
        ylabel('Migration distance')
        title('ESS migration strategy')
        set(gca,'YTick',[])
        set(gca,'XTick',[])
        colorbar
        
        subplot(2,2,2)
        imagesc(flipud(p_stage_history))
        xlabel('time')
        ylabel('Infection intensity stage')
        title('proportion')
        set(gca,'YTick',[])
        set(gca,'XTick',[])
        colorbar
        
        subplot(2,2,3)
        imagesc(flipud(p_distance_history))
        xlabel('time')
        ylabel('Migration distance')
        title('proportion')
        set(gca,'YTick',[])
        set(gca,'XTick',[])
        colorbar
        
        saveas(10,strcat(['fig_time_lambda=' num2str(lambda) '_betamin=' num2str(beta_min) '_gamma=' num2str(gamma) '_alphas=' num2str(alpha_s) '_mu=' num2str(mu)  '_psi=' num2str(psi) '_alphaf=' num2str(alpha_f)  '_phimax=' num2str(phi_max) '.jpg']))
    end
end
