text:
using Pkg
Pkg.activate("docs")
Pkg.develop(PackageSpec(path=pwd()))
Pkg.add(url="https://github.com/JuliaNLSolvers/LsqFit.jl", rev="master") # I implemented a bugfix (PR #182) that is merged, but not yet released. This line can be removed once a new version is released.
Pkg.instantiate()

using MRIgeneralizedBloch
using Documenter
using Literate
using Plots # to not capture precompilation output

# HTML Plotting Functionality
struct HTMLPlot
    p # :: Plots.Plot
end
const ROOT_DIR = joinpath(@__DIR__, "build")
const PLOT_DIR = joinpath(ROOT_DIR, "plots")
function Base.show(io::IO, ::MIME"text/html", p::HTMLPlot)
    mkpath(PLOT_DIR)
    path = joinpath(PLOT_DIR, string(UInt32(floor(rand()*1e9)), ".html"))
    Plots.savefig(p.p, path)
    if get(ENV, "CI", "false") == "true" # for prettyurl
        print(io, "<object type=\"text/html\" data=\"../../$(relpath(path, ROOT_DIR))\" style=\"width:100%;height:425px;\"></object>")
    else
        print(io, "<object type=\"text/html\" data=\"../$(relpath(path, ROOT_DIR))\" style=\"width:100%;height:425px;\"></object>")
    end
end

# Notebook hack to display inline math correctly
function notebook_filter(str)
    re = r"(?<!`)``(?!`)"  # Two backquotes not preceded by nor followed by another
    replace(str, re => "\$")
end

# Literate
OUTPUT = joinpath(@__DIR__, "src/build_literate")

files = [
    "tutorial.jl",
    "Greens_functions.jl",
    "Simulation_ContinuousWave.jl",
    "Simulation_Pulse.jl",
    "Analyze_NMR_IR_Data.jl",
    "Analyze_NMR_PreSat_Data.jl",
    "Linear_Approximation.jl",
]

for file in files
    file_path = joinpath(@__DIR__, "src/", file)
    Literate.markdown(file_path, OUTPUT)
    Literate.notebook(file_path, OUTPUT, preprocess=notebook_filter; execute=false)
    Literate.script(  file_path, OUTPUT)
end

DocMeta.setdocmeta!(MRIgeneralizedBloch, :DocTestSetup, :(using MRIgeneralizedBloch); recursive=true)

makedocs(;
    doctest = false,
    modules=[MRIgeneralizedBloch],
    authors="Jakob Asslaender <jakob.asslaender@nyumc.org> and contributors",
    repo="https://github.com/JakobAsslaender/MRIgeneralizedBloch.jl/blob/{commit}{path}#{line}",
    sitename="MRIgeneralizedBloch.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JakobAsslaender.github.io/MRIgeneralizedBloch.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Quick Start Tutorial" => "build_literate/tutorial.md",
        "Generalized Bloch Paper" => Any[
        "build_literate/Greens_functions.md",
        "build_literate/Simulation_ContinuousWave.md",
        "build_literate/Simulation_Pulse.md",
        "build_literate/Analyze_NMR_IR_Data.md",
        "build_literate/Analyze_NMR_PreSat_Data.md",
        "build_literate/Linear_Approximation.md",
        ],
        "API" => "api.md",
    ],
)

# Set dark theme as default independent of the OS's settings
run(`sed -i'.old' 's/var darkPreference = false/var darkPreference = true/g' docs/build/assets/themeswap.js`)

deploydocs(;
    repo="github.com/JakobAsslaender/MRIgeneralizedBloch.jl",
)


text:
using Random
using Test
using BiobakeryUtils
using BiobakeryUtils.Conda
using SparseArrays
using DelimitedFiles
using CSV

isdir(Conda.bin_dir(:BiobakeryUtils)) || BiobakeryUtils.install_deps()
ENV["PATH"] = ENV["PATH"] * ":" * Conda.bin_dir(:BiobakeryUtils)

@testset "CLI" begin
    @testset "Utilities" begin
        cmd = ["thing", "foo_bar"]
        cmd2 = copy(cmd)

        BiobakeryUtils.add_cli_kwargs!(cmd, Dict(:some_thing=> "foo", :bool=> true))
        @test all(cmd .== ["thing", "foo_bar", "--some_thing", "foo", "--bool"])
        BiobakeryUtils.add_cli_kwargs!(cmd2, Dict(:some_thing=> "foo", :bool=> true); optunderscores=false)
        @test all(cmd2 .== ["thing", "foo_bar", "--some-thing", "foo", "--bool"])
    end

    @testset "Metaphlan" begin
        @test BiobakeryUtils.check_for_install("metaphlan") |> isnothing
        @test BiobakeryUtils.check_for_install("merge_metaphlan_tables.py") |> isnothing

        @test metaphlan("", ""; help=true).exitcode == 0

        profiles = filter(f-> contains(f, "_profile.tsv"), readdir("files/metaphlan", join=true))
        @test metaphlan_merge(profiles, "files/metaphlan/merged_abundance_table.tsv").exitcode == 0
    end

    @testset "Humann" begin
        @test BiobakeryUtils.check_for_install("humann") |> isnothing
        @test BiobakeryUtils.check_for_install("humann_rename_table") |> isnothing
        @test BiobakeryUtils.check_for_install("humann_renorm_table") |> isnothing
        @test BiobakeryUtils.check_for_install("humann_join_tables") |> isnothing
        @test BiobakeryUtils.check_for_install("humann") |> isnothing
        @test humann("", ""; help=true).exitcode == 0

    end
end

@testset "Metaphlan" begin
    profile_1 = metaphlan_profile("files/metaphlan/SRS014464-Anterior_nares_profile.tsv";  sample="SRS014464")
    @test first(abundances(profile_1["Bacteria", "SRS014464"])) == 100.0
    @test first(abundances(profile_1["Pseudomonadales", "SRS014464"])) == 97.28734
    @test size(profile_1) == (13, 1)
    profile_2 = metaphlan_profile("files/metaphlan/SRS014459-Stool_profile.tsv", 3)
    @test size(profile_2) == (2, 1)
    @test first(abundances(profile_2["Firmicutes", "SRS014459-Stool_profile"])) == 68.90167
    profile_3 = metaphlan_profile("files/metaphlan/SRS014464-Anterior_nares_profile.tsv", :phylum)
    @test size(profile_3) == (2, 1)
    @test first(abundances(profile_3["Proteobacteria", 1])) == 97.28734

    merge_profile_1 = metaphlan_profiles("files/metaphlan/merged_abundance_table.tsv"; samplestart=3)
    @test size(merge_profile_1) == (62, 6)
    @test first(abundances(merge_profile_1["Moraxella", 5])) == 97.28734
    merge_profile_2 = metaphlan_profiles("files/metaphlan/merged_abundance_table.tsv", :family; samplestart=3)
    @test size(merge_profile_2) == (13, 6)
    @test first(abundances(merge_profile_2["Micrococcaceae", "SRS014464-Anterior_nares"])) == 0.0
    merge_profile_3 = metaphlan_profiles("files/metaphlan/merged_abundance_table.tsv", 7; samplestart=3)
    @test size(merge_profile_3) == (16, 6)
    @test first(abundances(merge_profile_3["Haemophilus_haemolyticus", 3])) == 1.35528
    CSV.write("files/metaphlan/merged_abundance_table2.csv", merge_profile_1)

    profiles = filter(f-> contains(f, "_profile.tsv"), readdir("files/metaphlan", join=true))
    @test_throws ArgumentError metaphlan_profiles(profiles; samples = ["sample1"])
    multi_profile_1 = metaphlan_profiles(profiles; samples=["sample$i" for i in 1:length(profiles)])
    @test abundances(multi_profile_1) == abundances(metaphlan_profiles(profiles))
    @test size(multi_profile_1) == (62, 6)
    @test first(abundances(multi_profile_1["Firmicutes", "sample1"])) == 68.90167
    multi_profile_2 = metaphlan_profiles(profiles, 3; samples=["sample$i" for i in 1:length(profiles)])
    @test abundances(multi_profile_2) == abundances(metaphlan_profiles(profiles, :class))
    @test size(multi_profile_2) == (6,6)
    @test first(abundances(multi_profile_2["Bacteroidia", "sample1"])) == 31.09833

    taxstring = "k__Archaea|p__Euryarchaeota|c__Methanobacteria|o__Methanobacteriales|f__Methanobacteriaceae|g__Methanobrevibacter|s__Methanobrevibacter_smithii"
    taxa = parsetaxa(taxstring)
    @test length(taxa) == 7
    @test parsetaxon(taxstring, 1) == Taxon("Archaea", :kingdom)
    @test parsetaxon(taxstring, :family) == Taxon("Methanobacteriaceae", :family)
    @test parsetaxon(taxstring) == Taxon("Methanobrevibacter_smithii", :species)
    @test_throws ArgumentError parsetaxon(taxstring, 8)

    @test parsetaxon("k__Archaea|p__Euryarchaeota|c__Methanobacteria", 2) == Taxon("Euryarchaeota", :phylum)
    @test parsetaxon("k__Archaea|p__Euryarchaeota|c__Methanobacteria") == Taxon("Methanobacteria", :class)
end

@testset "HUMAnN" begin
    p1 = humann_profile("files/humann/single_1.tsv")
    p2 = humann_profile("files/humann/single_2.tsv")
    @test p1 isa CommunityProfile
    @test size(p1) == (560, 1)
    @test samplenames(p1) == ["single_1"]
    @test samplenames(humann_profile("files/humann/single_1.tsv"; sample = "sample1")) == ["sample1"]
    @test samplenames(humann_profile("files/humann/single_1.tsv"; sample = MicrobiomeSample("sample1"))) == ["sample1"]

    @test all(f-> !hastaxon(f), features(p1)) # unstratified
    @test all(f-> !occursin('|', name(f)), features(p1))

    pj = humann_profiles("files/humann/joined.tsv")
    @test size(pj) == (560, 2)
    @test isempty(setdiff(features(pj), features(commjoin(p1, p2))))
    @test samplenames(pj) == samplenames(commjoin(p1, p2))

    pj_strat = humann_profiles("files/humann/joined.tsv"; stratified = true)
    @test size(pj_strat) == (1358, 2)
    @test !isempty(setdiff(features(pj_strat), features(pj)))
    @test isempty(setdiff(featurenames(pj_strat), featurenames(pj)))
    @test isempty(setdiff(features(filter(!hastaxon, pj_strat)), features(pj)))
    CSV.write("files/humann/joined_roundtrip.tsv", pj_strat; delim='\t')
end

ext: jl


text:
@recipe(TeXImg, tex) do scene
    merge(
        default_theme(scene),
        Attributes(
            color = Makie.automatic,
            implant = true,
            dpi = 72.0,
            align = (:center, :center),
            scale = 1.0,
            position = Point3{Float32}(0),
            rotation = 0f0,
        )
    )
end

function Makie.convert_arguments(::Type{<: TeXImg}, x::AbstractString)
    return (CachedTeX(implant_text(x)),)
end

function Makie.convert_arguments(::Type{<: TeXImg}, x::LaTeXString)
    if first(x) == "\$" && last(x) == "\$"
        return (CachedTeX(implant_math(x)),)
    else
        return (CachedTeX(implant_text(x)),)
    end
end

function Makie.convert_arguments(::Type{<: TeXImg}, doc::TeXDocument)
    return (CachedTeX(doc),)
end

function Makie.plot!(plot::T) where T <: TeXImg
    # image to draw
    # We always want to draw this at a 1:1 ratio, so increasing textsize or
    # changing dpi should rerender
    img = map(plot[1], plot.dpi, plot.scale) do texdoc, dpi, scale
        scale = scale
        rsvg2img(texdoc.handle, scale * dpi)
    end

    # Rect to draw in
    # This is mostly aligning
    xr = Observable(0.0..1.0)
    yr = Observable(0.0..1.0)
    lift(img, plot.position, plot.align) do img, pos, align
        halign, valign = align
        x, y = pos
        w, h = size(img)

        if halign == :left
            x -= 0
        elseif halign == :center
            x -= w / 2
        elseif halign == :right
            x -= w
        end

        if valign == :top
            y -= h
        elseif valign == :center
            y -= h / 2
        elseif valign == :bottom
            y -= 0
        end

        xr[] = x..x+w
        yr[] = y..y+h
        nothing
    end

    model = map(plot.model, plot.rotation, xr, yr) do model, angle, xr, yr
        x0 = xr.left; x1 = xr.right
        y0 = yr.left; y1 = yr.right
        model *
        Makie.translationmatrix(Vec3{Float32}(0.5(x1+x0), 0.5(y1+y0), 0)) *
        Makie.rotationmatrix_z(angle) *
        Makie.translationmatrix(- Vec3{Float32}(0.5(x1+x0), 0.5(y1+y0), 0))
    end

    image!(plot, xr, yr, img, model=model)
end

