#########################################################
#####           0. FUNCTIONS AND LIBRARIES          #####
#########################################################

# Libraries
require(mclust)
require(KernSmooth)
require(MASS)
require(R.utils)


# FUNCTION: 2D histogram of environmental space for numerical availability
Espace.numer<-function(maxF, maxT,steps, env1,env2)
{
  steps<-steps-1
  tempR<-c(0,maxT)
  foodR<-c(0,maxF)
  flab<-seq(foodR[1], foodR[2], len=steps)
  tlab<-seq(tempR[1], tempR[2], len=steps)
  fd<-rep(flab, steps)
  tp<-rep(tlab, 1, each=steps)
  avR<-rep(0, steps^2)
  for(jj in 1:d^2)
  {
    eu<-(fd-env1[jj])^2+(tp-env2[jj])^2
    eumin<-which(eu==min(eu))
    avR[eumin]<-avR[eumin]+1
  }
  avR<-matrix(avR/sum(avR), steps, steps)
  image(log(avR+0.001), xlab="Food", ylab="Temperature", main="e. Numerical availability", axes=FALSE, col=topo.colors(100))
  pos<-(1:steps)/steps
  axis(1, at = pos, labels=round(flab), tick=FALSE)
  axis(2, at = pos, labels=round(tlab), tick=FALSE)
}

# FUNCTION: 2D histogram of environmental space for approximated availability
Espace.approx<-function(maxF, maxT,steps, mus, sit, ws)
{
  steps<-steps-1
  L<-length(ws)
  tempR<-c(0,maxT)
  foodR<-c(0,maxF)
  flab<-seq(foodR[1], foodR[2], len=steps)
  tlab<-seq(tempR[1], tempR[2], len=steps)
  fd<-rep(flab, steps)
  tp<-rep(tlab, 1, each=steps)
  muf<-mus[1,]
  mut<-mus[2,]
  sif<-var[1]
  sit<-var[2]
  av<-exp(-1/2*((fd-muf[1])^2/sif+(tp-mut[1])^2/sit))
  av<-av/sum(av)*ws[1]
  for( l in 2:L)
  {
    avt<-exp(-1/2*((fd-muf[l])^2/sif+(tp-mut[l])^2/sit))
    av<-av+avt/sum(avt)*ws[l]
  }
  av<-matrix(av,steps,steps)
  image(log(av+0.001), xlab="Food", ylab="Temperature", main="f. Approximated availability", axes=FALSE, col=topo.colors(100))
  pos<-(1:steps)/steps
  axis(1, at = pos, labels=round(flab), tick=FALSE)
  axis(2, at = pos, labels=round(tlab), tick=FALSE)
}


# FUNCTION: Generates a random environmental layer in a square dxd arena using a total of x focal points
environ<-function(d,x,mx,bw,sl)
{
  ar<-array(0, dim=c(d,d))
  # Places seeds in arena
  slope<- sl # 0.05 # Introduces SW-to-NE gradient in resource
  cox<-cbind(1+(d-1)*runif(x, min=0, max=1)^slope, 1+(d-1)*runif(x, min=0, max=1)^slope)
  # Smooths seeds to create spatial autocorrelation
  sarx<-bkde2D(cox, bandwidth = c(bw,bw), gridsize=c(d,d),range.x=list(c(1,d),c(1,d)))
  ii<-rep(1:d,d)
  jj<-rep(1:d, each=d)
  
  sarx$fhat<-mx*(sarx$fhat/mean(sarx$fhat))
  return(sarx$fhat)
}

# FUNCTION: Prints out the confidence intervals from a distribution of residuals
resCI<-function(truth, est)
{res<-sort(truth-est)
 hi<-res[round(0.975*length(res))]
 lo<-res[round(0.025*length(res))]
 print(paste("(",lo,",",hi,")"))
}

########################################################
#####           1. POPULATION SIMULATION           #####
########################################################

# GLOBAL PARAMETERS & DATA STRUCTURES

# Environment-related
d<-50  #  Arena dimensions

# Population-related
ymax<-31 # Simulation years
yburnin<-20 #Initialisation years
tmax<-12 # "Months" in the year
p0<-100  # Initial population size
popmin<-5 # Extinction threshold
scmax<-20+100 # Number of environmental scenarios
scfit<-20 # Number of scenarios to be used for training

