# CSA (30th April 2020; Modifs du 25/03/2020 & 26/09/2023)

#!SPADR Class Specific Analysis
#==============================
# ===  Tab "Variables" definition ===
#SPAD.VARN STATUT = MCA; LABEL = active in MCA; ROLE = CAT
#SPAD.VAR1 STATUT = REF; LABEL = Reference; ROLE = CAT
#SPAD.VAR1 STATUT = GRP; LABEL = Class;     ROLE = CAT
#SPAD.VAR1 STATUT = WCASES; LABEL = Weights of cases; ROLE = CONT

#======================= 
# === Tab "Analyses" ===
#SPAD.TITRE[1]     LABEL = Analyses

#SPAD.SEP[1]       LABEL = Full path of the file of categories 
#SPAD.PAR.FICHIER  NOM = FICH ; LABEL = ;LABEL_AFTER = archived after MCA; DEFAUT = <> ; EXT = TXT

#SPAD.SEP[1]       LABEL =  
#SPAD.SEP[1]       LABEL = REFERENCE and CLASS definition (put & between categories to be grouped)

#SPAD.PAR.TEXTE[1] NOM = ACTIVE_LEVEL; LABEL = category(ies); LABEL_AFTER = defining the reference set; DEFAUT = <Name of category>
#SPAD.PAR.TEXTE[1] NOM = CLASS_LEVEL; LABEL = category(ies); LABEL_AFTER = defining the class; DEFAUT = <Name of the class>


#SPAD.SEP[1]        LABEL = CSA PARAMETERS
#SPAD.PAR.ENTIER[1] NOM = NB_AXES_CSA; LABEL = Number of axes of CSA; LABEL_AFTER= to be interpreted; MIN = 1; DEFAUT = 3
#SPAD.PAR.ENTIER[1] NOM = NB_AXES_REF; LABEL =Number of axes of MCA; LABEL_AFTER= retained for comparisons; MIN = 1; MAX = 50; DEFAUT = 10
#SPAD.PAR.LISTE[1]  NOM = CRITERION; LABEL = Criterion for contribution; STYLE = LARGE2; ITEMS = average#fixed; DEFAUT = 1
#SPAD.PAR.REEL[1]   NOM = CRITVALUE; LABEL = if fixed, the value is; LABEL_AFTER=%; MIN = 0; MAX = 100; DEFAUT = 0

#======================
# === Tab "Outputs" ===
#SPAD.TITRE[2]    LABEL = Outputs
#SPAD.SEP[2]      LABEL = Results for cases
#SPAD.PAR.BOOL[2] NOM   = COORDIND; LABEL = Coordinates;  DEFAUT = FALSE
#SPAD.SEP[2]      LABEL = Parameters for  graphs: scaling coefficient
#SPAD.PAR.REEL[2] NOM   = RAYON; LABEL=; LABEL_AFTER = for circles; MIN = 0.01; DEFAULT = 2
#SPAD.PAR.REEL[2] NOM   = TEXTSIZE; LABEL=; LABEL_AFTER = for font size of category labels; MIN = 0.01; DEFAULT = 0.9; MAX = 3.0
#SPAD.PAR.REEL[2] NOM   = TEXTWIDTH;LABEL=; LABEL_AFTER = for width of category labels; MIN = 0.01; DEFAULT = 0.2; MAX = 1

#### DATA AND PARAMETERS ####
baseMCA    = SPAD.getCategoryAsDataFrame("MCA")  # coord of MCA for reference set
reference  = SPAD.getCategoryAsDataFrame("REF")  # categorical
grp        = SPAD.getCategoryAsDataFrame("GRP")  # categorical
wcases     = SPAD.getCategoryAsDataFrame("WCASES")  # numerical
group_lbl  = SPAD.getInputParameter("CLASS_LEVEL")
active_lbl = SPAD.getInputParameter("ACTIVE_LEVEL")
filename   = SPAD.getInputParameter("FICH")
criterion  = SPAD.getInputParameter("CRITERION")
crit_value = SPAD.getInputParameter("CRITVALUE")
L_ref      = SPAD.getInputParameter("NB_AXES_REF")
L_class    = SPAD.getInputParameter("NB_AXES_CSA")
sortInd    = SPAD.getInputParameter("COORDIND")
echelle    = SPAD.getInputParameter("RAYON")
text_scale = SPAD.getInputParameter("TEXTSIZE")
label_width= SPAD.getInputParameter("TEXTWIDTH")
spad       = TRUE
sep_field  = ";"
dec_code   = "."
##############################

