text:
@testset "linear complementarity problem" begin
    M = [0  0 -1 -1 ;
         0  0  1 -2 ;
         1 -1  2 -2 ;
         1  2 -2  4 ]

    q = [2; 2; -2; -6]

    myfunc(x) = M*x + q

    n = 4
    lb = zeros(n)
    ub = 100*ones(n)

    status, z, f = solveLCP(myfunc, M, lb, ub)
    @show status
    @show z
    @show f
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
    @test status == :Solved

    status, z, f = solveLCP(myfunc, lb, ub)
    @show status
    @show z
    @show f
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
    @test status == :Solved

    println("-------------------------------------------------------")

    M = [0  0 -1 -1 ;
         0  0  1 -2 ;
         1 -1  2 -2 ;
         1  2 -2  4 ]

    q = [2; 2; -2; -6]

    myfunc(x) = M*x + q

    n = 4
    lb = zeros(n)
    ub = 100*ones(n)

    options(convergence_tolerance=1e-2, output=:no, time_limit=3600)

    status, z, f = solveLCP(myfunc, M, lb, ub)
    @show status
    @show z
    @show f
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
    @test status == :Solved

    status, z, f = solveLCP(myfunc, lb, ub)
    @show status
    @show z
    @show f
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
    @test status == :Solved

    println("-------------------------------------------------------")

    function elemfunc(x)
        val = similar(x)
        val[1] = -x[3]-x[4] + q[1]
        val[2] = x[3] -2x[4] + q[2]
        val[3] = x[1]-x[2]+2x[3]-2x[4] + q[3]
        val[4] = x[1]+2x[2]-2x[3]+4x[4] + q[4]
        return val
    end

    n = 4
    lb = zeros(n)
    ub = 100*ones(n)

    var_name = ["x1", "x2", "x3", "x4"]
    con_name = ["F1", "F2", "F3", "F4"]

    options(convergence_tolerance=1e-2, output=:yes, time_limit=3600)

    status, z, f = solveLCP(elemfunc, M, lb, ub)
    status, z, f = solveLCP(elemfunc, M, lb, ub, var_name)
    status, z, f = solveLCP(elemfunc, M, lb, ub, var_name, con_name)
    status, z, f = solveLCP(elemfunc, lb, ub)
    status, z, f = solveLCP(elemfunc, lb, ub, var_name)
    status, z, f = solveLCP(elemfunc, lb, ub, var_name, con_name)

    @show status
    @show z
    @show f

    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
    @test status == :Solved
end


text:
"""
    assert_allclose(actual, reference; [rtol], [atol], [err_msg], [error_location_fmt])

Raises an assertion exception when any entries in `actual` and `reference` are not equal up to the
desired tolerance. In other words, all entries of `actual` and `reference` must satisfy:
`abs(actual - reference) <= rtol * abs(reference) + atol`

This is inspired by numpy's assert_allclose function. However, we've changed some default values.
Additionally, in cases where actual and reference don't have equal values but do have the same
shape, the locations of the max errors are also printed.

# Keywords
- `rtol`: The desired relative tolerance
- `atol`: The desired absolute tolerance
- `err_msg`: An optional error message that can be printed when the arrays are not close
- `error_location_fmt`: An optional keyword used to specify a function that creates a string
  describing the location in the arrays where values are unequal, given a `CartesianIndex` object.
  This should not be specified in cases where actual and reference have different shapes.
"""
function assert_allclose(actual, reference; rtol = 1e-7, atol = 0.0, err_msg = nothing,
                         error_location_fmt = nothing)
    # make the initial check fast!
    if all( abs.(actual .- reference) .<= (rtol .* abs.(reference) .+ atol))
        return true
    end

    common_layout = size(actual) == size(reference)
    if isnothing(error_location_fmt)
        error_location_fmt = (arg) -> string( (arg isa Integer) ? arg : arg.I)
    elseif !common_layout
        throw("error_location_fmt must be nothing b/c actual and reference have different shapes")
    end

    # now let's provide a detailed error message (this can take longer):

    # determine what the max error is (and possibly where it happended)
    diff = abs.(actual .- reference)
    max_abs_diff_loc = argmax(diff)
    max_abs_diff = diff[max_abs_diff_loc]

    relative_diff = diff ./ abs.(reference)
    max_rel_error_magnitude_loc = argmax(relative_diff)
    max_rel_error_magnitude = relative_diff[max_rel_error_magnitude_loc]

    # format the message
    lines = isnothing(err_msg) ? [] : [err_msg]
    push!(lines, "The arrays aren't consistent to within atol = $atol, rtol = $rtol")

    if (rtol != 0.0) || (rtol == atol == 0.0)
        suffix = (common_layout) ? " at " * error_location_fmt(max_rel_error_magnitude_loc) : ""
        push!(lines, string("Max Rel Diff Magnitude: ",
                            isinf(max_rel_error_magnitude) ? "∞" : max_rel_error_magnitude,
                            suffix))
    end

    if (atol != 0.0) || (rtol == atol == 0.0)
        suffix = (common_layout) ? " at " * error_location_fmt(max_abs_diff_loc) : ""
        push!(lines, string("Max Abs Diff Magnitude: ", max_abs_diff, suffix))
    end

    throw(ErrorException(join(lines, "\n")))
    false
end

"""
    assert_allclose_grid(actual, reference, independent_vars; [rtol], [atol], [err_msg])

This function is a convenience wrapper around `assert_allclose` that provides slightly nicer error
messages based on the independent variables generating the grid.  It requires that `actual` and
`reference` have the same shape.

`independent_vars` should be a `ndims(actual)` element array. The ith element of `independent_vars`
should be 2-tuple or 3-tuple:
- `independent_vars[i][1]` should specify the name of the independent variable along the ith axis
- `independent_vars[i][2]` should specify the values of the independent variable along the ith axis
- `independent_vars[i][3]`, if present, should specify the units of the independent variable.
"""
function assert_allclose_grid(actual, reference, independent_vars; kwargs...)
    @assert length(independent_vars) == ndims(actual)

    function error_location_fmt(index)
        ax_strs = map(1:ndims(actual)) do ax
            if length(independent_vars[ax]) ∉ [2,3]
                error("independent_vars[$ax] should only have 2 or 3 elements")
            end

            name, vals = independent_vars[ax][1:2]
            unit_suffix = (length(independent_vars[ax]) == 3) ? " $(independent_vars[ax][3])" : ""

            if length(vals) != size(actual)[ax]
                error("length(independent_vars[$ax][2]) should be $(length(vals))")
            end

            "$name = $(vals[index[ax]])" * unit_suffix
        end
        join(ax_strs, ", ")
    end
    assert_allclose(actual, reference; error_location_fmt = error_location_fmt, kwargs...)
end

ext: jl


text:
export NCDstack, NCDarray

const NCD = NCDatasets

const UNNAMED_NCD_FILE_KEY = "unnamed"

const NCD_FILL_TYPES = (Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Float32,Float64,Char,String)

# CF standards don't enforce dimension names.
# But these are common, and should take care of most dims.
const NCD_DIMMAP = Dict(
    "lat" => Y,
    "latitude" => Y,
    "lon" => X,
    "long" => X,
    "longitude" => X,
    "time" => Ti,
    "lev" => Z,
    "mlev" => Z,
    "level" => Z,
    "vertical" => Z,
    "x" => X,
    "y" => Y,
    "z" => Z,
    "band" => Band,
)

haslayers(::Type{NCDfile}) = true
defaultcrs(::Type{NCDfile}) = EPSG(4326)
defaultmappedcrs(::Type{NCDfile}) = EPSG(4326)

# Raster ########################################################################

@deprecate NCDarray(args...; kw...) Raster(args...; source=NCDfile, kw...)

function Raster(ds::NCD.NCDataset, filename::AbstractString, key=nothing; kw...)
    key = _firstkey(ds, key)
    Raster(ds[key], filename, key; kw...)
end

_firstkey(ds::NCD.NCDataset, key::Nothing=nothing) = Symbol(first(layerkeys(ds)))
_firstkey(ds::NCD.NCDataset, key) = Symbol(key)

function FileArray(var::NCD.CFVariable, filename::AbstractString; kw...)
    da = RasterDiskArray{NCDfile}(var)
    size_ = size(da)
    eachchunk = DA.eachchunk(da)
    haschunks = DA.haschunks(da)
    T = eltype(var)
    N = length(size_)
    FileArray{NCDfile,T,N}(filename, size_; eachchunk, haschunks, kw...)
end

function Base.open(f::Function, A::FileArray{NCDfile}; write=A.write, kw...)
    _open(NCDfile, filename(A); key=key(A), write, kw...) do var
        f(RasterDiskArray{NCDfile}(var, DA.eachchunk(A), DA.haschunks(A)))
    end
end

"""
    Base.write(filename::AbstractString, ::Type{NCDfile}, s::AbstractRaster)

Write an NCDarray to a NetCDF file using NCDatasets.jl

Returns `filename`.
"""
function Base.write(filename::AbstractString, ::Type{NCDfile}, A::AbstractRaster)
    ds = NCD.Dataset(filename, "c"; attrib=_attribdict(metadata(A)))
    try
        _ncdwritevar!(ds, A)
    finally
        close(ds)
    end
    return filename
end

# Stack ########################################################################

@deprecate NCDstack(args...; kw...) RasterStack(args...; source=NCDfile, kw...)

"""
    Base.write(filename::AbstractString, ::Type{NCDfile}, s::AbstractRasterStack)

Write an NCDstack to a single netcdf file, using NCDatasets.jl.

Currently `Metadata` is not handled for dimensions, and `Metadata`
from other [`AbstractRaster`](@ref) @types is ignored.
"""
function Base.write(filename::AbstractString, ::Type{NCDfile}, s::AbstractRasterStack)
    ds = NCD.Dataset(filename, "c"; attrib=_attribdict(metadata(s)))
    try
        map(key -> _ncdwritevar!(ds, s[key]), keys(s))
    finally
        close(ds)
    end
    return filename
end

function create(filename, ::Type{NCDfile}, T::Union{Type,Tuple}, dims::DimTuple;
    name=:layer1, keys=(name,), layerdims=map(_->dims, keys), missingval=nothing, metadata=NoMetadata()
)
    types = T isa Tuple ? T : Ref(T)
    missingval = T isa Tuple ? missingval : Ref(missingval)
    # Create layers of zero arrays
    layers = map(layerdims, keys, types, missingval) do lds, key, t, mv
        A = FillArrays.Zeros{t}(map(length, lds))
        Raster(A, dims=lds; name=key, missingval=mv)
    end
    write(filename, NCDfile, Raster(first(layers)))
    return Raster(filename)
end

# DimensionalData methods for NCDatasets types ###############################

function DD.dims(ds::NCD.Dataset, crs=nothing, mappedcrs=nothing)
    map(_dimkeys(ds)) do key
        _ncddim(ds, key, crs, mappedcrs)
    end |> Tuple
end
function DD.dims(var::NCD.CFVariable, crs=nothing, mappedcrs=nothing)
    names = NCD.dimnames(var)
    map(names) do name
        _ncddim(var.var.ds, name, crs, mappedcrs)
    end |> Tuple
end

DD.metadata(ds::NCD.Dataset) = Metadata{NCDfile}(LA.metadatadict(ds.attrib))
DD.metadata(var::NCD.CFVariable) = Metadata{NCDfile}(LA.metadatadict(var.attrib))
DD.metadata(var::NCD.Variable) = Metadata{NCDfile}(LA.metadatadict(var.attrib))

function DD.layerdims(ds::NCD.Dataset)
    keys = Tuple(layerkeys(ds))
    dimtypes = map(keys) do key
        DD.layerdims(NCD.variable(ds, string(key)))
    end
    NamedTuple{map(Symbol, keys)}(dimtypes)
end
function DD.layerdims(var::NCD.Variable)
    map(NCD.dimnames(var)) do dimname
        _ncddimtype(dimname)()
    end
end

DD.layermetadata(ds::NCD.Dataset) = _layermetadata(ds, Tuple(layerkeys(ds)))
function _layermetadata(ds, keys)
    dimtypes = map(k -> DD.metadata(NCD.variable(ds, string(k))), keys)
    NamedTuple{map(Symbol, keys)}(dimtypes)
end

missingval(var::NCD.CFVariable) = missing
layermissingval(ds::NCD.Dataset) = missing

function layerkeys(ds::NCD.Dataset)
    dimkeys = _dimkeys(ds)
    toremove = if "bnds" in dimkeys
        dimkeys = setdiff(dimkeys, ("bnds",))
        boundskeys = String[]
        for k in dimkeys
            var = NCD.variable(ds, k)
            if haskey(var.attrib, "bounds")
                push!(boundskeys, var.attrib["bounds"])
            end
        end
        union(dimkeys, boundskeys)::Vector{String}
    else
        dimkeys::Vector{String}
    end
    return setdiff(keys(ds), toremove)
