科学技術計算用言語 Julia

Julia とは何か

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

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

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

という、大変多くの利点がある言語になっている．これからが楽しみな言語だ．

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

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

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

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

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

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

それぞれの学習環境での状況

今回は Julia およびその web インターフェイス jupyter が使えないといけない．詳細は下記に記載したが，なるべく自分の PC/Mac を用いるのが良いだろう．

環境	状況
阪大情報教育システム Windows10	Julia 1.3.1 がインストール済み．ただしfirewall 下にあるため, 様々なファイルを取得するための `git` コマンドや `curl` コマンドがそのままだとインターネットと通信できないなど，Julia を使うには現状ではやや不向きな環境だ．どうしてもという場合は Julia 利用のために: 0. 準備編. ファイヤーウォール中の場合の設定と， Julia 利用のために: 2. jupyter インストール編の作業をしてから使おう．なるべく，今回はこれを用いるのはやめておこう．
阪大情報教育システム CentOS7	firewall の問題もあるし，かなり遅いのでおすすめできない．どうしてもという場合は，Julia 利用のために: 0. 準備編. ファイヤーウォール中の場合の設定相当の作業を行った後に，Julia 利用のために: 1. Julia インストール編と Julia 利用のために: 2. jupyter インストール編の作業を行えば良いはず．なるべく，今回はこれを用いるのはやめておこう．
各自PC CentOS7 obtained at OSBoxes	Julia 利用のために: 1. Julia インストール編と Julia 利用のために: 2. jupyter インストール編の作業を行えば良い．今回はこれを用いても良い.
各自PC Windows	Julia 利用のために: 1. Julia インストール編と Julia 利用のために: 2. jupyter インストール編の作業を行えば良い．今回はこれを用いても良い．
各自PC Mac OS	Julia 利用のために: 1. Julia インストール編と Julia 利用のために: 2. jupyter インストール編の作業を行えば良い．今回はこれを用いても良い．

Julia and jupyter をインストールして使えるようにするには

まあ，下記を順番に読んでいけば良い．多くの人は準備編は関係ないだろう．

Julia 利用のために: 0. 準備編. ファイヤーウォール中の場合の設定

Julia のパッケージと呼ばれるライブラリ環境のインストールと更新は内部的に git, curl といった unix のネットワーク系コマンドを用いて行われる．そのため，firewall の中(たとえば阪大情報教育システムはそうである)では，適切な設定をしないとパッケージをうまく扱えない．

そこで，たとえば阪大情報教育システム(windows)上の Julia を用いるのであれば，(Julia インストール直後ぐらいに)次のような設定をしておく必要がある．

1．git コマンドのための設定: Windows のメニュー → “Command Prompt” として command prompt を起動し，次の3つのコマンドを実行する．すると，とある設定ファイルが作られて下記の内容が書き込まれ, 以降のgit コマンドは proxy サーバを介して外へつながるようになる．

1
2
3


git config --global http.proxy http://toyonakans.ecs.cmc.osaka-u.ac.jp:3128
git config --global https.proxy http://toyonakans.ecs.cmc.osaka-u.ac.jp:3128
git config --global url."https://".insteadOf git://

2．curl コマンドのための設定:

2-a．おそらく c:\Users(もしくはユーザー)\アカウント名\.julia というディレクトリがあるので，その中に config というディレクトリを作る．

2-b．上で作ったディレクトリに startup.jl というファイルを作り，以下の2行を書き込んでおく．

1
2


ENV["HTTP_PROXY"] = "http://toyonakans.ecs.cmc.osaka-u.ac.jp:3128"
ENV["HTTPS_PROXY"] = "http://toyonakans.ecs.cmc.osaka-u.ac.jp:3128"

Julia 利用のために: 1. Julia インストール編

CentOS 7 に Julia をインストールする場合

下記のようにすれば最新バージョンの julia がインストールできるはずだ．

1
2
3
4


sudo yum update
sudo reboot  ← 念の為のリブート
sudo yum-config-manager --add-repo https://copr.fedorainfracloud.org/coprs/nalimilan/julia/repo/epel-7/nalimilan-julia-epel-7.repo
sudo yum install julia

注: 上の 3行目でなぜかエラーになることがある．エラーになった場合は，リブートしてやりなおしたりなど何回か試してみると良い．

通常の Windows に Julia をインストールする場合

1．Windows のアカウント名に
半角アルファベット, 半角数字, 半角ハイフン, 半角ピリオド, 半角アンダースコア
以外の文字を使っていないか確認する(これ以外の文字を使っていると上記の jupyter 環境が使えない)．

