R/FitDoubleLogKlHeavy.R

In phenopix: Process Digital Images of a Vegetation Cover

FitDoubleLogKlHeavy <-
  function(
    ##title<<
    ## Fit a double logisitic function to a vector according to Elmore et al. (2012)
    ##description<<
    ## This function fits a double logistic curve to observed values using the function as described in Elmore et al. (2012) (equation 4).
    
    x,
    ### vector or time series to fit
    
    t = index(x),
    ### time steps
    
    tout = t,
    ### time steps of output (can be used for interpolation)
    
    max.iter=200,
    sf=quantile(x, probs=c(0.05, 0.95), na.rm=TRUE), 
    ### if TRUE the function returns parameters of the double logisitic fit, if FALSE it returns the fitted curve
    
    ### plot iterations for logistic fit?
    
    ...
    ### further arguments (currently not used)
    
    ##references<<
    ## Elmore, A.J., S.M. Guinn, B.J. Minsley and A.D. Richardson (2012): Landscape controls on the timing of spring, autumn, and growing season length in mid-Atlantic forests. - Global Change Biology 18, 656-674.
    
    ##seealso<<
    ## \code{\link{TSGFdoublelog}}, \code{\link{Phenology}}
    
  ) {
    .normalize <- function(x, sf) (x-sf[1])/(sf[2]-sf[1])
    .backnormalize <- function(x, sf) (x+sf[1]/(sf[2]-sf[1]))*(sf[2]-sf[1])    
    if (any(is.na(x))) stop('NA in the time series are not allowed: fill them with e.g. na.approx()')
    if (inherits(index(x), 'POSIXct')) {
      doy.vector <- as.numeric(format(index(x), '%j'))
      index(x) <- doy.vector
    }
    ## integrare da tmp.R per la funzione best.nls e tutte le sue dipendenze
    ## la teniamo strutturata cos e aggiustiamo i dati in modo che funzioni
    
    .rss <- function(fit, gcc) {
      sum.all <- sum(abs(predict(fit)-gcc))/sqrt(length(gcc))
    }
    
    .best.nls <- function(data, comb, max.iter=NULL, params, plot=FALSE) {
      the.funct <- formula(max.filtered~ (a1*t + b1) + (a2*t^2 + b2*t + c)*(1/(1+q1*exp(-B1*(t-m1)))^v1 - 1/(1+q2*exp(-B2*(t-m2)))^v2))
      
      #     offset + (upper/(1+q1*exp(-b1*(t-m1)))^v1 -
      #                                             upper/(1+q2*exp(-b2*(t-m2)))^v2))
      uniqued <- comb
      if (is.null(max.iter)) {counter <- dim(uniqued)[1]} else {
        counter <- max.iter
        sample.vect <- sample(1:dim(uniqued)[1], max.iter)
        uniqued <- uniqued[sample.vect, ]
      }
      all.rss <- numeric(counter)
      
      for (a in 1:counter) {
        row <- uniqued[a,]
        cfc <- params[which(names(params) %in% row[1]==TRUE)]
        fit.klos1 <- try(nls(the.funct, data=data, start=cfc), silent=TRUE)
        if (inherits(fit.klos1, 'try-error')) {
          rss.final <- 999 } else {
            rss1 <- .rss(fit.klos1, data$max.filtered)
            cfc <- coefficients(fit.klos1)
            cfc[2] <- params[which(names(params) %in% row[2]==TRUE)]
            names(cfc)[2] <- row[2]
            ##step 2 offset and b1
            fit.klos2 <- try(nls(the.funct, data=data, start=cfc), silent=TRUE)
            if (inherits(fit.klos2, 'try-error')) {
              last.fit <- fit.klos1
              rss.final <- rss1
            } else {
              rss2 <- .rss(fit.klos2, data$max.filtered)
              cfc <- coefficients(fit.klos2)
              cfc[3] <-  params[which(names(params) %in% row[3]==TRUE)]
              names(cfc)[3] <- row[3]
              ##step 3 offset b1 b2
              fit.klos3 <- try(nls(the.funct, data=data, start=cfc), silent=TRUE)
              if (inherits(fit.klos3, 'try-error')) {
                last.fit <- fit.klos2
                rss.final <- rss2
              } else {
                rss3 <- .rss(fit.klos3, data$max.filtered)
                cfc <- coefficients(fit.klos3)
                cfc[4] <- params[which(names(params) %in% row[4]==TRUE)]
                names(cfc)[4] <- row[4]
                ##step 4 offset b1 b2 q1
                fit.klos4 <- try(nls(the.funct, data=data, start=cfc), silent=TRUE)
                if (inherits(fit.klos4, 'try-error')) {
                  last.fit <- fit.klos3
                  rss.final <- rss3
                } else {
                  rss4 <- .rss(fit.klos4, data$max.filtered)
                  cfc <- coefficients(fit.klos4)
                  cfc[5] <- params[which(names(params) %in% row[5]==TRUE)]
                  names(cfc)[5] <- row[5]
                  ##step 5 offset b1 b2 q1 q2
                  fit.klos5 <- try(nls(the.funct, data=data, start=cfc), silent=TRUE)
                  if (inherits(fit.klos5, 'try-error')) {
                    last.fit <- fit.klos4
                    rss.final <- rss4
                  } else {
                    rss5 <- .rss(fit.klos5, data$max.filtered)
                    cfc <- coefficients(fit.klos5)
                    cfc[6] <- params[which(names(params) %in% row[6]==TRUE)]
                    names(cfc)[6] <- row[6]
                    ##step 6
                    fit.klos6 <- try(nls(the.funct, data=data,
                                         start=cfc), silent=TRUE)
                    if (inherits(fit.klos6, 'try-error')) {
                      last.fit <- fit.klos5
                      rss.final <- rss5
                    } else {
                      rss6 <- .rss(fit.klos6, data$max.filtered)
                      cfc <- coefficients(fit.klos6)
                      cfc[7] <- params[which(names(params) %in% row[7]==TRUE)]
                      names(cfc)[7] <- row[7]
                      ##step 6
                      fit.klos7 <- try(nls(the.funct, data=data,
                                           start=cfc), silent=TRUE)
                      if (inherits(fit.klos7, 'try-error')) {
                        last.fit <- fit.klos6
                        rss.final <- rss6
                      } else {
                        rss7 <- .rss(fit.klos7, data$max.filtered)
                        cfc <- coefficients(fit.klos7)
                        cfc[8] <- params[which(names(params) %in% row[8]==TRUE)]
                        names(cfc)[8] <- row[8]
                        ##step 6
                        fit.klos8 <- try(nls(the.funct, data=data,
                                             start=cfc), silent=TRUE)
                        if (inherits(fit.klos8, 'try-error')) {
                          last.fit <- fit.klos7
                          rss.final <- rss7
                        } else {
                          rss8 <- .rss(fit.klos8, data$max.filtered)
                          cfc <- coefficients(fit.klos8)
                          cfc[9] <- params[which(names(params) %in% row[9]==TRUE)]
                          names(cfc)[9] <- row[9]
                          ##step 6
                          fit.klos9 <- try(nls(the.funct, data=data,
                                               start=cfc), silent=TRUE)
                          if (inherits(fit.klos9, 'try-error')) {
                            last.fit <- fit.klos8
                            rss.final <- rss8
                          } else {
                            rss9 <- .rss(fit.klos9, data$max.filtered)
                            cfc <- coefficients(fit.klos9)
                            cfc[10] <- params[which(names(params) %in% row[10]==TRUE)]
                            names(cfc)[10] <- row[10]
                            ##step 6
                            fit.klos10 <- try(nls(the.funct, data=data,
                                                  start=cfc), silent=TRUE)
                            if (inherits(fit.klos10, 'try-error')) {                                                    last.fit <- fit.klos9
                            rss.final <- rss9
                            } else {
                              last.fit <- fit.klos10
                              rss.final <- .rss(fit.klos10, data$max.filtered)
                            }
                          }
                        }
                      }
                    }
                  }
                }
              }
            }
          }
        ## lines(t, predict(last.fit), col=colors()[a])
        all.rss[a] <- rss.final
        #     print(paste(round(a/counter*100), '% iteractions done', sep=''))
      }
      
      rss.min <- all.rss[which.min(all.rss)]
      min.pos <- which.min(all.rss)
      row <- uniqued[min.pos,]
      cfc <- params[which(names(params) %in% row[1]==TRUE)]
      fit.klos1 <- try(nls(the.funct, data=data,
                           start=cfc), silent=TRUE)
      if (inherits(fit.klos1, 'try-error')) {
        rss.final <- 999 } else {
          rss1 <- .rss(fit.klos1, gcc=data$max.filtered)
          cfc <- coefficients(fit.klos1)
          cfc[2] <- params[which(names(params) %in% row[2]==TRUE)]
          names(cfc)[2] <- row[2]
          ##step 2 offset and b1
          fit.klos2 <- try(nls(the.funct, data=data,
                               start=cfc), silent=TRUE)
          if (inherits(fit.klos2, 'try-error')) {
            last.fit <- fit.klos1
            rss.final <- rss1
          } else {
            rss2 <- .rss(fit.klos2, data$max.filtered)
            cfc <- coefficients(fit.klos2)
            cfc[3] <-  params[which(names(params) %in% row[3]==TRUE)]
            names(cfc)[3] <- row[3]
            ##step 3 offset b1 b2
            fit.klos3 <- try(nls(the.funct, data=data,
                                 start=cfc), silent=TRUE)
            if (inherits(fit.klos3, 'try-error')) {
              last.fit <- fit.klos2
              rss.final <- rss2
            } else {
              rss3 <- .rss(fit.klos3, data$max.filtered)
              cfc <- coefficients(fit.klos3)
              cfc[4] <- params[which(names(params) %in% row[4]==TRUE)]
              names(cfc)[4] <- row[4]
              ##step 4 offset b1 b2 q1
              fit.klos4 <- try(nls(the.funct, data=data,
                                   start=cfc), silent=TRUE)
              if (inherits(fit.klos4, 'try-error')) {
                last.fit <- fit.klos3
                rss.final <- rss3
              } else {
                rss4 <- .rss(fit.klos4, data$max.filtered)
                cfc <- coefficients(fit.klos4)
                cfc[5] <- params[which(names(params) %in% row[5]==TRUE)]
                names(cfc)[5] <- row[5]
                ##step 5 offset b1 b2 q1 q2
                fit.klos5 <- try(nls(the.funct, data=data,
                                     start=cfc), silent=TRUE)
                if (inherits(fit.klos5, 'try-error')) {
                  last.fit <- fit.klos4
                  rss.final <- rss4
                } else {
                  rss5 <- .rss(fit.klos5, data$max.filtered)
                  cfc <- coefficients(fit.klos5)
                  cfc[6] <- params[which(names(params) %in% row[6]==TRUE)]
                  names(cfc)[6] <- row[6]
                  ##step 6
                  fit.klos6 <- try(nls(the.funct, data=data,
                                       start=cfc), silent=TRUE)
                  if (inherits(fit.klos6, 'try-error')) {
                    last.fit <- fit.klos5
                    rss.final <- rss5
                  } else {
                    rss6 <- .rss(fit.klos6, data$max.filtered)
                    cfc <- coefficients(fit.klos6)
                    cfc[7] <- params[which(names(params) %in% row[7]==TRUE)]
                    names(cfc)[7] <- row[7]
                    ##step 6
                    fit.klos7 <- try(nls(the.funct, data=data,
                                         start=cfc), silent=TRUE)
                    if (inherits(fit.klos7, 'try-error')) {
                      last.fit <- fit.klos6
                      rss.final <- rss6
                    } else {
                      rss7 <- .rss(fit.klos7, data$max.filtered)
                      cfc <- coefficients(fit.klos7)
                      cfc[8] <- params[which(names(params) %in% row[8]==TRUE)]
                      names(cfc)[8] <- row[8]
                      ##step 6
                      fit.klos8 <- try(nls(the.funct, data=data,
                                           start=cfc), silent=TRUE)
                      if (inherits(fit.klos8, 'try-error')) {
                        last.fit <- fit.klos7
                        rss.final <- rss7
                      } else {
                        rss8 <- .rss(fit.klos8, data$max.filtered)
                        cfc <- coefficients(fit.klos8)
                        cfc[9] <- params[which(names(params) %in% row[9]==TRUE)]
                        names(cfc)[9] <- row[9]
                        ##step 6
                        fit.klos9 <- try(nls(the.funct, data=data,
                                             start=cfc), silent=TRUE)
                        if (inherits(fit.klos9, 'try-error')) {
                          last.fit <- fit.klos8
                          rss.final <- rss8
                        } else {
                          rss9 <- .rss(fit.klos9, data$max.filtered)
                          cfc <- coefficients(fit.klos9)
                          cfc[10] <- params[which(names(params) %in% row[10]==TRUE)]
                          names(cfc)[10] <- row[10]
                          ##step 6
                          fit.klos10 <- try(nls(the.funct, data=data,
                                                start=cfc), silent=TRUE)
                          if (inherits(fit.klos10, 'try-error')) {
                            last.fit <- fit.klos9
                            rss.final <- rss9
                          } else {
                            last.fit <- fit.klos10
                            rss.final <- .rss(fit.klos10, data$max.filtered)
                          }
                        }
                      }
                    }
                  }
                }
              }
            }
          }
        }
      output <- list(rss=rss.final, fit=last.fit, par=params)
      ## gestire l'estrazione dei massimi e minimi con local maxima e minima
    }
    
    ## format data
    
    n <- length(x)
    x <- .normalize(x, sf=sf)
    avg <- mean(x, na.rm=TRUE)
    mx <- max(x, na.rm=TRUE)
    mn <- min(x, na.rm=TRUE)
    ampl <- mx - mn
    
    ## get parameters
    doy <- quantile(t, c(0.25, 0.75), na.rm=TRUE)
    # doy2 <- diff(doy)
    a1 <- 0 #ok
    a2 <- 0 #ok
    b1 <- mn #ok
    b2 <- 0 #ok
    c <- 0.2*max(x) # ok 
    ## very slightly smoothed spline to get reliable maximum
    tmp <- smooth.spline(x, df=0.5*length(x))
    doy.max <- which.max(tmp$y)
    B1 <- 4/(doy.max-doy[1])
    B2 <- -4/(doy[2]-doy.max)
    m1 <- doy[1] + 0.5*(doy.max-doy[1])
    m2 <- doy.max + 0.5*(doy[2]-doy.max)
    m1.bis <- doy[1]
    m2.bis  <- doy[2]  
    q1 <- 0.5 #ok
    q2 <- 0.5 #ok
    v1 <- 2 # ok
    v2 <- 2 # ok
    
    # offset <- min(ts, na.rm=T)
    # upper <- diff(range(ts, na.rm=T))
    # ## offset <- 0.3
    # a1 <- 0.00003
    # b1 <- 0.05
    # a2 <- 0.001
    # b2 <- 0.08
    # c <- 0
    # q1 <- 1
    # q2 <- 1
    # m1 <- 130
    # m2 <- 250
    # v1 <- 1.3
    # v2 <- 1
    ## upper <- 0.05
    all.pars <- list(a1=a1, a2=a2, b1=b1, b2=b2, c=c, B1=B1, B2=B2, m1=m1, m2=m2, q1=q1,q2=q2, v1=v1, v2=v2)
    data <- data.frame(max.filtered=as.vector(x), t=index(x))
    ## days <- tt$days
    ## gcc <- tt$max.filtered
    
    par.names <- c('a1', 'a2', 'b1','b2', 'c', 'B1','B2', 'm1', 'm2', 'q1','q2', 'v1', 'v2')
    ## get a random sample of parameters (ten times larger than actual number of max.iter so that we make sure 
    ## we can sample at least max.iter of unique()d parameters order)
    
    reps <- max.iter *10
    
    df.pars <- data.frame(matrix(nrow=reps, ncol=length(par.names)))
    for (a in 1:reps) {
      df.pars[a, ] <- unlist(sample(par.names,length(par.names)))
    }
    ## remove replicates if any
    uniqued.tmp <- unique(df.pars)
    ## sample n=max.iter different combinations of params
    uniqued <- uniqued.tmp[sample(dim(uniqued.tmp)[1], 200), ]
    
    nls.best <- .best.nls(data, comb=uniqued, max.iter=max.iter, params=all.pars)
    predicted <- predict(nls.best$fit)
    predicted <- .backnormalize(predicted, sf=sf)
    predicted <- zoo(predicted, order.by=t)
    pars <- nls.best$par
    par.choosen <- pars
    for (a in 1:length(pars)) {
      act.name <- names(pars)[a]
      check <- act.name %in% names(coef(nls.best$fit))
      if (check==TRUE) {
        pos.in.fit <- names(coef(nls.best$fit)) %in% act.name 
        par.choosen[[a]]  <- coef(nls.best$fit)[pos.in.fit]
      }
    }
    the.names <- names(par.choosen)
    unlisted.pars <- as.vector(unlist(par.choosen))
    names(unlisted.pars) <- the.names
    fit.formula <- expression((a1*t + b1) + (a2*t^2 + b2*t + c)*(1/(1+q1*exp(-B1*(t-m1)))^v1 - 1/(1+q2*exp(-B2*(t-m2)))^v2))
    output <- list(predicted=predicted, params=unlisted.pars, formula=fit.formula, sf=sf)
    return(output)
  }