text:
import Distributions.sample
import Base.append!
using Flux: Zygote.@adjoint

function test(f, env, maxt=100)
    R = 0f0
    s = reset!(env)
    for _ in 1:maxt
        s, r, done, _ = step!(env, f(s))
        R += r
        done && break
    end
    return R
end

gather(y, a) = y[CartesianIndex.(a, axes(a, 1))]
gather(y, a::AbstractMatrix) = gather(y, vec(a))

@adjoint gather(y, a) = gather(y, a),
    ȳ -> begin
        o = zeros(eltype(ȳ), size(y))
        @inbounds for (i, j) in enumerate(a)
            o[j, i] = ȳ[i]
        end
        return (o, nothing)
    end

function softupdate!(target::T, model::T, τ=1f-2) where T
    for f in fieldnames(T)
        softupdate!(getfield(target, f), getfield(model, f), τ)
    end
end

function softupdate!(dst::A, src::A, τ=T(1f-2)) where {T, A<:AbstractArray{T}}
    dst .= τ .* src .+ (one(T) - τ) .* dst
end

mutable struct ReplayMemory{S<:Real, A<:Real, R<:Real}
    n::Int64
    maxsize::Int64
    batchsize::Int64
    s::Matrix{S}
    a::Matrix{A}
    r::Vector{R}
    d::Vector{Bool}
    function ReplayMemory{S, A, R}(obssize, actiondim, maxsize, batchsize) where {S<:Real, A<:Real, R<:Real}
        new(
            0,
            maxsize,
            batchsize,
            zeros(S, obssize, maxsize),
            zeros(A, actiondim, maxsize),
            zeros(R, maxsize),
            zeros(Bool, maxsize)
        )
    end
end

function append!(memory::ReplayMemory, s, a, r, d)
    ind = memory.n % memory.maxsize + 1
    memory.n += 1
    memory.s[:, ind] .= s
    memory.a[:, ind] .= a
    memory.r[ind] = r
    memory.d[ind] = d
end

function sample(memory::ReplayMemory)
    inds = sample(1:(min(memory.n, memory.maxsize)-1), min(memory.n-1, memory.batchsize), replace=false)
    return (
        memory.s[:, inds],
        memory.a[:, inds],
        memory.r[inds],
        memory.s[:, inds .+ 1],
        memory.d[inds]
    )
end


text:
"""
    ptrace(ρ, out)

Compute the partial trace of a linear `Operator` ρ by tracing out the subsystems specified by `out`. Multiple subsystems can be traced out by passing a sorted tuple of subsystem indices.

# Example
```jldoctest
Φ₊ = normalize!(basis(2,0)⊗basis(2,0) + basis(2,1)⊗basis(2,1)) # Bell pair
Ψ₊ = normalize!(basis(2,0)⊗basis(2,1) + basis(2,1)⊗basis(2,0)) # Bell pair
ρ  = 0.25 * Operator(Φ₊) + 0.75 * Operator(Ψ₊) # density matrix
ptrace(ρ,2) # trace out qubit 2
# output
2×2 Operator{Array{Float64,2},1} with dimensions 2
 0.5  0.0
 0.0  0.5
```
"""
function ptrace(σ::Operator, out)
    res = ptrace(Array(σ),out,dims(σ))
    return Operator(res[1],res[2])
end

"""
    ptrace(ψ, out)

Compute the partial trace of a state `Ket` or `Bra` ψ by tracing out the subsystems specified by `out`. Returns a density matrix. Multiple subsystems can be traced out by passing a sorted tuple of subsystem indices.

# Example
```jldoctest
julia> Φ₊ = normalize!(basis(2,0)⊗basis(2,0) + basis(2,1)⊗basis(2,1)) # Bell pair
4-d Ket{SparseVector{Float64,Int64},2} with dimensions 2⊗2
0.71∠0°|0,0⟩ + 0.71∠0°|1,1⟩

julia> ptrace(Φ₊,1) # trace out qubit 1
2×2 Operator{Array{Float64,2},1} with dimensions 2
 0.5  0.0
 0.0  0.5
```
"""
function ptrace(ψ::Ket, out)
    res = ptrace(Array(ψ),out,dims(ψ))
    return Operator(res[1], res[2])
