text:
printstyled("  units\n", color=:light_green)

# Cases from randSeisChannel
printstyled("    check that RandSeis uses valid UCUM units\n", color=:light_green)
redirect_stdout(out) do
  u = join(SeisIO.RandSeis.irregular_units, '.')
  @test(vucum(u, v=2))
  for u in SeisIO.RandSeis.irregular_units
    @test units2ucum(u) == u
  end
end
S = randSeisData()
@test validate_units(S) == trues(S.n)
S = randSeisData(2)
S.units = [S.units[1], "FOO"]
@test validate_units(S) == [true, false]

C = S[1]
@test validate_units(C)
@test_logs (:warn, "Error thrown for unit string: meters per second") SeisIO.vucum("meters per second")
@test vucum("m/s^2") == false

# Now some real units that show up constantly in web requests
printstyled("    test conversions to UCUM units\n", color=:light_green)

@test units2ucum("m/s^1")     == "m/s"
@test units2ucum("m/s**1")    == "m/s"
@test units2ucum("m s^-1")    == "m/s"
@test units2ucum("m*s^-1")    == "m/s"
@test units2ucum("m.s^-1")    == "m/s"
@test units2ucum("m⋅s^-1")    == "m/s"
@test units2ucum("m s**-1")   == "m/s"
@test units2ucum("m*s**-1")   == "m/s"
@test units2ucum("m.s**-1")   == "m/s"
@test units2ucum("m⋅s**-1")   == "m/s"
@test units2ucum("m × s^-1")  == "m/s"
@test units2ucum("m ⋅ s^-1")  == "m/s"

@test units2ucum("nm/s")      == "nm/s"
@test units2ucum("nm/s^1")    == "nm/s"
@test units2ucum("nm/s**1")   == "nm/s"
@test units2ucum("nm s^-1")   == "nm/s"
@test units2ucum("nm*s^-1")   == "nm/s"
@test units2ucum("nm.s^-1")   == "nm/s"
@test units2ucum("nm.s**-1")  == "nm/s"
@test units2ucum("nm*s**-1")  == "nm/s"
@test units2ucum("nm s**-1")  == "nm/s"

# Acceleration
@test units2ucum("m/s^2")   == "m/s2"
@test units2ucum("m/s**2")  == "m/s2"
@test units2ucum("m s^-2")  == "m/s2"
@test units2ucum("m*s^-2")  == "m/s2"
@test units2ucum("m.s^-2")  == "m/s2"
@test units2ucum("m⋅s^-2")  == "m/s2"
@test units2ucum("m s**-2") == "m/s2"
@test units2ucum("m*s**-2") == "m/s2"
@test units2ucum("m.s**-2") == "m/s2"
@test units2ucum("m⋅s**-2") == "m/s2"

# Fictitious cases; we're unlikely to see anything like these
@test units2ucum("s^43/s**3") == "s43/s3"
@test units2ucum("s**10/m^4") == "s10/m4"
@test units2ucum("s^43 m**-3") == "s43/m3"
@test units2ucum("s^43 cm^-3") == "s43/cm3"
@test units2ucum("s^43 cm^-3 g") == "s43/cm3.g"
@test units2ucum("s^43 cm**+3 g") == "s43.cm3.g"
@test units2ucum("m s^-1 K") == "m/s.K"
@test units2ucum("m s^-1*K") == "m/s.K"
@test units2ucum("s^43 cm^-3 V^-4") == "s43/cm3.V4"
@test units2ucum("s^43.cm^-3*V^-4") == "s43/cm3.V4"


text:
module DWDDynam

using Catlab.WiringDiagrams.DirectedWiringDiagrams
using  Catlab.CategoricalAlgebra
using  Catlab.CategoricalAlgebra.FinSets
import Catlab.WiringDiagrams: oapply

import ..UWDDynam: nstates, eval_dynamics, euler_approx

using OrdinaryDiffEq, DynamicalSystems
import OrdinaryDiffEq: ODEProblem
import DynamicalSystems: DiscreteDynamicalSystem

export AbstractMachine, ContinuousMachine, DiscreteMachine,
nstates, ninputs, noutputs, eval_dynamics, readout, euler_approx

using Base.Iterators
import Base: show, eltype

"""

A directed open dynamical system operating on information fo type `T`.
A machine  `m` has type signature  (`m.ninputs`, `m.outputs`).
"""
abstract type AbstractMachine{T} end

"""

An undirected open continuous system. The dynamics function `f` defines an ODE ``\\dot u(t) = f(u(t),x(t),p,t)`` where ``u`` is the state and ``x`` captures the exogenous variables.

The readout function may only depend on the state, so it must be of the form ``r(u(t))``.
"""
struct ContinuousMachine{T} <: AbstractMachine{T}
    ninputs::Int
    nstates::Int
    noutputs::Int
    dynamics::Function
    readout::Function
end

"""

A directed open discrete dynamical system. The dynamics function `f` defines a discrete update rule ``u_{n+1} = f(u_n, x_n, p, t)`` where ``u_n`` is the state and ``x_n`` is the value of the exogenous variables at the ``n``th time step.

The readout function may only depend on the state, so it must be of the form ``r(u_n)``.
"""
struct DiscreteMachine{T} <: AbstractMachine{T}
    ninputs::Int
    nstates::Int
    noutputs::Int
    dynamics::Function
    readout::Function
end

show(io::IO, vf::ContinuousMachine) = print("ContinuousMachine(ℝ^$(vf.nstates) × ℝ^$(vf.ninputs) → ℝ^$(vf.nstates))")
show(io::IO, vf::DiscreteMachine) = print("DiscreteMachine(ℝ^$(vf.nstates) × ℝ^$(vf.ninputs) → ℝ^$(vf.nstates))")
eltype(m::AbstractMachine{T}) where T = T

nstates(f::AbstractMachine) = f.nstates
ninputs(f::AbstractMachine) = f.ninputs
noutputs(f::AbstractMachine) = f.noutputs
readout(f::AbstractMachine, u::AbstractVector) = f.readout(u)

"""    eval_dynamics(m::AbstractMachine, u::AbstractVector, xs::AbstractVector, p, t)

Evaluates the dynamics of the machine `m` at state `u`, parameters `p`, and time `t`. The exogenous variables are set by `xs` which may either be a collection of functions ``x(t)`` or a collection of constant values.

The length of `xs` must equal the number of inputs to `m`.
"""
eval_dynamics(f::AbstractMachine, u::AbstractVector, xs::AbstractVector{T}, p, t::Real) where T <: Function =
    eval_dynamics(f, u, [x(t) for x in xs], p, t)

eval_dynamics(f::AbstractMachine, u::AbstractVector, xs::AbstractVector, p, t::Real)  = begin
    ninputs(f) == length(xs) || error("$xs must have length $(ninputs(f)) to set the exogenous variables.")
    f.dynamics(u, xs, p, t)
end

eval_dynamics(f::AbstractMachine, u::AbstractVector, xs::AbstractVector) =
    eval_dynamics(f, u, xs, nothing, 0)
eval_dynamics(f::AbstractMachine, u::AbstractVector, xs::AbstractVector, p) =
    eval_dynamics(f, u, xs, p, 0)

"""    euler_approx(m::ContinuousMachine, h)

Transforms a continuous machine `m` into a discrete machine via Euler's method with step size `h`. If the dynamics of `m` is given by ``\\dot{u}(t) = f(u(t),x(t),p,t)`` the the dynamics of the new discrete system is given by the update rule ``u_{n+1} = u_n + h f(u_n, x_n, p, t)``.
"""
euler_approx(f::ContinuousMachine{T}, h::Float64) where T = DiscreteMachine{T}(
    ninputs(f), nstates(f), noutputs(f),
    (u, x, p, t) -> u + h*eval_dynamics(f, u, x, p, t),
    f.readout
)

"""    euler_approx(m::ContinuousMachine)

Transforms a continuous machine `m` into a discrete machine via Euler's method where the step size is introduced as a new parameter, the last in the list of parameters.
"""
euler_approx(f::ContinuousMachine{T}) where T = DiscreteMachine{T}(
    ninputs(f), nstates(f), noutputs(f),
    (u, x, p, t) -> u + p[end]*eval_dynamics(f, u, x, p[1:end-1], t),
    f.readout
)
euler_approx(fs::Vector{ContinuousMachine{T}}, args...) where T =
    map(f->euler_approx(f,args...), fs)

euler_approx(fs::AbstractDict{S, ContinuousMachine{T}}, args...) where {S, T} =
    Dict(name => euler_approx(f, args...) for (name, f) in fs)

# Integration with ODEProblem in OrdinaryDiffEq.jl

"""    ODEProblem(m::ContinuousMachine, xs::Vector, u0::Vector, tspan)

Constructs an ODEProblem from the vector field defined by `(u,p,t) -> m.dynamics(u, x, p, t)`. The exogenous variables are determined by `xs`.
"""
ODEProblem(m::ContinuousMachine, u0::AbstractVector, xs::AbstractVector, tspan::Tuple{Real, Real}, p=nothing)  =
    ODEProblem((u,p,t) -> eval_dynamics(m, u, xs, p, t), u0, tspan, p)

ODEProblem(m::ContinuousMachine, u0::AbstractVector, x, tspan::Tuple{Real, Real}, p=nothing) =
    ODEProblem(m, u0, collect(repeated(x, ninputs(m))), tspan, p)

ODEProblem(m::ContinuousMachine{T}, u0::AbstractVector, tspan::Tuple{Real, Real}, p=nothing) where T =
    ODEProblem(m, u0, T[], tspan, p)

"""    DiscreteDynamicalSystem(m::DiscreteMachine, xs::Vector, u0::Vector, p)

Constructs an DiscreteDynamicalSystem from the equation of motion defined by
`(u,p,t) -> m.dynamics(u, x, p, t)`. The exogenous variables are determined by `xs`. Pass `nothing` in place of `p` if your system does not have parameters.
"""
DiscreteDynamicalSystem(m::DiscreteMachine, u0::AbstractVector, x,  p; t0::Int = 0) =
  DiscreteDynamicalSystem(m, u0, collect(repeated(x, ninputs(m))), p; t0=t0)

DiscreteDynamicalSystem(m::DiscreteMachine{T}, u0::AbstractVector, xs::AbstractVector, p; t0::Int = 0) where T = begin
  if nstates(m) == 1
    DiscreteDynamicalSystem1d(m, u0[1], xs, p; t0=t0)
  else
    !(T <: AbstractFloat) || error("Cannot construct a DiscreteDynamicalSystem if the type is a float")
    DiscreteDynamicalSystem(
        (u,p,t) -> SVector{nstates(m)}(eval_dynamics(m, u, xs, p, t)),
        u0, p; t0=t0)
  end
end

DiscreteDynamicalSystem(m::DiscreteMachine, u0::Real, xs::AbstractVector, p; t0::Int = 0) =
  DiscreteDynamicalSystem1D(m, u0, xs, p; t0=t0)

# if the system is 1D then the state must be represented by a number NOT by a 1D array
DiscreteDynamicalSystem1d(m::DiscreteMachine{T}, u0::Real, xs::AbstractVector, p; t0::Int = 0) where T = begin
  nstates(m) == 1 || error("The machine must have exactly 1 state")
  !(T <: AbstractFloat) || error("Cannot construct a DiscreteDynamicalSystem if the type is a float")
  DiscreteDynamicalSystem(
      (u,p,t) -> eval_dynamics(m, [u], xs, p, t)[1],
      u0, p; t0=t0)
end

DiscreteDynamicalSystem(m::DiscreteMachine, u0, p; t0::Int = 0) =
    DiscreteDynamicalSystem(m, u0, [], p; t0 = t0)

"""    fills(m::AbstractMachine, d::WiringDiagram, b::Int)

Checks if `m` is of the correct signature to fill box `b` of the  wiring diagram `d`.
"""
function fills(m::AbstractMachine, d::WiringDiagram, b::Int)
    b <= nboxes(d) || error("Trying to fill box $b, when $d has fewer than $b boxes")
    b = box_ids(d)[b]
    return ninputs(m) == length(input_ports(d,b)) && noutputs(m) == length(output_ports(d,b))
end

colimitmap!(f::Function, output, C::Colimit, input) = begin
    for (i,x) in enumerate(input)
        y = f(i, x)
        I = legs(C)[i](1:length(y))
        # length(I) == length(y) || error("colimitmap! attempting to fill $(length(I)) slots with $(length(y)) values")
        output[I] .= y
    end
    return output
end

@inline fillreadouts!(y, d, xs, Outputs, statefun) = colimitmap!(y, Outputs, xs) do i,x
    return x.readout(statefun(i))
end

@inline fillstates!(y, d, xs, States, statefun, inputfun, p, t) = colimitmap!(y, States, xs) do i, x
    return x.dynamics(statefun(i), inputfun(i), p, t)
end

@inline fillwire(w, d, readouts, Outputs) = readouts[legs(Outputs)[w.source.box - 2](w.source.port)] # FIX - box re-indexing

fillreadins!(readins, d, readouts, Outputs, Inputs, ins) = begin
    for (i,w) in enumerate(wires(d))
        if w.target.box == output_id(d)
            continue
        elseif w.source.box == input_id(d)
            readins[legs(Inputs)[w.target.box - 2](w.target.port)] += ins[w.source.port]
        else
            readins[legs(Inputs)[w.target.box - 2](w.target.port)] += fillwire(w, d, readouts, Outputs)
        end

    end
    return readins
