text:
"""
Provides implementation of Adaptive Stress Testing (AST) formulation of MDPs/POMDPs.
"""
module AST

using Random
using RLInterface
using Parameters
using Distributions
using DataStructures
using POMDPPolicies
using POMDPSimulators
using POMDPs

export
    BlackBox,
    GrayBox,

    ASTParams,
    ASTMDP,
    ASTState,
    ASTAction,
    ASTSeedAction,
    ASTSampleAction,
    ASTMetrics,
    actiontype,

    reward,
    action,
    initialstate,
    isterminal,
    hash_uint32,
    reset_metrics!,

    search!,
    playback,
    get_top_path,
    online_path

include("BlackBox.jl")
include("GrayBox.jl")
include("ast_types.jl")
include(joinpath("utils", "rand.jl"))
include(joinpath("utils", "seeding.jl"))
include(joinpath("utils", "hashing.jl"))
include(joinpath("utils", "recording.jl"))

# Julia 1.1 compat
if !@isdefined(hasproperty)
    hasproperty(x, s::Symbol) = s in propertynames(x)
end

# Julia 1.0 compat
if !@isdefined(isnothing)
    isnothing(x) = x === nothing
end

"""
    reward(mdp::ASTMDP, logprob::Real, isevent::Bool, isterminal::Bool, miss_distance::Real)::Float64

Reward function for the AST formulation. Defaults to:

    R_E       if isterminal and isevent (1)
    -d        if isterminal and !isevent (2)
    log(p)    otherwise (3)

    1) Terminates with event, collect reward_bonus (defaults to 0)
    2) Terminates without event, collect negative miss distance
    3) Each non-terminal step, accumulate reward correlated with the transition probability

For epsidic reward problems (i.e. rewards only at the end of an episode), set `mdp.params.episodic_rewards` to get:

    (log(p) - d)*R_E    if isterminal and isevent (1)
    log(p) - d          if isterminal and !isevent (2)
    0                   otherwise (3)

    1) Terminates with event, collect transition probability and miss distance with multiplicative reward bonus
    2) Terminates without event, collect transitions probability and miss distance
    3) Each non-terminal step, no intermediate reward (set `mdp.params.give_intermediate_reward` to use log transition probability)
"""
function POMDPs.reward(mdp::ASTMDP, logprob::Real, isevent::Bool, isterminal::Bool, miss_distance::Real)
    if mdp.params.episodic_rewards
        r = 0
        if isterminal
            r += logprob - miss_distance
            if isevent
                r *= mdp.params.reward_bonus # R_E (multiplicative)
            end
        else
            intermediate_reward = mdp.params.give_intermediate_reward ? logprob : 0
            r += intermediate_reward
        end
    else # Standard AST reward function
        r = logprob
        if isevent
            r += mdp.params.reward_bonus # R_E (additive)
        elseif isterminal
            r += -miss_distance # Only add miss distance cost if is terminal and not an event.
        end
    end

    if !mdp.params.episodic_rewards || (mdp.params.episodic_rewards && isterminal) || (mdp.params.episodic_rewards && mdp.params.give_intermediate_reward)
        record(mdp, prob=exp(logprob), logprob=logprob, miss_distance=miss_distance, reward=r, event=isevent)
    end

    return Float32(r)
end

"""
Initialize AST MDP state. Overridden from `POMDPs.initialstate` interface.
"""
function POMDPs.initialstate(mdp::ASTMDP, rng::AbstractRNG=Random.GLOBAL_RNG) # rng unused.
    mdp.t_index = 1
    BlackBox.initialize!(mdp.sim)

    if !isnothing(mdp.params.reset_seed)
        mdp.params.seed = mdp.params.reset_seed
    end

    if actiontype(mdp) == ASTSeedAction
        s = ASTState(t_index=mdp.t_index, action=ASTSeedAction(mdp.params.seed))
    elseif actiontype(mdp) == ASTSampleAction
        s = ASTState(t_index=mdp.t_index)
    end

    mdp.sim_hash = s.hash
    return s::ASTState
end

"""
Generate next state and reward for AST MDP (handles episodic reward problems). Overridden from `POMDPs.gen` interface.
"""
function POMDPs.gen(mdp::ASTMDP, s::ASTState, actions::Vector{Float32}, rng::AbstractRNG=Random.GLOBAL_RNG)
	a = ASTSampleAction()
	a.sample = GrayBox.pack(mdp.sim, Float64.(actions))
	return POMDPs.gen(mdp, s, a, rng)

end

function POMDPs.gen(mdp::ASTMDP, s::ASTState, a::ASTAction, rng::AbstractRNG=Random.GLOBAL_RNG)
    @assert mdp.sim_hash == s.hash
    mdp.t_index += 1
    isa(a, ASTSeedAction) ? set_global_seed(a) : nothing
    hasproperty(mdp.sim, :actions) ? push!(mdp.sim.actions, a) : nothing

    # Step black-box simulation
    if mdp.params.episodic_rewards
        # Do not evaluate when problem has episodic rewards
        if a isa ASTSeedAction
            if mdp.params.pass_seed_action
                logprob = GrayBox.transition!(mdp.sim, a.seed)
            else
                logprob = GrayBox.transition!(mdp.sim)
            end
        elseif a isa ASTSampleAction
            logprob = GrayBox.transition!(mdp.sim, a.sample)
        end
        isevent::Bool = false
        miss_distance::Float64 = NaN
    else
        if a isa ASTSeedAction
            if mdp.params.pass_seed_action
                (logprob, miss_distance, isevent) = BlackBox.evaluate!(mdp.sim, a.seed)
            else
                (logprob, miss_distance, isevent) = BlackBox.evaluate!(mdp.sim)
            end
        elseif a isa ASTSampleAction
            (logprob, miss_distance, isevent) = BlackBox.evaluate!(mdp.sim, a.sample)
        end
    end

    # Get simlation state
    state = GrayBox.state(mdp.sim)

    # Update MDP state
    sp = ASTState(t_index=mdp.t_index, parent=s, action=a, state=state)
    mdp.sim_hash = sp.hash
    sp.terminal = mdp.params.episodic_rewards ? false : BlackBox.isterminal(mdp.sim) # termination handled by end-of-rollout
    r::Float64 = reward(mdp, logprob, isevent, sp.terminal, miss_distance)
    sp.q_value = r

    # TODO: Optimize (debug?)
    record_trace(mdp, action_q_trace(sp)...)

    return (sp=sp, r=r)
end

"""
Determine if AST MDP is in a terminal state. Overridden from `POMDPs.isterminal` interface.
"""
function POMDPs.isterminal(mdp::ASTMDP, s::ASTState)
    @assert mdp.sim_hash == s.hash
    return BlackBox.isterminal(mdp.sim)
end

"""
AST problems are (generally) undiscounted to treat future reward equally. Overridden from `POMDPs.discount` interface.
"""
POMDPs.discount(mdp::ASTMDP) = mdp.params.discount

"""
Randomly select next action, independent of the state.
"""
random_action(mdp::ASTMDP, ::ASTState) = random_action(mdp)
function random_action(mdp::ASTMDP{ASTSeedAction})
    # Note, cannot use `rand` as that's controlled by the global seed.
    set_seed!(mdp)
    return ASTSeedAction(mdp.current_seed)
end

function random_action(mdp::ASTMDP{ASTSampleAction})
    environment::GrayBox.Environment = GrayBox.environment(mdp.sim)
    sample::GrayBox.EnvironmentSample = rand(environment)
    return GrayBox.unpack(mdp.sim, sample)
    #return ASTSampleAction(sample)
end

"""
Randomly select next action, independent of the state. Overridden from `POMDPs.action` interface.
"""
POMDPs.action(policy::RandomPolicy, s::ASTState) = random_action(policy.problem)
"""
Randomly select next actions, independent of the state. Overridden from `POMDPs.actions` interface.
"""
POMDPs.actions(mdp::ASTMDP) = [random_action(mdp)]

import Base.rand
rand(rng::AbstractRNG, s::ASTState) = s

"""
Used by the RLInterface to interact with deep RL solvers.
"""
function POMDPs.convert_s(::Type{Vector{Float32}}, s::ASTState, mdp::ASTMDP)
	if s.state == nothing
		s.state = GrayBox.state(mdp.sim)
	end
	return s.state
end

"""
Reset AST simulation to a given state; used by the MCTS DPWSolver as the `reset_callback` function.
"""
function go_to_state(mdp::ASTMDP, target_state::ASTState)
    s = initialstate(mdp, Random.GLOBAL_RNG)
    BlackBox.initialize!(mdp.sim)
    actions = action_trace(target_state)
    R = 0.0
    for a in actions
        s, r = @gen(:sp, :r)(mdp, s, a, Random.GLOBAL_RNG)
        R += r
    end
    @assert s == target_state

    # record_trace(mdp, actions, R)

    return (R, actions)
end

"""
Record the best paths from termination leaf node.
"""
function record_trace(mdp::ASTMDP, actions::Vector{ASTAction}, reward::Float64)
    if mdp.params.top_k > 0 && BlackBox.isterminal(mdp.sim)
        if !haskey(mdp.top_paths, actions)
            enqueue!(mdp.top_paths, actions, reward)
            while length(mdp.top_paths) > mdp.params.top_k
                dequeue!(mdp.top_paths)
            end
        end
    end
end

"""
Get k-th top path from the recorded `top_paths`.
"""
get_top_path(mdp::ASTMDP, k=min(mdp.params.top_k, length(mdp.top_paths))) = collect(keys(mdp.top_paths))[k]

"""
Rollout simulation for MCTS; used by the MCTS DPWSolver as the `estimate_value` function.
Custom rollout records action trace once the depth has been reached.
"""
function rollout(mdp::ASTMDP, s::ASTState, d::Int64)
    if d == 0 || isterminal(mdp, s)
        go_to_state(mdp, s) # Records trace through this call
        return 0.0
    else
        a::ASTAction = random_action(mdp)

        (sp, r) = @gen(:sp, :r)(mdp, s, a, Random.GLOBAL_RNG)
        q_value = r + discount(mdp)*rollout(mdp, sp, d-1)

        return q_value
    end
end

"""
Rollout to only execute SUT at end (`p` accounts for probabilities generated outside the rollout)

User defined:

`feed_gen` Function to feed best action, replaces call to `gen` when feeding

`feed_type` Indicate when to feed best action. Either at the start of the rollout `:start`, or mid-rollout `:mid`

`best_callback` Callback function to record best miss distance or reward for later feeding during rollout
"""
function rollout_end(mdp::ASTMDP, s::ASTState, d::Int64; max_depth=-1, feed_gen=missing, feed_type::Symbol=:none, best_callback::Function=(sm,r)->sm)
    sim = mdp.sim
    a::ASTAction = random_action(mdp)

    if d == 0 || isterminal(mdp, s)
        # End-of-rollout evaluations.

        isa(a, ASTSeedAction) ? set_global_seed(a) : nothing
        hasproperty(sim, :actions) ? push!(sim.actions, a) : nothing

        # Step black-box simulation
        if a isa ASTSeedAction
            if mdp.params.pass_seed_action
                (prob::Float64, miss_distance::Float64, isevent::Bool) = BlackBox.evaluate!(sim, a.seed)
            else
                (prob, miss_distance, isevent) = BlackBox.evaluate!(sim)
            end
        elseif a isa ASTSampleAction
            (prob, miss_distance, isevent) = BlackBox.evaluate!(sim, a.sample)
        end

        # Update state
        sp = ASTState(t_index=mdp.t_index, parent=s, action=a)
        mdp.sim_hash = sp.hash
        sp.terminal = BlackBox.isterminal(sim)
        r::Float64 = reward(mdp, prob, isevent, sp.terminal, miss_distance)
        sp.q_value = r

        best_callback(sim, miss_distance) # could use `r` instead

        return r
    else
        # Start of rollout best-action feeding OR mid-rollout best-action feeding.
        feed_available::Bool = !ismissing(feed_gen) && feed_type != :none
        start_of_rollout_feed::Bool = feed_type == :start && d == max_depth-1
        mid_rollout_feed::Bool = feed_type == :mid && d == div(max_depth,2)

        if feed_available && (start_of_rollout_feed || mid_rollout_feed)
            (sp, r) = feed_gen(mdp, s, a, Random.GLOBAL_RNG)
        else
            (sp, r) = @gen(:sp, :r)(mdp, s, a, Random.GLOBAL_RNG)
        end

        # Note, pass all keywords.
        q_value = r + discount(mdp)*rollout_end(mdp, sp, d-1; max_depth=max_depth, feed_gen=feed_gen, feed_type=feed_type, best_callback=best_callback)

        return q_value
    end
end

"""
Play back a given action trace from the `initialstate` of the MDP.
"""
playback(mdp::ASTMDP, actions::Nothing, func=nothing; kwargs...) = @warn("Action trace is `nothing`, please set mdp.params.top_k > 0.")
playback(planner, actions::Vector{ASTAction}, func=nothing; kwargs...) = playback(planner.mdp, actions, func; kwargs...)
function playback(mdp::ASTMDP, actions::Vector{ASTAction}, func=nothing; verbose=true, return_trace::Bool=false)
    rng = Random.GLOBAL_RNG # Not used.
    s = initialstate(mdp, rng)
    BlackBox.initialize!(mdp.sim)
    trace = []
    function trace_and_show(func)
        single_step = func(mdp.sim)
        push!(trace, single_step)
        println(single_step)
    end
    display_trace::Bool = verbose && !isnothing(func)
    display_trace ? trace_and_show(func) : nothing

    for a in actions
        isnothing(func) && verbose ? println(string(a)) : nothing
        (sp, r) = @gen(:sp, :r)(mdp, s, a, rng)
        s = sp
        display_trace ? trace_and_show(func) : nothing
    end

    if return_trace
        return trace::Vector # Returns trace output by `func`
    else
        return s::ASTState # Returns final state
    end
end

"""
Follow MCTS optimal path online calling `action` after each selected state.
"""
function online_path(mdp::MDP, planner::Policy, printstep=(sim, a)->println("Action: ", string(a));
                     verbose::Bool=false)
    s = initialstate(mdp, Random.GLOBAL_RNG)
    BlackBox.initialize!(mdp.sim)

    # First step
    a = action(planner, s)
    actions = ASTAction[a]
    verbose ? printstep(mdp.sim, a) : nothing
    (s, r) = @gen(:sp, :r)(mdp, s, a, Random.GLOBAL_RNG)

    while !BlackBox.isterminal(mdp.sim)
        a = action(planner, s)
        push!(actions, a)
        verbose ? printstep(mdp.sim, a) : nothing
        (s, r) = @gen(:sp, :r)(mdp, s, a, Random.GLOBAL_RNG)
    end

    # Last step: last action is NULL.
    ActionType::Type = actiontype(mdp)
    verbose ? printstep(mdp.sim, ActionType()) : nothing

    if BlackBox.isevent(mdp.sim)
        @info "Event found!"
    else
        @info "Hit terminal state."
    end

    return actions
