text:
module Zeptobot

using HTTP, JSON, GitHub, Dates, Sockets, Logging

USERNAME  = Ref{String}()
TOKEN     = Ref{String}()
LOGDIR    = Ref{String}()
REGISTRATION_OK_AGE = Day(3)

function __init__()
    if haskey(ENV, "GITHUB_USERNAME")
        USERNAME[] = ENV["GITHUB_USERNAME"]
    else
        print("Please enter a GitHub username: ")
        USERNAME[] = readline(stdin)
    end

    if haskey(ENV, "GITHUB_TOKEN")
        TOKEN[] = ENV["GITHUB_TOKEN"]
    else
        print("Please enter a GitHub token: ")
        TOKEN[] = readline(stdin)
    end

    if haskey(ENV, "ZEPTO_LOGDIR")
        LOGDIR[] = ENV["ZEPTO_LOGDIR"]
    else
        print("Please enter a path to the log directory: ")
        LOGDIR[] = readline(stdin)
    end
end

function parsepageinfo(s::AbstractString)
    ss = split(s, ',')
    rels_n_pages = map(t -> (strip(match(r"(?<=rel=).*", t).match, '\"'),
                             strip(match(r"(?<=<)(.*)(?=>)", t).match, '\"')),
                       ss)
    return Dict(rels_n_pages)
end

function _requestURL(s::AbstractString, type = "GET"; kw...)
    d = Dict(kw)
    r = HTTP.request(type,
                     s,
                     Dict("User-Agent" => USERNAME[],
                          "Authorization" => "token $(TOKEN[])",
                          "Accept" => "application/vnd.github.v3+json"),
                     # FixMe! We should probably strip all parameters from
                     # s as well.
                     query = isempty(d) ? nothing : d,
                     retry = true,
                     retries = 10,
                     status_exception = false)

    if r.status >= 300
        return r, false
    end

    return r, true
end

function getURLallpages(s::AbstractString; kw...)
    if haskey(kw, :page)
        # @error "Don't pass page argument since this function tranverses all pages"
        error("Don't pass page argument since this function tranverses all pages")
    end

    entries = []

    r, noerr = _requestURL(s, "GET", kw...)

    # Traverse the pages
    while true
        # For now just error if request didn't succeed
        if noerr
            # If success then parse the body
            append!(entries, JSON.parse(String(r.body)))

            # Parse the page info in the header
            d = Dict(r.headers)
            if haskey(d, "Link")
                pd = parsepageinfo(d["Link"])
            else
                @debug "Single page document"
                break
            end

            # If the header doesn't contain info about the next page then we are done
            if !haskey(pd, "next")
                @debug "No more pages to read"
                break
            else
                # else read the next page
                @debug "Reading page $(match(r"(?<=page=)[0-9]", pd["next"]).match)"
                r, noerr = _requestURL(pd["next"], "GET")
            end
        else
            @error "Request returned an error"
            break
        end
    end

    return entries
end

getPRs(;repo = "JuliaLang/METADATA.jl") = pull_requests(repo, params = Dict("state" => "open"))

getStatuses(pr::PullRequest) = getURLallpages(pr._links["statuses"]["href"])

getComments(pr::PullRequest) = getURLallpages(pr._links["comments"]["href"])

function checkComments(pr::PullRequest)
    comments = getComments(pr)
    # It's okay to merge if there are no comments or if the latest comment
    # is made by Attobot since it indicates that the PR has been updated
    return isempty(comments) || last(comments)["user"]["login"] == "attobot"
end

getLabels(pr::PullRequest) = getURLallpages(pr._links["issue"]["href"]*"/labels")

function checkLabels(pr)
    labels = getLabels(pr)
    # For now we abort on any label
    return isempty(labels)
end

function mergeable(pr::PullRequest)
    # Only do something with PRs not made by Attobot
    if pr.user.login != "attobot"
        @info "Skipping: PR not by Attobot"
        return false
    end

    # Check for comments
    @info "Checking for comments"
    if !checkComments(pr)
        @info "Skipping: PR contains comments"
        return false
    end

    # Check test status
    @info "Checking test status"
    statuses = map(t -> t["context"], filter(t -> t["state"] == "success", getStatuses(pr)))
    if !(("JuliaCIBot" ∈ statuses) &&
         ("continuous-integration/travis-ci/pr" ∈ statuses))
        @info "Skipping: Tests failed or still in progess"
        return false
    end

    # Check labels
    @info "Checking issue labels"
    if !checkLabels(pr)
        @info "Skipping: PR has labels attached"
        return false
    end

    # Handle releases (i.e. not new registrations)
    if occursin(r"Tag .*\.jl v[\d.\d.\d]", pr.title)

        @info "PR tries tag new release"

        return true

    elseif occursin(r"Register new package .*\.jl v[\d.\d.\d]", pr.title)

        @info "PR tries to register new package"

        if (now(UTC) - pr.created_at) >= REGISTRATION_OK_AGE
            return true
        else
            @info "Skipping: PR is younger than $REGISTRATION_OK_AGE"
            return false
        end
    else
        throw(ErrorException("this should never happen!"))
    end
end

function merge(pr::PullRequest;
    title = "$(pr.title) [$(match(r"(?<=\()(.*)(?=\))",pr.body).match)] (#$(pr.number))",
    comment = "Merged automatically by Zeptobot")

    r = HTTP.request(
        "PUT",
        "http://api.github.com/repos/$(pr.base.repo.full_name)/pulls/$(pr.number)/merge",
        Dict("User-Agent" => USERNAME[],
             "Authorization" => "token $(TOKEN[])",
             "Accept" => "application/vnd.github.v3+json"),
        JSON.json(
            Dict(
                "commit_title" => title,
                "commit_message" => comment,
                "merge_method" => "squash")
        ),
        retry = true,
        retries = 10,
        status_exception = false
    )

    return r
end

# FixMe! Not implemented yet
closeable(pr) = false

function process(;repo = "JuliaLang/METADATA.jl", prs::Vector{PullRequest} = pull_requests(repo, params = Dict("state" => "open"))[1], dryrun = false)

    if dryrun
        @warn "Running in dry-run mode. No actions taken."
    end

    merge_count         = 0
    merge_count_success = 0
    close_count         = 0
    close_count_success = 0

    @info "Traverse the PRs to determine actions"
    for pr in prs

        title = pr.title

        @info "Processing \"$title\""

        if mergeable(pr)

            @info "PR can be merged"
            merge_count += 1
            if !dryrun
                r = merge(pr)
                if r.status >= 300
                    @error "Merging failed with: $(String(r.body))"
                else
                    @info "Merge successful!"
                    merge_count_success += 1
                end
            end
            @info ""
            sleep(10)

        elseif closeable(pr)

            # FixME. Closing not implemented yet!
            @info "PR can be closed"
            close_count += 1
            if !dryrun
                r = close(pr)
                if r.status >= 300
                    @error "Closing PR failed with: $(String(r.body))"
                else
                    @info "PR closed successfully!"
                    close_count_success += 1
                end
            end

        else

            # do nothing
            @info ""

        end
    end

    @info "PR processing complete"
    @info "$merge_count_success out of $merge_count PRs merged successfully"
    @info "$close_count_success out of $close_count PRs closed successfully"

    return nothing
end

function restartPending(;repo = "JuliaLang/METADATA.jl", prs::Vector{PullRequest} = pull_requests(repo, params = Dict("state" => "open"))[1])

    for pr in prs
        # Check test status
        @info "Checking test status"
        statuses = map(t -> t["context"], filter(t -> t["state"] == "success", getStatuses(pr)))
        if !(("JuliaCIBot" ∈ statuses) &&
             ("continuous-integration/travis-ci/pr" ∈ statuses))
            @info "Skipping: Tests failed or still in progess"
            return false
        end
    end
end

# A small server to listen for METADATA events and run process if necessary
function serveMETADATA(;repo = "JuliaLang/METADATA.jl", port = 6384, auth = GitHub.authenticate(Zeptobot.TOKEN[]), secret = "")

    # Initialize "global" counter for the logs
    logcounter = 0

    # Setup EventListener
    e = GitHub.EventListener(auth = auth, secret = secret) do event

        # Create daily logging directory if needed
        todaydir = joinpath(LOGDIR[], replace(string(today()), "-" => ""))
        if !isdir(todaydir)
            mkdir(todaydir)
            logcounter = 1
        elseif isempty(todaydir)
            logcounter = 1
        else
            # Find the latest job id
            logcounter = parse(Int, match(r"[0-9]{4}", last(sort(readdir(todaydir)))).match) + 1
        end

        # Setup logfile environment
        prefix = joinpath(todaydir, "job$(lpad(logcounter, 4, '0'))")
        open(prefix*"stderr.log", "w") do f
            with_logger(SimpleLogger(f)) do

            @info "Received an event from GitHub. Processing!"
            if event.kind != "status"
                @error "event is not a status update"
                return HTTP.Response(500)
            else
                name = event.payload["name"]
                if name != repo
                    @error "Received an event from the wrong repo"
                    return HTTP.Response(500)
                end

                if event.payload["state"] == "success"
                    @info "Event status: success. Trying to merge PRs..."
                    process(repo = repo)
                    return HTTP.Response(200)
                else
                    @info "Status: pending"
                    return HTTP.Response(200)
                end
            end
        end
        end
    end

    if !isdir(LOGDIR[])
        @info "Creating log directory"
        mkdir(LOGDIR[])
    end

    GitHub.run(e, getipaddr(), port)
    return nothing
end

end


text:
abstract type AbstractSimulation{T <: SSDFloat} end

"""
    mutable struct Simulation{T <: SSDFloat} <: AbstractSimulation{T}

Collection of all parts of a Simulation of a Solid State Detector.
"""
mutable struct Simulation{T <: SSDFloat} <: AbstractSimulation{T}
    detector::Union{SolidStateDetector{T}, Missing}
    ρ::Union{ChargeDensity{T}, Missing}
    ρ_fix::Union{ChargeDensity{T}, Missing}
    ϵ::Union{DielectricDistribution{T}, Missing}
    point_types::Union{PointTypes{T}, Missing}
    electric_potential::Union{ElectricPotential{T}, Missing}
    weighting_potentials::Vector{Any}
    electric_field::Union{ElectricField{T}, Missing}

    charge_drift_model::Union{<:AbstractChargeDriftModel{T}, Missing}

    electron_drift_field::Union{ElectricField{T}, Missing}
    hole_drift_field::Union{ElectricField{T}, Missing}
end

function Simulation{T}() where {T <: SSDFloat}
    Simulation{T}(
        SolidStateDetector{T}(),
        missing,
        missing,
        missing,
        missing,
        missing,
        [missing],
        missing,
        ElectricFieldChargeDriftModel{T}(),
        missing,
        missing
    )
end

function get_precision_type(sim::Simulation{T}) where {T}
    typeof(sim).parameters[1]
end

function NamedTuple(sim::Simulation{T}) where {T <: SSDFloat}
    wps_strings = AbstractString[]
    for contact in sim.detector.contacts
        if !ismissing(sim.weighting_potentials[contact.id])
            push!(wps_strings, "WeightingPotential_$(contact.id)")
        end
    end
    wps_syms = Symbol.(wps_strings)
    nt = (
        detector_json_string = NamedTuple(sim.detector.config_dict),
        electric_potential = NamedTuple(sim.electric_potential),
        ρ = NamedTuple(sim.ρ),
        ρ_fix = NamedTuple(sim.ρ_fix),
        ϵ = NamedTuple(sim.ϵ),
        point_types = NamedTuple(sim.point_types),
        electric_field = NamedTuple(sim.electric_field),
        electron_drift_field = NamedTuple(sim.electron_drift_field),
        hole_drift_field = NamedTuple(sim.hole_drift_field)
    )
    if length(wps_strings) > 0
        nt = merge(nt, (weighting_potentials = NamedTuple{Tuple(wps_syms)}(NamedTuple.( skipmissing(sim.weighting_potentials))),))
    end
    return nt
end
Base.convert(T::Type{NamedTuple}, x::Simulation) = T(x)

function Simulation(nt::NamedTuple)
    ep = ElectricPotential(nt.electric_potential)
    T = eltype(ep.data)
    det = SolidStateDetector{T}( Dict(nt.detector_json_string) )
    sim = Simulation( det )
    if !ismissing(nt.electric_potential) sim.electric_potential = ep end
    if !ismissing(nt.ρ) sim.ρ = ChargeDensity(nt.ρ) end
    if !ismissing(nt.ρ_fix) sim.ρ_fix = ChargeDensity(nt.ρ_fix) end
    if !ismissing(nt.ϵ) sim.ϵ = DielectricDistribution(nt.ϵ) end
    if !ismissing(nt.point_types) sim.point_types = PointTypes(nt.point_types) end
    if !ismissing(nt.electric_field) sim.electric_field = ElectricField(nt.electric_field) end
    sim.weighting_potentials = [missing for contact in sim.detector.contacts]
    if !ismissing(nt.weighting_potentials)
        for contact in sim.detector.contacts
            if !ismissing(values(nt.weighting_potentials[contact.id])[1])
                sim.weighting_potentials[contact.id] = WeightingPotential(nt.weighting_potentials[contact.id])
            end
        end
    end
    if !ismissing(nt.electron_drift_field) sim.electron_drift_field = ElectricField(nt.electron_drift_field) end
    if !ismissing(nt.hole_drift_field) sim.hole_drift_field = ElectricField(nt.hole_drift_field) end
    sim.charge_drift_model = ADLChargeDriftModel(T=T)
    @info "I/O of charge drift model not yet supported. Loading default: ADLChargeDriftModel"
    return sim
