text:
"""
    AbstractRiemannianDiffBackend

An abstract type for diff backends. See [`RiemannianONBDiffBackend`](@ref) for
an example.
"""
abstract type AbstractRiemannianDiffBackend end

@doc raw"""
    differential(M::Manifold, f, t::Real, backend::AbstractDiffBackend = rdifferential_backend())

Compute the Riemannian differential of a curve $f: ℝ\to M$ on a manifold `M`
represented by function `f` at time `t` using the given backend.
It is calculated as the tangent vector equal to $\mathrm{d}f_t(t)[1]$.
"""
differential(::Manifold, ::Any, ::Real, ::AbstractRiemannianDiffBackend)

@doc raw"""
    gradient(M::Manifold, f, p, backend::AbstractRiemannianDiffBackend = rgradient_backend())

Compute the Riemannian gradient $∇f(p)$ of a real field on manifold `M` represented by
function `f` at point `p` using the given backend.
"""
gradient(::Manifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)

function differential!(M::Manifold, f::Any, X, t, backend::AbstractRiemannianDiffBackend)
    return copyto!(X, differential(M, f, t, backend))
end

function gradient!(M::Manifold, f, X, p, backend::AbstractRiemannianDiffBackend)
    return copyto!(X, gradient(M, f, p, backend))
end

differential(M::Manifold, f, p) = differential(M, f, p, rdifferential_backend())

differential!(M::Manifold, f, X, p) = differential!(M, f, X, p, rdifferential_backend())

gradient(M::Manifold, f, p) = gradient(M, f, p, rgradient_backend())

gradient!(M::Manifold, f, X, p) = gradient!(M, f, X, p, rgradient_backend())

"""
    RiemannianONBDiffBackend(
        diff_backend::AbstractDiffBackend
        retraction::AbstractRetractionMethod
        inverse_retraction::AbstractInverseRetractionMethod
        basis::Union{AbstractOrthonormalBasis,CachedBasis{<:AbstractOrthonormalBasis}},
    ) <: AbstractRiemannianDiffBackend

Riemannian differentiation based on differentiation in an [`AbstractOrthonormalBasis`](@ref)
`basis` using specified `retraction`, `inverse_retraction` and using backend `diff_backend`.
"""
struct RiemannianONBDiffBackend{
    TADBackend<:AbstractDiffBackend,
    TRetr<:AbstractRetractionMethod,
    TInvRetr<:AbstractInverseRetractionMethod,
    TBasis<:Union{
        AbstractOrthonormalBasis,
        CachedBasis{𝔽,<:AbstractOrthonormalBasis{𝔽}} where {𝔽},
    },
} <: AbstractRiemannianDiffBackend
    diff_backend::TADBackend
    retraction::TRetr
    inverse_retraction::TInvRetr
    basis::TBasis
end

function differential(M::Manifold, f, t::Real, backend::RiemannianONBDiffBackend)
    p = f(t)
    onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h
        return get_coordinates(
            M,
            p,
            inverse_retract(M, p, f(t + h), backend.inverse_retraction),
            backend.basis,
        )
    end
    return get_vector(M, p, onb_coords, backend.basis)
end

function differential!(M::Manifold, f, X, t::Real, backend::RiemannianONBDiffBackend)
    p = f(t)
    onb_coords = _derivative(zero(number_eltype(p)), backend.diff_backend) do h
        return get_coordinates(
            M,
            p,
            inverse_retract(M, p, f(t + h), backend.inverse_retraction),
            backend.basis,
        )
    end
    return get_vector!(M, X, p, onb_coords, backend.basis)
end

function gradient(M::Manifold, f, p, backend::RiemannianONBDiffBackend)
    X = get_coordinates(M, p, zero_tangent_vector(M, p), backend.basis)
    onb_coords = _gradient(X, backend.diff_backend) do Y
        return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction))
    end
    return get_vector(M, p, onb_coords, backend.basis)
end

function gradient!(M::Manifold, f, X, p, backend::RiemannianONBDiffBackend)
    X2 = get_coordinates(M, p, zero_tangent_vector(M, p), backend.basis)
    onb_coords = _gradient(X2, backend.diff_backend) do Y
        return f(retract(M, p, get_vector(M, p, Y, backend.basis), backend.retraction))
    end
    return get_vector!(M, X, p, onb_coords, backend.basis)
end

"""
    CurrentRiemannianDiffBackend(backend::AbstractRiemannianDiffBackend)

A mutable struct for storing the current Riemannian differentiation backend in global
constants [`Manifolds._current_rgradient_backend`](@ref) and
[`Manifolds._current_rdifferential_backend`](@ref).

# See also

[`AbstractRiemannianDiffBackend`](@ref), [`rdifferential_backend`](@ref),
[`rdifferential_backend!`](@ref)
"""
mutable struct CurrentRiemannianDiffBackend
    backend::AbstractRiemannianDiffBackend
end

"""
    _current_rgradient_backend

The instance of [`Manifolds.CurrentRiemannianDiffBackend`](@ref) that stores the
globally default differentiation backend for calculating gradients.

# See also

[`Manifolds.gradient(::Manifold, ::Any, ::Any, ::AbstractRiemannianDiffBackend)`](@ref)
"""
const _current_rgradient_backend = CurrentRiemannianDiffBackend(
    RiemannianONBDiffBackend(
        diff_backend(),
        ExponentialRetraction(),
        LogarithmicInverseRetraction(),
        DefaultOrthonormalBasis(),
    ),
)

"""
    _current_rdifferential_backend

The instance of [`Manifolds.CurrentRiemannianDiffBackend`](@ref) that stores the
globally default differentiation backend for calculating differentials.

# See also

[`Manifolds.differential`](@ref)
"""
const _current_rdifferential_backend = CurrentRiemannianDiffBackend(
    RiemannianONBDiffBackend(
        diff_backend(),
        ExponentialRetraction(),
        LogarithmicInverseRetraction(),
        DefaultOrthonormalBasis(),
    ),
)

"""
    rgradient_backend() -> AbstractRiemannianDiffBackend

Get the current differentiation backend for Riemannian gradients.
"""
rgradient_backend() = _current_rgradient_backend.backend

"""
    rgradient_backend!(backend::AbstractRiemannianDiffBackend)

Set current Riemannian gradient backend for differentiation to `backend`.
"""
function rgradient_backend!(backend::AbstractRiemannianDiffBackend)
    _current_rgradient_backend.backend = backend
    return backend
end

"""
    rdifferential_backend() -> AbstractRiemannianDiffBackend

Get the current differentiation backend for Riemannian differentials.
"""
rdifferential_backend() = _current_rdifferential_backend.backend

"""
    rdifferential_backend!(backend::AbstractRiemannianDiffBackend)

Set current Riemannian differential backend for differentiation to `backend`.
"""
function rdifferential_backend!(backend::AbstractRiemannianDiffBackend)
    _current_rdifferential_backend.backend = backend
    return backend
end

"""
    RiemannianProjectionGradientBackend(
        diff_backend::AbstractDiffBackend
    ) <: AbstractRiemannianDiffBackend

Riemannian differentiation based on differentiation in the ambient space and projection to
the given manifold. Differentiation in the ambient space is performed using
the backend `diff_backend`.

Only valid for manifolds that are embedded in a special way in the Euclidean space.
See [^Absil2008], Section 3.6.1 for details.

[^Absil2008]:
    > Absil, P. A., et al. Optimization Algorithms on Matrix Manifolds. 2008.
"""
struct RiemannianProjectionGradientBackend{TADBackend<:AbstractDiffBackend} <:
       AbstractRiemannianDiffBackend
    diff_backend::TADBackend
end

function gradient(M::Manifold, f, p, backend::RiemannianProjectionGradientBackend)
    amb_grad = _gradient(f, p, backend.diff_backend)
    return project(M, p, amb_grad)
end

function gradient!(M::Manifold, f, X, p, backend::RiemannianProjectionGradientBackend)
    amb_grad = embed(M, p, X)
    _gradient!(f, amb_grad, p, backend.diff_backend)
    return project!(M, X, p, amb_grad)
end


text:
module Impl

using Base: Dims
using Transducers
using Transducers: @return_if_reduced, next, complete

using ..Abstract

include("utils.jl")
include("vectors.jl")
include("arrays.jl")
include("reducibles.jl")

end  # module

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 51c330f6-7675-475e-b0eb-0fc8486fd094
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ 5c6ed844-279b-4576-a2fb-1c1d11fb62ee
using Arrowhead, LinearAlgebra

# ╔═╡ c3caaf95-bef6-446d-b197-23ad1a7da450
md"""
# Updating the SVD

In many applications which are based on the SVD, arrival of new data requires SVD of the new matrix. Instead of computing from scratch, existing SVD can be updated.

__Prerequisites__

The reader should be familiar with concepts of singular values and singular vectors, related perturbation theory, and algorithms.

__Competences__

The reader should be able to recognise applications where SVD updating can be sucessfully applied and apply it.
"""

# ╔═╡ e070e205-b8fb-4608-93c1-f5c5beb075c3
md"""
## Facts

For more details see
[M. Gu and S. C. Eisenstat, A Stable and Fast Algorithm for Updating the Singular Value Decomposition](http://www.cs.yale.edu/publications/techreports/tr966.pdf) and [M. Brand, Fast low-rank modifications of the thin singular value decomposition](http://www.sciencedirect.com/science/article/pii/S0024379505003812) and the references therein.

1. Let $A\in\mathbb{R}^{m\times n}$ with $m\geq n$ and $\mathop{\mathrm{rank}}(A)=n$, and  let $A=U\Sigma V^T$ be its SVD. Let $a\in\mathbb{R}^{n}$ be a vector, and let $\tilde A=\begin{bmatrix} A \\ a^T\end{bmatrix}$. Then

$$
\begin{bmatrix} A \\ a^T\end{bmatrix} =\begin{bmatrix} U & \\ & 1 \end{bmatrix}
\begin{bmatrix} \Sigma \\ a^TV \end{bmatrix}  V^T.$$

Let $\begin{bmatrix} \Sigma \\ a^T V \end{bmatrix} = \bar U \bar \Sigma \bar V^T$ be the SVD of the half-arrowhead matrix. _This SVD can be computed in $O(n^2)$ operations._ Then

$$
\begin{bmatrix} A \\ a^T\end{bmatrix} =
\begin{bmatrix} U & \\ & 1 \end{bmatrix} \bar U \bar\Sigma \bar V^T V^T \equiv
\tilde U \bar \Sigma \tilde V^T$$

is the SVD of $\tilde A$.

2. Direct computation of $\tilde U$ and $\tilde V$ requires $O(mn^2)$ and $O(n^3)$ operations. However, these multiplications can be performed using Fast Multipole Method. This is not (yet) implemented in Julia and is "not for the timid" (quote by Steven G. Johnson).

3. If $m<n$ and $\mathop{\mathrm{rank}}(A)=n$, then

$$
\begin{bmatrix} A \\ a^T\end{bmatrix} =\begin{bmatrix} U & \\ & 1 \end{bmatrix}
\begin{bmatrix} \Sigma & 0 \\ a^T V & \beta\end{bmatrix} \begin{bmatrix} V^T \\ v^T \end{bmatrix},$$

where $\beta=\sqrt{\|a\|_2^2-\|V^T a\|_2^2}$ and $v=(I-VV^T)a$. Notice that $V^Tv=0$ by construction. Let $\begin{bmatrix} \Sigma & 0 \\ a^T V &  \beta\end{bmatrix} = \bar U \bar \Sigma \bar V^T$ be the SVD of the half-arrowhead matrix. Then

$$
\begin{bmatrix} A \\ a^T\end{bmatrix} =
\begin{bmatrix} U & \\ & 1 \end{bmatrix} \bar U \bar\Sigma \bar V^T \begin{bmatrix} V^T \\ v^T \end{bmatrix}
\equiv \tilde U \bar \Sigma \tilde V^T$$

is the SVD of $\tilde A$.

4. Adding a column $a$ to $A$ is equivalent to adding a row $a^T$ to $A^T$.

5. If $\mathop{\mathrm{rank}}(A)<\min\{m,n\}$ or if we are using SVD approximation of rank $r$, and if we want to keep the rank of the approximation (this is the common case in practice), then the formulas in Fact 1 hold approximately. More precisely, the updated rank $r$ approximation is __not__ what we would get by computing the approximation of rank $r$ of the updated matrix, but is sufficient in many applications.
"""

