text:
using DataFrames
import MultivariateStats
using NearestNeighbors
using StatsBase: countmap, denserank
using Statistics

using Base.Threads

import CSV
import Distances

# Parse parameters

function default_param_value(param::Symbol, min_molecules_per_cell::Union{Int, Nothing};
                             n_molecules::Union{Int, Nothing}=nothing, n_genes::Union{Int, Nothing}=nothing)
    if min_molecules_per_cell === nothing
        error("Either `$param` or `min_molecules_per_cell` must be provided")
    end

    min_molecules_per_cell = max(min_molecules_per_cell, 3)

    if param == :min_molecules_per_segment
        return max(round(Int, min_molecules_per_cell / 4), 2)
    end

    if param == :confidence_nn_id
        return max(div(min_molecules_per_cell, 2) + 1, 5)
    end

    if param == :composition_neighborhood
        if n_genes === nothing
            return max(min_molecules_per_cell, 3)
        end

        return max(div(n_genes, 10), min_molecules_per_cell, 3)
    end

    if param == :n_gene_pcs
        (n_genes !== nothing) || error("Either `$param` or `n_genes` must be provided")
        return min(max(div(n_genes, 3), 30), 100, n_genes)
    end

    if param == :n_cells_init
        (n_molecules !== nothing) || error("Either `$param` or `n_molecules` must be provided")
        return div(n_molecules, min_molecules_per_cell) * 2
    end
end

default_if_not_provided(value::Union{<:Real, Nothing}, param_name::Symbol, args...; kwargs...) =
    (value === nothing) ? default_param_value(param_name, args...; kwargs...) : value

parse_scale_std(scale_std::Float64, ::Real) = scale_std
parse_scale_std(scale_std::Nothing, scale::Real) = 0.25 * scale
function parse_scale_std(scale_std::String, scale::Real)
    scale_std = strip(scale_std)
    if scale_std[end] == '%'
        return parse(Float64, scale_std[1:end-1]) * scale / 100.0
    end

    return parse(Float64, scale_std)
end

# Initialize cell positions

function position_data_by_assignment(pos_data::Matrix{Float64}, assignment::Vector{Int})
    filt_ids = findall(assignment .> 0)
    return [pos_data[:, ids] for ids in split(filt_ids, assignment[filt_ids])]
end

function covs_from_assignment(pos_data::Matrix{Float64}, assignment::Vector{Int}; min_size::Int=1)::Array{<:CovMat, 1}
    pos_data_by_assignment = position_data_by_assignment(pos_data, assignment)
    CM, MV = (size(pos_data, 1) == 2) ? (CovMat{2}, MeanVec{2}) : (CovMat{3}, MeanVec{3})

    covs = [estimate_sample_cov!(zeros(CM), pd) for pd in pos_data_by_assignment];
    mean_covs = MV(vec(median(hcat(Vector.(eigvals.(covs[size.(pos_data_by_assignment, 2) .>= min_size]))...), dims=2)))

    for i in findall(size.(pos_data_by_assignment, 2) .<= min_size)
        covs[i] = CM(diagm(0 => deepcopy(mean_covs)))
    end

    return covs
end

cell_centers_uniformly(spatial_df::DataFrame, args...; kwargs...) =
    cell_centers_uniformly(position_data(spatial_df), args..., (:confidence in propertynames(spatial_df)) ? spatial_df.confidence : nothing; kwargs...)