end

"""
Trace up the tree to get all ancestor states.
"""
function state_trace(s::ASTState)
    states::Vector{ASTState} = [s]
    while !isnothing(s.parent)
        prepend!(states, [s.parent])
        s = s.parent
    end
    return states::Vector{ASTState}
end

"""
Trace up the tree to get all ancestor actions.
"""
function action_trace(s::ASTState)
    actions::Vector{ASTAction} = ASTAction[]

    # Trace up the tree
    while !isnothing(s.parent)
        prepend!(actions, [s.action])
        s = s.parent
    end

    return actions
end

"""
Trace up the tree and accumulate Q-values.
"""
function q_trace(s::ASTState)
    q_value::Float64 = s.q_value
    while !isnothing(s.parent)
        q_value += s.parent.q_value
        s = s.parent
    end
    return q_value::Float64
end

"""
Trace up the tree to get all ancestor actions and summed Q-values.
"""
function action_q_trace(s::ASTState)
    actions::Vector{ASTAction} = ASTAction[]
    q_value::Float64 = s.q_value

    # Trace up the tree
    while !isnothing(s.parent)
        prepend!(actions, [s.action])
        q_value += s.parent.q_value
        s = s.parent
    end

    return actions::Vector, q_value::Float64
end

"""
Clear data stored in `mdp.metrics`.
"""
function reset_metrics!(mdp::ASTMDP)
    mdp.metrics = ASTMetrics()
end

"""
    search!(planner)

Search for failures given a `planner`. Implemented by each solver.
"""
function search! end

end  # module AST


text:
########
# Math #
########

@inline invlogit(x::T) where T<:Real = one(T) / (one(T) + exp(-x))
@inline logit(x::T) where T<:Real = log(x / (one(T) - x))

# More stable, faster version of rand(Categorical)
function randcat(p::Vector{Float64})
  # if(any(p .< 0)) error("Negative probabilities not allowed"); end
  r, s = rand(), one(Int)
  for j = 1:length(p)
    r -= p[j]
    if(r <= 0.0) s = j; break; end
  end

  s
end

struct NotImplementedException <: Exception end

# Numerically stable sum of values represented in log domain.
logsum(xs::Vector{T}) where T<:Real = begin
  largest = maximum(xs)
  ys = map(x -> exp.(x - largest), xs)

  log(sum(ys)) + largest
end

# KL-divergence
kl(p::Normal, q::Normal) = (log(q.σ / p.σ) + (p.σ^2 + (p.μ - q.μ)^2) / (2 * q.σ^2) - 0.5)

align_internal!(x,n) = begin
  m = length(x)
  resize!(x, n)
  x[m+1:end] .= zero(eltype(x))
end

align(x,y) = begin
  if length(x) < length(y)
    align_internal!(x, length(y))
  elseif length(x) > length(y)
    align_internal!(y, length(x))
  end

  (x,y)
end

#######
# I/O #
#######

macro sym_str(var)
  var_str = string(var)
  :(Symbol($var_str))
end

##########
# Helper #
##########

auto_tune_chunk_size!(mf::Function, rep_num=10) = begin
  dim = length(mf().vals)
  if dim > 8
    n = 1
    sz_cand = Int[]
    while dim / n > 8
      push!(sz_cand, ceil(Int, dim / n))
      n += 1
    end
    filter!(sz -> 8 < sz <= 50, sz_cand)
    sz_num = length(sz_cand)
    prof_log = Vector{Float64}(sz_num)
    for i = 1:sz_num
      println("[Turing] profiling chunk size = $(sz_cand[i])")
      setchunksize(sz_cand[i])
      prof_log[i] = @elapsed for _ = 1:rep_num mf() end
    end
    minval, minidx = findmin(prof_log)
    println("[Turing] final chunk size chosen = $(sz_cand[minidx])")
    setchunksize(sz_cand[minidx])
  else
    setchunksize(8)
  end
end

ext: jl


text:
# Algorithm 5.2
# Fourth order Runge-Kutta

function rk4(f, a, b, N, α)
	result = Tuple{Float64, Float64}[]
	k = Array{Float64}(undef, 4)
	h = (b - a)/N
	t = α
	ω = α
	push!(result, (t, ω))
	for i ∈ 1:N
		k[1] = h*f(t, ω)
		k[2] = h*f(t + h/2, ω + k[1]/2)
		k[3] = h*f(t + h/2, ω + k[2]/2)
		k[4] = h*f(t + h, ω + k[3])
		ω = ω +(k[1] + 2k[2] + 2k[3] + k[4])/6
		t = a + i*h
		push!(result, (t, ω))
	end
	return result
end

f(t, y) = y -t^2 + 1
println(rk4(f, 0, 2, 10, 0.5))


text:
# by JMW; taken from DataStructures.jl

macro delegate(source, targets)
    typename = esc(source.args[1])
    fieldname = source.args[2].value
    funcnames = targets.args
    n = length(funcnames)
    fdefs = Any[]
    for i in 1:n
        funcname = esc(funcnames[i])
        push!(fdefs, quote
            ($funcname)(a::($typename), args...) =
            ($funcname)(a.$fieldname, args...)
        end)
    end
    return Expr(:block, fdefs...)
end


text:
@testset "Active Set Core" begin

    # ActiveSetCore
    N = 10
    p = 3
    model = UnicycleGame(p=p)
    probsize = ProblemSize(N,model)
    ascore = ActiveSetCore(probsize)
    @test typeof(ascore) <: ActiveSetCore

end


text:
for oper in (:add, :sub, :mul, :div)
    fJ = Symbol(oper, :!)
    @eval begin
        export $fJ
    end

    for (specC, tJ, tC) in (("",    BigFloat,   Ref{BigFloat}),
                            ("_ui", CulongMax,  Culong),
                            ("_si", ClongMax,   Clong),
                            ("_d",  CdoubleMax, Cdouble),
                            ("_z",  BigInt,     Ref{BigInt}))

        @eval begin
            function ($fJ)(rop::BigFloat, op1::BigFloat, op2::$tJ, rnd=ROUNDING_MODE[])
                ccall(($(string("mpfr_", oper, specC)), "libmpfr"), Int32,
                      (Ref{BigFloat}, Ref{BigFloat}, $tC, MPFR.MPFRRoundingMode),
                      rop, op1, op2, rnd)
                return rop
            end
        end
    end
end

for oper in (:sub, :div)
    fJ = Symbol(oper, :!)
    for (specC, tJ, tC) in (("ui_", CulongMax,  Culong),
                            ("si_", ClongMax,   Clong),
                            ("d_",  CdoubleMax, Cdouble),
                            ("z_",  BigInt,     Ref{BigInt}))

        @eval begin
            function ($fJ)(rop::BigFloat, op1::$tJ, op2::BigFloat, rnd=ROUNDING_MODE[])
                ccall(($(string("mpfr_", specC,  oper)), "libmpfr"), Int32,
                      (Ref{BigFloat}, $tC, Ref{BigFloat}, MPFR.MPFRRoundingMode),
                      rop, op1, op2, rnd)
                return rop
            end
        end
    end
end

@doc """
    add!(rop, op1, op2, rnd)

Set `rop` to `op1+op2` rounded in the direction `rnd`. The IEEE 754
rules are used for signed zeros.

For types having no signed zeros, 0 is considered unsigned, i.e.,
`(+0)+0 = (+0)` and `(−0)+0 = (−0))`.
""" add!

@doc """
    sub!(rop, op1, op2, rnd)

Set `rop` to `op1-op2` rounded in the direction `rnd`. The IEEE 754
rules are used for signed zeros.

For types having no signed zeros, 0 is considered unsigned, i.e.,
`(+0)−0 = (+0)`, `(−0)−0 = (−0)`, `0−(+0) = (−0)` and `0−(−0) =
(+0))`.
""" sub!

@doc """
    mul!(rop, op1, op2, rnd)

Set `rop` to `op1×op2` rounded in the direction `rnd`. When a result
is zero, its sign is the product of the signs of the operands (for
types having no signed zeros, 0 is considered positive).
""" mul!

for f in (:sqr, :cbrt, :neg, :abs)
    fJ = Symbol(f, :!)
    @eval begin
        export $fJ
        function ($fJ)(rop::BigFloat, op::BigFloat, rnd=ROUNDING_MODE[])
            ccall(($(string("mpfr_", f)), "libmpfr"), Int32,
                  (Ref{BigFloat}, Ref{BigFloat}, MPFR.MPFRRoundingMode),
                  rop, op, rnd)
            return rop
        end
    end
end

for f in (:sqrt, :rec_sqrt)
    fJ = Symbol(f, :!)
    @eval begin
        export $fJ
        function ($fJ)(rop::BigFloat, op::BigFloat, rnd=ROUNDING_MODE[])
            op < 0 && throw(DomainError("Operand is negative"))
            ccall(($(string("mpfr_", f)), "libmpfr"), Int32,
                  (Ref{BigFloat}, Ref{BigFloat}, MPFR.MPFRRoundingMode),
                  rop, op, rnd)
            return rop
        end
    end
end

@doc """
    sqr!(rop, op, rnd)

Set `rop` to `op²` rounded in the direction `rnd`.
""" sqr!

@doc """
    sqrt!(rop, op, rnd)

Set `rop` to `√op` rounded in the direction `rnd`. Set `rop` to −0 if
`op` is −0, to be consistent with the IEEE 754 standard. Throw a
DomainError exception if `op` is negative.
""" sqrt!

@doc """
    rec_sqrt!(rop, op, rnd)

Set `rop` to 1/√op rounded in the direction `rnd`. Set `rop` to +Inf
if op is ±0,+0 if `op` is +Inf, and throw a DomainError exception if
`op` is negative. Warning!  Therefore the result on −0 is different
from the one of the `sqrt` function recommended by the IEEE 754-2008
standard (Section 9.2.1), which is −Inf instead of +Inf.
""" rec_sqrt!

@doc """
    cbrt!(rop, op, rnd)

Set `rop` to the cubic root of `op` rounded in the direction
`rnd`. For `op` negative (including −Inf), set `rop` to a negative
number. If `op` is zero, set `rop` to zero with the same sign as `op`.
""" cbrt!

export rootn!
"""
    rootn!(rop, op, k, rnd)

Set `rop` to the kth root of `op` rounded in the direction `rnd`. For
k=0, set `rop` to NaN. For k odd (resp. even) and `op` negative
(including −Inf), set `rop` to a negative number (resp. throw a
DomainError). If `op` is zero, set `rop` to zero with the sign
obtained by the usual limit rules, i.e., the same sign as `op` if `k`
is odd, and positive if `k` is even.

This function agrees with the `rootn` function of the IEEE 754-2008
standard (Section 9.2).
"""
function rootn!(rop::BigFloat, op::BigFloat, k::CulongMax, rnd=ROUNDING_MODE[])
    op < 0 && iseven(k) && throw(DomainError("Even root of negative number"))
    ccall((:mpfr_rootn_ui, "libmpfr"), Int32,
          (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFR.MPFRRoundingMode), rop, op, k, rnd)
    return rop
end

export pow!
for (specC, tJ, tC) in (("",    BigFloat,   Ref{BigFloat}),
                        ("_ui", CulongMax,  Culong),
                        ("_si", ClongMax,   Clong),
                        ("_z",  BigInt,     Ref{BigInt}))

    @eval begin
        function pow!(rop::BigFloat, op1::BigFloat, op2::$tJ, rnd=ROUNDING_MODE[])
            ccall(($(string("mpfr_pow", specC)), "libmpfr"), Int32,
                  (Ref{BigFloat}, Ref{BigFloat}, $tC, MPFR.MPFRRoundingMode),
                  rop, op1, op2, rnd)
            return rop
        end
    end
end

function pow!(rop::BigFloat, op1::CulongMax, op2::CulongMax, rnd=ROUNDING_MODE[])
    ccall((:mpfr_ui_pow_ui, "libmpfr"), Int32,
          (Ref{BigFloat}, Culong, Culong, MPFR.MPFRRoundingMode),
          rop, op1, op2, rnd)
    return rop
end

function pow!(rop::BigFloat, op1::CulongMax, op2::BigFloat, rnd=ROUNDING_MODE[])
    ccall((:mpfr_ui_pow, "libmpfr"), Int32,
          (Ref{BigFloat}, Culong, Ref{BigFloat}, MPFR.MPFRRoundingMode),
          rop, op1, op2, rnd)
    return rop
end

@doc """
    pow!(rop, op1, op2, rnd)

Set `rop` to op1^op2, rounded in the direction `rnd`.

Special values are handled as described inthe ISO C99 and IEEE
754-2008 standards for the `pow!` function.
""" pow!

@doc """
    neg!(rop, op, rnd)

Set `rop` to `−op`, rounded in the direction `rnd`. Just changes or
adjusts the sign if `rop` and `op` are the same variable, otherwise a
rounding might occur if the precision of `rop` is less than that of
`op`.
""" neg!

@doc """
    abs!(rop, op, rnd)

Set `rop` to the absolute value of `op`, rounded in the direction
`rnd`. Just changes or adjusts the sign if `rop` and `op` are the same
variable, otherwise a rounding might occur if the precision of `rop`
is less than that of `op`.
""" abs!

export dim!
"""
    dim!(rop, op1, op2, rnd)

Set `rop` to the positive difference of `op1` and `op2`,
i.e.,`op1−op2` rounded in the direction `rnd` if `op1>op2`, +0 if
`op1≤op2`, and throws a DomainError if `op1` or `op2` is NaN.
"""
function dim!(rop::BigFloat, op1::BigFloat, op2::BigFloat, rnd=ROUNDING_MODE[])
    isnan(op1) || isnan(op2) && throw(DomainError("Operands are NaNs"))
    ccall((:mpfr_dim, "libmpfr"), Int32,
          (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFR.MPFRRoundingMode),
          rop, op1, op2, rnd)
    return rop
end

for oper in (:mul_2, :div_2)
    fJ = Symbol(oper, :!)
    @eval begin
        export $fJ
    end
    for (specC, tJ, tC) in ((:ui, CulongMax,  Culong),
                            (:si, ClongMax,   Clong))

        @eval begin
            function ($fJ)(rop::BigFloat, op1::BigFloat, op2::$tJ, rnd=ROUNDING_MODE[])
                #ccall(($(string(:mpfr_, oper, specC)), "libmpfr"), Int32,
                ccall(($(string(:mpfr_, oper, specC)), "libmpfr"), Int32,
                      (Ref{BigFloat}, Ref{BigFloat}, $tC, MPFR.MPFRRoundingMode),
                      rop, op1, op2, rnd)
                return rop
            end
        end
    end
end

