いまから使えるJulia言語

Julia in Physics 2021 Online

本チュートリアルの目的

  • 「時代はJuliaかもしれない」という気にさせる
    • 速さだけがJuliaの良いところじゃない
    • 科学技術計算のためだけの言語じゃない
  • 全然Juliaを知らない人の「ちょっと試してみよう」のハードルを下げる
    • 事前アンケートによると約半数がJuliaを使ったことがない
    • まずはJuliaがどんな言語かをざっと紹介
    • 次にJuliaの構文や主な機能を解説
    • 最後に実用的なTipsを見る
  • 分量が多いので全て詳しくは説明できない
    • あとで参考にできるよう多めに書いてある
    • 詳しくは参考資料を参照
    • この資料は公開されます

宣伝

  • 1から始める Juliaプログラミング』発売中
  • 200ページ弱にJuliaのエッセンスをまとめた入門書
  • 科学技術計算を高速かつ手軽に行いたい学生や研究者向け
  • おかげさまでAmazonレビューも高評価
  • 出版社ウェブサイト

Juliaとは

現在までのJulia

  • 2009年から開発されている動的プログラミング言語
    • 2018年にバージョン1.0がリリース
    • 2021年9月現在の最新版はバージョン1.6.2
    • もうすぐバージョン1.7のリリース?
  • ここ数年で人気が急上昇

スクリプト言語として

  • 事前コンパイルの必要なし
  • julia script.jlで即実行
  • juliaコマンドで対話的実行
  • Jupyterもサポート
$ cat script.jl 
println("hello, world")
$ julia script.jl 
hello, world
$ julia -q
julia> A = randn(3, 3)
3×3 Matrix{Float64}:
 -0.378213   1.17271  0.646551
 -0.318493  -0.25825  0.121042
  0.822649   1.50609  1.01108

julia> A'A  # 行列積
3×3 Matrix{Float64}:
 0.921234  0.877702  0.548677
 0.877702  3.71027   2.24974
 0.548677  2.24974   1.45496

julia> 

特徴

  • 現代的な機能
    • 自動メモリ管理
    • ジェネリクス
    • メタプログラミング(マクロ)
  • MATLAB風の構文
    • endまでがブロック
    • 2xA'Bなどの数学記法
    • σなどユニコード変数をサポート
  • 高性能な実行時コンパイラ

# 1次元配列(Vector)は要素の型付き
typeof([1, 2, 3])        #> Vector{Int64}
typeof([1.0, 2.0, 3.0])  #> Vector{Float64}

# マクロによるコードの書き換え
x = 1
@assert x > 0
@assert x < 0  #! ERROR: AssertionError: x < 0

# if … end のブロック
if x > 0
    println("x is positive")
end

# ユニコード変数
μ, σ = -0.3, 1.2
x = μ + randn() * σ  # 正規分布 N(μ, σ^2)

クイックソートの実装例

  • どんなデータでも動く(generics)
    • 整数
    • 浮動小数点数
    • 文字列
    • etc.
  • データ型に応じて実行時コンパイル
    • 要素の型は入力から推論
    • isless関数が要素の型に応じて特殊化
  • マクロを通じた最適化ヒント
    • @inboundsは境界チェックを省略

quicksort(xs) = quicksort!(copy(xs))
quicksort!(xs) = quicksort!(xs, 1, length(xs))

function quicksort!(xs, lo, hi)
    if lo < hi
        p = partition!(xs, lo, hi)
        quicksort!(xs, lo, p - 1)
        quicksort!(xs, p + 1, hi)
    end
    return xs
end

function partition!(xs, lo, hi)
    pivot = div(lo + hi, 2)
    pvalue = xs[pivot]
    xs[pivot], xs[hi] = xs[hi], xs[pivot]
    j = lo
    @inbounds for i in lo:hi-1
        if isless(xs[i], pvalue)
            xs[i], xs[j] = xs[j], xs[i]
            j += 1
        end
    end
    xs[j], xs[hi] = xs[hi], xs[j]
    return j
end
                    

高機能な標準ライブラリ

  • 豊富な数値型
    • {8, 16, 32, 64, 128}bit ✕ {符号付き, 符号なし}の整数
    • 16, 32, 64bitの浮動小数点数
    • 任意精度整数と浮動小数点数
  • 線形代数
    • 多次元配列Array{T, n}
    • LinearAlgebra.jl: 内積・行列積・線形方程式・行列分解など
  • Pkg.jl
    • 高機能な標準パッケージマネジャ
    • 仮想環境の構築
    • Cライブラリのバイナリ管理

