何度も繰り返しているが、計算速度を欲する場合には Python ではなく C 言語などを使うべきである。 knight22a.py と同じ内容の処理を行う C 言語版の knight23a.c を用意した。これと setup.py を 同じディレクトリに置いて "pip3 install ." を実行することで、knight2 というモジュール名で C 言語版の knight23a がインストールされる。 knight22b.py に修正を加え、knight23a を使用するようにしたものを knight23b.py としよう。速さの違いを実感されよ。
また、ジョブ分配を 2 手目ではなく 3 手目で行うようにしたものが knight24a.c, setup.py, および knight24b.py である。
上に示した knight21b.py では、1 個の初期座標を 1 個のプロセスが担当していたが、この方法には少しばかり弱点がある。 一つは、演算装置がマス目の数よりも多いような場合に、全ての演算装置を使い切ることができない、という問題である。 たとえば、通常の PC を 10 台ほど連結した計算システムなどは、大学の研究室単位などで安価に運用することができる。 この場合、たとえば PC 1 台あたりの CPU コアが 6 個だとしても、全部で 60 コアとなる。 7 x 7 のチェス盤では、このうち 49 個のコアしか計算に使えないので、残り 11 個が暇になってしまうのである。 もう一つの問題は、計算終了時刻のバラツキが大きくなることである。 初期座標によって、すぐに計算が終わる場所と、多量の計算を要するために長く時間がかかる場所とがある。 そのため、全体としては「最も長く時間がかかる場所」の計算が終わるまで待たなければならないから、計算終了間際には実際に稼働している CPU の数が少なくなってしまう。 最悪のケースでは、「49 番目、つまり最後の初期座標をある CPU が計算開始した直後に、たまたま、他の全ての CPU が演算終了した」ということが起こりえる。 このとき、1 個の CPU だけが稼働する状態が続き、その CPU が仕事を終えるまで他の CPU はただ待ち続ける、という状況が生じる。
これらの問題を解決するために、初期座標ごとにプロセスを割り当てるのではなく、「初手と 2 手目の組み合わせ」ごとにプロセスを割り当てることにしよう。 こうすれば、たとえば 7 x 7 のチェス盤であっても 160 個の (2 手目で「行き止まり」になるものを除く。Python で数えてみるとよい。) 演算装置を使うことができる。 また、前述のような「最悪のケース」であっても、1 個の CPU のみが稼働している時間は最長で「1 通りの『初手と 2 手目の組みあわせ』を計算する時間」であるから、 knight21b.py の「1 通りの初手を計算する時間」に比べてマシである。 このような改修を考える。
knight21 で使っていた runsingle() は 1 通りの初手を計算するメソッドなので、そのままでは使えない。 1 通りの「初手と 2 手目の組み合わせ」を計算するメソッド runsecond() を knight21a.py に追加して knight22a.py を作ろう。 また、knight21a.init() で progress_total を計算する際に、ついでに「初手と 2 手目の組み合わせ」のリストを作ってしまうのがよい。 リストに要素を追加するのは
x1 = []
y1 = []
x2 = []
y2 = []
(中略)
x1.append(pos_x)
y1.append(pos_y)
x2.append(x)
y2.append(y)
といった具合で行うことができる。 このようにして作ったリストを用いて、knight22b.py ではジョブ分配を行うのがよいだろう。
今度は knight14 を MPI 化してみよう。 具体的には、次のような方針で knight21b.py を作ってみることにする。
なお、MPI で Send() や Recv() する際には、データの大きさに少々の注意を要する。 すなわち、numpy で dtype="int32" として MPI の Send() や Recv() で MPI.INT を指定する場合、数値を 32 ビットの整数で表現することになる。 つまり 2 進数で 32 桁までしか扱わないことを暗に前提としているので、32 桁を超える場合、つまり 4,294,967,296 を超えるような場合には、桁があふれてしまう。 (コンピューターに詳しい人は「アレ？符号の問題は？」と思ったかもしれないが、それは本質的な問題ではないのでここでは考えない。) そこで、ナイトの周遊のように、それなりに大きな数値を扱う場合には 64 ビットの整数値を使うのがよい。 これは、numpy で dtype="int64" とし、MPI の Send() や Recv() では MPI.LONG_LONG を指定すればよい。 世の中には 64 ビット整数でも足りないような問題があるのだが、それはここでは触れないことにしよう。
また、OpenMPI の仕様として、Barrier() や Recv() などで待機しているプロセスは、実際には何も有意義な計算をしていないにもかかわらず、CPU を占有するらしい。 我々のスクリプトでは、ランク 0 のプロセスで計算結果が返ってくるのを待つ部分が、これに該当する。そこで、Recv() の前に
while not mpi_comm.iprobe(source = self.rank):
time.sleep(1)
mpi_comm.Recv((buf_recv, 2, MPI.LONG_LONG), source = self.rank)
のような処理を入れよう。 iprobe() は、メッセージが送られてきているかどうかを確認するメソッドであって、何も送られて来ていなければ False を返す。 すなわち、メッセージの有無を 1 秒ごとに確認して sleep() しながら待つことにより、CPU を占有しないようにするわけである。 ただし、事情がよくわからないのだが、私の環境では 1 回目の iprobe() はメッセージの有無にかかわらず False を返すようなので、ここでは 2 連続で iprobe() を呼んでいる。 なお、sleep() を使うために time モジュールをインポートしておく必要がある。
ナイトの周遊問題に戻ろう。 純粋な Python スクリプトとしては knight14a.py と knight 14b.py が、C 言語との組み合わせでは knight15a.c と knight15b.py が、現時点での最新版である。 knight16 から knight 19 までには、本質的な部分での追加が行われていないので、knight14 と knight15 を元に、もう少し修正を加えよう。 まず、46 日目に紹介した numpy を使うよう、knight14a.py と knight14b.py を修正しよう。 numpy で二次元配列を作るには、たとえば
test = np.array([[0]*8]*3)
などとすればよい。 関連する部分を修正したものを knight20a.py, knight20b.py としよう。
筆者の環境では「python knight20b.py 5 6 6」の実行に実時間で 57 分、「python knight14b.py 5 6 6」は 31 分を要した。 すなわち、numpy を使わない knight14 の方が高速であった。 numpy が威力を発揮するのは、大きな行列や配列に対して複雑な演算をする場合であって、 今回の我々のように単にデータを格納するだけの配列ならば、numpy よりも通常のリストを使った方が良いものと思われる。
53-1.py ではデータ交換の前に Barrier() を呼んでいた。 これを省略すると、たとえば「上側からデータが来ると思って Recv() したら、先に右側からデータが来た」というような場合に問題が生じる。 Barrier() を呼べば、上下方向の通信が全部終わるまで左右方向の通信が始まらないため、そのようなエラーは生じない。 しかし頻回に Barrier() を呼ぶと、全体の処理効率が落ちてしまう。 従って、このような場合は Barrier() を呼ぶのではなく、Recv() で送信元を指定してしまう方がよい。 これは Recv() や Sendrecv() に "source =" を指定することで実現できる。
上に示した 54-1.py の例では、 プロセス間のデータやりとりは Transmit() 関数、計算終了後の集計 (データつなぎあわせ) は Integration() 関数で行うことにした。 細胞膜を介するイオンの出入りは Evolution() 関数である。 ここでは、10 x 10 の小さな系において、座標 (source_x, source_y) の膜電位を最初の 50 ステップは 100 に固定した。 1500 ステップ (150 x 10; resolution の値に注意) 計算した後の x 方向の膜電位の分布を 54-1.png に、3 次元プロットを 54-2.png に、 また初期状態で興奮している部分の隣にある心筋細胞の膜電位の経時的変化を 54-3.png に出力した。 この膜電位のグラフがそれらしい形になるよう、パラメーターを調整した。 一見、もっともらしいグラフが描けたようにみえるのだが、計算するステップ数を 10000 にすると、54-1.png から活動電位が消失する。 これは、既に興奮している細胞から流入するイオンの量に比べて、周囲に流出するイオンが多すぎるために、興奮が伝播する過程で徐々に活動電位が小さくなっているからである。 むろん、現実にはそのようなことは起こらない。我々のシミュレーションが悪いのである。 パラメーターをいじっても、なかなか、現実の興奮の伝播を再現することができない。 実際の興奮は、もっと複雑な機構でチャネルが開閉しているので、このような簡易なモデルで再現することは困難なのであろう。 もっと複雑なモデルを作っても良いのだが、Python の練習という意味では、これで十分である。
50 日目から 53 日目までにみてきた部分が、心筋の興奮伝播を二次元に拡張する際にキモとなる部分である。 これをふまえて、49-1.py を二次元化した 54-1.py を作ろう。
上下方向のデータ交換ができたので、今度は左方向とのデータ交換を行う。 基本的には上下方向と同様なのだが、データが連続していないので、for ループを使う必要がある。 53-1.py では、Transmit() の中に以下のように左側へデータを渡す部分を作ろう。
info.comm.Barrier()
# 左側にデータを渡す
target = info.GX * info.grid_y + info.grid_x - 1 # 送信先のランク番号
src = (info.size_x + 2) + 1 # 送信データの先頭位置
dest = (info.size_x + 2) + 1 + info.size_x # 受信データを格納する位置
buf_send = np.zeros(info.size_y, dtype="i") # 送信用バッファ
buf_recv = np.zeros(info.size_y, dtype="i") # 受信用バッファ
for loop_y in range(0, info.size_y):
buf_send[loop_y] = cells[src + (info.size_x + 2) * loop_y]
if info.grid_x == 0:
info.comm.Recv((buf_recv, info.size_y, MPI.INT))
for loop_y in range(0, info.size_y):
cells[dest + (info.size_x + 2) * loop_y] = buf_recv[loop_y]
elif info.grid_x == info.GX - 1:
info.comm.Send((buf_send, info.size_y, MPI.INT), target)
else:
info.comm.Sendrecv((buf_send, info.size_y, MPI.INT), target, recvbuf=(buf_recv, info.size_y, MPI.INT))
for loop_y in range(0, info.size_y):
cells[dest + (info.size_x + 2) * loop_y] = buf_recv[loop_y]
(中略)
for loop_y in range(1, info.size_y + 1):
pos = (info.size_x + 2) * loop_y + info.size_x
cells[pos] += cells[pos + 1] * 10
動作が確認できたら、同様に右側とデータ交換も行う部分も作り、53-2.py としよう。
今日は Transmit() 関数を実装しよう。 上下左右の四方向と順次データ交換を行う。まずは上方向とのデータ交換のみを実装して 52-1.py としよう。
info.comm.Barrier()
# 上側にデータを渡す
target = info.GX * (info.grid_y - 1) + info.grid_x # 送信先のランク番号
src = (info.size_x + 2) + 1 # 送信データの位置
dest = (info.size_x + 2) * (info.size_y + 1) + 1 # 受信データを格納する位置
buf_send = cells[src:src + info.size_x] # 送信用バッファ
buf_recv = np.zeros(info.size_x, dtype="i") # 受信用バッファ
if info.grid_y == 0:
info.comm.Recv((buf_recv, info.size_x, MPI.INT))
cells[dest:dest + info.size_x] = buf_recv
elif info.grid_y == info.GY - 1:
info.comm.Send((buf_send, info.size_x, MPI.INT), target)
else:
info.comm.Sendrecv((buf_send, info.size_x, MPI.INT), target, recvbuf=(buf_recv, info.size_x, MPI.INT))
cells[dest:dest + info.size_x] = buf_recv
info.comm.Barrier()
for loop_x in range(1, info.size_x + 1):
pos = (info.size_x + 2) + loop_x
cells[pos] += cells[pos - (info.size_x + 2)] * 10
ここでは Recv(), Send(), Sendrecv() と、大文字で始まるメソッドを使っている。 これは、複数のデータを一度に送受信するには recv(), send(), sendrecv() は使いにくいからである。 Recv() メソッドの第一引数には、受信したデータを入れるバッファ、データの数、各々のデータの型、をまとめたタプルを渡している。 MPI.INT というのは整数型の意味である。 また、np.zeros() に dtype="i" という引数が増えた。これは整数型として配列を作る、という意味である。 これまでは省略していたが、Send() や Recv() などと組み合わせる際には省略するとデータを正しく変換できないようである。 小数を扱いたいなら、dtype="float64" として、Recv などでは MPI.DOUBLE を指定すればよい。
この 52-1.py を -np 9 などで実行してみよう。意図した通りに表示されることを確認されたい。 うまくいったら、同様にして下方向についても実装し 52-2.py としよう。
二次元心筋塊について、キチンと意図した通りの通信ができるか確認するスクリプトを書いてみよう。 プロセス間で情報をやりとりする部分が、一次元の場合に比べて y 方向も通信する分だけ、ややこしくなる。 また gather した後の情報統合も、少しばかり複雑になる。 まずは gather 部分の動作確認スクリプト 51-1.py を次のように作ろう。 なお MPIinfo クラスの定義部分は重複になるので省略しているが、50-2.py と同様に作ればよい。 また Transmit() 関数ではプロセス間のデータやりとりを行う予定であるが、未実装である。 pass というのは、「文法的に何かを書かなければいけないが、実際には何もしない」という場合に使うコマンドである。
def Transmit(cells):
pass
def Integration():
all_cells = np.zeros((info.size_x * info.GX) * (info.size_y * info.GY))
for proc in range(0, info.size):
for loop_y in range(0, info.size_y):
dest = info.size_x * info.GX * (info.size_y * grid_y[proc] + loop_y) + info.size_x * grid_x[proc]
src = (info.size_x + 2) * loop_y + 1
all_cells[dest:dest + info.size_x] = data[proc][src:src + info.size_x]
View(all_cells)
def View(cells):
for loop_y in range(0, info.size_y * info.GY):
for loop_x in range(0, info.size_x * info.GX):
print("{0:3d}".format(int(cells[info.size_x * info.GX * loop_y + loop_x])), end="")
print("")
info = MPIinfo(20, 20)
my_cells = np.full((info.size_x + 2) * (info.size_y + 2), info.rank)
Transmit(my_cells)
data = info.comm.gather(my_cells, 0)
grid_x = info.comm.gather(info.grid_x, 0)
grid_y = info.comm.gather(info.grid_y, 0)
if info.rank == 0:
Integration()
my_cells には自分の担当領域の膜電位情報を格納する。all_cells は gather で集めた全プロセス分のデータを格納する。 all_cells は計算条件によっては非常に大きなサイズになるので、メモリを節約するためには、ランク 0 のみで実行される Integration() 関数の中で定義するべきである。 my_cells 自体は一次元の配列だが、データとしては二次元分の情報を格納している。 そのため、座標 (x, y) のデータを ((info.size_x + 2) * y + x) 番目に格納する、というように対応させている。 all_cells についても似たような対応関係を使っているが、info.GX なども絡むので少々、複雑にみえる。
これで mpirun -np 4 python 51-1.py などと実行して、表示を確かめてみよう。 -np に渡す値をイロイロといじって、領域の分割方法が変わる様子を確認するとよい。 なお、テスト目的で -np に CPU コア数より多い値を渡したい場合は --oversubscribe をつけて mpirun -np 10 --oversubscribe python 51-1.py などとするとよい。
テキスト (一週間でなれる! スパコンプログラマ) の Day 5 に入ろう。 前回までに、一次元の心筋における興奮の伝導をシミュレートした。 しかし一次元では、単に興奮が伝播していくだけなので、あまり面白くない。 そこで話を二次元に広げよう。 49-1.py で行ったのと同様の計算を、平面的に広がる心筋塊について、行ってみるのである。
二次元の場合、領域の分割がいささか面倒である。 たとえば 4 プロセスの並列計算であれば、全体を 2 x 2 の 4 領域に分割すれこともできるし、1 x 4 に分割することもできる。 どちらでも計算自体はできるのだが、プロセス間の情報伝達をなるべく少なくした方が、計算の効率は上がるであろう。 その意味で、領域同士が接している部分の長さは、できるだけ短い方がよい。 そのような分割方法は、MPI.Compute_dims() メソッドで得ることができる。 50-1.py を次のように作ろう。
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
division = MPI.Compute_dims(size, 2)
if rank == 0:
print(division)
Compute_dims() の第一引数にはプロセス数を、第二引数には次元を、それぞれ与える。 二次元の場合であれば、division には縦と横の最適な分割数が与えられる。 mpirun の -np に渡す値を変えて、50-1.py を実行してみるとよい。
たとえばプロセス数が 6 の場合、2 x 3 の領域に分割することになる。 この場合、たとえばランク 3 のプロセスは、いったい、どのプロセスと通信すればよいのか？などの情報がわかりにくいので、これをまとめて保存しておくと便利であろう。 そうした情報保存用のクラスを作ってしまおう。50-2.py は次のように作る。
from mpi4py import MPI
class MPIinfo():
rank = 0 # ランク番号
procs = 0 # 総プロセス数
GX = 0 # 横方向の分割数
GY = 0 # 縦方向の分割数
grid_x = 0 # 自分が担当する位置 (横方向)
grid_y = 0 # 自分が担当する位置 (縦方向)
size_x = 0 # 自分が担当する領域の横方向の大きさ (「のりしろ」部分を含まない)
size_y = 0 # 自分が担当する領域の縦方向の大きさ (「のりしろ」部分を含まない)
def __init__(self, total_size_x, total_size_y):
self.comm = MPI.COMM_WORLD
self.size = comm.Get_size()
self.rank = comm.Get_rank()
division = MPI.Compute_dims(size, 2)
self.GX = division[0]
self.GY = division[1]
self.grid_x = rank % GX
self.grid_y = rank // GX
self.size_x = total_size_x // GX
self.size_y = total_size_y // GY
領域分割による並列化を行うには、具体的にどのような情報をプロセス間でやりとりすればよいだろうか？ Evolution() 関数では、自分が計算を担当する領域の情報だけでなく、その両隣の細胞の情報も必要である。 そこで
mpi_comm = MPI.COMM_WORLD
mpi_size = mpi_comm.Get_size()
mpi_rank = mpi_comm.Get_rank()
my_cells = n_cells // mpi_size + 2
total_cells = (my_cells - 2) * mpi_size + 2
potential_old = np.full(mycells, E_init)
potential_new = np.full(mycells, E_init)
のようにしよう。各プロセスの担当する領域の長さは n_cells // mpi_size であるが、それよりも 2 長く potential_new を作っている。 これは potential_new[0] と potential_new[my_cells -1] に、自分の担当外領域の膜電位データを格納するためである。 この部分の値は、自分で計算するのではなく、sendrecv() などを使って隣のプロセスから受け取ることにする。 具体的には、Evolution() を呼ぶ前に
mpi_comm.Barrier()
# 左側にデータを渡す
if mpi_rank == 0:
potential_new[my_cells - 1] = mpi_comm.recv()
elif mpi_rank == mpi_size - 1:
mpi_comm.send(potential_new[1], mpi_rank - 1)
else:
potential_new[my_cells - 1] = mpi_comm.sendrecv(potential_new[1], mpi_rank - 1)
# 右側にデータを渡す
if mpi_rank == 0:
mpi_comm.send(potential_new[my_cells -2], mpi_rank + 1)
elif mpi_rank == mpi_size - 1:
potential_new[0] = mpi_comm.recv()
else:
potential_new[0] = mpi_comm.sendrecv(potential_new[my_cells - 2], mpi_rank + 1)
Evolution(loop)
と、すればよかろう。 このような方針で、47-1.py の並列版である 49-1.py を作ろう。
47-1.py の並列化を考えよう。 たとえば全体を 5 個の領域に分割し、各々の領域を 1 個のプロセスが担当する、という状況を考えよう。 興奮伝導の計算では、自分自身の状態と、隣の細胞の膜電位だけで、自分の膜電位の変化量を計算できる。 つまり、各領域について一番端の心筋細胞の膜電位変化を計算するには、隣の領域の一番端の細胞の膜電位情報が必要になる。 この場合、自明並列は使えず、こまめに他のプロセスと情報交換する必要がある。 具体的には、時間を 1 ステップ進めるごとに、隣の領域の一番端の膜電位情報を担当プロセスから受け取る必要がある。 このプロセス間の情報交換は、mpi4py の send(), recv(), sendrecv() を使って実現できる。 次のように 48-1.py を作ろう。
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
data = -1
if rank == 0:
comm.send(rank, 1)
elif rank == size - 1:
data = comm.recv()
else:
data = comm.sendrecv(rank, rank + 1)
for loop in range(0, size):
comm.Barrier()
if loop == rank:
print("{0:d}: {1:d}".format(rank, data), flush = True)
send() メソッドは、指定されたランクのプロセスに、指定されたデータを送信する。 send(rank, 1) であれば、ランク 1 のプロセスに対して rank をデータとして送信している。 recv() メソッドは、データを受信し、その中身を返す。 data = comm.recv() とすることで、送られてきたデータが data に格納されるわけである。 sendrecv() は、send() と recv() の作業を両方とも行う。
print() する部分に for ループを入れているのは、ランク順に表示するための工夫である。 print() 関数に flush = True を加えているが、これは「直ちに表示せよ」というオプションである。 これを入れない場合、print() を実行した際にすぐには画面に表示されず、忙しい処理が終わってからまとめて表示されることがある。 この場合、異なるプロセスから出力される内容については順序が保障されない。実際に、このオプションをナシにして何度かスクリプトを実行し、挙動を確かめてみるとよい。 なお、Cygwin のコマンドラインでは上下の矢印キーを使うことで、過去に入力したコマンドを参照できる。
計算が終了した時点では、各々のプロセスが、それぞれの担当領域の膜電位情報を持っている。 これを処理する方法はいくつか考えられるが、ここではランク 0 のプロセスに情報を集めることにしよう。 これは gather() で行うことができる。次のように 48-2.py を作ろう。
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
data = rank**2
result = comm.gather(data, 0)
if rank == 0:
print(result)
gather() メソッドも、send() と同様に、指定されたランクのプロセスに、指定されたデータを送信する。 しかし send() とは異なり recv() する必要がなく、gather() の返り値は配列になっている。 48-2.py ではまとめて print() しているが、例えばランク 3 から送られた値だけをみたいなら print(result[3]) とすればよい。
前回描いた膜電位はあまりに粗雑であるので、もう少し細かなシミュレーションを行うことにして 47-1.py を作ろう。 各パラメーターを次のように定める。なお、これらの値は最終的な活動電位をそれらしくするために設定したので、実際の細胞とは異なる。
resolution = 1 # 時間分解能: これを大きくすると 1 ステップあたりの変化量が小さくなる
n_cells = 100 # 計算する系に含まれている細胞数 (長さ)
time_max = 1000 * resolution # シミュレートする時間ステップ数
E_init = -85 # 膜電位の初期値
D_diffusion = 0.3 / resolution # 細胞間の膜電位伝播についての拡散係数
E_Na = 45 # ナトリウムの平衡電位
D_Na_static = 0.01 / resolution # 静止状態におけるナトリウムの拡散係数
D_Na_excited = 0.5 / resolution # 興奮状態におけるナトリウムの拡散係数
threashold = -60 # 興奮の閾値
duration_Na = 10 * resolution # 興奮した後、ナトリウムチャネルが開いている時間
E_Ca = 500 # カルシウムの平衡電位
D_Ca_static = 0.001 / resolution # 静止状態におけるカルシウムの拡散係数
D_Ca_excited = 0.1 / resolution # 興奮状態におけるカルシウムの拡散係数
duration_Ca = 300 * resolution # 興奮した後、カルシウムチャネルが開いている時間
E_K = -90 # カリウムの平衡電位
D_K_static = 0.5 / resolution # 静止状態におけるカリウムの拡散係数
D_K_excited = 0.55 / resolution # 興奮した後のカリウムの拡散係数
delay_K = 200 * resolution # 興奮してからカリウムチャネルが開くまでの待機時間
duration_K = 400 * resolution # 興奮した後に開いたカリウムチャネルが閉じるまでの時間 (興奮時を基準とする)
各細胞の膜電位変化は、次のように、細胞膜を介したイオンの出入り (拡散) を考慮するようにする。
for loop in range(1, n_cells + 1):
# 興奮する
if excited[loop] == 0 and potential_old[loop] > threashold:
excited[loop] = step
# 隣の心筋からの拡散
diff_left = potential_old[loop - 1] - potential_old[loop]
diff_right = potential_old[loop + 1] - potential_old[loop]
change_diffusion = (diff_left + diff_right) * D_diffusion
# ナトリウム
diff_Na = E_Na - potential_old[loop]
if excited[loop] == 0 or step - excited[loop] > duration_Na:
change_Na = diff_Na * D_Na_static
else:
change_Na = diff_Na * D_Na_excited
# カルシウム
diff_Ca = E_Ca - potential_old[loop]
if excited[loop] == 0 or step - excited[loop] > duration_Ca:
change_Ca = diff_Ca * D_Ca_static
else:
change_Ca = diff_Ca * D_Ca_excited
# カリウム
diff_K = E_K - potential_old[loop]
if excited[loop] == 0 or step - excited[loop] < delay_K or step - excited[loop] > duration_K:
change_K = diff_K * D_K_static
else:
change_K = diff_K * D_K_excited
potential_new[loop] = potential_old[loop] + change_diffusion + change_Na + change_Ca + change_K
これで、たとえば 50 番目の心筋細胞の膜電位の掲示変化を描くとそれらしい膜電位の図ができる。 また、resolution = 2, time_max = 250 * resolution として、計算終了時の各細胞の膜電位をプロットすると興奮が伝導している図ができる。 この図では、100 番目の心筋細胞は未だ興奮しておらず、90 番目あたりが興奮の最先端である。60 番目はプラトーであり、0 番目も未だ再分極はしていない。 time_max をもう少し大きくすれば、再分極している様子も描かれるであろう。
なお、実際のイオンチャンネルは、このように電位や時間に応じてカチッと開いたり閉じたりするものではない。 興奮は all or none の法則に従うが、チャンネルは確率的に開閉するのである。 従って、この膜電位のシミュレーションはかなりいい加減なものなのであるが、ここではあまり細かいことは気にせず、雰囲気だけ合っていればよしとしよう。
現時点での計算結果をグラフにしてみよう。 グラフの描画には matplotlib モジュールが便利である。Cygwin で「pip3 install matplotlib」を実行してインストールすればよいのだが、 筆者の環境では先に Cygwin に cmake, gcc-g++, ninja, および make をインストールしておく必要があった。 なお、グラフを描くツールとしては gnuplot も優秀であるし、筆者自身は gnuplot 派である。 しかしここでは、python のスクリプト内で描画を完結させるため、matplotlib を使うことにしよう。
matplotlib を使うために、スクリプト内で「from matplotlib import pyplot as plt」としてモジュールをインポートする。 最後の「as plot」というのは、以下「plt」という名前で pyplot クラスを呼び出しますよ、という意味である。省略して pyplot の名前のままで呼び出すことにしても、むろん、構わない。 最低限のプロットを行うためには、x 座標のデータと y 座標のデータをリストなどで作っておく必要がある。 たとえば、10 番目の心筋細胞の膜電位の経時的変化を描くのであれば、横軸は時刻なので「t = [i for i in range(0, time_max)]」などとすればよい。 縦軸は膜電位であるから、メインの for ループの中で「results[loop] = potential_new[10]」などとしておけばよかろう。 そして計算が終わったら、最後に
plt.plot(t, results)
plt.savefig("output.png")
とすることで、t と results の対応関係を PNG ファイルとして保存することができる。たとえば 46-1.png のようになる。
これは、心筋の膜電位としては、かなり変である。というのも、再分極の過程を実装していないので、やたらと高い膜電位が持続してしまっているのである。 再分極は、次のように表現することにしよう。 「興奮してから時間が 20 ステップ経過した後は、隣の心筋細胞との電位差に伴う膜電位変化は無視できるものとし、一方で細胞外への陽イオン流出による膜電位変化が起こるものとする。 この時、細胞外の膜電位は -0.1 として扱い、拡散係数 D は 0.2 とする。これは時間 50 ステップの間 (つまり興奮から 70 ステップ経過するまで) 継続する。」 これで計算してみると、我々が知っている心筋の活動電位とは異なるおかしなグラフが描かれるのであるが、とりあえずは気にせずに先に進むことにする。
ついでに、numpy モジュールを使うことにしよう。 これはベクトル計算を行うためのモジュールであって、リストの代わりに使えることがある。 大規模なデータを扱う際に、for ループなどを回してリストを処理するよりも高速に計算できる、らしい。 らしい、というのは、私は今のところ、自分でそれを確認していないからである。 とりあえず numpy を使うだけであれば、「import numpy as np」とした上で、リストの初期化部分を
t = np.arange(time_max)
potential_old = np.zeros(n_cells + 2)
potential_new = np.zeros(n_cells + 2)
excited = np.zeros(n_cells + 2)
results = np.zeros(time_max)
と、するだけである。np.zeros() は中身が全て 0 で、指定された要素数 (ベクトル的にいえば次元) のベクトルを作成する。 np.arange() は、0 から始まり 1 ずつ増えていく、という内容で、指定された要素数のベクトルを作成する。 これを 46-2.py としよう。
本日からはテキスト (一週間でなれる! スパコンプログラマ) の Day 4 である。 前回は自明並列、いわゆる馬鹿パラにより円周率を計算してみたが、本日は非自明並列計算を行う。 これは、プロセス間に情報をやりとりしながら並列計算を行うものであって、いわゆるスーパーコンピューターの本領が発揮される種類の計算である。 テキストでは拡散方程式を使った熱伝導のシミュレーションを例示している。 これは物理系の出身者にとっては馴染みやすい例であろうが、生物・医学系の人にとってはチンプンカンプンであろう。 そこで、ここでは心臓の刺激伝導系を簡略化したシミュレーションを考えよう。
とりあえずは刺激伝導系を一次元の、つまり一列に並んだ心筋細胞から成っているというシンプルなモデルで考える。 各々の心筋細胞は、実際よりも単純化して、次のような振る舞いをすることにしよう。 これらの心筋細胞はイオンチャネルでつながっているため、隣の細胞の膜電位が上昇すると、イオンの流れによって自分自身の膜電位も上がる。 そして膜電位が閾値を超えると、細胞外からイオンが流入することでさらに膜電位が上がる。これが興奮である。 興奮してから一定時間が経過すると、細胞外からのイオンの流入は途絶える。 この条件で、一番端の細胞 (洞房結節であろうか？) が初期状態で興奮しているとしたとき、どのように興奮が伝播していくかをシミュレートしてみよう。 この多数の心筋細胞を、いくつかの領域に分割し、各プロセスが一つの領域を担当して膜電位の変化を計算する、というのが領域分割による並列化である。 ある心筋細胞の膜電位の変化は、自分自身と、両隣にいる心筋細胞との膜電位の差によって決まる。 従って、各領域の端部にある心筋細胞の膜電位変化を計算するには、隣の領域の一番端にある心筋細胞の膜電位の情報が必要となる。 そのため、他のプロセスと情報をやりとりする必要が生じるのである。
まずは並列化を考えずに、単一プロセスでこの問題を解くことを考えよう。 各細胞の膜電位をリストとして保存しておき、微小な時間が経過した後の膜電位を次に示す 45-1.py のように計算することにしよう。興奮については未実装である。
n_cells = 10000
time_max = 100
D = 0.1
potential_old = [0.0 for i in range(0, n_cells + 2)]
potential_new = [0.0 for i in range(0, n_cells + 2)]
potential_new[0] = 1.0
potential_old[0] = 1.0
def Evolution():
global potential_old, potential_newv
potential_old, potential_new = potential_new, potential_old
for loop in range(1, n_cells + 1):
diff_left = potential_old[loop - 1] - potential_old[loop]
diff_right = potential_old[loop + 1] - potential_old[loop]
potential_new[loop] = potential_old[loop] + (diff_left + diff_right) * D
for loop in range(0, time_max):
print("{0:.2f} - {1:.2f} - {2:.2f} - {3:.2f} - {4:.2f}".format(potential_new[0], potential_new[1], potential_new[10], potential_new[11], potential_new[50]))
Evolution()
n_cells は細胞数を、time_max は何ステップの時間経過を計算するかを、それぞれ表す設定値である。 ある細胞の膜電位が、隣の細胞の膜電位よりも dE だけ低かったとしよう。このとき、時間を 1 ステップだけ進めると、 その細胞には隣の細胞からイオンが流入し、膜電位が少し高くなるであろう。 この膜電位の上昇量は、dE よりも小さく、また、だいたい dE に比例するであろう。 そこで、1 より小さい定数 D を用いて、膜電位の上昇量を D * dE と表現することができるであろう。 この D は、物理学の世界では拡散係数と呼ばれている。 膜電位は potential_old と potential_new の二つのリストに格納されている。 このリストは細胞数 n_cells よりも 2 だけ長くなっている。 これは、計算の便宜上、リストの両端部分については値を固定しておくためである。 具体的には、リストの左端は常に 1, 右端は常に 0 としている。このように、計算の前提条件として値が決められていることを「境界条件」と呼ぶ。
関数 Evolution() では、potential_new の値を potential_old に移した上で、potential_old の値から potential_new の新しい値を計算している。 各々の心筋細胞について、自分自身と両隣との間における膜電位の差に D を乗じた分だけ、膜電位を変化させているのである。 そして、for ループによって、Evolution() を 100 回、行っている。 途中経過を全部出力するとみにくいので、代表的な 5 点について、膜電位の変化を表示することにした。
この 45-1.py に、心筋細胞の興奮を実装しよう。 具体的には「ひとたび potential_old が 0.3 以上になると、その後 20 ステップの間、膜電位が 1 になる。その後は、二度と興奮しない。」ということにしよう。 これを 45-2.py とする。
43-1.py では、各々のプロセスが独自に円周率を推定した。これを集計するには、次のように allreduce() メソッドを使う。 Estimate() 関数をはじめとする前半部分は 43-1.py と同じなので省略する。
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
pi = Estimate(total)
#print("{0:d}: {1:f}".format(rank, pi))
comm.Barrier()
mean_pi = comm.allreduce(pi, MPI.SUM) / size
sum_pi2 = comm.allreduce(pi**2, MPI.SUM)
sigma2 = sum_pi2 / (size - 1) - mean_pi**2 * size / (size - 1)
sigma = math.sqrt(sigma2)
if rank == 0:
print(mean_pi, sigma)
allreduce() では、第一引数に与えられた変数について、第二引数に与えられた処理を全プロセスにわたって行う。 第二引数に MPI.SUM を渡せば、全プロセスにおける総和が計算されるわけである。 途中で comm.Barrier() メソッドが呼ばれているが、これは以前のマルチプロセスやマルチスレッドで使った join() と同様に、全てのプロセスが待ち合わせをする、というメソッドである。 最後の表示部分では、if 文を使わずに print してしまうと、全てのプロセスが同じ内容を表示することになってしまい、無駄である。 そのため、rank 0 のプロセスが代表して表示を行うことにしている。
テキスト (一週間でなれる! スパコンプログラマ) の Day 3 に入る。 ここでは自明並列、俗にいう「馬鹿パラ」を扱う。 馬鹿パラというのは「馬鹿でもできる並列 (パラレル) 計算」のことであるらしく、 複数の CPU が互いに独立に、少しパラメーターが違うだけの同じような計算を行う並列計算をいう。 プログラミングとしては単純なのだが、処理する問題の内容によっては、非常に強力である。 我々が以前、ナイトの周遊問題についてマルチプロセスで解を探索したのも、この自明並列の手法である。 我々がすでに扱ったマルチプロセスの並列処理と、MPI を使ったいわゆるスパコンプログラミングの間に、本質的な違いは存在しない。 ただ、MPI の場合、異なるコンピューターを組み合わせた並列計算が可能である、という特徴がある。 つまり、自宅のデスクトップと、持ち運び用のノートパソコンと、職場にある作業用パソコンとで協力して並列計算する、などということが MPI を使えば可能である。
ここではモンテカルロ法による円周率の推定を行ってみよう。 モンテカルロ法というのは、コンピューターで乱数 (ランダムな数) を生成し、それに応じてシミュレーションを多数回行い、統計的に何かを推定する、という手法である。 たとえば円周率について考える。 幾何学によれば、円とは「中心から等しい距離にある点の集合」であって、半径を r, 面積を S とすると、円周率 pi は pi = S / (r * r) と表される。 幾何学の教科書では、円周率は円周と直径の比として定義されることが多いかもしれないが、上述のように面積から円周率を定義しても同じことである。 さて、1 x 1 の正方形の中に、ランダムに 1 点を選ぶとする。この時、選ばれる点の位置が偏らないようにしよう。 この時、選ばれた点が、正方形の中心から距離 0.5 未満である確率 (頻度) を p としよう。 p は、この正方形の面積と、半径 0.5 の円の面積との比に一致するであろう。 正方形の面積は 1、半径 0.5 の円の面積は 0.25 * pi であることから、p = 0.25 pi となる。 そこで、p をモンテカルロ法により推定することができれば、円周率がわかる、というわけである。
まずは並列計算を行わずに、円周率を推定するスクリプト 43-1.py を書いてみよう。たとえば、次のような具合である。 なお、座標 (x, y) と原点との距離は math.sqrt(x**2 + y**2) で計算できる。これは幾何学の分野では三平方の定理と呼ばれている。
import sys
import random
import math
from mpi4py import MPI
def Estimate(n):
count = 0
for loop in range(0, n):
x = random.random() - 0.5
y = random.random() - 0.5
r = math.sqrt(x**2 + y**2)
if r < 0.5:
count += 1
return count / n * 4
if len(sys.argv) > 1:
total = int(sys.argv[1])
else:
total = 100000
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
pi = Estimate(total)
print("pi = {0:f} ({1:d}/{2:d})".format(pi, rank, size)
これを mpirun により並列計算してみるとよい。
今日は、特に作業はない。 テキスト (一週間でなれる! スパコンプログラマ) の Day 2 を読んで、 いわゆるスーパーコンピューターというのがどういうものなのか、おおまかに理解されよ。
昨日までで、Python の基本的な部分に一通り触れたように思う。 特に言及しなかった要素も少なくないが、適宜検索すれば容易に使えるようになるであろう。 ここからは、上級編として 一週間でなれる! スパコンプログラマ をなぞってみよう。 これは、いわゆるスーパーコンピューター向けのプログラミングテクニックを解説する初心者向けの文書である。 この文書は C/C++ 言語を基本として書かれているようだが、本質的な部分は Python であっても変わらないはずである。 以下、この文書を「テキスト」と表現することにする。 本文書ではテキストの要点のみをなぞり、詳細は省く部分も多いので、以下は適宜テキストを参照しながら読まれるとよい。
まずはテキストの Day 1 からみていこう。 スーパーコンピューターといっても、我々が使うパーソナルコンピューターと根本的には同じような機械である。 ただ、多数の演算装置を搭載しているために多数の計算を同時に進行できる上、一個一個の演算装置も高性能であるから、すごく速く計算を行うことができる。 計算を同時に行うようなプログラミングを、並列プログラミング、と呼んでいる。 この並列計算にあたっては、適切なプログラミングを行わないと効率が上がらず、スーパーコンピューターの性能を十分に引き出すことができない。 なお、計算速度を追求するならば Python を使うべきではないことはこれまでに何度も述べてきた。 従って、スーパーコンピューターで Python を使うなら、本当に Python が威力を発揮する部分に限定するべきであり、計算の主たる部分は C 言語などに任せるべきである。 とはいえ、本文書はあくまで入門書なので、スーパーコンピューターでのプログラミングの要領を理解することを主眼とし、Python と C 言語の使い分けは考えないことにしよう。
並列計算を実現する方法はいくつかあるが、ここでは MPI を使うことにする。 MPI というのは並列計算をサポートするライブラリの規格である。 Cygwin のインストーラーから openmpi と libopenmpi-devel をインストールし、さらに「pip3 install mpi4py」を実行し mpi4py モジュールをインストールするのがよかろう。 なお、著者の環境では mpi4py をインストールするためには、 libz-devel, libhwloc-devel, libevent-devel を Cygwin に追加インストールする必要があった。 諸君の環境では、また違うエラーが出てインストールできないかもしれない。このあたりのエラーに対処するには、ある程度の経験や C 言語の知識がないと厳しいであろう。 周囲のコンピューターオタクやインターネット上のお友達に相談して対処されたい。
スクリプトを並列実行するだけなら、特別な処理をスクリプト内に記載する必要はない。たとえば 41-1.py を
print("Hello, World!")
として、Cygwin から「mpirun -np 2 python 41-1.py」と実行してみよう。 mpirun は MPI で計算を実行するコマンドである。-np は、並列計算に使うプロセス数を指定するオプションである。 上述のコマンドを実行すると、「Hello, World!」が 2 回表示されたはずである。 すなわち、2 個のプロセスが生成され、諸君のコンピューターに搭載されている 2 個の CPU (実際には 1 個の CPU の中の 2 個のコアであったり、あるいは 1 個の CPU が一人二役をしているかもしれないが、そこは本質的な問題ではない) が、各々 41-1.py を実行したのである。実用上の意義はないが、一応、これも並列計算の一種である。
実用上は、プロセスごとに手分けして様々な計算を行うことで、全体の作業を速く進めたい、ということが多い。 従って、各々のプロセスには名前、あるいは識別番号のようなものをつけるのが便利である。MPI では、この識別番号を「ランク」と呼んでおり、次に示す 41-2.py のようにして確認できる。
from mpi4py import MPI
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
print("Hello, World! ({0:d}/{0:d})".format(rank, size))
Get_size() は全体のプロセス数を、Get_rank() は自分のランクを返すメソッドである。 諸君は、このランクの値に応じて各プロセスで異なる処理を行うこともできるし、同じ処理を行うこともできる。 諸君が何をしたいかによって、好きにすればよい。
なお、mpirun を実行する際に、-np に渡す数を大きくしすぎるとエラーになる。 大抵、コンピューターに搭載されている演算装置の数、つまり CPU コア数が上限である。 というのも、コンピューターに搭載されている演算装置より多くのプロセスを作って並列計算しても、実際上の処理速度は (ほとんど) 向上しないからである。 この制限を解除する方法もあるが、実用上の意義が乏しいので、ここでは触れない。 なお、諸君のコンピューターの論理プロセッサ数は、Windows のタスクマネージャーで確認することができる。
ここまで Python を使ってナイトの周遊問題の解を探索してきた。 しかしながら、何度も繰り返しているように、処理速度の観点からいって、こうした処理には python よりも C 言語などを使う方がはるかに適している。 Python の利点としてしばしば指摘されることの一つは、複雑な処理を比較的簡単な表現で書ける、という点である。 たとえば、変数 a に 10 を、変数 b に 20 を代入してから、a と b の中身を交換したい場合、C 言語であれば
int a, b, c;
a = 10;
b = 20;
c = a;
a = b;
b = c;
などと詳細に書かねばならないが、python であれば
a = 10
b = 20
a, b = b, a
と、比較的簡潔な表現で済む。これは、python は処理内容が曖昧でブラックボックス化しているという面もあるが、ちょっとしたスクリプトを簡単に書く際には便利な仕様である。
Python を使うことのもう一つの利点は、豊富なモジュールが多くの人によって作成されている点にある。 たとえば Excel ファイルを操作しようと思うなら、openpyxl というモジュールがある。 これをインストールするには、Cygwin で「pip3 install openpyxl」を実行すればよい。 pip3 は、このような python モジュールのインストール等を管理するソフトウェアである。
個々のモジュールの使い方を説明することは、この文書の趣旨を逸脱するので、ここでは行わない。 インターネット上には、これらのモジュールの使い方を説明している文書も多いので、敢えてここで繰り返す必要もなかろう。 試しに、knight17a.py, knight17b.py を修正し、発見した解を Excel ファイルとして出力するものを knight19a.py, knight19b としよう。
knight17a.py, knight17b.py, knight17c.py, knight17d.py を修正して、ファイルをバイナリ形式で読み書きするようにしたものを knight18a.py, knight18b.py, knight18c.py knight18d.py としよう。
knight17c.py および knight17d.py により、我々が発見した解は正しく、かつ重複がないことを確認することができた。 これで不都合はないのだが、コンピューターへの理解を深める意味で、もう少しいじってみよう。 knight17a.py では、解をテキストファイルで作成していたが、これをバイナリファイルで保存してみよう。 今日は、その準備である。
通常のコンピューターでは、テキストファイルにせよ画像ファイルにせよ、根本的には数値として扱っている。 それも、数値は 2 進法で取り扱うのが主流である。 2 進法による数の表現とは、たとえばコインを並べることで数を表現するようなものである。 コインの表を 1, 裏を 0 と決めておけば、たとえば 3 という数は 2 進法で 11 であるから、コインを「表表」と並べることで表される。 同様に、91 は「表裏表表裏表表」で表現できる。 これはコインの表裏でなくても、何か 2 つの状態をとるものであれば何を使って表現することもできる。 現代のコンピューターでは、磁気、つまり磁石の向きであるとか、電気、たとえばコンデンサーが充電されているか放電されているか、とか、電磁気を使うことが多い。 このような、二つの状態をとりえるものが表す情報量を 1 ビット (1 bit) と呼ぶ。 要するに 2 進数の一桁のことを 1 ビットと呼ぶ、と考えて差し支えない。
1 ビットで表されるのは「表か裏か」という限られた情報だけなのであまり多くの物事を表現できない。そこで実用上は、複数ビットをひとまとめにして扱うことで、豊富な情報を表現できるようにしるのが普通である。 このように、ひとまとめにして扱う複数ビットの単位をバイト (byte) と呼ぶ。 歴史的経緯から、現代のコンピューターでは、8 ビットを 1 バイトとして扱うのが普通である。 従って、1 バイトは 8 桁の 2 進数である、と考えてよい。 これは 00000000 から 11111111 までのパターンがあり、10 進数でいえば 0 から 255 である。 コンピューターオタクの人々が 256 や 1024 あるいは 65536 を「キリが良い数値」と考えるのは、それがちょうど 8 ビット、10 ビット、16 ビットの場合のパターン数に相当するからである。
1 バイトは 0 から 255 までの値をとりえるが、これを 16 進数で表すことがある。16 進数を使うと、1 バイトがちょうど 2 桁になり、表記しやすいのである。 普段我々が使う数字は 0 から 9 までしかないので、便宜上、10 から 15 までは A から F までのローマ字を使って表すことにしている。 たとえば 2 進数でいう 193 は 16 進数でいうと C1 になる。この場合は 16 の位が C なので 16 進数であるとすぐにわかるが、たとえば 82 とだけ表現されると、それが 10 進数なのか、あるいは 130 を表す 16 進数なのか、わかりにくい。 そのため、16 進数である場合には頭に 0x をつけて 0x82 などと表現することも多い。
コンピューター上では、文字もあくまで数値として扱われる。どの文字をどの数値で表現するのか、という対応関係を「文字コード」と呼んでいる。 以前は、言語ごとに異なる文字コードを使うのが主流であったが、近年では全ての言語を同一の文字コードで表現する Unicode が主流になりつつある。 Unicode というのは文字コードの規格であって、具体的な文字コードとしては UTF-8 などが使われる。 アラビア数字やローマ字は、ほとんどの文字コードで 1 バイトであるが、日本語の文字は 2 バイト以上であることが多い。
ここからが本題である。Python で 10 という値をファイルに書き込むとき具体的には、どのようなデータを我々は書き込むのだろうか？ これを「10」という文字を表す 2 バイトのデータとして書き込むのか、あるいは 0x0A の 1 バイトのデータとして書き込むのかは、諸君の自由である。 また 2 バイトの数値、たとえば 0x89F3 を書き込むとき、0x89, 0xF3 の順で書き込むのか、あるいは 0xF3, 0x89 の順で書き込むのかも、諸君の自由である。 前者はビッグエンディアン、後者はリトルエンディアンと呼ばれる方式であり、いずれも実際に用いられている。 いずれを使っても大きな優劣はないが、インターネット上の通信などでは、どちらの方式を使うのか予め決めておかないと、送信側と受信側で齟齬が生じる恐れがある。
Python では、ファイルを開くとき単に open("test.txt", "w") とすれば、テキストモードで開かれる。この場合、「10」を書き込むなら「10」という文字が書き込まれる。 これに対し open("test.bin", "wb") と "b" をつけて開けば、バイナリモードで開かれる。たとえば、次のようなスクリプトを実行してみるとよい。
file = open("test.bin", "wb")
n = 100
file.write(n.to_bytes(1, byteorder="little"))
file.close()
n は整数 (int) 型の変数として扱われている。 to_bytes は、整数型の変数をバイト型に変換するメソッドであり、第一引数は「何バイトのデータに変換するのか」を、第二引数は「ビッグエンディアン ("big") とリトルエンディアン ("little") のどちらを使うのか」を指定する。 このスクリプトによって作られた test.bin は、何バイトの大きさであろうか？
このようなファイルは、メモ帳やサクラエディタなどの普通のテキストエディタで開いても正しく読むことができない。 これらは「バイナリエディタ」と呼ばれるソフトウェアで適切に読んだり書いたりできるので、試してみられると良い。 Python でバイナリファイルを読む場合には、例えば次のように行う。
file = open("test.bin", "rb")
n = int.from_bytes(file.read(1), byteorder="little")
print(n)
file.close
次に解の重複がないことを確認する knight17d.py を作ろう。 6 x 6 のチェス盤では 6,637,920 通り、6 x 7 であれば 779,938,932 通りの解が発見されるのだが、これに重複がないことを確認するのは、一見、とてつもない作業である。 というのも、n 通りの解について、全ての組み合わせについて照合を行うならば、必要な検証の回数は n (n - 1) / 2 回となり、非常にたくさんの計算量が必要となるからである。 しかし幸い、解の重複がないことを確認するだけならば、全ての組み合わせを照合する必要はない。 我々のアルゴリズムでは、ランダムに解を発見するわけではなく、明確な順番に従って解を探索している。 すなわち、ある 2 個の解が示された場合、どちらの解を先に発見したのか、言い当てることができる。 これは次のような関数によって実現できよう。
Knight_x = [-2, -2, -1, -1, 1, 1, 2, 2]
Knight_y = [-1, 1, -2, 2, -2, 2, -1, 1]
def GetIndex(prev_x, prev_y, next_x, next_y):
for index in range(0, 8):
x = prev_x + Knight_x[index]
y = prev_y + Knight_y[index]
if x == next_x and y == next_y:
return index
else:
return -1
def Check(orig, prev, next):
if orig < 0:
return True
orig_x = orig % Size_x
orig_y = orig // Size_x
prev_x = prev % Size_x
prev_y = prev // Size_x
next_x = next % Size_x
next_y = next // Size_x
prev_index = GetIndex(orig_x, orig_y, prev_x, prev_y)
next_index = GetIndex(orig_x, orig_y, next_x, next_y)
if prev_index < next_index:
return True
else:
return False
Check() は、「n (任意の数) 番目に orig の位置におり、n + 1 番目に prev に移動した」という解と 「n 番目までは同じように動いていたが、n + 1 番目に next に移動した」という解を比較し、 前者の方が先に発見されたはずならば True を、そうでないならば False を返す関数である。 この関数を使い、solution.X.Y.txt の全ての連続する 2 行について、重複がなく、かつ発見順に並んでいることを確認すればよい。 連続する 2 行の順序が全て保たれているならば、任意の 2 行について重複していないことが担保されるからである。
Cygwin 上で 2 つのファイルの差分を表示するには diff コマンドを使うとよい。 たとえば「diff knight14b.py knight16b.py」を実行すると、knight14b.py と knight16b.py の差分が表示されるので、どこを修正したのかわかりやすい。 好みの問題だが、私は -u オプションをつけて「diff -u knight14b.py knight16b.py」とした方がみやすいと思う。
さて、次に solution.X.Y.txt の内容が全て正しい解であることを確認しよう。 これは
という手順を繰り返すことで実現できる。ファイルを一行だけ読み込むのは「input()」で実行できるので、これを反復して
while True:
line = input()
# 正しい解であるかどうか検証する
とすればよい。input() は 12 日目に扱った、標準入力から文字列を読み取る関数である。 Cygwin から knight17c.py を実行する際に「python knight17c.py < solution.X.Y.txt」とすれば、solution.X.Y.txt の内容が python の標準入力に渡されるので、input() で読むことができる。 さらに、「cat solution.*.*.txt |python knight17c.py」と実行すれば、sokution.*.*.txt に該当する全てのファイルを結合し、それをパイプで python に渡すことになるので、全ての solution.*.*.txt ファイルを一括で処理できる。 ここで * はワイルドカードと呼ばれる記号であり、任意の文字列として解釈され、該当する全てのファイル名として展開される。 また cat はファイルを結合するコマンドであって、引数に渡されたファイルを全て結合して標準出力に出力する。
上述のスクリプトを実行すると、ファイルの末尾でエラーが出てしまう。そこで、単に input() を実行するのではなく
while True:
try:
line = input()
except EOFError:
break
とするのがよい。これにより、input() を試しに実行し、EOFError のエラーが出たならば break する、という挙動が可能になる。 EOFError というのは、input() 関数がファイル末尾に到達した場合に生じるエラーである。
残りの部分は諸君の好きなように作ればよい。 作業をいくつかの単位に分割し、適宜、関数を作るとみやすいであろう。 私は、リスト形式の data から index 番目の移動先を探す関数 Search(data, index) と、 prev から next に移動できるか否かを検証する関数 Check(prev, next) を作った。
さて、我々のスクリプトは非常に高い効率でナイトの周遊問題の解を探索するが、これらの解は、本当に正しいのだろうか？ 我々が考えたアルゴリズムが正しく実装されているならば、このスクリプトは全ての解を網羅的に探すはずであるが、本当に我々は、意図したアルゴリズム通りにスクリプトを書けているのだろうか？ ひょっとすると、何らかのバグが生じていて、解ではないものを解として数えたり、あるいは同じ解を複数回数えたり、していないだろうか？ また、全ての解を網羅することができずに、見落としが発生してはいないだろうか。 次の要領で、計算結果の検証を行ってみよう。
これらを一度に行うのはたいへんなので、今日は knight17a.py と knight17b.py を作るだけにしよう。 ファイルの読み書きについては 13 日目に一度やっているので、忘れた方は復習されるとよい。
solution.txt の書式はいろいろ考えられるが、ここでは一つ一つのマスについて、そのマスを何番目に訪れるか、の数値を書き込むことにしよう。 6 x 6 のチェス盤であれば、一つの解につき 36 個の数値を書き込むことになる。各々の数値はスペースで区切ろう。
knight14b.py のようなマルチプロセスのスクリプトを実行中に、エラーなどに気が付いて Ctrl + C でスクリプトを終了したとする。 この場合、メインのプロセスは終了できるのだが、スクリプト中で生成したプロセスは終了されずに計算を続ける。 その様子はタスクマネージャーで、CPU 使用率が高いままであり、また「プロセス」をみると python3.9.exe などのプロセスが残っていることからも確認できる。 あるいは Cygwin 上で「ps」を実行すると、Cygwin から実行されている全てのプロセスを一覧表示することができる。 この一覧の中から python 関連のもののみを表示するためには「ps |grep python」と実行すればよい。 この「|」はパイプと呼ばれる記号であって、「ps」を実行して出力される内容を「grep python」の入力に渡す、という意味になる。 「grep」は、入力された内容のうち、引数で指定された文字列を含む行のみを表示する、というコマンドである。 すなわち「ps |grep python」と実行するとdi、python という文字列を含むプロセスのみが表示される。
Cygwin 上でプロセスを終了するには「kill」コマンドを使う。 ps では一列目にプロセスID (PID) が表示されるので、終了したいプロセスの PID を調べて kill の引数に渡せばよい。 たとえば PID 4129, 4130, 4131 の 3 個のプロセスを終了したいなら「kill 4129 4130 4131」と実行すればよい。
さらに手間を省くには「ps |grep python |awk '{print $1}' |xargs kill」を実行するという手もある。 awk は文字列を編集するエディタであって、引数に渡されたコマンドを実行する。「'{print $1}'」は「1 列目のみを表示する」というコマンドなので、これで PID のみが表示される。 xargs は、入力された内容を引数としてコマンドを実行する、というコマンドである。 つまり「4129 4130 4131」という内容を「xargs kill」に入力すれば「kill 4129 4130 4131」が実行される、というわけである。
knight14b.py を修正し、Ctrl + C で終了した際に全てのプロセスを終了するようにしよう。 「import signal」した上で、適切な場所に「signal.signal(signal.SIGINT, signal_hander)」を追加する。 これは、SIGINT のシグナル (python 上では signal.SIGINT) が送られてきた場合には signal_hander 関数で処理する、という意味である。 そして signal_handler 関数を新たに作り、生成されているプロセスを全て終了する処理を追加しよう。 「process.start()」で作成したプロセスを終了するには「process.terminate()」を行えばよい。 その後、念のため「process.join()」しておくとよいだろう。 これらの修正を加えたスクリプトを knight16b.py とする。
knight12a.py では、init() の中で progress_total を計算するにあたり、2 手目で手詰まりになることはない、と暗黙のうちに仮定していた。 しかし実際には、この仮定は正しくない。 たとえば 5 x 5 のチェス盤 (座標は (0, 0) から (4, 4) とする) で考えると、初手が (1, 2) の場合に問題が生じる。 2 手目で (2, 0) に移動すると、(0, 0) と (4, 0) が「行き止まり」になるため、手詰まりなのである。 move() 関数では、2 手目で手詰まりになる場合にはグローバル変数 search を触らずに、次の手の探索を行っている。 すなわち、進捗状況の表示に際しては、手詰まりになる 2 手目が分母には含まれているが分子には含まれていないため、進捗割合が過小評価されるのである。
このバグを修正する方法として、簡単なのは move() 関数において、2 手目手詰まりの場合でも search += 1 を実行する、という方法が考えられる。 たとえば
if dead(matrix, x, y):
if count == 2:
search += 1
matrix[x][y] = 0
continue
とすれば、一応、目的は達される。 しかし、これには 2 点、問題がある。 一つは、2 手目手詰まりの場合は、通常の 2 手目に比べてほとんど探索に時間をかけずに search += 1 が実行されてしまうため、進捗状況の評価としてあまり適切でないように思われる。 もう一つは、この if 文はかなり高い頻度で実行されるので、そこに「count == 2」の場合のみ実行されるような処理を挟むことは、無駄な計算を増やすことになり、効率が良くないように思われる。 実際には、この if 文による計算量の増加は全体からみてごく僅かであるから、あまり重要な問題ではない。とはいえ、無駄な分岐はなるべく減らした方がバグを生じさせないためにも良いであろう。
そこで上のリンク先に示した knight14a.py では、init() 関数の中で dead() 関数を呼んで行き止まり判定をすることにより、2 手目手詰まりを除外して progress_total を計算している。 これにより init() 関数が要する計算量が少し増えてしまうが、init() は一回しか実行されない関数なので、計算量が多少増えたところで問題はない。
何度も繰り返すが、計算量の多い処理は、Python で行うべきではない。 我々の例でいえば、knight14b.py の内容は計算量が少なく、しかしプロセスやスレッドの制御など複雑な処理を実行しているので、便利な (ただし遅い) モジュールが揃っている Python に向いた処理である。 一方、knight14a.py の内容は、スレッド制御などの複雑な仕事はしていないが、計算量が多いため、C 言語などに向いている。 こうした場合、knight14a.py の中身を C 言語で記述し、それを Python から呼び出すのが良い。 ここでは Python から C 言語のプログラムを呼び出す方法については記載しないが、 私が書いた C 言語のソースコードと、これを pip3 などでインストールするための setup.py を置いておこう。 これを呼び出す python スクリプトが knight15b.py である。
knight14b.py と knight15b.py は、行っている処理自体は同じである。計算の中心となる部分を Python で書かれた knight14a.py が担うか、C 言語の knight15a.c で行うかの違いである。 筆者の環境では、プロセス数 5 で 6 x 6 のチェス盤を計算した場合、knight14b.py では実時間で 30 分程度かかったが、knight15b.py では 25 秒であった。 C 言語の方が 60 倍、速かったのである。 なお 6 x 7 のチェス盤では、実時間 50 分で 779,938,932 個の解を発見した。
既にお気づきかもしれないが、この knight13b.py と knight12a.py の組み合わせにおける進捗状況の表示には、不正確な部分がある。 チェス盤の一辺の大きさが 5 以下の場合、進捗状況が過小評価されるのである。 たとえば、進捗表示を 1 秒間隔にして、プロセス数 1 で 5 x 6 のチェス盤について探索を行ってみるとよい。 進捗が 93% 程度になったところで、唐突に終了するはずである。これは、処理が途中で終わっているのではなく、進捗状況の表示が間違っているのである。 このように、設計意図とは異なるプログラムの挙動は「バグ (bug)」と呼ばれる。 また、こうしたバグを修正する作業を「デバッグ (debug)」と呼ぶ。 プログラミングにおいては、このデバッグに膨大な労力を要することが多い。従って、はじめから、バグが生じにくいような、わかりやすいスクリプトを書くことが重要である。 knight12a.py や knight13b.py にデバッグを行い、knight14a.py, knight14b.py を作ろう。
knight12a.py は、マス目の数だけプロセスを生成し、並列計算を行うことで高速化している。 しかしながら、コンピューターが同時に実行可能な演算の数より多くのプロセスを生成しても効率は変わらないので、生成するプロセスはもう少し少なくても構わない。 また、常時 CPU を 100% 使い切る状態になるため、他の作業を並行して行いたい場合には、いささか不便である。 そこで、スクリプトを実行する際に、生成するプロセスの数を指定できるようにしよう。 具体的には「python knight13a.py 10 6 6」などとすれば、10 個のプロセスを生成して 6 x 6 のチェス盤について解を探索するようにしよう。 これを実現するにはいくつかの方法が考えられるが、まずは以下のような方針を採用して knight13b.py を作ろう。 knigh12a.py には手を加える必要がない。
この時問題になるのが、Control.run() が新しい初期位置について探索を開始する際、どの初期位置について計算を行うべきか、という問題である。 単純に考えれば、たとえばグローバル変数 index に「次に探索すべき初期位置」を格納しておけばよい、と思われる。 しかし、複数のプロセスがほぼ同時に終了した場合、複数のスレッドがほぼ同時に index を参照したり更新したりすると、思わぬ間違いが生じるかもしれない。 そういう事態は稀であろうが、稀であっても起こりうる事態には備えておくべきである。そこで threading.Lock() を使う。これは、次のように使う。
lock = threading.Lock()
lock.acquire()
# 何かの処理
lock.release()
lock.acquire() は「ロックをかける」メソッドである。既にロックがかかっている場合、それが解除されるまで lock.acquire() は待機する。ロックの解除は lock.release() で行う。 つまり、あるスレッドが「# 何かの処理」を実行している間に別のスレッドも同じく「# 何かの処理」を実行しようとした場合、 最初のスレッドが lock.release() を行うまで、後のスレッドは待たされることになる。 これにより「# 何かの処理」は、同時に単一のスレッドによってのみ実行されることが保障される。
昨日作成した knight11a.py では、計算中に「発見済の解の数」が表示される。 しかし、現時点で全体の何割程度の探索を終えたのかはわからないので、計算終了までにあと何時間かかるのか、あるいは何日かかるのかを、推定することができない。 そこで今度は、全体のうちどの程度の探索が終了したのか、おおまかな見積もりを表示することにしよう。
進捗状況を見積もる方法はイロイロあるだろうが、この見積もりのためにあまり計算量を増やしたくはない。 そこで今回は、2 手目のうち何% を探索し終えたか、によって進捗状況を評価することにしよう。 例えば 7 x 7 のチェス盤の場合、初手は 49 通りであり、各々の初手に 8 通りの 2 手目があるのだから、2 手目は全部で 392 通りである。 正確にいえば、この 392 通りの中にはチェス盤の外に移動しようとするような、現実には不可能な 2 手目も含まれているが、ここでは一旦、それは気にしないことにする。 この 392 通りの 2 手目のうち、探索しつくしたものの割合を表示することにしよう。 これを実現するには、knight8a.py に手を加える必要がある。knight8a.py を knight12a.py としてコピーし、グローバル変数 search を加えよう。
results = 0
search = 0
そして move() 関数の中に「global search」を追加して、さらに for ループの最後の部分を
# 今の手を取り消して、続きを検索する
matrix[x][y] = 0
if count == 2:
search += 1
continue
とする。これにより、search 変数は、何通りの 2 手目を探索し終えたのかを保持することになる。 次に、可能な 2 手目が何通りあるのかも計算しておこう。上でも述べたように、チェス盤のマス目の数に 8 を乗じただけでは、不可能な 2 手目も含まれてしまうため、過大評価になる。 この可能な 2 手目を数学的に計算するのは少し大変なので、直接数えてしまおう。 グローバル変数「progress_total = 0」を加えたうえで、関数 init() に
for pos_x in range(0, Size_x):
for pos_y in range(0, Size_y):
for shift in range(0, 8):
x = pos_x + Knight_x[shift]
y = pos_y + Knight_y[shift]
if x < 0 or x >= Size_x or y < 0 or y >= Size_y:
continue
progress_total += 1
を加える。これにより、「全てのプロセスにおける search の和」を progress_total で除することで、全体の進捗状況を推定することができる。 そのような集計処理を knight11a.py の Knight クラスや Progress クラスなどに追加し、knight12b.py としよう。
knight9a.py を、単にマルチプロセス化したものが上に示した knight10a.py である。 これでかなり高速化することができ、実時間 20 分ほどで 6637920 個の解を発見することができた。
knight10a.py では、初期座標ごとに 1 個のプロセスを作成し、全てのプロセスを並列で計算している。 そして定期的に途中経過を表示しているのだが、ここでは「既に終了したプロセスが発見した解の数」のみが示されている。 そのため、最初のうちは発見済の解の数として 0 が表示され続け、ある程度処理が進むと、急速に「発見済の解の数」が増えているようにみえる。 これはあまり美しくない。既に実行中のプロセスが発見した解も含めて表示する方が良いであろう。 これを実現するために、knight10a.py に次のような修正を加えて knight11a.py としよう。
この場合、マス目の数 (つまりプロセスの数) だけ共有メモリが必要になる。 この場合、multiprocessing.Value() で共有メモリを確保するのは大変なので、multiprocessing.Array() を使うと良い。 これは、共有メモリを配列として確保するものである。配列というのは、リストと同じようなものであるが、長さ (要素の数) を後から変更することができない (これまで触れていなかったが、リストは、作成した後に長さを変えることができる)。 例えば「n = multiprocessing.Array('d', 3)」とすれば、n は 3 個の要素から成る配列となり、n[0], n[1], n[2] として各要素にアクセスできる。
knight9a.py を元に、マルチプロセスを使うことで高速化した knight10a.py を作ってみよう。
さて、ナイトの周遊の高速化に話を戻そう。 現代のコンピューターは、大抵、複数コアから成る CPU を搭載しているため、複数の演算を同時に行うことができる。 従って、knight9a.py を実行している最中であっても、諸君の CPU にはかなりの余裕があるはずである。 Windows であれば、Ctrl + Alt + Delete からタスクマネージャーを開くことで、CPU の使用状況を確認することができる。
一つのスクリプト内で複数の計算を並列に行うことで、全体としての処理時間を短縮することができる。 27-1.py を次のように作ろう。
import time
from multiprocessing import Process
def Run():
time.sleep(3)
print("Done.")
p = Process(target = Run)
p.start()
q = Process(target = Run)
q.start()
p.join()
q.join()
ここでは単に「import multiprocessing」とするのではなく「from multiprocessing import Process」としている。 これは multiprocessing モジュール全体を読み込むのではなく、そのうち Process クラスだけを使用する、という意味である。 multiprocessing モジュールには、他にも Pool クラスなどが含まれている。Process と Pool の違いなどの詳細は、ここでは触れない。
プロセスというのは、コンピューター用語で、現に実行されている一つのプログラムのことである。 基本的には、一つのプログラムを実行すれば一つのプロセスが作成され、終了すればそのプロセスは消滅する。 一つのプログラムを二重に起動すれば、二つのプロセスが作成される。 しかし、プログラム内で意図的に複数のプロセスを作成し、各々独立して動くように設計すれば、一つのプログラムから多数のプロセスが形成される。 こうした処理方法をマルチプロセスと呼ぶ。 このマルチプロセスを行うためのモジュールが multiprocessing であり、27-1.py にあるように、Process クラスを使うことで新しいプロセスを作ることができる。 実行する関数に引数を渡したい場合は「p = Process(target = Run, args=(0,)」などと、args を指定すればよい。 ここで引数 0 の後に , を入れているのは、これがタプル、すなわち複数の値のセットであることを示すためである。とりあえずは、そういうオマジナイだと思っても構わない。 新しいプロセスで実行する関数を target 引数として渡し、start() でプロセスを開始する。join() はプロセスの終了を待つ関数である。
このマルチプロセスは、以前に使ったマルチスレッドとよく似ている。 しかし python の場合、マルチスレッドでは基本的に並列計算は行われず、つまり複数のスレッドが同時に計算処理を行うわけではないので、高速化には使えない。 一方、マルチプロセスの場合、各々のプロセスは独立して動くので、高速化に使える一方、メモリの内容を共有はしないので、データのやりとりには多少の工夫が必要になる。 簡単なのは、共有メモリを使う方法である。27-2.py を次のように作ろう。
import time
import multiprocessing
def Run():
time.sleep(3)
count.value += 1
print("Done.")
count = multiprocessing.Value('i', 0)
p = multiprocessing.Process(target = Run)
p.start()
q = multiprocessing.Process(target = Run)
q.start()
p.join()
q.join()
print(count.value)
今回は、冒頭部分を「import multiprocessing」としたので、Process() ではなく multiprocessing.Process() としなければならない。 Value クラスを使うことで、全てのプロセスからアクセスできる変数 count を作ることができる。 引数として指定した 'i' は、この変数は整数である、という意味である。小数を使いたい場合は 'd' とすればよい。第 2 引数の 0 は初期値を指定している。
knight8a.txt に対して、マルチスレッドを使い、10 秒毎に進捗状況を表示する機能を追加してみよう。これを knight9a.py とする。 なお「datetime」モジュールをインポートすると「datetime.datetime.now()」で現在時刻を取得できるので、これも表示に含めると読みやすいであろう。
24-1.py のマルチスレッドに少し手を加えてみよう。 France クラスの __init__() に「self.event = threading.Event()」を加え、さらに run() を次のように変更して 25-1.py を作ろう。
if self.event.wait(timeout=5):
print("TRUE", self.location)
else:
print("FALSE", self.location)
そして marseille.start() の後に「paris.event.set()」を入れてみよう。
threading モジュールに含まれている Event クラスを使うことで、スレッドの挙動を外部から制御することができる。 France クラスでは、__init__() 内で、data attribute として event を作成している。 これに対し event.set() を実行すると、その event について「フラグが立った」状態になる。 なお、event.clear() を実行すれば、フラグを取り消すことができる。 event.wait() は「フラグが立つまで待機する」という関数である。timeout を引数で指定している場合には、その時間 (単位は秒) だけ待ち、 フラグが立てば真、フラグが立たずにタイムアウトすれば偽を返す。
25-1.py では、paris については paris.event.set() が実行されているので、すぐに "TRUE Paris" が表示される。 一方で marseille については event() を実行していないので、5 秒経ってから "FALSE marseille" が表示されるのである。
以前に書いた knight8a.txt は、たとえば 6 x 6 のチェス盤なら実行に一時間以上を要した。 解を探している間、スクリプトは何の出力もしないので、はたしてキチンと動作しているのか、不安になってしまう。 そこで、knight8a.txt を修正して、1 分ごとに進捗を報告するようにしよう。 これを実現するには、runsingle() 関数の中に処理を書き足すという方法もあるが、ここではマルチスレッドを使うことにしよう。 マルチスレッドというのは、一つのスクリプト内で、処理を並行して行うという意味である。 これまでのスクリプトでは、関数などを使ってアッチコッチに飛んではいるものの、常に処理は一箇所のみで行われていた。これを、複数箇所で同時進行するのである。 24-1.py を次のように作ろう。
import time, threading
class France(threading.Thread):
def __init__(self, city):
threading.Thread.__init__(self)
self.location = city
def run(self):
time.sleep(5)
print(self.location)
paris = France("Paris")
marseille = France("Mariselle")
paris.start()
marseille.start()
paris.join()
marseille.join()
まず最初に、time モジュールと threading モジュールがインポートされている。 この threading モジュールに、マルチスレッド処理に必要なものが含まれている。 France クラスを定義する際に「class France(threading.Thread)」としている。これは、threading モジュールに入っている Thread クラスを継承する、という意味である。 継承というのは、元となるクラス (ここでは threading.Thread) と基本的には同じ内容だが、少しだけ違う部分のあるクラスを作る、という意味である。 元のクラスと違う部分だけを、France クラスの定義として書くのである。
__init__() は、そのクラスのオブジェクトが作られた際に実行される method である。 threading.Thread クラスにも __init__() method は含まれているので、France クラスの __init__() では、threading.Thread.__init__(self) を実行した後に、 独自の操作として self.location に city を代入している。
run() は、新しくスレッドが作られた際に実行される method である。 time.sleep() は、引数として与えられた時間 (単位は秒) の間だけ処理を停止する、という関数である。 この場合、5 秒間待ってから self.loation を出力している。
paris と marseilel は、いずれも France クラスのオブジェクトとして作られているが、それぞれ異なる引数が渡されている。 ここで渡した引数が、France クラスの __init__() の第二引数として渡されるわけである。
paris.start() と marseille.start() は、各々、新しいスレッドを開始している。 開始されたスレッドは、run() の中の time.sleep() を実行している間止まっているのだが、本体の方は start() をすぐに終え、次の処理に移っている。 すなわち、paris のスレッドと marseille のスレッドはほぼ同時に開始されているので、終了も概ね同時である。
paris.join() と marseille.join() は、それぞれ paris や marseille のスレッドが終了するまで待機する。
knight8a.txt のさらなる高速化を行う前に、クラスを触ってみよう。 クラスというのは、関数や変数などをグループにしたものである。 23-1.py を次のように作ろう。なお、これまでスクリプトの拡張子を .txt にしてきたが、世間では python スクリプトには .py を用いるのが主流なので、 我々もこれからは .py を使うことにする。なお、ファイルの中身はただのテキストファイルなので、これらの .py ファイルをメモ帳あるいはサクラエディタなど 好きなソフトウェアを使って編集すればよい。
class France():
location = 0
def Paris(self):
self.location = "North"
def Marseille(self):
self.location = "South"
nation1 = France()
nation2 = France()
print(nation1.location)
nation1.Paris()
nation2.Marseille()
print(nation1.location)
print(nation2.location)
23-1.py では、まず France というクラスを定義している。 この France クラスには location という変数が含まれていて、初期値は 0 である。 このような、クラスに含まれている変数のことを python 用語では data attribute と呼ぶらしい。 次に、France クラスには Paris() と Marseille() という二つの関数が含まれている。このような、クラスに属している関数はメソッド method と呼ばれるらしい。
スクリプトの本体部分では、まず「nation1 = France()」としている。これは、France クラスのオブジェクトを新しく作り、nation1 にそれを代入する、ということである。 オブジェクトというのは、そのクラスの型を有する変数である、と思ってよかろう。 つまり、nation1 はその中に location という data attribute や Paris という method を保持しているし、 nation2 も同様に location や Paris, Marseille を保持している。 「nation1.location」という表現は「nation1 が保持している location」という意味である。 「nation1.paris()」を実行すれば「nation1.location」は "North" になるが、nation2 には何の変化も起こらない。
メソッド Paris() や Marseille() の定義には self という引数が含まれている。 これらのメソッドの中では、この self を使うことによって、そのメソッドが含まれているオブジェクトを参照することができる。 すなわち、nation1.Paris() の中では self は nation1 を指しているし、nation2.Marseille() の中では self は nation2 を指している。 これらのメソッドを呼び出す際には、関数に引数を与える必要はなく、python が自動的に補完してくれる。
このように、クラスを使えば data attribute や method をまとめて操作することができるが、それだけであれば、わざわざクラスを使わなくてもリストなどで事足りる。 クラスが便利に使われる例を、次回、紹介しよう。
前回 knight7a.py をモジュールとして作成したので、これを使いまわすことで、新しいスクリプトを作りやすくなる。 ただし、前回の knight7a.py では、Size_x などの値を設定する初期化処理と、解を探索する処理とを同じ runsingle() 関数で扱っているため、 解を探索する目的で runsingle を実行するたびに Size_x などを初期化するという無駄があった。 また、一つの関数内で複数の質的に異なる処理を実行すると、それが何をする関数なのかわかりにくくなり、スクリプト全体が読みにくくなる。 そこで、初期化処理は init() 関数で行い、runsingle() 関数では解の探索だけを行うように修正しよう。 また、print() 関数による出力は削除してしまおう。 これを knight8a.py とする。
この新しいモジュール knight8a.py を用いることで、解を探索する処理の部分を毎回コピー & ペーストすることなく、様々なスクリプトに使いまわすことができる。 今回は、「チェス盤の大きさを指定すると、ナイトの周遊問題の解の総数を計算する」というスクリプトを書いて knight8a.txt としよう。 なお、初期座標は指定せず、あらゆる初期座標についての解の総和を計算するものとする。 knight8a.py を使うことで、knight8a.txt 自体は 20 行足らずに収まるはずである。
筆者の環境では、チェス盤が 5 x 5 の大きさの場合には 1 秒で 1728 通りの解を発見して終了した。 5 x 6 の場合は 40 秒で 37568 通りであった。 6 x 6 になると、69 分で 6637920 通りの解を発見した。
このように、チェス盤が大きくなっていくと、全ての解を探索するのに要する時間は膨大になってしまう。 7 x 7 までなら、何日も費やせば計算できそうであるが、8 x 8 の場合の計算は、かなり大変そうである。 どうすれば、計算時間をさらに短縮できるであろうか？
前回作成した knight6a.txt を元に、任意の大きさのチェス盤において、ナイトの周遊問題について全ての解を検索するスクリプトを書いてみよう。 knight6a.txt に手を加える形でもよいのだが、今回はモジュールを作成してみよう。
まず knight6a.txt をコピーして knight7a.py としよう。 そして、メインの処理を runsingle() 関数の中に入れる。すなわち
def runsingle(size_x, size_y, n, initial_x, initial_y):
global Size_x, Size_y, N_max, n_square, square
Size_x = size_x
Size_y = size_y
N_max = n
if Size_x < 1 or Size_y < 1:
print ("Error: size_x and size_y must be higher than 0.")
return 1
<中略>
print(results, "solutions found.")
return 0
といった具合にする。ついでに、move() 関数内の sys.exit(0) をなくすために、move() 関数内の該当箇所を
if results == N_max:
print(results, "solutions found.")
view(matrix)
return -1
matrix[x][y] = 0
return 0
if (move(matrix, count + 1, x, y) < 0):
return -1
と、しておこう。その上で、knight7a.txt を
import sys
import knight7a
# 適切な引数が与えられていなければエラーメッセージを表示して終了
if len(sys.argv) != 6:
print ("Usage: python {0:s} size_x size_y results initial_x initial_y".format(sys.argv[0]))
sys.exit(1)
Size_x = int(sys.argv[1])
Size_y = int(sys.argv[2])
N_max = int(sys.argv[3])
initial_x = int(sys.argv[4])
initial_y = int(sys.argv[5])
knight7a.runsingle(Size_x, Size_y, N_max, initial_x, initial_y)
と、することで、knight7a.txt から knight7a.py を呼び出し、knight6a.txt と同じ処理を行うことができる。
以前にも書いたが、こうした膨大な計算量を必要とする問題は python スクリプトで処理するべきではない。 knight4a.txt は初期座標 (3, 3) で計算完了までに 1 分 40 秒程度を要するが、 同じ内容の処理を C 言語でプログラミングすると、5.5 秒で完了した。 C 言語の方が 20 倍近く速いのである。
さて、8 x 8 のチェス盤においてナイトの周遊問題を考えると、解が非常にたくさんあるため、その全てを網羅するのはたいへん長い時間を要する。 そこで、今度は 7 x 7 の小さな「チェス盤」の上でナイトの周遊問題を解くことを考えよう。 スクリプトの仕様を以下のように変更する。
具体的には「python knight6a.txt size_x size_y results initial_x initial_y」と実行することで、 横 size_x, 縦 size_y のチェス盤において、初期座標を (initial\x, initial_y) として、results 個の解を発見するまで解を探索することにしよう。 results に負値を与えた場合には、全ての解を発見するまで探索する。 knights4a.txt を元に、必要箇所を修正したスクリプトを knight6a.txt としよう。
knight4a.txt における dead() は、手詰まり状態の検出、という点では不完全である。 というのも「行き止まりマスが 2 個以上存在する」以外にも、手詰まり状態は存在するからである。 たとえば、あるマスから移動できるマスの中に、route() 関数が 2 を返すようなマスが 4 個以上存在する場合について考える。 この状況から、そのマスに移動して、さらに別のマスに移動した後の状態を考えると、必ず行き止まりマスが 2 個作られることになる。 すなわち、上述のような状況では、行き止まりマスが存在しなくても、既に手詰まりなのである。 これを検出する処理を dead() 関数に追加して knight5a.txt を作ってみよう。
スクリプト開始部分の square 変数を定義している部分の下に
# route() 関数の結果を格納する
routemap = [[0 for i in range(0,8)] for j in range(0,8)]
として、route() の計算結果を格納する変数を作っておこう。そして dead() を
def dead(matrix, cur_x, cur_y):
global routemap
count = 0
# routemap のアップデート
for x in range(0, 8):
for y in range(0, 8):
# 行き止まり判定
if matrix[x][y] > 0:
routemap[x][y] = 0
continue
n = route(matrix, x, y, cur_x, cur_y)
routemap[x][y] = n
if n == 0:
return 1
if n == 1:
count += 1
if count > 1:
return 1
for x in range(0, 8):
for y in range(0, 8):
# 潜在的行き止まり判定
count2 = 0
for shift in range(0, 8):
test_x = x + Knight_x[shift]
test_y = y + Knight_y[shift]
if test_x < 0 or test_x >= 8 or test_y < 0 or test_y >= 8:
continue
if routemap[test_x][test_y] == 2:
count2 += 1
if count2 > 2:
count += count2 - 2
if count > 1:
return 1
return 0
と、変更しよう。ここでは for の二重ループが 2 回、登場している。 一回目では、knight4a.txt と同様の行き止まり判定を行うと同時に、route() 関数の実行結果を routemap に格納している。 これは、route() 関数はそれなりに計算量が多い関数であるため、何度も route() を呼び出すことを避ける目的である。 二回目の for 二重ループでは、上述のような「潜在的行き止まり」の判定を行っている。 この「潜在的行き止まり」と明らかな「行き止まり」の合計が 2 以上である場合を、手詰まりと判定することにした。
この knight5a.txt を筆者の環境で実行すると、100000 通りの解を発見するのに 3 分 3 秒を要した。knight4a.txt に比べ、かなり時間がかかったのである。 これは、上述のような「潜在的行き止まり」を発見することで短縮できる計算時間に比べて、潜在的行き止まりの判定に要する計算時間の方がずっと長くなってしまったからである。 他にも手詰まりの状況はいろいろと考えられるが、そのような比較的頻度の低い潜在的行き止まりについては考慮せず knight4a.txt のような単純なアルゴリズムで処理する方が良さそうである。
昨日作った knight3a.txt では、一つの解を発見するのに要する時間が短すぎて、アルゴリズムの改善効果を評価しにくい。 そこで、解を 100000 個発見したら処理を終了する、というように変更したスクリプト kngiht4a.txt を作ろう。 knight3a.txt で解を発見した場合に sys.exit(0) を呼び出していた部分を、次のように変更する。
# これが 64 手目であれば、チェス盤を全て埋めたので 0 を返す
if count == 64:
results += 1
if results == 100000:
sys.exit(0)
matrix[x][y] = 0
return 0
またスクリプトの開始部分で
results = 0
と初期化しておく。 しかし、これだけではうまく動かない。というのも、関数 move() 内で使われる変数は、基本的には関数内だけで有効だからである。 つまり、move() の外で変数 results に 0 を代入したとしても、関数 move() 内で使われる results 変数は、関数外で使われた変数 result とは同名の別物として扱われる。 従って、move() 内で results を参照した際には「まだ値が何も代入されていない変数」ということになってしまう。 また move() 内で results に値を代入しても、もう一度 move() を呼び出した際には、その代入はなかったかのように扱われる。 このような、関数内でのみ有効な変数のことを、プログラミング業界では「ローカル変数」と呼んでいる。
そこで、また move() 関数のはじめの部分に
global results
と追加する。これは、変数 results はグローバル変数として扱う、という宣言である。 グローバル変数というのは、ローカルでない変数のことである。つまり、これを宣言しておけば、move() 外で 0 を代入された results 変数を move() 関数内で参照できるし、 また move() 関数内で results += 1 を実行すれば、それは move() 関数外や、次回 move() を呼び出した時にも反映されるわけである。 グローバル変数とローカル変数の違いについて、簡単なスクリプトを書いて挙動を確認してみるとよい。
以上のようにして、初期座標 (3, 3) で 100000 個の解を発見するのに要する時間を測定すると、筆者の環境では 1 分 42 秒であった。
「手詰まり」というのがどういう状況であろうか。 あるマスからみて、ナイトの移動可能なマスが 1 個しかない場合、そのマスは「行き止まり」である。 行き止まりのマスは、64 手目で到着するしかない。 従って、盤面に「行き止まり」のマスが複数存在する場合、手詰まりであるといえる。 この手詰まりになった時点で探索をやめる、という方針にスクリプトを変更してみよう。
knight2b.txt を基に、move() 関数の中に、手詰まりかどうかを判定する処理を追加して knight3a.txt を作ってみよう。 まず、あるマスから移動できるマスの数を判定する関数 route() を次のように作ろう。
def route(matrix, pos_x, pos_y, cur_x, cur_y):
count = 0
for shift in range(0, 8):
x = pos_x + Knight_x[shift]
y = pos_y + Knight_y[shift]
if x < 0 or x >= 8 or y < 0 or y >= 8:
continue
if x == cur_x and y == cur_y:
count += 1
continue
if matrix[x][y] == 0:
count += 1
return count
この route() 関数は、そのマスから移動できるマスの数を返す。 すなわち、この関数が 1 を返したら行き止まりであるし、もし 0 を返すなら、そのマスには既にアクセス不能である。 なお、cur_x, cur_y として現在位置を渡すこととして、現在位置も「そのマスから移動できるマス」に含めることにしている。 従って、route() 関数が 0 を返すようなマスが 1 個でも存在するか、あるいは route() 関数が 1 を返すマスが複数存在する場合には、手詰まりであるといえる。 次に、この route() 関数を使って、「現在手詰まりであるかどうか」を返す関数 dead() を次のように作ろう。
def dead(matrix, cur_x, cur_y):
count = 0
for x in range(0, 8):
for y in range(0, 8):
if matrix[x][y] > 0:
continue
n = route(matrix, x, y, cur_x, cur_y)
if n == 0:
return 1
if n == 1:
count += 1
if count > 1:
return 1
return 0
dead() 関数は、matrix で示される盤面が手詰まりであれば 1 を、そうでなければ 0 を返す。 この dead() 関数を、move() 関数内の次の位置で呼ぶことにしよう。
# 動けるマスであれば、そこに動く
matrix[x][y] = count
# 動いた結果として手詰まりになるなら次へ
if dead(matrix, x, y):
matrix[x][y] = 0
continue
# これが 64 手目であれば、チェス盤を全て埋めたので終了
これにより、著者の環境では knight2b.txt と同じ解を返すのに要する時間が、初期座標 (3, 3) で 0.05 秒以下、初期座標 (2, 2) でも 0.2 秒以下となった。
まず簡単な高速化として、for ループを改善しよう。 knight1b.txt では x 座標と y 座標についてそれぞれ for ループを回し、全 64 マスについて「既に訪れたマスであるか否か」「動けるマスであるか否か」を判定した。 しかしナイトが動けるマスは、実際には高々 8 マスであり、その座標もすぐにわかるのだから、64 マス全てを検索するのは非効率である。たとえば
Knight_x = [-2, -2, -1, -1, 1, 1, 2, 2]
Knight_y = [-1, 1, -2, 2, -2, 2, -1, 1]
for shift in range(0, 8):
x = pos_x + Knight_x[shift]
y = pos_y + Knight_y[shift]
if x < 0 or x >= 8 or y < 0 or y >= 8:
continue
# ここで処理を行う
とすれば、ナイトが動けるマスだけを効率的に調べることができる。 knight1b.txt に対し、この改良を加えたスクリプトを knight2a.txt としよう。 筆者の環境では、knight2a.txt は初期座標 (3, 3) なら 7 秒、初期座標 (2, 2) でも 2 分 15 秒で、一つの解を発見することができた。
さらなる高速化のために、なぜ、このスクリプトは解を発見するまでにこれほど長い時間を要するのか調べてみよう。 move() 関数の for ループの前に
print("DEBUG: {0:d} ({1:d}, {2:d})".format(count, pos_x, pos_y))
を加えて knight2b.txt としよう。これにより、move() 関数の実行状況が表示されるので、スクリプトの詳細な挙動を追跡することができる。 このスクリプトをそのまま実行すると、表示が多すぎてみにくいので、cyginw 上で「python knight2b.txt 3 3 > log.txt」として実行すると良い。 これにより、通常なら画面に表示される内容が「log.txt」に書き込まれるようになる。 筆者の環境では 183 MB の log.txt が生成された。 以下の例は筆者の書いた knight2b.txt を初期座標 (3, 3) で実行した場合の話である。諸君が筆者とは少し違うスクリプトを書いた場合には、 筆者とは違う log.txt が生成されることに留意されたい。
log.txt の冒頭部分をみてみよう。DEBUG: に続く値は「何手目を探索しているか」であり、括弧内は「現在どの座標にいるか」である。 2 手目から 41 手目までは順調に進んでいったが、「41 (7, 0)」の次が「41 (7, 4)」になっている。 つまり、41 手目として (7, 0) に進んだが、そこで動けるマスがなくなってしまったので、41 手目を (7, 4) としてやり直したわけである。 その後は 56 手目で (6, 0) に行って止まり、それからは 50 手目以降あたりをウロウロと模索している。 さらにみると、log.txt の 108 行目では 57 手目を探索しているのに、その後は 136 行目で 46 手目まで戻り、さらに紆余曲折の末、5976890 行目では 36 手目まで戻っている。 この「57 手目まで進んでから 36 手目まで戻る」という動きは、甚だ非効率であるようにみえる。 そこで状況を解析するために、任意の手数まで進んだ時点でスクリプトが止まるように変更しよう。 スクリプトの、コマンドライン引数を処理するあたりを
if len(sys.argv) != 4:
print ("Usage: python {0:s} initial_x initial_y moves".format(sys.argv[0]))
sys.exit(1)
initial_x = int(sys.argv[1])
initial_y = int(sys.argv[2])
moves = int(sys.argv[3])
として、move() を呼ぶ部分を
move(square, 2, initial_x, initial_y, moves)
とし、さらに結果表示の部分を
# これが 64 手目であれば、チェス盤を全て埋めたので終了
if count == 64 or count + 1 == moves:
view(matrix)
sys.exit(0)
これらの修正を加えたスクリプトをknight2c.txt とする。筆者が書いたものを参考にしてもよい。 これを「python knight2c.txt 3 3 56」として実行すると、56 手目を探索する段階、すなわち log.txt でいう 57 行目時点での状況が表示される。 それによると、盤面左下 (0, 6) の位置に「55」があるが、そこから移動できる場所に空き、すなわち 0 のマスが存在しない。 それ故、スクリプトは探索を断念して、もう少し手前の状態から再開したのである。 さらにみると、盤面の左下と右下に 0 すなわち未だ訪れていないマスが分断されて残っている。 このような状況では、左下と右下の両方を踏破することは不可能である。 そう考えると、55 手目まで進んで「56 手目を発見できないゾ」となってから探索を断念するのではなく、 0 のマスが複数個所に分断されてしまった時点で探索を諦める方が効率的であろう。 それを実現するようなアルゴリズムを考えてみよう。
前回作成した knight1a.txt をコピーして knight1b.txt とし、move() 関数の中身を次のように作ろう。 ついでに、解答を表示するための関数 view() も作っておく。 スクリプトの本体部分は前回と同じなので、再掲しない。
def view(matrix):
for x in range(0, 8):
for y in range(0, 8):
print("{0:2d}".format(matrix[x][y]), end=" ")
print()
return 0
def move(matrix, count, pos_x, pos_y):
"次の一手を探す関数"
for x in range(0,8):
for y in range(0, 8):
# 既に訪れたマスであれば次へ
if matrix[x][y] > 0:
continue
# 動けるマスであれば、そこに動く
if (abs(x - pos_x) == 1 and abs(y - pos_y) == 2) or (abs(x - pos_x) == 2 and abs(y - pos_y) == 1):
matrix[x][y] = count
# これが 64 手目であれば、チェス盤を全て埋めたので終了
if count == 64:
view(matrix)
sys.exit(0)
move(matrix, count + 1, x, y)
# 上の move() がチェス盤を埋め尽くして終了 (sys.exit(0)) した場合は、ここより下は実行されない
# 今の手を取り消して、続きを検索する
matrix[x][y] = 0
continue
return -1
# で始まる行は、コメントである。Python がスクリプトを実行する際にはコメント行は無視されるので、スクリプトをみやすくするために、適宜コメントを入れると良い。 knight1b.txt が何をやっているのかは、初心者にはわかりにくいかもしれない。 move() 関数には 4 個の引数が渡される。最初の matrix はチェス盤上のこれまでの移動の履歴を含む二次元リストであり、実際は square と同じである。 square と matrix は同じなのだから、わざわざ move() 関数に matrix として渡さなくても、直接 square を操作すれば良いのだが、 敢えてここで別の変数にしておいた方が後々役立ちそうな気がするので、分けておくことにする。このあたりは、経験と勘の世界である。 count は「現在、何手目を探索しているか」を示す値である。pox_x, pox_y は、それぞれナイトの現在位置の x 座標、y 座標である。
move() 関数には 2 個の for 文が含まれているが、これはチェス盤のマスを端から順番にみていくものである。 matrix[x][y] が 0 より大きいならば、そのマスは既に訪れたマスなので、continue することで次のマスを調べることになる。 abs() 関数は絶対値を返す関数である。この長い if 文は、「現在位置からナイトが移動できるかどうか」を評価しているわけである。 もし移動できるならば、そこに移動する、つまり matrix[x][y] に count を代入する。 これが 64 手目であるならば、view() 関数によって解を表示して、終了 (sys.exit(0)) する。 そうでなければ、count に 1 を足して move() を実行することによって、次の手を探索する。 このように move() の中で move() 自身が呼び出されているが、このような方法を再帰的、と呼ぶ。 move() 関数は解を発見したら sys.exit() するので、もし move() より後の行が実行されたならば、それは move() 関数が解を発見できなかっとことを意味する。 その場合、先ほどの移動を取り消し、つまり matrix[x][y] を 0 に戻して、次のマスを探索することになる。 view() 関数では、print() の中に「end=" "」が書かれている。これは、print() はデフォルトでは末尾に改行を付すが、これを「" "」に変更する、という意味である。
筆者の環境では、初期座標を (3, 3) にした場合、この knight1b.txt は 30 秒で 1 個の解を提示した。 しかし、これは運が良かっただけであり、初期座標を (2, 2) にすると 1 個の解を発見するのに 9 分 30 秒を要した。
この「ナイトの周遊」問題は、私が工学部 1 回生だか 2 回生だった頃 (つまり 2002 年か 2003 年頃) に、情報処理の授業で扱った問題である。 当時の講師は「人間が感覚的に解くならばそれほど難しくないが、コンピューターで解こうとするとかなり時間がかかる例」として、この問題を示していた。 実際、現在我々が書いているスクリプトを 20 年前のコンピューターで実行したなら、一つの解を示すだけでも、結構な時間がかかったであろう。 しかし当時の私は、その講師の主張を疑い、コンピューターの方が人間よりも速く解けるのではないかと考えた。 そこでアルゴリズムに工夫を凝らし、当時のコンピューターで一秒間に 1 万ほどの解を発見するプログラムを書くことに成功した。 それを同級生の某君に自慢したところ、その某君の友人である情報科の某学生が関心を示してくれたので、私はその学生氏に自分が書いたプログラムを送った。 その思い出のプログラムを、以後何回かかけて、python で復活させよう。
ここからは、より複雑なスクリプトを扱うことにしよう。 ナイトの周遊問題について考える。 チェスにおけるナイトは、上のリンク先の図に示されているように、将棋における桂馬と類似の動きをするが、桂馬とは異なり、全方向すなわち最大 8 箇所に移動できる。 チェス盤は 8 x 8 の大きさである。このチェス盤の任意の位置からスタートし、全てのマスを 1 回だけ通過する方法を示せ、というのがナイトの周遊問題である。 同じマスを 2 回以上通ってはならないし、一度も通らないマスがあってもいけない。 この問題の解答を全て検索するスクリプトを書いてみよう。
こういう問題を扱う場合、どのような戦略で解答を探していくかによって、必要な計算時間が大きく変わってくる。 その探索の戦略をアルゴリズムと呼ぶ。 まずは、効率が悪いことを承知の上で、次のような方針でスクリプト knight1.txt を書いてみよう。
この内容を一気に作るのは大変なので、まず、以下の内容で仮のスクリプト knight1a.txt を作ろう。 これは、ナイトの初期位置をコマンドライン引数から受け取り、それが有効な座標か否かを評価し、場合によってはエラーメッセージを表示して終了する、というものである。
import sys
def move(matrix, count, pos_x, pos_y):
"次の一手を探す関数"
return 0
square = [[0 for i in range(0,8)] for j in range(0,8)]
if len(sys.argv) != 3:
print ("Usage: python {0:s} initial_x initial_y".format(sys.argv[0]))
sys.exit(1)
initial_x = int(sys.argv[1])
initial_y = int(sys.argv[2])
if initial_x < 0 or initial_x >= 8:
print ("Error: initial_x must be between 0 and 7.")
sys.exit(1)
if initial_y < 0 or initial_y >= 8:
print ("Error: initial_y must be between 0 and 7.")
sys.exit(1)
print ("The initial location is ({0:d}, {1:d}).".format (initial_x, initial_y))
square[initial_x][initial_y] = 1
move(square, 2, initial_x, initial_y)
square が 8 x 8 の二次元リストとして定義されている。ここでは全ての要素を 0 で初期化しているが、いささか癖のある記法である。 こういう書き方があるのだ、とだけ認識していただければ十分である。 sys.exit() は、スクリプトを終了するための関数である。 この関数に 1 を渡しているのは、正常終了ではなくエラー終了であることを示すためである。 慣習的に、正常終了では 0 を、エラー終了では非 0 を、渡すことが多い。 knight1a.txt では、とりあえず中身がカラッポの状態で move() という関数が定義されている。次回は、この関数を実装していくことにしよう。
文字や数字を出力する際に、単純に print() を使うと表示が美しくならず、読みにくいことがある。 13-1.txt を次のように作ろう。
for n in range(0, 5):
print(n**2, n**3)
for n in range(0, 5):
print("n**2 = {0:2d}, n**3 = {1:02d} ({1:.3f})".format(n**2, n**3))
前半部分では単に print() しているので、桁数の具合によっては表示が乱れ、読みにくい。 それを整えるためには、後半部のように format を使うのである。 前半の「("n**2 = {0:2d}, n**3 = {1:02d} ({1:.3f})"」の部分が書式を表しており、それに続く .format() によって、 書式指定文字列中の「{0:2d}」などの部分が、書式を整えられた数値に置き換えられる。 たとえば「{0:2d}」についてみると、最初の「0」は、format() に渡された引数のうち 0 番目の値を表示する、という意味であり、具体的には「n**2」に置き換えられる。 「:」の後の部分が書式指定であって、「2d」の「d」は「整数として扱う」ということを、「2」は桁数を、それぞれ表す。 すなわち「0:2d」というのは、0 番目の値を 2 桁の整数として表示する、という意味である。 この場合、0 は「 0」というように、スペースを空けて 2 桁分の幅をもって表示される。 同様に、この例では「1:02d」は「n**3」を 2 桁の整数として表示することになる。 「f」は、小数として扱うことを示す。「.3f」とすれば、小数として、小数点以下 3 桁まで表示する、ということになる。 また、非常に大きな数や非常に小さな小数を扱う場合には「e」が便利である。 他にもいろいろな書式があるので、Python のライブラリなどを参考にするとよい。
Python でファイルを読み書きすることもできる。 13-2.txt を次のように作ろう。
text1 = "Isaac is an English scientist.\n"
text2 = "Pierre was French.\n"
file = open("test.txt", "w", encoding="utf-8")
file.write(text1)
file.write(text2)
file.close()
これを実行すると「test.txt」が作成されるはずであるので、中を確認してみよう。 open() は、ファイルを開く関数である。開かれたファイルは、以下、この「file」変数を用いて操作することができる。 open() に渡した最初の引数は、ファイル名である。 2 個目の引数はモードであって、"w" とすれば書き込みモード、"r" とすれば読み込みモードとなる。 書き込みモードで開いた場合、既存のファイル内容は破棄されるので、もし追記したい場合には "a" を指定して追記モードで開く必要がある。 3 個目の引数はテキストのエンコード指定であるが、現在では特別な事情がない限りは utf-8 を指定するのがよい。 file.write() により、文字列を file に書き込むことができる。 処理が終わったら file.close() 行う。これを忘れると、ファイルへの書き込み処理がキチンと行われない恐れがある。 というのも、ハードディスク等への書き込みは、効率化の関係上、file.write() が実行されたタイミングではなく、ある程度まとめてから実施される。 そのため、file.close() がないと、「後で書こうと思っていたものが、実際には書き込まれないまま終わる」という事件が起こるのである。
ファイルを読む際には、read() を用いる。13-3.txt を次のように作ろう。
file = open("test.txt", "r", encoding="utf-8")
text = file.readline()
file.close()
print(text)
file.readline() は、file から一行だけ読み込む、という関数である。 二行目以降も読みたければ、readline() を繰り返せばよい。
12-1.txt を次のように作ろう。
text = input()
print ("Your input is \"", text, "\"")
これを実行すると、Python は実行中に停止する。 そこで諸君が何か好きな文字を入力すると、Python はその内容を反復して終了する。 これを用いれば、諸君の入力に応じて異なる挙動をするスクリプトを書くことができる。
似たようなことは、次のようなスクリプトでも実現できる。これを 12-2.txt としよう。これを「python 12-2.txt Helsinki」などとして実行されたい。
import sys
print (sys.argv[0])
print (sys.argv[1])
Python の実行時に、スクリプト名の後に書いた内容は「sys.argv」リストに格納されている。sys.argv[0] は大抵、スクリプト名である。
ところで、スクリプトにランダムな挙動をさせたいことがあるかもしれない。その場合には、次に示す 12-3.txt のような方法でランダム (にみえる) 数を生成することができる。 これは 0.0 以上 1.0 未満の数を表示するスクリプトである。実行するたびに異なる結果になることを確認されたい。
import random
print (random.random())
これらの機能を用いて、次のようなゲームを作ってみよう。 「まずスクリプトは 0 以上 100 以下の整数一つを正解として決める。プレイヤーが整数を入力すると、スクリプトは、それが正解であるか間違いであるかを答える。 間違いである場合には、大きすぎるのか小さすぎるのかも答える。プレイヤーが 6 回以内に正解できれば、勝ちである。」 これを 12-4.txt とする。
関数の中で変数を取り扱う際には、少し気を付けなければならない点がある。11-1.txt を次のように作ろう。
a = 3
b = 4
def function(a):
"This is a test of function."
b = 2
print (a, b)
function(5)
print (a, b)
関数の中の一行目に書かれている文章は、スクリプトの実行内容には影響を与えない。 スクリプトをみやすくするために、これが何をする関数なのか簡潔に説明したコメントである。Python 用語では「ドキュメンテーション文字列 documentation string」と呼ばれる。
11-1.txt では、最初に a に 3 が代入されているが、function 関数には 5 が渡されている。 従って、関数の中では a には 5 が代入された状態になっている。さらに b に 2 を代入してから print() しているので「5, 2」と表示されるわけである。 一方、a に 5 が代入されたり b に 2 が代入されたりといった操作は、この関数の中でのみ有効である。 すなわち、関数が終わった後の最終行の print() の時点では、最初の二行で a には 3 を、b には 4 を代入したままの状態になっているので、今度は「3, 4」と表示される。 非常に紛らわしく、バグの温床にもなるので、このようにスクリプト全体にわたって有効な変数 (グローバル変数、と呼ばれる) と、 関数の中でのみ有効な変数 (ローカル変数、と呼ばれる) で同じ名前を使わない方がよい。 つまり、11-1.txt は、次のように修正した方が安全である。これを 11-2.txt としよう。
a = 3
b = 4
def function(n):
"This is a test of function."
m = 2
print (n, m)
function(5)
print (a, b)
スクリプトの挙動は 11-1.txt と 11-2.txt で全く同じであるが、11-2.txt の方が読みやすく、間違いが起こりにくい。
関数に渡す値 (引数、と呼ばれる。これを「いんすう」と発音する人もいるが、「ひきすう」と読むのが多数派である。) については、デフォルト値を設定することもできる。 11-3.txt を次のように作ろう。
def function(n = 2):
print (n)
function()
function(5)
function()
function を定義する際に「n = 2」などとしておくと、この関数を呼び出す際に引数を省略できる。 その場合、この予め指定された値がデフォルト値として指定される。
次のような問題を考えよう。 「10 から 100 までの整数について、それが素数であるならば、その旨を表示する。 また、それが偶数であり、かつ、それを 2 で割った数が素数であるならば、その旨を表示する。」
前回までに、ある整数が素数であるか否かを判定するスクリプトを書き、しかも高速化を行ってきた。 それを踏まえれば、この問題を解くこと自体は簡単である。 たとえば、9-1.txt を下地として、以下のような具合のスクリプトを書けばよい。これを 10-1.txt とする。
import math
for n in range (10, 101):
if n % 2 == 0:
m = n / 2
for t in range (2, int(math.sqrt(m)) + 1):
if m % t == 0:
break
else:
print ("A half of", n, "is a prime number.")
else:
for t in range (2, int(math.sqrt(n)) + 1):
if n % t == 0:
break
else:
print (n, "is a prime number.")
しかし、このスクリプトはあまり美しくない。 というのも、ある値が素数であるか判定するための、全く同じ for ループが 2 回登場しているからである。 これでは、スクリプトが読みにくいだけでなく、この判定部分に後から何か修正を加えようとしたとき、2 箇所で同じ修正作業を行わなければならない。 これは、面倒な上に、修正漏れによるバグが生じる温床となる。 同じ作業を行う部分は、一箇所にまとめて書いた方が、間違いが少ないのである。 そこで、関数と呼ばれる機能を使うことにする。上のスクリプトを修正して、次のように 10-2.txt を作ろう。
import math
def prime(n):
for t in range (2, int(math.sqrt(n)) + 1):
if n % t == 0:
return 0
else:
return 1
for n in range (10, 101):
if n % 2 == 0:
m = n / 2
if prime(m):
print ("A half of", n, "is a prime number.")
else:
if prime(n):
print (n, "is a prime number.")
「def」は、関数を定義する命令である。ここでは「prime」という関数を定義しており、この関数には引数として「n」を渡すことにしている。 この関数は、n が素数であれば 1 を返して終了し、素数でなければ 0 を返して終了する。
10-1.txt と 10-2.txt は同じことを行っているのだが、10-2.txt の方が、何をやっているのか多少みやすく、メンテナンスも行いやすい。
Cygwin 上で、プログラムの実行に要した時間を計測するには「time」コマンドを用いる。 たとえば「time python 9-2.txt」を実行すれば、「python 9-2.txt」を実行するのに要した時間が表示される。 これは、実世界での経過時間 (real) と、コンピューター上でプログラムの実行に要した時間 (user) とが区別される。 後者は、実際に CPU を占有していた時間、と考えてよい。 つまり、もし他のプログラムを動かしていて CPU がたいへん忙しいような状況の場合、real の時間は長くなるが、user は変わらない。
前回の 8-3.txt は、実行するのにたいへん時間がかかったであろう。著者の環境では 118 分を要した。 はじめのうちは、そこそこ軽快に処理が行われるのだが、やがて動きが鈍くなり、 90000 番目の素数を計算する頃には、プログラムがフリーズしたのではないかと諸君は心配したことであろう。 なぜ、こうなるのかというと、while 文の中で n が大きくなってくると、for 文の反復もそれに応じて多くなり、膨大な計算量が必要になるからである。
そこで今回は、8-3.txt を改良し、スクリプトを高速化してみよう。 8-3.txt では、t は 2 から n の範囲を動いていたが、よく考えると、これは無駄である。 ほんとうは、t は 2 から "n の平方根" までの範囲を動けば十分である。 というのも、もし n が "n の平方根よりも大きな整数" で割り切れるならば、その商は "n の平方根よりは小さな整数" のはずである。 従って、もし n が "n の平方根よりも小さな整数" で割り切れないならば、"n の平方根よりも大きな整数" でも割り切れないのである。 そこで、t の範囲を狭めて、9-1.txt を次のように作る。今回初登場の表現を太字にした。
p = 100000
import math
n = 2
count = 0
while 1:
for t in range (2, int(math.sqrt(n)) + 1):
if n % t == 0:
break
else:
count += 1
print (n, "is the ", count, "th prime number.")
break
n += 1
「import math」は、Python に「math」モジュールを読み込ませる命令である。 これは、複雑な数学的処理をするためのモジュールであって、次に出てくる「math.sqrt()」の実行に必要である。とりあえず、そういうオマジナイだと思っておいて構わない。 「math.sqrt(n)」は「n の平方根を計算する関数」である。そして「int()」は、小数を切り捨てにより整数に変換する関数である。 平方根は一般には小数になるが、range() 関数は整数しか扱えないので、小数を整数に変換する必要があるため、math.sqrt() の結果を int() に入れたのである。
筆者の環境では、9-1.txt の実行は 10 秒で済んだ。しかし、より大きな素数を探索するには、より長い時間を要する。 たとえば 100 万番目の素数を計算するのに、筆者は 5 分 45 秒かかった。 これでは、1 億番目の素数を計算するのはたいへんである。さらに高速化することを考えよう。
9-1.txt では、t を 2 から順に 1 ずつ大きくしているが、これは無駄である。 たとえば 2 で割り切れなかったなら、4 でも割り切れないのは当然である。 効率的に素数を探すなら、これまでに発見した素数で n を割り切れるかどうかだけを調べればよい。そこで、9-2.txt を次のように作ろう。
p = 1000000
import math
n = 2
count = 0
prime = [0] * p
while 1:
limit = int(math.sqrt(n))
t = 0
while t < count and prime[t] <= limit:
if n % prime[t] == 0:
break
t += 1
else:
prime[count] = n
count += 1
if count == p:
print (n, "is the ", count, "th prime number.")
break
n += 1
「prime = [0] * p」というのは、p 個の要素から成るリストとして変数 prime を作り、全ての要素に 0 を代入する、という意味である。 この prime 変数に、発見した素数を順次格納していく。 また、while 文の中で「and」が登場しているが、これを用いることで、「t < count」かつ「prime[t] <= limit」の場合に while ループを回す、ということができる。 同様に「or」を使えば、いずれかを満足する場合に while ループを回す、ということもできる。 この 9-2.txt は、筆者の環境では 2 分 40 秒で完了した。 このように、処理の内容は同じであっても、スクリプトの書き方次第で処理速度は大きく異なる。
なお、計算量が多く、実行に長い時間を要する処理を Python で行うのはお勧めしない。 9-2.txt と同じ内容のプログラムを C 言語で書いたところ、13 秒で実行完了した。C 言語の方が 10 倍、速かった、というわけである。
Python で不適切な while 文などを書いてしまうと、ループが永久に終わらず、コンピューターが延々と計算し続ける、という事態が生じてしまう。 たとえば、0 から 100 までの整数を表示させるつもりで、うっかり
loop = 0と書いてしまうと、延々と、いつまでも 0 を表示し続けることになる。こうした場合には、Ctrl + C で Python を強制終了させることができる。
さて、while 文や for 文の便利な使い方を紹介しよう。8-1.txt を次のように作る。
loop = 0
while 1:
print (loop)
loop += 1
if loop > 5:
break
for loop in range(0, 5):
print (loop)
while の条件式には、このように数値を入れることもできる。この場合、その数値が非 0 であれば、条件が満足されたと判断される。 また、for 文に range() を用いたが、これは、指定された範囲の数列を作成する関数である。 すなわち「range (0, 5)」は「[0, 1, 2, 3, 4]」と同じように扱われる。
これらを踏まえて、次の内容で 8-2.txt を作ろう。これは、1000 番目の素数を表示するものである。
p = 1000
n = 2
count = 0
while 1:
for t in range (2, n):
if n % t == 0:
break
else:
count += 1
if count == p:
print (n, "is the ", count, "th prime number.")
break
n += 1
一行目で p に代入する値を変えれば、何番目の素数でも計算することができる。 試しに、8-3.txt を次のように作って、10 万番目の素数を計算してみよう。 なお、進捗を表示するために、太字部分を追加した。
p = 100000
n = 2
count = 0
while 1:
for t in range (2, n):
if n % t == 0:
break
else:
count += 1
if count == p:
print (n, "is the ", count, "th prime number.")
break
elif count % 1000 == 0:
print (n, "is the ", count, "th prime number.")
n += 1
Cygwin でファイルやリンクを削除するコマンドは rm である。 たとえば、ホームディレクトリに users というファイルまたはリンクがあるとき「rm users」を実行すれば、users は削除される。 「本当に削除しますか？」というような確認は表示されないし「ゴミ箱」のようなものもないので、削除したファイルは基本的には元に戻せない。注意されたい。
これまでにみてきた構文を適宜用いて、次のような挙動をするスクリプトを書いてみよう。 「5 以上 100 以下の全ての整数 n について、それが 2 または 3 の倍数であれば、その数を表示する。そうでなければ、何もしない」これを 7-1.txt としよう。
次に、1048573 が素数であるかどうかを、Python を用いて検証してみよう。これを 7-2.txt とする。
Cygwin では「ln」コマンドでリンクを作ることができる。これは Windows でいうショートカットのようなものである。 たとえば、ホームディレクトリで「ln -s /cygdrive/c/Users/ users」と実行すれば、ホームディレクトリに users という名前のリンクが作成される。 そこで「cd users」と実行すれば、/cygdrive/c/Users/ に移動することができる。
前回は while 文と if 文を使った。今回は if 文に変化をつけてみよう。6-1.txt は次のように作る。
loop = 0
while loop < 10:
if loop < 3:
print (loop, "is smaller than 3.")
elif loop < 7:
print (loop, "is smaller than 7.")
else:
print (loop, "is not smaller than 7.")
loop += 1
if に続く条件が満足されていれば、それに続く命令が実行される。 if の条件が満足されていないならば、elif の条件が評価される。 elif の条件も満足されないならば、else に続く命令が実行される。 elif は何個続けても構わないし、if と else だけでも問題ない。むろん、if だけでも構わない。 これらの構文を用いることで、条件によって細かに異なる挙動をさせることができる。
以前に紹介したリストの各項目について処理を行うには、for 文が便利である。6-2.txt は次のように作ろう。
list = [1, 4, 9, 16, 25, 36]
for value in list:
print (value + 1)
while や for によるループは、途中で抜けることもできる。6-3.txt を次のように作ろう。
loop = 0
while loop < 10:
if loop % 2 == 0:
loop += 1
continue
elif loop > 5:
break
print (loop)
loop += 1
continue は、for や while の中で使われると、残り部分を無視して次のループを開始する。 break は、その時点でループを抜ける。
また、else は for や while と組み合わせることもできる。6-4.txt を次のように作ろう。
loop = 2
while loop < 5:
if loop > 5:
break
loop += 1
else:
print ("I am Isaac.")
loop = 2
while loop < 10:
if loop > 5:
break
loop += 1
else:
print ("I am Gottfried.")
この場合の else は、break されずにループが終了した場合にのみ実行されるのである。
Cygwin でホームディレクトリ (Cygwin を起動した時にいるディレクトリ) に移動するには、単に「cd」を実行すればよい。
前回までの内容は、「手作業でやった方が速そうな気がするが、Python で実行することもできる」ような作業であった。 今回は、「手作業でもできるが、とても面倒なので、Python にやらせた方が楽」な処理を扱おう。 5-1.txt は次のように作る。while の次の行の頭にはスペースがあいているが、これは tab を使うか、複数の半角スペースを並べる。
loop = 0
while loop < 10:
print (loop)
loop += 1
print ("END")
「while」文は、その直後に示されている条件が成立している限り、次行以降の tab でインデントされている内容を反復して実行する。 行末の「:」は、条件文の終わりを示す記号である。 5-1.txt の例では、loop が 10 未満である限り、loop の内容を表示して loop を 1 増やす、というわけである。 「loop += 1」は、「loop = loop + 1」と同じ意味である。 while 文の中で使える条件としては、不等号の「<」や「>」の他、等しいことを示す場合には「==」を使うことができる。 単に「=」とすると代入の意味になってしまうので、数学的にはキモチワルイ記号であるが、「=」を二つ重ねて「==」とするのである。 また、以上、以下を示すには「<=」や「>=」が、等しくないことを示すには「!=」を使うことができる。
while と似たような構文として if がある。これは、条件が満足されている場合に限り、次行以降の tab でインデントされている内容を実行する、というものである。 5-2.txt を次のように作ろう。
n = 2
if n < 3:
print ("2 < 3")
n = 5
if n < 3:
print ("5 < 3")
この while 文や if を使って、100 以下の素数を全て表示するスクリプトを書いて 5-3.txt としよう。 素数とは、2 以上の整数のうち、1 と自分自身以外では割り切れないような数のことである。 このスクリプトの例は、次回、示すことにしよう。
これは、プログラミング初心者には難しい問題かもしれない。イマイチ解決方法が思い浮かばないなら、とりあえず現時点では触らずにおいても構わない。
Cygwin で現在のディレクトリを確認したい場合には pwd コマンドを使うとよい。
複数の変数をひとまとめに扱うと便利なことがある。4-1.txt を次のように作ろう。
number = [1, 4, 9, 16, 25, 36]
print (number[0], number[1], number[3])
このリストは、連結したり、内容を書き換えたりすることができる。4-2.txt は次のようにしよう。
number = [1, 4, 9, 16, 25, 36]
append = [1, 8, 27, 64, 125, 216]
result = number + append
print (result[0], result[6], result[8])
result[6] = 0
print (result[0], result[6], result[8])
result[1:4] = [10, 11, 12]
print (result[0], result[2], result[5])
さて、4-3.txt を次のようしたとき、何が起こっているか、わかるだろうか？
number = [1, 4, 9, 16, 25, 36]
number [1:3] = [0, 0, 0]
print (number[4])
number = ["a", "b", "c", "d", "e", "f"]
number [1:4] = ["X", "X"]
print (number[0], number[2], number[4])
Cygwin で、たとえば現在 /cygdrive/c/Users/ にいるとして、/cygdrive/c/ に移動したい、という状況を考えよう。 この場合「cd /cygdrive/c/」を実行すればよいのだが、いささか入力が面倒である。 こうした場合「cd ../」とする方が簡便である。 「..」というのは「親ディレクトリ」を意味する表現であって、Windows 風に言えば「一つ上のフォルダ」のことである。
さて、本日は Python におけるテキストの取り扱いをやってみよう。 print() で文字列を表示できることは既にみてきたが、表現の仕方にはいくつかのバリエーションがあるので、好きなものを使えばよい。3-1.txt は次のように作ろう。
print ("This is a pen.")
print ('What is that?')
print ("It's an apple.")
print ("I don\'t see it.")
以前にも述べたが、「"」で括られた部分は一塊の文字列として取り扱われる。「'」を用いても同じことである。 表示したい文中に「'」や「"」が登場する場合には、いささかの注意が必要である。それらは文の一部に過ぎず、文字列の終わりではないことを Paython が理解できるように指示しなければならないからである。 具体的には、1 行目のように「"」と「'」で区別することができるし、あるいは 4 行目のように「'」の前に「\」を置くことで、「これは特別な意味のない文字ですよ」と明示することができる。
文字列を表示する際には、次のような表現もできる。3-2.txt を以下のように作ろう。
print ("Hello.\nAre you Tom?")
print ("""
No.
I am Tomson.
""")
print ("""\
Oh, I am sorry.
""")
print (r"\n\n\n")
まず 1 行目の例が示すように、「\n」は改行として扱われる。 改行を扱う方法としては 2-5 行目のように、文字列の始まりと終わりを「"""」にする、という方法もある。この場合、スクリプト中の改行が実際の出力に反映されるわけであるが、 6 行目のように行末に「\」を置いておけば、スクリプト中の改行を実際の処理上では無視させることもできる。 あるいは文字列の始まりの「"」の前に「r」を置いておけば、文字列中の「\n」を、改行ではなくそのまま「\n」という文字列として扱わせることもできる。 なお、フォントによっては「\」は円記号で表示されたり、バックスラッシュで表示されたりするが、意味は同じである。
文字列に対して、ちょっとした演算のようなことを行うこともできる。3-3.txt は次のようにしよう。
print ("Hello." + " " + "I am Tom.")
print ("Hello. " "How are you?")
print ("Hello. " * 3)
print (4 * "Hello. ")
text = "Hello. "
print (text + "Today is a fine day.")
文字列は「+」でつなぐこともできるし、2 行目のように連続して並べることで暗黙のうちに連結させることもできる。 ただし、最終行のように変数が関係する場合は「+」を省略することはできない (この「+」を削除して実行し、エラーになることを確認してみるとよい)。
逆に文字列の一部を取り出す場合には、次の 3-4.txt のようにすればよい。
text = "This is a pen. What is that?\nIt is an apple. Oh, I see."
print (text[0])
print (text[1])
print (text[2])
print (text[3])
print (text[10:13])
print (text[:14])
print (text[15:])
print (text[-10:])
print (len (text[:14]))
「text[n]」とすれば、text に収められている文字列のうち n 番目の文字を取り出すことができる。 「text[10:13]」とすれば、10 番目の文字から 12 番目の文字までが取り出される。13 番目までではないことに注意を要する。 「text[:14]」や「text[15:]」のように、範囲の始まりや終わりを省略すれば、最初や最後を示したものと解釈される。 「text[-9:]」のように負値が示された場合には、後ろから数えられる。 また、文字列の長さを調べたい時には「len()」関数を使えばよい。
2 日目に入る前に、Cygwin の使い方について少し補足をしておこう。 たとえば現在 /cygdrive/c/ というディレクトリにいて、次に /cygdrive/c/Users/ に移動したい、と仮定しよう。 このとき「cd /cygdrive/c/Users/」と入力すれば目的は達されるのだが、長くて入力が面倒であれば 単に「cd Users/」でも構わない。「cd」コマンドは、続くディレクトリ名が「/」で始まっていない場合、現在位置を基準とする相対的な位置 (相対パス、と呼ばれる) としてディレクトリを解釈するからである。
さて、Python では、基本的な演算を行うことができる。次のようなスクリプトを 2-1.txt として作成し、実行してみよう。 コンピューターの世界では、慣習的に、乗算記号として「*」を、除算記号として「/」を使っていることに注意されたい。 また、英文では、乗算記号として、しばしば「x」が使用されるが、その記法は Python では採用されていない。
print ("3 + 4 = ", 3 + 4)
print ("2 - 5 = ", 2 - 5)
print ("8 x 7 = ", 8 * 7)
print ("3 / 6 = ", 3 / 6)
さて、諸君の意図した通りに Python は計算を行っただろうか？ 不審な点がなかったなら、次は、先ほどの 2-1.txt とは少し異なる、以下のようなスクリプトを 2-1e.txt として作成し、実行してみよう。
print ("3 + 4 = " 3 + 4)
print ("2 - 5 = " 2 - 5)
print ("8 x 7 = " 8 * 7)
print ("3 / 6 = " 3 / 6)
2-1.txt と異なり、2-1e.txt では「,」を省略している。この場合、Python は「Syntax Error」を起こして適切に計算を行わない。なぜか。
2-1.txt では、print() コマンドの中に「,」で区切って二つの項目を指定していた。すなわち「"3 + 4 = " と表示し、その次に 3 + 4 の計算結果を表示せよ」と、諸君は命令したわけである。 ところが 2-1e.txt では、「"3 + 4 = "」と「3 + 4」が漠然と並べられており、両者の関係が不明瞭であるから、命令内容を Python は理解できない。 「二つの別々の項目である」ということを明示しなかった諸君の指示が悪いのである。 これを「融通が利かない、不便である」と考えるべきではない。曖昧な指示をコンピューターが勝手に解釈して処理してしまうと、予期せぬ事故につながるからである。
さて、計算の続きを試してみよう。2-2.txt は、次のような内容で作成する。
print ("10 / 3 = ", 10 / 3)
print ("10 / 3 = ", 10 // 3, "remainder", 10 % 3)
一行目の実行結果からわかるように、Python における除算は、基本的には小数として実行される。 二行目には見慣れぬ演算記号が登場している。この「//」は、除算の商の整数部分、「%」は「余り」を計算するものである。つまり「10 割る 3 は 3 余り 1」というわけである。
もう一つ、基本的な演算として冪乗も触っておこう。これば「べきじょう」と読む。「2 の 3 乗」とか「5 の 4 乗」とかいう、アレのことである。 コンピューターの世界では「2 の 3 乗」を「2^3」などと表現することもあるが、Python ではその記法を採用しておらず、「2 ** 3」と書く。 2-3.txt を次のように作成し、実行しよう。
print ("2^3 = ", 2 ** 3)
print ("2^1.5 = ", 2 ** 1.5)
print ("3 * 2 ** 2 = ", 3 * 2 ** 2)
この計算結果からわかるように、冪乗の計算は、乗算や除算よりも先に処理される。
最後に変数を扱っておこう。これは数学でいう変数とは意味が異なるし、記法も数学的ではないので、人によってはキモチワルイと感じるかもしれない。 2-4.txt は次のように作る。
x = 3
y = x ** 2
x = 1
answer = y + x
print (answer)
Python における「=」は「代入する」という意味であって、数学でいう等号とは意味が異なる。 すなわち、2-4.txt の内容は「x に 3 を代入せよ」「y に x の 2 乗を代入せよ」「x に 1 を代入せよ」「answer に y + x を代入せよ」「answer の中身を表示せよ」という意味である。 ここで x, y, z, answer は変数である。変数の名前は、一文字でも構わないし、もっと長くても問題ない。 さて、このスクリプトの実行結果は、諸君が予期した通りのものだっただろうか？
もしこれが数学だったなら、answer の中身は y + x であって、y は x^2 であり、また x = 1 なのであるから、つまり y = 1 であって answer = 2 である、と考えたくなる。 しかし Python では「y = x ** 2」が実行された時点で、y には「9」が代入されているのであって、「x ** 2」自体が代入されているわけではない。 すなわち、後に x に 1 が代入されたからといって、その時 y も一緒に変わることはない。 従って、最終的には変数 answer に 10 が代入されている状態になる。
諸般の事情で今年度中に Python を習得する必要が生じたので、この機会に、プログラミング初心者向けの Python 入門書を書こうと思う。 類似の文書は既に世の中にたくさんあるのかもしれないが、知らぬ。私は、私が欲しいものを、自分で作るだけのことである。
Python というのは、スクリプト言語の一つである。 スクリプト言語というのは、あるプログラムに入力する命令を記述する言語体系のことである。 すなわち、Python なる言語を用いて、Python の実行プログラムに命令を与えることで、実に様々な演算処理を行うことができる。
Python がプログラミング言語であるか否か、という点については、議論の余地がある。 かつては、人間が読めるスクリプトをそのまま実行する Python や Perl などのスクリプト言語はプログラミング言語に含めず、 C 言語などの「人間には読めないプログラムに変更する、業界用語ではコンパイルと呼ばれる処理を行った後に実行するもの」のみをプログラミング言語と呼ぶ、という流儀もあったが、現在では圧倒的少数派のようである。 とはいえ、たとえば Excel のマクロはプログラミングに含まれるのか、という点は意見が分かれるであろう。あるいは最大限に広く解釈すれば、コンピューターへの入力は全て広義のプログラミングにあたる、と、いえなくもない。 結局「プログラミング」という語を厳密にどう定義するか、というのは、ほとんど信仰の問題に過ぎない。 コンピューターの世界には、こういった信仰の違いによる不毛な争いが多く、たとえば Windows と Mac と Linux と *BSD のいずれが優れた OS であるか、といった論争は何十年も前から続けられている。 こうした紛争はコンピューター界における伝統であり尊重すべきであるが、プログラミング初心者にとっては理解の及ばぬ世界であろう。 従って、ここでは Python がプログラミング言語であるか否かについては論じない。 「Python はスクリプト言語である」としておけば、そこに異論を唱える人はほとんどいないであろう。
私が現在職場で使用しているコンピューターの OS は、残念ながら Windows である。残念ながら、というのは、私はむしろ NetBSD の方が好きだから、である。 NetBSD は、いわゆる unix 系 OS の一つである。 Linux の親戚である、と思ってもよいが、「Linux は unix 系 OS であるか否か」というのは論争になりやすいデリケートな信仰上の問題なので、NetBSD と Linux を比較するのはやめた方がよい。 要するに、NetBSD というのはオタク向け OS の一つであると思っておけばよい。
NetBSD をはじめとする unix 系 OS は、初心者には扱いにくいが、熟練してくると実に使いやすい。 そこで本文書では、まず Windows 上で unix 系 OS に近い操作方法を実現するための手段として、Cygwin をインストールしよう。 これは https://cygwin.com で配布されている。 Cygwin をインストールする際に、どのパッケージをインストールするかを選ぶ場面がある。ここで Python を加えることを忘れてはならぬ。 「View」で「Full」を選び、「Search」に「python」と入力すれば、下のリストに様々なパッケージが表示されるので、その中から python を選び、「skip」となっているものを適切に変更すればよい。
インストールされた Cygwin を起動すると、何やら「コンピューターみたいな」画面が表示される。ここで「exit」と入力すれば、Cygwin を終了できる。 もう一度 Cygwin を起動し、今度は「cd C:」と入れてみる。すると、C ドライブに「移動」することができる。この「cd」というのは、場所を移動するためのコマンドである。 ドライブだけでなく、「cd C:/Users/」などと入力することで、好きなディレクトリに移動することもできる。ディレクトリというのは、Windows では「フォルダ」と呼ばれているアレのことである。 ディレクトリの区切りには「/」を使うのがよい。 好きなディレクトリに移動して「ls」と入力すれば、ディレクトリの内容をリストアップすることができる。 Cygwin にはいろいろなコマンドがあるが、詳しくは google 先生にでも尋ねるとよい。Unix 系 OS で使えるコマンドの大半は、Cygwin でも使える。
ひとしきり Cygwin で遊んだら、話を python に戻そう。ここでは、まず、次のような内容のファイルを作る。ファイル名は何でも構わないが、ここでは 1.txt としよう。 ファイルを作るのは、Windows のメモ帳でも構わないし、何か他の好きなソフトウェアを使ってもよい。
print ("Hello, World!")
Cygwin で、この 1.txt が置かれているディレクトリに移動し「python 1.txt」というコマンドを実行しよう。 これは「python」のプログラムに「1.txt」という文字列を渡して実行する、という意味である。 python の場合、渡された文字列をファイル名として認識し、そのファイルに書かれている命令を実行することになっている。 すなわち、python は「print ("Hello, World!")」という命令を実行するのである。
python において、「print ()」という命令は、括弧内の文字列を表示せよ、という意味である。 「Hello, World!」をダブルクォーテーションで括っているのは、「これを文字列として解釈せよ」という意味であるが、今のところは、そういうオマジナイだと思っておいてもよい。
さて、諸君が「python 1.txt」を実行したとき、予期された現象が起こったであろうか？ 予期せぬ結果が生じた場合、まず確実に、諸君が何かを間違えたのである。コンピューターは、指示された内容を、忠実に実行するのみだからである。
