text:
import JSON
import UUIDs: UUID
import HTTP
import Sockets

import Base: endswith
function endswith(vec::Vector{T}, suffix::Vector{T}) where T
    local liv = lastindex(vec)
    local lis = lastindex(suffix)
    liv >= lis && (view(vec, (liv-lis + 1):liv) == suffix)
end

function set_ENV_defaults()
    for pair in ENV_DEFAULTS
        haskey(ENV, pair.first) || setindex!(ENV, pair.second, pair.first)
    end
end

const connectedclients = Dict{Symbol,Client}()
const notebooks = Dict{UUID,Notebook}()

"Send `messages` to all clients connected to the `notebook`."
function putnotebookupdates!(notebook::Notebook, messages::UpdateMessage...)
    listeners = filter(collect(values(connectedclients))) do c
        c.connected_notebook !== nothing &&
        c.connected_notebook.uuid == notebook.uuid
    end
    for next_to_send in messages, client in listeners
        put!(client.pendingupdates, next_to_send)
    end
    flushallclients(listeners)
    listeners
end

"Send `messages` to all connected clients."
function putplutoupdates!(messages::UpdateMessage...)
    listeners = collect(values(connectedclients))
    for next_to_send in messages, client in listeners
        put!(client.pendingupdates, next_to_send)
    end
    flushallclients(listeners)
    listeners
end

"Send `messages` to a `client`."
function putclientupdates!(client::Client, messages::UpdateMessage...)
    for next_to_send in messages
        put!(client.pendingupdates, next_to_send)
    end
    flushclient(client)
end

"Send `messages` to the `Client` who initiated."
function putclientupdates!(initiator::Initiator, messages::UpdateMessage...)
    putclientupdates!(connectedclients[initiator.clientID], messages...)
end

# https://github.com/JuliaWeb/HTTP.jl/issues/382
const flushtoken = Channel{Nothing}(1)
put!(flushtoken, nothing)

function flushclient(client::Client)
    while isready(client.pendingupdates)
        next_to_send = take!(client.pendingupdates)

        token = take!(flushtoken)
        try
            if client.stream !== nothing
                if isopen(client.stream)
                    client.stream.frame_type = HTTP.WebSockets.WS_BINARY
                    write(client.stream, serialize_message(next_to_send) |> codeunits)
                else
                    @info "Client $(client.id) stream closed."
                    put!(flushtoken, token)
                    return false
                end
            end
        catch ex
            bt = stacktrace(catch_backtrace())
            if ex isa Base.IOError
                # client socket closed, so we return false (about 5 lines after this one)
            else
                @warn "Failed to write to WebSocket of $(client.id) " exception=(ex,bt)
            end
            put!(flushtoken, token)
            return false
        end
        put!(flushtoken, token)
    end
    true
end

function flushallclients(subset::Union{Set{Client}, AbstractArray{Client}})
    disconnected = Set{Symbol}()
    for client in subset
        stillconnected = flushclient(client)
        if !stillconnected
            push!(disconnected, client.id)
        end
    end
    for to_deleteID in disconnected
        delete!(connectedclients, to_deleteID)
    end
end

function flushallclients()
    flushallclients(values(connectedclients))
end

"Will hold all 'response handlers': functions that respond to a WebSocket request from the client. These are defined in `src/webserver/Dynamic.jl`."
responses = Dict{Symbol,Function}()

const MSG_DELIM = "IUUQ.km jt ejggjdvmu vhi" # riddle me this, Julius

"""Start a Pluto server _synchronously_ (i.e. blocking call) on `http://localhost:[port]/`.

This will start the static HTTP server and a WebSocket server. Pluto Notebooks will be started dynamically (by user input)."""
function run(host, port::Integer; launchbrowser::Bool = false)
    set_ENV_defaults()
    hostIP = parse(Sockets.IPAddr, host)
    serversocket = Sockets.listen(hostIP, UInt16(port))
    servertask = @async HTTP.serve(hostIP, UInt16(port), stream = true, server = serversocket) do http::HTTP.Stream
        # messy messy code so that we can use the websocket on the same port as the HTTP server

        if HTTP.WebSockets.is_upgrade(http.message)
            try
                HTTP.WebSockets.upgrade(http) do clientstream
                    if !isopen(clientstream)
                        return
                    end
                    while !eof(clientstream)
                        # This stream contains data received over the WebSocket.
                        # It is formatted and JSON-encoded by send(...) in editor.html
                        try
                            parentbody = let
                                # For some reason, long (>256*512 bytes) WS messages get split up - `readavailable` only gives the first 256*512
                                data = ""
                                while !endswith(data, MSG_DELIM)
                                    if(eof(clientstream))
                                        if data == ""
                                            return
                                        end
                                        @warn "Unexpected eof after" data
                                        data = data * MSG_DELIM
                                        break
                                    end
                                    data = data * String(readavailable(clientstream))
                                end
                                JSON.parse(SubString(data, 1:(lastindex(data)-length(MSG_DELIM))))
                            end
                            process_ws_message(parentbody, clientstream)
                        catch ex
                            if ex isa InterruptException
                                rethrow(ex)
                            elseif ex isa HTTP.WebSockets.WebSocketError
                                # that's fine!
                            elseif ex isa InexactError
                                # that's fine! this is a (fixed) HTTP.jl bug: https://github.com/JuliaWeb/HTTP.jl/issues/471
                                # TODO: remove this switch
                            else
                                bt = stacktrace(catch_backtrace())
                                @warn "Reading WebSocket client stream failed for unknown reason:" exception=(ex,bt)
                            end
                        end
                    end
                end
            catch ex
                if ex isa InterruptException
                    rethrow(ex)
                elseif ex isa Base.IOError
                    # that's fine!
                elseif ex isa ArgumentError && occursin("stream is closed", ex.msg)
                    # that's fine!
                else
                    bt = stacktrace(catch_backtrace())
                    @warn "HTTP upgrade failed for unknown reason" exception=(ex,bt)
                end
            end
        else
            request::HTTP.Request = http.message
            request.body = read(http)
            HTTP.closeread(http)

            request_body = IOBuffer(HTTP.payload(request))
            if eof(request_body)
                # no request body
                response_body = HTTP.handle(PLUTOROUTER, request)
            else
                # there's a body, so pass it on to the handler we dispatch to
                response_body = HTTP.handle(PLUTOROUTER, request, JSON.parse(request_body))
            end

            request.response::HTTP.Response = response_body
            request.response.request = request
            try
                HTTP.startwrite(http)
                write(http, request.response.body)
                HTTP.closewrite(http)
            catch e
                if isa(e, Base.IOError) || isa(e, ArgumentError)
                    # @warn "Attempted to write to a closed stream at $(request.target)"
                else
                    rethrow(e)
                end
            end
        end
    end

    address = if ENV["PLUTO_ROOT_URL"] == "/"
        hostPretty = (hostStr = string(hostIP)) == "127.0.0.1" ? "localhost" : hostStr
        portPretty = Int(port)
        "http://$(hostPretty):$(portPretty)/"
    else
        ENV["PLUTO_ROOT_URL"]
    end
    println("Go to $address to start writing ~ have fun!")
    println()
    println("Press Ctrl+C in this terminal to stop Pluto")
    println()

    launchbrowser && @warn "Not implemented yet"

    # create blocking call:
    try
        wait(servertask)
    catch e
        if isa(e, InterruptException)
            println("\n\nClosing Pluto... Restart Julia for a fresh session. \n\nHave a nice day! 🎈")
            close(serversocket)
            # TODO: HTTP has a kill signal?
            # TODO: put do_work tokens back
            empty!(notebooks)
            for client in values(connectedclients)
                @async close(client.stream)
            end
            empty!(connectedclients)
            for (uuid, ws) in WorkspaceManager.workspaces
                WorkspaceManager.unmake_workspace(ws)
            end
        else
            rethrow(e)
        end
    end
end

run(port::Integer=1234; kwargs...) = run("127.0.0.1", port; kwargs...)

function process_ws_message(parentbody::Dict{String, Any}, clientstream::HTTP.WebSockets.WebSocket)
    client = let
        clientID = Symbol(parentbody["clientID"])
        Client(clientID, clientstream)
    end
    # add to our list of connected clients
    connectedclients[client.id] = client

    messagetype = Symbol(parentbody["type"])
    requestID = Symbol(parentbody["requestID"])

    if messagetype == :disconnect
        delete!(connectedclients, client.id)
        close(clientstream)
    else
        body = parentbody["body"]

        args = []
        if haskey(parentbody, "notebookID")
            notebook = let
                notebookID = UUID(parentbody["notebookID"])
                get(notebooks, notebookID, nothing)
            end
            if notebook === nothing
                messagetype === :connect || @warn "Remote notebook not found locally!"
            else
                client.connected_notebook = notebook
                push!(args, notebook)
            end
        end

        if haskey(parentbody, "cellID")
            cellID = UUID(parentbody["cellID"])
            index = cellindex_fromuuid(notebook, cellID)
            if index === nothing
                @warn "Remote cell not found locally!"
            else
                push!(args, notebook.cells[index])
            end
        end

        if haskey(responses, messagetype)
            responsefunc = responses[messagetype]
            try
                responsefunc(body, args..., initiator=Initiator(client.id, requestID))
            catch ex
                @warn "Response function to message of type $(messagetype) failed"
                rethrow(ex)
            end
        else
            @warn "Message of type $(messagetype) not recognised"
        end
    end
end


text:
include("../logistic-map.jl")

colors = map(RGB24, [0xa50026, 0xd73027, 0xf46d43, 0xfdae61, 0xfee090, 0xe0f3f8, 0xabd9e9, 0x74add1, 0x4575b4, 0x313695])
i=1

layers = [layer(
        y=p[2],
        x=[1:k],
        Geom.point,
        Geom.line,
        Theme(default_color=colors[p[1]])
    ) for p in enumerate(xs)]

gp = plot(layers...,
  Guide.Title("Logistic maps with varying R=[0.1, 4.0]"),
  Guide.YLabel("f(x)"),
  Guide.XLabel("Iteration"))

fpp = plot(x=Rs, y=fps,
  color=[1:length(Rs)],
  Scale.color_discrete_manual(colors...),
  Guide.Title("Reaced fixed point in n'th iteration"),
  Guide.YLabel("Iteration"),
  Guide.XLabel("R"))

draw(PNG("logistic-map-contour.png", 600px, 600px / golden), gp)
draw(PNG("logistic-map-fixedpoints.png", 600px, 600px / golden), fpp)

ext: jl


text:
#Naive convolution method using for loops with some parallel using
#@simd and @inbounds

function convolve!(cLayer::Conv1D, Ai::AbstractArray{T,3}) where {T}

    Aip = padding(cLayer, Ai)

    f_H = cLayer.f
    s_H = cLayer.s
    lastDim = ndims(Ai)
    n_H, c, m = size(cLayer.Z)
    W = cLayer.W
    B = cLayer.B
    @simd for mi = 1:m
        @simd for ci = 1:c
            @simd for hi = 1:n_H
                h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                @inbounds ai = Aip[h_start:h_end, :, mi]
                @inbounds cLayer.Z[hi, ci, mi] = W[:, :, ci] ⋅ ai
                # ai = nothing
                # Base.GC.gc()
            end #for mi=1:m, hi=1:n_H
        end #for ci=1:c
    end #for mi=1:m
    @inbounds cLayer.Z .+= B[:, 1, :]
    Z = cLayer.Z
    actFun = cLayer.actFun
    cLayer.A = eval(:($actFun($Z)))

    W = B = Z = nothing

    Ai = nothing

    # Base.GC.gc()

    return nothing
end #function convolve(cLayer::Conv1D

function convolve!(cLayer::Conv2D, Ai::AbstractArray{T,4}) where {T}

    Aip = padding(cLayer, Ai)

    f_H, f_W = cLayer.f
    s_H, s_W = cLayer.s
    lastDim = ndims(Ai)
    n_H, n_W, c, m = size(cLayer.Z)
    W = cLayer.W
    B = cLayer.B
    @simd for mi = 1:m
        @simd for ci = 1:c
            @simd for wi = 1:n_W
                @simd for hi = 1:n_H
                    h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                    h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                    w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                    w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                    @inbounds ai = Aip[h_start:h_end, w_start:w_end, :, mi]
                    @inbounds cLayer.Z[hi, wi, ci, mi] = W[:, :, :, ci] ⋅ ai
                end #for hi=1:n_H
            end #for mi=1:m, wi=1:n_W, hi=1:n_H
        end #for ci=1:c
    end #for mi=1:m
    @inbounds cLayer.Z .+= B[:, :, 1, :]
    Z = cLayer.Z
    actFun = cLayer.actFun
    cLayer.A = eval(:($actFun($Z)))

    # Base.GC.gc()

    return nothing
end #function convolve(cLayer::Conv2D

