text:
"""

    arm64
    ~/miniforge3/bin/python3

    https://betterprogramming.pub/using-python-and-r-with-julia-b7019a3d1420

Final list:
[x] Http2Stream can read buffer with given length

[ ] On error close the stream

[ ] inner process reads are blocking ...

## Check which process opened the port.
    sudo lsof -nP -i4TCP:50051 | grep LISTEN

from grpc.tools import command
    command.build_package_protos(self.distribution.package_dir[''])

    https://pypi.org/project/grpcio-tools/

[x] Working Python gRPC server/client example.

[ ] Investigate
    Julia Proto implementation
    abstract type ProtoRpcChannel end
    abstract type ProtoRpcController end
    abstract type AbstractProtoServiceStub{B} end
"""

using Distributed
using gRPC
using Nghttp2
using OpenSSL
using ProtoBuf
using PyCall
using ResumableFunctions
using Sockets
using Test

# Include protobuf codegen files.
include("proto/proto_jl_out/helloworld.jl")
include("proto/proto_jl_out/routeguide.jl")

# Include certificate helper functions.
include("cert_helpers.jl")

"""
    Installs python gRPC modules.

    Equivalent of running the commands:
    python -m pip install grpcio
    python -m pip install grpcio-tools
"""
function python_install_requirements()
    pip_os = pyimport("pip")
    pip_os.main(["list"])

    #pip_os.main(["uninstall", "grpcio-tools", "-y"])
    #pip_os.main(["uninstall", "grpcio", "-y"])

    pip_os.main(["install", "protobuf==3.18.1"])
    pip_os.main(["install", "grpcio==1.41.0"])
    pip_os.main(["install", "grpcio-tools==1.41.0"])
    return nothing
end

private_key_pem, public_key_pem = create_self_signed_certificate()

# Install gRPC modules
python_install_requirements()

"""
    Configures PyCall, adds a new Julia process to handle Python cals.
"""
function configure_pycall()
    println("Configure python [aths]")
    # Create a worker process, where we run python interpreter.
    addprocs(1)

    # Load python wrappers into the worker processes.
    @everywhere include("py_helpers.jl")
    #@everywhere include("test/py_helpers.jl")
end

# Configure PyCall.
configure_pycall()

"""
    Handler.
    gRPC test server implementation.
"""
module RouteGuideTestHandler
using gRPC
using ResumableFunctions
include("proto/proto_jl_out/routeguide.jl")

const GRPC_SERVER = Ref{gRPCServer}()

function GetFeature(point::routeguide.Point)
    println("[Server]->GetFeature")

    feature = routeguide.Feature()
    feature.name = "from_julia"
    feature.location = point
    return feature
end

@resumable function ListFeatures(rect::routeguide.Rectangle)
    println("[Server]->ListFeatures")

    feature = routeguide.Feature()
    feature.name = "enumerate_from_julia_1"
    @yield feature

    feature = routeguide.Feature()
    feature.name = "enumerate_from_julia_2"
    @yield feature

    feature = routeguide.Feature()
    feature.name = "enumerate_from_julia_3"
    @yield feature

    feature = routeguide.Feature()
    feature.name = "enumerate_from_julia_4"
    @yield feature
end

function RouteEcho(route_note::routeguide.RouteNote)
    println("[Server]->RouteEcho")

    res = routeguide.RouteNote()
    res.message = "from_julia"
    return res
end

@resumable function RouteChat(routes::DeserializeStream{routeguide.RouteNote})
    println("[Server]->RouteChat")

    for route in routes
        println("[Server]::RouteChat receving and sending route")
        @yield route
    end
end

function TerminateServer(empty::routeguide.Empty)
    println("[Server]->TerminateServer")
    GRPC_SERVER.x.is_running = false
    return routeguide.Empty()
end

end # module RouteGuideTestHandler

@resumable function ListRouteNotes()
    route_node = routeguide.RouteNote()
    route_node.message = "Julia Client 1"
    @yield route_node

    route_node = routeguide.RouteNote()
    route_node.message = "Julia Client 2"
    @yield route_node

    route_node = routeguide.RouteNote()
    route_node.message = "Julia Client 3"
    @yield route_node
end

function server_call(socket)
    println("[server_call]::[$(Threads.threadid())]")

    controller = gRPCController()

    server = gRPCServer(
        Dict("routeguide.RouteGuide" => RouteGuideTestHandler.routeguide.RouteGuide(RouteGuideTestHandler))
    )
    RouteGuideTestHandler.GRPC_SERVER.x = server

    accepted_socket = accept(socket)

    nghttp2_server_session = Nghttp2.from_accepted(accepted_socket)

    while server.is_running
        handle_request(nghttp2_server_session, controller, server)
    end

    return nothing
end

function helloworld_client_call()
    controller = gRPCController()

    socket = connect(50200)

    client_session = Nghttp2.open(socket)

    # Create gRPC channel.
    grpc_channel = gRPCChannel(client_session)

    greeterClient = helloworld.GreeterBlockingStub(grpc_channel)

    hello_request = helloworld.HelloRequest()
    hello_request.name = "Hello from Julia"

    hello_reply = helloworld.SayHello(greeterClient, controller, hello_request)
    @show hello_reply
end

function client_call(use_ssl::Bool)
    println("[client_call]: use_ssl:$use_ssl")
    controller = gRPCController()

    local socket::IO

    if use_ssl
        socket = connect(50400)
        ssl_ctx = OpenSSL.SSLContext(OpenSSL.TLSv12ClientMethod())
        #OpenSSL.ssl_set_min_protocol_version(ssl_ctx, OpenSSL.TLS1_2_VERSION)
        result = OpenSSL.ssl_set_options(ssl_ctx, OpenSSL.SSL_OP_NO_SSL_MASK)
        result = OpenSSL.ssl_set_alpn(ssl_ctx, OpenSSL.HTTP2_ALPN)

        socket = SSLStream(ssl_ctx, socket, socket)

        connect(socket)
    else
        socket = connect(50200)
    end

    client_session = Nghttp2.open(socket)

    # Create gRPC channel.
    grpc_channel = gRPCChannel(client_session)

    routeGuide = routeguide.RouteGuideBlockingStub(grpc_channel)

    # RouteChat.
    route_nodes = routeguide.RouteChat(routeGuide, controller, ListRouteNotes())
    received_count::Int = 0

    for route_node in route_nodes
        received_count = received_count + 1
    end

    @test received_count == length(collect(ListRouteNotes()))

    println("-> route notes.")

    # Get feature.
    println("=> client_call.GetFeature")
    in_point = routeguide.Point()
    in_point.latitude = 1
    in_point.longitude = 2

    for n in 1:10
        result = routeguide.GetFeature(routeGuide, controller, in_point)
    end

    # List features.
    println("=> client_call.ListFeatures")
    in_rect = routeguide.Rectangle()

    list_features = routeguide.ListFeatures(routeGuide, controller, in_rect)
    for feature in list_features
        println("feature.name: $(feature.name)")
    end

    # Terminate the server.
    println("=> client_call.TerminateServer")
    _ = routeguide.TerminateServer(routeGuide, controller, routeguide.Empty())

    return nothing
end

function server_call()
    socket = listen(50200)
    server_call(socket)
    return nothing
end

function wait_for_server(port::UInt16)
    connected = false
    while !connected
        connected = try
            connect(port)
            true
        catch e
            println("waiting")
            sleep(1)
            false
        end
    end
end

"""
    Use python client.
"""
function test1()
    socket = listen(50200)

    f2 = @spawnat 2 python_client()

    server_call(socket)

    fetch(f2)

    close(socket)

    return nothing
end

function test2()
    socket = listen(50200)

    # listen(), then pass the socket
    f1 = @spawnat 1 server_call(socket)

    client_call(false)

    fetch(f1)

    close(socket)
    return nothing
end

function test3()
    f2 = @spawnat 2 python_server(private_key_pem, public_key_pem)

    wait_for_server(UInt16(50400))

    client_call(true)

    fetch(f2)

    return nothing
end

function test4()
    f2 = @spawnat 2 python_server(private_key_pem, public_key_pem)

    wait_for_server(UInt16(50200))

    client_call(false)

    fetch(f2)

    return nothing