end

"""    oapply(d::WiringDiagram, ms::Vector)

Implements the operad algebras for directed composition of dynamical systems given a
composition pattern (implemented by a directed wiring diagram `d`)
and primitive systems (implemented by a collection of
machines `ms`).

Each box of the composition pattern `d` is filled by a machine with the
appropriate type signature. Returns the composite machine.
"""
function oapply(d::WiringDiagram, xs::Vector{Machine}) where {T, Machine <: AbstractMachine{T}}
    isempty(wires(d, input_id(d), output_id(d))) || error("d is not a valid composition syntax because it has pass wires")
    nboxes(d) == length(xs)  || error("there are $nboxes(d) boxes but $length(xs) machines")
    for box in 1:nboxes(d)
        fills(xs[box], d, box) || error("$xs[box] does not fill box $box")
    end

    S = coproduct((FinSet∘nstates).(xs))
    Inputs = coproduct((FinSet∘ninputs).(xs))
    Outputs = coproduct((FinSet∘noutputs).(xs))
    ys = zeros(T, length(apex(S)))

    states(u::AbstractVector, b::Int) = u[legs(S)[b](1:xs[b].nstates)]

    v = (u::AbstractVector, ins::AbstractVector, p, t::Real) -> begin
        readouts = zeros(T, length(apex(Outputs)))
        readins = zeros(T, length(apex(Inputs)))

        get_states(b) = states(u,b)
        get_inputs(b) = view(readins, legs(Inputs)[b](:))

        fillreadouts!(readouts, d, xs, Outputs, get_states)
        fillreadins!(readins, d, readouts, Outputs, Inputs, ins)
        fillstates!(ys, d, xs, S, get_states, get_inputs, p, t)
        return ys
    end

    function readout(u::AbstractVector)
        readouts = zeros(T, length(apex(Outputs)))
        get_states(b) = states(u,b)

        fillreadouts!(readouts, d, xs, Outputs, get_states)
        r = zeros(T, length(d.output_ports))
        for w in in_wires(d, output_id(d))
            r[w.target.port] += fillwire(w, d, readouts, Outputs)
        end
        return r
    end

    return Machine(length(d.input_ports), length(apex(S)), length(d.output_ports), v, readout)

end

"""    oapply(d::WiringDiagram, m::AbstractMachine)

A version of `oapply` where each box of `d` is filled with the machine `m`.
"""
function oapply(d::WiringDiagram, x::AbstractMachine)
    oapply(d, collect(repeated(x, nboxes(d))))
end

"""    oapply(d::WiringDiagram, generators::Dict)

A version of `oapply` where `generators` is a dictionary mapping the name of each box to its corresponding machine.
"""
function oapply(d::WiringDiagram, xs::AbstractDict)
    oapply(d, [xs[box.value] for box in boxes(d)])
end

end #module

ext: jl


text:
"""
Cache of coefficient arrays for the AIM.
To each AIM problem corresponds a cache. As long as the problem doesn't change, the cache can be reused.

# Members
- `icda::Array{T,1}`: Hold the initial c data, i.e., c^i_0.
- `ccda::Array{T,1}`: Hold the coefficients for the current aim step, c^i_n.
- `pcda::Array{T,1}`: Hold the coefficients for the previous aim step, c^i_{n-1}.
- `bcda::Array{T,1}`: The work buffer used to actually compute the c coefficients in parallel.
- `idda::Array{T,1}`: Hold the initial d data, i.e., c^i_0.
- `cdda::Array{T,1}`: Hold the coefficients for the current aim step, d^i_n.
- `pdda::Array{T,1}`: Hold the coefficients for the previous aim step, d^i_{n-1}.
- `bdda::Array{T,1}`: The work buffer used to actually compute the d coefficients in parallel.
- `size::N`: The size of the arrays in the cache.
"""
struct AIMCache{N <: Unsigned, T <: Any}
    icda::Array{T,1}
    ccda::Array{T,1}
    pcda::Array{T,1}
    bcda::Array{T,1}

    idda::Array{T,1}
    cdda::Array{T,1}
    pdda::Array{T,1}
    bdda::Array{T,1}

    size::N
end

"""
    AIMCache(p::QuadraticEigenvalueProblem{N,T}) where {N <: Unsigned, T <: Number}

Create an AIMCache object suitable for Quadratic Eigenvalue Problems.

# Input
- `p::QuadraticEigenvalueProblem`: The problem data.

# Output
An `AIMCache{N,Polynomial{T}}` object.
"""
function AIMCache(p::QuadraticEigenvalueProblem{N,T}) where {N <: Unsigned, T <: Number}
    P = Polynomial{T}
    size = get_niter(p) + one(N)

    # The initial data arrays are initialized using the Taylor coefficient formula f^(k)(x0)/k!
    # f^(k)(x0) is a polynomial constructed by createPoly
    icda = [createPoly(p, i, λ0) for i in zero(N):get_niter(p)]
    ccda = zeros(P, size)
    pcda = zeros(P, size)
    bcda = zeros(P, size)

    idda = [createPoly(p, i, S0) for i in zero(N):get_niter(p)]
    cdda = zeros(P, size)
    pdda = zeros(P, size)
    bdda = zeros(P, size)

    AIMCache{N,P}(icda, ccda, pcda, bcda, idda, cdda, pdda, bdda, size)
end

"""
    AIMCache(p::NumericAIMProblem{N,T}) where {N <: Unsigned, T <: Number}

Create an AIMCache object suitable for Numeric Eigenvalue Problems.

# Input
- `p::NumericAIMProblem`: The problem data.

# Output
An `AIMCache{N,T}` object.
"""
function AIMCache(p::NumericAIMProblem{N,T}) where {N <: Unsigned, T <: Number}
    size = get_niter(p) + one(N)

    icda = zeros(T, size)
    ccda = zeros(T, size)
    pcda = zeros(T, size)
    bcda = zeros(T, size)

    idda = zeros(T, size)
    cdda = zeros(T, size)
    pdda = zeros(T, size)
    bdda = zeros(T, size)

    AIMCache{N,T}(icda, ccda, pcda, bcda, idda, cdda, pdda, bdda, size)
end

"""
    recomputeInitials(p::NumericAIMProblem{N,T}, c::AIMCache{N,T}, ω::T) where {N <: Unsigned, T <: Number}

Reevaluate (in-place) the initial data arrays. The initial data array elements are the Taylor expansion
coefficients of λ0 and S0 in the ODE variable x around x0 of order get_niter(p) at a point ω.

# Input
- `p::NumericAIMProblem`: The problem data.
- `c::AIMCache`: The problem data associated cache.
- `ω::T`: The value of the eigenvalue to evaluate the arrays in.

# Output
    nothing
"""
function recomputeInitials!(p::NumericAIMProblem{N,T}, c::AIMCache{N,T}, ω::T) where {N <: Unsigned, T <: Number}
    t = get_x0(p) + Taylor1(T, convert(Int64, get_niter(p)))
    c1 = λ0(p)(t,ω).coeffs
    c2 = S0(p)(t,ω).coeffs

    copy!(c.icda, c1)
    copy!(c.idda, c2)

    return nothing
end