グルー言語として

  • 外部コマンドの実行
    • バッククォートでコマンドを作成
    • runeachlineで実行
  • C言語の呼出し
    • ccallでCの関数を実行
    • 構造体型なども扱える
  • PyCall.jlパッケージ

# 外部コマンドの実行
for path in eachline(`find . -name '*.txt'`)
    # gzip <"$path" >"$path.gz"
    run(pipeline(path, `gzip`, "$(path).gz"))
end

# C関数の呼出し
ccall(:sin, Cdouble, (Cdouble,), 1.0)

# Pythonの使用
using PyCall
optimize = pyimport("scipy.optimize")
optimize.brentq(x -> x^2 - 2, 0, 2)  # x^2 = 2 の解

ベンチマーク:数独ソルバー

  • 数独ソルバーでベンチマーク
    • アルゴリズムは単純なバックトラッキング
    • 簡単な問題と難しい問題の2種類で計測
    • ソースコード
  • 難しい問題ではJuliaが圧倒的に高速
    • easy.txt ⇨ 5倍遅い
    • hard.txt ⇨ 85倍速い
    • Juliaは起動は遅いが速度は速い
数独ベンチマーク結果
JuliaPython
easy.txt257.3 ㍉秒 52.4 ㍉秒
hard.txt 5.1  秒438.2  秒

# Julia版の数独ソルバー
function solve(puzzle)
    next(i, j) = i < 9 ? (i + 1, j) : (1, j + 1)
    function solved(sol, i, j)
        if j > 9
            return true
        elseif sol[i,j] > 0
            return solved(sol, next(i, j)...)
        end
        m, n = (i - 1) ÷ 3, (j - 1) ÷ 3
        for k in 1:9
            if all(sol[i,j] != k for i in 1:9) &&
               all(sol[i,j] != k for j in 1:9) &&
               all(sol[i,j] != k for j in 3n+1:3n+3,
                                     i in 3m+1:3m+3)
                sol[i,j] = k
                solved(sol, next(i, j)...) && return true
            end
        end
        sol[i,j] = 0
        return false
    end
    sol = copy(puzzle)
    return solved(sol, 1, 1) ? sol : nothing
end

セットアップ

Juliaのインストール

  • Windows, macOS場合
    1. 公式サイトからインストーラをダウンロード
    2. ダウンロードしたインストーラを実行
    3. 環境変数PATHjuliaコマンドのあるディレクトリを加える
  • Linuxの場合
    1. 公式サイトからバイナリをダウンロード
    2. ダウンロードしたTARボールを展開し、ファイルを適当な場所(~/.localなど)にコピー
    3. 環境変数PATHjuliaコマンドのあるディレクトリを加える

どのJuliaを入れる?

  • 公式がリリースするバイナリはいくつかある
    • Current stable release
    • Long-term support (LTS)
    • Upcoming release
    • Nightly build
  • 基本はCurrent stable release
    • 最新かつ安定版
    • 多くのパッケージの最新版もこれはサポートしている
    • LTSは安定しているが、機能が古くなりがち
    • Upcoming releaseは新機能のお試し
  • なるべく公式のバイナリがオススメ
    • トラブルが少なく、公式がサポートしている
    • DNFやAPTは公式バイナリでない
    • Homebrewは公式バイナリと同じ

REPLの使い方

  • juliaコマンドでREPLを起動
    • Ctrl+Dまたはexit()で終了
    • Ctrl+P/Ctrl+Nで前後のコマンド表示
  • 主なREPLのモード
    • julia> Juliaモード(Delキー)
    • help?> ヘルプモード(?キー)
    • shell> シェルモード(;キー)
    • pkg> パッケージ管理モード(]キー)

基本構文

数値リテラル

  • 小数点があれば浮動小数点数
    • 整数は64bit符号付きがデフォルト
    • 浮動小数点数は倍精度がデフォルト
  • 真偽値は二値論理で使用
  • nothingは意味のない値
    • println関数の返り値など
    • PythonのNoneに相当
    • missingは「値がない」ことを表す値

# 数値
42           # 整数
3.14         # 浮動小数点数(小数表記)
6.02e+23     # 浮動小数点数(指数表記)
1.2 - 0.9im  # 複素数
2//3         # 有理数

# 真偽値
true
false

# Nothing
nothing

# 欠損値
missing
                    

配列 1/2

  • [ … ]で配列作成
    • 1次元配列は列ベクトル
    • 2次元配列は行列
  • 添字は1始まりが標準
    • begin: 最初の要素
    • end: 最後の要素
配列操作関数
ndims次元
sizeサイズ
length要素数
zerosゼロ初期化配列
copy複製配列
sum総和
minimum最小値
vcat縦結合
hcat横結合
fill!全要素設定
push!要素追加
pop!要素削除
sort!整列

