performance - 优化种群中配子频率的计算

我需要优化种群中配子频率的计算。

我在每个人口中都有np人口和Ne人口。每个人由两个配子(雄性和雌性)组成。每个配子包含三个基因。每个gen可以是0或1。因此，每个人都是2x3矩阵。矩阵的每一行都是由 parent 之一给出的配子。每个人口中的个体集合可以是任意的(但始终为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)

不幸的是，代码仍然太慢。有人可以帮助我加快速度吗？

最佳答案

这段代码:

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]]行将所有种群中所有个体的 parent 配子从二进制立即转换为十进制，并加1，这样它们就成为索引，在该索引中可能产生的配子列表存储在给定个体的表t中。然后，我们一次对所有人口进行Extract所有编码，并使用Flatten和Partition恢复人口结构。然后，我们使用Tally计算频率，将频率为零的缺失配子追加(通过Join[#, Thread[{Complement[Range[0, 7], #[[All, 1]]], 0}]]行完成)，并使用Sort固定人口的每个频率列表。最后，我们提取频率并丢弃配子十进制索引。

自从对压缩数组执行操作以来，所有操作都非常快。加速是由于问题的矢量化表示以及压缩数组的使用。它还具有更多的内存-高效。

关于performance - 优化种群中配子频率的计算，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/7187329/