text:
using ForwardDiff

"Parameters for right-hand side function of fully implicit Runge-Kutta methods."
mutable struct ParametersFIRK{DT, TT, D, S, ET <: NamedTuple, FT, JT} <: Parameters{DT,TT}
    equs::ET
    tab::Tableau{TT}
    Δt::TT

    F::FT
    Jconfig::JT

    t::TT
    q::Vector{DT}

    function ParametersFIRK{DT,D}(equs::ET, tab::Tableau{TT}, Δt::TT) where {DT, TT, D, ET <: NamedTuple}
        F = (v,q) -> equs[:v](zero(TT), q, v)
        tq = zeros(DT,D)
        tv = zeros(DT,D)
        Jconfig = ForwardDiff.JacobianConfig(F, tv, tq)

        new{DT, TT, D, tab.s, ET, typeof(F), typeof(Jconfig)}(equs, tab, Δt, F, Jconfig, zero(TT), zeros(DT,D))
    end
end

@doc raw"""
Fully implicit Runge-Kutta integrator cache.

### Fields

* `q̃`: initial guess of solution
* `ṽ`: initial guess of vector field
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of solution
* `V`: internal stages of vector field
* `Y`: vector field of internal stages
"""
struct IntegratorCacheFIRK{DT,D,S} <: ODEIntegratorCache{DT,D}
    q̃::Vector{DT}
    ṽ::Vector{DT}
    s̃::Vector{DT}

    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    Y::Vector{Vector{DT}}
    J::Vector{Matrix{DT}}

    function IntegratorCacheFIRK{DT,D,S}() where {DT,D,S}
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        Y = create_internal_stage_vector(DT, D, S)
        J = [zeros(DT,D,D) for i in 1:S]
        new(zeros(DT,D), zeros(DT,D), zeros(DT,D), Q, V, Y, J)
    end
end

function IntegratorCache{ST}(params::ParametersFIRK{DT,TT,D,S}; kwargs...) where {ST,DT,TT,D,S}
    IntegratorCacheFIRK{ST,D,S}(; kwargs...)
end

@inline CacheType(ST, params::ParametersFIRK{DT,TT,D,S}) where {DT,TT,D,S} = IntegratorCacheFIRK{ST,D,S}

@doc raw"""
Fully implicit Runge-Kutta integrator solving the system
```math
\begin{aligned}
V_{n,i} &= v (Q_{n,i}, P_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} .
\end{aligned}
```
"""
struct IntegratorFIRK{DT, TT, D, S, PT <: ParametersFIRK{DT,TT},
                                    ST <: NonlinearSolver{DT},
                                    IT <: InitialGuessODE{TT}} <: AbstractIntegratorRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorFIRK(params::ParametersFIRK{DT,TT,D,S}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,ST,IT}
        new{DT, TT, D, S, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorFIRK{DT,D}(equations::NamedTuple, tableau::Tableau{TT}, Δt::TT; exact_jacobian=true) where {DT,TT,D}
        # get number of stages
        S = tableau.s

        # check if tableau is fully implicit
        if get_config(:verbosity) ≥ 1
            if isexplicit(tableau)
                @warn "Initializing IntegratorFIRK with explicit tableau $(tableau.name).\nYou might want to use IntegratorERK instead."
            elseif isdiagnonallyimplicit(tableau)
                @warn "Initializing IntegratorFIRK with diagonally implicit tableau $(tableau.name).\nYou might want to use IntegratorDIRK instead."
            end
        end

        # create params
        params = ParametersFIRK{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        if exact_jacobian
            solver = create_nonlinear_solver_with_jacobian(DT, D*S, params, caches)
        else
            solver = create_nonlinear_solver(DT, D*S, params, caches)
        end

        # create initial guess
        iguess = InitialGuessODE(get_config(:ig_extrapolation), equations[:v], Δt)

        # create integrator
        IntegratorFIRK(params, solver, iguess, caches)
    end

    function IntegratorFIRK{DT,D}(v::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorFIRK{DT,D}(NamedTuple{(:v,)}((v,)), tableau, Δt; kwargs...)
    end

    function IntegratorFIRK{DT,D}(v::Function, h::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorFIRK{DT,D}(NamedTuple{(:v,:h)}((v,h)), tableau, Δt; kwargs...)
    end

    function IntegratorFIRK(equation::ODE{DT}, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorFIRK{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

@inline Base.ndims(::IntegratorFIRK{DT,TT,D,S}) where {DT,TT,D,S} = D

function initialize!(int::IntegratorFIRK, sol::AtomicSolutionODE)
    sol.t̄ = sol.t - timestep(int)

    equations(int)[:v](sol.t, sol.q, sol.v)

    initialize!(int.iguess, sol.t, sol.q, sol.v,
                            sol.t̄, sol.q̄, sol.v̄)
end

function update_params!(int::IntegratorFIRK, sol::AtomicSolutionODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
end

function initial_guess!(int::IntegratorFIRK{DT}, sol::AtomicSolutionODE{DT},
                        cache::IntegratorCacheFIRK{DT}=int.caches[DT]) where {DT}

    local offset::Int

    # compute initial guess for internal stages
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q̄, sol.v̄, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).c[i])
    end
    for i in eachstage(int)
        offset = ndims(int)*(i-1)
        for k in eachdim(int)
            int.solver.x[offset+k] = 0
            for j in eachstage(int)
                int.solver.x[offset+k] += timestep(int) * tableau(int).a[i,j] * cache.V[j][k]
            end
        end
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
                         params::ParametersFIRK{DT,TT,D}) where {ST,DT,TT,D}

    local tᵢ::TT

    # copy x to Y and compute Q = q + Δt Y
    for i in eachindex(Q,Y)
        for k in eachindex(Q[i],Y[i])
            Y[i][k] = x[D*(i-1)+k]
            Q[i][k] = params.q[k] + Y[i][k]
        end
    end

    # compute V = v(Q)
    for i in eachindex(Q,V)
        tᵢ = params.t + params.Δt * params.tab.c[i]
        params.equs[:v](tᵢ, Q[i], V[i])
    end
end

# Compute stages of fully implicit Runge-Kutta methods.
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersFIRK{DT,TT,D},
                          caches::CacheDict) where {ST,DT,TT,D}
    # temporary variables
    local y1::ST
    local y2::ST

    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache.Q, cache.V, cache.Y, params)

    # compute b = - (Y-AV)
    for i in eachindex(cache.Y)
        for k in eachindex(cache.Y[i])
            y1 = 0
            y2 = 0
            for j in eachindex(cache.V)
                y1 += params.tab.a[i,j] * cache.V[j][k]
                y2 += params.tab.â[i,j] * cache.V[j][k]
            end
            b[D*(i-1)+k] = - cache.Y[i][k] + params.Δt * (y1 + y2)
        end
    end
end

# Compute Jacobian of fully implicit Runge-Kutta methods.
function jacobian!(x::Vector{DT}, jac::Matrix{DT}, cache::IntegratorCacheFIRK{DT,D,S},
                   params::ParametersFIRK{DT,TT,D,S}) where {DT,TT,D,S}
    local tᵢ::TT

    for i in 1:S
        for k in 1:D
            cache.Y[i][k] = x[D*(i-1)+k]
            cache.Q[i][k] = params.q[k] + cache.Y[i][k]
        end
        tᵢ = params.t + params.Δt * params.tab.c[i]
        F = (v,q) -> params.equs[:v](tᵢ, q, v)
        ForwardDiff.jacobian!(cache.J[i], params.F, cache.ṽ, cache.Q[i], params.Jconfig)
    end

    jac .= 0

    @inbounds for j in 1:S
        for i in 1:S
            for l in 1:D
                for k in 1:D
                    m = D*(i-1)+k
                    n = D*(j-1)+l

                    if i == j && k == l
                        jac[m,n] += -1
                    end

                    jac[m,n] += params.Δt * params.tab.a[i,j] * cache.J[j][k,l]
                    jac[m,n] += params.Δt * params.tab.â[i,j] * cache.J[j][k,l]
                end
            end
        end
    end
end

function integrate_step!(int::IntegratorFIRK{DT,TT}, sol::AtomicSolutionODE{DT,TT},
                         cache::IntegratorCacheFIRK{DT}=int.caches[DT]) where {DT,TT}

    # update nonlinear solver parameters from atomic solution
    update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset atomic solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector field at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.V, cache.Y, int.params)

    # compute final update
    update_solution!(sol.q, sol.q̃, cache.V, tableau(int).b, tableau(int).b̂, timestep(int))

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.v)
end


