ubercoldist

ubercoldist <- function(Qi, ni=c(1,2,2,4), w=0.2){

Qi <- Qi[, 1:length(ni)]
Qi <- log(Qi)


e <- setNames(w/sqrt(ni), names(Qi))

# prepare output
pairsid <- t(combn(nrow(Qi),2))

patch1 <- row.names(Qi)[pairsid[, 1]]
patch2 <- row.names(Qi)[pairsid[, 2]]

res <- data.frame(patch1, patch2, stringsAsFactors=F)
  


#############
# NUMERATOR #
#############

# all n-2 combinations (first part numerator)
n1combs <- combn(names(Qi),dim(Qi)[2]-2)
  # ADD IF STATEMENT HERE IN CASE dim <= 2

# get those combinations of ei and prod(ei)^2

num1 <- setNames(apply(n1combs, 2, function(x) prod(e[x])^2), apply(n1combs, 2, paste, collapse=""))
  # ADD IF STATEMENT HERE IN CASE dim <= 2

# remaining 2 combinations (second part numerator)
n2combs <- apply(n1combs, 2, function(x) names(Qi)[ !names(Qi) %in% x ] )

# f_d and f_e

deltaqiqj <- 
lapply(1:length(num1), function(y)
  fidiff <- do.call('rbind', apply(res, 1, function(x)
  (Qi[x[1], n2combs[,y]] - Qi[x[2], n2combs[,y]])
  ))
)

names(deltaqiqj) <- apply(n2combs, 2, paste, collapse='')

# (f_d-f_e)^2

num2 <- do.call('cbind',lapply(deltaqiqj, function(x) apply(x, 1, function(z) diff(z)^2)))

# (e_abc)^2*(f_d-f_e)^2

etimesq <- num2 %*% diag(num1)

# sum numerator

numerator <- rowSums(etimesq)

###############
# DENOMINATOR #
###############

# all n-1 combinations
dcombs <- combn(names(Qi),dim(Qi)[2]-1)

den <- setNames(apply(dcombs, 2, function(x) prod(e[x])^2), apply(dcombs, 2, paste, collapse=""))

denominator <- sum(den)

###########
# DELTA S #
###########

res$dS <- sqrt(numerator/denominator)

res
}



# RUN SOME TESTS
require(pavo)
data(sicalis)

test1 <- ubercoldist(vismodel(sicalis, relative=FALSE, achro='none'))
test2 <- coldist(vismodel(sicalis, relative=FALSE, achro='none'), achro=FALSE)

all.equal(test1, test2) # only attributes different, values the same

test3 <- ubercoldist(vismodel(sicalis, relative=FALSE, achro='none'), ni=c(1,2,2))
test4 <- coldist(vismodel(sicalis, relative=FALSE, achro='none'), achro=FALSE, vis='tri')

all.equal(test3, test4) # only attributes different, values the same


# HOW ABOUT SOME MANTIS SHRIMP ACTION BRO

mantis <- data.frame(matrix(rnorm(10*10, 10, 2), ncol=10))
names(mantis) <- LETTERS[1:10]

head(ubercoldist(mantis, ni=c(1,1,1,1,1,2,2,2,2,2)))