function get_ink_extents(surf::CairoSurface)
    dims = zeros(Float64, 4)

    ccall(
        (:cairo_recording_surface_ink_extents, CairoMakie.Cairo.libcairo),
        Cvoid,
        (Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
        surf.ptr, Ref(dims, 1), Ref(dims, 2), Ref(dims, 3), Ref(dims, 4)
    )

    return dims
end

function CairoMakie.draw_plot(scene::Scene, screen::CairoMakie.CairoScreen, img::T) where T <: MakieTeX.TeXImg

    bbox = Makie.boundingbox(img)

    ctx = screen.context
    tex = img[1][]
    halign, valign = img.align[]

    handle = svg2rsvg(tex.svg)
    dims = tex.raw_dims

    surf, rctx = rsvg2recordsurf(handle)
    x0, y0, w, h = get_ink_extents(surf)

    pos = CairoMakie.project_position(scene, :data, img.position[], img.model[])
    scale = img.scale[]
    _w = scale * w; _h = scale * h
    scale_factor = CairoMakie.project_scale(scene, :data, Vec2{Float32}(_w, _h), img.model[])

    pos = if halign == :left
        pos
    elseif halign == :center
        pos .- (scale_factor[1] / 2, 0)
    elseif halign == :right
        pos .- (scale_factor[1], 0)
    end

    pos = if valign == :top
        pos .+ (0, scale_factor[2])
    elseif valign == :center
        pos .+ (0, scale_factor[2] / 2)
    elseif valign == :bottom
        pos
    end

    Cairo.save(ctx)
    Cairo.translate(
        ctx,
        pos[1],
        pos[2] - (h + y0) * scale_factor[2] / h
    )
    Cairo.rotate(ctx, -img.rotation[])
    # Rotated center - normal center
    cx = 0.5scale_factor[1] * cos(img.rotation[]) - 0.5scale_factor[2] * sin(img.rotation[]) - 0.5scale_factor[1]
    cy = 0.5scale_factor[1] * sin(img.rotation[]) + 0.5scale_factor[2] * cos(img.rotation[]) - 0.5scale_factor[2]
    Cairo.translate(ctx, cx, cy)
    Cairo.scale(
        ctx,
        scale_factor[1] / w,
        scale_factor[2] / h
    )

    render_surface(ctx, surf)
    Cairo.restore(ctx)
end


text:
using GNSSTools
using Statistics
using PyPlot
pygui(true)

N = 10
# idx = (rand(1:N), rand(1:N))
idx = (4, 4)
background_val = 1
peak_val = 5
corr = fill(background_val, N, N)
corr[idx[1],idx[2]] = peak_val

# Get codes for plot
unknown_code = circshift(deepcopy(l1ca_codes[1]), 500)
prn_1_code = deepcopy(l1ca_codes[1])
prn_2_code = deepcopy(l1ca_codes[2])
# Make zeros -1
A = 2
unknown_code[unknown_code .== 0] .= -A
prn_1_code[prn_1_code .== 0] .= -A
prn_2_code[prn_2_code .== 0] .= -A
unknown_code[unknown_code .== 1] .= A
prn_1_code[prn_1_code .== 1] .= A
prn_2_code[prn_2_code .== 1] .= A
# Perform auto and cross correlations
auto_correlation = real.(fft_correlate(unknown_code, prn_1_code))
cross_correlation = real.(fft_correlate(unknown_code, prn_2_code))

xlims = [-50, length(auto_correlation)+50]
ylims = [mean(cross_correlation)-120, maximum(auto_correlation)+53]

# # Simulate L1 C/A signal and perform course acquisition
# signal = definesignal(define_l1ca_code_type(), 5e6, 1; prn=1, f_d=1000,
#                       code_start_idx=1000, include_phase_noise=true)
# generatesignal!(signal)
# replica = definesignal(define_l1ca_code_type(), 5e6, 1e-3)
# fd_est, n0_est, SNR_est, corr_result = courseacquisition(signal, replica, 1;
#                                                          return_corrresult=true,
#                                                          Δfd=1000)
# corr_crop = corr_result[:,n0_est-5:n0_est+5]
# idx_crop = argmax(corr_crop)
# box_x = [idx_crop[2]-0.5, idx_crop[2]-0.5, idx_crop[2]+0.5, idx_crop[2]+0.5, idx_crop[2]-0.5] .- 1
# box_y = [idx_crop[1]+0.5, idx_crop[1]-0.5, idx_crop[1]-0.5, idx_crop[1]+0.5, idx_crop[1]+0.5] .- 1

title_font_size = 9
axis_label_font_size = 8
ytick_vals = [0, 1023]
xtick_vals = [1, 501, 1023]
xtick_labels = [string(i) for i in [524, 1, 523]]

fig = figure(figsize=(6.5,1.5))
ax1 = fig.add_subplot(1,2,1)
ax1.plot(auto_correlation, "k")
xlim(xlims)
ylim(ylims)
ax1.set_yticks(ytick_vals)
ax1.set_xticks(xtick_vals)
ax1.set_xticklabels(xtick_labels)
yticks(fontsize=axis_label_font_size)
xticks(fontsize=axis_label_font_size)
xlabel("Chips", fontsize=axis_label_font_size)
ylabel("rⁱ[n]  ⋆  sⁱ[n]", fontsize=axis_label_font_size)
title("(a)", fontsize=title_font_size)
# xlabel("Time", fontsize=axis_label_font_size)
tick_params(left=true, bottom=true, labelleft=true, labelbottom=true)
ax2 = fig.add_subplot(1,2,2)
ax2.plot(cross_correlation, "k")
xlim(xlims)
ylim(ylims)
ax2.set_yticks(ytick_vals)
ax2.set_xticks(xtick_vals)
ax2.set_xticklabels(xtick_labels)
yticks(fontsize=axis_label_font_size)
xticks(fontsize=axis_label_font_size)
xlabel("Chips", fontsize=axis_label_font_size)
ylabel("rⁱ[n]  ⋆  sʲ[n]", fontsize=axis_label_font_size)
title("(b)", fontsize=title_font_size)
tick_params(left=true, bottom=true, labelleft=true, labelbottom=true)
# tick_params(left=true, bottom=false, labelleft=true, labelbottom=false)
# tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
subplots_adjust(hspace=0.5, wspace=0.4, top=0.85, left=0.1, right=0.93, bottom=0.3)
savefig("figures/ch3_autocorrelation_vs_cross_correlation.svg", dpi=300)

# Simulate L1 C/A signal and perform course acquisition
signal = definesignal(define_l1ca_code_type(), 5e6, 1; prn=1, f_d=1000,
                      code_start_idx=1000, include_phase_noise=true)
generatesignal!(signal)
replica = definesignal(define_l1ca_code_type(), 5e6, 1; include_phase_noise=false,
                       include_thermal_noise=false, include_adc=false,
                       )


text:
using D3Trees
using Colors
using Printf

"""
Return text to display below the node corresponding to state or action s
"""
node_tag(s) = string(s)

"""
Return text to display in the tooltip for the node corresponding to state or action s
"""
tooltip_tag(s) = node_tag(s)

"""
Creates a D3Tree instance for a given tree.

Use it as
  a, info = action_info(planner, state)
  D3Tree(info[:tree])
"""
function D3Trees.D3Tree(tree::PISTree; title="MCTS-PIS Tree", kwargs...)
    lens = length(tree.state_nodes)
    lensa = length(tree.state_action_nodes)
    len = lens + lensa
    children = Vector{Vector{Int}}(undef, len)
    text = Vector{String}(undef, len)
    tt = fill("", len)
    style = fill("", len)
    link_style = fill("", len)
    max_q = maximum(node.q for (k, node) in tree.state_action_nodes)
    min_q = minimum(node.q for (k, node) in tree.state_action_nodes)

    snode_dict = Dict((snode.id, snode) for (_, snode) in tree.state_nodes)
    for (id, snode) in sort(collect(snode_dict), by=x->x[1])
        s_label = snode.s_label

        children_ids = sizehint!(Int64[], length(snode.children))
        for a in snode.children
            sanode = tree.state_action_nodes[(s_label, a)]
            push!(children_ids, sanode.id)

            w = 20.0 * sqrt(sanode.n / snode.total_n)
            link_style[sanode.id + lens] = "stroke-width:$(w)px"
        end

        children[id] = children_ids
        text[id] = Printf.@sprintf("""
                                   %25s
                                   N: %6d
                                   """,
                                   node_tag(s_label),
                                   snode.total_n)
        tt[id] = """
                 $(tooltip_tag(s_label))
                 N: $(snode.total_n)
                 """
    end
    sanode_dict = Dict((sanode.id, sanode) for (_, sanode) in tree.state_action_nodes)
    for (id, sanode) in sort(collect(sanode_dict), by=x->x[1])
        a_label = sanode.a_label

        children[id + lens] = let
            sp_children = if !isempty(sanode.unique_transitions)
                sanode.unique_transitions
            else
                first.(sanode.transitions)
            end
            @assert length(sp_children) == sanode.n_a_children
            collect(tree.state_nodes[sp].id for sp in sp_children)
        end

        text[id + lens] = @sprintf("""
                                   %25s
                                   Q: %6.2f
                                   N: %6d
                                   """,
                                   node_tag(a_label),
                                   sanode.q,
                                   sanode.n)
        tt[id + lens] = """
                        $(tooltip_tag(a_label))
                        Q: $(sanode.q)
                        N: $(sanode.n)
                        """

        rel_q = (sanode.q - min_q) / (max_q - min_q)
        if isnan(rel_q)
            color = colorant"gray"
        else
            color = weighted_color_mean(rel_q, colorant"green", colorant"red")
        end
        style[id + lens] = "stroke:#$(hex(color))"
    end

    return D3Tree(children;
                  text=text,
                  tooltip=tt,
                  style=style,
                  link_style=link_style,
                  title=title,
                  kwargs...)
end


text:
# Number of points to generate.
N = 80
M = round(Int, N / 4)
Random.seed!(1234)

# Generate artificial data.
x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
xt1s = Array([[x1s[i] + 0.5; x2s[i] + 0.5] for i = 1:M])
x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
append!(xt1s, Array([[x1s[i] - 5; x2s[i] - 5] for i = 1:M]))

x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
xt0s = Array([[x1s[i] + 0.5; x2s[i] - 5] for i = 1:M])
x1s = rand(M) * 4.5; x2s = rand(M) * 4.5;
append!(xt0s, Array([[x1s[i] - 5; x2s[i] + 0.5] for i = 1:M]))

# Store all the data for later.
xs = [xt1s; xt0s]
ts = [ones(2*M); zeros(2*M)]

X = hcat(xs...)
Y = ts#onehotbatch(ts, 0:1)

# Plot data points.
function plot_data()
    x1 = map(e -> e[1], xt1s)
    y1 = map(e -> e[2], xt1s)
    x2 = map(e -> e[1], xt0s)
    y2 = map(e -> e[2], xt0s)

    Plots.scatter(x1,y1, color="red", clim = (0,1))
    Plots.scatter!(x2, y2, color="blue", clim = (0,1))
end


text:
using Documenter
using EzXML

makedocs(
    format=:html,
    sitename="EzXML.jl",
    modules=[EzXML],
    pages=["index.md", "manual.md", "reference.md", "devnotes.md"])

deploydocs(
    repo="github.com/bicycle1885/EzXML.jl.git",
    julia="0.7",
    target="build",
    deps=nothing,
    make=nothing)


text:
# Adjoint test for F and J
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall, PyPlot, JUDI.TimeModeling, Test, LinearAlgebra

## Set up model structure
n = (160, 170)	# (x,y,z) or (x,z)
d = (10.,10.)
o = (0.,0.)

# Velocity [km/s]
v = ones(Float32,n) .* 2.0f0
v[:,Int(round(end/3)):end] .= 4f0
v0 = smooth10(v,n)
rho = ones(Float32, n)
rho[:, Int(round(end/2)):end] .= 1.5f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = (1f0 ./ v0).^2
dm = vec(m - m0)

# Setup info and model structure
nsrc = 1
model = Model(n,d,o,m,rho=rho)
model0 = Model(n,d,o,m0,rho=rho)

## Set up receiver geometry
nxrec = 141
xrec = range(600f0,stop=1000f0,length=nxrec)
yrec = 0f0
zrec = range(100f0,stop=100f0,length=nxrec)

# receiver sampling and recording time
timeR = 800f0	# receiver recording time [ms]
dtR = calculate_dt(n,d,o,v,rho)    # receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR,nsrc=nsrc)

## Set up source geometry (cell array with source locations for each shot)
xsrc = 800f0
ysrc = 0f0
zsrc = 50f0

# source sampling and number of time steps
timeS = 800f0
dtS = calculate_dt(n,d,o,v,rho) # receiver sampling interval

# Set up source structure
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

# Info structure
ntComp = get_computational_nt(srcGeometry,recGeometry,model0)
info = Info(prod(n), nsrc, ntComp)

# setup wavelet
f0 = 0.008f0
wavelet = ricker_wavelet(timeS,dtS,f0)
wave_rand = wavelet.*rand(Float32,size(wavelet))

###################################################################################################

# Modeling operators
opt = Options(sum_padding=true)
F = judiModeling(info,model0,srcGeometry,recGeometry; options=opt)
q = judiVector(srcGeometry,wavelet)

# Nonlinear modeling
d_hat = F*q

# Generate random noise data vector with size of d_hat in the range of F
qr = judiVector(srcGeometry,wave_rand)
d1 = F*qr

# Adjoint computation
q_hat = adjoint(F)*d1

# Result F
a = dot(d1, d_hat)
b = dot(q, q_hat)
@test isapprox(a/b - 1, 0, atol=1f-4)

# Linearized modeling
J = judiJacobian(F,q)

dD_hat = J*dm
dm_hat = adjoint(J)*d_hat

c = dot(dD_hat, d_hat)
d = dot(dm, dm_hat)
@test isapprox(c/d - 1, 0, atol=1f-4)


text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random
path1=realpath(dirname(@__FILE__)*"/..")
include(string(path1,"/src/ANNBN.jl"))

##### PARTIAL DIFFERENTIAL EQUATIONS
# Laplace's Equation http://www.robots.ox.ac.uk/~jmb/lectures/pdelecture5.pdf
a=1;b=1;n=50;dx=a/n;x0=zeros(n); x=Array{Float64}(undef,0); for i=1:n-1 global x0,x,dx; x=[x;x0.+(i)*dx] end; x
ny=50;dy=b/n;y=convert(Array{Float64},range(dy,b-dy,length=ny)); y=repeat(y,n-1)
xx_train=[x y];i_train=size(xx_train,1);vars=2;
# boundaries
nof_b=20;xb1=zeros(nof_b,2);xb1[:,2]=range(0,1,length=nof_b);xb2=zeros(nof_b,2);xb2[:,1]=range(0,1,length=nof_b);
xb3=ones(nof_b,2);xb3[:,2]=range(0,1,length=nof_b);xb4=ones(nof_b,2);xb4[:,1]=range(0,1,length=nof_b)
xb=[xb1;xb2;xb3;xb4];yb=zeros(size(xb,1));yb[end-nof_b+1:end].=sin.(pi.*range(0,1,length=nof_b))
# source
y_source=zeros(i_train).+0.0(rand(i_train)) # you may adjust the level of noise here

##### Clustering
neurons=50;inds_all,n_per_part=ANNBN.___clustering(neurons,xx_train,200)

##### Solve PDE
# Use # 1 Gaussian for calc_phi,calc_phi_deriv,calc_phi_deriv_deriv
cc1=1.0
a_all,a_layer1,layer1=ANNBN.train_layer_1_rbf_pde(neurons,vars,i_train,n_per_part,inds_all,xx_train,cc1,xb,yb,y_source);
predl1,layer1_train=ANNBN.predict_new_rbf(a_all,a_layer1,xx_train,i_train,neurons,n_per_part,inds_all,xx_train,vars,cc1)
scatter3d(xx_train[1:1:end,1],xx_train[1:1:end,2],predl1[1:1:end],label="ANNBN solution",
    color=:black,markershape=:xcross,legend=:topleft,markersize=1)
# Analytical solution
an_sol=sinh.(pi*xx_train[:,2]./a).*sin.(pi*xx_train[:,1]./a)./sinh(pi*b/a)
scatter3d!(xx_train[:,1],xx_train[:,2],an_sol,label="Analytical solution",
    color=:black,markershape=:star5,legend=:topleft,markersize=1)
# savefig("Laplace.pdf")
err=an_sol .- predl1;mean(abs.(err))
# scatter3d(xx_train[:,1],xx_train[:,2],err,color=:black,markershape=:xcross,legend=false,markersize=1)
# savefig("Laplace_noise_err.pdf")

# Calculate Derivatives
predl1_deriv_deriv1=ANNBN.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,1)
predl1_deriv_deriv2=ANNBN.predict_new_rbf_deriv_deriv(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1,2)
err_2nd_derivatives=predl1_deriv_deriv1+predl1_deriv_deriv2.-0.0
mean(abs.(err_2nd_derivatives))
scatter3d(xx_train[:,1],xx_train[:,2],predl1_deriv_deriv1,label="Derivative fxx",
    color=:black,markershape=:xcross,legend=:topleft,markersize=1,camera=(70, 40)
    ,xlabel="x",ylabel="y")
