cluster.loadings <- 
function (keys, r.mat, correct = TRUE,SMC=TRUE) 
{  cl <- match.call()
    if (!is.matrix(keys)) {
        keys <- as.matrix(keys)}
    r.mat[is.na(r.mat)] <- -9999999 
    
 
   
    item.sd <- sqrt(diag(r.mat))
    item.covar <- r.mat %*% keys          #item by cluster covariances
    covar <- t(keys) %*% item.covar  #variance/covariance of clusters
   
    var <- diag(covar)
    sd.inv <- 1/sqrt(var)
    
    
  #items corrected for communality lead to the Guttman G6 estimate
    if(SMC)  {r.smc <- smc(r.mat)
              r.smc[r.smc < 0 ] <- 1  #for a very weird condition
            diag(r.mat) <- r.smc } else {
            diag(r.mat) <- 0
            item.max <- apply(r.mat,1,max)
            diag(r.mat) <- item.max}
    c.item.covar <- r.mat %*% keys #this uses the communality estimate and thus corrects for item overlap
    c.covar <- t(keys) %*% c.item.covar 
    c.var <- diag(c.covar)
    G6 <- c.var/var
    n.keys <- dim(keys)[2] 
    if(n.keys >1) {
    c.item.cor <- c.item.covar %*% sqrt(diag(1/c.var))/item.sd } else {c.item.cor <- c.item.covar/sqrt(c.var*item.sd) }
    
  
     
    key.count <- diag(t(keys) %*% keys)    #how many items in each cluster?
    if (correct) {
                  cluster.correct <- diag((key.count/(key.count - 1)))
                   for (i in 1:dim(keys)[2]) {
                         if (key.count[i]<2 ) {   #fix the case of 1 item keys
                   	     cluster.correct[i,i] <- 1
                   	} else { cluster.correct[i,i] <- key.count[i]/(key.count[i]-1)
                   	item.covar[,i] <- item.covar[,i] - keys[,i]} #subtract the variance of the item
                   }   #i loop
                   correction.factor <- keys %*% cluster.correct #put back average correlation for the item if it loads on the key
                   correction.factor[ correction.factor < 1] <- 1
                  item.covar <- item.covar * correction.factor
                  }
   
      
    ident.sd <- diag(sd.inv, ncol = length(sd.inv))
   
    c.correl <- ident.sd %*% covar %*% ident.sd
    p.loading <- try(c.item.cor %*% solve(c.correl))
    if(inherits(p.loading,"try-error")) {message('the correlation matrix was singular, pattern loadings not found, proceed with caution')
                         p.loading  <- c.item.cor}
    
     c.item.cor[abs(c.item.cor)  > 99999] <- NA
     c.correl[abs(c.correl) > 99999] <- NA
    
   
    
    key.alpha <- ((var - key.count)/var) * (key.count/(key.count - 1))
    key.alpha[is.nan(key.alpha)] <- 1
    key.alpha[!is.finite(key.alpha)] <- 1
    key.av.r <- key.alpha/(key.count - key.alpha*(key.count-1))  #alpha 1 = average r
    colnames(c.item.cor) <- colnames(keys)
    colnames(p.loading) <- colnames(keys)
    colnames(c.correl) <- colnames(keys)
    rownames(c.correl) <- colnames(keys)
    rownames(c.item.cor) <- rownames(r.mat)
    
    if( ncol(keys) >1)  {cluster.corrected <- correct.cor(c.correl, t(key.alpha))} else {cluster.corrected <- c.correl}
    
    results <- list(loadings=c.item.cor,pattern=p.loading, cor=c.correl,corrected=cluster.corrected, sd = sqrt(var), alpha = key.alpha,av.r = key.av.r,
            size = key.count,G6=G6,Call=cl)
    class(results) <- c("psych","cluster.loadings")
    return(results)
    }