# Fitness-related
a0<--80 # Intercept (must be < a2^2/(4*a3))
a1<-0.024 # Coefficient for food
a2<-4.8 # First order temperature coefficient
a3<--0.08 # Second order temperature coefficient (must be < 0)
attr<-10   # Coefficient of attrition between conspecifics sharing a grid cell
ss<-10 # Parameter determining uncertainty in detection of fitness
if(a0>a2^2/(4*a3)) print("Check the fitness coefficients")
bmax<- 1
b0<- -5.5
b1<-0.0032
smax<-(0.95)^(1/tmax) #Maximum "monthly" probability of survival - derived from a specified annual background mortality
s0<- -1
s1<-0.02
enI<-500

# Sampling-related
sampInt<-1
dat2<-c()
sc<-1

# SCENARIO LOOP
dat.Master<-data.frame("scenario"=c(), "food"=c(), "temp"=c(), "varf"=c(), "vart"=c(), "weights"=c() )
ydat<-data.frame("use"=c(), "food"=c(), "temp"=c(), "fit"=c(), "Efd"=c(), "Ecv"=c(), "Efd2"=c(), "Ecv2"=c(), "Epr"=c(), "offs"=c(), "N"=c())
dat.N<-data.frame("scenario"=c(), "N2"=c(), "N1"=c(), "AvFood"=c(), "AvTemp"=c(), "AvFood2"=c(), "AvTemp2"=c(), "Avfit"=c())