text:
function surfc(gr::Graph, z::Array{cmgl.Float}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surfc(gr.ptr, Data(z).ptr, Data(c).ptr, stl, opt)
end

function surfc(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surfc_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(c).ptr, stl, opt)
end

function surf3c(gr::Graph, v::cmgl.Float, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surf3c_val(gr.ptr, v, Data(a).ptr, Data(c).ptr, stl, opt)
end

function surf3c(gr::Graph, v::cmgl.Float, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surf3c_xyz_val(gr.ptr, v, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(a).ptr, Data(c).ptr, stl, opt)
end

function surf3c(gr::Graph, v::Array{cmgl.Float, 1}, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float, 3}; stl::String="", opt::String="")
    for val in v
        cmgl.surf3c_val(gr.ptr, val, Data(a).ptr, Data(c).ptr, stl, opt)
    end
end

function surf3c(gr::Graph, v::Array{cmgl.Float, 1}, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float, 3}; stl::String="", opt::String="")
    for val in v
        cmgl.surf3c_xyz_val(gr.ptr, val, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(a).ptr, Data(c).ptr, stl, opt)
    end
end

function surf3c(gr::Graph, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surf3c(gr.ptr, Data(a).ptr, Data(c).ptr, stl, opt)
end

function surf3c(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surf3c_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(a).ptr, Data(c).ptr, stl, opt)
end

function surfa(gr::Graph, z::Array{cmgl.Float}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surfa(gr.ptr, Data(z).ptr, Data(c).ptr, stl, opt)
end

function surfa(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surfa_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(c).ptr, stl, opt)
end

function surf3a(gr::Graph, v::cmgl.Float, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surf3a_val(gr.ptr, v, Data(a).ptr, Data(c).ptr, stl, opt)
end

function surf3a(gr::Graph, v::cmgl.Float, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surf3a_xyz_val(gr.ptr, v, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(a).ptr, Data(c).ptr, stl, opt)
end

function surf3a(gr::Graph, v::Array{cmgl.Float, 1}, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float, 3}; stl::String="", opt::String="")
    for val in v
        cmgl.surf3a_val(gr.ptr, val, Data(a).ptr, Data(c).ptr, stl, opt)
    end
end

function surf3a(gr::Graph, v::Array{cmgl.Float, 1}, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float, 3}; stl::String="", opt::String="")
    for val in v
        cmgl.surf3a_xyz_val(gr.ptr, val, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(a).ptr, Data(c).ptr, stl, opt)
    end
end

function surf3a(gr::Graph, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surf3a(gr.ptr, Data(a).ptr, Data(c).ptr, stl, opt)
end

function surf3a(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, a::Array{cmgl.Float, 3}, c::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.surf3a_xyz(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(a).ptr, Data(c).ptr, stl, opt)
end

function tiles(gr::Graph, z::Array{cmgl.Float}, r::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tiles(gr.ptr, Data(z).ptr, Data(r).ptr, stl, opt)
end

function tiles(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, z::Array{cmgl.Float}, r::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.tiles_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(z).ptr, Data(r).ptr, stl, opt)
end

function map(gr::Graph, ax::Array{cmgl.Float}, ay::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.map(gr.ptr, Data(ax).ptr, Data(ay).ptr, stl, opt)
end

function map(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, ax::Array{cmgl.Float}, ay::Array{cmgl.Float}; stl::String="", opt::String="")
    cmgl.map_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(ax).ptr, Data(ay).ptr, stl, opt)
end

function stfa(gr::Graph, real::Array{cmgl.Float}, imag::Array{cmgl.Float}, dn::Int; stl::String="", opt::String="")
    cmgl.stfa(gr.ptr, Data(real).ptr, Data(imag).ptr, dn, stl, opt)
end

function stfa(gr::Graph, x::Array{cmgl.Float}, y::Array{cmgl.Float}, real::Array{cmgl.Float}, imag::Array{cmgl.Float}, dn::Int; stl::String="", opt::String="")
    cmgl.stfa_xy(gr.ptr, Data(x).ptr, Data(y).ptr, Data(real).ptr, Data(imag).ptr, dn, stl, opt)
end

ext: jl


text:
# test interpolation of fields from one SBP operator to another

function test_interp()

  @testset "----- Testing Field Interpolation -----" begin
    fname = "input_vals_channel_dg_large.jl"

    opts = read_input(fname)
    opts2 = copy(opts)
    opts2["order"] += 1

    degree = opts["order"]

    # solve problem on both grids

    mesh, sbp, eqn, opts = createObjects(opts)
    mesh2, sbp2, eqn2, opts2 = createObjects(opts2)
    #array1DTo3D(mesh2, sbp2, eqn2, opts2, eqn2.q_vec, eqn2.q)
    q_vec = copy(eqn.q_vec)
    q_vec2 = copy(eqn2.q_vec)

    # put degree p field into eqn
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        y = mesh.coords[2, j, i]

        for k=1:mesh.numDofPerNode
          val = x^degree + 2*y^degree + k
          eqn.q[k, j, i] = val
          q_vec[mesh.dofs[k, j, i]] = val
        end
      end
    end

    interpField(sbp, eqn.q, sbp2, eqn2.q)
    interpField(mesh, sbp, q_vec, mesh2, sbp2, q_vec2)

    for i=1:mesh2.numEl
      for j=1:mesh2.numNodesPerElement
        x = mesh2.coords[1, j, i]
        y = mesh2.coords[2, j, i]

        for k=1:mesh2.numDofPerNode
          val = x^degree + 2*y^degree + k

          @test isapprox( abs(eqn2.q[k, j, i] - val), 0.0) atol=1e-13
          @test isapprox( abs(q_vec2[mesh2.dofs[k, j, i]] - val), 0.0) atol=1e-13
        end
      end
    end

  end  # end facts block

  return nothing
end

add_func1!(EulerTests, test_interp, [TAG_SHORTTEST])


text:
using Test
using RegisterUtilities

@testset "Counter test" begin
    for empty_gridsize in ((), (0,), (1,0), (1,-1), (1,0,1))
        for c in Counter(empty_gridsize)
            @test c == nothing
        end
    end
    for (gridsize, index_vec) in (((1,),[[1]]), ((2,3),[[1,1],[2,1],[1,2],[2,2],[1,3],[2,3]]),
            ((2,3,4), [[1, 1, 1],[2, 1, 1],[1, 2, 1],[2, 2, 1],[1, 3, 1],[2, 3, 1],
                       [1, 1, 2],[2, 1, 2],[1, 2, 2],[2, 2, 2],[1, 3, 2],[2, 3, 2],
                       [1, 1, 3],[2, 1, 3],[1, 2, 3],[2, 2, 3],[1, 3, 3],[2, 3, 3],
                       [1, 1, 4],[2, 1, 4],[1, 2, 4],[2, 2, 4],[1, 3, 4],[2, 3, 4]]))
        cnt = Array{Array{Int64,1},1}()
        for c in Counter(gridsize)
            push!(cnt,c)
        end
        @test cnt == index_vec
    end
end


text:
module Parse

export parse_slide
using CommonMark
using Pkg
using ..Types

const TOML_NAME = ["JuliaHuan.toml", "Huan.toml"]

function default_features()
    [
        FrontMatterRule(toml=Pkg.TOML.parse),
        FootnoteRule(),
        MathRule(),
        TableRule(),
        RawContentRule(),
        AttributeRule(),
        CitationRule(),
        AutoIdentifierRule(),
        AdmonitionRule(),
    ]
end

function create_parser()
    parser = Parser()
    enable!(parser, default_features())
    return parser
end

function segment(file)
    return open(file) do io
        raw = read(io, String)
        split(raw, "---\n")
    end
end

function default_meta()
    Dict(
        "highlight" => "monokai",
        "theme" => "black",
        "revealjs" => Dict(
            "version" => "4.0.2",
            "path" => "reveal.js",
        )
    )
end

function extract_meta(ast::CommonMark.Node)
    for (node, entering) in ast
        if node.t isa CommonMark.FrontMatter
            return node.t.data
        end
    end
    return
end

function parse_slide(path, parser=create_parser())
    if isfile(path)
        return parse_file(path, parser)
    elseif ispath(path)
        return parse_project(path, parser)
    else
        error("invalid path: $path")
    end
end

function parse_file(filename, parser)
    slides = _parse_file(filename, parser)
    m = merge(default_meta(), meta(first(slides)))
    return Presentation(m, filename, slides, true)
end

function parse_project(path, parser)
    toml = nothing
    for each in TOML_NAME
        if isfile(joinpath(path, each))
            toml = joinpath(path, each)
            break
        end
    end

    toml === nothing && error("please specify meta in $(join(TOML_NAME, " or "))")
    meta = merge(default_meta(), Pkg.TOML.parsefile(toml))

    slides = _parse_file(joinpath(path, "index.md"), parser)
    return Presentation(meta, path, slides, false)
end

function _parse_file(filename, parser)
    return map(segment(filename)) do s
        ast = parser(s)
        meta = extract_meta(ast)

        if meta !== nothing && haskey(meta, "include")
            file = joinpath(dirname(filename), meta["include"])
            isfile(file) || error("$file does not exist!")
            slides = _parse_file(file, parser)
            return SectionSlide(meta, filename, slides)
        else
            return BasicSlide(meta, filename, ast)
        end
    end
end

end


text:
#################
### User Info ###
#################

getUserInfo(c::Client) = getParsedResponse(User,c,Requests.get,"/user")

getUserInfo(c::Client, user::String) = getParsedResponse(User,c,Requests.get,"/users/$(user)")

##########################
### Application Tokens ###
##########################

listAccessTokens(c::Client, user::String,pass::String) = getParsedResponse(Vector{AccessToken},c,Requests.get,"/users/$(user)/tokens"; headers = Dict("Authorization" => "Basic $(base64encode("$(user):$(pass)"))"))

function createAccessToken(c::Client,user::String,pass::String,name::String)
	data = Dict{String,Any}("name"=>name)
	headers = Dict("Authorization" => "Basic $(base64encode("$(user):$(pass)"))")
	getParsedResponse(AccessToken,c,Requests.post,"/users/$(user)/tokens"; headers = headers, json = data)
end

####################
### User E-mails ###
####################

function listEmails(c::Client)
	getParsedResponse(Vector{Email},c,Requests.get,"/user/emails")
end

addEmail(c::Client, emails::String) = addEmail(c,[emails])

function addEmail(c::Client, emails::Vector{String})
	data = Dict{String,Any}("emails"=>emails)
	getParsedResponse(Vector{Email},c,Requests.post,"/user/emails";json = data)
end

deleteEmail(c::Client, emails::String) = deleteEmail(c,[emails])

function deleteEmail(c::Client, emails::Vector{String})
	data = Dict{String,Any}("emails"=>emails)
	getResponse(c,Requests.delete,"/user/emails";json = data)
	nothing
end

###################
### User Follow ###
###################

listMyFollowers(c::Client,page::Int) = getParsedResponse(Vector{User},c,Requests.get,"/user/followers?page=$(page)")

listFollowers(c::Client,user::String, page::Int) = getParsedResponse(Vector{User},c,Requests.get,"/users/$(users)/followers?page=$(page)")

listMyFollowing(c::Client,page::Int) = getParsedResponse(Vector{User},c,Requests.get,"/user/following?page=$(page)")

listFollowing(c::Client,user::String, page::Int) = getParsedResponse(Vector{User},c,Requests.get,"/user/$(users)/following?page=$(page)")

function isFollowing(c::Client,target::String)
	try
		getResponse(c,Requests.get,"/user/following/$(target)")
	catch
		return false
	end
	return true
end

function isUserFollowing(c::Client,user::String, target::String)
	try
		getResponse(c,Requests.get,"/user/$(users)/following/$(target)")
	catch
		return false
	end
	return true
end

follow(c::Client,target::String) = (getResponse(c, Requests.put, "/user/following/$(target)"); nothing)

unfollow(c::Client,target::String) = (getResponse(c, Requests.delete, "/user/following/$(target)"); nothing)

####################
### User GPG Key ###
####################

listGPGKeys(c::Client,user::String) = getParsedResponse(Vector{GPGKey},c,Requests.get,"/users/$(user)/gpg_keys")

listMyGPGKeys(c::Client) = getParsedResponse(Vector{GPGKey},c,Requests.get,"/user/gpg_keys")

getGPGKey(c::Client,id::Int64) = getParsedResponse(GPGKey,c,Requests.get,"/user/gpg_keys/$(id)")

function createPublicKey(c::Client,key::String)
	data = Dict{String,Any}("armored_public_key"=>key)
	return getParsedResponse(PublicKey,c,Requests.post,"/user/gpg_keys";json = data)
end

deleteGPGKey(c::Client,keyID::Int64) = (getResponse(c,Requests.delete,"/user/gpg_keys/$(keyID)"); nothing)

####################
### User SSH Key ###
####################

listPublicKeys(c::Client,user::String) = getParsedResponse(Vector{PublicKey},c,Requests.get,"/users/$(user)/keys")

listMyPublicKeys(c::Client) = getParsedResponse(Vector{PublicKey},c,Requests.get,"/user/keys")

getPublicKey(c::Client,keyID::Int64) = getParsedResponse(PublicKey,c,Requests.get,"/user/keys/$(keyID)")

function createPublicKey(c::Client,title::String,key::String)
	data = Dict{String,Any}("title"=>title,"key"=>key)
	return getParsedResponse(PublicKey,c,Requests.post,"/user/keys";json = data)
end

deletePublicKey(c::Client,keyID::Int64) = (getResponse(c,Requests.delete,"/user/keys/$(keyID)"); nothing)


text:
# TODO: This is very specific to UW student seminars. in long term, try to make it generic

module Talks

using Dates
using Base.Iterators
using TimeZones
using EnglishText
using SExpressions.Lists
using Remarkable.Articles
using Remarkable.Tags

function fromjson(obj, tagmatrix)
    time = obj[:time]
    if length(time) == 2
        time = time[1]  # TODO: support time interval
    end
    properties = ArticleMetadata(
        obj[:identifier],
        obj[:title],
        [obj[:speaker]],
        ZonedDateTime(DateTime(time, "y-m-dTH:M:S"), tz"America/Toronto"))
    tag!(properties, tagmatrix, obj[:tags])
    LivePerformance(properties, obj[:location])
end

url(t) = "archive/$(identifier(t))"
iscompleted(t) = Date(DateTime(datetime(t))) < Dates.today()

abstractpath(talk) = "abstract/$(identifier(talk))"
hasabstract(talk) = isfile(abstractpath(talk))

summarypath(talk) = "summary/$(identifier(talk))"
hassummary(talk) = isfile(summarypath(talk))

# TODO: this has nothing to do with talks, move it to document
function brief(t)
    if haskey(t, :type) && t[:type] == "reference"
        Dict(:title => t[:title],
             :url => "document/$(t[:id])",
             :summary => join(flatten([
                ["This is a reference document on $(t[:title])",
                 "by $(ItemList(t[:authors]))."],
                (haskey(t, :subsetof) && !isempty(t[:subsetof]) ?
                    ["It covers a subset of the material of",
                     "$(ItemList(t[:subsetof], Disjunction()))."] :
                    [])
             ]), " "))
    end
end

export abstractpath, hasabstract, summarypath, hassummary, iscompleted,
       valuate, brief, fromjson, url

end


text:
@testset "list" begin
    list(TEST_CONFIG)
end


text:
const CMDSTAN_HOME_VAR = "JULIA_CMDSTAN_HOME"

"""

Get the path to the `CMDSTAN_HOME` environment variable.

$(SIGNATURES)

# Extended help

Example: `get_cmdstan_home()`
"""
function get_cmdstan_home()
    get(ENV, CMDSTAN_HOME_VAR) do
        throw(ErrorException("The environment variable $CMDSTAN_HOME_VAR needs to be set."))
    end
end

"""

Set the path for `CMDSTAN_HOME`.

$(SIGNATURES)

# Extended help

Example: `set_cmdstan_home!(homedir() * "/Projects/Stan/cmdstan/")`
"""
set_cmdstan_home!(path) = global CMDSTAN_HOME = path


text:
using GRIB

keylist = Vector{String}()
grib_file = "/home/neo/Projekte/Predly/weather/gfs_workflow/download/data/gfs.20210109/06/gfs.t06z.pgrb2.0p25.f002"

parameters = [7,11,586,422,423,424,426,428,430,433,434,435,437,438,439,440,441,442,443,451,453,455,456,457,458,472,473,474,476,479,480,481,482,483,495,496,498,588]

unneeded_identifiers = ["values", "codedValues", "latLonValues", "latitudes", "longitudes", "distinctLatitudes", "distinctLongitudes"]

start_time = time()

data_dict = Dict()

GribFile(grib_file) do f
    for (index, message) in enumerate(f)

        band_time = time()
        data_array = []
        parameter_name = message["parameterName"]
        parameter_unit = message["parameterUnits"]
        parameter_ecmf = message["nameECMF"]
        units_ecmf = message["unitsECMF"]
        name = message["name"]
        units = message["units"]
        level = message["level"]
        type_of_level = message["typeOfLevel"]
        short_name = message["shortName"]
        data_dict[index] = Dict()

        # calculate x/y raster from points where x == 2y - 2
        # x is 1440 here, y is 721, 1440x721=1038240
        number_points = message["numberOfDataPoints"]
        y =  round(sqrt(number_points/2), digits=0, RoundUp)
        x =  2 * y - 2

        lat0 = message["latitudeOfFirstGridPointInDegrees"]
        lon0 = message["longitudeOfFirstGridPointInDegrees"]
        lat1 = message["latitudeOfLastGridPointInDegrees"]
        lon1 = message["longitudeOfLastGridPointInDegrees"]

        if index in parameters

            data_dict[index]["latitudes"] = message["latitudes"]
            data_dict[index]["longitudes"] = message["longitudes"]
            data_dict[index]["values"] = message["values"]

            println("band ", index ," parameter name: ", parameter_name, " unit ", parameter_unit,
                    " name ecmf ", parameter_ecmf, " units ecmf ", units_ecmf,
                    " name ", name, " units ", units, " level ", level,
                    " type of level ", type_of_level, " short name ", short_name)

            println(" x ",x, " y ", y, " num points ", number_points, " ", lat0 , " ", lon0, " ", lat1, " ", lon1)

            Nearest(message) do near

                lons, lats, values, distances = find(near, message, 8, 50)
                println(lons, lats, values, distances)

            end

            for key in keys(message)
                if key ∉ keylist
                    try
                       push!(keylist, key)
                    catch
                        print(key)
                    end
                end
            end

            println(length(data_array), " took ", time() - band_time)
        end
    end
    println("Overall time: ", time() - start_time)
end

println(Base.summarysize(data_dict)/1000000000, " Gigabytes ")
print(keylist)


text:
#=
Sampling functions to generate random galaxy consistent with provided NFW parameters
=#

#=
Questions from while coding
1) Why more complicated PhiR (genphi()) function in previous functions
2) Why x^2*vt in fprob()
3) Why introduce dummy normalization? in f()
=#

using Optim, GSL

#=
Define gravitationla potential function,
results are in units of PhiS
=#
function GravPotential(x::Real)
    return 1 - (log(1 + x) / x)
end

#=
Calculate escape velocity for a radius scaled by the
characteristic radius, i.e. x = r / rs
=#
function EscapeVelocity(x, rlim, rs, PhiS)
    xlim = rlim / rs
    return sqrt(2 * PhiS * (GravPotential(xlim) - GravPotential(x)))
end

#=
Calculate the  Navarro-Frenk-White density based on x, the radius
scaled by the characteristic radius, vr the radial component of velocity,
vt the tangental component of velocity, and the NFW parameters.  The density
is not normalized, the normalization constant is set to 1000 for conveniance.
=#
function NFWDensity{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1})
    if x < 0
        return 0
    end
    #=
    a = parameters[1]
    d = parameters[2]
    e = parameters[3]
    Ec = parameters[4]
    rlim = parameters[5]
    b = parameters[6]
    q = parameters[7]
    Jb = parameters[8]
    vmax = parameters[9]
    rmax = parameters[10]
    =#

    a, d, e, Ec, rlim, b, q, Jb, vmax, rmax = parameters

    rs = rmax / 2.16
    PhiS = (vmax / 0.465)^2
    PhiR = GravPotential(x) * PhiS

    J = abs(x * rs * vt) #vt = v * sin(theta)
    E = (vt * vt + vr * vr) / 2 + PhiR

    Ec = Ec * PhiS
    xlim = rlim / rs
    Philim = PhiS * GravPotential(xlim)
    Jb = Jb * rs * sqrt(PhiS)

    if b <= 0
        gJ = 1 / (1 + (J / Jb)^(-b))
    else
        gJ = 1 + (J / Jb)^b
    end

    N = 1000
    if E < Philim && E >= 0
        hE = N * E^a * (E^q + Ec^q)^(d/q) * (Philim - E)^e
    else
        hE = 0
    end
    return hE * gJ
end

#=
Scale NFWDensity function
=#
function NFWProb{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1})
    return x * x * vt * NFWDensity(x, vr, vt, parameters)
end

#=
Randomly a generate with samplesize = # stars and NFW parameters parameters.
Basic rejection sampler with a factor of 1.1 is used.  Results are scaled radius (x),
radial velocity (vr), and tangential velocity (vt).
=#
function GenGalaxy{G<:Real}(parameters::Array{G,1}, samplesize::Integer)
    rlim = parameters[5]
    vmax = parameters[9]
    rmax = parameters[10]
    rs = rmax / 2.16
    PhiS = (vmax / 0.465)^2

    function NFWProbOpt(x)
        return -NFWProb(x[1], x[2], x[3], parameters)
    end

    minProb = optimize(NFWProbOpt, [0.1, 1.0, 1.0], method = :nelder_mead)
    fmax = max(0,  -1.1 * minProb.f_minimum)
    x0 = 10.0^-8
    vmax0 = EscapeVelocity(x0, rlim, rs, PhiS)
    #print("fmax: ", fmax, " vmax: ", vmax0, "\n")

    sampledValues = Array{Array{Float64,1}}(samplesize)
    tested = 0
    accepted = 0
    while accepted < samplesize
        x = (rlim  / rs) * rand()
        vr = vmax0 * rand()
        vt = vmax0 * rand()
        u = rand()
        sampleDensity = NFWProb(x, vr, vt, parameters)
        tested += 1
        if sampleDensity >= (u * fmax)
            accepted += 1
            sampledValues[accepted] = [x*rs , vr, vt]
        end
    end
    #print("acceptance rate:", accepted / tested, "\n")
    return sampledValues
end

#=
Transform the results of GenGalaxy to Euclidean coordinates.
=#
function converttoeuclidean{G<:Real}(coordinate::Array{G,1})
    radius = coordinate[1]
    theta = acos(1.0 - 2.0 * rand())
    phi = 2.0 * pi * rand()

    velocity = sqrt(coordinate[2]^2 + coordinate[3]^2)
    velSign = sign(rand() - 0.5)
    velPhi = 2.0 * pi * rand()

    x = radius * sin(theta) * cos(phi)
    y = radius * sin(theta) * sin(phi)
    z = radius * cos(theta)

    vz2 = velSign * coordinate[2]
    vx2 = coordinate[3] * cos(velPhi)
    vy2 = coordinate[3] * sin(velPhi)

    vx = cos(theta) * cos(phi) * vx2 - sin(phi) * vy2 + sin(theta) * cos(phi) * vz2
    vy = cos(theta) * sin(phi) * vx2 + cos(phi) * vy2 + sin(theta) * sin(phi) * vz2
    vz = -sin(theta) * vx2 + cos(theta) * vz2
    return x, y, z, vx, vy, vz
end

#=
Simplified posterior functions
=#
function NFWhEDensity{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1})
    a = parameters[1]
    d = parameters[2]
    e = parameters[3]
    Ec = parameters[4]
    rlim = parameters[5]
    b = parameters[6]
    q = parameters[7]
    Jb = parameters[8]
    vmax = parameters[9]
    rmax = parameters[10]

    rs = rmax / 2.16
    PhiS = (vmax / 0.465)^2
    PhiR = GravPotential(x) * PhiS

    J = abs(x * rs * vt) #vt = v * sin(theta)
    E = (vt * vt + vr * vr) / 2 + PhiR

    Ec = Ec * PhiS
    xlim = rlim / rs
    Philim = PhiS * GravPotential(xlim)
    Jb = Jb * rs * sqrt(PhiS)

    N = 1000
    if E < Philim && E >= 0
        hE = N * E^a * (E^q + Ec^q)^(d/q) * (Philim - E)^e
    else
        hE = 0
    end
    return hE
