%GELAM.m
%**************************************************************%
% Golden Eagle Lead Abatement Model (GELAM) version 1          %
%                                                              %
%                                                              %
% Runs simulation to estimate acute mortality of golden        %
%  eagles due to ingesting lead from gut-piles of big game     %
%  animals and to assess possible mitigation strategies        %
%   to quantify mortality rates and uncertainty under          %
%   two scenarios:  non-lead ammunition and gut-pile retrieval %
%                                                              %
%  Written by:    Eric Lonsdorf and Carol A. Sanders-Reed      %
%  Last Updated:   April 10 2014                               %
%      (Updated only comments May 21 2014)                     %
%                                                              %
%  Input files:   location2West.xls                            %
%  Output files:                                               %
%      dynamic name: GELAM_date_time_number of iterations      %                             
%      example:  GELAM_10-Apr-2014_2155_5000                   %
%**************************************************************%

%**************************************************************%
%  Set up simulation run -- clear variables and plots          %
%                           start clock                        %
%  Seed random number generator so results can be replicated   %
% Note:  no significance to seed number of 123                 %
%**************************************************************%
startclock=clock;
tic;

close all 
clear all
rand('seed',123);

script_name = 'GELAM';

%**************************************************************%
% Set parameters for all iterations                            %
%**************************************************************%
% Number of iterations run 
num_iters = 5000;

% Range of mitigation effects  
% gutpile mitigation alpha1_levels = % gutpiles removed
alpha1_levels = 0:0.1:1;

% ammunition mitigation alpha2_levels % nonlead ammunition
alpha2_levels = 0:0.1:1;

% Maximum lead (pb) concentration levels (pb_levels)
pb_levels = 0:10:1000;
num_pb_levels = 101;

% Number of days to blood half-life (T-half in Equation 5)
max_T_half_blood_pb = 20;   
min_T_half_blood_pb = 10;

%% Cauchy (parameters in Equation 3)
min_E_mode   = 25;
max_E_mode   = 75;
gamma_Cauchy = 25;

% Scavenge parameters - note low and high refer to the mortality line
%  (parameters in Equation 1)
half_sat_gp_low = 5;
half_sat_gp_high = 10;
power_gp_low = 3;
power_gp_high = 1.5;
scalar_gp = 1;
max_gp_low = 1;
max_gp_high = 5;
 
% Mortality parameters (parameters in Equation 7)
half_sat_mort_low  = 150;
half_sat_mort_high = 700;
power_mort         = 2.5;

% Carcass parameters
 % max number of gut piles eaten
max_carcass = 5;
carcass_levels = 0:5;
num_carcass_levels = length(carcass_levels);

%**************************************************************%
% Load current location information for hunting units
%**************************************************************%
locationsraw= xlsread('location2West.xlsx');
% c1 = FID; c2 = number of 100km blocks; c3 = total shot 
% c4 = density eagles per km c5 = hunting area 

% Select only the hunting units with eagle density estimates in West data
locations = locationsraw(locationsraw(:,4)>0,:);

% 90% of animals shot result in gutpiles
locations(:,3) = locations(:,3) * 0.9;
num_loc = size(locations,1);

% clear parameters 
clear locationsraw

%*************************************************************%
% Pre-allocate variables to hold simulation results            %
%**************************************************************%

% Iteration values: 
prob_pb           = zeros(1,num_pb_levels);
prob_pb_temp      = zeros(1,num_pb_levels);
max_gp            = zeros(num_iters,1);
T_half_blood_pb   = zeros(num_iters,1);
L                 = zeros(num_iters,1);
E_mode            = zeros(num_iters,1);
half_sat_mort     = zeros(num_iters,1);
gp_no_pb          = zeros(num_iters,1);
prob_pb_no_pb     = zeros(num_iters,num_pb_levels,1);
m_gp_no_mit       = zeros(num_iters,num_loc,1);
prob_pb_no_pb_mit = zeros(num_iters,10,num_pb_levels,1);

% Mitigation values: 
m_gp              = zeros(num_iters,num_loc,length(alpha1_levels),1);
mort_ratea2       = zeros(num_iters,num_loc,length(alpha2_levels),1);  
mort_rate         = zeros(num_iters,num_loc,10,1);

%**************************************************************%
% Begin simulation runs                                        %
%**************************************************************%