end

@testset "Python client - Julia server" begin
    test1()
    @test true
end

@testset "Julia server and client" begin
    test2()
    @test true
end

@testset "Secure Python server - Julia client" begin
    test3()
    @test true
end

@testset "Insecure Python server - Julia client" begin
    test4()
    @test true
end

@resumable function enumerate_test_features()
    for i in 1:10
        feature = routeguide.Feature()
        feature.name = "enumerate_from_julia_$i"
        @yield feature
    end
end

@testset "SerializeStream" begin
    serialize_stream = SerializeStream(enumerate(enumerate_test_features()))

    deserialize_stream = DeserializeStream{routeguide.Feature}(serialize_stream)

    for (index, value) in enumerate(deserialize_stream)
        @show index, value
    end
end


text:
using ConcaveHull
using Test

# Triangle
p = [[-1.0,0.0],[1.0,0.0],[0.0,1.0]]
h = concave_hull(p)
@test Set(h.vertices) == Set(p)
@test h.converged
@test area(h) == 0.5*2

# Triangle enclosing point
p = [[-1.0,0.0],[1.0,0.0],[0.0,1.0],[0.0,0.5]]
h = concave_hull(p)
@test h.converged
@test Set(p[1:3]) == Set(h.vertices)

# Square
p = [[-1.0,0.0],[1.0,0.0],[1.0,2.0],[-1.0,2.0]]
h = concave_hull(p)
@test Set(h.vertices) == Set(p)
@test area(h) == 4.0
@test h.converged

# Square with multiple points on one side
p = [[-1.0,0.0], [1.0,0.0], [1.0,0.25], [1.0,0.5], [1.0,0.75],
     [1.0,1.0], [1.0,1.25], [1.0,1.5], [1.0,1.75], [1.0,2.0], [-1.0,2.0]]
h = concave_hull(p)
@test area(h) == 4.0
@test h.converged

# Square enclosing point
p = [[-1.0,0.0],[1.0,0.0],[1.0,2.0],[-1.0,2.0],[0.0,1.0]]
h = concave_hull(p)
@test Set(p[1:4]) == Set(h.vertices)
@test h.converged

# Square enclosing many points
p = [[-1.0,0.0],[1.0,0.0],[1.0,2.0],[-1.0,2.0]]
pr = [[-0.99,0.01] + 1.98*rand(2) for i=1:100]
h = concave_hull(vcat(p,pr))
@test issubset(Set(p), Set(h.vertices))
@test h.converged

# Test non-convergence case
p = [[-1.0, -1.0], [0.0, 0.0], [1.0, 1.0], [2.0, 2.0]]
h = concave_hull(p)
@test h.converged == false

ext: jl


text:
push!(LOAD_PATH,"../src/")

using Documenter, Hymod

pages = [
    "Home" => "index.md",
    "API" => "api.md",
]

makedocs(;
    modules = [Hymod],
    authors = "Kel Markert",
    repo = "https://github.com/KMarkert/Hymod.jl/blob/{commit}{path}#L{line}",
    sitename = "Hymod.jl",
    # format = Documenter.HTML(;
    #     prettyurls = get(ENV, "CI", "false") == "true",
    #     canonical = "https://deltares.github.io/Wflow.jl",
    #     assets = String[],
    # ),
    pages = pages,
)

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


text:
@testset "Systems" begin

    @testset "SPSystem" begin
        fs = equations(katsura(5))

        F = SPSystem(fs)

        @test size(F) == (6, 6)
        @test length(F) == 6
        x = rand(Complex{Float64}, 6)
        system_cache = cache(F, x)
        @test system_cache isa HomotopyContinuation.SystemNullCache
        u = zeros(Complex{Float64}, 6)
        evaluate!(u, F, x, system_cache)
        @test evaluate(F, x, system_cache) ≈ u
        @test evaluate(F, x) ≈ u
        @test F(x) ≈ u

        U = zeros(Complex{Float64}, 6, 6)
        jacobian!(U, F, x, system_cache)
        @test jacobian(F, x, system_cache) ≈ U
        @test jacobian(F, x) ≈ U

        evaluate_and_jacobian!(u, U, F, x, system_cache)
        @test all(evaluate_and_jacobian(F, x, system_cache) .≈ (u, U))

        @polyvar x y a b
        F = SPSystem(
            [x^2 + 2y * x + a * x + b + 3.0, b * x * y^3 + a + 2];
            parameters = [a, b],
        )
        @test size(F) == (2, 2)
        @test length(F) == 2
        x = rand(Complex{Float64}, 2)
        p = rand(Complex{Float64}, 2)
        system_cache = cache(F, x, p)
        @test system_cache isa HomotopyContinuation.SystemNullCache
        u = zeros(Complex{Float64}, 2)
        evaluate!(u, F, x, p, system_cache)
        @test evaluate(F, x, p, system_cache) ≈ u
        @test evaluate(F, x, p) ≈ u
        @test F(x, p) ≈ u

        U = zeros(Complex{Float64}, 2, 2)
        jacobian!(U, F, x, p, system_cache)
        @test jacobian(F, x, p, system_cache) ≈ U
        @test jacobian(F, x, p) ≈ U

        evaluate_and_jacobian!(u, U, F, x, p, system_cache)
        @test all(evaluate_and_jacobian(F, x, p, system_cache) .≈ (u, U))
        @test all(evaluate_and_jacobian(F, x, p) .≈ (u, U))
    end

    @testset "FixedHomotopy" begin
        f = SPSystem(equations(katsura(5)))
        g = SPSystem(equations(cyclic(6)))
        H = StraightLineHomotopy(f, g)

        F = FixedHomotopy(H, 0.23)
        x = rand(Complex{Float64}, 6)
        u = zeros(Complex{Float64}, 6)
        U = zeros(Complex{Float64}, 6, 6)

        system_cache = cache(F, x)
        @test system_cache isa HomotopyContinuation.FixedHomotopyCache
        u = zeros(Complex{Float64}, 6)
        evaluate!(u, F, x, system_cache)
        @test evaluate(F, x, system_cache) ≈ u

        U = zeros(Complex{Float64}, 6, 6)
        jacobian!(U, F, x, system_cache)
        @test jacobian(F, x, system_cache) ≈ U

        evaluate_and_jacobian!(u, U, F, x, system_cache)
        @test all(evaluate_and_jacobian(F, x, system_cache) .≈ (u, U))
    end

    @testset "FPSystem" begin
        F = FPSystem(equations(katsura(5)))
        x = rand(Complex{Float64}, 6)
        u = zeros(Complex{Float64}, 6)
        U = zeros(Complex{Float64}, 6, 6)

        system_cache = cache(F, x)
        @test system_cache isa HomotopyContinuation.FPSystemCache
        u = zeros(Complex{Float64}, 6)
        evaluate!(u, F, x, system_cache)
        @test evaluate(F, x, system_cache) ≈ u

        U = zeros(Complex{Float64}, 6, 6)
        jacobian!(U, F, x, system_cache)
        @test jacobian(F, x, system_cache) ≈ U

        evaluate_and_jacobian!(u, U, F, x, system_cache)
        @test all(evaluate_and_jacobian(F, x, system_cache) .≈ (u, U))
    end

    @testset "TotalDegreeSystem" begin
        @polyvar x y z

        F = SPSystem([x^4 - z^4, y^3 - z^3])
        w = ProjectiveVectors.PVector(rand(ComplexF64, 3))

        G = TotalDegreeSystem([x^4 - z^4, y^3 - z^3])
        @test size(G) == (2, 3)
        @test evaluate(F, w) ≈ evaluate(G, w)
        @test jacobian(F, w) ≈ jacobian(G, w)
        u, U = evaluate_and_jacobian(F, w)
        @test u ≈ evaluate(F, w)
        @test U ≈ jacobian(G, w)

        @polyvar x y z

        F_affine = SPSystem([x^4 - 1, y^3 - 1])
        w = rand(ComplexF64, 3)

        G_affine = TotalDegreeSystem([x^4 - 1, y^3 - 1]; affine = true)
        @test size(G_affine) == (2, 2)
        @test evaluate(F_affine, w) ≈ evaluate(G_affine, w)
        @test jacobian(F_affine, w) ≈ jacobian(G_affine, w)

        u, U = evaluate_and_jacobian(F_affine, w)
        @test u ≈ evaluate(F_affine, w)
        @test U ≈ jacobian(G_affine, w)
    end

    @testset "FixedParameterSystem" begin
        @polyvar x y a b
        f = SPSystem([x^2 + y - a, b * x + y * a^2]; parameters = [a, b])
        F = FixedParameterSystem(f, [2.0, 3.0])

        @test size(F) == (2, 2)
        @test length(F) == 2
        x = rand(Complex{Float64}, 2)
        system_cache = cache(F, x)
        @test system_cache isa HomotopyContinuation.FixedParameterSystemCache
        u = zeros(Complex{Float64}, 2)
        evaluate!(u, F, x, system_cache)
        @test evaluate(F, x, system_cache) ≈ u
        @test evaluate(F, x) ≈ u

        U = zeros(Complex{Float64}, 2, 2)
        jacobian!(U, F, x, system_cache)
        @test jacobian(F, x, system_cache) ≈ U
        @test jacobian(F, x) ≈ U

        evaluate_and_jacobian!(u, U, F, x, system_cache)
        @test all(evaluate_and_jacobian(F, x, system_cache) .≈ (u, U))

        set_parameters!(F, (4.0, 2.3))
        @test F.p == [4.0, 2.3]
    end

    @testset "SquaredUpSystem" begin
        # affine
        @polyvar x y z

        F = SPSystem([x^4 - 1, y^3 - 1, z^2 + x, x + y + z - 1, x^2 + z^2 - 3])
        w = rand(ComplexF64, 3)
        A = randn(ComplexF64, 3, 2)
        S = SquaredUpSystem(F, A)
        system_cache = cache(S, w)
        @test system_cache isa AbstractSystemCache

        @test size(S) == (3, 3)
        @test evaluate(S, w, system_cache) ≈ [LinearAlgebra.I A] * evaluate(F, w) atol = 1e-12
        @test jacobian(S, w, system_cache) ≈ [LinearAlgebra.I A] * jacobian(F, w) atol = 1e-12
        u, U = evaluate_and_jacobian(S, w, system_cache)
        @test u ≈ evaluate(S, w)
        @test U ≈ jacobian(S, w)

        # projective
        @polyvar x y z w

        A = randn(ComplexF64, 3, 2)
        F = [x^4 - 1, y^3 - 1, z^2 + x, x^2 + z^2 - 3, x + y + z - 1]
        G = SPSystem(homogenize([LinearAlgebra.I A] * F, w))
        F = SPSystem(homogenize(F, w))
        v = ProjectiveVectors.PVector(rand(ComplexF64, 4))

        S = SquaredUpSystem(F, A, [4, 3, 2, 2, 1])
        system_cache = cache(S, v)
        @test system_cache isa AbstractSystemCache

        @test size(S) == (3, 4)
        @test evaluate(S, v, system_cache) ≈ evaluate(G, v) atol = 1e-12
        @test jacobian(S, v, system_cache) ≈ jacobian(G, v) atol = 1e-12
        u, U = evaluate_and_jacobian(S, v, system_cache)
        @test u ≈ evaluate(S, v)
        @test U ≈ jacobian(S, v)

        @test HC.compute_numerically_degrees(S) == [4, 3, 2]
    end