end

function FileStack{NCDfile}(ds::NCD.Dataset, filename::AbstractString; write=false, keys)
    keys = map(Symbol, keys isa Nothing ? layerkeys(ds) : keys) |> Tuple
    type_size_ec_hc = map(keys) do key
        var = ds[string(key)]
        Union{Missing,eltype(var)}, size(var), _ncd_eachchunk(var), _ncd_haschunks(var)
    end
    layertypes = NamedTuple{keys}(map(x->x[1], type_size_ec_hc))
    layersizes = NamedTuple{keys}(map(x->x[2], type_size_ec_hc))
    eachchunk = NamedTuple{keys}(map(x->x[3], type_size_ec_hc))
    haschunks = NamedTuple{keys}(map(x->x[4], type_size_ec_hc))
    return FileStack{NCDfile}(filename, layertypes, layersizes, eachchunk, haschunks, write)
end

# Utils ########################################################################

function _open(f, ::Type{NCDfile}, filename::AbstractString; key=nothing, write=false)
    lookup = write ? "a" : "r"
    if key isa Nothing
        NCD.Dataset(cleanreturn ∘ f, filename, lookup)
    else
        NCD.Dataset(filename, lookup) do ds
            cleanreturn(f(ds[_firstkey(ds, key)]))
        end
    end
end

cleanreturn(A::NCD.CFVariable) = Array(A)

function _ncddim(ds, dimname::Key, crs=nothing, mappedcrs=nothing)
    if haskey(ds, dimname)
        D = _ncddimtype(dimname)
        lookup = _ncdlookup(ds, dimname, D, crs, mappedcrs)
        return D(lookup)
    else
        # The var doesn't exist. Maybe its `complex` or some other marker,
        # so make it a custom `Dim` with `NoLookup`
        len = _ncfinddimlen(ds, dimname)
        len === nothing && _unuseddimerror()
        return Dim{Symbol(dimname)}(NoLookup(Base.OneTo(len)))
    end
end

function _ncfinddimlen(ds, dimname)
    for key in keys(ds)
        var = NCD.variable(ds, key)
        dimnames = NCD.dimnames(var)
        if dimname in dimnames
            return size(var)[findfirst(==(dimname), dimnames)]
        end
    end
    return nothing
end

# Find the matching dimension constructor. If its an unknown name
# use the generic Dim with the dim name as type parameter
_ncddimtype(dimname) = haskey(NCD_DIMMAP, dimname) ? NCD_DIMMAP[dimname] : DD.basetypeof(DD.key2dim(Symbol(dimname)))

# _ncdlookup
# Generate a `LookupArray` from a netcdf dim.
function _ncdlookup(ds::NCD.Dataset, dimname, D, crs, mappedcrs)
    dvar = ds[dimname]
    index = dvar[:]
    metadata = Metadata{NCDfile}(LA.metadatadict(dvar.attrib))
    return _ncdlookup(ds, dimname, D, index, metadata, crs, mappedcrs)
end
# For unknown types we just make a Categorical lookup
function _ncdlookup(ds::NCD.Dataset, dimname, D, index::AbstractArray, metadata, crs, mappedcrs)
    Categorical(index; metadata=metadata)
end
# For Number and AbstractTime we generate order/span/sampling
function _ncdlookup(
    ds::NCD.Dataset, dimname, D, index::AbstractArray{<:Union{Number,Dates.AbstractTime}},
    metadata, crs, mappedcrs
)
    # Assume the locus is at the center of the cell if boundaries aren't provided.
    # http://cfconventions.org/cf-conventions/cf-conventions.html#cell-boundaries
    order = _ncdorder(index)
    var = NCD.variable(ds, dimname)
    if haskey(var.attrib, "bounds")
        boundskey = var.attrib["bounds"]
        boundsmatrix = Array(ds[boundskey])
        span, sampling = Explicit(boundsmatrix), Intervals(Center())
        return _ncdlookup(D, index, order, span, sampling, metadata, crs, mappedcrs)
    elseif eltype(index) <: Dates.AbstractTime
        span, sampling = _ncdperiod(index, metadata)
        return _ncdlookup(D, index, order, span, sampling, metadata, crs, mappedcrs)
    else
        span, sampling = _ncdspan(index, order), Points()
        return _ncdlookup(D, index, order, span, sampling, metadata, crs, mappedcrs)
    end
end
# For X and Y use a Mapped <: AbstractSampled lookup
function _ncdlookup(
    D::Type{<:Union{<:XDim,<:YDim}}, index, order::Order, span, sampling, metadata, crs, mappedcrs
)
    # If the index is regularly spaced and there is no crs
    # then there is probably just one crs - the mappedcrs
    crs = if crs isa Nothing && span isa Regular
        mappedcrs
    else
        crs
    end
    dim = DD.basetypeof(D)()
    return Mapped(index, order, span, sampling, metadata, crs, mappedcrs, dim)
end
# Band dims have a Categorical lookup, with order
function _ncdlookup(D::Type{<:Band}, index, order::Order, span, sampling, metadata, crs, mappedcrs)
    Categorical(index, order, metadata)
end
# Otherwise use a regular Sampled lookup
function _ncdlookup(D::Type, index, order::Order, span, sampling, metadata, crs, mappedcrs)
    Sampled(index, order, span, sampling, metadata)
end

function _ncdorder(index)
    index[end] > index[1] ? ForwardOrdered() : ReverseOrdered()
end

function _ncdspan(index, order)
    # Handle a length 1 index
    length(index) == 1 && return Regular(zero(eltype(index)))
    step = index[2] - index[1]
    for i in 2:length(index)-1
        # If any step sizes don't match, its Irregular
        if !(index[i+1] - index[i] ≈ step)
            bounds = if length(index) > 1
                beginhalfcell = abs((index[2] - index[1]) * 0.5)
                endhalfcell = abs((index[end] - index[end-1]) * 0.5)
                if LA.isrev(order)
                    index[end] - endhalfcell, index[1] + beginhalfcell
                else
                    index[1] - beginhalfcell, index[end] + endhalfcell
                end
            else
                index[1], index[1]
            end
            return Irregular(bounds)
        end
    end
    # Otherwise regular
    return Regular(step)
end

# delta_t and ave_period are not CF standards, but CDC
function _ncdperiod(index, metadata::Metadata{NCDfile})
    if haskey(metadata, :delta_t)
        period = _parse_period(metadata[:delta_t])
        period isa Nothing || return Regular(period), Points()
    elseif haskey(metadata, :avg_period)
        period = _parse_period(metadata[:avg_period])
        period isa Nothing || return Regular(period), Intervals(Center())
    end
    return sampling = Irregular((nothing, nothing)), Points()
end

function _parse_period(period_str::String)
    regex = r"(\d\d\d\d)-(\d\d)-(\d\d) (\d\d):(\d\d):(\d\d)"
    mtch = match(regex, period_str)
    if mtch === nothing
        return nothing
    else
        vals = Tuple(parse.(Int, mtch.captures))
        periods = (Year, Month, Day, Hour, Minute, Second)
        if length(vals) == length(periods)
            compound = sum(map((p, v) -> p(v), periods, vals))
            if length(compound.periods) == 1
                return compound.periods[1]
            else
                return compound
            end
        else
            return nothing
        end
    end
end

_attribdict(md::Metadata{NCDfile}) = Dict{String,Any}(string(k) => v for (k, v) in md)
_attribdict(md) = Dict{String,Any}()

_dimkeys(ds::NCD.Dataset) = keys(ds.dim)

# Add a var array to a dataset before writing it.
function _ncdwritevar!(ds::NCD.Dataset, A::AbstractRaster{T,N}) where {T,N}
    _def_dim_var!(ds, A)
    attrib = _attribdict(metadata(A))
    # Set _FillValue
    eltyp = _notmissingtype(Base.uniontypes(T)...)
    if ismissing(missingval(A))
        fillval = if haskey(attrib, "_FillValue") && attrib["_FillValue"] isa eltyp
            attrib["_FillValue"]
        else
            NCD.fillvalue(eltyp)
        end
        attrib["_FillValue"] = fillval
        A = replace_missing(A, fillval)
    elseif missingval(A) isa T
        attrib["_FillValue"] = missingval(A)
    else
        missingval(A) isa Nothing || @warn "`missingval` $(missingval(A)) is not the same type as your data $T."
    end

    key = if string(name(A)) == ""
        UNNAMED_NCD_FILE_KEY
    else
        string(name(A))
    end

    dimnames = lowercase.(string.(map(name, dims(A))))
    var = NCD.defVar(ds, key, eltyp, dimnames; attrib=attrib)
    # TODO do this with DiskArrays broadcast ??
    var[:] = parent(read(A))
end

_def_dim_var!(ds::NCD.Dataset, A) = map(d -> _def_dim_var!(ds, d), dims(A))
function _def_dim_var!(ds::NCD.Dataset, dim::Dimension)
    dimkey = lowercase(string(name(dim)))
    haskey(ds.dim, dimkey) && return nothing
    NCD.defDim(ds, dimkey, length(dim))
    lookup(dim) isa NoLookup && return nothing

    # Shift index before conversion to Mapped
    dim = _ncdshiftlocus(dim)
    if dim isa Y || dim isa X
        dim = convertlookup(Mapped, dim)
    end
    attrib = _attribdict(metadata(dim))
    if span(dim) isa Explicit
        bounds = val(span(dim))
        boundskey = get(metadata(dim), :bounds, string(dimkey, "_bnds"))
        push!(attrib, "bounds" => boundskey)
        NCD.defVar(ds, boundskey, bounds, ("bnds", dimkey))
    end
    NCD.defVar(ds, dimkey, Vector(index(dim)), (dimkey,); attrib=attrib)
    return nothing
end

_notmissingtype(::Type{Missing}, next...) = _notmissingtype(next...)
_notmissingtype(x::Type, next...) = x in NCD_FILL_TYPES ? x : _notmissingtype(next...)
_notmissingtype() = error("Your data is not a type that netcdf can store")

_ncdshiftlocus(dim::Dimension) = _ncdshiftlocus(lookup(dim), dim)
_ncdshiftlocus(::LookupArray, dim::Dimension) = dim
function _ncdshiftlocus(lookup::AbstractSampled, dim::Dimension)
    if span(lookup) isa Regular && sampling(lookup) isa Intervals
        # We cant easily shift a DateTime value
        if eltype(dim) isa Dates.AbstractDateTime
            if !(locus(dim) isa Center)
                @warn "To save to netcdf, DateTime values should be the interval Center, rather than the $(nameof(typeof(locus(dim))))"
            end
            dim
        else
            shiftlocus(Center(), dim)
        end
    else
        dim
    end
end

_unuseddimerror(dimname) = error("Netcdf contains unused dimension $dimname")

function _ncd_eachchunk(var)
    # chunklookup, chunkvec = NCDatasets.chunking(var)
    # chunksize = chunklookup == :chunked ? Tuple(chunkvec) :
    chunksize = size(var)
    DA.GridChunks(var, chunksize)
end

function _ncd_haschunks(var)
    # chunklookup, _ = NCDatasets.chunking(var)
    # chunklookup == :chunked ? DA.Chunked() :
    DA.Unchunked()
end

# precompilation

const _NCDVar = NCDatasets.CFVariable{Union{Missing, Float32}, 3, NCDatasets.Variable{Float32, 3, NCDatasets.NCDataset}, NCDatasets.Attributes{NCDatasets.NCDataset{Nothing}}, NamedTuple{(:fillvalue, :scale_factor, :add_offset, :calendar, :time_origin, :time_factor), Tuple{Float32, Nothing, Nothing, Nothing, Nothing, Nothing}}}
precompile(Rasters.FileArray, (_NCDVar, String))
precompile(layerkeys, (NCDatasets.NCDataset{Nothing},))
precompile(dims, (_NCDVar,Symbol))
precompile(dims, (_NCDVar,Symbol,Nothing,Nothing))
precompile(dims, (_NCDVar,Symbol,Nothing,EPSG))
precompile(dims, (_NCDVar,Symbol,EPSG,EPSG))
precompile(_firstkey, (NCDatasets.NCDataset{Nothing},))
precompile(_ncddim, (NCDatasets.NCDataset{Nothing}, Symbol, Nothing, Nothing))
precompile(_ncddim, (NCDatasets.NCDataset{Nothing}, Symbol, Nothing, EPSG))
precompile(_ncddim, (NCDatasets.NCDataset{Nothing}, Symbol, EPSG, EPSG))
precompile(Raster, (NCDatasets.NCDataset{Nothing}, String, Nothing))
precompile(Raster, (NCDatasets.NCDataset{Nothing}, String, Symbol))
precompile(Raster, (_NCDVar, String, Symbol))

