INApreliminary: Impact Network Analysis

Documented in ntsteps2

#' Changes management info status and bioentity establishment status at each node across multiple time steps
#'
#' This function gives the info status and establishment status at each node after multiple time steps. Its use follows use of the function \code{setup2}, or another source of node locations and initial conditions.  (Note that communication and dispersal status are given for the beginning of the time step (rather than for the end), so the corresponding output for these has length equal to the number of time steps plus 1.  In contrast, adoption and establishment output is of length equal to the number of time steps.) It uses the following functions, as needed: genmovnet, spreadstep, makedec, and estab. It draws on output from function setup2 in terms of whether communication occurs (estinfo), generated geographic locations of nodes as needed (genloc), and the initial locations for management information and the bioentity.

#'
#' Updated 2020-09-05

#' @param nsteps number of time steps to be evaluated

#' @param infon is T if estimated effect size exceeds threshold for communication (so communication can occur), F otherwise (can be evaluated in estinfo and/or setup2)

#' @param readseam3 if T, the communication adjacency matrix is read in rather than being generated by function genmovnet
#' @param readbpam3 if T, the dispersal adjacency matrix is read in rather than being generated by function genmovnet
#' @param seam3 communication adjacency matrix, read in if readseam3=T (rows as sources and columns as sinks)
#' @param bpam3 dispersal adjacency matrix, read in if readbpam3=T (rows as sources and columns as sinks)
#' @param geocoords3n matrix of x,y coordinates of nodes

#' @param seamdist3 com (genmovnet) the function of distance used to estimate movement probability - 'random' (not related to distance) or 'powerlaw' 
#' @param seampla3 com (genmovnet) inverse power law parameter a in ad^(-b)
#' @param seamplb3 com (genmovnet) inverse power law parameter b in ad^(-b)
#' @param seamrandp3 com (genmovnet) random case, probably of link

#' @param vect1cn com (spreadstep) status of nodes before spread

#' @param bpamdist3 disp (genmovnet) the function of distance used to estimate movement probability - 'random' (not related to distance) or 'powerlaw' 
#' @param bpampla3 disp (genmovnet) inverse power law parameter a in ad^(-b)
#' @param bpamplb3 disp (genmovnet) inverse power law parameter b in ad^(-b)
#' @param bpamrandp3 disp (genmovnet) random case, probably of link

#' @param vect1dn disp (spreadstep) status of nodes before spread

#' @param probadoptmean3 (makedec) mean probability of adopting management if informed
#' @param probadoptsd3 (makedec) sd in truncated normal distribution of probability of adoption
#' @param comvec vector of 1=info is present, 0=info is not present
#' @param probadoptvec3 vector of probabilities of adoption if informed for nodes in the network (\code{readprobadoptvec3} determines whether \code{probadoptvec3} is read in to \code{makedec} or generated within \code{makedec} based on \code{probadoptmean3} and \code{probadoptsd3})
#' @param readprobadoptvec3 if T, then a VECTOR of probadoptvec3 values is read in; if F, \code{probadoptmean3} and \code{probadoptsd3} are used to generate the vector

#' @param maneffdir3 if maneffdir3='decrease_estab', the management reduces the probability of establishment; if maneffdir3='increase_estab', the management reduces the probability that establishment does NOT occur

#' @param probestabmean3 (estab) mean probability of establishment (new or CONTINUED) in absence of management
#' @param probestabsd3 (estab) sd in probability of establishment in truncated normal distribution
#' @param maneffmean3n (estab) mean management effect (proportional reduction in estabp) 
#' @param maneffsd3n (estab) sd of management effect
#' @param probestabvec3 (estab) vector of probabilities of establishment (read in or generated when the function is run, depending on \code{readprobestabvec3} equal to T or F)
#' @param readprobestabvec3 (estab) if T, then a vector \code{probestabvec3} is read in; otherwise the vector is generated using \code{probestabmean3} and \code{probestabsd3}

#' @param plotmpn if T, plots of resulting presence of information and species are generated

#' @keywords simulation experiments
#' @export 
#' @import truncnorm

#' @examples

