text:
macro import_slider_huge()
    :(
        using ComplexVisual:
            CV_SliderContainer, cv_create_hslider, cv_setup_hslider,
            cv_anchor, CV_SliderCreateData
    )
end

"""
A (Container-)Canvas with a slider (a `CV_Math2DCanvas`) inside and its rulers.
"""
struct CV_SliderContainer{ccsT, dcbT} <: CV_2DContainer # {{{
    surface        :: Cairo.CairoSurfaceImage{UInt32}
    pixel_width    :: Int32
    pixel_height   :: Int32
    bounding_box   :: CV_Rectangle{Int32} # zero-based
    user_box       :: CV_Rectangle{Int32} # user-coordinates (result of layout)
                                          # typically nonzero-based
    can_slider     :: CV_Math2DCanvas
    can_slider_l   :: CV_2DLayoutPosition # postion of can_slider inside
                                          # container
    cc_can_slider  :: ccsT
    decorator_cb   :: dcbT
end

function cv_anchor(cont::CV_SliderContainer, name::Symbol)
    name_str = string(name)
    if startswith(name_str, "slider_")
        inner_sym = Symbol(name_str[8:end])
        return cv_local2global(cont, cont.can_slider_l,
            cv_anchor(cont.can_slider, inner_sym)...)
    else
        return cv_anchor(cont.bounding_box, name)
    end
end
# }}}

struct CV_SliderCreateData{ccsT, dcbT, contextT <: CV_2DCanvasContext
                            } <: CV_CreateData
    slider_container   :: CV_SliderContainer{ccsT, dcbT}
    container_context  :: contextT
    internal_layout    :: CV_2DLayout
end

function cv_create_hslider(pixel_width::Integer, pixel_height::Integer,
        value_min::Real, value_max::Real,
        rulers::NTuple{N, CV_Ruler};
        attach::CV_AttachType=cv_south,
        decoraction_with_layout_and_position_callback=nothing) where {N}  # {{{

    if pixel_width < 0  ||  pixel_height < 0
        cv_error("pixel_width and pixel_height must be postive")
    end
    if !isfinite(value_min)  ||  !isfinite(value_max)
        cv_error("value_min and value_max must be finite")
    end
    if !(value_min < value_max)
        cv_error("value_min must be smaller than value_max")
    end
    pixel_width, pixel_height = Int32(pixel_width), Int32(pixel_height)
    value_min, value_max = Float64(value_min), Float64(value_max)
    res = pixel_width / (value_max - value_min)
    math_height = pixel_height / res

    slider_layout = CV_2DLayout()
    can_slider = CV_Math2DCanvas(
        value_min + math_height*1im, value_max + 0.0*1im, res)
    can_slider_l = cv_add_canvas!(slider_layout, can_slider, (0,0), (0,0))
    cc_can_slider = cv_create_context(can_slider)

    decorator_cb = (decoraction_with_layout_and_position_callback !== nothing) ?
        decoraction_with_layout_and_position_callback(
            slider_layout, can_slider_l) : nothing

    if N > 0
        can_axis_l = cv_ticks_labels(slider_layout, can_slider_l, attach, rulers)
    end

    bb = slider_layout.seen_boxes.bounding_box
    slider_container = CV_SliderContainer(
        cv_create_cairo_image_surface(cv_width(bb), cv_height(bb)),
        cv_width(bb), cv_height(bb),
        CV_Rectangle(cv_height(bb), Int32(0), Int32(0), cv_width(bb)),
        bb, can_slider, can_slider_l, cc_can_slider, decorator_cb)

    cc_container = cv_create_context(slider_container; fill_with=cv_color(1, 1, 1))
    if N >0
        can_axis_l(cc_container)
    end

    return CV_SliderCreateData(slider_container, cc_container, slider_layout)
end # }}}

function cv_setup_hslider(setup::CV_SceneSetupChain,
        slider_data::CV_SliderCreateData,
        cont_l::CV_2DLayoutPosition,
        painter::CV_Painter, set_slider_value_func;
        react_to_actionpixel_update::Bool=true,
        react_to_statepixel_update::Bool=false) # {{{

    cont, cc_cont = slider_data.slider_container, slider_data.container_context

    reaction = (px, py, layout) -> begin
        resp = nothing
        can_layout = cv_get_can_layout(layout)
        cc_can_layout = cv_get_cc_can_layout(layout)
        lx, ly = cv_global2local(can_layout, cont_l, px, py)
        if 0 ≤ lx ≤ cv_width(cont_l.rectangle)  &&
                0 ≤ ly ≤ cv_height(cont_l.rectangle)
            x, y = cv_pixel2math(cont, cont.can_slider_l, lx, ly)
            resp = set_slider_value_func(x)
            cv_paint(cont.cc_can_slider, painter)
            cont.can_slider_l(cc_cont)
            cont_l(cc_can_layout)
        end
        return resp
    end

    redraw_func = layout -> begin
        cc_can_layout = cv_get_cc_can_layout(layout)
        cv_paint(cont.cc_can_slider, painter)
        cont.can_slider_l(cc_cont)
        cont_l(cc_can_layout)
        return nothing
    end

    draw_once_func = layout -> begin
        if cont.decorator_cb !== nothing
            cont.decorator_cb(cc_cont)
        end
        cv_paint(cont.cc_can_slider, painter)
        cont.can_slider_l(cc_cont)
        cc_can_layout = cv_get_cc_can_layout(layout)
        cont_l(cc_can_layout)
        return nothing
    end

    return cv_combine(setup; draw_once_func, redraw_func,
        actionpixel_update=react_to_actionpixel_update ? reaction : missing,
        statepixel_update=react_to_statepixel_update ? reaction : missing)
end # }}}

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
function getrewardtests(T::Type{Tac})
    tests = []
    for (idx,s) in enumerate(states(T))
        action_rewards = map(a->(a, get_reward(s, a)), getlegalactions(s))
        reward_tests = map(action_reward->begin
            a, r = action_reward
            RewardTest(
                idx,
                s,
                get_image(s),
                a,
                get_action_image(T, a),
                r
            )
        end, action_rewards)
        append!(tests, reward_tests)
    end
    RewardTests(tests)
end

function MuZero.makepolicytests(T::Type{Tac})
    value_fn = getvaluefn(T)
    actionvaluefn = getactionvaluefn(T; value_fn=value_fn)
    tests = []
    for (idx,s) in enumerate(getstates(T))
        action_values = map(a->begin
            value = actionvaluefn(s,a)
            (a.idx, value)
        end, getlegalactions(s))
        max_value = maximum(x->x[2], action_values)
        ideal_actions = map(first, filter(x->x[2] == max_value, action_values))
        pt = PolicyTest(
            idx,
            s,
            get_image(s),
            ideal_actions
        )
        push!(tests, pt)
    end
    PolicyTests(tests)
end

function getoptimalactionsfn(
    T::Type{Tac};
    actionvaluefn=getactionvaluefn(T)
    )
    v = Dict()
    for s in getstates(T)
        action_values = map(a->begin
            (a, actionvaluefn(s,a))
        end, getlegalactions(s))
        max_value = maximum(x->x[2], action_values)
        ideal_actions = map(first, filter(x->x[2] == max_value, action_values))
        v[s] = ideal_actions
    end
    v
end

function makerrvptests(T::Type{Tac})
    value_fn = getvaluefn(T)
    optimalactionsfn = getoptimalactionsfn(
        T;
        actionvaluefn=getactionvaluefn(T; value_fn=value_fn)
    )
    tests = []
    for (idx,s) in enumerate(getstates(T))
        for a in getlegalactions(s)
            reward_next, s_next = dynamics(s, a)
            if !isterminal(s_next)
                rrvp_test = RecurrentRVPTest(
                    "$idx-$(a.idx)",
                    s,
                    s_next,
                    get_image(s),
                    a,
                    get_action_image(T, a),
                    map(x->x.idx, optimalactionsfn[s_next]),
                    value_fn[s_next],
                    reward_next,
                )
                push!(tests, rrvp_test)
            end
        end
    end
    RecurrentRVPTests(tests)
end

# @todo - this should be a MuZero function
# (i.e. MuZero.makepolicyvaluetests)
function makepolicyvaluetests(T::Type{Tac})
    value_fn = getvaluefn(T)
    actionvaluefn = getactionvaluefn(T; value_fn=value_fn)
    tests = []
    for (idx,s) in enumerate(getstates(T))
        action_values = map(a->begin
            value = actionvaluefn(s,a)
            (a.idx, value)
        end, getlegalactions(s))
        max_value = maximum(x->x[2], action_values)
        ideal_actions = map(first, filter(x->x[2] == max_value, action_values))
        actual_value = value_fn[s]
        pt = PolicyValueTest(
            idx,
            s,
            get_image(s),
            ideal_actions,
            actual_value
        )
        push!(tests, pt)
    end
    PolicyValueTests(tests)
end

function getoptimalpolicyfn(T::Type{Tac})
    actionvaluefn = getactionvaluefn(T)
    (s)->begin
        action_values = map(a->begin
            value = actionvaluefn(s,a)
            (a.idx, value)
        end, getlegalactions(s))
        max_value = maximum(x->x[2], action_values)
        ideal_actions = map(first, filter(x->x[2] == max_value, action_values))
        idealpolicy = zeros(Float16, length(getactions(s)))
        idealpolicy[ideal_actions] .= 1/length(ideal_actions)
        idealpolicy
    end
end

function getactionvaluefn(T::Type{Tac}; value_fn=getvaluefn(T))
    action_values = Dict()
    stateactionhash(s,a) = hash("$(hash(s))-$(hash(a))")
    for s in getstates(T)
        for a in getlegalactions(s)
            reward, s_next = dynamics(s, a)
            value = reward + -1*value_fn[s_next]
            action_values[stateactionhash(s,a)] = value
        end
    end
    (s,a)->action_values[stateactionhash(s,a)]
end

function MuZero.getstates(T::Type{Tac})
    terminal(s) = !isnothing(getwinner(s)) || length(getlegalactions(s)) == 0
    nonterminal(s) = !isterminal(s)
    filter(nonterminal, getstatesplus(T))
end

using Combinatorics

function Chess()
    board = emptyboard(ChessState)
    board[1, 1:8] = getpieceidx(Chess, [
        :rook_white, :knight_white,
        :bishop_white, :king_white,
        :queen_white, :bishop_white,
        :knight_white, :rook_white,
    ])
    board[2, 1:8] = getpieceidx(Chess, [
        :pawn_white, :pawn_white,
        :pawn_white, :pawn_white,
        :pawn_white, :pawn_white,
        :pawn_white, :pawn_white,
    ])
    board[7, 1:8] = getpieceidx(Chess, [
        :pawn_black, :pawn_black,
        :pawn_black, :pawn_black,
        :pawn_black, :pawn_black,
        :pawn_black, :pawn_black,
    ])
    board[8, 1:8] = getpieceidx(Chess, [
        :rook_black, :knight_black,
        :bishop_black, :queen_black,
        :king_black, :bishop_black,
        :knight_black, :rook_black,
    ])
    player = getplayer(ChessState, :white)
    states = [ChessState(board, player)]
    Chess([],[],states)
end

function MuZero.getstatesplus(T::Type{Tac})
    toplay(pieces) = isodd(length(filter(isnothing, pieces))) ? getplayer(T, :X) : getplayer(T, :O)
    perms(pieces) = pieces |> permutations |> collect |> unique
    a0 = [nothing for _ in 1:9]
    a1 = [:X, [nothing for _ in 1:8]...]
    a2 = [:X, :O, [nothing for _ in 1:7]...]
    a3 = [:X, :O, :X, [nothing for _ in 1:6]...]
    a4 = [:X, :O, :X, :O, [nothing for _ in 1:5]...]
    a5 = [:X, :O, :X, :O, :X, [nothing for _ in 1:4]...]
    a6 = [:X, :O, :X, :O, :X, :O, [nothing for _ in 1:3]...]
    a7 = [:X, :O, :X, :O, :X, :O, :X, [nothing for _ in 1:2]...]
    a8 = [:X, :O, :X, :O, :X, :O, :X, :O, nothing]
    a9 = [:X, :O, :X, :O, :X, :O, :X, :O, :X]
    all = []
    for aX in [a0, a1, a2, a3, a4, a5, a6, a7, a8, a9]
        append!(all, perms(aX))
    end
    [TacState(getboard(T, pieces), toplay(pieces)) for pieces in all]
end