function convolve!(cLayer::Conv3D, Ai::AbstractArray{T,5}) where {T}

    Aip = padding(cLayer, Ai)

    f_H, f_W, f_D = cLayer.f
    s_H, s_W, s_D = cLayer.s
    lastDim = ndims(Ai)
    n_H, n_W, n_D, c, m = size(cLayer.Z)
    W = cLayer.W
    B = cLayer.B
    @simd for mi = 1:m
        @simd for ci = 1:c
            @simd for di = 1:n_D
                @simd for wi = 1:n_W
                    @simd for hi = 1:n_H
                        h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                        h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                        w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                        w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                        d_start = di * s_D - (s_D == 1 ? 0 : s_D - 1)
                        d_end = di * s_D - (s_D == 1 ? 0 : s_D - 1) + f_D - 1
                        @inbounds ai = Aip[
                            h_start:h_end,
                            w_start:w_end,
                            d_start:d_end,
                            :,
                            mi,
                        ]
                        @inbounds cLayer.Z[hi, wi, di, ci, mi] =
                            W[:, :, :, :, ci] ⋅ ai
                        # ai = nothing
                        # Base.GC.gc()
                    end #for hi=1:n_H
                end #for wi=1:n_W
            end #for di=1:n_D
        end #for ci=1:c
    end #for mi=1:m
    @inbounds cLayer.Z .+= B[:, :, :, 1, :]
    Z = cLayer.Z
    actFun = cLayer.actFun
    cLayer.A = eval(:($actFun($Z)))

    W = B = Z = nothing

    Ai = nothing

    # Base.GC.gc()

    return nothing
end #function convolve(cLayer::Conv3D

export convolve!

### dconvolve

#Naive convolution method using for loops with some parallel using
#@simd and @inbounds

function dconvolve!(
    cLayer::Conv1D,
    Ai::AbstractArray{T,3},
    dAi::AbstractArray{T,3},
    dZ::AbstractArray{T,3},
) where {T}

    Aip = padding(cLayer, Ai)
    padS = paddingSize(cLayer, Ai)
    dAip = similar(Aip)
    dAip .= 0

    f_H = cLayer.f
    s_H = cLayer.s
    lastDim = ndims(Ai)
    n_H, c, m = size(cLayer.Z)
    W = cLayer.W
    cLayer.dW = similar(W)
    cLayer.dW .= 0
    B = cLayer.B
    cLayer.dB = similar(B)
    cLayer.dB .= 0
    @simd for mi = 1:m
        @simd for ci = 1:c
            for hi = 1:n_H
                h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                @inbounds ai = Aip[h_start:h_end, :, mi]
                @inbounds dAip[h_start:h_end, :, mi] .+=
                    dZ[hi, ci, mi] .* W[:, :, ci]

                @inbounds cLayer.dW[:, :, ci] .+= (ai .* dZ[hi, ci, mi])
                @inbounds cLayer.dB[:, :, ci] .+= dZ[hi, ci, mi]
            end #for
        end #for ci=1:c
    end #for mi=1:m

    dAi .= dAip[1+padS[1]:end-padS[2], :, :]

    return dAi
end #function dconvolve(cLayer::Conv1D

function dconvolve!(
    cLayer::Conv2D,
    Ai::AbstractArray{T,4},
    dAi::AbstractArray{T,4},
    dZ::AbstractArray{T,4},
) where {T}

    Aip = padding(cLayer, Ai)
    padS = paddingSize(cLayer, Ai)
    dAip = similar(Aip)
    dAip .= 0

    f_H, f_W = cLayer.f
    s_H, s_W = cLayer.s
    lastDim = ndims(Ai)
    n_H, n_W, c, m = size(cLayer.Z)
    W = cLayer.W
    cLayer.dW = similar(W)
    cLayer.dW .= 0
    B = cLayer.B
    cLayer.dB = similar(B)
    cLayer.dB .= 0
    @simd for mi = 1:m
        @simd for ci = 1:c
            for wi = 1:n_W, hi = 1:n_H
                h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                @inbounds ai = Aip[h_start:h_end, w_start:w_end, :, mi]
                @inbounds dAip[h_start:h_end, w_start:w_end, :, mi] .+=
                    dZ[hi, wi, ci, mi] .* W[:, :, :, ci]

                @inbounds cLayer.dW[:, :, :, ci] .+= (ai .* dZ[hi, wi, ci, mi])
                @inbounds cLayer.dB[:, :, :, ci] .+= dZ[hi, wi, ci, mi]
            end #for
        end #for ci=1:c
    end #for mi=1:m

    dAi .= dAip[1+padS[1]:end-padS[2], 1+padS[3]:end-padS[4], :, :]

    return dAi
end #function dconvolve(cLayer::Conv2D

function dconvolve!(
    cLayer::Conv3D,
    Ai::AbstractArray{T,5},
    dAi::AbstractArray{T,5},
    dZ::AbstractArray{T,5},
) where {T}

    Aip = padding(cLayer, Ai)
    padS = paddingSize(cLayer, Ai)
    dAip = similar(Aip)
    dAip .= 0

    f_H, f_W, f_D = cLayer.f
    s_H, s_W, s_D = cLayer.s
    lastDim = ndims(Ai)
    n_H, n_W, n_D, c, m = size(cLayer.Z)
    W = cLayer.W
    cLayer.dW = similar(W)
    cLayer.dW .= 0
    B = cLayer.B
    cLayer.dB = similar(B)
    cLayer.dB .= 0
    @simd for mi = 1:m
        @simd for ci = 1:c
            for wi = 1:n_W, hi = 1:n_H, di = 1:n_D
                h_start = hi * s_H - (s_H == 1 ? 0 : s_H - 1)
                h_end = hi * s_H - (s_H == 1 ? 0 : s_H - 1) + f_H - 1
                w_start = wi * s_W - (s_W == 1 ? 0 : s_W - 1)
                w_end = wi * s_W - (s_W == 1 ? 0 : s_W - 1) + f_W - 1
                d_start = di * s_D - (s_D == 1 ? 0 : s_D - 1)
                d_end = di * s_D - (s_D == 1 ? 0 : s_D - 1) + f_D - 1
                @inbounds ai =
                    Aip[h_start:h_end, w_start:w_end, d_start:d_end, :, mi]
                @inbounds dAip[
                    h_start:h_end,
                    w_start:w_end,
                    d_start:d_end,
                    :,
                    mi,
                ] .+= dZ[hi, wi, di, ci, mi] .* W[:, :, :, :, ci]

                @inbounds cLayer.dW[:, :, :, :, ci] .+=
                    (ai .* dZ[hi, wi, ci, mi])
                @inbounds cLayer.dB[:, :, :, :, ci] .+= dZ[hi, wi, di, ci, mi]
            end #for
        end #for ci=1:c
    end #for mi=1:m

    dAi .= dAip[1+padS[1]:end-padS[2],
        1+padS[3]:end-padS[4],
        1+padS[5]:end-padS[6],
        :,
        :,
    ]

    return nothing
end #function dconvolve(cLayer::Conv3D

export dconvolve!


text:
using SemiIsomorphism, ProgressMeter, Permutations, SimpleDrawing
using LinearAlgebraX,
    SimpleGraphs, ChooseOptimizer, Gurobi, SimpleGraphAlgorithms, SimpleTools

use_Gurobi()

"""
    cycle_pair(n::Int)

Find the semi-isomorphism between 2C_n and C_{2n}.
"""
function cycle_pair(n::Int)
    A = adjacency(Cycle(n))
    A = dcat(A, A)
    d = SemiIsomorphism.make_allow(A)
    perms = RPG(2n, d)

    np = length(perms)
    @info "There are $np row permutations to consider"
    PM = Progress(np)

    good_p = Permutation(2n).data
    for p in perms
        B = A[:, p]
        if SemiIsomorphism.is_adjacency(B)
            H = SimpleGraph(B)
            if is_connected(H)
                good_p = p
                break
            end
        end
    end
    if good_p == Permutation(2n).data
        println("No semi-isormorphism found")
        return
    end
    println("Adjacency matrix for 2C_$n")
    display(A)
    println()
    println("Row permutation to make a $(2n) cycle: $(Permutation(good_p))")
    println()
    println("Resulting matrix")
    B = A[:, good_p]
    display(B)
    my_spy(B)
end

function ez_cycle_pair(n::Int)
    A = adjacency(Cycle(n))
    Z = zeros(Int, n, n)
    B = [Z A; A Z]
    G = SimpleGraph(B)
    embed(G, :combined)
    return G
end


text:
@testset "search breakpoints" begin
    Random.seed!(2019)
    p, d = 1000, 20
    H = rand(0.0:1.0, p, d)

    bkpt = p >>> 1 # 500
    X = Vector{Union{Float64, Missing}}(undef, p)
    mask = [rand() < 0.1 for i in 1:p]
    # 1 | 2
    # 1 | 3
    X .= H[:, 1] + [H[1:bkpt, 2]; H[bkpt+1:p, 3]]
    X[mask] .= missing
    bkpt_optim, err_optim = MendelImpute.search_breakpoint(X, H[:, 1], H[:, 2], H[:, 3])
    # @code_warntype MendelImpute.search_breakpoint(X, H[:, 1], H[:, 2], H[:, 3])
    @test bkpt_optim == 499
    @test err_optim == 0

    # 4 | 5
    # 2 | 5
    bkpt = 200
    X .= H[:, 5] + [H[1:bkpt, 4]; H[bkpt+1:p, 2]]
    X[mask] .= missing
    bkpt_optim, err_optim = MendelImpute.search_breakpoint(X, H[:, 5], H[:, 4], H[:, 2])
    @test bkpt_optim == 198
    @test err_optim == 0

    # 2 | 3
    # 5 | 4
    bkpt1 = 200
    bkpt2 = 500
    X .= [H[1:bkpt1, 2]; H[bkpt1+1:p, 5]] + [H[1:bkpt2, 3]; H[bkpt2+1:p, 4]]
    X[mask] .= missing
    bkpt_optim, err_optim = MendelImpute.search_breakpoint(X, H[:, 2], H[:, 5], H[:, 3], H[:, 4])
    @test bkpt_optim == (199, 499)
    @test err_optim == 0

    # 5 | 3
    # 15 | 14
    bkpt1 = 300
    bkpt2 = 700
    X .= [H[1:bkpt1, 5]; H[bkpt1+1:p, 15]] + [H[1:bkpt2, 3]; H[bkpt2+1:p, 14]]
    X[mask] .= missing
    bkpt_optim, err_optim = MendelImpute.search_breakpoint(X, H[:, 5], H[:, 15], H[:, 3], H[:, 14])
    @test bkpt_optim == (297, 699)
    @test err_optim == 0
end


text:
module ODEs

export euler, heun, rk3, rk4,  rk5

include("euler.jl")
include("heun.jl")
include("rk3.jl")
include("rk4.jl")
include("rk5.jl")

end


text:
txt = (readlines("day2.txt") .|> split)
hor = 0
ver = 0

for i in 1:size(txt, 1)
    X = txt[i][2] |> x -> parse(Int64,x)
    if txt[i][1] == "forward"
        hor += X
    else  ver -= (txt[i][1] == "down" ? -X : X)
    end
end

print(hor*ver)


text:
function func_from_include(s)
    length(s) + a
end


text:
function initialize!(solver::iLQRSolver)
	reset!(solver)
    set_verbosity!(solver)
    clear_cache!(solver)

    solver.ρ[1] = solver.opts.bp_reg_initial
    solver.dρ[1] = 0.0

    # Initial rollout
    rollout!(solver)
    TO.cost!(solver.obj, solver.Z)
end

# Generic solve methods
"iLQR solve method (non-allocating)"
function solve!(solver::iLQRSolver{T}) where T<:AbstractFloat
    initialize!(solver)

    Z = solver.Z; Z̄ = solver.Z̄;

    n,m,N = size(solver)
    J = Inf
    _J = TO.get_J(solver.obj)
    J_prev = sum(_J)

    grad_only = false
    for i = 1:solver.opts.iterations
        J = step!(solver, J_prev, grad_only)

        # check for a change in solver status
        status(solver) > SOLVE_SUCCEEDED && break

        copy_trajectories!(solver)

        dJ = abs(J - J_prev)
        J_prev = copy(J)
        gradient_todorov!(solver)

        # Record iteration and evaluate convergence
        record_iteration!(solver, J, dJ)
        exit = evaluate_convergence(solver, i)

        # check if the cost function is quadratic
        if !grad_only && solver.opts.reuse_jacobians && TO.is_quadratic(solver.E)
            @logmsg InnerLoop "Gradient-only LQR"
            grad_only = true  # skip updating feedback gain matrix (feedforward only)
        end

        # check for cost blow up

        # Print iteration
        if is_verbose(solver)
            print_level(InnerLoop, global_logger())
        end
        exit && break

        if J > solver.opts.max_cost_value
            # @warn "Cost exceeded maximum cost"
            solver.stats.status = MAXIMUM_COST
            break
        end
    end
    terminate!(solver)
    return solver
end

# function step!(solver::iLQRSolver, J)
#     TO.state_diff_jacobian!(solver.G, solver.model, solver.Z)
# 	TO.dynamics_expansion!(integration(solver), solver.D, solver.model, solver.Z)
# 	TO.error_expansion!(solver.D, solver.model, solver.G)
#     TO.cost_expansion!(solver.quad_obj, solver.obj, solver.Z, true, true)
#     TO.error_expansion!(solver.E, solver.quad_obj, solver.model, solver.Z, solver.G)
# 	if solver.opts.static_bp
#     	ΔV = static_backwardpass!(solver)
# 	else
# 		ΔV = backwardpass!(solver)
#     end
#     forwardpass!(solver, ΔV, J)
# end