precompile(geoarray, (String,))
precompile(Raster, (String,))


text:
# LPS25H pressure/temperature sensor
# See http://www.st.com/resource/en/datasheet/lps25h.pdf

"""
    LPS25H_pressure()

The atmospheric pressure (in millibars) from the LPS25H sensor.
"""
function LPS25H_pressure()
    setaddr(LPS25H_ADDRESS)
    CTRL_REG1 = 0x20
    CTRL_REG2 = 0x21

    smbus_write(CTRL_REG1, 0x00) # power down
    smbus_write(CTRL_REG1, 0x84) # power on, block update
    smbus_write(CTRL_REG2, 0x01) # one-shot aquisition

    while true
        sleep(0.025)
        smbus_read(CTRL_REG2) == 0 && break # check if finished
    end

    pressure = (UInt64.(smbus_read(0x2a)) << 16 |
                UInt64.(smbus_read(0x29)) << 8 |
                UInt64.(smbus_read(0x28)))/4096.0

    smbus_write(CTRL_REG1, 0x00) # power down
    return pressure
end

"""
    LPS25H_temperature()

The temperature (in °C) from the LPS25H sensor.
"""
function LPS25H_temperature()
    setaddr(LPS25H_ADDRESS)
    CTRL_REG1 = 0x20
    CTRL_REG2 = 0x21

    smbus_write(CTRL_REG1, 0x00) # power down
    smbus_write(CTRL_REG1, 0x84) # power on, block update
    smbus_write(CTRL_REG2, 0x01) # one-shot aquisition

    while true
        sleep(0.025)
        smbus_read(CTRL_REG2) == 0 && break # check if finished
    end

    temp = 42.5 + (Int16(smbus_read(0x2c)) << 8 | smbus_read(0x2b)) / 480

    smbus_write(CTRL_REG1, 0x00) # power down
    return temp
end


text:
"""
测试上VHS的kagome
"""

using Plots
using MARY_fRG.Drawers
using MARY_fRG.Fermi.Kagome
using MARY_fRG.Fermi.Surface
using MARY_fRG.FlowEquation

#阻止画图
ENV["GKSwstype"] = "100"

"""
获取极低温的时候的bubble
"""
function get_ult_bubb(Γ4; maxnum=1000000)
    Γ4_ult = fliter_away_surface(Γ4, 2.0)
    Γ4_ult = refine_to_surface(Γ4_ult)
    while length(Γ4_ult.ltris) < maxnum
        Γ4_ult = refine_to_surface(Γ4_ult)
    end
    plt = draw_points([tri.center for tri in Γ4_ult.ltris])
    @savepng plt "ult_surface"
    maxi, mini = engpeak_to_surface(Γ4_ult)
    bubb_pp_u, bubb_fs_u, bubb_ex_u = all_bubble_tf_ult_mt(Γ4_ult, maxi)
    return bubb_pp_u, bubb_fs_u, bubb_ex_u
end

function run_tf()
    model = upperband_kagome_lattice(-0.3)
    Γ4 = TFGamma4(
        model, 6.0, 24, 233
    )
    Γ4.V[1, 1, 1, 1, :, :, :] .+= get_kagome_U_mt(3.0, Γ4.patches[1])
    println(
        all(isfinite.(Γ4.V[1, 1, 1, 1, :, :, :]))
    )
    #plt = heatmap(Γ4.V[1, 1, 1, 1, :, :, 1])
    #png(plt, "Gamma41")
    #return
    lval = 0.
    lstep = 0.01
    bpu, bfu, beu = get_ult_bubb(Γ4)
    @info "ult结束"
    for idx in 1:1:5001
        bubb_pp, bubb_fs, bubb_ex = all_bubble_tf_mt(Γ4, lval)
        dl = dl_tf_mix_ult_mt(Γ4, bubb_pp, bubb_fs, bubb_ex, bpu, bfu, beu)
        Γ4.V .+= dl .* lstep
        if idx % 50 == 1 || idx > 2001
            plt = heatmap(Γ4.V[1, 1, 1, 1, :, :, 1])
            png(plt, "Gamma4"*string(idx))
        end
        lval += lstep
    end
end

run_tf()

"""
绘制patch和费米面
"""
function draw_brlu()
    model = upperband_kagome_lattice(-0.2)
    Γ4 = ECGamma4(
        model, 6.0, 24, 223
    )
    plt = draw_lines(Γ4.model.brillouin.edges)
    suf = const_energy_line(Γ4.ltris, Γ4.ladjs, 0., Γ4.model.dispersion[1])
    draw_lines!(plt, suf; lcidx=ones(Int64, length(suf)))
    draw_points!(plt, Γ4.patches[1], text=1:1:24)
    @savepng plt "kag_pat_surf"
end

#draw_brlu()


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    struct DensitySample

A weighted sample drawn according to an statistical density,
e.g. a [`BAT.AbstractDensity`](@ref).

Fields:
    * `params`: Multivariate parameter vector
    * `logdensity`: log of the value of the density at `params`
    * `weight`: Weight of the sample
    * `info`: Additional info on the provenance of the sample. Content depends
       on the sampling algorithm.
    * aux: Custom user-defined information attatched to the sample.

Constructors:

```julia
DensitySample(
    params::AbstractVector{<:Real},
    logdensity::Real,
    weight::Real,
    info::Any,
    aux::Any
)
```

Use [`DensitySampleVector`](@ref) to store vectors of multiple samples with
an efficient column-based memory layout.
"""
struct DensitySample{
    P<:Real,
    T<:Real,
    W<:Real,
    R,
    Q,
    PA<:AbstractVector{P}
}
    params::PA
    logdensity::T
    weight::W
    info::R
    aux::Q
end

export DensitySample

# DensitySample behaves as a scalar type under broadcasting:
@inline Base.Broadcast.broadcastable(shape::DensitySample) = Ref(shape)

import Base.==
function ==(A::DensitySample, B::DensitySample)
    A.params == B.params && A.logdensity == B.logdensity &&
        A.weight == B.weight && A.info == B.info && A.aux == B.aux
end

function Base.similar(s::DensitySample{P,T,W,R,Q}) where {P,T,W,R,Q}
    params = fill!(similar(s.params), oob(eltype(s.params)))
    logdensity = convert(T, NaN)
    weight = zero(W)
    info = R()
    aux = Q()
    PA = typeof(params)
    DensitySample{P,T,W,R,Q,PA}(params, logdensity, weight, info, aux)
end

nparams(s::DensitySample) = length(s.params)

function _apply_shape(shape::AbstractValueShape, s::DensitySample)
    (
        params = shape(s.params),
        logdensity = s.logdensity,
        weight = s.weight,
        info = s.info,
        aux = s.aux,
    )
end

@static if VERSION >= v"1.3"
    (shape::AbstractValueShape)(s::DensitySample) = _apply_shape(shape, s)
else
    (shape::ScalarShape)(s::DensitySample) = _apply_shape(shape, s)
    (shape::ArrayShape)(s::DensitySample) = _apply_shape(shape, s)
    (shape::ConstValueShape)(s::DensitySample) = _apply_shape(shape, s)
    (shape::NamedTupleShape)(s::DensitySample) = _apply_shape(shape, s)
end

"""
    DensitySampleVector

Type alias for `StructArrays.StructArray{<:DensitySample,...}`.

Constructor:

```julia
    DensitySampleVector(
        (
            params::AbstractVector{<:AbstractVector{<:Real}}
            logdensity::AbstractVector{<:Real}
            weight::AbstractVector{<:Real}
            info::AbstractVector{<:Any}
            aux::AbstractVector{<:Any}
        )
    )
```
"""
const DensitySampleVector{
    P<:Real,T<:AbstractFloat,W<:Real,R,Q,PA<:AbstractVector{P},
    PAV<:AbstractVector{<:AbstractVector{P}},TV<:AbstractVector{T},WV<:AbstractVector{W},RV<:AbstractVector{R},QV<:AbstractVector{Q}
} = StructArray{
    DensitySample{P,T,W,R,Q,PA},
    1,
    NamedTuple{(:params, :logdensity, :weight, :info, :aux), Tuple{PAV,TV,WV,RV,QV}}
}

export DensitySampleVector

function StructArray{DensitySample}(
    contents::Tuple{
        AbstractVector{<:AbstractVector{P}},
        AbstractVector{T},
        AbstractVector{W},
        AbstractVector{R},
        AbstractVector{Q},
    }
) where {P<:Real,T<:AbstractFloat,W<:Real,R,Q}
    params, logdensity, weight, info, aux = contents
    PA = eltype(params)
    StructArray{DensitySample{P,T,W,R,Q,PA}}(contents)
end

DensitySampleVector(contents::NTuple{5,Any}) = StructArray{DensitySample}(contents)

_create_undef_vector(::Type{T}, len::Integer) where T = Vector{T}(undef, len)

function DensitySampleVector{P,T,W,R,Q}(::UndefInitializer, len::Integer, npar::Integer) where {P<:Real,T<:AbstractFloat,W<:Real,R,Q}
    contents = (
        VectorOfSimilarVectors(ElasticArray{P}(undef, npar, len)),
        Vector{T}(undef, len),
        Vector{W}(undef, len),
        _create_undef_vector(R, len),
        _create_undef_vector(Q, len)
    )

    DensitySampleVector(contents)
end

DensitySampleVector(::Type{S}, nparams::Integer) where {P<:Real,T<:AbstractFloat,W<:Real,R,Q,S<:DensitySample{P,T,W,R,Q}} =
    DensitySampleVector{P,T,W,R,Q}(undef, 0, nparams)

# Specialize getindex to properly support ArraysOfArrays, preventing
# conversion to exact element type:
@inline Base.getindex(A::StructArray{<:DensitySample}, I::Int...) =
    DensitySample(A.params[I...], A.logdensity[I...], A.weight[I...], A.info[I...], A.aux[I...])

# Specialize IndexStyle, current default for StructArray seems to be IndexCartesian()
Base.IndexStyle(::StructArray{<:DensitySample, 1}) = IndexLinear()

# Specialize comparison, currently StructArray seems fall back to `(==)(A::AbstractArray, B::AbstractArray)`
import Base.==
function(==)(A::DensitySampleVector, B::DensitySampleVector)
    A.params == B.params &&
    A.logdensity == B.logdensity &&
    A.weight == B.weight &&
    A.info == B.info &&
    A.aux == B.aux
end

function Base.merge!(X::DensitySampleVector, Xs::DensitySampleVector...)
    for Y in Xs
        append!(X, Y)
    end
    X
end

Base.merge(X::DensitySampleVector, Xs::DensitySampleVector...) = merge!(deepcopy(X), Xs...)

function UnsafeArrays.uview(A::DensitySampleVector)
    DensitySampleVector((
        uview(A.params),
        uview(A.logdensity),
        uview(A.weight),
        uview(A.info),
        uview(A.aux)
    ))
end

Base.@propagate_inbounds function _bcasted_apply(shape::AbstractValueShape, A::DensitySampleVector)
    TypedTables.Table(
        params = shape.(A.params),
        logdensity = A.logdensity,
        weight = A.weight,
        info = A.info,
        aux = A.aux
    )
end

Base.copy(
    instance::Base.Broadcast.Broadcasted{
        <:Base.Broadcast.AbstractArrayStyle{1},
        <:Any,
        <:AbstractValueShape,
        <:Tuple{DensitySampleVector}
    }
) = _bcasted_apply(instance.f, instance.args[1])

"""
    bat_stats(samples::DensitySampleVector)

Calculated parameter statistics on `samples`. Returns a
`NamedTuple{(:mode,:mean,:cov,...)}`. Result properties not listed here
are not part of the stable BAT API and subject to change.
"""
function bat_stats(samples::DensitySampleVector)
    par_mean = mean(samples.params, FrequencyWeights(samples.weight))
    par_mode_idx = findmax(samples.logdensity)[2]
    par_mode = samples.params[par_mode_idx]
    par_cov = cov(samples.params, FrequencyWeights(samples.weight))

    (
        mode = par_mode,
        mean = par_mean,
        cov = par_cov
    )
end

export bat_stats

"""
    drop_low_weight_samples(
        samples::DensitySampleVector,
        fraction::Real = 10^-4
    )

*BAT-internal, not part of stable public API.*