text:
function OneVsAllPredict(Θ;μ=0,σ=1,ξ=0.5)

  P=function (X)

      m=size(X,1)
      s=zeros(m)
      X=(X.-μ)./σ
      X=[ones(m) X]
      p=Sigmoid(X*Θ')
      for i=1:m
           s[i]=argmax(p[i,:])
        end
        return s
    end

        return P
end


text:
module AnnealedIS

using Turing
using Distributions
using AdvancedMH
using AdvancedHMC; const AHMC = AdvancedHMC
using ReverseDiff
using ProgressMeter
using Memoization
using LinearAlgebra: I
using StatsFuns: logsumexp
using Random
using Bijectors

# Exports
export AnnealedISSampler,
    RejectionSampler,
    RejectionResample,
    AnIS,
    AnISHMC,
    AnnealedISSamplerHMC,
    sample,
    SimpleRejection,
    ais_sample,
    sample_from_prior,
    make_log_prior_density,
    make_log_joint_density,
    effective_sample_size

# TODO: Make it possible to parallelize sampling; difficulty is thread-safe random number generation.

abstract type RejectionSampler end

"""
Resample when a sample got rejected.
"""
struct RejectionResample <: RejectionSampler end

struct SimpleRejection <: RejectionSampler end

# TODO: Add types.
struct AnnealedISSampler{S<:RejectionSampler}
    prior_sampling::Function
    prior_density::Function
    joint_density::Function
    betas::Array{Float64}
    transition_kernels
    rejection_sampler::S
end

function Base.show(io::IO, anis::AnnealedISSampler)
    s = "AnnealedIS(num_betas=$(length(anis.betas)),rejection=$(anis.rejection_sampler))"
    print(io, s)
end

function AnnealedISSampler(
    prior_sampling::Function,
    prior_density::Function,
    joint_density::Function,
    transition_kernel_fn::Function,
    N::Int,
    rejection_sampler::T
) where {T<:RejectionSampler}
    betas = collect(range(0, 1, length=N+1))

    prior_sample = prior_sampling(Random.GLOBAL_RNG)
    transition_kernels = [transition_kernel_fn(prior_sample, i) for i in 1:(N-1)]

    return AnnealedISSampler(
        prior_sampling,
        prior_density,
        joint_density,
        betas,
        transition_kernels,
        rejection_sampler
    )
end

function AnnealedISSampler(
    prior_sampling::Function,
    prior_density::Function,
    joint_density::Function,
    transition_kernel_fn::Function,
    N::Int
)
    return AnnealedISSampler(
        prior_sampling,
        prior_density,
        joint_density,
        transition_kernel_fn,
        N,
        SimpleRejection()
    )
end

function AnnealedISSampler(
    prior_sampling,
    prior_density,
    joint_density,
    N::Int,
    rejection_sampler::T
) where {T<:RejectionSampler}
    return AnnealedISSampler(
        prior_sampling,
        prior_density,
        joint_density,
        get_normal_transition_kernel,
        N,
        rejection_sampler
    )
end

function AnnealedISSampler(prior_sampling, prior_density, joint_density, N::Int)
    return AnnealedISSampler(
        prior_sampling,
        prior_density,
        joint_density,
        get_normal_transition_kernel,
        N
    )
end

function AnnealedISSampler(
    model::Turing.Model,
    transition_kernel_fn::Function,
    N::Int
    )
    return AnnealedISSampler(
        rng -> sample_from_prior(rng, model),
        make_log_prior_density(model),
        make_log_joint_density(model),
        transition_kernel_fn,
        N
    )
end

function AnnealedISSampler(
    model::Turing.Model,
    transition_kernel_fn::Function,
    N::Int,
    rejection_sampler::T
) where {T<:RejectionSampler}
    return AnnealedISSampler(
        rng -> sample_from_prior(rng, model),
        make_log_prior_density(model),
        make_log_joint_density(model),
        transition_kernel_fn,
        N,
        rejection_sampler
    )
end

function AnnealedISSampler(model::Turing.Model, N::Int)
    return AnnealedISSampler(
        model,
        get_normal_transition_kernel,
        N
    )
end

function AnnealedISSampler(
    model::Turing.Model,
    N::Int,
    rejection_sampler::T
) where {T<:RejectionSampler}
    return AnnealedISSampler(
        model,
        get_normal_transition_kernel,
        N,
        rejection_sampler
    )
end

struct AnISModel <: AbstractMCMC.AbstractModel
    prior_sampling::Function
    prior_density::Function
    joint_density::Function
end

function AnISModel(m::Turing.Model)
    return AnISModel(
        rng -> sample_from_prior(rng, m),
        make_log_prior_density(m),
        make_log_joint_density(m)
    )
end

function AnISModel(m::Turing.Model, vi)
    return AnISModel(
        rng -> sample_from_prior(rng, m),
        make_log_prior_density(m, vi),
        make_log_joint_density(m, vi)
    )
end

struct AnIS{S<:RejectionSampler} <: Turing.InferenceAlgorithm
    betas::Array{Float64}
    transition_kernel_fn::Function
    rejection_sampler::S
end

function AnIS(
    transition_kernel_fn,
    N,
    rejection_sampler::S
) where {S<:RejectionSampler}
    betas = collect(range(0, 1, length=N+1))
    return AnIS(
        betas,
        transition_kernel_fn,
        rejection_sampler
    )
end

function AnIS(N)
    return AnIS(
        get_normal_transition_kernel,
        N,
        SimpleRejection()
    )
end

function Base.show(io::IO, anis::AnIS)
    s = "AnnealedIS(num_betas=$(length(anis.betas)),rejection=$(anis.rejection_sampler))"
    print(io, s)
end

function AnnealedISSampler(model::AnISModel, alg::AnIS)
    N = length(alg.betas)-1
    prior_sample = model.prior_sampling(Random.GLOBAL_RNG)
    transition_kernels = [
        alg.transition_kernel_fn(prior_sample, i) for i in 1:(N-1)
    ]
    return AnnealedISSampler(
        model.prior_sampling,
        model.prior_density,
        model.joint_density,
        alg.betas,
        transition_kernels,
        alg.rejection_sampler
    )
end

function AnnealedISSampler(model::Turing.Model, alg::AnIS)
    return AnnealedISSampler(AnISModel(model), alg)
end

function AnnealedISSampler(model::Turing.Model, alg::AnIS, vi)
    return AnnealedISSampler(AnISModel(model, vi), alg)
end

# HMC initialisation.

struct AnISHMC{AD,S<:RejectionSampler,P<:AHMC.AbstractProposal,M<:AHMC.AbstractMetric} <: Turing.Hamiltonian{AD}
    betas::Array{Float64} # Length N, betas[1] = 0, betas[N] = 1
    proposals::Array{P} # Length N-1
    metrics::Array{M} # Length N-1
    num_samples::Int # Number of samples for single transition kernel
    rejection_sampler::S
end

# Set ForwardDiffAD as default AD
AnISHMC(args...; kwargs...) = AnISHMC{Turing.Core.ForwardDiffAD{1}}(args...; kwargs...)
function AnISHMC{AD}(
    betas,
    proposal::P,
    num_samples,
    rejection_sampler::S
) where {AD, S<:RejectionSampler, P<:AHMC.AbstractProposal}
    # NOTE: For the unit euclidean metric AdvancedHMC.jl will automatically
    #       resize the metric to the appropriate dimension so we choose dim=1 as
    #       a default here.
    metric = AHMC.UnitEuclideanMetric(1)
    M = typeof(metric)
    return AnISHMC{AD}(
        betas, proposal, metric, num_samples, rejection_sampler
    )
end
function AnISHMC{AD}(
    betas,
    proposal::P,
    metric::M,
    num_samples,
    rejection_sampler::S
) where {AD, S<:RejectionSampler, P<:AHMC.AbstractProposal,M<:AHMC.AbstractMetric}
    num_betas = size(betas, 1)
    proposals = repeat([proposal], num_betas - 1)
    metrics = repeat([metric], num_betas - 1)
    return AnISHMC{AD,S,P,M}(
        betas, proposals, metrics, num_samples, rejection_sampler
    )
end
function AnISHMC{AD}(
    betas,
    proposals::Array{P,1},
    metrics::Array{M,1},
    num_samples,
    rejection_sampler::S
) where {AD, S<:RejectionSampler, P<:AHMC.AbstractProposal,M<:AHMC.AbstractMetric}
    return AnISHMC{AD,S,P,M}(
        betas, proposals, metrics, num_samples, rejection_sampler
    )
end

Turing.DynamicPPL.getspace(::AnISHMC) = ()
# NOTE: ReverseDiff will give wrong results because support for using contexts is broken.
#Turing.Core.getADbackend(alg::AnISHMC) = Turing.Core.ReverseDiffAD{true}()
#Turing.Core.getADbackend(alg::AnISHMC) = Turing.ForwardDiffAD{1}()
Turing.Core.getADbackend(alg::AnISHMC{AD,S,P,M}) where {AD,S,P,M} = AD()
Turing.Core.getchunksize(::Type{<:AnISHMC{AD,S,P,M}}) where {AD,S,P,M} = Turing.Core.getchunksize(AD)

struct AnnealedISSamplerHMC{AD,S<:RejectionSampler,P<:AHMC.AbstractProposal,M<:AHMC.AbstractMetric}
    prior_sampling::Function
    prior_density::Function
    joint_density::Function
    prior_grad::Function
    joint_grad::Function
    hamiltonians::Array{AHMC.Hamiltonian,1}
    alg::AnISHMC{AD,S,P,M}
end

function AnnealedISSamplerHMC(model::Turing.Model, alg::AnISHMC, spl)
    prior_sampling = gen_prior_sample(model)
    logprior = gen_logprior(spl.state.vi, model, spl)
    logjoint = gen_logjoint(spl.state.vi, model, spl)

    logprior_grad = gen_logprior_grad(spl.state.vi, model, spl)
    logjoint_grad = gen_logjoint_grad(spl.state.vi, model, spl)

    hamiltonians = make_hamiltonians(
        alg.betas,
        logprior,
        logjoint,
        logprior_grad,
        logjoint_grad,
        alg.metrics #length(spl.state.vi[spl])
    )
    return AnnealedISSamplerHMC(
        prior_sampling,
        logprior,
        logjoint,
        logprior_grad,
        logjoint_grad,
        hamiltonians,
        alg
    )
end

function make_hamiltonians(
    betas,
    logprior,
    logjoint,
    logprior_grad,
    logjoint_grad,
    metrics
)
    hamiltonians = Array{AHMC.Hamiltonian,1}(undef, length(betas)-1)
    # We don't need a hamiltonian for the first beta (=0).
    # TODO: We also don't need a hamiltonian for the last betas (=1).
    for (i, b) in enumerate(betas[2:end])
        density(x) = (1 - b) * logprior(x) + b * logjoint(x)
        density_grad(x) = (1 - b) .* logprior_grad(x) .+ b .* logjoint_grad(x)

        #metric = AHMC.UnitEuclideanMetric(dim)
        hamiltonians[i] = AHMC.Hamiltonian(metrics[i], density, density_grad)
    end

    return hamiltonians
end

"""
One importance sample with associated weight.
"""
# TODO: Type annotations for better performance.
# TODO: Possibly rename params into θ and add lp field to be consistent with Turing.
struct WeightedSample
    log_weight
    params
end

"""
Get the log density of the ith annealed distribution.
"""
function logdensity(sampler::AnnealedISSampler, i, x)
    prior_density = sampler.prior_density(x)
    joint_density = sampler.joint_density(x)
    beta = sampler.betas[i]

    # NOTE: if beta == 0 and joint_density == -Inf then we get NaN. That's why
    # we handle the special cases explicitly.
    if beta == 0.0
        return prior_density
    elseif beta == 1.0
        return joint_density
    else
        return (1-beta) * prior_density + beta * joint_density
    end
end

function logdensity(sampler::AnnealedISSamplerHMC, i, x)
    prior_density = sampler.prior_density(x)
    joint_density = sampler.joint_density(x)
    beta = sampler.alg.betas[i]

    # NOTE: if beta == 0 and joint_density == -Inf then we get NaN. That's why
    # we handle the special cases explicitly.
    if beta == 0.0
        return prior_density
    elseif beta == 1.0
        return joint_density
    else
        return (1-beta) * prior_density + beta * joint_density
    end
end

function logdensity_ratio(sampler::AnnealedISSampler, i, x)
    prior_density = sampler.prior_density(x)
    joint_density = sampler.joint_density(x)
    beta_i = sampler.betas[i]
    beta_ip1 = sampler.betas[i+1]

    return (beta_ip1 - beta_i) * joint_density - (beta_ip1 - beta_i) * prior_density
end

function logdensity_ratio(sampler::AnnealedISSamplerHMC, i, x)
    prior_density = sampler.prior_density(x)
    joint_density = sampler.joint_density(x)
    beta_i = sampler.alg.betas[i]
    beta_ip1 = sampler.alg.betas[i+1]

    return (beta_ip1 - beta_i) * joint_density - (beta_ip1 - beta_i) * prior_density
end

"""
Sample from the ith transition kernel.
"""
function transition_kernel(rng, sampler::AnnealedISSampler, i, x)
    # The base and target distribution do not have a transition kernel.
    @assert i > 1
    @assert i-1 <= length(sampler.transition_kernels)

    density(params) = logdensity(sampler, i, params)
    model = DensityModel(density)

    # Do five steps of Metropolis-Hastings
    # TODO: Possible use an iterator approach from AbstractMCMC to avoid saving
    # unused samples. NOTE: Requires AbstractMCMC 2
    spl = AdvancedMH.RWMH(sampler.transition_kernels[i-1])
    #num_steps = 20
    num_steps = 5
    #num_steps = i > 2 ? 20 : 1000
    samples = sample(model, spl, num_steps; progress=false, init_params=x)

    return samples[end].params
end

function transition_kernel(rng, sampler::AnnealedISSamplerHMC, i, x)
    hamiltonian = sampler.hamiltonians[i-1]
    proposal = sampler.alg.proposals[i-1]

    samples, stats = sample(
        hamiltonian, proposal, x, sampler.alg.num_samples; verbose=false)
    average_acceptance_rate = mean(map(s -> s.acceptance_rate, stats))
    return samples[end], average_acceptance_rate
end

function resample(rng, sampler::AnnealedISSampler{SimpleRejection}, num_rejected)
    sample = sampler.prior_sampling(rng)
    log_prior = logdensity(sampler, 1, sample)
    log_transition = logdensity(sampler, 2, sample)
    if (log_prior == -Inf || log_transition == -Inf)
        # Reject sample.
        return WeightedSample(-Inf, sample), true, num_rejected+1
    else
        log_weight = log_transition - log_prior
        return WeightedSample(log_weight, sample), false, num_rejected
    end
end

function resample(
    rng,
    sampler::AnnealedISSamplerHMC{AD,SimpleRejection,P},
    num_rejected
) where {AD,P<:AHMC.AbstractProposal}
    sample = sampler.prior_sampling()
    log_prior = logdensity(sampler, 1, sample)
    log_transition = logdensity(sampler, 2, sample)
    if (log_prior == -Inf || log_transition == -Inf)
        # Reject sample.
        return WeightedSample(-Inf, sample), true, num_rejected+1
    else
        log_weight = log_transition - log_prior
        return WeightedSample(log_weight, sample), false, num_rejected
    end
end

function resample(rng, sampler::AnnealedISSampler{RejectionResample}, num_rejected)
    sample = sampler.prior_sampling(rng)
    log_prior = logdensity(sampler, 1, sample)
    log_transition = logdensity(sampler, 2, sample)
    if (log_prior == -Inf || log_transition == -Inf)
        # Reject sample.
        return resample(rng, sampler, num_rejected+1)
    else
        log_weight = log_transition - log_prior
        return WeightedSample(log_weight, sample), false, num_rejected
    end
end

function single_sample(
    rng,
    sampler::AnnealedISSampler;
    store_intermediate_samples=false
)
    N = length(sampler.betas)
    num_samples = N-1

    # TODO: Make this a named tuple.
    diagnostics = Dict()

    weighted_sample, early_return, num_rejected = resample(rng, sampler, 0)
    diagnostics[:num_rejected] = num_rejected
    log_weight, sample = weighted_sample.log_weight, weighted_sample.params

    if store_intermediate_samples
        # FIXME: Number of annealing distributions has to be multiple of 10.
        # NOTE: We store the first 10 annealing samples and after that every 10th.
        diagnostics[:intermediate_samples] = Array{typeof(sample)}(
            undef,
            Int(num_samples/10)+9
        )
        diagnostics[:intermediate_log_weights] = Array{typeof(log_weight)}(
            undef,
            Int(num_samples/10)+9
        )
        diagnostics[:intermediate_samples][1] = sample
        diagnostics[:intermediate_log_weights][1] = log_weight
    end
    if early_return
        return weighted_sample, diagnostics
    end

    for i in 2:num_samples
        sample = transition_kernel(rng, sampler, i, sample)
        log_weight += logdensity_ratio(sampler, i, sample)
        if store_intermediate_samples && (i % 10 == 0 || i < 10)
            ix = i < 10 ? i : Int(i / 10) + 9
            diagnostics[:intermediate_samples][ix] = sample
            diagnostics[:intermediate_log_weights][ix] = log_weight
        end
    end

    ws = WeightedSample(log_weight, sample)
    return ws, diagnostics
end

function single_sample(
    rng,
    sampler::AnnealedISSamplerHMC;
    store_intermediate_samples=false
)
    N = length(sampler.alg.betas)
    num_samples = N-1

    # TODO: Make this a named tuple.
    diagnostics = Dict()

    weighted_sample, early_return, num_rejected = resample(rng, sampler, 0)
    diagnostics[:num_rejected] = num_rejected
    log_weight, sample = weighted_sample.log_weight, weighted_sample.params

    if store_intermediate_samples
        # FIXME: Number of annealing distributions has to be multiple of 10.
        # NOTE: We store the first 10 annealing samples and after that every 10th.
        diagnostics[:intermediate_samples] = Array{typeof(sample)}(
            undef,
            Int(num_samples/10)+9
        )
        diagnostics[:intermediate_log_weights] = Array{typeof(log_weight)}(
            undef,
            Int(num_samples/10)+9
        )
        diagnostics[:intermediate_acceptance_rate] = Array{Float64}(
            undef,
            Int(num_samples/10)+9
        )
        diagnostics[:intermediate_samples][1] = sample
        diagnostics[:intermediate_log_weights][1] = log_weight
        # TODO: If we use rejection sampling in first step it might not be 1.0.
        diagnostics[:intermediate_acceptance_rate][1] = 1.0
    end
    if early_return
        return weighted_sample, diagnostics
    end

    acceptance_rates = Array{Float64}(undef, num_samples-1)
    for i in 2:num_samples
        sample, acceptance_rate = transition_kernel(rng, sampler, i, sample)
        acceptance_rates[i-1] = acceptance_rate

        log_weight += logdensity_ratio(sampler, i, sample)
        if store_intermediate_samples && (i % 10 == 0 || i < 10)
            ix = i < 10 ? i : Int(i / 10) + 9
            diagnostics[:intermediate_samples][ix] = sample
            diagnostics[:intermediate_log_weights][ix] = log_weight
            diagnostics[:intermediate_acceptance_rate][ix] = acceptance_rate
        end
    end
    diagnostics[:acceptance_rate] = mean(acceptance_rates)

    ws = WeightedSample(log_weight, sample)
    return ws, diagnostics
end

function ais_sample(
    rng,
    sampler,
    num_samples;
    store_intermediate_samples=false,
    progress=false
)
    first_sample, first_diagnostics = single_sample(
        rng,
        sampler;
        store_intermediate_samples=store_intermediate_samples
    )
    samples = Array{typeof(first_sample)}(undef, num_samples)
    samples[1] = first_sample

    diagnostics = init_diagnostics(first_diagnostics, num_samples)

    if progress
        pmeter = Progress(num_samples)
    end
    spls = [deepcopy(sampler) for _ in 1:Threads.nthreads()]
    Threads.@threads for i = 2:num_samples
        tid = Threads.threadid()
        samples[i], d = single_sample(
            rng,
            spls[tid];
            store_intermediate_samples=store_intermediate_samples
        )
        add_diagnostics!(diagnostics, d, i)
        if progress
            next!(pmeter)
        end
    end

    diagnostics[:ess] = effective_sample_size(samples)

    return samples, diagnostics
end

function ais_sample(
    rng,
    model::AnISModel,
    sampler::AnIS,
    num_samples;
    kwargs...
)
    ais = AnnealedISSampler(model, sampler)
    return ais_sample(rng, ais, num_samples; kwargs...)
end

function estimate_expectation(samples::Array{WeightedSample}, f)
    weighted_sum = sum(samples) do weighted_sample
        exp(weighted_sample.log_weight) * f(weighted_sample.params)
    end
    sum_weights = sum(samples) do weighted_sample
        exp(weighted_sample.log_weight)
    end
    return weighted_sum / sum_weights
end

function effective_sample_size(samples::Array{WeightedSample}; normalise=true)
    log_weights = map(s -> s.log_weight, samples)
    denominator = logsumexp(2 * log_weights)
    numerator = 1
    if normalise
        numerator = 2 * logsumexp(log_weights)
    end

    return exp(numerator - denominator)
end

##############################
# Transition kernels
##############################

# TODO: Is this an appropriate type annotation?
function get_normal_transition_kernel(prior_sample::T, i) where {T<:Real}
    return Normal(0, 1)
end

function get_normal_transition_kernel(prior_sample::AbstractArray, i)
    return MvNormal(size(prior_sample, 1), 1)
end

function get_normal_transition_kernel(prior_sample::NamedTuple, i)
    return map(x -> get_normal_transition_kernel(x, i), prior_sample)
end

# TODO: Possibly add this to AdvancedMH.jl
AdvancedMH.RWMH(nt::NamedTuple) = MetropolisHastings(map(x -> RandomWalkProposal(x), nt))

##############################
# Utility functions
##############################

function add_diagnostics!(diagnostics, d, i)
    for k in keys(d)
        if !haskey(diagnostics, k)
            # Not all keys are initialized because some samples might have been
            # rejected.
            first_k = iterate(keys(diagnostics))[1]
            diagnostics[k] = Array{typeof(d[k])}(
                undef, size(diagnostics[first_k], 1))
        end
        diagnostics[k][i] = d[k]
    end
end

function init_diagnostics(first_diagnostics, num_samples)
    diagnostics = Dict()
    for k in keys(first_diagnostics)
        diagnostics[k] = Array{typeof(first_diagnostics[k])}(
        undef, num_samples)
        diagnostics[k][1] = first_diagnostics[k]
    end

    return diagnostics
end

include("turing_interop.jl")

####
#### Compiler interface, i.e. tilde operators.
####
#### Taken from https://github.com/TuringLang/Turing.jl/blob/master/src/inference/hmc.jl
function DynamicPPL.assume(
    rng,
    spl::Sampler{<:AnISHMC},
    dist::Distribution,
    vn::VarName,
    vi,
)
    updategid!(vi, vn, spl)
    r = vi[vn]
    # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))
    # r
    return r, Bijectors.logpdf_with_trans(dist, r, istrans(vi, vn))
end

function DynamicPPL.dot_assume(
    rng,
    spl::Sampler{<:AnISHMC},
    dist::MultivariateDistribution,
    vns::AbstractArray{<:VarName},
    var::AbstractMatrix,
    vi,
)
    @assert length(dist) == size(var, 1)
    updategid!.(Ref(vi), vns, Ref(spl))
    r = vi[vns]
    var .= r
    return var, sum(Bijectors.logpdf_with_trans(dist, r, istrans(vi, vns[1])))
end
function DynamicPPL.dot_assume(
    rng,
    spl::Sampler{<:AnISHMC},
    dists::Union{Distribution, AbstractArray{<:Distribution}},
    vns::AbstractArray{<:VarName},
    var::AbstractArray,
    vi,
)
    updategid!.(Ref(vi), vns, Ref(spl))
    r = reshape(vi[vec(vns)], size(var))
    var .= r
    return var, sum(Bijectors.logpdf_with_trans.(dists, r, istrans(vi, vns[1])))
end

function DynamicPPL.observe(
    spl::Sampler{<:AnISHMC},
    d::Distribution,
    value,
    vi,
)
    return DynamicPPL.observe(SampleFromPrior(), d, value, vi)
end

function DynamicPPL.dot_observe(
    spl::Sampler{<:AnISHMC},
    ds::Union{Distribution, AbstractArray{<:Distribution}},
    value::AbstractArray,
    vi,
)
    return DynamicPPL.dot_observe(SampleFromPrior(), ds, value, vi)
end

end


text:
# Run with: julia --color=yes -L ../../src/BlackBoxOptim.jl compare_optimizers.jl ...

using DataFrames
using BlackBoxOptim
using ArgParse
using CPUTime
using Compat
using Statistics
using Printf: @sprintf

global logfilehandle = nothing
function log(color::Symbol, str)
    printstyled(str, color=color)
    if logfilehandle != nothing
        print(logfilehandle, str)
        flush(logfilehandle)
    end
end
log(str::AbstractString) = log(:white, str)

function main(args)

    start_time = time()
    global logfilehandle
    logfilehandle = open(Libc.strftime("%Y%m%d_%H%M%S.log", start_time), "a+")

    s = ArgParseSettings(description = "comparing optimizers to benchmark runs",
        commands_are_required = true,
        version = "0.1.0",
        add_version = true
    )

    @add_arg_table s begin
        "update"
            action = :command
            help = "update benchmark runs"

        "list"
            action = :command
            help = "list info about benchmark runs"

        "compare"
            action = :command
            help = "compare current optimizer(s) to benchmark runs"
    end

    for cmd in ["update", "compare", "list"]
        @add_arg_table s[cmd] begin
            "--benchmarkfile", "-b"
                arg_type = AbstractString
                default = "benchmark_runs.csv"
                help = "name of benchmark runs db file"
        end
    end

    for cmd in ["update", "compare"]
        @add_arg_table s[cmd] begin
           "--problems"
                arg_type = AbstractString
                default = "all"
                help = "name of problem set"

            "--optimizers"
                arg_type = AbstractString
                default = "stable"
                help = "name of optimizer or optimizer set"

            "--numreps", "-n"
                arg_type = Int
                default = 10
                help = "number of repetitions per problem and optimizer"

            "--adaptive"
                action = :store_true
                help = "adaptive number of reps until we can prove there is a difference (with a max of 30 reps)"
        end
    end

    @add_arg_table s["compare"] begin
        "--comparisonfile"
            arg_type = AbstractString
            default = ""
            help = "name of comparison db file"
        end

    @add_arg_table s["list"] begin
        "--outfile", "-o"
            arg_type = AbstractString
            default = ""
            help = "name of csv file where result comparison table is saved"
        end

    pargs = parse_args(args, s)

    cmd = pargs["%COMMAND%"]

    benchmarkfile = pargs[cmd]["benchmarkfile"]

    local pset, optimizers, nreps

    if in(cmd, ["update", "compare"])
        if haskey(ProblemSets, pargs[cmd]["problems"])
            pset = ProblemSets[pargs[cmd]["problems"]]
        else
            ps = pargs[cmd]["problems"]
            throw("Unknown problem set $(ps), existing ones are: ", keys(ProblemSets))
        end

        if haskey(OptimizerSets, pargs[cmd]["optimizers"])
            optimizers = OptimizerSets[pargs[cmd]["optimizers"]]
        else
            os = pargs[cmd]["optimizers"]
            throw("Unknown optimizer set $(os), existing ones are: ", keys(OptimizerSets))
        end

        nreps = pargs[cmd]["numreps"]
    end

    if cmd == "update"
        db = read_benchmark_db(benchmarkfile)
        db = update_benchmarks(db, pset, optimizers, nreps)
        save_result_database(db, benchmarkfile)
    elseif cmd == "list"
        db = read_benchmark_db(benchmarkfile)
        list_benchmark_db(db, pargs[cmd]["outfile"])
    elseif cmd == "compare"
        compare_optimizers_to_benchmarks(benchmarkfile, pset, optimizers, nreps, 0.01,
                pargs[cmd]["comparisonfile"])
    else
        raise("Unknown command") # FIXME logfilehandle is not closed
    end

    close(logfilehandle)

    println("elapsed time:    ", format_time(time() - start_time))
end

function format_time(t)
    if t < 60.0
        @sprintf("%.1f secs", t)
    elseif t < 3600.0
        m = floor(Int, t/60.0)
        @sprintf("%d min%s ", m, (m>1) ? "s" : "") * format_time(t - 60*m)
    else
        h = floor(Int, t/3600.0)
        @sprintf("%d hour%s ", h, (h>1) ? "s" : "") * format_time(t - 3600*h)
    end
end

ProblemSets = Dict{String,Any}(
    "easy" => [
        # ProblemName, NumDims, PopSize, MaxFevals
        ("Sphere",        5, 20, 5e3),
        ("Sphere",       10, 20, 1e4),
        ("Sphere",       30, 20, 3e4),

        ("Schwefel2.22",  5, 20, 5e3),
        ("Schwefel2.22", 10, 20, 1e4),
        ("Schwefel2.22", 30, 20, 3e4),

        ("Schwefel2.21",  5, 20, 5e3),
        ("Schwefel2.21", 10, 20, 1e4),
        ("Schwefel2.21", 30, 20, 3e4)
    ],

    "harder" => [
        # Harder problems
        ("Schwefel1.2",   5, 20, 5e3),
        ("Schwefel1.2",  10, 50, 5e4),
        ("Schwefel1.2",  30, 50, 2e5),
        ("Schwefel1.2",  50, 50, 3e5),

        ("Rosenbrock",    5, 20, 1e4),
        ("Rosenbrock",   10, 50, 5e4),
        ("Rosenbrock",   30, 50, 2e5),
        ("Rosenbrock",   50, 40, 3e5),

        ("Rastrigin",    50, 50, 5e5),
        ("Rastrigin",   100, 90, 8e5),

        ("Ackley",       50, 50, 5e5),
        ("Ackley",      100, 90, 8e5),

        ("Griewank",     50, 50, 5e5),
        ("Griewank",    100, 90, 8e5),
    ],

    "lowdim" => [
        ("Schwefel1.2",   2, 25, 1e4),
        ("Rosenbrock",    2, 25, 1e4),
        ("Rastrigin",     2, 25, 1e4),
        ("Ackley",        2, 25, 1e4),
        ("Griewank",      2, 25, 1e4),
    ],

    "test" => [
        ("Rosenbrock",   30, 50, 2e5),
    ]
)
ProblemSets["all"] = vcat(ProblemSets["easy"], ProblemSets["harder"])

OptimizerSets = Dict{String,Any}(
    "de" => [:de_rand_1_bin, :de_rand_1_bin_radiuslimited, :adaptive_de_rand_1_bin, :adaptive_de_rand_1_bin_radiuslimited],
    "stable_non_de" => [:probabilistic_descent, :generating_set_search, :random_search],
    "nes" => [:xnes, :separable_nes, :dxnes],
    "test" => [:de_rand_1_bin],
)
OptimizerSets["all"] = unique(collect(keys(BlackBoxOptim.SingleObjectiveMethods)))
OptimizerSets["stable"] = vcat(OptimizerSets["de"], OptimizerSets["stable_non_de"])

# Keep track of how many times we have executed each problem for this run of the script,
# since we should never take time measurements before we have ensured everything has been compiled.
const RunsPerProblem = Dict{Any,Int}()
runs_per_problem(problemname, numdims) = get(RunsPerProblem, (problemname, numdims), 0)
function increase_runs_per_problem(problemname, numdims)
    k = (problemname, numdims)
    RunsPerProblem[k] = get(RunsPerProblem, k, 0) + 1
end

function fitness_for_opt(family::FunctionBasedProblemFamily, NumDimensions::Int,
                         PopulationSize::Int, MaxFuncEvals::Int,
                         Method::Symbol, TraceMode::Symbol = :verbose)

    problem = instantiate(family, NumDimensions)

    CPUtic()
    res = bboptimize(problem; Method = Method,
        NumDimensions = NumDimensions,
        PopulationSize = PopulationSize,
        TraceMode = TraceMode,
        MaxFuncEvals = MaxFuncEvals
    )

    best_fitness(res), CPUtoc()
end

function latest_git_id(git_repo = joinpath(@__DIR__, "../.."))
    strip(read(`git -C $git_repo log --format="%H" -n 1`, String))
end

# Test an optimizer on multiple problems and dimensions, return results in dict.
function multitest_opt(problemDescriptions, method; NumRepetitions = 3)
    dfs = DataFrame[]

    gitid = latest_git_id()

    for pd in problemDescriptions
        log("Running $(NumRepetitions) reps for $(method)\n")
        for rep in 1:NumRepetitions
            probname, numdims, popsize, numfevals = pd

            prob = BlackBoxOptim.example_problems[probname]

            # Ensure everything is compiled before we start measuring the 1st time.
            if runs_per_problem(probname, numdims) == 0
                fitness_for_opt(prob, numdims, popsize, 100, method, :silent)
                increase_runs_per_problem(probname, numdims)
            end

            df = DataFrame(Problem = probname, NumDims = numdims,
                    Method = string(method), PopulationSize = popsize,
                    NumFevals = numfevals,
                    GitID = gitid
            )

            log("\n$(probname), n = $(numdims), optimizer = $(string(method)), run $(rep) of $(NumRepetitions)\n")

            start_time = time()
            ftn, dur = fitness_for_opt(prob, numdims, popsize, ceil(Int, numfevals), method)
            df[:ElapsedTime] = dur
            df[:StartTime] = Libc.strftime("%Y-%m-%d %H:%M:%S", start_time)
            df[:Fitness] = ftn

            push!(dfs, df)
        end
    end

    vcat(dfs)
end

using CSV

function read_benchmark_db(filename)
    if isfile(filename)
        return CSV.read(filename)
    else
        @warn "Benchmark results file $filename not found, returning empty frame"
        return DataFrame()
    end
end

function add_rank_per_group(df, groupcols, rankcol, resultcol)
    by(df, groupcols) do subdf
        orderedsubdf = subdf[sortperm(subdf[:,rankcol]), :]
        orderedsubdf[resultcol] = collect(1:size(orderedsubdf,1))
        return orderedsubdf
    end
end

function list_benchmark_db(db, saveResultCsvFile = nothing)
    numrows = size(db, 1)
    println("Number of runs in db: ", numrows)

    if numrows > 1
        # Find min fitness per problem
        minfitpp = by(db, [:Problem, :NumDims]) do df
            DataFrame(
                MinFitness = minimum(df[:Fitness])
            )
        end

        # Add col with order of magnitude worse than min fitness for each run
        db = join(db, minfitpp, on = [:Problem, :NumDims])
        db[:LogTimesWorseFitness] = log10.(db[:Fitness] ./ db[:MinFitness])

        # Calc median fitness and time per problem and method.
        sumdf = by(db, [:Problem, :NumDims, :Method]) do df
            DataFrame(N = size(df, 1),
                      MedianFitness = median(df[:,:Fitness]),
                      MedianTime = median(df[:,:ElapsedTime]))
        end

        # Rank on median fitness and median time for each problem.
        sumdf = add_rank_per_group(sumdf, [:Problem, :NumDims], :MedianFitness, :RankFitness)
        sumdf = add_rank_per_group(sumdf, [:Problem, :NumDims], :MedianTime, :RankTime)

        # Get number of runs and median magnitude worse per method
        permethod = by(db, [:Method]) do df
            DataFrame(
                NumRuns = size(df, 1),
                MedianLogTimesWorseFitness = round(median(df[:, :LogTimesWorseFitness]), digits=1)
            )
        end

        # and merge with table with mean ranks of fitness and time.
        summarydf = by(sumdf, [:Method]) do df
            DataFrame(
                MeanRank = round(mean(df[:RankFitness]), digits=3),
                Num1sFitness = sum(r -> (r == 1) ? 1 : 0, df[:RankFitness]),
                MeanRankTime = round(mean(df[:RankTime]), digits=3),
                Num1sTime = sum(r -> (r == 1) ? 1 : 0, df[:RankTime]),
            )
        end
        df = join(summarydf, permethod, on = :Method)

        # Now sort and print
        sort!(df, [:MeanRank, :MeanRankTime, :Num1sFitness])
        println(df)
        if saveResultCsvFile !== nothing
            CSV.write(saveResultCsvFile, df)
            println("Results written to file: ", saveResultCsvFile)
        end
    end
end

function save_result_database(db, filename)
    CSV.write(filename, db)
end

function update_benchmarks(db, pset, optimizers, nreps = 10)
    for optmethod in optimizers
        @time res = multitest_opt(pset, optmethod; NumRepetitions = nreps)
        if size(db, 1) < 1
            db = res
        else
            db = vcat(db, res)
        end
        #save_result_database(db, "temp.csv")
        list_benchmark_db(db)
    end
    return db
end

using HypothesisTests

function compare_optimizers_to_benchmarks(benchmarkfile, pset, optimizers, nreps,
    significancelevel::Float64 = 0.05, comparisonfile = "")
    if comparisonfile == ""
        db = read_benchmark_db(benchmarkfile)
        dfs = DataFrame[]
        totalruns = length(optimizers) * nreps * length(pset)
        runnum = 0
        for optmethod in optimizers
            optsel = db[:Method] .== string(optmethod)
            for pd in pset
                probname, numdims, popsize, numfevals = pd
                psel = db[:Problem] .== probname
                dsel = db[:NumDims] .== numdims
                df = db[optsel .& psel .& dsel, :]
                benchfitnesses = convert(Vector{Float64}, df[:Fitness])
                benchtimes = convert(Vector{Float64}, df[:ElapsedTime])
                newfs = Float64[]
                newtimes = Float64[]
                prob = BlackBoxOptim.example_problems[probname]
                for r in 1:nreps
                    runnum += 1
                    log("\n$(probname), n = $(numdims), optimizer = $(string(optmethod)), run $(r) of $(nreps) ($(round(100.0 * runnum / totalruns, digits=2))% of total runs)\n")
                    ftn, dur = fitness_for_opt(prob, numdims, popsize, ceil(Int, numfevals), optmethod)
                    push!(newfs, ftn)
                    push!(newtimes, dur)
                end

                log(:blue, "$(probname)($numdims), $(optmethod):\n")

                ftn_oldmed = median(benchfitnesses)
                ftn_newmed = median(newfs)
                ftn_ord = (ftn_oldmed < ftn_newmed) ? "<" : ">"
                log(:white, "Fitness median $(ftn_oldmed)(old) $(ftn_ord) $(ftn_newmed)(new)\n")
                ftn_pval = pvalue(MannWhitneyUTest(benchfitnesses, newfs))
                if ftn_pval > significancelevel
                    log(:white, "  No statistically significant fitness difference in $nreps repetitions (P-value=$ftn_pval)!\n")
                else
                    log(ftn_ord == "<" ? :red : :green, "  Statistically significant fitness difference in $nreps repetitions (P-value=$ftn_pval)!\n")
                end

                time_oldmed = median(benchtimes)
                time_newmed = median(newtimes)
                time_ord = (time_oldmed < time_newmed) ? "<" : ">"
                log(:white, "Time median $(time_oldmed)(old) $(time_ord) $(time_newmed)(new)\n")
                time_pval = pvalue(MannWhitneyUTest(benchtimes, newtimes))
                if time_pval > significancelevel
                    log(:white, "  No statistically significant time difference in $nreps repetitions (P-value=$time_pval)!\n")
                else
                    log(time_ord == "<" ? :red : :green, "  Statistically significant time difference in $nreps repetitions (P-value=$time_pval)!\n")
                end

                push!(dfs, DataFrame(
                    Problem = probname, NumDims = numdims, Method = string(optmethod), PvalueThreshold = significancelevel,
                    FitnessOldMedian = ftn_oldmed, FitnessOrder = ftn_ord, FitnessNewMedian = ftn_newmed,
                    FitnessPvalue = ftn_pval, FitnessIsSignif = ftn_pval <= significancelevel,
                    TimeOldMedian = time_oldmed, TimeOrder = time_ord, TimeNewMedian = time_newmed,
                    TimePvalue = time_pval, TimeIsSignif = time_pval <= significancelevel)
                )
            end
        end
        df = vcat(dfs...)

        # Use Benjamini-Hochberg to judge which pvalues are significant given we did
        # many comparisons.
        ftn_pvs = convert(Vector{Float64}, df[:FitnessPvalue])
        df[:FitnessSignificantBH001] = benjamini_hochberg(ftn_pvs, 0.01)
        df[:FitnessSignificantBH005] = benjamini_hochberg(ftn_pvs, 0.05)
        df[:FitnessSignificantBH010] = benjamini_hochberg(ftn_pvs, 0.10)

        time_pvs = convert(Vector{Float64}, df[:TimePvalue])
        df[:TimeSignificantBH001] = benjamini_hochberg(time_pvs, 0.01)
        df[:TimeSignificantBH005] = benjamini_hochberg(time_pvs, 0.05)
        df[:TimeSignificantBH010] = benjamini_hochberg(time_pvs, 0.10)

        CSV.write(Libc.strftime("comparison_%Y%m%d_%H%M%S.csv", time()), df)
    else
        df = CSV.read(comparisonfile)
    end
    sort!(df, [:FitnessPvalue])
    report_below_pvalue(df, col_prefix="Fitness", pvalue=1.00)
    report_below_pvalue(df, col_prefix="Fitness", pvalue=0.05)
    report_below_pvalue(df, col_prefix="Fitness", pvalue=0.01)
    report_below_pvalue(df, col_prefix="Time", pvalue=1.00)
    report_below_pvalue(df, col_prefix="Time", pvalue=0.05)
    report_below_pvalue(df, col_prefix="Time", pvalue=0.01)

    # Report (in color) on number of significant differences after Benjamini-Hochberg
    # correction.
    n_ftn_reg = sum(df[:FitnessSignificantBH005] .& (df[:FitnessOrder] .== "<"))
    n_ftn_imp = sum(df[:FitnessSignificantBH005] .& (df[:FitnessOrder] .== ">"))
    printstyled("\n$n_ftn_reg significant fitness regressions at Benjamini-Hochberg 0.05 level\n",
                color=n_ftn_reg > 0 ? :red : :green)
    printstyled("\n$n_ftn_imp significant fitness improvments at Benjamini-Hochberg 0.05 level\n",
                color=n_ftn_imp > 0 ? :green : :white)
    n_time_reg = sum(df[:TimeSignificantBH005] .& (df[:TimeOrder] .== "<"))
    n_time_imp = sum(df[:TimeSignificantBH005] .& (df[:TimeOrder] .== ">"))
    printstyled("\n$n_time_reg significant time regressions at Benjamini-Hochberg 0.05 level\n",
                color=n_time_reg > 0 ? :red : :green)
    printstyled("\n$n_time_imp significant time improvments at Benjamini-Hochberg 0.05 level\n",
                color=n_time_imp > 0 ? :green : :white)
end

function benjamini_hochberg(pvals, alpha = 0.05)
    n = length(pvals)
    (n <= 1) && return pvals # no pvalue correction needed

    perm = sortperm(pvals)
    origperm = invperm(perm)

    thresholds = alpha .* (collect(1:n) ./ n)
    khat = findlast(i -> pvals[perm[i]] < thresholds[i], 1:n)
    (khat === nothing) && return fill(false, n) # There were no significant pvals so return all false

    significant = vcat(fill(true, khat), fill(false, n-khat))
    return significant[origperm]
end

function report_below_pvalue(df; col_prefix="Fitness", pvalue = 0.05)
    selection = df[Symbol(col_prefix, "Pvalue")] .< pvalue
    log("Num problems with $col_prefix p-values < $(pvalue): $(sum(selection))\n")
    # workaround sum(isequal, []) throws
    num_new_worse = sum(df[selection, Symbol(col_prefix, "Order")] .== "<")
    num_new_better = sum(df[selection, Symbol(col_prefix, "Order")] .== ">")
    log_color = num_new_better > num_new_worse ? :green : :red
    log(log_color, "  $col_prefix is better in current run in $(num_new_better) tests\n")
    log(log_color, "  $col_prefix is worse in current run in $(num_new_worse) tests\n")
    any(selection) && println(df[selection, :])
end

@CPUtime main(ARGS)


text:
module CodeTracking

using Base: PkgId
using Core: LineInfoNode
using Base.Meta: isexpr
using UUIDs
using InteractiveUtils

export whereis, definition, pkgfiles, signatures_at

include("pkgfiles.jl")
include("utils.jl")

### Global storage

# These values get populated by Revise

const method_info = IdDict{Type,Union{Missing,Tuple{LineNumberNode,Expr}}}()

const _pkgfiles = Dict{PkgId,PkgFiles}()

# Callback for method-lookup. `lookupfunc = method_lookup_callback[]` must have the form
#     ret = lookupfunc(method)
# where `ret` is either `nothing` or `(lnn, def)`. `lnn` is a LineNumberNode (or any valid
# input to `CodeTracking.fileline`) and `def` is the expression defining the method.
const method_lookup_callback = Ref{Any}(nothing)

# Callback for `signatures_at` (lookup by file/lineno). `lookupfunc = expressions_callback[]`
# must have the form
#    mod, exsigs = lookupfunc(id, relpath)
# where
#    id is the PkgId of the corresponding package
#    relpath is the path of the file from the basedir of `id`
#    mod is the "active" module at that point in the source
#    exsigs is a ex=>sigs dictionary, where `ex` is the source expression and `sigs`
#        a list of method-signatures defined by that expression.
const expressions_callback = Ref{Any}(nothing)

const juliabase = joinpath("julia", "base")
const juliastdlib = joinpath("julia", "stdlib", "v$(VERSION.major).$(VERSION.minor)")

### Public API

"""
    filepath, line = whereis(method::Method)

Return the file and line of the definition of `method`. `line`
is the first line of the method's body.
"""
function whereis(method::Method)
    file, line = String(method.file), method.line
    startswith(file, "REPL[") && return file, line
    lin = get(method_info, method.sig, nothing)
    if lin === nothing
        f = method_lookup_callback[]
        if f !== nothing
            try
                Base.invokelatest(f, method)
                lin = get(method_info, method.sig, nothing)
            catch
            end
        end
    end
    if lin === nothing || ismissing(lin)
    else
        file, line = fileline(lin[1])
    end
    file = maybe_fix_path(file)
    return file, line
end

"""
    loc = whereis(sf::StackFrame)

Return location information for a single frame of a stack trace.
If `sf` corresponds to a frame that was inlined, `loc` will be `nothing`.
Otherwise `loc` will be `(filepath, line)`.
"""
function whereis(sf::StackTraces.StackFrame)
    sf.linfo === nothing && return nothing
    return whereis(sf, sf.linfo.def)
end

"""
    filepath, line = whereis(lineinfo, method::Method)

Return the file and line number associated with a specific statement in `method`.
`lineinfo.line` should contain the line number of the statement at the time `method`
was compiled. The current location is returned.
"""
function whereis(lineinfo, method::Method)
    file, line1 = whereis(method)
    # We could be in an expanded macro. Apply the correction only if the filename checks out.
    # (We're not super-fastidious here because of symlinks and other path ambiguities)
    samefile = basename(file) == basename(String(lineinfo.file))
    if !samefile
        return maybe_fix_path(String(lineinfo.file)), lineinfo.line
    end
    return file, lineinfo.line - method.line + line1
end
function whereis(lineinfo::Core.LineInfoNode, method::Method)
    # With LineInfoNode we have certainty about whether we're in a macro expansion
    if lineinfo.method == Symbol("macro expansion")
        return maybe_fix_path(String(lineinfo.file)), lineinfo.line
    end
    file, line1 = whereis(method)
    return file, lineinfo.line - method.line + line1
end

"""
    sigs = signatures_at(filename, line)

Return the signatures of all methods whose definition spans the specified location.
`line` must correspond to a line in the method body (not the signature or final `end`).

Returns `nothing` if there are no methods at that location.
"""
function signatures_at(filename::AbstractString, line::Integer)
    if !startswith(filename, "REPL[")
        filename = abspath(filename)
    end
    if occursin(juliabase, filename)
        rpath = postpath(filename, juliabase)
        id = PkgId(Base)
        return signatures_at(id, rpath, line)
    elseif occursin(juliastdlib, filename)
        rpath = postpath(filename, juliastdlib)
        spath = splitpath(rpath)
        libname = spath[1]
        project = Base.active_project()
        id = PkgId(Base.project_deps_get(project, libname), libname)
        return signatures_at(id, joinpath(spath[2:end]...), line)
    end
    if startswith(filename, "REPL[")
        id = PkgId("@REPL")
        return signatures_at(id, filename, line)
    end
    for (id, pkgfls) in _pkgfiles
        if startswith(filename, basedir(pkgfls)) || id.name == "Main"
            bdir = basedir(pkgfls)
            rpath = isempty(bdir) ? filename : relpath(filename, bdir)
            if rpath ∈ pkgfls.files
                return signatures_at(id, rpath, line)
            end
        end
    end
    error("$filename not found in internal data, perhaps the package is not loaded (or not loaded with `includet`)")
end

"""
    sigs = signatures_at(mod::Module, relativepath, line)

For a package that defines module `mod`, return the signatures of all methods whose definition
spans the specified location. `relativepath` indicates the path of the file relative to
the packages top-level directory, e.g., `"src/utils.jl"`.
`line` must correspond to a line in the method body (not the signature or final `end`).

Returns `nothing` if there are no methods at that location.
"""
function signatures_at(mod::Module, relpath::AbstractString, line::Integer)
    id = PkgId(mod)
    return signatures_at(id, relpath, line)
end

function signatures_at(id::PkgId, relpath::AbstractString, line::Integer)
    expressions = expressions_callback[]
    expressions === nothing && error("cannot look up methods by line number, try `using Revise` before loading other packages")
    try
        for (mod, exsigs) in Base.invokelatest(expressions, id, relpath)
            for (ex, sigs) in exsigs
                lr = linerange(ex)
                lr === nothing && continue
                line ∈ lr && return sigs
            end
        end
    catch
    end
    return nothing
end

"""
    src, line1 = definition(String, method::Method)

Return a string with the code that defines `method`. Also return the first line of the
definition, including the signature (which may not be the same line number returned
by `whereis`). If the method can't be located (line number 0), then `definition`
instead returns `nothing.`

Note this may not be terribly useful for methods that are defined inside `@eval` statements;
see [`definition(Expr, method::Method)`](@ref) instead.
"""
function definition(::Type{String}, method::Method)
    file, line = whereis(method)
    line == 0 && return nothing
    src = src_from_file_or_REPL(file)
    eol = isequal('\n')
    linestarts = Int[]
    istart = 1
    for i = 1:line-1
        push!(linestarts, istart)
        istart = findnext(eol, src, istart) + 1
    end
    ex, iend = Meta.parse(src, istart)
    iend = prevind(src, iend)
    if isfuncexpr(ex)
        iend = min(iend, lastindex(src))
        return strip(src[istart:iend], '\n'), line
    end
    # The function declaration was presumably on a previous line
    lineindex = lastindex(linestarts)
    while !isfuncexpr(ex) && lineindex > 0
        istart = linestarts[lineindex]
        try
            ex, iend = Meta.parse(src, istart)
        catch
        end
        lineindex -= 1
        line -= 1
    end
    return chomp(src[istart:iend-1]), line
end

"""
    ex = definition(Expr, method::Method)
    ex = definition(method::Method)

Return an expression that defines `method`. If the definition can't be found,
returns `nothing`.
"""
function definition(::Type{Expr}, method::Method)
    file = String(method.file)
    def = startswith(file, "REPL[") ? nothing : get(method_info, method.sig, nothing)
    if def === nothing
        f = method_lookup_callback[]
        if f !== nothing
            try
                Base.invokelatest(f, method)
                def = get(method_info, method.sig, nothing)
            catch
            end
        end
    end
    return def === nothing || ismissing(def) ? nothing : copy(def[2])
end

definition(method::Method) = definition(Expr, method)

"""
    info = pkgfiles(name::AbstractString)
    info = pkgfiles(name::AbstractString, uuid::UUID)

Return a [`CodeTracking.PkgFiles`](@ref) structure with information about the files that
define the package specified by `name` and `uuid`.
Returns `nothing` if this package has not been loaded.
"""
pkgfiles(name::AbstractString, uuid::UUID) = pkgfiles(PkgId(uuid, name))
function pkgfiles(name::AbstractString)
    project = Base.active_project()
    # The value returned by Base.project_deps_get depends on the Julia version
    id = Base.project_deps_get(project, name)
    (id == false || id === nothing) && error("no package ", name, " recognized")
    return isa(id, PkgId) ? pkgfiles(id) : pkgfiles(name, id)
end
pkgfiles(id::PkgId) = get(_pkgfiles, id, nothing)

"""
    info = pkgfiles(mod::Module)

Return a [`CodeTracking.PkgFiles`](@ref) structure with information about the files that
were loaded to define the package that defined `mod`.
"""
pkgfiles(mod::Module) = pkgfiles(PkgId(mod))

if ccall(:jl_generating_output, Cint, ()) == 1
    precompile(Tuple{typeof(setindex!), Dict{PkgId,PkgFiles}, PkgFiles, PkgId})
end

end # module


text:
module Custom

using MLJ, MLJScikitLearnInterface, DataFrames

export MLJ, MLJScikitLearnInterface, transform

function transform(data, model)
    select!(data, Not(:Species))
end

end

max_stars_repo_name: amperie/user-models

max_issues_repo_name: amperie/user-models

max_forks_repo_name: amperie/user-models


text:
# Modelo de flash provis򱨯rio, pois o PP tem cૣlculo de flash mas
# esta fun袯⭠o n⭠o est� disponibilizada no plugin (esta fun袯⭠o seria mais
# eficiente nos cૣlculos, n⭠o teria cૣlculos repetitivos)
type ETA_CICLO
	ETA_CICLO()=begin
		new(
			Eficiencia(),
			Pot_sinal(Dict{Symbol,Any}(
				:Symbol=>"_{inB1}",
				:PosX=>0.2,
				:PosY=>0
			)),
			Pot_sinal(Dict{Symbol,Any}(
				:Symbol=>"_{inB2}",
				:PosX=>0.4,
				:PosY=>0
			)),
			Pot_sinal(Dict{Symbol,Any}(
				:Symbol=>"_{inB3}",
				:PosX=>0.6,
				:PosY=>0
			)),
			Pot_sinal(Dict{Symbol,Any}(
				:Symbol=>"_{inB4}",
				:PosX=>0.8,
				:PosY=>0
			)),
			Pot_sinal(Dict{Symbol,Any}(
				:Symbol=>"_{inB5}",
				:PosX=>0.25,
				:PosY=>1
			)),
			Pot_sinal(Dict{Symbol,Any}(
				:Symbol=>"_{inGE}",
				:PosX=>0.5,
				:PosY=>1
			)),
			Pot_sinal(Dict{Symbol,Any}(
				:Symbol=>"_{inGV}",
				:PosX=>0.75,
				:PosY=>1
			)),
			[
				:(EF_CIC * POT_GV = POT_GE - POT_BMB1 - POT_BMB2 - POT_BMB3 - POT_BMB4 - POT_BMB5),
			],
			[
				"Eficiencia do Ciclo",
			],
			[:EF_CIC,:POT_BMB1,:POT_BMB2,:POT_BMB3,:POT_BMB4,:POT_BMB5,:POT_GE,:POT_GV,]
		)
	end
	EF_CIC::Eficiencia
	POT_BMB1::Pot_sinal
	POT_BMB2::Pot_sinal
	POT_BMB3::Pot_sinal
	POT_BMB4::Pot_sinal
	POT_BMB5::Pot_sinal
	POT_GE::Pot_sinal
	POT_GV::Pot_sinal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export ETA_CICLO
function setEquationFlow(in::ETA_CICLO)
	addEquation(1)
end
function atributes(in::ETA_CICLO,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/eficiencia"
	drive!(fields,_)
	return fields
end
ETA_CICLO(_::Dict{Symbol,Any})=begin
	newModel=ETA_CICLO()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
include("../core/eq-core.jl")
include("../core/score-engine.jl")

function test_royal_flush()
    failing = 0
    cards = [
        Card(Ranks.ACE, Suits.CLUBS),
        Card(Ranks.ACE, Suits.DIAMONDS),
        Card(Ranks.ACE, Suits.HEARTS),
        Card(Ranks.ACE, Suits.SPADES),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.TEN, Suits.HEARTS)
    ]
    best_cards = [
        Card(Ranks.ACE, Suits.CLUBS),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.CLUBS),
    ]
    royal_flush = _royal_flush_draw(cards)
    if ! expect(royal_flush, best_cards, desc="ERROR: match royal flush")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.ACE, Suits.CLUBS), Card(Ranks.TEN, Suits.HEARTS)))
    board = Board(Card(Ranks.KING, Suits.CLUBS), Card(Ranks.ACE, Suits.HEARTS), Card(Ranks.JACK, Suits.CLUBS), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.TEN, Suits.CLUBS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.ROYAL, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score royal flush")
        failing += 1
    end
    return failing
end

function test_straight_flush()
    failing = 0
    cards = [
        Card(Ranks.ACE, Suits.SPADES),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.KING, Suits.DIAMONDS),
        Card(Ranks.KING, Suits.HEARTS),
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.NINE, Suits.CLUBS)
    ]
    best_cards = [
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.NINE, Suits.CLUBS)
    ]
    straight_flush = _straight_flush_draw(cards)
    if ! expect(straight_flush, best_cards, desc="ERROR: match straight flush")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.ACE, Suits.SPADES), Card(Ranks.TEN, Suits.CLUBS)))
    board = Board(Card(Ranks.KING, Suits.CLUBS), Card(Ranks.NINE, Suits.CLUBS), Card(Ranks.JACK, Suits.CLUBS), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.TEN, Suits.HEARTS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.SFLUSH, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score straight flush")
        failing += 1
    end
    return failing
end

function test_four_of_a_kind()
    failing = 0
    cards = [
        Card(Ranks.ACE, Suits.SPADES),
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.KING, Suits.HEARTS),
        Card(Ranks.KING, Suits.DIAMONDS),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.NINE, Suits.SPADES),
        Card(Ranks.TWO, Suits.CLUBS)
    ]
    best_cards = [
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.KING, Suits.HEARTS),
        Card(Ranks.KING, Suits.DIAMONDS),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.ACE, Suits.SPADES)
    ]
    quad = _quad_draw(cards)
    if ! expect(quad, best_cards, desc="ERROR: match four-of-a-kind")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.KING, Suits.SPADES), Card(Ranks.KING, Suits.HEARTS)))
    board = Board(Card(Ranks.KING, Suits.CLUBS), Card(Ranks.ACE, Suits.SPADES), Card(Ranks.JACK, Suits.CLUBS), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.KING, Suits.DIAMONDS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.QUAD, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score four-of-a-kind")
        failing += 1
    end

    return failing
end

function test_full_house()
    failing = 0
    cards = [
        Card(Ranks.ACE, Suits.SPADES),
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.KING, Suits.DIAMONDS),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.TWO, Suits.CLUBS),
        Card(Ranks.NINE, Suits.HEARTS)
    ]
    best_cards = [
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.KING, Suits.DIAMONDS),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.TEN, Suits.CLUBS)
    ]
    full_house = _full_draw(cards)
    if ! expect(full_house, best_cards, desc="ERROR: match full house")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.KING, Suits.SPADES), Card(Ranks.KING, Suits.CLUBS)))
    board = Board(Card(Ranks.KING, Suits.DIAMONDS), Card(Ranks.NINE, Suits.CLUBS), Card(Ranks.TEN, Suits.CLUBS), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.TEN, Suits.HEARTS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.FULL, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score full house")
        failing += 1
    end
    return failing
end

function test_flush()
    failing = 0
    cards = [
        Card(Ranks.ACE, Suits.SPADES),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.TEN, Suits.SPADES),
        Card(Ranks.TWO, Suits.CLUBS),
        Card(Ranks.NINE, Suits.HEARTS)
    ]
    best_cards = [
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.TWO, Suits.CLUBS)
    ]
    flush = _flush_draw(cards)
    if ! expect(flush, best_cards, desc="ERROR: match flush")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.TWO, Suits.CLUBS), Card(Ranks.KING, Suits.CLUBS)))
    board = Board(Card(Ranks.KING, Suits.DIAMONDS), Card(Ranks.NINE, Suits.HEARTS), Card(Ranks.TEN, Suits.CLUBS), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.JACK, Suits.CLUBS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.FLUSH, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score flush")
        failing += 1
    end
    return failing
end

function test_straight()
    failing = 0
    cards = [
        Card(Ranks.ACE, Suits.SPADES),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.SPADES),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.NINE, Suits.HEARTS)
    ]
    best_cards = [
        Card(Ranks.ACE, Suits.SPADES),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.SPADES),
    ]
    straight = _straight_draw(cards)
    if ! expect(straight, best_cards, desc="ERROR: match straight")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.ACE, Suits.SPADES), Card(Ranks.KING, Suits.CLUBS)))
    board = Board(Card(Ranks.TEN, Suits.SPADES), Card(Ranks.NINE, Suits.HEARTS), Card(Ranks.TEN, Suits.SPADES), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.JACK, Suits.CLUBS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.STRAIGHT, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score straight")
        failing += 1
    end
    return failing
end

function test_three_of_a_kind()
    failing = 0
    cards = [
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.SPADES),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.TWO, Suits.HEARTS)
    ]
    best_cards = [
        Card(Ranks.TEN, Suits.SPADES),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.TEN, Suits.CLUBS),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
    ]
    trip = _trip_draw(cards)
    if ! expect(trip, best_cards, desc="ERROR: match three-of-a-kind")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.TEN, Suits.CLUBS), Card(Ranks.KING, Suits.CLUBS)))
    board = Board(Card(Ranks.TEN, Suits.SPADES), Card(Ranks.TWO, Suits.HEARTS), Card(Ranks.TEN, Suits.HEARTS), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.JACK, Suits.CLUBS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.TRIP, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score three-of-a-kind")
        failing += 1
    end
    return failing
end

function test_two_pair()
    failing = 0
    cards = [
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.SPADES),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.TWO, Suits.HEARTS)
    ]
    best_cards = [
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.KING, Suits.CLUBS),
        Card(Ranks.TEN, Suits.SPADES),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.QUEEN, Suits.CLUBS),
    ]
    dubs = _dubs_draw(cards)
    if ! expect(dubs, best_cards, desc="ERROR: match two-pair")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.KING, Suits.SPADES), Card(Ranks.KING, Suits.CLUBS)))
    board = Board(Card(Ranks.TEN, Suits.SPADES), Card(Ranks.TWO, Suits.HEARTS), Card(Ranks.TEN, Suits.HEARTS), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.JACK, Suits.CLUBS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.DUBS, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score two-pair")
        failing += 1
    end
    return failing
end

function test_pair()
    failing = 0
    cards = [
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.SPADES),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.FOUR, Suits.DIAMONDS),
        Card(Ranks.TWO, Suits.HEARTS)
    ]
    best_cards = [
        Card(Ranks.TEN, Suits.SPADES),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS)
    ]
    pair = _pair_draw(cards)
    if ! expect(pair, best_cards, desc="ERROR: match pair")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.KING, Suits.SPADES), Card(Ranks.FOUR, Suits.DIAMONDS)))
    board = Board(Card(Ranks.TEN, Suits.SPADES), Card(Ranks.TWO, Suits.HEARTS), Card(Ranks.TEN, Suits.HEARTS), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.JACK, Suits.CLUBS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.PAIR, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score pair")
        failing += 1
    end
    return failing
end

function test_high_card()
    failing = 0
    cards = [
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.SEVEN, Suits.HEARTS),
        Card(Ranks.FOUR, Suits.DIAMONDS),
        Card(Ranks.TWO, Suits.HEARTS)
    ]
    best_cards = [
        Card(Ranks.KING, Suits.SPADES),
        Card(Ranks.QUEEN, Suits.CLUBS),
        Card(Ranks.JACK, Suits.CLUBS),
        Card(Ranks.TEN, Suits.HEARTS),
        Card(Ranks.SEVEN, Suits.HEARTS)
    ]
    high = _high_draw(cards)
    if ! expect(high, best_cards, desc="ERROR: match high card")
        failing += 1
    end
    player = Player(0, Hole(0, Card(Ranks.KING, Suits.SPADES), Card(Ranks.FOUR, Suits.DIAMONDS)))
    board = Board(Card(Ranks.SEVEN, Suits.HEARTS), Card(Ranks.TWO, Suits.HEARTS), Card(Ranks.TEN, Suits.HEARTS), Card(Ranks.QUEEN, Suits.CLUBS), Card(Ranks.JACK, Suits.CLUBS))
    hand = score(player, board)
    best_hand = Hand(0, Rankings.HIGH, best_cards...)
    if ! expect(hand, best_hand, desc="ERROR: score high card")
        failing += 1
    end

    return failing
