text:
"""
**QuadratureVI**

Variational Inference solver by approximating gradients via numerical integration via Quadrature

```julia
QuadratureVI(ϵ::T=1e-5,nGaussHermite::Integer=20,optimizer::Optimizer=Momentum(η=0.0001))
```

**Keyword arguments**

    - `ϵ::T` : convergence criteria
    - `nGaussHermite::Int` : Number of points for the integral estimation
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl) package. Default is `Momentum(η=0.0001)`
"""
mutable struct QuadratureVI{T<:Real,N} <: NumericalVI{T}
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    nPoints::Int64 #Number of points for the quadrature
    nodes::Vector{T}
    weights::Vector{T}
    clipping::T
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64 #Number of samples of the data
    nMinibatch::Int64 #Size of mini-batches
    ρ::T #Stochastic Coefficient
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    vi_opt::NTuple{N,NVIOptimizer}
    MBIndices::Vector #Indices of the minibatch
    xview::SubArray{T,2,Matrix{T}}
    yview::SubArray

    function QuadratureVI{T}(ϵ::T,nPoints::Integer,optimizer::Optimizer,Stochastic::Bool,clipping::Real,nMinibatch::Int) where {T}
        return new{T,1}(ϵ,0,nPoints,[],[],clipping,Stochastic,0,nMinibatch,1.0,true,(NVIOptimizer{T}(0,0,optimizer),))
    end

    function QuadratureVI{T,1}(ϵ::T,Stochastic::Bool,nPoints::Int,clipping::Real,nFeatures::Int,nSamples::Int,nMinibatch::Int,nLatent::Int,optimizer::Optimizer) where {T}
        gh = gausshermite(nPoints)
        vi_opts = ntuple(_->NVIOptimizer{T}(nFeatures,nMinibatch,optimizer),nLatent)
        new{T,nLatent}(ϵ,0,nPoints,gh[1].*sqrt2,gh[2]./sqrtπ,clipping,Stochastic,nSamples,nMinibatch,nSamples/nMinibatch,true,vi_opts,collect(1:nMinibatch))
    end
end

function QuadratureVI(;ϵ::T=1e-5,nGaussHermite::Integer=100,optimizer::Optimizer=Momentum(η=1e-5),clipping::Real=0.0) where {T<:Real}
    QuadratureVI{T}(ϵ,nGaussHermite,optimizer,false,clipping,1)
end

"""
**QuadratureSVI**

Stochastic Variational Inference solver by approximating gradients via numerical integration via Quadrature

```julia
QuadratureSVI(nMinibatch::Integer;ϵ::T=1e-5,nGaussHermite::Integer=20,optimizer::Optimizer=Adam(α=0.1))
```
    -`nMinibatch::Integer` : Number of samples per mini-batches

**Keyword arguments**

    - `ϵ::T` : convergence criteria, which can be user defined
    - `nGaussHermite::Int` : Number of points for the integral estimation (for the QuadratureVI)
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl) package. Default is `Momentum(η=0.001)`
"""
function QuadratureSVI(nMinibatch::Integer;ϵ::T=1e-5,nGaussHermite::Integer=100,optimizer::Optimizer=Momentum(η=1e-5),clipping::Real=0.0) where {T<:Real}
    QuadratureVI{T}(ϵ,nGaussHermite,optimizer,clipping,nMinibatch)
end

function tuple_inference(i::TInf,nLatent::Integer,nFeatures::Integer,nSamples::Integer,nMinibatch::Integer) where {TInf <: QuadratureVI}
    return TInf(i.ϵ,i.Stochastic,i.nPoints,i.clipping,nFeatures,nSamples,nMinibatch,nLatent,i.vi_opt[1].optimizer)
end

function expec_log_likelihood(model::VGP{T,L,<:QuadratureVI}) where {T,L}
    tot = 0.0
    for gp in model.f
        for i in 1:model.nSamples
            x = model.inference.nodes*sqrt(max(gp.Σ[i,i],0.0)) .+ gp.μ[i]
            tot += dot(model.inference.weights,logpdf.(model.likelihood,model.y[i],x))
        end
    end
    return tot
end

function expec_log_likelihood(l::Likelihood,i::QuadratureVI,y::AbstractVector,μ::AbstractVector,Σ::AbstractMatrix)

    # for j in 1:length(y)
        # nodes = i.nodes*sqrt(Σ[j,j]) .+ μ[j]
        # tot += dot(i.weights,logpdf.(l,y[j],nodes))
    # end
    # return tot
    @show "TEST5"
    sum(apply_quad.(y,μ,diag(Σ),i,l))
end

function apply_quad(y::Real,μ::Real,σ²::Real,i::QuadratureVI,l::Likelihood) where {T}
    x = i.nodes*sqrt(max(σ², zero(σ²))) .+ μ
    return dot(i.weights,x)
    # return dot(i.weights,logpdf.(l,y,x))
end

function expec_log_likelihood(model::SVGP{T,L,<:QuadratureVI}) where {T,L}
    tot = 0.0
    y = get_y(model)
    for gp in model.f
        μ = mean_f(gp)
        Σ = opt_diag(gp.κ*gp.Σ,gp.κ)
        for i in 1:model.inference.nMinibatch
            x = model.inference.nodes*sqrt(max(Σ[i],zero(T))) .+ μ[i]
            tot += dot(model.inference.weights,logpdf.(model.likelihood,y[i],x))
        end
    end
    return tot
end

function compute_grad_expectations!(model::AbstractGP{T,L,<:QuadratureVI}) where {T,L}
    y = get_y(model)
    for (gp,opt) in zip(model.f,model.inference.vi_opt)
        μ = mean_f(gp)
        Σ = diag_cov_f(gp)
        for i in 1:model.inference.nMinibatch
            opt.ν[i], opt.λ[i] = grad_quad(model.likelihood,y[i],μ[i],Σ[i],model.inference)
        end
    end
end

#Compute the first and second derivative of the log-likelihood using the quadrature nodes
function grad_quad(l::Likelihood{T},y::Real,μ::Real,σ²::Real,i::Inference) where {T<:Real}
    x = i.nodes*sqrt(max(σ²,zero(T))) .+ μ
    Edlogpdf = dot(i.weights,grad_logpdf.(l,y,x))
    Ed²logpdf = dot(i.weights,hessian_logpdf.(l,y,x))
    if i.clipping != 0
        return (abs(Edlogpdf) > i.clipping ? sign(Edlogpdf)*i.clipping : -Edlogpdf::T,
                abs(Ed²logpdf) > i.clipping ? sign(Ed²logpdf)*i.clipping : -Ed²logpdf::T)
    else
        return -Edlogpdf::T, Ed²logpdf::T
    end
end


text:
using Random

Random.rand!(A::MtlArray) = Random.rand!(GPUArrays.default_rng(MtlArray), A)
Random.randn!(A::MtlArray) = Random.randn!(GPUArrays.default_rng(MtlArray), A)

ext: jl


text:
module Counter

export create # API
export add_to_file, read_from_file, create_new_invoice_nbr # Infrastructure

include("Domain.jl"); using .Domain
include("API.jl"); using .API
include("Infrastructure.jl"); using .Infrastructure

end


text:
module LocalDataCommon

using Distributions
using LinearAlgebra: norm
using ProgressMeter
using BlackBoxOptim

export estimate_admitted_discharged
export estimate_discharged
export estimate_admitted_uncertainty
export build_carepaths_data

function estimate_admitted_discharged(data; use_rounding=false, optim_time=1.0)
	N, T = data.N, data.T

	function unpack_combined_params(i, params)
		_admitted = params[T+1:end]
		_discharged = params[1:T]

		_discharged = _discharged .* (data.initial[i] / sum(_discharged))

		if use_rounding
			_admitted = round.(Int, _admitted)
			_discharged = round.(Int, _discharged)
		end

		return (
			admitted = _admitted,
			discharged = _discharged,
		)
	end

	function combined_score_func(i, params; s_penalty=0.0)
		z = unpack_combined_params(i, params)
		_admitted, _discharged = z.admitted, z.discharged

		_active = compute_active_nosent(data.initial[i], _discharged, _admitted, data.los_dist, use_rounding=use_rounding)
		score = norm(_active - data.active[i,:], 2)
		smooth_score = norm(_active[2:end] - _active[1:end-1], 2)
		return score + (s_penalty * smooth_score)
	end

	function combined_param_bounds(i)
		_admitted_ub = fill(0.2 * maximum(data.active[i,:]), T)
		_discharged_ub = 1 .- cdf.(data.los_dist, 1:T)
		ub = vcat(_discharged_ub, _admitted_ub)
		bds = collect(zip(zeros(Float64, length(ub)), ub))
		return bds
	end

	admitted_sim = zeros(Float64, N, T)
	discharged_sim = zeros(Float64, N, T)
	@showprogress for i in 1:N
		if maximum(data.active[i,:]) == 0.0
			continue
		end

		local_score_func = params -> combined_score_func(i, params, s_penalty=0.0)
		local_param_bounds = combined_param_bounds(i)

		results_combined = bboptimize(
			local_score_func,
			SearchRange = local_param_bounds,
			Method = :adaptive_de_rand_1_bin_radiuslimited,
			TraceMode = :silent,
			MaxTime = optim_time,
			RandomizeRngSeed = false,
			RngSeed = 0,
		)

		_local_params = best_candidate(results_combined)
		local_params = unpack_combined_params(i, _local_params)
		local_admitted, local_discharged = local_params.admitted, local_params.discharged
		admitted_sim[i,:] = local_admitted
		discharged_sim[i,:] = local_discharged
	end

	if use_rounding
		admitted_sim = Int.(admitted_sim)
		discharged_sim = Int.(discharged_sim)
	end

	return (
		admitted = admitted_sim,
		discharged = discharged_sim,
	)
end

function estimate_discharged(data; use_rounding=false, optim_time=1.0)
	N, T = data.N, data.T

	function unpack_params(i, _discharged)
		_discharged = _discharged .* (data.initial[i] / sum(_discharged))
		if use_rounding
			_discharged = round.(Int, _discharged)
		end
		return _discharged
	end

	function score_func(i, params; s_penalty=0.0)
		_discharged = unpack_params(i, params)
		_active = compute_active_nosent(data.initial[i], _discharged, data.admitted[i,:], data.los_dist, use_rounding=use_rounding)
		score = norm(_active - data.active[i,:], 2)
		smooth_score = norm(_active[2:end] - _active[1:end-1], 2)
		return score + (s_penalty * smooth_score)
	end

	function param_bounds(i)
		ub = 1 .- cdf.(data.los_dist, 1:T)
		bds = collect(zip(zeros(Float64, length(ub)), ub))
		return bds
	end

	discharged_sim = zeros(Float64, N, T)
	@showprogress for i in 1:N
		if maximum(data.active[i,:]) == 0.0
			continue
		end

		local_score_func = params ->score_func(i, params, s_penalty=0.0)
		local_param_bounds = param_bounds(i)

		r = bboptimize(
			local_score_func,
			SearchRange = local_param_bounds,
			Method = :adaptive_de_rand_1_bin_radiuslimited,
			TraceMode = :silent,
			MaxTime = optim_time,
			RandomizeRngSeed = false,
			RngSeed = 0,
		)

		_local_params = best_candidate(r)
		discharged_sim[i,:] = unpack_params(i, _local_params)
	end

	if use_rounding
		discharged_sim = Int.(discharged_sim)
	end

	return discharged_sim
end

function compute_active_nosent(_initial::Array{<:Real,1}, _discharged::Array{<:Real,2},
		_admitted::Array{<:Real,2}, los_dist::Distribution; use_rounding=false)
	N, T = size(_admitted)
	_L = 1.0 .- cdf.(los_dist, 0:T)
	_active = [(
		_initial[i]
		- sum(_discharged[i,1:t])
		+ sum(_L[t-t₁+1] * _admitted[i,t₁] for t₁ in 1:t)
	) for i in 1:N, t in 1:T]
	if use_rounding
		_active = round.(Int, _active)
	end
	return _active
end

function compute_active_nosent(_initial::Real, _discharged::Array{<:Real,1},
		_admitted::Array{<:Real,1}, los_dist::Distribution; use_rounding=false)
	T = length(_admitted)
	_L = 1.0 .- cdf.(los_dist, 0:T)
	_active = [(
		_initial
		- sum(_discharged[1:t])
		+ sum(_L[t-t₁+1] * _admitted[t₁] for t₁ in 1:t)
	) for t in 1:T]
	if use_rounding
		_active = round.(Int, _active)
	end
	return _active
end