function getvaluefn(T::Type{Tac}; max_i=typemax(Int), θ=Float32(0.1))
    v = Dict()
    for s in getstatesplus(T)
        v[s]=Float32(0)
    end
    function player_action_value(state, a)
        reward, s_next = dynamics(state, a)
        reward + -1 * v[s_next]
    end
    for i in 1:max_i
        Δ, num_updated = Float32(0), 0
        for state in getstates(T)
            v_prev = v[state]
            legal = getlegalactions(state)
            state_action_values = map(a->player_action_value(state, a), legal)
            max_found = maximum(state_action_values)
            if v[state] != max_found
                v[state] = max_found
                num_updated += 1
            end
            Δ = max(Δ, abs(v_prev - v[state]))
        end
        vs = v |> values |> collect
        pos = vs[vs .> 0] |> length
        neg = vs[vs .< 0] |> length
        #println("Δ=$Δ($i), $pos, $neg (num_updated: $num_updated)")
        i += 1
        (Δ >= θ && i <= max_i) || break
    end
    v
end

ext: jl


text:
using LinearAlgebra: BlasInt, chkstride1
using LinearAlgebra.LAPACK: chklapackerror
using LinearAlgebra.BLAS: @blasfunc

"""
    dspev!(jobz::Symbol, uplo::Symbol, x::StridedVector{Float64})

Computes all the eigenvalues and optionally the eigenvectors of a real
symmetric `n×n` matrix `A` in packed storage. Will corrupt `x`.

Arguments:

  `jobz`: `:N` if only eigenvalues, `:V` if eigenvalues and eigenvectors

  `uplo`: `:L` if lower triangle of `A` is stored, `:U` if upper

  `x`: `A` represented as vector of the lower (upper) n*(n+1)/2 elements, packed columnwise.

Returns:

  `W,Z` if `jobz == :V` or: `W` if `jobz == :N` such that `A=Z*diagm(W)*Z'`
"""
function dspev!(jobz::Symbol, uplo::Symbol, A::StridedVector{Float64})
    chkstride1(A)
    vecN = length(A)
    n = try
        Int(sqrt(1/4+2*vecN)-1/2)
    catch
        throw(DimensionMismatch("A has length $vecN which is not N*(N+1)/2 for any integer N"))
    end
    W     = similar(A, Float64, n)
    Z     = similar(A, Float64, n, n)
    work  = Array{Float64}(undef, 1)
    lwork = BlasInt(3*n)
    info  = Ref{BlasInt}()
    work = Array{Float64}(undef, lwork)
    ccall((@blasfunc(dspev_), Base.liblapack_name), Cvoid,
          (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{Float64},
          Ptr{Float64}, Ptr{Float64}, Ptr{BlasInt}, Ptr{Float64}, Ptr{BlasInt}),
          jobz, uplo, Ref(n), A,
          W, Z, Ref(n), work, info)
    chklapackerror(info[])
    jobz == :V ? (W, Z) : W
end

function dspevV!(uplo::Symbol, A::StridedVector{Float64})
    jobz = :V
    chkstride1(A)
    vecN = length(A)
    n = try
        Int(sqrt(1/4+2*vecN)-1/2)
    catch
        throw(DimensionMismatch("A has length $vecN which is not N*(N+1)/2 for any integer N"))
    end
    W     = similar(A, Float64, n)
    Z     = similar(A, Float64, n, n)
    work  = Array{Float64}(undef, 1)
    lwork = BlasInt(3*n)
    info  = Ref{BlasInt}()
    work = Array{Float64}(undef, lwork)
    ccall((@blasfunc(dspev_), Base.liblapack_name), Cvoid,
          (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{Float64},
          Ptr{Float64}, Ptr{Float64}, Ptr{BlasInt}, Ptr{Float64}, Ptr{BlasInt}),
          jobz, uplo, Ref(n), A,
          W, Z, Ref(n), work, info)
    chklapackerror(info[])
    return W, Z
end


text:
using Test

using Qecsim.GenericModels
using Qecsim.Model
using Qecsim.BasicModels:FiveQubitCode, SteaneCode
using Qecsim.PauliTools:bsp, to_bsf

@testset "NaiveDecoder" begin
    decoder = NaiveDecoder()
    # label
    @test label(decoder) == "Naive"
    # decode
    code = FiveQubitCode()
    for error_pauli in ("IIIII", "IXIII", "IIYII", "IIIIZ", "IZYXI", "YZYXX")
        error = to_bsf(error_pauli)
        syndrome = bsp(stabilizers(code), error)
        result = decode(decoder, code, syndrome)
        recovery = result.recovery
        @test bsp(stabilizers(code), recovery) == syndrome
        @test !any(bsp(stabilizers(code), xor.(error, recovery)))
    end
    # max_qubits
    code = SteaneCode()
    error = to_bsf("IIIIIII")
    syndrome = bsp(stabilizers(code), error)
    @test_throws ArgumentError decode(NaiveDecoder(5), code, syndrome)  # max_qubits=5
    @test decode(NaiveDecoder(), code, syndrome).recovery == error  # max_qubits=10
    @test decode(NaiveDecoder(-1), code, syndrome).recovery == error  # max_qubits=unlimited
end


text:
function test_blas(AT)
    @testset "BLAS" begin
        @testset "matmul with element type $elty" for elty in (Float32, Float64, ComplexF32, ComplexF64)
            A = rand(elty, 5, 6)
            B = rand(elty, 6, 5)
            C = rand(elty, 5, 5)
            x = rand(elty, 5)
            y = rand(elty, 6)

            @test compare(*, AT, A, y)
            @test compare((t,s) -> transpose(t)*s, AT, A, x)
            @test compare((t,s) -> adjoint(t)*s  , AT, A, x)

            @test compare(*, AT, A, B)
            @test compare((t,s) -> transpose(t)*s           , AT, A, C)
            @test compare((t,s) -> t*transpose(s)           , AT, C, B)
            @test compare((t,s) -> transpose(t)*transpose(s), AT, A, B)
            @test compare((t,s) -> adjoint(t)*s             , AT, A, C)
            @test compare((t,s) -> t*adjoint(s)             , AT, C, B)
            @test compare((t,s) -> adjoint(t)*adjoint(s)    , AT, A, B)
            @test compare((t,s) -> transpose(t)*adjoint(s)  , AT, A, B)
            @test compare((t,s) -> adjoint(t)*transpose(s)  , AT, A, B)
        end

        for T in (ComplexF32, Float32)
            @testset "rmul! $T" begin
                @test compare(rmul!, AT, rand(T, 13, 23), Ref(77f0))
            end
        end

        @testset "gbmv" begin
            m, n = 10, 11
            A, x, y = randn(Float32, 3, n), randn(Float32, n), fill(0f0, m)

            Ag, xg, yg = AT(A), AT(x), AT(y)

            BLAS.gbmv!('N', m, 1, 1, 1f0, A, x, 0f0, y)
            BLAS.gbmv!('N', m, 1, 1, 1f0, Ag, xg, 0f0, yg)
            @test y ≈ Array(yg)
        end
    end
end


text:
export  excelInitPopReport,
        excelInitPopAgeReport

function excelInitPopReport( mpSim::MPsim;
    filename::AbstractString="initPopReport", overwrite::Bool=true,
    timeFactor::Real=12.0 )

    filename, overwrite = setupFile( filename, overwrite )
    tStart = now()
    report = initPopReport( mpSim )
    reportGenerationTime = ( now() - tStart ).value / 1000.0

    if isempty( report )
        return
    end  # if isempty( report )

    writeInitialPopData( report, mpSim, timeFactor, reportGenerationTime,
        filename, overwrite )

end   # excelInitPopReport( mpSim, filename, overwrite, timeFactor )

function excelInitPopAgeReport( mpSim::MPsim, ageRes::Real, ageType::Symbol;
    filename::AbstractString="initPopReport", overwrite::Bool=true,
    timeFactor::Real=12.0 )

    if !haskey( ageTypesS, ageType )
        return
    end  # if !haskey( ageTypesS, ageType )

    filename, overwrite = setupFile( filename, overwrite )
    tStart = now()
    report = initPopAgeReport( mpSim, ageRes, ageType )
    reportGenerationTime = ( now() - tStart ).value / 1000.0

    if isempty( report )
        return
    end  # if isempty( report )

    writeInitialPopAgeData( report, mpSim, ageType, timeFactor,
        reportGenerationTime, filename, overwrite )
    #=
    tStart = now()
    filename = string( filename, endswith( filename, ".xlsx" ) ? "" : ".xlsx" )

    XLSX.openxlsx( filename, mode = overwrite ? "w" : "rw" ) do xf
        dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
        ws = xf[1]
        nReports = 1

        # Prep the new sheet.
        if overwrite
            XLSX.rename!( ws, string( "Report 1  ", dateStr ) )
        else
            nReports = length( filter( name -> startswith( name, "Report" ),
                XLSX.sheetnames( xf ) ) ) + 1
            ws = XLSX.addsheet!( xf, string( "Report ", nReports, "  ",
                dateStr ) )
        end  # if overwrite

        # Sheet header.
        ws["A1"] = "Simulation length"
        ws["B1"] = mpSim.simLength / timeFactor
        ws["C1"] = "years"
        ws["A2"] = "Report on"
        ws["B2"] = ageTypesS[ageType]
        ws["A3"] = "Data generation time"
        ws["B3"] = reportGenerationTime
        ws["C3"] = "seconds"
        ws["A4"] = "Excel generation time"
        ws["C4"] = "seconds"

        if isempty( report )
            ws["A6"] = "No initial population in simulation"
            excelGenerationTime = ( now() - tStart ).value / 1000.0
            ws["B4"] = excelGenerationTime
            return
        end  # if isempty( report )

        ws["A5"] = "Initial population size"
        ws["B5"] = sum( report[1:(end-5), 2] )

        ws["A7"] = "Summary statistics"
        ws["A8"] = "Mean"
        ws["A9"] = "Median"
        ws["A10"] = "Std. Dev."
        ws["A11"] = "Minimum"
        ws["A12"] = "Maximum"

        for ii in 1:5
            val = report[end - (5-ii), 2]
            ws[ii + 7, 2] = val isa Missing ? "N/A" : val / timeFactor
        end  # for ii in 1:5

        ws["A14"] = "Age"
        ws["A15"] = "Amount"
        nBrackets = size( report, 1 ) - 5

        for ii in 1:nBrackets
            ws[14, ii + 1] = report[ii, 1] / timeFactor
            ws[15, ii + 1] = report[ii, 2]
        end  # for ii in 1:nBrackets

        excelGenerationTime = ( now() - tStart ).value / 1000.0
        ws["B4"] = excelGenerationTime
    end  # XLSX.openxlsx( ... ) do xf
=#
end  # excelInitPopAgeReport( mpSim, ageRes, ageType, filename, overwrite,
     #   timeFactor )

include( "private/initpop.jl" )


text:
using FDTD
using Test

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


text:
using GlobalEnergyGIS, Plots, StatsBase

CRF(r,T) = r / (1 - 1/(1+r)^T)
meandrop(x; dims=dims) = dropdims(mean(x, dims=dims), dims=dims)
sumdrop(x; dims=dims) = dropdims(sum(x, dims=dims), dims=dims)

function mean_skipNaN_dim12(xx)
	hours, regs, classes = size(xx)
	out = Vector{Float64}(undef,classes)
	for c = 1:classes
		hourtot = 0.0
		hourcount = 0
		for h = 1:hours
			regtot = 0.0
			regcount = 0
			for r = 1:regs
				x = xx[h,r,c]
				if !isnan(x)
					regtot += x
					regcount += 1
				end
			end
			if regcount != 0
				hourtot += regtot / regcount
				hourcount += 1
			end
		end
		out[c] = hourtot / hourcount
	end
	return out
end

# for plant_area in [0.03, 0.12], persons_per_km2 in [75, 150]

# end

function supplycurves_pv(reg, plant_area, persons_per_km2, minclasses, step)
	cfr,cfpva,cfpvb,cfcsa,cfcsb,pvr,pva,pvb,csa,csb =
		GISsolar(savetodisk=false, gisregion=reg, pvclasses_min=minclasses, pvclasses_max=(minclasses.+step),
			plant_area=plant_area, plant_persons_per_km2=persons_per_km2)

	# results=loadresults("regionset=eurasia21, carboncap=0",resultsfile="results_prev.jld2");
	# sum(results.params[:demand][1:8,:]), sum(results.params[:demand][16:21,:])
	regdemand = Dict("Europe8" => 4.4046006293422505e6, "China6" => 1.2380343090170678e7)		# GWh/year
	demand = regdemand[reg]
	meanCF_a = mean_skipNaN_dim12(cfpva)
	meanCF_b = mean_skipNaN_dim12(cfpvb)
	meanCF_r = mean_skipNaN_dim12(cfr)
	capacity_a = sumdrop(pva, dims=1)
	capacity_b = sumdrop(pvb, dims=1)
	capacity_r = sumdrop(pvr, dims=1)
	# investcost = Dict(:pv => 800, :wind => 1200)
	# fixedcost = Dict(:pv => 16, :wind => 36)
	# lifetime = Dict(:pv => 25, :wind => 25)
	totalcost_a = capacity_a .* (600 * CRF(0.07, 25) + 16)
	totalcost_b = capacity_b .* ((600+200) * CRF(0.07, 25) + 16)
	totalcost_r = capacity_r .* (900 * CRF(0.07, 25) + 20)
	annualenergy_a = capacity_a .* meanCF_a .* 8760
	annualenergy_b = capacity_b .* meanCF_b .* 8760
	annualenergy_r = capacity_r .* meanCF_r .* 8760
	lcoe_a = totalcost_a ./ annualenergy_a .* 1000
	lcoe_b = totalcost_b ./ annualenergy_b .* 1000
	lcoe_r = totalcost_r ./ annualenergy_r .* 1000
	demandshare_a = annualenergy_a / demand
	demandshare_b = annualenergy_b / demand
	demandshare_r = annualenergy_r / demand
	return demandshare_a, demandshare_b, demandshare_r, lcoe_a, lcoe_b, lcoe_r
