"sim.conn.GRG" <-
function(Eall,N,prop,rconn,beta,cc.AB,cc.A=cc.AB,cc.B=cc.AB,tol.interp=c(.01,.05,.1,.15,.2,.25),tol.simul=tol.interp,tries=10,skipSP=FALSE,verbose=FALSE)
  {
    ### Simulate connectivity according to input parameters, using geometric random graph
    
    nA <- as.integer(prop*N)
    nB <- N-nA

    denom <- beta*nA^2*rconn - beta*nA*rconn + beta*nA*nB + rconn*nA*nB - rconn*nB + nB^2
    pA <- rconn*beta*N*Eall/denom
    pB <- -N*Eall*(rconn-nA*rconn-nB+beta*nA)/(nB-1)/denom
    pAB <- beta*N*Eall/denom

    MBp <- new("MBparams")
    MBp <- addParams(x=MBp,N=N,m=2)
    MBp <- calcSim(x=MBp,func=sim.thresh,pNum=2)

    # Calculate GRG params for sets A, B, and AB
    if(verbose) cat("Calculating GRG parameters based on geom.interp matrices...\n")
    GRG.pars <- matrix(nrow=3,ncol=2)
    p <- c(pA,pB,pAB)
    cc <- c(cc.A,cc.B,cc.AB)
    
    p.min <- min(geom.interp.P,na.rm=TRUE)
    p.max <- max(geom.interp.P,na.rm=TRUE)
    cc.min <- min(geom.interp.CC,na.rm=TRUE)
    cc.max <- max(geom.interp.CC,na.rm=TRUE)

    row <- 0
    for(i in 1:length(p)) {
      row <- row+1
      if(p[i] < p.min) p[i] <- p.min
      else if(p[i] > p.max) p[i] <- p.max
      if(cc[i] < cc.min) cc[i] <- cc.min
      else if(cc[i] > cc.max) cc[i] <- cc.max
      for(tol.cc in tol.interp) {
        for(tol.p in tol.interp) {
          inds.mat <- which(abs(geom.interp.P-p[i]) <= (p[i]*tol.p) & abs(geom.interp.CC-cc[i]) <= (cc[i]*tol.cc), arr.ind=TRUE)
          if(!is.null(nrow(inds.mat))) {
            GRG.pars[row,] <- inds.mat[sample(nrow(inds.mat),1),]
            GRG.pars[row,1] <- as.numeric(rownames(geom.interp.P)[GRG.pars[row,1]])
            GRG.pars[row,2] <- as.numeric(colnames(geom.interp.P)[GRG.pars[row,2]])
            break
          }
        }
        if(!any(is.na(GRG.pars[row,]))) break
      }
    }
    
    if(any(is.na(GRG.pars))) {
      stop("No adequate GRG parameter values found in geom.interp matrices!")
    }
    if(verbose) cat("Success!\n")

    # Simulate connectivity until results are adequate
    done <- FALSE
    for(tol.sim in tol.simul) {
      for(try in 1:tries) {
        if(verbose && (try==1 || !(try %% 10))) cat(paste("Simulating connectivity at tolerance ",tol.sim,", try ",try,"...\n",sep=""))
        GRG.AB <- as.matrix(dist(cbind(runif(N),runif(N),runif(N))))

        GRG.A <- GRG.AB[1:nA,1:nA]
        qA.1 <- quantile(as.numeric(GRG.A),GRG.pars[1,1])
        GRG.A[GRG.A >= qA.1] <- 0
        GRG.A[GRG.A < qA.1 & GRG.A > 0] <- 1
        rsums.A <- rowSums(GRG.A)
        qA.2 <- quantile(rsums.A,GRG.pars[1,2])
        zero.A <- which(rsums.A <= qA.2)
        GRG.A[zero.A,] <- 0
        GRG.A[,zero.A] <- 0
	if(sum(GRG.A)==0) next

        GRG.B <- GRG.AB[(nA+1):N,(nA+1):N]
        qB.1 <- quantile(as.numeric(GRG.B),GRG.pars[2,1])
        GRG.B[GRG.B >= qB.1] <- 0
        GRG.B[GRG.B < qB.1 & GRG.B > 0] <- 1
        rsums.B <- rowSums(GRG.B)
        qB.2 <- quantile(rsums.B,GRG.pars[2,2])
        zero.B <- which(rsums.B <= qB.2)
        GRG.B[zero.B,] <- 0
        GRG.B[,zero.B] <- 0
	if(sum(GRG.B)==0) next

        qAB.1 <- quantile(as.numeric(GRG.AB),GRG.pars[3,1])
        GRG.AB[GRG.AB >= qAB.1] <- 0
        GRG.AB[GRG.AB < qAB.1 & GRG.AB > 0] <- 1
        rsums.AB <- rowSums(GRG.AB)
        qAB.2 <- quantile(rsums.AB,GRG.pars[3,2])
        zero.AB <- which(rsums.AB <= qAB.2)
        GRG.AB[zero.AB,] <- 0
        GRG.AB[,zero.AB] <- 0

        GRG.AB[1:nA,1:nA] <- GRG.A
        GRG.AB[(nA+1):N,(nA+1):N] <- GRG.B

        gpars <- calc.graph.pars.old(as.matrix.csr(GRG.AB),as.character(1:N),1:N,nA,skipSP=TRUE,verbose=FALSE)
	if(is.na(gpars$mcoefALL) || is.na(gpars$mcoefA) || is.na(gpars$mcoefB)) next
        if(abs(gpars$pAB-pAB) <= (pAB*tol.sim) && abs(gpars$rconn-rconn) <= (rconn*tol.sim) && abs(gpars$EA/gpars$EB-beta) <= (beta*tol.sim) &&
           abs(gpars$mcoefALL-cc.AB) <= (cc.AB*tol.sim) && abs(gpars$mcoefA-cc.A) <= (cc.A*tol.sim) && abs(gpars$mcoefB-cc.B) <= (cc.B*tol.sim)) {
	  gpars <- calc.graph.pars.old(as.matrix.csr(GRG.AB),as.character(1:N),1:N,nA,skipSP=skipSP,verbose=verbose)
          done <- TRUE
          break
        }
      }
      if(done) break
    }
    
    # Return results
    if(!done) {
      stop("Adequate simulated connectivity not obtained!")
    }
    if(verbose) cat("Success!\n")
    MBp@Sim[[1]] <- as.matrix.csr(GRG.AB)
    return(MBp)
  }