end
ptrace(ψ::Bra, out) = ptrace(ψ',out)

# Method for when a single index is passed
ptrace(x, out::Int, sysdims::NTuple) = ptrace(x,(out,),sysdims)

# Dense partial trace
function ptrace(A::AbstractArray{T}, out::NTuple{no,Int}, sysdims::NTuple{ns,Int}) where {T,no,ns}
    # Make sure the input arguments make sense
    issorted(out) || throw(ArgumentError("subsystem indices $out must be sorted"))
    no>ns && throw(ArgumentError("more subsystem indices ($no) than number of subsystems ($ns)"))
    for i = 1:no
        out[i]>ns && throw(ArgumentError("subsystem index $(out[i]) larger than number of subsystems"))
    end
    # Calculate dimensions
    keep  = sorted_setdiff(ntuple(identity,Val(ns)),out)
    keepdims = gettuple(sysdims,keep)
    outdims  = gettuple(sysdims,out)
    # Generate iterators to loop over subscripts
    R, S = tensored_iterator(keepdims), tensored_iterator(outdims)
    # Initialize vectors for indexing purpose & output matrix
    subs_ii = Vector{Int}(undef,ns); subs_jj = Vector{Int}(undef,ns)
    B = zeros(T,prod(keepdims),prod(keepdims))
    # Main loop
    @inbounds for r in R, q in R
        # For each element in the reduced matrix
        i, j = tensored_sub2ind(keepdims,q), tensored_sub2ind(keepdims,r)
        for k=1:(ns-no)
            subs_ii[keep[k]] = q[k]; subs_jj[keep[k]] = r[k]
        end
        for s in S
            # Sum over all the corresponding "diagonal" elements in the full density matrix
            for o=1:no
                subs_ii[out[o]] = s[o]; subs_jj[out[o]] = s[o]
            end
            ii, jj = tensored_sub2ind(sysdims,subs_ii), tensored_sub2ind(sysdims,subs_jj)
            B[i,j] += _ptrace_ii_jj(A,ii,jj)
        end
    end
    return B, keepdims
end

# Dense partial trace kernel
@inline _ptrace_ii_jj(A::Matrix,ii,jj) = A[ii,jj]
@inline _ptrace_ii_jj(x::Vector,ii,jj) = x[ii]*x[jj]'

ext: jl


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dcc_slider

"""
    dcc_slider(;kwargs...)

A Slider component.
A slider component with a single handle.
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `marks` (optional): Marks on the slider.
The key determines the position (a number),
and the value determines what will show.
If you want to set the style of a specific mark point,
the value should be an object which
contains style and label properties.. marks has the following type: Dict with Strings as keys and values of type String | lists containing elements 'label', 'style'.
Those elements have the following types:
  - `label` (String; optional)
  - `style` (Dict; optional)
- `value` (Real; optional): The value of the input
- `drag_value` (Real; optional): The value of the input during a drag
- `className` (String; optional): Additional CSS class for the root DOM node
- `disabled` (Bool; optional): If true, the handles can't be moved.
- `dots` (Bool; optional): When the step value is greater than 1,
you can set the dots to true if you want to
render the slider with dots.
- `included` (Bool; optional): If the value is true, it means a continuous
value is included. Otherwise, it is an independent value.
- `min` (Real; optional): Minimum allowed value of the slider
- `max` (Real; optional): Maximum allowed value of the slider
- `tooltip` (optional): Configuration for tooltips describing the current slider value. tooltip has the following type: lists containing elements 'always_visible', 'placement'.
Those elements have the following types:
  - `always_visible` (Bool; optional): Determines whether tooltips should always be visible
(as opposed to the default, visible on hover)
  - `placement` (a value equal to: 'left', 'right', 'top', 'bottom', 'topLeft', 'topRight', 'bottomLeft', 'bottomRight'; optional): Determines the placement of tooltips
See https://github.com/react-component/tooltip#api
top/bottom{*} sets the _origin_ of the tooltip, so e.g. `topLeft`
will in reality appear to be on the top right of the handle
- `step` (Real; optional): Value by which increments or decrements are made
- `vertical` (Bool; optional): If true, the slider will be vertical
- `verticalHeight` (Real; optional): The height, in px, of the slider if it is vertical.
- `updatemode` (a value equal to: 'mouseup', 'drag'; optional): Determines when the component should update its `value`
property. If `mouseup` (the default) then the slider
will only trigger its value when the user has finished
dragging the slider. If `drag`, then the slider will
update its value continuously as it is being dragged.
If you want different actions during and after drag,
leave `updatemode` as `mouseup` and use `drag_value`
for the continuously updating value.
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `persistence` (Bool | String | Real; optional): Used to allow user interactions in this component to be persisted when
the component - or the page - is refreshed. If `persisted` is truthy and
hasn't changed from its previous value, a `value` that the user has
changed while using the app will keep that change, as long as
the new `value` also matches what was given originally.
Used in conjunction with `persistence_type`.
- `persisted_props` (Array of a value equal to: 'value's; optional): Properties whose user interactions will persist after refreshing the
component or the page. Since only `value` is allowed this prop can
normally be ignored.
- `persistence_type` (a value equal to: 'local', 'session', 'memory'; optional): Where persisted user changes will be stored:
memory: only kept in memory, reset on page refresh.
local: window.localStorage, data is kept after the browser quit.
session: window.sessionStorage, data is cleared once the browser quit.
"""
function dcc_slider(; kwargs...)
        available_props = Symbol[:id, :marks, :value, :drag_value, :className, :disabled, :dots, :included, :min, :max, :tooltip, :step, :vertical, :verticalHeight, :updatemode, :loading_state, :persistence, :persisted_props, :persistence_type]
        wild_props = Symbol[]
        return Component("dcc_slider", "Slider", "dash_core_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
using AssemblyCalculus
using AssemblyCalculus: overlap, freeze_synapses!, synapse_noise!
using Graphs
using Plots
using StatsBase
using LinearAlgebra: norm

K = 20
N = K^2
MEAN_DEGREE = K
P = MEAN_DEGREE / N
BETA = 0.0
TIMESTEPS = 10_000

g = erdos_renyi(N, P, is_directed=true)
brain = Brain(g, [N], [K], [BETA])
synapse_noise!(brain)

zstim = zero_stim(brain[1])

_, sp1 = simulate!([zstim], TIMESTEPS)
_, sp2 = simulate!([zstim], TIMESTEPS)

spikes1 = sp1[1]
spikes2 = sp2[1]

function visit_prob(neurons::Array{Int, N}, max_idx::Int) where N
    neuron_counts = counts(reshape(neurons, :), max_idx)
    neruon_prob = neuron_counts ./ sum(neuron_counts)
    return neruon_prob
end

prob1 = visit_prob(spikes1, N)
prob2 = visit_prob(spikes2, N)
adj = adjacency_matrix(g)
normalized_row_sums = sum(adj', dims=2) ./ sum(adj)

p1 = bar(prob1, alpha=0.5, label="Run 1")
bar!(p1, prob2, alpha=0.5, label="Run 2")
bar!(p1, normalized_row_sums, alpha=0.5, label="Normalized In-Degree")
xlabel!(p1, "Neuron ID")
ylabel!(p1, "Visitation Probability")
title!(p1, "Neuron Visitation Probability\n Same Network, Different Initial Condition, 10,000 iterations.")

p2 = histogram(prob1, alpha=0.5, label="Run 1")
histogram!(p2, prob2, alpha=0.5, label="Run 2")
histogram!(p2, normalized_row_sums, alpha=0.5, label="Normalized In-Degree")
xlabel!(p2, "Visitation Probaility")
ylabel!(p2, "Number of Neurons")

plot(p1, p2; size=(1200, 400), layout=@layout [a; b])


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

@test !MPI.Initialized()
MPI.Init()
@test MPI.Initialized()

@test 0 <= MPI.Comm_rank(MPI.COMM_WORLD) < MPI.Comm_size(MPI.COMM_WORLD)

@test MPI.has_cuda() isa Bool

if ArrayType != Array
    @test MPI.has_cuda()
end

@test !MPI.Finalized()
MPI.Finalize()
@test MPI.Finalized()


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

using WriteVTK

using CommonUtils
using Basis1D
using Basis2DQuad
using UniformQuadMesh
using UnPack

using SetupDG

N = 2
K1D = 20
t=2.0
CFL=0.9
flname = "euler2Dconv"
Np = (N+1)*(N+1)

"Mesh related variables"
VX, VY, EToV = uniform_quad_mesh(Int(round(3/2*K1D)),K1D)
@. VX = VX*7.5+2.5
@. VY = VY*5

rd = init_reference_quad(N,gauss_lobatto_quad(0,0,N))
md = init_mesh((VX,VY),EToV,rd)
@unpack x,y = md;

@unpack VDM = rd
rp,sp = equi_nodes_2D(10)
Vp = vandermonde_2D(N,rp,sp)/VDM

x_vtk = []
y_vtk = []
for k = 1:Int(round(3/2*K1D))
    push!(x_vtk,x[1,k])
    push!(x_vtk,x[2,k])
    push!(x_vtk,x[3,k])
end
for k = 1:K1D
    push!(y_vtk,y[1,1+Int(round((k-1)*3/2*K1D))])
    push!(y_vtk,y[4,1+Int(round((k-1)*3/2*K1D))])
    push!(y_vtk,y[7,1+Int(round((k-1)*3/2*K1D))])
end
x_vtk = x_vtk .* ones((N+1)*K1D,1)'
y_vtk = y_vtk' .* ones((N+1)*Int(round(3/2*K1D)),1)
rho_vtk,rhou_vtk,rhov_vtk,E_vtk,sigmax_vtk,sigmay_vtk = (zeros(size(x_vtk)) for i = 1:6)

# vtk_grid("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH",x,y) do vtk
#     vtk["rho"] = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,xt=0.014317166411179642=$XLENGTH,rho,dmr.txt")
#     vtk["rhou"] = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rhou,dmr.txt")
#     vtk["rhov"] = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,rhov,dmr.txt")
#     vtk["E"] = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,x=$XLENGTH,E,dmr.txt")
# end

vtk_grid("N=$N,K1D=$K1D,t=$t,CFL=$CFL,$flname",x_vtk,y_vtk) do vtk
    rho = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,rho,$flname.txt")
    rhou = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,rhou,$flname.txt")
    rhov = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,rhov,$flname.txt")
    E = readdlm("N=$N,K1D=$K1D,t=$t,CFL=$CFL,E,$flname.txt")

    for k = 1:Int(round(3/2*K1D))*K1D
        kx = mod1(k,Int(round(3/2*K1D)))
        ky = div(k-1,Int(round(3/2*K1D)))+1

        rho_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = rho[:,k]
        rhou_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = rhou[:,k]
        rhov_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = rhov[:,k]
        E_vtk[(kx-1)*(N+1)+1:(kx-1)*(N+1)+N+1, (ky-1)*(N+1)+1:(ky-1)*(N+1)+N+1] = E[:,k]
    end

    vtk["rho"]  = rho_vtk
    vtk["rhou"] = rhou_vtk
    vtk["rhov"] = rhov_vtk
    vtk["E"]    = E_vtk
end

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
module HueLights

using HTTP, JSON, Colors

"""
    get_bridge_ip()

Get the IP address of the Hue bridge. First check the environment variables,
and if key `HUE_BRIDGE_IP` is not present, use the Philips Hue website. Return
the address as a string or throw an error if none is found.
"""
function get_bridge_ip()
	if haskey(ENV, "HUE_BRIDGE_IP")
		return ENV["HUE_BRIDGE_IP"]
	end
	r = HTTP.get("https://www.meethue.com/api/nupnp")
	if r.status != 200
		error("Did not get successful response from https://www.meethue.com/api/nupnp")
	end
	# body should be a 1-element vector,
	# something like [{"id":"012345abcde","internalipaddress":"192.168.0.102"}]
	js = r.body |> String |> JSON.parse
	if length(js) < 1
		error("Did not find any Hue bridges.")
	end
	return js[1]["internalipaddress"]
end

function hue_api_addr()
	return "http://" * ENV["HUE_BRIDGE_IP"] * "/api/" * ENV["HUE_USERNAME"]
end

const HUE = hue_api_addr()

function get_hue_lights()
	r = HTTP.get(HUE * "/lights")
	return get_hue_lights(String(r))
end
# get light names given the string describing all lights
function get_hue_lights(str::String)
	dict = JSON.parse(str)
	lightind = 1
	lightnames = String[]
	while haskey(dict, string(lightind))
		push!(lightnames, dict[string(lightind)]["name"])
		lightind += 1
	end
	return lightnames
end

function set_light_state!(light::Int, body::Dict)
	HTTP.put(HUE * "/lights/$(light)/state"; body = JSON.json(body))
end

# brightness false by default because it seems like not all color specs
# can express the same brightness range, e.g. RGB(0,0,1) gives a brightness Y of
# ~.07 on a [0,1] scale.
function set_light_color!(light, color; brightness = false)
	cie = convert(xyY, color)
	if cie.Y == 0 # no brightness, light should be off
		set_light_state!(light, Dict("on" => false))
	end
	state = Dict{String, Any}("xy" => [cie.x, cie.y], "on" => true)
	if brightness
		state["bri"] = round(Int, 254 * cie.Y)
	end
	set_light_state!(light, state)
end

function get_light_state(light::Int)
	d = JSON.parse(String(HTTP.get(HUE * "/lights/$(light)")))
	return d["state"]
end

macro huewarn(light, expr)
	quote
		local val = $(esc(expr))
		@schedule begin
			l = $(light)
			state = get_light_state(l)
			set_light_state!(l, Dict("on" => true, "bri" => 254))
			set_light_color!(l, RGB(0,1,0))
			sleep(1)
			set_light_color!(l, RGB(1,0,0))
			sleep(1)
			set_light_color!(l, RGB(0,1,0))
			sleep(1)
			set_light_state!(l, state)
		end
		return val
	end
end

end # module


text:
abstract type AbstractIndexing end
struct ZeroBasedIndexing <: AbstractIndexing end
struct OneBasedIndexing <: AbstractIndexing end

first_index(::ZeroBasedIndexing) = 0
first_index(::OneBasedIndexing) = 1
shift(x, ::ZeroBasedIndexing, ::ZeroBasedIndexing) = x
shift(x::Integer, ::ZeroBasedIndexing, ::OneBasedIndexing) = x + 1
shift(x::Array{<:Integer}, ::ZeroBasedIndexing, ::OneBasedIndexing) = x .+ 1
shift(x::Integer, ::OneBasedIndexing, ::ZeroBasedIndexing) = x - 1
shift(x, ::OneBasedIndexing, ::OneBasedIndexing) = x

mutable struct SparseMatrixCSRtoCSC{Tv,Ti<:Integer,I<:AbstractIndexing}
    indexing::I
    m::Int # Number of rows
    n::Int # Number of columns
    colptr::Vector{Ti}
    rowval::Vector{Ti}
    nzval::Vector{Tv}
    function SparseMatrixCSRtoCSC{Tv, Ti, I}(n) where {Tv, Ti<:Integer, I}
        A = new{Tv, Ti, I}()
        A.n = n
        A.colptr = zeros(Ti, n + 1)
        return A
    end
end
function allocate_nonzeros(A::SparseMatrixCSRtoCSC{Tv, Ti}) where {Tv, Ti}
    for i in 3:length(A.colptr)
        A.colptr[i] += A.colptr[i - 1]
    end
    A.rowval = Vector{Ti}(undef, A.colptr[end])
    A.nzval = Vector{Tv}(undef, A.colptr[end])
end
function final_touch(A::SparseMatrixCSRtoCSC)
    for i in length(A.colptr):-1:2
        A.colptr[i] = shift(A.colptr[i - 1], ZeroBasedIndexing(), A.indexing)
    end
    A.colptr[1] = first_index(A.indexing)
end
function _allocate_terms(colptr, indexmap, terms)
    for term in terms
        colptr[indexmap[term.scalar_term.variable].value + 1] += 1
    end
end
function allocate_terms(A::SparseMatrixCSRtoCSC, indexmap, func)
    _allocate_terms(A.colptr, indexmap, func.terms)
end
function _load_terms(colptr, rowval, nzval, indexmap, terms, offset)
    for term in terms
        ptr = colptr[indexmap[term.scalar_term.variable].value] += 1
        rowval[ptr] = offset + term.output_index
        nzval[ptr] = -term.scalar_term.coefficient
    end
end
function load_terms(A::SparseMatrixCSRtoCSC, indexmap, func, offset)
    _load_terms(A.colptr, A.rowval, A.nzval, indexmap, func.terms, shift(offset, OneBasedIndexing(), A.indexing))
end

"""
    Base.convert(::Type{SparseMatrixCSC{Tv, Ti}}, A::SparseMatrixCSRtoCSC{Tv, Ti, I}) where {Tv, Ti, I}

Converts `A` to a `SparseMatrixCSC`. Note that the field `A.nzval` is **not
copied** so if `A` is modified after the call of this function, it can still
affect the value returned. Moreover, if `I` is `OneBasedIndexing`, `colptr`
and `rowval` are not copied either, i.e., the conversion is allocation-free.
"""
function Base.convert(::Type{SparseMatrixCSC{Tv, Ti}}, A::SparseMatrixCSRtoCSC{Tv, Ti}) where {Tv, Ti}
    return SparseMatrixCSC{Tv, Ti}(
        A.m,
        A.n,
        shift(A.colptr, A.indexing, OneBasedIndexing()),
        shift(A.rowval, A.indexing, OneBasedIndexing()),
        A.nzval,
    )
end


text:
"""
    MLP(h::Int...;kwargs...)

Creates a multi layer perceptron according to given hidden states.
First hidden state is equal to input size and the last one equal to output size.

    (m::MLP)(x)

Runs MLP with given input `x`. `prob` is the dropout probability.

# Keywords

* `winit=xavier`: weight initialization distribution
* `bias=zeros`: bias initialization distribution
* `activation=ReLU()`: activation layer or function
* `atype=KnetLayers.arrtype` : array type for parameters.
   Default value is KnetArray{Float32} if you have gpu device. Otherwise it is Array{Float32}
*  pdrop=0: dropout probability between layers

"""
mutable struct MLP <: Layer
    layers::Tuple{Vararg{Linear}}
    activation::Activation
    drop::Dropout
end

function MLP(h::Int...; winit=xavier, binit=zeros, activation=ReLU(), atype=arrtype, pdrop=0)
    singlelayer(input,ouput) = Linear(input=input, output=ouput, winit=winit, binit=binit, atype=atype)
    MLP(singlelayer.(h[1:end-1], h[2:end]),activation,Dropout(pdrop))
end

function (m::MLP)(x)
    for layer in m.layers
        x = layer(m.drop(x))
        if layer !== m.layers[end]
            x = m.activation(x)
        end
    end
    return x
end


text:
module Datasets

include("io.jl")
include("MNIST.jl")
include("CIFAR10.jl")
include("CIFAR100.jl")
include("PTBLM.jl")

end


text:
# The returned result from a langevin kernel discrepancy

abstract type LangevinKernelBaseResult end

type LangevinKernelResult <: LangevinKernelBaseResult
    # n x p matrix of distinct sample points
    points::Array{Float64}
    # n x 1 vector of weights associated with each sample
    weights::Array{Float64}
    # the sum of the kernels (squared) that bound the IPM
    discrepancy2::Float64
    # Time taken to set up and solve each optimization program
    solvetime::Float64
end

type LangevinKernelCheckpointsResult <: LangevinKernelBaseResult
    # n x p matrix of distinct sample points
    points::Array{Float64}
    # n x 1 vector of weights associated with each sample
    weights::Array{Float64}
    # the sum of the kernels (squared) that bound the IPM by checkpoints
    discrepancy2::Array{Float64}
    #checkpoints passed into the function
    checkpoints::Array{Int}
    # Time taken to set up and solve each optimization program
    solvetime::Float64
end


text:
#md# Big words will be placed closer to the center
using WordCloud
wc = wordcloud(
    processtext(open(pkgdir(WordCloud) * "/res/alice.txt"), stopwords=WordCloud.stopwords_en ∪ ["said"]),
    angles=0, density=0.55,
    maskshape=squircle, rt=2.5 * rand(),
    run=initwords!)
placewords!(wc, style=:gathering, level=5, centerlargestword=true)
pin(wc, "Alice") do # keep "Alice" in the center
    generate!(wc, teleporting=0.7) # don't teleport largest 30% words
end
println("results are saved to gathering.svg")
paint(wc, "gathering.svg")
wc
#eval# runexample(:gathering)
#md# ![](gathering.svg)

max_stars_repo_name: guo-yong-zhi/WordCloud

max_issues_repo_name: guo-yong-zhi/WordCloud

max_forks_repo_name: guo-yong-zhi/WordCloud


text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# # 0105 Custom Types

# ## Some stuff to set-up the environment..

cd(@__DIR__)
using Pkg
Pkg.activate(".")
## If using a Julia version different than 1.7 please uncomment and run the following line (reproductibility guarantee will however be lost)
## Pkg.resolve()
## Pkg.instantiate() # run this if you didn't in Segment 01.01
using Random
Random.seed!(123)
using InteractiveUtils # loaded automatically when working... interactively

# ## "Type" of types

# In Julia primitive, composite and abstract types can all be defined by the user

primitive type APrimitiveType 819200 end # name and size in bit - multiple of 8 and below 8388608 (1MB)
primitive type APrimitiveType2 819200 end
819200/(8*1024)
struct ACompositeType end                # fields, constructors.. we'll see this later in details
abstract type AnAbstractType end         # no objects, no instantialisation of objects

# ## Composite types
mutable struct Foo # starting with a capital letter
    field1
    field2::String
    field3::ACompositeType
end

# mutable struct Foo # error! can't change a struct after I defined it
# end

fieldnames(Foo)

o = Foo(123,"aaa",ACompositeType()) # call the default constructor (available automatically) - order matters!
typeof(o)

# Outer constructor
function Foo(f2,f3=ACompositeType()) # "normal" functions just happens it has the name of the object to create
    if startswith(f2,"Booo")         # put whatever logic you wish
        return nothing
    end
    return Foo(123,f2,f3)            # call the default constructor
end

o = Foo(123,"aaa", ACompositeType()) # call the default constructor
o = Foo("blaaaaa")                   # call the outer constructor we defined

o.field1       # access fields
o.field1 = 321 # modify field (because type defined as "mutable" !!!)
o

function Base.show(io::IO, x::Foo)             # function(o) rather than o.function()
    println(io,"My custom representation for Foo objects")
    println(io,"Field1: $(o.field1)")
    println(io,"Field2: $(o.field2)")
end
o

# Inner constructor
mutable struct Foo2
    field1::Int64
    field2::String
    function Foo2(f1,f2,f3)
        ## ... logic
        return new(f1+f2,f3)
    end
end
Foo2(1,2,"aaa")

# !!! tip
#     If any inner constructor method is defined, no default constructor method is provided.

## Foo2(1,"aaa") # Error, no default constructor !

# ## Parametric types

struct Point{T<:Number} # T must be a child of type "Number"
   x::T
   y::T
end

o = Point(1,2)
Point(1.0,2.)
## Point(1,2.0) # error !

function Point(x::T, y::T=zero(T)) where {T}
    return Point(x,y)
end
Point(2)
Point(1.5)

abstract type Figure{T<:Number} end

a = Array{Int64,2}(undef,2,2) # Array is nothing else than a parametric type with 2 parameters
typeof(a)
eltype(a)

# As we see for arrays, parameters doesn't need to be _types_, but can be any value of a bits type (in practice an integer value) :

struct MyType{T,N}
  data::Array{T,N}
end

intMatrixInside = MyType([1 2 3; 4 5 6])
floatVectorInside = MyType([1 2 3])

function getPlane(o::MyType{T,N},dim,pos) where {T,N}
sizes = size(o.data)
if length(sizes) > N
  error("Dim over the dimensions of the data")
elseif sizes[dim] < pos
  error("Non enought elements in dimension $dim to cut at $pos")
end
return selectdim(o.data,dim,pos)
end

getPlane(intMatrixInside,1,2)

# A package where non-type parameters are emploied to boost speed is [StaticArray.jl](https://github.com/JuliaArrays/StaticArrays.jl) where one parameter is the _size_ of the array that hence become known at compile time

# ## Inheritance

abstract type MyOwnGenericAbstractType end                       # the highest-level
abstract type MyOwnAbstractType1 <: MyOwnGenericAbstractType end # child of MyOwnGenericAbstractType
abstract type MyOwnAbstractType2 <: MyOwnGenericAbstractType end # also child of MyOwnGenericAbstractType
mutable struct AConcreteTypeA <: MyOwnAbstractType1
  f1::Int64
  f2::Int64
end
mutable struct AConcreteTypeB <: MyOwnAbstractType1
  f1::Float64
end
mutable struct AConcreteTypeZ <: MyOwnAbstractType2
  f1::String
end
oA = AConcreteTypeA(2,10)
oB = AConcreteTypeB(1.5)
oZ = AConcreteTypeZ("aa")

supertype(AConcreteTypeA)
subtypes(MyOwnAbstractType1)

# !!! tip
#     When multiple methods are available for an object, function calls are dispatched to the most stricter method, i.e. the one defined over the exact parameter's type or their immediate parents

function foo(a :: MyOwnGenericAbstractType)                      # good for everyone
  println("Default implementation: $(a.f1)")
end
foo(oA) # Default implementation: 2
foo(oB) # Default implementation: 1.5
foo(oZ) # Default implementation: aa
function foo(a :: MyOwnAbstractType1)                            # specialisation for MyOwnAbstractType1
  println("A more specialised implementation: $(a.f1*4)")
end
foo(oA) # A more specialised implementation: 8
foo(oB) # A more specialised implementation: 6.0
foo(oZ) # Default implementation: aa                             # doesn't match the specialisation, default to foo(a :: MyOwnGenericAbstractType)
function foo(a :: AConcreteTypeA)
     println("A even more specialised implementation: $(a.f1 + a.f2)")
end
foo(oA) # A even more specialised implementation: 12
foo(oB) # A more specialised implementation: 6.0
foo(oZ) # Default implementation: aa

# !!! warning
#     Attention to the inheritance for parametric types. If it is true that `Vector{Int64} <: AbstractVector{Int64}` and `Int64 <: Number`, it is FALSE that `AbstractVector{Int64} <: AbstractVector{Number}`. If you want to allow a function parameter to be a vector of numbers, use instead templates explicitly, e.g. `foo(x::AbstractVector{T}) where {T<:Number} = return sum(x)`

Vector{Int64} <: AbstractVector{Int64}
Int64 <: Number
Vector{Int64} <: Vector{Number}
AbstractVector{Int64} <: AbstractVector{Number}

# ## Object-oriented model

# OO model based on _composition_

struct Shoes
   shoesType::String
   colour::String
end
struct Person
  myname::String
  age::Int64
end
struct Student
   p::Person        # by referencing a `Person`` object, we do not need to repeat its fields
   school::String
   shoes::Shoes     # same for `shoes`
end
struct Employee
   p::Person
   monthlyIncomes::Float64
   company::String
   shoes::Shoes
end
gymShoes = Shoes("gym","white")
proShoes = Shoes("classical","brown")
Marc     = Student(Person("Marc",15),"Divine School",gymShoes)
MrBrown  = Employee(Person("Brown",45),3200.0,"ABC Corporation Inc.", proShoes)

function printMyActivity(self::Student)
   println("Hi! I am $(self.p.myname), I study at $(self.school) school, and I wear $(self.shoes.colour) shoes") # I can use the dot operator chained...
end
function printMyActivity(self::Employee)
  println("Good day. My name is $(self.p.myname), I work at $(self.company) company and I wear $(self.shoes.colour) shoes")
end

printMyActivity(Marc)     # Hi! I am Marc, ...
printMyActivity(MrBrown)  # Good day. My name is MrBrown, ...

# OO models based on Specialisation (Person → Student) or Weack Relation (Person → Shoes) instead of Composition (Person → Arm) can be implemented using third party packages, like e.g. [SimpleTraits.jl](https://github.com/mauro3/SimpleTraits.jl) or [OOPMacro.jl](https://github.com/ipod825/OOPMacro.jl)

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
__precompile__(true)
module Setfield
using MacroTools
using MacroTools: isstructdef, splitstructdef, postwalk
using Requires: @require

if VERSION < v"1.1-"
    using Future: copy!
end

include("setindex.jl")
include("lens.jl")
include("sugar.jl")
include("functionlenses.jl")

# To correctly dispatch to `show(::IO, ::CustomLens)` when it is defined by a
# user, we avoid defining the generic `show(::IO, ::Lens)`.  This way, we can
# safely call `show` inside `ComposedLens` without worrying about the
# `StackOverflowError` that can be easily triggered in the previous approach.
# See also:
# * https://github.com/jw3126/Setfield.jl/pull/86
# * https://github.com/jw3126/Setfield.jl/pull/88
for n in names(Setfield, all=true)
    T = getproperty(Setfield, n)
    if T isa Type && T <: Lens && (T === ComposedLens || has_atlens_support(T))
        @eval Base.show(io::IO, l::$T) = _show(io, nothing, l)
    end
end

function __init__()
    @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" begin
        setindex(a::StaticArrays.StaticArray, args...) =
            Base.setindex(a, args...)
    end
end

end


text:
@generated function compute_integral(Is)
    SU = Ultraspherical(1, 0..1)^2
    Q  = DefiniteIntegral(SU)

    quote
        return Number( $Q * Fun($SU, ApproxFun.transform($SU, Is)) )
    end
end

@generated function compute_derivative1(q)
    Dx = Derivative(Chebyshev(0..1)^2, [1,0])
    CU = Conversion(Ultraspherical(1, 0..1) ⊗ Chebyshev(0..1), Ultraspherical(1, 0..1)^2)

    quote
        return values($CU * ($Dx * q))
    end
end

@generated function compute_derivative2(q)
    Dy = Derivative(Chebyshev(0..1)^2, [0,1])
    CU = Conversion(Chebyshev(0..1) ⊗ Ultraspherical(1, 0..1), Ultraspherical(1, 0..1)^2)

    quote
        return values($CU * ($Dy * q))
    end
end

@generated function compute_canonical_invariant(pinv::Union{PoincareInvariant2nd{DT,ND,NC,NV},PoincareInvariant2ndCanonical{DT,ND,NC,NV}}, q::AbstractArray{DT,2}, p::AbstractArray{DT,2}) where {DT,ND,NC,NV}
    SC = Chebyshev(0..1)^2
    SU = Ultraspherical(1, 0..1)^2

    Js = zeros(DT, NV)
    qσ = zeros(DT, ND, NV)
    qτ = zeros(DT, ND, NV)
    pσ = zeros(DT, ND, NV)
    pτ = zeros(DT, ND, NV)

    quote
        for j in 1:size(p,1)
            fq = Fun($SC, ApproxFun.transform($SC, q[j,:]))
            fp = Fun($SC, ApproxFun.transform($SC, p[j,:]))

            # compute derivatives of q and store values
            $qσ[j,:] .= compute_derivative1(fq)
            $qτ[j,:] .= compute_derivative2(fq)

            # compute derivatives of p and store values
            $pσ[j,:] .= compute_derivative1(fp)
            $pτ[j,:] .= compute_derivative2(fp)
        end

        # compute integrands of integral invariants at all points
        for k in 1:NV
            $Js[k] = $pσ[:,k]' * $qτ[:,k] - $qσ[:,k]' * $pτ[:,k]
        end

        # compute canonical integral invariant
        return compute_integral($Js)
    end
end

@generated function compute_noncanonical_invariant(pinv::PoincareInvariant2nd{DT,ND,NC,NV}, t::DT, q::AbstractArray{DT,2}) where {DT,ND,NC,NV}
    SC = Chebyshev(0..1)^2
    SU = Ultraspherical(1, 0..1)^2
    CU = Conversion(SC, SU)

    Is = zeros(DT, NV)
    qs = zeros(DT, ND, NV)
    qσ = zeros(DT, ND, NV)
    qτ = zeros(DT, ND, NV)
    Ω  = zeros(DT, ND, ND)

    quote
        for j in 1:size(q,1)
            fq = Fun($SC, ApproxFun.transform($SC, q[j,:]))

            # obtain values of q at the same points as the derivatives
            $qs[j,:] .= itransform($SU, resize!(($CU * fq).coefficients, NC))

            # compute derivatives of q and store values
            $qσ[j,:] .= compute_derivative1(fq)
            $qτ[j,:] .= compute_derivative2(fq)
        end

        # compute integrands of integral invariants at all points
        for k in 1:NV
            pinv.ω(t, $qs[:,k], $Ω)
            $Is[k] = vector_matrix_vector_product($qτ[:,k], $Ω, $qσ[:,k])
        end

        # compute noncanonical integral invariant
        return compute_integral($Is)
    end
end

@generated function compute_noncanonical_correction(pinv::PoincareInvariant2nd{DT,ND,NC,NV}, t::DT, q::AbstractArray{DT,2}, λ::AbstractArray{DT,2}) where {DT,ND,NC,NV}
    SC = Chebyshev(0..1)^2
    SU = Ultraspherical(1, 0..1)^2
    CU = Conversion(SC, SU)

    Ks = zeros(DT, NV)
    qs = zeros(DT, ND, NV)
    λs = zeros(DT, ND, NV)
    qσ = zeros(DT, ND, NV)
    qτ = zeros(DT, ND, NV)
    λσ = zeros(DT, ND, NV)
    λτ = zeros(DT, ND, NV)
    Ω  = zeros(DT, ND, ND)
    D²ϑ= zeros(DT, ND, ND)

    quote
        for j in 1:size(q,1)
            fq = Fun($SC, ApproxFun.transform($SC, q[j,:]))
            fλ = Fun($SC, ApproxFun.transform($SC, λ[j,:]))

            # obtain values of q at the same points as the derivatives
            $qs[j,:] .= itransform($SU, resize!(($CU * fq).coefficients, NC))
            $λs[j,:] .= itransform($SU, resize!(($CU * fλ).coefficients, NC))

            # compute derivatives of q and store values
            $qσ[j,:] .= compute_derivative1(fq)
            $qτ[j,:] .= compute_derivative2(fq)

            # compute derivatives of λ and store values
            $λσ[j,:] .= compute_derivative1(fλ)
            $λτ[j,:] .= compute_derivative2(fλ)
        end

        # compute integrands of integral invariants at all points
        for k in 1:NV
            pinv.ω(t, $qs[:,k], $Ω)
            $Ks[k] = vector_matrix_vector_product($λτ[:,k], $Ω, $λσ[:,k])

            for j in 1:length(pinv.D²ϑ)
                pinv.D²ϑ[j](t, $qs[:,k], $D²ϑ)
                $Ks[k] += 2 * $λs[j,k] * vector_matrix_vector_product($qτ[:,k], $D²ϑ, $λσ[:,k])
                $Ks[k] -= 2 * $λs[j,k] * vector_matrix_vector_product($qσ[:,k], $D²ϑ, $λτ[:,k])
            end
        end

        # compute correction to noncanonical integral invariant
        return compute_integral($Ks)
    end
end


text:
using Base.Test, Mehs

@testset "Mehs" begin
    # test promote rules
    @test promote_type(Meh, Meh) == Meh
    @test promote_type(Meh, Int) == Union{Meh, Int}
    @test promote_type(Int, Meh) == Union{Meh, Int}
    @test promote_type(Int, Any) == Any
    @test promote_type(Any, Any) == Any
    @test promote_type(Meh, Any) == Any
    @test promote_type(Any, Meh) == Any
    @test promote_type(Union{Int, Meh}, Meh) == Union{Int, Meh}
    @test promote_type(Meh, Union{Int, Meh}) == Union{Int, Meh}
    @test promote_type(Union{Int, Meh}, Int) == Union{Int, Meh}
    @test promote_type(Int, Union{Int, Meh}) == Union{Int, Meh}
    @test promote_type(Any, Union{Int, Meh}) == Any
    @test promote_type(Union{Int, Meh}, Union{Int, Meh}) == Union{Int, Meh}
    @test promote_type(Union{Float64, Meh}, Union{String, Meh}) == Any
    @test promote_type(Union{Float64, Meh}, Union{Int, Meh}) == Union{Float64, Meh}
    @test promote_type(Union{Void, Meh, Int}, Float64) == Any

    bit_operators = [&, |, ⊻]

    arithmetic_operators = [+, -, *, /, ^, Base.div, Base.mod, Base.fld, Base.rem]

    elementary_functions = [abs, abs2, sign,
                            acos, acosh, asin, asinh, atan, atanh, sin, sinh,
                            conj, cos, cosh, tan, tanh,
                            exp, exp2, expm1, log, log10, log1p, log2,
                            exponent, sqrt, gamma, lgamma,
                            identity, zero,
                            iseven, isodd, ispow2,
                            isfinite, isinf, isnan, iszero,
                            isinteger, isreal, isimag,
                            isempty, transpose, ctranspose, float]

    rounding_functions = [ceil, floor, round, trunc]

    # All unary operators return meh when evaluating meh
    for f in [!, +, -]
        @test ismeh(f(meh))
    end

    # All elementary functions return meh when evaluating meh
    for f in elementary_functions
        @test ismeh(f(meh))
    end

    # All rounding functions return meh when evaluating meh as first argument
    for f in rounding_functions
        @test ismeh(f(meh))
        @test ismeh(f(meh, 1))
        @test ismeh(f(meh, 1, 1))
        @test ismeh(f(Union{Int, Meh}, meh))
        @test_throws MehException f(Int, meh)
    end

    @test zero(Union{Int, Meh}) === 0
    @test zero(Union{Float64, Meh}) === 0.0

    for T in (subtypes(Dates.DatePeriod)...,
              subtypes(Dates.TimePeriod)...)
        @test zero(Union{T, Meh}) === T(0)
    end

    # Comparison operators
    @test (meh == meh) === meh
    @test (1 == meh) === meh
    @test (meh == 1) === meh
    @test (meh != meh) === meh
    @test (1 != meh) === meh
    @test (meh != 1) === meh
    @test isequal(meh, meh)
    @test !isequal(1, meh)
    @test !isequal(meh, 1)
    @test (meh < meh) === meh
    @test (meh < 1) === meh
    @test (1 < meh) === meh
    @test (meh <= meh) === meh
    @test (meh <= 1) === meh
    @test (1 <= meh) === meh
    @test !isless(meh, meh)
    @test !isless(meh, 1)
    @test isless(1, meh)

    # All arithmetic operators return meh when operating on two meh's
    # All arithmetic operators return meh when operating on a scalar and an meh
    # All arithmetic operators return meh when operating on an meh and a scalar
    for f in arithmetic_operators
        @test ismeh(f(meh, meh))
        @test ismeh(f(1, meh))
        @test ismeh(f(meh, 1))
    end

    # All bit operators return meh when operating on two meh's
    for f in bit_operators
        @test ismeh(f(meh, meh))
    end

    @test ismeh(meh & true)
    @test ismeh(true & meh)
    @test !(meh & false)
    @test !(false & meh)
    @test ismeh(meh | false)
    @test ismeh(false | meh)
    @test meh | true
    @test true | meh
    @test ismeh(xor(meh, true))
    @test ismeh(xor(true, meh))
    @test ismeh(xor(meh, false))
    @test ismeh(xor(false, meh))

    @test ismeh(meh & 1)
    @test ismeh(1 & meh)
    @test ismeh(meh | 1)
    @test ismeh(1 | meh)
    @test ismeh(xor(meh, 1))
    @test ismeh(xor(1, meh))

    @test ismeh("a" * meh)
    @test ismeh(meh * "a")

    @test sprint(show, meh) == "meh"
    @test sprint(showcompact, meh) == "meh"
    @test sprint(show, [meh]) == "$Meh[meh]"
    @test sprint(show, [1 meh]) == "$(Union{Int, Meh})[1 meh]"
    b = IOBuffer()
    display(TextDisplay(b), [meh])
    @test String(take!(b)) == "1-element Array{$Meh,1}:\n meh"
    b = IOBuffer()
    display(TextDisplay(b), [1 meh])
    @test String(take!(b)) == "1×2 Array{$(Union{Int, Meh}),2}:\n 1  meh"

    x = Mehs.replace([1, 2, meh, 4], 3)
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (4,)
    @test collect(x) == collect(1:4)
    @test collect(x) isa Vector{Int}
    x = Mehs.replace([1, 2, meh, 4], 3.0)
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (4,)
    @test collect(x) == collect(1:4)
    @test collect(x) isa Vector{Int}
    x = Mehs.replace([1 2; meh 4], 3)
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (2, 2)
    @test collect(x) == [1 2; 3 4]
    @test collect(x) isa Matrix{Int}
    x = Mehs.replace((v for v in [meh, 1, meh, 2, 4]), 0)
    @test length(x) == 5
    @test size(x) == (5,)
    @test eltype(x) === Any
    @test collect(x) == [0, 1, 0, 2, 4]
    @test collect(x) isa Vector{Int}

    x = Mehs.fail([1, 2, 3, 4])
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (4,)
    @test collect(x) == [1, 2, 3, 4]
    @test collect(x) isa Vector{Int}
    x = Mehs.fail([1 2; 3 4])
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (2, 2)
    @test collect(x) == [1 2; 3 4]
    @test collect(x) isa Matrix{Int}
    @test_throws MehException collect(Mehs.fail([1, 2, meh, 4]))
    x = Mehs.fail(v for v in [1, 2, 4])
    @test eltype(x) === Any
    @test length(x) == 3
    @test size(x) == (3,)
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}

    x = Mehs.skip([1, 2, meh, 4])
    @test eltype(x) === Int
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}
    x = Mehs.skip([1  2; meh 4])
    @test eltype(x) === Int
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}
    x = collect(Mehs.skip([meh]))
    @test eltype(x) === Union{}
    @test isempty(collect(x))
    @test collect(x) isa Vector{Union{}}
    x = collect(Mehs.skip(Union{Int, Meh}[]))
    @test eltype(x) === Int
    @test isempty(collect(x))
    @test collect(x) isa Vector{Int}
    x = Mehs.skip([meh, meh, 1, 2, meh, 4, meh, meh])
    @test eltype(x) === Int
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}
    x = Mehs.skip(v for v in [meh, 1, meh, 2, 4])
    @test eltype(x) === Any
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}

    @test Mehs.coalesce(meh, 1) === 1
    @test Mehs.coalesce(1, meh) === 1
    @test Mehs.coalesce(meh, meh) === meh
    @test Mehs.coalesce.([meh, 1, meh], 0) == [0, 1, 0]
    @test Mehs.coalesce.([meh, 1, meh], 0) isa Vector{Int}
    @test Mehs.coalesce.([meh, 1, meh], [0, 10, 5]) == [0, 1, 5]
    @test Mehs.coalesce.([meh, 1, meh], [0, 10, 5]) isa Vector{Int}
    @test isequal(Mehs.coalesce.([meh, 1, meh], [0, meh, meh]), [0, 1, meh])
    # Fails in Julia 0.6 and 0.7.0-DEV.1556
    @test_broken Mehs.coalesce.([meh, 1, meh], [0, meh, meh]) isa Vector{Union{Meh, Int}}

    @test levels(1:1) == levels([1]) == levels([1, meh]) == levels([meh, 1]) == [1]
    @test levels(2:-1:1) == levels([2, 1]) == levels([2, meh, 1]) == [1, 2]
    @test levels([meh, "a", "c", meh, "b"]) == ["a", "b", "c"]
    @test levels([Complex(0, 1), Complex(1, 0), meh]) == [Complex(0, 1), Complex(1, 0)]
    @test levels(sparse([0 3 2])) == [0, 2, 3]
    @test typeof(levels([1])) === typeof(levels([1, meh])) === Vector{Int}
    @test typeof(levels(["a"])) === typeof(levels(["a", meh])) === Vector{String}
    @test typeof(levels(sparse([1]))) === Vector{Int}
    @test isempty(levels([meh]))
    @test isempty(levels([]))

    x = convert(Vector{Union{Int, Meh}}, [1.0, meh])
    @test isa(x, Vector{Union{Int, Meh}})
    @test isequal(x, [1, meh])
    x = convert(Vector{Union{Int, Meh}}, [1.0])
    @test isa(x, Vector{Union{Int, Meh}})
    @test x == [1]
    x = convert(Vector{Union{Int, Meh}}, [meh])
    @test isa(x, Vector{Union{Int, Meh}})
    @test isequal(x, [meh])

    @test Mehs.T(Union{Int, Meh}) == Int
    @test Mehs.T(Any) == Any
    @test Mehs.T(Meh) == Union{}

    @test isequal(mehs(1), [meh])
    @test isequal(mehs(Int, 1), [meh])
    @test mehs(Int, 1) isa Vector{Union{Int, Meh}}
    @test isequal(mehs(Union{Int, Meh}, 1, 2), [meh meh])
    @test mehs(Union{Int, Meh}, 1, 2) isa Matrix{Union{Int, Meh}}
    @test Union{Int, Meh}[1,2,3] == (Union{Int, Meh})[1,2,3]

    @test convert(Union{Int, Meh}, 1.0) == 1

    # AbstractArray{>:Meh}

    @test ismeh([1, meh] == [1, meh])
    @test ismeh(["a", meh] == ["a", meh])
    @test ismeh(Any[1, meh] == Any[1, meh])
    @test ismeh(Any[meh] == Any[meh])
    @test ismeh([meh] == [meh])
    @test ismeh(Any[meh, 2] == Any[1, meh])
    @test ismeh([meh, false] == BitArray([true, false]))
    @test ismeh(Any[meh, false] == BitArray([true, false]))
    @test Union{Int, Meh}[1] == Union{Float64, Meh}[1.0]
    @test Union{Int, Meh}[1] == [1.0]
    @test Union{Bool, Meh}[true] == BitArray([true])
    @test !(Union{Int, Meh}[1] == [2])
    @test !([1] == Union{Int, Meh}[2])
    @test !(Union{Int, Meh}[1] == Union{Int, Meh}[2])

    @test ismeh([1, meh] != [1, meh])
    @test ismeh(["a", meh] != ["a", meh])
    @test ismeh(Any[1, meh] != Any[1, meh])
    @test ismeh(Any[meh] != Any[meh])
    @test ismeh([meh] != [meh])
    @test ismeh(Any[meh, 2] != Any[1, meh])
    @test ismeh([meh, false] != BitArray([true, false]))
    @test ismeh(Any[meh, false] != BitArray([true, false]))
    @test !(Union{Int, Meh}[1] != Union{Float64, Meh}[1.0])
    @test !(Union{Int, Meh}[1] != [1.0])
    @test !(Union{Bool, Meh}[true] != BitArray([true]))
    @test Union{Int, Meh}[1] != [2]
    @test [1] != Union{Int, Meh}[2]
    @test Union{Int, Meh}[1] != Union{Int, Meh}[2]

    @test any([true, meh])
    @test any(x -> x == 1, [1, meh])
    @test ismeh(any([false, meh]))
    @test ismeh(any(x -> x == 1, [2, meh]))
    @test ismeh(all([true, meh]))
    @test ismeh(all(x -> x == 1, [1, meh]))
    @test !all([false, meh])
    @test !all(x -> x == 1, [2, meh])
    @test 1 in [1, meh]
    @test ismeh(2 in [1, meh])
    @test ismeh(meh in [1, meh])

    @test isequal(float([1, meh]), [1, meh])
    @test float([1, meh]) isa Vector{Union{Float64, Meh}}
    @test isequal(float(Union{Int, Meh}[meh]), [meh])
    @test float(Union{Int, Meh}[meh]) isa Vector{Union{Float64, Meh}}
    @test float(Union{Int, Meh}[1]) == [1]
    @test float(Union{Int, Meh}[1]) isa Vector{Union{Float64, Meh}}
    @test isequal(float([meh]), [meh])
    @test float([meh]) isa Vector{Meh}

    # MehException
    @test sprint(showerror, MehException("test")) == "MehException: test"
