text:
function document(spec, p)
    @match p[:category] begin
        :function => document_function(spec, p)
        :struct => document_struct(spec, p)
    end
end

document(spec::SpecHandle) = document.(Ref(spec), add_constructors(spec))

backquoted(arg) = string('`', arg, '`')

function document_return_codes(spec::SpecFunc)
    res = ""
    if !isempty(spec.success_codes) || !isempty(spec.error_codes)
        res *= "Return codes:"
    end
    if !isempty(spec.success_codes) && must_return_success_code(spec)
        res *= string("\n- Success:\n  - ", join(backquoted.(spec.success_codes), "\n  - "))
    end
    if !isempty(spec.error_codes)
        res *= string("\n- Error:\n  - ", join(backquoted.(spec.error_codes), "\n  - "))
    end

    if !isempty(res)
        res *= '\n'^2
    end

    res
end

document_return_codes(spec::SpecHandle) = document_return_codes(add_constructor(spec))
document_return_codes(spec::SpecStruct) = ""

function document_function(spec::SpecHandle, p)
    docstring(reconstruct_call(p), string(' '^4, reconstruct_call(p), '\n'^2))
end

function append_to_argdoc!(argdocs, spec, str)
    i = findfirst(==(wrap_identifier(spec)), first.(argdocs))
    argdocs[i] = concat_right(argdocs[i], str)
end

function document_function(spec, p)
    params = children(spec)
    externsync_params = filter(x -> x.is_externsync, params)
    argdocs = document_arguments(p)
    if !isempty(externsync_params)
        for param ∈ externsync_params
            append_to_argdoc!(argdocs, param, " (externsync)")
        end
    end
    extra = ""
    if spec isa SpecFunc
        if any(is_data_with_retrievable_size, params)
            extra = """
            !!! warning
                The pointer returned by this function holds memory owned by Julia. It is therefore **your** responsibility to free it after use (e.g. with `Libc.free`)."""
        elseif any(is_opaque_data, params)
            opaque_params = params[findall(is_opaque_data, params)]
            map(opaque_params) do param
                append_to_argdoc!(
                    argdocs,
                    param,
                    " (must be a valid pointer with `$(wrap_identifier(len(param)))` bytes)",
                )
            end
        end
    end

    docstring(
        p[:name],
        string(
            ' '^4,
            reconstruct_call(p),
            '\n'^2,
            document_return_codes(spec),
            join(["Arguments:"; last.(argdocs)], "\n- "),
            isempty(extra) ? "" : '\n'^2 * extra,
            '\n',
        ),
    )
end

function document_arguments(p)
    args, kwargs = p[:args], p[:kwargs]
    if isempty(args) && isempty(kwargs)
        []
    else
        backquoted_args = map(x -> name(x) => backquoted(x), args)
        documented_kwargs = map(kwargs) do kwarg
            identifier, value = kwarg.args
            identifier => string(backquoted(identifier), ": defaults to ", backquoted(value))
        end
        vcat(backquoted_args, documented_kwargs)
    end
end

function document_struct(spec::SpecStruct, p)
    error("Not implemented")
end

function docstring(name, docstring)
    p_with_docs = Dict(:category => :doc, :docstring => string('\n', docstring, '\n'), :ex => name)

    reconstruct(p_with_docs)
end

concat_right(pair::Pair, val) = pair.first => (pair.second * val)

ext: jl


