chlorophyllCycling/createModelCalibrations_CCmodel.R at main · k-wheeler/chlorophyllCycling · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
##Create Calibration Fits for mean air temperature no trigger model
library(rjags)
library(runjags)
library(doParallel)
source('generalVariables.R')
source('runModelIterations.R')

createChlorophyllCyclingModelCalibration <- function(ns,sVals){
  print(ns)
  generalModel = "
model {
    ### Data Models for complete years
    for(yr in 1:(N)){
    for(i in 1:n){
    p[i,yr] ~ dnorm(x[i,yr],p.PC)
    }
    }

    #### Process Model
    for(yr in 1:(N)){
    for(i in 2:n){
    Tair[i,yr] ~ dnorm(TairMu[i,yr],TairPrec[i,yr])

    xmu[i,yr] <- max(min((x[(i-1),yr] + b4 * x[(i-1),yr]) + max(0,(b0 + b3 * Tair[i,yr] * D[i,yr])),x[1,yr]),0)
    x[i,yr] ~ dnorm(xmu[i,yr],p.proc)
    }
    }

    #### Priors
    for(yr in 1:N){ ##Initial Conditions
    x[1,yr] ~ dbeta(x1.a[yr],x1.b[yr]) I(0.001,0.999)
    }
    p.PC ~ dgamma(s1.PC,s2.PC)
    p.proc ~ dgamma(s1.proc,s2.proc)

    b0 ~ dunif(b0_lower,b0_upper)
    b3 ~ dunif(b3_lower,b3_upper)
    b4 ~ dunif(b4_lower,b4_upper)

  }"

  registerDoParallel(cores=n.cores)
  variableNames <- c("p.PC","x","p.proc","b0","b3","b4")

  foreach(s =sVals) %dopar% {
    siteName <- sites[s]
    print(siteName)
    yearRemoved <- yearsRemoved[s]
    load(paste0(dataDirectory,siteName,"_dataFinal.RData"))

    for(n in ns){
      print(n)
      outputFileName <- paste0(CCmodelOutputsFolder,siteName,"_",n,"_ccModel_forecast_calibration_varBurn.RData")
      dataFinal$p[(n+1):nrow(dataFinal$p),] <- NA

      #Remove year
      yearInt <- which(dataFinal$years==yearRemoved)
      dataFinal$p[,yearInt] <- NA

      ##Add priors
      dataFinal$s1.PC <- 1.56
      dataFinal$s2.PC <- 0.016
      dataFinal$s1.proc <- 1.56
      dataFinal$s2.proc <- 0.016
      dataFinal$b0_lower <- -1 #intercept of synthesis
      dataFinal$b0_upper <- 1 #intercept of synthesis

      dataFinal$b3_lower <- 0 #slope of synthesis
      dataFinal$b3_upper <- 1
      dataFinal$b4_lower <- -1 #Breakdown rate
      dataFinal$b4_upper <- 0 #Breakdown rate

      j.model <- try(jags.model(file = textConnection(generalModel),
                                data = dataFinal,
                                n.chains = nchain,
                                n.adapt = 3000))#Load Model Output
      if(inherits(j.model,"try-error")){
        next
      }

      out.burn <- try(runForecastIter(j.model=j.model,variableNames=variableNames,
                                      baseNum = 15000,iterSize = 5000,effSize = 5000, maxIter=1000000))
      if(inherits(out.burn,"try-error")){
        next
      }

      ##Thin the data:
      if(typeof(out.burn)!=typeof(FALSE)){
        out.mat <- as.matrix(out.burn$params)
        thinAmount <- round(nrow(out.mat)/5000,digits=0)
        out.burn2 <- list()
        out.burn2$params <- window(out.burn$params,thin=thinAmount)
        out.burn2$predict <- window(out.burn$predict,thin=thinAmount)
        out.burn <- out.burn2
      }
      save(out.burn,file = outputFileName)
      print(paste("saved:",outputFileName))
    }
  }

}

createChlorophyllCyclingModelCalibration(ns=seq(45,183),
                                         sVals=seq_along(sites)) #Change for number of included days (ns) and sites (sVals)