###############################################################################
# UAGRA R Scripts                                 Incremental_Kernel_Analysis.R
###############################################################################
# Incremental Kernel Analysis procedure
###############################################################################
#
# Version 1.1
#
# Author: Damiano G. Preatoni (prea@uinsubria.it)
#
# Description: estimates core areas accodring to Hanski et al. 2000.
#
# Usage: res <- corearea(loc, id, method = "kernel")
#
# Requires: adehabitat
#
# References:
# Hanski, I. K.; Stevens, P. C.; Ihalempia, P. & Selonen, V. (2000)
# Home-range size, movements, and nest-site use in the Siberian flyng squirrel,
# Pteromys volans Journal of Mammalogy, 2000, 81, 798-809 
#
###############################################################################
# created prea mer ago  1 11:26:45 CEST 2007
# updated prea mar ago  7 11:23:50 CEST 2007
# revision history:
#  1.1 - added clustering
#
#############################################################################
#    Copyright (C) 2007  Damiano G. Preatoni
#
#    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 3 of the License, or
#    (at your option) any later version.
#
#    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.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.
#############################################################################


if (!require(adehabitat)) {
  stop("This procedure requires adehabitat to be installed")
}

corearea <- function(xy, id = NULL, levels = seq(20,95,by = 5), method = "mcp", unin = "m", unout = "m", plot = TRUE) {
  alpha <- 0.05
  ret <- list()
  if(!method %in% c("clusthr","mcp","kernel")) {
    stop("method should be either \"clusthr\" or \"mcp\" or \"kernel\"") }
  if(method == "clusthr") {
    clust <- clusthr(xy, id)
    ret[['area']] <- clusthr.area(clust, percent = levels, unin = unin, unout = unout, plotit = FALSE)
    }
  if(method == "mcp") {
    ret[['area']] <- mcp.area(xy, id, percent = levels, unin = unin, unout = unout, plotit = FALSE)
    }
  if(method == "kernel") {
    ret[['area']] <- kernel.area(xy, id, h = "href", levels = levels, unin = unin, unout = unout)
    }
  #print(is(ret[['area']]))
  level.areas <- stack(as.data.frame(t(ret[['area']])))
  names(level.areas) <- c('area','coropleth')
  level.areas$coropleth <- ordered(level.areas$coropleth,levels=levels) # as factor!
  ret[['aov']] <- aov(area ~ coropleth, data = level.areas)
  # if coropleth factor is significant go on, else stop here
  if (!summary(ret[['aov']])[[1]][,5][1] <= alpha) {
    warning("No significant differences among coropleth areas.", call. = FALSE)
  }
  else {
    # do TukeyHSD
    ret[['tukey']] <- TukeyHSD(ret[['aov']],"coropleth",ordered = TRUE)
    # pick only consecutive pairs of coropleths
    n <- length(levels)
    i <- 1
    t <- vector()
    for(j in 1:n) {
      t[j] <- i
      i <- (i+n-j)
    }
    p <- ret[['tukey']]$coropleth[t[1:length(t)-1],4]
    if (length(p[p <= alpha]) == 0) {
      warning("No significant differences between consecutive coropleths found.", call. = FALSE)
    }
    else {
      val <- as.numeric(strsplit(names(p[p <= alpha]),"-")[[1]])
      cat("Tukey multiple comparison of means between consecutive levels\n\n")
      print(ret[['tukey']]$coropleth[t[1:length(t)-1],])
      cat("\n\nSignificant differences are present between levels",val[2],"and",val[1],"%\n")
      if (plot) {
        #@TODO(prea) make the plot using average +- standard deviation
        ret[['plot']] <- boxplot(area ~ coropleth,data=level.areas,xlab='% UD Coropleth',ylab=paste('Area (',unout, ")"))
        abline(v = match(val[2],levels),col="red")
      }
    }
  }
  invisible(ret)

}