text:
"""
    compute_calcium_scores(image, spacing, mask, min_vol=nothing, max_vol=nothing)

Calculates Agatston score and volume given an input mask

Ported from
(https://github.com/sgmvanvelzen/calcium-scoring/blob/f63629f322ec2a22b3e7cb4a7eb6e494eba09fc5/src/calciumscoring/scores.py#L9)
"""
function compute_calcium_scores(image, spacing, mask, min_vol=nothing, max_vol=nothing)
	binary_mask = mask
	voxel_volume = spacing[1] * spacing[2] * spacing[3]

	agatston_score = 0
	calcium_volume = 0

	# Find individual lesions (in 3D) so that we can discard too small or too large lesions
	lesion_map = Images.label_components(binary_mask)
	n_lesions = length(unique(lesion_map))

	for lesion = 1:(n_lesions - 1)
		lesion_mask = map(x -> x == lesion, lesion_map)

		# Ignore too small or too large lesions
		lesion_volume = count(lesion_mask) * voxel_volume
		if ((min_vol !== nothing) && (lesion_volume < min_vol))
			continue
		end
		if ((max_vol !== nothing) && (lesion_volume > max_vol))
			continue
		end

		calcium_volume += lesion_volume

		# Calculate Agatston score for this lesion
		slices = sort(unique(lesion_mask)) .+ 1
		for z = 1:slices[end]
			fragment_mask = lesion_mask[z,:,:]
			n_pixels = count(fragment_mask)
			maximum_intensity = maximum(image[z,:,:][fragment_mask])
            if maximum_intensity < 200
                coefficient = 1
			elseif maximum_intensity < 300
				coefficient = 2
			elseif maximum_intensity < 400
				coefficient = 3
			else
				coefficient = 4
			end
			agatston_score += coefficient * n_pixels
		end
	end
	agatston_score *= spacing[1] / 3.0 * spacing[2] * spacing[3]
	return agatston_score, calcium_volume
end


text:
module MATLAB
    using Compat, Compat.Libdl
    using Compat.Libdl: dlopen, dlsym

    # mxarray
    export MxArray, mxClassID, mxComplexity
    export mxclassid, data_ptr
    export classid, nrows, ncols, nelems, elsize

    export is_double, is_single
    export is_int8, is_uint8, is_int16, is_uint16
    export is_int32, is_uint32, is_int64, is_uint64
    export is_numeric, is_complex, is_sparse, is_empty
    export is_logical, is_char, is_struct, is_cell

    export mxarray, mxempty, mxsparse, delete, duplicate
    export mxcellarray, get_cell, set_cell
    export mxstruct, mxstructarray, nfields, get_fieldname, get_field, set_field
    export jvariable, jarray, jscalar, jvector, jmatrix, jsparse, jstring, jdict

    # mstatments
    export mstatement

    # engine & matfile
    export MSession, MatFile
    export get_default_msession, restart_default_msession, close_default_msession
    export eval_string, get_mvariable, get_variable, put_variable, put_variables
    export variable_names, read_matfile, write_matfile
    export mxcall
    export @mput, @mget, @matlab, @mat_str, @mat_mstr

    import Base.eltype, Base.close, Base.size, Base.copy, Base.ndims, Compat.unsafe_convert, Compat.ASCIIString, Compat.String

    include("exceptions.jl")
    include("mxbase.jl")
    include("mxarray.jl")
    include("matfile.jl")

    include("mstatements.jl")
    include("engine.jl")
    include("matstr.jl")
end # module

ext: jl


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ f409dbe0-17b3-11ec-28ce-e78732e98c73
begin
	import Pkg
	Pkg.activate(joinpath(@__DIR__))
	using SimLib, Plots, PlutoUI, Statistics, LinearAlgebra, LaTeXStrings, XXZNumerics, SpinSymmetry
	gr()
	LinearAlgebra.BLAS.set_num_threads(8)
	cd(joinpath(@__DIR__, ".."))
end

# ╔═╡ 55288508-762b-4087-bf3d-69f5ca43806c
RandomPositionsXXZWithXField(PositionDataDescriptor(:box_pbc, 1,8), PowerLaw(6), nothing, nothing, nothing)

# ╔═╡ dd6040ec-a6b2-425a-b75a-80bd9e4396e4
model = RandomPositionsXXZWithXField(PositionDataDescriptor(:box_pbc, 1,8), PowerLaw(6))

# ╔═╡ fbf8e37f-d687-4a20-8e12-1ef932d35483
eldata = load_el("xmag", model);

# ╔═╡ 7b111ece-3835-442c-aa72-940c3f4eef03
eldata.geometry

# ╔═╡ 08d7e059-dff7-460f-9926-ce1aa38f7915
md"""
#values: $(length(eldata.data))

#Inf   : $(sum(eldata.data .== -Inf))
"""

# ╔═╡ ae169b70-885b-4659-968c-4426f0535a54
size(eldata.data)

# ╔═╡ 989f3a23-8e50-4acc-8cda-1bdb8e79da64
means = meandrop(eldata.data; dims=(1,2,3));

# ╔═╡ 3daa32fc-78b3-43d5-bde6-0cfa24f444ed
md"""
Field: $(@bind fieldIndex Slider(1:length(eldata.fields)))

Density: $(@bind rhoIndex Slider(1:length(eldata.ρs)))
"""

# ╔═╡ 89ffd707-e5ef-485f-ba47-db07e3caa6f2
let p = plot(;legend=nothing,
		title="ρ=$(eldata.ρs[rhoIndex]), h=$(eldata.fields[fieldIndex])")
	histogram!(p, vec(eldata.data[30:end-30,1,:,fieldIndex,rhoIndex]); bins=40)
end

