R/weightmaxentropy.R

In embryogrowth: Tools to Analyze the Thermal Reaction Norm of Embryo Growth

#' Search for the weights of the nests which maximize the entropy of nest temperatures distribution
#' @title Search for the weights of the nests which maximize the entropy of nest temperatures distribution
#' @author Marc Girondot
#' @return A named vector of weights
#' @param temperatures Timeseries of temperatures formated using FormatNests()
#' @param entropy.method Entropy function, for example entropy::entropy.empirical. See package entropy for description
#' @param weight A named vector of the initial weight search for each nest for likelihood estimation
#' @param control_optim A list with control paramaters for optim function
#' @param control_plot A list with control paramaters for plot function
#' @param control_entropy A list with control paramaters for entropy function
#' @param plot Do the plot of temperatures before and after weight must be shown ? TRUE or FALSE
#' @param col Colors for unweighted and weighted distributions
#' @description Search for the weights of the nests which maximize the entropy of nest temperatures distribution. Entropy is measured by Shanon index.\cr
#' Entropy method must be entropy.empirical because it is the only method insensitive to scaling.\cr
#' If no weight is given, the initial weight is uniformly distributed.\cr
#' Use control_optim=list(trace=0) for not show progress of search report.
#' @examples
#' \dontrun{
#' library(embryogrowth)
#' data(nest)
#' formated <- FormatNests(nest)
#' w <- weightmaxentropy(temperatures=formated, control_plot=list(xlim=c(20,36)))
#' x <- structure(c(120.940334922916, 467.467455887442,  
#' 	306.176613681557, 117.857995419495),  
#' 	names = c("DHA", "DHH", "T12H", "Rho25"))
#' # pfixed <- c(K=82.33) or rK=82.33/39.33
#' pfixed <- c(rK=2.093313)
#' # K or rK are not used for dydt.linear or dydt.exponential
#' resultNest_4p_weight <- searchR(parameters=x,  
#' 	fixed.parameters=pfixed, temperatures=formated,  
#' 	integral=integral.Gompertz, M0=1.7, hatchling.metric=c(Mean=39.33, SD=1.92),  
#' 	method = "BFGS", weight=w)
#' data(resultNest_4p_weight)
#' plotR(resultNest_4p_weight, ylim=c(0,0.50), xlim=c(15, 35))
#' # Standard error of parameters can use the GRTRN_MHmcmc() function
#' }
#' @export