#' x2 <- ntsteps2(nsteps=3, infon=T, geocoords3n=matrix(c(1,1, 1,2, 1,3, 2,1, 2,2, 2,3),byrow=T,ncol=2), vect1cn=c(1,1,1,0,0,0), vect1dn=c(0,0,0,1,1,1), readseam3=F, seamdist3='powerlaw', seampla3=1, seamplb3=1, readbpam3=F, bpamdist3='powerlaw', bpampla3=1, bpamplb3=1, probadoptmean3=0.5, probadoptsd3=0.1, probestabmean3=0.5, probestabsd3=0.1, maneffdir3='decrease_estab', maneffmean3n=0.5, maneffsd3n=0.1,  readprobestabvec3=F, readprobadoptvec3=F, plotmpn=T)

#' x2z <- ntsteps2(nsteps=3, infon=T, geocoords3n=matrix(c(1,1, 1,2, 1,3, 2,1, 2,2, 2,3),byrow=T,ncol=2), vect1cn=c(1,1,1,0,0,0), vect1dn=c(0,0,0,1,1,1), readseam3=F, seamdist3='powerlaw', seampla3=0.5, seamplb3=2, readbpam3=F, bpamdist3='powerlaw', bpampla3=0.5, bpamplb3=2, probadoptmean3=0.5, probadoptsd3=0.1, probestabmean3=0.5, probestabsd3=0.1, maneffdir3='decrease_estab', maneffmean3n=0.5, maneffsd3n=0.1, readprobestabvec3=F, readprobadoptvec3=F, plotmpn=T)

#' x3 <- ntsteps2(nsteps=2, infon=T, geocoords3n=matrix(runif(n=100)*100,byrow=T,ncol=2), vect1cn=c(rep(1,10),rep(0,40)), vect1dn=sample(c(rep(1,10),rep(0,40))), readseam3=F, seamdist3='powerlaw', seampla3=1, seamplb3=1, readbpam3=F, bpamdist3='powerlaw', bpampla3=1, bpamplb3=1, probadoptmean3=0.5, probadoptsd3=0.1, probestabmean3=0.5, probestabsd3=0.1, maneffdir3='decrease_estab', maneffmean3n=0.5, maneffsd3n=0.1, readprobestabvec3=F, readprobadoptvec3=F, plotmpn=T)

#' x9 <- ntsteps2(nsteps=3, infon=T, geocoords3n=matrix(c(1,1, 1,2, 1,3, 2,1, 2,2, 2,3),byrow=T,ncol=2), vect1cn=c(1,1,1,0,0,0), vect1dn=c(0,0,0,1,1,1), readseam3=T, seam3=matrix(c(0,1,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1, 1,0,0,0,0,0),byrow=T,ncol=6), readbpam3=T, bpam3=matrix(c(0,0,0,0,0,0, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1, 1,0,0,0,0,0),byrow=T,ncol=6), probadoptmean3=0.3, probadoptsd3=0.1, probestabmean3=0.2, probestabsd3=0.1, maneffdir3='decrease_estab', maneffmean3n=0.5, maneffsd3n=0.1, readprobestabvec3=F, readprobadoptvec3=F, plotmpn=T)

#' x9p <- ntsteps2(nsteps=3, infon=T, geocoords3n=matrix(c(1,1, 1,2, 1,3, 2,1, 2,2, 2,3),byrow=T,ncol=2), vect1cn=c(1,1,1,0,0,0), vect1dn=c(0,0,0,1,1,1), readseam3=T, seam3=matrix(c(0.5,1,0.2,0.6,0.1,0.3, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1, 1,0,0,0,0,0),byrow=T,ncol=6), readbpam3=T, bpam3=matrix(c(0.2,0.2,0.2,0.2,0.2,0.2, 0,0,1,0,0,0, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1, 1,0,0,0,0,0),byrow=T,ncol=6), probadoptmean3=0.3, probadoptsd3=0.1, probestabmean3=0.2, probestabsd3=0.1, maneffdir3='decrease_estab', maneffmean3n=0.5, maneffsd3n=0.1, readprobestabvec3=F, readprobadoptvec3=F, plotmpn=T)

