Model parameters in Forest Fire Game (igraph)

library(igraph)
library(compiler)

ffg.optim<-function(g, sample.size=100, ambassadors=1, multicore=TRUE, init.fw=.5, init.bw=1, hess=TRUE, ...){
  #_g_ is an igraph 'graph' object.
    #If _g_ has multiple edges and/or loops, they will be removed.
  #_sample.size_ is how many simulations to perform for each optimization iteration
    #Larger values increase the accuracy, yet decrease the speed.
  #_ambassadors_ is the ambassadors parameter for forest.fire.game().
    #It must be an integer greater than or equal to 1.
    #optim has a hard time estimating it.
  #Set _multicore_ to TRUE if you have this package installed and would like a speed boost.
  #_init.fw_ and _init.bw_ are the initial guesses for the forward and backward burning parameters.
  #_hess_ is to set the hessian logical argument in optim. Defaults to TRUE.
  #_..._ are other arguments to pass to optim.  E.g., hessian=TRUE will provide a hessian matrix for the estimate.
  require(igraph)
  if (is.simple(g)==FALSE)
    g<-simplify(g)
  ambassadors<-round(ambassadors)
  if (ambassadors<1){
    ambassadors<-1
    cat("\n  WARNING: ambassadors must be greater than or equal to one. Setting ambassadors=1 to correct user error.\n")
    }
  obs.stats<-function(the.g){
    dens<-graph.density(the.g)
    apl<-average.path.length(the.g, directed=FALSE)
    gtrans<-transitivity(the.g)
    return(c(density=dens, average.path.length=apl, transitivity=gtrans))
    }
  serr<-function(test.vec, obs) #Squared error function
    return(((test.vec-obs)/obs)^2)
  obs<-obs.stats(g) #Target observations

  min.this<-function(test.mat){ #Sum of squared errors
    a<-sapply(1:length(obs), function(x) return(serr(test.mat[,x], obs[x])))
    return(sqrt(sum(a)))}
  n<-vcount(g)

  fun2optimize<-function(params){
    fw<-params["fw"]; bw<-params["bw"]
    test.vals.wrap<-function(){ #This wrapper saves memory, as the simulations are never stored.
      temp.g<-forest.fire.game(nodes=n, fw.prob=fw, bw.factor=bw, ambs=ambassadors, directed=is.directed(g))
      temp.g.obs<-obs.stats(temp.g)
      return(temp.g.obs)
      }
    if (multicore&("multicore"%in%installed.packages())){
      require(multicore)
      test.vals<-do.call(rbind, mclapply(1:sample.size, function(x) return(test.vals.wrap())))
      }
    else
      test.vals<-do.call(rbind, lapply(1:sample.size, function(x) return(test.vals.wrap())))
    return(min.this(test.vals))
    }
  out.optim<-optim(par=c(fw=init.fw, bw=init.bw), fn=fun2optimize, method="L-BFGS-B", lower=c(fw=.000000001, bw=.000000001), upper=c(fw=1, bw=Inf), hessian=hess, ...)
  return(out.optim)
  }; ffg.optim<-cmpfun(ffg.optim)
g<-forest.fire.game(100, .1, 3, ambs=1, directed=TRUE)
g.fit<-ffg.optim(g, 5000, init.fw=.3, init.bw=1, multicore=TRUE, control=list(maxit=10000, lmm=10))
#g.fit$par is usually quite off and the standard errors, sqrt(diag(solve(-g.fit$hessian))), produce NaNs due to non-invertibility of hessian matrix (positive values)
#The target stats are usually good, though.