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

name = "LibAMVW"
version = v"1.0.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://people.cs.kuleuven.be/~raf.vandebril/homepage/software/AMVW.tar", "37000a9a5a83677fc18203bc6ae81657d74103b11eeb2257d3e68802731a85ae"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
FC=gfortran
LBLAS="-lopenblas"
FFLAGS="-shared -fPIC -O3"

patch doubleshift/src/DAMVW.f90 < damvw.patch
patch singleshift/src/zamvw.f90 < zamvw.patch

$FC $FFLAGS \
    doubleshift/src/DAMVW.f90 \
    doubleshift/src/DCB.f90 \
    doubleshift/src/DCDB.f90 \
    doubleshift/src/DCFD.f90 \
    doubleshift/src/DCFT.f90 \
    doubleshift/src/DFCC.f90 \
    doubleshift/src/DFGR.f90 \
    doubleshift/src/DGR.f90 \
    doubleshift/src/DMQF.f90 \
    doubleshift/src/DNORMALPOLY.f90 \
    doubleshift/src/DRANDPOLYJT.f90 \
    doubleshift/src/RESCHECK.f90 \
    doubleshift/src/balance.f90 \
    doubleshift/src/init_random_seed.f90 \
    doubleshift/src/turnovers/DGTO2.f90 \
    -o ${libdir}/libamvwdouble.${dlext} ${LBLAS}

$FC $FFLAGS \
    singleshift/src/balance.f90 \
    singleshift/src/buildbulge.f90 \
    singleshift/src/chasebulge.f90 \
    singleshift/src/cnormalpoly.f90 \
    singleshift/src/crgivens.f90 \
    singleshift/src/deflation.f90 \
    singleshift/src/diagblock.f90 \
    singleshift/src/factor.f90 \
    singleshift/src/fuse.f90 \
    singleshift/src/init_random_seed.f90 \
    singleshift/src/modified_quadratic.f90 \
    singleshift/src/normalpoly.f90 \
    singleshift/src/rescheck.f90 \
    singleshift/src/throughdiag.f90 \
    singleshift/src/zamvw.f90 \
    singleshift/src/zamvw2.f90 \
    singleshift/src/turnovers/dto4.f90 \
    -o ${libdir}/libamvwsingle.${dlext} $LBLAS
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:x86_64, libc=:musl),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:x86_64)
]
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libamvwdouble", :libamvwdouble)
    LibraryProduct("libamvwsingle", :libamvwsingle)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using Flux, Zygote, ForwardDiff, Random, Optim
using DifferentialEquations, OrdinaryDiffEq
using Plots
using Flux.Optimise: update!
using BSON: @save, @load
using ProgressBars
using Printf
using LinearAlgebra
using Statistics
using DiffEqFlux

n_epoch = 1000
nplot = 100
npoint = 100
ns = 5
opt = ADAMW(0.005, (0.9, 0.999), 1.f-6)

function f(y, k, t)
    dydt = [-2 * k[1] * y[1]^2 - k[2] * y[1],
            k[1] * y[1]^2 - k[4] * y[2] * y[4],
            k[2] * y[1] - k[3] * y[3],
            k[3] * y[3] - k[4] * y[2] * y[4],
            k[4] * y[2] * y[4]]
    return dydt
end

u0 = Float32[1, 1, 0, 0, 0]
k = Float32[0.1, 0.2, 0.13, 0.3]
tspan = Float32[0, 40.0]
dt = tspan[2] / npoint
ts = tspan[1]:dt:tspan[2]

prob = ODEProblem(f, u0, tspan, k)
sol = solve(prob, saveat=ts);

plt = scatter(sol, yscale=:identity);
png(plt, "figs/sol");

hls = 8;
phi_ = Flux.Chain(Dense(1, hls, tanh),
                Dense(hls, hls, tanh),
                Dense(hls, hls, tanh),
                Dense(hls, hls, tanh),
                Dense(hls, ns))

θ, re = Flux.destructure(phi_)

function phi(t, θ)
    return u0 .+ t .* re(θ)([t])
end
phi(1, θ)

dfdx = (t, θ) -> ForwardDiff.derivative(t -> phi(t, θ), t)
dfdx(1, θ)

function inner_loss(t, θ)
    mean(abs2, dfdx(t, θ) - f(phi(t, θ), k, t))
end

loss(θ) = mean(abs2, [inner_loss(t, θ) for t in ts])

l_loss_train = []
l_loss_val = []
l_grad = []
iter = 1
cb = function (p, loss_train, loss_val, g_norm)
    global iter
    global l_loss_train, l_loss_val, l_grad, iter
    push!(l_loss_train, loss_train)
    push!(l_loss_val, loss_val)
    push!(l_grad, g_norm)

    if iter % nplot == 0
        pred = hcat([phi(t, θ) for t in ts]...)
        l_plt = []
        for i in 1:ns
            plt = scatter(ts, sol[i, :], label="label");
            plot!(plt, ts, pred[i, :], label="PINN")
            xlabel!(plt, "Time")
            ylabel!(plt, "x$i")
            push!(l_plt, plt)
        end
        plt_all = plot(l_plt..., legend=:best, size=(800, 800))
        png(plt_all, "figs/pinn")

        plt_loss = plot(l_loss_train, xscale=:identity, yscale=:log10, label="Training")
        plot!(plt_loss, l_loss_val, xscale=:identity, yscale=:log10, label="Validation")
        plt_grad = plot(l_grad, xscale=:identity, yscale=:log10, label="grad_norm")
        xlabel!(plt_loss, "Epoch")
        xlabel!(plt_grad, "Epoch")
        ylabel!(plt_loss, "Loss")
        ylabel!(plt_grad, "Gradient Norm")
        plt_all = plot([plt_loss, plt_grad]..., legend=:top)
        png(plt_all, "figs/loss_grad")
    end
    iter = iter + 1
end

epochs = ProgressBar(1:n_epoch);
for epoch in epochs
    global θ
    l = loss(θ)
    grad = ForwardDiff.gradient(x -> loss(x), θ)
    update!(opt, θ, grad)
    set_description(epochs, string(@sprintf("Loss %.4e", l)))
    cb(θ, l, l, norm(grad, 2))
end

ext: jl


text:
include("euler/euler.jl")
using .Calculus: is_pentagonal
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute()::Int
    pentagonals = [1]
    j = 1
    for i ∈ pentagonals
        for k ∈ pentagonals
            if is_pentagonal(i - k) && is_pentagonal(abs(i - 2k))
                return abs(i - 2k)
            end
        end
        j += 3
        append!(pentagonals, i + j)
    end
end

compute()

@benchmark compute()

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


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

using Markdown
using InteractiveUtils

# ╔═╡ 565f87b8-b04e-11eb-38ac-bfec45a98555
using Revise, ReproducePlotUtils, Plots, FileIO, StatsPlots

# ╔═╡ 321dfa1c-7cec-465c-abeb-e8d518487863
using Reproduce

# ╔═╡ 0ce48396-7229-4801-ba66-d9e1d05563c6
ic, diff_dict =
		ReproducePlotUtils.load_data("example_data/dir_tmaze_er_rnn_rmsprop_10/")

# ╔═╡ ef1c7adc-dafe-4543-83b0-1295b31d2b10
FileIO.load(joinpath(ic[1].folder_str, "results.jld2"))["results"]

# ╔═╡ 51c42faa-c44f-47f9-b021-c4fbc3f001df
# data_col = ReproducePlotUtils.get_line_data_for(
#                ic,
#                ["numhidden", "truncation", "cell", "eta"],
#                [];
#                comp=:max,
#                get_comp_data=(x)->ReproducePlotUtils.get_AUE(x, :successes),
#                get_data=(x)->ReproducePlotUtils.get_AUE(x, :successes))

# ╔═╡ 17614876-687b-4b29-9bb2-65736f837ca3
data_col = ReproducePlotUtils.get_line_data_for(
               ic,
               ["numhidden", "truncation", "cell"],
               ["eta"];
               comp=:max,
               get_comp_data=(x)->ReproducePlotUtils.get_MUE(x, :successes),
               get_data=(x)->ReproducePlotUtils.get_rolling_mean_line(x, :successes, 300))

# ╔═╡ 4a520cfd-ed7d-4bf8-8c26-587e99e998e2
plot(data_col, Dict("numhidden"=>10, "truncation"=>12), palette=ReproducePlotUtils.custom_colorant, legend=:topleft)

# ╔═╡ 9bd94861-8f2a-4949-b5b8-14ffdc760341
data_col_sens = ReproducePlotUtils.get_line_data_for(
               ic,
               ["numhidden", "truncation", "cell", "eta"],
               [];
               comp=:max,
               get_comp_data=(x)->ReproducePlotUtils.get_MUE(x, :successes),
               get_data=(x)->ReproducePlotUtils.get_MUE(x, :successes))

# ╔═╡ e6762f85-8922-4b65-97a2-4e83b17eb89a
let
	plot(data_col_sens,
	 	 Dict("numhidden"=>10, "truncation"=>12, "cell"=>"MAGRU");
		 sort_idx="eta", z=1.97, lw=2, palette=ReproducePlotUtils.custom_colorant)
	plot!(data_col_sens,
	 	  Dict("numhidden"=>17, "truncation"=>20, "cell"=>"GRU");
	 	  sort_idx="eta", z=1.97, lw=2, palette=ReproducePlotUtils.custom_colorant)
    plot!(data_col_sens,
	 	  Dict("numhidden"=>20, "truncation"=>12, "cell"=>"AAGRU");
	 	  sort_idx="eta", z=1.97, lw=2, palette=ReproducePlotUtils.custom_colorant)
end

# ╔═╡ ff833f3f-7feb-4f6d-a5ae-0abcc3de4c5f
data_col_eps = ReproducePlotUtils.get_line_data_for(
               ic,
               ["numhidden", "truncation", "cell"],
               ["eta"];
               comp=:max,
               get_comp_data=(x)->ReproducePlotUtils.get_MUE(x, :successes),
               get_data=(x)->ReproducePlotUtils.get_extended_line(x, :successes, :total_steps, n=500))

# ╔═╡ e5331c45-6171-45a9-a537-9e3c64c17fd6
plot(data_col_eps, Dict("numhidden"=>15, "truncation"=>20, "cell"=>"MAGRU"), palette=ReproducePlotUtils.custom_colorant, legend=:topleft, z=1.0)

# ╔═╡ 3a23313d-ed57-4faa-b2b0-bec1d549c255
data_col_box = ReproducePlotUtils.get_line_data_for(
               ic,
               ["numhidden", "truncation", "cell"],
               ["eta"];
               comp=:max,
               get_comp_data=(x)->ReproducePlotUtils.get_MUE(x, :successes),
               get_data=(x)->ReproducePlotUtils.get_MUE(x, :successes))

# ╔═╡ 5e461d4b-98a5-4305-a579-78e59c4c62e7
let
	args = [
		Dict("numhidden"=>20, "truncation"=>20, "cell"=>"GRU"),
		Dict("numhidden"=>20, "truncation"=>12, "cell"=>"AAGRU"),
		Dict("numhidden"=>10, "truncation"=>12, "cell"=>"MAGRU")
		]
	boxplot(data_col_box, args;
			label_idx="cell",
			legend=false,
			palette=ReproducePlotUtils.custom_colorant)
end

# ╔═╡ 78a8f6b0-6425-434b-b2be-a30df9c2df39
let
	params = Dict("numhidden"=>20, "truncation"=>20, "cell"=>"GRU")
	idx = findfirst(data_col_box.data) do ld
		all(ld.line_params[i] == params[i] for i in keys(params))
	end
	data_col_box[idx].data, data_col_box[idx].data_pms
end

# ╔═╡ 85c8a929-5c1d-4375-b2a0-a5f6f1db4e4b
let
	args = [
		Dict("numhidden"=>20, "truncation"=>12, "cell"=>"RNN"),
		Dict("numhidden"=>20, "truncation"=>12, "cell"=>"AARNN"),
		Dict("numhidden"=>10, "truncation"=>12, "cell"=>"MARNN"),
		Dict("numhidden"=>20, "truncation"=>12, "cell"=>"GRU"),
		Dict("numhidden"=>20, "truncation"=>12, "cell"=>"AAGRU"),
		Dict("numhidden"=>10, "truncation"=>12, "cell"=>"MAGRU")
		]
	# args = Dict("numhidden"=>10, "truncation"=>12)
	plt = boxplot(data_col_box, args;
			label_idx="cell",
			legend=false,
			palette=ReproducePlotUtils.custom_colorant,
			trim=true)
	# boxplot!(data_col_box, args;
	# 		label_idx="cell",
	# 		legend=false,
	# 		palette=ReproducePlotUtils.custom_colorant,
	# 		show_median=true, lw=2, outliers=true, range=1.5)
	savefig("dir_tmaze_box_plot.pdf")
	plt
end

# ╔═╡ e564bfe9-7c26-49f8-a956-42b11a706ba1
let
	args = [
		Dict("numhidden"=>20, "truncation"=>12, "cell"=>"GRU"),
		Dict("numhidden"=>20, "truncation"=>12, "cell"=>"AAGRU"),
		Dict("numhidden"=>10, "truncation"=>12, "cell"=>"MAGRU")
		]
	violin(data_col_box, args;
			label_idx="cell",
			legend=false,
			palette=ReproducePlotUtils.custom_colorant)
end