Drop `fraction` of the total probability mass from samples to filter out the
samples with the lowest weight.
"""
function drop_low_weight_samples(samples::DensitySampleVector, fraction::Real = 10^-5)
    W = float(samples.weight)
    if minimum(W) / maximum(W) > 10^-2
        samples
    else
        W_s = sort(W)
        Q = cumsum(W_s)
        Q ./= maximum(Q)
        @assert last(Q) ≈ 1
        thresh = W_s[searchsortedlast(Q, fraction)]
        idxs = findall(x -> x >= thresh, samples.weight)
        samples[idxs]
    end
end


text:
include("optimizers.jl")

function update!(opt, x, g)
    x .-= apply!(opt, x, g)
end

function update!(opt, gs::Zygote.Grads, ps::Params)
    for p in ps
        isnothing(gs[p]) && continue
        update!(opt, p, gs[p])
    end
end

make_d_tuple(d::Tuple) = d
make_d_tuple(d) = tuple(d)

function train!(N, loss, data, opt)
    ps = params(N)
    for d in data
        try
            g = gradient(ps) do
                loss(make_d_tuple(d)...)
            end
            update!(opt, g, ps)
        catch
            @warn "you can't train by one data because something wrong with the data!"
        end
    end
end

"""
    @epochs n ex
This macro cruns `ex` `n` times. Basically this is useful for learning NeuralNetwork.

# Example
julia> a = 1

julia> @epochs 1000 a+=1
progress:1000/1000
julia>a
1001
"""
macro epochs(n, ex)
    :(for i in 1 : $(esc(n))
        progress = "progress:"*string(i)*"/"*string($(esc(n)))*"\r"
        print(progress)
        $(esc(ex))
    end)
end


text:
include("modulos/testing_models.jl")
include("modulos/datasets.jl")
include("modulos/attributes_from_dataset.jl")

#=
hay diferentes extracciones:
    · :A1 (aproximación 1)
    · :A21 (aproximación 2 extracción 1)
    · :A22 (aproximación 2 extracción 2)
    etc...
=#
extraction = :A33;
dataset_name="datasets/faces.data";
change = false;
if extraction==:A1
    x = 6;
    y = 7;
elseif extraction==:A21
    x = 42;
    y = 43;
elseif (extraction==:A22 || extraction==:A23)
    x = 36;
    y = 37;
elseif extraction==:A31
    x = 48;
    y = 49;
elseif extraction==:A32
    x = 66;
    y = 67;
elseif extraction==:A33
    x = 72;
    y = 73;
elseif extraction==:A34
    x = 84;
    y = 85;
end;
if (!isfile(dataset_name) || change)
    createDataset(dataset_name,extraction);
end;

dataset = readdlm(dataset_name,',');
inputs = convert(Array{Float64,2}, dataset[:,1:x]);
targets = convert(Array{Any,1},dataset[:,y]);
seed!(1);
numFolds = 10;

#=
Posibles representaciones
    · ninguna o :All = saca todo por pantalla.
    · :AccStd = accuracy y su std.
    · :F1Std = F1 y su std.
    · :AccF1 = accuracy y F1.
    · :Plot3D = solo para svm con kernel polinómico.

Para sacar los gráficos 3D:
sudo apt-get install python3-matplotlib
Pkg.add("PyPlot") || Pkg.build("PyPlot")
=#

# Entrenamos knn
testingModels(:KNN, Dict("maxNeighbors" => 20), inputs, targets, numFolds; rep=:All);

# Entrenamos los arboles de decision
testingModels(:DecisionTree, Dict("maxDepth" => 20), inputs, targets, numFolds; rep=:All);

# Entrenamos svm
modelHyperparameters = Dict();
modelHyperparameters["kernel"] = "rbf";
modelHyperparameters["kernelDegree"] = 15;
modelHyperparameters["maxGamma"] = 10;
testingModels(:SVM, modelHyperparameters, inputs, targets, numFolds; rep=:All);

# Entrenamos las RR.NN.AA.
modelHyperparameters = Dict();
modelHyperparameters["fstNeuron"] = 1;
modelHyperparameters["maxNNxlayer"] = 16;
modelHyperparameters["layers"] = 1;
testingModels(:ANN, modelHyperparameters, inputs, targets, numFolds; rep=:All);

modelHyperparameters = Dict();
modelHyperparameters["fstNeuron"] = 6;
modelHyperparameters["maxNNxlayer"] = 12;
modelHyperparameters["layers"] = 2;
testingModels(:ANN, modelHyperparameters, inputs, targets, numFolds; rep=:All);


text:
# Analytic magnetic field topology in 2D.
#
# This demo shows the two types of classical magnetic field null point in 2D
# calculated from an analytic function.
# In 2D, the magnetic field potential can be written as
# Az(x,y) = B₀/2*(x^2/dx^2 ± y^2/dy^2)
# where
# -: X-point
# +: O-point
# and the corresponding field can be expressed as
# B(x,y,z) = ẑ × ∇Az + ẑ Bg
# where Bg is the guide field.
#
# Two different methods of computing B are provided:
# 1. hand-calculated form
# 2. auto-differentiation

using LinearAlgebra
using FieldTracer
using ForwardDiff
using StaticArrays
using PyPlot

"Hand-calculated o-type magnetic field from 2D magnetic potential."
function o_type(X, Y, B₀, Bg, dx, dy)
   Bx = [-B₀/dy^2 * y for y in Y, x in X]
   By = [B₀/dx^2 * x for y in Y, x in X]
   Bz = fill(Bg, size(Bx))

   Bx, By, Bz
end

"Hand-calculated x-type magnetic field from 2D magnetic potential."
function x_type(X, Y, B₀, Bg, dx, dy)
   Bx = [B₀/dy^2 * y for y in Y, x in X]
   By = [B₀/dx^2 * x for y in Y, x in X]
   Bz = fill(Bg, size(Bx))

   Bx, By, Bz
end

"Return o-type magnetic field from 2D potential using auto-differentiation."
function o_type_autodiff(X, Y, B₀, Bg, dx, dy)

   Az(x, y) = B₀/2*(x^2/dx^2 + y^2/dy^2)
   ∇Az(x, y, z) = ForwardDiff.gradient(p -> Az(p[1], p[2]), SVector(x,y,z))
   ẑ = SVector(0.0, 0.0, 1.0)
   B = Array{Float64,3}(undef, 3, length(X), length(Y))
   for j in axes(Y,1), i in axes(X,1)
      B[:,i,j] = ẑ × ∇Az(X[i], Y[j], 0.0)
      B[3,i,j] += Bg
   end
   B
end

"Return x-type magnetic field from 2D potential using auto-differentiation."
function x_type_autodiff(X, Y, B₀, Bg, dx, dy)

   Az(x, y) = B₀/2*(x^2/dx^2 - y^2/dy^2)
   ∇Az(x, y, z) = ForwardDiff.gradient(p -> Az(p[1], p[2]), SVector(x,y,z))
   ẑ = SVector(0.0, 0.0, 1.0)
   B = Array{Float64,3}(undef, 3, length(X), length(Y))
   for j in axes(Y,1), i in axes(X,1)
      B[:,i,j] = ẑ × ∇Az(X[i], Y[j], 0.0)
      B[3,i,j] += Bg
   end
   B
end

## Main

# Parameters
Bg = 1.0 # Guide field, [nT]
B₀ = 1.0

fig, axs = plt.subplots(1, 2, figsize=(12,5))

Δ = 0.1           # grid interval
X = -5:Δ:5        # grid x range
Y = -5:Δ:5        # grid y range
dx, dy = 1.0, 1.0 # scale ratio of the analytic field

@time Bx, By, Bz = o_type(X, Y, B₀, Bg, dx, dy)
axs[1].streamplot(collect(X), collect(Y), Bx, By, minlength=0.1)
#@time B = o_type_autodiff(X, Y, B₀, Bg, dx, dy)
#axs[1].streamplot(collect(X), collect(Y), B[1,:,:]', B[2,:,:]', minlength=0.2)

xstart = 1.0:1.0:3.0 # seed starting point x
ystart = 1.0:1.0:3.0 # seed starting point y

for i in axes(xstart,1)
   x1, y1 = trace(Bx, By, xstart[i], ystart[i], X, Y;
      maxstep=3000, ds=0.1, gridType="meshgrid")
   line = axs[1].plot(x1, y1, "r--")
   add_arrow(line[1])
end
axs[1].plot(0, 0, color="k", marker="o", label="o-type null point")
axs[1].legend()
axs[1].axis("equal")

X = -5:Δ:5        # grid x range
Y = -3:Δ:3        # grid x range
dx, dy = 1.0, 0.1 # scale ratio of the analytic field

#Bx, By, Bz = x_type_autodiff(X, Y, B₀, Bg, dx, dy)
@time B = x_type_autodiff(X, Y, B₀, Bg, dx, dy)
axs[2].streamplot(collect(X), collect(Y), B[1,:,:]', B[2,:,:]')

xstart = [1.0, 3.0, 1.0] # seed starting point x
ystart = [0.0, 0.0, 1.0] # seed starting point y

for i in axes(xstart,1)
   x1, y1 = trace(B[1,:,:], B[2,:,:], xstart[i], ystart[i], X, Y;
      maxstep=3000, ds=0.1, gridType="ndgrid")
   axs[2].plot(x1, y1, "r--")
end
axs[2].plot(0, 0, color="k", marker="o", label="x-type null point")
axs[2].legend()
axs[2].axis("equal")

savefig("x_o_point.png", bbox_inches="tight")


text:
# Regex to validate a Roman numeral
const VALID_ROMAN_PATTERN =
    r"""
    ^\s*                   # Skip leading whitespace
    (
    M*                     # Thousands
    (C{0,9}|CD|DC{0,4}|CM) # Hundreds
    (X{0,9}|XL|LX{0,4}|XC) # Tens
    (I{0,9}|IV|VI{0,4}|IX) # Ones
    )
    \s*$                   # Skip trailing whitespace
    """ix # Be case-insensitive and verbose

const NUMERAL_MAP = [
    (1000, "M")
    (900,  "CM")
    (500,  "D")
    (400,  "CD")
    (100,  "C")
    (90,   "XC")
    (50,   "L")
    (40,   "XL")
    (10,   "X")
    (9,    "IX")
    (5,    "V")
    (4,    "IV")
    (1,    "I")
]

function parseroman(str::ASCIIString)
    m = match(VALID_ROMAN_PATTERN, str)
    m == nothing && throw(InvalidRomanError(str))
    # Strip whitespace
    str = m.captures[1]
    # Make `str` uppercase
    if !isupper(str); str = uppercase(str); end
    i = 1
    val = 0
    strlen = length(str)
    for (num_val, numeral) in NUMERAL_MAP
        numlen = length(numeral)
        while i+numlen-1 <= strlen && str[i:i+numlen-1] == numeral
            val += num_val
            i += numlen
        end
    end
    val
end

function toroman(val::Integer)
    val <= 0 && throw(DomainError())
    val > 5000 && warn("Roman numerals do not handle large numbers well")

    str = ""
    for (num_val, numeral) in NUMERAL_MAP
        i = div(val, num_val)
        # Never concatenate an empty string to `str`
        i == 0 && continue
        str *= (numeral^i)
        val -= i*num_val
        # Stop when ready
        val == 0 && break
    end
    str
end


text:
using Spark
sc = SparkContext(master="local")
path = string("file:///", ENV["SPARK_HOME"], "/README.md")
txt = text_file(sc, path)
# Normally we would use a flatmap, but currently only has map_partitions


text:
module MyType

using Lib

using BigLib: thingx, thingy

import Base.show

export MyType, fuck

struct MyType
    x
end

bar(n) = 2n

fuck(mt::MyType) = bar(mt.x) + 9527

show(io::IO, mt::MyType) = print(io, "1314")

end

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
pattern = r"(\d+)-(\d+) (\w): (\w+)"
rows = open("password_list.txt") do f
    match.(pattern, readlines(f))
end

println(count(row -> (parse(Int, row[1]) <= count(row[3], row[4]) <= parse(Int, row[2])), rows))

println(count(row -> ((row[4][parse(Int, row[1])] == row[3][1]) ⊻ (row[4][parse(Int, row[2])] == row[3][1])), rows))


text:
@verilog function wa(x::Wire{1:0v})

  wirearray = Vector{Wire{1:0v}}(4)

  for idx = 1:4
    wirearray[idx] = x
  end

  result = Wire([wirearray[idx] for idx in 1:4]...)

end

@test wa() == """
module wa(
  input [1:0] x,
  output [7:0] result);

  wire [1:0] wirearray[3:0];

  assign wirearray[0] = x;
  assign wirearray[1] = x;
  assign wirearray[2] = x;
  assign wirearray[3] = x;
  assign result = {wirearray[0],wirearray[1],wirearray[2],wirearray[3]};