end


text:
"""
```
JPT{T} <: AbstractRepModel{T}
```

The `JPT` type defines the structure of JPT, an implementation
of the DSGE model in Justiniano et al. (2010)
"Investment Shocks and Business Cycles".

### 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::OrderedDict{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-readable names to row and
column indices in the matrix representations of of the measurement equation and
equilibrium conditions.

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

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

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

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

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

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

* `pseudo_observables::OrderedDict{Symbol,Int}`: Maps each pseudo-observable to
  a row in the model's pseudo-measurement equation matrices.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"m1002\", cached here for
  filepath computation.

* `subspec::String`: The model subspecification number, indicating that
  some parameters from the original model spec (\"ss10\") 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`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model
  units. 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`
  and `Observable` for further details.

* `pseudo_observable_mappings::OrderedDict{Symbol,PseudoObservable}`: A
  dictionary that stores names and transformations to/from model units. See
  `PseudoObservable` for further details.
"""
mutable struct JPT{T} <: AbstractRepModel{T}
    parameters::ParameterVector{T}                         # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                       # model steady-state values
    keys::OrderedDict{Symbol,Int}                          # human-readable names for all the model
                                                           # parameters and steady-states

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

    spec::String                                           # Model specification number (eg "m990")
    subspec::String                                        # 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

    observable_mappings::OrderedDict{Symbol, Observable}
    pseudo_observable_mappings::OrderedDict{Symbol, PseudoObservable}
end

description(m::JPT) = "Justiniano, Primiceri, and Tambalotti (2010), $(m.subspec)."

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

Arguments:
`m:: JPT`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::JPT)
    # Endogenous states
    endogenous_states = [:y_t, :k_t, :L_t, :Rk_t, :w_t, :π_t, :s_t, :λ_t,       # sticky prices
                         :c_t, :R_t, :u_t, :ϕ_t, :i_t, :kbar_t, :wgap_t,
                         :gdp_t, :z_t, :g_t, :μ_t, :λ_p_t, :λ_w_t,
                         :b_t, :mp_t, :Eπ_t, :Ec_t, :Eλ_t, :Eϕ_t, :ERk_t, :Ei_t,
                         :Ew_t, :gdp_t1, :c_t1, :i_t1, :w_t1, :λ_p_t1, :λ_w_t1,
                                                                                # CHECK BACK IF ANY CAN BE AUGMENTED VARIABLES
                         :y_f_t, :k_f_t, :L_f_t, :Rk_f_t, :w_f_t, :s_f_t,       # flexible prices
                         :λ_f_t, :c_f_t, :R_f_t, :u_f_t, :ϕ_f_t,
                         :i_f_t, :kbar_f_t, :wgap_f_t, :gdp_f_t, :Ec_f_t, :Eλ_f_t,
                         :Eϕ_f_t, :ERk_f_t, :Ei_f_t]

    # Exogenous shocks
    exogenous_shocks = [:mp_sh, :z_sh, :g_sh, :μ_sh, :λ_p_sh, :λ_w_sh, :b_sh]

    # Expectations shocks
    expected_shocks = [:Eπ_sh, :Ec_sh, :Eλ_sh, :Eϕ_sh, :ERk_sh, :Ei_sh, :Ew_sh, # sticky prices
                       :Ec_f_sh, :Eλ_f_sh, :Eϕ_f_sh, :ERk_f_sh, :Ei_f_sh]       # flexible prices

    # Equilibrium conditions
    equilibrium_conditions = [:eq_y, :eq_k, :eq_L, :eq_Rk, :eq_w, :eq_π, :eq_s, :eq_λ,     # sticky prices
                              :eq_c, :eq_R, :eq_u, :eq_ϕ, :eq_i, :eq_kbar, :eq_wgap,
                              :eq_gdp, :eq_z, :eq_g, :eq_μ, :eq_λ_p, :eq_λ_w,
                              :eq_b, :eq_mp, :eq_Eπ, :eq_Ec, :eq_Eλ, :eq_Eϕ,
                              :eq_ERk, :eq_Ei, :eq_Ew, :eq_gdp1, :eq_c1, :eq_i1, :eq_w1,    # CHECK BACK IF ANY CAN BE AUGMENTED VARIABLES
                              :eq_λ_p1, :eq_λ_w1,

                              :eq_y_f, :eq_k_f, :eq_L_f, :eq_Rk_f, :eq_w_f, :eq_s_f,       # flexible prices
                              :eq_λ_f, :eq_c_f, :eq_R_f, :eq_u_f, :eq_ϕ_f,
                              :eq_i_f, :eq_kbar_f, :eq_wgap_f, :eq_gdp_f, :eq_Ec_f, :eq_Eλ_f,
                              :eq_Eϕ_f, :eq_ERk_f, :eq_Ei_f]

    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = []

    # Observables
    observables = keys(m.observable_mappings)

    # Pseudo-observables
    pseudo_observables = keys(m.pseudo_observable_mappings)

    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
    for (i,k) in enumerate(pseudo_observables);          m.pseudo_observables[k]          = i end
