################################################################################
# File: SXvsKerSX.R
# Aim : Test if the two covariance matrices are equal
#-------------------------------------------------------------------------------
# Author: Yuan Huili (Peking University)
# Email : hlyuan@pku.edu.cn
# Date  : 05/29/2015
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
product<-function(X,Y)
{
m = ncol(X)
n = ncol(Y)
product=t(X)%*%Y
}
#-------------------------------------------------------------------------------
kproduct<-function(X,Y,coef)
{
m = ncol(X)
n = ncol(Y)
temp = matrix(0,m,n)
for (i in 1:m)
{    for (j in 1:n)
 {       temp[i,j] = sum((X[,i]-Y[,j])^2)
 }
}
temp = temp/coef^2
kproduct = exp(-temp)
}
#-------------------------------------------------------------------------------
#
A<-function(n,X)
{
answer=0
temp = product(X,X)
temp = temp - diag(diag(temp))
temp1 = temp^2
answer = sum(temp1)/(n*(n-1))
temp2 = temp%*%temp
temp2 = temp2 - diag(diag(temp2))
tempans = sum(temp2)
answer = answer - 2*tempans/(n*(n-1)*(n-2))
temp3 = kronecker(temp,temp)
tempans = sum(temp3)-4*tempans-2*sum(diag(temp3))
A = answer + tempans/(n*(n-1)*(n-2)*(n-3))
}
#-------------------------------------------------------------------------------
AA<-function(n,X,coef)
{
answer=0
temp = kproduct(X,X,coef)
temp = temp - diag(diag(temp))
temp1 = temp^2
answer = sum(temp1)/(n*(n-1))
temp2 = temp%*%temp
temp2 = temp2 - diag(diag(temp2))
tempans = sum(temp2)
answer = answer - 2*tempans/(n*(n-1)*(n-2))
temp3 = kronecker(temp,temp)
tempans = sum(temp3)-4*tempans-2*sum(diag(temp3))
AA = answer + tempans/(n*(n-1)*(n-2)*(n-3))
}
#-------------------------------------------------------------------------------
C<-function(n1,n2,X1,X2)
{
answer = 0
temp = product(X1,X2)
temp1 = temp^2
answer = sum(temp1)/(n1*n2)
temp2 = temp %*%t(temp)
temp2 = temp2 - diag(diag(temp2))
temp3 = t(temp) %*% temp
temp3 = temp3 - diag(diag(temp3))
answer = answer - sum(temp2)/(n1*n2*(n1-1))-sum(temp3)/(n1*n2*(n2-1))
temp4 = kronecker(temp,temp)
temp4 = sum(temp4)-sum(temp2)-sum(temp3)-sum(temp1)
C = answer + temp4/(n1*n2*(n1-1)*(n2-1))
}
#-------------------------------------------------------------------------------
CC<-function(n1,n2,X1,X2,coef)
{
answer = 0
temp = kproduct(X1,X2,coef)
temp1 = temp^2
answer = sum(temp1)/(n1*n2)
temp2 = temp %*%t(temp)
temp2 = temp2 - diag(diag(temp2))
temp3 = t(temp) %*% temp
temp3 = temp3 - diag(diag(temp3))
answer = answer - sum(temp2)/(n1*n2*(n1-1))-sum(temp3)/(n1*n2*(n2-1))
temp4 = kronecker(temp,temp)
temp4 = sum(temp4)-sum(temp2)-sum(temp3)-sum(temp1)
CC = answer + temp4/(n1*n2*(n1-1)*(n2-1))
}
#-------------------------------------------------------------------------------
#Songxi Chen's method
SXCHEN<-function(n1,n2,X1,X2)
{
a<-A(n1,X1)
b<-A(n2,X2)
SXCHEN<- a+b- 2*C(n1,n2,X1,X2)
}
#-------------------------------------------------------------------------------
#Songxi Chen's method under Gauss Kernel
GSXCHEN<-function(n1,n2,X1,X2,coef)
{
a<-AA(n1,X1,coef)
b<-AA(n2,X2,coef)
GSXCHEN<-a + b- 2*CC(n1,n2,X1,X2,coef)
}
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------