endmodule
"""

@verilog function wa2(x::Wire{1:0v})

  wirearray = Vector{Wire{1:0v}}(4)

  wirearray[1] = Wire(0x3, 2)

  for idx = 2:4
    wirearray[idx] = x
  end

  result = Wire([wirearray[idx] for idx in 1:4]...)

end

@test wa2() == """
module wa2(
  input [1:0] x,
  output [7:0] result);

  wire [1:0] wirearray[3:0];

  assign wirearray[0] = 2'b11;
  assign wirearray[1] = x;
  assign wirearray[2] = x;
  assign wirearray[3] = x;
  assign result = {wirearray[0],wirearray[1],wirearray[2],wirearray[3]};
endmodule
"""

@verilog function flipme(x::Wire{1:0v})
  flipped = Wire(x[0], x[1])
end

@verilog function wa3(x::Wire{1:0v})

  wirearray = Vector{Wire{1:0v}}(4)

  wirearray[1] = x
  wirearray[2] = flipme(wirearray[1])

  for idx = 3:4
    wirearray[idx] = x
  end

  result = Wire([wirearray[idx] for idx in 1:4]...)

end

@test wa3() == """
module wa3(
  input [1:0] x,
  output [7:0] result);

  wire [1:0] wirearray[3:0];

  flipme flipme_wirearray_1(
    .x (wirearray[0]),
    .flipped (wirearray[1]));

  assign wirearray[0] = x;
  assign wirearray[2] = x;
  assign wirearray[3] = x;
  assign result = {wirearray[0],wirearray[1],wirearray[2],wirearray[3]};
endmodule
"""

@verilog function twooutputs(inp::Wire{1:0v})
  andor = Wire(inp[1] & inp[0], inp[1] | inp[0])
  xandor = andor ^ inp
  (andor, xandor)
end

@verilog function wa4(x::Wire{1:0v})

  wirearray = Vector{Wire{1:0v}}(4)

  wirearray[1] = x
  (wirearray[2], wirearray[3]) = twooutputs(wirearray[1])
  wirearray[4] = x

  result = Wire([wirearray[idx] for idx in 1:4]...)

end

@test wa4() == """
module wa4(
  input [1:0] x,
  output [7:0] result);

  wire [1:0] wirearray[3:0];

  twooutputs twooutputs_wirearray_1_wirearray_2(
    .inp (wirearray[0]),
    .andor (wirearray[1]),
    .xandor (wirearray[2]));

  assign wirearray[0] = x;
  assign wirearray[3] = x;
  assign result = {wirearray[0],wirearray[1],wirearray[2],wirearray[3]};
endmodule
"""


text:
function addN(n::Number)::Function
  adder(x::Number) = n + x
  return adder
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "type_morespecific" begin
    @testset "subtype" begin
        @test type_morespecific(Int, Integer)
        @test !type_morespecific(Integer, Int)
        @test !type_morespecific(Int, Int)
    end

    @testset "number of arguments" begin
        @test type_morespecific(Tuple{Int,AbstractFloat}, Tuple{Integer})
    end

    @testset "ambiguous" begin
        @test !type_morespecific(Tuple{Int,AbstractFloat}, Tuple{Integer,Float64})
        @test !type_morespecific(Tuple{Integer,Float64}, Tuple{Int,AbstractFloat})
    end

    @testset "length" begin
        @test type_morespecific(Tuple{Int,Int}, Tuple{Int,Vararg{Int}})
    end

    @testset "union" begin
        @test type_morespecific(Tuple{Int}, Tuple{Union{Int,Float64}})
        @test type_morespecific(Tuple{Union{Int,Float64}}, Tuple{Integer})
    end

    @testset "union all" begin
        @test type_morespecific(Tuple{Int}, Tuple{T} where T <: Integer)

        @test !type_morespecific(Tuple{Integer}, Tuple{T} where T <: Integer)
        @test !type_morespecific(Tuple{T} where T <: Integer, Tuple{Integer})

        @test type_morespecific(Tuple{Type{Integer}}, Tuple{Type{T}} where T <: Integer)
        @test !type_morespecific(Tuple{Type{T}} where T <: Integer, Tuple{Type{Integer}})
    end
end

@testset "anonymous_signature" begin
    @testset "diagonal dispatch" begin
        diag(::T, ::T) where T <: Integer = nothing

        diag_method = first(methods(diag, (Int, Int)))
        diag_sig = Tuple{typeof(diag),T,T} where T <: Integer
        @test diag_method.sig == diag_sig

        @test anonymous_signature(Tuple{typeof(diag),Int,Int}) == Tuple{Int,Int}
        @test anonymous_signature(diag_sig) == Tuple{T,T} where T <: Integer
        @test anonymous_signature(diag_method) == Tuple{T,T} where T <:Integer
    end

    @testset "triangular dispatch" begin
        tri(::Vector{T}, ::S) where {T, S<:T} = nothing

        tri_method = first(methods(tri, (Vector, Int)))
        tri_sig = Tuple{typeof(tri),Vector{T},S} where {T, S<:T}
        @test tri_method.sig == tri_sig

        @test anonymous_signature(Tuple{typeof(tri),Vector,Int}) == Tuple{Vector,Int}
        @test anonymous_signature(tri_sig) == Tuple{Vector{T},S} where {T, S<:T}
        @test anonymous_signature(tri_method) == Tuple{Vector{T},S} where {T, S<:T}
    end
end

@testset "anon_morespecific" begin
    plus_method = first(methods(+, (Int, Int)))
    max_method = first(methods(max, (Int, Int)))

    @test !type_morespecific(plus_method.sig, max_method.sig)
    @test anon_morespecific(plus_method, max_method)
end

@testset "dispatch" begin
    @testset "no functions" begin
        m, f = dispatch(Function[])
        @test f === nothing
        @test m === nothing
    end

    @testset "more specific function" begin
        funcs = [
            (::Int) -> 1
            (::Integer) -> 2
        ]

        m, f = dispatch(funcs, zero(Int))
        @test f == funcs[1]
        @test m.sig == Tuple{typeof(funcs[1]), Int}

        m, f = dispatch(funcs, zero(UInt))
        @test f == funcs[2]
        @test m.sig == Tuple{typeof(funcs[2]), Integer}
    end

    @testset "more specific method" begin
        f1(::Number) = 1
        f1(::Int) = 2
        f2(::Any) = 3
        f2(::Integer) = 4
        funcs = [f1, f2]

        m, f = dispatch(funcs, zero(Int))
        @test f == funcs[1]
        @test m.sig == Tuple{typeof(funcs[1]), Int}

        m, f = dispatch(funcs, zero(UInt))
        @test f == funcs[2]
        @test m.sig == Tuple{typeof(funcs[2]), Integer}
    end

    @testset "ambiguous" begin
        funcs = [
            (::Integer) -> 1
            (::Integer) -> 2
        ]

        # In ambiguous cases the last function is preferred
        m, f = dispatch(funcs, 0)
        @test f == funcs[2]
        @test m.sig == Tuple{typeof(funcs[2]), Integer}

        m, f = dispatch(reverse(funcs), 0)
        @test f == funcs[1]
        @test m.sig == Tuple{typeof(funcs[1]), Integer}
    end

    @testset "optional" begin
        funcs = [
            (x, y=0) -> 1
            (x) -> 2
        ]

        m, f = dispatch(funcs, 0)
        @test f == funcs[2]
        @test m.sig == Tuple{typeof(funcs[2]), Any}

        m, f = dispatch(funcs, 0, 0)
        @test f == funcs[1]
        @test m.sig == Tuple{typeof(funcs[1]), Any, Any}

        m, f = dispatch(reverse(funcs), 0)
        @test f == funcs[1]
        @test m.sig == Tuple{typeof(funcs[1]), Any}
    end

    @testset "Type{T}" begin
        funcs = [
            (::Type{Integer}) -> 1
            (::Type{Int}) -> 2
        ]

        m, f = dispatch(funcs, 0)
        @test f === nothing
        @test m === nothing

        m, f = dispatch(funcs, Integer)
        @test f == funcs[1]
        @test m.sig == Tuple{typeof(funcs[1]), Type{Integer}}

        m, f = dispatch(reverse(funcs), Type{Integer})
        @test f === nothing
        @test m === nothing
    end

    @testset "Type{T} where T" begin
        funcs = [
            (::Type{Int}) -> 1
            (::Type{<:Integer}) -> 2
        ]

        m, f = dispatch(funcs, Int)
        @test f == funcs[1]
        @test m.sig == Tuple{typeof(funcs[1]), Type{Int}}

        m, f = dispatch(reverse(funcs), Integer)
        @test f == funcs[2]
        @test m.sig == Tuple{typeof(funcs[2]), Type{<:Integer}}
    end
end


text:
# ------------------------------------------------------------------------------
# lumped BDDC on macro elements example
# ------------------------------------------------------------------------------

using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
p = 1
numberelements1d = 4

# operator
diffusion = GalleryMacroElementOperator("diffusion", p + 1, p + 2, numberelements1d, mesh)

# BDDC smoother
bddc = DirichletBDDC(diffusion)

# compute operator symbols
A = computesymbols(bddc, [0.2], [π, π])
eigenvalues = real(eigvals(A))

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


text:
# ---
# title: 1637. Widest Vertical Area Between Two Points Containing No Points
# id: problem1637
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Sort
# link: <https://leetcode.com/problems/widest-vertical-area-between-two-points-containing-no-points/description/>
# hidden: true
# ---
#
# Given `n` `points` on a 2D plane where `points[i] = [xi, yi]`, Return _  the
# **widest vertical area** between two points such that no points are inside the
# area._
#
# A **vertical area** is an area of fixed-width extending infinitely along the
# y-axis (i.e., infinite height). The **widest vertical area** is the one with
# the maximum width.
#
# Note that points **on the edge** of a vertical area **are not** considered
# included in the area.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/19/points3.png)
#
#
#
#     Input: points = [[8,7],[9,9],[7,4],[9,7]]
#     Output: 1
#     Explanation: Both the red and the blue area are optimal.
#
#
# **Example 2:**
#
#
#
#     Input: points = [[3,1],[9,0],[1,0],[1,4],[5,3],[8,8]]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * `n == points.length`
#   * `2 <= n <= 105`
#   * `points[i].length == 2`
#   * `0 <= xi, yi <= 109`
#
#
## @lc code=start
using LeetCode

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


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ af957aaa-1cb9-11eb-1eb9-6bf4907a1628
using LinearAlgebra

# ╔═╡ 2572e992-1cba-11eb-3315-1f132089db17
A = [10 2;3 4]

# ╔═╡ 35f6a224-1cba-11eb-0e07-c9c4b43619b8
det(A)

# ╔═╡ 5ef57da4-1cba-11eb-12ae-8bd131d32935
tr(A)

# ╔═╡ 9bc6985a-1cba-11eb-3836-c314153cc1cd
inv(A)

# ╔═╡ ad81fcd8-1cba-11eb-2a91-af2cfbc9cc73
eigvals(A)

# ╔═╡ db9c32d2-1cba-11eb-36fd-4db5d1ed8a9c
eigvecs(A)

# ╔═╡ 97dfd928-1cbb-11eb-2858-2f5a079ee685
factorize(A)

# ╔═╡ d36a84b4-1cbb-11eb-1884-b7165e7741b3
Symmetric(A)

# ╔═╡ eba461d0-1cbb-11eb-042a-6f6d85b8f6e6
Diagonal(A)

# ╔═╡ f59c977a-1cbb-11eb-1a8a-17f77b281e40
Tridiagonal(A)

# ╔═╡ 1387ed2a-1cbc-11eb-1d9b-ed159167d81f
svd(A)

# ╔═╡ 6ffbd562-1cbc-11eb-37ab-d368cb70fe31
qr(A)

# ╔═╡ 86c0ffa2-1cbc-11eb-3f90-8fe517e4a514
eigen(A)

# ╔═╡ f3f9a022-1cbe-11eb-37b9-392a6bd21c8d
rank(A)

# ╔═╡ fd489688-1cbe-11eb-19f0-83313a3c4f0d
A*A

# ╔═╡ 3e1faae8-1cbf-11eb-11b7-a7f3ce3b43c1
A.*A

# ╔═╡ Cell order:
# ╠═af957aaa-1cb9-11eb-1eb9-6bf4907a1628
# ╠═2572e992-1cba-11eb-3315-1f132089db17
# ╠═35f6a224-1cba-11eb-0e07-c9c4b43619b8
# ╠═5ef57da4-1cba-11eb-12ae-8bd131d32935
# ╠═9bc6985a-1cba-11eb-3836-c314153cc1cd
# ╠═ad81fcd8-1cba-11eb-2a91-af2cfbc9cc73
# ╠═db9c32d2-1cba-11eb-36fd-4db5d1ed8a9c
# ╠═97dfd928-1cbb-11eb-2858-2f5a079ee685
# ╠═d36a84b4-1cbb-11eb-1884-b7165e7741b3
# ╠═eba461d0-1cbb-11eb-042a-6f6d85b8f6e6
# ╠═f59c977a-1cbb-11eb-1a8a-17f77b281e40
# ╠═1387ed2a-1cbc-11eb-1d9b-ed159167d81f
# ╠═6ffbd562-1cbc-11eb-37ab-d368cb70fe31
# ╠═86c0ffa2-1cbc-11eb-3f90-8fe517e4a514
# ╠═f3f9a022-1cbe-11eb-37b9-392a6bd21c8d
# ╠═fd489688-1cbe-11eb-19f0-83313a3c4f0d
# ╠═3e1faae8-1cbf-11eb-11b7-a7f3ce3b43c1


text:
# Generate scales.

import JSON

# This ID is incremented for every problem generated.  Barring changes
# to the generation, this should generate the same ids on subsequent
# runs.
cur_id = 1

function emit(question, answer)
    global cur_id
    println("INSERT INTO \"probs\" VALUES($cur_id, '$question', '$answer');")
    cur_id += 1
end

key_names = [ "C", "G", "D", "A", "E", "B", "F♯", "G♭", "D♭", "A♭", "E♭", "B♭", "F" ]

function gen_problem(desc)
    for key in key_names
        emit("$key $desc", "played")
    end
end

println("DELETE FROM probs;")

# Learn symmetric dominant over ii-V-I major.
gen_problem("sym dom over ii-V-I")


text:
@testset "Running WBS tests..." begin

    n = 1000;        # Number of samples
    λ = 100;         # Frequencey of changes

    ########################
    # Normal mean segments #
    ########################
    μ, σ = Normal(0.0, 10.0), 1.0
    sample, cps = @changepoint_sampler n λ Normal(μ, σ)
    seg_cost_CUSUM = Changepoints.CUSUM(sample);
    seg_cost_sSIC = Changepoints.sSIC(sample);

    @test_nowarn WBS(seg_cost_CUSUM, n);
    @test_nowarn get_WBS_changepoints(seg_cost_sSIC, WBS(seg_cost_CUSUM, n))
    @test_nowarn WBS(seg_cost_CUSUM, n, 1.3, 1.0, 5000, true , 1/sqrt(2));

end


text:
"""
```
impulse_responses(m, system)

impulse_responses(system, horizon)
```