end

function JPT(subspec::String="ss1";
                   custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                   testing = false)

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

    # Initialize empty model
    m = JPT{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(), OrderedDict{Symbol,Int}(),

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

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            OrderedDict{Symbol,Observable}(),
            OrderedDict{Symbol,PseudoObservable}())

    # Set settings
    model_settings!(m)
    DSGE.default_test_settings!(m)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Set observable and pseudo-observable transformations
    init_observable_mappings!(m)
    #init_pseudo_observable_mappings!(m)

    # Initialize parameters
    init_parameters!(m)

    init_model_indices!(m)
    init_subspec!(m)
    steadystate!(m)

    return m
end

"""
```
init_parameters!(m::JPT)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::JPT)
    # Calibrated parameters
    m <= parameter(:δ, 0.025, fixed=true,
                   description="δ: The capital depreciation rate.",
                   tex_label = "\\delta" )
    m <= parameter(:g_ss_ratio, 0.22, fixed=true,
                   description="g_ss_ratio: steady state government spending to GDP ratio.",
                   tex_label = "g^{ss}" )

    # Non standard devation parameters
    m <= parameter(:α, 0.25, (0., 1.), (1e-5, 0.999), ModelConstructors.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.21, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(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(:ι_w, 0.15, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.5, 0.15), fixed = false,
                   description="ι_w: The weight attributed to last period's wage in wage indexation. "
                    * "(1-ι_w) is the weight attributed to steady-state wages.",
                   tex_label = "\\iota_w")

    m <= parameter(:γ100, 0.53, (-5.0, 5.0), (-5., 5.), ModelConstructors.Untransformed(), Normal(0.5, 0.025), fixed = false,

                   description="γ: The log of the steady-state growth rate of technology.",
                   tex_label = "100\\gamma")

    m <= parameter(:h, 0.85, (-5.0, 5.0), (-5., 5.), ModelConstructors.Untransformed(), Beta(0.5, 0.1), fixed = false,
                   description="h: habit formation parameter.",
                   tex_label = "h")

    m <= parameter(:λ_p_ss, 0.25,(-1000.,1000.),(-1000., 1000.), ModelConstructors.Untransformed(), Normal(0.15, 0.05),  fixed = false,
                   description="λ_p_ss: The steady state net price markup.",
                   tex_label = "\\lambda_p")

    m <= parameter(:λ_w_ss, 0.12,(-1000.,1000.),(-1000., 1000.), ModelConstructors.Untransformed(), Normal(0.15, 0.05), fixed = false,
                   description = "λ_w_ss: The steady state net wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label = "\\lambda_w")

    m <= parameter(:L_ss, 0.5, (-1000., 1000.), (-1e3, 1e3), ModelConstructors.Untransformed(), Normal(0.0,0.5), fixed=false,
                   description="L_ss: The steady state for log hours.",
                   tex_label = "\\log(L)^{ss}")

    m <= parameter(:π_ss100, 0.75, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Normal(0.5, 0.1), fixed = false,

                   description="π_ss100: The steady-state rate of net inflation multiplied by 100.",
                   tex_label = "100 \\pi^{ss}")

    m <= parameter(:Fβ, 0.1, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Gamma(0.25, 0.1), fixed=false,
                   description = "Fβ: Discount rate transformed.",
                   tex_label = "100(\\beta^{-1} - 1)")

    m <= parameter(:ν, 4., (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), Gamma(2, 0.75), fixed = false,
                   description="ν_l: The inverse Frisch elasticity.",
                   tex_label = "\\nu")

    m <= parameter(:ξ_p, 0.85, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.66, 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 (π_ss).",
                   tex_label = "\\xi_p")

    m <= parameter(:ξ_w, 0.75, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.66, 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 = "\\xi_w")

    m <= parameter(:χ, 5., (0., 10.), (1e-5, 0.), ModelConstructors.Exponential(), Gamma(5., 1.), fixed = false,
                   description="χ: The elasticity of the capital utilization cost function.",
                   tex_label = "\\chi")

    m <= parameter(:S′′, 2.5, (-15., 15.), (-15., 15.), ModelConstructors.Untransformed(), Gamma(4., 1), fixed = false,
                   description="S′′: The investment adjust cost.",
                   tex_label = "S^{\\prime\\prime}")

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

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

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

    m <= parameter(:η_λ_p, 0.75, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.50, 0.20), fixed = false,
                   description="η_λ_p: Moving average component in the price markup shock.",
                   tex_label = "\\eta_{\\lambda_p}")

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

    # exogenous processes - autocorrelation
    m <= parameter(:ρ_R, 0.8, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_R: persistence in the monetary policy rule.",
                   tex_label = "\\rho_{R}")

    m <= parameter(:ρ_z, 0.25, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_{z}")

    m <= parameter(:ρ_g, 0.98, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")

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

    m <= parameter(:ρ_λ_p, 0.95, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.6, 0.2), fixed = false,
                   description="ρ_λ_p: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_p}")

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

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

    m <= parameter(:ρ_mp, 0.15, (0., 1.), (1e-5, 0.999), ModelConstructors.SquareRoot(), Beta(0.4, 0.2), fixed = false,
                   description="ρ_mp: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{mp}")

    # exogenous processes - standard deviation
    m <= parameter(:σ_mp, 0.2, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(2, 0.10), fixed = false,
                   description="σ_mp: The standard deviation of the monetary policy shock.",
                   tex_label = "\\sigma_{mp}")

    m <= parameter(:σ_z, 0.9, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_z: The standard deviation of the technology process.",
                   tex_label = "\\sigma_{z}")

    m <= parameter(:σ_g, 0.35, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.5, 1), fixed = false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label = "\\sigma_{g}")

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

    m <= parameter(:σ_λ_p, 0.15, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.1, 1), fixed = false,
                   description="σ_λ_p: The mean of the process that generates the price elasticity of the composite good. " *
                   "Specifically, the elasticity is (1+λ_{f,t})/(λ_{p_t}).",
                   tex_label = "\\sigma_{\\lambda_p}")

    m <= parameter(:σ_λ_w, 0.2, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.1, 1), fixed = false,
                   tex_label = "\\sigma_{\\lambda_w}")

    m <= parameter(:σ_b, 0.05, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), RootInverseGamma(0.1, 1), fixed = false,
                   description="σ_b: The standard deviation of the intertemporal preference shifter process.",
                   tex_label = "\\sigma_{b}")

    # steady states
    m <= SteadyStateParameter(:γ, NaN, tex_label = "\\gamma")
    m <= SteadyStateParameter(:β, NaN, tex_label = "\\beta")
    m <= SteadyStateParameter(:r_ss, NaN, tex_label = "r^{ss}")
    m <= SteadyStateParameter(:r_ss100, NaN, tex_label = "100 r^{ss}")
    m <= SteadyStateParameter(:π_ss, NaN, tex_label = "\\pi^{ss}")
    m <= SteadyStateParameter(:g_ss, NaN, tex_label = "G^{ss}")
    m <= SteadyStateParameter(:expL_ss, NaN, tex_label = "L^{ss}")
    m <= SteadyStateParameter(:Rk_ss, NaN, description = "Steady-state short-term rate of return on capital.",
                              tex_label = "R^{k, ss}")
    m <= SteadyStateParameter(:s_ss, NaN, description = "Steady-state marginal cost", tex_label = "s^{ss}")
    m <= SteadyStateParameter(:w_ss, NaN, description = "Steady-state wage", tex_label = "w^{ss}")
    m <= SteadyStateParameter(:kL_ss, NaN, tex_label = "k^{ss}/L^{ss}")
    m <= SteadyStateParameter(:FL_ss, NaN, description = "Steady-state ratio of fixed costs to labor", tex_label = "F^{ss}/L^{ss}")
    m <= SteadyStateParameter(:yL_ss, NaN, description = "Steady-state ratio of output to labor", tex_label = "y^{ss}/L^{ss}")
    m <= SteadyStateParameter(:k_ss, NaN, tex_label = "k^{ss}")
    m <= SteadyStateParameter(:i_ss, NaN, tex_label = "i^{ss}")
    m <= SteadyStateParameter(:F_ss, NaN, tex_label = "F^{ss}")
    m <= SteadyStateParameter(:y_ss, NaN, tex_label = "y^{ss}")
    m <= SteadyStateParameter(:c_ss, NaN, tex_label = "c^{ss}")

end

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

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

    m[:γ]        = m[:γ100] / 100.
    m[:β]        = 100. / (m[:Fβ] + 100.)
    m[:r_ss]     = exp(m[:γ]) / m[:β] - 1.
    m[:r_ss100]  = 100. * m[:r_ss]
    m[:π_ss]     = m[:π_ss100] / 100.
    m[:g_ss]     = 1. / (1. - m[:g_ss_ratio])

    m[:expL_ss]  = exp(m[:L_ss])
    m[:Rk_ss]    = exp(m[:γ]) / m[:β] - 1. + m[:δ]
    m[:s_ss]     = 1. / (1. + m[:λ_p_ss])
    m[:w_ss]     = (m[:s_ss] * ((1. - m[:α]) ^ (1. - m[:α])) / ((m[:α] ^ (-m[:α])) * m[:Rk_ss] ^ m[:α])) ^ (1. / (1. - m[:α]))
    m[:kL_ss]    = (m[:w_ss] / m[:Rk_ss]) * (m[:α] / (1. - m[:α]))
    m[:FL_ss]    = m[:kL_ss] ^ m[:α] - m[:Rk_ss] * m[:kL_ss] - m[:w_ss]
    m[:yL_ss]    = m[:kL_ss] ^ m[:α] - m[:FL_ss]
    m[:k_ss]     = m[:kL_ss] * m[:expL_ss]
    m[:i_ss]     = (1. - (1. - m[:δ]) * exp(-m[:γ])) * m[:k_ss] * exp(m[:γ])
    m[:F_ss]     = m[:FL_ss] * m[:expL_ss]
    m[:y_ss]     = m[:yL_ss] * m[:expL_ss]
    m[:c_ss]     = m[:y_ss] / m[:g_ss] -  m[:i_ss]

    return m
end

function model_settings!(m::JPT)

    DSGE.default_settings!(m)

    # Anticipated shocks
    m <= Setting(:n_anticipated_shocks, 0,
                 "Number of anticipated policy shocks")
    m <= Setting(:n_anticipated_shocks_padding, 0,
                 "Padding for anticipated policy shocks")

    # Data
    m <= Setting(:data_id, 3, "Dataset identifier")
    if get_setting(m, :cond_id) in collect(1:5)
        m <= Setting(:cond_full_names, [:obs_gdp, :obs_corepce, :obs_spread, :obs_nominalrate, :obs_longrate],
                     "Observables used in conditional forecasts")
    elseif get_setting(m, :cond_id) == 6
        m <= Setting(:cond_full_names, [:obs_gdp, :obs_corepce, :obs_spread, :obs_nominalrate, :obs_longrate,
                                        :obs_gdpdeflator], "Observables used in conditional forecasts")
    end
    m <= Setting(:cond_semi_names, [:obs_spread, :obs_nominalrate, :obs_longrate],
                 "Observables used in semiconditional forecasts")

    # Forecast
    m <= Setting(:use_population_forecast, true,
                 "Whether to use population forecasts as data")
    m <= Setting(:shockdec_startdate, Nullable(quartertodate("2007-Q1")),
                 "Date of start of shock decomposition output period. If null, then shockdec starts at date_mainsample_start")
end


text:
struct Atoms
  x :: Array{Float64}
  status :: Vector{Int64}
end
Atoms(n :: Int64) = Atoms(zeros(n,2),zeros(Int64,n))


text:
@symbol_func function phi_PF(cur_reactor::AbstractReactor)
  cur_phi = cur_reactor.R_0 ^ 2

  cur_phi -= ( R_CS(cur_reactor) + d(cur_reactor) ) ^ 2

  cur_phi *= cur_reactor.pi

  cur_phi *= B_V(cur_reactor)

  cur_phi
end


text:
export with_options

struct OptionsFactory{T}
	# See https://github.com/JuliaOpt/JuMP.jl/blob/94e2cbb3b1a079b327d2fb9b80bfd3b3afa3415d/src/JuMP.jl#L84-L121
	ObjectType::Type
	args::Tuple
	kwargs::Base.Iterators.Pairs
end

function with_options(object_type, args...; kwargs...)
	OptionsFactory{object_type}(object_type, args, kwargs)
end

function (options_factory::OptionsFactory{T})(args...; kwargs...) where {T}
    return options_factory.ObjectType(args..., options_factory.args...;
									kwargs..., options_factory.kwargs...)
end

"Absctract supertype for accelerator activation."
abstract type AbstractActivationReason end

"Activate accelerator immediately."
struct ImmediateActivation <: AbstractActivationReason end

"""
	COSMO.Settings{T}(; kwargs) where {T <: AbstractFloat}