end

function supplycurves_wind(reg, area_onshore, persons_per_km2, minclasses, step)
	cfa,cfb,cfoff,wa,wb,woff =
		GISwind(savetodisk=false, gisregion=reg, onshoreclasses_min=minclasses, onshoreclasses_max=(minclasses.+step),
			offshoreclasses_min=minclasses, offshoreclasses_max=(minclasses.+step),
			area_onshore=area_onshore, persons_per_km2=persons_per_km2)

	# results=loadresults("regionset=eurasia21, carboncap=0",resultsfile="results_prev.jld2");
	# sum(results.params[:demand][1:8,:]), sum(results.params[:demand][16:21,:])
	regdemand = Dict("Europe8" => 4.4046006293422505e6, "China6" => 1.2380343090170678e7)		# GWh/year
	demand = regdemand[reg]
	meanCF_a = mean_skipNaN_dim12(cfa)
	meanCF_b = mean_skipNaN_dim12(cfb)
	meanCF_off = mean_skipNaN_dim12(cfoff)
	capacity_a = sumdrop(wa, dims=1)
	capacity_b = sumdrop(wb, dims=1)
	capacity_off = sumdrop(woff, dims=1)
	totalcost_a = capacity_a .* (1200 * CRF(0.07, 25) + 43)
	totalcost_b = capacity_b .* ((1200+200) * CRF(0.07, 25) + 43)
	totalcost_off = capacity_off .* (2300 * CRF(0.07, 25) + 86)
	annualenergy_a = capacity_a .* meanCF_a .* 8760
	annualenergy_b = capacity_b .* meanCF_b .* 8760
	annualenergy_off = capacity_off .* meanCF_off .* 8760
	lcoe_a = totalcost_a ./ annualenergy_a .* 1000
	lcoe_b = totalcost_b ./ annualenergy_b .* 1000
	lcoe_off = totalcost_off ./ annualenergy_off .* 1000
	demandshare_a = annualenergy_a / demand
	demandshare_b = annualenergy_b / demand
	demandshare_off = annualenergy_off / demand
	return demandshare_a, demandshare_b, demandshare_off, lcoe_a, lcoe_b, lcoe_off
end

function plotsupply(capac, lcoe; hold=false, options...)
	s = sortperm(lcoe)
	x = cumsum(capac[s])
	y = lcoe[s]
	if hold
		plot!(x,y; options...)
	else
		plot(x,y; options...)
	end
end

function plottriplepv(cap, lc, area, pop)
	demandshare_a, demandshare_b, demandshare_r = cap[:pv,area,pop]
	lcoe_a, lcoe_b, lcoe_r = lc[:pv,area,pop]
	plotsupply(demandshare_a, lcoe_a, label="PV class A", linewidth=3)
	plotsupply(demandshare_b, lcoe_b, label="PV class B", linewidth=3, hold=true)
	display(plotsupply(demandshare_r, lcoe_r, label="PV rooftop", linewidth=3, hold=true,
					size=(1000,500), xlim=(0,1.2), ylim=(0,100), title="Solar PV",
					xlabel="Potential supply (relative to annual electricity demand in 2050)", ylabel="€/MWh",
					tickfont=16, legendfont=16, guidefont=16, titlefont=20))
end

function plottriplewind(cap, lc, area, pop)
	demandshare_wa, demandshare_wb, demandshare_off = cap[:wind,area,pop]
	lcoe_wa, lcoe_wb, lcoe_off = lc[:wind,area,pop]
	plotsupply(demandshare_wa, lcoe_wa, label="Onshore class A", linewidth=3)
	plotsupply(demandshare_wb, lcoe_wb, label="Onshore class B", linewidth=3, hold=true)
	display(plotsupply(demandshare_off, lcoe_off, label="Offshore", linewidth=3, hold=true,
					size=(1000,500), xlim=(0,1.2), ylim=(0,100), title="Wind",
					xlabel="Potential supply (relative to annual electricity demand in 2050)", ylabel="€/MWh",
					tickfont=16, legendfont=16, guidefont=16, titlefont=20))
end

function mergevars(cap, lc)
	mergecap, mergelc = Dict(), Dict()
	for tech = [:pv, :wind], area = [:low, :high], pop = [150]
		cc,ll = cap[tech,area,pop], lc[tech,area,pop]
		mergecap[tech,area,pop] = [cc[1]; cc[2]; cc[3]]
		mergelc[tech,area,pop] = [ll[1]; ll[2]; ll[3]]
	end
	return mergecap, mergelc
end

function plotmonopv(cap, lc)
	mcap, mlc = mergevars(cap, lc)
	plotsupply(mcap[:pv,:low,150], mlc[:pv,:low,150], label="area=5%", linewidth=3)
	display(plotsupply(mcap[:pv,:high,150], mlc[:pv,:high,150], label="area=10%", linewidth=3, hold=true,
					size=(1000,500), xlim=(0,1.2), ylim=(0,100), title="Solar PV (class A+B+rooftop)",
					xlabel="Potential supply (relative to annual electricity demand in 2050)", ylabel="€/MWh",
					tickfont=16, legendfont=16, guidefont=16, titlefont=20))
end

function plotmonowind(cap, lc)
	mcap, mlc = mergevars(cap, lc)
	plotsupply(mcap[:wind,:low,150], mlc[:wind,:low,150], label="area=8%", linewidth=3)
	display(plotsupply(mcap[:wind,:high,150], mlc[:wind,:high,150], label="area=16%", linewidth=3, hold=true,
					size=(1000,500), xlim=(0,1.2), ylim=(0,100), title="Wind (class A+B+offshore)",
					xlabel="Potential supply (relative to annual electricity demand in 2050)", ylabel="€/MWh",
					tickfont=16, legendfont=16, guidefont=16, titlefont=20))
end

function plotmono(cap, lc)
	mcap, mlc = mergevars(cap, lc)
	plotsupply(mcap[:pv,:low,150], mlc[:pv,:low,150], label="PV area=5%", line=(3, :solid, RGB([240,220,0]/255...)))
	plotsupply(mcap[:wind,:low,150], mlc[:wind,:low,150], label="wind area=8%", line=(3, :solid, RGB([149,179,215]/255...)), hold=true)
	plotsupply(mcap[:pv,:high,150], mlc[:pv,:high,150], label="PV area=10%", line=(3, :dash, RGB([240,220,0]/255...)), hold=true)
	display(plotsupply(mcap[:wind,:high,150], mlc[:wind,:high,150], label="wind area=16%", line=(3, :dash, RGB([149,179,215]/255...)), hold=true,
					size=(1000,500), xlim=(0,1.2), ylim=(0,100), title="Solar PV & wind (class A+B+rooftop+offshore)",
					xlabel="Potential supply (relative to annual electricity demand in 2050)", ylabel="€/MWh",
					tickfont=16, legendfont=16, guidefont=16, titlefont=20))
end

function calcvars(reg)
	cap, lc = Dict(), Dict()
	area = :low
	pop = 150

	step = 0.00005
	plantarea = area == :low ? 0.05 : 0.10
	demandshare_a, demandshare_b, demandshare_r, lcoe_a, lcoe_b, lcoe_r = supplycurves_pv(reg, plantarea, pop, [0:step:0.311-step;], step)
	cap[:pv,area,pop] = (demandshare_a, demandshare_b, demandshare_r)
	lc[:pv,area,pop] = (lcoe_a, lcoe_b, lcoe_r)

	onshorearea = area == :low ? 0.08 : 0.16
	step = 0.002
	demandshare_wa, demandshare_wb, demandshare_off, lcoe_wa, lcoe_wb, lcoe_off = supplycurves_wind(reg, onshorearea, pop, [0:step:25-step;], step)
	cap[:wind,area,pop] = (demandshare_wa, demandshare_wb, demandshare_off)
	lc[:wind,area,pop] = (lcoe_wa, lcoe_wb, lcoe_off)

	cap[:pv,:high,pop] = cap[:pv,:low,pop] .* 2
	lc[:pv,:high,pop] = lc[:pv,:low,pop]
	cap[:wind,:high,pop] = cap[:wind,:low,pop] .* 2
	lc[:wind,:high,pop] = lc[:wind,:low,pop]

	return cap, lc
end

function plotwindturbinecurve()
	speeds = powercurves[:,1]
	cfturbine = powercurves[:,2]
	cfpark = powercurves[:,3]
	display(plot(speeds, cfturbine, line=(3, :dash, RGB([149,179,215]/255...)), label="single turbine"))
	display(plot!(speeds, cfpark, line=(3, :solid, RGB([149,179,215]/255...)), label="park output",
					size=(1000,500), xlim=(0,32), ylim=(0,1.1), title="Wind turbine and park output",
					xlabel="Wind speed (m/s)", ylabel="Instantaneous capacity factor",
					tickfont=18, legendfont=18, guidefont=18, titlefont=20))
end

# Plot supply curves for a given model region
function plotsupplycurves(reg)
	plotly()

	cap, lc = calcvars(reg)

	plottriplepv(cap, lc, :low, 150)
	plottriplewind(cap, lc, :low, 150)
	plotmono(cap, lc)					# this is figure 3 in the model paper
end

# plotsupplycurves("China6")
# plotsupplycurves("Europe8")