Compute impulse responses for a single draw.

### Inputs

- `m::AbstractDSGEModel`: model object
- `system::System{S}`: state-space system matrices
- `horizon::Int`: number of periods ahead to forecast
- `flip_shocks::Bool`: Whether to compute IRFs in response to a positive shock (by default the shock magnitude is a negative 1 std. shock)

where `S<:AbstractFloat`

### Outputs

- `states::Array{S, 3}`: matrix of size `nstates` x `horizon` x `nshocks` of
  state impulse response functions
- `obs::Array{S, 3}`: matrix of size `nobs` x `horizon` x `nshocks` of
  observable impulse response functions
- `pseudo::Array{S, 3}`: matrix of size `npseudo` x `horizon` x `nshocks` of
  pseudo-observable impulse response functions
"""
function impulse_responses(m::AbstractRepModel,
                           system::Union{System{S}, RegimeSwitchingSystem{S}};
                           flip_shocks::Bool = false) where {S<:AbstractFloat}
    horizon = impulse_response_horizons(m)
    if typeof(system) == RegimeSwitchingSystem{S}
        last_reg = get_setting(m, :n_regimes)
        states, obs, pseudo = impulse_responses(system[last_reg], horizon, flip_shocks = flip_shocks)
    else
        states, obs, pseudo = impulse_responses(system, horizon, flip_shocks = flip_shocks)
    end

    return states, obs, pseudo
end

function impulse_responses_augmented(m::Union{AbstractHetModel,RepDSGEGovDebt}, system::System{S};
                                     flip_shocks::Bool = false,
                                     use_alternate_consumption = false) where {S<:AbstractFloat}
    horizon = impulse_response_horizons(m)

    TTT_jump, TTT_state, eu  = klein(m)
    if eu == -1
        throw(KleinError())
    end

    TTT, RRR = klein_transition_matrices(m, TTT_state, TTT_jump)
    CCC      = zeros(n_model_states(m))

    TTT, RRR, CCC = #if typeof(m) != RepDSGEGovDebt{Float64}
        augment_states(m, TTT, RRR, CCC)
   #= else
        TTT, RRR, CCC
    end=#

    measurement_equation = measurement(m, TTT, RRR, CCC)
    transition_equation = Transition(TTT, RRR, CCC)
    system = System(transition_equation, measurement_equation)

    states, obs, pseudo = impulse_responses(system, horizon, flip_shocks = flip_shocks)

    if typeof(m) == RepDSGEGovDebt{Float64}
        return states, obs, pseudo
    end

    state_indices_orig = stack_indices(m.endogenous_states_original, get_setting(m, :states))
    jump_indices_orig  = stack_indices(m.endogenous_states_original, get_setting(m, :jumps))

    Qx = get_setting(m, :Qx)
    Qy = get_setting(m, :Qy)

    # In this case, the length of state_indices and jump_indices seems to give the number of
    # UNNORMALIZED states/jumps however I'm not sure if this will always be the case/if this is
    # really what we want saved into these settings...so might need to modify in future.
    states_unnormalized = Array{Float64}(undef, length(state_indices_orig),
                                         horizon, size(states, 3))
    jumps_unnormalized = Array{Float64}(undef, length(jump_indices_orig), horizon, size(states, 3))

    model_states_unnormalized = Array{Float64}(undef, length(state_indices_orig) +
                                               length(jump_indices_orig) +
                                               length(m.endogenous_states_augmented), horizon,
                                               size(states, 3))
    for i in 1:size(states, 3)
        state_inds = 1:n_backward_looking_states(m)
        states_unnormalized[:, :, i] = Matrix(Qx')*states[state_inds, :, i]
        jump_inds = n_backward_looking_states(m)+1:n_backward_looking_states(m) + n_jumps(m)
        jumps_unnormalized[:, :, i] = Matrix(Qy')*states[jump_inds, :, i]

        augmented_states = states[n_backward_looking_states(m) + n_jumps(m) + 1 : end, :, i]
        model_states_unnormalized[:, :, i] = vcat(states_unnormalized[:, :, i],
                                                  jumps_unnormalized[:, :, i],
                                                  augmented_states)
    end

    if use_alternate_consumption
        endo = m.endogenous_states_original
        c_implied = (m[:ystar]/m[:g]) - m[:xstar]
        IRFC_implied = m[:ystar]/(c_implied*m[:g])*(model_states_unnormalized[endo[:y′_t], :, :] -
                                                    model_states_unnormalized[endo[:g′_t], :, :])-
                                                    (m[:xstar]/c_implied)*
        model_states_unnormalized[endo[:I′_t], :, :]

        model_states_unnormalized[1229, :, :] = IRFC_implied

        z_consumption = model_states_unnormalized[endo[:z′_t], :, :]
        g_lag = cat(0, model_states_unnormalized[endo[:g′_t], :, :], dims = 2)[:, 1:40, :]
        IRFC_implied_lag = m[:ystar]/(c_implied*m[:g])*(model_states_unnormalized[endo[:y′_t1], :, :]  - g_lag) - (m[:xstar]/c_implied)*model_states_unnormalized[endo[:I′_t1], :, :]
        IRF_observable_implied = IRFC_implied - IRFC_implied_lag .+ z_consumption
        obs[m.observables[:obs_consumption], :, :] = IRF_observable_implied
    end

    return model_states_unnormalized, obs, pseudo
end

function impulse_responses(m::Union{AbstractHetModel,RepDSGEGovDebt}, system::System{S};
                           flip_shocks::Bool = false,
                           use_augmented_states::Bool = true,
                           use_alternate_consumption::Bool = false) where {S<:AbstractFloat}
    if use_augmented_states
        return impulse_responses_augmented(m, system, flip_shocks = flip_shocks,
                                           use_alternate_consumption = use_alternate_consumption)
    end

    horizon = impulse_response_horizons(m)
    states, obs, pseudo = impulse_responses(system, horizon, flip_shocks = flip_shocks)

    if typeof(m) == RepDSGEGovDebt
        return states, obs, pseudo
    end

    Qx = get_setting(m, :Qx)
    Qy = get_setting(m, :Qy)

    state_indices_orig = stack_indices(m.endogenous_states_original, get_setting(m, :states))
    jump_indices_orig  = stack_indices(m.endogenous_states_original, get_setting(m, :jumps))

    # In this case, the length of state_indices and jump_indices seems to give the number of
    # UNNORMALIZED states/jumps however I'm not sure if this will always be the case/if this is
    # really what we want saved into these settings...so might need to modify in future.
    states_unnormalized = Array{Float64}(undef, length(state_indices_orig),
                                         horizon, size(states, 3))
    jumps_unnormalized = Array{Float64}(undef, length(jump_indices_orig), horizon, size(states, 3))
    model_states_unnormalized = Array{Float64}(undef, length(state_indices_orig) +
                                               length(jump_indices_orig), horizon, size(states, 3))
    for i in 1:size(states, 3)
        state_inds = 1:n_backward_looking_states(m)
        states_unnormalized[:, :, i] = Matrix(Qx')*states[state_inds, :, i]
        jump_inds = n_backward_looking_states(m)+1:n_backward_looking_states(m) + n_jumps(m)
        jumps_unnormalized[:, :, i] = Matrix(Qy')*states[jump_inds, :, i]

        model_states_unnormalized[:, :, i] = vcat(states_unnormalized[:, :, i],
                                                  jumps_unnormalized[:, :, i])
    end

    if use_alternate_consumption
        endo_orig = m.endogenous_states_original
        c_implied = (m[:ystar]/m[:g]) - m[:xstar]
        IRFC_implied = m[:ystar]/(c_implied*m[:g])*(model_states_unnormalized[endo_orig[:y′_t], :, :]  - model_states_unnormalized[endo_orig[:g′_t], :, :]) - (m[:xstar]/c_implied)*model_states_unnormalized[endo_orig[:I′_t], :, :]
        model_states_unnormalized = vcat(model_states_unnormalized, IRFC_implied)

        z_consumption = model_states_unnormalized[m.endogenous_states_original[:z′_t], :, :]
        IRFC_implied_lag = cat(0, IRFC_implied, dims = 2)[:, 1:40, :]
        IRF_observable_implied = IRFC_implied - IRFC_implied_lag .+ z_consumption
        obs[m.observables[:obs_consumption], :, :] = IRF_observable_implied
    end

    return model_states_unnormalized, obs, pseudo
end

function impulse_responses(system::System{S}, horizon::Int;
                           flip_shocks::Bool = false) where {S<:AbstractFloat}
    # Setup
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    for i = 1:nshocks
        # Isolate single shock
        shocks = zeros(S, nshocks, horizon)
        if flip_shocks
            shocks[i, 1] = sqrt(system[:QQ][i, i]) # a positive 1 s.d. shock
        else
            shocks[i, 1] = -sqrt(system[:QQ][i, i]) # a negative 1 s.d. shock
        end
        # Iterate state space forward
        states[:, :, i], obs[:, :, i], pseudo[:, :, i], _ = forecast(system, s_0, shocks)
    end

    return states, obs, pseudo
end

# Method for specifying the subset of shocks, and the size of each shock
function impulse_responses(m::AbstractDSGEModel, system::System{S},
                           horizon::Int, shock_names::Vector{Symbol},
                           shock_values::Vector{Float64}) where S<:AbstractFloat

    # Must provide a name and value for each shock
    @assert length(shock_names) == length(shock_values)

    # Setup
    exo          = m.exogenous_shocks
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    for (i, shock) in enumerate(shock_names)
        # Isolate single shock
        shocks = zeros(S, nshocks, horizon)
        shocks[exo[shock], 1] = shock_values[i]

        # Iterate state space forward
        states[:, :, exo[shock]], obs[:, :, exo[shock]], pseudo[:, :, exo[shock]], _ = forecast(system, s_0, shocks)
    end

    return states, obs, pseudo
end

# Method for picking a specific shock and the size of the desired shift in the state or
# observed variable using that shock. (Omitting the feature to do pseudo-observables
# now, for simplicity)
function impulse_responses(m::AbstractDSGEModel, system::System{S},
                           horizon::Int, shock_name::Symbol,
                           var_name::Symbol, var_value::Float64) where S<:AbstractFloat
    # Setup
    var_names, var_class =
    if var_name in keys(m.endogenous_states)
        m.endogenous_states, :states
    elseif var_name in keys(m.observables)
        m.observables, :obs
    end
    exo          = m.exogenous_shocks
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    states = zeros(S, nstates, horizon, nshocks)
    obs    = zeros(S, nobs,    horizon, nshocks)
    pseudo = zeros(S, npseudo, horizon, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    # Isolate single shock
    shocks = zeros(S, nshocks, horizon)
    if var_class == :states
        shocks[exo[shock_name], 1] = obtain_shock_from_desired_state_value(var_value,
                                                                      var_names[var_name],
                                                                      exo[shock_name],
                                                                      system[:RRR])
    else # == :obs
        shocks[exo[shock_name], 1] = obtain_shock_from_desired_obs_value(var_value,
                                                                    var_names[var_name],
                                                                    exo[shock_name],
                                                                    system[:ZZ],
                                                                    system[:RRR])
    end

    # Iterate state space forward
    states[:, :, exo[shock_name]], obs[:, :, exo[shock_name]], pseudo[:, :, exo[shock_name]], _ = forecast(system, s_0, shocks)

    return states, obs, pseudo
end

"""
```
function impulse_responses(system::System{S}, horizon::Int, permute_mat::AbstractMatrix{T},
                           shocks::AbstractVector{S} = Vector{S}(undef, 0);
                           restriction::Symbol = :short_run, flip_shocks::Bool = false,
                           get_shocks::Bool = false) where {S <: Real, T <: Number}
```
computes impulse responses using a Cholesky-identified shock
to observables, with the ability to apply a permutation matrix.

