text:
# 0 -> empty space, 1 -> O, 2 -> X
str = [' ', 'O', 'X']

function playTicTacToe(;ai::Int=1,method::String="minimax",abDepth::Int=4)
    # initialise board and display
    M = ttt_initialise()
    ttt_boardDisplay(M)

    # which player(s) are AI?
    if ai == 0
        isAI = [false false]
    elseif ai == 1
        isAI = rand()<0.5 ? [false true] : [true false]
    elseif ai == 2
        isAI = [true true]
    end

    # randomly pick starting player (sets symbols)
    player = rand((1,2))

    while true # take turns until the game ends
        # take a turn, either by human or AI
        if isAI[player]
            M = ttt_aiTurn!(M,player,method,abDepth)
        else
            M = ttt_takeTurn!(M,player)
        end

        ttt_boardDisplay(M) # show the board after the turn

        player = 3-player # swap player
        W = ttt_checkWinner(M) # check if someone has won

        if W!=0 # someone has won
            print("$(str[W+1])'s win\n")
            return nothing
        elseif ttt_gameOver(M) # board filled but no winner -> draw
            print("It's a draw\n")
            return nothing
        end
    end
end

function ttt_initialise()
    # empty board
    return zeros(Int8,3,3)
end

function ttt_boardDisplay(M)
    # Displays the board with symbols
    r1 = "         |   |   \n"
    r2(a) = "Row $a  $(str[M[a,1]+1]) | $(str[M[a,2]+1]) | $(str[M[a,3]+1]) \n"
    r3 = "      ___|___|___\n"
    print("Column 1   2   3\n",r1,r2(1),r3,
        r1,r2(2),r3,r1,r2(3),r1)
        print("Board eval: $(ttt_scoreHelper(M,1))\n")
end

function ttt_takeTurn!(M,player)
    # asks for player input and updates the game state
    coords = ttt_getMove(M)
    M[coords[1],coords[2]] = player
    return M
end

function ttt_getMove(M)
    # gets a move from a human player
    # TODO: More input validation

    while true # repeat until valid location selected
        print("Enter the row/column coordinates for your move separated by a space: ")
        # read a space-separated list of integers
        x = readdlm(IOBuffer(readline()),' ',Int64)
        if M[x[1],x[2]]==0 # must be an empty location
            return x
        end
    end
end

function ttt_aiTurn!(M,player,method,abDepth)
    # work out the best move and then plays it
    moves = 1:9
    if isequal(method,"minimaxAB")
        idx = getBestMoveAB(M,player,ttt_score,moves,ttt_movePossible,ttt_makeMove!,ttt_gameOver,abDepth)
    else
        idx = getBestMove(M,player,ttt_score,moves,ttt_movePossible,ttt_makeMove!,ttt_gameOver)
    end
    M[idx] = player
    return M
end

function ttt_score(M,player)
    win = ttt_checkWinner(M)
    if win!=0
        return win==player ? 100 : -100
    elseif ttt_gameOver(M)
        return 0
    else
        return ttt_scoreHelper(M,player)
    end

    # returns 0 if the game isn't over, 1 if the player wins,
    # -1 if the player loses
end

function ttt_scoreHelper(M,player)
    # See who has won the game

    # p1  and p2 are boolean arrays indicating where each player has played
    p = M.==player
    o = M.==3-player

    s = sum(filter(x -> x>=2 || x<=-2,
        [sum(p-o,dims=1) sum(p-o,dims=2)' sum(p[1:4:9]-o[1:4:9]) sum(p[3:2:7]-o[3:2:7])]))

    return s
end

function ttt_checkWinner(M)
    # See who has won the game

    # p1  and p2 are boolean arrays indicating where each player has played
    p1 = M.==1
    p2 = M.==2

    #   column of 3       or   row of 3         or main diagonal  or other diagonal
    if any(all(p1,dims=1)) | any(all(p1,dims=2)) | all(p1[1:4:9]) | all(p1[3:2:7])
        W = 1
    elseif any(all(p2,dims=1)) | any(all(p2,dims=2)) | all(p2[1:4:9]) | all(p2[3:2:7])
        W = 2
    else
        W = 0 # no winner or draw
    end

    return W
end

function ttt_makeMove!(N,i,player)
    N[i] = player
    return N
end

function ttt_movePossible(M,i)
    return M[i]==0
end

function ttt_gameOver(M)
    # if all the squares are taken (filled with 1s/2s) the game ends
    return all(M.>0) || ttt_checkWinner(M)!=0
end


text:
module NDTools
using Core: add_int
using Base.Iterators, PaddedViews, LinearAlgebra
using Statistics, OffsetArrays

include("offset_types.jl")
include("MutablePaddedViews.jl")
include("datatype_tools.jl")
include("size_tools.jl")
include("iteration_tools.jl")
include("selection_tools.jl")
include("generation_tools.jl")
include("reverse.jl")
include("base_extensions.jl")

end # module

ext: jl


text:
using Jchemo

n = 50 ; p = 7 ; q = 2
Xtrain = rand(n, p) ; Ytrain = rand(n, q)
ytrain = Ytrain[:, 1]
m = 3
Xtest = rand(m, p) ; Ytest = rand(m, q)
ytest = Ytest[:, 1]

########### PLSR

nlv = 0:2
pars = mpars(nlv = nlv)
gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = plskern, pars = pars)

# Faster
gridscorelv(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = plskern, nlv = nlv)

########### KPLSR

nlv = 0:2
pars = mpars(nlv = nlv, gamma = [1; 10])
#pars = mpars(nlv = 0:nlv, kern = ["kpol";], degree = 1:2)
gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = kplsr, pars = pars)

# Faster
pars = mpars(gamma = [1, 10])
gridscorelv(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = kplsr, nlv = nlv, pars = pars)

########### DKPLSR

nlv = 0:2
pars = mpars(nlv = nlv, gamma = [1; 10])
gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = dkplsr, pars = pars)

# Faster
pars = mpars(gamma = [1; 10])
gridscorelv(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = dkplsr, nlv = nlv, pars = pars)

############ RR

lb = [.01; .1]
pars = mpars(lb = lb)
gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = rr, pars = pars)

gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = rrchol, pars = pars)

# Faster (only for rr, not for rrchol)
gridscorelb(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = rr, lb = lb)

############ KRR

lb = [.01; .1]
pars = mpars(lb = lb, gamma = [1; 10])
gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = krr, pars = pars)

# Faster
pars = mpars(gamma = [1; 10])
gridscorelb(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = krr, lb = lb, pars = pars)

############ kNN-R

nlvdis = 5 ; metric = ["mahal"; "eucl"]
h = [1; 3] ; k = [20; 10]
pars = mpars(nlvdis = nlvdis, metric = metric, h = h, k = k)
gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = knnr, pars = pars, verbose = true)

############ kNN-LWPLSR

nlvdis = 5 ; metric = ["mahal";]
h = [1.; 3.] ; k = [100; 20]
nlv = 1:2
pars = mpars(nlv = nlv, nlvdis = nlvdis, metric = metric, h = h, k = k)
gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = lwplsr, pars = pars, verbose = true)

# Faster
pars = mpars(nlvdis = nlvdis, metric = metric, h = h, k = k)
gridscorelv(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = lwplsr, nlv = nlv, pars = pars, verbose = false)

########### PLSR-AGG

# Here there is no sense to use gridscorelv

pars = mpars(nlv = ["1:2"; "1:3"])
gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = plsr_agg, pars = pars)

############ kNN-LWPLSR-AGG

# Here there is no sense to use gridscorelv

nlvdis = 5 ; metric = ["mahal";] ;
h = [1.; 3.] ; k = [20; 10]
nlv = ["1:2"; "2:5"] ;
pars = mpars(nlv = nlv, nlvdis = nlvdis, metric = metric, h = h, k = k) ;
gridscore(Xtrain, Ytrain, Xtest, Ytest;
    score = msep, fun = lwplsr_agg, pars = pars, verbose = true)


text:
using Thebes, Luxor, ColorSchemes

function makecube()
    cube = [
        Point3D(1,   1, -1),
        Point3D(1,  -1, -1),
        Point3D(-1, -1, -1),
        Point3D(-1,  1, -1),
        Point3D(1,   1,  1),
        Point3D(1,  -1,  1),
        Point3D(-1, -1,  1),
        Point3D(-1,  1,  1)]
    r = Point3D[]

    for e in (
            [1, 2, 3, 4, 1],
            [5, 6, 7, 8, 5],
            [5, 1, 2, 6, 7],
            [7, 3, 4, 8, 5])
        append!(r, cube[e])
    end
    return r
end

function frame(scene, framenumber)
    eyepoint(200, 180, 100)
    background("lightblue1")
    sethue("brown")
    setlinejoin("bevel")
    carpet()
    axes3D()
    fontsize(30)
    p = rescale(framenumber, 1, scene.framerange.stop, 0, 1500)
    if p < 300
        p = 0
    end
    perspective(p)
    sethue("black")
    pin(10makecube(), gfunction = (p3, p2) -> poly(p2, :stroke))
    text(string(round(p, digits=0)), boxtopleft(BoundingBox() * 0.9))
end

width, height = (400, 400)
perspectivemovie = Movie(width, height, "perspectivemovie")

animate(perspectivemovie, [
    Scene(perspectivemovie, frame, 1:150)
    ],
    creategif=true,
    framerate=24,
    pathname="/tmp/perspectivemovie-animation.gif")


text:
"""
    UNetDynamic(backbone, inputsize[; kwargs...])

Create a U-Net model from convolutional `backbone` architecture. After every
downsampling layer (i.e. pooling or strided convolution), a skip connection and
an upsampling block are inserted, resulting in a convolutional network with
the same spatial output dimensions as its input.

## Keyword arguments

- `upsample`: A *constructor* for an upsampling block callable with `upsample(insize, k_out)`.
    If `insize` is `(h, w, k, b)`, then the output should have size `(2h, 2w, k_out, b)`.
    Defaults to [`FastAI.Models.upsample_block_small`](#).
- `agg`: Aggregation function for skip connection. Default concatenates in the
  channel dimension. Use `+` for summing and see [`Flux.SkipConnection`](#) for more
  details.
- `fdownsample = 0`: Number of upsampling steps to leave out. By default there will be one
    upsampling step for every downsampling step in `backbone`. Hence if the input spatial
    size is `(h, w)`, the output size will be `(h/2^fdownsample, w/2^fdownsample)`, i.e.
    to get outputs at half the resolution, set `fdownsample = 1`.
- `kwargs...`: Other keyword arguments are passed through to `upsample`.

## Examples

```julia
using FastAI, Metalhead

backbone = Metalhead.ResNet50(pretrain=true).layers[1:end-3]
unet = UNetDynamic(backbone, (256, 256, 3, 1); k_out = 10)
Flux.outputsize(unet, (256, 256, 3, 1)) == (256, 256, 10, 1)

unet = UNetDynamic(backbone, (256, 256, 3, 1); fdownscalk_out = 10)
Flux.outputsize(unet, (256, 256, 3, 1)) == (256, 256, 10, 1)
```
"""
function UNetDynamic(backbone, inputsize, final; kwargs...)
    backbonelayers = collect(iterlayers(backbone))
    unet = unet_from_layers(backbonelayers, inputsize; kwargs...)
    outsz = Flux.outputsize(unet, inputsize)
    return Chain(unet, final(outsz))
end

function UNetDynamic(backbone, inputsize, k_out::Int; kwargs...)
    final = insz -> convxlayer(insz[end-1], k_out; ks = 1)
    return UNetDynamic(backbone, inputsize, final; kwargs...)
end

function unet_from_layers(
        backbonelayers,
        insz;
        fdownscale = 0,
        upsample = upsample_block_small,
        agg = (mx, x) -> cat(mx, x; dims = length(insz)-1),
        kwargs...)
    layers = []
    channeldim = length(insz) - 1

    for (i, layer) in enumerate(backbonelayers)
        outsz = Flux.outputsize(layer, insz)
        if (insz[1] ÷ 2 == outsz[1])
            if fdownscale == 0
                child_unet = unet_from_layers(
                    backbonelayers[i+1:end],
                    outsz;
                    upsample = upsample,
                    fdownscale = fdownscale,
                    agg = agg)
                outsz = Flux.outputsize(child_unet, outsz)

                upsample_k_out = insz[channeldim]
                up = upsample(outsz, upsample_k_out; kwargs...)

                push!(layers, SkipConnection(Chain(layer, child_unet, up), agg))

                break
            else
                fdownscale -= 1
            end
        end
        push!(layers, layer)
        insz = outsz
    end

    unet = length(layers) == 1 ? only(layers) : Chain(layers...)
    return unet
end

iterlayers(m::Chain) = Iterators.flatten(iterlayers(l) for l in m.layers)
iterlayers(m) = (m,)