# ╔═╡ eb41c814-0e4b-4087-b35c-1a5626c87167
let
	get_data = (x)->ReproducePlotUtils.get_MUE(x, :successes)
	sub_ic = search(ic, Dict("numhidden"=>20, "truncation"=>20, "cell"=>"GRU", "eta"=>7.8125e-5))
	diff_dict = diff(sub_ic)
	tmp = get_data(FileIO.load(joinpath(sub_ic[1].folder_str, "results.jld2")))
	d = typeof(tmp)[]
	pms_template = NamedTuple(Symbol(k)=>diff_dict[k][1] for k ∈ keys(diff_dict))
	pms = typeof(pms_template)[]
	for (idx, item) ∈ enumerate(sub_ic)
		push!(d, get_data(FileIO.load(joinpath(item.folder_str, "results.jld2"))))
		push!(pms, NamedTuple(Symbol(k)=>item.parsed_args[k] for k ∈ keys(diff_dict)))
	end
	d, pms
end

# ╔═╡ f5524caa-825e-444f-9f4a-627f8c7b67d3

# ╔═╡ Cell order:
# ╠═565f87b8-b04e-11eb-38ac-bfec45a98555
# ╠═321dfa1c-7cec-465c-abeb-e8d518487863
# ╠═0ce48396-7229-4801-ba66-d9e1d05563c6
# ╠═ef1c7adc-dafe-4543-83b0-1295b31d2b10
# ╠═51c42faa-c44f-47f9-b021-c4fbc3f001df
# ╠═17614876-687b-4b29-9bb2-65736f837ca3
# ╠═4a520cfd-ed7d-4bf8-8c26-587e99e998e2
# ╠═9bd94861-8f2a-4949-b5b8-14ffdc760341
# ╠═e6762f85-8922-4b65-97a2-4e83b17eb89a
# ╠═ff833f3f-7feb-4f6d-a5ae-0abcc3de4c5f
# ╠═e5331c45-6171-45a9-a537-9e3c64c17fd6
# ╠═3a23313d-ed57-4faa-b2b0-bec1d549c255
# ╠═5e461d4b-98a5-4305-a579-78e59c4c62e7
# ╠═78a8f6b0-6425-434b-b2be-a30df9c2df39
# ╠═85c8a929-5c1d-4375-b2a0-a5f6f1db4e4b
# ╠═e564bfe9-7c26-49f8-a956-42b11a706ba1
# ╠═eb41c814-0e4b-4087-b35c-1a5626c87167
# ╠═f5524caa-825e-444f-9f4a-627f8c7b67d3


text:
@testset "Math to Symbolics conversions" begin
    @variables A B C D Time

    test = SBML.MathApply(
        "*",
        SBML.Math[
            SBML.MathApply(
                "+",
                SBML.Math[
                    SBML.MathApply(
                        "*",
                        SBML.Math[SBML.MathIdent("A"), SBML.MathIdent("B")],
                    ),
                    SBML.MathApply(
                        "-",
                        SBML.Math[SBML.MathApply(
                            "*",
                            SBML.Math[SBML.MathIdent("C"), SBML.MathIdent("D")],
                        )],
                    ),
                ],
            ),
            SBML.MathTime("Time"),
        ],
    )

    @test isequal(convert(Num, test), (A * B - C * D) * Time)

    test = SBML.MathApply(
        "piecewise",
        SBML.Math[
            SBML.MathVal(123),
            SBML.MathApply("lt", SBML.Math[SBML.MathVal(1), SBML.MathVal(0)]),
            SBML.MathVal(456),
        ],
    )

    @test isequal(convert(Num, test), 456)
end


text:
export monitor_functions, monitor_function, get_active, set_active!

"""
    $(TYPEDEF)

An individual monitor function.
"""
struct MonitorFunction{F}
    f::F
    initial_active::Base.RefValue{Bool}
end

# Slightly annoying special casing of the different combinations of inputs to account for
# the fact that Tuples are unwrapped if they are singletons

function (mfunc::MonitorFunction)(res, u::Tuple, data::Tuple, prob...)
    res[1] =
        mfunc.f(Base.tail(u), Base.tail(data), prob...) -
        (isempty(u[1]) ? data[1][] : u[1][1])
    return nothing
end

function (mfunc::MonitorFunction)(res, u::Tuple{<:Any,<:Any}, data::Tuple, prob...)
    res[1] = mfunc.f(u[2], Base.tail(data), prob...) - (isempty(u[1]) ? data[1][] : u[1][1])
    return nothing
end

function (mfunc::MonitorFunction)(res, u::Tuple, data::Tuple{<:Any,<:Any}, prob...)
    res[1] = mfunc.f(Base.tail(u), data[2], prob...) - (isempty(u[1]) ? data[1][] : u[1][1])
    return nothing
end

function (mfunc::MonitorFunction)(
    res, u::Tuple{<:Any,<:Any}, data::Tuple{<:Any,<:Any}, prob...
)
    res[1] = mfunc.f(u[2], data[2], prob...) - (isempty(u[1]) ? data[1][] : u[1][1])
    return nothing
end

function (mfunc::MonitorFunction)(res, u::Tuple, data::Tuple{<:Any}, prob...)
    res[1] = mfunc.f(Base.tail(u), prob...) - (isempty(u[1]) ? data[] : u[1][1])
    return nothing
end

function (mfunc::MonitorFunction)(res, u::Tuple{<:Any,<:Any}, data::Tuple{<:Any}, prob...)
    res[1] = mfunc.f(u[2], prob...) - (isempty(u[1]) ? data[] : u[1][1])
    return nothing
end

function monitor_function(
    name,
    f,
    var=(),
    data=();
    initial_value=nothing,
    active=false,
    group=:embedded,
    pass_problem=false,
    top_level=true,
)
    mvar = Var(name; initial_dim=1, initial_u=initial_value, top_level=top_level)
    mdata = Data("mfunc_data", Ref(initial_value))
    fullgroup = (group isa Symbol ? push! : append!)([:mfunc], group)
    return Func(
        name,
        MonitorFunction(f, Ref(active)),
        (mvar, var...),
        (mdata, data...);
        initial_dim=1,
        group=fullgroup,
        pass_problem=pass_problem,
    )
end

function parameter(name, var; active=false, top_level=true, index=1)
    local mfunc
    let index = index
        mfunc = monitor_function(
            name, u -> u[index], (var,); active=active, top_level=top_level
        )
    end
end

function parameters(names, var; kwargs...)
    return [parameter(name, var; index=i, kwargs...) for (i, name) in enumerate(names)]
end

"""
    $(SIGNATURES)

Set a monitor function active or inactive. Only used during the creation of the problem
structure; once the continuation problem is closed, changing this has no effect.
"""
function set_active!(mfunc::Func, active::Bool)
    if mfunc.func isa MonitorFunction
        mfunc.func.initial_active[] = active
    else
        throw(ArgumentError("Func provided is not a MonitorFunction"))
    end
end

"""
    $(SIGNATURES)

Return whether a monitor function is active (able to change value) or not.
"""
function get_active(mfunc::Func)
    if mfunc.func isa MonitorFunction
        return mfunc.func.initial_active[]
    else
        throw(ArgumentError("Func provided is not a MonitorFunction"))
    end
end

"""
    $(TYPEDEF)

A problem structure to contain monitor functions.
"""
struct MonitorFunctions <: ProblemOwner
    mfunc::Vector{MonitorFunction}  # MonitorFunction is an abstract type but that's fine since it's not used in a hot loop
    idx_var::Vector{Int64}
    idx_data::Vector{Int64}
    idx_func::Vector{Int64}
    name::Dict{String,Int64}
end
function MonitorFunctions()
    return MonitorFunctions(
        MonitorFunction[], Int64[], Int64[], Int64[], Dict{String,Int64}()
    )
end

function (mfuncs::MonitorFunctions)(::Signal{:post_correct}, problem)
    # Store the var value in data (if var is non-empty)
    chart = get_current_chart(problem)
    u = get_uview(chart)
    data = get_data(chart)
    for (idx_var, idx_data) in zip(mfuncs.idx_var, mfuncs.idx_data)
        if !isempty(u[idx_var])
            data[idx_data][] = u[idx_var]
        end
    end
    return nothing
end

function (mfuncs::MonitorFunctions)(::Signal{:initial_state}, problem, u, t, data)
    T = eltype(eltype(u))
    for (mfunc, idx_var) in zip(mfuncs.mfunc, mfuncs.idx_var)
        if !mfunc.initial_active[]
            u[idx_var] = T[]
            t[idx_var] = T[]
        end
    end
    return nothing
end

function init!(mfuncs::MonitorFunctions, problem)
    flat = get_flatproblem(problem)
    # Find all monitor functions and store their locations
    for (i, func) in pairs(flat.func)
        if func.func isa MonitorFunction
            idx_var = findfirst(==(func.var[1]), flat.var)
            idx_data = findfirst(==(func.data[1]), flat.data)
            push!(mfuncs.mfunc, func.func)
            push!(mfuncs.idx_var, idx_var)
            push!(mfuncs.idx_data, idx_data)
            push!(mfuncs.idx_func, i)
            mfuncs.name[first(k for (k, v) in flat.var_names if v == idx_var)] = lastindex(mfuncs.mfunc)
        end
    end
end

function exchange_pars(flat::FlatProblem) end

"""
    $(SIGNATURES)

Create a problem structure to contain monitor functions
"""
function monitor_functions(name="mfuncs")
    problem = Problem(name, MonitorFunctions())
    return problem
end


text:
using ArgCheck
using LinearAlgebra
using DocStringExtensions
using .Scattering
import QuadGK

"""
Contains parameters for running K-Matrix calculations.

Takes as input the mesh size, equivalent to the number of
iteration steps. The time required to solve KMatrix(N) goes as 𝒪(N²).
"""
struct KMatrix <: Method
    N::Int  # Number of iteration steps
end

"""
Compute the phase shift of S-waves by computing the reactance matrix.

$(SIGNATURES)
where `k₀` is the momentum for which to compute, `m` is the mass, and `P` the potential.

# Details

The goal is to obtain the S-wave phase shift of the potential P for the momentum `k`.
This is done by solving the matrix equation
             `K = A⁻¹V`
where K is the reactance matrix. V is the matrix of the potential evaluated at the mesh points,
and A describes the Lippmann-Schwinger equation for K. The phase shift is then related to
the last diagonal element of K as
             `K(k₀, k₀) = 1/mk₀ ⋅ tan(δₗ₌₀(k₀))`
"""
function (method::KMatrix)(k₀, m, P::Potential)
    A, V = createA(k₀, m, method.N, P)
    # Create the reactance matrix K using the matrix equation AK = V
    K = inv(A)*V
    @assert size(A) == size(V) == size(K) == (method.N+1, method.N+1)

    # Diagonal elements of K are related to the phase shifts δₗ
    # as K(k₀, k₀) = -tan(δₗ)/mk₀
    Kk0 = K[end, end]
    δᵢ = atan(-Kk0*m*k₀) / pi * 180
    #energy = k₀^2 / m * 197
    δᵢ
end

"""
Construct the A matrix involved in K-Matrix calculations.

$(SIGNATURES)

where `k₀` is the momentum for which to compute, `m` is the mass, `N` the
mesh size and `P` the potential.

# Details

The goal is the set up the matrices V and A of the matrix equation

    AK = V

In the first phase, the mesh points are constructed by Gaussian quadrature,
and V is constructed as the potential evaluated at P(kᵢ, kⱼ) ∀ i,j in the mesh.
The set up of A is more complicated, and is delegated to another function.
"""
function createA(k₀, m, N::Int, P::Potential)
    k, ω = QuadGK.gauss(N)
    transform!(k, ω)
    @assert k₀ ∉ k "k₀ can not be in the mesh points k"
    push!(k, k₀)
    V = [P(ki, kj) for ki in k, kj in k]
    A = createA(V, k, ω, m)
    return A, V
end

"""
Construct the A matrix involved in K-Matrix calculations.

$(SIGNATURES)

where `V` matrix of the potential evaluated at the mesh points, `k` is the momenta,
`ω` are the quadrature weights, and `m` is the mass.
# Details

The goal is the set up the matrices V and A of the matrix equation

    K = VA.

Here A is constructed as described in
https://manybodyphysics.github.io/NuclearForces/doc/pub/scatteringtheory/html/scatteringtheory.html
"""
function createA(V, k, ω, m)
    @argcheck length(k) - 1 == length(ω)
    @argcheck size(V, 1) == size(V, 2) == length(k)
    @argcheck m > 0

    N = length(k)-1
    k₀ = k[end]
    A = diagm(0 => repeat([1.0,], N+1))  # δᵢⱼ
    uⱼ = 0.0
    @inbounds for j in 1:N+1
        uⱼ = 0.0
        if j ≠ N+1
            uⱼ = 2/π * ω[j]*k[j]^2/((k₀^2 - k[j]^2)/m)
        else
            for n in 1:N
                uⱼ += ω[n]*k₀^2/((k₀^2 - k[n]^2)/m)
            end
            uⱼ *= -2/π
        end

        for i in 1:N+1
            A[i, j] -= V[i, j] * uⱼ
        end
    end
    return A
end

"""
Map the Gaussian quadrature points from [0, 1] to (-∞, ∞)

$(SIGNATURES)

"""
function transform!(xs, weights)
    # Times 197 for MeV
    @. weights = π/4 * weights / cos(π/4*(1.0 + xs))^2
    @. xs      = tan(π/4*(1.0 + xs))
end

max_stars_repo_name: ErlendLima/Scattering

max_issues_repo_name: ErlendLima/Scattering

max_forks_repo_name: ErlendLima/Scattering


text:
struct AutoZygote <: AbstractADType end

