text:
module CompositeStructures

using ..Prerequisites: rawtype
using ..TypeTraits: efficientoperations

export CompositeTuple, CompositeNTuple, CompositeVector, CompositeDict, NamedContainer

"""
    CompositeTuple{T<:Tuple}

A composite tuple is a tuple that is implemented by including an ordinary `Tuple` as one of its attributes with the name `:contents`.
"""
abstract type CompositeTuple{T<:Tuple} end
Base.length(::CompositeTuple{T}) where T = fieldcount(T)
Base.length(::Type{<:CompositeTuple{T}}) where T = fieldcount(T)
Base.eltype(::CompositeTuple{T}) where T = eltype(T)
Base.eltype(::Type{<:CompositeTuple{T}}) where T = eltype(T)
@generated function Base.hash(ct::CompositeTuple, h::UInt)
    ts = Expr(:tuple, [:(getfield(ct, $i)) for i = 1:fieldcount(ct)]...)
    return :(hash($ts, h))
end
Base.:(==)(ct1::CompositeTuple, ct2::CompositeTuple) = ==(efficientoperations, ct1, ct2)
Base.isequal(ct1::CompositeTuple, ct2::CompositeTuple) = isequal(efficientoperations, ct1, ct2)
Base.getindex(ct::CompositeTuple, i::Union{<:Integer, CartesianIndex}) = getfield(ct, :contents)[i]
@generated function Base.getindex(ct::CompositeTuple, inds)
    exprs = [(name == :contents) ? :(getfield(ct, :contents)[inds]) : :(getfield(ct, $i)) for (i, name) in enumerate(fieldnames(ct))]
    return :(rawtype(typeof(ct))($(exprs...)))
end
Base.lastindex(ct::CompositeTuple) = lastindex(getfield(ct, :contents))
Base.iterate(ct::CompositeTuple) = iterate(getfield(ct, :contents))
Base.iterate(ct::CompositeTuple, state) = iterate(getfield(ct, :contents), state)
Base.iterate(rv::Iterators.Reverse{<:CompositeTuple}, state=length(rv.itr)) = (state < 1) ? nothing : (rv.itr[state], state-1)
Base.keys(ct::CompositeTuple) = keys(getfield(ct, :contents))
Base.values(ct::CompositeTuple) = values(getfield(ct, :contents))
Base.pairs(ct::CompositeTuple) = pairs(getfield(ct, :contents))
@generated function Base.reverse(ct::CompositeTuple)
    exprs = [(name == :contents) ? :(reverse(getfield(ct, :contents))) : :(getfield(ct, $i)) for (i, name) in enumerate(fieldnames(ct))]
    return :(rawtype(typeof(ct))($(exprs...)))
end
Base.convert(::Type{Tuple}, ct::CompositeTuple) = getfield(ct, :contents)

"""
    CompositeNTuple{N, T}

A composite ntuple is a ntuple that is implemented by including an ordinary `NTuple` as one of its attributes with the name `:contents`.

Alias for `CompositeTuple{NTuple{N, T}}`.
"""
const CompositeNTuple{N, T} = CompositeTuple{NTuple{N, T}}

"""
    CompositeVector{T}

A composite vector is a vector that is implemented by including an ordinary `Vector` as one of its attributes with the name `:contents`.
"""
abstract type CompositeVector{T} <: AbstractVector{T} end
Base.size(cv::CompositeVector) = size(getfield(cv, :contents))
Base.size(cv::CompositeVector, i) = size(getfield(cv, :contents), i)
Base.length(cv::CompositeVector) = length(getfield(cv, :contents))
Base.:(==)(cv1::CompositeVector, cv2::CompositeVector) = ==(efficientoperations, cv1, cv2)
Base.isequal(cv1::CompositeVector, cv2::CompositeVector) = isequal(efficientoperations, cv1, cv2)
Base.getindex(cv::CompositeVector, i::Union{<:Integer, CartesianIndex}) = getfield(cv, :contents)[i]
@generated function Base.getindex(cv::CompositeVector, inds)
    exprs = [(name == :contents) ? :(getfield(cv, :contents)[inds]) : :(getfield(cv, $i)) for (i, name) in enumerate(fieldnames(cv))]
    return :(rawtype(typeof(cv))($(exprs...)))
end
Base.lastindex(cv::CompositeVector) = lastindex(getfield(cv, :contents))
Base.setindex!(cv::CompositeVector, value, inds) = (getfield(cv, :contents)[inds] = value)
Base.push!(cv::CompositeVector, values...) = (push!(getfield(cv, :contents), values...); cv)
Base.pushfirst!(cv::CompositeVector, values...) = (pushfirst!(getfield(cv, :contents), values...); cv)
Base.insert!(cv::CompositeVector, index::Integer, value) = (insert!(getfield(cv, :contents), index, value); cv)
Base.append!(cv::CompositeVector, values) = (append!(getfield(cv, :contents), values); cv)
Base.prepend!(cv::CompositeVector, values) = (prepend!(getfield(cv, :contents), values); cv)
Base.splice!(cv::CompositeVector, index::Integer, replacement=Base._default_splice) = splice!(getfield(cv, :contents), index, replacement)
@generated function Base.splice!(cv::CompositeVector, range::UnitRange{<:Integer}, replacement=Base._default_splice)
    exprs = [(name == :contents) ? :(splice!(getfield(cv, :contents), range, replacement)) : :(getfield(cv, $i)) for (i, name) in enumerate(fieldnames(cv))]
    return :(rawtype(typeof(cv))($(exprs...)))
end
Base.deleteat!(cv::CompositeVector, indices) = (deleteat!(getfield(cv, :contents), indices); cv)
Base.pop!(cv::CompositeVector) = pop!(getfield(cv, :contents))
Base.popfirst!(cv::CompositeVector) = popfirst!(getfield(cv, :contents))
Base.empty!(cv::CompositeVector) = (empty!(getfield(cv, :contents)); cv)
@generated function Base.empty(cv::CompositeVector)
    exprs = [(name == :contents) ? :(empty(getfield(cv, :contents))) : :(getfield(cv, $i)) for (i, name) in enumerate(fieldnames(cv))]
    return :(rawtype(typeof(cv))($(exprs...)))
end
@generated function Base.reverse(cv::CompositeVector)
    exprs = [(name == :contents) ? :(reverse(getfield(cv, :contents))) : :(getfield(cv, $i)) for (i, name) in enumerate(fieldnames(cv))]
    return :(rawtype(typeof(cv))($(exprs...)))
end
@generated function Base.similar(cv::CompositeVector, dtype::Type=eltype(cv), dims::Tuple{Vararg{Int}}=size(cv))
    exprs = [(name == :contents) ? :(similar(getfield(cv, :contents), dtype, dims)) : :(getfield(cv, $i)) for (i, name) in enumerate(fieldnames(cv))]
    return :(rawtype(typeof(cv))($(exprs...)))
end
Base.iterate(cv::CompositeVector, state=1) = iterate(getfield(cv, :contents), state)
Base.keys(cv::CompositeVector) = keys(getfield(cv, :contents))
Base.values(cv::CompositeVector) = values(getfield(cv, :contents))
Base.pairs(cv::CompositeVector) = pairs(getfield(cv, :contents))
Base.convert(::Type{Vector}, cv::CompositeVector) = getfield(cv, :contents)

"""
    CompositeDict{K, V}

A composite dict is a dict that is implemented by including an ordinary `Dict` as one of its attributes with the name `:contents`.
"""
abstract type CompositeDict{K, V} <: AbstractDict{K, V} end
Base.isempty(cd::CompositeDict) = isempty(getfield(cd, :contents))
Base.length(cd::CompositeDict) = length(getfield(cd, :contents))
Base.haskey(cd::CompositeDict, key) = haskey(getfield(cd, :contents), key)
Base.in(p::Pair, cd::CompositeDict, valcmp=(==)) = in(p, getfield(cd, :contents), valcmp)
Base.:(==)(cd1::CompositeDict, cd2::CompositeDict) = ==(efficientoperations, cd1, cd2)
Base.isequal(cd1::CompositeDict, cd2::CompositeDict) = isequal(efficientoperations, cd1, cd2)
Base.get(cd::CompositeDict, key, default) = get(getfield(cd, :contents), key, default)
Base.get(f::Base.Callable, cd::CompositeDict, key) = get(f, getfield(cd, :contents), key)
Base.getkey(cd::CompositeDict, key, default) = getkey(getfield(cd, :contents), key, default)
Base.getindex(cd::CompositeDict{K, V}, index::K) where {K, V} = getfield(cd, :contents)[index]
Base.push!(cd::CompositeDict, ps::Pair...) = (push!(getfield(cd, :contents), ps...); cd)
Base.get!(cd::CompositeDict, key, default) = get!(getfield(cd, :contents), key, default)
Base.get!(default::Union{Function, Type}, cd::CompositeDict, key) = get!(default, getfield(cd, :contents), key)
Base.setindex!(cd::CompositeDict{K, V}, value::V, index::K) where {K, V} = (getfield(cd, :contents)[index] = value)
Base.pop!(cd::CompositeDict) = pop!(getfield(cd, :contents))
Base.pop!(cd::CompositeDict, key) = pop!(getfield(cd, :contents), key)
Base.pop!(cd::CompositeDict, key, default) = pop!(getfield(cd, :contents), key, default)
Base.delete!(cd::CompositeDict, key) = (delete!(getfield(cd, :contents), key); cd)
Base.empty!(cd::CompositeDict) = (empty!(getfield(cd, :contents)); cd)
Base.merge(cd::CD, others::CD...) where CD <: CompositeDict = merge!(empty(cd), cd, others...)
Base.merge(combine::Function, cd::CD, others::CD...) where {CD<:CompositeDict} = merge!(combine, empty(cd), cd, others...)
@generated function Base.empty(cd::CompositeDict)
    exprs = [(name == :contents) ? :(Dict{keytype(cd), valtype(cd)}()) : :(getfield(cd, $i)) for (i, name) in enumerate(fieldnames(cd))]
    return :(rawtype(typeof(cd))($(exprs...)))
end
Base.iterate(cd::CompositeDict) = iterate(getfield(cd, :contents))
Base.iterate(cd::CompositeDict, state) = iterate(getfield(cd, :contents), state)
Base.keys(cd::CompositeDict) = keys(getfield(cd, :contents))
Base.values(cd::CompositeDict) = values(getfield(cd, :contents))
Base.pairs(cd::CompositeDict) = pairs(getfield(cd, :contents))
Base.convert(::Type{Dict}, cd::CompositeDict) = getfield(cd, :contents)

"""
    NamedContainer{T, Names} = NamedTuple{Names, <:Tuple{Vararg{T}}}

NamedContainer is just a wrapper of Julia NamedTuple, but not a composite type.
"""
const NamedContainer{T, Names} = NamedTuple{Names, <:Tuple{Vararg{T}}}

"""
    NamedContainer{Names}(contents) where Names -> NamedTuple{Names, typeof(contents)}

Construct a named container.
"""
@generated function NamedContainer{Names}(contents::Tuple) where Names
    @assert length(Names) == fieldcount(contents) "NamedContainer error: dismatched length between names and contents."
    fieldcount(contents) == 0 && return NamedTuple()
    return Expr(:tuple, [:($name = contents[$i]) for (i, name) in enumerate(Names)]...)
end

end #module


text:
# VectorLogger is a Logging destintation that's an im-memory vector.

module VectorLogging

include("vector_logger.jl")

end

ext: jl


text:
# generate the distribution plots found in
# https://arxiv.org/abs/2201.05967

using DyadicKDE
include("./plot_helpers.jl")

# specify parameters
n_evals = 100
degeneracies = ["total", "partial", "none"]
evals = collect(range(-2.0, stop=2.0, length=n_evals))
ps = Dict(
    "total" => [0.5, 0.0, 0.5],
    "partial" => [0.25, 0.0, 0.75],
    "none" => [0.2, 0.2, 0.6],
)

println("Making plots")
linewidth = 1.0
y_lim = [-0.01, 0.45]
handle_f = PyPlot.matplotlib.lines.Line2D([0], [0], color="k", lw=linewidth, label="\$f_W(w)\$")
handle_Var_f_given_A = PyPlot.matplotlib.lines.Line2D([0], [0], color="k", lw=linewidth,
    label="\$\\mathrm{Var}[f_{W \\mid A}(w \\mid A_i)]^{1/2}\$")
handles = [handle_f, handle_Var_f_given_A]

for degen in degeneracies

    p = ps[degen]
    f = generate_f(evals, p)
    Var_f_given_A = generate_Var_f_given_A(evals, p)

    # plot f
    fig, ax = plt.subplots(figsize=(4,4))
    ax.plot(evals, f, color = "black", linewidth=linewidth)

    # plot Var_f_given_A
    ax.plot(evals, Var_f_given_A.^(0.5),
        color = "black", linewidth=linewidth, linestyle=(0, (1,1)))

    # save plot
    PyPlot.xlabel("Evaluation point, \$w\$")
    plt.ylim(y_lim)
    x_min = minimum(evals)
    x_max = maximum(evals)
    plt.xlim((x_min, x_max))
    plt.yticks(range(0.0, stop=0.4, step=0.1))
    legend(handles=handles, loc="upper left")
    plt.ylabel("Density", labelpad=4.0)
    plt.tight_layout()
    PyPlot.savefig("distribution_plot_$degen.pdf")
    close("all")

end


text:
using Documenter, JuDoc, JuDocTemplates

makedocs(
    modules = [JuDoc, JuDocTemplates],
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        # custom CSS
        assets = ["assets/custom.css"]
        ),
    sitename = "JuDoc.jl",
    authors  = "Thibaut Lienart",
    pages    = [
        "Home" => "index.md",
        "Manual" => [
             "Workflow" => "man/workflow.md",
             "Syntax" => "man/syntax.md",
             "Templating" => "man/templating.md",
             "Themes" => "man/themes.md",
             "Contributing" => "man/contrib.md",
             "Troubleshooting" => "man/troubleshooting.md",
            ],
        "Library" => [
            "Design"    => "lib/design.md",
            "Public"    => "lib/public.md",
            "Internals" => "lib/internals.md",
            ],
        ], # end page
)

deploydocs(
    repo = "github.com/tlienart/JuDoc.jl.git",
)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PkgTestSuite.jl/blob/master/LICENSE

