科学技術計算用言語 Julia

Julia とは何か

ここで紹介する Julia とは、MIT とか NYU あたりの人たちが開発中の，2012年に公開が始まった，科学技術計算の専用言語のこと．本家 web にマニュアルやらチュートリアルやら動画やら、大変多くの様々な資料があるので、困ったらまずはここをみると良い．

そしてその言語の能力であるが、開発者らが欲張りだと自認するだけあって，後発の強みも活かして，

科学技術計算に便利な文法になっている
科学技術計算に便利なコマンドやライブラリが最初からある
*計算が速い
注: 実行時間は(JIT コンパイルにかかる時間は除いて)ほんとに高速で，C や Fortran とほぼ同等．julia の web の top にも表で載っているので見てみよう．
言語設計が近代的なので，プログラミングが随分楽

という、大変多くの利点がある言語になっている．開発版ですらまだバージョン 0.5 という未完成な状態ながら，これからが楽しみな言語だ．

今のところ？JIT compile, つまり「実行時にコンパイルして動かすスタイル」をとっているので，使い勝手的にはインタプリタと言えるプログラム開発中や，いろいろ試しながら動かすにも，これは便利．

よって，「たいへん気楽に使える高速な科学技術計算言語」として使える．また，ライセンスはMIT ライセンスなので安心して無料で使えるし，ソースも全部公開されているので中身をいじりたければ自由にいじれる．

ぜひ、この言語を自分のものとして、便利に使えるようになってもらいたい．

知っておくと混乱しないで済む？基本知識

少し、大雑把なことを書いておこう．

Julia は(上に書いたように) JIT コンパイラがその場でプログラムソースをコンパイルして動かす仕組み．
つまり，プログラムを「本当に」実行するマシン上には Julia 環境(Julia カーネル)がインストールされている必要がある．
プログラムやルーチンをそのカーネルセッションで「初めて動かすとき」は JIT コンパイルにかかる時間が余計にかかる．
「ちょっと 1回だけ触ってみた」という使い方だとこのせいで「Julia って遅くね？」と誤解するかもしれないので注意しよう.
Julia プログラムを動かす「カーネル」と，人間用インターフェイスは分離できるし，違うコンピュータに置くこともできる．
よって「使い勝手の良いインターフェイス」を自分で選ぶことが推奨されるし，インターフェイスは手元の PC で動かしつつ，実際にプログラムが動く環境は計算用の大きなサーバに置くということもできる．
インターフェイスとしては，以下の様なものが有名．好きなものが選べるが…
1. jupyter
  web browser で操作するインターフェイス．わかりやすいし， おすすめ ．
  これを使うと，Mathematica や Matlab のような見た目と使い勝手になる．
  後で紹介する JuliaBox もこれを使っている．
2. juno
  Julia 専用の統合開発環境．文法に沿ったいろいろなアシスト機能があって大変便利．ただ，うまく動かないときもあって，そうした場合のトラブル解決はやや面倒かも．
3. 文字端末ソフトで直接対話する REPL
  Windows だと Powershell で，Unix 系だと shell から Julia を立ち上げて，いわゆる「黒い画面」で操作するもので，マニュアルなんかでは read-eval-print loop の略で REPL と呼ばれている．
  悪くない．悪くないが，カーネルが落ちるとそれまでのやりとりも画面から消えてしまうので，長時間扱うならやっぱインターフェイスは別にしておくといいかも．
  ちなみに，REPL でもグラフや画像の表示は思ったよりきちんとできる．画像がぽこっと別ウィンドウや web browser 呼び出しで出てくる様子はちょっとかわいいので一度試してみるといいかも．
4. ESS (Emacs Speaks Statistics)
  Emacs で動く統合開発環境のようだ．
  ちなみに，Emacs の julia-mode 用の elisp もあって，syntax hilighting についての解説もあったりする．
5. Julia-vim
  vim の julia-mode らしい．なかなか良さそうだ．
グラフを描くライブラリや文法は自由に選べる．まあ，緩いデフォルト的な存在(Plots package)はあるけど．

ともかく試してみる

文法云々の堅い話はおいておいて、まあ、まずは触ってみよう．

実習
インストール不要, 無料で Julia をすぐ使える JuliaBox が公開されている
へ飛んで、jupyter の操作の練習を兼ねてごく簡単に Julia を触ってこよう．

文法等のマニュアル、解説

日本語で書かれた解説が載っている M. Hiroi’s Home Page: Julia Language Programming や、英語でよければ本家のマニュアルあたりをみると詳細な文法の解説があるので、困ったらこれらを読もう．

また、できれば本を一冊買っておくと良い．おすすめは書籍: Julia 1.0 Programming で、英語だけれども大変読み易い．

もう少し試してみる

実際に自分でプログラムを組む前に、もう少し試してみたほうが良いだろう．というわけで、下記の実習をやっておこう．

実習

自分でやってみる

さて、次に、ごく簡単な問題に実際に Julia で取り組んでみるのが良いだろう．というわけで、次の

Collatz の予想(角谷の問題)

正の整数 $n$ に対して、次のような写像 $f$ を考える．

