xdependence <- function(array3d,ts,u,fun) {

# Compute extreme dependence measures between a given three-dimensional array
# with first two space dimensions (e.g. longitude and latitude) and third
# time dimension, and a given time series with the same length as the time
# dimension of the three dimensional array
#
# Description:
#
#      Returns a map (matrix) of the extreme dependence measures described in
#      Coles et al. (1999) and in section 8.4 of Coles (2001)
#
# Usage:
#
#      xdependence(array3d,ts,u,fun)
#
# Input:
#
#      array3d: three-dimensional array with p longitude points and q latitude
#               points as the first two dimensions and n as the third time 
#               dimension 
#
#          ts: time series of data with the same length as the time dimension
#              of array3d for which extreme dependence with each grid point of
#              array3d will be computed
#
#           u: high threshold between 0 and 1 that will define the extreme 
#              level to compute dependence. Must be a high quantile (e.g. 0.95)
#
#         fun: String of characters that defines which extreme dependence 
#              measure is computed (e.g. "chi", "chibar" or "chicount").
#              If fun is "chi", computes chi(u) as in Eq.(3.2) of 
#              Coles et al. 1999 and in section 8.4 of Coles (2001). 
#              If fun is "chibar", computes chibar(u) as in section 3.3.2 of 
#              Coles et al. 1999 and in section 8.4 of Coles (2001).
#              If fun is "chicount", computes the conditional probability of
#              each grid point variable of array3d to be large (above a large 
#              threshold u) conditional on (given that) the variable ts is 
#              also large (with value above the same large threshold u)
#              by counting values above threshold and computing relative 
#              frequencies    
#          
# Output:
#
#      $out: an map (matrix) of the extreme dependence measure for each grid
#            point of array3d. 
#
# Authors:
#
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 21 Dec 2005
#      Christopher Ferro <c.a.t.ferro@reading.ac.uk>
#
# References: Coles, S., Heffernan, J., and Tawn, J., 1999: Dependence measures
#             of extreme value analyses. Extremes 2:4, 339-365.
#
#             Coles, S., 2001: An introduction to statistical modeling of 
#             extreme values. Spring series in statistics. 208pp.
#
# Examples:
#
#      x <- seq(-20, 20, 5)
#      y <- seq(30, 60, 5)
#      dim <- c(length(x), length(y), 1000)
#      data <- array(rnorm(prod(dim)), dim)
#      xdependence(data,data[4,5,],u=0.95,fun="chi")$out
#      xdependence(data,data[4,5,],u=0.95,fun="chibar")$out
#      xdependence(data,data[4,5,],u=0.95,fun="chicount")$out


out1<-aperm(apply(array3d,c(1,2),rank),c(2,3,1))/dim(array3d)[3]

if (fun == "chi"){
ts<-rank(ts)/length(ts)
chi <- function(y){
2-(log((sum(ts<=u & y<=u)+0.5)/(length(ts)+1))/(log((sum(ts<=u)+0.5)/(length(ts)+1))))
}
out<-apply(out1,c(1,2),chi)
out[out<0]<-0
}

if (fun == "chibar"){
ts<-rank(ts)/length(ts)
chibar <- function(y){
((2*log((sum(ts>=u)+0.5)/(length(ts)+1)))/log((sum(ts>=u & y>=u)+0.5)/(length(ts)+1)))-1
}
out<-apply(out1,c(1,2),chibar)
}

if (fun == "chicount"){
chicount <- function(y){
length(y[ts>=quantile(ts,u)][y[ts>=quantile(ts,u)]>=quantile(y,u)])/length(y[ts>=quantile(ts,u)])
}
out<-apply(array3d,c(1,2),chicount)
}

invisible(list(out=out))

}