R/activity_code.r

In activity: Animal Activity Statistics

#' Animal activity statistics
#'
#' Provides functions to estimate and compare activity parameters from sensor data.
#'
#' @details Sensors that record active animals (eg camera traps) build up a record of
#' the distribution of activity over the course of the day. Records are more frequent
#' when animals are more active, and less frequent or absent when animals are inactive.
#' The area under the distribution of records thus contains information on the overall
#' level of activity in a sampled population. This package provides tools for plotting
#' activity distributions, quantifying the overall level of activity with error, and
#' statistically comparing distributions through bootstrapping.
#'
#' The core function is \code{fitact}, which creates an \code{actmod} object containing
#' the circular kernel PDF, and the activity level estimate derived from this. The
#' generic plot function for \code{actmod} objects plots the distribution. Functions
#' starting with \code{compare} make statistical comparisons between distributions or
#' activity estimates. Note that all time or other circular data should be in radians
#' (in the range 0 to 2*pi).
#'
#' @references Rowcliffe, M., Kays, R., Kranstauber, B., Carbone, C., Jansen, P.A. (2014) Quantifying animal activity level using camera trap data. Methods in Ecology and Evolution.
#' @name activity
"_PACKAGE"

#' Animal record time of day data
#'
#' Barro Colorado Island 2008 data: times of day at which animal records occured
#' (\code{time}), together with species (\code{species}).
#'
#' @format A dataframe with 17820 observations and 2 variables.
#' @source http://dx.doi.org/10.6084/m9.figshare.1160536
#' @name BCItime
#' @docType data
NULL

#' Animal speed data
#'
#' Barro Colorado Island 2008 data: speeds of animal passages past camera traps
#' (\code{speed}), together with species (\code{species}) and time of day (\code{time})
#' for each record.
#'
#' @format A dataframe with 2204 observations and 3 variables.
#' @source http://dx.doi.org/10.6084/m9.figshare.1160536
#' @name BCIspeed
#' @docType data
NULL

#' Activity model class.
#'
#' An S4 class describing activity models fitted to time of observation data.
#'
#' @slot data Object of class \code{"numeric"}, the input data.
#' @slot wt Object of class \code{"numeric"}, weights applied to the data.
#' @slot bw Object of class \code{"numeric"}, kernel bandwidth.
#' @slot adj Object of class \code{"numeric"}, kernel bandwidth adjustment multiplier.
#' @slot pdf Object of class \code{"matrix"} describing fitted probability density function:
#'  Column 1: A regular sequence of radian times at which PDF evaluated; range is [0, 2*pi] if unbounded, and sequence steps are range difference divided by 512.
#'  Column 2: Corresponding circular kernel PDF values.
#' Additionally if errors bootstrapped:
#'  Column 3: PDF standard error.
#'  Column 4: PDF lower 95\% confidence limit. Column 5: PDF upper 95\% confidence limit.
#' @slot act Object of class \code{"numeric"} giving activity level estimate and, if errors boostrapped, standard error and 95 percent confidence limits.
#' @export
setClass("actmod",
         representation(data="numeric", wt="numeric", bw="numeric", adj="numeric",
                        pdf="matrix", act="numeric"))

#' An S4 class describing linear-circular relationships.
#'
#' @slot data Object of class \code{"data.frame"}, the input data, with columns
#' \code{lindat} (linear data) and \code{circdat} (circular data).
#' @slot fit Object of class \code{"data.frame"}, summary of the model fit, with columns:
#'  \code{x}: A regular ascending sequence from 0 to 2*pi at which other columns evaluated;
#'  \code{fit}: The linear fitted values;
#'  \code{p}: The two tailed probability of observing the fitted values under a random (null) circular distribution;
#'  \code{nullLCL}: The lower 95\% confidence limit of the null distribution;
#'  \code{nullUCL}: The upper 95\% confidence limit of the null distribution.
#' @export
setClass("lincircmod", representation(data="data.frame", fit="data.frame"))


#' von Mises density function
#'
#' Probability density function for the von Mises circular distribution.
#'
#' If more than one of x, mu and k have length > 1, values are recycled.
#'
#' @param x numeric angles (assumed to be radian).
#' @param mu numeric, the mean direction of the distribution.
#' @param k non-negative numeric, the concentration parameter distribution (kappa).
#' @param log if TRUE log probabilities are returned.
#' @return Probability density value(s).
#' @examples
#' dvonm(seq(0, 2*pi, len=10), pi, 1)
#' @export
dvonm <- function(x, mu, k, log=FALSE){
  if(any(k<0)) stop("The concentration parameter k must be non-negative")
  x <- x %% (2*pi)
  mu <- mu %% (2*pi)
  res <- (exp(cos(x-mu) - 1))^k / (2 * pi * besselI(k, 0, TRUE))
  if(log) res <- log(res)
  res
}

#' title trigonometric moment length
#'
#' Calculate trigonometric moment length
#'
#' @param x a vector of circular values, assumed to be radian.
#' @param p order of trigonometric moment to be computed.
#' @return Trigonometric moment length of the input.
#' @export
trigmolen <- function(x, p){
  n <- length(x)
  cmean <- atan2(sum(sin(x)), sum(cos(x)))
  sin.p <- sum(sin(p * (x - cmean)))/n
  cos.p <- sum(cos(p * (x - cmean)))/n
  sqrt(sin.p^2 + cos.p^2)
}

#' Index of overlap between circular distributions.
#'
#' Calculates Dhat4 overlap index (see reference) between two kernel distributions.
#'
#' Uses linear interpolation to impute values from kernel distributions.
#'
#' @param fit1,fit2 Fitted activity models of class actmod created using function fitact.
#' @return Scalar overlap index (specifically Dhat4).
#' @references Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337.
#' @examples
#' data(BCItime)
#' oceAct <- fitact(subset(BCItime, species=="ocelot")$time*2*pi)
#' broAct <- fitact(subset(BCItime, species=="brocket")$time*2*pi)
#' ovl4(oceAct, broAct)
#' @export
ovl4 <- function(fit1, fit2){
  f <- stats::approxfun(fit1@pdf[,1], fit1@pdf[,2])
  g <- stats::approxfun(fit2@pdf[,1], fit2@pdf[,2])
  fx <- f(fit1@data)
  gx <- g(fit1@data)
  fy <- f(fit2@data)
  gy <- g(fit2@data)
  xr <- gx/fx
  yr <- fy/gy
  (mean(ifelse(xr>1, 1, xr)) + mean(ifelse(yr>1, 1, yr))) / 2
}

#' Calculate circular kernel bandwidth.
#'
#' Uses an optimisation procedure to calculate the circular kernel bandwidth giving the best fit to the data.
#'
#' Mainly for internal use.
#'
#' @param dat Numeric data vector of radian times.
#' @param K Integer number of values of kappa over which to maximise (see references for details).
#' @return Single numeric bandwidth value.
#' @references Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337.
#' @export
bwcalc <- function(dat,K=3)
{  if(!all(dat>=0 & dat<=2*pi)) warning("some dat values are <0 or >2*pi, expecting radian data")
   if(max(dat)<1) warning("max(dat) < 1, expecting radian data")

   minfunc <- function(kap,k,dat){
     trigmom <- trigmolen(dat, k)
     (besselI(kap,k)/besselI(kap,0) - trigmom)^2
   }
   kapk.calc <- function(k,dat)
     stats::optimise(minfunc,c(0,100),k,dat)$minimum
   kap <- max(sapply(1:K, kapk.calc, dat))
   ((3*length(dat)*kap^2*besselI(2*kap,2)) / (4*pi^0.5*besselI(kap,0)^2))^(2/5)
}

