彭哲也 謝民主

(湖南師范大學(xué)物理與信息科學(xué)學(xué)院 湖南 長沙 410081)

一種基于Stirling圖枚舉算法的分球入盒問題求解

彭哲也 謝民主*

(湖南師范大學(xué)物理與信息科學(xué)學(xué)院 湖南 長沙 410081)

已有的分球入盒問題解法通常只關(guān)注分球的總方案數(shù)，目前尚沒有公開的計算機算法來枚舉出所有具體的分球方案，而方案的枚舉是生物信息學(xué)中一些分區(qū)優(yōu)化算法的基礎(chǔ)。受第二類Stirling數(shù)的遞推公式的啟發(fā)，提出一個新的數(shù)據(jù)結(jié)構(gòu)——Stirling圖。在此基礎(chǔ)上設(shè)計一個算法來枚舉p個不同球分配到q個相同盒子里的所有不同的方案。當(dāng)p和q較大，枚舉出所有的方案不可行時，設(shè)計另一個算法在整個方案空間實現(xiàn)均勻采樣，輸出指定個數(shù)的方案。測試結(jié)果表明，這些算法在內(nèi)存為8 GB的普通PC上可在合理的時間內(nèi)枚舉出上百萬組不同的方案。

分球入盒問題 第二類Stirling數(shù) 枚舉算法 Stirling圖 均勻采樣

0 引 言

分球入盒是組合數(shù)學(xué)中的經(jīng)典問題，在組合優(yōu)化中有廣泛的應(yīng)用。把p個球放入q個盒子里，根據(jù)球和盒子是否相同可以分為四種不同的模型：(1) 球相同盒子相同；(2) 球相同盒子不同；(3) 球不同盒子不同；(4) 球不同盒子相同。每種模型又進一步分為盒子是否允許為空這兩種情形。本文主要討論球不同盒子相同的這類模型。

將p個有區(qū)別的球放入q個相同的盒子中，要求所有盒子都非空的方案數(shù)為第二類Stirling數(shù)，記作S(p,q)[1]。例如，將編號為A、B、C、D的4個小球放入2個相同的盒子里的所有不同放法可以表示如下(盒子用集合表示)：[{A},{B,C,D}]，[{B},{A,C,D}]，[{C},{A,B,D}]，[{D},{A,B,C}]，[{A,B},{C,D}]，[{A,C},{B,D}]，[{A,D},{B,C}]。因此S(4,2)=7。

第二類Stirling數(shù)在數(shù)論、組合論、代數(shù)拓?fù)?、算法最?yōu)化等領(lǐng)域都有相關(guān)的應(yīng)用[2]。顯然，當(dāng)pq且q>1時，第二類Stirling可由以下遞歸算法求解:考慮第p個球，該球不是獨占一個盒子就是和其他球在同一個盒子中。第p個球獨占一個盒子：其他p-1個球分配到其他q-1個非空盒子中的所有方案數(shù)為S(p-1,q-1)；第p個球不獨占一個盒子的情況則是：其他p-1個球放入q個非空的盒子中，第p個球放入q個盒子的任意一個中，這樣共有qS(p-1,q)種方法。這樣就得到了第二類Stirling數(shù)的遞推公式：

S(p,q)=S(p-1,q-1)+qS(p-1,q)

(1)

將p個有區(qū)別的球放入q個相同的盒子中, 當(dāng)盒子允許為空時，令其不同的方法數(shù)為F(p,q)，顯然：

F(p,q)=S(p,1)+S(p,2)+…+S(p,q)

(2)

第二類Stirling數(shù)的計算已經(jīng)得到了廣泛深入的研究[3]，但是目前尚沒有公開的計算機算法來枚舉出所有具體的分球方案。枚舉是計算復(fù)雜性理論中的一類基本運算，是理論計算機科學(xué)中一個重要分支。自從1979 年 Valiant 首次提出枚舉運算以來[6]，人們對枚舉運算進行了大量的研究[7]，其中最典型的枚舉方法有枚舉問題實例的所有解或統(tǒng)計問題實例中解的個數(shù)[8]。

