# CIGDA_Geo-v1 (January 21, 2019)
# B. Le Roux, S. Bienaise, J.-L. Durand

#### SPAD Interface  ####
#!SPADR Geometric Typicality Test

#SPAD.VARN STATUT = SEL; LABEL = Selected; ROLE = C
#SPAD.VAR1 STATUT = GRP; LABEL = Factor defining the group; ROLE = N

#SPAD.TITRE[1]      LABEL = Analyses
#SPAD.SEP[1]        LABEL = Group definition
#SPAD.PAR.LISTE[1]  NOM = CASES; LABEL =; ITEMS = all cases#selection of categories (put & between categories to be grouped); STYLE = LARGE2; DEFAUT = 1
#SPAD.PAR.TEXTE[1]  NOM = GROUP_LEVEL; LABEL = Name of the group; DEFAUT = <Name of the group>
#SPAD.SEP[1]        LABEL = Reference point
#SPAD.PAR.LISTE[1]  NOM = REF; LABEL = ; ITEMS = origin#last line of the dataset#coordinates (with a semicolon as separator); STYLE = LARGE1; DEFAUT = 1
#SPAD.PAR.TEXTE[1]  NOM = REF_TXT; LABEL = coordinates; DEFAUT = <one value per selected variable>

#SPAD.SEP[1]        LABEL = Choice of analyses
#SPAD.PAR.LISTE[1]  NOM = SPACE_TYPE; LABEL =; ITEMS = one-dimensional#multidimensional#one & multidimensional; STYLE = LARGE2; DEFAUT = 2

#SPAD.SEP[1]        LABEL = Combinatorial procedures
#SPAD.PAR.LISTE[1]  NOM = PROC_TYPE; LABEL = ; ITEMS = test#compatibility region#test & region; STYLE = LARGE2; DEFAUT = 1


#SPAD.TITRE[2]      LABEL = Parameters
#SPAD.SEP[2]        LABEL = Notable limit
#SPAD.PAR.REEL[2]   NOM = NOTABLE; LABEL = M-distance/scaled deviation; DEFAUT = 0.4; MIN = 0

#SPAD.SEP[2]        LABEL = Alpha-levels
#SPAD.PAR.REEL[2]   NOM = ALPHA_T; LABEL = test; MIN = 0.001; MAX = 0.5; DEFAUT = 0.05
#SPAD.PAR.REEL[2]   NOM = ALPHA_C; LABEL = compatibility region; MIN = 0.001; MAX = 0.5; DEFAUT = 0.05

#SPAD.SEP[2]        LABEL = Permutations
#SPAD.PAR.ENTIER[2] NOM = MAX_NUMBER; LABEL = maximum number; MIN = 5000; DEFAUT = 100000

#SPAD.SEP[2]        LABEL = Monte Carlo
#SPAD.PAR.LISTE[2]  NOM = RANDOM_SEED; LABEL = seed; STYLE = LARGE2; ITEMS = random#fixed; DEFAUT = 1
#SPAD.PAR.ENTIER[2] NOM = SEED; LABEL = seed (if fixed); DEFAUT = 1234

#SPAD.SEP[2]        LABEL = Adjusted compatibility region
#SPAD.PAR.ENTIER[2] NOM = NDIR; LABEL = number of lines; MIN = 50; DEFAUT = 100; STEP = 10


#SPAD.TITRE[3]      LABEL = Outputs
#SPAD.SEP[3]        LABEL = Number of decimals
#SPAD.PAR.ENTIER[3] NOM = M_DIGITS; LABEL = mean; MIN = 0; MAX = 8; DEFAUT = 3
#SPAD.PAR.ENTIER[3] NOM = K_DIGITS; LABEL = M-distance & scale parameter; MIN = 1; MAX = 8; DEFAUT = 2
#SPAD.PAR.ENTIER[3] NOM = P_DIGITS; LABEL = p-value; MIN = 2; MAX = 8; DEFAUT = 3

#SPAD.SEP[3]        LABEL = Approximate results
#SPAD.PAR.LISTE[3]  NOM = APPROXIMATE; LABEL = ; ITEMS = yes#no; STYLE = LARGE2; DEFAUT = 2

#SPAD.SEP[3]        LABEL = Export of test statistic distribution
#SPAD.PAR.LISTE[3]  NOM = EXPORT; LABEL = ; STYLE = LARGE2; ITEMS = yes#no; DEFAUT = 2

