Оптимизация расчета частот гамет в популяциях - PullRequest
Новая
       12

Оптимизация расчета частот гамет в популяциях

5 голосов
/ 25 августа 2011

Мне нужно оптимизировать расчет частоты гамет в популяциях.

У меня есть np популяции и Ne особей в каждой популяции.Каждый индивид состоит из двух гамет (мужского и женского пола).Каждая гамета содержит три гена.Каждый ген может быть 0 или 1.Таким образом, каждый человек представляет собой матрицу 2х3.Каждый ряд матрицы представляет собой гамету, данную одним из родителей.Набор особей в каждой популяции может быть произвольным (но всегда длиной Ne).Для простоты начальные популяции с индивидуумами могут быть представлены как: 

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

Полный набор всех возможных гамет: 

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

Каждый индивид может сгенерировать гамету по 8 возможным путям с равной вероятностью.Эти гаметы: Tuples@Transpose@ind[[iPop, iInd]] (где iPop и iInd - индексы населения и отдельных в этой популяции).Мне нужно рассчитать частоты гамет, генерируемых индивидуумами для каждой популяции.

На данный момент мое решение заключается в следующем.

Сначала я конвертирую каждого индивида в гаметы, которые он может произвести: 

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

Но более эффективный способ сделать это: 

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

Во-вторых, я вычисляючастоты произведенных гамет, включая нулевые частоты для гамет, которые возможны, но отсутствуют среди населения: 

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

Более эффективная версия этого кода: 

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

К сожалению, код все еще слишкоммедленный.Кто-нибудь может помочь мне ускорить это?

1 Ответ

6 голосов
/ 25 августа 2011

Этот код: 

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)
    ]
]

использует преобразование в десятичные числа и упакованные массивы и ускоряет ваш код в 40 раз на моей машине. Тесты: 

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}}}

Обратите внимание, что в целом (для случайных групп населения) порядок частот в вашем и моих решениях по какой-то причине различен, и 

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

Теперь некоторые пояснения: сначала мы создаем таблицу t, которая построена следующим образом: каждой гамете присваивается число от 0 до 7, которое соответствует нулям и единицам в ней, которые рассматриваются как двоичные цифры. Затем в таблице содержатся возможные гаметы, произведенные индивидуумом, которые хранятся в позиции {i,j}, где i - десятичная дробь для, скажем, мамы матери, и j - для отцов, для этого индивида (каждый индивид уникален идентифицируется парой {i,j}). Гамет, произведенный индивидуумом, также конвертируется в десятичные. Вот как это выглядит: 

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}}}

Очень важным (решающим) шагом является преобразование этой таблицы в упакованный массив.

Строка Flatten[ind.(2^Range[0, 2]) + 1, 1]] преобразует гаметы родителей из двоичного в десятичное для всех индивидуумов сразу во всех популяциях и добавляет 1, так что они становятся индексами, при которых список возможных для производства гамет хранится в таблице t для данного человека. Затем мы Extract всех их сразу, для всех групп населения, и используем Flatten и Partition для восстановления структуры населения. Затем мы вычисляем частоты с Tally, добавляем отсутствующие гаметы с нулевыми частотами (по линии Join[#, Thread[{Complement[Range[0, 7], #[[All, 1]]], 0}]]) и Sort каждый список частот для фиксированной популяции. Наконец, мы извлекаем частоты и отбрасываем десятичный индекс гаметы.

Все операции выполняются довольно быстро, поскольку выполняются над упакованными массивами. Ускорение связано с векторизацией постановки задачи и использованием упакованных массивов. Это также намного более эффективно для памяти.

...