ROSALIND | Longest Increasing Subsequence (LGIS)

發表於 更新於

給一正整數列，求其最長的遞增與遞減子數列。例如數列 5, 1, 4, 2, 3 的最長遞增數列為 1, 2, 3，最長遞減數列則為 5, 4, 3。

Given: A positive integer n≤10000 followed by a permutation π of length n.

Return: A longest increasing subsequence of π , followed by a longest decreasing subsequence of π.

(https://rosalind.info/problems/lgis/)

知識點

Synteny block 是物種基因體間基因順序彼此一致的片段。這些片段在演化過程中沒有因為重排（genome arrangement）而分離或反轉，暗示著坐落其中的基因不只位置相鄰，功能也可能彼此關聯。

由於基因體動輒包含數億對鹼基，即使比較的對象是親緣關係相近的物種，仍有可觀的計算量與複雜性。相較之下，synteny block 的數量遠少於鹼基數，卻同樣能對基因體的功能和結構提供一些線索。因此，透過比較 synteny block 來研究推演物種的演化過程也不失為一種可行的方式。

舉例來說，假設甲物種與乙物種共享 5 個 synteny blocks，並且依照它們在甲基因體的排列順序命名為 1, 2, 3, 4, 5。接著標記乙基因體上對應的 synteny blocks，便能概括兩者基因體的差異程度，或是推測基因體的結構變化。

1
2
甲：1, 2, 3, 4, 5
乙：3, 2, 1, 4, 5  <-- 可能發生一次轉置

儘管以符號標記 synteny blocks 可直觀呈現基因體差異，不過隨著比較對象和內容增加，可以使用距離指標來簡化基因體差異分析。


（人類與老鼠的染色體 synteny blocks 分布圖，https://en.wikipedia.org/wiki/Synteny）

在 “Counting Point Mutations” 當中已介紹 Hamming distance 來評估發生點突變的序列間的差異。在這題則引進最長遞增/遞減子數列的概念，透過 synteny block 的順序來定量基因體間的差異。最長遞增/遞減子數列的長度越長，表示兩基因體越相似；長度越短，則表示可能發生較多次基因體重排等結構上的變化。

題解

由於找最長遞增子數列的方式也能用在找最長遞減子數列，所以我將專注討論最長遞增子數列的解法。

假設數列皆如 1, 2, 3, 4,… 嚴格遞增，那麼尋找最長遞增數列並不困難。然而，對於數列元素順序凌亂的狀況，如 2, 5, 3, 1, 4, …，隨著數列長度增加，會越來越難直觀看出答案。

造成困難的主因在於，在檢查某項數字時，我們難以確定它是延續某個既有的遞增子數列，還是作為新的遞增子數列的開端。在檢查完所有數字之前，也無從得知當前檢查的數字是否為最終答案的一部份。

以數列 2, 5, 3, 1, 4 為例，2 能與 5 或 3 組成遞增子數列，但該怎麼判斷要納入何者，以組成最長的遞增子數列？又或在數列 3, 5, 1, 2, 4 的案例中，該依據什麼資訊來決定放棄包含 3, 5 的子數列，而轉為建立新的子數列 1, 2, 4？

綜上所述，可以將讓我們左右為難的因素歸納如下，

  1. 當前數字大於先前數字時，不確定要忽略它還是將其延伸至既有的遞增子數列。
  2. 當前數字小於先前數字時，不確定要忽略它還是以它為起點建立新的遞增子數列。

為了解決這問題，可以使用當前數列的最長遞增子數列長度為依據，來衡量納入當前數字和建立新數列的效益，假如行動能增加遞增子數列的長度則為之。

由於每次檢查都可能推翻先前的答案，所以要比對整個數列後才能得知最終結果。以數列 3, 5, 1, 2, 4 為例，由左至右檢查數字來構築最長遞增子數列：

檢查數字 當前最長遞增子數列 遞增數列長度 判斷 行為
3 [3] 1 起始數字 納入當前數字
5 [3, 5] 2 5 > 3，延長遞增數列 納入當前數字
1 [3, 5] 2 1 < 5，不能延長現有遞增數列 忽略當前數字
2 [3, 5], [1, 2] 2 2 > 1，形成新的遞增數列，與現有數列等長 建立當前數列
4 [1, 2, 4] 3 4 > 2，延長新建立的遞增數列，長度大於 [3, 5] 與 [1, 2] 納入當前數字

為了實現這套算法，需要紀錄以下資訊

  • lis_len：紀錄截至當前數列位置，最長遞增子數列的長度，用於衡量是否要納入新數字或建立新數列。
  • prev：該數字在遞增子數列的前一項的位置，用於追蹤遞增子數列的元素。
1
2
lis_len = [1] * len(l)    # 因為每個數字自成數列，所以在開始計算前，遞增子數列的長度至少為 1。
prev = [-1] * len(l)      # 在開始計算前，還沒有遞增子數列的位置資訊，所以預設為 -1。

假設現在檢查第 i 個數字，那麼我們需要把它與數列之前所有數字（l[j], i > j）比較，才能涵蓋到所有可能性。

1
2
3
4
5
6
7
for i in range(1, len(l)):
    for j in range(i):
        if l[i] > l[j] and lis_len[j] + 1 > lis_len[i]:
            lis_len[i] = lis_len[j] + 1
            prev[i] = j
        else:
            pass

簡言之，在當前數字大於先前數字時（l[i] > l[j]），有兩種可能的行動：

  • 將當前數字 l[i] 納入以 l[j] 結尾的遞增子數列，形成長度為 lis_len[j] + 1 的遞增子數列。
  • 忽略先前的數字，不更動 l[i] 所屬的遞增子數列，讓長度維持 lis_len[i]

為了找到最長的遞增子數列，應該選擇能增加子數列長度的行動。換句話說，若 lis_len[j] + 1 > lis_len[i]，那麼就該把 l[i] 納入以 l[j] 結尾的遞增子數列，並且更新 lis_len[i] = lis_len[j] + 1。

與此同時，也透過 prev 紀錄當前數字 l[i] 在遞增子數列的前一項的位置 prev[i] = j，以利在檢查全部數字之後，從輸入的數列重建出最長的遞增子數列。

透過這樣的比較，可以確保 lis_len 每一項都記錄了截至當前位置能獲得的最常遞增子數列長度。因此，lis_len 的最大值便是整個數列的最長遞增子數列長度；該值在 lis_len 的位置，也對應到遞增子數列最後一項的位置。

1
final_pos = max(range(len(l)), key=lambda x: lis_len[x])

由於prev陣列記錄了當前位置的數字在最長遞增子數列中前一項的索引位置,因此我們需要依序讀取prev中的值l[prev[i]],並將讀取出的數列反轉,才能獲得正確的最長遞增子數列順序。

由於 prev 記錄了當前位置的數字在最長遞增子數列中前一項的位置，所以回溯 prev 得到的數列順序與正確結果相反。仍需將讀取出的數列反轉，才能得到最長遞增子數列。

1
2
3
4
5
output = []
while final_pos != -1:
    output.append(l[final_pos])
    final_pos = prev[final_pos]
return output[::-1]

完整的功能實現如下，

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def lgis(l):
    """Longest Increasing Subsequence
    """
    lis_len = [1] * len(l)
    prev = [-1] * len(l)
    for i in range(1, len(l)):
        for j in range(i):
            if l[i] > l[j] and lis_len[j] + 1 > lis_len[i]:
                lis_len[i] = lis_len[j] + 1
                prev[i] = j
            else:
                pass
    final_pos = max(range(len(l)), key=lambda x: lis_len[x]) 
    output = []
    while final_pos != -1:
        output.append(l[final_pos])
        final_pos = prev[final_pos]
    return output[::-1]

只要將數列倒序輸入這個副程式，再反轉回傳結果，即可獲得最長遞減子數列，因此不用額外寫一個專門的副程式。由於題目提供的輸入為純文字檔，所以在執行程式前要記得將文字轉換為整數，並在輸出答案時再將整數轉換回文字以利後續排版。若省略型態轉換的步驟，則程式會依照字母排序，結果自然會與預期有所出入。

1
2
3
4
5
def rosalind_lgis(path):
    with open(path, "r") as f:
        nums = [ int(n) for n in f.readlines()[1].strip().split() ]
        print(" ".join([ str(n) for n in lgis(nums) ]))
        print(" ".join([ str(n) for n in lgis(nums[::-1])[::-1] ]))

討論

這道題的解法是動態規劃的案例。對我來說，最困難的不是熟悉動態規劃的概念，而是判斷適用動態規劃的情境，還有選擇合適的資訊存儲和判斷式。雖然與資訊科學無關（但跟生物有點關係），練習題目時，一直想起高中生物老師耳提面命的話，

蓋棺才能論定

每個人多少有想追求的目標與想累積的經歷，當中免不了要面對兩難的人生抉擇：接受新的職位，還是待下去？繼續打拼，還是止損停利？與他人合作，還是獨行江湖？即使擁有預知能力，可以判斷當下決策帶來的影響，也很難保證長期而言，會帶來最高的效益。

於是，我們也只能依照當下的條件，選擇最契合個人方向與渴望的行動，直到人生告一段落再也沒有任何可能性的時候，才有機會判斷哪些舉動與決策對於當初訂立的目標最有幫助。只是知道了又如何？雖然其他人可以從這些經驗得到啟發，可是當事人也沒機會用到自己一生濃縮的最佳解。

何況每個階段的追求也可能大相逕庭，處心積慮規劃了老半天，才在某個陽光明媚的早晨發現，自己渴望 B 而不是 A 的情況應該也不少見。

因此，當初生物老師在聯考後第一堂課不檢討考卷，而花了老半天講些莫名其妙的話的時候，可能不是在勉勵我們「別因挫折而氣餒，能否達成目標還說不定呢！」（對，我們班上很多人沒考好）。他也許是期望我們就用力活著，去體驗挫敗、毫無意義、徒勞無功之類說不上是好是壞的情感，不用煩惱纏繞腦中的時程、途徑、機會、指標、最佳化等有的沒的東西。