const powercurves = [
         0         0         0
    0.1000         0    0.0000
    0.2000         0    0.0000
    0.3000         0    0.0000
    0.4000         0    0.0000
    0.5000         0    0.0001
    0.6000         0    0.0001
    0.7000         0    0.0001
    0.8000         0    0.0002
    0.9000         0    0.0003
    1.0000         0    0.0003
    1.1000         0    0.0005
    1.2000         0    0.0006
    1.3000         0    0.0008
    1.4000         0    0.0010
    1.5000         0    0.0013
    1.6000         0    0.0016
    1.7000         0    0.0020
    1.8000         0    0.0025
    1.9000         0    0.0030
    2.0000         0    0.0036
    2.1000         0    0.0044
    2.2000         0    0.0052
    2.3000         0    0.0062
    2.4000         0    0.0073
    2.5000         0    0.0085
    2.6000         0    0.0098
    2.7000         0    0.0114
    2.8000         0    0.0130
    2.9000         0    0.0148
    3.0000         0    0.0168
    3.1000    0.0106    0.0190
    3.2000    0.0147    0.0213
    3.3000    0.0187    0.0237
    3.4000    0.0229    0.0263
    3.5000    0.0270    0.0291
    3.6000    0.0312    0.0321
    3.7000    0.0354    0.0352
    3.8000    0.0397    0.0385
    3.9000    0.0440    0.0419
    4.0000    0.0484    0.0455
    4.1000    0.0529    0.0492
    4.2000    0.0574    0.0531
    4.3000    0.0621    0.0572
    4.4000    0.0668    0.0615
    4.5000    0.0716    0.0659
    4.6000    0.0766    0.0705
    4.7000    0.0816    0.0753
    4.8000    0.0868    0.0802
    4.9000    0.0921    0.0854
    5.0000    0.0976    0.0908
    5.1000    0.1032    0.0964
    5.2000    0.1090    0.1023
    5.3000    0.1150    0.1084
    5.4000    0.1214    0.1147
    5.5000    0.1281    0.1213
    5.6000    0.1353    0.1282
    5.7000    0.1429    0.1354
    5.8000    0.1510    0.1428
    5.9000    0.1597    0.1506
    6.0000    0.1691    0.1587
    6.1000    0.1791    0.1670
    6.2000    0.1898    0.1757
    6.3000    0.2011    0.1847
    6.4000    0.2129    0.1940
    6.5000    0.2253    0.2037
    6.6000    0.2380    0.2136
    6.7000    0.2512    0.2238
    6.8000    0.2647    0.2344
    6.9000    0.2786    0.2453
    7.0000    0.2927    0.2564
    7.1000    0.3070    0.2679
    7.2000    0.3216    0.2796
    7.3000    0.3363    0.2917
    7.4000    0.3514    0.3040
    7.5000    0.3666    0.3166
    7.6000    0.3822    0.3294
    7.7000    0.3980    0.3425
    7.8000    0.4141    0.3558
    7.9000    0.4305    0.3694
    8.0000    0.4473    0.3832
    8.1000    0.4643    0.3972
    8.2000    0.4817    0.4114
    8.3000    0.4994    0.4257
    8.4000    0.5175    0.4402
    8.5000    0.5360    0.4549
    8.6000    0.5548    0.4697
    8.7000    0.5740    0.4845
    8.8000    0.5937    0.4994
    8.9000    0.6137    0.5144
    9.0000    0.6341    0.5293
    9.1000    0.6550    0.5442
    9.2000    0.6762    0.5590
    9.3000    0.6975    0.5738
    9.4000    0.7188    0.5884
    9.5000    0.7400    0.6028
    9.6000    0.7609    0.6170
    9.7000    0.7814    0.6309
    9.8000    0.8014    0.6446
    9.9000    0.8206    0.6579
   10.0000    0.8390    0.6709
   10.1000    0.8565    0.6835
   10.2000    0.8729    0.6958
   10.3000    0.8884    0.7076
   10.4000    0.9028    0.7189
   10.5000    0.9161    0.7299
   10.6000    0.9283    0.7403
   10.7000    0.9394    0.7503
   10.8000    0.9494    0.7599
   10.9000    0.9582    0.7690
   11.0000    0.9659    0.7776
   11.1000    0.9723    0.7858
   11.2000    0.9777    0.7935
   11.3000    0.9821    0.8009
   11.4000    0.9856    0.8078
   11.5000    0.9884    0.8143
   11.6000    0.9906    0.8205
   11.7000    0.9922    0.8264
   11.8000    0.9934    0.8319
   11.9000    0.9944    0.8371
   12.0000    0.9951    0.8420
   12.1000    0.9958    0.8467
   12.2000    0.9964    0.8512
   12.3000    0.9970    0.8554
   12.4000    0.9975    0.8595
   12.5000    0.9980    0.8633
   12.6000    0.9984    0.8670
   12.7000    0.9988    0.8706
   12.8000    0.9992    0.8740
   12.9000    0.9996    0.8774
   13.0000    1.0000    0.8806
   13.1000    1.0000    0.8837
   13.2000    1.0000    0.8867
   13.3000    1.0000    0.8896
   13.4000    1.0000    0.8925
   13.5000    1.0000    0.8953
   13.6000    1.0000    0.8980
   13.7000    1.0000    0.9006
   13.8000    1.0000    0.9031
   13.9000    1.0000    0.9056
   14.0000    1.0000    0.9080
   14.1000    1.0000    0.9103
   14.2000    1.0000    0.9126
   14.3000    1.0000    0.9147
   14.4000    1.0000    0.9167
   14.5000    1.0000    0.9187
   14.6000    1.0000    0.9205
   14.7000    1.0000    0.9223
   14.8000    1.0000    0.9239
   14.9000    1.0000    0.9254
   15.0000    1.0000    0.9268
   15.1000    1.0000    0.9281
   15.2000    1.0000    0.9293
   15.3000    1.0000    0.9304
   15.4000    1.0000    0.9314
   15.5000    1.0000    0.9323
   15.6000    1.0000    0.9330
   15.7000    1.0000    0.9337
   15.8000    1.0000    0.9343
   15.9000    1.0000    0.9349
   16.0000    1.0000    0.9353
   16.1000    1.0000    0.9357
   16.2000    1.0000    0.9361
   16.3000    1.0000    0.9363
   16.4000    1.0000    0.9366
   16.5000    1.0000    0.9368
   16.6000    1.0000    0.9369
   16.7000    1.0000    0.9371
   16.8000    1.0000    0.9372
   16.9000    1.0000    0.9373
   17.0000    1.0000    0.9373
   17.1000    1.0000    0.9374
   17.2000    1.0000    0.9374
   17.3000    1.0000    0.9374
   17.4000    1.0000    0.9375
   17.5000    1.0000    0.9375
   17.6000    1.0000    0.9375
   17.7000    1.0000    0.9375
   17.8000    1.0000    0.9375
   17.9000    1.0000    0.9375
   18.0000    1.0000    0.9375
   18.1000    1.0000    0.9375
   18.2000    1.0000    0.9375
   18.3000    1.0000    0.9375
   18.4000    1.0000    0.9375
   18.5000    1.0000    0.9375
   18.6000    1.0000    0.9375
   18.7000    1.0000    0.9375
   18.8000    1.0000    0.9375
   18.9000    1.0000    0.9375
   19.0000    1.0000    0.9375
   19.1000    1.0000    0.9375
   19.2000    1.0000    0.9375
   19.3000    1.0000    0.9375
   19.4000    1.0000    0.9375
   19.5000    1.0000    0.9375
   19.6000    1.0000    0.9375
   19.7000    1.0000    0.9375
   19.8000    1.0000    0.9375
   19.9000    1.0000    0.9375
   20.0000    1.0000    0.9375
   20.1000    1.0000    0.9375
   20.2000    1.0000    0.9375
   20.3000    1.0000    0.9375
   20.4000    1.0000    0.9375
   20.5000    1.0000    0.9375
   20.6000    1.0000    0.9375
   20.7000    1.0000    0.9375
   20.8000    1.0000    0.9375
   20.9000    1.0000    0.9375
   21.0000    1.0000    0.9375
   21.1000    1.0000    0.9375
   21.2000    1.0000    0.9375
   21.3000    1.0000    0.9375
   21.4000    1.0000    0.9375
   21.5000    1.0000    0.9375
   21.6000    1.0000    0.9375
   21.7000    1.0000    0.9375
   21.8000    1.0000    0.9375
   21.9000    1.0000    0.9375
   22.0000    1.0000    0.9375
   22.1000    1.0000    0.9375
   22.2000    1.0000    0.9374
   22.3000    1.0000    0.9372
   22.4000    1.0000    0.9370
   22.5000    1.0000    0.9366
   22.6000    1.0000    0.9361
   22.7000    1.0000    0.9353
   22.8000    1.0000    0.9343
   22.9000    1.0000    0.9329
   23.0000    1.0000    0.9311
   23.1000    1.0000    0.9289
   23.2000    1.0000    0.9260
   23.3000    1.0000    0.9226
   23.4000    1.0000    0.9183
   23.5000    1.0000    0.9131
   23.6000    1.0000    0.9068
   23.7000    1.0000    0.8994
   23.8000    1.0000    0.8906
   23.9000    1.0000    0.8803
   24.0000    1.0000    0.8684
   24.1000    1.0000    0.8547
   24.2000    1.0000    0.8392
   24.3000    1.0000    0.8216
   24.4000    1.0000    0.8020
   24.5000    1.0000    0.7804
   24.6000    1.0000    0.7566
   24.7000    1.0000    0.7308
   24.8000    1.0000    0.7030
   24.9000    1.0000    0.6734
   25.0000    1.0000    0.6421
   25.1000    0.9000    0.6093
   25.2000    0.8000    0.5753
   25.3000    0.7000    0.5403
   25.4000    0.6000    0.5047
   25.5000    0.5000    0.4688
   25.6000    0.4000    0.4328
   25.7000    0.3000    0.3972
   25.8000    0.2000    0.3622
   25.9000    0.1000    0.3282
   26.0000         0    0.2954
   26.1000         0    0.2641
   26.2000         0    0.2345
   26.3000         0    0.2067
   26.4000         0    0.1809
   26.5000         0    0.1571
   26.6000         0    0.1355
   26.7000         0    0.1159
   26.8000         0    0.0983
   26.9000         0    0.0828
   27.0000         0    0.0691
   27.1000         0    0.0572
   27.2000         0    0.0469
   27.3000         0    0.0381
   27.4000         0    0.0307
   27.5000         0    0.0244
   27.6000         0    0.0192
   27.7000         0    0.0149
   27.8000         0    0.0115
   27.9000         0    0.0086
   28.0000         0    0.0064
   28.1000         0    0.0046
   28.2000         0    0.0032
   28.3000         0    0.0022
   28.4000         0    0.0014
   28.5000         0    0.0009
   28.6000         0    0.0005
   28.7000         0    0.0003
   28.8000         0    0.0001
   28.9000         0    0.0000
   29.0000         0         0
   29.1000         0         0
   29.2000         0         0
   29.3000         0         0
   29.4000         0         0
   29.5000         0         0
   29.6000         0         0
   29.7000         0         0
   29.8000         0         0
   29.9000         0         0
   30.0000         0         0
   30.1000         0         0
   30.2000         0         0
   30.3000         0         0
   30.4000         0         0
   30.5000         0         0
   30.6000         0         0
   30.7000         0         0
   30.8000         0         0
   30.9000         0         0
   31.0000         0         0
   31.1000         0         0
   31.2000         0         0
   31.3000         0         0
   31.4000         0         0
   31.5000         0         0
   31.6000         0         0
   31.7000         0         0
   31.8000         0         0
   31.9000         0         0
   32.0000         0         0
   32.1000         0         0
   32.2000         0         0
   32.3000         0         0
   32.4000         0         0
   32.5000         0         0
   32.6000         0         0
   32.7000         0         0
   32.8000         0         0
   32.9000         0         0
   33.0000         0         0
   33.1000         0         0
   33.2000         0         0
   33.3000         0         0
   33.4000         0         0
   33.5000         0         0
   33.6000         0         0
   33.7000         0         0
   33.8000         0         0
   33.9000         0         0
   34.0000         0         0
   34.1000         0         0
   34.2000         0         0
   34.3000         0         0
   34.4000         0         0
   34.5000         0         0
   34.6000         0         0
   34.7000         0         0
   34.8000         0         0
   34.9000         0         0
   35.0000         0         0
]

nothing


text:
#####
##### `AbstractDifferential`
#####

"""
The subtypes of `AbstractDifferential` define a custom \"algebra\" for chain
rule evaluation that attempts to factor various features like complex derivative
support, broadcast fusion, zero-elision, etc. into nicely separated parts.

All subtypes of `AbstractDifferential` implement the following operations:

`add(a, b)`: linearly combine differential `a` and differential `b`

`mul(a, b)`: multiply the differential `a` by the differential `b`

`Base.conj(x)`: complex conjugate of the differential `x`

`extern(x)`: convert `x` into an appropriate non-`AbstractDifferential` type for
use outside of `ChainContext`.

Valid arguments to these operations are `T` where `T<:AbstractDifferential`, or
where `T` has proper `+` and `*` implementations.

Additionally, all subtypes of `AbstractDifferential` support `Base.iterate` and
`Base.Broadcast.broadcastable(x)`.
"""
abstract type AbstractDifferential end

"""
    extern(x)

Return `x` converted to an appropriate non-`AbstractDifferential` type, for use
with external packages that might not handle `AbstractDifferential` types.

Note that this function may return an alias (not necessarily a copy) to data
wrapped by `x`, such that mutating `extern(x)` might mutate `x` itself.
"""
@inline extern(x) = x

@inline Base.conj(x::AbstractDifferential) = x

#=
This `AbstractDifferential` algebra has a monad-y "fallthrough" implementation;
each step handles an element of the algebra before dispatching to the next step.
This way, we don't need to implement promotion/conversion rules between subtypes
of `AbstractDifferential` to resolve potential ambiguities.
=#

const PRECEDENCE_LIST = [:wirtinger, :casted, :zero, :dne, :one, :thunk, :fallback]

global defs = Expr(:block)

let previous_add_name = :add, previous_mul_name = :mul
    for name in PRECEDENCE_LIST
        next_add_name = Symbol(string(:add_, name))
        next_mul_name = Symbol(string(:mul_, name))
        push!(defs.args, quote
            @inline $(previous_add_name)(a, b) = $(next_add_name)(a, b)
            @inline $(previous_mul_name)(a, b) = $(next_mul_name)(a, b)
        end)
        previous_add_name = next_add_name
        previous_mul_name = next_mul_name
    end
end

eval(defs)

@inline add_fallback(a, b) = a + b

@inline mul_fallback(a, b) = a * b

@inline add(x) = x

@inline mul(x) = x

#####
##### `Wirtinger`
#####