end


text: # TODO: Test ApproxQuantile


text:
# using Revise, Plots, Test
using BifurcationKit, LinearAlgebra, Setfield, SparseArrays, ForwardDiff, Parameters
const BK = BifurcationKit
norminf = x -> norm(x, Inf)

Fbp(x, p) = [p.a0 + x[1] * (3.23 .* p.μ - p.x2 * x[1] + p.x3 * 0.234 * x[1]^2) + x[2], -x[2]]
par = (a0 = 0., μ = -0.2, ν = 0, x2 = 1.12, x3 = 1.0)
####################################################################################################
opt_newton = NewtonPar(tol = 1e-9, maxIter = 20)
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01, pMax = 0.4, pMin = -0.5, detectBifurcation = 3, nev = 2, newtonOptions = opt_newton, maxSteps = 100, nInversion = 4, tolBisectionEigenvalue = 1e-8, dsminBisection = 1e-9)

br, = continuation(
	Fbp, [0.1, 0.1], par, (@lens _.μ),
	opts_br; plot = false, verbosity = 0, normC = norminf, recordFromSolution = (x, p) -> x[1])

####################################################################################################
# normal form computation
jet = BK.getJet(Fbp, (x, p) -> ForwardDiff.jacobian(z -> Fbp(z, p), x))