end
Base.convert(T::Type{Simulation}, x::NamedTuple) = T(x)

function println(io::IO, sim::Simulation{T}) where {T <: SSDFloat}
    println(typeof(sim), " - Coordinate system: ", get_coordinate_system(sim.detector))
    println("  Detector: $(sim.detector.name)")
    println("  Electric potential: ", !ismissing(sim.electric_potential) ? size(sim.electric_potential) : missing)
    println("  Charge density: ", !ismissing(sim.ρ) ? size(sim.ρ) : missing)
    println("  Fix Charge density: ", !ismissing(sim.ρ_fix) ? size(sim.ρ_fix) : missing)
    println("  Dielectric distribution: ", !ismissing(sim.ϵ) ? size(sim.ϵ) : missing)
    println("  Point types: ", !ismissing(sim.point_types) ? size(sim.point_types) : missing)
    println("  Electric field: ", !ismissing(sim.electric_field) ? size(sim.electric_field) : missing)
    println("  Weighting potentials: ")
    for contact in sim.detector.contacts
        print("    Contact $(contact.id): ")
        println(!ismissing(sim.weighting_potentials[contact.id]) ? size(sim.weighting_potentials[contact.id]) : missing)
    end
    println("  Charge drift model: ", !ismissing(sim.electric_field) ? typeof(sim.charge_drift_model) : missing)
    println("  Electron drift field: ", !ismissing(sim.electron_drift_field) ? size(sim.electron_drift_field) : missing)
    println("  Hole drift field: ", !ismissing(sim.hole_drift_field) ? size(sim.hole_drift_field) : missing)
end

function print(io::IO, sim::Simulation{T}) where {T <: SSDFloat}
    print(io, "Simulation{$T} - ", "$(sim.detector.name)")
end

function show(io::IO, sim::Simulation{T}) where {T <: SSDFloat} println(io, sim) end

function show(io::IO, ::MIME"text/plain", sim::Simulation{T}) where {T <: SSDFloat}
    show(io, sim)
end

function Simulation(detector::SolidStateDetector{T})::Simulation{T} where {T <: SSDFloat}
    sim::Simulation{T} = Simulation{T}()
    sim.detector = detector
    sim.weighting_potentials = Missing[ missing for i in 1:length(sim.detector.contacts)]
    return sim
end

function Simulation{T}(config_file::AbstractString)::Simulation{T} where{T <: SSDFloat}
    return Simulation( SolidStateDetector{T}(config_file) )
end
function Simulation(config_file::AbstractString)::Simulation{Float32}
    return Simulation{Float32}( config_file )
end

# Functions
"""
    function apply_initial_state!(sim::Simulation{T}, ::Type{ElectricPotential}, grid::Grid{T} = Grid(sim.detector))::Nothing

Applies the initial state of the electric potential calculation.
It overwrites `sim.electric_potential`, `sim.ρ`, `sim.ρ_fix`, `sim.ϵ` and `sim.point_types`.
"""
function apply_initial_state!(sim::Simulation{T}, ::Type{ElectricPotential}, grid::Grid{T} = Grid(sim.detector))::Nothing where {T <: SSDFloat}
    fssrb::PotentialSimulationSetupRB{T, 3, 4, get_coordinate_system(sim.detector)} =
        PotentialSimulationSetupRB(sim.detector, grid);

    sim.ρ = ChargeDensity(ChargeDensityArray(fssrb), grid)
    sim.ρ_fix = ChargeDensity(FixedChargeDensityArray(fssrb), grid)
    sim.ϵ = DielectricDistribution(DielektrikumDistributionArray(fssrb), grid)
    sim.point_types = PointTypes(PointTypeArray(fssrb), grid)
    sim.electric_potential = ElectricPotential(ElectricPotentialArray(fssrb), grid)
    nothing
end

"""
    function apply_initial_state!(sim::Simulation{T}, ::Type{WeightingPotential}, contact_id::Int, grid::Grid{T} = Grid(sim.detector))::Nothing

Applies the initial state of the weighting potential calculation for the contact with the id `contact_id`.
It overwrites `sim.weighting_potentials[contact_id]`.
"""
function apply_initial_state!(sim::Simulation{T}, ::Type{WeightingPotential}, contact_id::Int, grid::Grid{T} = Grid(sim.detector))::Nothing where {T <: SSDFloat}
    fssrb::PotentialSimulationSetupRB{T, 3, 4, get_coordinate_system(sim.detector)} =
        PotentialSimulationSetupRB(sim.detector, grid, weighting_potential_contact_id = contact_id);

    sim.weighting_potentials[contact_id] = WeightingPotential(ElectricPotentialArray(fssrb), grid)
    nothing
end

"""
    function update_till_convergence!( sim::Simulation{T} ::Type{ElectricPotential}, convergence_limit::Real; kwargs...)::T

Takes the current state of `sim.electric_potential` and updates it until it has converged.
"""
function update_till_convergence!( sim::Simulation{T},
                                   ::Type{ElectricPotential},
                                   convergence_limit::Real = 1e-7;
                                   n_iterations_between_checks::Int = 500,
                                   max_n_iterations::Int = -1,
                                   depletion_handling::Bool = false,
                                   use_nthreads::Int = Base.Threads.nthreads(),
                                   sor_consts::Union{Missing, T, NTuple{2, T}} = missing
                                    )::T where {T <: SSDFloat}
    CS = get_coordinate_system(sim.detector)
    if ismissing(sor_consts)
        sor_consts = CS == :cylindrical ? (T(1.4), T(1.85)) : T(1.4)
    elseif length(sor_consts) == 1 && CS == :cylindrical
        sor_consts = (T(sor_consts), T(sor_consts))
    elseif length(sor_consts) > 1 && CS == :cartesian
        sor_consts = T(sor_consts[1])
    end
    only_2d::Bool = length(sim.electric_potential.grid.axes[2]) == 1

    fssrb::PotentialSimulationSetupRB{T, 3, 4, get_coordinate_system(sim.electric_potential.grid)} =
        PotentialSimulationSetupRB(sim.detector, sim.electric_potential.grid, sim.electric_potential.data, sor_consts = T.(sor_consts))

    cf::T = _update_till_convergence!( fssrb, T(convergence_limit);
                                       only2d = Val{only_2d}(),
                                       depletion_handling = Val{depletion_handling}(),
                                       is_weighting_potential = Val{false}(),
                                       use_nthreads = use_nthreads,
                                       n_iterations_between_checks = n_iterations_between_checks,
                                       max_n_iterations = max_n_iterations )

    grid::Grid = Grid(fssrb)
    sim.ρ = ChargeDensity(ChargeDensityArray(fssrb), grid)
    sim.ρ_fix = ChargeDensity(FixedChargeDensityArray(fssrb), grid)
    sim.ϵ = DielectricDistribution(DielektrikumDistributionArray(fssrb), grid)
    sim.electric_potential = ElectricPotential(ElectricPotentialArray(fssrb), grid)
    sim.point_types = PointTypes(PointTypeArray(fssrb), grid)

    if depletion_handling == false
        update_again::Bool = false # With SOR-Constant = 1
        if fssrb.bulk_is_ptype
            @inbounds for i in eachindex(sim.electric_potential.data)
                if sim.electric_potential.data[i] < fssrb.minimum_applied_potential
                    sim.electric_potential.data[i] = fssrb.minimum_applied_potential
                    update_again = true
                end
            end
        else # ntype
            @inbounds for i in eachindex(sim.electric_potential.data)
                if sim.electric_potential.data[i] > fssrb.maximum_applied_potential
                    sim.electric_potential.data[i] = fssrb.maximum_applied_potential
                    update_again = true
                end
            end
        end
        if update_again
            fssrb.sor_const[:] .= T(1)
            cf = _update_till_convergence!( fssrb, T(convergence_limit);
                                            only2d = Val{only_2d}(),
                                            depletion_handling = Val{depletion_handling}(),
                                            is_weighting_potential = Val{false}(),
                                            use_nthreads = use_nthreads,
                                            n_iterations_between_checks = n_iterations_between_checks,
                                            max_n_iterations = max_n_iterations )
            sim.electric_potential = ElectricPotential(ElectricPotentialArray(fssrb), grid)

            if fssrb.bulk_is_ptype
                @inbounds for i in eachindex(sim.electric_potential.data)
                    if sim.electric_potential.data[i] < fssrb.minimum_applied_potential
                        @warn """Detector seems not to be fully depleted at a bias voltage of $(fssrb.bias_voltage) V.
                            At least one grid point has a smaller potential value ($(sim.electric_potential.data[i]) V)
                            than the minimum applied potential ($(fssrb.minimum_applied_potential) V). This should not be.
                            However, small overshoots could be due to numerical precision."""
                        break
                    end
                end
            else # ntype
                @inbounds for i in eachindex(sim.electric_potential.data)
                    if sim.electric_potential.data[i] > fssrb.maximum_applied_potential
                        @warn """Detector seems not to be not fully depleted at a bias voltage of $(fssrb.bias_voltage) V.
                            At least one grid point has a higher potential value ($(sim.electric_potential.data[i]) V)
                            than the maximum applied potential ($(fssrb.maximum_applied_potential) V). This should not be.
                            However, small overshoots could be due to numerical precision."""
                        break
                    end
                end
            end
        end
    end

    cf
end

"""
    function update_till_convergence!( sim::Simulation{T} ::Type{WeightingPotential}, contact_id::Int, convergence_limit::Real; kwargs...)::T

Takes the current state of `sim.weighting_potentials[contact_id]` and updates it until it has converged.
"""
function update_till_convergence!( sim::Simulation{T},
                                   ::Type{WeightingPotential},
                                   contact_id::Int,
                                   convergence_limit::Real;
                                   n_iterations_between_checks::Int = 500,
                                   max_n_iterations::Int = -1,
                                   depletion_handling::Bool = false,
                                   use_nthreads::Int = Base.Threads.nthreads(),
                                   sor_consts::Union{Missing, T, NTuple{2, T}} = missing
                                    )::T where {T <: SSDFloat}
    CS = get_coordinate_system(sim.detector)
    if ismissing(sor_consts)
        sor_consts = CS == :cylindrical ? (T(1.4), T(1.85)) : T(1.4)
    elseif length(sor_consts) == 1 && CS == :cylindrical
        sor_consts = (T(sor_consts), T(sor_consts))
    elseif length(sor_consts) > 1 && CS == :cartesian
        sor_consts = T(sor_consts[1])
    end

    only_2d::Bool = length(sim.weighting_potentials[contact_id].grid.axes[2]) == 1
    fssrb::PotentialSimulationSetupRB{T, 3, 4, get_coordinate_system(sim.weighting_potentials[contact_id].grid)} =
        PotentialSimulationSetupRB(sim.detector, sim.weighting_potentials[contact_id].grid, sim.weighting_potentials[contact_id].data,
                sor_consts = T.(sor_consts), weighting_potential_contact_id = contact_id)

    cf::T = _update_till_convergence!( fssrb, T(convergence_limit);
                                       only2d = Val{only_2d}(),
                                       depletion_handling = Val{depletion_handling}(),
                                       is_weighting_potential = Val{true}(),
                                       use_nthreads = use_nthreads,
                                       n_iterations_between_checks = n_iterations_between_checks,
                                       max_n_iterations = max_n_iterations )

    sim.weighting_potentials[contact_id] = WeightingPotential(ElectricPotentialArray(fssrb), sim.weighting_potentials[contact_id].grid)

    cf
end

"""
    function refine!(sim::Simulation{T}, ::Type{ElectricPotential}, max_diffs::Tuple{<:Real,<:Real,<:Real}, minimum_distances::Tuple{<:Real,<:Real,<:Real})

Takes the current state of `sim.electric_potential` and refines it with respect to the input arguments
`max_diffs` and `minimum_distances`.
"""
function refine!(sim::Simulation{T}, ::Type{ElectricPotential},
                    max_diffs::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0)),
                    minimum_distances::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0));
                    update_other_fields::Bool = false) where {T <: SSDFloat}
    sim.electric_potential = refine(sim.electric_potential, max_diffs, minimum_distances)

    if update_other_fields
        fssrb::PotentialSimulationSetupRB{T, 3, 4, get_coordinate_system(sim.electric_potential.grid)} =
            PotentialSimulationSetupRB(sim.detector, sim.electric_potential.grid, sim.electric_potential.data)

        sim.ρ = ChargeDensity(ChargeDensityArray(fssrb), sim.electric_potential.grid)
        sim.ρ_fix = ChargeDensity(FixedChargeDensityArray(fssrb), sim.electric_potential.grid)
        sim.ϵ = DielectricDistribution(DielektrikumDistributionArray(fssrb), sim.electric_potential.grid)
        sim.point_types = PointTypes(PointTypeArray(fssrb), sim.electric_potential.grid)
    end
    nothing
