text:
"""
```
prepare_forecast_inputs!(m, input_type, cond_type, output_vars;
    df = DataFrame(), verbose = :none)
```

Add required outputs using `add_requisite_output_vars` and load data if
necessary.

### Inputs

- `m::AbstractDSGEModel`: model object
- `input_type::Symbol`: See `?forecast_one`.
- `cond_type::Symbol`: See `?forecast_one`.
- `output_vars::Vector{Symbol}`: vector of desired output variables. See
  `?forecast_one_draw`

### Keyword Arguments

- `df::DataFrame`: historical data. If `cond_type in [:semi, :full]`, then the
   final row of `df` should be the period containing conditional data. If not
   provided, then `df` will be loaded using `load_data` with the appropriate
   `cond_type`
- `only_filter::Bool`: do not run the smoother and only run the filter. This limits the number of
  output variables which can be calculated.
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`

### Outputs

- `output_vars`
- `df`
"""
function prepare_forecast_inputs!(m::AbstractDSGEModel{S},
    input_type::Symbol, cond_type::Symbol, output_vars::Vector{Symbol};
    df::DataFrame = DataFrame(), subset_inds::AbstractRange{Int64} = 1:0,
    check_empty_columns::Bool = true, bdd_fcast::Bool = true, only_filter::Bool = false,
    verbose::Symbol = :none) where {S<:AbstractFloat}

    @assert cond_type in [:none, :semi, :full] "cond_type must be one of :none, :semi, or :full"

    if only_filter
        smooth_vars = [:histstates, :histpseudo, :histshocks, :histstdshocks, :forecastshocks, :forecaststdshocks,
                       :shockdecstates, :shockdecpseudo, :shockdecobs, :dettrendstates, :dettrendpseudo, :dettrendobs,
                       :shockdecseqstates, :shockdecseqpseudo, :shockseqdecobs]
        @assert isempty(intersect(output_vars, smooth_vars)) "Some output variables requested cannot be use with keyword only_filter"
    end

    implied_horizons = subtract_quarters(date_forecast_end(m), date_forecast_start(m)) + 1
    if implied_horizons != get_setting(m, :forecast_horizons)
        @warn "Setting :forecast_horizons does not match the number of periods implied by Settings :date_forecast_end" *
            " and :date_forecast_start. Updating :forecast_horizons to match"
        m <= Setting(:forecast_horizons, implied_horizons, "Number of periods to forecast ahead")
    end

    # Compute everything that will be needed to plot original output_vars
    output_vars = add_requisite_output_vars(output_vars, bdd_fcast = bdd_fcast)
    if input_type == :prior
        output_vars = setdiff(output_vars, [:bddforecastobs])
    end

    # Get products and classes computed
    output_prods   = unique(map(get_product, output_vars))
    output_classes = unique(map(get_class,   output_vars))

    # Throw error if input_type = :subset but no subset_inds provided
    if input_type == :subset && isempty(subset_inds)
        error("Must supply nonempty subset_inds if input_type = :subset")
    end

    # Determine if we are only running IRFs. If so, we won't need to load data below
    irfs_only = all(prod -> prod == :irf, output_prods)

    # Load data if not provided, else check data well-formed
    if !irfs_only
        if isempty(df)
            data_verbose = verbose == :none ? :none : :low
            df = load_data(m; cond_type = cond_type, try_disk = true, check_empty_columns = check_empty_columns, verbose = data_verbose)
        else
            @assert df[1, :date] == date_presample_start(m)
            @assert df[end, :date] == (cond_type == :none ? date_mainsample_end(m) : date_conditional_end(m))
        end
    end

    # If regime switching, check settings are consistent
    regswitch = haskey(get_settings(m), :regime_switching) ? get_setting(m, :regime_switching) : false
    if regswitch
        @assert get_setting(m, :n_regimes) == length(get_setting(m, :regime_dates)) "The number" *
            " in Setting :n_regimes ($(string(get_setting(m, :n_regimes))))" *
            " does not match the length of Setting :regime_dates ($(string(length(get_setting(m, :regime_dates)))))."
        @assert get_setting(m, :regime_dates)[1] == date_presample_start(m) "The first regime" *
            " date ($(string(get_setting(m, :regime_dates)[1]))) must match the presample start date " *
            "($(string(date_presample_start(m))))."
    end

    # Print out settings of temporary policies
    if verbose == :high && regswitch
        @info "Regime-Switching and Temporary Policy Settings"
        println("n_regimes = $(get_setting(m, :n_regimes))")
        println("reg_forecast_start = $(get_setting(m, :reg_forecast_start))")
        println("reg_post_conditional_end = $(get_setting(m, :reg_post_conditional_end))")
        println("time-varying information set: ", haskey(get_settings(m), :tvis_information_set) ?
                get_setting(m, :tvis_information_set) : false)

        println("gensys2: " * string(haskey(get_settings(m), :gensys2) ? get_setting(m, :gensys2) : false))
        replace_eqcond = haskey(get_settings(m), :replace_eqcond) ? get_setting(m, :replace_eqcond) : false
        println("replace_eqcond: " * string(replace_eqcond))
        if replace_eqcond
            replace_regs = sort!(collect(keys(get_setting(m, :regime_eqcond_info))))
            println("regime_eqcond_info: ", OrderedDict(k => get_setting(m, :regime_eqcond_info)[k] for k in replace_regs))
        end

        println("uncertain_altpolicy: " *
                string(haskey(get_settings(m), :uncertain_altpolicy) ? get_setting(m, :uncertain_altpolicy) : false))
        println("uncertain_temporary_altpolicy: " * string(haskey(get_settings(m), :uncertain_temporary_altpolicy) ? get_setting(m, :uncertain_temporary_altpolicy) : false))
        if replace_eqcond
            println("time-varying credibility: " * string(any(v -> !ismissing(v.weights), values(get_setting(m, :regime_eqcond_info)))))
        end
        if haskey(get_settings(m), :uncertain_altpolicy) ? get_setting(m, :uncertain_altpolicy) : false
            println("Desired policy rule: " *
                    string(haskey(get_settings(m), :regime_eqcond_info) ?
                           get_setting(m, :regime_eqcond_info)[maximum(collect(keys(get_setting(m, :regime_eqcond_info))))] :
                           alternative_policy(m).key))
            println("Other policy rules: " * join([string(x.key) for x in get_setting(m, :alternative_policies)], ", "))
            if haskey(get_settings(m), :imperfect_awareness_weights) # does not work currently, need to update
                println("Credibility weights: ", get_setting(m, :imperfect_awareness_weights))
            else
                # Only other way to specify imperfect awareness is regime_eqcond_info, so it must exist
                sorted_eqcond_info_regs = sort!(collect(keys(get_setting(m, :regime_eqcond_info))))
                println("Credbility weights: ", OrderedDict(reg => get_setting(m, :regime_eqcond_info)[reg].weights for reg in sorted_eqcond_info_regs))
            end
        end
        if haskey(get_settings(m), :temporary_altpolicy_length)
            println("temporary altpol length = ", get_setting(m, :temporary_altpolicy_length))
        end
        altkey = alternative_policy(m).key
        if altkey in [:smooth_ait_gdp_alt, :smooth_ait_gdp, :flexible_ait]
            println("skip_altpolicy_state_init: ", haskey(get_settings(m), :skip_altpolicy_state_init) ?
                get_setting(m, :skip_altpolicy_state_init) : false)
            println("add_initialize_pgap_ygap_pseudoobs: ", haskey(get_settings(m), :add_initialize_pgap_ygap_pseudoobs) ?
                get_setting(m, :add_initialize_pgap_ygap_pseudoobs) : false)
            println("φ_π = ", get_setting(m, Symbol(altkey, "_φ_π")))
            println("φ_y = ", get_setting(m, Symbol(altkey, "_φ_y")))
            println("ait_Thalf = ", get_setting(m, :ait_Thalf))
            println("gdp_Thalf = ", get_setting(m, :gdp_Thalf))
            println("ρ_smooth = ", get_setting(m, Symbol(altkey, "_ρ_smooth")))
            println("pgap type = ", get_setting(m, :pgap_type))
            println("ygap type = ", get_setting(m, :ygap_type))
            println("initial pgap = ", get_setting(m, :pgap_value))
            println("initial ygap = ", get_setting(m, :ygap_value))
        end

        if haskey(get_settings(m), :temporary_altpolicy_names)
            if :zero_rate in get_setting(m, :temporary_altpolicy_names)
                println("zero_rate_zlb_value = ", haskey(get_settings(m), :zero_rate_zlb_value) ?
                        get_setting(m, :zero_rate_zlb_value) : round(.1 / 4., digits = 3))
            end
            if :zlb_rule in get_setting(m, :temporary_altpolicy_names)
                if haskey(get_settings(m), :zlb_rule_value)
                    println("zlb_rule_value = ", get_setting(m, :zlb_rule_value))
                elseif haskey(get_settings(m), :zlb_rule_values)
                    println("zlb_rule_values = ", get_setting(m, :zlb_rule_values))
                else
                    println("zlb_rule_value = ", forecast_zlb_value(m))
                end
            end
        end

        println("") # add some space
    end

    return output_vars, df
end