# ╔═╡ eacbb47c-b86f-448c-b072-95e775f94e22
plot(eldata.fields, means; labels="ρ=" .* string.(eldata.ρs'))

# ╔═╡ 3eaa32b2-b41f-4609-8deb-9d69d80796f1
plot(eldata.ρs, means'; labels="h=" .* string.(eldata.fields'))

# ╔═╡ Cell order:
# ╠═f409dbe0-17b3-11ec-28ce-e78732e98c73
# ╠═55288508-762b-4087-bf3d-69f5ca43806c
# ╠═dd6040ec-a6b2-425a-b75a-80bd9e4396e4
# ╠═7b111ece-3835-442c-aa72-940c3f4eef03
# ╠═fbf8e37f-d687-4a20-8e12-1ef932d35483
# ╠═08d7e059-dff7-460f-9926-ce1aa38f7915
# ╠═ae169b70-885b-4659-968c-4426f0535a54
# ╠═989f3a23-8e50-4acc-8cda-1bdb8e79da64
# ╟─3daa32fc-78b3-43d5-bde6-0cfa24f444ed
# ╠═89ffd707-e5ef-485f-ba47-db07e3caa6f2
# ╠═eacbb47c-b86f-448c-b072-95e775f94e22
# ╠═3eaa32b2-b41f-4609-8deb-9d69d80796f1


text:
using MathOptInterfaceSCS
using Base.Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities

using MathOptInterfaceBridges
const MOIB = MathOptInterfaceBridges

MOIU.@model SCSModelData () (EqualTo, GreaterThan, LessThan) (Zeros, Nonnegatives, Nonpositives, SecondOrderCone, ExponentialCone, PositiveSemidefiniteConeTriangle) () (SingleVariable,) (ScalarAffineFunction,) (VectorOfVariables,) (VectorAffineFunction,)

MOIB.@bridge SplitInterval MOIB.SplitIntervalBridge () (Interval,) () () () (ScalarAffineFunction,) () ()
MOIB.@bridge LogDet MOIB.LogDetBridge () () (LogDetConeTriangle,) () () () (VectorOfVariables,) (VectorAffineFunction,)
MOIB.@bridge RootDet MOIB.RootDetBridge () () (RootDetConeTriangle,) () () () (VectorOfVariables,) (VectorAffineFunction,)

# linear9test needs 1e-3 with SCS < 2.0 and 5e-1 with SCS 2.0
# linear2test needs 1e-4
const config = MOIT.TestConfig(atol=1e-4, rtol=1e-4)

@testset "Continuous linear problems" begin
    # AlmostSuccess for linear9 with SCS 2
    MOIT.contlineartest(SplitInterval{Float64}(MOIU.CachingOptimizer(SCSModelData{Float64}(), SCSOptimizer())), config, ["linear9"])
end

@testset "Continuous conic problems" begin
    MOIT.contconictest(RootDet{Float64}(LogDet{Float64}(MOIU.CachingOptimizer(SCSModelData{Float64}(), SCSOptimizer()))), config, ["rsoc", "geomean", "psds", "rootdet", "logdets"])
end


text:
# "smart" summary and REPL printing

function Base.summary(sys::System)
    return "System (base power $(sys.basepower))"
end

function Base.show(io::IO, sys::System)
    println(io, "basepower=$(sys.basepower)")
    show(io, sys.data)
end

function Base.show(io::IO, ::MIME"text/plain", sys::System)
    println(io, "System")
    println(io, "======")
    println(io, "Base Power: $(sys.basepower)\n")
    show(io, MIME"text/plain"(), sys.data)
end

function Base.show(io::IO, ::MIME"text/html", sys::System)
    println(io, "
System
")
    println(io, "
Base Power: $(sys.basepower)
")
    show(io, MIME"text/html"(), sys.data)
end

function Base.summary(tech::TechnicalParams)
    return "$(typeof(tech))"
end

function Base.summary(arc::Arc)
    return "$(get_name(get_from(arc))) -> $(get_name(get_to(arc))): ($(typeof(arc)))"
end

function Base.show(io::IO, ::MIME"text/plain", data::PowerSystemTableData)
    println(io, "$(typeof(data)):")
    println(io, "  directory:  $(data.directory)")
    if !isnothing(data.timeseries_metadata_file)
        println(io, "  timeseries_metadata_file:  $(data.timeseries_metadata_file)")
    end
    println(io, "  basepower:  $(data.basepower)")
    for field in ("branch", "bus", "dcline", "gen", "load", "services")
        print(io, "  $field:  ")
        val = getfield(data, Symbol(field))
        if isnothing(val)
            println(io, "no data")
        else
            println(io, "$(summary(val))")
        end
    end
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#Octahemioctahedron prototype DEBUG
"""
Return the Lar.Struct of octahemioctahedron prototype
"""
function octahemioctahedron()
    a=0.5
    b=0.735
    V, (VV,EV,FV,CV) = Lar.simplex(3, true)
    tetra = [0 1 a a;0 0 a*-1 a; 0 0 b b], EV,FV,CV
    twotetra = Lar.Struct([ tetra ,Lar.r(π,0,0),tetra])
    fourtetra=Lar.Struct([ twotetra ,Lar.r(0,0,π), twotetra])
    Lar.Struct([ fourtetra ,Lar.r(0,0,π/2), fourtetra])
    V,EV,FV,CV = Lar.struct2lar(Lar.Struct([ fourtetra ,Lar.r(0,0,π/2), fourtetra]))
    return Lar.Struct([(V,collect(Set(EV)),FV,CV)])
end


text:
include("sand_03.jl")
using .Sand_03
using Plots
fig = Plots.plot()

NN = 100 # no of cells
tmax = 4.0 # max time
dt = 0.4 # how often to plot
tp = 0.0#5000 # when to start plotting

r= Sand_03.ring(NN; U=0.5)
for t in 0.0:dt:tmax
#for i in 1:10
    t = r.t
    if t ≥ tp
        x0_ = 0.5 * (r.X_[1:end-1] + r.X_[2:end])
        Plots.plot!(x0_[1:end],r.Σ_[1:end],label="t: $t",legend=false)#:bottomright)#,title="Σ")#;legend=false)
        display(fig)
#        Plots.plot!(x0_[1:end],r.J_[1:end],label=false)#,label="t: $t",legend=:topleft)#,title="J")#;legend=false)
#        display(fig)
    end
    Sand_03.tstep!(r,dt)
    #Sand_03.step!(r)
end

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
"""
    Base.show(io::IO, ::MIME"text/plain", m::Serialized{M}) where {M}

Show the [`Serialized`](@ref) model without unnecessarily loading it.
"""
function Base.show(io::IO, ::MIME"text/plain", m::Serialized{M}) where {M}
    print(
        io,
        "Serialized{$M} saved in \"$(m.filename)\" ($(isnothing(m.m) ? "not loaded" : "loaded"))",
    )
end


text:
using Lowranksvd
using Test

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


text:
module FirstOrderLogic

export PropositionalLogic

using Reexport

include("unicode.jl")
include("propositional.jl")

end # module

max_stars_repo_name: TotalVerb/FirstOrderLogic

max_issues_repo_name: TotalVerb/FirstOrderLogic

max_forks_repo_name: TotalVerb/FirstOrderLogic


text:
struct KernelRegression{R <: Union{Vector, Matrix}, T <: Real}
    x::R
    y::Vector{T}
    k::Function
    h::Float64
end

# Want a version for x as a vector (this) and X as a matrix (doesn't exist yet)
function kernelregression(x::Vector{R}, y::Vector{T}; kernel::Symbol = :epanechnikov, bandwidth::Real = NaN, getbandwidth::Function = kNN.bandwidth) where {R <: Any, T <: Real}
    if isnan(bandwidth)
        h_x = getbandwidth(x)
        h_y = getbandwidth(y)
        h = sqrt(h_x * h_y)
    else
        h = bandwidth
    end
    k = kernels[kernel]
    return KernelRegression(x, y, k, h)
end

function StatsBase.predict(model::KernelRegression, x::Real)
    y, n = 0.0, length(model.x)
    h, k = model.h, model.k
    normalizer = 0.0
    for i in 1:n
        # TODO: Get definitions and usage of bandwidth right
        # w_i = (1 / h) * k((x - model.x[i]) / h)
        w_i = h * k.(h * (x - model.x[i]))
        y += model.y[i] * w_i
        normalizer += w_i
    end
    return y / normalizer
end

function StatsBase.predict!(ys::Vector, model::KernelRegression, xs::AbstractVector{T}) where T <: Real
    n = length(xs)
    # @assert length(ys) == n
    for i in 1:n
        ys[i] = predict(model, xs[i])
    end
    return
end

function StatsBase.predict(model::KernelRegression, xs::AbstractVector{T}) where T <: Real
    ys = Array{T}(undef, length(xs))
    predict!(ys, model, xs)
    return ys
end


text:
using DSGE, ModelConstructors, Dates, Test, LinearAlgebra, FileIO, Random, JLD2

writing_output = false # Write output for tests which use random values
if VERSION < v"1.5"
    ver = "111"
else
    ver = "150"
end

m = AnSchorfheide()
system = compute_system(m)
Φ, Ψ, F_ϵ, F_u = DSGE.compute_system_function(system)
zero_sys = DSGE.zero_system_constants(system)

@testset "Check System, Transition, Measurement, and PseduoMeasurement access functions" begin
    @test typeof(system) == System{Float64}
    @test typeof(system[:transition]) == Transition{Float64}
    @test typeof(system[:measurement]) == Measurement{Float64}
    @test typeof(system[:pseudo_measurement]) == PseudoMeasurement{Float64}
    @test typeof(system[:transition][:TTT]) == Matrix{Float64}
    @test typeof(system[:transition][:RRR]) == Matrix{Float64}
    @test typeof(system[:transition][:CCC]) == Vector{Float64}
    @test typeof(system[:measurement][:ZZ]) == Matrix{Float64}
    @test typeof(system[:measurement][:DD]) == Vector{Float64}
    @test typeof(system[:measurement][:QQ]) == Matrix{Float64}
    @test typeof(system[:measurement][:EE]) == Matrix{Float64}
    @test typeof(system[:pseudo_measurement][:ZZ_pseudo]) == Matrix{Float64}
    @test typeof(system[:pseudo_measurement][:DD_pseudo]) == Vector{Float64}
end

@testset "Check miscellaneous functions acting on System types" begin
    @test sum(zero_sys[:CCC]) == 0.
    @test sum(zero_sys[:DD]) == 0.
    @test sum(zero_sys[:DD_pseudo]) == 0.
    @test Φ(ones(size(system[:transition][:TTT],2)), zeros(size(system[:transition][:RRR],2))) ==
        system[:transition][:TTT] * ones(size(system[:transition][:TTT],2))
    @test Ψ(ones(size(system[:transition][:TTT],2))) ==
        system[:measurement][:ZZ] * ones(size(system[:transition][:TTT],2)) + system[:measurement][:DD]
    @test F_ϵ.μ == zeros(3)
    @test F_ϵ.Σ.mat == system[:measurement][:QQ]
    @test F_u.μ == zeros(3)
    @test F_u.Σ.mat == system[:measurement][:EE]
end

@testset "Using compute_system to update an existing system" begin

    # Check updating
    system1 = compute_system(m, system)
    system2 = compute_system(m, system; observables = collect(keys(m.observables))[1:end - 1],
                                  pseudo_observables = collect(keys(m.pseudo_observables))[1:end - 1],
                                  shocks = collect(keys(m.exogenous_shocks))[1:end - 1],
                                  states = collect(keys(m.endogenous_states))[1:end - 1])
    system3 = compute_system(m, system; zero_DD = true)
    system4 = compute_system(m, system; zero_DD_pseudo = true)

    @test system1[:TTT] ≈ system[:TTT]
    @test system1[:RRR] ≈ system[:RRR]
    @test system1[:CCC] ≈ system[:CCC]
    @test system1[:ZZ] ≈ system[:ZZ]
    @test system1[:DD] ≈ system[:DD]
    @test system1[:QQ] ≈ system[:QQ]
    @test system1[:EE] ≈ zeros(size(system[:ZZ], 1), size(system[:ZZ], 1))
    @test system1[:ZZ_pseudo] ≈ system[:ZZ_pseudo]
    @test system1[:DD_pseudo] ≈ system[:DD_pseudo]

    @test system2[:TTT] ≈ system[:TTT][1:end - 1, 1:end - 1]
    @test system2[:RRR] ≈ system[:RRR][1:end - 1, 1:end - 1]
    @test sum(abs.(system2[:CCC])) ≈ sum(abs.(system[:CCC][1:end - 1])) ≈ 0.
    @test system2[:ZZ] ≈ system[:ZZ][1:end - 1, 1:end - 1]
    @test system2[:DD] ≈ system[:DD][1:end - 1]
    @test system2[:QQ] ≈ system[:QQ][1:end - 1, 1:end - 1]
    @test system2[:EE] ≈ zeros(size(system[:ZZ], 1) - 1, size(system[:ZZ], 1) - 1)
    @test system2[:ZZ_pseudo] ≈ system[:ZZ_pseudo][1:end - 1, 1:end - 1]
    @test system2[:DD_pseudo] ≈ system[:DD_pseudo][1:end - 1]

    @test system3[:TTT] ≈ system[:TTT]
    @test system3[:RRR] ≈ system[:RRR]
    @test system3[:CCC] ≈ system[:CCC]
    @test system3[:ZZ] ≈ system[:ZZ]
    @test system3[:DD] ≈ zeros(length(system[:DD]))
    @test system3[:QQ] ≈ system[:QQ]
    @test system3[:EE] ≈ zeros(size(system[:ZZ], 1), size(system[:ZZ], 1))
    @test system3[:ZZ_pseudo] ≈ system[:ZZ_pseudo]
    @test system3[:DD_pseudo] ≈ system[:DD_pseudo]

    @test system4[:TTT] ≈ system[:TTT]
    @test system4[:RRR] ≈ system[:RRR]
    @test system4[:CCC] ≈ system[:CCC]
    @test system4[:ZZ] ≈ system[:ZZ]
    @test system4[:DD] ≈ system[:DD]
    @test system4[:QQ] ≈ system[:QQ]
    @test system4[:EE] ≈ zeros(size(system[:ZZ], 1), size(system[:ZZ], 1))
    @test system4[:ZZ_pseudo] ≈ system[:ZZ_pseudo]
    @test system4[:DD_pseudo] ≈ zeros(length(system[:DD_pseudo]))

    # Check errors
    @test_throws ErrorException compute_system(m, system; observables = [:blah])
    @test_throws ErrorException compute_system(m, system; states = [:blah])
    @test_throws KeyError compute_system(m, system; shocks = [:blah])
    @test_throws ErrorException compute_system(m, system; pseudo_observables = [:blah])
end

@testset "VAR approximation of state space" begin
    m = Model1002("ss10"; custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true),
                                       :add_NominalWageGrowth =>
                                       Setting(:add_NominalWageGrowth, true)))
    system = compute_system(m)
    system = compute_system(m, system; observables = [:obs_hours, :obs_gdpdeflator,
                                                      :laborshare_t, :NominalWageGrowth],
                            shocks = collect(keys(m.exogenous_shocks)))
    yyyyd, xxyyd, xxxxd = DSGE.var_approx_state_space(system[:TTT], system[:RRR], system[:QQ],
                                                 system[:DD], system[:ZZ], system[:EE],
                                                 zeros(size(system[:ZZ], 1),
                                                       DSGE.n_shocks_exogenous(m)),
                                                 4; get_population_moments = true)
    yyyydc, xxyydc, xxxxdc = DSGE.var_approx_state_space(system[:TTT], system[:RRR], system[:QQ],
                                                         system[:DD], system[:ZZ], system[:EE],
                                                         zeros(size(system[:ZZ], 1),
                                                               DSGE.n_shocks_exogenous(m)),
                                                         4; get_population_moments = true,
                                                         use_intercept = true)
    β, Σ = DSGE.var_approx_state_space(system[:TTT], system[:RRR], system[:QQ], system[:DD],
                                  system[:ZZ], system[:EE],
                                  zeros(size(system[:ZZ], 1), DSGE.n_shocks_exogenous(m)),
                                  4; get_population_moments = false)
    βc, Σc = DSGE.var_approx_state_space(system[:TTT], system[:RRR], system[:QQ], system[:DD],
                                  system[:ZZ], system[:EE],
                                  zeros(size(system[:ZZ], 1), DSGE.n_shocks_exogenous(m)),
                                  4; get_population_moments = false, use_intercept = true)

    expmat = load(joinpath(dirname(@__FILE__), "../reference/exp_var_approx_state_space.jld2"))
    @test @test_matrix_approx_eq yyyyd expmat["yyyyd"]
    @test @test_matrix_approx_eq xxyyd expmat["xxyyd"][2:end, :]
    @test @test_matrix_approx_eq xxxxd expmat["xxxxd"][2:end, 2:end]
    @test @test_matrix_approx_eq yyyydc expmat["yyyyd"]
    @test @test_matrix_approx_eq xxyydc expmat["xxyyd"]
    @test @test_matrix_approx_eq xxxxdc expmat["xxxxd"]

    expβ = \(expmat["xxxxd"][2:end, 2:end], expmat["xxyyd"][2:end, :])
    expΣ = expmat["yyyyd"] - expmat["xxyyd"][2:end, :]' * expβ
    expΣ += expΣ'
    expΣ ./= 2
    @test @test_matrix_approx_eq β expβ
    @test @test_matrix_approx_eq Σ expΣ

    expβc = \(expmat["xxxxd"], expmat["xxyyd"])
    expΣc = expmat["yyyyd"] - expmat["xxyyd"]' * expβc
    expΣc += expΣc'
    expΣc ./= 2
    @test @test_matrix_approx_eq βc expβc
    @test @test_matrix_approx_eq Σc expΣc

    # Check DSGEVAR automates this properly
    dsgevar = DSGEVAR(m)
    DSGE.update!(dsgevar, shocks = collect(keys(m.exogenous_shocks)),
                 observables = [:obs_hours, :obs_gdpdeflator, :laborshare_t, :NominalWageGrowth],
                 lags = 4, λ = Inf)
    yyyyd, xxyyd, xxxxd = compute_system(dsgevar; get_population_moments = true)
    yyyydc, xxyydc, xxxxdc = compute_system(dsgevar; get_population_moments = true, use_intercept = true)
    β, Σ = compute_system(dsgevar)
    βc, Σc = compute_system(dsgevar; use_intercept = true)

    @test @test_matrix_approx_eq yyyyd expmat["yyyyd"]
    @test @test_matrix_approx_eq xxyyd expmat["xxyyd"][2:end, :]
    @test @test_matrix_approx_eq xxxxd expmat["xxxxd"][2:end, 2:end]
    @test @test_matrix_approx_eq yyyydc expmat["yyyyd"]
    @test @test_matrix_approx_eq xxyydc expmat["xxyyd"]
    @test @test_matrix_approx_eq xxxxdc expmat["xxxxd"]

    expβ = \(expmat["xxxxd"][2:end, 2:end], expmat["xxyyd"][2:end, :])
    expΣ = expmat["yyyyd"] - expmat["xxyyd"][2:end, :]' * expβ
    expΣ += expΣ'
    expΣ ./= 2
    @test @test_matrix_approx_eq β expβ
    @test @test_matrix_approx_eq Σ expΣ

    expβc = \(expmat["xxxxd"], expmat["xxyyd"])
    expΣc = expmat["yyyyd"] - expmat["xxyyd"]' * expβc
    expΣc += expΣc'
    expΣc ./= 2
    @test @test_matrix_approx_eq βc expβc
    @test @test_matrix_approx_eq Σc expΣc
end

@testset "VAR using DSGE as a prior" begin
    m = Model1002("ss10"; custom_settings =
                  Dict{Symbol,Setting}(:add_laborshare_measurement =>
                                       Setting(:add_laborshare_measurement, true),
                                       :add_NominalWageGrowth =>
                                       Setting(:add_NominalWageGrowth, true)))
    dsgevar = DSGEVAR(m)
    jlddata = load(joinpath(dirname(@__FILE__), "../reference/test_dsgevar_lambda_irfs.jld2"))
    DSGE.update!(dsgevar, shocks = collect(keys(m.exogenous_shocks)),
                 observables = [:obs_hours, :obs_gdpdeflator, :laborshare_t, :NominalWageGrowth],
                 lags = 4, λ = Inf)
    data = jlddata["data"]
    yyyydc1, xxyydc1, xxxxdc1 = compute_system(dsgevar; get_population_moments = true, use_intercept = true)
    βc1, Σc1 = compute_system(dsgevar; use_intercept = true)

    yyyyd2, xxyyd2, xxxxd2 = compute_system(dsgevar, data; get_population_moments = true)
    β2, Σ2 = compute_system(dsgevar, data)

    # Check when λ = Inf
    @test @test_matrix_approx_eq yyyydc1 yyyyd2
    @test @test_matrix_approx_eq xxyydc1 xxyyd2
    @test @test_matrix_approx_eq xxxxdc1 xxxxd2
    @test @test_matrix_approx_eq βc1 β2
    @test @test_matrix_approx_eq Σc1 Σ2

    # Check when λ is finite
    DSGE.update!(dsgevar, λ = 1.)
    Random.seed!(1793) # need to seed for this
    yyyyd, xxyyd, xxxxd = compute_system(dsgevar, data; get_population_moments = true)
    β, Σ = compute_system(dsgevar, data)

    if writing_output
        JLD2.jldopen(joinpath(dirname(@__FILE__), "../reference/test_dsgevar_lambda_irfs_statespace_output_version=" * ver * ".jld2"),
            true, true, true, IOStream) do file
            write(file, "exp_data_beta", β)
            write(file, "exp_data_sigma", Σ)
        end
    end

    file = jldopen(joinpath(dirname(@__FILE__), "../reference/test_dsgevar_lambda_irfs_statespace_output_version=" * ver * ".jld2"), "r")
    saved_β = read(file, "exp_data_beta")
    saved_Σ = read(file, "exp_data_sigma")
    close(file)

    @test @test_matrix_approx_eq saved_β β
    @test @test_matrix_approx_eq saved_Σ Σ
end

@testset "VECM approximation of state space" begin
    matdata = load(joinpath(dirname(@__FILE__), "../reference/vecm_approx_state_space.jld2"))
    nobs = Int(matdata["nvar"])
    p = Int(matdata["nlags"])
    coint = Int(matdata["coint"])
    TTT = matdata["TTT"]
    RRR = matdata["RRR"]
    ZZ = matdata["ZZ"]
    DD = vec(matdata["DD"])
    QQ = matdata["QQ"]
    EE = matdata["EE"]
    MM = matdata["MM"]
    yyyyd, xxyyd, xxxxd = DSGE.vecm_approx_state_space(TTT, RRR, QQ,
                                                       DD, ZZ, EE, MM, nobs,
                                                       p, coint; get_population_moments = true,
                                                       test_GA0 = matdata["GA0"])
    yyyydc, xxyydc, xxxxdc = DSGE.vecm_approx_state_space(TTT, RRR, QQ,
                                                          DD, ZZ, EE, MM, nobs,
                                                          p, coint; get_population_moments = true,
                                                          use_intercept = true, test_GA0 = matdata["GA0"])
    β, Σ = DSGE.vecm_approx_state_space(TTT, RRR, QQ, DD,
                                        ZZ, EE, MM, nobs, p, coint;
                                        get_population_moments = false,
                                        test_GA0 = matdata["GA0"])
    βc, Σc = DSGE.vecm_approx_state_space(TTT, RRR, QQ, DD,
                                          ZZ, EE, MM, nobs, p, coint;
                                          get_population_moments = false, use_intercept = true,
                                          test_GA0 = matdata["GA0"])

    no_int_inds = vcat(1:coint, coint + 2:size(matdata["xxyyd"], 1))
    @test @test_matrix_approx_eq yyyyd matdata["yyyyd"]
    @test @test_matrix_approx_eq xxyyd matdata["xxyyd"][no_int_inds, :]
    @test @test_matrix_approx_eq xxxxd matdata["xxxxd"][no_int_inds, no_int_inds]
    @test @test_matrix_approx_eq yyyydc matdata["yyyyd"]
    @test @test_matrix_approx_eq xxyydc matdata["xxyyd"]
    @test @test_matrix_approx_eq xxxxdc matdata["xxxxd"]

    expβ = \(matdata["xxxxd"][no_int_inds, no_int_inds], matdata["xxyyd"][no_int_inds, :])
    expΣ = matdata["yyyyd"] - matdata["xxyyd"][no_int_inds, :]' * expβ
    @test @test_matrix_approx_eq β expβ
    @test @test_matrix_approx_eq Σ expΣ

    expβc = matdata["xxxxd"] \ matdata["xxyyd"]
    expΣc = matdata["yyyyd"] - matdata["xxyyd"]' * expβc
    @test @test_matrix_approx_eq βc expβc
    @test @test_matrix_approx_eq Σc expΣc

    # Check DSGEVECM automates this properly. Need a proper model before being able to do this though.
    # m = Model1002("ss10")
    # dsgevecm = DSGEVECM(m)
    # DSGE.update!(dsgevecm, shocks = collect(keys(m.exogenous_shocks)),
    #              observables = [:obs_hours, :obs_gdpdeflator, :laborshare_t, :NominalWageGrowth],
    #              lags = 4, λ = Inf)
    # yyyyd, xxyyd, xxxxd = compute_system(dsgevecm; get_population_moments = true)
    # yyyydc, xxyydc, xxxxdc = compute_system(dsgevecm; get_population_moments = true, use_intercept = true)
    # β, Σ = compute_system(dsgevecm)
    # βc, Σc = compute_system(dsgevecm; use_intercept = true)

    # @test @test_matrix_approx_eq yyyyd matdata["yyyyd"]
    # @test @test_matrix_approx_eq xxyyd matdata["xxyyd"][2:end, :]
    # @test @test_matrix_approx_eq xxxxd matdata["xxxxd"][2:end, 2:end]
    # @test @test_matrix_approx_eq yyyydc matdata["yyyyd"]
    # @test @test_matrix_approx_eq xxyydc matdata["xxyyd"]
    # @test @test_matrix_approx_eq xxxxdc matdata["xxxxd"]

    # expβ = \(matdata["xxxxd"][2:end, 2:end], matdata["xxyyd"][2:end, :])
    # expΣ = matdata["yyyyd"] - matdata["xxyyd"][2:end, :]' * expβ
    # expΣ += expΣ'
    # expΣ ./= 2
    # @test @test_matrix_approx_eq β expβ
    # @test @test_matrix_approx_eq Σ expΣ

    # expβc = \(matdata["xxxxd"], matdata["xxyyd"])
    # expΣc = matdata["yyyyd"] - matdata["xxyyd"]' * expβc
    # expΣc += expΣc'
    # expΣc ./= 2
    # @test @test_matrix_approx_eq βc expβc
    # @test @test_matrix_approx_eq Σc expΣc
end

@testset "Updating a system for a DSGEVECM" begin
    dsge = AnSchorfheide()
    m = DSGEVECM(dsge)
    sys = compute_system(dsge)
    Dout1 = DSGE.compute_DD_coint_add(m, sys, [:obs_gdp, :obs_cpi])
    @info "The following 2 warnings about an empty vector are expected"
    sys_dsgevecm, Dout2 = compute_system(m, sys; get_DD_coint_add = true,
                                         cointegrating_add = [:obs_gdp, :obs_cpi])
    Dout3 = DSGE.compute_DD_coint_add(m, sys, Vector{Symbol}(undef, 0))
    _, Dout4 = compute_system(m, sys; get_DD_coint_add = true)

    # Check computing DD_coint_add
    @test Dout1 == sys[:DD][1:2]
    @test Dout2 == sys[:DD][1:2]
    @test isempty(Dout3)
    @test isempty(Dout4)

    # Check the system looks right
    @test @test_matrix_approx_eq sys[:TTT] sys_dsgevecm[:TTT]
    @test @test_matrix_approx_eq sys[:RRR] sys_dsgevecm[:RRR]
    @test @test_matrix_approx_eq sys[:CCC] sys_dsgevecm[:CCC]
    @test @test_matrix_approx_eq sys[:ZZ] sys_dsgevecm[:ZZ]
    @test @test_matrix_approx_eq sys[:DD] sys_dsgevecm[:DD]
    @test @test_matrix_approx_eq sys[:QQ] sys_dsgevecm[:QQ]
    @test @test_matrix_approx_eq sys_dsgevecm[:EE] zeros(size(sys_dsgevecm[:EE]))
end

@testset "RegimeSwitchingSystem" begin
    transitions         = Vector{Transition{Float64}}(undef, 3)
    measurements        = Vector{Measurement{Float64}}(undef, 3)
    pseudo_measurements = Vector{PseudoMeasurement{Float64}}(undef, 3)
    for i in 1:3
        transitions[i]         = system.transition
        measurements[i]        = system.measurement
        pseudo_measurements[i] = system.pseudo_measurement
    end

    regswitch_sys1 = RegimeSwitchingSystem(transitions, measurements, pseudo_measurements)
    regswitch_sys2 = RegimeSwitchingSystem(transitions, measurements)
    regswitch_sys3 = RegimeSwitchingSystem([system, system, system])
    regswitch_sys4 = copy(regswitch_sys1)
    regswitch_sys5 = deepcopy(regswitch_sys1)

    # Test access/utility functions first
    @test regswitch_sys1[:regimes]             == 1:3
    @test regswitch_sys1[:transitions]         == transitions
    @test regswitch_sys1[:measurements]        == measurements
    @test regswitch_sys1[:pseudo_measurements] == pseudo_measurements
    @test n_regimes(regswitch_sys1)            == 3
    for i in 1:3
        # Creating a System for a specific regime
        for tmpsys in [System(regswitch_sys1, i), System(regswitch_sys3, i),
                       System(regswitch_sys4, i), System(regswitch_sys5, i)]
            @test isa(tmpsys, System)
            @test @test_matrix_approx_eq tmpsys[:TTT] system[:TTT]
            @test @test_matrix_approx_eq tmpsys[:RRR] system[:RRR]
            @test @test_matrix_approx_eq tmpsys[:CCC] system[:CCC]
            @test @test_matrix_approx_eq tmpsys[:ZZ]  system[:ZZ]
            @test @test_matrix_approx_eq tmpsys[:DD]  system[:DD]
            @test @test_matrix_approx_eq tmpsys[:QQ]  system[:QQ]
            @test @test_matrix_approx_eq tmpsys[:EE]  system[:EE]
            @test @test_matrix_approx_eq tmpsys[:ZZ_pseudo]  system[:ZZ_pseudo]
            @test @test_matrix_approx_eq tmpsys[:DD_pseudo]  system[:DD_pseudo]
        end

        tmpsys2 = regswitch_sys2[i] # this creates a System underneath the hood
        @test isa(tmpsys2, System)
        @test @test_matrix_approx_eq tmpsys2[:TTT] system[:TTT]
        @test @test_matrix_approx_eq tmpsys2[:RRR] system[:RRR]
        @test @test_matrix_approx_eq tmpsys2[:CCC] system[:CCC]
        @test @test_matrix_approx_eq tmpsys2[:ZZ]  system[:ZZ]
        @test @test_matrix_approx_eq tmpsys2[:DD]  system[:DD]
        @test @test_matrix_approx_eq tmpsys2[:QQ]  system[:QQ]
        @test @test_matrix_approx_eq tmpsys2[:EE]  system[:EE]
        @test all(tmpsys2[:ZZ_pseudo] .== 0.)
        @test all(tmpsys2[:DD_pseudo] .== 0.)

        # Accessing specific data types or regime matrices
        @test @test_matrix_approx_eq regswitch_sys1[i, :transition][:TTT] transitions[i][:TTT]
        @test @test_matrix_approx_eq regswitch_sys1[i, :transition][:RRR] transitions[i][:RRR]
        @test @test_matrix_approx_eq regswitch_sys1[i, :transition][:CCC] transitions[i][:CCC]
        @test @test_matrix_approx_eq regswitch_sys1[i, :TTT]              transitions[i][:TTT]
        @test @test_matrix_approx_eq regswitch_sys1[i, :RRR]              transitions[i][:RRR]
        @test @test_matrix_approx_eq regswitch_sys1[i, :CCC]              transitions[i][:CCC]

        @test @test_matrix_approx_eq regswitch_sys1[i, :measurement][:ZZ] measurements[i][:ZZ]
        @test @test_matrix_approx_eq regswitch_sys1[i, :measurement][:DD] measurements[i][:DD]
        @test @test_matrix_approx_eq regswitch_sys1[i, :measurement][:QQ] measurements[i][:QQ]
        @test @test_matrix_approx_eq regswitch_sys1[i, :measurement][:EE] measurements[i][:EE]
        @test @test_matrix_approx_eq regswitch_sys1[i, :ZZ]               measurements[i][:ZZ]
        @test @test_matrix_approx_eq regswitch_sys1[i, :DD]               measurements[i][:DD]
        @test @test_matrix_approx_eq regswitch_sys1[i, :QQ]               measurements[i][:QQ]
        @test @test_matrix_approx_eq regswitch_sys1[i, :EE]               measurements[i][:EE]

        @test @test_matrix_approx_eq regswitch_sys1[i, :pseudo_measurement][:ZZ_pseudo] pseudo_measurements[i][:ZZ_pseudo]
        @test @test_matrix_approx_eq regswitch_sys1[i, :pseudo_measurement][:DD_pseudo] pseudo_measurements[i][:DD_pseudo]
        @test @test_matrix_approx_eq regswitch_sys1[i, :ZZ_pseudo]                      pseudo_measurements[i][:ZZ_pseudo]
        @test @test_matrix_approx_eq regswitch_sys1[i, :DD_pseudo]                      pseudo_measurements[i][:DD_pseudo]

        # Check errors when trying to access values
        @test_throws KeyError regswitch_sys1[1, :a]
        @test_throws KeyError regswitch_sys1[:a]
        @test_throws BoundsError regswitch_sys1[4, :transition]
        @test_throws BoundsError regswitch_sys1[4]
        @test_throws BoundsError System(regswitch_sys1, 4)

        # Check copying and deepcopying
        oldval = copy(regswitch_sys1[i, :TTT][1, 1])
        regswitch_sys1[i, :TTT][1, 1] = oldval + 1.
        @test @test_matrix_approx_eq regswitch_sys1[i, :TTT] regswitch_sys4[i, :TTT]
        @test !(regswitch_sys1[i, :TTT] ≈ regswitch_sys5[i, :TTT])
        regswitch_sys1[i, :TTT][1, 1] = oldval
    end
end

@testset "Implement alternative policy using regime_eqcond_info" begin
    output_vars = [:forecastobs, :histobs, :histpseudo, :forecastpseudo]

    m = Model1002("ss10", custom_settings = Dict{Symbol, Setting}(:add_altpolicy_pgap =>
                                                                  Setting(:add_altpolicy_pgap, true)))
    m <= Setting(:date_forecast_start, Date(2020, 6, 30))

    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                                  2 => Date(2020, 3, 31),
                                                  3 => Date(2020, 6, 30)))
    m = setup_regime_switching_inds!(m)
    m <= Setting(:forecast_horizons, 12)

    fp = dirname(@__FILE__)
    df = load(joinpath(fp, "../reference", "regime_switch_data.jld2"), "regime_switch_df_none")

    m <= Setting(:replace_eqcond, true)
    m <= Setting(:regime_eqcond_info, Dict{Int, DSGE.EqcondEntry}(
                                                                  3 => DSGE.EqcondEntry(DSGE.ngdp(), [1., 0.])))
    m <= Setting(:pgap_type, :ngdp)
    m <= Setting(:pgap_value, 12.0)

    m <= Setting(:gensys2, false) # don't treat as a temporary policy
    sys1 = compute_system(m)
    #fcast_altperm = DSGE.forecast_one_draw(m, :mode, :full, output_vars, map(x -> x.value, m.parameters),
    #                                       df; regime_switching = true, n_regimes = get_setting(m, :n_regimes))

    # Testing ore basic permanent NGDP
    m = Model1002("ss10"; custom_settings = Dict{Symbol, Setting}(:add_altpolicy_pgap =>
                                                                    Setting(:add_altpolicy_pgap, true)))
    m <= Setting(:forecast_horizons, 12)
    m <= Setting(:date_forecast_start, Date(2020, 6, 30))
    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                                2 => Date(2020, 3, 31),
                                                3 => Date(2020, 6, 30)))
    m = setup_regime_switching_inds!(m)

    m <= Setting(:pgap_value, 12.0)
    m <= Setting(:pgap_type, :ngdp)
    m <= Setting(:regime_eqcond_info, Dict{Int, DSGE.EqcondEntry}(i => DSGE.EqcondEntry(DSGE.ngdp(), [1., 0.]) for i in 1:3))
    m <= Setting(:replace_eqcond, true)
    sys2 = compute_system(m)

    for i in 1:3
        @test !(sys1[1, :TTT] ≈ sys2[i, :TTT])
        @test !(sys1[2, :TTT] ≈ sys2[i, :TTT])
        @test sys1[3, :TTT] ≈ sys2[i, :TTT]
    end
    @test sys1[1, :TTT] ≈ sys1[2, :TTT]