### Inputs
- `system::System{S}`: state-space system matrices
- `horizon::Int`: number of periods ahead to forecast
- `permute_mat::AbstractMatrix{T}`: permutation matrix
- `shocks::AbstractVector{S}`: vector of orthogonal shocks to observables.
    See the section Restriction Types below for a mathematical explanation.
    If no vector is provided, then we assume the user wants a shock vector
    whose first entry is one and all other entries are zero (a 1 standard deviation
    shock to the innovation affecting the first variable).
- `restriction::Symbol`: type of restriction for the Cholesky identification.
    Can be either `:short_run` or `:long_run`. We also let `:cholesky` refer
    specifically to the `:short_run` restriction and
    `:choleskyLR` or `:cholesky_long_run` refer to the `:long_run` restriction.
- `flip_shocks::Bool`: Whether to compute IRFs in response to a positive shock
    (by default the shock magnitude is a negative 1 std. shock)
- `get_shocks::Bool`: Whether to return the Cholesky-identified shocks and the
    underlying structural shocks

### Outputs

- `states::Matrix{S}`: matrix of size `nstates` x `horizon` x `nshocks` of
  state impulse response functions
- `obs::Matrix{S}`: matrix of size `nobs` x `horizon` x `nshocks` of
  observable impulse response functions
- `pseudo::Matrix{S}`: matrix of size `npseudo` x `horizon` x `nshocks` of
  pseudo-observable impulse response functions
- `cholesky_shock::Vector{S}`: Cholesky-identified shock on impact (only if `get_shocks = true`)
- `structural_shock::Vector{S}`: structural shocks causing the
    Cholesky-identified shock (only if `get_shocks = true`)

### Restriction Types
Consider a state space system
```
sₜ = TTT * sₜ₋₁ + RRR * sqrt(QQ) * ϵₜ
yₜ = ZZ * sₜ
```
where `ϵₜ ∼ 𝒩 (0, I)` and `QQ` is a diagonal matrix specifying the
variances of the structural shocks `ϵₜ`. Under this specification,
the units of `ϵₜ` are standard deviations.

A Cholesky identification searches for a lower triangular matrix `M` such that `M * M'` equals the
covariance of observables. The restriction is short-run if we take the short-run covariance and is long
run if we take the long-run covariance (i.e. the covariance of observables when shocks occur to the
stationary level of sₜ). Explicitly, the short- and long-run restrictions are, respectively,
```
Covₛᵣ = (ZZ * RRR) * QQ * (ZZ * RRR)',
Covₗᵣ = (ZZ * (I - TTT)⁻¹ * RRR) * QQ * (ZZ * (I - TTT)⁻¹ * RRR)'.
```
Taking the Cholesky decomposition of the appropriate covariance matrix yields the
lower triangular matrix `M`.

Using `M`, we may re-write the measurement equation of the state space as
```
yₜ = ZZ * TTT * sₜ₋₁ + M * uₜ,
```
where `uₜ ∼ 𝒩 (0, I)` are orthogonal shocks identified from the covariance of `yₜ`
using `M`. The orthogonal shocks `uₜ` are precisely the input argument `shocks`.
We can map `uₜ` to `ϵₜ`, i.e. identify structural shocks, by solving the linear system
```
(ZZ * RRR * sqrt(QQ)) * ϵₜ = M * uₜ
```
when using the short-run restriction and
```
(ZZ * (I - TTT)⁻¹ * RRR * sqrt(QQ)) * ϵₜ = M * uₜ
```
when using the long-run restriction.
"""
function impulse_responses(system::System{S}, horizon::Int, permute_mat::AbstractMatrix{T},
                           shocks::AbstractVector{S} = Vector{S}(undef, 0);
                           restriction::Symbol = :short_run, flip_shocks::Bool = false,
                           get_shocks::Bool = false) where {S <: Real, T <: Number}
    if !(restriction in [:short_run, :long_run, :cholesky, :choleskyLR, :cholesky_long_run])
        error("The restriction $restriction has not been implemented. "
              * "Choose either `:short_run` or `:long_run")
    end
    if restriction == :cholesky
        restriction = :short_run
    elseif restriction in [:choleskyLR, :cholesky_long_run]
        restriction = :long_run
    end

    # Permute matrices
    if !(permute_mat ≈ Matrix{eltype(permute_mat)}(I, size(permute_mat)...))
        do_permute = true
        ZZold = copy(system[:ZZ])
        system.measurement.ZZ = permute_mat * system[:ZZ]
    else
        do_permute = false
    end

    # Set up IRFs
    nshocks = size(system[:RRR], 2)
    nstates = size(system[:TTT], 1)
    nobs    = size(system[:ZZ], 1)
    npseudo = size(system[:ZZ_pseudo], 1)

    system = DSGE.zero_system_constants(system) # Set constant system matrices to 0
    s₀ = zeros(S, nstates)

    # Back out structural shocks corresponding to Cholesky shock to observables
    if isempty(shocks)
        shocks = vcat(1., zeros(nobs - 1)) # First orthogonal shock to observables
    end

    obs_std = if restriction == :short_run
        system[:ZZ] * system[:RRR] * sqrt.(system[:QQ])
    elseif restriction == :long_run
        system[:ZZ] * inv(Matrix{S}(I, nstates, nstates) - system[:TTT]) *
            system[:RRR] * sqrt.(system[:QQ])
    end

    obs_cov = obs_std * obs_std'
    cholmat = try
        cholesky(obs_cov).L
    catch e
        if isa(e, PosDefException)
            cholesky((obs_cov + obs_cov') ./ 2).L
        else
            rethrow(e)
        end
    end
    structural_shock = obs_std \ (cholmat * shocks)

    # Run IRFs
    shocks_augmented = zeros(S, nshocks, horizon)
    shocks_augmented[:, 1] = flip_shocks ? structural_shock : -structural_shock # negative shock by default
    states, obs, pseudo = forecast(system, s₀, shocks_augmented)

    # Restore system if permutation matrix is not the identity
    if do_permute
        system.measurement.ZZ = ZZold
    end

    if get_shocks
        return states, obs, pseudo, structural_shock
    else
        return states, obs, pseudo
    end
end

"""
```
function impulse_responses(system::System{S}, horizon::Int, frequency_band::Tuple{S,S},
                           n_obs_shock::Int; flip_shocks::Bool = false,
                           get_shocks::Bool = false) where {S<:Real}
```
computes impulse responses by identifying the shock which
maximizes the variance explained of a chosen observable within
a certain frequency band.

For typical business-cycle frequences, we recommend setting
`frequency_band = (2 * π / 32, 2 * π / 6)`.

### Inputs
- `system::System{S}`: state-space system matrices
- `horizon::Int`: number of periods ahead to forecast
- `frequency_band::Tuple{S, S}`: the frequencies at which we want to maximize the variance explained.
    The first frequency must be less than the second one.
- `n_obs_shock::Int`: the index of the observable whose variance we want to maximally explain
    in the state space model implied by `system`.
- `flip_shocks::Bool`: Whether to compute IRFs in response to a positive shock (by default the shock magnitude is a negative 1 std. shock)
- `get_shocks::Bool`: Whether to return the structural shocks.

where `S <: Real`

### Outputs
- `states::Matrix{S}`: matrix of size `nstates` x `horizon` x `nshocks` of
  state impulse response functions
- `obs::Matrix{S}`: matrix of size `nobs` x `horizon` x `nshocks` of
  observable impulse response functions
- `pseudo::Matrix{S}`: matrix of size `npseudo` x `horizon` x `nshocks` of
  pseudo-observable impulse response functions
- `structural_shock::Vector{S}`: structural shocks causing the
    Cholesky-identified shock (only if `get_shocks = true`)
