##########################################################################################
##                                                                                      ##
##                         Ana Flávia Delbem Vidigal Nazareth                           ##
##                                                                                      ##
##                            ALGORITMO DE AGRUPAMENTOS                                 ##
##                       ***        K-PROTÓTIPOS         ***                            ##
##           Funções modificadas a partir dos pacotes ade4, clustMixType e FactoExtra   ##
##                                                                                      ##
##########################################################################################

##########################################################################################
##                       FUNÇÕES E PACOTES A SEREM UTILIZADOS                           ##
##########################################################################################

##Download dos códigos fontes utilizados:
##Modificar o diretório de destino desejado: destdir="C:/"
#download.packages("clustMixType", destdir= "C:/Users/anafl/Desktop", type="source")
#download.packages("ade4", destdir= "C:/Users/anafl/Desktop", type="source")
#download.packages("factoextra", destdir= "C:/Users/anafl/Desktop", type="source")
#download.packages("cluster", destdir= "C:/Users/anafl/Desktop", type="source")
#download.packages("NbClust", destdir= "C:/Users/anafl/Desktop", type="source")

# Criação da função para o cálculo da distância euclidiana
cppFunction('NumericMatrix weighted_distance (NumericMatrix x, NumericMatrix y, NumericVector lambda){
            
            int n_x = x.nrow();
            int n_y = y.nrow();
            
            
            NumericMatrix out(n_x, n_y);
            
            //begin the loop
            
            for (int i = 0 ; i < n_x; i++){
            for (int j = 0  ; j < n_y ; j ++) {
            double d = sum(pow(x.row(i) - y.row(j), 2)*lambda);
            out(i,j) = d;
            }
            }
            return (out) ;
            }')  
  

##########################################################################################
##                   FUNÇÃO Kproto() MODIFICADA: kproto.modif()                         ##
##########################################################################################

##TESTE
#x <- dados1 #Data frame com ambas variáveis quantitativas e fatoriais.
#k = K  #Pode ser um número de grupos, um vetor especificando índices de protótipos iniciais ou um data frame de protótipos com as mesmas colunas que x.
#lambda = LAMB #Parâmetro > 0 para ponderar entre a distância euclidiana de variáveis quantitativas e a distância escolhida entre variáveis categóricas. Também é possível ser um vetor de pesos específicos para as variáveis onde o primeiro termo deve ser o peso dado às variáveis qualitativas e os demais termos correspondem à ordem das variáveis quantitativas nos dados. Nesse caso, cada distância entre cada variável quantitativa será multiplicada pelo valor lambda correspondente àquela variável e a distância entre as variáveis qualitativas será multiplicada pelo lambda das variáveis qualitativas único. Deste modo, deve ser um vetor de dimensão igual ao número de variáveis quantitativas mais um.
#iter.max = 100 #Número máximo de iterações, se não houver convergência antes.
#nstart = 1 #Se > 1 cálculos repetitivos com inicializações aleatórias são calculados e o resultado com um mínimo de tot.dist é retornado.
#na.rm = TRUE #Um valor lógico que indica se os valores ausentes (NAs) devem ser removidos antes que a computação prossiga.
#keep.data = FALSE #Um valor lógico que indica se o banco de dados original deve ser incluído no objeto retornado.
#verbose = TRUE #Um valor lógico que indica se as informações sobre o procedimento de agrupamentos devem ser dadas. Cuidado: Se verbose = FALSE, a redução do número de grupos não é mencionada.