end

function test_compare_hands()
    failing = 0
    board = Board(Card(Ranks.KING, Suits.HEARTS), Card(Ranks.QUEEN, Suits.HEARTS), Card(Ranks.JACK, Suits.HEARTS), Card(Ranks.TEN, Suits.HEARTS), Card(Ranks.TEN, Suits.CLUBS))
    p1 = Player(1, Hole(1, Card(Ranks.ACE, Suits.HEARTS), Card(Ranks.TWO, Suits.CLUBS)))
    p2 = Player(2, Hole(2, Card(Ranks.NINE, Suits.HEARTS), Card(Ranks.TWO, Suits.DIAMONDS)))
    p3 = Player(3, Hole(3, Card(Ranks.TEN, Suits.DIAMONDS), Card(Ranks.TEN, Suits.SPADES)))
    p4 = Player(4, Hole(4, Card(Ranks.KING, Suits.CLUBS), Card(Ranks.KING, Suits.SPADES)))
    p5 = Player(5, Hole(5, Card(Ranks.TWO, Suits.HEARTS), Card(Ranks.THREE, Suits.CLUBS)))
    p6 = Player(6, Hole(6, Card(Ranks.ACE, Suits.SPADES), Card(Ranks.THREE, Suits.SPADES)))
    p7 = Player(7, Hole(7, Card(Ranks.TEN, Suits.SPADES), Card(Ranks.FOUR, Suits.SPADES)))
    p8 = Player(8, Hole(8, Card(Ranks.JACK, Suits.SPADES), Card(Ranks.FOUR, Suits.CLUBS)))
    p9 = Player(9, Hole(9, Card(Ranks.EIGHT, Suits.CLUBS), Card(Ranks.SEVEN, Suits.CLUBS)))
    p10 = Player(10, Hole(10, Card(Ranks.FIVE, Suits.CLUBS), Card(Ranks.SIX, Suits.SPADES)))
    p1_hand = score(p1, board)
    p2_hand = score(p2, board)
    p3_hand = score(p3, board)
    p4_hand = score(p4, board)
    p5_hand = score(p5, board)
    p6_hand = score(p6, board)
    p7_hand = score(p7, board)
    p8_hand = score(p8, board)
    p9_hand = score(p9, board)
    p10_hand = score(p10, board)
    begin
        best_hands = compare_hands(p1_hand, p2_hand, p3_hand, p4_hand, p5_hand, p6_hand, p7_hand, p8_hand, p9_hand, p10_hand)
        if ! expect(best_hands, [p1_hand], desc="ERROR: royal flush is best hand")
            failing += 1
        end
    end
    begin
        best_hands = compare_hands(p2_hand, p3_hand, p4_hand, p5_hand, p6_hand, p7_hand, p8_hand, p9_hand, p10_hand)
        if ! expect(best_hands, [p2_hand], desc="ERROR: straight flush is best hand")
            failing += 1
        end
    end
    begin
        best_hands = compare_hands(p3_hand, p4_hand, p5_hand, p6_hand, p7_hand, p8_hand, p9_hand, p10_hand)
        if ! expect(best_hands, [p3_hand], desc="ERROR: four-of-a-kind is best hand")
            failing += 1
        end
    end
    begin
        best_hands = compare_hands(p4_hand, p5_hand, p6_hand, p7_hand, p8_hand, p9_hand, p10_hand)
        if ! expect(best_hands, [p4_hand], desc="ERROR: full-house is best hand")
            failing += 1
        end
    end
    begin
        best_hands = compare_hands(p5_hand, p6_hand, p7_hand, p8_hand, p9_hand, p10_hand)
        if ! expect(best_hands, [p5_hand], desc="ERROR: flush is best hand")
            failing += 1
        end
    end
    begin
        best_hands = compare_hands(p6_hand, p7_hand, p8_hand, p9_hand, p10_hand)
        if ! expect(best_hands, [p6_hand], desc="ERROR: straight is best hand")
            failing += 1
        end
    end
    begin
        best_hands = compare_hands(p7_hand, p8_hand, p9_hand, p10_hand)
        if ! expect(best_hands, [p7_hand], desc="ERROR: three-of-a-kind is best hand")
            failing += 1
        end
    end
    begin
        best_hands = compare_hands(p8_hand, p9_hand, p10_hand)
        if ! expect(best_hands, [p8_hand], desc="ERROR: two-pair is best hand")
            failing += 1
        end
    end
    begin
        best_hands = compare_hands(p9_hand, p10_hand)
        if ! expect(best_hands, [p9_hand, p10_hand], desc="ERROR: tie hands")
            failing += 1
        end
    end
    return failing