"""
    Wirtinger(primal::Union{Number,AbstractDifferential},
              conjugate::Union{Number,AbstractDifferential})

Returns a `Wirtinger` instance representing the complex differential:

```
df = ∂f/∂z * dz + ∂f/∂z̄ * dz̄
```

where `primal` corresponds to `∂f/∂z * dz` and `conjugate` corresponds to `∂f/∂z̄ * dz̄`.

The two fields of the returned instance can be accessed generically via the
[`wirtinger_primal`](@ref) and [`wirtinger_conjugate`](@ref) methods.
"""
struct Wirtinger{P,C} <: AbstractDifferential
    primal::P
    conjugate::C
    function Wirtinger(primal::Union{Number,AbstractDifferential},
                       conjugate::Union{Number,AbstractDifferential})
        return new{typeof(primal),typeof(conjugate)}(primal, conjugate)
    end
    function Wirtinger(primal, conjugate)
        error("`Wirtinger` only supports elements of type <: Union{Number,AbstractDifferential} for now")
    end
end

wirtinger_primal(x::Wirtinger) = x.primal
wirtinger_primal(x) = x

wirtinger_conjugate(x::Wirtinger) = x.conjugate
wirtinger_conjugate(::Any) = Zero()

extern(x::Wirtinger) = error("`Wirtinger` cannot be converted into an external type.")

Base.Broadcast.broadcastable(w::Wirtinger) = Wirtinger(broadcastable(w.primal),
                                                       broadcastable(w.conjugate))

Base.iterate(x::Wirtinger) = (x, nothing)
Base.iterate(::Wirtinger, ::Any) = nothing

Base.conj(x::Wirtinger) = error("`conj(::Wirtinger)` not yet defined")

function add_wirtinger(a::Wirtinger, b::Wirtinger)
    return Wirtinger(add(a.primal, b.primal), add(a.conjugate, b.conjugate))
end

add_wirtinger(a::Wirtinger, b) = add(a, Wirtinger(b, Zero()))
add_wirtinger(a, b::Wirtinger) = add(Wirtinger(a, Zero()), b)

function mul_wirtinger(a::Wirtinger, b::Wirtinger)
    error("""
          cannot multiply two Wirtinger objects; this error likely means a
          `WirtingerRule` was inappropriately defined somewhere. Multiplication
          of two Wirtinger objects is not defined because chain rule application
          often expands into a non-commutative operation in the Wirtinger
          calculus. To put it another way: simply given two Wirtinger objects
          and no other information, we can't know "locally" which components to
          conjugate in order to implement the chain rule. We could pick a
          convention; for example, we could define `a::Wirtinger * b::Wirtinger`
          such that we assume the chain rule application is of the form `f_a ∘ f_b`
          instead of `f_b ∘ f_a`. However, picking such a convention is likely to
          lead to silently incorrect derivatives due to commutativity assumptions
          in downstream generic code that deals with the reals. Thus, ChainRules
          makes this operation an error instead.
          """)
end

mul_wirtinger(a::Wirtinger, b) = Wirtinger(mul(a.primal, b), mul(a.conjugate, b))
mul_wirtinger(a, b::Wirtinger) = Wirtinger(mul(a, b.primal), mul(a, b.conjugate))

#####
##### `Casted`
#####

"""
TODO
"""
struct Casted{V} <: AbstractDifferential
    value::V
end

cast(x) = Casted(x)
cast(f, args...) = Casted(broadcasted(f, args...))

extern(x::Casted) = materialize(broadcasted(extern, x.value))

Base.Broadcast.broadcastable(x::Casted) = x.value

Base.iterate(x::Casted) = iterate(x.value)
Base.iterate(x::Casted, state) = iterate(x.value, state)

Base.conj(x::Casted) = cast(conj, x.value)

add_casted(a::Casted, b::Casted) = Casted(broadcasted(add, a.value, b.value))
add_casted(a::Casted, b) = Casted(broadcasted(add, a.value, b))
add_casted(a, b::Casted) = Casted(broadcasted(add, a, b.value))

mul_casted(a::Casted, b::Casted) = Casted(broadcasted(mul, a.value, b.value))
mul_casted(a::Casted, b) = Casted(broadcasted(mul, a.value, b))
mul_casted(a, b::Casted) = Casted(broadcasted(mul, a, b.value))

#####
##### `Zero`
#####

"""
TODO
"""
struct Zero <: AbstractDifferential end

extern(x::Zero) = false

Base.Broadcast.broadcastable(::Zero) = Ref(Zero())

Base.iterate(x::Zero) = (x, nothing)
Base.iterate(::Zero, ::Any) = nothing

add_zero(::Zero, ::Zero) = Zero()
add_zero(::Zero, b) = b
add_zero(a, ::Zero) = a

mul_zero(::Zero, ::Zero) = Zero()
mul_zero(::Zero, ::Any) = Zero()
mul_zero(::Any, ::Zero) = Zero()

#####
##### `DNE`
#####

"""
TODO
"""
struct DNE <: AbstractDifferential end

extern(x::DNE) = error("`DNE` cannot be converted into an external type.")

Base.Broadcast.broadcastable(::DNE) = Ref(DNE())

Base.iterate(x::DNE) = (x, nothing)
Base.iterate(::DNE, ::Any) = nothing

add_dne(::DNE, ::DNE) = DNE()
add_dne(::DNE, b) = b
add_dne(a, ::DNE) = a

mul_dne(::DNE, ::DNE) = DNE()
mul_dne(::DNE, ::Any) = DNE()
mul_dne(::Any, ::DNE) = DNE()

#####
##### `One`
#####

"""
TODO
"""
struct One <: AbstractDifferential end

extern(x::One) = true

Base.Broadcast.broadcastable(::One) = Ref(One())

Base.iterate(x::One) = (x, nothing)
Base.iterate(::One, ::Any) = nothing

add_one(a::One, b::One) = add(extern(a), extern(b))
add_one(a::One, b) = add(extern(a), b)
add_one(a, b::One) = add(a, extern(b))

mul_one(::One, ::One) = One()
mul_one(::One, b) = b
mul_one(a, ::One) = a

#####
##### `Thunk`
#####

"""
TODO
"""
struct Thunk{F} <: AbstractDifferential
    f::F
end

macro thunk(body)
    return :(Thunk(() -> $(esc(body))))
end

@inline extern(x::Thunk{F}) where {F} = x.f()

Base.Broadcast.broadcastable(x::Thunk) = broadcastable(extern(x))

@inline function Base.iterate(x::Thunk)
    externed = extern(x)
    element, state = iterate(externed)
    return element, (externed, state)
end

@inline function Base.iterate(::Thunk, (externed, state))
    element, new_state = iterate(externed, state)
    return element, (externed, new_state)
end

Base.conj(x::Thunk) = @thunk(conj(extern(x)))

add_thunk(a::Thunk, b::Thunk) = add(extern(a), extern(b))
add_thunk(a::Thunk, b) = add(extern(a), b)
add_thunk(a, b::Thunk) = add(a, extern(b))

mul_thunk(a::Thunk, b::Thunk) = mul(extern(a), extern(b))
mul_thunk(a::Thunk, b) = mul(extern(a), b)
mul_thunk(a, b::Thunk) = mul(a, extern(b))

#####
##### misc.
#####

"""
    Wirtinger(primal::Real, conjugate::Real)

Return `add(primal, conjugate)`.

Actually implementing the Wirtinger calculus generally requires that the
summed terms of the Wirtinger differential (`∂f/∂z * dz` and `∂f/∂z̄ * dz̄`) be
stored individually. However, if both of these terms are real-valued, then
downstream Wirtinger propagation mechanisms resolve to the same mechanisms as
real-valued calculus, so that the terms' sum can be eagerly computed and
propagated without requiring a special `Wirtinger` representation

This method primarily exists as an optimization.
"""
function Wirtinger(primal::Union{Real,DNE,Zero,One},
                   conjugate::Union{Real,DNE,Zero,One})
    return add(primal, conjugate)
end


text:
#!/usr/bin/env julia
######################################################################
# Overwrite REQUIRE using dependency information from Project.toml.
#
# Call from the root of the package repository.
#
# Has some basic sanity checks, but **use at your own risk**, `REQUIRE`
# will be overwritten.
#
# The purpose of this script is to appease attobot, until
# https://github.com/attobot/attobot/issues/50 is fixed.
######################################################################

@assert VERSION ≥ v"0.7"

import Pkg
const PT = Pkg.Types

Pkg.activate(pwd())             # current directory as the project
ctx = PT.Context()
pkg = ctx.env.pkg
if pkg ≡ nothing
    @error "Not in a package, I won't generate REQUIRE."
    exit(1)
else
    @info "found package" pkg = pkg
end

deps = PT.get_deps(ctx)
non_std_deps = sort(collect(setdiff(keys(deps), values(ctx.stdlibs))))

open("REQUIRE", "w") do io
    println(io, "julia 0.7")
    for d in non_std_deps
        println(io, d)
        @info "listing $d"
    end
end


text:
using CxxWrap: ConstCxxRef
using CxxWrap.StdLib: StdString, SharedPtr, UniquePtr, WeakPtr
using CxxWrap.CxxWrapCore:
    CxxBool,
    CxxChar,
    CxxWchar,
    CxxSignedChar,
    CxxUChar,
    CxxShort,
    CxxUShort,
    CxxInt,
    CxxUInt,
    CxxLong,
    CxxULong,
    CxxLongLong,
    CxxULongLong,
    SmartPointer

function Base.getproperty(x::Union{Contract, SmartPointer{Contract}}, name::Symbol)
    Base.@_inline_meta
    return if name === :conId
            Contract_get_conId(x)
        elseif name === :comboLegs
            Contract_get_comboLegs(x)
        elseif name === :comboLegsDescrip
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_comboLegsDescrip(x))
        elseif name === :secId
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_secId(x))
        elseif name === :secIdType
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_secIdType(x))
        elseif name === :includeExpired
            Contract_get_includeExpired(x)
        elseif name === :tradingClass
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_tradingClass(x))
        elseif name === :localSymbol
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_localSymbol(x))
        elseif name === :currency
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_currency(x))
        elseif name === :primaryExchange
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_primaryExchange(x))
        elseif name === :exchange
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_exchange(x))
        elseif name === :multiplier
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_multiplier(x))
        elseif name === :right
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_right(x))
        elseif name === :strike
            Contract_get_strike(x)
        elseif name === :lastTradeDateOrContractMonth
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_lastTradeDateOrContractMonth(x))
        elseif name === :secType
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_secType(x))
        elseif name === :symbol
            CxxWrap.CxxWrapCore.dereference_argument(Contract_get_symbol(x))
        else
            (Base).getfield(x, name)
        end
end
function Base.setproperty!(x::Union{Contract, SmartPointer{Contract}}, name::Symbol, value)
    Base.@_inline_meta
    if name === :conId
        Contract_set_conId(x, value)
    elseif name === :comboLegsDescrip
        Contract_set_comboLegsDescrip(x, value)
    elseif name === :secId
        Contract_set_secId(x, value)
    elseif name === :secIdType
        Contract_set_secIdType(x, value)
    elseif name === :includeExpired
        Contract_set_includeExpired(x, value)
    elseif name === :tradingClass
        Contract_set_tradingClass(x, value)
    elseif name === :localSymbol
        Contract_set_localSymbol(x, value)
    elseif name === :currency
        Contract_set_currency(x, value)
    elseif name === :primaryExchange
        Contract_set_primaryExchange(x, value)
    elseif name === :exchange
        Contract_set_exchange(x, value)
    elseif name === :multiplier
        Contract_set_multiplier(x, value)
    elseif name === :right
        Contract_set_right(x, value)
    elseif name === :strike
        Contract_set_strike(x, value)
    elseif name === :lastTradeDateOrContractMonth
        Contract_set_lastTradeDateOrContractMonth(x, value)
    elseif name === :secType
        Contract_set_secType(x, value)
    elseif name === :symbol
        Contract_set_symbol(x, value)
    else
        (Base).setfield!(x, name, value)
    end
    return x
end
function Base.propertynames(x::Union{Contract, SmartPointer{Contract}})
    Base.@_inline_meta
    return (:conId, :symbol, :secType, :lastTradeDateOrContractMonth, :strike, :right, :multiplier, :exchange, :primaryExchange, :currency, :localSymbol, :tradingClass, :includeExpired, :secIdType, :secId, :comboLegsDescrip, :comboLegs)
end


text:
include("_load_common.jl")
include("lp.jl")
include("lp_jad.jl")
include("lp_bench.jl")
include("lp_glpk.jl")

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
using CairoMakie, ProgressMeter, Statistics

