text:
using ArgParse
using DataFrames
using DataFramesMeta
using ProgressMeter
using Statistics

import Colors
import CSV
import Dates
import LibGit2
import Pkg
import Pkg.TOML
import Plots

## Common

parse_toml_config(config::T where T <: AbstractString) =
    parse_toml_config(TOML.parsefile(config))

function get_default_config()
    return deepcopy(Dict{String, Any}(
        "Data" => Dict{String, Any}(
            "x-column" => "x",
            "y-column" => "y",
            "gene-column" => "gene",
            "min-molecules-per-gene" => 1,
            "min-molecules-per-cell" => 3,
            "estimate-scale-from-centers" => true,
            "scale" => nothing,
            "scale-std" => "25%",
            "min-molecules-per-segment" => nothing
        ),
        "Sampling" => Dict{String, Any}(
            "new-component-weight" => 0.2,
            "new-component-fraction" => 0.3
        ),
        "Plotting" => Dict{String, Any}(
            "gene-composition-neigborhood" => nothing,
            "min-pixels-per-cell" => 15
        )
    ))
end

function parse_toml_config(config::Dict{AS, Any}) where AS <: AbstractString
    res_config = get_default_config()
    for (k,v) in config
        if !(k in keys(res_config))
            error("Unexpected value in the config: '$k'")
        end

        cur_def = res_config[k]

        for (k2,v2) in v
            if !(k2 in keys(cur_def))
                error("Unexpected value in the config: '$k' -> '$k2'")
            end

            cur_def[k2] = v2
        end
    end

    return res_config
end

function extend_params_with_config!(params::Dict, config::Dict)
    for sub_cfg in values(config)
        for (k, v) in sub_cfg
            if !(k in keys(params)) || params[k] === nothing
                params[k] = v
            end
        end
    end
end

load_df(args::Dict; kwargs...) = load_df(args["coordinates"]; x_col=args["x-column"], y_col=args["y-column"], gene_col=args["gene-column"], min_molecules_per_gene=args["min-molecules-per-gene"], kwargs...)

append_suffix(output::String, suffix) = "$(splitext(output)[1])_$suffix"

## Main run

function parse_commandline(args::Union{Nothing, Array{String, 1}}=nothing) # TODO: add verbosity level
    s = ArgParseSettings(prog="baysor run")
    @add_arg_table! s begin
        "--config", "-c"
            help = "TOML file with config"
        "--x-column", "-x"
            help = "Name of x column. Overrides the config value."
        "--y-column", "-y"
            help = "Name of gene column. Overrides the config value."
        "--gene-column", "-g"
            help = "Name of gene column. Overrides the config value."

        "--iters", "-i"
            help = "Number of iterations"
            arg_type = Int
            default = 500
        "--min-molecules-per-cell", "-m"
            help = "Minimal number of molecules for a cell to be considered as real. It's an important parameter, as it's used to infer several other parameters. Overrides the config value."
            arg_type = Int
        "--n-clusters"
            help = "Number of molecule clusters, i.e. major cell types. Depends on protocol resolution, but should not be too high. In most cases something between 3 and 15 should work well."
            arg_type = Int
            default=4
        "--n-frames", "-n"
            help = "Number of frames, which is the same as number of processes. Algorithm data is splitted by frames to allow parallel run over frames."
            arg_type = Int
            default=1
        "--num-cells-init"
            help = "Initial number of cells."
            arg_type = Int
        "--output", "-o"
            help = "Name of the output file or path to the output directory"
            default = "segmentation.csv"
        "--plot", "-p"
            help = "Save pdf with plot of the segmentation"
            action = :store_true
        "--save-polygons"
            help = "Save estimated cell boundary polygons to a file with a specified FORMAT. Only 'GeoJSON' format is currently supported. The option requires setting '-p' to work."
            arg_type = String
            range_tester = (x -> in(lowercase(x), ["geojson"]))
            metavar = "FORMAT"

        "--scale", "-s"
            help = "Scale parameter, which suggest approximate cell radius for the algorithm. Must be in the same units as 'x' and 'y' molecule coordinates. Overrides the config value. Sets 'estimate-scale-from-centers' to false."
            arg_type = Float64
        "--prior-segmentation-confidence"
            help = "Confidence of the `prior_segmentation` results. Value in [0; 1]. If you want the final segmentation not contradicting to prior_segmentation, set it to 1. Otherwise, if you assume errors in prior_segmentation, values in [0.2-0.7] allow flexibility for the algorithm."
            arg_type = Float64
            default = 0.2

        "coordinates"
            help = "CSV file with coordinates of molecules and gene type"
            required = true
        "prior_segmentation"
            help = "Image or a MAT file with segmentation mask (either boolean or component indexing) or CSV column with integer segmentation labels. If it's the column name, it should be preceded ':' symbol (e.g. :cell)"
    end

    return (args === nothing) ? parse_args(s) : parse_args(args, s)
end

function parse_configs(args::Union{Nothing, Array{String, 1}}=nothing)
    r = parse_commandline(args)
    if r["scale"] !== nothing
        r["estimate-scale-from-centers"] = false
    end

    if r["config"] !== nothing
        extend_params_with_config!(r, parse_toml_config(r["config"]))
    else
        @warn "No config file provided. Back-up to default parameters."
        extend_params_with_config!(r, get_default_config())
    end

    for k in ["gene-column", "x-column", "y-column"]
        if !(k in keys(r)) || (r[k] === nothing) # should never be the case as we have defaults
            error("$k must be specified")
        end
        r[k] = Symbol(r[k])
    end

    if r["prior_segmentation"] === nothing && r["scale"] === nothing
        println("Either `prior_segmentation` or `scale` must be provided.")
        exit(1)
    end

    if r["min-molecules-per-segment"] === nothing
        r["min-molecules-per-segment"] = default_param_value(:min_molecules_per_segment, r["min-molecules-per-cell"])
    end

    if isdir(r["output"]) || isdirpath(r["output"])
        r["output"] = joinpath(r["output"], "segmentation.csv")
    end

    if (r["save-polygons"] !== nothing) && (!r["plot"])
        @warn "--plot option is required for saving polygons (--save-polygons). The polygons will not be saved."
    end

    return r
end

