【Julia】基礎: for 文 と while 文

<!-- 【Julia】 基礎 2. for 文 と while 文 -->

# 0. はじめに
本稿では Julia の for 文と while 文の基礎を説明します．

出力は記事を書く便宜上 `print` を多用しますが，実際上はすべて `println` で問題ありません．

Julia は MATLAB に多大なる影響を受けており，その典型が for です．[^1]

以下反復の指数は `j, k, l,...` などを用います．整数を $j$ から始める主な理由は虚数単位 $i$ (または $\mathrm{i}$ ) と記号の衝突を避けること，タイピングで `J` がホームポジションで打てることです．
(虚数単位は Python では `1j`, Julia では `im` です．)

[^1]: [Why We Created Julia](https://julialang.org/blog/2012/02/why-we-created-julia/): "We are power Matlab users."

# 1. 反復の for

## 1.1 for j in ...
他のほとんどすべてのプログラミング言語と同様に Julia でも繰り返しの文が書けます．

所謂 **for loop** です．Julia の for loop は MATLAB のものと同じです．[^2]
[^2]: [for ループを指定した回数で繰り返す - MATLAB for - MathWorks 日本](https://jp.mathworks.com/help/matlab/ref/for.html)

公式には `for j =` が標準のようですがコーディングとして楽なので後の書き方と統一して `for j in ` と書いています．

```julia
N = 10
for j in 1:N
    print(j, " ")
end
println()
# 1 2 3 4 5 6 7 8 9 10 
```

Python とは異なり<ins>最後に `end` を書かないとエラー</ins>になります．また Fortran の `do...end do` と似ています．
C/C++ と同じように改行については規則がなく自由です．ただし，特に理由がない限り通常改行とタブは推奨されます．これは スニペット機能を使えば特に気にする必要はありません．

## 1.2 for j = ...

`for j =` と書いても同じです．

```julia
N = 10
for j = 1:N# equal
    print(j, " ")
end
println()
# 1 2 3 4 5 6 7 8 9 10 
```

## 1.3 `eachindex()`


<!-- <div style="padding: 10px; margin-bottom: 10px; border: 3px solid #333333;"> -->
:::note info
ベクトルや行列の添字に関して for 文を書きたい機会は頻繁に訪れます．
そんな時に便利なのが `eachindex` 関数です．
配列の添字を取得してくれます．
:::
<!-- </div> -->


例えば1次元イジングモデルのハミルトニアンを考えてみましょう．
イジングモデルでは隣り合ったスピンの積の和を計算する必要があるため添え字に関して for loop を書くことが必要です．
4行目で `eachindex` 関数がスライスと組み合せて使われています．
ここでは 強磁性 $J = -1$ をデフォルトにし，パラメータ変数 `args` に格納します．


```julia
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](https://github.com/Krypf/20230323Qiita/blob/main/eachindex.jl)

もう一つ考えうる，`length` 関数を使う方法は推奨されない場合があります．これについては以下をご参考にしてください．

[1:length(A) considered harmful — or how to make Julia code “generic safe” | by Roland Schätzle | Towards Data Science](https://towardsdatascience.com/1-length-a-considered-harmful-or-how-to-make-julia-code-generic-safe-ac7b39cfc2f0)

## 1.4 逆順ループ

逆順は `reverse()` で計算します．

```julia
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` 型の配列です．

```julia
xs = ["a", "b", "c"]
for x in xs
    print(x, " ")
end
println()
# a b c 
```

もちろん集合のタイプでも出力可能ですがもし何らか順序を合わせたい時には使えません．

```julia
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 を止めてみましょう．

```julia
#%% 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乗なら開平のアルゴリズムは筆算式で書けますから，その場合アルゴリズムとして冗長です)

```julia
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 を飛ばすと次のようになります．

```julia
#%% 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 
```

<font size = 2> `#%%` は VSCode 対応のタブ機能です． </font>

---

応用例として，外れ値をナイーブに除去し平均と標準偏差を更新するアルゴリズム `ExcludeOutlier` を考えてみましょう．
外れ値は偏差と標準偏差の比が大きいと考えられ，統計上含めないほうが良い可能性があります．
ですから，それが存在する時，**配列からそれを取り除いて新たな平均値と標準偏差を求める**ことにはナイーブには意義があります．

```julia
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)
```

<font size = 2> 平均値ではなく，中央値の偏差を用いてもよいでしょう．ただし，そうした場合標準偏差の意味との整合を考え直す必要はあります． </font>

# 4. 内包表記の for

Python ではしばしばリスト内包表記が用いられますがJuliaでも同様に改行なしの**内包表記 comprehension**ができます．

```julia
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})$ があったとします：

```math
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}.
```

添字集合を

```math
J = \{1,2,3,4\}
```

とおくのが便利です．
そして $\mathscr{A}\colon (j,k)\in J \times J \to \mathbb{Z} $,  
```math
\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}$ というわけです．


<!-- <div style="padding: 10px; margin-bottom: 10px; border: 3px solid #333333;"> -->
:::note info
このように行列を始めとする多次元配列の計算と多重ループは，本質を一にしています．
:::
<!-- </div>  -->


## 5.1 入れ子の2重ループ 

行列を2重ループで生成できます．
上の例をコードにしてみましょう．

```julia
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
```

<font size = 2> Julia では Matrix 型が用意されているのでここで使っている Vector の Vector は型としてはあくまで2次的なものです．</font>


## 5.2 1行立ての2重ループ

数学では $∀ j, k \in J$ という略記がありますが Julia では $∀ j \in J, k\in J$ という書き方ができます．
必ず覚えましょう．

```julia
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
```

数学の略記法の方は使えません．

```julia
# for j, k in 1:N
#     A[j][k] = 4 * (j - 1) + k
# end
# ERROR: syntax: invalid iteration specification
```

## 5.3 3重ループ

2重ループと同様に繰り返して下さい．

```julia
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$ の自然数解を求めるプログラムです．
総数は

```math
(N+2)(N+1)N / 6 = \begin{pmatrix}
n \\
3
\end{pmatrix}
```

と求まります．

```julia
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` で閉じるのでした．

```julia
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 を書くことができます．

```julia
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` 秒経過することを終了のトリガーにしてリミッター（安全装置）を設定した上で，無限ループを書いてみましょう．

```julia
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 予想です．
漸化式

```math
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 文で書くことができます．

```julia

# 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](https://en.wikipedia.org/wiki/Collatz_conjecture)

# 8. おわりに

アルゴリズムの基本である `for` 文と `while` 文について説明しました．

Julia では Python よりも多くの場面で `for` 文を使うことができますので必ず覚えることが推奨されます．