text:
# Helper functions for methods based on inner-outer iterations

using LinearAlgebra

export inner_solve;
export InnerSolver;
export NewtonInnerSolver
export PolyeigInnerSolver
export DefaultInnerSolver
export IARInnerSolver
export IARChebInnerSolver
export SGIterInnerSolver
export ContourBeynInnerSolver
export NleigsInnerSolver

export compute_rf

"""
    abstract type InnerSolver

Structs inheriting from this type are used to solve inner problems
in an inner-outer iteration.

The `InnerSolver` objects are passed to the NEP-algorithms,
which uses it to dispatch the correct version of the function [`inner_solve`](@ref).
Utilizes existing implementations of NEP-solvers and [`inner_solve`](@ref) acts as a wrapper
to these.

# Example

There is a [`DefaultInnerSolver`](@ref) that dispatches an inner solver based on the
provided NEP.
However, this example shows how you can force [`nlar`](@ref) to use the
[`IARInnerSolver`](@ref) for a PEP.
```julia-repl
julia> nep=nep_gallery("pep0", 100);
julia> λ,v = nlar(nep, inner_solver_method=NEPSolver.IARInnerSolver(), neigs=1, num_restart_ritz_vecs=1, maxit=70, tol=1e-8);
julia> norm(compute_Mlincomb(nep,λ[1],vec(v)))
8.68118417430353e-9
```

See also: [`inner_solve`](@ref), [`DefaultInnerSolver`](@ref), [`NewtonInnerSolver`](@ref),
[`PolyeigInnerSolver`](@ref), [`IARInnerSolver`](@ref), [`IARChebInnerSolver`](@ref),
[`SGIterInnerSolver`](@ref), [`ContourBeynInnerSolver`](@ref)
"""
abstract type InnerSolver end;

"""
    struct DefaultInnerSolver <: InnerSolver

Dispatches a version of [`inner_solve`](@ref) based on the
type of the NEP provided. This function tries to automatically
detect which solver is recommended.

See also: [`InnerSolver`](@ref), [`inner_solve`](@ref)
"""
struct DefaultInnerSolver <: InnerSolver end;

"""
    struct NewtonInnerSolver
    function NewtonInnerSolver(;tol=1e-13,maxit=80,starting_vector=:Vk,
                               newton_function=augnewton)

Uses a Newton-like method to solve the inner problem, with tolerance,
and maximum number of iterations given by
`tol` and `maxit`. The starting vector can be `:ones`,
`:randn`, or `:Vk`. The value `:Vk` specifies the use
of the outer NEP-solver keyword argument (`Vk`).
This is typically the previous iterate in the outer method.

The kwarg `newton_function`, specifies a `Function`
which is called. The type supports [`augnewton`](@ref),
[`newton`](@ref), [`resinv`](@ref)
[`quasinewton`](@ref), [`newtonqr`](@ref).
In principle it can be any
method which takes the same keyword arguments as these methods.

See also: [`InnerSolver`](@ref), [`inner_solve`](@ref)
"""
struct NewtonInnerSolver <: InnerSolver
    tol::Float64
    maxit::Int
    starting_vector::Symbol;
    newton_function::Function;
    function NewtonInnerSolver(;tol=1e-13,maxit=80,starting_vector=:Vk,
                               newton_function=augnewton)
        return new(tol,maxit,starting_vector,newton_function);
    end
end;

"""
    struct PolyeigInnerSolver <: InnerSolver
    function PolyeigInnerSolver()

Specifies the use of [`polyeig`](@ref) to solve the inner problem.
This is intended for NEPs of the type [`PEP`](@ref).

See also: [`InnerSolver`](@ref), [`inner_solve`](@ref)
"""
struct PolyeigInnerSolver <: InnerSolver end;

"""
    struct IARInnerSolver
    function IARInnerSolver(;tol=1e-13,maxit=80,
               starting_vector=:ones,normalize_DEPs=:auto,
               iar_function=iar)

Uses [`iar`](@ref), [`tiar`](@ref) or [`iar_chebyshev`](@ref)
to solve the inner problem, with tolerance,
and maximum number of iterations given by
`tol` and `maxit`. The starting vector can be `:ones` or
`:randn`.
The `iar_function` can be `iar`, `tiar` or `iar_chebyshev`
(or any function taking the same parameters as input).
`normalize_DEPs` determines if the we should carry out
precomputation and make sure the projection of a [`DEP`](@ref)
 is again a `DEP` (which can speed up performance).
It can take the value `true`, `false` or `:auto`. `:auto`
sets it to true if we use the `iar_chebyshev` solver.

The kwarg `iar_function`, specifies a `Function`
which is called. Examples of functions are `iar`
and `iar_chebyshev`. It can be any NEP-solver
which takes the same keyword arguments as these methods.

See also: [`InnerSolver`](@ref), [`inner_solve`](@ref)
"""
struct IARInnerSolver <: InnerSolver
    tol::Float64
    maxit::Int
    starting_vector::Symbol;
    normalize_DEPs::Bool;
    iar_function::Function;
    function IARInnerSolver(;tol=1e-13,maxit=80,starting_vector=:ones,
                            normalize_DEPs=:auto,iar_function=iar)
        # Default starting_vector = :ones since we would otherwise
        # get a place with difficult reproducability and hidden.

        if (normalize_DEPs == :auto)
            normalize_DEPs = (iar_function == iar_chebyshev)
        end

        return new(tol,maxit,starting_vector,normalize_DEPs,iar_function);
    end
end;

"""
    function IARChebInnerSolver(;tol=1e-13,maxit=80,starting_vector=:ones,
                                normalize_DEPs=true)

Uses [`iar_chebyshev`](@ref) to solve the inner problem.
See [`IARInnerSolver`](@ref) for keyword argument documentation.

See also: [`IARInnerSolver`](@ref), [`InnerSolver`](@ref), [`inner_solve`](@ref)
"""
function IARChebInnerSolver(;tol=1e-13,maxit=80,starting_vector=:ones,
                            normalize_DEPs=true)
    return IARInnerSolver(tol=tol,maxit=maxit,starting_vector=starting_vector,
                          normalize_DEPs=normalize_DEPs,
                          iar_function=iar_chebyshev)
end

"""
    struct SGIterInnerSolver <: InnerSolver

Uses [`sgiter`](@ref) to solve the inner problem.

See also: [`InnerSolver`](@ref), [`inner_solve`](@ref)
"""
struct SGIterInnerSolver <: InnerSolver end;

"""
    struct ContourBeynInnerSolver <: InnerSolver
    function ContourBeynInnerSolver(;tol=sqrt(eps(real(Float64))),
                                    radius=:auto,N=1000)

Uses [`contour_beyn`](@ref) to solve the inner problem, with radius and number
of quadrature nodes, given by `radius` and `n`. If the variable `radius` is set
to `:auto`, the integration radius will be automatically selected by
using the eigenvalue approximations specified by the outer solver.

See also: [`InnerSolver`](@ref), [`inner_solve`](@ref)
"""
struct ContourBeynInnerSolver <: InnerSolver
    tol::Float64
    radius::Union{Real,Tuple,Array,Symbol} # integration radius
    N::Integer;  # Nof quadrature nodes
    function ContourBeynInnerSolver(;tol=sqrt(eps(real(Float64))),radius= :auto,N=1000)
        return new(tol,radius,N);
    end
end;

"""
    struct NleigsInnerSolver <: InnerSolver
    function NleigsInnerSolver(;Σ= :auto,nodes =:auto, tol=1e-6 )

Uses [`nleigs`](@ref) to solve the inner problem, in the region `Σ` with shifts
`nodes` and with tolerance `tol`. If the variable `Σ` is set to `:auto`, the
region `Σ` will be set by using the eigenvalues approximations.
See [`nleigs`](@ref) for description of parameters.

See also: [`InnerSolver`](@ref), [`inner_solve`](@ref)
"""
struct NleigsInnerSolver <: InnerSolver
    Σ::Union{Vector,Symbol} #
    nodes::Union{Vector,Symbol}
    tol::Float64;
    function NleigsInnerSolver(;Σ= :auto,nodes =:auto, tol=1e-6 )
        return new(Σ,nodes,tol);
    end
end;

"""
    inner_solve(is::InnerSolver,T_arit,nep;kwargs...)

Solves the projected linear problem with solver specied with `is`. This is to be used
as an inner solver in an inner-outer iteration. T specifies which method
to use. The most common choice is [`DefaultInnerSolver`](@ref). The function returns
`(λv,V)` where `λv` is an array of eigenvalues and `V` a matrix with corresponding
vectors.
The struct `T_arit` defines the arithmetics used in the outer iteration and should prefereably
also be used in the inner iteration.

Different inner_solve methods take different kwargs.
These are standardized kwargs:
- `neigs`: Number of wanted eigenvalues (but less or more may be returned)
- `σ`: target specifying where eigenvalues
- `λv`, `V`: Vector/matrix of guesses to be used as starting values
- `j`: the jth eigenvalue in a min-max characterization
- `tol`: Termination tolarance for inner solver
- `inner_logger`: Determines how the inner solves are logged. See [`Logger`](@ref) for further references
"""
function inner_solve(is::DefaultInnerSolver,T_arit::Type,nep::NEPTypes.Proj_NEP;kwargs...)
    if isa(nep.orgnep,NEPTypes.PEP)
        return inner_solve(PolyeigInnerSolver(),T_arit,nep;kwargs...);
    elseif isa(nep.orgnep,NEPTypes.DEP)
        # Should be Cheb IAR
        return inner_solve(IARChebInnerSolver(),T_arit,nep;kwargs...);
    elseif isa(nep.orgnep,NEPTypes.SPMF_NEP) # Default to IAR for SPMF
        return inner_solve(IARInnerSolver(),T_arit,nep;kwargs...);
    else
        return inner_solve(NewtonInnerSolver(),T_arit,nep;kwargs...);
    end
end

function inner_solve(is::NewtonInnerSolver,T_arit::Type,nep::NEPTypes.Proj_NEP;
                     λv=zeros(T_arit,1),
                     V=Matrix{T_arit}(rand(size(nep,1),size(λv,1))),
                     tol=sqrt(eps(real(T_arit))),
                     inner_logger=0,
                     kwargs...)
    @parse_logger_param!(inner_logger)
    for k=1:size(λv,1)
        try
            if (is.starting_vector == :ones)
                v0=ones(size(nep,1));
            elseif (is.starting_vector == :randn)
                v0=randn(size(nep,1));
            elseif (is.starting_vector == :Vk)
                v0=V[:,k]; # Starting vector for projected problem
            end

            n=size(nep,1);

            projerrmeasure=ResidualErrmeasure(nep);
            # Compute a solution to projected problem with Newton's method
            λ1,vproj=is.newton_function(
                T_arit,nep,logger=inner_logger,λ=λv[k],
                v=v0,maxit=is.maxit,tol=is.tol,
                errmeasure=projerrmeasure);
            V[:,k]=vproj;
            λv[k]=λ1;
        catch e
            if (isa(e, NoConvergenceException))
                V[:,k] = e.v
                λv[k] = e.λ;
            else
                rethrow(e)
            end
        end
    end
    return λv,V
end

function inner_solve(is::PolyeigInnerSolver,T_arit::Type,nep::NEPTypes.Proj_NEP;kwargs...)
    if !isa(nep.orgnep,NEPTypes.PEP)
        error("Wrong type. PolyeigInnerSolver only handles the PEP type.");
    end
    pep=NEPTypes.PEP(NEPTypes.get_Av(nep.nep_proj))
    return polyeig(T_arit,pep);
end

function inner_solve(is::IARInnerSolver,T_arit::Type,nep::NEPTypes.Proj_NEP;σ=0,neigs=10,inner_logger=0,kwargs...)
    @parse_logger_param!(inner_logger)
    try
        if (isa(nep.orgnep, NEPTypes.DEP) && is.normalize_DEPs)
            # normalize DEPs: Make projected the λ-term in the DEP normalized
            AA = get_Av(nep)
            TT = eltype(AA);
            if (TT  <: SubArray)
                # If it's better to transform to store in Matrix instead
                TT=Matrix{eltype(AA[1])};
            end
            BB = Vector{TT}(undef, size(AA,1)-1)
            for i = 1:(size(AA,1)-1)
                BB[i] = AA[1]\AA[1+i]
            end
            nep = DEP(BB,nep.orgnep.tauv)
        end

        if (is.starting_vector == :ones)
            v0=ones(size(nep,1));
        else
            v0=randn(size(nep,1));
        end
        λ,V=is.iar_function(
            T_arit,nep,σ=σ,neigs=neigs,tol=is.tol,
            maxit=is.maxit,logger=inner_logger,v=v0);
        return λ,V
    catch e
        if (isa(e, NoConvergenceException))
            # If we fail to find the wanted number of eigvals, we can still
            # return the ones we found
            return e.λ,e.v
        else
            rethrow(e)
        end
    end
end

function inner_solve(is::SGIterInnerSolver,T_arit::Type,nep::NEPTypes.Proj_NEP;λv=[0],j=0,inner_logger=0,kwargs...)
    @parse_logger_param!(inner_logger)
    λ,V=sgiter(T_arit,nep,j,logger=inner_logger)
    return [λ],reshape(V,size(V,1),1);
end

function inner_solve(is::ContourBeynInnerSolver,T_arit::Type,nep::NEPTypes.Proj_NEP;σ=0,λv=[0,1],neigs=10,inner_logger=0,kwargs...)
    @parse_logger_param!(inner_logger)
    # Radius computed as the largest distance σ and λv and a litte more
    if is.radius == :auto
        radius = maximum(abs.(σ .- λv))*1.5
    else
        radius = is.radius
    end
    neigs = Int(min(neigs,size(nep,1)-1))
    λ,V = contour_beyn(T_arit,nep,neigs=neigs,σ=σ,radius=radius,N=is.N,logger=inner_logger,tol=is.tol,sanity_check=false)
    return λ,V
end

function inner_solve(is::NleigsInnerSolver,T_arit::Type,nep::NEPTypes.Proj_NEP;λv=[0,1],inner_logger=0,kwargs...)
    @parse_logger_param!(inner_logger)

    # Σ computed as tje smallest disk containing the estimated eigenvalues λv and a bit more
    if is.Σ == :auto
        σ=sum(λv)/length(λv)
        r=maximum(abs.(σ.-λv))*1.5
        θ=range(0,stop=2π,length=1000)
        Σ=σ.+r*cos.(θ) + 1im*r*sin.(θ)
    else
        Σ = is.Σ
    end

    if is.nodes == :auto
        nodes = [zero(Complex{Float64})]
    else
        nodes = is.nodes
    end
    tol = is.tol
    λ,V = nleigs(nep,Σ;nodes=nodes,tol=tol,static=true)
    return λ,V
end


text:
# include("plot/shader.jl")
# include("plot/camera.jl")
# include("plot/plot.jl")

function plot(boxset::BoxSet{<:AbstractBoxPartition{<:Box{N}}}; kwargs...) where N
    m = HyperRectangle(GeometryBasics.Vec3f0(0), GeometryBasics.Vec3f0(1))
    c = [box.center .- box.radius for box in boxset]
    r = [1.9 .* box.radius for box in boxset]
    fig, ax, ms = meshscatter(GeometryBasics.Vec{N, Float32}.(c), marker = m, markersize = r,
                            color =:red; kwargs...)
end

function plot(boxfun::BoxFun{<:AbstractBoxPartition{<:Box{1}}}; kwargs...)
    c, r, v = Float32[], Float32[], Float32[]

    for (key, value) in boxfun.dict
        box = key_to_box(boxfun.partition, key)
        push!(c, box.center[1])
        push!(r, 1.9 * box.radius[1])
        push!(v, value)
    end
    fig, ax, bp = barplot(c, v, width = r; kwargs ...)