"""
```
load_draws(m, input_type; subset_inds = 1:0, verbose = :low)

load_draws(m, input_type, block_inds; verbose = :low)
```

Load and return parameter draws from Metropolis-Hastings or SMC.

### Inputs

- `m::AbstractDSGEModel`: model object
- `input_type::Symbol`: one of the options for `input_type` described in the
  documentation for `forecast_one`
- `block_inds::AbstractRange{Int64}`: indices of the current block (already indexed by
  `jstep`) to be read in. Only used in second method

### Keyword Arguments

- `subset_inds::AbstractRange{Int64}`: indices specifying the subset of draws to be read
  in. Only used in first method
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`. If `:low` or greater, prints
  location of input file.

### Outputs

- `params`: first method returns a single parameter draw of type
  `Vector{Float64}`. Second method returns a `Vector{Vector{Float64}}` of
  parameter draws for this block.
"""
function load_draws(m::AbstractDSGEModel, input_type::Symbol;
                    subset_inds::AbstractRange{Int64} = 1:0,
                    verbose::Symbol = :low, filestring_addl::Vector{String} =
                    Vector{String}(undef, 0), use_highest_posterior_value::Bool = false,
                    input_file_name::String = "")

    if isempty(input_file_name)
        input_file_name = get_forecast_input_file(m, input_type, filestring_addl = filestring_addl)
    end

    # Load single draw
    if input_type in [:mean, :mode, :mode_draw_shocks]
        if get_setting(m, :sampling_method) == :MH
            params = convert(Vector{Float64}, h5read(input_file_name, "params"))
        elseif get_setting(m, :sampling_method) == :SMC
            if input_type in [:mode, :mode_draw_shocks]
                # If mode is in forecast overrides and is an extension, want to use the file specified
                if :mode in collect(keys(forecast_input_file_overrides(m))) &&
                    occursin(".h5", input_file_name)
                    params = convert(Vector{Float64}, h5read(input_file_name, "params"))
                    # If not, load it from the cloud
                elseif (occursin("smc_paramsmode", input_file_name) &&
                        !use_highest_posterior_value && input_type == :mode)
                    params = convert(Vector{Float64}, h5read(input_file_name, "params"))
                else
                    # Check that input_file_name is correct. Note that if
                    # it already has smc_cloud, then the following code block does not do anything.
                    if occursin("smc_paramsmode", input_file_name) || occursin(".h5", input_file_name)
                        input_file_name = replace(replace(input_file_name,
                                                          "smc_paramsmode" => "smc_cloud"),
                                                  ".h5" => ".jld2")
                        println(verbose, :low, "Switching estimation file of draws to $input_file_name")
                    end

                    cloud = load(input_file_name, "cloud")
                    params = if typeof(cloud) <: Union{DSGE.Cloud,SMC.Cloud}
                        SMC.get_highest_posterior_particle_value(SMC.Cloud(cloud))
                    else
                        cloud.particles[argmax(get_logpost(cloud))].value
                    end
                end
            else
                error("SMC mean not implemented yet")
            end
        else
            error("Invalid :sampling method specification. Change in setting :sampling_method")
        end
    # Load full distribution
    elseif input_type == :full

        if get_setting(m, :sampling_method) == :MH
            params = map(Float64, h5read(input_file_name, "mhparams"))
        elseif get_setting(m, :sampling_method) == :SMC
            input_file_name = replace(replace(input_file_name, ".h5" => ".jld2"), "smcsave" => "smc_cloud")
            cloud = load(input_file_name, "cloud")

            # Re-sample SMC draws according to their weights
            W = load(input_file_name, "W")
            weights = W[:, end]
            inds = SMC.resample(weights)

            params = if isa(cloud, ParticleCloud)
                # get_vals gets the particle values into n_parameters x n_particles,
                # we select the draws we want using inds, and finally
                # we transpose to get something with dimension n_particles x n_parameters
                Matrix{Float64}(get_vals(cloud)[:, inds]')
            else
                # Note cloud.particles has dimension n_particles x (n_parameters + 5), where the 5 is for metadata
                # when cloud is a `SMC.Cloud`
                cloud.particles[inds, 1:SMC.ind_para_end(size(cloud.particles, 2))]
            end
        else
            error("Invalid :sampling method specification. Change in setting :sampling_method")
        end

    # Load subset of full distribution
    elseif input_type == :subset

        if isempty(subset_inds)
            error("Must supply nonempty range of subset_inds if input_type == :subset")
        else
            if get_setting(m, :sampling_method) == :MH
                params = map(Float64, h5read(input_file_name, "mhparams", (subset_inds, :)))
            elseif get_setting(m, :sampling_method) == :SMC
                params = map(Float64, h5read(input_file_name, "smcparams", (subset_inds, :)))
            else
                throw("Invalid :sampling method specification. Change in setting :sampling_method")
            end
        end
    elseif input_type == :prior
        params = rand(m.parameters, n_forecast_draws(m, :prior))

    elseif input_type == :init || input_type == :init_draw_shocks
        # Return initial parameters of model object

#=        if m.spec == "het_dsge"
            init_parameters!(m, testing_gamma = false)
        else
            init_parameters!(m)
        end =#
        tmp = map(α -> α.value, m.parameters)
        params = convert(Vector{Float64}, tmp)

    end

    println(verbose, :low, "Loaded draws from $input_file_name")

    return params
end

function load_draws(m::AbstractDSGEModel, input_type::Symbol, block_inds::AbstractRange{Int64};
                    verbose::Symbol = :low,
                    filestring_addl::Vector{String} = Vector{String}(undef, 0),
                    input_file_name::String = "")

    if isempty(input_file_name)
        input_file_name = get_forecast_input_file(m, input_type, filestring_addl = filestring_addl)
    end
    println(verbose, :low, "Loading draws from $input_file_name")

    if input_type in [:full, :subset]
        if isempty(block_inds)
            error("Must supply nonempty range of block_inds for this load_draws method")
        else
            ndraws = length(block_inds)
            params = Vector{Vector{Float64}}(undef, ndraws)
            if get_setting(m, :sampling_method) == :SMC
                params_unweighted = similar(params)
            end
            for (i, j) in zip(1:ndraws, block_inds)
                if get_setting(m, :sampling_method) == :MH
                    params[i] = vec(map(Float64, h5read(input_file_name, "mhparams", (j, :))))
                elseif get_setting(m, :sampling_method) == :SMC
                    if input_file_name[end-1:end] == "h5"
                        params_unweighted[i] = vec(map(Float64, h5read(input_file_name, "smcparams", (j, :))))
                    else
                        dataset = jldopen(input_file_name, "r") do file
                            params_unweighted[i] = vec(map(Float64, file["cloud"].particles[j,1:end-5]))
                        end
                    end
                else
                    throw("Invalid :sampling_method setting specification.")
                end
            end

            # Re-sample draws according to weights if the sampling_method was SMC
            if get_setting(m, :sampling_method) == :SMC
                # Re-sample SMC draws according to their weights
                @load replace(replace(input_file_name, "smcsave" => "smc_cloud"), "h5" => "jld2") W
                weights = W[:, end][block_inds]
                inds = SMC.resample(weights)

                params = params_unweighted[inds]
            end

            return params
        end
    elseif input_type == :init_draw_shocks
        ndraws = length(block_inds)
        params = repeat([map(x -> x.value, m.parameters)], ndraws)
    elseif input_type == :mode_draw_shocks
        ndraws = length(block_inds)
        params = repeat([map(x -> x.value, m.parameters)], ndraws)
    elseif input_type == :prior
        ndraws = length(block_inds)
        params = Vector{Vector{Float64}}(undef, ndraws)
        for (i, j) in zip(1:ndraws, block_inds)
            try_again = true
            param_try = vec(rand(m.parameters, ndraws)[:, i])
            # Need to do this to ensure that draws from prior don't return -Inf likelihood
            while try_again
                param_try = vec(rand(m.parameters, ndraws)[:, i])
                try
                    update!(m, param_try)
                    likelihood(m, rand(length(m.observables), 100))
                    try_again = false
                catch
                    try_again = true
                end
            end
            params[i] = param_try
        end
        return params
    else
        error("This load_draws method can only be called with input_type in [:full, :subset]")
    end
end

function load_draws(m::AbstractDSGEVARModel, input_type::Symbol; subset_inds::AbstractRange{Int64} = 1:0,
                    verbose::Symbol = :low,
                    filestring_addl::Vector{String} = Vector{String}(undef, 0),
                    use_highest_posterior_value::Bool = false,
                    input_file_name::String = "")

    if isempty(input_file_name)
        input_file_name = get_forecast_input_file(m, input_type;
                                                  filestring_addl = filestring_addl)
    end

    return load_draws(get_dsge(m), input_type; subset_inds = subset_inds, verbose = verbose,
                      filestring_addl = filestring_addl, use_highest_posterior_value =
                      use_highest_posterior_value,
                      input_file_name = input_file_name)
end

"""
```
forecast_one(m, input_type, cond_type, output_vars; df = DataFrame(),
    subset_inds = 1:0, forecast_string = "", verbose = :low, ...)
```

Compute and save `output_vars` for input draws given by `input_type` and
conditional data case given by `cond_type`.

### Inputs

- `m::AbstractDSGEModel`: model object

- `input_type::Symbol`: one of:

```
  - `:mode`: forecast using the modal parameters only
  - `:mean`: forecast using the mean parameters only
  - `:init`: forecast using the initial parameter values only
  - `:full`: forecast using all parameters (full distribution)
  - `:subset`: forecast using a well-defined user-specified subset of draws
```

- `cond_type::Symbol`: one of:

```
  - `:none`: no conditional data
  - `:semi`: use \"semiconditional data\" - average of quarter-to-date
    observations for high frequency series
  - `:full`: use \"conditional data\" - semiconditional plus nowcasts for
    desired observables
```

- `output_vars::Vector{Symbol}`: vector of desired output variables. See
  `?forecast_one_draw`.

### Keyword Arguments

- `df::DataFrame`: Historical data. If `cond_type in [:semi, :full]`, then the
   final row of `df` should be the period containing conditional data. If not
   provided, will be loaded using `load_data` with the appropriate `cond_type`
- `subset_inds::AbstractRange{Int64}`: indices specifying the draws we want to use. If a
  more sophisticated selection criterion is desired, the user is responsible for
  determining the indices corresponding to that criterion. If `input_type` is
  not `subset`, `subset_inds` will be ignored
- `forecast_string::String`: short string identifying the subset to be
  appended to the output filenames. If `input_type = :subset` and
  `forecast_string` is empty, an error is thrown.
- `only_filter::Bool`: do not run the smoother and only run the filter. This limits the number of
  output variables which can be calculated.
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`.
- `check_empty_columns::Bool = true`: check empty columns or not when loading data (if `df` is empty)
- `bdd_fcast::Bool = true`: are we computing the bounded forecasts or not?
- `params::AbstractArray{Float64} = Vector{Float64}(undef, 0)`: parameter draws for the forecast.
     If empty, then we load draws from estimation files implied by the settings in `m`.
- `zlb_method::Symbol`: method for enforcing the zero lower bound. Defaults to `:shock`,
    meaning we use a monetary policy shock to enforce the ZLB.

  Other available methods:
  1. `:temporary_altpolicy` -> use a temporary alternative policy to enforce the ZLB.

- `rerun_smoother::Bool = false`: if true, rerun the conditional forecast when automatically enforcing
    the ZLB as a temporary alternative policy.
- `nan_endozlb_failures::Bool = false`: if true, failures of an endogenous ZLB (when implemented
    as a temporary policy) will be handled by throwing NaNs instead of using
    unanticipated monetary policy shocks.
- `set_regime_vals_altpolicy::Function`: `Function` that adds new regimes to parameters when
    using temporary alternative policies (if needed). Defaults to identity (which does nothing)
    This function should take as inputs the model object `m` and the total number of regimes
    (after adding the required temporary regimes), i.e. `set_regime_vals_altpolicy(m, n)`. It should then
    set up regime-switching parameters for these new additional regimes.
- `set_info_sets_altpolicy::Function = auto_temp_altpolicy_info_set`: `Function` that automatically updates
    the `tvis_information_set`, e.g. when `zlb_method = :temporary_altpolicy`.
- `update_regime_eqcond_info!::Function = (x1, x2, x3, x4) ->
    default_update_regime_eqcond_info(x1, x2, x3, x4, alternative_policy(m)`: `Function` that automatically
    updates the setting `:regime_eqcond_info`. The arguments of `update_regime_eqcond_info!` should be (in order)
    `m::AbstractDSGEModel`, `eqcond_dict::AbstractDict{Int64, EqcondEntry}`, `zlb_start_regime::Int64`,
    and `liftoff_regime::Int64`. The last two arguments are the regime numbers of the first regime for which
    the ZLB applies and the regime after the ZLB ends, respectively.
    The `eqcond_dict` argument should specify the `EqcondEntry` during the historical/conditional horizon regime
    (if it is desired) but can otherwise be empty. This function should then update `eqcond_dict`
    in place to implement a temporary ZLB and any other permanent alternative policies/regime-switching
    in the forecast horizon (after the conditional horizon).
    The user should also be careful and make sure `update_regime_eqcond_info!` handles imperfect awareness properly
    if they want to implement an imperfectly credible ZLB. For an example, see `?default_update_regime_eqcond_info`.
- `show_failed_percent::Bool = false`: prints out the number of failed forecasts, which are returned as NaNs.
    These may occur when the ZLB is not enforced, for example.
- `pegFFR::Bool = false`: peg the nominal FFR at the value specified by `FFRpeg`
- `FFRpeg::Float64 = -0.25/4`: value of the FFR peg
- `H::Int = 4`: number of horizons for which the FFR is pegged
- `testing_carer_kohn::Bool = false`: whether to create a file storing some property of Σ in Carter Kohn

### Outputs

None. Output is saved to files returned by
`get_forecast_output_files(m, input_type, cond_type, output_vars)`.
"""
function forecast_one(m::AbstractDSGEModel{Float64},
                      input_type::Symbol, cond_type::Symbol, output_vars::Vector{Symbol};
                      df::DataFrame = DataFrame(), subset_inds::AbstractRange{Int64} = 1:0,
                      forecast_string::String = "",
                      use_filtered_shocks_in_shockdec::Bool = false,
                      shock_name::Symbol = :none, shock_var_name::Symbol = :none,
                      shock_var_value::Float64 = 0.0, check_empty_columns = true,
                      bdd_fcast::Bool = true, params::AbstractArray{Float64} = Vector{Float64}(undef, 0),
                      zlb_method::Symbol = :shock, set_regime_vals_altpolicy::Function = identity,
                      set_info_sets_altpolicy::Function = auto_temp_altpolicy_info_set,
                      update_regime_eqcond_info!::Function =
                      (a, b, c, d) -> default_update_regime_eqcond_info!(a, b, c, d, alternative_policy(m)),
                      rerun_smoother::Bool = false, nan_endozlb_failures::Bool = false,
                      catch_smoother_lapack::Bool = false,
                      pegFFR::Bool = false, FFRpeg::Float64 = -0.25/4, H::Int = 4,
                      show_failed_percent::Bool = false, only_filter::Bool = false,
                      verbose::Symbol = :low, testing_carter_kohn::Bool = false,
                      trend_nostates_obs = Array{(0,0)}, trend_nostates_pseudo = Array{(0,0)},
                      full_shock_decomp::Bool = true, n_back::Int64 = 0, back_shocks::Vector{Symbol} = Symbol[])

    ### Common Setup

    # Add necessary output_vars and load data
    output_vars, df = prepare_forecast_inputs!(m, input_type, cond_type, output_vars;
                                               df = df, verbose = verbose, bdd_fcast = bdd_fcast,
                                               subset_inds = subset_inds,
                                               check_empty_columns = check_empty_columns, only_filter = only_filter)

    # Get output file names
    forecast_output = Dict{Symbol, Array{Float64}}()
    forecast_output_files = get_forecast_output_files(m, input_type, cond_type, output_vars;
                                                      forecast_string = forecast_string)
    output_dir = rawpath(m, "forecast")

    # Determine if we are regime_switching
    regime_switching = haskey(get_settings(m), :regime_switching) ?
        get_setting(m, :regime_switching) : false
    n_regimes        = regime_switching && haskey(get_settings(m), :n_regimes) ?
        get_setting(m, :n_regimes) : 1

    # Print
    info_print(verbose, :low, "Forecasting input_type = $input_type, cond_type = $cond_type...")
    println(verbose, :low, "Start time: $(now())")
    println(verbose, :low, "Forecast outputs will be saved in $output_dir")

    ### Single-Draw Forecasts

    if input_type in [:mode, :mean, :init]

        elapsed_time = @elapsed let
            if isempty(params)
                params = load_draws(m, input_type; verbose = verbose)
            end
            forecast_output = forecast_one_draw(m, input_type, cond_type, output_vars,
                                                params, df, verbose = verbose,
                                                shock_name = shock_name,
                                                shock_var_name = shock_var_name,
                                                shock_var_value = shock_var_value,
                                                regime_switching = regime_switching,
                                                n_regimes = n_regimes, zlb_method = zlb_method,
                                                set_regime_vals_altpolicy = set_regime_vals_altpolicy,
                                                set_info_sets_altpolicy = set_info_sets_altpolicy,
                                                update_regime_eqcond_info! = update_regime_eqcond_info!,
                                                pegFFR = pegFFR, FFRpeg = FFRpeg, H = H, only_filter = only_filter,
                                                rerun_smoother = rerun_smoother, nan_endozlb_failures = nan_endozlb_failures,
                                                catch_smoother_lapack = catch_smoother_lapack,
                                                testing_carter_kohn = testing_carter_kohn, trend_nostates_obs = trend_nostates_obs,
                                                trend_nostates_pseudo = trend_nostates_pseudo, full_shock_decomp = full_shock_decomp,
                                                n_back = n_back, back_shocks = back_shocks)

            write_forecast_outputs(m, input_type, output_vars, forecast_output_files,
                                   forecast_output; df = df, block_number = Nullable{Int64}(),
                                   verbose = verbose)
        end

        write_forecast_outputs(m, input_type, output_vars, forecast_output_files,
                               forecast_output; df = df, block_number = Nullable{Int64}(),
                               verbose = verbose)

        total_forecast_time     = elapsed_time
        total_forecast_time_min = total_forecast_time/60
        println(verbose, :low, "\nTotal time to forecast: $total_forecast_time_min minutes")

    ### Multiple-Draw Forecasts

    elseif input_type in [:full, :subset, :prior, :init_draw_shocks, :mode_draw_shocks]
        if get_setting(m, :sampling_method) == :SMC
            subset_cond = input_type == :subset && get_jstep(m, length(subset_inds)) > 1
            other_cond = input_type in [:full, :prior, :init_draw_shocks, :mode_draw_shocks] &&
                get_jstep(m, n_forecast_draws(m, :full)) > 1
            if subset_cond || other_cond
                answer = ""
                invalid_answer = true
                while invalid_answer
                    println("WARNING: The forecast draws are being thinned by $(get_jstep(m, n_forecast_draws(m,:full))). This is not recommended for forecasts using an SMC estimation. Do you want to continue? (y/n)")
                    answer = readline(stdin)
                    if (answer == "y" || answer == "n")
                        invalid_answer = false
                    else
                        println("Invalid answer! Must be `y` or `n`")
                    end
                end

                if answer == "n"
                    throw(error("Forecast halted by user."))
                end
            end
        end

        if !isempty(params) & isa(params, AbstractMatrix)
            m <= Setting(:forecast_ndraws, size(params, 1))
        elseif !isempty(params) & (input_type == :mode_draw_shocks) & (!haskey(get_settings(m), :forecast_ndraws))
            error("If using :mode_draw_shocks with a user-defined matrix of parameter draws, " *
                  "then the user must manually add the Setting :forecast_ndraws " *
                  "to indicate the number of times shocks will be drawn.")
        end

        # Block info
        block_inds, block_inds_thin = forecast_block_inds(m, input_type; subset_inds = subset_inds)
        nblocks = length(block_inds)
        start_block = forecast_start_block(m)

        # Info needed for printing progress
        total_forecast_time = 0.0
        block_verbose = verbose == :none ? :none : :low

        if testing_carter_kohn
            arred = zeros(block_inds[end][end], 246)
        end

        for block = start_block:nblocks
            println(verbose, :low, )
            info_print(verbose, :low, "Forecasting block $block of $nblocks...")
            begin_time = time_ns()

            # Get to work!
            if isempty(params)
                params_for_map = load_draws(m, input_type, block_inds[block]; verbose = verbose)
            elseif input_type == :mode_draw_shocks
                ndraws = length(block_inds[block])
                @assert isa(params, Vector) "To use mode_draw_shocks with params passed in as a keyword, params must be a Vector."
                params_for_map = Vector{Float64}[params for i in block_inds[block]]
            else
                params_for_map = Vector{Float64}[params[i, :] for i in block_inds[block]]
            end

            mapfcn = use_parallel_workers(m) ? pmap : map
            forecast_outputs = mapfcn(param -> forecast_one_draw(m, input_type, cond_type, output_vars,
                                                                 param, df, verbose = verbose,
                                                                 use_filtered_shocks_in_shockdec =
                                                                 use_filtered_shocks_in_shockdec,
                                                                 shock_name = shock_name,
                                                                 shock_var_name = shock_var_name,
                                                                 shock_var_value = shock_var_value,
                                                                 regime_switching = regime_switching,
                                                                 n_regimes = n_regimes, zlb_method = zlb_method,
                                                                 set_regime_vals_altpolicy = set_regime_vals_altpolicy,
                                                                 set_info_sets_altpolicy = set_info_sets_altpolicy,
                                                                 update_regime_eqcond_info! = update_regime_eqcond_info!,
                                                                 pegFFR = pegFFR, FFRpeg = FFRpeg, H = H, only_filter = only_filter,
                                                                 rerun_smoother = rerun_smoother,
                                                                 nan_endozlb_failures = nan_endozlb_failures,
                                                                 catch_smoother_lapack = catch_smoother_lapack,
                                                                 testing_carter_kohn = testing_carter_kohn,
                                                                 n_back = n_back, back_shocks = back_shocks),
                                      params_for_map)

            # Assemble outputs from this block and write to file
            forecast_outputs = convert(Vector{Dict{Symbol, Array{Float64}}}, forecast_outputs)

            if testing_carter_kohn && get_setting(m, :forecast_smoother) == :carter_kohn
                for i in 1:length(forecast_outputs)
                    arred[block_inds[block][i],1:length(forecast_outputs[i][:conded])] = forecast_outputs[i][:conded]
                    delete!(forecast_outputs[i], :conded)
                end
            end

            forecast_output = assemble_block_outputs(forecast_outputs; show_failed_percent = show_failed_percent)
            write_forecast_outputs(m, input_type, output_vars, forecast_output_files,
                                   forecast_output; df = df, block_number = Nullable(block),
                                   verbose = block_verbose, block_inds = block_inds_thin[block],
                                   subset_inds = subset_inds)
            GC.gc()

            # Calculate time to complete this block, average block time, and
            # expected time to completion
            block_time = (time_ns() - begin_time)/1e9
            total_forecast_time += block_time
            total_forecast_time_min     = total_forecast_time/60
            blocks_elapsed              = block - start_block + 1
            expected_time_remaining     = (total_forecast_time/blocks_elapsed)*(nblocks - block)
            expected_time_remaining_min = expected_time_remaining/60

            println(verbose, :low, "\nCompleted $block of $nblocks blocks.")
            println(verbose, :low, "Total time elapsed: $total_forecast_time_min minutes")
            println(verbose, :low, "Expected time remaining: $expected_time_remaining_min minutes")
        end # of loop through blocks

        if testing_carter_kohn && get_setting(m, :forecast_smoother) == :carter_kohn
            jldopen("conded_svd_diagm.jld2", true, true, true, IOStream) do file
                write(file, "conds", arred)
            end
        end
    end # of input_type
    combine_raw_forecast_output_and_metadata(m, forecast_output_files, verbose = verbose)

    println(verbose, :low, "\nForecast complete: $(now())")
end

"""
```
forecast_one_draw(m, input_type, cond_type, output_vars, params, df;
    verbose = :low, only_filter = false)
```

Compute `output_vars` for a single parameter draw, `params`. Called by
`forecast_one`.

### Inputs

- `m::AbstractDSGEModel{Float64}`: model object
- `input_type::Symbol`: See `?forecast_one`.
- `cond_type::Symbol`: See `?forecast_one`.
- `output_vars::Vector{Symbol}`: vector of desired output variables. See Outputs
  section
- `params::Vector{Float64}`: parameter vector
- `df::DataFrame`: historical data.
- `verbose::Symbol`: desired frequency of function progress messages printed to
  standard out. One of `:none`, `:low`, or `:high`.

### Output

- `forecast_outputs::Dict{Symbol, Array{Float64}}`: dictionary of forecast
  outputs. Keys are `output_vars`, which is some subset of:

```
  - `:histstates`: `Matrix{Float64}` of smoothed historical states
  - `:histobs`: `Matrix{Float64}` of smoothed historical data
  - `:histpseudo`: `Matrix{Float64}` of smoothed historical
    pseudo-observables
  - `:histshocks`: `Matrix{Float64}` of smoothed historical shocks
  - `:forecaststates`: `Matrix{Float64}` of forecasted states
  - `:forecastobs`: `Matrix{Float64}` of forecasted observables
  - `:forecastpseudo`: `Matrix{Float64}` of forecasted pseudo-observables
  - `:forecastshocks`: `Matrix{Float64}` of forecasted shocks
  - `:bddforecaststates`, `:bddforecastobs`, `:bddforecastpseudo`, and
    `:bddforecastshocks`: `Matrix{Float64}`s of forecasts where we enforce
    the zero lower bound to be `forecast_zlb_value(m)`
  - `:shockdecstates`: `Array{Float64, 3}` of state shock decompositions
  - `:shockdecobs`: `Array{Float64, 3}` of observable shock decompositions
  - `:shockdecpseudo`: `Array{Float64, 3}` of pseudo-observable shock
    decompositions
  - `:dettrendstates`: `Matrix{Float64}` of state deterministic trends
  - `:dettrendobs`: `Matrix{Float64}` of observable deterministic trends
  - `:dettrendpseudo`: `Matrix{Float64}` of pseudo-observable deterministic
    trends
  - `:trendstates`: `Vector{Float64}` of state trends, i.e. the `CCC` vector
  - `:trendobs`: `Vector{Float64}` of observable trends, i.e. the `DD` vector
  - `:trendpseudo`: `Vector{Float64}` of pseudo-observable trends, i.e. the
    `DD_pseudo` vector
  - `:irfstates`: `Array{Float64, 3}` of state impulse responses
  - `:irfobs`: `Array{Float64, 3}` of observable impulse responses
  - `:irfpseudo`: `Array{Float64, 3}` of pseudo-observable impulse responses
```
"""
function forecast_one_draw(m::AbstractDSGEModel{Float64}, input_type::Symbol, cond_type::Symbol,
                           output_vars::Vector{Symbol}, params::Vector{Float64}, df::DataFrame; verbose::Symbol = :low,
                           use_filtered_shocks_in_shockdec::Bool = false,
                           shock_name::Symbol = :none, shock_var_name::Symbol = :none,
                           shock_var_value::Float64 = 0.0, zlb_method::Symbol = :shock,
                           set_regime_vals_altpolicy::Function = identity,
                           set_info_sets_altpolicy::Function = auto_temp_altpolicy_info_set,
                           update_regime_eqcond_info!::Function =
                           (a, b, c, d) -> default_update_regime_eqcond_info!(a, b, c, d, alternative_policy(m)),
                           pegFFR::Bool = false, FFRpeg::Float64 = -0.25/4, H::Int = 4,
                           regime_switching::Bool = false, n_regimes::Int = 1, only_filter::Bool = false,
                           filter_smooth::Bool = false, rerun_smoother::Bool = false,
                           nan_endozlb_failures::Bool = false,
                           catch_smoother_lapack::Bool = false, testing_carter_kohn::Bool = false,
                           return_loglh::Bool = false, trend_nostates_obs = Array{(0,0)},
                           trend_nostates_pseudo = Array{(0,0)}, full_shock_decomp::Bool = false,
                           n_back::Int64 = 0, back_shocks::Vector{Symbol} = Symbol[])
    ### Setup

    # Re-initialize model indices if forecasting under an alternative policy
    # rule (in case new states or equations were added)
    if alternative_policy(m).key != :historical
        init_model_indices!(m)
    end

    # Time-Varying Information Set?
    tvis = haskey(get_settings(m), :tvis_information_set) ? !isempty(get_setting(m, :tvis_information_set)) : false

    # Are we only running IRFs?
    output_prods = map(get_product, output_vars)
    irfs_only = all(x -> x == :irf, output_prods)

    # Compute state space
    update!(m, params) # Note that params is a Vector{Float64}, not a ParameterVector. This `update!` infers if the forecast is regime-switching if length(params) > length(m.parameters)

    system = compute_system(m; tvis = tvis)
    if haskey(m.settings, :old_shock_decs) && get_setting(m, :old_shock_decs)
        # Must be using TV cred system
        m2 = deepcopy(m)
        m2 <= Setting(:regime_eqcond_info, get_setting(m2, :alternative_policies)[1].regime_eqcond_info)
        m2 <= Setting(:tvis_information_set, get_setting(m2, :alternative_policies)[1].infoset)
        m2 <= Setting(:identical_eqcond_regimes, get_setting(m2, :alternative_policies)[1].identical_eqcond_regimes)
        m2 <= Setting(:perfect_credibility_identical_transitions, get_setting(m2, :alternative_policies)[1].perfect_credibility_identical_transitions)
        m2 <= Setting(:n_regimes, get_setting(m2, :alternative_policies)[1].n_regimes)
        m2 <= Setting(:gensys2, get_setting(m2, :alternative_policies)[1].gensys2)
        delete!(m2.settings, :perfect_credibility_identical_transitions)
        for i in collect(keys(get_setting(m2, :regime_eqcond_info)))
            get_setting(m2, :regime_eqcond_info)[i].weights = [1.0, 0.0]
        end

        old_system = compute_system(m2; tvis=tvis)
    end

    # Initialize output dictionary
    forecast_output = Dict{Symbol, Array{Float64}}()

    # Decide whether to draw states/shocks in smoother/forecast
    uncertainty_override = forecast_uncertainty_override(m)
    uncertainty = if Nullables.isnull(uncertainty_override)
        if input_type in [:init, :mode, :mean, :prior]
            false
        elseif input_type in [:full, :subset, :mode_draw_shocks, :init_draw_shocks]
            true
        end
    else
        get(uncertainty_override)
    end

    ### 1. Smoothed Histories

    # Must run smoother for conditional data in addition to explicit cases
    hist_vars = [:histstates, :histpseudo, :histshocks, :histstdshocks]
    shockdec_vars = [:shockdecstates, :shockdecpseudo, :shockdecobs]
    shockdecseq_vars = [:shockdecseqstates, :shockdecseqpseudo, :shockdecseqobs]
    dettrend_vars = [:dettrendstates, :dettrendpseudo, :dettrendobs]
    smooth_vars = vcat(hist_vars, shockdec_vars, dettrend_vars, shockdecseq_vars)
    hists_to_compute = intersect(output_vars, smooth_vars)

    run_smoother = (!isempty(hists_to_compute) ||
        (cond_type in [:semi, :full] && !irfs_only)) && !only_filter
    if run_smoother
        # Call smoother
        if filter_smooth && (get_setting(m, :forecast_smoother) == :carter_kohn)
            kal = filter(m, df, system; cond_type = cond_type)
            histstates, histshocks, histpseudo, initial_states =
                smooth(m, df, system; cond_type = cond_type, draw_states = uncertainty,
                       s_pred = kal[:s_pred], P_pred = kal[:P_pred], s_filt = kal[:s_filt], P_filt = kal[:P_filt],
                       catch_smoother_lapack = catch_smoother_lapack)

            if haskey(m.settings, :old_shock_decs) && get_setting(m, :old_shock_decs)
                old_kal = filter(m2, df, old_system; cond_type = cond_type)
                old_histstates, old_histshocks, old_histpseudo, old_initial_states =
                    smooth(m2, df, old_system; cond_type = cond_type, draw_states = uncertainty,
                           s_pred = kal[:s_pred], P_pred = kal[:P_pred], s_filt = kal[:s_filt], P_filt = kal[:P_filt],
                           catch_smoother_lapack = catch_smoother_lapack)
            end
        else
            histstates, histshocks, histpseudo, initial_states =
                smooth(m, df, system; cond_type = cond_type, draw_states = uncertainty,
                       catch_smoother_lapack = catch_smoother_lapack)

            if haskey(m.settings, :old_shock_decs) && get_setting(m, :old_shock_decs)
                old_histstates, old_histshocks, old_histpseudo, old_initial_states =
                smooth(m2, df, old_system; cond_type = cond_type, draw_states = uncertainty,
                       catch_smoother_lapack = catch_smoother_lapack)
            end
        end

        # For conditional data, transplant the obs/state/pseudo vectors from hist to forecast
        if cond_type in [:full, :semi]
            T = n_mainsample_periods(m)

            forecast_output[:histstates] = transplant_history(histstates, T)
            forecast_output[:histshocks] = transplant_history(histshocks, T)
            forecast_output[:histpseudo] = transplant_history(histpseudo, T)
        else
            forecast_output[:histstates] = histstates
            forecast_output[:histshocks] = histshocks
            forecast_output[:histpseudo] = histpseudo
        end

        # Standardize shocks if desired
        if :histstdshocks in output_vars
            if regime_switching
                start_date = max(date_mainsample_start(m), df[1, :date]) # use mainsample b/c shocks only includes mainsample, no presample
                end_date   = prev_quarter(date_forecast_start(m))
                regime_inds = regime_indices(m, start_date, end_date)
                if regime_inds[1][1] < 1
                    regime_inds[1] = 1:regime_inds[1][end]
                end
                forecast_output[:histstdshocks] =
                    standardize_shocks(forecast_output[:histshocks],
                                       Matrix{eltype(system[1, :QQ])}[system[i, :QQ] for i in 1:n_regimes], regime_inds)
            else
                forecast_output[:histstdshocks] = standardize_shocks(forecast_output[:histshocks], system[:QQ])
            end
        end
    end

    ### 2. Forecasts

    unbddforecast_vars = [:forecaststates, :forecastobs, :forecastpseudo, :forecastshocks, :forecaststdshocks]
    bddforecast_vars = [:bddforecaststates, :bddforecastobs, :bddforecastpseudo, :bddforecastshocks, :bddforecaststdshocks]
    forecast_vars = vcat(unbddforecast_vars, bddforecast_vars)
    forecasts_to_compute = intersect(output_vars, forecast_vars)

    if !isempty(forecasts_to_compute)
        # Get initial forecast state vector s_T
        s_T = if run_smoother && !return_loglh # ONLY THIS BRANCH WORKS FOR REGIME SWITCHING
            # The last smoothed state is either s_{T|T} (if !uncertainty) or
            # drawn from N(s_{T|T}, P_{T|T}) (if uncertainty)
            histstates[:, end]
        else
            kal = Kalman(Vector{Float64}(undef,0), Matrix{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0, 0), Matrix{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0, 0), Vector{Float64}(undef, 0), Array{Float64}(undef, 0, 0, 0), Vector{Float64}(undef, 0), Array{Float64}(undef, 0, 0, 0))
            try
                kal = filter(m, df, system; cond_type = cond_type)
            catch err
                if isa(err, DomainError)
                    return Dict{Symbol, Array{Float64}}()
                else
                    rethrow(err)
                end
            end
            if uncertainty
                # If we want to draw s_T but haven't run the smoother, draw from
                # N(s_{T|T}, P_{T|T}) directly
                U, singular_values, _ = svd(kal[:P_T])
                dist = DegenerateMvNormal(kal[:s_T], U*Matrix(Diagonal(sqrt.(singular_values))))
                rand(dist)
            else
                # If we don't want to draw s_T, simply use the mean s_{T|T}
                kal[:s_T]
            end
        end

        # Re-solve model with alternative policy rule, if applicable
        apply_altpolicy = alternative_policy(m).solve != solve
        if apply_altpolicy

            if haskey(get_settings(m), :skip_altpolicy_state_init) ? !get_setting(m, :skip_altpolicy_state_init) : true
                # Adjust the initial state vector for pgap and ygap
                if haskey(m.settings, :pgap_type) && haskey(get_settings(m), :pgap_value)
                    update_hist = true
                    _pgap_type = get_setting(m, :pgap_type)
                    if _pgap_type == :ngdp
                        _, s_T = ngdp_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
                    elseif _pgap_type == :ait
                        _, s_T = ait_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
                    else
                        update_hist = false
                    end
                    if update_hist # Don't do this computation unless pgap_type in [:ngdp, :ait]
                        histstates[:, end] = s_T
                        ZZ_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :ZZ_pseudo] : system[:ZZ_pseudo] # currently assumes these are
                        DD_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :DD_pseudo] : system[:DD_pseudo] # the same across regimes
                        histpseudo[:, end] = ZZ_pseudo * s_T + DD_pseudo
                    end
                end

                if haskey(m.settings, :ygap_type) && haskey(get_settings(m), :ygap_value) &&
                    haskey(m.settings, :pgap_type) && haskey(get_settings(m), :pgap_value)
                    update_hist = true
                    _ygap_type = get_setting(m, :ygap_type)
                    _pgap_type = get_setting(m, :pgap_type)
                    if _ygap_type == _pgap_type == :smooth_ait_gdp
                        _, s_T = smooth_ait_gdp_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
                    elseif _ygap_type == _pgap_type == :smooth_ait_gdp_alt
                        _, s_T = smooth_ait_gdp_alt_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
                    elseif (_ygap_type == _pgap_type == :flexible_ait) &&
                        (haskey(m.settings, :flexible_ait_policy_change) ?
                         !get_setting(m, :flexible_ait_policy_change) : true) &&
                         (haskey(m.settings, :initialize_pgap_ygap) ? get_setting(m, :initialize_pgap_ygap) : false)
                        _, s_T = flexible_ait_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
                    else
                        update_hist = false
                    end
                    if update_hist # Don't do this computation unless pgap_type and ygap_type are okay
                        histstates[:, end] = s_T
                        ZZ_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :ZZ_pseudo] : system[:ZZ_pseudo] # currently assumes these are
                        DD_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :DD_pseudo] : system[:DD_pseudo] # the same across regimes
                        histpseudo[:, end] = ZZ_pseudo * s_T + DD_pseudo
                    end
                end

                if haskey(m.settings, :Rref_type)
                    if get_setting(m, :ygap_type) == :rw && get_setting(m, :pgap_type) == :rw &&
                        get_setting(m, :Rref_type) == :rw
                        _, s_T = rw_forecast_init(m, zeros(0, 0), s_T, cond_type = cond_type)
                        histstates[:, end] = s_T
                        ZZ_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :ZZ_pseudo] : system[:ZZ_pseudo] # currently assumes these are
                        DD_pseudo = isa(system, RegimeSwitchingSystem) ? system[1, :DD_pseudo] : system[:DD_pseudo] # the same across regimes
                        histpseudo[:, end] = ZZ_pseudo * s_T + DD_pseudo
                    end
                end
            end
        end

        # 2A. Unbounded forecasts
        if !isempty(intersect(output_vars, unbddforecast_vars))
            if pegFFR
                nshocks = size(system[:RRR], 2)
                nstates = size(system[:TTT], 1)
                shocks = m.exogenous_shocks
                PsiR1 = 0
                PsiR2 = zeros(nstates)
                PsiR2[m.endogenous_states[:R_t]] = 1.
                Rht = system[:RRR][:, vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")])]
                bb = zeros(H+1, 1)
                MH = zeros(H+1, H+1)
                for hh = 1:H+1
                    bb[hh, 1] = (FFRpeg - PsiR1 - PsiR2'*(system[:TTT])^hh*s_T)[1]
                    MH[hh, :] = PsiR2'*(system[:TTT])^(hh-1)*Rht
                end
                monshocks = MH\bb
                etpeg = zeros(nshocks, forecast_horizons(m))
                etpeg[vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")]), 1] = monshocks
                forecaststates, forecastobs, forecastpseudo, forecastshocks =
                    forecast(system, s_T, etpeg; cond_type = cond_type, enforce_zlb = false, draw_shocks = uncertainty)
                println("The forecasted interest rate path is $(forecastobs[m.observables[:obs_nominalrate], :])")
            else
                forecaststates, forecastobs, forecastpseudo, forecastshocks =
                    forecast(m, system, s_T;
                             cond_type = cond_type, enforce_zlb = false, draw_shocks = uncertainty)
            end

            # For conditional data, transplant the obs/state/pseudo vectors from hist to forecast
            # NOTE: ZZ REGIME SWITCHING NOT FULLY SUPPORTED, SO JUST TAKE THE LAST ZZ IN THE SYSTEM
            #       OR THE ZZ SPECIFIED BY THE FOLLOWING SETTING
            cond_ZZ_regime = if regime_switching
                if haskey(get_settings(m), :cond_ZZ_regime)
                    get_setting(m, :cond_ZZ_regime)
                else # an educated guess at the correct conditional regime
                    if only_filter
                        @warn "The setting :cond_ZZ_regime is not specified. Making an educated guess..."
                    end
                    if get_setting(m, :n_cond_regimes) > 0 # if regime-switching during conditional horizon, need to add them and
                        get_setting(m, :reg_forecast_start) + get_setting(m, :n_cond_regimes) - 1 # minus 1 for reg_forecast_start
                    else # otherwise no regime-switching during conditional horizon, so just use the same regime as forecast start date
                        get_setting(m, :reg_forecast_start)
                    end
                end
            else
                1
            end
            if cond_type in [:full, :semi] && !only_filter
                forecast_output[:forecaststates] = transplant_forecast(histstates, forecaststates, T)
                forecast_output[:forecastshocks] = transplant_forecast(histshocks, forecastshocks, T)
                forecast_output[:forecastpseudo] = transplant_forecast(histpseudo, forecastpseudo, T)
                forecast_output[:forecastobs]    =
                    transplant_forecast_observables(histstates, forecastobs,
                                                    regime_switching ? system[cond_ZZ_regime] : system, T)
            elseif cond_type in [:full, :semi] && only_filter
                if regime_switching
                    cond_obs                         = system[cond_ZZ_regime, :ZZ] * s_T + system[cond_ZZ_regime, :DD]
                    cond_pseudo                      = system[cond_ZZ_regime, :ZZ_pseudo] * s_T + system[cond_ZZ_regime, :DD_pseudo]
                    forecast_output[:forecaststates] = hcat(s_T,         forecaststates)
                    forecast_output[:forecastobs]    = hcat(cond_obs,    forecastobs)
                    forecast_output[:forecastpseudo] = hcat(cond_pseudo, forecastpseudo)
                else
                    cond_obs                         = system[:ZZ] * s_T + system[:DD]
                    cond_pseudo                      = system[:ZZ_pseudo] * s_T + system[:DD_pseudo]
                    forecast_output[:forecaststates] = hcat(s_T,         forecaststates)
                    forecast_output[:forecastobs]    = hcat(cond_obs,    forecastobs)
                    forecast_output[:forecastpseudo] = hcat(cond_pseudo, forecastpseudo)
                    # forecastshocks cannot be calculated w/out smoother
                end
            else
                forecast_output[:forecaststates] = forecaststates
                forecast_output[:forecastshocks] = forecastshocks
                forecast_output[:forecastpseudo] = forecastpseudo
                forecast_output[:forecastobs]    = forecastobs
            end

            # Standardize shocks if desired
            if :forecaststdshocks in output_vars
                forecast_output[:forecaststdshocks] =
                    standardize_shocks(forecast_output[:forecastshocks],
                                       regime_switching ? system[cond_ZZ_regime, :QQ] : system[:QQ])
            end
        end

        # 2B. Bounded forecasts

        if !isempty(intersect(output_vars, bddforecast_vars))
            if pegFFR
                nshocks = size(system[:RRR], 2)
                nstates = size(system[:TTT], 1)
                shocks = m.exogenous_shocks
                PsiR1 = 0
                PsiR2 = zeros(nstates)
                PsiR2[m.endogenous_states[:R_t]] = 1.
                Rht = system[:RRR][:, vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")])]
                bb = zeros(H+1, 1)
                MH = zeros(H+1, H+1)
                @show s_T[m.endogenous_states[:R_t]] + 400*log(m[:Rstarn])
                for hh = 1:H+1
                    bb[hh, 1] = (FFRpeg - PsiR1 - PsiR2'*(system[:TTT])^hh*s_T)[1]
                    MH[hh, :] = PsiR2'*(system[:TTT])^(hh-1)*Rht
                end
                monshocks = MH\bb
                etpeg = zeros(nshocks, forecast_horizons(m))
                etpeg[vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")]), 1] = monshocks
                forecaststates, forecastobs, forecastpseudo, forecastshocks = forecast(system, s_T, etpeg)
                @show forecaststates[m.endogenous_states[:R_t], :]
                @show forecastobs[m.observables[:obs_nominalrate], :]
            elseif zlb_method == :temporary_altpolicy
                # Run the unbounded forecast if they haven't already been computed
                if isempty(intersect(output_vars, unbddforecast_vars))
                    forecaststates, forecastobs, forecastpseudo, forecastshocks =
                        forecast(m, system, s_T;
                                 cond_type = cond_type, enforce_zlb = false, draw_shocks = uncertainty)
                end

                # Now run the ZLB enforcing forecast
                if rerun_smoother
                    zlb_enforced_output =
                        forecast(m, s_T, forecaststates, forecastobs, forecastpseudo, forecastshocks;
                                 cond_type = cond_type,
                                 set_zlb_regime_vals = set_regime_vals_altpolicy,
                                 set_info_sets_altpolicy = set_info_sets_altpolicy,
                                 update_regime_eqcond_info! = update_regime_eqcond_info!,
                                 rerun_smoother = rerun_smoother, df = df,
                                 nan_failures = nan_endozlb_failures,
                                 draw_states = uncertainty,
                                 histstates = histstates, histshocks = histshocks, histpseudo = histpseudo,
                                 initial_states = initial_states)
                    forecaststates, forecastobs, forecastpseudo, histstates, histshocks, histpseudo, initial_states =
                        zlb_enforced_output

                    if cond_type in [:full, :semi]
                        T = n_mainsample_periods(m)

                        forecast_output[:histstates] = transplant_history(histstates, T)
                        forecast_output[:histshocks] = transplant_history(histshocks, T)
                        forecast_output[:histpseudo] = transplant_history(histpseudo, T)
                    else
                        forecast_output[:histstates] = histstates
                        forecast_output[:histshocks] = histshocks
                        forecast_output[:histpseudo] = histpseudo
                    end

                    # Standardize shocks if desired
                    if :histstdshocks in output_vars
                        if regime_switching
                            start_date = max(date_mainsample_start(m), df[1, :date]) # use mainsample b/c shocks only includes mainsample, no presample
                            end_date   = prev_quarter(date_forecast_start(m))
                            regime_inds = regime_indices(m, start_date, end_date)
                            if regime_inds[1][1] < 1
                                regime_inds[1] = 1:regime_inds[1][end]
                            end
                            forecast_output[:histstdshocks] =
                                standardize_shocks(forecast_output[:histshocks],
                                                   Matrix{eltype(system[1, :QQ])}[system[i, :QQ] for i in 1:n_regimes], regime_inds)
                        else
                            forecast_output[:histstdshocks] = standardize_shocks(forecast_output[:histshocks], system[:QQ])
                        end
                    end
                else
                    forecaststates, forecastobs, forecastpseudo =
                        forecast(m, s_T, forecaststates, forecastobs, forecastpseudo, forecastshocks;
                                 cond_type = cond_type,
                                 set_zlb_regime_vals = set_regime_vals_altpolicy,
                                 set_info_sets_altpolicy = set_info_sets_altpolicy,
                                 update_regime_eqcond_info! = update_regime_eqcond_info!)
                end
            else
                forecaststates, forecastobs, forecastpseudo, forecastshocks =
                    forecast(m, system, s_T;
                             cond_type = cond_type, enforce_zlb = true, draw_shocks = uncertainty)
            end

            # For conditional data, transplant the obs/state/pseudo vectors from hist to forecast
            # NOTE: ZZ REGIME SWITCHING NOT FULLY SUPPORTED, SO JUST TAKE THE LAST ZZ IN THE SYSTEM
            #       OR THE ZZ SPECIFIED BY THE FOLLOWING SETTING
            cond_ZZ_regime = if regime_switching
                if haskey(get_settings(m), :cond_ZZ_regime)
                    get_setting(m, :cond_ZZ_regime)
                else # an educated guess at the correct conditional regime
                    if only_filter
                        @warn "The setting :cond_ZZ_regime is not specified. Making an educated guess..."
                    end
                    if get_setting(m, :n_cond_regimes) > 0 # if regime-switching during conditional horizon, need to add them and
                        get_setting(m, :reg_forecast_start) + get_setting(m, :n_cond_regimes) - 1 # minus 1 for reg_forecast_start
                    else # otherwise no regime-switching during conditional horizon, so just use the same regime as forecast start date
                        get_setting(m, :reg_forecast_start)
                    end
                end
            else
                1
            end
            if cond_type in [:full, :semi] && !only_filter
                forecast_output[:bddforecaststates] = transplant_forecast(histstates, forecaststates, T)
                forecast_output[:bddforecastshocks] = transplant_forecast(histshocks, forecastshocks, T)
                forecast_output[:bddforecastpseudo] = transplant_forecast(histpseudo, forecastpseudo, T)
                forecast_output[:bddforecastobs]    =
                    transplant_forecast_observables(histstates, forecastobs,
                                                    regime_switching ? system[cond_ZZ_regime] : system, T)
            elseif cond_type in [:full, :semi] && only_filter
                forecast_output[:bddforecaststates] = hcat(s_T,forecaststates)
                #forecast_output[:bddforecastshocks] = hcat(s_T,forecastshocks)
                #forecast_output[:bddforecastpseudo] = hcat(s_T,forecastpseudo)
                cond_var = system[cond_ZZ_regime][:ZZ]*s_T .+ system[cond_ZZ_regime][:DD]
                forecast_output[:bddforecastobs]    = hcat(cond_var,forecastobs)
            else
                forecast_output[:bddforecaststates] = forecaststates
                forecast_output[:bddforecastshocks] = forecastshocks
                forecast_output[:bddforecastpseudo] = forecastpseudo
                forecast_output[:bddforecastobs]    = forecastobs
            end

            # Standardize shocks if desired
            if :bddforecaststdshocks in output_vars
                forecast_output[:bddforecaststdshocks] =
                    standardize_shocks(forecast_output[:bddforecastshocks],
                                       regime_switching ? system[cond_ZZ_regime, :QQ] : system[:QQ])
            end
        end
    end

    if return_loglh
        return kal[:loglh]
    end

    ### 3. Shock Decompositions

    shockdecs_to_compute = intersect(output_vars, shockdec_vars)

    if !isempty(shockdecs_to_compute)
        histshocks_shockdec = if use_filtered_shocks_in_shockdec
            filter_shocks(m, df, system, cond_type = cond_type)
        else
            histshocks
        end

        if haskey(m.settings, :old_shock_decs) && get_setting(m, :old_shock_decs)
            old_histshocks_shockdec = if use_filtered_shocks_in_shockdec
                filter_shocks(m, df, old_system, cond_type = cond_type)
            else
                old_histshocks
            end
        end

        start_date = max(date_mainsample_start(m), df[1, :date]) # smooth doesn't return presample
        end_date   = if cond_type in [:semi, :full] # end date of histshocks includes conditional periods
            max(date_conditional_end(m), prev_quarter(date_forecast_start(m)))
        else
            prev_quarter(date_forecast_start(m)) # this is the end date of history period
        end

        if haskey(m.settings, :old_shock_decs) && get_setting(m, :old_shock_decs)
            shockdecstates, shockdecobs, shockdecpseudo = shock_decompositions(m, system, old_system, histshocks_shockdec, old_histshocks_shockdec, start_date, end_date, cond_type, full_shock_decomp = full_shock_decomp)
        else
            shockdecstates, shockdecobs, shockdecpseudo = isa(system, RegimeSwitchingSystem) ?
                shock_decompositions(m, system, histshocks_shockdec, start_date, end_date, cond_type) :
                shock_decompositions(m, system, histshocks_shockdec)
        end

        forecast_output[:shockdecstates] = shockdecstates
        forecast_output[:shockdecobs]    = shockdecobs
        forecast_output[:shockdecpseudo] = shockdecpseudo
    end

    ### 4. Trend

    trend_vars = [:trendstates, :trendobs, :trendpseudo]
    trends_to_compute = intersect(output_vars, trend_vars)

    if !isempty(trends_to_compute)
        trendstates, trendobs, trendpseudo = if regime_switching &&
            (haskey(get_settings(m), :time_varying_trends) ? get_setting(m, :time_varying_trends) : false)

            start_date = max(date_mainsample_start(m), df[1, :date]) # smooth doesn't return presample
            end_date   = if cond_type in [:semi, :full] # end date of histshocks includes conditional periods
                max(date_conditional_end(m), prev_quarter(date_forecast_start(m)), date_forecast_end(m))
            else
                max(prev_quarter(date_forecast_start(m)), date_forecast_end(m)) # this is the end date of history period
            end # need to add date_forecast_end to end_date to avoid dimension mismatch in write_forecast_output

            trends(m, system, start_date, end_date, cond_type)
        else
            trends(system)
        end

        forecast_output[:trendstates] = trendstates
        forecast_output[:trendobs]    = trendobs
        forecast_output[:trendpseudo] = trendpseudo

        if trend_nostates_obs != Array{(0,0)}
            forecast_output[:trendobs] = trend_nostates_obs
        end
        if trend_nostates_pseudo != Array{(0,0)}
            forecast_output[:trendpseudo] = trend_nostates_pseudo
        end
    end

    ### 5. Deterministic Trend

    dettrends_to_compute = intersect(output_vars, dettrend_vars)

    if !isempty(dettrends_to_compute)
        if isa(system, RegimeSwitchingSystem)
            start_date = max(date_mainsample_start(m), df[1, :date]) # smooth doesn't return presample
            end_date   = if cond_type in [:semi, :full] # end date of histshocks includes conditional periods
                max(date_conditional_end(m), prev_quarter(date_forecast_start(m)))
            else
                prev_quarter(date_forecast_start(m)) # this is the end date of history period
            end

            if haskey(m.settings, :old_shock_decs) && get_setting(m, :old_shock_decs) && full_shock_decomp
                dettrendstates, dettrendobs, dettrendpseudo =
                    deterministic_trends(m, system, old_system, initial_states, start_date, end_date, cond_type)
            else
                dettrendstates, dettrendobs, dettrendpseudo =
                    deterministic_trends(m, system, initial_states, start_date, end_date, cond_type)
            end
        else
            dettrendstates, dettrendobs, dettrendpseudo = deterministic_trends(m, system, initial_states)
        end

        forecast_output[:dettrendstates] = dettrendstates
        forecast_output[:dettrendobs]    = dettrendobs
        forecast_output[:dettrendpseudo] = dettrendpseudo
    end

    ### 6. Impulse Responses

    irf_vars = [:irfstates, :irfobs, :irfpseudo]
    irfs_to_compute = intersect(output_vars, irf_vars)

    if !isempty(irfs_to_compute)
        if shock_name!=:none
            irfstates, irfobs, irfpseudo = impulse_responses(m, system, impulse_response_horizons(m),
                                                             shock_name,
                                                             shock_var_name,
                                                             shock_var_value)
        else
            irfstates, irfobs, irfpseudo = impulse_responses(m, system)
        end
        forecast_output[:irfstates] = irfstates
        forecast_output[:irfobs] = irfobs
        forecast_output[:irfpseudo] = irfpseudo
    end

    # 7. Shockdec Sequence
    shockdecseqs_to_compute = intersect(output_vars, shockdecseq_vars)

    if !isempty(shockdecseqs_to_compute)
        histshocks_shockdecseq = if use_filtered_shocks_in_shockdec
            filter_shocks(m, df, system, cond_type = cond_type)
        else
            histshocks
        end
        if haskey(m.settings, :old_shock_decs) && get_setting(m, :old_shock_decs)
            old_histshocks_shockdecseq = if use_filtered_shocks_in_shockdec
                filter_shocks(m2, df, old_system, cond_type = cond_type)
            else
                old_histshocks
            end
        end

        start_date = max(date_mainsample_start(m), df[1, :date]) # smooth doesn't return presample
        end_date   = if cond_type in [:semi, :full] # end date of histshocks includes conditional periods
            max(date_conditional_end(m), prev_quarter(date_forecast_start(m)))
        else
            prev_quarter(date_forecast_start(m)) # this is the end date of history period
        end

        if haskey(m.settings, :old_shock_decs) && get_setting(m, :old_shock_decs)
            shockdecseqstates, shockdecseqobs, shockdecseqpseudo = shock_decompositions_sequence(m, system, old_system, histshocks_shockdecseq, old_histshocks_shockdecseq, start_date, end_date, cond_type, full_shock_decomp = full_shock_decomp, n_back = n_back, back_shocks = back_shocks)
        else
            shockdecseqstates, shockdecseqobs, shockdecseqpseudo = isa(system, RegimeSwitchingSystem) ?
                shock_decompositions_sequence(m, system, histshocks_shockdecseq, start_date, end_date, cond_type, n_back = n_back, back_shocks = back_shocks) :
                shock_decompositions_sequence(m, system, histshocks_shockdecseq, n_back = n_back, back_shocks = back_shocks)
        end

        forecast_output[:shockdecseqstates] = shockdecseqstates
        forecast_output[:shockdecseqobs]    = shockdecseqobs
        forecast_output[:shockdecseqpseudo] = shockdecseqpseudo
    end

    ### Return only desired output_vars

    for key in keys(forecast_output)
        if !(key in output_vars)
            delete!(forecast_output, key)
        end
    end

    if testing_carter_kohn && input_type == :full && get_setting(m, :forecast_smoother) == :carter_kohn
        forecast_output[:conded] = conded
    end
    return forecast_output
end


text:
module Nested

export nested

"""
    nested(T::Type, path, expr_builder, combiner=default_combiner)
Builds an arbitrary expression from each nested field in the passed in object.

Arguments:
- `T`: the type of the current object
- `P`: the type of the parent object
- `expr_builder`: function that returns an expression given T, path and a field name or index
- `combiner`: function that handles the collected expressions
"""
nested(T::Type, expr_builder, expr_combiner=default_combiner) =
    nested(T, Nothing, expr_builder, expr_combiner)
nested(T::Type, P::Type, expr_builder, expr_combiner) =
    expr_combiner(T, [Expr(:..., expr_builder(T, fn)) for fn in fieldnames(T)])

default_combiner(T, expressions) = Expr(:tuple, expressions...)

end # module

ext: jl


text:
using Printf
a, b = parse.(Int, split(readline()))
@printf("%d %d %0.8f\n", div(a,b), a%b, a/b)


text:
# # Example 2 : Lazily Read
#
# Two-dimensional fields are read from the netcdf file generated in `Example1`, and then re-written to a new netcdf file.

using NCTiles,NetCDF

inputs=NCTiles.NCTILES_TESTCASES
NCTiles.ensure_testcases_installed()

outputs = joinpath(tempdir(),"NCTILES_TESTCASES_OUTPUT/")
if ~ispath(outputs); mkpath(outputs); end

# Using NCDatasets backend
file_in=outputs*"ex1/ex1_NCDatasets.nc"
~isfile(file_in) ? error("Running Example1 first is needed to run Example2") : nothing

file_out=outputs*"ex2/ex2_NCDatasets.nc"
if ~ispath(outputs*"ex2/"); mkpath(outputs*"ex2/"); end

README = readlines(joinpath(inputs,"README"))

# Get all the metadata from the file and set up NCvars
ncvars,ncdims,fileatts = readncfile(file_in)

# Rewrite to a file
write(ncvars,file_out,README=README,globalattribs=fileatts)

# Using NetCDF backend
file_out=outputs*"ex2/ex2_NetCDF.nc"
ncvars,ncdims,fileatts = readncfile(file_in,NetCDF)
write(ncvars,file_out,README=README,globalattribs=fileatts)


text:
"""Enables deserialization of VariableCost. The default implementation can't figure out the
variable Union.
"""
function JSON2.read(io::IO, ::Type{VariableCost})
    data = JSON2.read(io)
    if data.cost isa Real
        return VariableCost(Float64(data.cost))
    elseif data.cost[1] isa Array
        variable = Vector{Tuple{Float64, Float64}}()
        for array in data.cost
            push!(variable, Tuple{Float64, Float64}(array))
        end
    else
        @assert data.cost isa Tuple || data.cost isa Array
        variable = Tuple{Float64, Float64}(data.cost)
    end

    return VariableCost(variable)
end


text:
@inline function data_matrix(points::Array{Float64, 2},
                             choices::Array{Int64, 2})
    points[:, choices]
end

@inline function construct_kernel(logσ², logϵ², logℓ)
    # x ∈ R^{parameters, latents}
    σ² = exp(logσ²)
    ϵ² = exp(logϵ²)
    ℓ  = exp.(-logℓ)
    k  = KernelFunctions.Matern52Kernel()
    t  = KernelFunctions.ARDTransform(ℓ)
    ϵ  = ϵ²*KernelFunctions.EyeKernel()
    k  = σ²*KernelFunctions.transform(k, t) + ϵ
    k
end

function compute_gram_matrix(X, logσ², logϵ², logℓ)
    kernel = construct_kernel(logσ², logϵ², logℓ)
    K = KernelFunctions.kernelmatrix(kernel, X, obsdim=2)
end


text:
@doc raw"""
read 1 byte as `UInt8` from `io`.

# Examples

```jldoctest
julia> Pickle.read_uint1(IOBuffer(b"\xff"))
0xff

```
"""
read_uint1(io::IO) = read(io, UInt8)

@doc raw"""
read 2 bytes as `UInt16` from `io`.

# Examples
```jldoctest
julia> Pickle.read_uint2(IOBuffer(b"\xff\x00"))
0x00ff

julia> Pickle.read_uint2(IOBuffer(b"\xff\xff"))
0xffff
```
"""
read_uint2(io::IO) = read(io, UInt16)

@doc raw"""
read 4 bytes as `Int32` from `io`.

# Examples
```jldoctest
julia> Pickle.read_int4(IOBuffer(b"\xff\x00\x00\x00"))
255

julia> Pickle.read_int4(IOBuffer(b"\x00\x00\x00\x80")) == -2^31
true
```
"""
read_int4(io::IO) = read(io, Int32)

@doc raw"""
read 4 bytes as `UInt32` from `io`.

# Examples
```jldoctest
julia> Pickle.read_uint4(IOBuffer(b"\xff\x00\x00\x00"))
0x000000ff

julia> Pickle.read_uint4(IOBuffer(b"\x00\x00\x00\x80")) == 2^31
true
```
"""
read_uint4(io::IO) = read(io, UInt32)

@doc raw"""
read 8 bytes as `UInt64` from `io`.

# Examples
```jldoctest
julia> Pickle.read_uint8(IOBuffer(b"\xff\x00\x00\x00\x00\x00\x00\x00"))
0x00000000000000ff

julia> Pickle.read_uint8(IOBuffer(b"\xff\xff\xff\xff\xff\xff\xff\xff")) == UInt64(2^64) - UInt64(1)
true
```
"""
read_uint8(io::IO) = read(io, UInt64)

@doc raw"""
read a string with quotes end with "\n" (newline) from `io`.

# Examples
```jldoctest
julia> Pickle.read_stringnl(IOBuffer(b"'abcd'\nefg\n"))
"abcd"

julia> Pickle.read_stringnl(IOBuffer(b"\n"))
ERROR: no string quotes around
[...]

julia> Pickle.read_stringnl(IOBuffer(b"\n"), stripquotes=false)
""

julia> Pickle.read_stringnl(IOBuffer(b"''\n"))
""

julia> Pickle.read_stringnl(IOBuffer(b"\\"abcd\\""))
ERROR: no newline found when trying to read stringnl
[...]
```
"""
function read_stringnl(io::IO; decode=true, stripquotes=true)
    data = readline(io; keep=true)
    data[end] != '\n' && error("no newline found when trying to read stringnl")

    if endswith(data, "\r\n") # handle windows newline
        nl_idx = lastindex(data) - 1
    else
        nl_idx = lastindex(data)
    end

    if stripquotes
        isquote(c) = isequal('\'')(c) || isequal('\"')(c)
        start = first(data)
        !isquote(start) && error("no string quotes around $data")

        endl = data[nl_idx-1]
        !isequal(start, endl) && error("string quote $start not found at both end of $data")

        data = data[2:nl_idx-2]
    else
        data = data[1:nl_idx-1]
    end

    if decode
        return unescape_string(data)
    end
    return data
end

read_stringnl_noescape(io::IO) = read_stringnl(io, stripquotes=false)

@doc raw"""
read a pair of unescape string from `io`.

# Examples
```jldoctest
julia> Pickle.read_stringnl_noescape_pair(IOBuffer(b"Queue\nEmpty\njunk"))
("Queue", "Empty")
```
"""
read_stringnl_noescape_pair(io::IO) = (read_stringnl_noescape(io), read_stringnl_noescape(io))

function read_multiple(io::IO, name::String, nf::Function)
    n = nf(io)
    data = read(io, n)
    length(data) != n && error("expected $n bytes in $name, but only $(length(data)) remain")
    return data
end
read_multiple(io::IO, T::Type, name::String, nf::Function) = T(read_multiple(io, name, nf))

@doc raw"""
read 1 byte as `UInt8` for length n and read n bytes as `String` from `io`.`

# Examples
```jldoctest
julia> Pickle.read_string1(IOBuffer(b"\x00"))
""

julia> Pickle.read_string1(IOBuffer(b"\x03abcdef"))
"abc"
```
"""
read_string1(io::IO) = read_multiple(io, String, "string1", read_uint1)

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_string4(IOBuffer(b"\x00\x00\x00\x00abc"))
""

julia> Pickle.read_string4(IOBuffer(b"\x03\x00\x00\x00abcdef"))
"abc"

julia> Pickle.read_string4(IOBuffer(b"\x00\x00\x00\x03abcdef"))
ERROR: expected 50331648 bytes in string4, but only 6 remain
[...]
```
"""
read_string4(io::IO) = read_multiple(io, String, "string4", read_int4)

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_bytes1(IOBuffer(b"\x00")) |> isempty # 0-element Vector{UInt8}
true

julia> Pickle.read_bytes1(IOBuffer(b"\x03abcdef"))
3-element Vector{UInt8}:
 0x61
 0x62
 0x63

julia> Pickle.read_bytes1(IOBuffer(b"\x03abcdef")) |> String
"abc"
```
"""
read_bytes1(io::IO) = read_multiple(io, "bytes1", read_uint1)

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_bytes4(IOBuffer(b"\x00\x00\x00\x00abc")) |> isempty # 0-element Vector{UInt8}
true

julia> Pickle.read_bytes4(IOBuffer(b"\x03\x00\x00\x00abcdef"))
3-element Vector{UInt8}:
 0x61
 0x62
 0x63

julia> Pickle.read_bytes4(IOBuffer(b"\x00\x00\x00\x03abcdef"))
ERROR: expected 50331648 bytes in bytes4, but only 6 remain
[...]
```
"""
read_bytes4(io::IO) = read_multiple(io, "bytes4", read_uint4)

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_bytes8(IOBuffer(b"\x00\x00\x00\x00\x00\x00\x00\x00abc")) |> isempty # 0-element Vector{UInt8}
true

julia> Pickle.read_bytes8(IOBuffer(b"\x03\x00\x00\x00\x00\x00\x00\x00abcdef"))
3-element Vector{UInt8}:
 0x61
 0x62
 0x63
```
"""
read_bytes8(io::IO) = read_multiple(io, "bytes8", read_uint8)

########################
# this function only make sense in python.

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_bytearray8(IOBuffer(b"\x00\x00\x00\x00\x00\x00\x00\x00abc")) |> isempty # 0-element Vector{UInt8}
true

julia> Pickle.read_bytearray8(IOBuffer(b"\x03\x00\x00\x00\x00\x00\x00\x00abcdef"))
3-element Vector{UInt8}:
 0x61
 0x62
 0x63
```
"""
read_bytearray8(io::IO) = read_multiple(io, "bytearray8", read_uint8)

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_unicodestringnl(IOBuffer(b"abc\\uabcd\njunk")) == "abc\uabcd"
true
```
"""
function read_unicodestringnl(io::IO)
    data = readline(io; keep=true)
    data[end] != '\n' && error("no newline found when trying to read unicodestringnl")

    if endswith(data, "\r\n") # handle windows newline
        nl_idx = lastindex(data) - 1
    else
        nl_idx = lastindex(data)
    end

    data = data[1:nl_idx-1]
    return unescape_string(data)
end

read_unicodestring1(io::IO)  = read_multiple(io, String, "unicodestring1", read_uint1)
read_unicodestring4(io::IO)  = read_multiple(io, String, "unicodestring4", read_int4)
read_unicodestring8(io::IO)  = read_multiple(io, String, "unicodestring8", read_uint8)
########################

function parseint(s; base=10)
    i = tryparse(Int, s; base=base)
    !isnothing(i) && return i
    bi = tryparse(BigInt, s; base=base)
    !isnothing(bi) && return bi
    return nothing
end

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_decimalnl_short(IOBuffer(b"1234\n56"))
1234

julia> Pickle.read_decimalnl_short(IOBuffer(b"1234L\n56"))
ERROR: invalid literal for Integer with base 10: 1234L
[...]
```
"""
function read_decimalnl_short(io::IO)
    s = read_stringnl(io, decode=false, stripquotes=false)

    if isequal("00")(s)
        return false
    elseif isequal("01")(s)
        return true
    end

    int = parseint(s)
    isnothing(int) &&
        error("invalid literal for Integer with base 10: $s")
    return int
end

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_decimalnl_long(IOBuffer(b"1234L\n56"))
1234

julia> Pickle.read_decimalnl_long(IOBuffer(b"123456789012345678901234L\n6"))
123456789012345678901234
```
"""
function read_decimalnl_long(io::IO)
    s = read_stringnl(io, decode=false, stripquotes=false)
    if s[end] == 'L'
        s = s[1:end-1]
    end

    int = parseint(s)
    isnothing(int) &&
        error("invalid literal for Integer with base 10: $s")
    return int
end

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_floatnl(IOBuffer(b"-1.25\n6"))
-1.25
```
"""
function read_floatnl(io::IO)
    s = read_stringnl(io, decode=false, stripquotes=false)
    return parse(Float64, s)
end

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_float8(IOBuffer(b"\xbf\xf4\x00\x00\x00\x00\x00\x00\n"))
-1.25
```
"""
read_float8(io::IO) = bswap(read(io, Float64))

@doc raw"""
convert a 2's complement byte array (Vector{UInt8}) into Signed Integer.

# Examples
```jldoctest
julia> Pickle.int_from_bytes(b"")
0

julia> Pickle.int_from_bytes(b"\xff\x00")
255

julia> Pickle.int_from_bytes(b"\xff\x7f")
32767

julia> Pickle.int_from_bytes(b"\x00\xff")
-256

julia> Pickle.int_from_bytes(b"\x00\x80")
-32768

julia> Pickle.int_from_bytes(b"\x80")
-128

julia> Pickle.int_from_bytes(b"\x7f")
127
```
"""
function int_from_bytes(x)
    isempty(x) && return 0
    islittle_endian() && (x = reverse(x))
    isneg = isone(first(x) >> 7)
    if !isneg
        s = Base.bytes2hex(x)
        i = parseint(s; base=16)
        isnothing(i) && error("2's complement bytes parse error")
        return i
    else
        s = Base.bytes2hex(x .⊻ 0xff)
        i = parseint(s; base=16)
        isnothing(i) && error("2's complement bytes parse error")
        return -i-1
    end
end

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_long1(IOBuffer(b"\x00"))
0

julia> Pickle.read_long1(IOBuffer(b"\x02\xff\x00"))
255

julia> Pickle.read_long1(IOBuffer(b"\x02\xff\x7f"))
32767

julia> Pickle.read_long1(IOBuffer(b"\x02\x00\xff"))
-256

julia> Pickle.read_long1(IOBuffer(b"\x02\x00\x80"))
-32768
```
"""
read_long1(io::IO) = int_from_bytes(read_multiple(io, "bytes1", read_uint1))

@doc raw"""
# Examples
```jldoctest
julia> Pickle.read_long4(IOBuffer(b"\x02\x00\x00\x00\xff\x00"))
255

julia> Pickle.read_long4(IOBuffer(b"\x02\x00\x00\x00\xff\x7f"))
32767

julia> Pickle.read_long4(IOBuffer(b"\x02\x00\x00\x00\x00\xff"))
-256

julia> Pickle.read_long4(IOBuffer(b"\x02\x00\x00\x00\x00\x80"))
-32768

julia> Pickle.read_long4(IOBuffer(b"\x00\x00\x00\x00"))
0
```
"""
read_long4(io::IO) = int_from_bytes(read_multiple(io, "bytes4", read_int4))


text:
#
# NuvuCameras.jl -
#
# Provide Julia interface to the Nüvü Camēras SDK.
#

__precompile__(true)

module NuvuCameras

# Low level interface.  All functions, types and constants are made available
# with the prefix `NC.`.
module NC
include("nctypes.jl")
include("ncconst.jl")
include("nccalls.jl")
end
import .NC

# Private module with stuff hidden to the end-user.  All methods, types and
# constants are made available with the prefix `Priv.`.
module Priv
include("private.jl")
end
import .Priv

# High level interface.
include("types.jl")
include("public.jl")

end


text:
"""
Invert Geotransform.

This function will invert a standard 3x2 set of GeoTransform coefficients. This
converts the equation from being pixel to geo to being geo to pixel.

### Parameters
* `gt_in`       Input geotransform (six doubles - unaltered).
* `gt_out`      Output geotransform (six doubles - updated).

### Returns
`gt_out`
"""
function invgeotransform!(gt_in::Vector{Cdouble}, gt_out::Vector{Cdouble})
    result = Bool(GDAL.invgeotransform(pointer(gt_in), pointer(gt_out)))
    result || error("Geotransform coefficients is uninvertable")
    gt_out
end

invgeotransform(gt_in::Vector{Cdouble}) =
    invgeotransform!(gt_in, Array{Cdouble}(undef, 6))

"""
Apply GeoTransform to x/y coordinate.

Applies the following computation, converting a (pixel,line) coordinate into a
georeferenced `(geo_x,geo_y)` location.
```
    geo_x = geotransform[1] + pixel*geotransform[2] + line*geotransform[3]
    geo_y = geotransform[4] + pixel*geotransform[5] + line*geotransform[6]
```
### Parameters
* `geotransform`  Six coefficient GeoTransform to apply.
* `pixel`           input pixel position.
* `line`            input line position.
"""
function applygeotransform(geotransform::Vector{Cdouble},
                           pixel::Cdouble,
                           line::Cdouble)
    geo_xy = Array{Cdouble}(undef, 2)
    geo_x = pointer(geo_xy)
    geo_y = geo_x + sizeof(Cdouble)
    GDAL.applygeotransform(pointer(geotransform), pixel, line, geo_x, geo_y)
    geo_xy
end

"""
Compose two geotransforms.

The resulting geotransform is the equivelent to `padfGT1` and then `padfGT2`
being applied to a point.

### Parameters
* `gt1`     the first geotransform, six values.
* `gt2`     the second geotransform, six values.
* `gtout`   the output geotransform, six values, may safely be the same
array as `gt1` or `gt2`.
"""
function composegeotransform!(gt1::Vector{Cdouble}, gt2::Vector{Cdouble},
                              gtout::Vector{Cdouble})
    GDAL.composegeotransform(pointer(gt1), pointer(gt2), pointer(gtout))
    gtout
end

composegeotransform(gt1::Vector{Cdouble}, gt2::Vector{Cdouble}) =
    composegeotransform!(gt1, gt2, Array{Cdouble}(undef, 6))


text:
abstract type AbstractFVMesh end
abstract type AbstractFVBoundaryCondition end
abstract type AbstractFVProbType end
struct GeneralProblem <: AbstractFVProbType end
struct ScalarProblem <: AbstractFVProbType end

# Reference for cell and node indexes
#   1   2   3          N-1  N
# |---|---|---|......|---|---|
# 1   2   3   4 ... N-1  N  N+1

# Wrapper type used to manage boundary conditions
struct FVMesh1D{PType <: AbstractFVProbType, MType, LBType,RBType} <: AbstractFVMesh
  mesh::MType
  left_boundary::LBType
  right_boundary::RBType
end

function mesh_setup(mesh, dbcs::Vector{T}, PType::AbstractFVProbType) where {T<:AbstractFVBoundaryCondition}
  dim = getdimension(mesh)
  if  dim == 1
    left_bd = nothing
    right_bd = nothing
    for bdc in dbcs
        if isLeftCondition(bdc)
            left_bd = bdc
        end
        if isRightCondition(bdc)
            right_bd = bdc
        end
    end
    return FVMesh1D{typeof(PType), typeof(mesh), typeof(left_bd), typeof(right_bd)}(mesh, left_bd, right_bd)
  else
    error("No methods available for mesh of dimension ", dim)
  end
end

@inline cell_indices(fvmesh::AbstractFVMesh) = 1:getncells(fvmesh.mesh)
@inline node_indices(fvmesh::AbstractFVMesh) = 1:getnnodes(fvmesh.mesh)
@inline cell_volume(fvmesh::AbstractFVMesh, cell_idx::Int) = cell_volume(fvmesh.mesh, cell_idx)
@inline value_at_cell(u, j,fvmesh::FVMesh1D{ScalarProblem}) = u[j]
@inline value_at_cell(u, j,fvmesh::FVMesh1D{GeneralProblem}) = u[:,j]

"""
    cellval_at_left(node::Int, A::AbstractArray{T,2}, mesh::AbstractFVMesh1D) where {T}

   cell values of variable `A` to the left of `node` in `mesh`.
"""
function cellval_at_left(node::Int, A, fvmesh::FVMesh1D{ScalarProblem})
    checkbounds(Bool, A, node-1) && return A[node-1]
    return A[cellidx_at_left(node, getncells(fvmesh.mesh), fvmesh.left_boundary)]
end

function cellval_at_left(node::Int, A, fvmesh::FVMesh1D{GeneralProblem})
    checkbounds(Bool, A, (:,node-1)...) && return A[:,node-1]
    return A[:,cellidx_at_left(node, getncells(fvmesh.mesh), fvmesh.left_boundary)]
end

"""
    cellval_at_right(edge::Int, A::AbstractArray{T,2}, mesh::AbstractFVMesh1D) where {T}

cell values of variable `A` to the right of `node` in `mesh`.
"""
function cellval_at_right(node::Int, A, fvmesh::FVMesh1D{ScalarProblem})
    checkbounds(Bool, A, node) && return A[node]
    return A[cellidx_at_right(node, getncells(fvmesh.mesh), fvmesh.right_boundary)]
end

function cellval_at_right(node::Int, A, fvmesh::FVMesh1D{GeneralProblem})
    checkbounds(Bool, A, (:,node)...) && return A[:,node]
    return A[:,cellidx_at_right(node, getncells(fvmesh.mesh), fvmesh.right_boundary)]
end

"""
    get_cellvals(A::AbstractArray{T,2}, idx..., mesh::AbstractFVMesh1D) where {T}
   cell values of variable `A` on cells `idx` of `mesh`.
"""
function get_cellvals(A, fvmesh::FVMesh1D{GeneralProblem}, idx...) where {T}
    checkbounds(Bool, A, idx...) && return A[idx...]
    if (minimum(idx[2]) < 1)
        return A[idx[1],getIndex(idx[2], size(A,2), fvmesh.left_boundary)]
    elseif (maximum(idx[2]) > getncells(fvmesh.mesh))
        return A[idx[1],getIndex(idx[2], size(A,2), fvmesh.right_boundary)]
    else
        error("unknown index")
    end
end

function get_cellvals(A, fvmesh::FVMesh1D{ScalarProblem}, idx...) where {T}
    checkbounds(Bool, A, idx...) && return A[idx...]
    if (minimum(idx[1]) < 1)
        return A[getIndex(idx[1], size(A,1), fvmesh.left_boundary)]
    elseif (maximum(idx[1]) > getncells(fvmesh.mesh))
        return A[getIndex(idx[1], size(A,1), fvmesh.right_boundary)]
    else
        error("unknown index")
    end
end

"""
    left_edge(cell::Int, mesh::AbstractFVMesh1D)

The index of the node to the left of `cell` in `mesh`.
"""
function left_node_idx(cell::Int, fvmesh::FVMesh1D)
    @boundscheck begin
        @assert (1 <= cell <= getncells(fvmesh.mesh))
    end
    cell
end

"""
    right_edge(cell::Int, mesh::AbstractFVMesh1D)

The index of the node to the right of `cell` in `mesh`.
"""
function right_node_idx(cell::Int, fvmesh::FVMesh1D)
    @boundscheck begin
        @assert (1 <= cell <= getncells(fvmesh.mesh))
    end
    cell+1
end

function left_cell_idx(node::Int, fvmesh::FVMesh1D)
    @boundscheck begin
        @assert (1 <= node <= getnnodes(fvmesh.mesh))
    end
    cellidx_at_left(node, getncells(fvmesh.mesh), fvmesh.left_boundary)
end
function right_cell_idx(node::Int, fvmesh::FVMesh1D)
    @boundscheck begin
        @assert (1 <= node <= getnnodes(fvmesh.mesh))
    end
    cellidx_at_right(node, getncells(fvmesh.mesh), fvmesh.right_boundary)
end

function apply_bc_in_fluxes!(fluxes, fvmesh)
    apply_lbc_in_fluxes!(fluxes, fvmesh, fvmesh.left_boundary)
    apply_rbc_in_fluxes!(fluxes, fvmesh, fvmesh.right_boundary)
end

function apply_bc_in_du!(du, fvmesh)
    apply_lbc_in_du!(du, fvmesh, fvmesh.left_boundary)
    apply_rbc_in_du!(du, fvmesh, fvmesh.right_boundary)
end

apply_lbc_in_fluxes!(fluxes, fvmesh, ::AbstractFVBoundaryCondition) = nothing
apply_rbc_in_fluxes!(fluxes, fvmesh, ::AbstractFVBoundaryCondition) = nothing

apply_lbc_in_du!(du, fvmesh, ::AbstractFVBoundaryCondition) = nothing
apply_rbc_in_du!(du, fvmesh, ::AbstractFVBoundaryCondition) = nothing

################################
# Periodic boundary conditions
###############################
struct Periodic <: AbstractFVBoundaryCondition
    axis::Integer
end
function Periodic(;axis = 1)
    Periodic(axis)
end
isLeftCondition(::Periodic) = true
isRightCondition(::Periodic) = true

function getIndex(I, n::Int, ::Periodic)
    if typeof(I) <: Int
      return mod1(I, n)
    else
      return [mod1(i, n) for i in I]
    end
end

cellidx_at_left(node::Int, n::Int, ::Periodic) = mod1(node-1, n)
cellidx_at_right(node::Int, n::Int, ::Periodic) =  mod1(node, n)

################################
# Zero flux boundary conditions
###############################
struct ZeroFlux <: AbstractFVBoundaryCondition
    side::Symbol
end
ZeroFlux(;side=:Both) = ZeroFlux(side)

function apply_lbc_in_fluxes!(fluxes, fvmesh::FVMesh1D{GeneralProblem}, ::ZeroFlux)
    fluxes[:,1] .= zero(eltype(fluxes))
end
function apply_lbc_in_fluxes!(fluxes, fvmesh::FVMesh1D{ScalarProblem}, ::ZeroFlux)
    fluxes[1] = zero(eltype(fluxes))
end
function apply_rbc_in_fluxes!(fluxes, fvmesh::FVMesh1D{GeneralProblem}, ::ZeroFlux)
    fluxes[:, getnnodes(fvmesh.mesh)] .= zero(eltype(fluxes))
end
function apply_rbc_in_fluxes!(fluxes, fvmesh::FVMesh1D{ScalarProblem}, ::ZeroFlux)
    fluxes[getnnodes(fvmesh.mesh)] = zero(eltype(fluxes))
end

################################
# Dirichlet boundary conditions
###############################
struct Dirichlet <: AbstractFVBoundaryCondition
    side::Symbol
end
Dirichlet(;side=:Both) = Dirichlet(side)
isLeftCondition(c::Union{ZeroFlux,Dirichlet}) = c.side==:Both || c.side == :Left
isRightCondition(c::Union{ZeroFlux,Dirichlet}) = c.side==:Both || c.side == :Right

function getIndex(I, n::Int, ::Union{ZeroFlux,Dirichlet})
    if typeof(I) <: Int
      return min(n,max(1,I))
    else
      return [min(n,max(1,i)) for i in I]
    end
end

cellidx_at_left(node::Int, n::Int, ::Union{ZeroFlux,Dirichlet}) = min(n,max(1,node-1))
cellidx_at_right(node::Int, n::Int, ::Union{ZeroFlux,Dirichlet}) =  min(n,max(1,node))

function apply_lbc_in_du!(du, fvmesh::FVMesh1D{GeneralProblem}, ::Dirichlet)
    du[:,1] .= zero(eltype(du))
end
function apply_lbc_in_du!(du, fvmesh::FVMesh1D{ScalarProblem}, ::Dirichlet)
    du[1] = zero(eltype(du))
end

function apply_rbc_in_du!(du, fvmesh::FVMesh1D{GeneralProblem}, ::Dirichlet)
    du[:,getncells(fvmesh.mesh)] .= zero(eltype(du))
end

function apply_rbc_in_du!(du, fvmesh::FVMesh1D{ScalarProblem}, ::Dirichlet)
    du[getncells(fvmesh.mesh)] = zero(eltype(du))
end


text:
function argmax(x::Array, dim::Int)
    _, index = findmax(x, dim)
    ind2sub(size(x), vec(index))[dim]
end

T = Float32
x = rand(T, 100, 100)
cux = CuArray{T}(x)
@test argmax(x) == Array(argmax(cux))