#### DATA AND PARAMETERS ####
base  = as.data.frame(SPAD.getCategoryAsDataFrame("SEL"))  # numerical
group = SPAD.getCategoryAsDataFrame("GRP")  # categorical

all_cases = (SPAD.getInputParameter("CASES") == 1) 
group_lbl = SPAD.getInputParameter("GROUP_LEVEL")
ref_zero  = (SPAD.getInputParameter("REF") == 1)
ref_last  = (SPAD.getInputParameter("REF") == 2)
ref_user  = (SPAD.getInputParameter("REF") == 3)
ref_txt   = SPAD.getInputParameter("REF_TXT")
test_dim  = SPAD.getInputParameter("SPACE_TYPE")
unidim    = (1 %in% test_dim | 3 %in% test_dim)
multidim  = (2 %in% test_dim | 3 %in% test_dim)
proc_type = SPAD.getInputParameter("PROC_TYPE")
test      = (1 %in% proc_type | 3 %in% proc_type)
region    = (2 %in% proc_type | 3 %in% proc_type)

notable_D = SPAD.getInputParameter("NOTABLE")
alpha_t = SPAD.getInputParameter("ALPHA_T")  # alpha level for test
alpha_c = SPAD.getInputParameter("ALPHA_C")  # alpha level for compatibility region
max_number = SPAD.getInputParameter("MAX_NUMBER")      # number of possible clouds
n_dir        = SPAD.getInputParameter("NDIR") # number of pseudo-random directions
random_seed = (1 %in% SPAD.getInputParameter("RANDOM_SEED"))  # TRUE if random seed is used for each test
seed      = SPAD.getInputParameter("SEED")
if (random_seed) {seed <- NULL}

m_digits    = SPAD.getInputParameter("M_DIGITS")    # Nb of decimal places for mean or sd
p_digits    = SPAD.getInputParameter("P_DIGITS")    # Nb of decimal places for p-values
k_digits    = SPAD.getInputParameter("K_DIGITS")  # Nb of decimal places for M-distance & kappa (scale parameter)
approximate = SPAD.getInputParameter("APPROXIMATE") == 1
export    = NULL  # dataframe for results export
export_TF = SPAD.getInputParameter("EXPORT") == 1

# *** The 3 following lines are temporary functions ***
# setwd("D:/CIGDA/data")
# save.image()
# load("D:/CIGDA/data/.RData")

# technical parameter
iota <- 1e-12 
time1 <- Sys.time()
continue <- TRUE

cat(" ******************************************************** \n",
    " *              GEOMETRIC TYPICALITY TEST               * \n", 
    " *        R script interfaced with Coheris SPAD         * \n", 
    " *                        * * *                         * \n",
    " *              See Chapter 4 of the book:              * \n",
    " * \"Combinatorial Inference in Geometric Data Analysis\" * \n", 
    " *         B. Le Roux, S. Bienaise, J.-L. Durand        * \n",
    " *               Chapman & Hall/CRC, 2019               * \n", 
    " ******************************************************** \n", 
    sep = "")

if(continue & any(is.na(base))) {
  cat("\n INPUT ERROR: there are missing data. \n")
  continue <- FALSE
}
if(continue & !all_cases & any(is.na(group))) {
  cat("INPUT ERROR: ",
      "you must choose Group ('Variables' dialog box). \n", sep = "")
  continue <- FALSE
}

# === group definition ===
txt_tmp <- group_lbl
# split chain into character vector (split = ",")
# and remove leading and trailing whitespace from character strings
selected_levels <- trimws(strsplit(x = txt_tmp, split = "&", fixed = TRUE)[[1]])
groupName <- paste0(selected_levels, collapse = " & ")
rm(txt_tmp)
tmp_levels <- setdiff(selected_levels, unlist(group))

#  ERROR MESSAGE PRINTING
if (continue & !all_cases) {
  if (length(tmp_levels) > 0) {
    cat("Selected level(s) for group \n", sep = "")
    
    if (length(tmp_levels) == 1) {
      cat("INPUT ERROR \n ", tmp_levels, ": category not observed. \n", sep ="")
    } else {
      cat("INPUT ERROR \n ", tmp_levels, ": categories not observed. \n", sep = "")
    } 
    continue <- FALSE
  }
  
  if (length(selected_levels) == 0 & continue) {
    cat("INPUT ERROR \n", 
        " No Group category is defined. \n",
        sep = "")
    continue <- FALSE
  }
  group_TF <- unlist(group) %in% selected_levels
  base <- base[group_TF, ]
}
# === End of group definition ===

