text:
#
# Code generation for boxing and unboxing.
#

emit_box!(cg::CodeCtx, x) = codegen!(cg, x)
emit_box!(cg::CodeCtx, ::Type{T}) where {T} = cg.datatype[T]

emit_map = Dict(
    Bool => :jl_box_bool,
    Int8 => :jl_box_int8,
    Int16 => :jl_box_int16,
    Int32 => :jl_box_int32,
    Int64 => :jl_box_int64,
    UInt8 => :jl_box_uint8,
    UInt16 => :jl_box_uint16,
    UInt32 => :jl_box_uint32,
    UInt64 => :jl_box_uint64,
    Float32 => :jl_box_float32,
    Float64 => :jl_box_float64)

bitsT = Union{Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, UInt64, Float32, Float64}
emit_box!(cg::CodeCtx, x::T) where {T <: bitsT}  = LLVM.call!(cg.builder, cg.extern[emit_map[T]], LLVM.Value[codegen!(cg, x)])
emit_box!(cg::CodeCtx, x::Bool) = LLVM.call!(cg.builder, cg.extern[:jl_box_bool], LLVM.Value[emit_val!(cg, codegen!(cg, x))])

function emit_box!(cg::CodeCtx, x::SlotNumber)
    slottype = cg.code_info.slottypes[x.id]
    @debug "$(cg.name): boxing Slot $(x.id)" slottype
    if haskey(emit_map, slottype)
        return LLVM.call!(cg.builder, cg.extern[emit_map[slottype]], LLVM.Value[codegen!(cg, x)])
    else
        return codegen!(cg, x)
    end
end

function emit_box!(cg::CodeCtx, x::SSAValue)
    ssatype = cg.code_info.ssavaluetypes[x.id+1] # I think the first position is a placeholder
    v = codegen!(cg, x)
    @debug "$(cg.name): boxing SSAValue $(x.id)" ssatype
    if isbits(ssatype)
        return LLVM.call!(cg.builder, cg.extern[emit_map[ssatype]], LLVM.Value[v])
    else
        return v
    end
end

emit_box!(cg::CodeCtx, e::Expr) = LLVM.call!(cg.builder, cg.extern[emit_map[e.typ]], LLVM.Value[codegen!(cg, e)])

function emit_box!(cg::CodeCtx, v::LLVM.Value)
    f!(s::Symbol) = LLVM.call!(cg.builder, cg.extern[s], LLVM.Value[v])
    t = LLVM.llvmtype(v)
    return t == int8_t ? f!(:jl_box_int8) :
           t == int16_t ? f!(:jl_box_int16) :
           t == int32_t ? f!(:jl_box_int32) :
           t == int64_t ? f!(:jl_box_int64) :
           t == float32_t ? f!(:jl_box_float32) :
           t == float64_t ? f!(:jl_box_float64) : error("not supported")
end

# function emit_box!(cg::CodeCtx, ::Type{T}, v) where T <: Tuple
# #   params = jl_svec(n_args, arg1_type, arg2_type, ...);
# #   tuple_type = jl_apply_tuple_type(params);
# #   result = jl_new_struct(tuple_type, jl_box...(arg1), jl_box...(arg2), ...);
#     if !isbits(T)
#         return v
#     end
#     println("Boxing a bits-type tuple")
#     n = length(T.parameters)
#     dt = LLVM.call!(cg.builder, cg.extern[:jl_apply_tuple_type],
#                     LLVM.call!(cg.builder, cg.extern[:jl_svec], codegen!(cg, Cint(n)),
#                                load_and_emit_datatype!.(cg, collect(T.parameters))...))
#     llvmargs = LLVM.Value[dt]
#     for i in 1:n
#         push!(llvmargs, codegen!(cg, x)) # need gep, load, store
#     end
#     return LLVM.call!(cg.builder, cg.extern[:jl_new_struct], llvmargs)
# end

# function emit_box!(cg::CodeCtx, @nospecialize(x::T)) where T
#     isa(x, Type)         && return cg.datatype[x]
#     v = codegen!(cg, x)
#     @debug "$(cg.name): boxing " x v
#     T == Any && return v
#     T <: Base.BitInteger && return emit_box!(cg, T, v)
#     T <: Base.IEEEFloat  && return emit_box!(cg, T, v)
#     T == Bool            && return emit_box!(cg, T, v)
#     # T == SSAValue        && return emit_box!(cg, T, codegen!(cg, x.id))
#     # T == SlotNumber      && return emit_box!(cg, T, codegen!(cg, x.id))
#     if T == SlotNumber
#         slottype = cg.code_info.slottypes[x.id]
#         # if isbits(slottype)
#             return emit_box!(cg, slottype, v)
#         # else
#         #     return v
#         # end
#     end
#     if T == SSAValue
#         ssatype = cg.code_info.ssavaluetypes[x.id+1] # I think the first position is a placeholder
#         @debug "$(cg.name): boxing SSAValue $(x.id)" ssatype
#         if isbits(ssatype)
#             return emit_box!(cg, ssatype, v)
#         else
#             return v
#         end
#     end
#     if T == Expr
#         # if isbits(x.typ)
#             return emit_box!(cg, x.typ, v)
#         # else
#         #     return v
#         # end
#     end
#     if T <: LLVM.Value   # BROKEN
#         return emit_box!(cg, LLVM.llvmtype(x), x)
#     end
#     return v
#     error("Boxing of $T not supported")
# end
# # emit_box!(cg::CodeCtx, v, t::LLVM.PointerType) = v

function emit_unbox!(cg::CodeCtx, v, ::Type{T}) where T
    t = LLVM.llvmtype(v)
    if t == jl_value_t_ptr
        T == Bool  && return LLVM.call!(cg.builder, cg.extern[:jl_unbox_bool], LLVM.Value[v])
        T == Int8  && return LLVM.call!(cg.builder, cg.extern[:jl_unbox_int8], LLVM.Value[v])
        T == Int16 && return LLVM.call!(cg.builder, cg.extern[:jl_unbox_int16], LLVM.Value[v])
        T == Int32 && return LLVM.call!(cg.builder, cg.extern[:jl_unbox_int32], LLVM.Value[v])
        T == Int64 && return LLVM.call!(cg.builder, cg.extern[:jl_unbox_int64], LLVM.Value[v])
        T == Float32 && return LLVM.call!(cg.builder, cg.extern[:jl_unbox_float32], LLVM.Value[v])
        T == Float64 && return LLVM.call!(cg.builder, cg.extern[:jl_unbox_float64], LLVM.Value[v])
        if T <: Tuple && isbits(T)
            bcp = LLVM.bitcast!(cg.builder, v, LLVM.PointerType(llvmtype(T)))
            return LLVM.load!(cg.builder, bcp)
        end
    end
    return v
end

Base.isbits(cg::CodeCtx, x) = false
Base.isbits(cg::CodeCtx, x::SlotNumber) = isbits(cg.code_info.slottypes[x.id])
Base.isbits(cg::CodeCtx, x::SSAValue) = isbits(cg.code_info.ssavaluetypes[x.id + 1])
Base.isbits(cg::CodeCtx, x::Expr) = isbits(x.typ)
Base.fieldnames(cg::CodeCtx, x::SlotNumber) = fieldnames(cg.code_info.slottypes[x.id])
Base.fieldnames(cg::CodeCtx, x::SSAValue) = fieldnames(cg.code_info.ssavaluetypes[x.id + 1])
Base.fieldnames(cg::CodeCtx, x::Expr) = fieldnames(x.typ)
_typeof(cg::CodeCtx, x::SlotNumber) = cg.code_info.slottypes[x.id]
_typeof(cg::CodeCtx, x::SSAValue) = cg.code_info.ssavaluetypes[x.id + 1]
_typeof(cg::CodeCtx, x::Expr) = x.typ
_typeof(cg::CodeCtx, x) = Base.typeof(x)
_typeof(cg::CodeCtx, x::GlobalRef) = Any


text:
module HTTP

export startwrite, startread, closewrite, closeread, stack, insert, insert_default!,
    remove_default!, AWS4AuthLayer, BasicAuthLayer, CanonicalizeLayer, ConnectionPoolLayer,
    ContentTypeDetectionLayer, DebugLayer, ExceptionLayer, MessageLayer, RedirectLayer,
    RetryLayer, StreamLayer, TimeoutLayer

const DEBUG_LEVEL = Ref(0)

Base.@deprecate escape escapeuri

using Base64, Sockets, Dates
using URIs

include("debug.jl")

include("Pairs.jl")                    ;using .Pairs
include("IOExtras.jl")                 ;using .IOExtras
include("Strings.jl")
include("sniff.jl")
include("multipart.jl")
include("Parsers.jl")                  ;import .Parsers: Headers, Header,
                                                         ParseError
include("ConnectionPool.jl")
include("Messages.jl")                 ;using .Messages
include("cookies.jl")                  ;using .Cookies
include("Streams.jl")                  ;using .Streams
include("layers.jl")                   ;using .Layers

const nobody = UInt8[]

"""

    HTTP.request(method, url [, headers [, body]]; <keyword arguments>]) -> HTTP.Response

Send a HTTP Request Message and receive a HTTP Response Message.

e.g.
```julia
r = HTTP.request("GET", "http://httpbin.org/ip")
println(r.status)
println(String(r.body))
```

`headers` can be any collection where
`[string(k) => string(v) for (k,v) in headers]` yields `Vector{Pair}`.
e.g. a `Dict()`, a `Vector{Tuple}`, a `Vector{Pair}` or an iterator.

`body` can take a number of forms:

 - a `String`, a `Vector{UInt8}` or any `T` accepted by `write(::IO, ::T)`
 - a collection of `String` or `AbstractVector{UInt8}` or `IO` streams
   or items of any type `T` accepted by `write(::IO, ::T...)`
 - a readable `IO` stream or any `IO`-like type `T` for which
   `eof(T)` and `readavailable(T)` are defined.

The `HTTP.Response` struct contains:

 - `status::Int16` e.g. `200`
 - `headers::Vector{Pair{String,String}}`
    e.g. ["Server" => "Apache", "Content-Type" => "text/html"]
 - `body::Vector{UInt8}`, the Response Body bytes
    (empty if a `response_stream` was specified in the `request`).

Functions `HTTP.get`, `HTTP.put`, `HTTP.post` and `HTTP.head` are defined as
shorthand for `HTTP.request("GET", ...)`, etc.

`HTTP.request` and `HTTP.open` also accept optional keyword parameters.

e.g.
```julia
HTTP.request("GET", "http://httpbin.org/ip"; retries=4, cookies=true)

HTTP.get("http://s3.us-east-1.amazonaws.com/"; aws_authorization=true)

conf = (readtimeout = 10,
        pipeline_limit = 4,
        retry = false,
        redirect = false)

HTTP.get("http://httpbin.org/ip"; conf..)
HTTP.put("http://httpbin.org/put", [], "Hello"; conf..)
```

URL options

 - `query = nothing`, replaces the query part of `url`.

Streaming options

 - `response_stream = nothing`, a writeable `IO` stream or any `IO`-like
    type `T` for which `write(T, AbstractVector{UInt8})` is defined.
 - `verbose = 0`, set to `1` or `2` for extra message logging.

Connection Pool options

 - `connection_limit = 8`, number of concurrent connections to each host:port.
 - `pipeline_limit = 16`, number of concurrent requests per connection.
 - `reuse_limit = nolimit`, number of times a connection is reused after the
                            first request.
 - `socket_type = TCPSocket`

Timeout options

 - `readtimeout = 0`, close the connection if no data is received for this many
   seconds. Use `readtimeout = 0` to disable.

Retry options

 - `retry = true`, retry idempotent requests in case of error.
 - `retries = 4`, number of times to retry.
 - `retry_non_idempotent = false`, retry non-idempotent requests too. e.g. POST.

Redirect options

 - `redirect = true`, follow 3xx redirect responses.
 - `redirect_limit = 3`, number of times to redirect.
 - `forwardheaders = true`, forward original headers on redirect.

Status Exception options

 - `status_exception = true`, throw `HTTP.StatusError` for response status >= 300.

SSLContext options

 - `require_ssl_verification = false`, pass `MBEDTLS_SSL_VERIFY_REQUIRED` to
   the mbed TLS library.
   ["... peer must present a valid certificate, handshake is aborted if
     verification failed."](https://tls.mbed.org/api/ssl_8h.html#a5695285c9dbfefec295012b566290f37)
 - `sslconfig = SSLConfig(require_ssl_verification)`

Basic Authentication options

 - Basic authentication is detected automatically from the provided url's `userinfo` (in the form `scheme://user:password@host`)
   and adds the `Authorization: Basic` header

AWS Authentication options

 - `aws_authorization = false`, enable AWS4 Authentication.
 - `aws_service = split(url.host, ".")[1]`
 - `aws_region = split(url.host, ".")[2]`
 - `aws_access_key_id = ENV["AWS_ACCESS_KEY_ID"]`
 - `aws_secret_access_key = ENV["AWS_SECRET_ACCESS_KEY"]`
 - `aws_session_token = get(ENV, "AWS_SESSION_TOKEN", "")`
 - `body_sha256 = digest(MD_SHA256, body)`,
 - `body_md5 = digest(MD_MD5, body)`,

Cookie options

 - `cookies::Union{Bool, Dict{<:AbstractString, <:AbstractString}} = false`, enable cookies, or alternatively,
        pass a `Dict{AbstractString, AbstractString}` of name-value pairs to manually pass cookies
 - `cookiejar::Dict{String, Set{Cookie}}=default_cookiejar`,

Canonicalization options

 - `canonicalize_headers = false`, rewrite request and response headers in
   Canonical-Camel-Dash-Format.

Proxy options

 - `proxy = proxyurl`, pass request through a proxy given as a url

Alternatively, HTTP.jl also respects the `http_proxy`, `https_proxy`, and `no_proxy`
environment variables; if set, they will be used automatically when making requests.

## Request Body Examples

String body:
```julia
HTTP.request("POST", "http://httpbin.org/post", [], "post body data")
```

Stream body from file:
```julia
io = open("post_data.txt", "r")
HTTP.request("POST", "http://httpbin.org/post", [], io)
```

Generator body:
```julia
chunks = ("chunk\$i" for i in 1:1000)
HTTP.request("POST", "http://httpbin.org/post", [], chunks)
```

Collection body:
```julia
chunks = [preamble_chunk, data_chunk, checksum(data_chunk)]
HTTP.request("POST", "http://httpbin.org/post", [], chunks)
```

`open() do io` body:
```julia
HTTP.open("POST", "http://httpbin.org/post") do io
    write(io, preamble_chunk)
    write(io, data_chunk)
    write(io, checksum(data_chunk))
end
```

## Response Body Examples

String body:
```julia
r = HTTP.request("GET", "http://httpbin.org/get")
println(String(r.body))
```

Stream body to file:
```julia
io = open("get_data.txt", "w")
r = HTTP.request("GET", "http://httpbin.org/get", response_stream=io)
close(io)
println(read("get_data.txt"))
```

Stream body through buffer:
```julia
io = Base.BufferStream()
@async while !eof(io)
    bytes = readavailable(io))
    println("GET data: \$bytes")
end
r = HTTP.request("GET", "http://httpbin.org/get", response_stream=io)
close(io)
```

Stream body through `open() do io`:
```julia
r = HTTP.open("GET", "http://httpbin.org/stream/10") do io
   while !eof(io)
       println(String(readavailable(io)))
   end
end

using HTTP.IOExtras

HTTP.open("GET", "https://tinyurl.com/bach-cello-suite-1-ogg") do http
    n = 0
    r = startread(http)
    l = parse(Int, HTTP.header(r, "Content-Length"))
    open(`vlc -q --play-and-exit --intf dummy -`, "w") do vlc
        while !eof(http)
            bytes = readavailable(http)
            write(vlc, bytes)
            n += length(bytes)
            println("streamed \$n-bytes \$((100*n)÷l)%\\u1b[1A")
        end
    end
end
```

## Request and Response Body Examples

String bodies:
```julia
r = HTTP.request("POST", "http://httpbin.org/post", [], "post body data")
println(String(r.body))
```

Interfacing with RESTful JSON APIs:
```julia
using JSON
params = Dict("user"=>"RAO...tjN", "token"=>"NzU...Wnp", "message"=>"Hello!")
base_url = "http://api.domain.com"
endpoint = "/1/messages.json"
url = base_url * endpoint
r = HTTP.request("POST", url,
             ["Content-Type" => "application/json"],
             JSON.json(params))
println(JSON.parse(String(r.body)))
```

Stream bodies from and to files:
```julia
in = open("foo.png", "r")
out = open("foo.jpg", "w")
HTTP.request("POST", "http://convert.com/png2jpg", [], in, response_stream=out)
```

Stream bodies through: `open() do io`:
```julia
using HTTP.IOExtras

HTTP.open("POST", "http://music.com/play") do io
    write(io, JSON.json([
        "auth" => "12345XXXX",
        "song_id" => 7,
    ]))
    r = startread(io)
    @show r.status
    while !eof(io)
        bytes = readavailable(io))
        play_audio(bytes)
    end
end
```
"""
function request(method, url, h=Header[], b=nobody;
                 headers=h, body=b, query=nothing, kw...)::Response
    return request(HTTP.stack(;kw...), string(method), request_uri(url, query), mkheaders(headers), body; kw...)
