#base.phen = phen.df
#base.nf = fam.nf
#Sigma_w = oualkacha_estimates$Sigmaw
#Sigma_b = oualkacha_estimates$Sigmab
#suffix = "1"

addPCLine <- function(PC, mu0_barra, out, color="black") {
  # reta do PC: y = ax+b, com a = PC_2/PC_1 e b = mu_2 - a mu_1
  slope_PC <- PC[2]/PC[1]
  intercept_PC <- mu0_barra[2] - slope_PC * mu0_barra[1]

  out <- out + geom_abline(intercept = intercept_PC, slope = slope_PC, colour=color)
  return(out)
}

geraCPsW <- function(base.phen, Sigma_w, nCPs) {
  base.phen <- as.matrix(base.phen)
  decompW <- eigen(Sigma_w)
  autoValoresW <- decompW$values
  autoVetoresW <- decompW$vectors
  autoValoresW[autoValoresW <= 0] <- 10e-4
  scoresW <- as.matrix(base.phen) %*% autoVetoresW[,1:nCPs]

  return (list(autoVetoresW=autoVetoresW, autoValoresW=autoValoresW, scoresW=scoresW))
}

geraCPsB <- function(base.phen, Sigma_b, nCPs) {
  base.phen <- as.matrix(base.phen)
  decompB <- eigen(Sigma_b)
  autoValoresB <- decompB$values
  autoVetoresB <- decompB$vectors
  scoresB <- base.phen %*% autoVetoresB[,1:nCPs]

  return (list(autoVetoresB=autoVetoresB, autoValoresB=autoValoresB, scoresB=scoresB))
}


geraCPsBWinv <- function(base.phen, Sigma_w, Sigma_b, nCPs) {
  base.phen <- as.matrix(base.phen)
  decompW <- eigen(Sigma_w)
  autoValoresW <- decompW$values
  autoVetoresW <- decompW$vectors
  autoValoresW[autoValoresW <= 0] <- 10e-4
  scoresW <- as.matrix(base.phen) %*% autoVetoresW[,1:nCPs]

  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 
  SigmaBWinv_equiv <- SigmaW_raiz_inv %*% Sigma_b %*% t(SigmaW_raiz_inv)

  decompSigmaBWinv_equiv <- eigen(SigmaBWinv_equiv)

  # Dos autovetores dessa matriz equivamente,
  # obtemos os autovetores de Sigma_b %*% Sigma_w^(-1)
  autoVetoresBWinv <- SigmaW_raiz_inv  %*% decompSigmaBWinv_equiv$vectors[,1:2]

  # Os autovalores sao os mesmos
  autoValoresBWinv <- decompSigmaBWinv_equiv$values

  # Obtendo os escores com nCPs
  scoresBWinv <- base.phen %*% autoVetoresBWinv[,1:nCPs]

  return (list(autoVetoresBWinv=autoVetoresBWinv, autoValoresBWinv=autoValoresBWinv, scoresBW=scoresBWinv))
}