end
"""
    function refine!(sim::Simulation{T}, ::Type{WeightingPotential}, max_diffs::Tuple{<:Real,<:Real,<:Real}, minimum_distances::Tuple{<:Real,<:Real,<:Real})

Takes the current state of `sim.weighting_potentials[contact_id]` and refines it with respect to the input arguments
`max_diffs` and `minimum_distances`.
"""
function refine!(sim::Simulation{T}, ::Type{WeightingPotential}, contact_id::Int,
                    max_diffs::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0)),
                    minimum_distances::Tuple{<:Real,<:Real,<:Real} = (T(0), T(0), T(0))) where {T <: SSDFloat}
    sim.weighting_potentials[contact_id] = refine(weighting_potentials[contact_id], max_diffs, minimum_distances)
    nothing
end

function _calculate_potential!( sim::Simulation{T}, potential_type::UnionAll, contact_id::Union{Missing, Int} = missing;
        init_grid_size::Union{Missing, NTuple{3, Int}} = missing,
        init_grid_spacing::Union{Missing, Tuple{<:Real,<:Real,<:Real}} = missing,
        grid::Union{Missing, Grid{T}} = missing,
        convergence_limit::Real = 1e-7,
        max_refinements::Int = 3,
        refinement_limits::Union{Missing, Tuple{<:Real,<:Real,<:Real}} = missing,
        min_grid_spacing::Union{Missing, Tuple{<:Real,<:Real,<:Real}} = missing,
        depletion_handling::Bool = false,
        use_nthreads::Int = Base.Threads.nthreads(),
        sor_consts::Union{Missing, <:Real, Tuple{<:Real,<:Real}} = missing,
        max_n_iterations::Int = 50000,
        verbose::Bool = true,
    )::Nothing where {T <: SSDFloat}

    begin # preperations
        convergence_limit::T = T(convergence_limit)
        isEP::Bool = potential_type == ElectricPotential
        isWP::Bool = !isEP
        if isWP depletion_handling = false end
        CS = get_coordinate_system(sim.detector)
        if ismissing(grid)
            grid = Grid(sim.detector, init_grid_size = init_grid_size, init_grid_spacing = init_grid_spacing, for_weighting_potential = isWP)
        end
        if ismissing(sor_consts)
            sor_consts = CS == :cylindrical ? (T(1.4), T(1.85)) : T(1.4)
        elseif length(sor_consts) == 1 && CS == :cylindrical
            sor_consts = (T(sor_consts), T(sor_consts))
        elseif length(sor_consts) > 1 && CS == :cartesian
            sor_consts = T(sor_consts[1])
        end
        if ismissing(refinement_limits)
            refinement_limits::NTuple{3, T} = CS == :cylindrical ? (T(1e-5), T(1e-5) / (0.25 * grid.axes[1][end]), T(1e-5)) : (T(1e-5), T(1e-5), T(1e-5))
        end
        if ismissing(min_grid_spacing)
            min_grid_spacing::NTuple{3, T} = CS == :cylindrical ? (T(1e-5), T(1e-5) / (0.25 * grid.axes[1][end]), T(1e-5)) : (T(1e-5), T(1e-5), T(1e-5))
        end
        n_iterations_between_checks::Int = div(10^7, length(grid))
        if n_iterations_between_checks < 20 n_iterations_between_checks = 20 end
        if use_nthreads > Base.Threads.nthreads()
            use_nthreads = Base.Threads.nthreads();
            @warn "`use_nthreads` was set to `1`. The environment variable `JULIA_NUM_THREADS` must be set appropriately before the julia session is started."
        end
        refine::Bool = max_refinements > 0 ? true : false
        only_2d::Bool = length(grid.axes[2]) == 1 ? true : false
        if CS == :cylindrical
            cyclic::T = grid.axes[2].interval.right - grid.axes[2].interval.left
            n_φ_sym::Int = only_2d ? 1 : Int(round(T(2π) / cyclic, digits = 3))
            n_φ_sym_info_txt = if only_2d
                "φ symmetry: Detector is φ-symmetric -> 2D computation."
            else
                "φ symmetry: calculating just 1/$(n_φ_sym) in φ of the detector."
            end
        end
        contact_potentials::Vector{T} = [contact.potential for contact in sim.detector.contacts]
        bias_voltage::T = (length(contact_potentials) > 0) ? (maximum(contact_potentials) - minimum(contact_potentials)) : T(0)
        if isWP bias_voltage = T(1) end
        if verbose
            println(
                "$(isEP ? "Electric" : "Weighting") Potential Calculation\n",
                if isEP "Bias voltage: $(bias_voltage) V\n" else "" end,
                if CS == :cylindrical "$n_φ_sym_info_txt\n" else "" end,
                "Precision: $T\n",
                "Convergence limit: $convergence_limit => $(round(abs(bias_voltage * convergence_limit), sigdigits=2)) V\n",
                "Threads: $use_nthreads\n",
                "Coordinate system: $(CS)\n",
                "Initial grid dimension: $(size(grid))\n",
                "Refine? -> $refine\n",
                "Refinement parameters:\n",
                "\tmaximum number of refinements: $(max_refinements)\n",
                "\tminimum grid spacing:\n",
                "\t\tr: $(min_grid_spacing[1]) m\n",
                "\t\tφ: $(min_grid_spacing[2]) rad\n",
                "\t\tz: $(min_grid_spacing[3]) m\n",
                "\tRefinement limits:\n",
                "\t\tr: $(refinement_limits[1]) -> $(round(abs(bias_voltage * refinement_limits[1]), sigdigits=2)) V\n",
                "\t\tφ: $(refinement_limits[2]) -> $(round(abs(bias_voltage * refinement_limits[2]), sigdigits=2)) V\n",
                "\t\tz: $(refinement_limits[3]) -> $(round(abs(bias_voltage * refinement_limits[3]), sigdigits=2)) V\n",
                ""
            )
        end
    end
    if isEP
        apply_initial_state!(sim, potential_type, grid)
        update_till_convergence!( sim, potential_type, convergence_limit,
                                  n_iterations_between_checks = n_iterations_between_checks,
                                  max_n_iterations = max_n_iterations,
                                  depletion_handling = depletion_handling,
                                  use_nthreads = use_nthreads,
                                  sor_consts = sor_consts )
    else
        apply_initial_state!(sim, potential_type, contact_id, grid)
        update_till_convergence!( sim, potential_type, contact_id, convergence_limit,
                                    n_iterations_between_checks = n_iterations_between_checks,
                                    max_n_iterations = max_n_iterations,
                                    depletion_handling = depletion_handling,
                                    use_nthreads = use_nthreads,
                                    sor_consts = sor_consts )
    end

    refinement_counter::Int = 1
    if refine
        max_diffs::NTuple{3, T} = abs.(refinement_limits .* bias_voltage)
        refine_at_inds::NTuple{3, Vector{Int}} = if isEP
            _get_refinement_inds(sim.electric_potential.data, sim.electric_potential.grid, max_diffs, min_grid_spacing)
        else
            _get_refinement_inds(sim.weighting_potentials[contact_id].data, sim.weighting_potentials[contact_id].grid, max_diffs, min_grid_spacing)
        end
        while any(!isempty, refine_at_inds[1]) && refinement_counter <= max_refinements
            if isEP
                sim.electric_potential = ElectricPotential(add_points_and_interpolate(
                        sim.electric_potential.data, sim.electric_potential.grid, refine_at_inds...)...)
                if verbose @info "New Grid Size = $(size(sim.electric_potential.grid))" end
            else
                sim.weighting_potentials[contact_id] = WeightingPotential(add_points_and_interpolate(
                    sim.weighting_potentials[contact_id].data, sim.weighting_potentials[contact_id].grid, refine_at_inds...)...)
                if verbose @info "New Grid Size = $(size(sim.weighting_potentials[contact_id].grid))" end
            end
            n_iterations_between_checks = div(10^7, length(isEP ? sim.electric_potential.grid : sim.weighting_potentials[contact_id].grid))
            if n_iterations_between_checks < 20 n_iterations_between_checks = 20 end

            if isEP
                update_till_convergence!( sim, potential_type, convergence_limit,
                    n_iterations_between_checks = n_iterations_between_checks,
                    max_n_iterations = max_n_iterations,
                    depletion_handling = depletion_handling,
                    use_nthreads = use_nthreads,
                    sor_consts = sor_consts )
            else
                update_till_convergence!( sim, potential_type, contact_id, convergence_limit,
                                    n_iterations_between_checks = n_iterations_between_checks,
                                    max_n_iterations = max_n_iterations,
                                    depletion_handling = depletion_handling,
                                    use_nthreads = use_nthreads,
                                    sor_consts = sor_consts )
            end

            refine_at_inds = if isEP
                _get_refinement_inds(sim.electric_potential.data, sim.electric_potential.grid, max_diffs, min_grid_spacing)
            else
                _get_refinement_inds(sim.weighting_potentials[contact_id].data, sim.weighting_potentials[contact_id].grid, max_diffs, min_grid_spacing)
            end
            refinement_counter += 1
        end
    end

    nothing
end

"""
    calculate_weighting_potential!(sim::Simulation{T}, contact_id::Int; kwargs...)::Nothing

Compute the weighting potential for the contact with id `contact_id`
for the given Simulation `sim` on an adaptive grid through successive over relaxation.

There are serveral `<keyword arguments>` which can be used to tune the computation:

# Keywords
- `convergence_limit::Real`: `convergence_limit` times the bias voltage sets the convergence limit of the relaxation.
    The convergence value is the absolute maximum difference of the potential between two iterations of all grid points.
    Default of `convergence_limit` is `2e-6` (times bias voltage).
- `max_refinements::Int`: Number of maximum refinements. Default is `2`. Set it to `0` to switch off refinement.
- `refinement_limits::Tuple{<:Real, <:Real, <:Real}`: Tuple of refinement limits for each dimension
    (in case of cylindrical coordinates the order is `r`, `φ`, `z`).
    A refinement limit (e.g. `refinement_limits[1]`) times the bias voltage of the detector `det` is the
    maximum allowed voltage difference between two neighbouring grid points in the respective dimension.
    When the difference is larger, new points are created inbetween. Default is `[1e-5, 1e-5, 1e-5]`.
- `init_grid_spacing::Tuple{<:Real, <:Real, <:Real}`: Tuple of the initial distances between two grid points for each dimension.
    For normal coordinates the unit is meter. For angular coordinates, the unit is radiance.
    It prevents the refinement to make the grid to fine. Default is `[0.005, 10.0, 0.005]``.
- `min_grid_spacing::Tuple{<:Real, <:Real, <:Real}`: Tuple of the mimimum allowed distance between two grid points for each dimension.
    For normal coordinates the unit is meter. For angular coordinates, the unit is radiance.
    It prevents the refinement to make the grid to fine. Default is [`1e-6`, `1e-6`, `1e-6`].
- `grid::Grid{T, N, S}`: Initial grid used to start the simulation. Default is `Grid(detector, init_grid_spacing=init_grid_spacing)`.
- `depletion_handling::Bool`: Enables the handling of undepleted regions. Default is false.
- `use_nthreads::Int`: Number of threads to use in the computation. Default is `Base.Threads.nthreads()`.
    The environment variable `JULIA_NUM_THREADS` must be set appropriately before the Julia session was
    started (e.g. `export JULIA_NUM_THREADS=8` in case of bash).
- `sor_consts::Union{<:Real, NTuple{2, <:Real}}`: Two element tuple in case of cylindrical coordinates.
    First element contains the SOR constant for `r` = 0.
    Second contains the constant at the outer most grid point in `r`. A linear scaling is applied in between.
    First element should be smaller than the second one and both should be ∈ [1.0, 2.0]. Default is [1.4, 1.85].
    In case of cartesian coordinates only one value is taken.
- `max_n_iterations::Int`: Set the maximum number of iterations which are performed after each grid refinement.
    Default is `10000`. If set to `-1` there will be no limit.
- `verbose::Bool=true`: Boolean whether info output is produced or not.
"""
function calculate_weighting_potential!(sim::Simulation{T}, contact_id::Int, args...; n_points_in_φ::Union{Missing, Int} = missing, kwargs...)::Nothing where {T <: SSDFloat}
    # S = get_coordinate_system(sim.detector)
    # periodicity::T = get_periodicity(sim.detector.world.intervals[2])
    # if S == :cylindrical && periodicity == T(0)
    #     if ismissing(n_points_in_φ)
    #         @info "\tIn weighing potential calculation: Keyword `n_points_in_φ` not set.\n\t\tDefault is `n_points_in_φ = 36`. 2D field will be extended to 36 points in φ."
    #         n_points_in_φ = 36
    #     else
    #         if !(n_points_in_φ > 1 && iseven(n_points_in_φ))
    #             @info "\tIn weighing potential calculation: Keyword `n_points_in_φ` is $(n_points_in_φ) but must be even and larger than 1.\n\t\t`n_points_in_φ` is now set to 36. 2D field will be extended to 36 points in φ."
    #             n_points_in_φ = 36
    #         end
    #     end
    # end
    # wps = calculate_weighting_potential(sim.detector, contact_id, args...; kwargs...)
    # if S == :cylindrical && size(wps.potential, 2) == 1 && !ismissing(n_points_in_φ)
    #     wp = WeightingPotential(wps, n_points_in_φ = n_points_in_φ)
    # else
    #     wp = WeightingPotential(wps)
    # end
    _calculate_potential!(sim, WeightingPotential, contact_id, args...; kwargs...)
    nothing
