calculaSigmaWSigmaB <- function(phen.data, kinshipList, base.nf) {
  base.f <- length(base.nf)
  base.n <- sum(base.nf)

  #--------------------------------------------------------------------
  # Cálculo de tau_a e tau_b por família.
  # Sao utilizados para calcular a matriz de covariancia entre famílias
  # -------------------------------------------------------------------

  tau_af <- NULL;
  tau_bf <- NULL;

  for (i in 1:base.f) {
    tau_af[[i]] <- 2*sum(diag(kinshipList[[i]]));
    tau_bf[[i]] <- 2*sum(kinshipList[[i]]);
  }


  # --------------------------------------------------------------------
  # Calculo de tau_a, tau_b e tau_c
  # --------------------------------------------------------------------

  tau_a <- 0;
  tau_b <- 0;
  tau_c <- 0;

  for (i in 1:base.f) {
    tau_a <- tau_a + tau_af[i];
    tau_b <- tau_b + tau_bf[i];
    tau_c <- tau_c + (tau_bf[i]/base.nf[i]);
  }



  # --------------------------------------------------------------------
  # Carrega a base de dados fenotípicos
  # --------------------------------------------------------------------

  # Converte para matriz
  phen.data <- as.matrix(phen.data);
  p <- ncol(phen.data);

  # Calcula o vetor medio geral da base excluindo a primeira coluna
  base.mu <- as.vector(colMeans(phen.data, na.rm=TRUE));

  # Substitui valores NA dos dados pela media
  if (any(is.na(phen.data))) {
    for (i in 1:ncol(phen.data)) 
      phen.data[,i][is.na(phen.data[,i])] <- base.mu[i];
  }

  # Calcula o vetor medio por família
  base.muf <- NULL;
  inicio <- 1;
  for (i in 1:base.f) {
    fim <- inicio + base.nf[i] - 1;
    
    if (inicio == fim) {
      dados <- phen.data[inicio,];
      base.muf[[i]] <- as.vector(dados);
    }
    else {
      dados <- phen.data[inicio:fim,];
      base.muf[[i]] <- as.vector(colMeans(dados));
    }
    
    inicio <- inicio + base.nf[i];
  }


  #----------------------------------------------------------
  # Calcula as matrizes Sb e Sw onde:
  # - Sb:  matriz da soma das variancias entre grupos (b = between)
  # - Sw:  matriz da soma das variancias dentro dos grupos (w = within)
  #------------------------------------------------------------

  # Inicializa as matrizes com zeros
  Sw <- matrix(rep(0, p^2), ncol=p);
  Sb <- matrix(rep(0, p^2), ncol=p);

  inicio <- 1;
  for (i in 1:base.f) {
    percentual <- round(i/base.f,digits=3)*100;
    print(paste("Calculando dados da familia ", i, " - ", percentual, "% completo...", sep=""));
    flush.console();
    
    fim <- inicio + base.nf[i] - 1;
    if (inicio == fim) {
      dados <- data.frame(phen.data[inicio,]);
    }
    else {
      dados <- data.frame(phen.data[inicio:fim,]);
    }
    inicio <- inicio + base.nf[i];
    
    # Media da família
    media <- base.muf[[i]];
    
    Sb <- Sb + base.nf[i] * tcrossprod(media - base.mu);
    
    for (j in 1:base.nf[i]) {
      dado <- unlist(dados[j,]);
      Sw <- Sw + tcrossprod(dado - media);
    }
  }

  #Sb
  #Sw

  #-------------------------------------------------------------------
  # Calcula as matrizes Sigma_b e Sigma_w em que:
  # Sigma_b: matriz de covariancia de Y dos efeitos entre grupos
  # Sigma_w: matriz da covariancia de Y dos efeitos aleatorios dentro do grupo
  # -------------------------------------------------------------------

  denb <- (tau_c - tau_b/base.n)/(base.f - 1) - (tau_a - tau_c)/(base.n - base.f);
  Sigma_b <- ((Sb/(base.f - 1) - Sw/(base.n - base.f))/denb);
  Sigma_b <- fixNonPositiveDefiniteMatrix(Sigma_b)
  Sigma_w <- Sw/(base.n-base.f) - Sigma_b*(tau_a - tau_c)/(base.n - base.f);
  Sigma_w <- fixNonPositiveDefiniteMatrix(Sigma_w)

  # Efeito de grupo na variancia
  efGrupo <- diag(Sigma_b)/diag(Sigma_w)
  #efGrupo

  return(list(Sigmaw=Sigma_w, Sigmab=Sigma_b, Sw=Sw, Sb=Sb))
}