function instantiate_function(f, x, ::AutoZygote, p, num_cons = 0)
    num_cons != 0 && error("AutoZygote does not currently support constraints")

    _f = (θ, args...) -> f(θ,p,args...)[1]
    if f.grad === nothing
        grad = (res, θ, args...) -> res isa DiffResults.DiffResult ? DiffResults.gradient!(res, Zygote.gradient(x -> _f(x, args...), θ)[1]) : res .= Zygote.gradient(x -> _f(x, args...), θ)[1]
    else
        grad = f.grad
    end

    if f.hess === nothing
        hess = function (res, θ, args...)
            if res isa DiffResults.DiffResult
                DiffResults.hessian!(res, ForwardDiff.jacobian(θ) do θ
                                                Zygote.gradient(x -> _f(x, args...), θ)[1]
                                            end)
            else
                res .=  ForwardDiff.jacobian(θ) do θ
                    Zygote.gradient(x ->_f(x, args...), θ)[1]
                end
            end
        end
    else
        hess = f.hess
    end

    if f.hv === nothing
        hv = function (H, θ, v, args...)
            _θ = ForwardDiff.Dual.(θ, v)
            res = DiffResults.GradientResult(_θ)
            grad(res, _θ, args...)
            H .= getindex.(ForwardDiff.partials.(DiffResults.gradient(res)),1)
        end
    else
        hv = f.hv
    end

    return OptimizationFunction{false,AutoZygote,typeof(f),typeof(grad),typeof(hess),typeof(hv),Nothing,Nothing,Nothing}(f,AutoZygote(),grad,hess,hv,nothing,nothing,nothing)
end


text:
using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
X = raviartthomas(Γ)

Δt, Nt = 0.6, 200

δ = timebasisdelta(Δt, Nt)
h = timebasisc0d1(Δt, Nt)
p = timebasiscxd0(Δt, Nt)

width, delay, scaling = 8.0, 12.0, 1.0
gaussian = creategaussian(width, delay, scaling)

direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, derive(gaussian), 1.0)

@assert p(-Δt/2) ≈ 1
@assert h(0) ≈ 1

rhsd = assemble(E, X⊗δ)
rhsp = assemble(E, X⊗p)
rhsh = assemble(E, X⊗h)

imax = floor(Int,delay/Δt)+5

plot()
plot!(rhsd[95,:])
scatter!(rhsp[95,:]/Δt^1, markersize=2.0)


text:
module Graph500

  # Dependencies
  using LightGraphs
  using StatsBase
  using ProgressMeter
  import LightGraphs.sample
  using Distributed: @everywhere, @spawnat, @distributed, nprocs, myid
  using Markdown: @doc_str
  using Random
  using Statistics: std

  export
  # kronecker_generator
  kronecker_generator,get_min_type,

  # kernel_1
  kernel_1, get_next_type,

  # kernel_2
  key_sampling, kernel_2,

  # bfs validation
  validate,

  # driver function
  driver,

  # print statistics
  output

  Graph500
  include("kronecker_generator.jl")
  include("kernel_1.jl")
  include("kernel_2.jl")
  include("validation.jl")
  include("driver.jl")
  include("output.jl")

end # module


text:
# missings go to the end
# fast algorithm for integers with few unique values forward direction
function ds_sort_int_missatright!(x, original_P, copy_P, where, lo, hi, rangelen, minval)
    offs = 1 - minval

    @inbounds for i in 1:rangelen+2
        where[i] = 0
    end

    # where = fill(0, rangelen+1)
    where[1] = 1
    @inbounds for i = lo:hi
        ismissing(x[i]) ? where[rangelen+2] += 1 : where[Int(x[i]) + offs + 1] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 2:rangelen+2
        where[i] += where[i-1]
    end

    @inbounds for i = lo:hi
        ismissing(x[i]) ? label = rangelen + 1 : label = Int(x[i]) + offs
        original_P[where[label] + lo - 1] = copy_P[i]
        where[label] += 1
    end
    # rearrange data
    @inbounds for i in 1:rangelen+1
        where[i] = 0
    end

    @inbounds for i = lo:hi
        ismissing(x[i]) ? where[rangelen+1] += 1 : where[Int(x[i]) + offs] += 1
        # where[x[i] + offs] += 1
    end
    f_indx = lo
    @inbounds for i = 1:rangelen
        l_indx = f_indx + where[i] - 1
        val = i-offs
        for j = f_indx:l_indx
            x[j] = val
        end
        f_indx = l_indx + 1
    end
    if f_indx <= hi
        for j = f_indx:hi
            x[j] = missing
        end
    end
end
function ds_sort_int_missatleft!(x, original_P, copy_P, where, lo, hi, rangelen, minval)
    offs = 1 - minval

    @inbounds for i in 1:rangelen+2
        where[i] = 0
    end

    # where = fill(0, rangelen+1)
    where[1] = 1
    @inbounds for i = lo:hi
        ismissing(x[i]) ? where[2] += 1 : where[Int(x[i]) + offs + 2] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 2:rangelen+2
        where[i] += where[i-1]
    end

    @inbounds for i = lo:hi
        ismissing(x[i]) ? label = 1 : label = Int(x[i]) + offs + 1
        original_P[where[label] + lo - 1] = copy_P[i]
        where[label] += 1
    end
    # rearrange data
    @inbounds for i in 1:rangelen+1
        where[i] = 0
    end

    @inbounds for i = lo:hi
        ismissing(x[i]) ? where[1] += 1 : where[Int(x[i]) + offs + 1] += 1
        # where[x[i] + offs] += 1
    end
    f_indx = lo + where[1]
    @inbounds for i = 2:rangelen+1
        l_indx = f_indx + where[i] - 1
        val = i - 1 - offs
        for j = f_indx:l_indx
            x[j] = val
        end
        f_indx = l_indx + 1
    end
    for j = lo:(lo + where[1] - 1)
        x[j] = missing
    end
end

# to simplify the problem we assume number_of_chunks is 2^n for some n
function _sort_chunks_int_right!(x, idx::Vector{<:Integer}, idx_cpy, where, number_of_chunks, rangelen, minval, o::Ordering)
    cz = div(length(x), number_of_chunks)
    en = length(x)
    Threads.@threads for i in 1:number_of_chunks
        ds_sort_int_missatright!(x, idx, idx_cpy, where[Threads.threadid()], (i-1)*cz+1,i*cz, rangelen, minval)
    end
    # take care of the last few observations
    if number_of_chunks*div(length(x), number_of_chunks) < en
        ds_sort_int_missatright!(x, idx, idx_cpy, where[Threads.threadid()],  number_of_chunks*div(length(x), number_of_chunks)+1, en, rangelen, minval)
    end
end

# to simplify the problem we assume number_of_chunks is 2^n for some n
function _sort_chunks_int_left!(x, idx::Vector{<:Integer}, idx_cpy, where, number_of_chunks, rangelen, minval, o::Ordering)
    cz = div(length(x), number_of_chunks)
    en = length(x)
    Threads.@threads for i in 1:number_of_chunks
        ds_sort_int_missatleft!(x, idx, idx_cpy, where[Threads.threadid()], (i-1)*cz+1,i*cz, rangelen, minval)
    end
    # take care of the last few observations
    if number_of_chunks*div(length(x), number_of_chunks) < en
        ds_sort_int_missatleft!(x, idx, idx_cpy, where[Threads.threadid()],  number_of_chunks*div(length(x), number_of_chunks)+1, en, rangelen, minval)
    end
end

# sorting a vector using parallel quick sort
# it uses a simple algorithm for doing this, and to make it even simpler the number of threads must be in the form of 2^n
function hp_ds_sort_int!(x, idx, idx_cpy, where, rangelen, minval, missatleft, a::Base.Sort.Algorithm, o::Ordering)
    cpucnt = Threads.nthreads()
    @assert cpucnt >= 2 "we need at least 2 cpus for parallel sorting"
    cpucnt = 2 ^ floor(Int, log2(cpucnt))
    if missatleft
        _sort_chunks_int_left!(x , idx, idx_cpy, where, cpucnt, rangelen, minval, o)
    else
        _sort_chunks_int_right!(x , idx, idx_cpy, where, cpucnt, rangelen, minval, o)
    end
    _sort_multi_sorted_chunk!(x, idx, cpucnt, a, o)
end

function _cleanup_starts!(starts, sz)
    curloc=2
    i = 1
    @inbounds while true
        if starts[curloc] == starts[i]
            curloc += 1
        else
            i += 1
            starts[i] = starts[curloc]

        end
        starts[i]>sz && break
    end
    return resize!(starts, i-1)
end

function _ds_sort_int_missatright_nopermx!(x, original_P, copy_P, lo, hi, rangelen, minval, ::Val{T}) where T
    offs = 1 - minval
    where = zeros(T, rangelen + 3)

    # where = fill(0, rangelen+1)
    where[1] = 1
    where[2] = 1
    @inbounds for i = lo:hi
        ismissing(x[i]) ? where[rangelen+3] += 1 : where[Int(x[i]) + offs + 2] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[i] += where[i-1]
    end

    @inbounds for i = lo:hi
        ismissing(x[i]) ? label = rangelen + 2 : label = Int(x[i]) + offs + 1
        original_P[where[label] + lo - 1] = copy_P[i]
        where[label] += 1
    end
    where
end

function _ds_sort_int_missatright_nopermx!(x, original_P, rangelen, minval, ::Val{T}) where T
    offs = 1 - minval
    where = zeros(T, rangelen + 3)

    # where = fill(0, rangelen+1)
    where[1] = 1
    where[2] = 1
    @inbounds for i = 1:length(x)
        ismissing(x[i]) ? where[rangelen+3] += 1 : where[Int(x[i]) + offs + 2] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[i] += where[i-1]
    end

    @inbounds for i = 1:length(x)
        ismissing(x[i]) ? label = rangelen + 2 : label = Int(x[i]) + offs + 1
        original_P[where[label]] = i
        where[label] += 1
    end
    where
end

function _ds_sort_int_missatleft_nopermx!(x, original_P, copy_P, lo, hi, rangelen, minval, ::Val{T}) where T
    offs = 1 - minval
    where = zeros(T, rangelen + 3)

    # where = fill(0, rangelen+1)
    where[1] = 1
    where[2] = 1
    @inbounds for i = lo:hi
        ismissing(x[i]) ? where[3] += 1 : where[Int(x[i]) + offs + 3] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[i] += where[i-1]
    end

    @inbounds for i = lo:hi
        ismissing(x[i]) ? label = 2 : label = Int(x[i]) + offs + 2
        original_P[where[label] + lo - 1] = copy_P[i]
        where[label] += 1
    end
    where
end

function _ds_sort_int_missatleft_nopermx!(x, original_P, rangelen, minval, ::Val{T}) where T
    offs = 1 - minval
    where = zeros(T, rangelen + 3)

    # where = fill(0, rangelen+1)
    where[1] = 1
    where[2] = 1
    @inbounds for i = 1:length(x)
        ismissing(x[i]) ? where[3] += 1 : where[Int(x[i]) + offs + 3] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[i] += where[i-1]
    end

    @inbounds for i = 1:length(x)
        ismissing(x[i]) ? label = 2 : label = Int(x[i]) + offs + 2
        original_P[where[label]] = i
        where[label] += 1
    end
    where
end

# threaded version

function _ds_sort_int_missatright_nopermx_threaded!(x, original_P, copy_P, lo, hi, rangelen, minval, ::Val{T}) where T
    offs = 1 - minval
    nt = Threads.nthreads()
    where = [zeros(T, rangelen + 3) for _ in 1:nt]

    # where = fill(0, rangelen+1)
    Threads.@threads for i in 1:nt
        where[i][1] = 1
        where[i][2] = 1
    end
    Threads.@threads for i = lo:hi
        @inbounds ismissing(x[i]) ? where[Threads.threadid()][rangelen+3] += 1 : where[Threads.threadid()][Int(x[i]) + offs + 2] += 1
    end
    for j in 3:length(where[1])
        for i in 2:nt
            where[1][j] += where[i][j]
        end
    end
    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[1][i] += where[1][i-1]
    end

    @sync for thid in 0:nt-1
        Threads.@spawn for i = lo:hi
            if ismissing(x[i])
                # missing values are handled by first thread
                if thid == 0
                    label = rangelen + 2
                    original_P[where[1][label] + lo - 1] = copy_P[i]
                    where[1][label] += 1
                end
            else
                label = Int(x[i]) + offs + 1
                if label % nt == thid
                    original_P[where[1][label] + lo - 1] = copy_P[i]
                    where[1][label] += 1
                end
            end
        end
    end
    where[1]
end

function _ds_sort_int_missatright_nopermx_threaded!(x, original_P, rangelen, minval, ::Val{T}) where T
    offs = 1 - minval
    nt = Threads.nthreads()
    where = [zeros(T, rangelen + 3) for _ in 1:nt]

    # where = fill(0, rangelen+1)
    Threads.@threads for i in 1:nt
        where[i][1] = 1
        where[i][2] = 1
    end
    Threads.@threads for i = 1:length(x)
        @inbounds ismissing(x[i]) ? where[Threads.threadid()][rangelen+3] += 1 : where[Threads.threadid()][Int(x[i]) + offs + 2] += 1
    end
    for j in 3:length(where[1])
        for i in 2:nt
            where[1][j] += where[i][j]
        end
    end
    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[1][i] += where[1][i-1]
    end
    @sync for thid in 0:nt-1
        Threads.@spawn for i = 1:length(x)
            if ismissing(x[i])
                if thid == 0
                    label = rangelen + 2
                    original_P[where[1][label]] = i
                    where[1][label] += 1
                end
            else
                label = Int(x[i]) + offs + 1
                if label % nt == thid
                    original_P[where[1][label]] = i
                    where[1][label] += 1
                end
            end
        end
    end
    where[1]