end

"""
    calculate_electric_potential!(sim::Simulation{T}; kwargs...)::Nothing

Compute the electric potential for the given Simulation `sim` on an adaptive grid
through successive over relaxation.

There are serveral `<keyword arguments>` which can be used to tune the computation:

# Keywords
- `convergence_limit::Real`: `convergence_limit` times the bias voltage sets the convergence limit of the relaxation.
    The convergence value is the absolute maximum difference of the potential between two iterations of all grid points.
    Default of `convergence_limit` is `2e-6` (times bias voltage).
- `max_refinements::Int`: Number of maximum refinements. Default is `2`. Set it to `0` to switch off refinement.
- `refinement_limits::Tuple{<:Real, <:Real, <:Real}`: Tuple of refinement limits for each dimension
    (in case of cylindrical coordinates the order is `r`, `φ`, `z`).
    A refinement limit (e.g. `refinement_limits[1]`) times the bias voltage of the detector `det` is the
    maximum allowed voltage difference between two neighbouring grid points in the respective dimension.
    When the difference is larger, new points are created inbetween. Default is `[1e-5, 1e-5, 1e-5]`.
- `init_grid_spacing::Tuple{<:Real, <:Real, <:Real}`: Tuple of the initial distances between two grid points for each dimension.
    For normal coordinates the unit is meter. For angular coordinates, the unit is radiance.
    It prevents the refinement to make the grid to fine. Default is `[0.005, 10.0, 0.005]``.
- `min_grid_spacing::Tuple{<:Real, <:Real, <:Real}`: Tuple of the mimimum allowed distance between two grid points for each dimension.
    For normal coordinates the unit is meter. For angular coordinates, the unit is radiance.
    It prevents the refinement to make the grid to fine. Default is [`1e-6`, `1e-6`, `1e-6`].
- `grid::Grid{T, N, S}`: Initial grid used to start the simulation. Default is `Grid(detector, init_grid_spacing=init_grid_spacing)`.
- `depletion_handling::Bool`: Enables the handling of undepleted regions. Default is false.
- `use_nthreads::Int`: Number of threads to use in the computation. Default is `Base.Threads.nthreads()`.
    The environment variable `JULIA_NUM_THREADS` must be set appropriately before the Julia session was
    started (e.g. `export JULIA_NUM_THREADS=8` in case of bash).
- `sor_consts::Union{<:Real, NTuple{2, <:Real}}`: Two element tuple in case of cylindrical coordinates.
    First element contains the SOR constant for `r` = 0.
    Second contains the constant at the outer most grid point in `r`. A linear scaling is applied in between.
    First element should be smaller than the second one and both should be ∈ [1.0, 2.0]. Default is [1.4, 1.85].
    In case of cartesian coordinates only one value is taken.
- `max_n_iterations::Int`: Set the maximum number of iterations which are performed after each grid refinement.
    Default is `10000`. If set to `-1` there will be no limit.
- `verbose::Bool=true`: Boolean whether info output is produced or not.
"""
function calculate_electric_potential!(sim::Simulation{T}, args...; kwargs...)::Nothing where {T <: SSDFloat}
    _calculate_potential!(sim, ElectricPotential, args...; kwargs...)
    nothing
end

"""
    calculate_electric_field!(sim::Simulation{T}, args...; n_points_in_φ::Union{Missing, Int} = missing, kwargs...)::Nothing

ToDo...
"""
function calculate_electric_field!(sim::Simulation{T}, args...; n_points_in_φ::Union{Missing, Int} = missing, kwargs...)::Nothing where {T <: SSDFloat}
    S = get_coordinate_system(sim.detector)
    periodicity::T = get_periodicity(sim.detector.world.intervals[2])
    e_pot, point_types = if S == :cylindrical && periodicity == T(0) # 2D, only one point in φ
        if ismissing(n_points_in_φ)
            @info "\tIn electric field calculation: Keyword `n_points_in_φ` not set.\n\t\tDefault is `n_points_in_φ = 36`. 2D field will be extended to 36 points in φ."
            n_points_in_φ = 36
        else
            if !(n_points_in_φ > 1 && iseven(n_points_in_φ))
                @info "\tIn electric field calculation: Keyword `n_points_in_φ` is $(n_points_in_φ) but must be even and larger than 1.\n\t\t`n_points_in_φ` is now set to 36. 2D field will be extended to 36 points in φ."
                n_points_in_φ = 36
            end
        end
        get_2π_potential(sim.electric_potential, n_points_in_φ = n_points_in_φ),
        get_2π_potential(sim.point_types,  n_points_in_φ = n_points_in_φ);
    elseif S == :cylindrical
        get_2π_potential(sim.electric_potential),
        get_2π_potential(sim.point_types)
    else
        sim.electric_potential,
        sim.point_types
    end
    sim.electric_field = get_electric_field_from_potential(e_pot, point_types);
    nothing
end

function set_charge_drift_model!(sim::Simulation{T}, charge_drift_model::AbstractChargeDriftModel{T})::Nothing where {T <: SSDFloat}
    sim.charge_drift_model = charge_drift_model
    nothing
end

function calculate_drift_fields!(sim::Simulation{T};
    use_nthreads::Int = Base.Threads.nthreads())::Nothing where {T <: SSDFloat}
    sim.electron_drift_field = ElectricField(get_electron_drift_field(sim.electric_field.data, sim.charge_drift_model, use_nthreads = use_nthreads), sim.electric_field.grid)
    sim.hole_drift_field = ElectricField(get_hole_drift_field(sim.electric_field.data, sim.charge_drift_model, use_nthreads = use_nthreads), sim.electric_field.grid)
    nothing
end
@deprecate apply_charge_drift_model!(args...; kwargs...) calculate_drift_fields!(args...; kwargs...)

function get_interpolated_drift_field(ef::ElectricField)
    get_interpolated_drift_field(ef.data, ef.grid)
end

function drift_charges( sim::Simulation{T}, starting_positions::Vector{CartesianPoint{T}};
                        Δt::RealQuantity = 5u"ns", max_nsteps::Int = 1000, verbose::Bool = true )::Vector{EHDriftPath{T}} where {T <: SSDFloat}
    return _drift_charges(   sim.detector, sim.point_types.grid, sim.point_types, starting_positions,
                             get_interpolated_drift_field(sim.electron_drift_field), get_interpolated_drift_field(sim.hole_drift_field),
                             Δt, max_nsteps = max_nsteps, verbose = verbose)::Vector{EHDriftPath{T}}
end

function get_signal(sim::Simulation{T}, drift_paths::Vector{EHDriftPath{T}}, energy_depositions::Vector{T}, contact_id::Int; Δt::TT = T(5) * u"ns") where {T <: SSDFloat, TT}
    dt::T = to_internal_units(internal_time_unit, Δt)
    wp::Interpolations.Extrapolation{T, 3} = interpolated_scalarfield(sim.weighting_potentials[contact_id])
    timestamps = _common_timestamps( drift_paths, dt )
    signal::Vector{T} = zeros(T, length(timestamps))
    add_signal!(signal, timestamps, drift_paths, energy_depositions, wp, Val(get_coordinate_system(sim.detector)))
    return RDWaveform( range(zero(T) * unit(Δt), step = T(ustrip(Δt)) * unit(Δt), length = length(signal)), signal )
end

"""
    function simulate!(sim::Simulation{T};  max_refinements::Int = 1, verbose::Bool = false,
                                        depletion_handling::Bool = false, convergence_limit::Real = 1e-5 ) where {T <: SSDFloat}

ToDo...
"""
function simulate!(sim::Simulation{T};  max_refinements::Int = 1, verbose::Bool = false,
                                        depletion_handling::Bool = false, convergence_limit::Real = 1e-7 ) where {T <: SSDFloat}
    calculate_electric_potential!(  sim,
                                    max_refinements = max_refinements,
                                    verbose = verbose,
                                    init_grid_size = (10,10,10),
                                    depletion_handling = depletion_handling,
                                    convergence_limit = convergence_limit )
    for contact in sim.detector.contacts
        calculate_weighting_potential!(sim, contact.id, max_refinements = max_refinements,
                init_grid_size = (10,10,10),
                verbose = verbose, convergence_limit = convergence_limit)
    end
    calculate_electric_field!(sim)
    set_charge_drift_model!(sim, sim.charge_drift_model)
    calculate_drift_fields!(sim)
    @info "Detector simulation done"
end

function _get_abs_bias_voltage(det::SolidStateDetector{T}) where {T <: SSDFloat}
    potentials::Vector{T} = map(c -> c.potential, det.contacts)
    return (maximum(potentials) - minimum(potentials)) * u"V"
end

calculate_stored_energy(sim::Simulation) =
    calculate_stored_energy(sim.electric_field, sim.ϵ)

function calculate_stored_energy(ef::ElectricField{T,3,S}, ϵ::DielectricDistribution{T,3,S}) where {T <: SSDFloat, S}
    W::T = 0

    cylindric::Bool = S == :cylindrical
    cartesian::Bool = !cylindric
    r0_handling::Bool = typeof(ef.grid.axes[1]).parameters[2] == :r0

    ax1::Vector{T} = collect(ef.grid[1])
    ax2::Vector{T} = collect(ef.grid[2])
    ax3::Vector{T} = collect(ef.grid[3])
    mp1::Vector{T} = midpoints(get_extended_ticks(ef.grid[1]))
    mp2::Vector{T} = midpoints(get_extended_ticks(ef.grid[2]))
    mp3::Vector{T} = midpoints(get_extended_ticks(ef.grid[3]))
    Δax1::Vector{T} = diff(ax1)
    Δax2::Vector{T} = diff(ax2)
    Δax3::Vector{T} = diff(ax3)
    Δmp1::Vector{T} = diff(mp1)
    Δmp2::Vector{T} = diff(mp2)
    Δmp3::Vector{T} = diff(mp3)

    w1r::Vector{T} = inv.(Δmp1) .* (mp1[2:end] .- ax1)
    w1l::Vector{T} = inv.(Δmp1) .* (ax1 - mp1[1:end-1])
    w2r::Vector{T} = inv.(Δmp2) .* (mp2[2:end] .- ax2)
    w2l::Vector{T} = inv.(Δmp2) .* (ax2 - mp2[1:end-1])
    w3r::Vector{T} = inv.(Δmp3) .* (mp3[2:end] .- ax3)
    w3l::Vector{T} = inv.(Δmp3) .* (ax3 - mp3[1:end-1])

    if cylindric
        mp1[1] = 0
        mp1[end] = ax1[end]
        Δmp1 = ((mp1[2:end].^2) .- (mp1[1:end-1].^2)) ./ 2
    end
    V::T = 0
    for i3 in 1:size(ϵ, 3)-1
        _Δmp3::T = Δmp3[i3]
        if (i3 == 1 || i3 == size(ϵ, 3)-1) _Δmp3 /= 2 end
        for i2 in 1:size(ϵ, 2)-1
            _Δmp2::T = Δmp2[i2]
            if (cartesian && (i2 == 1 || i2 == size(ϵ, 2)-1)) _Δmp2 /= 2 end
            for i1 in 1:size(ϵ, 1)-1
                _Δmp1::T = Δmp1[i1]
                if (cartesian && (i1 == 1 || i1 == size(ϵ, 1)-1)) _Δmp1 /= 2 end
                ev::SArray{Tuple{3},Float32,1,3} = ef.data[i1, i2, i3]
                dV::T = _Δmp3 * _Δmp2 * _Δmp1
                _ϵ::T = sum([
                    ϵ[i1, i2, i3]             * w1l[i1] * w2l[i2] * w3l[i3],
                    ϵ[i1 + 1, i2, i3]         * w1r[i1] * w2l[i2] * w3l[i3],
                    ϵ[i1, i2 + 1, i3]         * w1l[i1] * w2r[i2] * w3l[i3],
                    ϵ[i1, i2, i3 + 1]         * w1l[i1] * w2l[i2] * w3r[i3],
                    ϵ[i1 + 1, i2 + 1, i3]     * w1r[i1] * w2r[i2] * w3l[i3],
                    ϵ[i1 + 1, i2, i3 + 1]     * w1r[i1] * w2l[i2] * w3r[i3],
                    ϵ[i1, i2 + 1, i3 + 1]     * w1l[i1] * w2r[i2] * w3r[i3],
                    ϵ[i1 + 1, i2 + 1, i3 + 1] * w1r[i1] * w2r[i2] * w3r[i3]
                ])
                V += dV
                W += sum(ev.^2) * dV * _ϵ
            end
        end
    end
    E =  W * ϵ0 / 2 * u"J"
    if cylindric && (size(ϵ, 2) - 1 != size(ef, 2))
        E *= size(ef, 2) / (size(ϵ, 2) - 1)
    end
    return E