# 1次元配列(ベクトル)
[1, 2, 3]

# 2次元配列(行列)
[1.0 2.0 3.0
 0.0 4.0 5.0
 0.0 0.0 6.0]

# 要素の参照
x = [10, 20, 30]
x[1]      #> 10
x[3]      #> 30
x[begin]  #> 10
x[end]    #> 30
x[end-1]  #> 20

配列 2/2

  • 配列は要素の型をパラメタとして持つ
    • Vector{T}Tがパラメタ
    • Vector{T}Array{T, 1}の別名
    • Matrix{T}Array{T, 2}の別名
  • 配列要素の型は値から推定される
    • 空の配列のときは要素の型を指定可能
    • 空の配列ではデフォルトでAny

# 配列は要素の型を持つ
typeof([0])    #> Vector{Int64}
typeof([0.0])  #> Vector{Float64}

# 要素の型を返すeltype関数
eltype([0])    #> Int64
eltype([0.0])  #> Float64

# 空の配列
Float64[]   # 浮動小数点数の1次元配列
[]          # Any[] と同じ

タプル

  • 複数の値を組にして扱う
    • 値の型は任意
    • 個数も任意
  • 2種類のタプル
    • 名前無しタプル
    • 名前付きタプル
  • タプルは不変 vs 配列は可変
    • タプルの要素は入れ替えられない
    • タプルの長さは変えられない
    • 制約が強い代わりに軽量 ⇨ 関数から複数の値を返すときに使われる

# 名前無しタプル
t = ("Albert Einstein", 1879)
t[1]  #> "Albert Einstein"
t[2]  #> 1879

# 名前付きタプル
t = (name = "Albert Einstein", birth = 1879)
t.name   #> "Albert Einstein"
t.birth  #> 1879
t[1]     #> "Albert Einstein"
t[2]     #> 1879

# 最大値とその位置を返すfindmax関数
findmax([3.9, 5.2, 2.1, 4.0])  #> (5.2, 2)

基本演算

  • 基本四則演算はPythonなどと同様
    • /は浮動小数点数を返すので注意
    • 整数の商と剰余は÷%を使う
    • それぞれdivremと同じ
  • 論理演算
    • !: 否定
    • |: 論理和
    • &: 論理積
  • 基本的な演算子も関数
    • 2 + 3+(2, 3)も同じ
    • 演算子もオブジェクトとして受け渡し可能

# 四則演算
2 + 3  #> 5
2 - 3  #> -1
2 * 3  #> 6
2 / 3  #> 0.6666666666666666

# 商と剰余
9 ÷ 2  #> 4  = div(9, 2)
9 % 2  #> 1  = rem(9, 2)

# 論理演算
!true         #> false
true | false  #> true
true & false  #> false

# 畳み込み演算
foldl(+, [1, 2, 3], init = 0)  #> 6

分岐

  • 条件分岐はif … elseif … else
  • 条件式はBool型のみ
    • falsetrueのみが有効
    • 0""は使えないので注意
  • 2つの短絡評価
    • x && yは「xが真 iff yを評価」
    • x || yは「xが偽 iff yを評価」

# 条件分岐
if x == 0
    println("x is zero.")
elseif x > 0
    println("x is positive.")
else
    println("x is negative.")
end

# 条件演算子
println(x < 0 ? "negative" : "nonnegative")

# 短絡評価
x > 0 && x < 10  # 0 < x < 10 と同じ
x < 0 || x > 10
                    

反復

  • start:stopは範囲オブジェクト
    • start:step:stopで間隔を変更できる
    • stepを負の値にすれば逆順
  • forwhileブロック内ではbreakcontinueが使える
    • break: ループを抜ける
    • continue: ループの残りを飛ばす

# 範囲
1:9    # 1, 2, …, 9
1:2:9  # 1, 3, …, 9

# for文
for i in 1:9
    println(i)
end

# 2重ループ
for j in 1:9, i in 1:9
    println(i, j)
end

# while文
while true
    println("yes")
end

関数定義

  • 2種類の定義方法
    • 数学記法に近い定義(1行定義)
    • functionを使う定義
  • returnは省略可能
    • 1行定義だと省略する人が多い
  • 引数に型の制約をかけられる
    • 後述する多重ディスパッチと関係
  • x -> x + 1と書くと無名関数

# 1行定義
iseven(n) = n % 2 == 0

# 標準的な定義
function collatz(n::Integer)
    steps = 0
    while n != 1
        if iseven(n)
            n = n ÷ 2
        else
            n = 3n + 1
        end
        steps += 1
    end
    return steps
