text:
@inline @fastmath Ca(q::Float64, p::Float64, costh::Float64) =
    (
        4 *
        (1 - costh^2) *
        (
            p * q * costh * (p * q * costh - 6 * (p^2 + q^2)) +
            8 * p^2 * q^2 +
            3 * p^4 +
            3 * q^4
        )
    ) * Qmp(q, p, costh)^-1

Ca(q::Float64, p::Vector{Float64}, costh::Float64) = @avx (@. 4 *
         (1 - costh^2) *
         (
             p * q * costh * (p * q * costh - 6 * (p^2 + q^2)) +
             8 * p^2 * q^2 +
             3 * p^4 +
             3 * q^4
         )) .* Qmp(q, p, costh) .^ -1

@inline @fastmath Cg(q::Float64, costh::Float64) = q^2 * (1 - costh^2)

@inline @fastmath Ct(costh::Float64) = 2.0 * costh^2 - 14.0

@inline @fastmath Ccbc(q::Float64, p::Float64, costh::Float64) =
    p^2 * q^2 * (1 - costh^2) * Qmp(q, p, costh)^-1

Ccbc(q::Float64, p::Vector{Float64}, costh::Float64) =
    @avx p .^ 2 .* q^2 .* (1 - costh^2) .* Qmp(q, p, costh) .^ -1

@inline @fastmath Cpcbc(p::Float64, costh::Float64) = p^2 * (1 - costh^2)

Cpcbc(p::Vector{Float64}, costh::Float64) = @. p^2 * (1 - costh^2)

@inline @fastmath GGG1(q::Float64, p::Float64, x1::Float64, x2::Float64) =
    (
        p^3 * (
            2 * p * (-1 + x1^2 + x1 * sqrt(3 - 3 * x1^2) * x2) +
            q * (
                -(sqrt(3 - 3 * x1^2) * x2) +
                2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 -
                x1 * (-1 + x1^2) * (-1 + 3 * x2^2)
            )
        )
    ) / 4

GGG1(q::Float64, p::Vector{Float64}, x1::Float64, x2::Float64) = @avx @. (
    p^3 * (
        2 * p * (-1 + x1^2 + x1 * sqrt(3 - 3 * x1^2) * x2) +
        q * (
            -(sqrt(3 - 3 * x1^2) * x2) + 2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 -
            x1 * (-1 + x1^2) * (-1 + 3 * x2^2)
        )
    )
) / 4

@inline @fastmath GGG2(q::Float64, p::Float64, x1::Float64, x2::Float64) =
    -(
        p^3 *
        q *
        (
            p^2 * sqrt(3 - 3 * x1^2) * x2 -
            p * q * (-1 + x1^2) * (1 + 3 * x2^2) +
            q^2 * (
                x1 + sqrt(3 - 3 * x1^2) * x2 -
                2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 - 3 * x1 * x2^2 +
                x1^3 * (-1 + 3 * x2^2)
            )
        )
    ) / (4 * (p^2 + q^2 - 2 * p * q * x1))
GGG2(q::Float64, p::Vector{Float64}, x1::Float64, x2::Float64) = @avx @. -(
    p^3 *
    q *
    (
        p^2 * sqrt(3 - 3 * x1^2) * x2 - p * q * (-1 + x1^2) * (1 + 3 * x2^2) +
        q^2 * (
            x1 + sqrt(3 - 3 * x1^2) * x2 - 2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 - 3 * x1 * x2^2 +
            x1^3 * (-1 + 3 * x2^2)
        )
    )
) / (4 * (p^2 + q^2 - 2 * p * q * x1))

@inline @fastmath GGG3(q::Float64, p::Float64, x1::Float64, x2::Float64) =
    -(
        p^3 *
        q^2 *
        (
            p * (
                2 * x1 * sqrt(3 - 3 * x1^2) * x2 - 3 * x2^2 +
                3 * x1^2 * (1 + x2^2)
            ) +
            2 *
            q *
            (
                x1 + sqrt(3 - 3 * x1^2) * x2 -
                2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 - 3 * x1 * x2^2 +
                x1^3 * (-1 + 3 * x2^2)
            )
        )
    ) / (8 * (p^2 + q^2 - p * q * (x1 + sqrt(3 - 3 * x1^2) * x2)))
GGG3(q::Float64, p::Vector{Float64}, x1::Float64, x2::Float64) = @avx @. -(
    p^3 *
    q^2 *
    (
        p * (
            2 * x1 * sqrt(3 - 3 * x1^2) * x2 - 3 * x2^2 + 3 * x1^2 * (1 + x2^2)
        ) +
        2 *
        q *
        (
            x1 + sqrt(3 - 3 * x1^2) * x2 - 2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 - 3 * x1 * x2^2 +
            x1^3 * (-1 + 3 * x2^2)
        )
    )
) / (8 * (p^2 + q^2 - p * q * (x1 + sqrt(3 - 3 * x1^2) * x2)))

@inline @fastmath GGS1(q::Float64, p::Float64, x1::Float64, x2::Float64) = (
    (
        p^3 *
        q^2 *
        (
            3 *
            p^3 *
            (-4 + 3 * x1 * sqrt(3 - 3 * x1^2) * x2 + x2^2 - x1^2 * (-4 + x2^2)) +
            4 *
            q^3 *
            (
                -(sqrt(3 - 3 * x1^2) * x2) +
                2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 -
                x1 * (-1 + x1^2) * (-1 + 3 * x2^2)
            ) +
            2 *
            p *
            q^2 *
            (
                -9 +
                3 * x2^2 +
                x1^4 * (3 - 15 * x2^2) +
                x1^3 * sqrt(3 - 3 * x1^2) * x2 * (7 - 3 * x2^2) +
                6 * x1^2 * (1 + 2 * x2^2) +
                x1 * sqrt(3 - 3 * x1^2) * x2 * (-1 + 3 * x2^2)
            ) +
            p^2 *
            q *
            (
                x1^3 * (23 - 27 * x2^2) -
                6 * x1^2 * sqrt(3 - 3 * x1^2) * x2 * (-5 + x2^2) +
                sqrt(3 - 3 * x1^2) * x2 * (-13 + 6 * x2^2) +
                x1 * (-23 + 27 * x2^2)
            )
        )
    ) / (
        8 *
        (p^2 + q^2 + 2 * p * q * x1) *
        (p^2 + q^2 + p * q * (x1 + sqrt(3 - 3 * x1^2) * x2))
    )
)

GGS1(q::Float64, p::Vector{Float64}, x1::Float64, x2::Float64) = @avx @. (
    (
        p^3 *
        q^2 *
        (
            3 *
            p^3 *
            (-4 + 3 * x1 * sqrt(3 - 3 * x1^2) * x2 + x2^2 - x1^2 * (-4 + x2^2)) +
            4 *
            q^3 *
            (
                -(sqrt(3 - 3 * x1^2) * x2) +
                2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 -
                x1 * (-1 + x1^2) * (-1 + 3 * x2^2)
            ) +
            2 *
            p *
            q^2 *
            (
                -9 +
                3 * x2^2 +
                x1^4 * (3 - 15 * x2^2) +
                x1^3 * sqrt(3 - 3 * x1^2) * x2 * (7 - 3 * x2^2) +
                6 * x1^2 * (1 + 2 * x2^2) +
                x1 * sqrt(3 - 3 * x1^2) * x2 * (-1 + 3 * x2^2)
            ) +
            p^2 *
            q *
            (
                x1^3 * (23 - 27 * x2^2) -
                6 * x1^2 * sqrt(3 - 3 * x1^2) * x2 * (-5 + x2^2) +
                sqrt(3 - 3 * x1^2) * x2 * (-13 + 6 * x2^2) +
                x1 * (-23 + 27 * x2^2)
            )
        )
    ) / (
        8 *
        (p^2 + q^2 + 2 * p * q * x1) *
        (p^2 + q^2 + p * q * (x1 + sqrt(3 - 3 * x1^2) * x2))
    )
)

@inline @fastmath GGS2(q::Float64, p::Float64, x1::Float64, x2::Float64) =
    (
        p^3 * (
            p^3 * (-9 + 9 * x1^2 - 5 * x1 * sqrt(3 - 3 * x1^2) * x2) +
            4 *
            q^3 *
            (
                -(sqrt(3 - 3 * x1^2) * x2) +
                2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 -
                x1 * (-1 + x1^2) * (-1 + 3 * x2^2)
            ) -
            2 *
            p *
            q^2 *
            (
                7 +
                3 * x2^2 +
                x1 * sqrt(3 - 3 * x1^2) * x2 * (1 - 3 * x2^2) +
                x1^4 * (1 + 3 * x2^2) +
                x1^3 * sqrt(3 - 3 * x1^2) * x2 * (1 + 3 * x2^2) -
                2 * x1^2 * (4 + 3 * x2^2)
            ) +
            p^2 *
            q *
            (
                -11 * sqrt(3 - 3 * x1^2) * x2 +
                x1 * (
                    13 + 16 * x1 * sqrt(3 - 3 * x1^2) * x2 - 9 * x2^2 +
                    x1^2 * (-13 + 9 * x2^2)
                )
            )
        )
    ) / (8 * (q^2 + p * (p - q * x1 + q * sqrt(3 - 3 * x1^2) * x2)))

GGS2(q::Float64, p::Vector{Float64}, x1::Float64, x2::Float64) = @avx @. (
    p^3 * (
        p^3 * (-9 + 9 * x1^2 - 5 * x1 * sqrt(3 - 3 * x1^2) * x2) +
        4 *
        q^3 *
        (
            -(sqrt(3 - 3 * x1^2) * x2) + 2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 -
            x1 * (-1 + x1^2) * (-1 + 3 * x2^2)
        ) -
        2 *
        p *
        q^2 *
        (
            7 +
            3 * x2^2 +
            x1 * sqrt(3 - 3 * x1^2) * x2 * (1 - 3 * x2^2) +
            x1^4 * (1 + 3 * x2^2) +
            x1^3 * sqrt(3 - 3 * x1^2) * x2 * (1 + 3 * x2^2) -
            2 * x1^2 * (4 + 3 * x2^2)
        ) +
        p^2 *
        q *
        (
            -11 * sqrt(3 - 3 * x1^2) * x2 +
            x1 * (
                13 + 16 * x1 * sqrt(3 - 3 * x1^2) * x2 - 9 * x2^2 +
                x1^2 * (-13 + 9 * x2^2)
            )
        )
    )
) / (8 * (q^2 + p * (p - q * x1 + q * sqrt(3 - 3 * x1^2) * x2)))

@inline @fastmath GGS3(q::Float64, p::Float64, x1::Float64, x2::Float64) = (
    (
        p^3 * (
            8 *
            q^3 *
            (
                -(sqrt(3 - 3 * x1^2) * x2) +
                2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 -
                x1 * (-1 + x1^2) * (-1 + 3 * x2^2)
            ) +
            2 *
            p^2 *
            q *
            (
                sqrt(3 - 3 * x1^2) * x2 * (2 - 9 * x2^2) +
                4 * x1 * (1 + 3 * x2^2) - 4 * x1^3 * (1 + 3 * x2^2) +
                3 * x1^2 * sqrt(3 - 3 * x1^2) * x2 * (1 + 3 * x2^2)
            ) +
            4 *
            p *
            q^2 *
            (
                -5 + 12 * x2^2 + x1^2 * (1 - 15 * x2^2) -
                3 * x1 * sqrt(3 - 3 * x1^2) * x2 * (1 + x2^2) +
                x1^4 * (1 + 3 * x2^2) +
                x1^3 * sqrt(3 - 3 * x1^2) * x2 * (1 + 3 * x2^2)
            ) -
            p^3 * (
                6 + 4 * x1 * sqrt(3 - 3 * x1^2) * x2 - 15 * x2^2 +
                3 * x1^2 * (3 + 5 * x2^2)
            )
        )
    ) / (16 * (p^2 + q^2 + p * q * (x1 - sqrt(3 - 3 * x1^2) * x2)))
)

GGS3(q::Float64, p::Vector{Float64}, x1::Float64, x2::Float64) = @avx @.(
    (
        p^3 * (
            8 *
            q^3 *
            (
                -(sqrt(3 - 3 * x1^2) * x2) +
                2 * x1^2 * sqrt(3 - 3 * x1^2) * x2 -
                x1 * (-1 + x1^2) * (-1 + 3 * x2^2)
            ) +
            2 *
            p^2 *
            q *
            (
                sqrt(3 - 3 * x1^2) * x2 * (2 - 9 * x2^2) +
                4 * x1 * (1 + 3 * x2^2) - 4 * x1^3 * (1 + 3 * x2^2) +
                3 * x1^2 * sqrt(3 - 3 * x1^2) * x2 * (1 + 3 * x2^2)
            ) +
            4 *
            p *
            q^2 *
            (
                -5 + 12 * x2^2 + x1^2 * (1 - 15 * x2^2) -
                3 * x1 * sqrt(3 - 3 * x1^2) * x2 * (1 + x2^2) +
                x1^4 * (1 + 3 * x2^2) +
                x1^3 * sqrt(3 - 3 * x1^2) * x2 * (1 + 3 * x2^2)
            ) -
            p^3 * (
                6 + 4 * x1 * sqrt(3 - 3 * x1^2) * x2 - 15 * x2^2 +
                3 * x1^2 * (3 + 5 * x2^2)
            )
        )
    ) / (16 * (p^2 + q^2 + p * q * (x1 - sqrt(3 - 3 * x1^2) * x2)))
)

# Regulator and q-p def

# Rb(x) = x > 1e-6 ? x * (exp(1 * x^2) - 1)^-1 : 1 / x # when x is samll Rb~1/x
#
#
#
#
# Rbp(x) =
#     x > 1e-6 ? (-1 + exp(x^2))^-2 * (-1 + (1 + -2 * x^2) * exp(x^2)) : -1 / x^2
# # when x is samll Rbp~-1/x^2

@inline @fastmath Rb(x::Float64) = x > 1e-8 ? (-1 + x^-1) * (1 + exp((x - 1) * a^-1))^-1 : 1 / x