end

@testset "Calculating transition matrices/vectors for k-periods ahead expectations and expected sums" begin
    ## Test k-periods ahead expectations

    # Set up
    m = Model1002("ss10")
    sys_constant = compute_system(m)
    m <= Setting(:date_forecast_start, Date(2020, 9, 30))
    m <= Setting(:date_conditional_end, Date(2020, 9, 30))
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 9, 30),
                                                3 => Date(2020, 12, 31), 4 => Date(2021, 3, 31),
                                                5 => Date(2021, 6, 30), 6 => Date(2021, 9, 30), 7 => Date(2021, 12, 31)))
    m <= Setting(:replace_eqcond, true)
    m <= Setting(:gensys2, true)
    m <= Setting(:regime_switching, true)
    regime_eqcond_info = Dict{Int, DSGE.EqcondEntry}()
    for i in 2:6
        regime_eqcond_info[i] = DSGE.EqcondEntry(DSGE.zlb_rule(), [1., 0.])
    end
    regime_eqcond_info[7] = DSGE.EqcondEntry(AltPolicy(:historical, eqcond, solve), [1., 0.])
    m <= Setting(:regime_eqcond_info, regime_eqcond_info)
    m <= Setting(:temporary_altpolicy_names, [:zlb_rule])
    setup_regime_switching_inds!(m)
    m <= Setting(:tvis_regime_eqcond_info, [regime_eqcond_info])
    m <= Setting(:tvis_information_set, [1:1, 2:7, 3:7, 4:7, 5:7, 6:7, 7:7])
    m <= Setting(:tvis_select_system, ones(Int, 7))
    sys = compute_system(m)

    # Test k-periods ahead expectations with empty TTTs
    T_acc1, C_acc1 = DSGE.k_periods_ahead_expectations(sys_constant[:TTT], sys_constant[:CCC], typeof(sys_constant[:TTT])[],
                                                       typeof(sys_constant[:CCC])[], 1, 3)
    T_acc2, C_acc2 = DSGE.k_periods_ahead_expectations(sys_constant[:TTT], sys_constant[:CCC], typeof(sys_constant[:TTT])[],
                                                       typeof(sys_constant[:CCC])[], 1, 3, 3)
    TTT_constant³ = sys_constant[:TTT] ^ 3
    @test T_acc1 ≈ TTT_constant³
    @test T_acc2 ≈ TTT_constant³
    @test C_acc1 ≈ sys_constant[:CCC]
    @test C_acc2 ≈ sys_constant[:CCC]

    T_acc3, C_acc3 = DSGE.k_periods_ahead_expectations(sys_constant[:TTT], sys_constant[:CCC] .+ .1, typeof(sys_constant[:TTT])[],
                                                       typeof(sys_constant[:CCC])[], 1, 3)
    TTT_constant³sum = (I - sys_constant[:TTT]) \ (I - TTT_constant³)
    @test T_acc3 ≈ TTT_constant³
    @test C_acc3 ≈ TTT_constant³sum * (sys_constant[:CCC] .+ .1)
    @test !(all(C_acc3 .≈ 0.))

    # Test k-periods ahead expectations with TTTs and CCCs (but CCCs are possibly zero)
    TTTs = [sys[i, :TTT] for i in 1:7]
    CCCs = [sys[i, :CCC] for i in 1:7]
    T_acc1, C_acc1 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3)
    T_acc2, C_acc2 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3, 7)
    T_acc3, C_acc3 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 5)
    T_acc4, C_acc4 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 7)
    T_acc5, C_acc5 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 7, 3)

    @test T_acc1 ≈ sys[5, :TTT] * sys[4, :TTT] * sys[3, :TTT]
    @test T_acc2 ≈ sys[5, :TTT] * sys[4, :TTT] * sys[3, :TTT]
    @test T_acc3 ≈ sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :TTT] * sys[3, :TTT]
    @test T_acc4 ≈ sys[7, :TTT] * sys[7, :TTT] * sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :TTT] * sys[3, :TTT]
    @test T_acc5 ≈ sys[7, :TTT] ^ 3
    @test C_acc1 ≈ sys[5, :CCC] + sys[5, :TTT] * sys[4, :CCC] + sys[5, :TTT] * sys[4, :TTT] * sys[3, :CCC]
    @test C_acc2 ≈ sys[5, :CCC] + sys[5, :TTT] * sys[4, :CCC] + sys[5, :TTT] * sys[4, :TTT] * sys[3, :CCC]
    @test C_acc3 ≈ sys[7, :CCC] + sys[7, :TTT] * sys[6, :CCC] + sys[7, :TTT] * sys[6, :TTT] * sys[5, :CCC] +
        sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :CCC] + sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :TTT] * sys[3, :CCC]
    @test C_acc4 ≈ sys[7, :TTT]^2 * (sys[7, :CCC] + sys[7, :TTT] * sys[6, :CCC] + sys[7, :TTT] * sys[6, :TTT] * sys[5, :CCC] +
                                     sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :CCC] +
                                     sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :TTT] * sys[3, :CCC])
    @test all(C_acc5 .≈ 0.)

    # Test ForwardExpectationsMemo
    memo = DSGE.ForwardExpectationsMemo(TTTs, 3, 7, 7, 1, 3)
    @test !haskey(memo.time_varying_memo, 7)
    @test memo.time_varying_memo[6] ≈ TTTs[6]
    @test memo.time_varying_memo[5] ≈ TTTs[6] * TTTs[5]
    @test memo.time_varying_memo[4] ≈ TTTs[6] * TTTs[5] * TTTs[4]
    for k in 1:3
        @test memo.permanent_memo[k] ≈ TTTs[1] ^ k
    end
    memo = DSGE.ForwardExpectationsMemo(TTTs, 1, 1, 7)
    @test isempty(memo.permanent_memo)

    # Create memo for k_periods_ahead_expectations
    memo1 = DSGE.ForwardExpectationsMemo(TTTs, 2, 2 + 3, 7)
    memo2 = DSGE.ForwardExpectationsMemo(TTTs, 2, 6, 7, 1, 3)

    T_acc_memo1, C_acc_memo1 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3, memo = memo1)
    T_acc_memo2, C_acc_memo2 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3, 7, memo = memo1)
    T_acc_memo3, C_acc_memo3 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 5, memo = memo2)
    T_acc_memo4, C_acc_memo4 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 7, memo = memo2)
    T_acc_memo5, C_acc_memo5 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 7, 3, memo = memo2)

    @test T_acc_memo1 ≈ T_acc1
    @test C_acc_memo1 ≈ C_acc1
    @test T_acc_memo2 ≈ T_acc2
    @test C_acc_memo2 ≈ C_acc2
    @test T_acc_memo3 ≈ T_acc3
    @test C_acc_memo3 ≈ C_acc3
    @test T_acc_memo4 ≈ T_acc4
    @test C_acc_memo4 ≈ C_acc4
    @test T_acc_memo5 ≈ T_acc5
    @test C_acc_memo5 ≈ C_acc5

    # Test k-periods ahead expectations with all nonzero-CCCs
    CCCs = [sys[i, :CCC] .+ .1 for i in 1:7]
    T_acc1, C_acc1 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3)
    T_acc2, C_acc2 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3, 7)
    T_acc3, C_acc3 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 5)
    T_acc4, C_acc4 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 7)
    T_acc5, C_acc5 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 7, 3)

    @test T_acc1 ≈ sys[5, :TTT] * sys[4, :TTT] * sys[3, :TTT]
    @test T_acc2 ≈ sys[5, :TTT] * sys[4, :TTT] * sys[3, :TTT]
    @test T_acc3 ≈ sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :TTT] * sys[3, :TTT]
    @test T_acc4 ≈ sys[7, :TTT] * sys[7, :TTT] * sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :TTT] * sys[3, :TTT]
    @test T_acc5 ≈ sys[7, :TTT] ^ 3
    @test C_acc1 ≈ CCCs[5] + sys[5, :TTT] * CCCs[4] + sys[5, :TTT] * sys[4, :TTT] * CCCs[3]
    @test C_acc2 ≈ CCCs[5] + sys[5, :TTT] * CCCs[4] + sys[5, :TTT] * sys[4, :TTT] * CCCs[3]
    @test C_acc3 ≈ CCCs[7] + sys[7, :TTT] * CCCs[6] + sys[7, :TTT] * sys[6, :TTT] * CCCs[5] +
        sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * CCCs[4] + sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :TTT] * CCCs[3]
    @test C_acc4 ≈ (I + sys[7, :TTT]) * CCCs[7] + sys[7, :TTT]^2 * (CCCs[7] + sys[7, :TTT] * CCCs[6] + sys[7, :TTT] * sys[6, :TTT] * CCCs[5] +
                                                                    sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * CCCs[4] +
                                                                    sys[7, :TTT] * sys[6, :TTT] * sys[5, :TTT] * sys[4, :TTT] * CCCs[3])
    @test !all(C_acc5 .≈ 0.)

    # Test memo version
    T_acc_memo1, C_acc_memo1 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3; memo = memo1)
    T_acc_memo2, C_acc_memo2 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3, 7; memo = memo1)
    T_acc_memo3, C_acc_memo3 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 5; memo = memo2)
    T_acc_memo4, C_acc_memo4 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 7; memo = memo2)
    T_acc_memo5, C_acc_memo5 = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 7, 3; memo = memo2)

    @test T_acc_memo1 ≈ T_acc1
    @test C_acc_memo1 ≈ C_acc1
    @test T_acc_memo2 ≈ T_acc2
    @test C_acc_memo2 ≈ C_acc2
    @test T_acc_memo3 ≈ T_acc3
    @test C_acc_memo3 ≈ C_acc3
    @test T_acc_memo4 ≈ T_acc4
    @test C_acc_memo4 ≈ C_acc4
    @test T_acc_memo5 ≈ T_acc5
    @test C_acc_memo5 ≈ C_acc5

    ## Test k-periods ahead expected sum

    # Set up
    m = Model1002("ss10")
    sys_constant = compute_system(m)
    m <= Setting(:date_forecast_start, Date(2020, 9, 30))
    m <= Setting(:date_conditional_end, Date(2020, 9, 30))
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 9, 30),
                                                3 => Date(2020, 12, 31), 4 => Date(2021, 3, 31),
                                                5 => Date(2021, 6, 30), 6 => Date(2021, 9, 30), 7 => Date(2021, 12, 31)))
    m <= Setting(:replace_eqcond, true)
    m <= Setting(:gensys2, true)
    m <= Setting(:regime_switching, true)
    regime_eqcond_info = Dict()
    for i in 2:6
        regime_eqcond_info[i] = DSGE.EqcondEntry(DSGE.zlb_rule(), [1., 0.])
    end
    regime_eqcond_info[7] = DSGE.EqcondEntry(AltPolicy(:historical, eqcond, solve), [1., 0.])
    m <= Setting(:regime_eqcond_info, regime_eqcond_info)
    m <= Setting(:temporary_altpolicy_names, [:zlb_rule])
    setup_regime_switching_inds!(m)
    sys = compute_system(m)

    # Start with empty TTTs
    T_acc1, C_acc1 = DSGE.k_periods_ahead_expected_sums(sys_constant[:TTT], sys_constant[:CCC], typeof(sys_constant[:TTT])[],
                                                        typeof(sys_constant[:CCC])[], 1, 3)
    T_acc2, C_acc2 = DSGE.k_periods_ahead_expected_sums(sys_constant[:TTT], sys_constant[:CCC], typeof(sys_constant[:TTT])[],
                                                        typeof(sys_constant[:CCC])[], 1, 3, 3)
    T_accum = Vector{typeof(sys_constant[:TTT])}(undef, 3)
    C_accum = Vector{typeof(sys_constant[:CCC])}(undef, 3)
    for i in 1:3
        T_accum[i], C_accum[i] = DSGE.k_periods_ahead_expectations(sys_constant[:TTT], sys_constant[:CCC], typeof(sys_constant[:TTT])[],
                                                                   typeof(sys_constant[:CCC])[], 1, i)
    end
    T_accum = sum(T_accum)
    C_accum = sum(C_accum)

    @test T_acc1 ≈ T_accum
    @test T_acc2 ≈ T_accum
    @test C_acc1 ≈ C_accum
    @test C_acc2 ≈ C_accum

    T_acc3, C_acc3 = DSGE.k_periods_ahead_expected_sums(sys_constant[:TTT], sys_constant[:CCC] .+ .1, typeof(sys_constant[:TTT])[],
                                                        typeof(sys_constant[:CCC])[], 1, 3)
    C_accum = Vector{typeof(sys_constant[:CCC])}(undef, 3)
    for i in 1:3
        _, C_accum[i] = DSGE.k_periods_ahead_expectations(sys_constant[:TTT], sys_constant[:CCC] .+ .1, typeof(sys_constant[:TTT])[],
                                                                   typeof(sys_constant[:CCC])[], 1, i)
    end
    @test T_acc3 ≈ T_accum
    @test C_acc3 ≈ sum(C_accum)

    # Now do with TTTs and CCCs (some of which are possibly zeros)
    TTTs = [sys[i, :TTT] for i in 1:7]
    CCCs = [sys[i, :CCC] for i in 1:7]
    T_acc1, C_acc1 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3)
    T_acc2, C_acc2 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3, 7)
    T_acc3, C_acc3 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 5)
    T_acc4, C_acc4 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 7)
    T_acc5, C_acc5 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 7, 3)

    # Test ForwardMultipleExpectationsMemo
    memo = DSGE.ForwardMultipleExpectationsMemo(TTTs, 3, 4, 7, 7, 1, 3)
    @test all([!haskey(memo.time_varying_memo, i) for i in 1:3])
    for i in 4:6
        @test !haskey(memo.time_varying_memo[i], 7)
        for j in i:-1:4
            @test memo.time_varying_memo[i][j] ≈ prod([TTTs[l] for l in i:-1:j])
        end
    end
    @test !haskey(memo.time_varying_memo[7], 7)
    for j in 6:-1:4 # from 6 instead of 7 b/c 7 == first permanent period, so we stick the first power of the permanent
        @test memo.time_varying_memo[7][j] ≈ prod([TTTs[l] for l in 6:-1:j]) # transition matrix inside memo.permanent_memo
    end
    for k in 1:3
        @test memo.permanent_memo[k] ≈ TTTs[1] ^ k
    end
    @test !haskey(memo.permanent_memo, 4)
    memo = DSGE.ForwardMultipleExpectationsMemo(TTTs, 1, 1, 1, 7)
    @test isempty(memo.permanent_memo)

    # Now actually use the memo versions
    memo1 = DSGE.ForwardMultipleExpectationsMemo(TTTs, 2, 3, 7, 7)
    memo2 = DSGE.ForwardMultipleExpectationsMemo(TTTs, 2, 3, 7, 7, 1, 4)

    T_acc_memo1, C_acc_memo1 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3; memo = memo1)
    T_acc_memo2, C_acc_memo2 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 3, 7; memo = memo1)
    T_acc_memo3, C_acc_memo3 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 5; memo = memo2)
    T_acc_memo4, C_acc_memo4 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, 7; memo = memo2)
    T_acc_memo5, C_acc_memo5 = DSGE.k_periods_ahead_expected_sums(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 7, 3; memo = memo2)

    # Test output from above against doing it manually, as shown here
    T_accum = Vector{typeof(sys[1, :TTT])}(undef, 3)
    C_accum = Vector{typeof(sys[1, :CCC])}(undef, 3)
    for i in 1:3
        T_accum[i], C_accum[i] = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, i)
    end
    @test T_acc1 ≈ sum(T_accum)
    @test T_acc2 ≈ sum(T_accum)
    @test C_acc1 ≈ sum(C_accum)
    @test C_acc2 ≈ sum(C_accum)
    @test T_acc1 ≈ T_acc_memo1 # breaks
    @test T_acc2 ≈ T_acc_memo2
    @test C_acc1 ≈ C_acc_memo1
    @test C_acc2 ≈ C_acc_memo1

    T_accum = Vector{typeof(sys[1, :TTT])}(undef, 5)
    C_accum = Vector{typeof(sys[1, :CCC])}(undef, 5)
    for i in 1:5
        T_accum[i], C_accum[i] = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, i)
    end
    @test T_acc3 ≈ sum(T_accum)
    @test C_acc3 ≈ sum(C_accum)
    @test T_acc3 ≈ T_acc_memo3
    @test C_acc3 ≈ C_acc_memo3

    T_accum = Vector{typeof(sys[1, :TTT])}(undef, 7)
    C_accum = Vector{typeof(sys[1, :CCC])}(undef, 7)
    for i in 1:7
        T_accum[i], C_accum[i] = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 2, i)
    end
    @test T_acc4 ≈ sum(T_accum)
    @test C_acc4 ≈ sum(C_accum)
    @test T_acc4 ≈ T_acc_memo4
    @test C_acc4 ≈ C_acc_memo4

    T_accum = Vector{typeof(sys[1, :TTT])}(undef, 3)
    C_accum = Vector{typeof(sys[1, :CCC])}(undef, 3)
    for i in 1:3
        T_accum[i], C_accum[i] = DSGE.k_periods_ahead_expectations(sys[1, :TTT], sys[1, :CCC], TTTs, CCCs, 7, i)
    end
    @test T_acc5 ≈ sum(T_accum)
    @test C_acc5 ≈ sum(C_accum)
    @test T_acc5 ≈ T_acc_memo5
    @test C_acc5 ≈ C_acc_memo5