scatter3d!(xx_train[:,1],xx_train[:,2],predl1_deriv_deriv2,label="Derivative fyy",
    color=:black,markershape=:star5,legend=:topleft,markersize=1)
scatter3d!(xx_train[:,1],xx_train[:,2],predl1_deriv_deriv1+predl1_deriv_deriv2,label="fxx + fyy",
    color=:black,markershape=:circle,legend=:topleft,markersize=1)
# savefig("Laplace_fxx_fyy.pdf")


text:
using Plots

function khayyam_triangle(lentingle)
    output = []
    for i in 1:lentingle
        push!(output, [])
        for j in 1:i
            if j == 1
                append!(output[i], 1)
            elseif j == i
                append!(output[i], 1)
            else j != 1 && j != i
                add = output[i - 1][j] + output[i - 1][j - 1]
                append!(output[i], add)
            end
        end
    end
    return output
end

order = 2^11
triangle = khayyam_triangle(order)
x = []
y = []
for i in 1:order
    for j in 1:i
        if triangle[i][j]%2 != 0
            push!(x, j - i/2)
            push!(y, -i)
        end
    end
end
scatter(x,y,markersize = 0.000001, legend = false, border=:none, dpi=300)
savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\A-few-fractals-in-Julia\\Fractals-Fig\\SKT.png")

max_stars_repo_name: shahmari/A-few-fractals-in-Julia

max_issues_repo_name: shahmari/A-few-fractals-in-Julia

max_forks_repo_name: shahmari/A-few-fractals-in-Julia


text:
include("GEMM.jl")
include("HEMM.jl")
include("HER2K.jl")
include("HERK.jl")
include("SYMM.jl")
include("SYR2K.jl")
include("SYRK.jl")
include("TRMM.jl")
include("TRSM.jl")


text:
using BenchmarkTools
using Distances

##
function nu_L2!(nu, xa, xb, l)
  nu .= 0.0
  for ka = axes(xa, 2)
    @simd for kd = axes(xa, 1)
      @inbounds nu[ka] = nu[ka] + 0.5*(((xa[kd, ka] - xb[kd])/l[kd])^2)
    end
  end
end

function kern_sqexp!(out, nu)
  out[:] = exp.(-nu)
end

function build_K!(K, xa, xb, l, kern!, nufunc!)
  # Build a kernel matrix using a function `kern` to construct columns/rows
  nu = Array{Float64}(undef, size(xa, 2))
  for kb = axes(xb, 2)
    @views nufunc!(nu, xa, xb[:, kb], l)
    @views kern!(K[:, kb], nu)
  end
end

function kern_sqexp_elem(nu)
  exp(-nu)
end

@inline function nu_L2_elem(xa, xb, l)
  # mapreduce(k -> ((xa[k] - xb[k])/l[k])^2, +, axes(xa, 1)) / 2
  out = 0.0
  @inbounds @simd for k in axes(xa, 1)
    out += ((xa[k] - xb[k])/l[k])^2
  end
  out = 0.5*out
end

@inline function nu_L2_index(xa, xb, ka, kb, l)
  out = 0.0
  @inbounds @simd for k in axes(xa, 1)
    out += ((xa[k, ka] - xb[k, kb])/l[k])^2
  end
  out = 0.5*out
end

function build_K_elem!(K, xa, xb, l, kern::Function, nufunc::Function)
  # Build a kernel matrix using a function `kern` to construct elements
  @inbounds for kb = axes(xb, 2)
    for ka = axes(xa, 2)
      #@views nu = sqeuclidean(xa[:, ka], xb[:, kb])
      #@views nu = nu_L2_index(xa, xb, ka, kb, l)
      @views nu = nufunc(xa[:, ka], xb[:, kb], l)
      K[ka, kb] = kern(nu)
      #@views nu = nu_L2_elem(xa[:, ka], xb[:, kb], l)
      #K[ka, kb] = kern_sqexp_elem(nu)
    end
  end
end

nd = 4
na = 4096
nb = na

l = ones(nd)
xa = rand(nd, na)
xb = xa
K = zeros((na, nb))
##
@time build_K_elem!(K, xa, xb, l, kern_sqexp_elem, nu_L2_elem)
#@code_native build_K_elem!(K, xa, xb, l, kern_sqexp_elem, nu_L2_elem)
##
@time build_K!(K, xa, xb, l, kern_sqexp!, nu_L2!)
##
# @btime @views K = [kern_sqexp_elem(nu_L2_index(xa, xb, ka, kb, l)) for ka in axes(xa, 2), kb in axes(xb, 2)]
# print()
##
# using KernelFunctions

# k₁ = SqExponentialKernel()
# ##
# @time kernelmatrix!(K₁,k₁,xa,obsdim=2)
# print()

max_stars_repo_name: krystophny/profit

max_issues_repo_name: krystophny/profit

max_forks_repo_name: krystophny/profit


text:
"""
`module Build` : this is a very poor man's version of `ase.build`. At the moment
only a subset of the `bulk` functionality is supported.
"""
module Build

import JuLIP
using ..Chemistry
using JuLIP: JVec, JMat, JVecF, JMatF, JVecsF, mat,
      Atoms, cell, cell_vecs, positions, momenta, masses, numbers, pbc,
      chemical_symbols, set_cell!, set_pbc!, update_data!,
      set_defm!, defm

import Base: union

export repeat, bulk, cluster, autocell!, append

const _unit_cells = Dict(     # (positions, cell matrix, factor of a)
   :fcc => ( [ [0 0 0], ],
             [0 1 1; 1 0 1; 1 1 0],  0.5),
   :bcc => ( [ [0.0,0.0,0.0], ],
             [-1 1 1; 1 -1 1; 1 1 -1], 0.5),
   :diamond => ( [ [0.0, 0.0, 0.0], [0.5, 0.5, 0.5] ],
                 [0 1 1; 1 0 1; 1 1 0], 0.5)
)

const _cubic_cells = Dict(   # (positions, factor of a)
   :fcc => ( [ [0 0 0], [0 1 1], [1 0 1], [1 1 0] ], 0.5 ),
   :bcc => ( [ [0 0 0], [1 1 1] ], 0.5 ),
   :diamond => ( [ [0 0 0], [1 1 1], [0 2 2], [1 3 3], [2 0 2],
                   [3 1 3], [2 2 0], [3 3 1] ], 0.25 )
)

_simple_structures = [:fcc, :bcc, :diamond]

function _simple_bulk(sym::Symbol, cubic::Bool)
   if cubic
      X, scal = _cubic_cells[symmetry(sym)]
      C = eye(3) / scal
   else
      X, C, scal = _unit_cells[symmetry(sym)]
   end
   a = lattice_constant(sym)
   return [ JVecF(x) * a * scal  for x in X ], JMatF(C * a * scal)
end

function _bulk_hcp(sym::Symbol)
   D = Chemistry.refstate(sym)
   a = D["a"]
   c = a * D["c/a"]
   return [JVecF(0.0, 0.0, 0.0), JVecF(0.0, a / sqrt(3), c / 2)],
          JMatF( [a, -a / 2,  0.0, 0.0, a * sqrt(3)/2, 0.0, 0.0, 0.0, c] )
end

_convert_pbc(pbc::NTuple{3, Bool}) = pbc
_convert_pbc(pbc::Bool) = (pbc, pbc, pbc)

function bulk(sym::Symbol; cubic = false, pbc = (true,true,true))
   symm = symmetry(sym)
   if symm in _simple_structures
      X, C = _simple_bulk(sym, cubic)
   elseif symm == :hcp
      X, C = _bulk_hcp(sym)  # cubic parameter is irrelevant for hcp
   end
   m = atomic_mass(sym)
   z = atomic_number(sym)
   nat = length(X)
   return Atoms( X, fill(zero(JVecF), nat), fill(m, nat), fill(z, nat), C,
                 _convert_pbc(pbc) )
end

"""
auxiliary function to convert a general cell to a cubic one; this is a bit of
a hack, so to not make a mess of things, this will only work in very restrictive
circumstances and otherwise throw an error. But it could be revisited.
"""
function _cubic_cell(atu::Atoms)
   @assert length(atu) == 1
   ru = JuLIP.rmin(atu)
   at = bulk(chemical_symbol(atu.Z[1]), cubic=true)
   r = JuLIP.rmin(at)
   return set_defm!(at, (ru/r) * defm(at); updatepositions=true)
end