end

function NFWhEProb{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1})
    return x * x * vt * NFWhEDensity(x, vr, vt, parameters)
end

function GenGalaxySimplified{G<:Real}(parameters::Array{G,1}, samplesize::Integer)
    rlim = parameters[5]
    vmax = parameters[9]
    rmax = parameters[10]
    rs = rmax / 2.16
    PhiS = (vmax / 0.465)^2

    function NFWhEProbOpt(x)
        return -NFWhEProb(x[1], x[2], x[3], parameters)
    end

    minProb = optimize(NFWhEProbOpt, [0.1, 1.0, 1.0], method = :nelder_mead)
    fmax = max(0,  -1.1 * minProb.f_minimum)

    sampledValues = Array{Array{Float64,1}}(samplesize)
    tested = 0
    accepted = 0
    while accepted < samplesize
        x = (rlim  / rs) * rand()
        x0 = 10.0^-8
        vmax0 = EscapeVelocity(x0, rlim, rs, PhiS)
        vr = vmax0 * rand()
        vt = vmax0 * rand()
        u = rand()
        sampleDensity = NFWhEProb(x, vr, vt, parameters)
        tested += 1
        if sampleDensity >= (u * fmax)
            accepted += 1
            sampledValues[accepted] = [x*rs , vr, vt]
        end
    end
    #print("acceptance rate:", accepted / tested, "\n")
    return sampledValues
end

function convertgalaxy(galaxy::Array{Array{Float64, 1}, 1})
    n = length(galaxy)
    x = Array{Float64}(n)
    y = Array{Float64}(n)
    z = Array{Float64}(n)
    vx = Array{Float64}(n)
    vy = Array{Float64}(n)
    vz = Array{Float64}(n)
    for ii in 1:n
        x[ii], y[ii], z[ii], vx[ii], vy[ii], vz[ii] = converttoeuclidean(galaxy[ii])
    end
    return x, y, z, vx, vy, vz
end

function radiusvelocity(galaxy::Array{Array{Float64, 1}, 1})
    n = length(galaxy)
    r = Array{Float64}(n)
    v = Array{Float64}(n)
    for ii in 1:n
        r[ii] = galaxy[ii][1]
        v[ii] = sqrt(galaxy[ii][2]^2 + galaxy[ii][3]^2)
    end
    return r, v
end

function radiusvelocity(x::Array{Float64}, y::Array{Float64}, z::Array{Float64}, vx::Array{Float64}, vy::Array{Float64}, vz::Array{Float64})
    n = length(x)
    r = Array{Float64}(n)
    v = Array{Float64}(n)
    for ii in 1:n
        r[ii] = sqrt(x[ii]^2 + y[ii]^2 + z[ii]^2)
        v[ii] = sqrt(vx[ii]^2 + vy[ii]^2 + vz[ii]^2)
    end
    return r, v
end

function radius(x::Array{Float64}, y::Array{Float64})
    n = length(x)
    r = Array{Float64}(n)
    for ii in 1:n
        r[ii] = sqrt(x[ii]^2 + y[ii]^2)
    end
    return r
end

function rv2(r::Array{Float64, 1}, v::Array{Float64, 1})
    n = length(r)
    quantity = Array{Float64}(n)
    for ii in 1:n
        quantity[ii] = r[ii] * v[ii]^2
    end
    return quantity
end

#new functions

function genphi0_integral(x::Float64, alpha::Float64, beta::Float64, gamma::Float64, rhos::Float64, rs::Float64)
    PhiS = rhos * rs^2
    function integrand1(xi)
        return (xi ^ -gamma) * (1 + xi ^ alpha) ^ ((gamma - beta) / alpha) * xi * xi
    end
    function integrand2(xi)
        return (xi ^ -gamma) * (1 + xi ^ alpha) ^ ((gamma - beta) / alpha) * xi
    end
    phi1 = quadgk(integrand1, 0.0, x)[1]
    phi2 = quadgk(integrand2, x, Inf)[1]

    return PhiS * (1 - (phi1 / x + phi2))
end

function hypergeom_2F1_expanded(a::Float64, b::Float64, c::Float64, z::Float64)
    #if z >= -1.0 && z < 1.0
    #    return hypergeom([a, b], c, z)
    #end
    if z >= -1.0 && z < 0.5
        return hypergeom([a, b], c, z)
    end
    if z >= 0.5 && z < 1.0
        coef1 = gamma(c) * gamma(c - a - b) / (gamma(c - a) * gamma(c - b))
        hypergeom1 = hypergeom_2F1_expanded(a, b, a + b - c + 1, 1 - z)
        coef2 = (1 - z) ^ (c - a - b) * gamma(c) * gamma(a + b - c) / (gamma(a) * gamma(b))
        hypergeom2 = hypergeom_2F1_expanded(c - a, c - b, c - a - b + 1, 1 - z)
        return coef1 * hypergeom1 + coef2 * hypergeom2
    end
    if z >= 1.0
        coef1 = z ^ (-1.0 * a) * gamma(c) * gamma(c - a - b) / (gamma(c - a) * gamma(c - b))
        hypergeom1 = hypergeom_2F1_expanded(a, a - c + 1, a + b - c + 1, 1 - 1 / z)
        coef2 = z ^ (a - c) * (1 - z) ^ (c - a - b) * gamma(c) * gamma(a + b - c) / (gamma(a) * gamma(b))
        hypergeom2 = hypergeom_2F1_expanded(c - a, 1 - a, c - a - b + 1, 1 - 1 / z)
        return coef1 * hypergeom1 + coef2 * hypergeom2
    end
    if z < -1.0
        return (1 - z) ^ (-1.0 * a) * hypergeom_2F1_expanded(a, c - b, c, z / (z - 1))
    end
end

function genphi0(x::Float64, alpha::Float64, beta::Float64, gamma::Float64, rhos::Float64, rs::Float64)
    PhiS = rhos * rs^2
    x0 = 10.0 ^ -13
    p1a = 0
    p1b = hypergeom_2F1_expanded((3.0 - gamma) / alpha, (beta - gamma) / alpha, (3.0 + alpha - gamma) / alpha, -x ^ alpha)
    I1 = (x0 ^ (3.0 - gamma) * p1a - x ^ (3.0 - gamma) * p1b) / (x * (gamma - 3.0))
    p2 = hypergeom_2F1_expanded((beta - 2.0) / alpha, (beta - gamma) / alpha, (alpha + beta - 2.0) / alpha, -x ^ (-alpha))
    I2 = x ^ (2.0 - beta) * p2 / (beta - 2.0)
    return PhiS * (1 - (I1 + I2))
end

function genphi(x::Float64, alpha::Float64, beta::Float64, gamma::Float64, rhos::Float64, rs::Float64)
    x0 = 10.0 ^ -8
    return genphi0(x, alpha, beta, gamma, rhos, rs) -  genphi0(x0, alpha, beta, gamma, rhos, rs)
end

function genphi(x::Float64, alpha::Float64, beta::Float64, gamma::Float64, rhos::Float64, rs::Float64, Phi0::Float64)
    return genphi0(x, alpha, beta, gamma, rhos, rs) -  Phi0
end

#new vesc function line 385
function escape_velocity(x, rlim, rhos, rs, alpha, beta, gamma)
    xlim = rlim / rs
    PhiR = genphi(x, alpha, beta, gamma, rhos, rs)
    Philim = genphi(xlim, alpha, beta, gamma, rhos, rs)
    return sqrt(2 * (Philim - PhiR))
end