end
function request(stack::Type{<:Layer}, method, url, h=Header[], b=nobody;
                 headers=h, body=b, query=nothing, kw...)::Response
    return request(stack, string(method), request_uri(url, query), mkheaders(headers), body; kw...)
end

request(::Type{Union{}}, resp::Response) = resp

request_uri(url, query) = URI(URI(url); query=query)
request_uri(url, ::Nothing) = URI(url)

"""
    HTTP.open(method, url, [,headers]) do io
        write(io, body)
        [startread(io) -> HTTP.Response]
        while !eof(io)
            readavailable(io) -> AbstractVector{UInt8}
        end
    end -> HTTP.Response

The `HTTP.open` API allows the Request Body to be written to (and/or the
Response Body to be read from) an `IO` stream.

e.g. Streaming an audio file to the `vlc` player:
```julia
HTTP.open(:GET, "https://tinyurl.com/bach-cello-suite-1-ogg") do http
    open(`vlc -q --play-and-exit --intf dummy -`, "w") do vlc
        write(vlc, http)
    end
end
```
"""
open(f::Function, method::Union{String,Symbol}, url, headers=Header[]; kw...)::Response =
    request(string(method), url, headers, nothing; iofunction=f, kw...)

"""
    HTTP.openraw(method, url, [, headers])::Tuple{TCPSocket, Response, ByteView}

Open a raw socket that is unmanaged by HTTP.jl. Useful for doing HTTP upgrades
to other protocols.  Any bytes of the body read from the socket when reading
headers, is returned as excess bytes in the last tuple argument.

Example of a WebSocket upgrade:
```julia
headers = Dict(
    "Upgrade" => "websocket",
    "Connection" => "Upgrade",
    "Sec-WebSocket-Key" => "dGhlIHNhbXBsZSBub25jZQ==",
    "Sec-WebSocket-Version" => "13")

socket, response, excess = HTTP.openraw("GET", "ws://echo.websocket.org", headers)

# Write a WebSocket frame
frame = UInt8[0x81, 0x85, 0x37, 0xfa, 0x21, 0x3d, 0x7f, 0x9f, 0x4d, 0x51, 0x58]
write(socket, frame)
```
"""
function openraw(method::Union{String,Symbol}, url, headers=Header[]; kw...)::Tuple{IO, Response}
    socketready = Channel{Tuple{IO, Response}}(0)
    @async HTTP.open(method, url, headers; kw...) do http
        HTTP.startread(http)
        socket = http.stream
        put!(socketready, (socket, http.message))
        while(isopen(socket))
            Base.wait_close(socket)
        end
    end
    take!(socketready)
end

"""
    HTTP.get(url [, headers]; <keyword arguments>) -> HTTP.Response

Shorthand for `HTTP.request("GET", ...)`. See [`HTTP.request`](@ref).
"""
get(a...; kw...) = request("GET", a...; kw...)

"""
    HTTP.put(url, headers, body; <keyword arguments>) -> HTTP.Response

Shorthand for `HTTP.request("PUT", ...)`. See [`HTTP.request`](@ref).
"""
put(u, h=[], b=""; kw...) = request("PUT", u, h, b; kw...)

"""
    HTTP.post(url, headers, body; <keyword arguments>) -> HTTP.Response

Shorthand for `HTTP.request("POST", ...)`. See [`HTTP.request`](@ref).
"""
post(u, h=[], b=""; kw...) = request("POST", u, h, b; kw...)

"""
    HTTP.patch(url, headers, body; <keyword arguments>) -> HTTP.Response

Shorthand for `HTTP.request("PATCH", ...)`. See [`HTTP.request`](@ref).
"""
patch(u, h=[], b=""; kw...) = request("PATCH", u, h, b; kw...)

"""
    HTTP.head(url; <keyword arguments>) -> HTTP.Response

Shorthand for `HTTP.request("HEAD", ...)`. See [`HTTP.request`](@ref).
"""
head(u; kw...) = request("HEAD", u; kw...)

"""
    HTTP.delete(url [, headers]; <keyword arguments>) -> HTTP.Response

Shorthand for `HTTP.request("DELETE", ...)`. See [`HTTP.request`](@ref).
"""
delete(a...; kw...) = request("DELETE", a...; kw...)

include("RedirectRequest.jl");          using .RedirectRequest
include("BasicAuthRequest.jl");         using .BasicAuthRequest
include("AWS4AuthRequest.jl");          using .AWS4AuthRequest
include("CookieRequest.jl");            using .CookieRequest
include("CanonicalizeRequest.jl");      using .CanonicalizeRequest
include("TimeoutRequest.jl");           using .TimeoutRequest
include("MessageRequest.jl");           using .MessageRequest
include("ExceptionRequest.jl");         using .ExceptionRequest
                                        import .ExceptionRequest.StatusError
include("RetryRequest.jl");             using .RetryRequest
include("ConnectionRequest.jl");        using .ConnectionRequest
include("DebugRequest.jl");             using .DebugRequest
include("StreamRequest.jl");            using .StreamRequest
include("ContentTypeRequest.jl");       using .ContentTypeDetection
include("exceptions.jl")

"""
The `stack()` function returns the default HTTP Layer-stack type.
This type is passed as the first parameter to the [`HTTP.request`](@ref) function.

`stack()` accepts optional keyword arguments to enable/disable specific layers
in the stack:
`request(method, args...; kw...) request(stack(; kw...), args...; kw...)`

The minimal request execution stack is:

```julia
stack = MessageLayer{ConnectionPoolLayer{StreamLayer}}
```

The figure below illustrates the full request execution stack and its
relationship with [`HTTP.Response`](@ref), [`HTTP.Parsers`](@ref),
[`HTTP.Stream`](@ref) and the [`HTTP.ConnectionPool`](@ref).

```
 ┌────────────────────────────────────────────────────────────────────────────┐
 │                                            ┌───────────────────┐           │
 │  HTTP.jl Request Execution Stack           │ HTTP.ParsingError ├ ─ ─ ─ ─ ┐ │
 │                                            └───────────────────┘           │
 │                                            ┌───────────────────┐         │ │
 │                                            │ HTTP.IOError      ├ ─ ─ ─     │
 │                                            └───────────────────┘      │  │ │
 │                                            ┌───────────────────┐           │
 │                                            │ HTTP.StatusError  │─ ─   │  │ │
 │                                            └───────────────────┘   │       │
 │                                            ┌───────────────────┐      │  │ │
 │     request(method, url, headers, body) -> │ HTTP.Response     │   │       │
 │             ──────────────────────────     └─────────▲─────────┘      │  │ │
 │                           ║                          ║             │       │
 │   ┌────────────────────────────────────────────────────────────┐      │  │ │
 │   │ request(RedirectLayer,     method, ::URI, ::Headers, body) │   │       │
 │   ├────────────────────────────────────────────────────────────┤      │  │ │
 │   │ request(BasicAuthLayer,    method, ::URI, ::Headers, body) │   │       │
 │   ├────────────────────────────────────────────────────────────┤      │  │ │
 │   │ request(CookieLayer,       method, ::URI, ::Headers, body) │   │       │
 │   ├────────────────────────────────────────────────────────────┤      │  │ │
 │   │ request(CanonicalizeLayer, method, ::URI, ::Headers, body) │   │       │
 │   ├────────────────────────────────────────────────────────────┤      │  │ │
 │   │ request(MessageLayer,      method, ::URI, ::Headers, body) │   │       │
 │   ├────────────────────────────────────────────────────────────┤      │  │ │
 │   │ request(AWS4AuthLayer,             ::URI, ::Request, body) │   │       │
 │   ├────────────────────────────────────────────────────────────┤      │  │ │
 │   │ request(RetryLayer,                ::URI, ::Request, body) │   │       │
 │   ├────────────────────────────────────────────────────────────┤      │  │ │
 │   │ request(ExceptionLayer,            ::URI, ::Request, body) ├ ─ ┘       │
 │   ├────────────────────────────────────────────────────────────┤      │  │ │
┌┼───┤ request(ConnectionPoolLayer,       ::URI, ::Request, body) ├ ─ ─ ─     │
││   ├────────────────────────────────────────────────────────────┤         │ │
││   │ request(DebugLayer,                ::IO,  ::Request, body) │           │
││   ├────────────────────────────────────────────────────────────┤         │ │
││   │ request(TimeoutLayer,              ::IO,  ::Request, body) │           │
││   ├────────────────────────────────────────────────────────────┤         │ │
││   │ request(StreamLayer,               ::IO,  ::Request, body) │           │
││   └──────────────┬───────────────────┬─────────────────────────┘         │ │
│└──────────────────┼────────║──────────┼───────────────║─────────────────────┘
│                   │        ║          │               ║                   │
│┌──────────────────▼───────────────┐   │  ┌──────────────────────────────────┐
││ HTTP.Request                     │   │  │ HTTP.Response                  │ │
││                                  │   │  │                                  │
││ method::String                   ◀───┼──▶ status::Int                    │ │
││ target::String                   │   │  │ headers::Vector{Pair}            │
││ headers::Vector{Pair}            │   │  │ body::Vector{UInt8}            │ │
││ body::Vector{UInt8}              │   │  │                                  │
│└──────────────────▲───────────────┘   │  └───────────────▲────────────────┼─┘
│┌──────────────────┴────────║──────────▼───────────────║──┴──────────────────┐
││ HTTP.Stream <:IO          ║           ╔══════╗       ║                   │ │
││   ┌───────────────────────────┐       ║   ┌──▼─────────────────────────┐   │
││   │ startwrite(::Stream)      │       ║   │ startread(::Stream)        │ │ │
││   │ write(::Stream, body)     │       ║   │ read(::Stream) -> body     │   │
││   │ ...                       │       ║   │ ...                        │ │ │
││   │ closewrite(::Stream)      │       ║   │ closeread(::Stream)        │   │
││   └───────────────────────────┘       ║   └────────────────────────────┘ │ │
│└───────────────────────────║────────┬──║──────║───────║──┬──────────────────┘
│┌──────────────────────────────────┐ │  ║ ┌────▼───────║──▼────────────────┴─┐
││ HTTP.Messages                    │ │  ║ │ HTTP.Parsers                     │
││                                  │ │  ║ │                                  │
││ writestartline(::IO, ::Request)  │ │  ║ │ parse_status_line(bytes, ::Req') │
││ writeheaders(::IO, ::Request)    │ │  ║ │ parse_header_field(bytes, ::Req')│
│└──────────────────────────────────┘ │  ║ └──────────────────────────────────┘
│                            ║        │  ║
│┌───────────────────────────║────────┼──║────────────────────────────────────┐
└▶ HTTP.ConnectionPool       ║        │  ║                                    │
 │                     ┌──────────────▼────────┐ ┌───────────────────────┐    │
 │ getconnection() ->  │ HTTP.Transaction <:IO │ │ HTTP.Transaction <:IO │    │
 │                     └───────────────────────┘ └───────────────────────┘    │
 │                           ║    ╲│╱    ║                  ╲│╱               │
 │                           ║     │     ║                   │                │
 │                     ┌───────────▼───────────┐ ┌───────────▼───────────┐    │
 │              pool: [│ HTTP.Connection       │,│ HTTP.Connection       │...]│
 │                     └───────────┬───────────┘ └───────────┬───────────┘    │
 │                           ║     │     ║                   │                │
 │                     ┌───────────▼───────────┐ ┌───────────▼───────────┐    │
 │                     │ Base.TCPSocket <:IO   │ │MbedTLS.SSLContext <:IO│    │
 │                     └───────────────────────┘ └───────────┬───────────┘    │
 │                           ║           ║                   │                │
 │                           ║           ║       ┌───────────▼───────────┐    │
 │                           ║           ║       │ Base.TCPSocket <:IO   │    │
 │                           ║           ║       └───────────────────────┘    │
 └───────────────────────────║───────────║────────────────────────────────────┘
                             ║           ║
 ┌───────────────────────────║───────────║──────────────┐  ┏━━━━━━━━━━━━━━━━━━┓
 │ HTTP Server               ▼                          │  ┃ data flow: ════▶ ┃
 │                        Request     Response          │  ┃ reference: ────▶ ┃
 └──────────────────────────────────────────────────────┘  ┗━━━━━━━━━━━━━━━━━━┛
```
*See `docs/src/layers`[`.monopic`](http://monodraw.helftone.com).*
"""
function stack(;redirect=true,
                aws_authorization=false,
                cookies=false,
                canonicalize_headers=false,
                retry=true,
                status_exception=true,
                readtimeout=0,
                detect_content_type=false,
                verbose=0,
                kw...)

    NoLayer = Union

    stack = (redirect     ? RedirectLayer             : NoLayer){
                            BasicAuthLayer{
    (detect_content_type  ? ContentTypeDetectionLayer : NoLayer){
    (cookies === true || (cookies isa AbstractDict && !isempty(cookies)) ?
                            CookieLayer               : NoLayer){
    (canonicalize_headers ? CanonicalizeLayer         : NoLayer){
                            MessageLayer{
    (aws_authorization    ? AWS4AuthLayer             : NoLayer){
    (retry                ? RetryLayer                : NoLayer){
    (status_exception     ? ExceptionLayer            : NoLayer){
                            ConnectionPoolLayer{
    (verbose >= 3 ||
     DEBUG_LEVEL[] >= 3   ? DebugLayer                : NoLayer){
    (readtimeout > 0      ? TimeoutLayer              : NoLayer){
                            StreamLayer{Union{}}
    }}}}}}}}}}}}

    reduce(Layers.EXTRA_LAYERS; init=stack) do stack, (before, custom)
        insert(stack, before, custom)
    end
