#########################
#
#  CODE DEVELOPED BY JAMES THORSON FOR PAPER "The importance of spatial models for estimating the strength of density dependence"
#
# NOTES
#  1.  Recommend using Windows 64-bit with R v. 3.0.2
#  2.  TMB is downloaded from GitHub: https://github.com/kaskr/adcomp
#  3.  R-INLA is also downloaded separately: http://www.r-inla.org/download
#  5.  The actual TMB dynamically linked library ("spacetime_v7c.dll" and "time_v2b.dll") must be compiled by each user from included scripts ("spacetime_v7c.cpp" and "time_v2b.cpp"), as done at line 249
#
#########################

#########################
# Simulation evaluation of spatial and non-spatial Gompertz model
#########################

# File structure -- must be modified for each system
File = TmbFile = "C:/Users/James.Thorson/Desktop/UW Hideaway/Collaborations/2014 -- Spatial Gompertz model/Drafts/Ecology/R2/Supplement/"

# Date file
Date = Sys.Date()
  DateFile=paste(File,Date,"/",sep="")
  dir.create(DateFile)

# Random seed
RandomSeed = as.numeric(paste(na.omit(as.numeric(strsplit(as.character(Date),"")[[1]])),collapse=""))
  if( !exists("ThreadNum") ) ThreadNum = 1
  RepSet = 1:100
  RepSet = RepSet + max(RepSet)*(ThreadNum-1) 

# Libraries
library(INLA)
library(TMB)
  newtonOption(smartsearch=TRUE)
library(RandomFields)

# Settings
Scenario = 1  # 1=Fig_2 in main text ; 2=Fig_D1 in Appendix D; 3=Fig_3 in main text;  
  ConfigSet = list( 1:3, 4:7, 8:13)[[Scenario]]    
  RhoTrue_Set = c(0.5,0.5,0.0,  0.5,0.5,0.5,0.5,  0.5,0.5,0.5,0.5,0.5,0.5)
  AlphaTrue_Set = c(1.0,1.0,2.0,  1.0,1.01.0,1.0,  1.0,1.0,1.0,1.0,1.0,1.0)   # Average of null field; Equilibrum field = AlphaTrue/(1-RhoTrue)
  PhiTrue_Set = c(0.0,-2.0,0.0,  0.0,0.0,0.0,0.0,  0.0,0.0,0.0,0.0,0.0,0.0)  # log-ratio of median initial field and median equilibrium field
  n_stations_Set = c(100,100,100,  25,50,100,200,  100,100,100,100,100,100)
  n_years_Set = c(10,10,10,  10,10,10,  4,7,10,15,20,25)
  SD_E = 1   # marginal SD for Epsilon (annual innotation)
  SD_O = 1   # marginal SD for Omega (Gompertz field)
  SpatialScale = 0.25 # parameter input to RFsimulate; Range is derived from this
  GridPrediction = FALSE  # TRUE is much slower
  NegBin = 1  # Sampling distribution during estimation -- 0: Poisson; 1: Negative binomial
  MeshType = c("Minimal","Recommended")[1] # Minimal: Refine=FALSE; Recommended: much slower!
  Version_Spatial = "spacetime_v7d"
  Version_Nonspatial = "time_v2d"

# Save record of settings
RecordList = list("ConfigSet"=ConfigSet, "RhoTrue_Set"=RhoTrue_Set, "AlphaTrue_Set"=AlphaTrue_Set, 
  "PhiTrue_Set"=PhiTrue_Set, "n_stations_Set"=n_stations_Set, "n_years_Set"=n_years_Set, "SD_E"=SD_E, 
  "SD_O"=SD_O, "SpatialScale"=SpatialScale, "GridPrediction"=GridPrediction, 
  "NegBin"=NegBin, "MeshType"=MeshType, "Version_Spatial"=Version_Spatial, "Version_Nonspatial"=Version_Nonspatial)
save( RecordList, file=paste(DateFile,"RecordList.RData",sep=""))
  capture.output( RecordList, file=paste(DateFile,"RecordList.txt",sep=""))