end

export calculate_capacitance
"""
    calculate_capacitance(sim::Simulation{T})::T where {T <: SSDFloat}

Calculates the capacitance of an detector in Farad.
"""
function calculate_capacitance(sim::Simulation{T}) where {T <: SSDFloat}
    W = calculate_stored_energy(sim)
    return uconvert(u"pF", 2 * W / (_get_abs_bias_voltage(sim.detector)^2))
end

ext: jl


text:
#
# XPA.jl --
#
# Implement XPA communication via the dynamic library.
#
#------------------------------------------------------------------------------
#
# This file is part of XPA.jl released under the MIT "expat" license.
# Copyright (C) 2016-2019, Éric Thiébaut (https://github.com/emmt/XPA.jl).
#

module XPA

export
    XPA_VERSION

using Base: ENV, @propagate_inbounds

isfile(joinpath(@__DIR__, "..", "deps", "deps.jl")) ||
    error("XPA not properly installed.  Please run Pkg.build(\"XPA\")")
include(joinpath("..", "deps", "deps.jl"))
include("types.jl")
include("misc.jl")
include("client.jl")
include("server.jl")

end # module


text:
module Propensity

# Dependencies
using Statistics
using StatsModels
using GLM
using DataFrames
using Distributions
using Gadfly
using StatsBase:sample

# Module test
greet() = print("Hello World!")

# Core functionality (refactor)

function assign_formula(y, df)

    y = Symbol(y)
    fm = term(y) ~ sum(term.(Symbol.(names(df, Not(y)))))
    return fm
end

# split to minority and majority class
# record size of minority class (n)
function split_by_class(df, target)
    target = Symbol(target)
    df_1 = df[df[target] .== 1, :]
    df_0 = df[df[target] .== 0, :]
    min_cls_size = size(df_1)[1]
    return df_1, df_0, min_cls_size
end

# Generate n trainings sets
# concat minority and n subsampled majorities
function make_train_sets(df_0, df_1, n_sets)

    min_cls_size = size(df_1)[1]

    samples = []

    for i in 1:n_sets
        subsample = sample(1:nrow(df_0),
            min_cls_size,
            replace=false)

        df_0_sampled = df_0[subsample, :]
        df_train = [df_0_sampled; df_1]

        samples = vcat(samples, df_train)
    end
    return samples
end

# train logit functions for each n train set
function train_n_models(y, df_list)

    models = []
    for df in df_list
        logit = fit_logit(y, df)
        models = append!(models, [logit])
    end
    return models
end

function fit_logit(y, df)

    fm = assign_formula(y, df)

    logit = glm(fm, df, Binomial(), LogitLink())
    return logit
end

function assign_propensity_scores(df, logit)

    df[!,Symbol("Propensity_Scores")] = predict(logit, df)
    return df
end

# Matching functionality (refactor)

function bin_quartile(x, Q1, Q2, Q3, Q4)
    if x .<= Q1
        return "Q1"
    elseif x .<= Q2
        return "Q2"
    elseif x .<= Q3
        return "Q3"
    else
        return "Q4"
    end
end

function assign_quartile(df, col, new_col)

    df = sort(df, [order(Symbol(col), rev=false)])

    Q1, Q2, Q3, Q4 = quantile(df[!, Symbol(col)], [0.25 0.5 0.75 1.0])

    df[!, Symbol(new_col)] = bin_quartile.(df[:, Symbol(col)], Q1, Q2, Q3, Q4)

    return df
end

function sample_by_covariate_quartile()
end

function sample_by_propensity_score()
end

# Visualize/Analysis (refactor)

function plot_prop_by_factor(df, factor)
    set_default_plot_size(12cm, 8cm)

    fig = plot(
        df,
        x=Symbol(factor),
        y=:Propensity_Scores,
        color=Symbol(factor),
        Stat.x_jitter(range=0.5),
        Guide.title("Propensity Score")
    )
    return fig
end

function plot_prop_by_covariate(
        df,
        factor,
        covariate,
        covariate_quartiles
    )

    set_default_plot_size(21cm, 8cm)

    factor = Symbol(factor)
    covariate = Symbol(covariate)
    covariate_quartiles = Symbol(covariate_quartiles)

    covariate_by_prop = plot(
    df,
        x=covariate,
        y=:Propensity_Scores,
    color=factor,
        Stat.x_jitter(range=0.5),
    Guide.title("$covariate by Propensity Score")
    )
    covariate_quartile_by_prop = plot(
    df,
        x=covariate_quartiles,
        y=:Propensity_Scores,
    color=factor,
        Stat.x_jitter(range=0.5),
    Guide.title("$covariate Quartile by Propensity Score")
    )

    h = hstack(covariate_by_prop, covariate_quartile_by_prop)

    return h
end

end # module

max_stars_repo_name: pkmklong/propensity

max_issues_repo_name: pkmklong/propensity

max_forks_repo_name: pkmklong/propensity


text:
@testset "Recorded tests" begin
    include("recorded-tests/basic-read.jl")
end


text:
#
## Nil methods
#

function excl(::Sidef_Types_Nil_Nil)
    TRUE
end

function eq(::Sidef_Types_Nil_Nil, ::Sidef_Types_Nil_Nil)
    TRUE
end

function eq(::Sidef_Types_Nil_Nil, ::Sidef_Object)
    FALSE
end

function ne(::Sidef_Types_Nil_Nil, ::Sidef_Types_Nil_Nil)
    FALSE
end

function ne(::Sidef_Types_Nil_Nil, ::Sidef_Object)
    TRUE
end

max_stars_repo_name: trizen/Sidef-Deparse-Julia

max_issues_repo_name: trizen/Sidef-Deparse-Julia

max_forks_repo_name: trizen/Sidef-Deparse-Julia


text:
# iterative_methods.jl : convenience wrappers for smoothers in the IterativeSolvers package

abstract type Smoother end

struct Jacobi <: Smoother end

struct GaussSeidel <: Smoother end

struct RedBlackGaussSeidel <: Smoother end

smooth!(x::AbstractVector, A::AbstractMatrix, b::AbstractVector, n::Int, ::Jacobi) = jacobi!(x,A,b,maxiter=n)

smooth!(x::AbstractVector, A::AbstractMatrix, b::AbstractVector, n::Int, ::GaussSeidel) = gauss_seidel!(x,A,b,maxiter=n)


text:
#  type-based valuations

zero{G<:Grasp, R<:Real}(::Type{G},::Type{R}) = (G)(zero(R), zero(R))
 one{G<:Grasp, R<:Real}(::Type{G},::Type{R}) = (G)( one(R),  one(R))

# Inf{G<:Grasp, R<:Real}(::Type{G},::Type{R}) = (G)((R)(Inf), (R)(Inf))
# NaN{G<:Grasp, R<:Real}(::Type{G},::Type{R}) = (G)((R)(NaN), (R)(NaN))


text:
#=
Histogram plot for the paper.
left subplot: histograms of tapping, uwe, PGmusic
right: histogram of 1 playalong, with colors denoting instruments
=#

using DrWatson
quickactivate(@__DIR__, "VelocityMTD")
using MusicManipulations, PyPlot, StatsBase, Statistics
include(srcdir("extract_data.jl"))
include(srcdir("vel_normalization.jl"))
include(srcdir("style.jl"))

# Get velocities/drumnotes
datano = [1,1,1,7]
mtds, velstap = tapping_data()
tapping = velstap[datano[1]]
mtds, velsuwe = uwe_data(; subset="melody")
uwe = velsuwe[datano[2]]
mtds, velspg = pgmusic_data(;subset="melody")
pgmusic = velspg[datano[3]]

datanames = ["tapping",  "Uwe Meile", "PG music", "drums playalongs"]