# === control of base ===
if(continue) {
  K <- dim(base)[2]
  if (unidim & multidim & K == 1) { multidim <- FALSE}
  if (!unidim & multidim & K == 1) {
    cat("\n SETTING ERROR: only one variable is selected, multidimensional analysis is irrelevant. \n")
    continue <- FALSE
  }
}

# === Control of reference point ===
if (ref_user) {
  # remove spaces
  if (grepl(pattern= " ", x= ref_txt, fixed= TRUE)) {
    while (grepl(pattern= " ", x= ref_txt, fixed= TRUE)) {
      ref_txt <- sub(pattern= " ", replacement= "", x= ref_txt, fixed= TRUE)
    }
  }
  
  # split chain into character vector (split= ";")
  txt_tmp <- strsplit(x= ref_txt, split= ";", fixed= TRUE) [[1]]
  
  if (length(txt_tmp) == 0) {
    continue <- FALSE
    cat("ERROR \n ",
        "Reference mean point coordinate",
        ifelse (K == 1, " is", "s are"),
        " required. \n", sep= "")
  } else {
    for (i in 1:length(txt_tmp)) {
      if (is.na(as.numeric(txt_tmp[i]))) {
        continue <- FALSE
        cat("\n INPUT ERROR: ", txt_tmp[i], 
        " is not a valid coordinate for the reference point. \n", sep="")
      }
    }
    num_tmp <- as.numeric(txt_tmp)
  }
  
  if (continue) {
    if (length(num_tmp) == K) {
      ref.K <- num_tmp
    } else {
      continue <- FALSE
      cat("\n INPUT ERROR: the number of reference point coordinates \n", 
          "is not equal to the number of selected variables. \n")
    }
  }
}

