text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 0053bace-00d4-11eb-1072-cbc1284550c5
# 環境設定
begin
	using Printf
	using Dates
	using PlutoUI
	using Plots
	# using LinearAlgebra # Most of them are in Base
	using RowEchelon # 單元 9．Column Correspondence Theorem

	_date_ = today()
	_time_ = Dates.format(now(), "HH:MM")
	md"""
	## 環境配置
	今日： $(_date_) \
	此時： $(_time_) \
	"""
end
#=
Space before $ does matter ... _O_
[Line breaks in markdown]
(https://gist.github.com/shaunlebron/746476e6e7a4d698b373)
=#

# ╔═╡ 9d9a54ea-ff9c-11ea-1db7-1ba194b9fb3c
md"""
# 線性代數 + Julia + $$\LaTeX$$ 的學習筆記
版本日期：$_date_

整個學習過程將以如下「線性代數」課程為主軸：

### 線性代數 台灣大學電機系 蘇柏青

本課程是線性代數的入門課程。線性代數係以「向量空間」(Vector Space)為核心概念之數學工具，擁有極廣泛之應用，非常值得理工商管等科系大學部同學深入修習，作為日後專業應用之基礎。

課程來源：[http://ocw.aca.ntu.edu.tw/ntu-ocw/index.php/ocw/cou/102S207](http://ocw.aca.ntu.edu.tw/ntu-ocw/index.php/ocw/cou/102S207)

## 學習目標
如下為幾個學習的子目標：

### 學科
- 線性代數 - 重新學習，了解重要概念的中文及英文詞𢑥及應用。

### 工具

- Julia - 深入學習，了解重要套件的應用及使用。
- Pluto - 隨之成長，作為撰寫學習記錄的工具。
- LaTeX - 隨緣學習，作為撰寫學習記錄的工具。
- Markdown - 隨緣學習，作為撰寫學習記錄的工具。

### 服務
- GitHub - 隨緣學習，用來記錄整個學習歷程及分享學習內容。
"""

# ╔═╡ a40a1be4-ff9c-11ea-33d3-8bcf3b73e930
md"## 單元 1．Basic Concepts on Matrices and Vectors"

# ╔═╡ 71136776-fff2-11ea-0eb6-47a04c1c77d6
md"""
### Matrix
$$A=\begin{bmatrix} a_{11} & \dots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{m1} & \dots & a_{mn}\end{bmatrix}=\begin{bmatrix} a_{ij} \end{bmatrix}=M_{mn}$$
"""

# ╔═╡ a156894a-ffd4-11ea-1e22-29fb46b20097
md"""
### Matrix Addition
$$\begin{bmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{bmatrix} + \begin{bmatrix} 1 & 1 \\ 1 & 1 \\ 1 & 2 \end{bmatrix} = \begin{bmatrix} 2 & 3 \\ 4 & 5 \\ 6 & 8 \end{bmatrix}$$
"""

# ╔═╡ 06754e3a-ffe2-11ea-0220-9b066f8651e7
[1 2; 3 4; 5 6]+[1 1; 1 1; 1 2]

# ╔═╡ c50b8e9a-ffe2-11ea-2326-c3b69913b27c
md"""
### Scalar Multiplication
$$cA$$
$$3 \cdot \begin{bmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{bmatrix}$$
"""

# ╔═╡ 04823ca4-ffe3-11ea-3fab-0d8eaedb8c9a
3 * [1 2; 3 4; 5 6]

# ╔═╡ b7bcfc2c-ffe4-11ea-2862-432081314a58
3 .* [1 2; 3 4; 5 6]

# ╔═╡ f4077336-ffe4-11ea-2482-81f261e32798
md"""
### Transpose
$$C=\begin{bmatrix} 7 & 9 \\ 18 & 31 \\ 52 & 68 \end{bmatrix} \;\:\Rightarrow\;\; C^T = \begin{bmatrix} 7 & 18 & 52 \\ 9 & 31 & 68 \end{bmatrix}$$
"""

# ╔═╡ 047d2898-ffe5-11ea-241f-6d1064a665e8
let
	C=[7 9; 18 31; 52 68]
	C'
end

# ╔═╡ ea9c0a1c-ffe8-11ea-3df9-c5a82dc5c077
md"""
### Vectors

Row Vector:

$$\begin{bmatrix} 1 & 2 & 3 & 4 \end{bmatrix}$$

Column Vector:

$$\begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix}$$

$$\Downarrow$$

$$\begin{bmatrix} 1 & 2 & 3 & 4 \end{bmatrix}^T$$

The ith componet of $$\textbf{v}$$

$$v_i$$

"""

# ╔═╡ 5018d0f0-ffe9-11ea-0ba4-f3ca17a7eddb
[ 1 2 3 4]

# ╔═╡ 676ff59e-ffe9-11ea-22ec-31be78dfa211
[1; 2; 3; 4;]

# ╔═╡ 491d2fd0-ffeb-11ea-009f-4b10a02ba938
[ 1 2 3 4]'

# ╔═╡ 47e73308-fff9-11ea-2cee-7f050938fe19
md"""
### Linear Combination
A $$\textcolor{red}{linear\;combination}\;of\;vectors \;
$$\textbf{u}_1,\textbf{u}_2,\dots,\textbf{u}_k$$ is a vector of the form

$$c_1\textbf{u}_1+c_2\textbf{u}_2+\dots+c_k\textbf{u}_k$$

where $$c_1,c_2,\dots,c_k$$ are scalars. These scalars are called the $$\textcolor{red}{coefficients}$$ of the linear combination.
"""

# ╔═╡ 99e6bccc-fffa-11ea-223a-a3dfd65dee92
md"""
### Standard Vectors

The standard vectors of $$R^n$$ are defined as

$$e_1=\begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}, e_2=\begin{bmatrix} 0 \\ 1 \\ \vdots \\ 0 \end{bmatrix}, \dots , e_n=\begin{bmatrix} 0 \\ 0 \\ \vdots \\ 1 \end{bmatrix}$$

"""

# ╔═╡ 741db230-fffb-11ea-23aa-63d365767cbd
md"""
### Matrix-Vector Product
$$Av=v_1a_1+v_2a_2+\dots+v_na_n$$
"""

# ╔═╡ c497553e-fffb-11ea-3559-e59103e843e2
let
	A=[1 2; 3 4; 5 6]
	v=[7;8]
	A*v
end

# ╔═╡ aff2680a-fffe-11ea-0494-1d349105eff0
md"""
### Identity Matrix
$$I_3=\begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1\end{bmatrix}$$

"""

# ╔═╡ 42f11192-ffff-11ea-34ca-81545f4b5374
md"""
### Stochastic Matrix
$$A=\begin{bmatrix} 0.85 & 0.03 \\ 0.15 & 0.97 \end{bmatrix}$$
"""

# ╔═╡ 988b781e-026f-11eb-02b7-8d2de841d478
begin
	u01xslider = @bind u01x Slider(1:100; default=40, show_value=true)
	u01cslider = @bind u01c Slider(0:1200; default=1200, show_value=true)
	u01sslider = @bind u01s Slider(0:1200; default=0, show_value=true)
	md"""
	Slide to set number of **years**: $(u01xslider)

	Slide to set population of **city**: $(u01cslider)

	Slide to set population of **suburban**: $(u01sslider)
	"""
end

# ╔═╡ cb2c5e54-ffff-11ea-3e4e-a349d0a09341
let
	x=u01x # Number of Years (x)
	pc=u01c # Population of City
	ps=u01s # Population of Suburban
	A=[0.85 0.03; 0.15 0.97]
	#=
	# p0 Population in year 0
	p0=[500; 700]
	p1=A*p0
	p2=A*(p1)
	p3=A*(p2)
	p4=A*(p3)
	p5=A*(p4)
	x=0:5
	Y=hcat(p0, p1, p2, p3, p4, p5)
	plot(x, Y', title = "Population", label = ["City" "Suburban"])
	=#
	p=[pc; ps]
	Y=p
	for i in 1:x
		p=A*p
		Y=hcat(Y, p)
	end
	plot(0:x, Y', title = "Population Trend", label = ["City" "Suburban"])
end

# ╔═╡ 67ed416e-00d1-11eb-29cc-23eaa30c3446
html"""

"""

# ╔═╡ afcd29d0-ff9c-11ea-09bf-e5a0ed087c42
md"## 單元 2．System of Linear Equations"

# ╔═╡ 8c300b8a-ffd5-11ea-3ec9-27ccc82cea6e
# \backslash \setminus
md"""
### System of Linear Equations
$$A=\begin{bmatrix} 1 & -2 & -1 \\ 3 & -6 & -5 \\ 2 & -1 & 1 \end{bmatrix}
\;\;
b=\begin{bmatrix} 3 \\ 3\\ 0 \end{bmatrix}$$
$$Ax=b$$

Solves Ax = b by (essentially) Gaussian elimination (Julia \ Operator):

$$x=A\;\backslash\;b$$
"""

# ╔═╡ 70803966-ff9c-11ea-0a8a-d7860f4f87c0
# Solve System of Linear Equations
let
	with_terminal() do
		A=[1 -2 -1; 3 -6 -5; 2 -1 1]
		b=[3; 3; 0]
		x=A \ b
		println(x)
		println(A*x)
		println(A*x == b)
	end
end

# ╔═╡ 6a110912-0009-11eb-0a73-0b2e2ec0b634
md"""
### Row Echelon Form & Reduced Row Echelon Form
"""

# ╔═╡ 4b8f75e0-000a-11eb-175b-d7a9102f65ca
# Solve System of Linear Equations
let
	with_terminal() do
		A=[1 -3 0 2 0; 0 0 1 6 0; 0 0 0 0 1; 0 0 0 0 0]
		b=[7; 9; 2; 0]
		x=A \ b
		println(x)
		println(A*x)
		println(A*x == b)
	end
end

# ╔═╡ 1f9799ea-00d2-11eb-1b94-d7d3056a2c41
html"""

"""

# ╔═╡ ecfe3ef0-007a-11eb-3efb-7740876167a8
md"""
## 單元 3．Gaussian Elimination

實作參考：

[Gaussian-elimination.pdf](http://web.mit.edu/18.06/www/Spring17/Gaussian-elimination.pdf)

[Numerical Analysis by Julia Series 1 — Gauss Elimination | by Treee July | Medium]
(https://medium.com/@julytreee/numerical-analysis-by-julia-series-1-gauss-elimination-68e902a43c7e)
"""

# ╔═╡ b76720a6-00d2-11eb-19ae-d59e7204a005
md"""
對列及行的參照：
"""

# ╔═╡ 0c9dd736-00d8-11eb-077a-67cf34db152d
let
	A=[ 1 2 3; 4 5 6; 7 8 9]
end

# ╔═╡ 2fa8b248-00d6-11eb-0f85-871512ef16d2
let
	o=[]
	# Matrix
	A=[ 1 2 3; 4 5 6; 7 8 9]
	push!(o, @sprintf("A: %s", A))
	# Elements
	push!(o, @sprintf("A[1, 1]: %s", A[1, 1]))
	push!(o, @sprintf("A[end, end]: %s", A[end, end]))
	# Row
	r1=A[1,:]
	push!(o, @sprintf("r1: %s", r1))
	∑Ai=A[1,:]+A[2,:]+A[3,:]
	push!(o, @sprintf("∑Ai: %s", ∑Ai))
	# Column
	c1=A[:,1]
	push!(o, @sprintf("c1: %s", c1))
	∑Aj=A[:,1]+A[:,2]+A[:,3]
	push!(o, @sprintf("∑Aj: %s", ∑Aj))
	with_terminal(dump, o)
end

# ╔═╡ 479e95e2-0274-11eb-2c2a-993ba9119fba
let
	with_terminal() do
		println("👍 選定之輸出方案： 1) 容易以 do ... end 區塊包裝 3) 轉貼程式碼到他處不用修改\n")
		# Get Current Time
		command=`date`
		run(command)
		# Matrix
		A=[ 1 2 3; 4 5 6; 7 8 9]
		print("A: "); dump(A)
		# Elements
		print("A[1, 1]: "); dump(A[1, 1])
		print("A[end, end]: "); dump(A[end, end])
		# Row
		r1=A[1,:]
		print("r1: ");	dump(r1)
		∑Ai=A[1,:]+A[2,:]+A[3,:]
		print("∑Ai: "); dump(∑Ai)
		# Column
		c1=A[:,1]
		print("c1: ");	dump(c1)
		∑Aj=A[:,1]+A[:,2]+A[:,3]
		print("∑Aj: "); dump(∑Aj)
		println()
		run(`cal -h`)
	end
end

# ╔═╡ 816cc338-03a1-11eb-2872-c7643faab770
let
	Text() do io
		# Matrix
		A=[ 1 2 3; 4 5 6; 7 8 9]
		print(io, "A: "); dump(io, A)
		# Elements
		print(io, "A[1, 1]: "); dump(io, A[1, 1])
		print(io, "A[end, end]: "); dump(io, A[end, end])
		# Row
		r1=A[1,:]
		print(io, "r1: ");	dump(io, r1)
		∑Ai=A[1,:]+A[2,:]+A[3,:]
		print(io, "∑Ai: "); dump(io, ∑Ai)
		# Column
		c1=A[:,1]
		print(io, "c1: ");	dump(io, c1)
		∑Aj=A[:,1]+A[:,2]+A[:,3]
		print(io, "∑Aj: "); dump(io, ∑Aj)
	end
end

# ╔═╡ ee3ed802-03a3-11eb-025c-5bee178511af
let
		# Matrix
		A=[ 1 2 3; 4 5 6; 7 8 9]
		# Row
		r1=A[1,:]
		∑Ai=A[1,:]+A[2,:]+A[3,:]
		# Column
		c1=A[:,1]
		∑Aj=A[:,1]+A[:,2]+A[:,3]
	md"""
	A: $(Text(A))

	A[1, 1]: $(Text(A[1, 1]))

	A[end, end]: $(Text(A[end, end]))

	r1: $(Text(r1))

	∑Ai: $(Text(∑Ai))

	c1: $(Text(c1))

	∑Aj: $(Text(∑Aj))
	"""
end

# ╔═╡ f7de8296-0393-11eb-0782-57e9272ad999
md"""
## 單元 4．The language of set theory
"""

# ╔═╡ 35ec879e-03ae-11eb-15ce-016ae8a89168
md"""
#### Subset
$$Let\; S_1=\{a, b, c, d, e\},\; S_2={a, b, e}$$$$
S_2 ⊂ S_1\;means\;$$$$
∀x ∈ S_2,\:x\;is\;also ∈ S_1.$$
"""

# ╔═╡ cfaeab08-03b2-11eb-3930-ddb92e799789
let
	with_terminal() do
		s1=Set(["a", "b", "c", "d", "e"])
		println("🖐 Give Me Five, 原來集合在 Julia 裹頭是長這樣子喔！☯")
		dump(s1)
		println(s1)
		s2=Set(["a", "b", "e"])
		println(s2)
		# subset
		println(⊆(s2, s1))
		# union set
		println(∪(s1, s2))
		# intersection set
		println(∩(s1, s2))
		# difference set
		println(setdiff(s1, s2))
		println(setdiff(s2, s1))
	end
end

# ╔═╡ 0445c5ba-03cb-11eb-2093-b71bae438a19
md"""
## 單元 5．Span of a Set of Vectors
"""

# ╔═╡ 698fc53a-03d5-11eb-39b7-a776fabae2f2
md"""
實作參考：

[Linear Algebra – Quantitative Economics with Julia](https://julia.quantecon.org/tools_and_techniques/linear_algebra.html)
"""

# ╔═╡ 408f674e-03d5-11eb-3dc8-ed1addf002be
let
	x_vals = [0 0 0 ; 2 -3 -4]
	y_vals = [0 0 0 ; 4 3 -3.5]

	plot(x_vals, y_vals, arrow = true, color = :blue,
	     legend = :none, xlims = (-5, 5), ylims = (-5, 5),
	     annotations = [(2.2, 4.4, "[2, 4]"),
	                    (-3.3, 3.3, "[-3, 3]"),
	                    (-4.4, -3.85, "[-4, -3.5]")],
	     xticks = -5:1:5, yticks = -5:1:5,
	     framestyle = :origin)
end

# ╔═╡ 73360416-03da-11eb-1721-cf76aa462151
begin
	u05islider = @bind u05i Slider(1:5; default=1, show_value=true)
	u05jslider = @bind u05j Slider(1:5; default=5, show_value=true)
	md"""
	Slide to set **i**: $(u05islider)

	Slide to set **j**: $(u05jslider)
	"""
end

# ╔═╡ ac04f844-03d8-11eb-2db2-dbc19b343325
let
	i=u05i
	j=u05j
	# fixed linear function, to generate a plane
	f(x, y) = 0.2x + 0.1y

	# lines to vectors
	x_vec = [0 0 0 0 0; 3 3 -4 -4 3.5]
	y_vec = [0 0 0 0 0; 4 -4 -4 4 0]
	z_vec = [0 0 0 0 0; f(3, 4) f(3, -4) f(-4,-4) f(-4, 4) f(3.5, 0)]
	color = [:blue :green :red :brown :cyan]

	# draw the plane
	n = 20
	grid = range(-5, 5, length = n)
	z2 = [ f(grid[row], grid[col]) for row in 1:n, col in 1:n ]
	# wireframe(grid, grid, z2, fill = :blues, gridalpha =1 )
	plot(grid, grid, z2, fill = :blues, gridalpha = 1, lindwidth = 0.5, seriestype = :wireframe)
	# plot(grid, grid, z2, fill = :blues, gridalpha = 1, lindwidth = 0.5, seriestype = :surface)
	# Dots
	# plot!([0; 4; 4; -4; -4], [0; 4; -4; 4; -4], [-1.5; -1.5; -1.5; -1.5; -1.5], labels = "", seriestype = :scatter3d)
	p = [ 0 0 -1.5; 4 4 -1.5; 4 -4 -1.5; -4 4 -1.5; -4 -4 -1.5]' # Transpose
	plot!(p[1, i:j], p[2, i:j], p[3, i:j], labels = "", seriestype = :scatter3d)
	plot!(p[1, i:j], p[2, i:j], p[3, i:j], labels = "", seriestype = :path3d)
	# Vectors
	plot!(x_vec[:, i:j], y_vec[:, i:j], z_vec[:, i:j], color = color[:,i:j], linewidth = 3, xlabel = "x", ylabel = "y", zlabel = "z", labels = "", colorbar = false)
	# plot!(x_vec, y_vec, z_vec, color = color[:,i:j], linewidth = 3, xlabel = "x", ylabel = "y", zlabel = "z", labels = "", colorbar = false)
end

# ╔═╡ 8133c7e0-03d8-11eb-0697-afe2082ed481
md"""
### 單元 6．Linear Dependence and Linear Independence
"""

# ╔═╡ 82e1b674-03df-11eb-26e4-6f09a0a832d8
let
	with_terminal() do
		A=[1 2 -1; -1 1 -8; 2 -1 13; 1 -1 8]
		b=[0; 1; -2; 1]
		x=A \ b
		println(x)
		println(A*x)
		println(A*x == b)
	end
end

# ╔═╡ 55de62fc-03e8-11eb-2ade-6500eacc11c8
let
	with_terminal() do
		A=[1 -4 2 -1 2; 2 -8 3 2 1]
		b=[0; 0]
		x=A \b
		println(x)
		println(A*x)
		println(A*x == b)
	end
end

# ╔═╡ 69edd6d6-0469-11eb-1a3a-69b7f222541c
md"""
## 單元 7．Matrix Multiplication
"""

# ╔═╡ 4e3a788a-04b5-11eb-3574-e53d1db76173
md"""
### Matrix Multiplication
$$\,
\begin{align*}
Let\;v, x, y &∈ R^n. Suppose\;A\;and\;B\;are\;n × n\;matrices. \\
x &= Bv \\
y &= Ax \\
&\Downarrow\\
y &= Cv = A(Bv) = (AB)v
\end{align*}
\,$$
"""

# ╔═╡ c95c00b4-046b-11eb-09e0-69068c23a0e0
let
	with_terminal() do
		A=[1 2; 1 1]
		B = [ 1 1; 1 -1]
		dump(A * B)
		dump(B * A)
	end
end

# ╔═╡ a7edc48a-046e-11eb-2920-73dd6ab80b73
md"""
## 單元 8．Invertibility and Elmentary Matrices
"""

# ╔═╡ b6f85278-046e-11eb-25db-e71ea510fe6a
md"""
### Inverse
$$\,
\begin{align*}
&\text{An n × n matrix A is called invertible if there exists an n × n matrix B such
that }\\
\\
&AB = BA = I_n.\\
\\
&\text{In this case, B is called an inverse of A.}
\end{align*}
\,$$
"""

# ╔═╡ c42784e8-0483-11eb-2105-41108f06910e
# use of & to state the alignment points.
# Why does \\ not return a new line in an equation? - TeX - LaTeX Stack Exchange
# https://tex.stackexchange.com/questions/194236/why-does-not-return-a-new-line-in-an-equation
md"""
$$\,
\begin{matrix}
\begin{align*}
A&=\begin{bmatrix} 1 & 2 \\ 3 & 5 \end{bmatrix}
\\
B&=\begin{bmatrix} -5 & 2 \\ 3 & -1 \end{bmatrix}
\end{align*}
&
\Rightarrow
&
\begin{split}
AB=
\begin{bmatrix} 1 & 2 \\ 3 & 5 \end{bmatrix}
\begin{bmatrix} -5 & 2 \\ 3 & -1 \end{bmatrix}=\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}=I_2
\\
BA=
\begin{bmatrix} -5 & 2 \\ 3 & -1 \end{bmatrix}
\begin{bmatrix} 1 & 2 \\ 3 & 5 \end{bmatrix}=\begin{bmatrix} 1 & 0 \\ 0 & 1\end{bmatrix}=I_2
\end{split}
\end{matrix}
\,$$
"""

# ╔═╡ 6a0af07c-0482-11eb-287b-8d97045e4f36
let
	with_terminal() do
		A=[1 2; 3 5]
		B=inv(A)
		dump(A)
		dump(B)
		dump(round.(B))
		dump.(round.(A*B))
		dump.(round.(B*A))
	end
end

# ╔═╡ 7cadabce-0483-11eb-1a22-d7fcb40f6fe9
md"""
疑：為什麼只是近似的值？因為 3?
"""

# ╔═╡ 061630b8-0496-11eb-1646-137d406b7c37
md"""
## 單元 9．Column Correspondence Theorem
"""

# ╔═╡ d32ddc0a-0497-11eb-14e8-81331c56ebd8
md"""
實作參考：

[blegat/RowEchelon.jl: Small package containing the rref fonction for computing the reduced row echelon form of the matrix A]
(https://github.com/blegat/RowEchelon.jl)
> 可做為 module 及 RREF 計算實作參考。
> 提供的 function 為
> rref, rref!
> rref_with_pivots, rref_with_pivots!,
"""

# ╔═╡ 5f220836-049a-11eb-389c-b75fa5bb21d7
md"""
### Reduced Row Echelon Form (RREF)
"""

# ╔═╡ 86827bfe-049a-11eb-241c-cbbefed2bd08
let
	with_terminal() do
		println("Using RREF:")
		A=[ 1 2 -1 2 1 2; -1 -2 1 2 3 6; 2 4 -3 2 0 3; -3 -6 2 0 3 9]
		dump(A)
		println(Text(A))
		B=rref(A)
		dump(round.(B))
		println(Text(round.(B)))

		println("Using \\:")
		A1=A[:, 1:5]
		b1=A[:, 6]
		x=A1 \ b1
		dump(round.(x))
		println(Text(round.(x)))
		println(Text(round.(A1*x)))

		println("Using \\ with RREF:")
		A2=B[:, 1:5]
		b2=B[:, 6]
		y=A2 \ b2
		dump(round.(y))
		println(Text(round.(y)))
		println(Text(round.(A2*y)))
	end
end

# ╔═╡ 5db63852-04ac-11eb-09e2-b3a09a28fe5a
md"""
## 單元 10．The Inverse of a Matrix
"""

# ╔═╡ 68240db4-04ac-11eb-3e83-819f3efb8c66
md"""
### Matrix Inversion
"""

# ╔═╡ 2ae45a8a-04b1-11eb-195e-1399a9b27328
md"""
$$\,
\begin{bmatrix} \color{blue}{A} & I_3 \end{bmatrix}=
\begin{bmatrix}
\begin{array}{ccc|ccc}1 & 2 & 3 & \color{red}{1} & 0 & 0 \\ 2 & 5 & 6 & 0 & \color{red}{1} & 0\\ 3 & 4 & 8 & 0 & 0 & \color{red}{1}\end{array}
\end{bmatrix}
\rightarrow
\begin{bmatrix}
\begin{array}{ccc|ccc}\color{red}{1} & 0 & 0 & -16 & 4 & 3 \\ 0 & \color{red}{1} & 0 & -2 & 1 & 0\\ 0 & 0 & \color{red}{1} & 7 & -2 & -1\end{array}
\end{bmatrix}
=
\begin{bmatrix} I_3 & \color{blue}{B}\end{bmatrix}
\,$$
"""

# ╔═╡ 9642ea30-04ac-11eb-0a86-b7d63a7594d7
let
	with_terminal() do
		AI=[ 1 2 3 1 0 0; 2 5 6 0 1 0; 3 4 8 0 0 1]
		println(AI)
		IB=round.(rref(AI))
		println(IB)
		C=IB[:, 4:6]
		println(C)
		D=round.(inv(AI[:, 1:3]))
		println(D)
	end
end

# ╔═╡ 97f12ddc-00cf-11eb-173d-47a7931a9a08
md"""
### <<<
"""

# ╔═╡ 4a9f304c-0543-11eb-3b96-5fb9fb043444
md"""
## 附錄：Pluto 遊樂場/Playground
"""

# ╔═╡ c13b4d10-0545-11eb-178c-9fbfe174d4a6
md"""
### 注意 \$\$ 及 \$
"""

# ╔═╡ 5367d0b2-0543-11eb-3810-b31a6e903242
let
	x=1
	y=2
	z=3
	md"""
	#### No space before $
	x:$x \
	y:$y \
	z:$z
	"""
end

# ╔═╡ 1e622640-0545-11eb-20fc-29bed94aac0c
let
	x=1
	y=2
	z=3
	md"""
	#### With space before $
	x: $x \
	y: $y \
	z: $z
	"""
end

# ╔═╡ b3f4edae-054b-11eb-02b5-199ce15d9620
let
	_test_ = "x^2"
	# 在表格中沒有變數替換，有公式替換
	md"""
	### 測試 \$\$ 及 \$ 混用
	|變數|數值|
	|:---|:---|
	|今日| $(_date_) |
	|其他| $(_date_) |
	|函數| $$x^3$$ $(_test_)  |
	"""
	# 只好用這個版本，先公式替換然後變數替換
	md"""
	## 測試 \$\$ 及 \$ 混用
	今日： $(_date_) \
	其他： $(_date_) \
	函數： $$x^3$$ $(_test_)
	"""
end

# ╔═╡ e4b8f094-04c8-11eb-3bf7-93bd85db2dc9
md"""
## 附錄：$$\LaTeX$$ 遊樂場/Playground
"""

# ╔═╡ 5897b63a-04c9-11eb-34ab-b908c731d3d5
#=
\begin{equation*}
multline, gather, align, aligned[t|b|c], alignat{9}, align*, split
\end{equation*}
=#
md"""
實作參考：

[User’s Guide for the amsmath Package]
(https://www.latex-project.org/help/documentation/amsldoc.pdf)

$$\,
\begin{multline}
a+b+c+d+e+f\\
+i+j+k+l+m+n
\end{multline}
$$$$
\begin{gather}
a_1=b_1+c_1\\
a_2=b_2+c_2-d_2+e_2
\end{gather}
$$$$
\begin{align}
a_1& =b_1+c_1\\
a_2& =b_2+c_2-d_2+e_2
\end{align}
$$$$
\left.
\begin{aligned}[b]
&{\color{red} \text{Text in red}} \\
&{\color{blue} \text{Text in blue}}\\
&\text{Text with equation $a_{11}, a_{12},\dots$}\\
&A=πr^2\\
&c^2=a^2+b^2\\
&x
\end{aligned}
\big\}\quad\text{The formula}
\Bigg( 1+2 \Bigg]
F = G \left( \frac{m_1 m_2}{r^2} \right)
\right\} \text{This is it!}
$$$$
\begin{align}
x&=y & X&=Y & a&=b+c\\
x’&=y’ & X’&=Y’ & a’&=b\\
x+x’&=y+y’ & X+X’&=Y+Y’ & a’b&=c’b
\end{align}
$$$$
\begin{alignat}{2}
x& = y_1-y_2+y_3-y_5+y_8-\dots
&\quad& \text{by \eqref{eq:C}}\\
& = y’\circ y^* && \text{by \eqref{eq:D}}\\
& = y(0) y’ && \text {by Axiom 1.}
\end{alignat}
$$$$
\left.\begin{aligned}
B’&=-\partial\times E,\\
E’&=\partial\times B - 4\pi j,
\end{aligned}
\right\}
\qquad \text{Maxwell’s equations}
$$$$
P_{r-j}=\begin{cases}
0& \text{if $r-j$ is odd},\\
r!\,(-1)^{(r-j)/2}& \text{if $r-j$ is even}.
\end{cases}
\,$$
"""

# ╔═╡ 51c9222c-0518-11eb-3827-8fcf8a59b0fd
md"""
## 附錄：Markdown 遊樂場/Playground
"""

# ╔═╡ 602f56ce-0518-11eb-27d3-7fde99127b36
md"""
實作參考：

[Markdown · The Julia Language]
(https://docs.julialang.org/en/v1/stdlib/Markdown/)

```math
f(a) = \frac{1}{2\pi}\int_{0}^{2\pi} (\alpha+R\cos(\theta))d\theta
```

A paragraph containing some ``\LaTeX`` markup.
```
Ax=b
```
```math
Ax=b
```
Some Markdown text with some *blue* text.

註：[How to apply color in Markdown? - Stack Overflow]
(https://stackoverflow.com/questions/35465557/how-to-apply-color-in-markdown)
"""

# ╔═╡ 7dbac3ec-00d1-11eb-2c17-3bab40ffaa2e
html"""

"""

# ╔═╡ 770ab5ac-fff8-11ea-1ed1-87b0ae3aca70
md"""
## 參考資料
### Linear Algebra

[ ] [線性代數 - 臺大開放式課程 (NTU OpenCourseWare)]
(http://ocw.aca.ntu.edu.tw/ntu-ocw/index.php/ocw/cou/102S207/3)

[ ] [Introduction to Applied Linear Algebra – Vectors, Matrices, and Least Squares]
(http://vmls-book.stanford.edu/)
> [Julia language companion](http://vmls-book.stanford.edu/vmls-julia-companion.pdf)

### Julia

[ ] [Introduction to Julia]
(https://juliaacademy.com/courses/enrolled/375479)

> [ ] Advanced topics

[ ] [Julia for Data Science]
(https://juliaacademy.com/courses/enrolled/937702)

[ ] [18.S191 Introduction to Computational Thinking]
(https://computationalthinking.mit.edu/Fall20/)

[Linear Algebra – Quantitative Economics with Julia]
(https://julia.quantecon.org/tools_and_techniques/linear_algebra.html)

> [QuantEcon.cheatsheet/julia-cheatsheet.pdf](https://github.com/QuantEcon/QuantEcon.cheatsheet/blob/master/julia/julia-cheatsheet.pdf)

[cheatsheets/plotsjl-cheatsheet.pdf]
(https://github.com/sswatson/cheatsheets/blob/master/plotsjl-cheatsheet.pdf)

[Unicode Input · The Julia Language]
(https://docs.julialang.org/en/v1/manual/unicode-input/)

[Visualizing Graphs in Julia using Plots and PlotRecipes – Tom Breloff]
(http://www.breloff.com/Graphs/)

### Pluto
[fonsp/Pluto.jl: 🎈 Simple reactive notebooks for Julia]
(https://github.com/fonsp/Pluto.jl)

[Docstrings · PlutoUI.jl]
(https://juliahub.com/docs/PlutoUI/abXFp/0.6.3/autodocs/)

### $$\LaTeX$$

[LaTeX Documentation]
(https://www.latex-project.org/help/documentation/)

[User’s Guide for the amsmath Package]
(https://www.latex-project.org/help/documentation/amsldoc.pdf)

[LaTeX syntax · Documenter.jl]
(https://juliadocs.github.io/Documenter.jl/v0.7/man/latex.html)

[Documentation - Overleaf, Online LaTeX Editor]
(https://www.overleaf.com/learn/latex/Main_Page)

[LaTeX - Mathematical Python]
(https://www.math.ubc.ca/~pwalls/math-python/jupyter/latex/)

[LaTeX help 1.1 - Table of Contents]
(http://www.emerson.emory.edu/services/latex/latex_toc.html)

[List of mathematical symbols - Wikiwand]
(https://www.wikiwand.com/en/List_of_mathematical_symbols)

### Markdown
[Markdown Cheatsheet · adam-p/markdown-here Wiki]
(https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet)

[Markdown · The Julia Language]
(https://docs.julialang.org/en/v1/stdlib/Markdown/)

### GitHub

[ ] [Hello World · GitHub Guides]
(https://guides.github.com/activities/hello-world/)

### 其他
[三度辭典網 > 術語中英雙語詞典]
(https://www.3du.tw/term/)
"""

# ╔═╡ b5016eda-0522-11eb-1d8b-37135a22c7bf
md"""
## 感謝建議
"""

# ╔═╡ e9001c9a-0522-11eb-0d63-f718249e03ec

# ╔═╡ Cell order:
# ╟─9d9a54ea-ff9c-11ea-1db7-1ba194b9fb3c
# ╟─0053bace-00d4-11eb-1072-cbc1284550c5
# ╟─a40a1be4-ff9c-11ea-33d3-8bcf3b73e930
# ╟─71136776-fff2-11ea-0eb6-47a04c1c77d6
# ╟─a156894a-ffd4-11ea-1e22-29fb46b20097
# ╠═06754e3a-ffe2-11ea-0220-9b066f8651e7
# ╟─c50b8e9a-ffe2-11ea-2326-c3b69913b27c
# ╠═04823ca4-ffe3-11ea-3fab-0d8eaedb8c9a
# ╠═b7bcfc2c-ffe4-11ea-2862-432081314a58
# ╟─f4077336-ffe4-11ea-2482-81f261e32798
# ╠═047d2898-ffe5-11ea-241f-6d1064a665e8
# ╟─ea9c0a1c-ffe8-11ea-3df9-c5a82dc5c077
# ╠═5018d0f0-ffe9-11ea-0ba4-f3ca17a7eddb
# ╠═676ff59e-ffe9-11ea-22ec-31be78dfa211
# ╠═491d2fd0-ffeb-11ea-009f-4b10a02ba938
# ╟─47e73308-fff9-11ea-2cee-7f050938fe19
# ╟─99e6bccc-fffa-11ea-223a-a3dfd65dee92
# ╟─741db230-fffb-11ea-23aa-63d365767cbd
# ╠═c497553e-fffb-11ea-3559-e59103e843e2
# ╟─aff2680a-fffe-11ea-0494-1d349105eff0
# ╟─42f11192-ffff-11ea-34ca-81545f4b5374
# ╟─988b781e-026f-11eb-02b7-8d2de841d478
# ╠═cb2c5e54-ffff-11ea-3e4e-a349d0a09341
# ╟─67ed416e-00d1-11eb-29cc-23eaa30c3446
# ╟─afcd29d0-ff9c-11ea-09bf-e5a0ed087c42
# ╟─8c300b8a-ffd5-11ea-3ec9-27ccc82cea6e
# ╠═70803966-ff9c-11ea-0a8a-d7860f4f87c0
# ╟─6a110912-0009-11eb-0a73-0b2e2ec0b634
# ╠═4b8f75e0-000a-11eb-175b-d7a9102f65ca
# ╟─1f9799ea-00d2-11eb-1b94-d7d3056a2c41
# ╟─ecfe3ef0-007a-11eb-3efb-7740876167a8
# ╟─b76720a6-00d2-11eb-19ae-d59e7204a005
# ╠═0c9dd736-00d8-11eb-077a-67cf34db152d
# ╠═2fa8b248-00d6-11eb-0f85-871512ef16d2
# ╠═479e95e2-0274-11eb-2c2a-993ba9119fba
# ╠═816cc338-03a1-11eb-2872-c7643faab770
# ╠═ee3ed802-03a3-11eb-025c-5bee178511af
# ╟─f7de8296-0393-11eb-0782-57e9272ad999
# ╟─35ec879e-03ae-11eb-15ce-016ae8a89168
# ╠═cfaeab08-03b2-11eb-3930-ddb92e799789
# ╠═0445c5ba-03cb-11eb-2093-b71bae438a19
# ╠═698fc53a-03d5-11eb-39b7-a776fabae2f2
# ╠═408f674e-03d5-11eb-3dc8-ed1addf002be
# ╠═73360416-03da-11eb-1721-cf76aa462151
# ╠═ac04f844-03d8-11eb-2db2-dbc19b343325
# ╟─8133c7e0-03d8-11eb-0697-afe2082ed481
# ╠═82e1b674-03df-11eb-26e4-6f09a0a832d8
# ╠═55de62fc-03e8-11eb-2ade-6500eacc11c8
# ╟─69edd6d6-0469-11eb-1a3a-69b7f222541c
# ╠═4e3a788a-04b5-11eb-3574-e53d1db76173
# ╠═c95c00b4-046b-11eb-09e0-69068c23a0e0
# ╟─a7edc48a-046e-11eb-2920-73dd6ab80b73
# ╟─b6f85278-046e-11eb-25db-e71ea510fe6a
# ╟─c42784e8-0483-11eb-2105-41108f06910e
# ╠═6a0af07c-0482-11eb-287b-8d97045e4f36
# ╟─7cadabce-0483-11eb-1a22-d7fcb40f6fe9
# ╟─061630b8-0496-11eb-1646-137d406b7c37
# ╠═d32ddc0a-0497-11eb-14e8-81331c56ebd8
# ╟─5f220836-049a-11eb-389c-b75fa5bb21d7
# ╠═86827bfe-049a-11eb-241c-cbbefed2bd08
# ╟─5db63852-04ac-11eb-09e2-b3a09a28fe5a
# ╟─68240db4-04ac-11eb-3e83-819f3efb8c66
# ╟─2ae45a8a-04b1-11eb-195e-1399a9b27328
# ╠═9642ea30-04ac-11eb-0a86-b7d63a7594d7
# ╠═97f12ddc-00cf-11eb-173d-47a7931a9a08
# ╟─4a9f304c-0543-11eb-3b96-5fb9fb043444
# ╟─c13b4d10-0545-11eb-178c-9fbfe174d4a6
# ╠═5367d0b2-0543-11eb-3810-b31a6e903242
# ╠═1e622640-0545-11eb-20fc-29bed94aac0c
# ╠═b3f4edae-054b-11eb-02b5-199ce15d9620
# ╟─e4b8f094-04c8-11eb-3bf7-93bd85db2dc9
# ╠═5897b63a-04c9-11eb-34ab-b908c731d3d5
# ╟─51c9222c-0518-11eb-3827-8fcf8a59b0fd
# ╠═602f56ce-0518-11eb-27d3-7fde99127b36
# ╟─7dbac3ec-00d1-11eb-2c17-3bab40ffaa2e
# ╠═770ab5ac-fff8-11ea-1ed1-87b0ae3aca70
# ╟─b5016eda-0522-11eb-1d8b-37135a22c7bf
# ╠═e9001c9a-0522-11eb-0d63-f718249e03ec

max_stars_repo_name: ysaereve/LABoy

max_issues_repo_name: ysaereve/LABoy

max_forks_repo_name: ysaereve/LABoy


text:
@testset "target" begin
    @test_throws LLVMException Target("invalid")

    host_triple = triple()
    host_t = Target(host_triple)

    host_name = name(host_t)
    description(host_t)

    @test hasjit(host_t)
    @test hastargetmachine(host_t)
    @test hasasmparser(host_t)

    # target iteration
    let ts = targets()
        @test haskey(ts, host_name)
        @test ts[host_name] == host_t

        @test !haskey(ts, "invalid")
        @test get(ts, "invalid", "dummy") == "dummy"
        @test_throws KeyError ts["invalid"]

        @test eltype(ts) == Target

        first(ts)

        for t in ts
            # ...
        end

        @test any(t -> t == host_t, collect(ts))
    end
end

ext: jl


text:
struct Material
    Z::Int
    A::Int
    width::Float64
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
using EchelleTools
using Dates
(fns,zbsrv) = EchelleTools.HARPSN.get_harpsn_data_filenames_and_bsrv(date_start=Date(2015,1,1),date_stop=Date(2015,9,1))
data = EchelleTools.HARPSN.read_data(fns)
zlist = map(d->d[:z_bc],data.metadata_list)
blaze_list = read_blaze_list(data) # Need to inject lambda limits from somewhere into blaze_list

sm = SpectralModel(data, blaze_list=blaze_list)
normalizations =  EchelleTools.Blaze.calc_normalizations(blaze_list,data)
obsm = ObsSeriesModel(num_obs(data), z_bc=zlist, normalization=normalizations)
order = 30
obsid = 1
om = ObservationModel(obsm,obsid)
som = SpectralOrderModel(sm,order)

m11 = eval_model(som, om)[som.blaze.pixel_range]
m1 = eval_model(sm, om, log.(data.lambda[:,:,obsid]), )
m = eval_model(sm, obsm, log.(data.lambda))
#chisq =  eval_model_fit(data, sm, obsm) # WARNING: slow
num_pixels_summed = mapreduce(i->length(SpectralOrderModel(sm,i).blaze.pixel_range),+,1:num_orders(data))

using Plots
pyplot()
plot(data.lambda[:,order,obsid],m2[:,order,obsid]./obsm.normalization[obsid],markersize=1)
plot!(data.lambda[:,order,obsid+1],m2[:,order,obsid+1]./obsm.normalization[obsid+1],markersize=1)
plot!(data.lambda[:,order,obsid+2],m2[:,order,obsid+2]./obsm.normalization[obsid+2],markersize=1)
scatter!(data.lambda[:,order,obsid],data.flux[:,order,obsid]./obsm.normalization[obsid],markersize=1)
scatter!(data.lambda[:,order,obsid+1],data.flux[:,order,obsid+1]./obsm.normalization[obsid+1],markersize=1)
scatter!(data.lambda[:,order,obsid+2],data.flux[:,order,obsid+2]./obsm.normalization[obsid+2],markersize=1)

#=
flat_blaze_list = make_flat_blaze_polynomial_2d(data)
blaze_list = make_blaze_polynomial_2d(data)
sm = SpectralModel(data)
normalizations = calc_normalizations(blaze_list,data)
ave_normalization = mean(normalizations)
blaze_list = map(i->renormalize(blaze_list[i],ave_normalization), 1:length(blaze_list))
blaze_list = read_blaze_list()
normalizations ./= ave_normalization
=#

#=
using Statistics

order = 1
obsid = 1
obs_i = EchelleTools.ObservationModel( star_score=view(obsm.star_score,:,obsid), telluric_score=view(obsm.telluric_score,:,obsid), z_bc=obsm.z_bc[obsid], normalization=obsm.normalization[obsid])
llr = make_model_wavelength_range(EchelleObservation(data,obsid), order)
som = EchelleTools.make_spectral_order_model(data, llr, order)

llr = make_model_wavelength_range(data)
sm = make_spectral_model(data, z_bc=zlist)

evalm = eval_model(llr,sm[1],obs_i)

obsm = make_observations_model(data, sm, z_bc=zlist)
evalm = eval_model(llr,sm[1],1)
=#

#=
obsid = 1
obs_i = EchelleTools.ObservationModel( star_score=view(obsm.star_score,:,obsid), telluric_score=view(obsm.telluric_score,:,obsid), z_bc=obsm.z_bc[obsid], normalization=obsm.normalization[obsid])
=#

#=
using Plots
pyplot()
plot(exp.(llr),evalm)
scatter!(data.lambda[:,1,1],data.flux[:,1,1],markersize=1)
=#

#=
import EchelleTools.AbstractBlazeType
import EchelleTools.AbstractEchelleObservationSet
import EchelleTools.SpectralOrderModel
import EchelleTools.num_obs
import EchelleTools.num_orders
using Interpolations
import EchelleTools.calc_blaze
import EchelleTools.make_stellar_basis
import EchelleTools.default_telluric_template
import EchelleTools.make_flat_blaze

=#

#=
som = EchelleTools.SpectralOrderModel(llr)
obsm = EchelleTools.ObservationsModel(length(zlist))
(blaze1, templ1) = make_blaze_and_template_1d(harpsn_data,order=1,z_bc=zlist)

templ = EchelleTools.make_template_1d(harpsn_data,llr, 1)
starb = EchelleTools.default_stellar_template(length(llr))
telluricb = EchelleTools.default_telluric_template(length(llr))
pixrng = 1:length(llr)
blaz = make_flat_blaze(pixel_range=pixrng)
m = EchelleTools.SpectralOrderModel(llr, star_basis=starb, telluric_basis=telluricb, blaze=blaz )
=#


text:
module JuliaProbo

using Reexport
import RecipesBase: @recipe
@reexport using Plots
import Plots: Plot, plot, plot!, scatter!, annotate!, quiver!, text, @userplot, @series
import Random: rand
import Distributions: Exponential, Normal, Uniform, MvNormal, pdf, mean, cov, det
import LinearAlgebra: Diagonal, I, eigen, diagm
import StatsBase: sample, Weights
import StatsPlots: covellipse!
import StatsFuns: chisqinvcdf

export AbstractObject,
    AbstractWorld,
    AbstractSensor,
    AbstractRobot,
    AbstractAgent,
    AbstractEstimator,
    AbstractMDPAgent
export draw
export PoseUniform, uniform
export Agent, EstimatorAgent, LoggerAgent, decision
export IdealRobot, RealRobot, state_transition
export IdealCamera, RealCamera, PsiCamera, observations, observation_function
export Landmark, EstimatedLandmark, Map, World, push!, getindex
export Particle,
    copy,
    Mcl,
    KalmanFilter,
    KldMcl,
    ResetMcl,
    AMcl,
    motion_update,
    matQ,
    matH,
    matM,
    matA,
    matF
export MapParticle, FastSlam1, FastSlam2
export Goal, Puddle, PuddleWorld, puddle_depth, inside, update_status
export PuddleIgnoreAgent,
    PolicyEvaluator,
    init_value,
    init_policy,
    init_state_transition_probs,
    init_depth,
    action_value,
    policy_evaluation_sweep,
    value_iteration_sweep,
    policy_iteration_sweep,
    DynamicProgramming
export DpPolicyAgent,
    QMDPAgent,
    BeliefDP,
    init_motion_sigma_transition_probs,
    init_obs_sigma_transition_probs,
    init_expected_depths,
    AMDPPolicyAgent

# `include` order does matter(needs to be topologically sorted based on the type definition dependency)
include("types.jl")
include("util.jl")
include("world.jl")
include("agent.jl")
include("robot.jl")
include("sensor.jl")
include("mcl_estimator.jl")
include("kf_estimator.jl")
include("fastslam.jl")
include("mdp_world.jl")
include("mdp.jl")
include("pomdp.jl")

const DynamicProgramming = PolicyEvaluator

end # module

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    PolyhedronProjection

Functor object for using polyhedral projection in the prox step of a quasigradient algorithm. Create by supplying a [`Polyhedron`](@ref) object through `prox ` to `QuasiGradient.Optimizer` or by setting the [`Prox`](@ref) attribute.

...
# Parameters
- `penaltyterm::PenaltyTerm = Quadratic`: Specify penaltyterm variant ([`Quadratic`](@ref), [`InfNorm`](@ref), [`ManhattanNorm`][@ref])
...
"""
struct PolyhedronProjection{T <: AbstractFloat, PT <: AbstractPenaltyTerm} <: AbstractProximal
    penaltyterm::PT
    projection_targets::Vector{MOI.VariableIndex}
    ξ::Vector{Decision{T}}

    function PolyhedronProjection(penaltyterm::AbstractPenaltyTerm, ξ₀::AbstractVector, ::Type{T}) where T <: AbstractFloat
        PT = typeof(penaltyterm)
        projection_targets = Vector{MOI.VariableIndex}(undef, length(ξ₀))
        ξ = map(ξ₀) do val
            KnownDecision(val, T)
        end
        return new{T, PT}(penaltyterm, projection_targets, ξ)
    end
end

function initialize_prox!(quasigradient::AbstractQuasiGradient, polyhedron::PolyhedronProjection)
    # Add projection targets
    ξ = polyhedron.ξ
    decisions = get_decisions(quasigradient.structure.first_stage)::Decisions
    for i in eachindex(ξ)
        name = add_subscript(:ξ, i)
        set = SingleDecisionSet(1, ξ[i], NoSpecifiedConstraint(), false)
        var_index, _ = MOI.add_constrained_variable(quasigradient.master, set)
        set_decision!(decisions[1], var_index, ξ[i])
        MOI.set(quasigradient.master, MOI.VariableName(), var_index, name)
        polyhedron.projection_targets[i] = var_index
    end
    # Initialize penaltyterm
    F = MOI.ScalarAffineFunction{Float64}
    MOI.set(quasigradient.master, MOI.ObjectiveFunction{F}(), zero(F))
    initialize_penaltyterm!(polyhedron.penaltyterm,
                            quasigradient.master,
                            1.0,
                            all_decisions(decisions, 1),
                            polyhedron.projection_targets)
    return nothing
end

function restore_proximal_master!(quasigradient::AbstractQuasiGradient, polyhedron::PolyhedronProjection)
    # Delete penalty-term
    remove_penalty!(polyhedron.penaltyterm, quasigradient.master)
    # Delete projection targets
    for var in polyhedron.projection_targets
        MOI.delete(quasigradient.master, var)
    end
    empty!(polyhedron.projection_targets)
    return nothing
end

function prox!(quasigradient::AbstractQuasiGradient, polyhedron::PolyhedronProjection, x::AbstractVector, ∇f::AbstractVector, γ::AbstractFloat)
    decisions = get_decisions(quasigradient.structure.first_stage)::Decisions
    # Update projection targets
    for i in eachindex(polyhedron.ξ)
        polyhedron.ξ[i].value = x[i] - γ * ∇f[i]
    end
    # Update penaltyterm
    update_penaltyterm!(polyhedron.penaltyterm,
                        quasigradient.master,
                        1.0,
                        all_decisions(decisions, 1),
                        polyhedron.projection_targets)
    # Solve projection problem
    MOI.optimize!(quasigradient.master)
    # Get solution
    x .= MOI.get.(quasigradient.master, MOI.VariablePrimal(), all_decisions(decisions, 1))
    return nothing
end

# API
# ------------------------------------------------------------
"""
    Polyhedron

Factory object for [`PolyhedronProjection`](@ref). Pass to `prox` in `Quasigradient.Optimizer` or set the [`Prox`](@ref) attribute. See ?PolyhedronProjection for parameter descriptions.

"""
mutable struct Polyhedron <: AbstractProx
    penaltyterm::AbstractPenaltyTerm
end
Polyhedron(; penaltyterm = Quadratic()) = Polyhedron(penaltyterm)

function (polyhedron::Polyhedron)(structure::StageDecompositionStructure, x₀::AbstractVector, ::Type{T}) where T <: AbstractFloat
    return PolyhedronProjection(polyhedron.penaltyterm, x₀, T)
end

function str(::Polyhedron)
    return ""
end


text:
#! format: off

using SIIPExamples
using PowerSystems
using PowerSimulations
const PSI = PowerSimulations

using PowerSystemCaseBuilder
using Dates
using DataFrames
using TimeSeries

using HiGHS #solver

sys = build_system(PSITestSystems, "modified_RTS_GMLC_DA_sys")

MultiDay = collect(
    DateTime("2020-01-01T00:00:00"):Hour(1):(DateTime("2020-01-01T00:00:00") + Hour(8783)),
);

remove_time_series!(sys, DeterministicSingleTimeSeries)

for gen in get_components(ThermalGen, sys)
    varcost = get_operation_cost(gen)
    data = TimeArray(MultiDay, repeat([get_cost(get_variable(varcost))], 8784))
    _time_series = SingleTimeSeries("variable_cost", data)
    add_time_series!(sys, gen, _time_series)
    #set_variable_cost!(sys, gen, _time_series)
end

transform_single_time_series!(sys, 24, Dates.Hour(24))

uc_template = template_unit_commitment()

set_device_model!(uc_template, ThermalMultiStart, ThermalMultiStartUnitCommitment)

solver = optimizer_with_attributes(HiGHS.Optimizer, "mip_rel_gap" => 0.5)

problem = DecisionModel(uc_template, sys, horizon = 4, optimizer = solver)
build!(problem, output_dir = mktempdir())

solve!(problem)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
import Base.string
using Base.Test, MuKanren

include("funcTests.jl")
include("macroTests.jl")
include("macroTests2.jl")

FactCheck.exitstatus()


text:
# using Revise

using DistributedFactorGraphs
using IncrementalInference
using Test

@testset "testing compare functions for variables and factors..." begin

fg = initfg()

v1 = addVariable!(fg, :x0, ContinuousScalar)
f1 = addFactor!(fg, [:x0;], Prior(Normal()))

@test compareVariable(v1,v1)
@test compareFactor(f1,f1)

v2 = addVariable!(fg, :x1, ContinuousScalar)
f2 = addFactor!(fg, [:x0;:x1], LinearConditional(Normal(2.0, 0.1)))

fg2 = deepcopy(fg)

@test !compareVariable(v1,v2)

# not testing different factors in this way
# @test !compareFactor(f1,f2)

@test compareAllVariables(fg, fg)

@test compareAllVariables(fg, fg2)

@test compareSimilarVariables(fg, fg)
@test compareSimilarVariables(fg, fg2)

@test compareSimilarFactors(fg, fg)
@test compareSimilarFactors(fg, fg2)

@test compareFactorGraphs(fg, fg)
@test compareFactorGraphs(fg, fg2)

tree, smt, hist = solveTree!(fg)

x1a = getVariable(fg, :x0)
x1b = getVariable(fg2, :x0)

@test !compareVariable(x1a, x1b, skipsamples=false)

@test !compareSimilarVariables(fg, fg2, skipsamples=false)
@test !compareSimilarFactors(fg, fg2, skipsamples=false)

@test compareFactorGraphs(fg, fg)
@test !compareFactorGraphs(fg, fg2, skipsamples=false)

ensureAllInitialized!(fg2)

@test compareSimilarVariables(fg, fg2, skipsamples=true, skip=Symbol[:initialized;:inferdim;:estimateDict])
# fg2 has been solved, so it should fail on the estimate dictionary
@test !compareSimilarVariables(fg, fg2, skipsamples=true, skip=Symbol[:initialized;:inferdim])

tree = wipeBuildNewTree!(fg2)

@test compareSimilarFactors(fg, fg2, skipsamples=true, skipcompute=true)

@test !compareSimilarFactors(fg, fg2, skipsamples=true, skipcompute=false)

@test compareFactorGraphs(fg, fg2, skipsamples=true, skipcompute=true, skip=[:initialized;:inferdim;:estimateDict])

end

@testset "test subgraph functions..." begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal()))

addVariable!(fg, :x1, ContinuousScalar)
addFactor!(fg, [:x0;:x1], LinearConditional(Normal(2.0, 0.1)))

addVariable!(fg, :x2, ContinuousScalar)
addFactor!(fg, [:x1;:x2], LinearConditional(Normal(4.0, 0.1)))

addVariable!(fg, :l1, ContinuousScalar)
addFactor!(fg, [:x1;:l1], LinearConditional(Rayleigh()))

sfg = buildSubgraphFromLabels!(fg, [:x0;:x1])

@test compareFactorGraphs(fg, sfg, skip=[:labelDict;:addHistory;:logpath])

# drawGraph(sfg)

end

#


text:
# Lens implemented using global variables

"Default Lens: do nothing"
deflenscall(::Type{T}, x) where T = nothing

mutable struct Wrapper
  x
end

const glens = Wrapper(deflenscall)

function setglens!(newglens)
  glens.x = newglens
end

setglens!(pair::Pair) = setglens!(lmap(pair))
setglens!(tpl::Tuple) = setglens!(lmap(tpl...))

function lens(::Type{T}, x) where T
  glens.x(T, x)
end

function lenscall(newlens, f, arg)
  setglens!(newlens)
  res = f(arg)
  setglens!(deflenscall)
  res
end

@inline lenscall(pair::Pair, f, arg) = lenscall(lmap(pair), f, arg)
@inline lenscall(tpl::Tuple, f, arg) = lenscall(lmap(tpl...), f, arg)

"Lens call f(args...; kwargs)"
function lenscall(newlens, f, args...; kwargs...)
  setglens!(newlens)
  res = f(args...; kwargs...)
  setglens!(deflenscall)
  res
end

@inline lenscall(pair::Pair, f, arg...; kwargs...) = lenscall(lmap(pair), f, arg...; kwargs...)
@inline lenscall(tpl::Tuple, f, args...; kwargs...) = lenscall(lmap(tpl...), f, args...; kwargs...)

function lmap(l1::Type{T1}, f1) where {T1}
  g(_, x) = nothing
  g(::Type{T1}, x) = f1(x)
  g
end

function lmap(l1::Type{T1}, f1, l2::Type{T2}, f2) where {T1, T2}
  g(_, x) = nothing
  g(::Type{T1}, x) = f1(x)
  g(::Type{T2}, x) = f2(x)
  g
end

function lmap(l1::Type{T1}, f1, l2::Type{T2}, f2, l3::Type{T3}, f3) where {T1, T2, T3}
  g(_, x) = nothing
  g(::Type{T1}, x) = f1(x)
  g(::Type{T2}, x) = f2(x)
  g(::Type{T3}, x) = f3(x)
  g
end

function lmap(l1::Type{T1}, f1, l2::Type{T2}, f2, l3::Type{T3}, f3, l4::Type{T4}, f4) where {T1, T2, T3, T4}
  g(_, x) = nothing
  g(::Type{T1}, x) = f1(x)
  g(::Type{T2}, x) = f2(x)
  g(::Type{T3}, x) = f3(x)
  g(::Type{T4}, x) = f4(x)
  g
end

@inline inner(::Type{T}) where T = T
@inline lmap((l1, f1)::Pair) = lmap(inner(l1), f1)
@inline lmap((l1, f1)::Pair, (l2, f2)::Pair) = lmap(inner(l1), f1, inner(l2), f2)
@inline lmap((l1, f1)::Pair, (l2, f2)::Pair, (l3, f3)::Pair) = lmap(inner(l1), f1, inner(l2), f2, inner(l3), f3)
@inline lmap((l1, f1)::Pair, (l2, f2)::Pair, (l3, f3)::Pair, (l4, l4)::Pair) = lmap(inner(l1), f1, inner(l2), f2, inner(l3), f3, inner(l4), f4)

"""
Lensed eval

```julia
struct L end
function g(x, y)
  lens(L, (x = x, y = y))
  2x + y
end

  @leval L => println ∘ sum ∘ values g(1, 2)
```
"""
macro leval(lmap, expr)
  quote
    setglens!($(esc(lmap)))
    res = $(esc(expr))
    setglens!(deflenscall)
    res
  end
end


text:
using Test
using LinearAlgebra
using Random
using Strided
using Strided: StridedView, UnsafeStridedView

Random.seed!(1234)

println("Base.Threads.nthreads() =  $(Base.Threads.nthreads())")

println("Running tests single-threaded:")
Strided.disable_threads()
include("othertests.jl")
include("blasmultests.jl")

println("Running tests multi-threaded:")
Strided.enable_threads()
include("othertests.jl")
include("blasmultests.jl")

Strided.enable_threaded_mul()
include("blasmultests.jl")


text:
using NonconvexCore
using NonconvexCore: VectorOfFunctions, FunctionWrapper, CountingFunction,
                     getdim, getfunction, getfunctions, Objective, IneqConstraint,
                     EqConstraint, SDConstraint
using Test

@testset "CountingFunction" begin
    f = CountingFunction(FunctionWrapper(log, 1))
    for i in 1:3
        x = rand()
        @test log(x) == f(x)
        @test f.counter[] == i
    end
    @test getdim(f) == 1
    @test getfunction(f) isa FunctionWrapper
    @test getfunction(getfunction(f)) === log
    @test length(f) == getdim(f)
end

@testset "VectorOfFunctions" begin
    f1 = FunctionWrapper(log, 1)
    f2 = FunctionWrapper(exp, 1)
    fs = [f1, f2]
    f = VectorOfFunctions([f1, f2])
    for _ in 1:3
        x = rand()
        @test [log(x), exp(x)] == f(x)
    end
    @test getdim(f) == 2
    @test getfunctions(f) == fs
    @test getfunction(f, 1) == fs[1]
    @test getfunction(f, 2) == fs[2]
    newf = VectorOfFunctions(Any[f, FunctionWrapper(exp, 1)])
    @test getdim(newf) == 3
end

@testset "Objective" begin
    f = Objective(log)
    @test getdim(f) == 1
    for _ in 1:3
        x = rand()
        @test f(x) == log(x)
    end
    @test getfunction(f) === log
end

@testset "IneqConstraint" begin
    f = IneqConstraint(log, 1.0)
    for _ in 1:3
        x = rand()
        @test f(x) == log(x) - 1
    end
    @test getdim(f) == 1
    @test getfunction(f) === log
end

@testset "EqConstraint" begin
    f = EqConstraint(log, 1.0)
    for _ in 1:3
        x = rand()
        @test f(x) == log(x) - 1
    end
    @test getdim(f) == 1
    @test getfunction(f) === log
end

@testset "SDConstraint" begin
    f = SDConstraint(exp, 1)
    for _ in 1:3
        x = rand()
        @test f(x) == exp(x)
    end
    @test getdim(f) == 1
    @test getfunction(f) === exp

    f = SDConstraint(exp, 2)
    for _ in 1:3
        x = rand(2, 2)
        @test f(x) == exp(x)
    end
    @test getdim(f) == 2
    @test getfunction(f) === exp
end


text:
using RoME
using RoMEPlotting
using IncrementalInference
using DistributedFactorGraphs

# Find the filezzzzz.
fg_dir = "/home/tonio/repos/bayes-garden/data/manhattan/batch-liljon-1500"
fg_name = "fg-after-solve.tar.gz"
fg_file = joinpath(fg_dir, fg_name)

# Load the factor graph.
fg = initfg()
loadDFG(fg_file, Main, fg)

# Draw some pretty plots.
pl1 = drawPoses(fg, spscale=0.6, lbls=false)
pl1 |> PDF(joinpath(fg_dir, "poses.pdf"), 20cm, 10cm)

# Useful functions to deal with sorting symbols.

# For crying out loud, find a way to lexicographically sort the symbols.
function natural(x, y)
    k(x) = [occursin(r"\d+", s) ? parse(Int, s) : s
            for s in split(replace(x, r"\d+" => s->" $s "))]
    A = k(x); B= k(y)
    for (a, b) in zip(A, B)
        if !isequal(a, b)
            return typeof(a) <: typeof(b) ? isless(a, b) :
                   isa(a,Int) ? true : false
        end
    end
    return length(A) < length(B)
end

function getSortedSymbols(fg::AbstractDFG, lexsortfxn)
    # First transform the factor graph symbols into strings to sort them
    # lexicographically, in case numbers are not padded.
    stringsyms = []
    for sym in sort(ls(fg))
        push!(stringsyms, "$(sym)")
    end
    sort!(stringsyms, lt=lexsortfxn)

    # Now transform sorted string array into sorted symbols array.
    sortedsyms = []
    for sym in stringsyms
        push!(sortedsyms, Symbol(sym))
    end
    return sortedsyms
end

# Output ground truth to file.
sortedsyms = getSortedSymbols(fg, natural)
fid = open(joinpath(fg_dir, "ground-truth.txt"), "w")
for i in 1:length(sortedsyms)
    @show i
    @show sortedsyms[i]
    @show getVariablePPE(fg, sortedsyms[i]).suggested
    x_hat = getVariablePPE(fg, sortedsyms[i]).suggested
    quat = normalize([0; 0; sin(x_hat[3]/2); cos(x_hat[3]/2)])
    println(fid, "$(i-1) $(x_hat[1]) $(x_hat[2]) 0.0 $(quat[1]) $(quat[2]) $(quat[3]) $(quat[4])")
end
close(fid)

# Encapsulate all of this into a function:
# Where `lexsort` is a function to sort lexicographically (e.g., `natural`).
function exportTrajEstimatesToFile(fg_path::String, lexsort)
    # Load the factor graph.
    fg = initfg()
    loadDFG(fg_path, Main, fg)

    # Draw some pretty plots.
    pl1 = drawPoses(fg, spscale=0.6, lbls=false)
    pl1 |> PDF(joinpath(dirname(fg_path), "trajectory-estimates-$(basename(dirname(fg_path))).pdf"), 20cm, 10cm)

    # Output ground truth to file.
    sortedsyms = getSortedSymbols(fg, lexsort)
    fid = open(joinpath(dirname(fg_path), "trajectory-estimates-$(basename(dirname(fg_path))).txt"), "w")
    for i in 1:length(sortedsyms)
        x_hat = getVariablePPE(fg, sortedsyms[i]).suggested
        quat = normalize([0; 0; sin(x_hat[3]/2); cos(x_hat[3]/2)])
        println(fid, "$(i-1) $(x_hat[1]) $(x_hat[2]) 0.0 $(quat[1]) $(quat[2]) $(quat[3]) $(quat[4])")
    end
    close(fid)
end

# Now run this for all of our runs:
b100_new_liljon_1411 = "/home/tonio/repos/bayes-garden/data/manhattan/b100-new-liljon-1411/fg-after-solve1411.tar.gz"
exportTrajEstimatesToFile(b100_new_liljon_1411, natural)

b10_old_augustus_1571 = "/home/tonio/repos/bayes-garden/data/manhattan/b10-old-augustus-1571/fg-after-solve1571.tar.gz"
exportTrajEstimatesToFile(b10_old_augustus_1571, natural)

b30_old_liljon_1451 = "/home/tonio/repos/bayes-garden/data/manhattan/b30-old-liljon-1451/fg-after-solve1451.tar.gz"
exportTrajEstimatesToFile(b30_old_liljon_1451, natural)

incremental_liljon_1101 = "/home/tonio/repos/bayes-garden/data/manhattan/incremental-liljon-1101/fg-after-solve1101.tar.gz"
exportTrajEstimatesToFile(incremental_liljon_1101, natural)

batch_liljon_1500 = "/home/tonio/repos/bayes-garden/data/manhattan/batch-liljon-1500/fg-after-solve.tar.gz"
exportTrajEstimatesToFile(batch_liljon_1500, natural)

batch_liljon_1800 = "/home/tonio/repos/bayes-garden/data/manhattan/batch-liljon-1800/fg-after-solve.tar.gz"
exportTrajEstimatesToFile(batch_liljon_1800, natural)


text:
using DifferentialEquations
using Plots
using DiffEqBase
using PowerDynamics
using NetworkDynamics
using NetworkDynamics: collect_ve_info

const iipfunc = true # is in-place function

"""
The almost the same as the original NetworkDynamics-function. Only ODE is switched to DAE
"""
function dae_network_dynamics(vertices, edges, graph; x_prototype = zeros(1)) #where {T <: ODEVertex, U <: ODEEdge}
    v_dims, e_dims, symbols_v, symbols_e, mmv_array, mme_array = collect_ve_info(vertices, edges, graph)

    x_array = similar(x_prototype, sum(v_dims) + sum(e_dims))

    v_array = view(x_array, 1:sum(v_dims))
    e_array = view(x_array, sum(v_dims)+1:sum(v_dims)+sum(e_dims))

    symbols = vcat(symbols_v, symbols_e)

    graph_stucture = GraphStruct(graph, v_dims, e_dims, symbols_v, symbols_e)

    graph_data = GraphData(v_array, e_array, graph_stucture)

    nd! = nd_ODE_ODE(vertices, edges, graph, graph_stucture, graph_data)

    mass_matrix = construct_mass_matrix(mmv_array, graph_stucture)

    #DAEFunction{iipfunc,true}(nd!; mass_matrix = mass_matrix, syms=symbols)

    # mass_matrix is not a supported kwarg here... ill look into that later
    DAEFunction{iipfunc,true}(nd!;syms=symbols)
end

"""
Almost the same as the PowerDynamics function but calls dae_network_dynamics
"""
function dae_rhs(pg::PowerGrid)
    dae_network_dynamics(map(construct_vertex, pg.nodes), map(construct_edge, pg.lines), pg.graph)
end

"""
Almost the same as the PowerDynamics-function, Switched to a DEA Problem!
"""
function dae_solve(pg::PowerGrid, x0, timespan)
    problem = DAEProblem(dae_rhs(pg),zeros(systemsize(pg)),x0.vec,timespan)
    solution = solve(problem,DABDF2(),initializealg = ShampineCollocationInit())
    PowerGridSolution(solution, pg)
end

"""
Calls dae_solve
"""
function dae_simulate(p::Perturbation, powergrid, x0; timespan)
    dae_solve(powergrid, p(x0), timespan);
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
function func(a, b)
    a+b
end


text:
# Digital IO API

#------------------------------------ Control
function digitalIOReset(hdwf::Hdwf)
  succ = ccall((:FDwfDigitalIOReset,libdwf),Cint,(Hdwf,),
    hdwf)
  succ == 0 && error("Error calling FDwfDigitalIOReset.")
  return nothing
end

function digitalIOConfigure(hdwf::Hdwf)
  succ = ccall((:FDwfDigitalIOConfigure,libdwf),Cint,(Hdwf,),
    hdwf)
  succ == 0 && error("Error calling FDwfDigitalIOConfigure.")
  return nothing
end

function digitalIOStatus(hdwf::Hdwf)
  succ = ccall((:FDwfDigitalIOStatus,libdwf),Cint,(Hdwf,),
    hdwf)
  succ == 0 && error("Error calling FDwfDigitalIOStatus.")
  return nothing
end

#------------------------------------ Configuration
function digitalIOOutputEnableInfo(hdwf::Hdwf)
  pfsOutputEnableMask = Ref{Cuint}()
  succ = ccall((:FDwfDigitalIOOutputEnableInfo,libdwf),Cint,(Hdwf,Ptr{Cuint}),
    hdwf,pfsOutputEnableMask)
  succ == 0 && error("Error calling FDwfDigitalIOOutputEnableInfo.")
  return pfsOutputEnableMask[]
end

function digitalIOOutputEnableSet(hdwf::Hdwf,fsOutputEnable::UInt32)
  succ = ccall((:FDwfDigitalIOOutputEnableSet,libdwf),Cint,(Hdwf,Cuint),
    hdwf,fsOutputEnable)
  succ == 0 && error("Error calling FDwfDigitalIOOutputEnableSet.")
  return nothing
end

function digitalIOOutputEnableGet(hdwf::Hdwf)
  pfsOutputEnable = Ref{Cuint}()
  succ = ccall((:FDwfDigitalIOOutputEnableGet,libdwf),Cint,(Hdwf,Ptr{Cuint}),
    hdwf,pfsOutputEnable)
  succ == 0 && error("Error calling FDwfDigitalIOOutputEnableGet.")
  return pfsOutputEnable[]
end

function digitalIOOutputInfo(hdwf::Hdwf)
  pfsOutputMask = Ref{Cuint}()
  succ = ccall((:FDwfDigitalIOOutputInfo,libdwf),Cint,(Hdwf,Ptr{Cuint}),
    hdwf,pfsOutputMask)
  succ == 0 && error("Error calling FDwfDigitalIOOutputInfo.")
  return pfsOutputMask[]
end

function digitalIOOutputSet(hdwf::Hdwf,fsOutput::UInt32)
  succ = ccall((:FDwfDigitalIOOutputSet,libdwf),Cint,(Hdwf,Cuint),
    hdwf,fsOutput)
  succ == 0 && error("Error calling FDwfDigitalIOOutputSet.")
  return nothing
end

function digitalIOOutputGet(hdwf::Hdwf)
  pfsOutput = Ref{Cuint}()
  succ = ccall((:FDwfDigitalIOOutputGet,libdwf),Cint,(Hdwf,Ptr{Cuint}),
    hdwf,pfsOutput)
  succ == 0 && error("Error calling FDwfDigitalIOOutputGet.")
  return pfsOutput[]
end

function digitalIOInputInfo(hdwf::Hdwf)
  pfsInputMask = Ref{Cuint}()
  succ = ccall((:FDwfDigitalIOInputInfo,libdwf),Cint,(Hdwf,Ptr{Cuint}),
    hdwf,pfsInputMask)
  succ == 0 && error("Error calling FDwfDigitalIOInputInfo.")
  return pfsInputMask[]
end

function digitalIOInputStatus(hdwf::Hdwf)
  pfsInput = Ref{Cuint}()
  succ = ccall((:FDwfDigitalIOInputStatus,libdwf),Cint,(Hdwf,Ptr{Cuint}),
    hdwf,pfsInput)
  succ == 0 && error("Error calling FDwfDigitalIOInputStatus.")
  return pfsInput[]
end

function digitalIOInputStatus!(pfsInput::Ref{Cuint},hdwf::Hdwf)
  succ = ccall((:FDwfDigitalIOInputStatus,libdwf),Cint,(Hdwf,Ptr{Cuint}),
    hdwf,pfsInput)
  succ == 0 && error("Error calling FDwfDigitalIOInputStatus.")
  return pfsInput[]
end


text:
"""
Basic discrete multivariate parameter state type

# Constructors

## BasicDiscMuvParameterState{NI<:Integer, NR<:Real}(value::Vector{NI}, <optional arguments>)

Construct a basic discrete multivariate parameter state with some ``value``.

###  Optional arguments:

* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``::Type{NR}=Float64``: the element type of target-related fields.
* ``diagnosticvalues::Vector=Array(Any, length(diagnostickeys))``: the diagnostic values of the state.

## BasicDiscMuvParameterState{NI<:Integer, NR<:Real}(size::Integer, <optional arguments>)

Construct a basic discrete multivariate parameter state with with a ``value`` of specified ``size``.

###  Optional arguments:

* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``::Type{NI}=Integer``: the element type of the state value.
* ``::Type{NR}=Float64``: the element type of target-related fields.
* ``diagnosticvalues::Vector=Array(Any, length(diagnostickeys))``: the diagnostic values of the state.

# Examples

```julia
julia> state = BasicDiscMuvParameterState(Int64[0, 1], [:accept], Float64, [false])
Klara.BasicDiscMuvParameterState{Int64,Float64}([0,1],NaN,NaN,NaN,Bool[false],2,[:accept])

julia> state.value
2-element Array{Int64,1}:
 0
 1

julia> diagnostics(state)
Dict{Symbol,Bool} with 1 entry:
  :accept => false
```
"""
type BasicDiscMuvParameterState{NI<:Integer, NR<:Real} <: ParameterState{Discrete, Multivariate}
  "Vector value of basic discrete multivariate parameter state"
  value::Vector{NI}
  "Value of log-likelihood at the state's value"
  loglikelihood::NR
  "Value of log-prior at the state's value"
  logprior::NR
  "Value of log-target at the state's value"
  logtarget::NR
  "Diagnostic values associated with the sampling of the state"
  diagnosticvalues::Vector
  "Integer representing the length of vector value of basic discrete multivariate parameter state"
  size::Integer
  "Diagnostic keys associated with the sampling of the state"
  diagnostickeys::Vector{Symbol}
end

function BasicDiscMuvParameterState{NI<:Integer, NR<:Real}(
  value::Vector{NI},
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{NR}=Float64,
  diagnosticvalues::Vector=Array(Any, length(diagnostickeys))
)
  v = convert(NR, NaN)
  BasicDiscMuvParameterState{NI, NR}(value, v, v, v, diagnosticvalues, length(value), diagnostickeys)
end

BasicDiscMuvParameterState{NI<:Integer, NR<:Real}(
  size::Integer,
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{NI}=Int64,
  ::Type{NR}=Float64,
  diagnosticvalues::Vector=Array(Any, length(diagnostickeys))
) =
  BasicDiscMuvParameterState(Array(NI, size), diagnostickeys, NR, diagnosticvalues)

value_support{NI<:Integer, NR<:Real}(::Type{BasicDiscMuvParameterState{NI, NR}}) = Discrete
value_support(::BasicDiscMuvParameterState) = Discrete

variate_form{NI<:Integer, NR<:Real}(::Type{BasicDiscMuvParameterState{NI, NR}}) = Univariate
variate_form(::BasicDiscMuvParameterState) = Univariate

Base.eltype{NI<:Integer, NR<:Real}(::Type{BasicDiscMuvParameterState{NI, NR}}) = (NI, NR)
Base.eltype{NI<:Integer, NR<:Real}(::BasicDiscMuvParameterState{NI, NR}) = (NI, NR)

generate_empty(state::BasicDiscMuvParameterState) =
  BasicDiscMuvParameterState(state.size, state.diagnostickeys, eltype(state)...)


text:
####
#### Atom
####

"""
    Atom{T}

Represents a single atom and its spatial position of type `CoordT`.

# Examples
```jldoctest
julia> Atom(:Li, (0.0, 0.0, 1.4))
Atom{Float64}(:Li, (0.0, 0.0, 1.4))
```
"""
struct Atom{CoordT}
    species::Symbol
    coords::NTuple{3, CoordT}

    function Atom(species::Symbol, coords::NTuple{3, T}) where T
        haskey(PeriodicTable.elements, species) || throw(ErrorException("No element \"$species\" exists."))
        return new{T}(species, coords)
    end
end

####
#### Geometry
####

"""
    Geometry{CoordT}

Represents the geometry of a molecule, a `Vector` of `Atom{CoordT}`s, each of
which represents an atomic species and its position. `CoordT` is the type of the
coordinates.
"""
struct Geometry{CoordT}
    atoms::Vector{Atom{CoordT}}
end

Base.copy(g::Geometry) = Geometry(copy(g.atoms))
Base.:(==)(g1::Geometry, g2::Geometry) = g1.atoms == g2.atoms

"""
    Geometry()

Create an empty `Geometry` object with `Float64` coordinates.

# Examples
```jldoctest
julia> Geometry()
Geometry{Float64}(Atom{Float64}[])
```
"""
Geometry() = Geometry(Atom{Float64}[])

"""
    Geometry(atoms::Atom...)

Initialize a `Geometry` object with `atoms`.

# Examples
```jldoctest
julia> Geometry(Atom(:H, (0., 0., 0.)), Atom(:H, (0., 0., 0.7414)))
Geometry{Float64}(Atom{Float64}[Atom{Float64}(:H, (0.0, 0.0, 0.0)), Atom{Float64}(:H, (0.0, 0.0, 0.7414))])
```
"""
Geometry(atoms::Atom...) = Geometry([atoms...])

for f in (:push!, :length, :getindex, :iterate)
    @eval Base.$f(g::Geometry, args...) = $f(g.atoms, args...)
end

####
#### Geometry
####

"""
    struct MolecularSpec{CoordT}

Contains data specifying an electronic structure problem. The results of electronic
structure calculations, for example, integrals are *not* stored here. They are stored
in the [`MolecularData`](@ref) type.

# Keywords:
- `geometry::Geometry`
-  `multiplicity::Int = 1`
-  `charge::Int = 0`
-  `basis::String = "sto-3g"`

# Additional properties
-  `spin`

# Examples
```jldoctest
julia> geometry = Geometry(Atom(:H, (0., 0., 0.)), Atom(:H, (0., 0., 0.7414)));

julia> MolecularSpec(geometry=geometry)
MolecularSpec{Float64}(Geometry{Float64}(Atom{Float64}[Atom{Float64}(:H, (0.0, 0.0, 0.0)), Atom{Float64}(:H, (0.0, 0.0, 0.7414))]), 1, 0, "sto-3g")
```
"""
Base.@kwdef struct MolecularSpec{CoordT}
    geometry::Geometry{CoordT}  #  Should we use just `Vector{Atom{CoordT}}` ?
    multiplicity::Int = 1
    charge::Int = 0
    basis::String = "sto-3g"
end

Base.copy(spec::MolecularSpec) = MolecularSpec(copy(spec.geometry), spec.multiplicity, spec.charge, spec.basis)

function Base.:(==)(spec1::MolecularSpec, spec2::MolecularSpec)
    return spec1.geometry == spec2.geometry &&
        spec1.multiplicity == spec2.multiplicity &&
        spec1.charge == spec2.charge &&
        spec1.basis == spec2.basis
end

function Base.getproperty(mol::MolecularSpec, sym::Symbol)
    if sym == :spin
        return mol.multiplicity - 1  # from OpenFermion
    else
        return getfield(mol, sym)
    end
end


text:
# This module provides routines for writing PGM/PBM images in Julia
# Owain Kenway

module imageio

export writepgm, writepbm

# Function for writing a PGM image.
# d is the pixel value array.
# white is the max vale which corresponds to white.
# filename is the file to write to.
function writepgm(d::Array{Int64,2}, white::Int64, filename::AbstractString)

  x = size(d, 1)
  y = size(d, 2)

  f = open(filename, "w")

# Write header
  write(f, "P2\n")
  write(f, "# Written by pnmmodules (https://github.com/owainkenwayucl/pnmmodules).\n")
  write(f, "$(x) $(y)\n")
  write(f, "$(white)\n")

  for j = 1:y
    for i = 1:x
      write(f, "$(d[i,j])\n")
    end
  end

  close(f)

end

# Function for writing a PBM image.
# d is the pixel value array.
# threshold is the value at which 0 becomes 1
# filename is the file to write to.
function writepbm(d::Array{Int64,2}, threshold::Float64, filename::AbstractString)

  x = size(d, 1)
  y = size(d, 2)

  f = open(filename, "w")

# Write header
  write(f, "P1\n")
  write(f, "# Written by pnmmodules (https://github.com/owainkenwayucl/pnmmodules).\n")
  write(f, "$(x) $(y)\n")

  for j = 1:y
    for i = 1:x
      if Float64(d[i,j]) >= threshold
        write(f, "1\n")
      else
        write(f, "0\n")
      end
    end
  end

  close(f)

end

end

max_stars_repo_name: owainkenwayucl/pnmmodules

max_issues_repo_name: owainkenwayucl/pnmmodules

max_forks_repo_name: owainkenwayucl/pnmmodules


text:
MathKernel.area(a::BB)   where {BB <: BoundingBox} = 8 * dot(a.e, SVector{3,Float64}(a.e[2], a.e[3], a.e[1]))
MathKernel.volume(a::BB) where {BB <: BoundingBox} = 8 * prod(a.e)


text:
using ReTest # useful when setup.jl included separately

module Trace
export trace
using ReTest: @test

const RUN = []

function trace(x)
    push!(RUN, x)
    @test true
end
end # module Trace

# NOTE: all keywords have the same defaults as `retest`, except `marks`
# Also, not sure why ReTest.def has to be qualified, when we just import def from ReTest
# in this module and call it unqualified, then `def` is said to not exist when running
# the test suite, from within a call to @chapter
function check(x...; runtests=false, output::Union{Nothing,String}=nothing,
               verbose=ReTest.def(:verbose), stats=ReTest.def(:stats), dry=ReTest.def(:dry),
               strict::Bool=ReTest.def(:strict), recursive=ReTest.def(:recursive),
               static=ReTest.def(:static), id=ReTest.def(:id), load=ReTest.def(:load),
               marks::Bool=false, tag=ReTest.def(:tag), clear::Bool=ReTest.def(:clear))
    @assert !(runtests & (output !== nothing)) "unimplemented"
    args = x[1:end-1]
    expected = x[end]
    if expected isa AbstractString
        expected = split(x[end])
    end

    empty!(Trace.RUN)
    if runtests
        getfield(args[1], :runtests)(args[2:end]...; verbose=verbose, stats=stats, dry=dry,
                                     strict=strict, recursive=recursive, static=static,
                                     id=id, load=load, marks=marks, tag=tag, clear=clear)
    elseif output === nothing
        retest(args...; verbose=verbose, stats=stats, dry=dry, strict=strict,
               recursive=recursive, static=static, id=id, load=load, marks=marks,
               tag=tag, clear=clear)
    else
        mktemp() do path, io
            redirect_stdout(io) do
                retest(args...; verbose=verbose, stats=stats, dry=dry, strict=strict,
                       recursive=recursive, static=static, id=id, load=load, marks=marks,
                       tag=tag, clear=clear)
            end
            seekstart(io)
            printed = join(map(rstrip, split(readchomp(io), '\n')), '\n')
            output = join(map(rstrip, split(chomp(output), '\n')), '\n')
            @test printed == output
        end
    end
    @test Trace.RUN == expected
end

# @chapter title [wrap_in_testset::Bool=true] body
macro chapter(title, wrap_in_testset, body=nothing)
    title = string(title)
    if body === nothing
        body = wrap_in_testset
        wrap_in_testset = true
    else
        wrap_in_testset::Bool
    end

    if isempty(ARGS) || any(pat -> occursin(Regex(pat, "i"), title), ARGS)
        printstyled("\n\n", rpad("## $title #", 78, '#'), "\n\n", bold=true, color=:cyan)
        if wrap_in_testset
            quote
                Test.@testset $("$title") begin
                    $(esc(body))
                end
            end
        else
            esc(body)
        end
    end
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment, JSON, GZip, DataStructures

parse_case14() = JSON.parse(GZip.gzopen("../instances/test/case14.json.gz"),
                            dicttype=()->DefaultOrderedDict(nothing))

@testset "Validation" begin
    @testset "fix!" begin

        @testset "Cost curve should be convex" begin
            json = parse_case14()
            json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150, 200]
            json["Generators"]["g1"]["Production cost curve (\$)"] = [10, 25, 30]
            instance = UnitCommitment.from_json(json, fix=false)
            @test UnitCommitment.fix!(instance) == 4
        end

        @testset "Startup limit must be greater than Pmin" begin
            json = parse_case14()
            json["Generators"]["g1"]["Production cost curve (MW)"] = [100, 150]
            json["Generators"]["g1"]["Production cost curve (\$)"] = [100, 150]
            json["Generators"]["g1"]["Startup limit (MW)"] = 80
            instance = UnitCommitment.from_json(json, fix=false)
            @test UnitCommitment.fix!(instance) == 1
        end

        @testset "Startup costs and delays must be increasing" begin
            json = parse_case14()
            json["Generators"]["g1"]["Startup costs (\$)"] = [300, 200, 100]
            json["Generators"]["g1"]["Startup delays (h)"] = [8, 4, 2]
            instance = UnitCommitment.from_json(json, fix=false)
            @test UnitCommitment.fix!(instance) == 4
        end

    end
end


text:
using GlyphSheets
using Base.Test

library = create_library()
@test_throws ErrorException create_typeface(library, "NOT_FOUND.ttf")
close_library(library)


text:
module Test1

using Pmodules
@Parent

@P import Test1.C: value
@P import Test1: B, A

a_value = A.value
b_value = B.value
c_value = value

end # module


text:
is_scale(::Function) = true
is_scale(::AbstractArray) = true
is_scale(::Any) = false

function compute_attribute(scale::Function, el, rankdict)
    scale(el)
end
function compute_attribute(scale::AbstractArray, el, rankdict)
    scale[mod1(rankdict[el], length(scale))]
end

function rankdict(col::AbstractVector)
    s = collect(uniquesorted(col))
    return Dict(zip(s, 1:length(s)))
end

rankdicts(s) = map(rankdict, Tables.columntable(s))

max_stars_repo_name: piever/StatMakie

max_issues_repo_name: piever/StatMakie

max_forks_repo_name: piever/StatMakie


text:
# Find a Profile-most Probable k-mer in a String

struct ProfileColumn
    A::Float64
    C::Float64
    G::Float64
    T::Float64
end

convert_to_floatlist(line) = map(x -> parse(Float64, x), split(line))

function read_profile_matrix(lines, k)
    matrix = Array{Float64, 2}(undef, 4, k)
    for (ix, line) in enumerate(lines)
        matrix[ix, :] = convert_to_floatlist(line)
    end
    [ProfileColumn(column...) for column in eachcol(matrix)]
end

calc_kmer_probability(kmer, profile_matrix) =
    foldl((p, (nucleotide, column)) -> p * getfield(column, Symbol(nucleotide)), zip(kmer, profile_matrix); init=1.0)

function profile_most_probable_kmer(text, profile_matrix, k)
    max_probability = 0.0
    most_probable_kmer = text[1:k]
    for ix in 1:(length(text)-k+1)
        kmer = text[ix:ix+k-1]
        p = calc_kmer_probability(kmer, profile_matrix)
        if p > max_probability
            max_probability = p
            most_probable_kmer = kmer
        end
    end
    most_probable_kmer
end

function main()
    text = readline()
    k = parse(Int64, readline())
    lines = readlines()
    profile_matrix = read_profile_matrix(lines, k)
    result = profile_most_probable_kmer(text, profile_matrix, k)
    println(result)
end

main()


text:
using ModularForms

include("poly_to_power_series-test.jl")
include("delta-test.jl")
include("eisenstein-test.jl")
include("generators-test.jl")
include("victormiller-test.jl")
include("hecke-test.jl")
include("prime_range-test.jl")
include("eta_quotient-test.jl")
include("dimensions-test.jl")

function test_all()
   test_poly_to_power_series()
   test_delta()
   test_eis_series()
   test_generators()
   test_vm_basis()
   test_hecke()
   test_prime_range()
   test_eta_quotient()
   test_dimensions()
end


text:
"""
 LogWeight
 represents a function on [-1,1] weighted by log((1+x)^α*(1-x)^β)
"""
immutable LogWeight{S,DD} <: WeightSpace{S,RealBasis,DD,1}
    α::Float64
    β::Float64
    space::S
end
LogWeight(α,β,space)=LogWeight{typeof(space),typeof(domain(space))}(α,β,space)

spacescompatible(A::LogWeight,B::LogWeight)=A.α==B.α && A.β == B.β && spacescompatible(A.space,B.space)
canonicalspace(A::LogWeight)=A

logweight(α,β,x)=log((1+x).^α.*(1-x).^β)
weight(sp::LogWeight,x)=logweight(sp.α,sp.β,tocanonical(sp,x))

setdomain(sp::LogWeight,d::Domain)=LogWeight(sp.α,sp.β,setdomain(sp.space,d))

function coefficients(f::Vector,sp1::LogWeight,sp2::LogWeight)
    α,β=sp1.α,sp1.β
    c,d=sp2.α,sp2.β

    if isapprox(c,α) && isapprox(d,β)
        coefficients(f,sp1.space,sp2.space)
    else
        (Conversion(sp1,sp2)*f)
    end
end

for (OPrule,OP) in ((:maxspace_rule,:maxspace),(:union_rule,:union))
    @eval begin
        function $OPrule(A::LogWeight,B::LogWeight)
            if isapprox(A.α,B.α) && isapprox(A.β,B.β)
                LogWeight(A.α,A.β,$OP(A.space,B.space))
            else
                NoSpace()
            end
        end
        # there are no other comatible spaces
        # this is mostly to overcome typing bug
        # in spacestes
        $OPrule(::LogWeight,::JacobiWeight)=NoSpace()
        $OPrule(::JacobiWeight,::LogWeight)=NoSpace()
    end
end

# Same as JacobiWeight

# avoid redundency
function Multiplication{SS,LWS,DD<:IntervalDomain,T}(f::Fun{JacobiWeight{SS,DD},T},S::LogWeight{LWS,DD})
    M=Multiplication(Fun(f.coefficients,space(f).space),S)
    rsp=JacobiWeight(space(f).α,space(f).β,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end

function Multiplication(f::Fun,S::LogWeight)
    M=Multiplication(f,S.space)
    rsp=LogWeight(S.α,S.β,rangespace(M))
    MultiplicationWrapper(f,SpaceOperator(M,S,rsp))
end


text:
using CVChannel
using Test

"""
This script shows the Siddhu and dephrasure channels tensored
admit multiplicative communication value.
"""

print("\n Here we verify the multiplicativity of the Siddhu and dephrasure")
print(" channel when ran in parallel. Note as a special")
print(" case this considers the erasure channel.")
@testset "Verifying Multiplicativity of Siddhu Channel with Dephrasure Channel" begin
    sid_scan_range =[0:0.1:0.5;]
    dephrasure_scan_range = [0:0.1:1;]
    sid_ctr, p_ctr, q_ctr = 1,1,1
    max_val = 0
    is_mult = true
    for s_id in sid_scan_range
        println("---Now scanning for s =", s_id, ".-----")
        for q_id in dephrasure_scan_range
            println("Now scanning for q =", q_id,".")
            target_val = 2*(2-q_id)
            for p_id in dephrasure_scan_range
                dephrasurepq(X) = dephrasureChannel(X,p_id,q_id)
                sidchan(X) = siddhuChannel(X,s_id)

                dephr_chan = Choi(dephrasurepq,2,3)
                sid_chan= Choi(sidchan,3,3)

                par_cv, = pptCV(parChoi(sid_chan, dephr_chan), :dual)

                non_mult = par_cv - target_val
                !isapprox(non_mult,0,atol=4e-6) ? is_mult = false : nothing
                max_val < non_mult ? max_val = non_mult : nothing

                p_ctr += 1
            end
            q_ctr += 1
            p_ctr = 1
        end
        sid_ctr += 1
        q_ctr = 1
    end
    @test is_mult
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
import Core: typeof

struct Foo
  a::Int
  b::Int
end

function typeof(x::Foo)
  println("here")
end

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
using Random;

function sample(labeled_indices, unlabeled_indices, labeled_bs, unlabeled_bs)
    labeled_indices = shuffle(labeled_indices)
    i = 1
    batches = []
    while i<=length(labeled_indices)
        lb = i+labeled_bs-1
        if lb > length(labeled_indices)
            lb = length(labeled_indices)
        end
        push!(batches, (labeled_indices[i:lb], shuffle(unlabeled_indices)[1:unlabeled_bs]))
        i = i+labeled_bs
    end
    batches
end

max_stars_repo_name: skoc/julia-project

max_issues_repo_name: skoc/julia-project

max_forks_repo_name: skoc/julia-project


text:
#####
##### Shortcut for truncating uniform distributions.
#####

truncated(d::Uniform, l::T, u::T) where {T <: Real} = Uniform(promote(max(l, d.a), min(u, d.b))...)

truncated(d::Uniform, l::Integer, u::Integer) = truncated(d, float(l), float(u))