近年來，隨著DNA測序的快速發(fā)展，利用多倍體基因組的DNA測序片段如何組裝出構(gòu)成多倍體的多個單體型[11]和利用環(huán)境群體微生物的深度DNA測序序列重構(gòu)該群體的微生物基因組[12]等問題的研究成為生物信息學(xué)的研究熱點。這些問題的數(shù)學(xué)建模通常是把DNA測序序列分配到多個分區(qū)中以滿足特定的優(yōu)化目標(biāo)。因此，分球入盒方案的枚舉成為一些求解這類問題優(yōu)化算法的基本操作[11]。

本文受第二類Stirling數(shù)的遞推公式的啟發(fā)提出了一個新的數(shù)據(jù)結(jié)構(gòu)——Stirling圖。在此基礎(chǔ)上設(shè)計一個算法來枚舉p個不同球分配到q個相同盒子里的所有不同的方案；當(dāng)p和q較大時，不可能枚舉出所有的方案，該算法則在整個方案空間實現(xiàn)均勻采樣，輸出指定個數(shù)的方案。在內(nèi)存為8 GB、主頻為2.5 GHz的普通PC上的測試結(jié)果表明，該算法可以在幾秒內(nèi)枚舉出上百萬組的不同的方案。

1 Stirling圖

本文用一個劃分函數(shù)(Partition function)P(i)：{1,…,p}→{1,…,q}來表示p個不同的球分配到q個盒子中的一個具體方案。例如P(i)|i=1,2,3=1表示將編號為1，2，3的球放入第一個盒子里,P(4)=2表示將編號為4的球放入第二個盒子里。為了簡便，下面用P(1..k)=1表示P(i)|i=1,…,k=1;P(1..k)=1..k表示P(i)|i=1,…,k=i。令p個不同的球分配到q個非空的相同盒子中的所有不同的方案的集合為S(p,q)，集合S(p,q)中不同方案的個數(shù)|S(p,q)|=S(p,q)。

當(dāng)p0時，S(p,q)=1，即S(p,q)中只有一種分配方案，不同的球放到不同的盒子里，給球i所在盒子貼上標(biāo)簽i，則這種分配方案用P(1..p) = 1..p表示；當(dāng)p>q=1，這時S(p,q)中同樣只存在一種分配方案，即所有的球放在唯一的盒子中，給該盒子貼上標(biāo)簽1，則這種分配方案用P(1..p) = 1來表示；當(dāng)p>q>1時，所有不同分配方案的集合S(p,q)可以用式(1)的證明過程構(gòu)造出來：

S(p,q) ={P(p)=q·S(p-1,q-1)}∪

{S(p-1,q)·(P(p)=1…q)}

(3)

其中{(P(p)=q)·S(p-1,q-1)}表示將球p單獨放到一個盒子中，且把該盒子貼上標(biāo)簽q(用P(p)=q表示)，然后把其他p-1個球分配到其他q-1個相同非空的盒子中(用S(p-1,q-1)表示)的所有不同方案的集合；而{(P(p)=1..q)·S(p-1,q-1)}則表示將前p-1個球分配到q個相同的非空盒子中(用S(p-1,q)表示，每個盒子由此過程獲得了不同的標(biāo)簽)。然后將第p個球放到標(biāo)簽為i(i=1,…,q)的盒子中(用P(p)=1..q表示)的所有不同方案的集合。為了使用式(3)構(gòu)造出分配方案的集合S(p,q)，本文引入如圖1所示的Stirling圖。

圖1 Stirling圖

Stirling圖是一個有向無環(huán)圖，其結(jié)點用S(j,k)來標(biāo)記(j≥k≥1)。Stirling圖的結(jié)點分為兩類，一類是終端結(jié)點S(k,k)和S(j,1)，一類是非終端結(jié)點S(j,k)(j>k>1)。非終端結(jié)點S(j,k)結(jié)點有左右兩個孩子結(jié)點，其左孩子結(jié)點為S(j-1,k-1)，右孩子結(jié)點為S(j-1,k)。圖中用一條標(biāo)簽為P(j)=k的有向邊由非終端結(jié)點S(j,k)指向其左孩子結(jié)點，用標(biāo)簽分別為P(j)=1,P(j)=2,…,P(j)=k的k條有向邊指向其右孩子結(jié)點。對于兩個結(jié)點S和S′，如果存在一條路徑從S到S′，則S為S′的先驅(qū)結(jié)點，S′為S的后繼結(jié)點。

