text:
@testset "L-BFGS ($T)" for T in [Float32, Float64, Complex{Float32}, Complex{Float64}]
    using LinearAlgebra
    using ProximalAlgorithms: LBFGS, update!
    using RecursiveArrayTools: ArrayPartition, unpack

    Q = T[
        32.0000 13.1000 -4.9000 -3.0000 6.0000 2.2000 2.6000 3.4000 -1.9000 -7.5000
        13.1000 18.3000 -5.3000 -9.5000 3.0000 2.1000 3.9000 3.0000 -3.6000 -4.4000
        -4.9000 -5.3000 7.7000 2.1000 -0.4000 -3.4000 -0.8000 -3.0000 5.3000 5.5000
        -3.0000 -9.5000 2.1000 20.1000 1.1000 0.8000 -12.4000 -2.5000 5.5000 2.1000
        6.0000 3.0000 -0.4000 1.1000 3.8000 0.6000 0.5000 0.9000 -0.4000 -2.0000
        2.2000 2.1000 -3.4000 0.8000 0.6000 7.8000 2.9000 -1.3000 -4.3000 -5.1000
        2.6000 3.9000 -0.8000 -12.4000 0.5000 2.9000 14.5000 1.7000 -4.9000 1.2000
        3.4000 3.0000 -3.0000 -2.5000 0.9000 -1.3000 1.7000 6.6000 -0.8000 2.7000
        -1.9000 -3.6000 5.3000 5.5000 -0.4000 -4.3000 -4.9000 -0.8000 7.9000 5.7000
        -7.5000 -4.4000 5.5000 2.1000 -2.0000 -5.1000 1.2000 2.7000 5.7000 16.1000
    ]

    q = T[
        2.9000,
        0.8000,
        1.3000,
        -1.1000,
        -0.5000,
        -0.3000,
        1.0000,
        -0.3000,
        0.7000,
        -2.1000,
    ]

    xs = [
        T[1.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09],
        T[0.09, 1.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08],
        T[0.08, 0.09, 1.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07],
        T[0.07, 0.08, 0.09, 1.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06],
        T[0.06, 0.07, 0.08, 0.09, 1.0, 0.01, 0.02, 0.03, 0.04, 0.05],
    ]

    dirs_ref = [
        T[
            -3.476000000000000e+01,
            -1.367700000000000e+01,
            2.961000000000000e+00,
            3.756000000000000e+00,
            -5.618000000000001e+00,
            -1.571000000000000e+00,
            -4.121000000000000e+00,
            -3.709000000000000e+00,
            4.010000000000000e-01,
            7.639999999999999e+00,
        ],
        T[
            -6.861170733797231e-01,
            -1.661270665201917e+00,
            2.217225828759783e-01,
            5.615134140894827e-01,
            -1.922426760799171e-01,
            -8.961101045874649e-02,
            -3.044802963260585e-01,
            -1.996235459345302e-01,
            1.267604425710271e-01,
            3.360845247013288e-01,
        ],
        T[
            -1.621334774299757e-01,
            2.870743130038511e-01,
            -5.485761164147891e-01,
            9.992734938824949e-02,
            -1.332550298134261e-02,
            5.326252573648003e-02,
            -6.299408068289100e-02,
            1.525398352758626e-02,
            -7.776943954825602e-02,
            -2.335884953507600e-02,
        ],
        T[
            -2.008976150849174e-01,
            2.237224648542354e-01,
            4.811889625788801e-02,
            -6.855884193567087e-01,
            -2.729265954345345e-02,
            3.651730112313705e-02,
            6.325330777317102e-02,
            2.871281112230844e-02,
            -1.285590864125103e-01,
            -3.204963735369062e-03,
        ],
        T[
            -2.317011191832649e-01,
            2.980080835636926e-02,
            -1.267017945785352e-01,
            4.328230970765587e-02,
            -2.437461022925742e-01,
            1.349716200511426e-02,
            -7.155992987801297e-04,
            -3.513449694839536e-03,
            -5.603489763638488e-02,
            5.612114259243499e-02,
        ],
    ]

    @testset "Arrays" begin
        mem = 3
        x = zeros(T, 10)
        H = LBFGS(x, mem)

        x = xs[1]
        grad = Q * x + q
        dir = -(H * grad)

        @test dir ≈ dirs_ref[1]

        for i = 2:5
            x_prev = x
            grad_prev = grad

            x = xs[i]
            grad = Q * x + q

            update!(H, x - x_prev, grad - grad_prev)
            mul!(dir, H, -grad)

            @test dir ≈ dirs_ref[i]
        end
    end

    @testset "ArrayPartition" begin
        mem = 3
        x = ArrayPartition(zeros(T, 10), zeros(T, 10))
        H = LBFGS(x, mem)

        x = ArrayPartition(xs[1], xs[1])
        temp = Q * unpack(x, 1) + q
        grad = ArrayPartition(temp, temp)
        dir = -(H * grad)

        @test unpack(dir, 1) ≈ dirs_ref[1]
        @test unpack(dir, 2) ≈ dirs_ref[1]

        for i = 2:5
            x_prev = x
            grad_prev = grad

            x = ArrayPartition(xs[i], xs[i])
            temp = Q * unpack(x, 1) + q
            grad = ArrayPartition(temp, temp)

            update!(H, x - x_prev, grad - grad_prev)
            mul!(dir, H, -grad)

            @test unpack(dir, 1) ≈ dirs_ref[i]
            @test unpack(dir, 2) ≈ dirs_ref[i]
        end
    end
end


text:
using EvolutionaryModelingTools
using EvolutionaryModelingTools: sample
using Test

@testset "EvolutionaryModelingTools" begin
    @testset "Examples" begin
        # Examples are test with hand-code simulations
        @testset "Logistic growth with noise" begin
            include("example/logistic.jl")
        end
        @testset "Competitive Dynamics with noise and mutation" begin
            include("example/competition.jl")
        end
        @testset "SIR model with noise" begin
            include("example/sir.jl")
        end
    end

    @testset "Scalar" begin
        include("scalar.jl")
    end

    @testset "Sample" begin
        @test sample(1, 1) == CartesianIndex()
        @test sample(fill(1), 1) == CartesianIndex()
        @test sample([1, 2], 2) == CartesianIndex(2)
        @test sample([1 2], 2) == CartesianIndex(1, 2)
        @test sample(reshape(1:4, 2, 2), 1) == CartesianIndex(1, 1)
        @test sample(reshape(1:4, 2, 2), 2) == CartesianIndex(2, 1)
        @test sample(reshape(1:4, 2, 2), 4) == CartesianIndex(1, 2)
        @test sample(reshape(1:4, 2, 2), 7) == CartesianIndex(2, 2)
    end

    @testset "Tools" begin
        include("tools.jl")
    end
end

ext: jl


text:
# We support two-way mapping as in different problems both might be needed.
# In the matrix representation rows are vertices and columns are hyperedges
# n is number of vertices and k is number of hyperedges
# if he[vertex, hyperedge] returns nothing this means that the vertex
# is not present in the hyperedge
# if he[vertex, hyperedge] returns a real number it is a weight
# of the vertex in this hyperedge

"""
    Hypergraph{T} <: AbstractMatrix{Union{T, Nothing}}

A hypergraph storing information about vertices and hyperedges.

**Constructors**

    Hypergraph{T}(n,k) where {T<:Real}
    Hypergraph{T,V}(n, k) where {T<:Real, V}
    Hypergraph{T,V,E}(n, k) where {T<:Real, V, E}

Construct a hypergraph with a given number of vertices and hyperedges.
Optionally, values of type `V` can be stored at vertices and values of type `E`
can be stored at hyperedges.

    Hypergraph(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real}
    Hypergraph{V}(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real, V}
    Hypergraph{V, E}(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real, V, E}

Construct a hypergraph using its matrix representation.
In the matrix representation rows are vertices and columns are hyperedges.
Optionally, values of type `V` can be stored at vertices and values of type `E`
can be stored at hyperedges.

**Arguments**

* `T` : type of weight values stored in the hypergraph
* `V` : type of values stored in the vertices of the hypergraph
* `E` : type of values stored in the edges of the hypergraph
* `n` : number of vertices
* `k` : number of hyperedges
* `m` : a matrix representation rows are vertices and columns are hyperedges

"""
struct Hypergraph{T,V,E} <: AbstractMatrix{Union{T, Nothing}}
    v2he::Vector{Dict{Int,T}}
    he2v::Vector{Dict{Int,T}}
    v_meta::Vector{Union{V,Nothing}}
    he_meta::Vector{Union{E,Nothing}}
    Hypergraph{T,V,E}(n, k) where {T<:Real, V, E} =
        new{T,V,E}([Dict{Int,T}() for i in 1:n], [Dict{Int,T}() for i in 1:k],
                          Vector{Union{V,Nothing}}(nothing, n), Vector{Union{E,Nothing}}(nothing, k))
end

Hypergraph{T,V}(n, k) where {T<:Real, V} = Hypergraph{T,V,Nothing}(n, k)
Hypergraph{T}(n, k) where {T<:Real} =  Hypergraph{T,Nothing,Nothing}(n, k)

function Hypergraph{V, E}(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real, V, E}
    n, k = size(m)
    h = Hypergraph{T, V, E}(n,k)
    h .= m
    h
end

function Hypergraph{V}(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real, V}
    Hypergraph{V, Nothing}(m)
end

function Hypergraph(m::AbstractMatrix{Union{T, Nothing}}) where {T<:Real}
    Hypergraph{Nothing, Nothing}(m)
end

"""
    Base.size(h::Hypergraph)

Returns the size of Hypergraph `h`.
The result is a tuple of the number of vertices and the number of hyperedges

"""
Base.size(h::Hypergraph) = (length(h.v2he), length(h.he2v))

"""
    Base.getindex(h::Hypergraph, idx::Vararg{Int,2})

Returns a value for a given vertex-hyperedge pair `idx` for a hypergraph `h`.
If a vertex does not belong to a hyperedge `nothing` is returned.

"""
@inline function Base.getindex(h::Hypergraph, idx::Vararg{Int,2})
    @boundscheck checkbounds(h, idx...)
    get(h.v2he[idx[1]], idx[2], nothing)
end

"""
    Base.setindex!(h::Hypergraph, ::Nothing, idx::Vararg{Int,2})

Removes a vertex from a given hyperedge for a hypergraph `h` and a given vertex-hyperedge pair `idx`.
Note that trying to remove a vertex from a hyperedge when it is not present will not throw an error.

"""
@inline function Base.setindex!(h::Hypergraph, ::Nothing, idx::Vararg{Int,2})
    @boundscheck checkbounds(h, idx...)
    pop!(h.v2he[idx[1]], idx[2], nothing)
    pop!(h.he2v[idx[2]], idx[1], nothing)
    h
end