function step!(solver::iLQRSolver{<:Any,<:Any,L}, J, grad_only::Bool=false) where L
    to = solver.stats.to
    init = !solver.opts.reuse_jacobians  # force recalculation if not reusing
    @timeit_debug to "diff jac"     TO.state_diff_jacobian!(solver.G, solver.model, solver.Z)
    if !solver.opts.reuse_jacobians || !(L <: RD.LinearModel) || !grad_only
        @timeit_debug to "dynamics jac" TO.dynamics_expansion!(integration(solver), solver.D, solver.model, solver.Z, solver.cache)
    end
	@timeit_debug to "err jac"      TO.error_expansion!(solver.D, solver.model, solver.G)
    @timeit_debug to "cost exp"     TO.cost_expansion!(solver.quad_obj, solver.obj, solver.Z, solver.exp_cache, init=init, rezero=true)
    @timeit_debug to "cost err"     TO.error_expansion!(solver.E, solver.quad_obj, solver.model, solver.Z, solver.G)
	@timeit_debug to "backward pass" if solver.opts.static_bp
    	ΔV = static_backwardpass!(solver, grad_only)
	else
		ΔV = backwardpass!(solver)
    end
    @timeit_debug to "forward pass" forwardpass!(solver, ΔV, J)
end

"""
$(SIGNATURES)
Simulate the system forward using the optimal feedback gains from the backward pass,
projecting the system on the dynamically feasible subspace. Performs a line search to ensure
adequate progress on the nonlinear problem.
"""
function forwardpass!(solver::iLQRSolver, ΔV, J_prev)
    Z = solver.Z; Z̄ = solver.Z̄
    obj = solver.obj

    _J = TO.get_J(obj)
    J::Float64 = Inf
    α = 1.0
    iter = 0
    z = -1.0
    expected = 0.0
    flag = true

    while (z ≤ solver.opts.line_search_lower_bound || z > solver.opts.line_search_upper_bound) && J >= J_prev

        # Check that maximum number of line search decrements has not occured
        if iter > solver.opts.iterations_linesearch
            for k in eachindex(Z)
                Z̄[k].z = Z[k].z
            end
            TO.cost!(obj, Z̄)
            J = sum(_J)

            z = 0
            α = 0.0
            expected = 0.0

            regularization_update!(solver, :increase)
            solver.ρ[1] += solver.opts.bp_reg_fp
            break
        end

        # Otherwise, rollout a new trajectory for current alpha
        flag = rollout!(solver, α)

        # Check if rollout completed
        if ~flag
            # Reduce step size if rollout returns non-finite values (NaN or Inf)
            # @logmsg InnerIters "Non-finite values in rollout"
            iter += 1
            α /= 2.0
            continue
        end

        # Calcuate cost
        TO.cost!(obj, Z̄)
        J = sum(_J)

        expected::Float64 = -α*(ΔV[1] + α*ΔV[2])
        if expected > 0.0
            z::Float64  = (J_prev - J)/expected
        else
            z = -1.0
        end

        iter += 1
        α /= 2.0
    end

    if J > J_prev
        # error("Error: Cost increased during Forward Pass")
        solver.stats.status = COST_INCREASE
        return NaN
    end

    @logmsg InnerLoop :expected value=expected
    @logmsg InnerLoop :z value=z
    @logmsg InnerLoop :α value=2*α
    @logmsg InnerLoop :ρ value=solver.ρ[1]
    @logmsg InnerLoop :dJ_zero value=solver.opts.dJ_counter_limit

    return J

end

function copy_trajectories!(solver::iLQRSolver)
    for k = 1:solver.N
        solver.Z[k].z = solver.Z̄[k].z
    end
end

"""
Stash iteration statistics
"""
function record_iteration!(solver::iLQRSolver, J, dJ)
    gradient = mean(solver.grad)
    record_iteration!(solver.stats, cost=J, dJ=dJ, gradient=gradient)
    i = solver.stats.iterations::Int

    if dJ ≈ 0
        solver.stats.dJ_zero_counter += 1
    else
        solver.stats.dJ_zero_counter = 0
    end

    @logmsg InnerLoop :iter value=i
    @logmsg InnerLoop :cost value=J
    @logmsg InnerLoop :dJ   value=dJ
    @logmsg InnerLoop :grad value=gradient
    # @logmsg InnerLoop :zero_count value=solver.stats[:dJ_zero_counter][end]
    return nothing
end

"""
$(SIGNATURES)
    Calculate the problem gradient using heuristic from iLQG (Todorov) solver
"""
function gradient_todorov!(solver::iLQRSolver)
	tmp = solver.S[end].r
    for k in eachindex(solver.d)
		tmp .= abs.(solver.d[k])
		u = abs.(control(solver.Z[k])) .+ 1
		tmp ./= u
		solver.grad[k] = maximum(tmp)
    end
end

"""
$(SIGNATURES)
Check convergence conditions for iLQR
"""
function evaluate_convergence(solver::iLQRSolver, ilqr_iterations::Int)
    # Get current iterations
    i = solver.stats.iterations
    grad = solver.stats.gradient[i]
    dJ = solver.stats.dJ[i]

    # Check for cost convergence
    # must satisfy both
    if (0.0 <= dJ < solver.opts.cost_tolerance) && (grad < solver.opts.gradient_tolerance)
        @logmsg InnerLoop "Cost criteria satisfied."
        solver.stats.status = SOLVE_SUCCEEDED
        return true
    end

    # Check total iterations
    if i >= solver.opts.iterations
        @logmsg InnerLoop "Hit max iterations. Terminating."
        solver.stats.status = MAX_ITERATIONS
        return true
    end

    # Exit iLQR and continue solving if the maximum number of iLQR iterations
    # has been reached.
    if ilqr_iterations >= solver.opts.iterations_inner
        return true
    end

    # Outer loop update if forward pass is repeatedly unsuccessful
    if solver.stats.dJ_zero_counter > solver.opts.dJ_counter_limit
        @logmsg InnerLoop "dJ Counter hit max. Terminating."
        solver.stats.status = NO_PROGRESS
        return true
    end

    return false
end

"""
$(SIGNATURES)
Update the regularzation for the iLQR backward pass
"""
function regularization_update!(solver::iLQRSolver,status::Symbol=:increase)
    # println("reg $(status)")
    if status == :increase # increase regularization
        # @logmsg InnerLoop "Regularization Increased"
        solver.dρ[1] = max(solver.dρ[1]*solver.opts.bp_reg_increase_factor, solver.opts.bp_reg_increase_factor)
        solver.ρ[1] = max(solver.ρ[1]*solver.dρ[1], solver.opts.bp_reg_min)
        # if solver.ρ[1] > solver.opts.bp_reg_max
        #     @warn "Max regularization exceeded"
        # end
    elseif status == :decrease # decrease regularization
        # TODO: Avoid divides by storing the decrease factor (divides are 10x slower)
        solver.dρ[1] = min(solver.dρ[1]/solver.opts.bp_reg_increase_factor, 1.0/solver.opts.bp_reg_increase_factor)
        # solver.ρ[1] = solver.ρ[1]*solver.dρ[1]*(solver.ρ[1]*solver.dρ[1]>solver.opts.bp_reg_min)
        solver.ρ[1] = max(solver.opts.bp_reg_min, solver.ρ[1]*solver.dρ[1])
    end
end


text:
using Distributions, Plots, LaTeXStrings, Random; pyplot()

function tStat(mu0,mu,sig,n)
    sample = rand(Normal(mu,sig),n)
    xBar   = mean(sample)
    s      = std(sample)
    (xBar-mu0) / (s/sqrt(n))
end

function powerEstimate(mu0,mu1,sig,n,alpha,N)
    Random.seed!(0)
    sampleH1 = [tStat(mu0,mu1,sig,n) for _ in 1:N]
    critVal  = quantile(TDist(n-1),1-alpha)
    sum(sampleH1 .> critVal)/N
end

mu0 = 20
sig = 5
alpha = 0.05
N = 10^4
rangeMu1 = 16:0.1:30
nList = [5,10,20,30]

powerCurves = [powerEstimate.(mu0,rangeMu1,sig,n,alpha,N) for n in nList]

plot(rangeMu1,powerCurves[1],c=:blue, label="n = $(nList[1])")
plot!(rangeMu1,powerCurves[2],c=:red, label="n = $(nList[2])")
plot!(rangeMu1,powerCurves[3],c=:green, label="n = $(nList[3])")
plot!(rangeMu1,powerCurves[4],c=:purple, label="n = $(nList[4])",
	xlabel= L"\mu", ylabel="Power",
	xlims=(minimum(rangeMu1) ,maximum(rangeMu1)), ylims=(0,1))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
#=
    DenseKKTSystem
=#

"""
    DenseKKTSystem{T, VT, MT} <: AbstractReducedKKTSystem{T, MT}

Implement [`AbstractReducedKKTSystem`](@ref) with dense matrices.

Requires a dense linear solver to be factorized (otherwise an error is returned).

"""
struct DenseKKTSystem{T, VT, MT} <: AbstractReducedKKTSystem{T, MT}
    hess::MT
    jac::MT
    pr_diag::VT
    du_diag::VT
    diag_hess::VT
    # KKT system
    aug_com::MT
    # Info
    ind_ineq::Vector{Int}
    ind_fixed::Vector{Int}
    jacobian_scaling::VT
    # Buffers
    etc::Dict{Symbol, Any}
end

function DenseKKTSystem{T, VT, MT}(n, m, ind_ineq, ind_fixed) where {T, VT, MT}
    ns = length(ind_ineq)
    hess = MT(undef, n, n)
    jac = MT(undef, m, n+ns)
    pr_diag = VT(undef, n+ns)
    du_diag = VT(undef, m)
    diag_hess = VT(undef, n)

    # If the the problem is unconstrained, then KKT system is directly equal
    # to the Hessian (+ some regularization terms)
    aug_com = if (m == 0)
        hess
    else
        MT(undef, n+ns+m, n+ns+m)
    end

    jacobian_scaling = VT(undef, m)

    # Init!
    fill!(aug_com, zero(T))
    fill!(hess,    zero(T))
    fill!(jac,     zero(T))
    fill!(pr_diag, zero(T))
    fill!(du_diag, zero(T))
    fill!(diag_hess, zero(T))
    fill!(jacobian_scaling, one(T))

    return DenseKKTSystem{T, VT, MT}(
        hess, jac, pr_diag, du_diag, diag_hess, aug_com,
        ind_ineq, ind_fixed, jacobian_scaling, Dict{Symbol, Any}(),
    )
end

function DenseKKTSystem{T, VT, MT}(nlp::AbstractNLPModel, info_constraints=get_index_constraints(nlp)) where {T, VT, MT}
    return DenseKKTSystem{T, VT, MT}(
        get_nvar(nlp), get_ncon(nlp), info_constraints.ind_ineq, info_constraints.ind_fixed
    )
end

is_reduced(::DenseKKTSystem) = true
num_variables(kkt::DenseKKTSystem) = length(kkt.pr_diag)

# Special getters for Jacobian
function get_jacobian(kkt::DenseKKTSystem)
    n = size(kkt.hess, 1)
    ns = length(kkt.ind_ineq)
    return view(kkt.jac, :, 1:n)
end

get_raw_jacobian(kkt::DenseKKTSystem) = kkt.jac

nnz_jacobian(kkt::DenseKKTSystem) = length(kkt.jac)

function diag_add!(dest::AbstractMatrix, d1::AbstractVector, d2::AbstractVector)
    n = length(d1)
    @inbounds for i in 1:n
        dest[i, i] = d1[i] + d2[i]
    end
end

function _build_dense_kkt_system!(dest, hess, jac, pr_diag, du_diag, diag_hess, n, m, ns)
    # Transfer Hessian
    for i in 1:n, j in 1:i
        if i == j
            dest[i, i] = pr_diag[i] + diag_hess[i]
        else
            dest[i, j] = hess[i, j]
            dest[j, i] = hess[j, i]
        end
    end
    # Transfer slack diagonal
    for i in 1:ns
        dest[i+n, i+n] = pr_diag[i+n]
    end
    # Transfer Jacobian
    for i in 1:m, j in 1:(n+ns)
        dest[i + n + ns, j] = jac[i, j]
        dest[j, i + n + ns] = jac[i, j]
    end
    # Transfer dual regularization
    for i in 1:m
        dest[i + n + ns, i + n + ns] = du_diag[i]
    end
end

function build_kkt!(kkt::DenseKKTSystem{T, VT, MT}) where {T, VT, MT}
    n = size(kkt.hess, 1)
    m = size(kkt.jac, 1)
    ns = length(kkt.ind_ineq)
    if m == 0 # If problem is unconstrained, just need to update the diagonal
        diag_add!(kkt.aug_com, kkt.diag_hess, kkt.pr_diag)
    else # otherwise, we update the full matrix
        _build_dense_kkt_system!(kkt.aug_com, kkt.hess, kkt.jac, kkt.pr_diag, kkt.du_diag, kkt.diag_hess, n, m, ns)
    end
    treat_fixed_variable!(kkt)
end

function compress_jacobian!(kkt::DenseKKTSystem{T, VT, MT}) where {T, VT, MT}
    m = size(kkt.jac, 1)
    n = size(kkt.hess, 1)
    # Add slack indexes
    for i in kkt.ind_ineq
        kkt.jac[i, i+n] = -one(T)
    end
    # Scale
    kkt.jac .*= kkt.jacobian_scaling
    return