end

@testset "Time-Varying Information Set state space system" begin
    m = Model1002("ss10")
    system = compute_system(m)
    m <= Setting(:date_forecast_start, Date(2020, 9, 30))
    m <= Setting(:date_conditional_end, Date(2020, 9, 30))
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 9, 30),
                                                3 => Date(2020, 12, 31), 4 => Date(2021, 3, 31),
                                                5 => Date(2021, 6, 30), 6 => Date(2021, 9, 30), 7 => Date(2021, 12, 31)))
    setup_regime_switching_inds!(m)
    m <= Setting(:replace_eqcond, true)
    m <= Setting(:gensys2, true)
    m <= Setting(:regime_switching, true)

    # Check that using the same regime_eqcond_info but different tvis_select_system yields the same results
    regime_eqcond_info1 = Dict()
    for i in 2:6
        regime_eqcond_info1[i] = DSGE.EqcondEntry(DSGE.zlb_rule())#, [1., 0.])
    end
    regime_eqcond_info1[7] = DSGE.EqcondEntry(AltPolicy(:historical, eqcond, solve))#, [1., 0.])
    regime_eqcond_info2 = deepcopy(regime_eqcond_info1)
    m <= Setting(:tvis_regime_eqcond_info, [regime_eqcond_info1, regime_eqcond_info2])
    m <= Setting(:tvis_information_set, [1:1, 2:7, 3:7, 4:7, 5:7, 6:7, 7:7])
    m <= Setting(:temporary_altpolicy_names, [:zlb_rule])
    m <= Setting(:tvis_select_system, ones(Int, 7))
    sys1 = DSGE.compute_tvis_system(m)
    m <= Setting(:tvis_select_system, fill(2, 7))
    sys2 = DSGE.compute_tvis_system(m)

    @test isa(sys1[:transitions], Vector{Vector{Transition{Float64}}})
    @test isa(sys1[:measurements], Vector{Measurement{Float64}})
    @test isa(sys1[:pseudo_measurements], Vector{PseudoMeasurement{Float64}})
    @test sys1[:information_set] == sys2.information_set
    @test all(sys1[:select] .!= sys2.select)
    @test sys1[1, :measurement][:ZZ] == sys1.measurements[1][:ZZ]
    @test sys1[1, :pseudo_measurement][:ZZ_pseudo] == sys1.pseudo_measurements[1][:ZZ_pseudo]
    @test sys1[1, :ZZ] == sys1.measurements[1][:ZZ]
    @test sys1[1, :DD] == sys1.measurements[1][:DD]
    @test sys1[1, :ZZ_pseudo] == sys1.pseudo_measurements[1][:ZZ_pseudo]
    @test sys1[1, 2][:TTT] == sys1.transitions[1][2][:TTT]
    @test sys1[1, 2, :TTT] == sys1.transitions[1][2][:TTT]
    @test n_regimes(sys1) == 7
    @test sys1[:regimes] == 1:7
    for reg in 1:n_regimes(sys1)
        sys1[reg, :ZZ] == sys2[reg, :ZZ]
    end

    # Check using the different regime_eqcond_info
    regime_eqcond_info2[6] = DSGE.EqcondEntry(AltPolicy(:historical, eqcond, solve))#, [1., 0.]) # shorten the ZLB
    m <= Setting(:tvis_regime_eqcond_info, [regime_eqcond_info1, regime_eqcond_info2])
    m <= Setting(:tvis_information_set, [1:1, 2:7, 3:7, 4:7, 5:7, 6:7, 7:7])
    m <= Setting(:tvis_select_system, [1, 1, 1, 1, 2, 2, 2])
    sys1 = DSGE.compute_tvis_system(m) # Does not seem to be doing regime-switching
    m <= Setting(:regime_eqcond_info, regime_eqcond_info1)
    m <= Setting(:tvis_regime_eqcond_info, [regime_eqcond_info1])
    m <= Setting(:tvis_select_system, ones(Int, 7))
    sys2 = compute_system(m; tvis = true)
    m <= Setting(:regime_eqcond_info, regime_eqcond_info2)
    m <= Setting(:tvis_regime_eqcond_info, [regime_eqcond_info2])
    m <= Setting(:tvis_select_system, ones(Int, 7)) # ones here b/c there's only one system for the TVIS to use
    sys3 = compute_system(m; tvis = true)
    for reg in 1:4
        @test sys1[reg, :ZZ] ≈ sys2[reg, :ZZ]
    end
    for reg in 5:7
        @test sys1[reg, :ZZ] ≈ sys3[reg, :ZZ]
    end
    for reg in 2:5
        @test sys2[reg, :ZZ] != sys3[reg, :ZZ]
    end
    for reg in 6:7 # same for 6 and 7 b/c in both cases, next period expectation does not use zero rate rule
        @test sys2[reg, :ZZ] == sys3[reg, :ZZ]
    end

    m <= Setting(:tvis_regime_eqcond_info, [regime_eqcond_info1, regime_eqcond_info2])
    m <= Setting(:tvis_select_system, [1, 1, 1, 1, 2, 2, 2])
    sys4 = compute_system(m; tvis = true)
    for reg in 1:4
        @test sys4[reg, :ZZ] == sys1[reg, :ZZ]
        @test sys4[reg, :ZZ] == sys2[reg, :ZZ]
        @test sys4[reg, :TTT] == sys1[1, reg, :TTT]
        @test sys4[reg, :TTT] == sys2[reg, :TTT]
    end
    for reg in 5:7
        @test sys4[reg, :ZZ] == sys1[reg, :ZZ]
        @test sys4[reg, :ZZ] == sys3[reg, :ZZ]
        @test sys4[reg, :TTT] == sys1[2, reg, :TTT]
        @test sys4[reg, :TTT] == sys3[reg, :TTT]
    end
    for reg in 2:4
        @test sys4[reg, :ZZ] != sys3[reg, :ZZ]
    end
    for reg in 5:5 # only to 5 b/c the difference in info sets about ZLB length
        @test sys4[reg, :ZZ] != sys2[reg, :ZZ]
    end
