text:
export RandomizedSystem, square_up

"""
    RandomizedSystem(F::Union{System,AbstractSystem}, k::Integer;
            identity_block = true,
            compile = $(COMPILE_DEFAULT[]))

Given a ``n × N`` system `F`` this constructs the system
``\\mathfrak{R}(F; k)(x) = A⋅F(x)`` where ``A`` is a ``k × N`` matrix.
If `identity_block = true` then the first `k` columns of ``A`` form an identity matrix.
The other entries are random complex numbers. See Chapter 13.5 in [^SW05] for more details.

    RandomizedSystem(F::Union{System,AbstractSystem}, A::Matrix{ComplexF64};
        identity_block = true,
        compile = $(COMPILE_DEFAULT[]))

Explicitly provide the used randomization matrix. If `identity = true` then `A` has to be `k × (n - k)` matrix.
Otherwise `A` has to be a `k × n` matrix..

[^SW05]: Sommese, A. J., & Wampler, C. W. (2005). The Numerical Solution of Systems of Polynomials Arising in Engineering and Science. World Scientific.
"""
struct RandomizedSystem{S<:AbstractSystem} <: AbstractSystem
    system::S
    A::Matrix{ComplexF64}
    identity_block::Bool
    u::Vector{ComplexF64}
    ū::Vector{ComplexDF64}
    U::Matrix{ComplexF64}
    taylor_u::TaylorVector{4,ComplexF64}
end

function RandomizedSystem(
    F::Union{AbstractSystem,System},
    k::Integer;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    identity_block = true,
)
    n, N = size(F)
    if identity_block
        A = randn(ComplexF64, k, n - k)
    else
        # Create a random matrix unit upper triangular
        A = [LA.UnitUpperTriangular(rand_unitary_matrix(k)) randn(ComplexF64, k, n - k)]
    end
    RandomizedSystem(F, A; identity_block = identity_block, compile = compile)
end

RandomizedSystem(
    F::System,
    A::Matrix{ComplexF64};
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    kwargs...,
) = RandomizedSystem(fixed(F; compile = compile), A; kwargs...)

function RandomizedSystem(
    F::AbstractSystem,
    A::Matrix{ComplexF64};
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    identity_block = true,
)
    n, N = size(F)
    is_compatible = identity_block ? size(A, 1) + size(A, 2) == n : size(A, 2) == n
    if !is_compatible
        throw(ArgumentError("Randomization matrix has incompatible size."))
    end

    u = zeros(ComplexF64, n)
    ū = zeros(ComplexDF64, n)
    U = zeros(ComplexF64, n, N)
    taylor_u = TaylorVector{4}(ComplexF64, n)
    RandomizedSystem(F, A, identity_block, u, ū, U, taylor_u)
end

function Base.show(io::IO, mime::MIME"text/plain", F::RandomizedSystem)
    println(io, typeof(F), ":")
    println(io, "A:")
    show(io, mime, F.A)
    println(io, "\n\nF:")
    show(io, mime, F.system)
end

Base.size(F::RandomizedSystem) = (size(F.A, 1), last(size(F.system)))
ModelKit.variables(F::RandomizedSystem) = variables(F.system)
ModelKit.parameters(F::RandomizedSystem) = parameters(F.system)
ModelKit.variable_groups(F::RandomizedSystem) = variable_groups(F.system)

"""
    square_up(F::Union{System, AbstractSystem}; identity_block = true, compile = compile = $(COMPILE_DEFAULT[]))

Creates the [`RandomizedSystem`](@ref) ``\\mathfrak{R}(F(x); N)`` where ``N`` is the number
of variables of `F`.
"""
square_up(F::Union{AbstractSystem,System}; kwargs...) =
    RandomizedSystem(F, last(size(F)); kwargs...)

@inline function randomize!(u, A, v::AbstractVector; identity_block::Bool)
    n, m = size(A, 1), length(v)

    if identity_block
        for i = 1:n
            u[i] = v[i]
        end
        for j = 1:(m-n)
            vnj = v[n+j]
            for i = 1:n
                u[i] += A[i, j] * vnj
            end
        end
    else
        for i = 1:n
            u[i] = 0
        end
        for j = 1:m, i = 1:n
            u[i] += A[i, j] * v[j]
        end
    end

    u
end

@inline function randomize!(
    u::Union{AbstractMatrix,TaylorVector},
    A,
    v::Union{AbstractMatrix,TaylorVector},
    ncols::Int = size(A, 1);
    identity_block::Bool,
)
    n, m = size(A, 1), size(v, 1)

    if identity_block
        for j = 1:ncols, i = 1:n
            u[i, j] = v[i, j]
        end
        for j = 1:ncols, k = 1:(m-n), i = 1:n
            u[i, j] += A[i, k] * v[n+k, j]
        end
    else
        for j = 1:ncols, i = 1:n
            u[i, j] = 0
        end
        for j = 1:ncols, k = 1:m, i = 1:n
            u[i, j] += A[i, k] * v[k, j]
        end
    end

    u
end

function (F::RandomizedSystem)(x, p = nothing)
    if F.identity_block
        [LA.I F.A] * F.system(x, p)
    else
        F.A * F.system(x, p)
    end
end

function ModelKit.evaluate!(u, F::RandomizedSystem, x::Vector{ComplexDF64}, p = nothing)
    evaluate!(F.ū, F.system, x, p)
    randomize!(u, F.A, F.ū; identity_block = F.identity_block)
    u
end
function ModelKit.evaluate!(u, F::RandomizedSystem, x, p = nothing)
    evaluate!(F.u, F.system, x, p)
    randomize!(u, F.A, F.u; identity_block = F.identity_block)
    u
end

function ModelKit.evaluate_and_jacobian!(u, U, F::RandomizedSystem, x, p = nothing)
    evaluate_and_jacobian!(F.u, F.U, F.system, x, p)
    randomize!(u, F.A, F.u; identity_block = F.identity_block)
    randomize!(U, F.A, F.U; identity_block = F.identity_block)
end

function ModelKit.taylor!(
    u::AbstractVector,
    v::Val{N},
    F::RandomizedSystem,
    tx,
    p = nothing,
) where {N}
    taylor!(F.u, v, F.system, tx, p)
    randomize!(u, F.A, F.u; identity_block = F.identity_block)
end
function ModelKit.taylor!(
    u::TaylorVector,
    v::Val{N},
    F::RandomizedSystem,
    tx,
    p = nothing,
) where {N}
    taylor!(F.taylor_u, v, F.system, tx, p)
    randomize!(u, F.A, F.taylor_u, N + 1; identity_block = F.identity_block)
end


text:
using Base.Iterators: product
using Dates
using Documenter: doctest
using Intervals
using Test
using TimeZones

const BOUND_PERMUTATIONS = product((Closed, Open), (Closed, Open))

@testset "Intervals" begin
    include("inclusivity.jl")
    include("endpoint.jl")
    include("interval.jl")
    include("anchoredinterval.jl")
    include("comparisons.jl")
    include("plotting.jl")

    doctest(Intervals)
end

ext: jl


text:
using DSGE, ModelConstructors, HDF5, Test, Dates
path = dirname(@__FILE__)

### Model
model = Model1002("ss10")

### Parameters

@testset "Test parameter bounds checking" begin
    for θ in model.parameters
        @test isa(θ, AbstractParameter)
        if !θ.fixed
            (left, right) = θ.valuebounds
            @test left < θ.value < right
        end
    end
end

### Model indices

@testset "Checking model indices" begin
    # Endogenous states
    endo = model.endogenous_states
    @test length(endo) == 68
    @test endo[:Ez_t] == 60

    # Exogenous shocks
    exo = model.exogenous_shocks
    @test length(exo) == 24
    @test exo[:corepce_sh] == 16

    # Expectation shocks
    ex = model.expected_shocks
    @test length(ex) == 13
    @test ex[:EL_f_sh] == 12

    # Equations
    eq = model.equilibrium_conditions
    @test length(eq) == 68
    @test eq[:eq_Ez] == 59

    # Additional states
    endo_new = model.endogenous_states_augmented
    @test length(endo_new) == 16
    @test endo_new[:y_t1] == 69

    # Observables
    obs = model.observables
    @test length(obs) == 19
    @test obs[:obs_tfp] == 12
end

### Equilibrium conditions
Γ0, Γ1, C, Ψ, Π = eqcond(model)

# Transition and measurement equations
TTT, RRR, CCC = solve(model)
meas = measurement(model, TTT, RRR, CCC)

### Pseudo-measurement equation
pseudo_meas = pseudo_measurement(model, TTT, RRR, CCC)

@testset "Check eqcond and state-space dimensions" begin
    # Matrices are of expected dimensions
    @test size(Γ0) == (68, 68)
    @test size(Γ1) == (68, 68)
    @test size(C) == (68,)
    @test size(Ψ) == (68, 24)
    @test size(Π) == (68, 13)

    @test size(meas[:ZZ]) == (19,84)
    @test size(meas[:DD]) == (19,)
    @test size(meas[:QQ]) == (24,24)
    @test size(meas[:EE]) == (19,19)

    @test size(TTT) == (84,84)
    @test size(RRR) == (84,24)
    @test size(CCC) == (84,)
end

G0_exp, G1_exp, C_exp, PSI_exp, PIE_exp = h5open("$path/eqcond.h5") do file
    read(file, "G0"), read(file, "G1"), read(file, "C"),
    read(file, "PSI"), read(file, "PIE")
end
TTT_exp, RRR_exp, CCC_exp = h5open("$path/transition.h5") do file
    read(file, "TTT"), read(file, "RRR"), read(file, "CCC")
end
Z_exp, D_exp, Q_exp, E_exp = h5open("$path/measurement.h5") do file
    read(file, "ZZ"), read(file, "DD"), read(file, "QQ"), read(file, "EE")
end
Z_pseudo_exp, D_pseudo_exp = h5open("$path/pseudo_measurement.h5") do file
    read(file, "ZZ_pseudo"), read(file, "DD_pseudo")
end

@testset "Compare eqcond, transition, measurement, and pseudo-measurement matrices against reference matrices" begin
    @test @test_matrix_approx_eq Γ0 G0_exp
    @test @test_matrix_approx_eq Γ1 G1_exp
    @test @test_matrix_approx_eq C C_exp
    @test @test_matrix_approx_eq Ψ PSI_exp
    @test @test_matrix_approx_eq Π PIE_exp

    @test @test_matrix_approx_eq TTT TTT_exp
    @test @test_matrix_approx_eq RRR RRR_exp
    @test @test_matrix_approx_eq CCC CCC_exp

    @test @test_matrix_approx_eq meas[:ZZ] Z_exp
    @test @test_matrix_approx_eq meas[:DD] D_exp
    @test @test_matrix_approx_eq meas[:QQ] Q_exp
    @test @test_matrix_approx_eq meas[:EE] E_exp

    @test @test_matrix_approx_eq pseudo_meas[:ZZ_pseudo] Z_pseudo_exp[1:21, :] # Test files had been generated w/extra pseudo-obs
    @test @test_matrix_approx_eq pseudo_meas[:DD_pseudo] D_pseudo_exp[1:21, :] # Test files had been generated w/extra pseudo-obs
end

@testset "Check ss62 is properly implemented" begin
    m = Model1002("ss62")
    for (i, d) in enumerate([date_presample_start(m), Date(2020, 3, 31), Date(2020, 6, 30), Date(2020, 9, 30),
                             Date(2020, 12, 31), Date(2021, 3, 31)])
        @test get_setting(m, :regime_dates)[i] == d
    end
    sys = compute_system(m)
    exo = m.exogenous_shocks
    inds1 = 1:24 # standard shock indices
    inds2 = vcat(1:14, 17:24) # standard shock indices excluding measurement error
    inds3 = 25:30 # covid shock indices
    for k in [:σ_biidc, :σ_φ, :σ_ziid]
        @test haskey(m[k].regimes, :prior)
        @test get(regime_prior(m[k], 1)).ν == get(regime_prior(m[k], 2)).ν
        @test get(regime_prior(m[k], 1)).ν == get(regime_prior(m[k], 3)).ν
        @test get(regime_prior(m[k], 1)).τ == get(regime_prior(m[k], 2)).τ
        @test get(regime_prior(m[k], 1)).τ != get(regime_prior(m[k], 3)).τ
    end

    for i in 2:get_setting(m, :n_regimes)
        @test sys[1, :TTT] ≈ sys[i, :TTT]
        @test sys[1, :CCC] ≈ sys[i, :CCC]
        @test sys[1, :ZZ]  ≈ sys[i, :ZZ]
        @test sys[1, :DD]  ≈ sys[i, :DD]

        if i >= 5
            @test sys[1, :QQ][inds1, inds1] ≈ sys[i, :QQ][inds1, inds1]
        end
        if i >= 4
            @test sys[1, :QQ][inds2, inds2] ≈ sys[i, :QQ][inds2, inds2]
        end
        if 2 <= i <= 4
            @test sys[i, :QQ][exo[:biidc_sh], exo[:biidc_sh]] == 16
            @test sys[i, :QQ][exo[:ziid_sh], exo[:ziid_sh]] == 25.
            @test sys[i, :QQ][exo[:φ_sh], exo[:φ_sh]] == 160000.
            @test (sys[1, :QQ][exo[:corepce_sh], exo[:corepce_sh]] * 10^2) == sys[i, :QQ][exo[:corepce_sh], exo[:corepce_sh]]
            @test (sys[1, :QQ][exo[:gdpdef_sh], exo[:gdpdef_sh]] * 10^2) == sys[i, :QQ][exo[:gdpdef_sh], exo[:gdpdef_sh]]
        elseif i == 5
            @test sys[i, :QQ][exo[:biidc_sh], exo[:biidc_sh]] == .04^2
            @test sys[i, :QQ][exo[:ziid_sh], exo[:ziid_sh]] == .05^2
            @test sys[i, :QQ][exo[:φ_sh], exo[:φ_sh]] == 16.
        end

        if 2 <= i <= 3
            @test (sys[1, :QQ][exo[:rm_sh], exo[:rm_sh]] * 10^2) == sys[i, :QQ][exo[:rm_sh], exo[:rm_sh]]
            @test (sys[1, :QQ][exo[:b_sh], exo[:b_sh]] * .25^2) == sys[i, :QQ][exo[:b_sh], exo[:b_sh]]
            @test (sys[1, :QQ][exo[:g_sh], exo[:g_sh]] * .25^2) == sys[i, :QQ][exo[:g_sh], exo[:g_sh]]
        end

        @test sys[i, :QQ][exo[:ziid_shl1], exo[:ziid_shl1]] == 0.
        @test sys[i, :QQ][exo[:φ_shl1], exo[:φ_shl1]] == 0.
        if 2 <= i <= 3 || i >= 6
            @test sys[i, :QQ][exo[:biidc_shl1], exo[:biidc_shl1]] == 0.
        else
            @test sys[i, :QQ][exo[:biidc_shl1], exo[:biidc_shl1]] == 16.
        end
    end