Creates a COSMO settings object that is used to pass user settings to the solver.

Argument | Description | Values (default)
--- | :--- | :---
rho | ADMM rho step | 0.1
sigma | ADMM sigma step | 1e-6
alpha | Relaxation parameter | 1.6
eps_abs | Absolute residual tolerance | 1e-5
eps_rel | Relative residual tolerance | 1e-5
nearly_ratio | Residual tolerance ratio between `MOI.NEARLY_FEASIBLE_POINT` and ``MOI.FEASIBLE_POINT`` | 100
eps\\_prim\\_inf | Primal infeasibility tolerance | 1e-5
eps\\_dual\\_inf | Dual infeasibility tolerance | 1e-5
max_iter | Maximum number of iterations | 5000
verbose | Verbose printing | false
verbose_timing | Verbose timing | false
kkt_solver | Linear System solver | `QdldlKKTSolver`
check_termination | Check termination interval | 25
check_infeasibility | Check infeasibility interval | 40
scaling | Number of scaling iterations | 10
adaptive_rho | Automatic adaptation of step size parameter | true
adaptive_rho_max_adaptions | Max number of rho adaptions | typemax(Int64) (deactivated)
decompose | Activate to decompose chordal psd constraints | true
complete_dual | Activate to complete the dual variable after decomposition | false
merge_strategy | Choose a strategy for clique merging | `CliqueGraphMerge`
compact_transformation | Choose how a decomposed problem is transformed | true
time_limit | Set solver time limit in s | 0 (deactivated)
accelerator | Acceleration scheme | `AndersonAccelerator{T, Type2{QRDecomp}, RestartedMemory, NoRegularizer}`
accelerator_activation | Accelerator activation | `ImmediateActivation`
safeguard | Accelerator safeguarding | true
safeguard_tol | Safeguarding tolerance | 2.0