text:
module _SampEn2D
export SampEn2D
using Statistics: mean, std
    """
        SE2D, Phi1, Phi2 = SampEn2D(Mat)

    Returns the bidimensional sample entropy estimate (`SE2D`) and the number
    of matched sub-matricess (m:Phi1, m+1:Phi2) estimated for the data
    matrix (`Mat`) using the default parameters: time delay = 1,
    radius distance threshold = 0.2*SD(`Mat`), logarithm = natural
    matrix template size = [floor(H/10) floor(W/10)], (where H and W represent
    the height (rows) and width (columns) of the data matrix `Mat`) \n
    ** The minimum dimension size of Mat must be > 10.**

        SE2D, Phi1, Phi2 = SampEn2D(Mat::AbstractArray{T,2} where T<:Real; m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10),
                                        tau::Int=1, r::Real=0.2*std(Mat,corrected=false), Logx::Real=exp(1), Lock::Bool=true)

    Returns the bidimensional sample entropy (`SE2D`) estimates for the data
    matrix (`Mat`) using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Template submatrix dimensions, an integer scalar (i.e. the same
              height and width) or a two-element vector of integers
              [height, width] with a minimum value > 1.
              (default: [floor(H/10) floor(W/10)]) \n
    `tau`   - Time Delay, a positive integer       (default: 1) \n
    `r`     - Distance Threshold, a positive scalar  (default: 0.2*SD(Mat)) \n
    `Logx`  - Logarithm base, a positive scalar    (default: natural)  \n
    `Lock`  - By default, SampEn2D only permits matrices with a maximum
              size of 128 x 128 to prevent memory errors when storing data
              on RAM. e.g. For Mat = [200 x 200], m = 3, and tau = 1,
              SampEn2D creates a vector of 753049836 elements.
              To enable matrices greater than [128 x 128] elements,
              set `Lock` to false.  (default: true)  \n
              `WARNING: unlocking the permitted matrix size may cause your Julia
              IDE to crash.`

    #  See also `SampEn`, `FuzzEn2D`, `XSampEn`, `MSEn`

    # References:
        [1] Luiz Eduardo Virgili Silva, et al.,
             "Two-dimensional sample entropy: Assessing image texture
             through irregularity."
             Biomedical Physics & Engineering Express
             2.4 (2016): 045002.

    """
    function SampEn2D(Mat::AbstractArray{T,2} where T<:Real;
        m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10),
        tau::Int=1, r::Real=0.2*std(Mat,corrected=false), Logx::Real=exp(1), Lock::Bool=true)

    NL, NW = size(Mat)
    ((NL > 128 || NW > 128) && Lock) ?
        error("To prevent memory errors, matrix width & length must have <= 128 elements.
            To estimate SampEn2D  for the current matrix ($NL,$NW) change Lock to 'false'.
            Caution: unlocking the safe matrix size may cause the Julia IDE to crash.") :
            nothing

    length(m)==1 ?  (mL = m; mW = m) : (mL = m[1];  mW = m[2])
    (NL > 10 && NW > 10) ? nothing :
        error("Number of rows and columns in Mat must be > 10")
    (minimum(m)>1) ? nothing :
        error("m:     must be an integer > 1, or a 2 element tuple of integer values > 1")
    (tau > 0) ? nothing : error("tau:   must be an integer > 0")
    (r >= 0) ?  nothing : error("r:     must be a positive value")
    (Logx>0) ? nothing :  error("Logx:     must be a positive number > 0")

    NL = NL - mL*tau
    NW = NW - mW*tau
    X = zeros(mL+1,mW+1,NL*NW)
    p = 0
    for k = 1:NL
        for n = 1:NW
            p += 1
            X[:,:,p] = Mat[k:tau:(mL)*tau+k,n:tau:(mW)*tau+n]
        end
    end
    p = size(X,3)
    p != NL*NW ? @warn("Potential error with submatrix division.") : nothing
    Ny = p*(p-1)/2
    Ny > 300000000 ? @warn("Number of pairwise distance calculations is $Ny") : nothing

    Y1 = zeros(p-1)
    Y2 = zeros(p-1)
    for k = 1:p-1
        Temp = maximum(abs.(X[1:mL,1:mW,k+1:p] .- X[1:mL,1:mW,k]), dims=(1,2))[:] .< r
        Y1[k] = sum(Temp);    Temp = findall(Temp.>0) .+ k
        Y2[k] = sum(maximum(abs.(X[:,:,Temp] .- X[:,:,k]),dims=(1,2)) .< r)
    end
    Phi1 = sum(Y1)/Ny
    Phi2 = sum(Y2)/Ny
    SE2D = -log(Logx, Phi2/Phi1)

    return SE2D, Phi1, Phi2
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
export @parseParameters_str,@parseEquations_str

macro parseParameters_str(str)
  print(str,"\n")
  str=split(str,'\n')
  map(x->eval(parse(x)),str)
end

macro parseEquations_str(str)
  eqs=map(strip,split(str,'\n'));
  global vars=[]
  diffEqs=[]

  for k=1:length(eqs)
    if eqs[k][1]=='d'
      vars=[vars;match(r"(?<=d).*(?=/dt)",eqs[k]).match]
      str=vectorize(match(r"(?<==).*$",eqs[k]).match)
      diffEqs=[diffEqs;str]
    else
      str=vectorize(eqs[k])
      print(str)
      eval(parse(str));
    end
  end

  argument=join([vars;"t"],",");
  for k=1:length(vars)
    str=*("F",vars[k],"($argument)","=",diffEqs[k])
    print(str)
    eval(parse(str))
  end
  print("\n")
end


text:
module Bibliography

# BibInternal
import BibInternal
import BibInternal: AbstractEntry, Entry

# BibParser
import BibParser, BibParser.BibTeX

# Others
import DataStructures
import DataStructures.OrderedSet

export export_bibtex, import_bibtex
export export_web, bibtex_to_web
export select

include("select.jl")
include("bibtex.jl")
include("csl.jl")
include("staticweb.jl")

"""
    export_bibtex(e::Entry)
    export_bibtex(bibliography::DataStructures.OrderedDict{String,Entry})
    export_bibtex(target::String, bibliography::DataStructures.OrderedDict{String,Entry})
Export an entry or a bibliography to BibTeX format.
"""
function export_bibtex(target::String, bibliography::DataStructures.OrderedDict{String,Entry})
    data = export_bibtex(bibliography)
    if target != ""
        f = open(target, "w")
        write(f, data)
        close(f)
    end
    return data
end

"""
    bibtex_to_web(source::String)
Convert a BibTeX file to a web compatible format, specifically for the [StaticWebPages.jl](https://github.com/Azzaare/StaticWebPages.jl) pakcage.
"""
function bibtex_to_web(source::String)
    export_web(import_bibtex(source))
end

end # module


text:
using Test

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


text:
using PackageCompiler
create_sysimage(:CSV, sysimage_path="sys_mycsv.so", precompile_execution_file="precompile_csv.jl")

max_stars_repo_name: alsam/scriptics-samples

max_issues_repo_name: alsam/scriptics-samples

max_forks_repo_name: alsam/scriptics-samples


text:
"""
    readmps(filename)

    Parse an MPS file, and return a linear program in standard form.
"""
function readmps(file_name)
    # parse an MPS file
    section = ""
    d = Dict()

    row2idx = Dict()  # name <-> index correspondence for rows
    col2idx = Dict()  # name <-> index correspondence for columns

    senses = []  # sense of rows

    coeffs_row = Array{Int, 1}(0)
    coeffs_col = Array{Int, 1}(0)
    coeffs_val = Array{Float64, 1}(0)

    obj_col = Array{Int, 1}(0)
    obj_val = Array{Float64, 1}(0)

    rhs_row = Array{Int, 1}(0)
    rhs_val = Array{Float64, 1}(0)

    lb_col = Array{Int, 1}(0)
    lb_val = Array{Float64, 1}(0)
    ub_col = Array{Int, 1}(0)
    ub_val = Array{Float64, 1}(0)

    ranges_row = Array{Int, 1}(0)
    ranges_val = Array{Float64, 1}(0)

    d["nrows"] = 0
    d["ncols"] = 0
    d["rhs"] = ""
    d["bound"] = ""
    d["range"] = ""

    open(file_name) do f
        for ln in eachline(f)

            fields = parseline(ln)

            # pass empty lines or comments
            if length(fields) == 0
                continue
            end

            # check for indicators
            if ln[1] != ' '
                if fields[1] == "NAME"
                    # extract problem name
                    if length(fields) >= 2
                        # second field is name, rest is ignored
                        d["name"] = fields[2]
                    else
                        #  empty name
                        d["name"] = ""
                    end
                    continue
                end
                # update section value and pass to next line
                section = fields[1]
                continue
            end

            if section == "ROWS"
                parsemps_rows!(ln, fields, d, row2idx, senses)

            elseif section == "COLUMNS"
                parsemps_columns!(
                    ln, fields, d, col2idx, row2idx,
                    coeffs_row, coeffs_col, coeffs_val, obj_col, obj_val
                )

            elseif section == "RHS"
                parsemps_rhs!(ln, fields, d, row2idx, rhs_row, rhs_val)

            elseif section == "RANGES"
                parsemps_ranges!(ln, fields, d, row2idx, ranges_row, ranges_val)

            elseif section == "BOUNDS"
                parsemps_bounds!(ln, fields, d, col2idx, lb_col, lb_val, ub_col, ub_val)

            elseif section == "ENDATA"
                # end of file
            else
                warn("UNEXPECTED FORMAT PARSING\t$(ln)")
            end

        end

        # End of file reached
    end

    # extract relevant info
    ranges = sparsevec(ranges_row, ranges_val)

    # transform to standard form
    m, n, obj, rhs, coeffs, lb, ub = convert_to_standard_form(
        d,
        obj_col, obj_val,
        rhs_row, rhs_val, senses,
        coeffs_row, coeffs_col, coeffs_val,
        lb_col, lb_val, ub_col, ub_val,
        ranges_row, ranges_val
    )

    return m, n, obj, rhs, coeffs, lb, ub, ranges
end

function convert_to_standard_form(
    d,
    obj_col, obj_val,  # objective
    rhs_row, rhs_val, senses,  # right-hand side and senses
    coeffs_row, coeffs_col, coeffs_val,  # coefficients
    lb_col, lb_val, ub_col, ub_val,  # lower and uper bounds
    ranges_row, ranges_val
    )

    m = d["nrows"]  # number of constraints in original formulation
    n = d["ncols"]  # number of variables in original formulation

    # I. Tranform constraints to equality constraints
    for j=1:m
        if senses[j] == "L"
            # add slack variable
            n += 1
            push!(coeffs_row, j)
            push!(coeffs_col, n)
            push!(coeffs_val, 1.0)
            senses[j] = "E"

        elseif senses[j] == "G"
            # add surplus variable
            n += 1
            push!(coeffs_row, j)
            push!(coeffs_col, n)
            push!(coeffs_val, -1.0)
            senses[j] = "E"
        end
    end

    coeffs = sparse(coeffs_row, coeffs_col, coeffs_val, m, n)
    obj = Array(sparsevec(obj_col, obj_val, n))
    rhs = Array(sparsevec(rhs_row, rhs_val, m))

    # II. Form bounds
    lb = zeros(n)  # lower bound
    ub = Inf * ones(n)  # upper bound
    for i in lb_col
        lb[i] = lb_val[i]
    end
    for i in ub_col
        ub[i] = ub_val[i]
    end

    return m, n, obj, rhs, coeffs, lb, ub
end

"""
    parseline(ln)
"""
function parseline(ln)
    s = []
    if length(ln) == 0 || ln[1] == '*'
        # empty line or comment line
        return s
    elseif ln[1] != ' '
        # indicator row
        s = split(ln)
        return s
    end

    # fill-in with spaces to get a 80-long line
    if length(ln) < 80
        ln_ = ln * (" "^(80-length(ln)))
    end

    fields = [
        ln_[2:3],
        ln_[5:12],
        ln_[15:22],
        ln_[25:36],
        ln_[40:47],
        ln_[50:61]
    ]

    # parse each field from right to left
    # empty fields on the right are ignored
    for f in fields[end:-1:1]
        # remove left-trailing and right-trailing spaces
        f_ = strip(f)  # f_ may be empty

        if length(f_) == 0 && length(s) == 0
            # f_ is empty field at the right-end
            # do nothing
            continue
        # elseif contains(f_, "\$")
        #     # f contains the $ sign: remove everything after $
        #     s = [split(f_, '$')[1]]
        #     # rest of line is treated as comment
        #     continue
        else
            # f_ is empty but there is a non-empty field later
            unshift!(s, f_)
            continue
        end
    end

    return s
end

function parsemps_rows!(ln, fields, d, row2idx, senses)
    # parse a ROWS line of the MPS file
    # `ln` is a string that contains the current line
    # `fields` is an array such that fields == split(ln)

    # current line should contain exactly 2 fields
    if length(fields) != 2
        # input error, current line is ignored
        warn("INPUT ERROR:\t$(ln)")
        return nothing
    end

    # parse line
    # First field can be either of:
    #   "N" -> row corresponds to objective
    #   "E" -> equality constraint
    #   "L" -> less-or-equal inequality constraint
    #   "G" -> greater-or-equal inequality constraint
    if fields[1] == "N"
        # objective
        row2idx[fields[2]] = 0

    elseif fields[1] == "E" || fields[1] == "L" || fields[1] == "G"
        # constraint
        d["nrows"] += 1
        row2idx[fields[2]] = d["nrows"]  # name
        push!(senses, fields[1])  # sense (i.e., E, L, or G)

    else
        # input error, current line is ignored
        warn("INPUT ERROR:\t$(ln)")
        return nothing
    end

    return nothing
end

function parsemps_columns!(ln, fields, d, col2idx, row2idx,
    coeffs_row, coeffs_col, coeffs_val, obj_col, obj_val
)
    # parse a ROWS line of the MPS file
    # `ln` is a string that contains the current line
    # `fields` contains the 6 fields of that line

    # current line should contain at least three fields
    if length(fields) < 4
        warn("COLUMN LINE TOO SHORT|$(ln)")
        return nothing
    end

    # First field is empty, second field is variable's name
    cname = fields[2]
    if !haskey(col2idx, cname)
        d["ncols"] += 1
        col2idx[cname] = d["ncols"]  # Add column name to the pool
    end

    # Second and third fields are
    # the row's name, and the corresponding coefficient
    rname = fields[3]  # row name
    if !haskey(row2idx, rname)
        # current row not in list of rows, current line ignored
        warn("COL ERROR UNKNOWN ROW $(rname)|$(ln)")
        return nothing
    end
    coeff = parse(Float64, fields[4])  # coefficient value
    if row2idx[rname] == 0
        # objective coefficient
        push!(obj_col, col2idx[cname])
        push!(obj_val, coeff)
    else
        # constraint coefficient
        push!(coeffs_col, col2idx[cname])
        push!(coeffs_row, row2idx[rname])
        push!(coeffs_val, coeff)
    end

    # optional other fields
    if length(fields) >= 6
        rname = fields[5]  # row name
        if !haskey(row2idx, rname)
            # current row not in list of rows, current line ignored
            warn("COL ERROR UNKNOWN ROW $(rname)|$(ln)")
            return nothing
        end
        coeff = parse(Float64, fields[6])  # coefficient value
        if row2idx[rname] == 0
            # objective coefficient
            push!(obj_col, col2idx[cname])
            push!(obj_val, coeff)
        else
            # constraint coefficient
            push!(coeffs_col, col2idx[cname])
            push!(coeffs_row, row2idx[rname])
            push!(coeffs_val, coeff)
        end
    end

    return nothing
end

function parsemps_rhs!(ln, fields, d, row2idx, rhs_row, rhs_val)
    # parse a line of the RHS section
    # `ln` is the current line
    # `fields` contains the fields of that line

    # Check for too short lines
    if length(fields) < 4
        # input error, current line is ignored
        warn("RHS LINE TOO SHORT|$(ln)")
        return nothing
    end

    # First field is empty
    # Second field is RHS name (may be empty)
    rhsname = fields[2]
    if d["rhs"] == ""
        # first time RHS is read
        d["rhs"] = rhsname
    elseif d["rhs"] != rhsname
        # other RHS, current line is ignored
        return nothing
    end

    # parse line
    rname = fields[3]
    if !haskey(row2idx, rname)
        # current row not in list of rows, current line ignored
        warn("RHS ERROR UNKNOWN ROW $(rname)|$(ln)")
        return nothing
    end
    rval = parse(Float64, fields[4])
    # update index and value
    if row2idx[rname] == 0
        d["obj_offset"] = -rval
    else
        push!(rhs_row, row2idx[rname])
        push!(rhs_val, rval)
    end

    # optional fields
    if length(fields) >= 6
        rname = fields[5]  # row name
        if !haskey(row2idx, rname)
            # current row not in list of rows, current line ignored
            warn("RHS ERROR UNKNOWN ROW $(rname)|$(ln)")
            return nothing
        end
        rval = parse(Float64, fields[6])  # coefficient value
        # update index and value
        if row2idx[rname] == 0
            d["obj_offset"] = -rval
        else
            push!(rhs_row, row2idx[rname])
            push!(rhs_val, rval)
        end
    end

    return nothing
end

function parsemps_ranges!(ln, fields, d, row2idx, ranges_row, ranges_val)
    if length(fields) < 4
        # input error, current line is ignored
        warn("RNG LINE TOO SHORT|$(ln)")
        return nothing
    end

    rngname = fields[2]
    if d["range"] == ""
        d["range"] = rngname
    elseif d["range"] != rngname
        # other range vector, ignore
        return nothing
    end

    # parse line
    rname = fields[3]
    rngval = parse(Float64, fields[4])
    if !haskey(row2idx, rname)
        # unknown row
        warn("RNG ERROR UNKNOWN ROW $(rname)|$(ln)")
        return nothing
    end
    push!(ranges_row, row2idx[rname])
    push!(ranges_val, rngval)

    if length(fields) >=6
        rname = fields[5]
        rngval = parse(Float64, fields[6])
        if !haskey(row2idx, rname)
            # unknown row
            warn("RNG ERROR UNKNOWN ROW $(rname)|$(ln)")
            return nothing
        end
        push!(ranges_row, row2idx[rname])
        push!(ranges_val, rngval)
    end

    return nothing
end

function parsemps_bounds!(ln, fields, d, col2idx, lb_col, lb_val, ub_col, ub_val)
    # parse a line of the BOUNDS section

    if (length(fields) < 3
        || ((length(fields) < 4)
            && (fields[1] == "LO" || fields[1] == "UP" || fields[1] == "FX"))
    )
        warn("BND ERROR LINE TOO SHORT|$(ln)")
        return nothing
    end

    # first field should be bound type
    btype = fields[1]
    bname = fields[2]
    cname = fields[3]
    if !haskey(col2idx, cname)
        warn("BND ERROR UNKNOWN COL $(rname)|$(ln)")
        return nothing
    end

    if btype == "LO"
        push!(lb_col, col2idx[cname])
        push!(lb_val, parse(Float64, fields[4]))
    elseif btype == "UP"
        push!(ub_col, col2idx[cname])
        push!(ub_val, parse(Float64, fields[4]))
    elseif btype == "FR"
        push!(lb_col, col2idx[cname])
        push!(lb_val, -Inf)
        push!(ub_col, col2idx[cname])
        push!(ub_val, Inf)
    elseif btype == "FX"
        push!(lb_col, col2idx[cname])
        push!(lb_val, parse(Float64, fields[4]))
        push!(ub_col, col2idx[cname])
        push!(ub_val, parse(Float64, fields[4]))
    elseif btype == "MI"
        push!(lb_col, col2idx[cname])
        push!(lb_val, -Inf)
    elseif btype == "PL"
        push!(ub_col, col2idx[cname])
        push!(ub_val, Inf)
    else
        # error in bound type
        warn("BND ERROR WRONG BOUND TYPE $(btype)")
    end

    return nothing
end


text:
using Match

function solve_ship(input)
	dirs = [(1, 0), (0, 1), (-1, 0), (0, -1)]
	facing::Int64 = 1
	x = y = 0
	for instruction in input
		action = instruction[1]
		value::Int64 = parse(Int64, instruction[2:end])
		@match action begin
			'E' => begin
				x += value
			end
			'N' => begin
				y += value
			end
			'W' => begin
				x -= value
			end
			'S' => begin
				y -= value
			end
			'L' => begin
				facing = mod1(facing + value / 90, 4)
			end
			'R' => begin
				facing = mod1(facing - value / 90, 4)
			end
			'F' => begin
				x += dirs[facing][1] * value
				y += dirs[facing][2] * value
			end
		end
	end
	abs(x) + abs(y)
end

function solve_waypoint(input, waypoint = [10 1])
	rotations_l = [
		[ 0  1; -1  0],
		[-1  0;  0 -1],
		[ 0 -1;  1  0],
		[ 1  0;  0  1]
	]
	rotations_r = [
		[ 0 -1;  1  0],
		[-1  0;  0 -1],
		[ 0  1; -1  0],
		[ 1  0;  0  1]
	]
	x = y = 0
	waypoint_pos = copy(waypoint)
	for instruction in input
		action = instruction[1]
		value::Int64 = parse(Int64, instruction[2:end])
		@match action begin
			'E' => begin
				waypoint_pos[1] += value
			end
			'N' => begin
				waypoint_pos[2] += value
			end
			'W' => begin
				waypoint_pos[1] -= value
			end
			'S' => begin
				waypoint_pos[2] -= value
			end
			'L' => begin
				waypoint_pos *= rotations_l[mod1(value ÷ 90, 4)]
			end
			'R' => begin
				waypoint_pos *= rotations_r[mod1(value ÷ 90, 4)]
			end
			'F' => begin
				x += waypoint_pos[1] * value
				y += waypoint_pos[2] * value
			end
		end
	end
	abs(x) + abs(y)
end

@assert solve_ship(eachline("input/test/test_day12.txt")) == 25 ["test case 1 failed"]
@assert solve_waypoint(eachline("input/test/test_day12.txt")) == 286 ["test case 2 failed"]

println("First star: " * string(solve_ship(eachline("input/day12.txt"))))
println("Second star: " * string(solve_waypoint(eachline("input/day12.txt"))))


text:
using LinearAlgebra
P = [   0   1   0   0   0;
        1/4 0   1/4 1/4 1/4;
        0   1/2 0   0   1/2;
        0   1/2 0   0   1/2;
        0   0   0   0   1]
T = P[1:4,1:4]
p0 = [1 0 0 0]
for n in 1:10
    println(first(p0*sum([T^k for k in 0:n])*ones(4)))
end
println("Using inverse: ", first(p0*inv(I-T)*ones(4)))
println("Eigenvalues of T: ", sort(eigvals(T)))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
function create_data_for_test()
    data = zeros(Float32,2,1000)
    data[:,1:250] .= [-1,-1]
    data[:,251:500] .= [-1,1]
    data[:,501:750] .= [1,-1]
    data[:,751:1000] .= [1,1]
    return data
end

@testset "Testing Module (Determinstic)" begin
    data = create_data_for_test()
    labels,clusters,weights = fit(data,100.0, iters = 200, seed = 12345, burnout = 5)
    @test all(data[:,1:250] .== data[:,1])
    @test all(data[:,251:500] .== data[:,251])
    @test all(data[:,501:750] .== data[:,501])
    @test all(data[:,751:1000] .== data[:,751])
    @test length(clusters) == 4
    println(weights)
    @test all(weights .>= 0.15)
    labels_histogram = get_labels_histogram(labels)
    for (k,v) in labels_histogram
        @test v == 250
    end
end

@testset "Testing Module (Random mess)" begin
    @everywhere Random.seed!(12345)
    x,labels,clusters = generate_gaussian_data(10^5,3,10,100.0)

    hyper_params = DPMMSubClusters.niw_hyperparams(Float32(1.0),
               zeros(Float32,3),
               Float32(5),
               Matrix{Float32}(I, 3, 3)*1)

    dp = dp_parallel(x,hyper_params,Float32(1000000000000000000000.0), 100,1,nothing,true,false,15,labels)
    @test length(dp[1].group.local_clusters) > 1
end

@testset "Multinomial Module And save load" begin
    @everywhere Random.seed!(12345)
    x,labels,clusters = generate_mnmm_data(10^3,100,20,50)
    @test size(x,1) == 100
    @test size(x,2) == 10^3
    npzwrite("save_load_test/mnm_data.npy",x')

    dp = dp_parallel("save_load_test/multinomial_params.jl")
    @test length(dp[1].group.local_clusters) > 1
    dp = run_model_from_checkpoint("save_load_test/checkpoint_20.jld2")
    @test length(dp[1].group.local_clusters) > 1
end


text:
@testset "Curve" begin

    curve = Curve([3, 4, 5])
    @test curve == Curve([1, 2, 3], [3, 4, 5])

    curve = Curve([1.0, 2.0], [3.0, 4.0])

    # test basic operations
    @test length(curve) == 2
    @test curve[1] == (1.0, 3.0)
    @test curve[end] == (2.0, 4.0)
    @test [p for p in curve] == [(1.0, 3.0), (2.0, 4.0)] # test iteration
    @test [curve; curve] == Curve([1.0, 2.0, 1.0, 2.0], [3.0, 4.0, 3.0, 4.0])
    @test curve + curve.y == Curve([1.0, 2.0], [6.0, 8.0])
    @test minimum(curve) == (1.0, 3.0)
    @test maximum(curve) == (2.0, 4.0)
    @test curve + 2.0 == Curve([1.0, 2.0], [5.0, 6.0])
    @test curve - 2.0 == Curve([1.0, 2.0], [1.0, 2.0])
    @test curve * 2.0 == Curve([1.0, 2.0], [6.0, 8.0])
    @test curve / 2.0 == Curve([1.0, 2.0], [1.5, 2.0])
    @test 2.0 + curve == Curve([1.0, 2.0], [5.0, 6.0])
    @test 2.0 - curve == Curve([1.0, 2.0], [-1.0, -2.0])
    @test 2.0 * curve == Curve([1.0, 2.0], [6.0, 8.0])
    @test 2.0 / curve == Curve([1.0, 2.0], [2.0/3.0, 0.5])

    # x and y need to have same length
    @test_throws ArgumentError Curve([1.0], [4.0, 4.0])
    @test_throws ArgumentError Curve([1.0, 2.0], [4.0])
end;

@testset "test sort" begin
    c = Curve([2, 6, 1], [4, 12, 2])
    @test sort(c) == Curve([1, 2, 6], [2, 4, 12])
    c = Curve(2:-0.5:0, 4:-0.5:2)
    @test sort(c) == Curve([0.0, 0.5, 1.0, 1.5, 2.0], [2.0, 2.5, 3.0, 3.5, 4.0])
    ns = NoiseSample(c)
    @test sort(ns) == NoiseSample([0.0, 0.5, 1.0, 1.5, 2.0], [-1.0, -0.5, 0.0, 0.5, 1.0])
    seed!(1)
    uc = c ± 0.01
    @test sort(uc).x == [0.0, 0.5, 1.0, 1.5, 2.0]
    for (yᵢ, y) in zip(sort(uc).y, collect(2:0.5:4))
        @test yᵢ ≈ y
    end
end

@testset "create Curves from collections" begin
    x = 1:10
    y = 11:20
    c = Curve(x, y)
    @test eltype(c) == Int
    @test c.x == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    @test c.y == [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]
    c = Curve(1:3, [1. , 2., 3.])
    @test eltype(c) == Float64
    @test c.x == [1., 2., 3.]
    @test c.y == [1., 2., 3.]
    c = Curve([1. , 2., 3.], 1:3)
    @test eltype(c) == Float64
    @test c.x == [1., 2., 3.]
    @test c.y == [1., 2., 3.]
    # element types must be subtypes of Number
    @test_throws MethodError Curve(["a", "c"], ["b", "d"])
end

@testset "crop()" begin
    c = Curve(1:10, 11:20)
    @test crop(c, 4, 6) == Curve([4, 5, 6], [14, 15, 16])
    c = Curve(1:10, collect(Float64, 11:20))
    @test crop(c, 3.1, 6.9) == Curve([4., 5., 6.], [14., 15., 16.])
    @test crop(c, 3.99999999999999, 6.00000000000001) == Curve([4., 5., 6.], [14., 15., 16.])
    @test crop(c, 3.1, 6.9) == Curve([4., 5., 6.], [14., 15., 16.])
    @test crop(c, 3.1, 7.0) == Curve([4., 5., 6., 7.], [14., 15., 16., 17.])
    @test crop(c, 3.0, 7.0) == Curve([3., 4., 5., 6., 7.], [13., 14., 15., 16., 17.])
end


text:
import .SpecialFunctions
import .SpecialFunctions: erf, erfc, besselj0, besselj1, bessely0, bessely1,
    besselj, bessely, gamma, lgamma

erf(x::Float128) =
    Float128(@ccall(libquadmath.erfq(x::Cfloat128)::Cfloat128))
erfc(x::Float128) =
    Float128(@ccall(libquadmath.erfcq(x::Cfloat128)::Cfloat128))

besselj0(x::Float128) =
    Float128(@ccall(libquadmath.j0q(x::Cfloat128)::Cfloat128))
besselj1(x::Float128) =
    Float128(@ccall(libquadmath.j1q(x::Cfloat128)::Cfloat128))

bessely0(x::Float128) =
    Float128(@ccall(libquadmath.y0q(x::Cfloat128)::Cfloat128))
bessely1(x::Float128) =
    Float128(@ccall(libquadmath.y1q(x::Cfloat128)::Cfloat128))

besselj(n::Integer, x::Float128) =
    Float128(@ccall(libquadmath.jnq(n::Cint, x::Cfloat128)::Cfloat128))
bessely(n::Integer, x::Float128) =
    Float128(@ccall(libquadmath.ynq(n::Cint, x::Cfloat128)::Cfloat128))

gamma(x::Float128) =
    Float128(@ccall(libquadmath.tgammaq(x::Cfloat128)::Cfloat128))
lgamma(x::Float128) =
    Float128(@ccall(libquadmath.lgammaq(x::Cfloat128)::Cfloat128))


text:
# function that display the results

function q_show_rep_a(account_a, account_b)
    gens, n_para = size(account_b)
    mean_output = account_a[1]
    mean_consume = account_a[2]
    mean_capital = account_a[3]
    mean_s_rate = account_a[4]

    println("mean_output: ", mean_output[end])
    println("mean_consume: ", mean_consume[end])
    println("mean_capital: ", mean_capital[end])
    println("mean_s_rate: ", mean_s_rate[end])

    periods = length(mean_output)
    rr = floor(periods / gens)

    waste_b = 1
    waste_a =Int((waste_b-1) * rr + 1)

    p1 = plot(waste_a:periods, [mean_output[waste_a:end] mean_consume[waste_a:end]], title ="output & consumption",  label ="")
    p2 = plot(waste_a:periods, mean_capital[waste_a:end], title ="capital", label ="")
    p3 = plot(waste_a:periods, mean_s_rate[waste_a:end], title ="savings rate", label ="")
    if n_para == 1
        p4 = plot(waste_b:gens, account_b[:,1][waste_b:end] , title ="Para a", label ="")
        plot(p3, p2, p1, p4, layout=(2,2))
    elseif n_para == 3
        p4 = plot(waste_b:gens, [account_b[:,1][waste_b:end] account_b[:,2][waste_b:end] account_b[:,3][waste_b:end]], title ="Para a,b", label ="")
        plot(p3, p2, p1, p4, layout=(2,2))
    elseif n_para == 4
        p4 = plot(waste_b:gens, [account_b[:,1][waste_b:end] account_b[:,2][waste_b:end] account_b[:,3][waste_b:end] account_b[:,4][waste_b:end]], title ="Para a,b", label ="")
        plot(p3, p2, p1, p4, layout=(2,2))
    else
        p4 = plot(waste_b:gens, [account_b[:,1][waste_b:end] account_b[:,2][waste_b:end]] , title ="Para a,b", label ="")
        p5 = plot(waste_b:gens, [account_b[:,3][waste_b:end] account_b[:,4][waste_b:end]] , title ="Para c,d", label ="")
        p6 = plot(waste_b:gens, [account_b[:,5][waste_b:end] account_b[:,6][waste_b:end]],  title ="Para e,f", label ="")
        plot(p3, p2, p1, p4, p5, p6, layout=(3,2))
    end
end

function q_show_sim_rep(account_a)
    mean_output = account_a[1]
    mean_consume = account_a[2]
    mean_capital = account_a[3]
    mean_s_rate = account_a[4]
    periods = length(mean_output)

    println("mean_output: ", mean_output[end])
    println("mean_consume: ", mean_consume[end])
    println("mean_capital: ", mean_capital[end])
    println("mean_s_rate: ", mean_s_rate[end])

    yMin_s = minimum(mean_s_rate) * .9
    yMax_s = maximum(mean_s_rate) * 1.1

    yMin_c = minimum(mean_capital) * .9
    yMax_c = maximum(mean_capital) * 1.1

    waste_a = 1
    p1 = plot(waste_a:periods, [mean_output[waste_a:end] mean_consume[waste_a:end]], title ="Output & Consumption",  label ="")
    p2 = plot(waste_a:periods, mean_capital[waste_a:end], title ="Capital", ylim=(yMin_c, yMax_c), label ="")
    p3 = plot(waste_a:periods, mean_s_rate[waste_a:end], title ="Savings Rate", ylim=(yMin_s, yMax_s), label ="")
    plot(p1, p2, p3,  layout=(3,1))
end

function q_show_eval_r_plot(e_para, decide_set)
    scale =  decide_set["scale"]

    x=range(0,stop=1,length=100)
    y=range(0,stop=1,length=100)

    println()
    function q_eval_a(capital,productivity)
        shock_p = [sim_set["t_shock"] sim_set["shock_duration"] 1 sim_set["periods"]]
        input = hcat(capital, productivity)
        ss = q_action_gen_one(input, e_para, scale, shock_p)
        return ss[1]
    end

    function q_e_display()
        a_f_eva = zeros(3,3)
        e_input = ones(1,2)
        for e_cap in 1:3
            e_input[1] = e_cap
            for e_prod in 1:3
                e_input[2] = e_prod
                shock_p = [sim_set["t_shock"] sim_set["shock_duration"] 1 sim_set["periods"]]
                dd = q_action_gen_one(e_input, e_para, scale, shock_p)
                a_f_eva[e_cap,e_prod]= dd[1]
            end
        end
    end

    q_e_display()
    display(plot(x,y,q_eval_a,st=:surface,camera=(-10,30)))
    return
end

function q_show_eval_rep(para, state, sim_set, decide_set)

    println()
    println("capital min mean max ", minimum(state[1]), "  ", mean(state[1]),"  ", maximum(state[1]), " std ", std(state[1]) )
    println("productivity min mean max ", minimum(state[2]), "  ", mean(state[2]),"  ", maximum(state[2]), " std ", std(state[2]) )
    println()
    println("para : ", para)
    shock_p = [sim_set["t_shock"] sim_set["shock_duration"] 1 sim_set["periods"]]
    ca = maximum( [mean(state[1]) - std(state[1]), 0.0])
    cs = std(state[1])
    pa = maximum( [mean(state[2]) - std(state[2]), 0.0])
    ps = std(state[2])
    dd= zeros(3,3)
    for i in 1:3
        for j in 1:3
            e_inputs = [(ca + (i-1)*cs) (pa + (j-1)*ps) ]
            qd = q_action_gen_one(e_inputs, para, decide_set["scale"], shock_p)
            dd[i,j] = qd[1]
        end
    end
    println("row: capital, col:productivity")
    show(stdout, "text/plain", dd)
    println()
end

function q_show_results(state, s_para, econ_set, decide_set, sim_set, account_econ, account_para)
    type =  sim_set["type"]
    action_type =  sim_set["action_type"]
    market_type =  sim_set["market_type"]
    shock_type =  sim_set["shock_type"]
    shock_duration =  sim_set["shock_duration"]
    scale =  decide_set["scale"]
    println()
    println("Simulation Complete")
    println("Type: ", type)
    println("Decsion_function_Type: ", action_type )
    println("sim_set: ", sim_set)
    println("decide_set: ", decide_set)
    println("econ_set: ", econ_set)

    # visualize decsion function
    e_para = account_para[end:end,:]  # last mean parameters
    println("Representative Agent")
    println("e_para ", e_para)
    println("show a")
    oo = q_show_eval_rep(e_para, state, sim_set, decide_set)
    q_show_eval_r_plot(e_para, decide_set)
    q_show_rep_a(account_econ, account_para)
end


text:
using DSGE, FileIO, JLD2, ModelConstructors, Test, Random, Dates, HDF5
path = dirname(@__FILE__)

generate_regime_switch_tests = false # Set to true if you want to regenerate the jld2 files for testing

if VERSION < v"1.5"
    ver = "111"
elseif VERSION < v"1.6"
    ver = "150"
else
    ver = "160"
end

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:cond_id, 0)
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:date_conditional_end, quartertodate("2015-Q4"))
m <= Setting(:use_population_forecast, true)