kproto.modif <- function(
           x,
           k,
           lambda = NULL,
           iter.max = 100,
           nstart = 1,
           na.rm = TRUE,
           keep.data = TRUE,
           verbose = TRUE,
           ...) {
    
    # Checagem de erros inicial
    if (!is.data.frame(x))
      stop("x deve ser um data frame!")
    if (ncol(x) < 2)
      stop("Para agrupamentos x deve conter pelo menos duas variaveis!")
    if (iter.max < 1 |
        nstart < 1)
      stop("iter.max e nstart nao devem ser especificados como < 1!")
    if (!is.null(lambda)) {
      if (any(lambda < 0))
        stop("lambda deve ser especificado como >= 0!")
      if (!any(lambda > 0))
        stop("lambda deve ser especificado como > 0 para pelo menos uma variavel!")}
    
    # Verificação para variáveis quantitativas e fatoriais
    numvars <- sapply(x, is.numeric)
    anynum <- any(numvars)
    catvars <- sapply(x, is.factor)
    anyfact <- any(catvars)
    if (!anynum)
      stop("\n Nao ha variaveis quantitativas em x! Tente usar o metodo kmodes() do pacote klaR...\n\n")
    if (!anyfact)
      stop("\n Nao ha variaveis fatoriais em x! Tente usar o metodo kmeans()...\n\n")
    
    # Tratamento de valores ausentes
    NAcount <- apply(x, 2, function(z)
      sum(is.na(z)))
    if (verbose) {
      cat("# NAs nas variaveis:\n")
      print(NAcount)
      }
    if (any(NAcount == nrow(x)))
      stop(paste("Variavel(is) possuem apenas NAs, favor remove-las:",
                 names(NAcount)[NAcount == nrow(x)], "!"))
    if (na.rm) {
      miss <- apply(x, 1, function(z)
        any(is.na(z)))
      if (verbose) {
        cat(sum(miss), "Observacao(es) com NAs.\n")
        if (sum(miss) > 0)
          message("Observacoes com NAs sao removidas.\n")
        cat("\n")}
      x <- x[!miss,]} # Remoção de observações com valores ausentes
    
    if (!na.rm) {
      allNAs <- apply(x, 1, function(z)
        all(is.na(z)))
      if (sum(allNAs) > 0) {
        if (verbose)
          cat(sum(allNAs), "Observacao(es) onde todas as variaveis nao estao presentes.\n")
        warning("Nenhuma atribuicao de grupos significativa e possivel para observacoes onde todas as variaveis nao estao presentes.\n")
        if (verbose)
          cat("\n")}}
    
    if (nrow(x) == 1)
      stop("Apenas o agrupamento de uma observacao nao e significativo.")
    
    k_input <- k # Armazena o valor de k para nstart > 1 uma vez que os grupos podem ser fundidos
    
    # Inicialização dos protótipos
    if (!is.data.frame(k)) {
      if (length(k) == 1) {
        if (as.integer(k) != k) {
          k <- as.integer(k)
          warning(paste("k foi modificado para", k, "!"))}
        if (nrow(x) < k)
          stop("O data frame tem um numero de observacoes menor que o de grupos!")
        ids <- sample(nrow(x), k)
        protos <- x[ids,]}
      if (length(k) > 1) {
        if (nrow(x) < length(k))
          stop("O data frame tem um numero de observacoes menor que o de grupos!")
        ids <- k
        k <- length(ids)
        if (length(unique(ids)) != length(ids))
          stop("Se k e especificado como um vetor ele deve conter diferentes indices!")
        if (any(ids < 1) |
            any(ids > nrow(x)))
          stop("Se k e especificado como um vetor todos os elementos devem ser indices validos de x!")
        #Verificação de inteiros
        protos <- x[ids,]}
      rm(ids)}
    if (is.data.frame(k)) {
      if (nrow(x) < nrow(k))
        stop("O data frame tem um numero de observacoes menor que o de grupos!")
      if (length(names(k)) != length(names(x)))
        stop("k e x possuem um numero diferente de colunas!")
      if (any(names(k) != names(x)))
        stop("k e x possuem nomes de colunas diferentes!")
      if (anynum) {
        if (any(sapply(k, is.numeric) != numvars))
          stop("Variaveis quantitativas de k e x nao sao correspondentes!")}
      if (anyfact) {
        if (any(sapply(k, is.factor) != catvars))
          stop("Variaveis fatoriais de k e x nao sao correspondentes!")}
      protos <- k
      k <- nrow(protos)}
    if (k < 1)
      stop("O numero de grupos k nao deve ser menor que 1!")
    
    # Cálculo automático de lambda
    if (length(lambda) > 1) {
      if (length(lambda) != (sum(c(numvars))+1)) #!= (sum(c(numvars, catvars)))
        stop("se lambda e um vetor, seu tamanho deve ser de um mais a soma do numero de variaveis quantitativas do data frame!")}
    if (is.null(lambda)) {
      warning("Lambda precisa ser informado, como escalar ou vetor!")
    if (anynum & anyfact) {
       vnum <- mean(sapply(x[, numvars, drop = FALSE], var, na.rm = TRUE))
       vcat <- mean(sapply(x[, catvars, drop = FALSE], function(z) return(1 - sum((table(z) / sum(!is.na(z))) ^ 2))))
       if (vnum == 0) {
         if (verbose)
           warning("Todas as variaveis quantitativas possuem variancia zero.")
         anynum <- FALSE}
       if (vcat == 0) {
         if (verbose)
           warning("Todas as variaveis categoricas possuem variancia zero.")
         anyfact <- FALSE}
       if (anynum & anyfact) {
         lambda <- vnum / vcat
         if (verbose)
           cat("Lambda estimado:", lambda, "\n\n")}
         else{
         lambda <- 1}}}
    
    # Inicialização cálculos dos grupos
    clusters  <- numeric(nrow(x))
    tot.dists <- NULL
    moved   <- NULL
    iter <- 1
    
    # Verificação de protótipos iguais e redução do número de grupos em caso de ocorrência
    if (k > 1) {
      keep.protos <- rep(TRUE, k)
      for (l in 1:(k - 1)) {
        for (m in (l + 1):k) {
          d1 <- sum((protos[l, numvars, drop = FALSE] - protos[m, numvars, drop = FALSE]) ^ 2) # Euclidiana para variáveis quantitativas
          d2 <- sum(protos[l, catvars, drop = FALSE] != protos[m, catvars, drop = FALSE]) # Coeficiente de concordância simples ponderado para variáveis categóricas
          if ((d1 + d2) == 0)
            keep.protos[m] <- FALSE}}
      if (!all(keep.protos)) {
        protos <- protos[keep.protos,]
        k <- sum(keep.protos)
        if (verbose)
          message("Prototipos iguais fundidos. Numero de grupos reduzido para:", k, "\n\n")}}
    
    # Caso especial de apenas um grupo
    if (k == 1) {
      clusters <- rep(1, nrow(x))
      size  <- table(clusters)
      iter <- iter.max} # REM: o tamanho do vetor nomeado é necessário mais tarde...
    
    # Transformação e desmembramento do data frame de dados para os cálculos necessários quando K > 1
    x.num <- x[,numvars]
    x.num <- as.matrix(x.num)
    
    x.cat <- x[,catvars]
    x.cat <- data.frame(lapply(x.cat, as.character), stringsAsFactors=FALSE)
    x.cat <- data.frame(lapply(x.cat, as.numeric), stringsAsFactors=FALSE)
    x.cat <- as.matrix(1 * (x.cat > 0))
    
    if (length(lambda) > 1) {
      lambda.num <- lambda[-1]
      as.numeric(unlist(lambda.num))
      lambda.num <- as.matrix(lambda.num)
      lambda.num <- t(lambda.num)
      #lambda.num
      #mode(lambda.num)
      
      lambda.cat <- lambda[1]} #Lambda para variáveis categóricas é único para todas. Corresponde ao primeiro elemento do vetor.
    
    # Início das iterações para o caso padrão (que é k > 1)
    while (iter < iter.max){
      #a0 <- proc.time()[3]
      
      # Transformação e desmembramento do data frame dos protpótipos para os cálculos necessários
      protos.num <- protos[,numvars]
      protos.num <- as.matrix(protos.num)
      
      protos.cat <- protos[,catvars]
      protos.cat <- data.frame(lapply(protos.cat, as.character), stringsAsFactors=FALSE)
      protos.cat <- data.frame(lapply(protos.cat, as.numeric), stringsAsFactors=FALSE)
      protos.cat <- as.matrix(1 * (protos.cat > 0))
      
      # Cálculo das distâncias
      nrows <- nrow(x)
      dists <- matrix(NA, nrow = nrows, ncol = k)
      v.num <- ncol(x.num) # = ncol(protos.num) #Número total de variáveis quantitativas
      #v.cat <- length(catvars[catvars==TRUE])
      v.cat <- ncol(x.cat) #Número total de variáveis categóricas
      
      # Distância Euclidiana para variáveis quantitativas
      d1 <- matrix(NA, nrow = nrows, ncol = k)
      if (length(lambda) == 1) { 
        d1 <- outer(rowSums(x.num^2), rowSums(protos.num^2), '+') - tcrossprod(x.num, 2 * protos.num)}
      
      if (length(lambda) > 1){
        d1 <- weighted_distance(x.num, protos.num, lambda = lambda)}
      
      #a1 <- proc.time()[3]
      
      # Distância para variáveis binárias
      a <- x.cat %*% t(protos.cat)
      b <- x.cat %*% (1 - t(protos.cat))
      c <- (1 - x.cat) %*% t(protos.cat)
      d <- ncol(x.cat) - a - b - c
      
      ## Método de cálculo das distâncias para as variáveis binárias (nominais)
      #Escolher o índice desejado
      #Alterar também na "Atualização final dos protótipos e das distâncias"
      d2 <- a / (a + b + c) #JACCARD index (1901)
      #d2 <- (a + d) %/% v.cat #Coeficiente de Concordância Simples
      #d2 <- (a + d) %/% (a + b + c + d) #SOCKAL & MICHENER index (1958)
      #d2 <- a %/% (a + 2 %*% (b + c)) #SOCKAL & SNEATH(1963)
      #d2 <- (a + d) %/% (a + 2 %*% (b + c) + d) #ROGERS & TANIMOTO (1960)
      #d2 <- 2 %*% a %/% (2 %*% a + b + c) #CZEKANOWSKI (1913) or SORENSEN (1948)
      #d2 <- (a - (b + c) + d) %/% (a + b + c + d) #index of GOWER & LEGENDRE (1986)
      #d2 <- a %/% sqrt((a+b) %*% (a+c)) #OCHIAI (1957)
      #d2 <- a %*% d %/% sqrt((a + b) %*% (a + c) %*% (d + b) %*% (d + c)) #SOKAL & SNEATH (1963)
      #d2 <- (a %*% d - b %*% c) %/% sqrt((a + b) %*% (a + c) %*% (b + d) %*% (d + c)) #Phi of PEARSON
      #d2 <- a %/% (a + b + c + d) #coefficient of GOWER & LEGENDRE
      #      diag(d2) <- 1
      #d2 <- Inserir cálculo caso queira utilizar outro coeficiente de similaridade.
      #a3 <- proc.time()[3]
      
      #Cálculo da dissimilaridade (distância), a partir do coeficiente de similaridade
      d2 <- 1 - d2 #Apenas usar esta relação se d2 for uma medida de distância
      #d2 <- (b+c) %/% (a + b + c) #Distância de Jaccard.
      #d2 <- (v.cat-(a+d)) %/% v.cat #Distância a partir do Coeficiente de Concordância Simples
      #d2 <- Modificar caso queira utilizar uma relação entre d e s diferente.
      #d2 = Matriz de distâncias entre os n objetos do banco de dados e os k protótipos [n,k]
      
      d2[is.na(d2)] <- FALSE
      if (length(lambda) == 1){
        d2 <- lambda * d2}
      if (length(lambda) > 1){
        d2 <- lambda.cat * d2}
      
      #Matriz de distâncias combinadas entre os objetos e os k protótipos
      dists <- d1 + d2
      
      #cat(a1-a0, a2-a1, "\n")
      
      # Determinação dos grupos
      old.clusters <- clusters
      # clusters <- apply(dists, 1, function(z) which.min(z))
      clusters <- apply(dists, 1, function(z) {
        a <- which.min(z)
        if (length(a) > 1)
          a <- sample(a, 1)
        return(a)}) # Amostragem aleatória no caso de múltiplos mínimos
      size  <- table(clusters)
      min.dists <- apply(cbind(clusters, dists), 1, function(z) z[z[1] + 1])
      within <- as.numeric(by(min.dists, clusters, sum))
      tot.within <- sum(within)
      # Prevenção de classes vazias
      #tot.within    <- numeric(k)
      #totw.list     <- by(min.dists, clusters, sum)
      #tot.within[names(totw.list)] <- as.numeric(totw.list)
      
      # ...Verificação de grupos vazios e eventual redução do número de protótipos
      if (length(size) < k) {
        k <- length(size)
        protos <- protos[1:length(size),]
        if (verbose)
          cat("Ocorrencia de grupos vazios. Numero de grupos reduzido para:", k, "\n\n")}
      
      # Rastreamento
      tot.dists <- c(tot.dists, sum(tot.within))
      moved <- c(moved, sum(clusters != old.clusters))
      
      # Cálculo dos novos protótipos
      remids <- as.integer(names(size))
      for (i in remids) {
        protos[which(remids == i), numvars] <-
          sapply(x[clusters == i, numvars, drop = FALSE], mean, na.rm = TRUE)
        protos[which(remids == i), catvars] <-
          sapply(x[clusters == i, catvars, drop = FALSE], function(z)
            levels(z)[which.max(table(z))])}
      
      if (k == 1) {
        clusters <- rep(1, length(clusters))
        size <- table(clusters)
        iter <- iter.max
        break}
      
      # Verificação de protótipos iguais e eventual redução do número de protótipos
      if (iter == (iter.max - 1)) {
        # REM: para a última iteração protótipos iguais são permitidos, caso contrário, menos protótipos que grupos designados.
        keep.protos <- rep(TRUE, k)
        for (l in 1:(k - 1)) {
          for (m in (l + 1):k) {
            d1 <- sum((protos[l, numvars, drop = FALSE] - protos[m, numvars, drop = FALSE]) ^ 2) # euclidean for numerics
            d2 <- sum(protos[l, catvars, drop = FALSE] != protos[m, catvars, drop = FALSE]) # wtd simple matching for categorics
            if ((d1 + d2) == 0)
              keep.protos[m] <- FALSE}}
        if (!all(keep.protos)) {
          protos <- protos[keep.protos,]
          k <- sum(keep.protos)
          if (verbose)
            cat(" Prototipos iguais fundidos.Numero de grupos reduzido para:", k, "\n\n")}}
      
      # Adição de regras de parada
      if (moved[length(moved)] ==  0)
        break
      
      if (k == 1) {
        clusters <- rep(1, length(clusters))
        size <- table(clusters)
        iter <- iter.max
        break}
      
      #cat("iter", iter, "moved", moved[length(moved)], "tot.dists",tot.dists[length(tot.dists)],"\n" )
      iter <- iter + 1}
    
    
    ### Atualização final dos protótipos e das distâncias
    if (iter == iter.max) {
      # Caso contrário, não haverá mais movimentos e os protótipos finais não serão atualizados
      # Cálculo dos novos protótipos
      remids <- as.integer(names(size))
      for (i in remids) {
        protos[which(remids == i), numvars] <-
          sapply(x[clusters == i, numvars, drop = FALSE], mean, na.rm = TRUE)
        protos[which(remids == i), catvars] <-
          sapply(x[clusters == i, catvars, drop = FALSE], function(z)
            levels(z)[which.max(table(z))])}
      
      # Cálculo das distâncias
      
      protos.num <- protos[,numvars]
      protos.num <- as.matrix(protos.num)
      
      protos.cat <- protos[,catvars]
      protos.cat <- data.frame(lapply(protos.cat, as.character), stringsAsFactors=FALSE)
      protos.cat <- data.frame(lapply(protos.cat, as.numeric), stringsAsFactors=FALSE)
      protos.cat <- as.matrix(1 * (protos.cat > 0))
      
      nrows <- nrow(x)
      dists <- matrix(NA, nrow = nrows, ncol = k)
      
      # Distância Euclidiana para variáveis quantitativas
      if (length(lambda) == 1) { 
        d1 <- outer(rowSums(x.num^2), rowSums(protos.num^2), '+') - tcrossprod(x.num, 2 * protos.num)}
      
      if (length(lambda) > 1){
        d1 <- weighted_distance(x.num, protos.num, lambda = lambda)}
      #a1 <- proc.time()[3]
      
      # Distância para variáveis binárias
      a <- x.cat %*% t(protos.cat)
      b <- x.cat %*% (1 - t(protos.cat))
      c <- (1 - x.cat) %*% t(protos.cat)
      d <- ncol(x.cat) - a - b - c
      
      ## Método de cálculo das distâncias para as variáveis binárias (nominais)
      d2 <- a / (a + b + c) #JACCARD index (1901)
      #d2 <- (a + d) %/% v.cat #Coeficiente de Concordância Simples
      #d2 <- (a + d) %/% (a + b + c + d) #SOCKAL & MICHENER index (1958)
      #d2 <- a %/% (a + 2 %*% (b + c)) #SOCKAL & SNEATH(1963)
      #d2 <- (a + d) %/% (a + 2 %*% (b + c) + d) #ROGERS & TANIMOTO (1960)
      #d2 <- 2 %*% a %/% (2 %*% a + b + c) #CZEKANOWSKI (1913) or SORENSEN (1948)
      #d2 <- (a - (b + c) + d) %/% (a + b + c + d) #index of GOWER & LEGENDRE (1986)
      #d2 <- a %/% sqrt((a+b) %*% (a+c)) #OCHIAI (1957)
      #d2 <- a %*% d %/% sqrt((a + b) %*% (a + c) %*% (d + b) %*% (d + c)) #SOKAL & SNEATH (1963)
      #d2 <- (a %*% d - b %*% c) %/% sqrt((a + b) %*% (a + c) %*% (b + d) %*% (d + c)) #Phi of PEARSON
      #d2 <- a %/% (a + b + c + d) #coefficient of GOWER & LEGENDRE
      #      diag(d2) <- 1
      #d2 <- Inserir cálculo caso queira utilizar outro coeficiente de similaridade.
      #a3 <- proc.time()[3]
      
      #Cálculo da dissimilaridade (distância), a partir do coeficiente de similaridade
      d2 <- 1 - d2 #Apenas usar esta relação se d2 for uma medida de distância
     #d2 <- (b+c) %/% (a + b + c) #Distância de Jaccard.
     #d2 <- (v.cat-(a+d)) %/% v.cat #Distância a partir do Coeficiente de Concordância Simples
     #d2 <- Modificar caso queira utilizar uma relação entre d e s diferente.
      #d2 = Matriz de distâncias entre os n objetos do banco de dados e os k protótipos [n,k]
    
      d2[is.na(d2)] <- FALSE
      if (length(lambda) == 1){
        d2 <- lambda * d2}
      if (length(lambda) > 1){
        d2 <- lambda.cat * d2}
    
      #Matriz de distâncias combinadas entre os objetos e os k protótipos
     dists <- d1 + d2}
    
    size          <- table(clusters)
    min.dists     <- apply(cbind(clusters, dists), 1, function(z) z[z[1] + 1])
    within        <- as.numeric(by(min.dists, clusters, sum))
    tot.within    <- sum(within)
    
    names(clusters) <- row.names(dists) <- row.names(x)
    rownames(protos) <- NULL
    
    # Criação do resultado:
    res <- list(
      data = dados1,
      data.numeric = dados.num.n, # Modificado para a função fviz_cluster.modif
      data.categorical = dados.cat,
      cluster.n = clusters,
      cluster = as.factor(clusters), # Modificado para a função fviz_cluster.modif
      centers = protos,
      #K = K, #Não colocar para construir o screeplot
      lambda = lambda,
      size = size,
      withinss = within,
      tot.withinss = tot.within,
      dists = dists,
      d1 = d1,
      d2 = d2,
      #lambda.cat = lambda.cat,
      #lambda.num = lambda.num,
      #Inserir variável desejada para ser armazenada no resultado da função
      iter = iter,
      trace = list(tot.dists = tot.dists, moved = moved))
    
    # loop: se nstart > 1:
    if (nstart > 1)
      for (j in 2:nstart) {
        res.new <- kproto.modif(
          x = x,
          k = k_input,
          lambda = lambda,
          iter.max = iter.max,
          nstart = 1,
          verbose = verbose)
        if (res.new$tot.withinss < res$tot.withinss)
          res <- res.new}
    
    if (keep.data)
      res$data = x
    class(res) <- "kproto.modif"
    return(res)} ####Fim da função kproto.modif
  