end

function _ds_sort_int_missatleft_nopermx_threaded!(x, original_P, copy_P, lo, hi, rangelen, minval, ::Val{T}) where T
    offs = 1 - minval
    nt = Threads.nthreads()
    where = [zeros(T, rangelen + 3) for _ in 1:nt]

    # where = fill(0, rangelen+1)
    Threads.@threads for i in 1:nt
        where[i][1] = 1
        where[i][2] = 1
    end
    Threads.@threads for i = lo:hi
        @inbounds ismissing(x[i]) ? where[Threads.threadid()][3] += 1 : where[Threads.threadid()][Int(x[i]) + offs + 3] += 1
    end
    for j in 3:length(where[1])
        for i in 2:nt
            where[1][j] += where[i][j]
        end
    end
    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[1][i] += where[1][i-1]
    end

    @sync for thid in 0:nt-1
        Threads.@spawn for i = lo:hi
            if ismissing(x[i])
                if thid == 0
                    label = 2
                    original_P[where[1][label] + lo - 1] = copy_P[i]
                    where[1][label] += 1
                end
            else
                label = Int(x[i]) + offs + 2
                if label % nt == thid
                    original_P[where[1][label] + lo - 1] = copy_P[i]
                    where[1][label] += 1
                end
            end
        end
    end
    where[1]

end

function _ds_sort_int_missatleft_nopermx_threaded!(x, original_P, rangelen, minval, ::Val{T}) where T
    offs = 1 - minval
    nt = Threads.nthreads()
    where = [zeros(T, rangelen + 3) for _ in 1:nt]

    # where = fill(0, rangelen+1)
    Threads.@threads for i in 1:nt
        where[i][1] = 1
        where[i][2] = 1
    end
    Threads.@threads for i = 1:length(x)
        @inbounds ismissing(x[i]) ? where[Threads.threadid()][3] += 1 : where[Threads.threadid()][Int(x[i]) + offs + 3] += 1
    end
    for j in 3:length(where[1])
        for i in 2:nt
            where[1][j] += where[i][j]
        end
    end
    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[1][i] += where[1][i-1]
    end
    @sync for thid in 0:nt-1
        Threads.@spawn for i = 1:length(x)
            if ismissing(x[i])
                if thid == 0
                    label = 2
                    original_P[where[1][label]] = i
                    where[1][label] += 1
                end
            else
                label = Int(x[i]) + offs + 2
                if label % nt == thid
                    original_P[where[1][label]] = i
                    where[1][label] += 1
                end
            end
        end
    end
    where[1]
end

# low memory threaded
function _ds_sort_int_missatright_nopermx_threaded_lm!(x, original_P, copy_P, lo, hi, rangelen, minval, ::Val{T}) where T
    nt = Threads.nthreads()

    offs = 1 - minval
    where = zeros(T, rangelen + 3)

    # where = fill(0, rangelen+1)
    where[1] = 1
    where[2] = 1
    @inbounds for i = lo:hi
        ismissing(x[i]) ? where[rangelen+3] += 1 : where[Int(x[i]) + offs + 2] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[i] += where[i-1]
    end

    @sync for thid in 0:nt-1
        Threads.@spawn for i = lo:hi
            if ismissing(x[i])
                # missing values are handled by first thread
                if thid == 0
                    label = rangelen + 2
                    original_P[where[label] + lo - 1] = copy_P[i]
                    where[label] += 1
                end
            else
                label = Int(x[i]) + offs + 1
                if label % nt == thid
                    original_P[where[label] + lo - 1] = copy_P[i]
                    where[label] += 1
                end
            end
        end
    end
    where
end

function _ds_sort_int_missatright_nopermx_threaded_lm!(x, original_P, rangelen, minval, ::Val{T}) where T
    nt = Threads.nthreads()

    offs = 1 - minval
    where = zeros(T, rangelen + 3)

    # where = fill(0, rangelen+1)
    where[1] = 1
    where[2] = 1
    @inbounds for i = 1:length(x)
        ismissing(x[i]) ? where[rangelen+3] += 1 : where[Int(x[i]) + offs + 2] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[i] += where[i-1]
    end
    @sync for thid in 0:nt-1
        Threads.@spawn for i = 1:length(x)
            if ismissing(x[i])
                if thid == 0
                    label = rangelen + 2
                    original_P[where[label]] = i
                    where[label] += 1
                end
            else
                label = Int(x[i]) + offs + 1
                if label % nt == thid
                    original_P[where[label]] = i
                    where[label] += 1
                end
            end
        end
    end
    where
end

function _ds_sort_int_missatleft_nopermx_threaded_lm!(x, original_P, copy_P, lo, hi, rangelen, minval, ::Val{T}) where T
    nt = Threads.nthreads()

    offs = 1 - minval
    where = zeros(T, rangelen + 3)

    # where = fill(0, rangelen+1)
    where[1] = 1
    where[2] = 1
    @inbounds for i = lo:hi
        ismissing(x[i]) ? where[3] += 1 : where[Int(x[i]) + offs + 3] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[i] += where[i-1]
    end

    @sync for thid in 0:nt-1
        Threads.@spawn for i = lo:hi
            if ismissing(x[i])
                if thid == 0
                    label = 2
                    original_P[where[label] + lo - 1] = copy_P[i]
                    where[label] += 1
                end
            else
                label = Int(x[i]) + offs + 2
                if label % nt == thid
                    original_P[where[label] + lo - 1] = copy_P[i]
                    where[label] += 1
                end
            end
        end
    end
    where

end

function _ds_sort_int_missatleft_nopermx_threaded_lm!(x, original_P, rangelen, minval, ::Val{T}) where T
    nt = Threads.nthreads()

    offs = 1 - minval
    where = zeros(T, rangelen + 3)

    # where = fill(0, rangelen+1)
    where[1] = 1
    where[2] = 1
    @inbounds for i = 1:length(x)
        ismissing(x[i]) ? where[3] += 1 : where[Int(x[i]) + offs + 3] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 3:rangelen+3
        where[i] += where[i-1]
    end
    @sync for thid in 0:nt-1
        Threads.@spawn for i = 1:length(x)
            if ismissing(x[i])
                if thid == 0
                    label = 2
                    original_P[where[label]] = i
                    where[label] += 1
                end
            else
                label = Int(x[i]) + offs + 2
                if label % nt == thid
                    original_P[where[label]] = i
                    where[label] += 1
                end
            end
        end
    end
    where
end


text:
export LinearMap
export SplitLayer
export CatLayer

struct LinearMap
    W::AbstractArray
end

LinearMap(in::Int, out::Int, initW=Flux.glorot_uniform) = LinearMap(initW(out, in))
(m::LinearMap)(x::AbstractArray) = m.W * x
Flux.@functor LinearMap

"""
    CatLayer(layers::Tuple)

Concatenates output of multiple layers e.g.:
x -> vcat(layer[1](x), ..., layer[end](x))
"""
struct CatLayer
    layers::Tuple
end

function CatLayer(layers...)
    CatLayer(layers)
end

function (m::CatLayer)(x)
    y = [layer(x) for layer in m.layers]
    vcat(y...)
end

Flux.@functor CatLayer


text:
# function _dot_nodes_to_dict(pg)
#     nodes = collect(DOT.nodes(pg))
#     n = length(nodes)
#     # try to convert label in index
#     # so that  write("file",g) == read("file")
#     try
#         if all(v->isalpha(first(v)), nodes)
#             d = Dict(v => parse(Int, v[2:end]) for v in nodes)
#         else
#             d = Dict(nodes[i] => parse(Int, nodes[i]) for i=1:n)
#         end
#         @assert minimum(values(d)) == 1
#         @assert maximum(values(d)) == n
#         return d
#     catch
#         return Dict(zip(nodes, 1:n))
#     end
# end

# function _readdot{G<:AGraphOrDiGraph}(pg, ::Type{G})
#     n = length(DOT.nodes(pg))
#     nodedict = _dot_nodes_to_dict(pg)
#     g = G(n)
#     for es in DOT.edges(pg)
#         s = nodedict[es[1]]
#         d = nodedict[es[2]]
#         add_edge!(g, s, d)
#     end
#     return g
# end
#
# function readdot{G<:AGraphOrDiGraph}(io::IO, ::Type{G})
#     pg = first(DOT.parse_dot(readstring(io)))
#     H = pg.directed ? digraphtype(G) : graphtype(G)
#     return _readdot(pg, H)
# end

function writedot(io::IO, g::AGraphOrDiGraph)
    if is_directed(g)
        println(io, "strict digraph {")
        eop = "->"
    else
        println(io, "strict graph {")
        eop = "--"
    end

    for i=1:nv(g)
        if degree(g, i) == 0
            println(io, "\tn$i;")
        end
    end

    for e in edges(g)
        println(io, "\tn$(src(e)) $eop n$(dst(e));")
    end

    println(io, "}")
    return 1
end

filemap[:dot] = (NI, writedot, NI, NI)


text:
module BlockMaps

using LinearAlgebra
using SparseArrays
using BlockArrays
import BlockArrays: getblock, setblock!, findblockindex, block, blockindex

using LinearMaps
using LazyArrays
using FillArrays

include("block_maps.jl")
include("sparse_block_maps.jl")
include("linear_maps.jl")

export Block

end # module


text:
function coupling!(system::LegSystem,n::Int64,h::Float64,terms::Vector{Int64},
    uprox::Float64)
    # println("Distal coupling:")
    LegHemodynamics.coupledistal!(system,n,h,terms);
    # println("Proximal coupling:")
    LegHemodynamics.coupleproximal!(system,n,h,uprox);
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
# initialization case
function promote_type2(@nospecialize(T), @nospecialize(S))
    if T === Union{} || T === S
        return S
    elseif T === Missing || S === Missing
        return Union{T, S}
    elseif T === Int64
        return S === Float64 ? Float64 : String
    elseif T === Float64
        return S === Int64 ? Float64 : String
    elseif T === Date
        return S === DateTime ? DateTime : String
    elseif T === DateTime
        return S === Date ? DateTime : String
    elseif T === String
        return String
    else
        return S === String ? Union{String, Missing} : Union{promote_type2(Base.nonmissingtype(T), Base.nonmissingtype(S)), Missing}
    end
end

# providedtypes: Dict{String, Type}, Dict{Int, Type}, Vector{Type}
initialtype(allowmissing) = (allowmissing === :auto || allowmissing === :none) ? Union{} : Missing
initialtypes(T, ::Nothing, names) = Type[T for _ = 1:length(names)]
initialtypes(T, t::Dict{String, V}, names) where {V} = Type[get(t, string(nm), T) for nm in names]
initialtypes(T, t::Dict{Symbol, V}, names) where {V} = Type[get(t, nm, T) for nm in names]
initialtypes(T, t::Dict{Int, V}, names) where {V} = Type[get(t, i, T) for i = 1:length(names)]
initialtypes(T, t::Vector, names) = length(t) == length(names) ? collect(Type, t) : throw(ArgumentError("length of user provided types ($(length(t))) does not match length of header (?$(length(names)))"))

struct Fib{N}
    len::Int
end
Fib(len) = Fib{1024}(len)

@inline function Base.iterate(f::Fib)
    f.len == 0 && return nothing
    return 1, (1, 1, 1, 1)
end

@inline function Base.iterate(f::Fib{N}, (off, len, prevfib, fib)) where {N}
    (off + fib) > f.len && return nothing
    if rem(len, N) == 0
        prevfib, fib = fib, prevfib + fib
    end
    return off + fib, (off + fib, len + 1, prevfib, fib)
end

function detect(types, io, positions, parsinglayers, kwargs, typemap, categorical, transpose, ref, debug)
    len = transpose ? ref[] : length(positions)
    cols = length(types)

    # prep if categorical
    levels = categorical !== false ? Dict{String, Int}[Dict{String, Int}() for _ = 1:cols] : Dict{String, Int}[]

    lastcode = Base.RefValue(Parsers.SUCCESS)
    rows = 0
    for row in Fib(len)
        rows += 1
        !transpose && seek(io, positions[row])
        lastcode[] = Parsers.SUCCESS
        for col = 1:cols
            if !transpose && newline(lastcode[])
                types[col] = promote_type2(types[col], Missing)
                continue
            end
            transpose && seek(io, positions[col])
            # if debug
            #     pos = position(io)
            #     debug && println("pos=$pos")
            #     result = Parsers.parse(parsinglayers, io, String)
            #     seek(io, pos)
            #     debug && println("col: $col, detecting type for: '$(result.result)'")
            # end
            T = detecttype(Base.nonmissingtype(types[col]), io, parsinglayers, kwargs, levels, row, col, categorical, lastcode)
            debug && println("col: $col, detected type: $T")
            types[col] = promote_type2(types[col], T)
            debug && println("col: $col, promoted to: $(types[col])")
            transpose && setindex!(positions, position(io), col)
        end
    end
    debug && println("scanned $rows / $len = $((rows / len) * 100)% of file to infer types")
    debug && @show types

    if categorical !== false
        pools = Vector{CategoricalPool{String, UInt32, CatStr}}(undef, cols)
        for col = 1:cols
            T = types[col]
            # FIXME: the first condition could be moved second for performance,
            # but it triggers an inference bug in Julia 0.7 and 1.0
            if (categorical === true || length(levels[col]) / sum(values(levels[col])) < categorical) &&
               (T === String || T === Union{String, Missing})
                types[col] = substitute(T, CatStr)
                pools[col] = CategoricalPool{String, UInt32}(collect(keys(levels[col])))
            end
        end
    else
        pools = CategoricalPool{String, UInt32, CatStr}[]
    end
    if !isempty(typemap)
        for col = 1:cols
            types[col] = get(typemap, types[col], types[col])
        end
    end
    return types, pools