function cell_centers_uniformly(pos_data::Matrix{Float64}, n_clusters::Int,
        confidences::Union{Vector{Float64}, Nothing}=nothing; scale::Union{Float64, Nothing})
    n_clusters = min(n_clusters, size(pos_data, 2))

    cluster_centers = pos_data[:, select_ids_uniformly(pos_data', confidences; n=n_clusters, confidence_threshold=0.25)]
    cluster_labels = vcat(knn(KDTree(cluster_centers), pos_data, 1)[1]...)

    covs = nothing
    if scale === nothing
        covs = covs_from_assignment(pos_data, cluster_labels)
    elseif size(pos_data, 1) == 2
        scale = scale ^ 2
        covs = [CovMat{2}(diagm(0 => (ones(2) .* scale))) for i in 1:size(cluster_centers, 2)]
    elseif size(pos_data, 1) == 3
        scale = scale ^ 2
        covs = [CovMat{3}(diagm(0 => (ones(3) .* scale))) for i in 1:size(cluster_centers, 2)]
    else
        error("Unexpected number of dimensions: $(size(pos_data, 1))")
    end

    return InitialParams(copy(cluster_centers'), covs, cluster_labels)
end

# Build MRF

function convert_edge_list_to_adj_list(edge_list::Matrix{Int}, edge_weights::Union{Vector{Float64}, Nothing}=nothing; n_verts::Int=maximum(edge_list))
    res_ids = [vcat(v...) for v in zip(split(edge_list[2,:], edge_list[1,:], max_factor=n_verts),
                                       split(edge_list[1,:], edge_list[2,:], max_factor=n_verts))];

    if edge_weights === nothing
        return res_ids
    end

    res_weights = [vcat(v...) for v in zip(split(edge_weights, edge_list[1,:], max_factor=n_verts),
                                           split(edge_weights, edge_list[2,:], max_factor=n_verts))];

    return res_ids, res_weights
end

function estimate_local_composition_similarities(df_spatial::DataFrame, edge_list::Matrix{Int}; composition_neighborhood::Int, n_gene_pcs::Int, min_similarity::Float64=0.01)
    neighb_cm = neighborhood_count_matrix(df_spatial, composition_neighborhood);
    neighb_cm = MultivariateStats.transform(MultivariateStats.fit(MultivariateStats.PCA, neighb_cm, maxoutdim=n_gene_pcs), neighb_cm);

    return max.(map(i -> cor(neighb_cm[:, edge_list[1,i]], neighb_cm[:, edge_list[2,i]]), 1:size(edge_list, 2)), min_similarity);
end

function build_molecule_graph(df_spatial::DataFrame; min_edge_quant::Float64=0.3, use_local_gene_similarities::Bool=false, n_gene_pcs::Int=0, composition_neighborhood::Int=0, kwargs...)
    if use_local_gene_similarities && ((composition_neighborhood == 0) || (n_gene_pcs == 0))
        @warn "composition_neighborhood=$composition_neighborhood and n_gene_pcs=$n_gene_pcs, while use_local_gene_similarities=true. Force use_local_gene_similarities=false."
        use_local_gene_similarities = false
    end

    if use_local_gene_similarities && any(ismissing.(composition_data(df_spatial)))
        @warn "use_local_gene_similarities is not supported with missing gene values. Setting to false."
        use_local_gene_similarities = false
    end

    edge_list, adjacent_dists = adjacency_list(df_spatial; kwargs...);

    min_edge_length = quantile(adjacent_dists, min_edge_quant);

    adjacent_weights = min_edge_length ./ max.(adjacent_dists, min_edge_length)
    if use_local_gene_similarities
        adjacent_weights .*= estimate_local_composition_similarities(df_spatial, edge_list; composition_neighborhood=composition_neighborhood, n_gene_pcs=n_gene_pcs)
    end

    adjacent_points, adjacent_weights = convert_edge_list_to_adj_list(edge_list, adjacent_weights; n_verts=size(df_spatial, 1));
    return adjacent_points, adjacent_weights, adjacent_dists
end

function build_molecule_graph_normalized(df_spatial::DataFrame, vertex_weights::Union{Vector{Float64}, Symbol}; kwargs...)
    if typeof(vertex_weights) <: Symbol
        vertex_weights = df_spatial[!, vertex_weights]
    end

    adjacent_points, adjacent_weights = build_molecule_graph(df_spatial; kwargs...)
    for i in 1:length(adjacent_weights)
        cur_points = adjacent_points[i]
        cur_weights = adjacent_weights[i]
        for j in 1:length(cur_weights)
            cur_weights[j] *= vertex_weights[cur_points[j]]
        end
    end

    return adjacent_points, adjacent_weights
end

function match_gene_names(gene_masks::Vector{String}, gene_names::Vector{String})
    matches = Set{String}()
    missing_genes = String[]
    for gm in gene_masks
        g_ids = findall(match.(Regex(gm), gene_names) .!== nothing)
        if length(g_ids) == 0
            push!(missing_genes, gm)
        else
            union!(matches, gene_names[g_ids])
        end
    end

    (length(missing_genes) == 0) || @warn "Genes $(join(missing_genes, ',')) are missing from the data"
    return matches
end

function load_df(data_path::String; min_molecules_per_gene::Int=0, exclude_genes::Vector{String}=String[], kwargs...)
    df_spatial = read_spatial_df(data_path; kwargs...)

    gene_counts = StatsBase.countmap(df_spatial[!, :gene]);
    large_genes = Set{String}(collect(keys(gene_counts))[collect(values(gene_counts)) .>= min_molecules_per_gene]);
    df_spatial = df_spatial[in.(df_spatial.gene, Ref(large_genes)),:];

    if length(exclude_genes) > 0
        exclude_genes = match_gene_names(exclude_genes, unique(df_spatial.gene))
        df_spatial = df_spatial[.!in.(df_spatial.gene, Ref(exclude_genes)),:];
        @info "Excluding genes: " * join(sort(collect(exclude_genes)), ", ")
    end

    df_spatial[!, :x] = Array{Float64, 1}(df_spatial[!, :x])
    df_spatial[!, :y] = Array{Float64, 1}(df_spatial[!, :y])
    df_spatial[!, :gene], gene_names = encode_genes(df_spatial[!, :gene]);
    return df_spatial, gene_names;
end

# Initialize BmmData

"""
    main function for initialization of bm_data
"""
function initialize_bmm_data(df_spatial::DataFrame; min_molecules_per_cell::Int, scale::T where T<: Real, scale_std::Union{<:Real, String, Nothing}=nothing,
        n_cells_init::Union{Int, Nothing}=nothing, confidence_nn_id::Union{Int, Nothing}=nothing, composition_neighborhood::Union{Int, Nothing}=nothing,
        adjacent_points::Union{Vector{Vector{Int64}}, Nothing}=nothing, adjacent_weights::Union{Vector{Vector{Float64}}, Nothing}=nothing,
        use_local_gene_similarities::Bool=true, adjacency_type::Symbol=:triangulation, n_gene_pcs::Union{Int, Nothing}=nothing, prior_seg_confidence::Float64=0.5,
        verbose::Bool=true, kwargs...)::BmmData
    df_spatial = deepcopy(df_spatial)

    ## Parse parameters
    confidence_nn_id = default_if_not_provided(confidence_nn_id, :confidence_nn_id, min_molecules_per_cell)
    n_cells_init = default_if_not_provided(n_cells_init, :n_cells_init, min_molecules_per_cell, n_molecules=size(df_spatial, 1))
    scale_std = parse_scale_std(scale_std, scale)

    ## Estimate confidence
    if confidence_nn_id > 0
        verbose && @info "Estimate confidence per molecule"
        prior_segmentation = (:prior_segmentation in propertynames(df_spatial)) ? df_spatial.prior_segmentation : nothing
        append_confidence!(df_spatial, prior_segmentation; nn_id=confidence_nn_id, prior_confidence=prior_seg_confidence)
        verbose && @info "Done"
    end

    ## Initialize BmmData array
    verbose && @info "Initializing algorithm. Scale: $scale, scale std: $scale_std, initial #components: $n_cells_init, #molecules: $(size(df_spatial, 1))."
    size_prior = (:z in propertynames(df_spatial)) ?
        ShapePrior{3}(Float64[scale, scale, scale], Float64[scale_std, scale_std, scale_std], min_molecules_per_cell) :
        ShapePrior{2}(Float64[scale, scale], Float64[scale_std, scale_std], min_molecules_per_cell)

    init_params = cell_centers_uniformly(df_spatial, n_cells_init; scale=scale)
    if adjacent_points === nothing
        composition_neighborhood = default_if_not_provided(composition_neighborhood, :composition_neighborhood,
            min_molecules_per_cell, n_genes=maximum(df_spatial.gene))
        n_gene_pcs = default_if_not_provided(n_gene_pcs, :n_gene_pcs, min_molecules_per_cell, n_genes=maximum(df_spatial.gene))

        adjacent_points, adjacent_weights = build_molecule_graph(df_spatial; use_local_gene_similarities=use_local_gene_similarities,
            n_gene_pcs=n_gene_pcs, composition_neighborhood=composition_neighborhood, adjacency_type=adjacency_type)[1:2]
    end

    components, sampler, assignment = initial_distributions(df_spatial, init_params; size_prior=size_prior)

    return BmmData(components, df_spatial, adjacent_points, adjacent_weights, assignment, sampler;
        prior_seg_confidence=prior_seg_confidence, kwargs...)
end

function initial_distributions(df_spatial::DataFrame, initial_params::InitialParams; size_prior::ShapePrior{N},
                               gene_smooth::Real=1.0, gene_num::Int=maximum(skipmissing(composition_data(df_spatial)))) where N
    position_distrubutions = [MvNormalF(initial_params.centers[i,:], initial_params.covs[i]) for i in 1:size(initial_params.centers, 1)]
    gene_distributions = [CategoricalSmoothed(ones(Float64, gene_num), smooth=Float64(gene_smooth)) for i in 1:length(position_distrubutions)]

    components = [Component(pd, gd, shape_prior=deepcopy(size_prior)) for (pd, gd) in zip(position_distrubutions, gene_distributions)]

    gene_sampler = CategoricalSmoothed(ones(Float64, gene_num))
    sampler = Component(MvNormalF(zeros(N)), gene_sampler, shape_prior=deepcopy(size_prior)) # position_params are never used

    return components, sampler, initial_params.assignment
end

# Utils

function staining_value_per_transcript(df_spatial::DataFrame, staining::MT where MT <: AbstractMatrix{T}; quiet::Bool=false)::Vector{T} where T <: Real
    x_vals = round.(Int, df_spatial.x)
    y_vals = round.(Int, df_spatial.y)
    if !quiet && ((maximum(x_vals) > size(staining, 2)) || (maximum(y_vals) > size(staining, 1)))
        @warn "Maximum transcript coordinates are $((maximum(y_vals), maximum(x_vals))), which is larger than the DAPI size: $(size(staining)). Filling it with 0."
    end

    if !quiet && ((minimum(x_vals) < 1) || (minimum(y_vals) < 1))
        @warn "Minimum transcript coordinates are < 1: $((minimum(y_vals), minimum(x_vals))). Filling it with 0."
    end

    if !quiet && ((maximum(x_vals) < 0.5 * size(staining, 2)) || (maximum(y_vals) < 0.5 * size(staining, 1)))
        @warn "Maximum transcript coordinates are $((maximum(y_vals), maximum(x_vals))), which is much smaller than the DAPI size: $(size(staining)). May be result of an error."
    end

    inds = findall((x_vals .> 0) .& (x_vals .<= size(staining, 2)) .& (y_vals .> 0) .& (y_vals .<= size(staining, 1)))
    staining_vals = zeros(T, length(x_vals))
    for i in inds
        staining_vals[i] = staining[y_vals[i], x_vals[i]];
    end

    return staining_vals
end

function encode_genes(gene_list)
    gene_names = sort(unique(gene_list));
    gene_ids = Dict(zip(gene_names, 1:length(gene_names)))
    return [gene_ids[g] for g in gene_list], gene_names
end

# This function is a determenistic analogue of sampling. It picks points in a manner that preserves the distributions across x and y.
function select_ids_uniformly(vals::Union{Vector{<:Real}, <:AbstractMatrix{<:Real}}, confidence::Union{Vector{Float64}, Nothing}=nothing; n::Int, confidence_threshold::Float64=0.95)::Vector{Int}
    if n <= 1
        error("n must be > 1")
    end

    high_conf_ids = (confidence===nothing) ? (1:size(vals, 1)) : findall(confidence .>= confidence_threshold)
    if length(high_conf_ids) < n
        @warn "n=$n, which is > length(high_conf_ids) ($(length(high_conf_ids)))"
        return high_conf_ids
    end

    vals = sum(vals, dims=2)[high_conf_ids]
    return high_conf_ids[sortperm(vals)[unique(round.(Int, range(1, length(high_conf_ids), length=n)))]];
end

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
module Savefig
using Plots, PGFPlots, Colors, TikzPictures
using Plots: Subplot
color_palette = [
    RGB(0.1399999, 0.1399999, 0.4),
    RGB(1.0, 0.7075, 0.35),
    RGB(0.414999, 1.0, 1.0),
    RGB(0.6, 0.21, 0.534999),
    RGB(0,0.6,0),
]

default(grid=false, color_palette=color_palette)

function replace_a_bunch_of_stuff(filename, replacements)
    temppath, temp = mktemp()
    open(filename) do f
        text = read(f, String)
        for p in replacements
            text = replace(text, p)
        end
        print(temp, text)
    end
    close(temp)
    mv(temppath, filename, force=true)
end

tt(x) = permutedims(x,(2,1))
function set_positions(filename, layout)
    temppath, temp = mktemp()
    fignum = 0
    grid = size(layout)
    indices = tt(CartesianIndices(grid))
    named = Int[]
    for line in eachline(filename)
        if occursin("\\begin{axis}[", line)
            fignum += 1
            tt(layout.grid)[fignum] isa Subplot || (fignum += 1)
            println(temp, line)
            println(temp, "name = fig$(fignum),")
            push!(named, fignum)
            fignum == 1 && continue
            if (indices[fignum][1] == 1 && (fignum-1) ∈ named) || (fignum-grid[1]) ∉ named # The leftmost, anchor north #
                println(temp, "at=(fig$(fignum-1).right of south east),")
                println(temp, "anchor=left of south west,")
                println(temp, "xshift=3mm,")
            else # otherwise anchor west
                println(temp, "at=(fig$(fignum-grid[1]).below south west),")
                println(temp, "anchor=above north west,")
                println(temp, "yshift=-3mm,")
            end
            continue
        end
        occursin(Char(8), line) && continue
        println(temp, line)
    end

    close(temp)
    mv(temppath, filename, force=true)

end

replacements_pgfplots = [
r"axis background/.style={fill={rgb,1:red,[\d\.]+;green,[\d\.]+;blue,[\d\.]+}}\s*?," => "",
"rgb,1:red,0.00000000;green,0.00000000;blue,0.00000000" => "black",
"rgb,1:red,1.00000000;green,1.00000000;blue,1.00000000" => "white",
r"height = \{[\.\d]+\w*?\}" => "height = {\\figureheight}",
r"width = \{[\.\d]+\w*?\}" => "width = {\\figurewidth}",
r",[\w\s]+? style = \{font = \{\\fontsize\{[\w\s\d\.]*?\}\{[\w\s\d\.]*?\}\\selectfont\}, color = \{rgb,1:red,0.00000000;green,0.00000000;blue,0.00000000\}, draw opacity = 1.0, rotate = 0.0\}" => "",
r",[\w\s]+? style = {color = \{rgb,1:red,0.00000000;green,0.00000000;blue,0.00000000\},\n*?draw opacity = 1.0,\n*?line width = 1,\n*?solid(,fill = \{rgb,1:red,1.00000000;green,1.00000000;blue,1.00000000\},font = \{\\fontsize\{[\w\s\d\.]*?\}\{[\w\s\d\.]*?\}\\selectfont})?\}" => "",
r"xtick\s?=\s?\{.*?\}," => "",
r"ytick\s?=\s?\{.*?\}," => "",
r"xshift\s?=\s?\d+?\.\d+?[^\\]\w*?," => ",",
r"yshift\s?=\s?\d+?\.\d+?[^\\]\w*?," => ",",
r"xticklabels\s?=\s?\{.*?\}," => "",
r"yticklabels\s?=\s?\{.*?\}," => "",
]

function savefig_pgfplots(filename, axis="", fig = current())

    layout = size(fig.layout.grid)
    pgffig = fig.o
    subfig = Iterators.product(0:layout[2]-1, 0:layout[1]-1)
    for (f,sfig) in zip(pgffig, subfig)
        # stylevec = split(f.style, ",")
        # push!(stylevec, "xshift=$(1.05*sfig[1])\\figurewidth, yshift=$(1.05*sfig[2])\\figureheight")
        # f.style = join(stylevec, ',')
        f.style = "xshift=$(1.02*sfig[1])\\figurewidth, yshift=$(1.05*sfig[2])\\figureheight,"*axis
        f.height = "\\figureheight"
        f.width = "\\figurewidth"
    end
    PGFPlots.save(filename, pgffig, include_preamble=false)
    # replace_a_bunch_of_stuff(filename, replacements_pgfplots)
end

# function savefig_pgfplots(filename, axis="")
#     fig = current()
#     layout = size(fig.layout.grid)
#     pgffig = fig.o
#     subfig = Iterators.product(0:layout[2]-1, 0:layout[1]-1)
#     for (f,sfig) in zip(pgffig, subfig)
#         f.style = axis
#         f.height = "\\figureheight"
#         f.width = "\\figurewidth"
#     end
#     tikzfig = PGFPlots.plot(pgffig)
#     TikzPictures.save(TEX(filename, include_preamble=false), tikzfig) # Goingthrough TikzPictures manually is required. If saved using PGFPlots.save the subplots layout is completely screwed up
#     set_positions(filename, fig.layout)
# end

using PyCall
const tikzplotlib = PyNULL()

function __init__()
    @eval global tikzplotlib = pyimport("tikzplotlib")
end

replacements_pyplot = [
r"\\definecolor\{color\d+\}\{rgb\}\{.*?\}" => Char(8),
r"xtick\s?=\s?\{.*?\}," => Char(8),
r"ytick\s?=\s?\{.*?\}," => Char(8),
r"xticklabels\s?=\s?\{.*?\}," => Char(8),
r"yticklabels\s?=\s?\{.*?\}," => Char(8),
"axis line style={none}," => Char(8),
]
"""
`savefig_pyplot(filename, fig = current().o; extra::Vector{String})`
"""
function savefig_pyplot(filename, fig = current(); kwargs...)
    tikzplotlib.save(filename,fig.o; figureheight = "\\figureheight", figurewidth = "\\figurewidth", kwargs...)
    replace_a_bunch_of_stuff(filename, replacements_pyplot)
    set_positions(filename, fig.layout)
end

end # module

# cd(@__DIR__)
# plot(randn(100), label="label", title="title", xlabel="xlabel", ylabel="ylabel")
# savefig_pgfplots("test.tex")

ext: jl


text:
transduce_shfl(xf::Transducer, op, init, xs; kwargs...) =
    transduce_shfl(xf'(op), init, xs; kwargs...)

function transduce_shfl(op, init, xs;)
    xf0, coll = extract_transducer(xs)
    # TODO: more systematic approach to this (and also support product)
    if coll isa Iterators.Zip
        arrays = coll.is
        xf = xf0
    elseif coll isa Iterators.Pairs
        arrays = (keys(coll), values(coll))
        xf = xf0
    else
        arrays = (coll,)
        xf = opcompose(Map(first), xf0)
    end
    rf = _reducingfunction(xf, op; init = init)
    acc = transduce_shfl_impl(rf, init, arrays...)
    rf_dev = cudaconvert(rf)
    if rf_dev === rf
        result = complete(rf, acc)
    else
        result = complete_on_device(rf_dev, acc)
    end
    if unreduced(result) isa DefaultInitOf
        throw(EmptyResultError(rf))
    end
    return result
end

macro _inbounds(ex)
    ex = :($Base.@inbounds $ex)
    esc(ex)
end

function transduce_shfl_impl(rf::F, init, arrays...) where {F}
    ys, = (dest, buf) = transduce_shfl!(nothing, rf, init, arrays...)
    if buf === nothing
        # The accumulator is a singleton. Once we are finished with the
        # side-effects of the basecase, transduce is done:
        return ys
    end
    # @info "ys, = transduce_shfl!(nothing, rf, ...)" Text(summary(ys))
    # @info "ys, = transduce_shfl!(nothing, rf, ...)" collect(ys)
    length(ys) == 1 && return @allowscalar ys[1]
    rf2 = AlwaysCombine(rf)
    while true
        ys, = transduce_shfl!(buf, rf2, init, ys)
        # @info "ys, = transduce_shfl!(buf, rf2, ...)" Text(summary(ys))
        # @info "ys, = transduce_shfl!(buf, rf2, ...)" collect(ys)
        length(ys) == 1 && return @allowscalar ys[1]
        dest, buf = buf, dest
        # reusing buffer; is it useful?
    end
end

@inline function transduce_shfl!(buf, rf, init, arrays...)
    idx = eachindex(arrays...)
    n = Int(length(idx))  # e.g., `length(UInt64(0):UInt64(1))` is not an `Int`

    dev = device()
    wsize = warpsize(dev)
    WARP_SIZE = Val(wsize)

    acctype = if buf === nothing
        # global _ARGS = (rf, zip(arrays...), init)
        # @show fake_transduce(rf, zip(arrays...), init)
        fake_args = (cudaconvert(rf), zip(map(cudaconvert, arrays)...), cudaconvert(init), Val(true))
        fake_args_tt = Tuple{map(Typeof, fake_args)...}
        # global FAKE_ARGS = fake_args
        # global FAKE_ARGS_TT = fake_args_tt
        return_type(fake_transduce, fake_args_tt)
        # Note: the result of `return_type` is not observable by the
        # caller of the API `transduce_impl`
    else
        eltype(buf)
    end
    buf0 = if Base.issingletontype(acctype)
        nothing
    elseif buf === nothing
        # TODO: find a way to compute type for `cufunction` without
        # creating a dummy object.
        allocate_buffer(acctype, 0)
    else
        buf
    end
    args = (buf0, Val{acctype}(), WARP_SIZE, rf, init, 0, idx, arrays...)
    # global _KARGS = args
    kernel_tt = Tuple{map(x -> Typeof(cudaconvert(x)), args)...}
    # global KERNEL_TT = kernel_tt
    kernel = cufunction(transduce_shfl_kernel!, kernel_tt)
    effelsize = if isbitstype(acctype)
        sizeof(acctype)
    else
        sizeof(UnionArrays.buffereltypefor(acctype)) + sizeof(UInt8)
    end
    # @show acctype UnionArrays.buffereltypefor(acctype) effelsize
    kernel_config = launch_configuration(kernel.fun)
    # @show kernel_config
    threads = let wanted_threads = nextwarp(dev, n)
        given_threads = kernel_config.threads
        # @show wanted_threads
        if wanted_threads > given_threads
            prevwarp(dev, given_threads)
        else
            wanted_threads
        end
    end

    @assert threads <= wsize * wsize  # = 32 * 32 = 1024

    nwarps_per_block, _nwarps_rem = divrem(threads, wsize)
    # @show threads nwarps_per_block _nwarps_rem
    @assert _nwarps_rem % wsize == 0
    basesize = nextwarp(dev, cld(n, kernel_config.blocks * nwarps_per_block))
    blocks = cld(n, basesize * nwarps_per_block)
    @assert blocks <= kernel_config.blocks

    # @show threads, blocks, basesize, acctype
    if Base.issingletontype(acctype)
        @cuda(
            threads = threads,
            blocks = blocks,
            transduce_shfl_kernel!(
                nothing,
                Val{acctype}(),
                WARP_SIZE,
                rf,
                init,
                basesize,
                idx,
                arrays...,
            )
        )
        return acctype.instance, nothing
    end

    if buf === nothing
        dest_buf = allocate_buffer(acctype, blocks + cld(blocks, threads))
        dest = view(dest_buf, 1:blocks)
        buf = view(dest_buf, blocks+1:length(dest_buf))
    else
        dest = view(buf, 1:blocks)
    end

    @cuda(
        threads = threads,
        blocks = blocks,
        transduce_shfl_kernel!(
            dest,
            Val{acctype}(),
            WARP_SIZE,
            rf,
            init,
            basesize,
            idx,
            arrays...,
        )
    )

    return dest, buf
end

@inline function transduce_shfl_kernel!(
    dest::Union{AbstractArray,Nothing},
    ::Val{T},
    ::Val{WARP_SIZE},
    rf::F,
    init,
    basesize,  # length of the main loop per warp
    idx,
    arrays...,
) where {T,WARP_SIZE,F}

    @inline function _shfl_down(x, delta)
        if dest === nothing
            return x
        end
        uv = unionvalue(T, x)
        return interpret(shfl_down_sync(typemax(UInt32), uv, delta, WARP_SIZE))
    end

    nwarps_per_block, _nwarps_rem = divrem(blockDim().x, WARP_SIZE)
    @assert _nwarps_rem == 0
    warpIdx0, warp_offset = divrem(threadIdx().x - 1, WARP_SIZE)
    warpIdx = warpIdx0 + 1
    warp_leader = 1 + warpIdx0 * WARP_SIZE  # first thread of this warp
    @assert warp_leader + warp_offset == threadIdx().x

    main_offset = warpIdx0 + (blockIdx().x - 1) * nwarps_per_block
    main_bound0 = basesize * (main_offset + 1)  # bound for this warp in `eachindex(idx)`
    need_remainder = main_bound0 > lastindex(idx)  # `idx` too short for the last block
    main_bound = min(main_bound0, lastindex(idx)) - WARP_SIZE
    need_remainder && @assert blockIdx().x == gridDim().x

    # if warp_offset == 0
    #     @cuprintf("%03ld: warpIdx0 = %ld main_offset = %d\n", threadIdx().x, Int(warpIdx0), Int(main_offset))
    # end

    # Main O(N) loop:
    acc = start(rf, init)
    warp_leader_offset = basesize * main_offset  # offset for this warp
    while warp_leader_offset <= main_bound
        i = warp_leader_offset + warp_offset + 1
        # @cuprintf("%03ld: i = %d\n", threadIdx().x, Int(i))
        acc = @manual_union_split(
            acc isa ithtype(T, Val(1)),
            acc isa ithtype(T, Val(2)),
            acc isa ithtype(T, Val(3)),
            acc isa ithtype(T, Val(4)),
            acc isa ithtype(T, Val(5)),
            acc isa ithtype(T, Val(6)),
        ) do
            next(rf, acc, @_inbounds getvalues(idx[i], arrays...))
        end

        # Warp-wide merge:
        delta = 1
        while delta < WARP_SIZE
            acc = _combine(rf, acc, _shfl_down(acc, delta))
            delta <<= 1
        end
        if warp_offset != 0
            acc = start(rf, init)
        end
        warp_leader_offset += WARP_SIZE
    end

    # Remainder of the main loop:
    if need_remainder
        let i = warp_leader_offset + warp_offset + 1
            # @cuprintf("%03ld: (rem) i = %d\n", threadIdx().x, Int(i))
            if i <= lastindex(idx)
                acc = @manual_union_split(
                    acc isa ithtype(T, Val(1)),
                    acc isa ithtype(T, Val(2)),
                    acc isa ithtype(T, Val(3)),
                    acc isa ithtype(T, Val(4)),
                    acc isa ithtype(T, Val(5)),
                    acc isa ithtype(T, Val(6)),
                ) do
                    next(rf, acc, @_inbounds getvalues(idx[i], arrays...))
                end
            end

            delta = 1
            while delta < WARP_SIZE
                acc = _combine(rf, acc, _shfl_down(acc, delta))
                delta <<= 1
            end
        end
    end
    # @cuprintf("%03ld: acc = %f\n", threadIdx().x, acc)

    dest === nothing && return

    # Preparing for block-wide merge:
    @assert nwarps_per_block <= 32
    if isbitstype(T)
        shared = @cuStaticSharedMem(T, 32)
    else
        S = UnionArrays.buffereltypefor(T)
        data = @cuStaticSharedMem(S, 32)
        typeids = @cuStaticSharedMem(UInt8, 32)
        @assert UInt(pointer(data, length(data) + 1)) == UInt(pointer(typeids))
        shared = UnionVector(T, data, typeids)
    end
    if warp_offset == 0
        @_inbounds shared[warpIdx] = acc
    end

    shared_bound =
        let n = length(idx),
            nbasecases = cld(n, basesize),
            offsetb = (blockIdx().x - 1) * nwarps_per_block,
            input_bound = nbasecases - offsetb

            min(input_bound, nwarps_per_block)
        end

    # Block-wide merge:
    shared_delta = 1
    while shared_delta < nwarps_per_block

        # Gather `WARP_SIZE` elements:
        i = warp_leader + shared_delta * warp_offset
        acc = start(rf, init)

        sync_threads()
        if i <= shared_bound
            acc = @_inbounds shared[i]
        end

        # Warp-wide merge:
        delta = 1
        while delta < WARP_SIZE
            acc = _combine(rf, acc, _shfl_down(acc, delta))
            delta <<= 1
        end

        if warp_offset == 0 && threadIdx().x <= lastindex(shared)
            @_inbounds shared[threadIdx().x] = acc
        end

        shared_delta *= WARP_SIZE
    end

    if threadIdx().x == 1
        @_inbounds dest[blockIdx().x] = acc
    end

    return
end


text:
module LevyArea

# Imports
import LinearAlgebra: mul!
import Random: randn!, default_rng
import SpecialFunctions: trigamma

"""
    abstract type AbstractIteratedIntegralAlgorithm end

Abstract type for algorithms for the simulation of iterated integrals.

```jldoctest; setup=:(using InteractiveUtils; using LevyArea)
julia> subtypes(AbstractIteratedIntegralAlgorithm)
4-element Vector{Any}:
 Fourier
 Milstein
 MronRoe
 Wiktorsson
```
"""
abstract type AbstractIteratedIntegralAlgorithm end

# Define error criteria
include("error_norms.jl")
export MaxL2
export FrobeniusL2

# Simulate in terms of Levy Area
include("iterated_integrals.jl")
export iterated_integrals

# Properties of Levy Area algorithms
include("alg_utils.jl")
export terms_needed
export optimal_algorithm

# Levy Area algorithms
include("fourier.jl")
include("milstein.jl")
include("wiktorsson.jl")
include("mronroe.jl")
export AbstractIteratedIntegralAlgorithm
export Fourier
export Milstein
export Wiktorsson
export MronRoe

const ITER_INT_ALGS = [Fourier(),Milstein(),Wiktorsson(),MronRoe()]

end # module

max_stars_repo_name: mrbauff/MasterThesis

max_issues_repo_name: mrbauff/MasterThesis


text:
# This file tests if the number of f and g calls are correctly counted
@testset "f_g_counts" begin

    fcalls = Ref(0)
    gcalls = Ref(0)

    function f_counts_ref!(x::Vector, fvec::Vector, fcalls)
        fcalls[] += 1
        fvec[1] = 1 - x[1]
        fvec[2] = 10(x[2]-x[1]^2)
    end
    f_counts!(fvec, x) = f_counts_ref!(x, fvec, fcalls)

    function g_counts_ref!(x::Vector, fjac::Matrix, gcalls)
        gcalls[] += 1
        fjac[1,1] = -1
        fjac[1,2] = 0
        fjac[2,1] = -20x[1]
        fjac[2,2] = 10
    end
    g_counts!(fjac, x) = g_counts_ref!(x, fjac, gcalls)

    x0 = [-1.2; 1.]

    df = OnceDifferentiable(f_counts!, g_counts!, similar(x0), similar(x0))

    fcalls[] = 0
    gcalls[] = 0

    r = nlsolve(df, x0, method = :trust_region)
    @test r.f_calls == fcalls[]
    @test r.g_calls == gcalls[]

    df = OnceDifferentiable(f_counts!, g_counts!, similar(x0), similar(x0))

    fcalls[] = 0
    gcalls[] = 0

    r = nlsolve(df, x0, method = :newton)
    @test r.f_calls == fcalls[]
    @test r.g_calls == gcalls[]
end


text:
using DataFrames, InteractiveFixedEffectModels,Test, CSV
df = df = CSV.read(joinpath(dirname(pathof(InteractiveFixedEffectModels)), "../dataset/Cigar.csv"))

for method in [:gauss_seidel,  :dogleg, :levenberg_marquardt]
	println(method)
	model = @formula Sales ~ 0 + ife(State, Year, 1)
	result = regife(df, model, method = method, maxiter = 10_000)
	@test abs(result.augmentdf[1, :factors1]) ≈ 0.18662770198472406 atol = 1e-1
	@test abs(result.augmentdf[1, :loadings1]) ≈ 587.2272 atol = 1e-1
	@test result.augmentdf[1, :residuals] ≈ - 15.6928 atol = 1e-1

	model = @formula Sales ~ 0 + ife(State, Year, 2)
	result = regife(df, model, method = method, maxiter = 10_000)
	@test abs(result.augmentdf[1, :factors1])  ≈ 0.18662770198472406 atol = 1e-1
	@test abs(result.augmentdf[1, :loadings1])  ≈ 587.227 atol = 1e-1
	@test result.augmentdf[1, :residuals] ≈ 2.16611 atol = 1e-1

	model = @formula Sales ~ ife(State, Year, 1) + fe(State)
	result = regife(df, model, method = method, maxiter = 10_000)
	@test abs(result.augmentdf[1, :factors1])   ≈ 0.17636 atol = 1e-1
	@test abs(result.augmentdf[1, :loadings1])  ≈ 20.176432452716522 atol = 1e-1
	@test result.augmentdf[1, :residuals]  ≈ - 10.0181 atol = 1e-1
	@test result.augmentdf[1, :fe_State]  ≈ 107.4766 atol = 1e-1

	model = @formula Sales ~ ife(State, Year, 2) + fe(State)
	result = regife(df, model, method = method, maxiter = 10_000)
	@test abs(result.augmentdf[1, :factors2]) ≈ 0.244 atol = 1e-1
	@test abs(result.augmentdf[1, :loadings2]) ≈ 49.7943 atol = 1e-1
	@test result.augmentdf[1, :residuals]  ≈ 2.165319 atol = 1e-1
	@test result.augmentdf[1, :fe_State]  ≈ 107.47666 atol = 1e-1

	model = @formula Sales ~ ife(State, Year, 2) + fe(State)
	result = regife(df, model, subset = df.State .<= 30, method = method, maxiter = 10_000)
	@test size(result.augmentdf, 1) == size(df, 1)
	@test abs(result.augmentdf[1, :factors1])  ≈ 0.20215 atol = 1e-1
	@test abs(result.augmentdf[1, :loadings1]) ≈ 29.546 atol = 1e-1
	@test abs(result.augmentdf[1, :factors2])  ≈ 0.250768  atol = 1e-1
	@test abs(result.augmentdf[1, :loadings2])   ≈ 43.578 atol = 1e-1
	@test result.augmentdf[1, :residuals]  ≈ 2.9448  atol = 1e-1
	@test result.augmentdf[1, :fe_State] ≈ 107.4766 atol = 1e-1
end


text:
function apply_row(func::Function, w::Array{Float64,2})
      δ = zeros(size(w,1));
      for (index, value) in enumerate(w[:,1])
           δ[index] = func(w[index,:]);
      end
      return(δ);
end

using Optim
using Optim: converged, maximum, maximizer, minimizer, iterations #some extra functions
# Generate IV regression data
begin "Test Optim"
      n = 3000;
      sig = [01.0 0.5 0.5;
             0.0 0.5 0.0 ;
             0.0 0.0 1.0];

      mu = [1,9,0];
      mn = MvNormal(mu,sig * sig')
      W = rand(mn,n)'
      x = W[:,1];z=W[:,2];ϵ=W[:,3];
      y = x * 3 + ϵ;
      β_OLS = inv(x'*x)*(x'*y);
      β_IV = inv(z'*x)*(z'*y);
      w = hcat(x,y,z);
      obj = θ->(begin
            η = (w[:,1] .* θ - w[:,2]) .* w[:,3];
            η'*η
            end)
      result = optimize(obj,zeros(1),Newton(),autodiff = :forward)
      # @show converged(result) || error("Failed to converge in $(iterations(result)) iterations")
      @show xmin = result.minimizer
      @show result.minimum
end

@testset begin "GMM"
      @test β_OLS - result.minimizer[1] < 0.1

end


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

valshape(d::Distributions.UnivariateDistribution) = ScalarShape{Real}()
valshape(d::Distributions.MultivariateDistribution) = ArrayShape{Real}(size(d)...)
valshape(d::Distributions.MatrixDistribution) = ArrayShape{Real}(size(d)...)


text:
using Pkg
Pkg.activate(pwd())

# # Linear regression
# ## Loading and preparing data

using Plots
using StatsPlots
using RDatasets

iris = dataset("datasets", "iris")
iris[1:5,:]

# ### Exercise:
# We will simplify the goal and estimate the dependence of petal width on petal length. # Create the data $X$ (do not forget to add the bias) and the labels $y$.
#
# Make a graph of the dependence of petal width on petal length.
#
# ---
# ### Solution:

# ---
#
# ## Training the classifier
# ### Exercise:
# Use the closed-form formula to get the coefficients $w$ for the linear regression. Then
# use the `optim` method derived in the previous lecture to solve the optimization problem
# via gradient descent. The results should be identical.
#
#---
#### Solution:

abstract type Step end

struct GD <: Step
    α::Real
end

optim_step(s::GD, f, g, x) = -s.α*g(x)

function optim(f, g, x, s::Step; max_iter=100)
    for i in 1:max_iter
        x += optim_step(s, f, g, x)
    end
    return x
end

# ---
#
# ### Exercise:
# Write the dependence on the petal width on the petal length. Plot it in the previous
# graph.
#
# ---
# ### Solution:


text:
using FMM3D
using Random
using LinearAlgebra

#  This is a sample code to demonstrate how to use
#  the fmm libraries
#

# sample with one density, sources to sources,
# charge interactions, and potential only
#
n = 2000
nd = 1
sources = rand(3,n)
eps = 10.0^(-5)

zk = 1.1 + im*0
charges = rand(n) + im*rand(n)
vals = hfmm3d(eps,zk,sources,charges=charges,pg=1)

# sample with a vector of densities, sources to
# sources and targets, dipole interactions,
# potential and gradietns

nd = 3
zk = 1.1 + im*0
nt = 1870
targ = rand(3,nt)
dipvecs = rand(nd,3,n) + im*rand(nd,3,n)
vals = hfmm3d(eps,zk,sources,dipvecs=dipvecs,
             targets=targ,nd=nd,pg=2,pgt=2)


text:
#
# Constant diffusivity
#

"""
    ConstantIsotropicDiffusivity(T=Float64; ν=1e-6, κ=1e-7)

or

    MolecularDiffusivity(T=Float64; ν=1e-6, κ=1e-7)

Return a `ConstantIsotropicDiffusivity` closure object of type `T` with
viscosity `ν` and scalar diffusivity `κ`.
"""
Base.@kwdef struct ConstantIsotropicDiffusivity{T} <: IsotropicDiffusivity{T}
    ν :: T = 1e-6
    κ :: T = 1e-7
end

const MolecularDiffusivity = ConstantIsotropicDiffusivity

ConstantIsotropicDiffusivity(T; kwargs...) =
    typed_keyword_constructor(T, ConstantIsotropicDiffusivity; kwargs...)

κ_ccc(i, j, k, grid, closure::ConstantIsotropicDiffusivity, args...) = closure.κ

∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, eos, g, u, v, w, T, S) = (
      closure.ν * ∂x²_faa(i, j, k, grid, u)
    + closure.ν * ∂y²_aca(i, j, k, grid, u)
    + closure.ν * ∂z²_aac(i, j, k, grid, u)
    )

∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, eos, g, u, v, w, T, S) = (
      closure.ν * ∂x²_caa(i, j, k, grid, v)
    + closure.ν * ∂y²_afa(i, j, k, grid, v)
    + closure.ν * ∂z²_aac(i, j, k, grid, v)
    )

∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, eos, g, u, v, w, T, S) = (
      closure.ν * ∂x²_caa(i, j, k, grid, w)
    + closure.ν * ∂y²_aca(i, j, k, grid, w)
    + closure.ν * ∂z²_aaf(i, j, k, grid, w)
    )

#
# Constant anisotropic diffusivity (tensor diffusivity with heterogeneous
#                                       diagonal components)
#

"""
    ConstantAnisotropicDiffusivity(T=Float64; νh=1e-6, νv=1e-6, κh=1e-6, κv=1e-6)

Returns a ConstantAnisotropicDiffusivity object with horizontal viscosity and
diffusivity `νh` and `κh`, and vertical viscosity and diffusivity
`νv` and `κv`.
"""
Base.@kwdef struct ConstantAnisotropicDiffusivity{T} <: TensorDiffusivity{T}
    νh :: T = 1e-6
    νv :: T = 1e-6
    κh :: T = 1e-6
    κv :: T = 1e-6
end

ConstantAnisotropicDiffusivity(T; kwargs...) =
    typed_keyword_constructor(T, ConstantAnisotropicDiffusivity; kwargs...)

∇_κ_∇ϕ(i, j, k, grid, ϕ, closure::ConstantAnisotropicDiffusivity, args...) = (
      ∂x_κ_∂x_ϕ(i, j, k, grid, ϕ, closure.κh, closure, args...)
    + ∂y_κ_∂y_ϕ(i, j, k, grid, ϕ, closure.κh, closure, args...)
    + ∂z_κ_∂z_ϕ(i, j, k, grid, ϕ, closure.κv, closure, args...)
    )

∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure::ConstantAnisotropicDiffusivity, eos, g, u, v, w, T, S) = (
      closure.νh * ∂x²_faa(i, j, k, grid, u)
    + closure.νh * ∂y²_aca(i, j, k, grid, u)
    + closure.νv * ∂z²_aac(i, j, k, grid, u)
    )

∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure::ConstantAnisotropicDiffusivity, eos, g, u, v, w, T, S) = (
      closure.νh * ∂x²_caa(i, j, k, grid, v)
    + closure.νh * ∂y²_afa(i, j, k, grid, v)
    + closure.νv * ∂z²_aac(i, j, k, grid, v)
    )

∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure::ConstantAnisotropicDiffusivity, eos, g, u, v, w, T, S) = (
      closure.νh * ∂x²_caa(i, j, k, grid, w)
    + closure.νh * ∂y²_aca(i, j, k, grid, w)
    + closure.νv * ∂z²_aaf(i, j, k, grid, w)
    )


