# 실전 줄리아 첫걸음 연습문제

## 알고리즘

### 015. n번째로 큰 인덱스 ☆

`argmax`는 최대값의 위치를, `argmin`는 최소값의 위치를 리턴하는 함수다. 이와 유사하게, 부분순서배열 `A`에서 `n`번째로 큰 원소의 인덱스를 리턴하는 `argnth(n, A)`를 작성하라.

#### 검증예시

julia> X = [0, -1, 7, 3, 9, 2, 11];

julia> argnth(2, X)
5

julia> argnth(3, X)
3

#### 답안예시

argnth(n, A) = sortperm(A, rev = true)[n]

<!-- ### 001. 기수 정렬 ☆☆

#### 답안예시

function row_sort(M, i)
    m = M[i,:]
    idx = vcat([findall(m .== "$digit") for digit in 0:9]...)
    return M[:, idx]
end
function radix_sort(arry)
    k = Int64(maximum(ceil.(log10.(arry))))
    M = stack(split.(lpad.(arry, k, '0'), ""))

    for i in k:(-1):1
        M = row_sort(M, i)
    end
    return parse.(Int64, reduce.(*, eachcol(M)))
end

arry = [190, 812, 9, 77, 21, 413]
radix_sort(arry)

### 002. 가중치 샘플링 ☆

$\mathbb{R}_{-}^{c} = \left\{ r \in \mathbb{R} : r \ge 0 \right\}$ 이라고 하자. 임의의 $d$-튜플 $\mathbf{I}$ 와 영벡터가 아닌 벡터 $\mathbf{w} \in \left( \mathbb{R}_{-}^{c} \right)^{d}$ 를 입력받아서 $\frac{\left( \mathbf{w} \right)_{k}}{\sum_{k} \left( \mathbf{w} \right)_{k}}$ 의 확률로 $\mathbf{I}$ 의 $k$번째 원소 $\left( \mathbf{I} \right)_{k}$ 을 리턴하는 함수 `randw(item, weight)`를 작성하라.

#### 검증예시

julia> randw(['a', 'b', 'c'], [1,2])
ERROR: AssertionError: length mismatch

julia> randw(['a', 'b', 'c'], [0,0,0])
ERROR: AssertionError: all weights are zero

julia> randw(['a', 'b', 'c'], [1,2,-1])
ERROR: AssertionError: there is a negative weight

julia> randw(['a', 'b', 'c'], [1,2,3])
'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)

julia> abc = [randw(['a', 'b', 'c'], [1,2,3]) for _ in 1:6000];

julia> count.([abc .== 'a', abc .== 'b', abc .== 'c'])
3-element Vector{Int64}:
 1029
 1980
 2991

<!-- #### 답안예시

function randw(item, weight)
    @assert length(item) == length(weight) "length mismatch"
    @assert any(weight .> 0) "all weights are zero"
    @assert all(weight .≥ 0) "there is a negative weight"
    
    Σ = sum(weight)
    if !isone(Σ)
        weight = weight ./ Σ
    end
    return item[findfirst(rand() .≤ cumsum(weight))]
end
randw(['a', 'b', 'c'], [1,2])
randw(['a', 'b', 'c'], [0,0,0])
randw(['a', 'b', 'c'], [1,2,-1])
randw(['a', 'b', 'c'], [1,2,3])
abc = [randw(['a', 'b', 'c'], [1,2,3]) for _ in 1:6000]
count.([abc .== 'a', abc .== 'b', abc .== 'c'])

### 003. 비복원 샘플링 ☆☆

**선수문제:** 002. 가중치 샘플링

`randw` 함수를 사용해서 가중치 샘플링을 `n`번 반복하되, 비복원 샘플링(sample without replacement)을 할 수 있는 함수 `randw(item, weight, n; replace = true)`를 재정의하라. `replace` 키워드의 디폴트는 `true`로 주어지며, `false`일 때 비복원 샘플링이 수행된다.

#### 검증예시

julia> randw(1:10, 1:10, 5)
5-element Vector{Any}:
 9
 7
 6
 6
 6

