0. はじめに
本稿では Julia の for 文と while 文の基礎を説明します．
出力は記事を書く便宜上 print を多用しますが，実際上はすべて println で問題ありません．
Julia は MATLAB に多大なる影響を受けており，その典型が for です．1
以下反復の指数は j, k, l,... などを用います．整数を $j$ から始める主な理由は虚数単位 $i$ (または $\mathrm{i}$ ) と記号の衝突を避けること，タイピングで J がホームポジションで打てることです．
(虚数単位は Python では 1j, Julia では im です．)
1. 反復の for
1.1 for j in ...
他のほとんどすべてのプログラミング言語と同様に Julia でも繰り返しの文が書けます．
所謂 for loop です．Julia の for loop は MATLAB のものと同じです．2
公式には for j = が標準のようですがコーディングとして楽なので後の書き方と統一して for j in と書いています．
N = 10
for j in 1:N
print(j, " ")
end
println()
# 1 2 3 4 5 6 7 8 9 10
Python とは異なり最後に end を書かないとエラーになります．また Fortran の do...end do と似ています．
C/C++ と同じように改行については規則がなく自由です．ただし，特に理由がない限り通常改行とタブは推奨されます．これは スニペット機能を使えば特に気にする必要はありません．
1.2 for j = ...
for j = と書いても同じです．
N = 10
for j = 1:N# equal
print(j, " ")
end
println()
# 1 2 3 4 5 6 7 8 9 10
1.3 eachindex()
ベクトルや行列の添字に関して for 文を書きたい機会は頻繁に訪れます．
そんな時に便利なのが eachindex 関数です．
配列の添字を取得してくれます．
例えば1次元イジングモデルのハミルトニアンを考えてみましょう．
イジングモデルでは隣り合ったスピンの積の和を計算する必要があるため添え字に関して for loop を書くことが必要です．
4行目で eachindex 関数がスライスと組み合せて使われています．
ここでは 強磁性 $J = -1$ をデフォルトにし，パラメータ変数 args に格納します．
function ising_hamiltonian1D(A; args=(-1,))
J = args[1]
_ans = A[end] * A[1]
for j in eachindex(A)[1:end-1]
_ans += A[j] * A[j + 1]
# @show (j, _ans)
end
return J * _ans
end
ソースコード全体についてはこちらをご覧ください．
→ https://github.com/Krypf/20230323Qiita/blob/main/eachindex.jl
もう一つ考えうる，length 関数を使う方法は推奨されない場合があります．これについては以下をご参考にしてください．
1.4 逆順ループ
逆順は reverse() で計算します．
kaibun1 = "よのなかね，かおかおかねかなのよ．"
display(reverse(kaibun1))
for j in reverse(1:5)
print(j, " ")
end
print("Thunderbirds are go!\n")
"．よのなかねかおかおか，ねかなのよ"
5 4 3 2 1 Thunderbirds are go!
2. 集合の for
集合のすべての元についての操作，つまり数学の $∀ x ∈ X$ (for all $x$ in $X$) を書くには同様の for loop を使います．
プログラムでは集合の代わりにおおよそ順序付きの配列を用います．
代表的なものは Vector 型の配列です．
xs = ["a", "b", "c"]
for x in xs
print(x, " ")
end
println()
# a b c
もちろん集合のタイプでも出力可能ですがもし何らか順序を合わせたい時には使えません．
xs = Set(["a", "b", "c"])
for x in xs
print(x, " ")
end
println()
# c b a
3. 反復中の処理
3.1 break
ループを止めるには break を使います．MATLAB や Python と同じです．
例えば，j が4以上になったら for loop を止めてみましょう．
#%% stop if >= 4
# N = 10
for j in 1:N
if j >= 4
break
end
print(j, " ")
end
println()
# 1 2 3
このように，3までで止まっています．
次は冪が 閾値 を超える最小の自然数を求めるプログラムです．
素朴な総当りとして break の好例となっていると思います．
(2乗なら開平のアルゴリズムは筆算式で書けますから，その場合アルゴリズムとして冗長です)
function FindThresholdPower(threshold; M = 10, p = 2)
for j in 1:M
power_j_p = j^p
if power_j_p > threshold
println("j to the p is greater than the threshold when")
@show j, p, power_j_p, threshold
break
end
end
println("Done")
end
FindThresholdPower(50)
FindThresholdPower(5000, M = 100)
FindThresholdPower(Int(2e8); M = Int(1e3), p = 3)
FindThresholdPower(Int(1e15); M = 10, p = 17)
j to the p is greater than the threshold when
(j, p, power_j_p, threshold) = (8, 2, 64, 50)
Done
j to the p is greater than the threshold when
(j, p, power_j_p, threshold) = (71, 2, 5041, 5000)
Done
j to the p is greater than the threshold when
(j, p, power_j_p, threshold) = (585, 3, 200201625, 200000000)
Done
j to the p is greater than the threshold when
(j, p, power_j_p, threshold) = (8, 17, 2251799813685248, 1000000000000000)
Done
Python であれば速さのために array にして書いて下さい．
3.2 continue
ループで操作を飛ばすには continue を使います．
例えば，4 と 9 を飛ばすと次のようになります．
#%% omit 4 or 9 in 1:10
# N = 10
for j in 1:N
if j == 4 || j == 9
continue
end
print(j, " ")
end
println()
# 1 2 3 5 6 7 8 10
#%% は VSCode 対応のタブ機能です．
応用例として，外れ値をナイーブに除去し平均と標準偏差を更新するアルゴリズム ExcludeOutlier を考えてみましょう．
外れ値は偏差と標準偏差の比が大きいと考えられ，統計上含めないほうが良い可能性があります．
ですから，それが存在する時，配列からそれを取り除いて新たな平均値と標準偏差を求めることにはナイーブには意義があります．
function ExcludeOutlier(X, nsigma = 3)
_mean = mean(X)
_std = std(X)# standard deviation
println("original mean: ", _mean)
println("original stdev: ", _std)
n = length(X)::Int64
j = 0
_sum, _var = 0.0, 0.0
for x in X
a = abs(x - _mean) / _std# deviation
if a > nsigma# over threshold
j += 1
@show j,x,a
continue# skip
end
_sum += x
_var += (x)^2
end
_mean_new = _sum / (n - j)
_std_new = sqrt(( _var / (n - j) - _mean_new^2 ) * n / (n - 1))# convert the sample variance to the unbiased sample variance
return _mean_new, _std_new
end
# X = append!(collect(1:N),collect(1:N))
# X = append!(collect(1:N),10,collect(1:N))
# X = append!(collect(1:N),100,collect(1:N))
# X = append!(rand(N),10,(rand(N)))
X = append!(rand(N),10 * rand(),(rand(N)))
ExcludeOutlier(X)
今，$3 \sigma$ より偏差が大きければ統計上無視します．
コメントはテストコードです．（疑似）乱数を使っているので，出力は一例です．
original mean: 0.8943312769474329
original stdev: 1.6918161113329115
(j, x, a) = (1, 8.181866432588865, 4.307522021350114)
(0.5299545191653613, 0.27886564104290185)
平均値ではなく，中央値の偏差を用いてもよいでしょう．ただし，そうした場合標準偏差の意味との整合を考え直す必要はあります．
4. 内包表記の for
Python ではしばしばリスト内包表記が用いられますがJuliaでも同様に改行なしの内包表記 comprehensionができます．
N = Int(1e6)
Comprehension(N) = [j for j in 1:N]
@time x = Comprehension(N)
# 0.024938 seconds (46.11 k allocations: 10.067 MiB, 72.04% compilation time)
# 1000000-element Vector{Int64}:
# 1
# 2
# 3
# ⋮
# 999999
# 1000000
5. 多重ループ
複数の for loop は数値計算で頻出し重要です．数学的には直積集合の元から何らかの対応を取ることに相当します．
例えば 4次正方行列 $A \in \mathrm{M}(4,\mathbb{Z})$ があったとします：
A ≡ (A_{j,k})_{j,k} := \begin{pmatrix}
1 & 2 & 3 & 4 \\
5 & 6 & 7 & 8 \\
9 & 10 & 11 & 12 \\
13 & 14 & 15 & 16
\end{pmatrix}.
添字集合を
J = \{1,2,3,4\}
とおくのが便利です．
そして $\mathscr{A}\colon (j,k)\in J \times J \to \mathbb{Z} $,
\mathscr{A} (j,k) = 4(j - 1) + k
とするとこの場合 $\mathscr{A} (j,k) = A_{j,k}$ です．例えば $4(2-1)+3 = 7 = A_{2,3}$.
複数の添字 $j,k$ が走って $(j,k) = (2,3) \iff j = 2 \land k = 3$ のときの $\mathscr{A}$ の値が $A_{2,3}$ というわけです．
このように行列を始めとする多次元配列の計算と多重ループは，本質を一にしています．
5.1 入れ子の2重ループ
行列を2重ループで生成できます．
上の例をコードにしてみましょう．
N = 4
A = [[0 for _ in 1:N] for _ in 1:N]
for j in 1:N
for k in 1:N
A[j][k] = 4 * (j - 1) + k
end
end
display(A)
display(A[1])
4-element Vector{Vector{Int64}}:
[1, 2, 3, 4]
[5, 6, 7, 8]
[9, 10, 11, 12]
[13, 14, 15, 16]
4-element Vector{Int64}:
1
2
3
4
Julia では Matrix 型が用意されているのでここで使っている Vector の Vector は型としてはあくまで2次的なものです．
5.2 1行立ての2重ループ
数学では $∀ j, k \in J$ という略記がありますが Julia では $∀ j \in J, k\in J$ という書き方ができます．
必ず覚えましょう．
N = 4
A = [[0 for _ in 1:N] for _ in 1:N]
for j in 1:N, k in 1:N
A[j][k] = 4 * (j - 1) + k
end
display(A)
@show A[2][3]
4-element Vector{Vector{Int64}}:
[1, 2, 3, 4]
[5, 6, 7, 8]
[9, 10, 11, 12]
[13, 14, 15, 16]
(A[2])[3] = 7
7
数学の略記法の方は使えません．
# for j, k in 1:N
# A[j][k] = 4 * (j - 1) + k
# end
# ERROR: syntax: invalid iteration specification
5.3 3重ループ
2重ループと同様に繰り返して下さい．
function ShowAllIndices(N)
for j in 1:N, k in 1:N, l in 1:N
print((j,k,l), " ")
if k == N && l == N
println()
end
end
end
N = 3
ShowAllIndices(N)
# (1, 1, 1) (1, 1, 2) (1, 1, 3) (1, 2, 1) (1, 2, 2) (1, 2, 3) (1, 3, 1) (1, 3, 2) (1, 3, 3)
# (2, 1, 1) (2, 1, 2) (2, 1, 3) (2, 2, 1) (2, 2, 2) (2, 2, 3) (2, 3, 1) (2, 3, 2) (2, 3, 3)
# (3, 1, 1) (3, 1, 2) (3, 1, 3) (3, 2, 1) (3, 2, 2) (3, 2, 3) (3, 3, 1) (3, 3, 2) (3, 3, 3)
6. 条件付きの for
内包表記の for loop に条件をつけることができます．
Julia の場合 そのまま if 文を続ければ良いです．
$1\le j \le k \le l\le n$ の自然数解を求めるプログラムです．
総数は
(N+2)(N+1)N / 6 = \begin{pmatrix}
n \\
3
\end{pmatrix}
と求まります．
function AllCombination(N)
[(j,k,l) for j in 1:N, k in 1:N, l in 1:N if j <= k && k <= l]
end
N = 4
display(AllCombination(N))
# show(stdout, "text/plain", AllCombination(N))
# (N + 2) * (N + 1) * N // 6
20-element Vector{Tuple{Int64, Int64, Int64}}:
(1, 1, 1)
(1, 1, 2)
(1, 2, 2)
⋮
(3, 4, 4)
(4, 4, 4)
通常の for 文で書いてもほとんど同じです．
Julia では改行は自由で最後を end で閉じるのでした．
for j in 1:N, k in 1:N, l in 1:N if j <= k && k <= l
@show (j,k,l)
end
end
7. while
最後に while 文です．
7.1 while と for
while loop によって for loop を書くことができます．
function sum1toN(N)
_ans = 0
j = 1# count
while j <= N
_ans += j
j += 1
end
return _ans
end
function sum1toNfor(N)
_ans = 0
for j in 1:N
_ans += j
end
return _ans
end
N = 4
@show sum1toN(N), sum1toNfor(N)
# (sum1toN(N), sum1toNfor(N)) = (10, 10)
# (10, 10)
7.2 無限ループ
while loop では絶対に無限ループを避けるための注意を払いましょう．
無限ループが生じるとどういうことになるのか具体例で解説しましょう．
while 文で，条件が恒真であったり，再帰操作が終了するトリガーがなかったりすると，プログラムは終了しません．
このようなバグは絶対に避けるべきです．
最悪の場合メモリオーバーフローを起こしたり，人間がプログラムを止めようにも止められなくなったりするからです．
そこで今回は threshold 秒経過することを終了のトリガーにしてリミッター（安全装置）を設定した上で，無限ループを書いてみましょう．
using Dates
function endless_eight(threshold)
t0 = now()# Dates library. today's date with a unit millisecond.
j = 0# counts
while 0 < 1# infinite loop!
j += 1
t1 = now()
Δt = t1 - t0# elapsed time
# print(Δt, " ")# Julia takes much time to print = low speed
if Δt >= Dates.Second(threshold)
print("\n")
println(t0)
println(t1)
break
end
end
return "今回が$(j)回目に該当する"
end
threshold = 8
println(endless_eight(threshold))
2023-03-22T10:58:06.716
2023-03-22T10:58:14.716
今回が35454812回目に該当する
コンピュータは無謬の存在ではありません．桁落ち誤差や計算の有限性を避けては通れません．
しかし反復操作は人間より得意です．人間を超えた冷徹なまでの繰り返しが while loop にはその本質として立ち現れているのです．
7.3 Collatz 予想*
while loop の好例は Collatz 予想です．
漸化式
a_{n+1}=\begin{cases}a_{n}/2 & (n ≡ 0 \mod 2),\\
3a_{n}+1 & (n ≡ 1).\end{cases}
によって，いかなる正の整数を初項としても $∃ M \in \mathbb{N}, a_M = 1$ という主張が Collatz 予想です．これは未解決です．
アルゴリズムは単純です．プログラムは漸化式と while 文で書くことができます．
# The Collatz conjecture
# one step
function recurrence(n::Int64)# recurrence relation
if n % 2 == 0# n is even.
n ÷= 2
else# n is odd
n = 3 * n + 1
end
return n
end
# whole process
function collatz(n::Int64)
counts = 0# the number of calculations (steps)
while n > 1
n = recurrence(n)
counts += 1
print((counts, n), " ")
end
println()
return (counts, n)
end
N = Int(12345)
collatz(N)
(1, 37036) (2, 18518) (3, 9259) (4, 27778) (5, 13889) (6, 41668) (7, 20834) (8, 10417) (9, 31252) (10, 15626) (11, 7813) (12, 23440) (13, 11720) (14, 5860) (15, 2930) (16, 1465) (17, 4396) (18, 2198) (19, 1099) (20, 3298) (21, 1649) (22, 4948) (23, 2474) (24, 1237) (25, 3712) (26, 1856) (27, 928) (28, 464) (29, 232) (30, 116) (31, 58) (32, 29) (33, 88) (34, 44) (35, 22) (36, 11) (37, 34) (38, 17) (39, 52) (40, 26) (41, 13) (42, 40) (43, 20) (44, 10) (45, 5) (46, 16) (47, 8) (48, 4) (49, 2) (50, 1)
(50, 1)
→ Collatz conjecture - Wikipedia
8. おわりに
アルゴリズムの基本である for 文と while 文について説明しました．
Julia では Python よりも多くの場面で for 文を使うことができますので必ず覚えることが推奨されます．
-
Why We Created Julia: "We are power Matlab users." ↩