text:
function addsetv!(s::AbstractVector{T}, v::T) where {T}
    for sᵢ ∈ s
        sᵢ === v && return nothing
    end
    push!(s, v)
    nothing
end
function mergesetv!(s1::AbstractVector{T}, s2::AbstractVector{T}) where {T}
    for s ∈ s2
        addsetv!(s1, s)
    end
    nothing
end
function mergesetdiffv!(
    s1::AbstractVector{T},
    s2::AbstractVector{T},
    s3::AbstractVector{T}
) where {T}
    for s ∈ s2
        s ∉ s3 && addsetv!(s1, s)
    end
    nothing
end
# Everything in arg2 (s1) that isn't in arg3 (s2) is added to arg1 (s3)
function setdiffv!(s3::AbstractVector{T}, s1::AbstractVector{T}, s2::AbstractVector{T}) where {T}
    for s ∈ s1
        (s ∈ s2) || (s ∉ s3 && push!(s3, s))
    end
end
function setdiffv!(s4::AbstractVector{T}, s3::AbstractVector{T}, s1::AbstractVector{T}, s2::AbstractVector{T}) where {T}
    for s ∈ s1
        (s ∈ s2) ? (s ∉ s4 && push!(s4, s)) : (s ∉ s3 && push!(s3, s))
    end
end
function update_deps!(deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, parent::Operation)
    mergesetv!(deps, loopdependencies(parent))#, reduceddependencies(parent))
    if !(isload(parent) || isconstant(parent)) && !isreductcombineinstr(parent)
        mergesetv!(reduceddeps, reduceddependencies(parent))
    end
    nothing
end

function pushparent!(parents::Vector{Operation}, deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, parent::Operation)
    push!(parents, parent)
    update_deps!(deps, reduceddeps, parent)
end
# function pushparent!(mpref::ArrayReferenceMetaPosition, parent::Operation)
#     pushparent!(mpref.parents, mpref.loopdependencies, mpref.reduceddeps, parent)
# end
function add_parent!(
    vparents::Vector{Operation}, deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, ls::LoopSet, var, elementbytes::Int, position::Int
)
    parent = if var isa Symbol
        # if var === :kern_1_1
        #     @show operations(ls) ls.preamble_symsym
        # end
        opp = getop(ls, var, elementbytes)
        # if var === :kern_1_1
        #     @show operations(ls) ls.preamble_symsym
        # end
        # @show var opp first(operations(ls)) opp === first(operations(ls))
        if iscompute(opp) && instruction(opp).instr === :identity && length(loopdependencies(opp)) < position && isone(length(parents(opp))) && name(opp) === name(first(parents(opp)))
            first(parents(opp))
        else
            opp
        end
    elseif var isa Expr #CSE candidate
        add_operation!(ls, gensym(:temporary), var, elementbytes, position)
    else # assumed constant
        add_constant!(ls, var, elementbytes)
        # add_constant!(ls, var, deps, gensym(:loopredefconst), elementbytes)
    end
    pushparent!(vparents, deps, reduceddeps, parent)
end
# function add_reduction!(
#     vparents::Vector{Operation}, deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, ls::LoopSet, var::Symbol, elementbytes::Int
# )
#     get!(ls.opdict, var) do
#         add_constant!(ls, var, elementbytes)
#     end
# end
function search_tree(opv::Vector{Operation}, var::Symbol) # relies on cycles being forbidden
    for opp ∈ opv
        name(opp) === var && return true
        search_tree(parents(opp), var) && return true
    end
    false
end
function search_tree(opv::Vector{Operation}, var::Operation) # relies on cycles being forbidden
    for opp ∈ opv
        opp === var && return true
        search_tree(parents(opp), var) && return true
    end
    false
end
function update_reduction_status!(parentvec::Vector{Operation}, deps::Vector{Symbol}, parent::Symbol)
    for opp ∈ parentvec
        if name(opp) === parent
            mergesetv!(reducedchildren(opp), deps)
            break
        elseif search_tree(parents(opp), parent)
            mergesetv!(reducedchildren(opp), deps)
            update_reduction_status!(parents(opp), deps, parent)
            break
        end
    end
end
# function add_compute!(ls::LoopSet, op::Operation)
    # @assert iscompute(op)
    # pushop!(ls, child, name(op))
# end
# function isreductzero(op::Operation, ls::LoopSet, reduct_zero::Symbol)
#     isconstant(op) || return false
#     reduct_zero === op.instruction.mod && return true
#     if reduct_zero === :zero
#         iszero(ls, op) && return true
#     elseif reduct_zero === :one
#         isone(ls, op) && return true
#     end
#     false
# end
function add_reduced_deps!(op::Operation, reduceddeps::Vector{Symbol})
    # op.dependencies = copy(loopdependencies(op))
    # mergesetv!(loopdependencies(op), reduceddeps)
    reduceddepsop = reduceddependencies(op)
    if reduceddepsop === NODEPENDENCY
        op.reduced_deps = copy(reduceddeps)
    else
        mergesetv!(reduceddepsop, reduceddeps)
    end
    # reduceddepsop = reducedchildren(op)
    # if reduceddepsop === NODEPENDENCY
    #     op.reduced_children = copy(reduceddeps)
    # else
    #     mergesetv!(reduceddepsop, reduceddeps)
    # end
    nothing
end

function substitute_op_in_parents!(
    vparents::Vector{Operation}, replacer::Operation, replacee::Operation, reduceddeps::Vector{Symbol}
)
    found = false
    for i ∈ eachindex(vparents)
        opp = vparents[i]
        if opp === replacee
            vparents[i] = replacer
            found = true
        else
            fopp = substitute_op_in_parents!(parents(opp), replacer, replacee, reduceddeps)
            fopp && add_reduced_deps!(opp, reduceddeps)
            found |= fopp
        end
    end
    found
end

function add_reduction_update_parent!(
    vparents::Vector{Operation}, deps::Vector{Symbol}, reduceddeps::Vector{Symbol}, ls::LoopSet,
    parent::Operation, instr::Instruction, reduction_ind::Int, elementbytes::Int
)
    var = name(parent)
    isouterreduction = parent.instruction === LOOPCONSTANT
    instrclass = reduction_instruction_class(instr) # key allows for faster lookups
    reduct_zero = reduction_zero(instrclass)
    # if parent is not an outer reduction...
    # if !isouterreduction && !isreductzero(parent, ls, reduct_zero)
    add_reduct_instruct = !isouterreduction && !isconstant(parent)
    if add_reduct_instruct
        # We add
        reductcombine = reduction_scalar_combine(instrclass)
        # reductcombine = :identity
        reductsym = gensym(:reduction)
        reductinit = add_constant!(ls, gensym(:reductzero), loopdependencies(parent), reductsym, elementbytes, :numericconstant)
        if reduct_zero === :zero
            push!(ls.preamble_zeros, (identifier(reductinit), IntOrFloat))
        else
            push!(ls.preamble_funcofeltypes, (identifier(reductinit), reduct_zero))
        end
    else
        reductinit = parent
        reductsym = var
        reductcombine = :identity#Symbol("")
    end
    combineddeps = copy(deps); mergesetv!(combineddeps, reduceddeps)
    # directdependency && pushparent!(vparents, deps, reduceddeps, reductinit)#parent) # deps and reduced deps will not be disjoint
    if reduction_ind > 0 # if is directdependency
        insert!(vparents, reduction_ind, reductinit)
        if instr.instr ∈ (:-, :vsub!, :vsub, :/, :vfdiv!, :vfidiv!)
            update_deps!(deps, reduceddeps, reductinit)#parent) # deps and reduced deps will not be disjoint
        end
    elseif !isouterreduction && reductinit !== parent
        substitute_op_in_parents!(vparents, reductinit, parent, reduceddeps)
    end
    update_reduction_status!(vparents, reduceddeps, name(reductinit))
    # this is the op added by add_compute
    op = Operation(length(operations(ls)), reductsym, elementbytes, instr, compute, deps, reduceddeps, vparents)
    parent.instruction === LOOPCONSTANT && push!(ls.outer_reductions, identifier(op))
    opout = pushop!(ls, op, var) # note this overwrites the entry in the operations dict, but not the vector
    # isouterreduction || iszero(length(reduceddeps)) && return opout
    # return opout
    isouterreduction && return opout
    # create child op, which is the reduction combination
    childrdeps = Symbol[]; childparents = Operation[ op ]#, parent ]
    add_reduct_instruct && push!(childparents, parent)
    childdeps = loopdependencies(reductinit)
    setdiffv!(childrdeps, loopdependencies(op), childdeps)
    child = Operation(
        length(operations(ls)), name(parent), elementbytes, reductcombine, compute, childdeps, childrdeps, childparents
    )
    pushop!(ls, child, name(parent))
    opout
end

function add_compute!(
    ls::LoopSet, var::Symbol, ex::Expr, elementbytes::Int, position::Int,
    mpref::Union{Nothing,ArrayReferenceMetaPosition} = nothing
)
    @assert ex.head === :call
    # instr = instruction(first(ex.args))::Symbol
    instr = instruction!(ls, first(ex.args))::Instruction
    args = @view(ex.args[2:end])
    (instr.instr === :(^) && length(args) == 2 && (args[2] isa Number)) && return add_pow!(ls, var, args[1], args[2], elementbytes, position)
    vparents = Operation[]
    deps = Symbol[]
    reduceddeps = Symbol[]
    reduction_ind = 0
    # @show ex first(operations(ls)) === getop(ls, :kern_1_1, elementbytes) first(operations(ls)) getop(ls, :kern_1_1, elementbytes)
    for (ind,arg) ∈ enumerate(args)
        if var === arg
            reduction_ind = ind
            # add_reduction!(vparents, deps, reduceddeps, ls, arg, elementbytes)
            getop(ls, arg, elementbytes)
        elseif arg isa Expr
            isref, argref = tryrefconvert(ls, arg, elementbytes, varname(mpref))
            if isref
                if mpref == argref
                    if varname(mpref) === var
                        id = findfirst(r -> r == mpref.mref, ls.refs_aliasing_syms)
                        mpref.varname = var = isnothing(id) ? var : ls.syms_aliasing_refs[id]
                        reduction_ind = ind
                        mergesetv!(deps, loopdependencies(add_load!(ls, argref, elementbytes)))
                    else
                        pushparent!(vparents, deps, reduceddeps, add_load!(ls, argref, elementbytes))
                    end
                else
                    argref.varname = gensym(:tempload)
                    pushparent!(vparents, deps, reduceddeps, add_load!(ls, argref, elementbytes))
                end
            else
                add_parent!(vparents, deps, reduceddeps, ls, arg, elementbytes, position)
            end
        elseif arg ∈ ls.loopsymbols
            loopsymop = add_loopvalue!(ls, arg, elementbytes)
            pushparent!(vparents, deps, reduceddeps, loopsymop)
        else
            add_parent!(vparents, deps, reduceddeps, ls, arg, elementbytes, position)
        end
    end
    reduction = reduction_ind > 0
    loopnestview = view(ls.loopsymbols, 1:position)
    if iszero(length(deps)) && reduction
        append!(deps, loopnestview)
        append!(reduceddeps, loopnestview)
    else
        newloopdeps = Symbol[]; newreduceddeps = Symbol[];
        setdiffv!(newloopdeps, newreduceddeps, deps, loopnestview)
        mergesetv!(newreduceddeps, reduceddeps)
        deps = newloopdeps; reduceddeps = newreduceddeps
    end
    op = if reduction || search_tree(vparents, var)
        parent = ls.opdict[var]
        setdiffv!(reduceddeps, deps, loopdependencies(parent))
        # parent = getop(ls, var, elementbytes)
        # if length(reduceddeps) == 0
        if all(!in(deps), reduceddeps)
            insert!(vparents, reduction_ind, parent)
            mergesetv!(deps, loopdependencies(parent))
            op = Operation(length(operations(ls)), var, elementbytes, instr, compute, deps, reduceddeps, vparents)
            pushop!(ls, op, var)
        else
            add_reduction_update_parent!(vparents, deps, reduceddeps, ls, parent, instr, reduction_ind, elementbytes)
        end
    else
        op = Operation(length(operations(ls)), var, elementbytes, instr, compute, deps, reduceddeps, vparents)
        pushop!(ls, op, var)
    end
    # maybe_const_compute!(ls, op, elementbytes, position)
    op
end

function add_compute!(
    ls::LoopSet, LHS::Symbol, instr, vparents::Vector{Operation}, elementbytes
)
    deps = Symbol[]
    reduceddeps = Symbol[]
    foreach(parent -> update_deps!(deps, reduceddeps, parent), vparents)
    op = Operation(length(operations(ls)), LHS, elementbytes, instr, compute, deps, reduceddeps, vparents)
    pushop!(ls, op, LHS)
end

# adds x ^ (p::Real)
function add_pow!(
    ls::LoopSet, var::Symbol, x, p::Real, elementbytes::Int, position::Int
)
    xop::Operation = if x isa Expr
        add_operation!(ls, gensym(:xpow), x, elementbytes, position)
    elseif x isa Symbol
        if x ∈ ls.loopsymbols
            add_loopvalue!(ls, x, elementbytes)
        else
            xo = get(ls.opdict, x, nothing)
            if isnothing(xo)
                pushpreamble!(ls, Expr(:(=), var, Expr(:call, :(^), x, p)))
                return add_constant!(ls, var, elementbytes)
            end
            xo
        end
    elseif x isa Number
        pushpreamble!(ls, Expr(:(=), var, x ^ p))
        return add_constant!(ls, var, elementbytes)
    end
    pint = round(Int, p)
    if p != pint
        pop = add_constant!(ls, p, elementbytes)
        return add_compute!(ls, var, :^, [xop, pop], elementbytes)
    end
    if pint == -1
        return add_compute!(ls, var, :inv, [xop], elementbytes)
    elseif pint < 0
        xop = add_compute!(ls, gensym(:inverse), :inv, [xop], elementbytes)
        pint = - pint
    end
    if pint == 0
        op = Operation(length(operations(ls)), var, elementbytes, LOOPCONSTANT, constant, NODEPENDENCY, Symbol[], NOPARENTS)
        push!(ls.preamble_funcofeltypes, (identifier(op),:one))
        return pushop!(ls, op)
    elseif pint == 1
        return add_compute!(ls, var, :identity, [xop], elementbytes)
    elseif pint == 2
        return add_compute!(ls, var, :abs2, [xop], elementbytes)
    end

    # Implementation from https://github.com/JuliaLang/julia/blob/a965580ba7fd0e8314001521df254e30d686afbf/base/intfuncs.jl#L216
    t = trailing_zeros(pint) + 1
    pint >>= t
    while (t -= 1) > 0
        varname = (iszero(pint) && isone(t)) ? var : gensym(:pbs)
        xop = add_compute!(ls, varname, :abs2, [xop], elementbytes)
    end
    yop = xop
    while pint > 0
        t = trailing_zeros(pint) + 1
        pint >>= t
        while (t -= 1) >= 0
            xop = add_compute!(ls, gensym(:pbs), :abs2, [xop], elementbytes)
        end
        yop = add_compute!(ls, iszero(pint) ? var : gensym(:pbs), :(*), [xop, yop], elementbytes)
    end
    yop