text:
struct ERBFilterbank{C,G,T<:Real,U<:Real,V<:Real} <: Filterbank
    filters::Vector{SecondOrderSections{C,G}}
    ERB::Vector{T}
    center_frequencies::Vector{U}
    fs::V
end

function make_erb_filterbank(fs, num_channels, low_freq, EarQ = 9.26449, minBW = 24.7, order = 1)
    T = 1/fs
    if length(num_channels) == 1
        cf = erb_space(low_freq, fs/2, num_channels)
    else
        cf = num_channels
        if size(cf,2) > size(cf,1)
        	cf = cf'
        end
    end
    ERB = ((cf/EarQ).^order .+ minBW.^order).^(1/order)
    B = 1.019*2*pi*ERB
    B0 = T
    B2 = 0.0
    A0 = 1.0
    A1 = -2*cos.(2*cf*pi*T)./exp.(B*T)
    A2 = exp.(-2*B*T)

    B11 = -(2*T*cos.(2*cf*pi*T)./exp.(B*T) .+ 2*sqrt(3+2^1.5)*T*sin.(2*cf*pi*T)./exp.(B*T))/2
    B12 = -(2*T*cos.(2*cf*pi*T)./exp.(B*T) .- 2*sqrt(3+2^1.5)*T*sin.(2*cf*pi*T)./exp.(B*T))/2
    B13 = -(2*T*cos.(2*cf*pi*T)./exp.(B*T) .+ 2*sqrt(3-2^1.5)*T*sin.(2*cf*pi*T)./exp.(B*T))/2
    B14 = -(2*T*cos.(2*cf*pi*T)./exp.(B*T) .- 2*sqrt(3-2^1.5)*T*sin.(2*cf*pi*T)./exp.(B*T))/2

    gain = abs.((-2*exp.(4*im*cf*pi*T)*T .+ 2*exp.(-(B*T) .+
      2*im*cf*pi*T).*T.*(cos.(2*cf*pi*T) .- sqrt(3 - 2^(3/2))*
      sin.(2*cf*pi*T))) .* (-2*exp.(4*im*cf*pi*T)*T .+ 2*exp.(-(B*T) .+
      2*im*cf*pi*T).*T.* (cos.(2*cf*pi*T) .+ sqrt(3 - 2^(3/2)) *
      sin.(2*cf*pi*T))).* (-2*exp.(4*im*cf*pi*T)*T .+ 2*exp.(-(B*T) .+
      2*im*cf*pi*T).*T.* (cos.(2*cf*pi*T) - sqrt(3 +
      2^(3/2))*sin.(2*cf*pi*T))) .* (-2*exp.(4*im*cf*pi*T)*T .+
      2*exp.(-(B*T) + 2*im*cf*pi*T).*T.* (cos.(2*cf*pi*T) + sqrt(3 +
      2^(3/2))*sin.(2*cf*pi*T))) ./ (-2 ./ exp.(2*B*T) -
      2*exp.(4*im*cf*pi*T) + 2*(1 .+ exp.(4*im*cf*pi*T))./exp.(B*T)).^4)

    C = typeof(B0)
    filters = Array{SOSFilter{C,C}}(undef,num_channels)
    for ch=1:num_channels
        biquads = Array{BiquadFilter{C}}(undef,4)
        biquads[1] = BiquadFilter(B0, B11[ch], B2, A0, A1[ch], A2[ch])
        biquads[2] = BiquadFilter(B0, B12[ch], B2, A0, A1[ch], A2[ch])
        biquads[3] = BiquadFilter(B0, B13[ch], B2, A0, A1[ch], A2[ch])
        biquads[4] = BiquadFilter(B0, B14[ch], B2, A0, A1[ch], A2[ch])
        filters[ch] = SOSFilter(biquads, 1/gain[ch])
    end
    ERBFilterbank(filters, ERB, cf, fs)
