Corr2Omni accompanying code
Code for generating the Figures in “Optimizing the Induced Correlation in Omnibus Joint Graph Embeddings”
Background Code
bitflip <- function(G,P,r){
n<-dim(P)[1]
H<-matrix(0,n,n)
for(i in 1:(n-1)){
for(j in 2:n){
if(G[i,j]==1){H[i,j]<-rbinom(1,1, ( P[i,j]+r*(1-P[i,j]) )) }
if(G[i,j]==0){H[i,j]<-rbinom(1,1, r*(1-P[i,j]) ) }
}
}
H<-H+t(H)
diag(H)<-rep(0,n)
return(H)
}
fromCtoM <- function(A,C){
# A is list of matrices
# in WOMNI setting (more general is more
# complicated....)
m=length(A)
n=dim(A[[1]])[1]
M<-matrix(0,m*n,m*n)
for(i in 1:m){
for(j in 1:m){
ind1 <- c( ((i-1)*n+1):(i*n))
ind2 <- c( ((j-1)*n+1):(j*n))
M[ind1,ind2] <- ( C[i,j,i]*A[[i]] + C[i,j,j]*A[[j]] )/2
}
}
return(M)
}
ase <- function(A, dim) {
require(irlba)
#diag(A) <- rowSums(A) / (nrow(A)-1) # diagaug line
if(nrow(A) >= 1000){
A.svd <- irlba(A, nu = dim, nv = dim)
A.svd.values <- A.svd$d[1:dim]
A.svd.vectors <- A.svd$v[,1:dim,drop=F]
if(dim == 1)
A.coords <- sqrt(A.svd.values) * A.svd.vectors
else
A.coords <- A.svd.vectors %*% diag(sqrt(A.svd.values))
} else{
A.svd <- svd(A)
A.svd.values <- A.svd$d[1:dim]
if(dim == 1)
A.coords <- matrix(A.svd$v[,1,drop=F] * sqrt(A.svd$d[1]),ncol=dim)
else
A.coords <- A.svd$v[,1:dim] %*% diag(sqrt(A.svd$d[1:dim]))
}
return(list(X=A.coords,D=A.svd.values))
}
make_Omni<- function(A,w)
{
n<-dim(A[[1]])[1]
m<-length(A)
ind<-list()
for(i in 1:m){
ind[[i]]<-c( (n*(i-1)+1): (n*i))
}
O<-matrix(0,n*m,n*m)
for(i in 1:m){
for(j in 1:m){
O[ ind[[i]] , ind[[j]] ]<-as.matrix((w[i]*A[[i]]+w[j]*A[[j]])/(w[i]+w[j]))
}}
for(i in 1:m){O[ ind[[i]] , ind[[i]] ]<-as.matrix(A[[i]])}
return(O)
}
getElbows <- function(dat, n = 3, threshold = FALSE, plot = F) {
if (is.matrix(dat))
d <- sort(apply(dat,2,sd), decreasing=TRUE)
else
d <- sort(dat,decreasing=TRUE)
if (!is.logical(threshold))
d <- d[d > threshold]
p <- length(d)
if (p == 0)
stop(paste("d must have elements that are larger than the threshold ",
threshold), "!", sep="")
lq <- rep(0.0, p) # log likelihood, function of q
for (q in 1:p) {
mu1 <- mean(d[1:q])
mu2 <- mean(d[-(1:q)]) # = NaN when q = p
sigma2 <- (sum((d[1:q] - mu1)^2) + sum((d[-(1:q)] - mu2)^2)) / (p - 1 - (q < p))
lq[q] <- sum( dnorm( d[1:q ], mu1, sqrt(sigma2), log=TRUE) ) +
sum( dnorm(d[-(1:q)], mu2, sqrt(sigma2), log=TRUE) )
}
q <- which.max(lq)
if (n > 1 && q < p) {
q <- c(q, q + getElbows(d[(q+1):p], n=n-1, plot=FALSE))
}
if (plot==TRUE) {
if (is.matrix(dat)) {
sdv <- d # apply(dat,2,sd)
plot(sdv,type="b",xlab="dim",ylab="stdev")
points(q,sdv[q],col=2,pch=19)
} else {
plot(dat, type="b")
points(q,dat[q],col=2,pch=19)
}
}
return(q)
}Code to implement the Corr2Omni Algorithm
BTrans <- function(W, Z, D){
DD <- as.matrix(dist(Z))
B <- - W*D/DD
B[is.na(B)] <- 0
diag(B) <- -rowSums(B)
return(B)
}
stress <- function(D,W,X){
DD <- as.matrix(dist(X))
WD <- W*(D-DD)^2
sigma <- sum(WD)/2
return(sigma)
}
corr2omni <- function (W, R1, R2, m, maxiter){
library(Matrix)
# R1 is m x m LATENT correlation matrix with 1's on the diagonal (if unknown,
# default should be Identity)
# R2 is TARGET correlation (R2=R1 is often typical)
L <- t(chol(R1))
Linv <- solve(L)
# equality constrint is A1%*%vec(AL)=matrix(m,m,1)
A1 <- kronecker( diag(m) , t(Linv %*% matrix(1,m,1)) )
b1<-m*rep(1,m)
# for WOMNI, we have the Aij+Aji=1 for i neq j
A22 <- t(combn(m, 2))
A2 <- matrix(0,choose(m,2), m^2)
for(i in 1:choose(m,2)){
a2 <- A22[i,]
ai <- a2[1]
aj <- a2[2]
A2[i, (ai-1)*m++aj]<-1
A2[i, (aj-1)*m++ai]<-1
}
A2 <- A2 %*% kronecker( diag(m), t(Linv) )
b2 <- rep(1, choose(m,2))
A <- rbind(A1,A2)
b <- c(b1,b2)
# inequality constraints
G <- list()
E <- cbind(matrix(1,m-1,1),-diag(m-1))
G[[1]]<- E
for(i in 2:m){
EE <- E
EE[,i-1]<-E[,i]
EE[,i]<-E[,i-1]
E<-EE
G[[i]]<-E
}
B1 <- bdiag(G)
B2 <- kronecker( diag(m), t(Linv) )
BB <- as.matrix(rbind(B1%*%B2,B2))
h <- rep(0,dim(BB)[1])
# initialize with classic omni
X0 <- matrix(1/2,m,m)
diag(X0)<-diag(X0)+m/2
X<-X0%*%L
D <- sqrt(2*m^2*(1-R2))
s0 <- 0
iter <- 0
s1 <- stress(D,W,X)
V <- diag(rowSums(W))-W
VV <- kronecker(V,diag(m))
library(piqp)
Z <- X
# iterate
while(iter<maxiter){
B <- BTrans(W,Z,D)
bb <- c( t(B%*%Z ) )
model <- piqp(2*VV,-2*bb, A, b, -BB, h)
res <- model$solve()
Z <- res$x
Z <- matrix(Z,m,m,byrow=T)
s0 <- s1
# print(c("s0 =",s0)) this will print the stress values
s1 <- stress(D,W,Z)
# print(c("s1 =",s1)) this will print the stress values
iter <-iter+1
}
return(list(config=Z%*%Linv,stress=s1))
}
fromAtoC <- function(A){
# in WOMNI setting (more general is more
# complicated....)
m <- dim(A)[1]
C <- array(0,c(m,m,m))
for(i in 1:m){
for(j in 1:m){
if(j==i){C[i,i,i]<-1}else{
C[i,j,j]<-A[i,j]
C[i,j,i]<-1-A[i,j]}
}
}
return(C)
}Corr2Omni simulations
Simulations in Section 4.1
set.seed(7654)
R <- matrix(2/3,3,3)+diag(3)/3
W <- matrix(1,3,3)-diag(3)
CE <- corr2omni(W, diag(3), R, 3, 400)
round(CE$config,4)## [,1] [,2] [,3]
## [1,] 2 1 0
## [2,] 0 2 1
## [3,] 1 0 2m=4
set.seed(6543)
R <- matrix(2/3,m,m)+diag(m)/3
W <- matrix(1,m,m)-diag(m)
CE <- corr2omni(W, diag(m), R, m, 1000)
round(CE$config,4)## [,1] [,2] [,3] [,4]
## [1,] 2.4259 0.0000 1.0000 0.5741
## [2,] 1.0000 2.4259 0.0000 0.5741
## [3,] 0.0000 1.0000 2.4259 0.5741
## [4,] 0.4259 0.4259 0.4259 2.7222RR <- 1-as.matrix(dist(CE$config))^2/(2*m^2)
RR## 1 2 3 4
## 1 1.0000000 0.7213010 0.7213010 0.7148272
## 2 0.7213010 1.0000000 0.7213010 0.7148272
## 3 0.7213010 0.7213010 1.0000000 0.7148272
## 4 0.7148272 0.7148272 0.7148272 1.0000000m=5
set.seed(5432)
R <- matrix(2/3,m,m)+diag(m)/3
W <- matrix(1,m,m)-diag(m)
CE <- corr2omni(W, diag(m), R, m, 5000)
round(CE$config,4)## [,1] [,2] [,3] [,4] [,5]
## [1,] 3 0 0 1 1
## [2,] 1 3 0 0 1
## [3,] 1 1 3 0 0
## [4,] 0 1 1 3 0
## [5,] 0 0 1 1 3RR <- 1-as.matrix(dist(CE$config))^2/(2*m^2)
RR## 1 2 3 4 5
## 1 1.00 0.72 0.68 0.68 0.72
## 2 0.72 1.00 0.72 0.68 0.68
## 3 0.68 0.72 1.00 0.72 0.68
## 4 0.68 0.68 0.72 1.00 0.72
## 5 0.72 0.68 0.68 0.72 1.00m=5
set.seed(54321)
R <- matrix(0.72,m,m)+0.28*diag(m)
W <- matrix(1,m,m)-diag(m)
CE <- corr2omni(W, diag(m), R, m, 5000)
round(CE$config,4)## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.000 0.344 0.000 1.000 0.656
## [2,] 0.656 3.000 0.344 0.000 1.000
## [3,] 1.000 0.656 3.000 0.344 0.000
## [4,] 0.000 1.000 0.656 3.000 0.344
## [5,] 0.344 0.000 1.000 0.656 3.000RR <- 1-as.matrix(dist(CE$config))^2/(2*m^2)
RR## 1 2 3 4 5
## 1 1.0000000 0.7242928 0.7208418 0.7208418 0.7242928
## 2 0.7242928 1.0000000 0.7242928 0.7208418 0.7208418
## 3 0.7208418 0.7242928 1.0000000 0.7242928 0.7208418
## 4 0.7208418 0.7208418 0.7242928 1.0000000 0.7242928
## 5 0.7242928 0.7208418 0.7208418 0.7242928 1.0000000m=6
set.seed(6543)
R <- matrix(2/3,m,m)+diag(m)/3
W <- matrix(1,m,m)-diag(m)
CE <- corr2omni(W, diag(m), R, m, 5000)
round(CE$config,4)## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 3.4711 0.0000 0.0000 0.5289 1.0000 1.0000
## [2,] 1.0000 3.4323 0.0000 0.5677 1.0000 0.0000
## [3,] 1.0000 1.0000 3.4722 0.0000 0.0000 0.5278
## [4,] 0.4711 0.4323 1.0000 3.6121 0.4845 0.0000
## [5,] 0.0000 0.0000 1.0000 0.5155 3.4845 1.0000
## [6,] 0.0000 1.0000 0.4722 1.0000 0.0000 3.5278RR <- 1-as.matrix(dist(CE$config))^2/(2*m^2)
RR## 1 2 3 4 5 6
## 1 1.0000000 0.7376601 0.7129782 0.7089057 0.7330337 0.7099561
## 2 0.7376601 1.0000000 0.7281490 0.7248096 0.7089436 0.7115142
## 3 0.7129782 0.7281490 1.0000000 0.7184095 0.7119099 0.7222222
## 4 0.7089057 0.7248096 0.7184095 1.0000000 0.7222531 0.7176988
## 5 0.7330337 0.7089436 0.7119099 0.7222531 1.0000000 0.7216022
## 6 0.7099561 0.7115142 0.7222222 0.7176988 0.7216022 1.0000000Code for Figure 1
#########################################
###### single generator experiment ######
#########################################
library(iGraphMatch)
set.seed(123456)
X <-sample_dirichlet(500,c(1,1,1))
X<-X[c(1:2),]
P <- t(X)%*% X
D1<-matrix(0,200,2)
D1<-list(D1,D1,D1)
D2<-matrix(0,200,2)
D2<-list(D2,D2,D2)
r<-c(0,0.25,0.5)
B0<-sample_dot_product(X, directed = FALSE)
for(j in 1:3){
for(i in 1:200){
B0<-as.matrix(B0[])
B1<-bitflip(B0,P,r[j])
B2<-bitflip(B0,P,r[j])
B3<-bitflip(B0,P,r[j])
OO <-make_Omni(list(B1,B2,B3),c(1,1,1))
OC <-rbind(
cbind(B1,B1,B3),
cbind(B1,B2,B2),
cbind(B3,B2,B3)
)
X1<-ase(OO,2)$X
Y1<-ase(OC,2)$X
D1[[j]][i,]<- X1[1,]-X1[501,]
D2[[j]][i,]<- Y1[1,]-Y1[501,]
}
}
E1<-lapply(D1,"*",sqrt(500))
E2<-lapply(D2,"*",sqrt(500))
par(mfrow=c(1,3))
for(i in 1:3){
plot(E1[[i]][,1],E1[[i]][,2],xlim=c(-3,3),ylim=c(-3,3),xlab="X1",ylab="X2")
points(E2[[i]][,1],E2[[i]][,2],col='red')
}
set.seed(123)
m1<-Mclust(E1[[1]],G=1,modelNames="EEI")
m1$parameters$variance$Sigma## [,1] [,2]
## [1,] 0.2558875 0.0000000
## [2,] 0.0000000 0.3062538m2<-Mclust(E2[[1]],G=1,modelNames="EEI")
m2$parameters$variance$Sigma## [,1] [,2]
## [1,] 0.3486078 0.0000000
## [2,] 0.0000000 0.3996435m1<-Mclust(E1[[2]],G=1,modelNames="EEI")
m1$parameters$variance$Sigma## [,1] [,2]
## [1,] 0.1588564 0.0000000
## [2,] 0.0000000 0.9792145m2<-Mclust(E2[[2]],G=1,modelNames="EEI")
m2$parameters$variance$Sigma## [,1] [,2]
## [1,] 0.2197546 0.000000
## [2,] 0.0000000 1.287586m1<-Mclust(E1[[3]],G=1,modelNames="EEI")
m1$parameters$variance$Sigma## [,1] [,2]
## [1,] 0.1058826 0.000000
## [2,] 0.0000000 2.024476m2<-Mclust(E2[[3]],G=1,modelNames="EEI")
m2$parameters$variance$Sigma## [,1] [,2]
## [1,] 0.1552535 0.000000
## [2,] 0.0000000 3.056915Code to make the social media example figures in Section 4.3
########################################
############### social media example ###
########################################
# The social media dataset is available for download at
# http://multilayer.it.uu.se/datasets.html
# it is the ff-tw-yt dataset
# The data is from
# Celli, Fabio, et al. "Social network data and practices: the case of friendfeed." International Conference on Social Computing, Behavioral Modeling, and Prediction. Berlin, Heidelberg: Springer Berlin Heidelberg, 2010.
# Magnani, Matteo, and Luca Rossi. "The ml-model for multi-layer social networks." 2011 International conference on advances in social networks analysis and mining. IEEE, 2011.
social <- read.csv("~/Documents/social.csv", header=FALSE)
ind_tw <- which(social[,3]=="tw")
ind_ff <- which(social[,3]=="ff")
ind_yt <- which(social[,3]=="yt")
social_tw <- social[ind_tw, c(1:2) ]
g_tw <- graph_from_edgelist(cbind(social_tw$V1,social_tw$V2) , directed=FALSE)
v_tw <- names(V(g_tw))
social_ff <- social[ind_ff, c(1:2) ]
g_ff <- graph_from_edgelist(cbind(social_ff$V1,social_ff$V2), directed=FALSE)
v_ff <- names(V(g_ff))
social_yt <- social[ind_yt, c(1:2) ]
g_yt <- graph_from_edgelist(cbind(social_yt$V1,social_yt$V2), directed=FALSE)
v_yt <- names(V(g_yt))
cv <- intersect(intersect(v_tw, v_ff), v_yt)
n <- length(cv)
ind1<-rep(0,n)
ind2<-rep(0,n)
ind3<-rep(0,n)
for(i in 1:n){
ind1[i] <- which(v_yt == cv[i])
ind2[i] <- which(v_ff == cv[i])
ind3[i] <- which(v_tw == cv[i])
}
A <- g_yt[]
A <- A[ind1,ind1]
B <- g_ff[]
B <- B[ind2,ind2]
C <- g_tw[]
C <- C[ind3,ind3]
diag(A) <- 0
diag(B) <- 0
diag(C) <- 0
a<-1
b<-1
c<-1
while((a+b+c) > 0){
zA <- which( as.numeric(rowSums(A)) == 0)
zB <- which( as.numeric(rowSums(B)) == 0)
zC <- which( as.numeric(rowSums(C)) == 0)
zt <- union(union(zA, zB), zC)
a<- sum(rowSums(A[-zt,-zt])==0)
b<- sum(rowSums(B[-zt,-zt])==0)
c<- sum(rowSums(C[-zt,-zt])==0)
A <- A[-zt,-zt]
B <- B[-zt,-zt]
C <- C[-zt,-zt]}
diag(A)<-rowSums(A)/421
diag(B)<-rowSums(B)/421
diag(C)<-rowSums(C)/421
G<-list(A,B,C)
RO<-matrix(0,3,3)
for(i in 1:3){
for(j in 1:3){
E <- G[[i]]
GG <- G[[j]]
EH <- E-diag(diag(E))
GH <- GG-diag(diag(GG))
tc1 <- norm(E-GG,"F")^2
RO[i,j]<- 1-tc1/( norm(E,"F")^2+norm(GG,"F")^2-
-2*sum(diag(E))*sum(diag(GG))/422-
sum(EH)*sum(GH) /choose(422,2) )
}
}
WW <-matrix(1,3,3)-diag(3)
CB <- corr2omni(WW, RO, RO, 3, 1000)
CC <- fromAtoC(CB$config)
for(i in 1:3){G[[i]]<-as.matrix(G[[i]])}
OO<-make_Omni(G,rep(1,3))
M <- fromCtoM(G,CC)
RO # called R## [,1] [,2] [,3]
## [1,] 1.0000000 0.2772830 0.2808594
## [2,] 0.2772830 1.0000000 0.5516354
## [3,] 0.2808594 0.5516354 1.0000000L <-t(chol(RO))
RR <- 1-as.matrix(dist(CB$config%*%L))^2/(2*3^2)
RR # called \hat R_{c2O}## 1 2 3
## 1 1.0000000 0.6787925 0.7294055
## 2 0.6787925 1.0000000 0.9501817
## 3 0.7294055 0.9501817 1.0000000RC <- 1-as.matrix(dist(( diag(3)+matrix(0.5,3,3) )%*%L))^2/(2*3^2)
RC # called \hat R_O## 1 2 3
## 1 1.0000000 0.9196981 0.9200955
## 2 0.9196981 1.0000000 0.9501817
## 3 0.9200955 0.9501817 1.0000000norm(RR-RO,"F")## [1] 1.021022norm(RO-RC,"F")## [1] 1.400114set.seed(12345678)
W<-irlba(OO,30)
d<-getElbows(W$d)
Y<-ase(OO,d[2])$X
W<-irlba(M,30)
d<-getElbows(W$d)
YY<-ase(M,d[2])$X
Z<-list()
for(i in 1:3){
Z[[i]]<-Y[c( ((i-1)*422+1):(i*422)),]
}
DZ<-matrix(0,3,3)
for(i in 1:3){
for(j in 1:3){
DZ[i,j]<-norm(Z[[i]]-Z[[j]],"F")^2
}
}
ZZ<-list()
for(i in 1:3){
ZZ[[i]]<-YY[c( ((i-1)*422+1):(i*422)),]
}
DZZ<-matrix(0,3,3)
for(i in 1:3){
for(j in 1:3){
DZZ[i,j]<-norm(ZZ[[i]]-ZZ[[j]],"F")^2
}
}
par(mfrow=c(1,3))
par(mar = c(0.5,0.5,0.5,0.5))
# Figure 3
AA<-A
diag(AA)<-rep(0,422)
aa<-graph_from_adjacency_matrix(AA,mode="undirected")
DA<-distances(aa)
DA[DA==Inf]=0
image(DA[,422:1], col=gray.colors(200,start=0, end=1), axes=FALSE)
AA<-B
diag(AA)<-rep(0,422)
aa<-graph_from_adjacency_matrix(AA,mode="undirected")
DA<-distances(aa)
DA[DA==Inf]=0
image(DA[,422:1], col=gray.colors(200,start=0, end=1), axes=FALSE)
AA<-C
diag(AA)<-rep(0,422)
aa<-graph_from_adjacency_matrix(AA,mode="undirected")
DA<-distances(aa)
DA[DA==Inf]=0
image(DA[,422:1], col=gray.colors(200,start=0, end=1), axes=FALSE)
# Figure 2
image(as.matrix(dist(Y))[,1266:1], col=gray.colors(200,start=0, end=1), axes=FALSE)
image(as.matrix(dist(YY))[,1266:1], col=gray.colors(200,start=0, end=1), axes=FALSE)
Figures from Section 4.4
############################################
#### Brain graph example ###################
############################################
patnum <- 25426
BB <-list()
for( i in 1:10){
pat<-patnum+i
BB[[i]]<-list()
for(j in 1:3){
loc <- paste0("~/Downloads/HNU1_Connectomes/Desikan_space-MNI152NLin6_res-2x2x2/sub-00",
pat,"_ses-",j,"_connectome.csv")
X <- read_delim(loc,delim = " ", col_names = FALSE)
g <- graph_from_data_frame(X, directed=F, vertices=c(1:70))
g <- as_adjacency_matrix(g, attr="X3")
diag(g)<-rowSums(g)/69
BB[[i]][[j]]<-g[]
}
}
G<- unlist(BB, recursive=F)
R <- matrix(0.54,30,30)+0.46*diag(30)
set.seed(1234567)
WW <- matrix(1,30,30)-diag(30)
CB <- corr2omni(WW, diag(30), R, 30, 500)
CCCB <- fromAtoC(CB$config)
for(i in 1:30){G[[i]]<-as.matrix(G[[i]])}
OO<-make_Omni(G,rep(1,30))
M <- fromCtoM(G,CCCB)
CBD<-CB$config-diag(diag(CB$config))
image(CBD[,30:1], col=gray.colors(200,start=0, end=1),axes=FALSE)
CBCO <-matrix(0.5,30,30)-0.5*diag(30)
image(CBCO[,30:1], col=gray.colors(200,start=0, end=0.5),axes=FALSE)
set.seed(12345678)
W<-irlba(OO,30)
d<-getElbows(W$d)
Y<-ase(OO,d[2])$X
W<-irlba(M,30)
d<-getElbows(W$d)
YY<-ase(M,d[2])$X
Z<-list()
for(i in 1:30){
Z[[i]]<-Y[c( ((i-1)*70+1):(i*70)),]
}
DZ<-matrix(0,30,30)
for(i in 1:30){
for(j in 1:30){
DZ[i,j]<-norm(Z[[i]]-Z[[j]],"F")
}
}
set.seed(12345678)
W<-irlba(DZ,10)
d<-getElbows(W$d)
cz <- cmdscale(DZ,d[2])
ZZ<-list()
for(i in 1:30){
ZZ[[i]]<-YY[c( ((i-1)*70+1):(i*70)),]
}
DZZ<-matrix(0,30,30)
for(i in 1:30){
for(j in 1:30){
DZZ[i,j]<-norm(ZZ[[i]]-ZZ[[j]],"F")
}
}
require(dendextend)
par(mfrow=c(1,2))
par(mar = c(1,0,1,0))
# Figure 5
g10 <- DZ %>%
as.dist() %>%
hclust(method = "ward.D") %>%
as.dendrogram() %>% set("branches_lwd", 4)
g10 <- color_branches(g10, 10)
plot(g10,yaxt='n')
gg10 <- DZZ %>%
as.dist() %>%
hclust(method = "ward.D") %>%
as.dendrogram() %>% set("branches_lwd", 4)
gg10 <- color_branches(gg10, 10)
plot(gg10,yaxt='n')
ll<-sort(rep(c(1:10),3))
d10 <- DZ %>%
as.dist() %>%
hclust(method = "ward.D") %>%
cutree(10)
adjustedRandIndex(ll,d10)## [1] 0.8117994dd10 <- DZZ %>%
as.dist() %>%
hclust(method = "ward.D") %>%
cutree(10)
adjustedRandIndex(ll,dd10)## [1] 0.9118541y<-rbind(ll,d10,dd10)
y<-data.frame(y)
rownames(y)<-c("Truth",
"Classical OMNI",
"corr2Omni")
colnames(y)<-c("1a","1b","1c","2a","2b","2c","3a","3b","3c","4a","4b","4c","5a",
"5b","5c","6a","6b","6c","7a","7b","7c","8a","8b","8c","9a","9b",
"9c","10a","10b","10c")
y$ARI<-round(c(1,adjustedRandIndex(ll,d10),adjustedRandIndex(ll,dd10)),3)Aplysia example
#############################################################
######### Aplysia example ###################################
#############################################################
# load O1 from Importance of being Correlated paper
O1 <- readRDS("/Users/vincelyzinski/KPbrains.RDS")
RO<-matrix(0,24,24)
for(i in 1:24){
for(j in 1:24){
E <- O1[[i]]
G <- O1[[j]]
EH <- E-diag(diag(E))
GH <- G-diag(diag(G))
tc1 <- norm(E-G,"F")^2
RO[i,j]<- 1-tc1/( norm(E,"F")^2+norm(G,"F")^2-
-2*sum(diag(E))*sum(diag(G))/82-
sum(EH)*sum(GH) /choose(82,2) )
}
}
set.seed(12345)
W <- matrix(1,24,24)-diag(24)
C <- corr2omni(W, RO, RO, 24, 500)
CCC <- fromAtoC(C$config)
par(mfrow=c(1,4))
par(mar = c(1,0.2,1,0.2))
# Figure 6
L <-t(chol(RO))
RR <- 1-as.matrix(dist(C$config%*%L))^2/(2*24^2)
image(-RO[,24:1], col=gray.colors(200,start=0, end=1), axes=FALSE)
image(-RR[,24:1], col=gray.colors(200,start=0, end=1), axes=FALSE)
CmD <- C$config-diag(diag(C$config))
image(-CmD[,24:1], col=gray.colors(200,start=0, end=1),axes=FALSE)
COmD <- matrix(0.5,24,24)-0.5*diag(24)
image(-COmD[,24:1], col=gray.colors(200,start=0.5, end=1),axes=FALSE)
library(irlba)
set.seed(12345)
M <- fromCtoM(O1,CCC)
OO<-make_Omni(O1,rep(1,24))
W<-irlba(OO,30)
d<-getElbows(W$d)
Y<-ase(OO,d[2])$X
W<-irlba(M,30)
d<-getElbows(W$d)
YY<-ase(M,d[2])$X
Z<-list()
for(i in 1:24){
Z[[i]]<-Y[c( ((i-1)*82+1):(i*82)),]
}
DZ<-matrix(0,24,24)
for(i in 1:24){
for(j in 1:24){
DZ[i,j]<-norm(Z[[i]]-Z[[j]],"F")
}
}
ZZ<-list()
for(i in 1:24){
ZZ[[i]]<-YY[c( ((i-1)*82+1):(i*82)),]
}
DZZ<-matrix(0,24,24)
for(i in 1:24){
for(j in 1:24){
DZZ[i,j]<-norm(ZZ[[i]]-ZZ[[j]],"F")
}
}
require(dendextend)
par(mfrow=c(2,4))
par(mar = c(1,0,1,0))
# Figure 8 Dendrograms
g2 <- DZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
as.dendrogram() %>% set("branches_lwd", 4)
g2 <- color_branches(g2, 2)
attr(g2, "height")<-attr(g2, "height")/2
plot(g2,yaxt='n')
g3 <- DZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
as.dendrogram() %>% set("branches_lwd", 4)
g3 <- color_branches(g3, 3)
attr(g3, "height")<-attr(g3, "height")/2
plot(g3,yaxt='n')
g4 <- DZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
as.dendrogram() %>% set("branches_lwd", 4)
g4 <- color_branches(g4, 4)
attr(g4, "height")<-attr(g4, "height")/2
plot(g4,yaxt='n')
g5 <- DZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
as.dendrogram() %>% set("branches_lwd", 4)
g5 <- color_branches(g5, 5)
attr(g5, "height")<-attr(g5, "height")/2
plot(g5,yaxt='n')
f2 <- DZZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
as.dendrogram() %>% set("branches_lwd", 4)
f2 <- color_branches(f2, 2)
attr(f2, "height")<-attr(f2, "height")/10
plot(f2,yaxt='n')
f3 <- DZZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
as.dendrogram() %>% set("branches_lwd", 4)
f3 <- color_branches(f3, 3)
attr(f3, "height")<-attr(f3, "height")/10
plot(f3,yaxt='n')
f4 <- DZZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
as.dendrogram() %>% set("branches_lwd", 4)
f4 <- color_branches(f4, 4)
attr(f4, "height")<-attr(f4, "height")/10
plot(f4,yaxt='n')
f5 <- DZZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
as.dendrogram() %>% set("branches_lwd", 4)
f5 <- color_branches(f5, 5)
attr(f5, "height")<-attr(f5, "height")/10
plot(f5,yaxt='n')
d5<-DZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
cutree(4)
dd5<-DZZ %>%
as.dist() %>%
hclust(method = "ward.D2") %>%
cutree(4)
ll <-c(1,2,3,3,rep(4,20))
adjustedRandIndex(d5,ll)## [1] 0.03545317adjustedRandIndex(dd5,ll)## [1] 0.1951258set.seed(0987)
W1<-irlba(DZ,10)
d<-getElbows(W1$d)
dz<-data.frame(cmdscale(DZ,d[2]))
W2<-irlba(DZZ,10)
d<-getElbows(W2$d)
dzz<-data.frame(cmdscale(DZZ,d[2]))
Phase <- c("Resting","Stimulus",rep("Gallop",2),rep("Crawl",20))
dz$Phase<-Phase
dzz$Phase<-Phase
f1<-ggplot(data=dz,aes(x=X1,y=X2))
f1<-f1+
geom_path()+
geom_point(aes(col=Phase, shape=Phase),size=3)+
scale_shape_manual(values=c(15,16,17,7))
f2<-ggplot(data=dz,aes(x=X1,y=X3))
f2<-f2+
geom_path()+
geom_point(aes(col=Phase, shape=Phase),size=3)+
scale_shape_manual(values=c(15,16,17,7))
f3<-ggplot(data=dz,aes(x=X2,y=X3))
f3<-f3+
geom_path()+
geom_point(aes(col=Phase, shape=Phase),size=3)+
scale_shape_manual(values=c(15,16,17,7))
g1<-ggplot(data=dzz,aes(x=X1,y=X2))
g1<-g1+
geom_path()+
geom_point(aes(col=Phase, shape=Phase),size=3)+
scale_shape_manual(values=c(15,16,17,7))
g2<-ggplot(data=dzz,aes(x=X1,y=X3))
g2<-g2+
geom_path()+
geom_point(aes(col=Phase, shape=Phase),size=3)+
scale_shape_manual(values=c(15,16,17,7))
g3<-ggplot(data=dzz,aes(x=X2,y=X3))
g3<-g3+
geom_path()+
geom_point(aes(col=Phase, shape=Phase),size=3)+
scale_shape_manual(values=c(15,16,17,7))
e1 <- ggplot(data=data.frame(cbind(1:10,W1$d)),aes(x=X1,y=X2))+
geom_point()+labs(x="Index",y="Singular Value")
e2 <- ggplot(data=data.frame(cbind(1:10,W2$d)),aes(x=X1,y=X2))+
geom_point()+labs(x="Index",y="Singular Value")
# Figure 7
library(patchwork)## Warning: package 'patchwork' was built under R version 4.2.3(e1|f1|f2|f3)/(e2|g1|g2|g3)+plot_layout(guides = "collect")