"""
    upsample_block_small(insize, k_out)

An upsampling block that increases the spatial dimensions of the input by 2
using pixel-shuffle upsampling.
"""
function upsample_block_small(insize, k_out; ks = 3, kwargs...)
    return Chain(
        Flux.PixelShuffle(2),
        convxlayer(insize[end-1] ÷ 4, k_out; kwargs...)
    )
end

function conv_final(insize, k_out; ks = 1, kwargs...)
    return convxlayer(insize[end-1], k_out; ks = ks, kwargs...)
end


text:
include("../src/autodiff_module.jl")

module dtest
using ..autodiff
using PyPlot
PyPlot.pygui(true)

fcn0_1(x) = 5 * x + 10.0

fcn1(x) = begin
    c1 = x^3
    c2 = x^2
    c3 = x
    t1 = 1/3 * c1
    return t1 + 1/2 * c2 + 5 * c3 + 10.0
end
fcn2(x) = sin(x) * exp(x) + x^5
fcn3(x) = sin(x)

mode = autodiff.ReverseMode()

# x0 = 2.1 # 5.0 #-1.5 # 4.6 # 3.2
# der_f = autodiff.derivative(autodiff.ForwardMode(), fcn1, x0)
# der_r = autodiff.derivative(autodiff.ReverseMode(), fcn1, x0)

x0 = collect(-3.0:0.1:5)
L = length(x0)
der_f = Vector{Float64}(undef, L)
der_r = Vector{Float64}(undef, L)
delta_der = Vector{Float64}(undef, L)

# der_f = zeros(Float64, L)
# der_r = zeros(Float64, L)
# delta_der = zeros(Float64, L)

for ii = 1:L
    der_r[ii] = autodiff.derivative(autodiff.ReverseMode(), fcn3, x0[ii])
end

for ii = 1:L
    der_f[ii] = autodiff.derivative(autodiff.ForwardMode(), fcn3, x0[ii])
end

delta_der = der_f - der_r

PyPlot.figure()
PyPlot.plot(x0, der_f, marker = :+)
PyPlot.plot(x0, der_r, marker = :o)

PyPlot.figure()
PyPlot.plot(x0, delta_der)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
#=
    The file contains functions to operate modular reconstrction and CRT
=#

#------------------------------------------------------------------------------

