彭哲也 謝民主

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

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

彭哲也 謝民主*

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

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

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

0 引 言

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

將p個(gè)有區(qū)別的球放入q個(gè)相同的盒子中，要求所有盒子都非空的方案數(shù)為第二類Stirling數(shù)，記作S(p,q)[1]。例如，將編號(hào)為A、B、C、D的4個(gè)小球放入2個(gè)相同的盒子里的所有不同放法可以表示如下(盒子用集合表示)：[{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時(shí)，第二類Stirling可由以下遞歸算法求解:考慮第p個(gè)球，該球不是獨(dú)占一個(gè)盒子就是和其他球在同一個(gè)盒子中。第p個(gè)球獨(dú)占一個(gè)盒子：其他p-1個(gè)球分配到其他q-1個(gè)非空盒子中的所有方案數(shù)為S(p-1,q-1)；第p個(gè)球不獨(dú)占一個(gè)盒子的情況則是：其他p-1個(gè)球放入q個(gè)非空的盒子中，第p個(gè)球放入q個(gè)盒子的任意一個(gè)中，這樣共有qS(p-1,q)種方法。這樣就得到了第二類Stirling數(shù)的遞推公式：

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

(1)

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

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

(2)

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

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

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

1 Stirling圖

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

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

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

圖1 Stirling圖

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

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

對(duì)Stirling圖，本文定義如下操作：

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

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

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

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

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

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

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

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

相關(guān)操作實(shí)現(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)計(jì)算過，返回計(jì)算值

returnM[j][k].count;

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

//終端結(jié)點(diǎn)

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個(gè)球分配到標(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é)點(diǎn)，輸出P(1..k) = 1..k

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

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

return; }

if (k = 1) {

//終端結(jié)點(diǎn)，輸出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é)點(diǎn)，輸出P(1..k) = 1..k

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

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

if (k = 1) {

//終端結(jié)點(diǎn)，輸出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];

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

//點(diǎn)的

// k+1條有向邊的采樣個(gè)數(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){

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

floatrand = random();

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

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

n_left--;}

}

對(duì)于盒子允許為空的情況，利用Stirling圖枚舉出所有的分配方案和按照均勻分布隨機(jī)采樣n個(gè)不同的方案的算法也容易實(shí)現(xiàn)。在盒子允許為空的情況下，令p個(gè)不同的球分配到q個(gè)相同的盒子中的所有不同方案的集合為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)需要按照均勻分布隨機(jī)采樣輸出n個(gè)不同的方案時(shí)，可按照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 實(shí)驗(yàn)結(jié)果

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

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

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

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

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

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

當(dāng)方案總數(shù)過大，枚舉所有的方案不可行時(shí)，我們?cè)谡麄€(gè)空間進(jìn)行按照均勻分布隨機(jī)采樣100萬(wàn)個(gè)不同的方案進(jìn)行輸出，表2顯示對(duì)于不同的(p,q)，隨機(jī)采樣算法randomSample(p,q,1 000 000)的運(yùn)行時(shí)間。

表2 randomSample(p,q,1 000 000)運(yùn)行時(shí)間

4 結(jié) 語(yǔ)

本文提出一種新的數(shù)據(jù)結(jié)構(gòu)Stirling圖，實(shí)現(xiàn)了該數(shù)據(jù)結(jié)構(gòu)的操作算法，能枚舉出不同的球分配到相同盒子中的所有不同方案。當(dāng)方案總數(shù)過多無(wú)法枚舉時(shí)，本文實(shí)現(xiàn)了隨機(jī)采樣操作算法，該算法按照均勻分布隨機(jī)采樣輸出指定個(gè)數(shù)的不同的分配方案。測(cè)試結(jié)果表明，這些算法能有效解決球不同盒子相同的“分球入盒”方案枚舉問題和均勻采樣問題。實(shí)際應(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] 劉運(yùn)龍, 王建新. 3-維匹配問題的一種固定參數(shù)枚舉算法[J].計(jì)算機(jī)科學(xué),2010,37(5):210-213.

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

[10] 張景云. 改進(jìn)的堆的枚舉算法的研究[J]. 計(jì)算機(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。國(guó)家自然科學(xué)基金項(xiàng)目(61370172)。彭哲也，碩士生，主研領(lǐng)域：生物信息學(xué)。謝民主，教授。