end

function expect(cards::Array{Card}, exp::Array{Card}; desc="")
    if length(cards) != length(exp)
        println(desc)
        return false
    end
    for i in 1:length(cards)
        if ! match(cards[i], exp[i])
            println(desc)
            return false
        end
    end
    return true
end

function expect(cards::Nothing, exp::Array{Card}; desc="")
    println(desc)
    return false
end

function expect(hand::Hand, exp::Hand; desc="")
    if hand.playerID != exp.playerID
        println(desc)
        return false
    end
    if hand.ranking != exp.ranking
        println(desc)
        return false
    end
    if ! match(hand.card1, exp.card1)
        println(desc)
        return false
    end
    if ! match(hand.card2, exp.card2)
        println(desc)
        return false
    end
    if ! match(hand.card3, exp.card3)
        println(desc)
        return false
    end
    if ! match(hand.card4, exp.card4)
        println(desc)
        return false
    end
    if ! match(hand.card5, exp.card5)
        println(desc)
        return false
    end
    return true
end

function expect(hand::Nothing, exp::Hand; desc="")
    println(desc)
    return false
end

function expect(best_hands::Array{Hand}, exp::Array{Hand}; desc="")
    if length(best_hands) != length(exp)
        println(desc)
        return false
    end
    for i in 1:length(best_hands)
        if best_hands[i].playerID != exp[i].playerID
            println(desc)
            return false
        end
    end
    return true