#' Circular kernel probability density function.
#'
#' Optionally weighted Von Mises kernel probability densities.
#'
#' If \code{bw} not provided it is calculated internally using \code{bw.calc}. The \code{adj} argument is used to adjust \code{bw} to facilitate exploration of fit flexibility.
#'
#' @param x Numeric vector of radian times at which to evaluate the PDF.
#' @param dat Numeric vector of radian time data to which the PDF is fitted.
#' @param wt A numeric vector of weights for each \code{dat} value.
#' @param bw Numeric value for kernel bandwidth.
#' @param adj Numeric kernel bandwidth multiplier.
#' @return Numeric vector of probability densities evaluated at \code{x}.
#' @seealso \code{\link{bwcalc}}
#' @examples
#' #Example with made up input
#' tt <- runif(100,0,2*pi)
#' xx <- seq(0,2*pi, pi/256)
#' pdf <- dvmkern(xx, tt)
#' plot(xx, pdf, type="l")
#' @export
dvmkern <- function(x,dat,wt=NULL,bw=NULL,adj=1){
  if(!all(dat>=0 & dat<=2*pi)) warning("some dat values are <0 or >2*pi, expecting radian data")
  if(max(dat)<1) warning("max(dat) < 1, expecting radian data")
  if(!all(x>=0 & x<=2*pi)) warning("some x values are <0 or >2*pi, expecting radian values")
  if(!is.null(wt) & length(wt)!=length(dat)) stop("dat and wt have different lengths")

  if(is.null(bw)) bw <- bwcalc(dat)
  if(is.null(wt)) wt <- rep(1,length(dat))
  dx <- expand.grid(dat,x)
  dif <- abs(dx[,2]-dx[,1])
  i <- dif>pi
  dif[i] <- 2*pi-dif[i]
  prob <- dvonm(dif, 0, bw*adj)
  apply(matrix(prob*wt, nrow=length(dat)),2,sum)/sum(wt)
}

#' Random numbers from empirical distribution function.
#'
#' Random numbers drawn from an empirical distribution defined by paired values and probabilities.
#'
#' @details The distribution function is defined by \code{fit}, which must be a dataframe containing (at least) columns named:
#'  x: a regular sequence of values from which to draw;
#'  y: corresponding pdf values.
#' @param n Integer number of random numbers to return.
#' @param fit Data frame defining the emprical distribution (see details).
#' @return A numeric vector.
#' @examples
#' data(BCItime)
#' tm <- 2*pi*subset(BCItime, species=="paca")$time
#' mod <- fitact(tm)
#' rn <- redf(1000, as.data.frame(mod@pdf))
#' @export
redf <- function(n, fit){
  if(sum(c("x","y") %in% names(fit)) != 2) stop("fit must be a dataframe with (at least) columns named x and y")
  if(diff(range(diff(fit$x)))>0.0001) stop("x doesn't seem to be a regular sequence")

  df <- (fit$y[-1]+fit$y[-nrow(fit)])/2
  cdf <- c(0,cumsum(df)/sum(df))
  rn <- stats::runif(n)
  stats::approx(cdf, fit$x, rn)$y
}

#' Modified kernel density function
#'
#' Modifies \code{stats::density} by:
#'  Adding SE and 95\% confidence intervals for the density to the output; and
#'  Truncating calculation (not just reporting) of density  values on from and/or to.
#'
#' @details Truncation copes with cases where no data are available outside truncation points.
#' Truncation is achieved by fitting the density to the data augmented by reflecting it
#' across each bound using the optimal bandwidth for the unaugmented data, and returning
#' the resulting densities for the region between the bounds.
#'
#' @param x numeric data vector
#' @param reps bootstrap iterations for SE/interval calculation; set to NULL to suppress
#' @param ... Additional arguments passed to \code{stas::density}
#' @return A list with the same components as \code{stats::density} output plus:
#'  \code{se}: standard error of the density
#'  \code{lcl}, \code{ucl}: lower and upper 95\% confidence intervals of the density

#' @examples
#' data(BCItime)
#' tm <- subset(BCItime, species=="ocelot")$time
#' dens <- density2(tm, from=0.25, to=0.75)
#' plot(dens$x, dens$y, type="l")
#' @export
density2 <- function(x, reps=999, ...){
  prm <- list(...)
  prmnm <- names(prm)

  if(sum(c("from","to") %in% prmnm)==2)
    if(prm$from>prm$to) stop("When double-bounded, from must be less than to")

  warn <- FALSE
  dmult <- 1
  xx <- x
  wt <- prm$weights
  if("from" %in% prmnm){
    if(any(x<prm$from)){
      warn <- TRUE
      x <- x[x>prm$from]
      if(!is.null(wt)) prm$weights <- prm$weights[x>prm$from]
    }
    dmult <- dmult+1
    xx <- c(xx, 2*prm$from-x)
    if(!is.null(wt)) wt <- c(wt, prm$weights)
  }
  if("to" %in% prmnm){
    if(any(x>prm$to)){
      warn <- TRUE
      x <- x[x<prm$to]
      if(!is.null(prm$weights)) prm$weights <- prm$weights[x<prm$to]
    }
    dmult <- dmult+1
    xx <- c(xx, 2*prm$to-x)
    if(!is.null(wt)) wt <- c(wt, prm$weights)
    }
  if(warn) warning("Some x values outside bounds were removed")
  if(!is.null(wt)) wt <- wt/sum(wt)
  prm$weights <- wt
  if(!("bw" %in% prmnm & is.numeric(prm$bw))){
    bw <- stats::density(x, ...)$bw
    prm$bw <- bw
  }

  dens <- do.call(stats::density, append(list(x=xx), prm))
  dens$y <- dens$y*dmult
  if(!is.null(reps)){
    f <- function(){
      xs <- sample(xx, length(xx), replace=T)
      denss <- do.call(stats::density, append(list(x=xs), prm))
      dmult*denss$y
    }
    bsres <- replicate(reps, f())
    se <- t(apply(bsres, 1, stats::sd))
    ci <- t(apply(bsres, 1, stats::quantile, c(0.025,0.975)))
    dens <- c(dens, se=se, lcl=list(ci[,1]), ucl=list(ci[,2]))
  }
  dens
}

