優化羣體中配子頻率的計算

Question

我需要優化羣體中配子頻率的計算。

我有np人口和Ne個人在每個人口中。每個人由兩個配子（男性和女性）組成。每個配子包含三個基因。每個gen可能是0或1。所以每個人都是2x3矩陣。矩陣的每一行都是由其中一位父母給出的配子。每個人羣中的一組人可能是任意的（但總是長度爲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) 

不幸的是，代碼仍然太慢。任何人都可以幫我加快速度嗎？

Answer 1

此代碼：

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的數字，這個數字對應於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。最後，我們提取頻率並丟棄配子的小數點。

由於在打包數組上執行，所有操作都相當快。加速是由於問題的矢量化表達和使用打包數組。它的記憶效率也更高。