while(sc<=scmax)
{
  recorded<-FALSE
  burnin<-0
  # Scenario initialisation
  
  x1<-round(runif(1,100,100))  # Number of food distribution foci
  x2<-round(runif(1,5,10))   # Number of temperature distribution foci
  mx1<-runif(1,2500,2600)  # Average amount of food across arena
  mx2<-runif(1,26,27)   # Average value of temperature across arena
  bw1<-3+rexp(1,1/4) # Smoothing bandwidth for generation of environmental layers
  sl1<-runif(1,0.5,2) # Geographical gradient in distribution of food
  sl2<-runif(1,0.5,2) # Geographical gradient in distribution of temperature
  env1<-environ(d,x1,mx1,bw1,sl1)  # Generation of food layer
  env2<-environ(d,x1, mx2,bw1,sl2)  # Generation of temperature layer
  image(env1)
  image(env2)
  fitinit<-a0+a1*env1+a2*env2+a3*env2^2  # True map of fitness in unpopulated lanscape
  
  # Approximation of E-space densities
  edat<-cbind("food"=c(env1),"temp"=c(env2))  # Environmental dataframe
  eclust<-Mclust(edat, G=300, modelNames=c("EEI"))  # Application of Gaussian clustering with 300 kernels
  mus<-eclust$parameters$mean  # Extraction of estimated kernel means
  vars<-eclust$parameters$variance  # Extraction of estimated kernel variances
  var<-c(vars[5][[1]][1],vars[5][[1]][4])  # List of variances in the two environmental dimensions
  ws<-eclust$parameters$pro  # Extraction of estimated kernel weights
  
  print(paste("Now calculating scenario ", sc, "  With average food: ",
              round(mean(env1),5),"(max ",round(max(env1),2),") and temperature: ",
              round(mean(env2),5),"(max ",round(max(env2),2),")"))
  
  print(paste("Approximate expected carrying capacity:", sum(pmax(0,fitinit)/attr)))
  
  # POPULATION MODEL_____________________________________
  par(mfrow=c(4,5))
  y<-2 # Set year
  popc<-p0 # Set current population size to initial value
  pop<-rep(0, ymax) # Initialise population vector with single initial value
  pop[1]<-p0
  x<-matrix(round(runif(2*p0, 1, d)), byrow=T, nrow=p0, ncol=2)   # Initial positions for initial population
  en<-rnorm(p0, enI, 0.1*enI) # Initial condition for animals in population
  traj<-env1*0   # Tracks usage distribution 
  for(i in 1:p0) {traj[x[i,1],x[i,2]]<-traj[x[i,1],x[i,2]]+1} # Initialises usage distribution
  
  while (y<=ymax && popc>popmin)  # Year loop
  {
    print(paste(y-1,") P:",pop[y-1]))
    
    t<-0
    while (t<tmax && popc>popmin)  # "Month" loop
    {
      t<-t+1
      fit<-a0+a1*env1+a2*env2+a3*env2^2-traj*attr   #  Calculates true map of profitability
      dead<-c()
      for (i in 1:length(x[,1]))    #Individuals loop
      {
        fit<-a0+a1*env1+a2*env2+a3*env2^2-traj*attr   #  Calculates true map of profitability
        
        # Redistribution 
        xi<-x[i,1]
        yi<-x[i,2]
        fit[xi,yi]<-fit[xi,yi]+attr
        xu<-ifelse(xi==d, 1, xi+1)
        xl<-ifelse(xi==1, d, xi-1)
        yu<-ifelse(yi==d, 1, yi+1)
        yl<-ifelse(yi==1, d, yi-1)
        nx<-c(xi,xi, xu, xl, xi)
        ny<-c(yi,yu, yi, yi, yl)
        fits<-fit[cbind(nx,ny)]+rnorm(5, 0, ss)
        if(sum(exp(fits))>0)  pem<-exp(fits)/sum(exp(fits)) else pem<-rep(0.2,5) # Probability of moving to different cells in neighborhood
        # Chooses cell to move to
        vs<-rmultinom(1,1,pem) 
        jump<-which(vs==max(vs))
        # Updates position and energy for animal
        en[i]<-min(3000,en[i]+fit[xi,yi])
        traj[xi,yi]<-traj[xi,yi]-1
        x[i,1]<-xi<-nx[jump]
        x[i,2]<-yi<-ny[jump]
        traj[xi,yi]<-traj[xi,yi]+1
        
        # Mortality
        sur<-smax*exp(s0+s1*en[i])/(1+exp(s0+s1*en[i]))  # Calculates survival probability based on current energetic state
        alive<-rbinom(1,1,sur) # Decides if animal will live
        if(alive==0) dead<-c(dead, i)
      }
      
      if(length(dead)>1) 
      {
        x<-x[-dead,]
        en<-en[-dead]
      }
      
      if(recorded==FALSE)
      {
        x<-x[1:min(p0,length(en)),]
        en<-en[1:min(p0,length(en))]
      }
      popc<-length(x[,1])
      traj<-env1*0   # Current usage distribution 
      for(i in 1:popc) {traj[x[i,1],x[i,2]]<-traj[x[i,1],x[i,2]]+1} # Initialises usage distribution
    }
    
    pop[y]<-popc # Records population time series from pre-breeding time
    
    if(popc>1)
    {
      trajPre<-env1*0   # Pre-breeding usage distribution 
      for(i in 1:popc) {trajPre[x[i,1],x[i,2]]<-trajPre[x[i,1],x[i,2]]+1} # Initialises usage distribution
      image(trajPre, zlim=c(0,5))
      # Fecundity
      for (i in 1:popc)    #Individuals loop
      {
        b<-bmax*exp(b0+b1*en[i]) #/(1+exp(b0+b1*en[i]))
        if(b>0)
        {      
          births<-rpois(1,b)
          en[i]<-en[i]-births*enI
          popc<-popc+births
          x<-rbind(x, cbind(round(runif(births,max(1,x[i,1]-10),min(d,x[i,1]+10))),round(runif(births,max(1,x[i,2]-10),min(d,x[i,2]+10))))) # New births positioned near parent's location
          en<-c(en, rnorm(births, enI, 0.1*enI) )
        }
      }
      traj<-env1*0   # Current usage distribution 
      for(i in 1:popc) {traj[x[i,1],x[i,2]]<-traj[x[i,1],x[i,2]]+1} # Initialises usage distribution
      
    }
    
    
    
    if(burnin>yburnin) recorded<-TRUE else burnin<-burnin+1
    # Storing distribution and population data
    if(recorded==TRUE)
    {
      if(sc <=scfit) # Only records spatial data for fitting set of scenarios
      {
        ydat<-rbind(ydat, data.frame("sc"=rep(sc, d^2),"yr"=rep(y, d^2),"use"=c(trajPre), "food"=c(env1), "temp"=c(env2), "food2"=c(env1^2), "temp2"=c(env2^2),"fit"=c(fitinit), "Efd"=rep(mean(env1), d^2), "Ecv"=rep(mean(env2), d^2), "Efd2"=rep(mean(env1^2), d^2), "Ecv2"=rep(mean(env2^2), d^2), "Epr"=rep(mean(fitinit), d^2), "offs"=rep(log(popc), d^2), "N"=rep(popc, d^2)))
      }
      dat.N<-rbind(dat.N, data.frame("scenario"=sc, "N2"=pop[y], "N1"=pop[y-1], "AvFood"=mean(env1), "AvTemp"=mean(env2),"AvFood2"=mean(env1^2), "AvTemp2"=mean(env2^2), "Avfit"=mean(fitinit)))
      y<-y+1
    }
  } # End of years loop
  
  
  if(recorded==TRUE)
  {
    # Storing environmental makeup parameters
    dat.Master<-rbind(dat.Master, cbind("scenario"=rep(sc, length(ws)), t(mus), "varf"=rep(var[1], length(ws)), "vart"=rep(var[2], length(ws)), "weights"=ws) )
    sc<-sc+1
    
    # Graphical output for scenario  
    par(mfrow=c(3,2))
    image(env1, main="a. Food", axes=FALSE, xlab="Longitude", ylab="Latitude", col=terrain.colors(100), mgp=c(0,0,0))
    image(env2, main="b. Temperature", axes=FALSE, xlab="Longitude", ylab="Latitude", col=terrain.colors(100), mgp=c(0,0,0))
    image(a0+a1*env1+a2*env2+a3*env2^2, main="c. Fitness", axes=FALSE, xlab="Longitude", ylab="Latitude", col=terrain.colors(100), mgp=c(0,0,0))
    image(trajPre, main="d. Usage", axes=FALSE, xlab="Longitude", ylab="Latitude", col=terrain.colors(100), mgp=c(0,0,0))
    Espace.numer(8000, 100,100, env1,env2)
    Espace.approx(8000, 100, 100, mus, var, ws)
    par(mfrow=c(1,1))
  }
  
} # End of scenario loop
dat.N<-cbind(dat.N, "year"=rep(1:(ymax-1),scmax))