#### SCRIPT R ####
cat(" the working directory is ", getwd(), "\n")  ### modif 26/09/22
# function of binarization of categorical variables (package ade4)
disjonctif = function (df)
{
  df <- as.data.frame(df)
  acm.util.df <- function(i) {
    cl <- df[, i]
    cha <- names(df)[i]
    n <- length(cl)
    cl <- as.factor(cl)
    x <- matrix(0L, n, length(levels(cl)))
    x[(1:n) + n * (unclass(cl) - 1)] <- 1L
    dimnames(x) <- list(row.names(df), paste(cha, levels(cl), sep = ":"))
    return(x)
  }
  if (ncol(df) == 1)
    G <- acm.util.df(1)
  else {
    G <- lapply(1:ncol(df), acm.util.df)
  }
  G <- data.frame(G, check.names = FALSE)
  return(G)
}
##############################
# *** The 3 following lines are temporary functions ***
# setwd("D:/R_Scripts/CSA") 
# save.image()
# load("D:/R_Scripts/CSA/.RData")
############################## 
iota <- 1e-3                                      ### Modif 20/03/21
cat(" The working directory is ", getwd(), "\n") ### modif 26/09/22

#=== verification of the inputs ===
#===----------------------------===

message <- paste("\n\n === ERROR #1 ===\n", 
           "For MCA, you must select at least 2 categorical variables ('Variables' dialog box).\n\n",
             sep = "")
cond <- FALSE
if (!any(is.na(baseMCA))) { cond <- (dim(baseMCA)[2] < 2) }
if ( any(is.na(baseMCA)) | cond) {
  cat(message) ; stop(message)
}
if (any(is.na(grp))) {
  message <- paste("\n === ERROR #2 ===\n",
      "you must choose the variable that defines the class ('Variables' dialog box).\n\n", sep = "")
  cat(message) ; stop(message)
}
cat_group <- trimws(strsplit(x = group_lbl, split = "&", fixed = TRUE)[[1]])
grp_names <- levels(as.factor(grp[, 1]))
# ERROR MESSAGE PRINTING
i_false <- which(!cat_group %in% grp_names)
#names_valid <- paste(c(grp_names, ""), collapse=", ")
names_valid <- paste(grp_names, collapse=", ")
if(length(i_false) == 1){
  message <- paste("\n\n === ERROR #3 ===\n For the definition of class, '", cat_group[i_false],
       "' is not a valid category name\n", " Valid names are: ", names_valid, ".\n\n", sep="" )
  cat(message) ; stop(message)
}
if (length(i_false) > 1){
   names_false <- paste(cat_group[i_false], collapse=", ")
   message <- paste("\n\n === ERROR #3 ===\n For the definition of class, ", names_false,
                     "are not valid category names.\n", " Valid names are: ", names_valid, ".\n", sep="")
   cat(message) ; stop(message)
}
if (length(wcases) == 1 & any(is.na(wcases))) {          #les poids des indivdus sont égaux à 1
  wcases <- as.data.frame(rep(1L,nrow(baseMCA)))
}
if (length(reference) == 1 & any(is.na(reference))){     #all cases are in the reference set
  base <- baseMCA
  group <- as.matrix(grp, ncol=1)
  nI.I  <- wcases[ ,1]
  message <- "\n\n === WARNING=== \n No variable for reference, therefore all cases belong to the reference set.\n"
  warning(message, call. = FALSE)

} else{
  cat_active <- trimws(strsplit(x = active_lbl, split = "&", fixed = TRUE)[[1]])
  ref_names <- levels(as.factor(reference[ ,1]))
  i_false  <- which(!cat_active %in% ref_names)
  names_valid <- paste(ref_names, collapse=", ")
  if(length(i_false) == 1){
    message <- paste("\n\n === ERROR #4 ===\n For the definition of the reference set, '", cat_active[i_false], "' is not valid category name\n",
                     " Valid names are: ", names_valid, ".\n", sep="")
    cat(message) ; stop(message)
  }
  if(length(i_false) > 1){
    names_false <- paste(cat_active[i_false], collapse=", ")
    stop("\n\n === ERROR #4 ===\n For the definition of the reference set, ", names_false, "are not valid category names\n",
         " valid names are: ", names_valid, "\n" )
  }
#=== remove non-active cases from the database === 
  ind_active <- which(reference[ ,1] %in%  cat_active)
  base <- baseMCA[ind_active, ]    #data IxQ for reference set
  group <- matrix(grp[ind_active, ],ncol=1) 
  nI.I  <- wcases[ind_active, ]
}
rm(baseMCA, grp, reference, wcases)