"""
function impulse_responses(system::System{S}, horizon::Int, frequency_band::Tuple{S,S},
                           n_obs_shock::Int; flip_shocks::Bool = false,
                           get_shocks::Bool = false) where {S<:Real}
    if frequency_band[1] > frequency_band[2]
        @warn "First frequency in `frequency_band` is larger than second one. Switching the order of the frequencies."
        frequency_band = (frequency_band[2], frequency_band1[1])
    end

    # Setup
    nshocks = size(system[:RRR], 2)
    nstates = size(system[:TTT], 1)
    nobs    = size(system[:ZZ], 1)
    npseudo = size(system[:ZZ_pseudo], 1)
    system  = zero_system_constants(system)

    # Compute business-cycle variance maximizing shock
    increment = abs(frequency_band[1] - frequency_band[2]) / 200.
    variance = zeros(S, nshocks, nshocks)
    sd_shock = system[:RRR] * sqrt.(system[:QQ]) # one standard deviation shock
    I_state = Matrix{S}(I, nstates, nstates)
    for f = frequency_band[1]:increment:round(frequency_band[2], digits=10) # not rounding sometimes leads to one fewer loop than desired
        sumG = system[:TTT] .* exp(-im * f)
        invA = system[:ZZ] * ((I_state - sumG) \ sd_shock)
        variance += real.(invA[n_obs_shock,:] * invA[n_obs_shock,:]') .* increment ./
            abs(frequency_band[1] - frequency_band[2])
    end
    eigout = eigen(variance)
    q = eigout.vectors[:, argmax(eigout.values)]
    q .*= sign(q[n_obs_shock])

    # Compute IRFs
    states = zeros(S, nstates, horizon)
    shocks_augmented      = zeros(S, nshocks, horizon)
    shocks_augmented[:,1] = flip_shocks ? (sqrt.(system[:QQ]) * q) : -(sqrt.(system[:QQ]) * q) # q -> positive shock. To stick with convention, "flip_shock" should yield a positive shock.
    states, obs, pseudo = forecast(system, zeros(S, nstates), shocks_augmented)

    if get_shocks
        return states, obs, pseudo, vec(shocks_augmented[:,1])
    else
        return states, obs, pseudo
    end
end

function impulse_responses_peg(m::AbstractRepModel, system::System{S}, horizon::Int;
                           flip_shocks::Bool = false, H::Int = 0, peg::Symbol = :all_periods, real_rate = false) where {S<:AbstractFloat}
    ## H is peg horizon

    # Setup
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 1)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    shocks = m.exogenous_shocks

    states = zeros(S, nstates, horizon) #, nshocks)
    obs    = zeros(S, nobs,    horizon) #, nshocks)
    pseudo = zeros(S, npseudo, horizon) #, nshocks)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    s_0 = zeros(S, nstates)

    ### IRFs to FFR peg for H periods
    # Set anticipated shocks to peg FFR from t+1 to t+H
    FFRpeg = -0.25/4 #/400

    PsiR1 = 0 #constant
    # ZZ matrix
    PsiR2 = zeros(nstates)
    PsiR2[m.endogenous_states[:R_t]] = 1
    if real_rate
        PsiR2[m.endogenous_states[:Eπ_t]] = -1
    end
    if H == 0
        Rht = system[:RRR][:,shocks[:rm_sh]]
    else
        Rht = system[:RRR][:,vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")])] #% Columns of Rh referring to anticipated monetary shocks
    end
    bb = zeros(H+1,1)
    MH = zeros(H+1,H+1)
    for hh = 1:H+1
        if peg == :all_periods
           # @show size(FFRpeg), size(PsiR1), size(PsiR2), size(system[:TTT]), hh, size(s_0)
            #@show size(FFRpeg - PsiR1 - PsiR2'*(system[:TTT])^hh*s_0)
            bb[hh,1] = (FFRpeg - PsiR1 - PsiR2'*(system[:TTT])^hh*s_0)
            #@show bb[hh, 1]
            MH[hh,:] = PsiR2'*(system[:TTT])^(hh-1)*Rht
        elseif peg == :some_periods
            if hh < H+1
                FFRpeg_ = 0
            else
                FFRpeg_ = FFRpeg
            end
            bb[hh,1] = (FFRpeg_ - PsiR1 - PsiR2'*(system[:TTT])^hh*s_0)[1]
            MH[hh,:] = PsiR2'*(system[:TTT])^(hh-1)*Rht
        end

    end
    monshocks = MH\bb
    #@show size(MH), bb, monshocks
    #% Simulate the model
    etpeg = zeros(nshocks,horizon)
    if H == 0
        etpeg[shocks[:rm_sh],1] = monshocks[1]
    else
        etpeg[vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")]),1] = monshocks
    end
    states[:, :], obs[:, :], pseudo[:, :] = forecast(system, s_0, etpeg)
    return states, obs, pseudo
end

# Given a desired value, x, for some state s^i_1, back out the necessary value of
# ϵ^j_1 needed s.t. s^i_1 = x.
# E.g. if I want the impulse response of output (s^i_1) to be 50 basis points
# based on a shock to productivity (ϵ^j_1), how large does that shock
# to productivity need to be?
# The answer is ϵ^j_1 = x/R_{i,j}
function obtain_shock_from_desired_state_value(state_value::Float64, state_ind::Int,
                                               shock_ind::Int, RRR::Matrix{Float64})
    return state_value/RRR[state_ind, shock_ind]
end

# Similar principle to the function for states
# The answer is ϵ^j_1 = x/{Σ_{k=1}^n Z_{i, k} R_{k, j}},
# where n is the total number of states, i is the index of the observable of interest
# y^i_1, and j is the index of the shock of interest, ϵ^j_1
function obtain_shock_from_desired_obs_value(obs_value::Float64, obs_ind::Int,
                                             shock_ind::Int, ZZ::Matrix{Float64},
                                             RRR::Matrix{Float64})
    return obs_value/dot(ZZ[obs_ind, :], RRR[:, shock_ind])
end


text:
module RulesNormalMeanVarianceMeanTest

using Test
using ReactiveMP
using Random
using Distributions

import ReactiveMP: @test_rules

@testset "rules:NormalMeanVariance:mean" begin

    @testset "Belief Propagation: (m_out::PointMass, m_v::PointMass)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (m_out = PointMass(-1.0), m_v = PointMass(2.0)), output = NormalMeanVariance(-1.0, 2.0)),
            (input = (m_out = PointMass(1.0), m_v = PointMass(2.0)),  output = NormalMeanVariance(1.0, 2.0)),
            (input = (m_out = PointMass(2.0), m_v = PointMass(1.0)),  output = NormalMeanVariance(2.0, 1.0))
        ]

    end

    @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, m_v::PointMass)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (m_out = NormalMeanVariance(0.0, 1.0),  m_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_out = NormalMeanVariance(-1.0, 1.0), m_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_out = NormalMeanVariance(2.0, 0.5),  m_v = PointMass(1.0)), output = NormalMeanVariance(2.0, 1.5)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (m_out = NormalMeanPrecision(0.0, 1.0),  m_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_out = NormalMeanPrecision(-1.0, 1.0), m_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_out = NormalMeanPrecision(2.0, 0.5),  m_v = PointMass(1.0)), output = NormalMeanVariance(2.0, 3.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (m_out = NormalWeightedMeanPrecision(0.0, 1.0),  m_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_out = NormalWeightedMeanPrecision(-1.0, 1.0), m_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_out = NormalWeightedMeanPrecision(2.0, 0.5),  m_v = PointMass(1.0)), output = NormalMeanVariance(4.0, 3.0)),
        ]

    end

    @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, q_v::PointMass)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (m_out = NormalMeanVariance(0.0, 1.0),  q_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_out = NormalMeanVariance(-1.0, 1.0), q_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_out = NormalMeanVariance(2.0, 0.5),  q_v = PointMass(1.0)), output = NormalMeanVariance(2.0, 1.5)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (m_out = NormalMeanPrecision(0.0, 1.0),  q_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_out = NormalMeanPrecision(-1.0, 1.0), q_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_out = NormalMeanPrecision(2.0, 0.5),  q_v = PointMass(1.0)), output = NormalMeanVariance(2.0, 3.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (m_out = NormalWeightedMeanPrecision(0.0, 1.0),  q_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_out = NormalWeightedMeanPrecision(-1.0, 1.0), q_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_out = NormalWeightedMeanPrecision(2.0, 0.5),  q_v = PointMass(1.0)), output = NormalMeanVariance(4.0, 3.0)),
        ]

    end

    @testset "Variational: (q_out::Any, q_v::Any)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (q_out = PointMass(-1.0), q_v = PointMass(2.0)), output = NormalMeanVariance(-1.0, 2.0)),
            (input = (q_out = PointMass(1.0), q_v = PointMass(2.0)),  output = NormalMeanVariance(1.0, 2.0)),
            (input = (q_out = PointMass(2.0), q_v = PointMass(1.0)),  output = NormalMeanVariance(2.0, 1.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (q_out = NormalMeanVariance(-1.0, 2.0), q_v = PointMass(2.0)), output = NormalMeanVariance(-1.0, 2.0)),
            (input = (q_out = NormalMeanPrecision(1.0, 4.0), q_v = PointMass(3.0)),  output = NormalMeanVariance(1.0, 3.0)),
            (input = (q_out = NormalWeightedMeanPrecision(2.0, 4.0), q_v = PointMass(1.0)),  output = NormalMeanVariance(0.5, 1.0))
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:μ, Marginalisation) [
            (input = (q_out = PointMass(-1.0), q_v = InverseGamma(2.0, 1.0)), output = NormalMeanVariance(-1.0, 1.0)),
            (input = (q_out = PointMass(1.0), q_v = InverseGamma(4.0, 2.0)),  output = NormalMeanVariance(1.0, 2.0 / 3.0)),
            (input = (q_out = PointMass(2.0), q_v = InverseGamma(4.0, 6.0)),  output = NormalMeanVariance(2.0, 2.0))
        ]

    end

end

end


text:
"""
```
carter_kohn_smoother(y, T, R, C, Q, Z, D, E, s_0, P_0;
    Nt0 = 0, draw_states = true)

carter_kohn_smoother(regime_indices, y, Ts, Rs, Cs, Qs,
    Zs, Ds, Es, s_0, P_0; Nt0 = 0, draw_states = true)
```
This program is a simulation smoother based on Carter and Kohn's \"On Gibbs
Sampling for State Space Models\" (Biometrika, 1994). It recursively samples
from the conditional distribution of time-t states given the full set of
observables and states from time t+1 to time T. Unlike the Durbin-Koopman
simulation smoother, this one does rely on inverting potentially singular
matrices using the Moore-Penrose pseudoinverse.

The state space is augmented with shocks (see `?augment_states_with_shocks`),
and the augmented state space is Kalman filtered and smoothed. Finally, the
smoothed states and shocks are indexed out of the augmented state vectors.

The original state space (before augmenting with shocks) is given by:

```
s_{t+1} = C + T*s_t + R*ϵ_t    (transition equation)
y_t     = D + Z*s_t + u_t      (measurement equation)

ϵ_t ∼ N(0, Q)
u_t ∼ N(0, E)
Cov(ϵ_t, u_t) = 0
```

### Inputs

- `y`: `Ny` x `Nt` matrix containing data `y_1, ... , y_T`
- `s_0`: `Ns` x 1 initial state vector
- `P_0`: `Ns` x `Ns` initial state covariance matrix

**Method 1 only:** state-space system matrices `T`, `R`, `C`, `Q`, `Z`,
`D`, `E`. See `?kalman_filter`
pp
**Method 2 only:** `regime_indices` and system matrices for each regime `Ts`,
`Rs`, `Cs`, `Qs`, `Zs`, `Ds`, `Es`. See `?kalman_filter`

where:

- `Nt`: number of periods for which we have data
- `Ns`: number of states
- `Ne`: number of shocks
- `Ny`: number of observables

### Keyword Arguments

- `Nt0`: if greater than 0, the returned smoothed states and shocks matrices
  will be shorter by that number of columns (taken from the beginning)
- `draw_states`: indicates whether to draw smoothed states from the distribution
  `N(s_{t|T}, P_{t|T})` or to use the mean `s_{t|T}` (reducing to the Hamilton
  smoother)

### Outputs

- `s_smth`: `Ns` x `Nt` matrix of smoothed states `s_{t|T}`
- `ϵ_smth`: `Ne` x `Nt` matrix of smoothed shocks `ϵ_{t|T}`
"""
function carter_kohn_smoother(y::AbstractMatrix,
    T::Matrix{S}, R::Matrix{S}, C::Vector{S},
    Q::Matrix{S}, Z::Matrix{S}, D::Vector{S}, E::Matrix{S},
    s_0::Vector{S}, P_0::Matrix{S};
    Nt0::Int = 0, draw_states::Bool = true) where {S<:AbstractFloat}

    Nt = size(y, 2)
    carter_kohn_smoother(UnitRange{Int}[1:Nt], y, Matrix{S}[T], Matrix{S}[R], Vector{S}[C],
        Matrix{S}[Q], Matrix{S}[Z], Vector{S}[D], Matrix{S}[E], s_0, P_0;
        Nt0 = Nt0, draw_states = draw_states)
end

function carter_kohn_smoother(regime_indices::Vector{UnitRange{Int}}, y::AbstractMatrix,
    Ts::Vector{Matrix{S}}, Rs::Vector{Matrix{S}}, Cs::Vector{Vector{S}}, Qs::Vector{Matrix{S}},
    Zs::Vector{Matrix{S}}, Ds::Vector{Vector{S}}, Es::Vector{Matrix{S}},
    s_0::Vector{S}, P_0::Matrix{S};
    Nt0::Int = 0, draw_states::Bool = true) where {S<:AbstractFloat}

    # Dimensions
    Nt = size(y,     2) # number of periods of data
    Ns = size(Ts[1], 1) # number of states
    Ne = size(Rs[1], 2) # number of shocks

    # Augment state space with shocks
    Ts, Rs, Cs, Zs, s_0, P_0 =
        augment_states_with_shocks(regime_indices, Ts, Rs, Cs, Qs, Zs, s_0, P_0)

    # Kalman filter stacked states and shocks stil_t
    _, stil_pred, Ptil_pred, stil_filt, Ptil_filt, _, _, _, _ =
        kalman_filter(regime_indices, y, Ts, Rs, Cs, Qs, Zs, Ds, Es, s_0, P_0,
                      outputs = [:pred, :filt])

    # Smooth the stacked states and shocks recursively, starting at t = T-1 and
    # going backwards
    stil_smth = copy(stil_filt)

    for i = length(regime_indices):-1:1
        # Get state-space system matrices for this regime
        T = Ts[i]

        for t in reverse(regime_indices[i])
            if t == Nt
                μ = stil_filt[:, end]
                Σ = Ptil_filt[:, :, end]
            else
                J = Ptil_filt[:, :, t] * T' * pinv(Ptil_pred[:, :, t+1])
                μ = stil_filt[:, t] + J*(stil_smth[:, t+1] - stil_pred[:, t+1]) # stil_{t|T}
                Σ = Ptil_filt[:, :, t] - J*T*Ptil_filt[:, :, t]                 # Ptil_{t|T}
            end

            # Draw stil_t ∼ N(stil_{t|T}, Ptil_{t|T})
            stil_smth[:, t] = if draw_states
                U, eig, _ = svd(Σ)
                μ + U * diagm(0 => (sqrt.(eig))) * randn(Ns+Ne)
            else
                μ
            end
        end
    end

    # Index out states and shocks
    s_smth = stil_smth[1:Ns,     :] # s_{t|T}
    ϵ_smth = stil_smth[Ns+1:end, :] # ϵ_{t|T}

    # Trim the presample if needed
    if Nt0 > 0
        insample = Nt0+1:Nt
        s_smth = s_smth[:, insample]
        ϵ_smth = ϵ_smth[:, insample]
    end

    return s_smth, ϵ_smth
end