function estimate_admitted_uncertainty(data, admitted)
	admitted_uncertainty = 0.05 * maximum(admitted, dims=2)[:]
	admitted_uncertainty = repeat(admitted_uncertainty, 1, data.T)
	admitted_uncertainty = min.(admitted_uncertainty, admitted)
	return admitted_uncertainty
end

function roundint(x::Real, direction::Symbol=:nearest)::Int
	rounddir = RoundNearest
	if direction == :down
		rounddir = RoundDown
	elseif direction == :up
		rounddir = RoundUp
	end
	return round(Int, x, rounddir)
end

function build_carepaths_data(data)
	N, T = data.N, data.T

	stack0(A...) = permutedims(cat(A..., dims=ndims(A[1])+1), [ndims(A[1])+1; 1:ndims(A[1])])
	beds_grouped = stack0(data.icu.beds, data.ward.beds)
	admitted_grouped = stack0(zeros(Float64, N, T), data.icu.admitted, zeros(Float64, N, T), data.ward.admitted)
	admitted_uncertainty_grouped = stack0(zeros(Float64, N, T), data.icu.admitted_uncertainty, zeros(Float64, N, T), data.ward.admitted_uncertainty)
	initial_grouped = stack0(data.icu.initial, zeros(Float64, N), zeros(Float64, N), data.ward.initial)
	discharged_grouped = stack0(data.icu.discharged, zeros(Float64, N, T), zeros(Float64, N, T), data.ward.discharged)

	group_transfer_graph = BitArray([
		0 0 1 0
		1 0 0 0
		0 0 0 0
		0 0 0 0
	])

	_p = zeros(Float64, T+1)
	_p[2+1] = 1
	los_dist_toicu = DiscreteNonParametric(0:T, _p)
	_p = zeros(Float64, T+1)
	_p[3+1] = 1
	los_dist_fromicu = DiscreteNonParametric(0:T, _p)
	los_bygroup = [data.icu.los_dist, los_dist_toicu, los_dist_fromicu, data.ward.los_dist]

	carepath_data = merge(data, (
		G = 4,
		B = 2,

		beds = beds_grouped,
		admitted = admitted_grouped,
		admitted_uncertainty = admitted_uncertainty_grouped,
		initial = initial_grouped,
		discharged = discharged_grouped,

		groups = [:icu, :to_icu, :from_icu, :ward],
		bedtypes = [:icu, :ward],
		los_bygroup = los_bygroup,
		bedtype_bygroup = [1, 2, 2, 2],
		group_transfer_graph = group_transfer_graph,
	))
	return carepath_data
end

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
if Sys.islinux()
    root_dir = dirname(pwd())
    file_prefix = "COTAHIST_A"
    file_begin_counter = 1995
    file_end_counter = 2019

    path_sufix = "/Datasets/Dataset_Bovespa/"
	path_sender = "Zipped/"
	path_destiny = "Unzipped/"

	location = root_dir * path_sufix * path_sender
	destiny = root_dir * path_sufix * path_destiny

    println(location)
    println(destiny)

    if !isdir(destiny)
        mkdir(destiny)
    else
        rm(destiny, recursive=true)
        mkdir(destiny)
    end

    for year in range(file_begin_counter, stop=file_end_counter)
        path_sender_with_year = [location * file_prefix * string(year)]
        path_destiny_with_year = [destiny]
        full_command = `unzip $path_sender_with_year -d $path_destiny_with_year`
        run(full_command)
    end
else
    if Sys.iswindows()
        root_dir = dirname(pwd())
        file_prefix = "COTAHIST_A"
        file_begin_counter = 1995
        file_end_counter = 2019
        fileExtension = ".zip"

        path_sufix = "\\Datasets\\Dataset_Bovespa\\"
        path_sender = "Zipped\\"
        path_destiny = "Unzipped\\"

        location = root_dir * path_sufix * path_sender
        destiny = root_dir * path_sufix * path_destiny

        println(location)
        println(destiny)

        if !isdir(destiny)
            mkdir(destiny)
        else
            rm(destiny, recursive=true)
            mkdir(destiny)
        end

        #=
        powerShellCommand = 'Expand-Archive -Force '

	    for year in range(fileStartYear, fileEndYear + 1):
		    completeCommand = powerShellCommand + location + filePrefix + str(year) + fileExtension + ' ' + destiny
		    print(completeCommand)
		    result = subprocess.Popen(['powershell.exe',completeCommand])
		    print(result)
        print("Done.")
        =#

        #=
        for year in range(file_begin_counter, stop=file_end_counter)
            path_sender_with_year = [location * file_prefix * string(year)]
            path_destiny_with_year = [destiny]
            full_command = `unzip $path_sender_with_year -d $path_destiny_with_year`
            run(full_command)
        end
        =#
    else
        if Sys.isapple()
            println("Macbook")
            root_dir = dirname(pwd())
            file_prefix = "COTAHIST_A"
            file_begin_counter = 1995
            file_end_counter = 2019

            path_sufix = "/Datasets/Dataset_Bovespa/"
	        path_sender = "Zipped/"
	        path_destiny = "Unzipped/"

	        location = root_dir * path_sufix * path_sender
	        destiny = root_dir * path_sufix * path_destiny

            println(location)
            println(destiny)

            if !isdir(destiny)
                mkdir(destiny)
            else
                rm(destiny, recursive=true)
                mkdir(destiny)
            end
            #=
            for year in range(file_begin_counter, stop=file_end_counter)
                path_sender_with_year = [location * file_prefix * string(year)]
                path_destiny_with_year = [destiny]
                full_command = `unzip $path_sender_with_year -d $path_destiny_with_year`
                run(full_command)
            end
            =#
        else
            println("System not implemented")
        end
    end
end

println("Done.")


text:
@testset "$TEST $G" begin

g10 = CompleteGraph(10, G)
@test local_clustering_coefficient(g10) == ones(10)
@test global_clustering_coefficient(g10) == 1
@test local_clustering(g10) == (fill(36, 10), fill(36, 10))
@test triangles(g10) == fill(36, 10)
@test triangles(g10, 1) == 36

end # testset