end

function compress_hessian!(kkt::DenseKKTSystem)
    # Transfer diagonal term for future regularization
    diag!(kkt.diag_hess, kkt.hess)
end

function mul!(y::AbstractVector, kkt::DenseKKTSystem, x::AbstractVector)
    mul!(y, kkt.aug_com, x)
end

function jtprod!(y::AbstractVector, kkt::DenseKKTSystem, x::AbstractVector)
    mul!(y, kkt.jac', x)
end

function set_jacobian_scaling!(kkt::DenseKKTSystem, constraint_scaling::AbstractVector)
    copyto!(kkt.jacobian_scaling, constraint_scaling)
end


text:
# file compied from https://github.com/davidanthoff/WRS.jl/tree/master/src/WRS.jl

module WRS

using Distributions
using Distributed

export hd, pb2gen

# This implements (3.16) from (2017) p. 71-72
function hd(X, q=0.5, issorted=false)
    n = length(X)

    a = (n+1)*q
    b = (n+1)*(1-q)

    β = Beta(a, b)

    sorted_X = issorted ? X : sort(X)

    θ = 0.
    for i=1:n
        w = cdf(β, i/n) - cdf(β, (i-1)/n)
        θ += w*sorted_X[i]
    end

    return θ
end

function pb2gen(x, y; quantiles=[0.05, 0.25, 0.5, 0.75, 0.95], alpha=0.05, nboot=2000, estimator=hd)
    bootstrapped_diff_est = pmap(1:nboot) do i
        sampled_values_x = sample(x, length(x), replace=true)
        sampled_values_y = sample(y, length(y), replace=true)

        sort!(sampled_values_x)
        sort!(sampled_values_y)

        return [estimator(sampled_values_x, q, true) - estimator(sampled_values_y, q, true) for q in quantiles]
    end
    bootstrapped_diff_estimates = [[bootstrapped_diff_est[r_i][q_i] for r_i in 1:nboot] for q_i in 1:length(quantiles)]

    return map(enumerate(bootstrapped_diff_estimates)) do (q_i, bootstrapped_diff_est)
        sort!(bootstrapped_diff_est)

        low = round(Int, (alpha/2)*nboot)+1
        up = nboot-low

        q = quantiles[q_i]

        est_x = estimator(x, q)
        est_y = estimator(y, q)
        est_diff = est_x - est_y
        ci_low = bootstrapped_diff_est[low]
        ci_up = bootstrapped_diff_est[up]

        A = count(i->i<0,bootstrapped_diff_est)
        C = count(i->i==0,bootstrapped_diff_est)
        p_hat_star = A/nboot+0.5*C/nboot
        pvalue = 2*(min(p_hat_star,1-p_hat_star))

        se = var(bootstrapped_diff_est)

        return (quantile=q, est_x=est_x, est_y=est_y, est_diff=est_diff, ci_low=ci_low, ci_up=ci_up, pvalue=pvalue, se=se, signif=pvalue < alpha)
    end
end

end # module

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
using Printf, LinearAlgebra,  Plots
theme(:dark)

vhuman=6.0
vraptor0=10.0 # the slow raptor velocity in m/s
vraptor=15.0 #

raptor_distance = 20.0

 # a raptor within 20 cm can attack
raptor_min_distance = 0.2
tmax=10. # the maximum time in seconds
nsteps=1000

"""
This function will compute the derivatives of the
positions of the human and the raptors
"""
function compute_derivatives(angle,h,r0,r1,r2)
    dh = [cos(angle),sin(angle)]*vhuman
    dr0 = (h-r0)/norm(h-r0)*vraptor0
    dr1 = (h-r1)/norm(h-r1)*vraptor
    dr2 = (h-r2)/norm(h-r2)*vraptor
    return dh, dr0, dr1, dr2
end
"""
This function will use forward Euler to simulate the Raptors
"""
function simulate_raptors(angle; output::Bool = true)
    # initial positions
    h = [0.0,0.0]
    r0 = [1.0,0.0]*raptor_distance
    r1 = [-0.5,sqrt(3.)/2.]*raptor_distance
    r2 = [-0.5,-sqrt(3.)/2.]*raptor_distance

    # how much time el
    dt = tmax/nsteps
    t = 0.

    hhist = zeros(2,nsteps+1)
    r0hist = zeros(2,nsteps+1)
    r1hist = zeros(2,nsteps+2)
    r2hist = zeros(2,nsteps+2)

    hhist[:,1] = h
    r0hist[:,1] = r0
    r1hist[:,1] = r1
    r2hist[:,1] = r2

    for i=1:nsteps
        dh, dr0, dr1, dr2 = compute_derivatives(angle,h,r0,r1,r2)
        h += dh*dt
        r0 += dr0*dt
        r1 += dr1*dt
        r2 += dr2*dt
        t += dt

        hhist[:,i+1] = h
        r0hist[:,i+1] = r0
        r1hist[:,i+1] = r1
        r2hist[:,i+1] = r2

        if norm(r0-h) <= raptor_min_distance ||
            norm(r1-h) <= raptor_min_distance ||
            norm(r2-h) <= raptor_min_distance
            if output
                @printf("The raptors caught the human in %f seconds\n", t)
            end

            # truncate the history
            hhist = hhist[:,1:i+1]
            r0hist = r0hist[:,1:i+1]
            r1hist = r1hist[:,1:i+1]
            r2hist = r2hist[:,1:i+1]

            break
        end
    end
    return hhist, r0hist, r1hist, r2hist
end

angle = pi/4.
simulate_raptors(angle)

function show_raptors(angle; args...)
    hhist, r0hist, r1hist, r2hist = simulate_raptors(angle; args...)
    plot(vec(hhist[1,:]),vec(hhist[2,:]),linewidth=3)
    plot!(vec(r0hist[1,:]),vec(r0hist[2,:]),color=:red)
    plot!(vec(r1hist[1,:]),vec(r1hist[2,:]),color=:red)
    plot!(vec(r2hist[1,:]),vec(r2hist[2,:]),color=:red)
    #plot!(grid=false,xticks=false,yticks=false,legend=false,border=false)
    plot!(legend=false, framestyle=:none, size=(600,600))
    plot!(xlim=[-20.,20.],ylim=[-20.,20.])
    annotate!(r0hist[1,1], r0hist[2,1]-1,text("Wounded Raptor",:right,:white))
    annotate!(r1hist[1,1], r1hist[2,1],text("Raptor 1",:right,:white))
    annotate!(r2hist[1,1], r2hist[2,1],text("Raptor 2",:right,:white))
    annotate!(hhist[1,end], hhist[2,end],text("C R U N C H",10,:left,:white))
    title!(@sprintf("Survival time = %f sec",(-1+length(hhist[2,:]))*tmax/nsteps))
end
angle = pi/4
show_raptors(angle)

##
using Interact
ui = button()
settheme!(:nativehtml)
display(ui)

##
mp = @manipulate for θ in 0.0:2π/100:2π
    show_raptors(θ)
    plot!(size=(800,800))
end


text:
using Test
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit.BipartiteGraphs
using ModelingToolkit.StructuralTransformations

function check_contract(sys)
    sys = tearing_substitution(sys)
    state = TearingState(sys)
    fullvars = state.fullvars
    graph = state.structure.graph

    eqs = equations(sys)
    var2idx = Dict(enumerate(fullvars))
    for (i, eq) in enumerate(eqs)
        actual = union(ModelingToolkit.vars(eq.lhs), ModelingToolkit.vars(eq.rhs))
        actual = filter(!ModelingToolkit.isparameter, collect(actual))
        current = Set(fullvars[𝑠neighbors(graph, i)])
        @test isempty(setdiff(actual, current))
    end
end

include("../examples/rc_model.jl")

sys = structural_simplify(rc_model)
check_contract(sys)
@test !isempty(ModelingToolkit.defaults(sys))
u0 = [
      capacitor.v => 0.0
      capacitor.p.i => 0.0
      resistor.v => 0.0
     ]
prob = ODEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Rodas4())

@test sol[resistor.p.i] == sol[capacitor.p.i]
@test sol[resistor.n.i] == -sol[capacitor.p.i]
@test sol[capacitor.n.i] == -sol[capacitor.p.i]
@test iszero(sol[ground.g.i])
@test iszero(sol[ground.g.v])
@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]

# Outer/inner connections
function rc_component(;name)
    R = 1
    C = 1
    @named p = Pin()
    @named n = Pin()
    @named resistor = Resistor(R=R)
    @named capacitor = Capacitor(C=C)
    eqs = [
           connect(p, resistor.p);
           connect(resistor.n, capacitor.p);
           connect(capacitor.n, n);
          ]
    @named sys = ODESystem(eqs, t)
    compose(sys, [p, n, resistor, capacitor]; name=name)
end

@named ground = Ground()
@named source = ConstantVoltage(V=1)
@named rc_comp = rc_component()
eqs = [
       connect(source.p, rc_comp.p)
       connect(source.n, rc_comp.n)
       connect(source.n, ground.g)
      ]
@named sys′ = ODESystem(eqs, t)
@named sys_inner_outer = compose(sys′, [ground, source, rc_comp])
expand_connections(sys_inner_outer, debug=true)
sys_inner_outer = structural_simplify(sys_inner_outer)
@test !isempty(ModelingToolkit.defaults(sys_inner_outer))
u0 = [
      rc_comp.capacitor.v => 0.0
      rc_comp.capacitor.p.i => 0.0
      rc_comp.resistor.v => 0.0
     ]
prob = ODEProblem(sys_inner_outer, u0, (0, 10.0))
sol_inner_outer = solve(prob, Rodas4())
@test sol[capacitor.v] ≈ sol_inner_outer[rc_comp.capacitor.v]

u0 = [
      capacitor.v => 0.0
     ]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

@test sol[resistor.p.i] == sol[capacitor.p.i]
@test sol[resistor.n.i] == -sol[capacitor.p.i]
@test sol[capacitor.n.i] == -sol[capacitor.p.i]
@test iszero(sol[ground.g.i])
@test iszero(sol[ground.g.v])
@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]
#using Plots
#plot(sol)

include("../examples/serial_inductor.jl")
sys = structural_simplify(ll_model)
check_contract(sys)
u0 = [
      inductor1.i => 0.0
      inductor2.i => 0.0
      inductor2.v => 0.0
     ]
prob = ODEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Rodas4())

prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

@variables t x1(t) x2(t) x3(t) x4(t)
D = Differential(t)
@named sys1_inner = ODESystem([D(x1) ~ x1], t)
@named sys1_partial = compose(ODESystem([D(x2) ~ x2], t; name=:foo), sys1_inner)
@named sys1 = extend(ODESystem([D(x3) ~ x3], t; name=:foo), sys1_partial)
@named sys2 = compose(ODESystem([D(x4) ~ x4], t; name=:foo), sys1)
@test_nowarn sys2.sys1.sys1_inner.x1 # test the correct nesting

# compose tests
@parameters t

function record_fun(;name)
    pars = @parameters a=10 b=100
    ODESystem(Equation[], t, [], pars; name)
end

function first_model(;name)
    @named foo=record_fun()

    defs = Dict()
    defs[foo.a] = 3
    defs[foo.b] = 300
    pars = @parameters x=2 y=20
    compose(ODESystem(Equation[], t, [], pars; name, defaults=defs), foo)
end
@named goo = first_model()
@unpack foo = goo
@test ModelingToolkit.defaults(goo)[foo.a] == 3
@test ModelingToolkit.defaults(goo)[foo.b] == 300


text:
# if !isdefined(Main, :Revise) && !haskey(ENV, "JULIA_REVISE")
#     # This doesn't work. Either the ENV from this module is different from the
#     # Revise ENV, or the import happens before ENV can be set.
#     ENV["JULIA_REVISE"] = "manual"   # FIXME; that's still not great
# end
import Revise
using Revise: ModDict, parse_source, RelocatableExpr

# This was exported back when it was stand-alone functionality, but it's an implementation
# detail for TraceCalls
# export apply_code!, revert_code!, update_code_revertible, RevertibleCodeUpdate,
#     CodeUpdate, source
export source

################################################################################

function only(collection)
    @assert length(collection)==1 "only: `collection` was expected to contain one element; contains $(length(collection))"
    return first(collection)
end
only(collection::Base.Generator) = only(collect(collection)) # could be better...

function counter(seq)  # could use DataStructures.counter, but it's a big dependency
    di = Dict()
    for x in seq; di[x] = get(di, x, 0) + 1 end
    di
end

################################################################################

struct CodeUpdate
    md::ModDict
end
CodeUpdate() = CodeUpdate(ModDict())

""" `CodeUpdate(::ModDict)` supports `apply_code!(::CodeUpdate)`, and can be `merge`d
together. """
function Base.merge(cu1::CodeUpdate, cus::CodeUpdate...)
    md = ModDict()
    for cu in [cu1, cus...]
        for (mod::Module, set) in cu.md
            md[mod] = union(get(md, mod, OrderedSet{RelocatableExpr}()), set)
        end
    end
    return CodeUpdate(md)
end
# Base.getindex(cu::CodeUpdate, ind::UnitRange) = CodeUpdate(cu.ecs[ind])
Base.getindex(cu::CodeUpdate, ind::Module) = cu.md[ind]
Base.length(cu::CodeUpdate) = length(cu.md)
apply_code!(cu::CodeUpdate) = scrub_redefinition_warnings() do
    Revise.eval_revised(cu.md)