end

@inline function trytype(io, pos, layers, T, kwargs, lastcode)
    seek(io, pos)
    res = Parsers.parse(layers, io, T; kwargs...)
    lastcode[] = res.code
    return Parsers.ok(res.code) ? typeof(res.result) : nothing
end

function detecttype(@nospecialize(prevT), io, layers, kwargs, levels, row, col, categorical, lastcode)
    pos = position(io)
    if categorical !== false
        result = Parsers.parse(layers, io, Tuple{Ptr{UInt8}, Int})
        Parsers.ok(result.code) || throw(Error(Parsers.Error(io, result), row, col))
        lastcode[] = result.code
        # @debug "res = $result"
        result.result === missing && return Missing
        # update levels
        incr!(levels[col], result.result::Tuple{Ptr{UInt8}, Int})
    end
    if prevT === Union{} || prevT === Missing || prevT === Int64
        int = trytype(io, pos, layers, Int64, kwargs, lastcode)
        int !== nothing && return int
    end
    if prevT === Union{} || prevT === Missing || prevT === Int64 || prevT === Float64
        float = trytype(io, pos, layers, Float64, kwargs, lastcode)
        float !== nothing && return float
    end
    if prevT === Union{} || prevT === Missing || prevT === Date || prevT === DateTime
        if !haskey(kwargs, :dateformat)
            date = trytype(io, pos, layers, Date, kwargs, lastcode)
            date !== nothing && return date
            try
                datetime = trytype(io, pos, layers, DateTime, kwargs, lastcode)
                datetime !== nothing && return datetime
            catch e
            end
        else
            # use user-provided dateformat
            T = timetype(kwargs.dateformat)
            dt = trytype(io, pos, layers, T, kwargs, lastcode)
            dt !== nothing && return dt
        end
    end
    if prevT === Union{} || prevT === Missing || prevT === Bool
        bool = trytype(io, pos, layers, Bool, kwargs, lastcode)
        bool !== nothing && return bool
    end
    str = trytype(io, pos, layers, Tuple{Ptr{UInt8}, Int}, kwargs, lastcode)
    return ifelse(str === nothing, String, ifelse(str === Missing, Missing, String))
end

function timetype(df::Dates.DateFormat)
    date = false
    time = false
    for token in df.tokens
        T = typeof(token)
        if T == Dates.DatePart{'H'}
            time = true
        elseif T == Dates.DatePart{'y'} || T == Dates.DatePart{'Y'}
            date = true
        end
    end
    return ifelse(date & time, DateTime, ifelse(time, Time, Date))
end

function incr!(dict::Dict{String, Int}, key::Tuple{Ptr{UInt8}, Int})
    index = Base.ht_keyindex2!(dict, key)
    if index > 0
        @inbounds dict.vals[index] += 1
        return
    else
        kk::String = convert(String, key)
        @inbounds Base._setindex!(dict, 1, kk, -index)
        return
    end
end


text:
using LightStructArrays
using Test

# write your own tests here
@testset "index" begin
    a, b = [1 2; 3 4], [4 5; 6 7]
    t = StructArray((a = a, b = b))
    @test (@inferred t[2,2]) == (a = 4, b = 7)
    @test (@inferred t[2,1:2]) == StructArray((a = [3, 4], b = [6, 7]))
    @test_throws BoundsError t[3,3]
    @test (@inferred view(t, 2, 1:2)) == StructArray((a = view(a, 2, 1:2), b = view(b, 2, 1:2)))
end

@testset "fieldarrays" begin
    t = StructArray(a = 1:10, b = rand(Bool, 10))
    @test LightStructArrays.propertynames(t) == (:a, :b)
    @test LightStructArrays.propertynames(LightStructArrays.fieldarrays(t)) == (:a, :b)
end

@testset "similar" begin
    t = StructArray(a = rand(10), b = rand(Bool, 10))
    s = similar(t)
    @test eltype(s) == NamedTuple{(:a, :b), Tuple{Float64, Bool}}
    @test size(s) == (10,)
    t = StructArray(a = rand(10, 2), b = rand(Bool, 10, 2))
    s = similar(t, 3, 5)
    @test eltype(s) == NamedTuple{(:a, :b), Tuple{Float64, Bool}}
    @test size(s) == (3, 5)
    s = similar(t, (3, 5))
    @test eltype(s) == NamedTuple{(:a, :b), Tuple{Float64, Bool}}
    @test size(s) == (3, 5)
end

@testset "empty" begin
    s = StructVector(a = [1, 2, 3], b = ["a", "b", "c"])
    empty!(s)
    @test isempty(s.a)
    @test isempty(s.b)
end

@testset "constructor from existing array" begin
    v = rand(ComplexF64, 5, 3)
    t = @inferred StructArray(v)
    @test size(t) == (5, 3)
    @test t[2,2] == v[2,2]
    t2 = convert(StructArray, v)::StructArray
    @test t2 == t
    t3 = StructArray(t)::StructArray
    @test t3 == t
end

@testset "tuple case" begin
    s = StructArray(([1], ["test"],))
    @test s[1] == (1, "test")
    @test Base.getproperty(s, 1) == [1]
    @test Base.getproperty(s, 2) == ["test"]
    t = StructArray{Tuple{Int, Float64}}([1], [1.2])
    @test t[1] == (1, 1.2)

    t[1] = (2, 3)
    @test t[1] == (2, 3.0)
    push!(t, (1, 2))
    @test getproperty(t, 1) == [2, 1]
    @test getproperty(t, 2) == [3.0, 2.0]
end

@testset "kwargs constructor" begin
    a = [1.2]
    b = [2.3]
    @test StructArray(a=a, b=b) == StructArray((a=a, b=b))
    @test StructArray{ComplexF64}(re=a, im=b) == StructArray{ComplexF64}(a, b)
    f1() = StructArray(a=[1.2], b=["test"])
    f2() = StructArray{Pair}(first=[1.2], second=["test"])
    t1 = @inferred f1()
    t2 = @inferred f2()
    @test t1 == StructArray((a=[1.2], b=["test"]))
    @test t2 == StructArray{Pair}([1.2], ["test"])
end

@testset "complex" begin
    a, b = [1 2; 3 4], [4 5; 6 7]
    t = StructArray{ComplexF64}(a, b)
    @test t[2,2] == ComplexF64(4, 7)
    @test t[2,1:2] == StructArray{ComplexF64}([3, 4], [6, 7])
    @test view(t, 2, 1:2) == StructArray{ComplexF64}(view(a, 2, 1:2), view(b, 2, 1:2))
end

@testset "copy" begin
    a, b = [1 2; 3 4], [4 5; 6 7]
    t = StructArray{ComplexF64}(a, b)
    t2 = @inferred copy(t)
    @test t2[1,1] == 1.0 + im*4.0
    t2[1,1] = 2.0 + im*4.0
    # Test we actually did a copy
    @test t[1,1] == 1.0 + im*4.0
end

@testset "undef initializer" begin
    t = @inferred StructArray{ComplexF64}(undef, 5, 5)
    @test eltype(t) == ComplexF64
    @test size(t) == (5,5)
    c = 2 + im
    t[1,1] = c
    @test t[1,1] == c
end

@testset "resize!" begin
    t = StructArray{Pair}([3, 5], ["a", "b"])
    resize!(t, 5)
    @test length(t) == 5
    p = 1 => "c"
    t[5] = p
    @test t[5] == p
end

@testset "concat" begin
    t = StructArray{Pair}([3, 5], ["a", "b"])
    push!(t, (2 => "c"))
    @test t == StructArray{Pair}([3, 5, 2], ["a", "b", "c"])
    append!(t, t)
    @test t == StructArray{Pair}([3, 5, 2, 3, 5, 2], ["a", "b", "c", "a", "b", "c"])
    t = StructArray{Pair}([3, 5], ["a", "b"])
    t2 = StructArray{Pair}([1, 6], ["a", "b"])
    @test cat(t, t2; dims=1)::StructArray == StructArray{Pair}([3, 5, 1, 6], ["a", "b", "a", "b"]) == vcat(t, t2)
    @test vcat(t, t2) isa StructArray
    @test cat(t, t2; dims=2)::StructArray == StructArray{Pair}([3 1; 5 6], ["a" "a"; "b" "b"]) == hcat(t, t2)
    @test hcat(t, t2) isa StructArray
end

f_infer() = StructArray{ComplexF64}(rand(2,2), rand(2,2))

g_infer() = StructArray([(a=(b="1",), c=2)], unwrap = t -> t <: NamedTuple)
tup_infer() = StructArray([(1, 2), (3, 4)])

@testset "inferrability" begin
    @inferred f_infer()
    @inferred g_infer()
    @test g_infer().a.b == ["1"]
    s = @inferred tup_infer()
    @test s[1] == (1, 2)
    @test s[2] == (3, 4)
end

@testset "propertynames" begin
    a = StructArray{ComplexF64}(Float64[], Float64[])
    @test sort(collect(propertynames(a))) == [:im, :re]
end

struct S
    x::Int
    y::Float64
    S(x) = new(x, x)
end

LightStructArrays.SkipConstructor(::Type{<:S}) = true

@testset "inner" begin
    v = StructArray{S}([1], [1])
    @test v[1] == S(1)
    @test v[1].y isa Float64
end


text:
using GLVisualize, GeometryTypes, GLAbstraction, Colors

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0.0, 1.0, 360))
end
const N = 2048
function spiral(i, start_radius, offset)
	Point2f0(sin(i), cos(i)) * (start_radius + ((i/2pi)*offset))
end
# 2D particles
curve_data(i, N) = Point2f0[spiral(i+x/20f0, 1, (i/20)+1) for x=1:N]

t = const_lift(x-> (1f0-x)*100f0, timesignal)
color = map(RGBA{Float32}, colormap("Blues", N))
_view(visualize(const_lift(curve_data, t, N), :lines, color=color))

if !isdefined(:runtests)
    renderloop(window)
end


text:
#
#    Butterworth_filter_function_alt.m
#    ver 2.3   January 27, 2012
#    by Tom Irvine  Email: tomirvine@aol.com
#
#    Butterworth filter, sixth-order, infinite impulse response,
#    cascade with refiltering option for phase correction
#
#
#    iband:   1=lowpass  2=highpass  3=bandpass
#
#    iphase=1  for refiltering for phase correction
#          =2  for no refiltering
#
function Butterworth_filter_function_alt(y,dt,iband,fl,fh,iphase)
#
include("th_weighted_filter_coefficients.jl")
include("th_weighted_apply_filter.jl")
#
sr=1/dt;
n=length(y);
ns = n;
#
#
if(iband == 1)  # lowpass
          f=fl;
end
if(iband == 2)  # highpass
          f=fh;
end
#
iflag=1;
#
if(iband==1 || iband==2)
    if( f > 0.5*sr)
	    println("  error: cutoff frequency must be < Nyquist frequency ");
	    iflag=999;
    end
end
if(iband == 3)  # bandpass
          f=fh;
#
    if( fl > 0.5*sr || fh > 0.5*sr)
	    println("  error: cutoff frequency must be < Nyquist frequency ");
	    iflag=999;
    end
#
end
#
freq=f;
#
#****** calculate coefficients *******
#
if(iflag != 999 && iband !=3)
#
    (a,b,iflag) = th_weighted_filter_coefficients(f,dt,iband,iflag);
#
end
#
if(iflag < 900 )
#
    if(iband == 1 || iband ==2)  # lowpass or highpass
#
         if(iphase==1)   # refiltering
			y_temp=th_weighted_apply_filter(y,iphase,a,b);
			y_temp=th_weighted_apply_filter(y_temp,iphase,a,b);
	      else
			y_temp=th_weighted_apply_filter(y,iphase,a,b);
          end
#
    end
#
    if(iband == 3)  # bandpass
 #
		 f=fh;
		 freq=f;
#
##		 disp(" Step 1");
		 iband=2;
         (a,b,iflag) =th_weighted_filter_coefficients(f,dt,iband,iflag);
 #
         if(iphase==1) # refiltering
			y_temp=th_weighted_apply_filter(y,iphase,a,b);
			y_temp=th_weighted_apply_filter(y_temp,iphase,a,b);
	      else
			y_temp=th_weighted_apply_filter(y,iphase,a,b);
          end
#
		 f=fl;
		 freq=f;
#
##		 disp(" Step 2");
		 iband=1;
         (a,b,iflag) = th_weighted_filter_coefficients(f,dt,iband,iflag);
#
         if(iphase==1) # refiltering
			y_temp=th_weighted_apply_filter(y_temp,iphase,a,b);
			y_temp=th_weighted_apply_filter(y_temp,iphase,a,b);
	      else
			y_temp=th_weighted_apply_filter(y_temp,iphase,a,b);
          end
#
    end
#
    mu=mean(y_temp);
    sd=std(y_temp);
    mx=maximum(y_temp);
    mi=minimum(y_temp);
    rms=sqrt(sd^2+mu^2);
#    out1 = sprintf("\n Filtered Signal:  mean = #8.4g    std = #8.4g    rms = #8.4g ",mu,sd,rms);
#    disp(out1);
#
#
else
    println(" ")
    println("  Abnormal termination.  No output file generated. ");
end