end

function erb_space(low_freq, high_freq, num_channels, EarQ = 9.26449, minBW = 24.7, order = 1)
    # All of the following expressions are derived in Apple TR #35, "An
    # Efficient Implementation of the Patterson-Holdsworth Cochlear
    # Filter Bank."  See pages 33-34.
    space = (1:num_channels).*(-log(high_freq + EarQ*minBW) + log(low_freq + EarQ*minBW))/num_channels
    cfArray = -(EarQ*minBW) .+ exp.(space) * (high_freq + EarQ*minBW)
end


text:
Base.include(Main, joinpath(@__DIR__, "startup.jl"))

command_list = copy(ARGS)
if isempty(command_list)
    command_list = ["add"]
end
for (i, command) in enumerate(command_list)
    if command ∉ ["add", "build", "compile"]
        error("Unknown $i-th command: $command")
    end
end

@info """
Base.DEPOT_PATH is initialized to:
$(sprint(show, "text/plain", Base.DEPOT_PATH))
"""

ENV["PYTHON"] = ENV["JULIA_VENV_PYTHON"]
import Pkg

let venv_depot = joinpath(@__DIR__, "depot")
    for command in command_list
        if command == "add"
            if "PyCall" ∉ keys(Pkg.installed())
                @info "Installing PyCall.jl..."
                @assert Base.DEPOT_PATH[1] == venv_depot
                Pkg.add("PyCall")
            else
                @info "PyCall.jl is already installed"
            end
        elseif command == "build"
            Pkg.build("PyCall")
        elseif command == "compile"
            # Manually removing PyCall's precmopile file since it may
            # not be detected as "stale" by julia when it was compiled
            # in non-PyJulia -mode.
            pkg = Base.identify_package("PyCall")
            if pkg !== nothing
                cache_files = Base.find_all_in_cache_path(pkg)
                if !isempty(cache_files) &&
                        startswith(cache_files[1], venv_depot)
                    @info """
                    Removing old compilation cache for PyCall:
                    $(cache_files[1])
                    """
                    rm(cache_files[1], force = true)
                end
            end
            @info "Compiling PyCall..."
            JuliaVenv.precompile_package("PyCall")
        else
            error("Unknown command $command")
        end
    end
end

max_stars_repo_name: tkf/julia-venv

max_issues_repo_name: tkf/julia-venv

max_forks_repo_name: tkf/julia-venv


text:
function horizontal_average_is_correct(arch, FT)
    model = BasicModel(N = (16, 16, 16), L = (100, 100, 100), architecture=arch, float_type=FT)

    # Set a linear stably stratified temperature profile.
    T₀(x, y, z) = 20 + 0.01*z
    set!(model; T=T₀)

    T̅ = HorizontalAverage(model, model.tracers.T; interval=0.5second)
    push!(model.diagnostics, T̅)

    time_step!(model, 1, 1)
    correct_profile = @. 20 + 0.01 * collect(model.grid.zC)
    all(Array(T̅.profile[:][2:end-1]) ≈ correct_profile)
end

function product_profile_is_correct(arch, FT)
    model = BasicModel(N = (16, 16, 16), L = (100, 100, 100), architecture=arch, float_type=FT)

    # Set a linear stably stratified temperature profile and a sinusoidal u(z) profile.
    u₀(x, y, z) = sin(z)
    T₀(x, y, z) = 20 + 0.01*z
    set!(model; u=u₀, T=T₀)

    uT = HorizontalAverage(model, [model.velocities.u, model.tracers.T]; interval=0.5second)
    run_diagnostic(model, uT)

    correct_profile = @. sin.(model.grid.zC) * (20 + 0.01 * model.grid.zC)
    Array(uT.profile[:][2:end-1]) ≈ correct_profile
end

function nan_checker_aborts_simulation(arch, FT)
    model = BasicModel(N = (16, 16, 2), L = (1, 1, 1), architecture=arch, float_type=FT)

    # It checks for NaNs in w by default.
    nc = NaNChecker(model; frequency=1, fields=Dict(:w => model.velocities.w.data.parent))
    push!(model.diagnostics, nc)

    model.velocities.w[4, 3, 2] = NaN

    time_step!(model, 1, 1);
end

