library(gtools)
dyn.load("segregacao.so")

p <- .2 # frequencia alélica média
g <- 50 # número de genes
prole <- 10 # número de filhos por cruzamento
n <- 450 # tamanho da população
h2 <- .80 # herdabilidade do caráter
gmd <- 0 # grau médio de dominância
int <- 10 # número de indivíduos selecionados

pop.inicial<-function(p, g, n, h2, gmd){ # argumentos da função
  id <- matrix(0, ncol = 3, nrow = n) # cria matriz de identificação dos indivíduos da população
  gen <- matrix(0, ncol = g, nrow = n) # cria a matriz que armazena dos genótipos
  for(i in 1:n){ # 'laço' - atribui valores p/ cada loco p/ todos indivíduos
    f <- rbeta(g,  p * 10, 10 * (1 - p)) # dist. beta - das frequencias aléilicas p/ loco
    n01 <- runif(g, 0, 1) # dist unif - gerará os velores genotípicos por loco
    const <- numeric() # vetor qualquer - genótipos 'loco a loco' por indivíduos
    for(j in 1:g){ # 'laço' - aplica a regra para AA Aa ou aa
      # se unif < p^2, loco recebe 'r' (recessivo)
      # se unif entre p^2 e 2p(1-p) loco recebe 'h' (heterozigoto)
      # se unif > que (1-p)^2 loco recebe 'd' (dominante)
      if(n01[j] < f[j]^2) const[j] <- 'd' else
      if(n01[j] < (f[j]^2 + 2 * f[j] * (1 - f[j]))) const[j] <- 'h' else
      const[j] <- 'r'
    }
    gen[i,] <- const # armazena os valores genotípicos do indivíduo 'i' da população
  }
  vg <- apply(gen, 1, soma.vg, gmd = gmd) # calcula os valores genotípicos totais dos indivíduos
  varg <- var(vg) # calcula a variância genética
  fen <- vg + rnorm(n, 0, sqrt(varg * (1 - h2) / h2)) # atribui o desvio fenotípico
  resp <- data.frame(id, gen, vg, fen) # planilha com os individuos da população
  names(resp) <- c('fami','pai','mae',1:g,'vg','fen') # dá os nomes as colunas da planilha
  return(resp) # retorna a população
}

dialelo <- function(pais, populacao, prole, h2, gmd) { # argumentos da função
  combinacoes <- combinations(length(pais), 2, v = pais) # cria todas as combinações híbridas
  progenies <- vector('list', length = nrow(combinacoes)) # cria a estrutura das progênies
  id <- cbind(rep(c(1:nrow(combinacoes)), each = prole),
              rep(combinacoes[,1], each = prole),
              rep(combinacoes[,2], each = prole)) # cria a estrutura de identificação das progênies
  g <- ncol(populacao) - 5 # verifica o número de genes envolvidos

#  for(i in 1:2){
  
 for(i in 1:nrow(combinacoes)) { # 'laço' -  para executar todos cruzamentos
    cruz.i <- cruz(populacao[combinacoes[i,1],4:(3 + g)],
                   populacao[combinacoes[i,2],4:(3 + g)],prole) # aplica a função 'cruz'
    progenies[[i]] <- cruz.i # armazena a prole de cada cruzamento
  }
  dialelo <- as.data.frame(do.call(rbind, progenies)) # monta a planilha com o dialelo
  vg <- apply(dialelo, 1, soma.vg, gmd = gmd) # calcula os valores genotípicos totais por indivíduo
  varg <- var(vg) # calcula a variância genética
  fen <- vg + rnorm(nrow(dialelo), mean = 0, sqrt(varg * (1 - h2) / h2)) # atribui o desvio fenotípico
  resp <- cbind(id, dialelo, vg, fen) # monta a planilha com os resultados do dialelo
  names(resp) <- c('fami', 'pai', 'mae', 1:g, 'vg', 'fen') # atribui os nomes as colunas da planilha
  return(resp) # retorna o dialelo
}

