R/reedgraph-scratch.R

library(doParallel)
registerDoParallel(cores=(detectCores()-1))

###    Copyright (C) 2012 Martin J Reed              martin@reednet.org.uk
###    University of Essex, Colchester, UK
###
###    This program is free software; you can redistribute it and/or modify
###    it under the terms of the GNU General Public License as published by
###    the Free Software Foundation; either version 2 of the License, or
###    (at your option) any later version.
###
###    This program is distributed in the hope that it will be useful,
###    but WITHOUT ANY WARRANTY; without even the implied warranty of
###    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
###    GNU General Public License for more details.
###
###    You should have received a copy of the GNU General Public License along
###    with this program; if not, write to the Free Software Foundation, Inc.,
###    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

rg.generate.documentation <- function() {
    package.skeleton.dx(pkgdir="/Users/mjreed/main/R/reedgraph",
                        excludePattern=paste0("(reedgraph-scratch.R)|",
                            "(reedgraph-altbloom.R)|",
                            "(reedgraph-bloom.R)|",
                            "(ddlpn.R)"
                            ))
}

### NOTE THIS IS USING INLINEDOCS
### EXAMPLE
rg.inlinedocs.example <- function # Example function
### Here is some maths in the longer description \eqn{n} 
### which probably has more lines
##references<< Abelson, Hal; Jerry Sussman, and Julie
##Sussman. Structure and Interpretation of Computer
##Programs. Cambridge: MIT Press, 1984.
(n,    ##<< description of vairiables 
 times
### here is a different way: Number of Fermat tests to perform. More
### tests are more likely to give accurate results.
 ){
  if(times==0)TRUE
  ##seealso<< \code{\link{fermat.test}}
  else if(fermat.test(n)) is.pseudoprime(n,times-1)
  else FALSE
### the return value logical TRUE if n is probably prime.
}


### count number of packets to get one PPM packet from each router
rg.ppm.count.required <- function(N,p) {

  all <- FALSE
  count <- 0
  received_marks <- rep(FALSE,N)
  while(!all) {
    count <- count + 1
    mark <- NULL
    for(i in 1:N) {
      if(runif(1) < p) {
        mark <- i
      }
    }
    if(!is.null(mark)) {
      received_marks[mark] <- TRUE
    }
    if(sum(received_marks) == N) {
      break
    }
  }

  return(count)
}

rg.ppm.count.required.averaged <- function(N,p,averages) {

  sum <- 0
  for(i in 1:averages) {
    sum <- rg.ppm.count.required(N,p) + sum
  }
  return(sum/averages)
}

rg.ppm.savage <- function(N,p) {
  e <- (log(N)) / ( p*(1-p)^(N-1) )
  return(e)
}

rg.ppm.reed <- function(N,p) {
  sum <- 1
  for(i in (N-1):0) {
    print(i)
    sum <- sum  + 1/ ( p * (1-p)^i )
  }
  
  return(sum)
}

rg.ppm.test <- function(N,p) {

  prob <- p * (1-p) * (1-p)^2 * (1-p * (1-p))

  return(prob)
}
### Generate graph to simulate attack graph
### Warning this only produces a tree at the moment
rg.gen.attack.graph <- function(depth=4,mindeg=2,maxdeg=5,lower=0.3) {

  nodes <- list()
  alledges <- list()
  nodes[[1]] <- c("1")
  allnodes <- c(nodes[[1]])
  ## for each level
  for(i in 2:depth) {
    cat("creating nodes level ",i,"\n")
    count <- 1
    nodes[[i]] <- c()


    ## for each node at higher level
    cat("nodes[i-1]length=",length(nodes[[i-1]]),"\n")
    first <- TRUE
    for(j in 1:length(nodes[[i-1]])) {
      edges <- c()


      num <- floor(runif(1,mindeg,maxdeg+1))
      edges <- c(paste(i,".",as.character(count),sep=""))
      count <- count +1
      for(k in 2:num) {
        edges <- append(edges,paste(i,".",as.character(count),sep=""))
        count <- count +1
      }
      if(first) {
        nodes[[i]] <- edges
        first <- FALSE
      } else {
        nodes[[i]] <- append(nodes[[i]],edges)
      }

      ### this does not do anything yet. needs moving up to where each
      ### new node name is created and add an edge
      if(runif(1) < lower ) {
        prevd <- floor(runif(1,1,i+1))
        prev <- floor(runif(1,1,length(nodes[[prevd]])+1))
        print("i would connect to")
        print(nodes[[prevd]][prev])
      }

      alledges[[nodes[[i-1]][j]]]<- edges
      
    }
    print(nodes[[i]])
    allnodes <- append(allnodes,nodes[[i]])

  }
  allnodes <- append(allnodes,"t")
  for(j in 1:length(nodes[[depth]])) {
    alledges[[nodes[[depth]][j]]] <- "t"
  }

  
  print(alledges)
  g <- new("graphNEL",nodes=allnodes,edgeL=alledges,"directed")

  ## should not need to do this but some bug seems to be evident if we do not
  g <- addEdge("t","1",g,1)
  g <- removeEdge("t","1",g)

  return(g)
}

### THIS MAY HAVE BUGS!
### Not a full implementation! This uses the Fleischer maximum commodity
### flow algorithm to compute an approximate max-flow. THIS IS NOT THE
### BEST WAY TO DO THIS! but it is all I had "off the shelf"
### Demands must be less than
### a valid flow. Once the max-flow is computed it can determine the best
### s-t cut (the dual problem).
### g - graphNEL with capacities set
### s the starting node (character)
### t the finishing node (characeter)
### progress - printed progress bar if TRUE
### e - default=0.05, bound on approximate max-flow
### tolerance = 0.05 anything with less than this portion of the capacity
###                  in the flow will be assumed to be a congesting flow
###                  and will be used to form the cut
rg.st.cut <- function(g,s,t,progress=TRUE,e=0.05,tolerance=0.05,res=NULL) {

  ### need to do this better as this will be very slow
  validflow <- min(rg.edgeVector(g,"capacity"))
  demand <- rg.set.demand(validflow,c(s,t))
  if(is.null(res)) {
    res <- rg.fleischer.max.concurrent.flow(g,progress=progress,demands=demand,
                                            e=e)
  }
  flow <- as.double(edgeData(res$gflow,attr="weight"))
  cap <- as.double(edgeData(res$gflow,attr="capacity"))
  #cat("(cap-flow)/cap=",(cap-flow)/cap,"\n")
                   
  cuts <- NULL
  gcut <- g
  for(e in edgeNames(g)) {
    from <- strsplit(e,"~")[[1]][1]
    to <- strsplit(e,"~")[[1]][2]
    cap <- as.double(edgeData(g,from=from,to=to,attr="capacity"))
    flow <- as.double(edgeData(res$gflow,from=from,to=to,attr="weight"))
    if( (cap-flow)/cap <= tolerance ) {
      cuts <- c(cuts,from,to)
      gcut <- removeEdge(from,to,gcut)
    }
  }
  numcomp <- length(connComp(gcut))
  if(numcomp !=2 ) {
    cat("Error, Number of components=",numcomp," it should be 2\n")
    cat("Returning res, consider passing this back in with lower tolerance\n")
    return(res)
  } else {
    if("1" %in% connComp(gcut)[[1]] ) {
      if ("t" %in% connComp(gcut)[[2]] ) {

      } else {
        cat("Something is wrong expecting 1 and t to be in ",
            "in different components but they are not\n")
      }
    } else if ("t" %in% connComp(gcut)[[1]] ) {
      if ("1" %in% connComp(gcut)[[2]] ) {
        
      } else {
        cat("Something is wrong expecting 1 and t to be in ",
            "in different components but they are not\n")
      }

    } else {
      cat("Something is wrong expecting 1 and t to be in ",
          "in different components but they are not\n")
    }
  }
  
  return(cuts)
}


rg.analyse.tests <- function(res) {

  var <- data.frame()
  for(N in unique(res$N)) {
    var <- rbind(var,lapply(res,subset,res$N==N & res$L==res$N * res$N))
  }
  return(var)
}

rg.test.diameter <- function() {
    Size <- seq(100,800,100)
    Diameter <- c()
    for(N in Size) {
        g <- rg.generate.random.graph(N,2,25)
        Diameter <- c(Diameter,diameter(igraph.from.graphNEL(g)))
    }
    res <- data.frame(Size=Size,Diameter=Diameter)
    return(res)
}

rg.test.integer.plot <- function(res) {
  ## convert to long format
  resl <- melt(res,measure.vars=c("MF","FF"),variable.name="Algorithm")
  ## summarise
  ress <- summarySE(resl,"value",groupvars=c("Algorithm","Load"))
  ## plot!
  pd <- position_dodge(0.5) # offest so errorbars do not obscure
  ggplot(ress, aes(x=Load,y=value, colour=Algorithm)) + geom_errorbar(aes(ymin=value-ci, ymax=value+ci), width=1, position=pd) + geom_line(position=pd) + geom_point(position=pd, size=3, shape=21, fill="white") + ylab("Blocking Probability") + expand_limits(y=0)

  ggplot(res, aes(x=factor(Load))) + geom_boxplot(aes(y=FF),position=pd,width=0.25) + geom_boxplot(aes(y=MF),position=pd,width=0.25)
  
}

### NOT COMPLETE YET, the augmentation works but  it crashes in 
rg.test.integer <- function(emin,emax,estep,repeats,N=10,wavelengths=4,e=0.1) {
  g <- rg.generate.random.graph(N,2,25)
  M <- length(edgeMatrix(g))/2
  g <- rg.set.capacity(g,1.0)

  au <- rg.augment.graph.for.wavelengths(g,wavelengths)
  ga <- au$ga
  linkgroupmap <- au$linkgroupmap
  
  ga <- rg.set.weight(ga,0.0)
  sp.blocking <- c()
  flow.blocking <- c()
  flowint.blocking <- c()
  ff.blocking <- c()
  ff.gamma <- c()
  int.gamma <- c()
  int.mgamma <- c()
  ff.mgamma <- c()
  dcount <- c()
  int.results <- list()
  for(i in seq(emin,emax,estep) ) {
    for(j in 1:repeats) {
      demands <- rg.gen.demands(g,i)
      ## adapt the demands to the new source/sink nodes.
      for(j in names(demands) ){
        demands[[j]]$source <- paste0(demands[[j]]$source,"s")
        demands[[j]]$sink <- paste0(demands[[j]]$sink,"t")
      }
      numdems <- length(demands)
      dcount <- c(dcount, numdems)
      
      
      ff.results <- rg.sp.ff.max.concurrent.flow(ga,demands)
      accepted <- rg.count.accepted.demands(
        rg.check.blocked.demands(ff.results$demands,ga))
      ff.blocking <- c(ff.blocking,(numdems-accepted)/numdems)
      ff.gamma <- c(ff.gamma,ff.results$gamma)
      weights <- as.double(edgeData(ff.results$gflow,attr="weight"))
      capacities <- as.double(edgeData(ff.results$gflow,attr="capacity"))
      gamma <- (capacities - weights) / capacities
      mgamma <- mean(gamma)
      ff.mgamma <- c(ff.mgamma,mgamma)
      
      int.results <- rg.max.concurrent.flow.int(ga,demands,e=e)
      ch.demands <- rg.check.blocked.demands(int.results$demands,ga)
      accepted <- rg.count.accepted.demands(ch.demands)
      gtmp <- rg.max.concurrent.flow.graph(ga,ch.demands)
      weights <- as.double(edgeData(gtmp,attr="weight"))
      capacities <- as.double(edgeData(gtmp,attr="capacity"))
      gamma <- (capacities - weights) / capacities
      mgamma <- mean(gamma)
      int.mgamma <- c(int.mgamma, mgamma)
      flowint.blocking <- c(flowint.blocking,(numdems-accepted)/numdems)
      int.gamma <- c(int.gamma, min(gamma))
      
      results <- data.frame("Load"=dcount,
                            "FFB" = ff.blocking,
                            "IntB"=flowint.blocking)
      
      print(results[nrow(results),],digits=4)
    }
  }
  results <- data.frame("Load"=dcount,
                        "FF" = ff.blocking,
                        "MF"=flowint.blocking)
  return(results)
}




