text:
using ParameterHandling: Positive, Bounded
using ParameterHandling: vec_to_tril, tril_to_vec

mvnormal(args...) = MvNormal(args...)
pdiagmat(args...) = PDiagMat(args...)

@testset "parameters" begin
    @testset "postive" begin
        @testset "$val" for val in [5.0, 0.001f0, 1.0e-7]
            p = positive(val)
            test_parameter_interface(p)
            @test value(p) ≈ val
            @test typeof(value(p)) === typeof(val)
        end

        # Test edge cases around the size of the value relative to the error tol.
        @test_throws ArgumentError positive(-0.1)
        @test_throws ArgumentError positive(1e-11)
        @test value(positive(1e-11, Bijectors.Exp(), 1e-12)) ≈ 1e-11
    end

    @testset "bounded" begin
        @testset "$val" for val in [-0.05, -0.1 + 1e-12, 2.0 - 1e-11, 2.0 - 1e-12]
            p = bounded(val, -0.1, 2.0)
            test_parameter_interface(p)
            @test value(p) ≈ val
        end

        @test_throws ArgumentError bounded(-0.05, 0.0, 1.0)
    end

    @testset "fixed" begin
        val = (a=5.0, b=4.0)
        p = fixed(val)
        test_parameter_interface(p)
        @test value(p) == val
    end

    @testset "deferred" begin
        test_parameter_interface(deferred(sin, 0.5); check_inferred=tuple_infers)
        test_parameter_interface(deferred(sin, positive(0.5)); check_inferred=tuple_infers)
        test_parameter_interface(
            deferred(
                mvnormal, fixed(randn(5)), deferred(pdiagmat, positive.(rand(5) .+ 1e-1))
            );
            check_inferred=tuple_infers,
        )
    end

    @testset "orthogonal" begin
        is_almost_orthogonal(X::AbstractMatrix, tol) = norm(X'X - I) < tol

        @testset "nearest_orthogonal_matrix($T)" for T in [Float64, ComplexF64]
            X_orth = ParameterHandling.nearest_orthogonal_matrix(randn(T, 5, 4))
            @test is_almost_orthogonal(X_orth, 1e-9)
            X_orth_2 = ParameterHandling.nearest_orthogonal_matrix(X_orth)
            @test X_orth ≈ X_orth_2 # nearest_orthogonal_matrix is a projection.
        end

        X = orthogonal(randn(5, 4))
        @test X == X
        test_parameter_interface(X)
        @test is_almost_orthogonal(value(X), 1e-9)

        # We do not implement any custom rrules, so we only check that `Zygote` is able to
        # differentiate, and assume that the result is correct if it doesn't error.
        @testset "Zygote" begin
            _, pb = Zygote.pullback(X -> value(orthogonal(X)), randn(3, 2))
            @test only(pb(randn(3, 2))) isa Matrix{<:Real}
        end
    end

    @testset "positive_definite" begin
        @testset "vec_tril_conversion" begin
            X = tril!(rand(3, 3))
            @test vec_to_tril(tril_to_vec(X)) == X
            @test_throws ErrorException tril_to_vec(rand(4, 5))
        end
        X_mat = ParameterHandling.A_At(rand(3, 3)) # Create a positive definite object
        X = positive_definite(X_mat)
        @test X == X
        @test value(X) ≈ X_mat
        @test isposdef(value(X))
        @test vec_to_tril(X.L) ≈ cholesky(X_mat).L
        @test_throws ArgumentError positive_definite(rand(3, 3))
        test_parameter_interface(X)

        x, re = flatten(X)
        Δl = first(
            Zygote.gradient(x) do x
                X = re(x)
                return logdet(value(X))
            end,
        )
        ΔL = first(
            Zygote.gradient(vec_to_tril(X.L)) do L
                return logdet(L * L')
            end,
        )
        @test vec_to_tril(Δl) == tril(ΔL)
        ChainRulesTestUtils.test_rrule(vec_to_tril, x)
    end

    function objective_function(unflatten, flat_θ::Vector{<:Real})
        θ = value(unflatten(flat_θ))
        return abs2(θ.a) + abs2(θ.b)
    end

    # This is more of a worked example.
    @testset "Integration" begin
        θ0 = (a=5.0, b=4.0)
        flat_parameters, unflatten = flatten(θ0)

        results = Optim.optimize(
            θ -> objective_function(unflatten, θ),
            θ -> only(Zygote.gradient(θ -> objective_function(unflatten, θ), θ)),
            flat_parameters,
            BFGS(),
            Optim.Options();
            inplace=false,
        )

        # Check that it's successfully optimised.
        @test objective_function(unflatten, results.minimizer) < 1e-12
    end

    @testset "Other Integration" begin
        θ0 = (a=5.0, b=fixed(4.0))
        flat_parameters, unflatten = flatten(θ0)

        results = Optim.optimize(
            θ -> objective_function(unflatten, θ),
            θ -> only(Zygote.gradient(θ -> objective_function(unflatten, θ), θ)),
            flat_parameters,
            BFGS(),
            Optim.Options();
            inplace=false,
        )

        # Check that it's successfully optimised.
        @test value(unflatten(results.minimizer).b) == 4.0
    end

    function other_objective_function(unflatten, flat_θ::Vector{<:Real})
        X = value(unflatten(flat_θ))
        return -logpdf(X, -1.0) - logpdf(X, 1.0)
    end

    @testset "Normal" begin
        θ0 = deferred(Normal, randn(), positive(1.0))
        flat_parameters, unflatten = flatten(θ0)

        results = Optim.optimize(
            θ -> other_objective_function(unflatten, θ),
            θ -> only(Zygote.gradient(θ -> other_objective_function(unflatten, θ), θ)),
            flat_parameters,
            BFGS(),
            Optim.Options();
            inplace=false,
        )

        # Check that it's successfully optimised.
        @test mean(value(unflatten(results.minimizer))) ≈ 0 atol = 1e-7
    end

    @testset "value_flatten" begin
        x = (ones(3), fixed(5.0), (a=fixed(5.0), b=[6.0, 2.1]))
        v, unflatten = value_flatten(x)

        @test length(v) == 5
        @test unflatten(v) == ParameterHandling.value(x)
    end
end


text:
module DeepQLearning

using Random
using StatsBase
using Printf
using Parameters
using Flux
using BSON
using POMDPs
using POMDPModelTools
using POMDPPolicies
using POMDPLinter
using RLInterface
using LinearAlgebra
using TensorBoardLogger: TBLogger, log_value
using EllipsisNotation
using Distributions

export DeepQLearningSolver,
       AbstractNNPolicy,
       NNPolicy,
       getnetwork,
       resetstate!,
       DQExperience,
       restore_best_model,
       PrioritizedReplayBuffer,
       EpisodeReplayBuffer,

       # helpers
       flattenbatch,
       huber_loss,
       isrecurrent,
       batch_trajectories

include("helpers.jl")
include("policy.jl")
include("exploration_policy.jl")
include("evaluation_policy.jl")
include("prioritized_experience_replay.jl")
include("episode_replay.jl")
include("dueling.jl")
include("solver.jl")

end # module DeepQLearning

ext: jl


text:
richprint(io::IO, dom::Vector{Node}, x::ANY) = richprint_default(io, dom, x)

function richprint_default(io::IO, dom::Vector{Node}, x::ANY)
    t = typeof(x)
    isbits(t) && Base.show_circular(io, x) && return
    let id = startNode(dom, io, x, "")
        if which(Base.show, (typeof(io), typeof(x))).sig === Tuple{typeof(Base.show), IO, ANY}
            # override Base.default_show behavior
            nf = nfields(t)
            nb = sizeof(x)
            if nf != 0 || nb == 0
                Base.show(io, t)
                print(io, '(')
                recur_io = IOContext(io, :SHOWN_SET => x)
                for i in 1:nf
                    nextChild(dom, io, id)
                    f = fieldname(t, i)
                    if !isdefined(x, f)
                        print(io, undef_ref_str)
                    else
                        richprint(recur_io, dom, getfield(x, f))
                    end
                    if i < nf
                        print(io, ", ")
                    end
                end
                print(io, ')')
            else
                Base.show_default(io, x)
            end
        else
            Base.show(io, x)
        end
        endNode(dom, io, id)
    end
end

function richprint(io::IO, dom::Vector{Node}, x::Union)
    let id = startNode(dom, io, x, "Union")
        print(io, "Union")
        sorted_types = sort!(uniontypes(x); by=string)
        show_delim_array(io, dom, id, sorted_types, '{', ',', '}', false)
        endNode(dom, io, id)
    end
end

function richprint(io::IO, dom::Vector{Node}, x::UnionAll)
    let id = startNode(dom, io, x, "UnionAll")
        if Base.print_without_params(x)
            return Base.show(io, Base.unwrap_unionall(x).name)
        end
        richprint(IOContext(io, :unionall_env => x.var), dom, x.body)
        print(io, " where ")
        richprint(io, dom, x.var)
        endNode(dom, io, id)
    end
end

richprint(io::IO, dom::Vector{Node}, x::DataType) = richprint_datatype(io, dom, x)

function richprint_datatype(io::IO, dom::Vector{Node}, x::DataType)
    let id = startNode(dom, io, x, "DataType")
        istuple = x.name === Tuple.name
        if (!isempty(x.parameters) || istuple) && x !== Tuple
            n = length(x.parameters)

            # Print homogeneous tuples with more than 3 elements compactly as NTuple{N, T}
            if istuple && n > 3 && all(i -> (x.parameters[1] === i), x.parameters)
                print(io, "NTuple{", n, ", ")
                richprint(io, dom, x.parameters[1])
                print(io, "}")
            else
                Base.show(io, x.name)
                # Do not print the type parameters for the primary type if we are
                # printing a method signature or type parameter.
                # Always print the type parameter if we are printing the type directly
                # since this information is still useful.
                print(io, '{')
                for (i, p) in enumerate(x.parameters)
                    nextChild(dom, io, id)
                    richprint(io, dom, p)
                    i < n && print(io, ", ")
                end
                print(io, '}')
            end
        else
            Base.show(io, x.name)
        end
        endNode(dom, io, id)
    end
end

function richprint(io::IO, dom::Vector{Node}, p::Pair)
    let id = startNode(dom, io, p, "Pair")
        print(io, "(")
        nextChild(dom, io, id)
        richprint(io, dom, p.first)
        print(io, " => ")
        nextChild(dom, io, id)
        richprint(io, dom, p.second)
        print(io, ")")
        endNode(dom, io, id)
    end
end

function richprint(io::IO, dom::Vector{Node}, l::Core.MethodInstance)
    let id = startNode(dom, io, l, "MethodInstance")
        if isdefined(l, :def)
            if l.def.isstaged && l === l.def.generator
                print(io, "(@generated thunk)")
                richprint(io, dom, l.def)
            else
                richprint_tuple_as_call(io, dom, id, l.def.name, l.specTypes)
            end
        else
            print(io, "(Toplevel thunk)")
        end
        endNode(dom, io, id)
    end
end

function richprint(io::IO, dom::Vector{Node}, src::CodeInfo)
    # Fix slot names and types in function body
    lambda_io = IOContext(io, :SOURCEINFO => src)
    if src.slotnames !== nothing
        lambda_io = IOContext(io, :SOURCE_SLOTNAMES => Base.sourceinfo_slotnames(src))
    end
    let id = startNode(dom, io, src, "CodeInfo")
        print(io, "CodeInfo(")
        body = Expr(:body)
        body.args = src.code
        richprint(lambda_io, dom, body)
        print(io, ")")
        endNode(dom, io, id)
    end
    nothing
end

function show_delim_array(io::IO, dom::Vector{Node}, id::Int,
                          itr::Union{AbstractArray,SimpleVector}, op, delim, cl,
                          delim_one, i1=first(linearindices(itr)), l=last(linearindices(itr)))
    print(io, op)
    if !show_circular(io, itr)
        recur_io = IOContext(io, :SHOWN_SET => itr)
        if !haskey(io, :compact)
            recur_io = IOContext(recur_io, :compact => true)
        end
        first = true
        i = i1
        if l >= i1
            while true
                nextChild(dom, io, id)
                if !isassigned(itr, i)
                    print(io, undef_ref_str)
                else
                    x = itr[i]
                    richprint(recur_io, dom, x)
                end
                i += 1
                if i > l
                    delim_one && first && print(io, delim)
                    break
                end
                first = false
                print(io, delim)
                print(io, ' ')
            end
        end
    end
    print(io, cl)
end

function show_delim_array(io::IO, dom::Vector{Node}, id::Int, itr, op, delim, cl, delim_one, i1=1, n=typemax(Int))
    print(io, op)
    if !show_circular(io, itr)
        recur_io = IOContext(io, :SHOWN_SET => itr)
        state = start(itr)
        first = true
        while i1 > 1 && !done(itr, state)
            _, state = next(itr, state)
            i1 -= 1
        end
        if !done(itr, state)
            while true
                x, state = next(itr, state)
                nextChild(dom, io, id)
                richprint(recur_io, dom, x)
                i1 += 1
                if done(itr, state) || i1 > n
                    delim_one && first && print(io, delim)
                    break
                end
                first = false
                print(io, delim)
                print(io, ' ')
            end
        end
    end
    print(io, cl)
end

function richprint(io::IO, dom::Vector{Node}, t::Tuple)
    let id = startNode(dom, io, t, "Tuple")
        show_delim_array(io, dom, id, t, '(', ',', ')', true)
        endNode(dom, io, id)
    end
end
function richprint(io::IO, dom::Vector{Node}, v::SimpleVector)
    let id = startNode(dom, io, v, "SimpleVector")
        show_delim_array(io, dom, id, v, "svec(", ',', ')', false)
        endNode(dom, io, id)
    end
end

## Abstract Syntax Tree (AST) printing ##

const ExprNode = Base.ExprNode

richprint(io::IO, dom::Vector{Node}, ex::ExprNode) =
    let id = startNode(dom, io, ex, "ExprNode")
        richprint_unquoted_quote_expr(io, dom, ex, 0, -1)
        endNode(dom, io, id)
    end
richprint_unquoted(io::IO, dom::Vector{Node}, ex, ::Int=0, ::Int=0) = richprint(io, dom, ex)

## AST printing constants ##

import Base:
    indent_width, quoted_syms, uni_ops,
    expr_infix_wide, expr_infix, expr_infix_any,
    all_ops, expr_calls, expr_parens

# AST decoding helpers ##

import Base:
    is_id_start_char, is_id_char, isidentifier, isoperator,
    operator_precedence, prec_power, prec_decl,
    is_expr, is_linenumber, is_quoted, unquoted

## AST printing helpers ##

import Base: show_expr_type, emphasize, typeemphasize
#function show_expr_type(io::IO, ty::ANY, emph::Bool)
#    if ty === Function
#        print(io, "::F")
#    elseif ty === Core.IntrinsicFunction
#        print(io, "::I")
#    else
#        if emph && (!isleaftype(ty) || ty == Core.Box)
#            emphasize(io, "::$ty")
#        else
#            print(io, "::$ty")
#        end
#    end
#end

import Base: show_linenumber

# show a block, e g if/for/etc
function richprint_block(io::IO, dom::Vector{Node}, head, args::Vector, body, indent::Int)
    print(io, head, ' ')
    richprint_list(io, dom, args, ", ", indent)

    ind = head === :module || head === :baremodule ? indent : indent + indent_width
    exs = (is_expr(body, :block) || is_expr(body, :body)) ? body.args : Any[body]
    for ex in exs
        if !is_linenumber(ex)
            print(io, '\n', " "^ind)
        end
        richprint_unquoted(io, dom, ex, ind, -1)
    end
    print(io, '\n', " "^indent)
end
richprint_block(io::IO, dom::Vector{Node}, head, block, i::Int) = richprint_block(io, dom, head, [], block, i)
function richprint_block(io::IO, dom::Vector{Node}, head, arg, block, i::Int)
    if is_expr(arg, :block)
        richprint_block(io, dom, head, arg.args, block, i)
    else
        richprint_block(io, dom, head, Any[arg], block, i)
    end
end

# show an indented list
function richprint_list(io::IO, dom::Vector{Node}, items, sep, indent::Int, prec::Int=0, enclose_operators::Bool=false)
    n = length(items)
    if n == 0; return end
    indent += indent_width
    first = true
    for item in items
        !first && print(io, sep)
        parens = enclose_operators && isa(item,Symbol) && isoperator(item)
        parens && print(io, '(')
        richprint_unquoted(io, dom, item, indent, prec)
        parens && print(io, ')')
        first = false
    end
end
# show an indented list inside the parens (op, cl)
function richprint_enclosed_list(io::IO, dom::Vector{Node}, op, items, sep, cl, indent, prec=0, encl_ops=false)
    print(io, op)
    richprint_list(io, dom, items, sep, indent, prec, encl_ops)
    print(io, cl)
end

# show a normal (non-operator) function call, e.g. f(x, y) or A[z]
function richprint_call(io::IO, dom::Vector{Node}, head, func, func_args, indent)
    op, cl = expr_calls[head]
    if isa(func, Symbol) || (isa(func, Expr) &&
            (func.head == :. || func.head == :curly))
        richprint_unquoted(io, dom, func, indent)
    else
        print(io, '(')
        richprint_unquoted(io, dom, func, indent)
        print(io, ')')
    end
    if head == :(.)
        print(io, '.')
    end
    if !isempty(func_args) && isa(func_args[1], Expr) && func_args[1].head === :parameters
        print(io, op)
        richprint_list(io, dom, func_args[2:end], ", ", indent)
        print(io, "; ")
        richprint_list(io, dom, func_args[1].args, ", ", indent)
        print(io, cl)
    else
        richprint_enclosed_list(io, dom, op, func_args, ", ", cl, indent)
    end
end

### AST printing ##

richprint_unquoted(io::IO, dom::Vector{Node}, sym::Symbol, ::Int, ::Int) =
    let id = startNode(dom, io, sym, "Symbol")
        print(io, sym)
        endNode(dom, io, id)
    end
richprint_unquoted(io::IO, dom::Vector{Node}, ex::LineNumberNode, ::Int, ::Int) =
    let id = startNode(dom, io, ex, "LineNumberNode")
        show_linenumber(io, ex.line)
        endNode(dom, io, id)
    end
richprint_unquoted(io::IO, dom::Vector{Node}, ex::LabelNode, ::Int, ::Int) =
    let id = startNode(dom, io, ex, "LabelNode")
        print(io, ex.label, ": ")
        endNode(dom, io, id)
    end
richprint_unquoted(io::IO, dom::Vector{Node}, ex::GotoNode, ::Int, ::Int) =
    let id = startNode(dom, io, ex, "GotoNode")
        print(io, "goto ", ex.label)
        endNode(dom, io, id)
    end
richprint_unquoted(io::IO, dom::Vector{Node}, ex::GlobalRef, ::Int, ::Int) =
    let id = startNode(dom, io, ex, "GlobalRef")
        print(io, ex.mod, '.', ex.name)
        endNode(dom, io, id)
    end

function richprint_unquoted(io::IO, dom::Vector{Node}, ex::Slot, ::Int, ::Int)
    let id = startNode(dom, io, ex, "Slot")
        typ = isa(ex,TypedSlot) ? ex.typ : Any
        slotid = ex.id
        src = get(io, :SOURCEINFO, false)
        if isa(src, CodeInfo)
            slottypes = (src::CodeInfo).slottypes
            if isa(slottypes, Array) && slotid <= length(slottypes::Array)
                slottype = slottypes[slotid]
                # The Slot in assignment can somehow have an Any type
                if isa(slottype, Type) && isa(typ, Type) && slottype <: typ
                    typ = slottype
                end
            end
        end
        slotnames = get(io, :SOURCE_SLOTNAMES, false)
        if (isa(slotnames, Vector{String}) &&
            slotid <= length(slotnames::Vector{String}))
            print(io, (slotnames::Vector{String})[slotid])
        else
            print(io, "_", slotid)
        end
        emphstate = typeemphasize(io)
        if emphstate || (typ !== Any && isa(ex,TypedSlot))
            show_expr_type(io, typ, emphstate)
        end
        endNode(dom, io, id)
    end
end

function richprint_unquoted(io::IO, dom::Vector{Node}, ex::QuoteNode, indent::Int, prec::Int)
    let id = startNode(dom, io, ex, "QuoteNode")
        if isa(ex.value, Symbol)
            richprint_unquoted_quote_expr(io, dom, ex.value, indent, prec)
        else
            print(io, "\$(QuoteNode(")
            richprint(io, dom, ex.value)
            print(io, "))")
        end
        endNode(dom, io, id)
    end
end

function richprint_unquoted_quote_expr(io::IO, dom::Vector{Node}, value, indent::Int, prec::Int)
    if isa(value, Symbol) && !(value in quoted_syms)
        s = string(value)
        if isidentifier(s) || isoperator(value)
            print(io, ":")
            print(io, value)
        else
            print(io, "Symbol(", repr(s), ")")
        end
    else
        if isa(value, Expr) && value.head === :block
            richprint_block(io, dom, "quote", value, indent)
            print(io, "end")
        else
            print(io, ":(")
            richprint_unquoted(io, dom, value, indent+indent_width, -1)
            print(io, ")")
        end
    end
end

function richprint_generator(io::IO, dom::Vector{Node}, ex::Expr, indent::Int)
    if ex.head === :flatten
        fg = ex
        ranges = Any[]
        while isa(fg, Expr) && fg.head === :flatten
            push!(ranges, fg.args[1].args[2:end])
            fg = fg.args[1].args[1]
        end
        push!(ranges, fg.args[2:end])
        richprint_unquoted(io, dom, fg.args[1], indent)
        for r in ranges
            print(io, " for ")
            richprint_list(io, dom, r, ", ", indent)
        end
    else
        richprint_unquoted(io, dom, ex.args[1], indent)
        print(io, " for ")
        richprint_list(io, dom, ex.args[2:end], ", ", indent)
    end
end

function richprint_unquoted(io::IO, dom::Vector{Node}, ex::Expr, indent::Int, prec::Int)
    id = startNode(dom, io, ex, "Expr")
    head, args, nargs = ex.head, ex.args, length(ex.args)
    emphstate = typeemphasize(io)
    show_type = true
    if (ex.head == :(=) || ex.head == :line ||
        ex.head == :boundscheck ||
        ex.head == :gotoifnot ||
        ex.head == :return)
        show_type = false
    end
    if !emphstate && ex.typ === Any
        show_type = false
    end
    # dot (i.e. "x.y"), but not compact broadcast exps
    if head === :(.) && !is_expr(args[2], :tuple)
        richprint_unquoted(io, dom, args[1], indent + indent_width)
        print(io, '.')
        if is_quoted(args[2])
            richprint_unquoted(io, dom, unquoted(args[2]), indent + indent_width)
        else
            print(io, '(')
            richprint_unquoted(io, dom, args[2], indent + indent_width)
            print(io, ')')
        end

    # infix (i.e. "x <: y" or "x = y")
    elseif (head in expr_infix_any && nargs==2) || (head === :(:) && nargs==3)
        func_prec = operator_precedence(head)
        head_ = head in expr_infix_wide ? " $head " : head
        if func_prec <= prec
            richprint_enclosed_list(io, dom, '(', args, head_, ')', indent, func_prec, true)
        else
            richprint_list(io, dom, args, head_, indent, func_prec, true)
        end

    # list (i.e. "(1, 2, 3)" or "[1, 2, 3]")
    elseif haskey(expr_parens, head)               # :tuple/:vcat
        op, cl = expr_parens[head]
        if head === :vcat
            sep = "; "
        elseif head === :hcat || head === :row
            sep = " "
        else
            sep = ", "
        end
        head !== :row && print(io, op)
        richprint_list(io, dom, args, sep, indent)
        if nargs == 1
            if head === :tuple
                print(io, ',')
            elseif head === :vcat
                print(io, ';')
            end
        end
        head !== :row && print(io, cl)

    # function call
    elseif head === :call && nargs >= 1
        func = args[1]
        fname = isa(func,GlobalRef) ? func.name : func
        func_prec = operator_precedence(fname)
        if func_prec > 0 || fname in uni_ops
            func = fname
        end
        func_args = args[2:end]

        if (in(ex.args[1], (GlobalRef(Base, :bitcast), :throw)) ||
            ismodulecall(ex))
            show_type = false
        end
        if show_type
            prec = prec_decl
        end

        # scalar multiplication (i.e. "100x")
        if (func === :* &&
            length(func_args)==2 && isa(func_args[1], Real) && isa(func_args[2], Symbol))
            if func_prec <= prec
                richprint_enclosed_list(io, dom, '(', func_args, "", ')', indent, func_prec)
            else
                richprint_list(io, dom, func_args, "", indent, func_prec)
            end

        # unary operator (i.e. "!z")
        elseif isa(func,Symbol) && func in uni_ops && length(func_args) == 1
            richprint_unquoted(io, dom, func, indent)
            if isa(func_args[1], Expr) || func_args[1] in all_ops
                richprint_enclosed_list(io, dom, '(', func_args, ", ", ')', indent, func_prec)
            else
                richprint_unquoted(io, dom, func_args[1])
            end

        # binary operator (i.e. "x + y")
        elseif func_prec > 0 # is a binary operator
            na = length(func_args)
            if (na == 2 || (na > 2 && func in (:+, :++, :*))) &&
                    all(!isa(a, Expr) || a.head !== :... for a in func_args)
                sep = " $func "
                if func_prec <= prec
                    richprint_enclosed_list(io, dom, '(', func_args, sep, ')', indent, func_prec, true)
                else
                    richprint_list(io, dom, func_args, sep, indent, func_prec, true)
                end
            elseif na == 1
                # 1-argument call to normally-binary operator
                op, cl = expr_calls[head]
                print(io, "(")
                richprint_unquoted(io, dom, func, indent)
                print(io, ")")
                richprint_enclosed_list(io, dom, op, func_args, ", ", cl, indent)
            else
                richprint_call(io, dom, head, func, func_args, indent)
            end

        # normal function (i.e. "f(x,y)")
        else
            richprint_call(io, dom, head, func, func_args, indent)
        end

    # other call-like expressions ("A[1,2]", "T{X,Y}", "f.(X,Y)")
    elseif haskey(expr_calls, head) && nargs >= 1  # :ref/:curly/:calldecl/:(.)
        funcargslike = head == :(.) ? ex.args[2].args : ex.args[2:end]
        richprint_call(io, dom, head, ex.args[1], funcargslike, indent)

    # comprehensions
    elseif head === :typed_comprehension && length(args) == 2
        richprint_unquoted(io, dom, args[1], indent)
        print(io, '[')
        richprint_generator(io, dom, args[2], indent)
        print(io, ']')

    elseif head === :comprehension && length(args) == 1
        print(io, '[')
        richprint_generator(io, dom, args[1], indent)
        print(io, ']')

    elseif (head === :generator && length(args) >= 2) || (head === :flatten && length(args) == 1)
        print(io, '(')
        richprint_generator(io, dom, ex, indent)
        print(io, ')')

    elseif head === :filter && length(args) == 2
        richprint_unquoted(io, dom, args[2], indent)
        print(io, " if ")
        richprint_unquoted(io, dom, args[1], indent)

    # comparison (i.e. "x < y < z")
    elseif head === :comparison && nargs >= 3 && (nargs&1==1)
        comp_prec = minimum(operator_precedence, args[2:2:end])
        if comp_prec <= prec
            richprint_enclosed_list(io, dom, '(', args, " ", ')', indent, comp_prec)
        else
            richprint_list(io, dom, args, " ", indent, comp_prec)
        end

    # function calls need to transform the function from :call to :calldecl
    # so that operators are printed correctly
    elseif head === :function && nargs==2 && is_expr(args[1], :call)
        richprint_block(io, dom, head, Expr(:calldecl, args[1].args...), args[2], indent)
        print(io, "end")

    elseif head === :function && nargs == 1
        print(io, "function ", args[1], " end")

    # block with argument
    elseif head in (:for,:while,:function,:if) && nargs==2
        richprint_block(io, dom, head, args[1], args[2], indent)
        print(io, "end")

    elseif head === :module && nargs==3 && isa(args[1],Bool)
        richprint_block(io, dom, args[1] ? :module : :baremodule, args[2], args[3], indent)
        print(io, "end")

    # type declaration
    elseif head === :type && nargs==3
        richprint_block(io, dom, args[1] ? Symbol("mutable struct") : Symbol("struct"), args[2], args[3], indent)
        print(io, "end")

    elseif head === :bitstype && nargs == 2
        print(io, "primitive type ")
        richprint_list(io, dom, reverse(args), ' ', indent)
        print(io, " end")

    elseif head === :abstract && nargs == 1
        print(io, "abstract type ")
        richprint_list(io, dom, args, ' ', indent)
        print(io, " end")

    # empty return (i.e. "function f() return end")
    elseif head === :return && nargs == 1 && args[1] === nothing
        print(io, head)

    # type annotation (i.e. "::Int")
    elseif head === Symbol("::") && nargs == 1
        print(io, "::")
        richprint_unquoted(io, dom, args[1], indent)

    # var-arg declaration or expansion
    # (i.e. "function f(L...) end" or "f(B...)")
    elseif head === :(...) && nargs == 1
        richprint_unquoted(io, dom, args[1], indent)
        print(io, "...")

    elseif (nargs == 0 && head in (:break, :continue))
        print(io, head)

    elseif (nargs == 1 && head in (:return, :const)) ||
                          head in (:local,  :global, :export)
        print(io, head, ' ')
        richprint_list(io, dom, args, ", ", indent)

    elseif head === :macrocall && nargs >= 1
        # Use the functional syntax unless specifically designated with prec=-1
        if prec >= 0
            richprint_call(io, dom, :call, ex.args[1], ex.args[2:end], indent)
        else
            richprint_list(io, dom, args, ' ', indent)
        end

    elseif head === :line && 1 <= nargs <= 2
        show_linenumber(io, args...)

    elseif head === :if && nargs == 3     # if/else
        richprint_block(io, dom, "if",   args[1], args[2], indent)
        richprint_block(io, dom, "else", args[3], indent)
        print(io, "end")

    elseif head === :try && 3 <= nargs <= 4
        richprint_block(io, dom, "try", args[1], indent)
        if is_expr(args[3], :block)
            richprint_block(io, dom, "catch", args[2] === false ? Any[] : args[2], args[3], indent)
        end
        if nargs >= 4 && is_expr(args[4], :block)
            richprint_block(io, dom, "finally", Any[], args[4], indent)
        end
        print(io, "end")

    elseif head === :let && nargs >= 1
        richprint_block(io, dom, "let", args[2:end], args[1], indent); print(io, "end")

    elseif head === :block || head === :body
        richprint_block(io, dom, "begin", ex, indent); print(io, "end")

    elseif head === :quote && nargs == 1 && isa(args[1],Symbol)
        richprint_unquoted_quote_expr(io, dom, args[1], indent, 0)

    elseif head === :gotoifnot && nargs == 2
        print(io, "unless ")
        richprint_list(io, dom, args, " goto ", indent)

    elseif head === :string && nargs == 1 && isa(args[1], AbstractString)
        richprint(io, dom, args[1])

    elseif head === :null
        print(io, "nothing")

    elseif head === :kw && length(args)==2
        richprint_unquoted(io, dom, args[1], indent+indent_width)
        print(io, '=')
        richprint_unquoted(io, dom, args[2], indent+indent_width)

    elseif head === :string
        print(io, '"')
        for x in args
            if !isa(x,AbstractString)
                print(io, "\$(")
                if isa(x,Symbol) && !(x in quoted_syms)
                    print(io, x)
                else
                    richprint_unquoted(io, dom, x)
                end
                print(io, ")")
            else
                escape_string(io, x, "\"\$")
            end
        end
        print(io, '"')

    elseif (head === :&#= || head === :$=#) && length(args) == 1
        print(io, head)
        a1 = args[1]
        parens = (isa(a1,Expr) && a1.head !== :tuple) || (isa(a1,Symbol) && isoperator(a1))
        parens && print(io, "(")
        richprint_unquoted(io, dom, a1)
        parens && print(io, ")")

    # transpose
    elseif (head === Symbol('\'') || head === Symbol(".'")) && length(args) == 1
        if isa(args[1], Symbol)
            richprint_unquoted(io, dom, args[1])
        else
            print(io, "(")
            richprint_unquoted(io, dom, args[1])
            print(io, ")")
        end
        print(io, head)

    # `where` syntax
    elseif head === :where && length(args) > 1
        parens = 1 <= prec
        parens && print(io, "(")
        richprint_unquoted(io, dom, args[1], indent, operator_precedence(:(::)))
        print(io, " where ")
        if nargs == 2
            richprint_unquoted(io, dom, args[2], indent, 1)
        else
            print(io, "{")
            richprint_list(io, dom, args[2:end], ", ", indent)
            print(io, "}")
        end
        parens && print(io, ")")

    elseif head === :import || head === :importall || head === :using
        print(io, head)
        first = true
        for a = args
            if first
                print(io, ' ')
                first = false
            else
                print(io, '.')
            end
            if a !== :.
                print(io, a)
            end
        end
    elseif head === :meta && length(args) >= 2 && args[1] === :push_loc
        print(io, "# meta: location ", join(args[2:end], " "))
        show_type = false
    elseif head === :meta && length(args) == 1 && args[1] === :pop_loc
        print(io, "# meta: pop location")
        show_type = false
    # print anything else as "Expr(head, args...)"
    else
        show_type = false
        if emphstate && ex.head !== :lambda && ex.head !== :method
            io = IOContext(io, :TYPEEMPHASIZE => false)
            emphstate = false
        end
        print(io, "\$(Expr(")
        show(io, ex.head)
        for arg in args
            print(io, ", ")
            richprint(io, dom, arg)
        end
        print(io, "))")
    end
    show_type && show_expr_type(io, ex.typ, emphstate)
    endNode(dom, io, id)
end

import Base: unwrap_unionall
function richprint_tuple_as_call(io::IO, dom::Vector{Node}, id::Int, name::Symbol, sig::Type)
    # print a method signature tuple for a lambda definition
    if sig === Tuple
        Base.print_with_color(color, io, name, "(...)")
        return
    end
    sig = unwrap_unionall(sig).parameters
    ft = sig[1]
    uw = unwrap_unionall(ft)
    nextChild(dom, io, id)
    if ft <: Function && isa(uw,DataType) && isempty(uw.parameters) &&
            isdefined(uw.name.module, uw.name.mt.name) &&
            ft == typeof(getfield(uw.name.module, uw.name.mt.name))
        print(io, uw.name.mt.name)
    elseif isa(ft, DataType) && ft.name === Type.body.name && isleaftype(ft)
        f = ft.parameters[1]
        print(io, f)
    else
        print(io, "(::", ft, ")")
    end
    first = true
    print(io, "(")
    for i = 2:length(sig)
        nextChild(dom, io, id)
        first || print(io, ", ")
        first = false
        print(io, "::", sig[i])
    end
    print(io, ")")
    nothing
end

import Base: ismodulecall

function richprint(io::IO, dom::Vector{Node}, tv::TypeVar)
    # If we are in the `unionall_env`, the type-variable is bound
    # and the type constraints are already printed.
    # We don't need to print it again.
    # Otherwise, the lower bound should be printed if it is not `Bottom`
    # and the upper bound should be printed if it is not `Any`.
    let id = startNode(dom, io, tv, "TypeVar")
        in_env = (:unionall_env => tv) in io
        function show_bound(io::IO, b::ANY)
            parens = isa(b,UnionAll) && !print_without_params(b)
            parens && print(io, "(")
            richprint(io, dom, b)
            parens && print(io, ")")
        end
        lb, ub = tv.lb, tv.ub
        if !in_env && lb !== Union{}
            if ub === Any
                write(io, tv.name)
                print(io, ">:")
                show_bound(io, lb)
            else
                show_bound(io, lb)
                print(io, "<:")
                write(io, tv.name)
            end
        else
            write(io, tv.name)
        end
        if !in_env && ub !== Any
            print(io, "<:")
            show_bound(io, ub)
        end
        endNode(dom, io, id)
    end
end

import Base: alignment, undef_ref_str, undef_ref_alignment, replace_with_centered_mark, replace_in_print_matrix
import Base: dims2string, inds2string, summary, array_eltype_show_how

"""
`richprint_matrix_repr(io, X)` prints matrix X with opening and closing square brackets.
"""
function richprint_matrix_repr(io::IO, dom::Vector{Node}, id::Int, X::AbstractArray)
    compact, prefix = array_eltype_show_how(X)
    if compact && !haskey(io, :compact)
        io = IOContext(io, :compact => compact)
    end
    indr, indc = indices(X, 1), indices(X, 2)
    nr, nc = length(indr), length(indc)
    rr1, rr2 = UnitRange{Int}(indr), 1:0
    cr1, cr2 = UnitRange{Int}(indc), 1:0
    print(io, prefix, "[")
    for rr in (rr1, rr2)
        for i in rr
            for cr in (cr1, cr2)
                for j in cr
                    nextChild(dom, io, id)
                    j > first(cr) && print(io, " ")
                    if !isassigned(X, i, j)
                        print(io, undef_ref_str)
                    else
                        el = X[i, j]
                        richprint(io, dom, el)
                    end
                end
                if last(cr) == last(indc)
                    i < last(indr) && print(io, "; ")
                end
            end
        end
    end
    print(io, "]")
end

richprint(io::IO, dom::Vector{Node}, X::AbstractArray) = richprintarray(io, dom, X)
richprint(io::IO, dom::Vector{Node}, X::AbstractVector) = richprintvector(io, dom, X)

function richprintarray(io::IO, dom::Vector{Node}, X::AbstractArray)
    N = ndims(X)
    let id = startNode(dom, io, X, "Array-$N")
        if !haskey(io, :compact)
            io = IOContext(io, :compact => true)
        end
        if !isempty(X)
            if N == 0
                if isassigned(X)
                    return richprint(io, dom, X[])
                else
                    return print(io, undef_ref_str)
                end
            end
            if N <= 2
                richprint_matrix_repr(io, dom, id, X)
            else
                tailinds = tail(tail(indices(X)))
                nd = N - 2
                for I in CartesianRange(tailinds)
                    idxs = I.I
                    slice = view(X, indices(X, 1), indices(X, 2), idxs...)
                    richprint_matrix_repr(io, dom, id, slice)
                    print(io, idxs == map(last, tailinds) ? "" : "\n\n")
                end
            end
        else
            Base.repremptyarray(io, X)
        end
        endNode(dom, io, id)
    end
end

function richprintvector(io::IO, dom::Vector{Node}, v)
    let id = startNode(dom, io, v, "Vector")
        compact, prefix = array_eltype_show_how(v)
        if compact && !haskey(io, :compact)
            io = IOContext(io, :compact => compact)
        end
        print(io, prefix)
        show_delim_array(io, dom, id, v, '[', ',', ']', false)
        endNode(dom, io, id)
    end
end

max_stars_repo_name: vtjnash/Firelight


text:
using Dates
using MozillaCACerts_jll
using OpenSSL
using OpenSSL_jll
using Sockets
using Test

include("http_helpers.jl")

macro catch_exception_object(code)
    quote
        err = try
            $(esc(code))
            nothing
        catch e
            e
        end
        if err == nothing
            error("Expected exception, got $err.")
        end
        err
    end
end

# Verifies calling into OpenSSL library.
@testset "OpenSSL" begin
    @test OpenSSL.BIO_STREAM_CALLBACKS.x.on_bio_create_ptr != C_NULL
    @test OpenSSL.BIO_STREAM_CALLBACKS.x.on_bio_destroy_ptr != C_NULL
    @test OpenSSL.BIO_STREAM_CALLBACKS.x.on_bio_read_ptr != C_NULL
    @test OpenSSL.BIO_STREAM_CALLBACKS.x.on_bio_write_ptr != C_NULL
    @test OpenSSL.BIO_STREAM_CALLBACKS.x.on_bio_puts_ptr != C_NULL
    @test OpenSSL.BIO_STREAM_CALLBACKS.x.on_bio_ctrl_ptr != C_NULL
end

@testset "RandomBytes" begin
    random_data = random_bytes(64)

    @test length(random_data) == 64
end

@testset "BigNumbers" begin
    n1 = BigNum(0x4)
    n2 = BigNum(0x8)
    @test String(n1 + n2) == "0xC"

    #n4 = n3 - n1 - n2 - n3
    #n5 = BigNum(0x2)
    #@show n1, n2, n3, n4, n1 * n5

    n1 = BigNum(0x10)
    n2 = BigNum(0x4)
    @test String(n1 / n2) == "0x4"

    n1 = BigNum(0x11)
    @test String(n1 % n2) == "0x1"

    n1 = BigNum(0x3)
    @test String(n1 * n2) == "0xC"
end

@testset "Asn1Time" begin
    @test String(Asn1Time()) == "Jan  1 00:00:00 1970 GMT"
    @test String(Asn1Time(2)) == "Jan  1 00:00:02 1970 GMT"

    asn1_time = Asn1Time()
    OpenSSL.adjust(asn1_time, Dates.Second(4))

    OpenSSL.free(asn1_time)
    @test String(asn1_time) == "C_NULL"

    # double free
    OpenSSL.free(asn1_time)
    @test String(asn1_time) == "C_NULL"

    OpenSSL.adjust(asn1_time, Dates.Second(4))
    OpenSSL.adjust(asn1_time, Dates.Second(4))
    OpenSSL.adjust(asn1_time, Dates.Day(4))
    OpenSSL.adjust(asn1_time, Dates.Year(2))

    @show asn1_time
end

@testset "X509Name" begin
    x509_name_1 = X509Name()
    add_entry(x509_name_1, "C", "US")
    add_entry(x509_name_1, "ST", "Isles of Redmond")
    add_entry(x509_name_1, "CN", "www.redmond.com")

    x509_name_2 = X509Name()
    add_entry(x509_name_2, "C", "US")
    @test x509_name_1 != x509_name_2

    add_entry(x509_name_2, "ST", "Isles of Redmond")
    @test x509_name_1 != x509_name_2

    add_entry(x509_name_2, "CN", "www.redmond.com")
    @test x509_name_1 == x509_name_2
end

@testset "ReadPEMCert" begin
    file_handle = open(MozillaCACerts_jll.cacert)
    file_content = String(read(file_handle))
    close(file_handle)

    start_line = "==========\n"
    certs_pem = split(file_content, start_line; keepempty=false)

    cert = certs_pem[2]

    x509_cert = X509Certificate(cert)
    @test String(x509_cert.subject_name) == "/C=BE/O=GlobalSign nv-sa/OU=Root CA/CN=GlobalSign Root CA"
    @test String(x509_cert.issuer_name) == "/C=BE/O=GlobalSign nv-sa/OU=Root CA/CN=GlobalSign Root CA"
    @test String(x509_cert.time_not_before) == "Sep  1 12:00:00 1998 GMT"
    @test String(x509_cert.time_not_after) == "Jan 28 12:00:00 2028 GMT"

    # finalizer will cleanup
    #finalize(x509_cert)
end

@testset "StackOf{X509Certificate}" begin
    file_handle = open(MozillaCACerts_jll.cacert)
    file_content = String(read(file_handle))
    close(file_handle)

    start_line = "==========\n"
    certs_pem = split(file_content, start_line; keepempty=false)

    x509_certificates = StackOf{X509Certificate}()

    foreach(2:length(certs_pem)) do i
        x509_cert = X509Certificate(certs_pem[i])
        push!(x509_certificates, x509_cert)
        finalize(x509_cert)
        nothing
    end

    free(x509_certificates)
end

@testset "StackOf{BigNum}" begin
    n1 = BigNum(0x4)
    n2 = BigNum(0x8)

    big_nums = StackOf{BigNum}()
    push!(big_nums, n1)
    push!(big_nums, n2)

    _n1 = pop!(big_nums)
    _n2 = pop!(big_nums)

    @test _n1 == n2
    @test _n1 == n2
end

@testset "X509Store" begin
    file_handle = open(MozillaCACerts_jll.cacert)
    file_content = String(read(file_handle))
    close(file_handle)

    start_line = "==========\n"

    certs_pem = split(file_content, start_line; keepempty=false)

    # X509 store.
    x509_store = X509Store()

    foreach(2:length(certs_pem)) do i
        x509_cert = X509Certificate(certs_pem[i])
        add_cert(x509_store, x509_cert)
        free(x509_cert)
        nothing
    end

    free(x509_store)
end

@testset "HttpsConnect" begin
    tcp_stream = connect("www.nghttp2.org", 443)

    ssl_ctx = OpenSSL.SSLContext(OpenSSL.TLSv12ClientMethod())
    result = OpenSSL.ssl_set_options(ssl_ctx, OpenSSL.SSL_OP_NO_COMPRESSION)

    # Create SSL stream.
    ssl_stream = SSLStream(ssl_ctx, tcp_stream, tcp_stream)

    #TODO expose connect
    OpenSSL.connect(ssl_stream)

    x509_server_cert = OpenSSL.get_peer_certificate(ssl_stream)

    @test String(x509_server_cert.issuer_name) == "/C=US/O=Let's Encrypt/CN=R3"
    @test String(x509_server_cert.subject_name) == "/CN=nghttp2.org"

    request_str = "GET / HTTP/1.1\r\nHost: www.nghttp2.org\r\nUser-Agent: curl\r\nAccept: */*\r\n\r\n"

    written = unsafe_write(ssl_stream, pointer(request_str), length(request_str))

    response = read(ssl_stream)

    close(ssl_stream)
    finalize(ssl_ctx)
end

@testset "ClosedStream" begin
    tcp_stream = connect("www.nghttp2.org", 443)

    ssl_ctx = OpenSSL.SSLContext(OpenSSL.TLSv12ClientMethod())
    result = OpenSSL.ssl_set_options(ssl_ctx, OpenSSL.SSL_OP_NO_COMPRESSION)
    OpenSSL.ssl_set_ciphersuites(ssl_ctx, "TLS_AES_256_GCM_SHA384:TLS_CHACHA20_POLY1305_SHA256:TLS_AES_128_GCM_SHA256")

    ssl_stream = SSLStream(ssl_ctx, tcp_stream, tcp_stream)

    OpenSSL.connect(ssl_stream)

    # Close the ssl stream.
    close(ssl_stream)

    request_str = "GET / HTTP/1.1\r\nHost: www.nghttp2.org\r\nUser-Agent: curl\r\nAccept: */*\r\n\r\n"

    err = @catch_exception_object unsafe_write(ssl_stream, pointer(request_str), length(request_str))
    @test typeof(err) == Base.IOError

    finalize(ssl_ctx)
end

@testset "NoCloseStream" begin
    ssl_ctx = OpenSSL.SSLContext(OpenSSL.TLSv12ClientMethod())
    result = OpenSSL.ssl_set_options(ssl_ctx, OpenSSL.SSL_OP_NO_COMPRESSION)

    # Create SSL stream.
    tcp_stream = connect("www.nghttp2.org", 443)
    ssl_stream = SSLStream(ssl_ctx, tcp_stream, tcp_stream)
    OpenSSL.connect(ssl_stream)

    request_str = "GET / HTTP/1.1\r\nHost: www.nghttp2.org\r\nUser-Agent: curl\r\nAccept: */*\r\n\r\n"
    unsafe_write(ssl_stream, pointer(request_str), length(request_str))

    response = read(ssl_stream)
    @test contains(String(response), "HTTP/1.1 200 OK")

    # Do not close SSLStream, leave it to the finalizer.
    #close(ssl_stream)
    #finalize(ssl_ctx)
end

@testset "InvalidStream" begin
    ssl_ctx = OpenSSL.SSLContext(OpenSSL.TLSv12ClientMethod())
    result = OpenSSL.ssl_set_options(ssl_ctx, OpenSSL.SSL_OP_NO_COMPRESSION)

    # Create SSL stream.
    tcp_stream = connect("www.nghttp2.org", 443)
    ssl_stream = SSLStream(ssl_ctx, tcp_stream, tcp_stream)

    err = @catch_exception_object read(ssl_stream)
    @test typeof(err) == OpenSSL.OpenSSLError

    close(ssl_stream)
    free(ssl_ctx)
end

@testset "Hash" begin
    res = digest(EvpMD5(), IOBuffer("The quick brown fox jumps over the lazy dog"))
    @test res == UInt8[0x9e, 0x10, 0x7d, 0x9d, 0x37, 0x2b, 0xb6, 0x82, 0x6b, 0xd8, 0x1d, 0x35, 0x42, 0xa4, 0x19, 0xd6]
end

@testset "SelfSignedCertificate" begin
    x509_certificate = X509Certificate()

    evp_pkey = EvpPKey(rsa_generate_key())
    x509_certificate.public_key = evp_pkey

    x509_name = X509Name()
    add_entry(x509_name, "C", "US")
    add_entry(x509_name, "ST", "Isles of Redmond")
    add_entry(x509_name, "CN", "www.redmond.com")

    x509_certificate.subject_name = x509_name
    x509_certificate.issuer_name = x509_name

    adjust(x509_certificate.time_not_before, Second(0))
    adjust(x509_certificate.time_not_after, Year(1))

    add_extension(x509_certificate, X509Extension("basicConstraints", "CA:TRUE"))
    add_extension(x509_certificate, X509Extension("keyUsage", "keyCertSign"))

    sign_certificate(x509_certificate, evp_pkey)

    iob = IOBuffer()
    write(iob, x509_certificate)

    seek(iob, 0)
    cert_pem = String(read(iob))

    x509_certificate2 = X509Certificate(cert_pem)

    x509_string = String(x509_certificate)
    x509_string2 = String(x509_certificate2)

    public_key = x509_certificate.public_key

    @test x509_string == x509_string2

    p12_object = P12Object(evp_pkey, x509_certificate)

    OpenSSL.unpack(p12_object)

    #seek(io, 0)
    # @show "p12:" String(read(io))
end

@testset "SignCertCertificate" begin
    # Create a root certificate.
    x509_certificate = X509Certificate()

    evp_pkey_ca = EvpPKey(rsa_generate_key())
    x509_certificate.public_key = evp_pkey_ca

    x509_name = X509Name()
    add_entry(x509_name, "C", "US")
    add_entry(x509_name, "ST", "Isles of Redmond")
    add_entry(x509_name, "CN", "www.redmond.com")

    x509_certificate.subject_name = x509_name
    x509_certificate.issuer_name = x509_name

    adjust(x509_certificate.time_not_before, Second(0))
    adjust(x509_certificate.time_not_after, Year(1))

    add_extension(x509_certificate, X509Extension("basicConstraints", "CA:TRUE"))
    add_extension(x509_certificate, X509Extension("keyUsage", "keyCertSign"))

    sign_certificate(x509_certificate, evp_pkey_ca)

    root_certificate = x509_certificate

    # Create a certificate sign request.
    x509_request = X509Request()

    evp_pkey = EvpPKey(rsa_generate_key())

    x509_name = X509Name()
    add_entry(x509_name, "C", "US")
    add_entry(x509_name, "ST", "Isles of Redmond")
    add_entry(x509_name, "CN", "www.redmond.com")

    x509_request.subject_name = x509_name

    x509_exts = StackOf{X509Extension}()

    ext = X509Extension("subjectAltName", "DNS:localhost")
    push!(x509_exts, ext)
    add_extensions(x509_request, x509_exts)
    finalize(ext)

    finalize(x509_exts)

    sign_request(x509_request, evp_pkey)

    # Create a certificate.
    x509_certificate = X509Certificate()
    x509_certificate.version = 2

    # Set issuer and subject name of the cert from the req and CA.
    x509_certificate.subject_name = x509_request.subject_name
    x509_certificate.issuer_name = root_certificate.subject_name

    x509_exts = x509_request.extensions

    ext = pop!(x509_exts)

    add_extension(x509_certificate, ext)
    add_extension(x509_certificate, X509Extension("keyUsage", "digitalSignature, nonRepudiation, keyEncipherment"))
    add_extension(x509_certificate, X509Extension("basicConstraints", "CA:FALSE"))

    # Set public key
    x509_certificate.public_key = x509_request.public_key

    adjust(x509_certificate.time_not_before, Second(0))
    adjust(x509_certificate.time_not_after, Year(1))

    sign_certificate(x509_certificate, evp_pkey_ca)
end

@testset "ErrorTaskTLS" begin
    err_msg = OpenSSL.get_error()
    @test err_msg == ""

    ssl_ctx = OpenSSL.SSLContext(OpenSSL.TLSv12ServerMethod())

    # Make direct invalid call to OpenSSL
    invalid_cipher_suites = "TLS_AES_356_GCM_SHA384:TLS_CHACHA20_POLY1305_SHA256:TLS_AES_128_GCM_SHA256"
    result = ccall(
        (:SSL_CTX_set_ciphersuites, libssl),
        Cint,
        (OpenSSL.SSLContext, Cstring),
        ssl_ctx,
        invalid_cipher_suites)

    # Verify the error message.
    err_msg = OpenSSL.get_error()
    @test contains(err_msg, "no cipher match")

    # Ensure error queue is empty.
    err_msg = OpenSSL.get_error()
    @test err_msg == ""

    # Make invalid OpenSSL (with fail and OpenSSL updates internal error queue).
    result = ccall(
        (:SSL_CTX_set_ciphersuites, libssl),
        Cint,
        (OpenSSL.SSLContext, Cstring),
        ssl_ctx,
        invalid_cipher_suites)
    # Copy and clear OpenSSL error queue to task TLS.
    OpenSSL.update_tls_error_state()
    # OpenSSL queue should be empty right now.
    @test ccall((:ERR_peek_error, libcrypto), Culong, ()) == 0

    # Verify the error message, error message should be retrived from the task TLS.
    err_msg = OpenSSL.get_error()
    @test contains(err_msg, "no cipher match")

    free(ssl_ctx)
end

@testset "PKCS12" begin
    x509_certificate = X509Certificate()

    evp_pkey = EvpPKey(rsa_generate_key())
    x509_certificate.public_key = evp_pkey

    x509_name = X509Name()
    add_entry(x509_name, "C", "US")
    add_entry(x509_name, "ST", "Isles of Redmond")
    add_entry(x509_name, "CN", "www.redmond.com")

    x509_certificate.subject_name = x509_name
    x509_certificate.issuer_name = x509_name

    adjust(x509_certificate.time_not_before, Second(0))
    adjust(x509_certificate.time_not_after, Year(1))

    sign_certificate(x509_certificate, evp_pkey)

    p12_object = P12Object(evp_pkey, x509_certificate)

    _evp_pkey, _x509_certificate, _x509_ca_stack = unpack(p12_object)

    @test _evp_pkey == evp_pkey
    @test _evp_pkey.key_type == evp_pkey.key_type

    @test _x509_certificate == x509_certificate
    @test _x509_certificate.subject_name == x509_certificate.subject_name
    @test _x509_certificate.issuer_name == x509_certificate.issuer_name
end

@testset "Encrypt" begin
    evp_ciphers = [
        EvpEncNull(),
        EvpBlowFishCBC(),
        EvpBlowFishECB(),
        #EvpBlowFishCFB(),
        EvpBlowFishOFB(),
        EvpAES128CBC(),
        EvpAES128ECB(),
        #EvpAES128CFB(),
        EvpAES128OFB(),
    ]

    foreach(evp_ciphers) do evp_cipher
        sym_key = random_bytes(evp_cipher.key_length)
        init_vector = random_bytes(evp_cipher.init_vector_length)

        enc_evp_cipher_ctx = EvpCipherContext()
        encrypt_init(enc_evp_cipher_ctx, evp_cipher, sym_key, init_vector)

        dec_evp_cipher_ctx = EvpCipherContext()
        decrypt_init(dec_evp_cipher_ctx, evp_cipher, sym_key, init_vector)

        in_string = "OpenSSL Julia"
        in_data = IOBuffer(in_string)
        enc_data = IOBuffer()

        cipher(enc_evp_cipher_ctx, in_data, enc_data)
        seek(enc_data, 0)
        @show String(read(enc_data))
        seek(enc_data, 0)

        dec_data = IOBuffer()
        cipher(dec_evp_cipher_ctx, enc_data, dec_data)
        out_data = take!(dec_data)
        out_string = String(out_data)

        @test in_string == out_string
    end
end

@testset "StackOf{X509Extension}" begin
    ext1 = X509Extension("subjectAltName", "DNS:openssl.jl.com")
    ext2 = X509Extension("keyUsage", "digitalSignature, keyEncipherment, keyAgreement")
    ext3 = X509Extension("basicConstraints", "CA:FALSE")

    st = StackOf{X509Extension}()
    push!(st, ext1)
    push!(st, ext2)
    push!(st, ext3)

    @test String(ext1) == "DNS:openssl.jl.com"
    @test String(ext2) == "Digital Signature, Key Encipherment, Key Agreement"
    @test String(ext3) == "CA:FALSE"

    finalize(ext1)
    finalize(ext2)
    finalize(ext3)

    @test length(st) == 3

    ext_1 = pop!(st)
    ext_2 = pop!(st)
    ext_3 = pop!(st)

    @test length(st) == 0

    @test String(ext_1) == "CA:FALSE"
    @test String(ext_2) == "Digital Signature, Key Encipherment, Key Agreement"
    @test String(ext_3) == "DNS:openssl.jl.com"

    finalize(ext_1)
    finalize(ext_2)
    finalize(ext_3)

    finalize(st)
end

@testset "SerializePrivateKey" begin
    evp_pkey = EvpPKey(rsa_generate_key())

    io = IOBuffer()
    write(io, evp_pkey)

    pkey_pem = String(take!(io))

    @test startswith(pkey_pem, "-----BEGIN PRIVATE KEY-----")

    _evp_pkey = EvpPKey(pkey_pem)

    @test _evp_pkey == evp_pkey

    free(evp_pkey)
    free(_evp_pkey)
end

@testset "DSA" begin
    dsa = dsa_generate_key()
end

@testset "X509Attribute" begin
    attr = X509Attribute()
    free(attr)
end

@testset "SSLServer" begin
    f1 = @async test_server()
    f2 = @async test_client()

    fetch(f2)
    fetch(f1)
end

max_stars_repo_name: grlap/OpenSSL

max_issues_repo_name: grlap/OpenSSL

max_forks_repo_name: grlap/OpenSSL


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm

n = 3000 #1000 #20000
t = 0.15 #0.4 #0.15

function randomlines(n=1000,t=0.2)
	V = zeros(Float64,2,2*n)
	EV = [zeros(Int64,2) for k=1:n]
	for k=1:n
		v1 = rand(Float64,2)
		v2 = rand(Float64,2)
		vm = (v1+v2)/2
		transl = rand(Float64,2)
		V[:,k] = (v1-vm)*t + transl
		V[:,n+k] = (v2-vm)*t + transl
		EV[k] = [k,n+k]
	end
	V = Plasm.normalize(V)
	model = (V,EV)
	return model
end

model = randomlines(n,t)
Sigma = Lar.spaceindex(model)
Plasm.view(model)

W,EW = Lar.fragmentlines(model)
Plasm.viewexploded(W,EW)(1.2,1.2,1.2)

V,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells)) # 2-connected components (H & T)
hpcs = [ Plasm.lar2hpc(V,EVs[i]) for i=1:length(EVs) ]
Plasm.view([ Plasm.color(Plasm.colorkey[(k%12)==0 ? 12 : k%12])(hpcs[k]) for k=1:(length(hpcs)) ])


text:
#2-dim CAs

#Two dimensional Cellular Automaton
#TODO: Change naming to be consistent
mutable struct CA2d

    #User given values
    k::Int #Number of states
    r::Int #r-nearest neigbors

    #Internal values
    cells::Array{Int8, 3}

    function CA2d(B::Array{Int,1},
                  S::Array{Int,1},
                  init::Array{Int,2},
                  gen::Int,
                  k::Int=2,
                  r::Int=1)

        h, w = size(init)

        cells = Array{Int8}(undef, (h, w, gen)) #Syntax A(T, dims) is deprecated
        cells[:, :, 1] = Array{Int8}(init[:, :])

        for g = 2:gen
            for i = 1:h, j = 1:w
                cc = -cells[i, j, g-1]
                for p = (i-r):(i+r), q = (j-r):(j+r)

                    #Cyclic boundary conditions
                    if p < 1; p = h-p; end
                    if p > h; p = p-h; end
                    if q < 1; q = w-q; end
                    if q > w; q = q-w; end

                    cc += cells[p, q, g-1]
                end
                cells[i, j, g] = eval_rule(cc, cells[i, j, g-1], B, S)
            end #hw ij
        end #gen

        new(k, r, cells)
    end
end

#Evaluate life rules
#TODO: make this more general
function eval_rule(cc, olds, B, S, k=2)
    if cc in B || (olds == 1 && cc in S)
        return 1
    end

    return 0
end


text:
import Ipopt, NLopt
using Nonconvex, LinearAlgebra, Test

f(x::AbstractVector) = x[2] * sin(10x[2])
g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2]