@inline @fastmath Rbp(x::Float64) = (
    x > 1e-8 ?
        -(a + (a + x - x^2) * exp((x - 1) * a^-1)) *
    (a * (1 + exp((x - 1) * a^-1))^2 * x^2)^-1 :
        -1 / x^2
)

function Rb(x::Vector{Float64})
    temp = similar(x)
    @avx for i = 1:N_pgrid
        temp[i] =
            x[i] > 1e-8 ? (-1 + x[i]^-1) * (1 + exp((x[i] - 1) * a^-1))^-1 :
            1 / x[i]
    end
    return temp
end

function Rbp(x::Vector{Float64})
    temp = similar(x)
    @avx for i = 1:N_pgrid
        temp[i] = x[i] > 1e-8 ?
            -(a + (a + x[i] - x[i]^2) * exp((x[i] - 1) * a^-1)) *
        (a * (1 + exp((x[i] - 1) * a^-1))^2 * x[i]^2)^-1 :
            -1 / x[i]^2
    end
    return temp
end

# Rbp(x::Float64)=x > 1e-6 ? -(a+(a+x-x^2)*exp((x-1) *a^-1))*(a*(1+exp((x-1) *a^-1))^2 *x^2)^-1 : -1 / x^2

@inline @fastmath Qmp(q::Float64, p::Float64, costh::Float64) = (p^2 - 2p * q * costh + q^2)

Qmp(q::Float64, p::Vector{Float64}, costh::Float64) =
    @avx @. (p^2 - 2p * q * costh + q^2)

@inline @fastmath sqrtQmp(q::Float64, p::Float64, costh::Float64) =
    sqrt(p^2 - 2p * q * costh + q^2)
sqrtQmp(q::Float64, p::Vector{Float64}, costh::Float64) =
    @avx @. sqrt(p^2 - 2p * q * costh + q^2)

max_stars_repo_name: Yangyang-Tan/FRG-Code-Collection

max_issues_repo_name: Yangyang-Tan/FRG-Code-Collection

max_forks_repo_name: Yangyang-Tan/FRG-Code-Collection


text:
using Base.Test
using RigidBodyDynamics
using RigidBodyDynamics: Bounds
using StaticArrays: SVector
using Cbc: CbcSolver
using LCPSim
using MeshCat
using MeshCatMechanisms