end

@testset "Check ss63 is properly implemented" begin
    m = Model1002("ss63")
    for (i, d) in enumerate([date_presample_start(m), Date(2020, 3, 31), Date(2020, 6, 30), Date(2020, 9, 30),
                             Date(2020, 12, 31), Date(2021, 3, 31)])
        @test get_setting(m, :regime_dates)[i] == d
    end
    sys = compute_system(m)
    exo = m.exogenous_shocks
    inds1 = 1:24 # standard shock indices
    inds2 = vcat(1:14, 17:24) # standard shock indices excluding measurement error
    inds3 = 25:30 # covid shock indices
    for k in [:σ_biidc, :σ_φ, :σ_ziid]
        @test haskey(m[k].regimes, :prior)
        @test get(regime_prior(m[k], 1)).ν == get(regime_prior(m[k], 2)).ν
        @test get(regime_prior(m[k], 1)).ν == get(regime_prior(m[k], 3)).ν
        @test get(regime_prior(m[k], 1)).τ == get(regime_prior(m[k], 2)).τ
        @test get(regime_prior(m[k], 1)).τ != get(regime_prior(m[k], 3)).τ
    end

    for i in 2:get_setting(m, :n_regimes)
        @test sys[1, :TTT] ≈ sys[i, :TTT]
        @test sys[1, :CCC] ≈ sys[i, :CCC]
        @test sys[1, :ZZ]  ≈ sys[i, :ZZ]
        @test sys[1, :DD]  ≈ sys[i, :DD]

        if i >= 5
            @test sys[1, :QQ][inds1, inds1] ≈ sys[i, :QQ][inds1, inds1]
        end
        if i >= 4
            @test sys[1, :QQ][inds2, inds2] ≈ sys[i, :QQ][inds2, inds2]
        end
        if 2 <= i <= 4
            @test sys[i, :QQ][exo[:biidc_sh], exo[:biidc_sh]] == 16
            @test sys[i, :QQ][exo[:ziid_sh], exo[:ziid_sh]] == 25.
            @test sys[i, :QQ][exo[:φ_sh], exo[:φ_sh]] == 160000.
            @test sys[1, :QQ][exo[:corepce_sh], exo[:corepce_sh]] * 10^2 ≈ sys[i, :QQ][exo[:corepce_sh], exo[:corepce_sh]]
            @test sys[1, :QQ][exo[:gdpdef_sh], exo[:gdpdef_sh]] * 10^2 ≈ sys[i, :QQ][exo[:gdpdef_sh], exo[:gdpdef_sh]]
        elseif i == 5
            @test sys[i, :QQ][exo[:biidc_sh], exo[:biidc_sh]] == .04^2
            @test sys[i, :QQ][exo[:ziid_sh], exo[:ziid_sh]] == .05^2
            @test sys[i, :QQ][exo[:φ_sh], exo[:φ_sh]] == 16.
        end

        if 2 <= i <= 3
            @test sys[1, :QQ][exo[:rm_sh], exo[:rm_sh]] * 10^2 == sys[i, :QQ][exo[:rm_sh], exo[:rm_sh]]
            @test sys[1, :QQ][exo[:g_sh], exo[:g_sh]] == sys[i, :QQ][exo[:g_sh], exo[:g_sh]]
            @test sys[1, :QQ][exo[:b_sh], exo[:b_sh]] == sys[i, :QQ][exo[:b_sh], exo[:b_sh]]
        end

        @test sys[i, :QQ][exo[:ziid_shl1], exo[:ziid_shl1]] == 0.
        @test sys[i, :QQ][exo[:φ_shl1], exo[:φ_shl1]] == 0.
        if 2 <= i <= 3 || i >= 6
            @test sys[i, :QQ][exo[:biidc_shl1], exo[:biidc_shl1]] == 0.
        else
            @test sys[i, :QQ][exo[:biidc_shl1], exo[:biidc_shl1]] == 16.
        end
    end
end


text:
using ChainRulesCore, Test
using LinearAlgebra, SparseArrays
using OffsetArrays, BenchmarkTools

# Like ForwardDiff.jl's Dual
struct Dual{T<:Real} <: Real
    value::T
    partial::T
end
Base.real(x::Dual) = x
Base.float(x::Dual) = Dual(float(x.value), float(x.partial))
Base.zero(x::Dual) = Dual(zero(x.value), zero(x.partial))

# Trivial struct
struct NoSuperType end

