Rosalindを解く - フィボナッチ数列

Rosalindの紹介はこちらから。

www.kimoton.com

本日は自然界でも見られる数列、フィボナッチ数列についての実装を見ていきます。

rosalind.info

生物学的知識のおさらい

参考:相補性 (分子生物学) - Wikipedia

中世で最も才能があったと評価されるイタリアの数学者、レオナルド=フィボナッチは、ウサギの個体群の繁殖に関する問題を考案しました。

  1. 個体数は最初の月に生まれたばかりのウサギのペアが存在する時点からカウントする。
  2. ウサギは1ヶ月後に生殖年齢に達するとする。
  3. 生殖年齢のすべてのウサギは、生殖年齢の別のウサギと交尾する。
  4. 2匹のウサギが交尾してからちょうど1か月後に、オス1匹とメス1匹のウサギが生まれる。
  5. ウサギが死んだり、繁殖を停止したりすることは考慮しない。

このとき、1年間に何組のウサギが残るか。 上記の問題の答えとしては、1年後、ウサギの個体数は144ペアとなります。

ウサギの不死には異議を立てられてしまいますが、捕食者のいない環境での繁殖ではこのモデルは当てはまるとされており、現に花びらの数や花や実に現れる螺旋の数ははフィボナッチ数であることとなることが知られています。

どの月にも、前月に生きていたウサギと新しい子孫が含まれます。すなわち、任意の月における子孫の数が2か月前に生きていたウサギの数に等しくなります。

結果として、 $F_{n} = F_{n-1} + F_{n-2}$($F_1 = F_2 = 1$)という漸化式によって定義されるフィボナッチ数列を得ることになります。 なお、この数列にはフィボナッチの名前が付いていますが、2千年以上前にインドの数学者に知られていたようです。

問題

「$n\leq40$ の $n$、$k\leq5$ の $k$ が与えられる。1ペアから始めて、各世代において繁殖年齢のペアのウサギが $k$ ペアのウサギの子孫を残すとしたとき、$n$ か月後に存在するウサギのペアの総数を求めよ」

Given: Positive integers n≤40 and k≤5.
Return: The total number of rabbit pairs that will be present after n months, if we begin with 1 pair and in each generation, every pair of reproduction-age rabbits produces a litter of k rabbit pairs (instead of only 1 pair).

具体例と方針

具体例として、 $n=5$, $k=3$ の場合のウサギの個体数増加を下記に示します。

f:id:kimoppy126:20201127090952j:plain

任意の月における子孫の数が2か月前に生きていたウサギの数の $k$ 倍に等しくなります。

これは、漸化式で書くと下記のようになります。

$$ F_{n} = F_{n-1} + k \times F_{n-2} $$

この漸化式における $n$ 番目の値を求めればよいわけです。

フィボナッチ数列の算出方法には計算量の異なる様々な方法があります。今回は動的計画法を用いた実装を見ていきます。

解 (Julia)

Juliaでフィボナッチ数列を実装している良い例があったのでこちらを参考にさせて頂きました。

qiita.com

Juliaでは多重ディスパッチを定義することができます。

多重ディスパッチ(英: Multiple dispatch)またはマルチメソッド(英: Multimethods)は、多重定義された関数やメソッドからそこで呼び出されるべき1つの定義を選出し実行する(ディスパッチする)際に、2個以上の複数の引数が関与してどれかひとつを選ぶこと(特殊化)がおこなわれるものである。

ポイントは引数の判定が呼び出し時に行われる点です。
オーバーロード(継承)では、メソッド引数の定義型を基に静的にコンパイル時にメソッドの呼び出しが決定されますが、多重ディスパッチではこれが呼び出し時に行われるわけです。

この特性を活用したフィボナッチ配列の関数を作成します。

0. ディレクトリ構成

ディレクトリ構成は下記とします。

$ tree 
.
├── main.jl # 実行ファイル
├── data 
│   └── rosalind_fib.txt # 入力データ
└── utils
    └── file_handler.jl # 実行ファイルで使用するモジュール

1. ファイルを読むためのutil関数を作る

DeliminatedFiles パッケージを使用します。readdlm関数では、ファイル名、デリミタ、型の順で指定してあげます。

using DelimitedFiles

function read_file(file)
    contents = readdlm(file, ' ', Int)
    return contents
end

2. メイン関数を作成する

動的計画法で求めます。