text:
@testset "util.jl" begin

    @testset "gridmake" begin
        want = [1 4 6
                2 4 6
                3 4 6
                1 5 6
                2 5 6
                3 5 6
                1 4 7
                2 4 7
                3 4 7
                1 5 7
                2 5 7
                3 5 7]

        # test allocating version
        @test want == @inferred gridmake([1, 2, 3], [4, 5], [6, 7])

        # test non-allocating version
        out = zeros(Int, 12, 3)
        @inferred gridmake!(out, [1, 2, 3], [4, 5], [6, 7])
        @test out == want

        # test single array version returns matrix of same type
        @test gridmake([1, 2, 3]) == [1 2 3]'

        # make sure we get the error we were expecting
        @test_throws AssertionError gridmake!(zeros(Int, 2, 2), [1, 2, 3], [4, 5])
    end

    @testset "fix" begin
        @test 1 == @inferred QuantEcon.fix(1.2)
        @test [0, 2] == @inferred QuantEcon.fix([0.9, 2.1])
        @test [0, 2] == @inferred QuantEcon.fix([0, 2])

        out = [100, 100]
        @inferred QuantEcon.fix!([0, 2], out)
        @test [0, 2] == out

        @test [0 2; 2 0] == @inferred QuantEcon.fix([0.5 2.9999; 3-2eps() -0.9])
    end

    @test ([1 2; 1 2], [3 3; 4 4]) == @inferred meshgrid([1, 2], [3, 4])

    @testset "simplex_tools" begin
        grid_3_4 = [0  0  0  0  0  1  1  1  1  2  2  2  3  3  4
                    0  1  2  3  4  0  1  2  3  0  1  2  0  1  0
                    4  3  2  1  0  3  2  1  0  2  1  0  1  0  0]
        points = [0 1 4
                  0 1 0
                  4 2 0]
        idx = Array{Int}(3)

        for i in 1:3
            idx[i] = simplex_index(points[:, i], 3, 4)
        end

        @test all(simplex_grid(3, 4) .== grid_3_4)
        @test all(grid_3_4[:, idx] .== points)
        @test size(grid_3_4, 2) == num_compositions(3, 4)

        # Output from QuantEcon.py
        grid_5_4 =
            [0 0 0 0 4 0 0 0 1 3 0 0 0 2 2 0 0 0 3 1 0 0 0 4 0
             0 0 1 0 3 0 0 1 1 2 0 0 1 2 1 0 0 1 3 0 0 0 2 0 2
             0 0 2 1 1 0 0 2 2 0 0 0 3 0 1 0 0 3 1 0 0 0 4 0 0
             0 1 0 0 3 0 1 0 1 2 0 1 0 2 1 0 1 0 3 0 0 1 1 0 2
             0 1 1 1 1 0 1 1 2 0 0 1 2 0 1 0 1 2 1 0 0 1 3 0 0
             0 2 0 0 2 0 2 0 1 1 0 2 0 2 0 0 2 1 0 1 0 2 1 1 0
             0 2 2 0 0 0 3 0 0 1 0 3 0 1 0 0 3 1 0 0 0 4 0 0 0
             1 0 0 0 3 1 0 0 1 2 1 0 0 2 1 1 0 0 3 0 1 0 1 0 2
             1 0 1 1 1 1 0 1 2 0 1 0 2 0 1 1 0 2 1 0 1 0 3 0 0
             1 1 0 0 2 1 1 0 1 1 1 1 0 2 0 1 1 1 0 1 1 1 1 1 0
             1 1 2 0 0 1 2 0 0 1 1 2 0 1 0 1 2 1 0 0 1 3 0 0 0
             2 0 0 0 2 2 0 0 1 1 2 0 0 2 0 2 0 1 0 1 2 0 1 1 0
             2 0 2 0 0 2 1 0 0 1 2 1 0 1 0 2 1 1 0 0 2 2 0 0 0
             3 0 0 0 1 3 0 0 1 0 3 0 1 0 0 3 1 0 0 0 4 0 0 0 0]
        @test simplex_grid(5, 4) == reshape(grid_5_4', 5, 70)
    end

    @testset "next_k_array" begin

        k_arrays = [1  2  3;
                    1  2  4;
                    1  3  4;
                    2  3  4;
                    1  2  5;
                    1  3  5;
                    2  3  5;
                    1  4  5;
                    2  4  5;
                    3  4  5;
                    1  2  6;
                    1  3  6;
                    2  3  6;
                    1  4  6;
                    2  4  6;
                    3  4  6;
                    1  5  6;
                    2  5  6;
                    3  5  6;
                    4  5  6]
        L, k = size(k_arrays)

        k_arrays_computed = similar(k_arrays)
        k_arrays_computed[1, :] = collect(1:k)
        @test k_array_rank(k_arrays_computed[1, :]) == 1
        for i = 2:L
            k_arrays_computed[i, :] = next_k_array!(k_arrays_computed[i-1, :])
            @test k_array_rank(k_arrays_computed[i, :]) == i
        end

        @test k_arrays_computed == k_arrays

        # test InexactError when ranking is out of range
        n, k = 100, 50
        @test k_array_rank(BigInt, collect(n-k+1:n)) == binomial(BigInt(n), BigInt(k))
        @test_throws InexactError k_array_rank(collect(n-k+1:n))

        # test returning wrong value when ranking is out of range
        n, k = 68, 29
        @test k_array_rank(BigInt, collect(n-k+1:n)) == binomial(BigInt(n), BigInt(k))
        @test k_array_rank(collect(n-k+1:n)) != binomial(BigInt(n), BigInt(k))

    end

end


text:
module WeightedNetwork

using SparseArrays
using Graphs
using SimpleWeightedGraphs
using DataStructures
using LinearAlgebra
using Plots
using GraphPlot
using DataStructures
using ShowGraphviz
# using Cairo
# using Fontconfig
using Random
using Tokenize

export
    # generators
    rand_network, rand_directed_network, rand_dag, mk_connected_layered_dag,

    # graph properties
    children, parent, neighbors,

    # dataformat
    adjacency_matrix, AdjacencyList,

    # visualization mit GraphPlot:
    plot_network, plot_graph,

    # visualization mit Graphviz (for small graphs or use file_export:)
    plot_graphviz, write_dot_file, read_dot_file, AttributeDict, default_attributes,

    # algorithm:
    dfs_path, bfs_path, toplogicalsort, node_color, path_to_adj_matrix, shortest_path

include("./data_types/graphdata_types.jl")
include("./generators/random_networks.jl")
include("./visualization/gplots.jl")
include("./visualization/dots.jl")
include("./algorithm/traversel.jl")
include("./algorithm/graph_utils.jl")

end

max_stars_repo_name: tragisch/NetworkOptimization

max_issues_repo_name: tragisch/NetworkOptimization

max_forks_repo_name: tragisch/NetworkOptimization


text:
using Test: @testset, @test

using StatsBase

using ScottishTaxBenefitModel
using .FRSHouseholdGetter

using .ExampleHouseholdGetter
using .ModelHousehold:
      BenefitUnit,
      BenefitUnits,
      BUAllocation,
      Household,
      People_Dict,
      PeopleArray,
      Person,
      count,
      default_bu_allocation,
      get_benefit_units,
      get_head,
      get_spouse,
      has_carer_member,
      has_disabled_member,
      is_lone_parent,
      num_children,
      num_people,
      printpids

using .RunSettings: Settings
using .TimeSeriesUtils: fy_from_bits

using .Definitions

using .HistoricBenefits:
      RATIO_BENS,
      make_benefit_ratios!
      using .ExampleHelpers

using DataFrames

start_year=2015
num_households = 0
total_num_people = 0
nhh2 = 0
# pyplot()

rc = @timed begin
      num_households,total_num_people,nhh2 = FRSHouseholdGetter.initialise( Settings() )
end
println( "num_households=$num_households, num_people=$(total_num_people)")

@testset "ratio tests" begin
      hh = ExampleHouseholdGetter.get_household( "single_parent_1" )
      bu = get_benefit_units( hh )[1]
      head = get_head( bu )
      fy = fy_from_bits( hh.interview_year, hh.interview_month )
      @test fy == 2019

      # pension fy 2019 = 129.2
      # wid fy 2019 = 119.9
      empty!(head.income)
      empty!(head.benefit_ratios)
      head.income[state_pension] = 129.20
      make_benefit_ratios!(head, hh.interview_year, hh.interview_month )
      println( "head.benefit_ratios = $(head.benefit_ratios)")
      @test length( head.benefit_ratios ) == 1
      @test head.benefit_ratios[state_pension] ≈ 1

      head.income[bereavement_allowance_or_widowed_parents_allowance_or_bereavement] = 129.20
      make_benefit_ratios!(head, hh.interview_year, hh.interview_month )
      @test length( head.benefit_ratios ) == 2
      # 129.2/119.9 \approx 1.0775646371976646
      @test head.benefit_ratios[bereavement_allowance_or_widowed_parents_allowance_or_bereavement] ≈ 1.0775646371976646

      head.income[dlamobility] = 62.25
      head.income[attendance_allowance] = 89.15
      head.income[dlaself_care] = 89.15
      make_benefit_ratios!(head, hh.interview_year, hh.interview_month )
      @test head.dla_mobility_type == high
      @test head.attendance_allowance_type == high
      @test head.dla_self_care_type == high

      head.income[dlamobility] = 23.60
      head.income[attendance_allowance] = 59.7
      head.income[dlaself_care] = 23.60
      make_benefit_ratios!(head, hh.interview_year, hh.interview_month )
      @test head.dla_mobility_type == low
      @test head.attendance_allowance_type == low
      @test head.dla_self_care_type == low
      head.income[dlaself_care] = 59.7
      make_benefit_ratios!(head, hh.interview_year, hh.interview_month )
      @test head.dla_self_care_type == mid

      # overall test
      num_rats = Dict()
      exact = 0
      tot = 0
      for target in RATIO_BENS
            num_rats[target] = Vector{Float64}(undef,0)
      end
      for hhno in 1:num_households
            hh = FRSHouseholdGetter.get_household( hhno )
            fy = fy_from_bits( hh.interview_year, hh.interview_month )
            # println("people in HH")
            # printpids(hh.people)
            for (pid,pers) in hh.people
                  for target in RATIO_BENS
                        if haskey( pers.benefit_ratios, target )
                              @test pers.income[target] > 0
                              pts = pers.benefit_ratios[target]
                              if pts ≈ 1
                                    exact += 1
                              end
                              tot += 1
                              push!(num_rats[target], pts )
                              if target == bereavement_allowance_or_widowed_parents_allowance_or_bereavement
                                    @test pers.bereavement_type in [
                                          bereavement_allowance,
                                          widowed_parents,
                                          bereavement_support]
                              end
                        end
                  end
                  #=
                  newrats = make_benefit_ratios( fy, pers.income )
                  @test length(symdiff(keys( newrats ),keys( pers.benefit_ratios ))) == 0 # all same keys
                  for k in keys( newrats )
                        @test  newrats[k] ≈ pers.benefit_ratios[k]
                  end
                  =#
            end
      end
      for target in RATIO_BENS
            rat = num_rats[target]
            print(typeof(rat))
            n = length( rat )
            println( "ratios for $target n = $n")
            nb = target == state_pension ? 40 : 6
            cm = fit(Histogram,rat, nbins=nb)
            #=
            p = Plots.plot( cm, title="$target ratio of actual vs standard for year"  )
            println( "hist:\n $cm")
            fname="tmp/$(target)_hist.svg"
            Plots.svg( p, fname)
            =#
            println( "exacts = $exact tot = $tot")
      end
end

@testset "people search functions" begin
      @time names = ExampleHouseholdGetter.initialise( Settings() )
      hh = ExampleHouseholdGetter.get_household( "single_parent_1" )
      printpids( hh.people )
      @test num_people( hh ) == 3
      @test is_lone_parent( hh )
      bus = get_benefit_units( hh )
      @test size(bus)[1] == 1
      bu = bus[1]
      @test is_lone_parent( bu )
      @test num_people( bu ) == 3
      head = get_head( bu )
      @test ! has_disabled_member( bu )
      @test ! has_carer_member( bu )
      # FIXME we'll need to change this if disability no longer based on benefit receipt
      head.registered_blind = true
      @test has_disabled_member( bu )
      @test ! has_carer_member( bu )
      head.registered_blind = false
      head.income[carers_allowance] = 1.0
      @test has_carer_member( bu )
      delete!(head.income,carers_allowance)
      @test ! has_carer_member( bu )
      head.income[severe_disability_allowance] = 1.0
      @test has_disabled_member( bu )

end

@testset "benefit unit allocations" begin
      people_count = 0
      over_16_kids = 0
      @time for hhno in 1:num_households
            hh = FRSHouseholdGetter.get_household( hhno )
            # println("people in HH")
            # printpids(hh.people)

            hhsize = size(collect(keys(hh.people)))[1]
            people_count += hhsize
            bus = get_benefit_units( hh )
            buallocation = default_bu_allocation( hh )
            nbus = size( bus )[1]
            @test 0 < nbus < 10
            by_bu_people_count = 0
            by_bu_child_count = 0
            i = 0
            for bu in bus
                  # println( "from BU")
                  # printpids(bu)

                  i += 1
                  bu_people_count = 1
                  bu_children = 0
                  head = get_head( bu )
                  @test head.age >= 16
                  spouse = get_spouse( bu )
                  if spouse !== nothing
                        bu_people_count += 1
                        @test spouse.age >= 16
                  end
                  for chno in bu.children
                        child = bu.people[chno]
                        @test child.age <= 19
                        if child.age >16
                            over_16_kids += 1
                        end
                        bu_people_count += 1
                        bu_children += 1
                  end
                  @test num_children( bu ) == bu_children
                  by_bu_child_count += bu_children
                  by_bu_people_count += bu_people_count
                  @test bu_people_count == num_people( bu )
                  @test bu_people_count == size( buallocation[i])[1]
            end
            @test by_bu_people_count == hhsize
            @test num_children(hh) == by_bu_child_count

            #
            # Test splitting up PIDs into source, year, etc.
            #
            for (pid,pers) in hh.people
                  bits = from_pid( pid )
                  @test pid == pers.pid
                  @test bits.year == hh.data_year
                  @test bits.datasource == FRS
                  @test bits.hid == hh.hid
            end
            bua_people_count = 0
            for bua in buallocation
                  bua_people_count += size( bua )[1]
            end

            # println("BU Allocation")
            # printpids(buallocation)
            @test bua_people_count == hhsize
      end
      @test people_count == total_num_people
      @test over_16_kids > 0
      println( "num people $people_count over_16_kids=$over_16_kids" )
end

@testset "aggregate nmt bens" begin
      n = num_households
      cases = DataFrame(
            hid = zeros(BigInt, n ),
            year = zeros(Int,n),
            weight = zeros(n),

            pip_daily_living_enhanced = zeros(n),
            pip_daily_living_standard = zeros(n),
            pip_mobility_enhanced = zeros(n),
            pip_mobility_standard = zeros(n),
            carers_allowance = zeros(n),
            attendance_allowance_high = zeros(n),
            attendance_allowance_low = zeros(n),
            dla_care_high = zeros(n),
            dla_care_mid = zeros(n),
            dla_care_low = zeros(n),
            dla_mobility_high = zeros(n),
            dla_mobility_low  = zeros(n))
      r = 0
      @time for hhno in 1:num_households
            hh = FRSHouseholdGetter.get_household( hhno )
            r += 1
            row = cases[r,:]
            row.weight = hh.weight;
            row.hid = hh.hid
            row.year = hh.interview_year
            println( "on hh $(hh.hid)")
            for (pid, pers) in hh.people
                  print( "mobtype $(pers.dla_mobility_type)")
                  if pers.dla_self_care_type == low
                        row.dla_care_low += 1
                  elseif pers.dla_self_care_type == mid
                        row.dla_care_mid += 1
                  elseif pers.dla_self_care_type == high
                        row.dla_care_high += 1
                  end
                  if pers.dla_mobility_type == low
                        row.dla_mobility_low += 1
                  elseif pers.dla_mobility_type == mid
                        @assert 1==2 "mobility mid should never happen"
                  elseif pers.dla_mobility_type == high
                        row.dla_mobility_high += 1
                  end
                  if pers.attendance_allowance_type == low
                        row.attendance_allowance_low += 1
                  elseif pers.attendance_allowance_type == mid
                        @assert 3==4 "attendance mid should never happen"
                  elseif pers.attendance_allowance_type == high
                        row.attendance_allowance_high += 1
                  end
                  if pers.pip_daily_living_type == standard_pip
                        row.pip_daily_living_standard += 1
                  elseif pers.pip_daily_living_type == enhanced_pip
                        row.pip_daily_living_enhanced += 1
                  end
                  if pers.pip_mobility_type == standard_pip
                        row.pip_mobility_standard += 1
                  elseif pers.pip_mobility_type == enhanced_pip
                        row.pip_mobility_enhanced += 1
                  end
                  if haskey( pers.income, carers_allowance )
                        row.carers_allowance += 1
                  end
            end
      end

      pip_daily_living_enhanced_tot = cases.pip_daily_living_enhanced'cases.weight
      println( "pip_daily_living_enhanced total = $pip_daily_living_enhanced_tot" )

      pip_daily_living_standard_tot = cases.pip_daily_living_standard'cases.weight
      println( "pip_daily_living_standard total = $pip_daily_living_standard_tot" )

      pip_mobility_enhanced_tot = cases.pip_mobility_enhanced'cases.weight
      println( "pip_mobility_enhanced total = $pip_mobility_enhanced_tot" )

      pip_mobility_standard_tot = cases.pip_mobility_standard'cases.weight
      println( "pip_mobility_standard total = $pip_mobility_standard_tot" )

      carers_allowance_tot = cases.carers_allowance'cases.weight
      println( "carers_allowance total = $carers_allowance_tot" )

      attendance_allowance_high_tot = cases.attendance_allowance_high'cases.weight
      println( "attendance_allowance_high total = $attendance_allowance_high_tot" )

      attendance_allowance_low_tot = cases.attendance_allowance_low'cases.weight
      println( "attendance_allowance_low total = $attendance_allowance_low_tot" )

      dla_care_high_tot = cases.dla_care_high'cases.weight
      println( "dla_care_high total = $dla_care_high_tot" )

      dla_care_mid_tot = cases.dla_care_mid'cases.weight
      println( "dla_care_mid total = $dla_care_mid_tot" )

      dla_care_low_tot = cases.dla_care_low'cases.weight
      println( "dla_care_low total = $dla_care_low_tot" )

      dla_mobility_high_tot = cases.dla_mobility_high'cases.weight
      println( "dla_mobility_high total = $dla_mobility_high_tot" )

      dla_mobility_low_tot = cases.dla_mobility_low'cases.weight
      println( "dla_mobility_low  total = $dla_mobility_low_tot" )

      CSV.write( "bencounts.csv", cases )

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))