在Stirling圖中從非終端結(jié)點S(j,k)到其后繼終端結(jié)點的一條路徑代表把p個不同的球分配到q個非空的相同盒子的一種分配方案。如從S(3,2)到S(2,1)的路徑代表一個分配方案：P(1..2)=1,P(3)=2;即把前兩個球放到一個盒子里，第三個球放到單獨的盒子了。由式(2)可知，從非終端結(jié)點S(j,k)到其所有后繼終端結(jié)點的路徑就代表把j個不同的球分配到k個非空相同盒子的所有不同分配方案。

對Stirling圖，本文定義如下操作：

1)initiate(p,q):構(gòu)造出一個Stirling圖，包含結(jié)點S(p,q)及其所有后繼結(jié)點，及相應(yīng)的有向邊。

2)destroy():刪除Stirling圖，釋放內(nèi)存。

3)count(j,k): 計算出S(j,k)，即|S(p,q)|。

4)pathTransverse(j,k):對Stirling圖遍歷從結(jié)點S(j,k)開始以任意終端節(jié)點結(jié)束的所有路徑，即枚舉出集合S(j,k)中的所有方案。

5)randomSample(j,k,n): 從集合S(p,q)中按照均勻分布隨機選擇n個不同的方案輸出。

2 數(shù)據(jù)結(jié)構(gòu)及算法實現(xiàn)

本文用一個(p+1)×(q+1)的二維數(shù)組M(假設(shè)下標(biāo)是從0開始)來存儲Stirling圖的結(jié)點S(j,k):0

為了支持pathTransverse和count操作，元素Mj,k包含兩個整型變量index和count，其中count用來記錄S(j,k)的值，初始化為-1；index的取值范圍為0到k，0表示進行pathTransverse操作時，下次訪問的邊是指向左孩子的邊，index=i(1≤i≤k)表示下次訪問的邊是第i條指向其右孩子的邊。

相關(guān)操作實現(xiàn)的偽代碼如下所示:

structnode{

int index = 0;

int count = -1;

}

node[][] M;

void initiate(p, q){M = new node[p+1][q+1];}

void destroy() { delete[] M;}

intcount(int j, int k) {

if (M[j][k].count > 0)

//已經(jīng)計算過，返回計算值

returnM[j][k].count;

if (j = k || k = 1) {

//終端結(jié)點

M[j][k].count = 1; return 1; }

M[j][k].count=count(j-1,k-1)+k*count(j-1,k);

returnM[j][k].count; }

voidpathTransverse(int j,int k){

// P記錄球分配方案，第i個球分配到標(biāo)簽為P[i]的盒

//子中

int[] P = new int[j+1];

output_pathes (j, k, P, j); }

voidoutput_pathes (int j, int k, int[] P, int balls){

if (j = k) {

//終端結(jié)點，輸出P(1..k) = 1..k

for(i= 1 .. k) { P[i] = i; }

輸出分配方案P[1], …, P[balls];

return; }

if (k = 1) {

//終端結(jié)點，輸出P(1..k) = 1

for(i= 1 .. j) { P[i] = 1; }

輸出分配方案P[1], …, P[balls];

return; }

if (M[j][k].index = 0) {

P[j] = k;

//訪問指向左孩子的邊

M[j][k].index++;

output_pathes(j-1, k-1, P, balls); }

//訪問左孩子

while(M[j][k].index ≤ k) {

//訪問指向右孩子的邊

P[j] = M[j][k].index; M[j][k].index++;

output_pathes(j-1, k, P, balls); }

//訪問右孩子

M[j][k].index=0;

//恢復(fù)初始值

}

void randomSample(int j, int k, int n){

//n必須小于count(j, k)

int[] P = new int[j+1];

sample_pathes(j, k, P, j, n);

}

