Optimizar el cálculo de frecuencias de gametos en poblaciones

Question

Necesito optimizar el cálculo de las frecuencias de gametos en poblaciones.

Tengo np poblaciones y Ne personas en cada población. Cada individuo está formado por dos gametos (masculino y femenino). Cada gameto contiene tres genes. Cada gen puede ser 0 o 1. Entonces cada individuo es una matriz de 2x3. Cada fila de la matriz es un gameto dado por uno de los padres. El conjunto de individuos en cada población puede ser arbitrario (pero siempre de Ne de longitud). Por simplicidad poblaciones iniciales con individuos pueden ser dados como:

Ne = 300; np = 3^7; 
(*This table may be arbitrary with the same shape*) 
ind = Table[{{0, 0, 0}, {1, 1, 1}}, {np}, {Ne}] 

conjunto completo de todos los gametos posibles:

allGam = Tuples[{0, 1}, 3] 

cada individuo puede generar un gameto por 8 posibles formas con la misma probabilidad. Estos gametos son: Tuples@Transpose@ind[[iPop, iInd]] (donde iPop y iInd - índices de población y de individuos en esa población). Necesito calcular las frecuencias de gametos generadas por individuos para cada población.

En este momento mi solución es la siguiente.

Al principio, me convierto cada individuo en gametos se puede producir:

gamsInPop = Map[Sequence @@ Tuples@Transpose@# &, ind, {2}] 

Pero más eficiente manera de hacer esto es:

gamsInPop = 
Table[Join @@ Table[Tuples@Transpose@ind[[i, j]], {j, 1, Ne}], {i, 1, np}] 

En segundo lugar, puedo calcular las frecuencias de gametos producidos incluyendo frecuencias cero para los gametos que son posibles pero ausentes en la población:

gamFrq = Table[Count[pop, gam]/(8 Ne), {pop, gamInPop}, {gam, allGam}] 

Versión más eficiente de este código:

gamFrq = Total[ 
    Developer`ToPackedArray[ 
    gamInPop /. Table[ 
     allGam[[i]] -> Insert[{0, 0, 0, 0, 0, 0, 0}, 1, i], {i, 1, 
     8}]], {2}]/(8 Ne) 

Desafortunadamente, el código es aún demasiado lento. ¿Alguien puede ayudarme a acelerarlo?

Answer 1

este código:

Clear[getFrequencies]; 
Module[{t = 
    Developer`ToPackedArray[ 
    Table[FromDigits[#, 2] & /@ 
     Tuples[Transpose[{ 
      PadLeft[IntegerDigits[i, 2], 3], 
      PadLeft[IntegerDigits[j, 2], 3]}]], 
     {i, 0, 7}, {j, 0, 7}] 
    ]}, 
    getFrequencies[ind_] := 
    With[{extracted = 
     Partition[ 
      Flatten@Extract[t, Flatten[ind.(2^Range[0, 2]) + 1, 1]], 
      Ne*8]}, 
     Map[ 
     Sort@Join[#, Thread[{Complement[Range[0, 7], #[[All, 1]]], 0}]] &@Tally[#] &, 
     extracted 
     ][[All, All, 2]]/(Ne*8) 
    ] 
] 

utiliza la conversión de números decimales y las matrices envasados, y acelera su código por un factor de 40 en mi máquina. Los puntos de referencia:

In[372]:= Ne=300;np=3^7; 
(*This table may be arbitrary with the same shape*) 
inds=Table[{{0,0,0},{1,1,1}},{np},{Ne}]; 

In[374]:= 
getFrequencies[inds]//Short//Timing 
Out[374]= {0.282,{{1/8,1/8,1/8,1/8,1/8,1/8,1/8,1/8},<<2185>>, 
{1/8,1/8,1/8,1/8,1/8,1/8,1/8,1/8}}} 

In[375]:= 
Timing[ 
    gamsInPop=Table[Join@@Table[Tuples@Transpose@inds[[i,j]],{j,1,Ne}],{i,1,np}]; 
    gamFrq=Total[Developer`ToPackedArray[gamsInPop/.Table[allGam[[i]]-> 
     Insert[{0,0,0,0,0,0,0},1,i],{i,1,8}]],{2}]/(8 Ne)//Short] 

Out[375]= {10.563,{{1/8,1/8,1/8,1/8,1/8,1/8,1/8,1/8},<<2185>>, 
    {1/8,1/8,1/8,1/8,1/8,1/8,1/8,1/8}}} 

Tenga en cuenta que en general (para las poblaciones al azar), el orden de las frecuencias en sus y mis soluciones son por alguna razón diferente, y

In[393]:= fr[[All,{1,5,3,7,2,6,4,8}]] == gamFrq 
Out[393]= True 

Ahora, algunas explicaciones: en primer lugar, creamos una tabla t, que se construye de la siguiente manera: a cada gameto se le asigna un número del 0 al 7, que corresponde a los ceros y los que se tratan como dígitos binarios. La tabla tiene los posibles gametos producidos por un individuo, almacenados en una posición {i,j}, donde i es un decimal para gametos de la madre (por ejemplo) y j - para padres, para esa persona (cada individuo está identificado únicamente por un par {i,j}) . Los gametos producidos por individuo también se convierten a decimales.Aquí es como se ve:

In[396]:= t//Short[#,5]& 
Out[396]//Short= {{{0,0,0,0,0,0,0,0},{0,1,0,1,0,1,0,1},{0,0,2,2,0,0,2,2}, 
{0,1,2,3,0,1,2,3},{0,0,0,0,4,4,4,4},{0,1,0,1,4,5,4,5},{0,0,2,2,4,4,6,6}, 
{0,1,2,3,4,5,6,7}},<<6>>,{{7,6,5,4,3,2,1,0},{7,7,5,5,3,3,1,1},{7,6,7,6,3,2,3,2}, 
<<2>>,{7,7,5,5,7,7,5,5},{7,6,7,6,7,6,7,6},{7,7,7,7,7,7,7,7}}} 

A (importante) paso muy importante es convertir esta tabla a un conjunto empaquetado.

La línea Flatten[ind.(2^Range[0, 2]) + 1, 1]] convierte gametos de los padres de binario a decimal para todas las personas a la vez, en todas las poblaciones, y se suma 1 a fin de que éstas se convierten en índices en los que se almacena la lista de posible producir gametos en una tabla t para una individuo dado Luego, Extract todos a la vez, para todas las poblaciones, y usamos Flatten y Partition para recuperar la estructura de la población. Luego, calculamos frecuencias con Tally, anexamos gametos faltantes con frecuencias cero (hechas por Join[#, Thread[{Complement[Range[0, 7], #[[All, 1]]], 0}]] línea), y Sort cada lista de frecuencias para una población fija. Finalmente, extraemos las frecuencias y descartamos el índice decimal de gametos.

Todas las operaciones son bastante rápidas ya que se realizan en matrices empaquetadas. La aceleración se debe a la formulación vectorizada del problema y al uso de matrices empaquetadas. También es mucho más eficiente con la memoria.