text:
export ## Types
       ## Methods
       unitary_propagator

CompositeQSystem() = CompositeQSystem(QSystem[], Interaction[], ParametricInteraction[], Vector{Vector{Int}}[], Vector{Vector{Int}}[], Dissipation[])

function getindex(c::CompositeQSystem, key::String)
    for s in c.subSystems
        if label(s) == key
            return s
        end
    end
    throw(KeyError(key))
end

+(s1::QSystem, s2::QSystem) = (c = CompositeQSystem(); c += s1; c += s2; c)
+(s::QSystem, i::Interaction) = (c = CompositeQSystem(); c += s; c += i; c)

function +(c::CompositeQSystem, q::QSystem)
    append!(c.subSystems, [q])
    update_expansion_indices!(c)
    return c
end

function +(c::CompositeQSystem, i::Interaction)
    append!(c.interactions, [i])
    update_expansion_indices!(c)
    return c
end

function +(c::CompositeQSystem, pi::ParametricInteraction)
    append!(c.parametericInteractions, [pi])
    return c
end

function +(c::CompositeQSystem, d::Dissipation)
    append!(c.dissipators, [d])
end

function update_expansion_indices!(c::CompositeQSystem)
    c.subSystemExpansions = [Vector{Int}[] for _ = 1:length(c.subSystems)]
    for (ct, sys) in enumerate(c.subSystems)
        c.subSystemExpansions[ct] = expand_indices([ct], dims(c))
    end

    c.interactionExpansions = [Vector{Int}[] for _ = 1:length(c.interactions)]
    for (ct, i) in enumerate(c.interactions)
        c.interactionExpansions[ct] = expand_indices(find_subsystem_pos(c, i), dims(c))
    end
end

labels(c::CompositeQSystem) = [label(s) for s in c.subSystems]
dims(c::CompositeQSystem) = [dim(s) for s in c.subSystems]
dim(c::CompositeQSystem) = prod([dim(s) for s in c.subSystems])

function find_subsystem_pos(c::CompositeQSystem, s::QSystem)
    @assert s in c.subSystems "Oops! Subsystem not found."
    findin(c.subSystems, [s])
end

function find_subsystem_pos(c::CompositeQSystem, i::Interaction)
    #Field-system interactions are one-body terms
    if isa(i, SemiClassicalDipole)
        return find_subsystem_pos(c, i.system2)
    elseif isa(i, FluxTransmon)
        return find_subsystem_pos(c, i.transmon)
    else
        return [find_subsystem_pos(c, i.system1), find_subsystem_pos(c, i.system2)]
    end
end

function hamiltonian(c::CompositeQSystem, t::Float64)

    #Initialize Hamiltonian
    Htot = zeros(Complex128, dim(c), dim(c))

    #Add in all the terms
    hamiltonian_add!(Htot, c, t)

    return Htot
end
hamiltonian(c::CompositeQSystem) = hamiltonian(c, 0.0)

function hamiltonian_add!(Ham::Matrix{Complex128}, c::CompositeQSystem, t::Float64)
    #Fast system hamiltonian calculator with total Hamiltonian preallocated

    #Zero the Hamiltonian memory
    Ham[:] = 0.0

    #Update the subsystems with the parameteric interactions
    for pi in c.parametericInteractions
        update_params(c, pi, t)
    end

    #Add together subsystem Hamiltonians
    for (subsys, expander) in zip(c.subSystems, c.subSystemExpansions)
        expand_add!(Ham, hamiltonian(subsys, t), expander)
    end

    #Add interactions
    for (i, expander) in zip(c.interactions, c.interactionExpansions)
        expand_add!(Ham, hamiltonian(i,t), expander)
    end
end

function expand(m::Matrix, actingOn::Vector, dims::Vector)
    #Expand an operator onto a larger Hilbert space
    # m: matrix form of  operator
    # actingOn: array of which subsystem index the operator should be acting on
    # dims: array of dimensions of all the subsystems

    @assert size(m, 1) == prod(dims[actingOn]) "Oops! Dimensions of matrix do not match dims argument."

    #Create the large matrix by tensoring on identity
    l = length(dims)
    eyeIndices = filter(x->!(x in actingOn), 1:l)
    M = isempty(eyeIndices) ? m : kron(m, eye(eltype(m), prod(dims[eyeIndices])))

    #Reshape into multi-dimensional array given by subsystem dimensions
    #Since we have a matrix we repeat for rows then columns
    M = reshape(M, tuple([dims, dims]...))

    #Permute magic
    forwardPerm = [actingOn, eyeIndices]
    reversePerm = invperm(forwardPerm)
    #Handle the way tensor product indices work (last subsystem is fastest)
    reversePerm = reverse((l+1) .- reversePerm)
    M = permutedims(M, [reversePerm, reversePerm .+ l])

    #Reshape back
    return reshape(M, prod(dims), prod(dims))
end

function expand(m::Matrix, indices::Vector, sizeM::Int )
    M = zeros(eltype(m), (sizeM, sizeM))
    for (ct, inds) in enumerate(indices)
        M[inds] = m[ct]
    end
    return M
end

function expand_add!(M::Matrix, m::Matrix, indices::Vector )
    #Add to certain indices of M with terms from m according to expansion indices.
    for ct=1:length(indices)
        # M[indices[ct]] += m[ct]
        for idx = indices[ct]
            M[idx] += m[ct]
        end
    end
end

function expand_indices(actingOn::Vector, dims::Vector)
    #Calculate the indices for expansion
    actingOnDim = prod(dims[actingOn])
    sm = (actingOnDim, actingOnDim)
    lenm = actingOnDim^2;
    M = expand(reshape(1:lenm, sm), actingOn, dims)
    return [find(M .== x) for x in 1:lenm]
end


text:
module TestIO
    using Base.Test
    using DataFrames

    #test_group("Confirm that we can read various file types.")
    testdir = dirname(@__FILE__)

    filenames = ["$testdir/data/blanklines/blanklines.csv",
                 "$testdir/data/compressed/movies.csv.gz",
                 "$testdir/data/newlines/os9.csv",
                 "$testdir/data/newlines/osx.csv",
                 "$testdir/data/newlines/windows.csv",
                 "$testdir/data/newlines/embedded_os9.csv",
                 "$testdir/data/newlines/embedded_osx.csv",
                 "$testdir/data/newlines/embedded_windows.csv",
                 "$testdir/data/padding/space_after_delimiter.csv",
                 "$testdir/data/padding/space_around_delimiter.csv",
                 "$testdir/data/padding/space_before_delimiter.csv",
                 "$testdir/data/quoting/empty.csv",
                 "$testdir/data/quoting/escaping.csv",
                 "$testdir/data/quoting/quotedcommas.csv",
                 "$testdir/data/scaling/10000rows.csv",
                 "$testdir/data/scaling/movies.csv",
                 "$testdir/data/separators/sample_data.csv",
                 "$testdir/data/separators/sample_data.tsv",
                 "$testdir/data/separators/sample_data.wsv",
                 "$testdir/data/typeinference/bool.csv",
                 "$testdir/data/typeinference/standardtypes.csv",
                 "$testdir/data/utf8/corrupt_utf8.csv",
                 "$testdir/data/utf8/short_corrupt_utf8.csv",
                 "$testdir/data/utf8/utf8.csv"]

    for filename in filenames
        try
            df = readtable(filename)
        catch
            error(@sprintf "Failed to read %s\n" filename)
        end
    end

    # Spot check movies.csv file
    filename = "$testdir/data/scaling/movies.csv"
    df = readtable(filename)
    @test df[1, 1] === 1
    # TODO: Figure out why strict equality won't work here
    #       Doesn't seem to be UTF8String vs. ASCIIString
    @test df[1, 2] == "\$"
    @test df[1, 3] === 1971
    @test df[1, 4] === 121
    @test df[1, 5] === NA
    @test df[1, 6] === 6.4
    @test df[1, 7] === 348
    @test df[1, 8] === 4.5
    @test df[1, 9] === 4.5
    @test df[1, 10] === 4.5
    @test df[1, 11] === 4.5
    @test df[1, 12] === 14.5
    @test df[1, 13] === 24.5
    @test df[1, 14] === 24.5
    @test df[1, 15] === 14.5
    @test df[1, 16] === 4.5
    @test df[1, 17] === 4.5
    # TODO: Figure out why strict equality won't work here
    #       Doesn't seem to be UTF8String vs. ASCIIString
    @test df[1, 18] == ""
    @test df[1, 19] === 0
    @test df[1, 20] === 0
    @test df[1, 21] === 1
    @test df[1, 22] === 1
    @test df[1, 23] === 0
    @test df[1, 24] === 0
    @test df[1, 24] === 0

    @test df[end, 1] === 58788
    # TODO: Figure out why strict equality won't work here
    #       Doesn't seem to be UTF8String vs. ASCIIString
    @test df[end, 2] == "xXx: State of the Union"
    @test df[end, 3] === 2005
    @test df[end, 4] === 101
    @test df[end, 5] === 87000000
    @test df[end, 6] === 3.9
    @test df[end, 7] === 1584
    @test df[end, 8] === 24.5
    @test df[end, 9] === 4.5
    @test df[end, 10] === 4.5
    @test df[end, 11] === 4.5
    @test df[end, 12] === 4.5
    @test df[end, 13] === 14.5
    @test df[end, 14] === 4.5
    @test df[end, 15] === 4.5
    @test df[end, 16] === 4.5
    @test df[end, 17] === 14.5
    # TODO: Figure out why strict equality won't work here
    #       Doesn't seem to be UTF8String vs. ASCIIString
    @test df[end, 18] == "PG-13"
    @test df[end, 19] === 1
    @test df[end, 20] === 0
    @test df[end, 21] === 0
    @test df[end, 22] === 0
    @test df[end, 23] === 0
    @test df[end, 24] === 0
    @test df[end, 25] === 0

    df1 = readtable("$testdir/data/comments/before_after_data.csv", allowcomments = true)
    df2 = readtable("$testdir/data/comments/middata.csv", allowcomments = true)
    df3 = readtable("$testdir/data/skiplines/skipfront.csv", skipstart = 3)

    @test df1 == df2
    @test df2 == df3

    df1 = readtable("$testdir/data/separators/sample_data.csv")
    df2 = readtable("$testdir/data/separators/sample_data.tsv")
    df3 = readtable("$testdir/data/separators/sample_data.wsv")
    df4 = readtable("$testdir/data/separators/sample_data_white.txt", separator = ' ')

    @test df1 == df2
    @test df2 == df3
    @test df3 == df4

    readtable("$testdir/data/quoting/quotedwhitespace.txt", separator = ' ')

    # TODO: Implement skipping lines at specified row positions
    # readtable("$testdir/data/skiplines/skipbottom.csv", skiprows = [1, 2, 3])

    # TODO: Implement skipping lines at bottom
    # readtable("$testdir/data/skiplines/skipbottom.csv", skipstartlines = 4)

    #test_group("Confirm that we can read a large file.")

    df = DataFrame()

    nrows, ncols = 100_000, 10

    for j in 1:(ncols - 1)
        df[j] = randn(nrows)
    end
    df[ncols] = "A"

    filename = tempname()

    writetable(filename, df, separator = ',')

    df1 = readtable(filename, separator = ',')

    @test isequal(df, df1)

    rm(filename)

    #
    # Lots of rows
    #

    df = DataFrame()

    nrows, ncols = 1_000_000, 1

    for j in 1:ncols
        df[j] = randn(nrows)
    end

    filename = tempname()

    writetable(filename, df, separator = ',')

    df1 = readtable(filename, separator = ',')

    @test isequal(df, df1)

    rm(filename)

    #test_group("Properties of data frames returned by readtable method.")

    # Readtable with makefactors active should only make factors from columns
    # of strings.
    filename = "$testdir/data/factors/mixedvartypes.csv"
    df = readtable(filename, makefactors = true)

    @test typeof(df[:factorvar]) == PooledDataArray{UTF8String,Uint32,1}
    @test typeof(df[:floatvar]) == DataArray{Float64,1}

    # Readtable shouldn't silently drop data when reading highly compressed gz.
    df = readtable("$testdir/data/compressed/1000x2.csv.gz")
    @test size(df) == (1000, 2)

    # Readtable type inference
    filename = "$testdir/data/typeinference/bool.csv"
    df = readtable(filename)
    @test typeof(df[:Name]) == DataArray{UTF8String,1}
    @test typeof(df[:IsMale]) == DataArray{Bool,1}
    @test df[:IsMale][1] == true
    @test df[:IsMale][4] == false

    filename = "$testdir/data/typeinference/standardtypes.csv"
    df = readtable(filename)
    @test typeof(df[:IntColumn]) == DataArray{Int,1}
    @test typeof(df[:IntlikeColumn]) == DataArray{Float64,1}
    @test typeof(df[:FloatColumn]) == DataArray{Float64,1}
    @test typeof(df[:BoolColumn]) == DataArray{Bool,1}
    @test typeof(df[:StringColumn]) == DataArray{UTF8String,1}

    filename = "$testdir/data/typeinference/mixedtypes.csv"
    df = readtable(filename)
    @test typeof(df[:c1]) == DataArray{UTF8String,1}
    @test df[:c1][1] == "1"
    @test df[:c1][2] == "2.0"
    @test df[:c1][3] == "true"
    @test typeof(df[:c2]) == DataArray{Float64,1}
    @test df[:c2][1] == 1.0
    @test df[:c2][2] == 3.0
    @test df[:c2][3] == 4.5
    @test typeof(df[:c3]) == DataArray{UTF8String,1}
    @test df[:c3][1] == "0"
    @test df[:c3][2] == "1"
    @test df[:c3][3] == "f"
    @test typeof(df[:c4]) == DataArray{Bool,1}
    @test df[:c4][1] == true
    @test df[:c4][2] == false
    @test df[:c4][3] == true
    @test typeof(df[:c5]) == DataArray{UTF8String,1}
    @test df[:c5][1] == "False"
    @test df[:c5][2] == "true"
    @test df[:c5][3] == "true"

    # Readtable defining column types
    filename = "$testdir/data/definedtypes/mixedvartypes.csv"

    df = readtable(filename)
    @test typeof(df[:n]) == DataArray{Int,1}
    @test df[:n][1] == 1
    @test typeof(df[:s]) == DataArray{UTF8String,1}
    @test df[:s][1] == "text"
    @test typeof(df[:f]) == DataArray{Float64,1}
    @test df[:f][1] == 2.3
    @test typeof(df[:b]) == DataArray{Bool,1}
    @test df[:b][1] == true

    df = readtable(filename,eltypes=[Int64, UTF8String, Float64, Bool])
    @test typeof(df[:n]) == DataArray{Int64,1}
    @test df[:n][1] == 1
    @test typeof(df[:s]) == DataArray{UTF8String,1}
    @test df[:s][1] == "text"
    @test df[:s][4] == "text ole"
    @test typeof(df[:f]) == DataArray{Float64,1}
    @test df[:f][1] == 2.3
    @test typeof(df[:b]) == DataArray{Bool,1}
    @test df[:b][1] == true
    @test df[:b][2] == false

    df = readtable(filename,eltypes=[Int64, UTF8String, Float64, UTF8String])
    @test typeof(df[:n]) == DataArray{Int64,1}
    @test df[:n][1] == 1.0
    @test isna(df[:s][3])
    @test typeof(df[:f]) == DataArray{Float64,1}
    # Float are not converted to int
    @test df[:f][1] == 2.3
    @test df[:f][2] == 0.2
    @test df[:f][3] == 5.7
    @test typeof(df[:b]) == DataArray{UTF8String,1}
    @test df[:b][1] == "T"
    @test df[:b][2] == "FALSE"

    # Readtable ignorepadding argument
    io = IOBuffer("A , \tB  , C\n1 , \t2, 3\n")
    @test readtable(io, ignorepadding=true) == DataFrame(A=1, B=2, C=3)

    # Readtable name normalization
    abnormal = "\u212b"
    ns = [:Å, :_B_C_, :_end]
    @test !in(symbol(abnormal), ns)

    io = IOBuffer(abnormal*",%_B*\tC*,end\n1,2,3\n")
    @test names(readtable(io)) == ns