"""
    Base.setindex!(h::Hypergraph, v::Real, idx::Vararg{Int,2})

Adds a vertex to a hyperedge (represented by indices `idx`) and assigns value
`v` to be stored with that assignment.

"""
@inline function Base.setindex!(h::Hypergraph, v::Real, idx::Vararg{Int,2})
    @boundscheck checkbounds(h, idx...)
    h.v2he[idx[1]][idx[2]] = v
    h.he2v[idx[2]][idx[1]] = v
    h
end

"""
    getvertices(h::Hypergraph, he_id::Int)

Returns vertices from a hypergraph `a` for a given hyperedge `he_id`.

"""
@inline getvertices(h::Hypergraph, he_id::Int) = h.he2v[he_id]

"""
    gethyperedges(h::Hypergraph, v_id::Int)

Returns hyperedges for a given vertex `v_id` in a hypergraph `h`.

"""
@inline gethyperedges(h::Hypergraph, v_id::Int) = h.v2he[v_id]

"""
    add_vertex!(h::Hypergraph{T, V, E}; hyperedges::Dict{Int,T} = Dict{Int,T}(),
                vertex_meta::Union{V,Nothing} nothing undef) where {T <: Real, V, E}

Adds a vertex to a given hypergraph `h`. Optionally, the vertex can be added
to existing hyperedges. The `hyperedges` parameter presents a dictionary
of hyperedge identifiers and values stored at the hyperedges.
Additionally, a value can be stored with the vertex using the `vertex_meta` keyword parameter.

"""
function add_vertex!(h::Hypergraph{T, V, E}; hyperedges::Dict{Int,T} = Dict{Int,T}(),
                     vertex_meta::Union{V, Nothing} = nothing) where {T <: Real, V, E}
    @boundscheck (checkbounds(h,1,k) for k in keys(hyperedges))
    push!(h.v2he,hyperedges)
    ix = length(h.v2he)
    for k in keys(hyperedges)
        h[ix,k]=hyperedges[k]
    end
    push!(h.v_meta, vertex_meta)
    ix
end

"""
    add_hyperedge!(h::Hypergraph{T, V, E}; vertices::Dict{Int,T} = Dict{Int,T}(),
                   hyperedge_meta::Union{E,Nothing}=nothing) where {T <: Real, V, E}

Adds a hyperedge to a given hypergraph `h`. Optionally, existing vertices can be added
to the created hyperedge. The paramater `vertices` represents a dictionary
of vertex identifiers and values stored at the hyperedges.
Additionally, a value can be stored with the hyperedge using the `hyperedge_meta` keyword parameter.

"""
function add_hyperedge!(h::Hypergraph{T, V, E}; vertices::Dict{Int,T} = Dict{Int,T}(),
                        hyperedge_meta::Union{E,Nothing}=nothing) where {T <: Real, V, E}
    @boundscheck (checkbounds(h,k,1) for k in keys(vertices))
    push!(h.he2v,vertices)
    ix = length(h.he2v)
    for k in keys(vertices)
        h[k,ix]=vertices[k]
    end
    push!(h.he_meta, hyperedge_meta)
    ix
end

"""
    set_vertex_meta!(h::Hypergraph{T, V, E}, new_value::Union{V,Nothing}, id::Int)
        where {T <: Real, V, E}

Sets a new meta value `new_value` for the vertex `id` in the hypegraph `h`.

"""
function set_vertex_meta!(h::Hypergraph{T, V, E}, new_value::Union{V,Nothing}, id::Int) where {T <: Real, V, E}
    checkbounds(h.v_meta, id)
    h.v_meta[id] = new_value
    h.v_meta
end

"""
    get_vertex_meta(h::Hypergraph{T, V, E}, id::Int) where {T <: Real, V, E}

Returns a meta value stored at the vertex `id` in the hypergraph `h`.

"""
function get_vertex_meta(h::Hypergraph{T, V, E}, id::Int) where {T <: Real, V, E}
    checkbounds(h.v_meta, id)
    h.v_meta[id]
end

"""
    set_hyperedge_meta!(h::Hypergraph{T, V, E}, new_value::Union{E,Nothing}, id::Int)
        where {T <: Real, V, E}

Sets a new meta value `new_value` for the hyperedge `id` in the hypegraph `h`.

"""
function set_hyperedge_meta!(h::Hypergraph{T, V, E}, new_value::Union{E,Nothing}, id::Int) where {T <: Real, V, E}
    checkbounds(h.he_meta, id)
    h.he_meta[id] = new_value
    h.he_meta
end

"""
    get_hyperedge_meta(h::Hypergraph{T, V, E}, id::Int) where {T <: Real, V, E}

Returns a meta value stored at the hyperedge `id` in the hypergraph `h`.

"""
function get_hyperedge_meta(h::Hypergraph{T, V, E}, id::Int) where {T <: Real, V, E}
    checkbounds(h.he_meta, id)
    h.he_meta[id]
end

# TODO needs remove_vertex!(h::Hypergraph{T}, v_id::Int)

# TODO needs add_hyperedge!(h::Hypergraph{T}, he_id::Int)

# TODO needs validate_hypergraph!(h::Hypergraph{T})


text:
# @copyright (c) 2019-2020 RWTH Aachen. All rights reserved.
#
# @file src/Maxvol.jl
# @version 1.0.0
# @author Aleksandr Mikhalev
# @date 2020-02-06

module Maxvol

using LinearAlgebra: I, norm
using LinearAlgebra.LAPACK: getrf!
using LinearAlgebra.BLAS: trsm!, ger!, iamax

export maxvol_generic!, maxvol!, rect_maxvol_generic, rect_maxvol

"""
    maxvol_generic!(A, tol=1.05, maxiters=100)

Generic maxvol method, that does not use any **BLAS** or **LAPACK** calls.

Finds good square submatrix. Uses greedy iterative maximization of 1-volume to
find good `r`-by-`r` submatrix in a given `N`-by-`r` matrix `A` of rank `r`.
Returns good submatrix and coefficients of expansion (`N`-by-`r` matrix) of
rows of matrix `A` by rows of good submatrix.

Can be used for special arithmetics, like `BigFloat`, `Rational` or
`Complex{Float16}`.

# Arguments:
- `A::Matrix`: Input matrix on entry and output coefficients on exit.
- `tol::Float64`: Stop when determinant growth is less or equal to this.
- `maxiters::Int`: Maximum number of iterations.

# Returns:
- `piv::Vector{Int}`: Indexes of pivoted rows
- `niters::Int`: Number of performed swap iterations

# Example:
```julia
julia> using Random, LinearAlgebra, Maxvol
julia> rng = MersenneTwister(100);
julia> A = rand(rng, BigFloat, 1000, 100);
julia> C = copy(A);
julia> piv, niters = maxvol_generic!(C);
julia> norm(A-C*A[piv,:]) / norm(A)
2.030657951400512330834848952202721164346464876711701213634530270353170311161736e-76
julia> A = rand(rng, ComplexF64, 1000, 100);
julia> C = copy(A);
julia> piv, niters = maxvol_generic!(C);
julia> norm(A-C*A[piv,:]) / norm(A)
4.863490630095799e-15
```

# See also: [`maxvol!`](@ref)
"""
function maxvol_generic!(A::Matrix{T}, tol::Float64=1.05,
                         maxiters::Int=100) where T <: Number
    # Check input parameters
    if tol < 1
        throw(ArgumentError("Parameter `tol` must be at least 1.0"))
    end
    if maxiters < 0
        throw(ArgumentError("Parameter `maxiters` must be non-negative"))
    end
    # Compute pivoted LU at first
    N, r = size(A)
    ipiv = Vector{typeof(N)}(undef, r)
    basis::Vector{typeof(N)} = 1 : N
    # Cycle over all columns
    for j = 1 : r
        maxi = j
        maxv = abs(A[j, j])
        # Cycle over lower triangular elements in a given column
        # to find maximum absolute value
        for i = j+1 : N
            tmpv = abs(A[i, j])
            if maxv < tmpv
                maxi = i
                maxv = tmpv
            end
        end
        # Throw error if maximum value is 0. This happens only in case of
        # singular matrix
        if maxv == 0
            throw(ArgumentError("Input `A` was a singular matrix"))
        end
        # Swap j-th and maxi-th rows
        if maxi != j
            tmp = A[j, :]
            A[j, :] = A[maxi, :]
            A[maxi, :] = tmp
            tmp = basis[j]
            basis[j] = basis[maxi]
            basis[maxi] = tmp
        end
        # Update ipiv
        ipiv[j] = maxi
        # Normalize lower triangular part
        A[j+1:N, j] /= A[j, j]
        # Eliminate j-th row/column
        for i = j+1 : N
            for k = j+1 : r
                A[i, k] -= A[i, j] * A[j, k]
            end
        end
    end
    # Now solve C A[:r,:] = A[r+1:,:] and overwrite result to A[r+1:,:]
    # A is in LU form, so we need only one triangular solve
    # C L[:r,:] = L[r+1:,:]
    for i = r+1 : N
        # Since diagonal of L is an identity matrix, last column of C is equal
        # to last column of L[r+1:,:], so no need to update anything in that
        # column
        for j = r-1 : -1 : 1
            for k = j+1 : r
                A[i, j] -= A[i, k] * A[k, j]
            end
        end
    end
    # Now make A[:r,:] identity matrix, since A[:r,:]*inv(A[:r,:]) = I
    for i = 1 : r
        for j = 1 : r
            A[i, j] = 0
        end
        A[i, i] = 1
    end
    # Revert swaps, caused by LU
    for i = r : -1 : 1
        if ipiv[i] != i
            tmp = A[i, :]
            A[i, :] = A[ipiv[i], :]
            A[ipiv[i], :] = tmp
        end
    end
    # Exit if maxiters is 0
    if maxiters == 0
        return basis[1:r], 0
    end
    iter::Int = 0
    while true
        maxi, maxj, maxv = -1, -1, tol
        for i = 1 : N
            for j = 1 : r
                tmpv = abs(A[i, j])
                if tmpv > maxv
                    maxv = tmpv
                    maxi = i
                    maxj = j
                end
            end
        end
        if maxi != -1
            basis[maxj] = maxi
            tmp_row = A[maxi, :]
            tmp_row[maxj] -= 1
            tmp_col = A[:, maxj] / A[maxi, maxj]
            for i = 1 : N
                for j = 1 : r
                    A[i, j] -= tmp_row[j] * tmp_col[i]
                end
            end
        end
        iter += 1
        if (maxi == -1) || (iter == maxiters)
            return basis[1:r], iter
        end
    end
end

