Rao-Stirling diversity index I need help to calculate the Rao-Stirling diversity index. I tried it several times but cannot achieve equal results to R packages yet (e.g. diverse). 
I used the Rao-Stirling diversity index as suggested by several references 
$$ \Delta=\sum_{ij} p_j p_i d_{ij} $$
including the reciprocal cosine similarity for computing the disparity $$ d_{ij}=(1-'cosine') $$
When I would have the following simple matrix, please could one show me how to do it?
\begin{bmatrix}
0 & 10 & 20 \\
10 & 0 & 30 \\
20 & 30 & 0
\end{bmatrix}

Edit 02/01/20
delta <- function(matrix){ 
proportion.vector <- prop.vector(matrix) 
distance.matrix <- dist.sim.matrix(matrix) 
result <- matrix(rep(NA,nrow(matrix)*ncol(matrix)),nrow=nrow(matrix)) 
for (i in 1:nrow(matrix)){ 
    for (j in 1:ncol(matrix)){ 

result[i,j] <- proportion.vector[i]*proportion.vector[j]*distance.matrix[i,j] } } 

return(colSums(result))
}


Oh sorry, I missed the point. Here is the code that I have used. 
## Computing Rao-Stirling Diversity
#install.packages('diverse')
library(diverse)

rows <- c(1,1,1,2,2,2,3,3,3)
cols <- c(1,2,3,1,2,3,1,2,3)
values <- c(0,10,20,10,0,30,20,30,0)
my.dataframe <- data.frame(rows,cols,values)

RAO = diversity(data=my.dataframe, type="rao-stirling", method="cosine")

 A: I just had a cursory look at the literature. Here comes some code I just wrote:
    # computes cosine similarity between vec1 and vec2
    cos.sim <- function(vec1,vec2){
    return(sum(vec1*vec2)/sqrt(sum(vec1*vec1)*sum(vec2*vec2)))
    }

    # returns D_ij matrix
    dist.sim.matrix<- function(M){
         result <- matrix(rep(NA,nrow(M)^2),nrow=nrow(M))
         for (i in 1:nrow(M)){
            for (j in 1:nrow(M)){
                result[i,j] <- cos.sim(M[,i],M[,j])
            }
         }
         return(1-result)
     }

     # returns proportion vector
     prop.vector <- function(matrix){
          result <- numeric(ncol(matrix))
          total.sum <- sum(matrix)
          for (i in 1:ncol(matrix)){
             result[i] <- sum(matrix[,i])/total.sum
          }
          return(result)
       }
      # returns delta
     delta <- function(matrix){
        proportion.vector <- prop.vector(matrix)
        distance.matrix <- dist.sim.matrix(matrix)
        result <- matrix(rep(NA,nrow(matrix)^2),nrow=nrow(matrix))
        for (i in 1:nrow(matrix)){
           for (j in 1:nrow(matrix)){
             result[i,j] <-            proportion.vector[i]*proportion.vector[j]*distance.matrix[i,j]
           }
       }
      return(sum(result))
     }

     your.matrix <- matrix(c(0,10,20,10,0,30,20,30,0),nrow=3,byrow=TRUE)
     # Result for your matrix
     delta(your.matrix)

I hope this is correct.
