# R/twoHmax.R In adiv: Analysis of Diversity

```twoHmax <-
function (C, epsilon = 1e-8, smooth = TRUE, comment = FALSE)
{

nonnegativedefinite <- function(x){
eigenval <- eigen(x, only.values = TRUE)\$values
nonzeroeigenval <- eigenval[abs(eigenval)>epsilon]
return(!any(nonzeroeigenval<0))
}
if (!inherits(C, "matrix"))
stop("Matrix expected")
if (epsilon <= 0)
stop("epsilon must be positive")
if (!nonnegativedefinite(C))
stop("C must be nonnegative definite")
C2 <- C^2
n <- dim(C2)[1]
relax <- 0
Z <- diag(1/diag(C))%*%C2%*%diag(1/diag(C))
D2 <- max(Z) - Z
x0 <- apply(D2, 1, sum)/sum(D2)
objective0 <- t(x0) %*% D2 %*% x0
if (comment == TRUE)
print("evolution of the objective function:")
xk <- x0
repeat {
repeat {
maxi.temp <- t(xk) %*% D2 %*% xk
if (comment == TRUE)
print(as.character(maxi.temp))
deltaf <- (-2 * D2 %*% xk)
sature <- (abs(xk) < epsilon)
if (relax != 0) {
sature[relax] <- FALSE
relax <- 0
}
yk <- (-deltaf)
yk[sature] <- 0
yk[!(sature)] <- yk[!(sature)] - mean(yk[!(sature)])
if (max(abs(yk)) < epsilon) {
break
}
alpha.max <- as.vector(min(-xk[yk < 0]/yk[yk < 0]))
alpha.opt <- as.vector(-(t(xk) %*% D2 %*% yk)/(t(yk) %*%
D2 %*% yk))
if ((alpha.opt > alpha.max) | (alpha.opt < 0)) {
alpha <- alpha.max
}
else {
alpha <- alpha.opt
}
if (abs(maxi.temp - t(xk + alpha * yk) %*% D2 %*%
(xk + alpha * yk)) < epsilon) {
break
}
xk <- xk + alpha * yk
}
if (prod(!sature) == 1) {
if (comment == TRUE)
print("KT")
break
}
vectD2 <- D2 %*% xk
u <- 2 * (mean(vectD2[!sature]) - vectD2[sature])
if (min(u) >= 0) {
if (comment == TRUE)
print("KT")
break
}
else {
if (comment == TRUE)
print("relaxation")
satu <- (1:n)[sature]
relax <- satu[u == min(u)]
relax <- relax[1]
}
}
if (comment == TRUE)
print(list(objective.init = objective0, objective.final = maxi.temp))
result <- as.vector(xk, mode = "numeric")
result[result < epsilon] <- 0
if(smooth){
respos <- (1:n)[result>0]
ressmooth <- (solve(Z[respos, respos])%*%rep(1, length(respos))/sum(solve(Z[respos, respos])))[, 1]
if(!any(ressmooth < -epsilon)){
result <- rep(0, n)
result[respos] <- ressmooth
result[result < epsilon] <- 0
}
}
result <- result / diag(C)
result <- result / sum(result)
result <- as.data.frame(result)
names(result) <- "pmax"
result\$pmax[result\$pmax < epsilon] <- 0
restot <- list()
p <- result[, 1]
restot\$value <- as.vector((t(p)%*%diag(C))^2)/as.vector(t(p)%*%C2%*%p)
restot\$vector <- result
return(restot)

}
```