# drums notes
drums = begin
    pdir = datadir("exp_raw", "playalongs")
    file = readdir(pdir)[datano[4]]
    grid = file[3] == 'j' ? (0:1//3:1) : (0:1//4:1)
    midi = readMIDIFile(joinpath(pdir, file))
    drums = Drums.getnotes_td50(midi.tracks[2], midi.tpq)
    drums = drums_cleanup(drums)
    d = estimate_delay_recursive(drums, grid, 5)
    drums = translate(drums, -d)
end

fig = figure(figsize = (2figx, 2figx/3))
axp = subplot(121)
axd = subplot(122)
fig.tight_layout()

# %% plot smoothened piano histograms
axp.clear()
for (i, v) in enumerate((tapping, uwe, pgmusic))
    sh = smoothhist(velhist(v), 3)
    sh ./= sum(sh)
    x = 0:length(sh)-1
    axp.plot(x, sh, color = "C$(i)", label = datanames[i])
    axp.fill_between(x, sh, color = "C$(i)",
    alpha=0.25)
end
axp.legend()
axp.set_yticks([])
axp.set_ylabel("normalized density")
axp.set_xlabel("velocity")

# %%
axd.clear()
channeled = channel_direct(drums)
sepa = separatepitches(channeled)

for (key, val) in sepa
    key == INVALID && continue
    sh = smoothhist(val, 3)
    x = 0:length(sh)-1
    axd.plot(x, sh, color = "C$(key-1)", label = CHANNELNAMES[key])
    axd.fill_between(x, sh, color = "C$(key-1)",
    alpha=0.25)
end
axd.set_yticks([])
axd.set_xlim(0, 150)
axd.legend()
axd.set_ylabel("density")
axd.set_xlabel("velocity")

# %% final adjust
fig.tight_layout()
add_identifiers!(fig)
# fig.savefig(papersdir("figs", "histograms"))

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
using Base.Test
using EGR

println("TestVerify")

for (gradientOracle, numVars, numTP, csDataType, LossFunctionString, myOutputter, L2reg, thisDataName, thisProblem)  in Oracles

	srand(1)

	println(" $thisDataName $LossFunctionString L2reg = $L2reg")
	tol = 1e-6

	xs = Any[zeros(numVars,1),2*rand(numVars,1)-1, randn(numVars,1)]
	xtexts = ["zeros(numVars,1)","2*rand(numVars,1)-1", "randn(numVars,1)"]
	for i in 1:3
		println("  For x = $(xtexts[i])")
		x=xs[i]
		# println("  x = $(x)")
		est=zeros(numVars)
		for i=1:numTP
			(f,g1)= gradientOracle(x,i)
			est+=g1
		end
		g=est/numTP
		# println("  g = $(g)")
		gDiff = zeros(numVars)
		passedDiff = false
		for displacementAccuracy in 1:15
			for j in 1:numVars
				a=zeros(numVars,1)
				b=zeros(numVars,1)
				a[j] = 10.0^(-displacementAccuracy)
				b[j] = -10.0^(-displacementAccuracy)
				a+=x
				b+=x
				fa= gradientOracle(a)
				fb= gradientOracle(b)
				gDiff[j] = (fa-fb)/(2*10.0^(-displacementAccuracy))
			end

			# println("  gDiff = $(gDiff)")
			relError = norm(gDiff-g)/norm(g)
			println("   For displacementAccuracy = $displacementAccuracy, relError = $relError")
			if relError<=tol
				passedDiff = true
				break
			end
		end
		if passedDiff == false
			error("F values do not work with gradient!")
		end

	end
end

println("TestVerify successful!")
println()


text:
"""
    Ke_VL_model.jl

Written by Carl Whitfield, University of Manchester, 2021

Only to be included by definitions.jl
"""

#===================== VL trajectory params  =====================#

#====> Ke et al. 2021 (fitted) <====#
#log means of (1) peak viral load (copies/ml), (2) peak time (days from infection)
#(3) 1/r -- inverse growth rate (days) (4) 1/d -- inverse decay rate (days)
const Kmvμ = [17.546426, 1.392914, -1.197782, -0.671588]
#log covariance matrix
const KmvΣ = [0.909403 0.043808 -0.010543 0.080628;
              0.043808 0.032527 0.028866 0.004503;
             -0.010543 0.028866 0.028640 0.000383;
              0.080628 0.004503 0.000383 0.026590]
const Log_VLparams_dist = MvNormal(Kmvμ, KmvΣ)  #Distribution used to generate params

#=================================================================#

#=================== Infectivity vs. VL params ===================#

#====> Ke et al. 2021 (fitted) <====#
#log means of parameters J and h (Inf = J VL^h / (VL^h + Km^h))
const Kmvinfμ = [1.782, -0.108]
#log covars
const KmvinfΣ = [0.3387 0.0265; 0.0265 0.1270]
const Jscale = 4.188/exp(Kmvinfμ[1] + KmvinfΣ[1,1]/2)  #normalisation of mean peak inf (a.u.)
const Kinf_Km = 4.0e6
const Log_infparams_dist = MvNormal(Kmvinfμ, KmvinfΣ)  #Distribution used to generate params

# KmvΣ = reshape(KmvΣ,(4,4))
function generate_ke_VL_params()
    Ps = exp.(rand(MvNormal(Kmvμ, KmvΣ)))
    return log10(Ps[1]), Ps[2], 1/Ps[3], 1/Ps[4]
end

function generate_ke_inf_params()
    return exp.(rand(MvNormal(Kmvinfμ, KmvinfΣ)))
end

function infectivity_Ke(PVL::Float64, r::Float64, d::Float64, tp::Float64, J::Float64, h::Float64, T::Int64)
    #if viral load is piecewise linear, and infectivity is hill function, what is cumulative infectivity?
    m_up = r/log(10)
    m_down = d/log(10)

    cum_inf = zeros(T+1)
    t_inf = collect(-1:(T-1)) .+ 0.5
    t_inf[1] = 0
    t_inf[T+1] = T - 1.0

    V0 = (PVL - m_up*tp)
    K_m = 4.0e6^h
    a = h*r
    PInf = Jscale*J
    c1 = K_m * 10^(-V0*h)
    cond1 = (t_inf .> 0) .* (t_inf .<=  tp)
    cum_inf[cond1] = PInf .* (log.(exp.(a .* t_inf[cond1]) .+ c1)
                              .- log(c1 + 1.0)) ./ a
    ci_peak = PInf .* (log.(exp.(a * (tp)) .+ c1) .- log(c1 + 1)) ./ a
    cond2 = (t_inf .>  tp)
    b = -h*d
    c2 = K_m * 10^(-PVL*h)
    cum_inf[cond2] = ci_peak .+ PInf .* (log.(exp.(b .* (t_inf[cond2] .- tp)) .+ c2)
                              .- log(c2 + 1.0)) ./ b

    inf = (cum_inf[2:(T+1)] .- cum_inf[1:T]) ./ (t_inf[2:(T+1)] .- t_inf[1:T])

    return inf
end


text:
module DiscreteModelsTests

using Test
using Gridap
using Gridap.DiscreteModels: DiscreteModelFromData
using Gridap.CellValuesGallery
using Gridap.Grids: GridFromData
using Gridap.GridGraphs: FullGridGraphFromData

c2v = [[1,2,4,5],[2,3,5,6]]
v2c = [[1,],[1,2],[2,],[1,],[1,2],[2,]]
c2e = [[1,2,3,4],[5,6,4,7]]
e2c = [[1,],[1,],[1,],[1,2],[2,],[2,],[2,]]
e2v = [[1,2],[4,5],[1,4],[2,5],[2,3],[5,6],[3,7]]
v2e = [[1,3],[1,4,5],[5,7],[2,3],[2,4,6],[6,7]]
v2v = [[1,],[2,],[3,],[4,],[5,],[6,]]
e2e = [[1,],[2,],[3,],[4,],[5,],[6,],[7,]]
c2c = [[1,],[2,]]
v2x = Point{2,Float64}[(0.,0.),(0.,1.),(0.,2.),(1.,0.),(1.,1.),(1.,2.)]

c2v = CellValueFromArray(c2v)
v2c = CellValueFromArray(v2c)
c2e = CellValueFromArray(c2e)
e2c = CellValueFromArray(e2c)
e2v = CellValueFromArray(e2v)
v2e = CellValueFromArray(v2e)
v2v = CellValueFromArray(v2v)
e2e = CellValueFromArray(e2e)
c2c = CellValueFromArray(c2c)
v2x = CellValueFromArray(v2x)

data = Matrix{IndexCellVector}(undef,(3,3))
data[1,:] = [v2v,v2e,v2c]
data[2,:] = [e2v,e2e,e2c]
data[3,:] = [c2v,c2e,c2c]
graph = FullGridGraphFromData(data)

order = 1

nc = length(c2v)
t2 = (HEX_AXIS,HEX_AXIS)
c2t = ConstantCellValue(t2,nc)
c2o = ConstantCellValue(order,nc)
grid2 = GridFromData(v2x,c2v,c2t,c2o)

ne = length(e2v)
t1 = (HEX_AXIS,)
e2t = ConstantCellValue(t1,ne)
e2o = ConstantCellValue(order,ne)
grid1 = GridFromData(v2x,e2v,e2t,e2o)

nv = length(v2v)
t0 = ()
v2t = ConstantCellValue(t0,nv)
v2o = ConstantCellValue(order,nv)
grid0 = GridFromData(v2x,v2v,v2t,v2o)

v2l = CellValueFromArray([1,1,1,2,2,1])
e2l = CellValueFromArray([2,2,3,3,2,2,1])
c2l = CellValueFromArray([5,5])

dim_to_labels = [v2l, e2l, c2l]
tag1 = [1,3]
tag2 = [5,3,2]
tags = [tag1,tag2]
tag_to_name = ["tag1","tag2"]

facelabels = FaceLabels(dim_to_labels,tags,tag_to_name)

grids = Grid[grid0,grid1,grid2]

model = DiscreteModelFromData(grids,graph,facelabels)

@test string(model) == "DiscreteModelFromData object"

s = "DiscreteModelFromData object:\n celldim: 2\n 0-faces: 6\n 1-faces: 7\n 2-faces: 2\n tags: 2"

@test sprint(show,"text/plain",model) == s

test_discrete_model(model,2)

vgrid = Grid(model)

end # module


text:
using HTTP, JSON
using MAT
using CSV
using DataFrames
import Dates

% # ***************************************************************************
% #   ---------------------------------
% #   Written by: Vardhan Agnihotri
% #   - for -
% #   Mints: Multi-scale Integrated Sensing and Simulation
% #   ---------------------------------
% #   Date: July 8th, 2021
% #   ---------------------------------
% #   This module is written for generic implimentation of MINTS projects
% #   ---------------------------------------------------------------------
% #   https://github.com/mi3nts
% #
% #   Contact:
% #      email: 24agnihotriv@smtexas.org
% # ***********************************************************************
% #   DataFrames concatenation + EPA data analysis w/ modifiable function

function getDataFrame(bdate, edate, param, state, county, site)
    url = "https://aqs.epa.gov/data/api/sampleData/bySite?email=va648314@gmail.com&key=bayhare31&param=$param&bdate=$bdate&edate=$edate&state=$state&county=$county&site=$site"
    resp = HTTP.get(url)
    str = String(resp.body)
    jobj = JSON.parse(str)

    dateTimeList = []
    latList = []
    longList = []
    measureEPAList = []

    dfUTD = DataFrame(CSV.File("C:/Users/usr/PycharmProjects/mintsML/files/UTD2_5Data.csv"))
    measureUTDList = dfUTD.OPCN2_pm2_5

    for i in 1:length(jobj["Data"])
        dateStr = jobj["Data"][i]["date_local"]
        time = jobj["Data"][i]["time_local"]
        latitude = jobj["Data"][i]["latitude"]
        longitude = jobj["Data"][i]["longitude"]
        measurement = jobj["Data"][i]["sample_measurement"]

    dateTime = dateStr * " " * time

        if measurement == nothing
            continue
        else
            push!(dateTimeList, dateTime)
            push!(latList, latitude)
            push!(longList, longitude)
            push!(measureEPAList, measurement)
        end
    end

    reverseDateTime = reverse(dateTimeList)
    reverseLat = reverse(latList)
    reverseLong = reverse(longList)
    reverseEPAMeasure = reverse(measureEPAList)
    reverseUTDMeasure = reverse(measureUTDList)

    df = DataFrame(dateTime = reverseDateTime, latitude = reverseLat, longitude = reverseLong, measurementEPA = reverseEPAMeasure)
    println(df)
    # CSV.write("files/EPAdata.csv", df)

end

getDataFrame(20200101, 20201231, 88101, 48, 439, 1006)


text:
# -*- coding: utf-8 -*-
# + {}
module MonomialRegressionModels
import Base: show

using SemanticModels
using SemanticModels.ExprModels
using SemanticModels.ExprModels.Parsers
using SemanticModels.ExprModels.Transformations
import SemanticModels: model
import SemanticModels.ExprModels: isexpr

export MonomialRegressionModel, model

struct MonomialRegressionModel <: AbstractModel
    expr
    f
    coefficient
    objective
    interval
end

eval(m::MonomialRegressionModel) = eval(m.expr)

function model(::Type{MonomialRegressionModel}, ex::Expr)
    if ex.head == :block
        return model(MonomialRegressionModel, ex.args[2])
    end

    objective = callsites(ex, :optimize)[end].args[2]
    f = filter(isexpr, findfunc(ex, :f))[1]
    interval = findassign(ex, :a₀)[1]
    ldef = filter(isexpr, findfunc(ex, objective))
    coeff = funcarg(ldef[1])
    return MonomialRegressionModel(ex,
                              f,
                              coeff,
                              objective,
                              interval)
end

function show(io::IO, m::MonomialRegressionModel)
    write(io, "MonomialRegressionModel(\n  f=$(repr(m.f)),\n  objective=$(repr(m.objective)),\n  coefficient=$(repr(m.coefficient)),\n  interval=$(repr(m.interval)))")
end

function (t::Pow)(m::MonomialRegressionModel)
    x = m.f.args[1].args[3]
    replacer(a::Any) = a
    replacer(ex::Expr) = begin
        if ex.head == :call && ex.args[1] == :(^)
            # increment the power
            try
                ex.args[3]+=t.inc
            catch
                @warn "Possible invalid xform"
                @show ex
            end
        end
        return ex
    end
    m.f.args[2] = postwalk(replacer, m.f.args[2])
end

end


text:
push!(LOAD_PATH,"../src/")

using Documenter, DataStructuresAlgorithms

makedocs(sitename="Documentation for DAT038")


text:
using Interpolations

@testset "Interpolation" begin

@testset "Testing in-place interpolation" begin
    # make Gaussian line
    x1 = range(5434.0, 5435.0, step=0.01)
    y1 = GRASS.gaussian_line.(x1, mid=5434.5, width=0.1, depth=1.0)

    # new grid to re-sample on
    x2 = range(5433.5, 5435.5, step=0.05)

    # test that the interpolation runs
    itp1 = extrapolate(interpolate!(eltype(x1), (x1,), y1, Gridded(Linear())), 1.0)
    @test typeof(itp1) <: Interpolations.FilledExtrapolation

    # test boundary-conditions
    y2 = itp1.(x2)
    @test y2[1] == 1.0
    @test y2[end] == 1.0
end

end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
using CSV
using DataFrames
using StatsPlots

btc_df = CSV.File("./btc_jpy.csv") |> DataFrame
@df btc_df plot(:time, [:best_bid :best_ask :ltp], legend = :bottomright, title="BTC_JPY,Bitflyer.com 12/31/2020")

max_stars_repo_name: RockwallNest/RubyWithJulia

max_issues_repo_name: RockwallNest/RubyWithJulia

max_forks_repo_name: RockwallNest/RubyWithJulia


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using Documenter
# using DocumenterMarkdown
using GmshTools

# include("generate.jl")

makedocs(
    doctest=false,
    modules = [GmshTools],
    sitename = "GmshTools",
    # format = Markdown(),
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = [
        "Index" => "index.md",
        "Element Types" => "element_types.md",
    ],
)

deploydocs(
  repo = "github.com/shipengcheng1230/GmshTools.jl.git",
  # deps = Deps.pip("pymdown-extensions", "pygments", "mkdocs", "python-markdown-math", "mkdocs-material"),
  target = "build",
  # make = () -> run(`mkdocs build`),
)


text:
using Pkg;Pkg.activate(".")
using PairVelocities
using StaticArrays
using BenchmarkTools
using DelimitedFiles

rbins = LinRange(0.001,150.,150)
rbins_c = 0.5*(rbins[2:end] + rbins[1:end-1])

filename = "moments.csv"

positions, velocities, boxsize, redshift = read_abacus()
println(positions[:][1:20])
N = size(positions)[end]
println("N halos")
println(N)
positions = convert(Array{Float64}, positions)
velocities = convert(Array{Float64}, velocities)

moments = PairVelocities.compute_pairwise_velocity_moments(
    positions,
    velocities,
    rbins,
    boxsize,
)
DD = 2 .*moments[1][:]
bin_volume = 4/3 * pi .* diff(rbins.^3)
RR = N * (N-1)/boxsize^3 * bin_volume
println(DD)
println(RR)
xi = DD./RR .- 1.

open(filename; write=true) do f
    write(f, "r_c\txi\tv_r\tsigma_r\tsigma_t\tskewness_r\tskewness_rt\tkurtosis_r\tkurtosis_t\tkurtosis_rt\n")
    writedlm(f,
             zip(rbins_c, xi, moments[2][:],
            moments[3][:], moments[6][:],
            moments[4][:], moments[7][:],
            moments[5][:], moments[8][:], moments[9][:]))
    end

println("Wrote file !")


text:
using SparseMatrixVBCs
using Documenter

makedocs(;
    modules=[SparseMatrixVBCs],
    authors="Peter Ahrens",
    repo="https://github.com/peterahrens/SparseMatrixVBCs.jl/blob/{commit}{path}#L{line}",
    sitename="SparseMatrixVBCs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://peterahrens.github.io/SparseMatrixVBCs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/peterahrens/SparseMatrixVBCs.jl",
)


text:
using Test
using ClimateMachine.Microphysics_0M
using ClimateMachine.Microphysics
using ClimateMachine.Thermodynamics

using CLIMAParameters
using CLIMAParameters.Planet: ρ_cloud_liq, R_v, grav, R_d, molmass_ratio
using CLIMAParameters.Atmos.Microphysics
using CLIMAParameters.Atmos.Microphysics_0M

struct LiquidParameterSet <: AbstractLiquidParameterSet end
struct IceParameterSet <: AbstractIceParameterSet end
struct RainParameterSet <: AbstractRainParameterSet end
struct SnowParameterSet <: AbstractSnowParameterSet end

struct MicropysicsParameterSet{L, I, R, S} <: AbstractMicrophysicsParameterSet
    liquid::L
    ice::I
    rain::R
    snow::S
end

struct EarthParamSet{M} <: AbstractEarthParameterSet
    microphys_param_set::M
end

microphys_param_set = MicropysicsParameterSet(
    LiquidParameterSet(),
    IceParameterSet(),
    RainParameterSet(),
    SnowParameterSet(),
)

prs = EarthParamSet(microphys_param_set)
liq_prs = prs.microphys_param_set.liquid
ice_prs = prs.microphys_param_set.ice
rai_prs = prs.microphys_param_set.rain
sno_prs = prs.microphys_param_set.snow

@testset "τ_relax" begin

    @test τ_relax(liq_prs) ≈ 10
    @test τ_relax(ice_prs) ≈ 10

end

@testset "0M_microphysics" begin

    _τ_precip = τ_precip(prs)
    _qc_0 = qc_0(prs)
    _S_0 = S_0(prs)

    q_vap_sat = 10e-3
    qc = 3e-3
    q_tot = 13e-3
    frac = [0.0, 0.5, 1.0]

    # no rain if no cloud
    q = PhasePartition(q_tot, 0.0, 0.0)
    @test remove_precipitation(prs, q) ≈ 0
    @test remove_precipitation(prs, q, q_vap_sat) ≈ 0

    # rain based on qc threshold
    for lf in frac
        q_liq = qc * lf
        q_ice = (1 - lf) * qc

        q = PhasePartition(q_tot, q_liq, q_ice)

        @test remove_precipitation(prs, q) ≈
              -max(0, q_liq + q_ice - _qc_0) / _τ_precip
    end

    # rain based on supersaturation threshold
    for lf in frac
        q_liq = qc * lf
        q_ice = (1 - lf) * qc

        q = PhasePartition(q_tot, q_liq, q_ice)

        @test remove_precipitation(prs, q, q_vap_sat) ≈
              -max(0, q_liq + q_ice - _S_0 * q_vap_sat) / _τ_precip
    end
end

@testset "RainFallSpeed" begin
    # eq. 5d in Smolarkiewicz and Grabowski 1996
    # https://doi.org/10.1175/1520-0493(1996)124<0487:TTLSLM>2.0.CO;2
    function terminal_velocity_empir(
        q_rai::FT,
        q_tot::FT,
        ρ::FT,
        ρ_air_ground::FT,
    ) where {FT <: Real}
        rr = q_rai / (1 - q_tot)
        vel = FT(14.34) * ρ_air_ground^FT(0.5) * ρ^-FT(0.3654) * rr^FT(0.1346)
        return vel
    end

    # some example values
    q_rain_range = range(1e-8, stop = 5e-3, length = 10)
    ρ_air, q_tot, ρ_air_ground = 1.2, 20 * 1e-3, 1.22

    for q_rai in q_rain_range
        @test terminal_velocity(prs, rai_prs, ρ_air, q_rai) ≈
              terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground) atol =
            0.2 * terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground)

    end