"""
mutable struct Settings{T <: AbstractFloat}
	rho::T
	sigma::T
	alpha::T
	eps_abs::T
	eps_rel::T
	nearly_ratio::T
	eps_prim_inf::T
	eps_dual_inf::T
	max_iter::Int
	verbose::Bool
	kkt_solver::Union{Type{<: AbstractKKTSolver}, OptionsFactory{<: AbstractKKTSolver}}
	check_termination::Int
	check_infeasibility::Int
	scaling::Int
	MIN_SCALING::T
	MAX_SCALING::T
	adaptive_rho::Bool
	adaptive_rho_interval::Int
	adaptive_rho_tolerance::T
	adaptive_rho_fraction::T
	adaptive_rho_max_adaptions::Int64
	verbose_timing::Bool
	RHO_MIN::T
	RHO_MAX::T
	RHO_TOL::T
	RHO_EQ_OVER_RHO_INEQ::T
	COSMO_INFTY::T
  	decompose::Bool
  	complete_dual::Bool
	time_limit::T
	obj_true::T
	obj_true_tol::T
	merge_strategy::Union{Type{<: AbstractMergeStrategy}, OptionsFactory{<: AbstractMergeStrategy}}
	compact_transformation::Bool
	accelerator::Union{Type{<: AbstractAccelerator}, OptionsFactory{<: AbstractAccelerator}}
	safeguard::Bool
	safeguard_tol::T

	#constructor
	function Settings{T}(;
		rho::Real=T(0.1),
		sigma::Real=T(1e-6),
		alpha::Real=T(1.6),
		eps_abs::Real=T(1e-5),
		eps_rel::Real=T(1e-5),
		nearly_ratio::Real=T(100),
		eps_prim_inf::Real=T(1e-4),
		eps_dual_inf::Real=T(1e-4),
		max_iter::Integer=5000,
		verbose::Bool=false,
		kkt_solver=QdldlKKTSolver,
		check_termination::Int=25,
		check_infeasibility::Int=40,
		scaling::Integer=10,
		MIN_SCALING::Real = T(1e-4),
		MAX_SCALING::Real = T(1e4),
		adaptive_rho::Bool = true,
		adaptive_rho_interval::Int = 40,
		adaptive_rho_tolerance::Int = 5,
		adaptive_rho_fraction::Real = T(0.4),
		adaptive_rho_max_adaptions::Int = typemax(Int),
		verbose_timing::Bool = false,
		RHO_MIN::Real = T(1e-6),
		RHO_MAX::Real = T(1e6),
		RHO_TOL::Real = T(1e-4),
		RHO_EQ_OVER_RHO_INEQ::Real = T(1e3),
		COSMO_INFTY::Real = T(1e20),
		decompose::Bool = true,
    	complete_dual::Bool = false,
		time_limit::Real = zero(T),
		obj_true::Real = T(NaN),
		obj_true_tol::Real = T(1e-3),
		merge_strategy = CliqueGraphMerge,
		compact_transformation::Bool = true,
		accelerator = with_options(AndersonAccelerator{T, Type2{QRDecomp}, RestartedMemory, NoRegularizer}, mem = 15),
		safeguard::Bool = true,
		safeguard_tol::T = T(2)
		) where {T <: AbstractFloat}
		if !isa(kkt_solver, OptionsFactory)
			kkt_solver = with_options(kkt_solver)
		end

		if !isa(merge_strategy, OptionsFactory)
			merge_strategy = with_options(merge_strategy)
		end

		if !isa(accelerator, OptionsFactory)
			accelerator = with_options(accelerator)
		end

		new(rho, sigma, alpha, eps_abs, eps_rel, nearly_ratio, eps_prim_inf, eps_dual_inf, max_iter, verbose, kkt_solver, check_termination, check_infeasibility, scaling, MIN_SCALING, MAX_SCALING, adaptive_rho, adaptive_rho_interval, adaptive_rho_tolerance, adaptive_rho_fraction, adaptive_rho_max_adaptions, verbose_timing, RHO_MIN, RHO_MAX, RHO_TOL, RHO_EQ_OVER_RHO_INEQ, COSMO_INFTY, decompose, complete_dual, time_limit, obj_true, obj_true_tol, merge_strategy, compact_transformation, accelerator, safeguard, safeguard_tol)
	end
end

# The default case it to return a Settings{Float64} object
Settings(args...; kwargs...) = Settings{DefaultFloat}(args...; kwargs...)

function Base.show(io::IO, obj::COSMO.Settings{T}) where {T <: AbstractFloat}
	println(io,"A COSMO.Settings{$(T)} object. To list the available options type `?` and `help?>COSMO.Settings`.")
end
eltype(::Type{<:Settings{T}}) where {T} = T

function Settings(d::Dict)

	# converts strings to julia types by look-up
	string_settings_converter = Dict("QdldlKKTSolver" => COSMO.QdldlKKTSolver, "CholmodKKTSolver" => COSMO.CholmodKKTSolver, "NoMerge" => COSMO.NoMerge, "ParentChildMerge" => COSMO.ParentChildMerge, "CliqueGraphMerge" => COSMO.CliqueGraphMerge)

	settings = COSMO.Settings{DefaultFloat}()
	# convert pure dictionary settings to COSMO.Settings object
	settings = COSMO.Settings()
	for (key, val) in d
		if key in ["kkt_solver"; "merge_strategy"]
			setfield!(settings, Symbol(key), with_options(string_settings_converter[val]))
		else
			setfield!(settings, Symbol(key), val)
		end
	end
	return settings
end


text: include("smps/smps.jl")


text:
# Welcome to your new Toolips server!

using Main.Toolips

PUBLIC = "../public"

IP = "127.0.0.1"

PORT = 8000

function main()
        # Essentials
    routes = make_routes()
    server_template = ServerTemplate(IP, PORT, routes)
        # Fun stuff (examples !, you should probably delete these.)
    delayed = Route("/delay", fn(delay))
    suicide = Route("/suicide", fn(suicide_fn))
    arguments = Route("/args", fn(args))
    server_template.add(delayed)
    server_template.add(suicide)
    server_template.add(arguments)
    global TLSERVER = server_template.start()
    return(TLSERVER)
end

# Routes

function make_routes()
        # Pages
    four04 = html("
404, Page not found!
")
    index = html("
Hello world!

Not so exciting, is it?
     well, it is a work in progress :p.
")
        # Routes
    routes = []
    homeroute = Route("/", index)
    four04route = Route("404", four04)
    push!(routes, homeroute)
    push!(routes, four04route)
    routes
end

# Routes can either route to a function or a page. Using the html() method,
#   we have avoided making a page. This can be done for anything that is a func.
#   This includes fn, as we can see when these methods are referenced in

suicide_fn = http -> stop!(TLSERVER)

args = http -> string(getargs(http))

function delay(http::Any)

        for character in "Hello World!"

            write(http, string(character))

            sleep(1)

        end

end
main()


text:
@enum AgreementType begin
    require
    ensure
    loopInvariant
end

AgreementsEnabled = true

agreementBreachMessage = Dict(
    require => "Breach on Requirement Expression",
    ensure => "Breach on Ensure Expression",
    loopInvariant => "Breach on Loop Invariant Expression",
)

mutable struct Agreement
    type::AgreementType
    expressions::Array{Union{Expr,Symbol}}
    breachMessage::String
    processFunction::Any
    functionName::Union{String,Nothing}
end

function Agreement(
    agreementType::AgreementType,
    addChecksFunction,
    functionName::Union{String,Nothing},
)::Agreement
    Agreement(
        agreementType,
        [],
        agreementBreachMessage[agreementType],
        addChecksFunction,
        functionName,
    )
end

function (agreement::Agreement)(functionBody::Expr)
    AgreementsEnabled ? agreement.processFunction(functionBody, agreement) : functionBody
end

struct ContractBreachException <: Exception
    functionName::Union{String,Nothing}
    expression::String
    breachMessage::String
end

function Base.showerror(io::IO, e::ContractBreachException)
    print(io, "ContractBreachException: ", e.breachMessage, " '", e.expression, "'")
    if !isnothing(e.functionName)
        print(io, " in function '", e.functionName, "'")
    end
end

function agreementHolds(agreement::Agreement, expressionIndex::Int64)
    contractExpression = agreement.expressions[expressionIndex]
    stringContractExpression = string(contractExpression)
    exceptionThrownExpression = :(throw(
        ContractBreachException(
            $(agreement.functionName),
            $stringContractExpression,
            $(agreement.breachMessage),
        ),
    ))

    checkExpression = Expr(:block)
    checkExpression.head = :if
    checkExpression.args =
        [agreement.expressions[expressionIndex], :nothing, exceptionThrownExpression]
    return checkExpression
end

# Check Expressions list generation

function createCheckExpressions(agreement::Agreement)
    return [agreementHolds(agreement, i) for i = 1:length(agreement.expressions)]
end

function setAgreementEnabling(bool::Bool)
    global AgreementsEnabled
    AgreementsEnabled = bool
end


text:
"""
    f = initmetagauss(sx,st,nso,nto)