end

include("download.jl")
include("Servers.jl")                  ;using .Servers; using .Servers: listen
include("Handlers.jl")                 ;using .Handlers; using .Handlers: serve
include("parsemultipart.jl")
include("WebSockets.jl")               ;using .WebSockets

import .ConnectionPool: Transaction, Connection

function Base.parse(::Type{T}, str::AbstractString)::T where T <: Message
    buffer = Base.BufferStream()
    write(buffer, str)
    close(buffer)
    m = T()
    http = Stream(m, Transaction(Connection(buffer)))
    m.body = read(http)
    closeread(http)
    return m
end

end # module

ext: jl


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
using EntropyStableEuler
#using EntropyStableEuler.Fluxes1D

function wavespeed_1D(rho,rhou,E)
    p = pfun_nd(rho,(rhou,),E)
    cvel = @. sqrt(γ*p/rho)
    return @. abs(rhou/rho) + cvel
end
unorm(U) = sum(map((x->x.^2),U))
function pfun_nd(rho, rhoU, E)
    rhoUnorm2 = unorm(rhoU)./rho
    return @. (γ-1)*(E - .5*rhoUnorm2)
end

function primitive_to_conservative_hardcode(rho,U,p)
    rhoU = rho.*U
    Unorm = unorm(U)
    E = @. p/(γ-1) + .5*rho*Unorm
    return (rho,rhoU,E)
end

function euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
    rhoR,uR,betaR,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))

    unorm = (@. uL*uR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. f4aux*uavg)

    return (FxS1,FxS2,FxS3)
end

function euler_fluxes(rhoL,uL,betaL,rhoR,uR,betaR)
    rhologL,betalogL,rhologR,betalogR = map(x->log.(x),(rhoL,betaL,rhoR,betaR))
    return euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
                           rhoR,uR,betaR,rhologR,betalogR)
end

const TOL = 1e-16
"Approximation parameters"
N = 8 # The order of approximation
K = 10
T = 1.0
T = 6.0
#T = 0.0039

# # Sod shocktube
# const γ = 1.4
# const Bl = -0.5
# const Br = 0.5
# const rhoL = 1.0
# const rhoR = 0.125
# const pL = 1.0
# const pR = 0.1
# const xC = 0.0
# const GIFINTERVAL = 20
# T = 0.2

# Leblanc shocktube
const γ = 5/3
const Bl = 0.0
const Br = 9.0
const rhoL = 1.0
const rhoR = 0.001
const pL = 0.1
const pR = 1e-7
#const pR = 1e-15
const xC = 3.0
const GIFINTERVAL = 60
T = 6.0

"Mesh related variables"
VX = LinRange(Bl,Br,K+1)
EToV = transpose(reshape(sort([1:K; 2:K+1]),2,K))

"Initialize reference element"
r,_ = gauss_lobatto_quad(0,0,N)   # Reference nodes
VDM = vandermonde_1D(N,r)         # modal to nodal
Dr = grad_vandermonde_1D(N,r)/VDM # nodal differentiation
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1]) # nodal linear interpolation

Nq = N
rq,wq = gauss_lobatto_quad(0,0,Nq)
Vq = vandermonde_1D(N,rq)/VDM
M = Vq'*diagm(wq)*Vq
Mlump = zeros(size(M))
Mlump_inv = zeros(size(M))
for i = 1:Nq+1
    Mlump[i,i] = sum(M[i,:])
    Mlump_inv[i,i] = 1.0/Mlump[i,i]
end

# operators
Qr = M*Dr
B = zeros(N+1,N+1)
B[1,1] = -1
B[end,end] = 1
L = Array(spdiagm(0=>-2*ones(N+1), 1=>ones(N), -1=>ones(N)))
L[1,1] = -1
L[end,end] = -1
psi = pinv(L)*-1/2*B*ones(N+1)
S0 = zeros(N+1,N+1)
for i = 1:N+1
    for j = 1:N+1
        if L[i,j] != 0
            S0[i,j] = psi[j] - psi[i]
        end
    end
end
Qr0 = S0+1/2*B

