R/multiCWRsimSPEC.R

In DeadParrot69/CWERNI: Community Wide Evolutionary Rescue Neutral Implementation

#' Multiple simulations of the Community Wide Rescue model (CWR) saving specific timepoints and also saving the output matrix for every simulation.
#'
#'@param tmax Arbitrary units of time the simulation should be run
#' @param b1 Birthrate of the mutants
#' @param d1 Deathrate of the mutants
#' @param b2 Birthrate of the residents
#' @param d2 Death rate of the residents
#' @param m12 Mutation rate of mutants to residents
#' @param m21 Mutation rate of residents to mutants
#' @param k1 Carrying capacity
#' @param abun_original A vector of initial abundances for each species such as those created by generate_spat_abund
#' @param interval After how many evenst should the population state be saved to the output matrix
#' @param nsim The number of simulations to be performed
#' @param wantedtimes Specific timepoints for which the RAC is to be plotted
#' @return this function plots the RAC of the simulations for the times specified in Wantedtimes as .tiff files in the current working directory also writes a .csv file outputting the poplation matrix for every simulation performed
#'#'@author Timo van Eldijk
#' @return this function plots all the results of the simulations in a .tiff file in the current working directory
#' @export
#'
#' @examples
#' multiCWRsimSPEC(10,0.6 , 0.1, 0.05, 0.1, 0, 0.0005, 16000,
#' generate_spat_abund(theta = 200,Ivec = rep(40,1),Jvec = c(16000)),
#' 200, 2, c(15,30,50,75,100))


multiCWRsimSPEC=function (tmax, b1, d1, b2, d2, m12, m21, k1, abun_original, interval, nsim, wantedtimes){

  #Do the simulations and write them to .csv files
  counter=1
  while(counter<=nsim){
    print (counter)
    test=CWRsimSPEC(tmax, b1, d1, b2, d2, m12, m21, k1, abun_original, interval, wantedtimes)
    utils::write.csv(test, file = paste(counter, ".csv"))
    counter=counter+1

  }


  ##### Plot racs of wantedtimes
  timepointcounter=1
  while(timepointcounter<=length(wantedtimes)){ #Cycle through the timepoints for which you want plots (whole integers!)

    wantedtime=wantedtimes[timepointcounter]
    # 2 Go throught the .csv files one by one to make RAC through time plots
    racstor=data.frame()
    counter=1
    while(counter<=nsim){
      MyData <- utils::read.csv(file=paste(counter, ".csv"), header=TRUE, sep=",")
      readdata=as.matrix(MyData)
      readdata=readdata[,2:ncol(readdata)]
      readdata=apply(readdata, 2, as.numeric)

      timepoint=readdata[1, readdata[1,]<=wantedtime]
      timepoint=max(timepoint)


      racmut=readdata [(((nrow(readdata)-1)/2)+2):nrow(readdata),readdata[1,]==timepoint]
      racres=readdata [2:(((nrow(readdata)-1)/2)+1),readdata[1,]==timepoint]

      rac=sort((racres+racmut), decreasing = T)
      racstor=rbind(racstor, rac)

      print (paste(counter, "rac", wantedtime))

      counter=counter+1

    }




    grDevices::tiff(filename=paste(wantedtime,".tiff"), width = 1000, height = 1000, units = "px", pointsize = 32)
    graphics::par(mfrow=c(1,1))
    graphics::plot(0 , xlab="Species Rank", ylab="Nr of individuals", ylim=c(1,10000), xlim=c(0,(length(abun_original))),log="y", pch = 20, cex = .8 )

    newracstor=sapply(racstor, function(x){stats::quantile(x,probs=c(0.05,0.25, 0.5, 0.75,0.95), type = 1)})
    vec1=newracstor[1,]
    vec2=newracstor[2,]
    vec3=newracstor[3,]
    vec4=newracstor[4,]
    vec5=newracstor[5,]
    graphics::lines (1:length(vec1[vec1>0]), vec1[vec1>0], cex = .6, col="grey", lwd=2, lty=2)
    graphics::lines (1:length(vec5[vec5>0]), vec5[vec5>0], cex = .6, col="grey", lwd=2, lty=2)
    graphics::lines (1:length(vec2[vec2>0]), vec2[vec2>0], cex = .6, col="grey", lwd=2, lty=3)
    graphics::lines (1:length(vec4[vec4>0]), vec4[vec4>0], cex = .6, col="grey", lwd=2, lty=3)
    graphics::lines (1:length(vec3[vec3>0]),vec3[vec3>0], pch = 20, cex = .8, col="black", lwd=2)
    graphics::polygon(c((1:length(vec2[vec2>0])), rev(1:length(vec1[vec1>0]))), c(vec2[vec2>0], rev(vec1[vec1>0])),col=grDevices::gray(0.9))
    graphics::polygon(c((1:length(vec5[vec5>0])), rev(1:length(vec4[vec4>0]))), c(vec5[vec5>0], rev(vec4[vec4>0])),col=grDevices::gray(0.9))
    graphics::polygon(c((1:length(vec3[vec3>0])), rev(1:length(vec2[vec2>0]))), c(vec3[vec3>0], rev(vec2[vec2>0])),col='lightblue')
    graphics::polygon(c((1:length(vec4[vec4>0])), rev(1:length(vec3[vec3>0]))), c(vec4[vec4>0], rev(vec3[vec3>0])),col="lightblue")
    graphics::points (1:ncol(racstor), newracstor[3,], pch = 20, cex = 0.8, col="black", lwd=2)
    graphics::points(1:ncol(racstor), sort(abun_original, decreasing=T), col="Red", pch = 20, cex = .8)


    grDevices:: dev.off()

    timepointcounter=timepointcounter+1
  }


  #Do the poptrajectories

  grDevices::tiff(filename="Trajectory.tiff", width = 1000, height = 1000, units = "px", pointsize = 32)
  graphics::plot(1, type="n", xlab="time", ylab="Total community size", ylim=c(0,20000), xlim=c(0,tmax))

  counter=1 #cycle .csv files
  while(counter<=nsim){
    MyData <- utils::read.csv(file=paste(counter, ".csv"), header=TRUE, sep=",") #reading in .csv
    readdata=as.matrix(MyData)
    readdata=readdata[,2:ncol(readdata)]
    readdata=apply(readdata, 2, as.numeric)

    traj=cbind(readdata[1,], colSums(readdata [2:nrow(readdata), ])) #Calculte total community size trajectory over time


    graphics::points (traj[,1], traj[,2], pch = 20, cex = .8) #plot it all



    print(paste(counter, "trajectory"))
    counter=counter+1

  }
  grDevices:: dev.off()


}