Generate the spectral amplitudes `f` for a metagaussian field of size `nso * nso * nto`
with slopes `sx` and `st`.
"""
function initmetagauss(sx,st,ns,nt);

        sx=abs(sx);
        st=abs(st);

#       This spares a lot of memory!
#       [kx,ky,kt]=ndgrid([0:ns/2 -ns/2+1:-1],[0:ns/2 -ns/2+1:-1],[0:nt/2 -nt/2+1:-1]);
        kx=[collect(0:ns/2)' collect(-ns/2+1:-1)']';
        kx=repeat(kx*ones(1,ns),outer=[1,1,nt]);
# kky=permdute(kkx,[2 1 3]);

if(nt>1)
        kt=[collect(0:nt/2)' collect(-nt/2+1:-1)']';
        kt=repeat(kt*ones(1,ns),outer=[1 1 ns]);
        kt=permute(kt,[3 2 1]).^2;

        kt[:,:,1]=0.000001;
end

        kx=kx.^2+permutedims(kx,[2 1 3]).^2;

        kx[1,1,:].=0.000001;
# f = (kx.^(-(sx+1)/4)).*kt.^(-st/4);
    if(nt>1)
        kx = (kx.^(-(sx+1)/4)).*kt.^(-st/4);
    else
        kx = (kx.^(-(sx+1)/4));
    end
    kx[1,1,1]=0;
    kx[1,1,:].=0;
    if(nt>1)
       kx[:,:,1]=0;
    end
   #  kx=kx./sqrt(sum(lin(abs(kx).^2)))*ns*ns*nt;
   kx=kx./sqrt(sum(abs2.(kx[:])))*ns*ns*nt;
end


text:
@testset "Apply mutation" begin

    inpt(size, name="in") = InputSizeVertex(name, size)
    nt(name) = t -> NamedTrait(t, name)
    tf(name) = t -> nt(name)(SizeChangeValidation(t))
    av(outsize, in, name="av") = absorbvertex(identity, outsize, in, traitdecoration=tf(name))

    function mcv(nin, nout, in, name="mcv")
        mm = MatMul(nin, nout)
        return absorbvertex(mm, nout, in, traitdecoration=tf(name)), mm
    end

    function NaiveNASlib.Δnin(v::AbstractVertex, inds::Vector{<:Integer}...)
        function valfun(vv)
            vv ∉ inputs(v) && return ones(nout_org(vv))
            tosel = inds[vv .== inputs(v)]
            value = fill(-10, nout_org(vv))
            selinds = filter(i -> i > 0, inds[inputs(v) .== vv][1])
            value[selinds] .= 1000
            return value
        end

        Δ = max.(0, length.(inds) .- nin(v))
        Δnin(v, Δ...)
        Δoutputs(OutSelectRelaxed(), v, valfun)
    end

    function NaiveNASlib.Δnout(v::AbstractVertex, inds::Vector{<:Integer})
        function valfun(vv)
            vv != v && return ones(nout_org(vv))
            value = fill(-10, nout_org(vv))
            value[filter(i -> i > 0, inds)] .= 1000
            return value
        end

        Δ = max(0, length(inds) - nout(v))
        Δnout(v, Δ)
        Δoutputs(OutSelectRelaxed(), v, valfun)
    end

    @testset "AbsorbVertex mutation" begin
        nin1, nout1 = 3, 5
        nin2, nout2 = 5, 2
        iv1 = av(nin1, inpt(nin1), "iv1")
        mmv1, mm1 = mcv(nin1, nout1, iv1, "mmv1")
        mmv2, mm2 = mcv(nin2, nout2, mmv1, "mmv2")

        # Select subset
        Δnin(mmv2, Integer[1, 3, 4])
        apply_mutation.((mmv1, mmv2))
        @test mm1.W == [ 1 7 10; 2 8 11; 3 9 12]
        @test mm2.W == [ 1 6; 3 8; 4 9]

        # Increase size
        Δnout(mmv1, [collect(1:nout(mmv1))..., -1, -1])
        apply_mutation.((mmv1, mmv2))
        @test mm1.W == [ 1 7 10 0 0; 2 8 11 0 0; 3 9 12 0 0]
        @test mm2.W == [ 1 6; 3 8; 4 9; 0 0; 0 0]
    end

    @testset "StackingVertex mutation" begin

        mmv1, mm1 = mcv(2, 3, inpt(2, "in1"), "mmv1")
        mmv2, mm2 = mcv(3, 4, inpt(3, "in2"), "mmv2")
        cc1 = conc(mmv1, mmv2, dims=2, traitdecoration=tf("cc1"))
        mmv3, mm3 = mcv(7, 3, cc1, "mmv3")

        Δnout(cc1, Integer[1,3,5,7])
        apply_mutation.((cc1, mmv1, mmv2, mmv3))

        @test mm1.W == [1 5; 2 6]
        @test mm2.W == [4 10; 5 11; 6 12]
        @test mm3.W == [1 8 15; 3 10 17; 5 12 19; 7 14 21]

        Δnout(mmv2, [1, -1, -1])
        apply_mutation.((cc1, mmv1, mmv2, mmv3))

        @test mm1.W == [1 5; 2 6]
        @test mm2.W == [4 0 0; 5 0 0; 6 0 0]
        @test mm3.W == [1 8 15; 3 10 17; 5 12 19; 0 0 0; 0 0 0]

        Δnout(mmv1, [2, -1, -1, -1])
        apply_mutation.((cc1, mmv1, mmv2, mmv3))

        @test mm1.W == [5 0 0 0; 6 0 0 0]
        @test mm2.W == [4 0 0; 5 0 0; 6 0 0]
        @test mm3.W == [3 10 17; 0 0 0; ; 0 0 0; 0 0 0; 5 12 19; 0 0 0; 0 0 0]

        Δnin(mmv3, [1, 5])
        apply_mutation.((cc1, mmv1, mmv2, mmv3))

        @test mm1.W == reshape([5; 6], :, 1)
        @test mm2.W == reshape([4; 5; 6], :, 1)
        @test mm3.W == [3 10 17; 5 12 19]
    end

    @testset "InvariantVertex mutation" begin

        addsize = 7
        mmv1, mm1 = mcv(2, addsize, inpt(2, "in1"), "mmv1")
        mmv2, mm2 = mcv(3, addsize, inpt(3, "in2"), "mmv2")
        add = "add" >> mmv1 + mmv2
        mmv3, mm3 = mcv(addsize, 3, add, "mmv3")

        Δnout(add, Integer[1,3,5,7])
        apply_mutation.((add, mmv1, mmv2, mmv3))

        @test mm1.W == [1 5 9 13; 2 6 10 14]
        @test mm2.W == [1 7 13 19; 2 8 14 20; 3 9 15 21]
        @test mm3.W == [1 8 15; 3 10 17; 5 12 19; 7 14 21]

        Δnout(mmv2, [2, 3, -1, -1, -1])
        apply_mutation.((add, mmv1, mmv2, mmv3))

        @test mm1.W == [5 9 0 0 0; 6 10 0 0 0]
        @test mm2.W == [7 13 0 0 0; 8 14 0 0 0; 9 15 0 0 0]
        @test mm3.W == [3 10 17; 5 12 19; 0 0 0; 0 0 0; 0 0 0]
    end

    @testset "No apply to vertex with no inputs" begin
        mmv1, mm1 = mcv(2, 4, inpt(2, "in"), "mmv1")
        mmv2, mm2 = mcv(4, 3, mmv1, "mmv2")

        w1pre = mm1.W
        w2pre = mm2.W

        remove_edge!(mmv1, mmv2)
        apply_mutation.([mmv1, mmv2])

        @test mm1.W == w1pre
        @test mm2.W == w2pre
    end

    @testset "Size only mutation" begin
        mutable struct SizeProbe
            nin
            nout
        end
        NaiveNASlib.mutate_inputs(p::SizeProbe, newsize...) = p.nin = newsize[1]
        NaiveNASlib.mutate_outputs(p::SizeProbe, newsize) = p.nout = newsize
        NaiveNASlib.minΔninfactor(::SizeProbe) = 1
        NaiveNASlib.minΔnoutfactor(::SizeProbe) = 1

        p = SizeProbe(3,5)
        in = av(3, inpt(1, "in"))
        v = absorbvertex(p, p.nout, in, traitdecoration = tf("v"))

        Δnin(v, -1)
        apply_mutation(v)
        @test p.nin == [1,2]
        @test p.nout == [1,2,3,4,5]

        Δnout(v, 2)
        apply_mutation(v)
        @test p.nin == [1, 2]
        @test p.nout == [1,2,3,4,5,-1,-1]
    end

    @testset "Pass kwargs" begin
        mutable struct KwargsProbe
            kws
        end
        KwargsProbe() = KwargsProbe(nothing)
        NaiveNASlib.mutate_inputs(p::KwargsProbe, newsize...; kwargs...) = p.kws = kwargs
        NaiveNASlib.mutate_outputs(p::KwargsProbe, newsize; kwargs...) = p.kws = kwargs

        compgraph_apply(v; kwargs...) = apply_mutation(CompGraph(inputs(v), [v]); kwargs...)

        @testset "Kwargs pass $f" for f in (
            mutate_inputs,
            mutate_outputs,
            apply_mutation,
            compgraph_apply)
            p = KwargsProbe()
            iv = inputvertex("in", 1)
            v = absorbvertex(p, 1, iv)

            f(v; test=:pass)
            @test p.kws == pairs((test=:pass,))
        end
    end
end


text:
import StatsBase: sample

# Model selection

# Taken from https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2017.0009
"""
	$(SIGNATURES)