end

function plot(boxfun::BoxFun{<:AbstractBoxPartition{<:Box{3}}}; kwargs...) where N
    center, radius, color = GeometryBasics.Vec{3, Float32}[], GeometryBasics.Vec{3, Float32}[], Float32[]

    for (key, value) in boxfun.dict
        box = GAIO.key_to_box(boxfun.partition, key)
        push!(center, box.center)
        push!(radius, 1.9 .* box.radius)
        push!(color, value)
    end
    m = HyperRectangle(GeometryBasics.Vec3f0(0), GeometryBasics.Vec3f0(1))
    fig = Figure()
    fig, ax, ms = meshscatter(center, marker = m, markersize = radius,
                              color = color, colormap =:jet; kwargs...)
    Colorbar(fig[1,2], ms)
    fig
end

ext: jl


text:
using Convex, SCS
set_default_solver(SCSSolver(verbose = 0, max_iters = 1000000))

println(
    "small test to compile and see if we throw warnings when warmstarting what can't be warmstarted.",
)
n = 10
y = rand(n)
x = Variable(n)

lambda = [100]
problem = minimize(sumsquares(y - x) + lambda[1] * sumsquares(x - 10))
@time solve!(problem, warmstart = true)

println("\nnow try first solving and then warmstarting")
n = 1000
y = rand(n)
x = Variable(n)

lambda = [100]
problem = minimize(sumsquares(y - x) + lambda[1] * sumsquares(x - 10))
@time solve!(problem)
@show problem.optval
lambda[1] = 105
println(
    "this run should be faster than the last and have slightly higher optimal value",
)
@time solve!(problem, warmstart = true)
@show problem.optval

println("\nnow try with scalar variables")
lambda = 100
problem = minimize(sumsquares(y - x) + lambda * sumsquares(x - 10))
@time solve!(problem)
@show problem.optval
lambda = 105
println(
    "this run should be faster than the last and have slightly higher optimal value",
)
@time solve!(problem, warmstart = true)
@show problem.optval


text:
let
    fig,  ax = subplots(1,1)
    ax.plot(1:length(trace_log_p), trace_log_p)
    ax.set_yscale("log")
    ax.set_xscale("log")
    ax.grid(true)

end
let
    fig,  ax = subplots(1,1)
    -diff(trace_log_p) |> x -> ax.plot(1:length(x), x)
    ax.set_yscale("log")
    ax.set_xscale("log")
    ax.grid(true)
    ax.set_title("diff p")
end

plot_trace_value(p, ii) = begin
    fig,  ax = subplots(1,1)
    ax.plot(trace_p_cat[:k], trace_p_cat[p][ii,:])
    ax.plot(trace_p_cat[:k], p_true[p][ii]*ones(length(trace_p_cat[:k])), "-r")
    ax.grid(true)
    ax.set_xscale("log")
    ax.set_title("Trace $(p) index $(ii)")
end

plot_error(p, ii) = begin
    fig,  ax = subplots(1,1)
    ax.plot(trace_p_cat[:k], trace_p_cat[p][ii,:] .- p_true[p][ii] .|> abs)
    ax.grid(true)
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_title("Error $(p) index $(ii)")
end

parameters_2_plot =[
    (:T, CartesianIndex(1,1,1)),
    (:T, CartesianIndex(1,1,9)),
    (:T, CartesianIndex(1,2,9)),
    (:r, CartesianIndex(2,2)),
]
for (p, ii) in parameters_2_plot
    plot_error(p,ii)
    plot_trace_value(p, ii)
end


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

mutable struct ExpressRoutePeeringState <: SwaggerModel

    function ExpressRoutePeeringState(;)
        o = new()
        o
    end
end # type ExpressRoutePeeringState

const _property_map_ExpressRoutePeeringState = Dict{Symbol,Symbol}()
const _property_types_ExpressRoutePeeringState = Dict{Symbol,String}()
Base.propertynames(::Type{ ExpressRoutePeeringState }) = collect(keys(_property_map_ExpressRoutePeeringState))
Swagger.property_type(::Type{ ExpressRoutePeeringState }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ExpressRoutePeeringState[name]))}
Swagger.field_name(::Type{ ExpressRoutePeeringState }, property_name::Symbol) =  _property_map_ExpressRoutePeeringState[property_name]

function check_required(o::ExpressRoutePeeringState)
    true
end

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


text:
import DCELs._connect_vertices!, DCELs._split_face!
import DCELs._connect_manual!, DCELs._connect_automatic!

@testset "Connecting - primitive" begin
    u, v, w = Vertex(1, 0), Vertex(0, 1), Vertex(-1, 0)
    a1, a2, b1, b2, c1, c2 = [Edge() for _ in 1:6]
    f1, f2 = Face(), Face()

    u.edge, v.edge, w.edge = a1, b1, c1

    a1.dest, b1.dest, c1.dest = u, v, w
    a2.dest, b2.dest, c2.dest = w, u, v

    a1.twin, a2.twin = a2, a1
    b1.twin, b2.twin = b2, b1
    c1.twin, c2.twin = c2, c1

    a1.next, b1.next, c1.next = b1, c1, a1
    a1.prev, b1.prev, c1.prev = c1, a1, b1

    a2.next, b2.next, c2.next = c2, a2, b2
    a2.prev, b2.prev, c2.prev = b2, c2, a2

    a1.face = b1.face = c1.face = f1
    a2.face = b2.face = c2.face = f2

    f1.edge, f2.edge = a1, a2

    z = Vertex(2, 0)

    # No splitting
    h1 = _connect_vertices!(u, z)
    h2 = h1.twin

    @test h1.next == h2
    @test h1.prev == b2
    @test h2.next == a2
    @test h2.prev == h1

    @test h1.face == h2.face == f2

    # Splitting
    h3, f_del = _split_face!(z, v)
    h4 = h3.twin

    @test h1.next == h3
    @test h1.prev == b2
    @test h2.next == a2
    @test h2.prev == h4
    @test h3.next == b2
    @test h3.prev == h1
    @test h4.next == h2
    @test h4.prev == c2

    @test h1.face == h3.face == b2.face
    @test h4.face == h2.face == a2.face == c2.face

    @test f_del == f2
end

@testset "Connecting - manual" begin
    tess = Tessellation()

    u, v, w = Vertex(1, 0), Vertex(0, 1), Vertex(-1, 0)
    a1, a2, b1, b2, c1, c2 = [Edge() for _ in 1:6]
    f1, f2 = Face(), Face()

    push!(tess.vertices, u, v, w)
    push!(tess.edges, a1, a2, b1, b2, c1, c2)
    push!(tess.faces, f1, f2)

    u.edge, v.edge, w.edge = a1, b1, c1

    a1.dest, b1.dest, c1.dest = u, v, w
    a2.dest, b2.dest, c2.dest = w, u, v

    a1.twin, a2.twin = a2, a1
    b1.twin, b2.twin = b2, b1
    c1.twin, c2.twin = c2, c1

    a1.next, b1.next, c1.next = b1, c1, a1
    a1.prev, b1.prev, c1.prev = c1, a1, b1

    a2.next, b2.next, c2.next = c2, a2, b2
    a2.prev, b2.prev, c2.prev = b2, c2, a2

    a1.face = b1.face = c1.face = f1
    a2.face = b2.face = c2.face = f2

    f1.edge, f2.edge = a1, a2

    z = Vertex(2, 0)

    # No splitting
    h1 = _connect_manual!(tess, u, z, false)
    h2 = h1.twin

    @test h1.next == h2
    @test h1.prev == b2
    @test h2.next == a2
    @test h2.prev == h1

    @test h1.face == h2.face == f2

    # Splitting
    h3 = _connect_manual!(tess, z, v, true)
    h4 = h3.twin

    @test h1.next == h3
    @test h1.prev == b2
    @test h2.next == a2
    @test h2.prev == h4
    @test h3.next == b2
    @test h3.prev == h1
    @test h4.next == h2
    @test h4.prev == c2

    @test h1.face == h3.face == b2.face
    @test h4.face == h2.face == a2.face == c2.face
end

@testset "Connecting - automatic" begin
    tess = Tessellation()

    u, v, w = Vertex(1, 0), Vertex(0, 1), Vertex(-1, 0)
    a1, a2, b1, b2, c1, c2 = [Edge() for _ in 1:6]
    f1, f2 = Face(), Face()

    push!(tess.vertices, u, v, w)
    push!(tess.edges, a1, a2, b1, b2, c1, c2)
    push!(tess.faces, f1, f2)

    u.edge, v.edge, w.edge = a1, b1, c1

    a1.dest, b1.dest, c1.dest = u, v, w
    a2.dest, b2.dest, c2.dest = w, u, v

    a1.twin, a2.twin = a2, a1
    b1.twin, b2.twin = b2, b1
    c1.twin, c2.twin = c2, c1

    a1.next, b1.next, c1.next = b1, c1, a1
    a1.prev, b1.prev, c1.prev = c1, a1, b1

    a2.next, b2.next, c2.next = c2, a2, b2
    a2.prev, b2.prev, c2.prev = b2, c2, a2

    a1.face = b1.face = c1.face = f1
    a2.face = b2.face = c2.face = f2

    f1.edge, f2.edge = a1, a2

    z = Vertex(2, 0)

    # No splitting
    h1 = _connect_automatic!(tess, u, z)
    h2 = h1.twin

    @test h1.next == h2
    @test h1.prev == b2
    @test h2.next == a2
    @test h2.prev == h1

    @test h1.face == h2.face == f2

    # Splitting
    h3 = _connect_automatic!(tess, z, v)
    h4 = h3.twin

    @test h1.next == h3
    @test h1.prev == b2
    @test h2.next == a2
    @test h2.prev == h4
    @test h3.next == b2
    @test h3.prev == h1
    @test h4.next == h2
    @test h4.prev == c2

    @test h1.face == h3.face == b2.face
    @test h4.face == h2.face == a2.face == c2.face
end

@testset "Connecting - final" begin
    tess = Tessellation()

    u, v, w = Vertex(1, 0), Vertex(0, 1), Vertex(-1, 0)
    a1, a2, b1, b2, c1, c2 = [Edge() for _ in 1:6]
    f1, f2 = Face(), Face()

    push!(tess.vertices, u, v, w)
    push!(tess.edges, a1, a2, b1, b2, c1, c2)
    push!(tess.faces, f1, f2)

    u.edge, v.edge, w.edge = a1, b1, c1

    a1.dest, b1.dest, c1.dest = u, v, w
    a2.dest, b2.dest, c2.dest = w, u, v

    a1.twin, a2.twin = a2, a1
    b1.twin, b2.twin = b2, b1
    c1.twin, c2.twin = c2, c1

    a1.next, b1.next, c1.next = b1, c1, a1
    a1.prev, b1.prev, c1.prev = c1, a1, b1

    a2.next, b2.next, c2.next = c2, a2, b2
    a2.prev, b2.prev, c2.prev = b2, c2, a2

    a1.face = b1.face = c1.face = f1
    a2.face = b2.face = c2.face = f2

    f1.edge, f2.edge = a1, a2

    z = Vertex(2, 0)

    # No splitting
    h1 = connect(tess, u, z)
    h2 = h1.twin

    @test h1.next == h2
    @test h1.prev == b2
    @test h2.next == a2
    @test h2.prev == h1

    @test h1.face == h2.face == f2

    # Splitting
    h3 = connect(tess, z, v, split=true)
    h4 = h3.twin

    @test h1.next == h3
    @test h1.prev == b2
    @test h2.next == a2
    @test h2.prev == h4
    @test h3.next == b2
    @test h3.prev == h1
    @test h4.next == h2
    @test h4.prev == c2

    @test h1.face == h3.face == b2.face
    @test h4.face == h2.face == a2.face == c2.face
end


text:
push!( LOAD_PATH, "../src/" )

using Documenter, Touchstone

makedocs(
  format = :html,
  sitename = "Touchstone.jl",
  modules = [
    Touchstone
  ],
  pages = [
    "Home" => "index.md",
    "Library" => [
      "Public"    => "public.md",
      "Internals" => "internals.md"
    ]
  ]
)

deploydocs(
  repo   = "github.com/mpichl87/Touchstone.jl.git",
  target = "build",
  julia = "1.0",
  deps = nothing,
  make = nothing,
)


text:
using ITensors,
      ITensorsGPU,
      LinearAlgebra, # For tr()
      Combinatorics, # For permutations()
      CuArrays,
      Test

      # gpu tests!
@testset "cuITensor, Dense{$SType} storage" for SType ∈ (Float64, ComplexF64)
  mi,mj,mk,ml,ma = 2,3,4,5,6,7
  i = Index(mi,"i")
  j = Index(mj,"j")
  k = Index(mk,"k")
  l = Index(ml,"l")
  a = Index(ma,"a")
  @testset "Test add CuDense" begin
    A  = [SType(1.0) for ii in 1:dim(i), jj in 1:dim(j)]
    dA = ITensorsGPU.CuDense{SType, CuVector{SType}}(SType(1.0), dim(i)*dim(j))
    B  = [SType(2.0) for ii in 1:dim(i), jj in 1:dim(j)]
    dB = ITensorsGPU.CuDense{SType, CuVector{SType}}(SType(2.0), dim(i)*dim(j))
    dC = +(dA, IndexSet(i, j), dB, IndexSet(j, i))
    hC = collect(dC)
    @test collect(A + B) ≈ hC
  end
  @testset "Test subtract CuDense" begin
    A  = [SType(1.0) for ii in 1:dim(i), jj in 1:dim(j)]
    dA = ITensorsGPU.CuDense{SType, CuVector{SType}}(SType(1.0), dim(i)*dim(j))
    B  = [SType(2.0) for ii in 1:dim(i), jj in 1:dim(j)]
    dB = ITensorsGPU.CuDense{SType, CuVector{SType}}(SType(2.0), dim(i)*dim(j))
    dC = -(dA, IndexSet(i, j), dB, IndexSet(i, j))
    hC = collect(dC)
    @test A - B ≈ hC
  end
  @testset "Test permute CuDense" begin
    A  = [SType(ii*jj) for ii in 1:dim(i), jj in 1:dim(j)]
    dA = ITensorsGPU.CuDense{SType, CuVector{SType}}(NDTensors.Dense(vec(A)))
    B  = [SType(0.0) for ii in 1:dim(j), jj in 1:dim(j)]
    dB = ITensorsGPU.CuDense{SType, CuVector{SType}}(SType(0.0), dim(i)*dim(j))
    dC = permute!(dB, IndexSet(j, i), dA, IndexSet(i, j))
    hC = collect(dC)
    @test vec(transpose(A)) == hC
  end
  @testset "Test move CuDense on/off GPU" begin
    A  = [SType(1.0) for ii in 1:dim(i), jj in 1:dim(j)]
    dA = ITensorsGPU.CuDense{SType, CuVector{SType}}(NDTensors.Dense(vec(A)))
    dB = ITensorsGPU.Dense{SType, Vector{SType}}(dA)
    @test NDTensors.data(dB) == vec(A)
  end
  @testset "Test basic CuDense features" begin
    @test NDTensors.Dense{SType, CuVector{SType}}(10) isa ITensorsGPU.CuDense{SType}
    @test complex(NDTensors.Dense{SType, CuVector{SType}}) == NDTensors.Dense{complex(SType), CuVector{complex(SType), Nothing}}
  end
  if SType == Float64
      @testset "Test CuDense complex" begin
        A  = CuArrays.rand(SType, dim(i)*dim(j))
        dA = ITensorsGPU.CuDense{SType, CuVector{SType}}(A)
        dC = complex(dA)
        @test typeof(dC) !== typeof(dA)
        cdC = CuArray(dC)
        hC  = collect(cdC)
        ccA = complex.(A)
        @test hC == collect(ccA)
      end
  end
