migflow: Calculates the maximum flow through a network

Documented in randomSTAR

#' Generate a random network
#'
#' @description This function generates a random migratory network
#'
#' @param nbreeding number of breeding sites.
#' @param nwintering number of nonbreeding residency/wintering sites.
#' @param nstop number of sites used duting migratory stopover
#' @param toplot TRUE/FALSE to determine whether the output is plotted or not
#' @param pop population size flowing through network
#' @param minforks minimum number of forks in a river branch
#' @param maxforks maximum number of forks in a river branch
#' @param anadromous TRUE or FALSE. If TRUE, then population breeds at the tips of network branches, otheriwe, it breeding "at sea"
#'
#' @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.
#'
#' @examples
#' par(mfrow=c(1,1))
#' randomSTAR(nbreeding = 6, toplot=TRUE)
#'
#'
#'
#' @import igraph
#' @importFrom stats runif
#' @export
randomSTAR <- function(pop = 100000,
                       toplot= TRUE,
                       nbreeding = 8,
                       nwintering = 4,
                       nstop=40,
                       minforks = 2,
                       maxforks = 5,
                       anadromous = TRUE){

  # # For testing purposes
  # pop = 100000
  # toplot= TRUE
  # nbreeding = 50
  # # nwintering = 3
  # # nstop=100
  # minforks = 2
  # maxforks = 5
  # anadromous = TRUE
  # # anadromous = FALSE
  # nwintering=5
  # nstop = 200
  # # nbreeding=5

  nsites = nwintering + nbreeding + nstop
  branches = ifelse(anadromous == TRUE, nwintering, nbreeding)


  # randomly determine how many branches there will be
  samples <- runif(branches,0.5,0.8)
  samples <- samples/sum(samples)
  branch_lengths <- round(nsites *samples)
  # branch_lengths <- branch_lengths[!(branch_lengths<3)]

  # !KD
  # if(any(branch_lengths < 3)){
  #   branch_lengths[which(branch_lengths < 3)] <- (branch_lengths[which(branch_lengths < 3)]
  #                                                 + ( 3- branch_lengths[which(branch_lengths < 3)]))
  # }
  #
  branch_lengths

  #ensure that branch lengths are not bigger than the number of sites
  diff <- nsites-sum(branch_lengths)
  diff
  if(diff < 0) branch_lengths[order(branch_lengths, decreasing = TRUE)[1:abs(diff)]] <-  branch_lengths[order(branch_lengths, decreasing = TRUE)[1:abs(diff)]]-1

  #    if(length(max(branch_lengths))>=abs(diff)){
  #    branch_lengths[which(branch_lengths == max(branch_lengths))[1:abs(diff)]] <-  branch_lengths[which(branch_lengths == max(branch_lengths))[1: abs(1)]]-1
  # } }

  if(diff > 0) branch_lengths[order(branch_lengths)[1:abs(diff)]] <-  branch_lengths[order(branch_lengths)[1:abs(diff)]]+1

  # if(length(max(branch_lengths))>=abs(diff)){
  # branch_lengths[which(branch_lengths == min(branch_lengths))[1:diff]] <-  branch_lengths[which(branch_lengths == min(branch_lengths))[1:diff]]+1
  # }else{

  # } }
  diff <- nsites-sum(branch_lengths)
  diff
  # if(diff < 0) branch_lengths[which(branch_lengths == max(branch_lengths))] <-  branch_lengths[which(branch_lengths == max(branch_lengths))]+diff
  # if(diff > 0) branch_lengths[which(branch_lengths == min(branch_lengths))] <-  branch_lengths[which(branch_lengths == min(branch_lengths))]+diff


  branch_lengths
  site_counter = 1
  store = list()
  for (b in 1:branches){ # for every branch

    # Create a fake list of sites where animals were seen at, with latitude, longitude and number of anumals seen there
    site_list <- data.frame(Lat= runif(branch_lengths[b], min=-20, max=40),
                            Lon= runif(branch_lengths[b], min=-10, max=10),
                            Pop=runif(branch_lengths[b], min=500, max=10000))

    #sort according to latitude
    site_list = site_list[order(site_list$Lat, decreasing=T),]
    site_list$Site = site_counter:(site_counter + branch_lengths[b] - 1)

    if(anadromous == TRUE){
      site_list$NB <- 0
      site_list$NB[1]<- 1
    } else{
      site_list$B <- 0
      site_list$B[1]<- 1
    }

    # add a dummy breeding and wintering site
    site_list<- rbind(
      c(41,0,pop,0,0),
      site_list,
      c(-21,0,pop,9999,0))

    # create a distance matrix based on these data
    dist <- point2DIST(site_list)
    rownames(dist)[1] <- "supersource"
    colnames(dist)[1] <- "supersource"
    rownames(dist)[branch_lengths[b]+2] <- "supersink"
    colnames(dist)[branch_lengths[b]+2] <- "supersink"

    # Make the network directed
    network <- directedNET(max(dist[2:(branch_lengths[b]+1), 2:(branch_lengths[b]+1)]) - dist[2:(branch_lengths[b]+1), 2:(branch_lengths[b]+1)], include_diagonal = TRUE)
    network[network==0] <- NA

    if(branch_lengths[b] == 2){
      network[1,2] = 1
    }
    if(branch_lengths[b] == 3){
      network[1,2] = 1
      network[2,3] = 1
    }
    if(branch_lengths[b] > 3){

      i=1
      forks <- round(runif(1,1,2))
      sorted <- network[i,]
      sorted[which(!is.na(sorted))] <- sort(network[i,which(!is.na(sorted))],index.return = TRUE)$ix
      keep <- which (sorted >= (max(sorted,na.rm=TRUE) - (forks-1)))
      # val= runif(length(keep),0,1)
      network[i,] <- NA
      network[i,keep] <- 1 # val/(sum(val))

      i=2
      forks <- round(runif(1,minforks,maxforks))
      sorted <- network[i,]
      sorted[which(!is.na(sorted))] <- sort(network[i,which(!is.na(sorted))],index.return = TRUE)$ix
      keep <- which (sorted >= (max(sorted,na.rm=TRUE) - (forks-1)))
      already_inflowing <- which(!is.na(network[i-1,]))
      keep <- keep[which(!(keep %in% already_inflowing))]
      # val= runif(length(keep),0,1)
      network[i,] <- NA
      network[i,keep] <- 1# val/(sum(val))
      if (sum(network[,i], na.rm=TRUE) == 0){
        network[i-1,i]<- 1
      }


      suppressWarnings(for (i in 3:nrow(network)){
        forks <- round(runif(1,minforks,maxforks))
        sorted <- network[i,]
        sorted[which(!is.na(sorted))] <- sort(network[i,which(!is.na(sorted))],index.return = TRUE)$ix

        if(any(duplicated(sorted[!is.na(sorted)]))){
          idx <- which(any(duplicated(sorted[!is.na(sorted)])) == TRUE)
          sorted[!is.na(sorted)][idx:length(sorted)] <- sorted[!is.na(sorted)][idx:length(sorted)]+1
        }

        keep <- which (sorted >= (max(sorted,na.rm=TRUE) - (forks-1)))
        already_inflowing <- which(apply(network[1:i-1,],2,sum,na.rm=TRUE)>0)
        keep <- keep[which(!(keep %in% already_inflowing))]
        # val= runif(length(keep),0,1)
        network[i,] <- NA
        network[i,keep] <- 1# val/(sum(val))

        if (sum(network[,i], na.rm=TRUE) == 0){
          network[i-1,i]<- 1
        }
      })
    }
    site_counter <- site_counter + branch_lengths[b]
    store[[b]] <- list(network, site_list)
  }

  # Put all the tree branches together into a big network
  network <- matrix(0, ncol=nsites, nrow=nsites)
  for (b in 1:branches){
    network[store[[b]][[2]]$Site[2:(length(store[[b]][[2]]$Site)-1)],
            store[[b]][[2]]$Site[2:(length(store[[b]][[2]]$Site)-1)]] <- store[[b]][[1]]
  }
  colnames(network) = rownames(network) = 1:nsites


  # Put all the sites together into one big site list
  site_list <- unique(do.call(rbind,lapply(1:branches, function(x) store[[x]][[2]])))
  site_list <- site_list[order(site_list$Site),]
  site_list$Site[site_list$Site == 0]<- "supersource"
  site_list$Site[site_list$Site == 9999]<- "supersink"

  # specify the sinks based on how close sites are the sinks (i.e. the longest branches)
  dist <- point2DIST(site_list)
  rownames(dist)[1] <- "supersource"
  colnames(dist)[1] <- "supersource"
  rownames(dist)[branch_lengths[b]+2] <- "supersink"
  colnames(dist)[branch_lengths[b]+2] <- "supersink"
  sinks = which(apply(network,1, function(x) sum(which(x>0)))==0)
  if (anadromous==TRUE){
    # if(nbreeding != "ALL")
    sinks = sinks[which(sort(dist[sinks,"supersink"],index.return = TRUE)$ix <= nbreeding)]
    site_list$B <- 0
    site_list$B[sinks+1] <- 1
    network <- addSUPERNODE(network, sources= site_list$Site[site_list$NB == 1],sinks = sinks)
  }else{
    # if(nwintering != "ALL")
    sinks = sinks[which(sort(dist[sinks,"supersink"],index.return = TRUE)$ix <= nwintering)]
    site_list$NB <- 0
    site_list$NB[sinks+1] <- 1
    network <- addSUPERNODE(network, sources= site_list$Site[site_list$B == 1],sinks = sinks)
  }

  # if(nbreeding != "ALL") sinks = sinks[which(sort(dist[sinks,"supersink"],index.return = TRUE)$ix <= nbreeding)]
  #
  # #Add supersource and sink nodes
  # network <- addSUPERNODE(network, sources= site_list$Site[2],
  #                         sinks = sinks)

  network <- ifelse(network == Inf, pop, network)
  network[,"supersink"] <- ifelse(network[,"supersink"] == pop,
                                  max(dist[,"supersink"]) - dist[,"supersink"],
                                  network[,"supersink"])
  network[,"supersink"] <- network[,"supersink"] /sum(network[,"supersink"] )#*pop


  network["supersource",] <- ifelse(network["supersource",] == pop,
                                    max(dist["supersource",]) - dist["supersource",],
                                    network["supersource",])
  network["supersource",] <- network["supersource",] /sum(network["supersource",])

  network[is.na(network)] = 0
  network <- network * pop

  #-------------------------

  if(anadromous == TRUE){
    to_connect <- site_list$Site[site_list$NB==1]
  }else{
    to_connect <- site_list$Site[site_list$B==1]
  }

  network[to_connect,to_connect] <- rep(network[1,to_connect],each=length(to_connect))

  #-------------------------

  weight <- graph_from_adjacency_matrix(network,  mode="directed", weighted = TRUE)

  flow = max_flow(weight, source = V(weight)["supersource"],
                  target = V(weight)["supersink"], capacity = E(weight)$weight )

  nodes = get.edgelist(weight, names=TRUE)
  nodes = as.data.frame(nodes)
  nodes$flow = flow$flow
  nodes$Lat_from = unlist(lapply(1:nrow(nodes), function(i) as.numeric(site_list$Lat[site_list$Site %in% nodes[i,1]])))
  nodes$Lon_from = unlist(lapply(1:nrow(nodes), function(i) as.numeric(site_list$Lon[site_list$Site %in% nodes[i,1]])))
  nodes$Lat_to   = unlist(lapply(1:nrow(nodes), function(i) as.numeric(site_list$Lat[site_list$Site %in% nodes[i,2]])))
  nodes$Lon_to   = unlist(lapply(1:nrow(nodes), function(i) as.numeric(site_list$Lon[site_list$Site %in% nodes[i,2]])))


  # neti <- weight
  # E(neti)$weight <- flow$flow
  # network <- as.matrix(as_adjacency_matrix(neti, attr="weight"))
  # neti[neti== "."] <- flow$flow

  # if (toplot == TRUE){
  weight <- delete.vertices(weight, c(1, length(V(weight))))
  index = which(nodes$V1 != "supersource" & nodes$V2 != "supersink")
  nodeindex = which(nodes$V1 != "supersource")
  # nodes[nodeindex,]
  sizes <- unlist(lapply(1:nsites ,function(x) sum(nodes$flow[nodeindex][nodes$V1[nodeindex]==x])))


  site_list$NM <- site_list$SM <- 0
  site_list$NM[site_list$B == 0 & site_list$NB == 0] <- 1
  site_list$SM[site_list$B == 0 & site_list$NB == 0] <- 1

  # Bcols = c("black",ifelse(network[3:(nsites+1),"supersink"] > 0, "orange", "royalblue4"))
  lyt = layout_with_kk(weight)#*2

  # if (toplot == TRUE){
  #     plot(weight, layout= lyt, edge.width = ((flow$flow[index]/pop)*20), edge.arrow.mode=0,
  #          edge.color = "royalblue4",
  #          vertex.color =  Bcols,
  #          vertex.label="",
  #          vertex.size = ((sizes/pop)*20)+5)
  #   }

  site_list$Pop <- apply(network,1,sum)#c(pop,sizes, pop)
  site_list$Pop[length(site_list$Pop)] <- pop
  site_list$Lon[2:(nrow(site_list)-1)] <- lyt[,1]
  site_list$Lat[2:(nrow(site_list)-1)] <- lyt[,2]

  if (anadromous == TRUE){
    colnames(network)[1]<- rownames(network)[1] <- "supersink"
    colnames(network)[ncol(network)]<- rownames(network)[nrow(network)] <- "supersource"
    mirror_network <- network # matrix(rev(network),ncol=52)
    mirror_network <- mirror_network[, ncol(mirror_network):1]
    mirror_network <- mirror_network[nrow(mirror_network):1, ]
    mirror_network <- t(mirror_network)
    SMnet <- mirror_network
    NMnet <- network
  } else {
    mirror_network <- network # matrix(rev(network),ncol=52)
    mirror_network <- mirror_network[, ncol(mirror_network):1]
    mirror_network <- mirror_network[nrow(mirror_network):1, ]
    mirror_network <- t(mirror_network)
    SMnet <- network
    NMnet <- mirror_network
  }


  #---------------------------------------------
  # Join the two networks
  #---------------------------------------------

  topright=matrix(0,nrow(SMnet),ncol(NMnet)+1)
  colnames(topright) <- c("NB",colnames(NMnet))

  bottomleft=matrix(0, nrow(NMnet)+1, ncol(SMnet))
  rownames(bottomleft) <- c("NB",rownames(NMnet))

  bottomright = cbind(rep_len(0,nrow(NMnet)),NMnet)
  bottomright = rbind(rep_len(0,ncol(bottomright)),bottomright)
  # bottomright[1,1] = 1

  top = cbind(SMnet,topright)
  bottom = cbind(bottomleft,bottomright)

  network=rbind(top,bottom)

  colnames(network)=c(paste0("S",colnames(SMnet)),"NB",paste0("N",colnames(NMnet)))
  rownames(network)=c(paste0("S",rownames(SMnet)),"NB",paste0("N",rownames(NMnet)))

  network["Ssupersink","NB"]<- pop
  network["NB","Nsupersink"]<- pop


  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)["Ssupersource"],
                  target = V(weight)["Nsupersource"], capacity = E(weight)$weight)

  sites <-site_list

  # plot flow network
  nodes = get.edgelist(weight, names=TRUE)
  nodes = as.data.frame(nodes)
  nodes$flow = flow$flow
  nodes$V1 <- substring(nodes$V1, 2)
  nodes$V2 <- substring(nodes$V2, 2)

  nodes = nodes[nodes$V1 != "supersource" & nodes$V1 != "supersink" & nodes$V2 != "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]])))


  # library(shape)
  # par(mfrow=c(1,1))
  # par(mar=c(4,4,4,4))
  index=2:(nrow(sites)-1)
  if (toplot == TRUE){
    plot(sites$Lon[index], sites$Lat[index], pch=16,
         cex=0, xlab="", ylab="", xaxt="n", yaxt = "n",
         frame.plot=FALSE)


    index=1:nrow(nodes)
    segments(x0 = nodes$Lon_from[index],
             y0 = nodes$Lat_from[index],
             x1 = nodes$Lon_to[index],
             y1 = nodes$Lat_to[index],
             col= "black",
             lwd=(nodes$flow[index]/(max(nodes$flow)))*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
  colorz = ifelse(sites$B[index]==1,"royalblue",ifelse(sites$NB[index]==1,"orange","gray"))
  if (toplot == TRUE){
    points(sites$Lon[index],
           sites$Lat[index],
           pch=21,
           cex=(((nodeflowplot$x)/
                   as.numeric(max(nodeflowplot$x)))+0.4)*4,
           bg=colorz , col="black")
  }



  return(list( network = network,
               # tracks = tracks,
               sites = site_list))

}