bp = BK.computeNormalForm(jet..., br, 1; verbose=false)
@test BK.isTranscritical(bp) == true

jet = BK.getJet(Fbp, (x, p) -> BK.finiteDifferences(z -> Fbp(z, p), x))
bp = BK.computeNormalForm(jet..., br, 1; verbose=false)
@test BK.isTranscritical(bp) == true
show(bp)

# normal form
nf = bp.nf

@test norm(nf.a) < 1e-10
	@test norm(nf.b1 - 3.23) < 1e-10
	@test norm(nf.b2/2 - -1.12) < 1e-10
	@test norm(nf.b3/6 - 0.234) < 1e-10

####################################################################################################
# same but when the eigenvalues are not saved in the branch but computed on the fly
br_noev, = BK.continuation(
	Fbp, [0.1, 0.1], par, (@lens _.μ),
	recordFromSolution = (x, p) -> norminf(x),
	(@set opts_br.saveEigenvectors = false); plot = false, verbosity = 0, normC = norminf)
bp = BK.computeNormalForm(jet..., br_noev, 1; verbose=false)
nf = bp.nf
@test norm(nf[1]) < 1e-10
	@test norm(nf[2] - 3.23) < 1e-10
	@test norm(nf[3]/2 - -1.12) < 1e-10
	@test norm(nf[4]/6 - 0.234) < 1e-10
####################################################################################################
# Automatic branch switching
br2, = continuation(jet..., br, 1, setproperties(opts_br; pMax = 0.2, ds = 0.01, maxSteps = 14); recordFromSolution = (x, p) -> x[1], verbosity = 0)
@test br2 isa Branch
@test BK.haseigenvalues(br2) == true
@test BK.haseigenvector(br2) == true
BK.eigenvals(br2, 1, true)
BK.getfirstusertype(br2)
@test length(br2) == 12
# plot(br,br2)

br3, = continuation(jet..., br, 1, setproperties(opts_br; ds = -0.01); recordFromSolution = (x, p) -> x[1], verbosity = 0, usedeflation = true)
# plot(br,br3)

# automatic bifurcation diagram (Transcritical)
bdiag = bifurcationdiagram(jet..., [0.1, 0.1], par,  (@lens _.μ), 2,
	(args...) -> setproperties(opts_br; pMin = -1.0, pMax = .5, ds = 0.01, dsmax = 0.05, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 30, newtonOptions = (@set opt_newton.verbose=false), maxSteps = 15);
	recordFromSolution = (x, p) -> x[1],
	tangentAlgo = BorderedPred(),
	plot = false, verbosity = 0, normC = norminf)
####################################################################################################
# Case of the pitchfork
par_pf = setproperties(par; x2 = 0.0, x3 = -1.0)
brp, = BK.continuation(
	Fbp, [0.1, 0.1], par_pf, (@lens _.μ),
	recordFromSolution = (x, p) -> x[1],
	opts_br; plot = false, verbosity = 0, normC = norminf)
bpp = BK.computeNormalForm(jet..., brp, 1)
show(bpp)

nf = bpp.nf
@test norm(nf[1]) < 1e-9
	@test norm(nf[2] - 3.23) < 1e-9
	@test norm(nf[3]/2 - 0) < 1e-9
	@test norm(nf[4]/6 + 0.234) < 1e-9

# test automatic branch switching
br2, = continuation(jet..., brp, 1, setproperties(opts_br; maxSteps = 19, dsmax = 0.01, ds = 0.001, detectBifurcation = 2, newtonOptions = (@set opt_newton.verbose=false)); recordFromSolution = (x, p) -> x[1], verbosity = 0, ampfactor = 1)
	# plot(brp,br2, marker=:d)

# test methods for aBS
BK.from(br2) |> BK.type
BK.from(br2) |> BK.istranscritical
BK.type(nothing)
BK.show(stdout, br2)
BK.propertynames(br2)

# automatic bifurcation diagram (Pitchfork)
bdiag = bifurcationdiagram(jet..., [0.1, 0.1], par_pf,  (@lens _.μ), 2,
	(args...) -> setproperties(opts_br; pMin = -1.0, pMax = .5, ds = 0.01, dsmax = 0.05, nInversion = 6, detectBifurcation = 3, maxBisectionSteps = 30, newtonOptions = (@set opt_newton.verbose=false), maxSteps = 15);
	recordFromSolution = (x, p) -> x[1],
	# tangentAlgo = BorderedPred(),
	plot = false, verbosity = 0, normC = norminf)