end

function match(c1::Card, c2::Card)
    return c1.rank == c2.rank && c1.suit == c2.suit
end

function run_tests()
    failing = 0
    failing += test_royal_flush()
    failing += test_straight_flush()
    failing += test_four_of_a_kind()
    failing += test_full_house()
    failing += test_flush()
    failing += test_straight()
    failing += test_three_of_a_kind()
    failing += test_two_pair()
    failing += test_pair()
    failing += test_high_card()
    failing += test_compare_hands()
    println("Done!")
    if failing == 0
        println("All tests passing!")
    else
        println("$(failing) test faililng!")
    end
end

run_tests()

max_stars_repo_name: ttowncompiled/EQPoker

max_issues_repo_name: ttowncompiled/EQPoker

max_forks_repo_name: ttowncompiled/EQPoker


text:
using MLStyle

@testset "adt-type" begin

@data List{T} begin
    Nil{T}
    Cons{T}(head :: T, tail :: List{T})
end

lst = Cons(1, Nil{Int}())

let test =
    @match lst begin
    	   Cons(1, ::Cons) => nothing
           Cons(head, ::Nil{Int}) => head
	end
    @test test === 1
end

end


text:
function evolve()
sinfo = Params()

#== start template ==#
#fundamental constants/units
  ħ = 1.0545718e-34
  kB = 1.38064852e-23
  amu = 1.660339040e-27
  Bohr = 5.29e-11
