text:
using Eirene
using Pipe

#%%
# ==============================
#  ======== Tested code ========

#%%
# ================================
#  ======== Untested code ========

"""
    get_barcodes(results_eirene::Dict, max_dim::Integer; min_dim::Int=1)

Calls Eirene.barcode for 'dim' in range from `min_dim` up to 'max_dim' and
stack the resulting Arrays into a vector.

The returned value is a Vector of Arrays{Float64,2}. Each array is of
different size, because of different number of detected cycles. First column
of each array contains birth step, second column contains death step.

Arrays in returned vector correspond to Betti curve dimensions form range
`min_dim` up to 'max_dim'.

"""
function get_barcodes(results_eirene::Dict, max_dim::Integer; min_dim::Int = 1, sorted::Bool=false)
    barcodes = Matrix{Float64}[]
    for d = min_dim:max_dim
        result = barcode(results_eirene, dim = d)
        if isempty(result) && d > 1
            result = zeros(size(barcodes[d-1]))
        end
        if sorted
            result = result[sortperm(result[:,1]),:]
        end
        push!(barcodes, result)
    end
    return barcodes
end

"""
  plot_barcodes(barcodes; min_dim::Integer=1, betti_labels::Bool=true, default_labels::Bool=true kwargs...)

Creates a plot for set of barcodesstored in `barcodes` and return the
handler to the plot.

'kwargs' are plot parameters

Some of the possible 'kwargs' are:
  - title::String
  - legend:Bool
  - size::Tuple{T, T} where {T::Number}
  - lw::Integer or linewidth:Integer
(for more, see plots documentation):
"""
function plot_barcodes!(barcodes::Vector, plot_ref;
                        min_dim::Integer = 1,
                        barcodes_labels::Bool = true,
                        default_labels::Bool = true,
                        sort_by_birth::Bool = false,
                        alpha::Float64 = 1.,
                        kwargs...)#; plot_size = (width=1200, height=800),
    # TODO add change of x label based on x values- so it is either edge density for 0:1 range values or Filtration step otherwise
    # TODO add ordering of bars to firstly birth time, then death time
    # TODO dims should be given as range, not a single min dim

    # barcodes = all_barcodes_geom
    max_dim = size(barcodes, 1) - (1-min_dim) # TODO not sure if this is correct solution

    if min_dim > max_dim
        throw(DomainError(
            min_dim,
            "\'min_dim\' must be greater that maximal dimension in \'bettis\'",
        ))
    end

    lw_pos = findfirst(x -> x == :lw || x == :linewidth, keys(kwargs))
    if !isnothing(lw_pos)
        lw = kwargs[lw_pos]
    else
        lw = 8
    end

    colors_set = get_bettis_color_palete(min_dim=min_dim)

    dims_indices = 1:length(min_dim:max_dim)
    all_sizes = [size(barcodes[k],1) for k = dims_indices]
    ranges_sums = vcat(0,[sum(all_sizes[1:k]) for k = dims_indices])
    y_val_ranges = [ranges_sums[k]+1:ranges_sums[k+1] for k = dims_indices]

    if sort_by_birth
        sort_barcodes!(barcodes, min_dim, max_dim)
    end

    total_cycles = sum([size(barcodes[k],1) for (k, dim) in enumerate(min_dim:max_dim)])
    # for p = min_dim:(max_dim) #TODO ths can not be starting from min_dim, because it may be 0
    for (p, dim) = enumerate(min_dim:max_dim)# 1:(max_dim) #TODO ths can not be starting from min_dim, because it may be 0
        # @info p, dim
        args = (lc = colors_set[p], linewidth = lw)

        # b = barcodes[p][1:1:(end-1),:]
        b = barcodes[p][:,:]
        all_infs = findall(x->isinf(x), b)
        for k in all_infs
            b[k] = 1
        end

        if dim==0
            b = sort(b, dims=1)
        end

        total_bars = size(b,1)
        y_vals = [[k, k] for k in y_val_ranges[p]]
        lc = colors_set[p]
        for k = 1:total_bars
            # TODO change label to empty one
            plot!(b[k,:], y_vals[k]; label="", lc=lc, alpha=alpha)#; args...)
        end
        if false && betti_labels
            label = "β$(dim)"
        end
        # plot!(label=label)
    end

    # display(plot_ref)

    legend_pos = findfirst(x -> x == :legend, keys(kwargs))
    if !isnothing(legend_pos)
        plot!(legend = kwargs[legend_pos])
    else
        all_labels = reshape(["β$(k)" for k in 1:max_dim], (1,max_dim))
        plot!(label = all_labels)
        plot!(legend = true)
    end

    x_pos = findfirst(x -> x == :xlabel, keys(kwargs))
    y_pos = findfirst(x -> x == :ylabel, keys(kwargs))
    if !isnothing(x_pos)
        xlabel!(kwargs[x_pos])
    elseif default_labels
        xlabel!("Birth/Death")
    end
    if !isnothing(y_pos)
        ylabel!(kwargs[y_pos])
    elseif default_labels
        ylabel!("Cycle")
    end

    # ylims!(0, 2*total_cycles)
    ylims!(0, total_cycles+2)

    return plot_ref
end

function plot_barcodes(barcodes::Vector;
                        kwargs...)
    plot_ref = plot(; kwargs...)
    plot_barcodes!(barcodes, plot_ref;
                        kwargs...)

end

function sort_barcodes!(barcodes, min_dim, max_dim)
    for (dim_index, dim) = enumerate(min_dim:max_dim)
        if dim==0
            sort!(barcodes[dim_index], dims=1)
        else
            sort!(barcodes[dim_index], dims=1)
        end
    end
end
# TODO This has to be imported from other file
# function get_bettis_color_palete(; min_dim = 1, use_set::Integer = 1)
#     """
#     	function get_bettis_color_palete()
#
#     Generates vector with colours used for Betti plots. Designed for Betti plots consistency.
#     """
#     # TODO what does the number in the function below is used for?
#
#     if use_set == 1
#         cur_colors = [Gray(bw) for bw = 0.0:0.025:0.5]
#         if min_dim == 0
#             colors_set = [RGB(87 / 256, 158 / 256, 0 / 256)]
#         else
#             colors_set = RGB[]
#         end
#         colors_set = vcat(
#             colors_set,
#             [
#                 RGB(255 / 256, 206 / 256, 0 / 256),
#                 RGB(248 / 256, 23 / 256, 0 / 256),
#                 RGB(97 / 256, 169 / 256, 255 / 256),
#                 RGB(163 / 256, 0 / 256, 185 / 256),
#                 RGB(33 / 256, 96 / 256, 45 / 256),
#                 RGB(4 / 256, 0 / 256, 199 / 256),
#                 RGB(135 / 256, 88 / 256, 0 / 256),
#             ],
#             cur_colors,
#         )
#     else
#         use_set == 2
#         cur_colors = get_color_palette(:auto, 1)
#         cur_colors3 = get_color_palette(:lightrainbow, 1)
#         cur_colors2 = get_color_palette(:cyclic1, 1)
#         if min_dim == 0
#             # colors_set =  [cur_colors[3], cur_colors[5], [:red], cur_colors[1]] #cur_colors[7],
#             colors_set = [cur_colors3[3], cur_colors[5], cur_colors3[end], cur_colors[1]] #cur_colors[7],
#         else
#             colors_set = [cur_colors[5], cur_colors3[end], cur_colors[1]] #cur_colors[7],
#             # colors_set =  [cur_colors[5], [:red], cur_colors[1], cur_colors[14]]
#         end
#         # for c =  [collect(11:25);]
#         #     push!(colors_set, cur_colors2[c])
#         # end
#         colors_set = vcat(colors_set, [cur_colors2[c] for c in [collect(11:25);]])
#     end
#
#     return colors_set
# end

function get_birth_death_ratio(barcodes; max_dim::Integer=3)
    # birth_death_raio_π = [[all_barcodes_geom[k][m,2]/all_barcodes_geom[k][m,1] for m= 1:size(all_barcodes_geom[k],1)] for k in 1:max_dim]
    birth_death_ratio_π = [barcodes[k][:,2]./barcodes[k][:,1] for k in 1:max_dim]
    return birth_death_ratio_π
end

function get_barcode_lifetime(barcodes; max_dim::Integer=3)
    lifetime = [barcodes[k][:,2].-barcodes[k][:,1] for k in 1:max_dim]
    return lifetime
end

#%%
"""
    get_barcode_max_lifetime(lifetimes, min_dim, max_dim)

Returns the maximal life times of barcode for all dimensions.
"""
function get_barcode_max_lifetime(lifetimes)
    total_lifetimes = length(lifetimes)
    all_max_lifetimes = zeros(total_lifetimes, 1)

    for k in 1:total_lifetimes
        all_max_lifetimes[k] = findmax(lifetimes[k])[1]
    end

    return all_max_lifetimes
end

"""
    boxplot_birth_death(areas_matrix, min_dim::Integer, max_dim::Integer)

Plots the boxplot of area under betti curves.
"""
function boxplot_birth_death(birth_death_ratio_π, min_dim::Integer, max_dim::Integer)
    bplot = StatsPlots.boxplot()

    data_colors = get_bettis_color_palete()
    total_plots = size(birth_death_ratio_π,1)

    for k in 1:total_plots
        StatsPlots.boxplot!(bplot, birth_death_ratio_π[k], labels="β$(k)", color=data_colors[k])
        StatsPlots.dotplot!(bplot, birth_death_ratio_π[k],color=data_colors[k])
    end

    return bplot
end

"""
    boxplot_birth_death(areas_matrix, min_dim::Integer, max_dim::Integer)

Plots the boxplot of area under betti curves.
"""
function boxplot_lifetime(barcode_lifetime, min_dim::Integer, max_dim::Integer)
    bplot = StatsPlots.boxplot()

    data_colors = get_bettis_color_palete()
    total_plots = size(barcode_lifetime,1)

    for k in 1:total_plots
        StatsPlots.boxplot!(bplot, barcode_lifetime[k], labels="β$(k)", color=data_colors[k])
        StatsPlots.dotplot!(bplot, barcode_lifetime[k],color=data_colors[k])
    end

    return bplot
end

#%%
"""
    get_barcode_max_db_ratios(lifetimes, min_dim, max_dim)

Returns the maximal life times of barcode for all dimensions.
"""
function get_barcode_max_db_ratios(db_ratos)
    total_db = length(db_ratos)
    all_max_db = zeros(total_db, 1)

    for k in 1:total_db
        all_max_db[k] = findmax(db_ratos[k])[1]
    end

    return all_max_db
end

"""
    get_normalised_barcodes(barcodes::Vector, betti_numbers::Array)

Returns the barcodes which values are within [0,1] range.
1. get corresponding bettis
2. get max number of steps
3. divide all values by total number of steps.

"""
function get_normalised_barcodes(barcodes, betti_numbers::Array)
    if typeof(betti_numbers) == Vector
        total_steps = size(betti_numbers[1],1)
    else
        total_steps = size(betti_numbers,1)
    end

    return barcodes ./ total_steps
end

"""
    get_normalised_barcodes_collection(barcodes_collection, bettis_collection)

Applies get_normalised_barcodes to the collection of barcodes and corresponding
betti curves.
"""
function get_normalised_barcodes_collection(barcodes_collection, bettis_collection)
    if size(barcodes_collection,1) != size(bettis_collection,1)
        throw(BoundsError(barcodes_collection, bettis_collection,
            "Both collections must have same number of elements",
        ))
    else
        total_collections  = size(barcodes_collection,1)
    end
    normed_collection = deepcopy(barcodes_collection)

    for k = 1:total_collections
        normed_collection[k]=get_normalised_barcodes(barcodes_collection[k], bettis_collection[k])
    end
    return normed_collection
end