##########################################################################################
##             FUNÇÃO fviz_cluster() MODIFICADA: fviz_cluster.modif()                   ##
##########################################################################################

fviz_cluster.modif <- function(object, data = NULL, choose.vars = NULL, stand = TRUE, 
                               axes = c(1, 2),
                               geom = c("point", "text"), repel = FALSE,
                               show.clust.cent = TRUE,
                               ellipse = TRUE, ellipse.type = "convex", ellipse.level = 0.95,
                               ellipse.alpha = 0.2,
                               shape = NULL, pointsize = 1.5, labelsize = 12,
                               main = "Cluster plot",  xlab = NULL, ylab = NULL,
                               outlier.color = "black", outlier.shape = 19,
                               ggtheme = theme_grey(), ...){
  
  # Check axis lengths
  .check_axes <- function(axes, .length){
    if(length(axes) != .length) stop("axes should be of length ", 2)}
  
  # Deprecated arguments
  extra_args <- list(...)
  .check_axes(axes, .length = 2)
  
  if (!is.null(extra_args$jitter)) {
    warning("argument jitter is deprecated; please use repel = TRUE instead, to avoid overlapping of labels.", 
            call. = FALSE)
    if(!is.null(extra_args$jitter$width) | !is.null(extra_args$jitter$height) ) repel = TRUE}
  
  if(!is.null(extra_args$frame)) ellipse <- .facto_dep("frame", "ellipse", ellipse)
  if(!is.null(extra_args$frame.type)) ellipse.type <- .facto_dep("frame.type", "ellipse.type", extra_args$frame.type)
  if(!is.null(extra_args$frame.level)) ellipse.level <- .facto_dep("frame.level", "ellipse.level", extra_args$frame.level)
  if(!is.null(extra_args$frame.alpha)) ellipse.alpha <- .facto_dep("frame.alpha", "ellipse.alpha", extra_args$frame.alpha)
  if(!is.null(extra_args$title)) main <- .facto_dep("title", "main", extra_args$title)
  
  # object from cluster package
  if(inherits(object, c("partition", "hkmeans", "eclust"))) data <- object$data
  
  # Object from kmeans (stats package)
  else if((inherits(object, "kmeans") & !inherits(object, "eclust"))| inherits(object, "dbscan")){
    if(is.null(data)) stop("data is required for plotting kmeans/dbscan clusters")} 
  # Object from mclust package
  else if(inherits(object, "Mclust")) {
    object$cluster <- object$classification
    data <- object$data}
  # HCPC in FactoMineR
  else if(inherits(object, "HCPC")) {
    object$cluster <- object$call$X$clust
    data <- res.hcpc <- object
    stand <- FALSE # to avoid trying to standardize HCPC results
    #     data <- object$data.clust[, -ncol(object$data.clust), drop = FALSE]
    #     object$cluster <- as.vector(object$data.clust$clust)
    }
  else if(inherits(object, "hcut")){
    if(inherits(object$data, "dist")){
      if(is.null(data)) stop("The option 'data' is required for an object of class hcut." )}
    else data <- object$data}
  
  
  # Modificação para se adequar à função kproto.modif
  # Object from kproto.modif function
  else if(inherits(object, "kproto.modif")) {
    object$cluster <- object$cluster
    data <- object$data.numeric}
  
  
  # Any obects containing data and cluster elements
  else if(!is.null(object$data) & !is.null(object$cluster)){
    data <- object$data
    cluster <- object$cluster}
  else stop("Can't handle an object of class ", class(object))
  
  # Choose variables
  if(!is.null(choose.vars))
    data <- data[, choose.vars, drop = FALSE]
  if(stand) data <- scale(data)
  cluster <- as.factor(object$cluster)
  
  pca_performed <- FALSE
  
  # Prepare the data for plotting
  # ++++++++++++++++++++++++
  # PCA is performed depending on the number of variables
  if(inherits(data, c("matrix", "data.frame"))){
    # ncol(data) > 2 --> PCA
    if(ncol(data)>2){
      pca <- stats::prcomp(data, scale = FALSE, center = FALSE)
      ind <- facto_summarize(pca, element = "ind", result = "coord", axes = axes)
      eig <- get_eigenvalue(pca)[axes,2]
      if(is.null(xlab)) xlab = paste0("Dim", axes[1], " (", round(eig[1],1), "%)") 
      if(is.null(ylab)) ylab = paste0("Dim", axes[2], " (", round(eig[2], 1),"%)")}
    # PCA is not performed
    else if(ncol(data) == 2){
      ind <- as.data.frame(data)
      ind <- cbind.data.frame(name = rownames(ind), ind)
      if(is.null(xlab)) xlab <- colnames(data)[1]
      if(is.null(ylab)) ylab <- colnames(data)[2]
      
      if(xlab=="x") xlab <- "x value"
      if(ylab == "y") ylab <- "y value"}
    else{
      stop("The dimension of the data < 2! No plot.")}
    colnames(ind)[2:3] <-  c("x", "y")
    label_coord <- ind}
  else if(inherits(data, "HCPC")){
    ind <- res.hcpc$call$X[, c(axes, ncol(res.hcpc$call$X))]
    colnames(ind) <- c("Dim.1", "Dim.2", "clust")
    ind <- cbind.data.frame(name = rownames(ind), ind)
    colnames(ind)[2:3] <-  c("x", "y")
    label_coord <- ind
    eig <- get_eigenvalue(res.hcpc$call$t$res)[axes,2]
    if(is.null(xlab)) xlab = paste0("Dim", axes[1], " (", round(eig[1], 1), "%)") 
    if(is.null(ylab)) ylab = paste0("Dim", axes[2], " (", round(eig[2], 1),"%)")}
  else stop("A data of class ", class(data), " is not supported.")
  
  # Plot data and labels
  # ++++++++++++++++++++++++
  label = FALSE
  if("text" %in% geom) label = TRUE
  if(!("point" %in% geom)) pointsize = 0
  
  plot.data <- cbind.data.frame(ind, cluster = cluster)
  label_coord <- cbind.data.frame(label_coord, cluster = cluster)
  # Augment data
  if(inherits(object, "Mclust")){
    plot.data$uncertainty <- object$uncertainty
    label_coord$uncertainty <- object$uncertainty}
  
  # IF DBSCAN: cluster 0 is outliers. We don't want to make ellipse around
  # these observations. Let's remove them. They will be added to the plot later
  is_outliers = FALSE
  if(inherits(object, c("dbscan", "Mclust"))){
    outliers <- which(cluster == 0)
    if(length(outliers) > 0){
      is_outliers = TRUE
      outliers_data <- plot.data[outliers, , drop = FALSE]
      outliers_labs <- label_coord[outliers, , drop = FALSE]
      
      ind <- ind[-outliers, , drop = FALSE]
      cluster <- cluster[-outliers]
      plot.data <- plot.data[-outliers, , drop = FALSE]
      label_coord <- label_coord[-outliers, , drop = FALSE]}}
  
  # Plot
  # ++++++++++++++++++++++++
  lab <- NULL
  if("text" %in% geom) lab <- "name"
  if(is.null(shape)) shape <- "cluster"
  
  if(inherits(object, "partition") & missing(show.clust.cent))
    show.clust.cent <- FALSE # hide mean point for PAM, CLARA
  
  p <- ggpubr::ggscatter(plot.data, "x", "y",
                         color="cluster", shape = shape, size = pointsize,
                         point = "point" %in% geom, 
                         label = lab,
                         font.label = labelsize, repel = repel,
                         mean.point = show.clust.cent, 
                         ellipse = ellipse, ellipse.type = ellipse.type,
                         ellipse.alpha = ellipse.alpha, ellipse.level = ellipse.level,
                         main = main, xlab = xlab, ylab = ylab,
                         ggtheme = ggtheme, ...)
  
  # Add outliers (can exist only in dbscan)
  if(is_outliers)
    p <- .add_outliers(p, outliers_data, outliers_labs, outlier.color, outlier.shape,
                       pointsize, labelsize, geom, repel = repel)
  
  p}