if (filename == "<>" | filename == ""){
  message <- paste("\n === ERROR #5 ===\n No filename is given for category attributes ('Analyses' dialog box).\n\n")
  cat(message) ; stop(message)
}
# Study of the file of categories 
cat_table <-  read.table(filename, sep=";", header=TRUE)        ### Modif 21/03/20
if(dim(cat_table)[2] == 0){
  message <- paste("\n\n === ERROR #6 ===\n The file '",filename,
           "' is wrong: the field separator is not a semicolon.\n\n", sep="")
  cat(message); stop(message, call. = TRUE)
}

numcol_var <- 2
if( !(colnames(cat_table)[2] == "Variable") ) {
  message <- paste("\n\n === ERROR #7 ===", 
          "\n The 'Variable' checkbox \n of the file of category attributes (", filename,
          ") must be selected.\n\n")
  cat(message) ; stop(message)
}
numcol_cat <- 3
message <- paste("\n\n === ERROR #8 ===",
                 "\n The Category/Modalite checkbox \n of the file of category attributes (",
                 filename,") must be selected.\n\n")
if (dim(cat_table)[2] == 3) {cat(message) ; stop(message, call. = TRUE)}
if (dim(cat_table)[2] > 3  & !(colnames(cat_table)[4] %in% c("Status", "Statut"))){
   cat(message) ; stop(message, call. = TRUE)
}    
numcol_status <- 4
if(!(colnames(cat_table)[4] %in% c("Status", "Statut"))) {
  message <- paste("\n\n=== ERROR #9 ===\n Tthe 'Status/Statut' checkbox \n of the file of category attributes (",
                    filename,") must be selected.\n\n")
  cat(message) ; stop(message)
}
numcol_eff <- 5
if ( !(colnames(cat_table)[5] %in% c("Weight", "Poids" )) ) {
  message <- paste("\n\n=== ERROR #10 ===",
                   "\nthe 'Weight/Poids' checkbox \n of the file of category attributes  (",
                   filename,")  must be selected.\n\n")
  cat(message) ; stop(message)
}
numcol_weight <- 6
if ( !(colnames(cat_table)[6] %in% c("Relative.weight","Poids.relatif"))) {
  message <- paste("\n\n===  ERROR #11 ===",
                   "\nthe 'Relative weight/Poids relatif' checkbox \n of the file of category attributes  (", filename,") must be selected.\n\n")
  cat(message) ; stop(message)
}
#=== End of the verification of inputs ===
#========================================

# === Selection of the reference set ===

# Construction of category_table with only categories of active questions
# list_supp: list of categories of supplementary questions
list_supp <- which(cat_table[ ,numcol_status] %in% c("supplementary", "illustrative"))
list_K <- c(1:dim(cat_table)[1])
if (length(list_supp) > 0) {list_K <- list_K[-list_supp]}
category_table <- cat_table[list_K, ]
# re-defining the names of categories in category_table (first column)  #modif 25 Mars 2021
category_table[ , 1] <- paste(category_table[ ,2], category_table[ ,3], sep=":")