#' Fit activity model to time-of-day data
#'
#' Fits kernel density to radian time-of-day data and estimates activity level from this distribution.
#' Optionally: 1. bootstraps the distribution, in which case SEs and confidence limits are also
#' stored for activity level and PDF; 2. weights the distribution; 3. truncates the distribution at given times.
#'
#' @details When no \code{bounds} are given (default), a circular kernel distribution is fitted using \code{dvmkern}.
#' Otherwise, a normal kernel distribution is used, truncated at the values of \code{bounds}, using \code{density2}.
#'
#' The bandwidth adjustment multiplier \code{adj} is provided to allow
#' exploration of the effect of adjusting the internally calculated bandwidth on
#' accuracy of activity level estimates.
#'
#' The alternative bootstrapping methods defined by \code{sample} are:
#' \itemize{
#'  \item{\code{"none"}: no bootstrapping}
#'  \item{\code{"data"}: sample from the data}
#'  \item{\code{"model"}: sample from the fitted probability density distribution}
#'  }
#' It's generally better to sample from the data, but sampling from
#' the fitted distribution can sometimes provide more sensible confidence intervals when
#' the number of observations is very small.
#' @param dat A numeric vector of radian time-of-day data.
#' @param wt A numeric vector of weights for each \code{dat} value.
#' @param reps Number of bootstrap iterations to perform. Ignored if \code{sample=="none"}.
#' @param bw Numeric value for kernel bandwidth. If NULL, calculated internally.
#' @param adj Numeric bandwidth adjustment multiplier.
#' @param sample Character string defining sampling method for bootstrapping errors (see details).
#' @param bounds A two-element vector defining radian bounds at which to truncate.
#' @param show Logical whether or not to show a progress bar while bootstrapping.
#' @return An object of type \code{actmod}
#' @examples
#' #Fit without confidence limits
#' data(BCItime)
#' tm <- 2*pi*subset(BCItime, species=="brocket")$time
#' mod1 <- fitact(tm)
#' plot(mod1)
#'
#' #Fit with confidence limits (limited reps to speed up)
#' mod2 <- fitact(tm, sample="data", reps=10)
#' plot(mod2)
#'
#' #Fit weighted function to correct for detection radius 1.2 times higher
#' #by day than by night, assuming day between pi/2 (6 am) and pi*2/3 (6 pm)
#' weight <- 1/ifelse(tm>pi/2 & tm<pi*3/2, 1.2, 1)
#' mod3 <- fitact(tm, wt=weight)
#' plot(mod3)
#' #Overplot unweighted version for comparison
#' plot(mod1, add=TRUE, tline=list(col=2))
#'
#' #Fit truncated function to consider only night time records,
#' #assuming night between pi*3/2 (6 pm) and pi/3 (6 am)
#' mod4 <- fitact(tm, bounds=c(pi*3/2, pi/2))
#' plot(mod4, centre="night")
#' @export
fitact <- function(dat, wt=NULL, reps=999, bw=NULL, adj=1, sample=c("none","data","model"),
                   bounds=NULL, show=TRUE)
{ if(!is.null(wt) & length(wt)!=length(dat)) stop("dat and wt have different lengths")

  sample <- match.arg(sample)
  if(!is.null(wt)){
    if(any(wt<0)) stop("Weights must be non-negative")
    wt <- wt/sum(wt)
  }

  if(is.null(bounds)){ #circular kernel
    if(is.null(bw)) bw <- bwcalc(dat)
    x <- seq(0,2*pi,pi/256)
    pdf <- dvmkern(x, dat, wt, adj, bw)
    act <- 1/(2*pi*max(pdf))
  } else{ #truncated normal kernel
    if(length(bounds)!=2) stop("If provided, bounds must be a two-element vector")
    if(min(bounds)<0 | max(bounds)>2*pi) stop("bounds must be radian (between 0 and 2*pi)")
    bdiff <- diff(bounds)
    if(diff(bounds)<0){
      dat <- dat-ifelse(dat>pi, 2*pi, 0)
      bounds[1] <- bounds[1]-2*pi
      bdiff <- diff(bounds)
    }
    oob <- dat<bounds[1] | dat>bounds[2]
    if(sum(oob>0)){
      warning("Some x values outside bounds were removed")
      dat <- dat[!oob]
      if(!is.null(wt)){
        wt <- wt[!oob]
        wt <- wt/sum(wt)
      }
    }
    dens <- density2(dat, from=bounds[1], to=bounds[2], weights=wt, adjust=adj,
                     bw=if(is.null(bw)) "nrd0" else bw, reps=NULL)
    bw <- dens$bw
    x <- dens$x
    pdf <- dens$y
    act <- 1/(bdiff*max(pdf))
  }

  if(sample=="none")
    sepdf <- lclpdf <- uclpdf <- seact <- lclact <- uclact <- numeric(0) else{
      if(sample=="model")
        samp <- matrix(redf(reps*length(dat), data.frame(x=x,y=pdf)), ncol=reps) else
          samp <- matrix(sample(dat, reps*length(dat), replace=TRUE, prob=wt), ncol=reps)

        if(is.null(bounds)){
          if(show)
            pdfs <- pbapply::pbapply(samp, 2, function(dat) dvmkern(x,dat,wt,adj=adj)) else
              pdfs <- apply(samp, 2, function(dat) dvmkern(x,dat,wt,adj=adj))
        } else
          if(show)
            pdfs <- pbapply::pbapply(samp, 2, function(dat) density2(dat, from=bounds[1], to=bounds[2],
                                     weights=wt, adjust=adj, bw=if(is.null(bw)) "nrd0" else bw,
                                     reps=NULL)$y) else
              pdfs <- apply(samp, 2, function(dat) density2(dat, from=bounds[1], to=bounds[2], weights=wt,
                          adjust=adj, bw=if(is.null(bw)) "nrd0" else bw, reps=NULL)$y)

    sepdf <- apply(pdfs,1,stats::sd)
    lclpdf <- apply(pdfs,1,stats::quantile,probs=0.025)
    uclpdf <- apply(pdfs,1,stats::quantile,probs=0.975)
    if(is.null(bounds)) acts <- 1/(2*pi*apply(pdfs,2,max)) else
      acts <- 1/(bdiff*apply(pdfs,2,max))
    seact <- stats::sd(acts)
    lclact <- stats::quantile(acts,0.025)
    uclact <- stats::quantile(acts,0.975)
  }

  if(is.null(wt)) wt <- 1
  if(min(x)<0){
    dat <- dat+ifelse(dat<0, 2*pi, 0)
    x <- x+ifelse(x<0, 2*pi, 0)
  }
  pdftab <- cbind(x=x, y=pdf, se=sepdf, lcl=lclpdf, ucl=uclpdf)[order(x), ]

  methods::new("actmod", data=dat, wt=wt, bw=bw, adj=adj, pdf=pdftab,
      act=c(act=act, se=seact, lcl=lclact, ucl=uclact))
}