estroot = normpath(joinpath(dirname(@__FILE__), "..", "reference"))
overrides = forecast_input_file_overrides(m)
overrides[:mode] = joinpath(estroot, "optimize.h5")
overrides[:full] = joinpath(estroot, "metropolis_hastings.h5")

if haskey(ENV, "FRED_API_KEY") || isfile(joinpath(homedir(),".freddatarc"))
    df = load_data(m)
    skip_forecast_one_draw = false
else
    skip_forecast_one_draw = true
    @warn "Skipping forecast_one_draw tests because FRED_API_KEY not present"
end
# Make sure output_vars ignores the untransformed and 4Q things because they are
# computed in compute_meansbands
output_vars = add_requisite_output_vars([:histpseudo, :histobs, :histstdshocks,
                                         :histutpseudo, :histutobs,
                                         :hist4qpseudo, :hist4qobs,
                                         :forecaststates, :forecastpseudo, :forecastobs, :forecaststdshocks,
                                         :forecastutpseudo, :forecastutobs,
                                         :forecast4qpseudo, :forecast4qobs,
                                         :bddforecaststates, :bddforecastshocks, :bddforecastpseudo, :bddforecastobs,
                                         :shockdecpseudo, :shockdecobs,
                                         :trendstates, :trendobs, :trendpseudo,
                                         :dettrendstates, :dettrendobs, :dettrendpseudo,
                                         :irfstates, :irfpseudo, :irfobs])