end
function MakeRelocatableExpr(ex::Expr)
    # because the Revise constructor (via the convert method) is unsafe (!)
    rex = RelocatableExpr(ex.head, ex.args)
    rex.typ = ex.typ
    rex
end

to_expr(rex::RelocatableExpr) =
    # Necessary because sometimes rex.typ is #undef, and Revise.convert(::RelocatableExpr)
    # uses rex.typ. No idea why/when that's happening. Don't think it's on my end.
    Expr(rex.head, rex.args...)

empty_rex = MakeRelocatableExpr(:(identity(nothing))) # a dummy
function apply(fn::Function, rex::RelocatableExpr, mod::Module)
    r = fn(convert(Expr, rex), mod)
    r === nothing ? empty_rex : MakeRelocatableExpr(r)
end

is_empty_rex(rex::RelocatableExpr) = rex === empty_rex
Base.map(fn::Function, cu::CodeUpdate) =
    CodeUpdate(ModDict(mod=>filter(rex->!is_empty_rex(rex),
                                   OrderedSet{RelocatableExpr}(apply(fn, rex, mod)
                                                               for rex in set_rex))
                       for (mod, set_rex) in cu.md))
Base.filter(fn::Function, cu::CodeUpdate) =
    map((expr, mod)->fn(expr, mod) ? expr : nothing, cu) #a lazy & wasteful implementation

""" `RevertibleCodeUpdate(apply::CodeUpdate, revert::CodeUpdate)` contains code
to modify a module, and revert it back to its former state. Use `apply_code!` and
`revert_code!`, or `(::RevertibleCodeUpdate)() do ... end` to temporarily apply the code.
"""
struct RevertibleCodeUpdate
    apply::CodeUpdate
    revert::CodeUpdate
end
function RevertibleCodeUpdate(fn::Function, current::CodeUpdate)
    # Formerly defined as simply
    #    RevertibleCodeUpdate(map((expr, mod)->fn(expr), revert), revert)
    # but we want to get rid of the unmodified code in `set_revert`, so it got imperative.
    md_apply = ModDict()
    md_revert = ModDict()
    for (mod, set_rex) in current.md
        set_apply = OrderedSet{RelocatableExpr}()
        set_revert = OrderedSet{RelocatableExpr}()
        for rex in set_rex
            apply_res = fn(convert(Expr, rex))
            if apply_res !== nothing
                push!(set_apply, MakeRelocatableExpr(apply_res))
                push!(set_revert, rex)
            end
        end
        if !isempty(set_apply)
            md_apply[mod] = set_apply
            md_revert[mod] = set_revert
        end
    end
    return RevertibleCodeUpdate(CodeUpdate(md_apply), CodeUpdate(md_revert))
end

EmptyRevertibleCodeUpdate() = RevertibleCodeUpdate(CodeUpdate(), CodeUpdate())
Base.merge(rcu1::RevertibleCodeUpdate, rcus::RevertibleCodeUpdate...) =
    RevertibleCodeUpdate(merge((rcu.apply for rcu in (rcu1, rcus...))...),
                         merge((rcu.revert for rcu in (rcu1, rcus...))...))
apply_code!(rcu::RevertibleCodeUpdate) = apply_code!(rcu.apply)
revert_code!(rcu::RevertibleCodeUpdate) = apply_code!(rcu.revert)
function (rcu::RevertibleCodeUpdate)(body_fn::Function)
    try
        # It's safer to have the `apply_code!` inside the try, because we should be
        # able to assume that running `revert_code!` is harmless even if apply_code!
        # had an error half-way through.
        apply_code!(rcu)
        @eval $body_fn()   # necessary to @eval because of world age
    finally
        revert_code!(rcu)
    end
end

################################################################################
# These should go into MacroTools/ExprTools

function is_function_definition(expr::Expr)
    l = longdef1(expr)
    l.head == :function && length(l.args) > 1 # `function foo end` is not a definition
end
is_function_definition(::Any) = false
""" `generated2normal(expr::Expr)` makes a normal function out of the staged function.
This, of course, DOES NOT PRESERVE SEMANTICS. It's only useful to parse to `splitdef`
et al. """
function generated2normal(expr::Expr)
    @assert is_generated_function_definition(expr)
    Expr(:function, expr.args...)
end
function normal2generated(expr::Expr)
    @assert expr.head == :function
    Expr(:stagedfunction, expr.args...)
end

is_generated_function_definition(expr::Expr) = expr.head == :stagedfunction
is_generated_function_definition(::Any) = false

is_call_definition(fundef_di::Dict) = @capture(fundef_di[:name], (a_::b_) | (::b_))
is_call_definition(fundef) = is_call_definition(splitdef_memo(fundef))
""" `is_fancy_constructor_definition(fundef)` is true for constructors that have both
parameters and where-parameters (eg. `Vector{T}(x) where T = ...`) """
is_fancy_constructor_definition(fundef_di::Dict) =
    !isempty(get(fundef_di, :params, ())) && !isempty(get(fundef_di, :whereparams, ()))
is_fancy_constructor_definition(fundef) =
    is_fancy_constructor_definition(splitdef_memo(fundef))

const doc_mac = GlobalRef(Core, Symbol("@doc"))

strip_docstring(x) = x
function strip_docstring(x::Expr)
    if x.head == :macrocall && x.args[1] == doc_mac
        strip_docstring(x.args[3])
    else
        x
    end
end

################################################################################

function revertible_update_helper(fn)
    function (code)
        res = fn(code)
        if res === nothing
            (nothing, nothing)
        else
            (res, code)
        end
    end
end

# Return the list of files in `mod`
function module_files(mod::Module)
    if Base.find_source_file(string(mod)) isa Void
        error("TraceCalls error: cannot trace module $mod. It must be defined in a file called $mod.jl whose directory is in LOAD_PATH, and it must have been loaded with `using/import/require` (not `include`).")
    else
        if !haskey(Revise.module2files, Symbol(mod))
            # This will happen, for instance, for modules that were loaded before Revise
            # was loaded.
            Revise.parse_pkg_files(Symbol(mod))
        end
        return Revise.module2files[Symbol(mod)]
    end
end

is_macro_call(ex) = false
is_macro_call(ex::Expr) = ex.head==:macrocall

# Helper function, from Compat. Use Base.macroexpand in 0.7
macroexpandmodule(mod::Module, x::ANY) = eval(mod, :(macroexpand($(QuoteNode(x)))))

#dbg = []
function expand_macros(mod::Module, os::OrderedSet{RelocatableExpr})
    new_set = OrderedSet{RelocatableExpr}()
    trav(::Any) = nothing
    trav(rex::RelocatableExpr) = trav(convert(Expr, rex))
    function trav(expr::Expr)
        if Revise.isdocexpr(expr)
            trav(expr.args[Revise.nargs_docexpr])
        elseif is_macro_call(expr)
            #push!(dbg, (expr => macroexpandmodule(mod, expr)))
            # I wish I could just macroexpand the top-level. #21662
            trav(macroexpandmodule(mod, expr))
        elseif expr.head in (:begin, :block)
            foreach(trav, expr.args)
        else
            push!(new_set, MakeRelocatableExpr(expr))
        end
    end
    foreach(trav, os)
    new_set
end
expand_macros(md::ModDict) =
    Dict(mod=>expand_macros(mod, rexes) for (mod, rexes) in md)

code_of(mod::Module) =
    merge((code_of(mod, file) for file in module_files(mod))...)

function code_of(mod::Module, file::String)
    if haskey(Revise.file2modules, file)
        CodeUpdate(expand_macros(Revise.file2modules[file].md))
    elseif mod==Main
        # That's a bit of a hack, to support interactively-included files. I think it's
        # OK, but not 100% sure. See also Revise#40
        Revise.parse_source(file, Main, dirname(file))
        return (haskey(Revise.file2modules, file) ? code_of(mod, file) : CodeUpdate())
    else
        CodeUpdate()
    end
end
function code_of(included_file::String)::CodeUpdate
    parse_source(included_file, Main, pwd())
    code_of(Main, included_file)
end

method_file_counts(fn_to_change) =
    counter((mod, isa(file, AbstractString) ? realpath(abspath(file)): file)
            # The Set is so that we count methods that have the same file and line number.
            # (i.e. optional files, although it might catch macroexpansions too; not
            # sure if that's good or not)
            for (mod, file, line) in Set((m.module, functionloc(m)...)
                                         for m in methods(fn_to_change).ms))

struct UpdateInteractiveFailure
    fn::Union{Function, Type}
end
Base.show(io::IO, upd::UpdateInteractiveFailure) =
    write(io, "Cannot find source of some method of $(upd.fn). Perhaps it was defined interactively? Try adding `@traceable` in front of its definition.")

struct MissingMethodFailure
    count::Int
    correct_count::Int
    fn::Union{Function, Type}
    file::String
end
Base.show(io::IO, fail::MissingMethodFailure) =
    write(io, "Only $(fail.count)/$(fail.correct_count) methods of $(fail.fn) in $(fail.file) were found.")

""" `parse_mod!` fills up Revise.file2modules for that module, and returns `nothing` """
function parse_mod!(mod::Module)
    if module_parent(mod) !== Main
        parse_mod!(module_parent(mod))
    elseif mod == Base
        if !haskey(Revise.file2modules, Revise.sysimg_path)
            parse_source(Revise.sysimg_path, Main, dirname(Revise.sysimg_path))
        end
    elseif !haskey(Revise.file2modules, module_definition_file(mod))
        Revise.parse_pkg_files(Symbol(mod)) # it's a side-effect of this function...
        @assert haskey(Revise.file2modules, module_definition_file(mod))
    end
    nothing
end

""" `get_function(mod::Module, fundef::Expr)::Function` returns the `Function` which this
`fundef` is defining. This code works only when the Function already exists. """
function get_function_(mod::Module, fundef::Expr)::Union{Function, Type, Void}
    try
        eval(mod, splitdef_memo(fundef)[:name])
    catch
        # This can happen, for instance, in macro expansions that create functions
        # with `gensym` names.
        nothing
    end
end

""" `get_function(mod, def)` returns the ::Function corresponding to `def` if there is
one (and it's not a callable-object definition), otherwise returns `nothing` """
function get_function(mod::Module, def)
    expr = strip_docstring(def)
    if is_function_definition(expr) &&
        !is_call_definition(expr)
        return get_function_(mod, expr)
    elseif is_generated_function_definition(expr)
        return get_function_(mod, generated2normal(expr))
    else
        return nothing
    end
end

function code_of(fn::Union{Function, Type}; when_missing=warn)::CodeUpdate
    if when_missing in (false, nothing); when_missing = _->nothing end
    function process(mod, file, correct_count)
        if file === nothing
            when_missing(UpdateInteractiveFailure(fn))
            return CodeUpdate()
        end
        if mod !== Main && !haskey(Revise.file2modules, file)
            parse_mod!(mod)
            if !haskey(Revise.file2modules, file)
                # Should fail somehow?
                return CodeUpdate()
            end
        end
        count = 0  # how many methods were found
        rcu = filter((expr, expr_mod) -> (get_function(expr_mod, expr) == fn ?
                                          (count += 1; true) : false),
                     code_of(mod, file)::CodeUpdate)
        if count != correct_count
            when_missing(MissingMethodFailure(count, correct_count, fn, file))
        end
        rcu
    end
    return merge((process(mod, file, correct_count)
                  for ((mod, file), correct_count) in method_file_counts(fn))...)
end

"""
    update_code_revertible(new_code_fn::Function, obj::Union{Module, Function, String})

applies the source code transformation function `new_code_fn` to each expression in the
source code of `obj`, and returns a `RevertibleCodeUpdate` which can put into
effect/revert that new code. `obj` can be a module, a function (will transform each
method), or a Main-included ".jl" filename.

`update_code_revertible` itself is side-effect free (it neither modifies the source file,
nor the state of `Julia`). See the README for usage info.

IMPORTANT: if some expression `x` should not be modified, return `nothing` instead of `x`.
This will significantly improve performance. """
function update_code_revertible(fn::Function, mod::Module)
    if mod == Base; error("Cannot trace all of Base. Please trace specific functions (eg. Base.sparse), and beware that tracing ubiquitous, basic functions like `convert` or `length` will trap you in an infinite recursion.") end
    if mod == TraceCalls; error("Cannot trace all of TraceCalls - it yields a StackOverflow, but you may trace specific functions like TraceCalls.redgreen") end
    return RevertibleCodeUpdate(fn, code_of(mod))
end

update_code_revertible(new_code_fn::Function, file::String) =
    RevertibleCodeUpdate(new_code_fn, code_of(file))
update_code_revertible(new_code_fn::Function, fn_to_change::Union{Function, Type};
                       when_missing=warn) =
    RevertibleCodeUpdate(new_code_fn, code_of(fn_to_change; when_missing=when_missing))

""" `source(fn::Function, when_missing=warn)::Vector` returns a vector of the parsed
code corresponding to each method of `fn`. It can fail for any number of reasons,
and `when_missing` is a function that will be passed an exception when it cannot find
the code.
"""
function source(obj::Union{Module, Function}; kwargs...)
    res = []
    # Loop because flattening comprehensions are broken in 0.6
    for (mod2, rex_set) in code_of(obj).md
        for rex in rex_set
            push!(res, to_expr(rex))
        end
    end
    res