もしもこれらの文字をアカウント名の中で使っている人は，アカウント名をこの条件に沿う様に変えるか，この条件に沿ったアカウントを別に用意してそちらにログインし直し，下記の作業へ進もう．これらの文字をアカウント名の中に使っていない人は，そのまま次の作業へ進もう．

2．ver.1.4.2(for 64bit windows) をダウンロードしてダブルクリックしよう. 素の Julia がインストールできる．

通常の Mac に Julia をインストールする場合

Julia 1.4.2 for Mac 64bit をダウンロードし，help に従ってインストールしよう．

Julia 利用のために: 2. jupyter インストール編

jupyter インストールにはいくつかの方法があるが，ここでは最も簡単で，かつ，OS に依存しない方法を紹介しよう．次の手順で一回だけ作業を行えば良い．

1．まず，Julia を起動する．Julia がインストールされているならば，どの OS でもメニューのどこかに Julia があるはずなのでそれを使おう．すると，

といった感じの window が起動するだろう(これは julia の REPL インターフェイス)．ただし，window の色は環境依存なのでそこは気にしないように．

2．この Julia の window 中で ] キーを押す．すると，”julia>” というプロンプトが “(v1.4) pkg>” という感じのプロンプトに変わるだろう．これは Julia のパッケージ管理モードに入ったことを意味する．

3．このパッケージ管理モードのまま，

1

add IJulia

と入力する(大文字小文字の区別があるので注意)．これは jupyter を使うための “IJulia” というパッケージのインストール作業である．しばらくかかるので待とう．

4．上の作業が終わったら，backspace キーを押す．すると通常モードに戻り “julia>” というプロンプトに戻るはずだ．そうしてから，次の 2つのコマンドを Julia に実行させよう．

1
2


using IJulia
notebook()

すると上の2つ目のコマンド notebook() を実行した時に jupyter 関連のソフトウェア一式をインストールしてくれる．このときインストール方式等について尋ねてくるが，よくわからない場合は y キーを押しておけば良い．

しばらく待つと，作業が終わる．もしかしたら web browser が起動したりもするが，この段階ではいったんその web browser も閉じ，この julia 環境も終了させてしまおう．これで一応，Julia と jupyter が使えるようになったはずだ．

Julia & jupyter を起動/終了するには

起動の手順

Julia を jupyter インターフェイスで利用するには，基本的には

1．Julia を起動．
2．Julia の中で以下の2行を行う．

1
2


using IJulia
notebook()

だけで良い．こうするだけで web browser が呼び出され，そこに下のような jupyter インターフェイスが表れるはずだ(見えている細かい内容は環境によって異なる)．

これは jupyter のホーム画面とでも言うべき画面で，様々なファイル操作等はここを経由して行う．

次に Julia カーネルと接続して Julia を実際に使う画面を起動する．それには，この画面で右上に小さく見えている New ▼ というボタンを押して，

というメニューを出現させ，ここで “Julia 1.4.2” という部分を選べば良い．すると

といった画面が起動する．ここで入手力を行うことで，Julia を使う，Julia プログラムを動かす，といったことができる．

終了の手順

実際の使い方はこのあとすぐ説明するとして，次にここから終了する手順を示しておこう．

まず，Julia を使っている画面の左上を見て，

上図の赤い丸がついている部分のアイコンを押し，念の為に現状をセーブしておく．

次に，その上の “File” という部分を押すとメニューが出るので，その一番下の “Close and Halt” を選ぶ．

次に，jupyter のホーム画面に移り, 右上の “Quit” ボタンを押す．これで jupyter インターフェイスが終了する．あとは Julia 環境を終了させれば良い．

ともかく試してみる

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

実習
インストール不要, 無料で Julia をすぐ使える JuliaBox が公開されている
へ飛ぼう．書いてあることの半分ぐらいが今回の貴殿らの環境に合致するだろうから，そうした部分について、jupyter の操作の練習を兼ねてごく簡単に Julia を触ってこよう．
(注: 上の web で言及されている JuliaBox というサービスは現在は廃止されている)

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

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

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

もう少し試してみる

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

実習

自分でやってみる

さて、次に、ごく簡単な問題に実際に 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) を求めてグラフを描く

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

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

1

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

2．適当なグラフ描画パッケージ(デフォルトは 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)$ をおおよその数式にしてみよう！(学生のみんなへの課題だね)

レポート

以下の課題について能う限り賢明な調査と考察を行い，
2020-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]