# Check error handling for input_type = :subset
@testset "Ensure properly error handling for input_type = :subset" begin
    @test_throws ErrorException forecast_one(m, :subset, :none, output_vars,
                                             subset_inds = 1:10, forecast_string = "",
                                             verbose = :none)
    @test_throws ErrorException forecast_one(m, :subset, :none, output_vars,
                                             forecast_string = "test",
                                             verbose = :none)
end

# Run modal forecasts
out = Dict{Symbol, Dict{Symbol, Array{Float64}}}()
for cond_type in [:none, :semi, :full]
    forecast_one(m, :mode, cond_type, output_vars, verbose = :none)
    forecast_one(m, :init, cond_type, output_vars, verbose = :none)

    # Read output
    out[cond_type] = Dict{Symbol, Array{Float64}}()
    output_files = get_forecast_output_files(m, :mode, cond_type, output_vars)
    for var in keys(output_files)
        out[cond_type][var] = load(output_files[var], "arr")
    end
end

# Read expected output
exp_out = JLD2.jldopen("$path/../reference/forecast_one_out.jld2", "r") do file
    read(file, "exp_out")
end

# Test modal forecast outputs
specify_mode!(m, DSGE.get_forecast_input_file(m, :mode); verbose = :none)

@testset "Test modal forecast for all major output_vars" begin
    for cond_type in [:none, :semi, :full]
        # Histories
        @test @test_matrix_approx_eq exp_out[cond_type][:histpseudo]     out[cond_type][:histpseudo]

        # Forecasts
        @test @test_matrix_approx_eq exp_out[cond_type][:forecastobs]    out[cond_type][:forecastobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:forecastpseudo] out[cond_type][:forecastpseudo]

        # Shock decompositions, deterministic trends, trends
        @test @test_matrix_approx_eq exp_out[cond_type][:shockdecobs]    out[cond_type][:shockdecobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:shockdecpseudo] out[cond_type][:shockdecpseudo]
        @test @test_matrix_approx_eq exp_out[cond_type][:dettrendobs]    out[cond_type][:dettrendobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:dettrendpseudo] out[cond_type][:dettrendpseudo]
        @test @test_matrix_approx_eq exp_out[cond_type][:trendobs]       out[cond_type][:trendobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:trendpseudo]    out[cond_type][:trendpseudo]

        # IRFs
        @test @test_matrix_approx_eq exp_out[cond_type][:irfobs]         out[cond_type][:irfobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:irfpseudo]      out[cond_type][:irfpseudo]
    end
end

@testset "Test full-distribution forecasts run" begin
    m <= Setting(:forecast_block_size, 5)
    forecast_one(m, :prior, :none, output_vars, verbose = :none)
    for sampling_method in [:MH]
        m <= Setting(:sampling_method, sampling_method)
        for cond_type in [:none, :semi, :full]
            forecast_one(m, :full, cond_type, output_vars, verbose = :none)
            @test_throws ErrorException forecast_one(m, :subset, cond_type, output_vars, subset_inds = 1:10, verbose = :none)
            forecast_one(m, :subset, cond_type, output_vars, subset_inds = 1:10, forecast_string = "test", verbose = :none)
            forecast_one(m, :init_draw_shocks, cond_type, output_vars, verbose = :none)
            forecast_one(m, :mode_draw_shocks, cond_type, output_vars, verbose = :none)
        end
    end
end

# Test full-distribution blocking
m <= Setting(:forecast_block_size, 5)
forecast_one(m, :full, :none, output_vars, verbose = :none)
# Test read_forecast_output
@testset "Test full-distribution blocking" begin
    for input_type in [:mode, :full]
        output_files = get_forecast_output_files(m, input_type, :none, output_vars)
        @test ndims(DSGE.read_forecast_series(output_files[:trendobs], :trend, m.observables[:obs_gdp])) == 2
        @test ndims(DSGE.read_forecast_series(output_files[:forecastobs], :forecast, m.observables[:obs_gdp])) == 2
        @test ndims(DSGE.read_forecast_series(output_files[:irfobs], m.observables[:obs_gdp], m.exogenous_shocks[:rm_sh])) == 2
    end
end

@testset "Test forecast_one_draw" begin
    if !skip_forecast_one_draw
        m <= Setting(:regime_switching, true)
        for input_type in [:mode, :full]
            params = if input_type == :mode
                load_draws(m, input_type)
            else
                load_draws(m, input_type)[1, :]
            end
            setup_permanent_altpol!(m, AltPolicy(:historical, eqcond, solve))
            @test typeof(DSGE.forecast_one_draw(m, input_type, :none, output_vars, params, df;
                                                regime_switching = get_setting(m, :regime_switching),
                                                n_regimes = get_setting(m, :n_regimes))) == Dict{Symbol, Array{Float64}}

            # Test with alternative policy
            setup_permanent_altpol!(m, AltPolicy(:taylor93, eqcond, solve))
            @test typeof(DSGE.forecast_one_draw(m, input_type, :none, output_vars, params, df;
                                                regime_switching = get_setting(m, :regime_switching),
                                                n_regimes = get_setting(m, :n_regimes))) == Dict{Symbol, Array{Float64}}
        end
    end
end

## Now check for regime switching
# TODO: ADD TEST WHEN USING OVERRIDES W/REGIME-SWITCHING MODEL, also two more TODO below
custom_settings = Dict{Symbol, Setting}(
                                        :data_vintage             => Setting(:data_vintage, "160812"),
                                        :cond_vintage             => Setting(:cond_vintage, "160812"),
                                        :cond_id                  => Setting(:cond_id, 0),
                                        :use_population_forecast  => Setting(:use_population_forecast, true),
                                        :date_presample_start     => Setting(:date_presample_start, Date(1959, 9, 30)),
                                        :date_forecast_start      => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
                                        :date_conditional_end     => Setting(:date_conditional_end, DSGE.quartertodate("2016-Q3")),
                                        :n_mon_anticipated_shocks => Setting(:n_mon_anticipated_shocks, 6))
:forecast_horizons => Setting(:forecast_horizons, 12)
:impulse_response_horizons => Setting(:impulse_response_horizons, 12)

# Now check that regime switching works for different types of possible regimes
regime_dates_dicts = [Dict{Int, Date}(1 => DSGE.quartertodate("1959-Q3"),
                                      2 => DSGE.quartertodate("2010-Q1"),
                                      3 => DSGE.quartertodate("2012-Q4")),
                      Dict{Int, Date}(1 => DSGE.quartertodate("1959-Q3"),
                                      2 => DSGE.quartertodate("1980-Q2"),
                                      3 => DSGE.quartertodate("2012-Q4")),
                      Dict{Int, Date}(1 => DSGE.quartertodate("1959-Q3"),
                                      2 => DSGE.quartertodate("1980-Q2"),
                                      3 => DSGE.quartertodate("2003-Q4")),
                      Dict{Int, Date}(1 => DSGE.quartertodate("1959-Q3"),
                                      2 => DSGE.quartertodate("2008-Q4"),
                                      3 => DSGE.quartertodate("2012-Q4")),
                      Dict{Int, Date}(1 => DSGE.quartertodate("1959-Q3"),
                                      2 => DSGE.quartertodate("2000-Q2"),
                                      3 => DSGE.quartertodate("2008-Q4"))]

# Read in the expected outputs for all the regimes except when the first one,
# which we use to test for shockdec and irf outputs. Those cases are handled separately
exp_out_dict = JLD2.jldopen("$path/../reference/forecast_one_out_rs1_version=" * ver * ".jld2", "r") do file
    read(file, "exp_out_regime_switch_cases")
end

m = Model1002("ss10", custom_settings = custom_settings, testing = true)
m <= Setting(:rate_expectations_source, :ois)
dfs = Dict()
dfs[:none] = load("$path/../reference/regime_switch_data.jld2", "none")
dfs[:semi] = load("$path/../reference/regime_switch_data.jld2", "semi")
dfs[:full] = load("$path/../reference/regime_switch_data.jld2", "full")

if generate_regime_switch_tests
    exp_out_dict_new = exp_out_dict # We won't be testing, but we want to have the same structure as the existing one dict
end

