SocioInfo#20 講演「フラクタルへの誘い」の講演資料とコーディング例をまとめました．
次の手順でご利用ください．
[資料一式へのリンク]←ここをクリックしてデータ一式をダウンロードしてください。$ cd 〇〇〇〇〇〇/Code [RET]
以下の操作をして，このディレクトリで python 3.x を使えるようにしてください（システム全体で python 3.x を使えるようになっている場合は不要です）．
$ pyenv local (インストール済みpython 3.x のバージョン) [RET]
これで準備OKです．
次のようにしてコードを試してみることができます．
$ python circle.py [RET]
マージソートです。併合ソートともいいます。マージソートの基本的な考え方を次図に示します。分割して、順序よく並べ替えることでソートを実現します。
配列xには、[2, 18, 1, 9, 7, 13, 5, 8, 4, 15]という数が並んでいます。これを、前半部と後半部に分割し、[2, 18, 1, 9, 7]と[13, 5, 8, 4, 15]に分けます 。
これを、なんらかの方法でソートします。ソートした結果は[1, 2, 7, 9, 18]と[4, 5, 8, 13, 15]になっているでしょう。その結果をマージしていきます。両者を先頭から順番に睨んでいき、小さいほうを選んで並べていきます。両方とも既に個別にはソートされているので、左から順番に見比べていけばよいということはわかるでしょう。両者をマージし終わった段階で、全体の並べ替えが終了します。
問題は、「なんらかの方法でソート」というところでした。ここで、再帰的な思考が役立ちます。すなわち、半分にしたものもマージソートで並べ替えてしまえばよいだろうという発想です。分割を繰り返していけば、最後は要素が1個になるはずです。そのときは、並べ替えの必要はありません。その1個をそのまま返してあげればよいのです。その処理が再帰の終端条件になります。
マージソートは、分割統治法（divide and conquer algorithm）の一例です。分割統治法とは、問題を部分問題に分けて考えることで容易に解き、それらを統合することで全体を解決しようという考え方です。分割統治法では多くの場合、再帰的に分割統治を行うことで、問題を簡単に解決することを考えます。分割統治法の考え方のマージソートへの応用は、すでに整列している複数個の列を1個の列にマージする際に、小さいものから先に新しい列に並べれば、新しい列も整列されているだろうという考え方に相当します。
なお、マージソートの計算量は、O (n log n)です。データを半分にする処理はO (log n)で実現されるということはこれまで説明してきたとおりで、それぞれの段階でマージする処理は先頭から順番に見ていくだけなのでO (n)です。したがって、全体ではO (n log n)の手数がかかるということがわかります。
Pythonでマージソートを実装してみます。プログラムは次のとおりです。バブルソートやビンソートと比べると、多少、複雑になりました。ただし、再帰的定義をしていることに注意さえすれば、理解はさほど難しくないと考えますが、いかがでしょうか。
def merge_sort(x):
retary = []
if len(x) <= 1:
retary.extend(x)
else:
m = len(x) // 2
first = merge_sort(x[:m])
second = merge_sort(x[m:])
while len(first) > 0 and len(second) > 0:
retary.append(first.pop(0) \
if first[0] < second[0] else second.pop(0))
retary.extend(first if len(first) > 0 else second)
return retary
解説を加えます。まず、返り値を戻すための配列を空の配列retaryとして用意しておきます。引数で与えられる変数xの長さが1以下であれば、その配列をextend()でretaryにくっつけてしまいましょう。なお、この処理はlen(x) == 1:としてもほぼ問題ありません。しかし、その条件にするとxに空配列を与えたときにエラーが発生します。そのエラーを回避するために1以下としました。
さて、else: 以降が本体です。まず、対象の配列を真ん中で分割したいので、真ん中のインデクスを求めています。演算子「//」はPython特有のもので、「a // b」はaをbで割った商の整数値をとります。配列xの長さが奇数であるときを配慮したものです。
次のx[:m]とx[m:]も説明が必要でしょう。これはスライスを用いて配列を分割する処理です。前者はxの先頭からm個目までの要素を持つ配列を作り、後者はm番目以降最後までの要素を持つ配列を作ります。したがって、この表記で前半と後半に分割していることになるのです。
再帰的定義でそれらをソートした結果を、配列変数firstとsecondに格納しています。続くwhile文は、それぞれを先頭から順番に眺めて、小さいほうをpop(0)で取り出してretaryに追加していく処理です。どちらかの要素がなくなるまで その処理を続けます。
最後に、残りの要素を追加しています。どちらかの要素がなくなった時点で、片方には要素が残っているはずですが、それは最後に追加した要素よりも大きな要素のみが整列した状態で残っているはずです。したがって、配列の最後にextend()で追加してあげればよいでしょう。while文での処理とこの最後の処理に関しては、三項演算子の「〜 if 条件式 else〜」という記述になっているので少しわかりにくいかもしれませんが、丁寧に追いかけてみればそれほど難しくはないでしょう。
実行例を図に示します。適切に並べ替えが行われていることを確認してください。
標準的によく利用されているソート方法でありながら、その原理をきちんと理解するのはなかなか難しいクイックソート、そのアルゴリズムを解説します。
クイックソートは分割統治法の考え方で進められます。クイックソートの考え方を次図に示します。
まず、ピボット（pivot）となる数値を決めます。これはどの要素でも構いません。ここでは、並べ替える配列の最後の要素としています。図では8をピボットとしています。
ピボットを決めたら、残りの要素を先頭と最後から中心に向かって調べていきます。両者を比較し、前者がピボット以上の値、後者がピボットより小さい値のときは、その値を交換します。前者がピボットより小さい値ときはそのまま後ろに進みます。後者がピボット以上の値のときもそのまま前に進みます。
探索が進み、両者がぶつかったところで全体の探索は終了です。このとき、探索がぶつかったところより先はピボットより小さい値ばかりとなっているはずで、後ろはピボット以上の値になっているはずだということはおわかりでしょう。これが理解できれば、前と後ろ、それぞれの集合をクイックソートの再帰呼び出しで処理すればよいということがわかります。
再帰の終端条件も簡単です。引数として与える配列の要素が1個ないしは空配列のときは、再帰呼び出しを行わずにそれのみからなる配列を返せばよいのです。図では最初のピボット以下の左側だけを図示していますが、全く同じ方法で右側も処理できます。
Pythonによるクイックソートの実装例を次に示します。少し、複雑なプログラムになってしまいました。しかし、前述したクイックソートのアルゴリズムを理解していれば、このコードを理解するのはそれほど難しくはないでしょう。
def quick_sort(x):
retary = []
if len(x) > 1:
y = x[:]; first = []; second = []
pivot = y.pop()
while len(y) > 0:
top = y.pop(0)
if top < pivot: first.append(top)
else: second.insert(0, top)
if len(y) == 0: break
last = y.pop()
if last < pivot: first.append(last)
else: second.insert(0, last)
retary = quick_sort(first) + [pivot] + quick_sort(second)
else:
retary.extend(x)
return retary
さて、ソートの処理ですが、アルゴリズムの説明で示したように配列の添字を工夫して要素を入れ替えて、という処理は細かな条件判定がややこしくなりそうなので、配列を複製して先頭と末尾から要素を取り出しつつ進めていくという方針で実装しています。取り出した値を比較して、ピボットの左側（first）と右側（second）の配列に値を追加していくという方法で、先に示した交換方式と同じ結果となるようにしています。
y ＝ x[:]というコードは、スライスを用いて配列を複製するものです。その後、空の配列firstとsecondを用意し、複製したyの末尾からピボットを取り出しています。
while文は、複製したyが空になるまで繰り返します。先頭から取り出し、ピボットより小さければfirstに前から詰めていきます。ピボットより大きければsecondに後ろから詰めていきます。append()した結果が前から詰めることになる、あるいは、insert(0)した結果 が後ろから詰めることになるという処理に注意してください。
配列xが偶数個の要素を持つときは、複製したyからピボットを取り出した結果、yの要素数は奇数個になってしまいます。このとき、pop(0)としてtopを得た時点で配列yは空になってしまうので、そのときは繰り返し処理の途中で抜けなければいけません 。その対応が次のif文とbreak文の組み合わせです。
問題なければ後ろからlastを得て、同様にfirstまたはsecondに追加していきます。
繰り返しを終了した時点で配列firstにはpivotよりも小さい値のみが、配列secondにはpivot以上の値のみが詰められているはずです。したがって、再帰呼び出しでそれぞれをソートしたものを用意して、pivotを挟んでそれらを順番に並べてやれば 、全体がソートされているはず、という手順です。
上図にクイックソートの実行例を示します。配列xの内容が昇順に並べ替えられていることがわかるでしょう。
オリンピック準備の狂想曲が原因か，はたまた梅雨があけたとたんのこの暴力的な暑さが原因か，なにやら漢字変換がおかしい．たとえばこれ．
「構造化」という言葉を打ち込みたくて「こうぞうか」とタイプし漢字変換したら「Googleの」と出てきた．上の図はいったん削除して，再度「こうぞうか」と打ち込んだ状態でスナップショットを撮ってみたものだ．
そして極め付きはこれ．
「けいたいそかいせき」で変換したら「く，なんというか……」と変換された．なにそれ？
しかもそれしか選べない．どういうことでしょうかw
思うところありGoogle Trendsで「Python」と「AI」というキーワードの関係を調べてみた．結果はこのとおり．
青がPythonで赤がAIである．地域は日本，期間はGoogle Trendsが提供している期間のなかで最大にして2004年から，ということにした．
いやはや見事に相関していることがわかる．両者とも，2015年くらいからの伸びが顕著だ．AIというキーワードでの検索が増えていると同時に，Pythonというキーワードによる検索も急激に増えている．
データをダウンロードしてその相関係数をサクッと計算してみたところ，0.9117という値が出た．作為的なデータじゃないのに0.9を超える相関係数って滅多にみないよ？（でも，皆無というわけでもなくて，今回のようにたまに出くわして驚くことがある）．
いろいろな言語のトレンド見てたら島根県の特異な地域性が浮かび上がってきた．これについては何も申し上げますまい．
以前，JavaScriptでヒルベルト曲線を描くプログラムを紹介した．今回は，Pythonのタートルグラフィクスライブラリを使って同様のヒルベルト曲線を描くプログラムである．
ヒルベルト曲線に限らずフラクタルは再帰の考え方で描くことができるが，まさに下記の関数がそれを表している．こんなに簡単に表現することができるのは素晴らしい．NumPyを利用して行列演算をたんなる掛け算として表すことができる点も，単純化に寄与している要因である．
def hilbert(n, m):
if n > 1:
for m_ in tm: hilbert(n-1, m * m_)
else:
for q in [ m * p_ for p_ in p ]: draw_line_to(q)
なお，一時的に用意している「m_」や「p_」といった変数は，とくに名前を付ける必要はないのでアンダースコア1文字で次のように表現することも許されているらしい．しかし，ここまでやるのは少しやり過ぎな気もする．
def hilbert(n, m):
if n > 1:
for _ in tm: hilbert(n-1, m * _)
else:
for q in [ m * _ for _ in p ]: draw_line_to(q)
タートルグラフィクスなので，これだけのプログラムでアニメーション動作する．動画を用意した．ずっと眺めていても飽きないゾ？
プログラムの全体は以下のとおりである．
#!/usr/bin/env python
from turtle import *
import numpy as np
SIZE = 300
MARGIN = 50
penup_flag = True
settings = [ { 'color': x[0], 'width': x[1] } for x in [
['forestgreen', 5], ['navy', 4], ['purple', 3], ['brown', 2],
['red', 2], ['orange', 1], ['yellowgreeen', 1] ] ]
tm = [ np.matrix(x) for x in [
[ [0.0, -0.5, -0.5], [-0.5, 0.0, 0.5], [0.0, 0.0, 1.0] ],
[ [0.5, 0.0, -0.5], [ 0.0, 0.5, -0.5], [0.0, 0.0, 1.0] ],
[ [0.5, 0.0, 0.5], [ 0.0, 0.5, -0.5], [0.0, 0.0, 1.0] ],
[ [0.0, 0.5, 0.5], [ 0.5, 0.0, 0.5], [0.0, 0.0, 1.0] ] ] ]
e = np.eye(3)
p = [ np.matrix(x).T for x in [
[-0.5, 0.5, 1.0], [-0.5, -0.5, 1.0],
[ 0.5, -0.5, 1.0], [ 0.5, 0.5, 1.0] ] ]
def draw_line_to(x):
goto(x[0,0]*SIZE, x[1,0]*SIZE)
global penup_flag
if penup_flag:
pendown()
penup_flag = False
def reset():
penup()
global penup_flag
penup_flag = True
def hilbert(n, m):
if n > 1:
for m_ in tm: hilbert(n-1, m * m_)
else:
for q in [ m * p_ for p_ in p ]: draw_line_to(q)
setup(width = 2*SIZE+MARGIN, height = 2*SIZE+MARGIN)
for i in range(len(settings)):
reset()
color(settings[i]['color'], 'white')
width(settings[i]['width'])
hilbert(i+1, e)
onkey(exit, 'q')
listen()
mainloop()
梅雨も終盤にさしかかり，雨が降り続いている．晴耕雨読，そんなときは部屋の中で仕事をするに限るのかな，なんて考えていたら「そういえばここのところよく『線状降水帯』って言葉聞くけど，そういう言い方するようになったのは最近よね」と尋ねられた．
はて，そうだっけか．そういうときは，Google Trends に聞いてみよう．
https://trends.google.com にアクセス，「線状降水帯」というキーワードを入れて，期間は直近5年間にしてみた．
面白い．毎年，梅雨時になると「線状降水帯」というキーワードで皆さんがググっていることが明白である．昨年にぐっと人気(?)が出たということもわかる．今年も多めかな？
ところでWikipediaで「線状降水帯」の項を紐解くと「日本でこの用語が頻繁に用いられるようになったのは、2014年の平成26年8月豪雨による広島市の土砂災害以降とみられる」との記述がある．本当にそうだろうか．
Google Trends の期間を2011年7月から2021年7月までの10年間にして再検索してみた．うん，やはり2015年くらいから検索対称に現れてきているようだということがわかる．
なお，このGoogleトレンド，季節性のある単語を入れてみると，あーやはり季節のイベントっていうのは年周期でやってくるんだなあということがよく分かって面白い．以下の図は，5年間で「バレンダインデー，お花見，お盆休み」というキーワードの状況を見たものだ．いかがでしょう？
来月発売予定になっている書籍の校正をしている．オンライン授業を猫が邪魔する問題について触れている箇所ではショータとフーちゃん2匹を登場させた．本ニャンたちの許可は得ていないが，印税でちゅーるでも買ってやろうか．
なお，次の図は私が2000年4月に出版した「Linuxによる画像処理プログラミング」という本にこっそり掲載した図版である．身内の画像を使う悪い癖はその頃から．当然のことながら，こちらも本人の許可は得ていない．
Google My Maps に久しぶりにアクセスしたら，なんだかずいぶん進化していたので驚いた．ジオコーディングを自動でやってくれるので，試してみると，けっこうビックリする．学生にやらせたら「わーすごい」と素直にビックリしてた．まあ，Googleマップで検索するとピンポイントで探してくれるので，それをまとめてやってくれているというだけではあるが……
（こちらもどうぞ→「可視化は楽しい」）
今回，手元に用意したのは次のようなデータである．井の頭線の各駅と渋谷駅からの距離の情報を，CSV（Comma Separated Values）形式で記録したものである．距離はレイルラボのデータを参考にした．1行目にラベルとしてstations, distanceという名前が付けられている．これを，stations.csvというファイルに保存しておこう．
さて，Google Drive にアクセスする．ところで，アカウントの設定によっては Google My Maps にアクセスが許可されていないことがあるので気を付けられたい．大学のアカウントで試してみたら，アクセスすることができなかった．授業で学生にやらせてみようといそいそと準備したのに，「先生！なんかエラーメッセージが出てアクセスできません」が頻発してしまった．んー，残念．
左上の「＋ 新規」ボタンから，Googleマイマップを選ぶ．日本のユーザであれば，日本の地図が現れるはずだ．
「無題の地図」の「無題のレイヤ」にデータを投入しよう．「インポート」をクリックする．
受け付けるファイルはCSVの他，エクセルのデータや地理データなどである．ここに，先のCSVファイルをドラッグ＆ドロップする．
目印を配置する列を聞いてくるので，ここではstationsを選ぶ．駅の名前が並ぶカラムである．駅名をヒントにジオコーディングが行われ，地図上の位置が定まる．「続行」ボタンで続けよう．
続いてマーカーのタイトルとして使用するカラムを選ぶ．ここも stations としておこう．「完了」ボタンを押せばインポート作業終了である．
＼(^o^)／ 駅にマークが付いた！
渋谷駅からの距離に応じて表示を変えてみよう．「均一スタイル」のリンクをクリックする．
均一スタイルのプルダウンメニューから，「データ列別のスタイル：」→「distance」を選ぼう．
「範囲」にチェックを入れる．
グラデーションの色を選べるので，虹色に変えてみよう．グラデーションのプルダウンメニューで，一番上の派手なやつを選ぶ．
さらに，4段階のグラデーションを，12段階に変更する．
以上で，できあがり！