function plot_diagnostics_panel(df_res::DataFrame, assignment::Array{Int, 1}, tracer::Dict, args::Dict; margin=5*Plots.mm,
        max_diffs::Union{Vector{Float64}, Nothing}=nothing, change_fracs::Union{Vector{Float64}, Nothing}=nothing)
    @info "Plot diagnostics"
    open(append_suffix(args["output"], "diagnostics.html"), "w") do io
        println(io, "")
        # Molecule clustering convergence
        if max_diffs !== nothing
            p_mol_conv = Plots.plot(max_diffs, xlabel="Iteration", ylabel="Maximal probability change", title="Molcecule clustering convergence", label="Maximal change");
            if change_fracs !== nothing
                p_mol_conv = Plots.plot!(change_fracs, label="Fraction of molecules changed")
            end
            show(io, MIME("text/html"), p_mol_conv)
        end

        # Main algorithm convergence
        if (:n_components in keys(tracer)) && length(tracer[:n_components]) != 0
            p_conv = plot_num_of_cells_per_iterarion(tracer);
            show(io, MIME("text/html"), p_conv)
        end

        println(io, "

")

        # Confidence per molecule
        if :confidence in propertynames(df_res)
            bins = 0.0:0.025:1.0
            p_conf = Plots.histogram(df_res.confidence[assignment .!= 0], bins=bins, label="Assigned molecules",
                xlabel="Confidence", ylabel="#Molecules", title="Confidence per molecule", margin=margin)
            p_conf = Plots.histogram!(df_res.confidence[assignment .== 0], alpha=0.5, bins=bins, label="Noise molecules")
            show(io, MIME("text/html"), p_conf)
        end

        # Assignment confidence
        if :assignment_confidence in propertynames(df_res)
            p_conf = Plots.histogram(df_res.assignment_confidence[assignment .> 0], bins=50, legend=false)
            p_conf = Plots.vline!([0.95], xlabel="Assignment confidence", ylabel="#Molecules", xlims=(-0.01, 1.03),
                title="Assignment confidence per real molecules")
            show(io, MIME("text/html"), p_conf)
        end

        println(io, "

")

        # Num. of molecules per cell
        n_mols_per_cell = count_array(assignment, drop_zero=true)
        p_n_mols = Plots.histogram(n_mols_per_cell[(n_mols_per_cell .> 1) .& (n_mols_per_cell .< quantile(n_mols_per_cell, 0.99) / 0.99)],
            title="Num. molecules per cell", xlabel="Num. molecules per cell", ylabel="Num. cells", label=:none)
        show(io, MIME("text/html"), p_n_mols)
        println(io, "")
    end
end

function plot_transcript_assignment_panel(df_res::DataFrame, assignment::Vector{Int}, args::Dict;
        clusters::Vector{Int}, prior_polygons::Array{Matrix{Float64}, 1}, gene_colors::Union{Vector{<:Colors.ColorTypes.Color}, Nothing}=nothing)
    if gene_colors === nothing
        @info "Estimating local colors"
        gene_colors = gene_composition_colors(df_res, args["gene-composition-neigborhood"])
    end

    @info "Estimating boundary polygons" # For some parameters, this step takes a lot of time and memory
    grid_step = args["scale"] / args["min-pixels-per-cell"];
    polygons = boundary_polygons(df_res, assignment; grid_step=grid_step, bandwidth=args["scale"]/10);

    @info "Plot transcript assignment"
    gc_plot = plot_dataset_colors(df_res, gene_colors; polygons=polygons, prior_polygons=prior_polygons, min_molecules_per_cell=args["min-molecules-per-cell"],
        min_pixels_per_cell=args["min-pixels-per-cell"], title="Local expression similarity", alpha=0.5)

    clust_plot = nothing
    if !isempty(clusters)
        clust_plot = plot_dataset_colors(df_res, gene_colors; polygons=polygons, prior_polygons=prior_polygons, annotation=clusters, min_molecules_per_cell=args["min-molecules-per-cell"],
            min_pixels_per_cell=args["min-pixels-per-cell"], title="Molecule clustering", alpha=0.5)
    end

    open(append_suffix(args["output"], "borders.html"), "w") do io
        show(io, MIME("text/html"), gc_plot)

        if clust_plot !== nothing
            show(io, MIME("text/html"), clust_plot)
        end
    end

    return polygons
end

function parse_prior_assignment(df_spatial::DataFrame, args::Dict{String, Any})
    prior_col = Symbol(args["prior_segmentation"][2:end])
    prior_segmentation = df_spatial[!, prior_col]
    try
        prior_segmentation = Int.(prior_segmentation);
    catch
        error("The prior segmentation column '$prior_col' must be of integer type")
    end
    if minimum(prior_segmentation) < 0
        error("The prior segmentation column '$prior_col' must not contain negative numbers")
    end

    filter_segmentation_labels!(prior_segmentation, min_molecules_per_segment=args["min-molecules-per-segment"])
    if args["scale"] === nothing
        pos_data = position_data(df_spatial)
        scale, scale_std = estimate_scale_from_assignment(pos_data, prior_segmentation;
            min_mols_per_cell=args["min-molecules-per-cell"])
    else
        scale, scale_std = args["scale"], args["scale-std"]
    end

    return prior_segmentation, Matrix{Float64}[], scale, scale_std
end

function load_prior_segmentation(df_spatial::DataFrame, args::Dict{String, Any})
    if args["prior_segmentation"][1] == ':'
        return parse_prior_assignment(df_spatial, args)
    end

    @info "Loading segmentation mask..."

    prior_seg_labels = load_segmentation_mask(args["prior_segmentation"])
    prior_segmentation = Int.(staining_value_per_transcript(df_spatial, prior_seg_labels));
    filter_segmentation_labels!(prior_seg_labels, prior_segmentation; min_molecules_per_segment=args["min-molecules-per-segment"])
    GC.gc()

    scale, scale_std = args["scale"], args["scale-std"]
    if args["estimate-scale-from-centers"]
        scale, scale_std = estimate_scale_from_centers(prior_seg_labels)
    end
    @info "Done"

    @info "Estimating prior segmentation polygons..."
    prior_polygons = extract_polygons_from_label_grid(Matrix(prior_seg_labels[1:3:end, 1:3:end]); grid_step=3.0) # subset to save memory and time
    @info "Done"

    return prior_segmentation, prior_polygons, scale, scale_std
end

function get_baysor_run_str()::String
    pkg_info = Pkg.dependencies()[Base.UUID("cc9f9468-1fbe-11e9-0acf-e9460511877c")]
    pkg_str = "v$(pkg_info.version)"
    try
        repo = LibGit2.GitRepo(pkg_info.source)
        hash_str = "$(LibGit2.GitShortHash(LibGit2.peel(LibGit2.GitCommit, LibGit2.head(repo))))"
        pkg_str = "$pkg_str [$hash_str]"
    catch err
    end

    return "($(Dates.Date(Dates.now()))) Run Baysor $pkg_str"
end

function run_cli_main(args::Union{Nothing, Array{String, 1}}=nothing)
    args_str = join(args, " ")
    args = parse_configs(args)

    # Dump parameters

    if args["config"] !== nothing
        dump_dst = abspath(append_suffix(args["output"], "config.toml"))
        if abspath(args["config"]) != dump_dst
            cp(args["config"], dump_dst, force=true)
        end
    end

    open(append_suffix(args["output"], "params.dump"), "w") do f
        println(f, "# CLI params: `$args_str`")
        TOML.print(f, Dict(k => (v !== nothing) ? ((typeof(v) === Symbol) ? String(v) : v) : "" for (k,v) in args))
    end

    # Set up logger

    log_file = open(append_suffix(args["output"], "log.log"), "w")
    Base.CoreLogging.global_logger(DoubleLogger(log_file, stdout; force_flush=true))

    # Set up plotting

    ENV["GKSwstype"] = "100"; # Disable output device
    Plots.gr()

    # Run algorithm

    @info get_baysor_run_str()
    @info "Loading data..."
    df_spatial, gene_names = load_df(args, filter_cols=false)
    df_spatial[!, :molecule_id] = 1:size(df_spatial, 1)

    @info "Loaded $(size(df_spatial, 1)) transcripts"
    if size(df_spatial, 1) != size(unique(df_spatial), 1)
        @warn "$(size(df_spatial, 1) - size(unique(df_spatial), 1)) records are duplicates. You may need to filter them beforehand."
    end

    if args["gene-composition-neigborhood"] === nothing
        args["gene-composition-neigborhood"] = default_param_value(:composition_neighborhood, args["min-molecules-per-cell"], n_genes=length(gene_names))
    end

    prior_polygons = Matrix{Float64}[]
    if args["prior_segmentation"] !== nothing
        df_spatial[!, :prior_segmentation], prior_polygons, args["scale"], args["scale-std"] = load_prior_segmentation(df_spatial, args)
    end
    GC.gc()

    bm_data_arr = BmmData[]
    confidence_nn_id = default_param_value(:confidence_nn_id, args["min-molecules-per-cell"])

    @info "Estimating noise level"
    append_confidence!(df_spatial, (args["prior_segmentation"]===nothing ? nothing : df_spatial.prior_segmentation), nn_id=confidence_nn_id, prior_confidence=args["prior-segmentation-confidence"])
    @info "Done"

    max_diffs, change_fracs = nothing, nothing
    if args["n-clusters"] > 1
        @info "Clustering molecules..."
        # , adjacency_type=:both, k_adj=fmax(1, div(args["min-molecules-per-cell"], 2))
        adjacent_points, adjacent_weights = build_molecule_graph_normalized(df_spatial, :confidence, filter=false);

        mol_clusts = cluster_molecules_on_mrf(df_spatial, adjacent_points, adjacent_weights;
            n_clusters=args["n-clusters"], weights_pre_adjusted=false)

        df_spatial[!, :cluster] = mol_clusts.assignment;
        max_diffs, change_fracs = mol_clusts.diffs, mol_clusts.change_fracs

        # TODO: store mol_clusts.exprs at bm_data.misc?
        @info "Done"
    end

    bm_data_arr = initial_distribution_arr(df_spatial; n_frames=args["n-frames"], scale=args["scale"], scale_std=args["scale-std"],
            n_cells_init=args["num-cells-init"], prior_seg_confidence=args["prior-segmentation-confidence"],
            min_molecules_per_cell=args["min-molecules-per-cell"], confidence_nn_id=0);

    history_depth = round(Int, args["iters"] * 0.1)
    bm_data = run_bmm_parallel!(bm_data_arr, args["iters"], new_component_frac=args["new-component-fraction"], new_component_weight=args["new-component-weight"],
                                min_molecules_per_cell=args["min-molecules-per-cell"], assignment_history_depth=history_depth);

    @info "Processing complete."

    # Save results

    segmentated_df = get_segmentation_df(bm_data, gene_names)
    cell_stat_df = get_cell_stat_df(bm_data, segmentated_df; add_qc=true)

    @info "Estimating local colors"
    gene_colors = gene_composition_colors(bm_data.x, args["gene-composition-neigborhood"])
    segmentated_df[!, :ncv_color] = "#" .* Colors.hex.(gene_colors)

    @info "Saving results to $(args["output"])"
    CSV.write(args["output"], segmentated_df[sortperm(segmentated_df.molecule_id), :]);
    CSV.write(append_suffix(args["output"], "cell_stats.csv"), cell_stat_df);

    cm = convert_segmentation_to_counts(composition_data(bm_data), bm_data.assignment; gene_names=gene_names)
    count_str = join(names(cm), "\t") * "\n" * join([join(["$v" for v in r], '\t') for r in eachrow(cm)], '\n');
    open(append_suffix(args["output"], "counts.tsv"), "w") do f; print(f, count_str) end

    if args["plot"]
        plot_diagnostics_panel(segmentated_df, segmentated_df.cell, bm_data.tracer, args; max_diffs=max_diffs, change_fracs=change_fracs)
        polygons = plot_transcript_assignment_panel(bm_data.x, bm_data.assignment, args; clusters=bm_data.cluster_per_molecule, prior_polygons=prior_polygons,
            gene_colors=gene_colors)

        if args["save-polygons"] !== nothing
            if lowercase(args["save-polygons"]) == "geojson"
                save_polygons_to_geojson(polygons, append_suffix(args["output"], "polygons.json"))
            else
                @warn "Unknown polygon format: $(args["save-polygons"])"
            end
        end
    end

    @info "All done!"

    close(log_file)

    return 0
end

## Pre-view

function parse_preview_commandline(args::Union{Nothing, Array{String, 1}}=nothing)
    s = ArgParseSettings(prog="baysor preview")
    @add_arg_table! s begin
        "--config", "-c"
            help = "TOML file with config"
        "--x-column", "-x"
            help = "Name of x column. Overrides the config value."
        "--y-column", "-y"
            help = "Name of gene column. Overrides the config value."
        "--gene-column", "-g"
            help = "Name of gene column. Overrides the config value."
        "--min-molecules-per-cell", "-m"
            help = "Minimal number of molecules for a cell to be considered as real. It's an important parameter, as it's used to infer several other parameters. Overrides the config value."
            arg_type = Int
        "--min-pixels-per-cell"
            help = "Minimal number of pixels per cell. Used to estimate size of the dataset plot."
            arg_type = Int
            default = 15
        "--output", "-o"
            help = "Name of the output file or path to the output directory"
            default = "preview.html"
        "coordinates"
            help = "CSV file with coordinates of transcripts and gene type"
            required = true
        # TODO: add gene-composition-neighborhood
        # TODO: add verbosity level
    end

    return (args === nothing) ? parse_args(s) : parse_args(args, s)
end

function parse_preview_configs(args::Union{Nothing, Array{String, 1}}=nothing)
    r = parse_preview_commandline(args)
    if r["config"] !== nothing
        extend_params_with_config!(r, parse_toml_config(r["config"]))
    else
        extend_params_with_config!(r, get_default_config())
    end

    for k in ["gene-column", "x-column", "y-column"]
        r[k] = Symbol(r[k])
    end

    if isdir(r["output"]) || isdirpath(r["output"])
        r["output"] = joinpath(r["output"], "preview.html")
    end

    return r
end

function run_cli_preview(args::Union{Nothing, Array{String, 1}}=nothing)
    args = parse_preview_configs(args)

    # Set up logger

    log_file = open(append_suffix(args["output"], "preview_log.log"), "w")
    Base.CoreLogging.global_logger(DoubleLogger(log_file, stdout; force_flush=true))

    # Set up plotting

    ENV["GKSwstype"] = "100"; # Disable output device
    Plots.gr()

    # Run preview

    @info get_baysor_run_str()
    @info "Loading data..."
    args["min-molecules-per-gene"] = 0
    df_spatial, gene_names = load_df(args)

    @info "Loaded $(size(df_spatial, 1)) transcripts"
    if size(df_spatial, 1) != size(unique(df_spatial), 1)
        @warn "$(size(df_spatial, 1) - size(unique(df_spatial), 1)) records are duplicates. You may need to filter them beforehand."
    end

    @info "Estimating noise level"
    confidence_nn_id = default_param_value(:confidence_nn_id, args["min-molecules-per-cell"])
    edge_lengths, confidences, d1, d2 = append_confidence!(df_spatial, nn_id=confidence_nn_id) # TODO: use segmentation mask if available here
    @info "Done"

    @info "Estimating local neighborhoods"

    if args["gene-composition-neigborhood"] === nothing
        args["gene-composition-neigborhood"] = default_param_value(:composition_neighborhood, args["min-molecules-per-cell"], n_genes=length(gene_names))
    end

    @info "Estimating local colors"
    gene_colors = gene_composition_colors(df_spatial, args["gene-composition-neigborhood"])

    @info "Building transcript plots"
    gc_plot = plot_dataset_colors(df_spatial, gene_colors; min_molecules_per_cell=args["min-molecules-per-cell"],
        min_pixels_per_cell=args["min-pixels-per-cell"], title="Local expression similarity")

    conf_colors = map_to_colors(confidences, lims=(0.0, 1.0), palette=Colors.diverging_palette(10, 250, s=0.75, w=1.0));
    cc_plot = plot_dataset_colors(df_spatial, conf_colors[:colors]; min_molecules_per_cell=args["min-molecules-per-cell"],
        min_pixels_per_cell=args["min-pixels-per-cell"], title="Transcript confidence")

    @info "Building gene structure plot"
    gene_structure_plot = plot_gene_structure(df_spatial, gene_names)
    ## Plots

    n_tr_plot = plot_num_transcript_overview(df_spatial.gene, confidences, gene_names)
    noise_dist_plot = plot_noise_estimation_diagnostics(edge_lengths, confidences, d1, d2, confidence_nn_id=confidence_nn_id)

    @info "Plotting"

    open(args["output"], "w") do io
        print(io, """
            <!DOCTYPE html>

Content

Transcript plots

Noise level

Gene structure

            """)

        println(io, "
Transcript plots

")
        show(io, MIME("text/html"), gc_plot)
        println(io, "
")

        println(io, "
\n
Noise level

")
        show(io, MIME("text/html"), cc_plot)
        println(io, "
")
        show(io, MIME("text/html"), noise_dist_plot)
        println(io, "
")
        println(io, "Minimal noise level=$(round(100 * mean(confidences .< 0.01), sigdigits=3))%. ",
            "Expected noise level=$(round(100 * mean(1 .- confidences), sigdigits=2))%.")
        println(io, "
")

        println(io, "
\n
Gene structure

")
        show(io, MIME("text/html"), n_tr_plot)
        println(io, "
")
        show(io, MIME("text/html"), gene_structure_plot)
        println(io, "")
        println(io, "")
    end

    @info "All done!"

    close(log_file)

    return 0
end

## All

run_cli(args::String) = run_cli(String.(Base.split(args)))

function run_cli(args::Vector{String}=ARGS)::Cint
    help_message = "Usage: baysor <command> [options]\n\nCommands:\n\trun\t\trun segmentation of the dataset\n\tpreview\t\tgenerate preview diagnostics of the dataset\n"

    debug = false
    if "--debug" in args
        args = args[args .!= "--debug"]
        debug = true
    end

    try
        if (length(args) == 0) || (length(args) == 1) && (args[1] == "-h" || args[1] == "--help")
            println(help_message)
            return 0
        end

        if args[1] == "run"
            return run_cli_main(args[2:end])
        end

        if args[1] == "preview"
            return run_cli_preview(args[2:end])
        end

        @error "Can't parse argument $(args[1])"
        println(help_message)
        return 1
    catch err
        if debug
            rethrow()
        else
            @error("$err\n\n" * join(["$s" for s in stacktrace(catch_backtrace())], "\n"))
        end
    end

    return 2
end

julia_main()::Cint = run_cli()


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c8736c18-1f88-11eb-1025-91cf3d825edf
using Pkg

# ╔═╡ 490206e6-1f43-11eb-39db-05eb9f02b2b2
begin
	using TextAnalysis
	using CSV
	using DataFrames
	using Pipe: @pipe
	using Plots
end

# ╔═╡ 4d1d0a34-2478-11eb-3c8e-b1e62b58dcff
using TextAnalysis: NaiveBayesClassifier, fit!, predict

# ╔═╡ c251d572-247e-11eb-29d9-3d382c09f425
begin
	using HTTP
	using JSON3
end

# ╔═╡ 9ecf464a-1ff4-11eb-2c0a-dfa5efc8e305
md"""
# Sentiment analysis using TextAnalysis.jl

This notebook explores the [Twitter US Airline Sentiment](https://www.kaggle.com/crowdflower/twitter-airline-sentiment) data set from Kaggle.
We will use [TextAnalysis.jl](https://github.com/JuliaText/TextAnalysis.jl) as the primary tool for analyzing textual data.
"""

# ╔═╡ f500ca32-1f88-11eb-0f34-2d84f8316de7
dir = "/Users/tomkwong/Julia/HumansOfJulia-WeeklyContest/Week2-TextAnalysis.jl/tk3369"

# ╔═╡ cbed3fd6-1f88-11eb-1267-d172e19f441f
Pkg.activate(dir)

# ╔═╡ ff756d94-1ff4-11eb-25ae-5983ea00b7e3
md"""
## Loading data
"""

# ╔═╡ ffc6a542-1f88-11eb-284b-056a0af0a139
cd(dir)

# ╔═╡ f13c9a7a-1f88-11eb-3c7c-5371c8fe8203
df = DataFrame(CSV.File("data/Tweets.csv"))

# ╔═╡ e29c55f0-1ff1-11eb-2f0c-97d6c9fb83c3
md"""
## Data Wrangling

We will take a look at the data and get a little more understanding about what's going on in this data set.
"""

# ╔═╡ ea7dadbc-1ff1-11eb-2397-1fe07e2899be
describe(df, :eltype, :nmissing, :first => first)

# ╔═╡ 2fa37940-1ff3-11eb-01ef-813c561d5679
let x = combine(groupby(df, :airline_sentiment), nrow)
	bar(x.airline_sentiment, x.nrow;
		title = "Airline Sentiment",
		label = :none,
		legend = :topright)
end

# ╔═╡ 4797ae3c-1ff2-11eb-3f31-f32f365ea0ee
histogram(df.airline_sentiment_confidence;
	legend = nothing,
	title = "Airline Sentiment Confidence")

# ╔═╡ 747a707c-1ff4-11eb-2681-6b6787174d0e
let x = combine(groupby(dropmissing(df, :negativereason), :negativereason), nrow)
	bar(x.negativereason, x.nrow;
		title = "Negative Reasons", label = :none, xrotation = 45)
end

# ╔═╡ 5ca1e714-1f8b-11eb-30fd-dbdd6e574805
md"""
# Examining tweets

The CSV file contains over 14,000 tweets. Let's quickly examine some individual data.
"""

# ╔═╡ 097e9b2c-1ff4-11eb-255d-1d0816f3ddee
md"""
Before we go further, it would be nice to display a single record in table format. We can define a `table` function that converts an indexable object into Markdown format, which can be displayed in this Pluto notebook.
"""

# ╔═╡ 74192f10-1fe0-11eb-3eea-1d9cbc0c4c3b
function table(nt)
	io = IOBuffer()
	println(io, "|name|value|")
	println(io, "|---:|:----|")
	for k in keys(nt)
		println(io, "|`", k, "`|", nt[k], "|")
	end
	return Markdown.parse(String(take!(io)))
end;

# ╔═╡ 32974c34-1ff4-11eb-18ea-67999c24906a
md"""
Here, we will define a variable called `row` and bind it to a slider for quick experimentation.
"""

# ╔═╡ e97bc70a-2002-11eb-3f57-ab1838225b0e
@bind row html"""<input type="range" min="1" max="100" value="36"/>"""

# ╔═╡ db9272f6-1ff3-11eb-12e7-0fa2c8837975
"Current Record: $row"

# ╔═╡ cbcac7e6-1fe0-11eb-0820-9f4c6f3fa51f
table(df[row, :])

# ╔═╡ dd832f9c-1f8a-11eb-001b-3d8be911d80a
md"""
As an example, record #36 has the tweet text as:
```
Nice RT @VirginAmerica: Vibe with the moodlight from takeoff to touchdown. #MoodlitMonday #ScienceBehindTheExperience http://t.co/Y7O0uNxTQP
```

This is a tricky one because it contains all of the followings:
- mention (`@VirginAmeria`)
- hash tag (`#MoodlitMonday` and `#ScienceBehindTheExperience`)
- URL (`http://t.co/Y7O0uNxTQP`)

Technically `RT` is a shorthand for "retweet" so perhaps it should be expanded but let's not worry about that for now.
"""

# ╔═╡ b231f13c-1ff0-11eb-278d-43db738ca57f
md"""
## Handling mentions and hashtags
"""

# ╔═╡ b636fc32-1ff5-11eb-048e-cbc4a77c5f75
md"""
If we just ignore these problems then it can be a disaster.
"""

# ╔═╡ 6cb54136-1f8c-11eb-21d0-1542ddd12437
let s = df[36, :text]
	sd = StringDocument(lowercase(s))
	op = 0x00
	op |= strip_punctuation
	op |= strip_stopwords
	op |= strip_html_tags
	prepare!(sd, op)
	stem!(sd)
	table(ngrams(sd))
end

# ╔═╡ cd6619b0-1ff5-11eb-19f5-092f2f06b172
md"""
Right off the bat, I can see some problems here. It seems that when I stripped punctuations, it also took the `@` and `#` signs away. The URL also became weird.  Oh yeah, that's what stripping punctuation means, right? :-)
"""

# ╔═╡ 8a3e8bbc-1ff6-11eb-2a12-c3fba022267b
md"""
##### Extracting mentions, hash tags, and URL's.
This neat idea came from José Bayoán Santiago Calderón when I asked the question on Slack. Let's define some functions using regular expressions.
"""

# ╔═╡ c72afbc0-2474-11eb-3f52-e76f60ecd80d
const regexp = Dict(
	:mention => r"@\w+",
	:hashtag => r"#\w+",
	:url => r"http[s]?://(?:[a-zA-Z]|[0-9]|[$-_@.&+]|[!*\(\),]|(?:%[0-9a-fA-F][0-9a-fA-F]))+"
);

# ╔═╡ fd1a4328-1ff7-11eb-0b7f-2d62c0e9d0df
function extract_tokens(s, token_type)
	return collect(x.match for x in eachmatch(regexp[token_type], s))
end

# ╔═╡ b3b62480-1ff8-11eb-3ef9-7985882f19d4
function remove_tokens(s)
	for re in values(regexp)
		s = replace(s, re => "")
	end
	return s
end

# ╔═╡ d4ea52ea-24b0-11eb-11b1-f3079e437485
md"##### Create a new data frame with extracted and clean text fields"

# ╔═╡ 457104ae-1ff8-11eb-04d8-536205dc3fc9
begin
	df2 = DataFrame()
	df2.airline_sentiment = df.airline_sentiment
	df2.text = df.text
	df2.mentions = extract_tokens.(df.text, :mention)
	df2.hashtags = extract_tokens.(df.text, :hashtag)
	df2.urls = extract_tokens.(df.text, :url)
	df2.clean_text = lowercase.(remove_tokens.(df.text))
	df2
end;

# ╔═╡ 91e59d40-1ff8-11eb-0fab-df3732a771c2
table(df2[36, :])

# ╔═╡ de6d7008-215a-11eb-1082-dd3adfb676e6
md"""
As you can see, the mentions/hashtags/urls are extracted into separate columns in the data frame. The `clean_text` field contains the cleaned version of `text`.
"""

# ╔═╡ 5c30f168-2124-11eb-3343-79f2765bd382
md"""
## Using Naive Bayes Classifier
"""

# ╔═╡ 6ebd952a-2124-11eb-30ee-2b80714ccb11
md"""
In our data frame, we already have a column `x_string_doc` with `StringDocuments` values. So we can just fit them to the classifier.
"""

# ╔═╡ 238efc32-1ff9-11eb-3887-17016268dfe2
function create_string_doc(s)
	sd = StringDocument(s)
	op = 0x00
	op |= strip_punctuation
	op |= strip_stopwords
	op |= strip_html_tags
	prepare!(sd, op)
	stem!(sd)
	return sd
end

# ╔═╡ 76ab55ba-2124-11eb-3999-71ce72708dd3
model = let
	classes = unique(df2.airline_sentiment)
	nbc = NaiveBayesClassifier(classes)
	for (clean_text, class) in zip(df2.clean_text , df2.airline_sentiment)
		sd = create_string_doc(clean_text)
		fit!(nbc, sd, class)
	end
	nbc
end;

# ╔═╡ 983678c6-2125-11eb-07a7-29321450d287
md"Let's create a model test function and then try our predictor for a few simple test cases."

# ╔═╡ cdcc1764-212c-11eb-1a24-cb5f0aa323a5
function test_model(model, tweets)
	df = DataFrame(text = tweets)
	df.doc = TextAnalysis.text.(create_string_doc.(remove_tokens.(tweets)))
	df.analysis = predict.(Ref(model), df.doc)

	df.positive = getindex.(df.analysis, "positive")
	df.negative = getindex.(df.analysis, "negative")
	df.neutral = getindex.(df.analysis, "neutral")

	select!(df, Not(:analysis))

	return df
end;

# ╔═╡ 058b9f92-2486-11eb-3853-d7f83cf90198
let
	tweets = [
		"whatever airline sucks!",
		"i love @american service :-)",
		"just ok",
		"hello world"]
	test_model(model, tweets)
end

# ╔═╡ 02494270-2130-11eb-1342-390a8240d584
md"""
## Determining accuracy

How well does the Naive Bayes Classifier work?
"""

# ╔═╡ b4aa64ac-2136-11eb-142f-5dfb445e2e96
md"""
As the `predict` function returns a `Dict` object with the probabilities assigned to each class, we need to choose the best option. Let's define a function for that.
"""

# ╔═╡ 54a76b44-2130-11eb-2680-211869c363da
function predict_and_choose(c::NaiveBayesClassifier, sd::StringDocument)
	val = predict(c, sd)
	return argmax(val)
end;

# ╔═╡ 0d66ccd4-2137-11eb-1f8e-2d03c16e29d9
md"Now, make prediction over all 14K tweets."

# ╔═╡ 0fbd02e8-2130-11eb-1f17-d77c0e80fa2a
yhat = let sds = create_string_doc.(lowercase.(remove_tokens.(df2.text)))
	predict_and_choose.(Ref(model), sds)
end;

# ╔═╡ f27be3ec-2130-11eb-0114-33bf20368dba
hits = count(df2.airline_sentiment .== yhat)

# ╔═╡ 55f555f6-2137-11eb-2312-99fbac7dad6e
misses = length(yhat) - hits

# ╔═╡ 669e0fd8-2137-11eb-1eb4-01caefb25d55
wayoff = count(
			(df2.airline_sentiment .!== yhat) .&
			(df2.airline_sentiment .!== "neutral") .&
			(yhat .!== "neutral"))

# ╔═╡ c283de5e-2132-11eb-2cb6-df48c17efd1d
accuracy_percentage = hits / (hits + misses) * 100

# ╔═╡ ed660bee-2132-11eb-3217-a58156b889a9
slightly_off_percentage = (misses - wayoff) / (hits + misses) * 100

# ╔═╡ ffaba636-2132-11eb-1593-63aab54d6c98
way_off_percentage = wayoff / (hits + misses) * 100

# ╔═╡ 398a14bc-2135-11eb-0351-a314bfbb0db4
bar(["hits", "misses slightly", "misses way-off"],
	[accuracy_percentage, slightly_off_percentage , way_off_percentage];
	legend = :none)

# ╔═╡ b9da38a8-247e-11eb-1fcd-0b19be7ece16
md"# Analyze some random tweets"

# ╔═╡ c57bb7b8-247e-11eb-3055-c7c83a9bfe93
token = readline("/Users/tomkwong/.twitter-bearer");

# ╔═╡ 0863fb44-247f-11eb-3fc9-935c87fdff28
response = HTTP.get("https://api.twitter.com/1.1/search/tweets.json?q=lang%3Aen%20flight", ["authorization" => "Bearer $token"]);

# ╔═╡ 19c808ee-247f-11eb-2082-59aa94a3e77c
data = JSON3.read(response.body);

# ╔═╡ 5a873d72-2485-11eb-10e1-bb119c1e2795
data[:statuses][1]

# ╔═╡ 32fb7e76-2485-11eb-3f0b-6f779c314a81
tweets = [x.text for x in data.statuses]

# ╔═╡ de5600b6-2485-11eb-1111-a729caa3a913
result = test_model(model, tweets)

# ╔═╡ 82d9d400-24a4-11eb-0c28-950af58adee4
table(result[9,:])

# ╔═╡ Cell order:
# ╟─9ecf464a-1ff4-11eb-2c0a-dfa5efc8e305
# ╠═c8736c18-1f88-11eb-1025-91cf3d825edf
# ╠═f500ca32-1f88-11eb-0f34-2d84f8316de7
# ╠═cbed3fd6-1f88-11eb-1267-d172e19f441f
# ╠═490206e6-1f43-11eb-39db-05eb9f02b2b2
# ╟─ff756d94-1ff4-11eb-25ae-5983ea00b7e3
# ╠═ffc6a542-1f88-11eb-284b-056a0af0a139
# ╠═f13c9a7a-1f88-11eb-3c7c-5371c8fe8203
# ╟─e29c55f0-1ff1-11eb-2f0c-97d6c9fb83c3
# ╠═ea7dadbc-1ff1-11eb-2397-1fe07e2899be
# ╠═2fa37940-1ff3-11eb-01ef-813c561d5679
# ╠═4797ae3c-1ff2-11eb-3f31-f32f365ea0ee
# ╠═747a707c-1ff4-11eb-2681-6b6787174d0e
# ╟─5ca1e714-1f8b-11eb-30fd-dbdd6e574805
# ╟─097e9b2c-1ff4-11eb-255d-1d0816f3ddee
# ╠═74192f10-1fe0-11eb-3eea-1d9cbc0c4c3b
# ╟─32974c34-1ff4-11eb-18ea-67999c24906a
# ╠═e97bc70a-2002-11eb-3f57-ab1838225b0e
# ╠═db9272f6-1ff3-11eb-12e7-0fa2c8837975
# ╠═cbcac7e6-1fe0-11eb-0820-9f4c6f3fa51f
# ╟─dd832f9c-1f8a-11eb-001b-3d8be911d80a
# ╟─b231f13c-1ff0-11eb-278d-43db738ca57f
# ╟─b636fc32-1ff5-11eb-048e-cbc4a77c5f75
# ╠═6cb54136-1f8c-11eb-21d0-1542ddd12437
# ╟─cd6619b0-1ff5-11eb-19f5-092f2f06b172
# ╟─8a3e8bbc-1ff6-11eb-2a12-c3fba022267b
# ╠═c72afbc0-2474-11eb-3f52-e76f60ecd80d
# ╠═fd1a4328-1ff7-11eb-0b7f-2d62c0e9d0df
# ╠═b3b62480-1ff8-11eb-3ef9-7985882f19d4
# ╟─d4ea52ea-24b0-11eb-11b1-f3079e437485
# ╠═457104ae-1ff8-11eb-04d8-536205dc3fc9
# ╠═91e59d40-1ff8-11eb-0fab-df3732a771c2
# ╟─de6d7008-215a-11eb-1082-dd3adfb676e6
# ╟─5c30f168-2124-11eb-3343-79f2765bd382
# ╟─6ebd952a-2124-11eb-30ee-2b80714ccb11
# ╠═4d1d0a34-2478-11eb-3c8e-b1e62b58dcff
# ╠═238efc32-1ff9-11eb-3887-17016268dfe2
# ╠═76ab55ba-2124-11eb-3999-71ce72708dd3
# ╟─983678c6-2125-11eb-07a7-29321450d287
# ╠═cdcc1764-212c-11eb-1a24-cb5f0aa323a5
# ╠═058b9f92-2486-11eb-3853-d7f83cf90198
# ╟─02494270-2130-11eb-1342-390a8240d584
# ╟─b4aa64ac-2136-11eb-142f-5dfb445e2e96
# ╠═54a76b44-2130-11eb-2680-211869c363da
# ╟─0d66ccd4-2137-11eb-1f8e-2d03c16e29d9
# ╠═0fbd02e8-2130-11eb-1f17-d77c0e80fa2a
# ╠═f27be3ec-2130-11eb-0114-33bf20368dba
# ╠═55f555f6-2137-11eb-2312-99fbac7dad6e
# ╠═669e0fd8-2137-11eb-1eb4-01caefb25d55
# ╠═c283de5e-2132-11eb-2cb6-df48c17efd1d
# ╠═ed660bee-2132-11eb-3217-a58156b889a9
# ╠═ffaba636-2132-11eb-1593-63aab54d6c98
# ╠═398a14bc-2135-11eb-0351-a314bfbb0db4
# ╟─b9da38a8-247e-11eb-1fcd-0b19be7ece16
# ╠═c251d572-247e-11eb-29d9-3d382c09f425
# ╠═c57bb7b8-247e-11eb-3055-c7c83a9bfe93
# ╠═0863fb44-247f-11eb-3fc9-935c87fdff28
# ╠═19c808ee-247f-11eb-2082-59aa94a3e77c
# ╠═5a873d72-2485-11eb-10e1-bb119c1e2795
# ╠═32fb7e76-2485-11eb-3f0b-6f779c314a81
# ╠═de5600b6-2485-11eb-1111-a729caa3a913
# ╠═82d9d400-24a4-11eb-0c28-950af58adee4

ext: jl

max_stars_repo_name: Humans-of-Julia/Challenges

max_issues_repo_name: Humans-of-Julia/WeeklyContest

max_forks_repo_name: Humans-of-Julia/WeeklyContest


text:
# auxiliary functions
find_negatives(targets) = findall(vec(targets) .== 0)
find_positives(targets) = findall(vec(targets) .== 1)

@nograd find_negatives, find_positives

# surrogate functions
hinge(x, ϑ::Real = 1) = max(zero(x), 1 + ϑ * x)
quadratic(x, ϑ::Real = 1) = max(zero(x), 1 + ϑ * x)^2

# objectives
function fnr(targets, scores, t, surrogate = quadratic)
    return mean(surrogate.(t .- scores[find_positives(targets)]))
end

function fpr(targets, scores, t, surrogate = quadratic)
    return mean(surrogate.(scores[find_negatives(targets)] .- t))
end

# root finding
function find_root(f::Function, x_0::Real)
    try
        x = Roots.secant_method(f, x_0)
        isnan(x) && error("secant_method failed")
        return x
    catch
        try
            x = Roots.find_zero(f, (x_0, Inf))
            isnan(x) && error("find_zero failed")
            return x
        catch
            xrange = range(x_0; stop = 10*x_0, length = 100)
            xzero, ind  = findmin(abs.(f.(xrange)))
            return xrange[ind]
        end
    end
end


text:
using BosonSampling
using BenchmarkTools
using BenchmarkPlots, StatsPlots
using Permanents
using LaTeXStrings

BenchmarkTools.DEFAULT_PARAMETERS.samples = 100
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1
BenchmarkTools.DEFAULT_PARAMETERS.evals = 1

suite_permanent = BenchmarkGroup(["string"])
for n in 15:25
    A = RandHaar(n)
    suite_permanent["n=$n"] = @benchmarkable ryser($A.U)
end

res = run(suite_permanent, verbose=true, seconds=1)
data_mean = []

for i in 15:25
    mean_res = mean(res["n=$i"])
    push!(data_mean, (mean_res.time)/10^9)
end

f = open("data_python.txt")
f = readlines(f)
python_data = split(f[1], " ")
python_data = map(e->parse(Float64,e), python_data)

f = open("data_matlab.txt")
f = readlines(f)
matlab_data = split(f[1], " ")
matlab_data = map(e->parse(Float64,e), matlab_data)

perm_fig = plot(xlabel="n", ylabel="time (s)", legend=:topleft, dpi=300)
scatter!(15:25, data_mean, label="Julia", dpi=300)
scatter!(15:25, python_data, label="python", dpi=300)
scatter!(15:25, matlab_data, label="Matlab", dpi=300)
savefig(perm_fig, "docs/src/benchmarks/compute_perm.png")

suite_sampler = BenchmarkGroup()
suite_sampler["cliffords_sampler"] = BenchmarkGroup(["string"])
suite_sampler["noisy_sampler"] = BenchmarkGroup(["string"])

for n in 15:25
    interf = RandHaar(n)

    input = Input{Bosonic}(first_modes(n,n))
    suite_sampler["cliffords_sampler"]["n=$n"] = @benchmarkable cliffords_sampler(input=$input, interf=$interf)

    input = Input{OneParameterInterpolation}(first_modes(n,n), 0.5)
    suite_sampler["noisy_sampler"]["n=$n,reflec=0.25"] = @benchmarkable noisy_sampler(input=$input, reflectivity=0.25, interf=$interf)
    suite_sampler["noisy_sampler"]["n=$n,reflec=0.75"] = @benchmarkable noisy_sampler(input=$input, reflectivity=0.75, interf=$interf)
end

res = run(suite_sampler, verbose=true, seconds = 1)
data_bosonic = []
data_noisy_25 = []
data_noisy_75 = []

for j in 15:25
    mean_bosonic = mean(res["cliffords_sampler"]["n=$j"])
    push!(data_bosonic, (mean_bosonic.time)/10^9)
    mean_noisy = mean(res["noisy_sampler"]["n=$j,reflec=0.25"])
    push!(data_noisy_25, (mean_noisy.time)/10^9)
    mean_noisy = mean(res["noisy_sampler"]["n=$j,reflec=0.75"])
    push!(data_noisy_75, (mean_noisy.time)/10^9)
end

fig_samp = plot(xlabel="n", ylabel="time (s)", legend=:topleft, dpi=300)
scatter!(15:25, data_bosonic, label="cliffords sampler", dpi=300)
scatter!(15:25, data_noisy_75, label=L"noisy sampler, $η=0.75", dpi=300)
savefig(fig_samp, "docs/src/benchmarks/sampler.png")


text:
using Makie

function main()
    ts = Node(0.001)
    scene = Scene()
    v = range(0, stop=4pi, length=50)
    f(v, t) = sin(v + t) # some function
    s = lines!(
        scene,
        lift(t -> f.(v, t), ts),
    )[end];

    for i = 1:length(v)
        push!(ts, i)
        sleep(1/24)
    end

    r = 1:100
    for i = 1:length(r)
        scene[:markersize] = r[i]
        force_update!()
        sleep(1/24)
    end
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


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

using Lint
using Base.Test

results = lintpkg("FEMMaterials")
if !isempty(results)
    info("Lint.jl is a tool that uses static analysis to assist in the development process by detecting common bugs and potential issues.")
    info("For this package, Lint.jl report is following:")
    display(results)
    info("For more information, see https://lintjl.readthedocs.io/en/stable/")
    warn("Package syntax test has failed.")
    @test isempty(results)
else
    info("Lint.jl: syntax check pass.")
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Dex"
version = v"2.27.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/dexidp/dex.git", "0f9e2888ab65c5b18c4881eaeeb7e38d997e1d92"),
    DirectorySource("bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mkdir -p "${bindir}" "${prefix}/share"
cd dex/
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
install_license LICENSE
make
mkdir -p $bindir
mv bin/dex "$bindir/dex${exeext}"
tar -czvf $prefix/share/webtemplates.tar.gz -C ./web static templates themes
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="musl"),
    Platform("x86_64", "linux"; libc="glibc")
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("dex", :dex),
    FileProduct("share/webtemplates.tar.gz", :webtemplates),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers = [:go, :c])

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
coefs = Float64.([1 1 1])

thumbstack_stencil = @def_stencil_expression(
            @sum(i,-2,2,(
                @sum(j,-2,2, c[max(abs(i),abs(j))]*D[x+i,y+j,z])
            )) + @sum(i,1,2, c[i]*(D[x,y,z+i]+D[x,y,z-i]))
)

st_inst = NewStencilInstance(thumbstack_stencil, 2, coefs)

radius = 2
nx = 5
ny = 5
nz = 5
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 5
## Actual time steps
data = CreateData(dx,dy,dz,radius)
dout = ApplyStencil(st_inst, data, t_steps)
#heatmap(dout1[:,:,Int(round(dz/2))])
o = Array(dout)

function cpu_thumbstack_stencil(data, radius,coefs)
    out = zeros(Float32,size(data))
    for i = (radius+1):(size(data,1)-radius)
        for j = (radius+1):(size(data,2)-radius)
            for k = (radius+1):(size(data,3)-radius)
                c = 0#coefs[1]*data[i,j,k]
                for rx = -radius:radius, ry = -radius:radius
                    cr = max(abs.([rx,ry])...)
                    c += coefs[cr+1]*(data[i+rx,j+ry,k])
                end
                for rz = 1:radius
                    c += coefs[rz+1]*(data[i,j,k+rz] + data[i,j,k-rz])
                end
                out[i,j,k] = c
            end
        end
    end
    return out
end

cout = cpu_thumbstack_stencil(data, radius, coefs)
for i = 2:t_steps
    global cout
    cout = cpu_thumbstack_stencil(cout, radius, coefs)
end

@assert(sum(abs.(cout .- o))./(reduce(*, size(data))) < 0.1)
mean_err = sum(abs.(cout .- o))./(reduce(*, size(data)))
# using Plots
# shift = 6
# heatmap(cout[:,:,Int(round(dz/2))+shift])
# heatmap(o[:,:,Int(round(dz/2))+shift])
#
# heatmap(o[:,:,Int(round(dz/2))+shift] .- cout[:,:,Int(round(dz/2))+shift])
#
#
# diff =  abs.(cout .- o)#./(reduce(*, size(data)))
# findmax(diff)
#


text:
# This functionality is derived directly from the MATLAB interface of SuiteSparse:GraphBLAS
# by Tim Davis.

function argminmax(A::GBMatOrTranspose{T}, minmax, dim) where {T}
    nrows, ncols = size(A)
    n = (dim == 2) ? ncols : nrows
    y = GBVector{T}(n)
    y[:] = 1
    if minmax == min || min == BinaryOps.MIN
        rig = Semirings.MIN_FIRST
    else
        rig = Semirings.MAX_FIRST
    end
    dim == 2 ? desc = nothing : desc = T0
    x = mul(A, y, rig; desc=desc) # x = [min|max](A)
    D = Diagonal(x)
    if dim == 1
        G = mul(A, D, Semirings.ANY_EQ)
    else
        G = mul(D, A, Semirings.ANY_EQ)
    end
    select!(NONZERO, G, G)
    mul(G, y, Semirings.MIN_SECONDI1; desc)
end


text:
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::StridedMatrix, B::StridedMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::StridedMatrix, A::AbstractStrideMatrix, B::StridedMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::StridedMatrix, A::StridedMatrix, B::AbstractStrideMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::AbstractStrideMatrix, B::StridedMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::StridedMatrix, B::AbstractStrideMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::StridedMatrix, A::AbstractStrideMatrix, B::AbstractStrideMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)
@inline LinearAlgebra.mul!(C::AbstractStrideMatrix, A::AbstractStrideMatrix, B::AbstractStrideMatrix) = (@gc_preserve(matmul!(C, A, B)); return C)

@inline function stridematmul(
    A::AbstractStrideMatrix{<:Any,<:Any,TA},
    B::AbstractStrideMatrix{<:Any,<:Any,TB}
) where {TA<:Base.HWReal,TB<:Base.HWReal}
    M, KA = size(A)
    KB, N = size(B)
    @assert KA == KB "Size mismatch."
    K = Octavian._select(KA, KB)
    TC = promote_type(TA,TB)
    C = StrideArray{TC}(undef, (M, N))
    @gc_preserve matmul!(C, A, B, One(), Zero(), nothing, (M, K, N), ArrayInterface.contiguous_axis(C))
    return C
end
@inline Base.:*(A::AbstractStrideMatrix{<:Any,<:Any,TA}, B::StridedMatrix{TB}) where {TA<:Base.HWReal,TB<:Base.HWReal} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::AbstractStrideMatrix{<:Any,<:Any,TA}, B::StridedMatrix{TB}) where {TA<:Union{Float32,Float64},TB<:Union{Float32,Float64}} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::StridedMatrix{TA}, B::AbstractStrideMatrix{<:Any,<:Any,TB}) where {TA<:Base.HWReal,TB<:Base.HWReal} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::StridedMatrix{TA}, B::AbstractStrideMatrix{<:Any,<:Any,TB}) where {TA<:Union{Float32,Float64},TB<:Union{Float32,Float64}} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::AbstractStrideMatrix{<:Any,<:Any,TA}, B::AbstractStrideMatrix{<:Any,<:Any,TB}) where {TA<:Base.HWReal,TB<:Base.HWReal} = @gc_preserve(stridematmul(A, B))
@inline Base.:*(A::AbstractStrideMatrix{<:Any,<:Any,TA}, B::AbstractStrideMatrix{<:Any,<:Any,TB}) where {TA<:Union{Float32,Float64},TB<:Union{Float32,Float64}} = @gc_preserve(stridematmul(A, B))

@inline extract_λ(a) = a
@inline extract_λ(a::UniformScaling) = a.λ
@inline function Base.:*(A::AbstractStrideArray{S,D,T}, bλ::Union{Tb,UniformScaling{Tb}}) where {S,D,T<:VectorizationBase.NativeTypes,Tb <: Real}
    mv = similar(A)
    b = T(extract_λ(bλ))
    @avx for i ∈ eachindex(A)
        mv[i] = A[i] * b
    end
    mv
end

function LinearAlgebra.mul!(
    C::AbstractStrideMatrix{<:Any,<:Any,T},
    A::LinearAlgebra.Diagonal{T,<:AbstractStrideVector{<:Any,T}},
    B::StridedMatrix{T}
) where {T}
    M, K, N = matmul_axes(C, A, B)
    MandK = ArrayInterface._pick_range(M, K)
    vA = parent(A)
    @avx for n ∈ N, m ∈ 1:MandK
        C[m,n] = vA[m] * B[m,n]
    end
    C
end
@inline function Base.:*(
    A::LinearAlgebra.Diagonal{T,<:AbstractVector{T}},
    B::AbstractStrideMatrix{<:Any,<:Any,T}
) where {T}
    mul!(similar(B), A, B)
end


text:
module Harmonie

import JSONSchema, JSON, Glob
import Base.isvalid

export isvalid, diagnose

const moduledir=@__DIR__
const schemafile="$moduledir/../docs/harmonie.schema.json"
const testbeddir=joinpath(moduledir,"../test/testbed_configurations/")
const configdir=joinpath(moduledir, "../test/harmonie_configurations")
const configurations = getindex.(splitext.(readdir(configdir)),1)

"""
    __init__()

Initialize Harmonie.schema from Harmonie.schemafile.

"""
function __init__()
# Should we remove __init__ and read schema at compile time instead of runtime to make it const?
   global schema = JSONSchema.Schema(read(schemafile,String),parentFileDirectory="$moduledir/../docs")
   return nothing
end

"""
    isvalid(config)

Check that `config` is valid against Harmonie.schema
"""
isvalid(config::Dict)  = JSONSchema.isvalid(config,schema)

"""
    diagnose(config)

Check that `config` is valid against Harmonie.schema. If
valid return `nothing`, and if not, return a diagnostic String containing a
selection of one or more likely causes of failure.
"""
diagnose(config::Dict) = JSONSchema.diagnose(config,schema)

end


text:
using GrowableArrays, EllipsisNotation

using Test

const NUM_RUNS = 100
const PROBLEM_SIZE = 1000
function test1()
  u =    [1 2 3 4
          1 3 3 4
          1 5 6 3
          5 2 3 1]

  uFull = u
  for i = 1:PROBLEM_SIZE
    uFull = hcat(uFull,u)
  end
  uFull
end

function test2()
  u =    [1 2 3 4
          1 3 3 4
          1 5 6 3
          5 2 3 1]

  uFull = u

  for i = 1:PROBLEM_SIZE
    uFull = vcat(uFull,u)
  end
  uFull
end

function test3()
  u =    [1 2 3 4
          1 3 3 4
          1 5 6 3
          5 2 3 1]

  uFull = Vector{Int}(undef, 0)
  sizehint!(uFull,PROBLEM_SIZE*16)
  append!(uFull,vec(u))

  for i = 1:PROBLEM_SIZE
    append!(uFull,vec(u))
  end
  reshape(uFull,4,4,PROBLEM_SIZE+1)
  uFull
end

function test4()
  u =    [1 2 3 4
          1 3 3 4
          1 5 6 3
          5 2 3 1]

  uFull = Vector{Array{Int}}(undef, 0)
  push!(uFull,copy(u))

  for i = 1:PROBLEM_SIZE
    push!(uFull,copy(u))
  end
  uFull
end

function test5()
  u =    [1 2 3 4
          1 3 3 4
          1 5 6 3
          5 2 3 1]

  uFull = Vector{Array{Int,2}}(undef, 0)
  push!(uFull,copy(u))

  for i = 1:PROBLEM_SIZE
    push!(uFull,copy(u))
  end
  uFull
end

function test6()
  u =    [1 2 3 4
          1 3 3 4
          1 5 6 3
          5 2 3 1]

  uFull = Vector{typeof(u)}(undef, 0)
  push!(uFull,u)

  for i = 1:PROBLEM_SIZE
    push!(uFull,copy(u))
  end
  uFull
end

function test7()
  u =    [1 2 3 4
          1 3 3 4
          1 5 6 3
          5 2 3 1]

  uFull = GrowableArray(u)
  for i = 1:PROBLEM_SIZE
    push!(uFull,u)
  end
  uFull
end

function test8()
  u =    [1 2 3 4
          1 3 3 4
          1 5 6 3
          5 2 3 1]

  uFull = GrowableArray(u)
  sizehint!(uFull,PROBLEM_SIZE)
  for i = 1:PROBLEM_SIZE
    push!(uFull,u)
  end
  uFull
end

println("Run Benchmarks")
println("Pre-Compile")
#Compile Test Functions
test1()
test2()
test3()
test4()
test5()
test6()
test7()
test8()

println("Running Benchmarks")
t1 = @elapsed for i=1:NUM_RUNS test1() end
t2 = @elapsed for i=1:NUM_RUNS test2() end
t3 = @elapsed for i=1:NUM_RUNS test3() end
t4 = @elapsed for i=1:NUM_RUNS test4() end
t5 = @elapsed for i=1:NUM_RUNS test5() end
t6 = @elapsed for i=1:NUM_RUNS test6() end
t7 = @elapsed for i=1:NUM_RUNS test7() end
t8 = @elapsed for i=1:NUM_RUNS test8() end

println("Benchmark results: $t1 $t2 $t3 $t4 $t5 $t6 $t7 $t8")
println("Note the overhead of building the type before the loop for small runs")
A = [1 2; 3 4]
B = [1 2; 4 3]
G = GrowableArray(A)
push!(G,A)
push!(G,2*A)
push!(G,3*A)
@test begin
  Garr = copy(G) #Should change it to a regular array
  typeof(Garr)<:AbstractArray{Int,3}
end

# G2 = GrowableArray(1)
# push!(G2,4)
# G2[2] = 3
# G2[4,..] = 3
# @test G[4] == 12
# @test reshape(G[4,..],2,2) == B #This acts as a standard array
u =  Float64[1 2 3 4
             1 3 3 4
             1 5 6 3
             5 2 3 1]

uFull = Vector{Array{Float64,2}}(undef, 0)
push!(uFull,u)

let u = u
for i = 1:100
  u = 1.13 .* u
  push!(uFull,u)
end
end

S = StackedArray(uFull)
K = S[1,..] + S[3,..]
Sarr = copy(S)
@test typeof(Sarr)<:AbstractArray{Float64,3}

max_stars_repo_name: ChrisRackauckas/GrowableArrays

max_issues_repo_name: ChrisRackauckas/GrowableArrays

max_forks_repo_name: ChrisRackauckas/GrowableArrays


text:
using CUDA: i32

group = addgroup!(SUITE, "kernel")

dummy_kernel() = nothing
group["launch"] = @benchmarkable @cuda dummy_kernel()

wanted_threads = 10000
group["occupancy"] = @benchmarkable begin
    kernel = @cuda launch=false dummy_kernel()
    config = launch_configuration(kernel.fun)
    threads = min($wanted_threads, config.threads)
    blocks = cld($wanted_threads, threads)
end

src = CUDA.rand(Float32, 512, 1000)
dest = similar(src)
function indexing_kernel(dest, src)
    i = (blockIdx().x-1i32) * blockDim().x + threadIdx().x
    @inbounds dest[i] = src[i]
    return
end
group["indexing"] = @async_benchmarkable @cuda threads=size(src,1) blocks=size(src,2) $indexing_kernel($dest, $src)

function checked_indexing_kernel(dest, src)
    i = (blockIdx().x-1i32) * blockDim().x + threadIdx().x
    dest[i] = src[i]
    return
end
group["indexing_checked"] = @async_benchmarkable @cuda threads=size(src,1) blocks=size(src,2) $checked_indexing_kernel($dest, $src)

function rand_kernel(dest::AbstractArray{T}) where {T}
    i = (blockIdx().x-1i32) * blockDim().x + threadIdx().x
    dest[i] = rand(T)
    return
end
group["rand"] = @async_benchmarkable @cuda threads=size(src,1) blocks=size(src,2) $rand_kernel($dest)


text:
module TapirSchedulers

include("dummy_recalls.jl")

const Recalls = try
    Base.require(Base.PkgId(Base.UUID(0x30af7cf3eb4344c7afa733725b72a81e), "Recalls"))
catch
    DummyRecalls
end

include("utils.jl")
include("bounded_deque.jl")
include("bounded_scheduler.jl")
include("multi_queue.jl")
include("depth_first_scheduler.jl")

if isdefined(Base.Experimental, :Tapir)
    export ConstantPriorityTaskGroup,
        DepthFirstTaskGroup,
        RandomPriorityTaskGroup,
        WorkStealingTaskGroup,
        @sync_cp,
        @sync_df,
        @sync_rp,
        @sync_ws
    using Base.Experimental: Tapir
    include("tapir.jl")
end

function __init__()
    init_bounded_scheduler()
    init_depth_first_scheduler()
end

end


text:
import Pkg

if !isdir(joinpath(DEPOT_PATH[1], "registries", "General"))
    Pkg.Registry.add("General")
end


text:
function horner(x, p...)
    t = gensym(:t)
    ex = p[end]
    for i ∈ length(p)-1:-1:1
        ex = :(SIMDPirates.vmuladd($t, $ex, $(p[i])))
    end
    Expr(:block, :($t = $x), ex)
end

function _pirate(ex)
    postwalk(ex) do x
        # @show x
        # if @capture(x, SIMDPirates.vadd(SIMDPirates.vmul(a_, b_), c_)) || @capture(x, SIMDPirates.vadd(c_, SIMDPirates.vmul(a_, b_)))
        #     return :(SIMDPirates.vmuladd($a, $b, $c))
        # elseif @capture(x, SIMDPirates.vadd(SIMDPirates.vmul(a_, b_), SIMDPirates.vmul(c_, d_), e_)) || @capture(x, SIMDPirates.vadd(SIMDPirates.vmul(a_, b_), e_, SIMDPirates.vmul(c_, d_))) || @capture(x, SIMDPirates.vadd(e_, SIMDPirates.vmul(a_, b_), SIMDPirates.vmul(c_, d_)))
        #     return :(SIMDPirates.vmuladd($a, $b, SIMDPirates.vmuladd($c, $d, $e)))
        # elseif @capture(x, a_ += b_)
        if @capture(x, a_ += b_)
            return :($a = SIMDPirates.vadd($a, $b))
        elseif @capture(x, a_ -= b_)
            return :($a = SIMDPirates.vsub($a, $b))
        elseif @capture(x, a_ *= b_)
            return :($a = SIMDPirates.vmul($a, $b))
        elseif @capture(x, a_ /= b_)
            return :($a = SIMDPirates.vdiv($a, $b))
        elseif @capture(x, @horner a__)
            return horner(a...)
        elseif @capture(x, Base.Math.muladd(a_, b_, c_))
            return :( SIMDPirates.vmuladd($a, $b, $c) )
        elseif isa(x, Symbol) && !occursin("@", string(x))
            if x ∈ keys(VECTOR_SYMBOLS)
                return :(SIMDPirates.$(VECTOR_SYMBOLS[x]))
            else
                return x
            end
            # return get(VECTOR_SYMBOLS, x, x)
            # return :(extract_data($(get(VECTOR_SYMBOLS, x, x))))
        else
            return x
        end
    end |> esc
end

macro pirate(ex) _pirate(ex) end

struct vBitArray
    ptr::Ptr{UInt64}
end
@inline VectorizationBase.vectorizable(b::BitArray) = vBitArray(pointer(b.chunks))
@inline function Base.:+(i, b::vBitArray)
    vBitArray((i >> 3) + b.ptr)
end
@inline function Base.:+(b::vBitArray, i)
    vBitArray((i >> 3) + b.ptr)
end
@inline function Base.:-(b::vBitArray, i)
    vBitArray(b.ptr - (i >> 3))
end
function unsigned_type(N)
    if N <= 8
        return :UInt8
    elseif N <= 16
        return :UInt16
    elseif N <= 32
        return :UInt32
    elseif N <= 64
        return :UInt64
    elseif N <= 128
        return :UInt128
    else
        throw("Vectors of length $N > 128 not yet supported.")
    end
end
"""
The method on BitArrays is a hack, ignoring the type passed to Vec, to better support the LoopVectorization.jl implementation.
"""
@generated function vload(::Type{Vec{N,T}}, b::vBitArray, i::Integer) where {N,T}
    N < 8 && throw("Bit array vectors with $N < 8 not yet supported.")
    utype = unsigned_type(N)
    quote
        $(Expr(:meta, :inline))
        unsafe_load(Base.unsafe_convert(Ptr{$utype}, b.ptr + (i >> 3)))
    end
end
@generated function vload(::Type{Vec{N,T}}, b::vBitArray) where {N,T}
    N < 8 && throw("Bit array vectors with $N < 8 not yet supported.")
    utype = unsigned_type(N)
    quote
        $(Expr(:meta, :inline))
        unsafe_load(Base.unsafe_convert(Ptr{$utype}, b.ptr))
    end
end
@generated function vload(::Type{SVec{N,T}}, b::vBitArray) where {N,T}
    N < 8 && throw("Bit array vectors with $N < 8 not yet supported.")
    utype = unsigned_type(N)
    quote
        $(Expr(:meta, :inline))
        unsafe_load(Base.unsafe_convert(Ptr{$utype}, b.ptr))
    end
end
@generated function vload(::Type{SVec{N,T}}, b::vBitArray, i) where {N,T}
    N < 8 && throw("Bit array vectors with $N < 8 not yet supported.")
    utype = unsigned_type(N)
    quote
        $(Expr(:meta, :inline))
        unsafe_load(Base.unsafe_convert(Ptr{$utype}, b.ptr + (i >> 3)))
    end
end

"""
Masks on vectorizable bit arrays are currently ignored.
"""
@generated function vload(::Type{Vec{N,T}}, b::vBitArray, i::Integer, mask::Union{<:Unsigned,Vec{N,Bool}}) where {N,T}
    N < 8 && throw("Bit array vectors with $N < 8 not yet supported.")
    utype = unsigned_type(N)
    quote
        $(Expr(:meta, :inline))
        unsafe_load(Base.unsafe_convert(Ptr{$utype}, b.ptr + (i >> 3)))
    end
end
@generated function vload(::Type{Vec{N,T}}, b::vBitArray, mask::Union{<:Unsigned,Vec{N,Bool}}) where {N,T}
    N < 8 && throw("Bit array vectors with $N < 8 not yet supported.")
    utype = unsigned_type(N)
    quote
        $(Expr(:meta, :inline))
        unsafe_load(Base.unsafe_convert(Ptr{$utype}, b.ptr))
    end
end
@generated function vload(::Type{SVec{N,T}}, b::vBitArray, mask::Union{<:Unsigned,Vec{N,Bool}}) where {N,T}
    N < 8 && throw("Bit array vectors with $N < 8 not yet supported.")
    utype = unsigned_type(N)
    quote
        $(Expr(:meta, :inline))
        unsafe_load(Base.unsafe_convert(Ptr{$utype}, b.ptr))
    end
end
@generated function vload(::Type{SVec{N,T}}, b::vBitArray, i, mask::Union{<:Unsigned,Vec{N,Bool}}) where {N,T}
    N < 8 && throw("Bit array vectors with $N < 8 not yet supported.")
    utype = unsigned_type(N)
    quote
        $(Expr(:meta, :inline))
        unsafe_load(Base.unsafe_convert(Ptr{$utype}, b.ptr + (i >> 3)))
    end
end


text:
# ---
# title: 553. Optimal Division
# id: problem553
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, String
# link: <https://leetcode.com/problems/optimal-division/description/>
# hidden: true
# ---
#
# Given a list of **positive integers** , the adjacent integers will perform the
# float division. For example, [2,3,4] -> 2 / 3 / 4.
#
# However, you can add any number of parenthesis at any position to change the
# priority of operations. You should find out how to add parenthesis to get the
# **maximum** result, and return the corresponding expression in string format.
# **Your expression should NOT contain redundant parenthesis.**
#
# **Example:**
#
#
#
#     Input: [1000,100,10,2]
#     Output: "1000/(100/10/2)"
#     Explanation:
#     1000/(100/10/2) = 1000/((100/10)/2) = 200
#     However, the bold parenthesis in "1000/( **(** 100/10 **)** /2)" are redundant,
#     since they don't influence the operation priority. So you should return "1000/(100/10/2)".
#
#     Other cases:
#     1000/(100/10)/2 = 50
#     1000/(100/(10/2)) = 50
#     1000/100/10/2 = 0.5
#     1000/100/(10/2) = 2
#
#
# **Note:**
#
#   1. The length of the input array is [1, 10].
#   2. Elements in the given array will be in range [2, 1000].
#   3. There is only one optimal division for each test case.
#
#
## @lc code=start
using LeetCode

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


text:
st_fl = "3.14"
fl = parse(Float32, st_fl)
println(5fl)

st_in = "15"
in = parse(Int, st_in)
println(3in)

# c style printf, sprintf is deprecated
using Printf
@printf "fl = %0.2f\n" fl

@printf "a character: %c\n" 'α'

max_stars_repo_name: startprogress/Hello-Julia

max_issues_repo_name: startprogress/Hello-Julia

max_forks_repo_name: startprogress/Hello-Julia


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

using Documenter, QuasiGeostrophy

equations = ["Home" => "equations.md",
             "2-Layer QG" => "2_layers.md",
             "N-Layer QG" => "n_layers.md",
             ]

makedocs(
    modules = [QuasiGeostrophy],
    pages = [
        "Home" => "index.md",
        "Equations" => equations,
    ],
    sitename = "QuasiGeostrophy",
    format = Documenter.HTML(collapselevel = 1),
)

deploydocs(repo = "github.com/sandreza/QuasiGeostrophy.git")

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: upiterbarg/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
using DocStringExtensions

"""
Represents a rule scheduler for the equality saturation process

"""
abstract type AbstractScheduler end

"""
Should return `true` if the e-graph can be said to be saturated
```
cansaturate(s::AbstractScheduler)
```
"""
function cansaturate end

"""
Should return `false` if the rule `r` should be skipped
```
cansearch(s::AbstractScheduler, r::Rule)
```
"""
function cansearch end

"""
This function is called **after** pattern matching on the e-graph,
informs the scheduler about the yielded matches.
Returns `false` if the matches should not be yielded and ignored.
```
inform!(s::AbstractScheduler, iteration, r::Rule, matches)
```
"""
function inform! end

function setiter! end


text:
using Plots
a = [3,4,5,6,7,8]
plot(a)


text:
include("utils.jl")

@testset "CenteredMatrices" begin
    M = CenteredMatrices(3, 2)
    M_complex = CenteredMatrices(3, 2, ℂ)
    A = [1 2; 4 5; -5 -7]
    B = [1 2 3; 4 5 6; -5 -7 -9]    #wrong dimensions
    C = [-3 -im; 2 im; 1 0]         #complex
    D = [1 2; 3 4; 5 6]             #not centered
    types = [Matrix{Float64}]
    @testset "Real Centered Matrices Basics" begin
        @test repr(M) == "CenteredMatrices(3, 2, ℝ)"
        @test representation_size(M) == (3, 2)
        @test base_manifold(M) === M
        @test typeof(get_embedding(M)) === Euclidean{Tuple{3,2},ℝ}
        @test check_manifold_point(M, A) === nothing
        @test_throws DomainError is_manifold_point(M, B, true)
        @test_throws DomainError is_manifold_point(M, C, true)
        @test_throws DomainError is_manifold_point(M, D, true)
        @test check_tangent_vector(M, A, A) === nothing
        @test_throws DomainError is_tangent_vector(M, A, D, true)
        @test_throws DomainError is_tangent_vector(M, D, A, true)
        @test_throws DomainError is_tangent_vector(M, A, B, true)
        @test manifold_dimension(M) == 4
        @test A == project!(M, A, A)
        @test A == project(M, A, A)
        A2 = similar(A)
        embed!(M, A2, A)
        A3 = embed(M, A)
        @test A2 == A
        @test A3 == A
        E = [0 0; 1 -1; -1 1]
        F = [0.5 1; -1 -0.7; 0.5 -0.3]
        test_manifold(
            M,
            [A, E, F],
            test_injectivity_radius = false,
            test_reverse_diff = false,
            test_project_tangent = true,
            test_musical_isomorphisms = true,
            test_default_vector_transport = true,
            is_tangent_atol_multiplier = 1,
            is_point_atol_multiplier = 1,
        )
    end
    @testset "Complex Centered Matrices Basics" begin
        @test repr(M_complex) == "CenteredMatrices(3, 2, ℂ)"
        @test manifold_dimension(M_complex) == 8
        G = [1.0 1.0im; -1.0im 0.0; -1.0+1.0im -1.0im]
        H = [1.0im 0.0; -2.0im 1.0im; 1.0im -1.0im]
        test_manifold(
            M_complex,
            [C, G, H],
            test_injectivity_radius = false,
            test_reverse_diff = false,
            test_project_tangent = true,
            test_musical_isomorphisms = true,
            test_default_vector_transport = true,
            is_tangent_atol_multiplier = 1,
            is_point_atol_multiplier = 1,
        )
    end
end


text:
module Monadic
export @pure, @monadic, monadic
using Compat

"""
Simple helper type to mark pure code parts in monadic code block.

Typically you don't use this directly, but use macro `@pure` instead.
"""
struct PureCode
  code
end
"""
Mark code to contain non-monadic code.

This can be thought of as generalized version of `pure` typeclass.
"""
macro pure(expr)
  PureCode(expr)
end

# @pure is expanded within monadic
function _mergepure!(a::Int, b::Int, block::Vector{Any})
  for k ∈ a:b
    if block[k] isa PureCode
      block[k] = block[k].code
    end
  end
end

"""
    myflatmap(f, x) = Iterators.flatten(map(f, x))
    iteratorresult = @monadic map myflatmap begin
      x = 1:3
      y = [1, 6]
      @pure x + y
    end
    collect(iteratorresult)  # [2, 7, 3, 8, 4, 9]

The `@monadic` macro allows a syntax where containers and other contexts are treated rather as values, hiding the
respective well-defined side-effects.
Each line without @pure is regarded as a container, each line with @pure is treated as normal code which should be
inlined.

For the example above you see that the side-effect semantics of iterables are the same as for nested for loops. With the
crucial distinction, that the @monadic syntax has a return value.

----------------------

`@monadic` can take a third argument `wrapper` in order to first apply a custom function before executing the @monadic
code.

```
mywrapper(n::Int) = 1:n
mywrapper(any) = any
myflatmap(f, x) = Iterators.flatten(map(f, x))
iteratorresult = @monadic map myflatmap mywrapper begin
  x = 3
  y = [1, 6]
  @pure x + y
end
collect(iteratorresult)  # [2, 7, 3, 8, 4, 9]
```
"""
macro monadic(maplike, flatmaplike, expr)
  expr = macroexpand(__module__, expr)
  esc(monadic(maplike, flatmaplike, expr))
end

macro monadic(maplike, flatmaplike, wrapper, expr)
  expr = macroexpand(__module__, expr)
  esc(monadic(maplike, flatmaplike, wrapper, expr))
end

_ismonad(x) = !isa(x, Union{PureCode, LineNumberNode})  # all normal syntax should be treated as monadic

monadic(maplike, flatmaplike, expr) = monadic(maplike, flatmaplike, :identity, expr)
function monadic(maplike, flatmaplike, wrapper, expr)
  @assert(isa(expr, Expr) && expr.head == :block,
    "@monadic only supports plain :block expr, got instead $(isa(expr, Expr) ? expr.head : typeof(expr))")
  # @pure marked lines are not Expr but PureCode, hence everything which is a normal Expr is a Monad
  i = findfirst(_ismonad, expr.args)
  if isnothing(i)
    error("There need to be at least one non-@pure expression in a @monadic block")
  end
  # for everything before i we merge @pure expressions into normal code
  _mergepure!(1, i - 1, expr.args)
  _monadic(maplike, flatmaplike, wrapper, i, expr.args)
end

function _monadic(_, _, _, i::Nothing, block::Vector{Any})
  # we are checking for isnothing already beforehand
  error("this should never happen")
  # Expr(:block, block...)
end

function _monadic(maplike, flatmaplike, wrapper, i::Int, block::Vector{Any})
  wrap(expr) = wrapper === :identity ? expr : Expr(:call, wrapper, expr)

  expr = block[i]
  j = findnext(_ismonad, block, i+1)
  if isnothing(j) # last _monadic Expr is a special case
    # either i is the last entry at all, then this can be returned directly
    if i == length(block)
      block[i] = wrap(expr)
      Expr(:block, block...)
    # or this not the last entry, but @pure expressions may follow, then we construct a final fmap
    else
      _mergepure!(i + 1, length(block), block) # merge all left @pure
      lastblock = Expr(:block, block[i+1:end]...)

      callmap = if (isa(expr, Expr) && expr.head == :(=))
        subfunc = Expr(:->, Expr(:tuple, expr.args[1]), lastblock)  # we need to use :tuple wrapper to support desctructuring https://github.com/JuliaLang/julia/issues/6614
        Expr(:call, maplike, subfunc, wrap(expr.args[2]))
      else
        subfunc = Expr(:->, :_, lastblock)
        Expr(:call, maplike, subfunc, wrap(expr))
      end
      Expr(:block, block[1:i-1]..., callmap)
    end
  # if i is not the last _monadic Expr
  else
    _mergepure!(i + 1, j - 1, block) # merge all new @pure inbetween
    submonadic = _monadic(maplike, flatmaplike, wrapper, j - i, block[i+1:end])

    callflatmap = if (isa(expr, Expr) && expr.head == :(=))
      subfunc = Expr(:->, Expr(:tuple, expr.args[1]), submonadic)  # we need to use :tuple wrapper to support desctructuring https://github.com/JuliaLang/julia/issues/6614
      Expr(:call, flatmaplike, subfunc, wrap(expr.args[2]))
    else
      subfunc = Expr(:->, :_, submonadic)
      Expr(:call, flatmaplike, subfunc, wrap(expr))
    end
    Expr(:block, block[1:i-1]..., callflatmap)
  end
end

end # module


text:
struct DB
    nodes_dmp::String
    names_dmp::String
    parents::Dict{Int,Int}
    ranks::Dict{Int,Symbol}
    names::Dict{Int,String}
end

Base.show(io::IO, db::DB) = print(io, "Taxonomy.DB(\"$(db.nodes_dmp)\",\"$(db.names_dmp)\")")

"""
    Taxonomy.DB

# Constructors
```julia
DB(nodes_dmp::String, names_dmp::String)
DB(db_path::String, nodes_dmp::String, names_dmp::String)
```
Create DB(taxonomy database) object from nodes.dmp and names.dmp files.
You can specify the paths of the nodes.dmp and names.dmp files, or the directory where they exist and the names.
"""

function DB(nodes_dmp::String, names_dmp::String)
    @assert isfile(nodes_dmp)
    @assert isfile(names_dmp)

    parents, ranks = importnodes(nodes_dmp)
    namaes = importnames(names_dmp)

    return DB(nodes_dmp, names_dmp, parents, ranks, namaes)
end

function DB(db_path::String, nodes_dmp::String, names_dmp::String)
    @assert ispath(db_path)

    nodes_dmp_path = joinpath(db_path, nodes_dmp)
    names_dmp_path = joinpath(db_path, names_dmp)

    return DB(nodes_dmp_path, names_dmp_path)
end

function importnodes(nodes_dmp_path::String)
    parents = Dict{Int,Int}()
    ranks = Dict{Int,Symbol}()

    f = open(nodes_dmp_path, "r")
    for line in eachline(f)
        cols = split(line, "\t")
        taxid = parse(Int, cols[1])
        parent = parse(Int, cols[3])
        rank = Symbol(cols[5])

        parent != taxid || continue

        parents[taxid] = parent
        ranks[taxid] = rank
    end
    close(f)
    return parents, ranks
end

function importnames(names_dmp_path::String)
    namaes = Dict{Int,String}()
    f = open(names_dmp_path, "r")
    for line in eachline(f)
        cols = split(line, "\t")
        if cols[7] == "scientific name"
            taxid = parse(Int, cols[1])
            name = cols[3]
            namaes[taxid] = name
        end
    end
    close(f)
    return namaes
end


text:
export iauNutm80
"""
Form the matrix of nutation for a given date, IAU 1980 model.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2    double          TDB date (Note 1)

Returned:
   rmatn          double[3][3]    nutation matrix

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The matrix operates in the sense V(true) = rmatn * V(mean),
   where the p-vector V(true) is with respect to the true
   equatorial triad of date and the p-vector V(mean) is with
   respect to the mean equatorial triad of date.

Called:
   iauNut80     nutation, IAU 1980
   iauObl80     mean obliquity, IAU 1980
   iauNumat     form nutation matrix

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauNutm80(double date1, double date2, double rmatn[3][3])
function iauNutm80(date1::Real, date2::Real)
   # Allocate return values
   rmatn = zeros(Float64, 3, 3)

   ccall((:iauNutm80, libsofa_c), Cvoid,
        (Cdouble, Cdouble, Ptr{Cdouble},),
        convert(Float64, date1), convert(Float64, date2), rmatn)

   return SMatrix{3,3}(rmatn')
end


text:
module ColorPrinting

include("interpolation.jl")
include("parser.jl")
include("colorscheme.jl")

export @c_str, ColorScheme

end # module


text:
################################################################################
#       0. CONFIGURATION
################################################################################

# Setup initial configuration
begin
    num_workers = 4
    include("setup.jl")
end
# NOTE: first time running can take some minutes, packages installing

################################################################################
#       1. DATA GENERATION
################################################################################

# Fixed parameters
@everywhere begin
    repetitions = 3         # number of ecological networks
    vals = 2                # number of values for sweeping parameter
    N = 8                   # number of species
    M = (2^N)-1             # number of communities
end

# Select "Strength", "Connectivity", and "Universality" as example
pars = synthetic_data[1:3]
# NOTE: variable `synthetic_data` and values imported from module cNODE

# Generate GLV data
begin
    "generating synthethic data..." |> println
    generate_data(N,M,repetitions,vals,pars)
end

################################################################################
#       2. HYPERPARAMETER SEARCH
################################################################################

# Grid of values
@everywhere begin
    # Define potential hyperparameter values
    LearningRates = [[0.01,0.01],[0.01,0.05],[0.01,0.1]]
    Minibatches = [5,10]
    # Parameters for run
    max_epochs = 1000
    early_stoping = 100
    report = 10
    # Iterator over hyperparams, params and repetitions
    inx = collect(product(
                enumerate(LearningRates),
                enumerate(Minibatches),
                0:repetitions-1
            ))[:]
    # Run hyperparameter search using
    # the low value:
    v = vals[1]
    # of "Strength":
    DATA = synthetic_data[1]
    "search hyperparameters..." |> println
end

@sync @distributed for it in inx
    # Unpack hyperparameters and report
    ((j,lr),(k,mb),rep) = it
    "\t LR: $lr MB: $mb Rep: $rep" |> println
    # Load Data
    path = "./data/synthetic/$N/$DATA/P$(v)_$rep.csv"
    Z,P = import_data(path)
    ztrn,ptrn,zval,pval,ztst,ptst = split_data(Z,P)
    # Get cNODE
    cnode = getModel(N)
    # Minibatch size
    mb_size = mb > size(ztrn,2) ? size(ztrn,2) : mb
    # Train cNODE
    W, loss_train, loss_val, loss_test = train_reptile(
                                            cnode, max_epochs,
                                            mb_size, lr,
                                            ztrn, ptrn,
                                            zval, pval,
                                            ztst, ptst,
                                            report, early_stoping
                                            )
    # Save
    results = "./results/synthetic/$N/hyperparameters/"
    !ispath(results) && mkpath(results)
    writedlm(results* "train_loss_$(j)$(k)_$val$rep.csv", loss_train[end], ',')
    writedlm(results* "val_loss_$(j)$(k)_$val$rep.csv",   loss_val[end],   ',')
    writedlm(results* "test_loss_$(j)$(k)_$val$rep.csv",  loss_test[end],  ',')
end

################################################################################
#       3. HYPERPARAMETER SELECTION
################################################################################

begin
    "select hyperparameters..." |> println
    # Path to hyperparameter results
    root = "./results/synthetic/$N/hyperparameters/val_loss_"
    # Read every combination and average over repetitions
    val_hp = Array{Float64,2}(undef,length(LearningRates),length(Minibatches))
    for i in 1:length(LearningRates), j in 1:length(Minibatches)
        val =  [ readdlm(root * "$(i)$(j)_0$r.csv", ',',Float64,'\n')[end]
                    for r in 0:repetitions-1]
        val_hp[i,j] = sum(val)/length(val)
    end
    # Find minimum validation error
    best_hp = findmin(val_hp)|>last
    # Select optimum hyperparams in search
    @everywhere lr,mb = collect(product(LearningRates,Minibatches))[findmin(val_hp)|>last]
end

################################################################################
#       4. VALIDATION OF ALL DATASETS
################################################################################

@everywhere begin
    # Relative size of training dataset
    percentage = [.001,.01,0.1,0.25,0.5,0.75,1]
    # Iterator over distinct parameters, different values and repetitions
    inx = collect(product(pars,0:vals-1,0:repetitions-1))[:]
end

@sync @distributed for it in inx
    # Unpack parameters and report
    (DATA,i,j) = it
    "$DATA - \t LR: $lr MB: $mb Par: $i Rep: $j" |> println
    # Load Data
    path = "./data/synthetic/$N/$DATA/P$(i)_$(j).csv"
    Z, P = import_data(path)
    # Allocate for sync workers
    LossTrain = SharedArray{Float64}(length(percentage))
    LossVal = SharedArray{Float64}(length(percentage))
    LossTest = SharedArray{Float64}(length(percentage))
    # Results folders
    results = "./results/synthetic/$N/$DATA/"
    !ispath(results) && mkpath(results)
    !ispath(results*"loss_epochs/") && mkpath(results*"loss_epochs/")
    !ispath(results*"loss_percentage/") && mkpath(results*"loss_percentage/")
    # For different training dataset sizes
    @time for p in 1:length(percentage)
        # Load Data
        ztrn,ptrn,zval,pval,ztst,ptst = split_data(Z,P,percentage[p])
        # Get cNODE
        cnode = getModel(N)
        # Minibatches
        mb_size = mb > size(ztrn,2) ? size(ztrn,2) : mb_size
        # Train cNODE
        W, loss_train, loss_val, loss_test = train_reptile(
                                                cnode, max_epochs,
                                                mb_size, lr,
                                                ztrn, ptrn,
                                                zval, pval,
                                                ztst, ptst,
                                                report, early_stoping
                                                )
        # Save last value
        LossTrain[p] = loss_train[end]
        LossVal[p] = loss_val[end]
        LossTest[p] = loss_test[end]
        # Save last training
        if percentage[p]==percentage[end]
            writedlm(results*"loss_epochs/train_$i$j.csv", loss_train, ',')
            writedlm(results*"loss_epochs/val_$i$j.csv",   loss_val,   ',')
            writedlm(results*"loss_epochs/test_$i$j.csv",  loss_test,  ',')
        end
        "--- $(percentage[p]) ---" |> println
    end
    # Save
    writedlm(results*"loss_percentage/train_$i$j.csv", LossTrain, ',')
    writedlm(results*"loss_percentage/val_$i$j.csv",   LossVal,   ',')
    writedlm(results*"loss_percentage/test_$i$j.csv",  LossTest,  ',')
end


text:
using Random: AbstractRNG, GLOBAL_RNG

using Base:
    splat,
    Fix1

using Test: @test, @testset

import JLSO

import Dates

using Memoization: @memoize

include("InternalTestSet.jl")

ArgsDict = Dict{Symbol, Vector}

PredsAssoc = Vector{NamedTuple{(:pred, :desc, :args),
                               Tuple{Function, String, Vector{Symbol}}}}
"""
    Quickcheck

Contain a set of property to check through the generation of random input.

# Fields

- `description::AbstractString`: description for the instance.
- `rng::AbstractRNG`: PRNG used to generate inputs.
- `predicates::PredsAssoc`: predicates to check.
- `variables::ArgsDict`: Arguments used by the predicates.
- `n::Int`: Number of random inputs to generate.
- `serialize_fails::Bool`: If true, serialize failing inputs to a JLSO file.
"""
struct Quickcheck
    description::AbstractString

    ## PRNG used to generate inputs.
    rng::AbstractRNG

    ## Predicates and the names of their arguments.
    predicates::PredsAssoc

    ## Arguments, their type and a vector of randomly generated values.
    variables::ArgsDict

    ## Number of random inputs to generate.
    n::Int

    ## Should failing inputs be serialized?
    serialize_fails::Bool
end

"""
    Quickcheck(desc; rng=GLOBAL_RNG, n=100, serialize_fails=true)

Constructor for type [`Quickcheck`](@ref).

# Arguments

- `desc::AbstractString`: description for the instance.
- `rng::AbstractRNG`: PRNG used to generate inputs.
- `n::Int`: Number of random inputs to generate.
- `serialize_fails::Bool`: If true, serialize failing inputs to a JLSO file.

# Examples
```jldoctest
julia> qc = Quickcheck("A Test")
A Test: 0 predicate and 0 free variable.
```
"""
function Quickcheck(desc::AbstractString;
                    rng::AbstractRNG = GLOBAL_RNG,
                    n::Int = 100,
                    serialize_fails = true)
    Quickcheck(desc, rng, PredsAssoc(), ArgsDict(), n, serialize_fails)
end

function add_predicate(qc::Quickcheck,
                       desc::AbstractString,
                       args::Vector{Symbol},
                       types::Union{Vector{DataType},
                                    Vector{UnionAll},
                                    Vector{Union}},
                       pred::Function)
    ## Make sure that no predicate with the same description exists in
    ## `qc`.
    for (_, olddesc, _) ∈ qc.predicates
        olddesc === desc &&
            error("A predicate with the same description already exists")
    end

    for (arg, type) ∈ Iterators.zip(args, types)
        ## Make sure that suitable variables are available in `qc`. If
        ## `arg` is not present already, simply add it. Otherwise, we
        ## make sure types are the same.
        if !haskey(qc.variables, arg)
            qc.variables[arg] = generate(qc.rng, type, qc.n)
            continue
        end

        oldtype = eltype(qc.variables[arg])
        oldtype === type && continue

        error("A declaration for variable $arg already exists with type \
               $oldtype; please choose another name for $arg")
    end

    push!(qc.predicates, (pred = pred, desc = desc, args = [args...]))

    qc
end

"""
    @add_predicate qc desc pred

Add the predicate `pred` to the set of tests `qc` with description `desc`.

# Arguments

- `qc`: Object of type [`Quickcheck`](@ref Quickcheck).
- `desc`: String describing the predicate.
- `pred`: Predicate in the form of an [anonymous function](https://docs.julialang.org/en/v1/manual/functions/#man-anonymous-functions).

# Notes

The form of `pred` is very strict:
- It has to be an anonymous function. Formally, it should be an
  `Expr` of type `->`.
- The type of each argument appearing on the left-hand side of `pred` has
  to be specified with the `x::Type` syntax.
- The names of the arguments of `pred` matter! Specifically, in a given
  [`Quickcheck`](@ref Quickcheck) object, the type of every argument must
  be consistent across predicates (see examples).
- Each predicate stored in a given [`Quickcheck`](@ref Quickcheck) object
  must be given a distinct description.

# Examples
```jldoctest
julia> qc = Quickcheck("A Test")
A Test: 0 predicate and 0 free variable.

julia> @add_predicate qc "Identity" (x::Float64 -> x == x)
A Test: 1 predicate and 1 free variable.
x::Float64

julia> @add_predicate qc "Sum commute" ((n::Int, x::Float64) -> n + x == x + n)
A Test: 2 predicates and 2 free variables:
n::Int64
x::Float64

julia> @add_predicate qc "Is odd" isodd(x)
ERROR: Predicate declaration must have the form of an anonymous function (... -> ...)
[...]

julia> @add_predicate qc "Is odd" (x::Int -> is_odd(x))
ERROR: A declaration for variable x already exists with type Float64; please choose another name for x
[...]
```
"""
macro add_predicate(qc, desc, pred)
    ## Perform a bunch of checks.
    pred.head === :(->) ||
        return :(error("Predicate declaration must have the form of an \
                        anonymous function (... -> ...)"))

    local declarations = first(pred.args)

    if !(declarations isa Expr)
        var = Meta.quot(declarations)
        return :(error("A type must be specified for free variable \
                        $($var)"))
    end

    declarations.head ∈ [:(::), :tuple] ||
        return :(error("The rhs of the predicate must either be a single \
                        variable of a tuple of variables."))

    if declarations.head === :tuple
        for declaration ∈ declarations.args
            if !(declaration isa Expr)
                var = Meta.quot(declaration)
                return :(error("A type must be specified for free \
                                variable $($var)"))
            end
        end
    end

    ## Manage differences between unary and non unary predicates.
    local pred_arguments::Vector{Expr} =
        declarations.head === :tuple ?
        declarations.args : [declarations]

    local args, types = destructure_declaration(pred_arguments)

    :(add_predicate($(esc(qc)),
                    $desc,
                    [$args...],
                    [$types...],
                    $(esc(pred))))
end

function destructure_declaration(vardec::Union{Vector{Expr},
                                               Tuple{Vararg{Expr}}})
    args, types = vardec |>
        args -> Iterators.map(e -> [first(e.args), last(e.args)], args) |>
        splat(Iterators.zip)

    args, esc(Expr(:tuple, types...))
end

function Base.show(io::IO, qc::Quickcheck)
    nbpredicates = length(qc.predicates)
    nbvariables = length(qc.variables)

    predicate_string = "predicate" * (nbpredicates < 2 ? "" : 's')
    variable_string = "variable" * (nbvariables < 2 ? '.' : "s:")

    header = "$(qc.description): $nbpredicates $predicate_string and \
              $nbvariables free $variable_string"

    isempty(qc.variables) && return print(io, header)

    vars = String[]
    for (variable, entry) ∈ pairs(qc.variables)
        push!(vars, string(variable) * "::" * string(eltype(entry)))
    end

    print(io, header, "\n", join(vars, "\n"))
end

@memoize function evaluate_shrink(predicate::Function,
                                  valuation::Tuple,
                                  depth::Int)
    ## If `valuation` satisfies `predicate`, an empty tuple is returned.
    predicate(valuation...) && return valuation, zero(Int)

    ## Verify if at least one entry of valuation is shrinkable. If
    ## that's not the case, just return it.
    any(shrinkable, valuation) || return valuation, depth

    ## Let the shrinkage begin!
    candidates = valuation |>
        Fix1(Iterators.map, shrink) |>
        splat(Iterators.product)

    isempty(candidates) && return valuation, depth

    best_candidate, candidates = Iterators.peel(candidates)
    shrunk_best_candidates = evaluate_shrink(predicate, best_candidate, depth + 1)
    for candidate ∈ candidates
        shrunk_candidate = evaluate_shrink(predicate, candidate, depth + 1)

        if last(shrunk_candidate) > last(shrunk_best_candidates)
            shrunk_best_candidates = shrunk_candidate
        end
    end

    ## If no candidate falsifies the predicate, just return the
    ## original valuation.
    iszero(last(shrunk_best_candidates)) && return valuation, depth

    shrunk_best_candidates
end

evaluate_shrink(predicate::Function, valuation::Tuple) =
    first(evaluate_shrink(predicate, valuation, 1))

function quickcheck(qc::Quickcheck, file_id::AbstractString)
    ## Stores failed cases. The key is the test predicate's
    ## description (as a symbol). The value is a named tuple
    ## containing the actual predicate (predicate) and a vector of
    ## inputs for which the predicate evaluated to `false`
    ## (valuation).
    local failed = Dict{Symbol,
                        NamedTuple{(:predicate, :valuations),
                                   Tuple{Function, Vector{Tuple}}}}()

    @testset InternalTestSet "Test $desc" for
        (pred, desc, args) ∈ qc.predicates

        ## Flip to `false` if predicate evaluates to `false` for any
        ## valuation.
        local holds = true

        ## Special cases.
        specialcases_itr = args |>
            Fix1(Iterators.map,
                 specialcases ∘ eltype ∘ Fix1(getindex, qc.variables)) |>
            splat(Iterators.product)

        ## Random cases.
        randomcases_itr = args |>
            Fix1(Iterators.map, Fix1(getindex, qc.variables)) |>
            splat(Iterators.zip)

        for valuation ∈ Iterators.flatten((specialcases_itr,
                                           randomcases_itr))
            if !pred(valuation...)
                try
                    valuation = evaluate_shrink(pred, valuation)
                catch e
                    @warn "Could not shrink valuation, got the following error \
                           $e"
                end

                ex = Expr(:tuple,
                          (Expr(:(=), arg, val)
                           for (arg, val) ∈ zip(args, valuation))...)
                @warn "Predicate \"$desc\" does not hold \
                               for valuation $ex"

                if qc.serialize_fails
                    if !haskey(failed, Symbol(desc))
                        failed[Symbol(desc)] = (predicate = pred,
                                                valuations = Tuple[])
                    end

                    push!(failed[Symbol(desc)].valuations, valuation)
                end
                holds = false
            end
        end
        @test holds
    end

    if !isempty(failed) && qc.serialize_fails
        filename = "JCheck_" * file_id
        fileextension = ".jchk"

        ## Make sure that no 2 files with the same name are
        ## created. Most likely useless except for people with
        ## nothing better to do trying to break shit.
        if ispath(filename * fileextension)
            filename *= "--1"
        end

        while ispath(filename * fileextension)
            id_start = last(findfirst("--", filename)) + one(Int)

            idx = parse(Int, filename[id_start:end]) + 1

            filename = filename[begin:(id_start - 1)] * string(idx)
        end

        JLSO.save(filename * fileextension, failed)

        ## Warn the user that some predicates do not hold.
        message_color = Base.warn_color()

        println()
        printstyled("Some predicates do not hold for some valuations; \
                             they have been saved to \
                             $(filename * fileextension). Use function ",
                    color = message_color)
        printstyled("load", color = message_color, underline = true)
        printstyled(" and macro ", color = message_color)
        printstyled("@getcases",
                    color = message_color,
                    underline = true)
        printstyled(" to explore problematic cases.",
                    color = message_color)
        println("\n")
    end

    nothing
end

"""
    @quickcheck qc [file_id::AbstractString="yyyy-mm-dd_HH-MM-SS"]

Check the properties specified in object `qc` of type [`Quickcheck`](@ref).

If `qc.serialize_fails` is `true`, serialize the failing cases to
`JCheck_<file_id>.jchk`. Those can latter be analyzed using
[`load`](@ref) and [`@getcases`](@ref).

# Note

If no argument `file_id` is passed, defaults to current time.
"""
macro quickcheck(qc, file_id::AbstractString)
    _qc = esc(qc)

    :(quickcheck($_qc, $file_id))
end

macro quickcheck(qc)
    _qc = esc(qc)

    file_id = Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS")

    :(@quickcheck($_qc, $file_id))
end