geraCPs <- function(phen.data, Sigma_w, Sigma_b, suffix, base.nf) {
  phen.data <- as.matrix(phen.data)
  base.n <- dim(phen.data)[1]

  decompT <- eigen(Sigma_w + Sigma_b)
  autoValoresT <- decompT$values
  autoVetoresT <- decompT$vectors
  autoValoresT[autoValoresT <= 0] <- 1e-10
  scoresT <- phen.data %*% autoVetoresT[,1:2]


  decompW <- eigen(Sigma_w)
  autoValoresW <- decompW$values
  autoVetoresW <- decompW$vectors
  autoValoresW[autoValoresW <= 0] <- 1e-10
  scoresW <- phen.data %*% autoVetoresW[,1:2]

  colorsNames <- c("green", "blue", "red", "black", "orange", "magenta", "cyan")
  colors <- rep("darkgray", base.n)
  nColoredFam <- 5
  for (i in 1:nColoredFam) {
    colors[((i-1)*base.nf[i]+1):(i*base.nf[i])] <- colorsNames[i]
  }
  coloredids = which(colors != "darkgray")

  widths <- rep(1, base.n)
  widths[1:(nColoredFam*base.nf[i])] <- 3

  #Proporcao explicada pelos dois primeiros autovetores
  varexpl_w <- sum(autoValoresW[1:2])/sum(autoValoresW)


  #png(paste0("cps_sigma_w", suffix, ".png"), width=800, height=600)
  png(paste0("com_cps_sigma_w_b_wb", suffix, ".png"), width=1500, height=600)
  par(mfrow=c(1,3))  
  plot(scoresW[,1], scoresW[,2],
   main=paste0("Componentes Principais da Sigma_w\n",
      "Variancia explicada: ", signif(varexpl_w,4)),
   xlab="Componente Principal 1", ylab="Componente Principal 2",
   col=colors, lwd=widths, cex.lab=1.5, cex.main=1.7, cex=2)
  points(scoresW[coloredids,1], scoresW[coloredids,2], col=colors[coloredids], lwd=widths[coloredids])
  #dev.off()

  decompB <- eigen(Sigma_b)
  autoValoresB <- decompB$values
  autoVetoresB <- decompB$vectors
  scoresB <- phen.data %*% autoVetoresB[,1:2]

  #Proporcao explicada pelos dois primeiros autovetores
  varexpl_b <- sum(autoValoresB[1:2])/sum(autoValoresB)

  #png(paste0("cps_sigma_b", suffix, ".png"), width=800, height=600)
  plot(scoresB[,1], scoresB[,2],
   main=paste0("Componentes Principais da Sigma_b\n",
      "Variancia explicada: ", signif(varexpl_b,4)),
   xlab="Componente Principal 1", ylab="Componente Principal 2",
   col=colors, lwd=widths, cex.lab=1.5, cex.main=1.7, cex=2)
  points(scoresB[coloredids,1], scoresB[coloredids,2], col=colors[coloredids], lwd=widths[coloredids])
  #dev.off()


  SigmaW_raiz_inv <- autoVetoresW %*% diag(autoValoresW^(-0.5)) %*% t(autoVetoresW)

  # SigmaBW: Sigma_b %*% Sigma_w^(-1), mas usaremos uma que é garantidamente
  # simetrica, tem os mesmos autovalores e conseguimos obter os autovetores
  # a partir de uma correcao 
  SigmaBW_inv_equiv <- SigmaW_raiz_inv %*% Sigma_b %*% t(SigmaW_raiz_inv)

  decompSigmaBW_inv_equiv <- eigen(SigmaBW_inv_equiv)
  # Dos autovetores dessa matriz equivamente,
  # obtemos os autovetores de Sigma_b %*% Sigma_w^(-1)
  autoVetoresBW_inv <- SigmaW_raiz_inv  %*% decompSigmaBW_inv_equiv$vectors[,1:2]
  # Os autovalores sao os mesmo
  autoValoresBW_inv <- decompSigmaBW_inv_equiv$values
  # Obtendo os escores
  scoresBW_inv <- phen.data %*% autoVetoresBW_inv[,1:2]

  #Proporca£o explicada pelos dois primeiros autovetores
  varexpl_bw <- sum(autoValoresBW_inv[1:2])/sum(autoValoresBW_inv)

  #png(paste0("cps_sigma_bw", suffix, ".png"), width=800, height=600)
  plot(scoresBW_inv[,1], scoresBW_inv[,2],
   main=paste0("Componentes Principais da Sigma_b Sigma_w^(-1)\n",
      "Variancia explicada: ", signif(varexpl_bw,4)),
   xlab="Componente Principal 1", ylab="Componente Principal 2",
   col=colors, lwd=widths, cex.lab=1.5, cex.main=1.7, cex=2)
  points(scoresBW_inv[coloredids,1], scoresBW_inv[coloredids,2], col=colors[coloredids], lwd=widths[coloredids])
  dev.off()

  scores <- cbind(scoresBW_inv, scoresB, scoresW)

  n.phens <- dim(phen.data)[2]
  colnames(scores) <- c("CP1_BW", "CP2_BW", "CP1_B", "CP2_B", "CP1_W", "CP2_W") 

  return(list(eigvB=autoVetoresB[,1:2], eigvW=autoVetoresW[,1:2], eigvBW=autoVetoresBW_inv[,1:2], scores=scores))
}