#cardQ <- nlevels( factor(category_table$Variable) )
if (dim(base)[2] != nlevels( factor(category_table$Variable) )){
  message <- paste("\n\n===  ERROR #11 ===\nthe number of 'active questions' of MCA", 
                   "selected in the 'variables' checkbox ",
                   "\nis not equal to the one of the MCA. \n\n")
  cat(message) ; stop(message)
}
nntot = sum(category_table[ , numcol_eff])/ncol(base)
if( abs(nntot - sum(nI.I)) > iota ){                              ### Modif 20/03/21
  message <- paste("\n\n===  ERROR #12 ===\nthere are ", sum(nI.I), " active individuals in the CSA script and ", nntot,
                   "in the MCA: \nVERIFY the definition of the reference set ('Variables' and 'Analyses' dialog boxes)\n\n")
  cat(message); stop(message)
}
cat("List of active variables of MCA:\n")
cat(names(base), sep=" ; ")
  
#Z.IrefK <- as.matrix(disjonctif(base) )[, row.names(category_table)]
# row.names(disjonctif(base)) %in% category_table[ ,1]
Z.IrefK <- as.matrix(disjonctif(base) )[, category_table[ ,1]]     ### Modif 21/03/21
cat_status   <- category_table[ ,numcol_status]
label_status <- colnames(category_table)[numcol_status]
#typeof(Z.IrefK) ; object.size(Z.IrefK)

cardI 	<- nrow(base)         # number of cases in reference set
cardQ   <- ncol(base)         # number of active questions
cardK	  <- dim(Z.IrefK)[2]    # number of categories of active questions
ntot    <- sum(nI.I)          # total of absolute frequencies of active cases
rm(base)

cat(cbind("\n\n Number of reference cases:",cardI, " sum of weights: ", ntot))
cat(cbind("\n Number of active questions:",cardQ))
# categories management
list_passive <-which(cat_status=="passive")         #list of passive categories
cardKp  = length(list_passive)                                 
list_active <- c(1: cardK)                          # list of active categories
if (cardKp > 0 ){
  list_active <- list_active[-list_passive]
}

cat(cbind(" ; total number of categories:", cardK, " ; number of passive categories:", cardKp))
if (cardKp>0) {
  cat("\n list of passive categories: \n")
  print(colnames(Z.IrefK)[list_passive], sep =" ; ")
}

#computation of reference absolute (n_ref.K),  relative (freq_ref.K) frequencies and relative weights (weight.K)
#n_ref.K    <- as.integer(apply(diag(nI.I) %*% Z.IrefK,2,sum))
n_ref.K    <- apply(diag(nI.I) %*% Z.IrefK,2,sum)       ### Modif 21/03/21
freq_ref.K <- n_ref.K/ntot
weight.K  <- category_table[ ,numcol_weight]/100 
label_weight <- colnames(category_table)[numcol_weight]
if (cardKp > 0) {weight.K[list_passive] <- 0}
freq_q.K    <- weight.K/freq_ref.K                 #relative weigths from variables

#cat("\nrelative weigths from variables\n")
#cat(freq_q.K )

### selection of cases of the class
list_class = which(group[ ,1] %in% cat_group) #list of cases in CSA (CSA)
n_class.K <- as.integer(apply((diag(nI.I) %*% Z.IrefK)[list_class,],2,sum)) # counts of categories in the class
n_class   <- sum(n_class.K)/cardQ                          # sum of weights of cases in the class
{
  cat("\n\nCategories defining the class: ", cat_group)
  cat(cbind("\n Number of cases in the class:",length(list_class), " sum of weights: ", n_class))
}
#===----------------------------===
#=== Computation of CSA results ===
#===----------------------------===
# matrix of reference set to be diagonalised
S.KK <- diag(sqrt(freq_q.K/n_ref.K)) %*%  
  ( t(Z.IrefK) %*% diag(nI.I) %*% Z.IrefK - (n_ref.K %*% t(n_ref.K)/ntot) ) %*%
  diag(sqrt(freq_q.K/n_ref.K))