julia> randw(1:10, 1:10, 5, replace = false)
5-element Vector{Int64}:
  6
  4
  3
 10
  8

#### 답안예시

다음은 샘플 벡터가 제약조건을 만족시키지 못하면 폐기하고 원하는 결과가 나올때까지 반복하는 방법이다.

function randw(item, weight, n; replace = true)
    @assert n ≤ length(item) "n = $n is larger than the length of item = $(length(item))"

    while true
        sample = [randw(item, weight) for _ in 1:n]
        if allunique(sample) || replace
            return sample
        end
    end
end

다음은 샘플된 아이템을 하나씩 금지 시켜가며 가중치를 업데이트하는 방법이다. 이 방법의 결과로 얻는 샘플 벡터의 원소들이 적법한 순서를 가지는지 논하라.

function randw(item, weight, n; replace = true)
    @assert n ≤ length(item) "n = $n is larger than the length of item = $(length(item))"
    N = length(item)

    if replace
        sample = [randw(item, weight) for _ in 1:n]
    else
        weight = deepcopy(collect(weight))
        sample = []
        for _ in 1:n
            k = randw(1:N, weight)
            weight[k] = 0
            push!(sample, item[k])
        end
    end
    return sample # OK with this order?
end

## 네트워크

### 012. 길버트 모델 ☆