function escape_velocity(x, rlim, rhos, rs, alpha, beta, gamma, Philim)
    PhiR = genphi(x, alpha, beta, gamma, rhos, rs)
    return sqrt(2 * (Philim - PhiR))
end

function escape_velocity(x, rlim, rhos, rs, alpha, beta, gamma, Philim, Phi0)
    PhiR = genphi(x, alpha, beta, gamma, rhos, rs, Phi0)
    return sqrt(2 * (Philim - PhiR))
end

#new f function
function SFW_density{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1})
    if x < 0
        return 0
    end
    a, d, e, Ec, rlim, b, q, Jb, rhos, rs, alpha, beta, gamma = parameters

    PhiS = rhos * rs ^ 2
    PhiR = genphi(x, alpha, beta, gamma, rhos, rs)

    J = abs(x * rs * vt) #vt = v * sin(theta)
    E = (vt * vt + vr * vr) / 2.0 + PhiR

    Ec = Ec * PhiS
    xlim = rlim / rs
    Philim = genphi(xlim, alpha, beta, gamma, rhos, rs)
    Jb = Jb * rs * sqrt(PhiS)

    if b <= 0
        gJ = 1.0 / (1.0 + (J / Jb)^(-b))
    else
        gJ = 1.0 + (J / Jb)^b
    end

    N = 1000.0
    if E < Philim && E >= 0
        hE = N * E^a * (E^q + Ec^q)^(d/q) * (Philim - E)^e
    else
        hE = 0
    end
    return hE * gJ
end

function SFW_density{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1}, Philim::G)
    if x < 0
        return 0
    end
    a, d, e, Ec, rlim, b, q, Jb, rhos, rs, alpha, beta, gamma = parameters

    PhiS = rhos * rs ^ 2
    PhiR = genphi(x, alpha, beta, gamma, rhos, rs)

    J = abs(x * rs * vt) #vt = v * sin(theta)
    E = (vt * vt + vr * vr) / 2.0 + PhiR

    Ec = Ec * PhiS
    Jb = Jb * rs * sqrt(PhiS)

    if b <= 0
        gJ = 1.0 / (1.0 + (J / Jb)^(-b))
    else
        gJ = 1.0 + (J / Jb)^b
    end

    N = 1000.0
    if E < Philim && E >= 0
        hE = N * E^a * (E^q + Ec^q)^(d/q) * (Philim - E)^e
    else
        hE = 0
    end
    return hE * gJ
end

function SFW_density{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1}, Philim::G, Phi0::G)
    if x < 0
        return 0
    end
    a, d, e, Ec, rlim, b, q, Jb, rhos, rs, alpha, beta, gamma = parameters

    PhiS = rhos * rs ^ 2
    PhiR = genphi(x, alpha, beta, gamma, rhos, rs, Phi0)

    J = abs(x * rs * vt) #vt = v * sin(theta)
    E = (vt * vt + vr * vr) / 2.0 + PhiR

    Ec = Ec * PhiS
    Jb = Jb * rs * sqrt(PhiS)

    if b <= 0
        gJ = 1.0 / (1.0 + (J / Jb)^(-b))
    else
        gJ = 1.0 + (J / Jb)^b
    end

    N = 1000.0
    if E < Philim && E >= 0
        hE = N * E^a * (E^q + Ec^q)^(d/q) * (Philim - E)^e
    else
        hE = 0
    end
    return hE * gJ
end

#new fprob
function SFW_prob{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1})
    return x * x * vt * SFW_density(x, vr, vt, parameters)
end

function SFW_prob{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1}, Philim)
    return x * x * vt * SFW_density(x, vr, vt, parameters, Philim)
end

function SFW_prob{G<:Real}(x::G, vr::G, vt::G, parameters::Array{G,1}, Philim, Phi0)
    return x * x * vt * SFW_density(x, vr, vt, parameters, Philim, Phi0)
end

#new sample function
function general_gen_galaxy{G<:Real}(parameters::Array{G,1}, samplesize::Integer)
    a, d, e, Ec, rlim, b, q, Jb, rhos, rs, alpha, beta, gamma = parameters
    xlim = rlim / rs
    x0 = 10.0 ^ -8
    Phi0 = genphi0(x0, alpha, beta, gamma, rhos, rs)
    Philim = genphi(xlim, alpha, beta, gamma, rhos, rs, Phi0)
    function SFW_prob_opt(x)
        return -SFW_prob(x[1], x[2], x[3], parameters, Philim, Phi0)
    end

    minProb = optimize(SFW_prob_opt, [0.1, 1.0, 1.0], method = :nelder_mead)
    #minProb = optimize(SFW_prob_opt, [0.1, 1.0, 1.0], method = Optim.NelderMead())
    fmax = max(0,  -1.1 * minProb.f_minimum)
    x0 = 10.0^-8
    vmax0 = escape_velocity(x0, rlim, rhos, rs, alpha, beta, gamma, Philim, Phi0)
    #print("fmax: ", fmax, " vmax: ", vmax0, "\n")

    sampledValues = Array{Array{Float64,1}}(samplesize)
    tested = 0
    accepted = 0
    while accepted < samplesize
        x = (rlim  / rs) * rand()
        vr = vmax0 * rand()
        vt = vmax0 * rand()
        u = rand()
        sampleDensity = SFW_prob(x, vr, vt, parameters, Philim, Phi0)
        tested += 1
        if sampleDensity >= (u * fmax)
            accepted += 1
            sampledValues[accepted] = [x * rs , vr, vt]
        end
    end
    #print("acceptance rate:", accepted / tested, "\n")
    return sampledValues
end

#=
param=[2.0, -5.3, 2.5, 0.16, 1.5, -9.0, 6.9, 0.086, 4229.2, .694444, 1.01, 3.01, 1.01]
@time general_gen_galaxy(param, 1000); #7.5 seconds before Philim
@time general_gen_galaxy(param, 1000); #3.6 seconds before Phi0
@time general_gen_galaxy(param, 1000); #2.1 seconds with Philim and Phi0

#hypergeom([1.0, 2.0], 2.0, -0.98) #julia syntax
#ss.hyp2f1(1.0, 2.0, 2.0, -1.0) #python syntax
hypergeom()

@time genphi0(x, alpha, beta, gamma, rhos, rs)
@time genphi0_2(1.1, alpha, beta, gamma, rhos, rs)

@everywhere function f(n::Int64)
    x, vt, vr = general_gen_galaxy(param, 100)
    return mean(x)
end
=#


text:
module PaperData

import PowerModelsDistribution, PowerModelsDistributionStateEstimation

const _PMD = PowerModelsDistribution
const _PMDSE = PowerModelsDistributionStateEstimation

include("non_gaussian_se.jl")

export get_EU_LV_feeder_data, get_feeder_data

end

max_stars_repo_name: MartaVanin/PaperData

max_issues_repo_name: MartaVanin/PaperData

max_forks_repo_name: MartaVanin/PaperData


text:
import BioMetaheuristics.TestFunctions: Benchmark, evaluate

# Some useful constants
const ln2 = log(2.0)
const sqrtpi = √π

struct SimulatedAnnealing <: TrajectoryBase end

struct GeneralSimulatedAnnealing <: TrajectoryBase end

# To enable dispatch based on the type
function optimize(f, range, dim, iters, rng, ::SimulatedAnnealing; kwargs...)
    return SimulatedAnnealing(f, range..., dim, rng; low_temp=iters, kwargs...)
end

function optimize(f, range, dim, iters, rng, ::GeneralSimulatedAnnealing; kwargs...)
    return GeneralSimulatedAnnealing(f, range..., dim, rng; low_temp=iters, kwargs...)
end

"""
    SimulatedAnnealing(f::Function, a, b, dim::Integer, rng;
        low_temp=5000, t0=500.0
    ) -> OptimizationResults
    SimulatedAnnealing(f::Benchmark, a, b, dim::Integer, rng;
        low_temp=5000, t0=500.0
    ) -> OptimizationResults

Implementation of the *classical* version of simulated annealing.

This implementation uses a logarithmic cooling schedule and searches possible candidate solutions by sampling from an approximated Boltzmann distribution, drawn as a normal distribution.

Returns an `OptimizationResults` type with information relevant to the run executed, see [`OptimizationResults`](@ref).

# Arguments

- `f`: any user defined `Function` that can take `AbstractArray`'s.
- `a`: **lower** bound for the solution search space.
- `b`: **upper** bound for the solution search space.
- `dim`: dimension of the optimization problem.
- `rng`: an object of type `AbstractRNG`.

# Keyword arguments

_It is recommended to use the default values provided._

- `t0`: initial value for the *temperature* that is used. The default is an okay value, but should be changed depending on the optimization problem.
- `low_temp`: total number of iterations, short for *lowering temperature steps*. This also corresponds to the famous *Monte Carlo steps*, which are the total number of steps until the algorithm finishes.

# Examples

```julia
using BioMetaheuristics
using Random

rng = MersenneTwister()

# Define the 2D Rosenbrock function
rosenbrock2d(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

# Implement Simulated Annealing for a 2-dimensional Rosenbrock function, with
# 5000 iterations.
val = SimulatedAnnealing(rosenbrock2d, -5.0, 5.0, 2, rng; low_temp=5000)
```
"""
function SimulatedAnnealing(
    f::Function,
    a,
    b,
    dim::Integer,
    rng;
    t0=500.0,
    low_temp=5000
)
    # Enforce type stability
    a, b = promote(a, b)

    # Create the solution with random initial values within the
    # search space bounds
    x_solution = a .+ rand(rng, typeof(a), dim) * (b - a)

    # We create a temporal array for copying values
    xtmp = similar(x_solution)

    t = t0

    for j = 1:low_temp
        # Determine the temperature for current iteration
        x = _temperature!(t)

        # Create a neighbor solution
        _classical_visit(x_solution, xtmp, sqrt(x), rng)

        # Clip the resulting solution to the bounds
        _clip_trajectory!(x_solution, a, b)

        # Employ the Metropolis-Hastings algorithm
        _annealing!(f, x, x_solution, xtmp, rng)

        # Update the temperature with the classical cooling schedule
        t *= x * t0
    end

    res = OptimizationResults(x_solution, f(x_solution), "SimulatedAnnealing", low_temp)

    return res
end  # function SimulatedAnnealing

function SimulatedAnnealing(
    f::Benchmark,
    a,
    b,
    dim::Integer,
    rng;
    t0=500.0,
    low_temp=5000
)
    return SimulatedAnnealing(x -> evaluate(f, x), a, b, dim, rng;
        t0=t0, low_temp=low_temp)
end

@inline function _temperature!(x)
    # Avoid evaluating the logarithm at zero
    @assert x > -1.0

    return ln2 / log(1.0 + x)
end  # function temperature

function _classical_visit(x, xtmp, σ, rng)
    for i in eachindex(x)
        @inbounds xtmp[i] = x[i] + randn(rng, Float64) * σ
    end

    return nothing
end  # function _classical_visit!

function _annealing!(f::Function, t, x, xtmp, rng)
    # Evaluate the energies
    new_energy = f(xtmp)
    old_energy = f(x)

    # If the new energy is better, update it
    if new_energy <= old_energy
        copyto!(x, xtmp)
    # If not, we accept it with the Metropolis-Hastings algorithm
    else
        Δ = old_energy - new_energy
        # Compute acceptance probability
        acceptance = exp(Δ / t)
        # Metropolis-Hastings criterion
        if rand(rng) <= acceptance
            copyto!(x, xtmp)
        end
    end

    return nothing
end  # function _annealing

@doc raw"""
    GeneralSimulatedAnnealing(f::Function, a, b, dim::Integer, rng;
        low_temp=5_000, t0=500.0, qv=2.7, qa=-5.0
    ) -> OptimizationResults
    GeneralSimulatedAnnealing(f::Benchmark, a, b, dim::Integer, rng;
        low_temp=20_000, t0=500.0, qv=2.7, qa=-5.0
    ) -> OptimizationResults

Implementation of the *generalized* version of simulated annealing.

This implementation uses all the theory from Tsallis & Stariolo for the cooling schedule and the neighbor solution search.
See [`GeneralSimulatedAnnealing`](@ref generalized-sm) for the implementation details.

Returns an `OptimizationResults` type with information relevant to the run executed, see [`OptimizationResults`](@ref).

# Arguments

- `f`: any user defined `Function` that can take `AbstractArray`'s.
- `a`: **lower** bound for the solution search space.
- `b`: **upper** bound for the solution search space.
- `dim`: dimension of the optimization problem.
- `rng`: an object of type `AbstractRNG`.

# Keyword arguments
_It is recommended to use the default values provided._

- `t0`: initial value for the *temperature* that is used. The default is an okay value, but should be changed depending on the optimization problem.
- `low_temp`: total number of iterations, short for *lowering temperature steps*. This also corresponds to the famous *Monte Carlo steps*, which are the total number of steps until the algorithm finishes.
- `qv`: This is known as the *Tsallis parameter*; particularly this parameter controls the cooling schedule convergence and neighbor search. Positive values in the interval ``[1,5/3)`` are best because for values larger than 5/3 the neighbor search diverges.
- `qa`: Another *Tsallis parameter*; this particular parameter controls convergence for the Metropolis-Hastings algorithm and the acceptance probability involved. The more negative the value is, the better, but Tsallis & Stariolo report that the default value is best.

# Examples

```julia
using BioMetaheuristics
using Random

rng = MersenneTwister()
# Define the 2D Rosenbrock function
rosenbrock2d(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

# Implement Simulated Annealing for a 2-dimensional Rosenbrock function, with
# 15000 iterations and default values.
val = GeneralSimulatedAnnealing(rosenbrock2d, -5.0, 5.0, 2, rng; low_temp=15_000)
```
"""
function GeneralSimulatedAnnealing(
    f::Function,
    a,
    b,
    dim::Integer,
    rng;
    t0=500.0,
    low_temp=20_000,
    qv=2.7,
    qa=-5.0
)
    # Enforce type stability
    a, b = promote(a, b)

    # Create the solution with random initial values within the
    # search space bounds
    x_solution = a .+ rand(rng, typeof(a), dim) * (b - a)

    # We create a temporal array for copying values
    xtmp = similar(x_solution)

    t = t0

    for j = 1:low_temp
        # Determine the temperature for current iteration
        x = _general_temperature!(t, qv)

        # Create a neighbor solution
        _general_visit(x_solution, xtmp, x, qv, rng)

        # Clip the resulting solution to the bounds
        _clip_trajectory!(x_solution, a, b)

        # Employ the Metropolis-Hastings algorithm
        _general_annealing!(f, x, x_solution, xtmp, qa, rng)

        # Update the temperature with the generalized cooling schedule
        t *= x * t0
    end

    res = OptimizationResults(x_solution, f(x_solution), "GeneralizedSimulatedAnnealing", low_temp)

    return res
