function summary = rangeshift(rowmax, colmax, Ninit, Rmax, K, C, muRate_D, muSD_D, muRate_H, muSD_H, ... gmaxB, gmaxT, avshift, savefiles, blowup, modeltype, platmin, platmax, gmargin1, margsize) % INPUTS % rowmax = number of rows of patches in suitable habitat % colmax = number of columns of patches in suitable habitat % Ninit = initial population size (typically set to rowmax*colmax*K) % Rmax = mean maximum number of offspring per individual per generation % K = carrying capacity per patch % C = cost of dispersal (probability of death) % muRate_D = probability of mutation to dispersal probability % muSD_D = standard deviation of mutation to dispersal probability % muRate_H = probability of mutation to environmental optimum % muSD_H = standard deviation of mutation to environmental optimum % gmaxB = number of generations of equilibration before climate change starts % gmaxT = total number of generations % avshift = mean number of columns climate shifts per generation after generation gmaxB % savefiles = a switch: 1 means save data, figures, etc.; 0 means don't <1 or 0> % blowup = factor by which map images are 'blown up' (otherwise they are very small) % modeltype = a switch: 1 means there is no previous local adaptation, 2 or 3 means there is <1, 2, or 3> %% NOTE that only model types 1 and 2 are explored in Hargreaves et al 2015 % platmin = for modeltype 1 or 2, the minimum column where relative fitness equals 1 %% NOTE that platmin is equal to margin+1 in Hargreaves et al 2015 % platmax = for modeltype 1 or 2, the maximum column where relative fitness equals 1 % gmargin1 = for modeltype 3, half the range size +1 of an individual genotype. if less than platmin, will result in genotypes locally adapted to marginal habitat %% NOTE that gmargin1(from this code) is equal to g.margin+1 from Hargreaves et al 2015 % margsize = columns averaged across to get D lines in summary figure at the lagging and leading edges of the suitable habitat % % OUTPUT % summary = matrix with gmaxT + 1 rows and 3*colmax + 2 columns(description below) % ** Column Description: ** % column 1 = generation (0 to gmaxT) % columns 2 to colmax + 1 = average density by column of suitable range % columns colmax + 2 to 2*colmax + 1 = average dispersal probability (D) by column % columns 2*colmax + 2 to 3*colmax + 1 = average environmental optimum (H) by column % column 3*colmax + 2 = realized size of range shift % % NOTE ON SAVED FILES % If 'savefiles' is set to 1, the function saves files to sub-directories % of c:\rangeshift\ (~/Documents/MATLAB/rangeshift/ for macs) % _________________________________________________________________________ % (0) PREAMBLE % Update history % 2015-06-29 by Susan (cleaning up and further annotation for Dryad submission) % 2014-12-15 by Rob (so will work again when modeltype == 2) % 2014-11-29 by Anna (switches so code will work when modeltype==1) % 2014-11-24 altered H code (H genotypes start locally adapted) % 2014-11-21 mutation in H by Rob % 2014-03-16 annotation added by Susan % 2014-01-23 by Susan (annotated AH Jan 27) % 2013-09-24 by Susan % 2013-06-11 by Rob % 2013-04-16 by Rob % 2013-04-01 by Rob % 2013-03-12 by Rob % 2013-03-07 by Susan % 2013-03-05 by Rob % Example values for user-defined parameters % rowmax = 20; % number of rows in the climate window % colmax = 101; % number of cols in the climate window % Ninit = 80800; % initial population size (actually a variable, not a parameter) % Rmax = 5; % maximum number of offspring per individual % K = 40; % carrying capacity per cell % C = 0.1; % dispersal cost % muRate_D = 0.05; % mutation rate of dispersal probability % muSD_D = 0.05; % standard deviation of mutation size for dispersal probability % muRate_H = 0.005; % mutation rate of environmental optimum % muSD_H = 0.05; % standard deviation of mutation size for environmental optimum % gmaxB = 1000; % number of generations before climate change starts % gmaxT = 4000; % total number of generations investigated (must be greater than gmaxB) % avshift = 0.1; % average distance of shift during climate change % savefiles = 1; % save output files (1) or not (0) % blowup = 40; % factor by which map images are 'blown up' (otherwise they are very small) % modeltype = 1, 2, or 3 % model version (determines whether there is local adaptation or not, and if so, how it is implemented) % platmin = 21; % minimum column where relative fitness equals 1 % platmax = 81; % maximum column where relative fitness equals 1 % gmargin1 = 25; % half the width +1 of a genoytpe's range (typically <= margin) % margsize = 5; % columns averaged across to get D lines in summary fig % _________________________________________________________________________ % (1) CHECK INPUTS FOR ERRORS AND WARNINGS if mod(rowmax, 1) ~= 0 || rowmax < 1; error('rowmax must be a positive integer'); end if mod(colmax, 1) ~= 0 || colmax < 1; error('colmax must be a positive integer'); end if mod(Ninit, 1) ~= 0 || Ninit < 1; error('Ninit must be a positive integer'); end if Ninit ~= rowmax*colmax*K; warning('WarnTests:Ninit','warning: Ninit is typically set to rowmax*colmax*K'); end if Rmax <= 0; error('Rmax must be positive'); end if K <= 0; error('K must be positive'); end if C < 0 || C > 1; error('C must be between 0 and 1, inclusive'); end if muRate_D < 0 || muRate_D > 1; error('muRate_D must be between 0 and 1, inclusive'); end if muSD_D < 0; error('muSD_D must be non-negative (and should be small, preferably'); end if muRate_H < 0 || muRate_D > 1; error('muRate_H must be between 0 and 1, inclusive'); end if muSD_H < 0; error('muSD_H must be non-negative (and should be small, preferably'); end if mod(gmaxB, 1) ~= 0 || gmaxB < 1; error('gmaxB must be a positive integer'); end if mod(gmaxT, 1) ~= 0 || gmaxT < 1 || gmaxT < gmaxB; error('gmaxT must be a positive integer greater than or equal to gmaxB'); end if avshift < 0; error('avshift must be non-negative'); end if savefiles ~= 0 && savefiles ~= 1; error('savefiles must be 0 or 1'); end if mod(blowup, 1) ~=0 || blowup < 1; error('blowup must be a positive integer'); end if ~ismember(modeltype, [1 2 3]); error('modeltype must be 1, 2, or 3'); end if ~ismember(platmin, (1:colmax)) || platmin > platmax; error('platmin must be a positive integer between 1 and colmax, and must be less than or equal to platmax'); end if ~ismember(platmax, (1:colmax)) || platmax < platmin; error('platmax must be a positive integer between 1 and colmax, and must be greater than or equal to platmin'); end if ~ismember(gmargin1, (1:colmax)) || gmargin1 > platmin; error('gmargin1 must be a positive integer between 1 and colmax, and must be less than or equal to platmin'); end if ~ismember(margsize,(1:floor(colmax/2))); error('margsize must be a positive integer between 1 and colmax/2 (rounded down)'); end % _________________________________________________________________________ % (2) SET-UP AND INITIAL CONDITIONS % generate timestamp, unique for each individual run (don't bother if not saving files) if savefiles st = clock; year = num2str(st(1)); if st(2) >= 10; month = num2str(st(2)); else month = ['0', num2str(st(2))]; end if st(3) >= 10; day = num2str(st(3)); else day = ['0', num2str(st(3))]; end if st(4) >= 10; hour = num2str(st(4)); else hour = ['0', num2str(st(4))]; end if st(5) >= 10; minute = num2str(st(5)); else minute = ['0', num2str(st(5))]; end timestamp = [year, month, day, hour, minute]; end % create and simultaneously initialize N, inds, densmap, and Dmap N = Ninit; % create N, the current population size inds = zeros(N, 4); % create inds, which tracks the traits of the individiuals densmap = zeros(rowmax, colmax); % create densmap, a map of local population density Dmap = zeros(rowmax, colmax); % create Dmap, a map of mean dispersal probabilities for i = 1:N row = ceil(rand*rowmax); inds(i, 1) = row; % row of individual i col = ceil(rand*colmax); inds(i, 2) = col; % col of individual i D = rand; inds(i, 3) = D; % dispersal Pr of individual i if (modeltype == 2 || modeltype == 3) % start all individuals at their best possible location. If gmargin1 colmax - gmargin1 + 1 % col in far leading edge to which no genotype is optimally adapted H = colmax - gmargin1 + 1; else H = col; end inds(i, 4) = H; % optimum column of individual i Hmap(row, col) = Hmap(row, col) + H; %increase running total of H at pos. (row, col) Hmap = Hmap./densmap; % change Hmap to a mean by dividing by local density % if there are 0 individuals in a particular grid cell there is no mean % environmental optimum, so the entry in Hmap is 'NaN' (i.e., 0/0), not zero. end densmap(row, col) = densmap(row, col) + 1; % increment density at pos. (row, col) Dmap(row, col) = Dmap(row, col) + D; % increase running total of D at pos. (row, col) end Dmap = Dmap./densmap; % change Dmap to a mean by dividing by local density % if there are 0 individuals in a particular grid cell there is no mean % dispersal probility, so the entry in Dmap is 'NaN' (i.e., 0/0), not zero. % create and initialize densmapinit, Dmapinit and Hmapinit densmapinit = densmap; % create densmapinit (used to save image of distribution of original density) Dmapinit = Dmap; % create Dmapinit (used to save image of distribution of original dispersal probability) if (modeltype == 2 || modeltype == 3) Hmapinit = Hmap; % create Hmapinit (used to save image of distribution of original environmental optimum) end % determine column mean of Dmap and Hmap (weighted by density) meanDmap = zeros(1, colmax); % initialize meanDmap. if (modeltype == 2 || modeltype == 3) meanHmap = zeros(1, colmax); % initialize meanHmap. end for i = 1:N meanDmap(inds(i, 2)) = meanDmap(inds(i, 2)) + inds(i, 3); if (modeltype == 2 || modeltype == 3) meanHmap(inds(i, 2)) = meanHmap(inds(i, 2)) + inds(i, 4); end end for j = 1:colmax if sum(densmap(:, j)) == 0 % if there are no individuals in column j meanDmap(j) = NaN; if (modeltype == 2 || modeltype == 3) meanHmap(j) = NaN; end else meanDmap(j) = meanDmap(j)/sum(densmap(:, j)); if (modeltype == 2 || modeltype == 3) meanHmap(j) = meanHmap(j)/sum(densmap(:, j)); end end end % record individuals information at start of model run inds_init = inds; % create and initialize Hmapinit if (modeltype == 2 || modeltype == 3) Hmapinit = Hmap; % create Hmapinit (used to save image of distribution of original environmental optimum) % determine column mean of Hmap (weighted by density) meanHmap = zeros(1, colmax); % initialize meanHmap. for i = 1:N meanHmap(inds(i, 2)) = meanHmap(inds(i, 2)) + inds(i, 4); end for j = 1:colmax if sum(densmap(:, j)) == 0 % if there are no individuals in column j meanHmap(j) = NaN; else meanHmap(j) = meanHmap(j)/sum(densmap(:, j)); end end end % create summary matrix to record data if (modeltype == 1) summaryinit = [0, mean(densmap), meanDmap, 0]; % first row of summary (to be concatenated with summary at end)) summary = zeros(gmaxT, 2*colmax + 2); % records generation, mean density by row, mean dispersal Pr by row, mean H by row, size of range shift else summaryinit = [0, mean(densmap), meanDmap, meanHmap, 0]; % first row of summary (to be concatenated with summary at end)) summary = zeros(gmaxT, 3*colmax + 2); % records generation, mean density by row, mean dispersal Pr by row, mean H by row, size of range shift end % _________________________________________________________________________ % (3) THE MODEL for g = 1:gmaxT rs = 0; % the distance the range is shifted to the right (modified if climate change occurs) % climate change if g > gmaxB % climate chante occurs after transitory period of gmaxB generations % determine rs, the distance the range is shifted to the right (Poisson process) rs = poissrnd1(avshift); % poissrnd1 is a custom function (below) if rs > 0 % if a range shift occurs % range shift right of rs is equivalent to position shift left of rs of every individual inds(:, 2) = inds(:, 2) - rs; % determine who is still within a suitable climate indsnew = zeros(N, 4); % the new list of individuals remaining in the population j = 0; % number of individuals remaining in the population; also index of indsnew for i = 1:N col = inds(i, 2); if col >= 1 j = j + 1; % increment the number of individuals remaining in the population indsnew(j, :) = inds(i, :); % include the surviving individual to the new population list end end N = j; % update the population size inds = indsnew(1:N, :); % inds now represents the individuals that have survived climate change (truncated to the number that survived) % update densmap (because it is needed for calculating density-dependent reproduction) densmap = [densmap(:, (rs + 1):colmax), zeros(rowmax, rs)]; % no need to update Dmap or Hmap at this stage end end % create next generation offs = zeros(N*Rmax, 4); % traits of offspring (max num rows is N*Rmax; will be truncated) j = 0; % total number of new offspring; also index of offs for i = 1:N % determine the traits of individual i row = inds(i, 1); % row of individual i col = inds(i, 2); % col of individual i D = inds(i, 3); % D of individual i H = inds(i, 4); % H of individual i % determine w, the relative fitness of individual i if modeltype == 1 if col < platmin % marginal habitat near trailing edge w = col/platmin; elseif col > platmax % marginal habitat near leading edge w = (colmax - col + 1)/(colmax - platmax + 1); else % i.e., if platmin <= col <= platmax; optimal habitat near range centre w = 1; end elseif modeltype == 2 dist = abs(H - col); % number of columns away from optimum w = (gmargin1 - dist)/gmargin1; % if dist >= gmargin1 fitness = 0 if gmargin1 colmax - platmin + 1 % genotypes adapted to marginal habitat @ leading edge w = ((gmargin1 - dist)/gmargin1)*(1+(colmax - platmin + 1-H)/platmin); end end if w < 0 w = 0; % minimum relative fitness can be is 0 end elseif modeltype == 3 dist = abs(H - col); % number of columns away from optimum midpoint = (colmax + 1)/2; % midpoint of the distribution if H <= midpoint halfwidth = H; % half the width of the fitness distribution else halfwidth = colmax - H + 1; end w = (halfwidth - dist)/halfwidth; if w < 0 w = 0; % minimum relative fitness can be is 0 end end % determine k, the number of offspring of individual i (Poisson process) lambda = w*Rmax/(1 + densmap(row, col)/K*(Rmax - 1)); k = poissrnd1(lambda); % poissrnd1 is a custom function (below) % offspring are born for m = 1:k j = j + 1; % increment total number of offspring offs(j, 1) = row; % offspring retains parental row (for now) offs(j, 2) = col; % offspring retains parental col (for now) offs(j, 3) = D; % offspring retains parental D (for now) offs(j, 4) = H; % offspring retains parental H % mutation to disperal probability if rand < muRate_D newD = D + muSD_D*randn; % potential new D for offspring while newD < 0 || newD > 1 % throw out invalid values of D newD = D + muSD_D*randn; % potential new D for offspring end offs(j, 3) = newD; % offspring has new, mutated D end % mutation to environmental optimum if (modeltype==2 || modeltype==3) && rand < muRate_H newH = round(H + colmax*muSD_H*randn); % potential new H for offspring while newH < gmargin1 || newH > colmax - gmargin1 + 1; % throw out invalid values of H newH = round(H + colmax*muSD_H*randn); % potential new H for offspring end offs(j, 4) = newH; end end end % parents die, offspring remain N = j; % new population size is equal to number of offspring inds = offs(1:N, :); % offspring comprise the new population (truncate to N rows) % dispersal in a cylindrical world (i.e., top and bottom rows are adjacent) % simultaneously determine new densmap, Dmap, Hmap, and inds densmap = zeros(rowmax, colmax); Dmap = zeros(rowmax, colmax); if (modeltype==2 || modeltype==3) Hmap = zeros(rowmax, colmax); end j = 0; % count the number of individuals that do not disperse out of suitable range; also index of indsnew indsnew = zeros(N, 4); % keep track of the individuals that do not disperse out of suitable range (max num rows is N; will be truncated) for i = 1:N row = inds(i, 1); % row of individual i col = inds(i, 2); % col of individual i D = inds(i, 3); % dispersal Pr of i ddd = 0; % death due to dispersal (hence 'ddd'); specifically due to cost of dispersal, not leaving suitable range (altered if death occurs) if rand < D % if dispersal occurs dir = ceil(rand*4); % direction of dispersal is random (1 = U, 2 = D, 3 = L, 4 = R) if dir == 1 if row == 1 % if i is in the top row row = rowmax; % wrap i to rowmax; column stays the same %row = row; else % otherwise row = row - 1; % move i one row up %row = row; end elseif dir == 2 if row == rowmax % if i is in the bottom row row = 1; % wrap i to 1st row; column stays the same %row = row; else % otherwise row = row + 1; % move i one row down %row = row; end elseif dir == 3 col = col - 1; % move i one column left (incl. to col 0) %col = col; else col = col + 1; % move i one column right (incl. to col colmax + 1) %col = col; end if rand < C % if death occurs during dispersal ddd = 1; % set ddd to 1 end end inds(i, 1) = row; % update i's row inds(i, 2) = col; % update i's col if col >= 1 && col <= colmax && ddd == 0 % count individuals that have not dispersed out of suitable range, nor died during dispersal j = j + 1; % increment the number of individuals remaining indsnew(j, :) = inds(i, :); % record individuals that remain in suitable range densmap(row, col) = densmap(row, col) + 1; % increment density at pos. (row, col) Dmap(row, col) = Dmap(row, col) + D; % increase running total of D at pos. (row, col) if (modeltype==2 || modeltype==3) Hmap(row, col) = Hmap(row, col) + H; % increase running total of H at pos. (row, col) end end end Dmap = Dmap./densmap; % change Dmap to a mean by dividing by local density % if there are 0 individuals in a particular grid cell there is no mean % dispersal probility, so the entry in Dmap is 'NaN' (i.e., 0/0), not zero. if (modeltype==2 || modeltype==3) Hmap = Hmap./densmap; % change Hmap to a mean by dividing by local density end % if there are 0 individuals in a particular grid cell there is no mean % environmental optimum, so the entry in Hmap is 'NaN' (i.e., 0/0), not zero. % determine new population (i.e., those that haven't dispersed outside suitable range) N = j; % determine new population size inds = indsnew(1:N, :); % update inds (truncate to N rows) % determine column mean of Dmap and Hmap (weighted by density) meanDmap = zeros(1, colmax); % initialize meanDmap if (modeltype == 2 || modeltype == 3) meanHmap = zeros(1, colmax); % initialize meanHmap end for i = 1:N meanDmap(inds(i, 2)) = meanDmap(inds(i, 2)) + inds(i, 3); if (modeltype == 2 || modeltype == 3) meanHmap(inds(i, 2)) = meanHmap(inds(i, 2)) + inds(i, 4); end end for j = 1:colmax if sum(densmap(:, j)) == 0 % if there are no individuals in column j meanDmap(j) = NaN; if (modeltype == 2 || modeltype == 3) meanHmap(j) = NaN; end else meanDmap(j) = meanDmap(j)/sum(densmap(:, j)); if (modeltype == 2 || modeltype == 3) meanHmap(j) = meanHmap(j)/sum(densmap(:, j)); end end end % determine column mean of Hmap (weighted by density) for model2 & 3 if (modeltype==2 || modeltype==3) meanHmap = zeros(1, colmax); % initialize meanHmap for i = 1:N meanHmap(inds(i, 2)) = meanHmap(inds(i, 2)) + inds(i, 4); end for j = 1:colmax if sum(densmap(:, j)) == 0 % if there are no individuals in column j meanHmap(j) = NaN; else meanHmap(j) = meanHmap(j)/sum(densmap(:, j)); end end end % record data if (modeltype==1) summary(g, :) = [g, mean(densmap), meanDmap, rs]; % add row of summary else summary(g, :) = [g, mean(densmap), meanDmap, meanHmap, rs]; % add row of summary end % record state of 'maps' at critical generations (used to save 'snapshots' of the state of the population) if g == gmaxB densmapeq = densmap; % record densmapeq, the distribution of density at 'equilibrium' Dmapeq = Dmap; % record Dmapeq, the distribution of dispersal probability at 'equilibrium' if (modeltype==2 || modeltype==3) Hmapeq = Hmap; % record Hmapeq, the distribution of environmental optimum at 'equilibrium' end inds_eq = inds; % record inds_eq, information about all individuals at 'equilibrium' elseif g == gmaxT densmapfinal = densmap; % record densmapfinal, the distribution of density at the end of the model run Dmapfinal = Dmap; % record Dmapfinal, the distribution of dispersal probability at the end of the model run if (modeltype==2 || modeltype==3) Hmapfinal = Hmap; % record Hmapfinal, the distribution of environmental optimum at the end of the model run end inds_final = inds; % record inds_final, information about all individuals at the end of the model run end end summary = [summaryinit; summary]; % concatenate summaryinit and summary % _________________________________________________________________________ % (4) CREATE ON-SCREEN FIGURE DTS = summary(:, (colmax + 2):(2*colmax + 1)); % dispersal time series (DTS) DTSIm = zeros(gmaxT + 1, colmax, 3); % image of DTS for i = 1:gmaxT + 1 for j = 1:colmax if isnan(DTS(i, j)) DTSIm(i, j, :) = [1 0 0.5]; % cells with no individuals are pink else s = 1 - DTS(i, j); % shade of cell DTSIm(i, j, :) = [s s s]; end end end axes('position', [0.1 0.35 0.2 0.55]); image(DTSIm); hold on; plot([0; colmax], [gmaxB; gmaxB], '--', 'color', [0 0.5 0], 'linewidth', 1.5); hold off; ylabel('Model generations'); axis([1 colmax 0 gmaxT]); axes('position', [0.4 0.35 0.2 0.55]); plot([0; colmax], [gmaxB; gmaxB], '--', 'color', [0 0.5 0], 'linewidth', 1.5); hold on; for i = 0:(margsize - 1) plot(summary(:, 2*colmax + 1 - i), (0:gmaxT)', 'color', 'black', 'linewidth', 1.5); set(gca, 'YDir', 'reverse'); % leading edge (black) plot(summary(:, colmax + 2 + i), (0:gmaxT)', 'color', 'red', 'linewidth', 1.5); set(gca, 'YDir', 'reverse'); % trailing edge (red) end hold off; axis([0 1.05*max(max([summary(:, 2*colmax + 1), summary(:, colmax + 2)])) 0 gmaxT]); xlabel('Mean D'); if modeltype > 1 axes('position', [0.4 0.1 0.2 0.2]); plot((1:colmax)', 1-abs(summary(gmaxB + 1, (2*colmax + 2):(3*colmax + 1))-(1:colmax))/platmin, 'color', [0 0 1], 'linewidth', 1.5); hold on; plot((1:colmax)', 1-abs(summary(gmaxT + 1, (2*colmax + 2):(3*colmax + 1))-(1:colmax))/platmin, 'color', [1 0.5 0.5], 'linewidth', 1.5); hold off; xlabel('x position'); ylabel('1-abs(mean H - optimum)'); % i.e. how far from the optimum, on average, are the inhabitants of each column axis([1 colmax 0 1]) end axes('position', [0.1 0.1 0.2 0.2]); plot((1:colmax)', summary(gmaxB + 1, (colmax + 2):(2*colmax + 1)), 'color', [0 0 1], 'linewidth', 1.5); hold on; plot((1:colmax)', summary(gmaxT + 1, (colmax + 2):(2*colmax + 1)), 'color', [1 0.5 0.5], 'linewidth', 1.5); hold off; xlabel('x position'); ylabel('Mean D'); axis([1 colmax 0.95*min(min([summary(gmaxB + 1, (colmax + 2):(2*colmax + 1)), summary(gmaxT + 1, (colmax + 2):(2*colmax + 1))])) ... 1.05*max(max([summary(gmaxB + 1, (colmax + 2):(2*colmax + 1)), summary(gmaxT + 1, (colmax + 2):(2*colmax + 1))]))]); axes('position', [0.7 0.1 0.2 0.8]); axis off; text(0, 1.0, 'Parameter values', 'fontweight', 'bold'); text(0, 1-1/14, ['rowmax = ', num2str(rowmax)]); text(0, 1-2/14, ['colmax = ', num2str(colmax)]); text(0, 1-3/14, ['Ninit = ', num2str(Ninit)]); text(0, 1-4/14, ['Rmax = ', num2str(Rmax)]); text(0, 1-5/14, ['K = ', num2str(K)]); text(0, 1-6/14, ['C = ', num2str(C)]); text(0, 1-7/14, ['muRate_D = ', num2str(muRate_D)]); text(0, 1-8/14, ['muSD_D = ', num2str(muSD_D)]); text(0, 1-9/14, ['avshift = ', num2str(avshift)]); text(0, 1-10/14, ['platmin = ', num2str(platmin)]); text(0, 1-11/14, ['platmax = ', num2str(platmax)]); text(0, 1-12/14, ['modeltype = ', num2str(modeltype)]); if (modeltype==2 || modeltype==3) text(0, 1-13/14, ['muRate_H = ', num2str(muRate_H)]); text(0, 0, ['muSD_H = ', num2str(muSD_H)]); text(0, -1/14, ['gmargin1 = ', num2str(gmargin1)]); else text(0, 1-13/14, 'muRate_H = NA'); text(0, 0, 'muSD_H = NA'); text(0, -1/14, 'gmargin1 = NA'); end % _________________________________________________________________________ % (5) DATA PROCESSING AND FILE OUTPUT if savefiles % make directory in which to deposit data mkdir(['~/Documents/MATLAB/rangeshift/', timestamp]); % save txt file of parameter values params = [rowmax; colmax; Ninit; Rmax; K; C; muRate_D; muSD_D; muRate_H; muSD_H; gmaxB; gmaxT; avshift; modeltype; platmin; platmax]; filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 00. PARAMS.txt']; save(filename, 'params', '-ascii'); % save txt file of summary filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 01. SUMMARY.txt']; save(filename, 'summary', '-ascii'); % save .tif file of summary figure filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 02. SUMMARY.tif']; saveas(gcf, filename, 'tif'); % save .txt files of maps from various points of time filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 03. DENSMAP_INIT.txt']; save(filename, 'densmapinit', '-ascii'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 04. DENSMAP_EQ.txt']; save(filename, 'densmapeq', '-ascii'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 05. DENSMAP_FINAL.txt']; save(filename, 'densmapfinal', '-ascii'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 06. DMAP_INIT.txt']; save(filename, 'Dmapinit', '-ascii'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 07. DMAP_EQ.txt']; save(filename, 'Dmapeq', '-ascii'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 08. DMAP_FINAL.txt']; save(filename, 'Dmapfinal', '-ascii'); if(modeltype==2 || modeltype==3) filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 09. HMAP_INIT.txt']; save(filename, 'Hmapinit', '-ascii'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 10. HMAP_EQ.txt']; save(filename, 'Hmapeq', '-ascii'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 11. HMAP_FINAL.txt']; save(filename, 'Hmapfinal', '-ascii'); end %filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 21. INDS_INIT.txt']; save(filename, 'inds_init', '-ascii'); %filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 22. INDS_EQ.txt']; save(filename, 'inds_eq', '-ascii'); %filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 23. INDS_FINAL.txt']; save(filename, 'inds_final', '-ascii'); % save images of maps from various points of time densmapinitIm = zeros(rowmax*blowup, colmax*blowup, 3); % image of densmapinit densmapeqIm = zeros(rowmax*blowup, colmax*blowup, 3); % image of densmapeq densmapfinalIm = zeros(rowmax*blowup, colmax*blowup, 3); % image of densmapfinal DmapinitIm = zeros(rowmax*blowup, colmax*blowup, 3); % image of Dmapinit DmapeqIm = zeros(rowmax*blowup, colmax*blowup, 3); % image of Dmapeq DmapfinalIm = zeros(rowmax*blowup, colmax*blowup, 3); % image of Dmapfinal if(modeltype==2 || modeltype==3) HmapinitIm = zeros(rowmax*blowup, colmax*blowup, 3); % image of Hmapinit HmapeqIm = zeros(rowmax*blowup, colmax*blowup, 3); % image of Hmapeq HmapfinalIm = zeros(rowmax*blowup, colmax*blowup, 3); % image of Hmapfinal end for i = 1:rowmax for j = 1:colmax for ii = 1:blowup for jj = 1:blowup if densmapinit(i, j) < 0.5*K densmapinitIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [0 0 0]; % black elseif densmapinit(i, j) < 0.75*K densmapinitIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [0 0 1]; % blue elseif densmapinit(i, j) < 0.95*K densmapinitIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [0 1 0]; % green elseif densmapinit(i, j) < 1.05*K densmapinitIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 1 0]; % yellow elseif densmapinit(i, j) < 1.25*K densmapinitIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 0.5 0]; % orange elseif densmapinit(i, j) < 1.5*K densmapinitIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 0 0]; % red else densmapinitIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 1 1]; % white end if densmapeq(i, j) < 0.5*K densmapeqIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [0 0 0]; % black elseif densmapeq(i, j) < 0.75*K densmapeqIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [0 0 1]; % blue elseif densmapeq(i, j) < 0.95*K densmapeqIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [0 1 0]; % green elseif densmapeq(i, j) < 1.05*K densmapeqIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 1 0]; % yellow elseif densmapeq(i, j) < 1.25*K densmapeqIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 0.5 0]; % orange elseif densmapeq(i, j) < 1.5*K densmapeqIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 0 0]; % red else densmapeqIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 1 1]; % white end if densmapfinal(i, j) < 0.5*K densmapfinalIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [0 0 0]; % black elseif densmapfinal(i, j) < 0.75*K densmapfinalIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [0 0 1]; % blue elseif densmapfinal(i, j) < 0.95*K densmapfinalIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [0 1 0]; % green elseif densmapfinal(i, j) < 1.05*K densmapfinalIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 1 0]; % yellow elseif densmapfinal(i, j) < 1.25*K densmapfinalIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 0.5 0]; % orange elseif densmapfinal(i, j) < 1.5*K densmapfinalIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 0 0]; % red else densmapfinalIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [1 1 1]; % white end s = 1 - Dmapinit(i,j); DmapinitIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [s s s]; % grey scale s = 1 - Dmapeq(i,j); DmapeqIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [s s s]; % grey scale s = 1 - Dmapfinal(i,j); DmapfinalIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [s s s]; % grey scale if(modeltype==2 || modeltype==3) s = (colmax - Hmapinit(i,j))/(colmax - 1); HmapinitIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [s s s]; % grey scale s = (colmax - Hmapeq(i,j))/(colmax - 1); HmapeqIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [s s s]; % grey scale s = (colmax - Hmapfinal(i,j))/(colmax - 1); HmapfinalIm((i - 1)*blowup + ii, (j - 1)*blowup + jj, :) = [s s s]; % grey scale end end end end end filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 12. DENSMAP_INIT_IM.tif']; imwrite(densmapinitIm, filename, 'tif'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 13. DENSMAP_EQ_IM.tif']; imwrite(densmapeqIm, filename, 'tif'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 14. DENSMAP_FINAL_IM.tif']; imwrite(densmapfinalIm, filename, 'tif'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 15. DMAP_INIT_IM.tif']; imwrite(DmapinitIm, filename, 'tif'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 16. DMAP_EQ_IM.tif']; imwrite(DmapeqIm, filename, 'tif'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 17. DMAP_FINAL_IM.tif']; imwrite(DmapfinalIm, filename, 'tif'); if(modeltype==2 || modeltype==3) filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 18. HMAP_INIT_IM.tif']; imwrite(HmapinitIm, filename, 'tif'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 19. HMAP_EQ_IM.tif']; imwrite(HmapeqIm, filename, 'tif'); filename = ['~/Documents/MATLAB/rangeshift/', timestamp, '/', timestamp, ' 20. HMAP_FINAL_IM.tif']; imwrite(HmapfinalIm, filename, 'tif'); end end % _________________________________________________________________________ % (6) CUSTOM ANCILLARY FUNCTIONS function poissrnd1 = poissrnd1(lambda) % POISSRND1(LAMBDA) - alternative to poissrnd from statistics toolbox % % Draws a number from a Poisson distribution with mean of lambda (Knuth algorithm) k = 0; u = rand; while u >= exp(-lambda) u = u*rand; k = k + 1; end poissrnd1 = k; % _________________________________________________________________________ % (7) THE OFFICIAL MOTTO OF THE MODEL SHEEP % 0.05 < p < 'a stick in the eye'