plt.hist(data_nmiss.plant_age, color="blue", edgecolor="white", bins=100,
      density=true, alpha=0.5)
plt.axvline(mean_age, label = "Mean", color = "red")
plt.axvline(median_age, label = "Median")

plt.legend()

plt.xlim(0,)

savefig(joinpath(@OUTPUT, "D0-processing-g2.svg")) # hide


text:
module Particles

using FFTW
using OffsetArrays
using StatsBase

export Equation, Simulation
export ExplicitEuler
export PoissonSolver, solve!, eval_field

include("equation.jl")
include("integrators.jl")
include("simulation.jl")
include("poisson.jl")

end # module

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
import Pluto:throttled

@testset "Throttled" begin
    x = Ref(0)

    function f()
    	x[] += 1
    end
    f()
    # f was not throttled
    @test x[] == 1

    dt = 4 / 100
    ft, flush = throttled(f, dt)

    for x in 1:10
    	ft()
    end
    # we have an initial cooldown period in which f should not fire...
    # ...so x is still 1...
    @test x[] == 1
    sleep(2dt)
    # ...but after a delay, the call should go through.
    @test x[] == 2

    # sleep(0) ## ASYNC MAGIC :(

    # at this point, the *initial* cooldown period is over
    # and the cooldown period for the first throttled calls is over

    for x in 1:10
    	ft()
    end
    # we want to send plots to the user as soon as they are available,
    # so no leading timeout
    @test x[] == 3
    # the 2nd until 10th calls were still queued
    sleep(2dt)
    @test x[] == 4

    for x in 1:5
    	ft()
    	sleep(1.5dt)
    end
    @test x[] == 9
    sleep(2dt)
    @test x[] == 9

    ###

    # "call 1"
    ft()
    # no leading timeout, immediately set to 10
    @test x[] == 10

    sleep(.5dt)

    # "call 2"
    ft()
    # throttled
    @test x[] == 10

    sleep(.7dt)

    # we waited 1.2dt > dt seconds since "call 1", which should have started the dt cooldown. "call 2" landed during that calldown, and should have triggered by now
    @test x[] == 11

    sleep(2dt)
    @test x[] == 11

    ###

    ft()
    ft()
    @test x[] == 12
    flush()
    @test x[] == 13
    sleep(2dt)
    @test x[] == 13

end


text:
using FPGAUtils
using Test

@testset "FPGAUtils.jl" begin
    # Write your tests here.
end


text:
# SPDX-License-Identifier: X11
# 2020-11-20
# gacha (300pt)

function main()
  N = parse(Int, readline())
  d = Dict()
  for i ∈ 1:N
    d[readline()] = true
  end
  (println ∘ length)(d)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
# TODO(odow):
#
# This script can be used to build the C interface to Clp. However, it requires
# you to manually do the following steps first:
#
# 1) Copy Clp_C_Interface.h from clp into this /scripts directory
# 2) Copy Coin_C_defines.h from coinutils into this /scripts directory
#
# It should be possible to build the wrapper using the jll's, but I couldn't
# figure out how to do that, and I didn't have enough time to spend on it.

import Clang

LIBCLP_HEADERS = [
    joinpath(@__DIR__, "Coin_C_defines.h"),
    joinpath(@__DIR__, "Clp_C_Interface.h"),
]

wc = Clang.init(
    headers = LIBCLP_HEADERS,
    output_file = joinpath(@__DIR__, "..", "src", "libclp_api.jl"),
    common_file = joinpath(@__DIR__, "..", "src", "libclp_common.jl"),
    header_wrapped = (root, current) -> root == current,
    header_library = x -> "libClp",
    clang_diagnostics = true,
)

run(wc)


text:
"""
Submodule of Utilities that groups the utilities related to command-line
parsing.
"""
module Args

using DocStringExtensions # for TYPEDFIELDS
using ArgParse
export Arg, accepted_arg_list, throw_if_incompatible_options
export create_normalized_arg_subset

"""
An argument with `name`, `default` value, and `help` message.

$(TYPEDFIELDS)

"""
struct Arg{T}
  "The name with no initial double dashes but with dashes instead of spaces."
	name :: String
	"The default value if the argument is not specified. Use `false` for flags."
	default :: T
	"The help message to be displayed about this argument."
	help :: String
end

# The re-definition of Base.iterate is necessary to allow the flattening
# of arrays of arrays of Args.
Base.iterate(a :: Arg) = (a, nothing)
Base.iterate(a :: Arg, nothing) = nothing

"""
    accepted_arg_list(::Val{T}) :: Vector{Arg} where {T}

Generic error fallback. Any model or solver to be supported by
`GuillotineModules.SolversArgs.run` should implement their own version of this
method (replacing the `T` in `::Val{T}` by a Symbol identifying the solver
package or the name of the model). Look at module SolversArgs source for
examples of the implementation for a solver, and at PPG2KP module source for
examples for a model.

Gives the list of accepted options by some model or solver. The options have
either a boolean default (in this case they take no argument, the presence of
the option just flips the default value) or a non-boolean default (in this
case, if they are passed in the command-line, they must have a parameter; if
they are not passed, the default value is used).  Every solver implementation
must support the `no-output` and preferably the `raw-argument` option too. The
name of the supported solvers and implemented methods need to be passed to the
`run` method for them to be considered by it.
"""
function accepted_arg_list(::Val{T}) :: Vector{Arg} where {T}
	@error(
		"Solver " * string(T) * " is not supported (i.e., there is not" *
		" an implementation of Julia method `GuillotineModels.Utilities." *
		".Args.accepted_arg_list(::Val{Symbol(\"$(string(T))\")})` for it) or" *
		" the solver package was not imported before this method was called."
	)
end

"""
    throw_if_incompatible_options(::Val{T}, p_args) where {T}

Generic error fallback. Any model or solver to be supported by
`GuillotineModels.SolversArgs.run` should implement their own version of this
method (replacing the `T` in `::Val{T}` by a Symbol identifying the solver
package or the name of the model). Look at module SolversArgs source for
examples of the implementation for a solver, and at PPG2KP module source for
examples for a model. Every solver implementation must support the `no-output`
option, and there is no need to guarantee that conflicts with involving
the `raw-argument` option are detected.

Often is specialized to an empty method, as is not so often that solvers
or models have conflicting options.
"""
function throw_if_incompatible_options(
	::Val{T}, p_args
) where {T}
	@error (
		"A specialized method for $(T) should exist, but instead this " *
		" generic error fallback was called."
	)
end

"""
ArgParse.add_arg_table!(settings :: ArgParseSettings, arg :: Arg)

A specialization of `ArgParse.add_arg_table!` to transform `Arg` objects into
options of `ArgParseSettings`. Boolean arguments become `:store_{true|false}`
options (depending on the default) and non-boolean arguments become
`:store_arg` options with a default (and enforcing the same arg_type of the
default).
"""
function ArgParse.add_arg_table!(settings :: ArgParseSettings, arg :: Arg)
	if isa(arg.default, Bool)
		conf = Dict{Symbol, Any}(
			:help => arg.help,
			# default == true means that passing it as a flag turns it false
			:action => (arg.default ? :store_false : :store_true)
		)
	else
		conf = Dict{Symbol, Any}(
			:help => arg.help,
			:arg_type => typeof(arg.default),
			:default => arg.default,
			:action => :store_arg
		)
	end
	# only long options are accepted, enforced clarity
	ArgParse.add_arg_table!(settings, "--" * arg.name, conf)
end

"""
    create_normalized_arg_subset(p_args, selected :: Vector{Arg})

Create a new `Dict` with keys equal to `selected` names and values equal to
`selected` defaults, except that if `p_args` has a key with the same name
as an argument then the value in `p_args` is used instead the default.

It shows an warning for every key in `p_args` that is not a name in
`selected`.

If a `quiet` flag is recognized by selected and is true (either by explicit
definition or by default), then the warnings are not printed.
"""
function create_normalized_arg_subset(
	p_args, selected :: Vector{Arg}
)
	new_dict = empty(p_args)
	for arg in selected
		if haskey(p_args, arg.name)
			new_dict[arg.name] = p_args[arg.name]
		else
			new_dict[arg.name] = arg.default
		end
	end
	quiet = get(new_dict, "quiet", false)
	if !quiet
		sel_keys = sort(getfield.(selected, :name))
		old_keys = keys(p_args)
		for k in old_keys
			isempty(searchsorted(sel_keys, k)) && @warn(
				"Key $(k) is not recognized, it will not be used."
			)
		end
	end

	return new_dict
end
end # submodule Args


text:
# Autogenerated wrapper script for GnuPG_jll for x86_64-w64-mingw32
export dirmngr, dirmngr_client, gpg, gpg_agent, gpg_connect_agent, gpgconf, gpgscm, gpgsm, gpgtar, gpgv, kbxutil

using GnuTLS_jll
using Libksba_jll
using Libgcrypt_jll
using Libgpg_error_jll
using nPth_jll
using Zlib_jll
using Libassuan_jll
using OpenLDAPClient_jll
using Bzip2_jll
using SQLite_jll
using libusb_jll
using Nettle_jll
JLLWrappers.@generate_wrapper_header("GnuPG")
JLLWrappers.@declare_executable_product(dirmngr)
JLLWrappers.@declare_executable_product(dirmngr_client)
JLLWrappers.@declare_executable_product(gpg)
JLLWrappers.@declare_executable_product(gpg_agent)
JLLWrappers.@declare_executable_product(gpg_connect_agent)
JLLWrappers.@declare_executable_product(gpgconf)
JLLWrappers.@declare_executable_product(gpgscm)
JLLWrappers.@declare_executable_product(gpgsm)
JLLWrappers.@declare_executable_product(gpgtar)
JLLWrappers.@declare_executable_product(gpgv)
JLLWrappers.@declare_executable_product(kbxutil)
function __init__()
    JLLWrappers.@generate_init_header(GnuTLS_jll, Libksba_jll, Libgcrypt_jll, Libgpg_error_jll, nPth_jll, Zlib_jll, Libassuan_jll, OpenLDAPClient_jll, Bzip2_jll, SQLite_jll, libusb_jll, Nettle_jll)
    JLLWrappers.@init_executable_product(
        dirmngr,
        "bin\\dirmngr.exe",
    )

    JLLWrappers.@init_executable_product(
        dirmngr_client,
        "bin\\dirmngr-client.exe",
    )

    JLLWrappers.@init_executable_product(
        gpg,
        "bin\\gpg.exe",
    )

    JLLWrappers.@init_executable_product(
        gpg_agent,
        "bin\\gpg-agent.exe",
    )

    JLLWrappers.@init_executable_product(
        gpg_connect_agent,
        "bin\\gpg-connect-agent.exe",
    )

    JLLWrappers.@init_executable_product(
        gpgconf,
        "bin\\gpgconf.exe",
    )

    JLLWrappers.@init_executable_product(
        gpgscm,
        "bin\\gpgscm.exe",
    )

    JLLWrappers.@init_executable_product(
        gpgsm,
        "bin\\gpgsm.exe",
    )

    JLLWrappers.@init_executable_product(
        gpgtar,
        "bin\\gpgtar.exe",
    )

    JLLWrappers.@init_executable_product(
        gpgv,
        "bin\\gpgv.exe",
    )

    JLLWrappers.@init_executable_product(
        kbxutil,
        "bin\\kbxutil.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MicrobiomeAnalysis
# using SummarizedExperiments
using Test

@testset "MicrobiomeAnalysis.jl" begin

    t, Xapp = LVmodel()

    @test @isdefined t
    @test @isdefined Xapp

end


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

name = "P4est"
version = v"2.2.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://p4est.github.io/release/p4est-2.2.tar.gz", "1549cbeba29bee2c35e7cc50a90a04961da5f23b6eada9c8047f511b90a8e438"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/p4est-2.2/
if [[ "${target}" == *-freebsd* ]]; then
  export LIBS="-lm"
elif [[ "${target}" == x86_64-linux-musl ]]; then
    # We can't run Fortran programs for the native platform, so a check that the
    # Fortran compiler works would fail.  Small hack: swear that we're
    # cross-compiling.  See:
    # https://github.com/JuliaPackaging/BinaryBuilderBase.jl/issues/50.
    sed -i 's/cross_compiling=no/cross_compiling=yes/' configure
    sed -i 's/cross_compiling=no/cross_compiling=yes/' sc/configure
fi
export CPPFLAGS="-I${includedir}"
export LDFLAGS="-L${libdir}"
export BLAS_LIBS="${libdir}/libopenblas.${dlext}"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-static
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=Sys.iswindows)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2"))
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
version = v"6.2.1"