#Rb87 mass and scattering length
  m = 86.909180527*amu
  as = 100*Bohr
#trap frequencies
  ωx = 2π
  ωy = 4ωx
  ωz = 0.
#choice of time, length, energy units
  t0 = 1.0/ωy
  x0 = sqrt(ħ*t0/m)
  E0 = ħ/t0
#interactions
  #g  = (4*pi*ħ^2*as/m)*x0^3/E0 #dimensionless 3D
  g = 0.1 #test 2D
#damping parameters (dimensionless)
  γ  = 0.1
  ℳ  = 0.0
#chemical potential (dimensionless)
  μ  = 12.0
#time evolution parameters
  ti = 0.0
  tf = 1.0/γ  #evolve for 2 damping times
  Nt = 50
  t = collect(linspace(ti,tf,Nt))
  dt = 0.01π/μ #integrate step size [ - should have dt ≪ 2π/μ]
#== end template ==#

  @pack sinfo = ωx,ωy,ωz,γ,ℳ,g,t0,E0,x0,μ,ti,tf,Nt,t,dt

  #Initialize CField
  basis = "Hermite"
  ecut = 30 #units of ħ*ωy
  Ω = [ωx; ωy]*t0
  cinfo = makecinfo(ecut,Ω,basis)
  tinfo = maketinfo(cinfo,4)
  @unpack en, P, M = cinfo; Mx, My = M;
  #x,wx,Tx,y,wy,Ty = makealltrans(M,Ω,n=4,basis=basis)
  #W = wx.*wy'