end # End Dense storage test


text:
using WAV
using SpeechAugment: farfieldWav

# read a wav file
wav,fs = wavread("xxx.wav")
# simulate far-field effects by auto-gain-ctrl method.
# recommended maxvalue is ∈ [0.5,1.0]
y = farfieldWav(wav, fs; maxvalue=0.6)

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
"""
_clean_path_dots(dir)

FTP directories use hidden dot folders frequently, function removes dots for local use
"""
function _clean_path_dots(dir)
    return abspath(replace(dir, "." => "")) # Remove . from folders to un-hide them
end

"""
    Mission

Each mission must have some defined functions:

`_mission_name(<:Mission)` - returns the name of the mission

`_mission_master_url(<:Mission)` - returns the url to the HEASARC master table (.tdat.gz) file

`_mission_paths(<:Mission)` - returns a `NamedTuple` with entries corresponding to folder and file paths:
* download - path the observation folder will be downloaded to from the FTP
* jaxtam   - path JAXTAM-generated files will be saved to
* web      - path the web report files will be saved to
* rmf      - path to the **local** copy of the mission RMF file from caldb

`_obs_path_server(<:Mission, obs_row::Union{DataFrame,DataFrameRow{DataFrame}})` - returns the **server-side** path to an observation

`_mission_good_e_range(<:Mission)` - returns a `Tuple{Float64,Float64}` of the **low** and **high** good energy ranges, respectively

`_mission_instruments(<:Mission)` - returns a `Array{Symbol1,}` with each element as a symbol of the instrument name, e.g. [:FPMA, :FPMB] for NuSTAR

Look at /src/missions/nicer.jl for an example. Put any custom mission functions in /src/missions/custom/ folder.
"""
abstract type Mission end

[include("base/$mission")   for mission in readdir(joinpath(JAXTAM.__sourcedir__,   "src/missions/base"))]
[include("custom/$mission") for mission in readdir(joinpath(JAXTAM.__sourcedir__, "src/missions/custom"))]

_mission_log(log_path::String=joinpath(JAXTAM.__sourcedir__, "mission_paths.json")) = JSON.parsefile(log_path)

mission_paths(; log_path::String=joinpath(JAXTAM.__sourcedir__, "mission_paths.json")) = _mission_log(log_path)

function mission_paths(mission::Mission; log_path::String=joinpath(JAXTAM.__sourcedir__, "mission_paths.json"), overwrite::Bool=false)
    mission_name = _mission_name(mission)

    if !isfile(log_path)
        log_dict = Dict()
    else
        log_dict = JSON.parsefile(log_path)
    end

    if !haskey(log_dict, mission_name) || overwrite
        @warn "Mission not found in $log_path, please enter paths:"
        print("Download path: ");                 path_download = readline()
        print("JAXTAM (processed data) path: ");  path_jaxtam   = readline()
        print("Web (html reports) path: ");       path_web      = readline()
        print("RMF (caldb mission file) path: "); path_rmf      = readline()
        log_dict[mission_name] = (download=path_download, jaxtam=path_jaxtam, web=path_web, rmf=path_rmf)
        write(log_path, JSON.json(log_dict, 4))
        @info "Wrote to $log_path" "Add custom keys in to JSON file manually if required"
    end

    # Return named tuple to prevent user from attempting to change the paths
    nt_keys = Tuple(Symbol.(collect(keys(log_dict[mission_name]))))
    return NamedTuple{nt_keys}((values(log_dict[mission_name])))
end

function _obs_path_local(mission::Mission, obs_row::Union{DataFrames.DataFrame,DataFrameRow{DataFrame}}; kind::Symbol)
    return joinpath(mission_paths(mission)[kind], _clean_path_dots(_obs_path_server(mission, obs_row))[2:end])
end

function mission_summary(mission::Mission)
    name = string("Name: ", _mission_name(mission))
    master_url = string("URL: ", _mission_master_url(mission))
    paths = string("Paths: \n\t\t", join(["$k: $v" for (k,v) in pairs(mission_paths(mission))], "\n\t\t"))
    e_range = string("Good energy range: ", _mission_good_e_range(mission), " keV")
    instruments = string("Instruments: ", join([i for i in _mission_instruments(mission)], ", "))

    println("$name\n\t$master_url\n\t$paths\n\t$e_range\n\t$instruments")
end

show(mission::Mission) = JAXTAM.mission_sumary(mission)

function _mission_symbol_to_type(mission_symbol::Symbol)
    if isdefined(JAXTAM, mission_symbol)
        return getproperty(JAXTAM, mission_symbol)
    end
end


text:
module TestGenieAWS

using Genie, Logging, LoggingExtras

function main()
  Base.eval(Main, :(const UserApp = TestGenieAWS))

  Genie.genie(; context = @__MODULE__)

  Base.eval(Main, :(const Genie = TestGenieAWS.Genie))
  Base.eval(Main, :(using Genie))
end

end


text:
using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer
using Test

SDL2_pkg_dir = joinpath(@__DIR__, "..","..")

@testset "Init+Quit" begin
# Test that you can init and quit SDL_ttf multiple times correctly.
@test 0 == SDL2.Init(UInt32(SDL2.INIT_VIDEO))
@test 0 == SDL2.TTF_Init()
SDL2.TTF_Quit()
SDL2.Quit()

@test 0 == SDL2.Init(UInt32(SDL2.INIT_VIDEO))
@test 0 == SDL2.TTF_Init()
SDL2.TTF_Quit()
SDL2.Quit()
end

@testset "Text" begin
@test 0 == SDL2.Init(UInt32(SDL2.INIT_VIDEO))
@test 0 == SDL2.TTF_Init()

@testset "Simple text" begin
font = TTF_OpenFont(joinpath(SDL2_pkg_dir,"assets/fonts/FiraCode/ttf/FiraCode-Regular.ttf"), 14)
@test font != C_NULL
txt = "Hello World!"
text = TTF_RenderText_Blended(font, txt, SDL2.Color(20,20,20,255))
@test text != C_NULL

fx,fy = Int[1], Int[1]
@test 0 == TTF_SizeText(font, txt, pointer(fx), pointer(fy))
fx,fy = fx[1],fy[1]
@test fx > 0
@test fy > 0
end

SDL2.TTF_Quit()
SDL2.Quit()
end

@testset "Rendering" begin
@test 0 == SDL2.Init(UInt32(SDL2.INIT_VIDEO))
@test 0 == SDL2.TTF_Init()
win = SDL2.CreateWindow("Hello World!", Int32(100), Int32(100), Int32(300), Int32(400),
         UInt32(SDL2.WINDOW_SHOWN))
renderer = SDL2.CreateRenderer(win, Int32(-1),
             UInt32(SDL2.RENDERER_ACCELERATED | SDL2.RENDERER_PRESENTVSYNC))

# try rendering with complicated text
@testset "Special Characters text" begin
font = TTF_OpenFont(joinpath(SDL2_pkg_dir,"assets/fonts/FiraCode/ttf/FiraCode-Regular.ttf"), 14)
txt = "@BinDeps.install Dict([ (:glib, :libglib) ])"
text = TTF_RenderText_Blended(font, txt, SDL2.Color(20,20,20,255))
@test text != C_NULL
tex = SDL2.CreateTextureFromSurface(renderer,text)
@test tex != C_NULL

fx,fy = Int[1], Int[1]
@test 0 == TTF_SizeText(font, txt, pointer(fx), pointer(fy))
fx,fy = fx[1],fy[1]
@test fx > 0
@test fy > 0

@test 0 == SDL2.RenderCopy(renderer, tex, C_NULL, pointer_from_objref(SDL2.Rect(0,0,fx,fy)))
SDL2.RenderPresent(renderer)

end

SDL2.TTF_Quit()
SDL2.Quit()
end


text:
const BSON_TYPE_DOUBLE = 0x01
const BSON_TYPE_STRING = 0x02
const BSON_TYPE_DOCUMENT = 0x03
const BSON_TYPE_ARRAY = 0x04
const BSON_TYPE_BINARY = 0x05
const BSON_TYPE_UNDEFINED = 0x06
const BSON_TYPE_OBJECTID = 0x07
const BSON_TYPE_BOOL = 0x08
const BSON_TYPE_DATETIME = 0x09
const BSON_TYPE_NULL = 0x0A
const BSON_TYPE_REGEX = 0x0B
const BSON_TYPE_DB_POINTER = 0x0C
const BSON_TYPE_CODE = 0x0D
const BSON_TYPE_SYMBOL = 0x0E
const BSON_TYPE_CODE_WITH_SCOPE = 0x0F
const BSON_TYPE_INT32 = 0x10
const BSON_TYPE_TIMESTAMP = 0x11
const BSON_TYPE_INT64 = 0x12
const BSON_TYPE_DECIMAL128 = 0x13
const BSON_TYPE_MIN_KEY = 0xFF
const BSON_TYPE_MAX_KEY = 0x7F

const BSON_SUBTYPE_GENERIC_BINARY = 0x00
const BSON_SUBTYPE_FUNCTION = 0x01
const BSON_SUBTYPE_BINARY_OLD = 0x02
const BSON_SUBTYPE_UUID_OLD = 0x03
const BSON_SUBTYPE_UUID = 0x04
const BSON_SUBTYPE_MD5 = 0x05
const BSON_SUBTYPE_ENCRYPTED = 0x06

const UnsafeBSONString = WeakRefString{UInt8}

struct BSONTimestamp
    counter::UInt32
    time::UInt32
end

BSONTimestamp(x::UInt64) = BSONTimestamp(x % UInt32, (x >> 32) % UInt32)

struct BSONCode
    code::String
end

struct BSONCodeWithScope
    code::String
    mappings::Dict{String, Any}
end

struct BSONBinary
    data::Vector{UInt8}
    subtype::UInt8
end

BSONBinary(data::Vector{UInt8}) = BSONBinary(data, BSON_SUBTYPE_GENERIC_BINARY)

struct UnsafeBSONBinary
    data::UnsafeArray{UInt8, 1}
    subtype::UInt8
end

UnsafeBSONBinary(data::UnsafeArray{UInt8, 1}) = UnsafeBSONBinary(data, BSON_SUBTYPE_GENERIC_BINARY)

struct BSONRegex
    pattern::String
    options::String
end

struct BSONSymbol
    value::String
end

struct BSONDBPointer
    collection::String
    ref::BSONObjectId
end

struct BSONUUIDOld
    value::UUID
end

struct BSONMinKey end
struct BSONMaxKey end
struct BSONUndefined end

const ValueField = Union{
    Float64,
    Int64,
    Int32,
    Bool,
    DateTime,
    Dec128,
    UUID,
    String,
    UnsafeBSONString,
    Nothing,
    BSONTimestamp,
    BSONObjectId,
    BSONBinary,
    UnsafeBSONBinary,
    BSONRegex,
    BSONCode,
    BSONSymbol,
    BSONMinKey,
    BSONMaxKey,
    BSONUndefined,
    BSONDBPointer,
    BSONUUIDOld
}

bson_type_(::Type{Float64}) = BSON_TYPE_DOUBLE
bson_type_(::Type{Int64}) = BSON_TYPE_INT64
bson_type_(::Type{Int32}) = BSON_TYPE_INT32
bson_type_(::Type{Bool}) = BSON_TYPE_BOOL
bson_type_(::Type{Dec128}) = BSON_TYPE_DECIMAL128
bson_type_(::Type{UUID}) = BSON_TYPE_BINARY
bson_type_(::Type{BSONUUIDOld}) = BSON_TYPE_BINARY
bson_type_(::Type{DateTime}) = BSON_TYPE_DATETIME
bson_type_(::Type{Nothing}) = BSON_TYPE_NULL
bson_type_(::Type{String}) = BSON_TYPE_STRING
bson_type_(::Type{UnsafeBSONString}) = BSON_TYPE_STRING
bson_type_(::Type{BSONTimestamp}) = BSON_TYPE_TIMESTAMP
bson_type_(::Type{BSONBinary}) = BSON_TYPE_BINARY
bson_type_(::Type{UnsafeBSONBinary}) = BSON_TYPE_BINARY
bson_type_(::Type{BSONRegex}) = BSON_TYPE_REGEX
bson_type_(::Type{BSONCode}) = BSON_TYPE_CODE
bson_type_(::Type{BSONSymbol}) = BSON_TYPE_SYMBOL
bson_type_(::Type{BSONObjectId}) = BSON_TYPE_OBJECTID
bson_type_(::Type{BSONMinKey}) = BSON_TYPE_MIN_KEY
bson_type_(::Type{BSONMaxKey}) = BSON_TYPE_MAX_KEY
bson_type_(::Type{BSONUndefined}) = BSON_TYPE_UNDEFINED
bson_type_(::Type{BSONDBPointer}) = BSON_TYPE_DB_POINTER


text:
struct BoxConstraints{V1, V2}
    lb::V1
    ub::V2
    al::V1
    au::V2
end
function BoxConstraints(lb::V1, ub::V2) where {V1, V2}
    BoxConstraints(lb, ub,
                   convert(V1, [!isfinite(lb[i]) ? 1 : min((ub[i] - lb[i])/2, (1 + abs(lb[i]))/20)
                                for i in eachindex(lb)]),
                   convert(V2, [!isfinite(ub[i]) ? 1 : min((ub[i] - lb[i])/2, (1 + abs(ub[i]))/20)
                                for i in eachindex(ub)])
                  )
end
_constraints(::Nothing, ::Nothing) = nothing
_constraints(l::AbstractVector, ::Any) = BoxConstraints(l, fill(Inf, length(l)))
_constraints(::Any, u::AbstractVector) = BoxConstraints(fill(-Inf, length(u)), u)
_constraints(l::AbstractVector, u::AbstractVector) = BoxConstraints(l, u)
transform(c, x) = transform!(c, copy(x))
backtransform(c, x) = backtransform!(c, copy(x))
transform!(::Nothing, x) = x
backtransform!(::Nothing, x) = x
function _shift_or_mirror_into_invertible(x, lb, al, ub, au)
    if x < lb - 2 * al - (ub - lb) / 2 || x > ub + 2 * au + (ub - lb) / 2
        r = 2 * (ub - lb + al + au)  # period
        s = lb - 2 * al - (ub - lb) / 2  # start
        x -= r * floor((x - s) / r)  # shift
    end
    if x > ub + au
        x -= 2 * (x - ub - au)
    end
    if x < lb - al
        x += 2 * (lb - al - x)
    end
    x
end
function _linquad_transform(x, lb, al, ub, au)
    x = _shift_or_mirror_into_invertible(x, lb, al, ub, au)
    x < lb + al && return lb + (x - (lb - al))^2 / 4 / al
    x < ub - au && return x
    x < ub + 3au && return ub - (x - (ub + au))^2 / 4 / au
    ub + au - (x - (ub + au))
end
function _linquad_transform_inverse(y, lb, al, ub, au)
    lb <= y <= ub || error("$y is not within [$lb, $ub].")
    y < lb + al && return (lb - al) + 2*√(al * (y - lb))
    y < ub - au && return y
    (ub + au) - 2*√(au * (ub - y))
end
function transform!(b::BoxConstraints, x::AbstractVector)
    for i in eachindex(x)
        x[i] = _linquad_transform(x[i], b.lb[i], b.al[i], b.ub[i], b.au[i])
    end
    x
end
function backtransform!(b::BoxConstraints, y::AbstractVector)
    for i in eachindex(y)
        y[i] = _linquad_transform_inverse(y[i], b.lb[i], b.al[i], b.ub[i], b.au[i])
    end
    y