####################################################################################################
function Fbp2d(x, p)
	return [ p.α * x[1] * (3.23 .* p.μ - 0.123 * x[1]^2 - 0.234 * x[2]^2),
			 p.α * x[2] * (3.23 .* p.μ - 0.456 * x[1]^2 - 0.123 * x[2]^2),
			 -x[3]]
end

jet = BK.getJet(Fbp2d, (x, p) -> ForwardDiff.jacobian(z -> Fbp2d(z, p), x))

par = (μ = -0.2, ν = 0, α = -1)

for α in (1,1)
	global par = @set par.α = α
	br, = BK.continuation(
		Fbp2d, [0.01, 0.01, 0.01], par, (@lens _.μ),
		recordFromSolution = (x, p) -> norminf(x),
		setproperties(opts_br; nInversion = 2); plot = false, verbosity = 0, normC = norminf)
	# we have to be careful to have the same basis as for Fbp2d or the NF will not match Fbp2d
	bp2d = BK.computeNormalForm(jet..., br, 1; ζs = [[1, 0, 0.], [0, 1, 0.]]);
	show(bp2d)

	BK.nf(bp2d)
	length(bp2d)
	bp2d(rand(2), 0.2)
	bp2d(Val(:reducedForm), rand(2), 0.2)

	@test abs(bp2d.nf.b3[1,1,1,1] / 6 - -par.α * 0.123) < 1e-10
	@test abs(bp2d.nf.b3[1,1,2,2] / 2 - -par.α * 0.234) < 1e-10
	@test abs(bp2d.nf.b3[1,1,1,2] / 2 - -par.α * 0.0)   < 1e-10
	@test abs(bp2d.nf.b3[2,1,1,2] / 2 - -par.α * 0.456) < 1e-10
	@test norm(bp2d.nf.b2, Inf) < 3e-6
	@test norm(bp2d.nf.b1 - par.α * 3.23 * I, Inf) < 1e-9
	@test norm(bp2d.nf.a, Inf) < 1e-6
end

##############################
# same but when the eigenvalues are not saved in the branch but computed on the fly instead
br_noev, = BK.continuation(
	Fbp2d, [0.01, 0.01, 0.01], par, (@lens _.μ),
	recordFromSolution = (x, p) -> norminf(x),
	setproperties(opts_br; nInversion = 2, saveEigenvectors = false); plot = false, verbosity = 0, normC = norminf)
bp2d = BK.computeNormalForm(jet..., br_noev, 1; ζs = [[1, 0, 0.], [0, 1, 0.]]);
@test abs(bp2d.nf.b3[1,1,1,1] / 6 - -0.123) < 1e-15
@test abs(bp2d.nf.b3[1,1,2,2] / 2 - -0.234) < 1e-15
@test abs(bp2d.nf.b3[1,1,1,2] / 2 - -0.0)   < 1e-15
@test abs(bp2d.nf.b3[2,1,1,2] / 2 - -0.456) < 1e-15
@test norm(bp2d.nf.b2, Inf) < 3e-15
@test norm(bp2d.nf.b1 - 3.23 * I, Inf) < 1e-9
@test norm(bp2d.nf.a, Inf) < 1e-15
####################################################################################################
# vector field to test nearby secondary bifurcations
FbpSecBif(u, p) = @. -u * (p + u * (2-5u)) * (p -.15 - u * (2+20u))
dFbpSecBif(x,p)=  ForwardDiff.jacobian( z-> FbpSecBif(z,p), x)
jet = BK.getJet(FbpSecBif, dFbpSecBif)

br_snd1, = BK.continuation(
	FbpSecBif, [0.0], -0.2, (@lens _),
	recordFromSolution = (x, p) -> x[1],
	# tangentAlgo = BorderedPred(),
	setproperties(opts_br; pMin = -1.0, pMax = .3, ds = 0.001, dsmax = 0.005, nInversion = 8, detectBifurcation=3); plot = false, verbosity = 0, normC = norminf)

# plot(br_snd1)

br_snd2, = BK.continuation(
	jet..., br_snd1, 1,
	recordFromSolution = (x, p) -> x[1],
	setproperties(opts_br; pMin = -1.2, pMax = 0.2, ds = 0.001, detectBifurcation = 3, maxSteps=19, nInversion = 8, newtonOptions = NewtonPar(opts_br.newtonOptions),dsminBisection =1e-18, tolBisectionEigenvalue=1e-11, maxBisectionSteps=20); plot = false, verbosity = 0, normC = norminf,
	# tangentAlgo = BorderedPred(),
	# finaliseSolution = (z, tau, step, contResult) ->
	# 	(Base.display(contResult.eig[end].eigenvals) ;true)
	)

	# plot(plot(br_snd2.branch[1,:] |> diff, marker = :d),
	# plot(br_snd1,br_snd2, putbifptlegend=false))

bdiag = bifurcationdiagram(jet..., [0.0], -0.2, (@lens _), 2,
	(args...) -> setproperties(opts_br; pMin = -1.0, pMax = .3, ds = 0.001, dsmax = 0.005, nInversion = 8, detectBifurcation = 3,dsminBisection =1e-18, tolBisectionEigenvalue=1e-11, maxBisectionSteps=20, newtonOptions = (@set opt_newton.verbose=false));
	recordFromSolution = (x, p) -> x[1],
	# tangentAlgo = BorderedPred(),
	plot = false, verbosity = 0, normC = norminf)

# plot(bdiag; putbifptlegend=false, markersize=2, plotfold=false, title = "#branch = $(size(bdiag))")

# test calls for aBD
BK.level(bdiag)
BK.hasbranch(bdiag)
BK.from(bdiag.child[1].γ)
BK.getBranchesFromBP(bdiag, 2)
BK.getContResult(br)
BK.getContResult(getBranch(bdiag,(1,)).γ)
size(bdiag)
getBranch(bdiag, (1,))
show(stdout, bdiag)
####################################################################################################
# test of the pitchfork-D6 normal form
function FbpD6(x, p)
	return [ p.μ * x[1] + (p.a * x[2] * x[3] - p.b * x[1]^3 - p.c*(x[2]^2 + x[3]^2) * x[1]),
			 p.μ * x[2] + (p.a * x[1] * x[3] - p.b * x[2]^3 - p.c*(x[3]^2 + x[1]^2) * x[2]),
			 p.μ * x[3] + (p.a * x[1] * x[2] - p.b * x[3]^3 - p.c*(x[2]^2 + x[1]^2) * x[3])]
end
jet = BK.getJet(FbpD6, (x, p) -> ForwardDiff.jacobian(z -> FbpD6(z, p), x))

pard6 = (μ = -0.2, a = 0.3, b = 1.5, c = 2.9)

br, = BK.continuation(
	FbpD6, zeros(3), pard6, (@lens _.μ),
	recordFromSolution = (x, p) -> norminf(x),
	setproperties(opts_br; nInversion = 6, ds = 0.001); plot = false, verbosity = 0, normC = norminf)

# plot(br;  plotfold = false)
# we have to be careful to have the same basis as for Fbp2d or the NF will not match Fbp2d
bp2d = BK.computeNormalForm(jet..., br, 1; ζs = [[1, 0, 0.], [0, 1, 0.], [0, 0, 1.]])
BK.nf(bp2d)

@test abs(bp2d.nf.b3[1,1,1,1] / 6 - -pard6.b) < 1e-10
	@test abs(bp2d.nf.b3[1,1,2,2] / 2 - -pard6.c) < 1e-10
	@test abs(bp2d.nf.b2[1,2,3] - pard6.a)   < 1e-10

# test the evaluation of the normal form
x0 = rand(3); @test norm(FbpD6(x0, set(pard6, br.lens, 0.001))  - bp2d(Val(:reducedForm), x0, 0.001), Inf) < 1e-12

br1, = BK.continuation(
	jet..., br, 1,
	setproperties(opts_br; nInversion = 4, dsmax = 0.005, ds = 0.001, maxSteps = 300, pMax = 1.); plot = false, verbosity = 0, normC = norminf, recordFromSolution = (x, p) -> norminf(x))
	# plot(br1..., br, plotfold=false, putbifptlegend=false)

bdiag = bifurcationdiagram(jet..., zeros(3), pard6, (@lens _.μ), 3,
	(args...) -> setproperties(opts_br; pMin = -0.250, pMax = .4, ds = 0.001, dsmax = 0.005, nInversion = 4, detectBifurcation = 3, dsminBisection =1e-18, tolBisectionEigenvalue=1e-11, maxBisectionSteps=20, newtonOptions = (@set opt_newton.verbose=false));
	recordFromSolution = (x, p) -> norminf(x),
	# tangentAlgo = BorderedPred(),
	plot = false, verbosity = 0, normC = norminf)

# plot(bdiag; putspecialptlegend=false, markersize=2,plotfold=false);title!("#branch = $(size(bdiag))")

####################################################################################################
# test of the Hopf normal form
function Fsl2!(f, u, p, t)
	@unpack r, μ, ν, c3, c5 = p
	u1 = u[1]
	u2 = u[2]
	ua = u1^2 + u2^2

	f[1] = r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1
	f[2] = r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2

	return f
end

Fsl2(x, p) = Fsl2!(similar(x), x, p, 0.)
par_sl = (r = 0.5, μ = 0.132, ν = 1.0, c3 = 1.123, c5 = 0.2)
jet = BK.getJet(Fsl2, (x, p) -> ForwardDiff.jacobian(z -> Fsl2(z, p), x))

# detect hopf bifurcation
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.02, ds = 0.01, pMax = 0.1, pMin = -0.3, detectBifurcation = 3, nev = 2, newtonOptions = (@set opt_newton.verbose = false), maxSteps = 100)

br, = BK.continuation(
	jet[1], jet[2], [0.0, 0.0], (@set par_sl.r = -0.1), (@lens _.r),
	recordFromSolution = (x, p) -> norminf(x),
	opts_br; plot = false, verbosity = 0, normC = norminf)

hp = BK.computeNormalForm(jet..., br, 1)

nf = hp.nf
BK.type(hp)

@test abs(nf.a - 1) < 1e-9
@test abs(nf.b/2 - (-par_sl.c3 + im*par_sl.μ)) < 1e-14

####################################################################################################
# same but when the eigenvalues are not saved in the branch but computed on the fly instead
br, _ = BK.continuation(
	Fsl2, [0.0, 0.0], (@set par_sl.r = -0.1), (@lens _.r),
	recordFromSolution = (x, p) -> norminf(x),
	setproperties(opts_br, saveEigenvectors = false); plot = false, verbosity = 0, normC = norminf)

hp = BK.computeNormalForm(jet..., br, 1)

nf = hp.nf

@test abs(nf.a - 1) < 1e-9
@test abs(nf.b/2 - (-par_sl.c3 + im*par_sl.μ)) < 1e-14

show(hp)
####################################################################################################
# test for the Cusp normal form
Fcusp(x, p) = [p.β1 + p.β2 * x[1] + p.c * x[1]^3]
par = (β1 = 0.0, β2 = -0.01, c = 3.)
jet  = BK.getJet(Fcusp; matrixfree=false)
br, = continuation(jet[1], jet[2], [0.01], par, (@lens _.β1), opts_br;)

sn_codim2, = continuation(jet[1:2]..., br, 1, (@lens _.β2), ContinuationPar(opts_br, detectBifurcation = 1, saveSolEveryStep = 1, maxSteps = 40) ;
	updateMinAugEveryStep = 1,
	d2F = jet[3], d3F = jet[4],
	bdlinsolver = MatrixBLS()
	)
# find the cusp point
ind = findall(map(x->x.type == :cusp, sn_codim2.specialpoint))
cuspnf = computeNormalForm(jet..., sn_codim2, ind[1])
show(cuspnf)
BK.type(cuspnf)
@test cuspnf.nf.c == par.c
####################################################################################################
# test for the Bogdanov-Takens normal form
Fbt(x, p) = [x[2], p.β1 + p.β2 * x[2] + p.a * x[1]^2 + p.b * x[1] * x[2]]
par = (β1 = 0.01, β2 = -0.1, a = -1., b = 1.)
jet  = BK.getJet(Fbt; matrixfree=false)
opt_newton = NewtonPar(tol = 1e-9, maxIter = 40, verbose = false)
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.01, ds = 0.01, pMax = 0.5, pMin = -0.5, detectBifurcation = 3, nev = 2, newtonOptions = opt_newton, maxSteps = 100, nInversion = 8, tolBisectionEigenvalue = 1e-8, dsminBisection = 1e-9, saveSolEveryStep = 1)

br, = continuation(jet[1], jet[2], [0.01, 0.01], par, (@lens _.β1), opts_br; bothside = true)

sn_codim2, = continuation(jet[1:2]..., br, 1, (@lens _.β2), ContinuationPar(opts_br, detectBifurcation = 1, saveSolEveryStep = 1, maxSteps = 40) ;
	detectCodim2Bifurcation = 2,
	updateMinAugEveryStep = 1,
	d2F = jet[3], d3F = jet[4],
	bdlinsolver = MatrixBLS()
	)

hopf_codim2, = continuation(jet[1:2]..., br, 2, (@lens _.β2), ContinuationPar(opts_br, detectBifurcation = 1, saveSolEveryStep = 1, maxSteps = 40) ;
	detectCodim2Bifurcation = 2,
	updateMinAugEveryStep = 1,
	bothside = true,
	d2F = jet[3], d3F = jet[4],
	bdlinsolver = MatrixBLS(),
	)

# plot(sn_codim2, hopf_codim2, branchlabel = ["Fold", "Hopf"])

btpt = computeNormalForm(jet..., sn_codim2, 1; nev = 2)
show(btpt)
BK.type(btpt)
@test norm(btpt.nf.b * sign(sum(btpt.ζ[1])) - par.b, Inf) < 1e-5
@test norm(btpt.nf.a * sign(sum(btpt.ζ[1])) - par.a, Inf) < 1e-5
@test isapprox(abs.(btpt.ζ[1]), [1, 0])
@test isapprox(abs.(btpt.ζ[2]), [0, 1];rtol = 1e-6)
@test isapprox(abs.(btpt.ζstar[1]), [1, 0];rtol = 1e-6)