# Rational number reconstruction implementation borrowed from CLUE
# and modified a bit to suit the 'Modern Computer Algebra' definitions
# Returns a rational r // h of QQ field in a canonical form such that
#   r // h ≡ a (mod m)
#
# let n = max( λ(a), λ(m) ) , where λ(x) is a number of bits for x
# O(n^2)
function rational_reconstruction(a::I, m::I) where {I<:Union{Int, BigInt}}
    a = mod(a, m)
    if a == 0 || m == 0
        return QQ(0, 1)
    end
    if m < 0
        m = -m
    end
    if a < 0
        a = m - a
    end
    if a == 1
        return QQ(1, 1)
    end
    bnd = sqrt(float(m) / 2)

    # an issue: we can not make good use of staticarrays here
    # as long as BigInts are heap allocated anyways
    U = (I(1), I(0), m)
    V = (I(0), I(1), a)
    while abs(V[3]) >= bnd
        q = div(U[3], V[3])
        # TODO: use MutableArithmetics
        T = U .- q .* V
        U = V
        V = T
    end

    t = abs(V[2])
    r = V[3] * sign(V[2])
    # changed from `<= bnd` to `<= m / bnd`
    # we can speed up this !
    if t <= m / bnd && gcd(r, t) == 1
        return QQ(r, t)
    end

    # TODO: not needed
    throw(DomainError(
        :($a//$m), "rational reconstruction of $a (mod $m) does not exist"
    ))

    return QQ(0, 1)
end

function reconstruct_modulo(
        coeffs_zz::Vector{Vector{BigInt}},
        modulo::BigInt)

    coeffs_qq = Vector{Vector{Rational{BigInt}}}(undef, length(coeffs_zz))
    for i in 1:length(coeffs_zz)
        poly_zz = coeffs_zz[i]
        coeffs_qq[i] = Vector{Rational{BigInt}}(undef, length(poly_zz))
        for j in 1:length(poly_zz)
            coeffs_qq[i][j] = rational_reconstruction(poly_zz[j], modulo)
        end
    end

    coeffs_qq
end

#------------------------------------------------------------------------------

function CRT(
            a1::BigInt, m1::BigInt,
            a2::U, m2::BigInt) where {U<:Integer}

    CRT(a1, m1, BigInt(a2), m2)
end

function CRT(a1::BigInt, m1::BigInt, a2::BigInt, m2::BigInt)
    M = m1 * m2
    Ms1, Ms2 = m2, m1
    t1 = invmod(Ms1, m1)
    t2 = invmod(Ms2, m2)
    mod(a1*t1*Ms1 + a2*t2*Ms2, M)::BigInt
end

# TODO: move to coeffs.jl
function reconstruct_crt!(
            gbcoeffs_ff::Vector{Vector{UInt64}}, ch::UInt64)

    gbcoeffs_zz = Vector{Vector{BigInt}}(undef, length(gbcoeffs_ff))
    bigch = BigInt(ch)
    for i in 1:length(gbcoeffs_ff)
        poly_ff = gbcoeffs_ff[i]
        gbcoeffs_zz[i] = Vector{BigInt}(undef, length(poly_ff))
        for j in 1:length(poly_ff)
            gbcoeffs_zz[i][j] = BigInt(poly_ff[j])
        end
    end
    gbcoeffs_zz, bigch
end

function reconstruct_crt!(
            gbcoeffs_accum::Vector{Vector{BigInt}}, modulo::BigInt,
            gbcoeffs_ff::Vector{Vector{UInt64}}, ch::UInt64)

    if isempty(gbcoeffs_accum)
        return reconstruct_crt!(gbcoeffs_ff, ch)
    end

    gbcoeffs_zz = Vector{Vector{BigInt}}(undef, length(gbcoeffs_ff))
    bigch = BigInt(ch)

    for i in 1:length(gbcoeffs_ff)
        poly_zz = gbcoeffs_accum[i]
        poly_ff = gbcoeffs_ff[i]
        gbcoeffs_zz[i] = Vector{BigInt}(undef, length(poly_ff))
        for j in 1:length(poly_ff)
            gbcoeffs_zz[i][j] = CRT(poly_zz[j], modulo, poly_ff[j], bigch)
        end
    end
    # TODO: MutableArithmetics
    gbcoeffs_zz, modulo * bigch
end

#------------------------------------------------------------------------------


text:
# This file was generated, do not modify it. # hide
r = report(mtm)

md = r.parameter_values[:,1]
mcw = r.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(md, mcw, r.measurements)

xlabel("Maximum tree depth", fontsize=14)
ylabel("Minimum child weight", fontsize=14)
xticks(3:2:10, fontsize=12)
yticks(fontsize=12)

savefig("assets/literate/EX-crabs-xgb-heatmap.svg") # hide


text:
using Random
using Zygote: @adjoint
using Zygote
using ForwardDiff
using DSP

Random.seed!(1)

Float64(x::ForwardDiff.Dual)=ForwardDiff.value(x)

function dconv(x, a, b)
    n = length(size(a))
    x = reshape(vec(x), (size(a) .+ size(b) .- ones(Int, n))...)
    r = (
        DSP.conv(x, reverse(b))[(
            # DSP.xcorr(x, b, padmode = :none)[(
            i:j for (i, j) in zip(size(b), size(x))
        )...],
        DSP.conv(x, reverse(a))[(
            # DSP.xcorr(x, a, padmode = :none)[(
            i:j for (i, j) in zip(size(a), size(x))
        )...],
    )
    return r
end
# @adjoint DSP.conv(a, b) = DSP.conv(a, b), x -> dconv(x, a, b)
@adjoint DSP.conv(a, b) = DSP.conv(Float64.(a), Float64.(b)),
x ->
    dconv(Float64.(x), Float64.(a), Float64.(b))
# function frule(
#     (_, ΔA, ΔB),
#     ::typeof(DPS.conv),
#     A,
#     B,
# )
#     Ω = conv(A , B)
#     ∂Ω = conv(ΔA * B )+ A * ΔB
#     return (Ω, ∂Ω)
# end
@show jacobian(DSP.conv, [1, 2], [3, 4])

const tol = 1e-9

mutable struct LinearOperator
    i::Int
    o::Int
    li::Int
    lo::Int
    lr::Int

    dx::Float64
    rmax::Float64
    sz::Any

    params::AbstractVector

    radfunc::Function
    conv::Function

    H::AbstractArray
end

function remake!(m::LinearOperator; dx = m.dx, rmax = m.rmax, sz = m.sz)
    m.H = makefilters(m.lr, dx, rmax, sz, m.radfunc, m.params)
end

invr =
    (r; p = 1, rmax = 1.0e20, kwargs...) ->
        abs(r) < tol || r > rmax + tol ? 0.0 : 1 / r^p

const FD = [[1.0], [0.0, 0.5]]
function δ(x, dx; n = 0)
    i = round(Int, x / dx)
    if abs(x - i) < tol
        return i < length(FD[n+1]) ? FD[n+1][i+1] : 0.0
    end
    return 0.0
end

function LinearOperator(name, dx; kwargs...)
    LinearOperator(dx; name = name, kwargs...)
end
function LinearOperator(
    dx;
    i::Int = -1,
    o::Int = -1,
    li = -1,
    lo = -1,
    lr = -1,
    name = nothing,
    radfunc = g,
    params = ones(8),
    rmax = 1 / tol,
    sz = nothing,
)
    dV = dx^3

    if name == :potential
        li = lo = lr = 0
        params = []
        radfunc = invr
    elseif name == :field
        li = 0
        lr = lo = 1
        params = []
        radfunc = (r; kwargs...) -> invr(r; p = 2)
    elseif name == :grad
        li = 0
        lr = lo = 1
        params = []
        rmax = dx
        radfunc = (r; kwargs...) -> δ(r, dx; n = 1) / dV
    end

    if rmax !== nothing
        nrmax = round(Int, rmax / dx)
    end
    if sz === nothing
        sz = (1 + 2nrmax) * ones(Int, 3)
    end

    conv = Conv(li, lr, lo)

    m = LinearOperator(
        i,
        o,
        li,
        lo,
        lr,
        dx,
        rmax,
        sz,
        params,
        radfunc,
        conv,
        [],
    )
    remake!(m)
    return m
end

function (m::LinearOperator)(X; remake = false)
    if remake
        remake!(m)
    end
    m.conv(X, m.H)
end

"""
radial function
"""
function g(r; params = ones(16), dx = 1.0)
    k, σ = params[end-2+1:end]
    dV = dx^3
    funcs = [
        abs(r) < tol ? 1 / dV : 0.0,
        -invr(r; rmax = dx) / dV,
        invr(r),
        invr(r; p = 2),
        exp(-abs(k) * r),
        exp(-(r / σ)^2),
    ]
    dot(params[1:length(funcs)], funcs)
end

const s3 = sqrt(3)
function makefilters(lr, dx, rmax, sz, radfunc, params;)
    dV = dx^3
    nrmax = round(Int, rmax / dx)
    center = (sz .+ 1) .÷ 2

    V = 0
    Zygote.ignore() do
        V = [[x, y, z] .- center for x = 1:sz[1], y = 1:sz[2], z = 1:sz[3]] * dx
    end
    R0 = norm.(V)

    radfunc1 =
        r -> r > rmax + tol ? 0.0 : radfunc(r; params = params, dx = dx)
    R = radfunc1.(R0) * dV

    if lr == 0
        ret = [R]
    else
        Vhat = V ./ broadcast(x -> abs(x) < tol ? 1.0 : x, R0)
        X = getindex.(Vhat, 1)
        Y = getindex.(Vhat, 2)
        Z = getindex.(Vhat, 3)
        # println(X)
        # println(R)
        if lr == 1
            ret = [R .* X, R .* Y, R .* Z]
        elseif lr == 2
            ret = broadcast(
                x -> R .* x,
                [
                    s3 * X .* Z,
                    s3 * Y .* X,
                    Y .^ 2 - (X .^ 2 + Z .^ 2) / 2,
                    s3 * Y .* Z,
                    (Z .^ 2 - X .^ 2) * s3 / 2,
                ],
            )
        end
    end
    # [x ./ sum(abs.(x)) for x in ret]
    ret
end

function conv(x, f; crop = 1, periodic = false)
    sx = 0
    sf = 0
    ix = 0
    # starts = 0
    Zygote.ignore() do
        sx = size(x)
        sf = size(f)
        if crop==2
            sx,sf=sf,sx
        end
        ix = [((i-1)÷2+1):((i-1)÷2+j) for (i, j) in zip(sf, sx)]
    end

    y = DSP.conv(x, f)

    if periodic
        for i = 1:3

        end
    else
        return y[ix...]
    end
end

c0__ = (X, F; kwargs...) -> [conv(X[1], f; kwargs...) for f in F]
c111 =
    (X, Y; kwargs...) -> [
        conv(X[2], Y[3]; kwargs...) - conv(X[3], Y[2]; kwargs...),
        conv(X[3], Y[1]; kwargs...) - conv(X[1], Y[3]; kwargs...),
        conv(X[1], Y[2]; kwargs...) - conv(X[2], Y[1]; kwargs...),
    ]
c__0 = (X, F; kwargs...) -> sum(conv.(X, F; kwargs...))
c121 =
    (X, Y; kwargs...) -> [
        conv(X[1], -Y[3] / s3 .- Y[5]; kwargs...) +
        conv(X[2], Y[2]; kwargs...) +
        conv(X[3], Y[1]; kwargs...),
        conv(X[1], Y[2]; kwargs...) +
        conv(X[2], 2 / s3 * Y[3]; kwargs...) +
        conv(X[3], Y[4]; kwargs...),
        conv(X[1], Y[1]; kwargs...) +
        conv(X[2], Y[4]; kwargs...) +
        conv(X[3], -Y[3] / s3 .+ Y[5]; kwargs...),
    ]

function Conv(l1, l2, l; crop = 1, kwargs...)
    if l1 <= l2
        kwargs = (crop = crop,)
        pre = f -> ((X, Y) -> f(X, Y; kwargs...))
    else
        l1, l2 = l2, l1
        crop = [2, 1][crop]
        kwargs = (crop = crop,)
        pre = f -> ((X, Y) -> f(Y, X; kwargs...))
    end

    if l1 == 0
        return pre(c0__)
    elseif l1 == l2 && l == 0
        return pre(c__0)
    elseif l1 == l2 == l == 1
        return pre(c111)
    elseif (l1, l2, l) == (1, 2, 1)
        return pre(c121)
    end
end

function conv(x, y, l)
    Conv(length(x), length(y), l)(x, y)
end

function toArray(x::Vector)
    vcat(vec.(x)...)
    # cat(x...; dims = 4)
    # cat([reshape(x,(size(x)...,1)) for x in x]...;dims=4)
end

# using GLMakie
# using CairoMakie

function tfplot(X)
    # fig = Figure()
    fig = Figure(backgroundcolor = RGBf0(0.98, 0.98, 0.98),
    resolution = (1000, 700))
 for i = 1:length(X)
    volume(fig[1,i],X[i])
end
    # plots = [volume(X[i]) for i = 1:length(X)]
    # plot(plots..., layout = length(X))
end
# p

max_stars_repo_name: paulxshen/eno

max_issues_repo_name: paulxshen/eno

max_forks_repo_name: paulxshen/eno


text:
using HyperbolicPlane, Plots

"""
`complete_graph(n,rad=2)` draws the complete graph
on `n` vertices in the Hyperbolic plane (using straight line
segments for edges). The vertices are distance `rad` from the
origin.
"""
function complete_graph(n::Int=5, radius::Real = 2)

    if n <= 0
        throw(DomainError(n,"n must be positive"))
    end

    pts = [ HPoint(radius, 2pi*k/n) for k=1:n ]
    Plots.plot()
    for j=1:n-1
        for k=j+1:n
            L = HSegment(pts[j], pts[k])
            set_line_style(L,:dot)
            set_color(L,:darkslategray)
            set_thickness(L, 2)
            draw(L)
        end
    end

    for j=1:n
        P = pts[j]
        set_radius(P,2)
        draw(P)
    end
    draw(HPlane())
    finish()
end

complete_graph()

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# main.jl (v0.1)

using DndDB

function main()

	path::String = "D:\\git\\JuliaModules\\DndDB\\data\\db";      # Path to current directory
	lookup::String = "enemynames.json";    # Path to lookup for quicker db-indexing
	tablepath::String = "enemies.json"; # Path to the table of the database

	cd(path);

	uniqueenemies::Array{String, 1} = readfile(lookup);
	enemytable::Enemyarray = readfile(tablepath);
	tmptable = deepcopy(enemytable);
	push!(tmptable, deepcopy(tmptable[1]));

	deleteall!(x -> x["name"] == "Goblin", tmptable)

	# end of main
	writefile(lookup, uniqueenemies);
	writefile(tablepath, enemytable);
end

main();


text:
module ELM
using Reexport
@reexport using DataFrames
export HiddenLayer, ExtremeLearningMachine,SVDCM
export activations, fit!, predict, fitSVD!
include("CM_SVD.jl")
include("initializers.jl")
include("activation_fucntions.jl")
include("main.jl")
end


text:
export non0hist, renyi
using StatsBase

"""
```julia
non0hist(ε, dataset)
non0hist(ε, vectors...)
```
Partition a data-set into tabulated intervals (boxes) of
size `ε` and return the *sum-normalized* histogram in a 1D form, discarding all
non-zero elements. Use the `fit(Histogram, ...)` from package `StatsBase` if you
wish to keep information about the edges of the binning as well as the zero elements.
"""
function non0hist end
@inbounds function non0hist(ε, vectors::Vararg{AbstractVector{T}}) where {T<:Real}
  # Initialize:
  D = length(vectors); L = length(vectors[1])
  for i in 1:D
  length(vectors[i]) == L ||
  throw(ArgumentError("All vectors given to `non0hist` must be of equal length!"))
  end

  if T == BigFloat || T == BigInt
    ranges = StepRangeLen{T, T, T}[]
  else
    ranges = StepRangeLen{T, Base.TwicePrecision{T},Base.TwicePrecision{T}}[]
  end
  # Create ranges:
  for v in vectors
    push!(ranges, minimum(v):ε:maximum(v)+ε)  #be sure to have that +ε at the end!
  end
  # Perform Histogram:
  # Write my OWN Version, which would be many times over faster,
  # since I don't need the 0 elements.
  pks = fit(Histogram, vectors, (ranges...), closed=:left).weights/L
  return T[x for x in pks if x != 0]
end

"""
```julia
renyi(α, ε, dataset)
renyi(α, ε, vectors...)
```
Compute the `α` order generalized Rényi entropy [1] of a dataset,
by first partitioning them into boxes of of size `ε` (log base-e is used).

```julia
renyi(α, p::AbstractArray)
```
Compute the `α` order generalized Rényi entropy of an Array `p` directly,
assuming that `p` is sum-normalized.

The Rényi entropy `R_α(p) = 1/(1-α) * sum_i(p_i^α)` generalizes other known entropies,
like e.g. the Shannon entropy
(α = 1, see *the* Shannon paper [2]), the Hartley entropy (α = 0, also known as
maximum entropy), or the Collision (or correlation) entropy (α = 2).

[1] : A. Rényi, Proceedings of the fourth Berkeley Symposium on Mathematics,
Statistics and Probability, pp 547 (1960)

[2] : C. E. Shannon, Bell System Technical Journal **27**, pp 379 (1948)
"""
function renyi(α, ε, vectors::Vararg{AbstractVector{T}}) where {T<:Real}
  p = non0hist(ε, vectors...)
  renyi(α, p)
end
function renyi{T<:Real}(α::Real, p::AbstractArray{T})
  α < 0 && throw(ArgumentError("Order of Rényi entropy must be ≥ 0."))

  if α ≈ 0
    return log(length(p)) #Hartley entropy, max-entropy
  elseif α ≈ 1
    return -sum( x*log(x) for x in p ) #Shannon entropy, information to locate with ε
  elseif (isinf(α))
    return -log(maximum(p)) #Min entropy
  else
    return (1/(1-α))*log( sum(x^α for x in p) )
  end
end

# # Linear fit on dimension:
# path = "C:\\Users\\datseris\\PowerFolders\\MyScience\\Lectures\\Introduction to Physics of Complex Systems\\Exercises 2016\\10 - Generilized Dimensions And Synchronization\\data1.txt"
# f = readdlm(path)
# f1 = view(f, :, 1); f2 = view(f, :, 2)
#
# es = logspace(-1, -3, 6)
# D = zeros(es)
# for (i, ε) in enumerate(es)
#   println("ε = $ε")
#   D[i] = renyi(1, ε, f1, f2)
# end
#
# using LsqFit
# model(x, p) = p[1].*x .+ p[2]
#
# cf = curve_fit(model, -ln.(es), D, rand(2))
# plot(lnes, -cf.param[1] .* ln.(es) .+ cf.param[2])


text:
using DataDeps

register(DataDep(
  "lpga-2008",
  "http://www.stat.ufl.edu/~winner/data/lpga2008.txt",
  "http://www.stat.ufl.edu/~winner/data/lpga2008.dat",
  "3d631844fb4b13f7d1f9205005d5ef3687a383b515a39c027e41e5124883ed14",
  post_fetch_method=(path -> begin
    # IDs have spaces in them which is the delimiter
    lines = readlines(path)
    map!(x -> lstrip(x[30:end]), lines, lines)
    write(path, join(lines, '\n'))

    UCIData.process_dataset(path,
        target_index=8,
        feature_indices=1:6,
        id_indices=10,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# this page is going to hold all of the setting-up functions
# we will need to know how many neurons are in each layer and how many connections there are between each layer

#-------------------------------------------------------------------------------------------------------#
#										Base Redefinitions
#-------------------------------------------------------------------------------------------------------#

#				Connection
#===========================================#
Base.size(A::Connection) = A.dims

Base.getindex(A::Connection, I::Any) = A.data[I]
Base.getindex(A::Connection, I::Vararg{Any,N}) where {N} = A.data[I...]
Base.setindex!(A::Connection, filler, I::Any) = A.data[I...] .= filler
Base.setindex!(A::Connection, filler, I::Vararg{Any,N}) where {N} = A.data[I...] .= filler
Base.setindex!(A::Connection, filler, I::Vararg{Int,N}) where {N} = A.data[I...] .= filler
Base.display(A::Connection) = display(A.data)

Base.:(*)(A::T, B::NormProduct) where {T<:Connection} = T(A.data .* B)
Base.:(-)(A::T, B::NormProduct) where {T<:Connection} = T(A.data .- B)
Base.:(+)(A::T, B::NormProduct) where {T<:Connection} = T(A.data .+ B)
Base.:(/)(A::T, B::NormProduct) where {T<:Connection} = T(A.data ./ B)

Base.:(*)(A::T, B::T) where {T<:Connection} = T(A.data .* B.data)
Base.:(-)(A::T, B::T) where {T<:Connection} = T(A.data .- B.data)
Base.:(+)(A::T, B::T) where {T<:Connection} = T(A.data .+ B.data)
Base.:(/)(A::T, B::T) where {T<:Connection} = T(A.data ./ B.data)

#				SparseConnection
#===========================================#
SparseConnection(A::AbstractSparseArray{T,N}) where {T,N<:Int} = SparseConnection(A, size(A))
SparseConnection(A::SparseConnection{T,N}) where {T,N} = SparseConnection(A.data, size(A))

#		NeuronLayer & SynapseLayer
#===========================================#
Base.size(A::T) where {T<:ConnectionLayer} = A.dims
Base.getindex(A::NeuronLayer{T,N}, Ind::Int) where {T,N} = get(A.data, Ind, nothing)
Base.setindex!(A::NeuronLayer, filler, Ind::Int) where {T,N} = (A.data[Ind] = filler)
Base.zero(::Type{NeuronLayer{T,N}}) where {T<:Number,N} = NeuronLayer{T,N}([zero(T)],(1,))
 NeuronLayer(::Type{T}, len::Integer) where {T} = NeuronLayer{T,1}(Array{T,1}(undef, len), (len,))
NeuronLayer(filler::T, len::Integer) where T = NeuronLayer{T,1}(fill(filler,len), (len,))
Base.copy(A::ConnectionLayer) = typeof(A)(copy(A.data),A.dims)
Base.copy(A::BrainTypes) = typeof(A)([copy(l) for l in A.layers])
average(A::BrainTypes) = mean([average(lay) for lay in A.layers])
average(A::ConnectionLayer) = sum(A.data) == 0 ? 0 : sum(A.data)/length(A.data)
average(A::Array{<:BrainTypes,1}) = mean([average(elem) for elem in A])

Base.zero(::Type{SynapseLayer{T,N}}) where {T<:Number, N} = SynapseLayer{T,N}(zeros(Int,ntuple(x->1,N)),ntuple(x->1,N))
Base.getindex(A::ConnectionLayer, I::Vararg) = A.data[I...]
Base.setindex!(A::ConnectionLayer, filler, I::Vararg) = A.data[I...] = filler
Base.length(bt::BrainTypes) = length(bt.layers)

Base.one(::Type{SynapseLayer{T,N}}) where {T,N} = one(T)

SynapseLayer{T,N}(A::Array{T,N}) where {T,N} = SynapseLayer{T,N}(A,size(A))
SynapseLayer(A::Array{T,N}) where {T,N} = SynapseLayer{T,N}(A)
SynapseLayer(filler::T, dims::NTuple{N,Int}) where {T,N} = SynapseLayer{T,N}(fill(filler, dims))

MBONLayer(A::Array{T,N}) where {T,N} = MBONLayer(A, size(A))
#-------------------------------------------------------------------------------------------------------#
#										Initialising Functions
#-------------------------------------------------------------------------------------------------------#
function create_MBON_synapses(lyrsize::NTuple{2,Int}; weight=20.0)

	num_MBONS = lyrsize[2]
	synapses = [SynapseLayer(zeros(lyrsize...)) for _ in 1:num_MBONS]
	for l in 1:num_MBONS
		synapses[l][:,l] .= weight
	end
	return MBONLayer(synapses, size(synapses))
end

function create_MBON_synapses(nn::Iterable, syndens::Iterable, weight::Iterable)

	lyrsizes = [(nn[i],nn[i+1]) for i = 1:length(nn)-1]
		# big horrible list comprehension to get the normal synapses, then the weird synapses
	out = ([i==length(lyrsizes) ? create_MBON_synapses(lyr, weight=weight[i]) : create_synapses(SynapseLayer, lyr, syndens[i], weight=weight[i]) for (i, lyr) in enumerate(lyrsizes)])
	DiverseLayers(out)
end

function create_synapses(::Type{SynapseLayer}, lyrsize::NTuple{2,Int}, syndens::Float64; weight=20.0)
	syndens==1 && return SynapseLayer(ones(lyrsize...))
	syndens==0 && return SynapseLayer(zeros(lyrsize...))
	ns = Int(round(syndens*lyrsize[1]))
	return create_synapses(SynapseLayer, lyrsize, ns, weight=weight)
end

function create_synapses(::Type{SynapseLayer}, lyrsize::NTuple{2,Int}, syndens::Int; weight=20.0)

	syns = SynapseLayer(zeros(lyrsize...))
	cons = Array{Int,2}(undef,syndens,lyrsize[2])
	for i = 1:lyrsize[2]
		cons[:,i] .= shuffle(1:lyrsize[1])[1:syndens]
	end

	for postlyr = 1:lyrsize[2]
		syns[cons[:,postlyr],postlyr] .= weight
	end
	return syns
end

function create_synapses(::Type{SynapseLayers}, nn::Array, syndens; weight=2)

	lyrsizes = [(nn[i],nn[i+1]) for i = 1:length(nn)-1]
	SynapseLayers([create_synapses(SynapseLayer, lrs) for lrs in lyrsizes])
end

function create_synapses(::Type{SynapseLayers}, nn::Iterable, syndens::Iterable, weight::Iterable)
	lyrsizes = [(nn[i],nn[i+1]) for i = 1:length(nn)-1]
	out = Array{SynapseLayer,1}(undef,length(nn)-1)
	for (i, ls) in enumerate(lyrsizes)
		out[i] = create_synapses(SynapseLayer, ls, syndens[i], weight=weight[i])
	end
	return SynapseLayers{Float64,2}(out)
end

function clone_synapses(SLS::T; filler=1) where {T<:BrainTypes}

	#T([typeof(l)((l.data .> 0) .* filler) for l in SLS.layers])

	out = [clone_synapses(l, filler=filler) for l in SLS.layers]
	return T(out)
end

function clone_synapses(SLS::T; filler=1) where {T<:SynapseLayer}

	SynapseLayer((SLS.data .> 0) .* one(eltype(SLS)) * filler)
end

function clone_synapses(SLS::T; filler=1) where {T<:MBONLayer}

	# basically just clone the lower layers.
	MBONLayer([clone_synapses(l, filler=filler) for l in SLS.data])
end

function fill_synapses(::Type{SynapseLayers}, nn::Array, filler::T) where T
	lyrsizes = [(nn[i],nn[i+1]) for i = 1:length(nn)-1]
	SynapseLayers([SynapseLayer(filler, lrs) for lrs in lyrsizes])
end

function Base.fill!(a::BrainTypes, filler::Iterable)

	for (i,layer) in enumerate(a.layers)
		fill!(layer, filler[i])
	end
end

function Base.fill!(a::BrainTypes, filler)

	for layer in a.layers
		fill!(layer, filler)
	end
end
Base.fill!(a::ConnectionLayer, filler) = fill!(a.data, filler)

function fillentries!(sparray, filler)

	for l = 1:length(sparray)

		(i, j, v) = findnz(sparray[l])
		sparray[l] = sparse(i,j, v .* filler)

	end

	return sparray
end

function fillcells!(arraytofill, filler)

	for l = 1:length(arraytofill)

		arraytofill[l] .= filler

	end
end

function create_neurons(nn, filler)

	out = []

	if length(filler) < length(nn)
		filler = fill(filler,length(nn))
	end

	for layer = 1:length(nn)

		push!(out, fill(filler[layer], nn[layer]))

	end
	return out
end

#-------------------------------------------------------------------------------------------------------#


text:
module TikZ
using TikzPictures

function indent(s::String)
    map(split(s, "\n")) do line
        "  "*line
    end |> l -> join(l, "\n")
end

const TikZarg = Union{String,<:Pair{String,<:Any}}

tikz_arg(arg::String) = arg
tikz_arg(arg::Pair{String,T}) where T = "$(arg[1])=$(arg[2])"

function Base.convert(::Type{MIME"text/tikz"}, args::Vector{<:TikZarg})
    if !isempty(args)
        "[" * join(map(tikz_arg, args), ",\n") * "]"
    else
        ""
    end
end

function tikz_environment(fun::Function, environment::String, args::TikZarg...)
    args = convert(MIME"text/tikz", TikZarg[args...])
    "\\begin{$(environment)}$(args)\n"*indent(fun())*"\n\\end{$(environment)}"
end

mutable struct TikZnode{T}
    label::String
    id::String
    args::Vector{TikZarg}
    x::T
    y::T
end
TikZnode(label::String, args::TikZarg...; id="", x::T=0, y::T=0) where T =
    TikZnode{T}(label, id, TikZarg[args...], x, y)

function Base.convert(::Type{MIME"text/tikz"}, node::TikZnode{T}) where T
    args = convert(MIME"text/tikz", node.args)
    id = if !isempty(node.id)
        " ($(node.id)) "
    else
        ""
    end
    pos = if node.x != zero(T) || node.y != zero(T)
        " at ($(node.x),$(node.y)) "
    else
        ""
    end
    "\\node$(args)$(id)$(pos){$(node.label)};"
end

tikz_scope(fun::Function, args::TikZarg...) = tikz_environment(fun, "scope", args...)

export indent, TikZarg, tikz_arg, tikz_environment, TikZnode, tikz_scope
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2D.jl/blob/master/LICENSE

using Documenter, MortarContact2D

makedocs(modules=[MortarContact2D],
         format = Documenter.HTML(),
         checkdocs = :all,
         sitename = "MortarContact2D.jl",
         pages = ["index.md"]
        )


text:
using FixedPointNumbers, Test

@test isempty(detect_ambiguities(FixedPointNumbers, Base, Core))

@testset "normed" begin
    include("normed.jl")
end
@testset "fixed" begin
    include("fixed.jl")
end

@testset "traits" begin
    include("traits.jl")
end


text:
function conic_form!(x::AbstractVariable, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        add_to_id_to_variables!(unique_conic_forms, x)
        if vexity(x) == ConstVexity()
            # do exactly what we would for a constant
            objective = ConicObj()
            objective[objectid(:constant)] = (vec([real(evaluate(x));]),vec([imag(evaluate(x));]))
            cache_conic_form!(unique_conic_forms, x, objective)
        else
            objective = ConicObj()
            vec_size = length(x)

            objective[x.id_hash] = (real_conic_form(x), imag_conic_form(x))
            objective[objectid(:constant)] = (spzeros(vec_size, 1), spzeros(vec_size, 1))
            # placeholder values in unique constraints prevent infinite recursion depth
            cache_conic_form!(unique_conic_forms, x, objective)
            if !(sign(x) == NoSign() || sign(x) == ComplexSign())
                conic_form!(sign(x), x, unique_conic_forms)
            end

            # apply the constraints `x` itself carries
            for constraint in constraints(x)
                conic_form!(constraint, unique_conic_forms)
            end
        end
    end
    return get_conic_form(unique_conic_forms, x)
end


text:
"""
`function construct_fcast_and_hist_dfs(m, cond_type, vars, save_to_table, table_caption,
filename, savedir, forecast_string, include_T_in_df_forecast, use_4q, hist_start,
forecast_end)`

Construct and return two DataFrame objects that respectively contain the forecast and corresponding
history series for a set of specified variables. Alternatively, if `save_to_table` is
specified to be true with `table_caption`, `filename`, and `savedir` also specified, then
`construct_fcast_and_hist_dfs` will save the DataFrames as LaTeX tables in the specified directory.

### Arguments
- `m::AbstractDSGEModel`: The model object that was used for the forecast.
- `cond_type::Symbol`: The conditioning of the forecast.
- `vars::Vector{Symbol}`: The vector of forecasted variables (can be any combination of observables or pseudo-observables) from which to construct the table.

### Keyword Arguments
- `save_to_table::Bool`: Boolean indicator of whether or not to save the table. If set to be true, then one should also specify the `table_caption`, `filename`, and `savedir`.
- `forecast_string::String`: The forecast string (if relevant) for the given forecast.
- `include_T_in_df_forecast`: Whether or not the last historical period value should be included in the forecast table for reference.
- `use_4q::Bool`: Whether to pull the 4q output_var as opposed to the standard quarterly.
- `hist_start::Date`: The date from which to start the historical table.
- `forecast_end::Date`: The date from which to end the forecast table.
"""
function construct_fcast_and_hist_dfs(m::AbstractDSGEModel, cond_type::Symbol,
                                      vars::Vector{Symbol};
                                      save_to_table::Bool = false,
                                      table_caption::String = "",
                                      filename::String = "",
                                      savedir::String = "",
                                      forecast_string::String = "",
                                      include_T_in_df_forecast::Bool = true,
                                      use_4q::Bool = false,
                                      hist_start::Date = quartertodate("2008-Q4"),
                                      forecast_end::Date = iterate_quarters(date_forecast_start(m), 12))

    # Assert that if save_to_table is true, there is a caption, filename, and a savedir
    if save_to_table
        @assert !isempty(table_caption) && !isempty(filename) && !isempty(savedir)
    end

    if use_4q
        mb_histobs      = read_mb(m, :full, cond_type, :hist4qobs, forecast_string = forecast_string)
        mb_histpseudo   = read_mb(m, :full, cond_type, :hist4qpseudo, forecast_string = forecast_string)
        mb_forecastobs      = read_mb(m, :full, cond_type, :forecast4qobs, forecast_string = forecast_string)
        mb_forecastpseudo   = read_mb(m, :full, cond_type, :forecast4qpseudo, forecast_string = forecast_string)

        mb_histobs = create_q4q4_mb(mb_histobs)
        mb_histpseudo = create_q4q4_mb(mb_histpseudo)
        mb_forecastobs = create_q4q4_mb(mb_forecastobs)
        mb_forecastpseudo = create_q4q4_mb(mb_forecastpseudo)
    else
        mb_histobs      = read_mb(m, :full, cond_type, :histobs, forecast_string = forecast_string)
        mb_histpseudo   = read_mb(m, :full, cond_type, :histpseudo, forecast_string = forecast_string)
        mb_forecastobs      = read_mb(m, :full, cond_type, :forecastobs, forecast_string = forecast_string)
        mb_forecastpseudo   = read_mb(m, :full, cond_type, :forecastpseudo, forecast_string = forecast_string)
    end

    obs = intersect(vars, m.observables.keys)
    pseudo = intersect(vars, m.pseudo_observables.keys)

    df_histobs    = mb_histobs.means[:, vcat(:date, obs)]
    df_histpseudo = mb_histpseudo.means[:, vcat(:date, pseudo)]
    hist_start_ind = findfirst(x -> x == hist_start, df_histobs[!, :date])
    df_histobs    = df_histobs[hist_start_ind:end, :]
    df_histpseudo = df_histpseudo[hist_start_ind:end, :]

    df_forecastobs      = mb_forecastobs.means[:, vcat(:date, obs)]
    df_forecastpseudo   = mb_forecastpseudo.means[:, vcat(:date, pseudo)]
    if use_4q
        # If producing 4q figures, then the forecast_end date must be a Q4 date
        forecast_end = quartertodate(string(Dates.year(forecast_end))*"-Q4")
    end
    forecast_end_ind = findfirst(x -> x == forecast_end, df_forecastobs[!, :date])
    df_forecastobs      = df_forecastobs[1:forecast_end_ind, :]
    df_forecastpseudo   = df_forecastpseudo[1:forecast_end_ind, :]

    if include_T_in_df_forecast
        df_forecastobs = append!(DataFrame(df_histobs[end, :]), df_forecastobs)
        df_forecastpseudo = append!(DataFrame(df_histpseudo[end, :]), df_forecastpseudo)
    end

    if isdefined(DataFrames, :innerjoin)
        df_forecast = innerjoin(df_forecastobs, df_forecastpseudo, on = :date)
        df_hist     = innerjoin(df_histobs, df_histpseudo, on = :date)
    else
        df_forecast = join(df_forecastobs, df_forecastpseudo, on = :date)
        df_hist     = join(df_histobs, df_histpseudo, on = :date)
    end

    obs_keys = m.observables.keys
    pseudo_keys = m.pseudo_observables.keys

    # Constructs mapping from standard key names in observable/pseudo observable means-bands objects
    # to actual names of variables and the corresponding inverse mapping (for the purposes
    # of creating a mapping of actual names of variables to units)
    header_mappings = create_table_header_mappings(m, vars)
    units = unit_mappings(m, df_forecast, header_mappings, use_4q = use_4q)

    rename!(df_forecast, header_mappings)
    rename!(df_hist, header_mappings)

    if save_to_table
        fcast_table_caption =  table_caption*" Forecast"
        hist_table_caption  =  table_caption*" History"
        fcast_filename = filename*"_forecast"
        hist_filename  = filename*"_history"

        df_to_table(df_forecast, fcast_table_caption, fcast_filename, savedir, units)
        df_to_table(df_hist, hist_table_caption, hist_filename, savedir, units)
    else
        return df_forecast, df_hist
    end
end

# Implement a df_to_table like function that splits up the df/units into
# sub-DataFrames/Dictionaries of 3/4 variables and creates a single LaTeX
# document that has all of the variables in it

# Enforce that the first column of df is a date column named :date

"""
`function df_to_table(df, caption, filename, savedir, units)`

This is the low level function that is called by `construct_fcast_and_hist_dfs` if the
`save_to_table` kwarg for that function is set to be true.

Alternatively, this function can be called directly, provided with the relevant arguments
for naming, labeling, and saving the table.

### Arguments
- `df::DataFrame`: The DataFrame object that is storing the various series.
- `caption::String`: The title of the LaTeX table.
- `filename::String`: The name of the file.
- `savedir::String`: The filepath ending in the lowest level directory that should contain the table.
- `units::OrderedDict{Symbol, String}`: A dictionary that maps the column names of `df` to
the units of that particular series (this calculation is done automatically in if
`save_to_table` is set to be true.
"""
function df_to_table(df::DataFrame, caption::String, filename::String, savedir::String,
                     units::OrderedDict{Symbol, String})

    # Open the TeX file
    savedir = savedir[end] == "/" ? savedir : savedir*"/"
    if !ispath(savedir)
        println("The savedir path provided does not currently exist. Do you want to create the path '"*savedir*"'? y/n")
        answer = readline(stdin)
        if answer == "y"
            mkpath(savedir)
        else
            error("Create the proper savedir and call df_to_table again.")
        end
    else
        mkpath(savedir)
    end
    table_out = savedir*filename*".tex"
    fid = open(table_out, "w")

    write_table_preamble(fid)
    function write_single_table(fid::IOStream, df::DataFrame, units::OrderedDict{Symbol, String})

        # Write header
        n_columns = length(propertynames(df))
        col_str   = repeat("c", n_columns)

        @printf fid "%s%s%s" "\\begin{longtable}{" col_str "}\n"
        @printf fid "\\caption{%s}\n" caption
        @printf fid "\\\\ \\hline\n"

        # Write column names
        date_range = Vector(df[:, :date])

        column_keys = propertynames(df)

        @printf fid "%s " column_keys[1]
        for i in 2:n_columns
            column_key = string(column_keys[i])
            if occursin(r"_", column_key) # if the key has an underscore then replace it with the proper LaTeX syntax
                sub_strs = split(column_key, "_")
                column_key = sub_strs[1]*"\\_"*sub_strs[2]
            end
            column_entry = "\\parbox\\{0.3\\linewidth\\}\\{\\centering "*column_key*"\\}"
            if i != n_columns
                @printf fid "& %s " column_entry
            else
                @printf fid "& %s \\\\\n" column_entry
            end
        end

        # Write units
        for i in 2:n_columns
            if i != n_columns
                @printf fid "& %s " units[column_keys[i]]
            else
                @printf fid "& %s\n" units[column_keys[i]]
            end
        end
        @printf fid "\\\\ \\hline\n"
        @printf fid "\\endhead\n"

        for (i, date) in enumerate(date_range)
            for (j, key) in enumerate(column_keys)
                if j != length(column_keys)
                    if key == :date
                        @printf fid "%s " df[:, key][i]
                    else
                        @printf fid "& %.2f " df[:, key][i]
                    end
                else
                    @printf fid "& %.2f \\\\\n" df[:, key][i]
                end
            end
        end

        @printf fid "\\end{longtable}\n"
    end

    k = 1
    for i in 1:3:length(units)
        units_keys   = i+2 < length(units) ? units.keys[i:i+2] : units.keys[i:end]
        df_subset    = DataFrame()
        units_subset = OrderedDict{Symbol, String}()

        df_subset[!, :date] = df[!, :date]
        for unit in units_keys
            df_subset[!, unit] = df[!, unit]
            units_subset[unit] = units[unit]
        end
        write_single_table(fid, df_subset, units_subset)
        if k % 2 == 0
            @printf fid "\\clearpage\n" # every two tables, break the page
        else
            @printf fid "\\vspace*{.5cm}\n"
        end
        k += 1
    end

    @printf fid "\\end{document}"

    close(fid)
end

# Rename keys in the obs dictionaries
# So that the DataFrame has LaTeX conformant names
function create_table_header_mappings(m::AbstractDSGEModel, vars::Vector{Symbol})
    obs_keys = m.observables.keys
    header_mappings = OrderedDict{Symbol, Symbol}()
    for var in vars
        if var in obs_keys
            header_mappings[var] = Symbol(m.observable_mappings[var].name)
        else
            header_mappings[var] = detexify(Symbol(m.pseudo_observable_mappings[var].name))
        end
    end
    return header_mappings
end

# Defining the units for the variables included
function unit_mappings(m::AbstractDSGEModel, df::DataFrame,
                       header_mappings::OrderedDict{Symbol, Symbol}; use_4q::Bool = false)
    units = OrderedDict{Symbol, String}()
    quarter = use_4q ? "Q4" : "Q"
    obs_keys = m.observables.keys
    pseudo_keys = m.pseudo_observables.keys
    for key in propertynames(df)
        name = key != :date ? header_mappings[key] : continue
        if key in obs_keys
            if occursin(r"pct_annualized", string(m.observable_mappings[key].rev_transform))
                units[name] = "("*quarter*"/"*quarter*") \\% Annualized"
            elseif occursin(r"quartertoannual", string(m.observable_mappings[key].rev_transform))
                units[name] = quarter
            end
        elseif key in pseudo_keys
            if occursin(r"pct_annualized", string(m.pseudo_observable_mappings[key].rev_transform))
                units[name] = "("*quarter*"/"*quarter*") \\% Annualized"
            elseif occursin(r"quartertoannual", string(m.pseudo_observable_mappings[key].rev_transform))
                units[name] = quarter
            elseif occursin(r"identity", string(m.pseudo_observable_mappings[key].rev_transform))
                units[name] = quarter
            end
        end
    end
    return units
end


text:
#
# lists.jl -
#
# Management of Tcl lists of objects.
#

# Let Tcl list be iterable.
@static if VERSION < v"0.7"

    function Base.start(iter::TclObj{List})
        objc, objv = __getlistelements(__objptr(iter))
        return (0, objc, objv)
    end

    function Base.done(iter::TclObj{List}, state)
        state[1] ≥ state[2]
    end

    function Base.next(iter::TclObj{List}, state)
        i, n, objv = state
        i += 1
        item = __objptr_to(Any, unsafe_load(objv, i))
        return item, (i, n, objv)
    end

else

    function Base.iterate(iter::TclObj{List})
        objc, objv = __getlistelements(__objptr(iter))
        return iterate(iter, (0, objc, objv))
    end

    function Base.iterate(iter::TclObj{List}, state)
        i, objc, objv = state
        if i < objc
            i += 1
            item = __objptr_to(Any, unsafe_load(objv, i))
            return item, (i, objc, objv)
        else
            return nothing
        end
    end

end

# Let Tcl list be indexable.

Base.length(list::TclObj{List}) = llength(list)

Base.lastindex(list::TclObj{List}) = llength(list)

Base.getindex(list::TclObj{List}, i::Integer) = lindex(list, i)

function Base.getindex(list::TclObj{List}, msk::AbstractVector{Bool})
    n = count(!iszero, msk)
    v = Array{Any,1}(undef, n)
    if n < 1
        return v
    end
    len = length(list)
    i = 0
    for j in 1:length(msk)
        if msk[j]
            i += 1
            v[i] = (j ≤ len ? list[j] : nothing)
        end
    end
    return standardizetype(v)
end

Base.getindex(list::TclObj{List}, msk::AbstractArray{Bool}) =
    error("indexing a list by a multi-dimensional array of booleans is not possible")

function Base.getindex(list::TclObj{List}, J::AbstractArray{<:Integer})
    A = similar(J, Any)
    if length(J) < 1
        return A
    end
    len = length(list)
    for i in eachindex(A, J)
        j = J[i]
        A[i] = (1 ≤ j ≤ len ? list[j] : nothing)
    end
    return standardizetype(A)
end

# FIXME: call Tcl_ListObjReplace to implement Base.setindex!
#
# function Base.setindex!(list::TclObj{List}, value, i::Integer)
# end
#
# function Base.setindex!(list::TclObj{List}, value, r::UnitRange)
# end

#
# Implement lists objects (see objects.jl).
#
#     Iterables like vectors and tuples yield lists.  No arguments, yield empty
#     lists.
#

TclObj() = TclObj{List}(__newlistobj())

__objptr() = __newobj()

__newobj() = __newlistobj()

TclObj(itr::Iterables) = TclObj{List}(__newobj(itr))

__objptr(itr::Iterables) = __newobj(itr)

function __newobj(itr::Iterables) ::TclObjPtr
    listptr = __newlistobj()
    try
        for val in itr
            __lappend(listptr, val)
        end
    catch ex
        Tcl_DecrRefCount(listptr)
        rethrow(ex)
    end
    return listptr
end

"""
```julia
__newlistobj(itr)
```

yields a pointer to a new Tcl list object whose items are taken from
the iterable collection `itr`.

```julia
__newlistobj(args...; kwds...)
```

yields a pointer to a new Tcl list object whose leading items are taken from
`args...` and to which are appended the `(key,val)` pairs from `kwds...` so as
to mimic Tk options.

Beware that the returned object is not managed and has a zero reference count.
The caller is reponsible of taking care of that.

"""
__newlistobj() = Tcl_NewListObj(0, Ptr{TclObjPtr}(0))

function __newlistobj(args...; kwds...) ::TclObjPtr
    listptr = __newlistobj()
    try
        @__build_list listptr args kwds
    catch ex
        Tcl_DecrRefCount(listptr)
        rethrow(ex)
    end
    return listptr
end

@inline __objptr_to(::Type{Vector}, listptr::Ptr{Cvoid}) =
    __buildvector(__getlistelements(listptr)...)

@inline __buildvector(objc::Integer, objv::Ptr{TclObjPtr}) =
    buildvector(i -> __objptr_to(Any, unsafe_load(objv, i)), objc)

"""
```julia
buildvector(f, n)
```
"""
function buildvector(f::Function, n::Integer) :: Vector
    v = Vector{Any}(undef, n)
    if n ≥ 1
        for i in 1:n
            v[i] = f(i)
        end
        return standardizetype(v)
    end
    return v
end

"""
```julia
standardizetype(A)
```

if all elements of array `A` can be promoted to the same type `T`, returns `A`
converted as a `Array{T}`; otherwise, returns `A` unchanged.  Promotions rules
are a bit different than the ones in Julia, in the sense that the family type
of elements is preserved.  Families are: strings, integers and floats.

"""
function standardizetype(A::AbstractArray{Any,N}) where N
    if length(A) ≥ 1
        T = Ref{DataType}()
        first = true
        for i in eachindex(A)
            T[] = (first ? typeof(A[i]) :
                   __promote_elem_type(T[], typeof(A[i])))
            first = false
        end
        if T[] != Any
            # A common type has been found, promote the vector to this common
            # type.
            return convert(Array{T[],N}, A)
        end
    end
    return A
end

# Rules for combining list element types and find a more precise common type
# than just `Any`.  Combinations of integers are promoted to the largest
# integer type and similarly for floats but mixture of floats and integers
# yield `Any`.

__promote_elem_type(::DataType, ::DataType) = Any

for T in (Integer, AbstractFloat)
    @eval begin

        function __promote_elem_type(::Type{T1},
                                     ::Type{T2}) where {T1<:$T,T2<:$T}
            return promote_type(T1, T2)
        end

        function __promote_elem_type(::Type{Vector{T1}},
                                     ::Type{Vector{T2}}) where {T1<:$T,
                                                                T2<:$T}
            return Vector{promote_type(T1, T2)}
        end

    end
end

__promote_elem_type(::Type{String}, ::Type{String}) = String

__promote_elem_type(::Type{Vector{String}}, ::Type{Vector{String}}) =
    Vector{String}

#------------------------------------------------------------------------------

"""
```julia
list([interp,] args...; kwds...)
```

yields a list of Tcl objects consisting of the one object per argument
`args...` (in the same order as they appear) and then followed by two objects
per keyword, say `key=val`, in the form `-key`, `val` (note the hyphen in front
of the keyword name).  To allow for option names that are Julia keywords, a
leading underscore is stripped, if any, in `key`.

Optional argument `interp` is the Tcl interpreter to use for error messages
and callbacks.

Lists are iterable and indexable, as illustrated by the following examples:

``julia
lst = Tcl.list(π,1,"hello",2:6)
length(lst) # -> the number of items in the list
lst[1]      # -> 3.1415...
lst[end]    # -> [2,3,4,5,6]
lst[2:3]    # -> Any[1,"hello"]
lst[0]      # -> nothing
lst[end+1]  # -> nothing
for itm in lst
    println(itm)
end
sel = map(i -> isa(i, Number), lst) # -> [true,true,false,false]
lst[sel] # -> Any[3.14159,1]
```

You may note that, (i) like Tcl lists, getting an out of bound list item just
yields nothing; (ii) lists are retrieved as Julia arrays with, if possible,
homogeneous element type (otherwise `Any`).

You may sub-select list elements.  For instance to extract the numbers of a
list:

``julia
lst = Tcl.list(π,1,"hello",2:6)
sel = map(i -> isa(i, Number), lst) # -> [true,true,false,false]
lst[sel] # -> Any[3.14159,1]
```

Use `push!` (or [`Tcl.lappend!`](@ref]) to append elements to a list.  Use
[`Tcl.concat`](@ref) to concatenate lists.

See also: [`Tcl.concat`](@ref), [`Tcl.lindex`](@ref), [`Tcl.lappend!`](@ref),
          [`Tcl.exec`](@ref), .

"""
function list(interp::TclInterp, args...; kwds...)
    listptr = C_NULL
    __set_context(interp)
    try
        listptr = Tcl_IncrRefCount(__newlistobj())
        @__build_list listptr args kwds
        return TclObj{List}(listptr)
    finally
        listptr == C_NULL || Tcl_DecrRefCount(listptr)
         __reset_context()
    end
end

list(args...; kwds...) =
    list(getinterp(), args...; kwds...)

"""
```julia
Tcl.llength(lst)
```
or
```julia
length(lst)
```

yeild the number of element of the Tcl list `lst`.

See also: [`Tcl.list`](@ref).

"""
function llength(lst::TclObj{List}) :: Int
    status, length = Tcl_ListObjLength(C_NULL, __objptr(lst))
    status == TCL_OK || Tcl.error("failed to query length of list")
    return length
end

"""
```julia
Tcl.lappend!([interp,] lst, args...; kwds...)
```
or
```julia
push!(lst, args...; kwds...)
```

append to the list `lst` of Tcl objects one object per argument `args...` (in
the same order as they appear) and then followed by two objects per keyword,
say `key=val`, in the form `-key`, `val` (note the hyphen in front of the
keyword name).  To allow for option names that are Julia keywords, a leading
underscore is stripped, if any, in `key`; for instance:

```julia
Tcl.lappend!(lst, _in="something")
```

appends `"-in"` and `something` to the list `lst`.

Optional argument `interp` is the Tcl interpreter to use for error messages
and callbacks.

See also: [`Tcl.list`](@ref).

"""
function lappend!(interp::TclInterp, list::TclObj{List}, args...; kwds...)
    __set_context(interp)
    try
        listptr = __objptr(list)
        @__build_list listptr args kwds
    finally
        __reset_context()
    end
    return list
end

lappend!(list::TclObj{List}, args...; kwds...) =
    lappend!(getinterp(), list, args...; kwds...)

Base.push!(list::TclObj{List}, args...; kwds...) =
    lappend!(list, args...; kwds...)

Base.push!(interp::TclInterp, list::TclObj{List}, args...; kwds...) =
    lappend!(interp, list, args...; kwds...)

# Appending a new item to a list with Tcl_ListObjAppendElement increments the
# reference count of the item, this is the only side effect for the item.  That
# is to say, the appended item is not be duplicated, just shared.  So, to
# manage the memory associated with the item, we can increment its reference
# count before appending and decrement it after without effects on the
# performances.  This is not necessary for a managed object.

function __lappend(listptr::TclObjPtr, item)
    objptr = Tcl_IncrRefCount(__newobj(item))
    status = Tcl_ListObjAppendElement(__intptr(), listptr, objptr)
    Tcl_DecrRefCount(objptr)
    status == TCL_OK || __lappend_error()
end

function __lappend(listptr::TclObjPtr, obj::ManagedObject)
    if Tcl_ListObjAppendElement(__intptr(), listptr, __objptr(obj)) != TCL_OK
        __lappend_error()
    end
end

function __lappend(listptr::TclObjPtr, func::Function)
    intptr = __intptr()
    cmd = __createcallback(intptr, f)
    if Tcl_ListObjAppendElement(intptr, listptr, __objptr(cmd)) != TCL_OK
        __lappend_error()
    end
end

function __lappend(listptr::TclObjPtr, cmd::Callback)
    intptr = __intptr()
    if intptr != C_NULL && cmd.intptr == intptr
        # Create the command in the interpreter.
        cmd = __createcallback(intptr, string(cmd.obj), cmd.func)
    end
    if Tcl_ListObjAppendElement(intptr, listptr, __objptr(cmd)) != TCL_OK
        __lappend_error()
    end
end

__lappend_error() =
    __contextual_error("failed to append a new item to the Tcl list")

function __lappendoption(listptr::TclObjPtr, key::String, val)
    # First, append the key.
    intptr = __intptr()
    option = "-"*(length(key) ≥ 1 && key[1] == '_' ? key[2:end] : key)
    objptr = Tcl_IncrRefCount(__newobj(option))
    status = Tcl_ListObjAppendElement(intptr, listptr, objptr)
    Tcl_DecrRefCount(objptr)
    if status == TCL_OK
        # Second, append the value.
        objptr = Tcl_IncrRefCount(__objptr(val))
        status = Tcl_ListObjAppendElement(intptr, listptr, objptr)
        Tcl_DecrRefCount(objptr)
    end
    if status != TCL_OK
        Tcl.error("failed to append a new option to the Tcl list")
    end
    nothing
end

__lappendoption(listptr::TclObjPtr, key::Symbol, val) =
    __lappendoption(listptr, string(key), val)

"""
```julia
Tcl.concat([interp,]args...)
```

concatenates the specified arguments and yields a Tcl list.  Compared to
`Tcl.list` which considers that each argument correspond to a single item,
`Tcl.concat` flatten its arguments and does not accept keyword arguments.

Optional argument `interp` is the Tcl interpreter to use for error reporting
or to create callbacks.

See also: [`Tcl.list`](@ref), [`Tcl.eval`](@ref).

"""
function concat(interp::TclInterp, args...)
    listptr = C_NULL
    __set_context(interp)
    try
        listptr = Tcl_IncrRefCount(__newlistobj())
        @__concat_args listptr args
        return TclObj{List}(listptr)
    finally
        listptr == C_NULL || Tcl_DecrRefCount(listptr)
        __reset_context()
    end
end

concat(args...) = concat(getinterp(), args...)

# The basic functions used by most Tcl list manipulation functions are
# Tcl_ListObjGetElements, Tcl_ListObjReplace and Tcl_ListObjAppendElement.
# Tcl_ListObjReplace does not call Tcl_ListObjAppendElement.
#
# As a general rule, modifying a shared list is not allowed.  Thus
# Tcl_ListObjReplace and Tcl_ListObjAppendElement must not be applied to a
# shared list object.  This limits the risk of building circular lists.
#

# For atomic objects which are considered as single list element, __concat is
# equivalent to __lappend.
__concat(listptr::TclObjPtr, item::T) where T =
    __concat(listptr, AtomicType(T), item)

__concat(listptr::TclObjPtr, ::Atomic, item) =
    __lappend(listptr, item)

# Strings are iterables but we want that making a list out of string(s) yields
# a single element per string (not per character) so we have to short-circuit
# __concat(listptr, itr).  Note that `Number` are perfectly usable as iterables
# but we add them to the union below in order to use a faster method for them.
function __concat(listptr::TclObjPtr, ::NonAtomic, str::AbstractString)
    objptr = Tcl_IncrRefCount(__newobj(str))
    status = Tcl_ListObjAppendList(__intptr(), listptr, objptr)
    Tcl_DecrRefCount(objptr)
    if status != TCL_OK
        __contextual_error("failed to concatenate a string to a Tcl list")
    end
end

function __concat(listptr::TclObjPtr, ::NonAtomic, obj::ManagedObject)
    if Tcl_ListObjAppendList(__intptr(), listptr, __objptr(obj)) != TCL_OK
        __contextual_error("failed to concatenate Tcl lists")
    end
end

# Everything else is assumed to be an iterable.
function __concat(listptr::TclObjPtr, ::NonAtomic, itr) ::TclObjPtr
    for val in itr
        __concat(listptr, val)
    end
end

"""
```julia
Tcl.lindex([T,] [interp,] list, i)
```

yields the element at index `i` in Tcl list `list`.  An *empty* result is
returned if index is out of range.

If optional argument `T` is omitted, the type of the returned value reflects
that of the Tcl variable; otherwise, `T` can be `String` to get the string
representation of the value or `TclObj` to get a managed Tcl object.  The
latter type is more efficient if the returned item is intended to be put in a
Tcl list or to be an argument of a Tcl script or command.

Tcl interpreter `interp` may be provided to have more detailed error messages
in case of failure.

See also: [`Tcl.list`](@ref), [`Tcl.getvar`](@ref).

"""
lindex(list::TclObj{List}, i::Integer) =
    lindex(Any, list, i)

lindex(::Type{T}, list::TclObj{List}, i::Integer) where {T} =
    __itemptr_to(T, C_NULL, __lindex(list, i))

lindex(interp::TclInterp, list::TclObj{List}, i::Integer) =
    lindex(Any, interp, list, i)

lindex(::Type{T}, interp::TclInterp, list::TclObj{List}, i::Integer) where {T} =
    __itemptr_to(T, interp.ptr, __lindex(interp, list, i))

__itemptr_to(::Type{T}, intptr::TclInterpPtr, objptr::TclObjPtr) where {T} =
    (objptr == C_NULL ? __missing_item(T) : __objptr_to(T, objptr))

"""
```julia
___missing_item(T)
```

yields the value of missing list item of type `T`.  May throw an error if
missing items of such type are not allowed.

"""
__missing_item(::Type{String}) = ""
__missing_item(::Type{Any}) = nothing
__itemptr_item(::Type{TclObj}) = TclObj()
__itemptr_item(::Type{Vector}) = Any[]
__missing_item(::Type{T}) where {T<:Union{Integer,AbstractFloat}} = zero(T)
__missing_item(::Type{Char}) = '\0'

# Get a list item.
#
#     The convention of Tcl_ListObjIndex is to return TCL_ERROR if some error
#     occured and TCL_OK with a NULL pointer if index is out of range.

function __lindex(list::TclObj{List}, i::Integer)
    status, objptr = Tcl_ListObjIndex(C_NULL, __objptr(list), i)
    if status != TCL_OK
        Tcl.error("failed to get Tcl list element at index $i")
    end
    return objptr
end

function __lindex(interp::TclInterp, list::TclObj{List}, i::Integer)
    status, objptr = Tcl_ListObjIndex(interp.ptr, __objptr(list), i)
    if status != TCL_OK
        Tcl.error(interp)
    end
    return objptr
end

# Yields (objc, objv) do not free this buffer (see Tcl doc.)
function __getlistelements(listptr::TclObjPtr)
    status, objc, objv = Tcl_ListObjGetElements(__intptr(), listptr)
    if status != TCL_OK
        __contextual_error("failed to convert Tcl object into a list")
    end
    return objc, objv
end


text:
using DynamicHMC
using DiffWrappers
using ContinuousTransformations
using StochasticVolatility
import Distributions: Uniform, InverseGamma
using JLD

RNG = Base.Random.GLOBAL_RNG

ρ = 0.8
σ = 0.6
y = simulate_stochastic(ρ, σ, 10000)

# set up the model
model = StochasticVolatility(StochasticVolatility(y, Uniform(-1, 1), InverseGamma(1, 1), 10000)...)

# we start the estimation process from the true values
θ₀ = inverse(model.transformation, (ρ, σ))

# wrap for gradient calculations

fgw = ForwardGradientWrapper(model, θ₀)

# sample

sample, tuned_sampler = NUTS_tune_and_mcmc(RNG, fgw, 5000; q = θ₀)

# posterior analysis

posterior = variable_matrix(sample)

sample_ρ, sample_σ = map_by_row(model.transformation, posterior)

# save

save(StochasticVolatility.path("results", "results.jld"), Dict("sample" => sample, "posterior" => posterior, "sample_ρ" => sample_ρ, "sample_σ" => sample_σ, "ρ" => ρ, "σ" => σ))


text:
module Hyperspecialize

using InteractiveUtils
using Base.Iterators

export @concretize, @concretization, @widen, @replicable

postwalk(f, x) = f(x)
postwalk(f, x::Expr) = f(Expr(x.head, map(arg->postwalk(f, arg), x.args)...))

macro isdefined(var)
 quote
   try
     local _ = $(esc(var))
     true
   catch err
     isa(err, UndefVarError) ? false : rethrow(err)
   end
 end
end

"""
    concretesubtypes(t)

Return an `Array` containing all concrete subtypes of `t` at load time.

# Examples
```julia-repl
julia> Hyperspecialize.concretesubtypes(Real)
16-element Array{Any,1}:
 BigFloat
 Float16
 Float32
 Float64
 Bool
 BigInt
 Int128
 Int16
 Int32
 Int64
 Int8
 UInt128
 UInt16
 UInt32
 UInt64
 UInt8
```
"""
function concretesubtypes(t)
  if isconcretetype(t)
    return [t]
  else
    return vcat([concretesubtypes(s) for s in subtypes(t)]...)
  end
end

"""
    allsubtypes(t)

Return an `Array` containing all subtypes of `t` at load time.

# Examples
```julia-repl
julia> Hyperspecialize.allsubtypes(Real)
24-element Array{Type,1}:
 Real
 AbstractFloat
 BigFloat
 Float16
 Float32
 Float64
 AbstractIrrational
 Irrational
 Integer
 Bool
 Signed
 BigInt
 Int128
 Int16
 Int32
 Int64
 Int8
 Unsigned
 UInt128
 UInt16
 UInt32
 UInt64
 UInt8
 Rational
```
"""
function allsubtypes(t)
  return vcat([t], [allsubtypes(s) for s in subtypes(t)]...)
end

function parse_element(base_mod, T)
  if T isa Expr && T.head == :. && length(T.args) == 2
    (M, K) = T.args
  else
    M = base_mod
    K = QuoteNode(T)
  end
  return (M, K)
end

struct Replicable
  E::Any
  defined::Set{Tuple{Vararg{Type}}}
  type_tags::Vector{Tuple{Module, Symbol}}
end

struct Concretization
  base::Set{Type}
  contextual::Dict{Tuple{Module, Int}, Set{Type}}
  replicable_tags::Vector{Tuple{Module, Int}}
  Concretization(types::Set{Type}) = new(types, Dict{Tuple{Module, Int}, Set{Type}}(), Vector{Tuple{Module, Int}}())
end

struct State
  replicables::Vector{Replicable}
  concretizations::Dict{Symbol, Concretization}
  State() = new(Vector{Replicable}(), Dict{Symbol, Concretization}())
end

function _concretize(base_mod::Module, target_mod::Module, key::Symbol, types::Type)
  return _concretize(base_mod, target_mod, key, [types])
end

function _concretize(base_mod::Module, target_mod::Module, key::Symbol, types)
  return _concretize(base_mod, target_mod, key, Set{Type}(types))
end

function _concretize(base_mod::Module, target_mod::Module, key::Symbol, types::Set{Type})
  if base_mod == target_mod
    if !isdefined(base_mod, :__hyperspecialize__)
      Core.eval(base_mod, quote
        const global __hyperspecialize__ = Hyperspecialize.State()
      end)
    end
    if haskey(target_mod.__hyperspecialize__.concretizations, key)
      error("cannot reconcretize \"$key\" in module \"$target_mod\"")
    else
      target_mod.__hyperspecialize__.concretizations[key] = Concretization(types)
    end
  else
    error("cannot concretize \"$key\" in module \"$target_mod\" from module \"$base_mod\"")
  end
  return Set{Type}(target_mod.__hyperspecialize__.concretizations[key].base)
end

"""
    @concretize(typetag, types)

Define the set of types corresponding to a type tag as `types`, where `types`
is either a single type or any collection that may be passed to a set
constructor.  A type tag is a module-qualified symbol `mod.:Key` where mod
specifies a module and `:Key` is a symbol.  If just the `:Key` is given, then
the module is assumed to be the module in which the macro was expanded.

Note that you may not concretize a type in another module.

# Examples
```julia-repl
julia> @concretize Main.BestInts [Int32, Int64]
Set(Type[Int32, Int64])

julia> @concretize BestFloats Float64
Set(Type[Float64])

julia> @concretize BestStrings (String,)
Set(Type[String])

julia> @concretization BestInts
Set(Type[Int32, Int64])
```
"""
macro concretize(K, T)
  (M, K) = parse_element(__module__, K)
  return :(_concretize($(esc(__module__)), $(esc(M)), $(K), $(esc(T))))
end

function _concretize(base_mod::Module, target_mod::Module, key::Symbol)
  if _concretization(base_mod, target_mod, key) == nothing
    if isdefined(target_mod, key)
      if Core.eval(target_mod, key) isa Type
        _concretize(base_mod, target_mod, key, concretesubtypes(Core.eval(target_mod, key)))
      else
        error("cannot create default concretization from type tag $target_mod.$key: not a type")
      end
    else
      error("cannot create default concretization from type tag $target_mod.$key: not defined")
    end
  end
  return _concretization(base_mod, target_mod, key)
end

"""
    @concretize(typetag)

If no concretization exists, define the set of types corresponding to `typetag`
as the concrete subtypes of whatever type shares the name of `Key` at load
time.  `typetag` is a type tag, or a module-qualified symbol `mod.:Key` where
mod specifies a module and `:Key` is a symbol.  If just the `:Key` is given,
then the module is assumed to be the module in which the macro was expanded.

Note that you may not concretize a type tag in another module.

# Examples
```julia-repl
julia> @concretization(Main.Real)
nothing

julia> @concretize(Main.Real)
Set(Type[BigInt, Bool, UInt32, Float64, Float32, Int64, Int128, Float16, UInt128, UInt8, UInt16, BigFloat, Int8, UInt64, Int16, Int32])

julia> @concretize BestInts [Int32, Int64]
Set(Type[Int32, Int64])

julia> @concretization BestInts
Set(Type[Int32, Int64])

julia> @concretize NotDefinedHere
ERROR: cannot create default concretization from type tag Main.NotDefinedHere: not defined.
```
"""
macro concretize(K)
  (M, K) = parse_element(__module__, K)
  return :(_concretize($(esc(__module__)), $(esc(M)), $(K)))
end

function _concretization(base_mod::Module, target_mod::Module, key::Symbol)
  if !isdefined(target_mod, :__hyperspecialize__) || !haskey(target_mod.__hyperspecialize__.concretizations, key)
    return nothing
  end
  return Set{Type}(target_mod.__hyperspecialize__.concretizations[key].base)
end

"""
    @concretization(typetag)

Return the set of types corresponding to a type tag.  A type tag is a
module-qualified symbol `mod.:Key` where mod specifies a module and `:Key` is a
symbol.  If just the `:Key` is given, then the module is assumed to be the
module in which the macro was expanded.  If no concretization exists, return
nothing.

A concretization can be set and modified with `@concretize` and `@widen`

# Examples
```julia-repl
julia> @concretization(BestInts)
nothing

julia> @concretize BestInts [Int32, Int64]
Set(Type[Int32, Int64])

julia> @concretization BestInts
Set(Type[Int32, Int64])
```
"""
macro concretization(K)
  (M, K) = parse_element(__module__, K)
  return :(_concretization($(esc(__module__)), $(esc(M)), $(K)))
end

function _widen(base_mod::Module, target_mod::Module, key::Symbol, types::Type)
  return _widen(base_mod, target_mod, key, [types])
end

function _widen(base_mod::Module, target_mod::Module, key::Symbol, types)
  return _widen(base_mod, target_mod, key, Set{Type}(types))
end

function _widen(base_mod::Module, target_mod::Module, key::Symbol, types::Set{Type})
  _concretize(base_mod, target_mod, key)
  union!(target_mod.__hyperspecialize__.concretizations[key].base, types)
  map(((target_mod, num),) -> _define(base_mod, target_mod, target_mod.__hyperspecialize__.replicables[num]), target_mod.__hyperspecialize__.concretizations[key].replicable_tags)
  return Set{Type}(target_mod.__hyperspecialize__.concretizations[key].base)
end

"""
    @widen(typetag, types)

Expand the set of types corresponding to a type tag to include `types`, where
`types` is either a single type or any collection that may be passed to a set
constructor.  A type tag is a module-qualified symbol `mod.:Key` where mod
specifies a module and `:Key` is a symbol.  If just the `:Key` is given, then
the module is assumed to be the module in which the macro was expanded.  If no
concretization exists, create a default concretization consisting of the
conrete subtypes of whatever type shares the name of `Key` at load time.

If `@widen` is called for a type tag which has been referenced by a
`@replicable` code block, then that code block will be replicated even more to
reflect the new concretization.

# Examples
```julia-repl
julia> @concretize BestInts [Int32, Int64]
Set(Type[Int32, Int64])

julia> @replicable println(@hyperspecialize(BestInts))
Int32
Int64

julia> @widen BestInts (Bool, Int32, UInt128)
Bool
UInt128
Set(Type[Bool, UInt128, Int32, Int64])

julia> @concretization BestInts
Set(Type[Bool, Int8, Int32, Int64, UInt128])
```
"""
macro widen(K, T)
  (M, K) = parse_element(__module__, K)
  return :(_widen($(esc(__module__)), $(esc(M)), $(K), $(esc(T))))
end

_is_hyperspecialize(X) = false
function _is_hyperspecialize(X::Expr)
  return X.head == :macrocall &&
    length(X.args) == 3 &&
    X.args[1] == Symbol("@hyperspecialize")
end

function _get_hyperspecialize(X)
  @assert _is_hyperspecialize(X)
  return X.args[3]
end

function _define(base_mod::Module, target_mod::Module, replicable::Replicable)
  for types in product(map(type_tag -> _concretization(base_mod, type_tag...), replicable.type_tags)...)
    if !(types in replicable.defined)
      Core.eval(target_mod, postwalk(X -> begin
        if _is_hyperspecialize(X)
          types[_get_hyperspecialize(X)]
        else
          X
        end
      end, replicable.E))
      push!(replicable.defined, types)
    end
  end
end

function _replicable(base_mod::Module, E, type_tags::Tuple{Module, Symbol}...)
  if !isdefined(base_mod, :__hyperspecialize__)
    Core.eval(base_mod, quote
      const global __hyperspecialize__ = Hyperspecialize.State()
    end)
  end
  replicable = Replicable(E, Set{Any}(), [type_tags...])
  push!(base_mod.__hyperspecialize__.replicables, replicable)
  num = length(base_mod.__hyperspecialize__.replicables)
  for (target_mod, key) in type_tags
    _concretize(base_mod, target_mod, key)
    push!(target_mod.__hyperspecialize__.concretizations[key].replicable_tags, (base_mod, num))
  end
  _define(base_mod, base_mod, replicable)
end

"""
    @replicable block

Replicate the code in `block` where each type tag referred to by
`@hyperspecialize(typetag)` is replaced by an element in the concretization of
`typetag`.  `block` is replicated at global scope in the module where
`@replicable` was expanded once for each combination of types in the
concretization of each `typetag`.  A type tag is a module-qualified symbol
`mod.:Key` where mod specifies a module and `:Key` is a symbol.  If just the
`:Key` is given, then the module is assumed to be the module in which the macro
was expanded.  If no concretization exists for a type tag, create a default
concretization consisting of the conrete subtypes of whatever type shares the
name of `Key` at load time.

If `@widen` is called for a type tag which has been referenced by a
`@replicable` code block, then that code block will be replicated even more to
reflect the new concretization.

# Examples
```julia-repl
julia> @concretize BestInts [Int32, Int64]
Set(Type[Int32, Int64])

julia> @replicable println(@hyperspecialize(BestInts), @hyperspecialize(BestInts))
Int32Int32
Int32Int64
Int64Int32
Int64Int64

julia> @widen BestInts (Bool,)
BoolBool
BoolInt32
BoolInt64
Int32Bool
Int64Bool
Set(Type[Bool, Int32, Int64])
```
"""
macro replicable(E)
  elements = []
  count = 0
  E = postwalk(X -> begin
    if _is_hyperspecialize(X)
      (M, K) = parse_element(__module__, _get_hyperspecialize(X))
      push!(elements, :(($(esc(M)), $(K))))
      count += 1
      :(@hyperspecialize($count))
    else
      X
    end
  end, E)
  return :(_replicable($(esc(__module__)), $(QuoteNode(E)), $(elements...)))
end

end # module


text:
# Lightweight Representative Agent Version of HetDSGEGovDebt
mutable struct RepDSGEGovDebt{T} <: AbstractRepModel{T}
    parameters::ParameterVector{T}                # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}              # model steady-state values

    grids::OrderedDict{Symbol,Union{Grid, Array}} # Human-readable names for all the model
    keys::OrderedDict{Symbol,Int}                 # parameters and steady-states

    state_variables::Vector{Symbol}               # Vector of symbols of the state variables
    jump_variables::Vector{Symbol}                # Vector of symbols of the jump variables

    # Vector of unnormalized ranges of indices
    endogenous_states::OrderedDict{Symbol,Int}
    exogenous_shocks::OrderedDict{Symbol,Int}
    expected_shocks::OrderedDict{Symbol,Int}
    equilibrium_conditions::OrderedDict{Symbol,Int}
    endogenous_states_augmented::OrderedDict{Symbol, Int}
    observables::OrderedDict{Symbol,Int}

    spec::String                                  # Model specification number (eg "m990")
    subspec::String                               # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                # Settings/flags for computation
    testing::Bool                                 # Whether we are in testing mode or not
    observable_mappings::OrderedDict{Symbol, Observable}
end

description(m::RepDSGEGovDebt) = "RepDSGEGovDebt, $(m.subspec)"

function RepDSGEGovDebt(het::HetDSGEGovDebt;
                        custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}())

    @assert !isnan(het[:βstar].value) "Must have already called steadystate! on `het` before instantiating a RepDSGEGovDebt model"

    # Model-specific specifications
    spec               = "rep_dsge"
    subspec            = deepcopy(het.subspec)
    settings           = deepcopy(het.settings)

    # initialize empty model
    m = RepDSGEGovDebt{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(),
            Vector{Float64}(),

            # grids and keys
            OrderedDict{Symbol,Union{Grid, Array}}(), OrderedDict{Symbol,Int}(),

            # state_inds, jump_inds
            Vector{Symbol}(), Vector{Symbol}(),

            # model indices
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),

            spec,
            subspec,
            settings,
            deepcopy(het.testing),
            deepcopy(het.observable_mappings))

    # Set settings
    model_settings!(m, het)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Initialize model indices
    init_model_indices!(m)

    # Initialize parameters
    init_parameters!(m, het)

    # Initialize grids
    init_grids!(m, het)

    # initialize subspec from het
    init_subspec!(m, het)

    # Solve for the steady state
    steadystate!(m, het, same_beta = get_setting(m, :use_same_beta_as_hank))

    return m
end

function model_settings!(m::RepDSGEGovDebt, het::HetDSGEGovDebt)
    for (key, setting) in deepcopy(het.settings)
        m <= setting
    end

    # Total grid x*s
    # Not sure what to do with the below line. Leave commented for now
    # m <= Setting(:n, (get_setting(m, :nx1) +get_setting(m, :nx2)),
                 # "Total grid size, multiplying across grid dimensions.")
    # The rest of these lines were taken directly from the rank_irfs.jl file
    m <= Setting(:n_jumps, 15, "num scalar jumps")
    m <= Setting(:n_backward_looking_states, 14, "num scalar states")
    m <= Setting(:n_model_states, get_setting(m, :n_backward_looking_states) +
                 get_setting(m, :n_jumps),
                 "Number of 'states' in the state space model. Because backward and forward
                 looking variables need to be explicitly tracked for the Klein solution
                 method, we have n_states and n_jumps")
    m <= Setting(:n_model_states_augmented, get_setting(m, :n_model_states) +
                 length([:c_t1]),
                 "The number of 'states' in the augmented state space model.")

    m <= Setting(:use_same_beta_as_hank, false)
end

function init_model_indices!(m::RepDSGEGovDebt)

    # Endogenous states
    states = collect([:k′_t, :R′_t1, :i′_t1, :y′_t1, :w′_t1, :I′_t1, :bg′_t,
                      :b′_t, :g′_t, :z′_t, :μ′_t, :λ_w′_t, :λ_f′_t, :rm′_t])

    jumps = collect([:R′_t, :i′_t, :t′_t, :w′_t, :L′_t, :π′_t, :π_w′_t, :margutil′_t, :y′_t, :I′_t,
                     :mc′_t, :Q′_t, :capreturn′_t, :l′_t, :tg′_t])

    endogenous_states = collect(vcat(states, jumps))

    # Exogenous shocks
    exogenous_shocks = collect([:b_sh,:g_sh,:z_sh,:μ_sh,:λ_w_sh,:λ_f_sh,:rm_sh])

    # Equilibrium conditions
    equilibrium_conditions = collect([:eq_euler, :eq_market_clearing, :eq_lambda, :eq_transfers,
                                      :eq_investment, :eq_tobin_q, :eq_capital_accumulation,
                                      :eq_wage_phillips, :eq_price_phillips, :eq_marginal_cost,
                                      :eq_gdp, :eq_optimal_kl, :eq_taylor, :eq_fisher,
                                      :eq_nominal_wage_inflation, :eq_fiscal_rule,
                                      :eq_g_budget_constraint, :LR, :LI, :LY,
                                      :LW, :LX, :eq_b, :eq_g, :eq_z, :eq_μ, :eq_λ_w,
                                      :eq_λ_f, :eq_rm]) #, :eq_consumption])

    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = [:c_t1]

    observables = keys(m.observable_mappings)
    endo   = m.endogenous_states
    eqcond = equilibrium_conditions

    for (i,k) in enumerate(endogenous_states);m.endogenous_states[k] = i end
    for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i+length(endogenous_states) end
    for (i,k) in enumerate(exogenous_shocks); m.exogenous_shocks[k]  = i end
    for (i,k) in enumerate(observables);      m.observables[k]       = i end
    for (i,k) in enumerate(equilibrium_conditions);      m.equilibrium_conditions[k]      = i end
end

function init_parameters!(m::RepDSGEGovDebt, het::HetDSGEGovDebt)
    for param in deepcopy(het.parameters)
        m <= param
    end
end

function init_grids!(m::RepDSGEGovDebt, het::HetDSGEGovDebt)
    m.grids = deepcopy(het.grids)
end

function steadystate!(m::RepDSGEGovDebt)
    ns = get_setting(m, :ns)

    #m <= SteadyStateParameter(:βstar, het[:βstar].value, description = "Steady-state discount factor", tex_label = "\\beta_*")
    #m[:r] = exp(m[:γ])/m[:βstar]
    m <= SteadyStateParameter(:Rkstar, m[:r] + m[:δ],
                              description = "Rental rate on capital", tex_label = "Rk_*")
    m <= SteadyStateParameter(:ωstar,
                              (m[:α]^(m[:α]/(1-m[:α])))*(1-m[:α])*m[:Rkstar]^(-m[:α]/(1-m[:α])),
                              description = "Real wage", tex_label = "\\omega_*")
    m <= SteadyStateParameter(:klstar, (m[:α]/(1-m[:α]))*(m[:ωstar]/m[:Rkstar])*exp(m[:γ]),
                              description = "Capital/Labor ratio", tex_label = "kl_*")
    m <= SteadyStateParameter(:kstar, m[:klstar]*m[:H], description = "Capital",
                              tex_label = "k_*")
    m <= SteadyStateParameter(:ystar, (exp(-m[:α]*m[:γ])*(m[:kstar]^m[:α])*m[:H]^(1-m[:α])),
                              description = "GDP", tex_label = "y_*")
    m <= SteadyStateParameter(:xstar, (1-(1-m[:δ])*exp(-m[:γ]))*m[:kstar],
                              description = "Investment", tex_label = "x_*")
    m <= SteadyStateParameter(:cstar, m[:ystar]/m[:g] - m[:xstar],
                              description = "Steady-state consumption", tex_label = "c_*")
end

function steadystate!(m::RepDSGEGovDebt, het::HetDSGEGovDebt;
                      same_beta::Bool = false)
    ns = get_setting(m, :ns)

    if same_beta
        m <= SteadyStateParameter(:βstar, het[:βstar].value, description = "Steady-state discount factor", tex_label = "\\beta_*")
        m[:r] = exp(m[:γ])/m[:βstar]
        m <= SteadyStateParameter(:Rkstar, m[:r] + m[:δ],
                                  description = "Rental rate on capital", tex_label = "Rk_*")
        m <= SteadyStateParameter(:ωstar,
                                  (m[:α]^(m[:α]/(1-m[:α])))*(1-m[:α])*m[:Rkstar]^(-m[:α]/(1-m[:α])),
                                  description = "Real wage", tex_label = "\\omega_*")
        m <= SteadyStateParameter(:klstar, (m[:α]/(1-m[:α]))*(m[:ωstar]/m[:Rkstar])*exp(m[:γ]),
                                  description = "Capital/Labor ratio", tex_label = "kl_*")
        m <= SteadyStateParameter(:kstar, m[:klstar]*m[:H], description = "Capital",
                                  tex_label = "k_*")
        m <= SteadyStateParameter(:ystar, (exp(-m[:α]*m[:γ])*(m[:kstar]^m[:α])*m[:H]^(1-m[:α])),
                                  description = "GDP", tex_label = "y_*")
        m <= SteadyStateParameter(:xstar, (1-(1-m[:δ])*exp(-m[:γ]))*m[:kstar],
                                  description = "Investment", tex_label = "x_*")
        m <= SteadyStateParameter(:cstar, m[:ystar]/m[:g] - m[:xstar],
                                  description = "Steady-state consumption", tex_label = "c_*")
    else
        m <= SteadyStateParameter(:βstar, exp(m[:γ])/(1 + m[:r]), description = "Steady-state discount factor",
                                  tex_label = "\\beta_*")
        m <= SteadyStateParameter(:Rkstar, m[:r] + m[:δ],
                                  description = "Rental rate on capital", tex_label = "Rk_*")
        m <= SteadyStateParameter(:ωstar,
                                  (m[:α]^(m[:α]/(1-m[:α])))*(1-m[:α])*m[:Rkstar]^(-m[:α]/(1-m[:α])),
                                  description = "Real wage", tex_label = "\\omega_*")
        m <= SteadyStateParameter(:klstar, (m[:α]/(1-m[:α]))*(m[:ωstar]/m[:Rkstar])*exp(m[:γ]),
                                  description = "Capital/Labor ratio", tex_label = "kl_*")
        m <= SteadyStateParameter(:kstar, m[:klstar]*m[:H], description = "Capital",
                                  tex_label = "k_*")
        m <= SteadyStateParameter(:ystar, (exp(-m[:α]*m[:γ])*(m[:kstar]^m[:α])*m[:H]^(1-m[:α])),
                                  description = "GDP", tex_label = "y_*")
        m <= SteadyStateParameter(:xstar, (1-(1-m[:δ])*exp(-m[:γ]))*m[:kstar],
                                  description = "Investment", tex_label = "x_*")
        m <= SteadyStateParameter(:cstar, (m.grids[:weights_total].*het[:μstar].value)'*(min.(1.0 ./ het[:lstar].value, repeat(m.grids[:xgrid].points, ns) .+ m[:η])),
                                  description = "Steady-state consumption", tex_label = "c_*")
    end
    m <= SteadyStateParameter(:lstar, 1/m[:cstar],
                              description = "Steady-state expected discounted
                              marginal utility of consumption", tex_label = "l_*")
    m <= SteadyStateParameter(:bg, m[:BoverY]*m[:ystar]*exp(m[:γ]), description = "Govt Debt",
                              tex_label = "bg")
    m <= SteadyStateParameter(:Tg, (exp(-m[:γ]) - 1/(1+m[:r]))*m[:bg] + (1-(1/m[:g]))*m[:ystar],
                              description = "Net lump sump taxes", tex_label = "Tg")
    m <= SteadyStateParameter(:Tstar, m[:Rkstar]*m[:kstar]*exp(-m[:γ]) - m[:xstar] - m[:Tg],
                              description = "Net transfer to households", tex_label = "T_*")

    return m
end


text:
using reproSim
using Test

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