@testset "Test modal and full distribution forecasts with regime switching for all major output_vars" begin
    # Loop over different times at which the regime switches, e.g. does the regime before or after ZLB
    m.settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
    m.test_settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
    m.settings[:forecast_ndraws] = Setting(:forecast_ndraws, 4)
    m.test_settings[:forecast_ndraws] = Setting(:forecast_ndraws, 4)
    for (k, v) in enumerate(regime_dates_dicts)

        # pseudo regime switching (no values have second/third regimes, but we are still saying regime_switching = true)
        m_rs1 = Model1002("ss10", testing = true, custom_settings = custom_settings)
        m_rs1 <= Setting(:rate_expectations_source, :ois)
        m_rs1.settings[:regime_switching] = Setting(:regime_switching, true)
        m_rs1.settings[:n_regimes] = Setting(:n_regimes, 3)
        m_rs1.settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
        m_rs1.settings[:forecast_ndraws] = Setting(:forecast_ndraws, 4)
        m_rs1.settings[:regime_dates] = Setting(:regime_dates, v)
        m_rs1.settings[:time_varying_trends] = Setting(:time_varying_trends, false)
        m_rs1.test_settings[:regime_switching] = Setting(:regime_switching, true)
        m_rs1.test_settings[:n_regimes] = Setting(:n_regimes, 3)
        m_rs1.test_settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
        m_rs1.test_settings[:forecast_ndraws] = Setting(:forecast_ndraws, 4)
        m_rs1.test_settings[:regime_dates] = Setting(:regime_dates, v)
        m_rs1.test_settings[:time_varying_trends] = Setting(:time_varying_trends, false)
        m_rs1 = DSGE.setup_regime_switching_inds!(m_rs1)

        # pseudo regime switching (identical values for standard deviations)
        m_rs2 = Model1002("ss51v", testing = true, custom_settings = custom_settings)
        m_rs2 <= Setting(:rate_expectations_source, :ois)
        m_rs2.settings[:regime_switching] = Setting(:regime_switching, true)
        m_rs2.settings[:n_regimes] = Setting(:n_regimes, 3)
        m_rs2.settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
        m_rs2.settings[:forecast_ndraws] = Setting(:forecast_ndraws, 4)
        m_rs2.settings[:regime_dates] = Setting(:regime_dates, v)
        m_rs2.settings[:time_varying_trends] = Setting(:time_varying_trends, false)
        m_rs2.test_settings[:regime_switching] = Setting(:regime_switching, true)
        m_rs2.test_settings[:n_regimes] = Setting(:n_regimes, 3)
        m_rs2.test_settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
        m_rs2.test_settings[:forecast_ndraws] = Setting(:forecast_ndraws, 4)
        m_rs2.test_settings[:time_varying_trends] = Setting(:time_varying_trends, false)
        m_rs2.test_settings[:regime_dates] = Setting(:regime_dates, v)
        m_rs2 = DSGE.setup_regime_switching_inds!(m_rs2)

        # non-trivial regime switching: eqcond matrices different and QQ matrix is different
        # - Does not test for switching in ZZ or ZZ_pseudo
        m_rs3 = Model1002("ss51v", testing = true, custom_settings = custom_settings)
        m_rs3 <= Setting(:rate_expectations_source, :ois)
        m_rs3.settings[:regime_switching] = Setting(:regime_switching, true)
        m_rs3.settings[:n_regimes] = Setting(:n_regimes, 3)
        m_rs3.settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
        m_rs3.settings[:forecast_ndraws] = Setting(:forecast_ndraws, 4)
        m_rs3.settings[:regime_dates] = Setting(:regime_dates, v)
        m_rs3.settings[:time_varying_trends] = Setting(:time_varying_trends, false)
        m_rs3.test_settings[:regime_switching] = Setting(:regime_switching, true)
        m_rs3.test_settings[:n_regimes] = Setting(:n_regimes, 3)
        m_rs3.test_settings[:regime_switching_ndraws] = Setting(:regime_switching_ndraws, 4)
        m_rs3.test_settings[:forecast_ndraws] = Setting(:forecast_ndraws, 4)
        m_rs3.test_settings[:time_varying_trends] = Setting(:time_varying_trends, false)
        m_rs3.test_settings[:regime_dates] = Setting(:regime_dates, v)
        if haskey(m_rs3.settings, :model2para_regimes)
            delete!(m_rs3.settings, :model2para_regimes)
        end
        if haskey(m_rs3.test_settings, :model2para_regimes)
            delete!(m_rs3.test_settings, :model2para_regimes)
        end
        m_rs3 = DSGE.setup_regime_switching_inds!(m_rs3)

        # Need to set shocks for second and third regimes
        global prop = 1.   # Use the same value as the first regime
        global prop3 = .95 # Set the values to .95 of the first regime values
        for i in 1:3
            if i == 1
                oldprop3 = prop3
                prop3 = prop
            end
            ModelConstructors.set_regime_val!(m_rs2[:α], i, prop * m[:α].value)
            ModelConstructors.set_regime_val!(m_rs2[:σ_g], i, prop * m[:σ_g].value)
            ModelConstructors.set_regime_val!(m_rs2[:σ_b], i, prop * m[:σ_b].value)
            ModelConstructors.set_regime_val!(m_rs2[:σ_μ], i, prop * m[:σ_μ].value)
            ModelConstructors.set_regime_val!(m_rs2[:σ_ztil], i, prop * m[:σ_ztil].value)
            ModelConstructors.set_regime_val!(m_rs2[:σ_λ_f], i, prop * m[:σ_λ_f].value)
            ModelConstructors.set_regime_val!(m_rs2[:σ_λ_w], i, prop * m[:σ_λ_w].value)
            ModelConstructors.set_regime_val!(m_rs2[:σ_r_m], i, prop * m[:σ_r_m].value)
            ModelConstructors.set_regime_val!(m_rs2[:σ_σ_ω], i, prop * m[:σ_σ_ω].value)
            ModelConstructors.set_regime_val!(m_rs2[:σ_μ_e], i, prop * m[:σ_μ_e].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs2[:σ_γ], i, prop * m[:σ_γ].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs2[:σ_π_star], i, prop * m[:σ_π_star].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs2[:σ_lr], i, prop * m[:σ_lr].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs2[:σ_z_p], i, prop * m[:σ_z_p].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs2[:σ_tfp], i, prop * m[:σ_tfp].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs2[:σ_gdpdef], i, prop * m[:σ_gdpdef].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs2[:σ_corepce], i, prop * m[:σ_corepce].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs2[:σ_gdp], i, prop * m[:σ_gdp].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs2[:σ_gdi], i, prop * m[:σ_gdi].value, override_bounds = true)

            for j = 1:DSGE.n_mon_anticipated_shocks(m_rs2)
                ModelConstructors.set_regime_val!(m_rs2[Symbol("σ_r_m$(j)")], i, prop * m[Symbol("σ_r_m$(i)")])
            end

            ModelConstructors.set_regime_val!(m_rs3[:α], i, prop3 * m[:α].value) # Change α to change eqcond matrices
            ModelConstructors.set_regime_val!(m_rs3[:σ_g], i, prop3 * m[:σ_g].value)
            ModelConstructors.set_regime_val!(m_rs3[:σ_b], i, prop3 * m[:σ_b].value)
            ModelConstructors.set_regime_val!(m_rs3[:σ_μ], i, prop3 * m[:σ_μ].value)
            ModelConstructors.set_regime_val!(m_rs3[:σ_ztil], i, prop3 * m[:σ_ztil].value)
            ModelConstructors.set_regime_val!(m_rs3[:σ_λ_f], i, prop3 * m[:σ_λ_f].value)
            ModelConstructors.set_regime_val!(m_rs3[:σ_λ_w], i, prop3 * m[:σ_λ_w].value)
            ModelConstructors.set_regime_val!(m_rs3[:σ_r_m], i, prop3 * m[:σ_r_m].value)
            ModelConstructors.set_regime_val!(m_rs3[:σ_σ_ω], i, prop3 * m[:σ_σ_ω].value)
            ModelConstructors.set_regime_val!(m_rs3[:σ_μ_e], i, prop3 * m[:σ_μ_e].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs3[:σ_γ], i, prop3 * m[:σ_γ].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs3[:σ_π_star], i, prop3 * m[:σ_π_star].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs3[:σ_lr], i, prop3 * m[:σ_lr].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs3[:σ_z_p], i, prop3 * m[:σ_z_p].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs3[:σ_tfp], i, prop3 * m[:σ_tfp].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs3[:σ_gdpdef], i, prop3 * m[:σ_gdpdef].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs3[:σ_corepce], i, prop3 * m[:σ_corepce].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs3[:σ_gdp], i, prop3 * m[:σ_gdp].value, override_bounds = true)
            ModelConstructors.set_regime_val!(m_rs3[:σ_gdi], i, prop3 * m[:σ_gdi].value, override_bounds = true)

            for j = 1:DSGE.n_mon_anticipated_shocks(m_rs3)
                ModelConstructors.set_regime_val!(m_rs3[Symbol("σ_r_m$(j)")], i, prop3 * m[Symbol("σ_r_m$(i)")])
            end

            if i == 1
                global prop3 = oldprop3
            end
        end

        # Set up output vars, w/shockdecs and irfs tested only once b/c they take up a lot of space
        if k == 1
            output_vars = add_requisite_output_vars([:histpseudo, :histobs, :histstdshocks,
                                                     :histutpseudo, :histutobs,
                                                     :hist4qpseudo, :hist4qobs,
                                                     :forecaststates, :forecastpseudo, :forecastobs, :forecaststdshocks,
                                                     :forecastutpseudo, :forecastutobs,
                                                     :forecast4qpseudo, :forecast4qobs,
                                                     :bddforecaststates, :bddforecastshocks, :bddforecastpseudo, :bddforecastobs,
                                                     :shockdecpseudo, :shockdecobs,
                                                     :trendstates, :trendobs, :trendpseudo,
                                                     :dettrendstates, :dettrendobs, :dettrendpseudo,
                                                     :irfstates, :irfpseudo, :irfobs])
        else
            output_vars = add_requisite_output_vars([:histpseudo, :histobs, :histstdshocks,
                                                     :histutpseudo, :histutobs,
                                                     :hist4qpseudo, :hist4qobs,
                                                     :forecaststates, :forecastpseudo, :forecastobs, :forecaststdshocks,
                                                     :forecastutpseudo, :forecastutobs,
                                                     :forecast4qpseudo, :forecast4qobs,
                                                     :bddforecaststates, :bddforecastshocks, :bddforecastpseudo, :bddforecastobs])
        end

        # Save output in these dictionaries for different cases
        out = Dict{Symbol, Dict{Symbol, Array{Float64}}}()
        out_rs1 = Dict{Symbol, Dict{Symbol, Array{Float64}}}()
        out_rs2 = Dict{Symbol, Dict{Symbol, Array{Float64}}}()
        out_rs3 = Dict{Symbol, Dict{Symbol, Array{Float64}}}()
        params = map(x -> x.value, m.parameters)
        params_rs1 = map(x -> x.value, m_rs1.parameters)
        params_rs2 = map(x -> x.value, m_rs2.parameters)
        for para in m_rs3.parameters
            ModelConstructors.toggle_regime!(para, 3)
        end
        params_rs3 = map(x -> x.value, m_rs3.parameters)

        # Run modal forecasts w/different conditioning
        for cond_type in [:none, :semi, :full]
            setup_regime_switching_inds!(m_rs1; cond_type = cond_type)
            setup_regime_switching_inds!(m_rs2; cond_type = cond_type)
            setup_regime_switching_inds!(m_rs3; cond_type = cond_type)
            forecast_one(m, :mode, cond_type, output_vars, verbose = :none, params = params, df = dfs[cond_type])
            forecast_one(m, :init, cond_type, output_vars, verbose = :none, params = params, df = dfs[cond_type])
            forecast_one(m_rs1, :mode, cond_type, output_vars, verbose = :none, params = params_rs1, df = dfs[cond_type])
            forecast_one(m_rs1, :init, cond_type, output_vars, verbose = :none, params = params_rs1, df = dfs[cond_type])
            forecast_one(m_rs2, :mode, cond_type, output_vars, verbose = :none, params = params_rs2, df = dfs[cond_type])
            forecast_one(m_rs2, :init, cond_type, output_vars, verbose = :none, params = params_rs2, df = dfs[cond_type])
            forecast_one(m_rs3, :mode, cond_type, output_vars, verbose = :none, params = params_rs3, df = dfs[cond_type])
            forecast_one(m_rs3, :init, cond_type, output_vars, verbose = :none, params = params_rs3, df = dfs[cond_type])

            # Read output
            out[cond_type] = Dict{Symbol, Array{Float64}}()
            out_rs1[cond_type] = Dict{Symbol, Array{Float64}}()
            out_rs2[cond_type] = Dict{Symbol, Array{Float64}}()
            out_rs3[cond_type] = Dict{Symbol, Array{Float64}}()
            output_files = get_forecast_output_files(m, :mode, cond_type, output_vars)
            output_files_rs1 = get_forecast_output_files(m_rs1, :mode, cond_type, output_vars)
            output_files_rs2 = get_forecast_output_files(m_rs2, :mode, cond_type, output_vars)
            output_files_rs3 = get_forecast_output_files(m_rs3, :mode, cond_type, output_vars)
            for var in keys(output_files)
                if !(var in [:shockdecpseudo, :shockdecobs, :irfstates, :irfobs, :irfpseudo]) ||
                    (k == 1)
                    out[cond_type][var] = load(output_files[var], "arr")
                    out_rs1[cond_type][var] = load(output_files_rs1[var], "arr")
                    out_rs2[cond_type][var] = load(output_files_rs2[var], "arr")
                    out_rs3[cond_type][var] = load(output_files_rs3[var], "arr")
                end
            end
        end

        # Check parameters actually switched and can switch back
        for para in m_rs3.parameters
            ModelConstructors.toggle_regime!(para, 3)
        end
        @test m_rs3[:α].value != m[:α].value # check the regimes do not match
        for para in m_rs3.parameters
            ModelConstructors.toggle_regime!(para, 1)
        end
        @test m_rs3[:α].value == m[:α].value # check the regimes match after toggling

        if generate_regime_switch_tests
            if k == 1
                exp_out_regime_switch_new = out
                exp_out_true_regime_switch_new = out_rs3
            else
                exp_out_dict_new[k][:out] = out
                exp_out_dict_new[k][:out_rs3] = out_rs3
            end
        else
            # Read expected output
            if k == 1
                exp_out, exp_out_true = JLD2.jldopen("$path/../reference/forecast_one_out_rs2_version=" * ver * ".jld2", "r") do file
                    read(file, "exp_out_regime_switch"), read(file, "exp_out_true_regime_switch")
                end
            else
                exp_out      = exp_out_dict[k][:out]
                exp_out_true = exp_out_dict[k][:out_rs3]
            end
            for cond_type in [:none]#, :semi, :full]
                # Histories
                @test @test_matrix_approx_eq exp_out[cond_type][:histpseudo]          out[cond_type][:histpseudo]
                @test @test_matrix_approx_eq exp_out[cond_type][:histpseudo]          out_rs1[cond_type][:histpseudo]
                @test @test_matrix_approx_eq exp_out[cond_type][:histpseudo]          out_rs2[cond_type][:histpseudo]
                @test @test_matrix_approx_eq exp_out_true[cond_type][:histpseudo]     out_rs3[cond_type][:histpseudo]
                @test !(exp_out[cond_type][:histpseudo] ≈                             out_rs3[cond_type][:histpseudo])

                # Forecasts
                @test @test_matrix_approx_eq exp_out[cond_type][:forecastobs]             out[cond_type][:forecastobs]
                @test @test_matrix_approx_eq exp_out[cond_type][:forecastpseudo]          out[cond_type][:forecastpseudo]

                @test @test_matrix_approx_eq exp_out[cond_type][:forecastobs]             out_rs1[cond_type][:forecastobs]
                @test @test_matrix_approx_eq exp_out[cond_type][:forecastobs]             out_rs2[cond_type][:forecastobs]
                @test @test_matrix_approx_eq exp_out_true[cond_type][:forecastobs]        out_rs3[cond_type][:forecastobs]
                @test @test_matrix_approx_eq exp_out[cond_type][:forecastpseudo]          out_rs1[cond_type][:forecastpseudo]
                @test @test_matrix_approx_eq exp_out[cond_type][:forecastpseudo]          out_rs2[cond_type][:forecastpseudo]
                @test @test_matrix_approx_eq exp_out_true[cond_type][:forecastpseudo]     out_rs3[cond_type][:forecastpseudo]
                @test !(exp_out[cond_type][:forecastobs] ≈                                out_rs3[cond_type][:forecastobs])
                @test !(exp_out[cond_type][:forecastpseudo] ≈                             out_rs3[cond_type][:forecastpseudo])

                # B/c memory limits, only testing shock decs w/first test case
                if k == 1
                    # Shock decompositions, deterministic trends, trends
                    @test @test_matrix_approx_eq exp_out[cond_type][:shockdecobs]    out[cond_type][:shockdecobs]
                    @test @test_matrix_approx_eq exp_out[cond_type][:shockdecpseudo] out[cond_type][:shockdecpseudo]
                    @test @test_matrix_approx_eq exp_out[cond_type][:dettrendobs]    out[cond_type][:dettrendobs]
                    @test @test_matrix_approx_eq exp_out[cond_type][:dettrendpseudo] out[cond_type][:dettrendpseudo]
                    @test @test_matrix_approx_eq exp_out[cond_type][:trendobs]       out[cond_type][:trendobs]
                    @test @test_matrix_approx_eq exp_out[cond_type][:trendpseudo]    out[cond_type][:trendpseudo]

                    @test @test_matrix_approx_eq exp_out[cond_type][:shockdecobs]    out_rs1[cond_type][:shockdecobs]
                    @test @test_matrix_approx_eq exp_out[cond_type][:shockdecpseudo] out_rs1[cond_type][:shockdecpseudo]
                    @test @test_matrix_approx_eq exp_out[cond_type][:dettrendobs]    out_rs1[cond_type][:dettrendobs]
                    @test @test_matrix_approx_eq exp_out[cond_type][:dettrendpseudo] out_rs1[cond_type][:dettrendpseudo]
                    @test @test_matrix_approx_eq exp_out[cond_type][:trendobs]       out_rs1[cond_type][:trendobs][:, 1]
                    @test @test_matrix_approx_eq exp_out[cond_type][:trendpseudo]    out_rs1[cond_type][:trendpseudo][:, 1]
                    @test @test_matrix_approx_eq exp_out[cond_type][:trendobs]       out_rs1[cond_type][:trendobs][:, 2]
                    @test @test_matrix_approx_eq exp_out[cond_type][:trendpseudo]    out_rs1[cond_type][:trendpseudo][:, 2]
                    @test @test_matrix_approx_eq exp_out[cond_type][:trendobs]       out_rs1[cond_type][:trendobs][:, 3]
                    @test @test_matrix_approx_eq exp_out[cond_type][:trendpseudo]    out_rs1[cond_type][:trendpseudo][:, 3]

                    @test @test_matrix_approx_eq exp_out_true[cond_type][:shockdecobs]    out_rs3[cond_type][:shockdecobs]
                    @test @test_matrix_approx_eq exp_out_true[cond_type][:shockdecpseudo] out_rs3[cond_type][:shockdecpseudo]
                    @test @test_matrix_approx_eq exp_out_true[cond_type][:dettrendobs]    out_rs3[cond_type][:dettrendobs]
                    @test @test_matrix_approx_eq exp_out_true[cond_type][:dettrendpseudo] out_rs3[cond_type][:dettrendpseudo]
                    @test @test_matrix_approx_eq exp_out_true[cond_type][:trendobs]       out_rs3[cond_type][:trendobs]
                    @test @test_matrix_approx_eq exp_out_true[cond_type][:trendpseudo]    out_rs3[cond_type][:trendpseudo]

                    # IRFs
                    @test @test_matrix_approx_eq exp_out[cond_type][:irfobs]              out[cond_type][:irfobs]
                    @test @test_matrix_approx_eq exp_out[cond_type][:irfpseudo]           out[cond_type][:irfpseudo]
                    @test @test_matrix_approx_eq exp_out[cond_type][:irfobs]              out_rs1[cond_type][:irfobs]
                    @test @test_matrix_approx_eq exp_out[cond_type][:irfpseudo]           out_rs1[cond_type][:irfpseudo]
                    @test @test_matrix_approx_eq exp_out_true[cond_type][:irfobs]         out_rs3[cond_type][:irfobs]
                    @test @test_matrix_approx_eq exp_out_true[cond_type][:irfpseudo]      out_rs3[cond_type][:irfpseudo]
                end
            end
        end

        if k == 1 # only testing full distribution with the first case of regime switching

            # Construct fake matrix of parameters draws
            params_rs1 = repeat(params_rs1', get_setting(m_rs1, :regime_switching_ndraws))
            params_rs2 = repeat(params_rs2', get_setting(m_rs2, :regime_switching_ndraws))
            params_rs3 = repeat(params_rs3', get_setting(m_rs3, :regime_switching_ndraws))
            m <= Setting(:forecast_block_size, 2)
            m_rs1 <= Setting(:forecast_block_size, 2)
            m_rs2 <= Setting(:forecast_block_size, 2)
            m_rs3 <= Setting(:forecast_block_size, 2)
            m <= Setting(:forecast_jstep, 1)
            m_rs1 <= Setting(:forecast_jstep, 1)
            m_rs2 <= Setting(:forecast_jstep, 1)
            m_rs3 <= Setting(:forecast_jstep, 1)

            # @testset "Test full-distribution forecasts run" begin
            # TODO: UNCOMMENT THE FOLLOWING LINES WHEN LOADING DRAWS FROM A SAVED ESTIMATION HAS BEEN IMPLEMENTED
            # for sampling_method in [:MH]
            #     m <= Setting(:sampling_method, sampling_method)
            #     for cond_type in [:none, :semi, :full]
            #         forecast_one(m, :full, cond_type, output_vars, verbose = :none)
            #         @test_throws ErrorException forecast_one(m, :subset, cond_type, output_vars, subset_inds = 1:10, verbose = :none)
            #         forecast_one(m, :subset, cond_type, output_vars, subset_inds = 1:10, forecast_string = "test", verbose = :none)
            #         forecast_one(m, :init_draw_shocks, cond_type, output_vars, verbose = :none)
            #         forecast_one(m, :mode_draw_shocks, cond_type, output_vars, verbose = :none)
            #     end
            # end
            # end

            out_rs1 = Dict{Symbol, Dict{Symbol, Dict{Symbol, Array{Float64}}}}()
            out_rs2 = Dict{Symbol, Dict{Symbol, Dict{Symbol, Array{Float64}}}}()
            out_rs3 = Dict{Symbol, Dict{Symbol, Dict{Symbol, Array{Float64}}}}()

            for cond_type in [:none, :semi, :full]
                # out[cond_type] = Dict{Symbol, Dict{Symbol, Array{Float64}}}()
                out_rs1[cond_type] = Dict{Symbol, Dict{Symbol, Array{Float64}}}()
                out_rs2[cond_type] = Dict{Symbol, Dict{Symbol, Array{Float64}}}()
                out_rs3[cond_type] = Dict{Symbol, Dict{Symbol, Array{Float64}}}()

                # Forecast and read output from forecast
                for fcast_type in [:full, :subset, :init_draw_shocks, :mode_draw_shocks, :prior]
                    if fcast_type == :subset
                        for (model, para_rs) in zip([m_rs1, m_rs2, m_rs3], [params_rs1, params_rs2, params_rs3])
                            Random.seed!(1793)
                            forecast_one(model, :subset, cond_type, output_vars, subset_inds = 1:2, forecast_string = "test",
                                         verbose = :none, params = para_rs, df = dfs[cond_type])
                            @test_throws ErrorException forecast_one(model, :subset, cond_type, output_vars, subset_inds = 1:2,
                                                                     verbose = :none, params = para_rs, df = dfs[cond_type])
                        end
                    elseif cond_type == :none && fcast_type == :prior
                        # TODO: IMPLEMENT DRAWING FROM PRIOR FOR FORECAST WITH REGIME SWITCHING
                        # forecast_one(m_rs1, :prior, :none, output_vars, verbose = :none)
                        continue
                    elseif fcast_type == :prior
                        continue
                    elseif fcast_type == :mode_draw_shocks
                        for (model, para_rs) in zip([m_rs1, m_rs2, m_rs3], [params_rs1, params_rs2, params_rs3])
                            Random.seed!(1793)
                            forecast_one(model, fcast_type, cond_type, output_vars, verbose = :none, params = para_rs[1, :], df = dfs[cond_type])
                        end
                    else
                        for (model, para_rs) in zip([m_rs1, m_rs2, m_rs3], [params_rs1, params_rs2, params_rs3])
                            Random.seed!(1793)
                            forecast_one(model, fcast_type, cond_type, output_vars, verbose = :none, params = para_rs, df = dfs[cond_type])
                        end
                    end
                    out_rs1[cond_type][fcast_type] = Dict{Symbol, Array{Float64}}()
                    out_rs2[cond_type][fcast_type] = Dict{Symbol, Array{Float64}}()
                    out_rs3[cond_type][fcast_type] = Dict{Symbol, Array{Float64}}()
                    if fcast_type == :subset
                        output_files_rs1 = get_forecast_output_files(m_rs1, fcast_type, cond_type, output_vars; forecast_string = "test")
                        output_files_rs2 = get_forecast_output_files(m_rs2, fcast_type, cond_type, output_vars; forecast_string = "test")
                        output_files_rs3 = get_forecast_output_files(m_rs3, fcast_type, cond_type, output_vars; forecast_string = "test")
                    else
                        output_files_rs1 = get_forecast_output_files(m_rs1, fcast_type, cond_type, output_vars)
                        output_files_rs2 = get_forecast_output_files(m_rs2, fcast_type, cond_type, output_vars)
                        output_files_rs3 = get_forecast_output_files(m_rs3, fcast_type, cond_type, output_vars)
                    end
                    for var in keys(output_files_rs1)
                        if fcast_type != :prior || (cond_type == :none)
                            out_rs1[cond_type][fcast_type][var] = load(output_files_rs1[var], "arr")
                        end
                        if fcast_type != :prior
                            out_rs2[cond_type][fcast_type][var] = load(output_files_rs2[var], "arr")
                            out_rs3[cond_type][fcast_type][var] = load(output_files_rs3[var], "arr")
                        end
                    end
                end
            end

            if generate_regime_switch_tests
                exp_out_regime_switch_full_new = out_rs1
                exp_out_true_regime_switch_full_new = out_rs3

                JLD2.jldopen("$path/../reference/forecast_one_out_rs2_version=" * ver * ".jld2", true, true, true, IOStream) do file
                    # Now for the new additions
                    write(file, "exp_out_regime_switch", exp_out_regime_switch_new)
                    write(file, "exp_out_true_regime_switch", exp_out_true_regime_switch_new)
                    write(file, "exp_out_regime_switch_full", exp_out_regime_switch_full_new)
                end

                JLD2.jldopen("$path/../reference/forecast_one_out_rs3_version=" * ver * ".jld2", true, true, true, IOStream) do file
                    write(file, "exp_out_true_regime_switch_full", exp_out_true_regime_switch_full_new)
                end
            else
                exp_out = JLD2.jldopen("$path/../reference/forecast_one_out_rs2_version=" * ver * ".jld2", "r") do file
                    read(file, "exp_out_regime_switch_full")
                end
                exp_out_true = JLD2.jldopen("$path/../reference/forecast_one_out_rs3_version=" * ver * ".jld2", "r") do file
                    read(file, "exp_out_true_regime_switch_full")
                end

                for cond_type in [:none, :semi, :full]
                    for fcast_type in [:full, :subset, :init_draw_shocks, :mode_draw_shocks] # TODO: ADD PRIOR
                        # Histories
                        @test maximum(abs.(exp_out[cond_type][fcast_type][:histpseudo] - out_rs1[cond_type][fcast_type][:histpseudo])) < 1e-5
                        @test maximum(abs.(exp_out[cond_type][fcast_type][:histpseudo] - out_rs2[cond_type][fcast_type][:histpseudo])) < 1e-5
                        @test maximum(abs.(exp_out_true[cond_type][fcast_type][:histpseudo] - out_rs3[cond_type][fcast_type][:histpseudo])) < 1e-5
                        @test !(exp_out[cond_type][fcast_type][:histpseudo] ≈                             out_rs3[cond_type][fcast_type][:histpseudo])
                        # Forecasts
                        @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:forecastobs]             out_rs1[cond_type][fcast_type][:forecastobs]
                        @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:forecastobs]             out_rs2[cond_type][fcast_type][:forecastobs]
                        @test @test_matrix_approx_eq exp_out_true[cond_type][fcast_type][:forecastobs][:, vcat(1:9, 13), 1]        out_rs3[cond_type][fcast_type][:forecastobs][:, vcat(1:9, 13), 1]
                        @test @test_matrix_approx_eq exp_out_true[cond_type][fcast_type][:forecastobs][:, :, 2]        out_rs3[cond_type][fcast_type][:forecastobs][:, :, 2]
                        @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:forecastpseudo]          out_rs1[cond_type][fcast_type][:forecastpseudo]
                        @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:forecastpseudo]          out_rs2[cond_type][fcast_type][:forecastpseudo]
                        @test @test_matrix_approx_eq exp_out_true[cond_type][fcast_type][:forecastpseudo]     out_rs3[cond_type][fcast_type][:forecastpseudo]
                        @test !(exp_out[cond_type][fcast_type][:forecastobs] ≈                                out_rs3[cond_type][fcast_type][:forecastobs])
                        @test !(exp_out[cond_type][fcast_type][:forecastpseudo] ≈                             out_rs3[cond_type][fcast_type][:forecastpseudo])

                        # We do not compare distributional shocks decs or IRFs b/c memory limits, although we do check that the methods run
                        if k == 1
                            # # Shock decompositions, deterministic trends, trends
                            # @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:shockdecobs]    out[cond_type][fcast_type][:shockdecobs]
                            # @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:shockdecpseudo] out[cond_type][fcast_type][:shockdecpseudo]
                            # @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:dettrendobs]    out[cond_type][fcast_type][:dettrendobs]
                            # @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:dettrendpseudo] out[cond_type][fcast_type][:dettrendpseudo]
                            # @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:trendobs]       out[cond_type][fcast_type][:trendobs]
                            # @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:trendpseudo]    out[cond_type][fcast_type][:trendpseudo]

                            # # IRFs
                            # @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:irfobs]         out[cond_type][fcast_type][:irfobs]
                            # @test @test_matrix_approx_eq exp_out[cond_type][fcast_type][:irfpseudo]      out[cond_type][fcast_type][:irfpseudo]
                        end
                    end
                end
            end
        end
    end

    if generate_regime_switch_tests
        JLD2.jldopen("$path/../reference/forecast_one_out_rs1_version=" * ver * ".jld2", true, true, true, IOStream) do file

            # Now for the new additions
            write(file, "exp_out_regime_switch_cases", exp_out_dict_new)
        end
    end