@test isapprox(btpt.nfsupp.K2, [0, 0]; atol = 1e-5)
@test isapprox(btpt.nfsupp.d, 0; atol = 1e-3)
@test isapprox(btpt.nfsupp.e, 0; atol = 1e-3)
@test isapprox(btpt.nfsupp.a1, 0; atol = 1e-3)
@test isapprox(btpt.nfsupp.b1, 0; atol = 1e-3)

HC = BK.predictor(btpt, Val(:HopfCurve), 0.)
	HC.hopf(0.)
SN = BK.predictor(btpt, Val(:FoldCurve), 0.)
Hom = BK.predictor(btpt, Val(:HomoclinicCurve), 0.)
	Hom.orbit(0,0)

# plot(sn_codim2, branchlabel = ["Fold"], vars = (:β1, :β2))
# 	_S = LinRange(-0.06, 0.06, 1000)
# 	plot!([HC.hopf(s)[1] for s in _S], [HC.hopf(s)[2] for s in _S], linewidth=5, label = "Hpred")
# 	plot!([SN.fold(s)[1] for s in _S], [SN.fold(s)[2] for s in _S], linewidth=5, label = "SNpred")
# 	_S = LinRange(-0.25, 0.25, 1000)
# 	plot!([Hom.α(s)[1] for s in _S], [Hom.α(s)[2] for s in _S], linewidth=5, label = "Hom")
#
# 	plot!(hopf_codim2, branchlabel = ["Hopf"], vars = (:β1, :β2), color = :black)
# 	xlims!(-0.001, 0.05)

# plot of the homoclinic orbit
hom1 = [Hom.orbit(t,0.1)[1] for t in LinRange(-1000, 1000, 10000)]
hom2 = [Hom.orbit(t,0.1)[2] for t in LinRange(-1000, 1000, 10000)]
# plot(hom1, hom2)

# branch switching from BT from Fold
opt = sn_codim2.contparams
@set! opt.newtonOptions.verbose = false
@set! opt.maxSteps = 20
hp_fromBT, = continuation(jet..., sn_codim2, 1, opt;
	verbosity = 0, plot = false,
	δp = 1e-4,
	updateMinAugEveryStep = 1,
	)

####################################################################################################
# test of the Bautin normal form
function Fsl2!(f, u, p, t)
	@unpack r, μ, ν, c3, c5 = p
	u1, u2 = u
	ua = u1^2 + u2^2
	f[1] = r * u1 - ν * u2 + ua * (c3 * u1 - μ * u2) + c5 * ua^2 * u1
	f[2] = r * u2 + ν * u1 + ua * (c3 * u2 + μ * u1) + c5 * ua^2 * u2
	return f
end

Fsl2(x, p) = Fsl2!(similar(x), x, p, 0.)
par_sl = (r = -0.5, μ = 0., ν = 1.0, c3 = 0.1, c5 = 0.3)
jet = BK.getJet(Fsl2, matrixfree=false)

@set! opts_br.newtonOptions.verbose = false
@set! opts_br.newtonOptions.tol = 1e-12
opts_br = setproperties(opts_br;nInversion = 10, maxBisectionSteps = 25)

br, = continuation( jet[1], jet[2], [0.01, 0.01], par_sl, (@lens _.r), opts_br)

hopf_codim2, = continuation(jet[1:2]..., br, 1, (@lens _.c3), ContinuationPar(opts_br, detectBifurcation = 0, saveSolEveryStep = 1, maxSteps = 40, pMin = -2., pMax = 2., ds = -0.001) ;
	detectCodim2Bifurcation = 2,
	startWithEigen = true,
	updateMinAugEveryStep = 1,
	d2F = jet[3], d3F = jet[4],
	bdlinsolver = MatrixBLS(),
	)

bautin = BifurcationKit.computeNormalForm(jet..., hopf_codim2, 1; nev = 2)
show(bautin)
BK.type(bautin)

@test bautin.nf.l2 ≈ par_sl.c5 * 4


text:
using CSV
using DataFrames
using Colors
using PyPlot

const rexeltype = r"N0f8|N0f16|Float32"
const idxeltype = Dict("N0f8"=>1, "N0f16"=>2, "Float32"=>3)
const eltypeidx = Dict(v=>k for (k, v) in idxeltype)
const rexsize = r"_(\d+)$"
const special_replacements = ["components" => "label_components", "dblcone" => "distance_transform", "spiral" => "flood"]

randomtag(fn) = (occursin(r"rgb", fn), idxeltype[match(rexeltype, fn).match], parse(Int, only(match(rexsize, fn).captures)))
randomlabel(tag) = (tag[1] ? "RGB" : "Gray") * "{" * eltypeidx[tag[2]] * "}, $(tag[3])"

function imagelt(image1, image2)
    image1 ∈ ("2d", "3d") && image2 ∈ ("2d", "3d") && return image1 < image2

    getcolor(str) = match(r"(Gray|RGB)", str).captures[1]::AbstractString
    geteltype(str) = idxeltype[match(r"(N0f8|N0f16|Float32)", str).captures[1]]
    getsz(str) = parse(Int, match(r", (\d)+", str).captures[1])

    col1, col2 = getcolor(image1), getcolor(image2)
    col1 != col2 && return col1 < col2
    et1, et2 = geteltype(image1), geteltype(image2)
    et1 != et2 && return et1 < et2
    return getsz(image1) < getsz(image2)
end

hexhash(c) = "#"*hex(c)
const language_colors = (cols = hexhash.(distinguishable_colors(4, [colorant"white", colorant"black"]; dropseed=true));
                         Dict("python-skimage"=>cols[1], "Matlab"=>cols[2], "OpenCV"=>cols[3], "Julia"=>cols[4]))

function taskplots(jobtag::Function, jobname::Function, filelabel::Pair...; colordict=language_colors, tasknames=nothing, pervoxel::Bool=true)
    # Read all the data
    benchdata = Dict{String,GroupedDataFrame{DataFrame}}()
    for (file, label) in filelabel
        benchdata[label] = groupby(DataFrame(CSV.File(file; delim=',')), "Benchmark")
    end
    # Ensure all inputs have the same tasks
    ks = keys(first(benchdata).second)
    knames = Set(only(Tuple(k)) for k in ks)
    for (_, df) in benchdata
        dfknames = Set(only(Tuple(k)) for k in keys(df))
        # kdf ⊆ ks && ks ⊆ kdf || error("disjoint tasks")
        dfknames == knames || error("disjoint tasks")
    end
    # Determine the plotting order
    if tasknames === nothing
        tasknames = sort(collect(knames))
    end
    # Assign colors to each suite
    if colordict === nothing
        colordict = Dict(zip(keys(benchdata), hexhash.(distinguishable_colors(length(benchdata), [colorant"white"]; dropseed=true))))
    end
    ncols = floor(Int, sqrt(length(tasknames)))
    nrows = ceil(Int, length(tasknames)/ncols)
    fig, axes = plt.subplots(nrows=nrows, ncols=ncols, gridspec_kw=Dict("bottom"=>0.2))
    for i = length(tasknames)+1:length(axes)
        axes[i].set_axis_off()  # close()
    end
    legend_labels, hlines = String[], []
    for (axidx, ax, k) in zip(CartesianIndices(axes), axes, tasknames)
        xlbls = nothing
        for (label, data) in benchdata
            data = sort(data[(Benchmark=k,)], "File"; by=jobtag)
            if xlbls === nothing
                xlbls = jobname.(jobtag.(data."File"))
            end
            y = data."Time(s)"
            if pervoxel
                nvox = map(data."File") do fn
                    jobtag(fn)[end]^2
                end
                y ./= nvox
            end
            hline = ax.plot(y; color=colordict[label], marker="x")
            if length(legend_labels) < length(filelabel)
                push!(legend_labels, label)
                push!(hlines, first(hline))
            end
        end
        if axidx[1] == size(axes, 1) || LinearIndices(axes)[axidx] == length(tasknames)
            ax.set_xticks(0:length(xlbls)-1)
            ax.set_xticklabels(xlbls; rotation=90)
            ax.set_xlabel("Jobs")
        else
            ax.set_xticks([])
        end
        ax.set_yscale("log")
        if axidx[2] == 1
            ax.set_ylabel("Time (s)")
        end
        ax.set_title(k)
    end
    fig.legend(hlines, legend_labels)
    return fig
end

# taskplots(randomtag, randomlabel, "julia/julia_generics.csv" => "Julia",
#                                   "python-skimage/python_generics.csv" => "python-skimage",
#                                   "matlab/matlab_generics.csv" => "Matlab",
#                                   "opencv/opencv_generics.csv" => "OpenCV",
# )

function specialplots(filelabel::Pair...; colordict=language_colors)
    # Read all the data
    benchdata = Dict{String,DataFrame}()
    for (file, label) in filelabel
        benchdata[label] = DataFrame(CSV.File(file; delim=','))
    end
    # Ensure all inputs have the same tasks
    tasknames = sort(first(benchdata).second."Benchmark")
    all(df -> Set(df."Benchmark") == Set(tasknames), values(benchdata))
    # Assign colors to each suite
    if colordict === nothing
        colordict = Dict(zip(keys(benchdata), hexhash.(distinguishable_colors(length(benchdata), [colorant"white"]; dropseed=true))))
    end
    fig, ax = plt.subplots()
    x = 1:length(tasknames)
    w = 1/(2*length(benchdata))
    for (i, (label, df)) in enumerate(benchdata)
        p = Int.(indexin(df."Benchmark", tasknames))
        dfs = df[invperm(p),:]
        display(dfs)
        ax.bar(x .+ (i-1)*w .- w/2, dfs."Time(s)", w; label=label, color=colordict[label])
    end
    ax.set_xticks(x)
    repnames = map(tasknames) do tn
        foldl((str, rep) -> replace(str, rep), special_replacements; init=tn)
    end
    ax.set_xticklabels(repnames; rotation=90)
    ax.set_yscale("log")
    ax.set_ylabel("Time (s)")
    fig.legend()
    fig.tight_layout()
    return fig
end

# specialplots("julia/julia_special.csv" => "Julia",
#              "python-skimage/python_special.csv" => "python-skimage",
#              "matlab/matlab_special.csv" => "Matlab",
#              #= "opencv/opencv_special.csv" => "OpenCV", =#
# )

function combine_all(prs...; szs=(r"64_64", r"2048_2048"), eltypes=(r"gray_N0f8", r"gray_Float32", r"rgb_N0f8", r"rgb_Float32"))
    df = DataFrame("Framework"=>String[], "Operation"=>String[], "Image"=>String[], "Time"=>Float64[])
    for (file, framework) in prs
        for row in CSV.Rows(file)
            if length(row) == 3  # one of the generics
                any(occursin(row.File), szs) || continue
                any(occursin(row.File), eltypes) || continue
                tag = randomtag(row.File)
                t = getproperty(row, Symbol("Time(s)"))
                push!(df, (framework, row.Benchmark, randomlabel(tag), parse(Float64, t)))
            else  # one of the special
                imagetag = occursin("3d", row.Benchmark) ? "3d" : "2d"
                optag = foldl((str, rep) -> replace(str, rep), special_replacements; init=row.Benchmark[1:end-2])
                t = getproperty(row, Symbol("Time(s)"))
                push!(df, (framework, optag, imagetag, parse(Float64, t)))
            end
        end
    end
    return df
end

# df = combine_all("julia/julia_generics.csv" => "Julia",
#                  "python-skimage/python_generics.csv" => "python-skimage",
#                  "matlab/matlab_generics.csv" => "Matlab",
#                  "opencv/opencv_generics.csv" => "OpenCV",
#                  "julia/julia_special.csv" => "Julia",
#                  "python-skimage/python_special.csv" => "python-skimage",
#                  "matlab/matlab_special.csv" => "Matlab")

function ratio_df(df; operation_order=["complement", "mean", "histeq", "blur", "gradient", "label_components", "flood", "distance_transform"])
    jdata = filter(:Framework=>isequal("Julia"), df)
    jgroups = groupby(jdata, :Operation)
    rdf = DataFrame("Framework"=>String[], "Operation"=>String[], "Image"=>String[], "Ratio"=>Float64[])
    for f in ("python-skimage", "Matlab", "OpenCV")
        fdata = filter(:Framework=>isequal(f), df)
        fgroups = groupby(fdata, :Operation)
        for op in operation_order
            k = (Operation=op,)
            jref = sort(jgroups[k], :Image; lt=imagelt)
            try
                fval = sort(fgroups[k], :Image; lt=imagelt)
                for (jrow, frow) in zip(eachrow(jref), eachrow(fval))
                    @assert jrow.Image == frow.Image
                    push!(rdf, [f, op, frow.Image, frow.Time / jrow.Time])
                end
            catch
                for jrow in eachrow(jref)
                    push!(rdf, [f, op, jrow.Image, NaN])
                end
            end
        end
    end
    return rdf
end

function plot_ratio(rdf; colordict=language_colors)
    frameworks = unique(rdf.Framework)
    operations = unique(rdf.Operation)
    n = Int(nrow(rdf)/length(frameworks))  # deliberately errors if not an even divisor
    x = 1:n+length(operations)-1
    xlabels = String[]
    opgaps = Int[]
    gdf = groupby(rdf, :Framework)
    fig, axouter = plt.subplots(figsize=(7.5,4), subplot_kw=Dict("position" => [0.1, 0.05, 0.9, 0.95]))
    insetgap = 0.4
    ax = axouter.inset_axes([0, insetgap, 1, 1-insetgap])
    for k in keys(gdf)
        g = gdf[k]
        col = colordict[k.Framework]
        oldop = ""
        y = Float64[]
        for (i, (op, t)) in enumerate(zip(g.Operation, g.Ratio))
            if op != oldop && i != 1
                push!(y, NaN)
                if length(xlabels) < n && i != 1
                    # push!(xlabels, "")
                    push!(opgaps, length(y))
                end
            end
            oldop = op
            push!(y, t)
            if length(xlabels) < n
                push!(xlabels, g.Image[i])
            end
        end

        ax.plot(x, y, color=col, label=k.Framework, marker="x")
    end
    ax.plot(x, fill(1, length(x)), "k--", label=nothing)
    ax.set_xticks(setdiff(x, opgaps))
    ax.set_xticklabels(xlabels; rotation=90, fontsize=6)
    ax.set_yscale(:log)
    ax.set_ylim((1/3000, 3000))
    ax.set_yticks([1/1000, 1/100, 1/10, 1, 10, 100, 1000])
    ax.set_ylabel("Time relative to Julia")
    ax.legend(loc="lower left", fontsize=8)
    pushfirst!(opgaps, 0)
    push!(opgaps, length(x)+1)
    xc = (opgaps[1:end-1] .+ opgaps[2:end])./2
    axouter.set_xlim(ax.get_xlim())
    axouter.set_xticks(xc)
    axouter.set_xticklabels(operations; rotation=90, fontsize=9)
    for obj in axouter.get_xticklines()
        obj.set_visible(false)
    end
    axouter.get_yaxis().set_visible(false)
    axouter.set_frame_on(false)
    # for spkey in axouter.spines
    #     axouter.spines[spkey].set_visible(false)
    # end
    # for sp in ("left", "bottom", "right", "top")
    #     axouter.spines[sp].set_visible(false)
    # end
    fig.tight_layout()
    return fig