#' x12 <- ntsteps2(nsteps=3, infon=F, geocoords3n=matrix(c(1,1, 1,2, 1,3, 2,1, 2,2, 2,3),byrow=T,ncol=2), vect1cn=c(1,1,1,0,0,0), vect1dn=c(0,0,0,1,1,1), readseam3=F, seamdist3='powerlaw', seampla3=1, seamplb3=1, readbpam3=F, bpamdist='powerlaw', bpampla3=1, bpamplb3=1, probadoptmean3=0.3, probadoptsd3=0.1, probestabmean3=0.2, probestabsd3=0.1, maneffdir3='decrease_estab', maneffmean3n=0.5, maneffsd3n=0.1, readprobestabvec3=F, readprobadoptvec3=F, plotmpn=T)

#' x13 <- ntsteps2(nsteps=5, infon=T, geocoords3n=matrix(c(1,1, 1,2, 1,3, 2,1, 2,2, 2,3),byrow=T,ncol=2), vect1cn=c(1,1,1,0,0,0), vect1dn=c(0,0,0,1,1,1), readseam3=F, seamdist3='powerlaw', seampla3=1, seamplb3=1, readbpam3=F, bpamdist3='powerlaw', bpampla3=1, bpamplb3=1, probadoptmean3=0.3, probadoptsd3=0.1, probestabmean3=0.2, probestabsd3=0.1, maneffdir3='decrease_estab', maneffmean3n=0.5, maneffsd3n=0.1, readprobestabvec3=F, readprobadoptvec3=F, plotmpn=T)