end  # function GeneralSimulatedAnnealing

function GeneralSimulatedAnnealing(
    f::Benchmark,
    a,
    b,
    dim::Integer,
    rng;
    t0=500.0,
    low_temp=20_000,
    qv=2.7,
    qa=-5.0
)
    return GeneralSimulatedAnnealing(x -> evaluate(f, x), a, b, dim, rng;
        t0=t0, low_temp=low_temp, qv=qv, qa=qa)
end

"""
    This is the algorithm from Tsallis & Stariolo to sample the distribution
shown in their paper, by approximating their probability distribution as a
Lévy probability distribution.

See the appendix on the paper:
Tsallis, C. and Stariolo, D. A. (1996) ‘Generalized simulated annealing’, Physica A: Statistical Mechanics and its Applications, 233(1), pp. 395–406. doi: 10.1016/S0378-4371(96)00271-3.
"""
function _general_visit(
    x::AbstractArray{T},
    xtmp::AbstractArray{T},
    τ::T,
    q::T,
    rng
) where {T <: Real}
    factor_1 = exp(log(τ) / (q - 1.0))
    factor_2 = exp((4.0 - q) * log(q - 1.0))
    factor_3 = exp((2.0 - q) * ln2 / (q - 1.0))
    factor_4 = (sqrtpi * factor_1 * factor_2) / (factor_3 * (3.0 - q))
    factor_5 = (1.0 / (q - 1.0)) - 0.5
    factor_6 = π * (1.0 - factor_5) / sin(π * (1.0 - factor_5)) / exp(_gammaln(2.0 - factor_5))
    sigmax = exp(-(q - 1.0) * log(factor_6 / factor_4)) / (3.0 - q)

    @inbounds for i in eachindex(x)
        xx = sigmax .* randn(rng, eltype(x))
        y = randn(rng, eltype(x))
        den = exp((q - 1.0) * log(abs(y)) / (3.0 - q))
        xtmp[i] = x[i] + (xx / den)
    end

    return nothing
end

function _general_annealing!(
    f::Function,
    t::T,
    x::AbstractArray{T},
    xtmp::AbstractArray{T},
    q::T,
    rng
) where {T <: Real}

    # Evaluate the energies
    new_energy = f(xtmp)
    old_energy = f(x)

    # If the new energy is better, update it
    if new_energy < old_energy
        copyto!(x, xtmp)
    # If not, we accept it with the Metropolis-Hastings algorithm
    else
        Δ = new_energy - old_energy
        if q < 0
            acceptance = ((q - 1.0) * Δ) / t
            acceptance += 1.0
            if acceptance >= 0
                if rand(rng) <= 1.0 / acceptance
                    copyto!(x, xtmp)
                end
            end
        else
            power = 1.0 / (q - 1.0)
            # Compute acceptance probability
            acceptance = ((q - 1.0) * Δ) / t
            acceptance += 1.0
            acceptance = acceptance^power
            acceptance = 1.0 / acceptance
            # Metropolis-Hastings criterion
            if rand(rng) <= acceptance
                copyto!(x, xtmp)
            end
        end
    end

    return nothing
end  # function _general_annealing!

"""
    A generalized cooling schedule that should work for every possible type
of Simulated Annealing implementation. When `q` = 1, this reduces to the
classical version.
"""
function _general_temperature!(x, q)
    # Avoid evaluating the logarithm at zero
    @assert x > -1.0

    factor_1 = 2^(q - 1.0) - 1.0
    factor_2 = (1.0 + x)^(q - 1.0) - 1.0

    return factor_1 / factor_2
end


text:
using LightXML

# Use the files provided
filename = ARGS[1]
outputFile = ascii(ARGS[2] * ".xml")

#Read from file and remove trailing \n\r and spaces
function parseShit(filename::String)
str = strip(readall(open(filename)))
#parse file, does not check if the program is correct or not!
#Just linewise check.
i = start(str)
exprs = {}
while !done(str,i)
	# try
	if str[i] == '#'
		println("this is a comment, $i")
	end
	ex, i = parse(str,i, raise=true)
	push!(exprs, ex)
# catch y
# 	println("Error!")
# 	if isa(y, Exception)
# 		# check if this line is a comment
# 		if contains(str, '#')
# 			println("This is a comment!")
# 		end
# 		println(y)
# 	end
# end
end
exprs
end
function lowerFirstChar(x::String)
	string(lowercase(x[1])) * x[2:end]
end

#create XML
xdoc = XMLDocument()
xroot = create_root(xdoc, "juliaTree")

function exprToXML(x::Any, node)
	set_attribute(node, lowerFirstChar(string(typeof(x))), string(x))
end

function exprToXML(x::String, node)
	add_text(node, x)
end

function exprToXML(x::Expr, node)
	xNode = new_child(node, "expr")
	set_attribute(xNode, "head", string(x.head))
	for elem in x.args
		exprToXML(elem, new_child(xNode, "args"))
	end
end

for expr in exprs
	exprToXML(expr, xroot)
end

save_file(xdoc, outputFile)

max_stars_repo_name: dadarakt/jamon

max_issues_repo_name: dadarakt/jamon

max_forks_repo_name: dadarakt/jamon


text:
module SomervilleCouncilParser

using HTTP
using EzXML
using HttpCommon
using Logging
using Dates
using Legolas
using DataFrames
using ProgressMeter

include("meeting.jl")
export request_meetings, Meeting

include("agenda.jl")
export get_agenda_items, filter_agenda, search_agendas_for_content,
       display_items_by_meeting, Agenda

end # module


text:
using Knet, JLD

include("retinaface.jl")
include("../backbones/resnet.jl")
include("../backbones/fpn.jl")
include("../backbones/context_module.jl")
include("../backbones/ssh.jl")

function load_pth_model(model::RetinaFace, path; dtype=Array{Float32}, load_heads=true)
    c = JLD.jldopen(path, "r") do file
        read(file, "data")
    end
    data = Dict()
    for k in keys(c)
        if !isempty(k) data[k] = c[k] end
    end
    return set_pth_data(model, data; dtype=dtype, load_heads=load_heads)
end

function set_pth_data(model::RetinaFace, data; dtype=Array{Float32}, load_heads=true)
    conv_w, conv_b, bn_mom, bn_b, bn_mult = get_pth_resnet_params(data, dtype=dtype)
    model.backbone = load_mat_weights(model.backbone, nothing; pre_weights=[conv_w, conv_b, nothing, nothing, bn_mom, bn_b, bn_mult])
    model = set_pth_fpn_data(model, data, dtype=dtype)
    model = set_pth_context_module_data(model, data, dtype=dtype)
    if load_heads
        model = set_pth_head_getters_data(model, data, dtype=dtype)
    end

    for p in params(model)
        p = Param(convert(dtype, value(p)))
    end

    return model
end

function get_pth_resnet_params(data; dtype=Array{Float32})
    conv_w = []; conv_b = convert(dtype, zeros(1, 1, size(data["module.body.conv1.weight"], 1), 1));
    bn_mom = []; bn_b = []; bn_mult = [];

    push!(conv_w, convert(dtype, reverse(reverse(permutedims(data["module.body.conv1.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    push!(bn_mom, convert(dtype, reshape(data["module.body.bn1.running_mean"], (1, 1, size(data["module.body.bn1.running_mean"], 1), 1))))
    push!(bn_mom, convert(dtype, reshape(data["module.body.bn1.running_var"], (1, 1, size(data["module.body.bn1.running_var"], 1), 1))))
    push!(bn_mult, convert(dtype, vec(data["module.body.bn1.weight"])))
    push!(bn_b, convert(dtype, vec(data["module.body.bn1.bias"])))

    stages = ["layer1", "layer2", "layer3", "layer4"]

    for stage in stages
        ws, moms, bs, mults = get_pth_block_data(data, stage, dtype=dtype)
        conv_w = vcat(conv_w, ws);
        bn_mom = vcat(bn_mom, moms);  bn_mult = vcat(bn_mult, mults);  bn_b = vcat(bn_b, bs);
    end
    return conv_w, conv_b, bn_mom, bn_b, bn_mult
end

function set_pth_fpn_data(model::RetinaFace, data; dtype=Array{Float32})
    # ResNet-to-FPN connections
    model.fpn.o5.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data["module.fpn.output3.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.fpn.o5.bn.bn_params = Param(convert(dtype, vcat(
            vec(data["module.fpn.output3.1.weight"]),
            vec(data["module.fpn.output3.1.bias"]))))
    model.fpn.o5.bn.bn_moments = bnmoments(
        mean=convert(dtype, data["module.fpn.output3.1.running_mean"]),
        var=convert(dtype, data["module.fpn.output3.1.running_var"]))

    model.fpn.o4.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data["module.fpn.output2.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.fpn.o4.bn.bn_params = Param(convert(dtype, vcat(
            vec(data["module.fpn.output2.1.weight"]),
            vec(data["module.fpn.output2.1.bias"]))))
    model.fpn.o4.bn.bn_moments = bnmoments(
        mean=convert(dtype, data["module.fpn.output2.1.running_mean"]),
        var=convert(dtype, data["module.fpn.output2.1.running_var"]))

    model.fpn.o3.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data["module.fpn.output1.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.fpn.o3.bn.bn_params = Param(convert(dtype, vcat(
            vec(data["module.fpn.output1.1.weight"]),
            vec(data["module.fpn.output1.1.bias"]))))
    model.fpn.o3.bn.bn_moments = bnmoments(
        mean=convert(dtype, data["module.fpn.output1.1.running_mean"]),
        var=convert(dtype, data["module.fpn.output1.1.running_var"]))

    # Merge connections
    model.fpn.merge4.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data["module.fpn.merge2.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.fpn.merge4.bn.bn_params = Param(convert(dtype, vcat(
            vec(data["module.fpn.merge2.1.weight"]),
            vec(data["module.fpn.merge2.1.bias"]))))
    model.fpn.merge4.bn.bn_moments = bnmoments(
        mean=convert(dtype, data["module.fpn.merge2.1.running_mean"]),
        var=convert(dtype, data["module.fpn.merge2.1.running_var"]))

    model.fpn.merge3.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data["module.fpn.merge1.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.fpn.merge3.bn.bn_params = Param(convert(dtype, vcat(
            vec(data["module.fpn.merge1.1.weight"]),
            vec(data["module.fpn.merge1.1.bias"]))))
    model.fpn.merge3.bn.bn_moments = bnmoments(
        mean=convert(dtype, data["module.fpn.merge1.1.running_mean"]),
        var=convert(dtype, data["module.fpn.merge1.1.running_var"]))

    return model
end

function set_pth_context_module_data(model::RetinaFace, data; dtype=Array{Float32})
    model.context_module.ssh_p3 = set_pth_ssh_data(model.context_module.ssh_p3, "module.ssh1", data, dtype=dtype)
    model.context_module.ssh_p4 = set_pth_ssh_data(model.context_module.ssh_p4, "module.ssh2", data, dtype=dtype)
    model.context_module.ssh_p5 = set_pth_ssh_data(model.context_module.ssh_p5, "module.ssh3", data, dtype=dtype)
    return model
end

function set_pth_head_getters_data(model::RetinaFace, data; dtype=Array{Float32})
    model.cls_head2 = set_pth_head_data(model.cls_head2, "module.ClassHead", data, dtype=dtype)
    model.bbox_head2 = set_pth_head_data(model.bbox_head2, "module.BboxHead", data, dtype=dtype)
    model.lm_head2 = set_pth_head_data(model.lm_head2, "module.LandmarkHead", data, dtype=dtype)
    return model
end