end


text:
"""
MulAddAccumulate(μ, A, B)

represents the vector satisfying v[k+1] == A[k]*v[k]+B[k] with v[1] == μ
"""
mutable struct MulAddAccumulate{T} <: AbstractCachedVector{T}
    data::Vector{T}
    A::AbstractVector
    B::AbstractVector
    datasize::Tuple{Int}
    function MulAddAccumulate{T}(data::Vector{T}, A::AbstractVector, B::AbstractVector, datasize::Tuple{Int}) where T
        length(A) == length(B) || throw(ArgumentError("lengths must match"))
        new{T}(data, A, B, datasize)
    end
end

size(M::MulAddAccumulate) = size(M.A)

MulAddAccumulate(data::Vector{T}, A::AbstractVector, B::AbstractVector, datasize::Tuple{Int}) where T =
    MulAddAccumulate{T}(data, A, B, datasize)
MulAddAccumulate(μ, A, B) = MulAddAccumulate([μ], A, B, (1,))
MulAddAccumulate(A, B) = MulAddAccumulate(A[1]+B[1], A, B)

function LazyArrays.cache_filldata!(K::MulAddAccumulate, inds)
    A,B = K.A,K.B
    @inbounds for k in inds
        K.data[k] = muladd(A[k], K.data[k-1], B[k])
    end
end

"""
    divmul(A, B, C)

is equivalent to A \\ (B*C)
"""
function divmul(Q::SemiclassicalJacobi, D::Derivative, P::SemiclassicalJacobi)
    A,B,_ = recurrencecoefficients(P)
    α,β,_ = recurrencecoefficients(Q)

    d = AccumulateAbstractVector(*, A ./ Vcat(1,α))
    v1 = AccumulateAbstractVector(+, B ./ A)
    v2 = MulAddAccumulate(Vcat(0,0,α[2:∞]) ./ α, Vcat(0,β ./ α) ./ α);
    v3 = AccumulateAbstractVector(*, Vcat(A[1]A[2], A[3:∞] ./ α))
    _BandedMatrix(Vcat(((1:∞) .* d)', (((1:∞) .* (v1 .+ B[2:end]./A[2:end]) .- (2:∞) .* (α .* v2 .+ β ./ α)) .* v3)'), ∞, 2,-1)'
end

@simplify function *(D::Derivative, P::SemiclassicalJacobi)
    Q = SemiclassicalJacobi(P.t, P.a+1,P.b+1,P.c+1,P)
    Q * divmul(Q, D, P)
end

