function out=hedgeSim4(reps, herit, cMu, cVar, b,d, weather, varBool, displayBool)
% hedgeSim4 implements the stochastic model of a population of 100
% flies reproducing and dying over a breeding season, under alternate modes
% of behavioral heritability, and alternate simulated weather conditions.
%=========================================================================
% Usage:
% out=hedgeSim4(reps, herit, cMu, cVar, b,d, weather, varBool, displayBool)
% Inputs:
%       * reps = number of times to repeat the simulation
%       * herit = {0,1}. 0=bet-hedging, 1=adaptive tracking.
%       * cMu = mean light choice probability
%       * cVar = variance in light choice probability
%       * b = model birth rate (beta)
%       * d = model death rate (delta)
%       * weather = struct of weather data, e.g. wp in hedgeWeatherPack.mat
%       * varBool = string of 3 boolean digits indicating which weather
%           conditions to simulate. 
%               '001' indicates only the mean daily temperatures
%               '011' indicates mean daily temperatures plus random daily deviations
%               '111' indicates daily means plus deviations plus random daily cloud cover
%       * displayBool = string of 4 boolean digits indicating which weather
%           conditions to simulate. 
%               'X---' plot or not the temperature history vector
%               '-X--' plot or not the population histories across reps
%               '--X-' plot or not the mean phototactic preference over time across all reps
%               '---X' display to the command window the current rep being simulated
% Output, a structure with these objects:
%       * .finalPops = vector of the final population of each of the rep simulations
%       * .popHistories = reps x seasonlength matrix of the population size vs time, for each rep
%       * .prefHistories = reps x seasonlength matrix of the mean phototactic preference vs time, for each rep 
%       * .data = simulation data for all flies simulated in the final rep
%       * .prefProfile = reps x seasonlength x 3 matrix of the 25th 50th 
%            and 75th percentile of fly phototactic preferences.

% num of flies at beginning of simulation
numFliesInit=100;

% unpack varBool for weather elements to simulate
vb=varBool;
varBool=[str2double(vb(1)) str2double(vb(2)) str2double(vb(3))];
cloudBool=varBool(1);
devBool=varBool(2);
seasonBool=varBool(3);

% unpack displayBool for display options
db=num2str(displayBool);
displayBool=[str2double(db(1)) str2double(db(2)) str2double(db(3)) str2double(db(4))];
graphTempHistBool=displayBool(1);
graphPopHistBool=displayBool(2);
graphPrefTrendBool=displayBool(3);
displayTick=displayBool(4);

%implement shade in the model?
shadeBool=1;

%double season length? If this is desired, set XL to 1.
XL=0;
if XL==1; daysToSim=length(weather.flySeasonInterval)*2; end;

% extract basic parameters
daysToSim=length(weather.flySeasonInterval);

%arbitrary value for fly initial health
flyHP=150;

%age at which a fly becomes fertile (days)
matureAge=10;

%new names for birth and death parameters
dailyMateP=b;
deathRate=d;

%extract weather mean temperature and set shade/sun temp diff.
tempMean=mean(weather.aveT(weather.flySeasonInterval));
shadeDiff=7;

%new name for mu. 
prefMean=cMu;

% calculate A and B parameters of beta-function give mu and var of observed
% distribution of light-choice probabilities
A=cMu*((cMu*(1-cMu))/cVar - 1);
B=(1-cMu)*((cMu*(1-cMu))/cVar - 1);


%indices into fly data matrix for each of these per fly variables
pref=1;
HP=2;
alive=3;
matur=4;
eggSun=5;
born=6;
dead=7;
heritPref=8;
howDead=10;

%initialize history vectors
runHists0=zeros(reps,daysToSim);
prefHists=zeros(reps,daysToSim);

% pick a random slightly desaturated color for plotting this rep
randColor=[rand() rand() rand()]*0.8;

% determine autoregression parameters for historical daily temp deviation data
ARparams=30;
[ar_coeffs,nV,reflect_coeffs] = aryule(weather.dev,ARparams);

%initialize output matrix for the 25th 50th and 75th preference percentile values
prefProfiles=zeros(reps,daysToSim,3);