# Summarizes data.
## Gives count, mean, standard deviation, standard error of the mean,
##   and confidence interval (default 95%).  data: a data frame.
##   measurevar: the name of a column that contains the variable to be
##   summariezed groupvars: a vector containing names of columns that
##   contain grouping variables na.rm: a boolean that indicates
##   whether to ignore NA's conf.interval: the percent range of the
##   confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
    require(plyr)

    # New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }

    # This is does the summary; it's not easy to understand...
    datac <- ddply(data, groupvars, .drop=.drop,
                   .fun= function(xx, col, na.rm) {
                           c( N    = length2(xx[,col], na.rm=na.rm),
                              mean = mean   (xx[,col], na.rm=na.rm),
                              sd   = sd     (xx[,col], na.rm=na.rm)
                              )
                          },
                    measurevar,
                    na.rm
             )

    # Rename the "mean" column    
    datac <- rename(datac, c("mean"=measurevar))

    datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

    # Confidence interval multiplier for standard error
    # Calculate t-statistic for confidence interval: 
    # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
    ciMult <- qt(conf.interval/2 + .5, datac$N-1)
    datac$ci <- datac$se * ciMult

    return(datac)
}

rg.check.blocked.demands <- function # Checks demands in order to lower flow if blocked
### The results from a max-flow (min-congestion) algorithm may assign
### more traffic than is actually possible. This function will test
### each demand (and path in demand) in turn and only allocate the flow that is possible
### thus it is possible to see which demands are blocked by comparing flows against demands
### Note that trying the demands in order is not the optimum. It should be randomised
### rather trying it in order. Another possibility is ordering the demands/paths in
### order of length, shortest first. This way more demands are likely to be accepted.
(demands, ##<< demand list in standard format with paths
  g       ##<< graphNEL object with capacity attribute on edges
 ) {
  
  g <- rg.set.all.graph.edge.weights(g)
  for(nd in names(demands)) {
    d <- demands[[nd]]
    ##reassign the flow to the actual allocated
    flow <- 0
    for(p in names(d$paths)) {
      #cat(nd,":",p,"=",d$paths[[p]],"\n")
      pv <- as.vector(strsplit(p,"|",fixed=TRUE)[[1]])
      fromlist <- pv[1:length(pv)-1]
      tolist <- pv[2:length(pv)]
      
      weights <- as.double(edgeData(g,from=fromlist,
                                    to=tolist,att="weight"))
      capacities <- as.double(edgeData(g,from=fromlist,
                                       to=tolist,att="capacity"))

      minfree <- min(capacities-weights)

      free <- (capacities-weights)
      ## only for debugging
      ##minind <- which(free == minfree)
      ##for(m in minind) {
      ##  cat("Most constrained ",fromlist[m],"->",tolist[m],
      ##      weights[m],capacities[m]," ")
      ##}
      ## end only for debugging
      
      if (minfree < d$paths[[p]]) {
        ##cat("blocked")
        demands[[nd]]$paths[[p]] <- minfree
        flow <- flow + minfree
      } else {
        flow <- flow + d$paths[[p]]
      }
      ##cat("\n")
      weights <- weights + demands[[nd]]$paths[[p]]
      edgeData(g,from=fromlist,
               to=tolist,
               att="weight") <- weights
      
    }
    demands[[nd]]$flow <- flow
    
  }
  demands
  ### Demands with only the accepted traffic flow up to the limit of the graph capacities
  ### other demands will have their flow set to the maximum allowable.
}


rg.path.lengths <- function(demands) {
  lengths <- c()
  for(i in demands) {
    for(j in names(i$paths)) {
      lengths <- c(lengths,length(strsplit(j,"\\|")[[1]])-1)
    }
  }
  return(lengths)
}

## This was the routine used for the icc2012 paper
rg.run.tests <- function(min=8,max=20,averages=10,e=0.05,target=0.1,dirname=NULL) {
  vecN <- c()
  vecL <- c()
  vecSPgamma <- c()
  vecFLOWgamma <- c()
  vecINTgamma <- c()
  vecFLOWtime <- c()
  vecMAXmpls <- c()
  vecMAXbf <- c()
  run <- 1
  
  for(N in min:max ) {
    for(a in 1:averages) {
      cat("\nN=",N, "\n")

      scenario <- rg.test.int.versus.nonint.flow(N=N,L=NULL,target=target,e=e)
      if( ! is.null(dirname) ) {
        save(scenario,file=paste(dirname,"/",run,".robj",sep=""))
        cat("writing",paste(dirname,"/",run,".robj",sep=""),"\n")
      }
      vecN <- c(vecN,N)
      vecSPgamma <- c(vecSPgamma,scenario$sp.results$gamma)
      vecFLOWgamma <- c(vecFLOWgamma,scenario$flow.results$gamma)
      vecINTgamma <- c(vecINTgamma,scenario$int.results$bestgamma)
      vecFLOWtime <- c(vecFLOWtime,scenario$flow.results$runtime)
      vecMAXmpls <- c(vecMAXmpls,scenario$flow.results$max.mpls)
      vecMAXbf <- c(vecMAXbf,scenario$flow.results$max.bf)
      run <- run +1
    }
    
  }
  cat("\n")
  results <- data.frame(N=vecN,
                        SPgamma=vecSPgamma,
                        INTgamma=vecINTgamma,
                        FLOWgamma=vecFLOWgamma,
                        FLOWtime=vecFLOWtime,
                        MAXmpls=vecMAXmpls,
                        MAXbf=vecMAXbf
                        )
  
  return(results)
}

rg.run.tests.target <- function(min=0,max=1,step=0.1,averages=10,e=0.05,N=20,dirname=NULL) {
  vecTarget <- c()
  vecL <- c()
  vecSPgamma <- c()
  vecFLOWgamma <- c()
  vecINTgamma <- c()
  vecFLOWtime <- c()
  vecMAXmpls <- c()
  vecMAXbf <- c()
  run <- 1
  
  for(target in seq(min,max,step)) {
    for(a in 1:averages) {
      cat("\ntarget=",target, "\n")

      scenario <- rg.test.int.versus.nonint.flow(N=N,L=NULL,target=target,e=e)
      if( ! is.null(dirname) ) {
        save(scenario,file=paste(dirname,"/",run,".robj",sep=""))
        cat("writing",paste(dirname,"/",run,".robj",sep=""),"\n")
      }
      vecTarget <- c(vecTarget,target)
      vecSPgamma <- c(vecSPgamma,scenario$sp.results$gamma)
      vecFLOWgamma <- c(vecFLOWgamma,scenario$flow.results$gamma)
      vecINTgamma <- c(vecINTgamma,scenario$int.results$bestgamma)
      vecFLOWtime <- c(vecFLOWtime,scenario$flow.results$runtime)
      vecMAXmpls <- c(vecMAXmpls,scenario$flow.results$max.mpls)
      vecMAXbf <- c(vecMAXbf,scenario$flow.results$max.bf)
      run <- run +1
    }
    
  }
  cat("\n")
  results <- data.frame(target=vecTarget,
                        SPgamma=vecSPgamma,
                        INTgamma=vecINTgamma,
                        FLOWgamma=vecFLOWgamma,
                        FLOWtime=vecFLOWtime,
                        MAXmpls=vecMAXmpls,
                        MAXbf=vecMAXbf
                        )
  
  return(results)
}


rg.run.tests.process <- function(dirname) {
  vecN <- c()
  vecL <- c()
  vecSPgamma <- c()
  vecSPmaxflow <- c()
  vecFLOWgamma <- c()
  vecFLOWmaxflow <- c()
  vecINTgamma <- c()
  vecINTmaxflow <- c()
  vecFLOWtime <- c()
  vecMAXmpls <- c()
  vecMAXbf <- c()
  vecRatio <- c()
  run <- 1
  file <- paste(dirname,"/",run,".robj",sep="")
  while(file.exists(file)) {
    load(file)
    N <- length(nodes(scenario$g))
    vecN <- c(vecN,N)
    vecSPgamma <- c(vecSPgamma,scenario$sp.results$gamma)
    vecSPmaxflow <- c(vecSPmaxflow,
                      rg.calc.max.flow(scenario$sp.results$demands,scenario$sp.results$gamma))
    vecFLOWgamma <- c(vecFLOWgamma,scenario$flow.results$gamma)
    vecFLOWmaxflow <- c(vecFLOWmaxflow,
                        rg.calc.max.flow(scenario$flow.results$demands,scenario$flow.results$gamma))
#    vecFLOWmaxflow <- c(vecFLOWmaxflow,
#                        rg.calc.max.flow(scenario$flow.results$demands,0))
    vecINTgamma <- c(vecINTgamma,scenario$int.results$bestgamma)
    vecINTmaxflow <- c(vecINTmaxflow,
                       rg.calc.max.flow(scenario$int.results$intdemands,
                                        scenario$int.results$bestgamma))
#    vecINTmaxflow <- c(vecINTmaxflow,rg.calc.max.flow(scenario$int.results$intdemands,0))
    vecFLOWtime <- c(vecFLOWtime,scenario$flow.results$runtime)
    vecMAXmpls <- c(vecMAXmpls,scenario$flow.results$max.mpls)
    vecMAXbf <- c(vecMAXbf,scenario$flow.results$max.bf)
    vecRatio <- c(vecRatio,scenario$flow.results$ratio)
    run <- run+1
    file <- paste(dirname,"/",run,".robj",sep="")
  }
  print(length(vecSPmaxflow))
  print(length(vecFLOWmaxflow))
  print(length(vecINTmaxflow))
  results <- data.frame(N=vecN,
                        SPgamma=vecSPgamma,
                        SPmaxflow=vecSPmaxflow,
                        INTgamma=vecINTgamma,
                        INTmaxflow=vecINTmaxflow,
                        FLOWgamma=vecFLOWgamma,
                        FLOWmaxflow=vecFLOWmaxflow,
                        FLOWtime=vecFLOWtime,
                        MAXmpls=vecMAXmpls,
                        MAXbf=vecMAXbf,
                        Ratio=vecRatio
                        )
}

rg.calc.max.flow <- function(demands,gamma) {
  flow <- 0.0
  for(i in demands) {
    flow <- flow + i$flow / (1 - gamma)
  }
  return(flow)
}