return (y_temp,mu,sd,rms)
end


text:
const _ICD_ICD9_CCS_APPENDIXASINGLEDX_FILE_HAS_BEEN_PARSED = Ref(false)
const _ICD9_CODE_TO_SINGLE_LEVEL_DX_CCS = Ref(Dict{String, Int}())
const _SINGLE_LEVEL_DX_CCS_NUMBER_TO_NAME = Ref(Dict{Int, String}())
const _SINGLE_LEVEL_DX_CCS_TO_LIST_OF_ICD9_CODES = Ref(Dict{Int, Vector{String}}())

function remove_all_full_stops(x::AbstractString)::String
    result = replace(x, "." => "")
    return result
end

function parse_icd_icd9_ccs_appendixasingledx_file!()::Nothing
    global _ICD_ICD9_CCS_APPENDIXASINGLEDX_FILE_HAS_BEEN_PARSED
    global _ICD9_CODE_TO_SINGLE_LEVEL_DX_CCS
    global _SINGLE_LEVEL_DX_CCS_NUMBER_TO_NAME
    global _SINGLE_LEVEL_DX_CCS_TO_LIST_OF_ICD9_CODES
    if _ICD_ICD9_CCS_APPENDIXASINGLEDX_FILE_HAS_BEEN_PARSED[]
    else
        filename = package_directory(
            "assets",
            "icd",
            "icd9",
            "ccs",
            "AppendixASingleDX.txt"
            )
        file_sections = strip.(split(strip(read(filename, String)), "\n\n"))
        for section in file_sections
            if isempty(section)
            elseif startswith(section, "Appendix")
            elseif startswith(section, "Revised")
            else
                ccs_number::Int = parse(
                    Int,
                    strip.(
                        split(strip(section))
                        )[1],
                    )
                ccs_name::String = strip(
                    join(
                        strip.(
                            split(
                                strip(
                                    split(strip(section), "\n")[1]
                                    )
                                )[2:end]
                            ), " "
                        )
                    )
                _SINGLE_LEVEL_DX_CCS_NUMBER_TO_NAME[][ccs_number] = ccs_name
                icd9_code_list::Vector{String} = strip.(
                    remove_all_full_stops.(
                        split(
                            strip(
                                split(strip(section), "\n")[2]
                            )
                        )
                        )
                    )
                _SINGLE_LEVEL_DX_CCS_TO_LIST_OF_ICD9_CODES[][
                    ccs_number
                    ] = icd9_code_list
                for icd9_code in icd9_code_list
                    _ICD9_CODE_TO_SINGLE_LEVEL_DX_CCS[][icd9_code] = ccs_number
                end
            end
        end
        _ICD_ICD9_CCS_APPENDIXASINGLEDX_FILE_HAS_BEEN_PARSED[] = true
    end
    return nothing
end

function single_level_dx_ccs_number_to_name(
        ccs_number::Int,
        )::String
    parse_icd_icd9_ccs_appendixasingledx_file!()
    global _SINGLE_LEVEL_DX_CCS_NUMBER_TO_NAME
    result = _SINGLE_LEVEL_DX_CCS_NUMBER_TO_NAME[][ccs_number]
    return result
end

function single_level_dx_ccs_to_list_of_icd9_codes(
        ccs_number::Int,
        )::Vector{String}
    parse_icd_icd9_ccs_appendixasingledx_file!()
    global _SINGLE_LEVEL_DX_CCS_TO_LIST_OF_ICD9_CODES
    result = _SINGLE_LEVEL_DX_CCS_TO_LIST_OF_ICD9_CODES[][ccs_number]
    return result
end

"""
"""
function icd9_code_to_single_level_dx_ccs(
        icd9_code::AbstractString,
        )::Int
    parse_icd_icd9_ccs_appendixasingledx_file!()
    global _ICD9_CODE_TO_SINGLE_LEVEL_DX_CCS
    result = _ICD9_CODE_TO_SINGLE_LEVEL_DX_CCS[][
        remove_all_full_stops(string(icd9_code))
        ]
    return result
end


text:
immutable SparseMatrixBSC{Tv, Ti <: Integer}  <: AbstractBlockMatrix{Tv}
    R::Int                 # Block size in rows
    C::Int                 # Block size in columns
    m::Int                 # Number of rows
    n::Int                 # Number of columns
    colptr::Vector{Ti}     # Column i is in colptr[i]:(colptr[i+1]-1)
    rowval::Vector{Ti}     # Row values of blocks
    nzval::Array{Tv, 3}    # Nonzero values, one "matrix" per block, nzval[i, j, block]

    function SparseMatrixBSC(R::Integer, C::Integer,
                             m::Integer, n::Integer,
                             colptr::Vector{Ti}, rowval::Vector{Ti},  nzval::Array{Tv, 3})
      m < 0 && throw(ArgumentError("number of rows (m) must be ≥ 0, got $m"))
      n < 0 && throw(ArgumentError("number of columns (n) must be ≥ 0, got $n"))

      R < 0 && throw(ArgumentError("block size y must be ≥ 0, got $x_block_size"))
      C < 0 && throw(ArgumentError("block size x must be ≥ 0, got $y_block_size"))

      m % R != 0 && throw(ArgumentError("row block size: $(R) must evenly divide number of rows: $m"))
      n % C != 0 && throw(ArgumentError("column block size: $(C) must evenly divide number of rows: $n"))
      new(Int(R), Int(C), Int(m), Int(n), colptr, rowval, nzval)
  end
end

function SparseMatrixBSC{Tv, Ti <: Integer}(m::Integer, n::Integer, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Array{Tv, 3})
    R, C = size(nzval, 1), size(nzval, 2)
    SparseMatrixBSC{Tv, Ti}(R, C, m, n, colptr, rowval, nzval)
end

#####################
# Utility functions #
#####################

nblocks(A::SparseMatrixBSC) = (length(A.colptr) - 1, A.n ÷ A.C)
nblocks(A::SparseMatrixBSC, i::Int) = nblocks(A)[i]

blocksize(A::SparseMatrixBSC) = A.R, A.C
blocksize(A::SparseMatrixBSC, i::Int) = blocksize(A)[i]

nnzblocks(A::SparseMatrixBSC) = size(A.nzval, 1) * size(A.nzval, 2)
nzblockrange(A::SparseMatrixBSC, col::Integer) =  Int(A.colptr[col]):Int(A.colptr[col + 1] - 1)

Base.LinearIndexing(::Type{SparseMatrixBSC}) = Base.LinearSlow()
Base.size(A::SparseMatrixBSC) = (A.m, A.n)
Base.nnz(A::SparseMatrixBSC) = length(A.nzval)

function rowinblock(A, rowval, i)
    ((rowval - 1) * A.R + 1) <= i <= rowval * A.R
end

# Computes the blockindex which is the block I, J where the element with global index i, j
# would occupy
@inline blockindex(A::SparseMatrixBSC, i::Integer, j::Integer) = (i - 1) ÷ A.R + 1, (j - 1) ÷ A.C + 1

# Given a blockindex I, J and a global index i, j computes the offset α, β into the block such that
# A[i, j] = A.nzval[α, β, blockindex(I, J)]
@inline function blockoffsets(A::SparseMatrixBSC, I::Integer, J::Integer, i::Integer, j::Integer)
     α = i - A.R * (I - 1)
     β = j - A.C * (J - 1)
     return α, β
end

# Transforms from a global index i, j to a `BlockIndex`, see the BlockArrays package
function global2blockindex(A::SparseMatrixBSC, i::Integer, j::Integer)
    I, J = blockindex(A, i, j)
    α, β = blockoffsets(A, I, J, i, j)
    return BlockIndex((I, J), (α, β))
end

# Returns -1 if no block found else returns the block A[., ., block]
function findblock(block_array::SparseMatrixBSC, block_index::BlockIndex{2})
  nzrange = nzblockrange(A, block_index.I[2])
  # No block stored in this column
  if start(nzrange) > last(nzrange)
    return -1
  end

  block = searchsortedlast(A.rowval, block_index.I[1], start(nzrange), last(nzrange), Forward)
  if start(nzrange) > block
      return -1
  end

  if rowinblock(A, A.rowval[block], i)
      return block
  else
      return -1
  end
end

function Base.getindex(block_array::SparseMatrixBSC, block_index::BlockIndex{2})
  block = findblock(block_array, block_index)
  if block == -1
      return zero(T)
  else
    return A.nzval[block_index.α[1], block_index.α[2], block]
  end
end

#function Base.setindex!{T,N}(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N})
#    getblock(block_array, block_index.I...)[block_index.α...] = v
#end

function getblock{T}(A::SparseMatrixBSC{T}, I::Integer, J::Integer)
    @boundscheck checkbounds(A, I * blocksize(A, 1), J * blocksize(A, 2))
    getblock!(zeros(T, blocksize(A)), A, I, J)
end

function getblock!{T}(block::Matrix{T}, A::SparseMatrixBSC{T}, I::Integer, J::Integer)
    @boundscheck checkbounds(A, I * blocksize(A, 1), J * blocksize(A, 2))
    @boundscheck @assert size(block) == blocksize(A, 1)

    nzrange = nzblockrange(A, J)

    if start(nzrange) > last(nzrange)
       return fill!(block, 0.0)
    end
    block = searchsortedlast(A.rowval, I, start(nzrange), last(nzrange), Forward)

    if (start(nzrange) > block) || (A.rowval[block] != I)
         return fill!(block, 0.0)
    else
        @inbounds for i in blocksize(A, 1), j in blocksize(A, 2)
            block[i, j] = A.nzval[i, j, block]
        end
        return block
    end
end

############
# Printing #
############
function Base.showarray(io::IO, A::SparseMatrixBSC;
                   header::Bool=true, repr=false)
    print(io, A.m, "×", A.n, " sparse block matrix with ", nnzblocks(A), " ", eltype(A), " block entries of size ",
        blocksize(A, 1),"×",blocksize(A, 2))
end

function Base.getindex{T}(A::SparseMatrixBSC{T}, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    J = (j - 1) ÷ A.C + 1
    I = (i - 1) ÷ A.R + 1
    nzrange = nzblockrange(A, J)
    # No block stored in this column
    if start(nzrange) > last(nzrange)
      return zero(T)
    end

    block = searchsortedlast(A.rowval, I, start(nzrange), last(nzrange), Forward)
    if start(nzrange) > block
        return zero(T)
    end

    @inbounds if rowinblock(A, A.rowval[block], i)
        i_blk, j_blk = blockoffsets(A, I, J, i, j)
        return A.nzval[i_blk, j_blk, block]
    else
        return zero(T)
    end
end

###########################
# Conversions CSC <-> BSC #
###########################

# Strategy:
function SparseMatrixCSC{Tv, Ti <: Integer}(A::SparseMatrixBSC{Tv, Ti})
    if blocksize(A) == (1,1)
        return SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, vec(A.nzval))
    end
    rowval = zeros(Ti, nnz(A))
    colptr = zeros(Ti, size(A, 2) + 1)
    nzval = zeros(Tv, length(A.nzval))

    count_row = 1
    count_col = 2
    colptr[1] = 1
    @inbounds for col in 1:nblocks(A, 2)
        blockrange = nzblockrange(A, col)
        n_blocks_col = length(blockrange)
        nnz_values_col = n_blocks_col * blocksize(A, 1)
        for j_blk in 1:blocksize(A, 2)
            # The new colptr is the previous one plus the number of nonzero elements in this column.
            colptr[count_col] = colptr[count_col - 1] + nnz_values_col
            count_col += 1
            for block in blockrange
                i_offset = (A.rowval[block] - 1) * blocksize(A, 1)
                for i_blk in 1:blocksize(A, 1)
                    nzval[count_row] = A.nzval[i_blk, j_blk, block]
                    rowval[count_row] = i_blk + i_offset
                    count_row += 1
                end
            end
        end
    end
    return SparseMatrixCSC(A.m, A.n, colptr, rowval, nzval)
end