end

@testset "Using gensys2 in historical regimes" begin
    function set_regime_vals_fnct!(m, n)
        if n > 4
            for p in m.parameters
                if !isempty(p.regimes) && haskey(p.regimes, :value)
                    for i in 5:n
                        ModelConstructors.set_regime_val!(p, i, ModelConstructors.regime_val(p, 4))
                    end
                end
            end
        end
    end

    imperfect_cred_new = 1.
    imperfect_cred_old = 1. - imperfect_cred_new
    custom_set = Dict{Symbol,Setting}(:n_mon_anticipated_shocks =>
                                      Setting(:n_mon_anticipated_shocks, 6, "Number of anticipated policy shocks"),
                                      :imperfect_awareness_weights => Setting(:imperfect_awareness_weights,
                                                                              [imperfect_cred_new, imperfect_cred_old]),
                                      :alternative_policies => Setting(:alternative_policies, [DSGE.taylor_rule()]),
                                      :flexible_ait_policy_change => Setting(:flexible_ait_policy_change, true))

    m = Model1002("ss59", custom_settings = custom_set)
    usual_model_settings!(m, "201117", fcast_date = Date(2020, 9, 30))

    set_regime_vals_fnct!(m, 4 + 4)
    m <= Setting(:replace_eqcond, true)
    m <= Setting(:gensys2, true)
    reg_dates = deepcopy(get_setting(m, :regime_dates))
    regime_eqcond_info = Dict{Int, DSGE.EqcondEntry}()
    for (regind, date) in zip(4:(4 + 4), # 4 + 1 b/c zero for 4 periods and liftoff, add 3 to make correct regime
                              DSGE.quarter_range(reg_dates[4], DSGE.iterate_quarters(reg_dates[4], 4)))
        reg_dates[regind] = date
        if regind != 4 + 4
            regime_eqcond_info[regind] = DSGE.EqcondEntry(DSGE.zlb_rule(), [1., 0.])
        end
    end
    regime_eqcond_info[8] = DSGE.EqcondEntry(DSGE.flexible_ait(), [imperfect_cred_new, imperfect_cred_old])
    m <= Setting(:regime_dates, reg_dates)
    m <= Setting(:regime_eqcond_info, regime_eqcond_info)
    setup_regime_switching_inds!(m; cond_type = :full)
    m <= Setting(:tvis_information_set, [1:1, 2:2, 3:3, [i:get_setting(m, :n_regimes) for i in 4:get_setting(m, :n_regimes)]...])
    m <= Setting(:temporary_altpolicy_names, [:zlb_rule])

    sys_true = compute_system(m; tvis = true)
    m <= Setting(:uncertain_temporary_altpolicy, true)
    sys_true_unc = compute_system(m; tvis = true)
    m <= Setting(:uncertain_temporary_altpolicy, false)

    m <= Setting(:date_forecast_start, Date(2020, 12, 31))
    m <= Setting(:date_conditional_end, Date(2020, 12, 31))
    setup_regime_switching_inds!(m; cond_type = :full)
    sys = compute_system(m; tvis = true)
    for i in 1:get_setting(m, :n_regimes)
        @test sys_true[i, :TTT] ≈ sys[i, :TTT]
        @test sys_true[i, :ZZ] ≈ sys[i, :ZZ]
    end

    m <= Setting(:uncertain_temporary_altpolicy, true)
    sys = compute_system(m; tvis = true)
    for i in 1:get_setting(m, :n_regimes)
        @test sys_true_unc[i, :TTT] ≈ sys[i, :TTT]
        @test sys_true_unc[i, :ZZ] ≈ sys[i, :ZZ]
    end
end