@doc """
    mul_2!(rop, op1, op2, rnd)

Set `rop` to `op1×2^op2` rounded in the direction `rnd`. Just
increases the exponent by `op2` when `rop` and `op1` are identical.
""" mul_2!

@doc """
    div_2!(rop, op1, op2, rnd)

Set `rop` to `op1/2^op2` rounded in the direction `rnd`. Just
decreases the exponent by `op2` when `rop` and `op1` are identical.
""" div_2!


text:
module WaveletRegression

using Wavelets
using Statistics
using LinearAlgebra

export
    wavelm,
    AbstractWaveletTransform,
    MODWT,
    DWT,
    DWPT,
    WaveletRegressionModel,
    modwt_matrix,
    imodwt_matrix,
    dwt_matrix,
    idwt_matrix,
    predict,
    predict_wavespace,
    coef,
    ncoef,
    dof,
    nlevels,
    residuals,
    regression_stderr,
    coef_stderr,
    coef_shrunk,
    shrink

export wavelet, WT

abstract type AbstractWaveletTransform end

for transformtype in [:MODWT, :DWT, :WPT]
    @eval begin
        struct $transformtype <: AbstractWaveletTransform
            wavelet::DiscreteWavelet
            nlevels::Integer
        end
        Wavelets.wavelet(trans::$(transformtype)) = trans.wavelet
        nlevels(trans::$(transformtype)) = trans.nlevels
    end
end

struct WaveletRegressionModel
    X::AbstractArray
    Y::AbstractArray
    n::Integer
    mx::Integer
    my::Integer
    Xw::AbstractArray
    Yw::AbstractArray
    trans::AbstractWaveletTransform
    coefficients
end
Wavelets.wavelet(mod::WaveletRegressionModel) = wavelet(mod.trans)
nlevels(mod::WaveletRegressionModel) = nlevels(mod.trans)
coef(mod::WaveletRegressionModel) = mod.coefficients
ncoef(mod::WaveletRegressionModel) = sum(length(B) for B in mod.coefficients)
dof(mod::WaveletRegressionModel) = mod.n - ncoef(mod) - 1

function predict_wavespace(mod::WaveletRegressionModel, newXw::AbstractArray, shrink=false)
    n = size(newXw, 1)
    if shrink
        BB = coef_shrunk(mod)
    else
        BB = mod.coefficients
    end
    Yw = cat([reshape(newXw[:, j, :] * BB[j], (n, 1, mod.my))
        for j in 1:nlevels(mod)+1]..., dims=2)
    return Yw
end

function predict(mod::WaveletRegressionModel, newX::AbstractVecOrMat, shrink=false)
    @assert size(newX, 2) == size(mod.X, 2)
    newXw = modwt_matrix(newX, mod.trans)
    Yw = predict_wavespace(mod, newXw, shrink)
    Y = imodwt_matrix(Yw, mod.trans)
    return Y
end
predict(mod::WaveletRegressionModel, shrink=false) = predict(mod, mod.X, shrink)

function residuals(mod::WaveletRegressionModel, shrink=false)
    Yw_pred = predict_wavespace(mod, mod.Xw, shrink)
    return mod.Yw .- Yw_pred
end

function regression_stderr(mod::WaveletRegressionModel; unbiased=true, shrink=false)
    resid = residuals(mod, shrink)
    s2 = reshape(sum(abs2 , resid, dims=1), (mod.my, :))
    if unbiased
        return s2 / dof(mod)
    else
        return s2 / mod.n
    end
end