#' Compare circular distributions.
#'
#' Randomisation test for the probability that two sets of circular observations come from the same distribution.
#'
#' Calculates overlap index Dhat4 (see references) for the two fitted distributions, then generates a null distribution of overlap indices using data sampled randomly with replacement from the combined data.
#' This randomised distribution is then used to define an empirical probability distribution against which  the probability that the observed overlap arose by chance is judged.
#' When one or both fitted models use weighted distributions, sampling probabilities are taken from the weights. If both models are weighted, the weights must therefore be on the same scale.
#'
#' @param fit1,fit2 Fitted activity models of class actmod created using function fitact.
#' @param reps Number of bootstrap iterations.
#' @return A named 4-element vector: obs = observed overlap index; null = mean null overlap index; seNull = standard error of the null distribution; pNull = probability observed index arose by chance.
#' @references Ridout, M.S. & Linkie, M. (2009) Estimating overlap of daily activity patterns from camera trap data. Journal of Agricultural Biological and Environmental Statistics, 14, 322-337.
#' @examples
#' #Example with bootstrap reps limited to reduce run time
#' data(BCItime)
#' tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"]
#' tRat <- 2*pi*BCItime$time[BCItime$species=="rat"]
#' fPaca <- fitact(tPaca)
#' fRat <- fitact(tRat)
#' compareCkern(fPaca,fRat,reps=10)
#' @export
compareCkern <- function(fit1, fit2, reps=999){
  if(!inherits(fit1, "actmod") | !inherits(fit2, "actmod"))
    stop("Input must be fitted activity models (class actmod)")
  bnd <- range(fit1@pdf[,1])
  if(!all(bnd==range(fit2@pdf[,1])))
    stop("Distribution bounds are not identical")

  if(diff(bnd)==2*pi) bnd <- NULL
  olp <- ovl4(fit1,fit2)
  y1 <- fit1@data
  y2 <- fit2@data
  w1 <- fit1@wt
  w2 <- fit2@wt
  if(length(w1)==1) w1 <- rep(1, length(y1))
  if(length(w2)==1) w2 <- rep(1, length(y2))
  y <- c(y1,y2)
  w <- c(w1,w2)
  samp <- matrix(sample(1:length(y), reps*length(y), replace=T, prob=w), nrow=reps)
  f <- function(s){
    m1 <- fitact(y[s[1:length(y1)]], sample="n", bw=fit1@bw, adj=fit1@adj, bounds=bnd)
    m2 <- fitact(y[s[(length(y1)+1):length(y)]], sample="n", bw=fit1@bw, adj=fit1@adj, bounds=bnd)
    ovl4(m1,m2)
  }
  res <- pbapply::pbapply(samp, 1, f)
  fun <- stats::ecdf(res)
  c(obs=olp, null=mean(res), seNull=stats::sd(res), pNull=fun(olp))
}


#' Compare activity level estimates
#'
#' Wald test for the statistical difference between two or more activitiy level estimates.
#'
#' Uses a Wald test to ask whether the difference between estimates a1 and a2 is
#' significantly different from 0: statistic W = (a1-a2)^2 / (SE1^2+SE2^2) tested
#' on chi-sq distribution with 1 degree of freedom.
#'
#' @param fits A list of fitted \code{actmod} objects
#' @return A matrix with 4 columns: 1. differences between estimates; 2. SEs of the differences; 3. Wald statistics; 4. p-values (H0 is no difference between estimates). Matrix rows give all possible pairwise comparisons, numbered in the order in which they entered in the list \code{fits}.
#' @examples
#' #Test whether paca have a sigificantly different activity level from rat.
#' #Bootstrap reps limited to speed up example.
#' data(BCItime)
#' tPaca <- 2*pi*BCItime$time[BCItime$species=="ocelot"]
#' tRat <- 2*pi*BCItime$time[BCItime$species=="rat"]
#' fPaca <- fitact(tPaca, sample="data", reps=10)
#' fRat <- fitact(tRat, sample="data", reps=10)
#' fPaca@act
#' fRat@act
#' compareAct(list(fPaca,fRat))
#' @export
compareAct <- function(fits)
{ if(!inherits(fits, "list") | !all(unlist(lapply(fits,inherits,"actmod"))))
    stop("fits must be a list of actmod objects")
  if(min(unlist(lapply(fits, function(x) length(x@act))))==1)
    stop("all input model fits must be boostrapped")

  len <- length(fits)
  i <- rep(1:(len-1), (len-1):1)
  j <- unlist(sapply(2:len, function(i) i:len))
  acts <- unlist(lapply(fits, function(fit) fit@act[1]))
  seacts <- unlist(lapply(fits, function(fit) fit@act[2]))
  dif <- acts[i]-acts[j]
  vardif <- seacts[i]^2 + seacts[j]^2
  W <- dif^2/vardif
  prob <- 1-stats::pchisq(W,1)
  res <- cbind(Difference=dif, SE=sqrt(vardif), W=W, p=prob)
  dimnames(res)[[1]] <- paste(i,j,sep="v")
  res
}

#' Compare activity between times of day
#'
#' Uses a Wald test to statistically compare activity levels at given radian times of day for a fitted activity distribution.
#'
#' Bootrapping the activity model yields standard error estimates for the PDF. This function uses these SEs to compute a Wald statistic for the difference between PDF values (by inference activity levels) at given times of day: statistic W = (a1-a2)^2 / (SE1^2+SE2^2) tested on chi-sq distribution with 1 degree of freedom.
#'
#' @param fit Fitted \code{actmod} object with errors boostrapped (fit using \code{fitact} with \code{sample} argument != "none").
#' @param times Numeric vector of radian times of day at which to compare activity levels. All pairwise comparisons are made.
#' @return A matrix with 4 columns: 1. differences between PDF values; 2. SEs of the differences; 3. Wald statistics; 4. p-values (H0 is no difference between estimates). Matrix rows give all possible pairwise comparisons, numbered in the order in which they appear in vector \code{times}.
#' @examples
#' data(BCItime)
#' tPaca <- 2*pi*BCItime$time[BCItime$species=="paca"]
#' fPaca <- fitact(tPaca, sample="data", reps=10)
#' plot(fPaca)
#' compareTimes(fPaca, c(5.5,6,0.5,1))
#' @export
compareTimes <- function(fit, times)
{ if(!inherits(fit, "actmod")) stop("fit input must be an actmod object")
  if(!all(times>=0 & times<=2*pi)) stop("some times are <0 or >2*pi, expecting radian data")

  len <- length(times)
  i <- rep(1:(len-1), (len-1):1)
  j <- unlist(sapply(2:len, function(i) i:len))
  k <- findInterval(times, fit@pdf[,1])
  p <- (times-fit@pdf[,1][k]) / (fit@pdf[,1][k+1]-fit@pdf[,1][k])
  pdfs1 <- fit@pdf[,2][k]
  pdfs2 <- fit@pdf[,2][k+1]
  pdfs <- pdfs1 + p*(pdfs2-pdfs1)
  sepdfs1 <- fit@pdf[,3][k]
  sepdfs2 <- fit@pdf[,3][k+1]
  sepdfs <- sepdfs1 + p*(sepdfs2-sepdfs1)
  dif <- pdfs[i]-pdfs[j]
  vardif <- sepdfs[i]^2 + sepdfs[j]^2
  W <- dif^2/vardif
  prob <- 1-stats::pchisq(W,1)
  res <- cbind(Difference=dif, SE=sqrt(vardif), W=W, p=prob)
  dimnames(res)[[1]] <- paste(i,j,sep="v")
  res
}