var_ref <- sum(diag(S.KK)) # variance of reference cloud (MCA)
eig <- eigen(S.KK[list_active, list_active], symmetric = TRUE)
dimSpace <- sum(eig$value > 1e-8)
if (L_ref > dimSpace){
  L_ref =dimSpace
}
lambda.L <- eig$value[1:L_ref]        ; names(lambda.L) <- paste("Axis_MCA", 1:L_ref, sep=".")
Vec_ref.KL <- eig$vectors[ , 1:L_ref] ; colnames(Vec_ref.KL) <- names(lambda.L)

Yref <- Z.IrefK[ ,list_active]  %*%  diag( sqrt(freq_q.K/freq_ref.K)[list_active]) %*% Vec_ref.KL
Yref.IL <- apply(Yref, 2, function(x) x - weighted.mean(x,nI.I))
Yref.IL[1:10, ]

cat("\n\n")
print("Variances of reference axes (MCA)",quote=FALSE)
print(lambda.L,digits=5)
rm(S.KK, eig)

#matrix of class to be digonalised
S_class.KK <- diag(sqrt(freq_q.K/n_ref.K)) %*%  
  ( t(Z.IrefK[list_class, ]) %*% diag(nI.I[list_class]) %*% Z.IrefK[list_class, ]  - (n_class.K %*% t(n_class.K)/n_class) ) %*%
  diag(sqrt(freq_q.K/n_ref.K)) * ntot/n_class

var_class <- sum(diag(S_class.KK[list_active, list_active]))                 # variance du nuage de l'ACM
eig       <- eigen(S_class.KK[list_active, list_active], symmetric = TRUE)
dimSpace  <- sum(eig$value > 1e-8)
if (L_class > dimSpace){ L_class <- dimSpace }

mu.L      <- eig$value[1:L_class] ; names(mu.L) <- paste("Axis_CSA", 1:L_class, sep=".")
Vec_class.KL <- eig$vectors[ , 1:L_class] ; colnames(Vec_class.KL) <- names(mu.L)
cat("\n\n")
print("Variances of class axes (CSA)",quote=FALSE)
print(mu.L,digits=5)
rm(S_class.KK, eig)

# Coordinates of categories in CSA
Y_csa.KL <- diag(1/sqrt(weight.K[list_active])) %*% Vec_class.KL %*% diag(sqrt(mu.L))
rownames(Y_csa.KL) <- colnames(Z.IrefK)[list_active]  #row.names(category_table)[list_active]
colnames(Y_csa.KL) <- names(mu.L)

# Coordinates of individuals in CSA
if(sortInd){
  Yc <- Z.IrefK[list_class,list_active]  %*%  diag( sqrt(freq_q.K/freq_ref.K)[list_active]) %*% Vec_class.KL
  Yclass.IL <- apply(Yc,2,function(x) x - weighted.mean(x,nI.I[list_class]))
}

# Cosinus of angles between axes of ACM and of CSA
Cos.LL <- round(t(Vec_ref.KL) %*% Vec_class.KL, digits=10)


