# ******************************************************************
#  R Script "Geometric Typicality Test" (January 3, 2019)
#  see Chapter 4, §4.2.3, 4.2.4, 4.3.2 (Target) and 4.4.1 (Student),
#  in                                      
#    "Combinatorial Inference in Geometric Data Analysis"   
#           B. Le Roux, S. Bienaise, J.-L. Durand          
#                  Chapman & Hall/CRC, 2019                 
# *****************************************************************

# === Step 1. The data  ===
  
setwd("D:/path")                         # replace "D:/path" by the path of the data directory
base <-  read.table(file = "Target.txt", header= TRUE,
                    sep= ";", dec= ".", row.names= 1)

# base <-  read.table(file = "Target_X.txt", header= TRUE, sep= ";", dec= ".", row.names= 1)
# base <-  read.table(file = "Student.txt", header= TRUE, sep= ";", dec= ".", row.names= 1)

n     <- dim(base)[1]       # number of points
K     <- dim(base)[2]       # number of variables
X.IK  <- as.matrix(base, nrow= n, ncol= K) #matrix of coordinates
X_P   <- rep(0,K)            # coordinates of reference point

# === Step 2. === Parameters of the method
notable_D  <- 0.4            # notable limit for D
alpha      <- 0.05           # alpha level 
max_number <- 1e+06          # maximum number of samples
n_dir      <- 200            # number of lines for compatibility region
seed       <- NULL           # seed for random number generator

# === Step 3. Descriptive analysis ===
Mcov.KK <- cov.wt(X.IK, method = "ML")$cov
eig <- eigen(Mcov.KK, symmetric = TRUE)
L <- sum(eig$values > 1.5e-8) ; lambda.L <- eig$values[1:L] 
BasisChange.KL <- eig$vectors[ ,1:L] %*% diag(1/sqrt(lambda.L), nrow=L)
Z_GM.IL  <- sweep(X.IK, 2, colMeans(X.IK), "-") %*% BasisChange.KL
Z_GP.L   <-  t(BasisChange.KL) %*% t(t(colMeans(X.IK) - X_P))
norm2_PG <- sum(Z_GP.L^2)    # squared M-distance PG
cat(" M-distance D = ", round(sqrt(norm2_PG),3), "\n",
    " Descriptively, the deviation is ",
    ifelse(norm2_PG >= notable_D^2, "notable.",
           "small: there is no point in doing the test."), sep= "")

# === Step 4. The testing procedure ===
if (2^(n-1) <= max_number){
  cardJ <- 2^(n-1)                    
  Epsilon.JI <- matrix(1L, nrow= cardJ, ncol= n)
  for(i in 1:(n-1))
    Epsilon.JI[ , i+1] <- rep(c(1L,-1L), each=2^(n-i-1), times=2^(i-1))
  Epsilon.J <- t(t(rowSums(Epsilon.JI)))/n
  Z_GGj.JL  <- Epsilon.JI %*% Z_GM.IL/n
  rm(Epsilon.JI)
} else {
  cardJ     <-  max_number; set.seed(seed)
  Epsilon.J <- matrix(0L, nrow= cardJ, ncol= 1)
  Z_GGj.JL  <- matrix(0,  nrow= cardJ, ncol= L)
  for (j in 1:cardJ) {
    epsilon_j.I   <- t(sample(c(-1L, 1L), size= n, replace= TRUE))
    Epsilon.J[j]  <- sum(epsilon_j.I)/n
    Z_GGj.JL[j, ] <- epsilon_j.I %*% Z_GM.IL/n
  }
}
if (norm2_PG >= notable_D^2) {
  Z_PPj.JL <- Z_GGj.JL + Epsilon.J %*% t(Z_GP.L)
  norm2_PPj.J <- rowSums(Z_PPj.JL^2)
  d2P.J <- norm2_PPj.J - (Z_PPj.JL %*% Z_GP.L)^2 / (1 + norm2_PG)
  rm(Z_PPj.JL)
  d2P_obs     <- norm2_PG/(1 + norm2_PG)
  n_sup       <- sum(d2P.J >= d2P_obs * (1 - 1e-12))
  p_value     <- n_sup/cardJ
  if(K>1){
    cat(" p-value = ", 2*n_sup, "/", 2*cardJ, " = ",
        format(round(p_value, 3), nsmall= 3),"\n", sep= "")
   } else {
    cat(" p-value = ", n_sup, "/", 2*cardJ, " = ",
        format(round(p_value/2, 3), nsmall= 3)," (one-sided) \n", sep= "") 
  }
}

# === Step 5. The compatibility region ===
set.seed(seed); rank_inf <- trunc(alpha * cardJ) + 1
# cat(rank_inf)
C.J <- t(t(rowSums(Z_GGj.JL^2)))                 
Anul <- which(abs(C.J + Epsilon.J^2 - 1)  < 1e-12)
if ( K== 1) n_dir <- 1

kappa.D <- rep(0, 2 * n_dir)
for (d in 1: n_dir){
  xy.L <- runif(L, -1, 1) ; beta.L <- t(t(xy.L / sqrt(sum(xy.L^2))))
  U.J <- Z_GGj.JL %*% beta.L
  A.J <- C.J - U.J^2 + Epsilon.J^2- 1
  B.J <- -Epsilon.J * U.J
  DeltaP.J <- abs(B.J^2 - A.J*C.J) 
  X1.J <- (-B.J + sqrt(DeltaP.J))/A.J
  X2.J <- (-B.J - sqrt(DeltaP.J))/A.J
  X1.J[Anul] <- -Inf ; X2.J[Anul] <- Inf
  
  X1_sorted  <- sort(X1.J); X2_sorted  <- sort(X2.J)
  kappa.D[2*d - 1] <- abs(X1_sorted[rank_inf])
  kappa.D[2*d] <- X2_sorted[cardJ + 1 - rank_inf]
}
if (max(kappa.D) == Inf) {
  warning("Finite compatibility region is not accessible")
} else {
  if (L > 1){
    cat("\n Adjusted ", 100*(1 - alpha), "% compatibility region:\n", 
        "  principal ellipsoid of the cloud with scale parameter", 
        " kappa = ", round(mean(kappa.D), 3), " \n", 
        "  (mean of ", 2*max(L,n_dir)," kappa values in the range ", 
        floor(min(kappa.D)*1000)/1000, " to ", 
        ceiling(max(kappa.D)*1000)/1000, "). \n", sep= "")
  } else { 
    x1 <- colMeans(X.IK) -  beta.L*kappa.D[1]*sqrt(Mcov.KK[1,1])
    x2 <- colMeans(X.IK) +  beta.L*kappa.D[2]*sqrt(Mcov.KK[1,1])
    cat("\n ", 100*(1 - alpha), "% compatibility interval [",
        round( min(x1,x2),3)," ; ", round( max(x1,x2), 3),"]", sep= "")
  }
}
#==================== end of the script ====================