end
function transform!(b::BoxConstraints, x::AbstractMatrix)
    n = size(x, 1)
    for i in eachindex(x)
        j = (i - 1) % n + 1
        x[i] = _linquad_transform(x[i], b.lb[j], b.al[j], b.ub[j], b.au[j])
    end
    x
end
function backtransform!(b::BoxConstraints, y::AbstractMatrix)
    n = size(x, 1)
    for i in eachindex(y)
        j = (i - 1) % n + 1
        y[i] = _linquad_transform_inverse(y[i], b.lb[j], b.al[j], b.ub[j], b.au[j])
    end
    y
end

max_stars_repo_name: jbrea/CMAEvolutionStrategy

max_issues_repo_name: jbrea/CMAEvolutionStrategy

max_forks_repo_name: jbrea/CMAEvolutionStrategy


text:
# SPDX-License-Identifier: MIT

"""
    elbo(model, X[, detailed = false, stats_scale = 1, params = params])

Compute the Evidence Lower-BOund (ELBO) of the model. If `detailed` is
set to `true` returns a tuple `elbo, loglikelihood, KL`. If
`params` is provided the function will return the natural gradient
for those parameters.
"""
elbo

function elbo(model, args...; cache = Dict(), detailed = false, stats_scale = 1)
    llh = loglikelihood(model, args..., cache)
    T = eltype(llh)
    sllh = sum(llh)*T(stats_scale)
    cache[:stats_scale] = stats_scale

    params = filter(isbayesianparam, getparams(model))
    KL = sum([EFD.kldiv(param.posterior, param.prior, μ = param.μ)
              for param in params])
    detailed ? (sllh - KL, sllh, KL) : sllh - KL
end

function _diagonal(param)
    d = similar(EFD.realform(param))
    fill!(d, 1)
    Diagonal(d)
end

function ∇elbo(model, cache, params)
    grads_Tμ = ∇sum_loglikelihood(model, cache)
    grads = Dict()
    for param in params
        ηq = EFD.naturalform(param.posterior.param)
        ηp = EFD.naturalform(param.prior.param)
        ∂KL_∂Tμ = (ηq - ηp)
        ∂𝓛_∂Tμ = grads_Tμ[param] * cache[:stats_scale] - ∂KL_∂Tμ
        #J = EFD.jacobian(param.posterior.param)
        J = _diagonal(param.posterior.param)
        grads[param] = J * ∂𝓛_∂Tμ
    end
    grads
end

"""
    gradstep(param_grad; lrate)

Update the parameters' posterior by doing one natural gradient steps.
"""
function gradstep(param_grad; lrate::Real)
    for (param, ∇𝓛) in param_grad
        ξ = param.posterior.param.ξ
        ξ[:] = ξ + lrate*∇𝓛
        param.μ[:] = EFD.gradlognorm(param.posterior)
    end
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
using LogitNormals
using Distributions, DistributionFits
using Test

@testset "LogitNormals.jl" begin
    include("logitnormal.jl")
end


text:
## Doggo leg model
using ForwardDiff, DiffResults
using Plots
using LinearAlgebra
using BenchmarkTools
using PartedArrays
using MatrixCalculus
using Test

include("doggo_model.jl")

# model
g = 9.81

L1 = 0.5
L2 = 0.5
L3 = 1.5
L4 = L3

l1 = L1/2
l2 = L2/2
l3 = L3/2
l4 = L4/2

m1 = 1.0
m2 = 1
m3 = 1
m4 = 1

m = [m1, m2, m3, m4]

J1 = 1
J2 = 1
J3 = 1
J4 = 1

M1 = Diagonal([m1*ones(2);J1])
M2 = Diagonal([m2*ones(2);J2])
M3 = Diagonal([m3*ones(2);J3])
M4 = Diagonal([m4*ones(2);J4])

M = [M1, M2, M3, M4]

# M = cat(M1,M2,M3,M4,dims=(1,2))

# Create Doggo Model
doggo = Doggo(L1,L3, [m1,m2,m3,m4], [J1,J2,J3,J4])

# inputs
a = pi/3 + rand()/5
b = -pi/4 + rand()/5
q = [a;b]

ȧ = 1.0
ḃ = 1.0
q̇ = [ȧ;ḃ]

# Forward Kinematics
function theta(a,b)
    d = L1*sin(0.5*(a-b))
    h1 = L1*cos(0.5*(a-b))
    h2 = sqrt(L3^2 -d^2)
    h = h1 + h2
    atan(2*h*d,L1^2 + L3^2 -h^2)
end

function trigtheta(a,b)
    d = L1*sin(0.5*(a-b))
    h1 = L1*cos(0.5*(a-b))
    h2 = sqrt(L3^2 -d^2)
    h = h1 + h2
    hyp = sqrt((2*h*d)^2 + (L1^2 + L3^2 - h^2)^2)
    st = 2*h*d/hyp
    ct = (L1^2 + L3^2 - h^2)/hyp
    return st, ct
end

function trigtheta2(a,b,L1,L3)
    d = L1*sin(0.5*(a-b))
    h1 = L1*cos(0.5*(a-b))
    h2 = sqrt(L3^2 -d^2)
    h = h1 + h2
    st = h*d/(L1*L3)
    ct = (L1^2 + L3^2 - h^2)/(2*L1*L3)
    return st, ct
end

theta(q) = theta(q...)

r1(a,b) = [l1*sin(a), -l1*cos(a), a]

r2(a,b) = [l2*sin(b), -l2*cos(b), b]

function r3_(a,b)
    a2 = theta(a,b) + a - pi
    t = theta(a,b)
    [L1*sin(a) + l3*sin(a2), -L1*cos(a) + -l3*cos(a2), a + t - pi]
end

function r3(a,b)
    st,ct = trigtheta(a,b)
    t = theta(a,b)
    [(L1 - l3*ct)*sin(a) - l3*st*cos(a),
     (-L1 + l3*ct)*cos(a) - l3*st*sin(a),
     a + t - pi]
 end

function r4(a,b)
    b2 = pi - theta(a,b) + b
    t = theta(a,b)
    [L2*sin(b) + l4*sin(b2), -L2*cos(b) + -l4*cos(b2), b - t + pi]
end

function rt(a,b)
    st,ct = trigtheta(a,b)
    a2 = theta(a,b) + a - pi
    b2 = pi - theta(a,b) + b
    [L1*sin(a) + L3*sin(a2), L1*cos(a) + L3*cos(a2), ]
end

r1(q) = r1(q...)
r2(q) = r2(q...)
r3(q) = r3(q...)
r4(q) = r4(q...)
rt(q) = rt(q...)

# Test the FK derivatives
fk_fun = [r1,r2,r3,r4]
val  = [r(q) for r in fk_fun]
grad = [ForwardDiff.jacobian(r,q) for r in fk_fun]
hess = [ForwardDiff.jacobian(q->vec(ForwardDiff.jacobian(r,q)),q) for r in fk_fun]
fk(doggo, q) ≈ val
fk_hess(doggo, q)[1] ≈ grad
fk_hess(doggo, q)[2] ≈ hess
jacobian(doggo, q) ≈ grad
jac = jacobian(doggo, q)

# @btime [ForwardDiff.jacobian(q->vec(ForwardDiff.jacobian(r,q)),$q) for r in $fk]
# @btime fk_hess($doggo, $q)

# Dynamics
x = [q; q̇]
part = create_partition2((2,2),(:x,:v))

d,h = leg_dims(doggo, q)
@inferred leg_dims(doggo, q)
@inferred trigtheta(doggo, d, h)
@inferred theta(doggo, d, h)
@inferred fk(doggo, q)
@inferred jacobian(doggo, q)
@inferred fk_hess(doggo, q)
@inferred mass_matrix(doggo, jac)
@inferred mass_matrix(doggo, q)
@inferred lagrangian(doggo, q, q̇)
@inferred get_V(doggo, q)

# Check gradient of lagrangian
lagrangian(x) = lagrangian(doggo, x[1:2], x[3:4])
L_grad = ForwardDiff.gradient(lagrangian, x)
grad_L(doggo, q, q̇) ≈ L_grad

# Check potential energy gradient
get_V(q) = get_V(doggo, q)
ForwardDiff.gradient(get_V, q) ≈ get_∇V(doggo, jac)

# Test mass matrix
part = NamedTuple{(:xx,:xv,:vx,:vv)}(create_partition2((2,2)))
L_hess = BlockArray(ForwardDiff.hessian(lagrangian, x), part)
M_ = mass_matrix(doggo, q)
M_ ≈ L_hess.vv

# Test full euler lagrange equation (assuming some qdd)
qdd = rand(2)
qd = q̇
el = L_hess.vx*q̇ + L_hess.vv*qdd - L_grad[1:2]
euler_lagrange_auto(doggo, q, qd, qdd) == el
euler_lagrange(doggo, q, qd, qdd) ≈ el
@inferred euler_lagrange(doggo, q, qd, qdd)
@inferred euler_lagrange_auto(doggo, q, qd, qdd)
@btime euler_lagrange_auto($doggo, $q, $qd, $qdd)
@btime euler_lagrange($doggo, $q, $qd, $qdd)

@code_warntype create_partition2((2,2),(2,2),(:x,:v),(:x,:v))
names_all = create_partition2((2,2),(2,2),(:x,:v),(:x,:v))
@code_warntype PartedArrays.create_nt(names_all,part)

res = DiffResults.HessianResult(x)
ForwardDiff.hessian!(res, lagrangian, x)
DiffResults.gradient(res) ≈ L_grad

dynamics(doggo, q, qd, [0,0])
n,m = 4,2
f(x,u) = [x[3:4]; dynamics(doggo, x[1:2], x[3:4], u)]
fd = rk3(f, 0.1)
x = [deg2rad(45), deg2rad(-45), 0, 0]
u = zeros(2)
fd(x,u)

function simulate(x0,u,N, doplot=false)
    X = zeros(n,N)
    X[:,1] = x0
    for k = 1:N-1
        X[:,k+1] = fd(X[:,k],u)
        if doplot
            p = plot(doggo, X[:,k])
            display(p)
        end
    end
    return X
end
u = [1,2]*10
X = simulate(x,u,51, true)
plot(doggo, X[:,51])

plot(doggo, [deg2rad(50), deg2rad(-90)], 2)

# Second Derivatives
ForwardDiff.hessian(theta,q)
hess = ForwardDiff.hessian(sintheta,q)
hess = ForwardDiff.hessian(costheta,q)[1,1]
ForwardDiff.gradient(sintheta,q)
phi = (a-b)/2

dsdaa = 1/(4L3)*((-(d^2/h2^2 + 3)*h1*d/h2)*cos(phi) - (h-d^2/h2)*sin(phi) - 2h1*sin(phi) - d*cos(phi))
dcdaa = 1/4L3*(1/h2*(-d^2 + h*h1 + (h/h2 - 1)*d^2*h1/h2 + h*h2)*cos(phi) - d*(h/h2 + (1+h1/h2))*sin(phi))
dsdbb = dsdaa
dcdbb = dcdaa
dsdab = -dsdaa
dcdab = -dcdaa

dx3aa = -(L1-l3*ct)* sin(a) - l3*cos(a)*dcda +
             l3*st * cos(a) + l3*sin(a)*dsda -
             l3*dcda*cos(a) - l3*sin(a)*dcdaa +
             l3*dsda*sin(a) - l3*cos(a)*dsdaa

dx3ab = -l3*cos(a)*dcdb  + l3*sin(a)*dsdb -
         l3*sin(a)*dcdab - l3*cos(a)*dsdab

dx3bb = -l3*sin(a)*dcdbb - l3*cos(a)*dsdbb

rx3(q) = r3(q)[1]
ForwardDiff.hessian(rx3,q) ≈ [dx3aa dx3ab; dx3ab dx3bb]

dy3a = (L1-l3*ct)*sin(a) - l3*st*cos(a) + l3*cos(a)*dcda - l3*sin(a)*dsda
dy3b = l3*cos(a)*dcdb - l3*sin(a)*dsdb

dy3aa = (L1-l3*ct)*cos(a)  - l3*sin(a)*dcda +
            l3*st*sin(a)   - l3*cos(a)*dsda -
            l3*sin(a)*dcda + l3*cos(a)*dcdaa -
            l3*cos(a)*dsda - l3*sin(a)*dsdaa
dy3ab = -l3*sin(a)*dcdb - l3*cos(a)*dsdb +
         l3*cos(a)*dcdab - l3*sin(a)*dsdab
dy3bb = l3*cos(a)*dcdbb -l3*sin(a)*dsdbb

ry3(q) = r3(q)[2]
ForwardDiff.gradient(ry3,q) ≈ [dy3a,dy3b]
ForwardDiff.hessian(ry3,q) ≈ [dy3aa dy3ab; dy3ab dy3bb]

dx4a = -l4*sin(b)*dcda + l4*cos(b)*dsda
dx4b = (L2 - l4*ct)*cos(b) - l4*st*sin(b) - l4*sin(b)*dcdb + l4*cos(b)*dsdb

dx4aa = -l4*sin(b)*dcdaa + l4*cos(b)*dsdaa
dx4ab = -l4*dcda*cos(b) - l4*sin(b)*dcdab +
        -l4*dsda*sin(b) + l4*cos(b)*dsdab
dx4bb = -(L2-l3*ct)*sin(b) - l4*cos(b)*dcdb +
        -l4*st*cos(b)      - l4*sin(b)*dsdb +
        -l3*dcdb*cos(b)    - l3*sin(b)*dcdbb +
        -l4*dsdb*sin(b)    + l3*cos(b)*dsdbb

rx4(q) = r4(q)[1]
ForwardDiff.gradient(rx4,q) ≈ [dx4a,dx4b]
ForwardDiff.hessian(rx4,q) ≈ [dx4aa dx4ab; dx4ab dx4bb]

dy4a = l4*cos(b)*dcda + l4*sin(b)*dsda
dy4b = (L2 - l4*ct)*sin(b) + l4*st*cos(b) + l4*cos(b)*dcdb + l4*sin(b)*dsdb

dy4aa = l4*cos(b)*dcdaa + l4*sin(b)*dsdaa
dy4ab = -l4*dcda*sin(b) + l4*cos(b)*dcdab +
         l4*dsda*cos(b) + l4*sin(b)*dsdab
dy4bb = (L2-l4*ct)*cos(b) - l4*sin(b)*dcdb +
        -l4*st*sin(b)     + l4*cos(b)*dsdb +
        -l4*dcdb*sin(b)   + l4*cos(b)*dcdbb +
         l4*dsdb*cos(b)   + l4*sin(b)*dsdbb

ry4(q) = r4(q)[2]

println("passed")

## generate function for four bar mechanism joint positions as function of α, β
function gen_four_bar(L1::T,L3::T) where T
    L2 = L1
    L4 = L3
    l3 = L3/2
    l4 = L4/2

    function theta(a::T,b::T) where t
        d = L1*sin(0.5*(a-b))
        h1 = L1*cos(0.5*(a-b))
        h2 = sqrt(L3^2 -d^2)
        h = h1 + h2
        atan(2*h*d,L1^2 + L3^2 -h^2)
    end

    function pos(a::
        T,b::T) where T
        t = theta(a,b)

        a2 = t + a - pi
        b2 = pi - t + b

        # right elbow
        x1 = L1*sin(a)
        y1 = -L1*cos(a)

        #left elbow
        x2 = L2*sin(b)
        y2 = -L2*cos(b)

        # bottom joint
        x3 = L1*sin(a) + l3*sin(a2)
        y3 = -L1*cos(a) -l3*cos(a2)

        x4 = L2*sin(b) + l4*sin(b2)
        y4 = -L2*cos(b) - l4*cos(b2)

        # tip
        xt = L1*sin(a) + L3*sin(a2)
        yt = -L1*cos(a) -L3*cos(a2)

        [(0.0,0.0), (x1,y1), (x2,y2), (xt,yt)]
    end

    function pos(a::T,b::T,c::Bool) where T

        # right elbow
        x1 = L1*sin(a)
        y1 = -L1*cos(a)

        #left elbow
        x2 = L2*sin(b)
        y2 = -L2*cos(b)

        st, ct = trigtheta2(a,b,L1,L3)

        x3 = (L1 - l3*ct)*sin(a) - l3*st*cos(a)
        y3 = (-L1 + l3*ct)*cos(a) - l3*st*sin(a)

        x4 = (L2 - l4*ct)*sin(b) + l4*st*cos(b)
        y4 = (-L2 + l4*ct)*cos(b) + l4*st*sin(b)

        xt = (L2 - L4*ct)*sin(b) + L4*st*cos(b)
        yt = (-L2 + L4*ct)*cos(b) + L4*st*sin(b)

        [(0.0,0.0), (x1,y1), (x2,y2), (xt,yt)]
    end
    return pos