rm(Vec_ref.KL,Vec_class.KL) 
#--------------------------------------------
# Contributions des modalités actives et passives
Ctr.KL = diag(weight.K[list_active]) %*% Y_csa.KL^2 %*% diag(1/mu.L)
rownames(Ctr.KL) = rownames(Y_csa.KL) ; colnames(Ctr.KL) = c(paste("Ctr",c(1:L_class),sep="_"))
#  Sélection des modalités actives selon la contribution à l'axe (Ctr_axes=1 sinon 2)
crit	= crit_value/100                                     #critère de sélection des modalités
if (criterion == 1) {
  crit	<-  1./(cardK -cardKp)
  crit_value <- 100*crit
}
cat(c("\n criterion: ", round(crit*100, digits=2), "%\n"), sep="")
Ctr_axes.KL <- matrix("no", nrow=cardK-cardKp, ncol=L_class)
colnames(Ctr_axes.KL) <-  c(paste("Crit",c(1:L_class),sep="_"))
for (i in seq(1,L_class) ){
  Ctr_axes.KL[which(Ctr.KL[ ,i] >= crit),i] <- "yes"
}
#export <- cbind.data.frame(category_table[list_active ,numcol_status], weight.K[list_active], Y_csa.KL, Ctr.KL, Ctr_axes.KL)
export <- cbind.data.frame(weight.K[list_active], Y_csa.KL, Ctr.KL*100, Ctr_axes.KL)
colnames(export) <- c(label_weight, names(mu.L), colnames(Ctr.KL), colnames(Ctr_axes.KL))
#

# Exportation of the data to SPAD
if(spad){
  SPAD.setDataFrame(export)
} else {
  file0 <- paste(wdir,"/", doss_name,"_coord_categories.csv", sep="")
  write.table(export, file0, sep = sep_field, dec = dec_code, append= FALSE,
              row.names = TRUE, col.names=NA)
}

###-------------------------------------###
#      Excel report and text results
###------------------------------------####
if(!spad) {
  cat("\n === PRINTING THE RESULTS ===\n")
  sink(file = file_out, append = FALSE, split = FALSE)
}     
#--- Elementary statistics
#  toto <- cbind.data.frame(category_table[ ,numcol_status],100*freq_ref.K,100*n_class.K/n_class)
toto <- cbind.data.frame(cat_status,100*freq_ref.K,100*n_class.K/n_class)
colnames(toto) <- c(label_status, "reference","class")
rownames(toto) <- colnames(Z.IrefK)
texte <- paste(" number of reference cases: ", cardI, "; sum of weights:", ntot,
               "\n number of cases in the class:", length(list_class),  "; sum of weights:", n_class,
               "\n number of questions:", cardQ, "\n total number of categories:", cardK, "; number of passive categories:",cardKp)
if(spad){
  SPAD.report.newPage("Elementary statistics")
  SPAD.report.addText (texte)
  SPAD.report.addTable(toto, title = "Reference frequencies and class frequencies (in %) of categories" ) 
} else {
  cat("\nElementary Statistics\n=====================\n\n", texte,
      "\n\nReference frequencies and class frequencies (in %) of categories\n")
  print(toto, digits = 2)
}
#---Results of CSA 
eigenValues <- t(t(mu.L))
taux <- 100*eigenValues/var_class
Var_Rate <- cbind.data.frame(eigenValues,taux )
colnames(Var_Rate) <- c("variance of axis", "% of variance")
result_active <- cbind.data.frame(100*weight.K[list_active], Y_csa.KL, 100.*Ctr.KL)
colnames(result_active) <- c(label_weight, names(mu.L), colnames(Ctr.KL))
if(spad){
  SPAD.report.newPage("Results of CSA") 
  SPAD.report.addTable(Var_Rate, title = paste("Variance of subcloud: ", round(var_class, digits=3)))
  SPAD.report.addTable(result_active, title = "Active categories: relative weights (in %), coordinates, contributions (in %) to axes" )
} else{
  cat("\nResults of CSA\n==============\n\n ")
  cat(" Variance of cloud: ", round(var_class, digits=3),"\n")
  print(Var_Rate)
  cat("\n\nActive categories: relative weights (in %), coordinates, contributions (in %) to axes\n")
  results <- cbind(round(100*weight.K[list_active], digits=nb_dec), 
                   round(Y_csa.KL, digits=nb_dec), 
                   round(100.*Ctr.KL,digits=2))
  print(results)
}
#---Comparison between CSA and MCA axes
#angles between axes
angles <- round((acos(abs(Cos.LL))*180/pi), digits=0)
report_axes_1 <- data.frame(apply(angles[ ,1:ncol(angles)], 2, as.integer))
rownames(report_axes_1)= rownames(Cos.LL)
# angles between principal lines
angles_espace <- matrix(0, nrow=L_ref, ncol=L_class)
angles_espace[1,] <- round( acos(abs(Cos.LL[1, ]))*180/pi, digits=0 )
namerow <- c(1:L_ref)
for (j in seq(2,L_ref)){
  angles_espace[j, ] <- round( acos( sqrt(colSums(Cos.LL[1:j, ]^2)) )*180/pi, digits=0 )
  namerow[j] <- paste(1,"to", j, sep=" ")
}
report_axes_2 <- data.frame(apply(angles_espace[ ,1:ncol(angles_espace)], 2, as.integer))
colnames(report_axes_2) <- colnames(Y_csa.KL)
rownames(report_axes_2 ) <- namerow