end

ext: jl


text:
"Outer paving with a contractor C"

function outer_pave{N,T}(C, working::Vector{IntervalBox{N,T}}, ϵ)

    outer_list = SubPaving{N,T}()

    while !isempty(working)

        X = pop!(working)

        contracted = C(X)

        if isempty(contracted)
            continue
        end

        if diam(contracted) < ϵ
            push!(outer_list, contracted)

        else
            push!(working, bisect(contracted)...)
        end

    end

    return outer_list

end

outer_pave(C, X, ϵ) = outer_pave(C, [X], ϵ)


text:
function np_to_jl_type(o::PyObject)
    t = o.name
    if t == "float32"
        return Float32
    elseif t == "float64"
        return Float64
    elseif t == "complex64"
        return ComplexF32
    elseif t == "complex128"
        return Complexf64
    elseif t == "int32"
        return Int32
    elseif t == "int64"
        return Int64
    elseif t == "uint32"
        return UInt32
    elseif t == "uint64"
        return UInt64
    elseif t == "bool"
        return Bool
    else
        error("Unknown scalar type $t")
    end
end


text:
# import Flux
# import GPUArrays
# import NNlib
import ProgressMeter
import ValueHistories
# import Zygote


text:
include("config.jl")

println("Adding workers")
include("addprocs_external.jl")
include("test_workers.jl")

if isfile("compute.jl")
    include("compute.jl")
end


text:
function cudnnCreate()
    handle_ref = Ref{cudnnHandle_t}()
    res = @retry_reclaim err->isequal(err, CUDNN_STATUS_ALLOC_FAILED) ||
                              isequal(err, CUDNN_STATUS_NOT_INITIALIZED) ||
                              isequal(err, CUDNN_STATUS_INTERNAL_ERROR) begin
        unsafe_cudnnCreate(handle_ref)
    end
    if res != CUDNN_STATUS_SUCCESS
         throw_api_error(res)
    end
    return handle_ref[]
end

@memoize function cudnnGetProperty(property::CUDA.libraryPropertyType)
  value_ref = Ref{Cint}()
  cudnnGetProperty(property, value_ref)
  value_ref[]
end

@memoize function version()
  ver = cudnnGetVersion()
  major, ver = divrem(ver, 1000)
  minor, patch = divrem(ver, 10)

  VersionNumber(major, minor, patch)
end

@memoize function cuda_version()
  ver = cudnnGetCudartVersion()
  major, ver = divrem(ver, 1000)
  minor, patch = divrem(ver, 10)

  VersionNumber(major, minor, patch)
end


text:
using Images
using Plots
include("funcs.jl")

let

    cd("../../erickson/output_files/d4_p4_spinup")
    @printf "found directory...\n"
    slip_file = collect(eachline(open("slip.dat", "r")))
    stress_data = collect(eachline(open("stress.dat","r")))

    nn = Int64(get_line(slip_file,1)[1])
    y = get_line(slip_file, 2)[3:end]
    slip_data = @view slip_file[3:end , :]

    img1 = load(string("V_png/velocities.", lpad(0,4,"0"), ".png"))
    plt1 = plot(img1, yaxis=false, xaxis=false, ticks=false)
    img2 = load(string("U_png/displacements.", lpad(0,4,"0"), ".png"))
    plt2 = plot(img2,  yaxis=false, xaxis=false, ticks=false)
    plt3 = plot(legend=false, yflip = true, ylabel="Depth(Km)", xlabel="Slip(m)")
    plt4 = plot(legend=false, yflip = true, ylabel="Depth(Km)", xlabel="Stress(MPa)")

    @printf "intialized plotting...\n"
    l = @layout [
        a{.7w} b{.3w}
        c{.7w} d{.3w}
    ]

    offset = zeros(nn)
    inter_flag = 1
    break_count = 0
    inter_count = 0
    cycle_count = 1
    final_index = size(slip_data)[1]

    for i = 1:final_index
        @printf "\r%f" i/final_index
        stress_index = i - break_count

        if inter_flag == 1
            if slip_data[i] == "BREAK"
                inter_flag = 0
                break_count += 1
            else

                slip = get_line(slip_data, i)[3:end]
                stress = get_line(stress_data, stress_index)[3:end]
                plt3 = plot!(plt3, slip .- offset, y, linecolor=:blue, linewidth=.1, title=string("Cycle # ", cycle_count))
                plt4 = plot(stress, y,
                            legend=false,
                            linecolor=:black,
                            flip = true,
                            ylabel="Depth(Km)",
                            xlabel="Stress(MPa)")
                plot(plt1, plt3, plt2, plt4, layout=l)
                #gui()
                inter_count += 1
            end
        else

            if slip_data[i] == "BREAK"
                inter_flag = 1
                break_count += 1
                cycle_count += 1
                if break_count == 12
                    offset .= get_line(slip_data, i+1)[3:end]
                    plt3 = plot(legend=false, yflip = true, ylabel="Depth(Km)", xlabel="Slip(m)")
                end
            else

                slip = get_line(slip_data, i)
                V_max = slip[2]

                if V_max > log(.1)

                    particle_index = i - break_count - inter_count
                    img1 = load(string("V_png/velocities.", lpad(particle_index,4,"0"), ".png"))
                    plt1 = plot(img1, yaxis=false, xaxis=false, ticks=false)
                    img2 = load(string("U_png/displacements.", lpad(particle_index,4,"0"), ".png"))
                    plt2 = plot(img2, yaxis=false, xaxis=false, ticks=false)
                    slip = slip[3:end]
                    stress = get_line(stress_data, stress_index)[3:end]
                    plt3 = plot!(plt3, slip .- offset, y, linecolor=:red, linewidth=.1)
                    plt4 = plot(stress, y,
                                legend=false,
                                linecolor=:black,
                                flip = true,
                                ylabel="Depth(Km)",
                                xlabel="Stress(MPa)")

                    plot(plt1, plt3, plt2, plt4, layout=l)
                    #gui()
                end

            end
        end
    end

    gif(anim, "ruptures.gif", fps=16)
    cd("/home/tharvey2/Basin_Simulations/utils")
end


text:
# This is not run by runtests.jl. Only intended to be checked manually.
using Plots
const test_funcs = [
                        (x -> 2, 0, 1),
                        (x-> x, -1, 1),
                        (x-> 5x, -1, 1),
                        (x-> 1/x, 0, 1),
                        (x-> 1/x, -1, 1),
                        (x->tan(x), 0, 4),
                        (x-> -1/x, 0, 1),
                        (x-> 1/abs(x),-1,1),
                        (x-> 1e6x, -1, 1),
                        (x-> 1e50x, -1, 1),
                        (x-> log(1+sin(cos(x))), -6, 6),
                        (x-> sin(x^3)+cos(x^3), 0, 6.28),
                        (x-> sin(x), -5, 200),
                        (x-> sin(1/x), -2, 2),
                        (x-> sin(1/x), 0, 2),
                        (x-> sin(x^4), -4, 4),
                        (x-> sin(300x), -4, 4),
                        (x-> 1 + x^2 + 0.0125log(abs(1-3(x-1))), -2, 2),
                        (x-> sin(exp(x)), -6, 6),
                        (x-> 1/sin(x), -10, 10),
                        (x-> sin(x)/x, -6, 6),
                        (x-> tan(x^3-x+1) + 1/(x+3exp(x)), -2, 2),
                    ]
##
plots = [plot(func...) for func in test_funcs]
pitchfork = plot(x->sqrt(x),0,10)
plot!(pitchfork, x->-sqrt(x),0,10)
plot!(pitchfork, x->0,-10,0)
plot!(pitchfork, x->0,0,10, linestyle=:dash)
push!(plots, pitchfork)
display.(plots)


text:
# Types for GeMM

"""
A simple struct to store organism trait names and values.
"""
mutable struct Trait
    nameindex::UInt8
    value::Float64
end

"""
An abstraction layer to hide the implementation details of genes.
"""
abstract type AbstractGene end

"""
The standard gene struct; including the base sequence (represented in binary and
stored and as an integer), and an array of traits the gene codes for.
"""
mutable struct Gene <: AbstractGene
    sequence::Int
    codes::Array{Trait, 1}
end

"""
A variation of the standard Gene struct that allows for longer base sequence
lengths through its use of BigInt.

Used for the compatibility gene when settings["usebiggenes"] == true.
"""
mutable struct BigGene <: AbstractGene
    sequence::BigInt
    codes::Array{Trait, 1}
end

"""
This gene type cuts out the conversion of the sequence to an integer. (This
conversion was introduced to save memory, but with low population sizes
(as in Zosterops mode), saving computation cycles may be more important.)

Used when settings["compressgenes"] == false.
"""
mutable struct StringGene <: AbstractGene
    sequence::String
    codes::Array{Trait, 1}
end

"""
A struct that saves a collection of genes and a tag to denote whether they are
maternal or not.
"""
mutable struct Chromosome
    genes::Array{AbstractGene, 1}
    maternal::Bool # parental origin of chromosome
    #TODO keep track of the lineage origin to measure heterozygosity?