function SparseMatrixBSC{Tv, Ti <: Integer}(A::SparseMatrixCSC{Tv, Ti}, R::Integer, C::Integer)
    if (R, C) == (1,1)
        return SparseMatrixBSC(1, 1, A.m, A.n, A.colptr, A.rowval, reshape(A.nzval, length(A.nzval), 1, 1))
    end

    A.m % R != 0 && throw(ArgumentError("row block size: $(R) must evenly divide number of rows: $(A.m)"))
    A.n % C != 0 && throw(ArgumentError("column block size: $(C) must evenly divide number of rows: $(A.n)"))

    Anzval = A.nzval
    Arowval = A.rowval
    Acolptr = A.colptr

    # Upper bound of number of nonzero blocks is nnz(A).
    rowval = zeros(Ti, nnz(A))
    colptr = zeros(Ti, A.n ÷ C + 1)

    n_colblocks = div(A.n, C)

    colptr[1] = 1
    row_counter = 1
    rows = Int[]

    # Strategy to compute rowval and colptr:
    # For each column block accumulate all the rowvalues in the CSC matrix for that block.
    # Convert these to what rowblock they represent
    # Each unique rowblock should now be entered in order into rowval.
    # Colptr for this column block is incremented by the number of unique rowblock values.
    @inbounds for colblock in 1:n_colblocks
        j_offset = (colblock - 1) * C
        row_block_counter = 1

        # Count the number of non zero values for the columns in this column block
        nzvals_block = Acolptr[j_offset + C + 1] - Acolptr[j_offset + 1]
        if nzvals_block == 0 # No nz values in this column block, exit early
            colptr[colblock + 1] = row_counter
            continue
        end

        # Accumulate rowvals for this block
        resize!(rows, nzvals_block)
        for j_blk in 1:C
            col = j_offset + j_blk
            nz_range = Acolptr[col]:Acolptr[col + 1] - 1
            for r in nz_range
                rows[row_block_counter] = Arowval[r]
                row_block_counter += 1
            end
        end

        # Convert from row values -> block rows
        @simd for i in 1:length(rows)
            rows[i] = (rows[i] - 1) ÷ R + 1
        end

        # Pick out the unique block rows and put them into rowval
        sort!(rows) # <- A bit of a bottle enck, takes about 30% of tot time.
        rowval[row_counter] = rows[1] # We have at least one value in rows so this is ok
        row_counter += 1
        for i in 2:length(rows)
            if rows[i] > rows[i-1]
                rowval[row_counter] = rows[i]
                row_counter += 1
            end
        end
        colptr[colblock + 1] = row_counter
    end

    # We now know the true number of non zero blocks so we reshape rowval
    # and allocate the exact space we need for nzval
    deleteat!(rowval, row_counter:length(rowval))
    nzval = zeros(Tv, R, C, length(rowval))

    @inbounds for colblock in 1:n_colblocks
        j_offset = (colblock - 1) * C
        for j_blk in 1:C
            current_block = colptr[colblock]
            col = j_offset + j_blk
            for r in Acolptr[col]:Acolptr[col + 1] - 1
                row = Arowval[r]
                # Looking for the correct block for this column
                while row > rowval[current_block] * R
                    current_block += 1
                end
                i_blk = row - (rowval[current_block] - 1) * R
                nzval[i_blk, j_blk, current_block] = Anzval[r]
            end
        end
    end

    SparseMatrixBSC(A.m, A.n, colptr, rowval, nzval)
end

##########
# LinAlg #
##########

# We do dynamic dispatch here so that the size of the blocks are known at compile time
function A_mul_B!{Tv}(b::Vector{Tv}, A::SparseMatrixBSC{Tv}, x::Vector{Tv})
    if blocksize(A) == (1,1)
        A_mul_B!(b, SparseMatrixCSC(A), x)
    else
        _A_mul_B!(b, A, x, Val{A.C}, Val{A.R})
    end
end

function _A_mul_B!{Tv, C, R}(b::Vector{Tv}, A::SparseMatrixBSC{Tv}, x::Vector{Tv},
                             ::Type{Val{C}}, ::Type{Val{R}})
    fill!(b, 0.0)
    n_cb, n_rb = nblocks(A)
    for J in 1:n_cb
        j_offset = (J - 1)  * blocksize(A, 2)
        for r in nzblockrange(A, J)
            @inbounds i_offset = (A.rowval[r] - 1) * blocksize(A, 1)
            matvec_kernel!(b, A, x, r, i_offset, j_offset, Val{C}, Val{R})
        end
    end
    return b
end

# TODO: Possibly SIMD.jl and do a bunch of simd magic coolness for small mat * vec
@inline function matvec_kernel!{C, R, T}(b::Vector{T}, A::SparseMatrixBSC{T}, x::Vector{T}, r,
                                              i_offset, j_offset, ::Type{Val{C}}, ::Type{Val{R}})
    @inbounds for j in 1:R
        for i in 1:C
            b[i_offset + i] += A.nzval[i, j, r] * x[j_offset + j]
        end
    end
end

function Base.:*{Tv}(A::SparseMatrixBSC{Tv}, x::Vector{Tv})
    A_mul_B!(similar(x, A.m), A, x)
end


text:
using Pkg
Pkg.instantiate()

using JuliaFormatter
using Dates
using JSON

format_text("")

const logfile = open(joinpath(@__DIR__, "juliaformatter.log"), "w")

function log(msg; spacer = " ")
    write(logfile, "[$(Dates.now())]$spacer$msg\n")
    flush(logfile)
end

const CONFIG_FILE_NAME = JuliaFormatter.CONFIG_FILE_NAME

function find_config_file(dir)
    dir2config = Dict{String,Any}()
    next_dir = dirname(dir)
    config = if (next_dir == dir || # ensure to escape infinite recursion
                 isempty(dir)) # reached to the system root
        nothing
    elseif haskey(dir2config, dir)
        dir2config[dir]
    else
        path = joinpath(dir, CONFIG_FILE_NAME)
        isfile(path) ? JuliaFormatter.parse_config(path) : find_config_file(next_dir)
    end
    return dir2config[dir] = config
end

function main()
    server_state = "start"
    while server_state != "quit"
        text = readline(stdin)
        data = String(text)
        format_options = Dict{Symbol,Any}()
        if length(data) == 0
            continue
        end
        try
            data = JSON.parse(String(data))
        catch e
            iob = IOBuffer()
            showerror(iob, e, catch_backtrace())
            log("Unable to parse json: \"$data\". $(String(take!(iob)))")
            continue
        end
        if data["method"] == "exit"
            server_state = "quit"
        elseif data["method"] == "isconnectedcheck"
            log("Connected.")
        elseif data["method"] == "format"
            log("Setting up defaults ...")
            text = data["params"]["text"]
            options = data["params"]["options"]
            filepath = data["params"]["filepath"]
            style = pop!(options, "style", nothing)
            style = if style == "blue"
                BlueStyle()
            elseif style == "yas"
                YASStyle()
            elseif style == nothing || style == "default"
                DefaultStyle()
            else
                log("Unknown style option $style")
                DefaultStyle()
            end
            for (k, v) in options
                format_options[Symbol(k)] = v
            end
            format_options[:style] = style
            log("Searching for .JuliaFormatter.toml ...")
            dir = dirname(filepath)
            format_options = if (config = find_config_file(dir)) !== nothing
                log("Found .JuliaFormatter.toml with config: $config")
                merge(format_options, config)
            else
                log("No .JuliaFormatter.toml")
                format_options
            end
            log("Using options: $format_options with style: $(format_options[:style])")
            output = text
            indent = typemax(Int64)
            for line in text
                if length(line) > 0
                    indent = min(length(line) - length(lstrip(line)), indent)
                end
            end
            if indent == typemax(Int64)
                indent = 0
            end
            log("Formatting: ")
            log(join(text, "\n"), spacer = "\n")
            try
                output = format_text(join(text, "\n"); format_options...)
                data["status"] = "success"
                log("Success")
            catch e
                output = join(text, "\n")
                data["status"] = "error"
                iob = IOBuffer()
                showerror(iob, e, catch_backtrace())
                log(String(take!(iob)))
            end
            log("\n---------------------------------------------------------------------\n")
            log("Formatted: ")
            log(output, spacer = "\n")
            data["params"]["text"] = [rstrip(lpad(l, length(l) + indent)) for l in split(output, "\n")]
            println(stdout, JSON.json(data))
            log("Done.")
        end
    end

    log("exiting ...")
end

log("calling main ...")

try
    main()
catch e
    iob = IOBuffer()
    showerror(iob, e, catch_backtrace())
    log(String(take!(iob)))
end


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

"""
    Plane(p₀, v, w)

A plane coincident with point `p₀`, defined by non-parallel vectors `v` and `w`.
It can be called as p(s, t) with numeric parameters `s` and `t` to cast it at
`p₀ + s*v + t*w`.
"""
struct Plane{Dim,T} <: Primitive{Dim,T}
  pₒ::Point{Dim,T}
  v::Vec{Dim,T}
  w::Vec{Dim,T}
end

Plane(pₒ::Tuple, v::Tuple, w::Tuple) = Plane(Point(pₒ), Vec(v), Vec(w))

function (p::Plane)(s, t)
  p.pₒ + s*p.v + t*p.w
end

paramdim(::Type{<:Plane}) = 2


text:
using LocalMemoize

abstract type Light end

const Lights = Vector{Light}

struct LightSample
    light_vector::Vector3
    irradiance::ColorRGB
end
zero(::Type{LightSample}) = LightSample(zero(Vector3), Black)
InShadow() = zero(LightSample)
in_shadow(ls::LightSample) = ls.light_vector == zero(Vector3)

function shadow_test(𝐥::Vector3, position::Vector3, scene::Geometry, r::Float64=∞)
    shadow_ray = Ray3(position, 𝐥)
    shadow_result = shadow_ray ∩ scene
    hit(shadow_result) && shadow_result.value.distance <= r
end

struct DirectionalLight <: Light
    irradiance::ColorRGB
    direction::Vector3
    shadow::Bool
end
DirectionalLight(irradiance, direction) = DirectionalLight(irradiance, direction, true)
@memoize get_𝐥(dl::DirectionalLight) = -normalize(dl.direction)

function sample(@memo(dl::DirectionalLight), scene::Geometry, position::Vector3)
    𝐥 = get_𝐥(dl)
    if dl.shadow && shadow_test(𝐥, position, scene)
        return InShadow()
    end
    LightSample(𝐥, dl.irradiance)
end

struct PointLight <: Light
    intensity::ColorRGB
    position::Vector3
    shadow::Bool
end
PointLight(intensity, position) = PointLight(intensity, position, true)

function sample(pl::PointLight, scene::Geometry, position::Vector3)
    Δr = pl.position - position
    r² = norm²(Δr)
    r = √r²
    𝐥 = Δr / r
    if pl.shadow && shadow_test(𝐥, position, scene, r)
        return InShadow()
    end
    LightSample(𝐥, pl.intensity / r²)
end

struct SpotLight <: Light
    intensity::ColorRGB
    position::Vector3
    direction::Vector3
    θ::Float64
    ϕ::Float64
    falloff::Float64
    shadow::Bool
end
SpotLight(intensity, position, direction, θ, ϕ, falloff) = SpotLight(
    intensity, position, direction, θ, ϕ, falloff, true)

@memoize function sample(@memo(sl::SpotLight), scene::Geometry, position::Vector3)
    @memo 𝐬 = -normalize(sl.direction)
    @memo cosθ = cos(sl.θ / 2)
    @memo cosϕ = cos(sl.ϕ / 2)
    @memo base_multiplier = 1 / (cosθ - cosϕ)
    Δr = sl.position - position
    r² = norm²(Δr)
    r = √r²
    𝐥 = Δr / r
    if sl.shadow && shadow_test(𝐥, position, scene, r)
        return InShadow()
    end
    sdotl = 𝐬 ⋅ 𝐥
    if sdotl >= cosθ
        spot = 1.0
    elseif sdotl <= cosϕ
        spot = 0.0
    else
        spot = ((sdotl - cosϕ) * base_multiplier) ^ sl.falloff
    end
    LightSample(𝐥, sl.intensity * spot / r²)
end


text:
module WhiskerStructs
struct HydroWhisker
    λ::Float64
    γ::Float64
    Ac::Float64
    At::Float64
    ϕ::Float64
    ϵ::Float64
    T::Float64
end

struct TopoWhisker
    a::Float64
    b::Float64
    k::Float64
    l::Float64
    α::Float64
    β::Float64
    M::Float64
end

export HydroWhisker, TopoWhisker

end


text:
using LightXML

function iwslt2016_init()

    # Helper function for extracting nested archives
    function extract_lang_archive(src, dst, fn)
        # Unpack outer archive and construct path to the archive for the language pair
        unpack(fn)
        archivename = "2016-01"
        innerdir = "$(src)-$(dst)"
        langarchive = joinpath(archivename, "texts", "$(src)", "$(dst)", innerdir * ".tgz")

        # Unpack language pair archive
        unpack(langarchive)
        innerfiles = readdir(innerdir)
        mv.(joinpath.(innerdir, innerfiles), innerfiles)

        for f ∈ innerfiles
             if occursin(".xml", f)
                 clean_xml(f)
             elseif occursin(".tags", f)
                 clean_tag(f)
             end
        end
        rm(innerdir)
        rm(archivename; recursive=true)
    end

    archivehash = "425a3688e0faff00ed4d6d04f1664d1edbd9932e5b17a73680aa81a70f03e2d6"

    message = (src, dst) -> """
    The IWSLT 2016 TED talk translation task

    These are the data sets for the MT tasks of the evaluation campaigns of IWSLT. They are parallel data sets used for building and testing MT systems. They are publicly available through the WIT3 website wit3.fbk.eu, see release: 2016-01.

    Data are crawled from the TED website and carry the respective licensing conditions (for training, tuning and testing MT systems).
    Approximately, for each language pair, training sets include 2,000 talks, 200K sentences and 4M tokens per side, while each dev and test sets 10-15 talks, 1.0K-1.5K sentences and 20K-30K tokens per side. In each edition, the training sets of previous editions are re-used and updated with new talks added to the TED repository in the meanwhile.

    from $(src) to $(dst)
    """

    for lang ∈ ("ar", "cs", "fr", "de")
        register(DataDep(
            "IWSLT2016 $(lang)-en",
            message(lang, "en"),
            "https://drive.google.com/file/d/1l5y6Giag9aRPwGtuZHswh3w5v3qEz8D8/",
            archivehash;
            fetch_method=gdownload,
            post_fetch_method = fn -> extract_lang_archive(lang, "en", fn)
        ))

        register(DataDep(
            "IWSLT2016 en-$(lang)",
            message("en", lang),
            "https://drive.google.com/file/d/1l5y6Giag9aRPwGtuZHswh3w5v3qEz8D8/",
            archivehash;
            fetch_method=gdownload,
            post_fetch_method = fn -> extract_lang_archive("en", lang, fn)
        ))
    end
end

function clean_xml(f)
    nf = replace(f, "xml"=>"txt")
    xdoc = parse_file(f)
    xroot = root(xdoc)
    node = length(xroot["srcset"]) == 0 ? xroot["refset"][1] : xroot["srcset"][1]
    open(nf, "w+") do fw
        for doc in node["doc"]
            for seg in doc["seg"]
                println(fw, content(seg))
            end
        end
    end
    free(xdoc)
    nf