end

function plot_four_bar(f::Function,α::T,β::T) where T
    pos = f(α,β)
    p = plot(title="Doggo Leg alpha=$(round(α,digits=2)),beta=$(round(β,digits=2))",xlim=(-2,2),ylim=(-2,2),aspectratio=:equal,label="")
    plot!([pos[1][1],pos[2][1]],[pos[1][2],pos[2][2]],color=:black,linewidth=2,label="")
    plot!([pos[1][1],pos[3][1]],[pos[1][2],pos[3][2]],color=:black,linewidth=2,label="")
    plot!([pos[2][1],pos[4][1]],[pos[2][2],pos[4][2]],color=:black,linewidth=2,label="")
    plot!([pos[3][1],pos[4][1]],[pos[3][2],pos[4][2]],color=:black,linewidth=2,label="")
end

four_bar = gen_four_bar(rand()+1,rand()+2)
four_bar(q[1],q[2])[4] .≈ four_bar(q[1],q[2],true)[4]
@btime four_bar($q[1],$q[2])
@btime four_bar($q[1],$q[2],true)
plot_four_bar(four_bar,pi/4,-pi/3)


text:
module dSFMT

export DSFMT_state, dsfmt_get_min_array_size, dsfmt_get_idstring,
       dsfmt_init_gen_rand, dsfmt_gv_init_gen_rand,
       dsfmt_init_by_array, dsfmt_gv_init_by_array,
       dsfmt_genrand_close1_open2, dsfmt_gv_genrand_close1_open2,
       dsfmt_genrand_close_open, dsfmt_gv_genrand_close_open,
       dsfmt_genrand_uint32, dsfmt_gv_genrand_uint32,
       randmtzig_randn, randmtzig_exprnd,
       win32_SystemFunction036!

type DSFMT_state
    val::Vector{Int32}
    DSFMT_state() = new(Array(Int32, 770))
end

function dsfmt_get_idstring()
    idstring = ccall((:dsfmt_get_idstring,:libdSFMT),
                     Ptr{Uint8},
                     ())
    return bytestring(idstring)
end

function dsfmt_get_min_array_size()
    min_array_size = ccall((:dsfmt_get_min_array_size,:libdSFMT),
                           Int32,
                           ())
end

const dsfmt_min_array_size = dsfmt_get_min_array_size()

function dsfmt_init_gen_rand(s::DSFMT_state, seed::Uint32)
    ccall((:dsfmt_init_gen_rand,:libdSFMT),
          Void,
          (Ptr{Void}, Uint32,),
          s.val, seed)
end

function dsfmt_gv_init_gen_rand(seed::Uint32)
    ccall((:dsfmt_gv_init_gen_rand,:libdSFMT),
          Void,
          (Uint32,),
          seed)
end

function dsfmt_init_by_array(s::DSFMT_state, seed::Vector{Uint32})
    ccall((:dsfmt_init_by_array,:libdSFMT),
          Void,
          (Ptr{Void}, Ptr{Uint32}, Int32),
          s.val, seed, length(seed))
end

function dsfmt_gv_init_by_array(seed::Vector{Uint32})
    ccall((:dsfmt_gv_init_by_array,:libdSFMT),
        Void,
        (Ptr{Uint32}, Int32),
        seed, length(seed))
end

function dsfmt_genrand_close_open(s::DSFMT_state)
    ccall((:dsfmt_genrand_close_open, :libdSFMT),
    Float64,
    (Ptr{Void},),
    s.val)
end

function dsfmt_gv_genrand_close_open()
    ccall((:dsfmt_gv_genrand_close_open, :libdSFMT),
    Float64,
    ())
end

function dsfmt_genrand_close1_open2(s::DSFMT_state)
    ccall((:dsfmt_genrand_close1_open2, :libdSFMT),
    Float64,
    (Ptr{Void},),
    s.val)
end

function dsfmt_gv_genrand_close1_open2()
    ccall((:dsfmt_gv_genrand_close1_open2, :libdSFMT),
    Float64,
    ())
end

function dsfmt_genrand_uint32(s::DSFMT_state)
    ccall((:dsfmt_genrand_uint32,:libdSFMT),
          Uint32,
          (Ptr{Void},),
          s.val)
end

function dsfmt_gv_genrand_uint32()
    ccall((:dsfmt_gv_genrand_uint32,:libdSFMT),
          Uint32,
          ())
end

# randmtzig
## Tables for normal variates
const ki =
    Uint64[0x0007799ec012f7b3,                 0,0x0006045f4c7de363,0x0006d1aa7d5ec0a6,
           0x000728fb3f60f778,0x0007592af4e9fbc0,0x000777a5c0bf655d,0x00078ca3857d2256,
           0x00079bf6b0ffe58c,0x0007a7a34ab092ae,0x0007b0d2f20dd1cb,0x0007b83d3aa9cb52,
           0x0007be597614224e,0x0007c3788631abea,0x0007c7d32bc192ef,0x0007cb9263a6e86d,
           0x0007ced483edfa84,0x0007d1b07ac0fd39,0x0007d437ef2da5fc,0x0007d678b069aa6e,
           0x0007d87db38c5c87,0x0007da4fc6a9ba63,0x0007dbf611b37f3c,0x0007dd7674d0f286,
           0x0007ded5ce8205f7,0x0007e018307fb62c,0x0007e141081bd124,0x0007e2533d712de9,
           0x0007e3514bbd7718,0x0007e43d54944b52,0x0007e5192f25ef43,0x0007e5e67481118d,
           0x0007e6a6897c1ce2,0x0007e75aa6c7f64d,0x0007e803df8ee498,0x0007e8a326eb6273,
           0x0007e93954717a28,0x0007e9c727f8648f,0x0007ea4d4cc85a3d,0x0007eacc5c4907a9,
           0x0007eb44e0474cf6,0x0007ebb754e47419,0x0007ec242a3d8474,0x0007ec8bc5d69646,
           0x0007ecee83d3d6e9,0x0007ed4cb8082f46,0x0007eda6aee01710,0x0007edfcae2dfe68,
           0x0007ee4ef5dccd3f,0x0007ee9dc08c394f,0x0007eee9441a17c8,0x0007ef31b21b4fb2,
           0x0007ef773846a8a7,0x0007efba00d35a17,0x0007effa32ccf69f,0x0007f037f25e1279,
           0x0007f0736112d12c,0x0007f0ac9e145c26,0x0007f0e3c65e1fcc,0x0007f118f4ed8e55,
           0x0007f14c42ed0dc9,0x0007f17dc7daa0c3,0x0007f1ad99aac6a5,0x0007f1dbcce80015,
           0x0007f20874cf56bf,0x0007f233a36a3b9a,0x0007f25d69a604ae,0x0007f285d7694a92,
           0x0007f2acfba75e3c,0x0007f2d2e4720909,0x0007f2f79f09c345,0x0007f31b37ec883c,
           0x0007f33dbae36abc,0x0007f35f330f08d5,0x0007f37faaf2fa79,0x0007f39f2c805381,
           0x0007f3bdc11f4f1d,0x0007f3db71b83851,0x0007f3f846bba121,0x0007f4144829f846,
           0x0007f42f7d9a8b9e,0x0007f449ee420432,0x0007f463a0f8675e,0x0007f47c9c3ea77b,
           0x0007f494e643cd8e,0x0007f4ac84e9c475,0x0007f4c37dc9cd51,0x0007f4d9d638a433,
           0x0007f4ef934a5b6b,0x0007f504b9d5f33e,0x0007f5194e78b352,0x0007f52d55994a96,
           0x0007f540d36aba0d,0x0007f553cbef0e78,0x0007f56642f9ec90,0x0007f5783c32f31e,
           0x0007f589bb17f609,0x0007f59ac2ff1525,0x0007f5ab5718b15a,0x0007f5bb7a71427d,
           0x0007f5cb2ff3100a,0x0007f5da7a67cebe,0x0007f5e95c7a24e8,0x0007f5f7d8b7171f,
           0x0007f605f18f5ef4,0x0007f613a958ad0a,0x0007f621024ed7ea,0x0007f62dfe94f8cb,
           0x0007f63aa036777a,0x0007f646e928065a,0x0007f652db488f88,0x0007f65e786213ff,
           0x0007f669c22a7d8b,0x0007f674ba44645a,0x0007f67f623fc8dc,0x0007f689bb9ac294,
           0x0007f693c7c22482,0x0007f69d881217a6,0x0007f6a6fdd6ac37,0x0007f6b02a4c61ee,
           0x0007f6b90ea0a7f4,0x0007f6c1abf254c1,0x0007f6ca03521664,0x0007f6d215c2db82,
           0x0007f6d9e43a355a,0x0007f6e16fa0b329,0x0007f6e8b8d23729,0x0007f6efc09e4569,
           0x0007f6f687c84cc0,0x0007f6fd0f07ea0a,0x0007f703570925e2,0x0007f709606cad03,
           0x0007f70f2bc80370,0x0007f714b9a5b292,0x0007f71a0a85725e,0x0007f71f1edc4d9e,
           0x0007f723f714c17a,0x0007f728938ed843,0x0007f72cf4a03fa1,0x0007f7311a945a17,
           0x0007f73505ac4bf8,0x0007f738b61f03bd,0x0007f73c2c193dc1,0x0007f73f67bd835d,
           0x0007f74269242559,0x0007f745305b31a1,0x0007f747bd666429,0x0007f74a103f12ed,
           0x0007f74c28d414f6,0x0007f74e0709a42e,0x0007f74faab939fa,0x0007f75113b16657,
           0x0007f75241b5a156,0x0007f753347e16b9,0x0007f753ebb76b7d,0x0007f75467027d05,
           0x0007f754a5f4199e,0x0007f754a814b208,0x0007f7546ce003af,0x0007f753f3c4bb2a,
           0x0007f7533c240e92,0x0007f75245514f41,0x0007f7510e91726d,0x0007f74f971a9012,
           0x0007f74dde135797,0x0007f74be2927972,0x0007f749a39e051d,0x0007f747202aba8b,
           0x0007f744571b4e3d,0x0007f741473f9efe,0x0007f73def53dc44,0x0007f73a4dff9c00,
           0x0007f73661d4deaf,0x0007f732294f0040,0x0007f72da2d19444,0x0007f728cca72bdb,
           0x0007f723a5000367,0x0007f71e29f09628,0x0007f7185970156c,0x0007f7123156c102,
           0x0007f70baf5c1e2d,0x0007f704d1150a24,0x0007f6fd93f1a4e6,0x0007f6f5f53b10b6,
           0x0007f6edf211023f,0x0007f6e587671cea,0x0007f6dcb202167a,0x0007f6d36e749c65,
           0x0007f6c9b91bf4c6,0x0007f6bf8e1c541b,0x0007f6b4e95ce015,0x0007f6a9c6835700,
           0x0007f69e20ef5211,0x0007f691f3b517ec,0x0007f6853997f322,0x0007f677ed03ff1a,
           0x0007f66a08075bdc,0x0007f65b844ab75a,0x0007f64c5b091860,0x0007f63c8506d4bc,
           0x0007f62bfa8798ff,0x0007f61ab34364b1,0x0007f608a65a599a,0x0007f5f5ca4737e8,
           0x0007f5e214d05b49,0x0007f5cd7af7066e,0x0007f5b7f0e4c2a1,0x0007f5a169d68fcf,
           0x0007f589d80596a6,0x0007f5712c8d0174,0x0007f557574c912b,0x0007f53c46c77193,
           0x0007f51fe7feb9f3,0x0007f5022646ecfb,0x0007f4e2eb17ab1d,0x0007f4c21dd4a3d2,
           0x0007f49fa38ea394,0x0007f47b5ebb62eb,0x0007f4552ee27473,0x0007f42cf03d58f5,
           0x0007f4027b4854a0,0x0007f3d5a44119e0,0x0007f3a63a8fb553,0x0007f37408155101,
           0x0007f33ed05b55ec,0x0007f3064f9c183f,0x0007f2ca399c7ba1,0x0007f28a384bb940,
           0x0007f245ea1b7a2c,0x0007f1fcdffe8f1c,0x0007f1ae9af758ce,0x0007f15a8917f27f,
           0x0007f10001ccaaac,0x0007f09e413c418a,0x0007f034627733d8,0x0007efc15815b8d5,
           0x0007ef43e2bf7f55,0x0007eeba84e31dff,0x0007ee237294df89,0x0007ed7c7c170141,
           0x0007ecc2f0d95d3b,0x0007ebf377a46782,0x0007eb09d6deb286,0x0007ea00a4f17809,
           0x0007e8d0d3da63d6,0x0007e771023b0fcf,0x0007e5d46c2f08d9,0x0007e3e937669691,
           0x0007e195978f1176,0x0007deb2c0e05c1d,0x0007db0362002a1a,0x0007d6202c15143a,
           0x0007cf4b8f00a2cc,0x0007c4fd24520efe,0x0007b362fbf81816,0x00078d2d25998e25]