###############################################################
#####      2. GENERATION OF FITTING DATA FRAMES           #####
###############################################################

# Thinning by scenario
scmaxi<-20 #scfit # Must be no greater than scfit
ydatRed0<-ydat[ydat$sc<=scmaxi,]

# Thinning by years
stp<-10 # The step size through the recorded years. Values larger than 1 are used to thin the data
year<-rep(1:((ymax-1)*scmaxi), each=d^2)
ydatRedi<-cbind(ydatRed0, "year"=year)
yrsub<-seq(1,max(as.numeric(ydatRedi$year)),stp)
ydatRed<-ydatRedi[ydatRedi$year %in% yrsub,]
ydatRed$year<-as.factor(ydatRed$year)  
ydatRed1<-ydatRed

# Thinning by spatial cells
stpc<-20 # The step size through the cells in the grid. Values larger than 1 are used to thin the data
ydatRed<-ydatRed1[seq(1, length(ydatRed1[,1]),stpc),]


###################################################
#####      3. SPATIAL MODEL FITTING           #####
###################################################

# Fitting Generalised Function Response model (GFR-HSF)
gmod1<-glm(use~Efd+Ecv+Efd2+Ecv2+N+
             food+food:Efd+food:Ecv+food:Efd2+food:Ecv2+food:N+
             temp+temp:Efd+temp:Ecv+temp:Efd2+temp:Ecv2+temp:N+
             temp2+temp2:Efd+temp2:Ecv+temp2:Efd2+temp2:Ecv2+temp2:N+
             offset(offs), data=ydatRed, family=poisson("log") )

ga<-as.numeric(gmod1$coefficients) # Extracts all model coefficients for later use

# Fitting simple Habitat Selection Function, by pooling the data from different scenarios (simple-HSF)
gmod0<-glm(use~food+temp+temp2+offset(offs), data=ydatRed, family=poisson("log") )

AIC(gmod0, gmod1) # Comparison of the two HSFs



########################################################
#####          4. SOME EXPLORATORY PLOTS           #####
########################################################

# Visual comparison of simple simple HSF and GFR-HSF
rng<-c(0,max(max(ydatRed$use), fitted(gmod0), fitted(gmod1)))
par(mfrow=c(2,2))
plot(ydatRed$use, fitted(gmod0), xlim=rng, ylim=rng)
plot(ydatRed$use, fitted(gmod1), xlim=rng, ylim=rng)
smoothScatter(ydatRed$use, fitted(gmod1), xlim=rng, ylim=rng)
abline(0,1)
par(mfrow=c(1,1))