"""
Calculate shrinkage factor for regression coefficients.
J.B. Copas (1983). Regression, prediction, and shrinkage. J. Royal Statistical
Society B 45(3), 311-354.  https://www.jstor.org/stable/2345402
"""
function coef_shrunk(mod::WaveletRegressionModel)
    BBshrunk = deepcopy(coef(mod))
    for j in 1:mod.mx
        Xw = mod.Xw[:, j, :]
        n = mod.n
        p = mod.mx * mod.my
        v = n - p - 1
        σ2 = regression_stderr(mod)[:, j]
        B = coef(mod)[j]
        V = Xw'Xw/n
        M = cholesky(V).U

        shrinkage = zeros(1, mod.my)
        for k in 1:mod.my
            β = B[:, k]
            ξ = M * β
            shrinkage[k] = 1 - ((p-2) * σ2[k] * v) / (n * (v+2) * ξ'ξ)
        end
        BBshrunk[j] = B .* shrinkage
    end
    return BBshrunk
end

function shrink(mod::WaveletRegressionModel)
    BB = coef_shrunk(mod)
    return WaveletRegressionModel(mod.X, mod.Y, mod.n, mod.mx, mod.my,
        mod.Xw, mod.Yw, mod.trans, BB)
end

function coef_stderr(mod::WaveletRegressionModel; unbiased=true, shrink=false)
    # https://en.wikipedia.org/wiki/Ordinary_least_squares
    s2 = regression_stderr(mod, unbiased=unbiased, shrink=shrink)
    L = mod.trans.nlevels + 1
    Qxx =  [mod.Xw[:, j, :]' * mod.Xw[:, j, :] for j in 1:L]
    SE = [similar(B) for B in coef(mod)]
    for j in 1:L
        for i in 1:mod.my
            SE[j][:, i] .= sqrt.(diag(s2[i, j] .* inv(Qxx[j])))
        end
    end
    return SE
end

function modwt_matrix(X::AbstractVecOrMat, trans::MODWT)
    return cat([modwt(X[:, i], trans.wavelet, trans.nlevels) for i in 1:size(X, 2)]..., dims=3)
end

function imodwt_matrix(Xw::AbstractArray, trans::MODWT)
    n, nlevels, nvars = size(Xw)
    hcat([imodwt(Xw[:, :, i], trans.wavelet) for i in 1:nvars]...)
end

function wavelm(X::AbstractVecOrMat, Y::AbstractVecOrMat, trans::MODWT)
    Xw = modwt_matrix(X, trans)
    Yw = modwt_matrix(Y, trans)
    BB = [Xw[:, j, :] \ Yw[:, j, :] for j in 1:(trans.nlevels+1)]
    n = size(X, 1)
    mx = div(length(X), n)
    my = div(length(Y), n)
    return WaveletRegressionModel(X, Y, n, mx, my, Xw, Yw, trans, BB)
end

function dwt_matrix(X::AbstractVecOrMat, trans::DWT)
    return hcat([dwt(X[:, i], trans.wavelet, trans.nlevels) for i in 1:size(X, 2)]...)
end

function idwt_matrix(Xw::AbstractArray, trans::DWT)
    n, nvars = size(Xw)
    hcat([idwt(Xw[:, i], trans.wavelet, trans.nlevels) for i in 1:nvars]...)
end

function wavelm(X::AbstractVecOrMat, Y::AbstractVecOrMat, trans::DWT)
    Xw = dwt_matrix(X, trans)
    Yw = dwt_matrix(Y, trans)
    n = size(Xw, 1)
    indices = [detailrange(n, j) for j in 1:trans.nlevels]
    push!(indices, 1:detailindex(n, trans.nlevels, 1)-1) # scaling coef indices
    BB = [Xw[ii, :] \ Yw[ii, :] for ii in indices]
    n = size(X, 1)
    mx = div(length(X), n)
    my = div(length(Y), n)
    return WaveletRegressionModel(X, Y, n, mx, my, Xw, Yw, trans, BB)
end

end # module


text:
using Base.Collections
using Base.Order

type HierarchicalSoftmaxNode
	parent::Int32
	branch::Bool
end

type HierarchicalOutput
	code::Array{Int8}
	path::Array{Int}
end

import Base.length

length(out::HierarchicalOutput) = length(out.path)

function HierarchicalSoftmaxNode()
	return HierarchicalSoftmaxNode(Int32(0), false)
end

function softmax_path(nodes::Array{HierarchicalSoftmaxNode},
		V::Integer, id::Integer)
	function path()
		while true
			node = nodes[id]
			if node.parent == 0 break; end
			@assert node.parent > V
			produce((Int32(node.parent - V), convert(Float64, node.branch)))
			id = node.parent
		end
	end

	return Task(path)
end

function build_huffman_tree{Tf <: Number}(freqs::Array{Tf})
	V = length(freqs)
	nodes = Array(HierarchicalSoftmaxNode, V)
	for v in 1:V
		nodes[v] = HierarchicalSoftmaxNode()
	end

	freq_ord = By(wf -> wf[2])
	heap = heapify!([(nodes[v], freqs[v]) for v in 1:V], freq_ord)

	function pop_initialize!(parent::Int, branch::Bool)
		node = heappop!(heap, freq_ord)
		node[1].parent = Int32(parent)
		node[1].branch = branch
		return node[2]
	end

	L = V
	while length(heap) > 1
		L += 1
		node = HierarchicalSoftmaxNode()
		push!(nodes, node)

		freq = 1
		freq = pop_initialize!(L, true) + pop_initialize!(L, false)
		heappush!(heap, (node, freq), freq_ord)
	end

	@assert length(heap) == 1

	return nodes
end

function convert_huffman_tree(nodes::Array{HierarchicalSoftmaxNode}, V::Integer)
	outputs = Array(HierarchicalOutput, V)
	for v in 1:V
		code = Array(Int8, 0)
		path = Array(Int, 0)

		for (n, branch) in softmax_path(nodes, V, v)
			push!(code, round(UInt8, branch))
			push!(path, n)
		end

		outputs[v] = HierarchicalOutput(code, path)
	end

	return outputs
end

export HierarchicalSoftmaxNode, softmax_path

max_stars_repo_name: glicerico/AdaGram

max_issues_repo_name: glicerico/AdaGram

max_forks_repo_name: glicerico/AdaGram


text:
using PyPlot

AF = A - B1 * F1 - B2 * F2
n = 20
x = Array{Float64}(3, n)
x[:, 1] = [1 1 1]
for t in 1:n-1
    x[:, t+1] = AF * x[:, t]
end
q1 = x[2, :]
q2 = x[3, :]
q = q1 + q2         # Total output, MPE
p = a0 - a1 * q     # Price, MPE

fig, ax = subplots(figsize=(9, 5.8))
ax[:plot](q, "b-", lw=2, alpha=0.75, label="total output")
ax[:plot](p, "g-", lw=2, alpha=0.75, label="price")
ax[:set_title]("Output and prices, duopoly MPE")
ax[:legend](frameon=false)


text:
# test/runtests.jl

using MIRTio
using Test: @test, detect_ambiguities

@test read_rdb_hdr(:test)

@test length(detect_ambiguities(MIRTio)) == 0


text:
"""
Generation — WaylandScanner module for generating Julia expressions with given WL protocol definitions.

Can be used in two ways - one can either evaluate the expressions immediately (preferably in a way that can be done once during precompilation) or serialise the generated code into a Julia script file to "lock" the current state of protocols.

The generated code assumes the core types and message-passing methods to already exist.
"""
module Generation

export genprotocol, generate_script, output_script

using ..WaylandScanner
using ..Parsing
import MacroTools, Dates

function gendocexpr(desc::SDescription)
	docstring = string(desc)
	if !isempty(docstring)
		quote
			$"""
			$docstring
			"""
		end
	else
		quote end
	end
end
gendocexpr(::Nothing) = quote end

"""
    genprotocol(protocol::SProtocol, serverside::Bool = false)

Generate the protocol module expression from parsed `protocol`.

`serverside` will cause client-side code (if set to false, default), or server-side code (if set to true) to be generated.
"""
function genprotocol(protocol::SProtocol, serverside::Bool = false)
	name = tostructid(protocol.name)
	quote
		$(gendocexpr(protocol.description))
		module $name
		$(geninterface.(protocol.interfaces, serverside)...)
		end # module $name
	end
end

"""
	geninterface(interface::SInterface, serverside::Bool = false)

Generate code defining an interface.
"""
function geninterface(interface::SInterface, serverside::Bool = false)
	name = tostructid(interface.name)
	quote
		$(gendocexpr(interface.description))
		struct $name
			id::WlID
		end
		$(genenum.(interface.enums, name)...)
		$((serverside ? genrequest_server.(interface.requests, name) : genrequest_client.(interface.requests, name))...)
		$((serverside ? genevent_server.(interface.events, name) : genevent_client.(interface.events, name))...)
		# end interface $name
	end
end

"""
    genenum(enum::SEnum, ifname::Symbol = "")

Generate code defining a Wayland enum. Sadly, Wayland protocols currently don't declare what enums are used by what, so the best we can do is implement them as simple dictionaries.
"""
function genenum(enum::SEnum, ifname::Symbol = "")
	name = Symbol(string(ifname, tostructid(enum.name)))
	quote
		$(gendocexpr(enum.description))
		const $name = Dict($(genenumpair.(values(enum.entries))...))
	end
end

"""
	genenumpair(entry::SEnumEntry)

Get a `name => value` pair from the `entry`.
"""
genenumpair(entry::SEnumEntry) = entry.name => entry.value

"""
    genrequest(request::SRequest, ifname::Symbol)

Generate code defining a client-side Wayland request method.

The generated method can be used by a client to send a request.
"""
function genrequest_client(request::SRequest, ifname::Symbol)
	name = Symbol(request.name)
	hasargs = request.args != nothing
	args = hasargs ? genarg.(request.args) : []
	argnames = hasargs ? argname.(request.args) : []
	quote
		$(gendocexpr(request.description))
		function $name(obj::$ifname, $(args...))
			message(obj, $(argnames...))
		end
	end
end

argname(arg::SArgument) = Symbol(arg.name)
genarg(arg::SArgument) = :($(argname(arg))::$(arg.type))

"""
    genevent_client(event::SEvent, ifname::Symbol)

Generate code used to listen to a Wayland event.
"""
function genevent_client(event::SEvent, ifname::Symbol)
	name = Symbol(event.name)
	quote
		$(gendocexpr(event.description))
		$"# TODO - event-related code for $name"
	end
end

"""
	generate_script(protocol::SProtocol, serverside::Bool = false)

Serialise the expression into a string.

Useful for writing the definitions into a Julia script file.
"""
function generate_script(protocol::SProtocol, serverside::Bool = false)
	header = """
	# Automatically generated by WaylandScanner's Generation module on $(Dates.now())
	"""
	string(header, '\n', MacroTools.flatten(MacroTools.striplines(genprotocol(protocol, serverside))))
end

"""
	output_script(protocol::SProtocol, filepath::Union{AbstractString,Nothing} = nothing, serverside::Bool = false)

Write the definitions script to a file.
"""
function output_script(protocol::SProtocol, filepath::AbstractString = "./$(protocol.name).jl", serverside::Bool = false)
	open(filepath, "w") do file
		println(file, generate_script(protocol, serverside))
	end
end

function output_script(protocols::Set{SProtocol}, args...)
	for prot in protocols
		output_script(prot, args...)
	end
end

end	# module Generation


text:
#####################
## Sequence OBJECT ##
#####################
"""
    Sequence(GR)

Sequence object.

# Argument
 - `GR::Array{Grad,2}` := Gradient array, fist axis is the dimension (x,y,z) and the second one is time.

# Examples
```julia-repl
julia> A = Sequence([Grad(1,1) for i=1:2,j=1:2])
Sequence(Grad[Grad(1, 1) Grad(1, 1); Grad(1, 1) Grad(1, 1)])

julia> B = Sequence([Grad(2,2) for i=1:2,j=1:2])
Sequence(Grad[Grad(2, 2) Grad(2, 2); Grad(2, 2) Grad(2, 2)])

julia> C = A + B
Sequence(Grad[Grad(1, 1) Grad(1, 1) Grad(2, 2) Grad(2, 2); Grad(1, 1) Grad(1, 1) Grad(2, 2) Grad(2, 2)])
```
"""
mutable struct Sequence
	GR::Array{Grad,2}		  #Sequence in (X, Y and Z) and time
	RF::Array{RF,2}			  #RF pulses in coil and time
	ADC::Array{ADC,1}		  #ADC in time
	DUR::Vector				  #Duration of each block, this enables delays after RF pulses to satisfy ring-down times
	DEF::Dict{String,Any} 	  #Dictionary with information relevant to the reconstructor
end
#MAIN CONSTRUCTORS
function Sequence(GR)	  #If no RF is defined, just use a zero amplitude pulse
	M,N = size(GR)
	Sequence([i <= M ? GR[i,j] : Grad(0, 0) for i=1:3, j=1:N],
		reshape([RF(0, 0) for i = 1:N],1,:),
		[ADC(0, 0) for i = 1:N],
		GR.dur,
		Dict()
		)
end
function Sequence(GR,RF)	#If no ADC is DEFined, just use a ADC with 0 samples
	@assert size(GR,2) .== size(RF,2) "The number of Gradient, and RF objects must be the same."
	M,N = size(GR)
	Sequence([i <= M ? GR[i,j] : Grad(0, 0) for i=1:3, j=1:N],
		RF,
		[ADC(0, 0) for i = 1:N],
		maximum([GR.dur RF.dur],dims=2)[:],
		Dict()
		)
end
function Sequence(GR,RF,ADC)
	@assert size(GR,2) .== size(RF,2) .== length(ADC) "The number of Gradient, RF, and ADC objects must be the same."
	M,N = size(GR)
	Sequence([i <= M ? GR[i,j] : Grad(0, 0) for i=1:3, j=1:N],
		RF,
		ADC,
		maximum([GR.dur RF.dur ADC.dur],dims=2)[:],
		Dict())
end
function Sequence(GR,RF,ADC,DUR)
	@assert size(GR,2) .== size(RF,2) .== length(ADC) .== length(DUR) "The number of Gradient, RF, ADC, and DUR objects must be the same."
	M,N = size(GR)
	Sequence([i <= M ? GR[i,j] : Grad(0, GR[1,j].T, GR[1,j].rise, GR[1,j].fall, GR[1,j].delay) for i=1:3, j=1:N],
		RF,
		ADC,
		maximum([GR.dur RF.dur ADC.dur DUR],dims=2)[:],
		Dict())
end
#OTHER CONSTRUCTORS
Sequence(GR::Array{Grad,1}) = Sequence(reshape(GR,1,:))
Sequence(GR::Array{Grad,1}, RF::Array{RF,1}) = Sequence(
														reshape(GR,:,1),
														reshape(RF,1,:),
														[ADC(0,0) for i = 1:size(GR,2)]
														)
Sequence(GR::Array{Grad,1}, RF::Array{RF,1}, A::ADC, DUR, DEF) = Sequence(
																reshape(GR,:,1),
																reshape(RF,1,:),
																[A],
																[DUR],
																DEF
																)
Sequence() = Sequence([Grad(0,0)])
#Sequence operations
Base.length(x::Sequence) = length(x.DUR)
Base.iterate(x::Sequence) = (Sequence(x.GR[:,1], x.RF[:,1], x.ADC[1], x.DUR[1], x.DEF), 2)
Base.iterate(x::Sequence, i::Integer) = (i <= length(x)) ? (Sequence(x.GR[:,i], x.RF[:,i], x.ADC[i], x.DUR[i], x.DEF), i+1) : nothing
Base.getindex(x::Sequence, i::UnitRange{Int}) = Sequence(x.GR[:,i], x.RF[:,i], x.ADC[i], x.DUR[i], x.DEF)
Base.getindex(x::Sequence, i::Int) = Sequence(x.GR[:,i], x.RF[:,i], x.ADC[i], x.DUR[i], x.DEF)
Base.getindex(x::Sequence, i::BitArray{1}) = any(i) ? Sequence(x.GR[:,i], x.RF[:,i], x.ADC[i], x.DUR[i], x.DEF) : nothing
Base.getindex(x::Sequence, i::Array{Bool,1}) = any(i) ? Sequence(x.GR[:,i], x.RF[:,i], x.ADC[i], x.DUR[i], x.DEF) : nothing
Base.lastindex(x::Sequence) = length(x.DUR)
Base.copy(x::Sequence) where Sequence = Sequence([deepcopy(getfield(x, k)) for k ∈ fieldnames(Sequence)]...)
Base.show(io::IO,s::Sequence) = begin
	compact = get(io, :compact, false)
	if !compact
		nGRs = sum(is_Gx_on.(s)) + sum(is_Gy_on.(s)) + sum(is_Gz_on.(s))
		print(io, "Sequence[ τ = $(round(dur(s)*1e3;digits=3)) ms | blocks: $(length(s)) | ADC: $(sum(is_ADC_on.(s))) | GR: $nGRs | RF: $(sum(is_RF_on.(s))) | DEF: $(length(s.DEF)) ]")
	else
		print(io, "Sequence[τ = $(round(dur(s)*1e3;digits=3)) ms]")
	end
end
#Arithmetic operations
+(x::Sequence, y::Sequence) = Sequence([x.GR  y.GR], [x.RF y.RF], [x.ADC; y.ADC], [x.DUR; y.DUR], merge(x.DEF, y.DEF))
-(x::Sequence, y::Sequence) = Sequence([x.GR -y.GR], [x.RF y.RF], [x.ADC; y.ADC], [x.DUR; y.DUR], merge(x.DEF, y.DEF))
-(x::Sequence) = Sequence(-x.GR, x.RF, x.ADC, x.DUR, x.DEF)
*(x::Sequence, α::Real) = Sequence(α*x.GR, x.RF, x.ADC, x.DUR, x.DEF)
*(α::Real, x::Sequence) = Sequence(α*x.GR, x.RF, x.ADC, x.DUR, x.DEF)
*(x::Sequence, α::ComplexF64) = Sequence(x.GR, α.*x.RF, x.ADC, x.DUR, x.DEF)
*(α::ComplexF64, x::Sequence) = Sequence(x.GR, α.*x.RF, x.ADC, x.DUR, x.DEF)
*(x::Sequence, A::Matrix{Float64}) = Sequence(A*x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
*(A::Matrix{Float64}, x::Sequence) = Sequence(A*x.GR, x.RF, x.ADC, x.DUR, x.DEF) #TODO: change this, Rotation fo waveforms is broken
/(x::Sequence, α::Real) = Sequence(x.GR/α, x.RF, x.ADC, x.DUR, x.DEF)
#Sequence object functions
size(x::Sequence) = size(x.GR[1,:])
is_ADC_on(x::Sequence) = any(x.ADC.N .> 0)
"Tells if the sequence has elements with ADC on during t."
is_ADC_on(x::Sequence, t::Union{Array{Float64,1},Array{Float64,2}}) = begin
	N = length(x)
	ΔT = durs(x)
	ts = cumsum([0 ; ΔT[1:end-1]])
	delays = x.ADC.delay
	Ts = 	 x.ADC.dur #delat+T
	t0s = ts .+ delays
	tfs = ts .+ Ts
	# The following just checks the ADC
	# |___∿  |
	#     △
	#     Here
	activeADC = any([is_ADC_on(x[i]) && any(t0s[i] .<= t .< tfs[i]) for i=1:N])
	activeADC
end
"Tells if the sequence has elements with RF on."
is_RF_on(x::Sequence) = any([sum(abs.(r.A)) for r = x.RF] .> 0)
"Tells if the sequence has elements with RF on during t."
is_RF_on(x::Sequence, t::Vector{Float64}) = begin
	N = length(x)
	ΔT = durs(x)
	ts = cumsum([0 ; ΔT[1:end-1]])
	delays = x.RF.delay
	Ts = 	 x.RF.dur #dur = delat+T
	t0s = ts .+ delays
	tfs = ts .+ Ts
	# The following just checks the RF waveform
	# |___∿  |
	#     △
	#     Here
	activeRFs = any([is_RF_on(x[i]) && any(t0s[i] .<= t .<= tfs[i]) for i=1:N])
	activeRFs
end
"Tells if the sequence has elements with GR on."
is_GR_on(x::Sequence) = any([sum(abs.(r.A)) for r = x.GR] .> 0)
"Tells if the sequence has elements with GR on."
is_Gx_on(x::Sequence) = any([sum(abs.(r.A)) for r = x.GR.x] .> 0)
"Tells if the sequence has elements with GR on."
is_Gy_on(x::Sequence) = any([sum(abs.(r.A)) for r = x.GR.y] .> 0)
"Tells if the sequence has elements with GR on."
is_Gz_on(x::Sequence) = any([sum(abs.(r.A)) for r = x.GR.z] .> 0)
"Tells if the sequence is a delay"
is_Delay(x::Sequence) = !(is_GR_on(x) || is_RF_on(x) || is_ADC_on(x))
"Duration of sequence's blocks `ΔT::Array{Real,1}` [s]."
durs(x::Sequence) = begin
	# ΔT_GR  = x.GR.dur
	# ΔT_RF  = x.RF.dur
	# ΔT_ADC = x.ADC.dur
	# ΔT_DUR = x.DUR
	# ΔT = maximum([ΔT_GR ΔT_RF ΔT_ADC ΔT_DUR],dims=2)
	x.DUR
end
"Duration of the sequence `T` [s]."
dur(x::Sequence) = sum(durs(x))

#Trapezoidal Waveform
⏢(A,t,ΔT,ζ1,ζ2,delay) = begin
	if sum(abs.(A)) != 0 && ΔT+ζ1+ζ2 != 0 # If no event just ignore calculations
		#Getting amplitudes, onlu supports uniformly sampled waveforms for now
		if length(A) != 1
			grad_raster = ΔT / length(A)
			idx = ceil.(Int, (t .- delay .- ζ1) ./ grad_raster ) #Time to integer index
			valid = 1 .<= idx .<= length(A)
			idx[(!).(valid)] .= 1
			B = A[idx] .* valid
		else
			B = A
		end
		#Trapezoidal waveform
		aux = (ζ1    .<= t .- delay .< ζ1+ΔT) .* B
		if ζ1 != 0
			aux .+= (0     .<= t .- delay .< ζ1) .* A[1] .* (t .- delay)./ζ1
		end
		if ζ2 !=0
			aux .+= (ζ1+ΔT .<= t .- delay .< ζ1+ΔT+ζ2) .* A[end] .* (1 .- (t.-delay.-ζ1.-ΔT)./ζ2)
		end
	else
		aux = zeros(size(t))
	end
	aux
end
"""
	get_grad(seq,t)

Get gradient array from sequence `seq` evaluated in time points `t` for the dimension `dim`.

# Arguments
 - `seq::Sequence`:= Sequence object.
 - `t`            := Time vector to evaluate.

# Examples
```julia-repl
julia> seq = Sequence([Grad(1,1) Grad(-1,1); Delay(1) Delay(1)])
Sequence(Grad[Grad(1, 1) Grad(-1, 1); Grad(0, 1) Grad(0, 1)])

julia> t = 0:.5:2
0.0:0.5:2.0

julia> get_grad(seq,t)
5-element Array{Float64,1}:
  1.0
  1.0
  0.0
 -1.0
 -1.0
```
"""
function get_grads(seq, t::Vector)
    gx = get_theo_Gi(seq,1)
    gy = get_theo_Gi(seq,2)
    gz = get_theo_Gi(seq,3)
    Gx = LinearInterpolation(gx..., extrapolation_bc=0)(t)
    Gy = LinearInterpolation(gy..., extrapolation_bc=0)(t)
    Gz = LinearInterpolation(gz..., extrapolation_bc=0)(t)
    (Gx, Gy, Gz)
end
function get_grads(seq, t::Matrix)
	t_vec = t[:]
    gx = get_theo_Gi(seq,1)
    gy = get_theo_Gi(seq,2)
    gz = get_theo_Gi(seq,3)
    Gx = LinearInterpolation(gx..., extrapolation_bc=0)(t_vec)
    Gy = LinearInterpolation(gy..., extrapolation_bc=0)(t_vec)
    Gz = LinearInterpolation(gz..., extrapolation_bc=0)(t_vec)
    (Gx', Gy', Gz')
end
# get_grads(seq::Sequence,t) = begin
# 	#Amplitude
# 	A = seq.GR.A
# 	#Grad Timings
# 	T = seq.GR.T
# 	ζ1 = seq.GR.rise
# 	ζ2 = seq.GR.fall
# 	delay = seq.GR.delay
# 	#Sequence timings
# 	ΔT = durs(seq) #Duration of sequence block
# 	T0 = cumsum([0; ΔT[:]]) #Start time of each block
# 	#Waveforms
# 	(sum([⏢(A[j,i],t.-T0[i],T[j,i],ζ1[j,i],ζ2[j,i],delay[j,i]) for i=1:length(seq)]) for j=1:3)
# end

get_rfs(seq::Sequence,t) = begin
	#Amplitude
	A  = seq.RF.A
	Δf = seq.RF.Δf
	#Timings
	T = seq.RF.T
	delay = seq.RF.delay
	T0 = cumsum([0; durs(seq)], dims=1)
	(sum([⏢(A[1,i], t.-T0[i],sum(T[i]),0,0,delay[i]) for i=1:length(seq)]),
	 sum([⏢(Δf[1,i],t.-T0[i],sum(T[i]),0,0,delay[i]) for i=1:length(seq)])
	)
end
get_flip_angles(x::Sequence) = get_flip_angle.(x.RF)[:]

get_ADC_on(seq::Sequence,t::Array{Float64,1}) = begin
	Δt = t[2]-t[1]
	M, N = size(seq.ADC)
	ADC = seq.ADC.N .> 0
	T = [seq.ADC[i].T for i=1:N]; T = [0; T]
	⊓(t) = (0 .<= t .< 1)*1.;
	ADC_bool = sum([ADC[i]*⊓((t.-sum(T[1:i]))/(T[i+1]+Δt)) for i=1:N])
	ADC_bool = (ADC_bool.>0)
	circshift(convert.(Bool,	[:]),-1)
end
"Calculates the `b`-value, in the PGSE sequnce is b = (2πγGδ)²(Δ-δ/3) and the normalized diffusion vector `n` from a *diffusion sequence*."
get_bvalue(DIF::Sequence) = begin
	M, N = size(DIF.GR)
	G = getproperty.(DIF.GR,:A) #[DIF.GR[i,j].A for i=1:M,j=1:N] #Strength of pulse
	δ = getproperty.(DIF.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
	τ = dur(DIF) #End of sequence
	# B-value
	b = 0
	for k=1:M,i=1:N,j=1:N
		ij = max(i,j)
		α = (i==j) ? 2/3 : 1/2
		b += G[k,i]*G[k,j]*δ[i]*δ[j]*(τ-T[ij]-α*δ[ij])
	end
	b_value = (2π*γ)^2*b # Trace of B tensor
end
"Calculates the `b`-matrix, such as `b`-value = g' B g [s/mm2] with g [T/m]."
get_Bmatrix(DIF::Sequence) = begin
	M, N = size(DIF.GR)
	δ = getproperty.(DIF.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
	τ = dur(DIF) #End of sequence
	# B-value, slower way of doing it
	# b = zeros(M,N,N)
	# for k=1,i=1:N,j=1:N
	# 	ij = max(i,j)
	# 	α = (i==j) ? 2/3 : 1/2
	# 	b[k,i,j] = δ[i]*δ[j]*(τ-T[ij]-α*δ[ij])
	# end
	ij = [max(i,j) for i=1:N, j=1:N]
	α = [(i==j) ? 2/3 : 1/2 for i=1:N, j=1:N]
	b = (δ' .* δ) .* (τ .- T[ij] .- α .* δ[ij])
	b_value = (2π*γ)^2*1e-6*b # Trace of B tensor
	b_value
end

"Calculates the `b` value for PGSE and a moment nulled sequence."
get_bvalue(G,Δ,δ;null::Bool=false) = (null ? 1/9 : 1)*(2π*γ*G*δ)^2*(Δ-δ/3)
"Calculates `q` = γδG diffusion vector from a *diffusion sequence*."
get_qvector(DIF::Sequence;type::String="val") = begin
	M, N = size(DIF.GR)
	G = getproperty.(DIF.GR,:A) #[DIF.GR[i,j].A for i=1:M,j=1:N] #Strength of pulse
	δ = getproperty.(DIF.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
	τ = dur(DIF) #End of sequence
	qτ = γ*sum([G[i,j]*δ[j] for i=1:M,j=1:N],dims=2) #Final q-value
	norm(qτ) ≥ 1e-10 ? error("This is not a diffusion sequence, M0{G(t)}=∫G(t)dt!=0.") : 0
	# q-trajectory
	int_rect(t,A,δ) = (t.≥0) ? A*(t.-(t.≥δ).*(t.-δ)) : 0
	q(t) = γ*sum([int_rect(t.-T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
 	# Output q-vector and q-trajectory
	if type=="val"
		q(τ/2) #1/m
	elseif type=="traj"
		t = range(0,stop=τ,length=10^3)
		qt = [q(t)[i] for i=1:M, t=t]
		l = PlotlyJS.Layout(;title="q(t) trajectory", yaxis_title="q [1/mm]",
	    xaxis_title="t [ms]",height=300)
		p = [PlotlyJS.scatter() for j=1:M]
		idx = ["qx" "qy" "qz"]
		for j=1:M
			p[j] = PlotlyJS.scatter(x=t*1e3, y=qt[j,:]*1e-3,name=idx[j])
		end
		PlotlyJS.plot(p, l)
	end
end
get_M0_M1_M2(SEQ::Sequence) = begin
	M, N = size(SEQ.GR)
	G = [SEQ.GR[i,j].A for i=1:M,j=1:N] #Strength of pulse
	δ = [SEQ.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T] #Position of pulse
	τ = dur(SEQ) #End of sequence
	#M0
	M0i(t,T,A,δ) = (t.≥T) ? A*(t.-T.-(t.≥δ+T).*(t.-δ.-T)) : 0
	M0(t) = sum([M0i(t,T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
	#M1
	M1i(t,T,A,δ) = (t.≥T) ? A/2*(t.^2 .-T^2 .-(t.≥δ+T).*(t.^2 .-(δ+T).^2))./2 : 0
	M1(t) = sum([M1i(t,T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
	#M1
	M2i(t,T,A,δ) = (t.≥T) ? A/3*(t.^3 .-T^3 .-(t.≥δ+T).*(t.^3 .-(δ+T).^3))./3 : 0
	M2(t) = sum([M2i(t,T[j],G[i,j],δ[j]) for i=1:M,j=1:N],dims=2)
	M0, M1, M2
end
get_max_grad(x::Sequence) = begin
	G = [x.GR[i,j].A for i=1:size(x.GR,1),j=1:size(x.GR,2)]
	Gmax = maximum(sqrt.(sum(G.^2,dims=1)))
end
"""Get RF centers and types. Useful for k-space calculations."""
function get_RF_types(seq, t)
	α = get_flip_angles(seq)
	RF_mask = is_RF_on.(seq)
	RF_ex = (α .<= 90.01) .* RF_mask
	RF_rf = (α .>  90.01) .* RF_mask
	rf_idx = Int[]
	rf_type = Int[]
	T0 = cumsum([0; durs(seq)[:]])
	for i = 1:length(seq)
		if is_RF_on(seq[i])
			trf = get_RF_center(seq[i].RF[1]) + T0[i]
			append!(rf_idx, argmin(abs.(trf.-t)))
			if RF_ex[i]
				append!(rf_type, 0)
			elseif RF_rf[i]
				append!(rf_type, 1)
			end
		end
	end
	rf_idx, rf_type
end
"Outputs designed k-space trajectory from sequence object."
get_kspace(seq::Sequence; Δt=1) = begin
	t, Δt = get_uniform_times(seq, Δt)
	Gx, Gy, Gz = get_grads(seq, t)
	G = Dict(1=>[Gx; 0], 2=>[Gy; 0], 3=>[Gz; 0])
	#kspace
	Nt = length(t)
	k = zeros(Nt,3)
	#get_RF_center_breaks
	idx_rf, rf_type = get_RF_types(seq, t)
	parts = kfoldperm(Nt,1,type="ordered", breaks=idx_rf)
	for i = 1:3
		kf = 0
		for (rf, p) in enumerate(parts)
			k[p,i] = cumtrapz(Δt[p]', G[i][p[1]:p[end]+1]')[:] #This is the exact integral
			if rf > 1 #First part does not have RF
				k[p,i] .-= rf_type[rf-1] * kf
			end
			kf = k[p[end],i]
		end
	end
	kspace = γ * k #[m^-1]
	#Interp, as Gradients are piece-wise linear, the integral is piece-wise quadratic
	#Nevertheless, the integral is sampled at the ADC times so a linear interp is sufficient
	ts = t .+ Δt
	t_adc =  get_sample_times(seq)
	kx_adc = LinearInterpolation(ts,kspace[:,1])(t_adc)
	ky_adc = LinearInterpolation(ts,kspace[:,2])(t_adc)
	kz_adc = LinearInterpolation(ts,kspace[:,3])(t_adc)
	kspace_adc = [kx_adc ky_adc kz_adc]
	#Final
	kspace, kspace_adc
end

""""Calculates the normalized moments at the end of the sequence. """
function get_Mmatrix(seq::Sequence;order=0)
    M, N = size(seq.GR)
    τ = dur(seq)
    δ = getproperty.(seq.GR,:T)[1,:] #[DIF.GR[1,j].T for j=1:N]; #Duration of pulse
	T = [sum(δ[1:j]) for j=1:N]; T = [0; T[1:end-1]] #Position of pulse
	if order == 0
		M0 = δ/τ
		M1 = δ.*(T .+ δ/2)/τ^2
		M2 = δ.*(T.^2 .+ T.*δ .+ δ.^2/3)/τ^3
		M3 = δ.*(T.^3 .+ 3/2 * T.^2 .*δ .+ T.*δ.^2 .+ δ.^3/4)/τ^4
	elseif order == 1
		M0 = δ/τ
		M1 = δ.*T/τ^2
		M2 = δ.*(T.^2 .+ δ.^2/6)/τ^3
		M3 = δ.*(T.^3 .+ T.*δ.^2/2)/τ^4
	end
    [M0'; M1'; M2'; M3']
end

using LinearAlgebra: I, Bidiagonal, norm
"""Slew rate matrix: |SR*g| ≤ Smax."""
function get_SRmatrix(seq::Sequence)
    _, N = size(seq.GR)
    T = getproperty.(seq.GR[1,:],:T)
    Δt = [T[i] for i=1:N-1]
    dv = [Δt; T[end]]
    ev = Δt;
	SR = Array(Bidiagonal(-1 ./ dv, 1 ./ ev, :U))
	SR = [SR[1,:]'*0 ; SR]; SR[1,1] = T[1]
	SR
end

## TO SCANNER (Philips)
"""Duration in [s] => samples, with dwell-time of Δt = 6.4 μs."""
δ2N(δ) = floor(Int64, δ * 156250) + 2

function write_diff_fwf(DIF, idx180, Gmax, bmax; filename="./qte_vectors_input.txt", name="Maxwell2")
    open(filename, "w") do io
        G1 = DIF[1:idx180-1]
        G2 = DIF[idx180 + 1:end]
        t1 = range(0,dur(G1),length=δ2N(dur(G1)))
		t2 = range(0,dur(G2),length=δ2N(dur(G2)))
		Gx2, Gy2, Gz2 = get_grads_linear(G2,t2)
		Gx1, Gy1, Gz1 = get_grads_linear(G1,t1)
        write(io, "$name $(δ2N(dur(G1))) $(δ2N(dur(G2))) $(bmax)\n")
        for i = 1:length(t1)
            fx1, fy1, fz1 = (Gx1[i], Gy1[i], Gz1[i])./Gmax
            fx2, fy2, fz2 = i ≤ length(t2) ? (-Gx2[i], -Gy2[i], -Gz2[i])./Gmax : (0,0,0)
            line = @sprintf "% .4f % .4f % .4f % .4f % .4f % .4f\n" fx1 fy1 fz1 fx2 fy2 fz2
            write(io, line)
        end
    end
end

function read_diff_fwf(filename="./qte_vectors_input.txt")
    f = readlines(filename)
    G = zeros(size(f[2:end],1),6)
    n1, n2 = tryparse.(Int64, split(f[1]," ")[2:end])
    for (i, l) in enumerate(f[2:end])
        g = tryparse.(Float64, split(l," "))
        G[i,:] = g[g .≠ nothing]
    end
    G, n1, n2
end


text:
export mapflute

"""
 𝒍 = mapflute(𝒇, 𝒅)

Calculate locations of toneholes, and length of flute
 for frequencies 𝒇 (Hz, descending pitch) and hole diameters 𝒅 (mm, same order)
 in mm measured from embouchure hole
"""
function mapflute(𝒇, 𝒅)
  # TODO: assumes standard temperature ϑ, bore diameter ⌀, wall thickness ℎ
  # hole positions
  𝒍 = []
  # number of holes
  ħ = length(𝒇)-1
  # accumulate closed-hole correction
  𝛥ℓᵥ = 0.0
  for h in 1:ħ
    𝑓 = 𝒇[h]    # opened hole frequency
    𝑓ₜ = 𝒇[h+1] # closed hole frequency
    𝑑ₕ = 𝒅[h]   # hole diameter
    ℓₕ = toneholelength(𝑓, 𝑓ₜ, 𝑑ₕ; 𝛥ℓᵥ=𝛥ℓᵥ) # calculate hole position
    𝒍 = [𝒍; ℓₕ] # append to result list
    # closed hole correction effect on remaining holes
    𝛥ℓᵥ += closedholecorrection(𝑓, 𝑓ₜ, 𝑑ₕ; 𝛥ℓᵥ=𝛥ℓᵥ)
  end
  # append flute tube length to result
  ℓₜ = flutelength(𝒇[end]; 𝛥ℓᵥ=𝛥ℓᵥ)
  return [𝒍; ℓₜ]
end

max_stars_repo_name: isaacwhanson/scad-flute

max_issues_repo_name: isaacwhanson/scad-flute

max_forks_repo_name: isaacwhanson/scad-flute


text:
using SwitchTimeOpt

# Plotting
using PyPlot, PyCall
close("all")
# Import Seaborn for nice plots
@pyimport seaborn as sns
sns.set_palette("hls", 8)
sns.set_context("paper", font_scale=1.5)
sns.set_style("whitegrid")
# Use Latex Labels in Plots
plt[:rc]("text", usetex=true)
plt[:rc]("font", family="serif")   # Use Serif math

maxiter = 15
using Ipopt
solver = IpoptSolver(
  print_level=0,
  tol = 1e-10,  # Increase required precision to avoid convergence before maxiter iterations
  max_iter = maxiter)

### Define system parameters
#  Time vector
t0 = 0.0
tf = 10.0

# Integer inputs
uvec = repmat([1.0; 2.0], 8, 1)'

# Number of switches
N = size(uvec,2) -1

# Cost function matrix
K = 1.0
C = sqrt(K)*[0.0 1.0 -1.0]
Q = C'*C

# Initial state
x0 = [2.0; 2.0; 3.0]

### Define system dynamics
function nldyn(x::Array{Float64,1}, u::Array{Float64,1})
  n = length(x)
  f = zeros(n)

  f[1] = -sqrt(x[1]) + u[1]
  f[2] = sqrt(x[1]) - sqrt(x[2])
  f[3] = -0.05
  return f

end

function nldyn_deriv(x::Array{Float64,1}, u::Array{Float64,1})
  df = [-1/(2*sqrt(x[1]))        0                      0;
        1/(2*sqrt(x[1]))        -1/(2*sqrt(x[2]))       0;
        0                        0                      0]
end

# Generate and solve problems with different grid points
ngrid = [10; 30; 50; 100]

# Preallocate vectors for results
objode45 = Array{Float64}(length(ngrid))
objlin = Array{Float64}(length(ngrid))
objiterates = Array{Float64}(maxiter+1, length(ngrid))
cputime = Array{Float64}(length(ngrid))
nobjeval = Array{Int}(length(ngrid))
ngradeval = Array{Int}(length(ngrid))
nhesseval = Array{Int}(length(ngrid))
tauopt = Array{Float64}(N, length(ngrid))
xsim = Array{Float64}(3, 10^4, length(ngrid))
xlinsim = Array{Float64}(3, 10^4, length(ngrid))
usim = Array{Float64}(1, 10^4, length(ngrid))

# Initialize the problem first
m = stoproblem(x0, nldyn, nldyn_deriv, uvec)

for i = 1:length(ngrid) # Iterate over all grid points numbers
  m = stoproblem(
    x0,                 # Initial state
    nldyn,              # Nonlinear dynamics
    nldyn_deriv,        # Nonlinear dynamics derivative
    uvec,               # Vector of integer inputs
    ngrid = ngrid[i],   # Number of grid points
    t0 = t0,            # InitialtTime
    tf = tf,            # Final time
    Q = Q,              # Cost matrix
    solver = solver)

  # Solve optimization
  solve!(m)

  # Obtain values
  tauopt[:, i] = gettau(m)
  objlin[i] = getobjval(m)
  cputime[i] = getsoltime(m)
  nobjeval[i] = getnobjeval(m)
  ngradeval[i] = getngradeval(m)
  nhesseval[i] = getnhesseval(m)

  # Simulate System
  # Nonlinear Simulation
  xsim[:, :, i], xpts, objode45[i], t = simulate(m, tauopt[:, i])
  usim[:, :, i], _ = simulateinput(m, t)

  # Linearized simulation
  xlinsim[:, :, i], _, Jlinsim, _ = simulatelinearized(m, tauopt[:, i], t)

  # Save objective iterates
  objiterates[:, i] = m.STOev.obj[2:end]

end

###  Print results
@printf("RESULTS\n")
@printf("-------\n\n")

@printf("+==================================================================================+\n")
@printf("|  ngrid   |   objode45  |   objlin    | deltaobj [%%] |  nobjeval  |  cputime [s]  |\n")
@printf("+==================================================================================+\n")

for i = 1:length(ngrid)
  @printf("|  %7.i | %9.4f   | %9.4f   |   %6.3f     | %9.i  | %12.4f  | \n", ngrid[i], objode45[i], objlin[i], 100*norm(objode45[i]- objlin[i])/objode45[i], nobjeval[i], cputime[i])
  @printf("+----------------------------------------------------------------------------------+\n")
end

@printf("\nOptimal Switching Times for ngrid = 10\n")
@printf("--------------------------------------\n")
@printf("tauopt = "); show(round.(tauopt[:, 1],2)); @printf("\n")

# Generate plots for ngrid = 10
t = linspace(t0, tf, 10000)

figure()
subplot(3,1,1)
plot(t, xlinsim[1,:, 1], sns.xkcd_rgb["grass green"], linestyle = "dashdot")
plot(t, xsim[1,:, 1], sns.xkcd_rgb["denim blue"])
ylim(1.5, 4)
yticks([2; 3])
ylabel(L"x_1")

subplot(3,1,2)
plot(t, xsim[3,:, 1], sns.xkcd_rgb["black"], linestyle = "dotted")
plot(t, xlinsim[2,:, 1], sns.xkcd_rgb["grass green"], linestyle = "dashdot")
plot(t, xsim[2,:, 1], sns.xkcd_rgb["denim blue"])
ylabel(L"x_2")
ylim(1.5, 4.0)
yticks([2; 3])

subplot(3,1,3)
ax1 = plot(t, usim[1,:, 1], sns.xkcd_rgb["denim blue"])
ylim(0.8, 2.2)
yticks([1; 2])
ax = gca()
ax[:set_yticklabels]([L"u_{\mathrm{min}}", L"u_{\mathrm{max}}"])
ylabel(L"u")
xlabel(L"$\mathrm{Time}\; [s]$")

## Save figure
# tight_layout()
# savefig("tank_problem.pdf")

# Do not return anything
nothing


text:
@polly function kernel_trmm(alpha, A, B)
    m,n = size(B)
    for i = 1:m, j = 1:n
        for k = (i+1):m
            B[i,j] += A[k,i] * B[k,j]
        end
        B[i,j] = alpha * B[i,j]
    end
end

let
    m = 1000
    n = 1200

    A = zeros(Float32,m,m)
    B = zeros(Float32,m,n)

    alpha = Float32(1.5)
    for i = 1:m
        for j = 1:(i-1)
            A[i,j] = ((i+j) % m)/m
        end
        A[i,i] = 1.0
        for j = 1:n
            B[i,j] = ((n+(i-j)) % n)/n
        end
    end

    SUITE["trmm"] = @benchmarkable kernel_trmm($alpha, A, B) setup = (A = copy($A); B = copy($B))
end


text:
function createconstants!(dataobject::InfoType)
    dataobject.constants = createconstants()
    return dataobject
end

function createconstants()

    #---------------------------------------------------
    # define constants in cgs units
    #---------------------------------------------------
    # Sources:
    # http://www.astro.wisc.edu/~dolan/constants.html
    # IAU
    # RAMSES
    constants = PhysicalUnitsType() #zeros(Float64, 17)...)
    constants.Au = 149597870700e-13    # [cm] Astronomical unit -> from IAU
    constants.pc = 3.08567758128e18    # [cm] Parsec -> from IAU
    constants.kpc = constants.pc * 1e3
    constants.Mpc = constants.pc * 1e6
    constants.mpc = constants.pc * 1e-3

    constants.ly = 9.4607304725808e17  # [cm] Light year -> from IAU
    constants.Msol = 1.9891e33         # [g] Solar mass -> from IAU
    constants.Rsol = 6.96e10 #cm: Solar radius
    # Lsol = #erg s-2: Solar luminosity
    constants.Mearth = 5.9722e27       # [g]  Earth mass -> from IAU
    constants.Mjupiter = 1.89813e30    # [g]  Jupiter -> from IAU

    constants.me = 9.1093897e-28 #g: electron mass
    constants.mp = 1.6726231e-24 #g: proton mass
    constants.mn = 1.6749286e-24 #g: neutron mass
    constants.mH = 1.66e-24            # [g]   H-Atom mass -> from RAMSES
    constants.amu = 1.6605402e-24 #g: atomic mass unit
    constants.NA = 6.0221367e23 # Avagadro's number
    constants.c = 2.99792458e10 #cm s-1: speed of light in a vacuum
    # h = #erg s: Planck constant
    # hbar = #erg s
    constants.G  = 6.67259e-8 # cm3 g-1 s-2 Gravitational constant
    constants.kB = 1.3806200e-16       # [cm2 g s-2 K-1] Boltzmann constant -> cooling_module.f90 RAMSES

    constants.yr  = 3.15576e7           # [s]  Year -> from IAU
    constants.Myr = constants.yr *1e6
    constants.Gyr = constants.yr *1e9

    return constants
end

"""
### Create an object with predefined scale factors from code to pysical units
```julia
function createscales!(dataobject::InfoType)

return ScalesType
```
"""
function createscales!(dataobject::InfoType)
    dataobject.scale = createscales(dataobject)
    return dataobject
end

# create scales-field from existing InfoType
function createscales(dataobject::InfoType)
    unit_l = dataobject.unit_l
    unit_d = dataobject.unit_d
    unit_t = dataobject.unit_t
    unit_m = dataobject.unit_m
    constants = dataobject.constants
    return createscales(unit_l, unit_d, unit_t, unit_m, constants)
end

function createscales(unit_l::Float64, unit_d::Float64, unit_t::Float64, unit_m::Float64, constants::PhysicalUnitsType)
    #Initialize scale-object
    scale = ScalesType() #zeros(Float64, 32)...)

    # Conversion factors from user units to astronomical units
    mH      =   constants.mH        # [g]   H-Atom mass -> from RAMSES
    kB      =   constants.kB        # [cm2 g s-2 K-1] = [erg K-1] Boltzmann constant -> cooling_module.f90 RAMSES
    #Mpc     =   constants.pc /1e6   # [cm] MegaParsec -> from IAU
    #kpc     =   constants.pc /1e3   # [cm] KiloParsec -> from IAU
    pc      =   constants.pc        # [cm] Parsec -> from IAU
    #mpc     =   constants.pc *1e3   # [cm] MilliParsec -> from IAU
    Au      =   constants.Au        # [cm] Astronomical unit -> from IAU
    ly      =   constants.ly        # [cm] Light year -> from IAU
    Msol    =   constants.Msol      # [g] Solar mass -> from IAU
    Mearth  =   constants.Mearth    # [g]  Earth mass -> from IAU
    Mjupiter=   constants.Mjupiter  # [g]  Jupiter -> from IAU
    #Gyr     =   constants.yr /1e9   # [s]  GigaYear -> from IAU
    #Myr     =   constants.yr /1e6   # [s]  MegaYear -> from IAU
    yr      =   constants.yr        # [s]  Year -> from IAU
    X_frac  =   0.76                # Hydrogen fraction by mass -> cooling_module.f90 RAMSES

    scale.Mpc       = unit_l / pc / 1e6
    scale.kpc       = unit_l / pc / 1e3
    scale.pc        = unit_l / pc
    scale.mpc       = unit_l / pc * 1e3
    scale.ly        = unit_l / ly
    scale.Au        = unit_l / Au
    scale.km        = unit_l / 1.0e5
    scale.m         = unit_l / 1.0e2
    scale.cm        = unit_l
    scale.mm        = unit_l * 10.
    scale.μm        = unit_l * 1e4

    scale.Mpc3       = scale.Mpc^3
    scale.kpc3       = scale.kpc^3
    scale.pc3        = scale.pc^3
    scale.mpc3       = scale.mpc^3
    scale.ly3        = scale.ly^3
    scale.Au3        = scale.Au^3
    scale.km3        = scale.km^3
    scale.m3         = scale.m^3
    scale.cm3        = scale.cm^3
    scale.mm3        = scale.mm^3
    scale.μm3        = scale.μm^3

    scale.Msol_pc3  = unit_d * pc^3 / Msol
    scale.g_cm3     = unit_d

    scale.Msol_pc2  = unit_d * unit_l * pc^2 / Msol
    scale.g_cm2     = unit_d * unit_l

    scale.Gyr       = unit_t / yr / 1e9
    scale.Myr       = unit_t / yr / 1e6
    scale.yr        = unit_t / yr
    scale.s         = unit_t
    scale.ms        = unit_t * 1e3

    scale.Msol      = unit_d * unit_l^3 / Msol
    scale.Mearth    = unit_d * unit_l^3 / Mearth
    scale.Mjupiter  = unit_d * unit_l^3 / Mjupiter
    scale.g         = unit_d * unit_l^3
    scale.km_s      = unit_l / unit_t / 1e5
    scale.m_s       = unit_l / unit_t / 1e2
    scale.cm_s      = unit_l / unit_t

    scale.nH        = X_frac / mH * unit_d  # Hydrogen number density in [H/cc]
    scale.erg       = unit_m * (unit_l / unit_t)^2 # [g (cm/s)^2]
    scale.g_cms2    = unit_m / (unit_l * unit_t^2)

    scale.T_mu      = mH / kB * (unit_l / unit_t)^2 # T/mu [Kelvin]
    scale.Ba        = unit_m / unit_l / unit_t^2 # Barye (pressure) [cm-1 g s-2]

    return scale
end

"""
### Get a list of all exported Mera types and functions:
```julia
function viewmodule(modulename::Module)
```
"""
function viewmodule(modulename::Module)
    println()
    printstyled("[Mera]: Get a list of all exported Mera types and functions:\n", bold=true, color=:normal)
    printstyled("===============================================================\n", bold=true, color=:normal)
    module_list = names(modulename, all=false,imported= true)
    show(IOContext(stdout), "text/plain", module_list )
    return module_list
end

"""
### Convert a value to human-readable astrophysical units and round to ndigits
(pass the value in code units and the quantity specification (length, time) )
```julia
function humanize(value::Float64, scale::ScalesType, ndigits::Int, quantity::String)

return value, value_unit
```
"""
function humanize(value::Float64, scale::ScalesType, ndigits::Int, quantity::String)

    if quantity == ""
        round(value, digits=ndigits)

    elseif value == 0
        value_buffer = 0.
        value_unit = "x"
        return round(value_buffer, digits=ndigits), value_unit
    else

        if quantity == "length"
            sign_buffer = sign(value)
            value_buffer = value * scale.Mpc * sign_buffer
            value_unit = "Mpc"
            if value_buffer <= 1.
                value_buffer = value * scale.kpc * sign_buffer
                value_unit = "kpc"
                if value_buffer <= 1.
                    value_buffer = value * scale.pc * sign_buffer
                    value_unit = "pc"
                    if value_buffer <= 1.
                        value_buffer = value * scale.mpc * sign_buffer
                        value_unit = "mpc"
                        #if value_buffer < 1. #todo check
                        #    value_buffer = value * scale.au
                        #    value_unit = "au"
                            if value_buffer <= .1
                                value_buffer = value * scale.cm * sign_buffer
                                value_unit = "cm"
                            if value_buffer <= .1
                                value_buffer = value * scale.μm * sign_buffer
                                value_unit = "μm"
                                end
                            end
                        #end
                    end
                end
            end
            value_buffer = value_buffer * sign_buffer
        end

        if quantity == "time"
            sign_buffer = sign(value)
            value_buffer = value * scale.Gyr * sign_buffer
            value_unit = "Gyr"
            if value_buffer <= 1.
                value_buffer = value * scale.Myr * sign_buffer
                value_unit = "Myr"
                if value_buffer <= .1
                    value_buffer = value * scale.yr * sign_buffer
                    value_unit = "yr"
                    if value_buffer <= 1.
                        value_buffer = value * scale.s * sign_buffer
                        value_unit = "s"
                        if value_buffer <= 1.
                            value_buffer = value * scale.ms * sign_buffer
                            value_unit = "ms"
                        end
                    end
                end
            end

        value_buffer = value_buffer * sign_buffer

        end

        return round(value_buffer, digits=ndigits), value_unit
    end

end

function humanize(value::Float64, ndigits::Int, quantity::String)

    if quantity == ""
        round(value, digits=ndigits)

    else
        if quantity == "memory"
            value_buffer = value
            value_unit = "Bytes"
            if value_buffer > 1000.
                value_buffer = value_buffer / 1024.
                value_unit = "KB"
                if value_buffer > 1000.
                    value_buffer = value_buffer / 1024.
                    value_unit = "MB"
                    if value_buffer > 1000.
                        value_buffer = value_buffer / 1024.
                        value_unit = "GB"
                        if value_buffer > 1000.
                            value_buffer = value_buffer / 1024.
                            value_unit = "TB"
                        end
                    end
                end
            end
        end

        return round(value_buffer, digits=ndigits), value_unit
    end
end

#todo: define file type?
function skiplines(file, nlines::Int)
    for i=1:nlines
        read(file)
    end
    return
end

function getunit(dataobject, quantity::Symbol, vars::Array{Symbol,1}, units::Array{Symbol,1}; uname::Bool=false)
    idx = findall(x->x==quantity, vars)
    if length(idx) >= 1
        idx = idx[1]
        if  length(units) >= idx
            unit = units[idx]
        else
            unit = :standard
        end
    else
        unit = :standard
    end

    if unit == :standard
        if uname == false
            return 1.
        else
            return 1., unit
        end
    else
        if uname == false
            return getfield(dataobject.info.scale, unit)
        else
            return getfield(dataobject.info.scale, unit), unit
        end
    end

end

function getunit(dataobject::InfoType, unit::Symbol; uname::Bool=false)
    if unit == :standard
        if uname == false
            return 1.
        else
            return 1., unit
        end
    else
        if uname == false
            return getfield(dataobject.scale, unit)
        else
            return getfield(dataobject.scale, unit), unit
        end
    end
end

"""
### Create a New DataSetType from a Filtered Data Table

```julia
function construct_datatype(data::JuliaDB.AbstractIndexedTable, dataobject::HydroDataType)
return HydroDataType

function construct_datatype(data::JuliaDB.AbstractIndexedTable, dataobject::PartDataType)
return PartDataType

function construct_datatype(data::JuliaDB.AbstractIndexedTable, dataobject::ClumpDataType)
return ClumpDataType
```

### Example
```julia
# read simulation information
julia> info = getinfo(420)
julia> gas = gethydro(info)

# filter and create a new` data table
julia> density = 3. /gas.scale.Msol_pc3
julia> filtered_db = @filter gas.data :rho >= density

# construct a new HydroDataType
# (comparable to the object "gas" but only with filtered data)
julia> gas_new = construct_datatype(filtered_db, gas)
```
"""
function construct_datatype(data::JuliaDB.AbstractIndexedTable, dataobject::HydroDataType)
    hydrodata = HydroDataType()
    hydrodata.data = data
    hydrodata.info = dataobject.info
    hydrodata.lmin = dataobject.lmin
    hydrodata.lmax = dataobject.lmax
    hydrodata.boxlen = dataobject.boxlen
    hydrodata.ranges = dataobject.ranges
    hydrodata.selected_hydrovars = dataobject.selected_hydrovars
    hydrodata.used_descriptors = dataobject.used_descriptors
    hydrodata.smallr = dataobject.smallr
    hydrodata.smallc = dataobject.smallc
    hydrodata.scale = dataobject.scale
    return hydrodata
end

function construct_datatype(data::JuliaDB.AbstractIndexedTable, dataobject::PartDataType)
    partdata = PartDataType()
    partdata.data = data
    partdata.info = dataobject.info
    partdata.lmin = dataobject.lmin
    partdata.lmax = dataobject.lmax
    partdata.boxlen = dataobject.boxlen
    partdata.ranges = dataobject.ranges
    partdata.selected_partvars = dataobject.selected_partvars
    partdata.used_descriptors = dataobject.used_descriptors
    partdata.scale = dataobject.scale
    return partdata
end

# function construct_datatype(data::JuliaDB.AbstractIndexedTable, dataobject::GravDataType)
#     gravitydata = GravDataType()
#     gravitydata.data = data
#     gravitydata.info = dataobject.info
#     gravitydata.lmin = dataobject.lmin
#     gravitydata.lmax = dataobject.lmax
#     gravitydata.boxlen = dataobject.boxlen
#     gravitydata.ranges = dataobject.ranges
#     gravitydata.selected_gravvars = dataobject.selected_gravvars
#     gravitydata.scale = dataobject.scale
#     return gravitydata
# end

function construct_datatype(data::JuliaDB.AbstractIndexedTable, dataobject::ClumpDataType)
    clumpdata = ClumpDataType()
    clumpdata.data = data
    clumpdata.info = dataobject.info
    clumpdata.boxlen = dataobject.boxlen
    clumpdata.ranges = dataobject.ranges
    clumpdata.selected_clumpvars = dataobject.selected_clumpvars
    clumpdata.used_descriptors = dataobject.used_descriptors
    clumpdata.scale = dataobject.scale
    return clumpdata

end


text:
function ueberpruefeEindeutigkeit(f1,f2,f3,f4,f5,f6,f7,f8,f9,B,m_neu)
mneu = m_neu';
f = [f1; f2; f3; f4; f5; f6; f7; f8; f9];
fhorzvert = ones(size(B,2),1);
mhorzvert = 3*ones(size(B,1),1);
mMatrix = B*f;
e1 = [1;0;0;0;0;0;0;0;0];
e3 = [0;0;1;0;0;0;0;0;0];
e7 = [0;0;0;0;0;0;1;0;0];
e9 = [0;0;0;0;0;0;0;0;1];

display0 = 1;
display1 = 1;
display2 = 1;
display3 = 1;
display4 = 1;
display5 = 1;

for Zeile = 1:size(B,1)
    if sum(B[Zeile,:]) > sqrt(2)*3 || sum(B[Zeile,:]) <= 1 #Strahlen kürzer als 1cm oder länger als 3* Wurzel(2)
        if display0 == 1
            println("Schaut euch erneut an durch welche Pixel eure Strahlen wie lange laufen und überprüft,")
            println("ob ihr die Längen den richtigen Matrixeinträgen zugeordent habt.")
            display0 = display0+1;
        end
    else
        counter = 1;
        Nulleintraege = 1;
        for Spalte = 1:size(B,2)
            if B[Zeile,Spalte] == 0 #für nur einen Eintrag, dann 1,3,7,9
                Nulleintraege = Nulleintraege+1;
            end
            if B[Zeile,Spalte] < 0 || B[Zeile,Spalte] > 0 #für max 5 Einträge pro Zeile
                counter = counter+1;
            end
        end
        if counter >= 7 #Mehr als 5 Einträge pro Zeile
            if display1 == 1 && display0 == 1
                println("Schaut euch erneut an durch welche Pixel eure Strahlen wie lange laufen und überprüft,")
                println("ob ihr die Längen den richtigen Matrixeinträgen zugeordent habt.")
                display1 = display1+1;
            end
        else
            if B*fhorzvert == mhorzvert #keine neuen Strahlen verwendet
                if display2 == 1 && display1 == 1 && display0 == 1
                    println("Verwendet neue Strahlen und nicht noch einmal die Horizontalen oder Vertikalen, da so keine zusätlichen")
                    println("Informationen weiter gegeben werden.")
                    display2 = display2+1;
                end
            else
                if Nulleintraege-1 >= size(B,2)-1 #nur ein Eintrag falls 1,3,7,9
                    if sum(B[Zeile,:]) > B[Zeile,:]*e1 || sum(B[Zeile,:]) < B[Zeile,:]*e1 || sum(B[Zeile,:]) > B[Zeile,:]*e3 || sum(B[Zeile,:]) < B[Zeile,:]*e3 || sum(B[Zeile,:]) > B[Zeile,:]*e7 || sum(B[Zeile,:]) < B[Zeile,:]*e7 || sum(B[Zeile,:]) > B[Zeile,:]*e9 || sum(B[Zeile,:]) < B[Zeile,:]*e9 && display4 == 1
                        println("Schaut euch erneut an durch welche Pixel eure Strahlen wie lange laufen und überprüft,")
                        println("ob ihr die Pixel korret nummeriert und die richtigen Werte in die Matrix eingetragen habt.")
                    end
                end
            end
        end
    end
end
if display4 == 1 && display3 == 1 && display2 == 1 && display1 == 1 && display0 == 1 && round(mMatrix[1];digits=2) == round(mneu[1];digits=2) && round(mMatrix[2];digits=2) == round(mneu[2];digits=2) &&  round(mMatrix[3];digits=2) == round(mneu[3];digits=2) && round(mMatrix[4];digits=2) == round(mneu[4];digits=2) && round(mMatrix[5];digits=2) == round(mneu[5];digits=2) && round(mMatrix[6];digits=2) == round(mneu[6];digits=2) && round(mMatrix[7];digits=2) == round(mneu[7];digits=2) && round(mMatrix[8];digits=2) == round(mneu[8];digits=2) && round(mMatrix[9];digits=2) == round(mneu[9];digits=2)
    println("Eure Idee die Anzahl der Strahlen zu vergrößern ist korrekt! Leider ist diese Lösung noch immer nicht eindeutig. Es liegen noch nicht ausreichend viele Informationen vor, so dass weitere Strahlen hizugefügt werden müssen. Das heißt, dass so die Struktur des Objekts weiterhin nicht bestimmen kann.
Das hinzufügen von immer mehr Strahlen soll automatisiert werden, weshalb wir uns im nächsten Schritt überlegen müssen, wie die Strahlen abhängig vom Startpunkt und einem Winkel beschrieben werden können. So kann im Anschluss geschaut werden, dass ausreichend viele Strahlen betrachtet werden, um die Struktur des Objekts eindeutig zu bestimmen.")
elseif display4 == 1 && display3 == 1 && display2 == 1 && display1 == 1 && display0 == 1 && mMatrix[1] == mneu[1]
    println("Überprüft eure Berechnung von m erneut. Beachtet, dass auch hier die Länge der Pixel berücksichtigt werden müssen.")
end

end


text:
"""
    module CPython

This module provides a direct interface to the Python C API.
"""
module C

import Base: @kwdef
import CondaPkg
using Libdl, Requires, UnsafePointers, Serialization, ..Utils

include("consts.jl")
include("pointers.jl")
include("extras.jl")
include("context.jl")
include("gil.jl")
include("jlwrap.jl")

function __init__()
    init_context()
    with_gil() do
        init_jlwrap()
    end
end

end


text:
"""
    Web3

A module for Ethereum connectivity.

Parse ABI files, encode/decode ABI data, and make JSON-RPC calls.
"""
module Web3

#=
Parse ABI files and encode / decode ABI data

TUPLE AND FIXED LENGTH ARRAY: [fix1 ... fixN][var1 ... varN]
for static types, fixN is the encoded data and varN is empty
for dynamic types, fixN is the length as a uint256 and varN is the encoded data

DYNAMIC ARRAY: [LEN][ARRAY ENCODING...]

UINT<M> and INT<M>: encoded as 32-byte, big-endian values, i.e. as uint256 and int256 numbers

BYTES: [len as uint256][data]

STRING: a bytes representation of the UTF-8 data

ADDRESS: a uint160, encoded as a 32-byte, big-endian value

BOOL: a uint8, encoded as a 32-byte, big-endian value

BYTES<M>: the sequence of bytes, padded with trailing zeroes to a length of 32-bytes

FUNCTION CALL: [4-byte hash][args tuple]
hash is the first 4 high-order bytes of the hash of the signature name(type,type...)

RETURN VALUE: [return tuple]

EVENT: [contract address][topics[0]: signature][topics[n]: indexed args[n - 1]][data -- unindexed args]
topics are 32-bytes; dynamic values are represented as a hash and lose information
=#

####
#### CAVEATS
####
#### 1) THIS USES THE JSON API FOR KECCAK, WHICH IS HORRIBLY INEFFICIENT
####    THIS IS ONLY DONE WHILE PARSING ABIS BUT WE SHOULD CHANGE IT TO A LOCAL CALL
####    MAYBE LINK TO THE CORRUS VERSION: https://github.com/coruus/keccak-tiny
####    THE PROBLEM IS THAT IT RELIES ON MEMSET_S WHICH IS NOT CURRENTLY AVAILABLE IN LINUX
####    SURE WOULD BE NICE TO HAVE SHAKE256 IN PURE JULIA
####    MAYBE SOMEONE WILL PORT [jsSHA](https://caligatio.github.io/jsSHA)
####

using HTTP, JSON, Core

export readABI, Web3Connection, Contract, contracts
export FunctionCall, encodefunctioncall, decodefunctioncall
export FunctionResult, encodefunctionresult, decodefunctionresult
export Event, encodeevent, decodeevent
export clientversion, eth, utils

####################
# Web3
####################

"""
    Web3Connection

A JSON-RPC connection to an Ethereum node

```
    web3 = Web3Connection("http://localhost:8545")
    clientversion(web3) # -- return the client version
    jsonget("http://localhost:8545", :web3_clientVersion, []) # -- equivalent to clientversion(web3)
    rawjsonget("http://localhost:8545", :web3_clientVersion, []) # -- return full JSON object
```
"""
struct Web3Connection
    url
end

"""
    jsonget(url, methodname, params...)

Call a JSON-RPC method and return the result property of the JSON result
"""
jsonget(url, method, params...) = rawjsonget(url, method, params...)["result"]

"""
    rawjsonget(url, method, params...)

Call a JSON-RPC method and return the full JSON result
"""
function rawjsonget(url, method, params...)
    req = JSON.json(Dict([
        :jsonrpc => "2.0"
        :method => method
        :params => params
    ]))
    if verbose
        println("REQUEST: $(repr(req))")
    end
    resp = HTTP.request("POST", url, [], req)
    if resp.status == 200
        resultstr = String(resp.body)
        result = JSON.parse(resultstr)
        if verbose
            println("RESULT JSON: $resultstr")
            println("RESULT: $(repr(result))")
        end
        result
    else
        throw(ErrorException("Error, status = $(resp.status)"))
    end
end

"""
    apifunc

Create a call to an API func, given a function that converts inputs to JSON-ready inputs
"""
apifunc(apimethod, func) = (con::Web3Connection, args...; raw=false)-> (raw ? rawjsonget : jsonget)(con.url, apimethod, func(args...)...)

function hash end

const clientversion = apifunc(:web3_clientVersion, ()->[])
const eth = (
    gettransactioncount = apifunc(:eth_getTransactionCount, (addr, ctx)-> (addr, ctx)),
    gettransactionbyhash = apifunc(:eth_getTransactionByHash, (hash)-> (hash)))
const utils = (
    sha3 = hash,
    keccak = hash)

# This is a separate function so that test code can override it
function hash(con::Web3Connection, str::String)
    hex2bytes(jsonget(con.url, :web3_sha3, ("0x" * bytes2hex(Vector{UInt8}(str))))[3:end])
end

"""
    resultbytes

Convert a JSON-RPC API call result to bytes
"""
resultbytes(func) = (args...)-> hex2bytes(func(args...)[3:end])

####################
# ABI
####################

struct ABIContext
    contract
    connection
end

struct Int256
    big::Int128
    little::UInt128
    Int256(i::Integer) = new(i < 0 ? -1 : 0, i)
    Int256(big::Integer, little::Unsigned) = new(big, little)
end

struct UInt256
    big::UInt128
    little::UInt128
    UInt256(i::Unsigned) = new(0, i)
    UInt256(big::Unsigned, little::Unsigned) = new(big, little)
end

struct ABIType{T}
    ABIType(arg) = new{arg}()
end

struct Decl{ENCMODE, X, Y, Z}
    name::String
    typename::String
    components::Array
    indexed::Bool
end

struct ABIFunction
    constant
    hash
    inputs::Array{Decl}
    name
    outputs
    payable
    signature
    statemutability
end

struct ABIEvent
    name
    hash
    signature
    anonymous
    inputs::Array{Decl}
end

struct FunctionCall
    decl
    inputs::Array
end

struct FunctionResult
    decl
    result
end

struct Event
    decl
    parameters::Array
end

struct Contract
    functions
    events
    Contract() = new(Dict(), Dict())
end

const NumDecl = Union{Decl{T, :int} where T, Decl{T, :uint} where T}
const FunctionABI = Union{ABIType{:function}, ABIType{:constructor}, ABIType{:fallback}};

##############
# VARS
##############

"A dictionary of contract-address => Contract structures"
contracts = Dict()

fixedarraypattern = r".*\[([^\]]+)\]"
bitspattern = r"^([^[0-9]+)([0-9]*)"

##############
# ENCODING
##############

"""
    encodefunctioncall(io::IO, f::ABIFunction, inputs::Array)

Encode a call to a function
"""
function encodefunctioncall(io::IO, f::ABIFunction, inputs::Array)
    basicencodefunctioncall(io, f, inputs)
end

function encodefunctioncall(io::IOBuffer, f::ABIFunction, inputs::Array)
    basicencodefunctioncall(io, f, inputs)
    io.data
end

function basicencodefunctioncall(io::IO, f::ABIFunction, inputs::Array)
    write(io, f.hash)
    encode(io, f.inputs, inputs)
end

"""
    encodefunctionresult(io::IO, f::ABIFunction, outputs::Array)

Encode the results of a function
"""
function encodefunctionresult(io::IO, f::ABIFunction, outputs::Array)
    basicencodefunctionresult(io, f, outputs)
end

function encodefunctionresult(io::IOBuffer, f::ABIFunction, outputs::Array)
    basicencodefunctionresult(io, f, outputs)
    io.data
end

function basicencodefunctionresult(io::IO, f::ABIFunction, outputs::Array)
    write(io, f.hash)
    encode(io, f.outputs, outputs)
end

"""
    encodeevent(io::IO, e::ABIEvent, inputs::Array)

Encode an event
"""
function encodeevent(io::IO, e::ABIEvent, inputs::Array)
    basicencodeevent(io, e, inputs)
end

function encodeevent(io::IOBuffer, e::ABIEvent, inputs::Array)
    basicencodeevent(io, e, inputs)
    io.data
end

function basicencodeevent(io::IO, e::ABIEvent, inputs::Array)
    write(io, e.hash)
    encode(io, e.inputs, inputs)
end

function encode(io::IO, decl::Union{Decl, Array}, value)
    encodehead(io, decl, value)
    encodetail(io, decl, value)
end

# utilities
writeint(io, i::Signed) = writeints(io, Int128(i < 0 ? -1 : 0), Int128(i))
writeint(io, i::Unsigned) = writeints(io, UInt128(0), UInt128(i))
writeint(io, i::Union{Int256, UInt256}) = writeints(io, i.big, i.little)
writeints(io, ints...) = write(io, hton.(ints)...)
encodescalar(io::IO, decl::Decl{Any, :bool}, value) = writeint(io, value ? 1 : 0)
encodescalar(io::IO, decl::NumDecl, value) = writeint(io, value)

# scalar types
encodehead(io, decl::Decl{:scalar}, v) = encodescalar(io, decl, v)
encodetail(io, ::Decl{:scalar}, v) = nothing

# dynamic types
encodehead(io, ::Union{Decl{:string}, Decl{:bytes}, Decl{:dynamic}}, v) = writeuint(io, length(v))
encodetail(io, ::Union{Decl{:string}, Decl{:bytes}}, v) = write(io, v)
function encodetail(io::IO, decl::Decl{:dynamic}, values)
    for v in values
        encodescalar(decl, v)
    end
end

# array
function encodehead(io::IO, decl::Decl{:array, BASE, BITS, LENGTH}, values) where {BASE, BITS, LENGTH}
    t = arraycomptype(decl)
    for i in 1:LENGTH
        encodehead(io, t, values[i])
    end
end
function encodetail(io::IO, decl::Decl{:array}, values)
    t = arraycomptype(decl)
    for i in 1:LENGTH
        encodetail(io, t, values[i])
    end
end
function arraycomptype(decl::Decl{:array, BASE, BITS, LENGTH}) where {BASE, BITS, LENGTH}
    if decl.components != nothing
        Decl{:tuple, length(decl.components), :n, :n}(decl.name, "tuple", decl.components, false)
    else
        Decl{:scalar, BASE, BITS, :none}(decl.name, decl.typename, :nothing, false)
    end
end

# tuple
encodehead(io::IO, decl::Decl{:tuple}, values) = encodehead(io, decl.components, values)
encodetail(io::IO, decl::Decl{:tuple}, values) = encodetail(io, decl.components, values)
function encodehead(io::IO, decls::Array, values)
    for i in 1:length(decls)
        encodehead(io, decls[i], values[i])
    end
end
function encodetail(io::IO, decls::Array, values)
    for i in 1:length(decls)
        encodetail(io, decls[i], values[i])
    end
end

############
# DECODING
############

signedTypes = Dict([t.size => t for t in (Int8, Int16, Int32, Int64, Int128)])
unsignedTypes = Dict([t.size => t for t in (UInt8, UInt16, UInt32, UInt64, UInt128)])

function readint(io::IO)
    big = read(io, Int128)
    Int256(ntoh(big), ntoh(read(io, UInt128)))
end

function readuint(io::IO)
    big = read(io, UInt128)
    Int256(ntoh(big), ntoh(read(io, UInt128)))
end

readlength(io::IO) = (read(io, UInt128);read(io, UInt128))

"""
    decodefunctioncall(io::IO, con::Contract)

Decode a function call
"""
function decodefunctioncall(io::IO, con::Contract)
    decl = con.functions[read(io, 4)]
    FunctionCall(decl, decode(io, decl.inputs))
end

"""
    decodefunctionresult(io::IO, con::Contract)

Decode a function call result
"""
function decodefunctionresult(io::IO, con::Contract)
    hash = read(io, 4)
    f = con.functions[hash]
    FunctionResult(f, decode(io, f.outputs))
end

"""
    decodeevent(io::IO, con::Contract)

Decode an event in a transaction log
"""
function decodeevent(io::IO, con::Contract)
    decl = con.events[read(io, 4)]
    Event(decl, decode(io, decl.inputs))
end

# general
decode(io::IO, decl) = decodetail(io, decl, decodehead(io, decl))

# scalar types
decodehead(io::IO, decl::Decl{:scalar}) = decodescalar(io, decl)
decodetail(io::IO, ::Decl{:scalar}, head) = head
function decodescalar(io::IO, ::Decl{T, :int, SIZE}) where {T, SIZE}
    big = ntoh(read(io, Int128))
    little = ntoh(read(io, UInt128))
    SIZE <= 128 ? smallint(SIZE, big < 0 ? -little : little) : Int256(big, little)
end
function decodescalar(io::IO, ::Decl{T, :uint, SIZE}) where {T, SIZE}
    big = ntoh(read(io, UInt128))
    little = ntoh(read(io, UInt128))
    SIZE <= 128 ? smalluint(SIZE, little) : UInt256(big, little)
end

smallint(size, value) = signedTypes[Int(2^ceil(log2(floor((33 + 7) / 8))))](value)
smalulint(size, value) = unsignedTypes[Int(2^ceil(log2(floor((33 + 7) / 8))))](value)

function decode(io::IO, ::Decl{T, :bool}) where T
    big = read(io, UInt128)
    little = read(io, UInt128)
    litle == 0 ? false : true
end

# dynamic types
decodehead(io::IO, decl::Decl{:dynamic, :bytes}) = readlength(io)
decodetail(io::IO, decl::Decl{:dynamic, :bytes}, head) = readbytes(io, head)
readbytes(io::IO, len) = read(io, len)

decodehead(io::IO, decl::Decl{:dynamic, :string}) = readlength(io)
decodetail(io::IO, decl::Decl{:dynamic, :string}, head) = readbytes(io, head)
readstring(io::IO, len) = String(read(io, len))

# array types
function decode(io::IO, decl::Decl{:array, BASE, BITS, LENGTH}) where {BASE, BITS, LENGTH}
    [decodescalar(io, decl) for i in 1:LENGTH]
end
function decode(io::IO, ::Decl{:array, :bytes, BITS, LENGTH}) where {BITS, LENGTH}
    lens = [readlength(io) for i in 1:LENGTH]
    [readbytes(io, len) for len in lens]
end
function decode(io::IO, ::Decl{:array, :string, BITS, LENGTH}) where {BITS, LENGTH}
    lens = [readlength(io) for i in 1:LENGTH]
    [String(readbytes(io, len)) for len in lens]
end
function decode(io::IO, decl::Decl{:array, :tuple, BITS, LENGTH}) where {BITS, LENGTH}
    decodetail(io::IO, decl.components, decodehead(io::IO, decl.components))
end

# tuple
decode(io::IO, decl::Decl{:tuple}) = decodetail(io, decl.components)
decode(io::IO, decls::Array) = decodetail(io, decls, decodehead(io, decls))
decodehead(io::IO, decls::Array) = [decodehead(io, head) for head in decls]
decodetail(io::IO, decls::Array, heads) = [decodetail(io, decls[i], heads[i]) for i in 1:length(decls)]

####################
# DECL PARSING
####################

parseABI(con::ABIContext, json) = parseABI(con, ABIType(Symbol(get(json, "type", "function"))), json)

function parseABI(con::ABIContext, ::FunctionABI, func)
    name = func["name"]
    args = join((arg-> arg["type"]).(func["inputs"]), ",")
    sig = "$name($args)"
    ABIFunction(
        func["constant"],
        utils.sha3(con.connection, sig)[1:4],
        parseargs(func["inputs"]),
        name,
        haskey(func, "outputs") ? parseargs(func["outputs"]) : [],
        func["payable"],
        sig,
        func["stateMutability"]
    )
end

function parseABI(con::ABIContext, ::ABIType{:event}, evt)
    name = evt["name"]
    args = join((arg-> arg["type"]).(evt["inputs"]), ",")
    sig = "$name($args)"
    ABIEvent(evt["name"], utils.sha3(con.connection, sig)[1:4], sig, get(evt, "anonymous", false), parseargs(evt["inputs"]))
end

"""
    readABI(con::Web3Connection, contractname::String, stream::IO)

Read an ABI file for a contract
"""
function readABI(con::Union{Web3Connection, Dict}, contractname::String, stream::IO)
    contract = Contract()
    context = ABIContext(contract, con)
    d = JSON.parse(stream)
    close(stream)
    for json in d
        obj = parseABI(context, json)
        if isa(obj, ABIFunction)
            contract.functions[obj.name] = contract.functions[obj.hash] = obj
            if verbose
                println("$(bytes2hex(obj.hash)) $(obj.signature) $(repr(obj)) $(bytes2hex(utils.sha3(con, obj.signature)))")
            end
        elseif isa(obj, ABIEvent)
            contract.events[obj.name] = contract.events[obj.hash] = obj
            if verbose
                println("$(bytes2hex(obj.hash)) $(obj.signature) $(repr(obj)) $(bytes2hex(utils.sha3(con, obj.signature)))")
            end
        elseif verbose
            println(repr(obj))
        end
    end
    contracts[contractname] = contract
end

function basetypefor(typename)
    m = match(bitspattern, typename)
    (Symbol(m[1]), m[2] == "" ? 256 : parse(Int, m[2]))
end

function typefor(typename, arg)
    atype = arg["type"]
    if endswith(atype, "[]")
        (:dynamic, basetypefor(typename)..., :none)
    elseif endswith(atype, "]")
        m = match(fixedarraypattern, typename)
        (:array, basetypefor(typename)..., parse(Int, m[1]))
    elseif atype == "tuple"
        (:tuple, length(arg["components"]), :none, :none)
    elseif atype in ["string", "bytes"]
        (Symbol(atype), :none, :none, :none)
    else
        (:scalar, basetypefor(typename)..., :none)
    end
end

function parsearg(arg)
    typename = arg["type"]
    (enctype, atype, bits, len) = typefor(typename, arg)
    Decl{enctype, atype, bits, len}(
        arg["name"],
        typename,
        haskey(arg, "components") ? parseargs(arg["components"]) : [],
        get(arg, "indexed", false)
    )
end

parseargs(args) = parsearg.(args)

####################
# UTILS
####################

global verbose = false

function setverbose(v)
    global verbose
    verbose = v
end

end # module


text:
# # Welford's algorithm
#
# * Welford, B. P. (1962). "Note on a method for calculating corrected sums of squares and products". Technometrics. 4 (3): 419–420. doi:10.2307/1266577. JSTOR 1266577
# * [Welford's online algorithm - Algorithms for calculating variance - Wikipedia](https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm)

using Transducers

using Transducers: R_, start, next, complete, inner, xform, wrap, unwrap, wrapping

using Measurements

export Mean

struct Mean <: Transducer
end

function Transducers.start(rf::R_{Mean}, result)
    # @info "Transducers.start" result
    private_state = (0, 0.0, 0.0)
    return wrap(rf, private_state, start(inner(rf), result))
end

function Transducers.next(rf::R_{Mean}, result, input)
    # @info "Transducers.next" result input
    wrapping(rf, result) do st, iresult
        (n, μ, M2) = st
        n += 1
        δ = input - μ
        μ += δ/n
        δ2 = input - μ
        M2 += δ*δ2
        stderr = sqrt(M2)/n
        iinput = μ ± stderr
        iresult = next(inner(rf), iresult, iinput)
        return (n, μ, M2), iresult
    end
end

function Transducers.complete(rf::R_{Mean}, result)
    # @info "Transducers.complete" result
    _private_state, inner_result = unwrap(rf, result)
    return complete(inner(rf), inner_result)
end


text:
module MsgPackRpcServer

using MsgPack

run(port,m) = begin
  server = listen(port)
  while true
    conn = accept(server)
    line = ""
  #  @async begin
      try
        while true
          line = string(line, MsgPack.unpack(conn))
          if sizeof(line)>0
            args = eval(parse(line))
            evaluated = eval(m, Expr(:call, parse(args[3]), args[4]...))
            #println(evaluated)
            write(conn, MsgPack.pack([1,args[2],nothing,evaluated]))
            line = ""
          end
        end
      catch err
        #print("connection ended with error $err")
      end
  #  end
    close(conn)
  end
end

end # module

max_stars_repo_name: remore/MsgPackRpcServer

max_issues_repo_name: remore/MsgPackRpcServer

max_forks_repo_name: remore/MsgPackRpcServer


text:
using Memento

# Create our module level logger (this will get precompiled)
const LOGGER = getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via
# `getlogger(MyModule)`
# NOTE: If this line is not included then the precompiled `MyModule.LOGGER` won't be
# registered at runtime.
__init__() = Memento.register(LOGGER)

function setloglevel!(level::String)
    # Options correspond to Memento.jl's levels.
    # https://invenia.github.io/Memento.jl/latest/man/intro.html#Logging-levels-1
    Memento.config!(level; recursive = true)
end

setloglevel!("notice")