end

nothing


text:
#=
rabin_karp_string_search:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-09-20
=#

module RabinKarpSearch
    # parameters
    const cardinal = 101

    function get_hash(s, len)
        hash_value = 0
        pow_num = len
        for i in 1:len
            pow_num -= 1
            ascii_num = Int(s[i])
            hash_value += ascii_num * (cardinal^pow_num)
        end
        return hash_value
    end

    function rabin_karp_search(text, pattern)
        println("")
        println("Search $(pattern) from $(text)")

        len_txt = length(text)
        len_pat = length(pattern)

        # convert from string to hash value
        hash_txt = get_hash(text[1:len_pat], len_pat)
        hash_pat = get_hash(pattern[1:len_pat], len_pat)

        # search
        for i in 1:(len_txt)
            if hash_txt == hash_pat
                if text[i:i+len_pat-1] == pattern
                    println("Found at index:$(i)~$(i+len_pat-1)")
                    return i
                end
            end
            # prevent bounds error
            if (i + len_pat) > len_txt
                break
            end
            # calculate next hash value
            hash_txt = hash_txt - (Int(text[i]) * (cardinal^(len_pat-1)))
            hash_txt *= cardinal
            hash_txt = hash_txt + (Int(text[i+len_pat]) * (cardinal^0))
        end

        println("Pattern not found")
        return -1
    end

    function main()
        @time println(rabin_karp_search("Hello sunshine", "sun"))
        @time println(rabin_karp_search("GEEKS FOR GEEKS", "GEEKS"))
        @time println(rabin_karp_search("aaaaaaaaaaaaaaaaaaaaaaaaaaaa", "bbb"))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .RabinKarpSearch
    RabinKarpSearch.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