using PkgTestSuite: checkheader
using Base.Test

new_file_path = Pkg.dir("PkgTestSuite", "src", "new_file_with_wrong_license_header.jl")
isfile(new_file_path) && rm(new_file_path)

@test checkheader("PkgTestSuite")

test_file_contents = """
# This file is not a part of JuliaFEM.
# License is not MIT: don't see https://github.com/JuliaFEM/PkgTestSuite.jl/blob/master/LICENSE

function foo()
    return 1
end
"""
open(new_file_path, "w") do fid
    write(fid, test_file_contents)
end
@test_throws Exception checkheader("PkgTestSuite")
isfile(new_file_path) && rm(new_file_path)


text:
using Printf
using Plots; pyplot(fmt = :png)
using SolidStateDetectors
using Unitful
using StatsPlots

include("WaveformAnalysis.jl")

figSize = (1400, 1000)
extension = ".png"

function GetDetectorFacePosition(sim::Simulation)
    #thickness = sim.detector.semiconductors[1].geometry.h
    thickness = 2e-3
    #z = sim.detector.semiconductors[1].geometry.translate.z
    z = thickness/2
    zBottom = z - thickness/2
    zTop = z + thickness/2

    return (zBottom, zTop)
end

function PlotEvents(events::Vector{Event}, dir::String="", prefix::String="", suffix::String="")::Nothing
    if !isempty(events)
        try
            for i in 1:size(events[1].waveforms, 1) # Loop over contacts
                # Plot integrated charge
                pIC = plot(events[1].waveforms[1].time / 1u"ns", [x.waveforms[i].value for x in events],
                xlabel="Time [ns]", ylabel="Integrated charge [arb.]", linecolor="red",
                linealpha=1/sqrt(length(events)), legend=false)

                #savefig(dir * prefix * "integrated_charge_contact_$i" * suffix * extension)

                # Histogram rise times

                riseTimes = CollectRiseTimes(events)
                pRT = histogram(riseTimes / 1u"ns", bins=:scott, xlabel="Rise time [ns]", ylabel="Counts [arb.]")

                #savefig(dir * prefix * "risetime_hist_$i" * suffix * extension)

                # Histogram maximal integrated current

                maxAmplitude = [maximum(abs.(x.waveforms[i].value)) for x in events]
                # println(maxAmplitude)
                # println(length(riseTimes), length(maxAmplitude))
                pMA = histogram(maxAmplitude, bins=:scott, xlabel="Maximal amplitude", ylabel="Counts [arb.]")

                #savefig(dir * prefix * "amplitude_hist_$i" * suffix * extension)

                # 2D Amplitude - Rise time interpolated scatter plot
                pKDE = marginalkde(riseTimes / 1u"ns", maxAmplitude)

                #savefig(dir * prefix * "rise_time_amplitude_kde_$i" * suffix * extension)

                plot(pRT, pIC, pKDE, pMA, layout=(2, 2), size=figSize)

                savefig(dir * prefix * "event_contact_$i" * suffix * extension)
            end
        catch

        end
    end
end