"""
    plot_rhoUV(
        ρ,u,v,
        Δx,Δy,
        range,
        save_plot=true,save_video=false,dir_name=pwd(),file_name="0.png";
        theme="dark", resolution=(800,600),arrow_step=20,fps=40,kwargs...
    )

Return the figure, axes, heatmap and arrows associated to the density `ρ` and velocity field
`(u,v)`. The plot is saved in the directory `dir_name` in the file `file_name`. Two themes
are available, either `"dark"` (default) or `"light"`.
"""
function plot_rhoUV(
    ρ,u,v,
    Δx,Δy,
    range,
    save_plot=true,save_video=false,dir_name=pwd(),file_name="0.png";
    theme="dark", resolution=(800,600),arrow_step=20,fps=40,kwargs...
)
    if theme == "dark"
        set_theme!(theme_black())
        colormap = :linear_kryw_5_100_c67_n256
        arrowcolor = :black
    elseif theme == "light"
        set_theme!()
        colormap = Reverse(:linear_kryw_5_100_c67_n256)
        arrowcolor = :white
    else
        AttributeError("Theme not defined!")
    end
    ncellx = size(ρ)[1] - 2
    ncelly = size(ρ)[2] - 2
    Lx = ncellx*Δx
    Ly = ncelly*Δy
    x = Δx .* collect(0:(ncellx-1)) .+ Δx/2
    y = Δy .* collect(0:(ncelly-1)) .+ Δy/2
    rho = ρ[2:ncellx+1,2:ncelly+1]
    fig = Figure(resolution=resolution)
    ax = Axis(fig[1, 1], title="time=0.00")
    hm = heatmap!(ax,x,y,rho,colorrange=range,colormap=colormap)
    ax.aspect = AxisAspect(1)
    xlims!(ax,0,Lx)
    ylims!(ax,0,Ly)
    Colorbar(fig[1,2],hm)

    extract_x = floor(Int,arrow_step/2):arrow_step:ncellx
    extract_y = floor(Int,arrow_step/2):arrow_step:ncelly
    xx = x[extract_x]
    yy = y[extract_y]
    u_grid = u[2:ncellx+1,2:ncelly+1]
    v_grid = v[2:ncellx+1,2:ncelly+1]
    uu = u_grid[extract_x,extract_y]
    vv = v_grid[extract_x,extract_y]
    arrows = arrows!(ax,xx,yy,uu,vv)
    arrows.color = arrowcolor
    arrows.lengthscale = 2/3*arrow_step*Δx
    arrows.arrowsize = floor(Int,1/80*resolution[1])
    arrows.linewidth = floor(Int,3/800*resolution[1])
    arrows.origin = :center
    if save_plot
        save(joinpath(dir_name,file_name),fig)
    end
    if save_video
        stream = VideoStream(fig,framerate=fps)
        return fig, ax, hm, arrows, stream
    else
        return fig, ax, hm, arrows
    end
end

"""
    update_plot!(fig,ax,hm,arrows,ρ,u,v,title,dir_name,file_name,save_plot=true,stream=nothing)

Update a plot created with the function `plot_rhoUV` and save it.
"""
function update_plot!(fig,ax,hm,arrows,ρ,u,v,title,dir_name,file_name,save_plot=true,stream=nothing)
    ncellx = size(ρ)[1] - 2
    ncelly = size(ρ)[2] - 2
    rho = ρ[2:ncellx+1,2:ncelly+1]
    hm[3] = rho
    u_grid = u[2:ncellx+1,2:ncelly+1]
    v_grid = v[2:ncellx+1,2:ncelly+1]
    extract_x = 10:20:ncellx
    extract_y = 10:20:ncelly
    uu = u_grid[extract_x,extract_y]
    vv = v_grid[extract_x,extract_y]
    arrows[:directions] = vec(Vec2f.(uu,vv))
    ax.title = title
    if save_plot
        save(joinpath(dir_name,file_name),fig)
    end
    if !isnothing(stream)
        recordframe!(stream)
    end
end

"""
    save_data!(
        data,
        ρ,u,v,
        dir_name,file_name;key="data")

Update the dictionary `data` with the values `(ρ,u,v)` and save it.
"""
function save_data!(
    data,
    ρ,u,v,
    dir_name,file_name;key="data")
    data[:ρ] = ρ
    data[:u] = u
    data[:v] = v
    save(joinpath(dir_name,file_name),key,data)
end

"""
    radial_density(ρ,Δx,Δy,K=400)

Compute the radial density as the mean of `ρ` in `K` annuli of
equal lengths and evenly spread radii between 0 and the maximal
radius in the box. Return the vector of radii and the radial density.
"""
function radial_density(ρ,Δx,Δy,K=400)
    ncellx = size(ρ)[1] - 2
    ncelly = size(ρ)[2] - 2
    Lx = ncellx*Δx
    Ly = ncelly*Δy
    mx = Lx/2
    my = Ly/2
    rmat = zeros(ncellx,ncelly)
    for i in 1:ncellx
        for j in 1:ncellx
            xij = (i-1)*Δx + Δx/2
            yij = (j-1)*Δy + Δy/2
            rmat[i,j] = sqrt((xij-mx)^2 + (yij-my)^2)
        end
    end
    rvec = vec(rmat)
    rho = vec(ρ[2:ncellx+1,2:ncelly+1])
    r = LinRange(0.,maximum(rvec),K)
    dr = r[2]-r[1]
    rho_r = zeros(K)
    for k in 1:K
        rk = r[k]+dr
        rho_r[k] = mean(rho[abs.(rvec.-rk).<dr/2])
    end
    return r, rho_r
end


text:
module NeuralOperators
    using DataDeps
    using Fetch
    using MAT
    using StatsBase

    using Flux
    using FFTW
    using Tullio
    using CUDA
    using CUDAKernels
    using KernelAbstractions
    using Zygote
    using ChainRulesCore

    function __init__()
        register_datasets()
    end

    include("data.jl")
    include("fourier.jl")
    include("model.jl")
end


text:
@testset "zygote_adjoints" begin

    rng = MersenneTwister(123456)
    x = rand(rng, 5)
    y = rand(rng, 5)
    r = rand(rng, 5)

    gzeucl = gradient(:Zygote, [x,y]) do xy
        evaluate(Euclidean(), xy[1], xy[2])
    end
    gzsqeucl = gradient(:Zygote, [x,y]) do xy
        evaluate(SqEuclidean(), xy[1], xy[2])
    end
    gzdotprod = gradient(:Zygote, [x,y]) do xy
        evaluate(KernelFunctions.DotProduct(), xy[1], xy[2])
    end
    gzdelta = gradient(:Zygote, [x,y]) do xy
        evaluate(KernelFunctions.Delta(), xy[1], xy[2])
    end
    gzsinus = gradient(:Zygote, [x,y]) do xy
        evaluate(KernelFunctions.Sinus(r), xy[1], xy[2])
    end

    gfeucl = gradient(:FiniteDiff, [x,y]) do xy
        evaluate(Euclidean(), xy[1], xy[2])
    end
    gfsqeucl = gradient(:FiniteDiff, [x,y]) do xy
        evaluate(SqEuclidean(), xy[1], xy[2])
    end
    gfdotprod = gradient(:FiniteDiff, [x,y]) do xy
        evaluate(KernelFunctions.DotProduct(), xy[1], xy[2])
    end
    gfdelta = gradient(:FiniteDiff, [x,y]) do xy
        evaluate(KernelFunctions.Delta(), xy[1], xy[2])
    end
    gfsinus = gradient(:FiniteDiff, [x,y]) do xy
        evaluate(KernelFunctions.Sinus(r), xy[1], xy[2])
    end

    @test all(gzeucl .≈ gfeucl)
    @test all(gzsqeucl .≈ gfsqeucl)
    @test all(gzdotprod .≈ gfdotprod)
    @test all(gzdelta .≈ gfdelta)
    @test all(gzsinus .≈ gfsinus)
end


text:
# translate pyfpgrowth
using Combinatorics
using DataStructures
using DataFrames

const ROOT_STRING = "{}"                      # FIXME: hacky and fragile

type FPNode
    value::String
    count::Int
    parent::FPNode
    children::Array{FPNode,1}
    link

    function FPNode(value, count, parent)
        kids = Array{FPNode,1}(0)
        return new(value, count, parent, kids)
    end

    function FPNode(value, count)
        # This method is for root node initialization
        x = new(value, count)
        x.children = Array{FPNode,1}(0)
        x.parent = x
        return x
    end
end

type FPTree
    frequent::Dict{String,Int}
    headers::Dict{String,FPNode}
    root::FPNode

    # This methods is dispatched for the initial
    # creation of the FP tree.
    function FPTree(transactions, threshold)
        frequent = find_frequent_items(transactions, threshold)
        headers = Dict{String,FPNode}() #build_header_table(frequent)
        res = new(frequent, headers)
        res.root = build_fptree!(res,
                                 transactions,
                                 ROOT_STRING,
                                 0,
                                 frequent,
                                 headers)
        return res
    end

    function FPTree(tree, transactions, threshold, root_value, root_count)
        frequent = find_frequent_items(transactions, threshold)
        headers = Dict{String,FPNode}() #build_header_table(frequent)
        root = build_fptree!(tree,
                             transactions,
                             root_value,
                             root_count,
                             frequent,
                             headers)

        res = new(frequent, headers, root)
        return res
    end
end

function child_position(pnode::FPNode, value)
    res = 0                 # HACK: must be a better solution
    for i, child in enumerate(pnode.children)
        if child.value == value
            res = i
            break
        end
    end
    res
end

function increment_child!(pnode::FPNode, idx)
    pnode.children[idx].count += 1
end

function add_child!(pnode::FPNode, value)
    child = FPNode(value, 1, pnode)
    push!(pnode.children, child)
end

# Create Dict of items with occurrences above threshold.
function find_frequent_items(transactions, threshold)
    items = Dict{String,Int}()

    for transaction in transactions
        for item in transaction
            items[item] = get(items, item, 0) + 1
        end
    end

    for k in keys(items)
        if items[k] < threshold
            delete!(items, k)
        end
    end
    items
end

# function build_header_table(frequent)
#     headers = Dict()
#     for k in keys(frequent)
#         headers[k] = nothing
#     end
#     headers
# end

# This method is to be dispatched on creation of the actual
# root of the tree (i.e., the first node, which has no parent).
function build_fptree!(tree::FPTree, transactions, root_value, root_count, frequent, headers)
    root = FPNode(root_value, root_count)

    for transaction in transactions
        sorted_items = Array{String,1}(0)
        for x in transaction
            if x ∈ keys(frequent)
                push!(sorted_items, x)
            end
        end
        sort!(sorted_items, by = x -> frequent[x], rev = true)

        if length(sorted_items) > 0
            insert_tree!(tree, sorted_items, root, headers)
        end
    end
    root
end

function insert_tree!(tree::FPTree, items, node, headers)
    # recursively grow FP tree
    first_item = first(items)

    # child_idx will be 0 if not found, otherwise
    # it is the index of their position.
    child_idx = child_position(node, first_item)

    if child_idx ≠ 0
        increment_child!(node, child_idx)            # increment existing child counter
        child = node.children[child_idx]
    elseif child_idx == 0
        add_child!(node, first_item)                 # add new child
        child = last(node.children)                  # child we just inserted

        # Now we link the newly created node
        # to the header structure.
        if !haskey(headers, first_item)
            headers[first_item] = child
        else
            current = headers[first_item]
            while isdefined(current, :link)
                current = current.link
            end
            current.link = child
        end
    end

    # call function recursively
    remaining_items = items[2:end]
    if length(remaining_items) > 0
        insert_tree!(tree, remaining_items, child, headers)
    end
end

# This function returns a boolean indiecating whether or
# not this tree has a single path.
function tree_has_single_path(tree::FPTree, node)
    num_children = length(node.children)
    if num_children > 1
        return false
    elseif num_children == 0
        return true
    else
        return true && tree_has_single_path(tree, node.children[1])
    end
end

# Mine the constructed FP tree for patterns
function mine_patterns(tree, threshold)
    if tree_has_single_path(tree, tree.root)
        return generate_pattern_list(tree)
    else
        subtrees = mine_sub_trees(tree, threshold)
        return zip_patterns(tree, subtrees)
    end
end

# Append suffix to patterns in dictionary if
# we are in a conditional FP tree
function zip_patterns(tree, patterns)
    if isdefined(tree.root, :value)
        suffix = tree.root.value

        if suffix ≠ ROOT_STRING
            new_patterns = Dict{Array{String,1},Int}()

            for k in keys(patterns)
                xpattern = vcat(k, suffix)
                sort!(xpattern)
                new_patterns[xpattern] = patterns[k]
            end
            return new_patterns
        end
    end
    patterns
end

# Generate a list of patterns with support counts
function generate_pattern_list(tree)
    patterns = Dict{Array{String,1},Int}()
    items = collect(keys(tree.frequent))

    # If we are in a conditional tree,
    # the suffix is a pattern on its own.
    if !isdefined(tree.root, :value) || tree.root.value == ROOT_STRING
        suffix_value = Array{String,1}(0)
    else
        suffix_value = [tree.root.value]
        patterns[suffix_value] = tree.root.count
    end

    for i in 1:length(items)
        for subset in collect(combinations(items, i))
            pattern = sort(vcat(subset, suffix_value))
            patterns[pattern] = minimum([tree.frequent[x] for x in subset])
        end
    end
    return patterns