voidsample_pathes(intj, intk, int[] P, int b, int n){

if(n< 1){ return; }

if (j = k) {

//終端結(jié)點，輸出P(1..k) = 1..k

for(i = 1 .. k) { P[i] = i; }

輸出分配方案P[1], …, P[b]; return; }

if (k = 1) {

//終端結(jié)點，輸出P(1..k) = 1

for(i= 1 .. j) { P[i] = 1; }

輸出分配方案P[1], …, P[b]; return; }

intn_left = n;

int[] n_a = new int[k+1];

//通過均勻分布隨機采樣，獲得通過其指向左右孩子結(jié)

//點的

// k+1條有向邊的采樣個數(shù)

get_sample_number(j, k, n, n_a)

if (M[j][k].index = 0) { P[j] = k;

//訪問指向左孩子的邊

M[j][k].index++;

sample_pathes(j-1, k-1, P, b, n_a[0]); }

//訪問左孩子

while(M[j][k].index < k) {

//訪問指向右孩子的邊

P[j] = M[j][k].index;

M[j][k].index++;

sample_pathes(j-1, k, P, b, n_a[P[j]]); }

//訪問右孩子

M[j][k].index=0;

//恢復(fù)初始值

}

voidget_sample_number(intj, intk, intn,int[] n_a){

float r = n/count(j, k);

//采樣率

令count(j-1, k-1)*r的整數(shù)部分為a1，小數(shù)部分為e1，count(j, k-1)*r的整數(shù)部分為a2，小數(shù)部分為e2；

n_a[0] = a1;

floatpl = e1/(e1+k*e2), pr = e2/(e1+k*e2);

for(i= 1.. k) {

n_a[i] = a2; }

intn_left = n - a1 - k*a2;

while(n_left> 0){

//得到一個在[0, 1]之間均勻分布的隨機數(shù)rand

floatrand = random();

if(rand ≤ pl) {n_a[0]++; }

else { intind = (rand - pl)/pr +1; n_a[ind]++; }

n_left--;}

}

對于盒子允許為空的情況，利用Stirling圖枚舉出所有的分配方案和按照均勻分布隨機采樣n個不同的方案的算法也容易實現(xiàn)。在盒子允許為空的情況下，令p個不同的球分配到q個相同的盒子中的所有不同方案的集合為F(p,q)，顯然有：

F(p,q)=S(p,q)∪S(p,2)∪…∪S(p,q)

(4)

枚舉出F(p,q)的算法如下所示:

F(int p, int q){

for(i = 1…q) { pathTransverse(p, i); }

}

當(dāng)需要按照均勻分布隨機采樣輸出n個不同的方案時，可按照get_sample_number的思路把n按照S(p,1),S(p,2),…,S(p,q)集合的大小劃分成n1,n2,…,nq，然后再依次調(diào)用randomSample(p,1,n1),randomSample(p,2,n2),…,randomSample(p,q,nq)即可。

3 實驗結(jié)果

本文算法用Java編程實現(xiàn)，測試計算機內(nèi)存為8 GB，CPU為Intel酷睿i5 3210M 2.5 GHz的PC機。圖2-圖4分別表示的是F(3,2)、F(4,2)、F(4,3)在控制臺輸出的枚舉結(jié)果。

圖2 F(3,2)枚舉結(jié)果

圖3 F(4,2)枚舉結(jié)果

圖4 F(4,3)枚舉結(jié)果

由于目前沒有類似的算法可以橫向比較，我們通過幾個用例來測試本文算法的性能。如表1所示，當(dāng)p和q增大時，F(xiàn)(p,q)枚舉出所有分配方案所需的時間基本上與方案總數(shù)成線性增長；當(dāng)p=12,q=5時，F(xiàn)(p,q)可以在約1秒鐘的時間內(nèi)枚舉出2百余萬種不同的方案。

表1 F(p,q)枚舉時間

當(dāng)方案總數(shù)過大，枚舉所有的方案不可行時，我們在整個空間進行按照均勻分布隨機采樣100萬個不同的方案進行輸出，表2顯示對于不同的(p,q)，隨機采樣算法randomSample(p,q,1 000 000)的運行時間。

表2 randomSample(p,q,1 000 000)運行時間

4 結(jié) 語