function PlotGeometry(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    @info "Plotting detector geometry"
    plot(sim.detector, size=figSize)

    savefig(dir * prefix * "geometry" * suffix * extension)
end

function PlotElectricPotential(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, thickness = GetDetectorFacePosition(sim)
    zTop = zBottom + thickness
    if !ismissing(sim.electric_potential)
        plot(plot(sim.electric_potential, y=0, aspect_ratio = :none), plot(sim.electric_potential, x=0, aspect_ratio = :none),
        plot(sim.electric_potential, z=zBottom), plot(sim.electric_potential, z=zTop),
        layout=(2, 2), size=figSize)

        savefig(dir * prefix * "electric_potential" * suffix * extension)
    end
end

function Plot2DWeightingPotential(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    for i in 1:size(sim.weighting_potentials, 1)
        if !ismissing(sim.weighting_potentials[i])
            plot(plot(sim.weighting_potentials[i], y=0, aspect_ratio = :none), plot(sim.weighting_potentials[i], x=0, aspect_ratio = :none),
            plot(sim.weighting_potentials[i], z=zBottom), plot(sim.weighting_potentials[i], z=zTop),
            layout=(2, 2), size=figSize)

            p = @sprintf "weighting_potential_%d" i

            savefig(dir * prefix * p * suffix * extension)
        end
    end
end

function PlotElectricField(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    if !ismissing(sim.electric_field)
        p_xz = plot(sim.electric_field, y=0, aspect_ratio = :none)
        plot_electric_fieldlines!(sim)
        p_yz = plot(sim.electric_field, x=0, aspect_ratio = :none)
        plot_electric_fieldlines!(sim)
        p_xy_bottom = plot(sim.electric_field, z=zBottom)
        p_xy_top = plot(sim.electric_field, z=zTop)
        plot(p_xz, p_yz, p_xy_bottom, p_xy_top, layout=(2, 2), size=figSize)

        savefig(dir * prefix * "electric_field" * suffix * extension)
    end
end

function PlotFields(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    @info "Plotting fields"
    PlotElectricPotential(sim, dir, prefix, suffix)
    PlotWeightingPotential(sim, dir, prefix, suffix)
    PlotElectricField(sim, dir, prefix, suffix)
end

function PlotChargeDensity(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    if !ismissing(sim.ρ)
        plot(plot(sim.ρ, x=0, aspect_ratio = :none), plot(sim.ρ, y=0, aspect_ratio = :none), plot(sim.ρ, z=zBottom), plot(sim.ρ, z=zTop), layout=(2,2), size=figSize)
        savefig(dir * prefix * "rho" * suffix * extension)
    end
end

function PlotDielectric(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    if !ismissing(sim.ϵ)
        plot(plot(sim.ϵ, x=0, aspect_ratio = :none), plot(sim.ϵ, y=0, aspect_ratio = :none), plot(sim.ϵ, z=zBottom), plot(sim.ϵ, z=zTop), layout=(2,2), size=figSize)
        savefig(dir * prefix * "epsilon" * suffix * extension)
    end
end

function PlotPointType(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    zBottom, zTop = GetDetectorFacePosition(sim)
    if !ismissing(sim.point_types)
        plot(plot(sim.point_types, x=0, aspect_ratio = :none), plot(sim.point_types, y=0, aspect_ratio = :none), plot(sim.point_types, z=zBottom), plot(sim.point_types, z=zTop), layout=(2,2), size=figSize)
        savefig(dir * prefix * "point_types" * suffix * extension)
    end
end

function PlotMaterialProperties(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    @info "Plotting material properties"

    PlotChargeDensity(sim, dir, prefix, suffix)
    PlotDielectric(sim, dir, prefix, suffix)
    PlotPointType(sim, dir, prefix, suffix)

    zBottom, zTop = GetDetectorFacePosition(sim)

    if !(ismissing(sim.ρ) || ismissing(sim.ϵ) || ismissing(sim.point_types))
        plot(plot(sim.ρ, x=0), plot(sim.ϵ, x=0), plot(sim.point_types, x=0), layout=(1, 3), size=figSize)
        savefig(dir * prefix * "material_properties_yz" * suffix * extension)
        plot(plot(sim.ρ, y=0), plot(sim.ϵ, y=0), plot(sim.point_types, y=0), layout=(1, 3), size=figSize)
        savefig(dir * prefix * "material_properties_xz" * suffix * extension)
        plot(plot(sim.ρ, z=zBottom), plot(sim.ϵ, z=zBottom), plot(sim.point_types, z=zBottom), layout=(1, 3), size=figSize)
        savefig(dir * prefix * "material_properties_xy_bottom" * suffix * extension)
        plot(plot(sim.ρ, z=zTop), plot(sim.ϵ, z=zTop), plot(sim.point_types, z=zTop), layout=(1, 3), size=figSize)
        savefig(dir * prefix * "material_properties_xy_top" * suffix * extension)
    end
end

function PlotCV(sims::AbstractVector{Simulation}, dir::String="", prefix::String="", suffix::String="")::Nothing
    cv = zeros(size(sims, 1), 2)
    for i in 1:size(sims, 1)
        cv[i, 1] = sims[i].detector.contacts[1].potential
        cv[i, 2] = calculate_capacitance(sims[i]) / 1u"pF"
    end

    sort!(cv, dims=1)

    pCV = plot(cv[:, 1], cv[:, 2], xlabel="Voltage [V]", ylabel="Capacitance [pF]", size=figSize)

    savefig(dir * prefix * "cv" * suffix * extension)

    dCm2dV = zeros(size(sims, 1)-1, 2)
    for i in 1:size(dCm2dV, 1)
        dCm2dV[i, 1] = cv[i, 1]
        dCm2dV[i, 2] = ((cv[i+1, 2])^(-2)-cv[i, 2]^(-2))/(cv[i+1, 1] - cv[i, 1])
    end

    pDerivCV = plot(dCm2dV[:, 1], dCm2dV[:, 2], xlabel="Voltage [V]", ylabel="d(1/C^2)/dV", size=figSize)

    savefig(dir * prefix * "deriv_cv" * suffix * extension)

    plot(pCV, pDerivCV, layout=(1, 2), size=figSize)

    savefig(dir * prefix * "cv_deriv_cv" * suffix * extension)
end

function Plot1DWeightingPotential(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    # Edge and corner plots only work for hexagons
    @info "Plotting weighting potentials vs z"

    for i in 2:length(sim.weighting_potentials)
        grid = sim.weighting_potentials[i].grid
        gridX = collect(grid.x)
        gridY = collect(grid.y)
	center = sim.detector.contacts[i].geometry.translate
        R = sim.detector.contacts[i].geometry.rOuter
        r = sqrt(3)*R/2
        edge = (0.9.*r + center[1], 0.0 + center[2])
        corner = (0.9.*r + center[1], R/2 + center[2])

        thresholdX = maximum(gridX[2:end]-gridX[1:end-1])/2
        thresholdY = maximum(gridY[2:end]-gridY[1:end-1])/2

        #println(thresholdX, " ", thresholdY)
        x0 = findlast(x->abs(x-center[1]) < thresholdX, collect(grid.x))
        y0 = findlast(x->abs(x-center[2]) < thresholdY, collect(grid.y))
        xc = findlast(x->abs(x-corner[1]) < thresholdX, collect(gridX))
        yc = findlast(x->abs(x-corner[2]) < thresholdY, collect(gridY))
        xe = findlast(x->abs(x-edge[1]) < thresholdX, collect(gridX))
        ye = findlast(x->abs(x-edge[2]) < thresholdY, collect(gridY))

	println(center)
        println(grid[x0, y0, 1])
        println(grid[xe, ye, 1])
        println(grid[xc, yc, 1])

        wp0 = sim.weighting_potentials[i][x0, y0, :]
        wpe = sim.weighting_potentials[i][xe, ye, :]
        wpc = sim.weighting_potentials[i][xc, yc, :]

        plot(collect(grid.z), [wp0, wpe, wpc], labels = ["Center" "Edge" "Corner"],
        xlabel = "Z [m]", ylabel="Weighting potential", size=figSize)

        savefig(dir * prefix * "weighting_potential_vs_z_$i" * suffix * extension)
    end
end

function Plot1DWeightingPotentialRange(sim::Simulation, initPos::Real, finalPos::Real,length::Integer,contact::Integer=2, dir::String="", prefix::String="", suffix::String="")::Nothing

    @info "Plotting weighting potentials vs z"

    grid = sim.weighting_potentials[contact].grid
    gridX = collect(grid.x)
    gridY = collect(grid.y)

    edge = 6e-3
    tolerance = 5e-4

    thresholdX = maximum(gridX[2:end]-gridX[1:end-1])/2
    println(gridX[end])
    thresholdY = maximum(gridY[2:end]-gridY[1:end-1])/2
    y0 = findlast(x->x==0, collect(grid.y))

    wps = []
    posRange = range(initPos, stop=finalPos, length=length)
    for i in 1:size(posRange,1)
    	xi = posRange[i]
	println(xi)
	x0 = findlast(x->(abs(x-xi)<thresholdX && x<gridX[end]), collect(grid.x))

    	println(grid[x0, y0, 1])

        push!(wps,sim.weighting_potentials[contact][x0, y0, :])

    end

    posStr = map(x->string(x),posRange)
    poslabels = reshape(posStr, (1,size(posStr)[1]))

    plot(collect(grid.z), wps, labels = poslabels, xlabel = "Z [m]", ylabel="Weighting potential", size=figSize)

    savefig(dir * prefix * "weighting_potential_vs_z_range$contact" * suffix * extension)

end

function PlotIDP(sim::Simulation, dir::String="", prefix::String="", suffix::String="")::Nothing
    @info "Plotting IDP along z and top of detector"

    bulk = sim.detector.semiconductors[1]
    bulk_geom = bulk.geometry
    center = bulk_geom.translate === missing ? (0,0) : bulk_geom.translate
    R = bulk_geom.r_interval.right
    H = bulk_geom.z_interval.right

    zs = 0:1e-8:H
    charge_vecz = similar(zs)
    for iz in eachindex(zs)
        pt = CartesianPoint{T}(0,0,zs[iz])
        charge_vecz[iz] = pt in bulk ?
                            SolidStateDetectors.get_charge_density(bulk, pt) : 0
    end

    xs = -R:1e-7:R
    charge_vecx = similar(xs)
    charge_vecr = similar(xs)
    for ix in eachindex(xs)
        pt = CartesianPoint{T}(xs[ix]+center[1], 0, H-1e-6)
        charge_vecx[ix] = pt in bulk ?
                            SolidStateDetectors.get_charge_density(bulk, pt) : 0
        pt = CartesianPoint{T}(xs[ix] + center[1], 0, H/2)
        charge_vecr[ix] = pt in bulk ? SolidStateDetectors.get_charge_density(bulk, pt) : 0
    end

    plot(plot(zs, charge_vecz), plot(xs, charge_vecx), plot(xs, charge_vecr), layout=(1, 3), size=figSize)
    savefig(dir * prefix * "IDP_z_and_top_and_r" * suffix * extension)
end

function Plot1DElectricRange(sim::Simulation, initPos::Real, finalPos::Real,length::Integer, dir::String="", prefix::String="", suffix::String="")::Nothing

    @info "Plotting weighting potentials vs z"

    grid = sim.electric_field.grid
    gridX = collect(grid.x)
    gridY = collect(grid.y)
    gridZ = collect(grid.z)

    edge = 6e-3
    tolerance = 5e-4

    thresholdX = maximum(gridX[2:end]-gridX[1:end-1])/2
    println(gridX[end])
    thresholdY = maximum(gridY[2:end]-gridY[1:end-1])/2
    y0 = findlast(x->x==0, collect(grid.y))

    ezs = []
    posRange = range(initPos, stop=finalPos, length=length)
    for i in 1:size(posRange,1)
    	xi = posRange[i]
	println(xi)
	x0 = findlast(x->(abs(x-xi)<thresholdX && x<gridX[end]), collect(grid.x))

    	println(grid[x0, y0, 1])
        ez = []
	for z in 1:size(gridZ,1)
            push!(ez,sim.electric_field[x0, y0, z][3])
	end

	push!(ezs, ez)

    end

    posStr = map(x->string(x),posRange)
    poslabels = reshape(posStr, (1,size(posStr)[1]))

    plot(collect(grid.z), ezs, labels = poslabels, xlabel = "Z [m]", ylabel="Ez", size=figSize)

    savefig(dir * prefix * "electric_field_vs_z_range" * suffix * extension)

end

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
using Pkg.TOML

@testset "Quadrature.jl" begin
	params = TOML.parsefile("Quadrature.toml")

    magnetization = MPISimulations.initialize(MPISimulations.SeparableMagnetization,"SeparableMagnetization1",params,Float64)
    aob = magnetization.tracerDistribution
    gl1 = MPISimulations.QuadratureNodes(magnetization)
    gl2 = MPISimulations.QuadratureNodes(aob)

    @test gl1.nodesx == gl2.nodesx
    @test gl1.weightsx == gl2.weightsx
    @test gl1.nodesy == gl2.nodesy
    @test gl1.weightsy == gl2.weightsy
    @test gl1.nodesz == gl2.nodesz
    @test gl1.weightsz == gl2.weightsz
    @test gl1.n == gl2.n

    for density in [1/2e-2,1/2e-3,1/2e-4,1/2e-5]
        nodesAndWeights = MPISimulations.QuadratureNodes(aob, density)
        result = 0.0
        for (r,w) in nodesAndWeights
            result += w*MPISimulations.tracerConcentration(aob,r,0.0)
        end
        @test result ≈ 4.8e-8
    end

    magnetization = MPISimulations.initialize(MPISimulations.SeparableMagnetization,"SeparableMagnetization2",params,Float64)
    sphere = magnetization.tracerDistribution
    gl1 = MPISimulations.QuadratureNodes(magnetization)
    gl2 = MPISimulations.QuadratureNodes(sphere)
    @test gl1.nodesx == gl2.nodesx
    @test gl1.nodesy == gl2.nodesy
    @test gl1.nodesz == gl2.nodesz
    @test gl1.weight == gl2.weight
    @test gl1.n == gl2.n

    for (density,rtol) in [(1/2e-2,1e0),(1/2e-3,1e0),(1/2e-4,1e-1),(1/2e-5,1e-3)]
        nodesAndWeights = MPISimulations.QuadratureNodes(sphere, density)
        result = 0.0
        for (r,w) in nodesAndWeights
            result += w*MPISimulations.tracerConcentration(sphere,r,0.0)
        end
        @test result ≈ 4/3*π*(sphere.radius)^3 rtol=rtol
    end
end


text:
using NNlib: conv, ∇conv_filter, ∇conv_data, ∇maxpool, maxpool

@testset "conv2d" begin
    x = reshape(Float64[1:20;], 5, 4, 1, 1)
    w = reshape(Float64[1:4;], 2, 2, 1, 1)

    @test squeeze(conv(x, w),(3,4)) == [
        29 79 129;
        39 89 139;
        49 99 149;
        59 109 159.]

    @test squeeze(conv(x, w; stride=2),(3,4)) == [
        29 129;
        49 149.]

    @test squeeze(conv(x, w; pad=1),(3,4)) == [
        1.0   9.0   29.0   49.0   48.0;
        4.0  29.0   79.0  129.0  115.0;
        7.0  39.0   89.0  139.0  122.0;
        10.0  49.0   99.0  149.0  129.0;
        13.0  59.0  109.0  159.0  136.0;
        10.0  40.0   70.0  100.0   80.0
    ]

    # for gradients, check only size
    # correctness of gradients is cross-checked with CUDNN.jl
    # (it's assumed convolution code won't change often)

    @test size(∇conv_filter(reshape(rand(4,3), 4, 3, 1, 1), x, w)) == size(w)
    @test size(∇conv_data(reshape(rand(4,3), 4, 3, 1, 1), x, w)) == size(x)

    # Test that stride/pad work backward as well
    y = conv(x, w; stride=2, pad=1)
    @test size(y) == (3, 3, 1, 1)
    @test size(∇conv_filter(y, x, w; stride=2, pad=1)) == size(w)
    @test size(∇conv_data(y, x, w; stride=2, pad=1)) == size(x)
end

@testset "maxpool2d" begin

    x = reshape(Float64[1:20;], 5, 4, 1, 1)

    @test squeeze(maxpool(x, (2,2)), (3,4)) == [7 17; 9 19]
    @test squeeze(maxpool(x, (2,2); stride=(2,2)), (3,4)) == [7 17; 9 19]
    @test squeeze(maxpool(x, (2,2); pad=(1,1)), (3,4)) == [
        1.0  11.0  16.0;
        3.0  13.0  18.0;
        5.0  15.0  20.0;
    ]

    # for gradients, check only size
    # correctness of gradients is cross-checked with CUDNN.jl
    # (it's assumed maxpooling code won't change often)

    y = maxpool(x, (2,2))
    dy = reshape(rand(2,2), 2, 2, 1, 1)
    @test size(∇maxpool(dy, y, x, (2,2))) == size(x)

end

@testset "conv3d" begin

    x = reshape(Float64[1:60;], 5, 4, 3, 1, 1)
    w = reshape(Float64[1:8;], 2, 2, 2, 1, 1)
    res = zeros(4,3,2)
    res[:, :, 1] = [
        322.0  502.0  682.0;
        358.0  538.0  718.0;
        394.0  574.0  754.0;
        430.0  610.0  790.0]
    res[:, :, 2] = [
        1042.0  1222.0  1402.0;
        1078.0  1258.0  1438.0;
        1114.0  1294.0  1474.0;
        1150.0  1330.0  1510.0]
    @test squeeze(conv(x, w),(4,5)) == res

    @test squeeze(conv(x, w; stride=2),(3,4,5)) == [
        322.0 682.0;
        394.0 754.0]

    res = zeros(6,5,4)
    res[:, :, 1] = [
        1.0   9.0   29.0   49.0   48.0;
        4.0  29.0   79.0  129.0  115.0;
        7.0  39.0   89.0  139.0  122.0;
        10.0  49.0   99.0  149.0  129.0;
        13.0  59.0  109.0  159.0  136.0;
        10.0  40.0   70.0  100.0   80.0]
    res[:, :, 2] = [
        26.0  126.0  206.0  286.0  220.0;
        80.0  322.0  502.0  682.0  502.0;
        94.0  358.0  538.0  718.0  524.0;
        108.0  394.0  574.0  754.0  546.0;
        122.0  430.0  610.0  790.0  568.0;
        80.0  260.0  360.0  460.0  320.0]
    res[:, :, 3] = [
        146.0   446.0   526.0   606.0   420.0;
        360.0  1042.0  1222.0  1402.0   942.0;
        374.0  1078.0  1258.0  1438.0   964.0;
        388.0  1114.0  1294.0  1474.0   986.0;
        402.0  1150.0  1330.0  1510.0  1008.0;
        240.0   660.0   760.0   860.0   560.0]
    res[:, :, 4] = [
        205.0   517.0   577.0   637.0  392.0;
        456.0  1133.0  1263.0  1393.0  847.0;
        467.0  1159.0  1289.0  1419.0  862.0;
        478.0  1185.0  1315.0  1445.0  877.0;
        489.0  1211.0  1341.0  1471.0  892.0;
        270.0   660.0   730.0   800.0  480.0]
    @test squeeze(conv(x, w; pad=1),(4,5)) == res

    # for gradients, check only size
    # correctness of gradients is cross-checked with CUDNN.jl
    # (it's assumed convolution code won't change often)

    @test size(∇conv_filter(reshape(rand(4,3,2), 4, 3, 2, 1, 1), x, w)) == size(w)
    @test size(∇conv_data(reshape(rand(4,3,2), 4, 3, 2, 1, 1), x, w)) == size(x)

end

@testset "maxpool3d" begin

    x = reshape(Float64[1:60;], 5, 4, 3, 1, 1)

    @test squeeze(maxpool(x, (2,2,2)), (3,4,5)) == [27 37; 29 39.]
    @test squeeze(maxpool(x, (2,2,2); stride=(2,2,2)), (3,4,5)) == [27 37; 29 39.]
    res = zeros(3,3,2)
    res[:, :, 1] = [
        1.0  11.0  16.0;
        3.0  13.0  18.0;
        5.0  15.0  20.0]
    res[:, :, 2] = [
        41.0  51.0  56.0;
        43.0  53.0  58.0;
        45.0  55.0  60.0]
    @test squeeze(maxpool(x, (2,2,2), pad=(1,1,1)), (4,5)) == res

    # for gradients, check only size
    # correctness of gradients is cross-checked with CUDNN.jl
    # (it's assumed maxpooling code won't change often)

    y = maxpool(x, (2,2,2))
    dy = reshape(rand(2,2), 2, 2, 1, 1, 1)
    @test size(∇maxpool(dy, y, x, (2,2,2))) == size(x)

end


text:
using Test, GaussNewton

function rosenblock(r::AbstractArray{Float64}, x::AbstractArray{Float64})
    r[1] = 10 * (x[2] - x[1]^2)
    r[2] = 1 - x[1]
end

function brownbad(r::AbstractArray{T}, x::AbstractArray{T}) where {T}
    r[1] = x[1] - 1000000
    r[2] = x[2] - 2e-6
    r[3] = x[1] * x[2] - 2
end

function beale(r::AbstractArray{Float64}, x::AbstractArray{Float64})
    r[1] = 1.5 - x[1] * (1 - x[2])
    r[2] = 2.25 - x[1] * (1 - x[2]^2)
    r[3] = 2.625 - x[1] * (1 - x[2]^3)
end

@testset "GNRosenblock" begin
    x0 = [-1.2, 1.0]
    obj = rosenblock
    f0, info = optimize!(obj, x0, zeros(2))
    println(f0, " ", x0, " ", info)
    @test f0 == 0
    @test isapprox(1, x0[1], atol = 1e-7)
    @test isapprox(1, x0[2], atol = 1e-7)
    @test has_converged(info) && !is_fatal(info)
end

@testset "GNbrownbad" begin
    x0 = [1.0, 1.0]
    f0, info = optimize!(brownbad, x0, zeros(3))
    println(f0, " ", x0, " ", info)
    @test abs(f0) < eps(Float64)
    @test isapprox(1e6, x0[1], atol = 1e-8)
    @test isapprox(2e-6, x0[2], atol = 1e-8)
    @test has_converged(info) && !is_fatal(info)

    x0 = [1.0, 1.0]
    f0, info = optimize!(brownbad, x0, zeros(3), autodiff = :forward)
    println(f0, " ", x0, " ", info)
    @test abs(f0) < eps(Float64)
    @test isapprox(1e6, x0[1], atol = 1e-8)
    @test isapprox(2e-6, x0[2], atol = 1e-8)
    @test has_converged(info) && !is_fatal(info)

    x0 = [1.0, 1.0]
    f0, info = optimize!(brownbad, x0, zeros(3), autodiff = :central)
    println(f0, " ", x0, " ", info)
    @test abs(f0) < eps(Float64)
    @test isapprox(1e6, x0[1], atol = 1e-8)
    @test isapprox(2e-6, x0[2], atol = 1e-8)
    @test has_converged(info) && !is_fatal(info)

end

@testset "GNBeale" begin
    x0 = [1.0, 1.0]
    f0, info = optimize!(beale, x0, zeros(3))
    println(f0, " ", x0, " ", info)
    @test abs(f0) < sqrt(eps(Float64))
    @test isapprox(3.0, x0[1], atol = 1e-7)
    @test isapprox(0.5, x0[2], atol = 1e-7)
    @test has_converged(info) && !is_fatal(info)
end


text:
function isRooted(args...; phy = thetree)
	# if typeof(phy) != Phylo
	# 	 error("isRooted:wrongtpe\n", "# object phy is not class Phylo")
	# end
	if countmap(phy["edge"][:,1])[length(phy["tip_label"]) + 1] > 2
		false
	else true
	end
end


text:
# -------------------------------------------------------------------
# # MODEL METHODS

# *models* are externally defined objects with certain functionality
# that is exposed by overloading these methods:

err_train(model) =
    ArgumentError("`IterationControl.train!(model, ::Int)` "*
                  "not overloaded "*
                  "for `typeof(model) = $(typeof(model))`. ")

err_ingest(model) =
    ArgumentError("Cannot use `Data` control here, as "*
                  "`IterationControl.ingest!(model)` "*
                  "not overloaded "*
                  "for `typeof(model) = $(typeof(model))`. ")

# ## COMPULSORY

# train `model` for another `Δn` iterations:
train!(model, Δn) = throw(err_train(model))

# ## REQUIRED FOR SOME CONTROLS

# extract a single numerical estimate of `models`'s performance such
# as an out-of-sample loss; smaller is understood to be better:
loss(model) = nothing

# extract a vector of per-iteration "training" losses, accumulated in
# the last `train!(model, Δn)` call; this will have length `Δn`:
training_losses(model) = nothing

# ingest data:
ingest!(model, data) = throw(err_ingest(model))

# ## OPTIONAL

# so that user-specified functions (in eg `Info`) are applied to some
# internal part of the model, rather than `model` directly:
expose(model) = model

# switch for expose; not to be overloaded:
expose(model, raw) = expose(model, Val(raw))
expose(model, ::Val{true}) = model
expose(model, ::Val{false}) = expose(model)

# -------------------------------------------------------------------
# # CONTROL METHODS

# compulsory: `update!`
# optional: `done`, `takedown`, `needs_training_losses`

# called after first training event; returns initialized control
# "state":
update!(control, model, verbosity, n) = nothing

# called after all subsequent training events; returns new "state":
update!(control, model, verbosity, n, state) = state

# should we stop?
done(control, state) = false

# What to do after `control`, or some other control, has triggered
# a stop. Returns user-inspectable outcomes associated with the
# control's applications (separate from logging). This should be a
# named tuple, except for composite controls which return a tuple of
# named tuples (see composite_controls.jl):
takedown(control, verbosity, state) = NamedTuple()

# If it is nonsensical to apply `control` to any model for which
# `training_losses(model)` has not been overloaded and we want an
# error thrown when this is attempted, then overload this trait to
# take the value `true`. Otherwise `control` is applied anyway, and
# `training_losses`, if called, returns `nothing`.
needs_training_losses(control) = false

# If it is nonsensical to apply `control` to any model for which
# `loss(model)` has not been overloaded and we want an error thrown
# when this is attempted, then overload this trait to take the value
# `true`. Otherwise `control` is applied anyway, and
# `loss`, if called, returns `nothing`.
needs_loss(control) = false


text:
using DataStructures
using Distances

import Base.isequal
import Base.hash

"""
  DTBSplit

implements split of KD-Tree
"""
struct DTBSplit
    dim::Int64
    val::Float64
    sgn::Bool # 0: smaller, 1 : greater
end

"""
  MyE

implements undirected edge
"""
struct MyE
    a::Int64
    b::Int64
    MyE(pa::Int64,pb::Int64) = (pa<pb) ? new(pa,pb) : new(pb,pa)
end
function isequal(a::MyE,b::MyE)
    return a.a==b.a && a.b==b.b
end
function hash(a::MyE)
    return xor( hash(a.a) , hash(a.b) )
end
function to_matrix(ee::Vector{MyE})
    m::Array{Int64,2} = Array{Int64}(length(ee),2)
    for zi=1:length(ee); m[zi,1] = ee[zi].a; m[zi,2] = ee[zi].b; end
    return m
end

struct DTBProfiling
    n_comparisons::Int64
    n_break_boxdist::Int64
    n_break_onecomp::Int64
end
function DTBProfiling()
    return DTBProfiling(0,0,0)
end

"""
  DTBNode

implements node of a KD tree
"""
mutable struct DTBNode
    id::Int64 # identifies the node ; root has 00..00 , then
              # levels (right/left) are encoded as 01 or 10
              # i.e. node root->"right"->"right"->"left" ends with:
              # 10 01 01 00
              #  l r  r  R

    # data / indidces of this kd tree node
    data::Array{Float64,2}
    subset::Array{Int64}
    #splits::Array{DTBSplit,1}

    # kd tree node bounds
    box_lb::Array{Float64,1}
    box_ub::Array{Float64,1}

    # max_{v in this node} ( length of candidate edge of v)
    dQ::Float64

    left::DTBNode
    right::DTBNode

    #DTBNode(data::Array{Float64,2}) = (x = new(data); x.left=x; x.right=x)
    #function DTBNode(data::Array{Float64,2},subset::Array{Int64},splits::Array{DTBSplit,1},dQ::Float64)
    function DTBNode(id::Int64,data::Array{Float64,2},subset::Array{Int64},box_lb::Array{Float64,1},box_ub::Array{Float64,1},dQ::Float64)
        n = new(id,data,subset,box_lb,box_ub,dQ)
        n.left=n
        n.right=n
        return n
    end
end
function isequal(a::DTBNode,b::DTBNode)
    return a.id==b.id
end
function hash(a::DTBNode)
    #a.id
    hash(a.id) # we could also drop this hash.. should not matter..
end
function is_leaf(n::DTBNode)
    return n==n
    #return ( n.id & (1<<63) ) == (1<<63)
end

"""
  kdtree(xx::Array{Float64,2})

initializes root node of KD tree (but does not compute the splits)
To build the tree use:
  kdt_root = kdtree( data )
  kdtree_split!( kdt_root , 10 ) # 10 is the node size where we stop splitting
"""
function kdtree(xx::Array{Float64,2})
    #root = DTBNode(xx,-1,NaN,nothing,nothing)
    #root = DTBNode(Int64(0),xx,collect(1:size(xx,2)),fill!(ones(size(xx,1)),-Inf),fill!(ones(size(xx,1)),Inf),Inf)
    root = DTBNode(Int64(0),xx,collect(Int64(1):Int64(size(xx,2))),fill!(ones(size(xx,1)),-Inf),fill!(ones(size(xx,1)),Inf),Inf)
    return root
end

"""
  kdtree_split!(node::DTBNode, nmin::Int64)

computes splits of KD tree. nmin indicates the node size where we stop
splitting nodes.
"""
function kdtree_split!( node::DTBNode, nmin::Int64)

    # ensure that parent node is in the tree.
    #t[node.id] = node

    if( size(node.data,2) <= nmin )
        return node
    end

    mind = minimum(node.data,2)
    maxd = maximum(node.data,2)
    ds   = findmax( maxd-mind )
    ds   = ds[2]
    #println("$ds")
    vs   = median(node.data[ds,:])
    #println("$vs")
    bx   = node.data[ds,:] .<= vs

    range_a = node.subset[bx]
    range_b = node.subset[.~bx]

    data_a = node.data[:,bx]
    data_b = node.data[:,.~bx]

    box_lb_a = copy(node.box_lb)
    box_ub_a = copy(node.box_ub)
    box_ub_a[ds] = vs
    box_lb_b = copy(node.box_lb)
    box_lb_b[ds] = vs
    box_ub_b     = copy(node.box_ub)

    id::Int64 = node.id
    id_depth  = Int( ceil( (64-leading_zeros(id)) / 2) ) + 1 # in this cell we are, i.e. we have to shift id_depth times left by two bits..
    id_l = id   |  (1)<<(2*id_depth)
    #id_l = id_l |  (1)<<(63)
    id_r = id   |  (2)<<(2*id_depth)
    #id_r = id_r |  (1)<<(63)
    #println("id $(id) -> r $(id_l) , l $(id_r)")

    node_left  = kdtree_split!( DTBNode(id_l,data_a,range_a, box_lb_a   , box_ub_a  , Inf) , nmin )
    node_right = kdtree_split!( DTBNode(id_r,data_b,range_b, box_lb_b   , box_ub_b  , Inf) , nmin )

    node.left  = node_left
    node.right = node_right
    #println("SPLIT") #println(data_a) #println(data_b)
    #println(splits_l)
    return node
end

"""
  compute_emst(data::Array{Float64,2};nmin::Int64=64)

Computes EMST for the given data (where columns are samples).
nmin is the max number of elements in kd-tree node.
"""
function compute_emst(data::Array{Float64,2};nmin::Int=64)
    nmin64 = Int64(nmin)
    root = kdtree(data)
    kdtree_split!(root,nmin64)
    edges = dtb(root,IntDisjointSets(size(data,2)))
    return EMST.to_matrix( collect(edges) )
end

"""
  dtb(q::DTBNode,e::IntDisjointSets)

implements the dual-tree Boruvka algorithm which computes the EMST for the
given KD-tree.
"""
function dtb(q::DTBNode,e::IntDisjointSets)
    edges = Set{MyE}()

    #hashmap_QvsR = Dict{Tuple{DTBNode,DTBNode},Array{Float64,2}}()

    while(e.ngroups>1)
        ngroups = e.ngroups;
        println("--> ngroups: $ngroups")

        # prepare dicts for candidate edges..
        C_dcq = Dict{Int64,Float64}()
        C_e   = Dict{Int64,MyE}()
        # init themz
        roots = unique(e.parents)
        for ri in roots
            C_dcq[ri] = Inf
        end

        # compute..
        #fcn(q,q,e,C_dcq,C_e , hashmap_QvsR )
        fcn(q,q,e,C_dcq,C_e)

        # and now add the edges..
        for ne::MyE in values(C_e)
            union!(e,ne.a,ne.b)
            push!(edges,ne)
        end
    end
    return edges
end

"""
  fcn(q::DTBNode,r::DTBNode,e::IntDisjointSets , C_dcq::Dict{Int64,Float64} , C_e::Dict{Int64,MyE} )
Find Component Neighbors

Components are identified by the root of the component
C_dcq : component distances to candidate edges
C_e   : component candidate edges (comp i candidate edge is C_e[:,i])
"""
function fcn(q::DTBNode,r::DTBNode,e::IntDisjointSets , C_dcq::Dict{Int64,Float64} , C_e::Dict{Int64,MyE} )
#function fcn(q::DTBNode,r::DTBNode,e::IntDisjointSets , C_dcq::Dict{Int64,Float64} , C_e::Dict{Int64,Tuple{Int64,Int64}} )
#function fcn(q::DTBNode,r::DTBNode,e::IntDisjointSets , C_dcq::Dict{Int64,Float64} , C_e::Dict{Int64,MyE} , hashmap_QvsR::Dict{Tuple{DTBNode,DTBNode},Array{Float64,2}})

    # check all in same component
    onecomp::Bool  = true
    joined         = [q.subset;r.subset]
    #println(joined)
    for ji in joined
        if(~in_same_set(e,joined[1],ji))
            onecomp=false
            break
        end
    end
    if(onecomp)
        return
    end

    #check d(Q,R) > d(Q)
    dqr = computeDQR( q.box_lb,q.box_ub , r.box_lb,r.box_ub )
    #println("dq : $(dqr) , q.dQ : $(q.dQ)")
    if( dqr > q.dQ)
        #println(".")
        return
    end

    #check if R and Q in a leaf node
    if( q.left==q && r.left==r ) # that's not very elegant but probably what we have to do?!?
    #if( is_leaf(q) && is_leaf(r) )
        #println("LEAF NODES")
        #n_dQ = -Inf
        n_dQ::Float64 = q.dQ

        all_d_qr = Distances.pairwise(Euclidean(),q.data,r.data)

        # --> hashing these values does not really speed up things..
        #     and it uses tons of memory, such that fancier handling of
        #     hash map entries would be required
        #if(!haskey(hashmap_QvsR,(q,r)))
        #    if(!haskey(hashmap_QvsR,(r,q)))
        #        all_d_qr = Distances.pairwise(Euclidean(),q.data,r.data)
        #        #hashmap_QvsR[(q,r)] = all_d_qr
        #    else
        #        all_d_qr = transpose(hashmap_QvsR[(r,q)])
        #    end
        #    hashmap_QvsR[(r,q)] = transpose(all_d_qr)
        #else
        #    all_d_qr = hashmap_QvsR[(q,r)]
        #    #print('.')
        #end

        for iq = 1:size(q.subset,1)
            for ir = 1:size(r.subset,1)
                qq = q.subset[iq]
                rr = r.subset[ir]
                if( in_same_set(e , qq , rr ) )
                    continue
                end

                cq = find_root(e,rr) # compoment of q

                # check distance:
                # dist_qr::Float64 = sqrt( sum( (q.data[:,iq]-r.data[:,ir]).^2 ) )
                dist_qr = all_d_qr[iq,ir]

                if( dist_qr < C_dcq[ cq ] )
                    C_dcq[ cq ] = dist_qr
                    C_e[ cq ]   = MyE(qq,rr) #(qq,rr)
                    # and dQ !
                    n_dQ = max(n_dQ,dist_qr)
                    #n_dQ = dist_qr
                end
            end
        end
        q.dQ = n_dQ

        return
    end

    # do recursions..
    #fcn( q.left,r.left,e ,  C_dcq::Dict{Int64,Float64} , C_e::Dict{Int64,Tuple{Int64,Int64}} )
    #q_left  = is_leaf(q) ? q : t[q.left]
    #q_right = is_leaf(q) ? q : t[q.right]
    #r_left  = is_leaf(r) ? r : t[r.left]
    #r_right = is_leaf(r) ? r : t[r.right]
    fcn( q.left,r.left    , e ,  C_dcq , C_e ) #, hashmap_QvsR )
    fcn( q.right,r.left   , e ,  C_dcq , C_e ) #, hashmap_QvsR )
    fcn( q.left,r.right   , e ,  C_dcq , C_e ) #, hashmap_QvsR )
    fcn( q.right,r.right  , e ,  C_dcq , C_e ) #, hashmap_QvsR )
    q.dQ = max( q.left.dQ , q.right.dQ )
end

"""
  computeDQR( q_lb::Array{Float64,1} , q_ub::Array{Float64,1} , r_lb::Array{Float64,1} , r_ub::Array{Float64,1} )

compute min dist. between bounding boxes, i.e. between rectangular boxes Q/R with
bounds given by q_lb/q_ub and r_lb/r_ub
"""
function computeDQR( q_lb::Array{Float64,1} , q_ub::Array{Float64,1} , r_lb::Array{Float64,1} , r_ub::Array{Float64,1} )
    d::Int64  = length(q_lb)#size(q_lb,1)
    rdists::Array{Float64}    = zeros(d)
    for zd = 1:d
        rdists[zd] = max(  max( q_lb[zd] - r_ub[zd] , r_lb[zd] - q_ub[zd] ) , 0 )
    end
    dqr::Float64 = sqrt( sum( rdists.^2 ))
    return dqr
end


text:
export
    DriverModel,

    StaticDriver,

    get_name,
    action_type,
    set_desired_speed!,
    observe!,
    reset_hidden_state!,
    prime_with_history!

abstract type DriverModel{DriveAction} end

get_name(::DriverModel) = "???"
action_type(::DriverModel{A}) where {A} = A
set_desired_speed!(model::DriverModel, v_des::Float64) = model # do nothing by default
reset_hidden_state!(model::DriverModel) = model # do nothing by default
observe!(model::DriverModel, scene::EntityFrame{S,D,I}, roadway::R, egoid::Int) where {S,D,I,R} = model  # do nothing by default
Base.rand(model::DriverModel) = error("rand not implemented for model $model")
Distributions.pdf(model::DriverModel{A}, a::A) where {A} = error("pdf not implemented for model $model")
Distributions.logpdf(model::DriverModel{A}, a::A) where {A} = error("logpdf not implemented for model $model")

function prime_with_history!(
    model::DriverModel,
    trajdata::ListRecord{Entity{S,D,I}},
    roadway::R,
    frame_start::Int,
    frame_end::Int,
    egoid::I,
    scene::EntityFrame{S,D,I} = allocate_frame(trajdata),
    ) where {S,D,I,R}

    reset_hidden_state!(model)

    for frame in frame_start : frame_end
        get!(scene, trajdata, frame)
        observe!(model, scene, roadway, egoid)
    end

    return model
end
function prime_with_history!(model::DriverModel, rec::EntityQueueRecord{S,D,I}, roadway::R, egoid::I;
    pastframe_start::Int=1-nframes(rec),
    pastframe_end::Int=0,
    ) where {S,D,I,R}

    reset_hidden_state!(model)

    for pastframe in pastframe_start : pastframe_end
        scene = rec[pastframe]
        observe!(model, scene, roadway, egoid)
    end

    model
end

####

struct StaticDriver{A,P<:ContinuousMultivariateDistribution} <: DriverModel{A}
    distribution::P
end

get_name(::StaticDriver) = "StaticDriver"
function Base.rand(model::StaticDriver{A,P}) where {A,P}
    a = rand(model.distribution)
    return convert(A, a)
end
Distributions.pdf(model::StaticDriver{A}, a::A) where {A} = pdf(model.distribution, convert(Vector{Float64}, a))
Distributions.logpdf(model::StaticDriver{A}, a::A) where {A} = logpdf(model.distribution, convert(Vector{Float64}, a))

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
# high level user API

# two major visualization functions are used, namely:
# - Memory direction with visualize(::FactorGraph, ...)
# - Cloud access with visualize(robot::String, session::String, ...)

"""
    $(SIGNATURES)

Set the global transform for modifying the final draw step, for example Z up or Z down convention.

Example:
```julia
using Arena
# default is Z up

# make Z down
setGlobalDrawTransform!(quat=Quat(0.0,1.0,0.0,0.0))
```
"""
function setGlobalDrawTransform!(;trans=Translation(0.0,0.0,0.0), quat::Rotations.Quat=Quat(1.0,0.0,0.0,0.0))
  global drawtransform = trans ∘ LinearMap(quat)
end

"""
    $(SIGNATURES)

Toggle global boolean flag to terminate the visualization loop.
"""
function stopVis!()
  global loopvis
  loopvis = false
  nothing
end

function extractRobotSession!(rose_fgl::Union{FactorGraph, Tuple{<:AbstractString, <:AbstractString}},
                              params::Dict{Symbol, Any}  )

    config = nothing

    # robot and session defaults
    params[:robotId] = :Bot
    params[:sessionId] = :Session
    if isa(rose_fgl, FactorGraph)
        if rose_fgl.robotname != "" && rose_fgl.robotname != "NA"
            params[:robotId] = Symbol(rose_fgl.robotname)
        end
        if rose_fgl.sessionname != "" && rose_fgl.sessionname != "NA"
            params[:sessionId] = Symbol(rose_fgl.sessionname)
        end
    else
        # setup Graff config (iff a Graff visualization)
        getGraffConfig() == nothing && error("GraffSDK is not configured, please call loadGraffConfig() to set up Graff.")
        config = getGraffConfig()
        params[:robotId] = Symbol(rose_fgl[1])
        params[:sessionId] = Symbol(rose_fgl[2])
        config.robotId = string(params[:robotId])
        config.sessionId = string(params[:sessionId])
    end
    return config
end

"""
    $(SIGNATURES)

High level interface to launch webserver process that draws the factor graph contents using Three.js and MeshCat.jl.

Examples:
---------

```julia
# start webserver visible anywhere on the network on port 8000.
using GraffSDK, Arena
@async visualize( (robot, session), meanmax=:max )
```

Or view local memory `FactorGraph` version:
```julia
# also open a browser to connect to the three.js webserver
using Arena
@async visualize(fg, show=true, plugins=[myplugin1;])
```

See src/plugins/Template.jl for defining your own plugin to be rendered in the visualization loop.
```julia
function pluginCallback(vis::MeshCat.Visualizer, params::Dict{String, Any}, rose_fgl) ... end
```

**Note** plugins do not have to be defined in Arena, Main context or users repo is sufficient.
"""
function visualize(rose_fgl::Union{FactorGraph, Tuple{<:AbstractString, <:AbstractString}};
                   show::Bool=true,
                   meanmax=:max,
                   plugins::Vector=Function[]  )::Nothing
    #
    global loopvis
    global drawtransform

    loopvis = true

    # standard parameters dictionary
    params = Dict{Symbol, Any}()

    # set up basics in params and load Graff config if available
    extractRobotSession!(rose_fgl, params)

    # the visualizer object itself
    vis = startDefaultVisualization(show=show)

    # default variable caching format:  (softtype, already-drawn, mapEst)
    params[:cachevars] = Dict{Symbol, Tuple{Symbol, Vector{Bool}, Vector{Float64}}}()

    # Prepend default plugins to be executed
    plugins = union([cacheVariablePointEst!; visualizeVariableCache!], plugins)

    # run the visualization loop
    while loopvis
        # iterate through all listed plugin callbacks, such sa pointclouds / images / reprojections / etc.
        for callback in plugins
            try
                callback(vis, params, rose_fgl)
            catch ex
                io = IOBuffer()
                showerror(io, ex, catch_backtrace())
                err = String(take!(io))
                @error "Visualization plugin failure -- callback=$(string(callback)) errored: $ex"
                @error "Error! $err"
            end
        end

        # take a break and repeat
        sleep(1)
    end

    @info "visualize is finalizing."
    nothing
end

### Some old code -- will be deleted soon enough

# xx, ll = ls(fgl)
# for x in xx
#   X = getKDE(fgl, x)
#   xmx = meanmax == :max ? getKDEMax(X) : getKDEMean(X)
#
#   if !haskey(cachevars, x)
#     triad = Triad(1.0)
#     cachevars[x] = deepcopy(xmx)
#   end
#
#     setobject!(vis[:poses][x], triad)
#
#   cachevars[x][:] .= xmx
#   trans = Translation(xmx[1:2]..., 0.0) ∘ LinearMap(RotZ(xmx[3]))
#   settransform!(vis[:poses][x], drawtransform ∘ trans)
# end


text:
function write_stl(dataset::AbstractStaticVTKData, filepath_no_ext::String)
    polydata = PyVTK(trianulate(dataset))
    return write_stl(polydata, filepath_no_ext)
end

function write_stl(polydata::PyCall.PyObject, filepath_no_ext::String)
    writer = vtk.vtkSTLWriter()
    writer.SetFileName(string(filepath_no_ext, ".stl"))
    writer.SetInputData(polydata)
    writer.Update()
    return string(filepath_no_ext, ".stl")
end

function read_stl(filepath::String)
    if splitext(filepath)[2] == ".stl"
        reader = vtk.vtkSTLReader()
        reader.SetFileName(filepath)
        reader.Update()

        polydata = reader.GetOutput()
        return extract_simple_block(polydata)
    else
        throw("The file input is not an stl file.")
    end
end


text:
using DelimitedFiles

include("../../../src/transit_poly_struct.jl")

function profile_transit_poly!(trans,flux,b)
    t1 = time_ns()
    for k=1:10
        for i=1:length(b)
          trans.b = b[i]
          flux[i] = transit_poly!(trans)
        end
    end
    elapsed = 1e-10 * (time_ns() - t1)
    return elapsed
end

b = readdlm("b.txt")
u = readdlm("u.txt")
u_n = zeros(length(u))
for i = 1:length(u)
    u_n[i] = u[i]
end
trans = transit_init(0.1, 0.0, u_n, false)
flux = zeros(length(b))
println(profile_transit_poly!(trans, flux, b))
writedlm("flux.txt", flux)

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
module StructureLearning
    include("BaseStructureLearning.jl")
    include("HillClimb.jl")
end  # module StructureLearning


text:
cod2amino = Dict( "ATA" => 8, "ATC" => 8, "ATT"=>8, "ATG"=> 11,
        "ACA"=>17, "ACC"=>17, "ACG"=>17, "ACT"=> 17,
        "AAC"=>12, "AAT"=>12, "AAA"=>9, "AAG"=>9,
        "AGC"=>16, "AGT"=> 16, "AGA"=> 15, "AGG"=> 15,
        "CTA"=>10, "CTC"=>10, "CTG"=>10, "CTT"=>10,
        "CCA"=>13, "CCC"=>13, "CCG"=>13, "CCT"=>13,
        "CAC"=>7, "CAT"=>7, "CAA"=>14, "CAG"=>14,
        "CGA"=>15, "CGC"=>15, "CGG"=>15, "CGT"=>15,
        "GTA"=>18, "GTC"=>18, "GTG"=>18, "GTT"=>18,
        "GCA"=>1, "GCC"=>1, "GCG"=>1, "GCT"=>1,
        "GAC"=>3, "GAT"=>3, "GAA"=>4, "GAG"=>4,
        "GGA"=>6, "GGC"=>6, "GGG"=>6, "GGT"=>6,
        "TCA"=>16, "TCC"=>16, "TCG"=>16, "TCT"=>16,
        "TTC"=>5, "TTT"=>5, "TTA"=>10, "TTG"=>10,
        "TAC"=>20, "TAT"=>20, "TAA"=> 21, "TAG"=> 21,
        "TGC"=> 2, "TGT"=>2 , "TGA"=> 21, "TGG"=> 19,
        "---" => 21)

amino2cod = Dict()
for amino in 1:21
    codons = Vector{String}()
    for (key, val) in cod2amino
       val == amino && push!(codons, key)
    end
    amino2cod[amino] = codons
end

mutable struct SeqToEvolve
    Amino   :: Array{Int64}
    DNA :: Array{AbstractString}
end

function cond_proba_DNA_gibbs(k, amino_list, mutated_seq, h, J,N,  T = 1)
	prob = zeros(length(amino_list))
	for (index, q_k) in enumerate(amino_list)
		log_proba = h[q_k, k]
 		for i in 1:N
			log_proba += J[mutated_seq[i], q_k ,i, k]
		end
		prob[index] = exp(log_proba/T)
	end
	return normalize(prob,1)
end

function evol_seq_fix_steps_DNA_gibbs(ref_seq, MC_steps, h, J, N, T = 1)
	mutated_seq = deepcopy(ref_seq)
	non_gapped_pos = [pos for (pos, amino) in enumerate(ref_seq.Amino) if amino != 21]
	@inbounds for steps in 1: MC_steps
        pos_mut = rand(non_gapped_pos)

		old_codon = mutated_seq.DNA[pos_mut]

		amino_list, codon_list = get_accessible_muts_DNA_gibbs(old_codon)

		aa = sample(amino_list, ProbabilityWeights(cond_proba_DNA_gibbs(pos_mut, amino_list, mutated_seq.Amino, h, J,N,  T)))

		new_codon = rand(  [codon for codon in codon_list if (cod2amino[codon] == aa)]   )

		mutated_seq.DNA[pos_mut] = new_codon
		mutated_seq.Amino[pos_mut] = aa
	end
	return mutated_seq.Amino
end

function get_accessible_muts_DNA_gibbs(old_codon)
	old_codon = [string(old_codon[i]) for i in 1:3 ]
	codon_list = Vector{AbstractString}(undef, 12)
	for i in 1:3
		new_codon = deepcopy(old_codon)
		for (j, nucl) in enumerate(["A", "C", "G", "T"])
			new_codon[i] = nucl
			codon_list[ (i-1)*4 + j ] = join(new_codon)
		end
	end

	amino_list = get.(Ref(cod2amino),codon_list, 0)
	amino_list = unique!(filter!(aa->aa != 21, amino_list))

	return amino_list, codon_list
end

####
#MAIN FUNCTION
####

function evolMSA(output_path::AbstractString, params_path::AbstractString, wt_path::AbstractString; steps::Integer = 10, nseq::Integer = 100, T::Real = 1, wt_name::AbstractString = "unknown wt")
	for file in [params_path, wt_path]
	    !isfile(file) && error("Error: the file \"$(file)\" does not exist. Please check the spelling or the folder path.")
	end

	steps < 1 && throw(DomainError(steps, "'steps' must be a positive integer."))
	nseq < 1 && throw(DomainError(nseq, "'nseq' must be a positive integer."))
	T <= 0 && throw(DomainError(T, "'T' must be a positive real number."))

	h, J = extract_params(params_path)

	seq = join(readdlm(wt_path, skipstart = 1))
	L = Int64(length(seq)/3)
	DNA_seq = [seq[((i-1)*3 +1):(i*3)] for i in 1:L]
	amino_seq = [cod2amino[codon] for codon in DNA_seq]
	seed_seq = SeqToEvolve(amino_seq, DNA_seq)
	N = length(seed_seq.Amino)
	FastaWriter(output_path, "a") do file
		for i in 1:nseq
			seq = evol_seq_fix_steps_DNA_gibbs(seed_seq, steps, h, J, N, T)
			writeentry(file, "$i | original wt: $wt_name | $steps MC steps | T = $(T)", vec2string(seq))
		end
	end
end

function evolMSA(output_path::AbstractString, params::Tuple{Array{Float64, 2}, Array{Float64, 4}}, wt_path::AbstractString; steps::Integer = 10, nseq::Integer = 100, T::Real = 1, wt_name::AbstractString = "unknown wt")
	for file in [wt_path]
	    !isfile(file) && error("Error: the file \"$(file)\" does not exist. Please check the spelling or the folder path.")
	end

	steps < 1 && throw(DomainError(steps, "'steps' must be a positive integer."))
	nseq < 1 && throw(DomainError(nseq, "'nseq' must be a positive integer."))
	T <= 0 && throw(DomainError(T, "'T' must be a positive real number."))

	h, J = params

	seq = join(readdlm(wt_path, skipstart = 1))
	L = Int64(length(seq)/3)
	DNA_seq = [seq[((i-1)*3 +1):(i*3)] for i in 1:L]
	amino_seq = [cod2amino[codon] for codon in DNA_seq]
	seed_seq = SeqToEvolve(amino_seq, DNA_seq)
	N = length(seed_seq.Amino)
	FastaWriter(output_path, "a") do file
		for i in 1:nseq
			seq = evol_seq_fix_steps_DNA_gibbs(seed_seq, steps, h, J, N, T)
			writeentry(file, "$i | original wt: $wt_name | $steps MC steps | T = $(T)", vec2string(seq))
		end
	end
end

max_stars_repo_name: matteobisardi/SeqEvol

max_issues_repo_name: matteobisardi/SeqEvol

max_forks_repo_name: matteobisardi/SeqEvol


text:
"""
    domain_decomposition(N::Int64, N_workers::Int64)

Calculate relevant array slices for each worker.
Could be done better!
"""
function domain_decomposition(N::Integer, N_workers::Integer)

    batch = Array{typeof(1:2)}(undef, N_workers)
    size = Int(floor(N/N_workers))
    @inbounds for i = 1:N_workers-1
        batch[i] = 1+(i-1)*size:i*size
    end
    batch[N_workers] = 1 + (N_workers-1)*size:N

    return batch
end


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra
using OrdinaryDiffEq
using DiffEqFlux, Flux

using JLD2, FileIO
using Statistics, StatsBase
# Set a random seed for reproduceable behaviour
using StatsPlots
using Plots
gr()

function retrieve_results(i::Int64, fname::String = "Scenario_1_recovery_loop.jld2")
    try
        jldopen(fname, "r") do file
            res = file["$i"]["result"]
            losses = file["$i"]["losses"]
            return res, losses
        end
    catch
        @info "Failed on $i"
        return
    end
end

function check_eqs(Ψ::SparseIdentificationResult)
    sum(get_sparsity(Ψ)) != 2 && return false
    ps = parameters(Ψ.equations)
    xs = variables(Ψ.equations)
    return all(isequal.([ps[1]*xs[1]*xs[2]; ps[2]*xs[1]*xs[2]], [x.rhs for x in Ψ.equations.eqs]))
end

function collect_results(rng = 1:1:500)
    founds = zeros(Bool, length(rng))
    sparsities = zeros(Float64, length(rng))
    errors = zeros(Float64, length(rng))
    aiccs = zeros(Float64, length(rng))
    losses = Array{AbstractVector}(undef, length(rng))
    failed = zeros(Bool, length(rng))
    for (j,i) in enumerate(rng)
        exp_result = retrieve_results(i)
        if isnothing(exp_result)
            sparsities[j] = Inf
            aiccs[j] = Inf
            errors[j] = Inf
            losses[j] = [Inf]
            failed[j] = true
            continue
        end
        sparsities[j] = sum(get_sparsity(exp_result[1]))
        errors[j] = norm(get_error(exp_result[1]), 2)
        aiccs[j] = norm(get_aicc(exp_result[1]), 2)
        founds[j] = check_eqs(exp_result[1])
        losses[j] = exp_result[2]
    end
    return founds, sparsities, errors, aiccs, losses, failed
end

function return_trajectory_and_parameter(i::Int64, fname::String = "Scenario_1_recovery_loop.jld2")
    try
        jldopen(fname, "r") do file
            X = file["$i"]["X"]
            t = file["$i"]["t"]
            ps = file["$i"]["trained_parameters"]
            return X, t, ps
        end
    catch
        @info "Failed on $i"
        return
    end
end

rbf(x) = exp.(-(x.^2))

function build_estimate(i::Int64,fname::String = "Scenario_1_recovery_loop.jld2")
    X, t, ps = return_trajectory_and_parameter(i, fname)
    p_ = Float32[1.3, 0.9, 0.8, 1.8]
    U = FastChain(
        FastDense(2, 5, rbf),
        FastDense(5, 5, rbf),
        FastDense(5, 5, rbf),
        FastDense(5, 2),
    )

    # Define the hybrid model
    function ude_dynamics!(du, u, p, t)
        û = U(u, p) # Network prediction
        du[1] = p_[1] * u[1] + û[1]
        du[2] = -p_[4] * u[2] + û[2]
    end

    prob_nn = ODEProblem(ude_dynamics!, X[:, 1], (t[1], t[end]), ps)
    sol = solve(prob_nn, Tsit5(), saveat = t)
    return X, t, sol
end

## Collect for each noise lvl
res_1 = collect_results(1:1:100)
res_2 = collect_results(101:1:200)
res_3 = collect_results(201:1:300)
res_4 = collect_results(301:1:400)
res_5 = collect_results(401:1:500)
# As a matrix
founds = [res_1[1] res_2[1] res_3[1] res_4[1] res_5[1]]
sparsities = [res_1[2] res_2[2] res_3[2] res_4[2] res_5[2]]
errors = [res_1[3] res_2[3] res_3[3] res_4[3] res_5[3]]
aiccs = [res_1[4] res_2[4] res_3[4] res_4[4] res_5[4]]
losses =  [res_1[5] res_2[5] res_3[5] res_4[5] res_5[5]]
failed = [res_1[6] res_2[6] res_3[6] res_4[6] res_5[6]]

sum(failed)
## Get the colorpalette
cpal = palette(:Dark2_5)

## Found Equations plot
p_founds = bar([sum(fi) for fi in eachcol(founds)], ylim = (0, 100), label = nothing,
    xticks = (1:1:5, [1e-3, 5e-3, 1e-2, 2.5e-2, 5e-2]),
    #color = [:green, :seagreen, :limegreen, :springgreen, :aquamarine],
    color = [cpal[i] for i in 1:5],
    xlabel = "Noise Lvl.", ylabel = "Recovered Equations")

savefig(p_founds, joinpath(pwd(), "plots", "Found_Equations_Loop.pdf"))

[sum(fi) for fi in eachcol(founds)] / 100
mean([sum(fi) for fi in eachcol(founds)])
std([sum(fi) for fi in eachcol(founds)])
## Complexity Plot

#marginalkde(errors[founds], aiccs[founds], xlabel = "L2-Error", ylabel = "AICC",
#    label = ["Noise Lvl. $i" for i in [1e-3 5e-3 1e-2 2.5e-2 5e-2]],
#    #color = [:green :seagreen :limegreen :springgreen :aquamarine], alpha = 0.3,
#    palette = :Dark2_5, alpha = 0.5,)

e_s = [e[f] for (e,f) in zip(eachcol(errors), eachcol(founds))]
a_s = [e[f] for (e,f) in zip(eachcol(aiccs), eachcol(founds))]

scatter(e_s, a_s, xlabel = "L2-Error", ylabel = "AICC",
        label = ["Noise Lvl. $i" for i in [1e-3 5e-3 1e-2 2.5e-2 5e-2]],
        #color = [:green :seagreen :limegreen :springgreen :aquamarine], alpha = 0.3,
        palette = :Dark2_5, alpha = 0.5,)

boxplot(e_s, palette = :Dark2_5, label = nothing, xlabel = "Noise Lvl.", ylabel = "L2-Error",
    xticks = (1:1:5, [1e-3, 5e-3, 1e-2, 2.5e-2, 5e-2]),)
boxplot(a_s, palette = :Dark2_5, label = nothing, xlabel = "Noise Lvl.", ylabel = "AICC",
    xticks = (1:1:5, [1e-3, 5e-3, 1e-2, 2.5e-2, 5e-2]),)

## Loss Plot
function loss_plot(losses)
    p = plot(ylabel = "Loss", xlabel = "Iterations")
    vline!([200], color = :black, lw = 1, style = :dot, label = nothing)

    lsucc = false
    lfailed = false

    for (i,l) in enumerate(losses)
        if failed[i]
            continue
        end
        indicator = founds[i] ? :green : :red

        if !lsucc && founds[i]
            label = "Successful"
            lsucc = true
        elseif !lfailed && !founds[i]
            label = "Failed"
            lfailed = true
        else
            label = nothing
        end
        scatter!([length(l)], [l[end]], color = indicator, alpha = 0.5, label = label, msize = 4.)
        plot!(l, color = :blue, alpha = 0.1, yaxis = :log10, label = nothing)
    end
    finite_losses = [findmax(l) for l in losses if length(l) > 1]
    lmax, _ = findmax(finite_losses)
    #display(p)
    annotate!([(100, lmax[1]*1.5, text("ADAM", :center))])
    annotate!([(800, lmax[1]*1.5, text("BFGS", :center))])
    return p
end

plosses = loss_plot(losses)
minlength, _ = findmin(length.(l for l in losses[founds]))
l_mins = hcat([l[1:minlength] for l in losses[founds]])
l_mean = mean(l_mins, dims = 2)
plot!(l_mean, color = :red, label = "Mean Loss")
savefig(plosses, joinpath(pwd(), "plots", "Losses.pdf"))
## Sample some results
neg_samples = collect(1:500)[vcat(.!founds...)]
neg_res = sample(neg_samples, 8,replace = false)
function plot_examples(inds, samplesize = 8, idxs = 1:500)
   _samples = sample(collect(idxs)[vcat(inds...)], samplesize, replace = false, ordered = true)
   p = []
   for s in _samples
       r_  = build_estimate(s)
       if isnothing(r_)
           continue
       end
       X, t, sol = r_
       p_ = plot(sol, color = [:red :blue],label = nothing)
       scatter!(t, X', color = [:red :blue], label = nothing)
       annotate!([(2, 4.3, text("Sample $s", :left))])
       push!(p, p_)
   end
   return p
end

n_pos = plot_examples(.!founds[end-200:10:end], 8, 300:10:500)
n_examples = plot(n_pos..., layout = (4,2), size = (800, 1000))
savefig(n_examples, joinpath(pwd(), "plots", "Failed_Examples.pdf"))
p_pos = plot_examples(founds[end-300:end], 8, 200:500)
p_examples = plot(p_pos..., layout = (4,2), size = (800, 1000))
savefig(p_examples, joinpath(pwd(), "plots", "Sucessful_Examples.pdf"))


text:
# kuramoto_7-24-19.jl
#
# Back to work.

@time begin
    using Plots
    using DifferentialEquations

    function KuramotoModel(N, K, omega)
        @assert(N==size(omega)[1], "number of oscillators != N")

        function Kuramoto(dtheta, theta, p, t)
            for i in 1:N
                theta[i] = mod2pi(theta[i])

                stabilizers = map(theta_j -> sin(theta_j - theta[i]), theta)
                stabilizing_term = (K / N) * sum(stabilizers)

                dtheta[i] = omega[i] + stabilizing_term

                theta[i] = mod2pi(theta[i])
            end
        end

        return Kuramoto
    end
end
# This code takes about 30 seconds to run on my machine. Oof.

#######################################################################
# Simulate a two-oscillator Kuramoto model, over 60 seconds.

tspan = (0.0, 60.0)

K2 = KuramotoModel(2, 1.0, [1.0, 1.0])

theta0 = [0.0, pi]
prob = ODEProblem(K2, theta0, tspan)

@time sol = solve(prob, Euler(), dt=0.001) # Yes, Euler() is important.
@time plot(sol)

# Success!

# But how long does it take?
# sol = solve(...) takes around 4.2 seconds on my machine, and
# the plotting takes a full 17.7 seconds JESUS CHRIST

# I'm going to see if I can improve this a bit over in the
# techniques-to-increase-speed folder.

#######################################################################
# Simulate a 12-oscillator Kuramoto model, all with natural frequency 1
# radians / second, over 10 seconds.

tspan = (0.0, 10.0)

N = 12;
K = 1.0;

omega_n = Array{Float64}(undef, N)
for i in 1:N
    omega_n[i] = 1.0
end

theta_0 = Array{Float64}(undef, N)
for j in 1:N
    theta_0[j] = mod2pi(j)
        # This is just an easy hack to space these out. ∀ n_1, n_2 ∈ N :
        # n_1 != n_2 ⟹ mod2pi(n_1) != mod2pi(n_2). I got that from...
        # Somewhere. Trust me.
end

K12 = KuramotoModel(N, 1.0, omega_n)

prob = ODEProblem(K12, theta_0, tspan)

@time sol = solve(prob, Euler(), dt=0.0001) # Yes, Euler() is important.
@time plot(sol)
#######################################################################
# How about a *thousand*? 😉

tspan = (0.0, 10.0)

N = 1000;
K = 1.0;

omega_n = Array{Float64}(undef, N)
for i in 1:N
    omega_n[i] = 1.0
end

theta_0 = Array{Float64}(undef, N)
for j in 1:N
    theta_0[j] = mod2pi(j)
        # This is just an easy hack to space these out. ∀ n_1, n_2 ∈ N :
        # n_1 != n_2 ⟹ mod2pi(n_1) != mod2pi(n_2). I got that from...
        # Somewhere. Trust me.
end

K1000 = KuramotoModel(N, 1.0, omega_n)

prob = ODEProblem(K1000, theta_0, tspan)

sol = solve(prob, Euler(), dt=0.0001) # Yes, Euler() is important.
plot(sol)


text:
using Test
using BenchmarkTools

@testset "simpler implementations" begin

    # bisection
    xrt = Roots.bisection(sin, 3.0, 4.0)
    @test isapprox(xrt, pi)

    xrt = Roots.bisection(sin, 3.0, 4.0, xatol=1e-3)
    @test abs(sin(xrt)) >= 1e-7  # not to0 close

    xrt = Roots.bisection(sin, big(3.0), big(4.0))
    @test isapprox(xrt, pi)

    # secant_method
    fpoly = x -> x^5 - x - 1
    xrt = Roots.secant_method(fpoly, 1.0)
    @test abs(fpoly(xrt)) <= 1e-15

    xrt = Roots.secant_method(fpoly, (1, 2))
    @test abs(fpoly(xrt)) <= 1e-14

    # muller
    fpoly = x -> x^5 - x - 1
    xrt = Roots.muller(fpoly, 1.0)
    @test xrt isa Real
    @test abs(fpoly(xrt)) <= 1e-15

    @test_throws DomainError Roots.muller(fpoly, -1.0)
    xrt = Roots.muller(fpoly, -1.0+0im)
    @test xrt isa Complex
    @test abs(fpoly(xrt)) <= 1e-15

    @test Roots.muller(cos, 1.0) ≈ π/2
    expoly(z) = log(-z)*asin(z)/tanh(z)

    @test Roots.muller(expoly, -0.7-0.5im) ≈ -1.0

    # dfree
    fpoly = x -> x^5 - x - 1
    xrt = Roots.dfree(fpoly, 1.0)
    @test abs(fpoly(xrt)) <= 1e-14

    # newton
    @test Roots.newton((sin, cos), 3.0)  ≈ pi
    u = Roots.newton(x -> (sin(x), sin(x)/cos(x)), 3.0, xatol=1e-10, xrtol=1e-10)
    @test abs(u -pi) <= 1e-8
end


text:
using Benchmark

################################################################################
# A benchmark to show how broadcast_getindex performs with respect to plain
# loop. It turns out to be a bit faster. @inbounds annotation helps a little
# bit here.
#
# | Row | Function                    | Average     | Relative | Replications |
# |-----|-----------------------------|-------------|----------|--------------|
# | 1   | "mc_logistic_iter"          | 0.000492798 | 3.62449  | 100          |
# | 2   | "mc_logistic_iter_inbounds" | 0.000239364 | 1.7605   | 100          |
# | 3   | "mc_logistic_broadcast"     | 0.000135963 | 1.0      | 100          |
################################################################################

function mc_logistic_broadcast(pred, label)
  width, height, channels, num = size(pred)
  loss = sum(-log(broadcast_getindex(pred,
      reshape(1:width, (width,1,1,1)),
      reshape(1:height, (1,height,1,1)),
      int(label) + 1,
      reshape(1:num, (1,1,1,num)))))
  return loss / (width*height*num)
end

function mc_logistic_iter(pred, label)
  width, height, channels, num = size(pred)
  loss = 0
  for w = 1:width
    for h = 1:height
      for n = 1:num
        idx = int(label[w,h,1,n])+1
        loss += -log(pred[w,h,idx,n])
      end
    end
  end
  return loss / (width*height*num)
end

function mc_logistic_iter_inbounds(pred, label)
  width, height, channels, num = size(pred)
  loss = 0
  for w = 1:width
    for h = 1:height
      @simd for n = 1:num
        @inbounds loss += -log(pred[w,h,int(label[w,h,1,n])+1,n])
      end
    end
  end
  return loss / (width*height*num)
end

width, height, channels, num = (5, 6, 30, 128)

prob = abs(rand(width, height, channels, num))
label = abs(rand(Int, (width, height, 1, num))) % channels
label = convert(Array{Float64}, label)

mc_logistic_iter() = mc_logistic_iter(prob, label)
mc_logistic_broadcast() = mc_logistic_broadcast(prob, label)
mc_logistic_iter_inbounds() = mc_logistic_iter_inbounds(prob, label)

@assert abs(mc_logistic_iter() - mc_logistic_broadcast()) < 1e-10
@assert abs(mc_logistic_iter_inbounds() - mc_logistic_broadcast()) < 1e-10

df = compare([mc_logistic_iter, mc_logistic_iter_inbounds, mc_logistic_broadcast], 100)
println("$df")


text:
instructions = []
open("day08/day08.in") do io
    global instructions
    for line in readlines(io)
        line = split(line, " ")
        push!(instructions, [line[1], parse(Int64, line[2])])
    end
end

function code_runner(checkpoints, num)
    global instructions
    executed = zeros(Bool, length(instructions))
    exec_pos = 1
    accumulator = 0
    while exec_pos <= length(instructions) && !executed[exec_pos]
        executed[exec_pos] = true
        instruction = instructions[exec_pos]
        exec_pos += 1
        if instruction[1] == "acc"
            accumulator += instruction[2]
        elseif instruction[1] == "jmp"
            exec_pos += instruction[2]-1
            if num == -1
                push!(checkpoints, exec_pos)
            elseif checkpoints[num] == exec_pos
                exec_pos -= instruction[2]-1
            end
        elseif instruction[1] == "nop"
            if num == -1
                push!(checkpoints, exec_pos)
            elseif checkpoints[num] == exec_pos
                exec_pos += instruction[2]-1
            end
        end
    end
    if num == -1
        return checkpoints
    else
        return accumulator, exec_pos == length(instructions)+1
    end
end

checkpoints = code_runner([], -1)
checkpoint_pos = length(checkpoints)
succ = false
accumulator = 0
while !succ
    global checkpoints, succ, checkpoint_pos, accumulator
    accumulator, succ = code_runner(checkpoints, checkpoint_pos)
    checkpoint_pos -= 1
end

println(accumulator)

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
using JLD, FLoops

function calc(c)
    z = c
    for _ in 1:15000
        z = z^2 + c
        abs(z) > 2 && return true
    end
    false
end

function main(dim)::Vector{ComplexF64}
    @floop for b in -1.5:3/dim:1.5
        for a in -2.:3/dim:1.
            c = complex(a, b)
            calc(c) && @reduce(points = append!([], c))
        end
        println(b)
    end
    println(length(points))
    points
end

@time save("./out/20k.jld", "points", main(20000))

max_stars_repo_name: BoardTM/fractals

max_issues_repo_name: BoardTM/fractals

max_forks_repo_name: BoardTM/fractals


text:
### BasicContMuvParameter

type BasicContMuvParameter <: Parameter{Continuous, Multivariate}
  key::Symbol
  index::Integer
  pdf::Union{ContinuousMultivariateDistribution, Void}
  prior::Union{ContinuousMultivariateDistribution, Void}
  setpdf::Union{Function, Void}
  setprior::Union{Function, Void}
  loglikelihood!::Union{Function, Void}
  logprior!::Union{Function, Void}
  logtarget!::Union{Function, Void}
  gradloglikelihood!::Union{Function, Void}
  gradlogprior!::Union{Function, Void}
  gradlogtarget!::Union{Function, Void}
  tensorloglikelihood!::Union{Function, Void}
  tensorlogprior!::Union{Function, Void}
  tensorlogtarget!::Union{Function, Void}
  dtensorloglikelihood!::Union{Function, Void}
  dtensorlogprior!::Union{Function, Void}
  dtensorlogtarget!::Union{Function, Void}
  uptogradlogtarget!::Union{Function, Void}
  uptotensorlogtarget!::Union{Function, Void}
  uptodtensorlogtarget!::Union{Function, Void}
  states::VariableStateVector

  function BasicContMuvParameter(
    key::Symbol,
    index::Integer,
    pdf::Union{ContinuousMultivariateDistribution, Void},
    prior::Union{ContinuousMultivariateDistribution, Void},
    setpdf::Union{Function, Void},
    setprior::Union{Function, Void},
    ll::Union{Function, Void},
    lp::Union{Function, Void},
    lt::Union{Function, Void},
    gll::Union{Function, Void},
    glp::Union{Function, Void},
    glt::Union{Function, Void},
    tll::Union{Function, Void},
    tlp::Union{Function, Void},
    tlt::Union{Function, Void},
    dtll::Union{Function, Void},
    dtlp::Union{Function, Void},
    dtlt::Union{Function, Void},
    uptoglt::Union{Function, Void},
    uptotlt::Union{Function, Void},
    uptodtlt::Union{Function, Void},
    states::VariableStateVector
  )
    args = (setpdf, setprior, ll, lp, lt, gll, glp, glt, tll, tlp, tlt, dtll, dtlp, dtlt, uptoglt, uptotlt, uptodtlt)
    fnames = fieldnames(BasicContMuvParameter)[5:21]

    # Check that all generic functions have correct signature
    for i in 1:17
      if isa(args[i], Function) &&
        isgeneric(args[i]) &&
        !method_exists(args[i], (BasicContMuvParameterState, VariableStateVector))
        error("$(fnames[i]) has wrong signature")
      end
    end

    new(
      key,
      index,
      pdf,
      prior,
      setpdf,
      setprior,
      ll,
      lp,
      lt,
      gll,
      glp,
      glt,
      tll,
      tlp,
      tlt,
      dtll,
      dtlp,
      dtlt,
      uptoglt,
      uptotlt,
      uptodtlt,
      states
    )
  end
end

function BasicContMuvParameter!(
  parameter::BasicContMuvParameter,
  setpdf::Union{Function, Void},
  setprior::Union{Function, Void},
  ll::Union{Function, Void},
  lp::Union{Function, Void},
  lt::Union{Function, Void},
  gll::Union{Function, Void},
  glp::Union{Function, Void},
  glt::Union{Function, Void},
  tll::Union{Function, Void},
  tlp::Union{Function, Void},
  tlt::Union{Function, Void},
  dtll::Union{Function, Void},
  dtlp::Union{Function, Void},
  dtlt::Union{Function, Void},
  uptoglt::Union{Function, Void},
  uptotlt::Union{Function, Void},
  uptodtlt::Union{Function, Void}
)
  args = (setpdf, setprior, ll, lp, lt, gll, glp, glt, tll, tlp, tlt, dtll, dtlp, dtlt, uptoglt, uptotlt, uptodtlt)

  # Define setpdf (i = 1) and setprior (i = 2)
  for (i, setter, distribution) in ((1, :setpdf, :pdf), (2, :setprior, :prior))
    setfield!(
      parameter,
      setter,
      if isa(args[i], Function)
        eval(codegen_setfield(parameter, distribution, args[i]))
      else
        nothing
      end
    )
  end

  # Define loglikelihood! (i = 3) and gradloglikelihood! (i = 6)
  # plfield stands for parameter likelihood-related field respectively
  for (i, plfield) in ((3, :loglikelihood!), (6, :gradloglikelihood!))
    setfield!(
      parameter,
      plfield,
      if isa(args[i], Function)
        eval(codegen_closure(parameter, args[i]))
      else
        nothing
      end
    )
  end

  # Define logprior! (i = 4) and gradlogprior! (i = 7)
  # ppfield and spfield stand for parameter prior-related field and state prior-related field repsectively
  for (i , ppfield, spfield, f) in ((4, :logprior!, :logprior, logpdf), (7, :gradlogprior!, :gradlogprior, gradlogpdf))
    setfield!(
    parameter,
      ppfield,
      if isa(args[i], Function)
        eval(codegen_closure(parameter, args[i]))
      else
        if (
            isa(parameter.prior, ContinuousMultivariateDistribution) &&
            method_exists(f, (typeof(parameter.prior), Vector{eltype(parameter.prior)}))
          ) ||
          isa(args[2], Function)
          eval(codegen_target_closure_via_distribution(parameter, :prior, f, spfield))
        else
          nothing
        end
      end
    )
  end

  # Define logtarget! (i = 5) and gradlogtarget! (i = 8)
  # ptfield, plfield and ppfield stand for parameter target, likelihood and prior-related field respectively
  # stfield, slfield and spfield stand for state target, likelihood and prior-related field respectively
  for (i , ptfield, plfield, ppfield, stfield, slfield, spfield, f) in (
    (5, :logtarget!, :loglikelihood!, :logprior!, :logtarget, :loglikelihood, :logprior, logpdf),
    (8, :gradlogtarget!, :gradloglikelihood!, :gradlogprior!, :gradlogtarget, :gradloglikelihood, :gradlogprior, gradlogpdf)
  )
    setfield!(
      parameter,
      ptfield,
      if isa(args[i], Function)
        eval(codegen_closure(parameter, args[i]))
      else
        if isa(args[i-2], Function) && isa(getfield(parameter, ppfield), Function)
          eval(codegen_sumtarget_closure(parameter, plfield, ppfield, stfield, slfield, spfield))
        elseif (
            isa(parameter.pdf, ContinuousMultivariateDistribution) &&
            method_exists(f, (typeof(parameter.pdf), Vector{eltype(parameter.pdf)}))
          ) ||
          isa(args[1], Function)
          eval(codegen_target_closure_via_distribution(parameter, :pdf, f, stfield))
        else
          nothing
        end
      end
    )
  end

  # Define tensorloglikelihood! (i = 9) and dtensorloglikelihood! (i = 12)
  # plfield stands for parameter likelihood-related field respectively
  for (i, plfield) in ((9, :tensorloglikelihood!), (12, :dtensorloglikelihood!))
    setfield!(
      parameter,
      plfield,
      if isa(args[i], Function)
        eval(codegen_closure(parameter, args[i]))
      else
        nothing
      end
    )
  end

  # Define tensorlogprior! (i = 10) and dtensorlogprior! (i = 13)
  # ppfield stands for parameter prior-related field respectively
  for (i, ppfield) in ((10, :tensorlogprior!), (13, :dtensorlogprior!))
    setfield!(
      parameter,
      ppfield,
      if isa(args[i], Function)
        eval(codegen_closure(parameter, args[i]))
      else
        nothing
      end
    )
  end

  # Define tensorlogtarget! (i = 11) and dtensorlogtarget! (i = 14)
  for (i , ptfield, plfield, ppfield, stfield, slfield, spfield) in (
    (
      11,
      :tensorlogtarget!, :tensorloglikelihood!, :tensorlogprior!,
      :tensorlogtarget, :tensorloglikelihood, :tensorlogprior
    ),
    (
      14,
      :dtensorlogtarget!, :dtensorloglikelihood!, :dtensorlogprior!,
      :dtensorlogtarget, :dtensorloglikelihood, :dtensorlogprior
    )
  )
    setfield!(
      parameter,
      ptfield,
      if isa(args[i], Function)
        eval(codegen_closure(parameter, args[i]))
      else
        if isa(args[i-2], Function) && isa(args[i-1], Function)
          eval(codegen_sumtarget_closure(parameter, plfield, ppfield, stfield, slfield, spfield))
        else
          nothing
        end
      end
    )
  end

  # Define uptogradlogtarget!
  setfield!(
    parameter,
    :uptogradlogtarget!,
    if isa(args[15], Function)
      eval(codegen_closure(parameter, args[15]))
    else
      if isa(parameter.logtarget!, Function) && isa(parameter.gradlogtarget!, Function)
        eval(codegen_uptotarget_closures(parameter, [:logtarget!, :gradlogtarget!]))
      else
        nothing
      end
    end
  )

  # Define uptotensorlogtarget!
  setfield!(
    parameter,
    :uptotensorlogtarget!,
    if isa(args[16], Function)
      eval(codegen_closure(parameter, args[16]))
    else
      if isa(parameter.logtarget!, Function) &&
        isa(parameter.gradlogtarget!, Function) &&
        isa(parameter.tensorlogtarget!, Function)
        eval(codegen_uptotarget_closures(parameter, [:logtarget!, :gradlogtarget!, :tensorlogtarget!]))
      else
        nothing
      end
    end
  )

  # Define uptodtensorlogtarget!
  setfield!(
    parameter,
    :uptodtensorlogtarget!,
    if isa(args[17], Function)
      eval(codegen_closure(parameter, args[17]))
    else
      if isa(parameter.logtarget!, Function) &&
        isa(parameter.gradlogtarget!, Function) &&
        isa(parameter.tensorlogtarget!, Function) &&
        isa(parameter.dtensorlogtarget!, Function)
        eval(codegen_uptotarget_closures(parameter, [:logtarget!, :gradlogtarget!, :tensorlogtarget!, :dtensorlogtarget!]))
      else
        nothing
      end
    end
  )
end

BasicContMuvParameter(key::Symbol, index::Integer=0; signature::Symbol=:high, args...) =
  BasicContMuvParameter(key, Val{signature}, index; args...)

function BasicContMuvParameter(
  key::Symbol,
  ::Type{Val{:low}},
  index::Integer=0;
  pdf::Union{ContinuousMultivariateDistribution, Void}=nothing,
  prior::Union{ContinuousMultivariateDistribution, Void}=nothing,
  setpdf::Union{Function, Void}=nothing,
  setprior::Union{Function, Void}=nothing,
  loglikelihood::Union{Function, Void}=nothing,
  logprior::Union{Function, Void}=nothing,
  logtarget::Union{Function, Void}=nothing,
  gradloglikelihood::Union{Function, Void}=nothing,
  gradlogprior::Union{Function, Void}=nothing,
  gradlogtarget::Union{Function, Void}=nothing,
  tensorloglikelihood::Union{Function, Void}=nothing,
  tensorlogprior::Union{Function, Void}=nothing,
  tensorlogtarget::Union{Function, Void}=nothing,
  dtensorloglikelihood::Union{Function, Void}=nothing,
  dtensorlogprior::Union{Function, Void}=nothing,
  dtensorlogtarget::Union{Function, Void}=nothing,
  uptogradlogtarget::Union{Function, Void}=nothing,
  uptotensorlogtarget::Union{Function, Void}=nothing,
  uptodtensorlogtarget::Union{Function, Void}=nothing,
  states::VariableStateVector=VariableState[]
)
  parameter = BasicContMuvParameter(key, index, pdf, prior, fill(nothing, 17)..., states)

  BasicContMuvParameter!(
    parameter,
    setpdf,
    setprior,
    loglikelihood,
    logprior,
    logtarget,
    gradloglikelihood,
    gradlogprior,
    gradlogtarget,
    tensorloglikelihood,
    tensorlogprior,
    tensorlogtarget,
    dtensorloglikelihood,
    dtensorlogprior,
    dtensorlogtarget,
    uptogradlogtarget,
    uptotensorlogtarget,
    uptodtensorlogtarget
  )

  parameter
end

function BasicContMuvParameter(
  key::Symbol,
  ::Type{Val{:high}},
  index::Integer=0;
  pdf::Union{ContinuousMultivariateDistribution, Void}=nothing,
  prior::Union{ContinuousMultivariateDistribution, Void}=nothing,
  setpdf::Union{Function, Void}=nothing,
  setprior::Union{Function, Void}=nothing,
  loglikelihood::Union{Function, Expr, Void}=nothing,
  logprior::Union{Function, Expr, Void}=nothing,
  logtarget::Union{Function, Expr, Void}=nothing,
  gradloglikelihood::Union{Function, Void}=nothing,
  gradlogprior::Union{Function, Void}=nothing,
  gradlogtarget::Union{Function, Void}=nothing,
  tensorloglikelihood::Union{Function, Void}=nothing,
  tensorlogprior::Union{Function, Void}=nothing,
  tensorlogtarget::Union{Function, Void}=nothing,
  dtensorloglikelihood::Union{Function, Void}=nothing,
  dtensorlogprior::Union{Function, Void}=nothing,
  dtensorlogtarget::Union{Function, Void}=nothing,
  uptogradlogtarget::Union{Function, Void}=nothing,
  uptotensorlogtarget::Union{Function, Void}=nothing,
  uptodtensorlogtarget::Union{Function, Void}=nothing,
  states::VariableStateVector=VariableState[],
  nkeys::Integer=0,
  vfarg::Bool=false,
  autodiff::Symbol=:none,
  order::Integer=1,
  chunksize::Integer=0,
  init::Vector=fill(Any[], 3)
)
  inargs = (
    setpdf,
    setprior,
    loglikelihood,
    logprior,
    logtarget,
    gradloglikelihood,
    gradlogprior,
    gradlogtarget,
    tensorloglikelihood,
    tensorlogprior,
    tensorlogtarget,
    dtensorloglikelihood,
    dtensorlogprior,
    dtensorlogtarget,
    uptogradlogtarget,
    uptotensorlogtarget,
    uptodtensorlogtarget
  )

  fnames = Array(Any, 17)
  fnames[1:2] = fill(Symbol[], 2)
  fnames[3:14] = [Symbol[f] for f in fieldnames(BasicContMuvParameterState)[2:13]]
  for i in 1:3
    fnames[14+i] = Symbol[fnames[j][1] for j in 5:3:(5+i*3)]
  end

  for i in 3:5
    if isa(inargs[i], Expr) && autodiff != :reverse
      error("The only case $(fnames[i][1]) can be an expression is when used in conjunction with reverse mode autodiff")
    end
  end

  if nkeys > 0
    if (autodiff == :forward || autodiff == :reverse) && vfarg
      error("In the case of autodiff, if nkeys is not 0, then vfarg must be false")
    end
  elseif nkeys < 0
    "nkeys must be non-negative, got $nkeys"
  end

  if !in(autodiff, (:none, :forward, :reverse))
    error("autodiff must be :nore or :forward or :reverse, got $autodiff")
  end

  if order < 0 || order > 2
    error("Derivative order must be 0, 1 or 2, got $order")
  elseif autodiff != :reverse && order == 0
    error("Zero order can be used only with reverse mode autodiff")
  end

  @assert chunksize >= 0 "chunksize must be non-negative, got $chunksize"

  initarg = Array(Any, 3)
  initlen = length(init)

  if initlen == 1 || initlen == 2
    if autodiff != :reverse
      @assert all(isempty, init) "init option is used only for reverse mode autodiff"
    end

    for i in 1:3
      initarg[i] = (inargs[i+2] != nothing) ? init : Any[]
    end
  elseif initlen == 3
    if autodiff != :reverse
      @assert all(isempty, init) "init option is used only for reverse mode autodiff"
    end

    initarg = init
  else
    error("init must be a vector of length 1, 2 or 3, got vector of length $initlen")
  end

  parameter = BasicContMuvParameter(key, index, pdf, prior, fill(nothing, 17)..., states)

  outargs = Union{Function, Void}[nothing for i in 1:17]

  for i in 1:17
    if isa(inargs[i], Function)
      outargs[i] = eval(
        codegen_lowlevel_variable_method(inargs[i], :BasicContMuvParameterState, true, fnames[i], nkeys, vfarg)
      )
    end
  end

  if autodiff == :forward
    fadclosure = Array(Union{Function, Void}, 3)
    for i in 3:5
      fadclosure[i-2] =
        if isa(inargs[i], Function)
          nkeys == 0 ? inargs[i] : eval(codegen_internal_autodiff_closure(parameter, inargs[i], nkeys))
        else
          nothing
        end
    end

    for i in 6:8
      if !isa(inargs[i], Function) && isa(inargs[i-3], Function)
        outargs[i] = eval(codegen_lowlevel_variable_method(
          eval(codegen_forward_autodiff_function(Val{:gradient}, fadclosure[i-5], chunksize)),
          :BasicContMuvParameterState,
          true,
          fnames[i],
          0
        ))
      end
    end

    if !isa(inargs[15], Function) && isa(inargs[5], Function)
      outargs[15] = eval(codegen_lowlevel_variable_method(
        eval(codegen_forward_autodiff_uptofunction(Val{:gradient}, fadclosure[3], chunksize)),
        :BasicContMuvParameterState,
        true,
        fnames[15],
        0
      ))
    end

    if order >= 2
      for i in 9:11
        if !isa(inargs[i], Function) && isa(inargs[i-6], Function)
          outargs[i] = eval(codegen_lowlevel_variable_method(
            eval(codegen_forward_autodiff_target(:hessian, fadclosure[i-8], chunksize)),
            :BasicContMuvParameterState,
            true,
            fnames[i],
            0
          ))
        end
      end

      if !isa(inargs[16], Function) && isa(inargs[5], Function)
        outargs[16] = eval(codegen_lowlevel_variable_method(
          eval(codegen_forward_autodiff_uptotarget(:hessian, fadclosure[3], chunksize)),
          :BasicContMuvParameterState,
          true,
          fnames[16],
          0
        ))
      end
    end
  elseif autodiff == :reverse
    local f::Function

    for i in 3:5
      if isa(inargs[i], Expr)
        if nkeys == 0
          f = eval(codegen_reverse_autodiff_function(inargs[i], :Vector, initarg[i-2][1], 0, false))
        else
          f = eval(codegen_reverse_autodiff_function(inargs[i], :Vector, initarg[i-2], 0, false))
          f = eval(codegen_internal_autodiff_closure(parameter, f, nkeys))
        end

        outargs[i] = eval(codegen_lowlevel_variable_method(f, :BasicContMuvParameterState, true, fnames[i], 0))
      end
    end

    for i in 6:8
      if !isa(inargs[i], Function)
        if isa(inargs[i-3], Function)
          if nkeys == 0
            f = ReverseDiffSource.rdiff(inargs[i-3], (initarg[i-5][1][2],), order=1, allorders=false)
          else
            f = ReverseDiffSource.rdiff(
              inargs[i-3], (initarg[i-5][1][2], initarg[i-5][2][2]), ignore=[initarg[i-5][2][1]], order=1, allorders=false
            )
            f = eval(codegen_internal_autodiff_closure(parameter, f, nkeys))
          end

          outargs[i] = eval(codegen_lowlevel_variable_method(f, :BasicContMuvParameterState, true, fnames[i], 0))
        elseif isa(inargs[i-3], Expr)
          if nkeys == 0
            f = eval(codegen_reverse_autodiff_function(inargs[i-3], :Vector, initarg[i-5][1], 1, false))
          else
            f = eval(codegen_reverse_autodiff_function(inargs[i-3], :Vector, initarg[i-5], 1, false))
            f = eval(codegen_internal_autodiff_closure(parameter, f, nkeys))
          end

          outargs[i] = eval(codegen_lowlevel_variable_method(f, :BasicContMuvParameterState, true, fnames[i], 0))
        end
      end
    end

    if !isa(inargs[15], Function)
      if isa(inargs[5], Function)
        if nkeys == 0
          f = ReverseDiffSource.rdiff(inargs[5], (initarg[3][1][2],), order=1, allorders=true)
        else
          f = ReverseDiffSource.rdiff(
            inargs[5], (initarg[3][1][2], initarg[3][2][2]), ignore=[initarg[3][2][1]], order=1, allorders=true
          )
          f = eval(codegen_internal_autodiff_closure(parameter, f, nkeys))
        end

        outargs[15] = eval(codegen_lowlevel_variable_method(f, :BasicContMuvParameterState, true, fnames[15], 0))
      elseif isa(inargs[5], Expr)
        if nkeys == 0
          f = eval(codegen_reverse_autodiff_function(inargs[5], :Vector, initarg[3][1], 1, true))
        else
          f = eval(codegen_reverse_autodiff_function(inargs[5], :Vector, initarg[3], 1, true))
          f = eval(codegen_internal_autodiff_closure(parameter, f, nkeys))
        end

        outargs[15] = eval(codegen_lowlevel_variable_method(f, :BasicContMuvParameterState, true, fnames[15], 0))
      end
    end

    if order >= 2
      for i in 9:11
        if !isa(inargs[i], Function)
          if isa(inargs[i-6], Function) || isa(inargs[i-6], Expr)
            if nkeys == 0
              f = eval(codegen_reverse_autodiff_target(:hessian, inargs[i-6], :Vector, initarg[i-8][1]))
            else
              f = eval(codegen_reverse_autodiff_target(:hessian, inargs[i-6], :Vector, initarg[i-8]))
              f = eval(codegen_internal_autodiff_closure(parameter, f, nkeys))
            end

            outargs[i] = eval(codegen_lowlevel_variable_method(f, :BasicContMuvParameterState, true, fnames[i], 0))
          end
        end
      end

      if !isa(inargs[16], Function)
        if isa(inargs[5], Function) || isa(inargs[5], Expr)
          if nkeys == 0
            f = eval(codegen_reverse_autodiff_uptotarget(:hessian, inargs[5], :Vector, initarg[3][1]))
          else
            f = eval(codegen_reverse_autodiff_uptotarget(:hessian, inargs[5], :Vector, initarg[3]))
            f = eval(codegen_internal_autodiff_closure(parameter, f, nkeys))
          end

          outargs[16] = eval(codegen_lowlevel_variable_method(f, :BasicContMuvParameterState, true, fnames[16], 0))
        end
      end
    end
  end

  BasicContMuvParameter!(parameter, outargs...)

  parameter
end

function codegen_internal_autodiff_closure(parameter::BasicContMuvParameter, f::Function, nkeys::Integer)
  fstatesarg = [Expr(:ref, :Any, [:($(parameter).states[$i].value) for i in 1:nkeys]...)]

  @gensym internal_forward_autodiff_closure

  quote
    function $internal_forward_autodiff_closure(_x::Vector)
      $(f)(_x, $(fstatesarg...))
    end
  end
end

value_support(::Type{BasicContMuvParameter}) = Continuous
value_support(::BasicContMuvParameter) = Continuous

variate_form(::Type{BasicContMuvParameter}) = Multivariate
variate_form(::BasicContMuvParameter) = Multivariate

default_state_type(::BasicContMuvParameter) = BasicContMuvParameterState

default_state{N<:Real}(variable::BasicContMuvParameter, value::Vector{N}, outopts::Dict) =
  BasicContMuvParameterState(
    value,
    [getfield(variable, fieldnames(BasicContMuvParameter)[i]) == nothing ? false : true for i in 10:18],
    (haskey(outopts, :diagnostics) && in(:accept, outopts[:diagnostics])) ? [:accept] : Symbol[]
  )

Base.show(io::IO, ::Type{BasicContMuvParameter}) = print(io, "BasicContMuvParameter")
Base.writemime(io::IO, ::MIME"text/plain", t::Type{BasicContMuvParameter}) = show(io, t)


text:
module Gitea

# package code goes here
import Requests
import HttpCommon
import JSON

export Client

DEBUG = false

include("utils.jl")
include("errors.jl")
include("client.jl")

# Types
include("types.jl")

# Functions
include("user.jl")
include("repo.jl")
include("issue.jl")
include("pull.jl")
include("org.jl")
include("release.jl")
include("status.jl")
include("hook.jl")
include("miscellaneous.jl")

# include("admin.jl")

end # module