# Add outliers to cluster plot (for dbscan only)
.add_outliers <-function(p, outliers_data, outliers_labs, 
                         outlier.color = "black", outlier.shape = 19,
                         pointsize = 2, labelsize = 4, geom = c("point", "text"), repel = FALSE){
  if("point" %in% geom) 
    p <-  p + geom_point(data = outliers_data, 
                         aes_string('x', 'y'),
                         size = pointsize, color = outlier.color, shape = outlier.shape)
  if("text" %in% geom){
    if(repel)
      p <- p +ggrepel::geom_text_repel(data = outliers_labs, 
                                       aes_string('x', 'y', label = 'name'),
                                       size = labelsize, color = outlier.color)
    else
      p <- p + geom_text(data = outliers_labs, 
                         aes_string('x', 'y', label = 'name'),  
                         size = labelsize, vjust = -0.7, color = outlier.color)}
  
  return(p)}


##########################################################################################
##              FUNÇÃO clprofiles() MODIFICADA: clprofiles.modif()                      ##
##########################################################################################

## Função para visualização dos resultados do agrupamento modificada a partir da função "clprofiles" do pacote "clustMixType"
clprofiles.modif <- function(object, x, vars = NULL, col = NULL){
  if(length(object$cluster.n) != nrow(x)) stop("Size of x does not match cluster result!")
  if(is.null(vars)) vars <- 1:ncol(x)
  if(!is.numeric(vars)) vars <- sapply(vars, function(z) return(which(colnames(x)==z)))
  if(length(vars) < 1) stop("Specified variable names do not match x!")
  if(is.null(col)){
    k <- max(unique(object$cluster.n))
    if(k > 2)  col <- brewer.pal(k, "Set3")
    if(k == 2) col <- c("lightblue","orange")
    if(k == 1) col <- "lightblue"}
  #clusids <- as.numeric(names(object$size)) # object size not named for kmeans
  clusids <- sort(unique(object$cluster.n)) 
  if(length(col) != max(clusids)) warning("Length of col should match number of clusters!")
  
  par(ask=TRUE)
  for(i in vars){
    if(is.numeric(x[,i])){
      boxplot(x[,i]~object$cluster.n, col = col, main = colnames(x)[i])
      legend("topright", legend=clusids, fill = col)} 
    if(is.factor(x[,i])){
      tab <- table(x[,i], object$cluster.n)
      for(j in 1:length(object$size)) tab[,j] <- tab[,j]/object$size[j]
      barplot(t(tab), beside = TRUE, main = colnames(x)[i], col = col, ylim=c(0,1))}} #Adição dos limites do eixo Y.
  par(ask="FALSE")
  invisible()}


