我在每个人群中都有np种群和Ne个体.每个人由两个配子(男性和女性)组成.每个配子包含三个基因.每个gen可以是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)
不幸的是,代码仍然太慢.任何人都可以帮助我加快速度吗？
这段代码：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中.然后,我们立即为所有人群提取所有这些,并使用Flatten和Partition来恢复人口结构.然后,我们使用Tally计算频率,追加频率为零的缺失配子(由Join [#,Thread [{Complement [Range [0,7],#[[All,1]]],0}]] line完成,并为固定人口排序每个频率列表.最后,我们提取频率并丢弃配子十进制索引.
自从在打包数组上执行以来,所有操作都非常快.加速是由于问题的矢量化公式和打包阵列的使用.它的内存效率也更高.