if(spad){
  SPAD.report.newPage (paste("Comparison between CSA and MCA axes"))
  SPAD.report.addTable(as.data.frame(Cos.LL), title = "cosinus between axes of MCA and axes of CSA")
  SPAD.report.addTable(report_axes_1, title = "Angles (in degrees) between the principal lines of CSA and of MCA" )
  SPAD.report.addTable(report_axes_2, 
                       title = "Angles (in degrees) between CSA principal lines and MCA principal subspaces")
} else{
  cat("\n\nComparison between CSA and MCA axes\n")
  cat("\n\nAngles (in degrees) between the principal lines of CSA and of MCA\n")
  print(report_axes_1, digits=0)
  cat("\n\nAngles (in degrees) between CSA principal lines and MCA principal subspaces\n")
  print(report_axes_2, digits = 0)
}
#Interpretation of axes using the  method of contributions
if (!spad){
  graphe <- paste(mygraph,".pdf", sep="")
  pdf(graphe)
}
for (i in seq(1,L_class)){
  K_Ctr_P <- which( (Ctr.KL[,i] > crit) & (Y_csa.KL[,i]>0))
  K_Ctr_M <- which( (Ctr.KL[,i] > crit) & (Y_csa.KL[,i]<0))
  Ctr_sum <- 100*(sum(Ctr.KL[K_Ctr_P,i]) + sum(Ctr.KL[K_Ctr_M,i]))
  tab1 <- cbind.data.frame(weight.K[list_active[K_Ctr_P]], Y_csa.KL[K_Ctr_P,i], 100*Ctr.KL[K_Ctr_P,i])
  colnames(tab1) <- c("relative weights", "coordinates", "contributions (in%)")
  rownames (tab1) <- rownames(Y_csa.KL)[K_Ctr_P]
  tab3 <- cbind.data.frame(weight.K[list_active[K_Ctr_M]], Y_csa.KL[K_Ctr_M,i], 100*Ctr.KL[K_Ctr_M,i])
  colnames(tab3) <-  c(" CENTRAL ZONE", " ", " ")
  rownames(tab3) <- rownames(Y_csa.KL)[K_Ctr_M]
  
  
  ####  Study of angles 
  text2 <- paste("Angles (in degrees) between CSA principal line",i,"and \n the first", min(L_ref,10), 
                 "MCA principal lines (ranked in ascending order)")
  list_axes <-order(abs(Cos.LL[ ,i]), decreasing=TRUE)
  angles <- as.data.frame(as.integer(round( acos(abs(Cos.LL[list_axes[1:min(L_ref,10)] , i]))*180/pi, digits=0)))
  colnames(angles) <- "angles"; rownames(angles) <- rownames(Cos.LL)[list_axes[1:min(L_ref,10)]]
  if(spad){
    SPAD.report.newPage (paste("Interpreting axis", i, sep=""))
    #    SPAD.report.addText (paste("  Categories whose contributions are >", )
    SPAD.report.addText (paste("Categories whose contributions are >", round(crit_value, digits=2),"%"))
    SPAD.report.addTable(tab1)
    SPAD.report.addTable(tab3, columnNames = NULL)
    SPAD.report.addText (paste("Sum of contributions: ", round(Ctr_sum, digits=0),"%"))
    SPAD.report.addText(text2)
    SPAD.report.addTable(angles)
  }
  #Comparison of CSA axes with planes of MCA
  if(L_ref > 1){
    c2 <- Cos.LL[ , i]^2
    anglePlan <- combn(min(L_ref,15), 2,  function(x) acos(sqrt(sum(c2[x])))*180/pi)
    IndPetit <- which.min(anglePlan) ; ind <- combn(1:min(L_ref,15), 2)[, IndPetit] 
    text1 <- paste(" Considering the first ", min(L_ref,15), " axes, for Axis-CS#", i,", the closest principal plane of is the plane", 
                   paste(ind[1],"-", ind[2], sep=""),
                   "\n and the (dihedral) angle is equal to ", round(anglePlan[IndPetit], digits=0), " degrees.")
  }
  if(L_ref > 2){
    angle3 <- combn(min(L_ref,15), 3,  function(x) acos(sqrt(sum(c2[x])))*180/pi)
    IndPetit <- which.min(angle3); ind <- combn(1:min(L_ref,15), 3)[, IndPetit] 
    text2 <- paste(" Considering the first ", min(L_ref,15), " axes, for Axis-CSA#", i,", the closest principal 3-dimensional subspace is the subspace", 
                   paste(ind[1],"-", ind[2],"-", ind[3], sep=""),
                   "\n and the (dihedral) angle is equal to ", round(angle3[IndPetit], digits=0), " degrees.")
    if(spad){
      SPAD.report.addText(text1)
      SPAD.report.addText(text2)
    } else{
      cat("\n\n", text1,"\n", text2)
    }
  }
  # Graph of contributions
  nbre <- length(K_Ctr_P) +  length(K_Ctr_M)
  CoordAxis <- as.vector(c(Y_csa.KL[K_Ctr_P,i], Y_csa.KL[K_Ctr_M,i]))
  Rayons <- as.vector(c(weight.K[K_Ctr_P], weight.K[K_Ctr_M]))
  names(CoordAxis) <- rownames(Y_csa.KL)[c(K_Ctr_P,K_Ctr_M)]
  Ctr.axis <- as.vector(c(Ctr.KL[K_Ctr_P,i],Ctr.KL[K_Ctr_M,i]))*100
  if (nbre > 0){
    plot(x = rep(0, nbre), y = CoordAxis,
         main = paste("Interpretation of axis", i, sep=""),
         type = "n",
         cex.main = 1,
         xlab = "Contributions (%)",
         ylab = paste("Coordinates on axis ", i),
         xlim = max(Ctr.axis) * c(-label_width, 1),  # il faut augmenter pour des étiquettes plus longues
         cex.axis =0.7,
         cex.lab = 1.15, mgp=c(1.5, 0.25, 0), tck=0.01,
         lwd = 1)
    abline(v = 0, col = "grey")
    abline(v= crit_value, col = "red")
    abline(h = 0, col = "grey", lty = 2)
    segments(x0=rep(0,nbre), x1 = Ctr.axis,
             y0 = CoordAxis, y1 = CoordAxis,
             lwd = 3,
             col = "blue")
    symbols(rep(0, nbre), CoordAxis, circles= Rayons, inches=echelle/25.4,fg="blue", 
            bg="blue", add=TRUE)
    text(x=rep(0, nbre), y = CoordAxis, labels = names(CoordAxis), pos = 2, cex = text_scale)
  }
}                            # end of the loop on i (axes of CSA)
if(!spad) { 
  sink()
  dev.off()
}
# end of result writing

## Writing of cases coordinates
if(sortInd){
  toto <- cbind.data.frame(nI.I[list_class], Yclass.IL)
  colnames(toto)[1] <- "weight"
  if(spad){
    SPAD.report.newPage("Coordinates of individuals")
    SPAD.report.addTable(toto)
  } else{
    file1 <- paste(wdir,"/", doss_name,"_coord_ind.csv", sep="")
    write.table(toto, file1, sep = sep_field, dec = dec_code, append= FALSE,
                row.names = TRUE, col.names=NA)
  }
}