end


text:
# const Union{}
# function convert(type::Tuple{Points, Array{Cells, 1}})
# end

# export convert

# @everywhere begin
    function convert(::Type{Lar.Cell}, cell::Lar.Cells)
        return Lar.characteristicMatrix(cell)
    end

    function convert(::Type{Lar.ChainOp}, cells::Lar.Cells)
        return Lar.lar2cop(cells)
    end

    function convert(::Type{Lar.Cells}, chainop::Lar.ChainOp)
        return Lar.cop2lar(chainop)
    end
# end


text:
"""
    AbstractIndexedDiGraph{T} <: AbstractIndexedGraph{T}

Abstract type for representing directed graphs.
"""
abstract type AbstractIndexedDiGraph{T} <: AbstractIndexedGraph{T}; end

edges(g::AbstractIndexedDiGraph) = @inbounds (IndexedEdge{Int}(i, g.A.rowval[k], k) for i=1:size(g.A,2) for k=nzrange(g.A,i))
ne(g::AbstractIndexedDiGraph) = nnz(g.A)

is_directed(g::AbstractIndexedDiGraph) = true
is_directed(::Type{AbstractIndexedDiGraph{T}}) where T = true

"""
    outedges(g::AbstractIndexedDiGraph, i::Integer)

Return a lazy iterator to the edges outgoing from node `i` in `g`.
"""
outedges(g::AbstractIndexedDiGraph, i::Integer) = @inbounds (IndexedEdge{Int}(i, g.A.rowval[k], k) for k in nzrange(g.A, i))

edge_idx(g::AbstractIndexedDiGraph, src::Integer, dst::Integer) = nzindex(g.A, dst, src)
edge_src_dst(g::AbstractIndexedDiGraph, id::Integer) = reverse(nzindex(g.A, id))

"""
    get_edge(g::AbstractIndexedDiGraph, src::Integer, dst::Integer)
    get_edge(g::AbstractIndexedDiGraph, id::Integer)

Get edge given source and destination or given edge index.
"""
function get_edge(g::AbstractIndexedDiGraph, src::Integer, dst::Integer)
    id = edge_idx(g, src, dst)
    IndexedEdge(src, dst, id)
end
function get_edge(g::AbstractIndexedDiGraph, id::Integer)
    i, j = edge_src_dst(g, id)
    IndexedEdge(i, j, id)
end

# Returns sparse adj matrix. Elements default to Int (to match Graphs)
function adjacency_matrix(g::AbstractIndexedDiGraph, T::DataType=Int)
    M = sparse(transpose(g.A))
    SparseMatrixCSC(M.m, M.n, M.colptr, M.rowval, ones(T, nnz(M)))
end

"""
    IndexedDiGraph{T<:Integer} <: AbstractIndexedDiGraph{T}

A type representing a sparse directed graph with access only to outedges.

### FIELDS

- `A` -- square matrix filled with `NullNumber`s. `A[i,j]` corresponds to an edge `j=>i`
"""
struct IndexedDiGraph{T<:Integer} <: AbstractIndexedDiGraph{T}
    A :: SparseMatrixCSC{NullNumber, T}
    function IndexedDiGraph(A::SparseMatrixCSC{NullNumber, T}) where {T<:Integer}
        _checksquare(A)
        _check_selfloops(A)
        new{T}(A)
    end
end

"""
    IndexedDiGraph(A::AbstractMatrix)

Constructs a IndexedDiGraph from the adjacency matrix A.

`IndexedDiGraph` internally stores the transpose of `A`. To avoid overhead due
to the transposition, use `IndexedDiGraph(transpose(At))` where `At` is the
transpose of `A`.
"""
function IndexedDiGraph(A::AbstractMatrix)
    B = convert(SparseMatrixCSC, transpose(A))
    At = SparseMatrixCSC(B.m, B.n, B.colptr, B.rowval, fill(NullNumber(), nnz(B)))
    IndexedDiGraph(At)
end

"""
    IndexedDiGraph(A::AbstractSimpleGraph)

Construct an `IndexedDiGraph` from any `AbstractSimpleGraph` (Graphs.jl),
directed or otherwise.
"""
IndexedDiGraph(sg::AbstractSimpleGraph) = IndexedDiGraph(adjacency_matrix(sg))

# WARNING: very slow! Not recommended, if you need inneighbors, check out `IndexedBiDiGraph`
# here only to comply with the requirements for subtyping `AbstractGraph` from Graphs.jl
function inneighbors(g::IndexedDiGraph, i::Integer)
    X = sparse(transpose(g.A))
    @view X.rowval[nzrange(X,i)]
end

Base.zero(g::IndexedDiGraph) = IndexedDiGraph(zero(g.A))


text:
# Hash
# ====
#
# MurmurHash3 function of BioSequence.
#
# The hash function defined here cares about the starting position of a
# character sequence in the underlying data. That means, even if the starting
# positions of two sequences (`s1` and `s2`) are different in their `data`
# field, their hash values are identical if `s1 == s1` is true.  The
# implementation is based on the 128bit MurmurHash3 function, which was written
# by Austin Appleby, and the source code is distributed under the public domain:
# https://github.com/aappleby/smhasher/blob/master/src/MurmurHash3.cpp
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

@inline function rotl64(x::UInt64, r)
    return (x << r) | (x >> (64 - r))
end

@inline function fmix64(k::UInt64)
    k = k ⊻ k >> 33
    k *= 0xff51afd7ed558ccd
    k = k ⊻ k >> 33
    k *= 0xc4ceb9fe1a85ec53
    k = k ⊻ k >> 33
    return k
end

macro murmur1()
    esc(quote
        k1 *= c1
        k1 = rotl64(k1, 31)
        k1 *= c2
        h1 = h1 ⊻ k1
    end)
end

macro murmur2()
    esc(quote
        k2 *= c2
        k2 = rotl64(k2, 33)
        k2 *= c1
        h2 = h1 ⊻ k2
    end)
end

macro murmur()
    esc(quote
        @murmur1
        h1 = rotl64(h1, 27)
        h1 += h2
        h1 = h1 * 5 + 0x52dce729

        @murmur2
        h2 = rotl64(h2, 31)
        h2 += h1
        h2 = h2 * 5 + 0x38495ab5
    end)
end

# ref: MurmurHash3_x64_128
function Base.hash(seq::BioSequence, seed::UInt64)
    # Mix sequence length so that dna"A" and dna"AA"
    # return the different hash values.
    h1::UInt64 = h2::UInt64 = hash(length(seq), seed)
    c1 = 0x87c37b91114253d5
    c2 = 0x4cf5ad432745937f

    next = bitindex(seq, 1)
    last = bitindex(seq, lastindex(seq) + 1)

    k1::UInt64 = 0
    k2::UInt64 = 0

    # body
    r = offset(next)
    data = seq.data
    if last - next ≥ 128
        if r == 0
            @inbounds while last - next ≥ 128
                j = index(next)
                k1 = data[j]
                k2 = data[j+1]
                @murmur
                next += 128
            end
        else
            x = data[index(next)]
            @inbounds while last - next ≥ 128
                j = index(next)
                y = data[j+1]
                z = data[j+2]
                k1 = x >> r | y << (64 - r)
                k2 = y >> r | z << (64 - r)
                @murmur
                x = z
                next += 128
            end
        end
    end

    # tail
    k1 = 0
    k2 = 0
    if next < last
        x = data[index(next)]
        k1 |= x >> r
        m1 = mask(last - next)
        m2 = mask(max(last - (next + 64), 0))
        next += 64 - r
        if next < last
            y = data[index(next)]
            k1 |= y << (64 - r)
            k2 |= y >> r
            next += 64
            if next < last
                z = data[index(next)]
                k2 |= z << (64 - r)
            end
        end
        k1 &= m1
        k2 &= m2
        @murmur1
        @murmur2
    end

    # finalization
    h1 = h1 ⊻ length(seq)
    h2 = h2 ⊻ length(seq)
    h1 += h2
    h2 += h1
    h1 = fmix64(h1)
    h2 = fmix64(h2)
    h1 += h2
    h2 += h1

    return h1
end


text:
using Gadfly

set_default_plot_size(20cm, 15cm)

#     A spiking neuron model based on:
#     Izhikevich, E. M.
#     Simple Model of Spiking Neurons
#     IEEE TRANSACTIONS ON NEURAL NETWORKS, VOL. 14, NO. 6, NOVEMBER 2003

#         a, b, c, d are the parameters of this model.
#         a: the time scale of the recovery variable.
#         b: the sensitivity of the recovery variable.
#         c: the after-spike reset value of the membrane potential.
#         d: after-spike reset of the recovery variable.

type IZNeuron
    a::Float64
    b::Float64
    c::Float64
    d::Float64
    v::Float64 # membrane potential
    u::Float64 # membrane recovery variable
    has_fired::Bool
    bias::Float64
    current::Float64
    function IZNeuron(bias::Float64=0; a::Float64=0.02,b::Float64=0.2, c::Float64=-65.0, d::Float64=8.0)
       new(a, b, c, d, c, b*c, false, bias, bias)
    end
end

function update!(n::IZNeuron)
#         Advances time in 1 ms.
    n.v += 0.5 * (0.04 * n.v ^ 2. + 5. * n.v + 140. - n.u + n.current)
    n.v += 0.5 * (0.04 * n.v ^ 2. + 5. * n.v + 140. - n.u + n.current)
    n.u += n.a * (n.b * n.v - n.u)
    if n.v > 30
        n.has_fired = true
        n.v = n.c
        n.u += n.d
    else
        n.has_fired = false
        n.current = n.bias
    end
    return n.v
end

function reset!(n::IZNeuron)
#         'Resets all state variables.'
    n.v = n.c
    n.u = n.b * n.v
    n.has_fired = false
    n.current = n.bias
end

function getResponse(V::Vector{Float64}, s::IZNeuron)
    n = length(V)
    data = zeros(n)
    tm = zeros(Int64,n)
    for i = 1:n
        s.current = V[i]
        v = update!(s)
        data[i] = v
        tm[i] = i == 1? 1 : tm[i] = tm[i-1] + 1
    end
    return data,tm
end

# The following parameters produce some known spiking behaviors:
# Regular spiking:
s = IZNeuron(0.,a = 0.02, b = 0.2, c = -65., d = 8.)
# Intrinsically bursting:
s = IZNeuron(0.,a = 0.02, b = 0.2, c = -55., d = 4.)
# Chattering:
s = IZNeuron(0.,a = 0.02, b = 0.2, c = -50., d = 2.)
# Fast spiking:
s = IZNeuron(0.,a = 0.1, b = 0.2, c = -65., d = 2.)
# Thalamo-/cortical:
s = IZNeuron(0.,a = 0.02, b = 0.25, c = -65., d = 0.05)
# Resonator:
s = IZNeuron(0.,a = 0.1, b = 0.25, c = -65., d = 2.)
# Low-threshold spiking:
s = IZNeuron(0.,a = 0.02, b = 0.25, c = -65., d = 2.)

n = 2000
V = ones(n) * 5
response, tm = getResponse(V, s)
vstack(plot(x = tm, y=V, Geom.line), plot(x=tm, y=response, Geom.line))

response[1:10]
tm[1:10]

Ne = 800;                     Ni = 200;
re = rand(Ne);                ri = rand(Ni);
a = [0.02 * ones(Ne),         0.02 + 0.08 * ri]
b = [0.2  * ones(Ne),         0.25 - 0.05 * ri]
c = [-65. + 15.*re.^2,        -65. * ones(Ni)]
d = [  8. - 6. *re.^2,        2. * ones(Ni)]
S = hcat(0.5*rand(Ne+Ni,Ne),  -rand(Ne+Ni,Ni))

v =  -65*ones(Ne+Ni)
u =  b .* v

firings = zeros(Ne+Ni,1000)

for t =  1:1000
    I = [5.*randn(Ne), 2.*randn(Ni)]
    fired = find(x -> x >= 30,v)
    firings[fired,t] = v[fired]
    v[fired] = c[fired]
    u[fired] = u[fired] + d[fired]
    I = I + sum(S[:,fired],2)
    v = v .+ 0.5 * (0.04 * v .^ 2. + 5. .* v + 140. .- u .+ I)
    v = v .+ 0.5 * (0.04 * v .^ 2. + 5. .* v + 140. .- u .+ I)
    u = u .+ a .* (b .* v .- u)
end
firings
# plot(x=[1:1000], y=firings[3,:], Geom.line)
spy(firings)

# class Synapse:
#     """ A synapse indicates the connection strength between two neurons (or itself) """
#     def __init__(self, source, dest, weight):
#         self.__weight = weight
#         self.__source = source
#         self.__dest = dest

#     def advance(self):
#         'Advances time in 1 ms.'
#         if self.__source.has_fired:
#             self.__dest.current += self.__weight # dest.current or dest.__v ?


text:
# tests from ecto/test/ecto/repo_test.exs

module MySchema
using Ecto.Schema
schema("my_schema") do
    field(:x, :string)
    field(:y, :binary)
end
end

module MySchemaNoPK
using Ecto.Schema
primary_key(false)
schema("my_schema") do
    field(:x, :string)
end
end

using Ecto
TestRepo = Ecto.Repo.t()

using Base.Test