rg.run.tests.process.maxflow <- function(dirname) {
  vecN <- c()
  vecL <- c()
  vecSPgamma <- c()
  vecFLOWgamma <- c()
  vecINTgamma <- c()
  vecFLOWtime <- c()
  vecMAXmpls <- c()
  vecMAXbf <- c()
  vecRatio <- c()
  run <- 1
  file <- paste(dirname,"/",run,".robj",sep="")
  while(file.exists(file)) {
    load(file)
    N <- length(nodes(scenario$g))
    vecN <- c(vecN,N)
    vecSPgamma <- c(vecSPgamma,scenario$sp.results$gamma)
    vecFLOWgamma <- c(vecFLOWgamma,scenario$flow.results$gamma)
    vecINTgamma <- c(vecINTgamma,scenario$int.results$bestgamma)
    vecFLOWtime <- c(vecFLOWtime,scenario$flow.results$runtime)
    vecMAXmpls <- c(vecMAXmpls,scenario$flow.results$max.mpls)
    vecMAXbf <- c(vecMAXbf,scenario$flow.results$max.bf)
    vecRatio <- c(vecRatio,scenario$flow.results$ratio)
    run <- run+1
    file <- paste(dirname,"/",run,".robj",sep="")
  }
  results <- data.frame(N=vecN,
                        SPgamma=vecSPgamma,
                        INTgamma=vecINTgamma,
                        FLOWgamma=vecFLOWgamma,
                        FLOWtime=vecFLOWtime,
                        MAXmpls=vecMAXmpls,
                        MAXbf=vecMAXbf,
                        Ratio=vecRatio
                        )
}

rg.demands.to.csv <- function(demands,file,paths=TRUE) {
  cat("",file=file)
  if(paths) {
    for(i in demands) {
      for(j in names(i$paths)) {
        cat(i$paths[[j]],",",gsub("\\|",",",j),"\n",sep="",file=file,append=TRUE)
      }
    }
  } else {
    for(i in demands) {
        cat(i$demand,",",i$source,",",i$sink,"\n",sep="",file=file,append=TRUE)
    }
  }
}

rg.smart.round <- function(x) {
  y <- floor(x)
  indices <- tail(order(x-y), round(sum(x)) - sum(y))
  y[indices] <- y[indices] + 1
  y
}

rg.random.round <- function(demand,paths) {
  total <- 0
  newpaths <- paths
  for(path in names(paths)){
    newpaths[[path]] <- floor(paths[[path]])
    total <- total + newpaths[[path]] 
  }
  #return(NULL)
  while(total < demand) {
    # we want the order to be random, in case we do this again
    # otherwise earlier paths will be biased to get rounded before we finish
    for(path in sample(names(paths))){
      if(rbinom(1,1,paths[[path]] %% 1)) {
        newpaths[[path]] <- newpaths[[path]] + 1
        total <- total +1
        if(total >= demand)
          break
      }
    }
    
  }
  paths <- list()
  for(path in names(newpaths)){
    if(newpaths[[path]] > 0) {
      paths[[path]] <- newpaths[[path]]
    }
  }
 
  return(paths)
}



rg.karcius.run.all <- function(gnames=NULL,targets=NULL,zoo=zoo,dir="./",incomm=NULL,progress=FALSE,e=0.1) {
  g <- NULL
  for(gname in gnames) {
    if(gname == "lattice") {
      g <- NULL
      gi <- make_lattice(c(2,3))
      g <- igraph.to.graphNEL(as.directed(gi))
    } else {
      gi <- zoo[[gname]]
      g <- igraph.to.graphNEL(as.directed(simplify(gi)))
    }
    N <- length(nodes(g))
    if(is.null(incomm)) {
      comm <- rg.gen.demands(g,val=10*rlnorm(N*N,16,1)/exp(16.5))
    } else {
      comm <- incomm
    }
    foreach(target=targets) %dopar% {
      #for(target in targets){
      network.name <- paste0(dir,gname,"-",as.character(target))
      cat(network.name,"\n")
      results <- rg.karcius(g=g,comm=comm,network.name=network.name,target=target,e=e,progress=progress)
    }
  }
}


## returns the bounds on gamma as absolute gamma values in the form
## probs can be a vector if you want to calculate more than one value
# c(maxbound,firstprob, secondprob .....)
rg.karcius.determine.bounds <- function(flow.results,e,probs=c(0.05,0.01,0.0)) {
  gcount <- flow.results$gflow
  gcount <- rg.set.weight(gcount,0.0)

  for(d in flow.results$demands) {
    for(p in names(d$paths)) {
      pv <- as.vector(strsplit(p,"|",fixed=TRUE)[[1]])
      fromlist <- pv[1:length(pv)-1]
      tolist <- pv[2:length(pv)]
      
      weights <- as.double(edgeData(gcount,from=fromlist,
                                    to=tolist,att="weight"))
      weights <- weights + 1
      edgeData(gcount,from=fromlist,
               to=tolist,
               att="weight") <- weights

    }
  }
  # tricky this, the bound we are using is for opt being lambda' of the transformed problem
  # and a found value lambda (which is likely to be smaller than the optimum)
  # lambda' >= lambda >= lambda * (1+e)^3
  # but we have transformation 1/lambda=gamma, so we have for found gamma that the opt gamma'
  # gamma' <= 1 - (1-gamma)/1-e)^-3
  maxgammabound=1 - (1-flow.results$gamma)/(1-e)^-3
  mingammabound <- rep(Inf,length(probs))
  for(pi in 1:length(probs)) {
    for(ep in rg.edgeL(gcount)) {
      weight <- as.numeric(edgeData(flow.results$gflow,from=ep[1],to=ep[2],att="weight"))
      capacity <- as.numeric(edgeData(flow.results$gflow,from=ep[1],to=ep[2],att="capacity"))
      n <- as.numeric(edgeData(gcount,from=ep[1],to=ep[2],att="weight"))
      # the most we can round is by n units (where n is the number of flows)
      if(p == 0.0){
        gammaest <- (capacity - (weight + n)) / capacity
      }
      # this is an estimate of a bound, it might be pessimistic
      add <- n/2 *sqrt(-6 * log(probs[pi]) / n)
      # it certainly cannot be bigger than n
      if(add > n) add <- n
      gammaest <- (capacity - (weight + add))/capacity
      if (mingammabound[pi] > gammaest) {
        mingammabound[pi] <- gammaest
      }
    }
    
  }
  return(c(maxgammabound,mingammabound))
}

rg.karcius <- function(g=NULL,comm=NULL,dir="./",network.name="unknown-net",target=0.3,e=0.1,progress=FALSE) {
  N <- length(nodes(g))
  # note this was 10000 before Sept 2018
  g <- rg.set.capacity(g,500)
  #comm <- rg.gen.demands(g,val=10)
  g <- rg.set.weight(g,0.0)
  results <- rg.minimum.congestion.flow(g,comm,e=e,progress=progress)
  ## now work out how much to scale demand to meet a max-flow target
  multiplier <- (1-target)/(1-results$gamma)
  ## and now scale the results
  comm <- rg.rescale.demands(comm,multiplier,integer=TRUE)
  #cat("done rescale\n")
  flow.results <- rg.minimum.congestion.flow(g,comm,e=e,progress=progress)
  #cat("done rg.minimum.congestion.flow\n")
  int.results <- rg.max.concurrent.flow.int(g,flow.results$demands,e=e,progress=progress)
  #cat("done rg.max.concurrent.flow.int\n")

  sp.results <- rg.sp.max.concurrent.flow(g,comm)
  #cat("done rg.sp.max.concurrent.flow\n")
  flow.results.rounded <- flow.results
  flow.results.random <- flow.results
  #cat("rounding",length(flow.results$demands),"\n")
  for(i in 1:length(flow.results.rounded$demands)) {
    #cat("num paths",length(flow.results$demands[[i]]$paths),"\n")
    #print(as.numeric(flow.results$demands[[i]]$paths))
    random.paths <- rg.random.round(flow.results$demands[[i]]$demand,flow.results$demands[[i]]$paths)
    flow.results.random$demands[[i]]$paths <- random.paths

    ## this rounding needs putting inside the rg.smart.round function, leave it here for now
    rounded.flows <- rg.smart.round(as.numeric(flow.results.rounded$demands[[i]]$paths))
    if(sum(rounded.flows) > flow.results.rounded$demands[[i]]$demand ) {
      rounded.flows[which.max(rounded.flows)] <- rounded.flows[which.max(rounded.flows)] -
        sum(rounded.flows) + flow.results.rounded$demands[[i]]$demand
    } else if(sum(rounded.flows) < flow.results.rounded$demands[[i]]$demand ) {
      rounded.flows[which.max(rounded.flows)] <- rounded.flows[which.max(rounded.flows)] +
        sum(rounded.flows) - flow.results.rounded$demands[[i]]$demand
    }
    if(sum(rounded.flows) > flow.results.rounded$demands[[i]]$demand ) {
      cat(i,"too high","\n")
    } else if(sum(rounded.flows) > flow.results.rounded$demands[[i]]$demand ) {
      cat(i,"too low","\n")
    }
    count <- 1
    for(j in names(flow.results.rounded$demands[[i]]$paths)) {
      flow.results.rounded$demands[[i]]$paths[[j]] <- rounded.flows[count]
      count <- count+1
    }
  }
  #cat("finished rounding\n")

  flow.results.rounded$gflow <- rg.max.concurrent.flow.graph(flow.results.rounded$gflow,flow.results.rounded$demands)
  flow.results.random$gflow <- rg.max.concurrent.flow.graph(flow.results.random$gflow,
                                                            flow.results.random$demands)
  #cat("done flow results\n")
  flow.results.rounded$gamma <- rg.mcf.find.gamma(flow.results.rounded$gflow)
  flow.results.random$gamma <- rg.mcf.find.gamma(flow.results.random$gflow)
  #cat("building output\n")
  output <- list()
  output$g <- g
  output$sp.results <- sp.results
  output$flow.results <- flow.results
  output$flow.results.rounded <- flow.results.rounded
  output$flow.results.random <- flow.results.random
  output$int.results <- int.results
  output$comm <- comm
  network.name <- paste0(dir,network.name,"-")
  bounds <- rg.karcius.determine.bounds(flow.results,e=e,probs=c(0.05,0.01,0))
  #cat(network.name,"sp=",sp.results$gamma,", flow=",flow.results$gamma,", int=",int.results$gamma,"\n")
  cat(network.name," sp=",sp.results$gamma,", int=",int.results$gamma,"\n",
      ", flowrounded=",flow.results.rounded$gamma,
      ", flowrandom=",flow.results.random$gamma,
      ", flow=",flow.results$gamma, ", flowub=",bounds[1], ", flowlb0.05=",bounds[2],
      ", flowlb0.01=",bounds[3],", flowlb=",bounds[4],"\n")
  cat("sp=",sp.results$gamma,", int=",int.results$gamma,"\n",
      ", flowrounded=",flow.results.rounded$gamma,
      ", flowrandom=",flow.results.random$gamma,
      ", flow=",flow.results$gamma, ", flowub=",bounds[1], ", flowlb0.05=",bounds[2],
      ", flowlb0.01=",bounds[3],", flowlb=",bounds[4],
      file=paste0(network.name,"results.txt"))
  file <- paste0(network.name,"integer-flows.csv")
  rg.demands.to.csv(int.results$demands,file)
  file <- paste0(network.name,"sp-flows.csv")
  rg.demands.to.csv(sp.results$demands,file)
  file <- paste0(network.name,"split-flows.csv")
  rg.demands.to.csv(flow.results$demands,file)
  file <- paste0(network.name,"rounded-flows.csv")
  rg.demands.to.csv(flow.results.rounded$demands,file)
  file <- paste0(network.name,"random-flows.csv")
  rg.demands.to.csv(flow.results.random$demands,file)
  file <- paste0(network.name,"graph.csv")
  write.table(as_adjacency_matrix(igraph.from.graphNEL(g),sparse=FALSE,attr="capacity"),file=file)
  file <- paste0(network.name,"demands.csv")
  rg.demands.to.csv(comm,file=file,paths=FALSE)
  #cat("written output\n")
  return(output)
}