本文提出一種新的數(shù)據(jù)結(jié)構(gòu)Stirling圖，實現(xiàn)了該數(shù)據(jù)結(jié)構(gòu)的操作算法，能枚舉出不同的球分配到相同盒子中的所有不同方案。當(dāng)方案總數(shù)過多無法枚舉時，本文實現(xiàn)了隨機采樣操作算法，該算法按照均勻分布隨機采樣輸出指定個數(shù)的不同的分配方案。測試結(jié)果表明，這些算法能有效解決球不同盒子相同的“分球入盒”方案枚舉問題和均勻采樣問題。實際應(yīng)用中有很多優(yōu)化問題如多倍體單體型組裝[11]、宏基因組組裝[12]等均可轉(zhuǎn)化為“分球入盒”問題求解。本文提出的算法可以為這些問題的求解提供有效的分球入盒方案枚舉基本操作的支持。

[1] Wikipedia. Stirling numbers of the second kind[EB/OL].(2016-12-7)[2016-12-14].https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind.

[2] 王娟. 第二類Stirling數(shù)及其推廣[D]. 大連理工大學(xué),2009.

[3] Abramowitz M, Stegun I A. Handbook of mathematical functions with formulas, graphs, and mathematical tables (9th printing)[M]. New York: Dover, 1972:824-825.

[4] Boyadzhiev K N. Close Encounters with the Stirling Numbers of the Second Kind[J]. Mathematics Magazine,2012, 85(4):252-266.

[5] Bleick W W. Asymptotics of Stirling Numbers of the Second Kind[J]. Proceedings of the American Mathematical Society,1974,42(2):575-580.

[6] Valiant L.The complexity of enumeration and reliability problems[J].SIAM Journal on Computing, 1979,8(3):410-421.

[7] Flum J, Grohe M .The parameterized complexity of counting problems[J]. SIAM Journal on Computing, 2004,33(4):892-922.

[8] 劉運龍, 王建新. 3-維匹配問題的一種固定參數(shù)枚舉算法[J].計算機科學(xué),2010,37(5):210-213.

[9] 謝民主, 羅鋒, 唐烽. 單體型組裝最大片段割參數(shù)化精確算法[J].小型微型計算機系統(tǒng),2014,35(2):353-357.

[10] 張景云. 改進的堆的枚舉算法的研究[J]. 計算機應(yīng)用與軟件,2012,29(7):264-265,273.

[11] Xie Minzhu,Wu Qiong,Wang Jianxin, et al. H-PoP and H-PoPG: heuristic partitioning algorithms for single indivi-dual haplotyping of polyploids[J/OL]. Bioinformatics. [2016-08-16] Doi: 10.1093/bioinformatics/btw537.

[12] Rose R, Constantinides B, Tapinos A, et al. Challenges in the analysis of viral metagenomes[J/OL].Virus Evolution.[2016-08-03].DOI:http://dx.doi.org/10.1093/ve/vew022.

RESEARCHONDISTRIBUTINGBALLSINTOBOXESENUMERATIONALGORITHM

Peng Zheye Xie Minzhu*

(CollegeofPhysicsandInformationScience,HunanNormalUniversity,Changsha410081,Hunan,China)

The existing researches on the problem of distributing balls into boxes usually focus on the total number of different ways to distribute balls into boxes, but there are no public computer algorithms to enumerate them. However, enumerating them is the foundation to design some partition optimal algorithms in Bioinformatics. Inspired by the recursive formula of the Stirling numbers of the second kind, the paper proposes a new data structure—Stirling diagram, and based on the data structure, designs an algorithm to enumerate all different ways to distribute p different balls into q same boxes. When p and q are larger and none of the schemes is feasible, we design another algorithm to achieve uniform sampling of a given number of different ways. Test results show that these algorithms can enumerate millions of different distributing ways in a reasonable period of time on a PC with 8 GB memory.

Distributing balls into boxes problem Stirling numbers of the second kind Enumerating algorithm Stirling diagram Uniform sampling

TP306.1

A

10.3969/j.issn.1000-386x.2017.10.044

2016-12-15。國家自然科學(xué)基金項目(61370172)。彭哲也，碩士生，主研領(lǐng)域：生物信息學(xué)。謝民主，教授。