end

@testset "CloudLiquidCondEvap" begin

    q_liq_sat = 5e-3
    frac = [0.0, 0.5, 1.0, 1.5]

    _τ_cond_evap = τ_relax(liq_prs)

    for fr in frac
        q_liq = q_liq_sat * fr

        @test conv_q_vap_to_q_liq_ice(
            liq_prs,
            PhasePartition(0.0, q_liq_sat, 0.0),
            PhasePartition(0.0, q_liq, 0.0),
        ) ≈ (1 - fr) * q_liq_sat / _τ_cond_evap
    end
end

@testset "CloudIceCondEvap" begin

    q_ice_sat = 2e-3
    frac = [0.0, 0.5, 1.0, 1.5]

    _τ_cond_evap = τ_relax(ice_prs)

    for fr in frac
        q_ice = q_ice_sat * fr

        @test conv_q_vap_to_q_liq_ice(
            ice_prs,
            PhasePartition(0.0, 0.0, q_ice_sat),
            PhasePartition(0.0, 0.0, q_ice),
        ) ≈ (1 - fr) * q_ice_sat / _τ_cond_evap
    end
end

@testset "RainAutoconversion" begin

    _q_liq_threshold = q_liq_threshold(rai_prs)
    _τ_acnv = τ_acnv(rai_prs)

    q_liq_small = 0.5 * _q_liq_threshold
    @test conv_q_liq_to_q_rai(rai_prs, q_liq_small) == 0.0

    q_liq_big = 1.5 * _q_liq_threshold
    @test conv_q_liq_to_q_rai(rai_prs, q_liq_big) ==
          0.5 * _q_liq_threshold / _τ_acnv
end

@testset "SnowAutoconversion" begin

    ρ = 1.0

    # above freezing temperatures -> no snow
    q = PhasePartition(15e-3, 2e-3, 1e-3)
    T = 273.15 + 30
    @test conv_q_ice_to_q_sno(prs, ice_prs, q, ρ, T) == 0.0

    # no ice -> no snow
    q = PhasePartition(15e-3, 2e-3, 0.0)
    T = 273.15 - 30
    @test conv_q_ice_to_q_sno(prs, ice_prs, q, ρ, T) == 0.0

    # no supersaturation -> no snow
    T = 273.15 - 5
    q_sat_ice = q_vap_saturation_generic(prs, T, ρ, Ice())
    q = PhasePartition(q_sat_ice, 2e-3, 3e-3)
    @test conv_q_ice_to_q_sno(prs, ice_prs, q, ρ, T) == 0.0

    # TODO - coudnt find a plot of what it should be from the original paper
    # just chacking if the number stays the same
    T = 273.15 - 10
    q_vap = 1.02 * q_vap_saturation_generic(prs, T, ρ, Ice())
    q_liq = 0.0
    q_ice = 0.03 * q_vap
    q = PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice)
    @test conv_q_ice_to_q_sno(prs, ice_prs, q, ρ, T) ≈ 1.8512022335645584e-9

end

@testset "RainLiquidAccretion" begin

    # eq. 5b in Smolarkiewicz and Grabowski 1996
    # https://doi.org/10.1175/1520-0493(1996)124<0487:TTLSLM>2.0.CO;2
    function accretion_empir(q_rai::FT, q_liq::FT, q_tot::FT) where {FT <: Real}
        rr = q_rai / (FT(1) - q_tot)
        rl = q_liq / (FT(1) - q_tot)
        return FT(2.2) * rl * rr^FT(7 / 8)
    end

    # some example values
    q_rain_range = range(1e-8, stop = 5e-3, length = 10)
    ρ_air, q_liq, q_tot = 1.2, 5e-4, 20e-3

    for q_rai in q_rain_range
        @test accretion(prs, liq_prs, rai_prs, q_liq, q_rai, ρ_air) ≈
              accretion_empir(q_rai, q_liq, q_tot) atol =
            (0.1 * accretion_empir(q_rai, q_liq, q_tot))
    end
end

@testset "Accretion" begin
    # TODO - coudnt find a plot of what it should be from the original paper
    # just chacking if the number stays the same

    # some example values
    ρ = 1.2
    q_tot = 20e-3
    q_ice = 5e-4
    q_sno = 5e-4
    q_liq = 5e-4
    q_rai = 5e-4

    @test accretion(prs, liq_prs, rai_prs, q_liq, q_rai, ρ) ≈
          1.4150106417043544e-6
    @test accretion(prs, ice_prs, sno_prs, q_ice, q_sno, ρ) ≈
          2.453070979562392e-7
    @test accretion(prs, liq_prs, sno_prs, q_liq, q_sno, ρ) ≈
          2.453070979562392e-7
    @test accretion(prs, ice_prs, rai_prs, q_ice, q_rai, ρ) ≈
          1.768763302130443e-6

    @test accretion_rain_sink(prs, ice_prs, rai_prs, q_ice, q_rai, ρ) ≈
          3.085229094251214e-5

    @test accretion_snow_rain(prs, sno_prs, rai_prs, q_sno, q_rai, ρ) ≈
          2.1705865794293408e-4
    @test accretion_snow_rain(prs, rai_prs, sno_prs, q_rai, q_sno, ρ) ≈
          6.0118801860768854e-5
end

@testset "RainEvaporation" begin

    # eq. 5c in Smolarkiewicz and Grabowski 1996
    # https://doi.org/10.1175/1520-0493(1996)124<0487:TTLSLM>2.0.CO;2
    function rain_evap_empir(
        prs::AbstractParameterSet,
        q_rai::FT,
        q::PhasePartition,
        T::FT,
        p::FT,
        ρ::FT,
    ) where {FT <: Real}

        q_sat = q_vap_saturation_generic(prs, T, ρ, Liquid())
        q_vap = q.tot - q.liq
        rr = q_rai / (1 - q.tot)
        rv_sat = q_sat / (1 - q.tot)
        S = q_vap / q_sat - 1

        ag, bg = FT(5.4 * 1e2), FT(2.55 * 1e5)
        G = FT(1) / (ag + bg / p / rv_sat) / ρ

        av, bv = FT(1.6), FT(124.9)
        F =
            av * (ρ / FT(1e3))^FT(0.525) * rr^FT(0.525) +
            bv * (ρ / FT(1e3))^FT(0.7296) * rr^FT(0.7296)

        return 1 / (1 - q.tot) * S * F * G
    end

    # example values
    T, p = 273.15 + 15, 90000.0
    ϵ = 1.0 / molmass_ratio(prs)
    p_sat = saturation_vapor_pressure(prs, T, Liquid())
    q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0))
    q_rain_range = range(1e-8, stop = 5e-3, length = 10)
    q_tot = 15e-3
    q_vap = 0.15 * q_sat
    q_ice = 0.0
    q_liq = q_tot - q_vap - q_ice
    q = PhasePartition(q_tot, q_liq, q_ice)
    R = gas_constant_air(prs, q)
    ρ = p / R / T

    for q_rai in q_rain_range
        @test evaporation_sublimation(prs, rai_prs, q, q_rai, ρ, T) ≈
              rain_evap_empir(prs, q_rai, q, T, p, ρ) atol =
            -0.5 * rain_evap_empir(prs, q_rai, q, T, p, ρ)
    end

    # no condensational growth for rain
    T, p = 273.15 + 15, 90000.0
    ϵ = 1.0 / molmass_ratio(prs)
    p_sat = saturation_vapor_pressure(prs, T, Liquid())
    q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0))
    q_rai = 1e-4
    q_tot = 15e-3
    q_vap = 1.15 * q_sat
    q_ice = 0.0
    q_liq = q_tot - q_vap - q_ice
    q = PhasePartition(q_tot, q_liq, q_ice)
    R = gas_constant_air(prs, q)
    ρ = p / R / T

    @test evaporation_sublimation(prs, rai_prs, q, q_rai, ρ, T) ≈ 0.0

end

@testset "SnowSublimation" begin
    # TODO - coudnt find a plot of what it should be from the original paper
    # just chacking if the number stays the same

    cnt = 0
    ref_val = [
        -1.9756907119482267e-7,
        1.9751292385808357e-7,
        -1.6641552112891826e-7,
        1.663814937710236e-7,
    ]
    # some example values
    for T in [273.15 + 2, 273.15 - 2]
        p = 90000.0
        ϵ = 1.0 / molmass_ratio(prs)
        p_sat = saturation_vapor_pressure(prs, T, Ice())
        q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0))

        for eps in [0.95, 1.05]
            cnt += 1

            q_tot = eps * q_sat
            q_ice = 0.0
            q_liq = 0.0
            q = PhasePartition(q_tot, q_liq, q_ice)

            q_sno = 1e-4

            R = gas_constant_air(prs, q)
            ρ = p / R / T

            @test evaporation_sublimation(prs, sno_prs, q, q_sno, ρ, T) ≈
                  ref_val[cnt]

        end
    end
end

@testset "SnowMelt" begin

    # TODO - find a good reference to compare with
    T = 273.15 + 2
    ρ = 1.2
    q_sno = 1e-4
    @test snow_melt(prs, sno_prs, q_sno, ρ, T) ≈ 9.518235437405256e-6

    # no snow -> no snow melt
    T = 273.15 + 2
    ρ = 1.2
    q_sno = 0.0
    @test snow_melt(prs, sno_prs, q_sno, ρ, T) ≈ 0

    # T < T_freeze -> no snow melt
    T = 273.15 - 2
    ρ = 1.2
    q_sno = 1e-4
    @test snow_melt(prs, sno_prs, q_sno, ρ, T) ≈ 0

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# Autogenerated wrapper script for librsync_jll for i686-linux-gnu
export rdiff, librsync

