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

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

# ╔═╡ f8e2f7e1-d6a7-43ef-879c-db33c5f66924
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ c44c6c6a-3b8c-46e4-93ad-365d03082226
using SymPy, LinearAlgebra

# ╔═╡ c226e28d-3a8a-42f6-8a51-2ac418023a27
md"""
# Singular Value Decomposition - Definitions and Facts

__Prerequisites__

The reader should be familiar with basic linear algebra concepts and notebooks related to eigenvalue decomposition.

__Competences__

The reader should be able to undestand and check the facts about singular value decomposition.

__Selected references__

There are many excellent books on the subject. Here we list a few:

* J.W. Demmel, _Applied Numerical Linear Algebra_, SIAM, Philadelphia, 1997.
* G. H. Golub and C. F. Van Loan, _Matrix Computations_, 4th ed., The John Hopkins University Press, Baltimore, MD, 2013.
* N. Higham, _Accuracy and Stability of Numerical Algorithms_, SIAM, Philadelphia, 2nd ed., 2002.
* L. Hogben, ed., _Handbook of Linear Algebra_, CRC Press, Boca Raton, 2014.
* B. N. Parlett, _The Symmetric Eigenvalue Problem_, Prentice-Hall, Englewood Cliffs, NJ, 1980, also SIAM, Philadelphia, 1998.
* G. W. Stewart, _Matrix Algorithms, Vol. II: Eigensystems_, SIAM, Philadelphia, 2001.
* L. N. Trefethen and D. Bau, III, _Numerical Linear Algebra_, SIAM, Philadelphia, 1997.

"""

# ╔═╡ ddbd255b-7ddd-4fba-ac4d-722e88801e59
md"""
# Singular value problems

For more details and the proofs of the Facts below, see [R. C. Li, Matrix Perturbation Theory, pp. 21.6-21.8](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and [R. Mathias, Singular Values and Singular Value Inequalities,pp. 24.1-24.17](https://www.routledge.com/Handbook-of-Linear-Algebra/Hogben/p/book/9781138199897) and the references therein.

## Definitions

Let $A\in\mathbb{C}^{m\times n}$ and let $q=\min\{m,n\}$.

The __singular value decomposition__ (SVD) of $A$ is

$$A=U\Sigma V^*,$$

where $U\in\mathbb{C}^{m\times m}$ and $V\in\mathbb{C}^{n\times n}$ are unitary, and $\Sigma=\mathop{\mathrm{diag}}(\sigma_1,\sigma_2,\ldots)\in\mathbb{R}^{m\times n}$ with all $\sigma_j\geq 0$.

Here $\sigma_j$ is the __singular value__, $u_j\equiv U_{:,j}$ is the corresponding __left singular vector__, and $v_j\equiv V_{:,j}$ is the corresponding __right singular vector__.

The __set of singular values__ is $sv(A)=\{\sigma_1,\sigma_2,\ldots,\sigma_{q}\}$.

We assume that singular values are ordered, $\sigma_1\geq\sigma_2\geq\cdots\sigma_q\geq 0$.

The __Jordan-Wielandt__ matrix is the Hermitian matrix

$$
J=\begin{bmatrix}0 & A \\ A^* & 0 \end{bmatrix}
\in \mathbb{C}^{(m+n) \times (m+n)}.$$
"""

# ╔═╡ caa44e3e-3af5-4836-8ae3-dd3b9d68d222
md"""
## Facts

There are many facts related to the singular value problem for general matrices. We state some basic ones:

1. If $A\in\mathbb{R}^{m\times n}$, then $U$ and $V$ are real.

2. Singular values are unique (uniquely determined by the matrix).

3.  $\sigma_j(A^T)=\sigma_j(A^*)=\sigma_j(\bar A)=\sigma_j(A)$ for $j=1,2,\ldots,q$.

4.  $A v_j=\sigma_j u_{j}$ and $A^* u_{j}=\sigma_j v_{j}$ for $j=1,2,\ldots,q$.

5.  $A=\sigma_1 u_{1} v_{1}^* + \sigma_2 u_{2} v_{2}^* +\cdots + \sigma_q u_{q} v_{q}^*$.

6. __Unitary invariance.__ For any unitary $U\in\mathbb{C}^{m\times m}$ and $V\in\mathbb{C}^{n\times n}$, $sv(A)=sv(UAV)$.

7. There exist unitary matrices $U\in\mathbb{C}^{m\times m}$ and $V\in\mathbb{C}^{n\times n}$ such that $A=UBV$ if and only if  $sv(A)=sv(B)$.

8. SVD of $A$ is related to eigenvalue decompositions of Hermitian matrices $A^*A=V\Sigma^T\Sigma V^*$ and $AA^*=U\Sigma\Sigma^TU^*$. Thus, $\sigma_j^2(A)=\lambda_j(A^*A)=\lambda_j(AA^*)$ for $j=1,2,\ldots,q$.

9. The eigenvalues of Jordan-Wielandt matrix are $\pm \sigma_1(A), \pm\sigma_2(A), \cdots,\pm\sigma_q(A)$ together with $|m-n|$ zeros. The eigenvectors are obtained from an SVD of $A$. This relationship is used to deduce singular value results from the results for eigenvalues of Hermitian matrices.

10.  $\mathop{\mathrm{trace}}(|A|_{spr})=\sum_{i=1}^q \sigma_i$, where $|A|_{spr}=(A^*A)^{1/2}$.

11. If $A$ is square, then $|\det(A)|=\prod_{i=1}^n \sigma_i$.

12. If $A$ is square, then $A$ is singular $\Leftrightarrow$ $\sigma_j(A)=0$ for some $j$.

13. __Min-max Theorem.__ It holds:

$$
\begin{aligned}
\sigma_k& =\max_{\dim(W)=k}\min_{x\in W, \ \|x\|_2=1} \|Ax\|_2\\
& =\min_{\dim(W)=n-k+1}\max_{x\in W, \ \|x\|_2=1} \|Ax\|_2.
\end{aligned}$$

14.  $\|A\|_2=\sigma_1(A)$.

15. For $B\in\mathbb{C}^{m\times n}$,

$$
|\mathop{\mathrm{trace}}(AB^*)|\leq \sum_{j=1}^q \sigma_j(A)\sigma_j(B).$$

16. __Interlace Theorems.__ Let $B$ denote $A$ with the one of its rows _or_ columns deleted. Then

$$
\sigma_{j+1}(A)\leq \sigma_j(B)\leq \sigma_j(A),\quad j=1,\ldots,q-1.$$

Let $B$ denote $A$ with a row _and_ a column deleted. Then

$$
\sigma_{j+2}(A)\leq \sigma_j(B)\leq \sigma_j(A),\quad j=1,\ldots,q-2.$$

17. __Weyl Inequalities.__ For $B\in\mathbb{C}^{m\times n}$, it holds:

$$
\begin{aligned}
\sigma_{j+k-1}(A+B)& \leq \sigma_j(A)+\sigma_k(B), \quad  j+k\leq q+1,\\
\sum_{j=1}^k \sigma_j(A+B)& \leq \sum_{j=1}^k \sigma_j(A) + \sum_{j=1}^k \sigma_j(B),
\quad k=1,\ldots,q.
\end{aligned}$$

"""

# ╔═╡ 4b764569-a32d-4170-a11d-2072d971846a
md"""
## Examples

### Symbolic computation
"""

# ╔═╡ 772342da-a960-4aa7-b30b-7122b049eff9
A=[  3   2   1
 -5  -1  -4
  5   0   2]

# ╔═╡ 314c41e0-03e7-4f8e-851c-a4374e7f5092
x=symbols("x")

# ╔═╡ 4ab83338-0428-4942-979e-c1021b81c309
B=A'*A

# ╔═╡ bf02e40e-94f4-4c69-ab74-3496cbcd6cdf
# Characteristic polynomial p_B(λ)
p(x)=simplify(det(B-x*Matrix{Int}(I,3,3)))

# ╔═╡ 54264dd9-98cf-4c8f-b6f0-c6ee2ba11cdb
λ=solve(p(x),x)

# ╔═╡ 6cc74bec-e335-471d-87a0-9b9bad391a1d
map(Float64,λ)

# ╔═╡ 2d1ea310-24a8-448c-9a3d-0a4420bbc90c
begin
	V=Array{Any}(undef,3,3)
	for j=1:3
	    V[:,j]=nullspace(B-map(Float64,λ[j])*I)
	end
	V
end