@testset "needs schema with primary key field" begin
    schema = %(MySchemaNoPK, x= "abc")
    @test isa(schema, Schema.t)
    @test "abc" == schema.struct[:x]
    @test_throws Ecto.NoPrimaryKeyFieldError TestRepo.update!(schema |> Ecto.Changeset.change, force= true)
    @test_throws Ecto.NoPrimaryKeyFieldError TestRepo.delete!(schema)
end

@testset "works with primary key value" begin
    schema = %(MySchema, id= 1, x= "abc")
    @test isa(schema, Schema.t)
    @test 1 == schema.struct[:id]
    @test "abc" == schema.struct[:x]
    TestRepo.get(MySchema, 123)
    TestRepo.get_by(MySchema, x= "abc")
    TestRepo.update!(schema |> Ecto.Changeset.change, force= true)
    TestRepo.delete!(schema)
end

@testset "works with custom source schema" begin
    schema = %(MySchema, id= 1, x= "abc") |> put_meta(source= "custom_schema")
    @test isa(schema, Schema.t)
    TestRepo.update!(schema |> Ecto.Changeset.change, force= true)
    TestRepo.delete!(schema)
    to_insert = %(MySchema, x= "abc") |> put_meta(source= "custom_schema")
    @test isa(to_insert, Schema.t)
    TestRepo.insert!(to_insert)
end

@testset "fails without primary key value" begin
    schema = %(MySchema, x= "abc")
    @test_throws Ecto.NoPrimaryKeyValueError TestRepo.update!(schema |> Ecto.Changeset.change, force= true)
    @test_throws Ecto.NoPrimaryKeyValueError TestRepo.delete!(schema)
end

@testset "validates schema types" begin
    schema = %(MySchema, x= 123)
    @test_throws Ecto.ChangeError TestRepo.insert!(schema)
end

@testset "validates get" begin
    TestRepo.get(MySchema, 123)
    @test_throws ArgumentError TestRepo.get(MySchema, nothing)
    @test_throws Ecto.Query.CastError TestRepo.get(MySchema, :atom)
    @test_throws Ecto.QueryError TestRepo.get(%(Ecto.Query), :atom)
end

@testset "validates get_by" begin
    TestRepo.get_by(MySchema, id= 123)
    TestRepo.get_by(MySchema, Dict(:id=>123))
    @test_throws Ecto.Query.CastError TestRepo.get_by(MySchema, id= :atom)
end

@testset "validates update_all" begin
    TestRepo.update_all(MySchema, set= [:x "321"])

    e = in(MySchema)
    query = from(e, where= e.x == "123", update= [set [:x "321"]])
    @test isa(query, Query.t)
    TestRepo.update_all(query, [])

    @test_throws ArgumentError TestRepo.update_all(MySchema, [set [:x "321"]], returning= [])
    @test_throws Ecto.QueryError TestRepo.update_all(from(e, select= e), set= [:x "321"])
    @test_throws Ecto.QueryError TestRepo.update_all(from(e, order_by= e.x), set= [:x "321"])
end

@testset "validates delete_all" begin
    TestRepo.delete_all(MySchema)

    e = in(MySchema)
    query = from(e, where= e.x == "123")
    TestRepo.delete_all(query)

    @test_throws ArgumentError TestRepo.delete_all(MySchema, returning= [])
    @test_throws Ecto.QueryError TestRepo.delete_all(from(e, select= e))
    @test_throws Ecto.QueryError TestRepo.delete_all(from(e, order_by= e.x))
end

## Changesets
@testset "insert, update, insert_or_update and delete accepts changesets" begin
    valid = Ecto.Changeset.cast(%(MySchema, id= 1), Dict(), [])
    @test (:ok, %(MySchema)) == TestRepo.insert(valid)
    @test (:ok, %(MySchema)) == TestRepo.update(valid)
    @test (:ok, %(MySchema)) == TestRepo.insert_or_update(valid)
    @test (:ok, %(MySchema)) == TestRepo.delete(valid)
end

@testset "insert, update, insert_or_update and delete errors on invalid changeset" begin
    invalid = %(Ecto.Changeset.t, valid= false, data= %(MySchema))

    insert = %(invalid, action= :insert, repo= TestRepo.repo)
    @test (:error, insert) == TestRepo.insert(invalid)

    update = %(invalid, action= :update, repo= TestRepo.repo)
    @test (:error, update) == TestRepo.update(invalid)

    update = %(invalid, action= :insert, repo= TestRepo.repo)
    @test (:error, update) == TestRepo.insert_or_update(invalid)

    delete = %(invalid, action= :delete, repo= TestRepo.repo)
    @test (:error, delete) == TestRepo.delete(invalid)
end

@testset "insert!, update!, insert_or_update! and delete! fail on invalid changeset" begin
    invalid = %(Ecto.Changeset.t, valid= false, data= %(MySchema))

    @test_throws Ecto.InvalidChangesetError TestRepo.insert!(invalid)
    @test_throws Ecto.InvalidChangesetError TestRepo.update!(invalid)
    @test_throws Ecto.InvalidChangesetError TestRepo.insert_or_update!(invalid)
    @test_throws Ecto.InvalidChangesetError TestRepo.delete!(invalid)
end

@testset "insert!, update! and delete! fail on changeset without data" begin
    invalid = %(Ecto.Changeset.t, valid= true, data= nothing)

    @test_throws ArgumentError TestRepo.insert!(invalid)
    @test_throws ArgumentError TestRepo.update!(invalid)
    @test_throws ArgumentError TestRepo.delete!(invalid)
end

@testset "insert!, update!, insert_or_update! and delete! fail on changeset with wrong action" begin
    invalid = %(Ecto.Changeset.t, valid= true, data= %(MySchema), action= :other)

    @test_throws ArgumentError TestRepo.insert!(invalid)
    @test_throws ArgumentError TestRepo.update!(invalid)
    @test_throws ArgumentError TestRepo.insert_or_update!(invalid)
    @test_throws ArgumentError TestRepo.delete!(invalid)
end

@testset "insert_or_update uses the correct action" begin
    built  = Ecto.Changeset.cast(%(MySchema, y= "built"), Dict(), [])
    @test isa(built, Ecto.Changeset.t)
    loaded =
        %(MySchema, y= "loaded") |> TestRepo.insert! |> Ecto.Changeset.cast(Schema.Assoc(y= "updated"), [:y])
    #:insert

    TestRepo.insert_or_update(built)
    #:insert

    TestRepo.insert_or_update(loaded)
    #:update
end


text:
"""
objectivefunction(circuitfunc,measurements,frequencies)
Returns the Least squares objective function between the measurements and a circuit's simulated
impedance for a given set of circuit parameters at a specified range of experimental frequencies.
This function is used to optimize the circuit parameters.

"""
function objectivefunction(circuitfunc,measurements,frequencies)
    function objective(x)
        model_output = [circuitfunc(x,fr) for fr in frequencies]
        return  mean((abs.(measurements - model_output).^2)./(abs.(measurements).^2 .+ abs.(model_output).^2))
    end
    return objective
end

function objectivefunction1(circuitfunc,measurements,frequencies) #fitness measure from Ramos et al.
    function objective(x)
        model_output = [circuitfunc(x,fr) for fr in frequencies]
        return  mean((abs.(measurements - model_output).^2)./(abs.(model_output).^2))
    end
    return objective
end

#Some alternate objective functions.

function objectivefunction2(circuitfunc,measurements,frequencies)
    function objective(x)
        model_output = [circuitfunc(x,fr) for fr in frequencies]
        rme,ime,rmo,imo = real(measurements),imag(measurements),real(model_output),imag(model_output)
        return sum((rme.-rmo).^2 .+ (ime.-imo).^2)
    end
    return objective
end

function objectivefunction3(circuitfunc,measurements,frequencies) #fitness measure from Ramos et al.
    function objective(x)
        model_output = [circuitfunc(x,fr) for fr in frequencies]
        return  mean((abs.(measurements - model_output).^2)./(abs.(measurements).^2))
    end
    return objective
end

function objectivefunction4(circuitfunc,measurements,frequencies) #fitness measure from Ramos et al.
    function objective(x)
        model_output = [circuitfunc(x,fr) for fr in frequencies]
        return  mean((abs.(measurements - model_output).^2))
    end
    return objective
end

function objectivefunction5(circuitfunc,measurements,frequencies) #fitness measure from Lempka et al.
    function objective(x)
        model_output = [circuitfunc(x,fr) for fr in frequencies]
        rme,ime,rmo,imo = real(measurements),imag(measurements),real(model_output),imag(model_output)
        return sum(((rme.-rmo).^2)./(rmo.^2 .+rme.^2) .+ ((ime.-imo).^2)./(imo.^2 .+rme.^2))
    end
    return objective
end

function objectivefunction6(circuitfunc,measurements,frequencies) #fitness measure from Lempka et al.
    function objective(x)
        model_output = [circuitfunc(x,fr) for fr in frequencies]
        rme,ime,rmo,imo = real(measurements),imag(measurements),real(model_output),imag(model_output)
        return sum(((rme.-rmo).^2)./(rmo.^2) .+ ((ime.-imo).^2)./(imo.^2))
    end
    return objective
end

function objectivefunction7(circuitfunc,measurements,frequencies)
    function objective(x)
        model_output = [circuitfunc(x,fr) for fr in frequencies]
        return  mean((abs.(measurements - model_output).^2)./(abs.(measurements).^2 .+ abs.(model_output).^2))
    end
    return objective
end


text:
@everywhere using O_structs
@everywhere using B_structs


text:
"""
    ZeroKernel()

Zero kernel.

# Definition

For inputs ``x, x'``, the zero kernel is defined as
```math
k(x, x') = 0.
```
The output type depends on ``x`` and ``x'``.

See also: [`ConstantKernel`](@ref)
"""
struct ZeroKernel <: SimpleKernel end

# SimpleKernel interface
kappa(::ZeroKernel, ::Real) = false
metric(::ZeroKernel) = Delta()