# ╔═╡ 0ff9f1d6-1ebf-4a42-b51a-5a67f6192fe7
md"""

## Examples

### Adding row to a tall matrix

If $m\geq n$, adding row does not increase the size of $\Sigma$.
"""

# ╔═╡ 1f7f9f7a-df70-40a3-aa06-b5621b729a0c
function SVDaddrow(svdA::SVD,a::Vector)
    # Create the transposed half-arrowhead
    m,r,n=size(svdA.U,1),length(svdA.S),size(svdA.V,1)
    T=typeof(a[1])
    b=svdA.Vt*a
    if m>=n || r<m
        M=HalfArrow(svdA.S,b)
    else
        β=√(norm(a)^2-norm(b)^2)
        M=HalfArrow(svdA.S,[b;β])
    end
    # From Arrowhead package
    svdM,info=svd(M)
    # Return the updated SVD
    if m>=n || r<m
        return SVD([svdA.U zeros(T,m); zeros(T,1,r) one(T)]*svdM.V,
            svdM.S, adjoint(svdA.V*svdM.U))
    else
        # Need one more row of svdA.V - v is an orthogonal projection
        v=a-svdA.V*b
        normalize!(v)
        return SVD([svdA.U zeros(T,m); zeros(T,1,r) one(T)]*svdM.V,
            svdM.S, adjoint([svdA.V v]*svdM.U))
    end
end

# ╔═╡ 445cacab-3bb2-42f6-8d34-73c77ceecd0a
methods(SVD)

# ╔═╡ 2d0c21e5-dc68-4434-94ca-f433355f0e5c
begin
	import Random
	Random.seed!(421)
	A=rand(10,6)
	a=rand(6)
end

# ╔═╡ c034e218-09d3-4403-a413-98932865c481
S=svd(A)

# ╔═╡ 485e8222-2bb7-490e-9a51-675f0baef274
# Residual
norm(A*S.V-S.U*Diagonal(S.S))

# ╔═╡ 3dd2734d-db44-4226-9163-e6a55dcc7667
Sa=SVDaddrow(S,a)

# ╔═╡ 2fb74892-bade-47fe-bcaa-ef5bfdb19445
begin
	Aa=[A;transpose(a)]
	size(Aa),size(Sa.U),size(S.V)
end

# ╔═╡ c074206d-ed7d-44fb-a0ee-22e5cdd6e0f8
[svdvals(Aa) Sa.S]

