科学技術計算用言語 Julia

Julia とは何か

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

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

1. 科学技術計算に便利な文法になっている
2. 科学技術計算に便利なコマンドやライブラリが最初からある
3. 計算が速い ¹
4. 言語設計が近代的なので，プログラミングが随分楽

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

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

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

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

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

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

1. Julia は(上に書いたように) JIT コンパイラがその場でプログラムソースをコンパイルして動かす仕組み．
   つまり，プログラムを「本当に」実行するマシン上には Julia 環境(Julia カーネル)がインストールされている必要がある．
   
2. プログラムやルーチンをそのカーネルセッションで「初めて動かすとき」は JIT コンパイルにかかる時間が余計にかかる．
   「ちょっと 1回だけ触ってみた」という使い方だとこのせいで「Julia って遅くね？」と誤解するかもしれないので注意しよう⁴.
   
3. Julia プログラムを動かす「カーネル」と，人間用インターフェイスは分離できるし，違うコンピュータに置くこともできる．
   よって「使い勝手の良いインターフェイス」を自分で選ぶことが推奨されるし，インターフェイスは手元の PC で動かしつつ，実際にプログラムが動く環境は計算用の大きなサーバに置くということもできる．
   
4. インターフェイスとしては，以下の様なものが有名．好きなものが選べるが…
   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 らしい．なかなか良さそうだ．
      
5. グラフを描くライブラリや文法はデフォルトでは決まってない．
   インターフェイスに自由度があるので，グラフをどういう手段で描くのがベストか，決めようがないので、まあ仕方ないか． どういうものがあるかはいずれ解説するとして，とりあえず「そういうもんだ」と思っておきましょう．

ともかく試してみる

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

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

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

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

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

もう少し試してみる

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

  実習

1. 超入門
   
2. 魔法の加速法

自分でやってみる

さて、次に、ごく簡単な問題に実際に 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$ を関数としてプログラムする．

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

素朴な実装例

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

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

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

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

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

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

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

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

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

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 になることを、以下のように確認しておこう．

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

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

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

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

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

2. 適当なグラフ描画パッケージ(PyPlot や Plots など)を使ってグラフを描く
   PyPlot の場合:

   using PyPlot
   plot(y,"o")   # "o" をつけると、点グラフになる
   

   
   
   Plot (+GR) の場合:

using Plots
gr()

plot(y, st = :scatter )   # series type を点描画で、という指定

すると、

という感じのグラフが得られる．

なにか規則性があるような無いような？

レポート

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

課題

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