"""
  plot_bd_diagram(barcodes;
                        dims::Range,
                        use_js::Bool=false,
                        kwargs...)

Creates a birth/death diagram from `barcodes` and returns the handlers to
the plots.

By default, dims is set to range '1:length(barcodes)', which plots all of
the diagrams. If set to an integer, plots only 1 dimension.

If 'use_js' is set to true, plotly backend is used for plotting.

'kwargs' are plot parameters

Some of the possible 'kwargs' are:
  - title::String
  - legend:Bool
  - size::Tuple{T, T} where {T::Number}
  - lw::Integer or linewidth:Integer
(for more, see plots documentation):
TODO dims are defined as if the dimensions always starts at 1- this has to be changed
"""
function plot_bd_diagram(barcodes; dims=1:length(barcodes), use_js::Bool=false,
                                                class_sizes=[],
                                                class_labels=[],
                                                kwargs...)
    # TODO max min should be ready to use from input data- might be better to have better structures as an inupt
    # max_dim = size(barcodes, 1)
    # min_dim = findmin(dims)[1]
    min_dim = dims[1]
    max_dim = dims[end]

    if max_dim > length(barcodes)
        throw(DimensionMismatch("Can not have dims range larger than barcodes length"))
    end
    # all_dims = min_dim:max_dim

    # if findmax(dims)[1] > max_dim
    #     throw(DomainError(
    #         min_dim,
    #         "\'dims\' must be less than maximal dimension in \'bettis\'",
    #     ))
    # end

    # lw_pos = findfirst(x -> x == :lw || x == :linewidth, keys(kwargs))
    # if !isnothing(lw_pos)
    #     lw = kwargs[lw_pos]
    # else
    #     lw = 2
    # end

    colors_set = TopologyPreprocessing.get_bettis_color_palete(min_dim=min_dim)

    if use_js
        plotly()
    else
        gr()
    end

    plot_ref = plot(;xlims=(0,1), ylims=(0,1), kwargs...)

    for (p,dim) in enumerate(dims)
        # colors_set[p]
        my_vec = barcodes[p]

        # TODO class size is not a default and basic bd diagram property- should be factored out to other function
        if class_labels != [] && class_sizes != []
            labels = ["class/size $(class_labels[k])/$(class_sizes[class_labels[k]])" for k in 1:size(class_labels,1)]
        elseif  class_sizes == []
            labels = ["class: $(k)" for k in 1:size(my_vec,1)]
        else
            labels = ["class/size $(k)/$(class_sizes[k])" for k in 1:size(my_vec,1)]
        end

        args = (color = colors_set[p],
                # linewidth = lw,
                label="β$(dim)",
                aspect_ratio=1,
                size = (600,600),
                legend = :bottomright,
                # hover = labels,
                kwargs...)
        if class_labels != [] class_sizes != []
            for x = class_labels
                plot!(my_vec[Int(x), 1], my_vec[Int(x), 2], seriestype = :scatter; args...)
            end
        end

        plot!(my_vec[:, 1], my_vec[:, 2], seriestype = :scatter; args...)
    end

    # Add diagonal
    max_x = 1#  @pipe vcat([barcodes[d][:,1] for d in dims]...) |> findmax(_)[1]
    max_y = 1# @pipe vcat([barcodes[d][:,2] for d in dims]...) |> findmax(_)[1]
    plot!([0, max_x], [0, max_y], label = "")
    # xlims!(0,1)
    # ylims!(0,1)

    return plot_ref
end

#
"""
  plot_all_bd_diagrams(barcodes_collection;
                        min_dim::Integer=1,
                        betti_labels::Bool=true,
                        default_labels::Bool=true,
                        all_legend=false,
                        my_alpha=0.12,
                        aspect_ratio=1,
                        kwargs...)

Creates a set of birth/death diagrams from `barcodes_collection`
     and returns a dictionary with the handlers to the plots.

'kwargs' are plot parameters

Some of the possible 'kwargs' are:
  - title::String
  - legend:Bool
  - size::Tuple{T, T} where {T::Number}
  - lw::Integer or linewidth:Integer
(for more, see plots documentation):
"""
function plot_all_bd_diagrams(barcodes_collection;
                                min_dim::Integer = 1,
                                betti_labels::Bool = true,
                                default_labels::Bool = true,
                                all_legend=false,
                                my_alpha=0.12,
                                aspect_ratio=1,
                                base_w=600,
                                base_h=600,
                                kwargs...)
    total_dims = size(barcodes_collection[1],1)

    lw_pos = findfirst(x -> x == :lw || x == :linewidth, keys(kwargs))
    if !isnothing(lw_pos)
        lw = kwargs[lw_pos]
    else
        lw = 2
    end

    title_pos = findfirst(x -> x == :title, keys(kwargs))
    if !isnothing(title_pos)
        my_title = kwargs[title_pos]
    else
        my_title = "Birth death diagram"
    end

    colors_set = TopologyPreprocessing.get_bettis_color_palete(min_dim=min_dim)
    plot_dict = Dict()

    for b = 1:total_dims
        args = (lc = colors_set[b],
                linewidth = lw,
                label=false,
                aspect_ratio=aspect_ratio,
                size = (base_w,base_h),
                kwargs...)
        plot_dict["β$(b)"] = scatter(;xlims=(0,1), ylims=(0,1), dpi=300, args...)
        for bars = barcodes_collection
            barcode = bars[b]

            scatter!(barcode[:, 1], barcode[:, 2],
                    markeralpha=my_alpha,
                    markercolor = colors_set[b],
                    dpi=300)
        end
        plot!(legend=all_legend)
        plot!(title=(my_title*", β$(b)"))

        # legend_pos = findfirst(x -> x == :legend, keys(kwargs))
        # if !isnothing(legend_pos)
        #     plot!(legend = kwargs[legend_pos])
        # else
        #     plot!(legend = betti_labels)
        # end

        x_pos = findfirst(x -> x == :xlabel, keys(kwargs))
        y_pos = findfirst(x -> x == :ylabel, keys(kwargs))
        if !isnothing(x_pos)
            xlabel!(kwargs[x_pos])
        elseif default_labels
            xlabel!("Birth")
        end
        if !isnothing(y_pos)
            ylabel!(kwargs[y_pos])
        elseif default_labels
            ylabel!("Death")
        end
    end

    return plot_dict
end

## ===-
# Simpler plotting

"""
  plot_bd_diagram(barcodes;
                        dims::Range,
                        use_js::Bool=false,
                        kwargs...)

Creates a birth/death diagram from `barcodes` and returns the handlers to
the plots.

By default, dims is set to range '1:length(barcodes)', which plots all of
the diagrams. If set to an integer, plots only 1 dimension.

If 'use_js' is set to true, plotly backend is used for plotting.

'kwargs' are plot parameters

Some of the possible 'kwargs' are:
  - title::String
  - legend:Bool
  - size::Tuple{T, T} where {T::Number}
  - lw::Integer or linewidth:Integer
(for more, see plots documentation):
TODO dims are defined as if the dimensions always starts at 1- this has to be changed
"""
function plot_simple_bd_diagram(barcodes; dims=1:length(barcodes), max_bd=0, use_js::Bool=false,
                                                kwargs...)
    # TODO max min should be ready to use from input data- might be better to have better structures as an inupt
    max_dim = size(barcodes, 1)
    min_dim = findmin(dims)[1]
    # all_dims = min_dim:max_dim

    if findmax(dims)[1] > max_dim
        throw(DomainError(
            min_dim,
            "\'dims\' must be less than maximal dimension in \'bettis\'",
        ))
    end

    lw_pos = findfirst(x -> x == :lw || x == :linewidth, keys(kwargs))
    if !isnothing(lw_pos)
        lw = kwargs[lw_pos]
    else
        lw = 2
    end

    colors_set = TopologyPreprocessing.get_bettis_color_palete(min_dim=1)

    if use_js
        plotly()
    else
        gr()
    end

    plot_ref = plot(;kwargs...)

    for p in dims
        # colors_set[p]
        my_vec = barcodes[p]

        args = (color = colors_set[p],
                linewidth = lw,
                aspect_ratio=1,
                size = (600,600),
                legend = :bottomright,
                kwargs...)

        # scatter!(my_vec[:, 1], my_vec[:, 2], args...)
        plot!(my_vec[:, 1], my_vec[:, 2], seriestype = :scatter; args...)
    end

    # Add diagonal

    if max_bd>0
        max_x = max_bd
        max_y = max_bd
        plot!([0, max_y], [0, max_y], label = "")
    else
        all_births = vcat([barcodes[d][:,1] for d in dims]...)
        all_deaths = vcat([barcodes[d][:,2] for d in dims]...)
        max_x = findmax(all_births)[1]
        max_y = findmax(all_deaths)[1]
        plot!([0, max_y], [0, max_y], label = "")
    end

    return plot_ref
end


text:
module Assets

export
    @asset_str,
    copy_assets

macro asset_str(relpath)
    dir_path = string(__source__.file)
    parsed_relpath = esc(Meta.parse("\"$(escape_string(relpath))\""))
    return quote
        joinpath(assetpath(getpkgfolder($dir_path)), $parsed_relpath)
    end
end

"""
Returns the full path to the project's assets.

If called from a regular Julia session it will just return
the supplied argument augmented with "/assets/"

However, if called from a session compiled using PackageCompiler
it will return the path to the assets subfolder in the packaged
app.
"""
function assetpath(pkgfolder::AbstractString)
    pkgfolder = normpath(pkgfolder)
    if is_standard_julia_session()
        return joinpath(pkgfolder, "assets")
    else
        modulename = splitpath(pkgfolder)[end]
        return normpath(joinpath(Sys.BINDIR, "..", "assets", modulename))
    end
end

function getpkgfolder(path)
    folders = splitpath(path)
    while folders[end] != "src"
        pop!(folders)
    end
    pop!(folders)
    return joinpath(folders...)
end

"""
Returns true if called from a standard Julia session.
False if called from an app created using PackageCompiler.
"""
function is_standard_julia_session()
    julia_folder = splitpath(Sys.BINDIR)[end-1]
    return match(r"Julia.\d+\.\d+\.\d+", julia_folder) !== nothing
end

"""
Used in the build step of a new app created using PackageCompiler.
Copies all assets from pkg to the new app.
"""
function copy_assets(src_module::Module, appfolder)
    app_assetpath = joinpath(appfolder, "assets")
    mkpath(app_assetpath)
    dst = joinpath(pkgdir(src_module), "assets")
    src = joinpath(app_assetpath, string(src_module))
    cp(dst, src)
end

end

ext: jl

max_issues_repo_name: nielsls/Assets

max_forks_repo_name: nielsls/Assets


text:
@polyvar x z
@variable(model, t)
p = x^4 + x^2 - 3*x^2*z^2 + z^6 - t
@constraint(model, p >= 0)
@objective(model, Max, t)
optimize!(model)
println("Solution: $(value(t))")
# Returns the lower bound -.17700

max_stars_repo_name: sums-of-squares/tools

max_issues_repo_name: sums-of-squares/tools

max_forks_repo_name: sums-of-squares/tools


text:
function hfun_bar(vname)
  val = Meta.parse(vname[1])
  return round(sqrt(val), digits=2)
end

function hfun_m1fill(vname)
  var = vname[1]
  return pagevar("index", var)
end

