## BioCro/R/constrOpBioGro.R by Fernando Ezequiel Miguez Copyright (C)
## 2007-2009
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the Free
## Software Foundation; either version 2 or 3 of the License (at your option).
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
## more details.
##
## A copy of the GNU General Public License is available at
## http://www.r-project.org/Licenses/
## This will be hopefully a more clever way of approaching the Optimization of
## Biomass Growth The design idea is to optimize one phenological stage at a
## time It will be mostly written in R as I don't see big problems and most of
## the simulation time goes in running the BioGro function which is almost
## completely written in C
constrOpBioGro <- function(phen = 1, iCoef = NULL, WetDat, data, day1 = NULL, dayn = NULL,
timestep = 1, lat = 40, iRhizome = 7, irtl = 1e-04, canopyControl = list(), seneControl = list(),
photoControl = list(), phenoControl = list(), soilControl = list(), nitroControl = list(),
centuryControl = list(), verbose = FALSE, ...) {
## This seems like repeated code from BioGro but it is needed in order to
## automatically pick the right elements from the data
if (missing(day1)) {
half <- as.integer(dim(WetDat)[1]/2)
WetDat1 <- WetDat[1:half, c(2, 5)]
WetDat1s <- WetDat1[which(WetDat1[, 2] < 0), ]
day1 <- max(WetDat1s[, 1])
if (day1 < 90)
day1 <- 90
}
if (missing(dayn)) {
half <- as.integer(dim(WetDat)[1]/2)
WetDat1 <- WetDat[half:dim(WetDat)[1], c(2, 5)]
WetDat1s <- WetDat1[which(WetDat1[, 2] < 0), ]
dayn <- min(WetDat1s[, 1])
if (dayn > 330)
dayn <- 330
}
if ((day1 < 0) || (day1 > 365) || (dayn < 0) || (dayn > 365))
stop("day1 and dayn should be between 0 and 365")
if (day1 > dayn)
stop("day1 should be smaller than dayn")
if ((timestep < 1) || (24%%timestep != 0))
stop("timestep should be a divisor of 24 (e.g. 1,2,3,4,6,etc.)")
vecsize <- (dayn - (day1 - 1)) * 24
indes1 <- (day1 - 1) * 24
indesn <- (dayn) * 24
Temp <- WetDat[indes1:indesn, 5]
Temp12 <- ifelse(Temp >= 12, as.numeric(Temp - 12), 0)
cTT <- cumsum(Temp12/24)
# cTT <- cumsum(Temp/24)
indTmp <- BioCro:::indfun(data[, 1], cTT)
doy <- WetDat[indes1:indesn, 2]
doyTemp <- cbind(doy, cTT)
indDoy <- doyTemp[indTmp, ]
hour <- WetDat[indes1:indesn, 3]
solar <- WetDat[indes1:indesn, 4]
rh <- WetDat[indes1:indesn, 6]
WindS <- WetDat[indes1:indesn, 6]
## Getting the Parameters
canopyP <- canopyParms()
canopyP[names(canopyControl)] <- canopyControl
soilP <- soilParms()
soilP[names(soilControl)] <- soilControl
nitroP <- nitroParms()
nitroP[names(nitroControl)] <- nitroControl
phenoP <- phenoParms()
phenoP[names(phenoControl)] <- phenoControl
photoP <- photoParms()
photoP[names(photoControl)] <- photoControl
seneP <- seneParms()
seneP[names(seneControl)] <- seneControl
centuryP <- centuryParms()
centuryP[names(centuryControl)] <- centuryControl
if (missing(iCoef)) {
iCoef <- as.vector(unlist(phenoP[7:31]))
} else {
if (length(iCoef) != 25)
stop("iCoef should be of length 25")
}
ThermalP <- as.vector(unlist(phenoP)[1:6]) + 1
if (phen == 0)
pheno <- TRUE else pheno <- FALSE
if (phen == 0)
convs <- numeric(6)
if ((phen == 1) || pheno) {
dayn0 <- indDoy[indDoy[, 2] <= ThermalP[1], ]
if (class(dayn0) == "numeric") {
dayn0 <- matrix(dayn0, nrow = 1, ncol = 2)
}
dayn0 <- as.vector(dayn0[NROW(dayn0), 1]) + 1
ui <- matrix(c(-1, -1, 0, 0, 0, -1, 1, 0, 0, 0, 1, 0), ncol = 3, byrow = TRUE)
ci <- c(-1, 0, 0, 0)
opar <- constrOptim(iCoef[c(1:2, 4)], BioCro:::objFun2, phenStage = 1, grad = NULL,
ui = ui, ci = ci, iCoefs = iCoef, ThermalP = ThermalP, WetDat = WetDat,
day1 = day1, dayn = dayn0, lat = lat, irtl = irtl, iRhizome = iRhizome,
timestep = timestep, canopyP = canopyP, seneP = seneP, photoP = photoP,
phenoP = phenoP, soilP = soilP, nitroP = nitroP, data = data, verbose = verbose,
...)
iCoef[c(1:2, 4)] <- opar$par[1:3]
iCoef[3] <- 1 - sum(iCoef[1:2])
iCoef[5:25] <- iCoef[5:25]
if (pheno) {
convs[1] <- opar$convergence
cat(" Stage 1. Converged: ", ifelse(convs[1] == 0, "YES... \n", "NO... \n"))
}
}
if ((phen == 2) || pheno) {
dayn0 <- indDoy[indDoy[, 2] <= ThermalP[2], ]
dayn0 <- as.vector(dayn0[NROW(dayn0), 1]) + 1
ui <- matrix(c(-1, -1, 0, 0, 0, -1, 1, 0, 0, 0, 1, 0), ncol = 3, byrow = TRUE)
ci <- c(-1, 0, 0, 0)
opar <- constrOptim(iCoef[c(5, 6, 8)], BioCro:::objFun2, phenStage = 2, grad = NULL,
ui = ui, ci = ci, iCoefs = iCoef, ThermalP = ThermalP, WetDat = WetDat,
day1 = day1, dayn = dayn0, lat = lat, irtl = irtl, iRhizome = iRhizome,
timestep = timestep, canopyP = canopyP, seneP = seneP, photoP = photoP,
phenoP = phenoP, soilP = soilP, nitroP = nitroP, data = data, verbose = verbose,
...)
iCoef[1:4] <- iCoef[1:4]
iCoef[c(5, 6, 8)] <- opar$par[1:3]
iCoef[7] <- 1 - sum(opar$par[1:2])
iCoef[9:25] <- iCoef[9:25]
if (pheno) {
convs[2] <- opar$convergence
cat(" Stage 2. Converged: ", ifelse(convs[2] == 0, "YES... \n", "NO... \n"))
}
}
if ((phen == 3) || pheno) {
dayn0 <- indDoy[indDoy[, 2] <= ThermalP[3], ]
dayn0 <- as.vector(dayn0[NROW(dayn0), 1]) + 1
ui <- matrix(c(-1, -1, -1, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol = 3, byrow = TRUE)
ci <- c(-1, 0, 0, 0)
opar <- constrOptim(iCoef[9:11], BioCro:::objFun2, grad = NULL, ui = ui,
ci = ci, phenStage = 3, iCoefs = iCoef, WetDat = WetDat, day1 = day1,
dayn = dayn0, lat = lat, ThermalP = ThermalP, irtl = irtl, iRhizome = iRhizome,
timestep = timestep, canopyP = canopyP, seneP = seneP, photoP = photoP,
phenoP = phenoP, soilP = soilP, nitroP = nitroP, data = data, verbose = verbose,
...)
iCoef[1:8] <- iCoef[1:8]
iCoef[9:11] <- opar$par[1:3]
iCoef[12] <- 1 - sum(opar$par[1:3])
iCoef[13:25] <- iCoef[13:25]
if (pheno) {
convs[3] <- opar$convergence
cat(" Stage 3. Converged: ", ifelse(convs[3] == 0, "YES... \n", "NO... \n"))
}
}
if ((phen == 4) || pheno) {
dayn0 <- indDoy[indDoy[, 2] <= ThermalP[4], ]
dayn0 <- dayn0[NROW(dayn0), 1] + 1
ui <- matrix(c(-1, -1, -1, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol = 3, byrow = TRUE)
ci <- c(-1, 0, 0, 0)
opar <- constrOptim(iCoef[13:15], BioCro:::objFun2, grad = NULL, ui = ui,
ci = ci, phenStage = 4, iCoefs = iCoef, WetDat = WetDat, day1 = day1,
dayn = dayn0, lat = lat, ThermalP = ThermalP, irtl = irtl, iRhizome = iRhizome,
timestep = timestep, canopyP = canopyP, seneP = seneP, photoP = photoP,
phenoP = phenoP, soilP = soilP, nitroP = nitroP, data = data, verbose = verbose,
...)
iCoef[1:12] <- iCoef[1:12]
iCoef[13:15] <- opar$par[1:3]
iCoef[16] <- 1 - sum(iCoef[13:15])
iCoef[17:25] <- iCoef[17:25]
if (pheno) {
convs[4] <- opar$convergence
cat(" Stage 4. Converged: ", ifelse(convs[4] == 0, "YES... \n", "NO... \n"))
}
}
if ((phen == 5) || pheno) {
dayn0 <- indDoy[indDoy[, 2] <= ThermalP[5], ]
dayn0 <- as.vector(dayn0[NROW(dayn0), 1]) + 1
ui <- matrix(c(-1, -1, -1, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol = 3, byrow = TRUE)
ci <- c(-1, 0, 0, 0)
opar <- constrOptim(iCoef[17:19], BioCro:::objFun2, grad = NULL, ui = ui,
ci = ci, phenStage = 5, iCoefs = iCoef, WetDat = WetDat, day1 = day1,
dayn = dayn0, lat = lat, ThermalP = ThermalP, irtl = irtl, iRhizome = iRhizome,
timestep = timestep, canopyP = canopyP, seneP = seneP, photoP = photoP,
phenoP = phenoP, soilP = soilP, nitroP = nitroP, data = data, verbose = verbose,
...)
iCoef[1:16] <- iCoef[1:16]
iCoef[17:19] <- opar$par[1:3]
iCoef[20] <- 1 - sum(iCoef[17:19])
iCoef[21:25] <- iCoef[21:25]
if (pheno) {
convs[5] <- opar$convergence
cat(" Stage 5. Converged: ", ifelse(convs[5] == 0, "YES... \n", "NO... \n"))
}
}
if ((phen == 6) || pheno) {
dayn0 <- indDoy[indDoy[, 2] <= ThermalP[6], ]
dayn0 <- as.vector(dayn0[NROW(dayn0), 1]) + 1
ui <- matrix(c(-1, -1, -1, -1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,
0, 1), ncol = 4, byrow = TRUE)
ci <- c(-1 + -1e-07, 0, 0, 0, 0)
opar <- constrOptim(iCoef[21:24], BioCro:::objFun2, grad = NULL, ui = ui,
ci = ci, phenStage = 6, iCoefs = iCoef, WetDat = WetDat, day1 = day1,
dayn = dayn0, lat = lat, ThermalP = ThermalP, irtl = irtl, iRhizome = iRhizome,
timestep = timestep, canopyP = canopyP, seneP = seneP, photoP = photoP,
phenoP = phenoP, soilP = soilP, nitroP = nitroP, data = data, verbose = verbose,
...)
iCoef[1:20] <- iCoef[1:20]
iCoef[21:24] <- opar$par[1:4]
iCoef[25] <- 1 - sum(iCoef[21:24])
if (pheno) {
convs[6] <- opar$convergence
cat(" Stage 6. Converged: ", ifelse(convs[6] == 0, "YES. \n", "NO. \n"))
}
}
if (pheno) {
phen <- 1:6
opar$convergence <- convs
}
list1 <- list(ThermalP = ThermalP, WetDat = WetDat, lat = lat, day1 = day1, timestep = timestep,
dayn = dayn, iRhizome = iRhizome, irtl = irtl, indTmp = indTmp, canopyP = canopyP,
seneP = seneP, photoP = photoP, phenoP = phenoP, soilP = soilP, nitroP = nitroP,
centuryP = centuryP)
structure(list(coefs = iCoef, data = data, opar = opar, phen = phen, list1 = list1),
class = "OpBioGro")
}
objFun2 <- function(coefficients, phenStage, iCoefs, ThermalP, data, WetDat, day1,
dayn, timestep, lat, iRhizome, irtl, canopyP, seneP, photoP, phenoP, soilP, nitroP,
verbose, ...) {
Coefs <- numeric(25)
if (phenStage == 1) {
Coefs[c(1, 2, 4)] <- coefficients[1:3]
Coefs[3] <- 1 - sum(coefficients[1:2])
Coefs[5:25] <- iCoefs[5:25]
dat <- data[data[, 1] <= ThermalP[1], ]
nrdat <- nrow(dat)
} else if (phenStage == 2) {
Coefs[1:4] <- iCoefs[1:4]
Coefs[c(5, 6, 8)] <- coefficients[1:3]
Coefs[7] <- 1 - sum(coefficients[1:2])
Coefs[9:25] <- iCoefs[9:25]
dat <- data[data[, 1] <= ThermalP[2], ]
nrdat <- nrow(dat)
} else if (phenStage == 3) {
Coefs[1:12] <- iCoefs[1:12]
Coefs[9:11] <- coefficients[1:3]
Coefs[12] <- 1 - sum(coefficients[1:3])
Coefs[13:25] <- iCoefs[13:25]
dat <- data[data[, 1] <= ThermalP[3], ]
nrdat <- nrow(dat)
} else if (phenStage == 4) {
Coefs[1:12] <- iCoefs[1:12]
Coefs[13:15] <- coefficients[1:3]
Coefs[16] <- 1 - sum(coefficients[1:3])
Coefs[17:25] <- iCoefs[17:25]
dat <- data[data[, 1] <= ThermalP[4], ]
nrdat <- nrow(dat)
} else if (phenStage == 5) {
Coefs[1:16] <- iCoefs[1:16]
Coefs[17:19] <- coefficients[1:3]
Coefs[20] <- 1 - sum(coefficients[1:3])
Coefs[21:25] <- iCoefs[21:25]
dat <- data[data[, 1] <= ThermalP[5], ]
nrdat <- nrow(dat)
} else if (phenStage == 6) {
Coefs[1:20] <- iCoefs[1:20]
Coefs[21:24] <- coefficients[1:4]
Coefs[25] <- 1 - sum(coefficients[1:4])
dat <- data[data[, 1] <= ThermalP[6], ]
nrdat <- nrow(dat)
}
phenoP[7:31] <- Coefs[7:31 - 6]
res <- BioGro(WetDat, day1, dayn, timestep, lat, iRhizome, irtl, canopyControl = canopyP,
seneControl = seneP, photoControl = photoP, phenoControl = phenoP, soilControl = soilP,
nitroControl = nitroP)
rss <- RssBioGro(dat, res)
if (verbose) {
cat("rss : ", rss, "\n")
print(matrix(Coefs, ncol = 4, byrow = TRUE))
}
rss
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.