#test transform
  c0   = randn(Mx,My)+im*randn(Mx,My);
  c0   = P.*c0;   #Project
  ψ0   = c2x(c0,tinfo)
  ψ    = simlar(ψ0) #a field to write to in place if needed

c0    = randn(Mx,My) + im*randn(Mx,My); #create new random initial state
tspan = (t[1],t[end])
prob = ODEProblem(Lgp!,c0,tspan)
alg = DP5()
#abstol = 1e-3!
#reltol = 1e-3
println("Started evolution ...")
@time sol = solve(prob,alg,dt=dt,saveat=t);
println("... Finished.")
return sinfo,cinfo,sol
end


text:
module Walkforward

export xwalkforward

include("xwalkforward.jl")

end


text:
# Let's test latex commands to death... especially in light of #444

@testset "lx++1" begin
    # 444
    s = "\\newcommand{\\note}[1]{#1} \\note{A `B` C} D" |> fd2html_td
    @test isapproxstr(s, """
        A B C D
        """)
    s = raw"""
        \newcommand{\note}[1]{@@note #1 @@}
        \note{A}
        \note{A `B` C}
        \note{A @@cc B @@ D}
        \note{A @@cc B `D` E @@ F}
        """ |> fd2html_td
    @test isapproxstr(s, """

A

A B C

A
B
 D

A
B D E
 F

        """)
end


text:
type LSTM{T} <: Layer{T}
    name::Symbol
    W::Matrix{T} # 4H×(D+H)
        # consists of Wx and Wh
        #  Wx: Input-to-hidden weights, of shape (4H, D)
        #  Wh: Hidden-to-hidden weights, of shape (4H, H)
    b::Vector{T} # 4H
    forget_bias::T

    input::Vector{T} # input, [D]
    state::Vector{T} # state (2H), contains [cell state, hidden state]

    # preallocated memory
    c::Vector{T} # [H] cell state
    m::Vector{T} # [H] member state
    cell_inputs::Vector{T}  # [D+H]
    lstm_mult_res::Vector{T} # [4H]
    i::Vector{T} # [H] input gate
    j::Vector{T} # [H] new input
    f::Vector{T} # [H] forget gate
    o::Vector{T} # [H] output gate
end
name(a::LSTM) = a.name
output{T}(a::LSTM{T}) = a.m
zero!{T}(a::LSTM{T}) = fill!(a.state, zero(T))
function Base.push!{T}(net::ForwardNet{T}, ::Type{LSTM},
    name::Symbol,
    parent_index::NameOrIndex,
    H::Int; # hidden layer size
    forget_bias::T = one(T),
    )

    input = output(net[parent_index])::Vector{T}
    D = length(input)

    W = Array{T}(4H, D+H)
    b = Array{T}(4H)
    state = Array{T}(2H)

    c = Array{T}(H)
    m = Array{T}(H)
    cell_inputs = Array{T}(D+H)
    lstm_mult_res = Array{T}(4H)
    i = Array{T}(H)
    j = Array{T}(H)
    f = Array{T}(H)
    o = Array{T}(H)

    node = LSTM(name, W, b, forget_bias, input, state, c, m, cell_inputs, lstm_mult_res, i, j, f, o)
    push!(net, node, parent_index)
end
function forward!{T}(a::LSTM{T})

    H = length(a.c)
    D = length(a.input)

    copy!(a.c, 1, a.state, 1, H)
    copy!(a.m, 1, a.state, H+1, H)

    copy!(a.cell_inputs, 1, a.input, 1, D)
    copy!(a.cell_inputs, D+1, a.m, 1, H)

    # lstm_matrix = a.W * a.cell_inputs + a.b
    copy!(a.lstm_mult_res, a.b) # y ← b
    Base.LinAlg.BLAS.gemv!('N', one(T), a.W, a.cell_inputs, one(T), a.lstm_mult_res) # y ← W*x + y

    copy!(a.i, 1, a.lstm_mult_res,    1, H)
    copy!(a.j, 1, a.lstm_mult_res,  H+1, H)
    copy!(a.f, 1, a.lstm_mult_res, 2H+1, H)
    copy!(a.o, 1, a.lstm_mult_res, 3H+1, H)

    # compute c_next and m_next
    for k in 1 : H
        a.f[k] += a.forget_bias
        c = sigmoid(a.f[k])*a.c[k] + sigmoid(a.i[k])*tanh(a.j[k]) # c_next
        a.state[k] = a.c[k] = c
        a.state[k+H] = a.m[k] = sigmoid(a.o[k])*tanh(c) # m_next
    end

    a
end
function restore!{T}(a::LSTM{T}, filename_W::String, filename_b::String)
    H = length(a.c)
    D = length(a.input)

    vec = open(io->read_binary_vec(io, T), filename_W)
    a.W[:] = convert_to_column_major_array(vec, (D+H,4H))'

    copy!(a.b, open(read_binary_vec, filename_b))

    a
end

#=
u = f(W_xu * x + W_hu*h_{t-1} + b_u)
r = f(W_xr * x + W_hr*h_{t-1} + b_r)
c = tanh(W_xc * x + (W_hc*h_{t-1}) .* r + b_c)
h_t = (1-u) .* h_{t-1} + u .* c
=#
type GRU{T} <: Layer{T}
    name::Symbol
    W_x::Matrix{T} # 3H×D
    W_h::Matrix{T} # 3H×H
    b::Vector{T} # 3H

    input::Vector{T} # input, [D]

    # preallocated memory
    xb_ruc::Vector{T}  # [3H]
    xb_r::Vector{T} # [H]
    xb_u::Vector{T} # [H]
    xb_c::Vector{T} # [H]

    h_ruc::Vector{T} # [3H]
    h_r::Vector{T} # [H]
    h_u::Vector{T} # [H]
    h_c::Vector{T} # [H]

    r::Vector{T} # [H] input gate
    u::Vector{T} # [H] new input
    c::Vector{T} # [H] forget gate
    h::Vector{T} # [H] GRU state
    h_prev::Vector{T} # [H] previous state
end
name(a::GRU) = a.name
output{T}(a::GRU{T}) = a.h
zero!{T}(a::GRU{T}) = fill!(a.h_prev, zero(T))
function Base.push!{T}(net::ForwardNet{T}, ::Type{GRU},
    name::Symbol,
    parent_index::NameOrIndex,
    H::Int; # hidden layer size
    )

    input = output(net[parent_index])::Vector{T}
    D = length(input)

    W_x = Array{T}(3H, D)
    W_h = Array{T}(3H, H)
    b = Array{T}(3H)

    xb_ruc = Array{T}(3H)
    xb_r = Array{T}(H)
    xb_u = Array{T}(H)
    xb_c = Array{T}(H)

    h_ruc = Array{T}(3H)
    h_r = Array{T}(H)
    h_u = Array{T}(H)
    h_c = Array{T}(H)

    r = Array{T}(H)
    u = Array{T}(H)
    c = Array{T}(H)
    h = Array{T}(H)
    h_prev = Array{T}(H)

    node = GRU(name, W_x, W_h, b, input, xb_ruc, xb_r, xb_u, xb_c, h_ruc, h_r, h_u, h_c, r, u, c, h, h_prev)
    push!(net, node, parent_index)
end

function forward!{T}(a::GRU{T})

    H = length(a.h_prev)
    D = length(a.input)

    copy!(a.xb_ruc, a.b) # y ← b
    Base.LinAlg.BLAS.gemv!('N', one(T), a.W_x, a.input, one(T), a.xb_ruc) # y ← W*x + y

    copy!(a.h_ruc, Base.LinAlg.BLAS.gemv('N', one(T), a.W_h, a.h_prev)) # y ← W*x

    copy!(a.xb_r, 1, a.xb_ruc,    1, H)
    copy!(a.xb_u, 1, a.xb_ruc,  H+1, H)
    copy!(a.xb_c, 1, a.xb_ruc, 2H+1, H)

    copy!(a.h_r, 1, a.h_ruc,    1, H)
    copy!(a.h_u, 1, a.h_ruc,  H+1, H)
    copy!(a.h_c, 1, a.h_ruc, 2H+1, H)

    # compute u, r, c, and h
    for k in 1 : H
        a.r[k] = sigmoid(a.xb_r[k] + a.h_r[k])
        a.u[k] = sigmoid(a.xb_u[k] + a.h_u[k])
        a.c[k] = tanh(a.xb_c[k] + a.r[k]*a.h_c[k])
        h_old = a.h[k]
        a.h[k] = (1 - a.u[k])*a.h_prev[k] + a.u[k]*a.c[k]
        a.h_prev[k] = h_old
    end

    a
end
function restore!{T}(a::GRU{T}, filename_Wx::String, filename_Wh::String, filename_b::String)
    H = length(a.h_prev)
    D = length(a.input)

    vec = open(io->read_binary_vec(io, T), filename_Wx)
    a.W_x[:] = convert_to_column_major_array(vec, (D,3H))'

    vec = open(io->read_binary_vec(io, T), filename_Wh)
    a.W_h[:] = convert_to_column_major_array(vec, (H,3H))'

    copy!(a.b, open(read_binary_vec, filename_b))

    a
end


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

mutable struct IoK8sApiDiscoveryV1beta1EndpointConditions <: SwaggerModel
    ready::Any # spec type: Union{ Nothing, Bool } # spec name: ready

    function IoK8sApiDiscoveryV1beta1EndpointConditions(;ready=nothing)
        o = new()
        validate_property(IoK8sApiDiscoveryV1beta1EndpointConditions, Symbol("ready"), ready)
        setfield!(o, Symbol("ready"), ready)
        o
    end
end # type IoK8sApiDiscoveryV1beta1EndpointConditions

const _property_map_IoK8sApiDiscoveryV1beta1EndpointConditions = Dict{Symbol,Symbol}(Symbol("ready")=>Symbol("ready"))
const _property_types_IoK8sApiDiscoveryV1beta1EndpointConditions = Dict{Symbol,String}(Symbol("ready")=>"Bool")
Base.propertynames(::Type{ IoK8sApiDiscoveryV1beta1EndpointConditions }) = collect(keys(_property_map_IoK8sApiDiscoveryV1beta1EndpointConditions))
Swagger.property_type(::Type{ IoK8sApiDiscoveryV1beta1EndpointConditions }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiDiscoveryV1beta1EndpointConditions[name]))}
Swagger.field_name(::Type{ IoK8sApiDiscoveryV1beta1EndpointConditions }, property_name::Symbol) =  _property_map_IoK8sApiDiscoveryV1beta1EndpointConditions[property_name]

function check_required(o::IoK8sApiDiscoveryV1beta1EndpointConditions)
    true
end

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