"""
    maxvol!(A, tol=1.05, maxiters=100)

Finds good square submatrix. Uses greedy iterative maximization of 1-volume to
find good `r`-by-`r` submatrix in a given `N`-by-`r` matrix `A` of rank `r`.
Returns good submatrix and coefficients of expansion (`N`-by-`r` matrix) of
rows of matrix `A` by rows of good submatrix.

Uses vendor-optimized LAPACK, provided by Julia. Supports only 4 input types:
`Float32` (single), `Float64` (double), `ComplexF32` (single complex) and
`ComplexF64` (double complex).

# Arguments:
- `A::Matrix{T}`: Input matrix on entry and output coefficients on exit. `T`
    must be one of `Float32`, `Float64`, `ComplexF32` and `ComplexF64`.
- `tol::Float64`: Stop when determinant growth is less or equal to this.
- `maxiters::Int`: Maximum number of iterations.

# Returns:
- `piv::Vector{Int}`: Indexes of pivoted rows
- `niters::Int`: Number of performed swap iterations

# Example:
```julia
julia> using Random, LinearAlgebra, Maxvol
julia> rng = MersenneTwister(100);
julia> A = rand(rng, Float64, 1000, 100);
julia> C = copy(A);
julia> piv, niters = maxvol!(C);
julia> norm(A-C*A[piv,:]) / norm(A)
2.3975097489579994e-15
julia> A = rand(rng, ComplexF32, 1000, 100);
julia> C = copy(A);
julia> piv, niters = maxvol!(C);
julia> norm(A-C*A[piv,:]) / norm(A)
2.0852597f-6
```

# See also: [`maxvol_generic!`](@ref)
"""
function maxvol!(A::Matrix{T}, tol::Float64=1.05, maxiters::Int=100) where
        T <: Union{Float32, Float64, ComplexF32, ComplexF64}
    # Check input parameters
    if tol < 1
        throw(ArgumentError("Parameter `tol` must be at least 1.0"))
    end
    if maxiters < 0
        throw(ArgumentError("Parameter `maxiters` must be non-negative"))
    end
    # Compute pivoted LU at first
    N, r = size(A)
    basis::Vector{typeof(N)} = 1 : N
    A, ipiv, info = getrf!(A)
    # Throw error if info is not 0. This happens only in case of
    # singular matrix
    if info != 0
        throw(ArgumentError("Input `A` was a singular matrix"))
    end
    # Compute basis rows
    for i = 1 : r
        if ipiv[i] != i
            tmp = basis[i]
            basis[i] = basis[ipiv[i]]
            basis[ipiv[i]] = tmp
        end
    end
    # Now solve C A[:r,:] = A[r+1:,:] and overwrite result to A[r+1:,:]
    # A is in LU form, so we need only one triangular solve
    # C L[:r,:] = L[r+1:,:]
    trsm!('R', 'L', 'N', 'U', one(T), view(A, 1:r, :), view(A, r+1:N, :))
    # Now make A[:r,:] identity matrix, since A[:r,:]*inv(A[:r,:]) = I
    # and revert swaps, caused by LU
    for i = r : -1 : 1
        if ipiv[i] != i
            A[i, :] = A[ipiv[i], :]
        end
        A[ipiv[i], :] .= 0
        A[ipiv[i], i] = 1
    end
    # Exit if maxiters is 0
    if maxiters == 0
        return basis[1:r], 0
    end
    iter::Int = 0
    while true
        tmpi = iamax(A)
        maxj, maxi = fldmod1(tmpi, size(A, 1))
        maxv = abs(A[tmpi])
        if maxv > tol
            basis[maxj] = maxi
            tmp_row = adjoint.(A[maxi, :])
            tmp_row[maxj] -= one(T)
            tmp_col = A[:, maxj] / A[maxi, maxj]
            alpha::T = -one(T)
            ger!(alpha, tmp_col, tmp_row, A)
        end
        iter += 1
        if (maxv <= tol) || (iter == maxiters)
            return basis[1:r], iter
        end
    end
end