soma.vg <- function(indv, gmd){ # argumentos da função
  h <- gmd # atribui o desvio de dominância ao heterozigoto
  gen <- numeric() # cria um vetor qualquer - armazena os valores genotípicos
  for(i in 1:length(indv)){ # 'laço' - passa a regra p/ cada loco
    # se o loco for igual a 'r' (recessivo) valor genotípico recebe -1
    # se o loco for igual a 'h' (heterozigoto) valor genotípico recebe gmd
    # se o loco for igual a 'd' (dominante) valor genotípico recebe 1
    if(indv[i] == 'r') {
      gen[i] <- -1
    } else {
      if(indv[i] == 'h') gen[i] <- h else gen[i] <- 1
    }
  }
  return(sum(gen)) # retorna a soma dos valores genotípcos por loco
}

# segregacao <- function(pai, mae) { # argumentos da função
#   filho <- numeric() # cria um vetor qualquer - armazena os valores genotípicos do filho
#   for(i in 1:length(pai)) { # 'laço' - passa a regra de segregação p/ cada loco
# se ambos pais são homozigotos dominantes (d), filho é dominante (d)
# se um pai é dominante e o outro é heterozigoto (h) filho pode ser 'h' ou 'd' 50% p/ cada
# se ambos pais são heterozigoto segregação 1:2:1 'd', 'h', e 'r'
# se um pai é dominante e o outro é recessivo filho é heterozigoto
# se um pai é reccessivo e o outro é heterozigoto filho pode ser 'r' ou 'r' 50% p/ cada
# se ambos pais são recessivos filho é recessivo
# if(pai[i] == 'd' & mae[i] == 'd') filho[i] <- 'd'
#     if(pai[i] == 'd' & mae[i] == 'h' | pai[i] == 'h' & mae[i] == 'd'){
#       u <- runif(1)
#       if(u <= .5) filho[i] <- 'd' else filho[i] <- 'h'
#     }
#     if(pai[i] == 'd' & mae[i] == 'r' | pai[i] == 'r' & mae[i] == 'd'){
#       filho[i] <- 'h'
#     }
#     if(pai[i] == 'h' & mae[i] == 'h'){
#       u <- runif(1)
#       if(u <= .25) filho[i] <- 'd' else
#       if((u > .25) & (u <= .75)) filho[i] <- 'h' else
#       filho[i] <- 'r'
#     }
#     if(pai[i] == 'h' & mae[i] == 'r' | pai[i] == 'r' & mae[i] == 'h'){
#       u <- runif(1)
#       if(u <= .5) filho[i] <- 'h' else filho[i] <- 'r'
#     }
#     if(pai[i] == 'r' & mae[i] == 'r') filho[i] <- 'r'
#   }
#   return(filho) # retorna a constituição genética do filho do cruzamento
# }
# 


# _____________________________________________________________________
#_ OK
#_ Função para fazer cruzamentos
#_ Parâmetros:    pai e mae = indvíduos, 'pai' e 'mãe'
#                     prole = número de indivíduos na progênie    
#______________________________________________________________________

# cruz<-function(pai, mae, prole){ # argumentos da função
#   progenie <- matrix(NA, ncol = length(pai), nrow = prole) # cria matriz qualquer [prole,g] - armazena a prole
#   for(i in 1:prole) { #'laço' - cria prole filhos do cruzamento
#     progenie[i,] <- segregacao(pai, mae)
#   }
#   return(progenie) # retorna a prole do cruzamento
# }
# 

segregacao.fmc <- function(pai,mae){

	filho <- rep(0,length(pai))
	n <- length(pai)
	teste <- .C("segregacao", as.integer(pai), as.integer(mae), filho2=as.integer(filho), as.integer(n))$filho2
	teste[teste==1] <- 'd'
	teste[teste==2] <- 'h'
	teste[teste==3] <- 'r'
	return(teste)
}

cruz<-function(pai, mae, prole){ # argumentos da função
  progenie <- matrix(0, ncol = length(pai), nrow = prole) # cria matriz qualquer [prole,g] - armazena a prole
  prog2 <- apply(progenie,2,function(x) segregacao.fmc(pai,mae))
  return(prog2) # retorna a prole do cruzamento
}

p0 <- pop.inicial(p, g, n, h2, gmd)
pais.teste <- sort(sample(c(1:n), 2*int , replace = FALSE))
system.time(teste.prog <- dialelo(pais.teste, p0, prole, h2, gmd))