weightmaxentropy <- function(temperatures=stop('Temperature data must be provided !') , 
                             weight= NULL                                             , 
                             entropy.method=entropy::entropy.empirical                , 
                             plot=TRUE                                                , 
                             control_optim=list(trace=0, maxit=500)                   , 
                             control_plot=NULL                                        , 
                             control_entropy=NULL                                     , 
                             col=c("black", "red")                                    ) {
  
  # 	weight= NULL; entropy.method=entropy.empirical; plot=TRUE; control_optim=list(trace=0, maxit=500); control_plot=NULL; control_entropy=NULL
  
  if (!requireNamespace("entropy", quietly = TRUE)) {
    stop("entropy package is necessary for this function")
  }
  
  
  if (is.null(control_plot)) control_plot <- list(NULL)
  if (is.null(control_entropy)) control_entropy <- list(NULL)
  if (is.null(control_optim)) control_optim <- list(NULL)
  
  NbTS <- temperatures$IndiceT["NbTS"]
  
  
  if (is.null(weight)) {
    par <- rep(1, NbTS)
    names(par) <- temperatures$Names
  } else {
    
    if (any(is.na(weight))) {
      par <- rep(1, NbTS)
      names(par) <- temperatures$Names
    } else {
      
      if (is.list(weight)) weight <- weight$weight
      
      if (length(setdiff(temperatures$Names, names(weight)))==0) {
        par <- weight
      } else {
        print("Check the weights")
        return(invisible())
      }
    }
  }
  
  # j'ai le weight mais je dois le classer dans le mme ordre que les temperatures
  pec <- NULL
  for(i in temperatures$Names) {
    pec <- c(pec, par[i])
  }
  
  par <- pec[1:(NbTS-1)]/pec[NbTS]
  
  vm <- setdiff(temperatures$Names, names(par))
  
  if (length(vm) != 1) {
    print("Check the weights")
    return(invisible())
  }
  
  
  
  
  control_optim <- modifyList(control_optim, list(fnscale=-1))
  
  
  ##################
  # Calcul de la distribution sans weight
  ##################
  
  pi <- NULL
  
  vec.weight <- c(abs(par), 1)
  
  for(series in temperatures$Names) {
    
    nids <- temperatures$Nests[[series]]$data
    
    # Je cree un tableau avec les donnees heures par heure
    tl1 <- (0:(nids[,1][length(nids[,1])]%/%60-1))*60
    # je prends les vraies donnees
    tl2 <- nids[,1]
    
    
    it <- findInterval(tl1, tl2)
    tls <- nids[it, 2]
    
    ctls <- cut(tls, 0:45)
    
    pi <- rbind(pi, table(ctls))
    
  }
  
  rownames(pi) <- names(temperatures)[1:NbTS]
  
  nrow <- dim(pi)[1]
  ncol <- dim(pi)[2]
  
  pisansweight <- NULL
  for(i in 1:ncol) pisansweight <- c(pisansweight, sum(pi[,i]))
  
  
  ##################
  # Recherche des weights pour entropy maximale
  ##################
  
  
  
  print(paste("Initial entropy is", getFromNamespace(".weightmaxentropy_fit", ns="embryogrowth")(par, pi=pi, entropy.method=entropy.method, pentropy=control_entropy, vm=vm)))
  
  
  result  <- optim(par, getFromNamespace(".weightmaxentropy_fit", ns="embryogrowth"), pi=pi, entropy.method=entropy.method, pentropy=control_entropy, vm=vm, control=control_optim, method="BFGS", hessian=FALSE)
  
  if (result$convergence==1) {
    print("Maximum number of iterations has been reached. Try increase maxit in controls")
  }
  
  par <- setNames(abs(c(result$par, 1)), temperatures$Names)
  
  ##################
  # Calcul de la distribution avec weight
  ##################
  
  
  vec.weight <- par
  
  
  # juste pour reprendre la structure
  pi2 <- pi
  
  for(i in 1:nrow)
    pi2[i,] <- pi[i, 1:ncol]*(vec.weight[i]/sum(vec.weight))
  
  piavecweight <- NULL
  for(i in 1:ncol) piavecweight <- c(piavecweight, sum(pi2[,i]))
  
  if (plot) {
    
    L <- modifyList(list(x=0:44, y=pisansweight, type="h", 
                         xlim=c(0,44), lwd=2, bty="n", xlab="Temperatures", 
                         ylab="Relative frequency", col=col[1]), control_plot)
    
    do.call(what="plot", args=L)
    
    # if (par("yaxs")=="i") {
    #    y1 <- par("usr")[3]
    #    y2 <- par("usr")[4]
    #  } else {
    #    y2 <- (par("usr")[3]+par("usr")[4]*26)/27
    #    y1 <- y2*26-par("usr")[4]/0.04
    #  }
    
    par(new=TRUE)
    L <- modifyList(list(x=(0:44)+0.3, y=piavecweight, type="h", 
                         xlim=c(0,44), lwd=2, bty="n", xlab="", axes=FALSE,
                         ylab="", col=col[2]), control_plot)
    
    do.call(what="plot", args=L)
    
    legend("topleft", legend=c("Without weight", "With weight"), col=col, lty=1, lwd=2, bty="n")
    
  }
  
  
  print(paste("Maximal entropy is", result$value))
  
  return(list(weight=abs(par), freq.obs=pisansweight, freq.max.entropy=piavecweight))
  
  
}