# functional definitions

import PDESolver.createFunctional

@doc """
###AdvectionEquationMod.QfluxData

Data type for storing relevant information pertaining to an a functional or an
objective function.

**Members**

*  `bcnums` : boundary condition groups on which the functional is to be
                             computed.
*  `val` : Computed value of the functional
*  `target_qFlux` : Target value for the functional qFlux

"""->
mutable struct QfluxData{Topt} <: AbstractBoundaryFunctional{Topt}
  bcnums::Array{Int,1}
  val::Topt
  target_qflux::Topt
end # End type OfluxData

"""
  Constructor for QfluxData
"""
function QfluxDataConstructor(::Type{Topt}, mesh, sbp, eqn, opts, bcnums) where Topt

  val = zero(Topt)
  target_qflux = zero(Topt)
  return QfluxData{Topt}(bcnums, val, target_qflux)
end

"""
  Functional that integrates the solution q over the specified boundary(/ies)

  **Fields**

   * bcnums: the boundary condition groups the functional is computed over
   * val: the value of the functional, initially 0.0

"""
mutable struct IntegralQData{Topt} <: AbstractBoundaryFunctional{Topt}
  bcnums::Array{Int, 1}
  val::Topt
end

"""
  Constructor for IntegralQData
"""
function IntegralQDataConstructor(::Type{Topt}, mesh, sbp, eqn, opts,
                                  bcnums) where Topt

  val = 0.0
  return IntegralQData{Topt}(bcnums, val)
end

"""
  See the generic method for docs.
"""
function createFunctional(mesh::AbstractMesh,
                  sbp::AbstractOperator,
                  eqn::AdvectionData{Tsol}, opts,
                  functional_name::AbstractString,
                  functional_bcs::Vector{I}) where {Tsol, I<:Integer}

  func_constructor = FunctionalDict[functional_name]
  functional = func_constructor(Tsol, mesh, sbp, eqn, opts, functional_bcs)

  return functional
end

"""
  Maps names of functionals to their outer constructor functions

  See `FunctionalDict` in the Euler module for descriptions of the constructors
  and arguments.
"""
global const FunctionalDict = Dict{String, Function}(
"qflux" => QfluxDataConstructor,
"integralq" => IntegralQDataConstructor,
)


text:
struct CompatiblePaths
    S::States
    C::Vector{ChanceNode}
    Z::DecisionStrategy
    fixed::Dict{Node, State}
    function CompatiblePaths(S, C, Z, fixed)
        C_j = Set([c.j for c in C])
        if !all(k∈C_j for k in keys(fixed))
            throw(DomainError("You can only fix chance states."))
        end
        new(S, C, Z, fixed)
    end
end

"""
    CompatiblePaths(S::States, C::Vector{ChanceNode}, Z::DecisionStrategy)

Interface for iterating over paths that are compatible and active given influence diagram and decision strategy.

1) Initialize path `s` of length `n`
2) Fill chance states `s[C]` by generating subpaths `paths(C)`
3) Fill decision states `s[D]` by decision strategy `Z` and path `s`

# Examples
```julia
for s in CompatiblePaths(S, C, Z)
    ...
end
```
"""

function CompatiblePaths(S::States, C::Vector{ChanceNode}, Z::DecisionStrategy)
    CompatiblePaths(S, C, Z, Dict{Node, State}())
end

function compatible_path(S::States, C::Vector{ChanceNode}, Z::DecisionStrategy, s_C::Path)
    s = Array{Int}(undef, length(S))
    for (c, s_C_j) in zip(C, s_C)
        s[c.j] = s_C_j
    end
    for (d, Z_j) in zip(Z.D, Z.Z_j)
        s[d.j] = Z_j((s[d.I_j]...,))
    end
    return (s...,)
end

function Base.iterate(a::CompatiblePaths)
    C_j = [c.j for c in a.C]
    if isempty(a.fixed)
        iter = paths(a.S[C_j])
    else
        ks = sort(collect(keys(a.fixed)))
        fixed = Dict{Int, Int}(i => a.fixed[k] for (i, k) in enumerate(ks))
        iter = paths(a.S[C_j], fixed)
    end
    next = iterate(iter)
    if next !== nothing
        s_C, state = next
        return (compatible_path(a.S, a.C, a.Z, s_C), (iter, state))
    end
end

function Base.iterate(a::CompatiblePaths, gen)
    iter, state = gen
    next = iterate(iter, state)
    if next !== nothing
        s_C, state = next
        return (compatible_path(a.S, a.C, a.Z, s_C), (iter, state))
    end
end

Base.eltype(::Type{CompatiblePaths}) = Path
Base.length(a::CompatiblePaths) = prod(a.S[c.j] for c in a.C)

"""
     UtilityDistribution

UtilityDistribution type.

"""
struct UtilityDistribution
    u::Vector{Float64}
    p::Vector{Float64}
end

"""
    UtilityDistribution(S::States, P::AbstractPathProbability, U::AbstractPathUtility, Z::DecisionStrategy)

Constructs the probability mass function for path utilities on paths that are compatible and active.

# Examples
```julia
UtilityDistribution(S, P, U, Z)
```
"""
function UtilityDistribution(S::States, P::AbstractPathProbability, U::AbstractPathUtility, Z::DecisionStrategy)
    # Extract utilities and probabilities of active paths
    S_Z = CompatiblePaths(S, P.C, Z)
    utilities = Vector{Float64}(undef, length(S_Z))
    probabilities = Vector{Float64}(undef, length(S_Z))
    for (i, s) in enumerate(S_Z)
        utilities[i] = U(s)
        probabilities[i] = P(s)
    end

    # Filter zero probabilities
    nonzero = @. (!)(iszero(probabilities))
    utilities = utilities[nonzero]
    probabilities = probabilities[nonzero]

    # Sort by utilities
    perm = sortperm(utilities)
    u = utilities[perm]
    p = probabilities[perm]

    # Compute the probability mass function
    u2 = unique(u)
    p2 = similar(u2)
    j = 1
    p2[j] = p[1]
    for k in 2:length(u)
        if u[k] == u2[j]
            p2[j] += p[k]
        else
            j += 1
            p2[j] = p[k]
        end
    end

    UtilityDistribution(u2, p2)
end

"""
    StateProbabilities

StateProbabilities type.
"""
struct StateProbabilities
    probs::Dict{Node, Vector{Float64}}
    fixed::Dict{Node, State}