#' Linear-circular kernel fit
#'
#' Fits a Von Mises kernel distribution describing a linear variable as a function of a circular predictor.
#'
#' @param x Numeric vector of radian values at which to evaluate the distribution.
#' @param circdat Numeric vector of radian data matched with \code{lindat}.
#' @param lindat Numeric vector of linear data matched with \code{circdat}.
#' @return A numeric vector of fitted \code{lindat} values matched with \code{x}.
#' @references Xu, H., Nichols, K. & Schoenberg, F.P. (2011) Directional kernel regression for wind and fire data. Forest Science, 57, 343-352.
#' @examples
#' data(BCIspeed)
#' i <- BCIspeed$species=="ocelot"
#' log_speed <- log(BCIspeed$speed[i])
#' time <- BCIspeed$time[i]*2*pi
#' circseq <- seq(0,2*pi,pi/256)
#' trend <- lincircKern(circseq, time, log_speed)
#' plot(time, log_speed, xlim=c(0, 2*pi))
#' lines(circseq, trend)
#' @export
lincircKern <- function(x,circdat,lindat)
{ if(length(lindat)!=length(circdat))
  stop("lindat and circdat lengths are unequal")
  if(min(circdat)<0 | max(circdat)>2*pi)
    stop("circdat values not between 0 and 2*pi, expecting radian data")
  if(min(x)<0 | max(x)>2*pi)
    stop("x values not between 0 and 2*pi, expecting radian values")

  hs <- 1.06 * min(stats::sd(lindat), (stats::quantile(lindat,0.75)-stats::quantile(lindat,0.25))/1.34) *
    length(lindat)^-0.2
  bw <- 1/hs^2
  dx <- expand.grid(circdat,x)
  dif <- abs(dx[,2]-dx[,1])
  i <- dif>pi
  dif[i] <- 2*pi-dif[i]
  prob <- matrix(dvonm(dif, 0, bw), nrow=length(circdat))
  apply(prob,2,function(z) mean(z*lindat)/mean(z))
}

#' Linear-circular regression
#'
#' Fits a Von Mises kernel distribution describing a linear variable as a function
#' of a circular predictor, and boostraps the null distribution in order to evaluate
#' significance of radial variation in the linear variable.
#'
#' Deviation of \code{lindat} from the null expecation is assessed either visually
#' by the degree to which the fitted distribution departs from the null confidence
#' interval (use generic plot function), or quantitatively by column \code{p} of
#' slot \code{fit} in the resulting \code{lincircmod-class} object.
#'
#' @param circdat Numeric vector of radian data matched with \code{lindat}.
#' @param lindat Numeric vector of linear data matched with \code{circdat}.
#' @param pCI Single numeric value between 0 and 1 defining proportional confidence interval to return.
#' @param reps Integer number of bootstrap repetitions to perform.
#' @param res Resolution of fitted distribution and null confidence interval - specifically a single integer number of points on the circular scale at which to record distributions.
#' @return An object of type \code{\link{lincircmod-class}}
#' @references Xu, H., Nichols, K. & Schoenberg, F.P. (2011) Directional kernel regression for wind and fire data. Forest Science, 57, 343-352.
#' @examples
#' #Example with reps limited to increase speed
#' data(BCIspeed)
#' i <- BCIspeed$species=="ocelot"
#' sp <- log(BCIspeed$speed[i])
#' tm <- BCIspeed$time[i]*2*pi
#' mod <- fitlincirc(tm, sp, reps=50)
#' plot(mod, CircScale=24, xaxp=c(0,24,4), xlab="Time", ylab="log(speed)")
#' legend(8,-3, c("Fitted speed", "Null CI"), col=1:2, lty=1:2)
#' @export
fitlincirc <- function(circdat, lindat, pCI=0.95, reps=10, res=512)
{ if(length(lindat)!=length(circdat))
  stop("lindat and circdat lengths are unequal")
  if(min(circdat)<0 | max(circdat)>2*pi)
    stop("circdat values not between 0 and 2*pi, expecting radian data")
  if(max(circdat)<1)
    warning("max(circdat) < 1, expecting radian data")

  n <- length(circdat)
  x <- seq(0,2*pi,2*pi/res)

  bs <- pbapply::pbsapply(1:reps, function(i)
  {  j <- sample(1:n,n,TRUE)
     lincircKern(x,circdat[j],lindat)
  })
  nulllcl <- apply(bs,1,stats::quantile,(1-pCI)/2)
  nullucl <- apply(bs,1,stats::quantile,(1+pCI)/2)
  fit <- lincircKern(x,circdat,lindat)
  p <- sapply(1:(res+1), function(i)
  { f <- stats::ecdf(bs[i,])
    f(fit[i])
  })
  p[p>0.5] <- 1-p[p>0.5]
  p <- 2*p

  methods::new("lincircmod", data=data.frame(circdat=circdat, lindat=lindat),
      fit=data.frame(x=x, fit=fit, p=p, nullLCL=nulllcl, nullUCL=nullucl))
}