# Drop zeros
Qr = Matrix(droptol!(sparse(Qr),TOL))
Qr0 = Matrix(droptol!(sparse(Qr0),TOL))
B = Matrix(droptol!(sparse(B),TOL))
L = Matrix(droptol!(sparse(L),TOL))
S0 = Matrix(droptol!(sparse(S0),TOL))
S = Matrix(droptol!(sparse((Qr-Qr')/2),TOL))

rf = [-1.0;1.0]
nrJ = [-1.0;1.0]
Vf = vandermonde_1D(N,rf)/VDM

"""High order mesh"""
x = V1*VX[transpose(EToV)]
xf = Vf*x
mapM = reshape(1:2*K,2,K)
mapP = copy(mapM)
mapP[1,2:end] .= mapM[2,1:end-1]
mapP[2,1:end-1] .= mapM[1,2:end]

# """Periodic"""
# mapP[1] = mapM[end]
# mapP[end] = mapP[1]

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

# """Low order mesh"""
# VX0 = [x[1:end-1,:][:]; x[end]]
# #EToV = transpose(reshape(sort([1:K; 2:K+1]),2,K)
# EToV0 = transpose(reshape(sort([1:N*K; 2:N*K+1]),2,N*K))
# x0 = VX0[transpose(EToV0)]
# xf0 = x0
# mapM0 = reshape(1:2*N*K,2,N*K)
# mapP0 = copy(mapM0)
# mapP0[1,2:end] .= mapM0[2,1:end-1]
# mapP0[2,1:end-1] .= mapM0[1,2:end]

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

"""Initial condition"""
rho_x(x) = (x <= xC) ? rhoL : rhoR
u_x(x) = 0.0
p_x(x) = (x <= xC) ? pL : pR

rho = @. rho_x(x)
u = @. u_x(x)
p = @. p_x(x)
U = primitive_to_conservative_hardcode(rho,u,p)

function limiting_param(U_low, P_ij)
    l = 1.0
    # Limit density
    if U_low[1] + P_ij[1] < -TOL
        l = min(abs(U_low[1])/(abs(P_ij[1])+1e-14), 1.0)
    end

    # limiting internal energy (via quadratic function)
    a = P_ij[1]*P_ij[3]-1.0/2.0*P_ij[2]^2
    b = U_low[3]*P_ij[1]+U_low[1]*P_ij[3]-U_low[2]*P_ij[2]
    c = U_low[3]*U_low[1]-1.0/2.0*U_low[2]^2

    l_eps_ij = 1.0
    if b^2-4*a*c >= 0
        r1 = (-b+sqrt(b^2-4*a*c))/(2*a)
        r2 = (-b-sqrt(b^2-4*a*c))/(2*a)
        if r1 > TOL && r2 > TOL
            l_eps_ij = min(r1,r2)
        elseif r1 > TOL && r2 < -TOL
            l_eps_ij = r1
        elseif r2 > TOL && r1 < -TOL
            l_eps_ij = r2
        end
    end

    l = min(l,l_eps_ij)
    return l
end

function construct_artificial_wavespd(rhoL,mL,EL,rhoR,mR,ER,F1,F2,F3)
    dij = abs(F1/(rhoL+rhoR))
    rhosum = rhoL+rhoR
    msum = mL+mR
    Esum = EL+ER
    a = rhosum*Esum-1/2*msum^2
    b = rhosum*F3+Esum*F1-msum*F2
    c = F1*F3-1/2*F2^2

    if c >= 0
        return max(dij,0.0)
    else
        return max(dij,(b+sqrt(b^2-4*a*c))/(2*a))
    end
end

function flux_lowIDP(U_i,U_j,f_i,f_j,c_ij,wavespd)
    return c_ij*(f_i+f_j)-abs(c_ij)*wavespd*(U_j-U_i)
end

function flux_high(f_i,f_j,c_ij)
    return c_ij*(f_i+f_j)
end

function flux_ES(rho_i,u_i,beta_i,rho_j,u_j,beta_j,c_ij)
    return 2*c_ij.*euler_fluxes(rho_i,u_i,beta_i,rho_j,u_j,beta_j)
end

function rhs_IDP(U,K,N,wq,S,S0,Mlump_inv)
    p = pfun_nd.(U[1],U[2],U[3])
    flux = zero.(U)
    @. flux[1] = U[2]
    @. flux[2] = U[2]^2/U[1]+p
    @. flux[3] = U[3]*U[2]/U[1]+p*U[2]/U[1]

    Ub = zero.(U)
    @. Ub[1] = U[1]
    @. Ub[2] = U[2]/U[1]
    @. Ub[3] = U[1]/(2*p)

    J = (Br-Bl)/K/2 # assume uniform interval

    # Low order and high order algebraic fluxes
    F_low       = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_high      = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_low_P     = [zeros(2),zeros(2),zeros(2)] # 1: left boundary, 2: right boundary
    F_high_P    = [zeros(2),zeros(2),zeros(2)]
    #L           = zeros(N+1,N+1) # Array of limiting params
    L           = zeros(N+1)
    rhsU        = [zeros(N+1,K),zeros(N+1,K),zeros(N+1,K)]

    L_plot = zeros(K)
    L_plot2 = zeros(N+1,K)

    wavespd_arr = zeros(N+1,K)
    for k = 1:K
        for i = 1:N+1
            wavespd_arr[i,k] = wavespeed_1D(U[1][i,k],U[2][i,k],U[3][i,k])
        end
    end

    d_ii_arr = zeros(N+1,K)
    for k = 1:K
        for i = 1:N+1
            for j = 1:N+1
                if i != j
                    wavespd = max(wavespd_arr[i,k],wavespd_arr[j,k])
                    d_ij = wavespd*abs(S0[i,j])
                    d_ii_arr[i,k] -= d_ij
                end
            end
            if i == 1
                U_left   = (k == 1) ? [rhoL; 0.0; pL/(γ-1)]    : [U[1][end,k-1]; U[2][end,k-1]; U[3][end,k-1]]
                wavespd = max(wavespd_arr[1,k],wavespeed_1D(U_left[1],U_left[2],U_left[3]))
                d_ij = wavespd/2.0
                d_ii_arr[i,k] -= d_ij
            end
            if i == N+1
                U_right  = (k == K) ? [rhoR; 0.0; pR/(γ-1)]    : [U[1][1,k+1]; U[2][1,k+1]; U[3][1,k+1]]
                wavespd = max(wavespd_arr[end,k],wavespeed_1D(U_right[1],U_right[2],U_right[3]))
                d_ij = wavespd/2.0
                d_ii_arr[i,k] -= d_ij
            end
        end
    end

    dt = minimum(-J/2*M*(1 ./d_ii_arr))

    for k = 1:K
        # Assemble matrix of low and high order algebraic fluxes
        # interior of the element
        for i = 1:N+1
            for j = 1:N+1
                if i != j # skip diagonal
                    wavespd = max(wavespd_arr[i,k],wavespd_arr[j,k])
                    fluxS = flux_ES(Ub[1][i,k],Ub[2][i,k],Ub[3][i,k],Ub[1][j,k],Ub[2][j,k],Ub[3][j,k],S[i,j])
                    for c = 1:3
                        F_low[c][i,j]  = flux_lowIDP(U[c][i,k],U[c][j,k],flux[c][i,k],flux[c][j,k],S0[i,j],wavespd)
                        F_high[c][i,j] = fluxS[c]
                    end
                end
            end
        end

        # Assemble matrix of low and high order algebraic fluxes
        # interface of the element

        U_left   = (k == 1) ? [rhoL; 0.0; pL/(γ-1)]    : [U[1][end,k-1]; U[2][end,k-1]; U[3][end,k-1]]
        Ub_left  = (k == 1) ? [rhoL; 0.0; rhoL/(2*pL)] : [Ub[1][end,k-1]; Ub[2][end,k-1]; Ub[3][end,k-1]]
        f_left   = (k == 1) ? [0.0; pL; 0.0]           : [flux[1][end,k-1]; flux[2][end,k-1]; flux[3][end,k-1]]
        U_right  = (k == K) ? [rhoR; 0.0; pR/(γ-1)]    : [U[1][1,k+1]; U[2][1,k+1]; U[3][1,k+1]]
        Ub_right = (k == K) ? [rhoR; 0.0; rhoR/(2*pR)] : [Ub[1][1,k+1]; Ub[2][1,k+1]; Ub[3][1,k+1]]
        f_right  = (k == K) ? [0.0; pR; 0.0]           : [flux[1][1,k+1]; flux[2][1,k+1]; flux[3][1,k+1]]
        wavespd_l = max(wavespd_arr[1,k],wavespeed_1D(U_left[1],U_left[2],U_left[3]))
        wavespd_r = max(wavespd_arr[end,k],wavespeed_1D(U_right[1],U_right[2],U_right[3]))

        for c = 1:3
            F_low_P[c][1] = flux_lowIDP(U[c][1,k],U_left[c],flux[c][1,k],f_left[c],-0.5,wavespd_l)
            F_low_P[c][2] = flux_lowIDP(U[c][end,k],U_right[c],flux[c][end,k],f_right[c],0.5,wavespd_r)

            F_high_P[c][1] = F_low_P[c][1]
            F_high_P[c][2] = F_low_P[c][2]
        end

        # Calculate limiting parameters over interior of the element
        # TODO: redundant
        U_low =  [zeros(N+1,1),zeros(N+1,1),zeros(N+1,1)]
        for c = 1:3
            U_low[c] .= sum(-F_low[c],dims=2)
            U_low[c][1] -= F_low_P[c][1]
            U_low[c][end] -= F_low_P[c][2]
            U_low[c] .= U[c][:,k]+dt/J*Mlump_inv*U_low[c]
        end

        for i = 1:N+1
            m_i = J*wq[i]
            P_i = [0.0;0.0;0.0]
            for j = 1:N+1
                for c = 1:3
                    P_i[c] += dt/m_i*(F_low[c][i,j]-F_high[c][i,j])
                end
            end
            L[i] = limiting_param([U_low[1][i]; U_low[2][i]; U_low[3][i]],P_i)
        end
        l = minimum(L)
        #@show l
        if l < 1
            l = l - 1e-10
        end

        # construct rhs
        for c = 1:3
            # With limiting
            rhsU[c][:,k] = sum((l-1).*F_low[c] - l*F_high[c],dims=2)
            # for i = 1:N+1
            #     rhsU[c][i,k] = sum((L[i]-1)*F_low[c][i,:] - L[i]*F_high[c][i,:])
            # end
            rhsU[c][1,k] += -F_low_P[c][1]
            rhsU[c][N+1,k] += -F_low_P[c][2]

            rhsU[c][:,k] .= 1/J*Mlump_inv*rhsU[c][:,k]
        end

        L_plot[k] = sum(L)
        for i = 1:N+1
            L_plot2[i,k] = L[i]
        end
    end

    return rhsU,L_plot,L_plot2,dt
end

function rhs_IDPlow(U,K,N,Mlump_inv,p,flux,J)
    dfdx = (zeros(N+1,K),zeros(N+1,K),zeros(N+1,K))
    for i = 2:K*(N+1)-1
        for c = 1:3
            dfdx[c][i] = 1/2*(flux[c][mod1(i+1,K*(N+1))] - flux[c][mod1(i-1,K*(N+1))])
        end
    end
    dfdx[1][1] = 1/2*(flux[1][2] - 0.0)
    dfdx[2][1] = 1/2*(flux[2][2] - pL)
    dfdx[3][1] = 1/2*(flux[3][2] - 0.0)
    dfdx[1][end] = 1/2*(0.0 - flux[1][end-1])
    dfdx[2][end] = 1/2*(pR - flux[2][end-1])
    dfdx[3][end] = 1/2*(0.0 - flux[3][end-1])

    visc = (zeros(N+1,K),zeros(N+1,K),zeros(N+1,K))
    for i = 2:K*(N+1)-1
        wavespd_curr = wavespeed_1D(U[1][i],U[2][i],U[3][i])
        wavespd_R = wavespeed_1D(U[1][mod1(i+1,K*(N+1))],U[2][mod1(i+1,K*(N+1))],U[3][mod1(i+1,K*(N+1))])
        wavespd_L = wavespeed_1D(U[1][mod1(i-1,K*(N+1))],U[2][mod1(i-1,K*(N+1))],U[3][mod1(i-1,K*(N+1))])
        dL = 1/2*max(wavespd_curr,wavespd_L)
        dR = 1/2*max(wavespd_curr,wavespd_R)
        for c = 1:3
            visc[c][i] = dL*(U[c][mod1(i-1,K*(N+1))]-U[c][i]) + dR*(U[c][mod1(i+1,K*(N+1))]-U[c][i])
        end
    end

    # i = 1
    wavespd_curr = wavespeed_1D(U[1][1],U[2][1],U[3][1])
    wavespd_R = wavespeed_1D(U[1][2],U[2][2],U[3][2])
    wavespd_L = wavespeed_1D(rhoL,0.0,pL/(γ-1))
    dL = 1/2*max(wavespd_curr,wavespd_L)
    dR = 1/2*max(wavespd_curr,wavespd_R)
    visc[1][1] = dL*(rhoL-U[1][1]) + dR*(U[1][2]-U[1][1])
    visc[2][1] = dL*(0.0-U[2][1]) + dR*(U[2][2]-U[2][1])
    visc[3][1] = dL*(pL/(γ-1)-U[3][1]) + dR*(U[3][2]-U[3][1])

    # i = end
    wavespd_curr = wavespeed_1D(U[1][end],U[2][end],U[3][end])
    wavespd_R = wavespeed_1D(rhoR,0.0,pR/(γ-1))
    wavespd_L = wavespeed_1D(U[1][end-1],U[2][end-1],U[3][end-1])
    dL = 1/2*max(wavespd_curr,wavespd_L)
    dR = 1/2*max(wavespd_curr,wavespd_R)
    visc[1][end] = dL*(U[1][end-1]-U[1][end]) + dR*(rhoR-U[1][end])
    visc[2][end] = dL*(U[2][end-1]-U[2][end]) + dR*(0.0-U[2][end])
    visc[3][end] = dL*(U[3][end-1]-U[3][end]) + dR*(pR/(γ-1)-U[3][end])

    rhsU = (x->1/J*Mlump_inv*x).(.-dfdx.+visc)
    return rhsU
end

function rhs_high(U,K,N,Mlump_inv,S)
    p = pfun_nd.(U[1],U[2],U[3])
    flux = zero.(U)
    @. flux[1] = U[2]
    @. flux[2] = U[2]^2/U[1]+p
    @. flux[3] = U[3]*U[2]/U[1]+p*U[2]/U[1]

    Ub = zero.(U)
    @. Ub[1] = U[1]
    @. Ub[2] = U[2]/U[1]
    @. Ub[3] = U[1]/(2*p)

    J = (Br-Bl)/K/2 # assume uniform interval

    # Low order and high order algebraic fluxes
    F_high      = [zeros(N+1,N+1),zeros(N+1,N+1),zeros(N+1,N+1)]
    F_high_P    = [zeros(2),zeros(2),zeros(2)]
    rhsU        = [zeros(N+1,K),zeros(N+1,K),zeros(N+1,K)]

    # Determine dt
    wavespd_arr = zeros(N+1,K)
    for k = 1:K
        for i = 1:N+1
            wavespd_arr[i,k] = wavespeed_1D(U[1][i,k],U[2][i,k],U[3][i,k])
        end
    end

    for k = 1:K
        # Assemble matrix of low and high order algebraic fluxes
        # interior of the element
        for i = 1:N+1
            for j = 1:N+1
                if i != j # skip diagonal
                    # wavespd = max(wavespd_arr[i,k],wavespd_arr[j,k])
                    # for c = 1:3
                    #     F_high[c][i,j] = flux_ES(Ub[c][i,k],Ub[c][j,k],S[i,j])
                    # end
                    tmp = flux_ES(Ub[1][i,k],Ub[2][i,k],Ub[3][i,k],Ub[1][j,k],Ub[2][j,k],Ub[3][j,k],S[i,j])
                    for c = 1:3
                        F_high[c][i,j] = tmp[c]
                    end
                end
            end
        end

        # Assemble matrix of low and high order algebraic fluxes
        # interface of the element
        U_left   = (k == 1) ? [rhoL; 0.0; pL/(γ-1)]    : [U[1][end,k-1]; U[2][end,k-1]; U[3][end,k-1]]
        Ub_left  = (k == 1) ? [rhoL; 0.0; rhoL/(2*pL)] : [Ub[1][end,k-1]; Ub[2][end,k-1]; Ub[3][end,k-1]]
        f_left   = (k == 1) ? [0.0; pL; 0.0]           : [flux[1][end,k-1]; flux[2][end,k-1]; flux[3][end,k-1]]
        U_right  = (k == K) ? [rhoR; 0.0; pR/(γ-1)]    : [U[1][1,k+1]; U[2][1,k+1]; U[3][1,k+1]]
        Ub_right = (k == K) ? [rhoR; 0.0; rhoR/(2*pR)] : [Ub[1][1,k+1]; Ub[2][1,k+1]; Ub[3][1,k+1]]
        f_right  = (k == K) ? [0.0; pR; 0.0]           : [flux[1][1,k+1]; flux[2][1,k+1]; flux[3][1,k+1]]
        wavespd_l = max(wavespd_arr[1,k],wavespeed_1D(U_left[1],U_left[2],U_left[3]))
        wavespd_r = max(wavespd_arr[end,k],wavespeed_1D(U_right[1],U_right[2],U_right[3]))

        tmp = .-euler_fluxes(Ub[1][1,k],Ub[2][1,k],Ub[3][1,k],Ub_left[1],Ub_left[2],Ub_left[3])
        tmp2 = euler_fluxes(Ub[1][end,k],Ub[2][end,k],Ub[3][end,k],Ub_right[1],Ub_right[2],Ub_right[3])
        for c = 1:3
            # F_high_P[c][1] = -1/2*(flux[c][1,k]+f_left[c])-wavespd_l/2*(U_left[c]-U[c][1,k])#tmp[c]-wavespd_l/2*(U_left[c]-U[c][1,k])#-euler_fluxes(Ub[c][1,k],Ub_left[c])#flux_high(Ub[c][1,k],Ub_left[c],-0.5)-wavespd_l/2*(Ub_left[c]-Ub[c][1,k])
            # F_high_P[c][2] = 1/2*(flux[c][end,k]+f_right[c])-wavespd_r/2*(U_right[c]-U[c][end,k])#tmp2[c]-wavespd_r/2*(U_right[c]-U[c][end,k])#euler_fluxes(Ub[c][end,k],Ub_right[c])#flux_high(Ub[c][end,k],Ub_right[c],0.5)-wavespd_r/2*(Ub_right[c]-Ub[c][end,k])
            F_high_P[c][1] = tmp[c]-wavespd_l/2*(U_left[c]-U[c][1,k])#-euler_fluxes(Ub[c][1,k],Ub_left[c])#flux_high(Ub[c][1,k],Ub_left[c],-0.5)-wavespd_l/2*(Ub_left[c]-Ub[c][1,k])
            F_high_P[c][2] = tmp2[c]-wavespd_r/2*(U_right[c]-U[c][end,k])#euler_fluxes(Ub[c][end,k],Ub_right[c])#flux_high(Ub[c][end,k],Ub_right[c],0.5)-wavespd_r/2*(Ub_right[c]-Ub[c][end,k])
        end

        # construct rhs
        for c = 1:3
            # With limiting
            rhsU[c][:,k] = -sum(F_high[c],dims=2)
            rhsU[c][1,k] += -F_high_P[c][1]
            rhsU[c][N+1,k] += -F_high_P[c][2]
            rhsU[c][:,k] .= 1/J*Mlump_inv*rhsU[c][:,k]
        end

    end

    return rhsU
end

# Time stepping
"Time integration"
t = 0.0
U = collect(U)
resU = [zeros(size(x)),zeros(size(x)),zeros(size(x))]

Vp = vandermonde_1D(N,LinRange(-1,1,10))/VDM
gr(size=(300,300),ylims=(0,1.2),legend=false,markerstrokewidth=1,markersize=2)
plot()
dt_hist = []
anim = Animation()
i = 1

while t < T
    rhsU,L_plot,L_plot2,dt = rhs_IDP(U,K,N,wq,S,S0,Mlump_inv)
    @. U = U + dt*rhsU
    push!(dt_hist,dt)
    global t = t + dt
    global i = i + 1
    println("Current time $t with time step size $dt, and final time $T")
    if mod(i,10) == 1
        plot(Vp*x,Vp*U[1])
        ptL = Bl+(Br-Bl)/K/(N+1)/2
        ptR = Br-(Br-Bl)/K/(N+1)/2
        hplot = (Br-Bl)/K/(N+1)
        for k = 1:K
            plot!(ptL+(k-1)*hplot*(N+1):hplot:ptL+k*hplot*(N+1), 1 .-L_plot2[:,k],st=:bar,alpha=0.2)
        end
        frame(anim)
    end
end

gif(anim,"~/Desktop/tmp.gif",fps=15)

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
module MyExample

export NSKRRegressor, KKRegressor
export fit,predict,predict_LOO,fit_iteratively

export SparseTensor
export readtensor

export Kfold, CVestimate, CVscore, optimizeHyperParameters

export auc_,R2_

export SE, self, L2

#include("NSKRRegressor.jl")
include("ModelSelection.jl")
#include("")
#include("MatrixProductState.jl")

#moeten uiteindelijk niet echt geexporteerd worden:
export Mps
export optimalmul

end # module

max_stars_repo_name: padwulf/MyExample

max_issues_repo_name: padwulf/MyExample

max_forks_repo_name: padwulf/MyExample


text:
module Tokens

import Base.eof

export Token

include("token_kinds.jl")

iskeyword(k::Kind) = begin_keywords < k < end_keywords
isliteral(k::Kind) = begin_literal < k < end_literal
isoperator(k::Kind) = begin_ops < k < end_ops

iscontextualkeyword(k::Kind) = begin_contextual_keywords < k < end_contextual_keywords

function iswordoperator(k::Kind)
    # Keyword-like operators
    k == Tokens.IN ||
    k == Tokens.ISA ||
    k == Tokens.WHERE
end

# Create string => keyword kind
const KEYWORDS = Dict{String, Kind}()

function _add_kws()
    for k in instances(Kind)
        if iskeyword(k)
            KEYWORDS[lowercase(string(k))] = k
        end
    end
end
_add_kws()

# TODO: more
@enum(TokenError,
    NO_ERR,
    EOF_MULTICOMMENT,
    EOF_CHAR,
    INVALID_NUMERIC_CONSTANT,
    INVALID_OPERATOR,
    INVALID_INTERPOLATION_TERMINATOR,
    UNKNOWN,
)

# Error kind => description
TOKEN_ERROR_DESCRIPTION = Dict{TokenError, String}(
    EOF_MULTICOMMENT => "unterminated multi-line comment #= ... =#",
    EOF_CHAR => "unterminated character literal",
    INVALID_NUMERIC_CONSTANT => "invalid numeric constant",
    INVALID_OPERATOR => "invalid operator",
    INVALID_INTERPOLATION_TERMINATOR => "interpolated variable ends with invalid character; use `\$(...)` instead",
    UNKNOWN => "unknown",
)

abstract type AbstractToken end

struct Token <: AbstractToken
    kind::Kind
    # Offsets into a string or buffer
    startpos::Tuple{Int, Int} # row, col where token starts /end, col is a string index
    endpos::Tuple{Int, Int}
    startbyte::Int # The byte where the token start in the buffer
    endbyte::Int # The byte where the token ended in the buffer
    val::String # The actual string of the token
    token_error::TokenError
    dotop::Bool
    suffix::Bool
end
function Token(kind::Kind, startposition::Tuple{Int, Int}, endposition::Tuple{Int, Int},
    startbyte::Int, endbyte::Int, val::String)
Token(kind, startposition, endposition, startbyte, endbyte, val, NO_ERR, false, false)
end
Token() = Token(ERROR, (0,0), (0,0), 0, 0, "", UNKNOWN, false, false)

struct RawToken <: AbstractToken
    kind::Kind
    # Offsets into a string or buffer
    startpos::Tuple{Int, Int} # row, col where token starts /end, col is a string index
    endpos::Tuple{Int, Int}
    startbyte::Int # The byte where the token start in the buffer
    endbyte::Int # The byte where the token ended in the buffer
    token_error::TokenError
    dotop::Bool
    suffix::Bool
end
function RawToken(kind::Kind, startposition::Tuple{Int, Int}, endposition::Tuple{Int, Int},
    startbyte::Int, endbyte::Int)
RawToken(kind, startposition, endposition, startbyte, endbyte, NO_ERR, false, false)
end
RawToken() = RawToken(ERROR, (0,0), (0,0), 0, 0, UNKNOWN, false, false)

const _EMPTY_TOKEN = Token()
const _EMPTY_RAWTOKEN = RawToken()
EMPTY_TOKEN(::Type{Token}) = _EMPTY_TOKEN
EMPTY_TOKEN(::Type{RawToken}) = _EMPTY_RAWTOKEN

function kind(t::AbstractToken)
    isoperator(t.kind) && return OP
    iskeyword(t.kind) && return KEYWORD
    return t.kind
end
exactkind(t::AbstractToken) = t.kind
startpos(t::AbstractToken) = t.startpos
endpos(t::AbstractToken) = t.endpos
startbyte(t::AbstractToken) = t.startbyte
endbyte(t::AbstractToken) = t.endbyte
function untokenize(t::Token)
    if t.kind == IDENTIFIER || isliteral(t.kind) || t.kind == COMMENT || t.kind == WHITESPACE || t.kind == NEWLINE_WS || t.kind == ERROR
        return t.val
    elseif iskeyword(t.kind)
        return lowercase(string(t.kind))
    elseif isoperator(t.kind)
        if t.dotop
            str = string(".", UNICODE_OPS_REVERSE[t.kind])
        else
            str = string(UNICODE_OPS_REVERSE[t.kind])
        end
        return string(str, t.val)
    elseif t.kind == LPAREN
        return "("
    elseif t.kind == LSQUARE
        return "["
    elseif t.kind == LBRACE
        return "{"
    elseif t.kind == RPAREN
        return ")"
    elseif t.kind == RSQUARE
        return "]"
    elseif t.kind == RBRACE
        return "}"
    elseif t.kind == DQUOTE
        return "\""
    elseif t.kind == TRIPLE_DQUOTE
        return "\"\"\""
    elseif t.kind == BACKTICK
        return "`"
    elseif t.kind == TRIPLE_BACKTICK
        return "```"
    elseif t.kind == AT_SIGN
        return "@"
    elseif t.kind == COMMA
        return ","
    elseif t.kind == SEMICOLON
        return ";"
    else
        return ""
    end
end

function untokenize(t::RawToken, str::String)
    String(codeunits(str)[1 .+ (t.startbyte:t.endbyte)])
end

function untokenize(ts)
    if !(eltype(ts) <: AbstractToken)
        throw(ArgumentError("element type of iterator has to be Token"))
    end
    io = IOBuffer()
    for tok in ts
        write(io, untokenize(tok))
    end
    return String(take!(io))
end

function Base.show(io::IO, t::Token)
    start_r, start_c = startpos(t)
    end_r, end_c = endpos(t)
    str = kind(t) == ENDMARKER ? "" : escape_string(untokenize(t))
    print(io, rpad(string(start_r, ",", start_c, "-", end_r, ",", end_c), 17, " "))
    print(io, rpad(kind(t), 15, " "))
    print(io, "\"", str, "\"")
end

Base.print(io::IO, t::Token) = print(io, untokenize(t))

function Base.show(io::IO, t::RawToken)
    start_r, start_c = startpos(t)
    end_r, end_c = endpos(t)
    print(io, rpad(string(start_r, ",", start_c, "-", end_r, ",", end_c), 17, " "))
    print(io, rpad(kind(t), 15, " "))
end

end # module


text:
type NPLCacheVars{T}
    m::Vector{Vector{T}}   # index
    eᴾ::Vector{Vector{T}}  # expected value of shock given optimality
    P::Vector{Vector{T}}   # for use in policy step
    Pᵈ::Vector{Vector{T}}  # for use in ML step
    z̃::Vector{Matrix{T}}   # ML
    z̃ᵈ::Vector{Matrix{T}}  # used in gradient!
    Δz̃::Vector{Matrix{T}}  # policy
    Δz̃ᵈ::Vector{Matrix{T}} # ML
    ẽ::Vector{Vector{T}}   # unsure if needed
    Δẽ::Vector{Vector{T}}  # policy
    Δẽᵈ::Vector{Vector{T}} # ML
    Fᵁ::AbstractMatrix{T}  # Controlled transition probability matrix
    Wᴾ::AbstractMatrix{T}  # Valuated Zs and eᵖs
    Θᵏ::Vector{T}          # Parameter cache
end
function NPLCacheVars(nX, nXᵈ, nvar, J, F)
    m   = [zeros(nX) for j = 1:J-1] # index
    eᴾ  = [zeros(nX) for j = 1:J]   # expected value of shock given optimality
    P   = [zeros(nX) for j = 1:J]   # for use in policy step
    Pᵈ  = [zeros(nXᵈ)  for j = 1:J] # for use in ML step
    z̃   = [zeros(nX,  nvar)  for j = 1:J]   # ML
    z̃ᵈ  = [zeros(nXᵈ, nvar) for j = 1:J] # used in gradient!
    Δz̃  = [zeros(nX,  nvar) for j = 1:J-1]  # policy
    Δz̃ᵈ = [zeros(nXᵈ, nvar) for j = 1:J-1] # ML
    ẽ   = [zeros(nX) for i in 1:J] # unsure if needed
    Δẽ  = [zeros(nX) for j = 1:J-1]  # policy
    Δẽᵈ = [zeros(nXᵈ) for j = 1:J-1] # ML
    Fᵁ  = zeros(F[1])         # Controlled transition probability matrix
    Wᴾ  = zeros(nX, nvar+1) # Valuated Zs and eᵖs
    Θᵏ  = zeros(nvar)
    NPLCacheVars(m, eᴾ, P, Pᵈ,
                 z̃, z̃ᵈ, Δz̃, Δz̃ᵈ,
                 ẽ, Δẽ, Δẽᵈ,
                 Fᵁ, Wᴾ, Θᵏ)
end

type ConvergenceNPL <: ConvergenceInfo
    maxK::Int64
    # Outer loop
    flag::Bool
    outer::Int64
    norm_Θ::Float64
    # Inner loop
    ll::Float64
    norm_grad::Float64
    iter_maxlike::Int64
end
ConvergenceNPL() = ConvergenceNPL(0, false, 0, 0., 0., 0., 0)

inner_iterations(conv::ConvergenceNPL) = conv.iter_maxlike

type NPL <: EstimationMethod
  U::LinearUtility
  Θu::Matrix{Float64}
  Θy::Vector{Float64}
  Θf::Vector{Float64}
  ll::Vector{Float64}
  first_step::ConvergenceNPL
  K::Integer
  finite::Bool
  verbose::Bool
end

function NPL(U, initial; maxK::Integer = 0, verbose::Bool = false)
    K = maxK == 0 ? 80 : maxK
    NPL(U,
        #Θu
        zeros(K, size(U.Z[1],2)),
        #Θy
        Vector{Float64}(),
        #Θf
        Vector{Float64}(),
        # ll
        zeros(K),
        # first step results (should just be a P)
        initial,
        # If maxK is set, it will run exactly maxK times.
        K,
        # Is it NPL-K? (as opposed to NPL-∞)
        K > 0,
        # Print information
        verbose)
end

function Base.display{M<:NPL}(est_res::EstimationResults{M})
    @printf "Results of estimation\n"
    @printf " * Method: NPL\n"
    @printf " * loglikelihood: %s\n" est_res.loglikelihood
    if length(join(est_res.trace.Θ[end,:], ",")) < 40
       @printf " * Estimate: [%s]\n" join(est_res.trace.Θ[end,:], ",")
       @printf " * Std.err.: [%s]\n" join(stderr(est_res), ",")
   else
       @printf " * Estimate: [%s, ...]\n" join(est_res.trace.Θ[end,1:2], ",")
       @printf " * Std.err.: [%s, ...]\n" join(stderr(est_res)[1:2], ",")
   end
   @printf " Iterations\n"
   n_PI = size(est_res.trace.Θ, 1)
   if n_PI == 1
       @printf " * Policy iteration: %d\n" size(est_res.trace.Θ, 1)
   else
       @printf " * Policy iteration: %d\n" size(est_res.trace.Θ, 1)
   end
   @printf " * Maximum likelihood: %d\n" est_res.conv.iter_maxlike
end

type TraceNPL
    Θ::Matrix{Float64}
    norm_ΔΘ::Vector{Float64}
    ll::Vector{Float64}
    Δll::Vector{Float64}
    norm_g::Vector{Float64}
end

TraceNPL(K, n) = TraceNPL(fill(NaN, K, n),
                          fill(NaN, K),
                          fill(NaN, K),
                          fill(NaN, K),
                          fill(NaN, K))


text:
# Two Factor Designs
#
# author: yizhan miao
# email: yzmiao@protonmail.com
# last update: Oct 26 2018

"""Two-way ANOVA

Syntax:
- `Result = anova2(Y::Array{T, 1}, A::Array{D, 1}, B::Array{D, 1};names=(:A, :B), preview::Bool=true) where {T <: Number, D <: Integer}`

Arguments:
- `Y`: data set in 1d vector
- `A`: corresponding factor A marker, as Array{Integer, 1}
- `B`: corresponding factor B marker, as Array{Integer, 1}

Keywords:
- `names`: the name of the factor [default: `(:A, :B)`]
- `preview`: print anova table [default: `true`]

Return:
- ANOVA Result

"""
function anova2(Y::Array{T, 1}, A::Array{D, 1}, B::Array{D, 1};names=(:A, :B), preview::Bool=true) where {T <: Number, D <: Integer}
    # check balance
    _factor_A = unique(A)
    _n_factor_A = zeros(Int64, length(_factor_A))
    for (idx, each) in enumerate(_factor_A)
        _n_factor_A[idx] = A[A .== each] |> length
    end

    _factor_B = unique(B)
    _n_factor_B = zeros(Int64, length(_factor_B))
    for (idx, each) in enumerate(_factor_B)
        _n_factor_B[idx] = B[B .== each] |> length
    end

    _rep = length(Y) / length(_factor_A) / length(_factor_B)

    @assert(length(unique(_n_factor_A)) == 1, "input data is not balanced across factor A.")
    @assert(length(unique(_n_factor_B)) == 1, "input data is not balanced across factor B.")
    #TODO: support unbalanced ANOVA

    _mean_total = mean(Y)
    _mean_A = [mean(Y[A .== _A]) for _A in _factor_A]
    _mean_B = [mean(Y[B .== _B]) for _B in _factor_B]
    _mean_BA = [mean(Y[(A .== _A) .& (B .== _B)]) for _A in _factor_A, _B in _factor_B]

    _SS_total = (Y .- _mean_total).^2 |> sum
    _SS_A = sum((_mean_A .- _mean_total).^2) * _rep * length(_factor_B)
    _SS_B = sum((_mean_B .- _mean_total).^2) * _rep * length(_factor_A)
    _SS_S_BA = [sum((Y[(A .== _A) .& (B .== _B)] .- _mean_BA[idx_A, idx_B]).^2) for (idx_A, _A) in enumerate(_factor_A), (idx_B, _B) in enumerate(_factor_B)] |> sum
    _SS_BA = _SS_total - _SS_A - _SS_B - _SS_S_BA

    _df_total = length(Y) - 1
    _df_A = length(_factor_A) - 1
    _df_B = length(_factor_B) - 1
    _df_BA = _df_A * _df_B
    _df_S_BA = _df_total - _df_A - _df_B - _df_BA

    ##### rm anova #####
    _MS_A = _SS_A / _df_A
    _MS_B = _SS_B / _df_B
    _MS_BA = _SS_BA / _df_BA
    _MS_S_BA = _SS_S_BA / _df_S_BA

    _F_A = _MS_A / _MS_S_BA
    _F_B = _MS_B / _MS_S_BA
    _F_BA = _MS_BA / _MS_S_BA

    _p_A = ccdf(FDist(_df_A, _df_S_BA), _F_A)
    _p_B = ccdf(FDist(_df_B, _df_S_BA), _F_B)
    _p_BA = ccdf(FDist(_df_BA, _df_S_BA), _F_BA)

    TOTAL = anova_entry(:total, _SS_total, _df_total, NaN, NaN, NaN)
    _A = anova_entry(names[1], _SS_A, _df_A, _MS_A, _F_A, _p_A)
    _B = anova_entry(names[2], _SS_B, _df_B, _MS_B, _F_B, _p_B)
    _BA = anova_entry(Symbol("$(names[2])*$(names[1])"), _SS_BA, _df_BA, _MS_BA, _F_BA, _p_BA)
    _S_BA = anova_entry(:error, _SS_A, _df_A, _MS_A, NaN, NaN)

    RESULT = anova_result("Two-way ANOVA", [_A, _B, _BA, _S_BA, TOTAL])

    preview&&preview_anova_result(RESULT)

    RESULT
end


text:
using SMLMData

"""
    smld, smld_preclustered, smld_combined, params = frameconnect!(smld::SMLMData.SMLD2D;
        nnearestclusters::Int=2, nsigmadev::Float64=5.0, maxframegap::Int=5, nmaxnn::Int=2)

Connect repeated localizations of the same emitter in `smld`.

# Description
Repeated localizations of the same emitter present in `smld` are connected and
combined into higher precision localizations of that emitter.  This is done by
1) forming pre-clusters of localizations, 2) estimating rate parameters from
the pre-clusters, 3) solving a linear assignment problem for connecting
localizations in each pre-cluster, and 4) combining the connected localizations
using their MLE position estimate assuming Gaussian noise.

# Inputs
- `smld`: SMLD2D structure containing the localizations that should be
          connected (see SMLMData.SMLD2D for organization/fields).
- `nnearestclusters`: Number of nearest preclusters used for local density
                      estimates. (default = 2)(see estimatedensities())
- `nsigmadev`: Multiplier of localization errors that defines a pre-clustering
               distance threshold. (default = 5)(see precluster())(pixels)
- `maxframegap`: Maximum frame gap between temporally adjacent localizations in
                 a precluster. (default = 5)(see precluster())(frames)
- `nmaxnn`: Maximum number of nearest-neighbors inspected for precluster
            membership.  Ideally, this would be set to inf, but that's not
            feasible for most data. (default = 2)(see precluster())

# Outputs
- `smld`: Input `smld` with field connectID updated to reflect connected
          localizations (however localizations remain uncombined).
- `smld_preclustered`: Copy of the input `smld` with the field connectID
                       populated to reflect the results of pre-clustering.
- `smld_combined`: Final frame-connection result (i.e., `smld` with
                   localizations that seem to be from the same blinking event
                   combined into higher precision localizations).
- `params`: Structure of parameters used in the algorithm, with some copied
            directly from the option kwargs to this function, and others
            calculated internally (see SMLMFrameConnection.ParamStruct).
"""
function frameconnect!(smld::SMLMData.SMLD2D;
    nnearestclusters::Int = 2, nsigmadev::Float64 = 5.0,
    maxframegap::Int = 5, nmaxnn::Int = 2)

    # Prepare a ParamStruct to keep track of parameters used.
    params = ParamStruct()
    params.nnearestclusters = nnearestclusters
    params.nsigmadev = nsigmadev
    params.maxframegap = maxframegap
    params.nmaxnn = nmaxnn

    # Generate pre-clusters of localizations in `smld`.
    smld_preclustered = precluster(smld, params)
    clusterdata = organizeclusters(smld_preclustered)

    # Estimate rate parameters.
    params.k_on, params.k_off, params.k_bleach, params.p_miss =
        estimateparams(smld_preclustered, clusterdata)

    # Estimate the underlying density of emitters.
    params.initialdensity =
        estimatedensities(smld_preclustered, clusterdata, params)

    # Connect localizations in `smld` by solving the LAP.
    nframes = isempty(smld.nframes) ? maximum(smld.framenum) : smld.nframes
    smld.connectID = connectlocalizations(smld_preclustered.connectID,
        clusterdata, params, nframes)

    # Combine the connected localizations into higher precision localizations.
    smld_combined = combinelocalizations(smld)

    return smld, smld_preclustered, smld_combined, params
end

"""
    smld, smld_preclustered, smld_combined, params = frameconnect(smld::SMLMData.SMLD2D;
        nnearestclusters::Int=2, nsigmadev::Float64=5.0, maxframegap::Int=5, nmaxnn::Int=2)

Connect repeated localizations of the same emitter in `smld`.

# Description
Repeated localizations of the same emitter present in `smld` are connected and
combined into higher precision localizations of that emitter.  This is done by
1) forming pre-clusters of localizations, 2) estimating rate parameters from
the pre-clusters, 3) solving a linear assignment problem for connecting
localizations in each pre-cluster, and 4) combining the connected localizations
using their MLE position estimate assuming Gaussian noise.

# Inputs
-`smld`: SMLD2D structure containing the localizations that should be
         connected (see SMLMData.SMLD2D for organization/fields).
-`nnearestclusters`: Number of nearest preclusters used for local density
                     estimates. (default = 2)(see estimatedensities())
-`nsigmadev`: Multiplier of localization errors that defines a pre-clustering
              distance threshold. (default = 5)(see precluster())(pixels)
-`maxframegap`: Maximum frame gap between temporally adjacent localizations in
                a precluster. (default = 5)(see precluster())(frames)
-`nmaxnn`: Maximum number of nearest-neighbors inspected for precluster
           membership.  Ideally, this would be set to inf, but that's not
           feasible for most data. (default = 2)(see precluster())

# Outputs
-`smld`: Input `smld` with field connectID updated to reflect connected
         localizations (however localizations remain uncombined).
-`smld_preclustered`: Copy of the input `smld` with the field connectID
                      populated to reflect the results of pre-clustering.
-`smld_combined`: Final frame-connection result (i.e., `smld` with
                  localizations that seem to be from the same blinking event
                  combined into higher precision localizations).
-`params`: Structure of parameters used in the algorithm, with some copied
           directly from the option kwargs to this function, and others
           calculated internally (see SMLMFrameConnection.ParamStruct).
"""
function frameconnect(smld::SMLMData.SMLD2D;
    nnearestclusters::Int = 2, nsigmadev::Float64 = 5.0,
    maxframegap::Int = 5, nmaxnn::Int = 2)

    return frameconnect!(deepcopy(smld);
        nnearestclusters = nnearestclusters, nsigmadev = nsigmadev,
        maxframegap = maxframegap, nmaxnn = nmaxnn)
end


text:
function basis_m(n::Int64,m::Int64)
    len = binomial(n+m-1,m)
    base = spzeros(len,n)
    counter = 1
    index = Int64[]
    for i in 1:((m+1)^n-1)
        bin_vec = splitter(string(i,base=(m+1)),n)
        if sum(bin_vec) == m
            base[counter,:] = bin_vec
            append!(index,i)
            counter += 1
        end
    end
    return base, index
end

function splitter(x,n)
    vectorized = zeros(Int64,0)
    while length(x) < n
        x = "0"*x
    end
    for i in x
        str_to_int = parse(Int, i)
        append!(vectorized, str_to_int)
    end
    return vectorized
end

function myfind(c,j)
    a = similar(c, Int)
    count = 1
    @inbounds for i in eachindex(c)
        a[count] = i
        count += (c[i] == j)
    end
    return a[1:(count-1)]
end

function bdb(n::Int64, m::Int64, i::Int64, j::Int64)
    base, ind = basis_m(n,m)
    l = binomial(n+m-1,m)
    op = spzeros(l,l)
    if i==j
        for k in 1:l
            if base[k,:][j] != 0
                op[k,k] = base[k,:][j]
            end
        end
    else
        for k in 1:l
           if base[k,:][j] != 0
                ill = Int(ind[k]-(m+1)^(n-j)+(m+1)^(n-i))
                ill2 = myfind(ind, ill)[1]
                op[ill2,k] = 1
            end
        end
    end
    return op
end

function bbd(n::Int64, m::Int64, i::Int64, j::Int64)
    return bdb(n, m, j, i)
end

function bdb(base::SparseMatrixCSC{Float64,Int64}, ind::Array{Int64,1}, m::Int64, i::Int64, j::Int64)
    l = size(base)[1]
    n = size(base)[2]
    op = spzeros(l,l)
    if i==j
        for k in 1:l
            if base[k,:][j] != 0
                op[k,k] = base[k,:][j]
            end
        end
    else
        for k in 1:l
           if base[k,:][j] != 0
                ill = Int(ind[k]-(m+1)^(n-j)+(m+1)^(n-i))
                ill2 = myfind(ind, ill)[1]
                op[ill2,k] = 1
            end
        end
    end
    return op
end

function bbd(base::SparseMatrixCSC{Float64,Int64}, ind::Array{Int64,1}, m::Int64, i::Int64, j::Int64)
    return bdb(base, ind, m, j, i)
end

struct Op_fixed
    dim::Int
    nume::Int
    bdbtot::SparseMatrixCSC{Float64,Int64}
    le::Int
    basis::SparseMatrixCSC{Float64,Int64}
    Op_fixed(dim,nume) = new(dim, nume, operators_fixed(dim,nume)...)
end

dim(o::Op_fixed) = o.dim
nume(o::Op_fixed) = o.nume
bdbtot(o::Op_fixed) = o.bdbtot
le(o::Op_fixed) = o.le
basis(o::Op_fixed) = o.basis

bdb(o::Op_fixed,i,j) = bdbtot(o)[(1+(i-1)*le(o)):i*le(o),(1+(j-1)*le(o)):j*le(o)]
bbd(o::Op_fixed,i,j) = bdb(o,j,i)

# Here we create the full matrix with all
# fixed particle operators definitions
function operators_fixed(n::Int,m::Int)
    l1 = binomial(n+m-1,m)
    l2 = n*l1
    bas,ind = basis_m(n,m)
    z = spzeros(l2,l2)
    for i in 1:n
        for j in i:n
            if i==j
                z[(1+(i-1)*l1):i*l1,(1+(j-1)*l1):j*l1] = 1/2*bdb(bas,ind, m, i,j)
            else
                z[(1+(i-1)*l1):i*l1,(1+(j-1)*l1):j*l1] = bdb(bas,ind, m, i,j)
            end
        end
    end
    op_general = z+z'
    return op_general, l1, bas
end

function bd(n::Int64, m::Int64, i::Int64)
    base1, ind1 = basis_m(n,m)
    l2 = binomial(n+m-1,m)
    base2, ind2 = basis_m(n,m+1)
    l1 = binomial(n+m,m+1)
    op = spzeros(l1,l2)
    s2 = sparse([(m+2)^l for l in 0:(n-1)])
    for k in 1:l2
        j = myfind(ind2, (s2[end:-1:1]'*base1[k,:])+(m+2)^(n-i))[1]
        op[j,k] = 1
    end
    if m==0
        return op
    end
    return op
end

function b(n::Int64, m::Int64, i::Int64)
    op = bd(n, m-1, i)
    return op'
end

# This serves for applying several creations op
# at once. Modes will be a vector with the chosen
# modes
function bd(n::Int64, m::Int64, modes::Array{Int64,1})
    l = length(modes)
    op  = bd(n,m,modes[1])
    for k in 2:l
        op = bd(n,m+k-1,modes[k])*op
    end
    return op
end

function b(n::Int64, m::Int64, modes::Array{Int64,1})
    l = length(modes)
    op  = b(n,m,modes[1])
    for k in 2:l
        op = b(n,m-k+1,modes[k])*op
    end
    return op
end

function state_form(v::SparseMatrixCSC{Float64, Int64})
  i, j, w = findnz(v)
  l = size(v)[1]
  return sparsevec(i,w,l)
end


text:
using Plots

# Model
include_model("quadruped_v2")
model = free_time_model(model)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)
# render(vis)

# Horizon
T = 61
Tm = 31
T_fix = 0#5

# Time step
tf = 1.0
h = tf / (T - 1)

function ellipse_traj(x_start, x_goal, z, T)
	dist = x_goal - x_start
	a = 0.5 * dist
	b = z
	z̄ = 0.0
	# x = range(x_start, stop = x_goal, length = T)
	x = circular_projection_range(x_start, stop = x_goal, length = T)
	z = sqrt.(max.(0.0, (b^2) * (1.0 .- ((x .- (x_start + a)).^2.0) / (a^2.0))))
	return x, z
end

function circular_projection_range(start; stop=1.0, length=10)
	dist = stop - start
	θr = range(π, stop=0, length=length)
	r = start .+ dist * ((1 .+ cos.(θr))./2)
	return r
end

function initial_configuration(model::QuadrupedV2; offset = 0.025)
	q1 = zeros(model.nq)

	# position
	q1[1] = model.l_torso
	q1[3] = 0.2

	# orientation
	mrp = MRP(RotZ(0.0))
	q1[4:6] = [mrp.x; mrp.y; mrp.z]

	# feet positions (in body frame)
	q1[7:9] = [2 * model.l_torso + offset; model.w_torso; 0.0]
	q1[10:12] = [2 * model.l_torso - offset; -model.w_torso; 0.0]
	q1[13:15] = [-offset; model.w_torso; 0.0]
	q1[16:18] = [offset; -model.w_torso; 0.0]

	return q1
end

# function intermediate_configuration(model::QuadrupedV2; offset = 0.025)
# 	q1 = zeros(model.nq)
#
# 	# position
# 	q1[1] = model.l_torso
# 	q1[3] = 0.2
#
# 	# orientation
# 	mrp = MRP(RotZ(0.0))
# 	q1[4:6] = [mrp.x; mrp.y; mrp.z]
#
# 	# feet positions (in body frame)
# 	q1[7:9] = [2 * model.l_torso; model.w_torso; 0.0]
# 	q1[10:12] = [2 * model.l_torso - 2 * offset; -model.w_torso; 0.0]
# 	q1[13:15] = [-2 * offset; model.w_torso; 0.0]
# 	q1[16:18] = [0.0; -model.w_torso; 0.0]
#
# 	return q1
# end

function intermediate_configuration(model::QuadrupedV2; offset = 0.025)
	q1 = zeros(model.nq)

	# position
	q1[1] = 0.0
	q1[3] = 0.2 + model.l_torso

	# orientation
	mrp = MRP(RotY(-0.5 * π))
	q1[4:6] = [mrp.x; mrp.y; mrp.z]

	# feet positions (in body frame)
	q1[7:9] = [0.0; model.w_torso; q1[3]]
	q1[10:12] = [0.0; -model.w_torso; q1[3]]
	q1[13:15] = [-offset; model.w_torso; 0.0]
	q1[16:18] = [offset; -model.w_torso; 0.0]

	return q1
end

function final_configuration(model::QuadrupedV2; offset = 0.025)
	q1 = zeros(model.nq)

	# position
	q1[1] = 0.0
	q1[3] = 0.2 + model.l_torso

	# orientation
	mrp = MRP(RotY(-0.5 * π))
	q1[4:6] = [mrp.x; mrp.y; mrp.z]

	# feet positions (in body frame)
	q1[7:9] = [0.0; model.w_torso; q1[3]]
	q1[10:12] = [0.0; -model.w_torso; q1[3]]
	q1[13:15] = [-offset; model.w_torso; 0.0]
	q1[16:18] = [offset; -model.w_torso; 0.0]

	return q1
end

q1 = initial_configuration(model)
qT = intermediate_configuration(model)
# qT = final_configuration(model)

visualize!(vis, model, [q1])
# visualize!(vis, model, [qM])
visualize!(vis, model, [qT])

# feet positions
pf1 = q1[6 .+ (1:3)]
pf2 = q1[9 .+ (1:3)]

pf3 = q1[12 .+ (1:3)]
pf4 = q1[15 .+ (1:3)]

q_ref = [q1, linear_interpolation(q1, qT, T)...]
visualize!(vis, model, q_ref)

# zh = 0.025
# xf1_el, zf1_el = ellipse_traj(pf1[1], pf1[1] + strd, zh, Tm - T_fix)
# xf1 = [[pf1[1] for t = 1:Tm + T_fix]..., xf1_el[2:end]...]
# zf1 = [[pf1[3] for t = 1:Tm + T_fix]..., zf1_el[2:end]...]
# pf1_ref = [[xf1[t]; pf1[2];  zf1[t]] for t = 1:T]

# xf4_el, zf4_el = ellipse_traj(pf4[1], pf4[1] + strd, zh, Tm - T_fix)
# xf4 = [[pf4[1] for t = 1:Tm + T_fix]..., xf4_el[2:end]...]
# zf4 = [[pf4[3] for t = 1:Tm + T_fix]..., zf4_el[2:end]...]
# pf4_ref = [[xf4[t]; pf4[2]; zf4[t]] for t = 1:T]

# xf2_el, zf2_el = ellipse_traj(pf2[1], pf2[1] + strd, zh, Tm - T_fix)
# xf2 = [[xf2_el[1] for t = 1:T_fix]..., xf2_el..., [xf2_el[end] for t = 1:Tm-1 + T_fix]...]
# zf2 = [[zf2_el[1] for t = 1:T_fix]..., zf2_el..., [zf2_el[end] for t = 1:Tm-1 + T_fix]...]
# pf2_ref = [[xf2[t]; pf2[2]; zf2[t]] for t = 1:T]

# xf3_el, zf3_el = ellipse_traj(pf3[1], pf4[1], zh, Tm - T_fix)
# xf3 = [[xf3_el[1] for t = 1:T_fix]..., xf3_el..., [xf3_el[end] for t = 1:Tm-1]...]
# zf3 = [[zf3_el[1] for t = 1:T_fix]..., zf3_el..., [zf3_el[end] for t = 1:Tm-1]...]
# pf3_ref = [[xf3[t]; pf3[2]; zf3[t]] for t = 1:T]
#
# tr = range(0, stop = tf, length = T)
# plot(tr, hcat(pf1_ref...)')
# plot!(tr, hcat(pf4_ref...)')
#
# plot!(tr, hcat(pf2_ref...)')
# plot(tr, hcat(pf3_ref...)')

# Bounds

# control
# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
_uu[end] = 1.0 * h
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
_ul[end] = 1.0 * h
ul, uu = control_bounds(model, T, _ul, _uu)

xl, xu = state_bounds(model, T,
    x1 = [q1; q1],
	xT = [qT; qT])

# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])

# render(vis)
visualize!(vis, model, q_ref, Δt = h)
x_ref = configuration_to_state(q_ref)
x0 = copy(x_ref)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e4, model.m - 1)

# quadratic tracking objective
# Σ (x - xref)' Q (x - x_ref) + (u - u_ref)' R (u - u_ref)
obj_control = quadratic_time_tracking_objective(
    [1.0 * Diagonal(1.0e-1 * ones(model.n)) for t = 1:T],
    [1.0 * Diagonal([1.0e-5 * ones(model.nu)..., 1.0e-3 * ones(model.nc + model.nb)..., zeros(model.m - model.nu - model.nc - model.nb)...]) for t = 1:T-1],
    [x_ref[t] for t = 1:T],
    [zeros(model.m) for t = 1:T],
    1.0)

# quadratic velocity penalty
#Σ v' Q v
v_penalty = 1.0 * ones(model.nq)
obj_velocity = velocity_objective(
    [h * Diagonal(v_penalty) for t = 1:T-1],
    model.nq,
    h = h,
    idx_angle = collect([3, 4, 5]))

obj_ctrl_velocity = control_velocity_objective(Diagonal([1.0 * ones(model.nu)..., 1.0e-3 * ones(model.nc + model.nb)..., zeros(model.m - model.nu - model.nc - model.nb)...]))

function l_stage(x, u, t)
	q1 = view(x, 1:18)
	q2 = view(x, 18 .+ (1:18))
    J = 0.0

	# torso height
	# J += 1000.0 * (q2[t] - q_ref[t])^2.0

	# orientation
	# J += 1000.0 * sum((q2[4:6] - q_ref[t][4:6]).^2.0)

	# forward velocity
	# J += 1000.0 * ((q2[1] - q1[1]) / h - strd / tf)^2.0

    if true
		J += 1000.0 * sum((q_ref[t][1:6] - q2[1:6]).^2.0)

		p1 = q2[6 .+ (1:3)]
		p2 = q2[9 .+ (1:3)]
		p3 = q2[12 .+ (1:3)]
		p4 = q2[15 .+ (1:3)]

		J += 10000.0 * sum((q_ref[t][6 .+ (1:3)] - p1).^2.0)
		J += 10000.0 * sum((q_ref[t][9 .+ (1:3)] - p2).^2.0)
		J += 100.0 * sum((q_ref[t][12 .+ (1:3)] - p3).^2.0)
		J += 100.0 * sum((q_ref[t][15 .+ (1:3)] - p4).^2.0)
	end

    return J
end

l_terminal(x) = 0.0
obj_shaping = nonlinear_stage_objective(l_stage, l_terminal)

obj = MultiObjective([obj_penalty,
                      obj_control,
                      obj_velocity,
					  obj_shaping,
					  obj_ctrl_velocity])

# Constraints
include_constraints(["stage", "contact", "free_time", "loop"])
function pinned1!(c, x, u, t)
    q = view(x, 1:18)
    # c[1:3] = pf1_ref[t] - q[6 .+ (1:3)]
    c[1:3] = q_ref[t][15 .+ (1:3)] - q[15 .+ (1:3)]
    nothing
end

function pinned2!(c, x, u, t)
    q = view(x, 1:18)
	# c[1:3] = pf2_ref[t] - q[9 .+ (1:3)]
    c[1:3] = q_ref[t][12 .+ (1:3)] - q[12 .+ (1:3)]
    nothing
end

n_stage = 3
# t_idx1 = vcat([t for t = 1:Tm + T_fix])
# t_idx2 = vcat([t for t = 1:T_fix]..., [t for t = Tm:T]...)
t_idx1 = vcat([t for t = 1:T])
t_idx2 = vcat([t for t = 1:T])
# con_pinned1 = stage_constraints(pinned1!, n_stage, (1:0), t_idx1)
# con_pinned2 = stage_constraints(pinned2!, n_stage, (1:0), t_idx2)
con_pinned1 = stage_constraints(pinned1!, n_stage, (1:0), t_idx1)
con_pinned2 = stage_constraints(pinned2!, n_stage, (1:0), t_idx2)

con_contact = contact_constraints(model, T)
con_free_time = free_time_constraints(T)
# con_loop1 = loop_constraints(model, collect(1:2model.nq), 1, Tm, shift = [q_shift1; q_shift1])
# con_loop2 = loop_constraints(model, collect(1:2model.nq), Tm, T, shift = [q_shift2; q_shift2])

con = multiple_constraints([con_contact,
	# con_loop1, con_loop2,
    con_free_time, con_pinned1, con_pinned2])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# trajectory initialization
u0 = [[1.0e-3 * rand(model.m - 1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
@time z̄, info = solve(prob, copy(z0),
    nlp = :ipopt,
    tol = 1.0e-2, c_tol = 1.0e-2,
	max_iter = 2000,
    time_limit = 60 * 2, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)
_tf, _t, h̄ = get_time(ū)
@show h̄[1]

# vis = Visualizer()
# render(vis)
visualize!(vis, model, [[q̄[1] for t = 1:50]..., q̄..., [q̄[end] for t = 1:50]...], Δt = h̄[1])

# check foot trajectories
_pf1_ref = [q[6 .+ (1:3)] for q in q̄]
_pf2_ref = [q[9 .+ (1:3)]  for q in q̄]

_pf3_ref = [q[12 .+ (1:3)]  for q in q̄]
_pf4_ref = [q[15 .+ (1:3)]  for q in q̄]

plot(hcat([q[6 .+ (1:3)] for q in q_ref]...)', width = 2.0, color = :black)
plot!(hcat(_pf1_ref...)', color = :red)

plot(hcat([q[9 .+ (1:3)]  for q in q_ref]...)', width = 2.0, color = :black)
plot!(hcat(_pf2_ref...)', color = :red)

plot(hcat([q[12 .+ (1:3)]  for q in q_ref]...)', width = 2.0, color = :black)
plot!(hcat(_pf3_ref...)', color = :red)

plot(hcat([q[15 .+ (1:3)]  for q in q_ref]...)', width = 2.0, color = :black)
plot!(hcat(_pf4_ref...)', color = :red)

# check control trajectory
plot(hcat(ū...)[1:model.nu, :]', linetype = :steppost, label = "")

# unpack solution
q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
ψ = [u[model.idx_ψ] for u in ū]
η = [u[model.idx_η] for u in ū]
hm = mean(h̄)
μm = model.μ

qm = q; um = u; γm = γ; bm = b; ψm = ψ; ηm = η;
@save joinpath(@__DIR__, "quadruped_v2_mirror_gait.jld2") qm um γm bm ψm ηm μm hm

plot(hcat(q...)')
plot(hcat(q...)')

[norm(fd(model, [q[t+1]; q[t+2]], [q[t]; q[t+1]], [u[t]; γ[t]; b[t]; hm], zeros(model.d), hm, t)) for t = 1:T-1]

# function mirror_gait(q, u, γ, b, ψ, η, T)
# 	qm = [deepcopy(q)...]
# 	um = [deepcopy(u)...]
# 	γm = [deepcopy(γ)...]
# 	bm = [deepcopy(b)...]
# 	ψm = [deepcopy(ψ)...]
# 	ηm = [deepcopy(η)...]
#
# 	stride = zero(qm[1])
# 	@show stride[1] = q[T+1][1] - q[2][1]
# 	@show 0.5 * strd
#
# 	p10 = q[2][6 .+ (1:3)]
# 	p20 = q[2][9 .+ (1:3)]
# 	p30 = q[2][12 .+ (1:3)]
# 	p40 = q[2][15 .+ (1:3)]
#
# 	@show norm(fd(model, [qm[end-1]; qm[end]], [qm[end-2]; qm[end-1]], [um[end]; γm[end]; bm[end]], zeros(model.d), hm, 1))
#
# 	for t = 1:1#T-1
# 		# configuration
# 		pt = q[t+2][1:3]
# 		a = q[t+2][3 .+ (1:3)]
# 		p1 = q[t+2][6 .+ (1:3)]
# 		# p2 = q[t+2][9 .+ (1:3)]
# 		# p3 = q[t+2][12 .+ (1:3)]
# 		p4 = q[t+2][15 .+ (1:3)]
#
# 		p2_diff = q[t+2][9 .+ (1:3)] - p20
# 		p3_diff = q[t+2][12 .+ (1:3)] - p30
#
# 		push!(qm, [pt + [0.5 * strd; 0.0; 0.0]; a;
# 			p10[1] + p2_diff[1]; q[t+2][8]; p10[3] + p2_diff[3]
# 			q[end][9 .+ (1:3)];
# 			q[end][12 .+ (1:3)];
# 			p40[1] + p3_diff[1]; q[t+2][17]; p40[3] + p3_diff[3];
# 			])
# 		# push!(qm, q[t+2])
#
# 		# control
# 		u1 = u[t][1:3]
# 		u2 = u[t][3 .+ (1:3)]
# 		u3 = u[t][6 .+ (1:3)]
# 		u4 = u[t][9 .+ (1:3)]
# 		push!(um, [u2; u1; u4; u3])
#
# 		# impact
# 		γ1 = γ[t][1]
# 		γ2 = γ[t][2]
# 		γ3 = γ[t][3]
# 		γ4 = γ[t][4]
# 		push!(γm, [γ2; γ1; γ4; γ3])
# 		# push!(γm, γ[t])
#
# 		# friction
# 		b1 = b[t][1:4]
# 		b2 = b[t][4 .+ (1:4)]
# 		b3 = b[t][8 .+ (1:4)]
# 		b4 = b[t][12 .+ (1:4)]
# 		push!(bm, [b2; b1; b4; b3])
# 		# push!(bm, b[t])
#
# 		# dual
# 		ψ1 = ψ[t][1]
# 		ψ2 = ψ[t][2]
# 		ψ3 = ψ[t][3]
# 		ψ4 = ψ[t][4]
# 		push!(ψm, [ψ2; ψ1; ψ4; ψ3])
# 		# push!(ψm, ψ[t])
#
# 		# dual
# 		η1 = η[t][1:4]
# 		η2 = η[t][4 .+ (1:4)]
# 		η3 = η[t][8 .+ (1:4)]
# 		η4 = η[t][12 .+ (1:4)]
# 		push!(ηm, [η2; η1; η4; η3])
# 		# push!(ηm, η[t])
#
# 		@show norm(fd(model, [qm[end-1]; qm[end]], [qm[end-2]; qm[end-1]], [um[end]; γm[end]; bm[end]], zeros(model.d), hm, 1))
# 	end
#
# 	return qm, um, γm, bm, ψm, ηm
# end
#
# qm, um, γm, bm, ψm, ηm = mirror_gait(q, u, γ, b, ψ, η, T)
#
# # @save joinpath(@__DIR__, "quadruped_v2_mirror_gait.jld2") qm um γm bm ψm ηm μm hm
#
plot(hcat(q...)', color = :black, width = 2.0, label = "")
plot!(hcat(qm...)', color = :red, width = 1.0, label = "")
#
plot(hcat(u...)', color = :black, width = 2.0, label = "", linetype = :steppost)
# plot!(hcat(um...)', color = :red, width = 1.0, label = "", linetype = :steppost)
#
plot(hcat(γ...)', color = :black, width = 2.0, label = "", linetype = :steppost)
# plot!(hcat(γm...)', color = :red, width = 1.0, label = "", linetype = :steppost)
#
plot(hcat(b...)', color = :black, width = 2.0, label = "", linetype = :steppost)
# plot!(hcat(bm...)', color = :red, width = 1.0, label = "", linetype = :steppost)
# #
# function get_q_viz(q̄)
# 	q_viz = [q̄...]
# 	shift_vec = zeros(model.nq)
# 	shift_vec[1] = q̄[end][1]
# 	for i = 1:5
# 		q_update = [q + shift_vec for q in q̄[2:end]]
# 		push!(q_viz, q_update...)
# 		shift_vec[1] = q_update[end][1]
# 	end
#
# 	return q_viz
# end

# vis = Visualizer()
# render(vis)
visualize!(vis, model,
	qm,
	Δt = h̄[1])

settransform!(vis["/Cameras/default"],
        compose(Translation(0.0, -25.0, -1.0), LinearMap(RotY(0.0 * π) * RotZ(-π / 2.0))))
setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 20)
# open(vis)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    NeighborSearchMethod

A method for searching neighbors given a reference point.
"""
abstract type NeighborSearchMethod end

"""
    search!(neighbors, pₒ, method; mask=nothing)

Update `neighbors` of point `pₒ` using `method` and return
number of neighbors found. Optionally, specify a `mask` for
all indices of the domain.
"""
function search! end

"""
    search(pₒ, method, mask=nothing)

Return neighbors of point `pₒ` using `method`. Optionally,
specify a `mask` for all indices of the domain.
"""
function search end

"""
    BoundedNeighborSearchMethod

A method for searching neighbors with the property that the number of neighbors
is bounded above by a known constant (e.g. k-nearest neighbors).
"""
abstract type BoundedNeighborSearchMethod <: NeighborSearchMethod end

"""
    maxneighbors(method)

Return the maximum number of neighbors obtained with `method`.
"""
function maxneighbors end

# ----------
# FALLBACKS
# ----------

function search(pₒ::Point, method::BoundedNeighborSearchMethod; mask=nothing)
  neighbors = Vector{Int}(undef, maxneighbors(method))
  nneigh = search!(neighbors, pₒ, method; mask=mask)
  view(neighbors, 1:nneigh)
end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("neighborsearch/ball.jl")
include("neighborsearch/knearest.jl")
include("neighborsearch/kball.jl")
include("neighborsearch/bounded.jl")


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 377dacc0-a3c5-11eb-2818-ef550190e23d
begin
	using Plots, PlutoUI, CommonMark
	plotlyjs()
end

# ╔═╡ b615f2e2-a411-4f17-933c-9ee0a15adc6a
md"""
``V`` $(@bind V Slider(-10:0.1:10,show_value=true,default=1))

``x`` $(@bind xi Slider(0.1:0.01:2,show_value=true,default=0))
"""

# ╔═╡ 58de3ac6-30a0-4ace-bc83-11ed91009a54
p(x) = V/x

# ╔═╡ cc8973a6-669b-460f-b515-48548c8bbc38
begin
r0 = range(0.05,stop=2,length=1000)
plot(r0,p.(r0),label=false,xlims=[0,2],ylim=sort([sign(V)*20,0]))
y = p(xi)
vline!([xi],label=false)
hline!([y],label=false)
scatter!([(xi,y)],label=false)
end

# ╔═╡ bb59f36f-d6c3-4afa-b5a1-0069e916a739
md"""
``A`` $(@bind A Slider(-10:0.1:10,show_value=true,default=1))

``B`` $(@bind B Slider(-2:0.1:2,show_value=true,default=0))

``C`` $(@bind C Slider(-10:0.1:10,show_value=true,default=0))
"""

# ╔═╡ 5ceac717-6fa2-4d20-9384-f90d63c94e55
f(x) = A / (x-B) + C

# ╔═╡ 44076b9e-d547-4fdd-a2e6-eebcf37175c2
Markdown.parse("""
	\`\`\`math
	f(x) = \\frac{$A}{x-$B}+$C
	\`\`\`
	""")

# ╔═╡ 52ef0e0f-cc79-44f9-b3a0-8779b831cef3
begin
r = range(-5,stop=5,length=500)
plot(r,f.(r),label=false,xlim=[-4,4],ylim=[-10,10])
hline!([C],label="y=C")
vline!([B],label="x=B")
end

# ╔═╡ 229da4be-7ef8-457a-924c-5802516e4328
begin
a(x) = x + 2
den(x) = x + 3
end

# ╔═╡ 39b33678-ab2f-4236-b3ec-adc157850498
rango = 1:1000

# ╔═╡ 97ba7812-fc14-49df-be09-5b63cddeac3e
a.(rango)

# ╔═╡ 349850cf-ec77-49ca-bce0-7e63b842be14
den.(rango)

# ╔═╡ cb875d9c-f590-4d30-86b7-0475ea206195
h(x) = a(x)/den(x)

# ╔═╡ f18fcc1b-ab88-4ab4-920b-814a2a7d5f3a
h.(rango)

# ╔═╡ 59c64275-2fee-489a-893f-1b100b2c3749
plot([a,den,h])

# ╔═╡ Cell order:
# ╟─377dacc0-a3c5-11eb-2818-ef550190e23d
# ╟─b615f2e2-a411-4f17-933c-9ee0a15adc6a
# ╟─58de3ac6-30a0-4ace-bc83-11ed91009a54
# ╟─cc8973a6-669b-460f-b515-48548c8bbc38
# ╟─bb59f36f-d6c3-4afa-b5a1-0069e916a739
# ╟─5ceac717-6fa2-4d20-9384-f90d63c94e55
# ╟─44076b9e-d547-4fdd-a2e6-eebcf37175c2
# ╠═52ef0e0f-cc79-44f9-b3a0-8779b831cef3
# ╠═229da4be-7ef8-457a-924c-5802516e4328
# ╠═39b33678-ab2f-4236-b3ec-adc157850498
# ╠═97ba7812-fc14-49df-be09-5b63cddeac3e
# ╠═349850cf-ec77-49ca-bce0-7e63b842be14
# ╠═cb875d9c-f590-4d30-86b7-0475ea206195
# ╠═f18fcc1b-ab88-4ab4-920b-814a2a7d5f3a
# ╠═59c64275-2fee-489a-893f-1b100b2c3749

max_stars_repo_name: akielbowicz/cbc-analisis-ingenieria

max_issues_repo_name: akielbowicz/cbc-analisis-ingenieria

max_forks_repo_name: akielbowicz/cbc-analisis-ingenieria


text:
"""

Need to turn on compiled mode by `C` to run through DataFrames code.
Cannot use `u` when DataFrames is used in the code because it is super slow.
Use `bp add line` to set break point instead of `u line`.
"""
#
# Use ` to enter Julia mode to exame variables
# Use backspace to go back to debug mode
#
using Debugger
db = Debugger

"""
    foo(n)

Dummy function for learning Debugger.
"""
function foo(n)
    x = n+1
    m = (BigInt[1 1; 1 0])^x
    Ans = m[2,1]
    #@bp
    # bp add 22
    return Ans
end

#db.@run foo(3)
#db.@enter foo(3)

nothing


text:
using ImageFiltering: Kernel, imfilter, kernelfactors, centered

function nld_step_scalar(Ld, c, stepsize)
    dx = 0.5 * stepsize * (c[1:end, 1:end-1] + c[1:end, 2:end]) .* (Ld[1:end, 2:end] - Ld[1:end, 1:end-1])
    dy = 0.5 * stepsize * (c[1:end-1, 1:end] + c[2:end, 1:end]) .* (Ld[2:end, 1:end] - Ld[1:end-1, 1:end])

    Ld[1:end, 1:end-1] .+= dx
    Ld[1:end, 2:end] .-= dx
    Ld[1:end-1, 1:end] .+= dy
    Ld[2:end, 1:end] .-= dy
end

pm_g1_diffusivity(Lx, Ly, k) =
    calculate_diffusivity′(Lx, Ly, k) do dL
        exp(-dL)
    end

pm_g2_diffusivity(Lx, Ly, k) =
    calculate_diffusivity′(Lx, Ly, k) do dL
        1.0 / (1.0 + dL)
    end

weickert_diffusivity(Lx, Ly, k) =
    calculate_diffusivity′(Lx, Ly, k) do dL
        1.0 - exp(-3.315 / (dL * dL * dL * dL))
    end

charbonnier_diffusivity(Lx, Ly, k) =
    calculate_diffusivity′(Lx, Ly, k) do dL
        1.0 / sqrt(1.0 + dL)
    end

function calculate_diffusivity′(func, Lx, Ly, k)
    dst = zeros(size(Lx))
    invk2 = 1.0 / (k * k)
    nrow, ncol = size(Lx)
    @inbounds begin
        @simd for j = 1:nrow
            @simd for k = 1:ncol
                lx, ly = Lx[j, k], Ly[j, k]
                dL = (lx * lx + ly * ly) * invk2
                dst[j, k] = func(dL)
            end
        end
    end
    dst
end

function compute_k_percentile(img, perc; gscale = 1.0, nbins = 300)
    Lsmooth = imfilter(img, Kernel.gaussian(gscale))
    fy, fx = Kernel.scharr()
    fy .*=32
    fx .*=32
    # fx,fy = Kernel.ando3()

    w = (size(fx)[1] + 1) ÷ 2
    Lx = @view imfilter(Lsmooth, fx)[w:end-w+1, w:end-w+1]
    Ly = @view imfilter(Lsmooth, fy)[w:end-w+1, w:end-w+1]
    compute_k_percentile′(Lx, Ly, perc, nbins)
end

function compute_k_percentile′(Lx, Ly, perc, nbins = 300)
    hist = zeros(Int32, nbins)

    hmax = sqrt(maximum(y -> sum(x -> x^2, y), zip(Lx[:], Ly[:])))

    for (lx, ly) in zip(Lx[:], Ly[:])
        modg = sqrt(lx^2 + ly^2)
        if modg > 1e-10
            nbin = ceil(Int, nbins * (modg / hmax))  # limiting to nbins should not be necessary
            hist[nbin] += 1
        end
    end

    nthreshold = floor(Int, sum(hist) * perc)

    k = findfirst(hx -> hx > nthreshold, cumsum(hist))

    return if k == nothing
        0.03
    else
        (hmax * k) / nbins
    end
end

function demo_diffusivity_functions()
    imshow(Lsmooth)
    aa = pm_g1_diffusivity(Lx, Ly, 0.01)
    imshow(aa)
    aa = pm_g2_diffusivity(Lx, Ly, 0.01)
    imshow(aa)
    aa = weickert_diffusivity(Lx, Ly, 0.01)
    imshow(aa)
    aa = charbonnier_diffusivity(Lx, Ly, 0.01)
    imshow(aa)
end

function demo_k_percentile()
    gthresh = compute_k_percentile(img, 0.95)
    imshow(sqrt.(Lx .^ 2 + Ly .^ 2))
    imshow(sqrt.(Lx .^ 2 + Ly .^ 2) .< gthresh)
end

function demo_nld()
    Lflow = charbonnier_diffusivity(Lx, Ly, 0.01)

    imshow(copy(Lt))
    for _ = 1:3
        for _ = 1:100
            nld_step_scalar(Lt, Lflow, 0.05)
        end
        imshow(copy(Lt))
    end
end

function compute_derivative_kernels(scale)

    ksize = 3 + 2*(scale-1)

    ## The usual Scharr kernel
    if (scale == 1)
        # Kernel.scharr()
        fx, fy = Kernel.scharr()
        # (fx.*32, fy.*32)
        (fx, fy)
    end

    kx = centered(zeros(ksize))
    ky = centered(zeros(ksize))

    w = 10.0/3.0
    norm = 1.0/(2.0*scale*(w+2.0))

    # ky[1+0] = norm
    # ky[1+ksize÷2] = w*norm
    # ky[1+ksize-1] = norm

    # kx[1+0] = -1
    # kx[1+ksize÷2] = 0
    # kx[1+ksize-1] = 1
    ky[-ksize÷2] = norm
    ky[0] = w*norm
    ky[ksize÷2] = norm

    kx[-ksize÷2] = -1
    kx[0] = 0
    kx[ksize÷2] = 1
    # (kernelfactors((ky, kx)), kernelfactors((kx, ky)))
    (kernelfactors((ky, kx)), kernelfactors((kx, ky)))
end


text:
#!/usr/bin/env julia
using FASTX
using BioAlignments

# Read sequence dictionary for HLA
seqFile = "data/S2_callseqs.fa"
name2seq = Dict()
open(FASTA.Reader,seqFile) do reader
	for record in reader
		name2seq[FASTA.identifier(record)]=FASTA.sequence(record)
	end
end

target = "HLA:HLA00602" # >HLA:HLA00602 DQA1*01:02:01:01 6484 bp
chrStart=32628178 # DQA1
chrEnd=32637531
hlaseq = name2seq[target]

# Read chromosome 6
chr6File= "data/chr6_hg38.fa"
chr6record = FASTA.Record()
reader = open(FASTA.Reader,chr6File)
read!(reader,chr6record)
close(reader)

# Pull out region of chromosome6
chr6seq = FASTA.sequence(chr6record)
hlaref = chr6seq[chrStart:chrEnd]

# Align that region to big HLA sequence
scoremodel = AffineGapScoreModel(EDNAFULL, gap_open=-5, gap_extend=-1);
res = pairalign(LocalAlignment(),hlaref,hlaseq,scoremodel)

print(res)

# HLA-DQA1
# ./cutandalign.jl  2.58s user 0.86s system 101% cpu 3.398 total

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
using Insolation
using Plots
using Dates
using Statistics
using Formatting

using CLIMAParameters
using CLIMAParameters.Planet
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using CLIMAParameters: AbstractParameterSet
const APS = AbstractParameterSet

function diurnal_cycle(lat, lon, date, timezone, filename)
    nhours = 1000
    hours = collect(range(0, stop = 24, length = nhours))
    insol = zeros(nhours)
    sza = zeros(nhours)
    for (i,hr) in enumerate(hours)
        h = Int(round(hr + timezone))
        m = Int(round((hr + timezone - h) * 60))
        datetime = date + Dates.Hour(h) + Dates.Minute(m)
        S, mu = solar_flux_and_cos_sza(datetime, lon, lat, param_set)
        insol[i] = S*mu
        sza[i] = rad2deg(acos(mu))
    end
    plot(hours, insol, color=:blue, lw=2, label="",
            ylabel="Insolation [W/m2]", yguidefontcolor=:blue, dpi=200,
            left_margein=10Plots.mm, right_margin=15Plots.mm)
    plot!(twinx(), hours, sza, color=:red, lw=2, label="",
            ylabel="SZA [deg]", yguidefontcolor=:red)
    xlabel!("Local Time")
    savefig(filename)
end


text:
"""
    struct AggregationV{Tf<:Function} <: AbstractVApproximator{Int}
        table::Vector{Float64}
        f::Tf
    end

Using `a.f` to map a state `s` into an `Int`, then use `a.table` to check the corresponding state value.
"""
struct AggregationV{Tf<:Function} <: AbstractVApproximator{Int}
    table::Vector{Float64}
    f::Tf
end

(a::AggregationV)(s) = a.table[a.f(s)]

update!(a::AggregationV, s, e) = a.table[a.f(s)] += e


text:
struct IterativeSolvers_LOBPCG{Tv} <: AbstractEigenMethod{Tv}
    precond::Matrix{Tv}
    IterativeSolvers_LOBPCG{Tv}(h::AbstractMatrix{Tv}, nev = 1) where {Tv} =
        new(rand(Tv, size(h,1), nev))
end

IterativeSolvers_LOBPCG(h::AbstractMatrix{Tv}, nev = 1) where {Tv} =
    IterativeSolvers_LOBPCG{Tv}(h, nev)

function (d::Diagonalizer{<:IterativeSolvers_LOBPCG, Tv})(nev::Integer;
        side = upper, precond = true, kw...) where {Tv}
    if size(d.method.precond) != (size(d.matrix, 1), nev)
        d.method = IterativeSolvers_LOBPCG(d.matrix, nev) # reset preconditioner
    end
    largest = ifelse(isfinite(d.point), side === upper , d.point > 0)
    result = IterativeSolvers.lobpcg(d.lmap, largest, d.method.precond; kw...)
    λs, ϕs = result.λ, result.X
    isfinite(d.point) && (λs .= 1 ./ λs .+ d.point)
    precond && foreach(i -> (d.method.precond[i] = ϕs[i]), eachindex(d.method.precond))
    return Eigen(real(λs), ϕs)
end