Computes the Akaike Information Criterion (AIC) given the free parameters `k` for the data `X` and its
estimate `Y` of the model. `likelihood` can be any function of `X` and `Y`.
"""
function AIC(k::Int64, X::AbstractArray, Y::AbstractArray; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return 2*k - 2*log(likelihood(X, Y))
end

# Taken from https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2017.0009
"""
	$(SIGNATURES)

Computes the Akaike Information Criterion compensated for finite samples (AICC) given the free parameters `k` for the data `X` and its
estimate `Y` of the model. `likelihood` can be any function of `X` and `Y`.
"""
function AICC(k::Int64, X::AbstractMatrix, Y::AbstractMatrix; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return AIC(k, X, Y, likelihood = likelihood)+ 2*(k+1)*(k+2)/(size(X)[2]-k-2)
end

function AICC(k::Int64, X::AbstractVector, Y::AbstractVector; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return AIC(k, X, Y, likelihood = likelihood)+ 2*(k+1)*(k+2)/(length(X)-k-2)
end

# Double check on that
# Taken from https://www.immagic.com/eLibrary/ARCHIVES/GENERAL/WIKIPEDI/W120607B.pdf
"""
	$(SIGNATURES)

Computes Bayes Information Criterion (BIC) given the free parameters `k` for the data `X` and its
estimate `Y` of the model. `likelihood` can be any function of `X` and `Y`.
"""
function BIC(k::Int64, X::AbstractMatrix, Y::AbstractMatrix; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return - 2*log(likelihood(X, Y)) + k*log(size(X)[2])
end

function BIC(k::Int64, X::AbstractVector, Y::AbstractVector; likelihood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return - 2*log(likelihood(X, Y)) + k*log(length(X))
end

# Optimal Shrinkage for data in presence of white noise
# See D. L. Donoho and M. Gavish, "The Optimal Hard Threshold for Singular
# Values is 4/sqrt(3)", http://arxiv.org/abs/1305.5870
# Code taken from https://github.com/erichson/optht

function optimal_svht(m::Int64, n::Int64; known_noise::Bool = false)
    @assert m/n > 0
    @assert m/n <= 1

    β = m/n
    ω = (8*β) / (β+1+sqrt(β^2+14β+1))
    c = sqrt(2*(β+1)+ω)

    if known_noise
        return c
    else
        median = median_marcenko_pastur(β)
        return c / sqrt(median)
    end
end

function marcenko_pastur_density(t, lower, upper, beta)
    sqrt((upper-t).*(t-lower))./(2π*beta*t)
end

function incremental_marcenko_pastur(x, beta, gamma)
    @assert beta <= 1
    upper = (1+sqrt(beta))^2
    lower = (1-sqrt(beta))^2

    @inline marcenko_pastur(x) = begin
        if (upper-x)*(x-lower) > 0
            return marcenko_pastur_density(x, lower, upper, beta)
        else
            return zero(eltype(x))
        end
    end

    if gamma ≈ zero(eltype(gamma))
        i, ϵ = quadgk(x->(x^gamma)*marcenko_pastur(x), x, upper)
        return i
    else
        i, ϵ = quadgk(x->marcenko_pastur(x), x, upper)
        return i
    end
end

function median_marcenko_pastur(beta)
    @assert 0 < beta <= 1
    upper = (1+sqrt(beta))^2
    lower = (1-sqrt(beta))^2
    change = true
    x = ones(eltype(upper), 5)
    y = similar(x)
    while change && (upper - lower > 1e-5)
        x = range(lower, upper, length = 5)
        for (i,xi) in enumerate(x)
            y[i] = one(eltype(x)) - incremental_marcenko_pastur(xi, beta, 0)
        end
        any(y .< 0.5) ? lower = maximum(x[y .< 0.5]) : change = false
        any(y .> 0.5) ? upper = minimum(x[y .> 0.5]) : change = false
    end
    return (lower+upper)/2
end

"""
    $(SIGNATURES)

Compute a feature reduced version of the data array `X` via thresholding the
singular values by computing the [optimal threshold for singular values](http://arxiv.org/abs/1305.5870).
"""
function optimal_shrinkage(X::AbstractArray{T, 2}) where T <: Number
    m,n = minimum(size(X)), maximum(size(X))
    U, S, V = svd(X)
    τ = optimal_svht(m,n)
    inds = S .>= τ*median(S)
    return U[:, inds]*Diagonal(S[inds])*V[:, inds]'
end

"""
    $(SIGNATURES)

Compute a feature reduced version of the data array `X` inplace via thresholding the
singular values by computing the [optimal threshold for singular values](http://arxiv.org/abs/1305.5870).
"""
function optimal_shrinkage!(X::AbstractArray{T, 2}) where T <: Number
    m,n = minimum(size(X)), maximum(size(X))
    U, S, V = svd(X)
    τ = optimal_svht(m,n)
    inds = S .>= τ*median(S)
    X .= U[:, inds]*Diagonal(S[inds])*V[:, inds]'
    return
end


text:
function apply(bitmask, rval)
    b = collect(bitstring(parse(Int64, rval))[end - 35:end])
    for i in 1:36
        if bitmask[i] != 'X'
            b[i] = bitmask[i]
        end
    end
    parse(Int64, join(b), base = 2)
end

function solve_decoder_chip_v1(input)
    mem = zeros(Int64, 99999)
    bitmask, tail = Iterators.peel(input)
    bitmask = split(bitmask)[3]
    for line in tail
        lval, _, rval = split(line)
        if lval[4] == '['
            mem[parse(Int64, lval[5:findfirst(x -> x == ']', lval) - 1])] = apply(bitmask, rval)
        else
            bitmask = rval
        end
    end
    sum(mem)
end

function resolve_addresses(bitmask, address)
    b = collect(bitstring(address)[end - 35:end])
    addresses = [join([bitmask[i] == '0' ? b[i] : bitmask[i] for i in 1:36])]
    idxs = findall(x -> x == 'X', addresses[1])
    for i in 1:length(idxs)
        news = []
        for a in addresses
            a0 = collect(a)
            a0[idxs[i]] = '0'
            push!(news, join(a0))
            a1 = collect(a)
            a1[idxs[i]] = '1'
            push!(news, join(a1))
        end
        addresses = [addresses[1 + 2^(i - 1):end]; news]
    end
    parse.(Int64, addresses, base = 2)
end

function solve_decoder_chip_v2(input)
    mem = Dict{Int64, Int64}()
    bitmask, tail = Iterators.peel(input)
    bitmask = split(bitmask)[3]
    for line in tail
        lval, _, rval = split(line)
        if lval[4] == '['
            for a in resolve_addresses(bitmask, parse(Int64, lval[5:findfirst(x -> x == ']', lval) - 1]))
                mem[a] = parse(Int64, rval)
            end
        else
            bitmask = rval
        end
    end
    mem |> values |> sum
end

@assert solve_decoder_chip_v1(eachline("input/test/test_day14.txt")) == 165 ["test case 1 failed"]
@assert solve_decoder_chip_v2(eachline("input/test/test_day14a.txt")) == 208 ["test case 2 failed"]

println("First star: " * string(solve_decoder_chip_v1(readlines("input/day14.txt"))))
println("Second star: " * string(solve_decoder_chip_v2(readlines("input/day14.txt"))))