"""
    rect_maxvol_generic(A, tol=1.0, maxK=size(A,1), min_add_K=0,
        min_K=size(A,2), start_maxvol_iters=10, identity_submatrix=true)

Generic rect_maxvol method, that does not call **BLAS** or **LAPACK** directly.

Finds good rectangular submatrix. Uses greedy iterative maximization of
2-volume to find good `K`-by-`r` submatrix in a given `N`-by-`r` matrix `A` of
rank `r`. Returns good submatrix and least squares coefficients of expansion
(`N`-by-`K` matrix) of rows of matrix `A` by rows of good submatrix.

Can be used for special arithmetics, like `BigFloat`, `Rational` or
`Complex{Float16}`.

# Arguments:
- `A::Matrix`: Input matrix on entry.
- `tol::Float64`: Stop when volume growth is less or equal to this.
- `maxK::Int`: Maximum number of rows in good submatrix.
- `minK::Int`: Minimum number of rows in good submatrix.
- `min_add_K::Int`: Minimum number of rows to add to the square submatrix.
        Resulting good matrix will have minimum of `r+min_add_K` rows.
- `start_maxvol_iters::Int`: How many iterations of square maxvol (optimization
        of 1-volume) is required to be done before actual rectangular 2-volume
        maximization.
- `identity_submatrix::Bool`: Coefficients of expansions are computed as least
        squares solution. If `identity_submatrix` is True, returned matrix of
        coefficients will have submatrix, corresponding to good rows,
        set to identity.

# Returns:
- `piv::Vector{Int}`: Indexes of pivoted rows
- `C::Matrix`: N-by-size(piv) matrix of coefficients of expansion of all rows
        of A by ``good'' rows `piv`.

# Example:
```julia
julia> using Random, LinearAlgebra, Maxvol
julia> rng = MersenneTwister(100);
julia> A = rand(rng, BigFloat, 1000, 100);
julia> piv, C = rect_maxvol_generic(A);
julia> norm(A-C*A[piv,:]) / norm(A)
1.487413551535496062888389670275572427766752748825675538705145867354156169433079e-76
julia> A = rand(rng, ComplexF64, 1000, 100);
julia> piv, C = rect_maxvol_generic(A);
julia> norm(A-C*A[piv,:]) / norm(A)
3.5140903782712447e-15
```

# See also: [`rect_maxvol`](@ref)
"""
function rect_maxvol_generic(A::Matrix{T}, tol::Float64=1.0,
        maxK::Int=size(A,1), min_add_K::Int=0, minK::Int=size(A,2),
        start_maxvol_iters::Int=10, identity_submatrix::Bool=true) where
        T <: Number
    # Square of tolerance
    tol2 = tol^2
    # N - number of rows, r - number of columns of matrix A
    N, r = size(A)
    # some work on parameters
    if N ≤ r
        return Vector{Int}(1:N), Matrix{T}(I, N, N)
    end
    if maxK > N
        maxK = N
    end
    if maxK < r
        maxK = r
    end
    if minK < r
        minK = r
    end
    if minK > N
        minK = N
    end
    if min_add_K > 0
        minK = max(minK, r+min_add_K)
    end
    if minK > maxK
        minK = maxK
        #raise ValueError('minK value cannot be greater than maxK value')
    end
    # choose initial submatrix and coefficients according to maxvol
    # algorithm
    index = zeros(Int, N)
    chosen = falses(N)
    C = copy(A)
    tmp_index, _ = maxvol_generic!(C, 1.05, start_maxvol_iters)
    # Now matrix A contains coefficients of linear dependancy of old A on
    # pivoted A[tmp_index, :]
    index[1:r] = tmp_index
    chosen[tmp_index] .= true
    # compute square 2-norms of each row in coefficients matrix C
    row_norm_sqr = [chosen[i] ? 0 : norm(C[i, :])^2 for i = 1 : N]
    # find maximum value in row_norm_sqr
    i = argmax(row_norm_sqr)
    K = r
    # augment maxvol submatrix with each iteration
    while (row_norm_sqr[i] > tol2 && K < maxK) || K < minK
        # add i to index and recompute C and square norms of each row
        # by SVM-formula
        K += 1
        index[K] = i
        chosen[i] = true
        c = copy(C[i:i, :])
        #println(size(C), ' ', size(c), ' ', size(c'))
        v = C * (c')
        l::real(T) = 1.0 / (1+v[i]) # v[i] is always real
        c .*= l
        for j = 1 : N
            C[j, :] -= v[j] * c[1, :]
        end
        C = [C l*v]
        #C = [[(C[i, :]-v[i]*(l*c)) v[i]] for i = 1 : N]
        #ger(-l,v,c,a=C,overwrite_a=1)
        #C = np.hstack([C, l*v.reshape(-1,1)])
        row_norm_sqr -= [chosen[i] ? 0 : l*abs(v[i])^2 for i = 1 : N]
        row_norm_sqr[i] = 0
        #row_norm_sqr -= (l*v[:N]*v[:N].conj()).real
        #row_norm_sqr *= chosen
        # find maximum value in row_norm_sqr
        i = argmax(row_norm_sqr)
    end
    # parameter identity_submatrix is True, set submatrix,
    # corresponding to maxvol rows, equal to identity matrix
    if identity_submatrix
        C[index[1:K], :] = Matrix{T}(I, K, K)
    end
    return index[1:K], C
end

"""
    rect_maxvol(A, tol=1.0, maxK=size(A,1), min_add_K=0, min_K=size(A,2),
        start_maxvol_iters=10, identity_submatrix=true)

Finds good rectangular submatrix. Uses greedy iterative maximization of
2-volume to find good `K`-by-`r` submatrix in a given `N`-by-`r` matrix `A` of
rank `r`. Returns good submatrix and least squares coefficients of expansion
(`N`-by-`K` matrix) of rows of matrix `A` by rows of good submatrix.

Uses vendor-optimized LAPACK, provided by Julia. Supports only 4 input types:
`Float32` (single), `Float64` (double), `ComplexF32` (single complex) and
`ComplexF64` (double complex).

# Arguments:
- `A::Matrix`: Input matrix on entry.
- `tol::Float64`: Stop when volume growth is less or equal to this.
- `maxK::Int`: Maximum number of rows in good submatrix.
- `minK::Int`: Minimum number of rows in good submatrix.
- `min_add_K::Int`: Minimum number of rows to add to the square submatrix.
        Resulting good matrix will have minimum of `r+min_add_K` rows.
- `start_maxvol_iters::Int`: How many iterations of square maxvol (optimization
        of 1-volume) is required to be done before actual rectangular 2-volume
        maximization.
- `identity_submatrix::Bool`: Coefficients of expansions are computed as least
        squares solution. If `identity_submatrix` is True, returned matrix of
        coefficients will have submatrix, corresponding to good rows,
        set to identity.

# Returns:
- `piv::Vector{Int}`: Indexes of pivoted rows
- `C::Matrix`: N-by-size(piv) matrix of coefficients of expansion of all rows
        of A by ``good'' rows `piv`.

# Example:
```julia
julia> using Random, LinearAlgebra, Maxvol
julia> rng = MersenneTwister(100);
julia> A = rand(rng, Float64, 1000, 100);
julia> piv, C = rect_maxvol(A);
julia> norm(A-C*A[piv,:]) / norm(A)
1.8426360389674326e-15
julia> A = rand(rng, ComplexF32, 1000, 100);
julia> piv, C = rect_maxvol(C);
julia> norm(A-C*A[piv,:]) / norm(A)
1.8014168f-6
```

# See also: [`rect_maxvol_generic`](@ref)
"""
function rect_maxvol(A::Matrix{T}, tol::Float64=1.0,
        maxK::Int=size(A, 1), min_add_K::Int=0, minK::Int=size(A, 2),
        start_maxvol_iters::Int=10, identity_submatrix::Bool=true) where
        T <: Union{Float32, Float64, ComplexF32, ComplexF64}
    # Square of tolerance
    tol2 = tol^2
    # N - number of rows, r - number of columns of matrix A
    N, r = size(A)
    # some work on parameters
    if N ≤ r
        return Vector{Int}(1:N), Matrix{T}(I, N, N)
    end
    if maxK > N
        maxK = N
    end
    if maxK < r
        maxK = r
    end
    if minK < r
        minK = r
    end
    if minK > N
        minK = N
    end
    if min_add_K > 0
        minK = max(minK, r+min_add_K)
    end
    if minK > maxK
        minK = maxK
        #raise ValueError('minK value cannot be greater than maxK value')
    end
    # choose initial submatrix and coefficients according to maxvol
    # algorithm
    index = zeros(Int, N)
    chosen = falses(N)
    C = copy(A)
    tmp_index, _ = maxvol!(C, 1.05, start_maxvol_iters)
    # Now matrix A contains coefficients of linear dependancy of old A on
    # pivoted A[tmp_index, :]
    index[1:r] = tmp_index
    chosen[tmp_index] .= true
    # compute square 2-norms of each row in coefficients matrix C
    row_norm_sqr = Vector{real(T)}(undef, N)
    row_norm_sqr[:] = [chosen[i] ? 0 : norm(C[i, :])^2 for i = 1 : N]
    # find maximum value in row_norm_sqr
    i = argmax(row_norm_sqr)
    K = r
    # augment maxvol submatrix with each iteration
    while (row_norm_sqr[i] > tol2 && K < maxK) || K < minK
        # add i to index and recompute C and square norms of each row
        # by SVM-formula
        K += 1
        index[K] = i
        chosen[i] = true
        c = Vector{T}(undef, K-1)
        c[:] = C[i, :]'
        v = C * c
        l::real(T) = 1.0 / (1+v[i]) # v[i] is always real
        ger!(T(-l), v, c, C)
        C = [C l*v]
        #C = [[(C[i, :]-v[i]*(l*c)) v[i]] for i = 1 : N]
        #ger(-l,v,c,a=C,overwrite_a=1)
        #C = np.hstack([C, l*v.reshape(-1,1)])
        row_norm_sqr -= [chosen[i] ? 0 : l*abs(v[i])^2 for i = 1 : N]
        row_norm_sqr[i] = 0
        #row_norm_sqr -= (l*v[:N]*v[:N].conj()).real
        #row_norm_sqr *= chosen
        # find maximum value in row_norm_sqr
        i = argmax(row_norm_sqr)
    end
    # parameter identity_submatrix is True, set submatrix,
    # corresponding to maxvol rows, equal to identity matrix
    if identity_submatrix
        C[index[1:K], :] = Matrix{T}(I, K, K)
    end
    return index[1:K], C
end

end # module Maxvol


text:
# using Revise
using POMDPStressTesting
using Distributions
using Parameters

@with_kw mutable struct CategoricalWalk1DParams
    startx::Float64 = 0 # Starting x-position
    threshx::Float64 = 10 # +- boundary threshold
    endtime::Int64 = 30 # Simulate end time
end

# Implement abstract GrayBox.Simulation
@with_kw mutable struct CategoricalWalk1DSim <: GrayBox.Simulation
    params::CategoricalWalk1DParams = CategoricalWalk1DParams() # Parameters
    x::Float64 = 0 # Current x-position
    t::Int64 = 0 # Current time
    distribution::Distribution = Categorical([0.95,0.025,0.025]) # Transition distribution [1,2,3]
end

# Override from GrayBox
GrayBox.environment(sim::CategoricalWalk1DSim) = GrayBox.Environment(:x => sim.distribution)

# Override from GrayBox (NOTE: used with ASTSeedAction)
function GrayBox.transition!(sim::CategoricalWalk1DSim)
    # We sample the environment and apply the transition
    environment::GrayBox.Environment = GrayBox.environment(sim) # Get the environment distributions
    sample::GrayBox.EnvironmentSample = rand(environment) # Sample from the environment
    sim.t += 1 # Keep track of time
    sim.x += sample[:x].value # Move agent using sampled value from input
    return logpdf(sample)::Real # Summation handled by `logpdf()`
end

# Override from GrayBox (NOTE: used with ASTSampleAction)
function GrayBox.transition!(sim::CategoricalWalk1DSim, sample::GrayBox.EnvironmentSample)
    # The environment was sampled for us, and we just apply the transition
    sim.t += 1 # Keep track of time
    sim.x += sample[:x].value # Move agent using sampled value from input
    return logpdf(sample)::Real # Summation handled by `logpdf()`
end

# Override from BlackBox
function BlackBox.initialize!(sim::CategoricalWalk1DSim)
    sim.t = 0
    sim.x = sim.params.startx
end

# Override from BlackBox
BlackBox.distance(sim::CategoricalWalk1DSim) = max(sim.params.threshx - abs(sim.x), 0)

# Override from BlackBox
BlackBox.isevent(sim::CategoricalWalk1DSim) = abs(sim.x) >= sim.params.threshx

# Override from BlackBox
BlackBox.isterminal(sim::CategoricalWalk1DSim) = BlackBox.isevent(sim) || sim.t >= sim.params.endtime

# Override from BlackBox (NOTE: used with ASTSeedAction)
function BlackBox.evaluate!(sim::CategoricalWalk1DSim)
    logprob::Real = GrayBox.transition!(sim) # Step simulation
    d::Real = BlackBox.distance(sim) # Calculate miss distance
    event::Bool = BlackBox.isevent(sim) # Check event indication
    return (logprob::Real, d::Real, event::Bool)
end

# Override from BlackBox (NOTE: used with ASTSampleAction)
function BlackBox.evaluate!(sim::CategoricalWalk1DSim, sample::GrayBox.EnvironmentSample)
    logprob::Real = GrayBox.transition!(sim, sample) # Step simulation given input sample
    d::Real = BlackBox.distance(sim) # Calculate miss distance
    event::Bool = BlackBox.isevent(sim) # Check event indication
    return (logprob::Real, d::Real, event::Bool)
end

function setup_ast(seed=AST.DEFAULT_SEED; solver=PPOSolver)
    # Create gray-box simulation object
    sim::GrayBox.Simulation = CategoricalWalk1DSim()

    # AST MDP formulation object
    # NOTE: Use either {ASTSeedAction} or {ASTSampleAction} (when using TRPO/PPO/CEM, use ASTSampleAction)
    if solver in [TRPOSolver, PPOSolver, CEMSolver]
        mdp::ASTMDP = ASTMDP{ASTSampleAction}(sim)
    else
        mdp = ASTMDP{ASTSeedAction}(sim)
    end
    mdp.params.debug = true # record metrics
    mdp.params.top_k = 10 # record top k best trajectories
    mdp.params.seed = seed # set RNG seed for determinism
    n_iterations = 1000 # number of algorithm iterations

    # Choose a solver (examples of each)
    if solver == RandomSearchSolver
        solver = RandomSearchSolver(n_iterations=n_iterations,
                                    episode_length=sim.params.endtime)
    elseif solver == MCTSPWSolver
        solver = MCTSPWSolver(n_iterations=n_iterations,
                              exploration_constant=1.0, # UCT exploration
                              k_action=1.0, # action widening
                              alpha_action=0.5, # action widening
                              depth=sim.params.endtime) # tree depth (i.e. episode length)
    elseif solver == CEMSolver
        solver = CEMSolver(n_iterations=n_iterations,
                           episode_length=sim.params.endtime)
    elseif solver == TRPOSolver
        solver = TRPOSolver(num_episodes=n_iterations,
                            episode_length=sim.params.endtime,
                            policy_type=:discrete)
    elseif solver == PPOSolver
        solver = PPOSolver(num_episodes=n_iterations,
                           episode_length=sim.params.endtime,
                           policy_type=:discrete)
    end

    # Get online planner (no work done, yet)
    planner = solve(solver, mdp)

    return planner
end

function run_ast(seed=AST.DEFAULT_SEED; kwargs...)
    planner = setup_ast(seed; kwargs...)

    action_trace::Vector{ASTAction} = search!(planner) # work done here
    final_state::ASTState = playback(planner, action_trace, sim->sim.x)
    failure_rate::Float64 = print_metrics(planner)

    return planner, action_trace::Vector{ASTAction}, failure_rate::Float64
end

(planner, action_trace, failure_rate) = run_ast()

nothing # Suppress REPL


text:
function create_generator(typ::Int=CURAND_RNG_PSEUDO_DEFAULT)
    ptr = Ref{curandGenerator_t}()
    @check ccall((:curandCreateGenerator, libcurand),
                 curandStatus_t,
                 (Ptr{curandGenerator_t}, Cint), ptr, typ)
    r = RNG(ptr[], typ)
    finalizer(destroy_generator, r)
    return r
end

function destroy_generator(rng::RNG)
    @check ccall((:curandDestroyGenerator, libcurand),
                 curandStatus_t,
                 (curandGenerator_t,), rng)
end

function get_version()
    ver = Ref{Cint}()
    @check ccall((:curandGetVersion, libcurand),
                 curandStatus_t,
                 (Ref{Cint},), ver)
    return ver[]
end

# TODO: curandSetStream

function set_pseudo_random_generator_seed(rng::RNG, seed::Int64)
    @check ccall((:curandSetPseudoRandomGeneratorSeed, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Clonglong), rng, seed)
end

function set_generator_offset(rng::RNG, offset::Int64)
    @check ccall((:curandSetGeneratorOffset, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Clonglong), rng, offset)
end

function set_generator_ordering(rng::RNG, order::Int)
    @check ccall((:curandSetGeneratorOrdering, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Cint), rng, order)
end

function set_quasi_random_generator_dimensions(rng::RNG, num_dimensions::UInt)
    @check ccall((:curandSetQuasiRandomGeneratorDimensions, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Cuint),
                 rng, num_dimensions)
end

"""
Generate 64-bit quasirandom numbers.
"""
function generate(rng::RNG, arr::CuArray, n::UInt)
    @check ccall((:curandGenerate, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Ptr{UInt32}, Csize_t),
                 rng, arr, length(arr))
    return arr
end

"""
Generate uniformly distributed floats.

Valid RNG types are:
 - CURAND_RNG_QUASI_SOBOL64
 - CURAND_RNG_QUASI_SCRAMBLED_SOBOL64
"""
function generate_long_long(rng::RNG, arr::CuArray)
    @check ccall((:curandGenerateLongLong, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Ptr{Culonglong}, Csize_t),
                 rng, arr, length(arr))
    return arr
end

# uniform
function generate_uniform(rng::RNG, arr::CuArray)
    @check ccall((:curandGenerateUniform, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Ptr{Float32}, Csize_t),
                 rng, arr, length(arr))
    return arr
end

function generate_uniform_double(rng::RNG, arr::CuArray)
    @check ccall((:curandGenerateUniformDouble, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Ptr{Float64}, Csize_t),
                 rng, arr, length(arr))
    return arr
end

# normal
function generate_normal(rng::RNG, arr::CuArray, mean, stddev)
    @check ccall((:curandGenerateNormal, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Ptr{Cfloat}, Csize_t, Cfloat, Cfloat),
                 rng, arr, length(arr), mean, stddev)
    return arr
end

function generate_normal_double(rng::RNG, arr::CuArray, mean, stddev)
    @check ccall((:curandGenerateNormalDouble, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Ptr{Cdouble}, Csize_t, Cdouble, Cdouble),
                 rng, arr, length(arr), mean, stddev)
    return arr
end

# lognormal
function generate_log_normal(rng::RNG, arr::CuArray, mean, stddev)
    @check ccall((:curandGenerateLogNormal, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Ptr{Cfloat}, Csize_t, Cfloat, Cfloat),
                 rng, arr, length(arr), mean, stddev)
    return arr
end

function generate_log_normal_double(rng::RNG, arr::CuArray, mean, stddev)
    @check ccall((:curandGenerateLogNormalDouble, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Ptr{Cdouble}, Csize_t, Cdouble, Cdouble),
                 rng, arr, length(arr), mean, stddev)
    return arr
end

# Poisson
"""Construct the histogram array for a Poisson distribution."""
function create_poisson_distribtion(lambda)
    ptr = Ref{curandDiscreteDistribution_t}()
    @check ccall((:curandCreatePoissonDistribution, libcurand),
                 curandStatus_t,
                 (Cdouble, Ptr{curandDiscreteDistribution_t}),
                 lambda, ptr)
    dist = DiscreteDistribution(ptr[])
    finalizer(destroy_distribution, dist)
    return dist
end

"""Destroy the histogram array for a discrete distribution (e.g. Poisson)."""
function destroy_distribution(dist::DiscreteDistribution)
    @check ccall((:curandDestroyDistribution, libcurand),
                 curandStatus_t,
                 (curandDiscreteDistribution_t,),
                 dist)
end

"""Generate Poisson-distributed unsigned ints."""
function generate_poisson(rng::RNG, arr::CuArray, lambda)
    @check ccall((:curandGeneratePoisson, libcurand),
                 curandStatus_t,
                 (curandGenerator_t, Ptr{Cuint}, Csize_t, Cdouble),
                 rng, arr, length(arr), lambda)
    return arr
end

# seeds
"""Generate the starting state of the generator. """
function generate_seeds(rng::RNG)
    @check ccall((:curandGenerateSeeds, libcurand),
                 curandStatus_t,
                 (curandGenerator_t,), rng)
end

# TODO: curandGetDirectionVectors32
# TODO: curandGetScrambleConstants32
# TODO: curandGetDirectionVectors64
# TODO: curandGetScrambleConstants64


text:
"""
Plain wrapper.

Usage:

```julia
PyCall.PyObject(self::MyModule.MyType) = PyBase.Plain.wrap(self)
```
"""
module Plain

using ..PyBase
using ..PyBase: super
using ..JuliaAPI: JuliaObject, convert_pyindices

struct Shim{T} <: PyBase.AbstractShim
    self::T
end

"""
    PyBase.Plain.wrap(self) :: PyObject
"""
wrap(x::Any) = JuliaObject(Shim(x))

PyBase._unwrap(shim::Shim) = PyBase.super(shim)

end  # module


text:
struct Cache{T<:Number}
    elts::Vector{AbstractArray{T}}
end

"""
    is_in_cache(cache::Cache{T}, L::Int)

Check whether a vector of size L is stored in cache.
"""
function is_in_cache(cache::Cache{T}, L::Int) where T<:Number
    if any(length.(cache.elts) .== L)
        return findfirst(length.(cache.elts) .== L)
    else
        return 0
    end
end

"""
    update_cache!(cache::Cache{T}, M::AbstractArray{T}) where T<:Number

Update cache with new preallocated array M if no array with the same length is
already stored.
"""
function update_cache!(cache::Cache{T}, M::AbstractArray{T}) where T<:Number
    if is_in_cache(cache, length(M)) == 0
        push!(cache.elts, vec(M))
    end
    return cache
end


text:
@testset "1603.design-parking-system.jl" begin
    parking = ParkingSystem(1, 1, 0)
    @test add_car(parking, 1) == true
    @test add_car(parking, 2) == true
    @test add_car(parking, 3) == false
    @test add_car(parking, 1) == false
end


text:
module bfga

    import Base.isless, Base.isequal

    include("Types.jl")
    using .Types

    #using SSA
    #include("SSA.jl")
    #using .SSA

    include("Bf.jl")
    import .BfInterpreter

    mutable struct bfMonster <: Entity
        dna::Array
        fitness
        bonus :: Float64
        age :: Int
        m_length :: Int
        program :: String

        # = new(create_entity(1, nbGenes).dna, nothing, 1, nbGenes)
        bfMonster(dna :: Array) = new(dna, 0, 0,  1, length(dna), join(genesToBfInstructions(dna), ""))
        bfMonster(dna :: Array, age :: Int) = new(dna, 0, 0, age, length(dna), join(genesToBfInstructions(dna), ""))
        #bfMonster(age :: Int) = create_entity(age,nbGenes)
    end

    #nbGenes = 150

    function isless(lhs::bfMonster, rhs::bfMonster)
        #lhs.fitness < rhs.fitness || lhs.bonus > rhs.bonus
        #(lhs.fitness, rhs.bonus) < (rhs.fitness , lhs.bonus) # abs(lhs.fitness) > abs(rhs.fitness)
        (lhs.fitness, rhs.bonus) < (rhs.fitness , lhs.bonus) # abs(lhs.fitness) > abs(rhs.fitness)
    end

    function isequal(lhs::bfMonster, rhs::bfMonster)
        (lhs.fitness, rhs.bonus) == (rhs.fitness , lhs.bonus)
    end

    function create_entity(num, nbgenes)
        babybf = [rand() for i in 1:nbgenes]
        bfMonster(babybf, num)
    end

    function create_entity(dna :: Array)
        bfMonster(deepcopy(dna))
    end

    function isPresent(ent, ens)
        if isempty(ens)
            return false
        elseif ent.program == ens[1].program #ent.dna == ens[1].dna
            return true
        elseif length(ens)== 1
            return false
        else
            return isPresent(ent, ens[2,:])
        end
    end

    function d_aux(x::Array,y::Array)
        sum(x-y).^2
    end

    function distance(x :: bfMonster ,y ::bfMonster)
        l = length(x.dna)
        mini = l
        for i in -1:1
            mini = min(mini, d_aux(x.dna[2:l-1], y.dna[2+i:l-1+i]))
        end
        mini
    end

    function addNew(refSet, population)
        l = []
        k = length(population[1].dna)
        for e in population
            mini = k
            for s in refSet
                mini = min(mini, sum((e.dna-s.dna).^2))
                #mini = min(mini, distance(e,s))
            end
            push!(l, mini)
        end
        index = argmax(l)
    end

    #Constants
    global instructions = ['<', '>', '+', '-', '.', ',', '[', ']'] # , '$', '!', '*']

    global instructions2 = ['<', '>', '+', '-', '.', ',', '[', ']','$', '!', '*' , '#' , '@', '/' , '%']

    function floatToBfInstruction( number ::Float64 )
        @assert (number >= 0.0 && number <= 1.0)
        memoire = number*length(instructions)
        ent = trunc(Int, memoire)
        if ent == memoire
            return instructions[ent]
        else
            return instructions[ent+1]
        end
    end

    function entityToBfInstructions!( ent )
        ent.program = join(genesToBfInstructions(ent.dna), "")
    end

    function genesToBfInstructions( genes :: Array{Float64,1} )
        instructions = []
        for num in genes
            push!(instructions,floatToBfInstruction(num))
        end
        instructions
    end

    function getBfCode(ent)
        ent.program
    end

    function getInstructionsSet()
        mem = BfInterpreter.getInstructionsDict()
        @show mem
        mem
    end

    function crossover(group)
        mem = rand(1:2)
        mem2 = mem%2 +1
        parent1 = group[mem]
        parent2 = group[mem2]
        n1 = length(parent1.dna)
        n2 = length(parent2.dna)
        pos = trunc(Int, rand() * n1)
        minN = min(n1,n2)

        child1 = create_entity(parent1.age, minN)
        child2 = create_entity(parent1.age, minN)
        child1.age = parent1.age + 1
        child2.age = parent1.age + 1

        for i in 1:minN
            if i < pos
                child1.dna[i] = parent2.dna[i]
                child2.dna[i] = parent1.dna[i]
            else
                child2.dna[i] = parent2.dna[i]
                child1.dna[i] = parent1.dna[i]
            end
        end
        mutate(child1, 0.01)
        mutate(child2, 0.01)
        #clearCode!(child1)
        #clearCode!(child2)
        child1.m_length = minN #length(child1.dna)
        child2.m_length = minN #length(child2.dna)
        child1, child2
    end

    #=function crossover(group)
        parent1 = group[1]
        parent2 = group[2]
        n1 = length(parent1.dna)
        n2 = length(parent2.dna)
        pos = trunc(Int, rand() * n1)
        minN = min(n1,n2)

        child1 = create_entity(parent1.age, minN)
        child2 = create_entity(parent1.age, minN)
        child1.age = parent1.age + 1
        child2.age = parent1.age + 1

        for i in 1:minN
            mem = rand(1:2)
            mem2 = (mem%2) +1
            child1.dna[i] = group[mem].dna[i]
            child2.dna[i] = group[mem2].dna[i]
        end
        mutate(child1, 0.01)
        mutate(child2, 0.01)
        #clearCode!(child1)
        #clearCode!(child2)
        child1.m_length = minN #length(child1.dna)
        child2.m_length = minN #length(child2.dna)
        child1, child2
    end=#

    function crossoverGroup(group)
        nEnt = length(group)
        n = length(group[1].dna)
        child1 = child1 = create_entity(1, n)
        child2 = child1 = create_entity(1, n)
        hop  = div(n, nEnt)
        num = 1
        for i in 1:n
            child1.dna[i] = group[num].dna[i]
            child2.dna[i] = group[(nEnt - num +1)].dna[i]
            if i < nEnt && (i%hop) ==0
                num += 1
            end
        end
        mutate(child1, 0.1)
        mutate(child2, 0.1)
        child1, child2
    end

    function mutate(ent, mutationRate)
        m_length = length(ent.dna)
        # Go through each bit.
        for pos in 1:m_length

            # Should this bit mutate?
            if (rand() < mutationRate)
                # Select a mutation type.
                r = rand()
                if (r <= 0.25)
                    # Insertion mutation.
                    # Get shift index.
                    mutationIndex = pos;

                    # Make a copy of the current bit before we mutate it.
                    shiftBit = ent.dna[mutationIndex];

                    # Set random bit at mutation index.
                    ent.dna[mutationIndex] = rand();

                    # Bump bits up or down by 1.
                    up = rand() >= 0.5;
                    if up
                        # Bump bits up by 1.
                        for i in ((mutationIndex + 1 ):m_length)
                            nextShiftBit = ent.dna[i];

                            ent.dna[i] = shiftBit;

                            shiftBit = nextShiftBit;
                        end
                    else
                        # Bump bits down by 1.
                        for i in mutationIndex:-1:1
                            nextShiftBit = ent.dna[i];

                            ent.dna[i] = shiftBit;

                            shiftBit = nextShiftBit;
                        end
                    end
                elseif (r <= 0.5)
                    # Deletion mutation.
                    # Get deletion index.
                    mutationIndex = pos;

                    # Bump bits up or down by 1.
                    up = rand() >= 0.5;
                    if (up)
                        # Bump bits up by 1.
                        for i in mutationIndex:-1:2
                            ent.dna[i] = ent.dna[i - 1];
                        end

                        # Add a new mutation bit at front of genome to replace the deleted one.
                        ent.dna[1] = rand();
                    else
                        # Bump bits down by 1.
                        for i in mutationIndex:(m_length - 1)
                            ent.dna[i] = ent.dna[i + 1];
                        end

                        # Add a new mutation bit at end of genome to replace the deleted one.
                        ent.dna[m_length] = rand();
                    end
                elseif (r <= 0.75)
                    # Shift/rotation mutation.
                    # Bump bits up or down by 1.
                    up = rand() >= 0.5;
                    if (up)
                        # Bump bits up by 1. 1, 2, 3 => 3, 1, 2
                        shiftBit = ent.dna[1];

                        for i in 1:m_length
                            if (i > 1)
                                # Make a copy of the current bit.
                                temp = ent.dna[i];

                                # Set the current bit to the previous one.
                                ent.dna[i] = shiftBit;

                                # Select the next bit to be copied.
                                shiftBit = temp;
                            else
                                # Wrap last bit to front.
                                ent.dna[i] = ent.dna[m_length];
                            end
                        end
                    else
                        # Bump bits down by 1. 1, 2, 3 => 2, 3, 1
                        shiftBit = ent.dna[m_length ];

                        for i in m_length:-1:1
                            if (i < m_length)
                                # Make a copy of the current bit.
                                temp = ent.dna[i];

                                # Set the current bit to the previous one.
                                ent.dna[i] = shiftBit;

                                # Select the next bit to be copied.
                                shiftBit = temp;
                            else
                                # Wrap first bit to end.
                                ent.dna[i] = ent.dna[1];
                            end
                        end
                    end
                else
                    # Replacement mutation.
                    # Mutate bits.
                    mutation = rand();
                    ent.dna[pos] = mutation;
                end
            end
        end
    end

    function expand(ent , size :: Int)
        #println("EXPAND CAAALLLED")
        originalSize = length(ent.dna);
        difference = size - originalSize;

        # Resize the genome array.
        newGenes = []

        if (difference > 0)
            if (rand() < 0.5)
                # Extend at front.
                newGenes = [ [rand() for i in 1:difference]; deepcopy(ent.dna) ]
            else
                # Extend at back.
                newGenes = [ deepcopy(ent.dna); [rand() for i in 1:difference] ]
            end
            ent.dna = newGenes;
        else
            ent.dna = ent.dna[1:size]
            #Array.Resize(ref ent.dna, size);
        end
        ent.m_length = size;
    end

    function clearCode!(ent)
        indexopen = []
        indexclose =[]

        for i in 1:ent.m_length
            if (floatToBfInstruction(ent.dna[i]) == ']' && isempty(indexopen))
                r = rand()
                #=while floatToBfInstruction(r) == ']'
                    r = rand()
                end=#
                if floatToBfInstruction(r) == '['
                    push!(indexopen, i)
                end
                ent.dna[i] = r
            elseif floatToBfInstruction(ent.dna[i]) == ']'
                d = pop!(indexopen)
                if d == i-1
                    r = rand()
                    #=while floatToBfInstruction(r) == ']'
                        r = rand()
                    end=#
                    if floatToBfInstruction(r) == '['
                        push!(indexopen, i)
                    end
                    ent.dna[i] = r
                end
            elseif floatToBfInstruction(ent.dna[i]) == '['
                push!(indexopen, i)
            end
        end

        #=if !isempty(indexopen)
            for e in indexopen
                r = rand()
                while (floatToBfInstruction(r) == ']' ||  floatToBfInstruction(r) == '[')
                    r = rand()
                end
                ent.dna[e] = r
            end
        end=#
    end

end


text:
doc"""
    Logistic(μ,θ)

The *Logistic distribution* with location `μ` and scale `θ` has probability density function

$f(x; \mu, \theta) = \frac{1}{4 \theta} \mathrm{sech}^2
\left( \frac{x - \mu}{\theta} \right)$

```julia
Logistic()       # Logistic distribution with zero location and unit scale, i.e. Logistic(0.0, 1.0)
Logistic(u)      # Logistic distribution with location u and unit scale, i.e. Logistic(u, 1.0)
Logistic(u, b)   # Logistic distribution with location u ans scale b

params(d)       # Get the parameters, i.e. (u, b)
location(d)     # Get the location parameter, i.e. u
scale(d)        # Get the scale parameter, i.e. b
```

External links

* [Logistic distribution on Wikipedia](http://en.wikipedia.org/wiki/Logistic_distribution)

"""
immutable Logistic <: ContinuousUnivariateDistribution
    μ::Float64
    θ::Float64

    Logistic(μ::Real, θ::Real) = (@check_args(Logistic, θ > zero(θ)); new(μ, θ))
    Logistic(μ::Real) = new(μ, 1.0)
    Logistic() = new(0.0, 1.0)
end

@distr_support Logistic -Inf Inf

#### Parameters

location(d::Logistic) = d.μ
scale(d::Logistic) = d.θ

params(d::Logistic) = (d.μ, d.θ)

#### Statistics

mean(d::Logistic) = d.μ
median(d::Logistic) = d.μ
mode(d::Logistic) = d.μ

std(d::Logistic) = π * d.θ / sqrt3
var(d::Logistic) = (π * d.θ)^2 / 3.0
skewness(d::Logistic) = 0.0
kurtosis(d::Logistic) = 1.2

entropy(d::Logistic) = log(d.θ) + 2.0

#### Evaluation

zval(d::Logistic, x::Float64) = (x - d.μ) / d.θ
xval(d::Logistic, z::Float64) = d.μ + z * d.θ

pdf(d::Logistic, x::Float64) = (e = exp(-zval(d, x)); e / (d.θ * (1.0 + e)^2))
logpdf(d::Logistic, x::Float64) = (u = -abs(zval(d, x)); u - 2.0 * log1pexp(u) - log(d.θ))

cdf(d::Logistic, x::Float64) = logistic(zval(d, x))
ccdf(d::Logistic, x::Float64) = logistic(-zval(d, x))
logcdf(d::Logistic, x::Float64) = -log1pexp(-zval(d, x))
logccdf(d::Logistic, x::Float64) = -log1pexp(zval(d, x))

quantile(d::Logistic, p::Float64) = xval(d, logit(p))
cquantile(d::Logistic, p::Float64) = xval(d, -logit(p))
invlogcdf(d::Logistic, lp::Float64) = xval(d, -logexpm1(-lp))
invlogccdf(d::Logistic, lp::Float64) = xval(d, logexpm1(-lp))

function gradlogpdf(d::Logistic, x::Float64)
    e = exp(-zval(d, x))
    ((2.0 * e) / (1.0 + e) - 1.0) / d.θ
end

mgf(d::Logistic, t::Real) = exp(t * d.μ) / sinc(d.θ * t)

function cf(d::Logistic, t::Real)
    a = (π * t) * d.θ
    a == zero(a) ? complex(one(a)) : cis(t * d.μ) * (a / sinh(a))
end

#### Sampling

rand(d::Logistic) = quantile(d, rand())


text:
using Colors
using OpenCV
const cv = OpenCV

export  StereoMatcher, StereoBM,
        compute

"""
Base type for all stereo matching algorithms, or
stereo correspondence algorithms.
"""
abstract type StereoMatcher end

"""
Computing stereo correspondence using the block matching algorithm. This is a wrapper
arround the corresponding OpenCV algorithm.
"""
struct StereoBM <: StereoMatcher
   ptr::OpenCV.cv_PtrAllocated{OpenCV.StereoBM}

   function StereoBM(num_disparities::Integer, block_size::Integer)
       new(cv.StereoBM_create(
               numDisparities=Int32(num_disparities),
               blockSize=Int32(block_size)))
   end
end

"""
     StereoBM(;num_disparities=0, block_size=21)
Creates stereo matching algorithm based on block matching.

- `num_disparities` the disparity search range. For each pixel algorithm will find the best disparity
   from 0 (default minimum disparity) to `num_disparities`.
   The search range can then be shifted by changing the minimum disparity.

- `block_size` the linear size of the blocks compared by the algorithm.
   The size should be odd (as the block is centered at the current pixel).
   Larger block size implies smoother, though less accurate disparity map.
   Smaller block size gives more detailed disparity map,
   but there is higher chance for algorithm to find a wrong correspondence.
"""
function StereoBM(;num_disparities=0, block_size=21)
    StereoBM(num_disparities, block_size)
end

"""
    compute(matcher::StereoMatcher, left::Matrix{T}, right::Matrix{T}) -> right::Matrix{Gray}

Performs stereo matching between two images given by `left` and `right` where each element has to be
of a type `T` which is either a color such `Gray` or `RGB` or some kind of floating point of fixed point
type which can represent values between 0 and 1.

We get back a dispartity map, which is 2D matrix of `Gray` values.
"""
function compute(matcher::StereoMatcher, left::Matrix{T}, right::Matrix{T}) where T <: Union{Real, Color}
    imgL = image_to_opencv_input(left)
    imgR = image_to_opencv_input(right)

    disparity = cv.compute(matcher.ptr, imgL, imgR)

    opencv_disparity_map_to_image(disparity)
end


text:
#
# Do work here ...
#

s = "In Fortran We Trust"
println("Julia side: I will return a String: ", s)

return s

max_stars_repo_name: yundantianchang/dynload-julia

max_issues_repo_name: yundantianchang/dynload-julia

max_forks_repo_name: yundantianchang/dynload-julia


text:
# Autogenerated wrapper script for nghttp3_jll for x86_64-linux-musl
export libnghttp3

JLLWrappers.@generate_wrapper_header("nghttp3")
JLLWrappers.@declare_library_product(libnghttp3, "libnghttp3.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnghttp3,
        "lib/libnghttp3.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
data_populate(c = 0.5; half = "no", level = 3, nlmaxit=10)

This makes the files that plotnsold uses to make the figures
and tables in the paper.

If you set half="yes" be prepared to wait a couple weeks for the
files to come back.

I have generously precomputed all this stuff and put it in the
Data_From_Paper directory.

This function is part of the MPResults2019 module. Do not try to use it
without doing a

using MPResults

at the Julia prompt.
"""
function data_populate(c = 0.5; half = "no", level = 3, nlmaxit=10)
    #
    # Check the value of c against what the paper uses.
    #
    if c == 0.5
        locator = "=5"
    elseif c == 0.99
        locator = "=99"
    elseif c == 1
        locator = "=1"
    else
        println("Your choices are c=.5, .99, 1.")
        return
    end

    #
    # Get the directories organized. Make the base directory if
    # it's not there already.
    #
    if half == "no"
        BaseDirectory = string("Mixed_Precision_c", locator)
    else
        BaseDirectory = string("Mixed_Precision_c", locator, "/Float16")
    end
    try
        mkdir(BaseDirectory)
    catch
    end
    #
    # Get the file names.
    #
    savedat = filenames(BaseDirectory)
    #
    # Do the computations and put the data where it needs to be.
    #
    for nd = 1:level
        n = 512 * 2^nd

        if half == "no"
            fout64 = heqtest(n, c, "no"; jmaxit = nlmaxit)
            fout32 = heqtest(n, c, "no"; jprecision = Float32, jmaxit = nlmaxit)
            @save savedat[nd] fout64 fout32
        else
            fout16 = heqtest(n, c, "no";
                 jprecision = Float16, jmaxit = nlmaxit)
            @save savedat[nd] fout16
        end

    end

end

function filenames(BaseDirectory)
    savedat = Array{String,1}(undef, 5)
    for nd = 1:5
       # dimstring = string("/paper", string(512 * 2^nd), ".jld")
        dimstring=string("/paper",string(512*2^nd),".jld2")
        savedat[nd] = string(BaseDirectory, dimstring)
    end
    return savedat
end

max_stars_repo_name: ctkelley/MPResults


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "PThreadPool"
version = v"0.0.20200302"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/Maratyszcza/pthreadpool.git", "ebd50d0cfa3664d454ffdf246fcd228c3b370a11"),
    GitSource("https://github.com/Maratyszcza/FXdiv.git", "fd804a929fc64be9e40ee58bb51ed9f9cac98244"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pthreadpool
mkdir build
cd build
cmake \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DFXDIV_SOURCE_DIR=$WORKSPACE/srcdir/FXdiv \
    -DPTHREADPOOL_LIBRARY_TYPE=shared \
    -DPTHREADPOOL_BUILD_TESTS=OFF \
    -DPTHREADPOOL_BUILD_BENCHMARKS=OFF \
    ..
cmake --build . -- -j $nproc
make install
if [[ $target == *-w64-mingw32 ]]; then
    install -Dvm 755 libpthreadpool.dll "${libdir}/libpthreadpool.dll"
fi
install_license ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libpthreadpool", :libpthreadpool),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"5",
    julia_compat="1.6")


text:
module Day7

TESTINPUT = "16,1,2,0,4,2,7,1,2,14"

single_cost(pos::Int64, target::Int64) = abs(pos - target)
cost(positions::Vector{Int64}, target::Int64) = single_cost.(positions, target) |> sum
start(ps::Vector{Int64}) = (maximum(ps) - minimum(ps)) / 2 |> floor |> Int64

function gradient_descent(func, start::Int64)
    mid = func(start)
    pre = func(start - 1)
    post = func(start + 1)
    if (mid <= pre) && (mid <= post)
        return (start, mid)
    elseif pre < mid < post
        return gradient_descent(func, start - 1)
    elseif pre > mid > post
        return gradient_descent(func, start + 1)
    end
    throw("wat?")
end

run_part1(s) = joinpath(s, "day7.txt") |>
    f -> read(f, String) |>
    s -> parse.(Int64, split(s, ",")) |>
    ps -> (y -> cost(ps, y), start(ps)) |>
    x -> gradient_descent(x...) |>
    t -> "Optimal position is $(t[1]) with a cost of $(t[2])"

single_cost_2(pos::Int64, target::Int64) = (abs(pos - target) / 2) * (1 + abs(pos - target)) |> Int64
cost_2(positions::Vector{Int64}, target::Int64) = single_cost_2.(positions, target) |> sum

run_part2(s) = joinpath(s, "day7.txt") |>
    f -> read(f, String) |>
    s -> parse.(Int64, split(s, ",")) |>
    ps -> (y -> cost_2(ps, y), start(ps)) |>
    x -> gradient_descent(x...) |>
    t -> "Optimal position is $(t[1]) with a cost of $(t[2])"

end


text:
@doc raw"""
    PositiveNumbers <: Manifold{ℝ}

The hyperbolic manifold of positive numbers $H^1$ is a the hyperbolic manifold represented
by just positive numbers.

# Constructor

    PositiveNumbers()

Generate the `ℝ`-valued hyperbolic model represented by positive positive numbers.
To use this with arrays (1-element arrays),
please use [`SymmetricPositiveDefinite`](@ref)`(1)`.
"""
struct PositiveNumbers <: Manifold{ℝ} end

"""
    PositiveVectors(n)

Generate the manifold of vectors with positive entries.
This manifold is modeled as a [`PowerManifold`](@ref) of [`PositiveNumbers`](@ref).
"""
PositiveVectors(n::Integer) = PositiveNumbers()^n

"""
    PositiveMatrices(m,n)

Generate the manifold of matrices with positive entries.
This manifold is modeled as a [`PowerManifold`](@ref) of [`PositiveNumbers`](@ref).
"""
PositiveMatrices(n::Integer, m::Integer) = PositiveNumbers()^(n, m)

"""
    PositiveArrays(n₁,n₂,...,nᵢ)

Generate the manifold of `i`-dimensional arrays with positive entries.
This manifold is modeled as a [`PowerManifold`](@ref) of [`PositiveNumbers`](@ref).
"""
PositiveArrays(n::Vararg{Int,I}) where {I} = PositiveNumbers()^(n)

@doc raw"""
    check_manifold_point(M::PositiveNumbers, p)

Check whether `p` is a point on the [`PositiveNumbers`](@ref) `M`, i.e. $p>0$.
"""
function check_manifold_point(M::PositiveNumbers, p; kwargs...)
    if any(p .<= 0.0)
        return DomainError(
            p,
            "The point $(p) does not lie on $(M), since it is nonpositive.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M::PositiveNumbers, p, X; check_base_point, kwargs...)

Check whether `X` is a tangent vector in the tangent space of `p` on the
[`PositiveNumbers`](@ref) `M`.
For the real-valued case represented by positive numbers, all `X` are valid, since the tangent space is the whole real line.
For the complex-valued case `X` [...]

"""
function check_tangent_vector(M::PositiveNumbers, p, X; check_base_point = true, kwargs...)
    if check_base_point
        perr = check_manifold_point(M, p; kwargs...)
        return perr # if x is valid all v that are real numbers are valid
    end
    return nothing
end

@doc raw"""
    distance(M::PositiveNumbers, p, q)

Compute the distance on the [`PositiveNumbers`](@ref) `M`, which is

````math
d(p,q) = \Bigl\lvert \log \frac{p}{q} \Bigr\rvert = \lvert \log p - \log q\rvert.
````
"""
distance(::PositiveNumbers, ::Any, ::Any)

distance(::PositiveNumbers, p::Real, q::Real) = abs(log(p / q))
distance(::PositiveNumbers, p, q) = sum(abs.(log.(p ./ q)))

@doc raw"""
    exp(M::PositiveNumbers, p, X)

Compute the exponential map on the [`PositiveNumbers`](@ref) `M`.
```math
\exp_p X = p\operatorname{exp}(X/p).
```
"""
Base.exp(::PositiveNumbers, ::Any, ::Any)

Base.exp(::PositiveNumbers, p::Real, X::Real) = p * exp(X / p)

exp!(::PositiveNumbers, q, p, X) = (q .= p .* exp.(X ./ p))

flat(::PositiveNumbers, ::Real, X::TFVector) = FVector(CotangentSpace, X.data)

flat!(::PositiveNumbers, ξ::CoTFVector, p, X::TFVector) = copyto!(ξ, X)

@doc raw"""
    injectivity_radius(M::PositiveNumbers[, p])

Return the injectivity radius on the [`PositiveNumbers`](@ref) `M`, i.e. $\infty$.
"""
injectivity_radius(::PositiveNumbers) = Inf
injectivity_radius(::PositiveNumbers, ::ExponentialRetraction) = Inf
injectivity_radius(::PositiveNumbers, ::Any) = Inf
injectivity_radius(::PositiveNumbers, ::Any, ::ExponentialRetraction) = Inf
eval(
    quote
        @invoke_maker 1 Manifold injectivity_radius(
            M::PositiveNumbers,
            rm::AbstractRetractionMethod,
        )
    end,
)

@doc raw"""
    inner(M::PositiveNumbers, p, X, Y)

Compute the inner product of the two tangent vectors `X,Y` from the tangent plane at `p` on
the [`PositiveNumbers`](@ref) `M`, i.e.

````math
g_p(X,Y) = \frac{XY}{p^2}.
````
"""
inner(::PositiveNumbers, ::Any...)
@inline inner(::PositiveNumbers, p::Real, X::Real, Y::Real) = X * Y / p^2
@inline inner(::PositiveNumbers, p, X, Y) = sum((X .* Y) ./ (p .^ 2))
function inverse_retract(M::PositiveNumbers, x, y)
    return inverse_retract(M, x, y, LogarithmicInverseRetraction())
end
function inverse_retract(M::PositiveNumbers, x, y, ::LogarithmicInverseRetraction)
    return log(M, x, y)
end

@doc raw"""
    log(M::PositiveNumbers, p, q)

Compute the logarithmic map on the [`PositiveNumbers`](@ref) `M`.
````math
\log_p q = p\log\frac{q}{p}.
````
"""
Base.log(::PositiveNumbers, ::Any, ::Any)
Base.log(::PositiveNumbers, p::Real, q::Real) = p * log(q / p)

log!(::PositiveNumbers, X, p, q) = (X .= p .* log.(q ./ p))

@doc raw"""
    manifold_dimension(M::PositiveNumbers)

Return the dimension of the [`PositiveNumbers`](@ref) `M`,
i.e. of the 1-dimensional hyperbolic space,

````math
\dim(H^1) = 1
````
"""
manifold_dimension(::PositiveNumbers) = 1

mid_point(M::PositiveNumbers, p1, p2) = exp(M, p1, log(M, p1, p2) / 2)

@inline LinearAlgebra.norm(::PositiveNumbers, p, X) = sum(abs.(X ./ p))

@doc raw"""
    project(M::PositiveNumbers, p, X)

Project a value `X` onto the tangent space of the point `p` on the [`PositiveNumbers`](@ref) `M`,
which is just the identity.
"""
project(::PositiveNumbers, ::Any, ::Any)
project(::PositiveNumbers, ::Real, X::Real) = X

project!(::PositiveNumbers, Y, p, X) = (Y .= X)

retract(M::PositiveNumbers, p, q) = retract(M, p, q, ExponentialRetraction())
retract(M::PositiveNumbers, p, q, ::ExponentialRetraction) = exp(M, p, q)

representation_size(::PositiveNumbers) = ()

sharp(::PositiveNumbers, ::Real, ξ::CoTFVector) = FVector(TangentSpace, ξ.data)

sharp!(::PositiveNumbers, X::TFVector, p, ξ::CoTFVector) = copyto!(X, ξ)

Base.show(io::IO, ::PositiveNumbers) = print(io, "PositiveNumbers()")

function Base.show(
    io::IO,
    ::PowerManifold{ℝ,PositiveNumbers,TSize,ArrayPowerRepresentation},
) where {TSize}
    s = [TSize.parameters...]
    (length(s) == 1) && return print(io, "PositiveVectors($(s[1]))")
    (length(s) == 2) && return print(io, "PositiveMatrices($(s[1]), $(s[2]))")
    return print(io, "PositiveArrays($(join(s, ", ")))")
end

@doc raw"""
    vector_transport_to(M::PositiveNumbers, p, X, q, ::ParallelTransport)

Compute the parallel transport of `X` from the tangent space at `p` to the tangent space at
`q` on the [`PositiveNumbers`](@ref) `M`.

````math
\mathcal P_{q\gets p}(X) = X\cdot\frac{q}{p}.
````
"""
vector_transport_to(::PositiveNumbers, ::Any, ::Any, ::Any, ::ParallelTransport)
function vector_transport_to(
    ::PositiveNumbers,
    p::Real,
    X::Real,
    q::Real,
    ::ParallelTransport,
)
    return X * q / p
end

function vector_transport_to!(::PositiveNumbers, Y, p, X, q, ::ParallelTransport)
    return (Y .= X .* q ./ p)
end

function vector_transport_direction(
    M::PositiveNumbers,
    p,
    X,
    Y,
    m::AbstractVectorTransportMethod,
)
    q = exp(M, p, Y)
    return vector_transport_to(M, p, X, q, m)
end

zero_tangent_vector(::PositiveNumbers, p::Real) = zero(p)
zero_tangent_vector!(::PositiveNumbers, X, p) = fill!(X, 0)


text:
using StochasticDiffEq, DiffEqNoiseProcess, Test, Random
Random.seed!(100)
α=1
β=1
u₀=1/2
f(u,p,t) = α*u
g(u,p,t) = β*u
dt = 1//2^(4)
tspan = (0.0,1.0)
prob = SDEProblem(f,g,u₀,(0.0,1.0))
sol =solve(prob,EulerHeun(),dt=0.001,save_noise=true)
_sol = deepcopy(sol) # to make sure the plot is correct
W3 = NoiseGrid(reverse!(_sol.t),reverse!(_sol.W))
prob3 = SDEProblem(f,g,sol[end],(1.0,0.0),noise=W3)
sol2 = solve(prob3,EulerHeun(),dt=0.001)
@test sol.u ≈ reverse!(sol2.u) atol=1e-1


text:
# Converting cartesian coordinates to orbital elements

"""

A 3d point in cartesian space
"""
struct Point{T<:AbstractFloat}
    x::T
    y::T
    z::T
end

# Constructors
Point(x::AbstractVector) = Point(x...)
Point(x::AbstractMatrix) = [Point(x[:,i]) for i in eachindex(x[1,:])]
Point(x::Real) = Point(float(x),float(x),float(x))
unpack(x::Point) = (x.x,x.y,x.z)

# Overloads for products
LinearAlgebra.dot(x::Point,v::Point) = x.x*v.x + x.y*v.y + x.z*v.z
LinearAlgebra.dot(x::Point) = x.x*x.x + x.y*x.y + x.z*x.z
LinearAlgebra.cross(x::Point,v::Point) = Point(x.y*v.z - x.z*v.y,-(x.x*v.z - x.z*v.x),x.x*v.y - x.y*v.x)

""" Calculate the relative positions from the A-Matrix. """
function get_relative_positions(x,v,ic::InitialConditions)
    n = ic.nbody
    X = zeros(3,n)
    V = zeros(3,n)

    X .= permutedims(ic.amat*x')
    V .= permutedims(ic.amat*v')

    return Point(X), Point(V)
end

""" Get relative masses(*G) from initial conditions. """
function get_relative_masses(ic::InitialConditions)
    N = length(ic.m)
    M = zeros(N-1)
    G = 39.4845/(365.242 * 365.242) # AU^3 Msol^-1 Day^-2
    for i in 1:N-1
        for j in 1:N
            M[i] += abs(ic.ϵ[i,j])*ic.m[j]
        end
    end
    return G .* M
end

# Position and velocity magnitudes
mag(x) = sqrt(dot(x))
mag(x,v) = sqrt(dot(x,v))
Rdotmag(R,V,h) = sqrt(V^2 - (h/R)^2)

function hvec(r,rdot)
    hx,hy,hz = unpack(cross(r,rdot))
    hz >= 0.0 ? hy *= -1 : hx *= -1
    return Point(hx,hy,hz)
end

function calc_Ω(hx,hy,h,I)
    sinΩ = hx/(h*sin(I))
    cosΩ = hy/(h*sin(I))
    return atan(sinΩ,cosΩ)
end

function calc_ϖ(x,R,Rdot,I,Ω,a,e,h)
    # Find ω + f
    wpf = 0.0
    if I != 0.0
        sinwpf = x.z/(R*sin(I))
        coswpf = ((x.x/R) + sin(Ω)*sinwpf*cos(I))/cos(Ω)
        wpf = atan(sinwpf,coswpf)
    end
    # Find f
    sinf = a*Rdot*(1.0 - e^2)/(h*e)
    cosf = (a*(1.0-e^2)/R - 1.0)/e
    f = atan(sinf,cosf)

    # ϖ = Ω + ω
    return Ω + wpf - f
end

function convert_to_elements(x,v,M)
    R = mag(x)
    V = mag(v)
    h = mag(hvec(x,v))
    hx,hy,hz = unpack(hvec(x,v))
    Rdot = sign(dot(x,v))*Rdotmag(R,V,h)

    Gmm = M

    a = 1.0/((2.0/R) - (V*V)/(Gmm))
    e = sqrt(1.0 - (h*h/(Gmm*a)))
    I = acos(hz/h)

    # Make sure Ω is defined
    I != 0.0 ? Ω = calc_Ω(hx,hy,h,I) : Ω = 0.0

    ϖ = calc_ϖ(x,R,Rdot,I,Ω,a,e,h)
    P = (2.0*π)*sqrt(a*a*a/Gmm)

    return [P,0.0,e*cos(ϖ),e*sin(ϖ),I,Ω,a,e,ϖ]
end

function get_orbital_elements(s::State{T},ic::InitialConditions{T}) where T<:AbstractFloat
    elems = Elements{T}[]
    μ = get_relative_masses(ic)
    X,V = get_relative_positions(s.x,s.v,ic)
    N = ic.nbody

    push!(elems,Elements(m=ic.m[1]))

    i = 1; b = 0
    while i < N

        # Check if new binary
        if first(ic.ϵ[i,:]) == zero(T)
            b+=1
        end

        new_elems = convert_to_elements(X[i+b],V[i+b],μ[i+b])
        push!(elems,Elements(ic.m[i+1],new_elems...))

        # Compensate for new binary in step
        if b > 0
            b -= 2
        elseif b < 0
            i += 1
        end

    i+=1
    end
    return elems
end


text:
"""
    metabolite_fluxes(model::MetabolicModel, flux_dict::Dict{String, Float64})

Return two dictionaries of metabolite `id`s mapped to reactions that consume or
produce them, given the flux distribution supplied in `flux_dict`.
"""
function metabolite_fluxes(model::MetabolicModel, flux_dict::Dict{String,Float64})
    S = stoichiometry(model)
    rids = reactions(model)
    mids = metabolites(model)

    producing = Dict{String,Dict{String,Float64}}()
    consuming = Dict{String,Dict{String,Float64}}()
    for (row, mid) in enumerate(mids)
        for (col, rid) in enumerate(rids)
            mf = flux_dict[rid] * S[row, col]
            if mf < 0 # consuming rxn
                if haskey(consuming, mid)
                    consuming[mid][rid] = mf
                else
                    consuming[mid] = Dict(rid => mf)
                end
            elseif mf >= 0 # producing rxn
                if haskey(producing, mid)
                    producing[mid][rid] = mf
                else
                    producing[mid] = Dict(rid => mf)
                end
            end
        end
    end
    return consuming, producing
end

"""
    atom_fluxes(model::MetabolicModel, reaction_fluxes::Dict{String, Float64})

Return a dictionary mapping the flux of atoms across a flux solution given by
`reaction_fluxes` using the reactions in `model` to determine the appropriate stoichiometry.

Note, this function ignores metabolites with no formula assigned to them, no error message
will be displayed.

Note, if a model is mass balanced there should be not net flux of any atom. By removing
reactions from the flux_solution you can investigate how that impacts the mass balances.

# Example
```
# Find flux of Carbon through all metabolic reactions except the biomass reaction
delete!(fluxes, "BIOMASS_Ecoli_core_w_GAM")
atom_fluxes(model, fluxes)["C"]
```
"""
function atom_fluxes(model::MetabolicModel, reaction_fluxes::Dict{String,Float64})
    rids = reactions(model)
    atom_flux = Dict{String,Float64}()
    for (ridx, rid) in enumerate(rids)
        haskey(reaction_fluxes, rid) || continue
        rflux = reaction_fluxes[rid]
        for (mid, mstoi) in reaction_stoichiometry(model, rid)
            atoms = metabolite_formula(model, mid)
            isnothing(atoms) && continue # missing formulas are ignored
            for (atom, abundance) in atoms
                atom_flux[atom] = get(atom_flux, atom, 0.0) + rflux * mstoi * abundance
            end
        end
    end
    return atom_flux
end


text:
using LinearAlgebra;

θ, dnn = Flux.destructure(Chain(Dense(numInvariants, hiddenDim, tanh),
                                Dense(hiddenDim, numSymBasis)));

function lossAndManualGrads(x,
                            y,
                            θ,
                            tb,
                            hiddenDim,
                            numHiddenLayers,
                            σ)
    @assert (1==1)#length(x) == length(y) == length(tb)
    W = [zeros(
    W = fill(zeros(hiddenDim, hiddenDim), (numHiddenLayers));
    b = fill(zeros(hiddenDim), (numHiddenLayers))

end

max_stars_repo_name: cmhyett/MLCoarseGrainedVGT

max_issues_repo_name: cmhyett/MLCoarseGrainedVGT

max_forks_repo_name: cmhyett/MLCoarseGrainedVGT