## test used for PURSUIT TM theory and ICC paper
## N - the number of nodes in the graph
## L - number of demands to generate (NULL will do N(N-1) demands)
## target - the demands is automatically scaled to give at least this much
##          free capacity on the most congested edge for the fractional
##          multicommodity flow (try approx 0.3 for first attempt)
## e - optimisation parameter
## g - a graph to try with, if NULL generate a graph
rg.test.int.versus.nonint.flow <- function(N=8,L=NULL,target=0.0,e=0.1,g=NULL) {
  if(is.null(g)) {
    g <- rg.generate.random.graph(N,2,25)
    
  } else {
    N <- length(nodes(g))
  }
  # number of edges
  M <- length(edgeMatrix(g))/2
  if(is.null(L)) {
    ## generate using lognormal traffic as in
    ## Nucci, Antonio and Sridharan, Ashwin and Taft, Nina,
    ## The problem of synthetically generating IP traffic matrices: initial recommendations,
    ## SIGCOMM Comput. Commun. Rev., vol 35(3) July 2005 pp19-32
    ## this generates a mean of 1 with a variance that matches the paper
    comm <- rg.gen.demands(g,val=rlnorm(N*N,16,1)/exp(16.5))
  } else {
    comm <- rg.gen.demands(g,L,val=rlnorm(L,16,1)/exp(16.5))
  }

  ## set some random capacities (this should be an input parameter)
  g <- rg.set.capacity(g,runif(M,5,10))
  g <- rg.set.weight(g,0.0)
  ##results <- rg.sp.max.concurrent.flow(g,comm)
  results <- rg.minimum.congestion.flow(g,comm,e=e,progress=TRUE)
  ## now work out how much to scale demand to meet a max-flow target
  multiplier <- (1-target)/(1-results$gamma)
  ## and now scale the results
  comm <- rg.rescale.demands(comm,multiplier)
  ## calculate the flow when using strict shortest path routing
  sp.results <- rg.sp.max.concurrent.flow(g,comm)
  ## calculate the flow when using minimum congestion flow (fractional)
  ## but this time with the scaled demand
  flow.results <- rg.minimum.congestion.flow(g,comm,e=e,progress=TRUE)
  ## count how many flows we have, 
  gcount <- rg.count.edge.flows(g,flow.results$demands)
  flow.results$max.mpls <- max(as.double(edgeData(gcount,att="weight")))
  flow.results$max.bf <- max(as.double(lapply(graph::edges(g),length)))
  ## calculate the integer results
  int.results <- rg.max.concurrent.flow.int.c(g,flow.results$demands,
                                              e=e,progress=TRUE)

  output <- list()
  output$g <- g
  output$sp.results <- sp.results
  output$flow.results <- flow.results
  output$int.results <- int.results
  output$comm <- comm
  cat("sp=",sp.results$gamma,",flow=",flow.results$gamma,"int=",int.results$bestgamma,"\n")
  return(output)
}

rg.gen.all.pairs.demands <- function(g,val=1.0) {
  nodes <- nodes(g)
  comm <- list()
  n=1;

  k <- 1
  vallen <- length(val)
  for(i in nodes) {
    for(j in nodes) {
      if ( i != j ) {
        comm[[as.character(n)]] <- list(source=i,sink=j,demand=val[k])
        k <- k + 1
        if(k > vallen) k <- 1
        n <- n+1
      }
    }
  }
  return(comm)
    
}


## Generate a random graph and augmented graph representing the
## wavelength switched paths through the network.
## n - number of nodes
## mindeg - minimum node degree (default = 3)
## maxdeg - maximum node degree (default = 7)
## grid - wavelength numbers to be used on each span (ITU grid 1-73), special
##        value of zero means the augmented graph is the same as the original
##        graph
## lamdacap - the bandwidth of the wavelengths in Gb/s (default 1.0)
##
## output - graph$g the original graph
##          graph$lg the augmented lambda graph
##          graph$accessnodes the list of access nodes (all of the original graph)

rg.generate.random.lambda.graph <- function(g,n=10,mindeg=3,maxdeg=7,grid=c(0),lambdacap=1.0) {

  #g <- rg.generate.random.graph(n,mindeg,maxdeg)

  nodelist <- c()

  accessnodes <- c()
  for (n in nodes(g)) {
    accessnodes <- c(accessnodes,paste("A",n,sep="."))
  }
  nodelist <- accessnodes
  
  for (n in nodes(g)) {
    for(l in grid) {
      nodelist <- c(nodelist,paste("A",n,l,sep="."))
    }
  }

  
  for (e in rg.edgeL(g)) {
    for(l in grid) {
      s <- e[1]
      t <- e[2]
      nodelist <- c(nodelist,paste(s,t,l,sep="."))
    }
  }
  print(nodelist)

  ag <- new("graphNEL",nodes=nodelist,edgemode="directed")


  for(l in grid) {
    for(n in nodes(g)) {
      ag <- addEdge(paste("A",n,sep="."),
                    paste("A",n,l,sep="."),
                    ag,
                    c(0))
      for(u in graph::edges(g)[[n]] ) {
        if( u == n) next
        
        ag <- addEdge(paste("A",n,l,sep="."),
                      paste(n,u,l,sep="."),
                      ag,
                      c(0))
      }
    }
  }
  for(l in grid) {
    for (e in rg.edgeL(g)) {
      s <- e[1]
      t <- e[2]

      ag <- addEdge(paste(s,t,l,sep="."),
                    paste("A",t,sep="."),
                    ag,
                    c(0))
      
      for(u in graph::edges(g)[[t]] ) {
        if( u == t) next
        ## s.t to t.access
        
        ## s.t to t.u

        ag <- addEdge(paste(s,t,l,sep="."),
                      paste(t,u,l,sep="."),
                      ag,
                      c(0))
      }
    }
  }

  ag <- rg.set.capacity(ag,lambdacap)
  graph <- list()
  graph$g <- g
  graph$ag <- ag
  graph$anodes <- accessnodes
  return(graph)
}


rg.sp.lambda.max.concurrent.flow <- function(g,demands) {
  for(c in names(demands)) {
    demands[[c]]$flow <- 0
  }

  nodes <- nodes(g)
  gsol <- g
  punults.done <- rep(FALSE,length(nodes))

  rg.set.all.graph.edge.weights(gsol)
  
  results <- list()
  results$g <- g
  results$demands <- demands

  ccount <- 1
  for (i in demands) {
    startind=which(nodes==i$source)
    finind=which(nodes==i$sink)
    if(!punults.done[startind]) {
      g.sp[[startind]] <- dijkstra.sp(g,start=i$source)$penult
      punults.done[startind]=TRUE
    }
    path <- extractPath(startind,finind,g.sp[[startind]])
    print(path)
    gsol <- rg.addto.weight.on.path(gsol,path,i$demand)
    demands[[ccount]]$flow <- demands[[ccount]]$flow +i$demand
    ccount <- ccount + 1
  }
  f <- as.double(edgeData(gsol,attr="weight"))
  c <- as.double(edgeData(gsol,attr="capacity"))
  lambda <- min(c/f)

  gamma <- min((c-f)/c)
    
  #demands <- rg.max.concurrent.flow.rescale.demands.flow(demands,lambda)
    
  #gsol <- rg.max.concurrent.flow.graph(gsol,demands)
  
  results$demands <- demands
  results$gflow <- gsol
  results$lambda <- lambda
  results$gamma <- gamma
    
  return(results)
}

rg.dfs <- function(g) {
  data <- list()
  data$color <- rep(0,length(nodes(g))+1)
  data$pred <- rep(NULL,length(nodes(g)))
  for (u in nodes(g)) {
    cat("I am visiting", u,"\n")
    if(data$color[as.integer(u)] == 0 ) {
      data <- mydfs.visit(g,u,data)
    }
  }
}

rg.dfs.visit <- function(g,u,data) {
  data$color[as.integer(u)] <- 1
  print(data$color[as.integer(u)])
  for (v in graph::edges(g,u)[[1]]) {
    cat("testing edge:",v,":\n")
    if(data$color[as.integer(v)] == 0 ) {
      data$pred[as.integer(v)] <- as.integer(u)
      data <- mydfs.visit(g,v,data)
    }
  }
  
  data$color[as.integer(u)] <- 2
  print(data$color[as.integer(u)])
  cat("I have finished at ",u,"\n")
  return(data)
}


analyse.runs <- function(nums=NULL,maxattempts=1,progress=FALSE,e=0.005) {

  if(is.null(nums)) {
    files <- Sys.glob("run[0-9]*")
  } else {
    files <- c()
    for(num in nums) {
      files <- c(files,Sys.glob(paste("run",num,".*",sep="")))
    }
  }

  files <- files[order(as.numeric(gsub("^.*run([0-9]*).*$","\\1",files,files)))]
  foundratio <- as.double(c())
  countzero <- 0
  neverfound <- 0
  faillist <- c()
  numpaths <- 0
  countruns <- 0
  for(file in files) {
##  for(i in seq(1:50)) {
##    file <- paste("run",i,".robj",sep="")
    cat("----------------- Doing ",file,"\n")
    load(file)
    if(is.null(res$error)) {
      scenario <- res$scenario
      if(scenario$intgammalist[[1]] > 0.0) {

        attempts <- 0
        countgamma <- 0
        while(countgamma == 0 && attempts < maxattempts) {
          retlist <- analyse.inner(scenario,foundratio,file,progress=progress,e=e)
          countgamma <- retlist$res$countgamma
          attempts <- attempts + 1
          if(retlist$res$countgamma == 0)
            countzero <- countzero +1
        }
        if(retlist$res$countgamma == 0) {
          neverfound <- neverfound +1
          faillist <- c(file,faillist)
        }
        countruns <- countruns + 1
        numpaths <- numpaths + retlist$numpaths

        foundratio <- retlist$foundratio
      }
    } else {
      ##cat("was null\n")
    }
  }
  cat("numzero=",countzero,
      "min=",min(foundratio)," mean=",mean(foundratio)," max=",max(foundratio),
      "\n")
  cat("av numpaths=",numpaths/countruns,"\n")
  if(neverfound >0) {
    cat("Did not find",neverfound," for:\n")
    cat(gsub("[^0-9]+([0-9]+)[^0-9]+","\\1",faillist),"\n")
  }
  
  

}

analyse.inner.part1 <- function(nums,e=0.005) {

  if(is.null(nums)) {
    files <- Sys.glob("run[0-9]*")
  } else {
    files <- c()
    for(num in nums) {
      files <- c(files,Sys.glob(paste("run",num,".*",sep="")))
    }
  }
  
  files <- files[order(as.numeric(gsub("^.*run([0-9]*).*$","\\1",files,files)))]
  for(file in files) {
##  for(i in seq(1:50)) {
##    file <- paste("run",i,".robj",sep="")
    cat("----------------- Doing ",file,"\n")
    load(file)
    if(is.null(res$error)) {
      scenario <- res$scenario
      if(scenario$intgammalist[[1]] > 0.0) {
        results <- rg.minimum.congestion.flow(scenario$g,scenario$commr,e=e,progress=TRUE,permutation="lowest")

      }
    }
  }
  results$scenario <- scenario
  return(results)
}