"""
`cluster(args...; kwargs...) -> at::AbstractAtoms`

Produce a cluster of approximately radius R. The center
atom is always at index 1.

## Methods
```
cluster(atu::AbstractAtoms, R::Real)  # specify unit cell
cluster(sym::Symbol, R::Real)         # atu = bulk(sym; cubic=true)
```
Both methods assume that there is only a single species, and both require
the use of an orthorhombic unit cell (for now).

* `sym`: chemical symbol of the atomic species
* `R` : length-scale, meaning depends on shape of the cluster

## Keyword Arguments:
* `dims` : dimension into which the cluster is extended, typically
   `(1,2,3)` for 3D point defects and `(1,2)` for 2D dislocations, in the
   remaining dimension(s) the b.c. will be periodic.
* `shape` : shape of the cluster, at the moment only `:ball` is allowed

## TODO
 * lift the restriction of single species
 * allow other shapes
"""
function cluster{T}(atu::Atoms{T}, R::Real; dims = [1,2,3], shape = :ball, x0=nothing)
   sym = chemical_symbols(atu)[1]
   # check that the cell is orthorombic
   if !isdiag(cell(atu))
      if length(atu) > 1
         error("""`JuLIP.cluster` requires as argument either a cubic cell or a
                  one-atom cell.""")
      end
      atu = _cubic_cell(atu)
   end
   # # check that the first index is the centre
   # @assert norm(atu[1]) == 0.0
   # determine by how much to multiply in each direction
   Fu = cell(atu)'
   L = [ j ∈ dims ? 2 * (ceil(Int, R/Fu[j,j])+3) : 1    for j = 1:3]
   # multiply
   at = Atoms(atu) * L
   # find point closes to centre
   x̄ = mean( x[dims] for x in at.X)
   i0 = findmin( norm(x[dims] - x̄) for x in at.X )[2]
   if x0 == nothing
      x0 = at[i0]
   else
      x0 += at[i0]
   end
   # swap positions
   X = positions(at)
   X[1], X[i0] = X[i0], X[1]
   F = diagm([Fu[j,j]*L[j] for j = 1:3])
   # carve out a cluster with mini-buffer to account for round-off
   r = [ norm(x[dims] - x0[dims]) for x in X ]
   IR = find( r .<= R+sqrt(eps(T)) )
   # generate new positions
   Xcluster = X[IR]
   # generate the cluster
   at_cluster = Atoms(sym, Xcluster)
   set_cell!(at_cluster, F')
   # take open boundary in all directions specified in dims, periodic in the rest
   set_pbc!(at_cluster, [ !(j ∈ dims) for j = 1:3 ])
   # return the cluster and the index of the center atom
   return at_cluster
end

cluster(sym::Symbol, R::Real; kwargs...) =
      cluster( bulk(sym, cubic=true), R; kwargs... )

"""
```
repeat(at::Atoms, n::NTuple{3}) -> Atoms
repeat(at::Atoms, n::Integer) = repeat(at, (n,n,n))
```

Takes an `Atoms` configuration / cell and repeats it n_j times
into the j-th cell-vector direction. For example,
```
at = repeat(bulk("C"), (3,2,4))
```
creates 3 x 2 x 4 unit cells of carbon.

The same can be achieved by `*`:
```
at = bulk("C") * (3, 2, 4)
```
"""
function Base.repeat(at::Atoms, n::NTuple{3})
   C0 = cell(at)
   c1, c2, c3 = cell_vecs(at)
   X0, P0, M0, Z0 = positions(at), momenta(at), masses(at), numbers(at)
   nrep = n[1] * n[2] * n[3]
   nat0 = length(at)
   X = Base.repeat(X0, outer=nrep)
   P = Base.repeat(P0, outer=nrep)
   M = Base.repeat(M0, outer=nrep)
   Z = Base.repeat(Z0, outer=nrep)
   i = 0
   for a in CartesianRange( CartesianIndex(1,1,1), CartesianIndex(n...) )
      b = c1 * (a[1]-1) + c2 * (a[2]-1) + c3 * (a[3]-1)
      X[i+1:i+nat0] = [b+x for x in X0]
      i += nat0
   end
   return Atoms(X, P, M, Z, JMat(diagm([n...]) * cell(at)), pbc(at))
end

Base.repeat(at::Atoms, n::Integer) = repeat(at, (n,n,n))
Base.repeat(at::Atoms, n::AbstractArray) = repeat(at, (n...))

import Base.*
*(at::Atoms, n) = repeat(at, n)
*(n, at::Atoms) = repeat(at, n)

"""
`deleteat!(at::Atoms, n) -> at`:

returns the same atoms object `at`, but with the atom(s) specified by `n`
removed.
"""
function Base.deleteat!(at::Atoms, n)
   deleteat!(at.X, n)
   deleteat!(at.P, n)
   deleteat!(at.M, n)
   deleteat!(at.Z, n)
   update_data!(at, Inf)
   return at
end

"""
simple way to construct an atoms object from just positions
"""
Atoms(s::Symbol, X::Vector{JVecF}) =
      Atoms( X, fill(zero(JVecF), length(X)), fill(atomic_mass(s), length(X)),
             fill(atomic_number(s), length(X)), _autocell(X),
             (false, false, false) )

Atoms(s::Symbol, X::Matrix{Float64}) = Atoms(s, vecs(X))

function _autocell(X::Vector{JVec{T}}) where T
   ext = extrema(mat(X), 2)
   C = diagm([ e[2] - e[1] + 1.0  for e in ext ][:])
   return JMat{T}(C)
end

"""
`autocell!(at::Atoms) -> Atoms`

generates cell vectors that contain all atom positions + a small buffer,
sets the cell of `at` accordingly (in-place) and returns the same atoms objet
with the new cell.
"""
autocell!(at::Atoms) = set_cell!(at, _autocell(positions(at)))

union(at1::Atoms, at2::Atoms) =
   Atoms( X = union(at1.X, at2.X),
          P = union(at1.P, at2.P),
          M = union(at1.M, at2.M),
          Z = union(at1.Z, at2.Z),
          cell = cell(at1),
          pbc = pbc(at1) )

append(at::Atoms, X::JVecsF) =
   Atoms( X = union(at.X, X),
          P = union(at.P, zeros(JVecF, length(X))),
          M = union(at.M, zeros(length(X))),
          Z = union(at.Z, zeros(Int, length(X))),
          cell = cell(at),
          pbc = pbc(at) )

end


text:
using GeometryTypes, Colors

function loop{T}(point::T, radius, imax, i=1, points=Point3f0[], scales=Vec3f0[])
    i == imax && return points, scales
    push!(points, point); push!(scales, radius)
    for point in decompose(T, Sphere(point, radius), 4)
        loop(point, radius/2, imax, i+1, points, scales)
    end
    points, scales
end
const cmap2 = map(RGBA{Float32}, colormap("RdBu", 3))
function loop2{T}(point::T, normal, radius, imax, mesh=loadasset("cat.obj"), i=1, points=Point3f0[], scales=Vec3f0[], rotations=Vec3f0[], colors=RGBA{Float32}[])
    i > imax && return points, scales, rotations, colors
    push!(points, point);push!(scales, radius);push!(rotations, normal);push!(colors, cmap2[i])
    for (point, normal) in zip(decompose(T, mesh), decompose(Normal{3, Float32}, mesh))
        loop2(point, normal, radius/10f0, imax, mesh, i+1, points, scales, rotations, colors)
    end
    points, scales, rotations, colors
end
using GLVisualize
p, s, r, c = loop2(Point3f0(0), Normal{3,Float32}(0,0,1), 1f0, 2, loadasset("cat.obj"));

using GLPlot
glplot((loadasset("cat.obj"), p), scale = s, rotation = r, color = c)


text:
"""
    Edge

Directed edge from `parent` to `child`, both of type `Int` (for indices).
Field `length` can be missing.
"""
mutable struct Edge
    parent::Int
    child::Int
    length::Union{Missing,Float64}
end
# default edge length is missing:
Edge(parent, child) = Edge(parent, child, missing)

function Base.show(io::IO, edge::Edge)
    str  = "edge from $(edge.parent) to $(edge.child)"
    str *= (ismissing(edge.length) ? "" : ", length $(edge.length)")
    str *= "\n"
    print(io, str)
end

"""
    Tree{T}(edge, label, foo)

Tree (or directed graph) described by its list of `edge`s.
`label` is a dictionary whose values (node labels) are of type `T`,
such that labels could be strings, symbols, integers etc., but they
should be of the **same type** for all nodes in the tree.

Field `foo` is just to illustrate the difference between
`mutable struct` and `struct`.
"""
struct Tree{T}
    edge::Vector{Edge}
    label::Dict{Int64,T}
    foo::Int
end
# if no arguments at all: create empty tree with String labels, and foo=0
Tree() = Tree{String}(Vector{Edge}(), Dict{Int64,String}(), 0)
# given edges only: declare labels of type String, foo = 0
Tree(edge::Vector{Edge}) =  Tree{String}(edge, Dict{Int64,String}(), 0)
# given edges and labels: extract type of labels to get T and create the tree
Tree(edge::Vector{Edge}, label::AbstractDict, foo=0::Int) =
  Tree{eltype(values(label))}(edge, label, foo)

function Base.show(io::IO, tree::Tree{T}) where T
    str = "parent -> child:"
    for e in tree.edge
        str *= "\n$(e.parent) $(e.child)"
    end
    if !isempty(tree.label)
        str *= "\nlabels:"
        for (ind,lab) in tree.label
            str *= "\n$ind: $lab"
        end
    end
    print(io, str)
end


text:
"""
    NearestNeighborCluster(p, n)

Create a random cluster nearest-neighbour neutral landscape model with
values ranging 0-1. `p` sets the density of original clusters, and `n`
sets the neighborhood for clustering (see `?label` for neighborhood options)
"""
struct NearestNeighborCluster <: NeutralLandscapeMaker
    p::Float64
    n::Symbol
    function NearestNeighborCluster(p::Float64 = 0.5, n::Symbol = :rook)
        @assert p > 0
        @assert n ∈ (:rook, :queen, :diagonal)
        new(p,n)
    end
end

function _landscape!(mat, alg::NearestNeighborCluster)
    _landscape!(mat, NoGradient())
    classify!(mat, [alg.p, 1 - alg.p])
    replace!(mat, 2.0 => NaN)
    clusters, nClusters = label(mat, alg.n)
    coordinates = _coordinatematrix(clusters)
    sources = findall(!isnan, vec(clusters))
    tree = KDTree(coordinates[:,sources])
    clusters[:] .= clusters[sources[nn(tree, coordinates)[1]]]
    randvals = rand(nClusters)
    mat .= randvals[Int.(clusters)]
    mat
end


text:
tests = ["mathfuns", "functors", "rtypes"]

for t in tests
    fp = "$(t).jl"
    println("running $(fp) ...")
    include(fp)
end


text:
using TableReader
using Dates
using Test

@testset "readtsv" begin
    @testset "simple" begin
        # no data
        buffer = IOBuffer("""
        col1
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1]
        @test size(df) == (0, 1)

        buffer = IOBuffer("""
        col1\tcol2\tcol3
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test size(df) == (0, 3)

        # integers
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        1\t23\t456
        -10\t-99\t0
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1, -10]
        @test df[:col2] == [23, -99]
        @test df[:col3] == [456, 0]

        # floats
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        1.0\t1.1\t12.34
        -1.2\t0.0\t-9.
        .000\t.123\t100.000
        1e3\t1.E+123\t-8.2e-00
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1.0, -1.2, 0.000, 1e3]
        @test df[:col2] == [1.1, 0.0, 0.123, 1e123]
        @test df[:col3] == [12.34, -9.0, 100.000, -8.2]

        # bools
        buffer = IOBuffer("""
        col1\tcol2
        true\tfalse
        True\tFalse
        TrUe\tFaLsE
        TRUE\tFALSE
        t\tf
        T\tF
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2]
        @test df[:col1] == [true, true, true, true, true, true]
        @test df[:col2] == [false, false, false, false, false, false]

        # strings
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        a\tb\tc
        foo\tbar\tbaz
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == ["a", "foo"]
        @test df[:col2] == ["b", "bar"]
        @test df[:col3] == ["c", "baz"]

        # dates
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        2015-12-21\t2019-01-01\t1999-09-11
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [Date(2015, 12, 21)]
        @test df[:col2] == [Date(2019, 1, 1)]
        @test df[:col3] == [Date(1999, 9, 11)]

        # datetimes
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        2015-12-21T00:00:00\t2015-12-21T11:22:33\t2015-12-21T11:22:33.444
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [DateTime(2015, 12, 21, 0, 0, 0)]
        @test df[:col2] == [DateTime(2015, 12, 21, 11, 22, 33)]
        @test df[:col3] == [DateTime(2015, 12, 21, 11, 22, 33, 444)]
    end

    @testset "tricky float" begin
        # NaN
        buffer = IOBuffer("""
        col1\tcol2\tcol3\tcol4
        nan\tNan\tNaN\tNAN
        -nan\t+Nan\t-NaN\t+NAN
        """)
        df = readtsv(buffer)
        @test all(isnan.(df[:col1]))
        @test all(isnan.(df[:col2]))
        @test all(isnan.(df[:col3]))
        @test all(isnan.(df[:col4]))

        # Inf
        buffer = IOBuffer("""
        col1\tcol2\tcol3\tcol4
        inf\tInf\tInF\tINF
        -inf\t+Inf\t-InF\t+INF
        infinity\tInfinity\tInFiNiTy\tINFINITY
        -infinity\t+Infinity\t-InFiNiTy\t+INFINITY
        """)
        df = readtsv(buffer)
        @test all(isinf.(df[:col1]))
        @test all(isinf.(df[:col2]))
        @test all(isinf.(df[:col3]))
        @test all(isinf.(df[:col4]))
    end

    @testset "quotation" begin
        data = """
        "col1"\t"col2"\t"col3"
        "1"\t"23"\t"456"
        "-10"\t"-99"\t"0"
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, -10]
        @test df[:col2] == [23, -99]
        @test df[:col3] == [456, 0]

        # floats
        buffer = IOBuffer("""
        "col1"\t"col2"\t"col3"
        "1.0"\t"1.1"\t"12.34"
        "-1.2"\t"0.0"\t"-9."
        ".000"\t".123"\t"100.000"
        "1e3"\t"1.E+123"\t"-8.2e-00"
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1.0, -1.2, 0.000, 1e3]
        @test df[:col2] == [1.1, 0.0, 0.123, 1e123]
        @test df[:col3] == [12.34, -9.0, 100.000, -8.2]

        # strings
        buffer = IOBuffer("""
        col1\tcol2\tcol3
        "a"\t"b"\t"c"
        "foo"\t"bar"\t"baz"
        """)
        df = readtsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == ["a", "foo"]
        @test df[:col2] == ["b", "bar"]
        @test df[:col3] == ["c", "baz"]

        buffer = IOBuffer("""
        col1\tcol2\tcol3
        foo\tbar\tbaz
        foo\tbar\tbaz
        """)
        df = readtsv(buffer)
        @test df[:col1] == ["foo", "foo"]
        @test df[:col2] == ["bar", "bar"]
        @test df[:col3] == ["baz", "baz"]

        # quotation marks in a quoted string
        buffer = IOBuffer("""
        col1\tcol2
        " ""OK"" "\t"\""OK"\""
        """)
        df = readtsv(buffer)
        @test df[:col1] == [" \"OK\" "]
        @test df[:col2] == ["\"OK\""]
    end

    @testset "trimming" begin
        # trimming space
        data = """
        col1\tcol2\tcol3
        1   \t   2\t   3
           4\t   5\t   6
          7 \t 8  \t 9
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, 4, 7]
        @test df[:col2] == [2, 5, 8]
        @test df[:col3] == [3, 6, 9]

        df = readtsv(IOBuffer(data); trim = true)
        @test df[:col1] == [1, 4, 7]
        @test df[:col2] == [2, 5, 8]
        @test df[:col3] == [3, 6, 9]

        df = readtsv(IOBuffer(data); trim = false)
        @test df[:col1] == ["1   ", "   4", "  7 "]
        @test df[:col2] == ["   2", "   5", " 8  "]
        @test df[:col3] == ["   3", "   6", " 9  "]

        data = """
         col1  \t col2 \t col3
         foo   \t  b  \t baz
        """
        df = readtsv(IOBuffer(data); trim = true)
        @test df[:col1] == ["foo"]
        @test df[:col2] == ["b"]
        @test df[:col3] == ["baz"]

        df = readtsv(IOBuffer(data); trim = false)
        @test df[Symbol(" col1  ")] == [" foo   "]
        @test df[Symbol(" col2 ")] == ["  b  "]
        @test df[Symbol(" col3  ")] == [" baz"]
    end

    @testset "missing" begin
        data = """
        col1\tcol2\tcol3
        1\t2\t3
        \t5\t
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] isa Vector{Union{Int,Missing}}
        @test df[:col2] isa Vector{Int}
        @test df[:col3] isa Vector{Union{Int,Missing}}
        @test df[1,:col1] == 1
        @test ismissing(df[2,:col1])
        @test df[1,:col2] == 2
        @test df[2,:col2] == 5
        @test df[1,:col3] == 3
        @test ismissing(df[2,:col3])

        data = """
        col1\tcol2\tcol3
        1.0\t2.2\t-9.8
        \t10\t
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] isa Vector{Union{Float64,Missing}}
        @test df[:col2] isa Vector{Float64}
        @test df[:col3] isa Vector{Union{Float64,Missing}}
        @test df[1,:col1] == 1.0
        @test ismissing(df[2,:col1])
        @test df[1,:col2] == 2.2
        @test df[2,:col2] == 10.0
        @test df[1,:col3] == -9.8
        @test ismissing(df[2,:col3])

        data = """
        col1\tcol2\tcol3
        foo\t\tbar
        baz\tqux\t
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] isa Vector{String}
        @test df[:col2] isa Vector{Union{String,Missing}}
        @test df[:col3] isa Vector{Union{String,Missing}}
        @test df[1,:col1] == "foo"
        @test df[2,:col1] == "baz"
        @test ismissing(df[1,:col2])
        @test df[2,:col2] == "qux"
        @test df[1,:col3] == "bar"
        @test ismissing(df[2,:col3])
    end

    @testset "large integer" begin
        data = """
        col1\tcol2
        $(typemax(Int))\t$(typemin(Int))
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [typemax(Int)]
        @test df[:col2] == [typemin(Int)]

        # not supported
        data = """
        col1
        99999999999999999999999
        """
        @test_throws OverflowError readtsv(IOBuffer(data))
    end

    @testset "malformed data" begin
        # empty
        @test_throws TableReader.ReadError("found no column names in the header at line 1") readtsv(IOBuffer(""))

        # less columns than expected
        data = """
        col1\tcol2\tcol3
        foo\tbar
        """
        @test_throws TableReader.ReadError("unexpected number of columns at line 2") readtsv(IOBuffer(data))

        # more columns than expected
        data = """
        col1\tcol2\tcol3
        foo\tbar\tbaz\tqux\tquux
        """
        @test_throws TableReader.ReadError("unexpected number of columns at line 2") readtsv(IOBuffer(data))
    end

    @testset "UTF-8 strings" begin
        data = """
        col1\tcol2\tcol3
        甲\t乙\t丙
        👌\t😀😀😀\t🐸🐓
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == ["甲", "👌"]
        @test df[:col2] == ["乙", "😀😀😀"]
        @test df[:col3] == ["丙", "🐸🐓"]
    end

    @testset "CR+LF" begin
        data = """
        col1\tcol2\tcol3\r
        1\t2\t3\r
        4\t5\t6 \r
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1, 4]
        @test df[:col2] == [2, 5]
        @test df[:col3] == [3, 6]

        data = """
        col1\tcol2\tcol3\r
        1.0\t2.0\t3.0\r
        4.0\t5.0\t6.0 \r
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [1.0, 4.0]
        @test df[:col2] == [2.0, 5.0]
        @test df[:col3] == [3.0, 6.0]

        data = """
        col1\tcol2\tcol3\r
        foo\tbar\tbaz\r
        hoge\tfuga\tpiyo \r
        """
        df = readtsv(IOBuffer(data))
        @test df[:col1] == ["foo", "hoge"]
        @test df[:col2] == ["bar", "fuga"]
        @test df[:col3] == ["baz", "piyo"]
    end

    @testset "CR" begin
        data = """col1\tcol2\r123\t456\r"""
        df = readtsv(IOBuffer(data))
        @test df[:col1] == [123]
        @test df[:col2] == [456]
    end

    @testset "invalid argument" begin
        @test_throws ArgumentError readtsv(IOBuffer(""), chunksize = -1)
    end

    @testset "large data" begin
        buf = IOBuffer()
        m = 10000
        n = 2000
        println(buf, "name", '\t', join(("col$(j)" for j in 1:n), '\t'))
        for i in 1:m
            print(buf, "row$(i)")
            for j in 1:n
                print(buf, '\t', i)
            end
            println(buf)
        end
        data = take!(buf)

        # with chunking
        df = readtsv(IOBuffer(data))
        @test size(df) == (m, n + 1)
        @test df[:name] == ["row$(i)" for i in 1:m]
        @test df[:col1] == 1:m
        @test df[:col2] == 1:m

        # without chunking
        df = readtsv(IOBuffer(data), chunksize = 0)
        @test size(df) == (m, n + 1)
        @test df[:name] == ["row$(i)" for i in 1:m]
        @test df[:col1] == 1:m
        @test df[:col2] == 1:m
    end

    @testset "from file" begin
        df = readtsv(joinpath(@__DIR__, "test.tsv"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readtsv(joinpath(@__DIR__, "test.tsv.gz"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readtsv(joinpath(@__DIR__, "test.tsv.zst"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readtsv(joinpath(@__DIR__, "test.tsv.xz"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]
    end
end

@testset "readcsv" begin
    @testset "simple" begin
        # integers
        buffer = IOBuffer("""
        col1,col2,col3
        1,23,456
        -10,-99,0
        """)
        df = readcsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1, -10]
        @test df[:col2] == [23, -99]
        @test df[:col3] == [456, 0]

        # floats
        buffer = IOBuffer("""
        col1,col2,col3
        1.0,1.1,12.34
        -1.2,0.0,-9.
        .000,.123,100.000
        1e3,1.E+123,-8.2e-00
        """)
        df = readcsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == [1.0, -1.2, 0.000, 1e3]
        @test df[:col2] == [1.1, 0.0, 0.123, 1e123]
        @test df[:col3] == [12.34, -9.0, 100.000, -8.2]

        # strings
        buffer = IOBuffer("""
        col1,col2,col3
        a,b,c
        foo,bar,baz
        """)
        df = readcsv(buffer)
        @test eof(buffer)
        @test names(df) == [:col1, :col2, :col3]
        @test df[:col1] == ["a", "foo"]
        @test df[:col2] == ["b", "bar"]
        @test df[:col3] == ["c", "baz"]
    end

    @testset "header" begin
        data = """
        1,2,3
        4,5,6
        """
        df = readcsv(IOBuffer(data), header = ["col1", "col2", "col3"])
        @test df[:col1] == [1, 4]
        @test df[:col2] == [2, 5]
        @test df[:col3] == [3, 6]
    end

    @testset "quotation" begin
        data = """
        "col1","col2"
        "hi, there",","
        "1,2,3,4", ",,,"
        """
        df = readcsv(IOBuffer(data))
        @test df[:col1] == ["hi, there", "1,2,3,4"]
        @test df[:col2] == [",", ",,,"]
    end

    @testset "missing" begin
        data = """
        "col1","col2"
        "",""
        "1","2"
        """
        df = readcsv(IOBuffer(data))
        @test ismissing(df[1,:col1])
        @test df[2,:col1] == 1
        @test ismissing(df[1,:col2])
        @test df[2,:col2] == 2
    end

    @testset "multiline field" begin
        data = """
        col1,col2
        "oh,
        there
        there","
        multi
        line
        field
        "
        """
        df = readcsv(IOBuffer(data))
        @test df[:col1] == ["oh,\nthere\nthere"]
        @test df[:col2] == ["\nmulti\nline\nfield\n"]
    end

    @testset "EOF without newline" begin
        data = """
        col1,col2
        1,2"""
        df = readcsv(IOBuffer(data))
        @test df[:col1] == [1]
        @test df[:col2] == [2]
    end

    @testset "unnamed column" begin
        data = """
        col1,,col3,
        1,foo,3,bar
        """
        df = readcsv(IOBuffer(data))
        @test df[:col1] == [1]
        @test df[:UNNAMED_2] == ["foo"]
        @test df[:col3] == [3]
        @test df[:UNNAMED_4] == ["bar"]
    end

    @testset "implicit column" begin
        data = """
        col1,col2,col3
        foo,1,2,3
        """
        df = readcsv(IOBuffer(data))
        @test df[:UNNAMED_0] == ["foo"]
        @test df[:col1] == [1]
        @test df[:col2] == [2]
        @test df[:col3] == [3]
    end

    @testset "skip lines" begin
        data = """
        foobarbaz
        hogehogehoge
        col1,col2,col3
        1,2,3
        """
        df = readcsv(IOBuffer(data), skip = 2)
        @test df[:col1] == [1]
        @test df[:col2] == [2]
        @test df[:col3] == [3]
    end

    @testset "malformed data" begin
        # empty
        @test_throws TableReader.ReadError("found no column names in the header at line 1") readcsv(IOBuffer(""))
    end

    @testset "from file" begin
        df = readcsv(joinpath(@__DIR__, "test.csv"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readcsv(joinpath(@__DIR__, "test.csv.gz"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readcsv(joinpath(@__DIR__, "test.csv.zst"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]

        df = readcsv(joinpath(@__DIR__, "test.csv.xz"))
        @test df[:col1] == [1, 2]
        @test df[:col2] == [1.0, 2.0]
        @test df[:col3] == ["one", "two"]
    end

    @testset "from command" begin
        if Sys.which("echo") === nothing
            @info "skip tests: echo command is not found"
        else
            df = readcsv(`echo $("col1,col2\n1,2")`)
            @test df[:col1] == [1]
            @test df[:col2] == [2]
        end

        if Sys.which("cat") === nothing || Sys.which("gzip") === nothing
            @info "skip tests: cat/gzip commands are not found"
        else
            testfile = joinpath(@__DIR__, "test.csv")
            df = readcsv(pipeline(`cat $(testfile)`, `gzip`))
            @test df[:col1] == [1, 2]
            @test df[:col2] == [1.0, 2.0]
            @test df[:col3] == ["one", "two"]
        end
    end
end

@testset "readdlm" begin
    @testset "custom parser parameters" begin
        # delimtied by |
        data = """
        col1|col2|col3
        1|2|3
        """
        df = readdlm(IOBuffer(data), delim = '|')
        @test df[:col1] == [1]
        @test df[:col2] == [2]
        @test df[:col3] == [3]

        # quoted by `
        data = """
        col1,col2,col3
        `foo`,`bar`,`baz`
        """
        df = readdlm(IOBuffer(data), delim = ',', quot = '`')
        @test df[:col1] == ["foo"]
        @test df[:col2] == ["bar"]
        @test df[:col3] == ["baz"]

        # invalid chunksize
        @test_throws ArgumentError readdlm(IOBuffer(""), delim = ',', chunksize = -1)
        @test_throws ArgumentError readdlm(IOBuffer(""), delim = ',', chunksize = 2^63-1)
    end
end


text:
@testset "954.array-of-doubled-pairs.jl" begin
    @test can_reorder_doubled([3, 1, 3, 6]) == false
    @test can_reorder_doubled([2, 1, 2, 6]) == false
    @test can_reorder_doubled([4, -2, 2, -4]) == true
    @test can_reorder_doubled([1, 2, 3, 16, 8, 4]) == false
end


text:
# Use
#
#     DOCUMENTER_DEBUG=true julia --color=yes make.jl local [nonstrict] [fixdoctests]
#
# for local builds.

using Documenter
import Literate
using MGVI

SRC=joinpath(@__DIR__, "src")

GENERATED = joinpath(@__DIR__, "build")
GENERATED_SRC = joinpath(GENERATED, "src")
GENERATED_DOCS = joinpath(GENERATED, "docs")

EXECUTE_MD = true

mkdir(GENERATED)
cp(SRC, GENERATED_SRC)

function dir_replace(content)
    content = replace(content, "@__DIR__" => '"' * GENERATED_SRC * '"')
    return content
end

Literate.notebook(joinpath(GENERATED_SRC, "advanced_tutorial_lit.jl"),
                  GENERATED_SRC; name="advanced_tutorial", execute=false)
Literate.script(joinpath(GENERATED_SRC, "advanced_tutorial_lit.jl"),
                GENERATED_SRC; name="advanced_tutorial")
Literate.markdown(joinpath(GENERATED_SRC, "advanced_tutorial_lit.jl"),
                  GENERATED_SRC; name="advanced_tutorial", execute=EXECUTE_MD, preprocess=dir_replace)

Literate.notebook(joinpath(GENERATED_SRC, "tutorial_lit.jl"),
                  GENERATED_SRC; name="tutorial", execute=false)
Literate.script(joinpath(GENERATED_SRC, "tutorial_lit.jl"),
                GENERATED_SRC; name="tutorial")
Literate.markdown(joinpath(GENERATED_SRC, "tutorial_lit.jl"),
                  GENERATED_SRC; name="tutorial", execute=EXECUTE_MD)

makedocs(
    sitename = "MGVI",
    modules = [MGVI],
    root = GENERATED,
    source = "src",
    build = "docs",
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://bat.github.io/MGVI.jl/stable/"
    ),
    pages = [
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Advanced Tutorial" => "advanced_tutorial.md",
        "API" => "api.md",
        "LICENSE" => "LICENSE.md",
    ],
    doctest = ("fixdoctests" in ARGS) ? :fix : true,
    linkcheck = ("linkcheck" in ARGS),
    strict = !("nonstrict" in ARGS),
)

deploydocs(
    root = GENERATED,
    target = "docs",
    repo = "github.com/bat/MGVI.jl.git",
    forcepush = true,
    push_preview = true
)


text:
using ExactD, Test, LinearAlgebra

@testset "expectation values" begin
L = 5
N = 3
basis = make_LN_basis(L, N)
s = normalize(rand(length(basis)))
r = normalize(rand(length(basis)))

@testset "<s|n_p1*n_p2*...*n_pn|s>" begin
    p = 1
    @test expected_n(basis, s, p) ≈ s[1]^2+s[2]^2+s[3]^2+s[5]^2+s[6]^2+s[8]^2
    p = 4
    @test expected_n(basis, s, p) ≈ s[2]^2+s[3]^2+s[4]^2+s[8]^2+s[9]^2+s[10]^2

    p = [1, 3]
    @test expected_n(basis, s, p) ≈ s[1]^2+s[3]^2+s[6]^2
    p = [4, 2]
    @test expected_n(basis, s, p) ≈ s[2]^2+s[4]^2+s[9]^2

    p = [1, 3, 5]
    @test expected_n(basis, s, p) ≈ s[6]^2
    p = [4, 2, 3]
    @test expected_n(basis, s, p) ≈ s[4]^2

    p = [1, 2, 3, 5]
    @test expected_n(basis, s, p) ≈ 0. atol=1e-10
end

@testset "<r|n_p1*n_p2*...*n_pn|s>" begin
    p = 1
    @test expected_n(basis, r, s, p) ≈ s[1]*r[1]+s[2]*r[2]+s[3]*r[3]+s[5]*r[5]+s[6]*r[6]+s[8]*r[8]
    p = 4
    @test expected_n(basis, r, s, p) ≈ s[2]*r[2]+s[3]*r[3]+s[4]*r[4]+s[8]*r[8]+s[9]*r[9]+s[10]*r[10]

    p = [1, 3]
    @test expected_n(basis, r, s, p) ≈ s[1]*r[1]+s[3]*r[3]+s[6]*r[6]
    p = [4, 2]
    @test expected_n(basis, r, s, p) ≈ s[2]*r[2]+s[4]*r[4]+s[9]*r[9]

    p = [1, 3, 5]
    @test expected_n(basis, r, s, p) ≈ s[6]*r[6]
    p = [4, 2, 3]
    @test expected_n(basis, r, s, p) ≈ s[4]*r[4]

    p = [1, 2, 3, 5]
    @test expected_n(basis, r, s, p) ≈ 0. atol=1e-10
end

@testset "<s|b^+_p[1]*...*b^+_p[n]*b_q[1]*...*b_q[m]|s>" begin
    p = 1
    q = 1
    @test expected_pq(basis, s, p, q) ≈ s[1]^2+s[2]^2+s[3]^2+s[5]^2+s[6]^2+s[8]^2
    p = 4
    q = 4
    @test expected_pq(basis, s, p, q) ≈ s[2]^2+s[3]^2+s[4]^2+s[8]^2+s[9]^2+s[10]^2

    p = [1, 2]
    q = [1, 3]
    @test expected_pq(basis, s, p, q) ≈ s[2]*s[3] + s[5]*s[6]
    p = [4, 2]
    q = [1, 3]
    @test expected_pq(basis, s, p, q) ≈ s[6]*s[9]
end

@testset "<r|b^+_p[1]*...*b^+_p[n]*b_q[1]*...*b_q[m]|s>" begin
    p = 1
    q = 1
    @test expected_pq(basis, r, s, p, q) ≈ s[1]*r[1]+s[2]*r[2]+s[3]*r[3]+s[5]*r[5]+s[6]*r[6]+s[8]*r[8]
    p = 4
    q = 4
    @test expected_pq(basis, r, s, p, q) ≈ s[2]*r[2]+s[3]*r[3]+s[4]*r[4]+s[8]*r[8]+s[9]*r[9]+s[10]*r[10]

    p = [1, 2]
    q = [1, 3]
    @test expected_pq(basis, r, s, p, q) ≈ r[2]*s[3] + r[5]*s[6]
    p = [4, 2]
    q = [1, 3]
    @test expected_pq(basis, r, s, p, q) ≈ s[6]*r[9]
end
end


text:
module Citations

# include("bib.jl")

export @add_reference
export @bibliography
export cite

include("types.jl")

global_reference_list = Dict()

macro add_reference(references...)
    for ref in references
        r = eval(ref)
        global_reference_list[r.first] = r.second
    end
end

function parse_bib(bib_contents)
    # s = convert.(Ref(String), split(bib_contents, "\n"))
    # return s
    bib_contents
end

macro bibliography(bib_file)
    contents = readlines(bib_file)
    contents = join(contents, "\n")
    @show typeof(contents)
    citation_list = parse_bib(contents)
    @show typeof(citation_list)
    @show citation_list
    for (k,v) in citation_list
        global_reference_list[k] = v
    end
    return nothing
end

"""
    cite(ref)