for rr = 1:num_iters
   
    %**************************************************************%
    % Set parameters that vary by iteration                        %
    %**************************************************************%
    scavenge_rand = rand;
    
    % maximum number of gut piles ingested
    max_gp(rr,1) =1+4*scavenge_rand;  
    
    % determine days between eating gut piles  (Equation 4)
    % If ave gut piles eaten is equal to 1 set L to between 1 and 29
    if max_gp(rr,1) < 2
        L(rr,1) = round(3 + rand * 26);
    else
        L(rr,1) = round(3 + rand*(30/(max_gp(rr,1)-1)-3));
    end
                           
    % number of days before blood lead decays by 50%  T half in Eq 5
    T_half_blood_pb(rr,1) = round(min_T_half_blood_pb + ...
                rand*(max_T_half_blood_pb - min_T_half_blood_pb)); 
    
    % blood lead decay calculation (Equation 4)
    D = exp(log(0.5)/ T_half_blood_pb(rr,1));
                       
    % Lead per gut pile (D^L parameter in Equation 6)
    lpgp  = D^L(rr,1);
    
    clear decay_rate
    
    % Percent gut piles with no lead  Equation 3b
    gp_no_pb(rr,1)  = 0.1 + rand*(0.40);
    
    % determining lead from gut piles consumed     
    lead_consump = round(([0:max_carcass]')*(pb_levels)); 
     
    % mode of blood lead ingested per gp -- Cauchy distribution
    % E-Mode in Equation 3
    E_mode(rr,1) = min_E_mode + rand*(max_E_mode - min_E_mode);
 
    % Probability distribution for lead ingested per gp 
    %  Equation 6
    for i = 1:num_pb_levels
        for j = 0:max_carcass
            if j==0
               lead_consump(j+1,i) = 0;
            elseif j ==1
                lead_consump(j+1,i) = pb_levels(i);
            else
                lead_consump(j+1,i) = lead_consump(j,i)* lpgp + ...
                   pb_levels(i);
            end
        end
        % Equation 3 -- Cauchy pdf
        prob_pb_temp(i) = (gamma_Cauchy/pi)*(1/((pb_levels(i)-E_mode(rr,1))^2 + ...
                        gamma_Cauchy^2));            
    end
        
    % need to standardize so that probability distribution equals one
    % Equation 3b
    prob_pb = prob_pb_temp / ...
                sum(prob_pb_temp);
    
    
    %% probability distribution of lead per gp removing % gutpiles with
    %% no lead   Equation 3b
    prob_pb_no_pb(rr,1) = prob_pb(1)*(1 - gp_no_pb(rr,1))...
                                        + gp_no_pb(rr,1);
    prob_pb_no_pb(rr,2:num_pb_levels) = ...
                    prob_pb(1,2:num_pb_levels)*(1 - gp_no_pb(rr,1));
    
    % Clear variables                
    clear prob_pb prob_pb_temp
    
    %**************************************************************%
    % Mortality calculation                                        %
    %**************************************************************%
    mort_rand = rand;
     
    half_sat_mort(rr,1) = half_sat_mort_low + mort_rand * ...
                            (half_sat_mort_high - half_sat_mort_low);
    
    % Equation 7
    mort_prob = (lead_consump.^power_mort)./...
                   (lead_consump.^power_mort + ...
                    half_sat_mort(rr,1)^power_mort);
     
    prop_carc = zeros(1,num_carcass_levels);
            
    %**************************************************************%
    %  Location Loop
    %*************************************************************
    for gg = 1:num_loc
        % eagle density / km * 100 = eagles/100 km block
        ee = locations(gg,4)*100;
    
        %% Scenarios of mitigation
    
        %**************************************************************%
        %  mitigation loop
        %**************************************************************%
        for aa = 1:1:length(alpha1_levels)  %0:0.1:1 % carcass removal
        % gutpiles per 100 km2 blocks per eagle
        %  number of gutpiles / (number of 100 km blocks * # eagles/km *
        %  100)
        gp = (locations(gg,3)/(locations(gg,2)*ee))*(1-alpha1_levels(aa));
                
        % Equation 1 (when aa = 0) Equation 11 (when aa > 0)
        m_gp_low = max_gp_low * scalar_gp * gp ^ power_gp_low / ...
                 (scalar_gp*gp ^ power_gp_low +  ...
                    half_sat_gp_low ^ power_gp_low);
        m_gp_high = max_gp_high * scalar_gp * gp ^ power_gp_high / ...
                     (scalar_gp*gp ^ power_gp_high + ...
                     half_sat_gp_high ^ power_gp_high);

        %% determine ave number of gut piles
        if gp ==  0     
            m_gp(rr,gg,aa,1) = 0;
            prop_carc(:) = 0;
        else
            % choosing random value based on Equation 1
            m_gp(rr,gg,aa,1) =  ...
                   m_gp_low + scavenge_rand * (m_gp_high - m_gp_low); 
            % Equation 2
            if m_gp(rr,gg,aa,1) < max_carcass/2
                prop_carc(:) = ...
                    poisspdf(0:max_carcass,m_gp(rr,gg,aa,1));
            else
                prop_carc(:) =  ...
                    poisspdf(max_carcass - (0:max_carcass),max_carcass - ... 
                            m_gp(rr,gg,aa,1));
            end
            % Sum to 1
            temp_sum = sum(prop_carc(2:max_carcass+1)); 
            prop_carc_2(2:max_carcass+1) = ...
                        prop_carc(2:max_carcass+1)/temp_sum;
            temp_val = max([m_gp(rr,gg,aa,1), ...
                sum([1:max_carcass].*prop_carc_2(2:max_carcass+1))]);
            prop_carc_2(1) = 1-m_gp(rr,gg,aa,1)/temp_val;
            temp_sum2 = sum(prop_carc_2(:));
            prop_carc(:) = prop_carc_2(:)/temp_sum2;
        end
        
        %Equation 8
        prob_ingest = prop_carc'*squeeze(prob_pb_no_pb(rr,:));
                 
        %% estimate mortality  Equation 9
        mort_rate(rr,gg,aa) = sum(sum(prob_ingest.* mort_prob));
                  
        end    % for aa = 1:1:10
        
        % gutpiles per 100 km2 blocks per eagle without mitigation
        %  number of gutpiles / (number of 100 km blocks * # eagles/km *
        %  100)
        gp = locations(gg,3)/(locations(gg,2)*ee);
                
        % Equation 1
        m_gp_low = max_gp_low * scalar_gp * gp ^ power_gp_low / ...
                 (scalar_gp*gp ^ power_gp_low +  ...
                    half_sat_gp_low ^ power_gp_low);
        m_gp_high = max_gp_high * scalar_gp * gp ^ power_gp_high / ...
                     (scalar_gp*gp ^ power_gp_high + ...
                     half_sat_gp_high ^ power_gp_high);

        %% determine ave number of gut piles without mitigation
        if gp ==  0     
            m_gp_no_mit(rr,gg,1) = 0;
            prop_carc(:) = 0;
        else
            % Select random number based on Equation 1
            m_gp_no_mit(rr,gg,1) =  ...
                   m_gp_low + scavenge_rand * (m_gp_high - m_gp_low); 
            % Equation 2
            if m_gp_no_mit(rr,gg,1) < max_carcass/2
                prop_carc(:) = ...
                    poisspdf(0:max_carcass,m_gp_no_mit(rr,gg,1));
            else
                prop_carc(:) =  ...
                    poisspdf(max_carcass - (0:max_carcass),max_carcass - ... 
                            m_gp_no_mit(rr,gg,1));
            end
            %Sum to 1
            temp_sum = sum(prop_carc(2:max_carcass+1)); 
            prop_carc_2(2:max_carcass+1) = ...
                        prop_carc(2:max_carcass+1)/temp_sum;
            temp_val = max([m_gp_no_mit(rr,gg,1), ...
                sum([1:max_carcass].*prop_carc_2(2:max_carcass+1))]);
            prop_carc_2(1) = 1-m_gp_no_mit(rr,gg,1)/temp_val;
            temp_sum2 = sum(prop_carc_2(:));
            prop_carc(:) = prop_carc_2(:)/temp_sum2;
        end

        for aa2 = 1:1:length(alpha2_levels)
            
        %% probability distribution of lead per gp -- mitigation by a2
        %% amount   Equation 12
        prob_pb_no_pb_mit(rr,aa2,1) = prob_pb_no_pb(rr,1)*(1 - alpha2_levels(aa2))...
                                        + alpha2_levels(aa2);
        prob_pb_no_pb_mit(rr,aa2,2:num_pb_levels) = ...
                    prob_pb_no_pb(rr,2:num_pb_levels)*(1 - alpha2_levels(aa2));
        %Equation 8
        prob_ingest = prop_carc'*squeeze(prob_pb_no_pb_mit(rr,aa2,:))';
                 
        %% estimate mortality -- Equation 9
        mort_ratea2(rr,gg,aa2) = sum(sum(prob_ingest.* mort_prob));
        end  % aa2 end
            
    end  % for gg = 1:1:num_loc
end  % for rr = 1:1;num_iters

          
%**************************************************************%
%   Calculate total deaths for state
%*************************************************************
toteagle = locations(:,2).*locations(:,4)*100;
deadeaglea1 = zeros(num_iters,num_loc,length(alpha1_levels));
deadeaglea2 = zeros(num_iters,num_loc,length(alpha1_levels));
deadeaglesum = zeros(num_iters,length(alpha1_levels));
for rr=1:num_iters
    for aa = 1:length(alpha1_levels)
        deadeaglea1(rr,:,aa) = squeeze(mort_rate(rr,:,aa))'.* toteagle;
        deadeaglea2(rr,:,aa) = squeeze(mort_ratea2(rr,:,aa))'.* toteagle;
        deadeaglesuma1(rr,aa) = sum(deadeaglea1(rr,:,aa));
        deadeaglesuma2(rr,aa) = sum(deadeaglea2(rr,:,aa));
    end; 
    
end

% for entire state
state_mort_ratea1 = deadeaglesuma1/sum(toteagle);
state_mort_ratea2 = deadeaglesuma2/sum(toteagle);
state_prcta1 = prctile(state_mort_ratea1,[0 10 20 30 40 50 60 70 80 90 100]);
state_prcta2 = prctile(state_mort_ratea2,[0 10 20 30 40 50 60 70 80 90 100]);
state_meana1 = mean(state_mort_ratea1);
state_meana2 = mean(state_mort_ratea2);

% Calculate for Casper WY region
% these FID correspond to hunting units 22 34 66 67 88 89
selectFID = [57 70 129 88 84 82]; 
toteagleCasper = 0;
deadeaglea1Casper = zeros(num_iters,6,length(alpha1_levels));
deadeaglea2Casper = zeros(num_iters,6,length(alpha1_levels));
for i = 1:6 
    loc_index = find(locations(:,1) == selectFID(i));  
    eagle_loc =  locations(loc_index,2) * locations(loc_index,4) * 100;
    toteagleCasper = toteagleCasper + eagle_loc;                   
    for rr=1:num_iters
        for aa = 1:length(alpha1_levels)
            deadeaglea1Casper(rr,i,aa) =  ... 
                      squeeze(mort_rate(rr,loc_index,aa))'.* eagle_loc;
            deadeaglea2Casper(rr,i,aa) = ...
                      squeeze(mort_ratea2(rr,loc_index,aa))'.* eagle_loc;
        end; 
    end
end

for rr=1:num_iters
    for aa = 1:length(alpha1_levels)
        deadeaglesuma1Casper(rr,aa) = sum(deadeaglea1Casper(rr,:,aa));
        deadeaglesuma2Casper(rr,aa) = sum(deadeaglea2Casper(rr,:,aa));
    end
 end
 
 % for Casper
Casper_mort_ratea1 = deadeaglesuma1Casper/sum(toteagleCasper);
Casper_mort_ratea2 = deadeaglesuma2Casper/sum(toteagleCasper);
Casper_prcta1 = prctile(Casper_mort_ratea1,[0 10 20 30 40 50 60 70 80 90 100]);
Casper_prcta2 = prctile(Casper_mort_ratea2,[0 10 20 30 40 50 60 70 80 90 100]);
Casper_meana1 = mean(Casper_mort_ratea1);
Casper_meana2 = mean(Casper_mort_ratea2);

% save variables 
timesave = clock;
tsave = [num2str(timesave(4)), num2str(timesave(5)),'_', num2str(num_iters)];
savename = [script_name,'_',date,'_', tsave];
save(savename)

toc;