[길버트 모델](https://freshrimpsushi.github.io/posts/2118)을 통해 [에르되시-레니 그래프](https://freshrimpsushi.github.io/posts/2116)와 유사한 결과를 얻을 수 있다. 노드의 수 $n \in \mathbb{N}$ 과 확률 $p \in [0,1]$ 을 입력받아 길버트 모델에 따라 랜덤 그래프를 리턴하는 함수 `gilbert(n, p)`를 작성하라.

### 013. 유클리드 그래프 ☆☆

노드의 수 $n \in \mathbb{N}$ 과 컷오프 $\delta \ge 0$ 를 입력받아 [유클리드 그래프](https://freshrimpsushi.github.io/posts/2396)을 리턴하는 함수 `euclidean(n, δ)`를 작성하라.

### 004. 바라바시-알버트 모델 ☆☆☆

**선수문제:** 003. 비복원 샘플링

최종 노드의 수 $N \in \mathbb{N}$ 과 링크 파라미터 $m \in \mathbb{N}$ 를 입력받아 [바라바시-알버트 모델](https://freshrimpsushi.github.io/posts/2187)을 통해 [스케일 프리 네트워크](https://freshrimpsushi.github.io/posts/2183)를 리턴하는 함수 `barabasi_albert(N, m = 1)`를 작성하라.

#### 검증예시

julia> ba = barabasi_albert(1000, 3)
1000-element Vector{Vector{Int64}}:
 [2, 3, 4, 5, 6, 7, 9, 11, 15, 18  …  826, 860, 864, 872, 908, 913, 929, 942, 960, 961]
 [1, 3, 4, 5, 9, 17, 18, 38, 46, 55  …  704, 710, 720, 829, 837, 846, 922, 962, 971, 999]
 [1, 2, 4, 5, 6, 7, 8, 12, 14, 22  …  389, 402, 503, 551, 560, 565, 617, 735, 833, 897]
 [2, 3, 1, 8, 10, 11, 14, 19, 20, 22  …  726, 779, 797, 820, 866, 894, 900, 926, 937, 984]
 [1, 3, 2, 6, 7, 8, 9, 10, 13, 15  …  736, 762, 773, 848, 857, 922, 923, 926, 928, 979]
 ⋮
 [53, 739, 342]
 [83, 142, 87]
 [447, 9, 19]
 [495, 2, 236]
 [530, 190, 359]

julia> degree = length.(ba)
1000-element Vector{Int64}:
 96
 42
 36
 63
 72
  ⋮
  3
  3
  3
  3
  3
  
julia> histogram(degree, xscale = :log10, yscale = :log10, xticks = 10.0 .^(-0:0.5:2), xlims = extrema(degree))


#### 답안예시

function randw(item, weight)    
    Σ = sum(weight)
    if !isone(Σ)
        weight = weight ./ Σ
    end
    return item[findfirst(rand() .≤ cumsum(weight))]
end
function randw(item, weight, n; replace = true)
    while true
        sample = [randw(item, weight) for _ in 1:n]
        if (length(unique(sample)) == n) || replace
            return sample
        end
    end
end

function barabasi_albert(N, m = 1)
    G = [setdiff(1:m, k) for k in 1:m]
    for new in (m+1):N
        idxs = randw(1:(new-1), length.(G), m, replace = false)
        push!(G, idxs)
        push!.(G[idxs], new)
    end
    return G
end

ba = barabasi_albert(1000, 3)
degree = length.(ba)

using Plots
histogram(degree)
histogram(degree, xscale = :log10, yscale = :log10, xticks = 10.0 .^(-0:0.5:2), xlims = extrema(degree))

## 선형대수

### 007. 라플라스 전개 ☆

정사각행렬 $A = \left( a_{ij} \right) \in \mathbb{R}^{n \times n}$ 에 대해 $i$번째 행과 $j$번째 열을 제거한 행렬의 행렬식^determinant $M_{ij}$ 를 소행렬식^minor, $C_{ij} := (-1)^{i+j}M_{ij}$ 를 코팩터^cofactor라 한다.
$$
\det A = \sum_{j=1}^{n} a_{ij} C_{ij} \qquad \text{, for fixed row } i
$$

위의 공식을 라플라스 전개^{Laplace expansion}라 한다. 정사각행렬을 입력받아 라플라스 전개를 리턴하는 함수 `laplace(A)`를 작성하라.

#### 검증예시

julia> B = [2 3; 4 6]
2×2 Matrix{Int64}:
 2  3
 4  6

julia> laplace(B)
0

julia> C = [7 1; 9 3]
2×2 Matrix{Int64}:
 7  1
 9  3

julia> laplace(C)
12

julia> M = rand(3,3)
3×3 Matrix{Float64}:
 0.478437  0.779194  0.334746
 0.342975  0.867767  0.0211038
 0.733843  0.97752   0.331435

julia> laplace(M)
-0.04971335247679029

julia> using LinearAlgebra

julia> det(M)
-0.049713352476790304

julia> abs(det(M) - laplace(M))
1.3877787807814457e-17

<!-- #### 답안예시

function laplace(A::AbstractMatrix)
    n = size(A, 1)
    @assert n == size(A, 2) "matrix $A is not square"

    if n |> isone
        aijCij = A
    else
        aijCij = [A[1,j] * (-(-1)^j) * laplace(A[2:end, setdiff(1:n, j)]) for j in 1:n]
    end
    return sum(aijCij)
end
B = [2 3; 4 6]
laplace(B)
C = [7 1; 9 3]
laplace(C)

M = rand(3,3)
laplace(M)

using LinearAlgebra
det(M)
abs(det(M) - laplace(M))

### 008. 한켈 매트릭스 ☆

$$
H = \begin{bmatrix}
h_{11} & h_{12} & h_{13} & \cdots & h_{1n}
\\\ h_{21} & h_{22} & h_{23} & \cdots & h_{2n}
\\\ h_{31} & h_{32} & h_{33} & \cdots & h_{3n}
\\\ \vdots & \vdots & \vdots & \ddots & \vdots
\\\ h_{m1} & h_{m2} & h_{m3} & \cdots & h_{mn}
\end{bmatrix}
$$

주어진 행렬 $H = \left( h_{ij} \right) \in \mathbb{R}^{m \times n}$ 가 모든 $k = 2 , 3 , \cdots , m+n$ 에 대해 다음을 만족하면 **한켈 행렬**^{Hankel Matrix}라 한다.
$$
i_{1} + j_{1} = k = i_{2} + j_{2} \implies h_{i_{1} j_{1}} = h_{i_{2} j_{2}}
$$

다시 말해, 한켈 행렬이란 다음과 같이 반대각선의 성분이 모두 같은 행렬이다.
$$
H = \begin{bmatrix}
h_{2} & h_{3} & h_{4} & \cdots & h_{1+n}
\\\ h_{3} & h_{4} & h_{5} & \cdots & h_{2+n}
\\\ h_{4} & h_{5} & h_{6} & \cdots & h_{3+n}
\\\ \vdots & \vdots & \vdots & \ddots & \vdots
\\\ h_{m+1} & h_{m+2} & h_{m+3} & \cdots & h_{m+n}
\end{bmatrix}
$$

벡터 $\left( h_{2}, \cdots , h_{m+n} \right)$ 과 두 자연수 $m, n$ 을 입력받아 위와 같은 한켈 행렬 $H \in \mathbb{R}^{m \times n}$ 을 리턴하는 함수 `hankelizer(vec, m, n)`을 작성하라.

#### 검증예시

julia> elements = [8,2,0,6,5,1,5,4,0,1];

julia> hankelizer(elements, 3, 8)
3×8 Matrix{Int64}:
 8  2  0  6  5  1  5  4
 2  0  6  5  1  5  4  0
 0  6  5  1  5  4  0  1

julia> hankelizer(elements, 2, 9)
2×9 Matrix{Int64}:
 8  2  0  6  5  1  5  4  0
 2  0  6  5  1  5  4  0  1

<!--
#### 답변예시

function hankelizer(vec)
    @assert length(vec) % 2 == 1 "vector must have odd length"

    m = (length(vec) ÷ 2) + 1
    H = zeros(eltype(vec), m, m)
    for i ∈ 1:m, j ∈ 1:m
        H[i,j] = vec[i+j-1]
    end
    return H
end

elements = [8,2,0,6,5,1,5,4,0]
hankelizer(elements)


function hankelizer(vec, m, n)
    @assert length(vec) == (m+n-1) "wrong dimensions"

    H = zeros(eltype(vec), m, n)
    for i ∈ 1:m, j ∈ 1:n
        H[i,j] = vec[i+j-1]
    end
    return H
end

elements = [8,2,0,6,5,1,5,4,0,1];
hankelizer(elements, 3, 8)
hankelizer(elements, 2, 9)

-->
### 011. 역행렬 ☆☆

주어진 행렬 $M \in \mathbb{R}^{n \times n}$ 에 대해 가우스 소거법을 통해 역행렬 $M^{-1}$ 을 리턴하는 함수 `gauss(M)` 을 작성하라.

### 014. 특이값 분해 ☆☆

행렬 $A \in \mathbb{R}^{n \times n}$ 에 대해 $A^{T} A$ 가 $\operatorname{rank} A = r$ 이라 하자. $A$ 의 고유값이 $\lambda_{1} \ge \cdots \lambda_{n} \ge 0$ 과 같이 정렬 되어 있다고 하고, 각자에 대응되는 고유벡터를 $\mathbf{v}_{1} , \cdots , \mathbf{v}_{n}$ 이라고 할 때, 다음과 같은 계산을 할 수 있다.
$$
\begin{align}
\sigma_{j} :=& \sqrt{\lambda_{j}} & , j = 1 , \cdots , n \\
\mathbf{u}_{j} :=& \frac{1}{\sigma_{j}} A \mathbf{v}_{j} & , j \in 1 , \cdots , r
\end{align}
$$

두 직사각행렬 $\hat{U} := \left( \mathbf{u}_{1}^{T} , \cdots , \mathbf{u}_{r}^{T} \right) \in \mathbb{R}^{m \times r}$ 과 $\hat{V} := \left( \mathbf{v}_{1}^{T} , \cdots , \mathbf{v}_{r}^{T} \right) \in \mathbb{R}^{n \times r}$ 과 대각행렬 $\hat{\Sigma} := \operatorname{diag} \left( \sigma_{1} , \cdots , \sigma_{r} \right) \in \mathbb{R}^{r \times r}$ 에 대해
$$
A = \hat{U} \hat{\Sigma} \hat{V}^{T}
$$

가 성립하고, 이를 축소 특이값 분해^{reduced singular value decomposition}이라 한다. $A$ 의 모든 고유값과 고유벡터를 알고 있다고 가정할 수 있을 때, $A$ 의 특이값 분해를 리턴하는 함수 `mySVD(A)`를 작성하라.
- 문제의 풀이와 검증을 위해서는 `LinearAlgebra.jl` 라이브러리의 `rank`, `eigen`, `diagm`, `svd` 와 같은 함수를 사용해도 좋다.

#### 검증예시

julia> A = [1 2 3
            4 5 6
            7 8 7
            4 2 1];

julia> myU, myΣ, myVt = mySVD(A);

julia> myU*myΣ*myVt'
4×3 Matrix{Float64}:
 1.0  2.0  3.0
 4.0  5.0  6.0
 7.0  8.0  7.0
 4.0  2.0  1.0

julia> B = [1 2 3 4 5
            2 4 6 8 10
            7 8 7 6 5]
3×5 Matrix{Int64}:
 1  2  3  4   5
 2  4  6  8  10
 7  8  7  6   5

julia> myU, myΣ, myVt = mySVD(B);

julia> myU*myΣ*myVt'
3×5 Matrix{Float64}:
 1.0  2.0  3.0  4.0   5.0
 2.0  4.0  6.0  8.0  10.0
 7.0  8.0  7.0  6.0   5.0

SVD가 원래의 폼으로 되돌아 오면 분해 자체는 성공한 것이다. 표준 라이브러리인 `LinearAlgebra.jl`과 비교하면 다음과 같다. 차원에는 차이가 있으나 근본적인 수치들, 특히 고유값의 경우 정확히 같음을 확인하라.

julia> myU
3×2 Matrix{Float64}:
  0.295963  -0.33527
  0.591925  -0.670541
 -0.749687  -0.661792

julia> U
3×3 Matrix{Float64}:
 -0.33527   -0.295963  -0.894427
 -0.670541  -0.591925   0.447214
 -0.661792   0.749687   1.66533e-16

julia> myΣ
2×2 Matrix{Float64}:
 6.27906   0.0
 0.0      21.4143

julia> Σ
3-element Vector{Float64}:
 21.414326423487807
  6.2790623367117195
  0.0

julia> myVt
5×2 Matrix{Float64}:
 -0.600089  -0.294611
 -0.48381   -0.403797
 -0.128741  -0.451175
  0.226328  -0.498552
  0.581397  -0.54593

julia> Vt
5×3 adjoint(::Matrix{Float64}) with eltype Float64:
 -0.294611   0.600089  -0.594015
 -0.403797   0.48381    0.769135
 -0.451175   0.128741  -0.213693
 -0.498552  -0.226328  -0.0979743
 -0.54593   -0.581397   0.0177441

#### 답안예시

using LinearAlgebra

function mySVD(A::AbstractMatrix)
    r = rank(A)
    λ, V = eigen(A'A) # we can assume eigenpairs are given
    λ = λ[  (end-r+1):end]
    V = V[:,(end-r+1):end]

    Σ = diagm(sqrt.(λ))
    U = (A * V) / Σ
    return U, Σ, V
end

A = [1 2 3
     4 5 6
     7 8 7
     4 2 1];
myU, myΣ, myVt = mySVD(A);
myU*myΣ*myVt'
B = [1 2 3 4 5
     2 4 6 8 10
     7 8 7 6 5]
myU, myΣ, myVt = mySVD(B);
myU*myΣ*myVt'

U, Σ, Vt = svd(B);
U*diagm(Σ)*Vt'

myU
U
myΣ
Σ
myVt
Vt

## 해석학

### 009. 극값 ☆

주어진 시퀀스 $x_{t}$ 에 대해 $x_{t} \ge x_{t-1}, x_{t+1}$ 혹은 $x_{t} \le x_{t-1}, x_{t+1}$ 를 만족하면 $x_{t}$ 를 **극값**^{extreme value}이라 정의한다. 배열 `arr`$=x_{t}$ 를 입력받아 모든 극값의 위치를 리턴하는 함수 `argext(arr)` 를 작성하라.

<!--
#### 답안 예시

function argext(arr)
    idxmin = (circshift(arr, 1) .> arr) .&& (circshift(arr, -1) .> arr)
    idxmax = (circshift(arr, 1) .< arr) .&& (circshift(arr, -1) .< arr)
    return idxmin .|| idxmax
end

-->

### 010. 바이섹션 메소드 ☆

주어진 함수 $f : \mathbb{R} \to \mathbb{R}$ 가 $f(a) f(b) < 0$ 를 만족한다고 하자. [바이섹션 메소드](https://freshrimpsushi.github.io/posts/676)를 구현한 함수 `bisect(f,a,b;atol=eps())` 를 작성하라. `atol`은 절대허용오차(absolute tolerance)를 의미한다.

#### 검증예시

julia> bisect(x -> x^3 + 8, -7, 7)
-1.9999999850988388

julia> bisect(x -> sin(x), 3, 4)
3.141592025756836

<!-- #### 답안예시

function bisect(f,a,b; atol = 10^(-6))
    @assert f(a)*f(b) < 0 "f(a) and f(b) have same sign"

    c = (a+b)/2
    while abs(f(c)) > atol
        if f(b)*f(c) < 0
            a = c
        else
            b = c
        end
        c = (a+b)/2
    end
    return c
end
bisect(x -> x^3 + 8, -7, 7)
bisect(x -> sin(x), 3, 4)

## 다이내믹스

### 005. 로지스틱 맵의 트래젝터리 ☆

$$
x_{t+1} = f_{r} \left( x_{t} \right) = r x_{t} \left( 1 - x_{t} \right)
$$

위와 같이 로지스틱 패밀리^{logistic family}를 통해 정의된 다이내믹 시스템이 주어져 있다고 하자.
1. 파라미터 $r$ 과 $x_{t}$ 에 대해 $x_{t+1}$ 을 리턴하는 함수 `logistic(r, x)`을 작성하라.
2. 맵 $f_{r}$ 와 파라미터 $r$, 초기값 $x_{0}$, 반복횟수 $n$을 입력 받아 오빗^orbit $\left\{ x_{t} \right\}_{t=0}^{n}$ 을 리턴하는 함수 `iter(f, p, x0, n)`을 작성하라.

#### 검증예시

julia> logistic(4, 0.5)
1.0

julia> traj = iter(logistic, 3, 0.5, 100);

julia> using Plots

julia> timeevolution = plot(traj, color = :black, legend = :none)


<!-- #### 답안예시

function logistic(r, x)
    return r * x * (1 - x)
end
logistic(4, 0.5)

function iter(f, p, x0, n)
    x = [x0]
    for t in 1:n
        push!(x, f(p, x[end]))
    end
    return x
end
traj = iter(logistic, 3, 0.5, 100)

using Plots
plot(traj)

### 006. 로지스틱 맵의 바이퍼케이션 ☆☆

<!--
#### 답안예시

function bifurcation(h)
xaxis_ = []
yaxis_ = []
for r ∈ 1:h:4
preiterated = iter(logistic, r, 0.5, 200)[101:end]
xaxis_ = [xaxis_; fill(r, length(preiterated))]
yaxis_ = [yaxis_; preiterated]
end
return xaxis_, yaxis_
end

@time xaxis, yaxis = bifurcation(0.01)
diagram = scatter(xaxis, yaxis, markersize = 0.1, markeralpha = 0.5, legend = false, xlims = (1, 4), ylims = (0, 1))
@time xaxis, yaxis = bifurcation(0.001)
diagram = scatter(xaxis, yaxis, markersize = 0.1, markeralpha = 0.5, legend = false, xlims = (1, 4), ylims = (0, 1));
png(diagram, "bifurcation.png")

 -->
![](https://hackmd.io/_uploads/ryBWLc-6h.png)