analyse.inner <- function(scenario,foundratio,file,progress,e) {
  
  results <- rg.minimum.congestion.flow(scenario$g,scenario$commr,e=e,progress=progress,permutation="lowest")
  #res <- rg.try.single.demands(scenario$g,scenario$commr,e=0.005,permutation="lowest",
  #                             scenario=scenario,progress=progress)
  #numpaths <- rg.count.paths(res$demands)

  numpaths <- rg.count.paths(results$demands)

  res <- rg.max.concurrent.flow.int.c(
                                      scenario$g,
                                      results$demands,
                                      e=e,
                                      #eInternal=0.001,
                                      scenario=scenario,
                                      progress=progress,
                                      permutation="lowest"
                                      )
  res$scenario <- scenario
  foundratio <- c(foundratio,as.double(res$countgamma)/as.double(res$phases))
  cat(file,":",res$countgamma,"/",res$phases,
      " gamma=",res$scenario$intgammalist[[1]],
      "best found=",res$bestgamma,"\n")
  cat("pathdifcount=",res$pathdiffcount,"\n")
  cat("phasepathdiffcount=",res$phasepathdiffcount,"\n")
  retlist <- list()
  retlist$res <- res
  retlist$foundratio <- foundratio
  retlist$numpaths <- numpaths
  return(retlist)
}

rg.try.single.demands <- function(g,demands,e=0.1,progress=FALSE,permutation="fixed",scenario) {

  newdemands <- list()
  for(i in 1:length(demands)) {
    tmp <- list()
    tmp[["1"]] <- demands[[i]]
    results <- rg.minimum.congestion.flow(g,tmp,e=0.01,progress=progress,permutation=permutation)
    newdemands[[i]] <- results$demands[[1]]
  }
  res <- rg.max.concurrent.flow.int.c(g,newdemands,e=e,progress=progress,scenario=scenario,permutation=permutation)
  return(res)
}





rg.test.idea <- function(maxtime=1.0,N=9,L=10,filebase=NULL) {
  ## maximum time in hours for a run
  ## this is based upon 24e6 patch combinations taking 106 seconds
  hoursperpath <- 1470 / ( 3600 * 253e6  )
  maxpathcomb <- maxtime / hoursperpath

  filename <- paste(filebase,".robj",sep="")
  scenario <- rg.gen.random.scenario(N=N,L=L)

  starttime <- Sys.time()
  cat("starting at ",format(starttime),"\n")

  
  
  if(scenario$pathcomb > maxpathcomb) {
    cat("scenario =",scenario$pathcomb," too big\n")
#        as.double(maxpathcomb)/as.double(scenario$pathcomb) * maxtime,
    
    cat("it would have taken ",
        scenario$pathcomb * hoursperpath,
        "hours\n")
    res <- list()
    res$scenario <- scenario
    res$error <- "Scenario too big"
    if(!is.null(filebase)) {
      save(res,file=filename)
    }
    
    return(res)
  }
  estfinished <- starttime + scenario$pathcomb * hoursperpath * 3600
  cat("Estimate finish at ",format(estfinished),"\n")
  scenario <- rg.testeverypath.c(scenario)
  res <-
    rg.fleischer.max.concurrent.flow.restricted.c(
                                                  scenario$g,
                                                  scenario$results$demands,
                                                  e=0.01,
                                                  bestgamma=
                                                  scenario$intgammalist[[1]])
  res$scenario <- scenario
  cat("N=",N," L=",L," pathcomb=",scenario$pathcomb,
      "gammahits=",res$countgamma,"out of ",res$phases," phases\n")
  
  if(!is.null(filebase)) {
    save(res,file=filename)
  }
  
  return(res)
}


### Generate a random graphNEL and scaled commodities so that gamma is 0.9
### N    Number of nodes in graph
### L    Number of commodities
### returns output$pathcomb - number of combinations of paths
###         output$g  - graphNEL with capacity
###         output$results - from running rg.minimum.congestion.flow()
###         output$commr - the commodities (scaled to give gamma a reasonable value)
###         output$comm - the original commodities
rg.gen.random.scenario <- function(N=6,L=8,target=0.9,progress=FALSE,e=0.02) {
  g <- rg.generate.random.graph(N,3,6)
  M <- length(edgeMatrix(g))/2
  comm <- rg.gen.demands(g,L,runif(L,1,5))
  g <- rg.set.capacity(g,runif(M,5,10))
  g <- rg.set.weight(g,0.0)
  results <- rg.minimum.congestion.flow(g,comm,e=e,progress=progress)
  ## need to customize this
  multiplier <- (1-target)/(1-results$gamma)
  commr <- rg.rescale.demands(comm,multiplier)
  results <- rg.minimum.congestion.flow(g,commr,e=e,progress=progress)
  runtotal <- 1
  for(i in results$demands) {
    val <- length(names(i$paths))
    runtotal <- val * runtotal
  }
  output <- list()
  output$pathcomb <- runtotal
  output$g <- g
  output$results <- results
  output$comm <- comm
  output$commr <- commr
  return(output)
}

### Brute force search for the best path. It returns data with new values
### input data$g - graph graphNEL with capacity edge attribute set
###       data$results data$g data$comm data$commr data$pathcomb
###       this is as from rg.gen.random.scenario()
### returns data as above plus
###         data$intdemands
###         data$gintflow
###         WARNING: it will overboook
###         WARNING: it seems to keep the edges attribute from the graph
###                  it copied, this should be set to null.

rg.testeverypath.c <- function(data,progress=FALSE) {
  recordlen <- 100
  demands <- data$results$demands
  intdemands <-  demands

  g <- data$g
  vdemands <- as.double(lapply(demands,"[[","demand"))
  vcapacity <- as.double(edgeData(g,attr="capacity"))
  vdemandflow <- rep(0.0,length(vcapacity))
  L <- length(demands)
  
  edgeMap <- adjMatrix(g)
  demandpaths <- list()
  j <- 1
  for(d in demands) {
    demandpaths[[j]] <- list()
    k <- 1
    for(p in names(d$paths)) {
      pv <- as.vector(strsplit(p,"|",fixed=TRUE)[[1]])
      fromlist <- as.integer(pv[1:length(pv)-1])
      tolist <- as.integer(pv[2:length(pv)])
      me <- rbind(fromlist,tolist)
      pv <- c()
      for(i in 1:length(fromlist)) {
        v <- as.vector(me[,i])
        em <- edgeMap[v[1],v[2]]
        pv <- append(pv,em)
      }
      demandpaths[[j]][[k]] <- pv -1
      k <- k + 1
    }
    j <- j + 1
  }

  m <- length(rg.edgeL(g))

  if(progress != FALSE) {
    pb <- txtProgressBar(title = "progress bar", min = 0,
                         max = data$pathcomb, style=3)
  } else {
    pb <- NULL
  }

  retlist <- .Call("rg_test_every_path_inner",
                   demandpaths,
                   vdemands,
                   vcapacity,
                   progress,
                   pb,
                   recordlen,
                   environment()
                   );

  if( progress != FALSE) {
    close(pb)
    cat("Wrapping up and creating ",recordlen, "best graphs\n")
  }

  data$gintflowlist <- list()
  data$intdemandslist <- list()
  data$intgammalist <- retlist$gammas
  recorddemands <- retlist$pathptrs
  for(j in seq(1:recordlen) ) {
    g <- data$g
    g <- rg.set.weight(g,0)
    bestPathptr <- recorddemands[[j]] + 1
#    cat("test1\n")
#    print(bestPathptr[[i]])
#    print(demands[[i]])
#    cat("test2\n")
    for(i in seq(1:L)) {
      intdemands[[i]]$paths <- list()
      intdemands[[i]]$flow <- demands[[i]]$flow
      intdemands[[i]]$paths[[names(demands[[i]]$paths[bestPathptr[[i]]])]] <-
        demands[[i]]$flow
      path <- names(demands[[i]]$paths)[ bestPathptr[[i]] ]
      pv <- as.vector(strsplit(path,"|",fixed=TRUE)[[1]])
      
      val <- demands[[i]]$demand
      
      if(length(pv) == 1) {
        edgeData(g,as.character(pv[1]),as.character(pv[2]),"weight") <-
          edgeData(g,as.character(pv[1]),as.character(pv[2]),"weight") +val
      } else {
        
        fromlist <- pv[1:{length(pv)-1}]
        tolist <- pv[2:{length(pv)}]
        newvals <- as.double(edgeData(g,as.character(fromlist),as.character(tolist),"weight")) + val
        edgeData(g,as.character(fromlist),as.character(tolist),"weight") <- newvals
      }
      
    }

    data$gintflowlist[[j]] <- g

    data$intdemandslist[[j]] <- intdemands
  }
  data$gintflow <- data$gintflowlist[[1]]
  data$intdemands <- data$intdemandslist[[1]]
  return(data)


}


rg.testeverypath <- function(data,progress=FALSE) {

  ## record this many of the best, sorted by best gamma
  recordlen <- 100
  recordgamma <- rep(-Inf,recordlen)
  recorddemands <- rep(list(NULL),recordlen)

  L=length(data$results$demands)
  pathptr <- rep(1,times=L)
  pathsz <- c()
  demands <-  data$results$demands
  count <- 1
  for(i in demands) {
    pathsz[count] <- length(i$paths)
    count <- count +1
  }
  finished=FALSE
  count=1
  intdemands <-  demands
  for(d in names(intdemands)) {
    intdemands[[d]]$paths <- list()
    intdemands[[d]]$flow <- 0.0
  }
  intgraph <- data$g
  bestgamma <- -Inf
  bestPathptr <- c()
  count <- 0

  if(progress != FALSE) {
    pb <- txtProgressBar(min = 0,
                         max = data$pathcomb, style=3)
  } else {
    pb <- NULL
  }
  updatepb <- as.integer(ceiling( data$pathcomb / 100.0))

  while(!finished) {

    if(progress != FALSE && count %% updatepb == 0) {
      setTxtProgressBar(pb,count)
    }

    #print(pathptr)

    ## pathptr is counter to path selectors

    ## add load to each path
    ## for each demand (index)
    g <- data$g
    
    for(i in seq(1:L)) {
      path <- names(demands[[i]]$paths)[ pathptr[i] ]
      pv <- as.vector(strsplit(path,"|",fixed=TRUE)[[1]])

#     print(path)
      val <- demands[[i]]$demand
      
      if(length(pv) == 1) {
        edgeData(g,as.character(pv[1]),as.character(pv[2]),"weight") <-
          edgeData(g,as.character(pv[1]),as.character(pv[2]),"weight") +val
      } else {
        
        fromlist <- pv[1:{length(pv)-1}]
        tolist <- pv[2:{length(pv)}]
        newvals <- as.double(edgeData(g,as.character(fromlist),as.character(tolist),"weight")) + val
        edgeData(g,as.character(fromlist),as.character(tolist),"weight") <- newvals
      }
      
    }
    
    ## calc max load
    weights <- as.double(edgeData(g,attr="weight"))
    capacities <- as.double(edgeData(g,attr="capacity"))
    gamma <- min((capacities - weights)/capacities)
    ## record best gamma and pathptr so far
    if( gamma > bestgamma ) {
      bestgamma <- gamma
      bestPathptr <- pathptr
      cat("best gamma so far =",bestgamma,"\n")
    }
    if( gamma > recordgamma[recordlen] ) {
      recordgamma[recordlen] <- gamma
      recorddemands[[recordlen]] <- pathptr
      o <- order(recordgamma,decreasing=TRUE)
      recordgamma <- recordgamma[o]
      recorddemands <- recorddemands[o]
      
    }
    #cat("found=",gamma," best=",bestgamma,"\n")

    #cat(recordgamma,"\n")
    increment=TRUE
    for(i in seq(1:L)) {
      if(pathsz[i]>1) {
        if(pathptr[i]%%pathsz[i] == 0)
          incrementnext <- TRUE
        else
          incrementnext <- FALSE
        if(increment)
          pathptr[i] <- pathptr[i]%%pathsz[i]+1
        if(incrementnext && increment)
          increment <- TRUE
        else
          increment <- FALSE
      }
    }

    
    if(increment)
      finished <- TRUE
    count <-  count + 1
    # just for testing
    #if(count > 20)
     # finished=TRUE
  }

  if(progress != FALSE)
    close(pb)

  ## have found best, now fill integer flow graph and demands

  data$gintflowlist <- list()
  data$intdemandslist <- list()
  data$intgammalist <- recordgamma
#  cat("bestpathptr=",recorddemands[[1]],"\n")
  for(j in seq(1:recordlen) ) {
    g <- data$g

    bestPathptr <- recorddemands[[j]]
    for(i in seq(1:L)) {
      intdemands[[i]]$flow <- demands[[i]]$flow
      intdemands[[i]]$paths <- NULL
      intdemands[[i]]$paths[[names(demands[[i]]$paths[bestPathptr[i]])]] <-
        demands[[i]]$flow
      path <- names(demands[[i]]$paths)[ bestPathptr[i] ]
      pv <- as.vector(strsplit(path,"|",fixed=TRUE)[[1]])
      
      val <- demands[[i]]$demand
      
      if(length(pv) == 1) {
        edgeData(g,as.character(pv[1]),as.character(pv[2]),"weight") <-
          edgeData(g,as.character(pv[1]),as.character(pv[2]),"weight") +val
      } else {
        
        fromlist <- pv[1:{length(pv)-1}]
        tolist <- pv[2:{length(pv)}]
        newvals <- as.double(edgeData(g,as.character(fromlist),as.character(tolist),"weight")) + val
        edgeData(g,as.character(fromlist),as.character(tolist),"weight") <- newvals
      }
      
    }

    data$gintflowlist[[j]] <- g

    data$intdemandslist[[j]] <- intdemands
  }
  data$gintflow <- data$gintflowlist[[1]]
  data$intdemands <- data$intdemandslist[[1]]
  return(data)
}

