biocro: Simulation of Energycane and Sugarcane

## 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
}