@testset "Measurement equation of forward-looking variables" begin
    function set_regime_vals_fnct!(m, n)
        if n > 4
            for p in m.parameters
                if !isempty(p.regimes) && haskey(p.regimes, :value)
                    for i in 5:n
                        ModelConstructors.set_regime_val!(p, i, ModelConstructors.regime_val(p, 4))
                    end
                end
            end
        end
    end
    imperfect_cred_new = .5
    imperfect_cred_old = 1. - imperfect_cred_new
    custom_set = Dict{Symbol,Setting}(:n_mon_anticipated_shocks =>
                                      Setting(:n_mon_anticipated_shocks, 6, "Number of anticipated policy shocks"),
                                      :imperfect_awareness_weights => Setting(:imperfect_awareness_weights,
                                                                              [imperfect_cred_new, imperfect_cred_old]),
                                      :alternative_policies => Setting(:alternative_policies, [DSGE.taylor_rule()]),
                                      :flexible_ait_policy_change => Setting(:flexible_ait_policy_change, true))

    m = Model1002("ss59", custom_settings = custom_set)
    usual_model_settings!(m, "201117", fcast_date = Date(2020, 9, 30))
    m <= Setting(:flexible_ait_policy_change, false) # was true just so initialize a bunch of settings

    set_regime_vals_fnct!(m, 4 + 4)
    m <= Setting(:replace_eqcond, true)
    m <= Setting(:gensys2, true)
    reg_dates = deepcopy(get_setting(m, :regime_dates))
    regime_eqcond_info = Dict{Int, DSGE.EqcondEntry}()
    for (regind, date) in zip(4:(4 + 4), # 4 + 1 b/c zero for 4 periods and liftoff, add 3 to make correct regime
                              DSGE.quarter_range(reg_dates[4], DSGE.iterate_quarters(reg_dates[4], 4)))
        reg_dates[regind] = date
        if regind != 4 + 4
            regime_eqcond_info[regind] = DSGE.EqcondEntry(DSGE.zlb_rule(), [imperfect_cred_new, imperfect_cred_old])
        end
    end
    regime_eqcond_info[8] = DSGE.EqcondEntry(DSGE.flexible_ait(), [imperfect_cred_new, imperfect_cred_old])
    regime_eqcond_info_cp = deepcopy(regime_eqcond_info)
    m <= Setting(:regime_dates, reg_dates)
    m <= Setting(:regime_eqcond_info, regime_eqcond_info)
    setup_regime_switching_inds!(m; cond_type = :full)
    m <= Setting(:temporary_altpolicy_names, [:zlb_rule])

    m <= Setting(:regime_switching, false)
    m <= Setting(:uncertain_temporary_altpolicy, false)
    m <= Setting(:uncertain_altpolicy, false)
    delete!(m.settings, :regime_eqcond_info)
    m <= Setting(:gensys2, false)
    sys_taylor = compute_system(m)

    m <= Setting(:regime_switching, true)
    m <= Setting(:uncertain_temporary_altpolicy, false)
    m <= Setting(:regime_eqcond_info, deepcopy(regime_eqcond_info_cp))
    m <= Setting(:uncertain_altpolicy, false)
    m <= Setting(:gensys2, true)
    sys_perfcred = compute_system(m)

    m <= Setting(:regime_eqcond_info, deepcopy(regime_eqcond_info_cp))
    m <= Setting(:uncertain_temporary_altpolicy, true)
    m <= Setting(:uncertain_altpolicy, true)
    sys_unczlb_uncalt = compute_system(m)

    for i in sort!(collect(keys(get_setting(m, :regime_eqcond_info))))
        @test sys_unczlb_uncalt[i, :ZZ] ≈ imperfect_cred_new * sys_perfcred[i, :ZZ] +
            imperfect_cred_old * sys_taylor[:ZZ]
        @test sys_unczlb_uncalt[i, :ZZ_pseudo] ≈ imperfect_cred_new * sys_perfcred[i, :ZZ_pseudo] +
            imperfect_cred_old * sys_taylor[:ZZ_pseudo]
        @test sys_unczlb_uncalt[i, :DD] ≈ imperfect_cred_new * sys_perfcred[i, :DD] +
            imperfect_cred_old * sys_taylor[:DD]
        @test sys_unczlb_uncalt[i, :DD_pseudo] ≈ imperfect_cred_new * sys_perfcred[i, :DD_pseudo] +
            imperfect_cred_old * sys_taylor[:DD_pseudo]
    end
end

nothing


text:
module SetPropagation_FEM_Examples

# See docs in NeuralNetworkAnalysis.jl
macro modelpath(model_path::String, name::String)
    __source__.file === nothing && return nothing
    _dirname = dirname(String(__source__.file))
    dir = isempty(_dirname) ? pwd() : abspath(_dirname)
    return joinpath(dir, name)
end

end # module

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
module STBImage

include("LibSTBImage.jl")
using .LibSTBImage

end # module


text:
export graph_ps, graph_ps_degopt

"""
     (graph,crefs)=graph_ps(a; input=:A,
                          B_base=:B, C_base=:C, P_base=:P)

Generates the graph for the Paterson–Stockmeyer procedure with monomial basis
coefficieents. More precisely, it corresponds to evaluation of the polynomial

    p(s) = a[1] + a[2]*s + ... + a[n]*s^(n-1).

The code follows the description in the second of the papers referenced below.

**References**

1. M. Paterson and L. Stockmeyer. "On the number of nonscalar multiplications
   necessary to evaluate polynomials". SIAM Journal on Scientific Computing,
   2(1):60-66, 1973. DOI: [10.1137/0202007](https://doi.org/10.1137/0202007)

2. M. Fasi. "Optimality of the Paterson–Stockmeyer method for evaluating matrix
   polynomials and rational matrix functions". Linear Algebra and its
   Applications, 574, 2019.
   DOI: [10.1016/j.laa.2019.04.001](https://doi.org/10.1016/j.laa.2019.04.001)
"""
function graph_ps(a; input = :A, B_base = :B, C_base = :C, P_base = :P)

    # Initial setup
    n = length(a)
    cref = Vector{Tuple{Symbol,Int}}(undef, n)
    graph = Compgraph(eltype(a))

    k = n - 1 # Polynomial degree

    # Degenerate cases k = 0 and k = 1
    if k < 2
        if k == 0
            add_lincomb!(graph, :P0, a[1], :I, 0.0, :I)
            cref[1] = (:P0, 1)
        elseif k == 1
            add_lincomb!(graph, :P0, a[2], :A, a[1], :I)
            cref[1] = (:P0, 1)
            cref[2] = (:P0, 2)
        end
        add_output!(graph, :P0)
        return (graph, cref)
    else
        s = ceil(Int, sqrt(k))
        v = floor(Int, k / s)
    end

    # Create monomial basis
    prevkey = input
    # up to s-1 are needed for sub-polynomials, and s needed for outer monomials
    for i = 2:s
        key = Symbol("$input$i")
        add_mult!(graph, key, input, prevkey)
        prevkey = key
    end
    Amax = Symbol("$input$s")

    # Create intermediate coefficient-polynomials. (3) and following equation in
    # Fasi 2019
    for i = 0:(v-1)
        # The number of coefficients in the sub-coefficient-polynomial
        # Note, degree is one lower than the index.
        lower_idx = i * s + 1 #plus 1 for 1-indexing
        upper_idx = lower_idx + s - 1 # s number of coeffients, so degree s-1
        idx = lower_idx:upper_idx
        c = view(a, idx)
        ccref = view(cref, idx)

        # B-keys are cumulative sums of the coefficient-polynomials.
        Bkey = Symbol("$(B_base)_$(i)_1")
        add_lincomb!(graph, Bkey, c[1], :I, c[2], input)
        ccref[1] = (Bkey, 1)
        ccref[2] = (Bkey, 2)
        for k = 3:s
            #As above, degree is one lower than the index
            powkey = Symbol("$input$(k-1)")
            prevBkey = Bkey
            Bkey = Symbol("$(B_base)_$(i)_$(k-1)")
            add_lincomb!(graph, Bkey, c[k], powkey, 1.0, prevBkey)
            ccref[k] = (Bkey, 1)
        end
    end

    # Handle highest powers separately. i=v. Part of (3), the following
    # equation, and (4) in Fasi 2019.
    ks = rem(k, s) # Degree of remaining highest order term. We compute B_v_ks
    if ks == 0 # Only constant term
        Pkey = Symbol("$(P_base)$(v-1)")
        Bkey = Symbol("$(B_base)_$(v-1)_$(s-1)")
        # P_{v-1} = a[n]*A^s + B_[v-1}
        add_lincomb!(graph, Pkey, a[n], Amax, 1.0, Bkey)
        cref[n] = (Pkey, 1)
    else # ks > 0. At least constant and linear term.
        # Repeat above procedue with possibly lower-degree polynomial
        idx = (n-(ks+1)+1):n
        c = view(a, idx)
        ccref = view(cref, idx)

        # P_{v-1} = A^s*B_v + B_{v-1}, i.e., have to create B_v

        # B-keys are cumulative sums of the coefficient-polynomials.
        Bkey = Symbol("$(B_base)_$(v)_1")
        add_lincomb!(graph, Bkey, c[1], :I, c[2], input)
        ccref[1] = (Bkey, 1)
        ccref[2] = (Bkey, 2)
        for k = 3:ks+1
            # As above, degree is one lower than the index
            powkey = Symbol("$input$(k-1)")
            prevBkey = Bkey
            Bkey = Symbol("$(B_base)_$(v)_$(k-1)")
            add_lincomb!(graph, Bkey, c[k], powkey, 1.0, prevBkey)
            ccref[k] = (Bkey, 1)
        end
        # C_{v-1} = A^s * B_v
        Ckey = Symbol("$(C_base)$(v-1)")
        Bkey = Symbol("$(B_base)_$(v)_$(ks)")
        add_mult!(graph, Ckey, Bkey, Amax)
        # P_{v-1} = C_{v-1}+ B_{v-1}
        Pkey = Symbol("$(P_base)$(v-1)")
        Bkey = Symbol("$(B_base)_$(v-1)_$(s-1)")
        add_lincomb!(graph, Pkey, 1.0, Ckey, 1.0, Bkey)
    end

    prevPkey = Pkey

    # Evaluate outer-polynomial with Horner-type scheme. (4) in Fasi 2019.
    # Pout = :Nothing;
    Pout = Pkey
    for i in reverse(0:(v-2))
        # C_i = A^s * P_{i+1}
        Ckey = Symbol("$(C_base)$i")
        add_mult!(graph, Ckey, prevPkey, Amax)

        # P_i = C_i + B_i
        Bkey = Symbol("$(B_base)_$(i)_$(s-1)")
        Pkey = Symbol("$(P_base)$i")
        add_lincomb!(graph, Pkey, 1.0, Ckey, 1.0, Bkey)

        prevPkey = Pkey
        Pout = Pkey
    end

    add_output!(graph, Pout)
    return (graph, cref)
end

"""
     (graph,crefs)=graph_ps_degopt(a; input=:A)

Generates the same polynomial as [`graph_ps`](@ref), i.e., the graph for the
Paterson–Stockmeyer procedure with monomial basis coefficieents. However, it
does so by wrapping a call to [`graph_degopt`](@ref), resulting in more degrees
of freedom in `crefs`.
"""
function graph_ps_degopt(a; input = :A)

    # Initial setup
    n = length(a)
    T = eltype(a)
    graph = Compgraph(T)
    k = n - 1 # Polynomial degree

    if k <= 2 # Degenerate cases k = 0, k = 1, and k = 2
        if k == 0
            error("Does not implement degree-zero polynomial.")
        else
            # PS is not more efficient than monomial evaluation for low degrees,
            # but gives complications in code below.
            x = Vector{Tuple{Vector{T},Vector{T}}}(undef, k - 1)
            for i = 1:(k-1)
                x[i] = (
                    vcat(zero(T), one(T), zeros(i - 1)),
                    vcat(zeros(i), one(T)),
                )
            end
            return graph_degopt(x, a)
        end
    end

    s = ceil(Int, sqrt(k))
    v = floor(Int, k / s)
    ks = rem(k, s) # Degree of remaining highest order term.
    # Treat special if only constant. Save one multiplication
    if ks == 0
        i_adj = 1
    else
        i_adj = 0
    end

    x = Vector{Tuple{Vector{T},Vector{T}}}(undef, s + v - 1 - i_adj)
    # Create monomial basis (The first s-1 elements in the x-vector  <=>
    # The first s+1 slots in an x-component)
    for i = 1:(s-1)
        x[i] =
            (vcat(zeros(T, i), one(T)), vcat(zero(T), one(T), zeros(T, i - 1)))
    end

    # Evaluate outer-polynomial with Horner-type scheme, and at the same time
    # create the intermediate coefficient-polynomials
    if ks != 0 # If highest order term non-constant we create the polynomial
        idx = (n-(ks+1)+1):n
        diff = s - (ks + 1)
        c = view(a, idx)
        # A^s*P_{v-1}
        x[s] = (vcat(zeros(T, s), one(T)), vcat(c, zeros(T, diff), zero(T)))
    end

    for i in reverse(0:v-2)
        lower_idx = (i + 1) * s + 1
        upper_idx = lower_idx + s - 1
        idx = lower_idx:upper_idx
        c = view(a, idx)
        # Coeffs for A^s
        xl = vcat(zeros(T, s), one(T), zeros(T, v - 1 - i - i_adj))
        # If highest order term is constant we incorporate it here
        if (i == v - 2) && (ks == 0)
            xr = vcat(c, a[n]) # Coeffs for P_{v-i} + A^s*a[n]
        else
            # Coeffs for P_{v-i} + A^s*P_{v+1-i}
            xr = vcat(c, zero(T), zeros(T, v - 2 - i - i_adj), one(T))
        end
        x[(s-i_adj)+v-1-i] = (xl, xr) # A^s*( P_{v-i} + A^s*P_{v+1-i} )
    end

    z = vcat(a[1:s], zero(T), zeros(T, v - 1 - i_adj), one(T))
    return graph_degopt(x, z, input = input)
end


text:
abstract type AbstractComb{T} <: AbstractArray{T,1} end

"""
    Combination{T} <: AbstractComb{T}

Structure to generate all combinations of the given modules with given length
low_level structure that a user do not need to use, normaly

Fields:
    - `mod`: Array with modulels or where combination can be generated from
    - `len`: lenght of the made constructs
"""
struct Combination{T} <: AbstractComb{T}
    mod::Array{T}
    len::Int
end

"""
    length(c::Combination)
Returns the number of unique combinations that can be made
with a length k and n given modules.
The parameters k and are obtained from the `c` object.
"""
Base.length(c::Combination) = binomial(length(c.mod),c.len);

"""
    eltype(::Combination{T})

Returns the type of the iterator.
"""
Base.eltype(K::Combination{T}) where T = UnorderedConstruct{T} ;
"""
    Base.size(c::Combination)

size needed to display `AbstractArray{T,1}`, but has no furder use
"""
Base.size(c::Combination) = (length(c),1) ;
"""
    iterate(c::Combination{T} where T, state = [i for i in c.len:-1:1])

Iterater for Combinations, printed in lexicographical order.
"""
function Base.iterate(c::Combination{T} where T, state = [i for i in c.len:-1:1])
    n = length(c.mod)
    # check when done, last constructed is printed, stop Iterator
    if state[1] > n
        return nothing
    end

    # special case len = 1
    if c.len == 1
        return _len1(c,state)
    end

    # case len > 1

    # make the current construct based on given state
    construct = sum(c.mod[state])

    # generated the new state

    # i start at the last index, length of the construct
    i = c.len
    #evaluated if the difference between the last and the next is higher than 1 (than their is still a construct that can be made)
    # if not than lower the index with one and recheck if the difference is 1 . i>1 prevents error
    while i > 1  && state[i-1] - state[i] == 1
        i -=1
    end
    # upate the state at position i and reset all higher index to lowest values
    if i < c.len
        state = _restate(state,i)
    else
        state[i] +=1
    end
    return (construct,state)

end

"""
    _restate(state, i)

Update state at index i and reset all values behind to lowsted allow values.
function used in the Base.iterate(c::Combination{T} where T, state = [i for i in c.len:-1:1])
"""
function _restate(state, i)
    state[i] +=1;
    state[(i+1):end] = collect(length(state)-i:-1:1);
    return state