calcDexplicit <- function(gdual) {
    return(sum(as.double(edgeData(gdual,attr="capacity"))*
        as.double(edgeData(gdual,attr="weight"))))
  }



###system.time(test1(gdual))
###   user  system elapsed 
### 44.084   0.283  44.589 
test1 <- function(g) {

  var <- 1.0000001
  val <- 1.0
  em <- edgeMatrix(g)


  
  for(i in seq(1:10000)) {
    w <- as.double(edgeData(g,attr="weight"))
    w <- w*var
    g <- rg.set.weight(g,w)
  }

  return(g)


}

###system.time(test2(gdual))
###   user  system elapsed 
### 43.785   0.214  44.427 
test2 <- function(g) {

  var <- 1.0000001
  val <- 1.0
  em <- edgeMatrix(g)


  
  for(i in seq(1:10000)) {
    w <- as.double(edgeData(g,attr="weight"))
    w <- w*var
    em <- edgeMatrix(g)
    if (match("weight",names(edgeDataDefaults(g)),nomatch=1)) {
      edgeDataDefaults(g,"weight") <- 1.0
    }
    edgeData(g,
             from=as.character(em[1,]),
             to=as.character(em[2,]),
             attr="weight") <- val
  }
  
  return(g)


}


rg.test2 <- function() {
  pb <- txtProgressBar(min=0,max=5,style=3)
  for(i in 1:5) {
    Sys.sleep(1)
    setTxtProgressBar(pb,i)
    Sys.sleep(1)
    close(pb)
  }
}

rg.explore.start <- function(g,demands,e=0.1,progress=FALSE,intdemands=NULL) {

  initresults <- rg.fleischer.max.concurrent.flow.c(g,demands,e=e)
  cat("initresult$lambda=",initresults$lambda,
      "initresults$beta=",initresults$beta,"\n")
  initdemands <- initresults$demands
#  res <- rg.explore(g,initdemands,e=e,progress=progress,intdemands=intdemands)

  res <- rg.explore.outer(g,initdemands,e=e,progress=progress,intdemands=intdemands)

  ## fix flow to equal demand
  for(d in 1:length(res)) {
    res[[d]]$paths[1] <- res[[d]]$demand
    res[[d]]$flow <- res[[d]]$demand
  }

  gint <- rg.max.concurrent.flow.graph(g,intdemands)
  gexplore <- rg.max.concurrent.flow.graph(g,res)

  gammaint <- rg.mcf.find.gamma(gint)
  gammaexplore <- rg.mcf.find.gamma(gexplore)


  for(d in 1:length(res)) {
    cat(names(res[[d]]$paths[1]), ", ",
        names(intdemands[[d]]$paths[1]))
    if( identical (names(res[[d]]$paths[1]),
                   names(intdemands[[d]]$paths[1]))) {
      cat("  OK \n")
    } else {
      cat(" !!!!!! Different \n")
    }      
        

  }
  
  cat("Int gamma=",gammaint," gexplore gamma=",gammaexplore,"\n")
  return(res)

}
rg.explore.outer <- function(g,demands,e=0.1,progress=FALSE,intdemands=NULL) {

  res <- rg.explore(g,demands,e=e,progress=progress,intdemands=intdemands)

  finished <- TRUE
  for(d in 1:length(res)) {
    if (length(res[[d]]$paths) > 1 )
      finished <- FALSE
    cat("num paths in ",d," = ",length(res[[d]]$paths),"\n")
  }

  if( ! finished ) {
    res <- rg.explore.outer(g,res,e=e,progress=progress,intdemands=intdemands)
  }
  
  return(res)
}


rg.explore <- function(g,initdemands,e=0.1,progress=FALSE,intdemands=NULL) {
  ## calc split flows

  builtdemands <- initdemands
  
  ## this line not strictly needed
  #res <- rg.fleischer.max.concurrent.flow.restricted.c(g,initdemands,
  #                                                   updateflow=FALSE,e=e)

  #cat("res$beta=",res$beta," res$betar=",res$betar,
  #    " res$lambda=",res$lambda,"\n")

  ## try one path at a time in each demand and calculate lambda, record maximum

  overallMaxMetric <- -Inf
  overallMaxMetricD <- 0

  demandpaths <- list()
  demandpathsbest <- list()
  j <- 1
  for(d in initdemands) {
    demandpaths[[j]] <- list()
    demandpathsbest[[j]] <- 0
    k <- 1
    for(p in names(d$paths)) {
      demandpaths[[j]][[k]] <- list()
      demandpaths[[j]][[k]]$name <- p
      demandpaths[[j]][[k]]$lambda <- 0.0
      demandpaths[[j]][[k]]$beta <- 0.0
      k <- k + 1
    }
    j <- j + 1
  }
  
  for(j in 1:length(initdemands)) {
    maxmetric <- -Inf
    maxmetrici <- i
    cat("demand number",j," num paths ",length(initdemands[[j]]$paths),"---------------\n")
    if (length(initdemands[[j]]$paths) > 1 ) {
      for(i in 1:length(initdemands[[j]]$paths)) {
        d <- initdemands

        ## remove the current path
        cat("removing path ",names(d[[j]]$paths[i]))
        #cat(" flow= ",d[[j]]$paths[[i]])
        d[[j]]$paths[i] <- NULL
        ##cat(names(d[[j]]$paths))
        res <- rg.max.concurrent.flow.capacity.restricted.c(g,d,e=e)

        gamma <- rg.mcf.find.gamma(res$gflow,res$lambda)
        
        #metric <- res$lambda * res$betar / initdemands[[j]]$paths[[i]]
        metric <- gamma
        cat(" metric=",metric," res$betar=",res$betar,
            " res$lambda=",res$lambda,"\n")
        demandpaths[[j]][[i]]$lambda <- res$lambda
        demandpaths[[j]][[i]]$beta <- res$betar

        if ( metric > maxmetric ) {
          maxmetric <- metric
          maxmetrici <- i
          demandpathsbest[[j]] <- i
        }

      }
      if(maxmetric > overallMaxMetric ) {
        cat("better metric by ",overallMaxMetric - maxmetric,"\n")
        overallMaxMetric <- maxmetric
        overallMaxMetricD <- j
      } 
    } else {
      cat("  skipping, only one path",names(initdemands[[j]]$paths),"\n")
    }


    ## summarise results
    cat("max metric=",maxmetric," in path",
        names(initdemands[[j]]$paths[maxmetrici]),
        " path number ",maxmetrici)
    if( identical(names(initdemands[[j]]$paths[maxmetrici]),
                  names(intdemands[[j]]$paths[1]))) {
      cat(" !! in intdemands\n") }
    else cat("\n")
    
  }

  bestpathtoremove <- demandpathsbest[[overallMaxMetricD]]
  cat("----- summary-----------------------------------------------\n")
  cat("overall best metric is for demand ",overallMaxMetricD,", path ",
      bestpathtoremove,"\n")

  nametoremove <- names(builtdemands[[overallMaxMetricD]]$paths[bestpathtoremove])

  builtdemands[[overallMaxMetricD]]$paths[bestpathtoremove] <- NULL
  
  cat(" removing path ",nametoremove,"\n")
  cat("  integer demand is ",names(intdemands[[overallMaxMetricD]]$paths[1]),"\n")
  if( identical(names(intdemands[[overallMaxMetricD]]$paths[1]),
                nametoremove)) {
    cat("ERROR removing INTEGER DEMAND\n")
  }
  return(builtdemands)
}

## note demands must be the result of running
##     initresults <- rg.fleischer.max.concurrent.flow.c(g,demands,e=e)
## demands <- initresults$demands
rg.explore2 <- function(g,initdemands,e=0.1,progress=FALSE,intdemands=NULL) {
  ## calc split flows

  builtdemands <- initdemands
  
  ## this line not strictly needed
  #res <- rg.fleischer.max.concurrent.flow.restricted.c(g,initdemands,
  #                                                   updateflow=FALSE,e=e)

  #cat("res$beta=",res$beta," res$betar=",res$betar,
  #    " res$lambda=",res$lambda,"\n")

  ## try one path at a time in each demand and calculate lambda, record maximum

  overallMaxLambda <- -Inf
  overallMaxLambdaD <- 0
  overallMaxLambdaBeta <- -Inf
  bestSumBetaLambda <- -Inf
  bestSumBetaLambdaD <- 0
  worstFlow <- Inf
  worstFlowD <- 0
  worstFlowi <- 0

  demandpaths <- list()
  demandpathsbest <- list()
  j <- 1
  for(d in initdemands) {
    demandpaths[[j]] <- list()
    demandpathsbest[[j]] <- 0
    k <- 1
    for(p in names(d$paths)) {
      demandpaths[[j]][[k]] <- list()
      demandpaths[[j]][[k]]$name <- p
      demandpaths[[j]][[k]]$lambda <- 0.0
      demandpaths[[j]][[k]]$beta <- 0.0
      k <- k + 1
    }
    j <- j + 1
  }

  res <- rg.max.concurrent.flow.capacity.restricted.c(g,initdemands,e=e)
  cat(" res$betar=",res$betar,
      " res$lambda=",res$lambda,"\n")

  d <- res$demands
  
  for(j in 1:length(initdemands)) {
    minlambda <- Inf
    maxlambda <- -Inf
    betaAtMaxLambda <- -Inf
    maxbeta <- -Inf
    minbeta <- Inf
    maxlambdai <- 0
    maxbetai <- 0
    cat("demand number",j," num paths ",length(initdemands[[j]]$paths),"---------------\n")
    if (length(d[[j]]$paths) > 1 ) {
      for(i in 1:length(d[[j]]$paths)) {
        cat(names(d[[j]]$paths[i])," flow=",
            d[[j]]$paths[[i]],"\n")
        if (worstFlow > d[[j]]$paths[[i]] / d[[j]]$flow ) {
          worstFlow <- d[[j]]$paths[[i]] / d[[j]]$flow
          worstFlowD <- j
          worstFlowi <- i
          cat("NEW worst flow\n")
        }
      }
    } else {
      cat("  skipping, only one path",names(d[[j]]$paths),"\n")
    }
    
  }
  
  bestpathtoremove <- names(d[[worstFlowD]]$paths[worstFlowi])
  cat("----- summary-----------------------------------------------\n")
  cat("overall worst flow is for demand ",worstFlowD,", path ",
      bestpathtoremove,"\n")



  builtdemands[[worstFlowD]]$paths[bestpathtoremove] <- NULL
  
  cat(" removing path ",bestpathtoremove,"\n")
  cat("  integer demand is ",names(intdemands[[worstFlowD]]$paths[1]),"\n")
  if( identical(names(intdemands[[worstFlowD]]$paths[1]),
                bestpathtoremove)) {
    cat("ERROR removing INTEGER DEMAND\n")
  }
  return(builtdemands)
}