Cite reference in `global_reference_list`
"""
cite(ref) = global_reference_list[ref]

# function __init__()
#     println("Hello!")
#     refs = []
#     @add_reference(refs...)
# end

end # module


text:
export AbstractPipelineStage, EmptyPipelineStage, CompositePipelineStage, apply_pipeline_stage

import Base: +

## Abstract Custom Pipeline Stage

"""
    AbstractPipelineStage

An abstract type for all custom pipelines stages

See also: [`apply_pipeline_stage`](@ref), [`EmptyPipelineStage`](@ref), [`CompositePipelineStage`](@ref)
"""
abstract type AbstractPipelineStage end

"""
    apply_pipeline_stage(stage, factornode, tag, stream)

Applies a given pipeline stage to the `stream` argument given `factornode` and `tag` of an edge.

See also: [`AbstractPipelineStage`](@ref), [`EmptyPipelineStage`](@ref), [`CompositePipelineStage`](@ref)
"""
function apply_pipeline_stage end

## Default empty pipeline

"""
    EmptyPipelineStage <: AbstractPipelineStage

Dummy empty pipeline stage that does not modify the original pipeline.

See also: [`AbstractPipelineStage`](@ref), [`apply_pipeline_stage`](@ref), [`CompositePipelineStage`](@ref)
"""
struct EmptyPipelineStage <: AbstractPipelineStage end

apply_pipeline_stage(::EmptyPipelineStage, factornode, tag, stream) = stream

## Composite pipeline

"""
    CompositePipelineStage{T} <: AbstractPipelineStage