\begin{equation} f(n) = \left\{ \begin{array}{lcl} 3n + 1 & \mbox{　　　　　} & : n \mbox{ が奇数の場合} \newline n/2 & & : n \mbox{ が偶数の場合} \end{array} \right. \end{equation} このとき、任意の正の整数に対して、$f$ を複数回適用するといずれ必ず結果が 1 になる.

という、現在なお未解決の数学の予想がある．この問題を対象にして少しいろいろ練習してみよう．

まず最初は例を示そう．

上の写像 $f$ を関数としてプログラムする．

やり方はいろいろだ．例を示しておこう．

素朴な実装例

1
2
3
4
5
6
7


function f(n)
  if isodd(n)       # isodd は、中身が奇数かどうかを判定する命令．
    return 3n+1     # Julia は定数係数と変数の間の積の記号を省略できる．
  else
    return div(n,2) # n/2 の結果は浮動小数扱いになってしまうので、整数を返す div を使う．
  end
end

三項演算子を使ってシンプルに書いた例

1
2
3


function f(n)
  isodd(n) ? 3n+1 : div(n,2)
end

関数宣言もシンプルに書いた例

1

f(n) = isodd(n) ? 3n+1 : div(n,2)

実習
自力で、もしくは上のように $f$ を実装して、$f(3)$ が 10 になること、$f(10)$ が 5 になることを、以下のように確認しておこう．

1
2


f(3)  # 結果は 10 のはず．
f(10) # 結果は 5 のはず．

$f$ を何回適用すると結果が 1 になるかをカウントする関数 count を作る

例を一つだけ示しておこう．しかし, なるべく例を見ないで、マニュアルを読むなどして自力で組んでみよう．

1
2
3
4
5
6
7
8
9


function count(n)
    m = n
    number = 0
    while m > 1
        m = f(m)
        number += 1
    end
    return number
end

実習
自力で、もしくは上のように count を実装して、count(3) が 7 になること、count(17) が 12 になることを、以下のように確認しておこう．

1
2


count(3)  # 結果は 7 のはず．
count(17) # 結果は 12 のはず．

奇数 n = 3,5,7,… に対して、count(n) を求めてグラフを描く

実習
これも例を示しておくが、なるべく自力で頑張ってみよう．

まず、count(n) の結果の配列を作る

1

y = [count(n) for n in 3:2:2001] # 1000個ほどデータを用意する.

適当なグラフ描画パッケージ(デフォルトは Plots だ)を使ってグラフを描く.まあ，ここでは Plots を使ってみよう．

1
2
3
4


using Plots # セッション中で一回宣言すれば良い．

plot(y, st = :scatter, leg = false )   
# st = :scatter は点描画，leg = false は legend表示なし，の指定．

すると、
collatz points
という感じの，なにか規則性があるような無いようなグラフが見て取れる．この規則性，数式にできる？

n を大きくしていくと，あたりまえだけど「大雑把には」count(n) の最大値が大きくなる傾向があるね

$i = 3,5,7,…,n$ での count(i) の最大値を v(n) とすると(式だと次のようになるね)

\begin{equation} v(n) = \max_{i \leq n} \{ \mbox{count}(i) \} \end{equation}

この挙動はどういうものだろうか．

すぐに数式で求められそうにはないので，プログラムで求めてグラフを描いてみるのが良さそうだ．

例えば次のような感じにできるだろう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18


N = 200001

n_sq = [ 3 ]

max_c = count(3)
v_sq = [ max_c ]

for i in 5:2:N
    global max_c, n_sq, v_sq
    
    v = count(i)
    if v > max_c
        max_c = v
        push!(n_sq, i)
        push!(v_sq, v)
    end
    
end

とすると，n_sq に最大値 $v(n)$ が更新されたときの $n$ のリストが，そして v_sq にはその $v(n)$ のリストが計算されて記録される．だから，グラフを描いてみるとこんな感じ．

1
2


plot(n_sq, v_sq, marker = :auto, leg = false)
# marker = :auto はマークを自動的に決める指定．

う～む，どうも $\log$ 関数に似ているような気もするので，適当に並べて描いてみて比べてみよう．

1
2
3


plot(n_sq, v_sq, marker = :auto, leg = false)
plot!(n_sq, 30 .* log.(n_sq) )
# plot! はグラフを重ね描く．

う～むむむむ… 似ているような似ていないような… もっと N を大きくすると様子がわかるかな？

N を大きくして，x 軸の対数をとってグラフにしてみよう．

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20


N = 10000001

n_sq = [ 3 ]

max_c = count(3)
v_sq = [ max_c ]

for i in 5:2:N
    global max_c, n_sq, v_sq
    
    v = count(i)
    if v > max_c
        max_c = v
        push!(n_sq, i)
        push!(v_sq, v)
    end
    
end

plot(n_sq, v_sq, marker = :auto, leg = false, xaxis = :log)

う～むむむむむむむむ… 直線に近いがちょっとだけ強い…？

両対数だとどうだろうか．

1

plot(n_sq, v_sq, marker = :auto, leg = false, xaxis = :log, yaxis = :log)

おお！なにか掴んだ気がするぞ！ $n > 10^2$ ぐらいでの $v(n)$ をおおよその数式にしてみよう！(学生のみんなへの課題だね)

レポート

以下の課題について能う限り賢明な調査と考察を行い，
AppliedMath7-Report-13
という題名をつけて e-mail にて教官宛にレポートとして提出せよ. なお，レポートを e-mail の代わりに TeX で作成した書面にて提出してもよい.

課題

上の内容を、一回自力でやり通してみよう．
「正 $n$ 角形の周長 $L(n)$ を定規で測ったデータ」をプログラム計算で擬似的に生成して¹ 、$n \rightarrow \infty$ で確かに $L(n) \rightarrow 2\pi r$ (ただし$r$ は正 $n$ 角形の外接円の半径相当)となることを数値データで確かめよ．
問2 のデータに対して、Aitken 加速法を適用したらどうなるか試して、その状況を報告しよう．さらに、多重に適用してみよう．

三角関数などを使って良いということ．もっとはっきり書くと、$L(n) = 2rn\sin(\pi/n)$ としろ、ということになるw. ^[return]