end


text:
for (elty, ext) in ((:Float32, :s),
                    (:Float64, :d))

    for (matA, matb, sym) in ((:Matrix, :Matrix, "_"),
                              (:DistMatrix, :DistMatrix, "Dist_"),
                              (:SparseMatrix, :Matrix, "Sparse_"),
                              (:DistSparseMatrix, :DistMultiVec, "DistSparse_"))
        @eval begin

            function lav!(A::$matA{$elty}, b::$matb{$elty}, x::$matb{$elty})
                err = ccall(($(string("ElLAV", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, Ptr{Void}, Ptr{Void}),
                    A.obj, b.obj, x.obj)
                err == 0 || throw(ElError(err))
                return x
            end

            function lav!(A::$matA{$elty}, b::$matb{$elty}, x::$matb{$elty}, ctrl::LPAffineCtrl{$elty})
                err = ccall(($(string("ElLAVX", sym, ext)), libEl), Cuint,
                    (Ptr{Void}, Ptr{Void}, Ptr{Void}, LPAffineCtrl{$elty}),
                    A.obj, b.obj, x.obj, ctrl)
                err == 0 || throw(ElError(err))
                return x
            end
        end
    end
end

function lav{T<:Union{Float32,Float64}}(A::Matrix{T}, b::Matrix{T})
    x = Matrix(T)
    return lav!(A, b, x)
end
function lav{T<:Union{Float32,Float64}}(A::DistMatrix{T}, b::DistMatrix{T})
    x = DistMatrix(T, MC, MR, Grid(A))
    return lav!(A, b, x)
end
function lav{T<:Union{Float32,Float64}}(A::DistSparseMatrix{T}, b::DistMultiVec{T})
    x = DistMultiVec(T, comm(A))
    return lav!(A, b, x)
end

function lav{T<:Union{Float32,Float64}}(A::DistSparseMatrix{T}, b::DistMultiVec{T}, ctrl::LPAffineCtrl{T})
    x = DistMultiVec(T, comm(A))
    return lav!(A, b, x, ctrl)
end

ext: jl


text:
include("common.jl")

# bkamins: this function mostly reuses what is below; to be discussed
function calibrate(rs::RideShare, initial_solution::Vector{Vector{Int}})
    temperature = 1000.0
    iterations_per_temp = 1000
    cooling = 0.6
    accept_prob = 1.0
    current_solution = deepcopy(initial_solution)
    i, accepted, total, k = 0, 0.0, 0, 1

    current_score = rideshare_optimality(rs, initial_solution)

    # bkamins: I have change changed the logic of calculation of accept_prob
    # to make it more consistent
    while accept_prob > 0.6
        if i > iterations_per_temp
            accept_prob = total == 0 ? 1.0 : accepted / total
            temperature *= cooling
            accepted, total, i = 0.0, 0, 0
        end
        s_new = iterate_solution(current_solution)
        score_new = rideshare_optimality(rs, s_new)
        if score_new > current_score
            current_score = score_new
            current_solution = s_new
        else
            delta_score = score_new - current_score
            prob = 1 - exp(delta_score / (k * temperature))
            accepted += prob
            total = total + 1
            if prob < rand()
                current_score = score_new
                current_solution = s_new
            end
        end
        i = i + 1
    end
    println("Calibrated Temperature $temperature")
    temperature
end

function simulated_annealing(rs::RideShare;
                             final_temp=0.001, max_iterations=1000,
                             max_accept=2, alpha=0.90, threshold=1000)
    progress = Float64[]
    temp_progress = Float64[]
    current_progress = Float64[]
    solution = get_initial_solution(rs)
    current_temp = calibrate(rs, solution)
    current = rideshare_optimality(rs, solution)
    best = current
    best_schedule = solution
    i, accept, k = 0, 0.0, 1

    println("Initial Solution:")
    foreach(((x, s),) -> println("Driver #$x : $s"), enumerate(solution))
    println("Initial revenue: $current")
    waitbar = Progress(max_iterations, 1)
    while current_temp > final_temp && i < max_iterations && current < threshold
        if accept >= max_accept
            accept = 0
            current_temp *= alpha
        end
        next_schedule = iterate_solution(solution)
        next = rideshare_optimality(rs,next_schedule)
        if next > current
            accept += 1
            solution, current = next_schedule, next
        else
            delta = next - current
            if exp(delta / (k * current_temp)) > rand()
                accept += 1
                solution, current = next_schedule, next
            end
        end
        if current > best
            best, best_schedule = current, solution
        end
        i += 1
        push!(progress, best)
        push!(temp_progress, current_temp)
        push!(current_progress, current)
        next!(waitbar)
    end

    println("Iterations: $i")
    println("Best Solution")
    foreach(((x, s),) -> println("Driver #$x : $s"), enumerate(best_schedule))
    println("Best revenue: $best")
    return (current_progress, best_schedule)
end

using Random
Random.seed!(1234)

res = simulated_annealing(rs_samples[1]; max_iterations=10000, alpha=0.95)
res = simulated_annealing(rs_samples[2]; max_iterations=10000, alpha=0.95)


text:
using Flux
using Random

mutable struct ImageDQNAgent{M, TN, O, LU, AP<:AbstractValuePolicy, Φ, ERP, IMB} <: AbstractAgent
    model::M
    target_network::TN
    opt::O
    lu::LU
    ap::AP
    er::ExperienceReplay{ERP}
    image_buffer::IMB
    γ::Float32
    batch_size::Int64
    tn_counter_init::Int64
    target_network_counter::Int64
    action::Int64
    prev_s::Φ
end

ImageDQNAgent(model, target_network, opt, lu, ap, size_buffer, im, γ, batch_size, tn_counter_init, s) =
    ImageDQNAgent(model,
             target_network,
             opt,
             lu,
             ap,
             ExperienceReplay(size_buffer,
                              (Int64, Int64, Int64, Float32, Bool),
                              (:s, :a, :sp, :r, :t)),
             ImageBuffer(size_buffer+2, im, Array{Float32, 3}),
             γ,
             batch_size,
             tn_counter_init,
             tn_counter_init,
             0,
             0)

function JuliaRL.start!(agent::ImageDQNAgent, env_s_tp1, rng::AbstractRNG; kwargs...)
    # Start an Episode
    agent.action = sample(agent.ap,
                          agent.model(env_s_tp1),
                          rng)
    agent.prev_s = add!(agent.imb, env_s_tp1)
    return agent.action
end

function JuliaRL.step!(agent::ImageDQNAgent, env_s_tp1, r, terminal, rng::AbstractRNG; kwargs...)

    img_idx = add!(agent.imb, env_s_tp1)

    add!(agent.er, (agent.prev_s, agent.action, img_idx, r, terminal))

    if size(agent.er)[1] > 1000
        e = sample(agent.er, agent.batch_size; rng=rng)
        update_params!(agent, e)
    end

    agent.prev_s .= env_s_tp1
    agent.action = sample(agent.ap,
                          agent.model(agent.prev_s),
                          rng)

    return agent.action
end

function update_params!(agent::ImageDQNAgent, e)

    if agent.tn_counter_init > 0

        update!(agent.model, agent.lu, agent.opt, agent.image_buffer[e.s], e.a, agent.image_buffer[e.sp], e.r, e.t, agent.target_network)

        if agent.target_network_counter == 1
            agent.target_network_counter = agent.tn_counter_init
            agent.target_network = deepcopy(mapleaves(
                Flux.Tracker.data,
                agent.model))
        else
            agent.target_network_counter -= 1
        end
    else
        update!(agent.model, agent.lu, agent.opt, agent.image_buffer[e.s], e.a, agent.image_buffer[e.sp], e.r, e.t)
    end
    return nothing

end


text:
# This file is a part of RadiationDetectorSignals.jl, licensed under the MIT License (MIT).

import Test
Test.@testset "Package RadiationDetectorSignals" begin
    include("test_detector_hits.jl")
    include("test_detector_waveforms.jl")
    include("test_docs.jl")
end # testset


text:
function process(env)
    for (msg, key, res) in zip([Dict("cmd"=>"obs_dimensions"), Dict("cmd"=>"n_actions")],
                               ["obs_dim", "n_actions"],
                               [obs_dimensions(env), length(actions(env))])
        respmsg = RLInterface.process!(env, msg)
        @test respmsg[key] == res
    end
end


text:
using Test
using Mimi

@testset "Standard API" begin

# Test that the function does not error and returns a valid value
    scc1 = compute_scc(year=2020)
    @test scc1 isa Float64

# Test that a higher discount rate makes a lower scc value
    scc2 = compute_scc(year=2020, eta=0., prtp=0.03)
    @test scc2 < scc1

# Test with a modified model
    m = get_model()
    set_param!(m, :tcr_transientresponse, 3)
    scc3 = compute_scc(m, year=2020)
    @test scc3 > scc1

# Test get_marginal_model
    mm = get_marginal_model(year=2040)
    mm[:ClimateTemperature, :rt_g_globaltemperature]

# Test compute_scc_mm
    result = compute_scc_mm(year=2050)
    @test result.scc > scc1
    @test result.mm isa Mimi.MarginalModel

# Test Monte Carlo SCC support
    sccs1 = compute_scc(year=2020, n=10, seed=350)
    sccs2 = compute_scc(year=2020, n=10, seed=350)
    @test sccs1 == sccs2

    sccs3 = compute_scc(year=2020, n=10, seed=351)
    @test sccs3 != sccs1

    sccs4 = compute_scc(year=2020, prtp=0.02, eta=1.5, n=10, seed=350)
    @test sccs4 != sccs1    # test that the user specified discounting scheme overrides the default random variable values

# Test equity weighting options

# Test that regional equity weighting inreases the SCC
    @test compute_scc(year=2020, eta=1.5, prtp=0.01, equity_weighting=true) >
        compute_scc(year=2020, eta=1.5, prtp=0.01, equity_weighting=false)

# Test that when eta==0, equity_weighting does not change the SCC value
    @test compute_scc(year=2020, eta=0., prtp=0.03, equity_weighting=true) ≈
        compute_scc(year=2020, eta=0., prtp=0.03, equity_weighting=false) atol = 1e3

# Test Monte Carlo w/ and w/o equity weighting, with the same seed
    scc10a = compute_scc(year=2020, eta=1.5, prtp=0.01, equity_weighting=true, n=10, seed=350)
    scc10b = compute_scc(year=2020, eta=1.5, prtp=0.01, equity_weighting=false, n=10, seed=350)
    @test all(scc10a .> scc10b)

end


text:
@doc raw"""
    PSFModels.Gaussian([T=Float64]; fwhm, x=0, y=0, amp=1, maxsize=3, extent=maxsize .* fwhm)
    PSFModels.Gaussian([T=Float64]; fwhm, pos, amp=1, maxsize=3, extent=maxsize .* fwhm)
    PSFModels.Gaussian([T=Float64]; fwhm, r, theta, amp=1, maxsize=3, extent=maxsize .* fwhm)

An unnormalized bivariate Gaussian distribution. The position can be specified in `(x, y)` coordinates as a `Tuple`, `AbstractVector`, or as separate arguments. By default the model is placed at the origin. The position can also be given as a polar coordinate using `r`/`ρ` and `theta`/`θ`, optionally centered around `origin`.

The `fwhm` can be a scalar (isotropic) or a vector/tuple (diagonal). For efficient calculations, we recommend using [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl). Here, `maxsize` is a multiple of the fwhm, and can be given as a scalar or as a tuple for each axis. The `extent` defines the bounding box for the model and is used for the default rendering size.

# Functional form
```
f(x | x̂, FWHM) = exp[-4ln(2) * ||x - x̂|| / FWHM^2]
```
where `x̂` and `x` are position vectors (indices) `||⋅||` represents the square-distance, and `FWHM` is the full width at half-maximum. If `FWHM` is a scalar, the Gaussian distribution will be isotropic. If `FWHM` is a vector or tuple, the weighting is applied along each axis (diagonal).
"""
struct Gaussian{T,FT,VT<:AbstractVector,IT<:Tuple} <: PSFModel{T}
    fwhm::FT
    pos::VT
    amp::T
    indices::IT

    function Gaussian(fwhm::FT, pos::VT, amp::T, indices::IT) where {T,FT,VT<:AbstractVector,IT<:Tuple}
        new{T,FT,VT,IT}(fwhm, pos, amp, indices)
    end
end

# Alias Normal -> Gaussian
"""
    PSFModels.Normal

An alias for [`PSFModels.Gaussian`](@ref)
"""
const Normal = Gaussian

## constructors
# default type is Float64
Gaussian(; kwargs...) = Gaussian(Float64; kwargs...)
function Gaussian(T; fwhm, amp=one(T), maxsize=3, extent = maxsize .* fwhm, position...)
    # get the position from keyword distpatch
    pos = _position(; position...)
    return Gaussian(fwhm, pos, convert(T, amp), indices_from_extent(pos, extent))
end

Base.size(g::Gaussian) = map(length, g.indices)
Base.axes(g::Gaussian) = g.indices

# Gaussian pre-factor for normalizing the exponential
const GAUSS_PRE = -4 * log(2)

# covers isotropic case
function (g::Gaussian{T})(point::AbstractVector) where T
    Δ = sqeuclidean(point, g.pos)
    val = g.amp * exp(GAUSS_PRE * Δ / g.fwhm^2)
    return convert(T, val)
end
# covers vector case
function (g::Gaussian{T,<:Union{Tuple,AbstractVector}})(point::AbstractVector) where T
    weights = SA[1/first(g.fwhm)^2, 1/last(g.fwhm)^2] # manually invert
    Δ = wsqeuclidean(point, g.pos, weights)
    val = g.amp * exp(GAUSS_PRE * Δ)
    return convert(T, val)
end


text:
module FactorGraph

using SparseArrays, LinearAlgebra
using HDF5, XLSX
using Random
using PrettyTables, Printf

### Form a continuous factor graph and initialize messages and marginals
include("continuousgraph.jl")
export continuousModel, continuousTreeModel

### Form a discrete factor graph and initialize messages and marginals
include("discretegraph.jl")
export discreteTreeModel

### Factor graph manipulation
include("graphmanipulation.jl")
export freezeFactor!, defreezeFactor!, freezeVariable!, defreezeVariable!,
       freezeVariableFactor!, defreezeVariableFactor!, freezeFactorVariable!, defreezeFactorVariable!,
       hideFactor!, addFactors!

### Inference
include("inference.jl")
export marginal, marginalUnnormalized, dynamicFactor!, ageingVariance!, damping!

### Vanilla GBP algorithm
include("vanillaGBP.jl")
export messageFactorVariable, meanFactorVariable, varianceFactorVariable,
       messageDampFactorVariable, meanDampFactorVariable,
       messageVariableFactor, meanVariableFactor, varianceVariableFactor

### Broadcast GBP algorithm
include("broadcastGBP.jl")
export messageFactorVariableBroadcast, meanFactorVariableBroadcast, varianceFactorVariableBroadcast,
       messageDampFactorVariableBroadcast, meanDampFactorVariableBroadcast,
       messageVariableFactorBroadcast, meanVariableFactorBroadcast, varianceVariableFactorBroadcast

### Broadcast GBP with Kahan-Babuska algorithm
include("kahanGBP.jl")
export messageFactorVariableKahan, meanFactorVariableKahan, varianceFactorVariableKahan,
       messageDampFactorVariableKahan, meanDampFactorVariableKahan,
       messageVariableFactorKahan, meanVariableFactorKahan, varianceVariableFactorKahan

# Tree factor graph
include("treeGBP.jl")
export forwardVariableFactor, forwardFactorVariable, backwardVariableFactor, backwardFactorVariable, isTree

# Compute and show results
include("utility.jl")
export wls, errorMetric, displayData

end # FactorGraph


text:
using Test, SafeTestsets

@safetestset "Ridge" begin include("ridge.jl") end
@safetestset "Lasso" begin include("lasso.jl") end


text:
type XPTInfo
    sasversion::String
    operatingsystem::String
    dtcreated::DateTime
    dtmodified::DateTime
end
function Base.show(io::IO, x::XPTInfo)
    println(io, "     SAS Version: ", x.sasversion)
    println(io, "Operating system: ", x.operatingsystem)
    println(io, "         Created: ", x.dtcreated)
    print(io, "        Modified: ", x.dtmodified)
end

function XPTInfo(io::IO)
    seekstart(io)
    first_header_record  = readbytes(io, ROW_LENGTH)
    @assert first_header_record[1:48] == LIBRARY_HEADER_RECORD

    first_real_header = ascii(readbytes(io, ROW_LENGTH))
    second_real_header = ascii(readbytes(io, ROW_LENGTH))

    sasversion = strip(first_real_header[25:32])
    opsys = first_real_header[33:40]
    dtcreated = parsedatetime(first_real_header[end-15:end])
    dtmodified = parsedatetime(second_real_header[1:16])

    XPTInfo(sasversion, opsys, dtcreated, dtmodified)
end


text:
__precompile__()

module CUDAnative

using LLVM
using CUDAdrv
import CUDAdrv: debug, DEBUG, trace, TRACE

const ext = joinpath(@__DIR__, "..", "deps", "ext.jl")
isfile(ext) || error("Unable to load $ext\n\nPlease run Pkg.build(\"CUDAnative\") and restart Julia.")
include(ext)

include("jit.jl")
include("profile.jl")
include(joinpath("device", "array.jl"))
include(joinpath("device", "intrinsics.jl")) # these files contain generated functions,
include("execution.jl")                      # so should get loaded late (JuliaLang/julia#19942)
include("reflection.jl")

function __init__()
    if CUDAdrv.version() != cuda_version ||
        LLVM.version() != llvm_version ||
        VersionNumber(Base.libllvm_version) != julia_llvm_version
        error("Your set-up has changed. Please re-run Pkg.build(\"CUDAnative\") and restart Julia.")
    end
    init_jit()
end

end


text:
"""
    KAEx(device; groupsize, basesize)

Transducers.jl executor implemented using KernelAbstractions.jl.
"""
struct KAEx{D,K} <: Executor
    device::D
    kwargs::K
end

KAEx(device; kwargs...) = KAEx(device, (; kwargs...))

popsimd(; simd = nothing, kwargs...) = kwargs

Transducers.transduce(xf, rf::RF, init, xs, exc::KAEx) where {RF} =
    transduce_ka(xf, rf, init, xs, exc.device; popsimd(; exc.kwargs...)...)

# TODO: Once `groupsize` and `basesize` can be auto-tuned for non-CUDA GPUs,
# hook `KAEx` into the executor promotion mechanism.

function Base.show(io::IO, exc::KAEx)
    @nospecialize exc
    T = typeof(exc.kwargs)
    if (
        fieldtype(typeof(exc), 2) !== typeof(exc.kwargs) ||
        any(i -> typeof(exc.kwargs[i]) != fieldtype(T, i), 1:nfields(exc.kwargs))
    )
        return invoke(show, Tuple{IO,Any}, io, exc)
    end
    print(io, KAEx, "(::", typeof(exc.device), "; ")
    isfirst = true
    for (k, v) in pairs(exc.kwargs)
        isfirst || print(io, ", ")
        isfirst = false
        print(io, k, " = ", v)
    end
    print(io, ")")
    return
end


text:
const LOWERA = codepoint('a')
const LOWERZ = codepoint('z')
const UPPERA = codepoint('A')
const UPPERZ = codepoint('Z')

function rotate(n::Int, s::String)
    String(wraprot.(n,codeunits(s)))
end

# Prettier, but 25% slower for big strings:
# rotate(n::Int, s::String) = map(c->rotate(n,c),s)

rotate(n,c::Char) = Char(wraprot(n,codepoint(c)))

function wraprot(n::Int,c::T)::UInt8 where T <: Unsigned
    if LOWERA <= c <= LOWERZ
        rem(c - T(LOWERA) + T(n), T(26)) + T(LOWERA)
    elseif UPPERA <= c <= UPPERZ
        rem(c - T(UPPERA) + T(n), T(26)) + T(UPPERA)
    else
        c
    end
end

# Shamelessly cribbed from base/reducedim.jl
for n in 1:26
    name = Symbol("R$(n)_str")
    @eval begin
         macro $name(s)
            rotate($n,s)
        end
    end
end

# For testing:
# using BenchmarkTools
# @benchmark rotate(23,s) setup=(s=rand('a':'z',100_000)|>join)


text:
# Betweenness centrality measures
# TODO - weighted, separate unweighted, edge betweenness

@doc_str """
    betweenness_centrality(g[, vs])
    betweenness_centrality(g, k)
    parallel_betweenness_centrality(g[, vs])
    parallel_betweenness_centrality(g, k)
Calculate the [betweenness centrality](https://en.wikipedia.org/wiki/Centrality#Betweenness_centrality)
of a graph `g` across all vertices, a specified subset of vertices `vs`, or a random subset of `k`
vertices. Return a vector representing the centrality calculated for each node in `g`.
### Optional Arguments
- `normalize=true`: If true, normalize the betweenness values by the
total number of possible distinct paths between all pairsin the graphs.
For an undirected graph, this number is ``\\frac{(|V|-1)(|V|-2)}{2}``
and for a directed graph, ``{(|V|-1)(|V|-2)}``.
- `endpoints=false`: If true, include endpoints in the shortest path count.
Betweenness centrality is defined as:
``
bc(v) = \\frac{1}{\\mathcal{N}} \sum_{s \\neq t \\neq v}
\\frac{\\sigma_{st}(v)}{\\sigma_{st}}
``.
### References
- Brandes 2001 & Brandes 2008
"""
function betweenness_centrality(
    g::AbstractGraph,
    vs::AbstractVector = vertices(g);
    normalize=true,
    endpoints=false)

    n_v = nv(g)
    k = length(vs)
    isdir = is_directed(g)

    betweenness = zeros(n_v)
    for s in vs
        if degree(g, s) > 0  # this might be 1?
            state = dijkstra_shortest_paths(g, s; allpaths=true, trackvertices=true)
            if endpoints
                _accumulate_endpoints!(betweenness, state, g, s)
            else
                _accumulate_basic!(betweenness, state, g, s)
            end
        end
    end

    _rescale!(betweenness,
    n_v,
    normalize,
    isdir,
    k)

    return betweenness
end

betweenness_centrality(g::AbstractGraph, k::Integer; normalize=true, endpoints=false) =
    betweenness_centrality(g, sample(vertices(g), k); normalize=normalize, endpoints=endpoints)

function parallel_betweenness_centrality(
    g::AbstractGraph,
    vs::AbstractVector = vertices(g);
    normalize=true,
    endpoints=false)::Vector{Float64}

    n_v = nv(g)
    k = length(vs)
    isdir = is_directed(g)

    # Parallel reduction

    betweenness = @parallel (+) for s in vs
        temp_betweenness = zeros(n_v)
        if degree(g, s) > 0  # this might be 1?
            state = dijkstra_shortest_paths(g, s; allpaths=true, trackvertices=true)
            if endpoints
                _accumulate_endpoints!(temp_betweenness, state, g, s)
            else
                _accumulate_basic!(temp_betweenness, state, g, s)
            end
        end
        temp_betweenness
    end

    _rescale!(betweenness,
    n_v,
    normalize,
    isdir,
    k)

    return betweenness
end

parallel_betweenness_centrality(g::AbstractGraph, k::Integer; normalize=true, endpoints=false) =
    parallel_betweenness_centrality(g, sample(vertices(g), k); normalize=normalize, endpoints=endpoints)

function _accumulate_basic!(
    betweenness::Vector{Float64},
    state::DijkstraState,
    g::AbstractGraph,
    si::Integer
    )

    n_v = length(state.parents) # this is the ttl number of vertices
    δ = zeros(n_v)
    σ = state.pathcounts
    P = state.predecessors

    # make sure the source index has no parents.
    P[si] = []
    # we need to order the source vertices by decreasing distance for this to work.
    S = reverse(state.closest_vertices) #Replaced sortperm with this
    for w in S
        coeff = (1.0 + δ[w]) / σ[w]
        for v in P[w]
            if v > 0
                δ[v] += (σ[v] * coeff)
            end
        end
        if w != si
            betweenness[w] += δ[w]
        end
    end
    return nothing
end

function _accumulate_endpoints!(
    betweenness::Vector{Float64},
    state::DijkstraState,
    g::AbstractGraph,
    si::Integer
    )

    n_v = nv(g) # this is the ttl number of vertices
    δ = zeros(n_v)
    σ = state.pathcounts
    P = state.predecessors
    v1 = [1:n_v;]
    v2 = state.dists
    S = reverse(state.closest_vertices)
    s = vertices(g)[si]
    betweenness[s] += length(S) - 1    # 289

    for w in S
        coeff = (1.0 + δ[w]) / σ[w]
        for v in P[w]
            δ[v] += σ[v] * coeff
        end
        if w != si
            betweenness[w] += (δ[w] + 1)
        end
    end
    return nothing
end

function _rescale!(betweenness::Vector{Float64}, n::Integer, normalize::Bool, directed::Bool, k::Int)
    if normalize
        if n <= 2
            do_scale = false
        else
            do_scale = true
            scale = 1.0 / ((n - 1) * (n - 2))
        end
    else
        if !directed
            do_scale = true
            scale = 1.0 / 2.0
        else
            do_scale = false
        end
    end
    if do_scale
        if k > 0
            scale = scale * n / k
        end
        for v = 1:length(betweenness)
            betweenness[v] *= scale
        end
    end
    return nothing
end


text:
include("unit_test_helper.jl")

# Turing

using Turing

@model gdemo() = begin
    s ~ InverseGamma(2, 3)
    m ~ Normal(0,sqrt(s))
    1.5 ~ Normal(m, sqrt(s))
    2.0 ~ Normal(m, sqrt(s))
    return s, m
  end

# Plain Julia

using ReverseDiff: GradientTape, GradientConfig, gradient, gradient!, compile
using Turing: invlink, logpdf

θ_dim = 2
function lj_func(θ)
  _lj = zero(Real)

  d_s = InverseGamma(2, 3)
  s = invlink(d_s, θ[1])
  _lj += logpdf(d_s, s, true)
  m = θ[2]
  _lj += logpdf(Normal(0, sqrt(s)), m)

  _lj += logpdf(Normal(m, sqrt(s)), 1.5)
  _lj += logpdf(Normal(m, sqrt(s)), 2.0)

  return _lj
end

neg_lj_func(θ) = -lj_func(θ)
const f_tape = GradientTape(neg_lj_func, randn(θ_dim))
const compiled_f_tape = compile(f_tape)

function grad_func(θ)

  inputs = θ
  results = similar(θ)
  all_results = DiffResults.GradientResult(results)

  gradient!(all_results, compiled_f_tape, inputs)

  neg_lj = all_results.value
  grad, = all_results.derivs

  return -neg_lj, grad

end

# Unit test for gradient

test_grad(gdemo, grad_func; trans=[1])


text:
#
# Parrot.jl/src/layers.jl
#
# By Elijah Mathews (me@elijahmathews.com)
# openpgp4fpr:ac1d3fb1e8a5eb7d14bd587b2932c725055a90d8
#

"""
    Alsing(in::Integer, out::Integer)

Create a non-linear `Alsing` layer with trainable parameters `W`, `b`,
`α`, and `β`.

    y = (β .+ σ.(α .* (W*x .+ b)) .* (1 .- β)) .* (W*x .+ b)

The input `x` must be a vector of length `in`, or a batch of vectors
represented as an `in × N` matrix. The out `y` will be a vector or batch
of length `out`.

See [Alsing et al. (2020)](https://doi.org/10.3847/1538-4365/ab917f)
for more information about this activation function. Intended to be used for
emulating stellar population synthesis codes.

# Example
```
julia> a = Alsing(8, 3)
Alsing(8, 3)

julia> a(rand(8))
3-element Array{Float32,1}:
 0.9413336
 0.30788675
 0.5125884
```
"""
struct Alsing{S<:AbstractArray, T<:AbstractArray}
    W::S
    b::T
    α::T
    β::T
end

function Alsing(in::Integer, out::Integer;
                initW = Flux.glorot_uniform, initb = Flux.zeros,
                initα = Flux.glorot_uniform, initβ = Flux.glorot_uniform)
    return Alsing(initW(out,in), initb(out), initα(out), initβ(out))
end

Flux.@functor Alsing

function (a::Alsing)(x::AbstractArray)
    W, b, α, β = a.W, a.b, a.α, a.β
    (β .+ Flux.σ.(α .* (W*x .+ b)) .* (1 .- β)) .* (W*x .+ b)
end

function Base.show(io::IO, l::Alsing)
    print(io, "Alsing(", size(l.W, 2), ", ", size(l.W, 1), ")")
end

#
# Efficiency hack.
#
(a::Alsing{W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
    invoke(a, Tuple{AbstractArray}, x)

(a::Alsing{W})(x::AbstractArray{<:AbstractFloat}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
    a(T.(x))

"""
    TransformPCA(P::MultivariateStats.PCA{AbstractFloat})
    TransformPCA(P::OnlineStats.CCIPCA)

Create a `TransformPCA` layer that uses Principal Component Analysis (PCA)
to transform its input into the PCA basis using a given PCA object `P`.

The input `x` must have dimensions compatible with `P` and the layer's output
will be of the dimensions given by the PCA transformation. No parameters are
trainable.

See also: [`ReconstructPCA`](@ref)
"""
struct TransformPCA{S<:AbstractArray, T<:AbstractArray}
    Pt::S
    μ::T
end

function TransformPCA(P::MultivariateStats.PCA{S}) where {S<:AbstractFloat}
    return TransformPCA(permutedims(P.proj), P.mean)
end

function TransformPCA(P::OnlineStats.CCIPCA)
    return TransformPCA(permutedims(P.U), P.center)
end

Flux.@functor TransformPCA

Flux.trainable(t::TransformPCA) = ()

function (t::TransformPCA)(x::AbstractArray)
    Pt, μ = t.Pt, t.μ
    Pt * (x .- μ)
end

function Base.show(io::IO, l::TransformPCA)
    print(io, "TransformPCA(", size(l.Pt, 2), ", ", size(l.Pt, 1), ")")
end

#
# Efficiency hack.
#
(t::TransformPCA{W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
    invoke(t, Tuple{AbstractArray}, x)

(t::TransformPCA{W})(x::AbstractArray{<:AbstractFloat}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
    a(T.(x))

"""
    ReconstructPCA(P::MultivariateStats.PCA{AbstractFloat})
    ReconstructPCA(P::OnlineStats.CCIPCA)

Create a `ReconstructPCA` layer that uses Principal Component Analysis (PCA)
to reconstruct its input from the PCA basis using a given PCA object `P`.

The input `x` must have dimensions compatible with `P` and the layer's output
will be of the dimensions given by the PCA reconstruction. No parameters are
trainable.

See also: [`TransformPCA`](@ref)
"""
struct ReconstructPCA{S<:AbstractArray, T<:AbstractArray}
    P::S
    μ::T
end

function ReconstructPCA(P::MultivariateStats.PCA{S}) where {S<:AbstractFloat}
    return ReconstructPCA(P.proj, P.mean)
end

function ReconstructPCA(P::OnlineStats.CCIPCA)
    return ReconstructPCA(P.U, P.center)
end

Flux.@functor ReconstructPCA

Flux.trainable(r::ReconstructPCA) = ()

function (r::ReconstructPCA)(x::AbstractArray)
    P, μ = r.P, r.μ
    P * x .+ μ
end

function Base.show(io::IO, l::ReconstructPCA)
    print(io, "ReconstructPCA(", size(l.P, 2), ", ", size(l.P, 1), ")")
end

#
# Efficiency hack.
#
(r::ReconstructPCA{W})(x::AbstractArray{T}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
    invoke(r, Tuple{AbstractArray}, x)

(r::ReconstructPCA{W})(x::AbstractArray{<:AbstractFloat}) where {T <: Union{Float32,Float64}, W <: AbstractArray{T}} =
    r(T.(x))

"""
    Normalize(μ::AbstractArray, σ::AbstractArray)

Create a simple `Normalize` layer with parameters consisting of
some previously known mean `μ` and standard deviation `σ` that are used
to normalize its input.

    y = (x .- μ) ./ σ

The input `x` must be a vector of equal length to both `μ` and `σ` or
an array with dimensions such that the broadcasted functions can be used
with the given `μ` and `σ`. No parameters are trainable.

See also: [`Denormalize`](@ref)
"""
struct Normalize{S<:AbstractArray}
    μ::S
    σ::S
end

function Normalize(μ::AbstractArray, σ::AbstractArray)
    return Normalize(μ, σ)
end

Flux.@functor Normalize

Flux.trainable(n::Normalize) = ()

function (n::Normalize)(x::AbstractArray)
    μ, σ = n.μ, n.σ
    (x .- μ) ./ σ
end

function Base.show(io::IO, l::Normalize)
    print(io, "Normalize(", size(l.μ, 1), ", ", size(l.μ, 1), ")")
end

"""
    Denormalize(μ::AbstractArray, σ::AbstractArray)

Create a simple `Denormalize` layer with parameters consisting of
some previously known mean `μ` and standard deviation `σ` that are used
to restore its input from a normalized state.

    y = (σ .* x) .+ μ

The input `x` must be a vector of equal length to both `μ` and `σ` or
an array with dimensions such that the broadcasted functions can be used
with the given `μ` and `σ`. No parameters are trainable.

See also: [`Normalize`](@ref)
"""
struct Denormalize{S<:AbstractArray}
    μ::S
    σ::S
end

function Denormalize(μ::AbstractArray, σ::AbstractArray)
    return Denormalize(μ, σ)
end

Flux.@functor Denormalize

Flux.trainable(d::Denormalize) = ()

function (d::Denormalize)(x::AbstractArray)
    μ, σ = d.μ, d.σ
    (σ .* x) .+ μ
end

function Base.show(io::IO, l::Denormalize)
    print(io, "Denormalize(", size(l.μ, 1), ", ", size(l.μ, 1), ")")
end


text:
#*************************************************************************#
#
# Series Model Type
#
#*************************************************************************#
struct SeriesModel{M<:StaticModel} <: ParametricModel
    ninputs::Int
    noutputs::Int
    nparams::Int
    nmodels::Int
    mdls::Vector{M}

    # Internal buffers
    vec_z::Vector{AbstractVector{Float64}}
    vec_dx::Vector{AbstractMatrix{Float64}}
    vec_dΘ::Vector{AbstractMatrix{Float64}}
end

"""
    SeriesModel([mdl1, mdl2, mdl3,..., mdln])

Initialize `ParametricModel` representing composition of models
`z = mdln(... mdl3(mdl2(mdl1(x))))`.
"""
function SeriesModel{M<:StaticModel}(mdls::Vector{M})

    # Check coherence between models
    nmodels = length(mdls)
    for i=1:length(mdls)-1
        if mdls[i].noutputs != mdls[i+1].ninputs
            throw(DimensionMismatch(string("The number of outputs of a ",
                                           "model should match the number ",
                                           "of inputs of the next one.")))
        end
    end

    # Get input, outputs and number of parameters
    ninputs = mdls[1].ninputs
    noutputs = mdls[end].noutputs

    nparams = 0
    nparametricmodels = 0
    for mdl in mdls
        if isparametric(mdl)
            nparams += mdl.nparams
        end
    end

    # Check if there are nested Series Models
    expand_model = false
    for mdl in mdls
        if typeof(mdl) <: SeriesModel
            nmodels += mdl.nmodels-1
            expand_model = true
        end
    end

    # Expand Series Models if needed
    if expand_model
        new_mdls = Vector{StaticModel}(nmodels)

        k = 1
        for mdl in mdls
            if typeof(mdl) <: SeriesModel
                for internal_mdl in mdl.mdls
                    new_mdls[k] = internal_mdl
                    k += 1
                end
            else
                new_mdls[k] = mdl
                k += 1
            end
        end
        mdls = new_mdls
    end

    # ---
    # Allocate Memory
    # ---
    vec_z = Vector{AbstractVector{Float64}}(nmodels+1)
    vec_dx = Vector{AbstractMatrix{Float64}}(nmodels)
    vec_dΘ = Vector{AbstractMatrix{Float64}}(nmodels)

    vec_z[1] = Vector{Float64}(ninputs)
    for (i, mdl) in Enumerate(mdls)
        vec_z[i+1] = z_buffer(mdl)
        vec_dx[i] = dx_buffer(mdl)
        if isparametric(mdl)
            vec_dΘ[i] = dΘ_buffer(mdl)
        end
    end
    # ---

    SeriesModel(ninputs, noutputs, nparams, nmodels, mdls, vec_z, vec_dx, vec_dΘ)
end

#*************************************************************************#
#
# Compute Output and its Derivatives
#
#*************************************************************************#
function _evaluate!(seriesmdl::SeriesModel,
                    x::AbstractVector{<:Real},
                    Θ::AbstractVector{<:Real},
                    z::AbstractVector{Float64},
                    dx::AbstractMatrix{Float64},
                    dΘ::AbstractMatrix{Float64})

    # Get memory allocated for vec_z
    vec_z = seriesmdl.vec_z
    vec_z[1] .= x

    # Empty auxiliar matrix
    empty = Vector{AbstractMatrix{Float64}}(seriesmdl.nmodels)
    fill!(empty, (Matrix{Float64}(0, 0)))

    # Use memory allocated for vec_dx and vec_dΘ
    if isempty(dx) && isempty(dΘ)
        vec_dx = empty
        vec_dΘ = empty
    elseif !isempty(dΘ)
        vec_dx = seriesmdl.vec_dx
        vec_dΘ = seriesmdl.vec_dΘ
    else
        vec_dx[i] = seriesmdl.vec_dx
        vec_dΘ[i] = empty
    end

    # ---
    # Forward Stage
    # ---
    k = 1
    for (i, mdl) in Enumerate(seriesmdl.mdls)
        if isparametric(mdl)
            _evaluate!(mdl, vec_z[i], Θ[k:k+mdl.nparams-1],
                       vec_z[i+1], vec_dx[i], vec_dΘ[i])
            k += mdl.nparams
        else
            _evaluate!(mdl, vec_z[i], vec_z[i+1], vec_dx[i])
        end
    end

    z .= vec_z[end]
    # ---

    if isempty(dx) && isempty(dΘ)
        return
    end

    # ---
    # Backward Stage
    # ---
    node_derivative = Matrix{Float64} # TODO: reuse the same `node_derivative` buffer for all `i`
    for i = seriesmdl.nmodels:-1:1

        # Get static model
        mdl = seriesmdl.mdls[i]

        # compute dΘ
        if !isempty(dΘ) && isparametric(mdl)
            k -= mdl.nparams
            if i == seriesmdl.nmodels
                dΘ[1:end, k:k+mdl.nparams-1] .= vec_dΘ[i]
            else
                dΘ[1:end, k:k+mdl.nparams-1] .= node_derivative * vec_dΘ[i]
            end
        end

        # Propagate derivatives backwards
        if i == seriesmdl.nmodels
            node_derivative = vec_dx[i]
        elseif  i>1 || !isempty(dx)
            node_derivative = node_derivative * vec_dx[i]
        end
    end

    if  !isempty(dx)
        dx .= node_derivative
    end
    # ---
end

#*************************************************************************#
#
# Initializer
#
#*************************************************************************#
"""
    ComposedInitializer()

The `ComposedInitializer` will initialize the parameter vector
according to its subinitializers.

Valid only for series static model.
"""
struct ComposedInitializer{M<:Initializer} <: Initializer
    init_vec::Vector{M}
end

function initial_guess(seriesmdl::SeriesModel, init::ComposedInitializer)

    Θ = Vector{Float64}(seriesmdl.nparams)
    k = 1
    i = 1
    for mdl in seriesmdl.mdls
        if isparametric(mdl)
            Θ[k:k+mdl.nparams-1] = initial_guess(mdl, init.init_vec[i])
            k += mdl.nparams
            i += 1
        end
    end
    return Θ
end

function initial_guess(seriesmdl::SeriesModel, init::DefaultInitializer)

    # Use default initializer in each one of the submodels
    init_vec = fill(DefaultInitializer(), seriesmdl.nmodels)
    initial_guess(seriesmdl, ComposedInitializer(init_vec))
end


text:
_NAMES_TO_ACCEPTEDNAMES = Dict{Symbol,Array{Symbol}}(
    :FavettoSamson => [:favettosamson, :favetto_samson, :favettoandsamson, :favetto_and_samson],
    :FavettoSamsonAux => [:favettosamsonaux, :favetto_samson_aux, :favettoandsamsonaux, :favetto_and_samson_aux],
    :FitzHughNagumo => [:fitzhughnagumo, :fitz_hugh_nagumo],
    :FitzHughNagumoAlt => [:fitzhughnagumoalt, :fitz_hugh_nagumo_alt],
    :FitzHughNagumoConjug => [:fitzhughnagumoconjug, :fitz_hugh_nagumo_conjug],
    :FitzHughNagumoAux => [:fitzhughnagumoaux, :fitz_hugh_nagumo_aux],
    :FitzHughNagumoAltAuxSimple => [:fitzhughnagumoaltauxsimple, :fitz_hugh_nagumo_alt_aux_simple],
    :FitzHughNagumoAltAuxLin => [Symbol(s, e) for s in [:fitzhughnagumoaltaux, :fitz_hugh_nagumo_alt_aux_] for e in [:lin, :linear, :linearised]],
    :FitzHughNagumoConjugAuxSimple => [:fitzhughnagumoconjugauxsimple, :fitz_hugh_nagumo_conjug_aux_simple],
    :FitzHughNagumoConjugAuxLin => [Symbol(s, e) for s in [:fitzhughnagumoconjugaux, :fitz_hugh_nagumo_conjug_aux_] for e in [:lin, :linear, :linearised]],
    :JansenRit => [:jansenrit, :jansen_rit, :jansen_and_rit, :jansenandrit],
    :JansenRitAux => [:jansenritaux, :jansen_rit_aux, :jansen_and_rit_aux, :jansenandritaux],
    :Lorenz => [:lorenz],
    :LorenzAux => [:lorenzaux, :lorenz_aux],
    :Lorenz96 => [:lorenz96, :lorenz_96],
    :Lorenz96Aux => [:lorenz96aux, :lorenz_96_aux],
    :LotkaVolterra => [:lotka_volterra, :lotkavolterra],
    :LotkaVolterraAux => [:lotka_volterra_aux, :lotkavolterraaux],
    :Prokaryote => [:prokaryote, :prokaryotic],
    :ProkaryoteAux => [:prokaryoteaux, :prokaryoticaux, :prokaryote_aux, :prokaryotic_aux],
    :Sine => [:sine],
    :SineAux => [:sineaux, :sine_aux],
    :SIR => [:sir],
    :SIRAux => [:sir, :sir_aux],
)

_ADMISSIBLENAMES = collect(Iterators.flatten(values(_NAMES_TO_ACCEPTEDNAMES)))

_ACCEPTEDNAMES_TO_NAMES = Dict{Symbol,Symbol}()
for key in keys(_NAMES_TO_ACCEPTEDNAMES)
    for v in _NAMES_TO_ACCEPTEDNAMES[key]
        _ACCEPTEDNAMES_TO_NAMES[v] = key
    end
end

_NAMES_TO_PATH = Dict{Symbol, String}(
    :FavettoSamson => joinpath("favetto_samson","favetto_samson.jl"),
    :FavettoSamsonAux => joinpath("favetto_samson","favetto_samson_aux.jl"),
    :FitzHughNagumo => joinpath("fitz_hugh_nagumo", "fitz_hugh_nagumo.jl"),
    :FitzHughNagumoAlt => joinpath("fitz_hugh_nagumo", "fitz_hugh_nagumo_alt.jl"),
    :FitzHughNagumoConjug => joinpath("fitz_hugh_nagumo", "fitz_hugh_nagumo_conjug.jl"),
    :FitzHughNagumoAux => joinpath("fitz_hugh_nagumo", "fitz_hugh_nagumo_aux.jl"),
    :FitzHughNagumoAltAuxSimple => joinpath("fitz_hugh_nagumo", "fitz_hugh_nagumo_alt_aux_simple.jl"),
    :FitzHughNagumoAltAuxLin => joinpath("fitz_hugh_nagumo", "fitz_hugh_nagumo_alt_aux_linearised.jl"),
    :FitzHughNagumoConjugAuxSimple => joinpath("fitz_hugh_nagumo", "fitz_hugh_nagumo_conjug_aux_simple.jl"),
    :FitzHughNagumoConjugAuxLin => joinpath("fitz_hugh_nagumo", "fitz_hugh_nagumo_conjug_aux_linearised.jl"),
    :JansenRit => joinpath("jansen_rit","jansen_rit.jl"),
    :JansenRitAux => joinpath("jansen_rit","jansen_rit_aux.jl"),
    :Lorenz => joinpath("lorenz_system", "lorenz_system.jl"),
    :LorenzAux => joinpath("lorenz_system","lorenz_system_aux.jl"),
    :Lorenz96 => joinpath("lorenz_96","lorenz_96.jl"),
    :Lorenz96Aux => joinpath("lorenz_96","lorenz_96_aux.jl"),
    :LotkaVolterra => joinpath("lotka_volterra", "lotka_volterra.jl"),
    :LotkaVolterraAux => joinpath("lotka_volterra", "lotka_volterra_aux.jl"),
    :Prokaryote => joinpath("prokaryote","prokaryote.jl"),
    :ProkaryoteAux => joinpath("prokaryote","prokaryote_aux.jl"),
    :Sine => joinpath("sine", "sine.jl"),
    :SineAux => joinpath("sine", "sine_aux.jl"),
    :SIR => joinpath("sir", "sir.jl"),
    :SIRAux => joinpath("sir", "sir_aux.jl"),
)


text:
module test_millboard_strings

using Millboard
using Test

Millboard.set_table_mode(:markdown)

board = ["Lorem ipsum dolor sit amet"]
@test sprint(show, table(board)) == """
|   | 1                          |
|---|----------------------------|
| 1 | Lorem ipsum dolor sit amet |"""

board = ["Lorem ipsum\ndolor sit amet"]
@test sprint(show, table(board)) == """
|   | 1              |
|---|----------------|
|   | Lorem ipsum    |
| 1 | dolor sit amet |"""

@test sprint(show, table(board, rownames=["first"])) == """
|       | 1              |
|-------|----------------|
|       | Lorem ipsum    |
| first | dolor sit amet |"""

@test sprint(show, table(board, colnames=["first\ncolumn"], rownames=["first row"])) == """
|           | first          |
|           | column         |
|-----------|----------------|
|           | Lorem ipsum    |
| first row | dolor sit amet |"""

board = ["Lorem ipsum\ndolor sit amet" 42]
@test sprint(show, table(board, colnames=["first\ncolumn"], rownames=["first row"])) == """
|           | first          | 2  |
|           | column         |    |
|-----------|----------------|----|
|           | Lorem ipsum    | 42 |
| first row | dolor sit amet |    |"""

end # module test_millboard_strings


text:
@require Revise begin
    Revise.track(Plots, joinpath(Pkg.dir("Plots"), "src", "backends", "plotlyjs.jl"))
end

# https://github.com/spencerlyon2/PlotlyJS.jl

const _plotlyjs_attr        = _plotly_attr
const _plotlyjs_seriestype  = _plotly_seriestype
const _plotlyjs_style       = _plotly_style
const _plotlyjs_marker      = _plotly_marker
const _plotlyjs_scale       = _plotly_scale

# --------------------------------------------------------------------------------------

function add_backend_string(::PlotlyJSBackend)
    """
    if !Plots.is_installed("PlotlyJS")
        Pkg.add("PlotlyJS")
    end
    if !Plots.is_installed("Rsvg")
        Pkg.add("Rsvg")
    end
    import Blink
    Blink.AtomShell.install()
    """
end

function _initialize_backend(::PlotlyJSBackend; kw...)
    @eval begin
        import PlotlyJS
        export PlotlyJS
    end

    # # override IJulia inline display
    # if isijulia()
    #     IJulia.display_dict(plt::AbstractPlot{PlotlyJSBackend}) = IJulia.display_dict(plt.o)
    # end
end

# ---------------------------------------------------------------------------

function _create_backend_figure(plt::Plot{PlotlyJSBackend})
    if !isplotnull() && plt[:overwrite_figure] && isa(current().o, PlotlyJS.SyncPlot)
        PlotlyJS.SyncPlot(PlotlyJS.Plot(), current().o.view)
    else
        PlotlyJS.plot()
    end
end

function _series_added(plt::Plot{PlotlyJSBackend}, series::Series)
    syncplot = plt.o
    pdicts = plotly_series(plt, series)
    for pdict in pdicts
        typ = pop!(pdict, :type)
        gt = PlotlyJS.GenericTrace(typ; pdict...)
        PlotlyJS.addtraces!(syncplot, gt)
    end
end

function _series_updated(plt::Plot{PlotlyJSBackend}, series::Series)
    xsym, ysym = (ispolar(series) ? (:t,:r) : (:x,:y))
    kw = KW(xsym => (series.d[:x],), ysym => (series.d[:y],))
    z = series[:z]
    if z != nothing
        kw[:z] = (isa(z,Surface) ? transpose_z(series, series[:z].surf, false) : z,)
    end
    PlotlyJS.restyle!(
        plt.o,
        findfirst(plt.series_list, series),
        kw
    )
end

# ----------------------------------------------------------------

function _update_plot_object(plt::Plot{PlotlyJSBackend})
    pdict = plotly_layout(plt)
    syncplot = plt.o
    w,h = plt[:size]
    PlotlyJS.relayout!(syncplot, pdict, width = w, height = h)
end

# ----------------------------------------------------------------

function Base.show(io::IO, ::MIME"text/html", plt::Plot{PlotlyJSBackend})
    prepare_output(plt)
    if isijulia() && !_use_remote[]
        write(io, PlotlyJS.html_body(PlotlyJS.JupyterPlot(plt.o)))
    else
        show(io, MIME("text/html"), plt.o)
    end
end

function plotlyjs_save_hack(io::IO, plt::Plot{PlotlyJSBackend}, ext::String)
    tmpfn = tempname() * "." * ext
    PlotlyJS.savefig(plt.o, tmpfn)
    write(io, read(open(tmpfn)))
end
_show(io::IO, ::MIME"image/svg+xml", plt::Plot{PlotlyJSBackend}) = plotlyjs_save_hack(io, plt, "svg")
_show(io::IO, ::MIME"image/png", plt::Plot{PlotlyJSBackend}) = plotlyjs_save_hack(io, plt, "png")
_show(io::IO, ::MIME"application/pdf", plt::Plot{PlotlyJSBackend}) = plotlyjs_save_hack(io, plt, "pdf")
_show(io::IO, ::MIME"image/eps", plt::Plot{PlotlyJSBackend}) = plotlyjs_save_hack(io, plt, "eps")

function write_temp_html(plt::Plot{PlotlyJSBackend})
    filename = string(tempname(), ".html")
    savefig(plt, filename)
    filename
end

function _display(plt::Plot{PlotlyJSBackend})
    if get(ENV, "PLOTS_USE_ATOM_PLOTPANE", true) in (true, 1, "1", "true", "yes")
        display(plt.o)
    else
        standalone_html_window(plt)
    end
end

function closeall(::PlotlyJSBackend)
    if !isplotnull() && isa(current().o, PlotlyJS.SyncPlot)
        close(current().o)
    end
end


text:
struct GEQDSKFile
    file::String                    # Source file
    nw::Int                         # Number of horizontal R grid points
    nh::Int                         # Number of vertical Z grid points
    r::AbstractRange{Float64}       # R grid points
    z::AbstractRange{Float64}       # Z grid points
    rdim::Float64                   # Horizontal dimension in meter of computational box
    zdim::Float64                   # Vertical dimension in meter of computational box
    rleft::Float64                  # Minimum R in meter of rectangular computational box
    zmid::Float64                   # Z of center of computational box in meter
    nbbbs::Int                      # Number of boundary points
    rbbbs::Vector{Float64}          # R of boundary points in meter
    zbbbs::Vector{Float64}          # Z of boundary points in meter
    limitr::Int                     # Number of limiter points
    rlim::Vector{Float64}           # R of surrounding limiter contour in meter
    zlim::Vector{Float64}           # Z of surrounding limiter contour in meter
    rcentr::Float64                 # R in meter of vacuum toroidal magnetic field BCENTR
    bcentr::Float64                 # Vacuum toroidal magnetic field in Tesla at RCENTR
    rmaxis::Float64                 # R of magnetic axis in meter
    zmaxis::Float64                 # Z of magnetic axis in meter
    simag::Float64                  # Poloidal flux at magnetic axis in Weber/rad
    sibry::Float64                  # Poloidal flux at the plasma boundary in Weber/rad
    psi::AbstractRange{Float64}     # Poloidal flux points
    current::Float64                # Plasma current in Ampere
    fpol::Vector{Float64}           # Poloidal current function in m-T, F = RBT on flux grid
    pres::Vector{Float64}           # Plasma pressure in nt / m2 on uniform flux grid
    ffprim::Vector{Float64}         # FF'(psi) in (mT)2 / (Weber/rad) on uniform fl ux grid
    pprime::Vector{Float64}         # P'(psi) in (nt/m2) / (Weber/rad) on uniform flux grid
    qpsi::Vector{Float64}           # q values on uniform flu x grid from axis to boundary
    psirz::Matrix{Float64}          # Poloidal flux in Weber/rad on the rectangular grid points
end

function Base.show(io::IO, g::GEQDSKFile)
    print(io,"GEQDSKFile: \"",g.file,"\"")
end

Base.broadcastable(g::GEQDSKFile) = (g,)

function file_numbers(f)
    Channel(ctype=Float64) do c
        while ~eof(f)
            line = readline(f)
            pattern = r"[+-]?\d*[\.]?\d+(?:[Ee][+-]?\d+)?"
            for p in eachmatch(pattern,line)
                push!(c,parse(Float64,p.match))
            end
        end
    end
end

function read_array(t,n)
    data = zeros(n)
    for i=1:n
        data[i] = take!(t)
    end
    return data
end

function read_array2d(t,n,m)
    data = zeros(n,m)
    for i=1:n, j=1:m
        data[i,j] = take!(t)
    end
    return data'
end

function readg(gfile)

    isfile(gfile) || error("$(gfile) does not exist")

    f = open(gfile)

    desc = readline(f)
    s = split(desc)

    time = Meta.parse(s[end-3])
    idum = Meta.parse(s[end-2])
    nw = Meta.parse(s[end-1])
    nh = Meta.parse(s[end])

    token = file_numbers(f)

    rdim   = take!(token)
    zdim   = take!(token)
    rcentr = take!(token)
    rleft  = take!(token)
    zmid   = take!(token)

    rmaxis = take!(token)
    zmaxis = take!(token)
    simag  = take!(token)
    sibry  = take!(token)
    bcentr = take!(token)

    current = take!(token)
    simag  = take!(token)
    xdum    = take!(token)
    rmaxis  = take!(token)
    xdum    = take!(token)

    zmaxis = take!(token)
    xdum   = take!(token)
    sibry = take!(token)
    xdum   = take!(token)
    xdum   = take!(token)

    fpol = read_array(token,nw)
    pres = read_array(token,nw)
    ffprim = read_array(token,nw)
    pprime = read_array(token,nw)
    psirz = read_array2d(token,nh,nw)
    qpsi = read_array(token,nw)

    nbbbs = Int(take!(token))
    limitr = Int(take!(token))

    if nbbbs > 0
        rbbbs = zeros(nbbbs)
        zbbbs = zeros(nbbbs)
        for i=1:nbbbs
            rbbbs[i] = take!(token)
            zbbbs[i] = take!(token)
        end
    else
        rbbbs = [0.0]
        zbbbs = [0.0]
    end

    if limitr > 0
        rlim = zeros(limitr)
        zlim = zeros(limitr)
        for i=1:limitr
            rlim[i] = take!(token)
            zlim[i] = take!(token)
        end
    else
        rlim = [0.0]
        zlim = [0.0]
    end

    close(f)
    close(token)

    r = range(rleft, rleft + rdim, length=nw)
    z = range(zmid - 0.5*zdim, zmid + 0.5*zdim, length=nh)
    psi = range(simag, sibry, length=nw)

    g = GEQDSKFile(gfile, nw,nh,r,z,rdim,zdim,rleft,zmid,nbbbs,rbbbs,zbbbs,limitr,rlim,zlim,
                   rcentr,bcentr,rmaxis,zmaxis,simag,sibry,psi,current,fpol,pres,ffprim,pprime,
                   qpsi,psirz)

    return g
end


text:
## AS2008 model tests
@testset "AS2008 GMPE PGA" begin
  # init model parameters
  include("../examples/as2008.conf")
  # test at epicenter with M7.0, VS=350
  @test gmpe_as2008(eq_7,config_as2008_pga,grid_epicenter)[1].pga == 22.45
  # run PGA modeling on grid without minpga M6.0
  A = gmpe_as2008(eq_6,config_as2008_pga,grid)
  @test length(A) == TEST_GRID_SIZE
  @test round(sum([A[i].pga for i=1:length(A)]),digits=2) == 4.39
  # run PGA modeling on grid with minpga M6.0
  A = gmpe_as2008(eq_6,config_as2008_pga,grid,min_val=0.22)
  @test length(A) == WITH_MINPGA
  @test round(sum([A[i].pga for i=1:length(A)]),digits=2) == 2.86
  # run PGA modeling with M4.0
  A = gmpe_as2008(eq_4,config_as2008_pga,grid)
  @test length(A) == TEST_GRID_SIZE
  @test round(sum([A[i].pga for i=1:length(A)]),digits=2) == 0.21
  # run PGA modeling for plotting M6.0
  A = gmpe_as2008(eq_6,config_as2008_pga)
  @test length(A) == SIMULATION_ARRAY_SIZE
  @test round(sum(A),digits=2) == 573.5
  # run PGA Modeling M=4
  @test round(sum(gmpe_as2008(eq_4,config_as2008_pga)),digits=2) == 84.36
  # run PGA Modeling M=7
  @test round(sum(gmpe_as2008(eq_7,config_as2008_pga)),digits=2) == 1486.84
  # run PGA Modeling M=7 with VS30 = 900
  @test round(sum(gmpe_as2008(eq_7,config_as2008_pga,VS30=900)),digits=2) == 1031.1
  # run PGA Modeling M=7 with VS30 = 1600
  @test round(sum(gmpe_as2008(eq_7,config_as2008_pga,VS30=1600)),digits=2) == 817.23
end


text:
export hessian_dense
"""Compute the full hessian of an nlp (AbstractNLPModel). Should only be used
for lower dimensions problems"""
function hessian_dense(nlp, x)
    n = length(x)
    H = hess(nlp, x)
    tempH = (H + tril(H, -1)')
    H = Matrix(Symmetric(tempH, :L))
    return H
end


text:
if Sys.isunix() && !Sys.isbsd()
    function unpack(file,directory,extension,secondary_extension)
        if ((extension == ".gz" || extension == ".Z") && secondary_extension == ".tar") || extension == ".tgz"
            return (`tar xzf $file --directory=$directory`)
        elseif (extension == ".bz2" && secondary_extension == ".tar") || extension == ".tbz"
            return (`tar xjf $file --directory=$directory`)
        elseif extension == ".xz" && secondary_extension == ".tar"
            return pipeline(`unxz -c $file `, `tar xv --directory=$directory`)
        elseif extension == ".tar"
            return (`tar xf $file --directory=$directory`)
        elseif extension == ".zip"
            return (`unzip -x $file -d $directory`)
        elseif extension == ".gz"
            return pipeline(`mkdir $directory`, `cp $file $directory`, `gzip -d $directory/$file`)
        end
        error("I don't know how to unpack $file")
    end
end

if Sys.iswindows()
    const exe7z = joinpath(JULIA_HOME, "7z.exe")

    function unpack(file,directory,extension,secondary_extension)
        if ((extension == ".Z" || extension == ".gz" || extension == ".xz" || extension == ".bz2") &&
                secondary_extension == ".tar") || extension == ".tgz" || extension == ".tbz"
            return pipeline(`$exe7z x $file -y -so`, `$exe7z x -si -y -ttar -o$directory`)
        elseif (extension == ".zip" || extension == ".7z" || extension == ".tar" ||
                (extension == ".exe" && secondary_extension == ".7z"))
            return (`$exe7z x $file -y -o$directory`)
        end
        error("I don't know how to unpack $file")
    end
end


text:
struct RDateRange
    from::Dates.Date
    to::Union{Dates.Date, Nothing}
    period::RDate
    inc_from::Bool
    inc_to::Bool
    calendar_mgr::CalendarManager
end

function Base.iterate(iter::RDateRange, state=nothing)
    if state === nothing
        count = 0
        elem = apply(multiply_no_roll(iter.period, count), iter.from, iter.calendar_mgr)
        while elem > iter.from
            count -= 1
            elem = apply(multiply_no_roll(iter.period, count), iter.from, iter.calendar_mgr)
        end

        from_op = iter.inc_from ? Base.:< : Base.:<=
        while from_op(elem, iter.from)
            count += 1
            elem = apply(multiply_no_roll(iter.period, count), iter.from, iter.calendar_mgr)
        end

        state = (elem, count)
    end
    elem, count = state
    op = iter.inc_to ? Base.:> : Base.:>=
    if iter.to !== nothing && op(elem,iter.to)
        return nothing
    end

    return (elem, (apply(multiply_no_roll(iter.period, count+1), iter.from, iter.calendar_mgr), count+1))
end

Base.IteratorSize(::Type{RDateRange}) = Base.SizeUnknown()
Base.eltype(::Type{RDateRange}) = Dates.Date
function Base.range(from::Dates.Date, period::RDate; inc_from::Bool=true, inc_to::Bool=true, cal_mgr::Union{CalendarManager,Nothing}=nothing)
    return RDateRange(from, nothing, period, inc_from, inc_to, cal_mgr !== nothing ? cal_mgr : calendar_mgr)
end

function Base.range(from::Dates.Date, to::Dates.Date, period::RDate; inc_from::Bool=true, inc_to::Bool=true, cal_mgr::Union{CalendarManager,Nothing}=nothing)
    return RDateRange(from, to, period, inc_from, inc_to, cal_mgr !== nothing ? cal_mgr : calendar_mgr)
end


text:
using ParserCombinator

# subset of the grammar at http://www.tptp.org/TPTP/SyntaxBNF.html

global e = Env()

spc = Drop(Star(Space()))
@with_pre spc begin
  # useful sundry things
  lowerAlpha, upperAlpha, numeric = p"[a-z]", p"[A-Z]", p"[0-9]"
  alphaNumeric = lowerAlpha | upperAlpha | numeric | p"_"
  lowerWord = lowerAlpha + Star(alphaNumeric) |> x -> string(x...)
  upperWord = upperAlpha + Star(alphaNumeric) |> x -> string(x...)
  integer = PInt32()
  atomicWord = lowerWord
  name = atomicWord | integer

  # logic definitions
  functor = atomicWord > Functor
  constant = functor
  proposition = constant > ConstantTerm
  variable = upperWord > Variable
  fofFunctionTerm = Delayed()
  fofTerm = fofFunctionTerm | variable
  fofArguments = PlusList(fofTerm, E",")
  fofPlainTerm = (functor + P"\s*\(" + fofArguments + P"\s*\)" |>
                  x -> FunctionTerm(x[1],x[2:end])) | proposition
  fofFunctionTerm.matcher = fofPlainTerm
  fofPlainAtomicFormula = fofPlainTerm
  fofAtomicFormula = fofPlainAtomicFormula > x-> PredicateTerm(x.name,x.args)
  vLine = E"|"
  negatedAtomicFormula = E"~" + fofAtomicFormula |>
    x-> NegationTerm(PredicateTerm(x[1].name, x[1].args))
  dollarWord = p"\s*\$" + lowerWord |> x -> "\$" * x[2]
  atomicDefinedWord = dollarWord
  atomicSystemWord = p"\s*\$\$" + lowerWord
  definedFunctor = atomicDefinedWord
  definedConstant = definedFunctor
  systemFunctor = atomicSystemWord
  systemConstant = systemFunctor
  untypedAtom = constant | systemConstant
  atom = untypedAtom | definedConstant
  fofTerm = fofFunctionTerm | variable
  fofInfixUnary = fofTerm + P"\s*!=" + fofTerm |>
    x -> NegationTerm(equalTerm(x[1],x[2]))

  # CNF formulae
  literal = fofAtomicFormula | negatedAtomicFormula | fofInfixUnary
  disjunction = Delayed()
  disjunction.matcher = literal | (literal + vLine + disjunction)
  cnfFormula = (disjunction | (E"(" + disjunction + E")")) |> x -> OrTerm(x...)

  # FOF formulae
  unaryConnective = E"~"
  nonassocConnective = p"\s*<=>" | p"\s*=>" | p"\s*<=" | p"\s*<~>" | p"\s*~\|" |
    p"\s*~&"
  fofQuantifier = p"\s*!" | p"\s*\?"
  fofVariableList = PlusList(variable, E",")

  fofUnitaryFormula = Delayed()
  fofUnaryFormula = (unaryConnective + fofUnitaryFormula > x-> NegationTerm(x)) |
    fofInfixUnary
  fofUnitFormula = fofUnitaryFormula | fofUnaryFormula
  fofQuantifiedFormula = fofQuantifier + P"\s*\[" + fofVariableList + P"\s*\]\s*:" +
    fofUnitFormula |> x -> (x[1]=="!" ? AQuantifierTerm : EQuantifierTerm)(x[2:end]...)
  fofLogicFormula = Delayed()
  fofUnitaryFormula.matcher =  (P"\s*\(" + fofLogicFormula + P"\s*\)\s*") |
    (fofAtomicFormula + spc) | (fofQuantifiedFormula + spc)
  fofOrFormula = fofUnitFormula + vLine + PlusList(fofUnitFormula, vLine) |>
    x-> OrTerm(x...)
  fofAndFormula = fofUnitFormula + P"\s*&\s*" + PlusList(fofUnitFormula, P"\s*&\s*") |>
    x-> AndTerm(x...)
  fofBinaryAssoc = fofOrFormula | fofAndFormula
  fofBinaryNonassoc = fofUnitFormula + nonassocConnective + fofUnitFormula |>
    x-> begin
      if x[2] == "=>"
        OrTerm([NegationTerm(x[1]), x[3]])
      elseif x[2] == "<=>"
        AndTerm([OrTerm([NegationTerm(x[1]), x[3]]), OrTerm([NegationTerm(x[3]), x[1]])])
      elseif x[2] == "<="
        OrTerm([NegationTerm(x[3]), x[1]])
      else
        error("$(x[2]) not implemented")
      end
    end
  fofBinaryFormula = fofBinaryAssoc | fofBinaryNonassoc
  fofLogicFormula.matcher = fofUnitaryFormula | fofBinaryFormula | fofUnaryFormula
  fofFormula = fofLogicFormula

  # TFF: typed formula with FOOL extensions
  infixInequality = P"\s*!="
  definedInfixPred = P"\s*="
  untypedAtom = constant
  typeFunctor = atomicWord
  typeConstant = typeFunctor
  definedType = atomicDefinedWord
  tfxUnitaryFormula = variable
  tffAtomicType = typeConstant | definedType
  tffTypedVariable = variable + P"\s*:" + tffAtomicType |>
    x -> (e[x[1].name] = Environment.PrimitiveType(x[2]); x[1])
  tffVariable = tffTypedVariable | (variable > x -> (e[x.name] = Environment.InputType(); x))
  tffVariableList = PlusList(tffVariable, E",") |> x -> convert(Vector{Variable}, x)
  tffTerm = Delayed()
  tffArguments = PlusList(tffTerm, P"\s*,") |>
    x -> convert(Vector{typejoin(typeof.(x)...)}, x)
  tffPlainAtomic = (functor + P"\s*\(" + tffArguments + P"\s*\)\s*" |>
                    x -> begin
                    fType = get(e, x[1], ())
                    outputType = isempty(fType) ? Environment.BoolType() : fType[end]
                    if isempty(fType)
                    numArgs = length(x[2])
                    e[x[1].name] = ntuple(i -> i>numArgs ? Environment.BoolType() : outputType, numArgs+1)
                    end
                    (outputType != Environment.BoolType() ? FunctionTerm : PredicateTerm)(x[1], x[2])
                    end) |
      (constant > ConstantTerm)
  tffAtomicFormula = tffPlainAtomic
  tfxUnitaryFormula  = variable
  tffUnitaryFormula = Delayed()
  tffLogicFormula = Delayed()
  tffUnitaryTerm = tffAtomicFormula | variable |
    (P"\s*\(" + tffLogicFormula + P"\s*\)\s*")
  tffPreunitFormula = Delayed()
  tffPrefixUnary = unaryConnective + tffPreunitFormula > x -> NegationTerm(x)
  tffInfixUnary = tffUnitaryTerm  + infixInequality + tffUnitaryTerm |>
    x -> NegationTerm(equalTerm(x[1],x[2]))
  tffPreunitFormula.matcher = tffPrefixUnary | tffUnitaryFormula
  tffUnaryFormula = tffPrefixUnary | tffInfixUnary
  tffDefinedInfix = tffUnitaryTerm + definedInfixPred + tffUnitaryTerm |>
    x -> equalTerm(x[1],x[2])
  tffUnitFormula = tffUnitaryFormula | tffUnaryFormula | tffDefinedInfix
  tffQuantifiedFormula = fofQuantifier + P"\s*\[" + tffVariableList + P"\s*\]\s*:" + tffUnitFormula |>
    x -> (x[1]=="!" ? AQuantifierTerm : EQuantifierTerm)(x[2:end]...)
  tffUnitaryFormula.matcher = tffQuantifiedFormula | tffAtomicFormula |
    (P"\s*\(" + tffLogicFormula + P"\s*\)\s*") | tfxUnitaryFormula
  tffBinaryNonassoc = tffUnitFormula  + nonassocConnective + tffUnitFormula |>
    x -> begin
      if x[2] == "=>"
        OrTerm([NegationTerm(x[1]), x[3]])
      elseif x[2] == "<=>"
        AndTerm([OrTerm([NegationTerm(x[1]), x[3]]), OrTerm([NegationTerm(x[3]), x[1]])])
      elseif x[2] == "<="
        OrTerm([NegationTerm(x[3]), x[1]])
      else
        error("$(x[2]) not implemented")
      end
    end
  tffOrFormula = tffUnitFormula + vLine  + PlusList(tffUnitFormula, vLine) |>
    x -> OrTerm(x...)
  tffAndFormula = tffUnitFormula + P"\s*&"  + PlusList(tffUnitFormula, P"\s*&") |>
    x -> AndTerm(x...)
  tffBinaryAssoc = tffOrFormula | tffAndFormula

  tffUnitaryType = Delayed()
  tffXprodType = Delayed()
  tffXprodType.matcher = (tffUnitaryType + P"\s*\*" + tffAtomicType |> x -> x) |
    (tffXprodType + P"\s*\*" + tffAtomicType |> x -> [x[1]...,x[2]])
  tffUnitaryType.matcher = tffAtomicType | (P"\s*\(" + tffXprodType + P"\s*\)")
  tffMappingType = tffUnitaryType + P"\s*>" + tffAtomicType |>
    x -> isa(x[1], AbstractArray) ? [x[1]..., x[2]] : [x[1], x[2]]

  tffTopLevelType = Delayed()
  tffTopLevelType.matcher = tffMappingType | (P"\s*\(" + tffTopLevelType + P"\s*\)") |
    tffAtomicType

  tffAtomTyping = Delayed()
  tffAtomTyping.matcher = (untypedAtom + P"\s*:" + tffTopLevelType |>
                           x -> e[x[1].name] = Environment.PrimitiveType.(x[2]))|
                           (P"\s*\(" + tffAtomTyping + P"\s*\)")

  tffBinaryFormula = tffBinaryNonassoc | tffBinaryAssoc
  tffLogicFormula.matcher = tffUnitaryFormula | tffUnaryFormula | tffBinaryFormula |
    tffDefinedInfix
  tffTerm.matcher = tffLogicFormula
  tffSubtype = untypedAtom + P"\s*<<" + atom
  tffFormula = tffLogicFormula | tffAtomTyping | tffSubtype

  # top level
  formulaRole = e"axiom" | e"hypothesis" | e"definition" | e"assumption" |
    e"lemma" | e"theorem" | e"corollary" | e"conjecture" | e"negated_conjecture" |
    e"plain" | e"fi_domain" | e"fi_functors" | e"fi_predicates" | e"type" |
    e"unknown"
  cnfAnnotated = (P"\s*cnf\(" + name + P"\s*," + formulaRole + P"\s*," + cnfFormula + P"\s*\)\.\s*") |>
    x -> (Symbol(x[1]), x[3])
  fofAnnotated = (P"\s*fof\(" + name + P"\s*," + formulaRole + P"\s*," + fofFormula + P"\s*\)\.\s*") |>
    x -> (Symbol(x[1]), x[2] == "conjecture" ? NegationTerm(x[3]) : x[3])
  tffAnnotated = (P"\s*tff\(" + name + P"\s*," + formulaRole + P"\s*," + tffFormula + P"\s*\)\.\s*") |>
    x -> (Symbol(x[1]),  x[2] == "conjecture" ? NegationTerm(x[3]) : x[3])
  annotatedFormula = cnfAnnotated | fofAnnotated | tffAnnotated

  sqChar = p"([\40-\46\50-\133\135-\176]|[\\]['\\])"
  fileName = E"'" + Plus(sqChar) + E"'"
  includeRule = P"\s*include\(" + fileName + P"\s*\)\.\s*" > string
  commentLine = P"%([ -~])*"
  comment = commentLine
  tptpInput = annotatedFormula | (comment + spc) | includeRule

  tptpFile = Plus(tptpInput)
end

"""
    makeFOL!(s, e::Env)

After parsing an untyped tff input we convert it to standard first order format and modify
the environment accordingly. If `s` is already typed then do nothing.
"""
function makeFOL!(v::Variable, e::Env, outputIsBool=true)
  haskey(e, v.name) || (e[v.name] = (Environment.InputType(),))
  v
end
function makeFOL!(f::FunctionTerm, e::Env, outputIsBool=false)
  if !haskey(e,f.name.name)
    e[f.name.name] = ntuple(_ -> Environment.InputType(), length(f.args)+1)
  end
  FunctionTerm(f.name, makeFOL!.(f.args, Ref(e), false))
end
function makeFOL!(p::PredicateTerm, e::Env, outputIsBool=true)
  numArgs = length(p.args)
  if outputIsBool
    if !haskey(e, p.name.name)
      e[p.name.name] = ntuple(i -> i>numArgs ? Environment.BoolType() : Environment.InputType(), numArgs+1)
    end
    PredicateTerm(p.name, makeFOL!.(p.args, Ref(e), false))
  else
    if !haskey(e, p.name.name)
      e[p.name.name] = ntuple(_ -> Environment.InputType(), numArgs+1)
    end
    FunctionTerm(p.name, makeFOL!.(p.args, Ref(e), false))
  end
end
makeFOL!(n::NegationTerm, e::Env, outputIsBool=true) =
  NegationTerm(makeFOL!(n.scope, e, true))
makeFOL!(j::JunctionTerm, e::Env, outputIsBool=true) =
  rootType(j)(makeFOL!.(j.juncts, Ref(e), true))
makeFOL!(q::QuantifierTerm, e::Env, outputIsBool=true) =
  rootType(q)(q.variables, makeFOL!(q.scope, e, true))

"""
    parseTPTP(io::IO), parseTPTP(fileName::AbstractString)

Read in the TPTP description at the given IO stream
"""
function parseTPTP(io::IO)
  resetEnv!(FirstOrderLogic.e)  # create a fresh environment
  ret = parse_one(read(io, String), tptpFile + Eos())
  close(io)
  ret
end
parseTPTP(fileName::AbstractString) = parseTPTP(open(fileName,"r"))

"""
    @fol_str(s)

Build the data structure representing the first order formula stored in TPTP
format in the string s. Note that \$ should not be escaped inside the fol_str
macro.
"""
macro fol_str(s)
  resetEnv!(FirstOrderLogic.e)  # create a fresh environment
  makeFOL!(simplify(parse_one(s, tffFormula+Eos())[1]), FirstOrderLogic.e)
end

"""
    @tptp_str(s)

Build the data structure representing the TPTP directive(s) in the string s.
Note that \$ should not be escaped inside the fol_str macro
"""
macro tptp_str(s)
  resetEnv!(e)  # create a fresh environment
  parse_one(s, tptpFile+Eos())[1]
end


text:
function FactoryManagedFrame(arg0::Frame, arg1::TransformProvider, arg2::jboolean, arg3::Predefined)
    return FactoryManagedFrame((Frame, TransformProvider, jboolean, Predefined), arg0, arg1, arg2, arg3)
end

function get_factory_key(obj::FactoryManagedFrame)
    return jcall(obj, "getFactoryKey", Predefined, ())
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function add1!(x,y)
    x[y] .+= 1
end

function devec_add1!(x,y)
    for i=1:length(y)
        x[y[i]] += 1
    end
end

function devec_add1_logical!(x,y)
    for i=1:length(y)
        if y[i]
            x[i] += 1
        end
    end
end


text:
function t_roc()

	# Generate log-likelihood scores; classes fully overlap - the only one to guarantee that
	# the calculated measure actually work
	tar = collect(-1:0.01:1);
	non = -tar;

	x=(1+[tar;non])./2;
	X=[x';(1-x)'];
	y=[fill("a",length(tar));fill("b",length(non))];

	# Define functions to verify that ROC optimization works
	calculate_measure(y::AbstractVector{T}, yest::AbstractVector{T}, class::T, measure::j4pr.ROC.TPr) where T =
		sum((y.==class) .& (yest.==class)) ./ (sum(y.==class)+eps())

	calculate_measure(y::AbstractVector{T}, yest::AbstractVector{T}, class::T, measure::j4pr.ROC.TNr) where T =
		sum((y.!=class) .& (yest.!=class)) ./ (sum(y.!=class)+eps())

	calculate_measure(y::AbstractVector{T}, yest::AbstractVector{T}, class::T, measure::j4pr.ROC.FPr) where T =
		sum((y.!=class) .& (yest.==class)) ./ (sum(y.!=class)+eps())

	calculate_measure(y::AbstractVector{T}, yest::AbstractVector{T}, class::T, measure::j4pr.ROC.FNr) where T =
		sum((y.==class) .& (yest.!=class)) ./ (sum(y.==class)+eps())

	perfmetric = [j4pr.ROC.TPr(), j4pr.ROC.FPr(), j4pr.ROC.TNr(), j4pr.ROC.FNr()];
	desiredval = collect(0:0.1:1);
	tol = 0.01
	class = ["a","b"];
	methods=[:j4pr,:ra]
	yu = sort(unique(y))
	for c in class
		for pm in perfmetric
			for dv in desiredval
				for m in methods
					op = j4pr.ROC.findop(X, y, c, pm, dv, m, 10_000); # large number of ops for low tolerances
					Xop = op.weights*X;
					yest = yu[j4pr.targets(indmax,Xop)];
					ev = calculate_measure(y, yest, c, pm)
					Test.@test abs(dv-ev) <=tol
				end
			end
		end
	end

	op = j4pr.ROC.findop(X, y, "a", j4pr.ROC.FPr(), 0.1);

	# Test changing an op
	pos = 2
	j4pr.ROC.changeop!(op,pos)
	Test.@test op.pos == pos

	j4pr.ROC.changeop!(op,j4pr.ROC.FPr(),0.0)
	Test.@test op.pos == size(op.rocdata,1)

	# Test simple op
	sop = j4pr.ROC.simpleop(op); Test.@test sop.weights == op.weights
	sop = j4pr.ROC.simpleop([1,2,3]); Test.@test sop.weights == diagm([1.0,2.0,3.0])

	# Test for the FunctionCell interface
	D=j4pr.DataGenerator.iris();
	(tr,ts)=j4pr.splitobs(j4pr.shuffleobs(D),0.7);
	w=j4pr.libsvm();
	wt=w(tr);
	out=tr|>wt;

	# Test op search
	Test.@test try
		uop=j4pr.findop("virginica",j4pr.ROC.FPr(),0.3);
		cop = uop(out);
		pt=wt + cop;
		+D |> pt;
		true
	catch
		false
	end

	# Test changing an op
	pos = 2
	j4pr.changeop!(cop,pos)
	Test.@test cop.x.data.pos == pos

	j4pr.changeop!(cop,j4pr.ROC.FPr(),0.0)
	Test.@test cop.x.data.pos == size(cop.x.data.rocdata,1)

	# Test simple op
	uop=j4pr.findop("virginica",j4pr.ROC.FPr(),0.3);
	cop = uop(out)
	sop = j4pr.simpleop(cop)
	Test.@test sop.x.data.weights == cop.x.data.weights

	sop = j4pr.simpleop([1,2,3]);
	Test.@test sop.x.data.weights == diagm([1.0,2.0,3.0])

end


text:
# The normal fit api in Distributions.jl is:
#   fit(Distr, values)
# so we want to do a more general
#   fitall(values)
# that will try all distributions that can be fitted and return back the "best" one
# for relevant criteria of "best", typically BIC or AIC.

using Distributions

AllContinuousDistributions = [
    Arcsine, Beta, BetaPrime, Cauchy, Chi, Chisq, Erlang, Exponential,
    FDist, Frechet, Gamma, GeneralizedExtremeValue, GeneralizedPareto,
    Gumbel, InverseGamma, InverseGaussian, Laplace, Levy, LogNormal,
    Logistic, Normal, NormalInverseGaussian, Pareto, Rayleigh,
    SymTriangularDist, TDist, TriangularDist, Uniform, VonMises, Weibull
]

samples1 = rand(Normal(3.4, 1.2), 100)

ContinuousDistributionsThatCanFit = filter(AllContinuousDistributions) do D
    try
        fit_mle(D, samples1)
        true
    catch _err
        false
    end
end

# Start for continuous distributions
function fitalldistr{T<:Real}(x::AbstractArray{T}, sortby = :BIC)
    n = length(x)
    allfits = map(ContinuousDistributionsThatCanFit) do D
        d = fit_mle(D, x)
        loglik = loglikelihood(d, x)
        k = length(params(d))
        if sortby == :BIC
            criteria = -2 * loglik + k * log(n)
        elseif sortby == :AIC
            criteria = -2 * loglik + 2 * k
        elseif sortby == :AICc
            criteria = -2 * loglik + 2 * k + ((2*k*(k+1))/(n-k-1))
        elseif sortby == :NLogL
            criteria = -loglik
        end
        (d, criteria)
    end
    sort(allfits, by = t -> t[2])
end

function fitall{T<:Real}(x::AbstractArray{T}, sortby = :BIC)
    fitalldistr(x, sortby)[1][1]
end

samples1 = rand(Normal(3.4, 1.2), 100)
d1 = fitall(samples1)

samples2 = rand(Exponential(0.3), 100)
d2 = fitall(samples2)

samples3 = rand(LogNormal(-5.0, 2.0), 100)
d3s = fitalldistr(samples3)
d3 = fitall(samples3)

samples4 = rand(Gamma(6.20, 1.3), 200)
d4s = fitalldistr(samples4)
d4 = fitall(samples4)

# Should we maybe use a non-parametric test to select between them?


text:
using BenchmarkTools
using MRIgeneralizedBloch
using Test

##
τmin = 0
τmax = 100

##
Gfull = greens_superlorentzian

print("Time to pre-calculate the interpolation:")
Ginte = @btime interpolate_greens_function(greens_superlorentzian, τmin, τmax)

##
for i = 1:100
    τ = τmin + (τmax - τmin) * rand()
    @test Gfull(τ) ≈ Ginte(τ)
end

## benchmark
τ = τmin + (τmax - τmin) * rand()
print("Time to evaluate the original     super-Lorenzian Green's function:")
@btime Gfull($τ)

print("Time to evaluate the interpolated super-Lorenzian Green's function:")
@btime Ginte($τ)

## #####################################################
# Test derivatives
########################################################
t = 100e-6
τ = 0
T2s = 10e-6
δT2s = eps()

## Lorentzian
dGdT2s = dG_o_dT2s_x_T2s_lorentzian((t-τ)/T2s)/T2s
dGdT2s_fd = (greens_lorentzian((t-τ)/(T2s + δT2s)) - greens_lorentzian((t-τ)/T2s)) / δT2s
@test dGdT2s ≈ dGdT2s_fd rtol = 1e-6

## Gaussian
dGdT2s = dG_o_dT2s_x_T2s_gaussian((t-τ)/T2s)/T2s
dGdT2s_fd = (greens_gaussian((t-τ)/(T2s + δT2s)) - greens_gaussian((t-τ)/T2s)) / δT2s
@test dGdT2s ≈ dGdT2s_fd rtol = 1e-6

## super-Lorentzian
dGdT2s = dG_o_dT2s_x_T2s_superlorentzian((t-τ)/T2s)/T2s
dGdT2s_fd = (greens_superlorentzian((t-τ)/(T2s + δT2s)) - greens_superlorentzian((t-τ)/T2s)) / δT2s
@test dGdT2s ≈ dGdT2s_fd rtol = 1e-4


text:
module bad

import Base.show, Base.isless, Base.isequal, Base.-, Base.+, Base.*, Base.string,
    Base.convert, Base.<=, Base.>=, Base.|, Base.length, Base.size

import Printf.@printf, Printf.@sprintf

import ProgressMeter

import SpecialFunctions.gamma, SpecialFunctions.beta_inc

import Distributions, Distributions.cdf, Distributions.pdf

import QuadGK.quadgk, QuadGK.gauss

import Roots

using JuMP, GLPK, Ipopt

import DataFrames

import Statistics.mean

include("priors/Prior.jl")
export is_proper, condition, update, pdf, cdf, mean, expectation, bounds

include("util.jl")
export valid

include("designs.jl")
export Design, OptimalDesign,
    n1, n2, n, c2,
    early_futility, early_efficacy, continuation_region, futility_region,
    efficacy_region, early_stop_region,
    as_table, reject, sample_space

include("priors/JeffreysPrior.jl")
export JeffreysPrior

include("priors/GenericDistribution.jl")
export GenericDistribution

include("priors/Beta.jl")
export Beta, BetaMixture

include("priors/PointMass.jl")
export PointMass

include("pmf.jl")
export pmf, pmf_x1, pmf_x2_given_x1, cdf, cdf_x1, cdf_x2_given_x1

include("estimators/Estimator.jl")
export estimate, bias, mean_squared_error, mean_absolute_error

include("estimators/MaximumLikelihoodEstimator.jl")
export MaximumLikelihoodEstimator

include("estimators/PosteriorMean.jl")
export PosteriorMean, PosteriorMeanPrecalculated

include("estimators/RaoBlackwellEstimator.jl")
export RaoBlackwellEstimator

include("estimators/CompatibleMLE.jl")
export CompatibleMLE

include("orderings/Ordering.jl")
export EstimatorOrdering, more_extreme, p_value,
    compatible, mlecompatible

include("orderings/PValue.jl")
export PValue, evaluate

include("interval-estimators/IntervalEstimator.jl")
export IntervalEstimator, get_bounds, coverage_probability, mean_width,
    compatible

include("interval-estimators/ClopperPearsonInterval.jl")
export ClopperPearsonInterval

include("interval-estimators/PosteriorCredibleInterval.jl")
export PosteriorCredibleInterval

include("Score.jl")
export update!,
    SampleSize, Power, TypeOneErrorRate,
    CompositeScore

include("Problem.jl")
export Problem, optimise, adapt

include("objectives/Objective.jl")
export ScoreObjective, minimise, maximise

include("objectives/MiniMaxSampleSize.jl")
export MiniMaxSampleSize

include("constraints.jl")
export PowerConstraint, TypeOneErrorRateConstraint, conditional

end # module


text:
using LinearAlgebra

const libblas = Base.libblas_name
const liblapack = Base.liblapack_name
import LinearAlgebra.BLAS.@blasfunc
import LinearAlgebra: BlasFloat, BlasInt, LAPACKException,
    DimensionMismatch, chkstride1, Givens

include("preprocWrapper.jl")

function dtgsja!(A::AbstractMatrix{Float64}, B::AbstractMatrix{Float64},
    U::AbstractMatrix{Float64}, V::AbstractMatrix{Float64}, Q::AbstractMatrix{Float64},
    k::Integer, l::Integer)
    # chkstride1(A, B)
    jobu = 'U'
    jobv = 'V'
    jobq = 'Q'
    m, n = size(A)
    p = size(B,1)
    tola = max(m, n) * opnorm(A, 1) * eps(Float64)
    tolb = max(p, n) * opnorm(B, 1) * eps(Float64)
    lda = max(1, stride(A, 2))
    ldb = max(1, stride(B, 2))
    ldu = max(1, m)
    ldv = max(1, p)
    ldq = max(1, n)
    alpha = Vector{Float64}(undef, n)
    beta = Vector{Float64}(undef, n)
    work = Vector{Float64}(undef, 2*n)
    ncycle  = Ref{BlasInt}()
    info  = Ref{BlasInt}()

    # for i in 1:2  # first call returns lwork as work[1]
    ccall((@blasfunc(dtgsja_), liblapack), Cvoid,
    (Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
    Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
    Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
    Ref{Float64}, Ref{Float64},
    Ptr{Float64}, Ptr{Float64},
    Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
    Ptr{Float64}, Ref{BlasInt}, Ref{BlasInt}
    ),
    jobu, jobv, jobq,
    m, p, n, k, l,
    A, lda, B, ldb,
    tola, tolb,
    alpha, beta,
    U, ldu, V, ldv, Q, ldq,
    work, ncycle, info)
    # end
    return A, B, alpha, beta, U, V, Q, work, ncycle, info
end

function gsvdwrapper(A::AbstractMatrix{Float64}, B::AbstractMatrix{Float64})
    m, n = size(A)
    p = size(B)[1]
    # U, V, Q, k, l, A, B = dggsvp3!(A, B)
    # profiling
    ans1 = @timed dggsvp3!(A, B)
    U = ans1[1][1]
    V = ans1[1][2]
    Q = ans1[1][3]
    k = ans1[1][4]
    l = ans1[1][5]
    A = ans1[1][6]
    B = ans1[1][7]

    # A, B, alpha, beta, U, V, Q, work, ncycle, info = dtgsja!(A, B, U, V, Q, k, l)
    # profiling
    ans2 = @timed dtgsja!(A, B, U, V, Q, k, l)
    A = ans2[1][1]
    B = ans2[1][2]
    alpha = ans2[1][3]
    beta = ans2[1][4]
    U = ans2[1][5]
    V = ans2[1][6]
    Q = ans2[1][7]

    R0 = zeros(k+l, n)
    D1 = Matrix(1.0I, m, k+l)
    D2 = zeros(p, k+l)
    if m-k-l >= 0
        @views R0[:, n-k-l+1:n]  = A[1:k+l, n-k-l+1:n]
        @views D1[k+1:k+l, k+1:k+l] = Diagonal(alpha[k+1:k+l])
        @views D2[1:l, k+1:k+l] = Diagonal(beta[k+1:k+l])
    else
        # println("$k $l")
        @views R0[1:m, n-k-l+1:n] = A[1:m, n-k-l+1:n]
        @views R0[m+1:k+l, m+n-k-l+1:n] = B[m-k+1:l, n+m-k-l+1:n]
        @views D1[k+1:m, k+1:m] = Diagonal(alpha[k+1:m])
        @views D2[1:m-k, k+1:m] = Diagonal(beta[k+1:m])
        @views D2[m-k+1:l, m+1:k+l] = Matrix(1.0I, k+l-m, k+l-m)
    end
    # return U, V, Q, D1, D2, k, l, R0
    # profiling
    return ans1[2], ans2[2]
end

# function compute_dot(DX::Vector{Float64}, DY::Vector{Float64})
#     @assert length(DX) == length(DY)
#     n = length(DX)
#     incx = incy = 1
#     ccall((@blasfunc(ddot_), libblas), Float64,
#       (Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
#       Ptr{Float64}, Ref{BlasInt}),
#       n, DX, incx, DY, incy)
# end
#
# function compute_dot_julia(dx::Vector{Float64}, dy::Vector{Float64})
#     # ddot = 0.0
#     dtemp = 0.0
#     n = length(dx)
#     m = n % 5
#     if m != 0
#         for i = 1:m
#             dtemp = dtemp + dx[i]*dy[i]
#         end
#         if n < 5
#             # ddot = dtemp
#             return dtemp
#         end
#     end
#     mp1 = m + 1
#     for i = mp1:5:n-5
#         dtemp = dtemp + dx[i]*dy[i] + dx[i+1]*dy[i+1] +
#               dx[i+2]*dy[i+2] + dx[i+3]*dy[i+3] + dx[i+4]*dy[i+4]
#     end
#     # ddot = dtemp
#     # return ddot
#     return dtemp
# end


text:
# not used!

function draw_lorenz( pre::Vector, post::Vector )
    np = size(pre)[1]
    step = 1/np
    xr = 0:step:1
    # insert a zero element
    v1 = insert!(copy(pre),1,0)
    v2 = insert!(copy(post),1,0)

    layout = Layout(
        title="Lorenz Curve",
        xaxis_title="Population Share",
        yaxis_title="Net Income Share",
        xaxis_range=[0, 1],
        yaxis_range=[0, 1],
        legend=attr(x=0.01, y=0.95),
        width=700,
        height=700)

    pre = scatter(
        x=xr,
        y=v1,
        mode="line",
        name="Post" )

    post = scatter(
        x=xr,
        y=v2,
        mode="line",
        name="Post" )

    diag = scatter(y=[0,1], x=[0,1], showlegend=false, name="")
    ##return [pre, post, diag]
    #return Dict([:data=>[pre, post, diag],:layout=>layout])
    PlotlyJS.Plot( [pre, post, diag], layout)
end

function drawDeciles( pre::Vector, post :: Vector )
    v = pre-post;
    PlotlyJS.plot(
     bar(
        x=1:10,
        y=v )

end

function gain_lose_table( gl :: NamedTuple )
    lt = sum( gl.losers )
    gt = sum( gl.gainers )
    popn = sum( gl.popn )
    nct = sum( gl.nc )

    losers = md_format(lt)
    gainers = md_format(gt)
    nc = md_format(nct)

    losepct = md_format(100*lt/popn)
    gainpct = md_format(100*gt/popn)
    ncpct = md_format(100*nct/popn)

    tab = table(
        # header_values=["", "", "%"],
        cells_values=[
            ["Gainers","No Change","Losers"],
            [gainers,nc,losers],
            ["($(gainpct)%)","($(ncpct)%)","($(losepct))"]
        ]
    )
    return tab #PlotlyJS.plot( tab,
        # Layout(width=200, height=300))
end


text:
using Behavior
using Niancat
using Niancat.Gameface
using Niancat.Scores

column(t::Behavior.Gherkin.DataTable) = [row[1] for row in t]

@when("a solution board notification is sent") do context
    game = context[:game]
    round = gameround(game)
    instanceid = gameinstanceid(game)

    # Note that this step assumes that the Game instance has a field
    # `gameservice`, which isn't necessarily the case always.
    # Should we create a getter for gameservice?
    gameservice = game.gameservice

    notify!(gameservice, SolutionboardNotification(instanceid, round))
end

@when("a score board notification is sent") do context
    game = context[:game]
    round = gameround(game)
    instanceid = gameinstanceid(game)

    # Note that this step assumes that the Game instance has a field
    # `gameservice`, which isn't necessarily the case always.
    # Should we create a getter for gameservice?
    gameservice = game.gameservice

    notify!(gameservice, ScoreboardNotification(instanceid, round))
end

@then("a notification is sent, containing words") do context
    # This is a fake HTTP client, where we can have a look at what
    # is actually sent to the network.
    httpclient = context[:httpclient]

    # These are the words we expect in some notification
    words = column(context.datatable)

    # Collect all text bodies sent
    bodies = [body for (_uri, body) in httpclient.posts]

    # Check that the expected words are all in _one of the bodies_.
    found_words = false

    for body in bodies
        found_words = all(w -> contains(body, w), words) || found_words
    end

    if !found_words
        println("Expected words: $(words)")
        println("Actual notifications: $(bodies)")
        @fail "Could not find the expected words in any notification"
    end
end