end

function mine_sub_trees(tree, threshold)
    patterns = Dict{Array{String,1},Int}()
    xkeys = collect(keys(tree.frequent))
    mining_order = sort(xkeys, by = x -> tree.frequent[x])

    # Get items in tree in reverse order of occurrences.
    for item in mining_order
        suffixes = Array{FPNode,1}(0)
        conditional_tree_input = Array{Array{String,1},1}(0)
        node = tree.headers[item]

        persist = true
        while persist
            push!(suffixes, node)
            if isdefined(node, :link)
                node = node.link
            else
                # push!(suffixes, node)
                persist = false
            end
        end

        # Follow node links to get a list of
        # all occurences of a certain item
        for suffix in suffixes
            frequency = suffix.count
            path = Array{String,1}(0)
            parent = suffix.parent
            while parent ≠ parent.parent
                push!(path, parent.value)
                parent = parent.parent
            end
            for i = 1:frequency
                push!(conditional_tree_input, path)
            end
        end
        # Now we have the input for a subtree,
        # so construct it and grab the patterns.
        subtree = FPTree(tree,
                         conditional_tree_input,
                         threshold,
                         item,
                         tree.frequent[item])

        subtree_patterns = mine_patterns(subtree, threshold)

        # Insert subtree patterns into main patterns dictionary.
        for pattern in keys(subtree_patterns)
            patterns[pattern] = get(patterns, pattern, 0) + subtree_patterns[pattern]
        end

    end
    return patterns
end

# Given a set of transactions, find the patterns in it
# over the specified support threshold.
function find_frequent_patterns(transactions, support_threshold)
    tree = FPTree(transactions, support_threshold)
    return mine_patterns(tree, support_threshold)
end

function make_transactions(X::Array{Any,2})
    n, p = size(X)
    X = map(string, X)

    out = Array{Array{String, 1}, 1}(n)
    for i = 1:n
        out[i] = convert(Array{String, 1}, X[i, :])
    end
    return out
end

function make_transactions(X::DataFrame)
    n, p = size(X)
    Xstr = Array{String}(n, p)
    for j = 1:p
        Xstr[:, j] = map(string, X[:, j])
    end

    out = Array{Array{String, 1}, 1}(n)
    for i = 1:n
        out[i] = convert(Array{String, 1}, Xstr[i, :])
    end
    return out
end


text:
"""
# module kernels

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-06-07

# Examples

```jldoctest
julia>
```
"""

module LaminartKernels

export kern_d, kern_b

# lgn to l6/l4
using NNlib, ImageFiltering, Images, OffsetArrays