const wi =
    [1.7367254121602630e-15,9.5586603514556339e-17,1.2708704834810623e-16,
     1.4909740962495474e-16,1.6658733631586268e-16,1.8136120810119029e-16,
     1.9429720153135588e-16,2.0589500628482093e-16,2.1646860576895422e-16,
     2.2622940392218116e-16,2.3532718914045892e-16,2.4387234557428771e-16,
     2.5194879829274225e-16,2.5962199772528103e-16,2.6694407473648285e-16,
     2.7395729685142446e-16,2.8069646002484804e-16,2.8719058904113930e-16,
     2.9346417484728883e-16,2.9953809336782113e-16,3.0543030007192440e-16,
     3.1115636338921572e-16,3.1672988018581815e-16,3.2216280350549905e-16,
     3.2746570407939751e-16,3.3264798116841710e-16,3.3771803417353232e-16,
     3.4268340353119356e-16,3.4755088731729758e-16,3.5232663846002031e-16,
     3.5701624633953494e-16,3.6162480571598339e-16,3.6615697529653540e-16,
     3.7061702777236077e-16,3.7500889278747798e-16,3.7933619401549554e-16,
     3.8360228129677279e-16,3.8781025861250247e-16,3.9196300853257678e-16,
     3.9606321366256378e-16,4.0011337552546690e-16,4.0411583124143332e-16,
     4.0807276830960448e-16,4.1198623774807442e-16,4.1585816580828064e-16,
     4.1969036444740733e-16,4.2348454071520708e-16,4.2724230518899761e-16,
     4.3096517957162941e-16,4.3465460355128760e-16,4.3831194100854571e-16,
     4.4193848564470665e-16,4.4553546609579137e-16,4.4910405058828750e-16,
     4.5264535118571397e-16,4.5616042766900381e-16,4.5965029108849407e-16,
     4.6311590702081647e-16,4.6655819856008752e-16,4.6997804906941950e-16,
     4.7337630471583237e-16,4.7675377680908526e-16,4.8011124396270155e-16,
     4.8344945409350080e-16,4.8676912627422087e-16,4.9007095245229938e-16,
     4.9335559904654139e-16,4.9662370843221783e-16,4.9987590032409088e-16,
     5.0311277306593187e-16,5.0633490483427195e-16,5.0954285476338923e-16,
     5.1273716399787966e-16,5.1591835667857364e-16,5.1908694086703434e-16,
     5.2224340941340417e-16,5.2538824077194543e-16,5.2852189976823820e-16,
     5.3164483832166176e-16,5.3475749612647295e-16,5.3786030129452348e-16,
     5.4095367096239933e-16,5.4403801186554671e-16,5.4711372088173611e-16,
     5.5018118554603362e-16,5.5324078453927836e-16,5.5629288815190902e-16,
     5.5933785872484621e-16,5.6237605106900435e-16,5.6540781286489604e-16,
     5.6843348504368141e-16,5.7145340215092040e-16,5.7446789269419609e-16,
     5.7747727947569648e-16,5.8048187991076857e-16,5.8348200633338921e-16,
     5.8647796628943653e-16,5.8947006281858718e-16,5.9245859472561339e-16,
     5.9544385684180598e-16,5.9842614027720281e-16,6.0140573266426640e-16,
     6.0438291839361250e-16,6.0735797884236057e-16,6.1033119259564394e-16,
     6.1330283566179110e-16,6.1627318168165963e-16,6.1924250213258470e-16,
     6.2221106652737879e-16,6.2517914260879998e-16,6.2814699653988953e-16,
     6.3111489309056042e-16,6.3408309582080600e-16,6.3705186726088149e-16,
     6.4002146908880247e-16,6.4299216230548961e-16,6.4596420740788321e-16,
     6.4893786456033965e-16,6.5191339376461587e-16,6.5489105502874154e-16,
     6.5787110853507413e-16,6.6085381480782587e-16,6.6383943488035057e-16,
     6.6682823046247459e-16,6.6982046410815579e-16,6.7281639938375311e-16,
     6.7581630103719006e-16,6.7882043516829803e-16,6.8182906940062540e-16,
     6.8484247305500383e-16,6.8786091732516637e-16,6.9088467545571690e-16,
     6.9391402292275690e-16,6.9694923761748294e-16,6.9999060003307640e-16,
     7.0303839345521508e-16,7.0609290415654822e-16,7.0915442159548734e-16,
     7.1222323861967788e-16,7.1529965167453030e-16,7.1838396101720629e-16,
     7.2147647093647067e-16,7.2457748997883870e-16,7.2768733118146927e-16,
     7.3080631231227429e-16,7.3393475611774048e-16,7.3707299057898310e-16,
     7.4022134917657997e-16,7.4338017116476479e-16,7.4654980185558890e-16,
     7.4973059291369793e-16,7.5292290266240584e-16,7.5612709640179217e-16,
     7.5934354673958895e-16,7.6257263393567558e-16,7.6581474626104873e-16,
     7.6907028037219191e-16,7.7233964170182985e-16,7.7562324486711744e-16,
     7.7892151409638524e-16,7.8223488367564108e-16,7.8556379841610841e-16,
     7.8890871414417552e-16,7.9227009821522709e-16,7.9564843005293662e-16,
     7.9904420171571300e-16,8.0245791849212591e-16,8.0589009952726568e-16,
     8.0934127848215009e-16,8.1281200422845008e-16,8.1630284158098775e-16,
     8.1981437207065329e-16,8.2334719476060504e-16,8.2690192710884700e-16,
     8.3047920588053737e-16,8.3407968811366288e-16,8.3770405214202216e-16,
     8.4135299867980282e-16,8.4502725197240968e-16,8.4872756101861549e-16,
     8.5245470086955962e-16,8.5620947401062333e-16,8.5999271183276646e-16,
     8.6380527620052589e-16,8.6764806112455816e-16,8.7152199454736980e-16,
     8.7542804025171749e-16,8.7936719990210427e-16,8.8334051523084080e-16,
     8.8734907038131345e-16,8.9139399442240861e-16,8.9547646404950677e-16,
     8.9959770648910994e-16,9.0375900262601175e-16,9.0796169037400680e-16,
     9.1220716831348461e-16,9.1649689962191353e-16,9.2083241632623076e-16,
     9.2521532390956933e-16,9.2964730630864167e-16,9.3413013134252651e-16,
     9.3866565661866598e-16,9.4325583596767065e-16,9.4790272646517382e-16,
     9.5260849610662787e-16,9.5737543220974496e-16,9.6220595062948384e-16,
     9.6710260588230542e-16,9.7206810229016259e-16,9.7710530627072088e-16,
     9.8221725991905411e-16,9.8740719604806711e-16,9.9267855488079765e-16,
     9.9803500261836449e-16,1.0034804521436181e-15,1.0090190861637457e-15,
     1.0146553831467086e-15,1.0203941464683124e-15,1.0262405372613567e-15,
     1.0322001115486456e-15,1.0382788623515399e-15,1.0444832676000471e-15,
     1.0508203448355195e-15,1.0572977139009890e-15,1.0639236690676801e-15,
     1.0707072623632994e-15,1.0776584002668106e-15,1.0847879564403425e-15,
     1.0921079038149563e-15,1.0996314701785628e-15,1.1073733224935752e-15,
     1.1153497865853155e-15,1.1235791107110833e-15,1.1320817840164846e-15,
     1.1408809242582780e-15,1.1500027537839792e-15,1.1594771891449189e-15,
     1.1693385786910960e-15,1.1796266352955801e-15,1.1903876299282890e-15,
     1.2016759392543819e-15,1.2135560818666897e-15,1.2261054417450561e-15,
     1.2394179789163251e-15,1.2536093926602567e-15,1.2688244814255010e-15,
     1.2852479319096109e-15,1.3031206634689985e-15,1.3227655770195326e-15,
     1.3446300925011171e-15,1.3693606835128518e-15,1.3979436672775240e-15,
     1.4319989869661328e-15,1.4744848603597596e-15,1.5317872741611144e-15,
     1.6227698675312968e-15]
const fi =
    [1.0000000000000000e+00,9.7710170126767082e-01,9.5987909180010600e-01,
     9.4519895344229909e-01,9.3206007595922991e-01,9.1999150503934646e-01,
     9.0872644005213032e-01,8.9809592189834297e-01,8.8798466075583282e-01,
     8.7830965580891684e-01,8.6900868803685649e-01,8.6003362119633109e-01,
     8.5134625845867751e-01,8.4291565311220373e-01,8.3471629298688299e-01,
     8.2672683394622093e-01,8.1892919160370192e-01,8.1130787431265572e-01,
     8.0384948317096383e-01,7.9654233042295841e-01,7.8937614356602404e-01,
     7.8234183265480195e-01,7.7543130498118662e-01,7.6863731579848571e-01,
     7.6195334683679483e-01,7.5537350650709567e-01,7.4889244721915638e-01,
     7.4250529634015061e-01,7.3620759812686210e-01,7.2999526456147568e-01,
     7.2386453346862967e-01,7.1781193263072152e-01,7.1183424887824798e-01,
     7.0592850133275376e-01,7.0009191813651117e-01,6.9432191612611627e-01,
     6.8861608300467136e-01,6.8297216164499430e-01,6.7738803621877308e-01,
     6.7186171989708166e-01,6.6639134390874977e-01,6.6097514777666277e-01,
     6.5561147057969693e-01,6.5029874311081637e-01,6.4503548082082196e-01,
     6.3982027745305614e-01,6.3465179928762327e-01,6.2952877992483625e-01,
     6.2445001554702606e-01,6.1941436060583399e-01,6.1442072388891344e-01,
     6.0946806492577310e-01,6.0455539069746733e-01,5.9968175261912482e-01,
     5.9484624376798689e-01,5.9004799633282545e-01,5.8528617926337090e-01,
     5.8055999610079034e-01,5.7586868297235316e-01,5.7121150673525267e-01,
     5.6658776325616389e-01,5.6199677581452390e-01,5.5743789361876550e-01,
     5.5291049042583185e-01,5.4841396325526537e-01,5.4394773119002582e-01,
     5.3951123425695158e-01,5.3510393238045717e-01,5.3072530440366150e-01,
     5.2637484717168403e-01,5.2205207467232140e-01,5.1775651722975591e-01,
     5.1348772074732651e-01,5.0924524599574761e-01,5.0502866794346790e-01,
     5.0083757512614835e-01,4.9667156905248933e-01,4.9253026364386815e-01,
     4.8841328470545758e-01,4.8432026942668288e-01,4.8025086590904642e-01,
     4.7620473271950547e-01,4.7218153846772976e-01,4.6818096140569321e-01,
     4.6420268904817391e-01,4.6024641781284248e-01,4.5631185267871610e-01,
     4.5239870686184824e-01,4.4850670150720273e-01,4.4463556539573912e-01,
     4.4078503466580377e-01,4.3695485254798533e-01,4.3314476911265209e-01,
     4.2935454102944126e-01,4.2558393133802180e-01,4.2183270922949573e-01,
     4.1810064983784795e-01,4.1438753404089090e-01,4.1069314827018799e-01,
     4.0701728432947315e-01,4.0335973922111429e-01,3.9972031498019700e-01,
     3.9609881851583223e-01,3.9249506145931540e-01,3.8890886001878855e-01,
     3.8534003484007706e-01,3.8178841087339344e-01,3.7825381724561896e-01,
     3.7473608713789086e-01,3.7123505766823922e-01,3.6775056977903225e-01,
     3.6428246812900372e-01,3.6083060098964775e-01,3.5739482014578022e-01,
     3.5397498080007656e-01,3.5057094148140588e-01,3.4718256395679348e-01,
     3.4380971314685055e-01,3.4045225704452164e-01,3.3711006663700588e-01,
     3.3378301583071823e-01,3.3047098137916342e-01,3.2717384281360129e-01,
     3.2389148237639104e-01,3.2062378495690530e-01,3.1737063802991350e-01,
     3.1413193159633707e-01,3.1090755812628634e-01,3.0769741250429189e-01,
     3.0450139197664983e-01,3.0131939610080288e-01,2.9815132669668531e-01,
     2.9499708779996164e-01,2.9185658561709499e-01,2.8872972848218270e-01,
     2.8561642681550159e-01,2.8251659308370741e-01,2.7943014176163772e-01,
     2.7635698929566810e-01,2.7329705406857691e-01,2.7025025636587519e-01,
     2.6721651834356114e-01,2.6419576399726080e-01,2.6118791913272082e-01,
     2.5819291133761890e-01,2.5521066995466168e-01,2.5224112605594190e-01,
     2.4928421241852824e-01,2.4633986350126363e-01,2.4340801542275012e-01,
     2.4048860594050039e-01,2.3758157443123795e-01,2.3468686187232990e-01,
     2.3180441082433859e-01,2.2893416541468023e-01,2.2607607132238020e-01,
     2.2323007576391746e-01,2.2039612748015194e-01,2.1757417672433113e-01,
     2.1476417525117358e-01,2.1196607630703015e-01,2.0917983462112499e-01,
     2.0640540639788071e-01,2.0364274931033485e-01,2.0089182249465656e-01,
     1.9815258654577511e-01,1.9542500351413428e-01,1.9270903690358912e-01,
     1.9000465167046496e-01,1.8731181422380025e-01,1.8463049242679927e-01,
     1.8196065559952254e-01,1.7930227452284767e-01,1.7665532144373500e-01,
     1.7401977008183875e-01,1.7139559563750595e-01,1.6878277480121151e-01,
     1.6618128576448205e-01,1.6359110823236570e-01,1.6101222343751107e-01,
     1.5844461415592431e-01,1.5588826472447920e-01,1.5334316106026283e-01,
     1.5080929068184568e-01,1.4828664273257453e-01,1.4577520800599403e-01,
     1.4327497897351341e-01,1.4078594981444470e-01,1.3830811644855071e-01,
     1.3584147657125373e-01,1.3338602969166913e-01,1.3094177717364430e-01,
     1.2850872227999952e-01,1.2608687022018586e-01,1.2367622820159654e-01,
     1.2127680548479021e-01,1.1888861344290998e-01,1.1651166562561080e-01,
     1.1414597782783835e-01,1.1179156816383801e-01,1.0944845714681163e-01,
     1.0711666777468364e-01,1.0479622562248690e-01,1.0248715894193508e-01,
     1.0018949876880981e-01,9.7903279038862284e-02,9.5628536713008819e-02,
     9.3365311912690860e-02,9.1113648066373634e-02,8.8873592068275789e-02,
     8.6645194450557961e-02,8.4428509570353374e-02,8.2223595813202863e-02,
     8.0030515814663056e-02,7.7849336702096039e-02,7.5680130358927067e-02,
     7.3522973713981268e-02,7.1377949058890375e-02,6.9245144397006769e-02,
     6.7124653827788497e-02,6.5016577971242842e-02,6.2921024437758113e-02,
     6.0838108349539864e-02,5.8767952920933758e-02,5.6710690106202902e-02,
     5.4666461324888914e-02,5.2635418276792176e-02,5.0617723860947761e-02,
     4.8613553215868521e-02,4.6623094901930368e-02,4.4646552251294443e-02,
     4.2684144916474431e-02,4.0736110655940933e-02,3.8802707404526113e-02,
     3.6884215688567284e-02,3.4980941461716084e-02,3.3093219458578522e-02,
     3.1221417191920245e-02,2.9365939758133314e-02,2.7527235669603082e-02,
     2.5705804008548896e-02,2.3902203305795882e-02,2.2117062707308864e-02,
     2.0351096230044517e-02,1.8605121275724643e-02,1.6880083152543166e-02,
     1.5177088307935325e-02,1.3497450601739880e-02,1.1842757857907888e-02,
     1.0214971439701471e-02,8.6165827693987316e-03,7.0508754713732268e-03,
     5.5224032992509968e-03,4.0379725933630305e-03,2.6090727461021627e-03,
     1.2602859304985975e-03]

