[R-br] Res: Res: Implementação em C

[Fabio Mathias Corrêa]

Fernando,

Faltou apenas um pequeno exemplo de como inicia a simulação!

Valeu!!!

             Fábio Mathias Corrêa
        Departamento de Estatística
   Universidade Estadual de Santa Cruz




Tel.: 73-3680-5076
Cel.: 73-9991-8155




________________________________
De: Fernando Henrique Toledo <fernandohtoledo em gmail.com>
Para: r-br em listas.c3sl.ufpr.br
Enviadas: Segunda-feira, 21 de Março de 2011 8:43:36
Assunto: Re: [R-br] Res: Implementação em C

Senhores, bom dia.

A simulação que pretendo realizar implementei através das funções que estão 
listadas abaixo, a combinations é do pacote gtools, necessáriamente queria só 
melhorar o tempo de processamento do função dialelo, especificamente no for() 
que corre pela matriz de combinações. Fui rodando odas as etapas da função e 
detectei justamente nesse ponto o gargalo do calculo.

att,
Fernando H

# _____________________________________________________________________
#_ OK
#_ Função para cálculo do valor genotípico
#_ Parâmetros:  indv = indvíduo sob análise
#                gmd = grau médio de dominância
#______________________________________________________________________

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
}

# _____________________________________________________________________
#_ OK
#_ Função para regra da segregação
#_ Parâmetros:  pai e mãe = indvíduos sob análise (cruzando)
#______________________________________________________________________

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
}



>
>> # _____________________________________________________________________
>> #
>> #_ Função para fazer um dialelo
>> #_ Parâmetros:  pais = indivíduos selecionados na população anterior
>> #          população = população dos pais
>> #              prole = tamanho das progênies
>> #                 h2 = herdabilidade do caráter
>> #                gmd = grau médio de dominância
>> #______________________________________________________________________
>> 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: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
>> }
>


      
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20110321/5d5701c5/attachment.html>