# Spatial snapshots of usage and the corresponding fitted layers from the GFR
par(mfrow=c(5,6))
dj<-round((ymax-1)*scmaxi/(5*6*0.5))
for( i in seq(1,(ymax-1)*scmaxi,dj))
{
  
  real<-matrix(ydatRed1[(i-1)*d^2+(1:d^2),3], nrow=d, ncol=d)
  rang<-c(min(log(real+1)),max(log(real+1)))
  rang<-c(0,2)
  image(log(real+1), main=paste("Real",i), zlim=rang)
  fitd<-matrix(c(predict(gmod1, newdata=ydatRed1[(i-1)*d^2+(1:d^2),], type="response")), nrow=d, ncol=d)
  image(log(fitd+1), main=paste("Fitted",i), zlim=rang)
}
par(mfrow=c(1,1))


# Plot of all population trajectories
scns<-max(dat.N$scenario)

plot(dat.N[,2][dat.N$scenario==1], ylim=c(0,max(dat.N[,2])), col=grey(0.7),lwd=2, type="l", ylab="Population size", xlab="Year")
for(i in (scfit+1):scmax) {lines(dat.N[,2][dat.N$scenario==i], col=grey(0.7),lwd=2)}
for(i in 2:scfit) {lines(dat.N[,2][dat.N$scenario==i], col="red",lwd=2)}



##########################################################################################
#####     5. CALCULATION OF CONSTRUCTED COVARIATES FOR SPATIAL POPULATION MODEL      #####
##########################################################################################


nis<-10 # Number of samples used for parametric bootstrap of HSF parameter estimates

bar<-txtProgressBar(0, 100, style=3)

for(res in 1:nis) # Bootstrap loop
{
  setTxtProgressBar(bar,res)
  gas<-mvrnorm(1, ga, vcov(gmod1)) # Stochastic generation of parameters for bootstrap iteration
  F1<-c()
  F2<-c()
  F11<-c()
  F12<-c()
  F13<-c()
  F3<-c()
  bG<-c()
  bar<-txtProgressBar(0, 100, style=3)
  for(y in 1:length(dat.N[,1]))
  {
    
    scn<-dat.N[y,1]
    
    f.data<-dat.Master[dat.Master$scenario==scn,]
    muf<-f.data$food
    mut<-f.data$temp
    sif<-f.data$varf
    sit<-f.data$vart
    ws<-f.data$weights
    L<-length(muf)
    
    
    p1<-c()
    p21<-c()
    p22<-c()
    p23<-c()
    p2<-c()
    p3<-c()
    
    
    avF<-dat.N[dat.N$scenario==scn,][1,4]
    avT<-dat.N[dat.N$scenario==scn,][1,5]
    avF2<-dat.N[dat.N$scenario==scn,][1,6]
    avT2<-dat.N[dat.N$scenario==scn,][1,7]
    Ni<-dat.N[y,3]
    
    gaGFR<-c(
      gas[1]+gas[2]*avF+gas[3]*avT+gas[4]*avF2+gas[5]*avT2+gas[6]*Ni, 
      gas[7]+gas[10]*avF+gas[11]*avT+gas[12]*avF2+gas[13]*avT2+gas[14]*Ni,
      gas[8]+gas[15]*avF+gas[16]*avT+gas[17]*avF2+gas[18]*avT2+gas[19]*Ni,
      gas[9]+gas[20]*avF+gas[21]*avT+gas[22]*avF2+gas[23]*avT2+gas[24]*Ni
    )
    
    
    for( l in 1:L)
    {
      p1<-c(  p1,  sqrt(2*pi*sif[l])
              *exp(0.5*gaGFR[2]^2*sif[l]+muf[l]*gaGFR[2])
              *sqrt(2*pi*sit[l]/(1-2*gaGFR[4]*sit[l]))
              *exp((0.5*gaGFR[3]^2*sit[l]+gaGFR[3]*mut[l]+gaGFR[4]*mut[l]^2)/(1-2*gaGFR[4]*sit[l]))
      )
      
      p21<-c(  p21,  (gaGFR[2]*sif[l]+muf[l]))  
      
      p22<-c(  p22,  (gaGFR[3]*sit[l]+mut[l])/(1-2*gaGFR[4]*sit[l]))
      
      p23<-c(  p23,  sit[l]/(1-2*gaGFR[4]*sit[l])*(1+(gaGFR[3]*sit[l]+mut[l])^2/(sit[l]*(1-2*gaGFR[4]*sit[l]))))  
      
      p2<-c(  p2,  gaGFR[2]*(gaGFR[2]*sif[l]+muf[l])
              +gaGFR[3]*(gaGFR[3]*sit[l]+mut[l])/(1-2*gaGFR[4]*sit[l])
              +gaGFR[4]*sit[l]/(1-2*gaGFR[4]*sit[l])*(1+(gaGFR[3]*sit[l]+mut[l])^2/(sit[l]*(1-2*gaGFR[4]*sit[l])))  
      )
      p3<-c(  p3,  sqrt(2*pi*sif[l])
              *exp(2*gaGFR[2]^2*sif[l]+2*muf[l]*gaGFR[2])
              *sqrt(2*pi*sit[l]/(1-4*gaGFR[4]*sit[l]))
              *exp((2*gaGFR[3]^2*sit[l]+2*gaGFR[3]*mut[l]+2*gaGFR[4]*mut[l]^2)/(1-4*gaGFR[4]*sit[l]))
      )
      
      
      
    }
    
    F1<-c(F1, sum(ws*p1*p2))
    F11<-c(F11, sum(ws*p1*p21))
    F12<-c(F12, sum(ws*p1*p22))
    F13<-c(F13, sum(ws*p1*p23))
    F2<-c(F2, sum(ws*p1))
    F3<-c(F3, sum(ws*p3))
    bG<-c(bG, 2*pi*sqrt(sif[1]*sit[1]))
  }
  
  
  # Storage of generated covariates
  if(res==1)
  {
    All.dat<-data.frame(cbind(dat.N, F1, F2, F3, "fi1"=F11/F2, "fi2"=F12/F2, "fi3"=F13/F2, bG, "F"=F1/F2, "r"=dat.N$N2/dat.N$N1))
  }
  if(res>1)
  {
    All.dat<-rbind(All.dat,data.frame(cbind(dat.N, F1, F2, F3, "fi1"=F11/F2, "fi2"=F12/F2, "fi3"=F13/F2, bG, "F"=F1/F2, "r"=dat.N$N2/dat.N$N1)))
    
  }
}