### Same as rg.fleischer.max.concurrent.flow.restricted()
### except that it tests to see if the integerdemands are in
### the demand paths used in a single phase.
### When running this on a small example it gave the suprising result
### that the integer solution paths were never used in the phase.
### Out of the eight solution paths only seven were ever the same! This
### was for seven out of the many (40000?) phases.
rg.fleischer.max.concurrent.flow.restricted.test <- function(g,demands,e=0.1,updateflow=TRUE, progress=FALSE,integerdemands=NULL,bestgamma=NULL) {

  findshortestpath <- function(paths,vlength) {
    min <- Inf
    minp <- NULL
    minpcount <- NULL
    for(i in 1:length(paths)) {
      if ( sum(vlength[paths[[i]]]) < min) {
        minp <- paths[[i]]
        minpcount <- i
        min <- sum(vlength[paths[[i]]])
      }
    }
    
    return(minpcount)
  }
  
  savedemands <- demands
  
  vdemands <- as.double(lapply(demands,"[[","demand"))
  vcapacity <- as.double(edgeData(g,attr="capacity"))
  vweight <- rep(0.0,length(vcapacity))
  vlength <- rep(0.0,length(vcapacity))
  vdemandflow <- rep(0.0,length(vcapacity))
  L <- length(vlength)
  M <- length(vdemands)
  edgeMap <- adjMatrix(g)

  demandpaths <- list()
  demandpathflows <- list()
  integerdemandsindex <- list()
  j <- 1
  for(d in demands) {
    demandpaths[[j]] <- list()
    demandpathflows[[j]] <- list()
    k <- 1
    for(p in names(d$paths)) {
      pv <- as.vector(strsplit(p,"|",fixed=TRUE)[[1]])
      fromlist <- as.integer(pv[1:length(pv)-1])
      tolist <- as.integer(pv[2:length(pv)])
      me <- rbind(fromlist,tolist)
      pv <- c()
      for(i in 1:length(fromlist)) {
        v <- as.vector(me[,i])
        em <- edgeMap[v[1],v[2]]
        pv <- append(pv,em)
      }
      demandpaths[[j]][[k]] <- pv
      if( identical(p, names(integerdemands[[j]]$paths[1]))) {
        integerdemandsindex[[j]] <- k
      }
      demandpathflows[[j]][[k]] <- 0.0
      k <- k + 1
    }
    j <- j + 1
  }

  ## note this is not the dual value
  calcD <- function() {
    sum(vcapacity*vlength)
  }
  
  doubleCount <- 0
  doubleDemands <- function(demands) {
    vdemands <- vdemands * 2.0
    doubleCount <- doubleCount + 1
    vdemands
  }
  gdual <- g

  ## number of arcs
  m <- length(rg.edgeL(g))
  ## number of nodes
  N <- length(nodes(g))
  delta <- (m / (1-e)) ^ (-1/e)
  vlength <- delta / vcapacity

  for(c in names(demands)) {
    demands[[c]]$paths <- list()
    demands[[c]]$flow <- 0
  }

  
  doubreq <- 2 * ceiling(1/e * log(m/(1-e),base=(1+e)))
  ## updatepb used to measure progress as percent
  updatepb <- as.integer(ceiling(doubreq / 100.0))

  if(progress != FALSE)
    pb <- txtProgressBar(min = 0,
                        max = doubreq, style=3)
  phases <- 0
  totalphases <- 0
  countallint <- 0
  countgamma <- 0
  D <- calcD()

  pdemands <- c(length(vdemands),1:(length(vdemands) -1 ))
  demseq <- 1:length(vdemands)

  ## main algorithm
  while(D < 1 ) {
    if(progress != FALSE && totalphases %% updatepb == 0) {
      setTxtProgressBar(pb,totalphases)
    }
    if(phases > doubreq) {
      cat("doubling",doubreq,totalphases,"\n");
      vdemands <- doubleDemands(vdemands)
      phases <- 0
    }

    allinteger <- 0
    demseq <- demseq[pdemands]
    vcaptmp <- vcapacity
    vweights <- rep(0.0,L)
    underdemands <- TRUE

    tmppaths <- list()
    
    for(i in demseq) {
      demand <- vdemands[i]
      D <- calcD()
      while( D < 1 && demand > 0 ) {
        p <- findshortestpath(demandpaths[[i]],vlength)
        ##cat(integerdemandsindex[[i]]," p=",p,"\n")
        tmppaths[[i]] <- as.integer(p)
        if ( identical(as.integer(p),as.integer(integerdemandsindex[[i]])) &&
            allinteger >= 0 ) {
          allinteger <- allinteger + 1
        }
        vweights[ demandpaths[[i]][[p]] ] <- vweights[ demandpaths[[i]][[p]] ] +
          demand
        caponpath <- vcaptmp[ demandpaths[[i]][[p]] ]
        mincap <- min(demand,caponpath)
        if(mincap < demand)
         underdemands <- FALSE
        demand <- demand - mincap
        lengths <- vlength[ demandpaths[[i]][[p]] ]
        lengths <- lengths * (1 + (e*mincap) / vcapacity[ demandpaths[[i]][[p]] ])
        vlength[ demandpaths[[i]][[p]] ] <- lengths
        demandpathflows[[i]][[p]] <- demandpathflows[[i]][[p]] + mincap
        vdemandflow[i] <- vdemandflow[i] + mincap
        #vcaptmp[ demandpaths[[i]][[p]] ] <-
        #  vcaptmp[ demandpaths[[i]][[p]] ] - mincap
        
        
        D <-  calcD()
      }
      
    }
    gamma = min( 1 - vweights/vcapacity )
#    cat("gamma=",gamma,"\n")
    if(underdemands && D < 1) {
      if(allinteger == M) {
        countallint <- countallint +1
      }
      if(gamma >= bestgamma * 0.99999) {
        countgamma <- countgamma + 1

      }
    }
#    cat("NUM INTEGER PATHS",allinteger,"\n")

    phases <- phases +1
    totalphases <- totalphases + 1
  }
  ## end of main algorithm - now need to pack results

  ## If there had been demands doubling we need to
  ## fix things for the returned demands
  for(n in 1:length(demands)) {
    demands[[n]]$demand <- savedemands[[n]]$demand
    demands[[n]]$paths <- savedemands[[n]]$paths
    demands[[n]]$flow <- vdemandflow[n]
    for(i in 1:length(demands[[n]]$paths)) {
      demands[[n]]$paths[[i]] <- demandpathflows[[n]][[i]]
    }
  }
  gdual <- rg.set.weight(gdual,vlength)
  scalef <- 1 / log(1/delta,base=(1+e))
  demands <- rg.max.concurrent.flow.rescale.demands.flow(demands,scalef)

  beta <- calcBeta(demands,gdual)
  betar <- calcBetaRestricted(demands,gdual)
  lambda=NULL
  lambda <- calcLambda(demands)
  foundratio <- beta / lambda
  ratiobound <- (1-e)^-3
    
  gflow <- rg.max.concurrent.flow.graph(gdual,demands)
  retval <- list(demands=demands,gflow=gflow,gdual=gdual,beta=beta,betar=betar,
                 lambda=lambda,phases=totalphases,e=e,vlength=vlength,
                 countallint=countallint,countgamma=countgamma)

  if(progress != FALSE)
    close(pb)

  retval
}


rg.max.concurrent.flow.capacity.restricted.c <- function(g,demands,e=0.1,updateflow=TRUE,progress=FALSE) {

  ## note this is not the dual value
  calcD <- function() {
    sum(vcapacity*vlength)
  }

  savedemands <- demands
  
  vdemands <- as.double(lapply(demands,"[[","demand"))
  vcapacity <- as.double(edgeData(g,attr="capacity"))
  vlength <- rep(0.0,length(vcapacity))
  vdemandflow <- rep(0.0,length(vcapacity))
  L <- length(vlength)
  delta <- (L / (1-e)) ^ (-1/e)
  vlength <- delta / vcapacity
  edgeMap <- adjMatrix(g)

  demands.paths <- as.integer(c())

  ## code the paths as integers
  ## [demandno(e.g=1) numdempaths lengthpath1 path1[1] path1[2] ...
  ##  lengthpath2 path2[1] path2[2]....demandno(e.g=2)....]
  rdemandpaths <- list();
  for(d in 1:length(demands)) {
    rdemandpaths[[d]] <- list()
    demands.paths <- append(demands.paths,d-1)
    demands.paths <- append(demands.paths,length(demands[[d]]$paths))
    for(p in 1:length(demands[[d]]$paths)) {
      pathi <-
        as.integer(strsplit
                   (names(demands[[d]]$paths[p]),"|",fixed=TRUE)[[1]])
      pathi <- pathi -1
      rdemandpaths[[d]][[p]] <- pathi
      demands.paths <- append(demands.paths,length(pathi))
      demands.paths <- append(demands.paths,pathi)
    }

  }

  demandpaths <- list()
  j <- 1
  for(d in demands) {
    demandpaths[[j]] <- list()
    k <- 1
    for(p in names(d$paths)) {
      pv <- as.vector(strsplit(p,"|",fixed=TRUE)[[1]])
      fromlist <- as.integer(pv[1:length(pv)-1])
      tolist <- as.integer(pv[2:length(pv)])
      me <- rbind(fromlist,tolist)
      pv <- c()
      for(i in 1:length(fromlist)) {
        v <- as.vector(me[,i])
        em <- edgeMap[v[1],v[2]]
        pv <- append(pv,em)
      }
      demandpaths[[j]][[k]] <- pv -1
      k <- k + 1
    }
    j <- j + 1
  }

  
  
  m <- length(rg.edgeL(g))

  doubreq <- 2 * ceiling(1/e * log(m/(1-e),base=(1+e)))

  if(progress != FALSE) {
    pb <- txtProgressBar(title = "progress bar", min = 0,
                         max = doubreq, style=3)
  } else {
    pb <- NULL
  }

  retlist <- .Call("rg_max_concurrent_flow_capacity_restricted_c",
                   demandpaths,
                   vdemands,
                   vcapacity,
                   e,
                   progress,
                   pb,
                   environment()
                   );

  if( progress != FALSE) {
    close(pb)
  }

  ## now need to unpack results
  demands <- savedemands
  
  for(n in 1:length(demands)) {
    flow <- 0
    for(i in 1:length(demands[[n]]$paths)) {
      flow <- flow + retlist$pathflows[[n]][[i]]
      demands[[n]]$paths[[i]] <- retlist$pathflows[[n]][[i]]
    }
    demands[[n]]$flow <- flow
  }

  gdual <- g
  gdual <- rg.set.weight(gdual,retlist$vlengths)
  scalef <- 1 / log(1/delta,base=(1+e))
  demands <- rg.max.concurrent.flow.rescale.demands.flow(demands,scalef)
  beta <- calcBeta(demands,gdual)
  betar <- calcBetaRestricted(demands,gdual)
  lambda=NULL
  lambda <- calcLambda(demands)
  foundratio <- beta / lambda
  ratiobound <- (1-e)^-3
    
  gflow <- rg.max.concurrent.flow.graph(gdual,demands)
  retval <- list(demands=demands,gflow=gflow,gdual=gdual,beta=beta,betar=betar,
                 lambda=lambda,phases=retlist$totalphases,e=e,vlength=retlist$vlengths)

  
  return(retval)
}