include("../common.jl")

# Build the tarballs
build_tarballs(ARGS, configure(version)...;
               preferred_gcc_version=v"6", julia_compat="1.7")


text:
using MicroLogging

function checkDerivative(fullF::Function,x;doPlot::Bool=false, smoothOffset::Bool=false)
    df = fullF(x)[2]'
    f(x) = fullF(x)[1]
    return checkDerivative(f,df,x,doPlot=doPlot, smoothOffset=smoothOffset)
end
function checkDerivative(f::Function,df,x;doPlot::Bool=false, smoothOffset::Bool=false)
    N = 17
    h = exp10.(2:-1:-N)
    if smoothOffset
        v = smoothArray(randn(size(x)).-0.5,5,5.0)
    else
        v = randn(size(x)).-0.5
    end
    fx = f(x)
    errquad = zeros(length(h))
    errlin  = zeros(length(h))
    for i=1:length(h)
        errlin[i]  = norm(fx             - f(x+h[i]*v))
        errquad[i] = norm(fx + h[i]*df*v - f(x+h[i]*v))
        s = @sprintf(" h: %1.1e  ||f(x+h*v)||: %2.2e   elin: %2.2e   equad: %2.2e \n",h[i], norm(f(x+h[i]*v)), errlin[i],errquad[i])
        @info s
    end
    #
    if(doPlot)
      PyPlot.rc("legend",fontsize=10)
      PyPlot.rc("xtick", labelsize=10, color="black", direction="in")
      PyPlot.rc("ytick", labelsize=10, color="black", direction="in")
      #fig = PyPlot.figure(figsize=(6,5),facecolor="white")
      PyPlot.clf()
      PyPlot.plot(h,errlin,label="||f(x) - f(x+hv)||")
      PyPlot.plot(h,errquad,".-",label="||f(x) + h nabla f v - f(x+hv)||")
      PyPlot.yscale("log")
      PyPlot.xscale("log")
      PyPlot.xlabel("h")
      PyPlot.ylabel("error")
      PyPlot.legend(loc="upper left")
    end
    return errlin, errquad
end

function checkErrorDecay(errquad::Vector)
  errorReduction = log10.(errquad[1:end-1]./errquad[2:end])
  errorReduction[isinf.(errorReduction)]=1
  @debug "# errorReduction.>1.7 = ", sum(errorReduction.>1.7)
  @info "sum(errorReduction...) = ", sum(errorReduction[errorReduction.>1.7]), "   (needs >6)"
  return sum(errorReduction[errorReduction.>1.7])>6
end


text:
using MagneticReadHead: should_instrument, should_breakon, Rule, BreakpointRules
using Test

@testset "breakon rules" begin
    # Helper for defining rules during tests
    borule(args...) = BreakpointRules(Rule[], [Rule(args...)])

    @testset "module" begin
        # I don't really think this is useful, but it does work
        rules = borule(Iterators)
        @test should_breakon(rules, first(methods(Iterators.drop)), 1)
        @test !should_breakon(rules, first(methods(Iterators.drop)), 3)
        @test should_breakon(rules, first(methods(Iterators.flatten)), 1)
        @test !should_breakon(borule(Iterators), first(methods(Iterators.flatten)), 2)
    end

    @testset "Function" begin
        rules = borule(eps)
        for meth in methods(eps)
            @test should_breakon(rules, meth, 1)
            @test !should_breakon(rules, meth, 2)
            @test !should_breakon(rules, meth, 3)
        end
    end

    @testset "Method" begin
        meths  = collect(methods(pwd))
        rules = borule(meths[1])
        @test should_breakon(rules, meths[1], 1)
        @test !should_breakon(rules, meths[1], 2)
        for meth in meths[2:end]
            @test !should_breakon(rules, meth, 1)
        end
    end

    @testset "Method + statement number" begin
        # Not going to test ever other permutation of x + statement number
        # This is the one that matters, method + statement number
        # is what a line number breakpoint will become
        meths  = collect(methods(pwd))
        rules = borule(meths[1], 3)
        @test !should_breakon(rules, meths[1], 1)
        @test !should_breakon(rules, meths[1], 2)
        @test should_breakon(rules, meths[1], 3)
        @test !should_breakon(rules, meths[1], 4)

        for meth in meths[2:end]
            @test !should_breakon(rules, meth, 1)
            @test !should_breakon(rules, meth, 3)
        end
    end
end

@testset "default" begin
    rules = BreakpointRules()

    # should instrument everything
    @test should_instrument(rules, +)
    @test should_instrument(rules, sum)
    @test should_instrument(rules, Iterators.flatten)
    @test should_instrument(rules, Iterators.drop)

    @test !should_breakon(rules, first(methods(+)), 1)
    @test !should_breakon(rules, first(methods(+)), 3)
    @test !should_breakon(rules, first(methods(sum)), 1)
    @test !should_breakon(rules, first(methods(sum)), 4)
    @test !should_breakon(rules, first(methods(Iterators.flatten)), 1)
    @test !should_breakon(rules, first(methods(Iterators.flatten)), 5)
    @test !should_breakon(rules, first(methods(Iterators.drop)), 1)
    @test !should_breakon(rules, first(methods(Iterators.drop)), 5)
end

@testset "not instrumenting some modules" begin
    rules = BreakpointRules(
        [Rule(Iterators)],
        []
    )
    @test should_instrument(rules, +)
    @test should_instrument(rules, sum)
    @test !should_instrument(rules, Iterators.flatten)
    @test !should_instrument(rules, Iterators.drop)

    @test !should_breakon(rules, first(methods(+)), 1)
    @test !should_breakon(rules, first(methods(+)), 3)
    @test !should_breakon(rules, first(methods(sum)), 1)
    @test !should_breakon(rules, first(methods(sum)), 4)
    @test !should_breakon(rules, first(methods(Iterators.flatten)), 1)
    @test !should_breakon(rules, first(methods(Iterators.flatten)), 5)
    @test !should_breakon(rules, first(methods(Iterators.drop)), 1)
    @test !should_breakon(rules, first(methods(Iterators.drop)), 5)
end

@testset  "not instrumenting parent module should still instrument child module" begin
    #TODO: Decide if this is the correct behavour
    rules = BreakpointRules(
        [Rule(Base)],
        []
    )
    @test !should_instrument(rules, +)
    @test !should_instrument(rules, sum)
    @test should_instrument(rules, Iterators.flatten)
    @test should_instrument(rules, Iterators.drop)

    @test !should_breakon(rules, first(methods(+)), 1)
    @test !should_breakon(rules, first(methods(+)), 3)
    @test !should_breakon(rules, first(methods(sum)), 1)
    @test !should_breakon(rules, first(methods(sum)), 4)
    @test !should_breakon(rules, first(methods(Iterators.flatten)), 1)
    @test !should_breakon(rules, first(methods(Iterators.flatten)), 5)
    @test !should_breakon(rules, first(methods(Iterators.drop)), 1)
    @test !should_breakon(rules, first(methods(Iterators.drop)), 5)
end

@testset "not instrumenting modules but breaking on it" begin
    rules = BreakpointRules(
        [Rule(Iterators)],
        [Rule(Iterators.drop)]
    )
    @test should_instrument(rules, +)
    @test should_instrument(rules, sum)
    @test !should_instrument(rules, Iterators.flatten)
    @test should_instrument(rules, Iterators.drop)

    @test !should_breakon(rules, first(methods(+)), 1)
    @test !should_breakon(rules, first(methods(+)), 3)
    @test !should_breakon(rules, first(methods(sum)), 1)
    @test !should_breakon(rules, first(methods(sum)), 4)
    @test !should_breakon(rules, first(methods(Iterators.flatten)), 1)
    @test !should_breakon(rules, first(methods(Iterators.flatten)), 5)
    @test should_breakon(rules, first(methods(Iterators.drop)), 1)
    @test !should_breakon(rules, first(methods(Iterators.drop)), 5)
end


text:
using Mamba

## Data
rats = Dict{Symbol, Any}(
  :y =>
    [151, 199, 246, 283, 320,
     145, 199, 249, 293, 354,
     147, 214, 263, 312, 328,
     155, 200, 237, 272, 297,
     135, 188, 230, 280, 323,
     159, 210, 252, 298, 331,
     141, 189, 231, 275, 305,
     159, 201, 248, 297, 338,
     177, 236, 285, 350, 376,
     134, 182, 220, 260, 296,
     160, 208, 261, 313, 352,
     143, 188, 220, 273, 314,
     154, 200, 244, 289, 325,
     171, 221, 270, 326, 358,
     163, 216, 242, 281, 312,
     160, 207, 248, 288, 324,
     142, 187, 234, 280, 316,
     156, 203, 243, 283, 317,
     157, 212, 259, 307, 336,
     152, 203, 246, 286, 321,
     154, 205, 253, 298, 334,
     139, 190, 225, 267, 302,
     146, 191, 229, 272, 302,
     157, 211, 250, 285, 323,
     132, 185, 237, 286, 331,
     160, 207, 257, 303, 345,
     169, 216, 261, 295, 333,
     157, 205, 248, 289, 316,
     137, 180, 219, 258, 291,
     153, 200, 244, 286, 324],
  :x => [8.0, 15.0, 22.0, 29.0, 36.0]
)
rats[:xbar] = mean(rats[:x])
rats[:N] = size(rats[:y], 1)
rats[:T] = size(rats[:y], 2)

rats[:rat] = Int[div(i - 1, 5) + 1 for i in 1:150]
rats[:week] = Int[(i - 1) % 5 + 1 for i in 1:150]
rats[:X] = rats[:x][rats[:week]]
rats[:Xm] = rats[:X] - rats[:xbar]