fitting.dat<-All.dat[All.dat$scenario<=scmaxi,]
fitting.dat<-fitting.dat[fitting.dat$year %in% yrsub,]
predict.dat<-All.dat[All.dat$scenario>scfit,]

 
#######################################################
#####      6. FITTING THE POPULATION MODELS       #####
#######################################################

yrange<-c(min(-2,log(dat3$N2)-log(dat3$N1)),max(log(dat3$N2)-log(dat3$N1)))
dat3<-na.omit(fitting.dat[fitting.dat[,3]>100,])
dat3<-cbind(dat3, "DD"=1/(d^2)*dat3$bG*dat3$F3/dat3$F2^2)
yrange<-c(-1,1)
yrange<-c(min(log(dat3$N2)-log(dat3$N1)),max(log(dat3$N2)-log(dat3$N1)))
wei<-rep(1/nis, length(dat3[,1]))

par(mfrow=c(1,2))
spatial<-glm(N2~fi1+fi2+fi3+N1:DD, family=poisson(), offset=log(N1), data=dat3, weights=wei)
plot(log(dat3$N2)-log(dat3$N1), log(fitted(spatial))-log(dat3$N1), ylim=yrange, ylab="Predicted (spatial)", xlab="Observed log(r)")
abline(0,1)


mean_field<-glm(N2~AvFood+AvTemp+I(AvTemp^2)+N1, family=poisson(), offset=log(N1), data=dat3, weights=wei)
plot(log(dat3$N2)-log(dat3$N1), log(fitted(mean_field))-log(dat3$N1),  ylim=yrange, ylab="Predicted (mean field)", xlab="Observed log(r)")
abline(0,1)
par(mfrow=c(1,1))


AIC(spatial)
AIC(mean_field)



#######################################################
#####          7. GOODNESS-OF-FIT PLOTS           #####
#######################################################