end

@enum Sex hermaphrodite male female

"""
One of the core structs of the model, representing an individual organism.
"""
mutable struct Individual
    lineage::String
    genome::Array{Chromosome, 1}
    traits::Dict{String, Float64}
    marked::Bool  # indicator whether individual is new to a patch (after dispersal or birth)
    precadaptation::Float64 # adaption to precipitation
    tempadaptation::Float64 # adaption to temperature
    size::Float64 # body mass
    sex::Sex
    partner::Int  # ID of the partner individual (if applicable, default 0)
    id::Int
end

"""
One of the core structs of the model, representing a one-hectare patch of ground.
"""
mutable struct Patch
    id::Int
    location::Tuple{Int, Int}
    capacity::Float64   # biomass carrying capacity in g ("cellsize")
    temp::Float64   # temperature (physiologically important)
    prec::Float64   # precipitation (no physiological effect, just a generic niche)
    nicheb::Float64 # additional generic niche - currently not used
    community::Array{Individual, 1}
    seedbank::Array{Individual, 1}
    isisland::Bool  # island? (if false -> mainland)
    invasible::Bool # can exotics land here?
    isolated::Bool  # add a distance penalty when dispersing?
    initpop::Bool   # initialise with a population?
end

# constructors:
# XXX These default values should be defined elsewhere (defaults.jl, to be precise)
Patch(id, location, capacity) =
    Patch(id, location, capacity, 298, 5, 0, Individual[], Individual[],
          false, false, false, false)

max_stars_repo_name: lleiding/gemm

max_issues_repo_name: lleiding/gemm

max_forks_repo_name: lleiding/gemm


text:
function BOMCP.x2s(m::LQGPOMDP, x::Vector{Float64})
    s = x
    return s
end

function BOMCP.s2x(m::LQGPOMDP, s::Vector{Float64})
    x = s
    return x
end

function BOMCP.gen_A(m::LQGPOMDP, s::Vector{Float64}, a::Vector{Float64})
    A = zeros(Float64, 3, 3)
    A[1,1] = m.A_matrix[1,1]
    A[1,2] = m.A_matrix[1,2]
    A[2,1] = m.A_matrix[2,1]
    A[2,2] = m.A_matrix[2,2]
    A[3,3] = (s[end] + 1)/(s[end])
    return A
end

function BOMCP.gen_C(m::LQGPOMDP, s::Vector{Float64})
    C = zeros(Float64, 3, 3)
    C[1,1] = m.C_matrix[1,1]
    C[1,2] = m.C_matrix[1,2]
    C[2,1] = m.C_matrix[2,1]
    C[2,2] = m.C_matrix[2,2]
    C[3,3] = 1.0
    return C
end

function BOMCP.vectorize!(v, dims, x::MvNormal)
    v = copy(mean(x))
    return v
end

POMDPs.updater(p::LQGPolicy) = EKFUpdater(p.m, [0.1, 0.1, 1e-300].^2, [0.1, 0.1, 1e-300].^2)

function POMDPs.update(up::EKFUpdater, b::MvNormal, a::Vector{Float64}, o::Vector{Float64})
    μ = mean(b)
    n = length(μ)
    Σ = cov(b)
    s = x2s(up.m, μ)
    # Predict
    sp, z = POMDPs.gen(DDNOut(:sp, :o), up.m, s, a, Random.GLOBAL_RNG)
    xp = s2x(up.m, sp)

    A = gen_A(up.m, μ, a)
    C = gen_C(up.m, xp)

    Σ_hat = A*Σ*transpose(A) + up.Q

    # Update
    y = o - z

    S = C*Σ_hat*transpose(C) + up.R
    K = Σ_hat*transpose(C)/S

    μp = xp + K*y

    Σp = (Matrix{Float64}(I, n, n) - K*C)*Σ_hat

    Σp = round.(Σp, digits=5)
    Σp[end, end] = 1e-300
    bp = MvNormal(μp, Σp)
    return bp
end

function BasicPOMCP.extract_belief(up::EKFUpdater, bn)
    pomcpow_belief = belief(bn).sr_belief

    n_obs = length(pomcpow_belief.dist.cdf)
    X = zeros(Float64, 3, n_obs)
    for (i, item) in enumerate(pomcpow_belief.dist.items)
        X[:,i] = item[1]
    end
    weights = [pomcpow_belief.dist.cdf[1]]
    for i = 2:n_obs
        dw = pomcpow_belief.dist.cdf[i] - pomcpow_belief.dist.cdf[i-1]
        append!(weights, dw)
    end
    W = ProbabilityWeights(weights)

    if n_obs > 1
        μ = mean(X, W, dims=2)[:,1]
        σ = cov(X, W)
        display(σ)
    else
        μ = X[:,1]
        σ = [0.1, 0.1, 1e-300]
        σ = diagm(σ)
        # display(σ)
    end
    # display(μ)
    return MvNormal(μ, σ)
end


text:
export energy, meantime, rmsduration, meanfrequency, rmsbandwidth, ifrequency

"""
$(SIGNATURES)
Computes total signal energy.
"""
energy(s::AbstractVector; fs=framerate(s)) = sum(abs2, s)/inHz(fs)
energy(s::AbstractMatrix; fs=framerate(s)) = vec(sum(abs2, s; dims=1))./inHz(fs)

"""
$(SIGNATURES)
Computes mean time of the signal.
"""
meantime(s::SampleBuf) = wmean(domain(s), abs2.(samples(s)))
meantime(s; fs) = wmean((0:size(s,1)-1)/fs, abs2.(s))

"""
$(SIGNATURES)
Computes RMS duration of the signal.
"""
rmsduration(s::SampleBuf) = sqrt.(wmean(domain(s).^2, abs2.(samples(s))) .- meantime(s).^2)
rmsduration(s; fs) = sqrt.(wmean(((0:size(s,1)-1)/fs).^2, abs2.(s)) .- meantime(s; fs=fs).^2)

"""
$(SIGNATURES)
Computes instantaneous frequency of the signal.
"""
function ifrequency(s; fs=framerate(s))
  s1 = analytic(s)
  f1 = inHz(fs)/(2π) * diff(unwrap(angle.(s1); dims=1); dims=1)
  f2 = Array{Float32}(undef, size(s))
  f2[1,:] = f1[1,:]
  f2[2:end-1,:] .= (f1[1:end-1,:] .+ f1[2:end,:])./2
  f2[end,:] = f1[end,:]
  signal(f2, fs)
end

"""
$(SIGNATURES)
Computes mean frequency of a signal.
"""
function meanfrequency(s::AbstractMatrix; fs=framerate(s), nfft=1024, window=nothing)
  fs = inHz(fs)
  rv = mapslices(samples(s); dims=1) do s1
    p = welch_pgram(s1, min(nfft, length(s1)); fs=fs, window=window)
    f = freq(p)
    wmean(f, power(p))
  end
  dropdims(rv, dims=1)
end

function meanfrequency(s::AbstractVector; fs=framerate(s), nfft=1024, window=nothing)
  p = welch_pgram(s, min(nfft, length(s)); fs=inHz(fs), window=window)
  f = freq(p)
  wmean(f, power(p))
end

"""
$(SIGNATURES)
Computes RMS bandwidth of a signal.
"""
function rmsbandwidth(s::AbstractMatrix; fs=framerate(s), nfft=1024, window=nothing) where T
  fs = inHz(fs)
  rv = mapslices(samples(s); dims=1) do s1
    p = welch_pgram(s1, min(nfft, length(s1)); fs=fs, window=window)
    f = freq(p)
    f0 = wmean(f, power(p))
    sqrt.(wmean((f.-f0).^2, power(p)))
  end
  dropdims(rv, dims=1)
end

function rmsbandwidth(s::AbstractVector; fs=framerate(s), nfft=1024, window=nothing) where T
  p = welch_pgram(s, min(nfft, length(s)); fs=inHz(fs), window=window)
  f = freq(p)
  f0 = wmean(f, power(p))
  sqrt.(wmean((f.-f0).^2, power(p)))
end

### utility functions

