R/randomNET.R

In KiranLDA/migflow: Calculates the maximum flow through a network

#' Generate a random network
#'
#' @description This function generates a random migratory network
#'
#'
#' @param nsites number of sites to be generated in the network
#' @param nedges number of edges in the network. If all sites should be connected then use nedges="ALL"
#' @param pop population size flowing through network
#' @param mean_dist average distance the species can travel (this can be estimated from real data)
#' @param sd_dist standard deviation of the distances a species can travel (this can be estimated from real data)
#' @param Latrange geographic range of species, by default latitude restricted to c(-20,40),
#' @param Lonrange geographic range of species, by default longitude restricted to c(-10,10),
#' @param Poprange min and max population sizes to be randomly generated, by default c(100,10000)
#' @param toplot TRUE/FALSE to determine whether the output is plotted or not
#' @param nbreeding Number of breeding sites
#' @param nwintering Number of wintering sites
#'
#' @return a list containting the network which was randomly generated,
#' the tracks that were randomly generated, and the sites that were randomly generated for animals to use.
#'
#' @import igraph
#' @importFrom stats rnorm runif
#' @export
randomNET <- function(nsites=10,
                      nedges = "ALL",
                      nbreeding = 3,
                      nwintering = 4,
                      pop = 100000,
                      mean_dist = 2000,
                      sd_dist = 1000,
                      Latrange = c(-20,40),
                      Lonrange = c(-10,10),
                      Poprange = c(100,10000),
                      toplot = TRUE){

  # generate random tracks
  tracks <- rnorm(100, mean_dist, sd_dist)

  #Create a fake list of sites where animals were seen at, with latitude, longitude and number of anumals seen there
  sites <- data.frame(Lat= runif(nsites, min=Latrange[1], max =Latrange[2]),
                      Lon= runif(nsites, min=Lonrange[1], max=Lonrange[2]),
                      Pop = round(runif(nsites, min=Poprange[1], max=Poprange[2])))

  #sort according to latitude
  sites = sites[order(sites$Lat, decreasing=T),]
  sites$Site= 1:nsites

  # create a distance matrix based on these data
  dist <- point2DIST(sites)

  # calculate the probability of going between these sites given the distance the animal can travel
  Dist_P <- distPROB(tracks, dist, adjust=1, plot=F)

  # Calculate prioritisation of population using a site
  Pop_P <- nodePopPROP(sites, population = pop)

  #Calculate the azimuth angle
  Azi_P <- absAZIMUTH(dist, lonlats=sites )

  # make birds/animals prefer sites which a larger prioritisation of the population has been seen and where the distance is better
  network <- Dist_P * Pop_P * Azi_P

  # Make the network directed
  network <- directedNET(network, include_diagonal = TRUE)

  #estimate number of birds entering and exiting sites based on distance, population count and azimuth
  # network <- popPROP(network, population = pop)

  if(nedges != "ALL"){
    # only keep nodes generated by random graph generator
    graph <- as.matrix(as_adj(sample_gnm(n=nsites, m=nedges, directed=F)))
    network <- graph*network
    colnames(network)<-1:nsites
    rownames(network)<-1:nsites
  }

  #Add supersource and sink nodes
  network <- addSUPERNODE(network, sources= sites$Site[1:nbreeding], sinks = sites$Site[(nrow(sites)-nwintering+1):nrow(sites)])
  # network[1,(1:nbreeding)+1] = sites$Pop[1:nbreeding]/sum(sites$Pop[1:nbreeding])*pop
  # network[((nrow(sites)-nwintering+1):nrow(sites))+1,ncol(network)] = sites$Pop[(nrow(sites)-nwintering+1):nrow(sites)]/sum(sites$Pop[(nrow(sites)-nwintering+1):nrow(sites)])*pop
  #
  network[1,(1:nbreeding)+1] = sites$Pop[1:nbreeding]/sum(sites$Pop[1:nbreeding])
  network[((nrow(sites)-nwintering+1):nrow(sites))+1,ncol(network)] = sites$Pop[(nrow(sites)-nwintering+1):nrow(sites)]/sum(sites$Pop[(nrow(sites)-nwintering+1):nrow(sites)])

  network <- popPROP(network, population = pop)



  network[is.na(network)] = 0

  sites<- rbind(
    c(41,0,pop,"supersource"),
    sites,
    c(-21,0,pop,"supersink"))

  #created a weigted igraph network
  if (toplot == T){
    weight <- graph_from_adjacency_matrix(network,  mode="directed", weighted = TRUE)

    # run the population through the network a forst time
    flow = max_flow(weight, source = V(weight)["supersource"],
                    target = V(weight)["supersink"], capacity = E(weight)$weight )

    # plot flow network
    par(mfrow=c(1,1))
    par(mar=c(0,0,0,0))
    index=2:(nrow(sites)-1)
    plot(sites$Lon[index], sites$Lat[index], pch=16,
         cex=0)

    nodes = get.edgelist(weight, names=TRUE)
    nodes = as.data.frame(nodes)
    nodes$flow = flow$flow

    nodes = nodes[nodes$V1 != "supersource" & nodes$V2 != "supersink" ,]

    nodes$Lat_from = unlist(lapply(1:nrow(nodes), function(i) as.numeric(sites$Lat[sites$Site %in% nodes[i,1]])))
    nodes$Lon_from = unlist(lapply(1:nrow(nodes), function(i) as.numeric(sites$Lon[sites$Site %in% nodes[i,1]])))
    nodes$Lat_to   = unlist(lapply(1:nrow(nodes), function(i) as.numeric(sites$Lat[sites$Site %in% nodes[i,2]])))
    nodes$Lon_to   = unlist(lapply(1:nrow(nodes), function(i) as.numeric(sites$Lon[sites$Site %in% nodes[i,2]])))

    index=2:(nrow(nodes)-1)
    segments(x0 = nodes$Lon_from[index],
             y0 = nodes$Lat_from[index],
             x1 = nodes$Lon_to[index],
             y1 = nodes$Lat_to[index],
             lwd=(nodes$flow[index]/pop)*30)


    # sort sites by flow
    nodeflow = merge(aggregate(nodes$flow, by=list(Category=as.character(nodes$V1)), FUN=sum),
                     aggregate(nodes$flow, by=list(Category=as.character(nodes$V2)), FUN=sum), all=T)
    nodeflow$x = as.numeric(nodeflow$x)
    nodeflow = data.frame( unique(as.matrix(nodeflow[ , 1:2 ]) ))
    nodeflow$x = as.numeric(as.character(nodeflow$x))
    nodeflow = nodeflow[nodeflow$Category != "supersource" & nodeflow$Category != "supersink",]

    # make sure it is numeric
    nodeflow$Category = as.numeric(as.character(nodeflow$Category))

    # plot sites
    nodeflowplot = nodeflow[order(nodeflow$Category),]
    index=as.numeric(nodeflowplot$Category)+1
    points(sites$Lon[index],
           sites$Lat[index],
           pch=21,
           cex=(((nodeflowplot$x)/
                   as.numeric(max(nodeflowplot$x)))+0.4)*4,
           bg="orange", col="black")
  }
  return(list( network = network,
               tracks = tracks,
               sites = sites  ))

}