if (continue) {
  
  if ( ref_zero ) {   # O is the reference point P
    X_OM.IK <- as.matrix(base)
    X_OP.K <- matrix(0, ncol= 1, nrow= K)
  }
  if ( ref_last ) {   # the last line give coordinates of reference point P
    X_OM.IK <- as.matrix(base[-dim(base)[1], ])
    X_OP.K <- t(base[dim(base)[1], ])
  }
  if ( ref_user ) {   # coordinates of P are given by user
    X_OM.IK <- as.matrix(base)
    X_OP.K <- matrix(ref.K, ncol= 1)
  }
  n <- dim(X_OM.IK)[1]  # sample size
  
  # Choice of method: exhaustive or Monte Carlo
  cardJ <- min(2^(n - 1), max_number)        # card_J: number of permutation clouds
  if (2^(n - 1) <= max_number) { 
    MC <- FALSE
    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)))
  } else { MC <- TRUE }

  if ( unidim ) {
  	first_variable <- TRUE
  	export <- NULL
    
  	# Permutation table and test statistic
  	if(MC) {
  	  Epsilon.J <-  matrix(rep(0L, cardJ), ncol= 1)
  	  Epsilon.J[1, 1] <- n
  	  set.seed(seed)
  	  E_u.JK <- X.JK <- matrix(0, nrow= cardJ, ncol= K)
  	  for (j in 2:(cardJ)) {
  	    epsilon_j.I <- matrix(sample(x= c(-1L, 1L), size= n, replace= TRUE), nrow= 1, ncol= n)
  	    Epsilon.J[j, 1] <- sum(epsilon_j.I)
  	    X.JK[j, ] <- epsilon_j.I %*% X_OM.IK
  	  }
  	  Epsilon_bar.J <- Epsilon.J/n
  	  E_u.JK <- (X.JK - Epsilon.J %*% t(X_OP.K)) / n
  	} else {
  	  X.JK <- Epsilon.JI %*% X_OM.IK
  	  E_u.JK <- (X.JK - Epsilon.J %*% t(X_OP.K)) / n 
  	}
 
  	for ( icol in 1:K ) {
      # === DESCRIPTIVE ANALYSIS ===
  	  x.I <- as.vector(X_OM.IK[, icol])
  	  u <- X_OP.K[icol]
  		xbar <- mean(x.I)              # observed mean
  		v <- sum(x.I^2)/n - xbar^2     # variance
  		d_obs <- (xbar - u) / sqrt(v)  # scaled deviation
  		notable <- abs(d_obs) >= notable_D
  		
      # === INDUCTIVE ANALYSIS ===
      # Geometric typicality test
  		if (test & notable) {
  		  n_obs <- sum(abs(E_u.JK[, icol]) >= abs(xbar - u) * (1 - iota))
  		  p_value <- n_obs / (2*cardJ)  # one-sided
   		}
    
      # Compatibility interval
  		if (region) {
        u1.J <- as.vector((sum(x.I) - X.JK[ ,icol])/(n - Epsilon.J))
        u2.J <- as.vector((sum(x.I) + X.JK[ ,icol])/(n + Epsilon.J))
        u.2J <- c(u1.J, u2.J)
        u_sorted.2J <- sort(u.2J[-1])
        rank_inf <- trunc(cardJ * alpha_c) + 1  # trunc(2*cardJ * alpha_c/2)
        if (rank_inf < 2) {
          warning_CI <- TRUE
        } else {
          warning_CI <- FALSE
          u_min <- u_sorted.2J[rank_inf - 1]
          u_max <- u_sorted.2J[2*cardJ + 1 - rank_inf]
      		d_min <- (u_min - u) / sqrt(v) ; d_max <- (u_max - u) / sqrt(v)
    		}
  		}

  		if (approximate) {
  		  # (Z test)
    		z_obs <- sqrt(n) * (xbar - u) / sqrt(v + (xbar - u)^2)
    		app_p_value <- pnorm(-abs(z_obs))  # one-sided p-value
    		# compatibility interval
    		z_alpha <- qnorm( 1 - alpha_c / 2)
    		if (z_alpha^2 < n) {
    		  h_alpha <- z_alpha * sqrt(v / (n - z_alpha^2))
    		} else {
    		  h_alpha <- Inf
    		}
    		app_u_max <- xbar + h_alpha
    		app_u_min <- xbar - h_alpha
    		app_d_min <- (app_u_min - u) / sqrt(v)
    		app_d_max <- (app_u_max - u) / sqrt(v)
  		}
  
      # ===  PRINTING the RESULTS ===
  		if (first_variable) {
  		  cat("\n")
  		  cat("=========================================================== \n")
#  		  cat(" Comparison of the mean of a variable to a reference value  \n")
  		  cat(" COMPARISON OF THE MEAN OF A VARIABLE TO A REFERENCE VALUE  \n")
  		  cat("=========================================================== \n")
  		  if (K >= 2) cat("  ", K," variables are analysed separately. \n", sep= "")
  		}
  		cat("\n ********** VARIABLE: ", colnames(base)[icol], " **********\n", sep= "")
  		cat("  size of the group of observations")
  		if (!all_cases) {
  		  cat(" (", groupName, ")", sep = "")
  		}
      cat(" = ", n, "\n", sep= "")
  		cat("  reference value: ", u,"\n\n", sep= "")    
  		
  		cat("Descriptive analysis \n")
  		cat("-------------------- \n")
			cat(
				"  observed mean: ", format(round(xbar, m_digits), nsmall=m_digits),
				"; standard deviation: ", format(round(sqrt(v), m_digits), nsmall=m_digits),
				"\n", sep= ""
			)
			cat("  scaled deviation from reference value to observed mean: d_obs = ", 
			format(round(d_obs, k_digits), nsmall=k_digits), "\n", sep= "")
			cat("  |d_obs| ",
			    ifelse(abs(d_obs) >= notable_D, "> ", "< "),
			    notable_D, ": descriptively, the deviation is ",
			    ifelse(abs(d_obs) >= notable_D, "notable", "small"),
			    ". \n\n", sep= "")

			cat("Combinatorial inference \n")
			cat("----------------------- \n")
			if (MC) {
			      cat("  Number of possible clouds", sep = "")
			  if (2^n < Inf) {
			    if (log10(2^n) < 10) {
			      cat(": ", 2^n, ";\n", sep = "")
			    } else {
			      cat(": ", format(2^n, digits = 1), ";\n", sep = "")
			    }
			    cat("  MC method is performed with 2x", cardJ, "possible clouds, using symmetry property.\n")
			  } else {
			    cat(" > ", format(.Machine$double.xmax, digits = 1), ") \n", sep = "")
			  }
			} else {
			        cat("  Number of possible clouds: ", 2^n, ";\n", sep = "")
					    cat("  exhaustive method is performed.\n", sep = "") 
			}
			
      if (test) {
        cat("\n", "~ Geometric typicality test \n", sep= "")
  			if (!notable) {
  				cat("  Deviation is small: there is no point in performing the test. \n")
  			} else {
  			  cat("  test statistic: deviation from reference value to mean\n\n", sep= "")
  			  side <- ifelse(d_obs > 0, "high", "low")
  			  cat("  p-value: p = ", n_obs, "/", 2*cardJ, " = ", 
  			      format(round(p_value, p_digits), nsmall= p_digits), sep= "")
  			  if (p_value <= alpha_t / 2) {
  			    cat(" < ", alpha_t / 2, " (one-sided) \n", 
  			        "  The observed mean is atypical of the reference value on the side of ", 
  			        side, " values, at level ", alpha_t/2, ". \n", sep= "")
  			  } else {
  			    cat(" > ", alpha_t / 2, " (one-sided) \n", 
  			        "  The observed mean is not atypical of the reference value, at level ", 
  			        alpha_t, ". \n", sep= "")
  			  }
  			}
      }

			# Compatibility interval
			if (region) {
  			if (warning_CI) {
  			  cat("\n", "~ ", "No ", 100 * (1 - alpha_c), "% compatibility interval available \n", sep = "") # à améliorer
  			} else {
    			cat("\n", "~ ", 100 * (1 - alpha_c), "% compatibility interval \n",
    			    "    for Mean: [",
    			    format(round(u_min, m_digits), nsmall= m_digits), ", ",
    			    format(round(u_max, m_digits), nsmall= m_digits),  "] \n", 
    			    "    for Scaled deviation: [",
    			    format(round(d_min, k_digits), nsmall= k_digits), ", ",
    			    format(round(d_max, k_digits), nsmall= k_digits), "] \n", sep= "")
        }
			}
			
			# APPROXIMATE RESULTS
			if (approximate) {
			  if ((test & notable) | region) {
			    cat("-------------- \n", sep = "")
			  }
			  if (test & notable) {
			    cat("  approximate test: z = ", 
			        format(round(z_obs, 3), nsmall= 3), 
			        "; p = ", 
			        format(round(app_p_value, p_digits), nsmall= p_digits),
			        " (one-sided) \n", sep= "")
			  }
    		if (region) {
    			# Compatibility interval
    			cat("  approximate ", 100 * (1 - alpha_c), "% compatibility interval \n", sep= "")
    			cat("    for Mean: [", 
    			    format(round(app_u_min, m_digits), nsmall= m_digits), ", ",
    			    format(round(app_u_max, m_digits), nsmall= m_digits),  "] \n",
    			    "    for Scaled deviation: [",
    			    format(round(app_d_min, k_digits), nsmall= k_digits), ", ", 
    			    format(round(app_d_max, k_digits), nsmall= k_digits), "] \n", sep= "")
			  }
			}
			
			# List of the values of statistics E_u (for export tp SPAD)
			if (export_TF) {
				if (first_variable == TRUE) {
					export <- as.data.frame(c(E_u.JK[ , icol], rev(-E_u.JK[ , icol])))
				} else {
					export <- cbind(export, c(E_u.JK[ , icol], rev(-E_u.JK[ , icol])))
				}
			  colnames(export)[dim(export)[2]] <- paste("E_", 
			                                            colnames(base)[icol], sep= "")
			}
			if (K > 1) {first_variable <- FALSE}
  	}
    
  	# === EXPORT TO SPAD ===
  	if (export_TF && unidim && !(multidim)) {
  	  SPAD.setDataFrame(export, generateColId = FALSE)
  	  cat("\n", "Test statistic distribution",
      ifelse(K == 1, " is", "s are"),
      " exported to SPAD. \n", sep= "")
  	}
  }
  
  if ( multidim ) {
    
    # === DESCRIPTIVE ANALYSIS ===
    X_OG.K <- t(t(colMeans(X_OM.IK)))
  	V <- cov.wt(X_OM.IK, method="ML")$cov   # covariance matrix
  
    # Passage matrix from initial basis to principal orthocalibrated basis
  	E <- eigen(V, symmetric= TRUE)
  	L <- sum(E$values > iota)        # dimension of the cloud
  	BasisChange <- as.matrix(E$vectors[ ,1:L]) %*% diag(1/sqrt(E$values[1:L]))
  	# rm(E)
  
  	Z_GM.IL <- scale(X_OM.IK, center= TRUE, scale= FALSE) %*% BasisChange
  	Z_PG.L <- t(BasisChange) %*% (X_OG.K - X_OP.K)  # coord of PG
  	norm2_PG <- as.vector(t(Z_PG.L) %*% Z_PG.L) # squared M-norm of PG
  	notable <- norm2_PG >= notable_D^2
  	d2P_obs <- norm2_PG/(1 + norm2_PG)

  	# === INDUCTIVE ANALYSIS ===
  	if (MC) {
  	  Epsilon.J <- matrix(0, ncol= 1, nrow= cardJ)
  	  set.seed(seed)
  	  SZ_GGj.JL <- matrix(0, ncol= L, nrow= cardJ)
  	  for (j in 1:(cardJ)) {
  	    epsilon_j.I <- t(sample(c(-1L, 1L), size= n, replace= TRUE))
  	    SZ_GGj.JL[j, ] <- epsilon_j.I %*% Z_GM.IL
  	    Epsilon.J[j] <- sum(epsilon_j.I)
  	  }
  	} else {
  	  SZ_GGj.JL <- Epsilon.JI %*% Z_GM.IL
  	  rm(Epsilon.JI)
  	}
  	SZ_PPj.JL <- SZ_GGj.JL + Epsilon.J %*% t(Z_PG.L)
  	norm2_PPj.J <- rowSums(SZ_PPj.JL^2)/n^2
  	d2P.J <- norm2_PPj.J - (SZ_PPj.JL %*% Z_PG.L)^2 / (n^2 *(1 + norm2_PG))

  	if (test & notable) {
  		n_sup <- sum(d2P.J >= (d2P_obs * (1 - iota) - iota))
  		p_value <- n_sup/cardJ
  	}

		if (region) {
    	# pseudo-random directions in the space
  		C.J <- t(t(rowSums(SZ_GGj.JL^2)))          # Calcul de n^2 x |GGj|^2
  		Anul <- which(abs(C.J + Epsilon.J^2 - n^2)  < iota )
  		rank_inf <- trunc(alpha_c * cardJ) + 1
  	 	set.seed(seed)
    	kappa.D <- rep(0, 2 * max(L, n_dir))
      # L principal axes and n_dir - L random lines 
    	for (d in 1:max(L, n_dir)) {
    	  if (d <= L){
    	    SU.J <- SZ_GGj.JL[ ,d]
    	  } else {
    	    xy.L <- runif(L, -1, 1)
    	    beta.L <- t(t(xy.L / sqrt(sum(xy.L^2))))  # random unit vector
    	    SU.J <- SZ_GGj.JL %*% beta.L
    	  }
        # determines x1 and x2 in direction of GP 
    	  A.J <- C.J - SU.J^2 + Epsilon.J^2 - n^2
    		B.J <- -Epsilon.J * SU.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)
    		kappa.D[2*d - 1] <- abs(X1_sorted[rank_inf])
    		X2_sorted  <- sort(X2.J)
    		kappa.D[2*d] <- X2_sorted[cardJ + 1 - rank_inf]
    	}
    	kappa <- mean(kappa.D)
    	
    	# Compatibility ellipse
    	if ((L == 2) & (kappa != Inf)) {
    	  nbp <- 720
    	  lambda.L <- E$values[1:2]
    	  M.22 <- diag(sqrt(lambda.L), nrow = L) %*% t(E$vectors[ , 1:2])
    	  Circle <- matrix(0, nrow = nbp, ncol = 2)
    	  for (a in 1:nbp) {
    	    Circle[a, 1] <- cos(a*2*pi/nbp)
    	    Circle[a, 2] <- sin(a*2*pi/nbp)
    	  }
    	  A0 <- kappa * diag(c(-1, -1))
    	  B0 <- kappa * diag(c(1, 1))
    	  A <- t(t(A0 %*% M.22) + as.vector(X_OG.K))
    	  B <- t(t(B0 %*% M.22) + as.vector(X_OG.K))
    	  Ellipse <- sweep(kappa * Circle %*% M.22, X_OG.K, FUN = "+", MARGIN = 2)
    	  
    	  # range of figure
    	  rx <- range(c(base[, 1], Ellipse[ , 1]))
    	  ry <- range(c(base[, 2], Ellipse[ , 2]))
    	  xG <- X_OG.K[1,]
    	  yG <- X_OG.K[2,]
    	  
    	  # for pretty(): same range for x and y
    	  dx <- (rx[2] - rx[1])
    	  dy <- (ry[2] - ry[1])
    	  if(dx > dy) {
    	    dmax <- dx
    	    ry <- ry + c(-(dx - dy)/2, (dx - dy)/2)
    	  } else {
    	    dmax <- dy
    	    rx <- rx + c(-(dy - dx)/2, (dy - dx)/2)
    	  }
    	  
    	  # small expansion of range (20%)
    	  rx <- rx + c(-dmax/10, dmax/10)
    	  ry <- ry + c(-dmax/10, dmax/10)
    	  
    	  SPAD.report.newPage(title = "Compatibility region")
    	  SPAD.report.addText("    (the numerical results and conclusions of tests are in the file 'Edition') ")
    	  plot(base,
    	       xlim = rx, ylim = ry,
    	       asp = 1.18,
    	       col = "grey", cex = 1.5, lwd = 2,
    	       main = substitute(paste(a1, "% compatibility region adjusted by the ",
                    kappa, "-ellipse of observed cloud, with ", kappa, " = ", k1), 
                    list(a1 = 100 * (1 - alpha_c), k1 = round(kappa, k_digits))))
    	  abline(h = pretty(ry, 10), v = pretty(rx, 10), col = "lightgray", lty = 2)
    	  points(xG, yG, cex = 1.5, lwd = 3, col = "blue")  # mean point
    	  text("G", x = xG + (rx[2] - rx[1]) / 25, y = yG, cex = 1.5)
    	  points(Ellipse, asp = 1, cex = 0.1, col = "blue", type = "o")  # ellipse
    	  segments(x0= A[,1], y0= A[,2], x1= B[,1], y1= B[,2], col= "blue", lty= 2)
    	}
    }
  	
    if (approximate) {
    	# Approximate test
    	chi2_obs <- n * d2P_obs
    	app_p_value <- 1 - pchisq(chi2_obs, df= L)
    	# Approximate compatibility ellipsoid
    	x <- qchisq(alpha_c, df= L, lower.tail= FALSE)
    	if (x < n) {
    	  app_kappa <- sqrt(x / (n - x))
    	} else {
    	  app_kappa <- Inf
    	}
    }
  
  	# ===  PRINTING of the RESULTS===
  	{
  	  cat("\n")  
  	  cat("============================================================== \n")
#  	  cat(" Comparison of the mean point of a cloud to a reference point  \n")
  	  cat(" COMPARISON OF THE MEAN MOINT OF A CLOUD TO A REFERENCE POINT  \n")
  	  cat("============================================================== \n")
  	  cat("", K," variables (", paste(colnames(base), collapse = ", "),
  	      ") are analysed jointly; \n")
  	  cat(" cloud is in a  ",L,"-dimensional geometric space.\n", sep = "")
  	  cat("    size of the group of observations")
  	  if (!all_cases) {
  	    cat(" (", groupName, ")", sep = "")
  	  }
  	  cat(" = ", n, "\n", sep= "")
  	  cat("    reference point:     (", sep= "")
  	  for (i in 1: (K-1)) {cat(format(round(X_OP.K[i], m_digits), nsmall=m_digits), "; ", sep= "")}
  	  cat(format(round(X_OP.K[K], m_digits), nsmall=m_digits), ") \n", sep= "")

  	  cat("\nDescriptive analysis \n")
  	  cat("-------------------- \n")
  	  cat("  observed mean point: (", sep= "")
  	  for (i in 1: (K-1)) {cat(format(round(X_OG.K[i], m_digits), nsmall=m_digits), "; ", sep= "")}
      cat(format(round(X_OG.K[K], m_digits), nsmall=m_digits), ") \n", sep= "")
  	  
  	  
  		# conclusion about the magnitude of deviation
  		cat("  M-distance between observed mean point and reference point: ",
  		    format(round(sqrt(norm2_PG), m_digits), nsmall=m_digits), 
  		    ifelse(norm2_PG >= notable_D^2, " > ", " < "),
  		    notable_D, "\n", sep= "")
  		cat("  Descriptively, the deviation is ",
  		    ifelse(norm2_PG >= notable_D^2, "notable.", "small."), 
  		    "\n\n", sep= "")
  		
  		cat("Combinatorial inference \n", sep= "")
  		cat("----------------------- \n", sep= "")
  		
  		if (MC) {
  		  cat("  Number of possible clouds", sep = "")
  		  if (2^n < Inf) {
  		    if (log10(2^n) < 10) {
  		      cat(": ", 2^n, ";\n", sep = "")
  		    } else {
  		      cat(": ", format(2^n, digits = 1), ";\n", sep = "")
  		    }
  		  } else {
  		    cat(" > ", format(.Machine$double.xmax, digits = 1), ";\n", sep = "")
  		  }
  		   cat("  MC method is performed with 2x", cardJ, "possible clouds, using symmetry property.\n")
  		} else {
  		  cat("  Number of possible clouds: ", 2^n, ";\n", sep = "")
  		  cat("  exhaustive method is performed.\n", sep = "")
  		}
      
  		if (test) {
    		if (!notable) {
    			cat("\n", "Deviation is small: there is no point in doing the test. \n")
    		} else {
    		  cat("\n~ Geometric typicality test \n", sep= "")
    		  cat("  test statistic: squared generalized Mahalanobis distance with respect to reference point ",
    		      "between mean point and reference point", " (observed value: ",
    		      format(round(d2P_obs, k_digits), nsmall=k_digits), ")\n\n", sep = "")
      		cat("  p-value = ", 
      			2*n_sup, "/", 2*cardJ, " = ",
      			format(round(p_value, p_digits), nsmal= p_digits),
      			ifelse(p_value <= alpha_t, " < ", " > "), alpha_t,
      			"\n", sep= "")
      		
      		if (p_value <= alpha_t) {
      			cat("  The observed mean point is atypical ",
      				  "of the reference point, at level ", alpha_t, ". \n", sep= "")
      		} else {
      			cat("  The observed mean point is not atypical ",
      				  "of the reference point, at level ", alpha_t, ". \n", sep= "")    
      		}
    		}
  		}
  	  
  		# === COMPATIBILITY REGION ===
  		if (region) {
    		ellips = ifelse(L == 2, "ellipse", "ellipsoid")
    		cat("\n", "~ adjusted ", 100 * (1 - alpha_c), 
    		    "% compatibility region is defined by \n",
    		    "  the principal kappa-", ellips, " of observed cloud, with kappa = ", 
    		    format(round(mean(kappa.D), k_digits), nsmall= k_digits), sep = "")
    		cat("\n  (mean of ", 2*max(L,n_dir)," values in the range ", 
    		    format(round(min(kappa.D), k_digits), nsmall= k_digits), " to ",
    		    format(round(max(kappa.D), k_digits), nsmall= k_digits), ") \n", sep= "")
  		}
  		
  		# Approximate Results
  		if (approximate) {
  		  if ((test & notable) | region) {
  		    cat("-------------- \n", sep = "")
  		  }
  		  if (test & notable) {
    			app_p_value_txt <- format(round(app_p_value, p_digits), nsmal= p_digits)
    			cat("  approximate test: Chi-Square = ", round(chi2_obs, 3), 
    			    ", df = ", L, "; ", sep= "")  
    			cat("p = ", app_p_value_txt, "\n", sep= "")
  		  }
  		  if (region) {
      		cat("  approximate kappa: ", 
      		    format(round(app_kappa, k_digits), nsmall= k_digits), "\n", sep= "")
  		  }
  		}
  	}
  
  	# === EXPORT of DISTRIBUTIONS TO SPAD ===
  	if (export_TF && !(unidim) && multidim) {
  		export <- as.data.frame(c(d2P.J, rev(d2P.J)))
  		colnames(export) <- "D2"
  		SPAD.setDataFrame(export, generateColId = FALSE)
  		cat("\n", "Test statistic distribution is exported to SPAD. \n", sep= "")
  	}
  }
  
  # === Export of one- and multidimensional distributions to SPAD ==
  if ( export_TF && unidim && multidim ) {
  	export <- cbind(export, c(d2P.J, rev(d2P.J)))
  	colnames(export)[dim(export)[2]] <- "D2"
  	SPAD.setDataFrame(export, generateColId = FALSE)
  	cat("\n", "Test statistic distributions are exported to SPAD. \n", sep= "")
  }
  
  if (!export_TF) {cat("\n")}
  time2 <- Sys.time()
  elapsed_time <- round(as.numeric(difftime(time2, time1, units= "secs")), 2)
  sec_s <- if (elapsed_time < 2) {" sec"} else {" secs"}
  cat("Elapsed time: ", elapsed_time, sec_s, "\n", sep= "")
}