end

function clean_tag(f)
    nf = replace(f, ".tags"=>"")
    open(f) do fr
        open(nf, "w+") do fw
            for l ∈ eachline(fr)
                if !startswith(l, "<")
                    println(fw, l)
                end
            end
        end
    end
    nf
end

struct IWSLT2016 <: Dataset
    src
    ref
    function IWSLT2016(src, ref)
        provided = [:cs, :ar, :fr, :de]
        !(src == :en ? ref ∈ provided : ref == :en && src ∈ provided) && error("language not provided")
        new(src, ref)
    end
end

function trainfile(iw::IWSLT2016)
    p = "$(iw.src)-$(iw.ref)"
    src = @datadep_str "IWSLT2016 $p/train.$(p).$(iw.src)"
    ref = @datadep_str "IWSLT2016 $p/train.$(p).$(iw.ref)"
    src, ref
end

function tunefile(iw::IWSLT2016, dev, year; tedx = false)
    p = "$(iw.src)-$(iw.ref)"
    hX = tedx ? "X" : ""
    catl = dev ? "dev" : "tst"
    srcf = "IWSLT16.TED$(hX).$(catl)$(year).$(p).$(iw.src).txt"
    reff = "IWSLT16.TED$(hX).$(catl)$(year).$(p).$(iw.ref).txt"

    !(srcf ∈ readdir(@datadep_str "IWSLT2016 $p")) && error("""no such file: $srcf,
                                                             only have the following:
                                                             $(join(filter(x->occursin(".$(iw.src).txt", x), readdir(datadep"IWSLT2016 $p/")), "\n"))""")

    src = @datadep_str "IWSLT2016 $p/$srcf"
    ref = @datadep_str "IWSLT2016 $p/$reff"
    src, ref
end

devfile(iw::IWSLT2016, year=2010; tedx = false) = tunefile(iw, true, year; tedx=tedx)
testfile(iw::IWSLT2016, year=2010; tedx = false) = tunefile(iw, false, year; tedx=tedx)

function get_vocab(iw::IWSLT2016; mode = :vocab, min_freq = 7)
    if mode == :vocab
        sf, rf = trainfile(iw)
        vocab = token_freq(sf, rf; min_freq = min_freq)

        return vocab
    else
        error("mode should be one of :vocab")
    end
end


text:
module TestLeafVararg

using ComoniconTypes
using ComoniconTargetExpr
using ComoniconTargetExpr: emit_expr, emit_body, emit_norm_body, emit_dash_body
using Test

const test_args = Ref{Vector{Any}}()

function foo(a, b=2, c...)
    test_args[] = [a, b, c...]
end

cmd = Entry(;
    version=v"1.1.0",
    root=LeafCommand(;
        fn=foo,
        name="leaf",
        args=[
            Argument(;name="a", type=Int),
            Argument(;name="b", type=Int, require=false),
        ],
        vararg=Argument(;name="c", type=Float64, vararg=true),
    )
)

eval(emit_expr(cmd))

@testset "test leaf optional argument" begin
    @test command_main(["3", "2", "5.5", "6", "7"]) == 0
    @test test_args[] == Any[3, 2, 5.5, 6.0, 7.0]
    @test command_main(["3", "2"]) == 0
    @test test_args[] == Any[3, 2]
end

end


text:
using HackerNews, Compat
using Compat.Test

# testing maxitems
maxitem = HackerNews.HN(HackerNews.MaxItem)
result_maxitem = HackerNews.getinfo(maxitem)

@test hn_route_string(maxitem.story) == "maxitem"
@test maxitem.nposts == 1
@test typeof(result_maxitem[1].data) == Dict{String,Any}

# testing jobstories
jobstories = HackerNews.HN(HackerNews.JobStories,2)
result_jobstories = HackerNews.getinfo(jobstories)

@test hn_route_string(jobstories.story) == "jobstories"
@test jobstories.nposts == 2
@test typeof(result_jobstories[1].data) == Dict{String,Any}

# testing updates

## user updates
user_updates = HackerNews.HN(HackerNews.Updates, 1, true)
result_user_updates = HackerNews.getinfo(user_updates)

@test hn_route_string(user_updates.story) == "updates"
@test user_updates.nposts == 1
@test typeof(result_user_updates) == Array{HackerNews.HNUser,1}
@test typeof(result_user_updates[1].data) == Dict{String,Any}

## post/item updates
post_updates = HackerNews.HN(HackerNews.Updates, 1, false)
result_post_updates = HackerNews.getinfo(post_updates)

@test hn_route_string(post_updates.story) == "updates"
@test post_updates.nposts == 1
@test typeof(result_post_updates) == Array{HackerNews.HNPost,1}
@test typeof(result_post_updates[1].data) == Dict{String,Any}

# testing user information
user = HackerNews.getuser("pg")

@test typeof(user) == HackerNews.HNUser
@test typeof(user.data) == Dict{String,Any}
@test user.id == "pg"


text:
struct SupervisedGridSearch
    gamma_range
    C_range
end

function SupervisedGridSearch(;gamma_range, C_range)
    SupervisedGridSearch(gamma_range, C_range)
end

function Base.string(m::SupervisedGridSearch)
    return "SupervisedGridSearch_C-$(minimum(m.C_range))-$(maximum(m.C_range))_gamma-$(minimum(m.gamma_range))-$(maximum(m.gamma_range))"
end

function find_parameters(method::SupervisedGridSearch, data, labels, oracle, solver, qm, seed)
    model = VanillaSVDD(data)
    best_gamma = 1.0
    best_C = 1.0
    best_score = -Inf
    for γ in method.gamma_range
        SVDD.initialize!(model, SimpleCombinedStrategy(FixedGammaStrategy(MLKernels.GaussianKernel(γ)), FixedCStrategy(1.0)))
        set_adjust_K!(model, true)
        for C in method.C_range
            set_param!(model, Dict(:C => C))
            SVDD.fit!(model, solver)
            current_score = qm(ConfusionMatrix(SVDD.classify.(SVDD.predict(model, data)), labels))

            if current_score > best_score
                best_C = C
                best_gamma = γ
                best_score = current_score
            end
        end
    end
    return (best_gamma, best_C)
end

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama


text:
@testset "1491.average-salary-excluding-the-minimum-and-maximum-salary.jl" begin
    @test average_salary_excluding_mami([4000, 3000, 1000, 2000]) ≈ 2500.00000
    @test average_salary_excluding_mami([1000, 2000, 3000]) ≈ 2000.00000
    @test average_salary_excluding_mami([6000, 5000, 4000, 3000, 2000, 1000]) ≈ 3500.00000
    @test average_salary_excluding_mami([8000, 9000, 2000, 3000, 6000, 1000]) ≈ 4750.00000
end


text:
module LibCImPlot

using CImGui.LibCImGui.CImGui_jll
using CImGui.LibCImGui

using CEnum
# using CImGui: ImVec2, ImVec4, ImGuiMouseButton, ImGuiKeyModFlags
# import CImGui.LibCImGui:ImGuiContext

include(joinpath(@__DIR__, "..", "..", "gen", "libcimplot_common.jl"))
include(joinpath(@__DIR__, "..", "..", "gen", "libcimplot_api.jl"))

# for i in [ImPlotFlags_, ImPlotAxisFlags_]
#     for j in instances(i)
#         @eval export $(Symbol(j))
#     end
# end

# this should export all flags etc
foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), "ImPlot")
        @eval export $s
    end
end

export CreateContext, DestroyContext, SetImGuiContext

end # module


text:
using MotionCaptureJointCalibration
using DrakeVisualizer
using RigidBodyTreeInspector
using Interact

using MotionCaptureJointCalibration.SyntheticDataGeneration

using RigidBodyDynamics
using StaticArrays
using ValkyrieRobot
using ForwardDiff
using Ipopt
using NBInclude
using Test
using LinearAlgebra

using MotionCaptureJointCalibration: Point3DS, reconstruct!, deconstruct, _marker_residual, _∇marker_residual!
using Random: seed!, rand!

T = Float64
val = Valkyrie()
mechanism = val.mechanism
remove_fixed_tree_joints!(mechanism)
state = MechanismState{T}(mechanism)

markerbodies = findbody.(Ref(mechanism), ["leftFoot", "pelvis"])
seed!(1)
body_weights = Dict(b => rand() for b in markerbodies)
marker_options = MarkerPositionGenerationOptions()
pose_options = PoseDataGenerationOptions()
problem, groundtruth = generate_calibration_problem(state, body_weights, marker_options = marker_options, pose_options = pose_options)

@testset "deconstruct/reconstruct!" begin
    q = zeros(num_positions(mechanism))
    marker_positions_body = Dict(b => [Point3D(default_frame(b), 0., 0., 0.) for i = 1 : num_markers(problem, b)] for b in markerbodies)
    reconstruct!(markerbodies, q, marker_positions_body, deconstruct(markerbodies, configuration(state), groundtruth.marker_positions)...)
    @test all(q .== configuration(state))
    for (body, positions) in groundtruth.marker_positions
        @test all(marker_positions_body[body] .== positions)
    end
end

@testset "marker_residual gradient" begin
    data = problem.pose_data[1]
    marker_positions_world = data.marker_positions

    function marker_residual_inefficient(x::X...) where {X}
        M = eltype(mechanism)
        state = MechanismState{X}(mechanism)
        marker_positions_body = Dict{RigidBody{M}, Vector{Point3DS{X}}}(
            b => [Point3D(default_frame(b), zero(X), zero(X), zero(X)) for i = 1 : num_markers(problem, b)] for b in markerbodies
        )
        reconstruct!(markerbodies, configuration(state), marker_positions_body, x...)
        normalize_configuration!(state)
        setdirty!(state)
        _marker_residual(state, markerbodies, marker_positions_world, marker_positions_body, body_weights)
    end

    f(args) = [marker_residual_inefficient(args...)]

    for i = 1 : 100
        rand!(state)
        q = configuration(state)
        Jcheck = ForwardDiff.jacobian(f, deconstruct(markerbodies, q, groundtruth.marker_positions))

        set_configuration!(state, q)
        g = zeros(length(Jcheck))
        paths_to_root = Dict(b => RigidBodyDynamics.path(mechanism, root_body(mechanism), b) for b in markerbodies)
        jacobians = Dict(b => (p => geometric_jacobian(state, p)) for (b, p) in paths_to_root)
        _∇marker_residual!(g, state, markerbodies, marker_positions_world, groundtruth.marker_positions, body_weights, jacobians)
        J = g'

        @test isapprox(J, Jcheck, atol = 1e-14)
    end
end

@testset "problem" begin
    @test num_poses(problem) == pose_options.num_poses
    @test num_calibration_params(problem) == 6
    @test num_markers(problem) == marker_options.num_markers * length(markerbodies)
    @test num_bodies(problem) == length(markerbodies)
    show(devnull, problem)
end

@testset "solve" begin
    # NLopt SLSQP works well with up to 10 poses, free floating joint configurations and two unmeasured markers
    # solver = NLoptSolver(algorithm = :LD_SLSQP)

    solver = IpoptSolver(print_level = 4, max_iter = 10000, derivative_test = "first-order", check_derivatives_for_naninf = "yes", tol = 1e-10)
    # other useful options:
    # hessian_approximation = "limited-memory"

    result = solve(problem, solver)
    @test result.status == :Optimal

    @test num_poses(problem) == num_poses(result)
    @test num_calibration_params(problem) == num_calibration_params(result)
    @test num_markers(problem) == num_markers(result)
    @test num_bodies(problem) == num_bodies(result)

    # check calibration parameters
    for joint in calibration_joints(problem)
        @test isapprox(result.calibration_params[joint], groundtruth.calibration_params[joint]; atol = 1e-3)
    end

    # check configurations
    solutionstate = MechanismState{T}(mechanism)
    groundtruthstate = MechanismState{T}(mechanism)
    for i = 1 : num_poses(problem)
        set_configuration!(solutionstate, result.configurations[i])
        set_configuration!(groundtruthstate, groundtruth.configurations[i])
        for body in bodies(mechanism)
            @test isapprox(transform_to_root(solutionstate, body), transform_to_root(groundtruthstate, body); atol = 2e-3)
        end
    end

    # printing and visualization (just to make sure the code doesn't error)
    show(devnull, result)

    vis = Visualizer()[:valkyrie]
    geometry = visual_elements(mechanism, URDFVisuals(ValkyrieRobot.urdfpath(); package_path = [ValkyrieRobot.packagepath()]))
    setgeometry!(vis, mechanism, geometry)
    inspect!(state, vis, problem, result)
end

using RigidBodyTreeInspector

@testset "example notebooks" begin
    notebookdir = joinpath(@__DIR__, "..", "notebooks")
    excludedirs = [".ipynb_checkpoints"]
    excludefiles = String[]
    for (root, dir, files) in walkdir(notebookdir)
        basename(root) in excludedirs && continue
        for file in files
            file in excludefiles && continue
            name, ext = splitext(file)
            lowercase(ext) == ".ipynb" || continue
            path = joinpath(root, file)
            @eval module $(gensym()) # Each notebook is run in its own module.
            using Test
            using NBInclude
            @testset "Notebook: $($name)" begin
                # Note: use #NBSKIP in a cell to skip it during tests.
                @nbinclude($path; regex = r"^((?!\#NBSKIP).)*$"s)
            end
            end # module
        end
    end
end