Composite pipeline stage that consists of multiple inner pipeline stages

See also: [`AbstractPipelineStage`](@ref), [`apply_pipeline_stage`](@ref), [`EmptyPipelineStage`](@ref)
"""
struct CompositePipelineStage{T} <: AbstractPipelineStage
    stages :: T
end

apply_pipeline_stage(composite::CompositePipelineStage, factornode, tag, stream) = reduce((stream, stage) -> apply_pipeline_stage(stage, factornode, tag, stream), composite.stages, init = stream)

Base.:+(stage::AbstractPipelineStage) = stage

Base.:+(left::EmptyPipelineStage,     right::EmptyPipelineStage)     = EmptyPipelineStage()
Base.:+(left::EmptyPipelineStage,     right::AbstractPipelineStage)  = right
Base.:+(left::AbstractPipelineStage,  right::EmptyPipelineStage)     = left
Base.:+(left::AbstractPipelineStage,  right::AbstractPipelineStage)  = CompositePipelineStage((left, right))
Base.:+(left::AbstractPipelineStage,  right::CompositePipelineStage) = CompositePipelineStage((left, right.stages...))
Base.:+(left::CompositePipelineStage, right::AbstractPipelineStage)  = CompositePipelineStage((left.stages..., right))
Base.:+(left::CompositePipelineStage, right::CompositePipelineStage) = CompositePipelineStage((left.stages..., right.stages...))


text:
## Environment and packages
import Pkg
Pkg.activate(dirname(Base.current_project())) # activate project

# using Logging: global_logger
# using TerminalLoggers: TerminalLogger
# global_logger(TerminalLogger())

using Distributed
using ProgressMeter
const processes = 10

if nprocs() < processes
    addprocs(processes - nprocs(); exeflags="--project")
end

println("Number of cores: ", nprocs())
println("Number of workers: ", nworkers())

@everywhere begin
    import Pkg
    Pkg.activate(dirname(Base.current_project()))
    Pkg.precompile()
end

@everywhere begin
using Statistics
using LinearAlgebra
using Random
using HDF5
using JLD2
using OrdinaryDiffEq
using DiffEqFlux
using Flux
using Tullio
using RecursiveArrayTools
using Infiltrator
using Plots
using ProgressMeter

include("helpers/parameters.jl")

#######################################################################################################
#############################             Glacier Setup            ####################################
#######################################################################################################

# Load the HDF5 file with Harry's simulated data
cd(@__DIR__)
root_dir = cd(pwd, "../..")
argentiere_f = h5open(joinpath(root_dir, "data/Argentiere_2003-2100_aflow2e-16_50mres_rcp2.6.h5"), "r")

struct Glacier
    bed::Array{Float32}    # bedrock height
    thick::Array{Float32}  # ice thickness
    vel::Array{Float32}    # surface velocities
    MB::Array{Float32}     # surface mass balance
    lat::Float32
    lon::Float32
end

# Fill the Glacier structure with the retrieved data
argentiere = Glacier(HDF5.read(argentiere_f["bed"])[begin:end-2,:],
                     HDF5.read(argentiere_f["thick_hist"])[begin:end-2,:,2:end],
                     HDF5.read(argentiere_f["vel_hist"])[begin:end-2,:,2:end],
                     HDF5.read(argentiere_f["s_apply_hist"])[begin:end-2,:,2:end],
                     0, 0)

# Domain size
const nx = size(argentiere.bed)[1]
const ny = size(argentiere.bed)[2]

(@isdefined B) || (const B  = copy(argentiere.bed))
(@isdefined H₀) || (const H₀ = copy(argentiere.thick[:,:,1]))

# Spatial and temporal differentials
const Δx, Δy = 50, 50 #m (Δx = Δy)

end # @everywhere

# Include all functions
include("helpers/iceflow.jl")

#######################################################################################################
#############################           Reference dataset          ####################################
#######################################################################################################

(@isdefined temp_series) || (const temp_series, norm_temp_series = fake_temp_series(t₁))

if create_ref_dataset
    println("Generating reference dataset for training...")

    # Compute reference dataset in parallel
    @everywhere solver = Ralston()
    H_refs, V̄x_refs, V̄y_refs = generate_ref_dataset(temp_series, H₀)

    ## Save reference plots
    # for (temps, H) in zip(temp_series, H_refs)

    #     ###  Glacier ice thickness difference  ###
    #     lim = maximum( abs.(H[end] .- H₀) )
    #     hm = heatmap(H[end] .- H₀, c = cgrad(:balance,rev=true), aspect_ratio=:equal,
    #         clim = (-lim, lim),
    #         title="Variation in ice thickness")

    #     tempn = floor(mean(temps))
    #     println("Saving reference_$tempn.png...")
    #     savefig(hm,joinpath(root_dir,"plots/references","reference_$tempn.png"))
    # end

    println("Saving reference data")
    jldsave(joinpath(root_dir, "data/PDE_refs.jld2"); H_refs, V̄x_refs, V̄y_refs)

end

# Load stored PDE reference datasets
PDE_refs = load(joinpath(root_dir, "data/PDE_refs.jld2"))

#######################################################################################################
#############################             Train UDE            ########################################
#######################################################################################################

UA = FastChain(
        FastDense(1,3, x->softplus.(x)),
        FastDense(3,10, x->softplus.(x)),
        FastDense(10,3, x->softplus.(x)),
        FastDense(3,1, sigmoid_A)
    )

θ = initial_params(UA)
current_epoch = 1
batch_size = length(temp_series)

cd(@__DIR__)
const root_plots = cd(pwd, "../../plots")
# Train iceflow UDE in parallel
# First train with ADAM to move the parameters into a favourable space
@everywhere solver = ROCK4()
train_settings = (ADAM(0.05), 20) # optimizer, epochs
iceflow_trained = @time train_iceflow_UDE(H₀, UA, θ, train_settings, PDE_refs, temp_series)
θ_trained = iceflow_trained.minimizer

# Continue training with a smaller learning rate
# train_settings = (ADAM(0.001), 20) # optimizer, epochs
# iceflow_trained = @time train_iceflow_UDE(H₀, UA, θ, train_settings, H_refs, temp_series)
# θ_trained = iceflow_trained.minimizer

# Continue training with BFGS
train_settings = (BFGS(initial_stepnorm=0.02f0), 20) # optimizer, epochs
iceflow_trained = @time train_iceflow_UDE(H₀, UA, θ_trained, train_settings, PDE_refs, temp_series)
θ_trained = iceflow_trained.minimizer

# Save trained NN weights
save(joinpath(root_dir, "data/trained_weights.jld"), "θ_trained", θ_trained)

# Plot the final trained model
data_range = -20.0:0.0
pred_A = predict_A̅(UA, θ_trained, collect(data_range)')
pred_A = [pred_A...] # flatten
true_A = A_fake(data_range)

scatter(true_A, label="True A")
train_final = plot!(pred_A, label="Predicted A")
savefig(train_final,joinpath(root_plots,"training","final_model.png"))


text:
using GMT

# -------------------------------------------------------------------------
function blockmean()
	gmt("blockmean -R0/150/0/100 -I10 -V > V:/ave1.txt", rand(Float64,100,3)*150)
	B = gmt("blockmean -R0/150/0/100 -I10", rand(Float64,100,3)*100);
	gmt("write -Td V:/ave2.txt", B)
	return B
end

# -------------------------------------------------------------------------
function surface()
	gmt("surface -R0/150/0/100 -I1 -GV:/lixo.grd -V", rand(Float64,100,3)*150)
	G = gmt("surface -R0/150/0/100 -I1", rand(Float64,100,3)*150)
	return G
end

# -------------------------------------------------------------------------
function pscoast()
	gmt("pscoast -R-10/0/35/45 -B1 -W1 -Gbrown -JM14c -P -V > V:/lixo.ps")
end

# -------------------------------------------------------------------------
function gmtinfo()
	r = gmt("gmtinfo -C -V",ones(Float32,9,3)*5)
	info = gmt("blockmean -R0/5/0/5 -I1 -V", rand(Float32,2048,3)*5)
end

# -------------------------------------------------------------------------
function gmtread()
	G = gmt("gmtread -Tg V:/lixo.grd")
	return G
end

# -------------------------------------------------------------------------
function gmtwrite()
	G = gmt("gmtread -Tg V:/lixo.grd")
	gmt("write -Tg V:/crap.grd", G)
	return G
end

# -------------------------------------------------------------------------
function grdimage()
	G = gmt("surface -R0/150/0/100 -I1", rand(Float64,100,3)*150)
	gmt("grdimage -JX8c -Ba -P -Cblue,red > V:/crap_img.ps", G)
end

# -------------------------------------------------------------------------
function gridbox()
	grd = ones(Float32,9,9)*5
	hdr = [0.0,8,0,8,0,1,0,1,1]
	grd_box = array_container(9, 9, 1, pointer(grd), pointer(hdr))
	#gmt("grdinfo -L0 -V", grd_box)
	gmt("write -Tg V:/lixo.grd -Vl", grd_box)
end

# -------------------------------------------------------------------------
function grd2xyz_box()
	grd = reshape(Float32(1):12,3,4)
	hdr = [1.0,4,1,3,1,12,0,1,1]
	grd_box = array_container(4, 3, 1, pointer(grd), pointer(hdr))
	xyz = gmt("grd2xyz -V", grd_box)
	return xyz
end

# -------------------------------------------------------------------------
function grd2xyz()
	G = gmt("surface -R0/150/0/100 -I1", rand(Float64,100,3)*150)
	xyz = gmt("grd2xyz -V", G)
	return xyz
end

# -------------------------------------------------------------------------
function filter1d()
	t = zeros(100,2)
	t[:,1] = 1:100
	t[:,2] = rand(100,1) * 100
	gmt("filter1d -Fg10 -E > V:/filt1.txt", t)
	F = gmt("filter1d -Fg10 -E", t)
	gmt("write -Td V:/filt2.txt", F)
	return F
end

# -------------------------------------------------------------------------
function gmtsimplify()
	println("Test gmtsimplify")
	t = gmt("simplify -T0.2", rand(50,2))
end

# -------------------------------------------------------------------------
function makecpt()
	C = gmt("makecpt -Crainbow -T5/35/5")
	gmt("psscale -D4/4/4/0.5 -P -Baf > V:/t.ps", C)
end

# -------------------------------------------------------------------------
function grdinfo()
	G = gmt("surface -R0/150/0/150 -I1", rand(100,3) * 100)
	T = gmt("grdinfo -V", G)
end

# -------------------------------------------------------------------------
function pstext()
	lines = Any["5 6 Some label", "6 7 Another label"]
	gmt("pstext -R0/10/0/10 -JM6i -Bafg -F+f18p -P > V:/text.ps", lines)
end

# -------------------------------------------------------------------------
function jlogo(L)
	# Create the Julia "Terminator" 3 colored circles triangle
	# L is the length of the equilateral triangle
	W = 2 * L 					# Region width
	H = L * sind(60) 			# Triangle height
	s_size = 0.8 * L 			# Circle diameter
	l_thick = s_size * 0.06 	# Line thickness

	com = @sprintf("psxy -Sc%fc -G191/101/95  -Jx1 -R0/%f/0/%f -W%fc,171/43/33  -P -K >  V:/Jlogo.ps", s_size, W, W, l_thick)
	gmt(com, [L/2 L/2])
	com = @sprintf("psxy -Sc%fc -G158/122/190 -Jx1 -R0/%f/0/%f -W%fc,130/83/171 -O -K >> V:/Jlogo.ps", s_size, W, W, l_thick)
	gmt(com, [L+L/2 L/2])
	com = @sprintf("psxy -Sc%fc -G128/171/93  -Jx1 -R0/%f/0/%f -W%fc,81/143/24  -O    >> V:/Jlogo.ps", s_size, W, W, l_thick)
	gmt(com, [L L/2+H])
end

jlogo() = jlogo(5)
# If not provided, make a logo where the triangle has a side of 5 cm.


text:
"""
    apply_drawings!(g, drawings, figid, axidx)