TestModel(::GPU, FT, ν=1.0, Δx=0.5) = BasicModel(N=(16, 16, 16), L=(16*Δx, 16*Δx, 16*Δx),
                                                 architecture=GPU(), float_type=FT, ν=ν, κ=ν)

TestModel(::CPU, FT, ν=1.0, Δx=0.5) = BasicModel(N=(3, 3, 3), L=(3*Δx, 3*Δx, 3*Δx),
                                                 architecture=CPU(), float_type=FT, ν=ν, κ=ν)

function max_abs_field_diagnostic_is_correct(arch, FT)
    model = TestModel(arch, FT)
    set!(model.velocities.u, rand(size(model.grid)))
    u_max = FieldMaximum(abs, model.velocities.u)
    return u_max(model) == maximum(abs, model.velocities.u.data.parent)
end

function advective_cfl_diagnostic_is_correct(arch, FT)
    model = TestModel(arch, FT)

    Δt = FT(1.3e-6)
    Δx = FT(model.grid.Δx)
    u₀ = FT(1.2)
    CFL_by_hand = Δt * u₀ / Δx

    model.velocities.u.data.parent .= u₀
    cfl = AdvectiveCFL(FT(Δt))

    return cfl(model) ≈ CFL_by_hand
end

function diffusive_cfl_diagnostic_is_correct(arch, FT)
    Δt = FT(1.3e-6)
    Δx = FT(0.5)
    ν = FT(1.2)
    CFL_by_hand = Δt * ν / Δx^2

    model = TestModel(arch, FT, ν, Δx)
    cfl = DiffusiveCFL(FT(Δt))

    return cfl(model) ≈ CFL_by_hand
end

get_iteration(model) = model.clock.iteration
get_time(model) = model.clock.time

function timeseries_diagnostic_works(arch, FT)
    model = TestModel(arch, FT)
    iter_diag = Timeseries(get_iteration, model; frequency=1)
    push!(model.diagnostics, iter_diag)
    Δt = FT(1e-16)
    time_step!(model, 1, Δt)

    return iter_diag.time[end] == Δt && iter_diag.data[end] == 1
end

function timeseries_diagnostic_tuples(arch, FT)
    model = TestModel(arch, FT)
    timeseries = Timeseries((iters=get_iteration, itertimes=get_time), model; frequency=2)
    model.diagnostics[:timeseries] = timeseries
    Δt = FT(1e-16)
    time_step!(model, 2, Δt)
    return timeseries.iters[end] == 2 && timeseries.itertimes[end] == 2Δt
end

function diagnostics_getindex(arch, FT)
    model = TestModel(arch, FT)
    iter_timeseries = Timeseries(get_iteration, model)
    time_timeseries = Timeseries(get_time, model)
    model.diagnostics[:iters] = iter_timeseries
    model.diagnostics[:times] = time_timeseries
    return model.diagnostics[2] == time_timeseries
end

function diagnostics_setindex(arch, FT)
    model = TestModel(arch, FT)
    iter_timeseries = Timeseries(get_iteration, model)
    time_timeseries = Timeseries(get_time, model)
    max_abs_u_timeseries = Timeseries(FieldMaximum(abs, model.velocities.u), model; frequency=1)

    push!(model.diagnostics, iter_timeseries, time_timeseries)
    model.diagnostics[2] = max_abs_u_timeseries

    return model.diagnostics[:diag2] == max_abs_u_timeseries
end

@testset "Diagnostics" begin
    println("Testing diagnostics...")

    for arch in archs
        @testset "Horizontal average [$(typeof(arch))]" begin
            println("  Testing horizontal average [$(typeof(arch))]")
            for FT in float_types
                @test horizontal_average_is_correct(arch, FT)
                @test product_profile_is_correct(arch, FT)
            end
        end
    end

    for arch in archs
        @testset "NaN Checker [$(typeof(arch))]" begin
            println("  Testing NaN Checker [$(typeof(arch))]")
            for FT in float_types
                @test_throws ErrorException nan_checker_aborts_simulation(arch, FT)
            end
        end
    end

    for arch in archs
        @testset "Miscellaneous timeseries diagnostics [$(typeof(arch))]" begin
            println("  Testing miscellaneous timeseries diagnostics [$(typeof(arch))]")
            for FT in float_types
                @test diffusive_cfl_diagnostic_is_correct(arch, FT)
                @test advective_cfl_diagnostic_is_correct(arch, FT)
                @test max_abs_field_diagnostic_is_correct(arch, FT)
                @test timeseries_diagnostic_works(arch, FT)
                @test timeseries_diagnostic_tuples(arch, FT)
                @test diagnostics_getindex(arch, FT)
                @test diagnostics_setindex(arch, FT)
            end
        end
    end
end


text:
# Copyright 2017 Stephen L. Smith and Frank Imeson
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
module GLNS
export solver
using Random
include("src/utilities.jl")
include("src/parse_print.jl")
include("src/tour_optimizations.jl")
include("src/adaptive_powers.jl")
include("src/insertion_deletion.jl")
include("src/parameter_defaults.jl")

"""
Main GTSP solver, which takes as input a problem instance and
some optional arguments
"""
function solver(problem_instance; args...)
	###### Read problem data and solver settings ########
	num_vertices, num_sets, sets, dist, membership = read_file(problem_instance)
	param = parameter_settings(num_vertices, num_sets, sets, problem_instance, args)
	#####################################################
	init_time = time()
	count = Dict(:latest_improvement => 1,
	  			 :first_improvement => false,
	 		     :warm_trial => 0,
	  		     :cold_trial => 1,
				 :total_iter => 0,
				 :print_time => init_time)
	lowest = Tour(Int64[], typemax(Int64))
	start_time = time_ns()
	# compute set distances which will be helpful
	setdist = set_vertex_dist(dist, num_sets, membership)
	powers = initialize_powers(param)

	while count[:cold_trial] <= param[:cold_trials]
		# build tour from scratch on a cold restart
		best = initial_tour!(lowest, dist, sets, setdist, count[:cold_trial], param)
		# print_cold_trial(count, param, best)
		phase = :early

		if count[:cold_trial] == 1
			powers = initialize_powers(param)
		else
			power_update!(powers, param)
		end

		while count[:warm_trial] <= param[:warm_trials]
			iter_count = 1
			current = Tour(copy(best.tour), best.cost)
			temperature = 1.442 * param[:accept_percentage] * best.cost
			# accept a solution with 50% higher cost with 0.05% change after num_iterations.
			cooling_rate = ((0.0005 * lowest.cost)/(param[:accept_percentage] *
									current.cost))^(1/param[:num_iterations])

			if count[:warm_trial]  > 0	  # if warm restart, then use lower temperature
		        temperature *= cooling_rate^(param[:num_iterations]/2)
				phase = :late
			end
			while count[:latest_improvement] <= (count[:first_improvement] ?
				  param[:latest_improvement] : param[:first_improvement])

				if iter_count > param[:num_iterations]/2 && phase == :early
					phase = :mid  # move to mid phase after half iterations
				end
				trial = remove_insert(current, best, dist, membership, setdist, sets, powers, param, phase)

		        # decide whether or not to accept trial
				if accepttrial_noparam(trial.cost, current.cost, param[:prob_accept]) ||
				   accepttrial(trial.cost, current.cost, temperature)
					param[:mode] == "slow" && opt_cycle!(current, dist, sets, membership, param, setdist, "full")
				    current = trial
		        end
		        if current.cost < best.cost
					count[:latest_improvement] = 1
					count[:first_improvement] = true
					if count[:cold_trial] > 1 && count[:warm_trial] > 1
						count[:warm_trial] = 1
					end
					opt_cycle!(current, dist, sets, membership, param, setdist, "full")
					best = current
	        	else
					count[:latest_improvement] += 1
				end

				# if we've come in under budget, or we're out of time, then exit
			    if best.cost <= param[:budget] || time() - init_time > param[:max_time]
					param[:timeout] = (time() - init_time > param[:max_time])
					param[:budget_met] = (best.cost <= param[:budget])
					timer = (time_ns() - start_time)/1.0e9
					lowest.cost > best.cost && (lowest = best)
					print_best(count, param, best, lowest, init_time)
					print_summary(lowest, timer, membership, param)
					return
				end

		        temperature *= cooling_rate  # cool the temperature
				iter_count += 1
				count[:total_iter] += 1
				print_best(count, param, best, lowest, init_time)
			end
			print_warm_trial(count, param, best, iter_count)
			# on the first cold trial, we are just determining
			count[:warm_trial] += 1
			count[:latest_improvement] = 1
			count[:first_improvement] = false
		end
		lowest.cost > best.cost && (lowest = best)
		count[:warm_trial] = 0
		count[:cold_trial] += 1

		# print_powers(powers)

	end
	timer = (time_ns() - start_time)/1.0e9
	print_summary(lowest, timer, membership, param)
	return [lowest.cost, timer]
end
end

max_stars_repo_name: Azizimj/GTSPbnd

max_issues_repo_name: Azizimj/GTSPbnd

max_forks_repo_name: Azizimj/GTSPbnd


text:
using RefractiveIndex
using Test

function midrange(material)
    λmin, λmax = material.λrange
    return λmin + 0.5(λmax - λmin)
end

function testRM(material, n_ref)
    n = material(midrange(material))
    # Compare only fractional parts
    isapprox(n_ref % 1, n % 1, rtol=1e-3)
end

@testset "RefractiveIndex.jl" begin
    @testset "Dispersion formulas" begin
        # Sellmeier
        @test testRM(RefractiveMaterial("main", "Ar", "Grace-liquid-90K"), 1.2281)

        # Sellmeier-2
        @test testRM(RefractiveMaterial("main", "CdTe", "Marple"), 2.7273)

        # Polynomial
        @test testRM(RefractiveMaterial("other", "Mg-LiTaO3", "Moutzouris-o"), 2.1337)

        # RefractiveIndex.INFO
        @test testRM(RefractiveMaterial("main", "ZnTe", "Li"), 2.6605)

        # Cauchy
        @test testRM(RefractiveMaterial("main", "SF6", "Vukovic"), 1.00072071)

        # Gases
        @test testRM(RefractiveMaterial("main", "He", "Mansfield"), 1.000034724)

        # Herzberger
        @test testRM(RefractiveMaterial("main", "Si", "Edwards"), 3.4208)

        # Retro
        @test testRM(RefractiveMaterial("main", "AgBr", "Schröter"), 2.2600)

        # Exotic
        @test testRM(RefractiveMaterial("organic", "urea","Rosker-e"), 1.6000)
    end

    @testset "Tabular data" begin
        # RefractiveNK
        @test testRM(RefractiveMaterial("main", "ZnO", "Stelling"), 1.5970)
    end
end


text:
@testset "Test delay & refractoriness" begin
    config = """
    spike_duration: 4.99
    inputs:
    - id: i1
    outputs:
    - id: o1
    neurons:
    - id: n1
    synapses:
    - {id: syn1, source: i1, target: n1}
    - {id: syn2, source: n1, target: o1, spike_duration: 0}
    """

    (net,objects) = load_network(YAML_source=config)

    inp=[Event(:input_spikes, 0.0, 1.0, objects[:i1])]

    logger=simulate!(net, inp, 10.0; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:o1))
    # should spike once every 1ms for the whole spike duration of 5ms with 2ms delay
    # (due to the 2 synapses) after the input spike at 1ms, i.e. after 3ms
    @test length(logger.data.t)==5 && logger.data.t ≈ [3.0,4.0,5.0,6.0,7.0]

    # if the refractoriness = 5.0ms > the duration of a spike, the neuron should only spike once
    config = """
    spike_duration: 4.99
    inputs:
    - id: i1
    outputs:
    - id: o1
    neurons:
    - id: n1
      refractory_duration: 5.0 # <-----------------
    synapses:
    - {id: syn1, source: i1, target: n1}
    - {id: syn2, source: n1, target: o1, spike_duration: 0}
    """

    (net,objects) = load_network(YAML_source=config)
    inp=[Event(:input_spikes, 0.0, 1.0, objects[:i1])]
    logger=simulate!(net, inp, 10.0; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:o1))
    @test length(logger.data.t)==1 && logger.data.t ≈ [3.0]
end

@testset "Reset test" begin
    config = """
    refractory_duration: 10
    plateau_duration: 10
    spike_duration: 10
    inputs:
    - id: i1
    outputs:
    - id: o1
    neurons:
    - id: n1
      branches:
        - id: s1
          branches:
            - id: s2
              branches:
                - id: s3
    synapses:
    - {id: syn1, source: i1, target: n1}
    - {id: syn2, source: i1, target: s1}
    - {id: syn3, source: i1, target: s2}
    - {id: syn4, source: i1, target: s3}
    - {id: syn5, source: n1, target: o1, spike_duration: 0}
    """

    (net,objects) = load_network(YAML_source=config)
    inp=[Event(:input_spikes, 0.0, 1.0, objects[:i1])]
    logger=simulate!(net, inp, 10.0; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:o1))
    @test length(logger.data.t)==1 && logger.data.t ≈ [3.0]
    logger=simulate!(net, inp, 10.0; show_progress=false, reset=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:o1))
    @test isempty(logger.data.t)
    logger=simulate!(net, inp, 10.0; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:o1))
    @test length(logger.data.t)==1 && logger.data.t ≈ [3.0]
end