### Compute the number of combinations of selecting paths from two path sets
### input
### Given c ways of chosing alternative paths from another path set size m
### there are ( n )
###           ( m )    ways (Binomial coeficient)
### so for m = 0 ... n numbers chosen (choose none, then two etc) we have
### number of comb = Sum_0^n
rg.num.path.comb <- function(n) {
  sum <- 0
  for(m in 0:n)
    sum <- sum + choose(n,m)
  return(sum)
}

rg.dual.primal.ratio <- function(numedges,e,eInternal=NULL)  {
  if(is.null(eInternal))
    eInternal <- e
  
  delta <- (numedges / (1-e)) ^ (-1/e)

  ratio <- e * log(1/delta,1+e) / ( (1-e)*log((1-e)/(numedges*delta)))
  return(ratio)

}

rg.run.reed <- function(scenario) {
  g <- scenario$g
  demands <- scenario$commr
  results <- rg.minimum.congestion.flow(g,demands,e=0.02,progress=TRUE,permutation="lowest")
  res <- rg.max.concurrent.flow.int.c(g,results$demands,e=0.02,progress=TRUE,permutation="lowest")
  cat("Gamma = ",res$bestgamma,"\n")
}

rg.run.sp <- function(scenario) {
  g <- scenario$g
  demands <- scenario$commr
  res <- rg.sp.max.concurrent.flow(g,demands)
  cat("Gamma = ",res$gamma,"\n")
}

rg.run.alice <- function() {
  for(i in 10:20) {
    s <- rg.gen.random.scenario(N=i,L=2*i,target=0.2)
    rg.run.reed(s)
    rg.run.sp(s)
  }

}

rg.max.int.flow <- function(g,demands,e=0.1,progress=FALSE,
                            permutation="fixed",deltaf=1.0) {
  savedemands <- demands
  ## first obtain a reasonable feasible flow using
  ## a poor shortest path flow
  estimate <- rg.sp.max.concurrent.flow(g,demands)
  ## we know that estimate$lambda < Beta
  ## so scale demands by this much to make sure now
  ## that beta > 1
  estlambdascale <- estimate$lambda
  demands <- rg.rescale.demands(demands,estimate$lambda)

  ## Beta might be very high so obtain 2-approximate solution
  ## (1+w) = 2 = (1-e)^-3
  e2= 1- 2^(-1/3)

  if(progress != FALSE)
    progress <- "Calculating 2 opt"
  
  res.2opt <- rg.fleischer.max.concurrent.flow.c(g,demands,e=e2,
                                                 updateflow=FALSE,
                                                 progress=progress,permutation,
                                                 deltaf)
  
  ## now have 2-opt solution as beta_est so  beta < beta_est < 2*beta
  ## scaling by beta_est/2 give best demand for reduced running time
  scalef <- res.2opt$beta /2
  opt2scale <- scalef
  
  demands <- rg.rescale.demands(demands,scalef)


  if(progress != FALSE)
    progress <- "Main calculation"

  res <- rg.fleischer.max.concurrent.flow.with.int.c(g,demands,e=e,
                                                     progress=progress,updateflow=FALSE,
                                                     permutation,
                                                     deltaf)
  
  for(n in names(demands)) {
    res$demands[[n]]$demand <- savedemands[[n]]$demand
  }
  
  res$lambda <- calcLambda(res$demands)
    
  res$beta <- calcBeta(res$demands,res$gdual)
    
  res$estlambdascale <- estlambdascale
  res$opt2scale <- opt2scale
  res
}


rg.fleischer.max.concurrent.flow.with.int.c <- function(g,demands,e=0.1,updateflow=TRUE,progress=FALSE,permutation="lowest",deltaf=1.0) {


  em <- edgeMatrix(g)
  nN <- nodes(g)
  nv <- length(nN)
  
  ne <- ncol(em)
  eW <- unlist(edgeWeights(g))

  cap <- as.double(edgeData(g,attr="capacity"))

  demands.sources <- as.integer(lapply(demands,"[[","source"))
  demands.sinks <- as.integer(lapply(demands,"[[","sink"))
  demands.demand <- lapply(demands,"[[","demand")

  doubreq <- 2/e * log(ne/(1-e),base=(1+e))

  if(progress != FALSE) {
    pb <- txtProgressBar(title = "progress bar", min = 0,
                        max = doubreq, style=3)
  } else {
    pb <- NULL
  }
  if(is.numeric(permutation))
    permutation <- permutation - 1
  else if(permutation == "fixed")
    permutation <- 0: (length(demands) -1)
  else if(permutation == "random")
    permutation <- -1
  else if(permutation == "lowest")
    permutation <- -2
  
  #permutation <- seq(0,length(demands)-1)
  retlist <- .Call("rg_fleischer_max_concurrent_flow_with_int_c",
                   as.integer(nv),
                   as.integer(ne),
                   as.integer(em-1),
                   as.double(eW),
                   as.double(cap),
                   as.integer(length(demands)),
                   as.integer(demands.sources -1 ),
                   as.integer(demands.sinks -1),
                   as.double(demands.demand),
                   as.double(e),
                   as.logical(updateflow),
                   pb,
                   environment(),
                   progress,
                   permutation,
                   deltaf
                 )
  
  
  demflow <- retlist[[2]]
  for(c in names(demands)) {
    demands[[c]]$flow <- demflow[[as.integer(c)]]
    demands[[c]]$paths <- list()
  }


  if(retlist[[3]][[1]] > 0) {
    retdemkey <- retlist[[4]]
    retdemval <- retlist[[5]]
    i <- 1
    for(n in seq(1:retlist[[3]][[1]])) {
      demand <- retdemkey[[i]]
      i <- i+1
      path <- retdemkey[[i]]
      i <- i+1
      demands[[demand]]$paths[[path]] <- retdemval[[n]]
    }
  }
  delta <- deltaf * (ne / (1-e)) ^ (-1/e) 

  scalef <- 1 / log(1/delta,base=(1+e))

  
  gdual <- g

  fromlist <- edgeMatrix(g)[1,]
  tolist <- edgeMatrix(g)[2,]

  edgeData(gdual,from=as.character(fromlist),to=as.character(tolist),attr="weight") <- retlist[[1]]

  demands <- rg.max.concurrent.flow.rescale.demands.flow(demands,scalef)
  gflow <- rg.max.concurrent.flow.graph(gdual,demands)

  beta <- calcBeta(demands,gdual)
  
  lambda=NULL
  if(updateflow) {
    lambda <- calcLambda(demands)
    foundratio <- beta / lambda
    ratiobound <- (1-e)^-3
  }

  if( progress != FALSE) {
    close(pb)
  }

  retlist2 <- list(demands=demands,gflow=gflow,gdual=gdual,beta=beta,lambda=lambda,
                   phases=retlist[[3]][[2]],e=e)
}

### Calculates the number of edge disjoint paths
### between every pair of nodes
### Returns a matrix where element [u,v] is the
### number of paths between u and v where
### 0 means disconnected (or to itself)
### 1 for one path etc
### NOTE this is an R stub calling C++ code
### NOTE THIS IS INCORRECT!
rg.num.edge.disjoint.paths.c <- function(g) {
  em <- edgeMatrix(g)
  nN <- nodes(g)
  nv <- length(nN)

  ne <- ncol(em)
  eW <- unlist(edgeWeights(g))

  count <- .Call("rg_num_edge_disjoint_paths_c",
     as.integer(nv), as.integer(ne),
     as.integer(em-1), as.double(eW))
  count <- matrix(count,nrow=nv,ncol=nv)
  count
}

### Demo for students to simulate a "dice" loaded with exponential
### like distribution.
### input: n - number of throws to return
### rate - the parameter in the (negative expontential)
### returns vector of n throw results
loadeddice <- function(n,rate=0.7) {
  x <- floor(rexp(n,rate))+1
  x <- x[ 1 <= x & x <=6]
  return(x)
}

### Demo for students to simulate a "dice"
### input: n - number of throws to return
### returns vector of n throw results
dice <- function(n) {
  x <- floor(runif(n,min=1,max=7))
  x <- x[ 1 <= x & x <=6]
  return(x)
}

### testing prescaled rg.max.concurrent.flow.int Ran this 20150920 -
### it showed that the prescaled is much faster (mainly for high
### lambda) for very similar quality. Gamma was the same. They both
### had slightly different number of blocked demands when lambda was
### low (lambda<1.0). On average (over 50 results with blocking) the
### ratio of (prescaled-nonscaled)/prescaled was 0.014 (ie nonscaled
### performed very slightly better but there was only 1% in it,
### sometimes nonscaled was better), this was not a statistically
### significant result (p-value 0.08)
rg.max.concurrent.flow.int.test <- function() {
    Nrange <- seq(10,100,10)
    results <<- data.frame(N=numeric(),
                          Lambda=numeric(),
                          Blocked=numeric(),
                          Gamma=numeric(),
                          Runtime=numeric(),
                          Type=character())
    for(N in Nrange) {
        g <- rg.generate.random.graph(N)
        demands <- rg.gen.demands(g,val=round(runif(5,1,4)))
        numDemands <- length(demands)
        g <- rg.set.capacity(g,1)
        linres <- rg.minimum.congestion.flow(g,demands)
        ##linres <- rg.minimum.congestion.flow(g,demands)
        scales <- seq(0.1,2.0,0.2)*linres$lambda
        for(scale in scales) {
            scaled <- rg.rescale.demands(demands,scale)
            timing <- system.time(res <- rg.max.concurrent.flow.int(g,scaled))
            blocked <- numDemands - rg.count.accepted.demands(res$demands)
            interim <- data.frame(N=N,
                                  Lambda=res$lambda,
                                  Gamma=res$gamma,
                                  Blocked=blocked,
                                  Runtime=timing[["user.self"]],
                                  Type="Prescaled")
            print(interim)
            results <<- rbind(results,interim)
            timing <- system.time(res <- rg.max.concurrent.flow.int(g,scaled,
                                                                    prescaled=FALSE))
            blocked <- numDemands - rg.count.accepted.demands(res$demands)
            interim <- data.frame(N=N,
                                  Lambda=res$lambda,
                                  Gamma=res$gamma,
                                  Blocked=blocked,
                                  Runtime=timing[["user.self"]],
                                  Type="NoScaling")
            print(interim)
            results <<- rbind(results,interim)
        }
    }
    return(results)
}
martinjreed/reedgraph documentation built on May 21, 2019, 12:39 p.m.