動的計画法とは、メモ化を活用しながら計算することで処理時間を短縮するアルゴリズムを指します。

動的計画法は、計算機科学の分野において、アルゴリズムの分類の1つである。対象となる問題を複数の部分問題に分割し、部分問題の計算結果を記録しながら解いていく手法を総称してこう呼ぶ。

毎回計算するとフィボナッチ数列の演算は $O(2^N)$ かかってしまいますが、この動的計画法を用いることで計算量は$O(N)$となります。

#!/usr/bin/env julia

include("utils/file_handler.jl")

# 多層ディスパッチで使用するInteger型を引数にとるメソッドを定義
function fib(n::Integer, k::Integer)
    # メモ初期化
    d = Dict(1=>big"1", 2=>big"1")
    fib(n, d, k)
end

# 引数が3つの場合は下記の関数に渡される
function fib(n, d, k)
    # メモ内に値があれば返す
    if haskey(d, n)
        return d[n]
    end
    # n-1, n-2項目の値からn項目の値を求める
    result = fib(n-1, d, k) + k * fib(n-2, d, k)
    # メモ化
    d[n] = result
    return result
end

n, k = read_file("data/rosalind_fib.txt")
println(fib(n, k))

Tips
Juliaではmethodsという便利関数があり、これを使うことで定義されているすべてのメソッドとその引数タイプを表示してくれます。

> methods(fib)
# 2 methods for generic function "fib":
[1] fib(n::Integer) in Main at main.jl:4
[2] fib(n, d) in Main at main.jl:11

解 (Python)

いつも通り、Pythonの実装に関しても載せておきます。

0. ディレクトリ構成

ディレクトリ構成は下記とします。

$ tree 
.
├── main.py # 実行ファイル
└── data 
    └── rosalind_fib.txt # 入力データ

1. ファイルを読むためのutil関数を作る

with句を使ってファイルを読み書きすると安全にファイルを開けるため便利です。

#!/usr/bin/env python

import os
from typing import List

def read_file(file_name: str, dir_name: str="data") -> List:
    file_path = os.path.join(dir_name, file_name)
    with open(file_path, 'r') as rf:
        contents = rf.read()
    return contents

2. メイン関数を作成する

Juliaの場合と大枠同様です。Pythonでは外部ライブラリを使用しない限りは多層ディスパッチを用いることはできないため、 n=1, n=2のときの値は条件分岐で入れてあげる必要があります。

#!/usr/bin/env python

from utils.file_handler import read_file

# メモの初期化
memo = {}
def fib(n, k):
    # メモ内に値があれば返す
    if n in memo:
        return memo[n]
    if n == 1:
        # n=1, n=2 の場合は1を返す
        result = 1
    elif n == 2:
        # n=1, n=2 の場合は1を返す
        result = 1
    else:
        # それ以外の場合は漸化式に従い、n-1, n-2項目の値からn項目の値を求める
        result = fib(n-1, k) + k * fib(n-2, k)
    # メモ化
    memo[n] = result
    return result

n, k = read_file("rosalind_fib.txt").split()
print(fib(5, 3))

実行

$ cat data/rosalind_revc.txt
5 3

$ python main.py
## 19

コード例

下記で公開しているので参考にしてみて下さい。

repl.it

(おまけ)バイナリ公式を使用するver

通常のフィボナッチ数列を求める場合下記のような公式を使うことで計算量を短くすることができます。

$$ \begin{eqnarray*} F_{2n} &=& (2F_{n-1} + F_n) F_n \\ F_{2n+1} &=& F_{n}^2 + F_{n+1}^2 \end{eqnarray*} $$

こちらを使うと、計算量は$O(logN)$となります。

#!/usr/bin/env julia

# 多層ディスパッチで使用するInteger型のメソッドを定義
function fib(n::Integer)
    # メモ初期化
    d = Dict(zero(n)=>big"0", one(n)=>big"1")
    fib(n, d)
end

# 引数が2つの場合は下記の関数に渡される
function fib(n, d)
    # メモ内に値があれば返す
    if haskey(d, n)
        return d[n]
    end
    # バイナリ公式を利用
    m = n ÷ 2
    result = if iseven(n)
        (2 * fib(m - 1, d) + fib(m, d)) * fib(m, d)
    else
        fib(m, d) ^ 2 + fib(m + 1, d) ^ 2
    end
    # メモ化
    d[n] = result
    return result
end

println(fib(20))