## Model Specification
model = Model(

  y = Stochastic(1,
    (alpha, beta, rat, Xm, s2_c) ->
      begin
        mu = alpha[rat] + beta[rat] .* Xm
        MvNormal(mu, sqrt(s2_c))
      end,
    false
  ),

  alpha = Stochastic(1,
    (mu_alpha, s2_alpha) -> Normal(mu_alpha, sqrt(s2_alpha)),
    false
  ),

  alpha0 = Logical(
    (mu_alpha, xbar, mu_beta) -> mu_alpha - xbar * mu_beta
  ),

  mu_alpha = Stochastic(
    () -> Normal(0.0, 1000),
    false
  ),

  s2_alpha = Stochastic(
    () -> InverseGamma(0.001, 0.001),
    false
  ),

  beta = Stochastic(1,
    (mu_beta, s2_beta) -> Normal(mu_beta, sqrt(s2_beta)),
    false
  ),

  mu_beta = Stochastic(
    () -> Normal(0.0, 1000)
  ),

  s2_beta = Stochastic(
    () -> InverseGamma(0.001, 0.001),
    false
  ),

  s2_c = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

## Initial Values
inits = [
  Dict(:y => rats[:y], :alpha => fill(250, 30), :beta => fill(6, 30),
       :mu_alpha => 150, :mu_beta => 10, :s2_c => 1, :s2_alpha => 1,
       :s2_beta => 1),
  Dict(:y => rats[:y], :alpha => fill(20, 30), :beta => fill(0.6, 30),
       :mu_alpha => 15, :mu_beta => 1, :s2_c => 10, :s2_alpha => 10,
       :s2_beta => 10)
]

## Sampling Scheme
scheme = [Slice(:s2_c, 10.0),
          AMWG(:alpha, 100.0),
          Slice([:mu_alpha, :s2_alpha], [100.0, 10.0], Univariate),
          AMWG(:beta, 1.0),
          Slice([:mu_beta, :s2_beta], 1.0, Univariate)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, rats, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

"""
Formulation described in:

    Arroyo, J. M., & Conejo, A. J. (2000). Optimal response of a thermal unit
    to an electricity spot market. IEEE Transactions on power systems, 15(3),
    1098-1104. DOI: https://doi.org/10.1109/59.871739
"""
module ArrCon2000

import ..RampingFormulation

struct Ramping <: RampingFormulation end

end


text:
using RecipesBase

@recipe function f(g::Function, cols::Vector{DataColumn}; xtransform=identity, ztransform=identity, contrast=8)
    xs, ys, zs = datamap(g, cols)
    xs, ys = xs[1, :], ys[:, 1]
    xguide := get(cols[1].metadata, :xsymbol, nothing)
    yguide := cols[1].ysymbol
    xs = xtransform.(xs)
    zs = ztransform.(zs)

    xlims --> get(cols[1].metadata, :xlims, (xs[1], xs[end]))
    ylims --> get(cols[1].metadata, :ylims, length(ys) >= 2 ? (ys[1], ys[end]) : nothing)

    fzs = filter(!isnan, zs)
    maxz = maximum(abs.(fzs))
    stdz = std(fzs)
    clim = max(min(maxz, contrast*stdz), maxz / contrast)

    if all(x -> x >= 0, zs)
        clims --> (0.0, clim)
        seriescolor --> :bilbao
    else
        clims --> (-clim, clim)
        seriescolor --> :vik
    end

    # @series begin
    #     linewidth --> 1
    #     linecolor --> :black
    #     seriesalpha --> 0.3
    #     seriestype := :path
    #     primary := false
    #     markershape := :none

    #     (xs, mapreduce(col -> col.metadata.sim.eig.values, hcat, cols)')
    # end

    (xs, ys, zs)
end

@recipe function f(cols::Vector{DataColumn}; xtransform=identity, markerfunction=x->:black, step=1)
    xs, ys = ranges(cols)
    xguide --> get(cols[1].metadata, :xsymbol, nothing)
    xs = xtransform.(xs)
    eigvalues, energy_range = energy_bands(cols, step=step)

    xlims --> get(cols[1].metadata, :xlims, (xs[1], xs[end]))
    ylims --> get(cols[1].metadata, :ylims, length(ys) >= 2 ? (ys[1], ys[end]) : nothing)

    yticks --> energy_range

    linewidth --> 1
    linecolor --> :black
    seriesalpha --> 0.3
    legend --> false
    # markercolor --> mapreduce(col -> map(markerfunction, col.data), hcat, cols)
    # seriestype := :path
    # markershape := :none

    @series begin
        μs = mapreduce(col -> col.metadata.sim.μ, vcat, cols)
        linewidth --> 1
        linecolor := :red
        linestyle --> :dash
        (xs, μs)
    end

    (xs, eigvalues)
end

@recipe function f(g::Function, slices::Vector{DataSlice}; xtransform=identity)
    xs, ys = datamap(g, slices)
    xs = xtransform.(xs)
    xlims --> (xs[1], xs[end])
    ylim = maximum(abs.(ys))
    if all(y -> y >= 0, ys)
        ylims --> (0.0, ylim)
    else
        ylims --> (-ylim, ylim)
    end

    #linewidth --> 1

    (xs, ys)
end

@recipe function f(col::DataColumn; xtransform=identity, markerfunction=x->:black)
    xs = col.metadata.sim.eig.values
    xticks --> range(floor(Int, minimum(xs)), ceil(Int, maximum(xs)), step=1)

    xlims --> get(col.metadata, :xlims, (xs[1], xs[end]))

    linewidth --> 1
    linecolor --> :black
    legend --> false

    @series begin
        μ = col.metadata.sim.μ
        linewidth --> 1
        linecolor := :red
        linestyle --> :dash
        [μ]
    end

    xs
end

@recipe function f(g::Function, col::DataColumn; xtransform=identity, markerfunction=x->:black)
    xs = col.ys
    #xticks --> range(floor(Int, minimum(xs)), ceil(Int, maximum(xs)), step=1)

    ys = map(g, col.data)
    xguide := col.ysymbol

    xlims --> get(col.metadata, :ylims, length(ys) >= 2 ? (xs[1], xs[end]) : nothing)

    linewidth --> 1
    linecolor --> :black
    legend --> false

    xs, ys
end

@recipe function f(T::Type{DataScatter}, g::Function, cols::Vector{DataColumn}; xtransform=identity, ztransform=identity, contrast=8)
    xs, ys, zs = T(g, cols)
    xguide := get(cols[1].metadata, :xsymbol, nothing)
    yguide := cols[1].ysymbol
    xs = xtransform.(xs)
    zs = ztransform.(zs)

    xlims --> get(cols[1].metadata, :xlims, (xs[1], xs[end]))
    ylims --> get(cols[1].metadata, :ylims, length(ys) >= 2 ? (ys[1], ys[end]) : nothing)

    fzs = filter(!isnan, zs)
    maxz = maximum(abs.(fzs))
    stdz = std(fzs)
    clim = max(min(maxz, contrast*stdz), maxz / contrast)

    if all(x -> x >= 0, zs)
        clims --> (0.0, clim)
        color --> :bilbao
    else
        clims --> (-clim, clim)
        color --> :vik
    end

    zcolor --> zs
    ms --> max.(1.0, 20 .* sqrt.(abs.(zs)))
    α --> sqrt(0.5)

    (xs, ys)
end


text:
"""
    Hopper <: Environment

    planar hopping robot, based on https://gym.openai.com/envs/Hopper-v2/
"""
struct Hopper end

function hopper(;
    representation=:minimal,
    timestep=0.05,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=1.0,
    spring=10.0,
    damper=50.0,
    seed=1,
    contact_foot=true,
    contact_body=true,
    info=nothing,
    vis=Visualizer(),
    name=:robot,
    opts_step=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    opts_grad=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    T=Float64)

    mechanism = get_hopper(
        timestep=timestep,
        gravity=gravity,
        friction_coefficient=friction_coefficient,
        spring=spring,
        damper=damper,
        contact_foot=contact_foot,
        contact_body=contact_body)

    initialize_hopper!(mechanism)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
    end
    nu = 3
    no = nx - 1 # full_state is not being returned by default

    # values taken from Mujoco's model, combining the control range -1, 1 and the motor gears.
    aspace = BoxSpace(nu,
        low=(-Inf * ones(nu)),
        high=(Inf * ones(nu)))
    ospace = BoxSpace(no,
        low=(-Inf * ones(no)),
        high=(Inf * ones(no)))

    rng = MersenneTwister(seed)

    z = get_maximal_state(mechanism)
    x = representation == :minimal ? maximal_to_minimal(mechanism, z) : z

    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = zeros(nu)
    control_mask = [zeros(nu, 3) I(nu)]
    motor_gear = [200, 200, 200.]
    control_scaling = Diagonal(timestep * motor_gear)

    build_robot(mechanism, vis=vis, name=name)

    TYPES = [Hopper, T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    env = Environment{TYPES...}(mechanism, representation, aspace, ospace,
        x, fx, fu,
        u_prev, control_mask' * control_scaling,
        nx, nu, no,
        info,
        [rng], vis,
        opts_step, opts_grad)

    return env
end

function Base.reset(env::Environment{Hopper};
    x=nothing, reset_noise_scale=0.005)

    if x != nothing
        env.state .= x
    else
        # initialize above the ground to make sure that with random initialization we do not violate the ground constraint.
        initialize!(env.mechanism, :hopper,
            z=0.25)
        x0 = get_minimal_state(env.mechanism)
        nx = minimal_dimension(env.mechanism)

        low = -reset_noise_scale
        high = reset_noise_scale
        x = x0 + (high - low) .* rand(env.rng[1], nx) .+ low # we ignored the normal distribution on the velocities
        z = minimal_to_maximal(env.mechanism, x)
        set_maximal_state!(env.mechanism, z)
        if env.representation == :minimal
            env.state .= get_minimal_state(env.mechanism)
        elseif env.representation == :maximal
            env.state .= get_maximal_state(env.mechanism)
        end
        env.input_previous .= 0.0
    end
    return get_observation(env)
end

function get_observation(env::Environment{Hopper};
    full_state=false)

    full_state && (return env.state)

    nx = minimal_dimension(env.mechanism)
    if env.representation == :minimal
        o = env.state
    elseif env.representation == :maximal
        o = maximal_to_minimal(env.mechanism, env.state)
    end

    # clamp velocities and remove x position
    ind = velocity_index(env.mechanism)
    vel = clamp.(o[ind], -10, 10)
    o[ind] .= vel
    o = o[[1; 3:nx]]
    return o
end

function cost(env::Environment{Hopper}, x, u;
        alive_bonus=0.1)

    if env.representation == :minimal
        x_velocity = -x[5]
    else
        i_torso = findfirst(body -> body.name == "torso", collect(env.mechanism.bodies))
        z_torso = x[(i_torso-1)*13 .+ (1:13)]
        x_velocity = -z_torso[4]
    end
    # c = -x_velocity/10 # Gym reward
    c = -x_velocity
    c -= alive_bonus
    c += 1e-4 * u'*u
    return c
end

function is_done(env::Environment{Hopper}, x)
    nx = minimal_dimension(env.mechanism)
    if env.representation == :minimal
        x0 = x
    elseif env.representation == :maximal
        x0 = maximal_to_minimal(env.mechanism, x)
    end
    height = x0[1]
    ang = x0[3]
    done = !(
        all(isfinite.(x0)) &&
        all(abs.(x0[[1;3:nx]]) .< 100) &&
        (height > 0.7) &&
        (abs(ang) < 0.2)
        )
    return done
end


text:
using Revise
using PyCall
using PyPlot

function phicoef(x, y, a)
    mat = zeros(length(x), 3)
    mat[:, 1] = @. (x / a) * (9 * (x / a) / 8 - 3 / 4)
    mat[:, 2] = @. (1 - 3 * (x / a) / 2) * (1 + 3 * (x / a) / 2)
    mat[:, 3] = @. (x / a) * (9 * (x / a) / 8 + 3 / 4)
    return inv(mat) * y
end

function fig_3nq2linear()
    linewidth = 0.5
    fontsize = 6
    n = 1000
    a = 1
    x = LinRange(-a, a, n)

    # ϕ shape functions and coefficients
    phi1 = @. (x / a) * (9 * (x / a) / 8 - 3 / 4)
    phi2 = @. (1 - 3 * (x / a) / 2) * (1 + 3 * (x / a) / 2)
    phi3 = @. (x / a) * (9 * (x / a) / 8 + 3 / 4)
    xvals = [-1.0 ; 0.0 ; 1.0]
    yvals = @. [-1.0 ; 0.5 ; 2.0] + 1.25 # linear slip
    coef = phicoef(xvals, yvals, a)

    close("all")
    figure(figsize=(3, 3))
    plot(x, coef[1] * phi1, "--k", color="blue", linewidth=linewidth)
    plot(x, coef[2] * phi2, "--k", color="blue", linewidth=linewidth)
    plot(x, coef[3] * phi3, "--k", color="blue", linewidth=linewidth)
    plot(x, coef[1] * phi1 + coef[2] * phi2 + coef[3] * phi3, "-r", linewidth=linewidth)
    xticks([-1, 0, 1])
    yticks([-4, 0, 4])
    xlim([-1, 1])
    ylim([-4, 4])
    xlabel(L"$x$ (m)", fontsize=fontsize)
    ylabel(L"$u$ (m)", fontsize=fontsize)
    gca().tick_params(labelsize = fontsize)
    show()
end
fig_3nq2linear()


text:
function laplace_approximation(K::PDMats.PDMat,
                               choices::Array{<:Int, 2},
                               initial_latent::Vector{<:Real},
                               scale::Float64;
                               verbose::Bool=true)
#=
    Variant of the Newton's method based mode-locating algorithm (GPML, Algorithm 3.1)
    Utilizes the Woodburry identity for avoiding two cholesky factorizations
    per Newton iteration.

    Reduces the stepsize whenever the marginal likelhood gets stuck
    Algortihm 3.1 utilizes the fact that W is diagonal which is not for our case.

    Note: ( K^{-1} + W )^{-1} = K ( I - ( I + W K )^{-1} ) W K
=##
    max_iter = 20

    f            = initial_latent
    prev_f       = deepcopy(f)
    prev_mll     = -Inf

    WK  = nothing
    α   = 1.0
    a   = nothing
    B   = nothing
    for iteration = 1:max_iter
        btl_mat = logbtl_matrix(choices, f, scale)
        ∇ll     = ∇logbtl(btl_mat, f, choices, scale)
        W       = -∇²logbtl(btl_mat, ∇ll, choices, scale)

        WK  = W*K
        b   = W*f + ∇ll
        B   = I + WK
        Blu = lu(B)
        a   = (b - Blu \ (WK*b))

        f   = α*(K*a) + (1 - α)*prev_f
        mll = sum(btl_mat[:,1]) + dot(a, f)/-2 + logdet(Blu)/-2

        ∇mll = norm(∇ll - a)
        Δf   = norm(prev_f - f)

        if(∇mll < 1e-4 || Δf < 1e-4)
            break
        end

        if(mll <= prev_mll)
            α /= 2
        end

        if(verbose)
            @info "Laplace approximation stat" iteration ∇mll Δf
        end

        prev_f   = deepcopy(f)
        prev_mll = mll
    end
    μ  = f
    Σ  = K * (I - (I + WK) \ WK)
    # Due to numerical inaccuracies,
    # Σ often is not symmetric out of the box
    Σu = UpperTriangular(Σ)
    Σd = Diagonal(Σ)
    Σ  = Σu + Σu' - Σd
    return μ, Σ, a, B
end


text:
export clipdb, fold, mps, grpdelay;
include("./clipdb.jl");
include("./fold.jl");
include("./mps.jl");
include("./grpdelay.jl");


text:
using Test
using FlightMechanics
using FlightMechanics.Models
using FlightMechanics.Aircrafts

#  Stevens, B. L., Lewis, F. L., & Johnson, E. N. (2015). Aircraft control
#  and simulation: dynamics, controls design, and autonomous systems. John Wiley
#  & Sons. (page 205)
pos = PositionEarth(0, 0, 0)
env = Environment(pos, atmos="ISAF16", wind="NoWind", grav="const")

FlightMechanics.Aircrafts.get_empty_cg(ac::F16) = [-0.30 * get_chord(ac), 0, 0]

ac = F16()
fcs = F16FCS()

h = 0.0 * M2FT
psi = 0.234077 # rad
gamma = 0.0
turn_rate = 0.3  # rad/s
tas_fts = 502.0  # ft/s

tas =  tas_fts * FT2M

# Initial conditions for trimmer
exp_α    =  2.392628e-1  # rad
exp_β    =  5.061803e-4  # rad
exp_thtl =  0.8349601
exp_de   = -1.481766     # deg
exp_da   =  9.553108e-2  # deg
exp_dr   = -4.118124e-1  # deg
α0 = exp_α + 10*DEG2RAD
β0 = 0.
stick_lon0 = exp_de/25.0 + 0.5
stick_lat0 = exp_da/21.5 + 0.5
pedals0 = exp_dr/30.0 + 0.5
thtl0 = exp_thtl + 0.3

set_stick_lon(fcs, stick_lon0)
set_stick_lat(fcs, stick_lat0)
set_pedals(fcs, pedals0)
set_thtl(fcs, thtl0)

ac, aerostate, state, fcs = steady_state_trim(
    ac, fcs, env, tas, pos, psi, gamma, turn_rate, α0, β0, show_trace=false
);

linear_system = Linearize(ac, state, aerostate, fcs, env);

A = get_A(linear_system)
A_ = [
    # u      v      theta  q      w      phi    p      r
    A[1,1] A[1,2] A[1,8] A[1,5] A[1,3] A[1,9] A[1,4] A[1,5];
    A[2,1] A[2,2] A[2,8] A[2,5] A[2,3] A[2,9] A[2,4] A[2,5];
    A[8,1] A[8,2] A[8,8] A[8,5] A[8,3] A[8,9] A[8,4] A[8,5];
    A[5,1] A[5,2] A[5,8] A[5,5] A[5,3] A[5,9] A[5,4] A[5,5];
    A[3,1] A[3,2] A[3,8] A[3,5] A[3,3] A[3,9] A[3,4] A[3,5];
    A[9,1] A[9,2] A[9,8] A[9,5] A[9,3] A[9,9] A[9,4] A[9,5];
    A[4,1] A[4,2] A[4,8] A[4,5] A[4,3] A[4,9] A[4,4] A[4,5];
    A[5,1] A[5,2] A[5,8] A[5,5] A[5,3] A[5,9] A[5,4] A[5,5]
]

# from book
A_exp = [
    -.090 -169 -31.2 -7.75 31.4 -7.73 5e-4 2e-3;
    -5e-4 -1.05 .0151 .903 3e-4 -.0607 -5e-4 -1e-4;
    .0 .0 .0 .203 .0 -.300 .0 -.979;
    1e-3 1.26 0 -1.66 1e-3 .0 .0589 -.0157;
    -1e-4 1.3e-4 -0.0032 7e-6 -.322 .0130 .248 -.961;
    .0 .0 .300 .508 .0 .0 1.0 .0105;
    -3e-4 .0578 .0 -.0469 -59.4 .0 -3.19 1.64;
    5e-5 -.0617 .0 .0123 8.88 .0 -.299 -.564
]


text:
struct FISH <: AbstractShell
    path::AbstractString
    prompt::AbstractString
end

FISH(path::AbstractString) = FISH(path, "\\e[0;35m\\u@\\h:\\W (playground)> \\e[m")
FISH() = FISH(strip(readstring(`which fish`)))

function Base.run(shell::FISH, env::Environment)
    init(env)
    # Currently can't support sourcing a custom fish prompt on startup until
    # 2.7.0 is released with the `-C` flag.
    # https://github.com/fish-shell/fish-shell/milestone/17
    @mock run(`$(shell.path) -i`)
end


text:
module TestTaskGroups

using Julio
using Test

function test_trivial_spawns()
    ref1 = Ref(0)
    ref2 = Ref(0)
    Julio.withtaskgroup() do tg
        Julio.spawn!(tg) do
            ref1[] = 111
        end
        Julio.spawn!(tg) do
            ref2[] = 222
        end
    end
    @test (ref1[], ref2[]) == (111, 222)
end

function test_cancel_with_single_main_error()
    err_cause = ErrorException("error in main")
    err = try
        Julio.withtaskgroup() do tg
            _ie, oe = Julio.channel()
            Julio.spawn!(tg) do
                take!(oe)  # should be unblocked by cancellation
            end
            throw(err_cause)
        end
        nothing
    catch err
        err
    end
    @test err == err_cause
end

end  # module


text:
export modifiedMatsuno

"""
    modifiedMatsuno(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the modified Matsuno model

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- numerical parameters to construct the mesh of collocation points as `mesh = Mesh(param)`

## Optional keyword arguments
- `ν`: shallow/deep water multiplication factor. By default, `ν=1` if `μ≦1` and `ν=1/√μ` otherwise. Set the infinite-layer case if `ν=0` (or `μ=Inf`).
- `IL`: Set the infinite-layer case if `IL=true` (or `μ=Inf`, or `ν=0`), in which case `ϵ` is the steepness parameter. Default is `false`.
- `ktol`: tolerance of the low-pass Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `label`: a label for future references (default is `"modified Matsuno"`);
- `verbose`: prints information if `true` (default is `true`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `modifiedMatsuno.f!` to be called in explicit time-integration solvers;
2. a function `modifiedMatsuno.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `modifiedMatsuno.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v)`, where
    - `η` is the surface deformation;
    - `v` is the derivative of the trace of the velocity potential.

"""
mutable struct modifiedMatsuno <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function

    function modifiedMatsuno(param::NamedTuple;
							ν		= nothing,
							IL	    = false,
							ktol	= 0,
							dealias	= 0,
							label	= "modified Matsuno",
							verbose	= true )

		# Set up
		μ 	= param.μ
		ϵ 	= param.ϵ
		if ν == nothing
			if μ > 1
				ν = 1/sqrt(μ)
				nu = "1/√μ (deep water case)"
			else
				ν = 1
				nu = "1 (shallow water case)"
			end
		else
			nu = "$ν"
		end
		if μ == Inf || ν==0 # infinite layer case
			IL = true;  # IL (=Infinite layer) is a flag to be used therefafter
			μ = 1; ν = 1; # Then we should set μ=ν=1 in subsequent formula.
		else IL = false
		end

		if verbose  # Print information
			info = "Build the modified 'exponential' Matsuno model of Duchêne and Melinand.\n"
			if IL == true
				info *= "Steepness parameter ϵ=$ϵ (infinite depth case).\n"
			else
				info *= "Shallowness parameter μ=$μ, nonlinearity parameter ϵ=$ϵ, \
						scaling parameter ν=$nu.\n"
			end
			if dealias == true || dealias == 1
				info *= "Dealiasing with Orszag’s 3/2 rule. "
			else
				info *= "No dealiasing. "
			end
			if ktol == 0
				info *= "No Krasny filter. "
			else
				info *= "Krasny filter with tolerance $ktol."
			end
			info *= "\nShut me up with keyword argument `verbose = false`."
			@info info
			@warn "The velocity is consistent with the \
			derivative of the trace of the velocity potential \
			used for water waves only when they are null."
		end

		# Pre-allocate useful data
		mesh = Mesh(param)
		x = mesh.x
		k = mesh.k
		if dealias == 0
			Π⅔ 	= ones(size(k)) # no dealiasing (Π⅔=Id)
		else
			K = mesh.kmax * (1-dealias/(2+dealias))
			Π⅔ 	= abs.(k) .<= K # Dealiasing low-pass filter
		end
		if IL == true
			∂ₓF₀ 	= 1im * sign.(k)
			G₀ 	= abs.(k)
		else
        	∂ₓF₀ 	= 1im* sign.(k) .* tanh.(sqrt(μ)*abs.(k))
			G₀ 	= sqrt(μ)*abs.(k).*tanh.(sqrt(μ)*abs.(k))
		end
		∂ₓ	=  1im * sqrt(μ)* k            # Differentiation
		z = zeros(Complex{Float64}, mesh.N)
		η = copy(z) ; v = copy(z) ; F = copy(z) ;
		fftη = copy(z) ; fftv = copy(z) ; Q = copy(z) ; R = copy(z) ;

		# Build raw data from physical data.
		# Discrete Fourier transform with, possibly, dealiasing and Krasny filter.
		# This function is correct only when the initial data for v is zero.
		function mapto(data::InitialData)
			U = [Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]
			U[abs.(U).< ktol ].=0
			return U
		end

		# Return `(η,v)`, where
		# - `η` is the surface deformation;
		# - `v` is the velocity variable.
		# Inverse Fourier transform and takes the real part.
		function mapfro(U)
			real(ifft(U[:,1])),real(ifft(U[:,2]))
		end

		# Evolution equations are ∂t U = f(U)
		function f!(U)

			fftη .= U[:,1]
		    fftv .= U[:,2]
			η  .= ifft(U[:,1])
			v  .= ifft(U[:,2])
			Q .= -∂ₓF₀.*fftv
			Q += -ϵ *∂ₓ.*fft(η.*ifft(fftv)) .+ ϵ*∂ₓF₀.*fft(η.*ifft(G₀.*fftv))
			F .= exp.(-ϵ*ifft(G₀.*fftη))
			R .= -fft(F.*ifft(∂ₓ.*fftη))*ν
			R += ϵ/2*∂ₓ.*fft(-v.^2)

			U[:,1] .= Π⅔.*Q/sqrt(μ)/ν
		   	U[:,2] .= Π⅔.*R/sqrt(μ)/ν
			U[abs.(U).< ktol ].=0

		end

        new(label, f!, mapto, mapfro )
    end
end


text:
module MiscUtil

    using LinearAlgebra, Printf, Crayons, PyPlot

    import Distributions: AbstractMvNormal

    import Statistics: mean, cov

    # export calculaterbfbasis, calculaterbfbasis!

    export safewrapper

    export addjitter, addjitter!

    export makematrixsymmetric, makematrixsymmetric!

    export nearestposdef

    export woodbury_A_plus_UCV

    export tasksleep

    export sigmoid, invsigmoid, makepositive, invmakepositive, transformbetween, invtransformbetween

    export logrange, LogRange

    export bracket

    export colourprint

    export waitforkey

    export getcolourrange

    export plot_ellipse

    include("transformutil.jl")

    include("matricesutil.jl")

    include("sleeputil.jl")

    include("safewrapper.jl")

    include("logrange.jl")

    include("bracket.jl")

    include("colourprint.jl")

    include("waitforkey.jl")

    include("getcolourrange.jl")

    include("plot_ellipse.jl")

end


text:
module MusicProcessingTest

using MusicProcessing
using Test

tests = [
    "audio.jl",
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end

end


text:
## Default node and edge style definitions ##

constattr(v) = ConstantAttribute(v, ()->1)

default_node_style() = AttributeDict(
        :shape          => constattr(:circle),
        :size           => constattr(100),
        :color          => constattr((.7, .1, .2)),
        :border_color   => constattr((1., 1., 1.)),
        :border_width   => constattr(.5),
        :opacity        => constattr(.8),
        :label          => constattr(""),
        :label_color    => constattr((0., 0., 0.)))

default_edge_style() = AttributeDict(
        :shape          => constattr(:arrow),
        :width          => constattr(.5),
        :color          => constattr((.5, .5, .5)),
        :opacity        => constattr(.8),
        :curve          => constattr(0))

function setup_node_style(g::Graph; kvargs...)
    style = default_node_style()
    for (k, v) in g.nodeattrs       # incorporate node attributes
        if haskey(style, k)
            style[k] = v
        end
    end
    for (k, v) in kvargs            # incorporate node kvargs
        ks = String(k)
        if startswith(ks, "node_")
            k_ = Symbol(ks[6:end])
            if haskey(style, k_)
                style[k_] = attribute(v, ()->1)
            end
        end
    end
    return style
end

function setup_edge_style(g::Graph; kvargs...)
    style = default_edge_style()
    for (k, v) in g.edgeattrs       # incorporate edge attributes
        if haskey(style, k)
            style[k] = v
        end
    end
    for (k, v) in kvargs            # incorporate edge kvargs
        ks = String(k)
        if startswith(ks, "edge_")
            k_ = Symbol(ks[6:end])
            if haskey(style, k_)
                style[k_] = attribute(v, ()->1)
            end
        end
    end
    return style
end


text:
"""
Peeks `n`(can be positive or negative) steps from `pos` into list.
Returns nothing if out of bounds and jumps over nothing and LineNumberNode's
"""
function peek_n(n, list, pos)
    i = pos; step = sign(n)
    while abs(n) > 0
        (i<1 || i>length(list)) && return nothing
        elem = list[i]
        if !(isa(elem, LineNumberNode) || isa(elem, Void))
            n -= step
        end
        i += step
    end
    list[i]
end

"""
Determines if `elem` is a `unless`/`gotoifnot` node
"""
isa_unless(elem) = isa(elem, Expr) && elem.head == :gotoifnot

"""
Replaces `goto` statements in a loop body with continue and break.
"""
function replace_continue_break!(list, continue_label, break_label)
    map!(list) do elem
        if isa(elem, GotoNode) && elem.label == continue_label.label
            Expr(:continue)
        elseif isa(elem, GotoNode) && elem.label == break_label.label
            Expr(:break)
        else
            elem
        end
    end
end

"""
Removes all gotos, labels and nothing's inplace recursively in `Expr` `expr`.
"""
remove_labelgoto!(expr::Expr) = (_remove_labelgoto!(expr); nothing)
_remove_labelgoto!(x) = true
_remove_labelgoto!(x::Union{LabelNode, GotoNode, Void}) = false
function _remove_labelgoto!(x::Expr)
    filter!(_remove_labelgoto!, x.args)
    true
end

"""
Checks if nodes at `startpos` and `labelpos` define a loop body.
If so it will build up the matching loop expression and
`return true, LoopExpr, labelpos`.
Else it'll `return false, Expr(:block), labelpos`
"""
function handle_loop(list, startpos, labelpos, body)
    unless, breaklabel = list[startpos], list[labelpos]
    @assert unless.args[2] == breaklabel.label # we should have made this sure already
    looplabel = peek_n(-1, list, startpos)
    isa(looplabel, LabelNode) || return false, Expr(:block), labelpos
    goto = peek_n(-1, list, labelpos)
    isa(goto, GotoNode) || return false, Expr(:block), labelpos
    looplabel.label == goto.label || return false, Expr(:block), labelpos
    continuelabel = peek_n(-2, list, labelpos)
    isa(continuelabel, LabelNode) || return false, Expr(:block), labelpos
    # now that we have continue/break labels, we can replace the gotos
    replace_continue_break!(body.args, continuelabel, breaklabel)
    x = _remove_goto(body.args)
    true, Expr(:while, unless.args[1], x), labelpos+1
end

"""
Checks if nodes at `startpos` and `labelpos` define an if body.
If so it will build up the matching if(else) expression and
`return true, LoopExpr, labelpos`.
Else it'll `return false, Expr(:block), labelpos`
"""
function handle_if(list, startpos, labelpos, ifbody)
    unless = list[startpos]
    if_expr = Expr(:if, unless.args[1], _remove_goto(ifbody.args))
    goto = peek_n(-1, list, labelpos)
    if isa(goto, GotoNode) # that means we're dealing with else
        pos = labelpos; else_body = Expr(:block)
        for i=pos:length(list)
            elem = list[i]
            push!(else_body.args, elem)
            if isa(elem, LabelNode) && elem.label == goto.label # found else end
                push!(if_expr.args, _remove_goto(else_body.args))
                return true, if_expr, i
            end
        end
    else
        return true, if_expr, labelpos+1
    end
    error("If with goto found (else), but goto label not found. Unless: $unless \n ifbody: $ifbody")
end
"""
Handles a branch (`unless`/`gotoifnot`) at `unless_idx` in `list`.
It will return the corresponding controle flow `Expr` (e.g. `while`/`if`/`else`) and the
new position in the list.
"""
function handle_unless(list, unless_idx)
    unless_body = Expr(:block)
    unless_goto = list[unless_idx].args[2]
    pos = unless_idx + 1
    for i=pos:length(list)
        elem = list[i]
        if isa(elem, LabelNode) && elem.label == unless_goto
            isloop, loopexpr, i = handle_loop(list, pos-1, i, unless_body)
            isloop && return loopexpr, i
            isif, ifexpr, i = handle_if(list, pos-1, i, unless_body)
            isif && return ifexpr, i
        else
            push!(unless_body.args, elem)
        end
    end
end
"""
Replaces goto's with control flow statements while leaving labels/gotos in the AST.
This makes the different steps of processing easier.
"""
function _remove_goto(list)
    # (╯°□°）╯︵ ┻━┻
    result = Expr(:block); i = 1
    while i <= length(list)
        elem = list[i]
        expr = if isa_unless(elem)
            expr, i = handle_unless(list, i)
            expr
        else
            i += 1; elem
        end
        push!(result.args, expr)
    end
    result
end
"""
Returns a new expression block with all goto statements replaced by the corresponding
control flow `Expr` (e.g. `while`/`if`/`else`).
"""
function remove_goto(list)
    expr_nogoto = _remove_goto(list)
    remove_labelgoto!(expr_nogoto)
    expr_nogoto
end


text:
for frange in (mrange, srange)
    @testset "intersect-$frange" begin
        @test intersect(frange(1, 5), frange(2, 3)) == 2:3
        @test intersect(frange(-3, 5), frange(2, 8)) == 2:5
        @test intersect(frange(-8, -3), frange(-8, -3)) == -8:-3
        @test intersect(frange(1, 5), frange(5, 13)) == 5:5
        @test isempty(intersect(frange(-8, -3), frange(-2, 2)))

        @test isempty(intersect(frange(-3, 7), frange(2, 1)))
        @test intersect(frange(1, 11), frange(-2, step=3, stop=15)) == 1:3:10
        @test intersect(frange(1, 11), frange(-2, step=2, stop=15)) == 2:2:10
        @test intersect(frange(1, 11), frange(-2, step=1, stop=15)) == 1:11
        @test intersect(frange(1, 11), frange(15, step=-1, stop=-2)) == 1:11
        @test intersect(frange(1, 11), frange(15, step=-4, stop=-2)) == 3:4:11
        @test intersect(frange(-20, -5), frange(-10, step=3, stop=-2)) == -10:3:-7
        @test isempty(intersect(frange(-5, 5), frange(-6, step=13, stop=20)))
        @test isempty(intersect(frange(1, 11), frange(15, step=4, stop=-2)))
        @test isempty(intersect(frange(11, 1), frange(15, step=-4, stop=-2)))
        #@test intersect(-5:5, 1+0*(1:3)) == 1:1
        #@test isempty(intersect(-5:5, 6+0*(1:3)))
        @test intersect(frange(-15, step=4, stop=7), frange(-10, -2)) == -7:4:-3
        @test intersect(frange(13, step=-2, stop=1), frange(-2, 8)) == 7:-2:1
        @test isempty(intersect(frange(13, step=2, stop=1), frange(-2, 8)))
        @test isempty(intersect(frange(13, step=-2, stop=1), frange(8, -2)))
        #@test intersect(5+0*(1:4), 2:8) == 5+0*(1:4)
        #@test isempty(intersect(5+0*(1:4), -7:3))
        @test intersect(frange(0, step=3, stop=24), frange(0, step=4, stop=24)) == 0:12:24
        @test intersect(frange(0, step=4, stop=24), frange(0, step=3, stop=24)) == 0:12:24
        @test intersect(frange(0, step=3, stop=24), frange(24, step=-4, stop=0)) == 0:12:24
        @test intersect(frange(24, step=-3, stop=0), frange(0, step=4, stop=24)) == 24:-12:0
        @test intersect(frange(24, step=-3, stop=0), frange(24, step=-4, stop=0)) == 24:-12:0
        @test intersect(frange(1, step=3, stop=24), frange(0, step=4, stop=24)) == 4:12:16
        @test intersect(frange(0, step=6, stop=24), frange(0, step=4, stop=24)) == 0:12:24
        @test isempty(intersect(frange(1, step=6, stop=2400), frange(0, step=4, stop=2400)))
        @test intersect(frange(-51, step=5, stop=100), frange(-33, step=7, stop=125)) == -26:35:79
        @test intersect(frange(-51, step=5, stop=100), frange(-32, step=7, stop=125)) == -11:35:94
        #@test intersect(0:6:24, 6+0*(0:4:24)) == 6:6:6
        #@test intersect(12+0*(0:6:24), 0:4:24) == AbstractRange(12, 0, 5)
        #@test isempty(intersect(6+0*(0:6:24), 0:4:24))
        @test intersect(frange(2, step=1, stop=1), frange(2, step=1, stop=1)) == intersect(range(2, step=1, stop=1), range(2, step=1, stop=1))
        @test intersect(frange(-10, step=3, stop=24), frange(-10, step=3, stop=24)) == -10:3:23
        @test isempty(intersect(frange(-11, step=3, stop=24), frange(-10, step=3, stop=24)))
        @test intersect(frange(typemin(Int), step=2, stop=typemax(Int)), 1:10) == 2:2:10
        @test intersect(1:10, frange(typemin(Int), step=2, stop=typemax(Int))) == 2:2:10

        @test intersect(reverse(typemin(Int):2:typemax(Int)), frange(typemin(Int), step=2, stop=typemax(Int))) == reverse(typemin(Int):2:typemax(Int))
        @test intersect(typemin(Int):2:typemax(Int), reverse(frange(typemin(Int), step=2, stop=typemax(Int)))) == typemin(Int):2:typemax(Int)

#            @test intersect(UnitRange(1,2),3) == UnitRange(3,2)
#            @test intersect(UnitRange(1,2), UnitRange(1,5), UnitRange(3,7), UnitRange(4,6)) == UnitRange(4,3)

        @test intersect(frange(1, 3), 2) == intersect(2, frange(1, 3)) == frange(2, 2)
        @test intersect(frange(1.0, 3.0), 2) == intersect(2, frange(1.0, 3.0)) == [2.0]

        if VERSION > v"1.2"
            @testset "Support StepRange with a non-numeric step" begin
                start = Date(1914, 7, 28)
                stop = Date(1918, 11, 11)

                @test intersect(frange(start, step=Day(1), stop=stop), start:Day(1):stop) == start:Day(1):stop
                @test intersect(start:Day(1):stop, start:Day(5):stop) == start:Day(5):stop
                @test intersect(start-Day(10):Day(1):stop-Day(10), start:Day(5):stop) ==
                                start:Day(5):stop-Day(10)-mod(stop-start, Day(5))
            end
        end
    end
end


text:
module day10
using StatsBase # countmap
using Pipe: @pipe

include("inputs.jl")

ex0 = inputs.ex0_d10()
ex1 = inputs.ex1_d10()
input = inputs.input_d10()

function part1(input)
    adapter = sort(input)
    adapter = vcat(0, adapter, adapter[end]+3)
    diff_counts = countmap(diff(adapter))
    return diff_counts[1] * diff_counts[3]
end

@assert part1(ex1) == 220
@assert part1(input) == 2176
@info "Day 10, Part 1 answer: $(part1(input))"

# this one sucked. alg dev by trial and error.
function part2(input)
    adapter = sort(input)
    adapter = vcat(0, adapter, adapter[end]+3)
    diffs = diff(adapter) #

    arrangements = @pipe join(diffs) |>
        split(_, '3', keepempty=false) |> # group the "1's" in arrays
        length.(_) |> # count the number of 1's
        binomial.(_, 2) .+ 1 |> # arrangment is num of 1's-choose-2 (plus the 3's)
        reduce(*, _) # all multiplied together. phew.
end

@assert part2(ex0) == 8
@assert part2(ex1) == 19208
@assert part2(input) == 18512297918464
@info "Day 10, Part 2 answer: $(part2(input))"
end


text:
# Autogenerated wrapper script for ColPack_jll for x86_64-w64-mingw32-cxx11
export libcolpack

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ColPack")
JLLWrappers.@declare_library_product(libcolpack, "libColPack-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libcolpack,
        "bin\\libColPack-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()