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

彭哲也 謝民主

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

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

彭哲也 謝民主*

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

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

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

0 引 言

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

將p個有區別的球放入q個相同的盒子中，要求所有盒子都非空的方案數為第二類Stirling數，記作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數在數論、組合論、代數拓撲、算法最優化等領域都有相關的應用[2]。顯然，當pq且q>1時，第二類Stirling可由以下遞歸算法求解:考慮第p個球，該球不是獨占一個盒子就是和其他球在同一個盒子中。第p個球獨占一個盒子：其他p-1個球分配到其他q-1個非空盒子中的所有方案數為S(p-1,q-1)；第p個球不獨占一個盒子的情況則是：其他p-1個球放入q個非空的盒子中，第p個球放入q個盒子的任意一個中，這樣共有qS(p-1,q)種方法。這樣就得到了第二類Stirling數的遞推公式：

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

(1)

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

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

(2)

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

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

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

1 Stirling圖

本文用一個劃分函數(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)中不同方案的個數|S(p,q)|=S(p,q)。

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

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單獨放到一個盒子中，且把該盒子貼上標簽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)表示，每個盒子由此過程獲得了不同的標簽)。然后將第p個球放到標簽為i(i=1,…,q)的盒子中(用P(p)=1..q表示)的所有不同方案的集合。為了使用式(3)構造出分配方案的集合S(p,q)，本文引入如圖1所示的Stirling圖。

圖1 Stirling圖

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

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

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

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

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

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

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

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

2 數據結構及算法實現

本文用一個(p+1)×(q+1)的二維數組M(假設下標是從0開始)來存儲Stirling圖的結點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條指向其右孩子的邊。

相關操作實現的偽代碼如下所示:

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)

//已經計算過，返回計算值

returnM[j][k].count;

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

//終端結點

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個球分配到標簽為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) {

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

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

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

return; }

if (k = 1) {

//終端結點，輸出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;

//恢復初始值

}

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) {

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

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

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

if (k = 1) {

//終端結點，輸出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];

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

//點的

// k+1條有向邊的采樣個數

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;

//恢復初始值

}

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

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

//采樣率

令count(j-1, k-1)*r的整數部分為a1，小數部分為e1，count(j, k-1)*r的整數部分為a2，小數部分為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]之間均勻分布的隨機數rand

floatrand = random();

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

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

n_left--;}

}

對于盒子允許為空的情況，利用Stirling圖枚舉出所有的分配方案和按照均勻分布隨機采樣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); }

}

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

3 實驗結果

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

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

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

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

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

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

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

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

4 結 語

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

[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數及其推廣[D]. 大連理工大學,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-維匹配問題的一種固定參數枚舉算法[J].計算機科學,2010,37(5):210-213.

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

[10] 張景云. 改進的堆的枚舉算法的研究[J]. 計算機應用與軟件,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。國家自然科學基金項目(61370172)。彭哲也，碩士生，主研領域：生物信息學。謝民主，教授。