#' Plot activity distribution
#'
#' Plot an activity probability distribution from a fitted \code{actmod} object.
#'
#' @details When xunit=="clock", The underlying numeric range of the x-axis is [0,24] if centre=="day",
#' or [-12,12] if centre=="night".
#' @param x Object of class \code{actmod}.
#' @param xunit Character string defining x-axis unit.
#' @param yunit Character string defining y-axis unit.
#' @param data Character string defining whether to plot the data distribution and if so which style to use.
#' @param centre Character string defining whether to centre the plot on midday or midnight.
#' @param dline List of plotting parameters for data lines.
#' @param tline List of plotting parameters for trend line.
#' @param cline List of plotting parameters for trend confidence interval lines.
#' @param add Logical defining whether to create a new plot (default) or add to an existing plot.
#' @param xaxis List of plotting parameters to pass to axis command for x-axis plot (see axis for arguments).
#' @param ... Additional arguments passed to internal plot call affecting only the plot frame and y axis. Modify x axis through xaxis.
#' @return No return value, called to create a plot visualising an activity model.
#' @examples
#' data(BCItime)
#' otm <- 2*pi*subset(BCItime, species=="ocelot")$time
#' btm <- 2*pi*subset(BCItime, species=="brocket")$time
#' omod <- fitact(otm)
#' bmod <- fitact(btm)
#' plot(omod, yunit="density", data="none")
#' plot(bmod, yunit="density", data="none", add=TRUE, tline=list(col="red"))
#' legend("topleft", c("Ocelot", "Brocket deer"), col=1:2, lty=1)
#'
#' mod <- fitact(otm, sample="data", reps=10)
#' plot(mod, dline=list(col="grey"),
#'           tline=list(col="red", lwd=2),
#'           cline=list(col="red", lty=3))
#'
#' mod2 <- fitact(otm, bounds=c(pi*3/2, pi/2))
#' plot(mod2, centre="night")
#' plot(mod2, centre="night", xlim=c(-6,6), xaxis=list(at=seq(-6,6,2)))
#' @export
plot.actmod <- function(x, xunit=c("clock","hours","radians"), yunit=c("frequency","density"),
                        data=c("histogram","rug","both","none"), centre=c("day","night"),
                        dline=list(lwd=ifelse(data=="rug", 0.1, 1)), tline=NULL, cline=list(lty=2),
                        add=FALSE, xaxis=NULL, ...){

  #function sets up the plot parameters with appropriate scales and labels with flexibility for modifcation using ...
  setup <- function(...){
    pprm <- list(...)
    pprm <- append(list(x=xlim, y=ylim, type="n", xaxt="n"), pprm)
    if(!("ylim" %in% names(pprm))) pprm <- append(pprm, list(ylim=ylim))
    if(!("ylab" %in% names(pprm)))
      pprm <- switch(yunit,
                     "frequency"=append(pprm, list(ylab="Frequency")),
                     "density"=append(pprm, list(ylab="Density"))
      )
    if(!("xlab" %in% names(pprm)))
      pprm <- switch(xunit,
                     "clock"=append(pprm, list(xlab="Time")),
                     "hours"=append(pprm, list(xlab="Hours")),
                     "radians"=append(pprm, list(xlab="Radians"))
      )
    do.call(graphics::plot, pprm)
    xaxis <- append(list(side=1), xaxis[!(names(xaxis)=="side")])
    if(xunit=="clock"){
      if(!("at" %in% names(xaxis))){
        at <- switch(centre,
                     "day"=seq(0, 24, 6),
                     "night"=seq(-12, 12, 6)
        )
        xaxis <- append(xaxis, list(at=at))
      }
      if(!("labels" %in% names(xaxis))){
        at <- xaxis$at
        hh <- trunc(at)
        mm <- round(60*(at-trunc(at)), 0)
        hh <- hh + ifelse(hh<0,24,0) - ifelse(mm<0, 1, 0)
        mm <- mm + ifelse(mm<0, 60, 0)
        lab <- paste0(ifelse(hh<10,"0",""), hh, ":", ifelse(mm<10,"0",""), mm)
        xaxis <- append(xaxis, list(labels=lab))
      }
    }
    do.call(graphics::axis, xaxis)
  }
  ### end setup function ###

  data <- match.arg(data)
  xunit <- match.arg(xunit)
  yunit <- match.arg(yunit)
  centre <- match.arg(centre)
  if(!"lty" %in% names(cline)) cline <- append(cline, list(lty=2))

  fit <- x
  pdf <- fit@pdf
  fdata <- fit@data

  if(centre=="night"){
    if(diff(range(pdf[,1]))==2*pi | min(pdf[,1]>=0)){
      if(diff(range(pdf[,1]))==2*pi) pdf <- pdf[-1,]
      pdf[,1] <- pdf[,1] - ifelse(pdf[,1]>pi, 2*pi, 0)
      pdf <- pdf[order(pdf[,1]), ]
    }
    fdata <- fdata-ifelse(fdata>pi,2*pi,0)
    xlim <- c(-pi,pi)
  }else xlim <- c(0,2*pi)

  x <- pdf[,"x"]
  y <- pdf[,"y"]
  if(ncol(pdf)==5)
  { lcl <- pdf[,"lcl"]
  ucl <- pdf[,"ucl"]
  } else lcl <- ucl <- numeric(0)

  if(xunit=="radians") maxbrk <- 2*pi else {
    xlim <- xlim*12/pi
    x <- x*12/pi
    fdata <- fdata*12/pi
    maxbrk <- 24
  }

  if(yunit=="frequency"){
    y <- y*length(fdata)*pi/12
    lcl <- lcl*length(fdata)*pi/12
    ucl <- ucl*length(fdata)*pi/12
  }else{
    if(xunit %in% c("clock","hours")){
      y <- y*pi/12
      lcl <- lcl*pi/12
      ucl <- ucl*pi/12
    }
  }

  if(data %in% c("histogram","both")){
    h <- switch(centre,
                "day"=graphics::hist(fdata, breaks=seq(0,maxbrk,maxbrk/24), plot=F),
                "night"=graphics::hist(fdata, breaks=seq(-maxbrk/2,maxbrk/2,maxbrk/24), plot=F))
    d <- switch(yunit,
                "frequency"=h$counts,
                "density"=h$density)
    ylim <- c(0,max(y,d,ucl,na.rm=T))
  } else ylim <- c(0,max(y,ucl,na.rm=T))

  #Plot data
  if(!add) setup(...)
  if(data %in% c("histogram","both")){
    if(!"lwd" %in% names(dline)) dline <- append(dline, list(lwd=1))
    do.call(graphics::lines, append(list(x=h$breaks, y=c(d,d[1]), type="s"), dline))
  }
  if(data %in% c("rug","both")){
    if(!"lwd" %in% names(dline)) dline <- append(dline, list(lwd=0.1))
    for(i in 1:length(fdata))
      do.call(graphics::lines, append(list(x=rep(fdata[i],2), y=max(y,na.rm=T)*-c(0,0.03)), dline))
  }

  #Plot trend
  i <- x < switch(centre, "day"=ifelse(xunit=="radians", pi, 12), 0)
  do.call(graphics::lines, append(list(x=x[i], y=y[i]), tline))
  do.call(graphics::lines, append(list(x=x[!i], y=y[!i]), tline))

  #Plot conf intervals
  if(length(lcl)>0)
  { do.call(graphics::lines, append(list(x=x[i], y=lcl[i]), cline))
    do.call(graphics::lines, append(list(x=x[i], y=ucl[i]), cline))
    do.call(graphics::lines, append(list(x=x[!i], y=lcl[!i]), cline))
    do.call(graphics::lines, append(list(x=x[!i], y=ucl[!i]), cline))
  }
}

#' Plot linear-circular relationship
#'
#' Plot linear against circular data along with the fitted and null confidence limit distributions from a fitted \code{lincircmod} object.
#'
#' @param x Object of class \code{lincircmod}.
#' @param CircScale Single numeric value defining the plotting maximum of the circular scale.
#' @param tlim Numeric vector with two elements >=0 and <=1 defining the lower and upper limits at which to plot distributions; default plots the full range.
#' @param fcol,flty,ncol,nlty Define line colour (\code{col}) and type (\code{lty}) for fitted (\code{f}) and null (\code{n}) distributions; input types as for \code{col} and \code{lty}, see \code{\link{par}}.
#' @param ... Additional arguments passed to the inital plot construction, affecting axes and data plot symbols.
#' @return No return value, called to create a plot visualising a linear-circular relationship.
#' @export
plot.lincircmod <- function(x, CircScale=2*pi, tlim=c(0,1), fcol="black", flty=1, ncol="red", nlty=2, ...){

  if(min(tlim)<0 | max(tlim)>1 | length(tlim)!=2) stop("tlim should contain two values >=0 and <=1")

  fit <- x
  fit <- x@fit
  dat <- x@data
  xx <- fit$x*CircScale/(2*pi)
  LinearData <- dat$lindat
  CircularData <- dat$circdat*CircScale/(2*pi)
  range <- tlim*CircScale
  graphics::plot(CircularData, LinearData, ...)
  if(range[1]<range[2])
  { i <- xx>=range[1] & xx<=range[2]
    graphics::lines(xx[i], fit$fit[i], col=fcol, lty=flty)
    graphics::lines(xx[i], fit$nullLCL[i], col=ncol, lty=nlty)
    graphics::lines(xx[i], fit$nullUCL[i], col=ncol, lty=nlty)
  } else
  {  i <- xx>=range[1]
     graphics::lines(xx[i], fit$fit[i], col=fcol, lty=flty)
     graphics::lines(xx[i], fit$nullLCL[i], col=ncol, lty=nlty)
     graphics::lines(xx[i], fit$nullUCL[i], col=ncol, lty=nlty)
     i <- xx<=range[2]
     graphics::lines(xx[i], fit$fit[i], col=fcol, lty=flty)
     graphics::lines(xx[i], fit$nullLCL[i], col=ncol, lty=nlty)
     graphics::lines(xx[i], fit$nullUCL[i], col=ncol, lty=nlty)
  }
}