end

"""
    _len1(c,state)

part of Base.iterate(c::Combination{T} where T, state = [i for i in c.len:-1:1]),
used to generated all possibilities when length of the allowed constructs = 1.

"""
function _len1(c,state)
        construct = [c.mod[state[1]]]
        state[1] +=1
        return (UnorderedConstruct(construct), state)
end

"""
    function getindex(c::Combination,pos::Int)

Get the index of `c` at position `pos`.

Index based on rank Combination.
https://en.wikipedia.org/wiki/Combinatorial_number_system
"""
function getindex(c::Combination,pos::Int)

    # spacial case of lenght is 1
    if c.len == 1
        return UnorderedConstruct([c.mod[pos]])
    end
    # if not 1
    com = Array{Int}(undef, c.len, 1)
    for i in c.len:-1:1
        n = i-1
        b_val = binomial(n,i)
        while b_val < pos
            n += 1
            b_val = binomial(n,i)
        end
        pos -= binomial(n-1,i)
        com[c.len-i+1] = n-1
    end
    index =  com .+=1
    return (sum(c.mod[index]))
end;

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
using Test
using RDates
using Dates
using ParserCombinator

cal_mgr = SimpleCalendarManager()
setcalendar!(cal_mgr, "WEEKEND", WeekendCalendar())
setcalendar!(cal_mgr, "WEEK/END", calendar(cal_mgr, "WEEKEND"))
setcalendar!(cal_mgr, "WEEK-END", calendar(cal_mgr, "WEEKEND"))
setcalendar!(cal_mgr, "WEEK END", calendar(cal_mgr, "WEEKEND"))
setcachedcalendar!(cal_mgr, "CACHED WEEKEND", calendar(cal_mgr, "WEEKEND"))