## Tables for exponential variates
const ke =
Uint64[0x000e290a13924be4,0                 ,0x0009beadebce18c0,0x000c377ac71f9e08,
       0x000d4ddb99075857,0x000de893fb8ca23e,0x000e4a8e87c4328e,0x000e8dff16ae1cba,
       0x000ebf2deab58c5a,0x000ee49a6e8b9639,0x000f0204efd64ee5,0x000f19bdb8ea3c1c,
       0x000f2d458bbe5bd2,0x000f3da104b78236,0x000f4b86d784571f,0x000f577ad8a7784f,
       0x000f61de83da32ac,0x000f6afb7843cce7,0x000f730a57372b44,0x000f7a37651b0e68,
       0x000f80a5bb6eea52,0x000f867189d3cb5c,0x000f8bb1b4f8fbbd,0x000f9079062292b9,
       0x000f94d70ca8d43a,0x000f98d8c7dcaa99,0x000f9c8928abe083,0x000f9ff175b734a6,
       0x000fa319996bc47e,0x000fa6085f8e9d08,0x000fa8c3a62e1991,0x000fab5084e1f660,
       0x000fadb36c84cccb,0x000faff041086846,0x000fb20a6ea22bb9,0x000fb404fb42cb3d,
       0x000fb5e295158174,0x000fb7a59e99727a,0x000fb95038c8789d,0x000fbae44ba684ec,
       0x000fbc638d822e60,0x000fbdcf89209ffb,0x000fbf29a303cfc5,0x000fc0731df1089d,
       0x000fc1ad1ed6c8b1,0x000fc2d8b02b5c8a,0x000fc3f6c4d92131,0x000fc5083ac9ba7d,
       0x000fc60ddd1e9cd7,0x000fc7086622e825,0x000fc7f881009f0c,0x000fc8decb41ac71,
       0x000fc9bbd623d7ec,0x000fca9027c5b26e,0x000fcb5c3c319c49,0x000fcc20864b4449,
       0x000fccdd70a35d41,0x000fcd935e34bf80,0x000fce42ab0db8bd,0x000fceebace7ec02,
       0x000fcf8eb3b0d0e7,0x000fd02c0a049b60,0x000fd0c3f59d199d,0x000fd156b7b5e27e,
       0x000fd1e48d670342,0x000fd26daff73552,0x000fd2f2552684bf,0x000fd372af7233c2,
       0x000fd3eeee528f62,0x000fd4673e73543b,0x000fd4dbc9e72ff8,0x000fd54cb856dc2c,
       0x000fd5ba2f2c4119,0x000fd62451ba02c3,0x000fd68b415fcff5,0x000fd6ef1dabc161,
       0x000fd75004790eb7,0x000fd7ae120c583f,0x000fd809612dbd0a,0x000fd8620b40effa,
       0x000fd8b8285b78fe,0x000fd90bcf594b1d,0x000fd95d15efd426,0x000fd9ac10bfa70c,
       0x000fd9f8d364df06,0x000fda437086566c,0x000fda8bf9e3c9ff,0x000fdad28062fed5,
       0x000fdb17141bff2d,0x000fdb59c4648085,0x000fdb9a9fda83cd,0x000fdbd9b46e3ed4,
       0x000fdc170f6b5d05,0x000fdc52bd81a3fb,0x000fdc8ccacd07ba,0x000fdcc542dd3902,
       0x000fdcfc30bcb794,0x000fdd319ef77143,0x000fdd6597a0f60c,0x000fdd98245a48a3,
       0x000fddc94e575271,0x000fddf91e64014f,0x000fde279ce914cb,0x000fde54d1f0a06b,
       0x000fde80c52a47d0,0x000fdeab7def394e,0x000fded50345eb36,0x000fdefd5be59fa1,
       0x000fdf248e39b26f,0x000fdf4aa064b4b0,0x000fdf6f98435894,0x000fdf937b6f30bb,
       0x000fdfb64f414572,0x000fdfd818d48262,0x000fdff8dd07fed8,0x000fe018a08122c5,
       0x000fe03767adaa5a,0x000fe05536c58a14,0x000fe07211ccb4c5,0x000fe08dfc94c532,
       0x000fe0a8fabe8ca2,0x000fe0c30fbb87a6,0x000fe0dc3ecf3a5a,0x000fe0f48b107522,
       0x000fe10bf76a82ef,0x000fe122869e4200,0x000fe1383b4327e1,0x000fe14d17c83188,
       0x000fe1611e74c023,0x000fe1745169635a,0x000fe186b2a09177,0x000fe19843ef4e08,
       0x000fe1a90705bf64,0x000fe1b8fd6fb37c,0x000fe1c828951444,0x000fe1d689ba4bfd,
       0x000fe1e4220099a5,0x000fe1f0f26655a0,0x000fe1fcfbc726d4,0x000fe2083edc2831,
       0x000fe212bc3bfeb4,0x000fe21c745adfe4,0x000fe225678a8895,0x000fe22d95fa23f4,
       0x000fe234ffb62282,0x000fe23ba4a800d9,0x000fe2418495fddd,0x000fe2469f22bffb,
       0x000fe24af3cce90e,0x000fe24e81ee9859,0x000fe25148bcda1a,0x000fe253474703fe,
       0x000fe2547c75fdc6,0x000fe254e70b7550,0x000fe25485a0fd1a,0x000fe25356a71450,
       0x000fe2515864173b,0x000fe24e88f316f2,0x000fe24ae64296fb,0x000fe2466e132f61,
       0x000fe2411df611bd,0x000fe23af34b6f73,0x000fe233eb40bf42,0x000fe22c02cee01c,
       0x000fe22336b81711,0x000fe2198385e5cd,0x000fe20ee586b707,0x000fe20358cb5dfc,
       0x000fe1f6d92465b1,0x000fe1e9621f2c9e,0x000fe1daef02c8da,0x000fe1cb7accb0a6,
       0x000fe1bb002d22ca,0x000fe1a9798349b9,0x000fe196e0d9140d,0x000fe1832fdebc44,
       0x000fe16e5fe5f932,0x000fe15869dccfcf,0x000fe1414647fe78,0x000fe128ed3cf8b2,
       0x000fe10f565b69cf,0x000fe0f478c633ab,0x000fe0d84b1bdd9e,0x000fe0bac36e6688,
       0x000fe09bd73a6b5c,0x000fe07b7b5d920b,0x000fe059a40c26d2,0x000fe03644c5d7f9,
       0x000fe011504979b3,0x000fdfeab887b95d,0x000fdfc26e94a448,0x000fdf986297e306,
       0x000fdf6c83bb8663,0x000fdf3ec0193eee,0x000fdf0f04a5d30a,0x000fdedd3d1aa204,
       0x000fdea953dcfc13,0x000fde7331e3100e,0x000fde3abe9626f3,0x000fddffdfb1dbd5,
       0x000fddc2791ff351,0x000fdd826cd068c7,0x000fdd3f9a8d3857,0x000fdcf9dfc95b0d,
       0x000fdcb1176a55fe,0x000fdc65198ba50c,0x000fdc15bb3b2daa,0x000fdbc2ce2dc4ae,
       0x000fdb6c206aaaca,0x000fdb117becb4a2,0x000fdab2a6379bf1,0x000fda4f5fdfb4e9,
       0x000fd9e76401f3a4,0x000fd97a67a9ce20,0x000fd9081922142a,0x000fd8901f2d4b02,
       0x000fd812182170e1,0x000fd78d98e23cd4,0x000fd7022bb3f083,0x000fd66f4edf96ba,
       0x000fd5d473200306,0x000fd530f9ccff94,0x000fd48432b7b351,0x000fd3cd59a8469f,
       0x000fd30b9368f90a,0x000fd23dea45f500,0x000fd16349e2e04b,0x000fd07a7a3ef98b,
       0x000fcf8219b5df06,0x000fce7895bcfcdf,0x000fcd5c220ad5e3,0x000fcc2aadbc17dd,
       0x000fcae1d5e81fbd,0x000fc97ed4e778fa,0x000fc7fe6d4d720f,0x000fc65ccf39c2fc,
       0x000fc4957623cb04,0x000fc2a2fc826dc8,0x000fc07ee19b01ce,0x000fbe213c1cf493,
       0x000fbb8051ac1566,0x000fb890078d120e,0x000fb5411a5b9a96,0x000fb18000547134,
       0x000fad334827f1e2,0x000fa839276708b9,0x000fa263b32e37ee,0x000f9b72d1c52cd1,
       0x000f930a1a281a05,0x000f889f023d820a,0x000f7b577d2be5f4,0x000f69c650c40a8f,
       0x000f51530f0916d9,0x000f2cb0e3c5933e,0x000eeefb15d605d9,0x000e6da6ecf27460]

const we =
    [1.9311480126418366e-15,1.4178028487910829e-17,2.3278824993382448e-17,
     3.0487830247064320e-17,3.6665697714474878e-17,4.2179302189289733e-17,
     4.7222561556862764e-17,5.1911915446217879e-17,5.6323471083955047e-17,
     6.0510082606427647e-17,6.4510165096727506e-17,6.8352646803700541e-17,
     7.2059939574689050e-17,7.5649815537392981e-17,7.9136643961951065e-17,
     8.2532235563518929e-17,8.5846436168850513e-17,8.9087554865647428e-17,
     9.2262679629663719e-17,9.5377914505292719e-17,9.8438560874559257e-17,
     1.0144925809006294e-16,1.0441409405585343e-16,1.0733669323436384e-16,
     1.1022028745670189e-16,1.1306777346479334e-16,1.1588176009705533e-16,
     1.1866460730417886e-16,1.2141845865694359e-16,1.2414526862326387e-16,
     1.2684682560606153e-16,1.2952477151912284e-16,1.3218061851538810e-16,
     1.3481576335745444e-16,1.3743149982367625e-16,1.4002902946807859e-16,
     1.4260947099321287e-16,1.4517386844829297e-16,1.4772319842763584e-16,
     1.5025837641447456e-16,1.5278026239101652e-16,1.5528966581595696e-16,
     1.5778735005459581e-16,1.6027403633350909e-16,1.6275040728083524e-16,
     1.6521711010420076e-16,1.6767475945078279e-16,1.7012393998770646e-16,
     1.7256520873568226e-16,1.7499909718432365e-16,1.7742611321380505e-16,
     1.7984674284430714e-16,1.8226145183195818e-16,1.8467068712763576e-16,
     1.8707487821298258e-16,1.8947443832625899e-16,1.9186976558915995e-16,
     1.9426124404443042e-16,1.9664924461299023e-16,1.9903412597830144e-16,
     2.0141623540485899e-16,2.0379590949693882e-16,2.0617347490308439e-16,
     2.0854924897123771e-16,2.1092354035891528e-16,2.1329664960238294e-16,
     2.1566886964838970e-16,2.1804048635167009e-16,2.2041177894111562e-16,
     2.2278302045723950e-16,2.2515447816331350e-16,2.2752641393233694e-16,
     2.2989908461180186e-16,2.3227274236804366e-16,2.3464763501180916e-16,
     2.3702400630653389e-16,2.3940209626069303e-16,2.4178214140547710e-16,
     2.4416437505894123e-16,2.4654902757768304e-16,2.4893632659702250e-16,
     2.5132649726057970e-16,2.5371976244007951e-16,2.5611634294614988e-16,
     2.5851645773082391e-16,2.6092032408240577e-16,2.6332815781331452e-16,
     2.6574017344147618e-16,2.6815658436579989e-16,2.7057760303623509e-16,
     2.7300344111887955e-16,2.7543430965657619e-16,2.7787041922541278e-16,
     2.8031198008751431e-16,2.8275920234049704e-16,2.8521229606393309e-16,
     2.8767147146315804e-16,2.9013693901073754e-16,2.9260890958589514e-16,
     2.9508759461219033e-16,2.9757320619372521e-16,3.0006595725014739e-16,
     3.0256606165070789e-16,3.0507373434762511e-16,3.0758919150899939e-16,
     3.1011265065151543e-16,3.1264433077316750e-16,3.1518445248623523e-16,
     3.1773323815073683e-16,3.2029091200858335e-16,3.2285770031865573e-16,
     3.2543383149302610e-16,3.2801953623454359e-16,3.3061504767600738e-16,
     3.3322060152114841e-16,3.3583643618764577e-16,3.3846279295240445e-16,
     3.4109991609932597e-16,3.4374805306980633e-16,3.4640745461620167e-16,
     3.4907837495850680e-16,3.5176107194449828e-16,3.5445580721360130e-16,
     3.5716284636474652e-16,3.5988245912849274e-16,3.6261491954370031e-16,
     3.6536050613905045e-16,3.6811950211971757e-16,3.7089219555951389e-16,
     3.7367887959883854e-16,3.7647985264877841e-16,3.7929541860172334e-16,
     3.8212588704887531e-16,3.8497157350504876e-16,3.8783279964117988e-16,
     3.9070989352498183e-16,3.9360318987020748e-16,3.9651303029500381e-16,
     3.9943976358986842e-16,4.0238374599574693e-16,4.0534534149283966e-16,
     4.0832492210071775e-16,4.1132286819038357e-16,4.1433956880894741e-16,
     4.1737542201763194e-16,4.2043083524385856e-16,4.2350622564821518e-16,
     4.2660202050715582e-16,4.2971865761233266e-16,4.3285658568752094e-16,
     4.3601626482415681e-16,4.3919816693657415e-16,4.4240277623809919e-16,
     4.4563058973923611e-16,4.4888211776926172e-16,4.5215788452263475e-16,
     4.5545842863172421e-16,4.5878430376746227e-16,4.6213607926964266e-16,
     4.6551434080870692e-16,4.6891969108099157e-16,4.7235275053955480e-16,
     4.7581415816285534e-16,4.7930457226372470e-16,4.8282467134125866e-16,
     4.8637515497845119e-16,4.8995674478861404e-16,4.9357018541385775e-16,
     4.9721624557917034e-16,5.0089571920591141e-16,5.0460942658884340e-16,
     5.0835821564116245e-16,5.1214296321235415e-16,5.1596457648410618e-16,
     5.1982399444994938e-16,5.2372218948478484e-16,5.2766016901098856e-16,
     5.3163897726836902e-16,5.3565969719590503e-16,5.3972345243389779e-16,
     5.4383140945596370e-16,5.4798477984116296e-16,5.5218482269752343e-16,
     5.5643284724928722e-16,5.6073021560139669e-16,5.6507834569605064e-16,
     5.6947871447763482e-16,5.7393286128396354e-16,5.7844239148359912e-16,
     5.8300898038105864e-16,5.8763437741400573e-16,5.9232041066909314e-16,
     5.9706899174600906e-16,6.0188212100252363e-16,6.0676189321700068e-16,
     6.1171050370897217e-16,6.1673025496306200e-16,6.2182356380685327e-16,
     6.2699296919933262e-16,6.3224114069342115e-16,6.3757088764394262e-16,
     6.4298516924135947e-16,6.4848710546189033e-16,6.5407998903644809e-16,
     6.5976729855445663e-16,6.6555271283433428e-16,6.7144012671064882e-16,
     6.7743366840910103e-16,6.8353771870512740e-16,6.8975693209068478e-16,
     6.9609626020748846e-16,7.0256097784459588e-16,7.0915671184495837e-16,
     7.1588947332085531e-16,7.2276569364381212e-16,7.2979226475290851e-16,
     7.3697658441912426e-16,7.4432660721604146e-16,7.5185090208325131e-16,
     7.5955871753377488e-16,7.6746005575784274e-16,7.7556575712157906e-16,
     7.8388759686228577e-16,7.9243839615735500e-16,8.0123215021130834e-16,
     8.1028417659131464e-16,8.1961128778061250e-16,8.2923199285818092e-16,
     8.3916673441467979e-16,8.4943816836487701e-16,8.6007149633349414e-16,
     8.7109486293879040e-16,8.8253983380721398e-16,8.9444197485198646e-16,
     9.0684155971316690e-16,9.1978444098118649e-16,9.3332313294229516e-16,
     9.4751817065249841e-16,9.6243983456584759e-16,9.7817036547844198e-16,
     9.9480684723838795e-16,1.0124650144288319e-15,1.0312843657756166e-15,
     1.0514351604044550e-15,1.0731281954224043e-15,1.0966288068517408e-15,
     1.1222774909350319e-15,1.1505212963006663e-15,1.1819635283304206e-15,
     1.2174462832361815e-15,1.2581958069755114e-15,1.3060984107128082e-15,
     1.3642786158057857e-15,1.4384889932178723e-15,1.5412190700064194e-15,
     1.7091034077168055e-15]