# ╔═╡ 89990d20-4752-4558-a328-a4ae6a428906
# Residual and orthogonality
norm(Aa*Sa.V-Sa.U*Diagonal(Sa.S)),
norm(Sa.U'*Sa.U-I), norm(Sa.Vt*Sa.V-I)

# ╔═╡ bd35d5e4-b454-4a02-a026-ded3bdf76cdf
md"""
### Adding row to a flat matrix
"""

# ╔═╡ 8f31ef4a-a5da-4837-9320-ec9578c5a277
begin
	# Now flat matrix
	Random.seed!(421)
	A₁=rand(6,10)
	a₁=rand(10)
	S₁=svd(A₁);
end

# ╔═╡ e2ff1d44-7f3c-4100-80d3-df97579bc56c
A₁

# ╔═╡ 2b055add-b46f-4a2e-8ae1-89e59f12c21c
Aa₁=[A₁;transpose(a₁)]

# ╔═╡ e68d3450-da63-4b3f-9244-632a4eba5aff
Sa₁=SVDaddrow(S₁,a₁);

# ╔═╡ fac291b0-cf08-46d6-a233-3e1e8675c7e0
size(Aa₁),size(Sa₁.U),size(S₁.V)

# ╔═╡ 9f2d9dcc-ee9c-49e5-a839-3dadfed371e8
# Residual and orthogonality
norm(Aa₁*Sa₁.V-Sa₁.U*Diagonal(Sa₁.S)),
norm(Sa₁.U'*Sa₁.U-I), norm(Sa₁.Vt*Sa₁.V-I)

# ╔═╡ 9727df4c-52f5-4527-84e6-dfc02820d583
md"""
### Adding columns

This can be viewed as adding rows to the transposed matrix, an elegant one-liner in Julia.
"""

# ╔═╡ 3b34587c-fc0c-4677-91a7-4efb8a9abdba
function SVDaddcol(svdA::SVD,a::Vector)
    X=SVDaddrow(SVD(svdA.V,svdA.S,adjoint(svdA.U)),a)
    SVD(X.V,X.S,adjoint(X.U))
end

# ╔═╡ d8807efc-0cff-4aee-8bc5-4a5403c092eb
begin
	# Tall matrix
	Random.seed!(897)
	A₂=rand(10,6)
	a₂=rand(10)
	S₂=svd(A₂)
	Sa₂=SVDaddcol(S₂,a₂)
end

# ╔═╡ 8ce644c3-f71a-4904-bf15-d3f602f44c05
Aa₂=[A₂ a₂];

# ╔═╡ ca5fd757-5ab8-4826-8636-59a29a8c5843
# Residual and orthogonality
norm(Aa₂*Sa₂.V-Sa₂.U*Diagonal(Sa₂.S)),
norm(Sa₂.U'*Sa₂.U-I), norm(Sa₂.Vt*Sa₂.V-I)

# ╔═╡ 5f5e6293-17b5-405d-96c6-d661151796dc
begin
	# Flat matrix
	Random.seed!(332)
	A₃=rand(6,10)
	a₃=rand(6)
	S₃=svd(A₃)
	Sa₃=SVDaddcol(S₃,a₃)
end

# ╔═╡ 3c5d82e1-81f4-4cff-9c2a-6cc9bfbefd58
Aa₃=[A₃ a₃]

# ╔═╡ a2f5ab12-5c09-4f2b-9665-d2d3d44abd17
# Residual and orthogonality
norm(Aa₃*Sa₃.V-Sa₃.U*Diagonal(Sa₃.S)),
norm(Sa₃.U'*Sa₃.U-I), norm(Sa₃.Vt*Sa₃.V-I)

# ╔═╡ 6563cdc2-8577-4fac-997d-413edc56e2ac
begin
	# Square matrix
	A₄=rand(10,10)
	a₄=rand(10)
	S₄=svd(A₄)
end

# ╔═╡ 0f656a69-2243-4273-8600-83a12b694497
begin
	Sa₄=SVDaddrow(S₄,a₄)
	Aa₄=[A₄;transpose(a₄)]
	norm(Aa₄*Sa₄.V-Sa₄.U*Diagonal(Sa₄.S)),
	norm(Sa₄.U'*Sa₄.U-I), norm(Sa₄.Vt*Sa₄.V-I)
end

# ╔═╡ 35b4c90e-830f-487e-b683-988952aa12af
begin
	Sa₅=SVDaddcol(S₄,a₄)
	Aa₅=[A₄ a₄]
	norm(Aa₅*Sa₅.V-Sa₅.U*Diagonal(Sa₅.S)),
	 norm(Sa₅.U'*Sa₅.U-I), norm(Sa₅.Vt*Sa₅.V-I)
end

# ╔═╡ 22b4dea8-79a0-438a-a3e4-d7c725d47051
md"""
### Updating a low rank approximation

"""

# ╔═╡ 0a5298eb-ed33-443b-b9e4-c0a32b408605
begin
	# Adding row to a tall matrix
	A₆=rand(10,6)
	S₆=svd(A₆)
	a₆=rand(6)
	# Rank of the approximation
	r=4
end

# ╔═╡ 99b2a284-4a24-492d-b677-4b08430a4399
# Low-rank approximation
Sr=SVD(S₆.U[:,1:r], S₆.S[1:r],adjoint(S₆.V[:,1:r]));

# ╔═╡ 3677f364-e828-4a42-84b0-e86f0e3fff3b
begin
	# Eckart, Young, Mirsky
	Ar=Sr.U*Diagonal(Sr.S)*Sr.Vt
	Δ=Ar-A₆
	opnorm(Δ),svdvals(A₆)[5]
end

# ╔═╡ 32880612-cdd1-4e16-8a05-025eb35110f9
Sa₆=SVDaddrow(Sr,a₆);

# ╔═╡ 65c83ee2-36ad-4388-b9b3-05830ab3ee9a
Aa₆=[A₆; transpose(a₆)];

# ╔═╡ 13deeb53-9670-4351-85b3-a366690a94a8
svdvals(Aa₆),svdvals([Ar;transpose(a₆)]),Sa₆.S

# ╔═╡ a02b7022-b775-4dbd-bb83-0d6f88a47312
begin
	# Adding row to a flat matrix
	A₇=rand(6,10)
	S₇=svd(A₇)
	a₇=rand(10)
	# Rank of the approximation
	r₇=4
end

# ╔═╡ 2f22e978-2fb4-4ec4-80af-519384e98f6c
begin
	Sr₇=SVD(S₇.U[:,1:r₇], S₇.S[1:r₇],adjoint(S₇.V[:,1:r₇]))
	Sa₇=SVDaddrow(Sr₇,a₇);
end

# ╔═╡ 2beba2fd-cb70-4da2-a9f3-d4c39815ff18
begin
	Ar₇=Sr₇.U*Diagonal(Sr₇.S)*Sr₇.Vt
	svdvals(Sa₇),svdvals([Ar₇;transpose(a₇)]),Sa₇.S
end

# ╔═╡ Cell order:
# ╟─51c330f6-7675-475e-b0eb-0fc8486fd094
# ╟─c3caaf95-bef6-446d-b197-23ad1a7da450
# ╟─e070e205-b8fb-4608-93c1-f5c5beb075c3
# ╟─0ff9f1d6-1ebf-4a42-b51a-5a67f6192fe7
# ╠═5c6ed844-279b-4576-a2fb-1c1d11fb62ee
# ╠═1f7f9f7a-df70-40a3-aa06-b5621b729a0c
# ╠═445cacab-3bb2-42f6-8d34-73c77ceecd0a
# ╠═2d0c21e5-dc68-4434-94ca-f433355f0e5c
# ╠═c034e218-09d3-4403-a413-98932865c481
# ╠═485e8222-2bb7-490e-9a51-675f0baef274
# ╠═3dd2734d-db44-4226-9163-e6a55dcc7667
# ╠═2fb74892-bade-47fe-bcaa-ef5bfdb19445
# ╠═c074206d-ed7d-44fb-a0ee-22e5cdd6e0f8
# ╠═89990d20-4752-4558-a328-a4ae6a428906
# ╟─bd35d5e4-b454-4a02-a026-ded3bdf76cdf
# ╠═8f31ef4a-a5da-4837-9320-ec9578c5a277
# ╠═e2ff1d44-7f3c-4100-80d3-df97579bc56c
# ╠═2b055add-b46f-4a2e-8ae1-89e59f12c21c
# ╠═e68d3450-da63-4b3f-9244-632a4eba5aff
# ╠═fac291b0-cf08-46d6-a233-3e1e8675c7e0
# ╠═9f2d9dcc-ee9c-49e5-a839-3dadfed371e8
# ╟─9727df4c-52f5-4527-84e6-dfc02820d583
# ╠═3b34587c-fc0c-4677-91a7-4efb8a9abdba
# ╠═d8807efc-0cff-4aee-8bc5-4a5403c092eb
# ╠═8ce644c3-f71a-4904-bf15-d3f602f44c05
# ╠═ca5fd757-5ab8-4826-8636-59a29a8c5843
# ╠═5f5e6293-17b5-405d-96c6-d661151796dc
# ╠═3c5d82e1-81f4-4cff-9c2a-6cc9bfbefd58
# ╠═a2f5ab12-5c09-4f2b-9665-d2d3d44abd17
# ╠═6563cdc2-8577-4fac-997d-413edc56e2ac
# ╠═0f656a69-2243-4273-8600-83a12b694497
# ╠═35b4c90e-830f-487e-b683-988952aa12af
# ╟─22b4dea8-79a0-438a-a3e4-d7c725d47051
# ╠═0a5298eb-ed33-443b-b9e4-c0a32b408605
# ╠═99b2a284-4a24-492d-b677-4b08430a4399
# ╠═3677f364-e828-4a42-84b0-e86f0e3fff3b
# ╠═32880612-cdd1-4e16-8a05-025eb35110f9
# ╠═65c83ee2-36ad-4388-b9b3-05830ab3ee9a
# ╠═13deeb53-9670-4351-85b3-a366690a94a8
# ╠═a02b7022-b775-4dbd-bb83-0d6f88a47312
# ╠═2f22e978-2fb4-4ec4-80af-519384e98f6c
# ╠═2beba2fd-cb70-4da2-a9f3-d4c39815ff18

max_stars_repo_path: src/Module B - Eigenvalue and Singular Value Decompositions/L8 Updating the SVD.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:
function randDirichlet(α::T...) where T<:Real
    p = SVector((αi -> Gamma(αi,one(T))).(α))
    p ./ sum(p)
end
randDirichlet(α::Integer...) = randDirichlet(Float64.(α)...)
randDirichlet(n::Integer, α::Real) = randDirichlet((α for _ in 1:n)...)
randDirichlet(n::Integer, α::Integer) = randDirichlet(n, Float64(α))

function logpdfDirichlet(x, α::Real...)
    a, b = sum(u -> SVector(u,loggamma(u)), α)
    s = sum(((u,v) -> (u-one(u))log(v)).(α, x))
    s - b + lgamma(a)
end

function _logpdfDirichlet(x, α::Real...)
    a, b = sum(u -> SVector(u,lgamma(u)), α)
    s = sum(((u,v) -> (u-one(u))CUDAnative.log(v)).(α, x))
    s - b + CUDAnative.lgamma(a)
end


text:
"""
    AbstractConstraints{T}

Abstract type for all possible constructs.
"""
abstract type AbstractConstraints{T} end

"""
    ElementConstraints{T}

Abstract Constraints with effect on the element that are allow in the constructs.
"""
abstract type ElementConstraints{T} <: AbstractConstraints{T} end

"""
    NoConstraint{T}

If no constraint is given.
"""
struct NoConstraint{T} <: ElementConstraints{T}
    con :: T
end

"""
    ConstructConstraints{T}

 If a particular construct is allowed is checked on the construct level,
 using fitler that block all unwanted constructs.
 This implies that first all constructed are generated and afterwards filterd,
 which is not efficient for large space with a limited number of constraints.
`filterconstrain` than filters the removes the unwanted constructs
 For more information see [`filterconstrain`](@ref)
"""
abstract type ConstructConstraints{T} <: ElementConstraints{T} end

"""
    SingleConstructConstraints{T}

 Type used to as input for a single constraint:
 two types are currently implemented:
 `UnOrderedConstraint` and `OrderedConstraint`.
"""
abstract type SingleConstructConstraints{T} <: ConstructConstraints{T} end

"""
    UnOrderedConstraint{T<:Mod}

an unorderedConstraint is a group of modules that cannot co-occur in one construct.
These group is given as:  ``Array{T,1} where T <: Mod``
This constraint is unorder because the location of these models in the constructed
is not considered

For more information see [`filterconstrain`]@ref
"""
struct UnOrderedConstraint{T<:Mod} <: SingleConstructConstraints{T}
    combination::Array{T}
end

make_unorderconstraint(c::Array) = UnOrderedConstraint([Mod(m) for m in c])

"""
    OrderedConstraint{T<:Mod}

Order Constraints a group of modules that can not co-occur in one construct if
these modules are in a specific position.
This type has two inputs both vector of the same length
The first vector contains the indexes,
the second contains the forbidden module at the correspoding possiton
The contraint is order because the location of these modules
in the constructed is used in the
evaluation of the constraint.

The corresponding positions are ``Array{::Int}``, these are the forbidden indexes for the corresponding module.
These modules are given as ``Array{T} where T <: Mod``

For more information see [`filterconstraint`]@ref.
or read the documentation regareding the use of constraints.

"""
struct OrderedConstraint{T<:Mod} <: SingleConstructConstraints{T}
    pos::Array{Int}
    combination::Array{T}
end

make_orderconstraint(pos::Array{Int},c::Array) = OrderedConstraint(pos,[Mod(m) for m in c])

"""
   ComposedConstructConstraints{T}

Constraint that are applicable on a generated construct.
These can only be used in a filter that is not very efficient for large design
spaces.
"""
struct ComposedConstructConstraints{T <: SingleConstructConstraints} <: ConstructConstraints{T}
    constructcon :: Array{T}
end

Base.eltype(K::ComposedConstructConstraints{T}) where {T} = T

# concentration of multiple Constraints can be done using + sign

"""
    +(con1::SingleConstructConstraints,con2::SingleConstructConstraints)

Returns a `ComposedConstructConstraints` containing both Constraints `con1` and `con2`.
"""
Base.:+(con1::SingleConstructConstraints,con2::SingleConstructConstraints) = ComposedConstructConstraints([con1,con2])

"""
    Base.:+(con1::ComposedConstructConstraints{T} where T, con2::N where N <: SingleConstructConstraints)

Returns a `ComposedConstructConstraints` containing where `con2` is add to the `con1`
"""

Base.:+(con1::ComposedConstructConstraints{T} where T, con2::N where N <: SingleConstructConstraints) = ComposedConstructConstraints([con1.constructcon;con2])

"""
    Base.:+(con1::SingleConstructConstraints,con2::ComposedConstructConstraints)

Returns a `ComposedConstructConstraints` containing where `con1` is add to the `con2`
"""
Base.:+(con1::SingleConstructConstraints,con2::ComposedConstructConstraints) = +(con2,con1)

"""
    Base.:+(con1::ComposedConstructConstraints{T} where T, con2::ComposedConstructConstraints{N} where N)

Returns a `ComposedConstructConstraints` concatenation of `con1` and `con2`
"""
Base.:+(con1::ComposedConstructConstraints{T} where T, con2::ComposedConstructConstraints{N} where N) = ComposedConstructConstraints([con1.constructcon;con2.constructcon;])

##################################
# Code below has currently no application in the package, this will probably not change shortly
# The idea is that for some constraints, we can prevent that unwanted constructs are generated, instead of generational filtering afterwards.
# To allow this a special type was constructed and see how it good be implemented
# No need to check,
#################################

#=
"""
space contrains has effect on how the space is repesanted without real evaluation of the constructs
"""

abstract type Space_Constraints{T} <: ElementConstraints{T} end

"""
Single possition constrain. If allow = 1 it indicaces that only these construcst are allow.
if allow = 0is false thes models aren't allowed on this possition, for non position specific constrain pos = 0
"""

struct Compose_Space_Constrain{Ts <: Space_Constraints{T} where T} <: Space_Constraints{Ts}
    space_con::AbstractArray{Ts}
end

abstract type Single_Space_Constraints{T} <: Space_Constraints{T} end

struct Possition_Constrain{T} <: Single_Space_Constraints{T}
    pos::Int
    mod::T
end

# add two  space Constraints
Base.:+(con1::N  ,con2::N ) where N <: Single_Space_Constraints{T} where T = Compose_Space_Constrain([con1,con2])
# add to compose Constraintse
Base.:+(con1::Compose_Space_Constrain,con2::Single_Space_Constraints) = Compose_Space_Constrain([con1.space_con;con2])
Base.:+(con1::Single_Space_Constraints,con2::Compose_Space_Constrain) = Compose_Space_Constrain([con2.space_con;con1])

Base.:+(con1::Compose_Space_Constrain,con2::Compose_Space_Constrain)  = Compose_space_Constraints([con1.space_con;con2.space_con])

Base.:+(con1::Element_Constraints,con2::Element_Constraints...) = +(con1,+(con2...))
"""
Group all elementConstraints
"""

struct Multipe_constrain{T<:Any, Ts<:Space_Constraints, Tc<:Construct_Constraints}  <: Element_Constraints{T}
    space_con:: Ts
    construct_con:: Tc
    function Multipe_constrain(space_con::Space_Constraints{T}, construct_con::Construct_Constraints{V}) where {T, V}
        return new{promote_type(T, V),typeof(space_con),typeof(construct_con)}(space_con,construct_con)
    end
end

Base.:+(con1:: Space_Constraints, con2:: Construct_Constraints) = Multipe_constrain(con1,con2)
Base.:+(con1:: Construct_Constraints, con2:: Space_Constraints) = Multipe_constrain(con2,con1)

Base.:+(con1::Multipe_constrain,con2::Construct_Constraints) = Multipe_constrain(con1.space_con, +(con1.construct_con,con2))
Base.:+(con1::Multipe_constrain,con2::Space_Constraints) = Multipe_constrain(+(con1.space_con,con2),con1.construct_con)

Base.:+(con1::Construct_Constraints,con2::Multipe_constrain) = Multipe_constrain(con2,con1)
Base.:+(con1::Space_Constraints,con2::Multipe_constrain) =  Multipe_constrain(con2,con1)

Base.:+(con1::Multipe_constrain,con2::Multipe_constrain) =  Multipe_constrain(con1.space_con + con2.space_con, con1.construct_con + con2.construct_con)

# fucntion(con1::Multipe_constrain,con2::Space_Constraints)

"""
Group length and Constraints to one object
"""
abstract type Full_Constrain{T} <: AbstractConstraints{T} end

"""
Gives an specific lenth to given Constraints
"""

struct Len_Constrain{T} <: Full_Constrain{T}
    len::Int
    el_con::Element_Constraints{T}
end

struct Group_Len_Constrain{T} <: Full_Constrain{T}
    group::Tuple{T}
end

group_len_constrain(con1::Len_Constrain,con2::Len_Constrain) = Group_Len_Constrain(Tuple(con1,con2))
group_len_constrain(con1::Len_Constrain,con2::Len_Constrain...) = Group_Len_Constrain(Tuple(con1,Tuple(con2...)))
=#

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


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

# # Tips and tricks

# **Originally Contributed by**: Arpit Bhatia

# !!! tip
#     A good source of tips is the [Mosek Modeling Cookbook](https://docs.mosek.com/modeling-cookbook/mio.html).

# This tutorial collates some tips and tricks you can use when formulating
# mixed-integer programs. It uses the following packages:

using JuMP

# ## Boolean Operators on Binary Variables

# Binary variables can be used to construct logical operators. Here are some
# example.

# ### OR

# $$x_3 = x_1 \lor x_2$$

model = Model()
@variable(model, x[1:3], Bin)
@constraints(model, begin
    x[1] <= x[3]
    x[2] <= x[3]
end)

# ### And

# $$x_3 = x_1 \land x_2$$

model = Model()
@variable(model, x[1:3], Bin)
@constraints(model, begin
    x[3] <= x[1]
    x[3] <= x[2]
    x[3] >= x[1] + x[2] - 1
end)

# ### Not

# $$x_1 \neg x_2$$

model = Model()
@variable(model, x[1:2], Bin)
@constraint(model, x[1] == 1 - x[2])

# ### Implies

# $$x_1 \implies x_2$$

model = Model()
@variable(model, x[1:2], Bin)
@constraint(model, x[1] <= x[2])

# ## Big-M Disjunctive Constraints (OR)

# **Problem** Suppose that we have two constraints $a^\top x \leq b$ and
# $c^\top x \leq d$, and we want at least one to hold.

# **Trick** Introduce a "big-M" multiplied by a binary variable to relax one of
# the constraints.

# **Example** Either $x_1 \leq 1$ and/or $x_2 \leq 2$.

model = Model()
@variable(model, x[1:2])
@variable(model, y, Bin)
M = 100
@constraint(model, x[1] <= 1 + M * y)
@constraint(model, x[2] <= 2 + M * (1 - y))

# !!! warning
#     If `M` is too small, the solution may be suboptimal. If `M` is too big,
#     the solver may encounter numerical issues. Try to use domain knowledge to
#     choose an `M` that is just right. Gurobi has a [good documentation section](https://www.gurobi.com/documentation/9.1/refman/dealing_with_big_m_constra.html)
#     on this topic.

# ## Indicator Constraints ($\implies$)

# ### Problem

# Suppose we want to model that a certain linear inequality must be satisfied
# when some other event occurs, i.e., for a binary variable $z$, we want to
# model the implication:

# $$z = 1 \implies a^Tx \leq b$$

# ### Trick 1

# Some solvers have native support for indicator constraints.

# **Example** $x_1 + x_2 \leq 1$ if $z = 1$.

model = Model()
@variable(model, x[1:2])
@variable(model, z, Bin)
@constraint(model, z => {sum(x) <= 1})

# **Example** $x_1 + x_2 \leq 1$ if $z = 0$.

model = Model()
@variable(model, x[1:2])
@variable(model, z, Bin)
@constraint(model, !z => {sum(x) <= 1})

# ### Trick 2

# If the solver doesn't support indicator constraints, you an use the big-M
# trick.

# **Example** $x_1 + x_2 \leq 1$ if $z = 1$.

model = Model()
@variable(model, x[1:2])
@variable(model, z, Bin)
M = 100
@constraint(model, sum(x) <=  1 + M * (1 - z))

# **Example** $x_1 + x_2 \leq 1$ if $z = 0$.

model = Model()
@variable(model, x[1:2])
@variable(model, z, Bin)
M = 100
@constraint(model, sum(x) <=  1 + M * z)

# ## Semi-Continuous Variables

# A semi-continuous variable is a continuous variable between bounds $[l,u]$
# that also can assume the value zero. ie.
# $$x \in \{0\} \cup [l,u].$$

# **Example** $$x \in \{0\}\cup [1, 2]$$

model = Model()
@variable(model, x in MOI.Semicontinuous(1.0, 2.0))

# ## Semi-Integer Variables

# A semi-integer variable is a variable which asummes integer values between
# bounds $[l,u]$ and can also assume the value zero:
# $$x \in \{0\} \cup [l, u] \cap \mathbb{Z}.$$

model = Model()
@variable(model, x in MOI.Semiinteger(5.0, 10.0))

# ## Special Ordered Sets of Type I (SOS1)

# A Special Ordered Set of Type I is a set of variables, at most one of which
# can take a non-zero value, all others being at 0.
#
# They most frequently apply where a set of variables are actually binary
# variables. In other words, we have to choose at most one from a set of
# possibilities.

model = Model()
@variable(model, x[1:3], Bin)
@constraint(model, x in SOS1())

# You can optionally pass `SOS1` a weight vector like

@constraint(model, x in SOS1([0.2, 0.5, 0.3]))

# If the decision variables are related and have a physical ordering, then the
# weight vector, although not used directly in the constraint, can help the
# solver make a better decision in the solution process.

# ## [Special Ordered Sets of Type II (SOS2)](@id tip_sos2)

# A Special Ordered Set of type 2 is a set of non-negative variables, of which
# at most two can be non-zero, and if two are non-zero these must be consecutive
# in their ordering.

model = Model()
@variable(model, x[1:3])
@constraint(model, x in MOI.SOS2([3.0, 1.0, 2.0]))

# The ordering provided by the weight vector is more important in this case as
# the variables need to be consecutive according to the ordering.
# For example, in the above constraint, the possible pairs are:
# * Consecutive
#   * (`x[1]` and `x[3]`) as they correspond to 3 and 2 resp. and thus can be non-zero
#   * (`x[2]` and `x[3]`) as they correspond to 1 and 2 resp. and thus can be non-zero
# * Non-consecutive
#   * (`x[1]` and `x[2]`) as they correspond to 3 and 1 resp. and thus cannot be non-zero

# ## Piecewise linear approximations

# [SOSII constraints](@ref tip_sos2) are most often used to form piecewise
# linear approximations of a function.

# Given a set of points for `x`:
x̂ = -1:0.5:2
# and a set of corresponding points for `y`:
ŷ = x̂.^2
# the piecewise linear approximation is constructed by reprsenting `x` and `y`
# as convex combinations of `x̂` and `ŷ`.

N = length(x̂)
model = Model()
@variable(model, -1 <= x <= 2)
@variable(model, y)
@variable(model, 0 <= λ[1:N] <= 1)
@objective(model, Max, y)
@constraints(model, begin
    x == sum(x̂[i] * λ[i] for i = 1:N)
    y == sum(ŷ[i] * λ[i] for i = 1:N)
    sum(λ) == 1
    λ in SOS2()
end)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function a()
    x = 0
	repeat
        x = + x 1
        print(x)
    until
        >= x 10

end

max_stars_repo_name: tylerholmes/Interpreter

max_issues_repo_name: tylerholmes/Interpreter

max_forks_repo_name: tylerholmes/Interpreter


text:
mutable struct SpanningTreeLightGraphsPrimSolver <: SpanningTreeSolver
  graph::SimpleGraph
  weights_matrix::Matrix{Float64}

  function SpanningTreeLightGraphsPrimSolver()
    return new(Graph(0), zeros(0, 0))
  end
end

function build!(solver::SpanningTreeLightGraphsPrimSolver, graph::SimpleGraph)
  n = nv(graph)
  solver.graph = graph
  solver.weights_matrix = zeros(n, n)
end

function solve_linear(solver::SpanningTreeLightGraphsPrimSolver, reward::Dict{Tuple{Int, Int}, Float64})
  # Make up a rewards matrix by copying the input dictionary into the right data structure for the shortest path computations.
  fill!(solver.weights_matrix, 0.0)
  for (k, v) in reward
    solver.weights_matrix[k[1], k[2]] = v
  end
  solver.weights_matrix = maximum(solver.weights_matrix) .- solver.weights_matrix

  # Compute the maximum spanning tree.
  mst = prim_mst(solver.graph, solver.weights_matrix)

  # Transform edges into tuples.
  return _mst_solution_normalise(reward, [(src(e), dst(e)) for e in mst])
end

has_lp_formulation(::SpanningTreeLightGraphsPrimSolver) = false


text:
# import ThreeBodyProblem.Constants
# include("constants.jl")
"""
    Body(m, R, a, T, name, color)
    Body(m, R, a, T, name)
    Body(m, R, a, T)

A planet, moon, or other gravitationally significant object
m::Float64      # mass {kg}
R::Float64      # mean radius {km}
a::Float64      # mean semimajor axis about parent body {km}
T::Float64      # sidereal orbital period {s}
name::String    # name of body (e.g. "Earth")
color::Symbol   # color of the body for plotting
"""
struct Body
    m::Float64      # mass {kg}
    R::Float64      # mean radius {km}
    a::Float64      # mean semimajor axis about parent body {km}
    T::Float64      # sidereal orbital period {s}
    name::String    # name of body (e.g. "Earth")
    color::Symbol   # color of the body for plotting
    # parent::Body    # parent body (e.g. SUN for EARTH)
end
Body(m::Float64, R::Float64, a::Float64, T::Float64, name::String) = Body(m, R, a, T, name, :blue)
Body(m::Float64, R::Float64, a::Float64, T::Float64) = Body(m, R, a, T, "NewPlanet", :blue)

"""
    System(prim, sec, μ₁, μ₂, μ, d, R₁, R₂, T, RUNIT, VUNIT, TUNIT, name)
    System(prim::Body, sec::Body)

A Circular Restricted Three-Body Problem System defined by primary and secondary bodies (prim and sec)
prim::Body      # Primary body
sec::Body       # Secondary body
μ₁::Float64     # {km^3/s^2} gravitational parameter of primary body
μ₂::Float64     # {km^3/s^2} gravitational parameter of secondary body
μ::Float64      # {} mass parameter
d::Float64      # {km} average distance between two primaries
R₁::Float64     # {km} Radius of primary body
R₂::Float64     # {km} Radius of secondary Body
T::Float64      # {s} sidereal orbital period
RUNIT::Float64  # {km} distance normalizing parameter
VUNIT::Float64  # {km/s} velocity normalizing parameter
TUNIT::Float64  # {s} time normalizing parameter
name::String    # name of system (e.g. "Earth/Moon")
"""
struct System
    prim::Body      # Primary body
    sec::Body       # Secondary body
    μ₁::Float64     # {km^3/s^2} gravitational parameter of primary body
    μ₂::Float64     # {km^3/s^2} gravitational parameter of secondary body
    μ::Float64      # {} mass parameter
    d::Float64      # {km} average distance between two primaries
    R₁::Float64     # {km} Radius of primary Body
    R₂::Float64     # {km} Radius of Secondary Body
    T::Float64      # {s} sidereal orbital period
    RUNIT::Float64  # {km} distance normalizing parameter
    TUNIT::Float64  # {s} time normalizing parameter
    VUNIT::Float64  # {km/s} velocity normalizing parameter
    AUNIT::Float64  # {km^2/s} acceleration normalizing parameter
    name::String    # name of system (e.g. "Earth/Moon")

    # Added these since it is unwise to require unicode operators from the user
    mu1::Float64     # {km^3/s^2} gravitational parameter of primary body
    mu2::Float64     # {km^3/s^2} gravitational parameter of secondary body
    mu::Float64      # {} mass parameter
    R1::Float64     # {km} Radius of primary Body
    R2::Float64     # {km} Radius of Secondary Body
end
System(prim::Body, sec::Body) = System(prim, sec, prim.m*G, sec.m*G, sec.m/(prim.m+sec.m), sec.a, prim.R, sec.R, sec.T, sec.a, sec.T/2π, sec.a/(sec.T/2π), sec.a/(sec.T/2π)^2, string(prim.name,"/",sec.name),
                                                  prim.m*G, sec.m*G, sec.m/(prim.m+sec.m),        prim.R, sec.R)

"""
    BicircularSystem(prim, sec, μ₁, μ₂, μ, d, R₁, R₂, T, RUNIT, VUNIT, TUNIT, name)
    BicircularSystem(prim::Body, sec::Body)

A BicircularSystem defined by primary, secondary, and tertiary bodies (prim, sec, and ter)
prim::Body      # Primary body
sec::Body       # Secondary body
ter::Body       # Tertiary body
μ₁::Float64     # {km^3/s^2} gravitational parameter of primary body
μ₂::Float64     # {km^3/s^2} gravitational parameter of secondary body
μ₃::Float64     # {km^3/s^2} gravitational parameter of tertiary body
μ::Float64      # {} mass parameter, μ₂/(μ₁+μ₂)
d::Float64      # {km} average distance between two primaries
R₁::Float64     # {km} Radius of primary body
R₂::Float64     # {km} Radius of secondary body
R₃::Float64     # {km} Radius of tertiary body
T::Float64      # {s} sidereal orbital period
RUNIT::Float64  # {km} distance normalizing parameter
VUNIT::Float64  # {km/s} velocity normalizing parameter
TUNIT::Float64  # {s} time normalizing parameter
name::String    # name of system (e.g. "Earth/Moon/Sun")
"""
struct BicircularSystem
    prim::Body      # Primary body
    sec::Body       # Secondary body
    ter::Body       # Tertiary body
    mu1::Float64    # {km^3/s^2} gravitational parameter of primary body
    mu2::Float64    # {km^3/s^2} gravitational parameter of secondary body
    mu3::Float64    # {km^3/s^2} gravitational parameter of tertiary body
    mu::Float64     # {} mass parameter, μ₂/(μ₁+μ₂)
    m3::Float64     # {} normalized mass of tertiary body
    n3::Float64     # {} normalized mean motion of tertiary body
    μ₁::Float64     # {km^3/s^2} gravitational parameter of primary body
    μ₂::Float64     # {km^3/s^2} gravitational parameter of secondary body
    μ₃::Float64     # {km^3/s^2} gravitational parameter of tertiary body
    μ::Float64      # {} mass parameter, μ₂/(μ₁+μ₂)
    m₃::Float64     # {} normalized mass of tertiary body
    n₃::Float64     # {} normalized mean motion of tertiary body
    d::Float64      # {km} average distance between two primaries
    R1::Float64     # {km} Radius of primary Body
    R2::Float64     # {km} Radius of Secondary Body
    R3::Float64     # {km} Radius of tertiary body
    T::Float64      # {s} sidereal orbital period of prim and sec about each other
    RUNIT::Float64  # {km} distance normalizing parameter
    TUNIT::Float64  # {s} time normalizing parameter
    VUNIT::Float64  # {km/s} velocity normalizing parameter
    AUNIT::Float64  # {km^2/s} acceleration normalizing parameter
    name::String    # name of system (e.g. "Earth/Moon/SUN")
end
BicircularSystem(prim::Body, sec::Body, ter::Body) = BicircularSystem(prim, sec, ter,
    prim.m*G, sec.m*G, ter.m*G, sec.m/(prim.m+sec.m), ter.m/(prim.m+sec.m), sec.T/prim.T,
    prim.m*G, sec.m*G, ter.m*G, sec.m/(prim.m+sec.m), ter.m/(prim.m+sec.m), sec.T/prim.T,
    sec.a, prim.R, sec.R, ter.R, sec.T, sec.a, sec.T/2π, sec.a/(sec.T/2π), sec.a/(sec.T/2π)^2,
    string(prim.name,"/",sec.name,"/",ter.name))

### Setting Parameters for Various Bodies and Systems
# Values from wikipedia
SUN         = Body(1.98847e30,  695700,     0.0,            0.0,            "Sun",  :yellow) # The sun does not orbit a central body

## Planets
# Values from:
# m, R, T - https://ssd.jpl.nasa.gov/?planet_phys_par,
# a - https://ssd.jpl.nasa.gov/txt/p_elem_t2.txt
MERCURY     = Body(0.330114e24,             2439.4,     0.38709843*AU,  0.2408467*JY,   "Mercury",  :grey)
VENUS       = Body(4.86747e24,              6051.8,     0.72332102*AU,  0.61519726*JY,  "Venus",    :orange)
EARTH       = Body(5.97237e24,              6371.0084,  1.00000018*AU,  1.0000174*JY,   "Earth",    :blue)
MARS        = Body(0.641712e24,             3389.50,    1.52371243*AU,  1.8808476*JY,   "Mars",     :red)
JUPITER     = Body(1898.187e24,             69911,      5.20248019*AU,  11.862615*JY,   "Jupiter",  :red)
SATURN      = Body(568.3174e24,             58232,      9.54149883*AU,  29.447498*JY,   "Saturn",   :yellow4)
URANUS      = Body(86.8127e24,              25362,      19.18797948*AU, 84.016846*JY,   "Uranus",   :cyan)
NEPTUNE     = Body(102.4126e24,             24622,      30.06952752*AU, 164.79132*JY,   "Neptune",  :blue4)
PLUTO       = Body(0.013030e24,             1188.3,     39.48686035*AU, 247.92065*JY,   "Pluto",    :grey4)

## Satellites
# Values for R, a, and T for all satellites from https://ssd.jpl.nasa.gov/?sat_phys_par, https://ssd.jpl.nasa.gov/?sat_elem
MOON        = Body(4902.801/G,              1737.5,     384400,         27.322*JD,      "Moon",     :grey)
# Values for m from ftp://ssd.jpl.nasa.gov/pub/eph/satellites/nio/LINUX_PC/mar097.2100-2600.txt,
PHOBOS      = Body(7.087546066894452E-04/G, 11.1,       9376.,          0.3189*JD,      "Phobos",   :grey)
DEIMOS      = Body(9.615569648120313E-05/G, 6.2,        23458.,         1.2624*JD,      "Deimos",   :grey)
# Values for m from ftp://ssd.jpl.nasa.gov/pub/eph/satellites/nio/LINUX_PC/jup310xl.txt
IO          = Body(5.959924010272514E+03/G, 1821.6,     421800.,        1.769*JD,       "Io",       :grey)
EUROPA      = Body(3.202739815114734E+03/G, 1560.8,     671100.,        3.551*JD,       "Europa",   :grey)
GANYMEDE    = Body(9.887819980080976E+03/G, 2631.2,     1070400.,       7.155*JD,       "Ganymede", :grey)
CALLISTO    = Body(7.179304867611079E+03/G, 2410.3,     1882700.,       16.69*JD,       "Callisto", :grey)
AMALTHEA    = Body(1.487604677404272E-01/G, 83.45,      181400.,        0.498*JD,       "Amalthea", :grey)
# Values for m from ftp://ssd.jpl.nasa.gov/pub/eph/satellites/nio/LINUX_PC/sat427l.txt
MIMAS       = Body(2.503617062809250E+00/G, 198.20,     185539.,        0.942*JD,       "Mimas",    :grey)
ENCELADUS   = Body(7.210497553340731E+00/G, 252.10,     238042.,        1.370*JD,       "Enceladus",:grey)
TETHYS      = Body(4.121405263872402E+01/G, 533.00,     294672.,        1.888*JD,       "Tethys",   :grey)
DIONE       = Body(7.311617801921636E+01/G, 561.70,     377415.,        2.737*JD,       "Dione",    :grey)
RHEA        = Body(1.539409077211430E+02/G, 764.30,     527068.,        4.518*JD,       "Rhea",     :grey)
TITAN       = Body(8.978137369591670E+03/G, 2574.73,    1221865.,       15.95*JD,       "Titan",    :grey)
HYPERION    = Body(3.704182596063880E-01/G, 135.00,     1500933.,       21.28*JD,       "Hyperion", :grey)
IAPETUS     = Body(1.205081845217891E+02/G, 735.60,     3560854.,       79.33*JD,       "Iapetus",  :grey)
PHOEBE      = Body(5.581081743011904E-01/G, 106.50,     12947918.,      548.02*JD,      "Phoebe",   :grey)
# Values for m from ftp://ssd.jpl.nasa.gov/pub/eph/satellites/nio/LINUX_PC/ura111.xl.txt
ARIEL       = Body(8.346344431770477E+01/G, 578.9,      190900.,        2.520*JD,       "Ariel",    :grey)
UMBRIEL     = Body(8.509338094489388E+01/G, 584.7,      266000.,        4.144*JD,       "Umbriel",  :grey)
TITANIA     = Body(2.269437003741248E+02/G, 788.9,      436300.,        8.706*JD,       "Titania",  :grey)
OBERON      = Body(2.053234302535623E+02/G, 761.4,      583500.,        13.46*JD,       "Oberon",   :grey)
MIRANDA     = Body(4.319516899232100E+00/G, 235.8,      129900.,        1.413*JD,       "Miranda",  :grey)
# Values for m from ftp://ssd.jpl.nasa.gov/pub/eph/satellites/nio/LINUX_PC/nep081xl.txt
TRITON      = Body(1.427598140725034E+03/G, 1353.4,     354759.,        5.877*JD,       "Triton",   :grey)
NEREID      = Body(2.060000000000000E+03/G, 170,        5513818.,       360.13*JD,      "Nereid",   :grey)
PROTEUS     = Body(2.588564729289845E+00/G, 210,        117646.,        1.122*JD,       "Proteus",  :grey)
# Values for m from ftp://ssd.jpl.nasa.gov/pub/eph/satellites/nio/LINUX_PC/plu055l.txt
CHARON      = Body(1.062509269522026E+02/G, 603.6,      19591.,         6.387*JD,       "Charon",   :grey)
NIX         = Body(2.150552267969335E-03/G, 23.0,       48671.,         24.85*JD,       "Nix",      :grey)
HYDRA       = Body(3.663917107480563E-03/G, 30.5,       48671.,         38.20*JD,       "Hydra",    :grey)
KERBEROS    = Body(4.540734312735987E-04/G, 14.0,       64698.,         32.17*JD,       "Kerberos", :grey)
STYX        = Body(2.000000000000000E-20/G, 10.0,       57729.,         20.16*JD,       "Styx",     :grey)

"""
    sun_mercury()
"""
sun_mercury() = System(SUN, MERCURY)
"""
    sun_venus()
"""
sun_venus() = System(SUN, VENUS)
"""
    sun_earth()
"""
sun_earth() = System(SUN, EARTH)
"""
    earth_moon()
"""
earth_moon() = System(EARTH, MOON)
"""
    sun_mars()
"""
sun_mars() = System(SUN, MARS)
"""
    sun_jupiter()
"""
sun_jupiter() = System(SUN, JUPITER)
"""
    sun_saturn()
"""
sun_saturn() = System(SUN, SATURN)
"""
    sun_uranus()
"""
sun_uranus() = System(SUN, URANUS)
"""
    sun_neptune()
"""
sun_neptune() = System(SUN, NEPTUNE)
"""
    jupiter_europa()
"""
jupiter_europa() = System(JUPITER, EUROPA)
"""
    saturn_enceladus()
"""
saturn_enceladus() = System(SATURN, ENCELADUS)

"""
    sun_earth_moon()
"""
earth_moon_sun() = BicircularSystem(EARTH, MOON, SUN)

## Extra values from ephemeris files
# MARS        = Body(4.282837362069909E+04/G,
# SYSTEM      = Body(4.282837442560939E+04/G,
#
# JUPITER     = Body(1.266865341960128E+08/G,
# SYSTEM      = Body(1.267127641334463E+08/G,
#
# SATURN      = Body(3.793120615901047E+07/G,
# SYSTEM      = Body(3.794058484179918E+07/G,
# HELENE      = Body(4.551624250415933E-04/G,
# TELESTO     = Body(0.000000000000000E+00/G,
# CALYPSO     = Body(0.000000000000000E+00/G,
# POLYDEUCES  = Body(0.000000000000000E+00/G,
# METHONE     = Body(0.000000000000000E+00/G,
#
# URANUS      = Body(5.793951322279009E+06/G,
# SYSTEM      = Body(5.794556465751799E+06/G,
#
# NEPTUNE     = Body(6.835099502439672E+06/G,
# SYSTEM      = Body(6.836527100580397E+06/G,
# NAIAD       = Body(8.028519741022843E-03/G,
# THALASSA    = Body(2.358873197992170E-02/G,
# DESPINA     = Body(1.179477866876596E-01/G,
# GALATEA     = Body(1.905015637145090E-01/G,
# LARISSA     = Body(2.551116091803392E-01/G,
#
# PLUTO       = Body(8.693390780381926E+02/G,
# SYSTEM      = Body(9.755962735332020E+02/G,


text:
# Download BedMachine
using Pkg
Pkg.activate(".")
using NCDatasets, Plots, Parameters

!ispath("../data") && mkdir("../data")
!ispath("../output") && mkdir("../output")

# smooth_surface = 2000.0 # smoothing radius, in meters
nsm_s  = 2 # surface
nsm_b  = 0 # bed

run_id = "run_ds_20" # DEBUG: define the output directory within ../output/

!ispath("../output/$run_id") && mkdir("../output/$run_id")

const rw = 1000.0
const ri = 910.0

url = "file://" * homedir() * "/itet-stor/glazioarch/GlacioData/BedMachine_Antarctica/168596330/BedMachineAntarctica_2019-11-05_v01.nc"
# url = "file://" * "/Volumes/luraess/glazioarch/GlacioData/BedMachine_Antarctica/168596330/BedMachineAntarctica_2019-11-05_v01.nc"

filename = if !isfile(joinpath("../data", splitdir(url)[2]))
    print("Downloading data ...")
    download(url, joinpath("../data", splitdir(url)[2]))
    println("done.")
else
    joinpath("../data", splitdir(url)[2])
end

# if !(@isdefined ds) # DEBUG: should be done differently to allow for changing ds while the file already downloaded
    print("Loading NC file ... ")
    ds = NCDataset(filename)

    # downscale = 2 finest possible on 16GB
    # downscale = 1 ok on 32GB ram (needs around 16GB)
    downscale = 30

    x  = ds["x"][1:downscale:end]
    y  = ds["y"][1:downscale:end]
    x  = x[1]:x[2]-x[1]:x[end]
    y  = y[1]:y[2]-y[1]:y[end]
    mask_   = ds["mask"][1:downscale:end, 1:downscale:end]
    mask    = (mask_.==2) .| (mask_.==4) .| (mask_.==1) # route water over(under): grounded ice or lake Vostok or ice-free-land
    surface = ds["surface"][1:downscale:end, 1:downscale:end]
    bed     = convert(Matrix{Float32}, ds["bed"][1:downscale:end, 1:downscale:end])

    # Pad the data to make it GPU compatible with good performance # DEBUG: do this until ParallelStencil is fixed: see https://github.com/samo-lin/ParallelStencil.jl/issues/42
    olx = rem(length(x), 32)
    oly = rem(length(y), 8 )
    x1=1; xE=0; if (olx!=0) x1 = Int(ceil(olx/2)); xE = olx-x1+1; end
    y1=1; yE=0; if (oly!=0) y1 = Int(ceil(oly/2)); yE = oly-y1+1; end

    x  = x[x1:end-xE]
    y  = y[y1:end-yE]
    mask    = mask[x1:end-xE, y1:end-yE]
    surface = surface[x1:end-xE, y1:end-yE]
    bed     = bed[x1:end-xE, y1:end-yE]

    println("done")
# end
println("grid size: nx=$(length(x)), ny=$(length(y))")

# if smooth_surface>0
#     dn = round(Int, smooth_surface/step(x))
#     if dn>0
#         #smooth
#         @warn "todo"
#     end
# end

@views function smooth!(A)
    A[2:end-1,2:end-1] .= A[2:end-1,2:end-1] .+ 1.0./6.1.*(diff(diff(A[:,2:end-1],dims=1),dims=1) .+ diff(diff(A[2:end-1,:],dims=2),dims=2))
    return
end

print("Smoothing DEM [surface $nsm_s step(s), bed $nsm_b step(s)] ... ")
for ism = 1:nsm_s
    smooth!(surface)
end
for ism = 1:nsm_b
    smooth!(bed)
end
println("done")

surf = copy(surface)
surf[.!mask] .= NaN

## Now run flow routing with surface and bed

# using WhereTheWaterFlows, PyPlot
# const WWF = WhereTheWaterFlows
# inds = [45:55,45:55]
# mask = mask[inds...]
# hydpot = (rw*bed .+ (surf.-bed)*ri)[inds...]

# hydpot = (rw*bed .+ (surf.-bed)*ri)
# @time area, slen, dir, nout, nin, pits, c, bnds = waterflows( hydpot, drain_pits=true, bnd_as_pits=true);
# area[.!mask] .= NaN;
# WWF.plotarea(1:size(area,1), 1:size(area,2), area, pits)
# if !all(getindex.(Ref(mask), pits).==0)
#     @warn "Not all pits were removed"
# end
#

using WhereTheWaterAlsoFlows
#const WWAF = WhereTheWaterAlsoFlows

@time D, F, xcp, ycp = flow_routing2D(x, reverse(y), bed, surface.-bed, mask, mask*1, outdir="../output/$run_id")

# visu
F[.!mask] .= NaN

plt = heatmap(xcp, ycp, log10.(F'), aspect_ratio=1, xlims=(xcp[1], xcp[end]), ylims=(ycp[1], ycp[end]), c=:viridis, framestyle=:box, title="log10 ||Flux||")

# savefig(plot(plt, dpi=300), joinpath(@__DIR__, "../output/$run_id/o$(length(x))_final.png"))

display(plt)

nothing


text:
@inline function oracle_type(::Lob{O}) :: OraOracleTypeNum where {O}
    return O
end

@inline function native_type(::Lob) :: OraNativeTypeNum
    return ORA_NATIVE_TYPE_LOB
end

function close(lob::Lob)
    if lob.is_open
        result = dpiLob_close(lob.handle)
        error_check(context(lob), result)
        lob.is_open = false
    end

    nothing
end

function add_ref(lob::Lob)
    result = dpiLob_addRef(lob.handle)
    error_check(context(lob), result)
end

@inline function check_valid_temp_lob_oracle_type_num(t::OraOracleTypeNum)
    is_valid = ( t == ORA_ORACLE_TYPE_CLOB
              || t == ORA_ORACLE_TYPE_NCLOB
              || t == ORA_ORACLE_TYPE_BLOB)

    @assert is_valid "Oracle type num $t is not allowed for a Temp Lob."
end

@inline function check_valid_lob_oracle_type_num(t::OraOracleTypeNum)
    is_valid = ( t == ORA_ORACLE_TYPE_CLOB
              || t == ORA_ORACLE_TYPE_NCLOB
              || t == ORA_ORACLE_TYPE_BLOB
              || t == ORA_ORACLE_TYPE_BFILE)

    @assert is_valid "Oracle type num $t is not allowed for Lob."
end

function Lob(conn::Connection, lob_type::OraOracleTypeNum)
    check_valid_temp_lob_oracle_type_num(lob_type)

    lob_handle_ref = Ref{Ptr{Cvoid}}()
    result = dpiConn_newTempLob(conn.handle, lob_type, lob_handle_ref)
    error_check(context(conn), result)
    return Lob(conn, lob_handle_ref[], lob_type)
end

# Reading and writing to the Lob in multiples of this size will improve performance.
function chunk_size(lob::Lob) :: UInt32
    chunk_size_ref = Ref{UInt32}()
    result = dpiLob_getChunkSize(lob.handle, chunk_size_ref)
    error_check(context(lob), result)
    return chunk_size_ref[]
end

# Returns the size of the data stored in the Lob. For character Lobs the size is in characters; for binary Lobs the size is in bytes.
function _lob_size(lob::Lob)
    lob_size_ref = Ref{UInt64}()
    result = dpiLob_getSize(lob.handle, lob_size_ref)
    error_check(context(lob), result)
    return lob_size_ref[]
end

function size_in_bytes(lob::Lob)
    @assert !is_character_data(lob) "This Lob holds character data. It's size is in Chars. Use `size_in_chars`."
    return _lob_size(lob)
end

function size_in_chars(lob::Lob)
    @assert is_character_data(lob) "This Lob holds byte data. It's size is in Bytes. Use `size_in_bytes`."
    return _lob_size(lob)
end

@inline function is_character_data(::T) :: Bool where {T<:Lob}
    return is_character_data(T)
end

function is_character_data(::Type{Lob{ORATYPE, T}}) :: Bool where {ORATYPE, T}
    if ORATYPE == ORA_ORACLE_TYPE_BFILE || ORATYPE == ORA_ORACLE_TYPE_BLOB
        return false
    else
        @assert ORATYPE == ORA_ORACLE_TYPE_CLOB || ORATYPE == ORA_ORACLE_TYPE_NCLOB
        return true
    end
end

# based on connection encoding/nencoding,
# returns the maximum number of bytes per char
# for this LOB. (CLOB -> encoding, NCLOB -> nencoding)
@inline function max_bytes_per_char(lob::Lob, size_in_chars::Integer=1) :: UInt64
    size_in_bytes_ref = Ref{UInt64}()
    result = dpiLob_getBufferSize(lob.handle, UInt64(size_in_chars), size_in_bytes_ref)
    error_check(context(lob), result)
    return size_in_bytes_ref[]
end

function Base.write(lob::Lob, data::Ptr{UInt8}, data_length::Integer)
    result = dpiLob_setFromBytes(lob.handle, data, UInt64(data_length))
    error_check(context(lob), result)
    nothing
end

Base.write(lob::Lob, data::Vector{UInt8}) = write(lob, pointer(data), length(data))
Base.write(lob::Lob, data::String) = write(lob, pointer(data), sizeof(data))

function Base.read(blob::Lob{ORA_ORACLE_TYPE_BLOB}) :: Vector{UInt8}
    blob_size = size_in_bytes(blob)
    result = undef_vector(UInt8, blob_size)

    open(blob, "r") do io
        i = 0
        while !eof(io)
            i += 1
            result[i] = read(io, UInt8)
        end
        @assert i == blob_size
    end

    return result
end

for ora_type in (ORA_ORACLE_TYPE_CLOB, ORA_ORACLE_TYPE_NCLOB)
    @eval begin
        function Base.read(clob::Lob{$ora_type}) :: String
            local result::String

            open(clob, "r") do io
                result = read(io, String)
            end

            return result
        end
    end
end

abstract type LobDataType end

struct BinaryLob <: LobDataType end
struct CharacterLob <: LobDataType end

abstract type AbstractLobIO{T<:LobDataType,L<:Lob,RW} <: IO end

#
# BinaryLobIO
#

mutable struct BinaryLobIO{L<:Lob,RW} <: AbstractLobIO{BinaryLob,L,RW}
    lob::L
    buffer::Vector{UInt8}
    buffer_start_pos::UInt64 # uninitialized implies 0. First position is 1.
    buffer_end_pos::UInt64
    pos::UInt64              # unread implies pos=0.
    mark::Int                # io.jl mark interface. Unmarked implies mark = -1.

    function BinaryLobIO(lob::Lob{O,P}, rw::Symbol, buffer::Vector{UInt8}) where {O,P}
        @assert !is_character_data(lob)
        check_lob_open_mode(rw)
        return new{Lob{O}, rw}(lob, buffer, 0, 0, 0, -1)
    end
end

mutable struct CharacterLobIO{L<:Lob,RW} <: AbstractLobIO{CharacterLob,L,RW}
    lob::L
    transfer_buffer::Vector{UInt8}
    transfer_buffer_capacity_in_chars::UInt64
    buffer::String
    buffer_start_pos::UInt64 # in chars
    buffer_end_pos::UInt64   # in chars
    pos::UInt64              # in chars
    buffer_pos_index::UInt64 # buffer string index relative to current position
    mark::Int

    function CharacterLobIO(lob::Lob{O,P}, rw::Symbol, transfer_buffer::Vector{UInt8}) where {O,P}
        @assert is_character_data(lob)
        check_lob_open_mode(rw)

        transfer_buffer_size_in_bytes = length(transfer_buffer)
        local transfer_buffer_capacity_in_chars::UInt64
        let
            bytes_per_char = max_bytes_per_char(lob)
            transfer_buffer_capacity_in_chars = UInt64(div(transfer_buffer_size_in_bytes, bytes_per_char))
        end

        return new{Lob{O}, rw}(lob, transfer_buffer, transfer_buffer_capacity_in_chars, "", 0, 0, 0, 0, -1)
    end
end

@inline function check_lob_open_mode(rw::Symbol)
    # TODO
    #@assert rw == :r || rw == :w || rw == :rw "Invalid mode for Lob IO: $rw."
    @assert rw == :r "Mode not supported: $rw."
end

function new_LobIO(lob::Lob, mode;
               buffer::Union{Nothing, Vector{UInt8}}=nothing,
               buffer_size::Union{Nothing, Integer}=nothing)

    buf = new_LobIO_buffer(lob, buffer, buffer_size)

    if is_character_data(lob)
        return CharacterLobIO(lob, Symbol(mode), buf)
    else
        return BinaryLobIO(lob, Symbol(mode), buf)
    end
end

function Base.open(f::Function, lob::Lob, mode;
              buffer::Union{Nothing, Vector{UInt8}}=nothing,
              buffer_size::Union{Nothing, Integer}=nothing)

    local io::AbstractLobIO = new_LobIO(lob, mode, buffer=buffer, buffer_size=buffer_size)

    try
        f(io)
    finally
    end
end

const APPROX_DEFAULT_CHUNK_SIZE_IN_BYTES = UInt64(1E6) # 1MB

function default_buffer_size_in_bytes(l::Lob)
    chunk_sz = chunk_size(l)
    num_chunks = div(APPROX_DEFAULT_CHUNK_SIZE_IN_BYTES, chunk_sz) + 1
    return chunk_sz*num_chunks
end

function new_LobIO_buffer(lob::Lob, buffer, capacity) :: Vector{UInt8}

    function alloc_buffer(cap) :: Vector{UInt8}
        buf = undef_vector(UInt8, cap)
        sizehint!(buf, cap)
        return buf
    end

    if buffer != nothing && capacity != nothing
        @assert length(buffer) == capacity "Provided capacity ($(Int(capacity))) argument diverges from provided buffer size ($(length(buffer)))."
        return buffer
    end

    if buffer != nothing
        return buffer
    end

    if capacity != nothing
        return alloc_buffer(capacity)
    end

    # both arguments were not given
    return alloc_buffer(default_buffer_size_in_bytes(lob))
end

@inline buffer_capacity_in_bytes(io::BinaryLobIO) = UInt64(length(io.buffer))
@inline buffer_capacity_in_bytes(io::CharacterLobIO) = UInt64(length(io.transfer_buffer))
@inline buffer_capacity_in_chars(io::CharacterLobIO) = io.transfer_buffer_capacity_in_chars

function read_chunk_to_buffer!(io::BinaryLobIO, pos::UInt64)
    @assert pos > 0 "Invalid Lob position: $(Int(pos))."
    max_bytes = buffer_capacity_in_bytes(io)
    buffer_len_ref = Ref{UInt64}(max_bytes)
    result = dpiLob_readBytes(io.lob.handle, pos, max_bytes, pointer(io.buffer), buffer_len_ref)
    error_check(context(io.lob), result)
    actual_bytes_read = buffer_len_ref[]

    @assert actual_bytes_read > 0 "Lob has no data at position $pos."

    io.pos = pos - 1 # put position right before the start of the new data

    io.buffer_start_pos = pos
    io.buffer_end_pos = pos + actual_bytes_read - 1

    nothing
end

function read_chunk_to_buffer!(io::CharacterLobIO, pos::UInt64)
    @assert pos > 0 "Invalid Lob position: $(Int(pos))."

    max_chars = buffer_capacity_in_chars(io)
    max_bytes = buffer_capacity_in_bytes(io)

    buffer_len_ref = Ref{UInt64}(max_bytes)
    result = dpiLob_readBytes(io.lob.handle, pos, max_chars, pointer(io.transfer_buffer), buffer_len_ref)
    error_check(context(io.lob), result)
    actual_bytes_read = buffer_len_ref[]

    @assert actual_bytes_read > 0 "Lob has no data at position $pos."

    # assign new data to buffer string, without truncating transfer_buffer
    io.buffer = unsafe_string(pointer(io.transfer_buffer), actual_bytes_read)

    io.pos = pos - 1 # put position right before the start of the new data
    io.buffer_start_pos = pos
    io.buffer_end_pos = pos + length(io.buffer) - 1 # length(io.buffer) is the actual number of characters in the string
    io.buffer_pos_index = 0

    nothing
end

@inline read_chunk_to_buffer!(io::AbstractLobIO, pos::Integer) = read_chunk_to_buffer!(io, UInt64(pos))

function read_next_chunk_to_buffer_if_exhausted!(io::AbstractLobIO)
    if buffer_data_available(io) == 0
        if io.buffer_start_pos == 0
            # this is the first time we read from the stream
            read_chunk_to_buffer!(io, 1) # read from first position
        else
            read_chunk_to_buffer!(io, io.buffer_end_pos + 1)
        end
    end
end

# Returns the number of chars or bytes available in the IO buffer.
function buffer_data_available(io::AbstractLobIO) :: UInt64

    if is_buffer_uninitialized(io)
        return UInt64(0)
    end

    if is_inside_buffer_position(io)
        # bytes available are the distance between end and current positions
        return io.buffer_end_pos - io.pos
    else
        # buffer is initialized, but nothing was read from it.
        return io.buffer_end_pos - io.buffer_start_pos + 1
    end
end

@inline function is_buffer_uninitialized(io::AbstractLobIO) :: Bool
    if io.buffer_start_pos == 0
        return true
    else
        # buffer must be initialized
        @assert io.buffer_end_pos > 0 "AbstractLobIO in inconsistent state."
        return false
    end
end

@inline function is_inside_buffer_position(io::AbstractLobIO, pos::Integer=io.pos) :: Bool
    return io.buffer_start_pos <= pos && pos <= io.buffer_end_pos
end

@inline function check_inside_buffer_position(io::AbstractLobIO)
    @assert is_inside_buffer_position(io) "AbstractLobIO in inconsistent state."
end

# we reach end-of-file if current position equals the size of the lob.
function Base.eof(io::AbstractLobIO) :: Bool
    if buffer_data_available(io) != 0
        return false
    else
        sz = _lob_size(io.lob) # number of chars or bytes, depending on Lob type
        return sz == 0 || sz == io.pos
    end
end

function prepare_read_next!(io::AbstractLobIO)
    if eof(io)
        throw(EOFError())
    end

    read_next_chunk_to_buffer_if_exhausted!(io)

    # inc position
    io.pos += 1
    check_inside_buffer_position(io)
    nothing
end

function Base.read(io::BinaryLobIO, ::Type{UInt8}) :: UInt8
    prepare_read_next!(io)
    @inbounds return io.buffer[io.pos - io.buffer_start_pos + 1]
end

function Base.read(io::CharacterLobIO, ::Type{Char}) :: Char
    prepare_read_next!(io)
    io.buffer_pos_index = nextind(io.buffer, io.buffer_pos_index)
    return io.buffer[io.buffer_pos_index]
end

function Base.read(io::CharacterLobIO, ::Type{String}) :: String
    prepare_read_next!(io)
    io.buffer_pos_index = nextind(io.buffer, io.buffer_pos_index)

    result = io.buffer[io.buffer_pos_index:end]
    io.pos = io.buffer_end_pos

    # continues till the end of the Lob
    while !eof(io)
        io.buffer_pos_index = 0
        read_chunk_to_buffer!(io, io.buffer_end_pos + 1)
        result *= io.buffer
        io.pos = io.buffer_end_pos
    end

    return result
end

function Base.position(io::AbstractLobIO) :: UInt64
    return io.pos
end

function Base.seek(io::AbstractLobIO, pos::Integer)
    if is_inside_buffer_position(io, pos)
        io.pos = pos # next read position will be pos + 1

        if isa(io, CharacterLobIO)
            char_number_in_buffer = io.pos - io.buffer_start_pos + 1
            char_index_in_buffer = firstindex(io.buffer)

            if char_number_in_buffer > 1
                for i in 2:char_number_in_buffer
                    char_index_in_buffer = nextind(io.buffer, char_index_in_buffer)
                end
            end

            io.buffer_pos_index = char_index_in_buffer
        end
    else
        read_chunk_to_buffer!(io, pos + 1)
    end

    nothing
end


text:
module RequiredKeywordArguments

"Error to be thrown if a required keyword argument is missing"
type RequiredKeywordArgumentError <: Exception
    varname
    reason
end
function Base.showerror(io::IO, e::RequiredKeywordArgumentError)
    text = "RequiredKeywordArgumentError: `$(e.varname)` is a required keyword argument, please provide `$(e.varname) = ...`."
    if !isempty(e.reason)
        text *= " $(e.reason)"
    end
    print(io, text)
end

"Function to be used to declare an argument as a required keyword argument."
RequiredKeywordArgument(varname, reason) = throw(RequiredKeywordArgumentError(varname, reason))
RequiredKeywordArgument(varname; reason="") = RequiredKeywordArgument(varname, reason)
RKA = RequiredKeywordArgument

"Convert a functions symbol argument to the corresponding required keyword argument."
function symbol2kw(sym::Symbol)
    Expr(:kw, sym, Expr(:call, RequiredKeywordArgument, string(sym)))
end
symbol2kw(arg) = arg

"""
With @rka, for a function you can define a required keyword argument by simply
not providing a default value.

This means that

    @rka f(;y) = ...

is equivalent to

    f(;y=RequiredKeywordArgument("y")) = ...

"""
macro rka(func)
    params = func.args[1].args[2]
    params.args = map(symbol2kw, params.args)
    return esc(func)
end

export RequiredKeywordArgument, RKA, @rka, RequiredKeywordArgumentError

end # module RequiredKeywordArguments


text:
module Solve_TF
    include("gausslegendre.jl")
    include("solve_tf_module.jl")
    using LinearAlgebra
    using Match
    using MKL
    using .GaussLegendre
    using .Solve_TF_module

    const BETA0 = 1000.0

    beta(x, xarray, beta_array) = let
        klo = 1
        max = length(xarray)
        khi = max

        # 表の中の正しい位置を二分探索で求める
        @inbounds while khi - klo > 1
            k = (khi + klo) >> 1

            if xarray[k] > x
                khi = k
            else
                klo = k
            end
        end

        # yvec_[i] = f(xvec_[i]), yvec_[i + 1] = f(xvec_[i + 1])の二点を通る直線を代入
        return (beta_array[khi] - beta_array[klo]) / (xarray[khi] - xarray[klo]) * (x - xarray[klo]) + beta_array[klo]
    end

    function boundary_conditions!(solve_tf_val, solve_tf_param)
        a = solve_tf_param.Y0
        solve_tf_val.vec_b_glo[1] = a
        solve_tf_val.vec_b_glo[2] -= a * solve_tf_val.tmp[1]

        b = solve_tf_param.YMAX
        solve_tf_val.vec_b_glo[solve_tf_param.NODE_TOTAL] = b
        solve_tf_val.vec_b_glo[solve_tf_param.NODE_TOTAL - 1] -= b * solve_tf_val.tmp[2]
    end

    function construct(data, xarray, yarray)
        solve_tf_param = Solve_TF_module.Solve_TF_param(
            data.grid_num,
            data.gauss_legendre_integ_num,
            data.grid_num + 1,
            data.usethread,
            data.xmax,
            data.xmin,
            yarray[1],
            yarray[end])

        solve_tf_val = Solve_TF_module.Solve_TF_variables(
            Array{Float64}(undef, solve_tf_param.ELE_TOTAL),
            SymTridiagonal(Array{Float64}(undef, solve_tf_param.NODE_TOTAL), Array{Float64}(undef, solve_tf_param.NODE_TOTAL - 1)),
            Array{Int64, 2}(undef, solve_tf_param.ELE_TOTAL, 2),
            Array{Float64, 2}(undef, solve_tf_param.ELE_TOTAL, 2),
            xarray,
            Array{Float64}(undef, 2),
            Array{Float64}(undef, solve_tf_param.NODE_TOTAL),
            Array{Float64}(undef, solve_tf_param.ELE_TOTAL, 2),
            Array{Float64}(undef, solve_tf_param.NODE_TOTAL),
            Array{Float64}(undef, solve_tf_param.INTEGTABLENUM),
            Array{Float64}(undef, solve_tf_param.INTEGTABLENUM))

        solve_tf_val.x, solve_tf_val.w = GaussLegendre.gausslegendre(solve_tf_param.INTEGTABLENUM)

        make_data!(solve_tf_param, solve_tf_val)
        make_global_matrix!(solve_tf_param, solve_tf_val)

        return solve_tf_param, solve_tf_val
    end

    make_beta(solve_tf_val, yarray) = let
        beta_array = Vector{Float64}(undef, length(solve_tf_val.node_x_glo))
        beta_array[1] = BETA0
        for i in 2:length(solve_tf_val.node_x_glo)
            beta_array[i] = yarray[i] * sqrt(yarray[i] / solve_tf_val.node_x_glo[i])
        end

        return beta_array
    end

    function make_data!(solve_tf_param, solve_tf_val)
        @inbounds for e = 1:solve_tf_param.ELE_TOTAL
            solve_tf_val.node_num_seg[e, 1] = e
            solve_tf_val.node_num_seg[e, 2] = e + 1
        end

        @inbounds for e = 1:solve_tf_param.ELE_TOTAL
            for i = 1:2
                solve_tf_val.node_x_ele[e, i] = solve_tf_val.node_x_glo[solve_tf_val.node_num_seg[e, i]]
            end
        end

        # 各線分要素の長さを計算
        @inbounds for e = 1:solve_tf_param.ELE_TOTAL
            solve_tf_val.length[e] = abs(solve_tf_val.node_x_ele[e, 2] - solve_tf_val.node_x_ele[e, 1])
        end
    end

    function make_element_vector!(solve_tf_param, solve_tf_val, yarray)
        beta_array = make_beta(solve_tf_val, yarray)

        if solve_tf_param.USETHREAD
            # Local節点ベクトルの各成分を計算
            Threads.@threads for e = 1:solve_tf_param.ELE_TOTAL
                for i = 1:2
                    solve_tf_val.vec_b_ele[e, i] =
                        @match i begin
                            1 => GaussLegendre.gl_integ(x -> -beta(x, solve_tf_val.node_x_glo, beta_array) * (solve_tf_val.node_x_ele[e, 2] - x) / solve_tf_val.length[e],
                                                        solve_tf_val.node_x_ele[e, 1],
                                                        solve_tf_val.node_x_ele[e, 2],
                                                        solve_tf_val)

                            2 => GaussLegendre.gl_integ(x -> -beta(x, solve_tf_val.node_x_glo, beta_array) * (x - solve_tf_val.node_x_ele[e, 1]) / solve_tf_val.length[e],
                                                        solve_tf_val.node_x_ele[e, 1],
                                                        solve_tf_val.node_x_ele[e, 2],
                                                        solve_tf_val)

                            _ => 0.0
                        end
                end
            end
        else
            # Local節点ベクトルの各成分を計算
            @inbounds for e = 1:solve_tf_param.ELE_TOTAL
                for i = 1:2
                    solve_tf_val.vec_b_ele[e, i] =
                        @match i begin
                            1 => GaussLegendre.gl_integ(x -> -beta(x, solve_tf_val.node_x_glo, beta_array) * (solve_tf_val.node_x_ele[e, 2] - x) / solve_tf_val.length[e],
                                                        solve_tf_val.node_x_ele[e, 1],
                                                        solve_tf_val.node_x_ele[e, 2],
                                                        solve_tf_val)

                            2 => GaussLegendre.gl_integ(x -> -beta(x, solve_tf_val.node_x_glo, beta_array) * (x - solve_tf_val.node_x_ele[e, 1]) / solve_tf_val.length[e],
                                                        solve_tf_val.node_x_ele[e, 1],
                                                        solve_tf_val.node_x_ele[e, 2],
                                                        solve_tf_val)

                            _ => 0.0
                        end
                end
            end
        end
    end

    function make_global_matrix!(solve_tf_param, solve_tf_val)
        # 要素行列
        mat_A_ele = Array{Float64}(undef, solve_tf_param.ELE_TOTAL, 2, 2)

        # 要素行列の各成分を計算
        @inbounds for e = 1:solve_tf_param.ELE_TOTAL
            for i = 1:2
                for j = 1:2
                    mat_A_ele[e, i, j] = (-1) ^ i * (-1) ^ j / solve_tf_val.length[e]
                end
            end
        end

        tmp_dv = zeros(solve_tf_param.NODE_TOTAL)
        tmp_ev = zeros(solve_tf_param.NODE_TOTAL - 1)

        # 全体行列を生成
        @inbounds for e = 1:solve_tf_param.ELE_TOTAL
            for i = 1:2
                for j = 1:2
                    if solve_tf_val.node_num_seg[e, i] == solve_tf_val.node_num_seg[e, j]
                        tmp_dv[solve_tf_val.node_num_seg[e, i]] += mat_A_ele[e, i, j]
                    elseif solve_tf_val.node_num_seg[e, i] + 1 == solve_tf_val.node_num_seg[e, j]
                        tmp_ev[solve_tf_val.node_num_seg[e, i]] += mat_A_ele[e, i, j]
                    end
                end
            end
        end

        # 全体ベクトルの境界条件処理のために保管しておく
        solve_tf_val.tmp[1] = tmp_ev[1]
        solve_tf_val.tmp[2] = tmp_ev[solve_tf_param.NODE_TOTAL - 1]

        # 全体行列の境界条件処理
        tmp_dv[1] = 1.0
        tmp_ev[1] = 0.0
        tmp_dv[solve_tf_param.NODE_TOTAL] = 1.0
        tmp_ev[solve_tf_param.NODE_TOTAL - 1] = 0.0

        solve_tf_val.mat_A_glo = SymTridiagonal(tmp_dv, tmp_ev)
    end

    function make_global_vector!(solve_tf_param, solve_tf_val)
        solve_tf_val.vec_b_glo = zeros(solve_tf_param.NODE_TOTAL)

        # 全体行列と全体ベクトルを生成
        @inbounds for e = 1:solve_tf_param.ELE_TOTAL
            for i = 1:2
                solve_tf_val.vec_b_glo[solve_tf_val.node_num_seg[e, i]] += solve_tf_val.vec_b_ele[e, i]
            end
        end
    end

    function solvetf!(solve_tf_param, solve_tf_val, yarray)
        # Local節点ベクトルを生成
        make_element_vector!(solve_tf_param, solve_tf_val, yarray)

        # 全体ベクトルを生成
        make_global_vector!(solve_tf_param, solve_tf_val)

        # 境界条件処理
        boundary_conditions!(solve_tf_val, solve_tf_param)

        # 連立方程式を解く
        solve_tf_val.ug = solve_tf_val.mat_A_glo \ solve_tf_val.vec_b_glo

        return solve_tf_val.ug
    end
end


text:
using Documenter, FITSIO

makedocs(
    modules = [FITSIO],
    format = :html,
    sitename = "FITSIO.jl",
    pages    = Any[
        "Introduction" => "index.md",
        "API Reference" => "api.md",
        "Libcfitsio Submodule" => "libcfitsio.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaAstro/FITSIO.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
    julia  = "0.6",
)