function set_pth_head_data(model::HeadGetter, start_txt, data; dtype=Array{Float32})
    model.layers[1].w = Param(convert(dtype, reverse(reverse(permutedims(data[start_txt * ".0.conv1x1.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.layers[1].b = Param(convert(dtype,
        reshape(data[start_txt * ".0.conv1x1.bias"], (1, 1, size(data[start_txt * ".0.conv1x1.bias"], 1), 1))))

    model.layers[2].w = Param(convert(dtype, reverse(reverse(permutedims(data[start_txt * ".1.conv1x1.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.layers[2].b = Param(convert(dtype,
        reshape(data[start_txt * ".1.conv1x1.bias"], (1, 1, size(data[start_txt * ".1.conv1x1.bias"], 1), 1))))

    model.layers[3].w = Param(convert(dtype, reverse(reverse(permutedims(data[start_txt * ".2.conv1x1.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.layers[3].b = Param(convert(dtype,
        reshape(data[start_txt * ".2.conv1x1.bias"], (1, 1, size(data[start_txt * ".2.conv1x1.bias"], 1), 1))))

    if start_txt == "module.ClassHead"
        model.task_len = 2
    elseif start_txt == "module.BboxHead"
        model.task_len = 4
    else
        model.task_len = 10
    end

    return model
end

function set_pth_ssh_data(model::SSH, start_txt, data; dtype=Array{Float32})
    model.conv128.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data[start_txt * ".conv3X3.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.conv128.bn.bn_params = Param(convert(dtype, vcat(
        vec(data[start_txt * ".conv3X3.1.weight"]), vec(data[start_txt * ".conv3X3.1.bias"]))))
    model.conv128.bn.bn_moments = bnmoments(
        mean=convert(dtype, data[start_txt * ".conv3X3.1.running_mean"]),
        var=convert(dtype, data[start_txt * ".conv3X3.1.running_var"]))

    model.conv64_1_1.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data[start_txt * ".conv5X5_1.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.conv64_1_1.bn.bn_params = Param(convert(dtype, vcat(
        vec(data[start_txt * ".conv5X5_1.1.weight"]), vec(data[start_txt * ".conv5X5_1.1.bias"]))))
    model.conv64_1_1.bn.bn_moments = bnmoments(
        mean=convert(dtype, data[start_txt * ".conv5X5_1.1.running_mean"]),
        var=convert(dtype, data[start_txt * ".conv5X5_1.1.running_var"]))

    model.conv64_1_2.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data[start_txt * ".conv5X5_2.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.conv64_1_2.bn.bn_params = Param(convert(dtype, vcat(
        vec(data[start_txt * ".conv5X5_2.1.weight"]), vec(data[start_txt * ".conv5X5_2.1.bias"]))))
    model.conv64_1_2.bn.bn_moments = bnmoments(
        mean=convert(dtype, data[start_txt * ".conv5X5_2.1.running_mean"]),
        var=convert(dtype, data[start_txt * ".conv5X5_2.1.running_var"]))

    model.conv64_2_1.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data[start_txt * ".conv7X7_2.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.conv64_2_1.bn.bn_params = Param(convert(dtype, vcat(
        vec(data[start_txt * ".conv7X7_2.1.weight"]), vec(data[start_txt * ".conv7X7_2.1.bias"]))))
    model.conv64_2_1.bn.bn_moments = bnmoments(
        mean=convert(dtype, data[start_txt * ".conv7X7_2.1.running_mean"]),
        var=convert(dtype, data[start_txt * ".conv7X7_2.1.running_var"]))

    model.conv64_2_2.conv.w = Param(convert(dtype, reverse(reverse(permutedims(data[start_txt * ".conv7x7_3.0.weight"], (4, 3, 2, 1)), dims=1), dims=2)))
    model.conv64_2_2.bn.bn_params = Param(convert(dtype, vcat(
        vec(data[start_txt * ".conv7x7_3.1.weight"]), vec(data[start_txt * ".conv7x7_3.1.bias"]))))
    model.conv64_2_2.bn.bn_moments = bnmoments(
        mean=convert(dtype, data[start_txt * ".conv7x7_3.1.running_mean"]),
        var=convert(dtype, data[start_txt * ".conv7x7_3.1.running_var"]))

    return model
end

function get_pth_block_data(data, stagename; dtype=Array{Float32})
    conv_w = []; bn_mom = []; bn_b = []; bn_mult = [];
    conv_ds = nothing; ds_bn_b = nothing; ds_bn_mult = nothing; ds_bn_moms = [];

    kvals = []
    for k in keys(data)
        if !isempty(k)
            push!(kvals, string(k))
        end
    end

    for k in sort(kvals)
        if typeof(data[k]) == Int64
            continue
        end

        val = data[k]
        # print(k, " --> ", size(val), "\n")

        if startswith(k, "module.body." * stagename) && occursin("conv", k) && endswith(k, "weight")
            # residual weight
            push!(conv_w, convert(dtype, reverse(reverse(permutedims(val, (4, 3, 2, 1)), dims=1), dims=2)))
        elseif startswith(k, "module.body." * stagename) && endswith(k, "downsample.0.weight")
            # downsample weight
            conv_ds = convert(dtype, reverse(reverse(permutedims(val, (4, 3, 2, 1)), dims=1), dims=2))
        else
            # downsample parameters
            val = convert(dtype, val)
            if startswith(k, "module.body." * stagename) && endswith(k, "downsample.1.bias")
                ds_bn_b = vec(val)
            elseif startswith(k, "module.body." * stagename) && endswith(k, "downsample.1.weight")
                ds_bn_mult = vec(val)
            elseif startswith(k, "module.body." * stagename) && endswith(k, "downsample.1.running_mean")
                push!(ds_bn_moms, reshape(val, (1, 1, size(val, 1), 1)))
            elseif startswith(k, "module.body." * stagename) && endswith(k, "downsample.1.running_var")
                push!(ds_bn_moms, reshape(val, (1, 1, size(val, 1), 1)))
            # residual parameters
            elseif startswith(k, "module.body." * stagename) && endswith(k, "bias")
                push!(bn_b, vec(val))
            elseif startswith(k, "module.body." * stagename) && endswith(k, "weight")
                push!(bn_mult, vec(val))
            elseif startswith(k, "module.body." * stagename) && occursin("running_", k)
                push!(bn_mom, reshape(val, (1, 1, size(val, 1), 1)))
            end
        end
    end

    conv_w = [conv_ds, conv_w...]
    bn_mom = [ds_bn_moms..., bn_mom...]
    bn_b = [ds_bn_b, bn_b...]
    bn_mult = [ds_bn_mult, bn_mult...]
    return conv_w, bn_mom, bn_b, bn_mult
end

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
using Distributed
using Random
using Logging

workers = 8

if nprocs() <= workers
    addprocs(workers + 1 - nprocs())
end

@everywhere include("models/Vibron.jl")
@everywhere include("modules/ClassicalDynamics.jl")

@everywhere disable_logging(Logging.Info)

""" Calculates Poincaré sections with Lyapunov exponents for various energies
    Parallel calculation, takes usually days to finish
"""

@everywhere function SolveItem(energy, parameters, dimension; file="Vibron.txt", path="", alreadySolved=[], sectionPlane=2)
    for (E, λ) in alreadySolved
        if isapprox(E, energy) && isapprox(parameters[1], λ)
            println("Skipped $parameters, E=$energy, dim=$dimension (already calculated)")
            return false
        end
    end

    time = @elapsed averageLyapunov, freg, trajectories, lyapunovs = SolveEnergy(energy, parameters, dimension, savePath=path, timeout=7200, showFigures=false, randomize=true)

    chaos = 0
    total = 0
    error = 0
    meanLyapunov = 0
    meanLyapunovChaos = 0

    for x in averageLyapunov
        if x > 0.0
            total += 1
            meanLyapunov += x
        end
        if x == 0.0
            error += 1
        end
        if x > 0.01
            chaos += 1
            meanLyapunovChaos += x
        end
    end

    varλ = 0
    if length(lyapunovs) == 0
        maxλ = 0
        meanλ = 0
    else
        maxλ = maximum(lyapunovs)
        meanλ = mean(lyapunovs)
        if length(lyapunovs) > 1
            varλ = var(lyapunovs)
        end
    end

    result = [energy, parameters[1], total, chaos, error, freg, total > 0 ? meanLyapunov / total : 0, chaos > 0 ? meanLyapunovChaos / chaos : 0, length(lyapunovs), maxλ, meanλ, varλ, myid(), time, trajectories]

    open(path * file, "a") do io
        println(io, result)
    end

    return true
end

function RunMapC(; C=0.2, path="", dimension=101, step=0.1, sectionPlane=3)
    path *= "Vibron_"

    file = "Map_dim=$(dimension)_C=$C.txt"
    alreadySolved = ReadMap(path * file)

    input = shuffle([(energy, [0.5 * (1 - A), -A, C], dimension) for energy in -1.2:step:1, A in 0:step:1])

    println("To be calculated $(length(input)).")
    println("Already calculated $(length(alreadySolved)) points.")

    pmap((args)->SolveItem(args...; file=file, path=path, alreadySolved=alreadySolved, sectionPlane=sectionPlane), input)

    return
end

""" Calculates freg for one given A and C """
function RunAC(; C=0.2, A=0.4, path="", dimension=101, step=0.01, sectionPlane=3)
    path *= "Vibron_"

    file = "Energy_dim=$(dimension)_$([A, C]).txt"
    alreadySolved = ReadMap(path * file)

    input = [(energy, [0.5 * (1 - A), -A, C], dimension) for energy in -1.2:step:1]

    println("To be calculated $(length(input)).")
    println("Already calculated $(length(alreadySolved)) points.")

    pmap((args)->SolveItem(args...; file=file, path=path, alreadySolved=alreadySolved, sectionPlane=sectionPlane), input)

    return
end

function ReadMap(file="")
    Eλ = []

    try
        for line in eachline(file)
            line = replace(line, "[" => "")
            line = replace(line, "]" => "")
            line = replace(line, "," => "")
            elements = split(line)

            append!(Eλ, [(parse(Float64, elements[1]), parse(Float64, elements[2]))])
        end
    catch x
    end

    return Eλ
end

max_stars_repo_name: PavelStransky/ClassicalDynamics

max_issues_repo_name: PavelStransky/ClassicalDynamics

max_forks_repo_name: PavelStransky/ClassicalDynamics


text:
"""
    compute_col_dry_kernel!(
        col_dry,
        p_lev,
        mol_m_dry,
        mol_m_h2o,
        avogadro
        helmert1,
        vmr_h2o
        lat,
        glaycol,
    )

This function computes the column amounts of dry or moist air.
"""
function compute_col_dry_kernel!(
    col_dry::AbstractArray{FT, 2},
    p_lev::AbstractArray{FT, 2},
    mol_m_dry::FT,
    mol_m_h2o::FT,
    avogadro::FT,
    helmert1::FT,
    vmr_h2o::Union{AbstractArray{FT, 2}, Nothing},
    lat::Union{AbstractArray{FT, 1}, Nothing},
    glaycol::Tuple{Int, Int},
) where {FT <: AbstractFloat}
    glay, gcol = glaycol[1], glaycol[2]  # global col & lay ids
    helmert2 = FT(0.02586)  # second constant of Helmert formula
    m2_to_cm2 = FT(100 * 100) # m^2 to cm^2
    if lat isa AbstractArray
        g0 = helmert1 - helmert2 * cos(FT(2) * FT(π) * lat[gcol] / FT(180)) # acceleration due to gravity [m/s^2]
    else
        g0 = helmert1
    end
    Δp = p_lev[glay, gcol] - p_lev[glay + 1, gcol]
    if vmr_h2o isa AbstractArray
        fact = FT(1) / (FT(1) + vmr_h2o[glay, gcol])
        # Get average mass of moist air per mole of moist air
        m_air = (mol_m_dry + mol_m_h2o * vmr_h2o[glay, gcol]) * fact
        # Hydrostatic equation
        col_dry[glay, gcol] = (Δp * avogadro / (m2_to_cm2 * m_air * g0)) * fact # molecules/cm^2
    else
        col_dry[glay, gcol] = (Δp * avogadro / (m2_to_cm2 * mol_m_air * g0)) # molecules/cm^2
    end
end
"""
    compute_interp_fractions(
        lkp::AbstractLookUp{I,FT},
        vmr,
        p_lay,
        t_lay,
        tropo,
        ibnd,
        glaycol,
    ) where {I<:Int,FT<:AbstractFloat}

compute interpolation fractions for binary species parameter, pressure and temperature.
"""
@inline function compute_interp_fractions(
    lkp::AbstractLookUp{I, FT},
    vmr,
    p_lay,
    t_lay,
    tropo,
    ibnd,
    glaycol,
) where {I <: Int, FT <: AbstractFloat}
    jftemp = compute_interp_frac_temp(lkp, t_lay, glaycol...)
    jfpress = compute_interp_frac_press(lkp, p_lay, tropo, glaycol...)
    jfη, col_mix = compute_interp_frac_η(lkp, vmr, tropo, jftemp[1], ibnd, glaycol...)
    return (jftemp, jfpress, jfη, col_mix)
end

"""
    compute_interp_frac_temp(
        lkp::AbstractLookUp{I,FT},
        t_lay,
        glay,
        gcol,
    ) where {I<:Int,FT<:AbstractFloat}

compute interpolation fraction for temperature.
"""
@inline function compute_interp_frac_temp(
    lkp::AbstractLookUp{I, FT},
    t_lay,
    glay,
    gcol,
) where {I <: Int, FT <: AbstractFloat}
    (; Δ_t_ref, n_t_ref, t_ref) = lkp

    @inbounds jtemp = loc_lower(t_lay, Δ_t_ref, n_t_ref, t_ref)
    @inbounds ftemp = (t_lay - t_ref[jtemp]) / Δ_t_ref

    return (jtemp, ftemp)
end

"""
    compute_interp_frac_press(
        lkp::AbstractLookUp,
        p_lay,
        tropo,
        glay,
        gcol,
    )

Compute interpolation fraction for pressure.
"""
@inline function compute_interp_frac_press(lkp::AbstractLookUp, p_lay, tropo, glay, gcol)
    (; Δ_ln_p_ref, p_ref, n_p_ref) = lkp

    @inbounds jpress = Int(min(max(fld(log(p_ref[1]) - log(p_lay), Δ_ln_p_ref) + 1, 1), n_p_ref - 1) + 1)

    @inbounds fpress = (log(p_ref[jpress - 1]) - log(p_lay)) / Δ_ln_p_ref
    jpress = jpress + tropo - 1

    return (jpress, fpress)
end

"""
    compute_interp_frac_η(
        lkp::AbstractLookUp{I,FT},
        vmr,
        tropo,
        jtemp,
        ibnd,
        glay,
        gcol,
    ) where {FT<:AbstractFloat,I<:Int}

Compute interpolation fraction for binary species parameter.
"""
@inline function compute_interp_frac_η(
    lkp::AbstractLookUp{I, FT},
    vmr,
    tropo,
    jtemp,
    ibnd,
    glay,
    gcol,
) where {FT <: AbstractFloat, I <: Int}
    (; n_η, key_species, vmr_ref) = lkp
    ig = view(key_species, :, tropo, ibnd)

    vmr1 = get_vmr(vmr, ig[1], glay, gcol)
    vmr2 = get_vmr(vmr, ig[2], glay, gcol)

    itemp = 1
    @inbounds η_half = vmr_ref[tropo, ig[1] + 1, jtemp + itemp - 1] / vmr_ref[tropo, ig[2] + 1, jtemp + itemp - 1]
    col_mix1 = vmr1 + η_half * vmr2
    η = col_mix1 ≥ eps(FT) * 2 ? vmr1 / col_mix1 : FT(0.5)
    loc_η = FT(η * (n_η - 1))
    jη1 = min(unsafe_trunc(I, loc_η) + 1, n_η - 1)
    #fη1 = loc_η % FT(1) # TODO: "%: operator seems unstable on GPU
    #fη1 = FT(loc_η % 1) # to be revisited
    fη1 = loc_η - fld(loc_η, FT(1))

    itemp = 2
    @inbounds η_half = vmr_ref[tropo, ig[1] + 1, jtemp + itemp - 1] / vmr_ref[tropo, ig[2] + 1, jtemp + itemp - 1]
    col_mix2 = vmr1 + η_half * vmr2
    η = col_mix2 ≥ eps(FT) * 2 ? vmr1 / col_mix2 : FT(0.5)
    loc_η = FT(η * (n_η - 1))
    jη2 = min(unsafe_trunc(I, loc_η) + 1, n_η - 1)
    #fη2 = loc_η % FT(1) # TODO: "%" operator seems unstable on GPU
    #fη2 = FT(loc_η % 1) # to be revisited
    fη2 = loc_η - fld(loc_η, FT(1))

    return ((jη1, jη2, fη1, fη2), (col_mix1, col_mix2))#nothing
end

"""
    compute_τ_ssa_lw_src!(
        lkp::AbstractLookUp{I,FT},
        vmr,
        col_dry,
        igpt,
        ibnd,
        p_lay::FT,
        t_lay,
        glaycol,
        src_args...,
    ) where {FT<:AbstractFloat,I<:Int}

Compute optical thickness, single scattering albedo, asymmetry parameter
and longwave sources whenever applicable.
"""
@inline function compute_τ_ssa_lw_src!(
    lkp::AbstractLookUp{I, FT},
    vmr,
    col_dry,
    igpt,
    ibnd,
    p_lay::FT,
    t_lay,
    glaycol,
    src_args...,
) where {FT <: AbstractFloat, I <: Int}
    # upper/lower troposphere
    tropo = p_lay > lkp.p_ref_tropo ? 1 : 2
    # volume mixing ration of h2o
    vmr_h2o = get_vmr(vmr, lkp.idx_h2o, glaycol...)

    jftemp, jfpress, jfη, col_mix = compute_interp_fractions(lkp, vmr, p_lay, t_lay, tropo, ibnd, glaycol)

    # computing τ_major
    τ_major = interp3d(jfη..., jftemp..., jfpress..., lkp.kmajor, igpt, col_mix...) * col_dry
    # computing τ_minor
    τ_minor =
        compute_τ_minor(lkp, tropo, vmr, vmr_h2o, col_dry, p_lay, t_lay, jftemp..., jfη..., igpt, ibnd, glaycol...)
    # compute τ_Rayleigh
    τ_ray = compute_τ_rayleigh(lkp, tropo, col_dry, vmr_h2o, jftemp..., jfη..., igpt)
    τ = τ_major + τ_minor + τ_ray
    ssa = FT(0)
    if τ > 2 * eps(FT) # single scattering albedo
        ssa = τ_ray / τ
    end
    # computing Planck sources for longwave problem
    compute_lw_planck_src!(lkp, jfη..., jfpress..., jftemp..., t_lay, igpt, ibnd, glaycol..., src_args...)
    return (τ, ssa)
end

"""
    compute_τ_minor(
        lkp::AbstractLookUp,
        tropo::I,
        vmr,
        vmr_h2o::FT,
        col_dry,
        p_lay::FT,
        t_lay::FT,
        jtemp::I,
        ftemp::FT,
        jη1::I,
        jη2::I,
        fη1::FT,
        fη2::FT,
        igpt,
        ibnd,
        glay,
        gcol,
    ) where {FT<:AbstractFloat,I<:Int}

Compute optical thickness contributions from minor gases.
"""
@inline function compute_τ_minor(
    lkp::AbstractLookUp,
    tropo::I,
    vmr,
    vmr_h2o::FT,
    col_dry,
    p_lay::FT,
    t_lay::FT,
    jtemp::I,
    ftemp::FT,
    jη1::I,
    jη2::I,
    fη1::FT,
    fη2::FT,
    igpt,
    ibnd,
    glay,
    gcol,
) where {FT <: AbstractFloat, I <: Int}

    if tropo == 1 # in lower atmosphere
        minor_bnd_st = lkp.minor_lower_bnd_st
        idx_gases_minor = lkp.idx_gases_minor_lower
        minor_scales_with_density = lkp.minor_lower_scales_with_density
        idx_scaling_gas = lkp.idx_scaling_gas_lower
        scale_by_complement = lkp.lower_scale_by_complement
        minor_gpt_sh = lkp.minor_lower_gpt_sh
        kminor = lkp.kminor_lower
    else # in upper atmosphere
        minor_bnd_st = lkp.minor_upper_bnd_st
        idx_gases_minor = lkp.idx_gases_minor_upper
        minor_scales_with_density = lkp.minor_upper_scales_with_density
        idx_scaling_gas = lkp.idx_scaling_gas_upper
        scale_by_complement = lkp.upper_scale_by_complement
        minor_gpt_sh = lkp.minor_upper_gpt_sh
        kminor = lkp.kminor_upper
    end

    τ_minor = FT(0)
    pa2hpa = FT(0.01) # pascals to hectopascals
    dry_fact = FT(1) / (FT(1) + vmr_h2o)

    @inbounds loc_in_bnd = igpt - (lkp.bnd_lims_gpt[1, ibnd] - 1)

    @inbounds for i in minor_bnd_st[ibnd]:(minor_bnd_st[ibnd + 1] - 1)
        vmr_imnr = get_vmr(vmr, idx_gases_minor[i], glay, gcol)
        if vmr_imnr > eps(FT) * 2
            scaling = vmr_imnr * col_dry

            if minor_scales_with_density[i] == 1
                scaling *= (pa2hpa * p_lay / t_lay)
                sgas = idx_scaling_gas[i]
                if sgas > 0
                    if scale_by_complement[i] == 1
                        scaling *= (FT(1) - get_vmr(vmr, sgas, glay, gcol) * dry_fact)
                    else
                        scaling *= get_vmr(vmr, sgas, glay, gcol) * dry_fact
                    end
                end
            end
            k_loc = minor_gpt_sh[i] + loc_in_bnd
            τ_minor += interp2d(fη1, fη2, ftemp, kminor, k_loc, jη1, jη2, jtemp) * scaling
        end
    end

    return τ_minor
end

"""
    compute_τ_rayleigh(
        lkp::LookUpSW,
        tropo::I,
        col_dry::FT,
        vmr_h2o::FT,
        jtemp::I,
        ftemp::FT,
        jη1::I,
        jη2::I,
        fη1::FT,
        fη2::FT,
        igpt::I,
    ) where {FT<:AbstractFloat,I<:Int}

Compute Rayleigh scattering optical depths for shortwave problem
"""
@inline function compute_τ_rayleigh(
    lkp::LookUpSW,
    tropo::I,
    col_dry::FT,
    vmr_h2o::FT,
    jtemp::I,
    ftemp::FT,
    jη1::I,
    jη2::I,
    fη1::FT,
    fη2::FT,
    igpt::I,
) where {FT <: AbstractFloat, I <: Int}
    if tropo == 1
        τ_ray = interp2d(fη1, fη2, ftemp, lkp.rayl_lower, igpt, jη1, jη2, jtemp) * (vmr_h2o + FT(1)) * col_dry
    else
        τ_ray = interp2d(fη1, fη2, ftemp, lkp.rayl_upper, igpt, jη1, jη2, jtemp) * (vmr_h2o + FT(1)) * col_dry
    end

    return τ_ray
end

@inline function compute_τ_rayleigh(lkp::LookUpLW{I, FT}, args...) where {FT <: AbstractFloat, I <: Int}
    return FT(0)
end

"""
    compute_lw_planck_src!(
        lkp::LookUpLW,
        jη1,
        jη2,
        fη1,
        fη2,
        jpresst,
        fpress,
        jtemp,
        ftemp,
        t_lay,
        igpt,
        ibnd,
        glay,
        gcol,
        sf,
        t_lev,
        t_sfc,
    )

Computes Planck sources for the longwave problem.

"""
@inline function compute_lw_planck_src!(
    lkp::LookUpLW,
    jη1,
    jη2,
    fη1,
    fη2,
    jpresst,
    fpress,
    jtemp,
    ftemp,
    t_lay,
    igpt,
    ibnd,
    glay,
    gcol,
    sf,
    t_lev,
    t_sfc,
)
    (; planck_fraction, t_planck, n_t_plnk, totplnk) = lkp
    (; lay_source, lev_source_inc, lev_source_dec, sfc_source) = sf
    # compute Planck fraction
    p_frac = interp3d(jη1, jη2, fη1, fη2, jtemp, ftemp, jpresst, fpress, planck_fraction, igpt)

    planck_args = (t_planck, totplnk, ibnd)
    # computing lay_source
    @inbounds lay_source[glay, gcol] = interp1d(t_lay, planck_args...) * p_frac
    # computing lev_source_inc
    @inbounds lev_source_inc[glay, gcol] = interp1d(t_lev[glay + 1, gcol], planck_args...) * p_frac
    # computing lev_source_dec
    @inbounds lev_source_dec[glay, gcol] = interp1d(t_lev[glay, gcol], planck_args...) * p_frac
    if glay == 1 # computing sfc_source
        @inbounds sfc_source[gcol] = interp1d(t_sfc, planck_args...) * p_frac
    end
    return nothing
end

@inline function compute_lw_planck_src!(lkp::LookUpSW, args...)
    return nothing
end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
""" Provide functions to communicate with the module `GeneticAlgorithms` """

# Public functions :
# - GAINTERFACE_exportProb : Convert `BiOKP` to `_MOMKP`
# - GAINTERFACE_exportLB : Convert `LowerBound` into a `front` for `GeneticAlgorithms`
# - GAINTERFACE_importLB : Convert a `front` for `GeneticAlgorithms` to a `LowerBound`

# Export functions convert types of the code to the types of the module GeneticAlgorithms
# Import functions convert types of the module GeneticAlgorithms to the types of the code

#-------------------------------------------------------------------#
#------------------------ PUBLIC FUNCTIONS -------------------------#
#-------------------------------------------------------------------#

"""
Convert `BiOKP` to `_MOMKP`
"""
function GAINTERFACE_exportProb(prob::BiOKP)
    @assert prob.isInteger "A _MOMKP object can only be integer"

    function weights_int(weights::Vector{Float64})
        W = Matrix{Int}(undef, 1, length(weights))
        for i in 1:length(weights)
            W[1,i] = Int(weights[i])
        end
        return W
    end

    return GeneticAlgorithms._MOMKP("myMOMKP",
                prob.nbVar,
                1,
                round.(Int, prob.profits),
                weights_int(prob.weights),
                [Int(prob.maxWeight)]
            )
end

"""
Convert `LowerBound` into a `front` for `GeneticAlgorithms`
"""
function GAINTERFACE_exportLB(prob::BiOKP, LB::LowerBound)
    n = length(LB.sols)
    LB_vec = Vector{GeneticAlgorithms.Sol}(undef, n)

    LB_aux = LB.sols # (LinkedList{Sol})
    iter = 1
    while LB_aux.tail != nil(Sol)
        # fill the LB (array)
        LB_vec[iter] = GAINTERFACE_exportSol(prob, LB_aux.head)

        LB_aux = LB_aux.tail # go next
        iter += 1
    end
    LB_vec[iter] = GAINTERFACE_exportSol(prob, LB_aux.head)

    return LB_vec
end

"""
Convert a `front` for `GeneticAlgorithms` to a `LowerBound`
"""
function GAINTERFACE_importLB(prob::BiOKP, front::Vector{GeneticAlgorithms.Sol})
    sort!(front, by = sol -> sol.z[1], rev = true) # the construction of the LB will reverse it again

    LB = LowerBound()

    LB.sols = cons(GAINTERFACE_importSol(front[1]), LB.sols)

    for iter in 2:length(front)
        LB.sols = cons(GAINTERFACE_importSol(front[iter]), LB.sols)
    end

    CONFIG.debug && DEBUG_LB(prob, LB)

    return LB
end

#-------------------------------------------------------------------#
#----------------------- PRIVATE FUNCTIONS -------------------------#
#-------------------------------------------------------------------#

"""
Convert type `Sol` to type `GeneticAlgorithms.Sol`
"""
function GAINTERFACE_exportSol(prob::BiOKP, sol::Sol)
    @assert sol.isBinary "GeneticAlgorithms module doesn't support non binary solutions"
    return GeneticAlgorithms.Sol(GAINTERFACE_exportProb(prob), round.(Int, sol.x), round.(Int, sol.y), [round(Int, sol.w)])
end

"""
Convert type `GeneticAlgorithms.Sol` to type `Sol`
"""
function GAINTERFACE_importSol(sol::GeneticAlgorithms.Sol)
    return Sol(sol.x, sol.z, sol.w[1], true)
end


text:
using NonconvexCore, LinearAlgebra, Test, Zygote

@testset "OrderedDict of reals" begin
    f(x) = sqrt(x[:b])
    g(x, a, b) = (a*x[:a] + b)^3 - x[:b]

    m = DictModel(f)
    addvar!(m, :a, 0.0, 10.0)
    addvar!(m, :b, 0.0, 10.0)
    add_ineq_constraint!(m, x -> g(x, 2, 0))
    add_ineq_constraint!(m, x -> g(x, -1, 1))

    @test getmin(m) == OrderedDict(:a => 0.0, :b => 0.0)
    @test getmax(m) == OrderedDict(:a => 10.0, :b => 10.0)
    @test getmin(m, :a) == 0.0
    @test getmax(m, :a) == 10.0
    @test getmin(m, :b) == 0.0
    @test getmax(m, :b) == 10.0

    p = OrderedDict(:a => 1.234, :b => 2.345)

    vec_model, _, un = NonconvexCore.tovecmodel(m)
    pvec = NonconvexCore.flatten(p)[1]
    val0, grad0 = NonconvexCore.value_gradient(NonconvexCore.getobjective(vec_model), pvec)
    @test val0 == f(un(pvec))
    @test grad0 == Zygote.gradient(x -> f(un(x)), pvec)[1]

    val1, grad1 = NonconvexCore.value_gradient(NonconvexCore.getineqconstraint(vec_model, 1), pvec)
    @test val1 == g(un(pvec), 2, 0)
    @test grad1 == Zygote.gradient(pvec -> g(un(pvec), 2, 0), pvec)[1]

    val2, grad2 = NonconvexCore.value_gradient(NonconvexCore.getineqconstraint(vec_model, 2), pvec)
    @test val2 == g(un(pvec), -1, 1)
    @test grad2 == Zygote.gradient(pvec -> g(un(pvec), -1, 1), pvec)[1]

    vals, jac = NonconvexCore.value_jacobian(NonconvexCore.getineqconstraints(vec_model), pvec)
    @test [val1, val2] == vals
    @test [grad1 grad2]' == jac

    vals, jac = NonconvexCore.value_jacobian(NonconvexCore.getobjectiveconstraints(vec_model), pvec)
    @test [val0, val1, val2] == vals
    @test [grad0 grad1 grad2]' == jac
end

struct S
    a
    b
end
@testset "Vector of structs" begin
    f(x) = sqrt(x[1].b)
    g(x, a, b) = (a*x[1].a + b)^3 - x[1].b

    m = Model(f)
    addvar!(m, S(0.0, 0.0), S(10.0, 10.0))
    add_ineq_constraint!(m, x -> g(x, 2, 0))
    add_ineq_constraint!(m, x -> g(x, -1, 1))

    p = [S(1.234, 2.345)]

    vec_model, _, un = NonconvexCore.tovecmodel(m)
    pvec, un = NonconvexCore.flatten(p)
    val0, grad0 = NonconvexCore.value_gradient(NonconvexCore.getobjective(vec_model), pvec)
    @test val0 == f(un(pvec))
    @test grad0 == Zygote.gradient(x -> f(un(x)), pvec)[1]

    val1, grad1 = NonconvexCore.value_gradient(NonconvexCore.getineqconstraint(vec_model, 1), pvec)
    @test val1 == g(un(pvec), 2, 0)
    @test grad1 == Zygote.gradient(pvec -> g(un(pvec), 2, 0), pvec)[1]

    val2, grad2 = NonconvexCore.value_gradient(NonconvexCore.getineqconstraint(vec_model, 2), pvec)
    @test val2 == g(un(pvec), -1, 1)
    @test grad2 == Zygote.gradient(pvec -> g(un(pvec), -1, 1), pvec)[1]

    vals, jac = NonconvexCore.value_jacobian(NonconvexCore.getineqconstraints(vec_model), pvec)
    @test [val1, val2] == vals
    @test [grad1 grad2]' == jac

    vals, jac = NonconvexCore.value_jacobian(NonconvexCore.getobjectiveconstraints(vec_model), pvec)
    @test [val0, val1, val2] == vals
    @test [grad0 grad1 grad2]' == jac
end


text:
# This file is for types in chp.c

type CHPState{N}<:AbstractQuState{N}
    X::AbstractMatrix # (2n+1)*n matrix for stabilizer/destabilizer x bits
    Z::AbstractMatrix # (2n+1)*n matrix for z bits
    R::AbstractVector # phase bits: 0 for +1, 1 for i, 2 for -1, 3 for -i. Normally either 0 or 2
end

CHPState(X::BitMatrix,Z::BitMatrix,R::BitVector) = CHPState{length(R)/2}(X,Z,R)

function CHPState(n::Int)
    X = spzeros(Bool,2n,n)
    Z = spzeros(Bool,2n,n)
    R = spzeros(Bool,2n)

    X[1:n,:] = speye(n)
    Z[n+1:2n,:] = speye(n)

    return CHPState{n}(X,Z,R)
end

function copy!(A::CHPState,B::CHPState)
    copy!(A.X,B.X)
    copy!(A.Z,B.Z)
    copy!(A.R,B.R)

    return A
end

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
#=
Solves the price function for the Lucas tree in a continuous state
setting, using piecewise linear approximation for the sequence of
candidate price functions.  The consumption endownment follows the
log linear AR(1) process

    log y' = alpha log y + sigma epsilon

where y' is a next period y and epsilon is an iid standard normal
shock. Hence

    y' = y^alpha * xi   where xi = e^(sigma * epsilon)

The distribution phi of xi is

    phi = LN(0, sigma^2) where LN means lognormal

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date : 2014-07-05

References
----------

http://quant-econ.net/jl/markov_asset.html

TODO: refactor. Python is much cleaner.
=#
using Distributions
using Grid: CoordInterpGrid, BCnan, BCnearest, InterpLinear

"""
The Lucas asset pricing model

##### Fields

-  `gam::Real` : coefficient of risk aversion in the CRRA utility function
-  `bet::Real` : Discount factor in (0, 1)
-  `alpha::Real` : Correlation coefficient in the shock process
-  `sigma::Real` : Volatility of shock process
-  `phi::Distribution` : Distribution for shock process
-  `grid::AbstractVector` : Grid of points on which to evaluate the prices. Each
point should be non-negative
-  `grid_min::Real` : Lower bound on grid
-  `grid_max::Real` : Upper bound on grid
-  `grid_size::Int` : Number of points in the grid
- `quad_nodes::Vector` : Quadrature nodes for integrating over the shock
- `quad_weights::Vector` : Quadrature weights for integrating over the shock
-  `h::Vector` : Storage array for the `h` vector in the lucas operator
"""
type LucasTree
    gam::Real
    bet::Real
    alpha::Real
    sigma::Real
    phi::Distribution
    grid::AbstractVector
    grid_min::Real
    grid_max::Real
    grid_size::Int
    quad_nodes::Vector
    quad_weights::Vector
    h::Vector

    # this needs to be an internal constructor because we need to incompletely
    # initialize the object before we can compute h.
    """
    Constructor for LucasTree

    ##### Arguments

    -  `gam::Real` : coefficient of risk aversion in the CRRA utility function
    -  `bet::Real` : Discount factor in (0, 1)
    -  `alpha::Real` : Correlation coefficient in the shock process
    -  `sigma::Real` : Volatility of shock process

    ##### Notes

    All other fields of the type are instantiated within the constructor
    """
    function LucasTree(gam::Real, bet::Real, alpha::Real, sigma::Real)
        phi = LogNormal(0.0, sigma)
        grid = make_grid(alpha, sigma)
        grid_min, grid_max, grid_size = minimum(grid), maximum(grid), length(grid)
        _int_min, _int_max = exp(-4 * sigma), exp(4 * sigma)
        n, w = qnwlege(21, _int_min, _int_max)

        # create lt object without h
        lt = new(gam, bet, alpha, sigma, phi, grid, grid_min, grid_max, grid_size,
                 n, w)

        # initialize h
        h = Array(Float64, grid_size)

        for (i, y) in enumerate(grid)
            integrand(z) = (y^alpha * z).^(1 - gam)
            h[i] = bet .* integrate(lt, integrand)
        end

        # now add h to it
        lt.h = h
        lt
    end
end

function make_grid(alpha, sigma)
    grid_size = 100
    if abs(alpha) >= 1
        grid_min, grid_max = 0.0, 10.
    else
        # Set the grid interval to contain most of the mass of the
        # stationary distribution of the consumption endowment
        ssd = sigma / sqrt(1 - alpha^2)
        grid_min, grid_max = exp(-4 * ssd), exp(4 * ssd)
    end
    return linspace_range(grid_min, grid_max, grid_size)
end

function integrate(lt::LucasTree, g::Function, qn::Array=lt.quad_nodes,
                   qw::Array=lt.quad_weights)
    int_func(x) = g(x) .* pdf(lt.phi, x)
    return do_quad(int_func, qn, qw)
end

"""
The approximate Lucas operator, which computes and returns the updated function
Tf on the grid points.

##### Arguments

- `lt::LucasTree` : An instance of the `LucasTree` type
- `f::Vector{Float64}` : A candidate function on R_+ represented as points on a
grid. It should be the same size as `lt.grid`

##### Returns

- `Tf::Vector{Float64}` : The updated function Tf

"""
function lucas_operator(lt::LucasTree, f::AbstractVector)
    grid, h, alpha, bet = lt.grid, lt.h, lt.alpha, lt.bet

    Tf = similar(f)
    Af = CoordInterpGrid(grid, f, BCnearest, InterpLinear)

    for (i, y) in enumerate(grid)
        to_integrate(z) = Af[y^alpha * z]
        Tf[i] = h[i] + bet * integrate(lt, to_integrate)
    end
    return Tf
end

"""
Compute the equilibrium price function associated with Lucas tree `lt`

##### Arguments

- `lt::LucasTree` : An instance of the `LucasTree` type
- `;kwargs...` : other arguments to be passed to `compute_fixed_point`

##### Returns

- `price::Vector{Float64}` : The price at each point in `lt.grid`

"""
function compute_lt_price(lt::LucasTree; kwargs...)
    # Simplify names, set up distribution phi
    grid, grid_size, gam = lt.grid, lt.grid_size, lt.gam

    f_init = zeros(grid)  # Initial condition
    func(x) = lucas_operator(lt, x)
    f = compute_fixed_point(func, f_init; max_iter=500, verbose=false, kwargs...)

    # p(y) = f(y) / u'(y) = f(y) * y^gamma
    price = f .* grid.^(gam)

    return price
end


text:
function swapmatrix(N::Integer)
   Dim  =   N*N;
   W    =   zeros(Dim,Dim);
   for k = 1:Dim
       j = k % N;
       if j == 0
           j = N;
       end
       i = convert(Int64,(k - j) / N + 1);
       W[(j-1)*N+i,k] = 1;
    end
    return W;
end

max_stars_repo_name: askuyue/wjOperators

max_issues_repo_name: askuyue/wjOperators

max_forks_repo_name: askuyue/wjOperators


text:
export ScoreActor

using Rocket

mutable struct ScoreActorProps
    count       :: Int
    is_released :: Bool
end

struct ScoreActor{L} <: Rocket.Actor{L}
    score :: Vector{L}
    inds  :: Vector{Int}
    props :: ScoreActorProps
end

ScoreActor()                              = ScoreActor(Real)
ScoreActor(::Type{L}) where { L <: Real } = ScoreActor{L}(Vector{L}(), Vector{Int}(), ScoreActorProps(0, true))

Base.show(io::IO, ::ScoreActor) = print(io, "ScoreActor()")

function Rocket.on_next!(actor::ScoreActor{L}, data::L) where L
    actor.props.count += 1
    actor.props.is_released = false
    push!(actor.score, data)
end

function Rocket.on_error!(actor::ScoreActor, err)
    error(err)
end

function Rocket.on_complete!(actor::ScoreActor)
    Rocket.release!(actor)
    nothing
end

function Rocket.release!(actor::ScoreActor)
    if !actor.props.is_released
        push!(actor.inds, actor.props.count)
        actor.props.is_released = true
    end
end

function Rocket.getvalues(actor::ScoreActor)
    if length(actor.inds) <= 1
        return actor.score
    end

    n        = length(actor.inds)
    indices  = Iterators.flatten((0, actor.inds))

    # Allocation free pairwise
    pairwise = zip(Iterators.take(indices, n), Iterators.take(Iterators.drop(indices, 1), n))

    columns = map(pairwise) do (left, right)
        return view(actor.score, (left+1):right)
    end

    maxlength = maximum(length.(columns))

    result = zeros(eltype(actor.score), maxlength)

    foreach(columns) do column
        broadcast!(+, result, result, Iterators.flatten((column, Iterators.repeated(last(column), maxlength - length(column)))))
    end

    return result
end