end

構造体型

  • 構造体型は複数のフィールドをまとめるデータ型
    • フィールド数は0個以上
    • フィールドには型を指定可能
  • デフォルトの構造体は不変
    • mutable structに変えると可変な構造体
    • 通常は不変な構造体を推奨
  • コンストラクタは定義しなくてもOK ⇨ デフォルトコンストラクタ

# 2次元極座標の構造体型
struct PolarCoords
    # フィールド
    r::Float64  # radius
    φ::Float64  # angle

    # ユーザ定義コンストラクタ(オブジェクトを作る関数)
    function PolarCoords(r, φ)
        r ≥ 0 || error("radius must be non-negative")
        0 ≤ φ < 2π || error("angle must be in [0, 2π)")
        return new(r, φ)
    end
end

# PolarCoordsオブジェクトの作成
p = PolarCoords(1.0, 2π/3)
p.r, p.φ  #> (1.0, 2.0943951023931953)

パラメトリック構造体型

  • フィールドの型はパラメータ化可能
    • 右の例のTがパラメータ
    • Tはどんな型にもなれる
  • 性能を損なわずデータ型の柔軟性が向上
    • 複素数型:Complex{T}
    • 配列型:Array{T, n}
    • 辞書型:Dict{K, V}
    • など

# 2次元極座標の構造体型(パラメトリック版)
struct PolarCoords{T}
    # フィールド
    r::T        # radius
    φ::T        # angle

    # ユーザ定義コンストラクタ(オブジェクトを作る関数)
    function PolarCoords(r::T, φ::T) where T
        r ≥ 0 || error("radius must be non-negative")
        0 ≤ φ < 2π || error("angle must be in [0, 2π)")
        return new{T}(r, φ)
    end
    PolarCoords(r, φ) = PolarCoords(promote(r, φ)...)
end

# PolarCoordsオブジェクトの作成
p = PolarCoords(1//2, 2//3)

基本機能

型システム

抽象型(abstract type)
他の型をグルーピングする型
具象型(concrete type)
実際にオブジェクトを作る型
  • 抽象型は下位の型を持てる
  • 抽象型と具象型で木構造を取る
    • 内部節が抽象型
    • 葉が具象型
    • 根はAny抽象型

Integer  # 抽象型
Int64    # 具象型
UInt8    # 具象型
Float64  # 具象型

# 階層関係
Int64   <: Integer  #> true
UInt8   <: Integer  #> true
Float64 <: Integer  #> false

42  isa Int64    #> true
42  isa Integer  #> true
0.5 isa Float64  #> true
0.5 isa Integer  #> false
                    

例:Number抽象型以下の階層

角丸の節が抽象型

多重ディスパッチ

  • 関数は複数のメソッドを持つ
    • 右の関数fには3つのメソッドがある
    • Pythonでは「メソッド → クラス」なのに対し、Juliaでは「メソッド → 関数」の関係
  • 引数の型に応じて実行されるメソッドが変わるのが多重ディスパッチ
    • ユビキタスな機能(*関数さえも)
    • 多重とは、すべての引数がメソッド決定に関与するという意味

# 関数fに3つのメソッドを定義
f(x::Int64)   = "Int64"
f(x::Float64) = "Float64"
f(x::Integer) = "Integer"

# 関数呼出し
f(42)    #> "Int64"
f(0.5)   #> "Float64"
f(true)  #> "Integer"
f('a')   #! MethodError

# Int * Intメソッドが呼び出される
2 * 3  #> 6

# Int * Vector{Int}メソッドが呼び出される
2 * [4, 5, 6]  #> [8, 10, 12]

多重ディスパッチの実用例

  • 極座標に演算を定義してみる
    • 新しい関数を定義するときは普通に定義
    • 標準ライブラリBaseの関数を拡張するときはBase.{関数名}で定義
  • 既存の関数をどんどん拡張してよい
    • ただし、意味をあまり変えすぎないように
    • type piracyもダメ

# 極座標の回転(新しい関数の定義)
rotate(p::PolarCoords, θ::Real) =
    PolarCoords(p.r, mod(p.φ + θ, 2π))

# 極座標のスカラ倍(関数の拡張)
Base.:*(α::Real, p::PolarCoords) =
    PolarCoords(α * p.r, p.φ)
Base.:*(p::PolarCoords, α::Real) = α * p

# 極座標の和(練習問題)
Base.:+(p1::PolarCoords,
        p2::PolarCoords) = …

多重ディスパッチの応用例

問題:fit関数を新しいデータ型MyDataに拡張したい

Pythonの場合


# Python
class Model:
    def fit(self, data):
        …

fit関数の実装を書き換える


def fit(self, data):
    if isinstance(data, MyData):
        …

Juliaの場合


# Julia
struct Model … end

function fit(mod::Model, data) … end

fit関数に新しいメソッドを加える


function fit(mod::Model, data::MyData)
    …
end

パッケージ管理

Pkg.jl

  • 標準のパッケージ管理ツール
    • パッケージの追加・更新・削除など
    • ~/.julia以下にファイルを配置
  • ]キーを押下するとパッケージ管理モードに移行
    • status: 状態を表示
    • add: パッケージを追加
    • remove: パッケージを削除