@testset "Deflated Ipopt" begin
    @testset "First order - $first_order" for first_order in [true, false]
        options = DeflatedOptions(
            ndeflations = 1, sub_options = IpoptOptions(),
        )
        m = Model(f)
        addvar!(m, [0.0, 0.0], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = DeflatedAlg(IpoptAlg())
        x0 = [1.234, 2.345]
        r = Nonconvex.optimize(m, alg, x0, options = options)
        @test abs(r.solutions[1][1][2] - 2.9879) <= 1e-4
        @test abs(r.solutions[2][1][2] - 4.2435) <= 1e-4
    end
end

@testset "Deflated NLopt MMA" begin
    options = DeflatedOptions(
        ndeflations = 1, sub_options = NLoptOptions(),
    )
    m = Model(f)
    addvar!(m, [0.0, 0.0], [10.0, 10.0])
    add_ineq_constraint!(m, x -> g(x, 2, 0))
    add_ineq_constraint!(m, x -> g(x, -1, 1))

    alg = DeflatedAlg(NLoptAlg(:LD_MMA))
    x0 = [1.234, 2.345]
    r = Nonconvex.optimize(m, alg, x0, options = options)
    @test abs(r.solutions[1][1][2] - 5.4996) <= 1e-4
    @test abs(r.solutions[2][1][2] - 6.7559) <= 1e-4
end


text:
alphabet = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"

function break_caesar_cipher(ciphertext, known_word)
    plaintext = ''

    # Your code goes here!

    return plaintext
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
using Random

"Where to download gutenberg from"
gutenbergurl = "http://www.gutenberg.org/files"

"Where to download gutenberg to"
gutenbergdir = joinpath(@__DIR__, "gutenberg")

"Download text from Project Gutenberg and return contents as String."
function gutenberg(name, fdir=nothing)
    (fdir == nothing) && (fdir = name)
    fname = split(fdir, "/")[end]
    isdir(gutenbergdir) || mkpath(gutenbergdir)
    path = joinpath(gutenbergdir, "$fname.txt")
    if !isfile(path)
        @info("Downloading Gutenberg $fname")
        url = "$gutenbergurl/$name/$fdir.txt"
        download(url,path)
    end
    return read(path, String)
end

"Filter and split Shakespeare text."
function shakespeare()
    global _shakespeare_trn, _shakespeare_tst, _shakespeare_chars
    if !(@isdefined _shakespeare_trn)
        s = gutenberg(100, "old/shaks12")
        a = []
        pos1 = 1
        while true
            pos2 = first(something(findnext("<<THIS", s, pos1), 0:-1)) - 1
            if pos2 == -1
                push!(a, s[pos1:end])
                break
            end
            push!(a, s[pos1:pos2])
            pos1 = last(something(findnext("SHIP.>>", s, pos2), 0:-1)) + 1
        end
        # 218 shakespeare texts
        a = a[3:end-1]
        # split trn, tst
        a = shuffle(MersenneTwister(42),a)
        trn = string(a[1:200]...)
        tst = string(a[201:end]...)
        # construct char vocab
        h = Dict{Char,Int}()
        for txt in (trn,tst), c in txt
            h[c] = 1 + get(h,c,0)
        end
        _shakespeare_chars = sort(collect(keys(h)),by=(x->h[x]),rev=true)
        _shakespeare_trn = UInt8[ something(findfirst(isequal(c), _shakespeare_chars), 0) for c in trn];
        _shakespeare_tst = UInt8[ something(findfirst(isequal(c), _shakespeare_chars), 0) for c in tst];
    end
    return _shakespeare_trn, _shakespeare_tst, _shakespeare_chars
end

nothing


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-beta2
#     language: julia
#     name: julia-1.7
# ---

# %%
VERSION

# %%
f(x, y) = (1 - x)^2 + 100(y - x^2)^2
X = range(-5, 5, 501)[1:end-1]
Y = range(0, 5, 501)

# %%
findmin(sin, Iterators.reverse(collect(X)))

# %%
findmin(sin∘last, enumerate(X))

# %%
findmin(sin∘first, zip(X, Y))

# %%
findmin(sin, Iterators.accumulate(+, X))

# %%
findmin(sin, Iterators.take(X, 486))

# %%
findmin(v -> f(v...), Iterators.product(X, Y))

# %%
findmin(sin, Iterators.flatten((X, Y)))

# %%
Base.keys(rv::Iterators.Reverse) = keys(reverse(rv.itr))
Base.keys(en::Iterators.Enumerate) = keys(en.itr)
Base.keys(::Iterators.Zip) = Iterators.countfrom()
Base.keys(::Iterators.Filter) = Iterators.countfrom()
Base.keys(ac::Iterators.Accumulate) = keys(ac.itr)
Base.keys(::Iterators.Rest) = Iterators.countfrom()
Base.keys(::Iterators.Count) = Iterators.countfrom()
Base.keys(::Iterators.Take) = Iterators.countfrom()
Base.keys(::Iterators.Drop) = Iterators.countfrom()
Base.keys(::Iterators.TakeWhile) = Iterators.countfrom()
Base.keys(::Iterators.DropWhile) = Iterators.countfrom()
Base.keys(::Iterators.Cycle) = Iterators.countfrom()
Base.keys(::Iterators.Repeated) = Iterators.countfrom()
Base.keys(pr::Iterators.ProductIterator) = CartesianIndices(axes(pr))
Base.keys(::Iterators.Flatten) = Iterators.countfrom()
Base.keys(::Iterators.PartitionIterator) = Iterators.countfrom()

# %%
findmin(sin, Iterators.reverse(collect(X)))

# %%
findmin(sin∘last, enumerate(X))

# %%
findmin(sin∘first, zip(X, Y))

# %%
findmin(sin, Iterators.accumulate(+, X))

# %%
findmin(sin, Iterators.take(X, 486))

# %%
findmin(v -> f(v...), Iterators.product(X, Y))

# %%
findmin(sin, Iterators.flatten((X, Y)))

# %%

max_stars_repo_path: 0008/definition of keys methods of iterators.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
struct SAFTVRMieParam <: EoSParam
    segment::SingleParam{Float64}
    sigma::PairParam{Float64}
    lambda_a::PairParam{Float64}
    lambda_r::PairParam{Float64}
    epsilon::PairParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
    Mw::SingleParam{Float64}
end

abstract type SAFTVRMieModel <: SAFTModel end
@newmodel SAFTVRMie SAFTVRMieModel SAFTVRMieParam

export SAFTVRMie
function SAFTVRMie(components; idealmodel=BasicIdeal, userlocations=String[], ideal_userlocations=String[], verbose=false)
    params,sites = getparams(components, ["SAFT/SAFTVRMie", "properties/molarmass.csv"]; userlocations=userlocations, verbose=verbose)

    params["Mw"].values .*= 1E-3
    Mw = params["Mw"]
    segment = params["m"]
    params["sigma"].values .*= 1E-10
    sigma = sigma_LorentzBerthelot(params["sigma"])
    epsilon = epsilon_HudsenMcCoubrey(params["epsilon"], sigma)
    lambda_a = lambda_LorentzBerthelot(params["lambda_a"])
    lambda_r = lambda_LorentzBerthelot(params["lambda_r"])
    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bondvol"]

    packagedparams = SAFTVRMieParam(segment, sigma, lambda_a, lambda_r, epsilon, epsilon_assoc, bondvol, Mw)
    references = ["10.1063/1.4819786", "10.1080/00268976.2015.1029027"]

    model = SAFTVRMie(packagedparams, sites, idealmodel; ideal_userlocations=ideal_userlocations, references=references, verbose=verbose)
    return model
end

function a_res(model::SAFTVRMieModel, V, T, z)
    return @f(a_mono) + @f(a_chain) + @f(a_assoc)
end

function a_mono(model::SAFTVRMieModel, V, T, z)
    return @f(a_hs)+@f(a_disp)
end
function a_disp(model::SAFTVRMieModel, V, T, z)
    return @f(a_1)+@f(a_2)+@f(a_3)
end

function a_hs(model::SAFTVRMieModel, V, T, z)
    ζ0 = @f(ζn, 0)
    ζ1 = @f(ζn, 1)
    ζ2 = @f(ζn, 2)
    ζ3 = @f(ζn, 3)
    N = N_A*∑(z)
    return 6*V/π/N*(3ζ1*ζ2/(1-ζ3) + ζ2^3/(ζ3*(1-ζ3)^2) + (ζ2^3/ζ3^2-ζ0)*log(1-ζ3))
end

function ζn(model::SAFTVRMieModel, V, T, z, n)
    return π/6*@f(ρ_S) * ∑(@f(x_S,i)*@f(d,i)^n for i ∈ @comps)
end

function ρ_S(model::SAFTVRMieModel, V, T, z)
    ∑z = ∑(z)
    N = N_A*∑z
    x = z/∑z
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    return N/V*m̄
end

function x_S(model::SAFTVRMieModel, V, T, z, i)
    x = z/∑(z)
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    return x[i]*m[i]/m̄
end

function d(model::SAFTVRMieModel, V, T, z, i)
    ϵ = model.params.epsilon.diagvalues
    σ = model.params.sigma.diagvalues
    λr = model.params.lambda_r.diagvalues
    λa = model.params.lambda_a.diagvalues
    u = SAFTVRMieconsts.u
    w = SAFTVRMieconsts.w
    θ = @f(C,i,i)*ϵ[i]/T
    return σ[i]*(1-∑(w[j]*(θ/(θ+u[j]))^(1/λr[i])*(exp(θ*(1/(θ/(θ+u[j]))^(λa[i]/λr[i])-1))/(u[j]+θ)/λr[i]) for j ∈ 1:5))
end

function d(model::SAFTVRMieModel, V, T, z, i, j)
    return (@f(d,i)+@f(d,j))/2
end

function C(model::SAFTVRMieModel, V, T, z, i, j)
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    return (λr[i,j]/(λr[i,j]-λa[i,j]))*(λr[i,j]/λa[i,j])^(λa[i,j]/(λr[i,j]-λa[i,j]))
end

function ζ_X(model::SAFTVRMieModel, V, T, z)
    comps = @comps
    return π/6*@f(ρ_S)*∑(@f(x_S,i)*@f(x_S,j)*(@f(d,i)+@f(d,j))^3/8 for i ∈ comps for j ∈ comps)
end

function a_1(model::SAFTVRMieModel, V, T, z)
    comps = @comps
    x = z/∑(z)
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    return m̄/T*∑(@f(x_S,i)*@f(x_S,j)*@f(a_1,i,j) for i ∈ comps for j ∈ comps)
end

function a_1(model::SAFTVRMieModel, V, T, z, i, j)
    ϵ = model.params.epsilon.values
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    x_0ij = @f(x_0,i,j)
    return 2*π*ϵ[i,j]*@f(d,i,j)^3*@f(C,i,j)*@f(ρ_S) *
        ( x_0ij^λa[i,j]*(@f(aS_1,λa[i,j])+@f(B,λa[i,j],x_0ij)) -
            x_0ij^λr[i,j]*(@f(aS_1,λr[i,j])+@f(B,λr[i,j],x_0ij)) )
end

function aS_1(model::SAFTVRMieModel, V, T, z, λ)
    ζeff_ = @f(ζeff,λ)
    return -1/(λ-3)*(1-ζeff_/2)/(1-ζeff_)^3
end

function ζeff(model::SAFTVRMieModel, V, T, z, λ)
    A = SAFTγMieconsts.A
    ζ_X_ = @f(ζ_X)
    return A * [1; 1/λ; 1/λ^2; 1/λ^3] ⋅ [ζ_X_; ζ_X_^2; ζ_X_^3; ζ_X_^4]
end

function B(model::SAFTVRMieModel, V, T, z, λ, x_0)
    I = (1-x_0^(3-λ))/(λ-3)
    J = (1-(λ-3)*x_0^(4-λ)+(λ-4)*x_0^(3-λ))/((λ-3)*(λ-4))
    ζ_X_ = @f(ζ_X)
    return I*(1-ζ_X_/2)/(1-ζ_X_)^3-9*J*ζ_X_*(ζ_X_+1)/(2*(1-ζ_X_)^3)
end

function x_0(model::SAFTVRMieModel, V, T, z, i, j)
    σ = model.params.sigma.values
    return σ[i,j]/@f(d,i,j)
end

function a_2(model::SAFTVRMieModel, V, T, z)
    comps = @comps
    x = z/∑(z)
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    return m̄/T^2*∑(@f(x_S,i)*@f(x_S,j)*@f(a_2,i,j) for i ∈ comps for j ∈ comps)
end

function a_2(model::SAFTVRMieModel, V, T, z,i,j)
    ϵ = model.params.epsilon.values
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    x_0ij = @f(x_0,i,j)
    ζ_X_ = @f(ζ_X)
    return π*@f(KHS)*(1+@f(χ,i,j))*@f(ρ_S)*ϵ[i,j]^2*@f(d,i,j)^3*@f(C,i,j)^2 *
        (x_0ij^(2*λa[i,j])*(@f(aS_1,2*λa[i,j])+@f(B,2*λa[i,j],x_0ij))
        - 2*x_0ij^(λa[i,j]+λr[i,j])*(@f(aS_1,λa[i,j]+λr[i,j])+@f(B,λa[i,j]+λr[i,j],x_0ij))
        + x_0ij^(2*λr[i,j])*(@f(aS_1,2*λr[i,j])+@f(B,2*λr[i,j],x_0ij)))
end

function KHS(model::SAFTVRMieModel, V, T, z)
    ζ_X_ = @f(ζ_X)
    return (1-ζ_X_)^4/(1+4ζ_X_+4ζ_X_^2-4ζ_X_^3+ζ_X_^4)
end

function χ(model::SAFTVRMieModel, V, T, z,i,j)
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    ζst_ = @f(ζst)
    α = @f(C,i,j)*(1/(λa[i,j]-3)-1/(λr[i,j]-3))
    return @f(f,α,1)*ζst_+@f(f,α,2)*ζst_^5+@f(f,α,3)*ζst_^8
end

function f(model::SAFTVRMieModel, V, T, z, α, m)
    ϕ = SAFTVRMieconsts.ϕ
    return ∑(ϕ[i+1][m]*α^i for i ∈ 0:3)/(1+∑(ϕ[i+1][m]*α^(i-3) for i ∈ 4:6))
end

function ζst(model::SAFTVRMieModel, V, T, z)
    comps = @comps
    σ = model.params.sigma.values
    return @f(ρ_S)*π/6*∑(@f(x_S,i)*@f(x_S,j)*σ[i,j]^3 for i ∈ comps for j ∈ comps)
end

function a_3(model::SAFTVRMieModel, V, T, z)
    comps = @comps
    x = z/∑(z)
    m = model.params.segment.values
    m̄ = ∑(x .* m)
    return m̄/T^3*∑(@f(x_S,i)*@f(x_S,j)*@f(a_3,i,j) for i ∈ comps for j ∈ comps)
end

function a_3(model::SAFTVRMieModel, V, T, z, i, j)
    ϵ = model.params.epsilon.values
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    ζst_ = ζst(model, V, T, z)
    α = @f(C,i,j)*(1/(λa[i,j]-3)-1/(λr[i,j]-3))
    return -ϵ[i,j]^3*@f(f,α,4)*ζst_ * exp(@f(f,α,5)*ζst_+@f(f,α,6)*ζst_^2)
end

function a_chain(model::SAFTVRMieModel, V, T, z)
    x = z/∑(z)
    m = model.params.segment.values
    return -∑(x[i]*(log(@f(g_Mie,i))*(m[i]-1)) for i ∈ @comps)
end

function g_Mie(model::SAFTVRMieModel, V, T, z, i)
    ϵ = model.params.epsilon.diagvalues
    g_HSi = @f(g_HS,i)
    return g_HSi*exp(ϵ[i]/T*@f(g_1,i)/g_HSi+(ϵ[i]/T)^2*@f(g_2,i)/g_HSi);
end

function g_HS(model::SAFTVRMieModel, V, T, z, i)
    x_0ij = @f(x_0,i,i)
    ζ_X_ = @f(ζ_X)
    k_0 = -log(1-ζ_X_)+(42ζ_X_-39ζ_X_^2+9ζ_X_^3-2ζ_X_^4)/(6*(1-ζ_X_)^3)
    k_1 = (ζ_X_^4+6*ζ_X_^2-12*ζ_X_)/(2*(1-ζ_X_)^3)
    k_2 = -3*ζ_X_^2/(8*(1-ζ_X_)^2)
    k_3 = (-ζ_X_^4+3*ζ_X_^2+3*ζ_X_)/(6*(1-ζ_X_)^3)
    return exp(k_0+x_0ij*k_1+x_0ij^2*k_2+x_0ij^3*k_3)
end

function g_1(model::SAFTVRMieModel, V, T, z, i)
    λr = model.params.lambda_r.diagvalues
    λa = model.params.lambda_a.diagvalues
    x_0ij = @f(x_0,i,i)
    return 3*@f(∂a_1╱∂ρ_S,i)-@f(C,i,i)*(λa[i]*x_0ij^λa[i]*(@f(aS_1,λa[i])+@f(B,λa[i],x_0ij))-λr[i]*x_0ij^λr[i]*(@f(aS_1,λr[i])+@f(B,λr[i],x_0ij)))
end

function ∂a_1╱∂ρ_S(model::SAFTVRMieModel, V, T, z, i)
    λr  = model.params.lambda_r.diagvalues
    λa  = model.params.lambda_a.diagvalues
    x_0ij = @f(x_0,i,i)
    return @f(C,i,i)*(x_0ij^λa[i]*(@f(∂aS_1╱∂ρ_S,λa[i])+@f(∂B╱∂ρ_S,λa[i],x_0ij))
                      - x_0ij^λr[i]*(@f(∂aS_1╱∂ρ_S,λr[i])+@f(∂B╱∂ρ_S,λr[i],x_0ij)))
end

function ∂aS_1╱∂ρ_S(model::SAFTVRMieModel, V, T, z, λ)
    A  = SAFTVRMieconsts.A
    ζ_X_ = @f(ζ_X)
    ∂ζeff╱∂ρ_S = A * [1; 1/λ; 1/λ^2; 1/λ^3] ⋅ [1; 2ζ_X_; 3ζ_X_^2; 4ζ_X_^3] * ζ_X_/@f(ρ_S)
    ζeff_  = @f(ζeff,λ)
    return -1/(λ-3)*((1-ζeff_/2)/(1-ζeff_)^3
                     + @f(ρ_S)*((3*(1-ζeff_/2)*(1-ζeff_)^2
                                 - 0.5*(1-ζeff_)^3)/(1-ζeff_)^6)*∂ζeff╱∂ρ_S);
end

function ∂B╱∂ρ_S(model::SAFTVRMieModel, V, T, z, λ, x_0)
    I = (1-x_0^(3-λ))/(λ-3)
    J = (1-(λ-3)*x_0^(4-λ)+(λ-4)*x_0^(3-λ))/((λ-3)*(λ-4))
    ζ_X_ = @f(ζ_X)
    ρ_S_ = @f(ρ_S)
    return ( ((1-ζ_X_/2)*I/(1-ζ_X_)^3-9*ζ_X_*(1+ζ_X_)*J/(2*(1-ζ_X_)^3))
            + ζ_X_*( (3*(1-ζ_X_/2)*(1-ζ_X_)^2
                - 0.5*(1-ζ_X_)^3)*I/(1-ζ_X_)^6
                - 9*J*((1+2*ζ_X_)*(1-ζ_X_)^3
                + ζ_X_*(1+ζ_X_)*3*(1-ζ_X_)^2)/(2*(1-ζ_X_)^6) ) );
end

function g_2(model::SAFTVRMieModel,V, T, z, i)
    return (1+@f(γ_c,i))*@f(gMCA_2,i)
end

function γ_c(model::SAFTVRMieModel,V, T, z, i)
    ϕ = SAFTVRMieconsts.ϕ
    ϵ = model.params.epsilon.diagvalues
    λr = model.params.lambda_r.diagvalues
    λa = model.params.lambda_a.diagvalues
    ζst_ = @f(ζst)
    α = @f(C,i,i)*(1/(λa[i]-3)-1/(λr[i]-3))
    return ϕ[1][7]*(1-tanh(ϕ[2][7]*(ϕ[3][7]-α)))*ζst_*(exp(ϵ[i]/T)-1)*exp(ϕ[4][7]*ζst_+ϕ[5][7]*ζst_^2)
end

function gMCA_2(model::SAFTVRMieModel, V, T, z, i)
    λr  = model.params.lambda_r.diagvalues
    λa  = model.params.lambda_a.diagvalues
    x_0ij = @f(x_0,i,i)
    ζ_X_  = @f(ζ_X)
    return 3*@f(∂a_2╱∂ρ_S,i)-@f(KHS)*@f(C,i,i)^2 *
    ( λr[i]*x_0ij^(2*λr[i])*(@f(aS_1,2*λr[i])+@f(B,2*λr[i],x_0ij))-
        (λa[i]+λr[i])*x_0ij^(λa[i]+λr[i])*(@f(aS_1,λa[i]+λr[i])+@f(B,λa[i]+λr[i],x_0ij))+
        λa[i]*x_0ij^(2*λa[i])*(@f(aS_1,2*λa[i])+@f(B,2*λa[i],x_0ij)))
end

function ∂a_2╱∂ρ_S(model::SAFTVRMieModel,V, T, z, i)
    λr = model.params.lambda_r.diagvalues
    λa = model.params.lambda_a.diagvalues
    x_0ij = @f(x_0,i,i)
    ζ_X_ = @f(ζ_X)
    ρ_S_ = @f(ρ_S)
    ∂KHS╱∂ρ_S = -ζ_X_/ρ_S_ *
    ( (4*(1-ζ_X_)^3*(1+4*ζ_X_+4*ζ_X_^2-4*ζ_X_^3+ζ_X_^4)
        + (1-ζ_X_)^4*(4+8*ζ_X_-12*ζ_X_^2+4*ζ_X_^3))/(1+4*ζ_X_+4*ζ_X_^2-4*ζ_X_^3+ζ_X_^4)^2 )
    return 0.5*@f(C,i,i)^2 *
    (@f(ρ_S)*∂KHS╱∂ρ_S*(x_0ij^(2*λa[i])*(@f(aS_1,2*λa[i])+@f(B,2*λa[i],x_0ij))
                         - 2*x_0ij^(λa[i]+λr[i])*(@f(aS_1,λa[i]+λr[i])+@f(B,λa[i]+λr[i],x_0ij))
                         + x_0ij^(2*λr[i])*(@f(aS_1,2*λr[i])+@f(B,2*λr[i],x_0ij)))
        + @f(KHS)*(x_0ij^(2*λa[i])*(@f(∂aS_1╱∂ρ_S,2*λa[i])+@f(∂B╱∂ρ_S,2*λa[i],x_0ij))
              - 2*x_0ij^(λa[i]+λr[i])*(@f(∂aS_1╱∂ρ_S,λa[i]+λr[i])+@f(∂B╱∂ρ_S,λa[i]+λr[i],x_0ij))
              + x_0ij^(2*λr[i])*(@f(∂aS_1╱∂ρ_S,2*λr[i])+@f(∂B╱∂ρ_S,2*λr[i],x_0ij))))
end

function a_assoc(model::SAFTVRMieModel, V, T, z)
    x = z/∑(z)
    X_ = @f(X)
    n = model.sites.n_sites
    return ∑(x[i]*∑(n[i][a]*(log(X_[i][a])+(1-X_[i][a])/2) for a ∈ @sites(i)) for i ∈ @comps)
end

function X(model::SAFTVRMieModel, V, T, z)
    _1 = one(V+T+first(z))
    ∑z = ∑(z)
    x = z/∑z
    ρ = N_A*∑z/V
    n = model.sites.n_sites
    itermax = 500
    dampingfactor = 0.5
    error = 1.
    tol = model.absolutetolerance
    iter = 1
    X_ = [[_1 for a ∈ @sites(i)] for i ∈ @comps]
    X_old = deepcopy(X_)
    while error > tol
        iter > itermax && throw("X has failed to converge after " * string(itermax) * " iterations")
        for i ∈ @comps, a ∈ @sites(i)
            rhs = 1/(1+∑(ρ*x[j]*∑(n[j][b]*X_old[j][b]*@f(Δ,i,j,a,b) for b ∈ @sites(j)) for j ∈ @comps))
            X_[i][a] = (1-dampingfactor)*X_old[i][a] + dampingfactor*rhs
        end
        error = sqrt(∑(∑((X_[i][a] - X_old[i][a])^2 for a ∈ @sites(i)) for i ∈ @comps))
        for i = 1:length(X_)
            X_old[i] .= X_[i]
        end
        iter += 1
    end
    return X_
end

function Δ(model::SAFTVRMieModel, V, T, z, i, j, a, b)
    σ = model.params.sigma.values
    ϵ = model.params.epsilon.values
    σ3_x = ∑(∑(@f(x_S,i)*@f(x_S,j)*σ[i,j]^3 for j ∈ @comps) for i ∈ @comps)
    ρR = @f(ρ_S)*σ3_x
    TR = T/ϵ[i,j]
    c  = SAFTVRMieconsts.c
    I = ∑(∑(c[n+1,m+1]*ρR^n*TR^m for m ∈ 0:(10-n)) for n ∈ 0:10)

    ϵ_assoc = model.params.epsilon_assoc.values
    K = model.params.bondvol.values
    F = (exp(ϵ_assoc[i,j][a,b]/T)-1)
    return F*K[i,j][a,b]*I
end

const SAFTVRMieconsts = (
    u = [0.26356031971814109102031,1.41340305910651679221800,3.59642577104072208122300,7.08581000585883755692200,12.6408008442757826594300],
    w = [0.5217556105828086524759,0.3986668110831759274500,7.5942449681707595390e-2,3.6117586799220484545e-3,2.3369972385776227891e-5],

    A = [0.81096   1.7888  -37.578   92.284;
    1.02050  -19.341   151.26  -463.50;
    -1.90570   22.845  -228.14   973.92;
    1.08850  -6.1962   106.98  -677.64],

    ϕ = [[7.5365557, -359.440,  1550.9, -1.199320, -1911.2800,  9236.9,  10.0],
        [-37.604630,  1825.60, -5070.1,  9.063632,  21390.175, -129430,  10.0],
        [71.745953, -3168.00,  6534.6, -17.94820, -51320.700,  357230,  0.57],
        [-46.835520,  1884.20, -3288.7,  11.34027,  37064.540, -315530, -6.70],
        [-2.4679820,- 0.82376, -2.7171,  20.52142,  1103.7420,  1390.2, -8.00],
        [-0.5027200, -3.19350,  2.0883, -56.63770, -3264.6100, -4518.2,   NaN],
        [8.0956883,  3.70900,  0.0000,  40.53683,  2556.1810,  4241.6,   NaN]],

    c  = [0.0756425183020431	-0.128667137050961	 0.128350632316055	-0.0725321780970292	   0.0257782547511452  -0.00601170055221687	  0.000933363147191978  -9.55607377143667e-05  6.19576039900837e-06 -2.30466608213628e-07 3.74605718435540e-09
          0.134228218276565	    -0.182682168504886 	 0.0771662412959262	-0.000717458641164565 -0.00872427344283170	0.00297971836051287	 -0.000484863997651451	 4.35262491516424e-05 -2.07789181640066e-06	4.13749349344802e-08 0
         -0.565116428942893	     1.00930692226792   -0.660166945915607	 0.214492212294301	  -0.0388462990166792	0.00406016982985030	 -0.000239515566373142	 7.25488368831468e-06 -8.58904640281928e-08	0	                 0
         -0.387336382687019	    -0.211614570109503	 0.450442894490509	-0.176931752538907	   0.0317171522104923  -0.00291368915845693	  0.000130193710011706  -2.14505500786531e-06  0	                0	                 0
          2.13713180911797	    -2.02798460133021 	 0.336709255682693	 0.00118106507393722  -0.00600058423301506	0.000626343952584415 -2.03636395699819e-05	 0	                   0	                0	                 0
         -0.300527494795524	     2.89920714512243   -0.567134839686498	 0.0518085125423494	  -0.00239326776760414	4.15107362643844e-05  0	                     0	                   0	                0                    0
         -6.21028065719194	    -1.92883360342573	 0.284109761066570	-0.0157606767372364	   0.000368599073256615	0 	                  0	                     0	                   0	                0	                 0
          11.6083532818029	     0.742215544511197  -0.0823976531246117	 0.00186167650098254   0	                0	                  0	                     0	                   0	                0	                 0
         -10.2632535542427	    -0.125035689035085	 0.0114299144831867	 0	                   0	                0	                  0	                     0	                   0	                0	                 0
          4.65297446837297	    -0.00192518067137033 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0
         -0.867296219639940	     0	                 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0],
)


text:
#test fast transform via mapped Chebyschev
using ApproxFun, FastTransforms, BenchmarkTools

import FastTransforms: Butterfly

M=128

# Get the actual matrix from the plan in ApproxFun.plan_transform(::Hermite(),rand(256)), say.

PlanHerm = ApproxFun.plan_transform(Hermite(),rand(M));
x,w = PlanHerm.plan;

V=ApproxFun.hermitep(0:M-1,x)';
nrm=(V.^2)*w;
P= Diagonal(inv.(nrm))*V*Diagonal(w)

#P=Q[1:N,:]
BF = Butterfly(P, floor(Int, log2(size(P, 1))-6); sketch = :none, atol = 1e-14);

r = FastTransforms.allranks(BF) # just diagnostic info.

mean(r), std(r)
y = zeros(M)+im*zeros(M)
x = randn(M)+im*randn(M);
zr = zeros(M);
zi = zeros(M)
z = zeros(x)
@time A_mul_B!(y, P, x)

@time A_mul_B!(zr, BF, real(x)); A_mul_B!(zi, BF, imag(x));z=zr+im*zi

norm(y-z) # Did it work?

#test for transform to quadrature grid
N=256
x,wx,Tx = nfieldtrans("Hermite",N,2)

BFT = Butterfly(Tx, floor(Int, log2(size(Tx, 1))-6); sketch = :none, atol = 1e-32)

r = FastTransforms.allranks(BFT) # just diagnostic info.

mean(r), std(r)

M=N
x0 = randn(N)
x = randn(N)+im*randn(N)
x1 = randn(M)+im*randn(M)
y = randn(M)
zr = randn(N)
zi = randn(N)
z = randn(N)+im*randn(N)

@time A_mul_B!(x1, Tx, x)

@time A_mul_B!(zr, BFT, real(x)); A_mul_B!(zi, BFT, imag(x));z=zr+im*zi

norm(x-z) # Did it work?


text:
function declare_old_library_product(product_name, product_soname)
    return esc(quote

    end)
end

function declare_new_library_product(product_name)
    handle_name = Symbol(string(product_name, "_handle"))
    get_path_name = Symbol(string("get_", product_name, "_path"))
    path_name = Symbol(string(product_name, "_path"))
end

macro declare_library_product(product_name, product_soname)
    handle_name = Symbol(string(product_name, "_handle"))
    get_path_name = Symbol(string("get_", product_name, "_path"))
    path_name = Symbol(string(product_name, "_path"))
    @static if VERSION < v"1.6.0-DEV"
        lib_declaration = quote
            # On Julia 1.5-, this must be `const` and must be the SONAME
            const $(product_name) = $(product_soname)
        end
    else
        lib_declaration = quote
            # On Julia 1.6+, this doesn't have to be `const`!  Thanks Jeff!
            $(product_name) = ""
        end
    end

    return excat(
        quote
            # These will be filled in by init_library_product()
            $(handle_name) = C_NULL
            $(path_name) = ""
            function $(get_path_name)()
                return $(path_name)::String
            end
        end,
        lib_declaration,
    )
end

function init_new_library_product(product_name)
    @static if VERSION < v"1.6.0-DEV"
        return nothing
    else
        return quote
            # Initialize non-const variable export with the path to this product
            global $(product_name) = $(Symbol(string(product_name, "_path")))
        end
    end
end

macro init_library_product(product_name, product_path, dlopen_flags)
    handle_name = Symbol(string(product_name, "_handle"))
    path_name = Symbol(string(product_name, "_path"))
    return excat(quote
            global $(path_name) = joinpath(artifact_dir, $(product_path))
            # Manually `dlopen()` this right now so that future invocations
            # of `ccall` with its path/SONAME will find this path immediately.
            # dlopen_flags === nothing means to not dlopen the library.
            if $(dlopen_flags) !== nothing
                global $(handle_name) = dlopen($(path_name), $(dlopen_flags))
                push!(LIBPATH_list, joinpath(artifact_dir, $(dirname(product_path))))
            end
        end,
        init_new_library_product(product_name),
    )
end


text:
function description(asset::AbstractAsset)::String

    # get some data for this asset -
    ticker = asset.ticker
    sense = asset.sense

    if (isa(asset, CallContract) == true)

        strike_price = asset.K
        return "$(ticker) CALL@$(strike_price) $(sense)"

    elseif (isa(asset, PutContract) == true)

        strike_price = asset.K
        return "$(ticker) PUT@$(strike_price) $(sense)"

    elseif (isa(asset, Equity) == true)

        C = asset.C
        return "$(ticker) EQUITY@$(C) $(sense)"
    end
end

function payoff(contract::PutContract,
    underlying_price::Float64)::Float64

    # get data from the contract -
    K = contract.K # contract strike price -

    # compute the intrinsic value for a put -
    return max(0.0, K - underlying_price)
end

function payoff(contract::CallContract,
    underlying_price::Float64)::Float64

    # get data from the contract -
    K = contract.K # contract strike price -

    # compute the intrinsic value for a put -
    return max(0.0, underlying_price - K)
end

function payoff(contract::Equity,
    underlying_price::Float64)::Float64

    # get data from this asset -
    C = contract.C

    # return -
    return (C - underlying_price)
end

function compute_profit_loss_at_expiration(assets::Array{AbstractAsset,1},
    underlying_prices::Array{Float64,1})::DataFrame

    # initialize -
    𝒜 = length(assets)
    𝒩 = length(underlying_prices)
    data_array = zeros(𝒩,(𝒜+2))
    col_name_array = Array{String,1}()

    # compute the profit and loss -
    for (asset_index, asset) ∈ enumerate(assets)

        # what is the sense of this asset?
        asset_sense = asset.sense
        C = asset.C  # cost of asset -

        for (price_index, price) ∈ enumerate(underlying_prices)

            # compute the intrinsic value of this contract -
            payoff_value = payoff(asset, price)

            plvalue = 0.0
            if asset_sense == :BUY
                plvalue = payoff_value - C # we purchased this asset/contract
            elseif asset_sense == :SELL
                plvalue = C - payoff_value # we sold this asset/contract
            end

            # capture -
            data_array[price_index,1] = price
            data_array[price_index, asset_index + 1] = plvalue
        end
    end

    for price_index ∈ 1:𝒩
        value = sum(data_array[price_index,2:end-1])
        data_array[price_index,end] = value
    end

    # construct the col names -
    push!(col_name_array,"P")
    for asset ∈ assets
        col_name = build_option_ticker_symbol(asset)
        push!(col_name_array, col_name)
    end
    push!(col_name_array,"Σ")

    # build the data frame -
    df = DataFrame(data_array,col_name_array)

    # return -
    return df
end

max_stars_repo_name: Pooksoft/Serenity

max_issues_repo_name: Alba-Intelligence/Serenity

max_forks_repo_name: Alba-Intelligence/Serenity


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

BLAS.set_num_threads(1)

mm = pyimport("mc_hf.mc_hf")

eri_ep, eri_ep_ao_ints, eia, ejb, e_nocc, e_nvir, p_nocc, p_nvir, mo_coeff, e_nao, p_nao = mm.get_mp2_quantities()

const nao = e_nao + p_nao
const co_e = mo_coeff[ :, 1:e_nocc]
const cv_e = mo_coeff[ :, 1+e_nocc:e_nocc+e_nvir]
const co_n = mo_coeff[ :, 1+e_nocc+e_nvir:e_nocc+e_nvir+p_nocc]
const cv_n = mo_coeff[ :, 1+e_nocc+e_nvir+p_nocc:e_nocc+e_nvir+p_nocc+p_nvir]

function index_2(i,j)
  if(i>j)
    return div(i*(i-1),2)+j
  else
    return div(j*(j-1),2)+i
  end
end

function ep_energy(eri_ep, eia, ejb, p_nvir, p_nocc, e_nvir, e_nocc)
  emp2_ep = 0.0
  for a = 1:p_nvir, i = 1:p_nocc, b = 1:e_nvir, j= 1:e_nocc
           @inbounds emp2_ep += eri_ep[j,b,i,a]^2 / (eia[j,b] + ejb[i,a])
  end
  emp2_ep = 2.0 * emp2_ep
  return emp2_ep
end

function half_loop1!(tmp, i, j,  e_nao, p_nao, eri_ep_ao_ints, x1, y1, z1, mo3, mo4, indx3, indx4)
  ij = index_2(i,j)
  kl=0
   for k = 1:p_nao
    for l = 1:k
      kl = index_2(e_nao+k,e_nao+l)
      ijkl = index_2(ij,kl)
      @inbounds x1[k,l] = eri_ep_ao_ints[ijkl]
      @inbounds x1[l,k] = x1[k,l]
    end
  end

  @inbounds mul!(y1, mo3,x1)
  @inbounds mul!(z1, y1, mo4)
  @inbounds tmp[:,ij] = vec(z1)

end

function half_loop2!(tei_mo, k, l, e_nao, mo1, mo2, indx1, indx2, indx3, indx4, x2,y2,z2, tmp)
  kl = (k-1)*indx4 + l
  ij = 0
   for i in 1:e_nao
     for j in 1:i
      ij = ij + 1
      @inbounds x2[i,j] = tmp[kl,ij]
      @inbounds x2[j,i] = x2[i,j]
    end
    end

  @inbounds mul!(y2, mo1, x2)
  @inbounds mul!(z2, y2, mo2)
  @inbounds tei_mo[:,:,k,l] = z2

end

function fast_conversion_general(nao, e_nao, p_nao, e_nocc, e_nvir, p_nocc, p_nvir, mo1, mo2, mo3, mo4, eri_ep_ao_ints )

  x1 = zeros(p_nao,p_nao,nthreads())
  y1 = zeros(p_nocc,p_nao,nthreads())
  z1 = zeros(p_nocc,p_nvir,nthreads())

  x2 = zeros(e_nao,e_nao,nthreads())
  y2 = zeros(e_nocc,e_nao,nthreads())
  z2 = zeros(e_nocc,e_nvir,nthreads())

  indx1 = size(mo1)[2]
  indx2 = size(mo2)[2]
  indx3 = size(mo3)[2]
  indx4 = size(mo4)[2]

  tei_mo = zeros(indx1, indx2, indx3, indx4)

  mo1_subset = transpose(mo1[1:e_nao,:])
  mo2_subset = mo2[1:e_nao,:]
  mo3_subset = transpose(mo3[1+e_nao:nao,:])
  mo4_subset = mo4[1+e_nao:nao,:]

  value1 = indx3*indx4
  value2 = div(nao*(nao+1),2)

  tmp = zeros(value1,value2)

  local n = threadid()

  @sync for i = 1:nao
     Threads.@spawn for j = 1:i
      half_loop1!(tmp, i, j, e_nao, p_nao, eri_ep_ao_ints, x1[ :, :, n], y1[ :, :,n], z1[ :,:, n], mo3_subset, mo4_subset, indx3, indx4)
     end
  end

  @sync for k in 1:indx3
    Threads.@spawn for l in 1:indx4
      half_loop2!(tei_mo, k, l, e_nao, mo1_subset, mo2_subset, indx1, indx2, indx3, indx4, x2[:, :,n], y2[:,:,n], z2[:,:,n], tmp)
    end
  end

  return tei_mo

end

@btime mo_int4 = fast_conversion_general(nao, e_nao, p_nao, e_nocc, e_nvir, p_nocc, p_nvir, co_e, cv_e, co_n, cv_n, eri_ep_ao_ints)
mo_int4 = fast_conversion_general(nao, e_nao, p_nao, e_nocc, e_nvir, p_nocc, p_nvir, co_e, cv_e, co_n, cv_n, eri_ep_ao_ints)

mp2_ep = ep_energy(mo_int4, eia, ejb, p_nvir, p_nocc, e_nvir, e_nocc)
println(mp2_ep)


text:
"""
Internals, run all callbacks
"""
function _run_callbacks(callbacks::C, scenes::Vector{Scene{Entity{S,D,I}}}, actions::Union{Nothing, Vector{Scene{A}}}, roadway::R, models::Dict{I,M}, tick::Int) where {S,D,I,A<:EntityAction,R,M<:DriverModel,C<:Tuple{Vararg{Any}}}
    isdone = false
    for callback in callbacks
        isdone |= run_callback(callback, scenes, actions, roadway, models, tick)
    end
    return isdone
end

"""
    run_callback(callback, scenes::Vector{EntityScene}, actions::Union{Nothing, Vector{A}}, roadway::Roadway, models::Dict{I, DriverModel}, tick::Int64)

Given a callback type, `run_callback` will be run at every step of a simulation run using `simulate`.
By overloading the `run_callback` method for a custom callback type one can log information or interrupt a simulation.
The `run_callback` function is expected to return a boolean. If `true` the simulation is stopped.

# Inputs:
- `callback` the custom callback type used for dispatch
- `scenes` where the simulation data is stored, note that it is only filled up to the current time step (`scenes[1:tick+1]`)
- `actions` where the actions are stored, it is only filled up to `actions[tick]`
- `roadway` the roadway where entities are moving
- `models` a dictionary mapping entity IDs to driver models
- `tick` the index of the current time step
"""
function run_callback end

## Implementations of useful callbacks

"""
    CollisionCallback

Terminates the simulation once a collision occurs
"""
@with_kw struct CollisionCallback
    mem::CPAMemory=CPAMemory()
end

function run_callback(
    callback::CollisionCallback,
    scenes::Vector{Scene{E}},
    actions::Union{Nothing, Vector{Scene{A}}},
    roadway::R,
    models::Dict{I,M},
    tick::Int
) where {E<:Entity,A<:EntityAction,R,I,M<:DriverModel}
    return !is_collision_free(scenes[tick], callback.mem)
end


text:
module PkgOptionsLoader

using Base: PkgId, package_slug, slug, _crc32c
using Pkg: TOML

macro return_if_nothing(ex)
    quote
        ans = $(esc(ex))
        ans === nothing && return nothing
        ans
    end
end

project_slug() = slug(_crc32c(@return_if_nothing Base.active_project()), 5)

optionsdir() =
    joinpath(Base.DEPOT_PATH[1], "options", @return_if_nothing project_slug())

optionspath(pkg::PkgId) =
    joinpath(@return_if_nothing(optionsdir()),
             package_slug(pkg.uuid),
             pkg.name * ".toml")

function load(mod)
    pkg = PkgId(mod)
    if pkg.uuid === nothing
        error("Cannot load package options in a non-package module ", mod)
    end

    tomlpath = optionspath(pkg)
    if tomlpath === nothing || !isfile(tomlpath)
        return Dict{String,Any}()
    end

    # Editing this TOML file should invalidate the cache:
    Base.include_dependency(tomlpath)

    return TOML.parsefile(tomlpath)
end

"""
    PkgOptionsLoader.@load

Load package options for current module at _expansion time_.  It
evaluates to a Dict.  It is recommended to load it at the top level of
your package module:

```julia
const pkgoptions = PkgOptionsLoader.@load
```

It also is a good idea to verify options and then assign it to a `const`:

```julia
get!(pkgoptions, "FEATURE_FLAG", false)  # default

if !(pkgoptions["FEATURE_FLAG"] isa Bool)
    error("some helpful error message")
end

const FEATURE_FLAG::Bool = pkgoptions["FEATURE_FLAG"]
```
"""
macro load()
    load(__module__)
end

end # module


text:
module MusicManipulations

using Reexport
@reexport using MIDI
import MIDI: Notes, Note

include("general.jl")
include("midifiles.jl")
include("quantize.jl")
include("scale.jl")
include("data_handling/data_extraction.jl")
include("data_handling/timeseries.jl")

include("specific_modules/jazz.jl")
include("specific_modules/drums.jl")

include("motifs/notes_with_motifs.jl")

include("staff_printing/musescore.jl")
import .MuseScore: musescore;

export MuseScore, musescore

# using Requires
# @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" begin
#     include("visualize.jl")
# end

export Jazz, Drums

end


text:
module PlantBiophysics

# For model parameters (efficient and still mutable!)
using MutableNamedTuples

# For IO:
using YAML
using CSV
using OrderedCollections

using Dates

import DataFrames.DataFrame # For convenience transformations
import DataFrames.Not
import DataFrames.rename!
import DataFrames: select!, select
import DataFrames.dropmissing!
import Base.show
import Base.getindex
import Impute.locf # For filling missing values in comments in Walz to identify curves
import LsqFit: curve_fit
using RecipesBase
import Statistics: mean

# Generic structures:
include("structs/Abstract_model_structs.jl")
include("structs/concrete_component_structs.jl")

# Physical constants:
include("structs/constants.jl")

# Atmosphere computations (vapor pressure...)
include("structs/atmosphere.jl")

# Models helpers:
include("models_helpers.jl")

# Conversions
include("conversions.jl")

# Light interception
include("light_interception/generic_structs.jl")
include("light_interception/Ignore.jl")
include("light_interception/Translucent.jl")

# Photosynthesis related files:
include("photosynthesis/photosynthesis.jl")
include("photosynthesis/constantA.jl")
include("photosynthesis/FvCB.jl")
include("photosynthesis/FvCBIter.jl")
include("photosynthesis/FvCBRaw.jl")
include("photosynthesis/temperature-dependence.jl")

# Stomatal conductance related files:
include("conductances/stomatal/stomatal_conductance.jl")
include("conductances/stomatal/constantGs.jl")
include("conductances/stomatal/medlyn.jl")

# Boundary layer conductance:
include("conductances/boundary/gb.jl")

# Energy balance
include("energy/longwave_energy.jl")
include("energy/energy_balance.jl")
include("energy/Missing.jl")
include("energy/Monteith.jl")

# File IO
include("io/read_model.jl")
include("io/read_weather.jl")
include("io/read_walz.jl")

# Parameters optimization
include("fitting/fit.jl")
include("fitting/fit_FvCB.jl")

# Model evaluation
include("evaluation/statistics.jl")

# Depreciations
include("depreciations/models.jl")

# File IO:
export read_model
export is_model
export read_weather
export read_walz

# Physical constants
export Constants

# Status
export init_status!

# Get models informations
export variables, variables_typed
export inputs, outputs, defaults
export to_initialise, is_initialised

# Atmosphere
export vapor_pressure
export e_sat
export e_sat_slope
export air_density
export Atmosphere
export Weather

# Conversions
export rh_from_vpd
export rh_from_e
export ms_to_mol
export mol_to_ms
export gbh_to_gbw
export gbw_to_gbh
export gsw_to_gsc
export gsc_to_gsw

# Energy balance
export AbstractEnergyModel
export black_body
export grey_body
export psychrometer_constant
export net_longwave_radiation
export energy_balance # main interface to user
export energy_balance! # main interface to user
export Monteith       # a struct to hold the values for the model of Monteith and Unsworth (2013)
export latent_heat
export sensible_heat
export γ_star
export latent_heat_vaporization

# structure for light interception
export Translucent
export Ignore
export OpticalProperties
export σ
export AbstractInterceptionModel

# Photosynthesis
export AbstractAModel
export ConstantA
export Fvcb # Parameters for the coupled Farquhar et al. (1980) model
export FvcbIter # To update...
export FvcbRaw # Parameters for the original Farquhar et al. (1980) model
export Constants
export photosynthesis!
export photosynthesis
export photosynthesis!_

# Conductances
export AbstractGsModel
export gbh_to_gbw
export gbₕ_free
export gbₕ_forced
export gs, gs!
export Medlyn
export ConstantGs

# Temporary structures (to move to another package)
export Translucent
export Ignore
export get_km, Γ_star, arrhenius, get_J, gs_closure, get_Cᵢⱼ, get_Cᵢᵥ, get_Dₕ
export init_variables_manual, init_variables, Fvcb_net_assimiliation
export get_component_type, get_process, get_model, instantiate, get_component_type

export AbstractModel

# Components (structures that hold models)
export AbstractComponentModel
export LeafModels

# Parameters optimization
export fit

# Convenience functions
export DataFrame, copy, getindex

# Model evaluation
export EF, RMSE, dr

end