@testset "Test cascades" begin
  config = """
  spike_duration: 0.99
  plateau_duration: 5.0
  inputs:
  - id: i1
  - id: i2
  - id: i3
  - id: i4
  outputs:
  - id: o1
  neurons:
  - id: n1
    branches:
      - id: s1
        branches:
          - id: s2
            branches:
              - id: s3
  synapses:
  - {id: syn1, source: i1, target: n1}
  - {id: syn2, source: i2, target: s1}
  - {id: syn3, source: i3, target: s2}
  - {id: syn4, source: i4, target: s3}
  - {id: syn5, source: n1, target: o1, spike_duration: 0}
  """

  (net,objects) = load_network(YAML_source=config)

  cascade_trials=[
    (
      summary="Slow sequence in wrong direction",
      input=[
        Event(:input_spikes, 0.0, 1.0, objects[:i1]), # should do nothing
        Event(:input_spikes, 0.0, 2.0, objects[:i2]), # should do nothing
        Event(:input_spikes, 0.0, 3.0, objects[:i3]), # should do nothing
        Event(:input_spikes, 0.0, 4.0, objects[:i4]), # should trigger a plateau on s3, but no spike!
      ],
      output=[]
    ),
    (
      summary="Re-triggering the neuron",
      input=[
        Event(:input_spikes, 0.0, 0.0, objects[:i1]), # should do nothing
        Event(:input_spikes, 0.0, 0.2, objects[:i1]), # should do nothing
        Event(:input_spikes, 0.0, 0.4, objects[:i1]), # should do nothing
      ],
      output=[]
    ),
    (
      summary="Quick sequence in wrong direction",
      input=[
        Event(:input_spikes, 0.0, 0.0, objects[:i1]), # should do nothing
        Event(:input_spikes, 0.0, 0.2, objects[:i2]), # should do nothing
        Event(:input_spikes, 0.0, 0.4, objects[:i3]), # should do nothing
        Event(:input_spikes, 0.0, 0.6, objects[:i4]), # triggers a spike at 1.6 -> spike on syn5 at 2.6
      ],
      output=[2.6]
    ),
    (
      summary="Slow sequence in correct direction",
      input=[
        Event(:input_spikes, 0.0, 0.0, objects[:i4]), # should do nothing
        Event(:input_spikes, 0.0, 2.0, objects[:i3]), # should do nothing
        Event(:input_spikes, 0.0, 4.0, objects[:i2]), # should do nothing
        Event(:input_spikes, 0.0, 6.0, objects[:i1]), # trigger a spike at 7.0 -> spike on syn5 at 8.0
        ],
      output=[8.0]
    ),
    (
      summary="Too slow sequence in correct direction",
      input=[
        Event(:input_spikes, 0.0, 0.0, objects[:i4]), # should do nothing
        Event(:input_spikes, 0.0, 2.0, objects[:i3]), # should do nothing
        Event(:input_spikes, 0.0, 4.0, objects[:i2]), # should do nothing
        Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should do nothing
        ],
      output=[]
    )
  ]

  @testset "Test cascade: $(trial.summary)" for trial in cascade_trials
    logger=simulate!(net, trial.input; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:o1))
    @test length(logger.data.t) == length(trial.output) && logger.data.t ≈ Vector{Float64}(trial.output)
  end
end

@testset "Test multi-branch" begin
    config = """
    refractory_duration: 10
    inputs:
    - id: i1
    - id: i2
    - id: i3
    outputs:
    - id: o1
    neurons:
    - id: n1
      branches:
        - id: s1
        - id: s2
    synapses:
    - {id: syn1, source: i1, target: n1}
    - {id: syn2, source: i2, target: s1}
    - {id: syn3, source: i3, target: s2}
    - {id: syn4, source: n1, target: o1, spike_duration: 0}
    """

    (net,objects) = load_network(YAML_source=config)

    trials=[
      (
        summary="Branch 1",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i2]), # should do nothing
          Event(:input_spikes, 0.0, 1.0, objects[:i1]), # should trigger a spike!
        ],
        output=[3.0]
      ),
      (
        summary="Branch 2",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # should do nothing
          Event(:input_spikes, 0.0, 1.0, objects[:i1]), # should trigger a spike!
        ],
        output=[3.0]
      ),
      (
        summary="No branch",
        input=[
          Event(:input_spikes, 0.0, 1.0, objects[:i3]), # should do nothing
          Event(:input_spikes, 0.0, 2.0, objects[:i2]), # should not trigger a spike!
        ],
        output=[]
      ),
    ]

    @testset "Test cascade: $(trial.summary)" for trial in trials
        logger=simulate!(net, trial.input; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:o1))
        @test length(logger.data.t) == length(trial.output) && logger.data.t ≈ Vector{Float64}(trial.output)
    end
end

@testset "Test inhibition at soma" begin
    config = """
    spike_duration: 4.99
    inputs:
    - id: i1
    - id: i2
    - id: i3
    outputs:
    - id: o1
    neurons:
    - id: n1
      refractory_duration: 5.0 # <-----------------
    synapses:
    - {id: syn1, source: i1, target: n1}
    - {id: syn2, source: i2, target: n1, weight: -1}
    - {id: syn3, source: i3, target: n1, weight: -1, spike_duration: 10}
    - {id: syn4, source: n1, target: o1, spike_duration: 0}
    """

    (net,objects) = load_network(YAML_source=config)

    trials=[
      (
        summary="Inhibit after spike",
        input=[
          Event(:input_spikes, 0.0, 1.0, objects[:i1]),
          Event(:input_spikes, 0.0, 1.5, objects[:i2])
        ],
        output=[3.0]
      ),
      (
        summary="Inhibit before spike",
        input=[
          Event(:input_spikes, 0.0, 1.0, objects[:i2]),
          Event(:input_spikes, 0.0, 1.5, objects[:i1])
        ],
        output=[7.99]
      ),
      (
        summary="Inhibit with long EPSP",
        input=[
          Event(:input_spikes, 0.0, 1.0, objects[:i3]),
          Event(:input_spikes, 0.0, 1.5, objects[:i1])
        ],
        output=[]
      ),
      (
        summary="Inhibition can prevent rebound spike",
        input=[
          Event(:input_spikes, 0.0, 1.0, objects[:i3]),
          Event(:input_spikes, 0.0, 1.0, objects[:i2]),
          Event(:input_spikes, 0.0, 1.5, objects[:i1])
        ],
        output=[]
      ),
    ]
    @testset "$(trial.summary)" for trial in trials
      logger=simulate!(net, trial.input; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:o1))
      @test length(logger.data.t) == length(trial.output) && logger.data.t ≈ Vector{Float64}(trial.output)
    end
end

#####

@testset "Test inhibition on multiple branches" begin
    config = """
    refractory_duration: 10
    inputs:
    - id: i1
    - id: i2
    - id: i3
    - id: i4
    - id: i5
    - id: i6
    - id: i7
    - id: i8
    outputs:
    - id: o1
    neurons:
    - id: n1
      branches:
        - id: s0
          branches:
            - id: s1
            - id: s2
    synapses:
    - {id: syn1, source: i1, target: n1}
    - {id: syn2, source: i2, target: n1, weight: -1, spike_duration: 10}
    - {id: syn3, source: i3, target: s1}
    - {id: syn4, source: i4, target: s1, weight: -1}
    - {id: syn5, source: i5, target: s2}
    - {id: syn6, source: i6, target: s2, weight: -1}
    - {id: syn7, source: i7, target: s0}
    - {id: syn8, source: i8, target: s0, weight: -1}
    - {id: syn9, source: n1, target: o1, spike_duration: 0}
    """

    (net,objects) = load_network(YAML_source=config)

    trials=[
      (
        summary="Cascade",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # trigger plateau on s1
          Event(:input_spikes, 0.0, 7.5, objects[:i7]), # trigger plateau on s0
          Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should trigger a spike!
        ],
        output=[12.0]
      ),
      (
        summary="Cascade interrupted 1",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # trigger plateau on s1
          Event(:input_spikes, 0.0, 5.0, objects[:i4]), # end plateau on s1
          Event(:input_spikes, 0.0, 7.5, objects[:i7]), # fails to trigger plateau on s0
          Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should not trigger a spike!
        ],
        output=[]
      ),
      (
        summary="Cascade not interrupted 2",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # trigger plateau in s1
          Event(:input_spikes, 0.0, 5.0, objects[:i7]), # trigger plateau in s0
          Event(:input_spikes, 0.0, 7.5, objects[:i4]), # end plateau in top s1 --> doesn't turns off s0
          Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should trigger a spike!
        ],
        output=[12.0]
      ),
      (
        summary="Cascade interrupted 3",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # trigger plateau
          Event(:input_spikes, 0.0, 5.0, objects[:i7]), # trigger plateau
          Event(:input_spikes, 0.0, 7.5, objects[:i8]), # end plateau in lower segment
          Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should not trigger a spike!
        ],
        output=[]
      ),
      (
        summary="Cascade interrupted 4",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # trigger plateau
          Event(:input_spikes, 0.0, 5.0, objects[:i7]), # trigger plateau
          Event(:input_spikes, 0.0, 7.5, objects[:i2]), # inhibit soma
          Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should not trigger a spike!
        ],
        output=[]
      ),
      (
        summary="Cross-branch uninterrupted",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # trigger plateau on s1
          Event(:input_spikes, 0.0, 1.0, objects[:i5]), # trigger plateau on s2
          Event(:input_spikes, 0.0, 2.0, objects[:i7]), # trigger plateau on s0
          Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should trigger a spike!
        ],
        output=[12.0]
      ),
      (
        summary="Cross-branch interrupted",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # trigger plateau on s1
          Event(:input_spikes, 0.0, 1.0, objects[:i5]), # trigger plateau on s2
          Event(:input_spikes, 0.0, 2.0, objects[:i7]), # trigger plateau on s0
          Event(:input_spikes, 0.0, 3.0, objects[:i8]), # end plateau on s0
          Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should not trigger a spike!
        ],
        output=[]
      ),
      (
        summary="Cross-branch not interrupted 1",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # trigger plateau on s1
          Event(:input_spikes, 0.0, 1.0, objects[:i5]), # trigger plateau on s2
          Event(:input_spikes, 0.0, 2.0, objects[:i7]), # trigger plateau on s0
          Event(:input_spikes, 0.0, 3.0, objects[:i4]), # end plateau on s1
          Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should still trigger a spike!
        ],
        output=[12.0]
      ),
      (
        summary="Cross-branch not interrupted 2",
        input=[
          Event(:input_spikes, 0.0, 0.0, objects[:i3]), # trigger plateau on s1
          Event(:input_spikes, 0.0, 1.0, objects[:i5]), # trigger plateau on s2
          Event(:input_spikes, 0.0, 2.0, objects[:i7]), # trigger plateau on s0
          Event(:input_spikes, 0.0, 3.0, objects[:i4]), # end plateau on s1
          Event(:input_spikes, 0.0, 5.0, objects[:i6]), # end plateau on s2
          Event(:input_spikes, 0.0, 10.0, objects[:i1]), # should not trigger a spike!
        ],
        output=[12.0]
      )
    ]

    @testset "Test cascade: $(trial.summary)" for trial in trials
        logger=simulate!(net, trial.input; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:o1))
        @test length(logger.data.t) == length(trial.output) && logger.data.t ≈ Vector{Float64}(trial.output)
    end
end

@testset "Don't forget reset" begin
    #This test for a bug, where a plateau-reset would not happen if plateau conditions are still met at the time where
    #the plateau should end. No new end would be scheduled, therefore the segment would stay on indefinitely.

    config = """
    plateau_duration: 100
    refractory_duration: 5.01
    inputs:
    - id: i1
    - id: i2
    neurons:
    - id: n
      branches:
        - id: seg
    synapses:
    - {id: syn1, source: i1, target: seg}
    - {id: syn2, source: i2, target: n}
    """

    (net,objects) = load_network(YAML_source=config)
    input=[
      Event(:input_spikes, 0.0, 0.0, objects[:i1]), # trigger plateau on seg
      Event(:input_spikes, 0.0, 99.0, objects[:i1]), # satisfy plateau conditions at time of plateau-end
      Event(:input_spikes, 0.0, 1000.0, objects[:i2]), # this would trigger a spike much later
    ]
    logger=simulate!(net, input; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->tp==:spikes))
    @test isempty(logger.data)
end

@testset "Prevent multiple triggers" begin
    config = """
    refractory_duration: 100.01
    inputs:
    - id: i1
    - id: i2
    outputs:
    - id: o
    neurons:
    - id: n
      θ_syn: 1
      θ_seg: 2
      branches:
        - id: seg1
        - id: seg2
    synapses:
    - {id: syn1, source: i1, target: seg1}
    - {id: syn2, source: i2, target: seg2}
    - {id: syn3, source: i1, target: n}
    - {id: syn4, source: i2, target: n}
    - {id: syn_out, source: n, target: o, spike_duration: 0, delay: 0}
    """

    (net,objects) = load_network(YAML_source=config)

    input=[
        Event(:input_spikes, 0.0,  0.0, objects[:i1]),
        Event(:input_spikes, 0.0,  2.0, objects[:i2]),
    ]

    logger=simulate!(net, input, 10.0; show_progress=false, logger! = Logger(net; filter=(t,tp,id,x)->id==:n))
    n = filter(x->(x.object==:n && x.event==:spikes), logger.data)[!, :t]
    @test length(logger.data.t)==1 && n ≈ [3.0]
end


text:
function value(p::InitialCondition{Float64})
    return p.value
end

function value(p::InitialCondition{PJ.ParameterRef})
    return PJ.value(p.value)
end


text:
"""
    AnalogInChannelCount(dev::DWFDevice)

Return the number of analog in channels.
"""
function AnalogInChannelCount(dev::DWFDevice)
    cChannel = Ref(Cint(0))
    FDwfAnalogInChannelCount(dev.hdwf[], cChannel)
    return cChannel[]