# ╔═╡ 92b68122-5e9a-443a-aaa3-ad301cb1878a
begin
	U=Array{Any}(undef,3,3)
	for j=1:3
	    U[:,j]=nullspace(A*A'-map(Float64,λ[j])*I)
	end
	U
end

# ╔═╡ edaaa469-6a73-43e6-a8b3-b365bf265e04
σ=sqrt.(λ)

# ╔═╡ c1ac5d19-4971-461f-b73d-10aabb4eb56f
# Diagonalizing
V'*A'*A*V

# ╔═╡ 72352834-0bf5-41b7-86e2-0b2cb733d7f4
# Residual -what is wrong?
A-U*Diagonal(float(σ))*V'

# ╔═╡ ace071e5-764a-4c04-8916-7287ed53a01e
A-U*Diagonal(float(σ))*(V*Diagonal(sign.(U'*A*V)))'

# ╔═╡ 8d911867-9b90-4bff-9d0c-7358dfcc04ec
#?svd

# ╔═╡ 1b1b7ddd-f50d-4d65-8cc8-398955dcce34
# Using the structure SVD
S=svd(A)

# ╔═╡ b6663c47-1293-43ba-8c13-20858054c975
S.Vt

# ╔═╡ 74e584f4-a233-4f18-b04c-dfd1710515e1
# Residual
A-S.U*Diagonal(S.S)*S.Vt

# ╔═╡ 6ba22886-faa3-4f30-af96-366d9c001c94
md"""
Observe the signs of the columns!
"""

# ╔═╡ 80d341c5-43ae-4568-bcf3-ced1947662ce
S.U

# ╔═╡ 1245fe85-25a8-44cb-bbd7-3bc477859d04
U

# ╔═╡ abe3e648-65e6-423c-83de-673e96013670
md"
### Random complex matrix
"

# ╔═╡ 64535706-eaba-4155-b593-12cba7366127
begin
	import Random
	Random.seed!(421)
	m=5
	n=3
	q=min(m,n)
	A₁=rand(ComplexF64,m,n)
end

# ╔═╡ bd88592a-781f-46e7-a175-342b78906c45
S₁=svd(A₁,full=true)

# ╔═╡ f54a7c12-d175-45cf-b996-d13024910135
norm(A₁-S₁.U[:,1:q]*Diagonal(S₁.S)*S₁.Vt), norm(S₁.U'*S₁.U-I), norm(S₁.Vt*S₁.V-I)

# ╔═╡ b39c5987-d320-499a-816c-31608eb1d856
@bind k Slider(1:q,show_value=true)

# ╔═╡ e7aa4857-76c0-42da-964a-cecd428ec70e
# Fact 4
norm(A₁*S₁.V[:,k]-S₁.S[k]*S₁.U[:,k]), norm(A₁'*S₁.U[:,k]-S₁.S[k]*S₁.V[:,k])

# ╔═╡ da2e1129-474e-4a6a-9d66-7224f1f25b34
λ₁,V₁=eigen(A₁'*A₁)

# ╔═╡ b740803d-7bb2-4ce2-8338-4a53e333f9f9
sqrt.(λ₁)

# ╔═╡ 43316c8b-149b-4a52-90e9-ccb2c2ca4d67
S₁.V

# ╔═╡ 8a06f9fc-62d0-4394-8f5c-b861a63da947
abs.(S₁.Vt*V₁)

# ╔═╡ 879285b0-8fc9-43a2-b06d-30fd95556d61
md"""
__Explain non-uniqueness of $U$ and $V$!__
"""

# ╔═╡ ff0ea61d-7944-4316-a224-72461a808cdc
# Jordan-Wielandt matrix
J=[zero(A*A') A; A' zero(A'*A)]

# ╔═╡ 5d941c94-b186-49ba-9f03-1613a638ecb4
λJ,UJ=eigen(J)

# ╔═╡ 9909c421-607d-4883-9b35-705290ce8bf3
λJ

# ╔═╡ dda804cf-71c1-4735-bfc7-a97453698c8f
 float(σ)

# ╔═╡ f4b2ce9d-4783-420d-bc87-2ce27671c7e4
md"
### Random real matrix
"

# ╔═╡ 87de7452-77af-466b-920c-1a562ee05ea2
begin
	m₂=8
	n₂=5
	q₂=min(m₂,n₂)
	A₂=rand(-9:9,m₂,n₂)
end

# ╔═╡ 27f76447-12cf-4096-a1ab-ad997367b336
S₂=svd(A₂)

# ╔═╡ 63af09b0-1b32-478e-868f-95ddf2e951e8
# Fact 10
tr(sqrt(A₂'*A₂)), sum(S₂.S)

# ╔═╡ c30340fd-cc4e-4cbe-b68a-fb0b01f383fa
begin
	# Fact 11
	B₁=rand(n₂,n₂)
	det(B₁), prod(svdvals(B₁))
end

# ╔═╡ ecae308e-7444-4883-84bf-8106bb7106e9
# Fact 14
norm(A₂), opnorm(A₂), S₂.S[1]

# ╔═╡ eb15023f-f126-4177-bc2c-057543b8fcc1
begin
	# Fact 15
	B₂=rand(m₂,n₂)
	abs(tr(A₂*B₂')), svdvals(A₂)⋅svdvals(B₂)
end

# ╔═╡ 07f30928-36aa-4d59-bb53-fa8df06fdad0
@bind j Slider(1:q₂,show_value=true)

# ╔═╡ bf28d15d-c243-4229-8d41-0945c63e3baa
# Interlace Theorems (repeat several times)
S₂.S, svdvals(A₂[[1:j-1;j+1:m₂],:]), svdvals(A₂[:,[1:j-1;j+1:n₂]])

# ╔═╡ 69470862-8321-4953-b324-3a5b4e93d9c8
# Weyl Inequalities
# B=rand(m,n)
[svdvals(A₂+B₂) S₂.S svdvals(B₂)]

# ╔═╡ 35fead5b-4425-40b3-a64d-9cf7ebab18d2
@bind l Slider(1:q₂,show_value=true)

# ╔═╡ 61ecca0f-1e5e-4dec-90d3-924961b5f9c8
sum(svdvals(A₂+B₂)[1:l]),sum(svdvals(A₂)[1:l].+sum(svdvals(B₂)[1:l]))

# ╔═╡ 2d804bbd-c6ce-49ce-9197-d7850ccb7db5
md"""
# Matrix approximation

Let $A=U\Sigma V^*$, let $\tilde \Sigma$ be equal to $\Sigma$ except that $\tilde \Sigma_{jj}=0$ for $j>k$, and let $\tilde A=U\tilde \Sigma V^*$. Then $\mathop{\mathrm{rank}}(\tilde A)\leq k$ and

$$
\begin{aligned}
\min\{\|A-B\|_2: \mathop{\mathrm{rank}}(B)\leq k\} & =\|A-\tilde A\|_2=\sigma_{k+1}(A)\\
\min\{\|A-B\|_F: \mathop{\mathrm{rank}}(B)\leq k\} & =\|A-\tilde A\|_F=
\bigg(\sum_{j=k+1}^{q}\sigma_{j}^2(A)\bigg)^{1/2}.
\end{aligned}$$

This is the __Eckart-Young-Mirsky Theorem__.
"""

# ╔═╡ d7221da7-78fc-4d86-a693-95bab6d12d5a
A₂

# ╔═╡ 4d9640da-d85c-41ec-afcf-fa9548c25804
S₂.S

# ╔═╡ 997d0c9c-22a4-4a60-9730-0668eda26af4
@bind k₀ Slider(1:q₂,show_value=true)

# ╔═╡ 76bc9a58-c163-4df2-a315-5f9474e746aa
B₀=S₂.U*Diagonal([S₂.S[1:k₀];zeros(q₂-k₀)])*S₂.Vt

# ╔═╡ 36fb7249-ecc5-432a-9842-1c9cbf17c30a
opnorm(A₂-B₀), S₂.S[k₀+1]

# ╔═╡ 42ec2898-6676-4254-b553-215da8ec7490
norm(A₂-B₀),norm(S₂.S[k₀+1:q₂])

# ╔═╡ Cell order:
# ╟─f8e2f7e1-d6a7-43ef-879c-db33c5f66924
# ╟─c226e28d-3a8a-42f6-8a51-2ac418023a27
# ╟─ddbd255b-7ddd-4fba-ac4d-722e88801e59
# ╟─caa44e3e-3af5-4836-8ae3-dd3b9d68d222
# ╟─4b764569-a32d-4170-a11d-2072d971846a
# ╠═c44c6c6a-3b8c-46e4-93ad-365d03082226
# ╠═772342da-a960-4aa7-b30b-7122b049eff9
# ╠═314c41e0-03e7-4f8e-851c-a4374e7f5092
# ╠═4ab83338-0428-4942-979e-c1021b81c309
# ╠═bf02e40e-94f4-4c69-ab74-3496cbcd6cdf
# ╠═54264dd9-98cf-4c8f-b6f0-c6ee2ba11cdb
# ╠═6cc74bec-e335-471d-87a0-9b9bad391a1d
# ╠═2d1ea310-24a8-448c-9a3d-0a4420bbc90c
# ╠═92b68122-5e9a-443a-aaa3-ad301cb1878a
# ╠═edaaa469-6a73-43e6-a8b3-b365bf265e04
# ╠═c1ac5d19-4971-461f-b73d-10aabb4eb56f
# ╠═72352834-0bf5-41b7-86e2-0b2cb733d7f4
# ╠═ace071e5-764a-4c04-8916-7287ed53a01e
# ╠═8d911867-9b90-4bff-9d0c-7358dfcc04ec
# ╠═1b1b7ddd-f50d-4d65-8cc8-398955dcce34
# ╠═b6663c47-1293-43ba-8c13-20858054c975
# ╠═74e584f4-a233-4f18-b04c-dfd1710515e1
# ╟─6ba22886-faa3-4f30-af96-366d9c001c94
# ╠═80d341c5-43ae-4568-bcf3-ced1947662ce
# ╠═1245fe85-25a8-44cb-bbd7-3bc477859d04
# ╟─abe3e648-65e6-423c-83de-673e96013670
# ╠═64535706-eaba-4155-b593-12cba7366127
# ╠═bd88592a-781f-46e7-a175-342b78906c45
# ╠═f54a7c12-d175-45cf-b996-d13024910135
# ╠═b39c5987-d320-499a-816c-31608eb1d856
# ╠═e7aa4857-76c0-42da-964a-cecd428ec70e
# ╠═da2e1129-474e-4a6a-9d66-7224f1f25b34
# ╠═b740803d-7bb2-4ce2-8338-4a53e333f9f9
# ╠═43316c8b-149b-4a52-90e9-ccb2c2ca4d67
# ╠═8a06f9fc-62d0-4394-8f5c-b861a63da947
# ╟─879285b0-8fc9-43a2-b06d-30fd95556d61
# ╠═ff0ea61d-7944-4316-a224-72461a808cdc
# ╠═5d941c94-b186-49ba-9f03-1613a638ecb4
# ╠═9909c421-607d-4883-9b35-705290ce8bf3
# ╠═dda804cf-71c1-4735-bfc7-a97453698c8f
# ╟─f4b2ce9d-4783-420d-bc87-2ce27671c7e4
# ╠═87de7452-77af-466b-920c-1a562ee05ea2
# ╠═27f76447-12cf-4096-a1ab-ad997367b336
# ╠═63af09b0-1b32-478e-868f-95ddf2e951e8
# ╠═c30340fd-cc4e-4cbe-b68a-fb0b01f383fa
# ╠═ecae308e-7444-4883-84bf-8106bb7106e9
# ╠═eb15023f-f126-4177-bc2c-057543b8fcc1
# ╠═07f30928-36aa-4d59-bb53-fa8df06fdad0
# ╠═bf28d15d-c243-4229-8d41-0945c63e3baa
# ╠═69470862-8321-4953-b324-3a5b4e93d9c8
# ╠═35fead5b-4425-40b3-a64d-9cf7ebab18d2
# ╠═61ecca0f-1e5e-4dec-90d3-924961b5f9c8
# ╟─2d804bbd-c6ce-49ce-9197-d7850ccb7db5
# ╠═d7221da7-78fc-4d86-a693-95bab6d12d5a
# ╠═4d9640da-d85c-41ec-afcf-fa9548c25804
# ╠═997d0c9c-22a4-4a60-9730-0668eda26af4
# ╠═76bc9a58-c163-4df2-a315-5f9474e746aa
# ╠═36fb7249-ecc5-432a-9842-1c9cbf17c30a
# ╠═42ec2898-6676-4254-b553-215da8ec7490

max_stars_repo_path: src/Module B - Eigenvalue and Singular Value Decompositions/L5a Singular Value Decomposition - Definitions and Facts .jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
include("Yfind.jl")
include("is_it_exit2.jl")
include("is_it.jl")
@everywhere function where_did_it_come_from_parallel(K)
    Q_end=1.5
    Q_start=-Q_end
    P_end=3
    P_start=-P_end

    n_iter_P=8001
    n_iter_Q=8001

    delta_P=(P_end-P_start)/(n_iter_P)
    delta_Q=(Q_end-Q_start)/(n_iter_Q)

    ArrP= P_start: delta_P: P_end
    ArrQ= Q_start: delta_Q: Q_end

    Q=ArrQ[Int(ceil(K/n_iter_Q))]
    P=ArrP[mod(K,n_iter_P)+1]

    max_hit=3
    Energy=.25
    t_end=1e8
    tol_dist=1e-7
    H=(2*Energy)^2
    #
    # Q_0=0
    # P_0=0
    #
    # Q_1=0
    # P_1=0
    #
    # Q_2=0
    # P_2=0
    #
    # Q_3=0
    # P_3=0
    #

    Y=Yfind(Q,P,H)
    # YB=YfindB(Q,P,H)

    if ~isempty(Y)
        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y[1] #Y
        u0[5]=0 #Y
        # counter variable
        #are the initial conditions right?

        prob_test= ODEProblem(Eq_of_M,u0,(0., 10.))
        sol_test=solve(prob_test, RK4(),maxiters=1e10 )
        S=sign(sol_test.u[2][1])
        if true#check if upcrossing
            prob_forward = ODEProblem(Eq_of_M,u0,(0., t_end))
            sol_forward=solve(prob_forward, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb, save_start=true,save_end=true,save_everystep=true)
            B=hcat(sol_forward.u[end-5],sol_forward.u[end])'
            if is_it_exit2(B,H,tol_dist) && is_it_exit(sol_forward.u[end],H,tol_dist) && sol_forward.u[end][5]==0
                prob_backward = ODEProblem(Eq_of_M,u0,(0.,-t_end))
                sol_backward=solve(prob_backward, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_backward,save_start=true,save_end=true,save_everystep=false)

                uf=zeros(5)
                uf[1]=sol_backward.u[end][1] #X
                uf[2]=sol_backward.u[end][2] #P
                uf[3]=sol_backward.u[end][3] #Q
                uf[4]=sol_backward.u[end][4] #Y
                dH=abs(H_test(uf)-H)
                # uf[1]=0
                # prob_test2= ODEProblem(Eq_of_M,u0,(0., -10.))
                # sol_test2=solve(prob_test2, RK4(),maxiters=1e10 )
                # S2=sign(sol_test.u[2][1])
                if  dH<1e-8 #&& S2<0
                    # prob_check= ODEProblem(Eq_of_M,uf,(0., t_end))
                    # sol_check=solve(prob_check, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_forward,save_start=true,save_end=true,save_everystep=false)
                    #return sol_backward.u[end][3],sol_backward.u[end][2],sol_backward.u[end][1],sol_backward.u[end][4],sol_backward.u[end][5],sol_check.u[end][3],sol_check.u[end][2],sol_check.u[end][1],sol_check.u[end][4],Y,sol_check.u[end][5]
                    # return Q,P,sol_backward.u[end][1],Y,sol_backward.u[end][5],sol_check.u[end][3],sol_check.u[end][2],sol_check.u[end][1],sol_check.u[end][4],Yi,sol_check.u[end][5]
                    return sol_backward.u[end][3],sol_backward.u[end][2],sol_backward.u[end][1],sol_backward.u[end][4],sol_backward.u[end][5],Y[1]
                end
            end
        end
    end
end

ext: jl


text:
abstract NLregModF{T<:FP}

## default methods for all NLregModF objects
Base.size(m::NLregModF) = size(tgrad(m))
Base.size(m::NLregModF,args...) = size(tgrad(m),args...)
StatsBase.model_response(m::NLregModF) = m.y
StatsBase.nobs(m::NLregModF) = length(model_response(m))
StatsBase.residuals(m::NLregModF) = m.resid
covariatemat(m::NLregModF) = m.x
expctd(m::NLregModF) = m.mu
mufunc(m::NLregModF) = m.f
npars(m::NLregModF) = size(tgrad(m),1)
tgrad(m::NLregModF) = m.tgrad # transposed gradient matrix
unscaledvcov(m::NLregModF) = inv(cholfact(tgrad(m) * tgrad(m)'))

## update the expected response and residuals, returning the sum of squared residuals
function updtmu!(m::NLregModF, pars::Vector)
    length(pars) == npars(m) || throw(DimensionMismatch(""))
    f  = mufunc(m)
    mu = expctd(mu)
    r  = residuals(m)
    tg = tgrad(m)
    x  = covariatemat(m)
    y  = model_response(m)
    for i in 1:length(y)
        mu[i] = f(pars,view(x,:,i),view(tgrad,:,i)) # pass subarrays by reference, not copies
        r[i] = y[i] - mu[i]
        rss += abs2(r[i])
    end
    rss
end

function updtmu!(m::NLregModF, pars::Matrix, inds::Vector)
    (n == nobs(m)) == length(inds) && size(pars,1) == npars(m) || throw(DimensionMismatch(""))
    f  = mufunc(m)
    mu = expctd(mu)
    r  = residuals(m)
    tg = tgrad(m)
    x  = covariatemat(m)
    y  = model_response(m)
    ii = 0; rss = zero(eltype(mu))
    for i in 1:n
        if ii != inds[i]
            ii = inds[i]
            pp = view(pars,:,ii)
        end
        mu[i] = f(pp,view(x,:,i),view(tgrad,:,i))
        r[i]  = y[i] - mu[i]
        rss  += abs2(r[i])
    end
    rss
end

abstract PLregModF{T<:FP} <: NLregModF{T}

## default methods for all PLregModF objects
Base.size(m::PLregModF) = size(mmjac(m))
Base.size(m::PLregModF,args...) = size(mmjac(m),args...)
StatsBase.model_response(m::PLregModF) = m.y
mmfunc(m::PLregModF) = m.mmf # model-matrix update function
mmjac(m::PLregModF) = m.MMD # model-matrix derivative (Jacobian)

## update the (transposed) model matrix for the linear parameters and the gradient MMD
function updtMM!(m::PLregModF,nlpars)
    mmd = mmjac(m)
    mmf = mmfunc(m)
    tg  = tgrad(m)
    x   = covariatemat(m)
    y   = model_response(m)
    for i in 1:size(x,2)
        mmf(nlpars,view(x,:,i),view(tg,:,i),view(mmd,:,:,i))
    end
    tg[1:size(m,2),:]
end

## update mu and resid from full parameter vector returning rss
function updtmu!(m::PLregModF,pars::Vector)
    length(pars) == npars(m) || throw(DimensionMismatch(""))
    mmd = mmjac(m)
    mmf = mmfunc(m)
    mu  = expctd(m)
    r   = residuals(m)
    tg  = tgrad(m)
    x   = covariatemat(m)
    y   = model_response(m)
    nnl,nl,n = size(m.MMD); lind = 1:nl; nlind = nl + (1:nnl)
    nlpars = view(pars,nlind); lpars = view(pars,lind)
    rss = zero(eltype(mu))
    for i in 1:n
        mmdi = view(mmd,:,:,i)
        mmf(nlpars,view(x,:,i),view(tg,:,i),mmdi)
        tg[nlind,i] = mmdi * lpars
        r[i] = y[i] - (mu[i] = dot(view(tg,lind,i),lpars))
        rss += abs2(r[i])
    end
    rss
end

## update mu and resid from parameter matrix and indicators returning rss
function updtmu!(m::PLregModF,pars::Matrix,inds::Vector)
    nnl,nl,n = size(m); lind = 1:nl; nlind = nl + (1:nnl)
    length(inds) == n && size(pars,1) == npars(m) || throw(DimensionMismatch(""))
    mmd = mmjac(m)
    mmf = mmfunc(m)
    mu  = expctd(m)
    r   = residuals(m)
    tg  = tgrad(m)
    x   = covariatemat(m)
    y   = model_response(m)
    nnl,nl,n = size(m); lind = 1:nl; nlind = nl + (1:nnl)
    ii  = 0
    rss = zero(eltype(mu))
    for i in 1:n
        if ii != inds[i]
            ii = inds[i]
            nlp = view(pars,nlind,ii)
            lp = view(pars,lind,ii)
        end
        mmdi = view(mmd,:,:,i)
        mmf(nlp,view(x,:,i),view(tg,:,i),mmdi)
        tg[nlind,i] = mmdi * lp
        r[i] = y[i] - (mu[i] = dot(view(tg,lind,i),lp))
        rss += abs2(r[i])
    end
    rss
end


text:
using GadgetUnits, GadgetIO, Test, Unitful, UnitfulAstro, Downloads, Cosmology

Downloads.download("http://www.usm.uni-muenchen.de/~lboess/GadgetIO/snap_144.0", "./snap_144.0")

@testset "GadgetUnits" begin

    @testset "Unit Structs" begin

        @testset "GadgetPhysicalUnits" begin

            @testset "Float64" begin
                GU = GadgetPhysicalUnits(Float64)

                @test GU.t_s ≈ 3.085678e16u"s"
                @test GU.E_cgs ≈ 1.989e53u"erg"
            end

            @testset "Float32" begin
                GU = GadgetPhysicalUnits(Float32)
                @test GU.t_s        ≈ 3.0856778f16u"s"
            end
        end

        @testset "GadgetPhysical" begin

            @testset "Float64" begin
                GU = GadgetPhysical(Float64)

                @test GU.t_s   ≈ 3.085678e16
                @test GU.E_cgs ≈ 1.989e53
            end

            @testset "Float32" begin
                GU = GadgetPhysical(Float32)

                @test GU.t_s        ≈ 3.085678e16
                @test GU.x_physical ≈ 1.0
            end
        end

    end # Conversion Structs

    @testset "Cosmology" begin

        h = read_header("snap_144.0")
        h.z = 0.5

        @testset "Setup" begin
            c = cosmology(h)
            @test c.h     ≈ h.h0
            @test c.Ω_m   ≈ h.omega_0
            #Cosmology.FlatLCDM{Float64}(0.704, 0.7279156596296679, 0.272, 8.434037033212373e-5)
        end

        @testset "Scale" begin
            c = cosmology(h)

            # tested against script by Klaus Dolag
            @test arcmin_to_kpc(1.0, h)                  ≈ 367.7593389798557
            @test arcmin_to_kpc(c, 1.0, h.z)             ≈ 367.7593389798557
            @test arcmin_to_kpc(1.0u"arcminute", h)      ≈ 367.7593389798557u"kpc"
            @test arcmin_to_kpc(c, 1.0u"arcminute", h.z) ≈ 367.7593389798557u"kpc"
        end

        @testset "Radiation" begin
            c = cosmology(h)

            # tested against script by Klaus Dolag
            @test mJy_to_W(1.0, h)            ≈ 9.681690084092813e23
            @test mJy_to_W(c, 1.0, h.z)       ≈ 9.681690084092813e23
            @test mJy_to_W(1.0u"mJy", h)      ≈ 9.681690084092813e23u"W/Hz"
            @test mJy_to_W(c, 1.0u"mJy", h.z) ≈ 9.681690084092813e23u"W/Hz"
        end

        @testset "Time" begin
            @test age(h)        ≈ 8.69582554048774u"Gyr"
            @test age(h, false) ≈ 8.69582554048774

            @test redshift(8.69582554048774u"Gyr", h) ≈ 0.5
            @test redshift(8.69582554048774, h)       ≈ 0.5

            c = cosmology(h)
            @test redshift(8.69582554048774u"Gyr", c) ≈ 0.5
            @test redshift(8.69582554048774, c)       ≈ 0.5
        end
    end

end

rm("snap_144.0")


text:
"""
    PeriodicKernel(r::AbstractVector)
    PeriodicKernel(dims::Int)
    PeriodicKernel(T::DataType, dims::Int)

Periodic Kernel as described in http://www.inference.org.uk/mackay/gpB.pdf eq. 47.
```
    κ(x,y) = exp( - 0.5 sum_i(sin (π(x_i - y_i))/r_i))
```
"""
struct PeriodicKernel{T} <: SimpleKernel
    r::Vector{T}
    function PeriodicKernel(; r::AbstractVector{T} = ones(Float64, 1)) where {T<:Real}
        @assert all(r .> 0)
        new{T}(r)
    end
end

PeriodicKernel(dims::Int) = PeriodicKernel(Float64, dims)

PeriodicKernel(T::DataType, dims::Int = 1) = PeriodicKernel(r = ones(T, dims))

metric(κ::PeriodicKernel) = Sinus(κ.r)

kappa(κ::PeriodicKernel, d::Real) = exp(- 0.5d)

Base.show(io::IO, κ::PeriodicKernel) = print(io, "Periodic Kernel, length(r) = ", length(κ.r), ")")


text:
using MultiObjectiveProblems
using Test

@testset "MultiObjectiveProblems.jl" begin
    T1 = MHT1();
    @test num_objectives(T1) == num_vars(T1) == 2
    ps = get_pareto_set( T1 );
    ps_p1 = get_points( ps, 10 );
    ps_p2 = get_points( ps, 10; method = :random)
    @test length(ps_p1) == length(ps_p2)
    @test all( [sum(p) ≈ 10 for p ∈ ps_p1] )
    @test all( [sum(p) ≈ 10 for p ∈ ps_p2] )
    @test get_starting_point(T1,1) == [-34.24; 47.06]
end


text:
export VectorModel, fit, vectorize, TfidfModel, TfModel, IdfModel, FreqModel, prune, prune_select_top

"""
    abstract type Model

An abstract type that represents a weighting model
"""
abstract type Model end

abstract type TfidfModel end
abstract type TfModel end
abstract type IdfModel end
abstract type FreqModel end

struct IdFreq
    id::Int32
    freq::Int32
end

const Vocabulary = Dict{Symbol, IdFreq}

"""
    mutable struct VectorModel

Models a text through a vector space
"""
mutable struct VectorModel <: Model
    config::TextConfig
    tokens::Vocabulary
    id2token::Dict{Int,Symbol}
    # id2token::Dict{Int,Symbol}
    maxfreq::Int
    m::Int  # vocsize
    n::Int  # collection size
end

"""
    fit(::Type{VectorModel}, config::TextConfig, corpus::AbstractVector)

Trains a vector model using the text preprocessing configuration `config` and the input corpus.
"""
function fit(::Type{VectorModel}, config::TextConfig, corpus::AbstractVector)
    bow = BOW()
    n = 0
    maxfreq = 0

    for data in corpus
        n += 1
        _, _maxfreq = compute_bow(tokenize(config, data), bow)
        maxfreq = max(maxfreq, _maxfreq)
    end

    tokens = Vocabulary()
    id2token = Dict{Int,Symbol}()
    i = 0
    for (t, freq) in bow
        i += 1
        id2token[i] = t
        tokens[t] = IdFreq(i, freq)
    end

    VectorModel(config, tokens, id2token, Int(maxfreq), length(tokens), n)
end

"""
    prune(model::VectorModel, minfreq, rank)

Cuts the vocabulary by frequency using lower and higher filter;
All tokens with frequency below `minfreq` are ignored; top `rank` tokens are also removed.
"""
function prune(model::VectorModel, minfreq::Integer, rank::Integer)
    W = [(token, idfreq) for (token, idfreq) in model.tokens if idfreq.freq >= minfreq]
    sort!(W, by=x->x[2])
    tokens = Vocabulary()
    for i in 1:length(W)-rank+1
        w = W[i]
        tokens[w[1]] = w[2]
    end

    id2token = Dict(idfreq.id => t for (t, idfreq) in tokens)
    VectorModel(model.config, tokens, id2token, model.maxfreq, model.m, model.n)
end

"""
    prune_select_top(model::VectorModel, k::Integer, kind::Type{T}=IdfModel)
    prune_select_top(model::VectorModel, ratio::AbstractFloat, kind::Type{T}=IdfModel)

Creates a new model with the best `k` tokens from `model` based on the `kind` scheme; kind must be either `IdfModel` or `FreqModel`.
`ratio` is a floating point between 0 and 1 indicating the ratio of the vocabulary to be kept
"""
function prune_select_top(model::VectorModel, k::Integer, kind::Type{T}=IdfModel) where T <: Union{IdfModel,FreqModel}
    tokens = Vocabulary()
    maxfreq = 0
    if kind == IdfModel
        X = [(t, idfreq, _weight(kind, 0, 0, model.n, idfreq.freq)) for (t, idfreq) in model.tokens]
        sort!(X, by=x->x[end], rev=true)
        for i in 1:k
            t, idfreq, w = X[i]
            tokens[t] = idfreq
            maxfreq = max(maxfreq, idfreq.freq)
        end

    else kind == FreqModel
        X = [(t, idfreq) for (t, idfreq) in model.tokens]
        sort!(X, by=x->x[end].freq, rev=true)
        for i in 1:k
            t, idfreq = X[i]
            tokens[t] = idfreq
            maxfreq = max(maxfreq, idfreq.freq)
        end
    end

    id2token = Dict(idfreq.id => t for (t, idfreq) in tokens)
    VectorModel(model.config, tokens, id2token, maxfreq, model.m, model.n)
end

prune_select_top(model::VectorModel, ratio::AbstractFloat, kind=IdfModel) = prune_select_top(model, floor(Int, length(model.tokens) * ratio), kind)

"""
    update!(a::VectorModel, b::VectorModel)

Updates `a` with `b` inplace; returns `a`.
"""
function update!(a::VectorModel, b::VectorModel)
    m = a.m
    for (token, idfreq) in b.tokens
        x = get(a.tokens, token, nothing)
        if x === nothing
            m += 1
            a.tokens[token] = IdFreq(m, x.freq)
        else
            a.tokens[token] = IdFreq(idfreq.id, idfreq.freq + x.freq)
        end
    end

    a.maxfreq = max(a.maxfreq, b.maxfreq)
    a.n += b.n
    a.m = m
    a.id2token = Dict(idfreq.id => t for (t, idfreq) in a.tokens)
    a
end

"""
    vectorize(model::VectorModel, weighting::Type, data; normalize=true)::Dict{Int, Float64}

Computes `data`'s weighted bow of words using the given model and weighting scheme;
the vector is normalized to the unit normed vector if normalize is true
"""
function vectorize(model::VectorModel, weighting::Type, data::DataType; normalize=true) where DataType <: Union{AbstractString, AbstractVector{S}} where S <: AbstractString
    bow, maxfreq = compute_bow(tokenize(model.config, data))
    vectorize(model, weighting, bow, maxfreq, normalize=normalize)
end

"""
    vectorize(model::VectorModel, weighting::Type, bow::BOW, maxfreq=0; normalize=true) where Tv<:Real

Computes a weighted vector using the given bag of words and the specified weighting scheme.
"""
function vectorize(model::VectorModel, weighting::Type, bow::BOW, maxfreq=0; normalize=true)
    if maxfreq == 0
        for v in values(bow)
            maxfreq = max(maxfreq, v)
        end
    end

    vec = SVEC()
    for (token, freq) in bow
        t = get(model.tokens, token, nothing)

        if t === nothing
            continue
        end

        w = _weight(weighting, freq, maxfreq, model.n, t.freq)

        if w > 1e-6
            vec[t.id] = w
        end
    end

    normalize && normalize!(vec)
    vec
end

"""
    vectorize(model::VectorModel, data; normalize=true)

Computes the vector of data using TfidfModel as default
"""
vectorize(model::VectorModel, data; normalize=true) = vectorize(model, TfidfModel, data, normalize=normalize)

function broadcastable(model::VectorModel)
    (model,)
end

"""
    _weight(::Type{T}, freq::Integer, maxfreq::Integer, n::Integer, global_freq::Integer)::Float64

Computes a weight for the given stats using scheme T
"""
function _weight(::Type{TfidfModel}, freq::Real, maxfreq::Real, n::Real, global_freq::Real)::Float64
    (freq / maxfreq) * log(2, 1 + n / global_freq)
end

function _weight(::Type{TfModel}, freq::Real, maxfreq::Real, n::Real, global_freq::Real)::Float64
    freq / maxfreq
end

function _weight(::Type{IdfModel}, freq::Real, maxfreq::Real, n::Real, global_freq::Real)::Float64
    log(2, n / global_freq)
end

function _weight(::Type{FreqModel}, freq::Real, maxfreq::Real, n::Real, global_freq::Real)::Float64
    freq
end


text:
export rats, pumps, seeds, coin

coin = @model flips begin
    pHeads ~ Beta(1,1)
    flips ~ Bernoulli(pHeads) |> iid(20)
end

export lda
lda = @model (α, η, K, V, N) begin
    M = length(N)
    β ~ Dirichlet(repeat([η],V)) |> iid(K)
    θ ~ Dirichlet(repeat([α],K)) |> iid(M)
    z ~ For(1:M) do m
            Categorical(θ[m]) |> iid(N[m])
        end
    w ~ For(1:M) do m
            For(1:N[m]) do n
                Categorical(β[z[m][n]])
            end
        end
end

export hello
hello = @model μ begin
    σ ~ HalfCauchy()
    x ~ Normal(μ,σ)
end

export normalModel
normalModel = @model x begin
    μ ~ Normal(0,5)
    σ ~ HalfCauchy(3)
    x ~ Normal(μ,σ) |> iid(10)
end

# Is this right?
export tdist
tdist = @model ν begin
    w ~ InverseGamma(ν/2,ν/2)
    x ~ Normal(0, w^2)
    return x
end

export nested
nested = @model x begin
    μ ~ simpleModel(s = 2.0)
    σ ~ HalfCauchy(3)
    # x ~ Normal(μ.z,σ) |> iid
end

# export nested2
# nested2 = @model x begin
#     M ~ @model begin
#         z ~ Normal(0,2.0)
#     end
#     σ ~ HalfCauchy(3)
#     x ~ Normal(M.z, σ)
# end

export simpleModel
simpleModel = @model s begin
    z ~ Normal(0,s)
end

export mix
mix = @model (K,α) begin
    p ~ Dirichlet(repeat([α],K))
    μ ~ Cauchy(0,5) |> iid(K)
    σ ~ HalfCauchy(3) |> iid(K)
    components = Normal.(μ,σ)
    N ~ Poisson(100)
    x ~ MixtureModel(components, p) |> iid(N)
end

export linReg1D

linReg1D_local(ŷ, σ) = n -> Normal(ŷ[n])
linReg1D = @model (x, y) begin
        α ~ Cauchy(0, 10)
        β ~ Cauchy(0, 2.5)
        σ ~ HalfCauchy(3)
        ŷ = α .+ β .* x
        N = length(x)
        y ~ For(linReg1D_local(ŷ, σ), 1:N)
    end

# From OpenBUGS and section 6 of Gelfand et al.
rats = @model x begin
    μdist = Normal(0,1000)
    μα ~ μdist
    μβ ~ μdist

    σ2dist = InverseGamma(0.001,0.001)
    σ2α ~ σ2dist
    σ2β ~ σ2dist
    σ2c ~ σ2dist

    α ~ Normal(μα, sqrt(σ2α)) |> iid(30)
    β ~ Normal(μβ, sqrt(σ2β)) |> iid(30)

    x̄ = mean(x)
    y ~ For([(i,j) for i in 1:30, j in 1:5]) do (i,j)
            Normal(α[i] + β[i]*(x[j]-x̄), sqrt(σ2c))
        end
end

pumps = @model t begin
    n = length(t)
    α ~ Gamma(1,1)
    β ~ Gamma(0.1,1)
    θ ~ For(1:n) do i
            Gamma(α, 1/β)
        end
    y ~ For(1:n) do i
            Poisson(θ[i] * t[i])
        end
end

# dogs = @model
#     lo = typemin(Float64) |> nextfloat
#     negReals = Uniform(lo,-1e-5)
#     α ~ negReals
#     β ~ negReals
#     y ~ For()
#     Bernoulli
#     α
# end

seeds = @model (n,x) begin
    σ2 ~ InverseGamma(0.001,0.001)
    σ = sqrt(σ2)
    αDist = Normal(0,1000)
    α0 ~ αDist
    α1 ~ αDist
    α2 ~ αDist
    α12 ~ αDist
    b ~ Normal(0,σ) |> iid(21)
    y = α0 .+ α1 .* x[1,:] .+ α2 .* x[2,:] .+ α12 .* x[1,:] .* x[2,:] .+ b
    r ~ For(1:21) do i
            LogisticBinomial(n[i],y[i])
        end
end

stacks = @model begin
    σ2 ~ InverseGamma(0.001,0.001)
    βDist = Normal(0,1000)
    β0 ~ βDist
    β1 ~ βDist
    β2 ~ βDist
    β3 ~ βDist
    μ = β0 + β1*z[1,:] + β2*z[2,:] + β3*z[3,:]
    y ~ For(1:21) do i
            Laplace(μ[i],σ2)
        end
end

y = [28.,  8., -3.,  7., -1.,  1., 18., 12.]
σ = [15., 10., 16., 11.,  9., 11., 10., 18.]

export school8
school8 = @model (y, σ) begin
  μ ~ Flat()
  τ ~ HalfFlat()
  η ~ Normal() |> iid(8)
  y ~ For(1:8) do j
      Normal(μ + τ * η[j], σ[j])
  end
end

export funnel
funnel = @model begin
    y ~ Normal(0, 3);
    x ~ Normal(0, exp(y/2))
end


text:
const MAX_SIZEOF_MANTISSA = 20

struct OraNumber <: Real
    len::UInt8
    ex::UInt8
    mantissa::NTuple{MAX_SIZEOF_MANTISSA, UInt8}
end


text:
include("../../../src/Dexter.jl")
using Dexter

using LightGraphs
using Gadfly
using GraphPlot
using JuMP
using Clp

m = Model(solver = ClpSolver())
@variable(m, 0 <= x <= 5)
@variable(m, 0 <= y <= 5)
@variable(m, 0 <= z <= 5)

@objective(m, Max, 5x + 3y + 2z)
# @NLconstraint(m, cos(x)/y * x * e^(-x^2) <= 3.0 )
# @NLconstraint(m, cos(x) * e^(-x^2) + y - x*x + 1 <= 5.0 - x^2)
# @NLconstraint(m, cos(x) * e^(-x^2) + y - x*x + -(x/2 + 2*y) + 1 <= x^2 + 5)
# @NLconstraint(m, 2x <= 3.0 )
# @NLconstraint(m, x + y + z <= 8)
# @NLconstraint(m, - 2 * x * e^(-x) - 5 * y * e^(-y^2) - 3 * z * e^(-z) <= 30)
@NLconstraint(m, 2 * x * e^(-x) + 5 * y * e^(-y^2) + 3 * z * e^(-z) <= 3)

print(m)

# TODO: Right now we are assuming <=. Consider >= case.

# Compute evaluations of parts of additively separable constraints
evals = Dexter.evaluate_separable_constraints(m)

# TODO: Transition function should support dynamic ordering

# Set ordering and objective
n = MathProgBase.numvar(m)
ordering = Dexter.NoOrdering()

# Create specs for decision diagram
eval = evals[1]
initial_state = Dexter.create_constraint_initial_state(eval)
domains = Dexter.create_constraint_domain_function(m)
transition = Dexter.create_constraint_transition_function(eval, domains, ordering)
mathprog_specs = Dexter.ProblemSpecs(initial_state, transition, domains)

# Construct decision diagram
dd = Dexter.construct_DD(n, mathprog_specs, ordering=ordering)

# Plot decision diagram
locs_x = [Dexter.get_node_idlayer(dd, i) * 1.0 for i in 1:nv(dd.graph)]
locs_y = [Dexter.get_node_layer(dd, i) / n for i in 1:nv(dd.graph)]
gplot(dd.graph, locs_x, locs_y) |> PNG("test_mp.png", 16cm, 16cm)

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
@testset "Thermodynamics" begin
    @test true == true
end

max_stars_repo_name: VictorVeraFrazao/OptStoic

max_issues_repo_name: VictorVeraFrazao/OptStoic

max_forks_repo_name: VictorVeraFrazao/OptStoic


text:
function print_node(io::IO, b::LazyBridgeOptimizer, node::VariableNode)
    S, = b.variable_types[node.index]
    MOIU.print_with_acronym(io, "[$(node.index)] constrained variables in `$S` are")
end
function print_node(io::IO, b::LazyBridgeOptimizer, node::ConstraintNode)
    F, S = b.constraint_types[node.index]
    MOIU.print_with_acronym(io, "($(node.index)) `$F`-in-`$S` constraints are")
end
function print_node(io::IO, b::LazyBridgeOptimizer, node::ObjectiveNode)
    F, = b.objective_types[node.index]
    MOIU.print_with_acronym(io, "|$(node.index)| objective function of type `$F` is")
end
function print_node_info(io::IO, b::LazyBridgeOptimizer, node::AbstractNode)
    print(io, " ")
    print_node(io, b, node)
    d = _dist(b.graph, node)
    if d == INFINITY
        print(io, " not supported\n")
    else
        index = bridge_index(b.graph, node)
        if iszero(index) || (node isa VariableNode && !is_variable_edge_best(b.graph, node))
            @assert node isa VariableNode
            println(io, " supported (distance $d) by adding free variables and then constrain them, see ($(b.graph.variable_constraint_node[node.index].index)).")
        else
            print(io, " bridged (distance $d) by ")
            MOIU.print_with_acronym(io, string(_bridge_type(b, node, index)))
            println(io, ".")
        end
    end
end
print_graph(b::LazyBridgeOptimizer) = print_graph(Base.stdout, b)
function print_graph(io::IO, b::LazyBridgeOptimizer)
    println(io, b.graph)
    for node in variable_nodes(b.graph)
        print_node_info(io, b, node)
    end
    for node in constraint_nodes(b.graph)
        print_node_info(io, b, node)
    end
    for node in objective_nodes(b.graph)
        print_node_info(io, b, node)
    end
end

function print_unsupported(io::IO, b::LazyBridgeOptimizer, node::AbstractNode)
    if _dist(b.graph, node) != INFINITY
        return
    end
    print(io, "  ")
    print_node(io, b, node)
    print(io, " not supported\n")
end
function print_unsupported(io::IO, b::LazyBridgeOptimizer, edge::Edge)
    for node in edge.added_variables
        print_unsupported(io, b, node)
    end
    for node in edge.added_constraints
        print_unsupported(io, b, node)
    end
end
function print_unsupported(io::IO, b::LazyBridgeOptimizer, edge::ObjectiveEdge)
    for node in edge.added_variables
        print_unsupported(io, b, node)
    end
    for node in edge.added_constraints
        print_unsupported(io, b, node)
    end
    print_unsupported(io, b, edge.added_objective)
end
function print_unsupported(io::IO, b::LazyBridgeOptimizer, edges::Vector, index_to_bridge::Function)
    no_bridge = true
    for edge in edges
        if no_bridge
            println(io, ":")
            no_bridge = false
        end
        MOIU.print_with_acronym(io, "  Cannot use `$(index_to_bridge(edge.bridge_index))` because:\n")
        print_unsupported(io, b, edge)
    end
    if no_bridge
        println(io, " no added bridge supports bridging it.")
    end
end
function _bridge_type(b::LazyBridgeOptimizer, node::VariableNode, bridge_index::Int)
    return Variable.concrete_bridge_type(b.variable_bridge_types[bridge_index], b.variable_types[node.index]...)
end
function _bridge_type(b::LazyBridgeOptimizer, node::ConstraintNode, bridge_index::Int)
    return Constraint.concrete_bridge_type(b.constraint_bridge_types[bridge_index], b.constraint_types[node.index]...)
end
function _bridge_type(b::LazyBridgeOptimizer, node::ObjectiveNode, bridge_index::Int)
    return Objective.concrete_bridge_type(b.objective_bridge_types[bridge_index], b.objective_types[node.index]...)
end
function print_unsupported(io::IO, b::LazyBridgeOptimizer, variables, constraints, objectives)
    for node in variables
        print_node(io, b, node)
        print(io, " not supported because")
        print_unsupported(io, b, b.graph.variable_edges[node.index],
                          bridge_index -> _bridge_type(b, node, bridge_index))
        print(io, "  Cannot add free variables and then constrain them because")
        constraint_node = b.graph.variable_constraint_node[node.index]
        if constraint_node.index == -1
            println(io, " free variables are bridged but no functionize bridge was added.")
        else
            println(io, ":")
            print_unsupported(io, b, constraint_node)
        end
    end
    for node in constraints
        print_node(io, b, node)
        print(io, " not supported because")
        print_unsupported(io, b, b.graph.constraint_edges[node.index],
                          bridge_index -> _bridge_type(b, node, bridge_index))
    end
    for node in objectives
        print_node(io, b, node)
        print(io, " not supported because")
        print_unsupported(io, b, b.graph.objective_edges[node.index],
                          bridge_index -> _bridge_type(b, node, bridge_index))
    end
end

function add_unsupported(graph::Graph, edges::Vector{Edge},
                         variables, constraints, objectives)
    for edge in edges
        for node in edge.added_variables
            add_unsupported(graph, node, variables, constraints, objectives)
        end
        for node in edge.added_constraints
            add_unsupported(graph, node, variables, constraints, objectives)
        end
    end
end
function add_unsupported(graph::Graph, edges::Vector{ObjectiveEdge},
                         variables, constraints, objectives)
    for edge in edges
        for node in edge.added_variables
            add_unsupported(graph, node, variables, constraints, objectives)
        end
        for node in edge.added_constraints
            add_unsupported(graph, node, variables, constraints, objectives)
        end
        add_unsupported(graph, edge.added_objective, variables, constraints, objectives)
    end
end
function add_unsupported(graph::Graph, node::VariableNode,
                         variables, constraints, objectives)
    if _dist(graph, node) != INFINITY
        return
    end
    push!(variables, node)
    add_unsupported(graph, graph.variable_edges[node.index], variables, constraints, objectives)
    constraint_node = graph.variable_constraint_node[node.index]
    if constraint_node.index != -1
        add_unsupported(graph, graph.variable_constraint_node[node.index], variables, constraints, objectives)
    end
end
function add_unsupported(graph::Graph, node::ConstraintNode,
                         variables, constraints, objectives)
    if _dist(graph, node) != INFINITY
        return
    end
    push!(constraints, node)
    add_unsupported(graph, graph.constraint_edges[node.index], variables, constraints, objectives)
end
function add_unsupported(graph::Graph, node::ObjectiveNode,
                         variables, constraints, objectives)
    if _dist(graph, node) != INFINITY
        return
    end
    push!(objectives, node)
    add_unsupported(graph, graph.objective_edges[node.index], variables, constraints, objectives)
end
function _sort(nodes::Set)
    vector = collect(nodes)
    sort!(vector, by = node -> node.index)
    return vector
end
function debug_unsupported(io::IO, b::LazyBridgeOptimizer, node::AbstractNode)
    variables = Set{VariableNode}()
    constraints = Set{ConstraintNode}()
    objectives = Set{ObjectiveNode}()
    add_unsupported(b.graph, node, variables, constraints, objectives)
    print_unsupported(io, b, _sort(variables), _sort(constraints), _sort(objectives))
end
function debug(b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
               S::Type{<:MOI.AbstractSet}; io::IO = Base.stdout)
    if MOI.supports_constraint(b, F, S)
        if F == MOIU.variable_function_type(S)
            # This may be thanks to a variable bridge so `F`-in-`S` constraints
            # are maybe not supported but constrained variables in `S` are
            # definitely supported.
            MOIU.print_with_acronym(io, "Constrained variables in `$S` are supported.\n")
        else
            MOIU.print_with_acronym(io, "`$F`-in-`$S` constraints are supported.\n")
        end
    else
        message = " are not supported and cannot be bridged into supported" *
            " constrained variables and constraints. See details below:"
        if F == MOIU.variable_function_type(S)
            MOIU.print_with_acronym(io, "Constrained variables in `$S`")
            println(io, message)
            debug_unsupported(io, b, node(b, S))
        else
            MOIU.print_with_acronym(io, "`$F`-in-`$S` constraints")
            println(io, message)
            debug_unsupported(io, b, node(b, F, S))
        end
    end
end

"""
    debug_supports_constraint(
        b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
        S::Type{<:MOI.AbstractSet}; io::IO = Base.stdout)

Prints to `io` explanations for the value of [`MOI.supports_constraint`](@ref)
with the same arguments.
"""
function debug_supports_constraint(
    b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet}; kws...)
    debug(b, F, S; kws...)
end

function debug(b::LazyBridgeOptimizer,
               F::Type{<:MOI.AbstractScalarFunction}; io::IO = Base.stdout)
    MOIU.print_with_acronym(io, "Objective function of type `$F` is")
    if supports_bridging_objective_function(b, F)
        MOIU.print_with_acronym(io, " supported.\n")
    else
        MOIU.print_with_acronym(io, " not supported and cannot be bridged" *
            " into a supported objective function by adding only supported" *
            " constrained variables and constraints. See details below:\n")
        debug_unsupported(io, b, node(b, F))
    end
end

"""
    debug_supports(b::LazyBridgeOptimizer, ::MOI.ObjectiveFunction{F}; io::IO = Base.stdout) where F

Prints to `io` explanations for the value of [`MOI.supports`](@ref) with the
same arguments.
"""
function debug_supports(
    b::LazyBridgeOptimizer, ::MOI.ObjectiveFunction{F}; kws...) where F
    debug(b, F; kws...)
end


text:
function test_pentagon_ispentagonal()
    @test find([ispentagonal(i) for i in 1:1000]) == [1, 5, 12, 22, 35, 51, 70, 92, 117, 145, 176, 210, 247, 287, 330, 376, 425, 477, 532, 590, 651, 715, 782, 852, 925]
end

function test_pentagon_nthpentagonal()
    @test nthpentagonal(1) == 1
    @test nthpentagonal(2) == 5
    @test nthpentagonal(10) == 145
end

function test_pentagon_npentagonal()
    @test npentagonal(10) == [1, 5, 12, 22, 35, 51, 70, 92, 117, 145]
    @test eltype(npentagonal(10)) == Int
    @test eltype(npentagonal(10, BigInt)) == BigInt
end

function test_pentagon_somepentagonal()
    @test collect(somepentagonal(91)) == [1, 5, 12, 22, 35, 51, 70]
    @test collect(somepentagonal(92)) == [1, 5, 12, 22, 35, 51, 70, 92]
    @test eltype(somepentagonal(92)) == Int
    @test eltype(somepentagonal(big(92))) == BigInt
end

function test_pentagon_exactpentagonal()
    @test collect(exactpentagonal(10)) == [1, 5, 12, 22, 35, 51, 70, 92, 117, 145]
    @test eltype(exactpentagonal(10)) == Int
    @test eltype(exactpentagonal(10, BigInt)) == BigInt
end

function test_pentagon_all()
    print(rpad("Math.NumberTheory.Pentagon...", 50, ' '))

    test_pentagon_ispentagonal()
    test_pentagon_nthpentagonal()
    test_pentagon_npentagonal()
    test_pentagon_somepentagonal()
    test_pentagon_exactpentagonal()

    println("PASS")
end


text:
#!/usr/bin/env julia
#=
Package Gnar
del2z <delta.z@aliyun.com>
=#

function loadwiki(fname::String)
    fin = open(fname, "r")
    doclist = Vector{String}[]
    count = 0
    doc = Vector{String}()
    IsKeep = false
    for line in eachline(fin)
        line = strip(line)
        if line == ""
            continue
        elseif startswith(line, "<doc id")
            IsKeep = true
            title = match(r"title=\"(.*)\"", line).captures[1]
            if startswith(title, "Wikipedia:") || startswith(title, "Help:")
                isKeep = false
                println("Exclude content of \"$title\"")
            end
        elseif startswith(line, "</doc>")
            if IsKeep && length(doc) > 0
                push!(doclist, copy(doc))
            end
            IsKeep = false
            empty!(doc)
        elseif IsKeep && isvalid(line)
            if startswith(line, "，") || startswith(line, "。") || (length(doc) > 0 && endswith(doc[end], "，"))
                doc[end] = doc[end] * clean(line)
            else
                push!(doc, lstrip(clean(line), ['，', '。']))
            end
        end
        count += 1
        if count > 5000
            break
        end
    end
    close(fin)
    doclist
end

function isvalid(str::AbstractString)
    reZh = r"[\u4e00-\u9fd5\uff0c]{8}"
    return match(reZh, str) != nothing && !startswith(str, "[[Category:") && !startswith(str, "[[File:")
end

function clean(str::AbstractString)
    str = replace(str, "（，" => "（")
    str = replace(str, "，）" => "）")
    str = replace(str, "「" => "“")
    str = replace(str, "」" => "”")
    str = replace(str, "“”" => "")
    str = replace(str, "-{}-" => "")
    str = replace(str, "《》" => "")
    str = replace(str, "（）" => "")
    str = replace(str, "()" => "")
    reCh = "\\u4e00-\\u9fd5a-zA-Z0-9·"
    str = replace(str, Regex("-{([$reCh]+)}-") => s"\1")
    if match(Regex("-{[$reCh;:\\- ]*zh-cn:[$reCh;:\\- ]*}-"), str) != nothing
        str = replace(str, Regex("-{[$reCh;:\\- ]*zh-cn:([$reCh]+)[$reCh;:\\- ]*}-") => s"\1")
    end
    if match(Regex("-{[$reCh;:\\- ]*zh-hans:[$reCh;:\\- ]*}-"), str) != nothing
        str = replace(str, Regex("-{[$reCh;:\\- ]*zh-hans:([$reCh]+)[$reCh;:\\- ]*}-") => s"\1")
    end
    if match(Regex("\\[\\[[$reCh]+\\]\\]"), str) != nothing
        str = replace(str, Regex("\\[\\[([$reCh]+)\\]\\]") => s"\1")
    end
    if match(Regex("\\[\\[[$reCh\\(\\) ]+\\|[$reCh\\(\\) ]+\\]\\]"), str) != nothing
        str = replace(str, Regex("\\[\\[[$reCh\\(\\) ]+\\|([$reCh\\(\\) ]+)\\]\\]") => s"\1")
    end
    str = replace(str, r"[。]{2,}" => "。")
    str = replace(str, r"[；]{2,}" => "；")
    # codice_, formula_
    str = lstrip(str, ['：', '？'])
    str = rstrip(str, ['：'])
    str
end

doclist = loadwiki("../../../../Resources/Corpus/zh_wiki.txt")
for doc in doclist
    println([x for x in doc if occursin("深化合作，研究和制定粤港澳大湾区", x)])
end


text:
export iauUt1tai
"""
Time scale transformation:  Universal Time, UT1, to International
Atomic Time, TAI.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   ut11,ut12  double    UT1 as a 2-part Julian Date
   dta        double    UT1-TAI in seconds

Returned:
   tai1,tai2  double    TAI as a 2-part Julian Date

Returned (function value):
               int       status:  0 = OK

Notes:

   1. ut11+ut12 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where ut11 is the Julian
   Day Number and ut12 is the fraction of a day.  The returned
   tai1,tai2 follow suit.

   2. The argument dta, i.e. UT1-TAI, is an observed quantity, and is
   available from IERS tabulations.

Reference:

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauUt1tai(double ut11, double ut12, double dta,
#               double *tai1, double *tai2)

function iauUt1tai(ut11::Real, ut12::Real, dta::Real)
   ref_tai1 = Ref{Float64}(0.0)
   ref_tai2 = Ref{Float64}(0.0)

   status = ccall((:iauUt1tai, libsofa_c), Cint,
                   (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, ut11), convert(Float64, ut12), convert(Float64, dta),
                   ref_tai1, ref_tai2)

   return status, ref_tai1[], ref_tai2[]
end


text:
module Butterfly

end # module


text:
using DiffEqFlux, OrdinaryDiffEq, Flux, CUDA
CUDA.allowscalar(false) # Makes sure no slow operations are occuring

# Generate Data
u0 = Float32[2.0; 0.0]
datasize = 30
tspan = (0.0f0, 1.5f0)
tsteps = range(tspan[1], tspan[2], length = datasize)
function trueODEfunc(du, u, p, t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
prob_trueode = ODEProblem(trueODEfunc, u0, tspan)
# Make the data into a GPU-based array if the user has a GPU
ode_data = gpu(solve(prob_trueode, Tsit5(), saveat = tsteps))

dudt2 = FastChain((x, p) -> x.^3,
                  FastDense(2, 50, tanh),
                  FastDense(50, 2))
u0 = Float32[2.0; 0.0] |> gpu
p = initial_params(dudt2) |> gpu
prob_neuralode = NeuralODE(dudt2, tspan, Tsit5(), saveat = tsteps)

function predict_neuralode(p)
  gpu(prob_neuralode(u0,p))
end
function loss_neuralode(p)
    pred = predict_neuralode(p)
    loss = sum(abs2, ode_data .- pred)
    return loss, pred
end
# Callback function to observe training
list_plots = []
iter = 0
callback = function (p, l, pred; doplot = false)
  global list_plots, iter
  if iter == 0
    list_plots = []
  end
  iter += 1
  display(l)
  # plot current prediction against data
  plt = scatter(tsteps, Array(ode_data[1,:]), label = "data")
  scatter!(plt, tsteps, Array(pred[1,:]), label = "prediction")
  push!(list_plots, plt)
  if doplot
    display(plot(plt))
  end
  return false
end
result_neuralode = DiffEqFlux.sciml_train(loss_neuralode, p,
                                          ADAM(0.05), cb = callback,
                                          maxiters = 300)


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

V,(VV,EV,FV,CV) = Lar.cuboidGrid([10,10,5],true)

GL.VIEW([
    GL.GLGrid(V,FV,GL.Point4d(1,1,1,0.1))
	GL.GLAxis(GL.Point3d(-1,-1,-1),GL.Point3d(1,1,1))
]);

GL.VIEW([
	GL.GLGrid(V,VV,GL.COLORS[4])
	GL.GLGrid(V,EV,GL.COLORS[2])
	GL.GLGrid(V,FV,GL.Point4d(0,1,0,0.1))
	GL.GLAxis(GL.Point3d(-1,-1,-1),GL.Point3d(1,1,1))
]);

# visualization of a simplicial 3-complex (CV are all tetrahedra)
V,CV = Lar.simplexGrid([4,4,4])
FV = Lar.simplexFacets(CV)
GL.VIEW([
	GL.GLGrid(V,FV,GL.Point4d(1,1,1,0.1))
	GL.GLAxis(GL.Point3d(0,0,0),GL.Point3d(1,1,1))
]);

# visualization of a simplicial 2-complex (TV are all triangles)
V,CV = Lar.simplexGrid([4,4,4])
FV = Lar.simplexFacets(CV)
TV = Lar.simplexBoundary_3(CV,FV)
EV = Lar.simplexFacets(TV)

GL.VIEW([
	GL.GLGrid(V,TV,GL.Point4d(1,1,1,0.1))
	GL.GLGrid(V,EV,GL.Point4d(1,1,1,0.25))
	GL.GLAxis(GL.Point3d(-1,-1,-1),GL.Point3d(1,1,1))
]);


text:
# Contains FIS and Rule types
# ---------------------------
#
'''
    # A Rule connecting input and output membership function
    #
    # Properties
    # ----------
    # `input_mf_names` is a Vector of AbstractString (in order of inputs) which contains name of membership functions
    # `output_mf` is either an AbstractString containing the connected output membership function name with `input_mf_names`
    # 		or is a Vector of AbstractFloat containing consequence parameters in case of Sugeno fis
'''
struct Rule
    input_mf_names::Vector{AbstractString}
    output_mf::Union{Vector{Float64},AbstractString}
    firing_method::AbstractString

    function Rule(input_mf_names, output_mf, firing_method="MIN")
        new(input_mf_names, output_mf, firing_method)
    end
end

'''
    # Fuzzy Inference System of Mamdani type
    #
    # Properties
    # ----------
    # `input_mfs_dicts` is a Vector of Dict (one for each input) which contains membership function
    # 		associated with its name
    # 		e.g. "small" => TriangularMF(1, 2, 3)
    # `output_mfs_dict` is a Dict containing membership function of output associated with names
    # 		e.g. "small" => GaussianMF(4, 2)
    # `rules` is a Vector of Rule that form the inference system
'''
struct FISMamdani
    input_mfs_dicts::Vector{Dict{AbstractString,MF}}
    output_mfs_dict::Dict{AbstractString,MF}
    rules::Vector{Rule}
end

'''
    # Fuzzy Inference System of Mamdani type
    #
    # Properties
    # ----------
    # `input_mfs_dicts` is a Vector of Dict (one for each input) which contains membership function
    # 		associated with its name
    # 		e.g. "small" => TriangularMF(1, 2, 3)
    # `rules` is a Vector of Rule (of sugeno type) that form the inference system
'''
struct FISSugeno
    input_mfs_dicts::Vector{Dict{AbstractString,MF}}
    rules::Vector{Rule}
end


text:
push!(LOAD_PATH, ".")

module GeneticAlgorithm

using Population
using Individual
using BinaryIndividual
using IntegerIndividual
using RealIndividual
using Fitness
using GAInitialization
using GASelection
using GACrossover
using GAMutation
using Debug
using Statistics
using Plots

include("utility.jl")

abstract type AbstractGeneticAlgorithm end

mutable struct GeneticAlgorithmType
  pop::Population.PopulationType{<: Individual.AbstractIndividual}
  pop_size::Integer
  best_solution::Vector{Population.IndFitType{<: Individual.AbstractIndividual}}
  elite_size::Integer
  init_args::Tuple
  sel_args::Tuple
  cross_args::Tuple
  mut_args::Tuple
  function GeneticAlgorithmType(x...;
                        init_args::Tuple = (GAInitialization.InitRandom,),
                        sel_args::Tuple = (GASelection.Tournament, 2),
                        cross_args::Tuple = (GACrossover.PointCrossover, 0.95),
                        mut_args::Tuple = (GAMutation.BitFlipMutation, 0.01))
    args = build(x..., init_args, sel_args, cross_args, mut_args)
    new(args[1], args[2], args[3], args[4], args[5], args[6], args[7], args[8])
  end
end

function build(ind_args::Tuple, fit::Fitness.AbstractFitness, ps::Integer, es::Integer,
               i_a::Tuple, s_a::Tuple, c_a::Tuple, m_a::Tuple)
  ps > 1 || error("GA: Population size must be greater than 1")
  return Population.PopulationType{ind_args[1]}(ind_args[2:end], fit), ps, Vector{Population.IndFitType{ind_args[1]}}(), es, i_a, s_a, c_a, m_a
end

function getBestSolution(this::GeneticAlgorithmType)
  return this.best_solution
end

function evolveSO!(this::GeneticAlgorithmType, num_it::Integer, log::Integer = 0)
  best_fitness = zeros(Float64, num_it)
  mean_fitness = zeros(Float64, num_it)
  if Debug.ga_plot
    gr()
    plt = plot(zeros(0, 2), zeros(0, 2), xlims=(1, num_it), label=["Best" "Mean"], xlabel="Generation", ylabel="Fitness", legend=:bottomright, show=true)
  end

  # Create initial population
  GAInitialization.initializePopulation!(this.pop, this.pop_size, this.init_args...)

  # Main loop
  @time for it = 1 : num_it
    if Debug.ga_debug
      println("----- Generation ", it, " -----\n")
      show(this.pop)
    end

    ind_args = Population.getIndArgs(this.pop)
    pop_fit = Population.getFitFunction(this.pop)
    new_pop = typeof(this.pop)(ind_args, pop_fit)

    # Population sorting
    cur_individuals = collect(this.pop)
    sort!(cur_individuals; rev = true)

    # Selection
    child_num = this.pop_size - this.elite_size
    parents_group = GASelection.selectParents(cur_individuals, child_num, this.sel_args...)

    # Crossover
    childs = typeof(this.pop).parameters[1][]
    for parents in parents_group
      res_ch = GACrossover.crossover(parents..., this.cross_args...)
      push!(childs, res_ch...)
      length(childs) < child_num || break
    end
    childs = childs[1:child_num]

    # Mutation
    for child in childs
      Population.insertIndividual!(new_pop, GAMutation.mutate!(child, this.mut_args...))
    end

    # Elitism
    for (ind, fit) in @view cur_individuals[1:this.elite_size]
      Population.insertIndividual!(new_pop, ind, fit)
    end

    # New generation
    Population.evalFitness!(new_pop)
    this.pop = new_pop
    cur_best_solution = maximum(this.pop)
    if isempty(this.best_solution) || cur_best_solution > this.best_solution[1]
      this.best_solution = [deepcopy(cur_best_solution)]
    end
    if log > 0 && (it - 1) % log == 0
      println(it, " -> Fitness: (", cur_best_solution[2], ", ", this.best_solution[1][2], ")")
    end

    # Statistics
    best_fitness[it] = cur_best_solution[2][1]
    mean_fitness[it] = mean([fit[1] for (ind, fit) in this.pop])
    if Debug.ga_plot
      push!(plt, [it], [best_fitness[it], mean_fitness[it]])
      gui()
    end
  end
  return best_fitness, mean_fitness
end

function evolveMO!(this::GeneticAlgorithmType, num_it::Integer, log::Integer = 0)
  # Create initial population
  GAInitialization.initializePopulation!(this.pop, this.pop_size, this.init_args...)

  # Main loop
  @time for it = 1 : num_it
    if Debug.ga_debug
      println("----- Generation ", it, " -----\n")
      show(this.pop)
    end

    ind_args = Population.getIndArgs(this.pop)
    pop_fit = Population.getFitFunction(this.pop)
    new_pop = typeof(this.pop)(ind_args, pop_fit)

    # Population sorting
    cur_individuals = [(this.pop[i], i) for i in eachindex(this.pop)]
    sorted_individuals = Array{Int64, 1}[]
    while !isempty(cur_individuals)
      dominated = fill(false, length(cur_individuals))
      frontier = Int64[]
      for i in eachindex(cur_individuals)
        for j in eachindex(cur_individuals)
          if cur_individuals[j][1] > cur_individuals[i][1]
            dominated[i] = true
            break
          end
        end
        if !dominated[i]
          push!(frontier, cur_individuals[i][2])
        end
      end
      cur_individuals = cur_individuals[filter(x -> dominated[x], eachindex(cur_individuals))]
      push!(sorted_individuals, frontier)
    end
    cur_individuals = [(this.pop[ind][1], (Float64(-i),))
                       for i in eachindex(sorted_individuals) for ind in sorted_individuals[i]]

    # Selection
    child_num = this.pop_size - this.elite_size
    parents_group = GASelection.selectParents(cur_individuals, child_num, this.sel_args...)

    # Crossover
    childs = typeof(this.pop).parameters[1][]
    for parents in parents_group
      res_ch = GACrossover.crossover(parents..., this.cross_args...)
      push!(childs, res_ch...)
      length(childs) < child_num || break
    end
    childs = childs[1:child_num]

    # Mutation
    for child in childs
      Population.insertIndividual!(new_pop, GAMutation.mutate!(child, this.mut_args...))
    end

    # Elitism
    for ind in @view vcat(sorted_individuals...)[1:this.elite_size]
      Population.insertIndividual!(new_pop, this.pop[ind]...)
    end

    # New generation
    Population.evalFitness!(new_pop)
    this.pop = new_pop
    cur_best_solution = typeof(this.best_solution)()
    for ind in this.pop
      dominated = false
      for ind2 in this.pop
        if ind2 > ind
          dominated = true
          break
        end
      end
      if !dominated
        push!(cur_best_solution, ind)
      end
    end
    this.best_solution = dominationIntersection(this.best_solution, cur_best_solution, Population.getPopSize(this.pop))
    if log > 0 && (it - 1) % log == 0
      println(it, " -> Fitness: (", length(this.best_solution), ")")
    end
  end
  return this.best_solution
end

function dominationIntersection(x::Array{Population.IndFitType{IndType}}, y::Array{Population.IndFitType{IndType}}, max_size::Integer) where {IndType}
  x_dominated = fill(false, length(x))
  y_dominated = fill(false, length(y))

  for indx in eachindex(x), indy in eachindex(y)
    if x[indx] > y[indy]
      y_dominated[indy] = true
    elseif y[indy] > x[indx]
      x_dominated[indx] = true
    end
  end

  x_non_dominated = x[filter(x -> !x_dominated[x], eachindex(x))]
  y_non_dominated = y[filter(x -> !y_dominated[x], eachindex(y))]
  res = [x_non_dominated..., y_non_dominated...]
  return res[1:min(length(res), max_size)]
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
using Qaintessent

function Qaintessent.backward(g::RNearbyValuesMixerGate, Δ::AbstractMatrix)
    delta = 1e-10

    # uses conjugated gradient matrix
    U_rnv1 = Qaintessent.matrix(RNearbyValuesMixerGate(-(g.β[] - delta/2), g.r, g.d))
    U_rnv2 = Qaintessent.matrix(RNearbyValuesMixerGate(-(g.β[] + delta/2), g.r, g.d))

    U_deriv = (U_rnv2 - U_rnv1) / delta

    return RNearbyValuesMixerGate(2 * sum(real(U_deriv .* Δ)), g.r, g.d)
end

function Qaintessent.backward(g::ParityRingMixerGate, Δ::AbstractMatrix)
    delta = 1e-10

    U_parity1 = Qaintessent.matrix(ParityRingMixerGate(-(g.β[] - delta/2), g.d))
    U_parity2 = Qaintessent.matrix(ParityRingMixerGate(-(g.β[] + delta/2), g.d))

    U_deriv = (U_parity2 - U_parity1) / delta

    return ParityRingMixerGate(2 * sum(real(U_deriv .* Δ)), g.d)
end

function Qaintessent.backward(g::PartitionMixerGate, Δ::AbstractMatrix)
    delta = 1e-10

    U_partition1 = partition_mixer_gate_matrix(g, -(g.β[] - delta/2))
    U_partition2 = partition_mixer_gate_matrix(g, -(g.β[] + delta/2))

    U_deriv = (U_partition2 - U_partition1) / delta

    return PartitionMixerGate(2 * sum(real(U_deriv .* Δ)), g.d, g.partition)
end

function Qaintessent.backward(g::MaxKColSubgraphPhaseSeparationGate, Δ::AbstractMatrix)
    H_P_enc = max_k_col_subgraph_phase_separation_hamiltonian(g.graph, g.κ)

    # we can exploit that H_P_enc is diagonal
    # uses conjugated gradient matrix
    U_deriv = im * H_P_enc * exp(im * g.γ[] * H_P_enc)

    return MaxKColSubgraphPhaseSeparationGate(2 * sum(real(U_deriv .* Δ)), g.κ[], g.graph)
end


text:
insts = OrderedDict()
insts["minimal"] = [
    ((3, false, false),),
    ((3, false, true),),
    ((3, true, false),),
    ((3, true, true),),
    ]
insts["fast"] = [
    ((20, false, false),),
    ((20, true, false),),
    ((8, false, true),),
    ((5, true, true),),
    ((50, false, false),),
    ((50, true, false),),
    ((12, false, true),),
    ((8, true, true),),
    ]
insts["various"] = [
    ((100, false, false),),
    ((100, true, false),),
    ((12, false, true),),
    ((8, true, true),),
    ((200, false, false),),
    ((200, true, false),),
    ((15, false, true), nothing, (default_tol_relax = 100,)),
    ((300, false, false),),
    ((250, true, false),),
    ((17, false, true), nothing, (default_tol_relax = 1000,)),
    ]
return (ClassicalQuantum, insts)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""TokenRequestSpec contains client provided parameters of a token request.

    IoK8sApiAuthenticationV1TokenRequestSpec(;
        audiences=nothing,
        boundObjectRef=nothing,
        expirationSeconds=nothing,
    )

    - audiences::Vector{String} : Audiences are the intendend audiences of the token. A recipient of a token must identitfy themself with an identifier in the list of audiences of the token, and otherwise should reject the token. A token issued for multiple audiences may be used to authenticate against any of the audiences listed but implies a high degree of trust between the target audiences.
    - boundObjectRef::IoK8sApiAuthenticationV1BoundObjectReference : BoundObjectRef is a reference to an object that the token will be bound to. The token will only be valid for as long as the bound object exists. NOTE: The API server's TokenReview endpoint will validate the BoundObjectRef, but other audiences may not. Keep ExpirationSeconds small if you want prompt revocation.
    - expirationSeconds::Int64 : ExpirationSeconds is the requested duration of validity of the request. The token issuer may return a token with a different validity duration so a client needs to check the 'expiration' field in a response.
"""
mutable struct IoK8sApiAuthenticationV1TokenRequestSpec <: SwaggerModel
    audiences::Any # spec type: Union{ Nothing, Vector{String} } # spec name: audiences
    boundObjectRef::Any # spec type: Union{ Nothing, IoK8sApiAuthenticationV1BoundObjectReference } # spec name: boundObjectRef
    expirationSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: expirationSeconds

    function IoK8sApiAuthenticationV1TokenRequestSpec(;audiences=nothing, boundObjectRef=nothing, expirationSeconds=nothing)
        o = new()
        validate_property(IoK8sApiAuthenticationV1TokenRequestSpec, Symbol("audiences"), audiences)
        setfield!(o, Symbol("audiences"), audiences)
        validate_property(IoK8sApiAuthenticationV1TokenRequestSpec, Symbol("boundObjectRef"), boundObjectRef)
        setfield!(o, Symbol("boundObjectRef"), boundObjectRef)
        validate_property(IoK8sApiAuthenticationV1TokenRequestSpec, Symbol("expirationSeconds"), expirationSeconds)
        setfield!(o, Symbol("expirationSeconds"), expirationSeconds)
        o
    end
end # type IoK8sApiAuthenticationV1TokenRequestSpec

const _property_map_IoK8sApiAuthenticationV1TokenRequestSpec = Dict{Symbol,Symbol}(Symbol("audiences")=>Symbol("audiences"), Symbol("boundObjectRef")=>Symbol("boundObjectRef"), Symbol("expirationSeconds")=>Symbol("expirationSeconds"))
const _property_types_IoK8sApiAuthenticationV1TokenRequestSpec = Dict{Symbol,String}(Symbol("audiences")=>"Vector{String}", Symbol("boundObjectRef")=>"IoK8sApiAuthenticationV1BoundObjectReference", Symbol("expirationSeconds")=>"Int64")
Base.propertynames(::Type{ IoK8sApiAuthenticationV1TokenRequestSpec }) = collect(keys(_property_map_IoK8sApiAuthenticationV1TokenRequestSpec))
Swagger.property_type(::Type{ IoK8sApiAuthenticationV1TokenRequestSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthenticationV1TokenRequestSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthenticationV1TokenRequestSpec }, property_name::Symbol) =  _property_map_IoK8sApiAuthenticationV1TokenRequestSpec[property_name]

function check_required(o::IoK8sApiAuthenticationV1TokenRequestSpec)
    (getproperty(o, Symbol("audiences")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAuthenticationV1TokenRequestSpec }, name::Symbol, val)
end


text:
#
# Protein - Water
#

using PDBTools
using MDDF
using Plots

# Here we use the PDBTools package to read the pdb file (from http://github.com/m3g/PDBTools)
atoms = PDBTools.readPDB("./structure.pdb")

# The solute is a single protein molecule (infinte dilution case). In this case,
# use the option nmols=1
solute_indexes = PDBTools.select(atoms,"protein")
solute = MDDF.Selection( solute_indexes, nmols=1 )

# The solvent is TMAO, which has 14 atoms. Use the natomspermol to indicate how many
# atoms each molecule has, such that there is no ambiguity on how to split the coordinates
# of the selection into individual molecules.
solvent_indexes = PDBTools.select(atoms,"resname TIP3")
solvent = MDDF.Selection( solvent_indexes, natomspermol=3 )

# Initialize trajectroy data structure and open input stream
trajectory = MDDF.Trajectory("./trajectory.dcd",solute,solvent)

# Input options for the calcualtion
options = MDDF.Options(binstep=0.2,n_random_samples=10)

# Run MDDF calculation, and get the resutls in the R structure
@time R = MDDF.mddf_linkedcells_parallel(trajectory,options)

plot(layout=(4,1))

x = [0,10]
y = [1, 1]

sp=1
plot!(ylabel="MDDF or RDF",subplot=sp)
plot!(R.d,R.mddf,subplot=sp,label="mddf")
plot!(R.d,R.rdf,subplot=sp,label="rdf")
plot!(x,y)
plot!(legend=:topright,subplot=sp)

sp=2
plot!(ylabel="KB",subplot=sp)
plot!(R.d,R.kb,subplot=sp,label="mddf")
plot!(R.d,R.kb_rdf,subplot=sp,label="rdf")
plot!(legend=:topright,subplot=sp)

sp=3
plot!(ylabel="atom contrib",subplot=sp)
for i in 1:solvent.natomspermol
  plot!(R.d,R.solvent_atom[:,i],subplot=sp,label="",linewidth=2)
end
plot!(legend=:topright,subplot=sp)

sp=4
plot!(ylabel="count",subplot=sp)
plot!(R.d,R.md_count,subplot=sp,label="R md")
plot!(R.d,R.md_count_random,subplot=sp,label="R rand")
plot!(legend=:topleft,subplot=sp)

plot!(size=(600,800))
savefig("./example5.pdf")
savefig("./plots.pdf")


text:
using Clang
using Test

@testset "macro" begin
    ctx = DefaultContext()
    parse_header!(ctx, joinpath(@__DIR__, "c", "enum.h"))
    trans_unit = ctx.trans_units[1]
    root_cursor = getcursor(trans_unit)
    cursors = children(root_cursor)
    wrap!(ctx, cursors[1])
    io = IOBuffer()
    print_buffer(io, dump_to_buffer(ctx.common_buffer))
    @test String(take!(io)) == "\n@cenum(_foo{Int32},\n    BAR = 0,\n    BAZ = -9,\n)\n"
end


text:
function build_A_sp(kx,U,E,Dy,Dy2,y,parms)
    U1 = U
    h1 = E.+parms.H

    # Define Differentiation Matrices for y: FD second order

    Ny  = length(U1)-1
    g1  = parms.g1

    #  Define Matrices to set up EVP

    dU1 = Diagonal(Dy*U1);
    dh1 = Diagonal(Dy*h1);

    F    = Diagonal(parms.f0.+0*y);
    Id    = Matrix{Float64}(I, Ny+1, Ny+1);
    h1f  = Diagonal(h1);

    kx2  = kx^2;
    ikx  = 1im*kx;

    # Form 1L-RSW Matrix
    #   [u1, v1, h1]

    A = [[Diagonal(U1),      -F[:,2:Ny].+dU1[:,2:Ny],                g1*I],
        [-F[2:Ny,:]/kx2,         Diagonal(U1[2:Ny]),  -g1/kx2*Dy[2:Ny,:]],
        [Diagonal(h1),       Dy*h1f[:,2:Ny],                    Diagonal(U1)]];
    return sparse(A)
end


text:
##############################################################################
#
# REFERENCES: Forbes et al. "Statistical Distributions"
#
##############################################################################

immutable BetaPrime <: ContinuousUnivariateDistribution
    alpha::Float64
    beta::Float64
    function BetaPrime(a::Real, b::Real)
        if a > 0.0 && b > 0.0
            new(float64(a), float64(b))
        else
            error("Both alpha and beta must be positive")
        end
    end
end

BetaPrime() = BetaPrime(1.0, 1.0)

cdf(d::BetaPrime, q::Real) = inc_beta(q / 1.0 + q, d.alpha, d.beta)

insupport(d::BetaPrime, x::Number) = isreal(x) && x > 0.0 ? true : false

function mean(d::BetaPrime)
    if d.beta > 1.0
        return d.alpha / (d.beta + 1.0)
    else
        error("mean not defined when beta <= 1")
    end
end

function modes(d::BetaPrime)
    if d.alpha > 1.0
        return [(d.alpha - 1.0) / (d.beta + 1.0)]
    else
        return [0.0]
    end
end

function pdf(d::BetaPrime, x::Real)
    a, b = d.alpha, d.beta
    return (x^(a - 1.0) * (1.0 + x)^(-(a + b))) / beta(a, b)
end

rand(d::BetaPrime) = 1.0 / rand(Beta(d.alpha, d.beta))

function skewness(d::BetaPrime)
    a, b = d.alpha, d.beta
    if b > 3.0
        return (2.0 * (2.0 * a + b - 1)) / (b - 3.0) *
               sqrt((b - 2.0) / (a * (a + b - 1.0)))
    else
        error("skewness not defined when beta <= 3")
    end
end

function var(d::BetaPrime)
    a, b = d.alpha, d.beta
    if b > 2.0
        return (a * (a + b - 1.0)) / ((b - 2.0) * (b - 1.0)^2)
    else
        error("var not defined when beta <= 2")
    end
end


text:
module NameResolution
using PrettyPrint
export Scope, LocalVar, GlobalVar
export is_global!, is_local!, enter!, require!, abs_interp_on_scopes, VarMap, run_analyzer
export child_analyzer!, new_analyzer, top_analyzer, new_scope, is_top_analyzer

include("Variable.jl")
VarMap = Dict{Symbol, LocalVar}
include("Analyzer.jl")

function is_global!(analyzer :: Analyzer, sym :: Symbol)
    push!(analyzer.globals, sym) ;nothing
end
function is_global!(:: Nothing, sym :: Symbol) end

function is_local!(analyzer :: Analyzer, sym :: Symbol)
    push!(analyzer.locals, sym) ;nothing
end
function is_local!(:: Nothing, sym :: Symbol) end

function enter!(analyzer :: Analyzer, sym::Symbol)
    if haskey(analyzer.entered, sym)
        analyzer.entered[sym] = true
    else
        analyzer.entered[sym] = false # has mutating sites
    end ;nothing
end
function enter!(:: Nothing, sym::Symbol) end

function require!(analyzer :: Analyzer, sym::Symbol)
    push!(analyzer.required, sym) ;nothing
end
function require!(:: Nothing, sym::Symbol) end

function request_freevar!(ana::Analyzer, var :: LocalVar)

    scope = ana.solved
    sym = var.sym
    bound = get(scope.bounds, var.sym, nothing)
    if bound === nothing
        getter = if ana.is_physical_scope
            var.is_shared.x = true
            get!
        else
            get
        end
        bound = getter(scope.freevars, sym) do
            # if 'var' is in free variables,
            request_freevar!(ana.parent, var)
            var
        end
    end
    @assert bound === var
    nothing
end

function abs_interp_on_scopes(analyzer::Analyzer, inherited::D) where {
    D <: AbstractDict{Symbol, LocalVar}
}
    scope    = analyzer.solved
    freevars = scope.freevars
    bounds   = scope.bounds

    entered  = analyzer.entered
    required = analyzer.required
    globals  = analyzer.globals
    locals   = analyzer.locals
    is_physical_scope = analyzer.is_physical_scope

    both_local_and_global = intersect(globals, locals)
    if !isempty(both_local_and_global)
        vars = join([string(x) for x in both_local_and_global], ", ")
        error("syntax: variable \"$vars\" declared both local and global")
    end

    for (sym, assign_twice) in entered
        if sym in globals
            @goto when_marked_global
        elseif sym in locals
            @goto when_bound
        elseif haskey(inherited, sym)
            @goto when_avaiable_outside
        else
            @goto when_bound
        end
        @label when_marked_global
            continue
        @label when_avaiable_outside
            var = inherited[sym]
            var.is_mutable.x = true
            request_freevar!(analyzer, var)
            continue
        @label when_bound
            var = get(bounds, sym, nothing)
            if var === nothing
                bounds[sym] = LocalVar(Ref(assign_twice), Ref(false), sym)
            else
                bounds[sym].is_mutable.x = true
            end
        continue # for good-looking
    end

    for sym in required
        if sym in globals
            continue
        elseif sym in locals
            var = get(bounds, sym, nothing)
            if var === nothing
                # uninitialized variables
                # TODO: add facilities to analyse the uninitialized variables
                bounds[sym] = LocalVar(Ref(false), Ref(false), sym)
            end
            continue
        elseif haskey(inherited, sym)
            # case for free variable
        elseif sym in globals && haskey(entered, sym)
            error("Writing global variable $sym without specifying the scope explicitly.\n"*
                  "Try to add 'global $sym' in that scope.")
        else
            continue
        end
        # is free variable
        var = inherited[sym]
        request_freevar!(analyzer, var)
    end
    inherited = VarMap(inherited..., bounds...)
    for child in analyzer.children
        abs_interp_on_scopes(child, inherited)
    end
end

run_analyzer(ana :: Analyzer) = abs_interp_on_scopes(ana, VarMap())
include("Pretty.jl")

end # module


text:
const dttmfmt = dateformat"YYYY.mm.dd-HH.MM.SS.sss"
const dtfmt = dateformat"YYYY.mm.dd"

struct OptionsTest
    d::Date
    dt::DateTime
end

JSON2.@format OptionsTest dateformat=dtfmt datetimeformat=dttmfmt

@testset "Read/write options" begin
    d = Date(1, 2, 3)
    dt = DateTime(1, 2, 3, 4, 5, 6, 7)

    @testset "Read" begin
        dstr = repr(Dates.format(d, ISODateFormat))
        dtstr = repr(Dates.format(dt, ISODateTimeFormat))

        @test JSON2.read(dstr, Date) == d
        @test JSON2.read(dtstr, DateTime) == dt

        @test JSON2.read(repr(Dates.format(d, dtfmt)), Date; dateformat=dtfmt) == d
        @test JSON2.read(repr(Dates.format(dt, dttmfmt)), DateTime; dateformat=dttmfmt) == dt
    end

    @testset "Write" begin
        @test JSON2.write(d) == repr(Dates.format(d, ISODateFormat))
        @test JSON2.write(dt) == repr(Dates.format(dt, ISODateTimeFormat))

        @test JSON2.write(d; dateformat=dtfmt) == repr(Dates.format(d, dtfmt))
        @test JSON2.write(dt; dateformat=dttmfmt) == repr(Dates.format(dt, dttmfmt))
    end

    @testset "Structs" begin
        opts = OptionsTest(d, dt)

        s = JSON2.write(opts; datetimeformat=dttmfmt, dateformat=dtfmt)
        @test JSON2.read(s, OptionsTest) == opts

        s = JSON2.write(opts; dateformat=Dates.ISODateFormat, datetimeformat=Dates.ISODateTimeFormat)
        @test JSON2.read(s, OptionsTest; dateformat=Dates.ISODateFormat, datetimeformat=Dates.ISODateTimeFormat) == opts
    end
end


text:
using Test
using MPI

MPI.Init()

parent_comm = MPI.Comm_get_parent()
world_comm = MPI.Intercomm_merge(parent_comm, true)
rank = MPI.Comm_rank(world_comm)

@test rank != 0

size = MPI.Comm_size(world_comm)
rank = MPI.Comm_rank(world_comm)
root = 0

val = rank == root ? size : nothing
@test MPI.Reduce(1, +, root, world_comm) == val

MPI.Comm_free(world_comm)
MPI.Comm_free(parent_comm)
MPI.Finalize()


text:
module Preprocessing
    using DataFrames
    using Statistics
    using SingularSpectrumAnalysis
    using ..ModelSelection

    export input, PREPROCESSING_EXTRAKEY

    include("const.jl")
    include("strings.jl")
    include("utils.jl")
    include("core.jl")
end


text:
"""
    preview_demos(demo_path; kwargs...)

Generate a docs preview for a single demo card.

# Return values

* `index_file`: absolute path to the index file of the demo page. You can open this in your browser
  and see the generated demos.

# Arguments

* `demo_path`: path to your demo file or folder. It can be path to the demo file, the folder of
  multiple scripts. If standard demo page folder structure is detected, use it, and otherwise will
  create a preview version of it.

# Parameters

* `theme`: the card theme you want to use in the preview. By default, it infers from your page
  config file. To not generate index page, you could pass `nothing` to it. See also
  [`cardtheme`](@ref DemoCards.cardtheme) for theme choices.
* `assets = String[]`: this is passed to `Documenter.HTML`
* `edit_branch = "master"`: same to that in `makedemos`
* `credit = true`: same to that in `makedemos`
* `throw_error = false`: same to that in `makedemos`
* `require_html = true`: if it needs to trigger the Documenter workflow and generate HTML preview.
  If set to `false`, the return value is then the path to the generated `index.md` file. This is an
  experimental keyword and should not be considered stable.
* `clean = true`: whether you need to first clean up the existing sandbox building dir.

"""
function preview_demos(demo_path::String;
                      theme = missing,
                      edit_branch = "master",
                      assets = String[],
                      credit = true,
                      throw_error = false,
                      require_html = true,
                      clean = true)
    # hard code these args in a sandbox environment -- there's no need for customization
    build = "build"
    src = "src"

    demo_path = abspath(rstrip(demo_path, ['/', '\\']))
    ispath(demo_path) || throw(ArgumentError("demo path does not exists: $demo_path"))

    build_dir = preview_build_dir()

    if clean
        # Ref: https://discourse.julialang.org/t/find-what-has-locked-held-a-file/23278/2
        Base.Sys.iswindows() && GC.gc()
        rm(build_dir; force=true, recursive=true)
        mkpath(build_dir)
    end

    page_dir = generate_or_copy_pagedir(demo_path, build_dir)
    copy_assets_and_configs(page_dir, build_dir)

    cd(build_dir) do
        page = @suppress_err DemoPage(page_dir)

        if theme === missing
            theme = page.theme
        else
            if page.theme != theme
                @info "Detect themes in both page config file and keyword `theme`, use \"$(string(theme))\"" config=string(page.theme) keyword=string(theme)
            end
        end

        if isnothing(theme)
            card_templates = nothing
        else
            card_templates, card_theme = cardtheme(theme; root = build_dir)
            push!(assets, card_theme)
        end

        page_index = joinpath(src, "index.md")
        mkpath(src)
        open(page_index, "w") do io
            println(io, "# Index page")
        end

        demos, demos_cb = makedemos(
            basename(page_dir), card_templates;
            root = build_dir,
            src = src,
            build = build,
            edit_branch = edit_branch,
            credit = credit,
            throw_error = throw_error,
        )

        # In cases that addtional Documenter pipeline is not needed
        # This is mostly for test usage right now and might be changed if there is better solution
        require_html || return abspath(build_dir, src, basename(page_dir))

        format = Documenter.HTML(
            edit_link = edit_branch,
            prettyurls = get(ENV, "CI", nothing) == "true",
            assets = assets
        )
        makedocs(
            root=build_dir,
            source=src,
            format=format,
            pages=[
                "Index" => relpath(page_index, src),
                demos
            ],
            sitename="DemoCards Preview"
        )

        demos_cb()

        return abspath(build, "index.html")
    end
end

# A common workflow when previewing demos is to directly refresh the browser this requires us to use
# the same build dir, otherwise, user might be confused why the contents are not updated.
const _preview_build_dir = String[]
function preview_build_dir()
    if isempty(_preview_build_dir)
        build_dir = mktempdir()
        push!(_preview_build_dir, build_dir)
    end
    return _preview_build_dir[1]
end

"""
    generate_or_copy_pagedir(src_path, build_dir)

copy the page folder structure from `src_path` to `build_dir/\$(basename(src_path)`. If `src_path`
does not live in a standard demo page folder structure, generate a preview version.
"""
function generate_or_copy_pagedir(src_path, build_dir)
    ispath(src_path) || throw(ArgumentError("$src_path does not exists"))

    if is_democard(src_path)
        card_path = src_path
        page_dir = infer_pagedir(card_path)
        if isnothing(page_dir)
            page_dir = abspath(build_dir, "preview_page")
            dst_card_path = abspath(page_dir, "preview_section", basename(card_path))
        else
            dst_card_path = abspath(build_dir, relpath(card_path, dirname(page_dir)))
        end
        mkpath(dirname(dst_card_path))
        cp(card_path, dst_card_path; force=true)
    elseif is_demosection(src_path)
        sec_dir = src_path
        page_dir = infer_pagedir(sec_dir)
        if isnothing(page_dir)
            page_dir = abspath(build_dir, "preview_page")
            dst_sec_dir = abspath(page_dir, basename(sec_dir))
        else
            dst_sec_dir = abspath(build_dir, relpath(sec_dir, dirname(page_dir)))
        end
        mkpath(dirname(dst_sec_dir))
        cp(sec_dir, dst_sec_dir; force=true)
    elseif is_demopage(src_path)
        page_dir = src_path
        cp(page_dir, abspath(build_dir, basename(page_dir)); force=true)
    else
        throw(ArgumentError("failed to parse demo page structure from path: $src_path. There might be some invalid demo files."))
    end

    return page_dir
end

"""
    copy_assets_and_configs(src_page_dir, dst_build_dir=pwd())

copy only assets, configs and templates from `src_page_dir` to `dst_build_dir`. The folder structure
is preserved under `dst_build_dir/\$(basename(src_page_dir))`

`order` in config files are modified accordingly.
"""
function copy_assets_and_configs(src_page_dir, dst_build_dir=pwd())
    (isabspath(src_page_dir) && isdir(src_page_dir)) || throw(ArgumentError("src_page_dir is expected to be absolute folder path: $src_page_dir"))

    for (root, dirs, files) in walkdir(src_page_dir)
        assets_dirs = intersect(ignored_dirnames, dirs)
        if !isempty(assets_dirs)
            for assets in assets_dirs
                src_assets = abspath(root, assets)
                dst_assets_path = abspath(dst_build_dir, relpath(src_assets, dirname(src_page_dir)))
                dst_sec_dir = dirname(dst_assets_path)

                # only copy assets when there are section or card entries
                if isdir(dst_sec_dir)
                    entry_names = filter(readdir(dst_sec_dir)) do x
                        x = joinpath(root, x)
                        is_democard(x) || is_demosection(x)
                    end
                    if !isempty(entry_names)
                        ispath(dst_assets_path) || cp(src_assets, dst_assets_path; force=true)
                    end
                end
            end
        end

        if template_filename in files
            src_template_path = abspath(root, template_filename)
            dst_template_path = abspath(dst_build_dir, relpath(src_template_path, dirname(src_page_dir)))
            mkpath(dirname(dst_template_path))
            ispath(dst_template_path) || cp(src_template_path, dst_template_path)
        end
    end

    # modify and copy config.json after the folder structure is already set up
    for (root, dirs, files) in walkdir(src_page_dir)
        if config_filename in files
            src_config_path = abspath(root, config_filename)
            dst_config_path = abspath(dst_build_dir, relpath(src_config_path, dirname(src_page_dir)))
            mkpath(dirname(dst_config_path))

            config = JSON.parsefile(src_config_path)
            config_dir = dirname(dst_config_path)
            if haskey(config, "order")
                order =  [x for x in config["order"] if x in readdir(config_dir)]
                if isempty(order)
                    delete!(config, "order")
                else
                    config["order"] = order
                end
            end

            if !isempty(config)
                # Ref: https://discourse.julialang.org/t/find-what-has-locked-held-a-file/23278/2
                Base.Sys.iswindows() && GC.gc()
                open(dst_config_path, "w") do io
                    JSON.print(io, config)
                end
            end
        end
    end
end


text:
using Images

fpart(x) = mod(x, one(x))
rfpart(x) = one(x) - fpart(x)

function drawline!(img::Matrix{Gray{N0f8}}, x0::Integer, y0::Integer, x1::Integer, y1::Integer)
    steep = abs(y1 - y0) > abs(x1 - x0)

    if steep
        x0, y0 = y0, x0
        x1, y1 = y1, x1
    end
    if x0 > x1
        x0, x1 = x1, x0
        y0, y1 = y1, y0
    end

    dx = x1 - x0
    dy = y1 - y0
    grad = dy / dx

    if iszero(dx)
        grad = oftype(grad, 1.0)
    end

    # handle first endpoint
    xend = round(Int, x0)
    yend = y0 + grad * (xend - x0)
    xgap = rfpart(x0 + 0.5)
    xpxl1 = xend
    ypxl1 = floor(Int, yend)

    if steep
        img[ypxl1,   xpxl1] = rfpart(yend) * xgap
        img[ypxl1+1, xpxl1] =  fpart(yend) * xgap
    else
        img[xpxl1, ypxl1  ] = rfpart(yend) * xgap
        img[xpxl1, ypxl1+1] =  fpart(yend) * xgap
    end
    intery = yend + grad # first y-intersection for the main loop

    # handle second endpoint
    xend = round(Int, x1)
    yend = y1 + grad * (xend - x1)
    xgap = fpart(x1 + 0.5)
    xpxl2 = xend
    ypxl2 = floor(Int, yend)
    if steep
        img[ypxl2,   xpxl2] = rfpart(yend) * xgap
        img[ypxl2+1, xpxl2] =  fpart(yend) * xgap
    else
        img[xpxl2, ypxl2  ] = rfpart(yend) * xgap
        img[xpxl2, ypxl2+1] =  fpart(yend) * xgap
    end

    # main loop
    if steep
        for x in xpxl1+1:xpxl2-1
            img[floor(Int, intery),   x] = rfpart(intery)
            img[floor(Int, intery)+1, x] =  fpart(intery)
            intery += grad
        end
    else
        for x in xpxl1+1:xpxl2-1
            img[x, floor(Int, intery)  ] = rfpart(intery)
            img[x, floor(Int, intery)+1] =  fpart(intery)
            intery += grad
        end
    end

    return img
end

img = fill(Gray(1.0N0f8), 250, 250);
drawline!(img, 8, 8, 192, 154)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module metrics
export eval_metrics

using CSV: File
using DataFrames: DataFrame
using Munkres: munkres

function eval_metrics(biclusters::Vector, input_path::String, ground_truth::Vector)

    ground_truth = deepcopy(ground_truth)

    for bclr in ground_truth
        bclr["cols"] .+= 1
        bclr["rows"] .+= 1
    end

    dataset = DataFrame(File(input_path))
    nrows = size(dataset, 1)
    ncols = size(dataset, 2) - 1 # omit column with g0, g1, ...

    relevance = prelic_relevance(biclusters, ground_truth)
    recovery = prelic_recovery(biclusters, ground_truth)
    ce = clustering_error(biclusters, ground_truth, nrows, ncols)

    return relevance, recovery, ce
end

function prelic_relevance(predicted_biclusters, reference_biclusters)::Float64
    """
    Based on: https://github.com/padilha/biclustlib/blob/master/biclustlib/evaluation/prelic.py
    """
    col_score = _match_score(predicted_biclusters, reference_biclusters, "cols")
    row_score = _match_score(predicted_biclusters, reference_biclusters, "rows")

    return sqrt(col_score * row_score)
end

function prelic_recovery(predicted_biclusters, reference_biclusters)::Float64
    return prelic_relevance(reference_biclusters, predicted_biclusters)
end

function _match_score(predicted_biclusters, reference_biclusters, attr)::Float64
    isempty(predicted_biclusters) && isempty(reference_biclusters) && return 1
    isempty(predicted_biclusters) || isempty(reference_biclusters) && return 0

    return sum([
        maximum([
            length(intersect(bp[attr], br[attr])) / length(union(bp[attr], br[attr]))
            for br in reference_biclusters
        ]) for bp in predicted_biclusters
    ]) / length(predicted_biclusters)
end

function clustering_error(predicted_biclusters, reference_biclusters, nrows, ncols)::Float64
    """
    Based on: https://github.com/padilha/biclustlib/blob/master/biclustlib/evaluation/subspace.py
    """
    isempty(predicted_biclusters) && isempty(reference_biclusters) && return 1
    isempty(predicted_biclusters) || isempty(reference_biclusters) && return 0

    union_size =
        _calc_size(predicted_biclusters, reference_biclusters, nrows, nrows, "union")
    dmax = _calc_dmax(predicted_biclusters, reference_biclusters)

    return dmax / union_size
end

function _calc_size(predicted_biclusters, reference_biclusters, nrows, ncols, operation)::Float64
    pred_count = _count_biclustering(predicted_biclusters, nrows, ncols)
    true_count = _count_biclustering(reference_biclusters, nrows, ncols)

    if operation == "union"
        return sum(broadcast(max, pred_count, true_count))
    elseif operation == "intersection"
        return sum(broadcast(min, pred_count, true_count))
    end

    throw(ArgumentError("Incorrect argument: $operation"))
end

function _count_biclustering(biclusters, nrows, ncols)
    count = zeros(Int, nrows, ncols)

    for b in biclusters
        for (x, y) in Base.product(b["rows"], b["cols"])
            count[x, y] += 1
        end
    end
    return count
end

function _calc_dmax(predicted_biclusters, reference_biclusters)::Float64
    pred_sets = _bic2sets(predicted_biclusters)
    true_sets = _bic2sets(reference_biclusters)

    cost_matrix = [
        maxintfloat() - length(intersect(b, g))
        for g in true_sets, b in pred_sets
    ]
    indices = munkres(cost_matrix)

    return sum([maxintfloat() - cost_matrix[x, y] for (x, y) in enumerate(indices)])
end

function _bic2sets(biclusters)
    return [Set(Base.product(b["rows"], b["cols"])) for b in biclusters]
end

end


text:
@safetestset "Binomial Model" begin
    using DifferentialEvolutionMCMC, Test, Random, Turing, Parameters, Distributions
    Random.seed!(29542)
    N = 10
    k = rand(Binomial(N, .5))
    data = (N = N,k = k)

    prior_loglike(θ) = logpdf(Beta(1, 1), θ)

    sample_prior() = rand(Beta(1, 1))

    bounds = ((0,1),)
    names = (:θ,)

    function loglike(data, θ)
        return logpdf(Binomial(data.N, θ), data.k)
    end

    model = DEModel(;
        sample_prior,
        prior_loglike,
        loglike,
        data,
        names
    )

    de = DE(;sample_prior, bounds, burnin=1500, Np=3)

    n_iter = 3000
    chains = sample(model, de, n_iter)
    μθ = describe(chains)[1][:,:mean][1]
    σθ = describe(chains)[1][:,:std][1]
    rhat = describe(chains)[1][:,:rhat][1]
    solution = Beta(k + 1, N - k + 1)
    @test μθ ≈ mean(solution) rtol = .02
    @test σθ ≈ std(solution) rtol = .02
    @test rhat ≈ 1.0 atol = .01
end