Internal function to apply a vector of `Drawing` objects contained in an
`Axes` container in a GLE context. The `axidx` and `figid` help keep track
of where individual drawings belong to which is useful when writing auxiliary
files containing the drawing data.
"""
function apply_drawings!(
            g::GS,
            drawings::Vector{<:Drawing},
            figid::String,
            axidx::Int
        )::Nothing

    # the element counter helps indicate which column in auxilliary files
    # to look into. A drawing application increments it depending on
    # how many columns it looks at.
    el_cntr = 1
    for drawing ∈ drawings
        el_cntr = apply_drawing!(g, drawing, el_cntr, figid, axidx)
    end
    return
end

"""
    auxpath(h, fid, axid)

Internal function to generate a path to an auxiliary file storing drawing
data for the current axes of the current figure.
The argument `h` is the hash of what has to be written (so that if the exact
same command is issued, the file is not re-written). The `axidx` and `figid`
are there to make the auxiliary file unique and easy to delete after use.
"""
auxpath(h::UInt, fid::String, axid::Int) =
    GP_ENV["TMP_PATH"] / "$(fid)_$(axid)_$h.csv"

"""
    auxdata(o, fid, axid)

Write the auxilliary data to a Drawing to file taking into account the figure
and axes id to avoid replacing an-already existing file.
"""
function auxdata(
            o::Drawing,
            fid::String,
            axid::Int;
            #
            nomiss = false
        )::String

    # write data to a temporary CSV file
    faux = auxpath(hash(o.data), fid, axid)
    # don't rewrite if it's the exact same zipper
    isfile(faux) || csv_writer(faux, o.data, !nomiss && o.hasmissing)
    return faux
end

##########################################################################
##########################################################################

"""
    apply_drawing!(g, scatter, ...)

Apply Scatter2D plot. Missing values are allowed.
"""
function apply_drawing!(
            g::GS,
            scatter::Scatter2D,
            el_cntr::Int,
            figid::String,
            axidx::Int
        )::Int

    faux = auxdata(scatter, figid, axidx)
    #
    # GLE syntax is:
    #
    #   data datafile.dat d1=c1,c2
    #   d1 line blue ...
    #
    # we do that for each object in the scatter.
    #
    for k ∈ 1:scatter.nobj

        # (1) indicate what data to read
        "\n\tdata \"$faux\" d$(el_cntr)=c1,c$(k+1)" |> g

        # if no color has been specified, assign one according to the palette
        if !isdef(scatter.linestyles[k].color)
            scatter.linestyles[k].color = palette(el_cntr)
        end

        # (2) Line and Marker description
        # (2.A) - there is a line
        if scatter.linestyles[k].lstyle != -1
            # Line
            "\n\td$el_cntr line" |> g
            apply_linestyle!(g, scatter.linestyles[k])

            # Marker
            mcol   = false
            mscale = 1.0
            # > if a marker color is specified and is different than the
            # > current line color, we need to add a dedicated GLE routine
            # > that can handle this on the relevant figure.
            if isdef(scatter.markerstyles[k].color) &&
                    (scatter.markerstyles[k].color != scatter.linestyles[k].color)
                mcol = true
                f    = Figure(figid; _noreset=true)
                add_sub_marker!(f, scatter.markerstyles[k])
                mscale = f.textstyle.hei
            end
            apply_markerstyle!(g, scatter.markerstyles[k]; mcol, mscale)

        # (2.B) - only markers
        else
            # Scatter plot; if there's no specified marker color,
            # take the default line color
            if !isdef(scatter.markerstyles[k].color)
                scatter.markerstyles[k].color = scatter.linestyles[k].color
            end
            "\n\td$el_cntr" |> g
            apply_markerstyle!(g, scatter.markerstyles[k])
        end

        el_cntr += 1
    end
    return el_cntr
end

"""
    apply_drawing!(g, fill, ...)

Apply Fill2D plot. Missing values are NOT allowed.
"""
function apply_drawing!(
            g::GS,
            fill::Fill2D,
            el_cntr::Int,
            figid::String,
            axidx::Int
        )::Int

    faux = auxdata(fill, figid, axidx; nomiss=true)
    #
    # GLE syntax is:
    #
    #   data datafile.dat d1=c1,c2 d2=c1,c3
    #   fill d1,d2 color rgb(1,1,1) xmin 0 xmax 1
    #
    "\n\tdata \"$faux\" d$(el_cntr)=c1,c2 d$(el_cntr+1)=c1,c3" |> g

    "\n\tfill d$(el_cntr),d$(el_cntr+1)"    |> g
    "color $(col2str(fill.fillstyle.fill))" |> g
    add(g, fill, :xmin, :xmax)

    el_cntr += 2
    return el_cntr
end

"""
    apply_drawing!(g, hist, ...)

Apply Hist2D plot. Missing values are allowed.
"""
function apply_drawing!(
            g::GS,
            hist::Hist2D,
            el_cntr::Int,
            figid::String,
            axidx::Int
        )::Int

    faux = auxdata(hist, figid, axidx)
    #
    # GLE syntax is:
    #
    #   data datafile.dat d1
    #   let d2 = hist d1 from xmin to xmax bins nbins
    #   bar d2 width width$ fill col$ color col2$ pattern pat$ horiz
    #

    # if no color has been specified, assign one according to the PALETTE
    if !isdef(hist.barstyle.color)
        if hist.barstyle.fill == c"white"
            cc = mod(el_cntr, GP_ENV["SZ_PALETTE"])
            (cc == 0) && (cc = GP_ENV["SZ_PALETTE"])
            hist.barstyle.color = GP_ENV["PALETTE"][cc]
        else
            hist.barstyle.color = c"white" # looks nicer than black
        end
    end

    # (1) indicate what data to read
    "\n\tdata \"$faux\" d$(el_cntr)" |> g

    # (2) hist description  | let d(k+1) = hist dk from min to max
    minx, maxx = hist.range
    "\n\tlet d$(el_cntr+1) = hist d$(el_cntr)" |> g
    "from $minx to $maxx" |> g
    el_cntr += 1

    # number of bins
    nobs = hist.nobs
    nbauto = (nobs == 0 ? 1 : ceil(Integer, log2(nobs))+1) # sturges, see StatsBase
    bins = isdef(hist.bins) ? hist.bins : nbauto
    add(g, hist, :bins)

    # (3) compute appropriate scaling
    width   = (maxx - minx) / bins
    scaling = 1.0
    hist.scaling == "probability" && (scaling /= nobs)
    hist.scaling == "pdf"         && (scaling /= (nobs * width))
    "\n\tlet d$(el_cntr) = d$(el_cntr)*$scaling" |> g

    # (4) apply histogram
    "\n\tbar d$(el_cntr) width $width" |> g

    # apply styling
    apply_barstyle!(g, hist.barstyle)
    add(g, hist, :horiz)

    return el_cntr + 1
end

"""
    apply_drawing!(g, bar, ...)

Apply Bar2D plot. Missing values are allowed.
"""
function apply_drawing!(
            g::GS,
            bar::Bar2D,
            el_cntr::Int,
            figid::String,
            axidx::Int
        )::Int

    faux = auxdata(bar, figid, axidx)
    #
    # general GLE syntax is:
    # (1) data datafile.dat d1 d2 ...
    # (2) bar d1,d2,d3 fill color_ color color_
    # (STACKED, NO GROUP)
    # bar d1 ...
    # bar d2 from d1 ...
    #
    # XXX (STACKED, GROUP) ---> see first whether this has a use case...
    # could solve it by specifying what gets stacked stack=[(1, 2), ...] but
    # clunky (and tbf not great visualisation anyway...)
    # bar d1,d2 ...
    # bar d3,d4 from d1,d2 ...
    #

    # if no color has been specified, assign with palette
    for c ∈ eachindex(bar.barstyles)
        if !isdef(bar.barstyles[c].color)
            if bar.barstyles[c].fill == c"white"
                bar.barstyles[c].color = palette(el_cntr + c - 1)
            else
                bar.barstyles[c].color = c"white"
            end
        end
    end

    nbars = bar.nobj

    # (1) indicate what data to read "data file d1 d2 d3..."
    "\n\tdata \"$faux\""      |> g
    dlist(el_cntr .+ 1:nbars) |> g

    # (2) non stacked (or single barset)
    if nbars==1 || !bar.stacked
        # bar d1,d2,d3
        "\n\tbar"                          |> g
        dlist(el_cntr .+ 1:nbars, sep=",") |> g
        # styling
        isdef(bar.bwidth) && "width $(bar.bwidth)" |> g
        apply_barstyles_nostack!(g, bar.barstyles)
        add(g, bar, :horiz)

    # (2) stacked
    else
        # first base bar
        "\n\tbar d$(el_cntr)"                      |> g
        isdef(bar.bwidth) && "width $(bar.bwidth)" |> g
        apply_barstyle!(g, bar.barstyles[1])
        add(g, bar, :horiz)
        # bars stacked on top
        for i ∈ 2:nbars
            "\n\tbar d$(el_cntr+i-1) from d$(el_cntr+i-2)" |> g
            apply_barstyle!(g, bar.barstyles[i])
        end
    end

    return el_cntr + nbars
end

"""
    apply_drawing!(g, boxplot, ...)