end

"""
    AnalogInFrequencyInfo(dev::DWFDevice)

Return the maximum sampling frequency of the analog in channels.
"""
function AnalogInFrequencyInfo(dev::DWFDevice)
    hzFreq = Ref(Cdouble(0))
    FDwfAnalogInFrequencyInfo(dev.hdwf[], C_NULL, hzFreq)
    return hzFreq[]
end


text:
using LegendGeSim

##
# import Pkg
# Pkg.add(url="https://github.com/legend-exp/LegendTestData.jl.git")
# Pkg.build("LegendTestData")
using LegendTestData
testdata_path = joinpath(LegendTestData.legend_test_data_path(), "data", "ldsim")

##

# PET (position-energy-time) information of the simulated energy depositions
pet_file = joinpath(testdata_path, "single-invcoax-th228-geant4.csv")
# detector
det_metadata = joinpath(testdata_path, "invcoax-metadata.json")
# simulation settings
sim_config_name = "configs/SSD_NoiseSim.json"
#sim_config_name = "configs/siggen_NoiseSim.json"

##
raw_table = LegendGeSim.simulate_raw(pet_file, det_metadata, sim_config_name)

##
using HDF5
using LegendHDF5IO

##
mkdir("output")
HDF5.h5open("output/my_name_raw.h5", "w") do f
    LegendHDF5IO.writedata(f, "raw", raw_table)
end

##
using Plots

##
plot_wf = plot(raw_table.waveform[1:5])
# png(plot_wf, "raw.png")

##
stp_table = LegendGeSim.pet_to_stp(pet_file, det_metadata)

##
mkdir("cache")
HDF5.h5open("cache/my_name_stp.h5", "w") do f
    LegendHDF5IO.writedata(f, "stp", stp_table)
end

##
pss_table, pss_truth = LegendGeSim.stp_to_pss("cache/my_name_stp.h5", det_metadata, sim_config_name)

##
pss_table, pss_truth = LegendGeSim.stp_to_pss(stp_table, det_metadata, sim_config_name)

##
using Plots

##
plot_pss = plot(pss_table.waveform[1:5])
# png(plot_pss, "plots_md/pss.png")

##
h5open("cache/my_name_pss.h5", "w") do f
    LegendHDF5IO.writedata(f, "pss/pss", pss_table)
    LegendHDF5IO.writedata(f, "pss/truth", pss_truth)
end

##
raw_table = LegendGeSim.pss_to_raw("cache/my_name_pss.h5", det_metadata, sim_config_name)

##
raw_table = LegendGeSim.pss_to_raw(pss_table, pss_truth, det_metadata, sim_config_name)

##
pss_table, pss_truth = LegendGeSim.pet_to_pss(pet_filename, det_metadata, sim_config_name)


text:
#InspectDR: Glyph definitions
#-------------------------------------------------------------------------------

#=TODO:
   -deprecate :diagcross?
   -add :asterisk?
=#

#==Constants
===============================================================================#
const SUPPORTED_GLYPHS = Symbol[
	:+, :cross, :x, :xcross, :diagcross,
	:o, :circle,
	:*,
	:square, :diamond,
	:utriangle, :dtriangle, :ltriangle, :rtriangle,
	:uarrow, :darrow, :larrow, :rarrow, #Currently implements triangles
	:star, :star4, :star5, :star6, :star7, :star8,
	:pentagon, :hexagon, :heptagon, :octagon,
	:hline, :vline,
]

#==Main types
===============================================================================#
abstract type Glyph end

struct GlyphCircle <: Glyph
	radius::DReal
end
struct GlyphPolyline <: Glyph
	#TODO: Should we just store deltas? would be more efficient!
	x::Vector{DReal}
	y::Vector{DReal}
	closepath::Bool
end
GlyphPolyline(x, y; closepath::Bool=true, scale::DReal=1.0) =
	GlyphPolyline(x*scale, y*scale, closepath)
struct GlyphLineSegments <: Glyph
	x1::Vector{DReal}
	x2::Vector{DReal}
	y1::Vector{DReal}
	y2::Vector{DReal}
end
#NOTE: defining x1::Vector so that we can add "scale" keyword without re-definition issues.
GlyphLineSegments(x1::Vector, x2, y1, y2; scale::DReal=1.0) =
	GlyphLineSegments(x1*scale, x2*scale, y1*scale, y2*scale)

#==Helper functions
===============================================================================#
isglyph(::Glyph) = true
isglyph(s::Symbol) = (s != :none) #Not necessarily a valid glyph id

#==Generator functions
===============================================================================#
function _stargenerator(n::Int, webscale::DReal; scale::DReal=1.0)
	x = Vector{DReal}(undef, 2*n)
	y = Vector{DReal}(undef, 2*n)
	ϕ = range(0, stop=2*pi, length=2*n+1)

	for i in 1:2:length(x)
		x[i] = sin(ϕ[i])
		y[i] = cos(ϕ[i])
		x[i+1] = webscale*sin(ϕ[i+1])
		y[i+1] = webscale*cos(ϕ[i+1])
	end

	return GlyphPolyline(x, y, scale=scale)
end

function _polygongenerator(n::Int; scale::DReal=1.0)
	x = Vector{DReal}(undef, n)
	y = Vector{DReal}(undef, n)
	ϕ = range(0, stop=2*pi, length=n+1)

	for i in 1:length(x)
		x[i] = sin(ϕ[i])
		y[i] = cos(ϕ[i])
	end

	return GlyphPolyline(x, y, scale=scale)
end

function ptmap_rotate(g::GlyphPolyline, Θ::Float64)
	cosΘ = cos(Θ); sinΘ = sin(Θ) #WANTCONST
	x = similar(g.x); y = similar(g.y)
	for i in 1:length(g.x)
		xi = g.x[i]; yi = g.y[i]
		x[i] = xi*cosΘ - yi*sinΘ
		y[i] = xi*sinΘ + yi*cosΘ
	end
	return GlyphPolyline(x, y)
end

#==Predefined glyphs (implicitly defines relative sizes)
===============================================================================#
const GLYPH_CIRCLE = GlyphCircle(0.5)
const GLYPH_SQUARE = GlyphPolyline(
	[-1, 1, 1,-1], #x
	[ 1, 1,-1,-1], #y
	scale = 0.5
)
const GLYPH_DIAMOND = GlyphPolyline(
	[0, 1, 0,-1], #x
	[1, 0,-1, 0], #y
	scale = 0.6
)
const GLYPH_UTRIANGLE = GlyphPolyline(
	[ 0, 1, -1], #x
	[ 1,-1, -1], #y
	scale = 0.5
)
const GLYPH_DTRIANGLE = GlyphPolyline(
	[ 0, 1, -1], #x
	[-1, 1,  1], #y
	scale = 0.5
)
const GLYPH_LTRIANGLE = GlyphPolyline(
	[-1, 1,  1], #x
	[ 0, 1, -1], #y
	scale = 0.5
)
const GLYPH_RTRIANGLE = GlyphPolyline(
	[1, -1, -1], #x
	[0,  1, -1], #y
	scale = 0.5
)
const GLYPH_CROSS = GlyphPolyline(
	[-1, +1, +1, +4, +4, +1, +1, -1, -1, -4, -4, -1], #x
	[+4, +4, +1, +1, -1, -1, -4, -4, -1, -1, +1, +1], #y
	scale = 1/8
)
const GLYPH_LINECROSS = GlyphLineSegments(
	[-1, 0], [ 1, 0], #x1, x2
	[ 0,-1], [ 0, 1], #y1, y2
	scale = sqrt(2)/2
)
const GLYPH_XCROSS = ptmap_rotate(GLYPH_CROSS, pi/4)
const GLYPH_LINEXCROSS = GlyphLineSegments(
	[-1,-1], [ 1, 1], #x1, x2
	[ 1,-1], [-1, 1], #y1, y2
	scale = 0.5
)
const GLYPH_ASTERISK = GlyphLineSegments(
	vcat(GLYPH_LINECROSS.x1, GLYPH_LINEXCROSS.x1), vcat(GLYPH_LINECROSS.x2, GLYPH_LINEXCROSS.x2),
	vcat(GLYPH_LINECROSS.y1, GLYPH_LINEXCROSS.y1), vcat(GLYPH_LINECROSS.y2, GLYPH_LINEXCROSS.y2),
)
const GLYPH_HLINE = GlyphLineSegments(
	[-1], [ 1], #x1, x2
	[ 0], [ 0], #y1, y2
	scale = sqrt(2)/2
)
const GLYPH_VLINE = GlyphLineSegments(
	[ 0], [0], #x1, x2
	[-1], [1], #y1, y2
	scale = sqrt(2)/2
)
const GLYPH_STAR4 = _stargenerator(4, .5, scale=0.7)
const GLYPH_STAR5 = _stargenerator(5, .5, scale=0.7)
const GLYPH_STAR6 = _stargenerator(6, .5, scale=0.7)
const GLYPH_STAR7 = _stargenerator(7, .5, scale=0.7)
const GLYPH_STAR8 = _stargenerator(8, .5, scale=0.7)
const GLYPH_PENTAGON = _polygongenerator(5, scale=0.6)
const GLYPH_HEXAGON = _polygongenerator(6, scale=0.6)
const GLYPH_HEPTAGON = _polygongenerator(7, scale=0.6)
const GLYPH_OCTAGON = _polygongenerator(8, scale=0.6)

#==Scaling functions
===============================================================================#
#TODO: Could scale once for all points (more efficient)
#_scale(p::GlyphPolyline, _size::DReal) =
#	GlyphPolyline(p.x*_size, p.y*_size, p.closepath)

#==Accessors
===============================================================================#
Glyph(g::Glyph) = g #Passthrough.
function Glyph(s::Symbol)
	#TODO: use Map instead?  precompile issues?
	s == :o && return GLYPH_CIRCLE
	s == :+ && return GLYPH_LINECROSS
	s == :x && return GLYPH_LINEXCROSS
	s == :* && return GLYPH_ASTERISK

	s == :circle && return GLYPH_CIRCLE
	s == :square && return GLYPH_SQUARE
	s == :diamond && return GLYPH_DIAMOND

	s == :utriangle && return GLYPH_UTRIANGLE
	s == :dtriangle && return GLYPH_DTRIANGLE
	s == :ltriangle && return GLYPH_LTRIANGLE
	s == :rtriangle && return GLYPH_RTRIANGLE
	s == :uarrow && return GLYPH_UTRIANGLE
	s == :darrow && return GLYPH_DTRIANGLE
	s == :larrow && return GLYPH_LTRIANGLE
	s == :rarrow && return GLYPH_RTRIANGLE

	s == :cross && return GLYPH_CROSS
	s == :xcross && return GLYPH_XCROSS
	s == :diagcross && return GLYPH_XCROSS

	s == :star && return GLYPH_STAR5
	s == :star4 && return GLYPH_STAR4
	s == :star5 && return GLYPH_STAR5
	s == :star6 && return GLYPH_STAR6
	s == :star7 && return GLYPH_STAR7
	s == :star8 && return GLYPH_STAR8

	s == :pentagon && return GLYPH_PENTAGON
	s == :hexagon && return GLYPH_HEXAGON
	s == :heptagon && return GLYPH_HEPTAGON
	s == :octagon && return GLYPH_OCTAGON

	s == :hline && return GLYPH_HLINE
	s == :vline && return GLYPH_VLINE

	if isglyph(s) #Glyph requested, but not supported
		@warn("Glyph shape not supported: $s")
	end
	return nothing
end

#Last line


text:
abstract type SaturationGoal end

reached(g::EGraph, goal::Nothing) = false
reached(g::EGraph, goal::SaturationGoal) = false

"""
This goal is reached when the `exprs` list of expressions are in the
same equivalence class.
"""
struct EqualityGoal <: SaturationGoal
    exprs::Vector{Any}
    ids::Vector{EClassId}
    function EqualityGoal(exprs, eclasses)
        @assert length(exprs) == length(eclasses) && length(exprs) != 0
        new(exprs, eclasses)
    end
end

function reached(g::EGraph, goal::EqualityGoal)
    all(x -> in_same_class(g, goal.ids[1], x), @view goal.ids[2:end])
end

"""
Boolean valued function as an arbitrary saturation goal.
User supplied function must take an [`EGraph`](@ref) as the only parameter.
"""
struct FunctionGoal <: SaturationGoal
    fun::Function
end

function reached(g::EGraph, goal::FunctionGoal)::Bool
    fun(g)
end


text:
# The code here is mostly copied from the flip_sign.jl code for FlipSignBridge and GreaterToLessBridge

"""
    AbstractToIntervalBridge{T, S1, F}

Bridge a `F`-in-`Interval` constraint into an `F`-in-`Interval{T}` constraint where we have either:
* `S1 = MOI.GreaterThan{T}`
* `S1 = MOI.LessThan{T}`

The `F`-in-`Interval{T}` constraint is stored in the `constraint`
field by convention.
It is required that T be a AbstractFloat type because otherwise
typemin and typemax would either be not implemented (e.g. BigInt)
or would not give infinite value (e.g. Int).
"""
abstract type AbstractToIntervalBridge{
    T<:AbstractFloat,
    S1<:MOI.AbstractScalarSet,
    F<:MOI.AbstractScalarFunction,
} <: SetMapBridge{T,MOI.Interval{T},S1,F,F} end

# The function map is the identity. It is also an involution, symmetric, and a symmetric involution.
map_function(::Type{<:AbstractToIntervalBridge{T}}, func) where {T} = func
inverse_map_function(::Type{<:AbstractToIntervalBridge}, func) = func
adjoint_map_function(::Type{<:AbstractToIntervalBridge}, func) = func
inverse_adjoint_map_function(::Type{<:AbstractToIntervalBridge}, func) = func