file.copy( from=paste(TmbFile,Version_Spatial,".cpp",sep=""), to=paste(DateFile,Version_Spatial,".cpp",sep=""), overwrite=TRUE)  
  file.copy( from=paste(TmbFile,Version_Nonspatial,".cpp",sep=""), to=paste(DateFile,Version_Nonspatial,".cpp",sep=""), overwrite=TRUE)  

# Loop across all configurations in a given scenario, and all replicates within a configuration
for(ConfigI in 1:length(ConfigSet)){
  for(RepI in RepSet){
    set.seed( RandomSeed+RepI )
    
    # Settings
    RhoTrue = RhoTrue_Set[ConfigSet[ConfigI]]
    AlphaTrue = AlphaTrue_Set[ConfigSet[ConfigI]]
    PhiTrue = PhiTrue_Set[ConfigSet[ConfigI]]
    n_stations = n_stations_Set[ConfigSet[ConfigI]]
    n_years = n_years_Set[ConfigSet[ConfigI]]
    
    # Save file
    ConfigFile = paste(DateFile,"Config=",ConfigSet[ConfigI],"/",sep="")
      dir.create(ConfigFile)
    RepFile = paste(ConfigFile,"Rep=",RepI,"/",sep="")
      dir.create(RepFile)
    
    # Define fields
    model_O <- RMgauss(var=SD_O^2*(1-RhoTrue^2), scale=SpatialScale)
    model_E <- RMgauss(var=SD_E^2, scale=SpatialScale)
  
    # Save record
    Settings = list("AlphaTrue"=AlphaTrue, "PhiTrue"=PhiTrue, "RhoTrue"=RhoTrue, "RangeTrue"=RangeTrue, "SpatialScale"=SpatialScale, "n_stations"=n_stations, "n_years"=n_years, "NegBin"=NegBin, "MeshType"=MeshType, "RandomSeed"=RandomSeed, "SD_E"=SD_E, "SD_O"=SD_O)
      capture.output(Settings, file=paste(ConfigFile,"Settings_Thread=",ThreadNum,".txt",sep=""))
      save(Settings, file=paste(ConfigFile,"Settings_Thread=",ThreadNum,".RData",sep=""))
  
    # Simulated data using fixed-station design                                   
      # Dynamics:  B(t+1) = ro*B(t) + B_null
    # Sampling locations
      loc_stations = cbind( "x"=runif(n_stations, min=0,max=1), "y"=runif(n_stations, min=0,max=1) )
      loc_grid = expand.grid( "x"=seq(0,1,length=21), "y"=seq(0,1,length=21) )
      loc_samp = rbind(loc_stations, loc_grid)
      n_grid = nrow(loc_grid)
      n_samp = nrow(loc_samp)
   
    # Simulate null field
      Omega_samp = RFsimulate(model = model_O, x=loc_samp[,'x'], y=loc_samp[,'y'])@data[,1] + AlphaTrue
      Omega_grid = matrix(Omega_samp[-c(1:n_stations)],21,21)
      Omega_stations = Omega_samp[1:n_stations]
      image(Omega_grid)
      Epsilon_samp = array(NA, dim=c(n_years, n_samp))
      Epsilon_grid = array(NA, dim=c(n_years, 21, 21))
      Epsilon_stations = array(NA, dim=c(n_years, n_stations))
    
    # Simulate density fields
      lnB_samp = array(NA, dim=c(n_years,nrow(loc_samp)))
      lnB_grid = array(NA, dim=c(n_years,dim(Omega_grid)))
      # Simulate B
      Epsilon_samp[1,] = RFsimulate(model_E, x=loc_samp[,'x'], y=loc_samp[,'y'])@data[,1]
      Epsilon_grid[1,,] = Epsilon_samp[1,-c(1:n_stations)]
      Epsilon_stations[1,] = Epsilon_samp[1,1:n_stations]
      lnB_samp[1,] = (Omega_samp)/(1-RhoTrue) + PhiTrue + Epsilon_samp[1,]
      lnB_grid[1,,] = lnB_samp[1,-c(1:n_stations)]
      for(i in 2:n_years){
        Epsilon_samp[i,] = RFsimulate(model_E, x=loc_samp[,'x'], y=loc_samp[,'y'])@data[,1]
        Epsilon_grid[i,,] = Epsilon_samp[i,-c(1:n_stations)]
        Epsilon_stations[i,] = Epsilon_samp[i,1:n_stations]
        lnB_samp[i,] = RhoTrue * lnB_samp[i-1,] + Omega_samp + Epsilon_samp[i,]
        lnB_grid[i,,] = lnB_samp[i,-c(1:n_stations)]
      }
      lnB_stations = lnB_samp[,1:n_stations]
  
      # simulate trawl data
      Y = rpois(n=n_stations*n_years, lambda=exp(lnB_samp[cbind(rep(1:n_years,each=n_stations),rep(1:n_stations,n_years))]) )
      Y_stations = matrix(Y, nrow=n_years, byrow=TRUE)
  
    # Make INLA inputs for TMB
      if(MeshType=="Recommended"){
        mesh_stations = inla.mesh.create( loc_stations, plot.delay=NULL, extend=list(n=8,offset=-0.15), refine=list(min.angle=26,max.edge.data=0.08,max.edge.extra=0.2) )  # loc_samp
        mesh_samp = inla.mesh.create( loc_samp, plot.delay=NULL, extend=list(n=8,offset=-0.15), refine=list(min.angle=26,max.edge.data=0.08,max.edge.extra=0.2) )  # loc_samp
      }
      if(MeshType=="Minimal"){
        mesh_stations = inla.mesh.create( loc_stations, plot.delay=NULL, extend=list(n=8,offset=-0.15), refine=F )  # loc_samp
        mesh_samp = inla.mesh.create( loc_samp, plot.delay=NULL, extend=list(n=8,offset=-0.15), refine=F )  # loc_samp
      }
      spde_stations = inla.spde2.matern(mesh_stations, alpha=2)  # Given 2D field, alpha=2 <=> Matern Nu=1 (Lindgren and Rue 2013, between Eq. 1&2)
      spde_samp = inla.spde2.matern(mesh_samp, alpha=2)
      plot(mesh_stations)
      Y_samp = cbind(Y_stations, matrix(NA, nrow=n_years, ncol=nrow(loc_grid)))

    # Assemble data 
      # TMB 
      NAind_stations = array( as.integer(ifelse(is.na(Y_stations),1,0)), dim=dim(Y_stations))
      NAind_samp = array( as.integer(ifelse(is.na(Y_samp),1,0)), dim=dim(Y_samp))
      X_stations = cbind( rep(1,n_stations*n_years) )
      X_samp = cbind( rep(1,n_samp*n_years) )
    
    # Run TMB
      # Recompile .DLL files from .CPP TMB code
      compile( paste0(Version_Spatial,".cpp") )
      compile( paste0(Version_Nonspatial,".cpp") )
      
      # Run spatial model
      dyn.load( paste0(TmbFile,dynlib(Version_Spatial)) )
      if(GridPrediction==TRUE | RepI==RepSet[1]){
        # Predict to grid and samples
        Data_spatial = list(n_data=n_samp*n_years, n_p=ncol(X_samp), X=X_samp, Y=as.vector(t(Y_samp)), NAind=as.vector(t(NAind_samp)), NegBin=as.integer(NegBin), n_stations=n_samp, meshidxloc=mesh_samp$idx$loc-1, n_years=n_years, G0=spde_samp$param.inla$M0, G1=spde_samp$param.inla$M1, G2=spde_samp$param.inla$M2)
        Parameters_spatial = list(alpha=c(0.0), phi=0.0, log_tau_E=0.0, log_tau_O=0.0, log_kappa=0.0,	rho=0.5, ln_VarInfl=c(0.0,0.0), Epsilon_input=matrix(0,spde_samp$n.spde,n_years), Omega_input=rep(0,spde_samp$n.spde))
      }else{
        # Predict to samples only (not grid)
        Data_spatial = list(n_data=n_stations*n_years, n_p=ncol(X_stations), X=X_stations, Y=as.vector(t(Y_stations)), NAind=as.vector(t(NAind_stations)), NegBin=as.integer(NegBin), n_stations=n_stations, meshidxloc=mesh_stations$idx$loc-1, n_years=n_years, G0=spde_stations$param.inla$M0, G1=spde_stations$param.inla$M1, G2=spde_stations$param.inla$M2)
        Parameters_spatial = list(alpha=c(0.0), phi=0.0, log_tau_E=0.0, log_tau_O=0.0, log_kappa=0.0,	rho=0.5, ln_VarInfl=c(0.0,0.0), Epsilon_input=matrix(0,spde_stations$n.spde,n_years), Omega_input=rep(0,spde_stations$n.spde))
      }
      obj_spatial <- MakeADFun(data=Data_spatial, parameters=Parameters_spatial, random=c("Epsilon_input","Omega_input"), hessian=FALSE)
      # Settings
      obj_spatial$control <- list(trace=1,parscale=rep(1,13),REPORT=1,reltol=1e-12,maxit=100)
      obj_spatial$fn(obj_spatial$par)
      # Run optimizer
      opt_spatial = nlminb(start=obj_spatial$par, objective=obj_spatial$fn, gradient=obj_spatial$gr, lower=c(rep(-20,2),rep(-10,3),-0.999,rep(-10,2)), upper=c(rep(20,2),rep(10,3),0.999,rep(10,2)), control=list(eval.max=1e4, iter.max=1e4))
        Report_spatial = obj_spatial$report()
      if( opt_spatial$par['rho'] > -0.99 & opt_spatial$par['rho'] < 0.99 ){
        sdreport_spatial = try( sdreport(obj_spatial) ) 
        if( !("condition" %in% names(attributes(sdreport_spatial))) ) opt_spatial[["summary"]] = summary(sdreport_spatial)
      }
  
      # Run non-spatial model
      dyn.load( paste(TmbFile,dynlib(Version_Nonspatial),sep="") )
      Data_nonspatial = list(n_data=n_stations*n_years, n_p=ncol(X_stations), X=X_stations, Y=as.vector(t(Y_stations)), NAind=as.vector(t(NAind_stations)), NegBin=as.integer(NegBin), n_stations=n_stations, n_years=n_years)
      Parameters_nonspatial = list(alpha=c(0.0), phi=0.0, log_tau=0.0, rho=0.5, ln_VarInfl=c(0.0,0.0), Epsilon_input=rep(0.0,n_years))
      obj_nonspatial <- MakeADFun(data=Data_nonspatial, parameters=Parameters_nonspatial, random=c("Epsilon_input"))
      # Settings
      obj_nonspatial$control <- list(trace=1,parscale=rep(1,13),REPORT=1,reltol=1e-12,maxit=100)
      obj_nonspatial$hessian <- FALSE
      obj_nonspatial$fn(obj_nonspatial$par)
      # Run optimizer
      opt_nonspatial = nlminb(obj_nonspatial$par, obj_nonspatial$fn, obj_nonspatial$gr, lower=c(rep(-20,2),rep(-10,1),-0.999,rep(-10,2)), upper=c(rep(20,2),rep(10,1),0.999,rep(10,2)), control=list(eval.max=1e4, iter.max=1e4) )
        Report_nonspatial = obj_nonspatial$report()
      if( opt_nonspatial$par['rho'] > -0.99 & opt_nonspatial$par['rho'] < 0.99 ){
        sdreport_nonspatial = try( sdreport(obj_nonspatial) ) 
        if( !("condition" %in% names(attributes(sdreport_nonspatial))) ) opt_nonspatial[["summary"]] = summary(sdreport_nonspatial)
      }
    
    # Save record
    Save = list("Omega_grid"=Omega_grid, "lnB_grid"=lnB_grid, "Epsilon_grid"=Epsilon_grid, "Omega_stations"=Omega_stations, "lnB_stations"=lnB_stations, "Epsilon_stations"=Epsilon_stations) 
      if( exists("opt_spatial") ) Save[["opt_spatial"]] = opt_spatial
      if( exists("opt_nonspatial") ) Save[["opt_nonspatial"]] = opt_nonspatial
      if( exists("Report_spatial") ) Save[["Report_spatial"]] = Report_spatial
      if( exists("Report_nonspatial") ) Save[["Report_nonspatial"]] = Report_nonspatial
      save(Save, file=paste(RepFile,"Save.RData",sep=""))
    
    # Visualize predicted surface on spatial grid (only do if GridPrediction=TRUE or RepI=1, and if model converged)
    if( exists("sdreport_spatial") ){
    if( !("condition" %in% names(attributes(sdreport_spatial))) & (GridPrediction==TRUE | RepI==RepSet[1]) ){
      pow = function(Num1,Num2) Num1^Num2
      REs = sdreport_spatial$par.random		 
      Report = obj_spatial$report()
      # Omega
      Omega_grid_hat = matrix( Report$Omega[mesh_samp$idx$loc][-c(1:n_stations)], ncol=21, nrow=21)
        Equil_grid_hat = matrix( Report$Equil[mesh_samp$idx$loc][-c(1:n_stations)], ncol=21, nrow=21)
      # Epsilon_grid + lnB_grid
      Epsilon_grid_hat = lnB_grid_hat = array(NA, dim=dim(lnB_grid))
      for(YearI in 1:n_years){
        lnB_grid_hat[YearI,,] = matrix( Report$Dji[mesh_samp$idx$loc[-c(1:n_stations)],YearI], ncol=21, nrow=21)
        Epsilon_grid_hat[YearI,,] = matrix( Report$Epsilon[mesh_samp$idx$loc[-c(1:n_stations)],YearI], ncol=21, nrow=21)
      }
      plot(x=lnB_grid_hat, y=lnB_grid)
      abline(a=0,b=1)
      # plot Epsilon_grid_hat
      Nrow = ceiling(sqrt(n_years)); Ncol=ceiling(n_years/Nrow)
      png(file=paste(RepFile,"Epsilon_grid_hat.png",sep=""), width=Ncol*2, height=Nrow*2, res=200, units="in")
        par(mfrow=c(Nrow,Ncol), mar=c(3,3,2,0), mgp=c(2,0.5,0))
        for(i in 1:n_years) image(Epsilon_grid_hat[i,,])
      dev.off()
      # plot Epsilon_grid
      Nrow = ceiling(sqrt(n_years)); Ncol=ceiling(n_years/Nrow)
      png(file=paste(RepFile,"Epsilon_grid.png",sep=""), width=Ncol*2, height=Nrow*2, res=200, units="in")
        par(mfrow=c(Nrow,Ncol), mar=c(3,3,2,0), mgp=c(2,0.5,0))
        for(i in 1:n_years) image(Epsilon_grid[i,,])
      dev.off()
       # plot lnB_grid_hat
      Nrow = ceiling(sqrt(n_years)); Ncol=ceiling(n_years/Nrow)
      png(file=paste(RepFile,"lnB_grid_hat.png",sep=""), width=Ncol*2, height=Nrow*2, res=200, units="in")
        par(mfrow=c(Nrow,Ncol), mar=c(3,3,2,0), mgp=c(2,0.5,0))
        for(i in 1:n_years) image(lnB_grid_hat[i,,])
      dev.off()
      # plot lnB_grid vs. lnB_grid_hat
      png(file=paste(RepFile,"lnB_true_vs_est.png",sep=""), width=5, height=5, res=200, units="in")
        par(mar=c(3,3,2,0), mgp=c(2,0.5,0), tck=-0.02)
        plot( x=lnB_grid_hat, y=lnB_grid )
        abline(a=0,b=1)
      dev.off()
      # Save estimates
      Save_Spatial = list("lnB_grid_hat"=lnB_grid_hat, "Omega_grid_hat"=Omega_grid_hat, "Equil_grid_hat"=Equil_grid_hat)
      save( Save_Spatial, file=paste(RepFile,"Save_Spatial.RData",sep=""))
    }}
    
  } # RepI 
} # ConfigI