#' Convert time of day data to numeric
#'
#' Accepts data of class POSIXct, POSIXlt or character and returns the  time of day element as numeric (any date element is ignored).
#'
#' @param x A vector of POSIXct, POSIXlt or character format time data to convert.
#' @param scale The scale on which to return times (see Value for options).
#' @param ... arguments passed to as.POSIXlt
#' @param tryFormats formats to try when converting date from character, passed to as.POSIXlt
#' @return A vector of numeric times of day in units defined by the \code{scale} argument:
#' radian, on the range [0, 2*pi];
#' hours, on the range [0, 24];
#' proportion, on the range [0, 1].
#' @seealso \code{\link{strptime}}
#' @examples
#' data(BCItime)
#' rtime <- gettime(BCItime$date)
#' htime <- gettime(BCItime$date, "hour")
#' ptime <- gettime(BCItime$date, "proportion")
#' summary(rtime)
#' summary(htime)
#' summary(ptime)
#' @export
gettime <- function(x, scale=c("radian","hour","proportion"), ...,
                    tryFormats=c("%Y-%m-%d %H:%M:%OS",
                                 "%Y/%m/%d %H:%M:%OS",
                                 "%Y:%m:%d %H:%M:%OS",
                                 "%Y-%m-%d %H:%M",
                                 "%Y/%m/%d %H:%M",
                                 "%Y:%m:%d %H:%M",
                                 "%Y-%m-%d",
                                 "%Y/%m/%d",
                                 "%Y:%m:%d")){
  scale <- match.arg(scale)
  x <- as.POSIXlt(x, tryFormats=tryFormats, ...)
  res <- x$hour + x$min/60 + x$sec/3600
  if(scale=="radian") res <- res*pi/12
  if(scale=="proportion") res <- res/24
  if(all(res==0, na.rm=T)) warning("All times are 0: may be just strptime default?")
  res
}


#' Wraps data on a given range.
#'
#' Input data outside the given bounds (default radian [0, 2*pi]) are wrapped to appear within the range.
#'
#' @details As an example of wrapping, on bounds [0, 1], a value of 1.2 will be converted to 0.2, while a value of -0.2 will be converted to 0.8.
#' @param x A vector of numeric data.
#' @param bounds The range within which to wrap \code{x} values
#' @return A vector of numeric values within the limits defined by \code{bounds}
#' @examples
#' data(BCItime)
#' adjtime <- BCItime$time + 1/24
#' summary(adjtime)
#' adjtime <- wrap(adjtime, c(0,1))
#' summary(adjtime)
#' @export
wrap <- function(x, bounds=c(0,2*pi)){
  bounds[1] + (x-bounds[1]) %% diff(bounds)
}


#' Circular mean
#'
#' Calculates the average direction of a set of radian circular values.
#'
#' @details The \code{base::mean} function is use internally, and additional arguments, e.g for missing data handling, are passed to this.
#' @param x A vector of radian values.
#' @param ... Arguments passed to \code{mean}.
#' @return A radian value giving mean direction.
#' @seealso \code{\link{mean}}
#' @examples
#' data(BCItime)
#' times <- subset(BCItime, species=="ocelot")$time*2*pi
#' cmean(times)
#' @export
cmean <- function(x, ...){
  X <- mean(cos(x), ...)
  Y <- mean(sin(x), ...)
  wrap(atan(Y/X) + ifelse(X<0,pi,0))
}

#' Calculates solar event times
#'
#' Calculates approximate times of sunrise and sunset and day lengths
#' for given dates at given locations.
#'
#' @details Function adapted from https://www.r-bloggers.com/2014/09/seeing-the-daylight-with-r/
#' @param date character, POSIX or Date format date/time value(s)
#' @param lat,lon latitude and longitude in decimal degrees
#' @param offset the time offset in hours relative to UTC (GMT) for results
#' @param ... arguments passed to as.POSIXlt
#' @param tryFormats formats to try when converting date from character, passed to as.POSIXlt
#' @return A dataframe with columns sunrise and sunset (given in the timezone defined by offset) and daylength, all expressed in hours.
#' @references Teets, D.A. 2003. Predicting sunrise and sunset times. The College Mathematics Journal 34(4):317-321.
#' @examples
#' data(BCItime)
#' dat <- subset(BCItime, species=="ocelot")$date
#' get_suntimes(dat, 9.156335, -79.847682, -5)
#' @export
get_suntimes <- function(date, lat, lon, offset, ...,
                         tryFormats=c("%Y-%m-%d %H:%M:%OS",
                                      "%Y/%m/%d %H:%M:%OS",
                                      "%Y:%m:%d %H:%M:%OS",
                                      "%Y-%m-%d %H:%M",
                                      "%Y/%m/%d %H:%M",
                                      "%Y:%m:%d %H:%M",
                                      "%Y-%m-%d",
                                      "%Y/%m/%d",
                                      "%Y:%m:%d")){
  nlat <- length(lat)
  nlon <- length(lon)
  ndat <- length(date)
  if((nlat>1 & nlat!=ndat) | (nlon>1 & nlon!=ndat))
    stop("lat and lon must have length 1 or the same length as date")
  # Day of year
  d <- as.POSIXlt(date, tryFormats=tryFormats, ...)$yday + 1
  # Radius of the earth (km)
  R <- 6378
  # Radians between the xy-plane and the ecliptic plane
  epsilon <- 23.45 * pi / 180
  # Convert observer's latitude to radians
  L <- lat * pi / 180
  # Calculate offset of sunrise based on longitude (min)
  # If lon is negative, then the mod represents degrees West of
  # a standard time meridian, so timing of sunrise and sunset should
  # be made later.
  lon_h <- 24 * lon / 360
  # The earth's mean distance from the sun (km)
  r <- 149598000
  theta <- 2 * pi / 365.25 * (d - 80)
  zs <- r * sin(theta) * sin(epsilon)
  rp <- sqrt(r^2 - zs^2)
  acos_term <- (R-zs * sin(L)) / (rp * cos(L))
  t0 <- suppressWarnings(1440 / (2*pi) * acos(acos_term))
  # A kludge adjustment for the radius of the sun
  that <- t0+5
  # Adjust "noon" for the fact that the earth's orbit is not circular:
  n <- 720 - 10 * sin(4 * pi * (d-80) / 365.25) + 8 * sin(2*pi*d / 365.25)
  ## now sunrise and sunset are:
  sunrise <- (n - that) / 60 - lon_h + offset
  sunset <- (n + that) / 60 - lon_h + offset
  daylength <- ifelse(acos_term < -1, 24,
                      ifelse(acos_term > 1, 0,
                             sunset-sunrise))
  data.frame(sunrise=sunrise, sunset=sunset, daylength=daylength)
}