@testset "RDates" verbose=true begin
    @testset "Calendar Holidays" begin
        @test holidays(calendar(cal_mgr, "WEEKEND"), Date(2021,7,8), Date(2021,7,11)) == Bool[0, 0, 1, 1]
        @test holidays(calendar(cal_mgr, "WEEKEND"), Date(2021,7,8), Date(2021,7,9)) == Bool[0, 0]
        @test holidaycount(calendar(cal_mgr, "WEEKEND"), Date(2021,7,4), Date(2021,7,20)) == 5
        @test bizdaycount(calendar(cal_mgr, "WEEKEND"), Date(2021,7,4), Date(2021,7,20)) == 12

        @test holidays(calendar(cal_mgr, "CACHED WEEKEND"), Date(2021,7,8), Date(2021,7,11)) == Bool[0, 0, 1, 1]
        @test holidays(calendar(cal_mgr, "CACHED WEEKEND"), Date(2021,7,8), Date(2021,7,9)) == Bool[0, 0]
        @test holidaycount(calendar(cal_mgr, "CACHED WEEKEND"), Date(2021,7,4), Date(2021,7,20)) == 5
        @test bizdaycount(calendar(cal_mgr, "CACHED WEEKEND"), Date(2021,7,4), Date(2021,7,20)) == 12
        @test holidays(calendar(cal_mgr, "CACHED WEEKEND"), Date(2031,1,30), Date(2031,2,3)) == Bool[0, 0, 1, 1, 0]
        @test holidays(calendar(cal_mgr, "CACHED WEEKEND"), Date(2011,2,3), Date(2011,2,6)) == Bool[0, 0, 1, 1]
    end

    @testset "Next and Previous" begin
        @test rd"Next(0E, 1E)" + Date(2021,1,1) == Date(2021,4,4)
        @test rd"Next(0E, 1E)" + Date(2021,4,3) == Date(2021,4,4)
        @test rd"Next(0E, 1E)" + Date(2021,4,4) == Date(2022,4,17)
        @test rd"Next!(0E, 1E)" + Date(2021,4,4) == Date(2021,4,4)
        @test rd"Next!(0E, 1E)" + Date(2021,4,5) == Date(2022,4,17)
        @test rd"Next(0E, -1E)" + Date(2021,1,1) == Date(2021,4,4)
        @test_throws ErrorException rd"Next(0E, -1E)" + Date(2021,12,31)

        @test rd"Previous(0E, -1E)" + Date(2021,1,1) == Date(2020,4,12)
        @test rd"Previous(0E, -1E)" + Date(2020,4,13) == Date(2020,4,12)
        @test rd"Previous(0E, -1E)" + Date(2020,4,12) == Date(2019,4,21)
        @test rd"Previous!(0E, -1E)" + Date(2020,4,12) == Date(2020,4,12)
        @test rd"Previous!(0E, -1E)" + Date(2020,4,11) == Date(2019,4,21)
        @test rd"Previous(0E, 1E)" + Date(2021,12,31) == Date(2021,4,4)
        @test_throws ErrorException rd"Previous(0E, 1E)" + Date(2021,1,1)

        @test -rd"Next(1d)" == rd"Previous(-1d)"
        @test rd"-Next(1d)" == rd"Previous(-1d)"
        @test -2*rd"Next(1d)" == rd"Previous(-2d)"
        @test rd"-2*Next(1d)" == rd"Previous(-2d)"

        @test -rd"Previous(-1d)" == rd"Next(1d)"
        @test rd"-Previous(-1d)" == rd"Next(1d)"
        @test -2*rd"Previous(-1d)" == rd"Next(2d)"
        @test rd"-2*Previous(-1d)" == rd"Next(2d)"
    end

    @testset "Rounding Conventions" begin
        @test apply(rd"0d@WEEKEND[NBD]", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0d@WEEKEND[NBD]", Date(2021,7,11), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0d@WEEKEND[NBDSM]", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0d@WEEKEND[NBDSM]", Date(2021,7,11), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0d@WEEKEND[PBD]", Date(2021,7,10), cal_mgr) == Date(2021,7,9)
        @test apply(rd"0d@WEEKEND[PBD]", Date(2021,7,11), cal_mgr) == Date(2021,7,9)
        @test apply(rd"0d@WEEKEND[PBDSM]", Date(2021,7,10), cal_mgr) == Date(2021,7,9)
        @test apply(rd"0d@WEEKEND[PBDSM]", Date(2021,7,11), cal_mgr) == Date(2021,7,9)
        @test apply(rd"0d@WEEKEND[NR]", Date(2021,7,10), cal_mgr) == Date(2021,7,10)
        @test apply(rd"0d@WEEKEND[NR]", Date(2021,7,11), cal_mgr) == Date(2021,7,11)
        @test apply(rd"0d@WEEKEND[NEAR]", Date(2021,7,10), cal_mgr) == Date(2021,7,9)
        @test apply(rd"0d@WEEKEND[NEAR]", Date(2021,7,11), cal_mgr) == Date(2021,7,12)
        # same month special casing
        @test apply(rd"0d@WEEKEND[NBDSM]", Date(2021,7,31), cal_mgr) == Date(2021,7,30)
        @test apply(rd"0d@WEEKEND[NBDSM]", Date(2021,8,1), cal_mgr) == Date(2021,8,2)
        @test apply(rd"0d@WEEKEND[PBDSM]", Date(2021,7,31), cal_mgr) == Date(2021,7,30)
        @test apply(rd"0d@WEEKEND[PBDSM]", Date(2021,8,1), cal_mgr) == Date(2021,8,2)
    end

    @testset "Calendar Formats" begin
        @test apply(rd"0d@WEEKEND[NBD]", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0d@WEEK/END[NBD]", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0d@WEEK END[NBD]", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0d@WEEK-END[NBD]", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0b@WEEKEND", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0b@WEEK/END", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0b@WEEK END", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
        @test apply(rd"0b@WEEK-END", Date(2021,7,10), cal_mgr) == Date(2021,7,12)
    end

    @testset "Calendar Adjustments" begin
        cal_mgr = RDates.SimpleCalendarManager(Dict("WEEKEND" => RDates.WeekendCalendar()))
        @test apply(RDates.CalendarAdj(["WEEKEND"], rd"1d", RDates.HolidayRoundingNBD()), Date(2019,4,16), cal_mgr) == Date(2019,4,17)
        @test apply(RDates.CalendarAdj(["WEEKEND"], rd"1d", RDates.HolidayRoundingNBD()), Date(2019,9,27), cal_mgr) == Date(2019,9,30)
        @test apply(RDates.CalendarAdj(["WEEKEND"], rd"1d", RDates.HolidayRoundingNBD()), Date(2019,11,29), cal_mgr) == Date(2019,12,2)
        @test apply(RDates.CalendarAdj(["WEEKEND"], rd"1d", RDates.HolidayRoundingNBDSM()), Date(2019,11,29), cal_mgr) == Date(2019,11,29)
    end

    @testset "Business Days" begin
        @test apply(rd"0b@WEEKEND", Date(2019,9,28), cal_mgr) == Date(2019,9,30)
        @test apply(rd"-0b@WEEKEND", Date(2019,9,28), cal_mgr) == Date(2019,9,27)
        @test apply(rd"1b@WEEKEND", Date(2019,9,28), cal_mgr) == Date(2019,10,1)
        @test apply(rd"-1b@WEEKEND", Date(2019,9,28), cal_mgr) == Date(2019,9,26)
        @test apply(rd"2b@WEEKEND", Date(2019,9,28), cal_mgr) == Date(2019,10,2)
    end

    @testset "Add Ordering" begin
        @test rd"1d" + Date(2019,4,16) == Date(2019,4,17)
        @test Date(2019,4,16) + rd"1d" == Date(2019,4,17)
        @test rd"1d+2d" == rd"1d" + rd"2d"
        @test rd"3*1d" == rd"3d"
    end

    @testset "Month and Year Conventions" begin
        @test RDates.Year(1) + Date(2016,2,29) == Date(2017,2,28)
        @test Dates.Year(1) + Date(2016,2,29) == Date(2017,2,28)
        @test rd"1y" + Date(2016,2,29) == Date(2017,2,28)
        @test RDates.Year(1, RDates.InvalidDayLDOM(), RDates.MonthIncrementPDOM()) + Date(2016,2,29) == Date(2017,2,28)
        @test rd"1y[LDOM;PDOM]" + Date(2016,2,29) == Date(2017,2,28)
        @test RDates.Year(1, RDates.InvalidDayFDONM(), RDates.MonthIncrementPDOM()) + Date(2016,2,29) == Date(2017,3,1)
        @test rd"1y[FDONM;PDOM]" + Date(2016,2,29) == Date(2017,3,1)

        @test RDates.Month(1) + Date(2019,1,31) == Date(2019,2,28)
        @test Dates.Month(1) + Date(2019,1,31) == Date(2019,2,28)
        @test rd"1m" + Date(2019,1,31) == Date(2019,2,28)
        @test RDates.Month(1, RDates.InvalidDayNDONM(), RDates.MonthIncrementPDOM()) + Date(2019,1,31) == Date(2019,3,3)
        @test rd"1m[NDONM;PDOM]" + Date(2019,1,31) == Date(2019,3,3)
        @test RDates.Month(1, RDates.InvalidDayFDONM(), RDates.MonthIncrementPDOM()) + Date(2019,1,31) == Date(2019,3,1)
        @test rd"1m[FDONM;PDOM]" + Date(2019,1,31) == Date(2019,3,1)
        @test RDates.Month(1, RDates.InvalidDayNDONM(), RDates.MonthIncrementPDOMEOM()) + Date(2019,1,31) == Date(2019,2,28)
        @test rd"1m[NDONM;PDOMEOM]" + Date(2019,1,31) == Date(2019,2,28)
        @test RDates.Month(1, RDates.InvalidDayNDONM(), RDates.MonthIncrementPDOMEOM()) + Date(2019,1,30) == Date(2019,3,2)
        @test rd"1m[NDONM;PDOMEOM]" + Date(2019,1,30) == Date(2019,3,2)

        # PDOMEOM can also support calendars to work off the last business day of the month
        @test apply(rd"3m[LDOM;PDOMEOM@WEEKEND]", Date(2019,8,30), cal_mgr) == Date(2019,11,29)
    end

    @testset "Parsing Whitespace" begin
        @test rd"  1d" == rd"1d"
        @test rd"1d  + 3d" == rd"1d+3d"
        @test rd"--1d" == rd"1d"
    end

    @testset "Days" begin
        @test rd"1d" + Date(2019,4,16) == Date(2019,4,17)
        @test rd"1d" + Date(2019,4,30) == Date(2019,5,1)
        @test rd"0d" + Date(2015,3,23) == Date(2015,3,23)
        @test rd"7d" + Date(2017,10,25) == Date(2017,11,1)
        @test rd"-1d" + Date(2014,1,1) == Date(2013,12,31)
    end

    @testset "Weeks" begin
        @test rd"1w" + Date(2019,4,16) == Date(2019,4,23)
        @test rd"1w" + Date(2019,4,30) == Date(2019,5,7)
        @test rd"0w" + Date(2015,3,23) == Date(2015,3,23)
        @test rd"7w" + Date(2017,10,25) == Date(2017,12,13)
        @test rd"-1w" + Date(2014,1,1) == Date(2013,12,25)
    end

    @testset "Months" begin
        @test rd"1m" + Date(2019,4,16) == Date(2019,5,16)
        @test rd"1m" + Date(2019,4,30) == Date(2019,5,30)
        @test rd"0m" + Date(2015,3,23) == Date(2015,3,23)
        @test rd"12m" + Date(2017,10,25) == Date(2018,10,25)
        @test rd"-1m" + Date(2014,1,1) == Date(2013,12,1)
    end

    @testset "Years" begin
        @test rd"1y" + Date(2019,4,16) == Date(2020,4,16)
        @test rd"1y" + Date(2019,4,30) == Date(2020,4,30)
        @test rd"0m" + Date(2015,3,23) == Date(2015,3,23)
        @test rd"12y" + Date(2017,10,25) == Date(2029,10,25)
        @test rd"-1y" + Date(2014,1,1) == Date(2013,1,1)
    end

    @testset "Day Months" begin
        @test rd"12MAR" + Date(2018,3,3) == Date(2018,3,12)
        @test rd"1JAN" + Date(2019,12,31) == Date(2019,1,1)
        @test rd"1JAN" + Date(2020,1,1) == Date(2020,1,1)
    end

    @testset "Easters" begin
        @test rd"0E" + Date(2018,3,3) == Date(2018,4,1)
        @test rd"0E" + Date(2018,12,3) == Date(2018,4,1)
        @test rd"1E" + Date(2018,3,3) == Date(2019,4,21)
        @test rd"-1E" + Date(2018,3,3) == Date(2017,4,16)
    end

    @testset "Weekdays" begin
        @test rd"1MON" + Date(2017,10,25) == Date(2017,10,30)
        @test rd"10SAT" + Date(2017,10,25) == Date(2017,12,30)
        @test rd"-1WED" + Date(2017,10,25) == Date(2017,10,18)
        @test rd"-10FRI" + Date(2017,10,25) == Date(2017,8,18)
        @test rd"-1TUE" + Date(2017,10,25) == Date(2017,10,24)
    end

    @testset "Dates" begin
        @test rd"1JAN2019" + Date(2017,10,25) == Date(2019,1,1)
        @test RDates.Date(Date(2019,1,1)) + Date(2017,10,25) == Date(2019,1,1)
    end

    @testset "Non Supported Negations" begin
        @test rd"1d" - rd"1st MON" == rd"1d" + rd"1st MON"
    end

    @testset "Nth Weekdays" begin
        @test rd"1st MON" + Date(2017,10,25) == Date(2017,10,2)
        @test rd"2nd FRI" + Date(2017,10,25) == Date(2017,10,13)
        @test rd"4th SAT" + Date(2017,11,25) == Date(2017,11,25)
        @test rd"5th SUN" + Date(2017,12,25) == Date(2017,12,31)
    end

    @testset "Nth Last Weekdays" begin
        @test rd"Last MON" + Date(2017,10,24) == Date(2017,10,30)
        @test rd"2nd Last FRI" + Date(2017,10,24) == Date(2017,10,20)
        @test rd"5th Last SUN" + Date(2017,12,24) == Date(2017,12,3)
    end

    @testset "Bad Nth Weekdays" begin
        @test_throws ArgumentError rd"5th WED" + Date(2017,10,25)
        @test_throws ArgumentError rd"5th MON" + Date(2017,6,1)
    end

    @testset "First/Last Day of Month" begin
        @test rd"FDOM" + Date(2017,10,25) == Date(2017,10,1)
        @test rd"LDOM" + Date(2017,10,25) == Date(2017,10,31)

        # FDOM and LDOM can support calendars as well.
        @test apply(rd"FDOM@WEEKEND", Date(2019,9,25), cal_mgr) == Date(2019,9,2)
        @test apply(rd"LDOM@WEEKEND", Date(2019,8,25), cal_mgr) == Date(2019,8,30)
    end

    @testset "Basic Compounds" begin
        @test rd"1d+1d" + Date(2017,10,26) == Date(2017,10,28)
        @test rd"2*1d" + Date(2017,10,26) == Date(2017,10,28)
        @test rd"1d+1d+1d+1d" + Date(2017,10,26) == Date(2017,10,30)
        @test rd"4*1d" + Date(2017,10,26) == Date(2017,10,30)
        @test rd"2*2d" + Date(2017,10,26) == Date(2017,10,30)
        @test rd"1d-1d+1d-1d" + Date(2017,10,26) == Date(2017,10,26)
        @test rd"2*3d+1d" + Date(2019,4,16) == Date(2019,4,23)
        @test rd"2*(3d+1d)" + Date(2019,4,16) == Date(2019,4,24)
        @test rd"2d-1E" + Date(2019,4,16) == Date(2018,4,1)
        @test rd"1d" - rd"2*1d" + Date(2019,5,1) == Date(2019,4,30)
        @test RDates.multiply(rd"1d", 3) + Date(2017,4,14) == Date(2017,4,17)
    end

    @testset "Roll vs No-Roll Multiplication" begin
        # If we apply each 1m addition individually, then invalid days in Feb kicks in.
        @test rd"1m+1m" + Date(2019,1,31) == Date(2019,3,28)
        @test rd"2*Repeat(1m)" + Date(2019,1,31) == Date(2019,3,28)
        # Normal multiplication will embed this in the period instead.
        @test rd"2m" + Date(2019,1,31) == Date(2019,3,31)
        @test rd"2*1m" + Date(2019,1,31) == Date(2019,3,31)
    end

    @testset "Ranges" begin
        @test collect(range(Date(2017,1,1), Date(2017,1,3), rd"1d")) == [Date(2017,1,1), Date(2017,1,2), Date(2017,1,3)]
        @test collect(range(Date(2017,1,1), Date(2017,1,3), rd"2d")) == [Date(2017,1,1), Date(2017,1,3)]
        @test collect(range(Date(2017,1,1), Date(2017,1,18), rd"1d+1w")) == [Date(2017,1,1), Date(2017,1,9), Date(2017,1,17)]
        @test collect(range(Date(2017,1,1), Date(2017,1,3), rd"1d", inc_from=false)) == [Date(2017,1,2), Date(2017,1,3)]
        @test collect(range(Date(2017,1,1), Date(2017,1,3), rd"1d", inc_to=false)) == [Date(2017,1,1), Date(2017,1,2)]
        @test collect(range(Date(2017,1,1), Date(2017,1,3), rd"1d", inc_from=false, inc_to=false)) == [Date(2017,1,2)]
        @test collect(range(Date(2017,1,1), Date(2018,1,1), rd"1DEC2017 + 3m + 3rd WED")) == [Date(2017,3,15), Date(2017,6,21), Date(2017,9,20), Date(2017,12,20)]
    end

    @testset "Parsing Methods" begin
        @test rdate("1d") == rd"1d"
        @test rdate("3*2d") == rd"3*2d"
        @test rd"3*2d" == 3*rd"2d"
    end

    @testset "Failed Parsing" begin
        @test_throws ParserCombinator.ParserException rdate("1*2")
        @test_throws ParserCombinator.ParserException rdate("1dw")
        @test_throws ParserCombinator.ParserException rdate("+2d")
        @test_throws ParserCombinator.ParserException rdate("2d+")
        @test_throws ParserCombinator.ParserException rdate("d")
    end

    @testset "String Forms" begin
        @test string(rd"1d") == "1d"
        @test string(rd"3d + 2w") == "17d"
        @test string(rd"2*(1m + 1y)") == "2m[LDOM;PDOM]+2y[LDOM;PDOM]"
        @test string(rd"2*Repeat(1m)") == "2*Repeat(1m[LDOM;PDOM])"
        @test string(rd"1m[LDOM;PDOMEOM@A]@B[NBD]") == "(1m[LDOM;PDOMEOM@A])@B[NBD]"
        @test string(rd"(1m+2d)@A[PBD]") == "(1m[LDOM;PDOM]+2d)@A[PBD]"
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c2e4a0d0-5cae-11eb-199f-3b8bc7ca4867
begin
	using ReinforcementLearning
	using Flux
	using Statistics
	using Plots
	using SparseArrays
end

# ╔═╡ 65614b3c-5cae-11eb-0613-93393f90cfc0
md"""
# Chapter 10 Mountain Car
"""

# ╔═╡ 054805aa-5caf-11eb-1f4a-cd04a590a023
md"""
The `MountainCarEnv` is already provided in `ReinforcementLearning.jl`. So we can use it directly here. Note that by default this environment will terminate at the maximum step of `200`. While in the example on the book there's no such restriction.
"""

# ╔═╡ db653a64-5cae-11eb-0a6d-cd071416fc50
env = MountainCarEnv()

# ╔═╡ dc826048-5cae-11eb-3141-1b589c25852a
S = state_space(env)

# ╔═╡ 02252420-5caf-11eb-2352-251cf8b60517
md"""
First let's define a `Tiling` structure to encode the state.
"""

# ╔═╡ e1e987aa-5cae-11eb-1f24-69fefb87d4e0
begin
	struct Tiling{N,Tr<:AbstractRange}
		ranges::NTuple{N,Tr}
		inds::LinearIndices{N,NTuple{N,Base.OneTo{Int}}}
	end

	Tiling(ranges::AbstractRange...) =Tiling(
		ranges,
		LinearIndices(Tuple(length(r) - 1 for r in ranges))
	)

	Base.length(t::Tiling) = reduce(*, (length(r) - 1 for r in t.ranges))

	function Base.:-(t::Tiling, xs)
		Tiling((r .- x for (r, x) in zip(t.ranges, xs))...)
	end

	encode(range::AbstractRange, x) = floor(Int, div(x - range[1], step(range)) + 1)

	encode(t::Tiling, xs) = t.inds[CartesianIndex(Tuple(map(encode,  t.ranges, xs)))]
end

# ╔═╡ a0f8a23c-5caf-11eb-01d2-1706e563dac1
begin
	ntilings = 8
	ntiles = 8
	tiling = Tiling(
		(
			range(r.left, step=(r.right-r.left)/ntiles, length=ntiles+2)
			for r in S
		)...
	)
	offset = map(x-> x.right - x.left, S) ./ (ntiles * ntilings)
	tilings = [tiling - offset .* (i-1) for i in 1:ntilings]
end

# ╔═╡ d62c2d40-5caf-11eb-36e4-5bb8ee95c24e
md"""
The rest parts are simple, we initialize `agent` and `env`, then roll out experiments:
"""

# ╔═╡ 214a577a-5cb0-11eb-356b-ab06ad454eef
function create_env_agent(α=2e-4, n=0)
    env = StateOverriddenEnv(
        MountainCarEnv(;max_steps=10000),
        s -> sparse(map(t -> encode(t, s), tilings), 1:8, ones(8), 81, 8) |> vec
    )

    agent = Agent(
        policy=QBasedPolicy(
            learner=TDLearner(
                approximator=LinearQApproximator(
                    n_state=81*8,
                    n_action=3,
                    opt = Descent(α)
                    ),
                method=:SARSA,
                n=n
                ),
            explorer=GreedyExplorer()
            ),
        trajectory=VectorSARTTrajectory(;state=Vector{Int})
    )

    env, agent
end

# ╔═╡ 3901d7b4-5cb0-11eb-2cc3-9ff34896f465
X = range(S[1].left, stop=S[1].right, length=40)

# ╔═╡ 3c0e7bea-5cb0-11eb-3190-598e56c887ff
Y = range(S[2].left, stop=S[2].right, length=40)

# ╔═╡ 230ee9ae-5cb0-11eb-0fe5-6fd8f5482af4
function show_approximation(n)
    env, agent = create_env_agent()
    run(agent, env, StopAfterEpisode(n))
    [
		agent.policy.learner.approximator(env.f([p, v])) |> maximum
        for p in X, v in Y
	]
end

# ╔═╡ 4ba11496-5cb0-11eb-202c-7f2baba059bf
n = 10

# ╔═╡ 4905a328-5cb0-11eb-03f1-0f6f27350ade
plot(X, Y, -show_approximation(n), linetype=:wireframe)

# ╔═╡ 68779ca2-5cb0-11eb-11e4-9fdb4ea0989f
begin
	fig_10_2 = plot(legend=:topright)
	n_runs = 5  # quite slow here, need revisit
	for α in [0.1/8, 0.2/8, 0.5/8]
		avg_steps_per_episode = zeros(500)
		for _ in 1:n_runs
			env, agent = create_env_agent(α)
			hook = StepsPerEpisode()
			run(agent, env, StopAfterEpisode(500; is_show_progress=false),hook)
			avg_steps_per_episode .+= hook.steps
		end
		plot!(fig_10_2, avg_steps_per_episode ./ n_runs, yscale=:log10, label="α=$α")
	end
	fig_10_2
end

# ╔═╡ 79e4a002-5cb0-11eb-049a-db8707d6c400
begin
	function run_once(α, n;is_reduce=true, n_episode=50)
		env, agent = create_env_agent(α, n)
		hook = StepsPerEpisode()
		run(agent, env, StopAfterEpisode(n_episode; is_show_progress=false),hook)
		is_reduce ? mean(hook.steps) : hook.steps
	end
	fig_10_3 = plot()
	plot!(fig_10_3, mean(run_once(0.5/8, 1; is_reduce=false, n_episode=500) for _ in 1:10), yscale=:log10)
	plot!(fig_10_3, mean(run_once(0.3/8, 8; is_reduce=false, n_episode=500) for _ in 1:10), yscale=:log10)
	fig_10_3
end

# ╔═╡ 832c5ce0-5cb0-11eb-06f6-4f9ad5b707a2
begin
	fig_10_4 = plot(legend=:topright)
	for (A, n) in [(0.4:0.1:1.7, 1), (0.3:0.1:1.6, 2), (0.2:0.1:1.4, 4), (0.2:0.1:0.9, 8), (0.2:0.1:0.7, 16)]
		plot!(fig_10_4, A, [mean(run_once(α/8, n) for _ in 1:5) for α in A], label="n = $n")
	end
	fig_10_4
end

# ╔═╡ Cell order:
# ╟─65614b3c-5cae-11eb-0613-93393f90cfc0
# ╠═c2e4a0d0-5cae-11eb-199f-3b8bc7ca4867
# ╟─054805aa-5caf-11eb-1f4a-cd04a590a023
# ╠═db653a64-5cae-11eb-0a6d-cd071416fc50
# ╠═dc826048-5cae-11eb-3141-1b589c25852a
# ╟─02252420-5caf-11eb-2352-251cf8b60517
# ╠═e1e987aa-5cae-11eb-1f24-69fefb87d4e0
# ╠═a0f8a23c-5caf-11eb-01d2-1706e563dac1
# ╟─d62c2d40-5caf-11eb-36e4-5bb8ee95c24e
# ╠═214a577a-5cb0-11eb-356b-ab06ad454eef
# ╠═3901d7b4-5cb0-11eb-2cc3-9ff34896f465
# ╠═3c0e7bea-5cb0-11eb-3190-598e56c887ff
# ╠═230ee9ae-5cb0-11eb-0fe5-6fd8f5482af4
# ╠═4ba11496-5cb0-11eb-202c-7f2baba059bf
# ╠═4905a328-5cb0-11eb-03f1-0f6f27350ade
# ╠═68779ca2-5cb0-11eb-11e4-9fdb4ea0989f
# ╠═79e4a002-5cb0-11eb-049a-db8707d6c400
# ╠═832c5ce0-5cb0-11eb-06f6-4f9ad5b707a2


text:
using Documenter, Sherlock

makedocs(;
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/caseykneale/Sherlock.jl/blob/{commit}{path}#L{line}",
    sitename="Sherlock.jl",
    authors="Casey Kneale",
    assets=String[],
)

deploydocs(;
    repo="github.com/caseykneale/Sherlock.jl.git",
    devbranch = "master",
)