wmean(x::AbstractVector, w::AbstractVector) = (x'w) / sum(w)
wmean(x, w::AbstractVector) = sum(x.*w) ./ sum(w)
wmean(x, w::AbstractMatrix) = dropdims(sum(x.*w; dims=1) ./ sum(w; dims=1), dims=1)


text:
using Chain: @chain
using Dates
using JSON3: JSON3

using Stonks: JSONContent, APIResponseError, ContentParserError
using Stonks.Models:
  AssetPrice,
  AssetInfo,
  ExchangeRate,
  IncomeStatement,
  BalanceSheet,
  CashflowStatement,
  Earnings

function parse_yahoo_info(
  content::AbstractString; kwargs...
)::Union{Vector{AssetInfo},Exception}
  maybe_js = @chain begin
    content
    validate_yahoo_response
    unpack_quote_summary_response
  end
  typeof(maybe_js) <: Exception && return maybe_js
  js = maybe_js
  key_check = setdiff([:quoteType, :price], keys(js))
  if !isempty(key_check)
    return ContentParserError("Missing keys: $(join(key_check, ","))")
  end
  res = (
    assetProfile=get(js, "assetProfile", Dict()),
    quoteType=get(js, "quoteType", Dict()),
    price=get(js, "price", Dict()),
  )
  return [
    AssetInfo(;
      symbol=res.quoteType["symbol"],
      currency=res.price["currency"],
      name=get(res.quoteType, "longName", missing),
      type=get(res.quoteType, "quoteType", missing),
      exchange=get(res.quoteType, "exchange", missing),
      country=get(res.assetProfile, "country", missing),
      industry=get(res.assetProfile, "industry", missing),
      sector=get(res.assetProfile, "sector", missing),
      timezone=get(res.quoteType, "timeZoneFullName", missing),
      employees=get(res.assetProfile, "fullTimeEmployees", missing),
    ),
  ]
end

function parse_yahoo_price(
  content::AbstractString; kwargs...
)::Union{Vector{AssetPrice},Exception}
  maybe_js = validate_yahoo_response(content)
  typeof(maybe_js) <: Exception && return maybe_js
  js = maybe_js
  args = Dict(kwargs)
  from, to = get(args, :from, missing), get(args, :to, missing)
  prices = @chain js begin
    [parse_price_record(v) for (k, v) in _]
    filter(x -> isa(x, Vector{AssetPrice}), _)
  end
  if isempty(prices)
    return ContentParserError("Couldn't constrct Vector{AssetPrice} from the json object")
  end
  return @chain prices begin
    isa(from, Date) ? map(price -> filter(row -> row.date >= from, price), _) : _
    isa(to, Date) ? map(price -> filter(row -> row.date <= to, price), _) : _
    vcat(_...)
    unique
  end
end

function parse_yahoo_exchange_rate(
  content::AbstractString; kwargs...
)::Union{Vector{ExchangeRate},Exception}
  maybe_js = validate_yahoo_response(content)
  typeof(maybe_js) <: Exception && return maybe_js
  js, args = maybe_js, Dict(kwargs)
  from, to = get(args, :from, missing), get(args, :to, missing)
  rates = @chain js begin
    [parse_exchange_record(v) for (k, v) in _]
    filter(x -> isa(x, Vector{ExchangeRate}), _)
  end
  if isempty(rates)
    return ContentParserError("Couldn't constrct Vector{ExchangeRate} from the json object")
  end
  return @chain rates begin
    isa(from, Date) ? map(xrate -> filter(row -> row.date >= from, xrate), _) : _
    isa(to, Date) ? map(xrate -> filter(row -> row.date <= to, xrate), _) : _
    vcat(_...)
    unique
  end
end

function parse_yahoo_income_statement(
  content::AbstractString; kwargs...
)::Union{Vector{IncomeStatement},Exception}
  maybe_js = @chain begin
    content
    validate_yahoo_response
    unpack_quote_summary_response
  end
  typeof(maybe_js) <: Exception && return maybe_js
  js = maybe_js
  from, to = [get(kwargs, arg, missing) for arg in [:from, :to]]
  frequency = get(kwargs, :frequency, missing)
  symbol = get(kwargs, :symbol, missing)
  keys_default = [:incomeStatementHistory, :incomeStatementHistoryQuarterly]
  currency = extract_currency(js)
  keys_exp = (
    if ismissing(frequency)
      keys_default
    elseif frequency in ["annual", "annualy", "year", "yearly"]
      [:incomeStatementHistory]
    elseif frequency in ["quarter", "quarterly"]
      [:incomeStatementHistoryQuarterly]
    else
      keys_default
    end
  )
  key_check = setdiff(keys_exp, keys(js))
  if !isempty(key_check)
    return ContentParserError("Missing keys: $(join(key_check, ","))")
  end
  for key in keys_exp
    if !haskey(js[key], :incomeStatementHistory)
      return ContentParserError("Missing keys js[$key][:incomeStatementHistory]")
    end
  end
  remaps = Dict(
    :selling_general_and_administrative => :selling_general_administrative,
    :research_and_development => :research_development,
    :net_income_common_shares => :net_income_applicable_to_common_shares,
  )
  data = IncomeStatement[]
  for key in keys_exp
    js_vals = js[key][:incomeStatementHistory]
    dvals = map(x -> begin
      dval = js_to_dict(x; to_snake_case=true)
      dval[:date] = Date(unix2datetime(dval[:end_date]))
      dval[:interest_expense] = abs(get(dval, :interest_expense, missing))
      if !ismissing(currency)
        dval[:currency] = currency
      end
      dval
    end, js_vals)
    freq = key == :incomeStatementHistory ? "yearly" : "quarterly"
    is = map(
      obj -> tryparse_js(
        IncomeStatement, obj;
        fixed=Dict(:symbol => symbol, :frequency => freq),
        remaps=remaps,
      ),
      dvals,
    )
    append!(data, is)
  end
  original_len, latest_date = length(data), maximum(map(x -> x.date, data))
  res = apply_filters(data, "date"; from=from, to=to)
  if isempty(res)
    @warn """No datapoints between '$from' and '$to' after filtering.
             Original length: $original_len. Latest date: $latest_date"""
    return IncomeStatement[]
  end
  return res
end

function parse_yahoo_balance_sheet(
  content::AbstractString; kwargs...
)::Union{Vector{BalanceSheet},Exception}
  maybe_js = @chain begin
    content
    validate_yahoo_response
    unpack_quote_summary_response
  end
  typeof(maybe_js) <: Exception && return maybe_js
  js = maybe_js
  from, to = [get(kwargs, arg, missing) for arg in [:from, :to]]
  frequency = get(kwargs, :frequency, missing)
  symbol = get(kwargs, :symbol, missing)
  keys_default = [:balanceSheetHistory, :balanceSheetHistoryQuarterly]
  currency = extract_currency(js)
  keys_exp = (
    if ismissing(frequency)
      keys_default
    elseif frequency in ["annual", "annualy", "year", "yearly"]
      [:balanceSheetHistory]
    elseif frequency in ["quarter", "quarterly"]
      [:balanceSheetHistoryQuarterly]
    else
      keys_default
    end
  )
  key_check = setdiff(keys_exp, keys(js))
  if !isempty(key_check)
    return ContentParserError("Missing keys: $(join(key_check, ","))")
  end
  for key in keys_exp
    if !haskey(js[key], :balanceSheetStatements)
      return ContentParserError("Missing keys js[$key][:balanceSheetStatements]")
    end
  end
  remaps = Dict(
    :cash_and_equivalents => :cash,
    :current_net_receivables => :net_receivables,
    :current_accounts_payable => :accounts_payable,
    :other_current_liabilities => :other_current_liab,
    :total_liabilities => :total_liab,
    :total_shareholder_equity => :total_stockholder_equity,
  )
  data = BalanceSheet[]
  for key in keys_exp
    js_vals = js[key][:balanceSheetStatements]
    dvals = map(
      x -> begin
        dval = js_to_dict(x; to_snake_case=true)
        ks = [
          :intangible_assets,
          :goodwill,
          :treasury_stock,
          :other_stockholder_equity,
          :long_term_debt,
          :long_term_debt_noncurrent,
        ]
        for k in ks
          if haskey(dval, k)
            dval[k] = !isa(dval[k], Int64) ? tryparse(Int64, dval[k]) : dval[k]
          else
            dval[k] = missing
          end
        end
        dval[:date] = Date(unix2datetime(dval[:end_date]))
        dval[:currency] = currency
        dval[:intangible_assets] = dval[:intangible_assets] + dval[:goodwill]
        dval[:treasury_stock] =
          abs(dval[:treasury_stock]) - abs(dval[:other_stockholder_equity])
        dval[:long_term_debt_noncurrent] = dval[:long_term_debt]
        delete!(dval, :long_term_debt)
        [delete!(dval, k) for (k, v) in dval if ismissing(v)]
        dval
      end,
      js_vals,
    )
    freq = key == :balanceSheetHistory ? "yearly" : "quarterly"
    bs = map(
      obj -> tryparse_js(
        BalanceSheet,obj;
        fixed=Dict(:symbol => symbol, :frequency => freq),
        remaps=remaps,
      ),
      dvals,
    )
    append!(data, bs)
  end
  original_len, latest_date = length(data), maximum(map(x -> x.date, data))
  res = apply_filters(data, "date"; from=from, to=to)
  if isempty(res)
    @warn """No datapoints between '$from' and '$to' after filtering.
             Original length: $original_len. Latest date: $latest_date"""
    return BalanceSheet[]
  end
  return res
end

function parse_yahoo_cashflow_statement(
  content::AbstractString; kwargs...
)::Union{Vector{CashflowStatement},Exception}
  maybe_js = @chain begin
    content
    validate_yahoo_response
    unpack_quote_summary_response
  end
  typeof(maybe_js) <: Exception && return maybe_js
  js = maybe_js
  from, to = [get(kwargs, arg, missing) for arg in [:from, :to]]
  frequency = get(kwargs, :frequency, missing)
  symbol = get(kwargs, :symbol, missing)
  keys_default = [:cashflowStatementHistory, :cashflowStatementHistoryQuarterly]
  currency = extract_currency(js)
  keys_exp = (
    if ismissing(frequency)
      keys_default
    elseif frequency in ["annual", "annualy", "year", "yearly"]
      [first(keys_default)]
    elseif frequency in ["quarter", "quarterly"]
      [last(keys_default)]
    else
      keys_default
    end
  )
  key_check = setdiff(keys_exp, keys(js))
  if !isempty(key_check)
    return ContentParserError("Missing keys: $(join(key_check, ","))")
  end
  for key in keys_exp
    if !haskey(js[key], :cashflowStatements)
      return ContentParserError("Missing keys js[$key][:cashflowStatements]")
    end
  end
  remaps = Dict(
    :operating_cashflow => :total_cash_from_operating_activities,
    :change_operating_liabilities => :change_to_liabilities,
    :depreciation_and_amortization => :depreciation,
    :change_receivables => :change_to_account_receivables,
    :change_inventory => :change_to_inventory,
    :cashflow_investment => :total_cashflows_from_investing_activities,
    :cashflow_financing => :total_cash_from_financing_activities,
    :dividend_payout => :dividends_paid,
    :stock_repurchase => :repurchase_of_stock,
  )
  data = CashflowStatement[]
  for key in keys_exp
    js_vals = js[key][:cashflowStatements]
    dvals = map(x -> begin
      dval = js_to_dict(x; to_snake_case=true)
      ks = [
        :change_in_cash,
        :effect_of_exchange_rate,
        :capital_expenditures,
        :change_to_account_receivables,
        :change_to_inventory,
        :dividends_paid
      ]
      for k in ks
        if haskey(dval, k)
          dval[k] = !isa(dval[k], Int64) ? tryparse(Int64, dval[k]) : dval[k]
        else
          dval[k] = missing
        end
      end
      dval[:date] = Date(unix2datetime(dval[:end_date]))
      dval[:currency] = currency
      dval[:change_cash_and_equivalents] = dval[ks[1]] - dval[ks[2]]
      inv_sign = [:capital_expenditures, :change_to_account_receivables,
                  :change_to_inventory, :dividends_paid]
      [dval[k] = -dval[k] for k in inv_sign]
      dval
    end, js_vals)
    freq = key == first(keys_default) ? "yearly" : "quarterly"
    items = map(
      obj -> tryparse_js(
        CashflowStatement, obj;
        fixed=Dict(:symbol => symbol, :frequency => freq),
        remaps=remaps,
      ),
    dvals,
    )
    append!(data, items)
  end
  original_len, latest_date = length(data), maximum(map(x -> x.date, data))
  res = apply_filters(data, "date"; from=from, to=to)
  if isempty(res)
    @warn """No datapoints between '$from' and '$to' after filtering.
             Original length: $original_len. Latest date: $latest_date"""
    return CashflowStatement[]
  end
  return res
end

function parse_yahoo_earnings(
  content::AbstractString; kwargs...
)::Union{Vector{Earnings},Exception}
  maybe_js = @chain begin
    content
    validate_yahoo_response
    unpack_quote_summary_response
  end
  typeof(maybe_js) <: Exception && return maybe_js
  js = maybe_js
  from, to = [get(kwargs, arg, missing) for arg in [:from, :to]]
  frequency = get(kwargs, :frequency, missing)
  if !ismissing(frequency) && frequency in ["yearly", "annualy", "year", "annual"]
    @warn "'$frequency' frequency not available for this resource, only 'quarterly'."
    return Earnings[]
  end
  symbol = get(kwargs, :symbol, missing)
  keys_exp = [:earningsHistory]
  currency = extract_currency(js)
  key_check = setdiff(keys_exp, keys(js))
  if !isempty(key_check)
    return ContentParserError("Missing keys: $(join(key_check, ","))")
  end
  for key in keys_exp
    if !haskey(js[key], :history)
      return ContentParserError("Missing keys js[$key][:history]")
    end
  end
  remaps = Dict(
    :actual => :eps_actual,
    :estimate => :eps_estimate,
  )
  data = Earnings[]
  for key in keys_exp
    js_vals = js[key][:history]
    dvals = map(x -> begin
      dval = js_to_dict(x; to_snake_case=true)
      dval[:date] = Date(unix2datetime(dval[:quarter]))
      dval[:currency] = currency
      dval
    end, js_vals)
    items = map(
      obj -> tryparse_js(
        Earnings, obj;
        fixed=Dict(:symbol => symbol, :frequency => "quarterly"),
        remaps=remaps,
      ),
    dvals,
    )
    append!(data, items)
  end
  original_len, latest_date = length(data), maximum(map(x -> x.date, data))
  res = apply_filters(data, "date"; from=from, to=to)
  if isempty(res)
    @warn """No datapoints between '$from' and '$to' after filtering.
             Original length: $original_len. Latest date: $latest_date"""
    return Earnings[]
  end
  return res
end

function extract_currency(js::JSONContent)
  if haskey(js, :price)
    return (haskey(js[:price], :currency) ? js[:price][:currency] : missing)
  else
    return missing
  end
end

function parse_price_record(js_value::JSONContent)::Union{Vector{AssetPrice},Nothing}
  _keys = [String(k) for k in keys(js_value)]
  if !isempty(setdiff(["symbol", "timestamp", "close"], _keys))
    return nothing
  end
  isempty(js_value["timestamp"]) && return nothing
  nrows = length(js_value["timestamp"])
  ticker = js_value["symbol"]
  return [
    AssetPrice(;
      symbol=ticker,
      date=Date(unix2datetime(js_value["timestamp"][i])),
      close=Float64(js_value["close"][i]),
    ) for i in 1:nrows
  ]
end

function parse_exchange_record(js_value::JSONContent)::Union{Vector{ExchangeRate},Nothing}
  _keys = [String(k) for k in keys(js_value)]
  if !isempty(setdiff(["symbol", "timestamp", "close"], _keys))
    return nothing
  end
  isempty(js_value["timestamp"]) && return nothing
  nrows = length(js_value["timestamp"])
  base, target = @chain js_value begin
    replace(_["symbol"], "=X" => "")
    (_[1:3], _[4:length(_)])
  end
  return [
    ExchangeRate(;
      base=base,
      target=target,
      date=Date(unix2datetime(js_value["timestamp"][i])),
      rate=Float64(js_value["close"][i]),
    ) for i in 1:nrows
  ]
end

function validate_yahoo_response(content::AbstractString)::Union{JSONContent,Exception}
  maybe_js = JSON3.read(content)
  maybe_js === nothing && return error("Content could not be parsed as JSON")
  js_keys = keys(maybe_js)
  js = length(js_keys) == 1 ? maybe_js[first(js_keys)] : maybe_js
  error_idx = findfirst(x -> contains(lowercase(x), "error"), [String(k) for k in keys(js)])
  error_in_response = error_idx !== nothing ? isa(js["error"], JSON3.Object) : false
  if error_in_response
    error_msg = "Response contains an error"
    if haskey(js, :error)
      if isa(js[:error], JSON3.Object)
        error_msg = js[:error][first(keys(js[:error]))]
      elseif isa(js[:error], String)
        error_msg = js[:error]
      end
    end
    return APIResponseError(error_msg)
  end
  # In case the response has only 1 key, return the original response, not the one inside the key
  return maybe_js
end

function unpack_quote_summary_response(js::Union{JSON3.Object,Exception})
  !isa(js, JSON3.Object) && return js
  !in("quoteSummary", keys(js)) &&
    return ContentParserError("expected key 'quoteSummary' not found in API response")
  js["quoteSummary"]["error"] !== nothing &&
    return APIResponseError("API response contains error")
  res = js["quoteSummary"]["result"]
  return length(res) == 1 ? first(res) : res
end


text:
#function that returns array with the bags
#that contain the color we are searching for
function contain_color(color, bags)

   bags_containing_color = ""

   for bag in bags
      if occursin(color, bag)

        #get what is before "contain"
        m = match(r".+(?= contain)",bag)
        before = m.match

        if !occursin(color,before)

           #get what is after "contain"
           #m = match(r"(?<=contain ).+",bag)
           #after = m.match

           #append the name of bag containing color
           bags_containing_color= bags_containing_color * before
        end

      end
   end
   return bags_containing_color
end

#read file
f = open("input_day7","r")

#file has 594 lines
lines = Array{String}(undef,594)

c = 1
for line in readlines(f)
   lines[c] = line
   global c += 1
end

#search for shiny gold
first_level = contain_color("shiny gold", lines)
let first_colors = split(first_level," bags")
println(first_colors)
total = length(first_colors) - 1 #account for the ""
println(total)

old_colors = first_colors

i = 1
while i < 10

   second_level_ts = ""

   #search for the ones containing the prev ones
   for color in first_colors
      if color != ""
         second_level = contain_color(color,lines)
         second_level_ts *= second_level
      end
   end

   second_colors_ts = split(second_level_ts," bags")
   second_colors =unique([s for s in second_colors_ts if !(s in old_colors)])
   println(second_colors)
   total = total + length(second_colors)
   println(total)
   first_colors = second_colors
   old_colors = vcat(old_colors,second_colors)

   i+= 1
end

end #let


text:
using IteratedMaps
using Base.Test

# A logistic map with parameter 2 converges to 0.5
lm = LogisticMap(2)

for i in take(lm(0.5),100)
    @test i == 0.5
end

# Test the alternative syntax for map generation
for i in 1:100
    @test collect(take(lm(0.2),100))[i] == collect(take(LogisticMap(2,0.2),100))[i]
end

# Check the Henon map calculates correctly
@test collect(take(HenonMap([1.4,0.3],[0.2,0.3]),2)) == Any[[0.2,0.3],[1.244,0.06]]

# Same with the Duffing map
@test collect(take(DuffingMap([2.75,0.2],[1.0,0.0]),2)) == Any[[1.0,0.0],[0.0,-0.2]]

# And the Tent map

@test collect(take(TentMap(2,0.2),2)) == [0.2,0.4]

max_stars_repo_name: wkearn/IteratedMaps

max_issues_repo_name: wkearn/IteratedMaps

max_forks_repo_name: wkearn/IteratedMaps


text:
"""
```
Model990{T} <: AbstractModel{T}
```

The `Model990` type defines the structure of the FRBNY DSGE model.

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model parameters.

* `steady_state::Vector{AbstractParameter}`: Model steady-state values, computed as a function of elements of
  `parameters`.

* `keys::Dict{Symbol,Int}`: Maps human-readable names for all model parameters and
  steady-states to their indices in `parameters` and `steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readible names to row and column
indices in the matrix representations of of the measurement equation and equilibrium
conditions.

* `endogenous_states::Dict{Symbol,Int}`: Maps each state to a column in the measurement and
  equilibrium condition matrices.

* `exogenous_shocks::Dict{Symbol,Int}`: Maps each shock to a column in the measurement and
  equilibrium condition matrices.

* `expected_shocks::Dict{Symbol,Int}`: Maps each expected shock to a column in the
  measurement and equilibrium condition matrices.

* `equilibrium_conditions::Dict{Symbol,Int}`: Maps each equlibrium condition to a row in the
  model's equilibrium condition matrices.

* `endogenous_states_augmented::Dict{Symbol,Int}`: Maps lagged states to their columns in
  the measurement and equilibrium condition equations. These are added after Gensys solves the
  model.

* `observables::Dict{Symbol,Int}`: Maps each observable to a row in the model's measurement
  equation matrices.

#### Model Specifications and Settings

* `spec::AbstractString`: The model specification identifier, "m990", cached here for
  filepath computation.

* `subspec::AbstractString`: The model subspecification number, indicating that some
  parameters from the original model spec ("ss0") are initialized differently. Cached here for
  filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect computation without changing
  the economic or mathematical setup of the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`, settings from
  `m.test_settings` are used in place of those in `m.settings`.

* `data_series::Dict{Symbol,Vector{Symbol}}`: A dictionary that
  stores data sources (keys) and lists of series mnemonics
  (values). DSGE.jl will fetch data from the Federal Reserve Bank of
  St. Louis's FRED database; all other data must be downloaded by the
  user. See `load_data` for further details.
"""
type Model990{T} <: AbstractModel{T}
    parameters::ParameterVector{T}                  # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                # model steady-state values
    keys::Dict{Symbol,Int}                          # human-readable names for all the model
                                                    # parameters and steady-n_states

    endogenous_states::Dict{Symbol,Int}             # these fields used to create matrices in the
    exogenous_shocks::Dict{Symbol,Int}              # measurement and equilibrium condition equations.
    expected_shocks::Dict{Symbol,Int}               #
    equilibrium_conditions::Dict{Symbol,Int}        #
    endogenous_states_augmented::Dict{Symbol,Int}   #
    observables::Dict{Symbol,Int}                   #

    spec::ASCIIString                               # Model specification number (eg "m990")
    subspec::ASCIIString                            # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                  # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}             # Settings/flags for testing mode
    rng::MersenneTwister                            # Random number generator
    testing::Bool                                   # Whether we are in testing mode or not

    data_series::Dict{Symbol,Vector{Symbol}}       # Keys = data sources, values = vector of series mnemonics
    data_transforms::OrderedDict{Symbol,Function}  # functions to transform raw data into input matrix
end

description(m::Model990) = "FRBNY DSGE Model m990, $(m.subspec)"

"""
`init_model_indices!(m::Model990)`

Arguments:
`m:: Model990`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::Model990)
    # Endogenous states
    endogenous_states = [[
        :y_t, :c_t, :i_t, :qk_t, :k_t, :kbar_t, :u_t, :rk_t, :Rtil_k_t, :n_t, :mc_t,
        :π_t, :μ_ω_t, :w_t, :L_t, :R_t, :g_t, :b_t, :μ_t, :z_t, :λ_f_t, :λ_f_t1,
        :λ_w_t, :λ_w_t1, :rm_t, :σ_ω_t, :μ_e_t, :γ_t, :π_star_t, :Ec_t, :Eqk_t, :Ei_t,
        :Eπ_t, :EL_t, :Erk_t, :Ew_t, :ERtil_k_t, :y_f_t, :c_f_t, :i_f_t, :qk_f_t, :k_f_t,
        :kbar_f_t, :u_f_t, :rk_f_t, :w_f_t, :L_f_t, :r_f_t, :Ec_f_t, :Eqk_f_t, :Ei_f_t,
        :EL_f_t, :Erk_f_t, :ztil_t, :π_t1, :π_t2, :π_a_t, :R_t1, :zp_t, :Ez_t];
        [symbol("rm_tl$i") for i = 1:n_anticipated_shocks(m)]]

    # Exogenous shocks
    exogenous_shocks = [[
        :g_sh, :b_sh, :μ_sh, :z_sh, :λ_f_sh, :λ_w_sh, :rm_sh, :σ_ω_sh, :μ_e_sh,
        :γ_sh, :π_star_sh, :lr_sh, :zp_sh, :tfp_sh, :gdpdef_sh, :corepce_sh];
        [symbol("rm_shl$i") for i = 1:n_anticipated_shocks(m)]]

    # Expectations shocks
    expected_shocks = [
        :Ec_sh, :Eqk_sh, :Ei_sh, :Eπ_sh, :EL_sh, :Erk_sh, :Ew_sh, :ERktil_sh, :Ec_f_sh,
        :Eqk_f_sh, :Ei_f_sh, :EL_f_sh, :Erk_f_sh]

    # Equilibrium conditions
    equilibrium_conditions = [[
        :eq_euler, :eq_inv, :eq_capval, :eq_spread, :eq_nevol, :eq_output, :eq_caputl, :eq_capsrv, :eq_capev,
        :eq_mkupp, :eq_phlps, :eq_caprnt, :eq_msub, :eq_wage, :eq_mp, :eq_res, :eq_g, :eq_b, :eq_μ, :eq_z,
        :eq_λ_f, :eq_λ_w, :eq_rm, :eq_σ_ω, :eq_μ_e, :eq_γ, :eq_λ_f1, :eq_λ_w1, :eq_Ec,
        :eq_Eqk, :eq_Ei, :eq_Eπ, :eq_EL, :eq_Erk, :eq_Ew, :eq_ERktil, :eq_euler_f, :eq_inv_f,
        :eq_capval_f, :eq_output_f, :eq_caputl_f, :eq_capsrv_f, :eq_capev_f, :eq_mkupp_f,
        :eq_caprnt_f, :eq_msub_f, :eq_res_f, :eq_Ec_f, :eq_Eqk_f, :eq_Ei_f, :eq_EL_f, :eq_Erk_f,
        :eq_ztil, :eq_π_star, :eq_π1, :eq_π2, :eq_π_a, :eq_Rt1, :eq_zp, :eq_Ez];
        [symbol("eq_rml$i") for i=1:n_anticipated_shocks(m)]]

    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = [
        :y_t1, :c_t1, :i_t1, :w_t1, :π_t1, :L_t1, :Et_π_t, :lr_t, :tfp_t, :e_gdpdef_t,
        :e_corepce_t, :u_t1]

    # Measurement equation observables
    observables = [[
        :obs_gdp,              # quarterly output growth
        :obs_hours,            # aggregate hours growth
        :obs_wages,            # real wage growth
        :obs_gdpdeflator,      # inflation (GDP deflator)
        :obs_corepce,          # inflation (core PCE)
        :obs_nominalrate,      # nominal interest rate
        :obs_consumption,      # consumption growth
        :obs_investment,       # investment growth
        :obs_spread,           # spreads
        :obs_longinflation,    # 10-year inflation expectation
        :obs_longrate,         # long-term rate
        :obs_tfp];             # total factor productivity
        [symbol("obs_nominalrate$i") for i=1:n_anticipated_shocks(m)]] # compounded nominal rates

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

function Model990(subspec::AbstractString="ss2")

    # Model-specific specifications
    spec               = split(basename(@__FILE__),'.')[1]
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister()
    testing            = false

    # Set up data sources and series
    fred_series        = [:GDP, :GDPCTPI, :PCE, :FPI, :CNP16OV, :CE16OV, :PRS85006013,
                          :UNRATE, :AWHNONAG, :DFF, :BAA, :GS10, :PRS85006063, :CES0500000030, :CLF16OV,
                          :PCEPILFE, :COMPNFB, :THREEFYTP10]
    spf_series         = [:ASACX10]
    fernald_series     = [:TFPJQ, :TFPKQ]
    longrate_series    = [:FYCZZA]
    # ois data taken care of in load_data

    data_series = Dict{Symbol,Vector{Symbol}}(:fred => fred_series, :spf => spf_series,
                                              :fernald => fernald_series, :longrate => longrate_series)

    # set up data transformations
    data_transforms = OrderedDict{Symbol,Function}()

    # initialize empty model
    m = Model990{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(), Dict{Symbol,Int}(),

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

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            data_series,
            data_transforms)

    # Set settings
    settings_m990!(m)
    default_test_settings!(m)

    # Set data transformations
    init_data_transforms!(m)

    # Initialize parameters
    m <= parameter(:α,      0.1596, (1e-5, 0.999), (1e-5, 0.999),   DSGE.SquareRoot(),     Normal(0.30, 0.05),         fixed=false,
                   description="α: Capital elasticity in the intermediate goods sector's production function (also known as the capital share).",
                   tex_label="\\alpha")

    m <= parameter(:ζ_p,   0.8940, (1e-5, 0.999), (1e-5, 0.999),   DSGE.SquareRoot(),     BetaAlt(0.5, 0.1),          fixed=false,
                   description="ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")

    m <= parameter(:ι_p,   0.1865, (1e-5, 0.999), (1e-5, 0.999),   DSGE.SquareRoot(),     BetaAlt(0.5, 0.15),         fixed=false,
                   description="ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",
                   tex_label="\\iota_p")

    m <= parameter(:δ,      0.025,  fixed=true,
                   description="δ: The capital depreciation rate.", tex_label="\\delta" )

    m <= parameter(:Upsilon,  1.000,  (0., 10.),     (1e-5, 0.),      DSGE.Exponential(),    GammaAlt(1., 0.5),          fixed=true,
                   description="Υ: The trend evolution of the price of investment goods relative to consumption goods. Set equal to 1.",
                   tex_label="\\mathcal{\\Upsilon}")

    m <= parameter(:Φ,   1.1066, (1., 10.),     (1.00, 10.00),   DSGE.Exponential(),    Normal(1.25, 0.12),         fixed=false,
                   description="Φ: Fixed costs.",
                   tex_label="\\Phi")

    m <= parameter(:S′′,       2.7314, (-15., 15.),   (-15., 15.),     DSGE.Untransformed(),  Normal(4., 1.5),            fixed=false,
                   description="S'': The second derivative of households' cost of adjusting investment.", tex_label="S\\prime\\prime")

    m <= parameter(:h,        0.5347, (1e-5, 0.999), (1e-5, 0.999),   DSGE.SquareRoot(),     BetaAlt(0.7, 0.1),          fixed=false,
                   description="h: Consumption habit persistence.", tex_label="h")

    m <= parameter(:ppsi,     0.6862, (1e-5, 0.999), (1e-5, 0.999),   DSGE.SquareRoot(),     BetaAlt(0.5, 0.15),         fixed=false,
                   description="ppsi: Utilization costs.", tex_label="ppsi")

    m <= parameter(:ν_l,     2.5975, (1e-5, 10.),   (1e-5, 10.),     DSGE.Exponential(),    Normal(2, 0.75),            fixed=false,
                   description="ν_l: The coefficient of relative risk aversion on the labor
                   term of households' utility function.", tex_label="\\nu_l")

    m <= parameter(:ζ_w,   0.9291, (1e-5, 0.999), (1e-5, 0.999),   DSGE.SquareRoot(),     BetaAlt(0.5, 0.1),          fixed=false,
                   description="ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label="\\zeta_w")

    m <= parameter(:ι_w,   0.2992, (1e-5, 0.999), (1e-5, 0.999),   DSGE.SquareRoot(),     BetaAlt(0.5, 0.15),         fixed=false,
                   description="ι_w: No description available.",
                   tex_label="\\iota_w")

    m <= parameter(:λ_w,      1.5000,                                                                               fixed=true,
                   description="λ_w: The wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label="\\lambda_w")

    m <= parameter(:β,      0.1402, (1e-5, 10.),   (1e-5, 10.),     DSGE.Exponential(),    GammaAlt(0.25, 0.1),        fixed=false,  scaling = x -> 1/(1 + x/100),
                   description="β: Discount rate.",
                   tex_label="\\beta ")

    m <= parameter(:ψ1,     1.3679, (1e-5, 10.),   (1e-5, 10.00),   DSGE.Exponential(),    Normal(1.5, 0.25),          fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ψ2,     0.0388, (-0.5, 0.5),   (-0.5, 0.5),     DSGE.Untransformed(),  Normal(0.12, 0.05),         fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3,     0.2464, (-0.5, 0.5),   (-0.5, 0.5),     DSGE.Untransformed(),  Normal(0.12, 0.05),         fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")

    m <= parameter(:π_star,   0.5000, (1e-5, 10.),   (1e-5, 10.),     DSGE.Exponential(),    GammaAlt(0.75, 0.4),        fixed=true,  scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")

    m <= parameter(:σ_c,   0.8719, (1e-5, 10.),   (1e-5, 10.),     DSGE.Exponential(),    Normal(1.5, 0.37),          fixed=false,
                   description="σ_c: No description available.",
                   tex_label="\\sigma_{c}")

    m <= parameter(:ρ,      0.7126, (1e-5, 0.999), (1e-5, 0.999),   DSGE.SquareRoot(),     BetaAlt(0.75, 0.10),        fixed=false,
                   description="ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho")

    m <= parameter(:ϵ_p,     10.000,                                                                               fixed=true,
                   description="ϵ_p: No description available.",
                   tex_label="\\varepsilon_{p}")

    m <= parameter(:ϵ_w,     10.000,                                                                               fixed=true,
                   description="ϵ_w: No description available.",
                   tex_label="\\varepsilon_{w}")

    # financial frictions parameters
    m <= parameter(:Fω,      0.0300, (1e-5, 0.99999), (1e-5, 0.99),  DSGE.SquareRoot(),    BetaAlt(0.03, 0.01),         fixed=true,    scaling = x -> 1 - (1-x)^0.25,
                   description="F(ω): The cumulative distribution function of ω (idiosyncratic iid shock that increases or decreases entrepreneurs' capital).",
                   tex_label="F(\\omega)")

    m <= parameter(:spr,     1.7444, (0., 100.),      (1e-5, 0.),    DSGE.Exponential(),   GammaAlt(2., 0.1),           fixed=false,  scaling = x -> (1 + x/100)^0.25,
                   description="spr_*: Steady-state spreads.",
                   tex_label="spr_*")

    m <= parameter(:ζ_spb, 0.0559, (1e-5, 0.99999), (1e-5, 0.99),  DSGE.SquareRoot(),    BetaAlt(0.05, 0.005),        fixed=false,
                   description="ζ_spb: The elasticity of the expected exess return on capital (or 'spread') with respect to leverage.",
                   tex_label="\\zeta_{spb}")

    m <= parameter(:γ_star, 0.9900, (1e-5, 0.99999), (1e-5, 0.99),  DSGE.SquareRoot(),    BetaAlt(0.99, 0.002),        fixed=true,
                   description="γ_star: No description available.",
                   tex_label="\\gamma_*")

    # exogenous processes - level
    m <= parameter(:γ,      0.3673, (-5.0, 5.0),     (-5., 5.),     DSGE.Untransformed(), Normal(0.4, 0.1),            fixed=false, scaling = x -> x/100,
                   description="γ: The log of the steady-state growth rate of technology.",
                   tex_label="\\gamma")

    m <= parameter(:Lmean,  -45.9364, (-1000., 1000.), (-1e3, 1e3),   DSGE.Untransformed(), Normal(-45, 5),              fixed=false,
                   description="Lmean: No description available.",
                   tex_label="Lmean")

    m <= parameter(:g_star,    0.1800,                                                                               fixed=true,
                   description="g_star: No description available.",
                   tex_label="g_*")

    # exogenous processes - autocorrelation
    m <= parameter(:ρ_g,      0.9863, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label="\\rho_g")

    m <= parameter(:ρ_b,      0.9410, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label="\\rho_b")

    m <= parameter(:ρ_μ,     0.8735, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label="\\rho_{\\mu}")

    m <= parameter(:ρ_z,      0.9446, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")

    m <= parameter(:ρ_λ_f,    0.8827, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label="\\rho_{\\lambda_f}")

    m <= parameter(:ρ_λ_w,    0.3884, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label="\\rho_{\\lambda_w}")

    #monetary policy shock - see eqcond
    m <= parameter(:ρ_rm,     0.2135, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label="\\rho_{rm}")

    m <= parameter(:ρ_σ_w,   0.9898, (1e-5, 0.99999), (1e-5, 0.99),  DSGE.SquareRoot(),    BetaAlt(0.75, 0.15),         fixed=false,
                   description="ρ_σ_w: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with mean ρ_σ_w. Innovations to the process are called _spread shocks_.",
                   tex_label="\\rho_{\\sigma_\\omega}")

    m <= parameter(:ρ_μ_e,    0.7500, (1e-5, 0.99999), (1e-5, 0.99),  DSGE.SquareRoot(),    BetaAlt(0.75, 0.15),         fixed=true,
                   description="ρ_μ_e: Verification costs are a fraction μ_e of the amount the bank extracts from an entrepreneur in case of bankruptcy???? This doesn't seem right because μ_e isn't a process (p12 of PDF)",
                   tex_label="\\rho_{\\mu_e}")

    m <= parameter(:ρ_γ,   0.7500, (1e-5, 0.99999), (1e-5, 0.99),  DSGE.SquareRoot(), BetaAlt(0.75, 0.15),         fixed=true,  description="ρ_γ: Autocorrelation coefficient on the γ shock.",              tex_label="\\rho_{\\gamma}")
    m <= parameter(:ρ_π_star,   0.9900, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(), BetaAlt(0.5, 0.2),           fixed=true,  description="ρ_π_star: Autocorrelation coefficient on the π_star shock.",  tex_label="\\rho_{\\pi^*}")
    m <= parameter(:ρ_lr,     0.6936, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false, description="ρ_lr: Autocorrelation coefficient on the lr shock.", tex_label="\\rho_{lr}")
    m <= parameter(:ρ_z_p,     0.8910, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false, description="ρ_z_p: Autocorrelation coefficient on the z_p shock.", tex_label="\\rho_{z^p}")
    m <= parameter(:ρ_tfp,    0.1953, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false, description="ρ_tfp: No description available.",            tex_label="\\rho_{tfp}")
    m <= parameter(:ρ_gdpdef, 0.5379, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false, description="ρ_gdpdef: GDP deflator.",                            tex_label="\\rho_{gdpdef}")
    m <= parameter(:ρ_corepce,    0.2320, (1e-5, 0.999),   (1e-5, 0.999), DSGE.SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false, description="ρ_corepce: No description available.",            tex_label="\\rho_{corepce}")

    # exogenous processes - standard deviation
    m <= parameter(:σ_g,      2.5230, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_b,      0.0292, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_b: No description available.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_μ,     0.4559, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_z,      0.6742, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_z: No description available.",
                   tex_label="\\sigma_{z}")

    m <= parameter(:σ_λ_f,    0.1314, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good.  Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w,    0.3864, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_λ_w: No description available.",
                   tex_label="\\sigma_{\\lambda_w}")

    m <= parameter(:σ_r_m,     0.2380, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_r_m: No description available.",
                   tex_label="\\sigma_{rm}")

    m <= parameter(:σ_σ_ω,   0.0428, (1e-7,100.),     (1e-5, 0.),    DSGE.Exponential(),   RootInverseGamma(4., 0.05),  fixed=false,
                   description="σ_σ_ω: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with standard deviation σ_σ_ω.",
                   tex_label="\\sigma_{\\sigma_\\omega}")

    m <= parameter(:σ_μ_e,    0.0000, (1e-7,100.),     (1e-5, 0.),    DSGE.Exponential(), RootInverseGamma(4., 0.05),  fixed=true,  description="σ_μ_e: Standard deviation of the μ_e shock.",         tex_label="\\sigma_{\\mu_e}")
    m <= parameter(:σ_γ,   0.0000, (1e-7,100.),     (1e-5, 0.),    DSGE.Exponential(),   RootInverseGamma(4., 0.01),  fixed=true,  description="σ_γ: Standard deviation of the γ shock.",              tex_label="\\sigma_{\\gamma}")
    m <= parameter(:σ_π_star,   0.0269, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(), RootInverseGamma(6., 0.03),  fixed=false, description="σ_π_star: Standard deviation of the π_star shock.", tex_label="\\sigma_{\\pi^*}")
    m <= parameter(:σ_lr,     0.1766, (1e-8,10.),      (1e-8, 5.),    DSGE.Exponential(), RootInverseGamma(2., 0.75),  fixed=false, description="σ_lr: Standard deviation of the lr shock.",           tex_label="\\sigma_{lr}")
    m <= parameter(:σ_z_p,     0.1662, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false, description="σ_z_p: Standard deviation of the z_p shock.",      tex_label="\\sigma_{z^p}")
    m <= parameter(:σ_tfp,    0.9391, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false, description="σ_tfp: Standard deviation of the σ_tfp shock.",     tex_label="\\sigma_{tfp}")
    m <= parameter(:σ_gdpdef, 0.1575, (1e-8, 5.), (1e-8, 5.),DSGE.Exponential(),
         RootInverseGamma(2., 0.10), fixed=false,
         description="σ_gdpdef: Standard deviation of the gdpdef shock.",
         tex_label="\\sigma_{gdpdef}")

    m <= parameter(:σ_corepce, 0.0999, (1e-8, 5.),(1e-8, 5.),DSGE.Exponential(),RootInverseGamma(2., 0.10),
                   fixed=false,
                   description="σ_corepce: Standard deviation of the pce shock.",
                   tex_label="\\sigma_{corepce}")

    # standard deviations of the anticipated policy shocks
    for i = 1:n_anticipated_shocks_padding(m)
        if i < 13
            m <= parameter(symbol("σ_r_m$i"), .2, (1e-7, 100.), (1e-5, 0.), DSGE.Exponential(),
                           RootInverseGamma(4., .2), fixed=false,
                           description="σ_r_m$i: Standard deviation of the $i-period-ahead anticipated policy shock.",
                           tex_label=@sprintf("\\sigma_{ant%d}",i))
        else
            m <= parameter(symbol("σ_r_m$i"), .0, (1e-7, 100.), (1e-5, 0.),
                           DSGE.Exponential(), RootInverseGamma(4., .2), fixed=true,
                           description="σ_r_m$i: Standard deviation of the $i-period-ahead anticipated policy shock.",
                           tex_label=@sprintf("\\sigma_{ant%d}",i))
        end
    end

    m <= parameter(:η_gz,       0.8400, (1e-5, 0.999), (1e-5, 0.999), DSGE.SquareRoot(),
                   BetaAlt(0.50, 0.20), fixed=false,
                   description="η_gz: Correlate g and z shocks.",
                   tex_label="\\eta_{gz}")

    m <= parameter(:η_λ_f,      0.7892, (1e-5, 0.999), (1e-5, 0.999), DSGE.SquareRoot(),
                   BetaAlt(0.50, 0.20),         fixed=false,
                   description="η_λ_f: Moving average component in the price markup shock.",
                   tex_label="\\eta_{\\lambda_f}")

    m <= parameter(:η_λ_w,      0.4226, (1e-5, 0.999), (1e-5, 0.999), DSGE.SquareRoot(),
                   BetaAlt(0.50, 0.20),         fixed=false,
                   description="η_λ_w: Moving average component in the wage markup shock.",
                   tex_label="\\eta_{\\lambda_w}")

    m <= parameter(:Iendoα, 0.0000, (0.000, 1.000), (0., 0.), DSGE.Untransformed(),
                   BetaAlt(0.50, 0.20), fixed=true,
                   description="Iendoα: Indicates whether to use the model's endogenous α in the capacity utilization adjustment of total factor productivity.",
                   tex_label="I\\{\\alpha^{model}\\}")

    m <= parameter(:Γ_gdpdef,  1.0354, (-10., 10.), (-10., -10.),  DSGE.Untransformed(),
                   Normal(1.00, 2.), fixed=false,
                   description="Γ_gdpdef: No description available.",
                   tex_label="\\Gamma_{gdpdef}")

    m <= parameter(:δ_gdpdef,   0.0181, (-9.1, 9.1), (-10., -10.),  DSGE.Untransformed(),
                   Normal(0.00, 2.),            fixed=false,
                   description="δ_gdpdef: No description available.",
                   tex_label="\\delta_{gdpdef}")

    # steady states
    m <= SteadyStateParameter(:z_star,  NaN, description="No description available.", tex_label="\\z_*")
    m <= SteadyStateParameter(:rstar,   NaN, description="No description available.", tex_label="\\r_*")
    m <= SteadyStateParameter(:Rstarn,  NaN, description="No description available.", tex_label="\\R_*_n")
    m <= SteadyStateParameter(:r_k_star,  NaN, description="No description available.", tex_label="\\r^k_*")
    m <= SteadyStateParameter(:wstar,   NaN, description="No description available.", tex_label="\\w_*")
    m <= SteadyStateParameter(:Lstar,   NaN, description="No description available.", tex_label="\\L_*")
    m <= SteadyStateParameter(:kstar,   NaN, description="Effective capital that households rent to firms in the steady state.", tex_label="\\k_*")
    m <= SteadyStateParameter(:kbarstar, NaN, description="Total capital owned by households in the steady state.", tex_label="\\bar{k}_*")
    m <= SteadyStateParameter(:istar,  NaN, description="Detrended steady-state investment", tex_label="\\i_*")
    m <= SteadyStateParameter(:ystar,  NaN, description="No description available.", tex_label="\\y_*")
    m <= SteadyStateParameter(:cstar,  NaN, description="No description available.", tex_label="\\c_*")
    m <= SteadyStateParameter(:wl_c,   NaN, description="No description available.", tex_label="\\wl_c")
    m <= SteadyStateParameter(:nstar,  NaN, description="No description available.", tex_label="\\n_*")
    m <= SteadyStateParameter(:vstar,  NaN, description="No description available.", tex_label="\\v_*")
    m <= SteadyStateParameter(:ζ_spσ_ω,  NaN, description="No description available.", tex_label="\\zeta_{sp_{\\sigma_\\omega}}")
    m <= SteadyStateParameter(:ζ_spμ_e,   NaN, description="No description available.", tex_label="\\zeta_{sp_{\\mu_e}}")
    m <= SteadyStateParameter(:ζ_nRk,     NaN, description="No description available.", tex_label="\\zeta_{n_R_k}")
    m <= SteadyStateParameter(:ζ_nR,      NaN, description="No description available.", tex_label="\\zeta_{n_R}")
    m <= SteadyStateParameter(:ζ_nqk,     NaN, description="No description available.", tex_label="\\zeta_{n_q_k}")
    m <= SteadyStateParameter(:ζ_nn,      NaN, description="No description available.", tex_label="\\zeta_{nn}")
    m <= SteadyStateParameter(:ζ_nμ_e,    NaN, description="No description available.", tex_label="\\zeta_{n_{\\mu_e}}")
    m <= SteadyStateParameter(:ζ_nσ_ω,   NaN, description="No description available.", tex_label="\\zeta_{n_{\\sigma_\\omega}}")

    init_model_indices!(m)
    init_subspec!(m)
    steadystate!(m)
    return m
end

# functions that are used to compute financial frictions
# steady-state values from parameter values
@inline function ζ_spb_fn(z, σ, spr)
    zetaratio = ζ_bω_fn(z, σ, spr)/ζ_zω_fn(z, σ, spr)
    nk = nk_fn(z, σ, spr)
    return -zetaratio/(1-zetaratio)*nk/(1-nk)
end

@inline function ζ_bω_fn(z, σ, spr)
    nk          = nk_fn(z, σ, spr)
    μstar       = μ_fn(z, σ, spr)
    ω_star       = ω_fn(z, σ)
    Γstar       = Γ_fn(z, σ)
    Gstar       = G_fn(z, σ)
    dΓ_dω_star   = dΓ_dω_fn(z)
    dG_dω_star   = dG_dω_fn(z, σ)
    d2Γ_dω2star = d2Γ_dω2_fn(z, σ)
    d2G_dω2star = d2G_dω2_fn(z, σ)
    return ω_star*μstar*nk*(d2Γ_dω2star*dG_dω_star - d2G_dω2star*dΓ_dω_star)/
        (dΓ_dω_star - μstar*dG_dω_star)^2/spr/(1 - Γstar + dΓ_dω_star*(Γstar - μstar*Gstar)/
            (dΓ_dω_star - μstar*dG_dω_star))
end

@inline function ζ_zω_fn(z, σ, spr)
    μstar = μ_fn(z, σ, spr)
    return ω_fn(z, σ)*(dΓ_dω_fn(z) - μstar*dG_dω_fn(z, σ))/
        (Γ_fn(z, σ) - μstar*G_fn(z, σ))
end

nk_fn(z, σ, spr) = 1 - (Γ_fn(z, σ) - μ_fn(z, σ, spr)*G_fn(z, σ))*spr
μ_fn(z, σ, spr)  =
    (1 - 1/spr)/(dG_dω_fn(z, σ)/dΓ_dω_fn(z)*(1 - Γ_fn(z, σ)) + G_fn(z, σ))
ω_fn(z, σ)        = exp(σ*z - σ^2/2)
G_fn(z, σ)        = cdf(Normal(), z-σ)
Γ_fn(z, σ)        = ω_fn(z, σ)*(1 - cdf(Normal(), z)) + cdf(Normal(), z-σ)
dG_dω_fn(z, σ)    = pdf(Normal(), z)/σ
d2G_dω2_fn(z, σ)  = -z*pdf(Normal(), z)/ω_fn(z, σ)/σ^2
dΓ_dω_fn(z)       = 1 - cdf(Normal(), z)
d2Γ_dω2_fn(z, σ)  = -pdf(Normal(), z)/ω_fn(z, σ)/σ
dG_dσ_fn(z, σ)    = -z*pdf(Normal(), z-σ)/σ
d2G_dωdσ_fn(z, σ) = -pdf(Normal(), z)*(1 - z*(z-σ))/σ^2
dΓ_dσ_fn(z, σ)    = -pdf(Normal(), z-σ)
d2Γ_dωdσ_fn(z, σ) = (z/σ-1)*pdf(Normal(), z)

"""
```
steadystate!(m::Model990)
```

Calculates the model's steady-state values. `steadystate!(m)` must be called whenever
the parameters of `m` are updated.
"""
function steadystate!(m::Model990)
    SIGWSTAR_ZERO = 0.5

    m[:z_star]    = log(1+m[:γ]) + m[:α]/(1-m[:α])*log(m[:Upsilon])
    m[:rstar]    = exp(m[:σ_c]*m[:z_star]) / m[:β]
    m[:Rstarn]   = 100*(m[:rstar]*m[:π_star] - 1)
    m[:r_k_star]   = m[:spr]*m[:rstar]*m[:Upsilon] - (1-m[:δ])
    m[:wstar]    = (m[:α]^m[:α] * (1-m[:α])^(1-m[:α]) * m[:r_k_star]^(-m[:α]) / m[:Φ])^(1/(1-m[:α]))
    m[:Lstar]    = 1.
    m[:kstar]    = (m[:α]/(1-m[:α])) * m[:wstar] * m[:Lstar] / m[:r_k_star]
    m[:kbarstar] = m[:kstar] * (1+m[:γ]) * m[:Upsilon]^(1 / (1-m[:α]))
    m[:istar]    = m[:kbarstar] * (1-((1-m[:δ])/((1+m[:γ]) * m[:Upsilon]^(1/(1-m[:α])))))
    m[:ystar]    = (m[:kstar]^m[:α]) * (m[:Lstar]^(1-m[:α])) / m[:Φ]
    m[:cstar]    = (1-m[:g_star])*m[:ystar] - m[:istar]
    m[:wl_c]     = (m[:wstar]*m[:Lstar])/(m[:cstar]*m[:λ_w])

    # FINANCIAL FRICTIONS ADDITIONS
    # solve for σ_ω_star and zω_star
    zω_star = quantile(Normal(), m[:Fω].scaledvalue)
    σ_ω_star = SIGWSTAR_ZERO
    try
        σ_ω_star = fzero(sigma -> ζ_spb_fn(zω_star, sigma, m[:spr]) - m[:ζ_spb], 0.5)
    catch
        σ_ω_star = SIGWSTAR_ZERO
    end
    # evaluate ωbarstar
    ωbarstar = ω_fn(zω_star, σ_ω_star)

    # evaluate all BGG function elasticities
    Gstar                   = G_fn(zω_star, σ_ω_star)
    Γstar               = Γ_fn(zω_star, σ_ω_star)
    dGdω_star            = dG_dω_fn(zω_star, σ_ω_star)
    d2Gdω2star          = d2G_dω2_fn(zω_star, σ_ω_star)
    dΓdω_star        = dΓ_dω_fn(zω_star)
    d2Γdω2star      = d2Γ_dω2_fn(zω_star, σ_ω_star)
    dGdσstar            = dG_dσ_fn(zω_star, σ_ω_star)
    d2Gdωdσstar     = d2G_dωdσ_fn(zω_star, σ_ω_star)
    dΓdσstar        = dΓ_dσ_fn(zω_star, σ_ω_star)
    d2Γdωdσstar = d2Γ_dωdσ_fn(zω_star, σ_ω_star)

    # evaluate μ, nk, and Rhostar
    μ_estar       = μ_fn(zω_star, σ_ω_star, m[:spr])
    nkstar        = nk_fn(zω_star, σ_ω_star, m[:spr])
    Rhostar       = 1/nkstar - 1

    # evaluate wekstar and vkstar
    wekstar       = (1-m[:γ_star]/m[:β])*nkstar - m[:γ_star]/m[:β]*(m[:spr]*(1-μ_estar*Gstar) - 1)
    vkstar        = (nkstar-wekstar)/m[:γ_star]

    # evaluate nstar and vstar
    m[:nstar]       = nkstar*m[:kstar]
    m[:vstar]       = vkstar*m[:kstar]

    # a couple of combinations
    ΓμG      = Γstar - μ_estar*Gstar
    ΓμGprime = dΓdω_star - μ_estar*dGdω_star

    # elasticities wrt ωbar
    ζ_bw       = ζ_bω_fn(zω_star, σ_ω_star, m[:spr])
    ζ_zw       = ζ_zω_fn(zω_star, σ_ω_star, m[:spr])
    ζ_bw_zw    = ζ_bw/ζ_zw

    # elasticities wrt σ_ω
    ζ_bσ_ω      = σ_ω_star * (((1 - μ_estar*dGdσstar/dΓdσstar) /
        (1 - μ_estar*dGdω_star/dΓdω_star) - 1)*dΓdσstar*m[:spr] + μ_estar*nkstar*
            (dGdω_star*d2Γdωdσstar - dΓdω_star*d2Gdωdσstar)/ΓμGprime^2) /
                ((1 - Γstar)*m[:spr] + dΓdω_star/ΓμGprime*(1-nkstar))
    ζ_zσ_ω      = σ_ω_star * (dΓdσstar - μ_estar*dGdσstar) / ΓμG
    m[:ζ_spσ_ω] = (ζ_bw_zw*ζ_zσ_ω - ζ_bσ_ω) / (1-ζ_bw_zw)

    # elasticities wrt μ_e
    ζ_bμ_e      = μ_estar * (nkstar*dΓdω_star*dGdω_star/ΓμGprime+dΓdω_star*Gstar*m[:spr]) /
        ((1-Γstar)*ΓμGprime*m[:spr] + dΓdω_star*(1-nkstar))
    ζ_zμ_e      = -μ_estar*Gstar/ΓμG
    m[:ζ_spμ_e] = (ζ_bw_zw*ζ_zμ_e - ζ_bμ_e) / (1-ζ_bw_zw)

    # some ratios/elasticities
    Rkstar        = m[:spr]*m[:π_star]*m[:rstar] # (r_k_star+1-δ)/Upsilon*π_star
    ζ_gw       = dGdω_star/Gstar*ωbarstar
    ζ_Gσ_ω    = dGdσstar/Gstar*σ_ω_star

    # elasticities for the net worth evolution
    m[:ζ_nRk]    = m[:γ_star]*Rkstar/m[:π_star]/exp(m[:z_star])*(1+Rhostar)*(1 - μ_estar*Gstar*(1 - ζ_gw/ζ_zw))
    m[:ζ_nR]     = m[:γ_star]/m[:β]*(1+Rhostar)*(1 - nkstar + μ_estar*Gstar*m[:spr]*ζ_gw/ζ_zw)
    m[:ζ_nqk]    = m[:γ_star]*Rkstar/m[:π_star]/exp(m[:z_star])*(1+Rhostar)*(1 - μ_estar*Gstar*(1+ζ_gw/ζ_zw/Rhostar)) - m[:γ_star]/m[:β]*(1+Rhostar)
    m[:ζ_nn]     = m[:γ_star]/m[:β] + m[:γ_star]*Rkstar/m[:π_star]/exp(m[:z_star])*(1+Rhostar)*μ_estar*Gstar*ζ_gw/ζ_zw/Rhostar
    m[:ζ_nμ_e]   = m[:γ_star]*Rkstar/m[:π_star]/exp(m[:z_star])*(1+Rhostar)*μ_estar*Gstar*(1 - ζ_gw*ζ_zμ_e/ζ_zw)
    m[:ζ_nσ_ω]  = m[:γ_star]*Rkstar/m[:π_star]/exp(m[:z_star])*(1+Rhostar)*μ_estar*Gstar*(ζ_Gσ_ω-ζ_gw/ζ_zw*ζ_zσ_ω)

    return m
end

function settings_m990!(m::Model990)
    default_settings!(m)
end

"""
```
init_data_transforms!(m::Model990)
```

This function initializes a dictionary of functions that map series read in in levels to the
appropriate transformed value. At the time that the functions are initialized, data is not
itself in memory. These functions are model-specific because they assume that certain series
are available. The keys of data transforms should match exactly the keys of `m.observables`.
"""
function init_data_transforms!(m::Model990)

    # 1. Output growth, per-capita
    m.data_transforms[:obs_gdp] = function (levels)
        # FROM: Level of nominal GDP (FRED :GDP series)
        # TO:   Quarter-to-quarter percent change of real, per-capita GDP, adjusted for population smoothing

        levels[:temp] = percapita(m, :GDP, levels)
        gdp = 1000 * nominal_to_real(:temp, levels)
        hpadjust(oneqtrpctchange(gdp), levels)
    end

    # 2. Aggregate hours, per-capita
    m.data_transforms[:obs_hours] = function (levels)
        # FROM: Average weekly hours (AWHNONAG) & civilian employment (CE16OV)
        # TO:   log (3 * aggregregate weekly hours / 100), per-capita
        # Note: Not sure why the 3 is there.

        aggregateweeklyhours = levels[:AWHNONAG] .* levels[:CE16OV]
        100*(log(3 * aggregateweeklyhours / 100) - log(levels[:filtered_population]))
    end

    # 3. Real wage growth
    m.data_transforms[:obs_wages] = function (levels)
        # FROM: Nominal compensation per hour (:COMPNFB from FRED)
        # TO: quarter to quarter percent change of real compensation (using GDP deflator)

        oneqtrpctchange(nominal_to_real(:COMPNFB, levels))
    end

    # 4. GDP deflator
    m.data_transforms[:obs_gdpdeflator] = function (levels)
        # FROM: GDP deflator (index)
        # TO:   Approximate quarter-to-quarter percent change of gdp deflator,
        #       i.e.  quarterly gdp deflator inflation

        oneqtrpctchange(levels[:GDPCTPI])
    end

    # 5. Core PCE inflation
    m.data_transforms[:obs_corepce] = function (levels)
        # FROM: Core PCE index
        # INTO: Approximate quarter-to-quarter percent change of Core PCE,
        # i.e. quarterly core pce inflation

        oneqtrpctchange(levels[:PCEPILFE])
    end

    # 6. Nominal short-term interest rate (3 months)
    m.data_transforms[:obs_nominalrate] = function (levels)
        # FROM: Nominal effective federal funds rate (aggregate daily data at a
        #       quarterly frequency at an annual rate)
        # TO:   Nominal effective fed funds rate, at a quarterly rate

        annualtoquarter(levels[:DFF])

    end

    # 7. Consumption growth, per-capita
    m.data_transforms[:obs_consumption] = function (levels)
        # FROM: Nominal consumption
        # TO:   Real consumption, approximate quarter-to-quarter percent change,
        #       per capita, adjusted for population filtering

        levels[:temp] = percapita(m, :PCE, levels)
        cons = 1000 * nominal_to_real(:temp, levels)
        hpadjust(oneqtrpctchange(cons), levels)
    end

    # 8. Investment growth, per-capita
    m.data_transforms[:obs_investment] = function (levels)
        # FROM: Nominal investment
        # INTO: Real investment, approximate quarter-to-quarter percent change,
        #       per capita, adjusted for population filtering

        levels[:temp] = percapita(m, :FPI, levels)
        inv = 10000 * nominal_to_real(:temp, levels)
        hpadjust(oneqtrpctchange(inv), levels)
    end

    # 9. Spread: BAA-10yr TBill
    m.data_transforms[:obs_spread] = function (levels)
        # FROM: Baa corporate bond yield (percent annualized), and 10-year
        #       treasury note yield (percent annualized)
        # TO:   Baa yield - 10T yield spread at a quarterly rate
        # Note: Moody's corporate bond yields on the H15 are based on corporate
        #       bonds with remaining maturities of at least 20 years.

        annualtoquarter(levels[:BAA] - levels[:GS10])
    end

    # 10. Long term inflation expectations
    m.data_transforms[:obs_longinflation] = function (levels)
        # FROM: SPF: 10-Year average yr/yr CPI inflation expectations (annual percent)
        # TO:   FROM, less 0.5
        # Note: We subtract 0.5 because 0.5% inflation corresponds to
        #       the assumed long-term rate of 2 percent inflation, but the
        #       data are measuring expectations of actual inflation.

        annualtoquarter(levels[:ASACX10]  .- 0.5)
    end

    # 11. Long rate (10-year, zero-coupon)
    m.data_transforms[:obs_longrate] = function (levels)
        # FROM: pre-computed long rate at an annual rate
        # TO:   10T yield - 10T term premium at a quarterly rate

        annualtoquarter(levels[:FYCZZA] - levels[:THREEFYTP10])
    end

    # 12. Fernald TFP
    m.data_transforms[:obs_tfp] = function (levels)
        # FROM: Fernald's unadjusted TFP series
        # TO:   De-meaned unadjusted TFP series, adjusted by Fernald's
        #       estimated alpha

        tfp_unadj      = levels[:TFPKQ]
        tfp_unadj_mean = mean(tfp_unadj[!isnan(tfp_unadj)])
        (tfp_unadj - tfp_unadj_mean) ./ (4*(1 - levels[:TFPJQ]))
    end

    # Columns 13 - 13 + n_anticipated_shocks
    for i = 1:n_anticipated_shocks(m)
        # FROM: OIS expectations of $i-period-ahead interest rates at a quarterly rate
        # TO:   Same

        m.data_transforms[symbol("obs_ois$i")] = function (levels)
            levels[:, symbol("ant$i")]
        end
    end
end


text:
include("model.jl")

using Knet: Param, @diff, value, grad

const seq_len   = 10
const hm_data   = 20
const hm_epochs = 10

const lr = .01

model = make()
data = [[randn(1,in_size) for _ in 1:seq_len] for __ in 1:hm_data]

runner() =
for i in 1:hm_epochs
    l = 0
    for seq in data
        input = seq[1:end-1]
        label = seq[2:end]
        result = @diff begin
            response, state, memory = prop(model, input)
            loss(response, label)
        end
        l += value(result)
        for layer in model
            for param in fieldnames(Layer)
                setfield!(layer, param, Param(getfield(layer, param) - grad(result, getfield(layer, param)) * lr))
            end
        end
    end

    @show i, l
end ; runner()

# initial test code.

# in_size = 10
# l_size = 10
#
# seq_len = 10
# hm_data = 20
#
# lstm = Layer(in_size,l_size)
#
# state = zeros(1,l_size)
#
# memory = zeros(1,memory_size)
#
#
# data = [[randn(1,in_size) for _ in 1:seq_len] for __ in 1:hm_data]
#
#
#
# main(model, state, memory, data) =
# begin
#
#     for datapoint in data
#
#         in_data = datapoint[1:end-1]
#         out_data = datapoint[2:end]
#
#
#         g = @diff begin
#
#             outs = []
#             for timestep in in_data
#
#                 out, state, memory = lstm(timestep, state, memory)
#                 push!(outs, out)
#
#             end
#
#         sum(sum([(e1-e2).^2 for (e1,e2) in zip(outs, out_data)]))
#         end
#
#         # for param in params(lstm)
#         #     param -= .01 .* grad(g, param)
#         # end
#
#         for field in fieldnames(Layer)
#             setfield!(model, field, Param(getfield(model, field) - .01 .* grad(g, getfield(model, field))))
#         end
#
#     end
#
# end
#
#
# test(model, state, memory, data) =
#
#     for timestep in data[1]
#
#         out, state, memory = lstm(timestep, state, memory)
#
#         # @show out
#         # @show state
#         @show memory
#         println(" ")
#
#     end
#
#
#
# test(lstm, state, memory, data)
#
#
# println("---")
#
# main(lstm, state, memory, data)
#
# println("---")
#
#
# test(lstm, state, memory, data)


text:
using DelimitedFiles
using StatsBase

function answer1()
    input = "input"
    initial_fish = parse.(Int,split(readline(input),","))
    fishes_now = copy(initial_fish)
    new_to_add = 0
    for day in 1:80
        # @show fishes_now
        new_to_add = 0
        for i in  eachindex(fishes_now)
            if fishes_now[i] == 0
                new_to_add +=1
                fishes_now[i] = 6
            else
                fishes_now[i] -=1
            end
        end
        append!(fishes_now, repeat([8],new_to_add))
    end
    fishes_now
end

@show answer1()

"Do not try to bruteforce part 2..."
function answer2()
    input = "input"
    initial_fish = parse.(Int,split(readline(input),","))
    fishes_now = zeros(Int,9)
    days = 256
    for fish in initial_fish
        fishes_now[fish+1] += 1
    end
    for day in 1:1:days
        fishes_to_add = fishes_now[1]
        for i in 2:1:9
            fishes_now[i-1] = copy(fishes_now[i])
        end
        fishes_now[7] += fishes_to_add
        fishes_now[9] = fishes_to_add
    end
    sum(fishes_now)
end

@show answer2()