@testset "projection" begin

    #####
    ##### `Base`
    #####

    @testset "Base: numbers" begin
        # real / complex
        @test ProjectTo(1.0)(2.0 + 3im) === 2.0
        @test ProjectTo(1.0 + 2.0im)(3.0) === 3.0 + 0.0im
        @test ProjectTo(2.0+3.0im)(1+1im) === 1.0+1.0im
        @test ProjectTo(2.0)(1+1im) === 1.0

        # storage
        @test ProjectTo(1)(pi) === pi
        @test ProjectTo(1 + im)(pi) === ComplexF64(pi)
        @test ProjectTo(1//2)(3//4) === 3//4
        @test ProjectTo(1.0f0)(1 / 2) === 0.5f0
        @test ProjectTo(1.0f0 + 2im)(3) === 3.0f0 + 0im
        @test ProjectTo(big(1.0))(2) === 2
        @test ProjectTo(1.0)(2) === 2.0

        # Tangents
        ProjectTo(1.0f0 + 2im)(Tangent{ComplexF64}(re=1, im=NoTangent())) === 1.0f0 + 0.0f0im
    end

    @testset "Dual" begin # some weird Real subtype that we should basically leave alone
        @test ProjectTo(1.0)(Dual(1.0, 2.0)) isa Dual
        @test ProjectTo(1.0)(Dual(1, 2)) isa Dual

        # real & complex
        @test ProjectTo(1.0 + 1im)(Dual(1.0, 2.0)) isa Complex{<:Dual}
        @test ProjectTo(1.0 + 1im)(
            Complex(Dual(1.0, 2.0), Dual(1.0, 2.0))
        ) isa Complex{<:Dual}
        @test ProjectTo(1.0)(Complex(Dual(1.0, 2.0), Dual(1.0, 2.0))) isa Dual

        # Tangent
        @test ProjectTo(Dual(1.0, 2.0))(Tangent{Dual}(; value=1.0)) isa Tangent
    end

    @testset "Base: arrays of numbers" begin
        pvec3 = ProjectTo([1, 2, 3])
        @test pvec3(1.0:3.0) === 1.0:3.0
        @test pvec3(1:3) == 1.0:3.0  # would prefer ===, map(Float64, dx) would do that, not important
        @test pvec3([1, 2, 3 + 4im]) == 1:3
        @test eltype(pvec3([1, 2, 3.0f0])) === Float64

        # reshape
        @test pvec3(reshape([1, 2, 3], 3, 1)) isa Vector
        @test_throws DimensionMismatch pvec3(reshape([1, 2, 3], 1, 3))
        @test_throws DimensionMismatch pvec3([1, 2, 3, 4])

        pmat = ProjectTo(rand(2, 2) .+ im)
        @test pmat([1 2; 3 4.0+5im]') isa Adjoint     # pass-through
        @test pmat([1 2; 3 4]') isa Matrix{ComplexF64}  # broadcast type change

        pmat2 = ProjectTo(rand(2, 2)')
        @test pmat2([1 2; 3 4.0+5im]) isa Matrix  # adjoint matrices are not re-created

        prow = ProjectTo([1im 2 3im])
        @test prow(transpose([1, 2, 3 + 4.0im])) == [1 2 3 + 4im]
        @test prow(transpose([1, 2, 3 + 4.0im])) isa Matrix  # row vectors may not pass through
        @test prow(adjoint([1, 2, 3 + 5im])) == [1 2 3 - 5im]
        @test prow(adjoint([1, 2, 3])) isa Matrix

        # some bugs
        @test pvec3(fill(NoTangent(), 3)) === NoTangent()  #410, was an array of such
        @test ProjectTo([pi])([1]) isa Vector{Int}  #423, was Irrational -> Bool -> NoTangent
    end

    @testset "Base: arrays of arrays, etc" begin
        pvecvec = ProjectTo([[1, 2], [3, 4, 5]])
        @test pvecvec([1:2, 3:5])[1] == 1:2
        @test pvecvec([[1, 2 + 3im], [4 + 5im, 6, 7]])[2] == [4, 6, 7]
        @test pvecvec(hcat([1:2, hcat(3:5)]))[2] isa Vector  # reshape inner & outer

        pvecvec2 = ProjectTo(reshape(Any[[1 2], [3 4 5]], 1, 2))  # a row of rows
        y1 = pvecvec2([[1, 2], [3, 4, 5]]')
        @test y1[1] == [1 2]
        @test !(y1 isa Adjoint) && !(y1[1] isa Adjoint)

        # arrays of other things
        @test ProjectTo([:x, :y]) isa ProjectTo{NoTangent}
        @test ProjectTo(Any['x', "y"]) isa ProjectTo{NoTangent}
        @test ProjectTo([(1,2), (3,4), (5,6)]) isa ProjectTo{AbstractArray}

        @test ProjectTo(Any[1, 2])(1:2) == [1.0, 2.0]  # projects each number.
        @test Tuple(ProjectTo(Any[1, 2 + 3im])(1:2)) === (1.0, 2.0 + 0.0im)
        @test ProjectTo(Any[true, false]) isa ProjectTo{NoTangent}

        # empty arrays
        @test isempty(ProjectTo([])(1:0))
        @test_throws DimensionMismatch ProjectTo(Int[])([2])
        @test ProjectTo(Bool[]) isa ProjectTo{NoTangent}
    end

    @testset "Base: zero-arrays" begin
        pzed = ProjectTo(fill(1.0))
        @test pzed(fill(3.14)) == fill(3.14)  # easy
        @test pzed(fill(3)) == fill(3.0)      # broadcast type change must not produce number
        @test pzed(hcat(3.14)) == fill(3.14)  # reshape of length-1 array
        @test pzed(3.14 + im) == fill(3.14)   # re-wrap of a scalar number

        @test_throws DimensionMismatch ProjectTo([1])(3.14 + im) # other array projectors don't accept numbers
        @test_throws DimensionMismatch ProjectTo(hcat([1, 2]))(3.14)
        @test pzed isa ProjectTo{AbstractArray}
    end

    @testset "Base: Ref" begin
        pref = ProjectTo(Ref(2.0))
        @test pref(Ref(3 + im)).x === 3.0
        @test pref(Tangent{Base.RefValue}(x = 3 + im)).x === 3.0
        @test pref(4).x === 4.0  # also re-wraps scalars
        @test pref(Ref{Any}(5.0)) isa Tangent{<:Base.RefValue}

        pref2 = ProjectTo(Ref{Any}(6 + 7im))
        @test pref2(Ref(8)).x === 8.0 + 0.0im
        @test pref2(Tangent{Base.RefValue}(x = 8)).x === 8.0 + 0.0im

        prefvec = ProjectTo(Ref([1, 2, 3 + 4im]))  # recurses into contents
        @test prefvec(Ref(1:3)).x isa Vector{ComplexF64}
        @test prefvec(Tangent{Base.RefValue}(x = 1:3)).x isa Vector{ComplexF64}
        @test_skip @test_throws DimensionMismatch prefvec(Tangent{Base.RefValue}(x = 1:5))

        @test ProjectTo(Ref(true)) isa ProjectTo{NoTangent}
        @test ProjectTo(Ref([false]')) isa ProjectTo{NoTangent}
    end

    @testset "Base: non-diff" begin
        @test ProjectTo(:a)(1) == NoTangent()
        @test ProjectTo('b')(2) == NoTangent()
        @test ProjectTo("cde")(345) == NoTangent()
    end

    #####
    ##### `LinearAlgebra`
    #####

    @testset "LinearAlgebra: $adj vectors" for adj in [transpose, adjoint]
        # adjoint vectors
        padj = ProjectTo(adj([1, 2, 3]))
        adjT = typeof(adj([1, 2, 3.0]))
        @test padj(transpose(1:3)) isa adjT
        @test padj([4 5 6 + 7im]) isa adjT
        @test padj([4.0 5.0 6.0]) isa adjT

        @test_throws DimensionMismatch padj([1, 2, 3])
        @test_throws DimensionMismatch padj([1 2 3]')
        @test_throws DimensionMismatch padj([1 2 3 4])

        padj_complex = ProjectTo(adj([1, 2, 3 + 4im]))
        @test padj_complex([4 5 6 + 7im]) == [4 5 6 + 7im]
        @test padj_complex(transpose([4, 5, 6 + 7im])) == [4 5 6 + 7im]
        @test padj_complex(adjoint([4, 5, 6 + 7im])) == [4 5 6 - 7im]

        # evil test case
        if VERSION >= v"1.7-"  # up to 1.6  Vector[[1,2,3]]'  is an error, not sure why it's called
            xs = adj(Any[Any[1, 2, 3], Any[4 + im, 5 - im, 6 + im, 7 - im]])
            pvecvec3 = ProjectTo(xs)
            @test pvecvec3(xs)[1] == [1 2 3]
            @test pvecvec3(xs)[2] == adj.([4 + im 5 - im 6 + im 7 - im])
            @test pvecvec3(xs)[2] isa LinearAlgebra.AdjOrTransAbsMat{ComplexF64,<:Vector}
            @test pvecvec3(collect(xs))[1] == [1 2 3]
            ys = permutedims([[1 2 3 + im], Any[4 5 6 7 + 8im]])
            @test pvecvec3(ys)[1] == [1 2 3]
            @test pvecvec3(ys)[2] == [4 5 6 7 + 8im]
            @test pvecvec3(xs)[2] isa LinearAlgebra.AdjOrTransAbsMat{ComplexF64,<:Vector}
            @test pvecvec3(ys) isa LinearAlgebra.AdjOrTransAbsVec

            zs = adj([[1 2; 3 4], [5 6; 7 8+im]'])
            pvecmat = ProjectTo(zs)
            @test pvecmat(zs) == zs
            @test pvecmat(collect.(zs)) == zs
            @test pvecmat(collect.(zs)) isa LinearAlgebra.AdjOrTransAbsVec
        end

        # issue #410
        @test padj([NoTangent() NoTangent() NoTangent()]) === NoTangent()

        @test ProjectTo(adj([true, false]))([1 2]) isa AbstractZero
        @test ProjectTo(adj([[true], [false]])) isa ProjectTo{<:AbstractZero}
    end

    @testset "LinearAlgebra: dense structured matrices" begin
        psymm = ProjectTo(Symmetric(rand(3, 3)))
        @test psymm(reshape(1:9, 3, 3)) == [1.0 3.0 5.0; 3.0 5.0 7.0; 5.0 7.0 9.0]
        @test psymm(psymm(reshape(1:9, 3, 3))) == psymm(reshape(1:9, 3, 3))
        @test psymm(rand(ComplexF32, 3, 3, 1)) isa Symmetric{Float64}
        @test ProjectTo(Symmetric(randn(3, 3) .> 0))(randn(3, 3)) == NoTangent() # Bool

        pherm = ProjectTo(Hermitian(rand(3, 3) .+ im, :L))
        # NB, projection onto Hermitian subspace, not application of Hermitian constructor
        @test pherm(reshape(1:9, 3, 3) .+ im) == [1.0 3.0 5.0; 3.0 5.0 7.0; 5.0 7.0 9.0]
        @test pherm(pherm(reshape(1:9, 3, 3))) == pherm(reshape(1:9, 3, 3))
        @test pherm(rand(ComplexF32, 3, 3, 1)) isa Hermitian{ComplexF64}

        pupp = ProjectTo(UpperTriangular(rand(3, 3)))
        @test pupp(reshape(1:9, 3, 3)) == [1.0 4.0 7.0; 0.0 5.0 8.0; 0.0 0.0 9.0]
        @test pupp(pupp(reshape(1:9, 3, 3))) == pupp(reshape(1:9, 3, 3))
        @test pupp(rand(ComplexF32, 3, 3, 1)) isa UpperTriangular{Float64}
        @test ProjectTo(UpperTriangular(randn(3, 3) .> 0))(randn(3, 3)) == NoTangent()

        # an experiment with allowing subspaces which aren't subtypes
        @test psymm(Diagonal([1, 2, 3])) isa Diagonal{Float64}
        @test pupp(Diagonal([1, 2, 3 + 4im])) isa Diagonal{Float64}
    end

    @testset "LinearAlgebra: sparse structured matrices" begin
        pdiag = ProjectTo(Diagonal(1:3))
        @test pdiag(reshape(1:9, 3, 3)) == Diagonal([1, 5, 9])
        @test pdiag(pdiag(reshape(1:9, 3, 3))) == pdiag(reshape(1:9, 3, 3))
        @test pdiag(rand(ComplexF32, 3, 3)) isa Diagonal{Float64}
        @test pdiag(Diagonal(1.0:3.0)) === Diagonal(1.0:3.0)
        @test ProjectTo(Diagonal(randn(3) .> 0))(randn(3, 3)) == NoTangent()
        @test ProjectTo(Diagonal(randn(3) .> 0))(Diagonal(rand(3))) == NoTangent()

        pbi = ProjectTo(Bidiagonal(rand(3, 3), :L))
        @test pbi(reshape(1:9, 3, 3)) == [1.0 0.0 0.0; 2.0 5.0 0.0; 0.0 6.0 9.0]
        @test pbi(pbi(reshape(1:9, 3, 3))) == pbi(reshape(1:9, 3, 3))
        @test pbi(rand(ComplexF32, 3, 3)) isa Bidiagonal{Float64}
        bi = Bidiagonal(rand(3, 3) .+ im, :L)
        @test pbi(bi) == real(bi)  # reconstruct via generic_projector
        bu = Bidiagonal(rand(3, 3) .+ im, :U)  # differs but uplo, not type
        @test pbi(bu) == diagm(0 => diag(real(bu)))
        @test_throws DimensionMismatch pbi(rand(ComplexF32, 3, 2))

        pstri = ProjectTo(SymTridiagonal(Symmetric(rand(3, 3))))
        @test pstri(reshape(1:9, 3, 3)) == [1.0 3.0 0.0; 3.0 5.0 7.0; 0.0 7.0 9.0]
        @test pstri(pstri(reshape(1:9, 3, 3))) == pstri(reshape(1:9, 3, 3))
        @test pstri(rand(ComplexF32, 3, 3)) isa SymTridiagonal{Float64}
        stri = SymTridiagonal(Symmetric(rand(3, 3) .+ im))
        @test pstri(stri) == real(stri)
        @test_throws DimensionMismatch pstri(rand(ComplexF32, 3, 2))

        ptri = ProjectTo(Tridiagonal(rand(3, 3)))
        @test ptri(reshape(1:9, 3, 3)) == [1.0 4.0 0.0; 2.0 5.0 8.0; 0.0 6.0 9.0]
        @test ptri(ptri(reshape(1:9, 3, 3))) == ptri(reshape(1:9, 3, 3))
        @test ptri(rand(ComplexF32, 3, 3)) isa Tridiagonal{Float64}
        @test_throws DimensionMismatch ptri(rand(ComplexF32, 3, 2))
    end

    #####
    ##### `SparseArrays`
    #####

    @testset "SparseArrays" begin
        # vector
        v = sprand(30, 0.3)
        pv = ProjectTo(v)

        @test pv(v) == v
        @test pv(v .* (1 + im)) ≈ v  # same nonzero elements

        o = pv(ones(Int, 30, 1))  # dense array
        @test nnz(o) == nnz(v)

        v2 = sprand(30, 0.7)  # different nonzero elements
        @test pv(v2) == pv(collect(v2))

        # matrix
        m = sprand(10, 10, 0.3)
        pm = ProjectTo(m)

        @test pm(m) == m
        @test pm(m .* (1 + im)) ≈ m
        om = pm(ones(Int, 10, 10))
        @test nnz(om) == nnz(m)

        m2 = sprand(10, 10, 0.5)
        @test pm(m2) == pm(collect(m2))

        @test_throws DimensionMismatch pv(ones(Int, 1, 30))
        @test_throws DimensionMismatch pm(ones(Int, 5, 20))
    end

    #####
    ##### `OffsetArrays`
    #####

    @testset "OffsetArrays" begin
        # While there is no code for this, the rule that it checks axes(x) == axes(dx) else
        # reshape means that it restores offsets. (It throws an error on nontrivial size mismatch.)

        poffv = ProjectTo(OffsetArray(rand(3), 0:2))
        @test axes(poffv([1, 2, 3])) == (0:2,)
        @test axes(poffv(hcat([1, 2, 3]))) == (0:2,)

        @test axes(poffv(OffsetArray(rand(3), 0:2))) == (0:2,)
        @test axes(poffv(OffsetArray(rand(3, 1), 0:2, 0:0))) == (0:2,)
    end

    #####
    ##### `ChainRulesCore`
    #####

    @testset "pass-through" begin
        @test ProjectTo(NoSuperType()) === identity
    end

    @testset "AbstractZero" begin
        pz = ProjectTo(ZeroTangent())
        pz(0) == NoTangent()
        @test pz(ZeroTangent()) === ZeroTangent()  # not sure how NB this is to preserve
        @test pz(NoTangent()) === NoTangent()

        pb = ProjectTo(true) # Bool is categorical
        @test pb(2) === NoTangent()
        @test pb(ZeroTangent()) isa AbstractZero  # was a method ambiguity!

        # all projectors preserve Zero, and specific type, via one fallback method:
        @test ProjectTo(pi)(ZeroTangent()) === ZeroTangent()
        @test ProjectTo(pi)(NoTangent()) === NoTangent()
        pv = ProjectTo(sprand(30, 0.3))
        @test pv(ZeroTangent()) === ZeroTangent()
        @test pv(NoTangent()) === NoTangent()
    end

    @testset "Thunk" begin
        th = @thunk 1 + 2 + 3
        pth = ProjectTo(4 + 5im)(th)
        @test pth isa Thunk
        @test unthunk(pth) === 6.0 + 0.0im
    end

    @testset "Tangent" begin
        x = 1:3.0
        dx = Tangent{typeof(x)}(; step=0.1, ref=NoTangent());
        @test ProjectTo(x)(dx) isa Tangent
        @test ProjectTo(x)(dx).step === 0.1
        @test ProjectTo(x)(dx).offset isa AbstractZero

        pref = ProjectTo(Ref(2.0))
        dy = Tangent{typeof(Ref(2.0))}(x = 3+4im)
        @test pref(dy) isa Tangent{<:Base.RefValue}
        @test pref(dy).x === 3.0
    end

    @testset "display" begin
        @test repr(ProjectTo(1.1)) == "ProjectTo{Float64}()"
        @test occursin("ProjectTo{AbstractArray}(element", repr(ProjectTo([1, 2, 3])))
        str = repr(ProjectTo([1, 2, 3]'))
        @test eval(Meta.parse(str))(ones(1, 3)) isa Adjoint{Float64,Vector{Float64}}
    end

    VERSION > v"1.1" && @testset "allocation tests" begin
        # For sure these fail on Julia 1.0, not sure about 1.3 etc.
        # We only really care about current stable anyway
        # Each "@test 33 > ..." is zero on nightly, 32 on 1.5.

        pvec = ProjectTo(rand(10^3))
        @test 0 == @ballocated $pvec(dx) setup=(dx = rand(10^3))    # pass through
        @test 90 > @ballocated $pvec(dx) setup=(dx = rand(10^3, 1))  # reshape

        @test 33 > @ballocated ProjectTo(x)(dx) setup = (x = rand(10^3); dx = rand(10^3)) # including construction

        padj = ProjectTo(adjoint(rand(10^3)))
        @test 0 == @ballocated $padj(dx) setup=(dx = adjoint(rand(10^3)))
        @test 0 == @ballocated $padj(dx) setup=(dx = transpose(rand(10^3)))

        @test 33 > @ballocated ProjectTo(x')(dx') setup = (x = rand(10^3); dx = rand(10^3))

        pdiag = ProjectTo(Diagonal(rand(10^3)))
        @test 0 == @ballocated $pdiag(dx) setup=(dx = Diagonal(rand(10^3)))

        psymm = ProjectTo(Symmetric(rand(10^3, 10^3)))
        @test_broken 0 == @ballocated $psymm(dx) setup=(dx = Symmetric(rand(10^3, 10^3)))  # 64
    end
end


text:
CnMatrix = AbstractMatrix{Int};
FdMatrix = AbstractMatrix{Float64};
CnVector = AbstractVector{Int};
FdVector = AbstractVector{Float64};

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
function DateComponents(arg0::DateComponents, arg1::jint)
    return DateComponents((DateComponents, jint), arg0, arg1)
end

function DateComponents(arg0::jint)
    return DateComponents((jint,), arg0)
end

function DateComponents(arg0::jint, arg1::Month, arg2::jint)
    return DateComponents((jint, Month, jint), arg0, arg1, arg2)
end

function DateComponents(arg0::jint, arg1::jint)
    return DateComponents((jint, jint), arg0, arg1)
end

function DateComponents(arg0::jint, arg1::jint, arg2::jint)
    return DateComponents((jint, jint, jint), arg0, arg1, arg2)
end

function compare_to(obj::DateComponents, arg0::DateComponents)
    return jcall(obj, "compareTo", jint, (DateComponents,), arg0)
end

function create_from_week_components(::Type{DateComponents}, arg0::jint, arg1::jint, arg2::jint)
    return jcall(DateComponents, "createFromWeekComponents", DateComponents, (jint, jint, jint), arg0, arg1, arg2)
end

function equals(obj::DateComponents, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_calendar_week(obj::DateComponents)
    return jcall(obj, "getCalendarWeek", jint, ())
end

function get_day(obj::DateComponents)
    return jcall(obj, "getDay", jint, ())
end

function get_day_of_week(obj::DateComponents)
    return jcall(obj, "getDayOfWeek", jint, ())
end

function get_day_of_year(obj::DateComponents)
    return jcall(obj, "getDayOfYear", jint, ())
end

function get_j2000_day(obj::DateComponents)
    return jcall(obj, "getJ2000Day", jint, ())
end

function get_mjd(obj::DateComponents)
    return jcall(obj, "getMJD", jint, ())
end

function get_month(obj::DateComponents)
    return jcall(obj, "getMonth", jint, ())
end

function get_month_enum(obj::DateComponents)
    return jcall(obj, "getMonthEnum", Month, ())
end

function get_year(obj::DateComponents)
    return jcall(obj, "getYear", jint, ())
end

function hash_code(obj::DateComponents)
    return jcall(obj, "hashCode", jint, ())
end

function parse_date(::Type{DateComponents}, arg0::JString)
    return jcall(DateComponents, "parseDate", DateComponents, (JString,), arg0)
end

function to_string(obj::DateComponents)
    return jcall(obj, "toString", JString, ())
end


text:
__precompile__()

module PolicyGradient

using mjWrap
using LearningStrategies
import LearningStrategies: setup!, update!, hook, cleanup!, finished
export setup!, update!, hook, cleanup!, finished

using Policy
using Common
using ExpFunctions: FunctionSet

using LinearAlgebra

export PolicyGradModel
export NPGStrategy, PolicyRollout

struct PolicyGradModel
   theta::Vector     # policy parameters
   mjsys::mjWrap.mjSet
   policy::Policy.AbstractPolicy

   # scratch stuff
   samples::mjWrap.TrajSamples
   meansamples::mjWrap.TrajSamples

   trace::Dict{Symbol,Vector{Float64}}

   niter::Int

   function PolicyGradModel(theta,
                            mjsys, policy,
                            T, K, meanT, meanK, niter)
      ns = mjsys.ns
      samples     = mjw.allocateTrajSamples(mjsys, T, K, ns)
      meansamples = mjw.allocateTrajSamples(mjsys, meanT, meanK, ns)
      return new(theta, mjsys, policy,
                 samples, meansamples,
                 Dict(:stocR => Vector{Float64}(),
                      :meanR => Vector{Float64}(),
                      :evalscore => Vector{Float64}() ), niter)
   end
end

include("npg.jl")

################## utilities

using Common: makesavedir, save
export PolicyExp

# Policy Gradient experiment specific things
mutable struct PolicyExp <: LearningStrategy
   prefix::String
   expfile::String
   expname::String
   meshfile::String
   dir::String
   statefile::String
   meanfile::String
   policyfile::String
   expresults::String

   function PolicyExp(prefix::String, expfile::String,
                      expname::String, meshfile::String)
      return new(prefix, expfile, expname, meshfile, "", "", "", "")
   end
end

function setup!(s::PolicyExp, model)
   @info("SETTING UP NPG EXPERIMENT")

   poltype = Base.datatype_name(typeof(model.policy))
   dir = "$(s.prefix)/$(s.expname)_$poltype"
   expdir = makesavedir(dir, s.expfile, model.mjsys.name, s.meshfile)
   s.dir        = "$(expdir)"
   s.statefile  = "$(expdir)/data.jld2"
   s.meanfile   = "$(expdir)/mean.jld2"
   s.policyfile = "$(expdir)/policy.jld2"
   s.expresults = "$(expdir)/expmt.jld2"
end

function savepolicygrad(s::PolicyExp, model)
   @info("\tSaving traj for experiment ", s.dir)
   mjWrap.save_traj(model.mjsys, model.samples, s.statefile)
   mjWrap.save_traj(model.mjsys, model.meansamples, s.meanfile)
   Policy.save_pol(model.policy,
                   model.mjsys.skip,
                   basename(model.mjsys.name),
                   s.policyfile) # needs skip and model for policy playback (FIX LATER)
   save(s.expresults, model.trace)
end

function hook(s::PolicyExp, model, i)
   #if i==1 || model.trace[:stocR][i] >= maximum(model.trace[:stocR][1:i-1])
   N = 100
   if mod1(i, N) == N
      savepolicygrad(s, model)
   end
end

function cleanup!(s::PolicyExp, model)
   savepolicygrad(s, model)  # TODO hack may not want to do this...
   @info("Saved Experiment to ", s.dir)
   # email myself when done
   #run(`bash ./plot2html.sh $(scores[end]) $(expfile)`)
   #if cluster == 1
   #   println(workers())
   #   for w in workers()
   #      rmprocs(w)
   #   end
   #end
end

# TODO for parallel code
# setup! starts up remote do-work functions
# update! does the baseline sync and distributed CG, same as old parallelepoch

end


text:
# simple keyword args case
kwf1(ones; tens=0, hundreds=0) = ones + 10*tens + 100*hundreds
@test kwf1(2) == 2
@test kwf1(2, hundreds=6) == 602
@test kwf1(2, tens=6) == 62
@test kwf1(1, hundreds=2, tens=7) == 271
@test kwf1(3, tens=7, hundreds=2) == 273

@test_throws kwf1()        # no method, too few args
@test_throws kwf1(1, z=0)  # unsupported keyword
@test_throws kwf1(1, 2)    # no method, too many positional args

# keyword args plus varargs
kwf2(x, rest...; y=1) = (x, y, rest)
@test isequal(kwf2(0), (0, 1, ()))
@test isequal(kwf2(0,1,2), (0, 1, (1,2)))
@test isequal(kwf2(0,1,2,y=88), (0, 88, (1,2)))
@test isequal(kwf2(0,y=88,1,2), (0, 88, (1,2)))
@test_throws kwf2(0, z=1)
@test_throws kwf2(y=1)

# keyword arg with declared type
kwf3(x; y::Float64 = 1.0) = x + y
@test kwf3(2) == 3.0
@test kwf3(2, y=3.0) == 5.0
@test_throws kwf3(2, y=3)  # wrong type keyword

# function with only keyword args
kwf4(;a=1,b=2) = (a,b)
@test isequal(kwf4(), (1,2))
@test isequal(kwf4(b=10), (1,10))

# in-order evaluation of keyword args
kwf9(;read=true,write=!read) = (read,write)
@test kwf9() === (true,false)
@test kwf9(read=false) === (false,true)
@test kwf9(write=true) === (true,true)

# rest keywords
kwdelegator(ones;kw...) = kwf1(ones;kw...)
@test kwdelegator(4,hundreds=8) == 804

# optional positional args
opaf1(a,b=1,c=2,d=3) = (a,b,c,d)
@test isequal(opaf1(0), (0,1,2,3))
@test isequal(opaf1(0,2), (0,2,2,3))
@test isequal(opaf1(0,2,4), (0,2,4,3))
@test isequal(opaf1(0,2,4,6), (0,2,4,6))
@test_throws opaf1()
@test_throws opaf1(0,1,2,3,4)

# optional positional plus varargs
opaf2(a=1,rest...) = (a,rest)
@test isequal(opaf2(), (1,()))
@test isequal(opaf2(2), (2,()))
@test isequal(opaf2(2,3), (2,(3,)))

# optional positional plus keyword args
opkwf1(a=0,b=1;k=2) = (a,b,k)
@test isequal(opkwf1(), (0,1,2))
@test isequal(opkwf1(10), (10,1,2))
@test isequal(opkwf1(10,20), (10,20,2))
@test_throws opkwf1(10,20,30)
@test isequal(opkwf1(10,20,k=8), (10,20,8))
@test isequal(opkwf1(11;k=8),    (11, 1,8))
@test isequal(opkwf1(k=8),       ( 0, 1,8))

# dictionaries as keywords
@test kwf1(4; {:hundreds=>9, :tens=>5}...) == 954

# with inner function
let
    function kwf_maker()
        f(;k=0) = 2k+1
    end
    kwf5 = kwf_maker()
    @test kwf5() == 1
    @test kwf5(k=2) == 5
    @test_throws kwf5(1)
end

# with every feature!
extravagant_args(x,y=0,rest...;color="blue",kw...) =
    (x,y,rest,color,kwf1(6;tens=8,kw...))

@test isequal(extravagant_args(1), (1,0,(),"blue",86))
@test isequal(extravagant_args(1;hundreds=7), (1,0,(),"blue",786))
@test isequal(extravagant_args(1,2,3;{:color=>"red", :hundreds=>3}...),
              (1,2,(3,),"red",386))

# passing empty kw container to function with no kwargs
@test sin(1.0) == sin(1.0;{}...)

# passing junk kw container
@test_throws extravagant_args(1; {[]}...)

# keyword args with static parmeters
kwf6{T}(x; k::T=1) = T
@test kwf6(1) === Int
@test kwf6(1;k=2.5) === Float64

kwf7{T}(x::T; k::T=1) = T
@test kwf7(2) === Int
@test kwf7(1.5;k=2.5) === Float64
@test_throws kwf7(1.5)
@test_throws kwf7(1.5;k=2)

# try to confuse it with quoted symbol
kwf8{T}(x::MIME{:T};k::T=0) = 0
@test kwf8(MIME{:T}()) === 0

# issue #4538
macro TEST4538()
    quote
        function $(esc(:test4538))(x=1)
            return x
        end
    end
end
@TEST4538
@test test4538() == 1
@test test4538(2) == 2

macro TEST4538_2()
    quote
        function $(esc(:test4538_2))(;x=1)
            return x
        end
    end
end
@TEST4538_2
@test test4538_2() == 1
@test_throws test4538_2(2)
@test test4538_2(x=2) == 2

# that, but in a module
module Foo4538
macro TEST()
    quote
        function $(esc(:test4538_foo_2))(;x=1)
            return x
        end
    end
end
end
@Foo4538.TEST()
@test test4538_foo_2() == 1
@test test4538_foo_2(x=2) == 2

f4538_3(;x=1) = x
macro TEST4538_3()
    quote
        x = 2
        f4538_3(x=3)
    end
end
@test (@TEST4538_3) == 3

# issue #4801
type T4801{X}
    T4801(;k=0) = new()
end
@test isa(T4801{Any}(k=0), T4801{Any})

# issue #4974
function Foo4974(f; kwargs...)
    (f(), kwargs)
end

function test4974(;kwargs...)
    Foo4974(;kwargs...) do
        2
    end
end

@test test4974(a=1) == (2, {(:a,1)})

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
# Shared Memory (part of B.2)

export @cuStaticSharedMem, @cuDynamicSharedMem, CuStaticSharedArray, CuDynamicSharedArray

"""
    CuStaticSharedArray(T::Type, dims) -> CuDeviceArray{T,AS.Shared}

Get an array of type `T` and dimensions `dims` (either an integer length or tuple shape)
pointing to a statically-allocated piece of shared memory. The type should be statically
inferable and the dimensions should be constant, or an error will be thrown and the
generator function will be called dynamically.
"""
@inline function CuStaticSharedArray(::Type{T}, dims) where {T}
    len = prod(dims)
    # NOTE: this relies on const-prop to forward the literal length to the generator.
    #       maybe we should include the size in the type, like StaticArrays does?
    ptr = emit_shmem(T, Val(len))
    CuDeviceArray(dims, ptr)
end

macro cuStaticSharedMem(T, dims)
    Base.depwarn("@cuStaticSharedMem is deprecated, please use the CuStaticSharedArray function", :CuStaticSharedArray)
    quote
        CuStaticSharedArray($(esc(T)), $(esc(dims)))
    end
end

"""
    CuDynamicSharedArray(T::Type, dims, offset::Integer=0) -> CuDeviceArray{T,AS.Shared}

Get an array of type `T` and dimensions `dims` (either an integer length or tuple shape)
pointing to a dynamically-allocated piece of shared memory. The type should be statically
inferable or an error will be thrown and the generator function will be called dynamically.

Note that the amount of dynamic shared memory needs to specified when launching the kernel.

Optionally, an offset parameter indicating how many bytes to add to the base shared memory
pointer can be specified. This is useful when dealing with a heterogeneous buffer of dynamic
shared memory; in the case of a homogeneous multi-part buffer it is preferred to use `view`.
"""
@inline function CuDynamicSharedArray(::Type{T}, dims, offset=0) where {T}
    len = prod(dims)
    @boundscheck if offset+len > dynamic_smem_size()
        throw(BoundsError())
    end
    ptr = emit_shmem(T) + offset
    CuDeviceArray(dims, ptr)
end

macro cuDynamicSharedMem(T, dims, offset=0)
    Base.depwarn("@cuDynamicSharedMem is deprecated, please use the CuDynamicSharedArray function", :CuStaticSharedArray)
    quote
        CuDynamicSharedArray($(esc(T)), $(esc(dims)), $(esc(offset)))
    end
end

dynamic_smem_size() = @asmcall("mov.u32 \$0, %dynamic_smem_size;", "=r", true, UInt32, Tuple{})

# get a pointer to shared memory, with known (static) or zero length (dynamic shared memory)
@generated function emit_shmem(::Type{T}, ::Val{len}=Val(0)) where {T,len}
    Context() do ctx
        eltyp = convert(LLVMType, T; ctx)
        T_ptr = convert(LLVMType, LLVMPtr{T,AS.Shared}; ctx)

        # create a function
        llvm_f, _ = create_function(T_ptr)

        # create the global variable
        mod = LLVM.parent(llvm_f)
        gv_typ = LLVM.ArrayType(eltyp, len)
        gv = GlobalVariable(mod, gv_typ, "shmem", AS.Shared)
        if len > 0
            # static shared memory should be demoted to local variables, whenever possible.
            # this is done by the NVPTX ASM printer:
            # > Find out if a global variable can be demoted to local scope.
            # > Currently, this is valid for CUDA shared variables, which have local
            # > scope and global lifetime. So the conditions to check are :
            # > 1. Is the global variable in shared address space?
            # > 2. Does it have internal linkage?
            # > 3. Is the global variable referenced only in one function?
            linkage!(gv, LLVM.API.LLVMInternalLinkage)
            initializer!(gv, null(gv_typ))
        end
        # by requesting a larger-than-datatype alignment, we might be able to vectorize.
        # we pick 32 bytes here, since WMMA instructions require 32-byte alignment.
        # TODO: Make the alignment configurable
        alignment!(gv, max(32, Base.datatype_alignment(T)))

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            ptr = gep!(builder, gv, [ConstantInt(0; ctx), ConstantInt(0; ctx)])

            untyped_ptr = bitcast!(builder, ptr, T_ptr)

            ret!(builder, untyped_ptr)
        end

        call_function(llvm_f, LLVMPtr{T,AS.Shared})
    end
end


text:
using Oceananigans.Coriolis: Ω_Earth
using Oceananigans.Advection: EnergyConservingScheme, EnstrophyConservingScheme

function instantiate_fplane_1(FT)
    coriolis = FPlane(FT, f=π)
    return coriolis.f == FT(π)
end

function instantiate_fplane_2(FT)
    coriolis = FPlane(FT, rotation_rate=2, latitude=30)
    return coriolis.f == FT(2)
end

function instantiate_constant_coriolis_1(FT)
    coriolis = ConstantCartesianCoriolis(FT, f=1, rotation_axis=[0, cosd(45), sind(45)])
    @test coriolis.fy == FT(cosd(45))
    @test coriolis.fz == FT(sind(45))
end

function instantiate_constant_coriolis_2(FT)
    coriolis = ConstantCartesianCoriolis(FT, f=10, rotation_axis=[√(1/3),√(1/3),√(1/3)])
    @test coriolis.fx == FT(10*√(1/3))
    @test coriolis.fy == FT(10*√(1/3))
    @test coriolis.fz == FT(10*√(1/3))
end

function instantiate_betaplane_1(FT)
    coriolis = BetaPlane(FT, f₀=π, β=2π)
    @test coriolis.f₀ == FT(π)
    @test coriolis.β  == FT(2π)
end

function instantiate_betaplane_2(FT)
    coriolis = BetaPlane(FT, latitude=70, radius=2π, rotation_rate=3π)
    @test coriolis.f₀ == FT(6π * sind(70))
    @test coriolis.β  == FT(6π * cosd(70) / 2π)
end

function instantiate_ntbetaplane_1(FT)
    coriolis = NonTraditionalBetaPlane(FT, fz=π, fy=ℯ, β=1//7, γ=5)
    @test coriolis.fz == FT(π)
    @test coriolis.fy == FT(ℯ)
    @test coriolis.β  == FT(1//7)
    @test coriolis.γ  == FT(5)
end

function instantiate_ntbetaplane_2(FT)
    Ω, φ, R = π, 17, ℯ
    coriolis = NonTraditionalBetaPlane(FT, rotation_rate=Ω, latitude=φ, radius=R)
    @test coriolis.fz == FT(+2Ω*sind(φ))
    @test coriolis.fy == FT(+2Ω*cosd(φ))
    @test coriolis.β  == FT(+2Ω*cosd(φ)/R)
    @test coriolis.γ  == FT(-4Ω*sind(φ)/R)
end

function instantiate_hydrostatic_spherical_coriolis1(FT)
    coriolis = HydrostaticSphericalCoriolis(FT, scheme=EnergyConservingScheme())
    @test coriolis.rotation_rate == FT(Ω_Earth)
    @test coriolis.scheme isa EnergyConservingScheme

    coriolis = HydrostaticSphericalCoriolis(FT, scheme=EnstrophyConservingScheme())
    @test coriolis.rotation_rate == FT(Ω_Earth)
    @test coriolis.scheme isa EnstrophyConservingScheme
end

function instantiate_hydrostatic_spherical_coriolis2(FT)
    coriolis = HydrostaticSphericalCoriolis(FT, rotation_rate=π)
    @test coriolis.rotation_rate == FT(π)
    @test coriolis.scheme isa EnergyConservingScheme # default
end

@testset "Coriolis" begin
    @info "Testing Coriolis..."

    @testset "Coriolis" begin
        for FT in float_types
            @test instantiate_fplane_1(FT)
            @test instantiate_fplane_2(FT)

            instantiate_constant_coriolis_1(FT)
            instantiate_constant_coriolis_2(FT)
            instantiate_betaplane_1(FT)
            instantiate_betaplane_2(FT)
            instantiate_hydrostatic_spherical_coriolis1(FT)
            instantiate_hydrostatic_spherical_coriolis2(FT)

            # Test that FPlane throws an ArgumentError
            @test_throws ArgumentError FPlane(FT)
            @test_throws ArgumentError FPlane(FT, rotation_rate=7e-5)
            @test_throws ArgumentError FPlane(FT, f=1, latitude=40)
            @test_throws ArgumentError FPlane(FT, f=1, rotation_rate=7e-5, latitude=40)

            # Test that ConstantCartesianCoriolis throws an ArgumentError
            @test_throws ArgumentError ConstantCartesianCoriolis(FT, rotation_axis=[0,1,1])
            @test_throws ArgumentError ConstantCartesianCoriolis(FT, f=1, latitude=45,)
            @test_throws ArgumentError ConstantCartesianCoriolis(FT, fx=1, latitude=45)
            @test_throws ArgumentError ConstantCartesianCoriolis(FT, fx=1, f=1)

            # Non-exhaustively test that BetaPlane throws an ArgumentError
            @test_throws ArgumentError BetaPlane(FT)
            @test_throws ArgumentError BetaPlane(FT, f₀=1)
            @test_throws ArgumentError BetaPlane(FT, β=1)
            @test_throws ArgumentError BetaPlane(FT, f₀=1e-4, β=1e-11, latitude=70)

            # Test that NonTraditionalBetaPlane throws an ArgumentError
            @test_throws ArgumentError NonTraditionalBetaPlane(FT)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, rotation_rate=7e-5)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fy=1, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fy=1, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, β=3, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, β=3, rotation_rate=7e-5, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, β=3, γ=4, latitude=40)
            @test_throws ArgumentError NonTraditionalBetaPlane(FT, fz=1, fy=2, β=3, γ=4, rotation_rate=7e-5, latitude=40)

            # Test show functions
            ✈ = FPlane(FT, latitude=45)
            show(✈); println()
            @test ✈ isa FPlane{FT}

            ✈ = ConstantCartesianCoriolis(FT, f=1e-4)
            show(✈); println()
            @test ✈ isa ConstantCartesianCoriolis{FT}

            ✈ = BetaPlane(FT, latitude=45)
            show(✈); println()
            @test ✈ isa BetaPlane{FT}

            ✈ = NonTraditionalBetaPlane(FT, latitude=45)
            show(✈); println()
            @test ✈ isa NonTraditionalBetaPlane{FT}
        end
    end
end


text:
export Tetrahedron, Dodecahedron, Icosahedron, Octahedron

function add_edge_matrix!(G::SimpleGraph, edges::Array{Int,2})
    ne = size(edges,1)
    for j=1:ne
        e = edges[j,:]
        add!(G,e[1],e[2])
    end
    return
end

"""
`Tetrahedron()` creates the tetrahedron `SimpleGraph`.
"""
Tetrahedron() = Complete(4)

"""
`Dodecahedron()` creates the dodecahedron `SimpleGraph`.
"""
function Dodecahedron()
    G = IntGraph()
    edges = [
             1 2
             1 11
             1 20
             2 3
             2 9
             3 4
             3 7
             4 5
             4 20
             5 6
             5 18
             6 7
             6 16
             7 8
             8 9
             8 15
             9 10
             10 11
             10 14
             11 12
             12 13
             12 19
             13 14
             13 17
             14 15
             15 16
             16 17
             17 18
             18 19
             19 20
             ]
    add_edge_matrix!(G,edges)
    return G
end

"""
`Icosahedron()` creates the icosahedron `SimpleGraph`.
"""
function Icosahedron()
    G = IntGraph()
    edges = [
             1 2
             1 6
             1 8
             1 9
             1 12
             2 3
             2 6
             2 7
             2 9
             3 4
             3 7
             3 9
             3 10
             4 5
             4 7
             4 10
             4 11
             5 6
             5 7
             5 11
             5 12
             6 7
             6 12
             8 9
             8 10
             8 11
             8 12
             9 10
             10 11
             11 12
             ]
    add_edge_matrix!(G, edges)
    return G
end

"""
`Octahedron()` creates the octaahedron `SimpleGraph`.
"""
function Octahedron()
    G = IntGraph()
    edges = [
             1 2
             1 3
             1 4
             1 5
             2 3
             2 4
             2 6
             3 5
             3 6
             4 5
             4 6
             5 6
             ]

    add_edge_matrix!(G, edges)
    return G
end


text:
module JlrsReflect

export reflect, renamestruct!, renamefields!

import Base: show, getindex, write

abstract type Wrapper end

struct StructParameter
    name::Symbol
    elide::Bool
end

struct TypeParameter
    name::Symbol
    value
end

struct GenericWrapper <: Wrapper
    name::Symbol
end

mutable struct StructField
    name::Symbol
    rsname::String
    fieldtype::Wrapper
    typeparams::Vector{TypeParameter}
    referenced::Set{TypeVar}
    framelifetime::Bool
    datalifetime::Bool
end

struct BitsUnionWrapper <: Wrapper
    union_of::Union
    typeparams::Vector{StructParameter}
    framelifetime::Bool
    datalifetime::Bool
    BitsUnionWrapper(union_of::Union) = new(union_of, [], false, false)
end

mutable struct StructWrapper <: Wrapper
    name::Symbol
    typename::Core.TypeName
    rsname::String
    fields::Vector{StructField}
    typeparams::Vector{StructParameter}
    framelifetime::Bool
    datalifetime::Bool
end

struct TupleField
    fieldtype::Wrapper
    typeparams::Vector{TypeParameter}
    framelifetime::Bool
    datalifetime::Bool
end

struct TupleWrapper <: Wrapper
    rsname::String
    fields::Vector{TupleField}
    framelifetime::Bool
    datalifetime::Bool
    TupleWrapper(fields::Vector{TupleField}, framelifetime::Bool, datalifetime::Bool) = new(string("::jlrs::wrappers::inline::tuple::Tuple", length(fields)), fields, framelifetime, datalifetime)
end

struct BuiltinWrapper <: Wrapper
    rsname::String
    typeparams::Vector{StructParameter}
    framelifetime::Bool
    datalifetime::Bool
end

struct Wrappers
    dict::Dict{Type,Wrapper}
end

struct StringWrappers
    dict::Dict{Type,String}
end

function StringWrappers(wrappers::Wrappers)
    strwrappers = Dict{Type,String}()

    for name in keys(wrappers.dict)
        rustimpl = strwrapper(wrappers.dict[name], wrappers.dict)
        if rustimpl !== nothing
            strwrappers[name] = rustimpl
        end
    end

    StringWrappers(strwrappers)
end

function getindex(sb::StringWrappers, els...)
    sb.dict[els...]
end

function show(io::IO, wrappers::Wrappers)
    rustimpls = []
    names = []

    for name in keys(wrappers.dict)
        push!(names, name)
    end

    for name in sort(names, lt=(a, b) -> string(a) < string(b))
        rustimpl = strwrapper(wrappers.dict[name], wrappers.dict)
        if rustimpl !== nothing
            push!(rustimpls, rustimpl)
        end
    end

    print(io, join(rustimpls, "\n\n"))
end

function write(io::IO, wrappers::Wrappers)
    rustimpls = ["use jlrs::prelude::*;"]
    names = []

    for name in keys(wrappers.dict)
        push!(names, name)
    end

    for name in sort(names, lt=(a, b) -> string(a) < string(b))
        rustimpl = strwrapper(wrappers.dict[name], wrappers.dict)
        if rustimpl !== nothing
            push!(rustimpls, rustimpl)
        end
    end

    write(io, join(rustimpls, "\n\n"), "\n")
end

function insertbuiltins!(wrappers::Dict{Type,Wrapper})::Nothing
    wrappers[UInt8] = BuiltinWrapper("u8", [], false, false)
    wrappers[UInt16] = BuiltinWrapper("u16", [], false, false)
    wrappers[UInt32] = BuiltinWrapper("u32", [], false, false)
    wrappers[UInt64] = BuiltinWrapper("u64", [], false, false)
    wrappers[Int8] = BuiltinWrapper("i8", [], false, false)
    wrappers[Int16] = BuiltinWrapper("i16", [], false, false)
    wrappers[Int32] = BuiltinWrapper("i32", [], false, false)
    wrappers[Int64] = BuiltinWrapper("i64", [], false, false)
    wrappers[Float16] = BuiltinWrapper("::half::f16", [], false, false)
    wrappers[Float32] = BuiltinWrapper("f32", [], false, false)
    wrappers[Float64] = BuiltinWrapper("f64", [], false, false)
    wrappers[Bool] = BuiltinWrapper("::jlrs::wrappers::inline::bool::Bool", [], false, false)
    wrappers[Char] = BuiltinWrapper("::jlrs::wrappers::inline::char::Char", [], false, false)
    wrappers[Core.SSAValue] = BuiltinWrapper("::jlrs::wrappers::inline::ssa_value::SSAValue", [], false, false)
    wrappers[Union{}] = BuiltinWrapper("::jlrs::wrappers::inline::union::EmptyUnion", [], false, false)

    wrappers[Any] = BuiltinWrapper("::jlrs::wrappers::ptr::ValueRef", [], true, true)
    wrappers[basetype(Array)] = BuiltinWrapper("::jlrs::wrappers::ptr::ArrayRef", [StructParameter(:T, true), StructParameter(:N, true)], true, true)
    wrappers[Core.CodeInstance] = BuiltinWrapper("::jlrs::wrappers::ptr::CodeInstanceRef", [], true, false)
    wrappers[DataType] = BuiltinWrapper("::jlrs::wrappers::ptr::DataTypeRef", [], true, false)
    wrappers[Expr] = BuiltinWrapper("::jlrs::wrappers::ptr::ExprRef", [], true, false)
    wrappers[Method] = BuiltinWrapper("::jlrs::wrappers::ptr::MethodRef", [], true, false)
    wrappers[Core.MethodInstance] = BuiltinWrapper("::jlrs::wrappers::ptr::MethodInstanceRef", [], true, false)
    wrappers[Core.MethodMatch] = BuiltinWrapper("::jlrs::wrappers::ptr::MethodMatchRef", [], true, false)
    wrappers[Core.MethodTable] = BuiltinWrapper("::jlrs::wrappers::ptr::MethodTableRef", [], true, false)
    wrappers[Core.Method] = BuiltinWrapper("::jlrs::wrappers::ptr::MethodRef", [], true, false)
    wrappers[Module] = BuiltinWrapper("::jlrs::wrappers::ptr::ModuleRef", [], true, false)
    wrappers[Core.SimpleVector] = BuiltinWrapper("::jlrs::wrappers::ptr::SimpleVectorRef", [], true, false)
    wrappers[String] = BuiltinWrapper("::jlrs::wrappers::ptr::StringRef", [], true, false)
    wrappers[Symbol] = BuiltinWrapper("::jlrs::wrappers::ptr::SymbolRef", [], true, false)
    wrappers[Task] = BuiltinWrapper("::jlrs::wrappers::ptr::TaskRef", [], true, false)
    wrappers[Core.TypeName] = BuiltinWrapper("::jlrs::wrappers::ptr::TypeNameRef", [], true, false)
    wrappers[TypeVar] = BuiltinWrapper("::jlrs::wrappers::ptr::TypeVarRef", [], true, false)
    wrappers[Core.TypeMapEntry] = BuiltinWrapper("::jlrs::wrappers::ptr::TypeMapEntryRef", [], true, false)
    wrappers[Core.TypeMapLevel] = BuiltinWrapper("::jlrs::wrappers::ptr::TypeMapLevelRef", [], true, false)
    wrappers[Union] = BuiltinWrapper("::jlrs::wrappers::ptr::UnionRef", [], true, false)
    wrappers[UnionAll] = BuiltinWrapper("::jlrs::wrappers::ptr::UnionAllRef", [], true, false)
    wrappers[WeakRef] = BuiltinWrapper("::jlrs::wrappers::ptr::WeakRefRef", [], true, false)

    nothing
end

function toposort!(data::Dict{DataType,Set{DataType}})::Vector{Type}
    for (k, v) in data
        delete!(v, k)
    end

    for item in setdiff(reduce(∪, values(data)), keys(data))
        data[item] = Set()
    end

    rst = Vector()
    while true
        ordered = Set(item for (item, dep) in data if isempty(dep))
        if isempty(ordered) break end
        append!(rst, ordered)
        data = Dict(item => setdiff(dep, ordered) for (item, dep) in data if item ∉ ordered)
    end

    @assert isempty(data) "a cyclic dependency exists amongst $(keys(data))"
    rst
end

function partialtype(type::UnionAll)::DataType
    t = type

    while t.body isa UnionAll
        t = t.body
    end

    return t.body
end

function partialtype(type::DataType)::DataType
    return type
end

function basetype(type::DataType)::DataType
    definition = getproperty(type.name.module, type.name.name)
    partialtype(definition)
end

function basetype(type::UnionAll)::DataType
    partial = partialtype(type)
    definition = getproperty(partial.name.module, partial.name.name)
    partialtype(definition)
end

const BUILTINS = begin
    d = Dict{Type,Wrapper}()
    insertbuiltins!(d)
    d
end

function isnonparametric(type::Union)::Bool
    for utype in Base.uniontypes(type)
        if utype isa DataType
            if findfirst(utype.parameters) do p p isa TypeVar end !== nothing
                return false
            end

            continue
        elseif utype isa TypeVar
            return false
        end
    end

    true
end

function extracttupledeps!(acc::Dict{DataType,Set{DataType}}, type::DataType)::Nothing
    for ttype in type.types
        extractdeps!(acc, ttype)
    end

    nothing
end

function extracttupledeps!(acc::Dict{DataType,Set{DataType}}, key::DataType, type::DataType)::Nothing
    for ttype in type.types
        if ttype isa DataType
            if ttype <: Tuple
                if !isconcretetype(ttype)
                    extracttupledeps!(acc, ttype)
                else
                    extracttupledeps!(acc, key, ttype)
                end
            else
                tbase = basetype(ttype)
                push!(acc[key], tbase)
                extractdeps!(acc, ttype)
            end
        elseif ttype isa Union
            extractdeps!(acc, ttype)
        end
    end

    nothing
end

function extractdeps!(acc::Dict{DataType,Set{DataType}}, type::Type)::Nothing
    if type isa DataType
        if type <: Tuple
            return extracttupledeps!(acc, type)
        elseif isabstracttype(type)
            return
        end

        partial = partialtype(type)
        base = basetype(type)

        if !(base in keys(acc)) && !(base in keys(BUILTINS) )
            acc[base] = Set()

            for btype in base.types
                if btype isa DataType
                    if btype <: Tuple
                        if findfirst(btype.parameters) do p p isa TypeVar end !== nothing
                            error("Tuple fields with type parameters are not supported")
                        elseif !isconcretetype(btype)
                            extracttupledeps!(acc, btype)
                        else
                            extracttupledeps!(acc, type, btype)
                        end
                    else
                        bbase = basetype(btype)
                        push!(acc[base], bbase)
                        extractdeps!(acc, btype)
                    end
                elseif btype isa UnionAll
                    extractdeps!(acc, btype)
                elseif btype isa Union
                    if !isnonparametric(btype)
                        error("Unions with type parameters are not supported")
                    end
                    extractdeps!(acc, btype)
                end
            end
        end

        btypes = base.parameters
        ptypes = partial.parameters

        for i in 1:length(btypes)
            btype = btypes[i]
            ptype = ptypes[i]
            if btype isa TypeVar && ptype isa Type
                extractdeps!(acc, ptype)
            end
        end
    elseif type isa UnionAll
        extractdeps!(acc, partialtype(type))
    elseif type isa Union
        for uniontype in Base.uniontypes(type)
            if uniontype isa TypeVar
                error("Unions with type parameters are not supported")
            end

            extractdeps!(acc, uniontype)
        end
    end

    nothing
end

function extractparams(ty::Type, wrappers::Dict{Type,Wrapper})::Set{TypeVar}
    out = Set()
    if ty <: Tuple
        for elty in ty.parameters
            union!(out, extractparams(elty, wrappers))
        end

        return out
    elseif ty isa Union
        return out
    elseif isabstracttype(ty)
        return out
    end

    partial = partialtype(ty)
    base = basetype(ty)

    wrapper = wrappers[base]

    if !hasproperty(partial, :parameters)
        return out
    end

    for (name, param) in zip(wrapper.typeparams, partial.parameters)
        if !name.elide
            if param isa TypeVar
                idx = findfirst(t -> t.name == name.name, wrapper.typeparams)
                if idx !== nothing
                    push!(out, param)
                end
            elseif param isa Type
                union!(out, extractparams(param, wrappers))
            end
        end
    end

    out
end

function concretetuplefield(tuple::Type, wrappers::Dict{Type,Wrapper})::TupleWrapper
    framelifetime = false
    datalifetime = false
    fieldwrappers::Vector{TupleField} = []

    for ty in tuple.types
        fieldwrapper = if ty isa DataType
            if Base.uniontype_layout(ty)[1]
                if ty <: Tuple
                    b = concretetuplefield(ty, wrappers)
                    framelifetime |= b.framelifetime
                    datalifetime |= b.datalifetime
                    TupleField(b, [], b.framelifetime, b.datalifetime)
                else
                    bty = basetype(ty)
                    b = wrappers[bty]
                    tparams = map(a -> TypeParameter(a[1].name, a[2]), zip(bty.parameters, ty.parameters))
                    framelifetime |= b.framelifetime
                    datalifetime |= b.datalifetime
                    TupleField(b, tparams, b.framelifetime, b.datalifetime)
                end
            elseif ty in keys(wrappers)
                b = wrappers[ty]
                if b isa BuiltinWrapper
                    framelifetime |= b.framelifetime
                    datalifetime |= b.datalifetime
                    TupleField(b, [], b.framelifetime, b.datalifetime)
                else
                    framelifetime = true
                    datalifetime = true
                    TupleField(wrappers[Any], [], true, true)
                end
            else
                framelifetime = true
                datalifetime = true
                TupleField(wrappers[Any], [], true, true)
            end
        else
            error("Invalid type")
        end

        push!(fieldwrappers, fieldwrapper)
    end

    TupleWrapper(fieldwrappers, framelifetime, datalifetime)
end

function structfield(fieldname::Symbol, fieldtype::Union{Type,TypeVar}, wrappers::Dict{Type,Wrapper})::StructField
    if fieldtype isa TypeVar
        StructField(fieldname, string(fieldname), GenericWrapper(fieldtype.name), [TypeParameter(fieldtype.name, fieldtype)], Set([fieldtype]), false, false)
    elseif fieldtype isa UnionAll
        bt = basetype(fieldtype)

        if bt isa Union
            error("Unions with type parameters are not supported")
        elseif bt.name.name == :Array
            fieldwrapper = wrappers[bt]
            tparams = map(a -> TypeParameter(a[1].name, a[2]), zip(bt.parameters, bt.parameters))
            references = extractparams(bt, wrappers)
            StructField(fieldname, string(fieldname), fieldwrapper, tparams, references, fieldwrapper.framelifetime, fieldwrapper.datalifetime)
        else
            StructField(fieldname, string(fieldname), wrappers[Any], [], Set(), true, true)
        end
    elseif fieldtype isa Union
        if Base.isbitsunion(fieldtype)
            StructField(fieldname, string(fieldname), BitsUnionWrapper(fieldtype), [], Set(), false, false)
        else
            StructField(fieldname, string(fieldname), wrappers[Any], [], Set(), true, true)
        end
    elseif fieldtype == Union{}
        StructField(fieldname, string(fieldname), wrappers[Union{}], [], Set(), false, false)
    elseif fieldtype <: Tuple
        params = extractparams(fieldtype, wrappers)
        if length(params) > 0
            error("Tuples with type parameters are not supported")
        elseif isconcretetype(fieldtype)
            wrapper = concretetuplefield(fieldtype, wrappers)
            StructField(fieldname, string(fieldname), wrapper, [], Set(), wrapper.framelifetime, wrapper.datalifetime)
        else
            StructField(fieldname, string(fieldname), wrappers[Any], [], Set(), true, true)
        end
    elseif fieldtype isa DataType
        bt = basetype(fieldtype)
        if bt in keys(wrappers)
            fieldwrapper = wrappers[bt]
            tparams = map(a -> TypeParameter(a[1].name, a[2]), zip(bt.parameters, fieldtype.parameters))
            references = extractparams(fieldtype, wrappers)
            StructField(fieldname, string(fieldname), fieldwrapper, tparams, references, fieldwrapper.framelifetime, fieldwrapper.datalifetime)
        elseif Base.uniontype_layout(fieldtype)[1]
            StructField(fieldname, string(fieldname), wrappers[Any], [], Set(), true, true)
        else
            error("Cannot create field wrapper")
        end
    else
        error("Unknown field type")
    end
end

function createwrapper!(wrappers::Dict{Type,Wrapper}, type::Type)::Nothing
    bt = basetype(type)

    if isdefined(Core, :OpaqueClosure) && bt <: Core.OpaqueClosure
        error("Core.OpaqueClosure is not supported")
    end

    if bt in keys(wrappers) return end

    if isabstracttype(bt)
        wrappers[bt] = wrappers[Any]
        return
    end

    fields = []
    framelifetime = false
    datalifetime = false
    typevars = Set()
    for (name, ty) in zip(fieldnames(bt), fieldtypes(bt))
        field = structfield(name, ty, wrappers)
        framelifetime |= field.framelifetime
        datalifetime |= field.datalifetime
        union!(typevars, field.referenced)
        push!(fields, field)
    end

    params = map(a -> StructParameter(a.name, !(a in typevars)), bt.parameters)
    wrappers[bt] = StructWrapper(type.name.name, type.name, string(type.name.name), fields, params, framelifetime, datalifetime)
    nothing
end

function haslifetimes(ty::Type, wrappers::Dict{Type,JlrsReflect.Wrapper})::Tuple{Bool,Bool}
    framelifetime = false

    if ty <: Tuple
        if isconcretetype(ty)
            for fty in ty.types
                framelt, datalt = haslifetimes(fty, wrappers)
                if datalt
                    return (true, true)
                end

                framelifetime |= framelt
            end
        else
            return (true, true)
        end
    else
        bt = basetype(ty)
        wrapper = wrappers[bt]

        if wrapper.datalifetime
            return (true, true)
        end

        framelifetime |= wrapper.framelifetime

        if wrapper isa StructWrapper
            for param in ty.parameters
                if param isa Type
                    framelt, datalt = haslifetimes(param, wrappers)
                    if datalt
                        return (true, true)
                    end

                    framelifetime |= framelt
                end
            end
        end
    end

    (framelifetime, false)
end

function setparamlifetimes!(wrappers::Dict{Type,JlrsReflect.Wrapper})::Nothing
    for (ty, wrapper) in wrappers
        if wrapper isa StructWrapper
            framelifetime = wrapper.framelifetime
            datalifetime = wrapper.datalifetime

            if datalifetime
                continue
            end

            for field in wrapper.fields
                for param in field.typeparams
                    if param.value !== nothing && !(param.value isa TypeVar)
                        framelt, datalt = haslifetimes(param.value, wrappers)
                        if datalt
                            framelifetime = true
                            datalifetime = true
                            break
                        end

                        framelifetime |= framelt
                    end
                end

                if wrapper.datalifetime
                    break
                end
            end

            wrapper.framelifetime = framelifetime
            wrapper.datalifetime = datalifetime
        end
    end

    nothing
end

"""
    reflect(types::Vector{<:Type})::Wrappers

Generate Rust wrappers for all types in `types` and their dependencies. The only requirement is
that these types must not contain any union or tuple fields that depend on a type parameter.
Wrappers are generated for the most general case by erasing the contents of all provided type
parameters, so you can't avoid this restriction by explicitly providing a more qualified type.
The only effect qualifying types has, is that wrappers for the used parameters will also be
generated. The wrappers will derive `Unbox` and `ValidLayout`, and `IntoJulia` if it's a
bits-type with no type parameters.

The result of this function can be written to a file, its contents will normally be a valid Rust
module.

When you use these wrappers with jlrs, these types must be available with the same path. For
example, if you generate wrappers for `Main.Bar.Baz`, this type must be available through that
exact path and not some other path like `Main.Foo.Bar.Baz`.

# Example
```jldoctest
julia> using JlrsReflect

julia> reflect([Complex])
#[repr(C)]
#[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
#[jlrs(julia_type = "Base.Complex")]
pub struct Complex<T>
where
    T: ::jlrs::layout::valid_layout::ValidLayout + Clone,
{
    pub re: T,
    pub im: T,
}
```
"""
function reflect(types::Vector{<:Type})::Wrappers
    deps = Dict{DataType,Set{DataType}}()
    for ty in types
        extractdeps!(deps, ty)
    end

    wrappers = Dict{Type,Wrapper}()
    insertbuiltins!(wrappers)

    for ty in toposort!(deps)
        createwrapper!(wrappers, ty)
    end

    setparamlifetimes!(wrappers)
    Wrappers(wrappers)
end

function strgenerics(wrapper::StructWrapper)::Union{Nothing,String}
    generics = []
    wheres = []
    if wrapper.framelifetime
        push!(generics, "'frame")
    end

    if wrapper.datalifetime
        push!(generics, "'data")
    end

    for param in wrapper.typeparams
        if !param.elide
            push!(generics, string(param.name))
            push!(wheres, string("    ", param.name, ": ::jlrs::layout::valid_layout::ValidLayout + Clone,"))
        end
    end

    if length(generics) > 0
        wh = if length(wheres) > 0
            string("\nwhere\n", join(wheres), "\n")
        else
            " "
        end
        string("<", join(generics, ", "), ">", wh)
    end
end

function strsignature(ty::DataType, wrappers::Dict{Type,Wrapper})::String
    if ty <: Tuple
        generics = []

        for ty in ty.types
            push!(generics, strsignature(ty, wrappers))
        end

        name = string("::jlrs::wrappers::inline::tuple::Tuple", length(generics))

        if length(generics) > 0
            return string(name, "<", join(generics, ", "), ">")
        else
            return name
        end
    end

    base = basetype(ty)
    wrapper = wrappers[base]

    name = wrapper.rsname

    generics = []
    if wrapper.framelifetime
        push!(generics, "'frame")
    end

    if wrapper.datalifetime
        push!(generics, "'data")
    end

    for (tparam, param) in zip(wrapper.typeparams, ty.parameters)
        if !tparam.elide
            if param isa TypeVar
                idx = findfirst(a -> a.name == param.name, wrapper.typeparams)
                if idx !== nothing
                    push!(generics, string(param.name))
                end
            elseif param isa DataType
                push!(generics, strsignature(param, wrappers))
            end
        end
    end

    if length(generics) > 0
        string(name, "<", join(generics, ", "), ">")
    else
        name
    end
end

function strsignature(wrapper::StructWrapper, field::Union{StructField,TupleField}, wrappers::Dict{Type,Wrapper})::String
    if field.fieldtype isa GenericWrapper
        return string(field.fieldtype.name)
    elseif field.fieldtype isa TupleWrapper
        return strtuplesignature(wrapper, field, wrappers)
    end

    generics = []

    if field.framelifetime
        push!(generics, "'frame")
    end

    if field.datalifetime
        push!(generics, "'data")
    end

    for (sparam, tparam) in zip(field.fieldtype.typeparams, field.typeparams)
        if !sparam.elide
            if tparam.value isa TypeVar
                idx = findfirst(a -> a.name == tparam.value.name, wrapper.typeparams)
                if idx !== nothing
                    push!(generics, string(tparam.value.name))
                end
            elseif tparam.value isa DataType
                push!(generics, strsignature(tparam.value, wrappers))
            end
        end
    end

    if length(generics) > 0
        string(field.fieldtype.rsname, "<", join(generics, ", "), ">")
    else
        field.fieldtype.rsname
    end
end

function strtuplesignature(wrapper::StructWrapper, field::Union{StructField,TupleField}, wrappers::Dict{Type,Wrapper})::String
    generics = []

    for fieldwrapper in field.fieldtype.fields
        push!(generics, strsignature(wrapper, fieldwrapper, wrappers))
    end

    if length(generics) > 0
        string(field.fieldtype.rsname, "<", join(generics, ", "), ">")
    else
        field.fieldtype.rsname
    end
end

function strstructname(wrapper::StructWrapper)::String
    generics = strgenerics(wrapper)
    if generics !== nothing
        string(wrapper.rsname, generics)
    else
        string(wrapper.rsname, " ")
    end
end

function strstructfield(wrapper::StructWrapper, field::StructField, wrappers::Dict{Type,Wrapper})::String
    if field.fieldtype isa BitsUnionWrapper
        align_field_name = string("_", field.rsname, "_align")
        flag_field_name = string(field.rsname, "_flag")

        ibu, sz, al = Base.uniontype_layout(field.fieldtype.union_of)
        @assert ibu "Not a bits union. This should never happen, please file a bug report."

        alignment = if al == 1
            "::jlrs::wrappers::inline::union::Align1"
        elseif al == 2
            "::jlrs::wrappers::inline::union::Align2"
        elseif al == 4
            "::jlrs::wrappers::inline::union::Align4"
        elseif al == 8
            "::jlrs::wrappers::inline::union::Align8"
        elseif al == 16
            "::jlrs::wrappers::inline::union::Align16"
        else
            error("Unsupported alignment")
        end

        string(
            "    #[jlrs(bits_union_align)]\n",
            "    ", align_field_name, ": ", alignment, ",\n",
            "    #[jlrs(bits_union)]\n",
            "    pub ", field.rsname, ": ::jlrs::wrappers::inline::union::BitsUnion<", sz, ">,\n",
            "    #[jlrs(bits_union_flag)]\n",
            "    pub ", flag_field_name, ": u8,",
        )
    else
        sig = strsignature(wrapper, field, wrappers)
        string("    pub ", field.rsname, ": ", sig, ",")
    end
end

strwrapper(::BuiltinWrapper, ::Dict{Type,Wrapper})::Union{Nothing,String} = nothing

function strwrapper(wrapper::StructWrapper, wrappers::Dict{Type,Wrapper})::Union{Nothing,String}
    ty = getproperty(wrapper.typename.module, wrapper.typename.name)
    isbits = ty isa DataType && findfirst(ty.parameters) do p p isa TypeVar end === nothing && isbitstype(ty)
    intojulia = isbits ? ", IntoJulia" : ""
    zst = isbits && ty.size == 0 ? ", zero_sized_type" : ""

    modname = string(wrapper.typename.module)
    if startswith(modname, "Main.__doctest__")
        modname = "Main"
    end

    parts = [
        "#[repr(C)]",
        string("#[derive(Clone, Debug, Unbox, ValidLayout, Typecheck", intojulia, ")]"),
        string("#[jlrs(julia_type = \"", modname, ".", wrapper.typename.name, "\"", zst, ")]"),
        string("pub struct ", strstructname(wrapper), "{")
    ]
    for field in wrapper.fields
        push!(parts, strstructfield(wrapper, field, wrappers))
    end
    push!(parts, "}")
    join(parts, "\n")
end

"""
    renamestruct!(wrappers::Wrappers, type::Type, rename::String)

Change a struct's name. This can be useful if the name of a struct results in invalid Rust code or
causes warnings.

# Example
```jldoctest
julia> using JlrsReflect

julia> struct Foo end

julia> wrappers = reflect([Foo]);

julia> renamestruct!(wrappers, Foo, "Bar")

julia> wrappers
#[repr(C)]
#[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
#[jlrs(julia_type = "Main.Foo", zero_sized_type)]
pub struct Bar {
}
```
"""
function renamestruct!(wrappers::Wrappers, type::Type, rename::String)::Nothing
    btype::DataType = basetype(type)
    wrappers.dict[btype].rsname = rename

    nothing
end

"""
    renamefields!(wrappers::Wrappers, type::Type, rename::Dict{Symbol,String})
    renamefields!(wrappers::Wrappers, type::Type, rename::Vector{Pair{Symbol,String})

Change some field names of a struct. This can be useful if the name of a struct results in invalid
Rust code or causes warnings.

# Example
```jldoctest
julia> using JlrsReflect

julia> struct Food burger::Bool end

julia> wrappers = reflect([Food]);

julia> renamefields!(wrappers, Food, [:burger => "hamburger"])

julia> wrappers
#[repr(C)]
#[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
#[jlrs(julia_type = "Main.Food")]
pub struct Food {
    pub hamburger: ::jlrs::wrappers::inline::bool::Bool,
}

```
"""
function renamefields! end

function renamefields!(wrappers::Wrappers, type::Type, rename::Dict{Symbol,String})::Nothing
    btype::DataType = basetype(type)
    for field in wrappers.dict[btype].fields
        if field.name in keys(rename)
            field.rsname = rename[field.name]
        end
    end

    nothing
end

function renamefields!(wrappers::Wrappers, type::Type, rename::Vector{Pair{Symbol,String}})::Nothing
    renamefields!(wrappers, type, Dict(rename))
end
end


text:
# Function for running Hamiltonian Monte Carlo (HMC)
function hmc(model::Model, opts::HmcOpts)
  mcmc = Array(Float64, opts.mcmc.nPostBurnin, model.nPars)
  z = Array(Float64, opts.mcmc.nPostBurnin, model.nPars)
  proposed, accepted = 0., 0.

  println("Running burn-in iterations...")

  currentPars = model.randPrior()
  currentLogPosterior = model.logPosterior(currentPars)
  currentGradLogPosterior = model.gradLogPosterior(currentPars)
  invMass = inv(opts.mass)
  leapStep = opts.setLeapStep(1, 0., opts.mcmc.n, opts.mcmc.nBurnin, 0.)

  for i = 1:opts.mcmc.n
    proposed += 1

    proposedPars = copy(currentPars)
    momentum = opts.mass'*randn(model.nPars)
    currentHamiltonian = -currentLogPosterior+(momentum'*(invMass*momentum))/2
    randomSteps = ceil(rand()*opts.nLeaps)
    proposedGradLogPosterior = copy(currentGradLogPosterior)
    for j = 1:randomSteps
      momentum = momentum+(leapStep/2)*proposedGradLogPosterior
      proposedPars = proposedPars+leapStep*momentum
      proposedGradLogPosterior = model.gradLogPosterior(proposedPars)
      momentum = momentum+(leapStep/2)*proposedGradLogPosterior
    end
    proposedLogPosterior = model.logPosterior(proposedPars)
    proposedHamiltonian = -proposedLogPosterior+(momentum'*(invMass*momentum))/2

    ratio = (currentHamiltonian-proposedHamiltonian)[1]

    if ratio > 0 || (ratio > log(rand()))
      accepted += 1

      currentPars = copy(proposedPars)
      currentLogPosterior = copy(proposedLogPosterior)
      currentGradLogPosterior = copy(proposedGradLogPosterior)
    end

    if i > opts.mcmc.nBurnin
      mcmc[i-opts.mcmc.nBurnin, :] = currentPars
      z[i-opts.mcmc.nBurnin, :] = -currentGradLogPosterior/2
    end

    if mod(i, opts.mcmc.monitorRate) == 0
      acceptanceRatio = accepted/proposed

    leapStep = opts.setLeapStep(i, acceptanceRatio, opts.mcmc.n,
      opts.mcmc.nBurnin, leapStep)

      println("Iteration $i of $(opts.mcmc.n): ", round(100*acceptanceRatio, 2),
        " % acceptance ratio")

      proposed, accepted = 0., 0.
    end

    if i == opts.mcmc.nBurnin
      println("Burn-in completed...\n\nRunning post burn-in MCMC...");
    end
  end

  return mcmc, z
end


text:
#module Plot

import Gadfly
import DataFrames: DataFrame

export plot

function plot(df::DataFrame, xs::Symbol, ys::Array{Symbol, 1}; scale = :linear, filename=None, height=3, width=6)
    dflong = vcat(map(l->stack(df,l,xs),ys)...)
    if scale ==:log
        p = Gadfly.plot(dflong,x=xs,y=:value,color=:variable,Gadfly.Scale.y_log10)
    else
        p = Gadfly.plot(dflong,x=xs,y=:value,color=:variable)
    end
    if !(filename==None)
        println("saving figure in $filename")
        Gadfly.draw(Gadfly.PDF(filename, width*Gadfly.inch, height*Gadfly.inch), p)
    end
    return p
end

#end # module


text:
@testset "Parametrized Algebraic Decision Diagrams" begin
    ADD = AlgebraicDecisionDiagrams # Alias
    MLE = MultilinearExpression # Alias
    @show x1 = indicator(MLE,1) # equivalent to
    #@show x1 = Node(1, MLE(0.0), MLE(1.0))
    @show x2 = indicator(MLE,2) # equivalent to
    #@show x2 = Node(2, MLE(0.0), MLE(1.0))
    @show x3 = x1 + x2
    @show d1 = Node(1, 1 + MLE(-1.0,1), MLE(1,1))
    @show m1 = x1 + d1
    @show d2 = Node(2, 1 + MLE(-1.0,2), MLE(1,2))
    @show m2 = x2 + d2
    @show m3 = d1 * d2
    @show convert(MLE,1)
    @show promote(MLE(1.0,2),2)
    println()
    @info "Example from the paper"
    # Dictionary of variables to indices (use negative to avoid collision with automatic variable creation)
    @show v =  Dict("a" => -1, "b" => -2)
    # Maps each expression to a fresh monomial
    let
        cache = Dict{MLE,MLE}()
        global new_monomial
        function new_monomial(e::MLE)
            get!(cache,e,MLE(1.0,length(cache)+1))
        end
    end
    # Print out constraints
    dummy = MLE(1)
    function print_constraint(e1::MLE, e2::MLE)
        println(e1, " = ", e2)
        dummy
    end
    @show ϕA = ADD.reduce(
        Node(0,
            Node(5,
                1+MLE(-1.0,v["a"]),
                MLE(1.0,v["a"])
            ),
            Node(1,
                Node(5,
                    1+MLE(-1.0,v["a"]),
                    MLE(1.0,v["a"])
                ),
                Node(6,
                    1+MLE(-1.0,v["a"]),
                    MLE(1.0,v["a"])
                ),
            )
        )
    )
    @show ϕ5 = Node(5, MLE(0.6), MLE(0.4))
    @show μ1 = marginalize(ϕA * ϕ5, 5)
    @show μ1p = apply(new_monomial,μ1)
    # equivalent to:
    # @show μ1p = ADD.reduce(
    #     Node(0,
    #         MLE(1.0,3]),
    #         Node(1,
    #             MLE(1.0,3),
    #             Node(6,
    #                 MLE(1.0,4),
    #                 MLE(1.0,5)
    #             )
    #         )
    #     )
    # )
    @show ϕ6 = Node(6, MLE(0.1), MLE(0.9))
    @show μ2 = marginalize(μ1p * ϕ6, 6)
    @show μ2p = apply(new_monomial,μ2)
    @show ϕ1 = Node(1, MLE(0.625), MLE(0.375))
    @show μ3 = marginalize(μ2p * ϕ1, 1)
    @show μ3p = apply(new_monomial,μ3)
    @show ϕB = ADD.reduce(
        Node(0,
            Node(7,
                1+MLE(-1.0,v["b"]),
                MLE(1.0,v["b"])
            ),
            Node(8,
                1+MLE(-1.0,v["b"]),
                MLE(1.0,v["b"])
            ),
        )
    )
    @show ϕ8 = Node(8, MLE(0.8), MLE(0.2))
    @show μ4 = marginalize(μ3p * ϕB * ϕ8, 8)
    @show μ4p = apply(new_monomial,μ4)
    @show ϕ0 = Node(0, MLE(0.2), MLE(0.8))
    @show μ5 = marginalize(μ4p * ϕ0, 0)
    @show μ5p = apply(new_monomial,μ5)
    @show ϕ7 = Node(7, MLE(0.3), MLE(0.7))
    @show μ6 = marginalize(μ5p * ϕ7, 7) # = 0.7x12 + 0.3x11
    @test value(μ6) == MLE(0.7,12) + MLE(0.3,11)
    println()
    @info "Bilinear program"
    println("Maximize: $(value(μ6))")
    println("Subject to:")
    apply(print_constraint,μ1p,μ1)
    apply(print_constraint,μ2p,μ2)
    apply(print_constraint,μ3p,μ3)
    apply(print_constraint,μ4p,μ4)
    apply(print_constraint,μ5p,μ5);

end # testset


text:
"""
    perm_rows(A::AbstractMatrix, p::Vector{Int64})::AbstractMatrix

Apply a permutation to the rows of a given Matrix.
Used to create corresponding permutation matrices with `A = identity`.
"""
function perm_rows(A::AbstractMatrix, p::Vector{Int64})::AbstractMatrix
    if !isperm(p)
        throw("p is not a permutation.")
        return nothing
    end
    B = copy(A)
    for i in 1:length(p)
        B[i,:] = A[p[i],:]
    end
    return B
end

"""
    perm_cols(A::AbstractMatrix, p::Vector{Int64})::AbstractMatrix

Apply a permutation to the columns of a given Matrix.
Used to create corresponding inverse permutation matrices with `A = identity`.
"""
function perm_cols(A::AbstractMatrix, p::Vector{Int64})::AbstractMatrix
    if !isperm(p)
        throw("p is not a permutation.")
        return nothing
    end
    B = copy(A)
    for i in 1:length(p)
        B[:,i] = A[:,p[i]]
    end
    return B
end

"""
    apply_permutation(A::AbstractMatrix, p::Vector{Int64})::AbstractMatrix

Applies a given permutation-vector `p` to the given matrix `A` using the corresponding permutation-matrices.
"""
function apply_permutation(A::AbstractMatrix, p::Vector{Int64})::AbstractMatrix
    if !isperm(p)
        throw("p is not a permutation.")
        return nothing
    end
    n = length(p)
    if size(A, 1) != n || size(A,2) != n
        throw("Only square matrices allowed with the same number of rows/columns as the permutation.")
        return nothing
    end
    I = id(n)
    P = perm_rows(I, p)
    P_inv = perm_cols(I,p)
    return P * A * P_inv
end

"""
    print_permutation_cycles(labelled_list::Bijection{String, T}, permutations::Vector{Vector{Int64}}, filename_identifier::String; filename_prefix="permutation_cycles_") where T

Applies all given permutations to all given labelled elements and prints out the permutation cycles using the labels.
The output is redirected and appended to a file `filename_prefix * filename_identifier * ".txt"`, e.g. `"permutation_cycles_1.txt"`
"""
function print_permutation_cycles(labelled_list::Bijection{String, T}, permutations::Vector{Vector{Int64}}, filename_identifier::String; filename_prefix="permutation_cycles_") where T
    # construct the filename
    filename = filename_prefix * filename_identifier * ".txt"
    n = length(permutations[1])
    # start output in "append" mode
    open(filename, "a") do io
        for k in sort(String.(keys(labelled_list)))
            for i in 1:length(permutations)
                el_prev = labelled_list[k]
                found = String[]
                cond = true
                I = id(n)
                p = permutations[i]
                P = perm_rows(I, p)
                P_inv = perm_cols(I,p)

                print(io, "p[$i]: " * k)
                push!(found, k)
                while cond
                    el = P * el_prev * P_inv
                    lbl = labelled_list(el)
                    print(io, " → " * lbl)
                    if lbl in found
                        cond = false
                    end
                    push!(found, lbl)
                    el_prev = el
                end
                println(io,"")
            end
        end
    end
end


text:
module kNN
    export knn, kernelregression

    using StatsBase
    using Distances
    using NearestNeighbors
	using SmoothingKernels
    using LinearAlgebra
    using Statistics

    include("bandwidth.jl")
    include("majority_vote.jl")
    include("classifier.jl")
    include("regress.jl")
    include("regression.jl")
    include("anomaly.jl")
end


text:
# CUCB, "combinatorial upper confidence bound".
# Based on http://proceedings.mlr.press/v28/chen13a.html.

mutable struct CUCB <: Policy end

struct CUCBDetails <: PolicyDetails
  solver_time::Float64
end

function choose_action(instance::CombinatorialInstance{T}, policy::CUCB, state::State{T}; with_trace::Bool=false) where T
  if any(v == 0 for v in values(state.arm_counts))
    # There is at least one arm that has never been tried: maximise the arms that have never been tested.
    weights = Dict(arm => (state.arm_counts[arm] == 0.0) ? 1.0 : 0.0 for arm in keys(state.arm_counts))
  else
    # All arms have been seen, thus this formula makes sense (no zero-valued arm count).
    weights = Dict(arm => state.arm_average_reward[arm] + sqrt((3 * log(state.round + 1)) / (2 * state.arm_counts[arm])) for arm in keys(state.arm_counts))
  end

  t0 = time_ns()
  sol = solve_linear(instance, weights)
  t1 = time_ns()

  if with_trace
    return sol, CUCBDetails((t1 - t0) / 1_000_000_000)
  else
    return sol
  end
end


text:
include("SupportFunctions.jl")

#Define the p in the Lp minimization, previous work suggests p=0.5 results in
#equivalent results to 0 with signifigantly less time to convergence

#eps is a function that starts at 1 and converges to 0 as x goes from 1->inf
#this is needed by the IRLS algorithm to converge
function IRLS(MeasurementMatrix::Array{Float64,2},MeasuredOutput::Array{Float64,1};
              verbose=false,maxiter=1000,p=.5,eps=x->1/x^3,threshold=1e-7,debug=false)

    #identify the size of the input
    m=size(MeasurementMatrix,2)

    ######### First we need to calculate a valid solution ###############

    GuessedInput = InitSolution(MeasurementMatrix,MeasuredOutput,m)

    ###### Initialize some values used in the IRLS algorithm ############
    #Construct the weight matrix, used in the ridge regression
    wn = zeros(m,m)

    #Save a copy of the "previous" guess, which in this case is the original guess
    PrevGuess=GuessedInput

    #transpose the sampling matrix, this is used every iteration so it
    #is far more efficient to calculate it before hand.
    tMeasurementMatrix=MeasurementMatrix'

    #set the distance between iterations to infinite
    PrevDist = fill(Inf,convert.(Int, maxiter/100)+1)

    #start on iteration 1
    iteration=1

    if verbose
        print("Iteration: \n",iteration)
    end

    #assume we are converging, this will be set to false if neccessary
    converges=true

    #Begin Iterating
    while PrevDist[convert.(Int,round(iteration/100)+1)]>threshold

        #Calculate the diagonal weights for the ridge regression
        for j in 1:m
            wn[j,j]=1. / (GuessedInput[j]^2. + eps(iteration))^(p/2. - 1.)
        end

        #Record Previous guess to convergence test later
        PrevGuess=GuessedInput

        #IRLS step
        GuessedInput=wn*tMeasurementMatrix*
            pinv(MeasurementMatrix*wn*tMeasurementMatrix)*
            MeasuredOutput

        #CONVERGENCE TEST
        #Every 1% of maxiter, see if we are approaching convergence
        if mod(iteration,convert.(Int, maxiter/100))==0
            #Measure convergence as euclidean distance between answers

            PrevDist[convert.(Int,round(iteration/100)+1)] = sqrt(sum((PrevGuess-GuessedInput).^2.))/m

            #Print progress if so desired
            if verbose
                print("\r"*string(iteration, "  Euclidean Distance between steps: "
                    ,PrevDist[convert.(Int,round(iteration/100)+1)]))
            end

            #most involved convergence test, see if the distance increases from one iteration
            #to the next, and if it does, break, converges=false since we didn't hit threshold
            if (convert.(Int,round(iteration/100)+1)>2) &&
                    (PrevDist[convert.(Int,round(iteration/100)+1)] > PrevDist[max(convert.(Int, round(iteration/100)),1)])
                converges=false
                break
            end
        end

        #if we pass maxiter iterations, give up
        if iteration>=maxiter
            converges=false
            break
        end

        iteration+=1
    end

    if verbose
        print("\nDone\n")
    end

    #if debugging, return additional information, such as distance to convergence and iteration number
    if debug
        return (GuessedInput,converges,iteration-1,PrevDist[1:(convert.(Int, iteration/100)+1)])
    else
        return GuessedInput
    end
end

IRLS(MeasurementMatrix,MeasuredOutput;x...)=
            IRLS(convert.(Array{Float64,2},MeasurementMatrix),
                convert.(Array{Float64,1},MeasuredOutput);x...)


text:
using Mustache

function test_data(str, iter)
    [Dict(str => x) for x in iter]
end

failing_template = mt"{{#report}}[{{#times1}}{{x}},{{/times1}}],[{{#times2}}{{y}},{{/times2}}]{{/report}}"
test_func = () -> render(failing_template,
                         Dict("report" => [ Dict("times1" => test_data("x", 1:4), "times2" => test_data("y",4:8))]))

res = test_func()
@assert res == "[1,2,3,4,],[4,5,6,7,8,]"

failing_template = mt"{{#report}}[{{#times1}}{{x}},[{{#times2}}{{y}},{{/times2}}],{{/times1}}]{{/report}}"
test_func = () -> render(failing_template,
        Dict("report" => [Dict("times1" => test_data("x", 1:3), "times2" => test_data("y",4:8))]))
res = test_func()
@assert res == "[1,[4,5,6,7,8,],2,[4,5,6,7,8,],3,[4,5,6,7,8,],]"

failing_template = mt"{{#report}}{{/report}}[{{#times1}}{{x}},[{{#times2}}{{x}},{{/times2}}],{{/times1}}]"
test_func = () -> render(failing_template,
                          Dict("report" => [],
	                              "times1" => test_data("x", 1:3),
	                              "times2" => test_data("x", 4:8)))
                                      res = test_func()
@assert res == "[1,[4,5,6,7,8,],2,[4,5,6,7,8,],3,[4,5,6,7,8,],]"

## issue #31 -- nested sections
mutable struct Location
  lat::Float64
  lon::Float64
end
mutable struct Thing
  location::Location
  name::AbstractString
end
x = [Thing(Location(1.,2.), "name"), Thing(Location(3.,4.), "nombre")]

tpl = """
{{#:x}}
{{#.}}
{{#location}}
{{lat}}--{{lon}}
{{/location}}
{{name}}
{{/.}}
{{/:x}}
"""

Mustache.render(tpl, x=x)  == "\n\n\n1.0--2.0\n\nname\n\n\n\n\n3.0--4.0\n\nnombre\n\n\n"