##########################################################################################
##           FUNÇÃO summary.kproto() MODIFICADA: summary.kproto.modif()                 ##
##########################################################################################

summary.kproto.modif <- function(object, data = NULL, pct.dig = 3, ...){
  if(class(object) != "kproto.modif") stop("object must be of class kproto.modif!")
  
  if(is.null(data)){
    data <- object$data
    cluster <- as.numeric(object$cluster)}
  #if(!is.null(data)) cluster <- predict(object, data)$cluster
  
  numvars <- sapply(data, is.numeric)
  #anynum <- any(numvars)
  catvars <- sapply(data, is.factor)
  #anyfact <- any(catvars)
  
  res <- NULL
  for(i in 1:ncol(data)){
    cat(names(data)[i],"\n")
    if(numvars[i]){
      resi <- by(data[,i], cluster, summary, ...)
      res[[i]] <- matrix(unlist(resi), nrow = length(unique(cluster)), byrow=TRUE)
      colnames(res[[i]]) <- names(resi[[1]])
      rownames(res[[i]]) <- sort(unique(cluster))
    }
    if(catvars[i])  res[[i]] <- round(prop.table(table(cluster, data[,i]),1), digits = pct.dig)
    print(res[[i]])
    cat("\n-----------------------------------------------------------------\n")}
  names(res) <- names(data)  
  
  #return(res)
  invisible(res)}