# Optimizations
(::ZeroKernel)(x, y) = false
kernelmatrix(::ZeroKernel, x::AbstractVector) = Falses(length(x), length(x))
function kernelmatrix(::ZeroKernel, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    return Falses(length(x), length(y))
end
function kernelmatrix!(K::AbstractMatrix, ::ZeroKernel, x::AbstractVector)
    validate_inplace_dims(K, x)
    return fill!(K, zero(eltype(K)))
end
function kernelmatrix!(
    K::AbstractMatrix, ::ZeroKernel, x::AbstractVector, y::AbstractVector
)
    validate_inplace_dims(K, x, y)
    return fill!(K, zero(eltype(K)))
end
kernelmatrix_diag(::ZeroKernel, x::AbstractVector) = Falses(length(x))
function kernelmatrix_diag(::ZeroKernel, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    return Falses(length(x))
end
function kernelmatrix_diag!(K::AbstractVector, ::ZeroKernel, x::AbstractVector)
    validate_inplace_dims(K, x)
    return fill!(K, zero(eltype(K)))
end
function kernelmatrix_diag!(
    K::AbstractVector, ::ZeroKernel, x::AbstractVector, y::AbstractVector
)
    validate_inplace_dims(K, x, y)
    return fill!(K, zero(eltype(K)))
end

Base.show(io::IO, ::ZeroKernel) = print(io, "Zero Kernel")

"""
    WhiteKernel()

White noise kernel.

# Definition

For inputs ``x, x'``, the white noise kernel is defined as
```math
k(x, x') = \\delta(x, x').
```
"""
struct WhiteKernel <: SimpleKernel end

"""
    EyeKernel()

Alias of [`WhiteKernel`](@ref).
"""
const EyeKernel = WhiteKernel

kappa(κ::WhiteKernel, δₓₓ::Real) = δₓₓ

metric(::WhiteKernel) = Delta()

Base.show(io::IO, ::WhiteKernel) = print(io, "White Kernel")

"""
    ConstantKernel(; c::Real=1.0)

Kernel of constant value `c`.

# Definition

For inputs ``x, x'``, the kernel of constant value ``c \\geq 0`` is defined as
```math
k(x, x') = c.
```

See also: [`ZeroKernel`](@ref)
"""
struct ConstantKernel{Tc<:Real} <: SimpleKernel
    c::Vector{Tc}

    function ConstantKernel(; c::Real=1.0)
        @check_args(ConstantKernel, c, c >= zero(c), "c ≥ 0")
        return new{typeof(c)}([c])
    end
end

@functor ConstantKernel

# SimpleKernel interface
kappa(κ::ConstantKernel, ::Real) = only(κ.c)
metric(::ConstantKernel) = Delta()

# Optimizations
(k::ConstantKernel)(x, y) = only(k.c)
kernelmatrix(k::ConstantKernel, x::AbstractVector) = Fill(only(k.c), length(x), length(x))
function kernelmatrix(k::ConstantKernel, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    return Fill(only(k.c), length(x), length(y))
end
function kernelmatrix!(K::AbstractMatrix, k::ConstantKernel, x::AbstractVector)
    validate_inplace_dims(K, x)
    return fill!(K, only(k.c))
end
function kernelmatrix!(
    K::AbstractMatrix, k::ConstantKernel, x::AbstractVector, y::AbstractVector
)
    validate_inplace_dims(K, x, y)
    return fill!(K, only(k.c))
end
kernelmatrix_diag(k::ConstantKernel, x::AbstractVector) = Fill(only(k.c), length(x))
function kernelmatrix_diag(k::ConstantKernel, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    return Fill(only(k.c), length(x))
end
function kernelmatrix_diag!(K::AbstractVector, k::ConstantKernel, x::AbstractVector)
    validate_inplace_dims(K, x)
    return fill!(K, only(k.c))
end
function kernelmatrix_diag!(
    K::AbstractVector, k::ConstantKernel, x::AbstractVector, y::AbstractVector
)
    validate_inplace_dims(K, x, y)
    return fill!(K, only(k.c))
end

Base.show(io::IO, κ::ConstantKernel) = print(io, "Constant Kernel (c = ", only(κ.c), ")")


text:
using .Tokenize.Tokens: Kind

include("token_kinds.jl")

"""
    K"s"

The full kind of a string "s".  For example, K")" is the kind of the
right parenthesis token.

Naming rules:
* Kinds which correspond to exactly one textural form are represented with that
  text. This includes keywords like K"for" and operators like K"*".
* Kinds which represent many textural forms have UpperCamelCase names. This
  includes kinds like K"Identifier" and K"Comment".
* Kinds which exist merely as delimiters are all uppercase
"""
macro K_str(str)
    get(_str_to_kind, str) do
        error("unknown token kind K$(repr(str))")
    end
end

"""
A set of kinds which can be used with the `in` operator.  For example

    k in KSet`+ - *`
"""
macro KSet_cmd(str)
    kinds = [get(_str_to_kind, s) do
         error("unknown token kind KSet`$(repr(str)[2:end-1])`")
    end
    for s in split(str)]

    quote
        ($(kinds...),)
    end
end

kind(k::Kind) = k
kind(raw::TzTokens.RawToken) = TzTokens.exactkind(raw)

# Some renaming for naming consistency
is_literal(k) = TzTokens.isliteral(kind(k))
is_keyword(k) = TzTokens.iskeyword(kind(k))
is_contextual_keyword(k) = TzTokens.iscontextualkeyword(kind(k))
is_operator(k) = TzTokens.isoperator(kind(k))
is_word_operator(k) = TzTokens.iswordoperator(kind(k))

# Predicates for operator precedence
# FIXME: Review how precedence depends on dottedness, eg
# https://github.com/JuliaLang/julia/pull/36725
is_prec_assignment(t)  = K"BEGIN_ASSIGNMENTS" < kind(t) < K"END_ASSIGNMENTS"
is_prec_pair(t)        = K"BEGIN_PAIRARROW"   < kind(t) < K"END_PAIRARROW"
is_prec_conditional(t) = K"BEGIN_CONDITIONAL" < kind(t) < K"END_CONDITIONAL"
is_prec_arrow(t)       = K"BEGIN_ARROW"       < kind(t) < K"END_ARROW"
is_prec_lazy_or(t)     = K"BEGIN_LAZYOR"      < kind(t) < K"END_LAZYOR"
is_prec_lazy_and(t)    = K"BEGIN_LAZYAND"     < kind(t) < K"END_LAZYAND"
is_prec_comparison(t)  = K"BEGIN_COMPARISON"  < kind(t) < K"END_COMPARISON"
is_prec_pipe(t)        = K"BEGIN_PIPE"        < kind(t) < K"END_PIPE"
is_prec_colon(t)       = K"BEGIN_COLON"       < kind(t) < K"END_COLON"
is_prec_plus(t)        = K"BEGIN_PLUS"        < kind(t) < K"END_PLUS"
is_prec_bitshift(t)    = K"BEGIN_BITSHIFTS"   < kind(t) < K"END_BITSHIFTS"
is_prec_times(t)       = K"BEGIN_TIMES"       < kind(t) < K"END_TIMES"
is_prec_rational(t)    = K"BEGIN_RATIONAL"    < kind(t) < K"END_RATIONAL"
is_prec_power(t)       = K"BEGIN_POWER"       < kind(t) < K"END_POWER"
is_prec_decl(t)        = K"BEGIN_DECL"        < kind(t) < K"END_DECL"
is_prec_where(t)       = K"BEGIN_WHERE"       < kind(t) < K"END_WHERE"
is_prec_dot(t)         = K"BEGIN_DOT"         < kind(t) < K"END_DOT"
is_prec_unicode_ops(t) = K"BEGIN_UNICODE_OPS" < kind(t) < K"END_UNICODE_OPS"

is_prec_pipe_lt(t)     = kind(t) == K"<|"
is_prec_pipe_gt(t)     = kind(t) == K"|>"

#=
# Sholuld we optimize membership a bit by unrolling?
@generated function in(k::Kind, t::NTuple{N,Kind}) where {N}
    ex = :(k === t[1])
    for i = 2:N
        ex = :($ex || k === t[$i])
    end
    quote
        $ex
    end
end
=#

function is_syntax_kind(t)
    K"BEGIN_SYNTAX_KINDS" < kind(t) < K"END_SYNTAX_KINDS"
end

function is_identifier(k)
    kind(k) == K"Identifier"
end

function is_macro_name(k)
    K"BEGIN_MACRO_NAMES" < kind(k) < K"END_MACRO_NAMES"
end

function is_number(t)
    kind(t) in (K"Integer", K"BinInt", K"HexInt", K"OctInt", K"Float")
end

function is_string_delim(t)
    kind(t) in (K"\"", K"\"\"\"")
end

function is_radical_op(t)
    kind(t) in (K"√", K"∛", K"∜")
end

function is_whitespace(t)
    kind(t) in (K"Whitespace", K"NewlineWs")
end

"""
Return the string representation of a token kind, or `nothing` if the kind
represents a class of tokens like K"Identifier".

When `unique=true` only return a string when the kind uniquely defines the
corresponding input token, otherwise return `nothing`.  When `unique=false`,
return the name of the kind.

TODO: Replace `untokenize()` with `Base.string()`?
"""
function untokenize(k::Kind; unique=true)
    get(unique ? _kind_to_str_unique : _kind_to_str, k, nothing)
end


text:
# -------------------------------------------------------------------
#* Distillation Column model with:
#*
#*	- NumberOfTrays like tray;
#*	- a vessel in the bottom of column;
#*	- a splitter who separate the bottom product and the stream to reboiler;
#*	- steady state reboiler (thermosyphon);
#*	- a steady state condenser with subcooling;
#*	- a vessel drum (layed cilinder);
#*	- a splitter which separate reflux and distillate;
#*	- a pump in reflux stream.
#*
#* ------------------------------------------------------------------
type User_Distillation_thermosyphon_subcooling
	User_Distillation_thermosyphon_subcooling()=begin
		new(
			User_Section_ColumnBasic(),
			DanaSwitcher(Dict{Symbol,Any}(
				:Valid=>["on", "off"],
				:Hidden=>true,
				:Default=>"on"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour Outlet in the section",
				:PosX=>1,
				:PosY=>0.41,
				:Protected=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid Outlet in the section",
				:PosX=>1,
				:PosY=>0.515,
				:Protected=>true
			)),
			condenserSteady(),
			tank_cylindrical(),
			splitter(),
			pump(),
			reboilerSteady(),
			tank(),
			splitter(),
			DanaReal(),
			energy_stream(Dict{Symbol,Any}(
				:Brief=>"Heat supplied to Condenser",
				:Hidden=>true
			)),
			energy_stream(Dict{Symbol,Any}(
				:Brief=>"Heat supplied to Reboiler",
				:Hidden=>true
			)),
			energy_stream(Dict{Symbol,Any}(
				:Brief=>"Heat supplied to Bottom Vessel",
				:Hidden=>true
			)),
			energy_stream(Dict{Symbol,Any}(
				:Brief=>"Heat supplied to Top Vessel",
				:Hidden=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream From Top Splitter",
				:PosX=>1,
				:PosY=>0.30
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Connector for Liquid outlet stream From Top Splitter",
				:Hidden=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream From Bottom Splitter",
				:PosX=>1,
				:PosY=>1
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Connector for Liquid outlet stream From Reboiler",
				:Hidden=>true
			)),
			[
				:(ConnectorSplitterTop.T = LiquidDistillate.T),
				:(ConnectorSplitterTop.P = LiquidDistillate.P),
				:(ConnectorSplitterTop.F = LiquidDistillate.F),
				:(ConnectorSplitterTop.z = LiquidDistillate.z),
				:(ConnectorSplitterBottom.T = BottomProduct.T),
				:(ConnectorSplitterBottom.P = BottomProduct.P),
				:(ConnectorSplitterBottom.F = BottomProduct.F),
				:(ConnectorSplitterBottom.z = BottomProduct.z),
				:(VapourDrawOff.F*_base_1.VapSideTrayIndex= trays.VapourSideStream.F),
				:(VapourDrawOff.T = trays(_base_1.VapourSideStreamLocation).VapourSideStream.T),
				:(VapourDrawOff.P = trays(_base_1.VapourSideStreamLocation).VapourSideStream.P),
				:(VapourDrawOff.z = trays(_base_1.VapourSideStreamLocation).VapourSideStream.z),
				:(LiquidDrawOff.F*_base_1.LiqSideTrayIndex= trays.LiquidSideStream.F),
				:(LiquidDrawOff.T = trays(_base_1.LiquidSideStreamLocation).LiquidSideStream.T),
				:(LiquidDrawOff.P = trays(_base_1.LiquidSideStreamLocation).LiquidSideStream.P),
				:(LiquidDrawOff.z = trays(_base_1.LiquidSideStreamLocation).LiquidSideStream.z),
				:(VapourDrawOffFlow = VapourDrawOff.F),
				:(LiquidDrawOffFlow = LiquidDrawOff.F),
				:(CondenserUnity.InletV.F*trays(1).vV = alfaTopo * _base_1.Ah * sqrt(2*(trays(1).OutletV.P - CondenserUnity.OutletL.P + 1e-8 * "atm") / (_base_1.alfa*trays(1).rhoV))),
				:(CondenserUnity.InletV.F = 0 * "mol/s"),
			],
			[
				"","","","","","","","","","","","","","","","","","","","",
			],
			[:CondenserVapourFlow,],
			[:VapourDrawOff,:LiquidDrawOff,:CondenserUnity,:TopVessel,:TopSplitter,:PumpUnity,:ReboilerUnity,:BottomVessel,:BottomSplitter,:alfaTopo,:HeatToCondenser,:HeatToReboiler,:HeatToBottomVessel,:HeatToTopVessel,:LiquidDistillate,:ConnectorSplitterTop,:BottomProduct,:ConnectorSplitterBottom,]
		)
	end
	_base_1::User_Section_ColumnBasic
	CondenserVapourFlow::DanaSwitcher
	VapourDrawOff::vapour_stream
	LiquidDrawOff::liquid_stream
	CondenserUnity::condenserSteady
	TopVessel::tank_cylindrical
	TopSplitter::splitter
	PumpUnity::pump
	ReboilerUnity::reboilerSteady
	BottomVessel::tank
	BottomSplitter::splitter
	alfaTopo::DanaReal
	HeatToCondenser::energy_stream
	HeatToReboiler::energy_stream
	HeatToBottomVessel::energy_stream
	HeatToTopVessel::energy_stream
	LiquidDistillate::liquid_stream
	ConnectorSplitterTop::stream
	BottomProduct::liquid_stream
	ConnectorSplitterBottom::stream
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export User_Distillation_thermosyphon_subcooling
function setEquationFlow(in::User_Distillation_thermosyphon_subcooling)
	# Top Splitter Connector Equations
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	# Bottom Splitter Connector Equations
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	let switch=CondenserVapourFlow
		if CondenserVapourFlow==CondenserUnity.InletV.F < 1e-6 * "kmol/h"
			set(switch,"off")
		end
		if CondenserVapourFlow==trays(1).OutletV.P > CondenserUnity.OutletL.P + 1e-1 * "atm"
			set(switch,"on")
		end
		if switch=="on"
			addEquation(19)
		elseif switch=="off"
			addEquation(20)
		end
	end
end
function atributes(in::User_Distillation_thermosyphon_subcooling,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/DistillationThermosyphonSubcooling"
	fields[:Brief]="Model of a distillation column with steady condenser and steady reboiler."
	fields[:Info]="== Specify ==
* the feed stream of each tray (Inlet);
* the Murphree eficiency for each tray Emv;
* the pump head;
* the condenser pressure drop;
* the heat supllied in top and bottom tanks;
* the heat supllied in condenser and reboiler;
* the Outlet1 flow in the bottom splitter (spbottom.Outlet1.F) that corresponds to the bottom product;
* both  top splitter outlet flows OR one of the splitter outlet flows and the splitter frac.

== Initial Conditions ==
* the trays temperature (OutletL.T);
* the trays liquid level (Level) OR the trays liquid flow (OutletL.F);
* (NoComps - 1) OutletL (OR OutletV) compositions for each tray;

* the top tank temperature (OutletL.T);
* the top tank liquid level (Level);
* (NoComps - 1) OutletL (OR OutletV) compositions;

* the bottom tank temperature (OutletL.T);
* the bottom tank liquid level (Level);
* (NoComps - 1) OutletL (OR OutletV) compositions.
"
	drive!(fields,_)
	return fields
end
User_Distillation_thermosyphon_subcooling(_::Dict{Symbol,Any})=begin
	newModel=User_Distillation_thermosyphon_subcooling()
	newModel.attributes=atributes(newModel,_)
	newModel
end