# FIXME are these modify functions necessary?
function MOI.modify(
    model::MOI.ModelLike,
    bridge::AbstractToIntervalBridge,
    change::MOI.ScalarCoefficientChange,
)
    return MOI.modify(model, bridge.constraint, change)
end
function MOI.modify(
    model::MOI.ModelLike,
    bridge::AbstractToIntervalBridge,
    change::MOI.MultirowChange{T},
) where {T}
    return MOI.modify(model, bridge.constraint, change)
end

"""
    GreaterToIntervalBridge{T, F<:MOI.AbstractScalarFunction} <:
        AbstractToIntervalBridge{T, MOI.GreaterThan{T}, F}

Transforms a `F`-in-`GreaterThan{T}` constraint into an `F`-in-`Interval{T}`
constraint.
"""
struct GreaterToIntervalBridge{T,F<:MOI.AbstractScalarFunction} <:
       AbstractToIntervalBridge{T,MOI.GreaterThan{T},F}
    constraint::CI{F,MOI.Interval{T}}
end
function map_set(::Type{<:GreaterToIntervalBridge}, set::MOI.GreaterThan)
    return MOI.Interval(set.lower, typemax(set.lower))
end
function inverse_map_set(::Type{<:GreaterToIntervalBridge}, set::MOI.Interval)
    return MOI.GreaterThan(set.lower)
end
function concrete_bridge_type(
    ::Type{<:GreaterToIntervalBridge{T}},
    F::Type{<:MOI.AbstractScalarFunction},
    ::Type{MOI.GreaterThan{T}},
) where {T}
    return GreaterToIntervalBridge{T,F}
end

"""
    LessToIntervalBridge{T, F<:MOI.AbstractScalarFunction} <:
        AbstractToIntervalBridge{T, MOI.LessThan{T}, F}

Transforms a `F`-in-`LessThan{T}` constraint into an `F`-in-`Interval{T}`
constraint.
"""
struct LessToIntervalBridge{T,F<:MOI.AbstractScalarFunction} <:
       AbstractToIntervalBridge{T,MOI.LessThan{T},F}
    constraint::CI{F,MOI.Interval{T}}
end
function map_set(::Type{<:LessToIntervalBridge}, set::MOI.LessThan)
    return MOI.Interval(typemin(set.upper), set.upper)
end
function inverse_map_set(::Type{<:LessToIntervalBridge}, set::MOI.Interval)
    return MOI.LessThan(set.upper)
end
function concrete_bridge_type(
    ::Type{<:LessToIntervalBridge{T}},
    F::Type{<:MOI.AbstractScalarFunction},
    ::Type{MOI.LessThan{T}},
) where {T}
    return LessToIntervalBridge{T,F}
end


text:
import PowerModels
import JSON
import PyCall

function _load_case(casefile::String)
    data = PowerModels.parse_file(casefile)
    net = JSON.json(data)
    return net
end

return load_case = PyCall.pyfunctionret(_load_case, String, String)


text:
using ImageMagick

info("ImageMagick version ", ImageMagick.libversion())

include("constructed_images.jl")
include("readremote.jl")
include("badimages.jl")

workdir = joinpath(tempdir(), "Images")
try
    rm(workdir, recursive=true)
end


text:
using EDF
using EDF: TimestampedAnnotationList, PatientID, RecordingID, SignalHeader,
           Signal, AnnotationsSignal
using Dates
using Test

#####
##### Testing utilities
#####

function deep_equal(a::T, b::T) where T
    nfields = fieldcount(T)
    if nfields == 0
        return isequal(a, b)  # Use `isequal` instead of `==` to handle `missing`
    else
        for i = 1:nfields
            typeof(getfield(a, i)) <: IO && continue # Two different files will have different IO sources
            isdefined(a, i) || return !isdefined(b, i)  # Call two undefs equal
            deep_equal(getfield(a, i), getfield(b, i)) || return false
        end
    end
    return true
end

function deep_equal(a::T, b::T) where T<:AbstractArray
    length(a) == length(b) || return false
    for (x, y) in zip(a, b)
        deep_equal(x, y) || return false
    end
    return true
end

deep_equal(a::T, b::S) where {T,S} = false

#####
##### Actual tests
#####

const DATADIR = joinpath(@__DIR__, "data")

@testset "Just Do It" begin
    # test EDF.read(::AbstractString)
    edf = EDF.read(joinpath(DATADIR, "test.edf"))
    @test sprint(show, edf) == "EDF.File with 140 16-bit-encoded signals"
    @test edf.header.version == "0"
    @test edf.header.patient == PatientID(missing, missing, missing, missing)
    @test edf.header.recording == RecordingID(Date(2014, 4, 29), missing, missing, missing)
    @test edf.header.is_contiguous
    @test edf.header.start == DateTime(2014, 4, 29, 22, 19, 44)
    @test edf.header.record_count == 6
    @test edf.header.seconds_per_record == 1.0
    @test edf.signals isa Vector{Union{Signal{Int16},AnnotationsSignal}}
    @test length(edf.signals) == 140
    for signal in edf.signals
        if signal isa EDF.Signal
            @test length(signal.samples) == signal.header.samples_per_record * edf.header.record_count
        else
            @test length(signal.records) == edf.header.record_count
            # XXX seems like this test file actually contains nonsensical onset timestamps...
            # according to the EDF+ specification, onsets should be relative to the start time of
            # the entire file, but it seems like whoever wrote these onsets might have used values
            # that were relative to the start of the surrounding data record
            expected = [[TimestampedAnnotationList(0.0, nothing, String[""]), TimestampedAnnotationList(0.0, nothing, ["start"])],
                        [TimestampedAnnotationList(1.0, nothing, String[""]), TimestampedAnnotationList(0.1344, 0.256, ["type A"])],
                        [TimestampedAnnotationList(2.0, nothing, String[""]), TimestampedAnnotationList(0.3904, 1.0, ["type A"])],
                        [TimestampedAnnotationList(3.0, nothing, String[""]), TimestampedAnnotationList(2.0, nothing, ["type B"])],
                        [TimestampedAnnotationList(4.0, nothing, String[""]), TimestampedAnnotationList(2.5, 2.5, ["type A"])],
                        [TimestampedAnnotationList(5.0, nothing, String[""])]]
            @test all(signal.records .== expected)
            @test AnnotationsSignal(signal.records).samples_per_record == 16
        end
    end

    # test EDF.write(::IO, ::EDF.File)
    io = IOBuffer()
    EDF.write(io, edf)
    seekstart(io)
    file = EDF.File(io)
    @test deep_equal(edf.header, file.header)
    @test all(isempty(s isa Signal ? s.samples : s.records) for s in file.signals)
    EDF.read!(file)
    @test deep_equal(edf.signals, file.signals)
    @test eof(io)
    # ensure that multiple `EDF.read!` calls don't error and have no effect by
    # simply rerunning the exact same test as above
    EDF.read!(file)
    @test deep_equal(edf.signals, file.signals)

    # test that EDF.write(::IO, ::EDF.File) errors if file is
    # discontiguous w/o an AnnotationsSignal present
    bad_file = EDF.File(IOBuffer(),
                        EDF.FileHeader(file.header.version,
                                       file.header.patient,
                                       file.header.recording,
                                       file.header.start,
                                       false, # is_contiguous
                                       file.header.record_count,
                                       file.header.seconds_per_record),
                        filter(s -> !(s isa AnnotationsSignal), file.signals))
    @test_throws ArgumentError EDF.write(IOBuffer(), bad_file)

    # test EDF.write(::AbstractString, ::EDF.File)
    mktempdir() do dir
        path = joinpath(dir, "tmp.edf")
        EDF.write(path, edf)
        file = EDF.File(open(path, "r"))
        @test deep_equal(edf.header, file.header)
        @test all(isempty(s isa Signal ? s.samples : s.records) for s in file.signals)
        EDF.read!(file)
        @test deep_equal(edf.signals, file.signals)
        @test eof(io)
    end

    @test EDF._edf_repr(EDF._nearest_representable_edf_time_value(-0.0023405432)) == "-0.00234"
    @test EDF._edf_repr(EDF._nearest_representable_edf_time_value(0.0023405432)) == "0.002340"
    @test EDF._edf_repr(EDF._nearest_representable_edf_time_value(1.002343)) == "1.002343"
    @test EDF._edf_repr(EDF._nearest_representable_edf_time_value(1011.05432)) == "1011.054"
    @test EDF._edf_repr(EDF._nearest_representable_edf_time_value(-1011.05432)) == "-1011.05"
    @test EDF._edf_repr(EDF._nearest_representable_edf_time_value(-1013441.5)) == "-1013442"
    @test EDF._edf_repr(EDF._nearest_representable_edf_time_value(-1013441.3)) == "-1013441"
    @test EDF._edf_repr(34577777) == "34577777"
    @test EDF._edf_repr(0.0345) == "0.034500"
    @test EDF._edf_repr(-0.02) == "-0.02000"
    @test EDF._edf_repr(-187.74445) == "-187.744"
    @test_throws ErrorException EDF._edf_repr(123456789)
    @test_throws ErrorException EDF._edf_repr(-12345678)
    @test_throws ErrorException EDF._edf_repr(0.00000000024)
    @test_throws ErrorException EDF.edf_write(IOBuffer(), "hahahahaha", 4)

    uneven = EDF.read(joinpath(DATADIR, "test_uneven_samp.edf"))
    @test sprint(show, uneven) == "EDF.File with 2 16-bit-encoded signals"
    @test uneven.header.version == "0"
    @test uneven.header.patient == "A 3Hz sinewave and a 0.2Hz block signal, both starting in their positive phase"
    @test uneven.header.recording == "110 seconds from 13-JUL-2000 12.05.48hr."
    @test uneven.header.is_contiguous
    @test uneven.header.start == DateTime(2000, 1, 31, 23, 0, 59)
    @test uneven.header.record_count == 11
    @test uneven.header.seconds_per_record == 10.0
    @test uneven.signals[1].header.samples_per_record != uneven.signals[2].header.samples_per_record
    @test length(uneven.signals) == 2

    nonint = EDF.read(joinpath(DATADIR, "test_float_extrema.edf"))
    signal = nonint.signals[1]
    @test signal.header.physical_minimum ≈ -29483.1f0
    @test signal.header.physical_maximum ≈ 29483.12f0
    @test signal.header.digital_minimum ≈ -32767.0f0
    @test signal.header.digital_maximum ≈ 32767.0f0

    # Python code for generating the comparison values used here:
    # ```
    # import mne
    # edf = mne.io.read_raw_edf("test/data/test_float_extrema.edf")
    # signal = edf.get_data()[0] * 1e6  # The 1e6 converts volts back to microvolts
    # with open("test/data/mne_values.csv", "w") as f:
    #     for x in signal:
    #         f.write("%s\n" % x)
    # ```
    mne = map(line->parse(Float32, line), eachline(joinpath(DATADIR, "mne_values.csv")))
    for (a, b) in zip(EDF.decode(signal), mne)
        @test a ≈ b atol=0.01
    end

    # Truncated files
    mktempdir() do dir
        # note that this tests a truncated final record, not an incorrect number of records
        full_edf_bytes = read(joinpath(DATADIR, "test.edf"))
        write(joinpath(dir, "test_truncated.edf"), full_edf_bytes[1:(end - 1)])
        @test_logs((:warn, "Number of data records in file header does not match " *
                    "file size. Skipping 1 truncated data record(s)."),
                   EDF.read(joinpath(dir, "test_truncated.edf")))
        truncated_edf = EDF.read(joinpath(dir, "test_truncated.edf"))
        for field in fieldnames(EDF.FileHeader)
            a = getfield(edf.header, field)
            b = getfield(truncated_edf.header, field)
            if field === :record_count
                @test b == a - 1
            else
                @test a == b
            end
        end
        for i in 1:length(edf.signals)
            good = edf.signals[i]
            bad = truncated_edf.signals[i]
            if good isa EDF.Signal
                @test deep_equal(good.header, bad.header)
                @test good.samples[1:(end - good.header.samples_per_record)] == bad.samples
            else
                @test good.samples_per_record == bad.samples_per_record
            end
        end
        @test deep_equal(edf.signals[end].records[1:(edf.header.record_count - 1)],
                         truncated_edf.signals[end].records)
    end

    @test EDF._size(IOBuffer("memes")) == 5
    @test EDF._size(Base.DevNull()) == -1

    @testset "BDF Files" begin
        # The corresponding EDF file was exported by 3rd party software based on the BDF,
        # so some differences are inevitable, but we just want to check that the values
        # are Close Enough™.
        bdf = EDF.read(joinpath(DATADIR, "bdf_test.bdf"))
        comp = EDF.read(joinpath(DATADIR, "bdf_test.edf"))
        for i in 1:8
            bdf_values = EDF.decode(bdf.signals[i])
            comp_values = EDF.decode(comp.signals[i])
            @test bdf_values ≈ comp_values rtol=0.01
        end
        # Ensure that BDF files can also be round-tripped
        mktempdir() do dir
            path = joinpath(dir, "tmp.bdf")
            EDF.write(path, bdf)
            file = EDF.read(path)
            @test deep_equal(bdf, file)
        end
        @test EDF.sample_type(bdf) == EDF.Int24
        @test EDF.sample_type(comp) == Int16
        @test EDF.is_bdf(bdf)
        @test !EDF.is_bdf(comp)
        @test sprint(show, bdf) == "EDF.File with 8 24-bit-encoded signals"
    end
end