#=
"""
`care(A, B, Q, R)`

Compute 'X', the solution to the continuous-time algebraic Riccati equation,
defined as A'X + XA - (XB)R^-1(B'X) + Q = 0, where R is non-singular.

Algorithm taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf
"""
function care(A, B, Q, R)
    G = try
        B*inv(R)*B'
    catch
        error("R must be non-singular.")
    end

    Z = [A  -G;
        -Q  -A']

    S = schurfact(Z)
    S = ordschur(S, real(S.values).<0)
    U = S.Z

    (m, n) = size(U)
    U11 = U[1:div(m, 2), 1:div(n,2)]
    U21 = U[div(m,2)+1:m, 1:div(n,2)]
    return U21/U11
end

"""
`lqr(A, B, Q, R)`

Calculate the optimal gain matrix `K` for the state-feedback law `u = K*x` that
minimizes the cost function:

J = integral(x'Qx + u'Ru, 0, inf).

For the continuous time model `dx = Ax + Bu`.

`lqr(sys, Q, R)`

Solve the LQR problem for state-space system `sys`. Works for both discrete
and continuous time systems.

See also `LQG`

Usage example:
```julia
A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = eye(2)
R = eye(1)
L = lqr(sys,Q,R)

u(t,x) = -L*x # Form control law,
t=0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0)
plot(t,x, lab=["Position", "Velocity"]', xlabel="Time [s]")
```
"""
function lqr(A, B, Q, R)
    S = care(A, B, Q, R)
    K = R\B'*S
    return K
end

function dynamics(mechanism, x, u)
    state = MechanismState{eltype(x)}(mechanism)
    result = DynamicsResult{eltype(x)}(mechanism)
    ẋ = zeros(eltype(x), 4)
    dynamics!(ẋ, result, state, x, u)
    ẋ
end

A = ForwardDiff.jacobian(x -> dynamics(mechanism, x, zeros(eltype(x), 2)), [π, 0, 0, 0])
B = ForwardDiff.jacobian(u -> dynamics(mechanism, eltype(u)[π, 0, 0, 0], [0, u[2]]), [0., 0])
Q = diagm([10, 10, 1, 1])
R = diagm([1, 1])
K = lqr(A, B, Q, R)
=#

@testset "acrobot" begin
    K = [0.0 0.0 0.0 0.0; -266.087 -107.349 -114.95 -54.3175]
    urdf = joinpath(@__DIR__, "..", "examples", "urdf", "Acrobot.urdf")
    mechanism = parse_urdf(Float64, urdf)
    world = root_body(mechanism)
    env = Environment{Float64}([])

    x0 = MechanismState{Float64}(mechanism)
    set_velocity!(x0, zeros(num_velocities(x0)))
    set_configuration!(x0, findjoint(mechanism, "shoulder"), [π])
    set_configuration!(x0, findjoint(mechanism, "elbow"), [-0.1])
    q0 = copy(configuration(x0))
    v0 = copy(velocity(x0))

    controller! = let x0 = vcat(q0, v0), K = K
        (τ, t, x) -> begin
            x̂ = vcat(configuration(x), velocity(x)) - x0
            τ .= -K * x̂
        end
    end

    Δt = 0.05
    N = 100
    results = LCPSim.simulate(x0, controller!, env, Δt, N, CbcSolver());
    @test length(results) == N
    @test norm(velocity(results[end].state)) < 0.03
    @test norm(configuration(results[end].state) - q0) < 1.2

    vis = Visualizer()
    if !haskey(ENV, "CI")
        open(vis)
        wait(vis)
    end

    mvis = MechanismVisualizer(mechanism, URDFVisuals(urdf), vis)
    for r in results
        sleep(Δt)
        set_configuration!(mvis, configuration(r.state))
    end
end

ext: jl


text:
module SimpleI18n

export get_language, set_language, set_language!,
    i18n, @i_str

include("utils.jl")
include("locales.jl")
include("types.jl")
include("fallback.jl")
include("setup.jl")
include("main.jl")

end # module


text:
function Base.iterate(traj::Trajectory, state=1)
	state > length(traj.f) ? nothing : ((traj.t[state], traj.f[state]), state + 1)
end

Base.eltype(::Type{Trajectory{T}}) where T = Tuple{T, Float64}
Base.length(T::Trajectory) = length(T.f)

Base.in(f, fb::FrequencyBin) = fb.f - fb.df < f < fb.f + fb.df
Base.isless(x::FrequencyBin, y::FrequencyBin) = isless(x.f, y.f)

@recipe function f(traj::Trajectory)
	label --> ""
	return traj.t, traj.f
end
@recipe function f(traj::Trajectory, fb::FrequencyBin)
	if haskey(traj.timetofreq, fb)
		if traj.final_state == :lost
			linecolor --> :blue
		elseif traj.final_state == :fixed
			linecolor --> :red
		else
			linecolor --> :black
		end
		label --> ""
		t0 = traj.timetofreq[fb]
		return traj.t .- t0, traj.f
	else
		label --> ""
		return ()
	end
end


text:
mutable struct Maxout <: Block
    w::VarOrNil # input to middle hidden weights
    b::VarOrNil # bias of middle hidden units
    h::Int
    k::Int
    function Maxout(inputSize::Int, hiddenSize::Int; k::Int=2, type::Type=Array{Float32})
        @assert (k>=2) "# of Affine layers should no less than 2"
        T = eltype(type)
        d = hiddenSize * k
        w = randn(T, d, inputSize) .* sqrt( T(1/d) )
        b = zeros(T, d, 1)
        new(Variable{type}(w,true,true,true),
            Variable{type}(b,true,true,true), hiddenSize, k)
    end
    function Maxout(hiddenSize::Int; k::Int=2)
        @assert (k>=2) "# of Affine layers should no less than 2"
        new(nothing, nothing, hiddenSize, k)
    end
end

function clone(this::Maxout; type::Type=Array{Float32})
    cloned = Maxout(this.h, k=this.k)
    cloned.w = clone(this.w, type=type)
    cloned.b = clone(this.b, type=type)
    return cloned
end

# pretty show
function Base.show(io::IO, m::Maxout)
    SIZE = size(m.w)
    TYPE = typeof(m.w.value)
    maxk = m.k
    print(io, "Maxout($(SIZE[2]), $(SIZE[1]÷maxk); k=$maxk, type=$TYPE)")
end

"""
    unbiasedof(m::Maxout)

unbiased weights of `Maxout` block
"""
function unbiasedof(m::Maxout)
    weights = Vector(undef, 1)
    weights[1] = m.w.value
    return weights
end

function paramsof(m::Maxout)
    params = Vector{Variable}(undef,2)
    params[1] = m.w
    params[2] = m.b
    return params
end

function xparamsof(m::Maxout)
    xparams = Vector{XVariable}(undef,2)
    xparams[1] = ('w', m.w)
    xparams[2] = ('b', m.b)
    return xparams
end

function nparamsof(m::Maxout)
    lw = length(m.w)
    lb = length(m.b)
    return (lw + lb)
end

function bytesof(model::Maxout, unit::String="MB")
    n = nparamsof(model) * elsizeof(model.w)
    return blocksize(n, uppercase(unit))
end

function forward(model::Maxout, x::Variable{T}) where T
    h = model.h
    k = model.k
    w = model.w
    b = model.b
    c = size(x, 2)
    x = matAddVec(w * x, b)         # dim=(h*k, c)
    temp = reshape(x.value, h,k,c)  # dim=(h,k,c)
    maxv = maximum(temp, dims=2)    # dim=(h,1,c)
    mask = temp .== maxv            # dim=(h,k,c)
    out  = Variable{T}(reshape(maxv, h,c), x.backprop)
    if x.backprop
        function maxoutBackward()
            if need2computeδ!(x)
                x.delta += reshape(mask .* reshape(out.delta, h,1,c), h*k,c)
            end
            ifNotKeepδThenFreeδ!(out);
        end
        push!(graph.backward, maxoutBackward)
    end
    return out
end

function predict(model::Maxout, x::AbstractArray)
    h = model.h
    k = model.k
    w = model.w.value
    b = model.b.value
    c = size(x, 2)
    x = w * x .+ b                  # dim=(h*k, c)
    temp = reshape(x, h,k,c)        # dim=(h,k,c)
    maxv = maximum(temp, dims=2)    # dim=(h,1,c)
    out  = reshape(maxv, h,c)       # dim=(h,  c)
    return out
end

function to(type::Type, m::Maxout)
    m.w = to(type, m.w)
    m.b = to(type, m.b)
    return m
end

function to!(type::Type, m::Maxout)
    m.w = to(type, m.w)
    m.b = to(type, m.b)
    return nothing
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
####################################
# AbstractArray eltype replacement #
####################################

function eltype_param_number{T<:AbstractArray}(::Type{T})
    if T.name.name == :AbstractArray
        return 1
    else
        T_super = supertype(T)
        param_number = eltype_param_number(T_super)
        tv = T_super.parameters[param_number]
        for i = 1:T.parameters.length
            if tv == T.parameters[i]
                return i
            end
        end
    end
end

@generated function replace_eltype{T,S}(x::AbstractArray{T}, ::Type{S})
    if x.mutable
        pnum = eltype_param_number(x.name.primary)
        tname = x.name.primary
        tparams = collect(x.parameters)
        tparams[pnum] = S
        newtype = quote
            $(Expr(:meta, :inline))
            $(tname){$(tparams...)}
        end
    else
        newtype = quote
            $(Expr(:meta, :inline))
            typeof(similar(x, S, size(x)))
        end
    end
    return newtype
end

#######################################
# caching for Jacobians and gradients #
#######################################

const JACOBIAN_CACHE = Dict{Tuple{Int,Int,DataType,Bool},Any}()

immutable JacobianCache{N,T,D}
    duals::D
    seeds::NTuple{N,Partials{N,T}}
end

function JacobianCache{N}(x, chunk::Chunk{N})
    T = eltype(x)
    duals = similar(x, Dual{N,T}, size(x))
    seeds = construct_seeds(T, chunk)
    return JacobianCache{N,T,typeof(duals)}(duals, seeds)
end

@inline jacobian_dual_type{T,M,N}(arr::AbstractArray{T,M}, ::Chunk{N}) = replace_eltype(arr, Dual{N,T})

Base.copy(cache::JacobianCache) = JacobianCache(copy(cache.duals), cache.seeds)

@eval function multithread_jacobian_cachefetch!{X,N}(x::X, chunk::Chunk{N}, usecache::Bool,
                                                   alt::Bool = false)
    T, xlen = eltype(x), length(x)
    if usecache
        result = get!(JACOBIAN_CACHE, (xlen, N, X, alt)) do
            construct_jacobian_caches(x, chunk)
        end
    else
        result = construct_jacobian_caches(x, chunk)
    end
    return result::NTuple{$NTHREADS,JacobianCache{N,T,jacobian_dual_type(x, chunk)}}
end

jacobian_cachefetch!(args...) = multithread_jacobian_cachefetch!(args...)[compat_threadid()]

########################
# caching for Hessians #
########################

# only used for vector mode, so we can assume that N == length(x)
const HESSIAN_CACHE = Dict{Tuple{Int,DataType},Any}()

immutable HessianCache{N,T,D}
    duals::D
    inseeds::NTuple{N,Partials{N,T}}
    outseeds::NTuple{N,Partials{N,Dual{N,T}}}
end

function HessianCache{N}(x, chunk::Chunk{N})
    T = eltype(x)
    duals = similar(x, Dual{N,Dual{N,T}}, size(x))
    inseeds = construct_seeds(T, chunk)
    outseeds = construct_seeds(Dual{N,T}, chunk)
    return HessianCache{N,T,typeof(duals)}(duals, inseeds, outseeds)
end

@inline hessian_dual_type{T,M,N}(arr::AbstractArray{T,M}, ::Chunk{N}) = replace_eltype(arr, Dual{N,Dual{N,T}})

Base.copy(cache::HessianCache) = HessianCache(copy(cache.duals), cache.inseeds, cache.outseeds)

@eval function multithread_hessian_cachefetch!{X,N}(x::X, chunk::Chunk{N}, usecache::Bool)
    T = eltype(x)
    if usecache
        result = get!(HESSIAN_CACHE, (N, X)) do
            construct_hessian_caches(x, chunk)
        end
    else
        result = construct_hessian_caches(x, chunk)
    end
    return result::NTuple{$NTHREADS,HessianCache{N,T,hessian_dual_type(x, chunk)}}
end

hessian_cachefetch!(args...) = multithread_hessian_cachefetch!(args...)[compat_threadid()]

#################
# Partial seeds #
#################

function seedall!{N,T}(xdual, x, seed::Partials{N,T})
    for i in eachindex(xdual)
        xdual[i] = Dual{N,T}(x[i], seed)
    end
    return xdual
end

function seed!{N,T}(xdual, x, seed::Partials{N,T}, index)
    offset = index - 1
    for i in 1:N
        j = i + offset
        xdual[j] = Dual{N,T}(x[j], seed)
    end
    return xdual
end

function seed!{N,T}(xdual, x,seeds::NTuple{N,Partials{N,T}}, index, chunksize = N)
    offset = index - 1
    for i in 1:chunksize
        j = i + offset
        xdual[j] = Dual{N,T}(x[j], seeds[i])
    end
    return xdual
end

function seedhess!{N,T}(xdual, x, inseeds::NTuple{N,Partials{N,T}},
                        outseeds::NTuple{N,Partials{N,Dual{N,T}}})
    for i in 1:N
        xdual[i] = Dual{N,Dual{N,T}}(Dual{N,T}(x[i], inseeds[i]), outseeds[i])
    end
    return xdual
end

#####################
# @eval'd functions #
#####################

@eval function construct_jacobian_caches{N}(x, chunk::Chunk{N})
    result = JacobianCache(x, chunk)
    return $(Expr(:tuple, :result, [:(copy(result)) for i in 2:NTHREADS]...))
end

@eval function construct_hessian_caches{N}(x, chunk::Chunk{N})
    result = HessianCache(x, chunk)
    return $(Expr(:tuple, :result, [:(copy(result)) for i in 2:NTHREADS]...))
end

for N in 1:MAX_CHUNK_SIZE
    ex = Expr(:tuple, [:(setindex(zero_partials, seed_unit, $i)) for i in 1:N]...)
    @eval function construct_seeds{T}(::Type{T}, ::Chunk{$N})
        seed_unit = one(T)
        zero_partials = zero(Partials{$N,T})
        return $ex
    end
end


text:
# Copyright 2017 Stephen L. Smith and Frank Imeson
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

#########  Input Parsing ############################

"""
parse the problem instance file given as a text file with the following format:
N: int
M: int
Symmetric: true/false
Triangle: true/false
1 5
5 1 8 9 10 11
8 2 3 4 6 7 12 13 14
DISTANCE_MATRIX

Note: the line 5 1 8 9 10 11 shows that the second set contains 5 vertices:
1,8,9,10,11.

TSPLIB Parser defined by:
  http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/TSPFAQ.html
"""
function read_file(filename)
	  if !isfile(filename)
		    error("the problem instance  ", filename, " does not exist")
	  end

    # Setup
    INF = 9999
    RE_INT = r"-?\d+"
    RE_DEC = r"-?\d+\.\d+"
    RE_NUMBER = r"-?\d+\.?\d*"

    parse_state = "UNKNOWN_FORMAT"
    data_type = ""
    data_format = ""
    dist = zeros(Int64, 0, 0)
    sets = Any[]
    vid00 = vid01 = 1
    coords = Any[]
    set_data = Int64[]
    num_vertices = -1
    num_sets = -1

    # parse
    s = open(filename, "r")
    for line in readlines(s)
        line = strip(line)
        # debug
#        println(line)
#        println(parse_state)

        # auto format select
        if parse_state == "UNKNOWN_FORMAT"
            if occursin(r"^\s*NAME\s*:\s*\w+", uppercase(line))
                parse_state = "TSPLIB_HEADER"
            elseif occursin(r"^\s*N\s*:\s*\w+", uppercase(line))
                parse_state = "SIMPLE_HEADER"
            end
        end

        # Parse Setup
        if parse_state == "TSPLIB_HEADER"
            value = strip(split(line,":")[end])
            if occursin(r"^\s*NAME\s*:\s*\w+", uppercase(line))
            elseif occursin(r"^\s*TYPE\s*:\s*\w+\s*$", uppercase(line))
            elseif occursin(r"^\s*DIMENSION\s*:\s*\d+\s*$", uppercase(line))
                num_vertices = parse(Int64, value)
                dist = zeros(Int64, num_vertices, num_vertices)
            elseif occursin(r"^\s*GTSP_SETS\s*:\s*\d+\s*$", uppercase(line))
                num_sets = parse(Int64, value)
            elseif occursin(r"^\s*EDGE_WEIGHT_TYPE\s*:\s*\w+\s*$", uppercase(line))
                data_type = value
                data_format = value
            elseif occursin(r"^\s*EDGE_WEIGHT_FORMAT\s*:\s*\w+\s*$", uppercase(line))
                if data_type == "EXPLICIT"
                    data_format = value
                end
            elseif occursin(r"^\s*EDGE_WEIGHT_SECTION\s*:?\s*$", uppercase(line))
                parse_state = "TSPLIB_MATRIX_DATA"
            elseif occursin(r"^\s*NODE_COORD_SECTION\s*:?\s*$", uppercase(line))
                parse_state = "TSPLIB_COORD_DATA"
            end

        # Parse matrix data
        elseif parse_state == "TSPLIB_MATRIX_DATA"
            if occursin(r"^[\d\se+-\.]+$", line)
                for x in split(line)
                    cost = parse(Int64, x)
                    # tested
                    if data_format == "FULL_MATRIX"
                        dist[vid00, vid01] = cost
                        vid01 += 1
                        if vid01 > num_vertices
                            vid00 += 1
                            vid01  = 1
                        end
                    # tested
                    elseif data_format == "LOWER_DIAG_ROW"
                        dist[vid00, vid01] = cost
                        dist[vid01, vid00] = cost
                        vid01 += 1
                        if vid01 > vid00
                            vid00 += 1
                            vid01  = 1
                        end
                    # not tested
                    elseif data_format == "LOWER_ROW"
                        println("WARNING: Not tested")

                        if vid00 == 0 && vid01 == 0
                            vid00 = 1
                        end
                        dist[vid00, vid01] = cost
                        dist[vid01, vid00] = cost
                        vid01 += 1
                        if vid01 >= vid00
                            vid00 += 1
                            vid01  = 0
                        end
                    # tested
                    elseif data_format == "UPPER_DIAG_ROW"
                        dist[vid00, vid01] = cost
                        dist[vid01, vid00] = cost
                        vid01 += 1
                        if vid01 > num_vertices
                            vid00 += 1
                            vid01  = vid00
                        end
                    # tested
                    elseif data_format == "UPPER_ROW"
                        if vid00 == 1 && vid01 == 1
                            vid01 = 2
                        end
                        dist[vid00, vid01] = cost
                        dist[vid01, vid00] = cost
                        vid01 += 1
                        if vid01 > num_vertices
                            vid00 += 1
                            vid01  = vid00+1
                        end
                    end
                end
            elseif occursin(r"^\s*DISPLAY_DATA_SECTION\s*:?\s*$", uppercase(line))
                parse_state = "TSPLIB_DISPLAY_DATA"
            elseif occursin(r"^\s*GTSP_SET_SECTION\s*:?\s*$", uppercase(line))
                parse_state = "TSPLIB_SET_DATA"
            end

        # Parse display data
        elseif parse_state == "TSPLIB_DISPLAY_DATA"
            if occursin(r"^\s*GTSP_SET_SECTION\s*:?\s*$", uppercase(line))
                parse_state = "TSPLIB_SET_DATA"
            end

        # Parse coord data
        elseif parse_state == "TSPLIB_COORD_DATA"
            if occursin(r"\s*\d+\s*", uppercase(line))
                coord = [parse(Float64, x) for x = split(line)[2:end]]
                push!(coords, coord)
            elseif occursin(r"^\s*GTSP_SET_SECTION\s*:?\s*$", uppercase(line))
                parse_state = "TSPLIB_SET_DATA"
            end

        # Parse set data
        elseif parse_state == "TSPLIB_SET_DATA"
            if occursin(r"\d+", uppercase(line))
                for x = split(line)
                    push!(set_data, parse(Int64,x))
                end
            elseif occursin(r"^\s*EOF\s*$", uppercase(line))
                parse_state = "TSPLIB"
            end

        # Parse header (simple)
        elseif parse_state == "SIMPLE_HEADER"
            if occursin(r"^\s*N\s*:\s*\w+", uppercase(line))
                value = strip(split(strip(line),":")[end])
                num_vertices = parse(Int64, value)
                dist = zeros(Int64, num_vertices, num_vertices)
            elseif occursin(r"^\s*M\s*:\s*\d+\s*$", uppercase(line))
                value = strip(split(strip(line),":")[end])
                num_sets = parse(Int64, value)
                parse_state = "SIMPLE_SETS"
            end

        # Parse set data (simple)
        elseif parse_state == "SIMPLE_SETS"
            if occursin(r"^[\d\se+-\.]+$", line)
                sid = parse(Int64, split(line)[1])
                set = [parse(Int64, x) for x in split(line)[2:end]]
                push!(sets, set)
                if sid == num_sets
                    parse_state = "SIMPLE_MATRIX"
                end
            end

        # Parse set data (simple)
        elseif parse_state == "SIMPLE_MATRIX"
            if occursin(r"^[\d\se+-\.]+$", line)
                for x in split(line)
                    cost = parse(Int64, x)
                    dist[vid00, vid01] = cost
                    vid01 += 1
                    if vid01 > num_vertices
                        vid00 += 1
                        vid01  = 1
                    end
                end
            else
                parse_state = "SIMPLE"
            end
        end

    end
    close(s)
    if occursin(r"TSPLIB", parse_state)
        parse_state = "TSPLIB"
    end

    # Convert coordinate data to matrix data
    if parse_state == "TSPLIB" && data_type != "EXPLICIT"
        # tested
        if data_format == "EUC_2D"
            for vid00 in 1:num_vertices
                for vid01 in 1:num_vertices
                    if vid00 == vid01
                        dist[vid00, vid01] = INF
                    else
                        dx = coords[vid00][1] - coords[vid01][1]
                        dy = coords[vid00][2] - coords[vid01][2]
                        cost = sqrt(dx^2 + dy^2)
                        dist[vid00, vid01] = nint(cost)
                    end
                end
            end
        # not tested
        elseif data_format == "MAN_2D"
            println("Warning: MAN_2D not tested")
            for vid00 in 1:num_vertices
                for vid01 in 1:num_vertices
                    if vid00 == vid01
                        dist[vid00, vid01] = INF
                    else
                        dx = coords[vid00][1] - coords[vid01][1]
                        dy = coords[vid00][2] - coords[vid01][2]
                        cost = abs(dx) + abs(dy)
                        dist[vid00, vid01] = nint(cost)
                    end
                end
            end
        # not working...
        elseif data_format == "GEO"
            RRR = 6378.388
            PI = 3.141592
            DEBUG01 = false
            TSPLIB_GEO = true
            for vid00 in 1:num_vertices
                d, m    = degree_minutes(coords[vid00][1])
                lat00   = PI * (d + 5.0 * m / 3.0) / 180.0
                d, m    = degree_minutes(coords[vid00][2])
                long00  = PI * (d + 5.0 * m / 3.0) / 180.0
                for vid01 in 1:num_vertices
                    d, m    = degree_minutes(coords[vid01][1])
                    lat01   = PI * (d + 5.0 * m / 3.0) / 180.0
                    d, m    = degree_minutes(coords[vid01][2])
                    long01  = PI * (d + 5.0 * m / 3.0) / 180.0
                    if vid00 == vid01
                        dist[vid00, vid01] = INF
                    else
                        if TSPLIB_GEO
                            q1      = cos(long00 - long01)
                            q2      = cos(lat00 - lat01)
                            q3      = cos(lat00 + lat01)
                            cost    = RRR * acos(0.5 * ((1.0 + q1) * q2 - (1.0 - q1) * q3)) + 1.0
                            dist[vid00, vid01] = floor(Int64, cost)
                        else
                            # http://andrew.hedges.name/experiments/haversine/
                            R       = 6373
                            dlat    = abs(coords[vid01][1] - coords[vid00][1])
                            dlong   = abs(coords[vid01][2] - coords[vid00][2])
                            a       = (sind(dlat/2))^2 + cosd(coords[vid00][1])*cosd(coords[vid01][1])*(sind(dlong/2))^2
                            c       = 2*atan2(sqrt(a), sqrt(1-a))
                            cost    = R * c
                            dist[vid00, vid01] = floor(Int64, cost)
                        end
                        if DEBUG01
                            println("lat00 = ", coords[vid00][1], ", long00 = ", coords[vid00][2], ", lat01 = ", coords[vid01][1], ", long01 = ", coords[vid01][2], ", dist = ", floor(Int64, cost))
                            DEBUG01 = false
                        end
                    end
                end
            end
        # tested
        elseif data_format == "ATT"
            for vid00 in 1:num_vertices
                for vid01 in 1:num_vertices
                    if vid00 == vid01
                        dist[vid00, vid01] = INF
                    else
                        dx = coords[vid00][1] - coords[vid01][1]
                        dy = coords[vid00][2] - coords[vid01][2]
                        r  = sqrt((dx^2 + dy^2)/10.0)
                        cost = ceil(r)
                        dist[vid00, vid01] = nint(cost)
                    end
                end
            end
        # tested: not sure if working or not
        elseif data_format == "CEIL_2D"
            for vid00 in 1:num_vertices
                for vid01 in 1:num_vertices
                    if vid00 == vid01
                        dist[vid00, vid01] = INF
                    else
                        dx = coords[vid00][1] - coords[vid01][1]
                        dy = coords[vid00][2] - coords[vid01][2]
                        cost = sqrt(dx^2 + dy^2)
                        dist[vid00, vid01] = ceil(cost)
                    end
                end
            end
        else
            error("coordinate type $data_format not supported")
        end
    end

    # construct sets
    if parse_state == "TSPLIB"
        i = 1
        sid00 = 1
        set = Int64[]
        set_data = set_data[2:end]
        while i <= length(set_data)
            x = set_data[i]
            if x == -1
                push!(sets, set)
                set = Int64[]
                i += 1 # skip set id
            else
                push!(set, x)
            end
            i += 1
        end
        if num_sets != length(sets)
            error("number of sets doesn't match set size")
        end
    end

    if length(sets) <= 1
        error("must have more than 1 set")
    end

	  membership = findmember(num_vertices, sets)

    return num_vertices, num_sets, sets, dist, membership
end

""" Computing degrees and minutes for GEO instances """
function degree_minutes(num)
	if num > 0
		deg = floor(Int64, num)
		return deg, num - deg
	else
		deg = ceil(Int64, num)
		return deg, num - deg
	end
end

#####################################################
#####################################################
######### Output Printing ###########################

""" print the main parameter settings """
function print_params(param::Dict{Symbol,Any})
	if param[:print_output] > 0
		println("\n", "--------- Problem Data ------------")
		println("Instance Name      : ", param[:problem_instance])
	    println("Number of Vertices : ", param[:num_vertices])
	    println("Number of Sets     : ", param[:num_sets])
		println("Initial Tour       : ", (param[:init_tour] == "rand" ?
				"Random" : "Random Insertion"))
		println("Maximum Removals   : ", param[:max_removals])
		println("Trials             : ", param[:cold_trials])
		println("Restart Attempts   : ", param[:warm_trials])
		println("Rate of Adaptation : ", param[:epsilon])
		println("Prob of Reopt      : ", param[:prob_reopt])
		println("Maximum Time       : ", param[:max_time])
		println("Tour Budget        : ", (param[:budget] == typemin(Int64) ?
				"None" : param[:budget]))
		println("-----------------------------------\n")
	end
end

"""
Print the powers in an easy-to-read format for debugging the adaptive weights
"""
function print_powers(powers)
	println("--  Printing powers -- ")
	println("Insertions:")
	for power in powers["insertions"]
		println(power.name, " ", power.value, ": ", power.weight)
	end

	println("\n Removals:")
	for power in powers["removals"]
		println(power.name, " ", power.value, ": ", power.weight)
	end
	print("\n")

	println("\n Noises:")
	for power in powers["noise"]
		println(power.name, " ", power.value, ": ", power.weight)
	end
	print("\n")
end

"""print statement at the beginning of a cold trial"""
function print_cold_trial(count::Dict{Symbol,Real}, param::Dict{Symbol,Any}, best::Tour)
	if param[:print_output] == 2
		println("\n||--- trial ", count[:cold_trial],
		" --- initial cost ", best.cost, " ---||")
	end
end

"""print details at the end of each warm trial"""
function print_warm_trial(count::Dict{Symbol,Real}, param::Dict{Symbol,Any},
							best::Tour, iter_count::Int)
	if param[:print_output] == 2
		println("-- ", count[:cold_trial], ".", count[:warm_trial],
		" - iterations ", iter_count, ":  ", "cost ", best.cost)
	end
end

""" print best cost so far """
function print_best(count::Dict{Symbol,Real}, param::Dict{Symbol,Any},
							best::Tour, lowest::Tour, init_time::Float64)
	if param[:print_output] == 1 && time() - count[:print_time] > param[:print_time]
		count[:print_time] = time()
		println("-- trial ", count[:cold_trial], ".", count[:warm_trial], ":",
				"  Cost = ", min(best.cost, lowest.cost),
				"  Time = ", round(count[:print_time] - init_time, digits=1), " sec")

	elseif (param[:print_output] == 3 && time() - count[:print_time] > 0.5) ||
		param[:budget_met] || param[:timeout]

		count[:print_time] = time()
		if param[:warm_trials] > 0
			progress = (count[:cold_trial] - 1)/param[:cold_trials] +
	 	   			 (count[:warm_trial])/param[:warm_trials]/param[:cold_trials]
		else
			progress = (count[:cold_trial] - 1)/param[:cold_trials]
		end
		tcurr = round(count[:print_time] - init_time, digits=1)
		cost = min(best.cost, lowest.cost)
		progress_bar(param[:cold_trials], progress, cost, tcurr)
	end
end

""" a string representing the progress bar """
function progress_bar(trials, progress, cost, time)
	ticks, trials_per_bar, total_length = 6, 5, 31
	progress == 1.0 && (progress -= 0.0001)
	n = floor(Int64, progress * trials/trials_per_bar)
	start_number = n * trials_per_bar
	trials_in_bar = min(trials_per_bar, trials - start_number)

	progress_in_bar = (progress * trials - start_number)/trials_in_bar
	bar_length = min(total_length - 1, (trials - start_number) * ticks)

	progress_bar = "|"
	for i=1:total_length
		if i == bar_length + 1
			progress_bar *= "|"
		elseif i > bar_length + 1
			progress_bar *= " "
		elseif i % ticks == 1
			progress_bar *= string(start_number + ceil(Int64, i / ticks))
		elseif i <= ceil(Int64, bar_length * progress_in_bar)
			progress_bar *= "="
		else
			progress_bar *= " "
		end
	end
	print(" ", progress_bar, "  Cost = ", cost, "  Time = ", time, " sec      \r")
end

"""print tour summary at end of execution"""
function print_summary(lowest::Tour, timer::Float64, member::Array{Int64,1},
						param::Dict{Symbol,Any})
	if param[:print_output] == 3 && !param[:timeout] && !param[:budget_met]
		progress_bar(param[:cold_trials], 1.0, lowest.cost, round(timer, digits=1))
	end
	if param[:print_output] > -1
		if (param[:print_output] > 0 || param[:output_file] == "None")
			println("\n\n", "--------- Tour Summary ------------")
			println("Cost              : ", lowest.cost)
			println("Total Time        : ", round(timer, digits=2), " sec")
			println("Solver Timeout?   : ", param[:timeout])
			println("Tour is Feasible? : ", tour_feasibility(lowest.tour, member,
																	param[:num_sets]))
			order_to_print = (param[:output_file] == "None" ?
					lowest.tour : "printed to " * param[:output_file])
			println("Output File       : ",  param[:output_file])
			println("Tour Ordering     : ",  order_to_print)
			println("-----------------------------------")
		end
		if param[:output_file] != "None"
			s = open(param[:output_file], "w")
			write(s, "Problem Instance : ", param[:problem_instance], "\n")
			write(s, "Vertices         : ", string(param[:num_vertices]), "\n")
			write(s, "Sets             : ", string(param[:num_sets]), "\n")
			write(s, "Comment          : To avoid ~0.5sec startup time, use the Julia REPL\n")
			write(s, "Host Computer    : ", gethostname(), "\n")
			write(s, "Solver Time      : ", string(round(timer, digits=3)), " sec\n")
			write(s, "Tour Cost        : ", string(lowest.cost), "\n")
			write(s, "Tour             : ", string(lowest.tour))
			close(s)
		end
	end
end

"""
init function defined by TSPLIB
  http://comopt.ifi.uni-heidelberg.de/software/TSPLIB95/TSPFAQ.html
"""
function nint(x::Float64)
    return floor(Int64, x + 0.5)
end


text:
# This file includes mapping of interpolation interval to different subintervals.

using FractalTools
using Plots

# Construct interpolation data
xi, dx, xf = 0, 0.1, 1
x = collect(xi : dx : xf)
y = sin.(2π * x)
z = cos.(2π * x)
d = 0.9 * ones(length(x) - 1)
h = 0.9 * ones(length(x) - 1)
l = 0.9 * ones(length(x) - 1)
m = 0.9 * ones(length(x) - 1)

# Interpolate the data
@time hiddeninterp = hiddenfif(x, y, z, d, h, l, m)

# Get transformations
ws = hiddeninterp.ifs.ws

# Get an example transformation
pts = hcat.(eachcol([x'; y'; z']))

# Plots
plt = plot(layout=(4,1))
plot!(x, sin.(2π*x), subplot=1, label="")
vline!(x, subplot=1, label="", ls=:dot)
for (i, w) in enumerate(ws)
    v = getindex.(w.(pts), 1)
    plot!(v, ones(length(v)) * i, label="", subplot=2)
end
for (i, w) in enumerate(ws)
    v = getindex.(w.(pts), 2)
    plot!(ones(length(v)) * i, v, label="", subplot=3, lw=3)
end
for (i, w) in enumerate(ws)
    v = getindex.(w.(pts), 3)
    plot!(ones(length(v)) * i, v, label="", subplot=4)
end
display(plt)

# # Checks
# ifs = hiddeninterp.ifs

# # Analytical
# ei = x[1 : end - 1]
# ai = x[2 : end] - ei
# fi = y[1 : end - 1] - h
# ci = y[2 : end] - h - fi
# gi = z[1 : end - 1] - h
# ki = z[2 : end] - h - gi

# # Computed
# e = [w.b[1] for w in ws]
# a = [w.A[1, 1] for w in ws]
# f = [w.b[2] for w in ws]
# c = [w.A[2, 1] for w in ws]
# g = [w.b[3] for w in ws]
# k = [w.A[3, 1] for w in ws]

# # Compare
# e .== ei
# a .== ai
# f .== fi
# c .== ci
# g .== gi
# k .≈ ki


text:
abstract type AbstractEBayesMethod end

broadcastable(eb_method::AbstractEBayesMethod) = Ref(eb_method)


text:
module ElideableMacros

using StaticArrays

export @elidableassert, @elidableclamp, @elidableenv, @elidablenanzeroer
export @elidetimeandfilepath

"""
Show the environment variables
"""
macro elidableenv()
  quote
    filter(x -> occursin("ELIDE", first(x)), $(esc(ENV)))
  end
end

"""
The Julia @assert macro, except it can be compiled by one of two ways:
1) Set `ENV["ELIDE_ASSERTS"]` to one of `"yes"`, `"true"`, `"1"`, or `"on"`
or
2) Put the function `elideasserts() = true` in the module that calls @elidableassert
but not both
"""
macro elidableassert(assertion, messages...)
  local elide = haskey(ENV, "ELIDE_ASSERTS") &&
    ENV["ELIDE_ASSERTS"] ∈ ("yes", "true", "1", "on");
  elide |= @isdefined(elideasserts) ? elideasserts() : false;
  if @isdefined(elideasserts) && haskey(ENV, "ELIDE_ASSERTS") &&
    ENV["ELIDE_ASSERTS"] ∈ ("yes", "true", "1", "on");
    error("Set only ENV[\"ELIDE_ASSERTS\"] or elideasserts(), not both.");
  end
  local message = isempty(messages) ? assertion : messages[1]
  if !elide
    quote
      @assert $(esc(assertion)) $(esc(message))
    end
  end
end

"""
Replace NaNs with zeros
Compile out this macro by one of two ways:
1) Set `ENV["ELIDE_NANZEROER"]` to one of `"yes"`, `"true"`, `"1"`, or `"on"`
or
2) Put the function `elidenanzeroer() = false` in the module that calls @elidablenanzeroer
but not both
"""
macro elidablenanzeroer(value)
  local elide = haskey(ENV, "ELIDE_NANZEROER") &&
    ENV["ELIDE_NANZEROER"] ∈ ("yes", "true", "1", "on");
  elide |= @isdefined(elidenanzeroer) ? elidenanzeroer() : false;
  if @isdefined(elidenanzeroer) && haskey(ENV, "ELIDE_NANZEROER") &&
    ENV["ELIDE_NANZEROER"] ∈ ("no", "false", "0", "off");
    error("Set only ENV[\"ELIDE_NANZEROER\"] or elidenanzeroer(), not both.");
  end;
  if !elide
    return quote
      @inline _replace(x::Real) = ifelse(isnan(x), zero(x), x);
      @inline _replace(x::T) where {T<:Complex} = T(_replace(real(x)),
                                                    _replace(imag(x)))
      _replace.($(esc(value)))
    end
  else
    return quote $(esc(value)) end
  end
end

"""
Clamp values to within ± value or between lower and upper bounds
Compile out this macro by one of two ways:
1) Set `ENV["ELIDE_CLAMP"]` to one of `"yes"`, `"true"`, `"1"`, or `"on"`
or
2) Put the function `elideclamp() = false` in the module that calls @elideclamp
but not both
"""
macro elidableclamp(value, valueorlower, nothingorupper=nothing)
  isnothing(nothingorupper) && @assert valueorlower >= 0
  local lower = isnothing(nothingorupper) ? -valueorlower : valueorlower
  local upper = isnothing(nothingorupper) ? valueorlower : nothingorupper
  local elide = haskey(ENV, "ELIDE_CLAMP") &&
    ENV["ELIDE_CLAMP"] ∈ ("yes", "true", "1", "on");
  elide |= @isdefined(elideclamp) ? elideclamp() : false;
  if @isdefined(elideclamp) && haskey(ENV, "ELIDE_CLAMP") &&
    ENV["ELIDE_CLAMP"] ∈ ("no", "false", "0", "off");
    error("Set only ENV[\"ELIDE_CLAMP\"] or elideclamp(), not both.");
  end;
  if !elide
    return quote
      _replace(x::T) where {T<:Real} = clamp(x, $(lower), $(upper));
      _replace(x::T) where {T<:Complex} = T(_replace(real(x)),
                                            _replace(imag(x)))
      _replace.($(esc(value)))
    end
  else
    return quote $(esc(value)) end
  end
end

"""
Display the time, and the file path
Compile out this macro by one of two ways:
1) Set `ENV["ELIDE_TIMEANDFILEPATH"]` to one of `"yes"`, `"true"`, `"1"`, or `"on"`
or
2) Put the function `elidetimeandfilepath() = false` in the module that calls @elidetimeandfilepath
but not both
"""
macro timeandfilepath()
  local elide = haskey(ENV, "ELIDE_TIMEANDFILEPATH") &&
    ENV["ELIDE_TIMEANDFILEPATH"] ∈ ("yes", "true", "1", "on");
  elide |= @isdefined(elideclamp) ? elideclamp() : false;
  if @isdefined(elidetimeandfilepath) &&
      haskey(ENV, "ELIDE_TIMEANDFILEPATH") &&
      ENV["ELIDE_TIMEANDFILEPATH"] ∈ ("no", "false", "0", "off");
    error("Set only ENV[\"ELIDE_TIMEANDFILEPATH\"] or elidetimeandfilepath(), not both.");
  end;
  if !elide
    quote
      using Dates; println("$(now()) $(@__FILE__)")
    end
  end
end

end


text:
using OnlineStats, CSV, ProgressMeter, IterTools, DataFrames

function get_empirical_teams(filename::String; league = "great")
    data = CSV.read(filename, DataFrame)
    numEmpiricalTeams = nrow(data)
    data = hcat(silph_to_pvpoke.(data[:, 1:3]), data[:, 4])
    rankings = get_rankings(league)
    empiricalTeams = Array{StaticTeam}(undef, numEmpiricalTeams)
    for i = 1:numEmpiricalTeams
        if get_rankings_mon_id(data[i, 1], league = league) != 0 &&
           get_rankings_mon_id(data[i, 2], league = league) != 0 &&
           get_rankings_mon_id(data[i, 3], league = league) != 0
            empiricalTeams[i] = StaticTeam(
                [data[i, 1] data[i, 2] data[i, 3]],
                league = league,
            )
        else
            @warn "Ignoring team $(data[i, 1]), $(data[i, 2]), $(data[i, 3])"
        end
    end
    weights = data[:, 4]
    i = 1
    while true
        if i > length(empiricalTeams)
            break
        end
        if !isassigned(empiricalTeams, i)
            deleteat!(empiricalTeams, i)
            deleteat!(weights, i)
            i -= 1
        end
        i += 1
    end
    return empiricalTeams, weights
end;

function get_theoretical_teams(numMons::Int64; league = "great")
    theoreticalMons = Array{String}(undef, numMons)
    rankings = get_rankings(league)
    for i = 1:numMons
        theoreticalMons[i] = rankings[i]["speciesId"]
    end
    dexDict = Dict{Int16,Array{Int64,1}}()
    for i = 1:numMons
        @inbounds dex = gamemaster["pokemon"][
            get_gamemaster_mon_id(theoreticalMons[i])]["dex"]
        if haskey(dexDict, dex)
            push!(dexDict[dex], i)
        else
            dexDict[dex] = [i]
        end
    end
    dexKeys = collect(keys(dexDict))
    index = 0
    for i = 1:length(dexKeys)
        for j in Iterators.flatten((1:(i-1), (i+1):length(dexKeys)))
            for k = (j+1):length(dexKeys)
                if i != k
                    @inbounds toAdd = length(dexDict[dexKeys[i]]) *
                                        length(dexDict[dexKeys[j]]) *
                                        length(dexDict[dexKeys[k]])
                    index += toAdd
                end
            end
        end
    end
    numTheoreticalTeams = index
    theoreticalTeams = Array{StaticTeam}(undef, numTheoreticalTeams)
    index = 1
    dexes(x) = dexDict[dexKeys[x]]
    for i = 1:length(dexKeys)
        for j in Iterators.flatten((1:(i-1), (i+1):length(dexKeys)))
            for k = (j+1):length(dexKeys)
                if i != k
                    @inbounds toAdd = length(dexDict[dexKeys[i]]) *
                                        length(dexDict[dexKeys[j]]) *
                                        length(dexDict[dexKeys[k]])
                    @inbounds theoreticalTeams[index:(index+toAdd-1)] = [
                        StaticTeam(
                            [l, m, n], league = league
                        ) for l in dexes(i), m in dexes(j), n in dexes(k)]
                    index += toAdd
                end
            end
        end
    end
    return theoreticalTeams, theoreticalMons
end;

function run_empirical_teams(
    theoreticalTeam::StaticTeam,
    empiricalTeams::Array{StaticTeam},
    weights::Array{Int64},
)
    histogram = Hist(0.0:0.025:1.0)
    @simd for i = 1:length(empiricalTeams)
        @simd for j = 1:weights[i]
            static_state = StaticState(theoreticalTeam, empiricalTeams[i])
            dynamic_state = DynamicState(static_state)
            fit!(
                histogram,
                play_battle(dynamic_state, static_state),
            )
        end
    end
    return histogram
end;

function run_theoretical_teams(
    theoreticalTeams::Array{StaticTeam},
    empiricalTeams::Array{StaticTeam},
    weights,
)
    histograms = Array{Hist}(undef, length(theoreticalTeams))
    @showprogress for i = 1:length(theoreticalTeams)
        @inbounds histograms[i] = run_empirical_teams(
            theoreticalTeams[i],
            empiricalTeams,
            weights,
        )
    end
    expected_wins = get_expected_win.(histograms, sum(weights))
    expected_battle_score = mean.(histograms)
    return (histograms, expected_wins, expected_battle_score)
end;

function get_expected_win(histogram::Hist, numEmpiricalTeams::Int64)
    sum(histogram.counts[21:40]) / (numEmpiricalTeams / 5)
end;

function get_mon_name(mon::StaticPokemon, mon_list::Array{String})
    i = mon_list[findfirst(x -> StaticPokemon(x) == mon, mon_list)]
end

function get_summary_stats(
    histograms,
    expected_wins,
    expected_battle_score,
    theoreticalTeams,
    theoreticalMons,
)
    summaryStats = Array{Any}(undef, length(theoreticalTeams), 5)

    for i = 1:length(theoreticalTeams)
        @inbounds summaryStats[i, :] = [expected_wins[i]
            expected_battle_score[i] get_mon_name(theoreticalTeams[i].mons[1],
            theoreticalMons) get_mon_name(theoreticalTeams[i].mons[2],
            theoreticalMons) get_mon_name(theoreticalTeams[i].mons[3],
            theoreticalMons)]
    end
    return summaryStats
end;

function rank(numMons, indigo_file, outfile; league = "great")
    println("Constructing Empirical Teams...")
    empiricalTeams, weights = get_empirical_teams(indigo_file, league = league)
    println("Constructing Theoretical Teams...")
    theoreticalTeams, theoreticalMons = get_theoretical_teams(numMons,
        league = league)
    println("Running Battles...")
    histograms, expected_wins, expected_battle_score = run_theoretical_teams(
        theoreticalTeams,
        empiricalTeams,
        weights,
    )
    println("Sorting and Saving...")
    summaryStats = get_summary_stats(
        histograms,
        expected_wins,
        expected_battle_score,
        theoreticalTeams,
        theoreticalMons,
    )
    summaryStats = sortslices(
        summaryStats,
        by = x -> x[1],
        dims = 1,
        rev = true,
    )
    df = DataFrame(summaryStats)
    CSV.write(outfile, df)
    println("Done.")
    return df
end;


text:
function get_inverse(obj::FieldDecompositionSolver)
    return jcall(obj, "getInverse", FieldMatrix, ())
end

function is_non_singular(obj::FieldDecompositionSolver)
    return jcall(obj, "isNonSingular", jboolean, ())
end

function solve(obj::FieldDecompositionSolver, arg0::FieldMatrix)
    return jcall(obj, "solve", FieldMatrix, (FieldMatrix,), arg0)
end

function solve(obj::FieldDecompositionSolver, arg0::FieldVector)
    return jcall(obj, "solve", FieldVector, (FieldVector,), arg0)
end


text:
include(joinpath(module_dir(), "src", "dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# qet hopper simulation
s = get_simulation("hopper_2D", "flat_2D_lc", "flat")
model = s.model
env = s.env
nq = s.model.dim.q
nu = model.dim.u
nc = model.dim.c
nb = nc * friction_dim(env)
nz = num_var(model, env)
nθ = num_data(model)

H = 92
# h = 0.04
# h = 0.02
h = 0.01
κ = 1.0e-8

# Design open-loop control trajectory
q0_ref = SVector{nq,T}([0.0, 0.50, 0.0, 0.5])
q1_ref = SVector{nq,T}([0.0, 0.50, 0.0, 0.5])

α = 0.0077 # N=46 h=0.01

u_ref = []
push!(u_ref,  [[0.0,  5.0*α*model.g*(model.mb+model.ml)/2] for k=1:2*6]...);
push!(u_ref,  [[0.0, -0.60*α*model.g*(model.mb+model.ml)/2] for k=1:2*10]...);
push!(u_ref,  [[0.0,  0.14*α*model.g*(model.mb+model.ml)/2] for k=1:2*15]...);
push!(u_ref,  [[0.0,  2.19*α*model.g*(model.mb+model.ml)/2] for k=1:H-length(u_ref)]...);

# Simulate
sim = simulator(s, q0_ref, q1_ref, h, H;
    p = open_loop_policy(u_ref),
    d = no_disturbances(s.model),
    ip_opts = InteriorPointOptions(r_tol=1e-9, κ_tol=2κ),
    sim_opts = SimulatorOptions())

simulate!(sim)

plot(hcat(Vector.(sim.traj.u)...)')
plot(hcat(Vector.(sim.traj.q)...)')
visualize_robot!(vis, s.model, sim.traj.q)

for t = 1:10
    @show norm(residual(s.model, s.env, sim.traj.z[t], sim.traj.θ[t], [κ]))
end

# Check ~perfect loop
@show round.(sim.traj.q[end-2] - q0_ref, digits=3)
@show round.(sim.traj.q[end-1] - q1_ref, digits=3)

# Save trajectory
traj = deepcopy(sim.traj)
gait_path = joinpath(@__DIR__, "..", "dynamics", "hopper_2D", "gaits", "gait_in_place.jld2")
@save gait_path traj

# Reload trajectory
res = JLD2.jldopen(gait_path)
loaded_traj = res["traj"]

traj = get_trajectory(s.model, s.env,
    joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/gait_in_place.jld2"),
    load_type = :joint_traj)
plot(hcat(Vector.(traj.q)...)')


text:
function graphml_read_one_graph!(el::EzXML.Node, ::Type{G}) where G
    nodes = Dict{String,Int}()
    E = Edge{Int}
    edges = Vector{E}()

    nodeid = 1
    for f in eachelement(el)
        if nodename(f) == "node"
            nodes[f["id"]] = nodeid
            nodeid += 1
        elseif nodename(f) == "edge"
            n1 = f["source"]
            n2 = f["target"]
            push!(edges, E(nodes[n1], nodes[n2]))
        else
            warn("Skipping unknown node '$(nodename(f))'")
        end
    end
    g = G(length(nodes))
    for edge in edges
        add_edge!(g, edge)
    end
    return g
end

graphmlparse(T, x::String) = parse(T, x)
graphmlparse(::Type{String}, x::String) = x
graphmlparse(::Type{Vector{T}}, x::String) where {T} = parse.(T, split(x, ','))

function graphml_read_one_net!(xg::EzXML.Node, ::Type{G},
                        gpropkeys, vpropkeys, epropkeys) where G
    nodes = Dict{String,Int}()
    nodeid = 1
    # traverse the tree to map id to 1:n
    for f in eachelement(xg)
        nodename(f) != "node" && continue
        nodes[f["id"]] = nodeid
        nodeid += 1
    end

    g = G(length(nodes))
    for (pname,T) in values(vpropkeys); vprop!(g, pname, T); end
    for (pname,T) in values(epropkeys); eprop!(g, pname, T); end

    for f in eachelement(xg)
        if nodename(f) == "node"
            for el in eachelement(f)
                nodename(el) != "data" && continue
                i = nodes[f["id"]]
                pname, T = vpropkeys[el["key"]]
                m = vprop(g, pname)
                m[i] = graphmlparse(T, nodecontent(el))
            end
        elseif nodename(f) == "edge"
            n1 = f["source"]
            n2 = f["target"]
            ok, e = add_edge!(g, nodes[n1], nodes[n2]) #TODO
            for el in eachelement(f)
                nodename(el) != "data" && continue
                pname, T = epropkeys[el["key"]]
                m = eprop(g, pname)
                m[e] = graphmlparse(T, nodecontent(el))
            end
        elseif nodename(f) == "data"
            pname, T = gpropkeys[f["key"]]
            gprop!(g, pname, graphmlparse(T, nodecontent(f)))
        else
            warn("Skipping unknown xml-node '$(nodename(f))'")
        end
    end
    return g
end

function readgraphml(io::IO, ::Type{G}) where G<:AGraphOrDiGraph
    xdoc = parsexml(read(io, String))
    xroot = root(xdoc)  # an instance of XMLElement
    nodename(xroot) == "graphml" || error("Not a GraphML file")
    xg = getchild(xroot, "graph")
    isdir = false
    if haskey(xg, "edgedefault")
        isdir = xg["edgedefault"] == "directed"  ? true  : false
    end
    H = isdir ? digraphtype(G) : graphtype(G)
    return graphml_read_one_graph!(xg, H)
end

function readnetgraphml(io::IO, ::Type{G}) where G<:AGraphOrDiGraph
    xdoc = parsexml(read(io, String))
    xroot = root(xdoc)  # an instance of XMLElement
    nodename(xroot) == "graphml" || error("Not a GraphML file")
    xg = getchild(xroot, "graph")
    isdir = false
    if haskey(xg, "edgedefault")
        isdir = xg["edgedefault"] == "directed"  ? true  : false
    end

    gpropkeys=Dict{String, Tuple{String, DataType}}();
    vpropkeys=Dict{String, Tuple{String, DataType}}();
    epropkeys=Dict{String, Tuple{String, DataType}}();
    for el in eachelement(xroot)
        nodename(el) != "key" && continue
        if el["for"] == "graph"
            gpropkeys[el["id"]] = (el["attr.name"], graphml_types_rev[el["attr.type"]])
        elseif el["for"] == "node"
            vpropkeys[el["id"]] = (el["attr.name"], graphml_types_rev[el["attr.type"]])
        elseif el["for"] == "edge"
            epropkeys[el["id"]] = (el["attr.name"], graphml_types_rev[el["attr.type"]])
        end
    end

    H = isdir ? digraphtype(G) : graphtype(G)
    return graphml_read_one_net!(xg, H, gpropkeys, vpropkeys, epropkeys)
end

function writegraphml(io::IO, g::AGraphOrDiGraph)
    xdoc = XMLDocument()
    xroot = setroot!(xdoc, ElementNode("graphml"))
    xroot["xmlns"] = "http://graphml.graphdrawing.org/xmlns"
    xroot["xmlns:xsi"] = "http://www.w3.org/2001/XMLSchema-instance"
    xroot["xsi:schemaLocation"] = "http://graphml.graphdrawing.org/xmlns/1.0/graphml.xsd"

    xg = addelement!(xroot, "graph")
    xg["edgedefault"] = is_directed(g) ? "directed" : "undirected"
    for i in 1:nv(g)
        xv = addelement!(xg, "node")
        xv["id"] = "n$(i-1)"
    end
    for e in edges(g)
        xe = addelement!(xg, "edge")
        xe["source"] = "n$(src(e)-1)"
        xe["target"] = "n$(dst(e)-1)"
    end
    prettyprint(io, xdoc)
    return 1
end

const graphml_types = Dict( Int32 => "int",
                            Int64 => "int",
                            Bool => "boolean",
                            Float32 => "double",
                            Float64 => "double",
                            String => "string",
                            Vector{Float32} => "vector_double",
                            Vector{Float64} => "vector_double"
                        )

const graphml_types_rev = Dict("int"  =>  Int,
                            "long"    =>  Int,
                            "boolean" =>  Bool,
                            "float"   =>  Float32,
                            "double"  =>  Float64,
                            "string"  =>  String,
                            "vector_float" => Vector{Float64},
                            "vector_double" => Vector{Float64}
                        )

graphmlstring(x) = string(x)
graphmlstring(v::Vector) = join((Printf.@sprintf("%.10g",x) for x in v), ", ")

function writenetgraphml(io::IO, g::ANetOrDiNet)
    xdoc = XMLDocument()
    xroot = setroot!(xdoc, ElementNode("graphml"))
    xroot["xmlns"] = "http://graphml.graphdrawing.org/xmlns"
    xroot["xmlns:xsi"] = "http://www.w3.org/2001/XMLSchema-instance"
    xroot["xsi:schemaLocation"] = "http://graphml.graphdrawing.org/xmlns/1.0/graphml.xsd"

    nprop = 0
    gpropkey = Dict{String,String}()
    for (pname, p) in gprop(g)
        xp = addelement!(xroot, "key")
        xp["id"] = "key$nprop"
        gpropkey[pname] = "key$nprop"
        xp["for"] = "graph"
        xp["attr.name"] = pname
        xp["attr.type"] = graphml_types[typeof(p)]
        nprop+=1
    end
    vpropkey = Dict{String,String}()
    for (pname, p) in vprop(g)
        xp = addelement!(xroot, "key")
        xp["id"] = "key$nprop"
        vpropkey[pname] = "key$nprop"
        xp["for"] = "node"
        xp["attr.name"] = pname
        xp["attr.type"] = graphml_types[valtype(p)]
        nprop+=1
    end
    epropkey = Dict{String,String}()
    for (pname, p) in eprop(g)
        xp = addelement!(xroot, "key")
        xp["id"] = "key$nprop"
        epropkey[pname] = "key$nprop"
        xp["for"] = "edge"
        xp["attr.name"] = pname
        xp["attr.type"] = graphml_types[valtype(p)]
        nprop+=1
    end

    xg = addelement!(xroot, "graph")
    xg["edgedefault"] = is_directed(g) ? "directed" : "undirected"
    for (pname, p) in gprop(g)
        xp = addelement!(xg, "data")
        xp["key"] = gpropkey[pname]
        setnodecontent!(xp, graphmlstring(p))
    end

    for i in 1:nv(g)
        xv = addelement!(xg, "node")
        xv["id"] = "n$(i-1)"
        for (pname,p) in vprop(g)
            xp = addelement!(xv, "data")
            xp["key"] = vpropkey[pname]
            setnodecontent!(xp, graphmlstring(p[i])) #TODO check has key
        end
    end

    for e in edges(g)
        xe = addelement!(xg, "edge")
        xe["source"] = "n$(src(e)-1)"
        xe["target"] = "n$(dst(e)-1)"
        for (pname, p) in eprop(g)
            xp = addelement!(xe, "data")
            xp["key"] = epropkey[pname]
            setnodecontent!(xp, graphmlstring(p[e])) #TODO check has key
        end
    end
    prettyprint(io, xdoc)
    return 1
end

filemap[:graphml] = (readgraphml, writegraphml, readnetgraphml, writenetgraphml)


text:
#= #TODO: find the following Items
abstract GtkWidgetI <: GObjectI
abstract GtkContainerI <: GtkWidgetI
abstract GtkBinI <: GtkContainerI
#@@ -9,65 +8,7 @@ abstract GtkBoxI <: GtkContainerI
abstract GtkMenuShellI <: GtkContainerI
abstract GtkMenuItemI <: GtkMenuShellI
=#

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
@not_implemented(::BitbucketAPI, ::typeof(get_branch), ::String, ::String, ::String)
@not_implemented(::BitbucketAPI, ::typeof(get_branches), ::String, ::String)
@not_implemented(::BitbucketAPI, ::typeof(delete_branch), ::String, ::String, ::String)


text:
function getAllAdaptersInBag(input)
    allAdapters::Vector{Int} = []
    for line in readlines(input)
        push!(allAdapters, parse(Int,line))
    end
    return sort(allAdapters)
end

function calculateAllOutlets(allAdapters)
    plane::Vector{Int} = [0]
    device::Vector{Int} = [findmax(allAdapters)[1] + 3]
    return vcat(plane,allAdapters,device)
end

function solveProblemA()
    allAdapters = getAllAdaptersInBag("data.txt")
    allOutlets = calculateAllOutlets(allAdapters)
    diffs::Vector{Int} = [0,0,0]
    for i in 2:length(allOutlets)
        diff = allOutlets[i] - allOutlets[i-1]
        diffs[diff] +=1
    end
    println(diffs[1]*diffs[3])
end

solveProblemA()


text:
"""
`trymkpath(path; kwargs...)` make path just like `mkpath(path; kwargs...)`, but does nothing rather than throwing an error if the directory already exists.
"""
function trymkpath(path; kwargs...)
    try
        mkpath(path; kwargs...)
    catch
        @info """
            "$path" already exist; do nothing.
        """
    end
end


text:
for p in ("Knet","ArgParse")
    Pkg.installed(p) == nothing && Pkg.add(p)
end
using Knet
!isdefined(:MNIST) && (local lo=isdefined(:load_only); load_only=true; include(Knet.dir("examples","mnist.jl")); load_only=lo)

"""

This example learns to classify hand-written digits from the
[MNIST](http://yann.lecun.com/exdb/mnist) dataset.  There are 60000
training and 10000 test examples. Each input x consists of 784 pixels
representing a 28x28 image.  The pixel values are normalized to
[0,1]. Each output y is converted to a ten-dimensional one-hot vector
(a vector that has a single non-zero component) indicating the correct
class (0-9) for a given image.  10 is used to represent 0.

You can run the demo using `julia lenet.jl`.  Use `julia lenet.jl
--help` for a list of options.  The dataset will be automatically
downloaded.  By default the [LeNet](http://yann.lecun.com/exdb/lenet)
convolutional neural network model will be trained for 10 epochs.  The
accuracy for the training and test sets will be printed at every epoch
and optimized parameters will be returned.

"""
module LeNet
using Knet,ArgParse
using Main.MNIST: minibatch, accuracy, xtrn, ytrn, xtst, ytst

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="lenet.jl (c) Deniz Yuret, 2016. The LeNet model on the MNIST handwritten digit recognition problem from http://yann.lecun.com/exdb/mnist."
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
        ("--batchsize"; arg_type=Int; default=100; help="minibatch size")
        ("--lr"; arg_type=Float64; default=0.1; help="learning rate")
        ("--fast"; action=:store_true; help="skip loss printing for faster run")
        ("--epochs"; arg_type=Int; default=3; help="number of epochs for training")
        ("--gcheck"; arg_type=Int; default=0; help="check N random gradients per parameter")
    end
    println(s.description)
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    println("opts=",[(k,v) for (k,v) in o]...)
    o[:seed] > 0 && srand(o[:seed])
    gpu() >= 0 || error("LeNet only works on GPU machines.")

    dtrn = minibatch4(xtrn, ytrn, o[:batchsize])
    dtst = minibatch4(xtst, ytst, o[:batchsize])
    w = weights()
    report(epoch)=println((:epoch,epoch,:trn,accuracy(w,dtrn,predict),:tst,accuracy(w,dtst,predict)))

    if o[:fast]
        @time (train(w, dtrn; lr=o[:lr], epochs=o[:epochs]); gpu()>=0 && Knet.cudaDeviceSynchronize())
    else
        report(0)
        @time for epoch=1:o[:epochs]
            train(w, dtrn; lr=o[:lr], epochs=1)
            report(epoch)
            if o[:gcheck] > 0
                gradcheck(loss, w, first(dtrn)...; gcheck=o[:gcheck])
            end
        end
    end
    return w
end

function train(w, data; lr=.1, epochs=20, nxy=0)
    for epoch=1:epochs
        for (x,y) in data
            g = lossgradient(w, x, y)
            for i in 1:length(w)
                # w[i] -= lr * g[i]
                axpy!(-lr, g[i], w[i])
            end
        end
    end
    return w
end

function predict(w,x,n=length(w)-4)
    for i=1:2:n
        x = pool(relu(conv4(w[i],x;padding=0) .+ w[i+1]))
    end
    x = mat(x)
    for i=n+1:2:length(w)-2
        x = relu(w[i]*x .+ w[i+1])
    end
    return w[end-1]*x .+ w[end]
end

function loss(w,x,ygold)
    ypred = predict(w,x)
    ynorm = logp(ypred,1)  # ypred .- log(sum(exp(ypred),1))
    -sum(ygold .* ynorm) / size(ygold,2)
end

lossgradient = grad(loss)

function weights(;ftype=Float32,atype=KnetArray)
    w = Array(Any,8)
    w[1] = xavier(Float32,5,5,1,20)
    w[2] = zeros(Float32,1,1,20,1)
    w[3] = xavier(Float32,5,5,20,50)
    w[4] = zeros(Float32,1,1,50,1)
    w[5] = xavier(Float32,500,800)
    w[6] = zeros(Float32,500,1)
    w[7] = xavier(Float32,10,500)
    w[8] = zeros(Float32,10,1)
    return map(a->convert(atype,a), w)
end

function minibatch4(x, y, batchsize; atype=KnetArray{Float32})
    data = minibatch(x,y,batchsize; atype=atype)
    for i=1:length(data)
        (x,y) = data[i]
        data[i] = (reshape(x, (28,28,1,batchsize)), y)
    end
    return data
end

function xavier(a...)
    w = rand(a...)
     # The old implementation was not right for fully connected layers:
     # (fanin = length(y) / (size(y)[end]); scale = sqrt(3 / fanin); axpb!(rand!(y); a=2*scale, b=-scale)) :
    if ndims(w) < 2
        error("ndims=$(ndims(w)) in xavier")
    elseif ndims(w) == 2
        fanout = size(w,1)
        fanin = size(w,2)
    else
        fanout = size(w, ndims(w)) # Caffe disagrees: http://caffe.berkeleyvision.org/doxygen/classcaffe_1_1XavierFiller.html#details
        fanin = div(length(w), fanout)
    end
    # See: http://jmlr.org/proceedings/papers/v9/glorot10a/glorot10a.pdf
    s = sqrt(2 / (fanin + fanout))
    w = 2s*w-s
end

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia lenet.jl --epochs 10
# julia> LeNet.main("--epochs 10")
!isinteractive() && (!isdefined(Main,:load_only) || !Main.load_only) && main(ARGS)

end # module

# SAMPLE RUN 65f57ff+ Wed Sep 14 10:02:30 EEST 2016
#
# lenet.jl (c) Deniz Yuret, 2016. The LeNet model on the MNIST handwritten digit recognition problem from http://yann.lecun.com/exdb/mnist.
# opts=(:seed,-1)(:batchsize,100)(:epochs,3)(:lr,0.1)(:gcheck,0)(:fast,true)
# ..................
# 9.319163 seconds (5.84 M allocations: 277.927 MB, 7.37% gc time)


text:
using Genie, App, App.Todos

### Your tests here
@test 1 == 1

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
#=

  Run the Project Euler programs (# 1:50) in Julia.

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using Distributed
include("Euler.jl")
total = 0
@time @sync @distributed for p in 1:50
    global total
    s = "euler$p.jl"
    if isfile(s)
        # include(s)
        t = @timed include(s)
        total += t.time
    end
end
println("Run total: $(total)s")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module TreeSitterHighlight

include("./libtree_sitter_highlight.jl")

import .LibTreeSitterHighlight as LTSH

function maybe_throw_ts_error(maybe_err)
    err(name) = error("Got error of type $name")

    if maybe_err == LibTreeSitterHighlight.TSHighlightOk
        # pass
    elseif maybe_err == LibTreeSitterHighlight.TSHighlightUnknownScope
        err("TSHighlightError::TSHighlightUnknownScope")
    elseif maybe_err == LibTreeSitterHighlight.TSHighlightTimeout
        err("TSHighlightError::TSHighlightTimeout")
    elseif maybe_err == LibTreeSitterHighlight.TSHighlightInvalidLanguage
        err("TSHighlightError::TSHighlightInvalidLanguage")
    elseif maybe_err == LibTreeSitterHighlight.TSHighlightInvalidUtf8
        err("TSHighlightError::TSHighlightInvalidUtf8")
    elseif maybe_err == LibTreeSitterHighlight.TSHighlightInvalidRegex
        err("TSHighlightError::TSHighlightInvalidRegex")
    elseif maybe_err == LibTreeSitterHighlight.TSHighlightInvalidQuery
        err("TSHighlightError::TSHighlightInvalidQuery")
    end
end

macro tscall(call)
    quote
        local res = $(esc(call))
        maybe_throw_ts_error(res)
    end
end

struct Language
    name::Symbol
    language::Ptr{Nothing}
end

mutable struct Highlighter
    # Keep these strings around
    names::Vector{String}
    attributes::Vector{String}

    ptr::Ptr{LTSH.TSHighlighter}
    languages::Vector{Language}

    function Highlighter(names, attributes, ptr)
        finalizer(LTSH.ts_highlighter_delete, new(names, attributes, ptr, Language[]))
    end
end

"""
    HighlightBuffer(names_attributes::Dict{String,String})

An HighlightBuffer is the object needed to perform syntax highlighting. The names
parameter represent the capture groups that should be wrapped in spans. The attributes
element corresponds to the corresponding HTML attributes that should be inserted in
the span of each of this spans.

```julia
julia> highlighter = Highlighter(Dict(
           "keyword" => "class=keyword",
       ))
Highlighter(Language[])
```
"""
function Highlighter(names_attributes::Dict{String,String})
    names = collect(keys(names_attributes))
    attributes = collect(values(names_attributes))
    Highlighter(names, attributes)
end

function Highlighter(names, attributes)
    @assert length(names) == length(attributes) "There should be an attribute for each name (got $(length(names)) names and $(length(attributes))) attributes)."
    res = LTSH.ts_highlighter_new(names, attributes, length(names))
    if res == C_NULL
        error("Failed to create an Highlighter")
    end
    Highlighter(names, attributes, res)
end

function Base.show(io::IO, highlighter::Highlighter)
    write(io, "Highlighter(")
    show(io, highlighter.languages)
    write(io, ")")
end

Base.unsafe_convert(::Type{Ptr{LTSH.TSHighlighter}}, highlighter::Highlighter) =
    highlighter.ptr

"""
    add_language!(
        highlighter::Highlighter,
        language::Language;
        scope::String,
        highlights_query::Union{Nothing,String}=nothing,
        injections_query::Union{Nothing,String}=nothing,
        locals_query::Union{Nothing,String} = nothing,
        injection_regex::String = string("^", language.name),
    )

Adds a language definition to the given highlighter.
"""
function add_language!(
    highlighter::Highlighter,
    language::Language;
    scope::String,
    highlights_query::Union{Nothing,String} = nothing,
    injections_query::Union{Nothing,String} = nothing,
    locals_query::Union{Nothing,String} = nothing,
    injection_regex::String = string("^", language.name),
)
    push!(highlighter.languages, language)
    # lang = Ref{TSLanguage}(language.language)

    @tscall LTSH.ts_highlighter_add_language(
        highlighter,
        scope,
        injection_regex,
        language.language,
        something(highlights_query, C_NULL),
        something(injections_query, C_NULL),
        something(locals_query, C_NULL),
        highlights_query === nothing ? 0 : length(highlights_query),
        injections_query === nothing ? 0 : length(injections_query),
        locals_query === nothing ? 0 : length(locals_query),
    )

    highlighter
end

mutable struct HighlightBuffer
    ptr::Ptr{LTSH.TSHighlightBuffer}

    function HighlightBuffer(ptr)
        finalizer(LTSH.ts_highlight_buffer_delete, new(ptr))
    end
end
function HighlightBuffer()
    res = LTSH.ts_highlight_buffer_new()
    if res == C_NULL
        error("Could not create HighlightBuffer")
    end
    HighlightBuffer(res)
end

Base.unsafe_convert(::Type{Ptr{LTSH.TSHighlightBuffer}}, buffer::HighlightBuffer) =
    buffer.ptr
function Base.convert(::Type{String}, buffer::HighlightBuffer)
    ptr = convert(Ptr{Cchar}, LTSH.ts_highlight_buffer_content(buffer))
    len = LTSH.ts_highlight_buffer_len(buffer)
    unsafe_string(ptr, len)
end

"""
    highlight(highlighter::Highlighter, source_code::String; scope::String)::String

Highlights the given source code in scope and returns an HTML string.
"""
function highlight(highlighter, source_code; scope)
    buffer = HighlightBuffer()
    @tscall LTSH.ts_highlighter_highlight(
        highlighter,
        scope,
        source_code,
        length(source_code),
        buffer,
        C_NULL,
    )

    # output_line_count = LTSH.ts_highlight_buffer_line_count(buffer)
    # output_line_offsets = LTSH.ts_highlight_buffer_line_offsets(buffer)
    # output_line_offsets = Base.unsafe_wrap(Array{Cuint}, output_line_offsets, output_line_count)

    # TODO: look into string indexing
    # lines = [
    #      output_string[start:end_]
    #      for (start, end_) in zip(output_line_offsets[begin:end-1], output_line_offsets[begin+1:end])
    # ]

    # for (i, line) in enumerate(lines)
    #     println("line $i #", line)
    # end

    convert(String, buffer)
end

export add_language!, Highlighter, highlight, Language

end # module


text:
############################################################################################
#Import
include("tuning/tuning.jl")
include("utility.jl")
include("particles.jl")
include("trajectories.jl")
include("algorithm.jl")
include("postprocessing.jl")

############################################################################################
#export


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    GraphMol, QueryMol,
    graphmol, remapnodes,
    querymol,
    SDFile, SMILES, SMARTS,
    getatom, getbond, hasbond,
    setatom!, setbond!, addatom!, addbond!,
    atomcount, bondcount

struct GraphMol{A<:Atom,B<:Bond} <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{A}
    edgeattrs::Vector{B}
    cache::Dict{Symbol,Any}
    attributes::Dict{Symbol,Any}
end

"""
    graphmol() -> GraphMol

Generate empty `GraphMol`.
"""
graphmol(::Type{A}, ::Type{B}
    ) where {A<:Atom,B<:Bond} = GraphMol{A,B}([], [], [], [], Dict(), Dict())

"""
    graphmol(atoms::Vector{Atom}, bonds::Vector{Bond}) -> GraphMol

Generate `GraphMol` that has the given atom objects and edge objects.
"""
function graphmol(
        edges, atoms::Vector{A}, bonds::Vector{B}) where {A<:Atom,B<:Bond}
    nbrmap = [Dict{Int,Int}() for i in 1:length(atoms)]
    edges = collect(edges)
    for (i, (u, v)) in enumerate(edges)
        nbrmap[u][i] = v
        nbrmap[v][i] = u
    end
    return GraphMol(
        nbrmap, edges, atoms, bonds, Dict{Symbol,Any}(), Dict{Symbol,Any}())
end

"""
    graphmol(mol::GraphMol; clearcache=true) -> GraphMol

Copy `GraphMol`. Cached properties will be removed.
"""
function graphmol(graph::GraphMol{A,B}; clearcache=true
        ) where {A<:Atom,B<:Bond}
    newg = graphmol(A, B)
    for nbr in graph.neighbormap
        push!(newg.neighbormap, copy(nbr))
    end
    append!(newg.edges, graph.edges)
    append!(newg.nodeattrs, graph.nodeattrs)
    append!(newg.edgeattrs, graph.edgeattrs)
    clearcache || merge!(newg.cache, graph.cache)
    merge!(newg.attributes, graph.attributes)
    return newg
end

Graph.clone(mol::GraphMol) = graphmol(mol, clearcache=false)

"""
    graphmol(mol::SubgraphView{GraphMol}) -> GraphMol

Generate a new `GraphMol` from a substructure view.

Graph property caches and attributes are not inherited.
"""
function graphmol(view::SubgraphView)
    newg = graphmol(nodeattrtype(view), edgeattrtype(view))
    nkeys = sort(collect(nodeset(view)))
    ekeys = sort(collect(edgeset(view)))
    nmap = Dict{Int,Int}()
    for (i, n) in enumerate(nkeys)
        nmap[n] = i
        push!(newg.nodeattrs, nodeattr(view, n))
        push!(newg.neighbormap, Dict())
    end
    for (i, e) in enumerate(ekeys)
        (oldu, oldv) = getedge(view, e)
        u = nmap[oldu]
        v = nmap[oldv]
        push!(newg.edges, (u, v))
        push!(newg.edgeattrs, edgeattr(view, e))
        newg.neighbormap[u][i] = v
        newg.neighbormap[v][i] = u
    end
    return newg
end

function remapnodes(graph::GraphMol, mapping::Dict{Int,Int})
    newg = graphmol(graph)
    for i in 1:nodecount(graph)
        newg.nodeattrs[mapping[i]] = graph.nodeattrs[i]
        newg.neighbormap[mapping[i]] = Dict(
            k => mapping[v] for (k, v) in graph.neighbormap[i])
    end
    for i in 1:edgecount(graph)
        u, v = graph.edges[i]
        newg.edges[i] = (mapping[u], mapping[v])
        newg.edgeattrs[i] = graph.edgeattrs[i]
    end
    return newg
end

struct QueryMol{A<:QueryAtom,B<:QueryBond} <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{A}
    edgeattrs::Vector{B}
    cache::Dict{Symbol,Any}
    attributes::Dict{Symbol,Any}
    connectivity::Vector{Vector{Int}}
end

"""
    querymol() -> QueryMol

Generate empty `QueryMol`.
"""
querymol(
    ::Type{A}, ::Type{B}
) where {A<:QueryAtom,B<:QueryBond} = QueryMol{A,B}(
    [], [], [], [], Dict(), Dict(), [])

"""
    querymol(atoms::Vector{Atom}, bonds::Vector{Bond}) -> GraphMol

Generate `QueryMol` that has the given atom objects and edge objects.
"""
function querymol(edges, atoms::Vector{A}, bonds::Vector{B},
        connectivity::Vector{Vector{Int}}) where {A<:QueryAtom,B<:QueryBond}
    nbrmap = [Dict{Int,Int}() for i in 1:length(atoms)]
    edges = collect(edges)
    for (i, (u, v)) in enumerate(edges)
        nbrmap[u][i] = v
        nbrmap[v][i] = u
    end
    return QueryMol(
        nbrmap, edges, atoms, bonds,
        Dict{Symbol,Any}(), Dict{Symbol,Any}(), connectivity)
end

# Aliases

SDFile = GraphMol{SDFileAtom,SDFileBond}
SMILES = GraphMol{SmilesAtom,SmilesBond}
SMARTS = QueryMol{SmartsAtom,SmartsBond}

function Base.getindex(graph::GraphMol, sym::Symbol)
    if haskey(graph.cache, sym)
        return graph.cache[sym]
    else
        return eval(Expr(:call, sym, graph))
    end
end

function Base.getindex(view::SubgraphView, sym::Symbol)
    if haskey(view.graph.cache, sym)
        return view.graph.cache[sym]
    else
        return eval(Expr(:call, sym, view.graph))
    end
end

Molecule = Union{GraphMol,QueryMol}
IAtom = Union{Atom,QueryAtom}
IBond = Union{Bond,QueryBond}

# Aliases

getatom(mol::Molecule, i::Int) = nodeattr(mol, i)
getbond(mol::Molecule, i::Int) = edgeattr(mol, i)
getbond(mol::Molecule, u::Int, v::Int) = edgeattr(mol, u, v)
hasbond(mol::Molecule, u::Int, v::Int) = hasedge(mol, u, v)

setatom!(mol::Molecule, i::Int, attr::IAtom) = setnodeattr!(mol, i, attr)
setbond!(mol::Molecule, i::Int, attr::IBond) = setedgeattr!(mol, i, attr)
addatom!(mol::Molecule, attr::IAtom) = addnode!(mol, attr)
addbond!(mol::Molecule, u::Int, v::Int, attr::IBond) = addedge!(mol, u, v, attr)

atomcount(mol::Molecule) = nodecount(mol)
bondcount(mol::Molecule) = edgecount(mol)


text:
module Jtalk
# julia 1.3

  struct  HTSVoiceList
    succ::Ptr{HTSVoiceList}
    path::Cstring
    name::Cstring
  end

  const dll = Sys.iswindows() ? "jtalk" : "libjtalk"

  struct JTalk

    handle::Ptr{Cvoid}
    voicesArray::Array{Dict{Symbol,String}}

    function JTalk(;voicePath = "",dicPath = "",voiceDirPath = "")
      h = ccall((:openjtalk_initialize,dll),Ptr{Cvoid},(Cstring,Cstring,Cstring,),voicePath,dicPath,voiceDirPath)
      list = ccall((:openjtalk_getHTSVoiceList,dll),Ptr{HTSVoiceList},(Ptr{Cvoid},),h)
      voicelist = list
      a = []
      while voicelist != C_NULL
        temp = Dict{Symbol,String}()
        data = unsafe_load(voicelist)
        temp[:path] = unsafe_string(data.path)
        temp[:name] = unsafe_string(data.name)
        push!(a,temp)
        if data.succ == C_NULL
          break
        end
        voicelist = data.succ
      end
      ccall((:openjtalk_clearHTSVoiceList,dll),Cvoid,(Ptr{Cvoid},Ptr{HTSVoiceList},),h,list)
      new(h,a)
    end
  end

  ### voices
  function getVoices(j::JTalk)
    return j.voicesArray
  end

  ### samplingFrequency
  function getSamplingFrequency(j::JTalk)
    return ccall((:openjtalk_getSamplingFrequency,dll),UInt,(Ptr{Cvoid},),j.handle)
  end

  function setSamplingFrequency(j::JTalk,value::UInt)
    if value < 1
      throw("sampling frequency の範囲は1以上の整数です: " * value)
    end
    ccall((:openjtalk_setSamplingFrequency,dll),Cvoid,(Ptr{Cvoid},UInt,),j.handle,value)
  end

  function getS(j::JTalk)
    return getSamplingFrequency(j)
  end

  function setS(j::JTalk,value::UInt)
    setSamplingFrequency(j,value)
  end

  ### fPeriod
  function getFperiod(j::JTalk)
    return ccall((:openjtalk_getFperiod,dll),UInt,(Ptr{Cvoid},),j.handle)
  end

  function setFperiod(j::JTalk,value::UInt)
    if value < 1
      throw("frame period の範囲は1以上の整数です: " * value)
    end
    ccall((:openjtalk_setFperiod,dll),Cvoid,(Ptr{Cvoid},UInt,),j.handle,value)
  end

  function getP(j::JTalk)
    return getFperiod(j)
  end

  function setP(j::JTalk,value::UInt)
    setFperiod(j,value)
  end

  ### alpha
  function getAlpha(j::JTalk)
    return ccall((:openjtalk_getAlpha,dll),Float64,(Ptr{Cvoid},),j.handle)
  end

  function setAlpha(j::JTalk,value::Float64)
    if value < 0.0 || value > 1.0
      throw("all-pass constant の範囲は0と1の間の浮動小数点数です: " * value)
    end
    ccall((:openjtalk_setAlpha,dll),Cvoid,(Ptr{Cvoid},Float64,),j.handle,value)
  end

  function getA(j::JTalk)
    return getAlpha(j)
  end

  function setA(j::JTalk,value::Float64)
    setAlpha(j,value)
  end

  ### beta
  function getBeta(j::JTalk)
    return ccall((:openjtalk_getBeta,dll),Float64,(Ptr{Cvoid},),j.handle)
  end

  function setBeta(j::JTalk,value::Float64)
    if value < 0.0 || value > 1.0
      throw("postfiltering coefficient の範囲は0と1の間の浮動小数点数です: " * value)
    end
    ccall((:openjtalk_setAlpha,dll),Cvoid,(Ptr{Cvoid},Float64,),j.handle,value)
  end

  function getB(j::JTalk)
    return getBeta(j)
  end

  function setB(j::JTalk,value::Float64)
    setBeta(j,value)
  end

  ### speed
  function getSpeed(j::JTalk)
    return ccall((:openjtalk_getSpeed,dll),Float64,(Ptr{Cvoid},),j.handle)
  end

  function setSpeed(j::JTalk,value::Float64)
    if value < 0.0
      throw("speech speed rate の範囲は0以上の浮動小数点数です: " * value)
    end
    ccall((:openjtalk_setSpeed,dll),Cvoid,(Ptr{Cvoid},Float64,),j.handle,value)
  end

  function getR(j::JTalk)
    return getSpeed(j)
  end

  function setR(j::JTalk,value::Float64)
    setSpeed(j,value)
  end

  ### additionalHalfTone
  function getAdditionalHalfTone(j::JTalk)
    return ccall((:openjtalk_getAdditionalHalfTone,dll),Float64,(Ptr{Cvoid},),j.handle)
  end

  function setAdditionalHalfTone(j::JTalk,value::Float64)
    ccall((:openjtalk_setAlpha,dll),Cvoid,(Ptr{Cvoid},Float64,),j.handle,value)
  end

  function getFm(j::JTalk)
    return getAdditionalHalfTone(j)
  end

  function setFm(j::JTalk,value)
    setAdditionalHalfTone(value::Float64)
  end

  ### msdThreshold
  function getMsdThreshold(j::JTalk)
    return ccall((:openjtalk_getMsdThreshold,dll),Float64,(Ptr{Cvoid},),j.handle)
  end

  function setMsdThreshold(j::JTalk,value::Float64)
    if value < 0.0 || value > 1.0
      throw("voiced/unvoiced threshold の範囲は0と1の間の浮動小数点数です: " * value)
    end
    ccall((:openjtalk_setAlpha,dll),Cvoid,(Ptr{Cvoid},Float64,),j.handle,value)
  end

  function getU(j::JTalk)
    return getMsdThreshold(j)
  end

  function setU(j::JTalk,value::Float64)
    setMsdThreshold(j,value)
  end

  ### gvWeightForSpectrum
  function getGvWeightForSpectrum(j::JTalk)
    return ccall((:openjtalk_getGvWeightForSpectrum,dll),Float64,(Ptr{Cvoid},),j.handle)
  end

  function setGvWeightForSpectrum(j::JTalk,value::Float64)
    if value < 0.0
      throw("weight of GV for spectrum の範囲は0以上の浮動小数点数です: " * value)
    end
    ccall((:openjtalk_setAlpha,dll),Cvoid,(Ptr{Cvoid},Float64,),j.handle,value)
  end

  function getJm(j::JTalk)
    return getGvWeightForSpectrum(j)
  end

  function setJm(j::JTalk,value::Float64)
    setGvWeightForSpectrum(j,value)
  end

  ### gvWeightForLogF0
  function getGvWeightForLogF0(j::JTalk)
    return ccall((:openjtalk_getGvWeightForLogF0,dll),Float64,(Ptr{Cvoid},),j.handle)
  end

  function setGvWeightForLogF0(j::JTalk,value::Float64)
    if value < 0.0
      throw("weight of GV for LogF0 の範囲は0以上の浮動小数点数です: " * value)
    end
    ccall((:openjtalk_setGvWeightForLogF0,dll),Cvoid,(Ptr{Cvoid},Float64,),j.handle,value)
  end

  function getJf(j::JTalk)
    return getGvWeightForLogF0(j)
  end

  function setJf(j::JTalk,value::Float64)
    setGvWeightForLogF0(j,value)
  end

  ### value
  function getVolume(j::JTalk)
    return ccall((:openjtalk_getVolume,dll),Float64,(Ptr{Cvoid},),j.handle)
  end

  function setVolume(j::JTalk,value::Float64)
    ccall((:openjtalk_setAlpha,dll),Cvoid,(Ptr{Cvoid},Float64,),j.handle,value)
  end

  function getG(j::JTalk)
    return getVolume(j)
  end

  function setG(j::JTalk,value::Float64)
    setVolume(j,value)
  end

  ### dic
  function getDic(j::JTalk)
    buffer = Vector{UInt8}(undef, 260)
    res = ccall((:openjtalk_getDic,dll),Cstring,(Ptr{Cvoid},Ptr{UInt8},),j.handle,buffer)
    return unsafe_string(res)
  end

  function setDic(j::JTalk,value)
    if value == ""
      throw("辞書フォルダを示す文字列が空です")
    end
    if !isfile(value)
      throw("辞書フォルダが見つかりません: " * value)
    end
    res = ccall((:openjtalk_setDic,dll),UInt8,(Ptr{Cvoid},Ptr{UInt8},),j.handle,value)
    if res==0
      throw("辞書フォルダを設定できません。UTF-8向けの辞書ではないかもしれません: " * value)
    end
  end

  ### voiceDir
  function getVoiceDir(j::JTalk)
    buffer = Vector{UInt8}(undef, 260)
    res = ccall((:openjtalk_getVoiceDir,dll),Cstring,(Ptr{Cvoid},Ptr{UInt8},),j.handle,buffer)
    return unsafe_string(res)
  end

  function setVoiceDir(j::JTalk,value)
    if value == ""
      throw("音響モデルフォルダを示す文字列が空です")
    end
    if !isfile(value)
      throw("音響モデルフォルダが見つかりません: " * value)
    end
    res = ccall((:openjtalk_setVoiceDir,dll),Cstring,(Ptr{Cvoid},Ptr{UInt8},),j.handle,value)
    if res==0
      throw("音響モデルフォルダを設定できません: " * value)
    end
    return res
  end

  ### voice
  function getVoice(j::JTalk)
    buffer = Vector{UInt8}(undef, 260)
    res = ccall((:openjtalk_getVoice,dll),Cstring,(Ptr{Cvoid},Ptr{UInt8},),j.handle,buffer)
    result = Dict()
    result[:path] = unsafe_string(res)
    result[:name] = splitext(basename(result[:path]))[1]
    return result
  end

  function setVoice(j::JTalk,value)
    if typeof(value)==Dict{Symbol,String} && haskey(value,:path)
      v = value[:path]
    else
      v = value
    end
    if typeof(v)!=String
      throw("voiceに代入できるのは、文字列かキーpathのある Hash だけです。")
    end
    if v == ""
      throw("音響モデルを示す文字列が空です")
    end
    res = ccall((:openjtalk_setVoice,dll),UInt8,(Ptr{Cvoid},Ptr{UInt8},),j.handle,v)
    if res==0
      throw("音響モデルを設定できません: " * v)
    end
  end

  ### Speak
  function speakAsync(j::JTalk,text)
    ccall((:openjtalk_speakAsync,dll),Cvoid,(Ptr{Cvoid},Ptr{UInt8},),j.handle,text)
  end

  function speakSync(j::JTalk,text)
    ccall((:openjtalk_speakSync,dll),Cvoid,(Ptr{Cvoid},Ptr{UInt8},),j.handle,text)
  end

  function speakToFile(j::JTalk,text,path)
    ccall((:openjtalk_speakToFile,dll),Cvoid,(Ptr{Cvoid},Ptr{UInt8},Ptr{UInt8},),j.handle,text,path)
  end

  ### pause
  function pause(j::JTalk)
    ccall((:openjtalk_pause,dll),Cvoid,(Ptr{Cvoid},),j.handle)
  end

  ### resume
  function resume(j::JTalk)
    ccall((:openjtalk_resume,dll),Cvoid,(Ptr{Cvoid},),j.handle)
  end

  ### Stop
  function stop(j::JTalk)
    ccall((:openjtalk_stop,dll),Cvoid,(Ptr{Cvoid},),j.handle)
  end

  ### isSpeaking
  function isSpeaking(j::JTalk)
    ccall((:openjtalk_isSpeaking,dll),Cint,(Ptr{Cvoid},),j.handle)
  end

  ### isPaused
  function isPaused(j::JTalk)
    ccall((:openjtalk_isPaused,dll),Cint,(Ptr{Cvoid},),j.handle)
  end

  ### isFinished
  function isFinished(j::JTalk)
    ccall((:openjtalk_isFinished,dll),Cint,(Ptr{Cvoid},),j.handle)
  end

  ### Wait
  function waitUntilDone(j::JTalk)
    ccall((:openjtalk_waitUntilDone,dll),Cvoid,(Ptr{Cvoid},),j.handle)
  end

  function wait(j::JTalk,value::UInt=0)
    if value==0
      ccall((:openjtalk_wait,dll),Cvoid,(Ptr{Cvoid},UInt64,),j.handle,value)
    else
      ccall((:openjtalk_waitUntilDone,dll),Cvoid,(Ptr{Cvoid},),j.handle)
    end
  end
end


text:
using Test, Astro

print("12.a Sidereal time (mean)")
# Meeus page 88
# Find mean and apparent sidereal  time at Greenwich on 1987 April 10 at 0 UT:
# Mean sidereal time should be 13h10m46.3668s
# Nutation in right ascension should be -0.2317 degrees
# Apparent sidereal time is 13h10m46.1351s
jd = cal_to_jd(1987, 4, 10) + hms_to_fday(0,0,0)
theta0 = mean_sidereal_time_greenwich(jd)
# theta0 is mean sidereal time in radians (2pi radians = 24 hrs)

theta0_st = fday_to_hms(theta0 / (pi * 2)) # is a tuple!
theta_fday = hms_to_fday(theta0_st...)     # so splat required... :)
# mean sidereal time
@test isapprox(theta0 / (pi * 2), theta_fday, atol=0.1)
@test(theta0_st[1] == 13)
@test(theta0_st[2] == 10)
@test isapprox(theta0_st[3], 46.3668, atol=0.1)
println(" passed")

print("12.a Sidereal time (apparent)")
jd = cal_to_jd(1987, 4, 10) + hms_to_fday(0,0,0)
theta0 = apparent_sidereal_time_greenwich(jd)
# apparent sidereal time
ast_hms = fday_to_hms(theta0/(2 * pi))
@test(ast_hms[1] == 13)
@test(ast_hms[2] == 10)
@test isapprox(ast_hms[3], 46.1351, atol=0.1)
println(" passed")

print("12.b Sidereal time (mean)")
# find the mean sidereal time at Greenwich on 1987 April 10 at 19h21m00s
# jd = cal_to_jd(1987, 4, 10) + hms_to_fday(19,21,0)
msd_midnight = mean_sidereal_time_greenwich(cal_to_jd(1987, 4, 10))
println(msd_midnight)

@test isapprox(theta0 / (pi * 2), 128.7378734 / 360, atol=0.2)
println(" passed")

println(" Mean and apparent sidereal time")
mst_g = fday_to_hms(mean_sidereal_time_greenwich(cal_to_jd(1987, 4, 10)) / (2 * pi))
@test(mst_g[1] == 13)
@test(mst_g[2] == 10)
@test isapprox(mst_g[3] , 46.3668, atol=0.01)

ast_g = fday_to_hms(apparent_sidereal_time_greenwich(cal_to_jd(1987, 4, 10)) / (2 * pi))
@test(ast_g[1] == 13)
@test(ast_g[2] == 10)
@test ast_g[3] ≈  46.1351 atol= 0.3#  good enough?

# mean sidereal time and gmt/ut should be the same around September 21/22
mst = fday_to_hms(mod2pi(mean_sidereal_time_greenwich(cal_to_jd(2014,9,21)))/(2 * pi))
@test(mst[1] == 23)
@test(mst[2] == 59)
println(" passed")

# Meeus Example 12.a
# mean_sidereal_time_greenwich(Dates.datetime2julian(DateTime(1987, 4, 10, 0, 0, 0))) |> radianstime_to_fday |> fday_to_hms
# (13,10,46.366821793410054)

#=

function gmst(dt::DateTime)
    jd = Dates.datetime2julian(dt)
    theta0 = mean_sidereal_time_greenwich(jd) # theta0 is mean sidereal time in radians (2pi radians = 24 hrs)
    theta0_st = fday_to_hms(theta0 / (pi * 2)) # is a tuple!
end

# from https://www.cv.nrao.edu/~rfisher/Ephemerides/times.html

function gmst1(dt::DateTime)
    jd  = Dates.datetime2julian(dt)
    d =   jd - 2451545.0
    T =  d/36525
    gmst = 24110.54841 + 8640184.812866T + 0.093104(T^2) - 0.0000062(T^3)
    # difference from UT1
    gmst += (Dates.hour(dt) * 60 * 60) + (Dates.minute(dt) * 60) + Dates.second(dt)
    return fday_to_hms((gmst % 86400) / 86400)
end

=#