####This R code is to define CMDCH with Gussian kernel
###main function 
library(kernlab) #calculate kernel matrxi



# inner product of matrix
u.inner.H <- function(X, Y){
  n <- dim(X)[1]
  ip <- sum(X * Y) / n / (n-3) 
  return(ip)
}


# inner product of matrix Y and conditional set X1
uYX1.center.H <- function(Y,X1) {
  n <- dim(Y)[1]
  A1 <- as.matrix((Y%*%t(Y)))
  rbf1 <- function(x,y){
    return(exp(-0.5*t(x-y)%*%(x-y))) #guassian kernel with bandwidth 2
  }
  class(rbf1)="kernel"
  A2 <- kernelMatrix(rbf1,X1)
  A <- A1*A2
  diag(A) <- 0
  R <- rowSums(A)
  C <- colSums(A)
  S <- sum(A)
  r <- matrix(rep(R, n), n, n) / (n-2)
  c <- t(matrix(rep(C, n), n, n)) / (n-2 )
  t <- matrix(S / (n-1 ) / (n-2 ), n, n)
  UA <- A - r - c + t
  diag(UA) <- 0
  return(UA)
}

#matrix of screening set X2 
uX2.center.H <- function(X2) {
  
  n <- dim(X2)[1]
  rbf1 <- function(x,y){
    return(exp(-0.5*t(x-y)%*%(x-y))) 
  }
  A <- kernelMatrix(rbf1,X2)
  diag(A) <- 0
  R <- rowSums(A)
  C <- colSums(A)
  S <- sum(A)
  r <- matrix(rep(R, n), n, n) / (n-2)
  c <- t(matrix(rep(C, n), n, n)) / (n-2 )
  t <- matrix(S / (n-1) / (n-2), n, n)
  UA <- A - r - c + t
  diag(UA) <- 0
  return(UA)
}



CMDCH<- function(Y,X1,X2) {
  X1 <- as.matrix(X1)
  X2 <- as.matrix(X2)
  Y  <- as.matrix(Y)
  n <- nrow(X1)
  if (n != nrow(Y)) {
    stop("The dimensions of X and Y do not agree.")
  }
  A <- uYX1.center.H(Y,X1)
  B <- uX2.center.H(X2)
  CMDC.H <- u.inner.H(A, B) / sqrt(u.inner.H(A, A) * u.inner.H(B, B))
  return(CMDC.H)
}





#############
####This R code is to calculate MCH with bandwidth of sample variance of (X) multiply by 2

library(kernlab)#calculate kernel matrix

# inner product of matrix Y and X
u.inner <- function(X, Y){
  n <- dim(X)[1]
  ip <- sum(X * Y) / n / (n - 3) 
  return(ip)
}

#  matrix for Y
uY.center <- function(Y) {
  n <- dim(Y)[1]
  A <- as.matrix((Y%*%t(Y)))
  diag(A) <- 0
  R <- rowSums(A)
  C <- colSums(A)
  S <- sum(A)
  r <- matrix(rep(R, n), n, n) / (n - 2)
  c <- t(matrix(rep(C, n), n, n)) / (n - 2)
  t <- matrix(S / (n - 1) / (n - 2), n, n)
  UA <- A - r - c + t
  diag(UA) <- 0
  return(UA)
}


#matrix for X
uX.center <- function(X) {
  n <- length(X)
  sigma1<-cov(X)
  rbf1<-function(x,y){
    return(exp(-0.5*t(x-y)%*% solve(sigma1)%*%(x-y)))
  }
  class(rbf1)<-"kernel"
  A <- kernelMatrix(rbf1,X)
  diag(A)=0
  R <- rowSums(A)
  C <- colSums(A)
  S <- sum(A)
  r <- matrix(rep(R, n), n, n) / (n - 2)
  c <- t(matrix(rep(C, n), n, n)) / (n - 2)
  t <- matrix(S / (n - 1) / (n - 2), n, n)
  UA <- A - r - c + t
  diag(UA) <- 0
  return(UA)
}


MCH <- function(X, Y) {
  X <- as.matrix(X)
  Y <- as.matrix(Y)
  n <- nrow(X)
  if (n != nrow(Y)) {
    stop("The dimensions of X and Y do not agree.")
  }
  A <- uX.center(X)
  B <- uY.center(Y)
  MC_H <- u.inner(A, B) / sqrt(u.inner(A, A) * u.inner(B, B))
  return(MC_H)
}