# Plotting predictions for population growth rate
yrange<-c(min(log(dat3$N2)-log(dat3$N1)),max(log(dat3$N2)-log(dat3$N1)))
par(mfrow=c(1,2))
smoothScatter(log(dat3$N2)-log(dat3$N1), log(fitted(mean_field))-log(dat3$N1),  nrpoints=500, bandwidth=0.02, ylim=yrange, xlim=yrange, xlab="True growth rate log(r)", ylab="log(r) predicted (mean field)")
abline(0,1)
smoothScatter(log(dat3$N2)-log(dat3$N1), log(fitted(spatial))-log(dat3$N1),  nrpoints=500, bandwidth=0.02, ylim=yrange, xlim=yrange, xlab="True growth rate log(r)", ylab="log(r) predicted (spatial)")
abline(0,1)
par(mfrow=c(1,1))

# Plotting predictions for population carrying capacity
dat4<-dat3[dat3$year>10,]
KKmean<-(cbind(dat4[,1]*0+1, dat4$AvFood, dat4$AvTemp, dat4$AvTemp^2)%*%mean_field$coefficients[1:4])/(-mean_field$coefficients[5])
KKspatial<-(cbind(dat4[,1]*0+1, dat4$fi1, dat4$fi2, dat4$fi3)%*%spatial$coefficients[1:4])/(-spatial$coefficients[5]*dat4$DD)
yrange<-c(0.9*min(dat4$N2),1.1*max(dat4$N2))
par(mfrow=c(1,2))
smoothScatter(dat4$N2, KKmean, nrpoints=500, bandwidth=40, ylim=yrange, xlim=yrange, xlab="True carrying capacity (K)", ylab="K predicted (mean field)")
abline(0,1)
smoothScatter(dat4$N2, KKspatial,  nrpoints=500, bandwidth=40, ylim=yrange, xlim=yrange, xlab="True carrying capacity (K)", ylab="K predicted (spatial)")
abline(0,1)
par(mfrow=c(1,1))


#########################################################
#####              8. PREDICTION PLOTS              #####
#########################################################
dat5<-na.omit(predict.dat[predict.dat[,3]>200,])
dat5<-cbind(dat5, "DD"=1/(d^2)*dat5$bG*dat5$F3/dat5$F2^2)

# Plotting predictions for population growth rate
yrange<-c(min(log(dat5$N2)-log(dat5$N1)),max(log(dat5$N2)-log(dat5$N1)))
RRtrue<-log(dat5$N2)-log(dat5$N1)
RRmean<-predict(mean_field, newdat=dat5)-log(dat5$N1)
RRspatial<-predict(spatial, newdat=dat5)-log(dat5$N1)
resCI(RRtrue, RRmean)
resCI(RRtrue, RRspatial)
  
par(mfrow=c(1,2))
smoothScatter(RRtrue, RRmean ,  nrpoints=2000, bandwidth=0.02, ylim=yrange, xlim=yrange, xlab="True growth rate log(r)", ylab="log(r) predicted (mean field)")
abline(0,1)
smoothScatter(RRtrue, RRspatial,  nrpoints=2000, bandwidth=0.02, ylim=yrange, xlim=yrange, xlab="True growth rate log(r)", ylab="log(r) predicted (spatial)")
abline(0,1)
par(mfrow=c(1,1))

# Plotting predictions for population carrying capacity
dat6<-dat5[dat5$year>15,]
KKtrue<-dat6$N2
KKmean<-(cbind(dat6[,1]*0+1, dat6$AvFood, dat6$AvTemp, dat6$AvTemp^2)%*%mean_field$coefficients[1:4])/(-mean_field$coefficients[5])
KKspatial<-(cbind(dat6[,1]*0+1, dat6$fi1, dat6$fi2, dat6$fi3)%*%spatial$coefficients[1:4])/(-spatial$coefficients[5]*dat6$DD)
resCI(KKtrue, KKmean)
resCI(KKtrue, KKspatial)

yrange<-c(0.9*min(dat6$N2),1.1*max(dat6$N2))
par(mfrow=c(1,2))
smoothScatter(KKtrue, KKmean, nrpoints=2000, bandwidth=50, ylim=yrange, xlim=yrange, xlab="True carrying capacity (K)", ylab="K predicted (mean field)")
abline(0,1)
smoothScatter(KKtrue, KKspatial,  nrpoints=2000, bandwidth=50, ylim=yrange, xlim=yrange, xlab="True carrying capacity (K)", ylab="K predicted (spatial)")
abline(0,1)
par(mfrow=c(1,1))