(@v1.6) pkg> status
      Status `~/.julia/environments/v1.6/Project.toml`
  [6e4b80f9] BenchmarkTools v1.1.1
  [944b1d66] CodecZlib v0.7.0
  [31c24e10] Distributions v0.25.11
  …

パッケージ管理ファイル

ペアで使われる2つのパッケージ管理ファイル:

プロジェクトファイルProject.toml
プロジェクトのメタ情報(バージョンや依存パッケージ等)を管理
マニフェストファイルManifest.toml
使用するパッケージのバージョンやインストール位置などを記録
  • ~/.julia/environments/v1.6/{Project, Manifest}.tomlがデフォルト(ユーザ固有)
  • パッケージ管理コマンドで自動的に更新

プロジェクトのパッケージ管理

  • プロジェクトのパッケージ管理
    1. プロジェクトのディレクトリに移動
    2. REPLを立ち上げてパッケージ管理モードに移行
    3. activate .コマンドを実行して環境を有効化
    4. パッケージを追加・削除など
  • プロジェクトのパッケージ環境有効化
    • julia --projectでJuliaを実行 OR
    • JULIA_PROJECT環境変数を@.に設定
  • プロジェクト固有の環境を作ってみる

Tips

LaTeX記法入力

  • ギリシャ文字などはLaTeX記法で入力可能
    • REPLで\sigmaの直後にTabキーを押すと、σに変換される
    • 意味が定義されている記号もある(例えばisapprox関数の別名)
  • VS Code, Vim, Emacs, Jupyterなどでも

julia> \sigma  # ← ここでTabキーを押下

julia> σ       # ← このようになる

julia> 1.0 ≈ 1.00000001  # isapprox
true
主なLaTeXコマンド
\div÷
\approx
\equiv
\in
\notin
\subseteq
\cap
\cup

型推論と最適化

  • 多重ディスパッチは型推論で決まる
    • 関数に渡された引数から式の型を推論
    • うまく書けばディスパッチのコストは無視できる
    • 型推論がうまくいくコードを書くことが性能を上げるためには重要
  • 例:ループ内のisless関数呼出し
    • isless関数は型に応じてディスパッチ
    • 周囲の文脈から引数の型は推論される

function partition!(xs, lo, hi)
    pivot = div(lo + hi, 2)
    pvalue = xs[pivot]
    xs[pivot], xs[hi] = xs[hi], xs[pivot]
    j = lo
    @inbounds for i in lo:hi-1
        if isless(xs[i], pvalue)
            xs[i], xs[j] = xs[j], xs[i]
            j += 1
        end
    end
    xs[j], xs[hi] = xs[hi], xs[j]
    return j
end

# 型推論結果を表示
code_typed(
    partition!,
    (Vector{Float32}, Int, Int),
    optimize = false)

マルチスレッド計算

  • マルチスレッドの並列化
    • julia -t4で4スレッドを起動
    • タスク(Task型)というコルーチンを並列実行できる
  • Threadsモジュール
    • @spawn: 並列実行タスクを起動
    • @sync: タスクを同期
    • @threads: forループを並列化

using .Threads: @threads, @spawn

# fooとbarを並列実行
@sync begin
    @spawn foo()
    @spawn bar()
end

# for文を並列実行
@threads for i in 1:100
    baz(i)
end

まとめ

  • Juliaはモダンな動的プログラミング言語
    • スクリプトを手軽に実行
    • REPLやJupyterで対話的実行
    • 充実した標準ライブラリ
  • 多重ディスパッチを中心パラダイムとした機能
    • 関数は型に応じた複数のメソッド(実装)を持つ
    • 四則演算から複雑な関数まで多重ディスパッチ
    • ユーザ定義型への拡張も多重ディスパッチ
    • 型が推論できれば多重ディスパッチのコストはゼロ
  • 優秀なパッケージ管理ツールを標準搭載
    • パッケージの追加・削除・更新もREPLから
    • プロジェクト固有の仮想環境も構築可能

参考資料