% repeat simulation reps times
for num=1:reps
    if displayTick==1; disp(num); end; %if display rep number, do so
    
    tempHist=zeros(daysToSim,1)+tempMean; %set the season temperature history to a constant value, by default
    
    if seasonBool==1 %if seasonality is to be simulated, set the daily temperatures to the normals
        tempHist=weather.normal(weather.flySeasonInterval);
        if XL==1; tempHist=[tempHist;tempHist]; end; %concatenate if this is a 2x season simulation
    end
    
    if devBool==1  % if daily temperature deviations are to be simulated, calculate them
        devHist=filter(1,ar_coeffs,randn(daysToSim+100,1)*nV); %filter gaussian white noise with auto-regression parameters
        devHist=devHist(101:end);                              %throw out burn-in time series data
        tempHist=tempHist+devHist*weather.devStd/std(devHist); %add these deviations to the tempHistory vector
    end
    
    if cloudBool==1 % if clouds are to be simulated (this option was not used in the final paper)
        cloudHist=hedgeMakeCloud2(weather.cloud(weather.flySeasonInterval),0:10); % call cloud-making function
        if XL==1; cloudHist=[cloudHist;cloudHist]; end;  %concatenate if this is a 2x season simulation
    else
        cloudHist=zeros(daysToSim,1)+1; %otherwise, cloud cover vector is ones.
        if XL==1; cloudHist=[cloudHist;cloudHist]; end;
    end
    
    %if this is the first simulation, and the populations are to be
    %displayed, start the figure
    if num==1; if graphTempHistBool==1; figure; hold on; plot(tempHist); plot(tempHist+shadeDiff*cloudHist,'r'); drawnow; figure; hold on; end;  end;
    
    %initialize a 100x10 matri for data characterizing the seed population
    flyData=zeros(numFliesInit,10);
    
    %initialize the seed population
    for i=1:numFliesInit
        randHP=flyHP*rand(); %random possible age, gets max'ed in two lines against matureAge, so that all flies are post-eclosion at the start of the season, given that they over winter as adults
        prefTemp=betarnd(A,B); %draw from a fit beta distribution a random preference for this fly.
        % these values initialize the fly's photopref, hit points, living
        % status, age, whether it was laid as an egg in the sun, birth date
        % (NA, since seed), death date (TBD), photopref inherited from
        % parent (NA, encoded as own pref), index, and deadness.
        % these correspond to the indices listed above.
        flyData(i,:)=[prefTemp, randHP, 1, max([matureAge flyHP-randHP]), rand()<prefMean, NaN, NaN, prefTemp, i, 0]; % use min([matureAge flyHP-randHP]) for overwintering as eggs and/or pupae
    end
    
    %initialize history data for this rep
    popHist=zeros(daysToSim,1)+NaN;
    prefHist=zeros(daysToSim,1)+NaN;
    
    %run simulation
    for i=1:daysToSim
        
        %if the population has exploded, abort
        if sum(flyData(:,alive))>700
            break;
        end;

        %for all the flies in the simulation
        numFlies=size(flyData,1);
        
        for j=1:numFlies
            if flyData(j,alive)==1                       % if fly under consideration alive
                if rand() < deathRate                    % if it to be killed, update the data matrix
                    flyData(j,alive)=0;
                    flyData(j,dead)=i;
                    flyData(j,howDead)=1;
                else                                        % otherwise
                    if flyData(j,matur) < matureAge         % is the animal eclosed?
                        if flyData(j,eggSun) == 0           % if no, calculate temperature experienced according to egg/larva/pupal sun position
                            x=tempHist(i);
                        else
                            x=tempHist(i)+shadeDiff*cloudHist(i)*shadeBool;
                        end
                        %this temp contributes to its maturation age
                        maturHit=matureAge/(0.2306*x^2-11.828*x+158.34+1);
                        % update maturation age in data matrix.
                        flyData(j,matur)=flyData(j,matur)+maturHit;       
                    else %if eclosed already
                         % calculate temp perceived, based on photopreference.
                        x=tempHist(i)+flyData(j,pref)*shadeDiff*cloudHist(i)*shadeBool;
                        
                        if rand()<dailyMateP                        % mate today?
                             %determine photopreference of progeny animal
                            progenyPref=betarnd(A,B)*(1-herit)+flyData(j,heritPref)*herit;
                            % add a row to the data matrix for the new fly,
                            % appropriately initialized.
                            flyData=[flyData;[progenyPref flyHP 1 0 rand()<flyData(j,pref) i NaN flyData(j,heritPref) j 0]];
                        end
                    end
                    
                    %for all flies, eclosed or not, subtract HP for
                    %temp-dependent aging rate
                    HPHit=flyHP/(0.4074*x^2-28.356*x+506.2);
                    flyData(j,HP)=flyData(j,HP)-HPHit;        % update HP
                    if flyData(j,HP) <0                      % if negative HP
                        flyData(j,alive)=0;                     % mark as dead
                        flyData(j,dead)=i;
                        flyData(j,howDead)=2;
                    end
                    
                end
            end
        end
        %save stats on this day's flies
        popHist(i)=sum(flyData(:,alive));
        prefHist(i)=mean(flyData(flyData(:,alive)==1,pref));
        prefProfiles(num,i,:)=prctile(flyData(flyData(:,alive)==1,pref),[25 50 75]);
    end
    
    %record this season's data into output matrices.
    runHists0(num,:)=popHist;
    prefHists(num,:)=prefHist;
    
    % if season populations are being plotted, plot now.
    if graphPopHistBool==1; plot(popHist,'Color', randColor); ylim([0 400]); drawnow; end;
    
end

% output final variables
out.finalPops=runHists0(:,end);
out.popHistories=runHists0;
out.prefHistories=prefHists;
out.data=flyData;
out.prefProfile=reshape(mean(prefProfiles,1),daysToSim,3);

% if photo preference trend is to be plotted, plot it here
if graphPrefTrendBool==1; figure; errorbar(nanmean(prefHists),nanstd(prefHists)/sqrt(reps),'g'); end;