using Zlib_jll
using Bzip2_jll
using Popt_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `rdiff`
const rdiff_splitpath = ["bin", "rdiff"]

# This will be filled out by __init__() for all products, as it must be done at runtime
rdiff_path = ""

# rdiff-specific global declaration
function rdiff(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(rdiff_path)
    end
end

# Relative path to `librsync`
const librsync_splitpath = ["lib", "librsync.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
librsync_path = ""

# librsync-specific global declaration
# This will be filled out by __init__()
librsync_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const librsync = "librsync.so.2"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"librsync")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Zlib_jll.PATH_list, Bzip2_jll.PATH_list, Popt_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Zlib_jll.LIBPATH_list, Bzip2_jll.LIBPATH_list, Popt_jll.LIBPATH_list,))

    global rdiff_path = normpath(joinpath(artifact_dir, rdiff_splitpath...))

    push!(PATH_list, dirname(rdiff_path))
    global librsync_path = normpath(joinpath(artifact_dir, librsync_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global librsync_handle = dlopen(librsync_path)
    push!(LIBPATH_list, dirname(librsync_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
# Objective

"""
    AugLag2Obj

The objective function of the augmented Lagrangian model, `model`. The following are its fields:
 - `model`: the original model optimized
 - `x`: the current primal solution
 - `λ`: the current dual solution
 - `quadweight`: the current quadratic penalty
 - `f`: the current original objective value
 - `g`: the current constraint function value
"""
@params struct AugLag2Obj <: AbstractFunction
    model::Model
    x::AbstractVector # primal solution
    λ::AbstractVector # dual solution
    quadweight::Base.RefValue # quadratic penalty factor
    f::Base.RefValue
    g::AbstractVector
end

"""
    AugLag2Obj(model::Model)

Constructs the objective function from the model, `model`.
"""
function AugLag2Obj(
    model::Model;
    linweights = ones(getnineqconstraints(model)),
    quadweight = 1e-6,
    kwarg...,
)
    x = copy(getmin(model))
    f = Ref(Inf)
    g = similar(linweights)
    return AugLag2Obj(model, x, linweights, Ref(quadweight), f, g)
end

getdim(::AugLag2Obj) = 1

"""
    (obj::AugLag2Obj)(x::AbstractVector, λ::AbstractVector)

Evaluates the augmented Lagrangian at the primal solution `x` and the dual solution `λ`.
"""
function (obj::AugLag2Obj)(x::AbstractVector, λ::AbstractVector; penalise = true)
    if debugging[]
        #@show x
    end
    parent = getparent(obj)
    origobjval = getobjective(parent)(x)
    origconstrval = getineqconstraints(parent)(x)
    savexλ!(obj, x, λ)
    savefg!(obj, origobjval, origconstrval)

    linepenalty = dot(λ, origconstrval)
    quadpenalty = penalise ? NonNegSumOfSquares(obj.quadweight[])(origconstrval) : zero(linepenalty)

    return origobjval + linepenalty + quadpenalty
end

"""
    (obj::AugLag2Obj)(primaloptimizer::Function, λ::AbstractVector)

Evaluates the augmented Lagrangian at the dual solution `λ` where `primaloptimizer(λ)` returns a tuple of:
 - The optimal primal solution `x`,
 - The original objective value at the optimal `x` and the current `λ`, and
 - The original constraint values at the optimal `x` and the current `λ`,
"""
function (obj::AugLag2Obj)(primaloptimizer::Function, λ::AbstractVector; penalise = true)
    parent = getparent(obj)
    x, origobjval, origconstrval = nogradcall(primaloptimizer, λ)
    savexλ!(obj, x, λ)
    savefg!(obj, origobjval, origconstrval)

    linepenalty = dot(λ, origconstrval)
    quadpenalty = penalise ? NonNegSumOfSquares(obj.quadweight[])(origconstrval) : zero(linepenalty)

    return origobjval + linepenalty + quadpenalty
end

nogradcall(f::Function, arg) = f(arg)
ChainRulesCore.@non_differentiable nogradcall(f::Function, arg)

function savefg!(obj::AugLag2Obj, f::Real, g::AbstractVector)
    setorigobjval!(obj, f)
    setorigconstrval!(obj, g)
    return obj
end
ChainRulesCore.@non_differentiable savefg!(obj::AugLag2Obj, f::Real, g::AbstractVector)

function savexλ!(obj::AugLag2Obj, x::AbstractVector, λ::AbstractVector)
    setx!(obj, x)
    setλ!(obj, λ)
    return obj
end
ChainRulesCore.@non_differentiable savexλ!(obj::AugLag2Obj, x::AbstractVector, λ::AbstractVector)

getparent(f::AugLag2Obj) = f.model

# To be minimized
#getprimalobjective(f::AugLag2Obj) = FunctionWrapper(x -> (out = f(x, f.λ); @show out; out), 1)
getprimalobjective(f::AugLag2Obj; kwargs...) = FunctionWrapper(x -> f(x, f.λ; kwargs...), 1)

# To be maximized
getdualobjective(f::AugLag2Obj, primaloptimizer::Function) = FunctionWrapper(λ -> begin
    if debugging[]
        #@show λ
    end
    f(primaloptimizer, λ)
end, 1)

getorigobjval(f::AugLag2Obj) = f.f[]

getorigconstrval(f::AugLag2Obj) = f.g

getlinweights(f::AugLag2Obj) = f.λ

getquadweight(f::AugLag2Obj) = f.quadweight[]

getx(f::AugLag2Obj) = f.x

getλ(f::AugLag2Obj) = f.λ

function setlinweights!(f::AugLag2Obj, λ::AbstractVector)
    getlinweights(f) .= λ
    return f
end

function setquadweight!(f::AugLag2Obj, quadweight::Number)
    f.quadweight[] = quadweight
    return f
end

function setx!(f::AugLag2Obj, x::AbstractVector)
    f.x .= x
    return f
end

function setλ!(f::AugLag2Obj, λ::AbstractVector)
    f.λ .= λ
    return f
end

function setorigobjval!(f::AugLag2Obj, val::Real)
    f.f[] = val
    return f
end

function setorigconstrval!(f::AugLag2Obj, g::AbstractVector)
    f.g .= g
    return f
end

# Model

@params struct AugLag2Model <: AbstractModel
    parent::Model
    objective::AugLag2Obj
end
function AugLag2Model(
    model::Model;
    kwargs...,
)
    return AugLag2Model(model, AugLag2Obj(model; kwargs...,))
end

getparent(model::AugLag2Model) = model.parent

getmin(model::AugLag2Model) = getmin(getparent(model))

getmax(model::AugLag2Model) = getmax(getparent(model))

getobjective(model::AugLag2Model) = model.objective

getineqconstraints(::AugLag2Model) = throw("`getineqconstraints` is not defined for `AugLag2Model`.")

geteqconstraints(::AugLag2Model) = throw("`geteqconstraints` is not defined for `AugLag2Model`.")

getobjectiveconstraints(::AugLag2Model) = throw("`getobjectiveconstraints` is not defined for `AugLag2Model`.")

getlinweights(model::AugLag2Model) = getlinweights(getobjective(model))

getx(model::AugLag2Model) = getx(getobjective(model))

getλ(model::AugLag2Model) = getλ(getobjective(model))

getorigobjval(model::AugLag2Model) = getorigobjval(getobjective(model))

getorigconstrval(model::AugLag2Model) = getorigconstrval(getobjective(model))

function setlinweights!(model::AugLag2Model, λ::AbstractVector)
    setlinweights!(getobjective(model), λ)
    return model
end

getquadweight(model::AugLag2Model) = getquadweight(getobjective(model))

function setquadweight!(model::AugLag2Model, quadweight::Number)
    setquadweight!(getobjective(model), quadweight)
    return model
end

function setx!(model::AugLag2Model, x::AbstractVector)
    setx!(getobjective(model), x)
    return model
end

function setλ!(model::AugLag2Model, λ::AbstractVector)
    setλ!(getobjective(model), λ)
    return model
end

function setorigobjval!(model::AugLag2Model, val::Real)
    setorigobjval!(getobjective(model), val)
    return model
end

function setorigconstrval!(model::AugLag2Model, g::AbstractVector)
    setorigconstrval!(getobjective(model), g)
    return model
end

function getresiduals(solution::Solution, model::AugLag2Model, ::GenericCriteria)
    @unpack prevx, x, prevf, f, g = solution
    Δx = maximum(abs(x[j] - prevx[j]) for j in 1:length(x))
    Δf = abs(f - prevf)
    infeas = max(0, maximum(g))
    return Δx, Δf, infeas
end
function getresiduals(solution::Solution, model::AugLag2Model, ::KKTCriteria)
    @unpack g, λ, x = solution
    xmin, xmax = getmin(model), getmax(model)
    T = eltype(x)
    auglag = getobjective(model)
    _, ∇L = value_gradient(getprimalobjective(auglag, penalise = true), solution.x)
    res = maximum(1:length(x)) do j
        temp = ∇L[j]
        if xmin[j] >= x[j]
            return abs(min(0, temp))
        elseif x[j] >= xmax[j]
            return max(0, temp)
        else
            return abs(temp)
        end
    end
    if debugging[]
        #@show λ, g
    end
    res = max(res, maximum(abs.(λ .* g)))
    if debugging[]
        @show maximum(abs, g)
        @show maximum(abs, λ)
        @show maximum(x)
    end
    infeas = max(maximum(g), 0)
    return res, infeas
end
function getresiduals(solution::Solution, model::AugLag2Model, ::IpoptCriteria)
    @unpack g, λ, x = solution
    if debugging[]
        #@show x
    end
    xmin, xmax = getmin(model), getmax(model)
    T = eltype(x)
    n, m, s = length(x), length(λ), zero(T)
    auglag = getobjective(model)
    _, ∇L = value_gradient(getprimalobjective(auglag, penalise = true), solution.x)
    res = maximum(1:n) do j
        temp = ∇L[j]
        if xmin[j] >= x[j]
            dj = temp
            s += max(dj, 0)
            return abs(min(0, dj))
        elseif x[j] >= xmax[j]
            yj = -temp
            s += max(yj, 0)
            return abs(min(0, yj))
        else
            return abs(temp)
        end
    end
    sd = max(100, (sum(abs, λ) + s) / (n + m)) / 100
    res = max(res, maximum(abs.(λ .* g)))
    res = res / sd
    infeas = max(maximum(g), 0)
    if debugging[]
        println("Agg infeas = ", infeas)
    end
    return res, infeas
end


text:
# Create Napari.view_* and Napari.add_* that alias into the Python versions
# This allows Julia code to take advantage of multiple dispatch by overloading
# these functions for specific types.

# The code below generates the add_* and view_* functions in the comments below
for layer in layers
    local view_name = Symbol("view_", layer)
    local add_name = Symbol("add_", layer)
    @eval begin
        $view_name( data, args...; kwargs... ) = napari.$view_name( data, args...; kwargs... )
        $add_name( viewer, data, args...; kwargs... ) = viewer.$add_name( data, args...; kwargs... )
    end
end

# view_image(   data, args...; kwargs... ) = napari.view_image(   data, args...; kwargs... )
# view_points(  data, args...; kwargs... ) = napari.view_points(  data, args...; kwargs... )
# view_labels(  data, args...; kwargs... ) = napari.view_labels(  data, args...; kwargs... )
# view_shapes(  data, args...; kwargs... ) = napari.view_shapes(  data, args...; kwargs... )
# view_surface( data, args...; kwargs... ) = napari.view_surface( data, args...; kwargs... )
# view_vectors( data, args...; kwargs... ) = napari.view_vectors( data, args...; kwargs... )
# view_tracks(  data, args...; kwargs... ) = napari.view_tracks(  data, args...; kwargs... )

# add_image(   viewer, data, args...; kwargs... ) = napari.add_image(   viewer, data, args...; kwargs... )
# add_points(  viewer, data, args...; kwargs... ) = napari.add_points(  viewer, data, args...; kwargs... )
# add_labels(  viewer, data, args...; kwargs... ) = napari.add_labels(  viewer, data, args...; kwargs... )
# add_shapes(  viewer, data, args...; kwargs... ) = napari.add_shapes(  viewer, data, args...; kwargs... )
# add_surface( viewer, data, args...; kwargs... ) = napari.add_surface( viewer, data, args...; kwargs... )
# add_vectors( viewer, data, args...; kwargs... ) = napari.add_vectors( viewer, data, args...; kwargs... )
# add_tracks(  viewer, data, args...; kwargs... ) = napari.add_tracks(  viewer, data, args...; kwargs... )

# Place holder modules mirroring upstream
# I'm not really sure if this is actually going to be useful
module Layers

    # module Image
    include("image.jl")

    module Labels
    end
    module Points
    end
    module Shapes
    end
    module Surface
    end
    module Tracks
    end
    module Utils
    end
    module Vectors
    end
end