function lx_baz(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return uppercase(brace_content)
end

"""
    {{news}}

Plug in the list of news contained in the `/news/` folder.
"""
function hfun_news()
    author_badges = globvar("author_badges")

    # Changing page width in "News" section
    io = IOBuffer()
    write(io, """""")

    # Add news
    curyear = year(Dates.today())
    write(io, """
""")
    for year in curyear:-1:2022
        ys = "$year"
        year < curyear
        for month in 12:-1:1
            ms = "0"^(month < 10) * "$month"
            base = joinpath("news", ys, ms)
            isdir(base) || continue
            posts = filter!(p -> endswith(p, ".md"), readdir(base))
            days  = zeros(Int, length(posts))
            lines = Vector{String}(undef, length(posts))
            for (i, post) in enumerate(posts)
                ps  = splitext(post)[1]
                url = "/news/$ys/$ms/$ps/"
                surl = strip(url, '/')
                post_title = pagevar(surl, :post_title)
                short_descr = pagevar(surl, :short_descr)
                post_author = pagevar(surl, :post_author)
                width = rand(250:290)
                pubdate = pagevar(surl, :published)
                if isnothing(pubdate)
                    date    = "$ys-$ms-01"
                    days[i] = 1
                else
                    date    = Date(pubdate, dateformat"d U Y")
                    days[i] = day(date)
                end
                lines[i] = """"""
                lines[i] *= """
"""
                lines[i] *= """$post_title $date
"""
                lines[i] *= """
$short_descr
"""
                lines[i] *= """$post_author"""
                for b in author_badges[post_author]
                      lines[i] *= """"""
                end
                lines[i] *= """
"""
            end
            # sort by day
            foreach(line -> write(io, line), lines[sortperm(days, rev=true)])
        end
    end
    write(io, """
""")
    return String(take!(io))
end

"""
    {{ addcomments }}

Add a comment widget, managed by utterances <https://utteranc.es>.
"""
function hfun_addcomments()
    html_str = """

    """
    return html_str
end


text:
###### This file is part of the MendelImpute.jl package.
###### It contains relevant code to impute and output a genotype matrix

"""
    write(outfile, X, compressed_haplotypes, X_sampleID, XtoH_idx)

Writes imputed `X` into `outfile`. All genotypes in `outfile` are non-missing.

# Inputs
- `outfile`: Output file name (ending in `.vcf.gz` or `.vcf`)
- `X`: Imputed matrix. If `X` is a matrix, output VCF will be unphased. Otherwise if `X` is a Tuple
    of matrix (i.e. `X = X1 + X2`), then all output will be phased.
- `compressed_haplotypes`: A `CompressedHaplotypes` object
- `X_sampleID`: Sample ID of imputation target
- `snp_score`: Imputation score for each typed SNP. `1` is best, `0` is worse
- `XtoH_idx`: Position of typed SNPs in complete set of SNPs
- `impute`: Boolean indicating whether to impute untyped SNPs (default true)

# Notes
Here the writing routine is emulating `write_dlm` in Base at
https://github.com/JuliaLang/julia/blob/3608c84e6093594fe86923339fc315231492484c/stdlib/DelimitedFiles/src/DelimitedFiles.jl#L736
"""
function Base.write(
    outfile::AbstractString,
    X::Union{AbstractMatrix, Tuple{AbstractMatrix, AbstractMatrix}},
    compressed_haplotypes::CompressedHaplotypes,
    X_sampleID::AbstractVector,
    snp_score::AbstractVector,
    XtoH_idx::AbstractVector,
    impute::Bool = true
    )
    threads = Threads.nthreads()
    snps = typeof(X) <: AbstractMatrix ? size(X, 1) : size(X[1], 1)
    len = div(snps, threads)
    files = ["tmp$i.vcf.gz" for i in 1:threads]
    typed = falses(snps)
    impute ? (typed[XtoH_idx] .= true) : (typed .= true)

    # retrieve reference file information
    chr = impute ? compressed_haplotypes.chr : @view(compressed_haplotypes.chr[XtoH_idx])
    pos = impute ? compressed_haplotypes.pos : @view(compressed_haplotypes.pos[XtoH_idx])
    ids = impute ? compressed_haplotypes.SNPid : @view(compressed_haplotypes.SNPid[XtoH_idx])
    ref = impute ? compressed_haplotypes.refallele : @view(compressed_haplotypes.refallele[XtoH_idx])
    alt = impute ? compressed_haplotypes.altallele : @view(compressed_haplotypes.altallele[XtoH_idx])

    # write minimal meta information to outfile
    io = [openvcf(files[i], "w") for i in 1:threads]
    pb = [PipeBuffer() for _ in 1:threads]
    print(pb[1], "##fileformat=VCFv4.2\n")
    print(pb[1], "##source=MendelImpute\n")
    print(pb[1], "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n")
    print(pb[1], "##INFO=<ID=IMPQ,Number=0,Type=Float,Description=" *
        "\"Quality of marker. 0 is best, larger is worse.\">\n")
    print(pb[1], "##INFO=<ID=IMP,Number=0,Type=Flag,Description=\"Imputed marker\">\n")

    # header line should match reffile (i.e. sample ID's should match)
    print(pb[1], "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT")
    for id in X_sampleID
        print(pb[1], "\t", id)
    end
    print(pb[1], "\n")
    bytesavailable(pb[1]) > 1048576 && write(io[1], take!(pb[1]))

    # each thread writes `len` SNPs
    pmeter = Progress(snps, 5, "Writing to file...")
    Threads.@threads for t in 1:threads
        id = Threads.threadid()
        cur_ranges = (id == threads ?
            ((threads-1)*len+1:snps) : (1:len) .+ (t-1)*len)

        @inbounds for i in cur_ranges
            # write meta info (chrom/pos/snpid/ref/alt/imputation-quality)
            print(pb[id], chr[i], "\t", string(pos[i]), "\t", ids[i][1], "\t",
                ref[i], "\t", alt[i][1], "\t.\tPASS\t")
            print(pb[id], "IMPQ=", round(snp_score[i], digits=3))
            typed[i] ? print(pb[id], "\tGT") : print(pb[id], ";IMP\tGT")
            # print ith record
            write_snp!(pb[id], X, i)
            bytesavailable(pb[id]) > 1048576 && write(io[id], take!(pb[id]))
            next!(pmeter)
        end
        write(io[id], take!(pb[id]))
    end
    close.(io); close.(pb) # close io and buffer

    # concatenate all files into 1 VCF file
    run(pipeline(`cat $files`, stdout=outfile))

    # delete intermediate files
    for i in 1:threads
        rm("tmp$i.vcf.gz", force=true)
    end

    return nothing
end

"""
Helper function for saving a record (SNP), not tracking phase information.
"""
function write_snp!(pb::IOBuffer, X::AbstractMatrix, i::Int)
    x = @view(X[i, :]) # current record
    n = length(x)
    @inbounds for j in 1:n
        if x[j] == 0
            print(pb, "\t0/0")
        elseif x[j] == 1
            print(pb, "\t1/0")
        elseif x[j] == 2
            print(pb, "\t1/1")
        else
            error("imputed genotypes can only be 0, 1, 2 but got $(x[j])")
        end
    end
    print(pb, "\n")
    nothing
end

"""
Helper function for saving a record (SNP), tracking phase information.
Here `X = X1 + X2`.
"""
function write_snp!(pb::IOBuffer, X::Tuple{AbstractMatrix, AbstractMatrix}, i::Int)
    X1, X2 = X[1], X[2]
    x1 = @view(X1[i, :])
    x2 = @view(X2[i, :])

    n = length(x1)
    @assert n == length(x2)
    @inbounds for j in 1:n
        if x1[j] == x2[j] == 0
            print(pb, "\t0|0")
        elseif x1[j] == 0 && x2[j] == 1
            print(pb, "\t0|1")
        elseif x1[j] == 1 && x2[j] == 0
            print(pb, "\t1|0")
        elseif x1[j] == 1 && x2[j] == 1
            print(pb, "\t1|1")
        else
            error("phased genotypes can only be 0|0, 0|1, 1|0 or 1|1 but
                got $(x1[j])|$(x2[j])")
        end
    end
    print(pb, "\n")
    nothing
end

"""
    impute!(X1, X2, H, phase, [impute_untyped])

Imputes `X = X1 + X2` completely using haplotype segments of `H`, where segments
information are stored in `phase`. `X1` is strand1 and `X2` is strand 2. If
`impute_untyped=false`, untyped SNPs will not be imputed.

Non-missing entries in `X` can be different after imputation.
"""
function impute!(
    X1::AbstractMatrix,
    X2::AbstractMatrix,
    compressed_Hunique::CompressedHaplotypes,
    phase::Vector{HaplotypeMosaicPair},
    impute_untyped::Bool=true
    )

    fill!(X1, 0)
    fill!(X2, 0)

    # loop over individuals
    @inbounds for i in 1:size(X1, 2)
        # strand 1
        for segment in 1:length(phase[i].strand1.start)
            Xrange, haplotype = _get_impute_ranges(segment, phase[i].strand1,
                compressed_Hunique, impute_untyped=impute_untyped)
            X1[Xrange, i] .= haplotype
        end

        # strand 2
        for segment in 1:length(phase[i].strand2.start)
            Xrange, haplotype = _get_impute_ranges(segment, phase[i].strand2,
                compressed_Hunique, impute_untyped=impute_untyped)
            X2[Xrange, i] .= haplotype
        end
    end
end

"""
    _get_impute_ranges(segment, strand, compressed_Hunique, impute_untyped,
        num_typed_snps)

Helper function for impute! that computes the range of `X` and `H` in window `w`

# Inputs:
- `segment`: Current segment of a `HaplotypeMosaic`
- `strand`: A `HaplotypeMosaic` (1 strand of haplotype for a person, contains
    many segments)
- `compressed_Hunique`: A `CompressedHaplotypes` object
- `impute_untyped`: If `true`, all untyped SNPs will be considered in indexing.
- `haplabel_typedonly`: If `true`, `strand.haplotypelabel` stores haplotype
    index with respect to `compressed_Hunique.CW`. Otherwise index is with
    respect to `compressed_Hunique.CW_typed`
"""
function _get_impute_ranges(
    segment::Int,
    strand::HaplotypeMosaic,
    compressed_Hunique::CompressedHaplotypes;
    impute_untyped::Bool = true,
    haplabel_typedonly::Bool = false
    )
    window = strand.window[segment]
    is_last_segment = segment == length(strand.start)
    total_snps = impute_untyped ? strand.length :
        last(compressed_Hunique.X_window_range[end])

    # get X range
    Xstart = strand.start[segment]
    X_end = is_last_segment ? total_snps : strand.start[segment + 1] - 1
    Xrange = Xstart:X_end

    # get haplotype vector
    H = impute_untyped ? compressed_Hunique.CW[window].uniqueH :
                         compressed_Hunique.CW_typed[window].uniqueH
    offset = impute_untyped ? compressed_Hunique.Hstart[window] :
        first(compressed_Hunique.X_window_range[window])
    Hstart = strand.start[segment] - offset + 1
    Hrange = Hstart:(Hstart + length(Xrange) - 1)
    Hlabel = strand.haplotypelabel[segment]
    haplabel_typedonly && (Hlabel = unique_all_idx_to_unique_typed_idx(Hlabel,
        window, compressed_Hunique))
    haplotype = @view(H[Hrange, Hlabel])

    return Xrange, haplotype
end

function impute!(
    X1::AbstractMatrix,
    X2::AbstractMatrix,
    H::AbstractMatrix,
    phase::Vector{HaplotypeMosaicPair},
    )

    fill!(X1, 0)
    fill!(X2, 0)

    # loop over individuals
    for i in 1:size(X1, 2)
        for s in 1:(length(phase[i].strand1.start) - 1)
            idx = phase[i].strand1.start[s]:(phase[i].strand1.start[s + 1] - 1)
            X1[idx, i] = H[idx, phase[i].strand1.haplotypelabel[s]]
        end
        idx = phase[i].strand1.start[end]:phase[i].strand1.length
        X1[idx, i] = H[idx, phase[i].strand1.haplotypelabel[end]]
        for s in 1:(length(phase[i].strand2.start) - 1)
            idx = phase[i].strand2.start[s]:(phase[i].strand2.start[s + 1] - 1)
            X2[idx, i] += H[idx, phase[i].strand2.haplotypelabel[s]]
        end
        idx = phase[i].strand2.start[end]:phase[i].strand2.length
        X2[idx, i] += H[idx, phase[i].strand2.haplotypelabel[end]]
    end
end

"""
    impute_discard_phase!(X, H, phase, [impute_untyped])

Imputes missing entries of `X` using corresponding haplotypes `H` via `phase`
information. Non-missing entries in `X` will not change, but X and H has to be
aligned. This does NOT preserve phase information.
"""
function impute_discard_phase!(
    X::AbstractMatrix,
    compressed_Hunique::CompressedHaplotypes,
    phase::Vector{HaplotypeMosaicPair},
    impute_untyped::Bool=true
    )

    p, n = size(X)

    # for person in 1:n
    Threads.@threads for person in 1:n
        @inbounds for snp in 1:p
            if ismissing(X[snp, person])
                #find which segment the snp is located
                hap1_segment = searchsortedlast(phase[person].strand1.start,snp)
                hap2_segment = searchsortedlast(phase[person].strand2.start,snp)

                #find haplotype pair in corresponding window for this segment
                h1 = phase[person].strand1.haplotypelabel[hap1_segment]
                h2 = phase[person].strand2.haplotypelabel[hap2_segment]
                w1 = phase[person].strand1.window[hap1_segment]
                w2 = phase[person].strand2.window[hap2_segment]
                i1 = snp - compressed_Hunique.Hstart[w1] + 1
                i2 = snp - compressed_Hunique.Hstart[w2] + 1

                # imputation step
                H1 = impute_untyped ? compressed_Hunique.CW[w1].uniqueH :
                    compressed_Hunique.CW_typed[w1].uniqueH
                H2 = impute_untyped ? compressed_Hunique.CW[w2].uniqueH :
                    compressed_Hunique.CW_typed[w2].uniqueH
                X[snp, person] = H1[i1, h1] + H2[i2, h2]
            end
        end
    end

    return nothing
end

"""
    update_marker_position!(phaseinfo, tgtfile, reffile)

Converts `phaseinfo`'s strand1 and strand2's starting position in
terms of matrix rows of `X` to starting position in terms matrix
rows in `H`.
"""
function update_marker_position!(
    phaseinfo::Vector{HaplotypeMosaicPair},
    XtoH_idx::AbstractVector,
    )
    people = length(phaseinfo)

    @inbounds for j in 1:people
        # update strand1's starting position
        for (i, idx) in enumerate(phaseinfo[j].strand1.start)
            phaseinfo[j].strand1.start[i] = XtoH_idx[idx]
        end
        # update strand2's starting position
        for (i, idx) in enumerate(phaseinfo[j].strand2.start)
            phaseinfo[j].strand2.start[i] = XtoH_idx[idx]
        end
    end

    # update first starting position
    @inbounds for j in 1:people
        phaseinfo[j].strand1.start[1] = 1
        phaseinfo[j].strand2.start[1] = 1
    end

    return nothing
end

"""
    untyped_snpscore(total_snps, typed_snp_scores, typed_index)

For each untyped SNP, average the nearest 2 typed SNP's quality score and
return a vector of quality scores for all SNPs, typed and untyped. 0 is best,
larger is worse.

# Inputs
- `total_snps`: Total number of SNPs, typed and untyped
- `typed_snp_scores`: Vector of scores for each typed SNP
- `typed_index`: Typed SNP's indices
"""
function untyped_snpscore(
    total_snps::Int,
    typed_snp_scores::AbstractVector,
    typed_index::AbstractVector
    )

    # copy typed SNPs' quality score into vector of complete SNPs
    complete_snpscore = Vector{eltype(typed_snp_scores)}(undef, total_snps)
    copyto!(@view(complete_snpscore[typed_index]), typed_snp_scores)

    # all untyped SNPs before first typed SNPs gets same quality score
    cur_range = 1:(typed_index[1] - 1)
    complete_snpscore[cur_range] .= typed_snp_scores[1]

    # loop through every segments of untyped SNPs
    for i in 2:length(typed_index)
        cur_range = (typed_index[i - 1] + 1):(typed_index[i] - 1)
        avg_score = (typed_snp_scores[i - 1] + typed_snp_scores[i]) / 2
        complete_snpscore[cur_range] .= avg_score
    end

    # last segment of untyped SNPs
    cur_range = (typed_index[end] + 1):total_snps
    complete_snpscore[cur_range] .= typed_snp_scores[end]

    return complete_snpscore
end

"""
    typed_snpscore(snps, typed_snp_scores, typed_index)

For each typed SNP, compute its average least square error over all samples.
Also for each sample, compute its total least squares error. In both cases
0 is best and larger is worse.

# Inputs
- `X`: Genotype matrix. Each row is a typed SNP and each column is a person.
- `phase`: A vector of `HaplotypeMosaicPair` keeping track of each person's
    phase information. `Haplotypelabels` point to complete haplotype set
- `compressed_haplotypes`: A `CompressedHaplotypes` object
"""
function typed_snpscore(
    X::AbstractMatrix,
    phase::Vector{HaplotypeMosaicPair},
    compressed_haplotypes::CompressedHaplotypes,
    )
    snps, samples = size(X)
    scores = zeros(snps)
    Ximp = zeros(eltype(X), snps)
    Xerr = zeros(Float64, snps)
    snp_error = zeros(Float64, snps)
    sample_error = zeros(Float64, samples)
    missing_counter = zeros(Int, snps)

    @inbounds for i in 1:samples
        fill!(Ximp, 0)
        # add strand1 to Ximp
        for segment in 1:length(phase[i].strand1.start)
            Xrange, haplotype = _get_impute_ranges(segment, phase[i].strand1,
                compressed_haplotypes, impute_untyped=false,
                haplabel_typedonly=true)
            Ximp[Xrange] .= haplotype
        end

        # add strand2 to Ximp
        for segment in 1:length(phase[i].strand2.start)
            Xrange, haplotype = _get_impute_ranges(segment, phase[i].strand2,
                compressed_haplotypes, impute_untyped=false,
                haplabel_typedonly=true)
            Ximp[Xrange] += haplotype
        end

        # accumulate error
        sample_missingness = 0
        for j in 1:snps
            if ismissing(X[j, i])
                missing_counter[j] += 1
                sample_missingness += 1
                continue
            end
            err = abs2(X[j, i] - Ximp[j])
            snp_error[j] += err
            sample_error[i] += err
        end

        # normalize per-sample error
        sample_error[i] = 1 - 0.25(sample_error[i] / (snps - sample_missingness))
    end

    for i in 1:snps
        @inbounds snp_error[i] = 1 - 0.25(snp_error[i] /
            (samples - missing_counter[i]))
    end

    return snp_error, sample_error
end

function write_sample_error(
    outfile::AbstractString,
    samplescore::AbstractVector,
    sampleID::AbstractVector
    )
    # write per-sample error to output file
    if endswith(outfile, ".jlso")
        strip_chr = 4
    elseif endswith(outfile, ".vcf")
        strip_chr = 3
    else # .vcf.gz
        strip_chr = 6
    end
    error_filename = outfile[1:end-strip_chr]

    open(error_filename * "sample.error", "w") do io
        println(io, "sampleID,error")
        for i in eachindex(sampleID)
            @inbounds println(io, sampleID[i], ",", samplescore[i])
        end
    end

    return nothing
end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
@testset "Vax Float D" begin
    d8_vax = [  0x0000000000004080,
                0x000000000000c080,
                0x0000000000004160,
                0x000000000000c160,
                0x68c0a2210fda4149,
                0x68c0a2210fdac149,
                0x48d81abbbdc27df0,
                0x48d81abbbdc2fdf0,
                0x5c7814541cea0308,
                0x5c7814541cea8308,
                0xcee814620652409e,
                0xcee814620652c09e]

    d8_ieee = Array{Float64}([  one(Float64),
                                -one(Float64),
                                3.5,
                                -3.5,
                                Float64(pi),
                                -Float64(pi),
                                1.0e37,
                                -1.0e37,
                                9.9999999999999999999999999e-38,
                                -9.9999999999999999999999999e-38,
                                1.2345678901234500000000000000,
                                -1.2345678901234500000000000000 ])

    @testset "Conversion..." begin
        for (vax, ieee) in zip(d8_vax, d8_ieee)
            @test VaxFloatD(vax) == VaxFloatD(ieee)
            @test convert(Float64, VaxFloatD(vax)) == ieee
        end
    end

    @testset "Promotion..." begin
        for t in [subtypes(VaxInt); subtypes(VaxFloat); Int8; Int16; Int32; Int64; Int128; Float16; Float32; Float64]
            @test isa(one(t)*VaxFloatD(1), Float64)
        end
        @test isa(one(BigInt)*VaxFloatD(1), BigFloat)
        @test isa(one(BigFloat)*VaxFloatD(1), BigFloat)

        un = one(VaxFloatD)
        @test promote(un, un, un) == (1.0, 1.0, 1.0)
        @test promote(un, un, un, un) == (1.0, 1.0, 1.0, 1.0)
    end

    @testset "Number definitions" begin
        @test floatmax(VaxFloatD) == typemax(VaxFloatD)
        @test floatmin(VaxFloatD) == typemin(VaxFloatD)

        @test zero(VaxFloatD) == 0
        @test one(VaxFloatD) == 1
    end

    @testset "Edge cases" begin
        # Reserved Vax floating point operand
        @test_throws InexactError convert(Float64, VaxFloatD(UInt64(0x8000)))

        # Inf and NaN should error too
        @test_throws InexactError VaxFloatD(Inf64)
        @test_throws InexactError VaxFloatD(-Inf64)
        @test_throws InexactError VaxFloatD(NaN64)

        # Both IEEE zeros should be converted to Vax true zero
        @test VaxFloatD(-0.0) === VaxFloatD(0.0) === zero(VaxFloatD)

        # Dirty zero
        @test convert(Float64, VaxFloatD(UInt64(0x08))) === zero(Float64)

        # Numbers smaller than floatmin(VaxFloatD) should underflow
        @test VaxFloatD(prevfloat(convert(Float64, floatmin(VaxFloatD)))) === zero(VaxFloatD)
        @test VaxFloatD(convert(Float64, floatmin(VaxFloatD))) === floatmin(VaxFloatD)

        # Numbers larger than floatmax(VaxFloatD) should error
        @test_throws InexactError VaxFloatD(nextfloat(convert(Float64, floatmax(VaxFloatD))))

        # Because the D Float as more precision, the conversion to Float64 and back to D Float will not be circular
        # @test VaxFloatD(convert(Float64, floatmax(VaxFloatD))) === floatmax(VaxFloatD)
    end

    @testset "IO" begin
        io = IOBuffer(reinterpret(UInt8, ones(VaxFloatD, 4)))
        @test read(io, VaxFloatD) === one(VaxFloatD)
        @test read!(io, Vector{VaxFloatD}(undef, 3)) == ones(VaxFloatD, 3)
    end
end


text:
#### Similarity and dissimilarity measures between images ####

"""
`hausdorff_distance(imgA, imgB)` is the modified Hausdorff distance
between binary images (or point sets).

References
----------
Dubuisson, M-P; Jain, A. K., 1994. A Modified Hausdorff Distance for
Object-Matching.
"""
function hausdorff_distance(imgA::AbstractArray, imgB::AbstractArray)
  @assert size(imgA) == size(imgB) "Images must have the same dimensions"

  # trivial case
  isequal(imgA, imgB) && return 0.

  A, B = find(imgA), find(imgB)

  # return if there is no object to match
  (isempty(A) || isempty(B)) && return Inf

  # grid coordinates (ndims by npoints)
  A = hcat([Float64[ind2sub(size(imgA), a)...] for a in A]...)
  B = hcat([Float64[ind2sub(size(imgB), b)...] for b in B]...)

  m = size(A, 2); n = size(B, 2)

  D = zeros(m, n)
  for j=1:n, i=1:m
    @inbounds D[i,j] = norm(A[:,i] - B[:,j])
  end

  dAB = mean(minimum(D, 2))
  dBA = mean(minimum(D, 1))

  max(dAB, dBA)
end


text:
#####
##### `AbstractDifferential`
#####

"""
The subtypes of `AbstractDifferential` define a custom \"algebra\" for chain
rule evaluation that attempts to factor various features like complex derivative
support, broadcast fusion, zero-elision, etc. into nicely separated parts.

All subtypes of `AbstractDifferential` implement the following operations:

`+(a, b)`: linearly combine differential `a` and differential `b`

`*(a, b)`: multiply the differential `a` by the differential `b`

`Base.conj(x)`: complex conjugate of the differential `x`

`extern(x)`: convert `x` into an appropriate non-`AbstractDifferential` type for
use outside of `ChainContext`.

Valid arguments to these operations are `T` where `T<:AbstractDifferential`, or
where `T` has proper `+` and `*` implementations.

Additionally, all subtypes of `AbstractDifferential` support `Base.iterate` and
`Base.Broadcast.broadcastable(x)`.
"""
abstract type AbstractDifferential end

Base.:+(x::AbstractDifferential) = x

"""
    extern(x)

Return `x` converted to an appropriate non-`AbstractDifferential` type, for use
with external packages that might not handle `AbstractDifferential` types.

Note that this function may return an alias (not necessarily a copy) to data
wrapped by `x`, such that mutating `extern(x)` might mutate `x` itself.
"""
@inline extern(x) = x

@inline Base.conj(x::AbstractDifferential) = x

#####
##### `Wirtinger`
#####

"""
    Wirtinger(primal::Union{Number,AbstractDifferential},
              conjugate::Union{Number,AbstractDifferential})

Returns a `Wirtinger` instance representing the complex differential:

```
df = ∂f/∂z * dz + ∂f/∂z̄ * dz̄
```

where `primal` corresponds to `∂f/∂z * dz` and `conjugate` corresponds to `∂f/∂z̄ * dz̄`.

The two fields of the returned instance can be accessed generically via the
[`wirtinger_primal`](@ref) and [`wirtinger_conjugate`](@ref) methods.
"""
struct Wirtinger{P,C} <: AbstractDifferential
    primal::P
    conjugate::C
    function Wirtinger(primal::Union{Number,AbstractDifferential},
                       conjugate::Union{Number,AbstractDifferential})
        return new{typeof(primal),typeof(conjugate)}(primal, conjugate)
    end
end

wirtinger_primal(x::Wirtinger) = x.primal
wirtinger_primal(x) = x

wirtinger_conjugate(x::Wirtinger) = x.conjugate
wirtinger_conjugate(::Any) = Zero()

extern(x::Wirtinger) = throw(ArgumentError("`Wirtinger` cannot be converted to an external type."))

Base.Broadcast.broadcastable(w::Wirtinger) = Wirtinger(broadcastable(w.primal),
                                                       broadcastable(w.conjugate))

Base.iterate(x::Wirtinger) = (x, nothing)
Base.iterate(::Wirtinger, ::Any) = nothing

# TODO: define `conj` for` `Wirtinger`
Base.conj(x::Wirtinger) = throw(MethodError(conj, x))

#####
##### `Zero`
#####

"""
    Zero()
The additive identity for differentials.
This is basically the same as `0`.
"""
struct Zero <: AbstractDifferential end

extern(x::Zero) = false  # false is a strong 0. E.g. `false * NaN = 0.0`

Base.Broadcast.broadcastable(::Zero) = Ref(Zero())

Base.iterate(x::Zero) = (x, nothing)
Base.iterate(::Zero, ::Any) = nothing

#####
##### `DoesNotExist`
#####

"""
    DoesNotExist()

This differential indicates that the derivative Does Not Exist (D.N.E).
This is not the cast that it is not implemented, but rather that it mathematically
is not defined.
"""
struct DoesNotExist <: AbstractDifferential end

function extern(x::DoesNotExist)
    throw(ArgumentError("Derivative does not exit. Cannot be converted to an external type."))
end

Base.Broadcast.broadcastable(::DoesNotExist) = Ref(DoesNotExist())

Base.iterate(x::DoesNotExist) = (x, nothing)
Base.iterate(::DoesNotExist, ::Any) = nothing

#####
##### `One`
#####

"""
     One()
The Differential which is the multiplicative identity.
Basically, this represents `1`.
"""
struct One <: AbstractDifferential end

extern(x::One) = true  # true is a strong 1.

Base.Broadcast.broadcastable(::One) = Ref(One())

Base.iterate(x::One) = (x, nothing)
Base.iterate(::One, ::Any) = nothing

#####
##### `AbstractThunk
#####
abstract type AbstractThunk <: AbstractDifferential end

Base.Broadcast.broadcastable(x::AbstractThunk) = broadcastable(extern(x))

@inline function Base.iterate(x::AbstractThunk)
    externed = extern(x)
    element, state = iterate(externed)
    return element, (externed, state)
end

@inline function Base.iterate(::AbstractThunk, (externed, state))
    element, new_state = iterate(externed, state)
    return element, (externed, new_state)
end

#####
##### `Thunk`
#####

"""
    Thunk(()->v)
A thunk is a deferred computation.
It wraps a zero argument closure that when invoked returns a differential.
`@thunk(v)` is a macro that expands into `Thunk(()->v)`.

Calling a thunk, calls the wrapped closure.
`extern`ing thunks applies recursively, it also externs the differial that the closure returns.
If you do not want that, then simply call the thunk

```
julia> t = @thunk(@thunk(3))
Thunk(var"##7#9"())

julia> extern(t)
3

julia> t()
Thunk(var"##8#10"())

julia> t()()
3
```

### When to `@thunk`?
When writing `rrule`s (and to a lesser exent `frule`s), it is important to `@thunk`
appropriately.
Propagation rule's that return multiple derivatives are not able to do all the computing themselves.
 By `@thunk`ing the work required for each, they then compute only what is needed.

#### So why not thunk everything?
`@thunk` creates a closure over the expression, which (effectively) creates a `struct`
with a field for each variable used in the expression, and call overloaded.

Do not use `@thunk` if this would be equal or more work than actually evaluating the expression itself. Examples being:
- The expression wrapping something in a `struct`, such as `Adjoint(x)` or `Diagonal(x)`
- The expression being a constant
- The expression being itself a `thunk`
- The expression being from another `rrule` or `frule` (it would be `@thunk`ed if required by the defining rule already)
"""
struct Thunk{F} <: AbstractThunk
    f::F
end

"""
    @thunk expr

Define a [`Thunk`](@ref) wrapping the `expr`, to lazily defer its evaluation.
"""
macro thunk(body)
    # Basically `:(Thunk(() -> $(esc(body))))` but use the location where it is defined.
    # so we get useful stack traces if it errors.
    func = Expr(:->, Expr(:tuple), Expr(:block, __source__, body))
    return :(Thunk($(esc(func))))
end

# have to define this here after `@thunk` and `Thunk` is defined
Base.conj(x::AbstractThunk) = @thunk(conj(extern(x)))

(x::Thunk)() = x.f()
@inline extern(x::Thunk) = extern(x())

Base.show(io::IO, x::Thunk) = println(io, "Thunk($(repr(x.f)))")

"""
    InplaceableThunk(val::Thunk, add!::Function)

A wrapper for a `Thunk`, that allows it to define an inplace `add!` function,
which is used internally in `accumulate!(Δ, ::InplaceableThunk)`.

`add!` should be defined such that: `ithunk.add!(Δ) = Δ .+= ithunk.val`
but it should do this more efficently than simply doing this directly.
(Otherwise one can just use a normal `Thunk`).

Most operations on an `InplaceableThunk` treat it just like a normal `Thunk`;
and destroy its inplacability.
"""
struct InplaceableThunk{T<:Thunk, F} <: AbstractThunk
    val::T
    add!::F
end

(x::InplaceableThunk)() = x.val()
@inline extern(x::InplaceableThunk) = extern(x.val)

function Base.show(io::IO, x::InplaceableThunk)
    println(io, "InplaceableThunk($(repr(x.val)), $(repr(x.add!)))")
end

# The real reason we have this:
accumulate!(Δ, ∂::InplaceableThunk) = ∂.add!(Δ)
store!(Δ, ∂::InplaceableThunk) = ∂.add!((Δ.*=false))  # zero it, then add to it.

"""
    NO_FIELDS

Constant for the reverse-mode derivative with respect to a structure that has no fields.
The most notable use for this is for the reverse-mode derivative with respect to the
function itself, when that function is not a closure.
"""
const NO_FIELDS = DoesNotExist()

"""
    refine_differential(𝒟::Type, der)

Converts, if required, a differential object `der`
(e.g. a `Number`, `AbstractDifferential`, `Matrix`, etc.),
to another  differential that is more suited for the domain given by the type 𝒟.
Often this will behave as the identity function on `der`.
"""
function refine_differential(::Type{<:Union{<:Real, AbstractArray{<:Real}}}, w::Wirtinger)
    return wirtinger_primal(w) + wirtinger_conjugate(w)
end
refine_differential(::Any, der) = der  # most of the time leave it alone.


text:
# create a dummy basis that is easy to track
f11(theta) = theta[1]^2 + 5
f12(theta) = theta[1] - 9
f21(theta) = theta[1]^2 * (theta[2]^3 + 4)
f22(theta) = theta[2]
f23(theta) = theta[1] * theta[2]
f31(theta) = theta[1] * theta[2] * theta[3]

struct TestBasis <: OrderedBasis
  component::Int
end

function (tuple::TestBasis)(theta::AbstractArray{<:Real})
  c = tuple.component
  funcs = c == 1 ? [f11, f12] : (c == 2 ? [f21, f22, f23] : [f31])
  map(f -> f(theta), funcs)
end

Base.length(tuple::TestBasis) = tuple.component == 1 ? 2 : (tuple.component == 2 ? 3 : 1)

@testset "TransportFamily linearity" begin
  T = TransportFamily(map(TestBasis, 1:3))
  for _=1:100
    gam1 = [randn(2), randn(3), randn(1)]
    gam2 = [randn(2), randn(3), randn(1)]
    scale = randn()
    gam3 = @. gam1 + scale * gam2
    for _=1:100
      theta = randn(3)
      @test T(gam1)(theta) + scale * T(gam2)(theta) ≈ T(gam3)(theta)
    end
  end
end

@testset "TransportFamily typing" begin
  T = TransportFamily(map(TestBasis, 1:3))
  @test isa(T, TransportFamily)
  @test isa(T([zeros(2), zeros(3), zeros(1)]), TransportMap)
end

include("polynomials.jl")

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
import RestrictedBoltzmannMachines as RBMs
using RestrictedBoltzmannMachines: energy
using Test: @test, @testset, @inferred
using WhitenedRBMs: shift_fields, shift_fields!

function energy_shift(offset::AbstractArray, x::AbstractArray)
    @assert size(offset) == size(x)[1:ndims(offset)]
    if ndims(offset) == ndims(x)
        return -sum(offset .* x)
    elseif ndims(offset) < ndims(x)
        ΔE = -sum(offset .* x; dims=1:ndims(offset))
        return reshape(ΔE, size(x)[(ndims(offset) + 1):end])
    end
end

@testset "shift_fields" begin
    N = (3, 4)
    layers = (
        RBMs.Binary(randn(N...)),
        RBMs.Spin(randn(N...)),
        RBMs.Potts(randn(N...)),
        RBMs.Gaussian(randn(N...), rand(N...)),
        RBMs.ReLU(randn(N...), rand(N...)),
        RBMs.dReLU(randn(N...), randn(N...), rand(N...), rand(N...)),
        RBMs.pReLU(randn(N...), rand(N...), randn(N...), rand(N...) .- 0.5),
        RBMs.xReLU(randn(N...), rand(N...), randn(N...), randn(N...)),
    )
    for layer in layers
        offset = randn(size(layer)...)
        x = RBMs.sample_from_inputs(layer, randn(size(layer)..., 2, 3))
        layer_shifted = @inferred shift_fields(layer, offset)
        @test energy(layer_shifted, x) ≈ energy(layer, x) + energy_shift(offset, x)
    end
end

@testset "shift_fields!" begin
    N = (3, 4)
    layers = (
        RBMs.Binary(randn(N...)),
        RBMs.Spin(randn(N...)),
        RBMs.Potts(randn(N...)),
        RBMs.Gaussian(randn(N...), rand(N...)),
        RBMs.ReLU(randn(N...), rand(N...)),
        RBMs.dReLU(randn(N...), randn(N...), rand(N...), rand(N...)),
        RBMs.pReLU(randn(N...), rand(N...), randn(N...), rand(N...) .- 0.5),
        RBMs.xReLU(randn(N...), rand(N...), randn(N...), randn(N...)),
    )
    for layer in layers
        offset = randn(size(layer)...)
        x = RBMs.sample_from_inputs(layer, randn(size(layer)..., 2, 3))
        E = energy(layer, x)
        @inferred shift_fields!(layer, offset)
        @test energy(layer, x) ≈ E + energy_shift(offset, x)
    end
end


text:
"""
    zygote_friendly_map(f, x)

This version of map is a bit weird. It makes slightly stronger assumptions about the nature
of what you're allowed to pass in to it than `Base.map` does and, in return, you get much
improved performance when used in conjunction with `Zygote`.

# Assumptions.
- No globals are used in `f`. This means that `TemporalGPs.NoContext` can be employed.
- `f` has no fields. If you've got data to share across elements, use a `Fill`.
- Similarly, `f` has no mutable state (follows from the above).
- `f` doesn't mutate its argument.
"""
zygote_friendly_map(f, x) = dense_zygote_friendly_map(f, x)

function dense_zygote_friendly_map(f, x)

    # Perform first iteration.
    y_1 = f(_getindex(x, 1))

    # Allocate for outputs.
    ys = Array{typeof(y_1)}(undef, size(x))
    ys[1] = y_1

    # Perform remainder of iterations.
    for n in 2:length(x)
        ys[n] = f(_getindex(x, n))
    end

    return ys
end

function Zygote._pullback(::AContext, ::typeof(dense_zygote_friendly_map), f, x)

    # Perform first iteration.
    y_1, pb_1 = Zygote._pullback(NoContext(), f, _getindex(x, 1))

    # Allocate for outputs.
    ys = Array{typeof(y_1)}(undef, size(x))
    ys[1] = y_1

    # Allocate for pullbacks.
    pbs = Array{typeof(pb_1)}(undef, size(x))
    pbs[1] = pb_1

    for n in 2:length(x)
        y, pb = Zygote._pullback(NoContext(), f, _getindex(x, n))
        ys[n] = y
        pbs[n] = pb
    end

    function zygote_friendly_map_pullback(Δ)
        Δ === nothing && return

        # Do first iteration.
        Δx_1 = pbs[1](Δ[1])

        # Allocate for cotangents.
        Δxs = get_adjoint_storage(x, 1, Δx_1[2])

        for n in 2:length(x)
            Δx = pbs[n](Δ[n])
            Δxs = _accum_at(Δxs, n, Δx[2])
        end

        return nothing, nothing, Δxs
    end

    return ys, zygote_friendly_map_pullback
end

zygote_friendly_map(f, x::Fill) = map(f, x)

function zygote_friendly_map(f, x::Base.Iterators.Zip{<:Tuple{Vararg{Fill, N}}}) where {N}
    return zygote_friendly_map(f, Fill(map(first, x.is), length(x)))
end


text:
# module for calculating particle motion
# particle has a map data for slam

using Distributions, LinearAlgebra, PDMats

include(joinpath(split(@__FILE__, "src")[1], "src/common/state_transition/state_transition.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/observation_function/observation_function.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/model/map/map.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/slam/estimated_object.jl"))

mutable struct MapParticle
  pose
  weight
  object_num
  map

  # init
  function MapParticle(init_pose, weight, object_num)
    self = new()
    self.pose = init_pose
    self.weight = weight
    self.object_num = object_num
    self.map = Map()
    for i in 1:object_num
      add_object(self.map, EstimatedObject(0.0, 0.0, id=i))
    end
    return self
  end
end

function motion_update(self::MapParticle, speed, yaw_rate,
                       time_interval, noise_rate_pdf)
  ns = rand(noise_rate_pdf) # [nn, no, on, oo]
  noised_spd = speed + ns[1]*sqrt(abs(speed)/time_interval) + ns[2]*sqrt(abs(yaw_rate)/time_interval)
  noised_yr = yaw_rate + ns[3]*sqrt(abs(speed)/time_interval) + ns[4]*sqrt(abs(yaw_rate)/time_interval)
  self.pose = state_transition(noised_spd, noised_yr, time_interval, self.pose)
end

function mat_H_m(mu_pose, obj_pose)
  obj_x, obj_y = obj_pose[1], obj_pose[2]
  mu_x, mu_y = mu_pose[1], mu_pose[2]
  mu_l = sqrt((mu_x - obj_x)^2 + (mu_y - obj_y)^2)
  return [(obj_x - mu_x)/mu_l (obj_y - mu_y)/mu_l;
          (mu_y - obj_y)/(mu_l^2) (obj_x - mu_x)/(mu_l^2)]
end

function mat_H_x(mu_pose, obj_pose)
  obj_x, obj_y = obj_pose[1], obj_pose[2]
  mu_x, mu_y = mu_pose[1], mu_pose[2]
  mu_l = sqrt((mu_x - obj_x)^2 + (mu_y - obj_y)^2)
  return [(mu_x - obj_x)/mu_l (mu_y - obj_y)/mu_l 0.0;
          (obj_y - mu_y)/(mu_l^2) (mu_x - obj_x)/(mu_l^2) -1.0]
end

function mat_Q(dist_dev, dir_dev)
  return [dist_dev^2 0.0;
          0.0 dir_dev^2]
end

function mat_M(speed, yaw_rate, time, stds)
  return diagm(0 => [stds["nn"]^2*abs(speed)/time + stds["no"]^2*abs(yaw_rate)/time,
                     stds["on"]^2*abs(speed)/time + stds["oo"]^2*abs(yaw_rate)/time])
end

function mat_A(speed, yaw_rate, time, theta)
  st, ct = sin(theta), cos(theta)
  stw, ctw = sin(theta + yaw_rate * time), cos(theta + yaw_rate * time)
  return [(stw - st)/yaw_rate  -speed/(yaw_rate^2)*(stw - st) + speed/yaw_rate*time*ctw;
          (-ctw + ct)/yaw_rate -speed/(yaw_rate^2)*(-ctw + ct) + speed/yaw_rate*time*stw;
          0                    time]
end

function params_for_drawing(pred_pose, landmark, dist_dev, dir_dev)
  # linearize observation function
  diff_x = pred_pose[1] - landmark.pose[1]
  diff_y = pred_pose[2] - landmark.pose[2]
  dist = hypot(diff_x, diff_y)
  Q = mat_Q(dist_dev * dist, dir_dev)
  pred_z = observation_function(pred_pose, landmark.pose)
  H_m = mat_H_m(pred_pose, landmark.pose) # 2x2
  H_x = mat_H_x(pred_pose, landmark.pose) # 2x3

  # calculate covariance matrix of sensor data
  # based on particle's pose and map
  Q_z = H_m * landmark.cov * H_m' + Q

  return pred_z, Q_z, H_x
end

function gauss_for_drawing(pred_pose, pred_cov, z,
                           landmark, dist_dev, dir_dev)
  pred_z, Q_z, H_x = params_for_drawing(pred_pose, landmark,
                                        dist_dev, dir_dev)
  K = pred_cov * H_x' * inv(Q_z + H_x*pred_cov*H_x')
  est_pose = K * (z - pred_z) + pred_pose
  est_cov = (Matrix{Float64}(I, 3, 3) - K*H_x) * pred_cov
  return est_pose, est_cov
end

function motion_update_2(self::MapParticle, speed, yaw_rate,
                         time_interval, motion_noise_stds,
                         observation, dist_dev, dir_dev)
  # distribution after movement
  M = mat_M(speed, yaw_rate, time_interval, motion_noise_stds)
  A = mat_A(speed, yaw_rate, time_interval, self.pose[3])
  est_cov = A * M * A'
  est_pose = state_transition(speed, yaw_rate, time_interval, self.pose)

  # update weight of particle
  for obs in observation
    pred_z, Q_z, H_x = params_for_drawing(est_pose,
                                          self.map.objects[obs[2]],
                                          dist_dev, dir_dev)
    sigma_z = H_x*est_cov*H_x' + Q_z
    self.weight *= pdf(MvNormal(pred_z, Symmetric(sigma_z)), obs[1])
  end

  # update pose and covariance of particle
  for obs in observation
    est_pose, est_cov = gauss_for_drawing(est_pose, est_cov, obs[1],
                                          self.map.objects[obs[2]],
                                          dist_dev, dir_dev)
  end
  est_cov += Matrix{Float64}(I,3,3).*0.00001
  self.pose = rand(MvNormal(est_pose, Symmetric(est_cov)))
end

function init_landmark_estimation(self::MapParticle, landmark, obs_pose,
                                  dist_dev_rate, dir_dev)
  lx = obs_pose[1] * cos(self.pose[3] + obs_pose[2]) + self.pose[1]
  ly = obs_pose[1] * sin(self.pose[3] + obs_pose[2]) + self.pose[2]
  landmark.pose = [lx, ly]
  H = mat_H_m(self.pose, landmark.pose)
  Q = mat_Q(dist_dev_rate * obs_pose[1], dir_dev)
  landmark.cov = inv(H' * inv(Q) * H)
end

function observation_update_landmark(self::MapParticle, landmark, obs_pose,
                                     dist_dev_rate, dir_dev, is_fastslam_1)
  est_obs_pose = observation_function(self.pose, landmark.pose)
  # not calculate when distance is too close
  if est_obs_pose[1] > 0.01
    # calculate kalman gain
    H = mat_H_m(self.pose, landmark.pose)
    Q = mat_Q(dist_dev_rate * est_obs_pose[1], dir_dev)
    K = landmark.cov * H' * inv(Q + H*landmark.cov*H')

    # update weight of particle
    # this process is for only fastslam 1.0
    if is_fastslam_1 == true
      Q_z = H * landmark.cov * H' + Q
      self.weight *= pdf(MvNormal(est_obs_pose, Symmetric(Q_z)), obs_pose)
    end

    # update landmark estimation
    landmark.pose += K * (obs_pose - est_obs_pose)
    landmark.cov = (Matrix{Float64}(I, 2, 2) - K*H) * landmark.cov
  end
end

function observation_update(self::MapParticle, observation,
                            dist_dev_rate, dir_dev, is_fastslam_1)
  for obs in observation
    obs_pose = obs[1] # [distance, direction]
    obs_id = obs[2]
    landmark = self.map.objects[obs_id]

    if landmark.cov === nothing
      init_landmark_estimation(self, landmark, obs_pose,
                               dist_dev_rate, dir_dev)
    else
      observation_update_landmark(self, landmark, obs_pose,
                                  dist_dev_rate, dir_dev, is_fastslam_1)
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# This does ML estimation of Negative Binomial model using the meps1996.data

# The MEPS data

# Define dep and expl vbls
#
# 	The dep. vbls, in corresponding column
# 	Office based doctor visits	1
# 	Outpatient doctor visits	2
# 	Emergency room visits		3
# 	Inpatient visits		4
# 	Dental visits			5
# 	Prescriptions			6
#

## load the data
using Econometrics, DelimitedFiles
which_dep = 1
if (which_dep == 1) dep = "OBDV"   end
if (which_dep == 2) dep = "OPV"    end
if (which_dep == 3) dep = "IPV"    end
if (which_dep == 4) dep = "ERV"    end
if (which_dep == 5) dep = "DV"     end
if (which_dep == 6) dep = "PRESCR" end
cd(@__DIR__)
data = readdlm("../Data/meps1996.data")

## create the regressors and dep var
y = data[:,which_dep]
x = data[:,7:12]
n = size(x,1)
x = [ones(n,1) x]

## to check effects of poor scaling, try commenting the next line
x[:,end] = x[:,end]/1000.0 # put income in thousands

## do ML
nbtype = 2
names = ["constant", "pub. ins.","priv. ins.", "sex", "age","edu","inc", "log(α)"]
title = "Negative Binomial-$nbtype model, "*dep*",  MEPS 1996 full data set"
model = θ  -> poisson(θ, y, x)
βhat, junk = mle(model, zeros(size(x,2)))
model = θ  -> negbin(θ, y, x, nbtype)
θstart = [βhat; 0.] # start values for estimation
# try adding the option vc=2 for "Hessian" or vc=3 for OPG
θhat, objvalue, V, converged = mleresults(model, θstart, title, names, vc=1);

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
# test typed version
using Plots, Revise, VortexDistributions

# include("typeversion.jl")
Exact()

# unitary?
v0 = [.2 .4 1]
w0 = PointVortex(v0)
@test rawData(w0) == v0

v1 = [.2 .4 1;0.7 1.5 -1;-.3 1.2 1]
w1 = PointVortex(v1)
@test rawData(w1) == v1

w3 = randPointVortex(1000)
v3 = rawData(w3)
@test v3 == rawData(w3)

# fast... really, really fast!
@time v3 = rawData(w3)
@time w3 = PointVortex(v3)

# Single vortex tests are reliable and don't require sorting
# (Phase cross talk will cause construction to wander)

@test foundNear(1)
@test foundNear(10)
@test foundNear(500)

# more tests below?
psi,vort = randVortexField(30); heatmap(psi.x,psi.y,abs2.(psi.ψ))

# Real data with masking examples
@load "./test/one_frame.jld2"
x0,y0 = 70,200

# circular
msk(x,y) = sqrt(x^2 + y^2) > 100
@. ψ1[msk(x-x0,y'-y0) ] = false
heatmap(x,y,abs2.(ψ1),transpose=true)
psi3 = Torus(ψ1,x,y)
vort3 = findvortices(psi3)
ψ1

# square
Lx,Ly = 100,100

xm = x[@. abs(x-x0)<Lx/2]
ym = y[@. abs(y-y0)<Ly/2]
psim = @. ψ1[abs(x-x0)<Lx/2,abs(y-y0)<Ly/2]

heatmap(xm,ym,angle.(psim),transpose=true)
psi2 = Torus(psim,xm,ym)
vort = findvortices(psi2)

scatter!([vort[1].xv],[vort[1].yv],alpha=0.4)
scatter!([vort[2].xv],[vort[2].yv],alpha=0.4)


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module RangeTest
using ParallelAccelerator

#ParallelAccelerator.DomainIR.set_debug_level(4)
#ParallelAccelerator.ParallelIR.set_debug_level(4)
#ParallelAccelerator.CGen.set_debug_level(4)
#ParallelAccelerator.set_debug_level(4)

@acc function reduce_col(col::Int)
    A = rand(10^5,10);
    x = rand(10^5);
    sum((A[:, col] .- x[:]) .* (A[:, col] .- x[:]))
end

@acc function reduce_col2(col::Int)
    A = rand(10^5,10);
    x = rand(10^5);
    r = 1:size(A, 1)
    sum((A[r, col] .- x[r]) .* (A[r, col] .- x[r]))
end

function test1()
    return reduce_col(3)
end

function test2()
    return reduce_col2(3)
end

end

using Base.Test

println("Testing ranges...")

@test RangeTest.test1() > 1.66e3
@test RangeTest.test2() > 1.66e3

println("Done testing ranges.")


text:
function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this        = parse(Int64, this)

#
# Important directories
simdir          = "/remote/ceph2/group/gerda/data/mpik/gerda-simulations/gerda-gems-db/pss-dev"
simdir_old      = "/remote/ceph/group/gerda/data/simulation/gerda-mage-sim"
mapping_file    = simdir_old*"/UTILS/det-data/ged-mapping.json"
parameters_file = simdir_old*"/UTILS/det-data/ged-parameters.json"
config_dir      = "../../../2020-02-06_8380701d_st_ffischer/res/Impurity_Scan/config-dep/"

#
# Load packages and functions
include("../src/init.jl")
include("../src/fct.jl")
include("../src/worker_fct.jl")

#
# Output paths and filter settings
plots_base_path = "../../waveforms/sim/plots/raw/"
base_path_raw   = "../../waveforms/sim/raw/"
base_path       = "../../waveforms/sim/raw_wf/"
n_sim_events    = 10000;
sampling_time   = 1u"ns";
generate_wf     = false;
generate_cl_wf  = true;
no_new_files    = 0

channels = []
for c in 0:1:36
    files = glob(base_path_raw * "raw*/" * lpad(c, 2, "0") * "-" * channel_to_name[c] * ".h5");
    if length(files) > 0
        push!(channels, c)
    end
end
channels = get_share_for_worker(channels, num_workers, this)

while no_new_files <= 3
    for ch in channels
        @info("Start Ch"*lpad(ch, 2, "0")*" | "*channel_to_name[ch])
        @info(">---------------------------------------------------------------<")
        files = glob(base_path_raw*"*/*-"*channel_to_name[ch]*".h5")

        # Read the SSD simulation file of depleted detector
        config_file = glob(config_dir * channel_to_name[ch] * "*.config")[1]
        sim_output  = glob(config_dir * "../output/"*channel_to_name[ch]*"*/" )[1]
        @info("Load SSD simulation")
        @time sim = @suppress readcreateh5(config_file, sim_output);

        for file in files
            output_dir = base_path * lpad(ch, 2, "0") * "-" * channel_to_name[ch] * "/cl-wf/"
            output_dir *= split(file, "/")[end-1] * "/"
            @time events = h5open(file, "r") do h5f
                LegendHDF5IO.readdata(h5f, "data")
            end;
            last_file = glob(output_dir * "*.h5")
            if length(last_file) == 0
                existance_check = false
            else
                last_file       = last_file[end]
                last_event      = parse(Int64,split(split(last_file,"-")[end], ".h5")[1])
                existance_check = isfile(last_file) && last_event == size(events,1) && stat(last_file).size > 0
            end
            if existance_check
                @info("This file has already been simluated")
            else
                no_new_files = 0
                isdir(output_dir) ? "Directory exists" : output_dir
                file_id = findall(x->x == file, files)[1]

                if generate_wf
                    # IJulia.clear_output(true)
                    Base.run(`clear`)
                    @info("File $file_id | " * string(length(files)))
                    @info("Start simulating waveforms for ch" * lpad(ch, 2, "0") * " | " * channel_to_name[ch])
                    @info(">---------------------------------------------------------------<")
                    output_basename  = lpad(ch, 2, "0") * "-" * channel_to_name[ch]*"-wf-"*split(file, "/")[end-1]
                    t_with_waveforms = @suppress SSD.simulate_waveforms(events, sim, output_dir, output_basename, chunk_n_physics_events = n_sim_events, Δt = sampling_time);
                end

                if generate_cl_wf
                    # IJulia.clear_output(true)
                    Base.run(`clear`)
                    @info("File $file_id | "*string(length(files)))
                    @info("Start simulating clustered waveforms for ch" * lpad(ch, 2, "0") * " | " * channel_to_name[ch])
                    @info(">---------------------------------------------------------------<")
                    @info("$(sum(length.(events.edep))) hits before clustering")
                    events_clustered = SSD.cluster_detector_hits(events, 0.2u"mm")
                    @info("$(sum(length.(events_clustered.edep))) hits after clustering")
                    output_basename            = lpad(ch, 2, "0") * "-" * channel_to_name[ch] * "-cl-wf-"*split(file, "/")[end-1]
                    t_clustered_with_waveforms = SSD.simulate_waveforms(events_clustered, sim, output_dir, output_basename, chunk_n_physics_events = n_sim_events, Δt = sampling_time);
                end
            end
        end
    end
    no_new_files += 1
    sleep(300)
end


text:
# used by OTS models
function check_br_status(sol)
    for (i,branch) in sol["branch"]
        @test branch["br_status"] == 0.0 || branch["br_status"] == 1.0
    end
end

if (Pkg.installed("AmplNLWriter") != nothing && Pkg.installed("CoinOptServices") != nothing)
    @testset "test ac ots" begin
        #Omitting this test, until bugs can be resolved, bonmin does not report a integral solution
        #@testset "3-bus case" begin
        #    result = run_ots("../test/data/case3.m", ACPPowerModel, BonminNLSolver(["bonmin.bb_log_level=0", "bonmin.nlp_log_level=0"]))

        #    check_br_status(result["solution"])

        #    @test result["status"] == :LocalOptimal
        #    @test isapprox(result["objective"], 5812; atol = 1e0)
        #end
        @testset "5-bus case" begin
            result = run_ots("../test/data/case5.m", ACPPowerModel, BonminNLSolver(["bonmin.bb_log_level=0", "bonmin.nlp_log_level=0"]))

            check_br_status(result["solution"])

            @test result["status"] == :LocalOptimal
            # NOTE this objective value is out of date, and this test will not pass
            @test isapprox(result["objective"], 15174; atol = 1e0)
        end
        #Omitting this test, returns local infeasible
        #@testset "6-bus case" begin
        #    result = run_ots("../test/data/case6.m", ACPPowerModel, BonminNLSolver(["bonmin.bb_log_level=0", "bonmin.nlp_log_level=0"]))

        #    check_br_status(result["solution"])

        #    @test result["status"] == :LocalOptimal
        #    println(result["objective"])
        #    @test isapprox(result["objective"], 15174; atol = 1e0)
        #end
    end
end

@testset "test dc ots" begin
    @testset "3-bus case" begin
        result = run_ots("../test/data/case3.m", DCPPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 5782.0; atol = 1e0)
    end
    @testset "5-bus case" begin
        result = run_ots("../test/data/case5.m", DCPPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 14991.2; atol = 1e0)
    end
    @testset "6-bus case" begin
        result = run_ots("../test/data/case6.m", DCPPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 11396.3; atol = 1e0)
    end
end

@testset "test dc+ll ots" begin
    @testset "3-bus case" begin
        result = run_ots("../test/data/case3.m", DCPLLPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 5885.2; atol = 1e0)
    end
    @testset "5-bus case" begin
        result = run_ots("../test/data/case5.m", DCPLLPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 15275.2; atol = 1e0)
    end
    @testset "6-bus case" begin
        result = run_ots("../test/data/case6.m", DCPLLPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 11515.6; atol = 1e0)
    end
end

@testset "test soc ots" begin
    @testset "3-bus case" begin
        result = run_ots("../test/data/case3.m", SOCWRPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 5746.7; atol = 1e0)
    end
    @testset "5-bus case" begin
        result = run_ots("../test/data/case5.m", SOCWRPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 15051.4; atol = 1e0)
    end
    @testset "6-bus case" begin
        result = run_ots("../test/data/case6.m", SOCWRPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 11559.8; atol = 1e0)
    end
end

@testset "test qc ots" begin
    @testset "3-bus case" begin
        result = run_ots("../test/data/case3.m", QCWRPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 5746.7; atol = 1e0)
    end
    @testset "5-bus case" begin
        result = run_ots("../test/data/case5.m", QCWRPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 15051.4; atol = 1e0)
    end
    @testset "5-bus asymmetric case" begin
        result = run_ots("../test/data/case5_asym.m", QCWRPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 14999.7; atol = 1e0)
    end
    @testset "6-bus case" begin
        result = run_ots("../test/data/case6.m", QCWRPowerModel, pajarito_solver)

        check_br_status(result["solution"])

        @test result["status"] == :Optimal
        @test isapprox(result["objective"], 11567.1; atol = 1e0)
    end
end


text:
push!(LOAD_PATH, "src")
using Documenter, SubmodularUtilities

makedocs(sitename="SubmodularUtilities.jl", doctest = false,
    pages = ["Home" => "index.md"]
)
# makedocs(sitename="SubmodularUtilities.jl",
# 	 repo="https://github.com/lchen91/Submodular_Utilities/blob/master{path}#{line}")

# deploydocs(
# 	repo = "github.com/lchen91/Submodular_Utilities.git",
# )

deploydocs(
    repo   = "github.com/lchen91/SubmodularUtilities.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
)


text:
# The jackknife is a resampling technique especially useful for variance and bias estimation.
"""
    jackknife_avg(`scores`)

Apply jackknife on the input list of `scores`
"""
function jackknife_avg(scores)
    if length(collect(Set(scores))) == 1
        #= In case the elements of the array are all equal=#
        return scores[1]
    else
        #=store the maximum scores
        from the m different sets of m-1 scores.
        such that m is the len(score_list)=#
        average = []
        for i in scores
            # dummy : list a particular combo of m-1 scores
            dummy = [j for j in scores if i != j]
            append!(average, max(dummy...))
            end

        return sum(average)/length(average)
    end
end

"""
    weighted_lcs(X, Y, weight_score::Bool, returns_string::Bool, weigthing_function::Function)

Compute the Weighted Longest Common Subsequence of X and Y.
"""
function weighted_lcs(X, Y, weighted = true, returns_string = false, f = sqrt)
    m, n = length(X), length(Y)
    c_table = [zeros(n+1) for i in 1:m+1]
    w_table = [zeros(n+1) for i in 1:m+1]

    for i in 2:(m+1)
        for j in 2:(n+1)
            if X[i-1] == Y[j-1]
                k = w_table[i-1][j-1]
                if weighted == true
                    increment = (f(k+1)) - (f(k))
                else
                    increment = 1
                end
                c_table[i][j] = c_table[i-1][j-1] + increment
                w_table[i][j] = k + 1
            else
                if c_table[i-1][j] > c_table[i][j-1]
                    c_table[i][j] = c_table[i-1][j]
                    w_table[i][j] = 0  # no match at i,j
                else
                    c_table[i][j] = c_table[i][j-1]
                    w_table[i][j] = 0  # no match at i,j
                end
            end
        end
    end

    lcs_length = (c_table[m+1][n+1])
    if returns_string == false
        return lcs_length
    end

    if weighted == true
        lcs_length = c_table[m][n]^(2)
    end

    lcs_length = round(Int64, lcs_length)
    lcs_length = convert(Int64, lcs_length)
    lcs = ["" for i in  1:(lcs_length+1)]
    lcs[lcs_length+1] = ""
    i = m + 1
    j = n + 1

    while i>1 && j>1
        if X[i-1] == Y[j-1]
            lcs[lcs_length+1] = X[i-1]
            i -= 1
            j -= 1
            lcs_length -= 1
        elseif c_table[i-1][j] > c_table[i][j-1]
            i -= 1
        else
            j -= 1
        end
    end

    return (join(lcs, " "))  # the lcs string
end

"""
    fmeasure_lcs(RLCS, PLCS, β)

Compute the F-measure based on WLCS.

# Arguments

- `RLCS` - Recall Factor
- `PLCS` - Precision Factor
- `β` - Parameter
"""
function fmeasure_lcs(RLCS, PLCS, β=1)
    try
        return ((1 + β ^ 2) * RLCS * PLCS) / (RLCS + (β ^ 2) * PLCS)
    catch ex
        if ex isa DivideError
            return 0
        else
            rethrow(ex)
        end
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    KrigingEstimator

A Kriging estimator (e.g. Simple Kriging).
"""
abstract type KrigingEstimator end

"""
    KrigingState(X, z, LHS, RHS)

A Kriging state stores information needed
to perform estimation at any given location.
"""
mutable struct KrigingState{T<:Real,V,
                            A<:AbstractMatrix{T},
                            B<:AbstractVector{V},
                            F<:Factorization,R}
  X::A
  z::B
  LHS::F
  RHS::Vector{R}
end

"""
    KrigingWeights(λ, ν)

An object storing Kriging weights `λ` and Lagrange multipliers `ν`.
"""
struct KrigingWeights{T<:Real,A<:AbstractVector{T}}
  λ::A
  ν::A
end

"""
    FittedKriging(estimator, state)

An object that can be used for making predictions using the
parameters in `estimator` and the current Kriging `state`.
"""
struct FittedKriging{E<:KrigingEstimator,S<:KrigingState}
  estimator::E
  state::S
end

"""
    status(fittedkrig)

Return the status of the `fittedkrig` object, meaning
the factorization of the Kriging system was successful.
"""
status(fittedkrig::FittedKriging) = issuccess(fittedkrig.state.LHS)

#--------------
# FITTING STEP
#--------------

"""
    fit(estimator, X, z)

Build Kriging system from coordinates `X` and
values `z` and return a fitted estimator.
"""
function fit(estimator::KrigingEstimator, X::AbstractMatrix, z::AbstractVector)
  # build Kriging system
  LHS = lhs(estimator, X)
  RHS = Vector{eltype(LHS)}(undef, size(LHS,1))

  # factorize LHS
  FLHS = factorize(estimator, LHS)

  # record Kriging state
  state = KrigingState(X, z, FLHS, RHS)

  # return fitted estimator
  FittedKriging(estimator, state)
end

"""
    lhs(estimator, X)

Return LHS of Kriging system using spatial configuration `X`.
"""
function lhs(estimator::KrigingEstimator, X::AbstractMatrix)
  γ = estimator.γ
  nobs = size(X, 2)
  ncons = nconstraints(estimator)

  # pre-allocate memory for LHS
  x = view(X,:,1)
  T = Variography.result_type(γ, x, x)
  m = nobs + ncons
  LHS = Matrix{T}(undef, m, m)

  # set variogram/covariance block
  pairwise!(LHS, γ, X)
  if isstationary(γ)
    for j=1:nobs, i=1:nobs
      @inbounds LHS[i,j] = sill(γ) - LHS[i,j]
    end
  end

  # set blocks of constraints
  set_constraints_lhs!(estimator, LHS, X)

  LHS
end

"""
    nconstraints(estimator)

Return number of constraints for `estimator`.
"""
nconstraints(estimator::KrigingEstimator) = error("not implemented")

"""
    set_constraints_lhs!(estimator, LHS, X)

Set constraints in LHS of Kriging system.
"""
set_constraints_lhs!(estimator::KrigingEstimator,
                     LHS::AbstractMatrix, X::AbstractMatrix) = error("not implemented")

"""
    factorize(estimator, LHS)

Factorize LHS of Kriging system with appropriate factorization method.
"""
factorize(estimator::KrigingEstimator, LHS::AbstractMatrix) = error("not implemented")

#-----------------
# PREDICTION STEP
#-----------------

"""
    predict(estimator, xₒ)

Compute mean and variance for the `estimator` at coordinates `xₒ`.
"""
predict(estimator::FittedKriging, xₒ::AbstractVector) =
  combine(estimator, weights(estimator, xₒ), estimator.state.z)

"""
    weights(estimator, xₒ)

Compute the weights λ (and Lagrange multipliers ν) for the
`estimator` at coordinates `xₒ`.
"""
function weights(estimator::FittedKriging, xₒ::AbstractVector)
  nobs = size(estimator.state.X, 2)

  set_rhs!(estimator, xₒ)

  # solve Kriging system
  x = estimator.state.LHS \ estimator.state.RHS

  λ = view(x,1:nobs)
  ν = view(x,nobs+1:length(x))

  KrigingWeights(λ, ν)
end

"""
    set_rhs!(estimator, xₒ)

Set RHS of Kriging system at coodinates `xₒ`.
"""
function set_rhs!(estimator::FittedKriging, xₒ::AbstractVector)
  γ = estimator.estimator.γ
  X = estimator.state.X
  RHS = estimator.state.RHS

  # RHS variogram/covariance
  @inbounds for j in 1:size(X, 2)
    xj = view(X,:,j)
    RHS[j] = isstationary(γ) ? sill(γ) - γ(xj, xₒ) : γ(xj, xₒ)
  end

  set_constraints_rhs!(estimator, xₒ)
end

"""
    set_constraints_rhs!(estimator, xₒ)

Set constraints in RHS of Kriging system.
"""
set_constraints_rhs!(estimator::FittedKriging, xₒ::AbstractVector) = error("not implemented")

"""
    combine(estimator, weights, z)

Combine `weights` with values `z` to produce mean and variance.
"""
function combine(estimator::FittedKriging, weights::KrigingWeights, z::AbstractVector)
  γ = estimator.estimator.γ
  b = estimator.state.RHS
  λ = weights.λ
  ν = weights.ν

  # compute b⋅[λ;ν]
  nobs  = length(λ)
  c₁ = view(b,1:nobs)⋅λ
  c₂ = view(b,nobs+1:length(b))⋅ν
  c = c₁ + c₂

  if isstationary(γ)
    z⋅λ, sill(γ) - c
  else
    z⋅λ, c
  end
end

#-----------------
# IMPLEMENTATIONS
#-----------------
include("estimators/simple_kriging.jl")
include("estimators/ordinary_kriging.jl")
include("estimators/universal_kriging.jl")
include("estimators/external_drift_kriging.jl")


text:
# ---
# title: 310. Minimum Height Trees
# id: problem310
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Breadth-first Search, Graph
# link: <https://leetcode.com/problems/minimum-height-trees/description/>
# hidden: true
# ---
#
# A tree is an undirected graph in which any two vertices are connected by
# _exactly_  one path. In other words, any connected graph without simple cycles
# is a tree.
#
# Given a tree of `n` nodes labelled from `0` to `n - 1`, and an array of `n -
# 1` `edges` where `edges[i] = [ai, bi]` indicates that there is an undirected
# edge between the two nodes `ai` and `bi` in the tree, you can choose any node
# of the tree as the root. When you select a node `x` as the root, the result
# tree has height `h`. Among all possible rooted trees, those with minimum
# height (i.e. `min(h)`)  are called **minimum height trees** (MHTs).
#
# Return _a list of all **MHTs '** root labels_. You can return the answer in
# **any order**.
#
# The **height** of a rooted tree is the number of edges on the longest downward
# path between the root and a leaf.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/01/e1.jpg)
#
#
#
#     Input: n = 4, edges = [[1,0],[1,2],[1,3]]
#     Output: [1]
#     Explanation: As shown, the height of the tree is 1 when the root is the node with label 1 which is the only MHT.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/01/e2.jpg)
#
#
#
#     Input: n = 6, edges = [[3,0],[3,1],[3,2],[3,4],[5,4]]
#     Output: [3,4]
#
#
# **Example 3:**
#
#
#
#     Input: n = 1, edges = []
#     Output: [0]
#
#
# **Example 4:**
#
#
#
#     Input: n = 2, edges = [[0,1]]
#     Output: [0,1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 2 * 104`
#   * `edges.length == n - 1`
#   * `0 <= ai, bi < n`
#   * `ai != bi`
#   * All the pairs `(ai, bi)` are distinct.
#   * The given input is **guaranteed** to be a tree and there will be **no repeated** edges.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
""" Parse Julia programs into morphisms represented as wiring diagrams.
"""
module ParseJuliaPrograms
export @program, parse_wiring_diagram

using Base: invokelatest
using Compat
using Match

using ...Catlab
import ...Meta: Expr0
using ...Doctrines: ObExpr, HomExpr
using ...WiringDiagrams
using ..GenerateJuliaPrograms: make_return_value

""" Parse a wiring diagram from a Julia program.

For the most part, this is standard Julia but we take a few liberties with the
syntax. Products are represented as tuples. So if `x` and `y` are variables of
type \$X\$ and \$Y\$, then `(x,y)` has type \$X \\otimes Y\$. Also, both `()`
and `nothing` are interpreted as the monoidal unit \$I\$.

Unlike in standard Julia, the call expressions `f(x,y)` and `f((x,y))` are
equivalent. Consequently, given morphisms \$f: W \\to X \\otimes Y\$ and \$g: X
\\otimes Y \\to Z\$, the code

```julia
x, y = f(w)
g(x,y)
```

is equivalent to `g(f(w))`. In standard Julia, at most one of these calls to `g`
would be valid.

The diagonals (copying and deleting) are implicit in the Julia syntax: copying
is variable reuse and deleting is variable non-use. For the codiagonals (merging
and creating), a special syntax is provided, reinterpreting Julia's vector
literals. The merge of `x1` and `x2` is represented by the vector `[x1,x2]` and
creation by the empty vector `[]`. For example, `f([x1,x2])` translates to
`compose(mmerge(X),f)`.

This macro is a wrapper around [`parse_wiring_diagram`](@ref).
"""
macro program(pres, expr)
  Expr(:call, GlobalRef(ParseJuliaPrograms, :parse_wiring_diagram),
       esc(pres), esc(Expr(:quote, expr)))
end
macro program(pres, call, body)
  Expr(:call, GlobalRef(ParseJuliaPrograms, :parse_wiring_diagram),
       esc(pres), esc(Expr(:quote, call)), esc(Expr(:quote, body)))
end

""" Parse a wiring diagram from a Julia function expression.

The macro version of this function is [`@program`](@ref).
"""
function parse_wiring_diagram(pres::Presentation, expr::Expr)::WiringDiagram
  @match expr begin
    Expr(:function, [call, body]) => parse_wiring_diagram(pres, call, body)
    Expr(:->, [call, body]) => parse_wiring_diagram(pres, call, body)
    _ => error("Not a function or lambda expression")
  end
end

function parse_wiring_diagram(pres::Presentation, call::Expr0, body::Expr)::WiringDiagram
  # FIXME: Presentations should be uniquely associated with syntax systems.
  syntax_module = Syntax.syntax_module(first(pres.generators))

  # Parse argument names and types from call expression.
  call_args = @match call begin
    Expr(:call, [name, args...]) => args
    Expr(:tuple, args) => args
    Expr(:(::), _) => [call]
    _::Symbol => [call]
    _ => error("Invalid function signature: $call")
  end
  parsed_args = map(call_args) do arg
    @match arg begin
      Expr(:(::), [name::Symbol, type_expr::Expr0]) =>
        (name, eval_type_expr(pres, syntax_module, type_expr))
      _ => error("Argument $arg is missing name or type")
    end
  end

  # Compile...
  args = first.(parsed_args)
  kwargs = make_lookup_table(pres, syntax_module, unique_symbols(body))
  func_expr = compile_recording_expr(body, args,
    kwargs = sort!(collect(keys(kwargs))))
  func = parentmodule(syntax_module).eval(func_expr)

  # ...and then evaluate function that records the function calls.
  arg_obs = syntax_module.Ob[ last(arg) for arg in parsed_args ]
  arg_blocks = Int[ length(to_wiring_diagram(ob)) for ob in arg_obs ]
  inputs = to_wiring_diagram(otimes(arg_obs))
  diagram = WiringDiagram(inputs, munit(typeof(inputs)))
  v_in, v_out = input_id(diagram), output_id(diagram)
  arg_ports = [ Tuple(Port(v_in, OutputPort, i) for i in (stop-len+1):stop)
                for (len, stop) in zip(arg_blocks, cumsum(arg_blocks)) ]
  recorder = f -> (args...) -> record_call!(diagram, f, args...)
  value = invokelatest(func, recorder, arg_ports...; kwargs...)

  # Add outgoing wires for return values.
  out_ports = normalize_arguments((value,))
  diagram.output_ports = [
    # XXX: Inferring the output port types is not reliable.
    port_value(diagram, first(ports)) for ports in out_ports
  ]
  add_wires!(diagram, [
    port => Port(v_out, InputPort, i)
    for (i, ports) in enumerate(out_ports) for port in ports
  ])
  substitute(diagram)
end

""" Make a lookup table assigning names to generators or term constructors.
"""
function make_lookup_table(pres::Presentation, syntax_module::Module, names)
  theory = syntax_module.theory().class().theory
  terms = Set([ term.name for term in theory.terms ])

  table = Dict{Symbol,Any}()
  for name in names
    if has_generator(pres, name)
      table[name] = generator(pres, name)
    elseif name in terms
      table[name] = (args...) -> invoke_term(syntax_module, name, args...)
    end
  end
  table
end

""" Evaluate pseudo-Julia type expression, such as `X` or `otimes{X,Y}`.
"""
function eval_type_expr(pres::Presentation, syntax_module::Module, expr::Expr0)
  function eval(expr)
    @match expr begin
      Expr(:curly, [name, args...]) =>
        invoke_term(syntax_module, name, map(eval, args)...)
      name::Symbol => generator(pres, name)
      _ => error("Invalid type expression $expr")
    end
  end
  eval(expr)
end

""" Generate a Julia function expression that will record function calls.

Rewrites the function body so that:

  1. Ordinary function calls are mapped to recorded calls, e.g.,
     `f(x,y)` becomes `recorder(f,x,y)`
  2. "Curly" function calls are mapped to ordinary function calls, e.g.,
     `f{X,Y}` becomes `f(X,Y)`
"""
function compile_recording_expr(body::Expr, args::Vector{Symbol};
    kwargs::Vector{Symbol}=Symbol[],
    recorder::Symbol=Symbol("##recorder"))::Expr
  function rewrite(expr)
    @match expr begin
      Expr(:call, [f, args...]) =>
        Expr(:call, Expr(:call, recorder, rewrite(f)), map(rewrite, args)...)
      Expr(:curly, [f, args...]) =>
        Expr(:call, rewrite(f), map(rewrite, args)...)
      Expr(head, args) => Expr(head, map(rewrite, args)...)
      _ => expr
    end
  end
  Expr(:function,
    Expr(:tuple,
      Expr(:parameters, (Expr(:kw, kw, nothing) for kw in kwargs)...),
      recorder, args...),
    rewrite(body))
end

""" Record a Julia function call as a box in a wiring diagram.
"""
function record_call!(diagram::WiringDiagram, f::HomExpr, args...)
  # Add a new box, itself a wiring diagram, for the call.
  subdiagram = to_wiring_diagram(f)
  v = add_box!(diagram, subdiagram)

  # Adding incoming wires.
  inputs = input_ports(subdiagram)
  arg_ports = normalize_arguments(Tuple(args))
  @assert length(arg_ports) == length(inputs)
  add_wires!(diagram, [
    Wire(port => Port(v, InputPort, i))
    for (i, ports) in enumerate(arg_ports) for port in ports
  ])

  # Return output ports.
  outputs = output_ports(subdiagram)
  return_ports = [ Port(v, OutputPort, i) for i in eachindex(outputs) ]
  make_return_value(return_ports)
end

""" Normalize arguments given as (possibly nested) tuples or vectors of values.
"""
function normalize_arguments(xs::Tuple)
  mapreduce(normalize_arguments, (xs,ys) -> (xs..., ys...), xs; init=())
end
function normalize_arguments(xs::Vector)
  xss = map(normalize_arguments, flatten_vec(xs)) # Vector of lists of vectors
  if isempty(xss)
    ([],) # Degenerate case
  else
    Tuple(reduce(vcat, xs) for xs in zip(xss...))
  end
end
normalize_arguments(::Nothing) = ()
normalize_arguments(x) = ([x],)
flatten_vec(xs::Vector) = mapreduce(flatten_vec, vcat, xs; init=[])
flatten_vec(x) = [x]

""" Set of all symbols occuring in a Julia expression.
"""
unique_symbols(expr::Expr) =
  reduce(union!, map(unique_symbols, expr.args); init=Set{Symbol}())
unique_symbols(x::Symbol) = Set([x])
unique_symbols(x) = Set{Symbol}()

end


text:
module ISP
include(joinpath(dirname(@__FILE__), "globSettings.jl"))

# Include common types
include(joinpath(dirname(@__FILE__), "ISPPDEtypes.jl"))

# Include functionals we have defined so far.
include(joinpath(dirname(@__FILE__), "Functionals.jl"))

end

max_stars_repo_name: jgoldfar/isp-julia-mwe

max_issues_repo_name: jgoldfar/isp-julia-mwe

max_forks_repo_name: jgoldfar/isp-julia-mwe


text:
include("TriPMatrices.jl")
n = 6
D = [linspace(1,n,n);] + 5.0
Dl = [linspace(1,n,n);] + 1.0
Du = [linspace(1,n,n);] + 2.0

M = TriPMatrices.TridiagonalP(Dl, D, Du)
Mf = full(M)

S = TriPMatrices.SymTridiagonalP(D, Du)
Sf = full(S)