In ToonVanDaele/metapop:

knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)

Een Individual based model

Dit script voegt een aantal R-scripts, R-functies en ook een Java methode samen om zo vanuit de input van een aantal parameters als output het gedrag van een metapopulatie te simuleren.

Input

In deze R-code pas je de populatie-parameters aan indien nodig.

seed = 420                        #Een seed voor random number generators
tmax = 1000                       #Total number of timesteps
nIters = 200                      #Total number of iterations of a whole loop
burnin.time = 15                  #Amount of time the population evolves before information is measured
nPop = 6                          #Number of populations
indiv.per.ha = 5                  #Amount of individuals that can live on 1 ha
nStages = 2                       #Amount of stages
stages.names = c("Juvenile","Adult")
maxAge = 20                       #Maximum age an individual can reach
quasiExtinctionBound = 10         #Below this bound (on total population) the population is seen a extinct
initial.agedistr = c(0,5)         #Age ditribution for the initial population
survivalChance = c(0.5,0.8)       #Survival chances for each stage
survivalVar = c(0.05,0)           #Variance of the survival chances
reproductiveChance = 0.5          #Probability that an idividual reproduces
reproductiveVar = 0               #Variance on reproductive probability
reproductive.age = 2              #Age/Stagenumber at which an individual can reproduce
nOffspring = 3                    # Mean nestsize (will be Poisson ditributed)
stochasticity = c(T,F,F)          #For which probabilities should stochastic values be calculated (survivalchances,reproduction)?
correlation = FALSE               #Should the parameters be correlated?      
correlation.onlyStoch = FALSE     #Should only the stochastic parameters be correlated or all parameters?

In dit stuk R-code kan je enkele parameters aanpassen die invloed hebben op de patches

mapsize = 1000    #Total size of the area in which the patches are created
dist_min = 200    #Minimum distance between the centre of two patches
area_M = 2.5      #Mean area of the patches
area_SD = 0.8     #Standard deviance of the patch areas
disp = 500        #Mean distance over which individuals are able to disperse
plotPatch = TRUE  #Should the patch configuration be plotted?
stay = 100        #If this number is small the individuals have a small migration probability

Hier kan je namen van files aanpassen waaring tussentijdse of finale informatie zal opgeslagen worden.

  output.patch = "patches.txt"                          #Where the patch info will be written to
  output.migration = "dispersion.txt"                   #Where the dispersion matrix will be written to
  output.stoch = "stochSurvival.txt"                    #Where the stochastic parameters will be written to
  output.population.evolution = "Evolution.txt"         #Where the results of a loop will be written to
  output.population.extinction = "ExtinctionTimes.txt"  #Where the extinction times will be written to

Geef hier aan wat voor run moet gebeuren: 1 run van r tmax tijdsstappen of 500 runs van r tmax tijdsstappen zodat ook extinctiekansen kunnen berekend worden.

oneLoop <- F
extinctionLoop <- T

Resultaten

Input genereren

Dit script genereert patches, de migratiekansen en stochastische waarden voor de survival- en reproductieparameters indien nodig.

source("../../R/ibm_input.R")

outputPatch(seed = seed, mapsize = mapsize, dist_min = dist_min, areaM = area_M, areaSD = area_SD,
            Npatch = nPop, disp = disp, plotG = plotPatch, stay = stay, output.patch = output.patch,
            output.migration = output.migration)
if (correlation.onlyStoch) {
  corrMat = corrMatrix(npar = length(survivalChance) + 1 , npop = nPop, correlation = correlation, stochasticity = stochasticity)

  }else{
  corrMat = corrMatrix(npar = length(survivalChance) + 1 , npop = nPop, correlation = correlation)
}
stochMatrix(seed = seed, tmax = tmax, Npop = nPop, stages = nStages, stages.names = stages.names,
            stochasticity = stochasticity, survivalM = survivalChance, survivalVAR = survivalVar,
            reproductionM = reproductiveChance, reproductionVAR = reproductiveVar,
            corrMat = corrMat, output.stoch = output.stoch)

Java Scripts runnen

vector.to.String <- function(v, sep=";"){
  vString = paste0("\"",v[1])
  for (i in 2:length(v)) {
    vString <- paste(vString,v[i],sep = sep)
  }
  vString <- paste0(vString,"\"")
  return(vString)
}

getwd()
setwd("../java/MetaPopulation")
system("javac -cp lib/jdistlib-0.4.4-bin.jar src/metapopulation/*.java -d build")

tmpfile = tempfile()

# Default argument
command <- paste("java -cp bin;lib/jdistlib-0.4.4-bin.jar metapopulation.MetaPopulation")
command
system(command)

# Textfile with arguments
command <- paste("java -cp bin;lib/jdistlib-0.4.4-bin.jar metapopulation.MetaPopulation",
                 "input.txt", 
                 ">", tmpfile)
command
system(command)

# Full argument list
command <- paste("java -cp bin;lib/jdistlib-0.4.4-bin.jar metapopulation.MetaPopulation", 
                 seed, tmax, nIters, burnin.time, oneLoop, extinctionLoop,
                 quasiExtinctionBound, nOffspring, reproductive.age, maxAge,
                 indiv.per.ha, reproductiveChance, vector.to.String(survivalChance),
                 vector.to.String(initial.agedistr), output.patch, output.migration,
                 output.stoch, output.population.evolution, output.population.extinction,
                 ">", tmpfile)
command
system(command)

cc <- t(as.vector(read.table(tmpfile, header = FALSE)))

Genereer wat output om resultaten te onderzoeken

source(paste0(path,"MetaPopulation/R/ibm_output.R"))
if (oneLoop) {
  EvolutionOutput(tmax = tmax, nPop = nPop, java.out = output.population.evolution, carying.cap = cc)
}
if (extinctionLoop) {
  EvolutionOutput(tmax = tmax, nPop = nPop, java.out = output.population.evolution, oneLoop = F, carying.cap = cc)
  ExtinctionOutput(nrep = nIters, tmax = tmax, java.out = output.population.extinction)
}

source(paste0(path,"MetaPopulation/R/ibm_output.R"))
  extinctionLoop = TRUE
correlation = F
  command <- paste("java -cp .:classes:jdistlib-0.4.4-bin.jar metapopulation.MetaPopulation", seed, tmax, nIters,
                   burnin.time, oneLoop, extinctionLoop, quasiExtinctionBound, nOffspring, reproductive.age, maxAge,
                   indiv.per.ha, reproductiveChance, vector.to.String(survivalChance), vector.to.String(initial.agedistr),
                   output.patch, output.migration, output.stoch, output.population.evolution,
                   output.population.extinction,">",tmpfile)
system(command)
cc<-t(as.vector(read.table(tmpfile, header = FALSE)))
ExtinctionOutput(nrep = nIters, tmax=tmax, java.out = output.population.extinction)
correlation = T 

correlation.onlyStoch = T

command <- paste("java -cp .:classes:jdistlib-0.4.4-bin.jar metapopulation.MetaPopulation", seed, tmax, nIters,
                   burnin.time, oneLoop, extinctionLoop, quasiExtinctionBound, nOffspring, reproductive.age, maxAge,
                   indiv.per.ha, reproductiveChance, vector.to.String(survivalChance), vector.to.String(initial.agedistr),
                   output.patch, output.migration, output.stoch, output.population.evolution,
                   output.population.extinction,">",tmpfile)
system(command)
ExtinctionOutput(nrep = nIters, tmax=tmax, java.out=output.population.extinction, lines=TRUE)
legend("topleft",legend = c("No correlation","Correlation"),col = 1:2, fill = 1:2)