ntsteps2 <- function(nsteps, geocoords3n, infon, vect1cn, vect1dn, readseam3, seam3, seamdist3, seampla3, seamplb3, seamrandp3, readbpam3, bpam3, bpamdist3, bpampla3, bpamplb3, bpamrandp3, readprobadoptvec3, probadoptvec3, probadoptmean3, probadoptsd3, maneffdir3, maneffmean3n, maneffsd3n, readprobestabvec3, probestabvec3, probestabmean3, probestabsd3, plotmpn=F){

# Numbering of steps refers to the numbering used in the INA User's Manual (1-5)

nodlen <- nrow(geocoords3n) # number of nodes

# 5. If the observed management effect exceeds the threshold for communication, or if the 'science of science' threshold is not used - the following commands are implemented
# (If these criteria are not met, there is no communication)

if(infon){ 

  # 3.1. communication adjacency matrix - underlying probabilities

  # if generating the communication adj mat (not reading it in)

  if (!readseam3){ 
    seam3 <- genmovnet(geocoords4n=geocoords3n, amdist4=seamdist3, amrandp4=seamrandp3, ampla4=seampla3, amplb4=seamplb3)
  }

  # 4.1. decision making - underlying probabilities of adoption
  
  # if generating the vector of probabilities of adoption for each node (not reading it in)

################ check whether this is actually done in function makedec

  if (!readprobadoptvec3){
    probadoptvec3 <- truncnorm::rtruncnorm(n=nodlen, a=0, b=1, mean = probadoptmean3, sd=probadoptsd3)
  }

} else if (!infon){ # 5. if there is no information spread

  if (!readseam3){seam3 <- 0}

  probadoptvec3 <- 0
}

### end first infon-determined portion for socioeconomic network


# 3.2. bioentity dispersal adjacency matrix - underlying probabilities

# if generating the bioentity dispersal network (rather than reading it in)

if (!readbpam3){ 
  bpam3 <- genmovnet(geocoords4n=geocoords3n, amdist4=bpamdist3, amrandp4=bpamrandp3, ampla4=bpampla3, amplb4=bpamplb3)  
}

# 4.2. establishment - underlying probabilities of establishment
  
  # if generating the vector of probabilities of establishment for each node (not reading it in)

################ but does estab do the same thing as the following?

  if (!readprobestabvec3){

    probestabvec3 <- truncnorm::rtruncnorm(n=nodlen, a=0, b=1, mean = probestabmean3, sd=probestabsd3)
  }


# prepare output lists for timesteps

vect1cL <- as.list(1:(nsteps+1))
# communication status vector at beginning of time step
vect1cL[[1]] <- vect1cn 
vect1dL <- as.list(1:(nsteps+1))
# dispersal status at beginning of time step
vect1dL[[1]] <- vect1dn

estabvec2 <- vect1dn

decvecL <- as.list(1:nsteps)
estabvecL <- as.list(1:nsteps)

# if there is no communication, 
#   the decision is always not to adopt
if(!infon){
  vect2c <- 0*(1:nodlen)
  decvec2 <- vect2c
  probadoptvec3 <- vect2c
}

### each time step

for(j in 1:nsteps) {

if(infon){ # if estimated management effect exceeds threshold for communication, or no threshold for communication is used

  # communication outcome

  # generate an adjacency matrix for this time step, tseam3

  # note that if seam3 is type pa (presence/absence) then the 
  #   following command has no effect
  # if seam3 is type pr (probability) then the structure can
  #   change in each time step


  tseam3 <- seam3 > matrix(runif(n=nodlen*nodlen),ncol=nodlen)

# the adjacency matrix tseam3 is composed of 1s and 0s

# the diagonal is set to 1 so that information is not lost 
#    in this step
  
  diag(tseam3) <- 1


  vect2c <- spreadstep(amat=tseam3, vect1=vect1cL[[j]]) # communication

  vect1cL[[j+1]] <- vect2c

  # outcome for adoption of management

  decvec2 <- makedec(comvec=vect2c, geocoords4n=geocoords3n, probadoptvec4=probadoptvec3, readprobadoptvec4=T, plotmp=F)

  decvecL[[j]] <- decvec2

} ### end of second infon-determined actions

# dispersal outcome 

# which links exist this time step in biophysical 
#   adjacency matrix

tbpam3 <- bpam3 > matrix(runif(n=nodlen*nodlen),ncol=nodlen)

# set the diagonal equal to zero so that the species will not 
#   be lost in these steps - whether it is lost is 
#   determined in the establishment steps

diag(tbpam3) <- 1

# note that bioentity dispersal and communication are handled differently - only locations where the bioentity has established can be a source of spread

# the adjacency matrix tbpam3 is composed of 1s and 0s

vect2d <- spreadstep(amat=tbpam3, vect1=estabvec2) # dispersal

vect1dL[[j+1]] <- vect2d

### outcome for establishment

######################### why is readprobestabvec4 set to hard T while at the same time... maybe it needs a different name for estab?? to avoid confusion ....

estabvec2 <- estab(decvec=decvec2, dispvec=vect2d, probestabvec4=probestabvec3, readprobestabvec4=T, maneffdir4=maneffdir3, maneffmean4n=maneffmean3n, maneffsd4n=maneffsd3n, geocoords4n=geocoords3, plotmp=F)

estabvecL[[j]] <- estabvec2

# plots, if called for

if(plotmpn){
    plot(geocoords3n, main=paste('Time',j,', 1 Information, shaded = yes'))
    if(!is.null(dim(geocoords3n[vect2c==1,]))){
      points(geocoords3n[vect2c==1,], pch=16)
    } 
    else if(is.null(dim(geocoords3n[vect2c==1,]))){
      points(geocoords3n[vect2c==1,][1], geocoords3n[vect2c==1,][2], pch=16)
    }

   plot(geocoords3n, main=paste('Time',j,', 2 Adoption, shaded = yes'))
    if(!is.null(dim(geocoords3n[decvec2==1,]))){
      points(geocoords3n[decvec2==1,], pch=16)
    } 
    else if(is.null(dim(geocoords3n[decvec2==1,]))){
      points(geocoords3n[decvec2==1,][1], geocoords3n[decvec2==1,][2], pch=16)
    }

    plot(geocoords3n, main=paste('Time',j,', 3 Sp dispersal, shaded = yes'))
    if(!is.null(dim(geocoords3n[vect2d==1,]))){
      points(geocoords3n[vect2d==1,], pch=16)
    } 
    else if(is.null(dim(geocoords3n[vect2d==1,]))){
      points(geocoords3n[vect2d==1,][1], geocoords3n[vect2d==1,][2], pch=16)
    }

    plot(geocoords3n, main=paste('Time',j,', 4 Sp establishment, shaded = yes'))
    if(!is.null(dim(geocoords3n[estabvec2==1,]))){
      points(geocoords3n[estabvec2==1,], pch=16)
    } 
    else if(is.null(dim(geocoords3n[estabvec2==1,]))){
      points(geocoords3n[estabvec2==1,][1], geocoords3n[estabvec2==1,][2], pch=16)
    }
}

} #### end of time steps

list(seam3=seam3, probadoptvec3=probadoptvec3, bpam3=bpam3, probestabvec3=probestabvec3, vect1cL=vect1cL, vect1dL=vect1dL, decvecL=decvecL, estabvecL=estabvecL)
}