#########################################
# Analyze results
#########################################

  DateFile = "C:/Users/James.Thorson/Desktop/UW Hideaway/Collaborations/2014 -- Spatial Gompertz model/2014-07-08-1/"
  RepSet = 1:100
  ConfigSet = (list( 1:3, 4:7, 8:13 )[[1]])
  ConfigNames = list( 1:3, c(25,50,100,200), c(4,7,10,13), c(4,7,10,15,20,25) )[[1]]
  ParamResults = array(NA, dim=c(2,length(ConfigSet),length(RepSet),7,5), dimnames=list(c("Spatial","Nonspatial"),paste("Config=",ConfigSet,sep=""),paste("Rep=",RepSet,sep=""),c("Rho","Phi","Alpha","meanK","SigmaE","SigmaO","Change_Error"),c("10%","50%","90%","Mean","SD")))

  # Load previous runs
  Change = function(Vec){ rev(Vec)[1]/Vec[1] }
  for(ConfigI in 1:length(ConfigSet)){
    ConfigFile = paste(DateFile,"Config=",ConfigSet[ConfigI],"/",sep="")
    load( paste(ConfigFile,"Settings_Thread=1.RData",sep="") )
    for(LoopI in 1:length(RepSet)){
      RepI = RepSet[LoopI]
      RepFile = paste(ConfigFile,"Rep=",RepI,"/",sep="")
      if( "Save.RData" %in% list.files(RepFile) ){
        load( paste(RepFile,"Save.RData",sep="") )
        if( !is.null(Save$opt_spatial$summary) ){
          ParamResults['Spatial',ConfigI,LoopI,c("Rho","Phi","Alpha","SigmaE","SigmaO"),c("10%","50%","90%")] = Save$opt_spatial$summary[c('rho','phi','alpha','SigmaO','SigmaE'),'Estimate']%o%rep(1,3) + Save$opt_spatial$summary[c('rho','phi','alpha','SigmaO','SigmaE'),'Std. Error']%o%qnorm(c(0.1,0.5,0.9))
          ParamResults['Spatial',ConfigI,LoopI,c("Rho","Phi","Alpha","SigmaE","SigmaO"),c("Mean","SD")] = cbind( Save$opt_spatial$summary[c('rho','phi','alpha','SigmaO','SigmaE'),'Estimate'], Save$opt_spatial$summary[c('rho','phi','alpha','SigmaO','SigmaE'),'Std. Error'] )
          ParamResults['Spatial',ConfigI,LoopI,c("Change_Error"),c("10%","50%","90%","Mean","SD")] = (Save$opt_spatial$summary[c('Change'),'Estimate']-Change(apply( exp(Save$lnB_grid),MARGIN=1,FUN=mean)))%o%c(1,1,1,1,0) + Save$opt_spatial$summary[c('Change'),'Std. Error']%o%c(qnorm(c(0.1,0.5,0.9)),0,1)
        }
        if( !is.null(Save$opt_nonspatial$summary) ){
          ParamResults['Nonspatial',ConfigI,LoopI,c("Rho","Phi","Alpha"),c("10%","50%","90%")] = Save$opt_nonspatial$summary[c('rho','phi','alpha'),'Estimate']%o%rep(1,3) + Save$opt_nonspatial$summary[c('rho','phi','alpha'),'Std. Error']%o%qnorm(c(0.1,0.5,0.9))
          ParamResults['Nonspatial',ConfigI,LoopI,c("Rho","Phi","Alpha"),c("Mean","SD")] = cbind( Save$opt_nonspatial$summary[c('rho','phi','alpha'),'Estimate'], Save$opt_nonspatial$summary[c('rho','phi','alpha'),'Std. Error'] )
          ParamResults['Nonspatial',ConfigI,LoopI,c("Change_Error"),c("10%","50%","90%","Mean","SD")] = (Save$opt_nonspatial$summary[c('Change'),'Estimate']-Change(apply( exp(Save$lnB_grid),MARGIN=1,FUN=mean)))%o%c(1,1,1,1,0) + Save$opt_nonspatial$summary[c('Change'),'Std. Error']%o%c(qnorm(c(0.1,0.5,0.9)),0,1)
        }
      }
    }
  }
  
  # Fig. 3/4 -- Plot posteriors for estimates -- boxplots
  BoxFn = function(Mat, At, Col, Xlim, Ylim, New, Lwd, ...){
    if(New==TRUE) plot( 1, type="n", xlim=Xlim, ylim=Ylim, ...)
    for(i in 1:ncol(Mat)){
      segments( x0=At[i], x1=At[i], y0=quantile(Mat[,i],prob=0.025,na.rm=TRUE), y1=quantile(Mat[,i],prob=0.25,na.rm=TRUE), lwd=Lwd )
      segments( x0=At[i], x1=At[i], y0=quantile(Mat[,i],prob=0.75,na.rm=TRUE), y1=quantile(Mat[,i],prob=0.975,na.rm=TRUE), lwd=Lwd )
      polygon( x=At[i]+c(-0.5,0.5)[c(1,2,2,1,1)], y=quantile(Mat[,i],prob=c(0.25,0.75),na.rm=TRUE)[c(1,1,2,2,1)], col=Col, lwd=Lwd)
      segments( x0=At[i]-0.5, x1=At[i]+0.5, y0=median(Mat[,i],na.rm=TRUE), y1=median(Mat[,i],na.rm=TRUE), lwd=Lwd)
    }
  }
  cMx=function(Input){as.matrix(Input)}
  rMx=function(Input){
    if(is.vector(Input)){Output<-t(as.matrix(Input))}
    if(!is.vector(Input)){Output<-as.matrix(Input)}
    Output
  }
  for(resI in 1:2){
    png( file=paste(DateFile,"Param_boxplots_",c("LO","HI")[resI],".png",sep=""), width=2*2.5, height=2*2.5, res=c(200,1200)[resI], units="in")
      par( mfrow=c(2,2), mar=c(0,1.5,0,0), mgp=c(1.5,0.25,0), tck=-0.02, oma=c(3,2,2,0), yaxs="i", xaxs="i")
      for(i in 1:2){
        Ylim1 = list( c(-1,1), c(-4,4), c(0,3), c(0,3), c(-1,1) )[[i]]
        Ylim2 = list( c(0,1), c(0,3), c(0,1), c(0,1), c(0,1) )[[i]]
        At = (1:length(ConfigSet)-1)*2 + 1
        # Plot means
        BoxFn( t(rMx(ParamResults['Spatial',,,c('Rho','Phi','SigmaO','SigmaE','Change_Error')[i],'Mean'])), Xlim=range(At)+c(-1,2), Ylim=Ylim1, Col=rgb(0,0,0,0.3), At=At, New=TRUE, xaxt="n", Lwd=1.5, xlab="", ylab="", main="", yaxt="n" )
        BoxFn( t(rMx(ParamResults['Nonspatial',,,c('Rho','Phi','SigmaO','SigmaE','Change_Error')[i],'Mean'])), Xlim=range(At)+c(-1,2), Ylim=Ylim1, Col=rgb(1,1,1,0.3), At=At+1, New=FALSE, Lwd=1.5 )
        for(ConfigI in 1:length(ConfigSet)){
          ConfigFile = paste(DateFile,"Config=",ConfigSet[ConfigI],"/",sep="")
          load( paste(ConfigFile,"Settings_Thread=1.RData",sep="") )
          if(any(i==1:5)){
            True = c(Settings[c("RhoTrue","PhiTrue","SD_O","SD_E")],0)[[i]]
            segments( x0=(ConfigI-1)*2+0.5, x1=(ConfigI-1)*2+2.5, y0=True, y1=True, lwd=3, lty="dotted")
          }
        }
        if(i==2){
          #axis( side=1, at=0, labels="Scenario=", tick=FALSE )
          axis( side=1, at=At+0.5, labels=ConfigNames, tick=FALSE )
          axis( side=1, at=At+0.5, labels=rep("",length(ConfigSet)*1), tick=TRUE )
          if(all(c(1,4,5,6)%in%ConfigSet)) axis( side=1, at=At[1]-1, labels=expression(paste(n[stations]," =",sep="")), tick=FALSE)
          if(all(c(1,8,9)%in%ConfigSet)) axis( side=1, at=At[1]-1, labels=expression(paste(n[years]," =",sep="")), tick=FALSE)
        }
        if(i%in%c(1,3,4,5)) axis( side=2)
        if(i==2) axis( side=2, at=rev(rev(axTicks(2))[-1]), labels=rev(rev(axTicks(2))[-1]))
        if(i%in%c(1:5)) mtext(side=2, outer=FALSE, line=2, text=list(expression(italic(rho)),expression(italic(phi)),"SigmaO","SigmaE",expression(Delta))[[i]])
        if(i==1) mtext(side=3, outer=FALSE, line=0.5, text="Estimated value")
        # Plot SDs
        BoxFn( t(rMx(ParamResults['Spatial',,,c('Rho','Phi','SigmaO','SigmaE','Change_Error')[i],'SD'])), Xlim=range(At)+c(-1,2), Ylim=Ylim2, Col=rgb(0,0,0,0.3), At=At, New=TRUE, xaxt="n", Lwd=1.5, xlab="", ylab="", main="", yaxt="n" )
        BoxFn( t(rMx(ParamResults['Nonspatial',,,c('Rho','Phi','SigmaO','SigmaE','Change_Error')[i],'SD'])), Xlim=range(At)+c(-1,2), Ylim=Ylim2, Col=rgb(1,1,1,0.3), At=At+1, New=FALSE, Lwd=1.5 )
        if(i==2){
          axis( side=1, at=At+0.5, labels=ConfigNames, tick=FALSE )
          axis( side=1, at=At+0.5, labels=rep("",length(ConfigSet)*1), tick=TRUE )
        }
        if(i%in%c(1,3,4,5)) axis( side=2)
        if(i==2) axis( side=2, at=rev(rev(axTicks(2))[-1]), labels=rev(rev(axTicks(2))[-1]))
        if(i==1) mtext(side=3, outer=FALSE, line=0.5, text="Standard error")
        mtext(side=1, outer=TRUE, line=1.25, text="Scenario")
      }
    dev.off()
  }
  
  ### Example plots for spatial fields
  Col_spatial = colorRampPalette(colors=c("white","grey","black"))
  Col_true = colorRampPalette(colors=c("white","grey","black"))
  ConfigI = 1
  for(ConfigI in 1:length(ConfigSet)){
    # Fig. 2 combined 
    for(resI in 1:2){
      png( file=paste(ConfigFile,"Fig_2comb_Example_fields_",c("LO","HI")[resI],".png",sep=""), width=2*1.5, height=1.5*(length(YearSet)+1), res=c(200,1200)[resI], units="in")
        par( mfrow=c(length(YearSet)+1,2), mar=c(0,0,0,0), mgp=c(1.5,0.25,0), tck=-0.02, oma=c(2,2,2,1.5)+0.1)
        RepFile = paste(ConfigFile,"/Rep=",RepSet[1],"/",sep="")
        load( paste(ConfigFile,"Settings_Thread=1.RData",sep="") )
        load( paste(RepFile,"Save.RData",sep=""))
        load( paste(RepFile,"Save_Spatial.RData",sep=""))
        # Productivity
        Omega_grid_trans = Save$Omega_grid - Settings[['AlphaTrue']]
        #Omega_grid_hat_trans = Save_Spatial$Omega_grid_hat * Save$opt_spatial
        Zlim = range( range(Omega_grid_trans), range(Save_Spatial$Omega_grid_hat) )
        # True
        image(Omega_grid_trans, main="", xaxt="n", yaxt="n", zlim=Zlim, col=Col_true(11))
        mtext(side=3, text="True", line=0.5)
        # Label
        mtext(side=2, text=expression(bold(Omega)), line=0.5)
        # Est
        image(Save_Spatial$Omega_grid_hat, main="", xaxt="n", yaxt="n", zlim=Zlim, col=Col_spatial(11))
        mtext(side=3, text="Estimated", line=0.5)
        # Biomass 
        Zlim = range( range(Save$lnB_grid), range(Save_Spatial$lnB_grid_hat) )
        for(i in YearSet){
          # True
          image(Save$lnB_grid[i,,], main="", xaxt="n", yaxt="n", zlim=Zlim, col=Col_true(11))
          # Label
          if(i==1) mtext(side=2, text=expression(paste("log(",D[1],")",sep="")), line=0.5)   # paste("t =",i)
          if(i==2) mtext(side=2, text=expression(paste("log(",D[2],")",sep="")), line=0.5)   # paste("t =",i)
          if(i==3) mtext(side=2, text=expression(paste("log(",D[3],")",sep="")), line=0.5)   # paste("t =",i)
          if(i==10) mtext(side=2, text=expression(paste("log(",D[10],")",sep="")), line=0.5)   # paste("t =",i)
          # Est
          image(Save_Spatial$lnB_grid_hat[i,,], main="", xaxt="n", yaxt="n", zlim=Zlim, col=Col_spatial(11))
        }
        #mtext(side=3, outer=TRUE, line=0, text=expression(paste("True                                    Est", cex=2)
        mtext(side=1, outer=TRUE, line=0.5, text="Eastings")
        mtext(side=4, outer=TRUE, line=0.5, text="Northings")      
      dev.off()
    }
  }
}