end

"""
    function StateProbabilities(S::States, P::AbstractPathProbability, Z::DecisionStrategy, node::Node, state::State, prev::StateProbabilities)

Associates each node with array of conditional probabilities for each of its states occuring in active paths given fixed states and prior probability.

# Examples
```julia
# Prior probabilities
prev = StateProbabilities(S, P, Z)

# Select node and fix its state
node = 1
state = 2
StateProbabilities(S, P, Z, node, state, prev)
```
"""
function StateProbabilities(S::States, P::AbstractPathProbability, Z::DecisionStrategy, node::Node, state::State, prev::StateProbabilities)
    prior = prev.probs[node][state]
    fixed = prev.fixed
    push!(fixed, node => state)
    probs = Dict(i => zeros(S[i]) for i in 1:length(S))
    for s in CompatiblePaths(S, P.C, Z, fixed), i in 1:length(S)
        probs[i][s[i]] += P(s) / prior
    end
    StateProbabilities(probs, fixed)
end

"""
    function StateProbabilities(S::States, P::AbstractPathProbability, Z::DecisionStrategy)

Associates each node with array of probabilities for each of its states occuring in active paths.

# Examples
```julia
StateProbabilities(S, P, Z)
```
"""
function StateProbabilities(S::States, P::AbstractPathProbability, Z::DecisionStrategy)
    probs = Dict(i => zeros(S[i]) for i in 1:length(S))
    for s in CompatiblePaths(S, P.C, Z), i in 1:length(S)
        probs[i][s[i]] += P(s)
    end
    StateProbabilities(probs, Dict{Node, State}())
end

"""
    function value_at_risk(u::Vector{Float64}, p::Vector{Float64}, α::Float64)

Value-at-risk.
"""
function value_at_risk(u::Vector{Float64}, p::Vector{Float64}, α::Float64)
    @assert 0 ≤ α ≤ 1 "We should have 0 ≤ α ≤ 1."
    i = sortperm(u)
    u, p = u[i], p[i]
    index = findfirst(x -> x≥α, cumsum(p))
    return if index === nothing; u[end] else u[index] end
end

"""
    function conditional_value_at_risk(u::Vector{Float64}, p::Vector{Float64}, α::Float64)

Conditional value-at-risk.
"""
function conditional_value_at_risk(u::Vector{Float64}, p::Vector{Float64}, α::Float64)
    x_α = value_at_risk(u, p, α)
    if iszero(α)
        return x_α
    else
        tail = u .≤ x_α
        return (sum(u[tail] .* p[tail]) - (sum(p[tail]) - α) * x_α) / α
    end
end


text:
biinsert(arr::VecI, val::T) where T = biinsert(arr, val, 1, length(arr)) # @code_warntype ✓
function biinsert(arr::VecI, val::T, lx::Int, rx::Int) where T           # @code_warntype ✓
    lx ≥ rx && return lx
    ub = rx # upper bound
    while lx < rx
        mx = (lx + rx) >> 1                                    # midpoint (binary search)
        @inbounds isless(val, arr[mx]) ? rx = mx : lx = mx + 1 # arr[mx].f == val in this case
    end
    @inbounds lx == ub && !isless(val, arr[lx]) && (lx += 1)   # lx = upper bound && arr[lx] ≤ val
    return lx
end

biinsertsort!(arr::VecI) = biinsertsort!(arr, 1, length(arr)) # @code_warntype ✓
function biinsertsort!(arr::VecI, lx::Int, rx::Int)           # @code_warntype ✓
    for ix in lx+1:rx
        @inbounds val = arr[ix]
        jx = ix
        lc = biinsert(arr, val, lx, ix) # location
        while jx > lc
            swap!(arr, jx, jx - 1)
            jx -= 1
        end
    end
end


text:
include("kepler_eqn.jl")

function calc_true_anom(ecc_anom::Real, e::Real)
	true_anom = 2*atan(sqrt((1+e)/(1-e))*tan(ecc_anom/2))
end

begin
	""" Calculate RV from t, P, K, e, ω and M0	"""
	function calc_rv_keplerian end
	calc_rv_keplerian(t, p::Vector) = calc_rv_keplerian(t, p...)
	function calc_rv_keplerian(t, P,K,e,ω,M0)
		mean_anom = t*2π/P-M0
		ecc_anom = calc_ecc_anom(mean_anom,e)
		true_anom = calc_true_anom(ecc_anom,e)
		rv = K/sqrt((1-e)*(1+e))*(cos(ω+true_anom)+e*cos(ω))
	end
end


text:
using Test, LinearAlgebra
using QBase

@testset "/src/information.jl" begin

@testset "shannon_entropy()" begin
    @test shannon_entropy(Probabilities([1,0,0,0])) == 0
    @test shannon_entropy([1,0,0,0]) == 0

    @test shannon_entropy(Probabilities(fill(1/8,8))) == 3
    @test shannon_entropy(fill(1/8,8)) == 3

    @test shannon_entropy(Probabilities([1/3,1/3,1/3])) == log2(3)
    @test shannon_entropy([1/3,1/3,1/3]) == log2(3)

    @test shannon_entropy(Probabilities([0,0.5,0.5])) == 1
    @test shannon_entropy(Probabilities([0.7,0.3])) ≈ 0.88129089
    @test shannon_entropy(Probabilities([1/2,1/4,1/8,1/8])) == 7/4

    @test_throws DomainError shannon_entropy([-0.1,1.1])
end

@testset "von_neumann_entropy()" begin
    @test von_neumann_entropy(State(diagm(0 => fill(1/4,4)))) == 2
    @test von_neumann_entropy(State([1 0;0 0])) == 0

    @test von_neumann_entropy([0.5 0.5;0.5 0.5]) == 0

    @test von_neumann_entropy.(State.([[1 0;0 0],[0.5 0;0 0.5]])) == [0,1]

    @test_throws DomainError von_neumann_entropy([1 0;0 1])
end

@testset "holevo_bound()" begin
    priors = Probabilities([0.5,0.5])
    ρ_set = State.([[1 0;0 0],[0 0;0 1]])
    @test holevo_bound(priors,ρ_set) == 1
    @test holevo_bound([0.5,0.5],[[1 0;0 0],[0 0;0 1]]) == 1
    @test holevo_bound([0.5,0.5], ρ_set) == 1

    priors = Probabilities([1/3,1/3,1/3])
    ρ_set = mirror_symmetric_qubit_states(0)
    @test holevo_bound(priors,ρ_set) ≈ 0

    priors = Probabilities([1/3,1/3,1/3])
    ρ_set = mirror_symmetric_qubit_states(π/2)
    @test isapprox(holevo_bound(priors,ρ_set), 2/3*log2(3/2)+1/3*log2(3),atol=1e-6)

    priors = Probabilities([1/3,1/3,1/3])
    ρ_set = mirror_symmetric_qubit_states(π/3)
    @test holevo_bound(priors,ρ_set) ≈ 1

    priors = Probabilities([1/3,1/3,1/3])
    ρ_set = mirror_symmetric_qubit_states(5π/12)
    @test isapprox(holevo_bound(priors,ρ_set), 0.957,atol=1e-3)

    qb1 = bloch_qubit_ket(2*π/3,π/2)
    qb2 = bloch_qubit_ket(2*π/3,3π/2)
    qb3 = bloch_qubit_ket(0,0)
    ρ_set = pure_state.([qb1,qb2,qb3])
    @test holevo_bound(priors,ρ_set) ≈ 1

    qb1 = bloch_qubit_state(2*π/3,0)
    qb2 = bloch_qubit_state(5*π/12,π)
    qb3 = bloch_qubit_state(0,0)
    @test holevo_bound(priors,[qb1,qb2,qb3]) ≈ 0.952526321

    states = [State([1 0;0 0]), State([1 0;0 0]), State([0 0;0 1])]
    @test holevo_bound(priors, states) ≈ 0.91829583

    priors = Probabilities(ones(4)/4)
    @test holevo_bound(priors, sic_qubit_states()) == 1
end

@testset "mutual_information()" begin
    priors = Probabilities([1/3,1/3,1/3])
    ρ_set = trine_qubit_states()
    povm = sqrt_povm(priors, ρ_set)
    @test mutual_information(priors, ρ_set, povm) ≈ 1/3

    @test mutual_information([0.5,0.5],[[1 0;0 0],[0 0;0 1]],[[1 0;0 0],[0 0;0 1]]) == 1

    rot_y = qubit_rotation(π,axis="y")

    rot_povm = POVM(map(el -> rot_y*el*rot_y', povm))

    @test mutual_information(priors, ρ_set, rot_povm) ≈ 0.5849625

    orth_povm = POVM([[0.5 0;0 0],[0.5 0;0 0],[0 0;0 1]])
    @test mutual_information(priors, ρ_set, orth_povm) ≈ 0.459147917
end

@testset "joint_entropy()" begin
    priors = Probabilities([0.5,0.5])
    conditionals = Conditionals([0.5 0.5;0.5 0.5])

    @test joint_entropy(priors, conditionals) == 2
    @test joint_entropy(priors, [0.5 0.5;0.5 0.5]) == 2

    priors = Probabilities([1/3,1/3,1/3])
    conditionals = Conditionals([2/3 1/6 1/6;1/6 2/3 1/6;1/6 1/6 2/3])
    @test joint_entropy(priors, conditionals) ≈ 2.8365917
end

@testset "mutual_information()" begin
    priors = Probabilities([1/3,1/3,1/3])
    conditionals = Conditionals([2/3 1/6 1/6;1/6 2/3 1/6;1/6 1/6 2/3])
    @test mutual_information(priors,conditionals) ≈ 1/3
    @test mutual_information([1/3,1/3,1/3],conditionals) ≈ 1/3

    priors = Probabilities([1/3,1/3,1/3])
    conditionals = Conditionals([0 1/2 1/2;1/2 0 1/2;1/2 1/2 0])
    @test mutual_information(priors,conditionals) ≈ 0.5849625

    conditionals = Conditionals([1 0 1;0 1 0;0 0 0])
    @test mutual_information(priors,conditionals) ≈ 2/3*log2(3/2)+1/3*log2(3)

    conditionals = Conditionals([1 0 0;0 1 0;0 0 1])
    @test mutual_information(priors,conditionals) ≈ log2(3)

    conditionals = Conditionals([1/2 1/8 1/8;1/2 1/8 1/8;0 3/4 3/4])
    @test mutual_information(priors,conditionals) ≈ 0.459147917
end

@testset "conditional_entropy()" begin
    priors = Probabilities([1/3, 1/3, 1/3])
    conditionals = Conditionals([2/3 1/6 1/6;1/6 2/3 1/6;1/6 1/6 2/3])
    @test conditional_entropy(priors, conditionals) ≈ 1.251629167

    cond_ent = conditional_entropy(priors, conditionals)
    mut_inf = mutual_information(priors,conditionals)
    @test cond_ent + mut_inf ≈ shannon_entropy(priors)

    priors = Probabilities([1/3, 1/3, 1/3])
    conditionals = Conditionals([3/4 1/8 1/8;1/8 3/4 1/8;1/8 1/8 3/4])
    conditional_entropy(priors, conditionals)

    cond_ent = conditional_entropy(priors, conditionals)
    mut_inf = mutual_information(priors,conditionals)
    @test cond_ent + mut_inf ≈ shannon_entropy(priors)
end

@testset "success_probability()" begin
    @testset "trivial case" begin
        Π = POVM([[1 0;0 0],[0 0;0 1]])
        ρ_states = State.([[1 0;0 0],[0 0;0 1]])
        priors = Probabilities([1/2,1/2])

        @test success_probability(priors, ρ_states, Π) == 1
        @test success_probability([1/2,1/2], [[1 0;0 0],[0 0;0 1]], Π) == 1
    end
end

@testset "error_probability()" begin
    @testset "trivial case" begin
        Π = POVM([[1 0;0 0],[0 0;0 1]])
        ρ_states = State.([[1 0;0 0],[0 0;0 1]])
        priors = Probabilities([1/2,1/2])

        @test error_probability(priors, ρ_states, Π) == 0
    end
end

end


text:
import brainflow

# specify logging library to use
brainflow.enable_dev_brainflow_logger(Integer(brainflow.BOARD_CONTROLLER))

params = brainflow.BrainFlowInputParams()
board_shim = brainflow.BoardShim(Integer(brainflow.SYNTHETIC_BOARD), params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(32, board_shim)
brainflow.release_session(board_shim)

brainflow.write_file(data, "test.csv", "w")
restored_data = brainflow.read_file("test.csv")

println("Original Data")
println(data)
println("Restored Data")
println(restored_data)

max_stars_repo_name: Cognixion-inc/brainflow

max_issues_repo_name: Cognixion-inc/brainflow

max_forks_repo_name: Cognixion-inc/brainflow