function divmul(wP::Weighted{<:Any,<:SemiclassicalJacobi}, D::Derivative, wQ::Weighted{<:Any,<:SemiclassicalJacobi})
    Q,P = wQ.P,wP.P
    ((-sum(orthogonalityweight(Q))/sum(orthogonalityweight(P))) * (Q \ (D * P))')
end

@simplify function *(D::Derivative, wQ::Weighted{<:Any,<:SemiclassicalJacobi})
    wP = Weighted(SemiclassicalJacobi(wQ.P.t, wQ.P.a-1,wQ.P.b-1,wQ.P.c-1))
    wP * divmul(wP, D, wQ)
end

##
# One-Weighted
##

function divmul(HQ::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi}, D::Derivative, HP::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi})
    Q = HQ.P
    P = HP.P
    t = P.t
    a = Q.a
    A,B,C = recurrencecoefficients(P)
    α,β,γ = recurrencecoefficients(Q)
    d = AccumulateAbstractVector(*, A ./ α)
    v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
    v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))

    _BandedMatrix(
        Vcat(
        ((a:∞) .* v2 .- ((a+1):∞) .* Vcat(1,v1[2:end] .* d))',
        (((a+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
end

function divmul(HQ::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi}, D::Derivative, HP::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi})
    Q = HQ.P
    P = HP.P
    t = P.t
    b = Q.b
    A,B,C = recurrencecoefficients(P)
    α,β,γ = recurrencecoefficients(Q)
    d = AccumulateAbstractVector(*, A ./ α)
    d2 = AccumulateAbstractVector(*, A ./ Vcat(1,α))
    v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
    v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))

    _BandedMatrix(
        Vcat(
        (-(b:∞) .* v2 .+ ((b+1):∞) .* Vcat(1,v1[2:end] .* d) .+ Vcat(0,(1:∞) .* d2))',
        (-((b+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
end

function divmul(HQ::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi}, D::Derivative, HP::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi})
    Q = HQ.P
    P = HP.P
    t = P.t
    c = Q.c
    A,B,C = recurrencecoefficients(P)
    α,β,γ = recurrencecoefficients(Q)
    d = AccumulateAbstractVector(*, A ./ α)
    d2 = AccumulateAbstractVector(*, A ./ Vcat(1,α))
    v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
    v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
    _BandedMatrix(
        Vcat(
        (-(c:∞) .* v2 .+ ((c+1):∞) .* Vcat(1,v1[2:end] .* d) .+ Vcat(0,(t:t:∞) .* d2))',
        (-((c+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
end

@simplify function *(D::Derivative, HP::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi})
    P = HP.P
    t = P.t
    HQ = HalfWeighted{:a}(SemiclassicalJacobi(t, P.a-1, P.b+1, P.c+1))
    HQ * divmul(HQ, D, HP)
end

@simplify function *(D::Derivative, HP::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi})
    P = HP.P
    t = P.t
    HQ = HalfWeighted{:b}(SemiclassicalJacobi(t, P.a+1, P.b-1, P.c+1))
    HQ * divmul(HQ, D, HP)
end

@simplify function *(D::Derivative, HP::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi})
    P = HP.P
    t = P.t
    HQ = HalfWeighted{:c}(SemiclassicalJacobi(t, P.a+1, P.b+1, P.c-1))
    HQ * divmul(HQ, D, HP)
end

##
# Double-Weighted
##

function divmul(HQ::HalfWeighted{:ab}, D::Derivative, HP::HalfWeighted{:ab})
    Q = HQ.P
    P = HP.P
    A,B,_ = recurrencecoefficients(P)
    α,β,_ = recurrencecoefficients(Q)
    a,b = Q.a,Q.b

    d = AccumulateAbstractVector(*, Vcat(1,A) ./ α)
    e = AccumulateAbstractVector(*, Vcat(1,A ./ α))
    f = MulAddAccumulate(Vcat(0,0,A[2:end] ./ α[2:end]), Vcat(0, (B./ α) .* e))
    g = cumsum(β ./ α)
    _BandedMatrix(Vcat((((a+1):∞) .* e .- ((b+a+1):∞).*f .+ ((a+b+2):∞) .* e .* g )',
                           (-((a+b+2):∞)  .* d)'),ℵ₀,1,0)
end

function divmul(HQ::HalfWeighted{:bc}, D::Derivative, HP::HalfWeighted{:bc})
    Q = HQ.P
    P = HP.P
    t = P.t
    A,B,_ = recurrencecoefficients(P)
    α,β,_ = recurrencecoefficients(Q)
    b,c = Q.b,Q.c

    d = AccumulateAbstractVector(*, Vcat(1,A) ./ α)
    e = AccumulateAbstractVector(*, Vcat(1,A ./ α))
    f = MulAddAccumulate(Vcat(0,0,A[2:end] ./ α[2:end]), Vcat(0, (B./ α) .* e))
    g = cumsum(β ./ α)
    _BandedMatrix(Vcat((-((t+1)* (0:∞) .+ (t*(b+1) + c+1)) .* e .+ ((c+b+1):∞).*f .- ((b+c+2):∞) .* e .* g )',
                        (((b+c+2):∞)  .* d)'),ℵ₀,1,0)
end

function divmul(HQ::HalfWeighted{:ac}, D::Derivative, HP::HalfWeighted{:ac})
    Q = HQ.P
    P = HP.P
    t = P.t
    A,B,_ = recurrencecoefficients(P)
    α,β,_ = recurrencecoefficients(Q)
    a,c = Q.a,Q.c

    d = AccumulateAbstractVector(*, Vcat(1,A) ./ α)
    e = AccumulateAbstractVector(*, Vcat(1,A ./ α))
    f = MulAddAccumulate(Vcat(0,0,A[2:end] ./ α[2:end]), Vcat(0, (B./ α) .* e))
    g = cumsum(β ./ α)
    _BandedMatrix(Vcat((t* ((a+1):∞) .* e .- ((c+a+1):∞).*f .+ ((a+c+2):∞) .* e .* g )',
            (-((a+c+2):∞)  .* d)'),ℵ₀,1,0)
end

@simplify function *(D::Derivative, HP::HalfWeighted{:ab,<:Any,<:SemiclassicalJacobi})
    P = HP.P
    t = P.t
    HQ = HalfWeighted{:ab}(SemiclassicalJacobi(t, P.a-1,P.b-1,P.c+1))
    HQ * divmul(HQ, D, HP)
end

@simplify function *(D::Derivative, HP::HalfWeighted{:ac,<:Any,<:SemiclassicalJacobi})
    P = HP.P
    t = P.t
    HQ = HalfWeighted{:ac}(SemiclassicalJacobi(t, P.a-1,P.b+1,P.c-1))
    HQ * divmul(HQ, D, HP)
end

@simplify function *(D::Derivative, HP::HalfWeighted{:bc,<:Any,<:SemiclassicalJacobi})
    P = HP.P
    t = P.t
    HQ = HalfWeighted{:bc}(SemiclassicalJacobi(t, P.a+1,P.b-1,P.c-1))
    HQ * divmul(HQ, D, HP)
end


text:
function has_ext_finalization()
    @use rt = Runtime() begin
        supported = Ref{Cint}(0)

        # query runtime for support for the extension
        err = ccall(
            (:hsa_system_extension_supported,libhsa),
            hsa_status_t,
            (UInt16,UInt16,UInt16,Ptr{Cint}),
            ExtensionFinalizer,
            UInt16(1),
            UInt16(0),
            supported)

        if test_status(err, warn_only = true) && supported.x != 0
            return true
        else
            return false
        end
    end
end


text:
using MLStyle
using SimpleGraphs
using SimplePosets

expr(x) = :(identity($x))

# like `something`, but doesn't throw an error
maybesomething() = nothing
maybesomething(x::Nothing, y...) = maybesomething(y...)
maybesomething(x::Some, y...) = x.value
maybesomething(x::Any, y...) = x

export argtuple
argtuple(m) = arguments(m) |> astuple

astuple(x) = Expr(:tuple,x...)
astuple(x::Symbol) = Expr(:tuple,x)

export arguments
arguments(m::Model) = m.args
arguments(d::JointDistribution) = d.args

export sampled
sampled(m::Model) = keys(m.dists) |> collect

export assigned
assigned(m::Model) = keys(m.vals) |> collect

export parameters
parameters(m::Model) = union(assigned(m), sampled(m))

export variables
variables(m::Model) = union(arguments(m), parameters(m))

function variables(expr :: Expr)
    leaf(x::Symbol) = begin
        [x]
    end
    leaf(x) = []
    branch(head, newargs) = begin
        union(newargs...)
    end
    foldall(leaf, branch)(expr)
end

variables(s::Symbol) = [s]
variables(x) = []

for f in [:arguments, :assigned, :sampled, :parameters, :variables]
    @eval function $f(m::Model, nt::NamedTuple)
        vs = $f(m)
        isempty(vs) && return NamedTuple()
        return select(nt, $f(m))
    end
end

export foldall
function foldall(leaf, branch; kwargs...)
    function go(ast)
        MLStyle.@match ast begin
            Expr(head, args...) => branch(head, map(go, args); kwargs...)
            x                   => leaf(x; kwargs...)
        end
    end

    return go
end

export foldall1
function foldall1(leaf, branch; kwargs...)
    function go(ast)
        MLStyle.@match ast begin
            Expr(head, arg1, args...) => branch(head, arg1, map(go, args); kwargs...)
            x                         => leaf(x; kwargs...)
        end
    end

    return go
end

import MacroTools: striplines, @q

# function arguments(model::Model)
#     model.args
# end

allequal(xs) = all(xs[1] .== xs)

# # fold example usage:
# # ------------------
# # function leafCount(ast)
# #     leaf(x) = 1
# #     expr(head, arg1, newargs) = sum(newargs)
# #     fold(leaf, expr)(ast)
# # end

# # leaves = begin
# #     leaf(x) = [x]
# #     expr(head, arg1, newargs) = union(newargs...)
# #     fold(leaf, expr)
# # end

# # ast = :(f(x + 3y))

# # leaves(ast)

# # Example of Tamas Papp's `as` combinator:
# # julia> as((;s=as(Array, as𝕀,4), a=asℝ))(randn(5))
# # (s = [0.545324, 0.281332, 0.418541, 0.485946], a = 2.217762640580984)

function buildSource(m, proc, wrap=identity; kwargs...)

    kernel = @q begin end

    for st in map(v -> findStatement(m,v), toposort(m))
        ex = proc(m, st; kwargs...)
        isnothing(ex) || push!(kernel.args, ex)
    end

    isnothing(m.retn) || push!(kernel.args, proc(m, Return(m.retn); kwargs...))
    # args = argtuple(m)

    # body = @q begin
    #     function $basename($args; kwargs...)
    #         $(wrap(kernel))
    #     end
    # end

    wrap(kernel) |> flatten
    # flatten(body)
end

# From https://github.com/thautwarm/MLStyle.jl/issues/66
@active LamExpr(x) begin
           @match x begin
               :($a -> begin $(bs...) end) =>
                 let exprs = filter(x -> !(x isa LineNumberNode), bs)
                   if length(exprs) == 1
                     (a, exprs[1])
                   else
                     (a, Expr(:block, bs...))
                     end
               end
                _  => nothing
           end
       end

# using BenchmarkTools
# f(;kwargs...) = kwargs[:a] + kwargs[:b]

# @btime invokefrozen(f, Int; a=3,b=4)  # 3.466 ns (0 allocations: 0 bytes)
# @btime f(;a=3,b=4)                    # 1.152 ns (0 allocations: 0 bytes)

# @isdefined
# Base.@locals
# @__MODULE__
# names

# getprototype(::Type{NamedTuple{(),Tuple{}}}) = NamedTuple()
getprototype(::Type{NamedTuple{N,T} where {T <: Tuple} } ) where {N} = NamedTuple{N}
getprototype(::NamedTuple{N,T} where {T<: Tuple} ) where N = NamedTuple{N}

function loadvals(argstype, datatype)
    args = getntkeys(argstype)
    data = getntkeys(datatype)
    loader = @q begin
    end

    for k in args
        push!(loader.args, :($k = _args.$k))
    end
    for k in data
        push!(loader.args, :($k = _data.$k))
    end

    src -> (@q begin
        $loader
        $src
    end) |> flatten
end

function loadvals(argstype, datatype, parstype)
    args = getntkeys(argstype)
    data = getntkeys(datatype)
    pars = getntkeys(parstype)

    loader = @q begin

    end

    for k in args
        push!(loader.args, :($k = _args.$k))
    end
    for k in data
        push!(loader.args, :($k = _data.$k))
    end

    for k in pars
        push!(loader.args, :($k = _pars.$k))
    end

    src -> (@q begin
        $loader
        $src
    end) |> flatten
end

getntkeys(::NamedTuple{A,B}) where {A,B} = A
getntkeys(::Type{NamedTuple{A,B}}) where {A,B} = A

# These macros quickly define additional methods for when you get tired of typing `NamedTuple()`
macro tuple3args(f)
    quote
        $f(m::Model, (), data) = $f(m::Model, NamedTuple(), data)
        $f(m::Model, args, ()) = $f(m::Model, args, NamedTuple())
        $f(m::Model, (), ())   = $f(m::Model, NamedTuple(), NamedTuple())
    end
end

macro tuple2args(f)
    quote
        $f(m::Model, ()) = $f(m::Model, NamedTuple())
    end
end

# This is just handy for REPLing, no direct connection to Soss

# julia> tower(Int)
# 6-element Array{DataType,1}:
#  Int64
#  Signed
#  Integer
#  Real
#  Number
#  Any

export tower

function tower(x)
    t0 = typeof(x)
    result = [t0]
    t1 = supertype(t0)
    while t1 ≠ t0
        push!(result, t1)
        t0, t1 = t1, supertype(t1)
    end
    return result
end

const TypeLevel = GeneralizedGenerated.TypeLevel

unVal(::Type{Val{T}}) where {T} = T
unVal(::Val{T}) where {T} = T


text:
using Gen

using Cairo
using FileIO
using ImageMagick
using Compat
using Compat.Base64
using Images
using ImageFiltering

const letters = String["A", "B", "C"]

###########################
# primitive distributions #
###########################

struct NoisyMatrix <: Distribution{Matrix{Float64}} end
const noisy_matrix = NoisyMatrix()

function Gen.logpdf(::NoisyMatrix, x::Matrix{Float64}, mu::Matrix{Float64}, noise::Float64)
    var = noise * noise
    diff = x - mu
    vec = diff[:]
    -(vec' * vec)/ (2.0 * var) - 0.5 * log(2.0 * pi * var)
end

function Gen.random(::NoisyMatrix, mu::Matrix{Float64}, noise::Float64)
    mat = copy(mu)
    (w, h) = size(mu)
    for i=1:w
        for j=1:h
            mat[i, j] = mu[i, j] + randn() * noise
        end
    end
    mat
end

struct Object
    x::Float64
    y::Float64
    angle::Float64
    fontsize::Int
    letter::String
end

const width = 40
const height = 40
const min_size = 5
const max_size = 10

function render(obj::Object)
    canvas = CairoRGBSurface(width, height)
    cr = CairoContext(canvas)
    Cairo.save(cr)

    # set background color to white
    set_source_rgb(cr, 1.0, 1.0, 1.0)
    rectangle(cr, 0.0, 0.0, width, height)
    fill(cr)
    restore(cr)
    Cairo.save(cr)

    # write some letters
    set_font_face(cr, "Sans $(obj.fontsize)")
    text(cr, obj.x, obj.y, obj.letter, angle=obj.angle)

    # convert to matrix of color values
    buf = IOBuffer()
    write_to_png(canvas, buf)
    Images.Gray.(readblob(take!(buf)))
end

@gen function model()

    # object prior
    x = @addr(uniform_continuous(0, 1), "x")
    y = @addr(uniform_continuous(0, 1), "y")
    size = @addr(uniform_continuous(0, 1), "size")
    letter = letters[@addr(uniform_discrete(1, length(letters)), "letter")]
    angle = @addr(uniform_continuous(-1, 1), "angle")
    fontsize = min_size + Int(floor((max_size - min_size + 1) * size))
    object = Object(height * x, width * y, angle * 45, fontsize, letter)

    # render
    image = render(object)

    # blur it
    blur_amount = 1
    blurred = imfilter(image, Kernel.gaussian(blur_amount))

    # add speckle
    mat = convert(Matrix{Float64}, blurred)
    noise = 0.1
    @addr(noisy_matrix(mat, noise), "image")
end

const num_input = width * height
const num_hidden1 = 100
const num_hidden2 = 100
const num_output = 11

relu(x) = x .* (x .> 0)

function hidden_layer(W, b, input)
    relu(W * input + b)
end
output_layer(W, b, input) = W * input + b

function inference_network(image, W1, b1, W2, b2, W3, b3)
    h1 = hidden_layer(W1, b1, image)
    h2 = hidden_layer(W2, b2, h1)
    output_layer(W3, b3, h2)
end

@gen function proposal(prev_choices)
    @param W1::Matrix{Float64}
    @param b1::Vector{Float64}
    @param W2::Matrix{Float64}
    @param b2::Vector{Float64}
    @param W3::Matrix{Float64}
    @param b3::Vector{Float64}

    image = prev_choices["image"][:]

    # inference network
    outputs = inference_network(image, W1, b1, W2, b2, W3, b3)
    x_mu = outputs[1]
    x_std = exp.(outputs[2])
    y_mu = outputs[3]
    y_std = exp.(outputs[4])
    r_mu = exp.(outputs[5])
    r_std = exp.(outputs[6])

    # distribution over sizes of the letters
    size_alpha = exp(outputs[7])
    size_beta = exp(outputs[8])

    # distribution over 3 possible letters
    log_letter_dist = outputs[9:9 + length(letters)-1]
    letter_dist = exp.(log_letter_dist)
    letter_dist = letter_dist / sum(letter_dist)

    @addr(normal(x_mu, x_std), "x")
    @addr(normal(y_mu, y_std), "y")
    @addr(normal(r_mu, r_std), "angle")
    @addr(beta(size_alpha, size_beta), "size")
    @addr(categorical(letter_dist), "letter")
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
export
    OdrManagerLite,
    Position,

    loadfile,
    printdata,
    create_position,
    activate_position,

    get_trackpos,
    get_lanepos,
    get_inertialpos,
    get_footpoint,

    set_pos,
    set_trackpos,
    set_lanepos,
    set_inertialpos,

    track2inertial,
    inertial2track,
    lane2inertial,
    inertial2lane,
    footpoint2inertial,

    print_odrmanagerlite,
    get_curvature,
    get_track_len,
    get_lane_width

type OdrManagerLite
    ptr::Ptr{Void}
    has_activated_position::Bool

    function OdrManagerLite()
        ptr = ccall( (:createOdrManagerLite, LIB_ODRMGR), Ptr{Void}, () )
        odrmanager = new(ptr, false)
        finalizer(odrmanager, obj -> begin
            ccall( (:free_OdrManagerLite, LIB_ODRMGR), Void, (Ptr{Void},), obj.ptr )
        end)
        odrmanager
    end
end

type Position
    ptr::Ptr{Void}
    Position(ptr::Ptr{Void}) = new(ptr)
end

loadfile(mgr::OdrManagerLite, name::String) =
    ccall((:odr_manager_loadFile, LIB_ODRMGR), Bool, (Ptr{Void}, Ptr{UInt8}), mgr.ptr, name)

printdata(mgr::OdrManagerLite) =
    ccall((:odr_manager_printData, LIB_ODRMGR), Void, (Ptr{Void},), mgr.ptr )

function create_position(mgr::OdrManagerLite)
    mgr.has_activated_position = true
    ptr = ccall((:odr_manager_createPosition, LIB_ODRMGR), Ptr{Void}, (Ptr{Void},), mgr.ptr)
    Position(ptr)
end

function activate_position(mgr::OdrManagerLite, pos::Position)
    mgr.has_activated_position = true
    ccall((:odr_manager_activatePosition, LIB_ODRMGR), Void, (Ptr{Void}, Ptr{Position}), mgr.ptr, pos.ptr)
end

function get_trackpos(mgr::OdrManagerLite)
    if !mgr.has_activated_position
        warn("OdrManagerLite does not have an activated position")
    else
        coord = TrackCoord()
        ccall((:odr_manager_getTrackPos, LIB_ODRMGR), Void, (Ptr{Void},Ptr{Void}), mgr.ptr, pointer_from_objref(coord))
        return coord
    end
end
function get_lanepos(mgr::OdrManagerLite)
    if !mgr.has_activated_position
        warn("OdrManagerLite does not have an activated position")
    else
        coord = LaneCoord()
        ccall((:odr_manager_getLanePos, LIB_ODRMGR), Void, (Ptr{Void},Ptr{Void}), mgr.ptr, pointer_from_objref(coord))
        return coord
    end
end
function get_inertialpos(mgr::OdrManagerLite)
    if !mgr.has_activated_position
        warn("OdrManagerLite does not have an activated position")
    else
        coord = Coord()
        ccall((:odr_manager_getInertialPos, LIB_ODRMGR), Void, (Ptr{Void},Ptr{Void}), mgr.ptr, pointer_from_objref(coord))
        return coord
    end
end
function get_footpoint(mgr::OdrManagerLite)
    if !mgr.has_activated_position
        warn("OdrManagerLite does not have an activated position")
    else
        coord = Coord()
        ccall((:odr_manager_getFootPoint, LIB_ODRMGR), Void, (Ptr{Void},Ptr{Void}), mgr.ptr, pointer_from_objref(coord))
        return coord
    end
end

set_pos(mgr::OdrManagerLite, value::TrackCoord) =
    ccall((:odr_manager_setpos_track_coord, LIB_ODRMGR), Void, (Ptr{Void}, Ptr{TrackCoord}), mgr.ptr, pointer_from_objref(value))
set_pos(mgr::OdrManagerLite, value::LaneCoord) =
    ccall((:odr_manager_setpos_lane_coord, LIB_ODRMGR), Void, (Ptr{Void}, Ptr{LaneCoord}), mgr.ptr, pointer_from_objref(value))
set_pos(mgr::OdrManagerLite, value::Coord) =
    ccall((:odr_manager_setpos_coord, LIB_ODRMGR), Void, (Ptr{Void}, Ptr{Coord}), mgr.ptr, pointer_from_objref(value))

set_trackpos(mgr::OdrManagerLite, value::TrackCoord) =
    ccall((:odr_manager_set_track_pos_track_coord, LIB_ODRMGR), Void, (Ptr{Void}, Ptr{TrackCoord}), mgr.ptr, pointer_from_objref(value))
set_trackpos(mgr::OdrManagerLite, id::Integer, s::Real, t::Real=0.0) =
    ccall((:odr_manager_set_track_pos_s_t, LIB_ODRMGR), Void, (Ptr{Void}, Cint, Cdouble, Cdouble), mgr.ptr, id, s, t)

set_lanepos(mgr::OdrManagerLite, value::LaneCoord) =
    ccall((:odr_manager_setLanePos_with_lanecoord, LIB_ODRMGR), Void, (Ptr{Void}, Ptr{LaneCoord}), mgr.ptr, pointer_from_objref(value))
set_lanepos(mgr::OdrManagerLite, trackId::Integer, laneId::Integer, s::Real, offset::Real=0.0) =
    ccall((:odr_manager_setLanePos, LIB_ODRMGR), Void, (Ptr{Void}, Cint, Cint, Cdouble, Cdouble), mgr.ptr, trackId, laneId, s, offset)

set_inertialpos(mgr::OdrManagerLite, x::Real, y::Real, z::Real) =
    ccall((:odr_manager_setInertialPos, LIB_ODRMGR), Void, (Ptr{Void}, Cdouble, Cdouble, Cdouble), mgr.ptr, x, y, z)

track2inertial(mgr::OdrManagerLite) =
     ccall((:odr_manager_track2inertial, LIB_ODRMGR), Bool, (Ptr{Void},), mgr.ptr)

inertial2track(mgr::OdrManagerLite) =
     ccall((:odr_manager_inertial2track, LIB_ODRMGR), Bool, (Ptr{Void},), mgr.ptr)

lane2inertial(mgr::OdrManagerLite) =
     ccall((:odr_manager_lane2inertial, LIB_ODRMGR), Bool, (Ptr{Void},), mgr.ptr)

inertial2lane(mgr::OdrManagerLite) =
     ccall((:odr_manager_inertial2lane, LIB_ODRMGR), Bool, (Ptr{Void},), mgr.ptr)

footpoint2inertial(mgr::OdrManagerLite) =
    ccall((:odr_manager_footPoint2inertial, LIB_ODRMGR), Void, (Ptr{Void},), mgr.ptr)

print_odrmanagerlite(mgr::OdrManagerLite, indent::Integer=0) =
    ccall((:odr_manager_print, LIB_ODRMGR), Void, (Ptr{Void},Cint), mgr.ptr, indent)

get_curvature(mgr::OdrManagerLite) =
     ccall((:odr_manager_getCurvature, LIB_ODRMGR), Cdouble, (Ptr{Void},), mgr.ptr)

get_track_len(mgr::OdrManagerLite, trackId::Integer) =
    ccall((:odr_manager_getTrackLen, LIB_ODRMGR), Cdouble, (Ptr{Void}, Cint), mgr.ptr, trackId)

get_lane_width(mgr::OdrManagerLite) =
    ccall((:odr_manager_getLaneWidth, LIB_ODRMGR), Cdouble, (Ptr{Void},), mgr.ptr)

max_stars_repo_name: sisl/OpenDriveManagerLite

max_issues_repo_name: sisl/OpenDriveManagerLite

max_forks_repo_name: sisl/OpenDriveManagerLite