#' Transforms clock time to solar time anchored to sun rise and sunset times for a given location.
#'
#' This is a wrapper for \code{transtime} that takes non-numeric date-time input together with latitude and longitude to calculate mean average sunrise and sunset times, which are then used to anchor the transformation using average anchoring.
#'
#' @details Time zone \code{tz} should be expressed in numeric hours relative to UTC (GMT).
#' @param dat A vector of character, POSIXct or POSIXlt date-time values.
#' @param lat,lon Single numeric values or numeric vectors the same length as \code{dat} giving site latitude and longitude in decimal format.
#' @param tz A single numeric value or numeric vector same length as \code{dat} giving time zone (see Details).
#' @param ... arguments passed to as.POSIXlt
#' @param tryFormats formats to try when converting date from character, passed to as.POSIXlt
#' @return A list with elements:
#' @return \code{input}: event input dates-times in POSIXlt format.
#' @return \code{clock}: radian clock time data.
#' @return \code{solar}: radian solar time data anchored to average sun rise and sun set times.
#' @references Vazquez, C., Rowcliffe, J.M., Spoelstra, K. and Jansen, P.A. in press. Comparing diel activity patterns of wildlife across latitudes and seasons: time transformation using day length. Methods in Ecology and Evolution.
#' @seealso \code{\link{strptime}, \link{transtime}}
#' @examples
#' data(BCItime)
#' subdat <- subset(BCItime, species=="ocelot")
#' times <- solartime(subdat$date, 9.156335, -79.847682, -5)
#' rawAct <- fitact(times$clock)
#' avgAct <- fitact(times$solar)
#' plot(rawAct)
#' plot(avgAct, add=TRUE, data="n", tline=list(col="cyan"))
#' @export
solartime <- function(dat, lat, lon, tz, ...,
                      tryFormats=c("%Y-%m-%d %H:%M:%OS",
                                   "%Y/%m/%d %H:%M:%OS",
                                   "%Y:%m:%d %H:%M:%OS",
                                   "%Y-%m-%d %H:%M",
                                   "%Y/%m/%d %H:%M",
                                   "%Y:%m:%d %H:%M",
                                   "%Y-%m-%d",
                                   "%Y/%m/%d",
                                   "%Y:%m:%d")){
  dat <- as.POSIXlt(dat, tryFormats=tryFormats, ...)
  posdat <- list(lat,lon,tz)
  if(!all(unlist(lapply(posdat, inherits, "numeric"))) |
     any(unlist(lapply(posdat, length)) != length(dat) &
         unlist(lapply(posdat, length)) != 1))
    stop("lat, lon and tz must all be numeric scalars, or vectors the same length as dat")

  suntimes <- wrap(get_suntimes(dat, lat, lon, tz)[,-3] * pi/12)
  tm <- gettime(dat)
  list(input=dat, clock=tm, solar=transtime(tm, suntimes))
}


#' Transforms clock time to solar times.
#'
#' Transforms time expressed relative to either the time of a single solar event (anchor times - Nouvellet et al. 2012), or two solar events (such as sun rise and sun set - Vazquez et al. in press).
#'
#' @details If double anchoring is requested (i.e. \code{type} is equinoctial
#' or average), the \code{anchor} argument requires a two-column matrix,
#' otherwise a vector. The argument \code{mnanchor} can usually be left at
#' its default \code{NULL} value. In this case, the mean anchors are set to
#' \code{c(pi/2, pi*3/2)} when \code{type}=="equinoctial", otherwise the
#' \code{anchor} mean(s).
#'
#' Although the anchors for transformation are usually likely to be solar
#' events (e.g. sun rise and/or sunset), they could be other celestial
#' (e.g. lunar) or human-related (e.g. timing of artificial lighting) events.
#' @param dat A vector of radian event clock times.
#' @param anchor A vector or matrix matched with \code{dat} containing radian anchor times on the day of each event (see Details).
#' @param mnanchor A scalar or two-element vector of numeric radian mean anchor times (see Details).
#' @param type The type of transformation to use (see Details).
#' @return  A vector of radian transformed times.
#' @references Vazquez, C., Rowcliffe, J.M., Spoelstra, K. and Jansen, P.A. in press. Comparing diel activity patterns of wildlife across latitudes and seasons: time transformation using day length. Methods in Ecology and Evolution.
#' @references Nouvellet, P., Rasmussen, G.S.A., Macdonald, D.W. and Courchamp, F. 2012. Noisy clocks and silent sunrises: measurement methods of daily activity pattern. Journal of Zoology 286: 179-184.
#' @examples
#' data(BCItime)
#' subdat <- subset(BCItime, species=="ocelot")
#' suntimes <- pi/12 * get_suntimes(subdat$date, 9.156335, -79.847682, -5)[, -3]
#' rawtimes <- subdat$time*2*pi
#' avgtimes <- transtime(rawtimes, suntimes)
#' eqntimes <- transtime(rawtimes, suntimes, type="equinoctial")
#' sngtimes <- transtime(rawtimes, suntimes[,1], type="single")
#' rawAct <- fitact(rawtimes)
#' avgAct <- fitact(avgtimes)
#' eqnAct <- fitact(eqntimes)
#' sngAct <- fitact(sngtimes)
#' plot(rawAct)
#' plot(avgAct, add=TRUE, data="n", tline=list(col="magenta"))
#' plot(eqnAct, add=TRUE, data="n", tline=list(col="orange"))
#' plot(sngAct, add=TRUE, data="n", tline=list(col="cyan"))
#' @export
transtime <- function(dat, anchor, mnanchor=NULL, type=c("average", "equinoctial", "single")){
  if(!all(dat>=0 & dat<=2*pi, na.rm=TRUE)) warning("some dat values are <0 or >2*pi, expecting radian data")
  if(max(dat, na.rm=TRUE)<1) warning("max(dat) < 1, expecting radian data")
  if(!all(anchor>=0 & anchor<=2*pi, na.rm=TRUE)) warning("some anchor values are <0 or >2*pi, expecting radian values")
  if(is.null(ncol(anchor))) anchor <- matrix(anchor, ncol=1)
  if(!all(apply(anchor, 2, is.numeric))) stop("anchor must be a numeric vector, matrix or data.frame")
  nr <- nrow(anchor)
  if(length(dat) != nr) stop("dat and anchor have different lengths")
  type <- match.arg(type)
  nc <- ncol(anchor)
  if(is.null(mnanchor)) mnanchor <- apply(anchor, 2, cmean, na.rm=TRUE)
  if(type=="single"){
    if(nc>1) warning("only one column needed for anchor; additional columns ignored")
  } else{
    if(nc==1) stop("double anchoring requires a two-column matrix or data.frame for anchor")
    if(!is.vector(mnanchor) | length(mnanchor)!=2) stop("if provided, mnanchor must be a 2-element vector for double anchoring")
    if(nc>2) warning("only two columns needed for anchor; additional columns ignored")
  }

  if(type=="single"){
    res <- wrap(mnanchor[1] + dat - anchor[,1])
  } else{
    difs <- wrap(cbind(dat,dat)-anchor)
    flip <- difs[,1]>difs[,2]
    a1 <- ifelse(flip, anchor[,2], anchor[,1])
    a2 <- ifelse(flip, anchor[,1], anchor[,2])
    relpos <- wrap(dat-a1) / wrap(a2-a1)
    interval <- switch(type,
                       "equinoctial"=pi,
                       "average"=ifelse(flip, wrap(mnanchor[1]-mnanchor[2]), wrap(mnanchor[2]-mnanchor[1]))
    )
    baseline <- switch(type,
                       "equinoctial"=ifelse(flip, pi*3/2, pi/2),
                       "average"=ifelse(flip, mnanchor[2], mnanchor[1])
    )
    res <- wrap(baseline + interval * relpos)
  }
  res
}