const fe =
    [1.0000000000000000e+00,9.3814368086217470e-01,9.0046992992574648e-01,
     8.7170433238120359e-01,8.4778550062398961e-01,8.2699329664305032e-01,
     8.0842165152300838e-01,7.9152763697249562e-01,7.7595685204011555e-01,
     7.6146338884989628e-01,7.4786862198519510e-01,7.3503809243142348e-01,
     7.2286765959357202e-01,7.1127476080507601e-01,7.0019265508278816e-01,
     6.8956649611707799e-01,6.7935057226476536e-01,6.6950631673192473e-01,
     6.6000084107899970e-01,6.5080583341457110e-01,6.4189671642726609e-01,
     6.3325199421436607e-01,6.2485273870366598e-01,6.1668218091520766e-01,
     6.0872538207962201e-01,6.0096896636523223e-01,5.9340090169173343e-01,
     5.8601031847726803e-01,5.7878735860284503e-01,5.7172304866482582e-01,
     5.6480919291240017e-01,5.5803828226258745e-01,5.5140341654064129e-01,
     5.4489823767243961e-01,5.3851687200286191e-01,5.3225388026304332e-01,
     5.2610421398361973e-01,5.2006317736823360e-01,5.1412639381474856e-01,
     5.0828977641064288e-01,5.0254950184134772e-01,4.9690198724154955e-01,
     4.9134386959403253e-01,4.8587198734188491e-01,4.8048336393045421e-01,
     4.7517519303737737e-01,4.6994482528395998e-01,4.6478975625042618e-01,
     4.5970761564213769e-01,4.5469615747461550e-01,4.4975325116275500e-01,
     4.4487687341454851e-01,4.4006510084235390e-01,4.3531610321563657e-01,
     4.3062813728845883e-01,4.2599954114303434e-01,4.2142872899761658e-01,
     4.1691418643300288e-01,4.1245446599716118e-01,4.0804818315203240e-01,
     4.0369401253053028e-01,3.9939068447523107e-01,3.9513698183329016e-01,
     3.9093173698479711e-01,3.8677382908413765e-01,3.8266218149600983e-01,
     3.7859575940958079e-01,3.7457356761590216e-01,3.7059464843514600e-01,
     3.6665807978151416e-01,3.6276297335481777e-01,3.5890847294874978e-01,
     3.5509375286678746e-01,3.5131801643748334e-01,3.4758049462163698e-01,
     3.4388044470450241e-01,3.4021714906678002e-01,3.3658991402867761e-01,
     3.3299806876180899e-01,3.2944096426413633e-01,3.2591797239355619e-01,
     3.2242848495608917e-01,3.1897191284495724e-01,3.1554768522712895e-01,
     3.1215524877417955e-01,3.0879406693456019e-01,3.0546361924459026e-01,
     3.0216340067569353e-01,2.9889292101558179e-01,2.9565170428126120e-01,
     2.9243928816189257e-01,2.8925522348967775e-01,2.8609907373707683e-01,
     2.8297041453878075e-01,2.7986883323697292e-01,2.7679392844851736e-01,
     2.7374530965280297e-01,2.7072259679906002e-01,2.6772541993204479e-01,
     2.6475341883506220e-01,2.6180624268936298e-01,2.5888354974901623e-01,
     2.5598500703041538e-01,2.5311029001562946e-01,2.5025908236886230e-01,
     2.4743107566532763e-01,2.4462596913189211e-01,2.4184346939887721e-01,
     2.3908329026244918e-01,2.3634515245705964e-01,2.3362878343743335e-01,
     2.3093391716962741e-01,2.2826029393071670e-01,2.2560766011668407e-01,
     2.2297576805812019e-01,2.2036437584335949e-01,2.1777324714870053e-01,
     2.1520215107537868e-01,2.1265086199297828e-01,2.1011915938898826e-01,
     2.0760682772422204e-01,2.0511365629383771e-01,2.0263943909370902e-01,
     2.0018397469191127e-01,1.9774706610509887e-01,1.9532852067956322e-01,
     1.9292814997677135e-01,1.9054576966319539e-01,1.8818119940425432e-01,
     1.8583426276219711e-01,1.8350478709776746e-01,1.8119260347549629e-01,
     1.7889754657247831e-01,1.7661945459049488e-01,1.7435816917135349e-01,
     1.7211353531532006e-01,1.6988540130252766e-01,1.6767361861725019e-01,
     1.6547804187493600e-01,1.6329852875190182e-01,1.6113493991759203e-01,
     1.5898713896931421e-01,1.5685499236936523e-01,1.5473836938446808e-01,
     1.5263714202744286e-01,1.5055118500103989e-01,1.4848037564386679e-01,
     1.4642459387834494e-01,1.4438372216063478e-01,1.4235764543247220e-01,
     1.4034625107486245e-01,1.3834942886358020e-01,1.3636707092642886e-01,
     1.3439907170221363e-01,1.3244532790138752e-01,1.3050573846833077e-01,
     1.2858020454522817e-01,1.2666862943751067e-01,1.2477091858083096e-01,
     1.2288697950954514e-01,1.2101672182667483e-01,1.1916005717532768e-01,
     1.1731689921155557e-01,1.1548716357863353e-01,1.1367076788274431e-01,
     1.1186763167005630e-01,1.1007767640518538e-01,1.0830082545103380e-01,
     1.0653700405000166e-01,1.0478613930657017e-01,1.0304816017125772e-01,
     1.0132299742595363e-01,9.9610583670637132e-02,9.7910853311492199e-02,
     9.6223742550432798e-02,9.4549189376055859e-02,9.2887133556043541e-02,
     9.1237516631040155e-02,8.9600281910032858e-02,8.7975374467270218e-02,
     8.6362741140756913e-02,8.4762330532368119e-02,8.3174093009632383e-02,
     8.1597980709237419e-02,8.0033947542319905e-02,7.8481949201606421e-02,
     7.6941943170480503e-02,7.5413888734058410e-02,7.3897746992364746e-02,
     7.2393480875708738e-02,7.0901055162371829e-02,6.9420436498728755e-02,
     6.7951593421936601e-02,6.6494496385339774e-02,6.5049117786753749e-02,
     6.3615431999807334e-02,6.2193415408540995e-02,6.0783046445479633e-02,
     5.9384305633420266e-02,5.7997175631200659e-02,5.6621641283742877e-02,
     5.5257689676697037e-02,5.3905310196046087e-02,5.2564494593071692e-02,
     5.1235237055126281e-02,4.9917534282706372e-02,4.8611385573379497e-02,
     4.7316792913181548e-02,4.6033761076175170e-02,4.4762297732943282e-02,
     4.3502413568888183e-02,4.2254122413316234e-02,4.1017441380414819e-02,
     3.9792391023374125e-02,3.8578995503074857e-02,3.7377282772959361e-02,
     3.6187284781931423e-02,3.5009037697397410e-02,3.3842582150874330e-02,
     3.2687963508959535e-02,3.1545232172893609e-02,3.0414443910466604e-02,
     2.9295660224637393e-02,2.8188948763978636e-02,2.7094383780955800e-02,
     2.6012046645134217e-02,2.4942026419731783e-02,2.3884420511558171e-02,
     2.2839335406385240e-02,2.1806887504283581e-02,2.0787204072578117e-02,
     1.9780424338009743e-02,1.8786700744696030e-02,1.7806200410911362e-02,
     1.6839106826039948e-02,1.5885621839973163e-02,1.4945968011691148e-02,
     1.4020391403181938e-02,1.3109164931254991e-02,1.2212592426255381e-02,
     1.1331013597834597e-02,1.0464810181029979e-02,9.6144136425022099e-03,
     8.7803149858089753e-03,7.9630774380170400e-03,7.1633531836349839e-03,
     6.3819059373191791e-03,5.6196422072054830e-03,4.8776559835423923e-03,
     4.1572951208337953e-03,3.4602647778369040e-03,2.7887987935740761e-03,
     2.1459677437189063e-03,1.5362997803015724e-03,9.6726928232717454e-04,
     4.5413435384149677e-04]

ziggurat_nor_r      = 3.6541528853610087963519472518
ziggurat_nor_inv_r  = inv(ziggurat_nor_r)
ziggurat_exp_r      = 7.6971174701310497140446280481

randi() = reinterpret(Uint64,dsfmt_gv_genrand_close1_open2()) & 0x000fffffffffffff
randi(state::DSFMT_state) = reinterpret(Uint64,dsfmt_genrand_close1_open2(state)) & 0x000fffffffffffff
for (lhs, rhs) in (([], []),
                  ([:(state::DSFMT_state)], [:state]))
    @eval begin
        function randmtzig_randn($(lhs...))
            @inbounds begin
                while true
                    r = randi($(rhs...))
                    rabs = int64(r>>1) # One bit for the sign
                    idx = rabs & 0xFF
                    x = (r&1 != 0x000000000 ? -rabs : rabs)*wi[idx+1]
                    if rabs < ki[idx+1]
                        return x # 99.3% of the time we return here 1st try
                    elseif idx == 0
                        while true
                            xx = -$(ziggurat_nor_inv_r)*log(rand($(rhs...)))
                            yy = -log(rand($(rhs...)))
                            if yy+yy > xx*xx
                                return (rabs & 0x100) != 0x000000000 ? -$(ziggurat_nor_r)-xx : $(ziggurat_nor_r)+xx
                            end
                        end
                    elseif (fi[idx] - fi[idx+1])*rand($(rhs...)) + fi[idx+1] < exp(-0.5*x*x)
                        return x # return from the triangular area
                    end
                end
            end
        end

        function randmtzig_exprnd($(lhs...))
            @inbounds begin
                while true
                    ri = randi($(rhs...))
                    idx = ri & 0xFF
                    x = ri*we[idx+1]
                    if ri < ke[idx+1]
                        return x # 98.9% of the time we return here 1st try
                    elseif idx == 0
                        return $(ziggurat_exp_r) - log(rand($(rhs...)))
                    elseif (fe[idx] - fe[idx+1])*rand($(rhs...)) + fe[idx+1] < exp(-x)
                        return x # return from the triangular area
                    end
                end
            end
        end
    end
end
## Windows entropy

@windows_only begin
    function win32_SystemFunction036!(a::Array{Uint32})
        ccall((:SystemFunction036,:Advapi32),stdcall,Uint8,(Ptr{Void},Uint32),a,length(a)*sizeof(eltype(a)))
    end
end

end # module

max_stars_repo_name: alanedelman/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
function arc_solution()
  cur_T = 14.0

  cur_reactor = arc_deck(cur_T)

  cur_reactor.n_bar = 1.3
  cur_reactor.I_P = 7.8

  cur_reactor.R_0 = 3.3
  cur_reactor.B_0 = 9.2

  cur_reactor.P_F = 525

  cur_reactor.f_BS = 0.63
  cur_reactor.f_CD = 1.00 - cur_reactor.f_BS
  cur_reactor.f_IN = 0.00

  cur_reactor.beta_N = cur_reactor.max_beta_N
  cur_reactor.q_95 = 7.2
  cur_reactor.P_W = 2.5

  cur_reactor.volume = 141

  cur_reactor
end


text:
__precompile__()
module Corona

using Dates
using TimeSeries
using DataStructures

include("utils.jl")
include("DA.jl")
include("dataio.jl")

end # module Corona


text:
using BEAST
using CompScienceMeshes

using Base.Test

Γ = meshrectangle(1.0,1.0,1.0)
γ = meshsegment(1.0,1.0,3)
X = raviartthomas(Γ, γ)

x = divergence(X)

@test numfunctions(X) == 2
@test numfunctions(x) == numfunctions(X)
for (f,g) in zip(x.fns, X.fns)
    @test length(f) == length(g)
    if length(f) == 1
        c = f[1].cellid
        s = chart(Γ, Γ.faces[c])
        @test f[1].coeff == 1 / volume(s)
    else
        @test f[1].coeff + f[2].coeff == 0
    end
end

edges = skeleton(Γ, 1)
vps = cellpairs(Γ, edges)
rwg = raviartthomas(Γ, vps)
dvg = divergence(rwg)

D = spzeros(Int,numcells(edges), numcells(Γ))
for (m,f) in enumerate(dvg.fns)
    for s in f
        e = s.cellid
        D[m,e] = s.coeff
    end
end

C = connectivity(edges, Γ)
@test size(C) == size(D')

## Test the curl of the linear Lagrange elements
using CompScienceMeshes
using BEAST
using Base.Test

m = meshrectangle(1.0, 1.0, 0.25)
X = duallagrangec0d1(m)
@test numfunctions(X) == numcells(m)

Y = curl(X)
@test numfunctions(X) == numfunctions(Y)
@test isa(Y, BEAST.RTBasis)
@test isa(refspace(Y), BEAST.BEAST.RTRefSpace)

# test the chain property
Z = divergence(Y)
for fn ∈ Z.fns
    total = zero(scalartype(Z))
    for sh ∈ fn
        total += sh.coeff
    end
    @test total ≈ 0
end


text:
floats = [
    5e-324,
    1.7976931348623157e308,
    4294967272.0,
    4.1855804968213567e298,
    5.5626846462680035e-309,
    2147483648.0,
    3.5844466002796428e+298,
    reinterpret(Float64, 0x0010000000000000),
    reinterpret(Float64, 0x000FFFFFFFFFFFFF),
    Float32(1e-45),
    3.4028234f38,
    Float32(4294967272.0),
    Float32(3.32306998946228968226e+35),
    Float32(1.2341e-41),
    Float32(3.3554432e7),
    Float32(3.26494756798464e14),
    Float32(3.91132223637771935344e37),
    reinterpret(Float32, 0x00800000),
    reinterpret(Float32, 0x007FFFFF),
    1.0,
    1.5,
    # 3.3161339052167390562200598e-237,
    7.9885183916008099497815232e+191,
    Float64(0xFFFFFFFF),
    1e21,
    999999999999999868928.00,
    6.9999999999999989514240000e+21,
    1.55,
    1.00000001,
    0.1,
    0.01,
    0.001,
    0.0001,
    0.00001,
    0.000001,
    0.0000001,
    0.00000001,
    0.000000001,
    0.0000000001,
    0.00000000001,
    0.000000000001,
    0.0000000000001,
    0.00000000000001,
    0.000000000000001,
    0.0000000000000001,
    0.00000000000000001,
    0.000000000000000001,
    0.0000000000000000001,
    0.00000000000000000001,
    0.10000000004,
    0.00010000004,
    0.00001000004,
    0.00000100004,
    0.00000010004,
    0.00000001004,
    0.00000000104,
    0.0000000001000004,
    0.0000000000100004,
    0.0000000000010004,
    0.0000000000001004,
    0.0000000000000104,
    0.6,
    0.96,
    0.996,
    0.9996,
    0.99996,
    0.999996,
    0.9999996,
    0.99999996,
    0.999999996,
    0.9999999996,
    0.99999999996,
    0.999999999996,
    0.9999999999996,
    0.99999999999996,
    0.999999999999996,
    0.9999999999999996,
    0.00999999999999996,
    0.000999999999999996,
    0.0000999999999999996,
    0.00000999999999999996,
    0.000000999999999999996,
    323423.234234,
    12345678.901234,
    98765.432109,
    42,
    0.5,
    1e-23,
    1e-123,
    1e-123,
    1e-23,
    1e-21,
    1e-22,
    6e-21,
    9.1193616301674545152000000e+19,
    4.8184662102767651659096515e-04,
    1.9023164229540652612705182e-23,
    1000000000000000128.0,
    floatmin(Float16),
    floatmax(Float16),
    Float16(1.0),
    Float16(1.5),
    Float16(1.55),
    Float16(1.00000001),
    Float16(0.1),
    Float16(0.01),
    Float16(0.001),
    Float16(0.0001),
    Float16(0.00001),
    Float16(0.000001),
    Float16(0.0000001),
    Float16(0.6),
    Float16(0.96),
    Float16(0.996),
    Float16(0.9996),
    Float16(0.99996),
    Float16(0.999996),
    Float16(0.9999996),
    Float16(0.99999996),
    Float16(42),
    Float16(0.5),
]
for f in floats
    try
        @test JSON2.read(JSON2.write(f), typeof(f)) == f
    catch
        println("parsing $(typeof(f)) '$f' failed")
    end
end

for F in (Float16, Float32, Float64)
    for f in rand(F, 100)
        try
            @test JSON2.read(JSON2.write(f), F) == f
        catch
            println("parsing $(typeof(f)) '$f' failed")
        end
    end
end