Apply Boxplot.
"""
function apply_drawing!(
            g::GS,
            bp::Boxplot,
            el_cntr::Int,
            figid::String,
            axidx::Int
        )::Int

    # 1. add boxplot subroutines (vertical and horizontal) if not there already
    subname = ifelse(bp.horiz, "bp_horiz", "bp_vert")
    f       = Figure(figid; _noreset=true)
    if subname ∉ keys(f.subroutines)
        f.subroutines[subname] = GLE_DRAW_SUB[subname]
    end

    # draw the boxplots one by one
    for k ∈ 1:bp.nobj
        # 1. retrieve the statistics and map to spaced string
        stats_str = vec2str(bp.stats[k, :])

        # 2. call the subroutine and apply the style
        s = bp.boxstyles[k]
        "\n\tdraw $subname $k $stats_str" |> g
        apply_boxplotstyle!(g, s, f)

        el_cntr += 1
    end
    return el_cntr
end

"""
    apply_drawing!(g, hm, ...)

Apply Heatmap.
"""
function apply_drawing!(
            g::GS,
            hm::Heatmap,
            el_cntr::Int,
            figid::String,
            axidx::Int
        )::Int

    # 1. apply the subroutine with a hash depending on figid and axidx
    hashid = hash((figid, axidx))
    add_sub_heatmap!(Figure(figid; _noreset=true), hm, hashid)

    # 2. write the zfile
    faux = auxpath(hash(hm.data), figid, axidx)
    isfile(faux) || csv_writer(faux, hm.transpose ? hm.data' : hm.data, false)

    # infer bw/bh based on ncols/nrows
    nrows, ncols = size(hm.data)
    bw, bh       = 1.0 ./ (ncols, nrows)
    hm.transpose && (tmp = bw; bw = bh; bh = bw)

    # number of columns of the final heatmap
    nct = ifelse(hm.transpose, nrows, ncols)

    # 3. load data
    vs = join(("d$j=c0,c$j" for j ∈ 1:nct), " ")
    "\n\tdata \"$faux\" $vs" |> g

    # 4. go over all the columns and draw the boxes (scales linearly with # cols)
    for j ∈ 1:nct
        "\n\tdraw hm_$hashid $j d$j $bw $bh" |> g
    end
    return el_cntr + 1
end


text:
# Wrapper for level-2 BLIS routines.
#

# Level2 API of common forms are put together.
#
macro blis_group_level2_her2(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          uploa,  BliUpLo,
                          conjx,  BliConj,
                          conjy,  BliConj,
                          m,      BliDim,
                          α,      Ptr{xType},
                          x,      Ptr{xType}, incx, BliInc,
                          y,      Ptr{xType}, incy, BliInc,
                          a,      Ptr{xType}, rsa,  BliInc, csa, BliInc)
    end
end

@blis_group_level2_her2 her2
@blis_group_level2_her2 syr2

macro blis_group_level2_hemv(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          uploa,  BliUpLo,
                          conja,  BliConj,
                          conjx,  BliConj,
                          m,      BliDim,
                          α,      Ptr{xType},
                          a,      Ptr{xType}, rsa,  BliInc, csa, BliInc,
                          x,      Ptr{xType}, incx, BliInc,
                          β,      Ptr{xType},
                          y,      Ptr{xType}, incy, BliInc)
    end
end

@blis_group_level2_hemv hemv
@blis_group_level2_hemv symv

macro blis_group_level2_tr(funcname)
    return quote
        @blis_ccall_group($funcname,
                          Cvoid,
                          uploa,  BliUpLo,
                          transa, BliTrans,
                          diaga,  BliDiag,
                          m,      BliDim,
                          α,      Ptr{xType},
                          a,      Ptr{xType}, rsa,  BliInc, csa, BliInc,
                          x,      Ptr{xType}, incx, BliInc)
    end
end

@blis_group_level2_tr trmv
@blis_group_level2_tr trsv

# Forms with only 1 member are defined directly.
#
@blis_ccall_group(gemv,
                  Cvoid,
                  transa, BliTrans,
                  conjx,  BliConj,
                  m,      BliDim,
                  n,      BliDim,
                  α,      Ptr{xType},
                  a,      Ptr{xType}, rsa,  BliInc, csa, BliInc,
                  x,      Ptr{xType}, incx, BliInc,
                  β,      Ptr{xType},
                  y,      Ptr{xType}, incy, BliInc)
@blis_ccall_group(ger,
                  Cvoid,
                  conjx,  BliConj,
                  conjy,  BliConj,
                  m,      BliDim,
                  n,      BliDim,
                  α,      Ptr{xType},
                  x,      Ptr{xType}, incx, BliInc,
                  y,      Ptr{xType}, incy, BliInc,
                  a,      Ptr{xType}, rsa,  BliInc, csa, BliInc)
@blis_ccall_group(her,
                  Cvoid,
                  uploa,  BliUpLo,
                  conjx,  BliConj,
                  m,      BliDim,
                  α,      Ptr{rType},
                  x,      Ptr{xType}, incx, BliInc,
                  a,      Ptr{xType}, rsa,  BliInc, csa, BliInc)
@blis_ccall_group(syr,
                  Cvoid,
                  uploa,  BliUpLo,
                  conjx,  BliConj,
                  m,      BliDim,
                  α,      Ptr{xType},
                  x,      Ptr{xType}, incx, BliInc,
                  a,      Ptr{xType}, rsa,  BliInc, csa, BliInc)


text:
module bib

export Parser

include("parser.jl")

end # module


text:
for i = 1:10
    if i % 3 != 0
        continue
    end
    print(i)
    print(" ")
end
# expected: 3 6 9

print("\n")
for j in 1:1000
    print(j)
    print(" ")
    if j == 5
        break
    end
end
# expected: 1 2 3 4 5

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
using LightGraphs
function bench_maxflow(n)
    p = 8.0/n
    A = sprand(n,n,p)
    g = DiGraph(A)
    cap = round(A*100)
    @time maximum_flow(g, 1, n, cap)
end

for n in 3:14
    bench_maxflow(2^n)
end


text:
module Neat

export sigmoid
export Connection
export Node, NodeType, INPUT
export Network

include("activation.jl")
include("connection.jl")
include("node.jl")
include("network.jl")

end # module


text:
# Halfway Bounce Back rules
struct BounceBack{D <: Direction, Ints} <: BoundaryCondition
    direction::D
    xs::Ints
    ys::Ints
end

function apply!(bc::BounceBack{<:North}, q::Quadrature, f_new, f_old)
    nx, ny, nf = size(f_new)
    for f_idx in 1:size(f_new, 3)
        # Bounce back the unkown populations
        opposite_f_idx = opposite(q, f_idx)

        for y_idx in bc.ys
            if y_idx + q.abscissae[2, opposite_f_idx] <= ny
                continue
            end
            for x_idx in bc.xs
                f_new[x_idx, y_idx, f_idx] = f_old[x_idx, y_idx, opposite_f_idx]
            end
        end
    end
end

function apply!(bc::BounceBack{<:South}, q::Quadrature, f_new, f_old)
    nx, ny, nf = size(f_new)
    for f_idx in 1:size(f_new, 3)
        opposite_f_idx = opposite(q, f_idx)

        for y_idx in bc.ys
            if y_idx + q.abscissae[2, opposite_f_idx] > 0
                continue
            end
            for x_idx in bc.xs
                f_new[x_idx, y_idx, f_idx] = f_old[x_idx, y_idx, opposite_f_idx]
            end
        end
    end
end

function apply!(bc::BounceBack{<:East}, q::Quadrature, f_new, f_old)
    nx, ny, nf = size(f_new)
    for f_idx in 1:size(f_new, 3)
        # Bounce back the unkown populations
        opposite_f_idx = opposite(q, f_idx)

        for x_idx in bc.xs
            if x_idx + q.abscissae[1, opposite_f_idx] <= nx
                continue
            end
            for y_idx in bc.ys
                f_new[x_idx, y_idx, f_idx] = f_old[x_idx, y_idx, opposite_f_idx]
            end
        end
    end
end
function apply!(bc::BounceBack{<:West}, q::Quadrature, f_new, f_old)
    nx, ny, nf = size(f_new)
    for f_idx in 1:size(f_new, 3)
        # Bounce back the unkown populations
        opposite_f_idx = opposite(q, f_idx)

        for x_idx in bc.xs
            if x_idx + q.abscissae[1, opposite_f_idx] > 0
                continue
            end
            for y_idx in bc.ys
                f_new[x_idx, y_idx, f_idx] = f_old[x_idx, y_idx, opposite_f_idx]
            end
        end
    end
end


text:
@inline function ExaTronKernels.dcopy(n::Int,dx::CuDeviceArray{Float64,1},incx::Int,
                       dy::CuDeviceArray{Float64,1},incy::Int)
    tx = threadIdx().x

    # Ignore incx and incy for now.
    if tx <= n
        @inbounds dy[tx] = dx[tx]
    end
    CUDA.sync_threads()

    return
end


text:
@with_kw mutable struct LShapedData{T <: Real}
    Q::T = 1e10
    θ::T = -1e10
    iterations::Int = 0
end

@with_kw mutable struct LShapedParameters{T <: Real}
    τ::T = 1e-6
    bundle::Int = 1
    log::Bool = true
end

"""
    LShaped

Functor object for the L-shaped algorithm. Create by supplying `:ls` to the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

...
# Algorithm parameters
- `τ::Real = 1e-6`: Relative tolerance for convergence checks.
- `bundle::Int = 1`: Amount of cutting planes in bundle. A value of 1 corresponds to a multicut algorithm and a value of at least the number of scenarios yields the classical L-shaped algorithm.
- `log::Bool = true`: Specifices if L-shaped procedure should be logged on standard output or not.
...
"""
struct LShaped{F, T <: Real, A <: AbstractVector, SP <: StochasticProgram, M <: LQSolver, S <: LQSolver} <: AbstractLShapedSolver{F,T,A,M,S}
    stochasticprogram::SP
    solverdata::LShapedData{T}

    # Master
    mastersolver::M
    mastervector::A
    c::A
    x::A
    Q_history::A

    # Subproblems
    nscenarios::Int
    subproblems::Vector{SubProblem{F,T,A,S}}
    subobjectives::A

    # Cuts
    θs::A
    cuts::Vector{SparseHyperPlane{T}}
    θ_history::A

    # Params
    parameters::LShapedParameters{T}
    progress::ProgressThresh{T}

    function (::Type{LShaped})(stochasticprogram::StochasticProgram, x₀::AbstractVector, mastersolver::MPB.AbstractMathProgSolver, subsolver::MPB.AbstractMathProgSolver, F::Bool; kw...)
        if nworkers() > 1
            @warn "There are worker processes, consider using distributed version of algorithm"
        end
        first_stage = StochasticPrograms.get_stage_one(stochasticprogram)
        length(x₀) != first_stage.numCols && error("Incorrect length of starting guess, has ", length(x₀), " should be ", first_stage.numCols)

        T = promote_type(eltype(x₀), Float32)
        c_ = convert(AbstractVector{T}, JuMP.prepAffObjective(first_stage))
        c_ *= first_stage.objSense == :Min ? 1 : -1
        x₀_ = convert(AbstractVector{T}, copy(x₀))
        mastervector = convert(AbstractVector{T}, copy(x₀))
        A = typeof(x₀_)
        SP = typeof(stochasticprogram)
        msolver = LQSolver(first_stage, mastersolver)
        M = typeof(msolver)
        S = LQSolver{typeof(MPB.LinearQuadraticModel(subsolver)),typeof(subsolver)}
        n = StochasticPrograms.nscenarios(stochasticprogram)

        lshaped = new{F,T,A,SP,M,S}(stochasticprogram,
                                    LShapedData{T}(),
                                    msolver,
                                    mastervector,
                                    c_,
                                    x₀_,
                                    A(),
                                    n,
                                    Vector{SubProblem{F,T,A,S}}(),
                                    A(),
                                    A(),
                                    Vector{SparseHyperPlane{T}}(),
                                    A(),
                                    LShapedParameters{T}(;kw...),
                                    ProgressThresh(1.0, "L-Shaped Gap "))
        # Initialize solver
        init!(lshaped, subsolver)
        return lshaped
    end
end
LShaped(stochasticprogram::StochasticProgram, mastersolver::MPB.AbstractMathProgSolver, subsolver::MPB.AbstractMathProgSolver, checkfeas::Bool; kw...) = LShaped(stochasticprogram, rand(decision_length(stochasticprogram)), mastersolver, subsolver, checkfeas; kw...)

function (lshaped::LShaped)()
    # Reset timer
    lshaped.progress.tfirst = lshaped.progress.tlast = time()
    # Start procedure
    while true
        status = iterate!(lshaped)
        if status != :Valid
            return status
        end
    end
end


text:
mutable struct MixedDotAcc{T}
    s :: T
end

mixedDotAcc(::Type{Float32})                   = MixedDotAcc(Float64(0))
mixedDotAcc(::Type{Vec{N, Float32}}) where {N} = MixedDotAcc(vzero(Vec{N, Float64}))

function add!(acc::MixedDotAcc, x, y)
    SIMDops.@fusible acc.s += SIMDops.f64(x) * SIMDops.f64(y)
end

function add!(acc::MixedDotAcc{T}, x::MixedDotAcc{T}) where {T}
    SIMDops.@fusible acc.s += x.s
end

Base.sum(acc::MixedDotAcc{T}) where {T<:Vec}  = MixedDotAcc(vsum(acc.s))
Base.sum(acc::MixedDotAcc{T}) where {T<:Real} = acc.s