function kern_d_ph(σ::Real, θ::Real, l = 4 * ceil(Int,σ) + 1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l>>1
    g = σ == 0 ? [1f0] : [exp(-x*cos(θ) / (2f0 * σ)) for x=-w:w]
    centered(g/sum(g))
end

function kern_d_pv(σ::Real, θ::Real, l = 4*ceil(Int,σ)+1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l>>1
    g = σ == 0 ? [1f0] : [exp(-x*sin(θ)/(2f0*σ)) for x=-w:w]
    centered(g/sum(g))
end

function kern_d_mh(σ::Real, θ::Real, l = 4*ceil(Int,σ)+1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l>>1
    g = σ == 0 ? [1f0] : [exp(x*cos(θ)/(2f0*σ)) for x=-w:w]
    centered(g/sum(g))
end

function kern_d_mv(σ::Real, θ::Real, l = 4*ceil(Int,σ)+1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l>>1
    g = σ == 0 ? [1f0] : [exp(x*sin(θ)/(2f0*σ)) for x=-w:w]
    centered(g/sum(g))
end

function kern_d_p(σ::Real, θ::Real, l = 4*ceil(Int,σ)+1)
    kern_d_pv(σ, θ, l) .* transpose(kern_d_ph(σ, θ, l))
end

function kern_d_m(σ::Real, θ::Real, l = 4*ceil(Int,σ)+1)
    kern_d_mv(σ, θ, l) .* transpose(kern_d_mh(σ, θ, l))
end

function kern_d(σ::Real, θ::Real, l = 4*ceil(Int,σ)+1)
    kern_d_p(σ, θ, l).-kern_d_m(σ, θ, l)
end

# function kern_d(σ::Real, θ::Real, l = 4*ceil(Int,σ)+1)
#     relu.(kern_d(σ, θ, l)) .- relu.(-1 .*(kern_d(σ, θ, l)))
# end

function kern_b(σ::Real, θ::Real, l = 4*ceil(Int,σ)+1)
    relu.(kern_d(σ, θ, l)) .+ relu.(-1f0 .*(kern_d(σ, θ, l)))
end

function gaussian_rot(σ_x::Real, σ_y::Real, θ::Real, l = 4*ceil(Int, max(σ_a,σ_b))+1)
    isodd(l) || throw(ArgumentError("length must be odd"))
    w = l>>1
    g = OffsetArray(fill(0.0f0, l, l), -w:w, -w:w)
    #todo add when σ_x or/and σ_y == 0
    for x ∈ -w:w, y ∈ -w:w
        g[x,y] = exp(-1f0/2f0*((((x*cos(θ)-y*sin(θ))/σ_x)^2f0)+(((x*sin(θ)+y*cos(θ))/σ_y)^2f0)))
    end
    centered(g/sum(g))
end

# Rotating function
function fun_R(x,y,θ)
    x*cos(θ)-y*sin(θ), x*sin(θ)+y*cos(θ)
end

end


text:
# export
#     cost

############################################################################################
#                              COST METHODS                                                #
############################################################################################

"""
	stage_cost(cost::CostFunction, z::AbstractKnotPoint)

Evaluate the cost at a knot point, and automatically handle terminal
knot point, multiplying by dt as necessary."""
function stage_cost(cost::CostFunction, z::AbstractKnotPoint)
    if is_terminal(z)
        stage_cost(cost, state(z))
    else
        stage_cost(cost, state(z), control(z))*z.dt
    end
end

"""
	cost(obj::Objective, Z::Traj)
	cost(obj::Objective, dyn_con::DynamicsConstraint{Q}, Z::Traj)

Evaluate the cost for a trajectory. If a dynamics constraint is given,
    use the appropriate integration rule, if defined.
"""
function cost(obj::AbstractObjective, Z::AbstractTrajectory{<:Any,<:Any,<:AbstractFloat})
    cost!(obj, Z)
    J = get_J(obj)
    return sum(J)
end

# ForwardDiff-able method
function cost(obj::AbstractObjective, Z::AbstractTrajectory{<:Any,<:Any,T}) where T
    J = zero(T)
    for k = 1:length(obj)
        J += stage_cost(obj[k], Z[k])
    end
    return J
end

# Default to no integration
cost(obj, dyn_con::DynamicsConstraint{Q}, Z) where Q<:QuadratureRule = cost(obj, Z)

"Evaluate the cost for a trajectory (non-allocating)"
@inline function cost!(obj::Objective, Z::AbstractTrajectory)
    map!(stage_cost, obj.J, obj.cost, Z)
end

############################################################################################
#                              COST EXPANSIONS                                             #
############################################################################################

"""
	cost_gradient!(E::Objective, obj::Objective, Z, init)

Evaluate the cost gradient along the entire tracjectory `Z`, storing the result in `E`.

If `init == true`, all gradients will be evaluated, even if they are constant.
"""
function cost_gradient!(E::Objective, obj::Objective, Z::AbstractTrajectory, init::Bool=false)
    is_const = E.const_grad
    N = length(Z)
    for k in eachindex(Z)
        if init || !is_const[k]
			if is_terminal(Z[k])
            	is_const[k] = gradient!(E.cost[k], obj.cost[k], state(Z[k]))
			else
            	is_const[k] = gradient!(E.cost[k], obj.cost[k], state(Z[k]), control(Z[k]))
			end
            dt_x = k < N ? Z[k].dt :  one(Z[k].dt)
            dt_u = k < N ? Z[k].dt : zero(Z[k].dt)
            E[k].q .*= dt_x
            E[k].r .*= dt_u
        end
    end
end

"""
	cost_hessian!(E::Objective, obj::Objective, Z, init)

Evaluate the cost hessian along the entire tracjectory `Z`, storing the result in `E`.

If `init == true`, all hessian will be evaluated, even if they are constant. If false,
they will only be evaluated if they are not constant.
"""
function cost_hessian!(E::Objective, obj::Objective, Z::AbstractTrajectory, init::Bool=false, rezero::Bool=false)
	is_const = E.const_hess
	if !init && all(is_const)
		return
	end
    N = length(Z)
    for k in eachindex(Z)
        if init || !is_const[k]
			if rezero
				E[k].Q .= 0
				E[k].R .= 0
				!is_blockdiag(E[k]) && (E[k].H .= 0)
			end
			if is_terminal(Z[k])
            	is_const[k] = hessian!(E.cost[k], obj.cost[k], state(Z[k]))
			else
				is_const[k] = hessian!(E.cost[k], obj.cost[k], state(Z[k]), control(Z[k]))
			end
            dt_x = k < N ? Z[k].dt :  one(Z[k].dt)
            dt_u = k < N ? Z[k].dt : zero(Z[k].dt)
            E[k].Q .*= dt_x
            E[k].R .*= dt_u
            E[k].H .*= dt_u
        end
    end
end

"""
	cost_expansion!(E::Objective, obj::Objective, Z, [init, rezero])

Evaluate the 2nd order Taylor expansion of the objective `obj` along the trajectory `Z`,
storing the result in `E`.

If `init == false`, the expansions will only be evaluated if they are not constant.

If `rezero == true`, all expansions will be multiplied by zero before taking the expansion.
"""
function cost_expansion!(E::Objective, obj::Objective, Z::Traj, init::Bool=false, rezero::Bool=false)
    cost_gradient!(E, obj, Z, init)
    cost_hessian!(E, obj, Z, init, rezero)
    return nothing
end

function error_expansion!(E::Objective, Jexp::Objective, model::AbstractModel, Z::Traj, G, tmp=G[end])
    @assert E === Jexp "E and Jexp should be the same object for AbstractModel"
    return nothing
end

function error_expansion!(E::Objective, Jexp::Objective, model::LieGroupModel, Z::Traj, G, tmp=G[end])
    for k in eachindex(E.cost)
        error_expansion!(E.cost[k], Jexp.cost[k], model, Z[k], G[k], tmp)
	end
	E.const_hess .= false   # hessian will always be dependent on the state
end

function error_expansion!(E::QuadraticCost, cost::QuadraticCost, model, z::AbstractKnotPoint,
        G, tmp)
	E.Q .= 0
    RobotDynamics.∇²differential!(E.Q, model, state(z), cost.q)
    if size(model)[1] < 15
        G = SMatrix(G)
        E.H .= SMatrix(cost.H) * G
        E.q .= G'cost.q
        E.Q .+= G'cost.Q*G
    else
        mul!(E.H, cost.H, G)
        mul!(E.q, Transpose(G), cost.q)
        mul!(tmp, cost.Q, G)
        mul!(E.Q, Transpose(G), tmp, 1.0, 1.0)
	end
    return nothing
end

struct StaticExpansion{T,N,M,NN,MM,NM} <: AbstractExpansion{T}
	x::SVector{N,T}
	xx::SMatrix{N,N,T,NN}
	u::SVector{M,T}
	uu::SMatrix{M,M,T,MM}
	ux::SMatrix{M,N,T,NM}
end

function StaticExpansion(E::AbstractExpansion)
	StaticExpansion(SVector(E.x), SMatrix(E.xx),
		SVector(E.u), SMatrix(E.uu), SMatrix(E.ux))
end

function StaticExpansion(x,xx,u,uu,ux)
	StaticExpansion(SVector(x), SMatrix(xx), SVector(u), SMatrix(uu), SMatrix(ux))
end

function static_expansion(cost::QuadraticCost)
	StaticExpansion(cost.q, cost.Q, cost.r, cost.R, cost.H)
end

"""
	dgrad(E::QuadraticExpansion, dZ::Traj)

Calculate the derivative of the cost in the direction of `dZ`, where `E` is the current
quadratic expansion of the cost.
"""
function dgrad(E::QuadraticExpansion{n,m,T}, dZ::Traj)::T where {n,m,T}
	g = zero(T)
	N = length(E)
	for k = 1:N-1
		g += dot(E[k].q, state(dZ[k])) + dot(E[k].r, control(dZ[k]))
	end
	g += dot(E[N].q, state(dZ[N]))
	return g
end

"""
	dhess(E::QuadraticCost, dZ::Traj)

Calculate the scalar 0.5*dZ'G*dZ where G is the hessian of cost
"""
function dhess(E::QuadraticExpansion{n,m,T}, dZ::Traj)::T where {n,m,T}
	h = zero(T)
	N = length(E)
	for k = 1:N-1
		x = state(dZ[k])
		u = control(dZ[k])
		h += dot(x, E[k].Q, x) + dot(u, E[k].R, u)
	end
	x = state(dZ[N])
	h += dot(x, E[N].Q, x)
	return 0.5*h
end

"""
	norm_grad(E::QuadraticExpansion, p=2)

Norm of the cost gradient
"""
function norm_grad(E::QuadraticExpansion{n,m,T}, p=2)::T where {n,m,T}
	J = get_J(E)
	for (k,cost) in enumerate(E)
		J[k] = norm(cost.q, p)
		if !cost.terminal
			J[k] = norm(SA[norm(cost.r, p) J[k]], p)
		end
	end
	return norm(J, p)
end


text:
# ---
# title: 195. Tenth Line
# id: problem195
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/tenth-line/description/>
# hidden: true
# ---
#
# Given a text file `file.txt`, print just the 10th line of the file.
#
# **Example:**
#
# Assume that `file.txt` has the following content:
#
#
#
#     Line 1
#     Line 2
#     Line 3
#     Line 4
#     Line 5
#     Line 6
#     Line 7
#     Line 8
#     Line 9
#     Line 10
#
#
# Your script should output the tenth line, which is:
#
#
#
#     Line 10
#
#
# **Note:**
# 1\. If the file contains less than 10 lines, what should you output?
# 2\. There's at least three different solutions. Try to explore all
# possibilities.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# ---
# title: 1521. Find a Value of a Mysterious Function Closest to Target
# id: problem1521
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Search, Bit Manipulation, Segment Tree
# link: <https://leetcode.com/problems/find-a-value-of-a-mysterious-function-closest-to-target/description/>
# hidden: true
# ---
#
# ![](https://assets.leetcode.com/uploads/2020/07/09/change.png)
#
# Winston was given the above mysterious function `func`. He has an integer
# array `arr` and an integer `target` and he wants to find the values `l` and
# `r` that make the value `|func(arr, l, r) - target|` minimum possible.
#
# Return _the minimum possible value_ of `|func(arr, l, r) - target|`.
#
# Notice that `func` should be called with the values `l` and `r` where `0 <= l,
# r < arr.length`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [9,12,3,7,15], target = 5
#     Output: 2
#     Explanation: Calling func with all the pairs of [l,r] = [[0,0],[1,1],[2,2],[3,3],[4,4],[0,1],[1,2],[2,3],[3,4],[0,2],[1,3],[2,4],[0,3],[1,4],[0,4]], Winston got the following results [9,12,3,7,15,8,0,3,7,0,0,3,0,0,0]. The value closest to 5 is 7 and 3, thus the minimum difference is 2.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1000000,1000000,1000000], target = 1
#     Output: 999999
#     Explanation: Winston called the func with all possible values of [l,r] and he always got 1000000, thus the min difference is 999999.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [1,2,4,8,16], target = 0
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^5`
#   * `1 <= arr[i] <= 10^6`
#   * `0 <= target <= 10^7`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function score(x, y)
    distance = (x^2 + y^2)^0.5
    if distance > 10
        return 0
    elseif distance > 5
        return 1
    elseif distance > 1
        return 5
    end

return 10

end

max_stars_repo_name: jfcliang/jl-exercism-solutions

max_issues_repo_name: jfcliang/jl-exercism-solutions

max_forks_repo_name: jfcliang/jl-exercism-solutions


text:
module PkgTest

include("ordmnfit_test.jl")
include("ordmntest_test.jl")
# include("../perf/polr_perf.jl")

end # PkgTest module


text:
import JSON, JSON2, PlotlyBase

JSON2.write(io::IO, p::PlotlyBase.Plot; kwargs...) = write(io, JSON.json(p))


text:
function add_export_callbacks(b::Gtk.GtkBuilder,handles::Tracker_Handles)

    signal_connect(export_button_cb,b["export_button"],"clicked",Void,(),false,(handles,))

    nothing
end

function export_button_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    setproperty!(han.b["export_button"],:label,"Exporting to MAT file...")

    e_angle = getproperty(han.b["angle_export_button"],:active,Bool)
    e_curve = getproperty(han.b["curvature_export_button"],:active,Bool)
    e_phase = getproperty(han.b["phase_export"],:active,Bool)
    e_whisker = getproperty(han.b["whisker_export"],:active,Bool)

    #Convert to Janelia
    (my_whiskers,mytracked)=convert_discrete_to_janelia(han.nn.predicted,han.nn.confidence_thres,han.wt.pad_pos);

    face_axis_num=getproperty(han.b["face_axis_combo"],:active,Int64)

    if face_axis_num == 0
        (mycurv,myangles)=get_curv_and_angle(my_whiskers,mytracked,han.wt.pad_pos);
    else
        (mycurv,myangles)=get_curv_and_angle(my_whiskers,mytracked,han.wt.pad_pos,face_axis='y');
    end

    #Whisking phase calculation
    phase_low_band=8.0 #Band pass filter lower band
    phase_high_band=30.0 #Band pass filter upper band
    myphase=WhiskerTracking.get_phase(myangles,bp_l=phase_low_band,bp_h=phase_high_band);

    filepath=string(han.wt.data_path,"output.mat")
    file=matopen(filepath,"w")
        if e_angle
            write(file,"Angles",myangles)
        end
        if e_curve
            write(file,"Curvature",mycurv)
        end
        if e_phase
            write(file,"Phase",myphase)
        end
        if e_whisker
            write(file,"Whisker",han.nn.predicted)
        end
    write(file,"Tracked",mytracked)
    close(file)

    setproperty!(han.b["export_button"],:label,"Export")
    println("Export complete")

    nothing
end


text:
"""
    Dataset(::PyObject)
    Dataset(; data_vars = nothing, coords = Dict(), attrs = Dict())

Loose wrapper around `xarray.Dataset`, mostly used for dispatch.

# Keywords

  - `data_vars::Dict{String,Any}`: Dict mapping variable names to

      + `Vector`: Data vector. Single dimension is named after variable.
      + `Tuple{String,Vector}`: Dimension name and data vector.
      + `Tuple{NTuple{N,String},Array{T,N}} where {N,T}`: Dimension names and data array.

  - `coords::Dict{String,Any}`: Dict mapping dimension names to index names. Possible
    arguments has same form as `data_vars`.
  - `attrs::Dict{String,Any}`: Global attributes to save on this dataset.

In most cases, use [`convert_to_dataset`](@ref) or [`convert_to_constant_dataset`](@ref) or
to create a `Dataset` instead of directly using a constructor.
"""
struct Dataset
    o::PyObject

    function Dataset(o::PyObject)
        if !pyisinstance(o, xarray.Dataset)
            throw(ArgumentError("$o is not an `xarray.Dataset`."))
        end
        return new(o)
    end
end

Dataset(; kwargs...) = Dataset(xarray.Dataset(; kwargs...))
@inline Dataset(data::Dataset) = data

@inline PyObject(data::Dataset) = getfield(data, :o)

Base.convert(::Type{Dataset}, obj::PyObject) = Dataset(obj)
Base.convert(::Type{Dataset}, obj::Dataset) = obj
Base.convert(::Type{Dataset}, obj) = convert_to_dataset(obj)

Base.hash(data::Dataset) = hash(PyObject(data))

Base.propertynames(data::Dataset) = propertynames(PyObject(data))

function Base.getproperty(data::Dataset, name::Symbol)
    o = PyObject(data)
    name === :o && return o
    return getproperty(o, name)
end

Base.getindex(data::Dataset, k) = py"$(data)[$k]"

function Base.show(io::IO, data::Dataset)
    out = pycall(pybuiltin("str"), String, data)
    out = replace(out, "<xarray.Dataset>" => "Dataset (xarray.Dataset)")
    print(io, out)
    return nothing
end
function Base.show(io::IO, ::MIME"text/html", data::Dataset)
    obj = PyObject(data)
    (:_repr_html_ in propertynames(obj)) || return show(io, data)
    out = obj._repr_html_()
    out = replace(out, r"(<|<)?xarray.Dataset(>|>)?" => "Dataset (xarray.Dataset)")
    print(io, out)
    return nothing
end

attributes(data::Dataset) = getproperty(PyObject(data), :_attrs)

function setattribute!(data::Dataset, key, value)
    attrs = merge(attributes(data), Dict(key => value))
    setproperty!(PyObject(data), :_attrs, attrs)
    return attrs
end

@doc doc"""
    convert_to_dataset(obj; group = :posterior, kwargs...) -> Dataset

Convert a supported object to a `Dataset`.

In most cases, this function calls [`convert_to_inference_data`](@ref) and returns the
corresponding `group`.
"""
convert_to_dataset

function convert_to_dataset(obj; group=:posterior, kwargs...)
    group = Symbol(group)
    idata = convert_to_inference_data(obj; group=group, kwargs...)
    dataset = getproperty(idata, group)
    return dataset
end
convert_to_dataset(data::Dataset; kwargs...) = data

@doc doc"""
    convert_to_constant_dataset(obj::Dict; kwargs...) -> Dataset
    convert_to_constant_dataset(obj::NamedTuple; kwargs...) -> Dataset

Convert `obj` into a `Dataset`.

Unlike [`convert_to_dataset`](@ref), this is intended for containing constant parameters
such as observed data and constant data, and the first two dimensions are not required to be
the number of chains and draws.

# Keywords

- `coords::Dict{String,Vector}`: Map from named dimension to index names
- `dims::Dict{String,Vector{String}}`: Map from variable name to names of its dimensions
- `library::Any`: A library associated with the data to add to `attrs`.
- `attrs::Dict{String,Any}`: Global attributes to save on this dataset.
"""
convert_to_constant_dataset

function convert_to_constant_dataset(
    obj; coords=Dict(), dims=Dict(), library=nothing, attrs=Dict()
)
    base = arviz.data.base

    obj = _asstringkeydict(obj)
    coords = _asstringkeydict(coords)
    dims = _asstringkeydict(dims)
    attrs = _asstringkeydict(attrs)

    data = Dict{String,PyObject}()
    for (key, vals) in obj
        vals = _asarray(vals)
        val_dims = get(dims, key, nothing)
        (val_dims, val_coords) = base.generate_dims_coords(
            size(vals), key; dims=val_dims, coords=coords
        )
        data[key] = xarray.DataArray(vals; dims=val_dims, coords=val_coords)
    end

    default_attrs = base.make_attrs()
    if library !== nothing
        default_attrs = merge(default_attrs, Dict("inference_library" => string(library)))
    end
    attrs = merge(default_attrs, attrs)
    return Dataset(; data_vars=data, coords=coords, attrs=attrs)
end

@doc doc"""
    dict_to_dataset(data::Dict{String,Array}; kwargs...) -> Dataset

Convert a dictionary with data and keys as variable names to a [`Dataset`](@ref).

# Keywords

- `attrs::Dict{String,Any}`: Json serializable metadata to attach to the dataset, in
    addition to defaults.
- `library::String`: Name of library used for performing inference. Will be attached to the
    `attrs` metadata.
- `coords::Dict{String,Array}`: Coordinates for the dataset
- `dims::Dict{String,Vector{String}}`: Dimensions of each variable. The keys are variable
    names, values are vectors of coordinates.

# Examples

```@example
using ArviZ
ArviZ.dict_to_dataset(Dict("x" => randn(4, 100), "y" => randn(4, 100)))
```
"""
dict_to_dataset

function dict_to_dataset(data; library=nothing, attrs=Dict(), kwargs...)
    if library !== nothing
        attrs = merge(attrs, Dict("inference_library" => string(library)))
    end
    return arviz.dict_to_dataset(data; attrs=attrs, kwargs...)
end

@doc doc"""
    dataset_to_dict(ds::Dataset) -> Tuple{Dict{String,Array},NamedTuple}

Convert a `Dataset` to a dictionary of `Array`s. The function also returns keyword arguments
to [`dict_to_dataset`](@ref).
"""
dataset_to_dict

function dataset_to_dict(ds::Dataset)
    ds_dict = ds.to_dict()
    data_vars = ds_dict["data_vars"]
    attrs = ds_dict["attrs"]

    coords = ds_dict["coords"]
    delete!(coords, "chain")
    delete!(coords, "draw")
    coords = Dict(k => v["data"] for (k, v) in coords)

    data = Dict{String,Array}()
    dims = Dict{String,Vector{String}}()
    for (k, v) in data_vars
        data[k] = v["data"]
        dim = v["dims"][3:end]
        if !isempty(dim)
            dims[k] = [dim...]
        end
    end

    return data, (attrs=attrs, coords=coords, dims=dims)
end