text:
module HTTPC

using Compat
using LibCURL
using LibCURL.Mime_ext

import Base.convert, Base.show, Base.get, Base.trace

export init, cleanup, get, put, post, trace, delete, head, options
export RequestOptions, Response

def_rto = 0.0

##############################
# Struct definitions
##############################

type RequestOptions
    blocking::Bool
    query_params::Vector{Tuple}
    request_timeout::Float64
    callback::Union{Function,Bool}
    content_type::AbstractString
    headers::Vector{Tuple}
    ostream::Union{IO, AbstractString, Void}
    auto_content_type::Bool

    RequestOptions(; blocking=true, query_params=Array(Tuple,0), request_timeout=def_rto, callback=null_cb, content_type="", headers=Array(Tuple,0), ostream=nothing, auto_content_type=true) =
    new(blocking, query_params, request_timeout, callback, content_type, headers, ostream, auto_content_type)
end

type Response
    body
    headers :: Dict{AbstractString, Vector{AbstractString}}
    http_code
    total_time
    bytes_recd::Integer

    Response() = new(nothing, Dict{AbstractString, Vector{AbstractString}}(), 0, 0.0, 0)
end

function show(io::IO, o::Response)
    println(io, "HTTP Code   :", o.http_code)
    println(io, "RequestTime :", o.total_time)
    println(io, "Headers     :")
    for (k,vs) in o.headers
        for v in vs
            println(io, "    $k : $v")
        end
    end

    println(io, "Length of body : ", o.bytes_recd)
end

type ReadData
    typ::Symbol
    src::Any
    str::AbstractString
    offset::Csize_t
    sz::Csize_t

    ReadData() = new(:undefined, false, "", 0, 0)
end

type ConnContext
    curl::Ptr{CURL}
    url::AbstractString
    slist::Ptr{Void}
    rd::ReadData
    resp::Response
    options::RequestOptions
    close_ostream::Bool

    ConnContext(options::RequestOptions) = new(C_NULL, "", C_NULL, ReadData(), Response(), options, false)
end

immutable CURLMsg2
  msg::CURLMSG
  easy_handle::Ptr{CURL}
  data::Ptr{Any}
end

type MultiCtxt
    s::curl_socket_t    # Socket
    chk_read::Bool
    chk_write::Bool
    timeout::Float64

    MultiCtxt() = new(0,false,false,0.0)
end

##############################
# Callbacks
##############################

function write_cb(buff::Ptr{UInt8}, sz::Csize_t, n::Csize_t, p_ctxt::Ptr{Void})
#    println("@write_cb")
    ctxt = unsafe_pointer_to_objref(p_ctxt)
    nbytes = sz * n
    unsafe_write(ctxt.resp.body, buff, nbytes)
    ctxt.resp.bytes_recd = ctxt.resp.bytes_recd + nbytes

    nbytes::Csize_t
end

c_write_cb = cfunction(write_cb, Csize_t, (Ptr{UInt8}, Csize_t, Csize_t, Ptr{Void}))

function header_cb(buff::Ptr{UInt8}, sz::Csize_t, n::Csize_t, p_ctxt::Ptr{Void})
#    println("@header_cb")
    ctxt = unsafe_pointer_to_objref(p_ctxt)
    hdrlines = split(unsafe_string(buff, convert(Int, sz * n)), "\r\n")

#    println(hdrlines)
    for e in hdrlines
        m = match(r"^\s*([\w\-\_]+)\s*\:(.+)", e)
        if (m != nothing)
            k = strip(m.captures[1])
            v = strip(m.captures[2])
            if haskey(ctxt.resp.headers, k)
                push!(ctxt.resp.headers[k], v)
            else
                ctxt.resp.headers[k] = (AbstractString)[v]
            end
        end
    end
    (sz*n)::Csize_t
end

c_header_cb = cfunction(header_cb, Csize_t, (Ptr{UInt8}, Csize_t, Csize_t, Ptr{Void}))

function curl_read_cb(out::Ptr{Void}, s::Csize_t, n::Csize_t, p_ctxt::Ptr{Void})
#    println("@curl_read_cb")

    ctxt = unsafe_pointer_to_objref(p_ctxt)
    bavail::Csize_t = s * n
    breq::Csize_t = ctxt.rd.sz - ctxt.rd.offset
    b2copy = bavail > breq ? breq : bavail

    if (ctxt.rd.typ == :buffer)
        ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt64),
                out, convert(Ptr{UInt8}, pointer(ctxt.rd.str)) + ctxt.rd.offset, b2copy)
    elseif (ctxt.rd.typ == :io)
        b_read = read(ctxt.rd.src, UInt8, b2copy)
        ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, UInt64), out, b_read, b2copy)
    end
    ctxt.rd.offset = ctxt.rd.offset + b2copy

    r = convert(Csize_t, b2copy)
    r::Csize_t
end

c_curl_read_cb = cfunction(curl_read_cb, Csize_t, (Ptr{Void}, Csize_t, Csize_t, Ptr{Void}))

function curl_socket_cb(curl::Ptr{Void}, s::Cint, action::Cint, p_muctxt::Ptr{Void}, sctxt::Ptr{Void})
    if action != CURL_POLL_REMOVE
        muctxt = unsafe_pointer_to_objref(p_muctxt)

        muctxt.s = s
        muctxt.chk_read = false
        muctxt.chk_write = false

        if action == CURL_POLL_IN
            muctxt.chk_read = true

        elseif action == CURL_POLL_OUT
            muctxt.chk_write = true

        elseif action == CURL_POLL_INOUT
            muctxt.chk_read = true
            muctxt.chk_write = true
        end
    end

    # NOTE: Out-of-order socket fds cause problems in the case of HTTP redirects, hence ignoring CURL_POLL_REMOVE
    ret = convert(Cint, 0)
    ret::Cint
end

c_curl_socket_cb = cfunction(curl_socket_cb, Cint, (Ptr{Void}, Cint, Cint, Ptr{Void}, Ptr{Void}))

function curl_multi_timer_cb(curlm::Ptr{Void}, timeout_ms::Clong, p_muctxt::Ptr{Void})
    muctxt = unsafe_pointer_to_objref(p_muctxt)
    muctxt.timeout = timeout_ms / 1000.0

#    println("Requested timeout value : " * string(muctxt.timeout))

    ret = convert(Cint, 0)
    ret::Cint
end

c_curl_multi_timer_cb = cfunction(curl_multi_timer_cb, Cint, (Ptr{Void}, Clong, Ptr{Void}))

##############################
# Utility functions
##############################

macro ce_curl(f, args...)
    quote
        cc = CURLE_OK
        cc = $(esc(f))(ctxt.curl, $(args...))

        if(cc != CURLE_OK)
            error(string($f) * "() failed: " * unsafe_string(curl_easy_strerror(cc)))
        end
    end
end

macro ce_curlm(f, args...)
    quote
        cc = CURLM_OK
        cc = $(esc(f))(curlm, $(args...))

        if(cc != CURLM_OK)
            error(string($f) * "() failed: " * unsafe_string(curl_multi_strerror(cc)))
        end
    end
end

null_cb(curl) = return nothing

function set_opt_blocking(options::RequestOptions)
        o2 = RequestOptions()
        for n in filter(x -> !(x in [:ostream, :blocking]),fieldnames(o2))
            setfield!(o2, n, deepcopy(getfield(options, n)))
        end
        o2.blocking = true
        o2.ostream = options.ostream
        return o2
end

function get_ct_from_ext(filename)
    fparts = split(basename(filename), ".")
    if (length(fparts) > 1)
        if haskey(MimeExt, fparts[end]) return MimeExt[fparts[end]] end
    end
    return false
end

function setup_easy_handle(url, options::RequestOptions)
    ctxt = ConnContext(options)

    curl = curl_easy_init()
    if (curl == C_NULL) throw("curl_easy_init() failed") end

    ctxt.curl = curl

    @ce_curl curl_easy_setopt CURLOPT_FOLLOWLOCATION 1

    @ce_curl curl_easy_setopt CURLOPT_MAXREDIRS 5

    if length(options.query_params) > 0
        qp = urlencode_query_params(curl, options.query_params)
        url = url * "?" * qp
    end

    ctxt.url = url

    @ce_curl curl_easy_setopt CURLOPT_URL url
    @ce_curl curl_easy_setopt CURLOPT_WRITEFUNCTION c_write_cb

    p_ctxt = pointer_from_objref(ctxt)

    @ce_curl curl_easy_setopt CURLOPT_WRITEDATA p_ctxt

    @ce_curl curl_easy_setopt CURLOPT_HEADERFUNCTION c_header_cb
    @ce_curl curl_easy_setopt CURLOPT_HEADERDATA p_ctxt

    if options.content_type != ""
        ct = "Content-Type: " * options.content_type
        ctxt.slist = curl_slist_append(ctxt.slist, ct)
    else
        # Disable libCURL automatically setting the content type
        ctxt.slist = curl_slist_append(ctxt.slist, "Content-Type:")
    end

    for hdr in options.headers
        hdr_str = hdr[1] * ":" * hdr[2]
        ctxt.slist = curl_slist_append(ctxt.slist, hdr_str)
    end

    # Disabling the Expect header since some webservers don't handle this properly
    ctxt.slist = curl_slist_append(ctxt.slist, "Expect:")
    @ce_curl curl_easy_setopt CURLOPT_HTTPHEADER ctxt.slist

    if isa(options.ostream, AbstractString)
        ctxt.resp.body = open(options.ostream, "w+")
        ctxt.close_ostream = true
    elseif isa(options.ostream, IO)
        ctxt.resp.body = options.ostream
    else
        ctxt.resp.body = IOBuffer()
    end

    ctxt
end

function cleanup_easy_context(ctxt::Union{ConnContext,Bool})
    if isa(ctxt, ConnContext)
        if (ctxt.slist != C_NULL)
            curl_slist_free_all(ctxt.slist)
            ctxt.slist = C_NULL
        end

        if (ctxt.curl != C_NULL)
            curl_easy_cleanup(ctxt.curl)
            ctxt.curl = C_NULL
        end

        if ctxt.close_ostream
            close(ctxt.resp.body)
            ctxt.resp.body = nothing
            ctxt.close_ostream = false
        end
    end
end

function process_response(ctxt)
    http_code = Array(Clong,1)
    @ce_curl curl_easy_getinfo CURLINFO_RESPONSE_CODE http_code

    total_time = Array(Cdouble,1)
    @ce_curl curl_easy_getinfo CURLINFO_TOTAL_TIME total_time

    ctxt.resp.http_code = http_code[1]
    ctxt.resp.total_time = total_time[1]
end

# function blocking_get (url)
#     try
#         ctxt=nothing
#         ctxt = setup_easy_handle(url)
#         curl = ctxt.curl
#
#         @ce_curl curl_easy_perform
#
#         process_response(ctxt)
#
#         return ctxt.resp
#     finally
#         if isa(ctxt, ConnContext) && (ctxt.curl != 0)
#             curl_easy_cleanup(ctxt.curl)
#         end
#     end
# end

##############################
# Library initializations
##############################

init() = curl_global_init(CURL_GLOBAL_ALL)
cleanup() = curl_global_cleanup()

##############################
# GET
##############################

function get(url::AbstractString, options::RequestOptions=RequestOptions())
    if (options.blocking)
        ctxt = false
        try
            ctxt = setup_easy_handle(url, options)

            @ce_curl curl_easy_setopt CURLOPT_HTTPGET 1

            return exec_as_multi(ctxt)
        finally
            cleanup_easy_context(ctxt)
        end
    else
        return remotecall(get, myid(), url, set_opt_blocking(options))
    end
end

##############################
# POST & PUT
##############################

function post(url::AbstractString, data, options::RequestOptions=RequestOptions())
    if (options.blocking)
        return put_post(url, data, :post, options)
    else
        return remotecall(post, myid(), url, data, set_opt_blocking(options))
    end
end

function put(url::AbstractString, data, options::RequestOptions=RequestOptions())
    if (options.blocking)
        return put_post(url, data, :put, options)
    else
        return remotecall(put, myid(), url, data, set_opt_blocking(options))
    end
end

function put_post(url::AbstractString, data, putorpost::Symbol, options::RequestOptions)
    rd::ReadData = ReadData()

    if isa(data, AbstractString)
        rd.typ = :buffer
        rd.src = false
        rd.str = data
        rd.sz = length(data)

    elseif isa(data, Dict) || (isa(data, Vector) && issubtype(eltype(data), Tuple))
        arr_data = isa(data, Dict) ? Array{Tuple,1}(map((d) -> (d[1], d[2]), collect(data))) : data
        rd.str = urlencode_query_params(arr_data)  # Not very optimal since it creates another curl handle, but it is clean...
        rd.sz = length(rd.str)
        rd.typ = :buffer
        rd.src = arr_data
        if ((options.content_type == "") && (options.auto_content_type))
            options.content_type = "application/x-www-form-urlencoded"
        end

    elseif isa(data, IO)
        rd.typ = :io
        rd.src = data
        seekend(data)
        rd.sz = position(data)
        seekstart(data)
        if ((options.content_type == "") && (options.auto_content_type))
            options.content_type = "application/octet-stream"
        end

    elseif isa(data, Tuple)
        (typsym, filename) = data
        if (typsym != :file) error("Unsupported data datatype") end

        rd.typ = :io
        rd.src = open(filename)
        rd.sz = filesize(filename)

        try
            if ((options.content_type == "") && (options.auto_content_type))
                options.content_type = get_ct_from_ext(filename)
            end
            return _put_post(url, putorpost, options, rd)
        finally
            close(rd.src)
        end

    else
        error("Unsupported data datatype")
    end

    return _put_post(url, putorpost, options, rd)
end

function _put_post(url::AbstractString, putorpost::Symbol, options::RequestOptions, rd::ReadData)
    ctxt = false
    try
        ctxt = setup_easy_handle(url, options)
        ctxt.rd = rd

        if (putorpost == :post)
            @ce_curl curl_easy_setopt CURLOPT_POST 1
            @ce_curl curl_easy_setopt CURLOPT_POSTFIELDSIZE rd.sz
        elseif (putorpost == :put)
            @ce_curl curl_easy_setopt CURLOPT_UPLOAD 1
            @ce_curl curl_easy_setopt CURLOPT_INFILESIZE rd.sz
        end

        if (rd.typ == :io) || (putorpost == :put)
            p_ctxt = pointer_from_objref(ctxt)
            @ce_curl curl_easy_setopt CURLOPT_READDATA p_ctxt

            @ce_curl curl_easy_setopt CURLOPT_READFUNCTION c_curl_read_cb
        else
            ppostdata = pointer(rd.str)
            @ce_curl curl_easy_setopt CURLOPT_COPYPOSTFIELDS ppostdata
        end

        return exec_as_multi(ctxt)
    finally
        cleanup_easy_context(ctxt)
    end
end

##############################
# HEAD, DELETE and TRACE
##############################
function head(url::AbstractString, options::RequestOptions=RequestOptions())
    if (options.blocking)
        ctxt = false
        try
            ctxt = setup_easy_handle(url, options)

            @ce_curl curl_easy_setopt CURLOPT_NOBODY 1

            return exec_as_multi(ctxt)
        finally
            cleanup_easy_context(ctxt)
        end
    else
        return remotecall(head, myid(), url, set_opt_blocking(options))
    end

end

delete(url::AbstractString, options::RequestOptions=RequestOptions()) = custom(url, "DELETE", options)
trace(url::AbstractString, options::RequestOptions=RequestOptions()) = custom(url, "TRACE", options)
options(url::AbstractString, options::RequestOptions=RequestOptions()) = custom(url, "OPTIONS", options)

for f in (:get, :head, :delete, :trace, :options)
    @eval $(f)(url::AbstractString; kwargs...) = $(f)(url, RequestOptions(; kwargs...))
end

# put(url::AbstractString, data::AbstractString; kwargs...) = put(url, data, options=RequestOptions(; kwargs...))
# post(url::AbstractString, data::AbstractString; kwargs...) = post(url, data, options=RequestOptions(; kwargs...))

for f in (:put, :post)
    @eval $(f)(url::AbstractString, data::AbstractString; kwargs...) = $(f)(url, data, RequestOptions(; kwargs...))
end

function custom(url::AbstractString, verb::AbstractString, options::RequestOptions)
    if (options.blocking)
        ctxt = false
        try
            ctxt = setup_easy_handle(url, options)

            @ce_curl curl_easy_setopt CURLOPT_CUSTOMREQUEST verb

            return exec_as_multi(ctxt)
        finally
            cleanup_easy_context(ctxt)
        end
    else
        return remotecall(custom, myid(), url, verb, set_opt_blocking(options))
    end
end

##############################
# EXPORTED UTILS
##############################

function urlencode_query_params{T<:Tuple}(params::Vector{T})
    curl = curl_easy_init()
    if (curl == C_NULL) throw("curl_easy_init() failed") end

    querystr = urlencode_query_params(curl, params)

    curl_easy_cleanup(curl)

    return querystr
end
export urlencode_query_params

function urlencode_query_params{T<:Tuple}(curl, params::Vector{T})
    querystr = ""
    for x in params
        k,v = x
        if (v != "")
            ep = urlencode(curl, string(k)) * "=" * urlencode(curl, string(v))
        else
            ep = urlencode(curl, string(k))
        end

        if querystr == ""
            querystr = ep
        else
            querystr = querystr * "&" * ep
        end

    end
    return querystr
end

function urlencode(curl, s::AbstractString)
    b_arr = curl_easy_escape(curl, s, sizeof(s))
    esc_s = unsafe_string(b_arr)
    curl_free(b_arr)
    return esc_s
end

function urlencode(s::AbstractString)
    curl = curl_easy_init()
    if (curl == C_NULL) throw("curl_easy_init() failed") end

    esc_s = urlencode(curl, s)
    curl_easy_cleanup(curl)
    return esc_s

end

urlencode(s::SubString) = urlencode(unsafe_string(s))

export urlencode

function exec_as_multi(ctxt)
    curl = ctxt.curl
    curlm = curl_multi_init()

    if (curlm == C_NULL) error("Unable to initialize curl_multi_init()") end

    try
        if isa(ctxt.options.callback, Function) ctxt.options.callback(curl) end

        @ce_curlm curl_multi_add_handle curl

        n_active = Array(Cint,1)
        n_active[1] = 1

        no_to = 30 * 24 * 3600.0
        request_timeout = 0.001 + (ctxt.options.request_timeout == 0.0 ? no_to : ctxt.options.request_timeout)

        started_at = time()
        time_left = request_timeout

    # poll_fd is unreliable when multiple parallel fds are active, hence using curl_multi_perform

# START curl_multi_socket_action  mode

#         @ce_curlm curl_multi_setopt CURLMOPT_SOCKETFUNCTION c_curl_socket_cb
#         @ce_curlm curl_multi_setopt CURLMOPT_TIMERFUNCTION c_curl_multi_timer_cb
#
#         muctxt = MultiCtxt()
#         p_muctxt = pointer_from_objref(muctxt)
#
#         @ce_curlm curl_multi_setopt CURLMOPT_SOCKETDATA p_muctxt
#         @ce_curlm curl_multi_setopt CURLMOPT_TIMERDATA p_muctxt
#
#
#         @ce_curlm curl_multi_socket_action CURL_SOCKET_TIMEOUT 0 n_active
#
#         while (n_active[1] > 0) && (time_left > 0)
#             evt_got = 0
#             if (muctxt.chk_read || muctxt.chk_write)
#                 t1 = int64(time() * 1000)
#
#                 poll_to = min(muctxt.timeout < 0.0 ? no_to : muctxt.timeout, time_left)
#                 pfd_ret = poll_fd(RawFD(muctxt.s), poll_to, readable=muctxt.chk_read, writable=muctxt.chk_write)
#
#                 evt_got = (isreadable(pfd_ret) ? CURL_CSELECT_IN : 0) | (iswritable(pfd_ret) ? CURL_CSELECT_OUT : 0)
#             else
#                 break
#             end
#
#             if (evt_got == 0)
#                 @ce_curlm curl_multi_socket_action CURL_SOCKET_TIMEOUT 0 n_active
#             else
#                 @ce_curlm curl_multi_socket_action muctxt.s evt_got n_active
#             end
#
#             time_left = request_timeout - (time() - started_at)
#         end

# END curl_multi_socket_action  mode

# START curl_multi_perform  mode

        cmc = curl_multi_perform(curlm, n_active);
        while (n_active[1] > 0) &&  (time_left > 0)
            nb1 = ctxt.resp.bytes_recd
            cmc = curl_multi_perform(curlm, n_active);
            if(cmc != CURLM_OK) error("curl_multi_perform() failed: " * unsafe_string(curl_multi_strerror(cmc))) end

            nb2 = ctxt.resp.bytes_recd

            if (nb2 > nb1)
                yield() # Just yield to other tasks
            else
                sleep(0.005) # Just to prevent unnecessary CPU spinning
            end

            time_left = request_timeout - (time() - started_at)
        end

# END OF curl_multi_perform

        if (n_active[1] == 0)
            msgs_in_queue = Array(Cint,1)
            p_msg::Ptr{CURLMsg2} = curl_multi_info_read(curlm, msgs_in_queue)

            while (p_msg != C_NULL)
#                println("Messages left in Q : " * string(msgs_in_queue[1]))
                msg = unsafe_load(p_msg)

                if (msg.msg == CURLMSG_DONE)
                    ec = convert(Int, msg.data)
                    if (ec != CURLE_OK)
#                        println("Result of transfer: " * string(msg.data))
                        throw("Error executing request : " * unsafe_string(curl_easy_strerror(ec)))
                    else
                        process_response(ctxt)
                    end
                end

                p_msg = curl_multi_info_read(curlm, msgs_in_queue)
            end
        else
            error("request timed out")
        end

    finally
        curl_multi_remove_handle(curlm, curl)
        curl_multi_cleanup(curlm)
    end

    ctxt.resp
end

end


text:
@testset "ModelKit integration" begin
    @testset "total degree" begin
        @var x y
        f = [2 * x^2 + 3y + 5, x * y + 2 * y^2 + 2x + 3]
        @test_throws ArgumentError solve(f)
        res = solve(f; variable_ordering = [x, y])
        @test nnonsingular(res) == 4

        res = solve(System(f, [x, y]))
        @test nnonsingular(res) == 4

        # overdetermined
        g = [f; sum(f)]
        res = solve(System(g, [x, y]))
        @test nnonsingular(res) == 4

        @var x y z
        @test_throws ArgumentError problem_startsolutions(
            [x - 2y, y^2 + 3 * x^2, x^3 + y^3],
            variable_ordering = [x, y],
        )

        prob, starts = problem_startsolutions(
            [x^2 - 2 * z^2, y^2 + 3 * z^2, z^2 + x^2, y * z + x^2],
            variable_ordering = [x, y, z],
        )
        @test prob isa HC.OverdeterminedProblem{HC.ProjectiveTracking}
        @test starts isa HC.TotalDegreeSolutionIterator
        @test starts.degrees == [2, 2]
    end

    @testset "polyhedral" begin
        # converts into polynomial in the background
        @var x y
        result = solve(
            [2y + 3 * y^2 - x * y^3, x + 4 * x^2 - 2 * x^3 * y];
            variable_ordering = [x, y], seed = 32241, start_system = :polyhedral,
        )
        @test nsolutions(result) == 6

        result = solve(
            System([2y + 3 * y^2 - x * y^3, x + 4 * x^2 - 2 * x^3 * y], [x, y]),
            seed = 32241,
            start_system = :polyhedral,
        )
        @test nsolutions(result) == 6
    end

    @testset "parameter homotopy" begin
        @var x y z a b

        f = System([x^2 - a, y^2 - b], [x, y], [a, b])
        res = solve(
            f,
            [[1, -1], [1, 1], [-1, 1], [-1, -1]];
            start_parameters = [1, 1], target_parameters = [3, 4],
        )
        @test nnonsingular(res) == 4

        f = System([x^2 - z^2 * a, y^2 - z^2 * b], [x, y, z], [a, b])
        res = solve(
            f,
            [[1, -1], [1, 1], [-1, 1], [-1, -1]];
            start_parameters = [1, 1], target_parameters = [3, 4],
        )
        @test nnonsingular(res) == 4

        f = System([x^2 - a, y^2 - b], [x, y], [a, b])
        res = solve(
            f,
            [[1, -1], [1, 1], [-1, 1], [-1, -1]];
            start_parameters = [1, 1],
            target_parameters = [3, 4],
            γ₁ = cis(2.132im),
            γ₀ = cis(0.412im),
        )
        @test nnonsingular(res) == 4
    end

    @testset "Monodromy" begin
        @var x y
        f = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
        # define new variables u₁, u₂ and λ₁
        @var u[1:2] λ[1:1]
        # define the jacobian of F
        J = differentiate([f], [x, y])
        # J' defines the transpose of J
        C = [[x, y] - u - J' * λ; f]
        F = System(C, [x; y; λ], u)

        s₀ = [
            0.3785092033935262 + 0.20148856470590054im,
            0.008762115289894847 + 1.0025348274195496im,
            -0.08344446740253277 + 0.01612603621102832im,
        ]
        u₀ = [
            0.3373116901371924 + 0.12535355543080084im,
            0.27964298581644154 + 1.9612473908330865im,
        ]
        @test_throws ArgumentError monodromy_solve(C, s₀, u₀)
        @test_throws ArgumentError monodromy_solve(C, s₀, u₀; variable_ordering = [x; y; λ])

        res = monodromy_solve(
            C,
            s₀,
            u₀;
            variable_ordering = [x; y; λ],
            parameters = u,
            seed = 1312,
            target_solutions_count = 36,
        )
        @test nsolutions(res) == 36
        res = monodromy_solve(
            F,
            s₀,
            u₀;
            threading = false, target_solutions_count = 36, seed = 1312,
        )
        @test nsolutions(res) == 36

        # automatic start pair
        @var x y a b
        f = [x^2 - a, y^2 - b]
        res =
            monodromy_solve(f; parameters = [a, b], variable_ordering = [x, y], seed = 1312)
        @test nsolutions(res) == 4
        F = System(f, [x, y], [a, b])
        res = monodromy_solve(F; seed = 1312)
        @test nsolutions(res) == 4
    end
end

ext: jl


text:
module panel
import Term:
    split_lines,
    reshape_text,
    do_by_line,
    join_lines,
    truncate,
    textlen,
    fillin

import ..consoles: console_width, console_height
import ..measure: Measure
import ..renderables: AbstractRenderable, RenderablesUnion, Renderable, RenderableText
import ..segment: Segment
using ..box
import ..layout: pad, vstack, Padding
import ..style: apply_style

export Panel, TextBox

abstract type AbstractPanel <: AbstractRenderable end

# ---------------------------------------------------------------------------- #
#                                     PANEL                                    #
# ---------------------------------------------------------------------------- #
"""
    Panel

`Renderable` with a panel surrounding some content:
        ╭──────────╮
        │ my panel │
        ╰──────────╯
"""
mutable struct Panel <: AbstractPanel
    segments::Vector
    measure::Measure
    title::Union{Nothing,String}
    title_style::Union{String,Nothing}
    style::Union{String,Nothing}
end

"""
    Panel(
        content::RenderablesUnion;
        title::Union{Nothing,String} = nothing,
        title_style::Union{String,Nothing} = nothing,
        title_justify::Symbol = :left,
        subtitle::Union{String,Nothing} = nothing,
        subtitle_style::Union{String,Nothing} = nothing,
        subtitle_justify::Symbol = :left,
        style::Union{String,Nothing} = "default",
        box::Symbol = :ROUNDED,
        width::Int = 88,
        height::Union{Nothing,Int} = nothing,
        fit::Symbol=:nofit,
        justify = :left,
    )

`Panel` constructor to fit a panel to a piece of (renderable) content.

`title` can be used to specify a title to be addded to the top row and
`title_style` and `title_justify` set its appearance and position.
Same for `subtitle` but for the panel's bottom row.
`box` accepts one of the named `Term.box.Box` object and sets
which box should be used.

`width` and `height` are used to set the `Panel`'s size.
If `fit=true` is passed the width and height arguments are
ignored and the panel's box is fitted to the content's size
`justify` (:left, :center, :right) defines how the content should
be justified withing the panel.
`padding` specifies the ammount of padding between the box and the content.
"""
function Panel(
    content::Union{Nothing, RenderablesUnion};
    title::Union{Nothing,String} = nothing,
    title_style::Union{String,Nothing} = nothing,
    title_justify::Symbol = :left,
    subtitle::Union{String,Nothing} = nothing,
    subtitle_style::Union{String,Nothing} = nothing,
    subtitle_justify::Symbol = :left,
    style::Union{String,Nothing} = "default",
    box::Symbol = :ROUNDED,
    width::Int = 88,
    height::Union{Nothing,Int} = nothing,
    fit::Bool=false,
    justify::Symbol = :left,
    padding::Union{Vector, Padding, NTuple} = Padding(2, 2, 0, 0)

)
    box = eval(box)  # get box object from symbol

    # if content is text, ensure all lines have same width
    content = content isa AbstractString ? fillin(content) : content

    # get measure
    WIDTH = console_width(stdout)
    padding = padding isa Padding ? padding : Padding(padding...)
    Δw = padding.left + padding.right + 2
    Δh = padding.top + padding.bottom

    # define convenience function
    function resize_text(content, _width)
        if content isa AbstractString || content isa RenderableText
            content = RenderableText(content; width=_width)
            content_measure = content.measure
            return content, content.measure
        else
            return content, content_measure
        end
    end

    # get measure of panel's box and optionally resize text.
    if isnothing(content)
        if fit
            panel_measure = Measure(3, 2)
        else
            height = isnothing(height) ? 2 : height
            panel_measure = Measure(width, height)
        end
        # panel_measure = Measure(2, 2)
        content = ""
        content_measure = Measure(0, 0)
    else
        content_measure = Measure(content)
        if fit
            # if content width too large, resize content if its text
            if content_measure.w > WIDTH - Δw
                width = WIDTH-Δw
                content, content_measure = resize_text(content, width-Δw+2)
                panel_measure = Measure(width+2, content_measure.h+Δh+2)
            else
                panel_measure = Measure(content_measure.w+Δw, content_measure.h+Δh+2)
            end
        end

        if !fit
            # check that the content fits within the given width
            if content isa AbstractString || content isa RenderableText
                width = min(width, WIDTH)
                if content_measure.w > width-Δw
                    content, content_measure = resize_text(content, width-Δw)
                end
            else
                # if width too small for content, try to enlarge
                width = width < content_measure.w+Δw ?  min(content_measure.w+Δw, WIDTH-Δw) : width
            end

            # get target height
            _h = content_measure.h + Δh + 2
            height = isnothing(height) ? _h  : max(height, _h)
            panel_measure = Measure(width, height)
        end
    end

    # style stuff
    title_style = isnothing(title_style) ? style : title_style
    σ(s) = "[$style]$s[/$style]"  # applies the main style markup to a string to make a segment

    # create segments
    segments::Vector{Segment} = []

    # create top/bottom rows with titles
    top = get_title_row(
        :top,
        box,
        title;
        width = panel_measure.w-2,
        style = style,
        title_style = title_style,
        justify = title_justify,
    )

    bottom = get_title_row(
        :bottom,
        box,
        subtitle;
        width = panel_measure.w-2,
        style = style,
        title_style = subtitle_style,
        justify = subtitle_justify,
    )

    # add a panel row for each content row
    push!(segments, top)

    # add padding lines at top
    left, right = σ(string(box.mid.left)), σ(string(box.mid.right))

    addempty() = push!(segments, Segment(left * " "^(panel_measure.w-2) * right))
    for i in 1:padding.top
        addempty()
    end

    content_lines = split_lines(content)

    for n in 1:(content_measure.h)
        # apply style and pad
        line = pad(apply_style(content_lines[n]), panel_measure.w-Δw, justify)
        line = pad(line, padding.left, padding.right)

        # make line
        segment = Segment(left * line * right)

        push!(segments, segment)
    end

    # add empty lines to ensure target height is reached
    if content_measure.h < panel_measure.h - 2 - Δh
        for i in 1:(panel_measure.h-2 - content_measure.h)
            addempty()
        end
    end

    # add padding at the bottom
    for i in 1:padding.bottom
        addempty()
    end

    push!(segments, bottom)

    return Panel(
        segments, panel_measure, title, title_style, style
    )
end

"""
    Panel(renderables; kwargs...)

`Panel` constructor for creating a panel out of multiple renderables at once.
"""
function Panel(content, renderables...; kwargs...)
    content = vstack(content, Renderable.(renderables)...)
    return Panel(content; kwargs...)
end

Panel(content::AbstractVector; kwargs...) = Panel(content...; kwargs...)
Panel(; kwargs...) = Panel(nothing; kwargs...)

# ---------------------------------------------------------------------------- #
#                                    TextBox                                   #
# ---------------------------------------------------------------------------- #

"""
    TextBox

Creates a `Panel` and fits input text to it.
The pannel is hidden so that the result is just a text box.
"""
mutable struct TextBox <: AbstractPanel
    segments::Vector
    measure::Measure
end

"""
    TextBox(
        text::Union{Vector, AbstractString};
        width::Union{Nothing, Int}=nothing,
        title::Union{Nothing, String}=nothing,
        title_style::Union{String, Nothing}="default",
        title_justify::Symbol=:left,
        subtitle::Union{String, Nothing}=nothing,
        subtitle_style::Union{String, Nothing}="default",
        subtitle_justify::Symbol=:left,
        justify::Symbol=:left,
        fit::Symbol=:fit,
        )

Creates an hidden `Panel` with `text` in it.

If a `width` is passed, the input `text` is reshaped to have
that size, unless `fit=:truncate` in which case it's cut to size.
If no `width` is passed and `fit=true` the `TextBox`'s size
matches the size of the input `text`.
Other arguments behave like `Panel`.

See also [`Panel`](@ref).
"""
function TextBox(
    text::Union{AbstractString};
    width::Int = 88,
    title::Union{Nothing,String} = nothing,
    title_style::Union{String,Nothing} = "default",
    title_justify::Symbol = :left,
    subtitle::Union{String,Nothing} = nothing,
    subtitle_style::Union{String,Nothing} = "default",
    subtitle_justify::Symbol = :left,
    justify::Symbol = :left,
    fit::Symbol = :nofit,
)

    # fit text or get width
    if fit==:fit
        # the box's size depends on the text's size
        width = Measure(text).w + 4

        # too large, fit to console
        if width > console_width(stdout)
            width=console_width(stdout)
            fit = true
        end
    else
        width = width > console_width(stdout) ? console_width(stdout) - 4 : width
    end

    # truncate or reshape text
    if fit == :truncate
        # truncate the text to fit the given width
        text = do_by_line(ln -> truncate(ln, width - 7), text)
    else
        text = do_by_line((ln) -> reshape_text(ln, width - 4), text)
    end

    # create panel with text inside
    panel = Panel(
        text;
        style = "hidden",
        title = title,
        title_style = title_style,
        title_justify = title_justify,
        subtitle = subtitle,
        subtitle_style = subtitle_style,
        subtitle_justify = subtitle_justify,
        justify = justify,
        width = width,
        fit=false
    )

    return TextBox(panel.segments, Measure(panel.measure.w, panel.measure.h))
end

TextBox(texts...; kwargs...) = TextBox(join_lines(texts); kwargs...)

end


text:
using Makie, GLVisualize
using GeometryTypes, Colors
using Bridge: _viridis

skippoints = 1
extra = 50
d2 = length(valtype(V))
scene = Scene(resolution = (400, 400))

sphere2 = Sphere(Point2f0(0), 1.0f0)
lines((@view Xtrue.yy[1:skippoints:end]), linewidth = 0.2, color = :darkblue)
#lines((@view Xtrue.yy[1:skippoints:end]), linewidth = 0.2, color = [:darkblue,:darkred][1+(map(x->dot(R{3}(1,1,0),x), Xtrue.yy[1:skippoints:end]) .> 0) ])

if d2 == 3
    viri = _viridis[round.(Int,range(1, stop=256, length=length(V)))]
    scatter(V.yy, marker=sphere2, markersize = 0.2,  strokewidth = 0.02, strokecolor = :white, color = map(x->RGBA(Float32.(x)..., 0.9f0), viri))
elseif d2 == 2
 #   viri = _viridis[round.(Int,linspace(1,256, 2*length(V)))]
  #  linesegment(map(v->(Point3f0(-10, v...), Point3f0(10, v...)), V.yy), color = map(x->RGBA(Float32.(x)..., 0.9f0), viri))
    viri = _viridis[round.(Int,range(1, stop=256, length=length(V)))]
    scatter(map(v->Point3f0(0, v...), V.yy), marker=sphere2, markersize = 0.2,  strokewidth = 0.02, strokecolor = :white, color = map(x->RGBA(Float32.(x)..., 0.9f0), viri))
end
#for pt in Pt
#    lines((@view pt.Y.yy[1:skippoints:end]), linewidth = 0.2, color = :darkred)
#end

Nu = vcat([Point3f0.(Pᵒ[i].V[1:end-1]) for i in 1:m]...)
lines((@view Nu[1:skippoints:end]), linewidth = 0.2, color = :darkred)

lines((@view Xmean[1:skippoints:end]), linewidth = 0.2, color = :darkgreen)

visualize_uncertainty(scene, (Xmean, Xrot, Xscal), extra*skippoints; color = RGBA(0.2f0, 0.4f0, 0.9f0, 0.1f0))

sphere = Sphere(Point3f0(0,0,0), 1.0f0)

if d2 == 3
    axis(map(x -> UnitRange(round.(Int,x)...), extrema(Bridge.mat(V.yy), 2))...)
else
    #axis(map(x -> UnitRange(round.(Int,x)...), extrema(Bridge.mat(Xtrue.yy), 2))...)
    ma = maximum(Bridge.mat(X2.yy), 2)
    axis(map(x -> UnitRange(0, round.(Int,x)), ma)..., showaxis = (true,true,true), showgrid = (false,false,false))

    for (x, c) in zip(ℝ{3}[(ma[1], 0, 0), (0, ma[2], 0),(0, 0, ma[3]), (0,0,0)], ['x', 'y', 'z', '•'])
        scatter([x], marker = c, markersize = 1., color = :black)
    end
end

scatter(collect(Models.foci(P)), marker=Sphere(Point2f0(0), 1.0f0), markersize = 1.,  strokewidth = 0.1, strokecolor = :white, color = :green)

center!(scene)


text:
##############################################################################
#
# test_regularisation_functions.jl
#
# Test CVortex's regularisation methods
#
# Copyright 2020 HJA Bird
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.
###############################################################################


text:
using Printf
using Statistics
using Oceananigans
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom, mask_immersed_field!

@inline bottom_drag_func(x, y, t, u, Cᴰ) = - Cᴰ * u^2

function hilly_simulation(;
                          Nx = 64,
                          Nz = Nx,
                          ϵ = 0.1,
                          Re = 1e4,
                          N² = 1e-2,
                          bottom_drag = false,
                          stop_time = 1,
                          save_interval = 0.1,
                          architecture = CPU(),
                          name = "flow_over_hills")

    underlying_grid = RectilinearGrid(architecture,
                                      size = (Nx, Nz),
                                      halo = (3, 3),
                                      x = (0, 2π),
                                      z = (0, 2π),
                                      topology = (Periodic, Flat, Bounded))

    if ϵ > 0
        x, y, z = nodes((Center, Center, Center), underlying_grid, reshape=true)
        hills = @. ϵ * (1 + sin(x)) / 2
        grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(hills))
    else # no hills
        grid = underlying_grid
    end

    closure = isfinite(Re) ? ScalarDiffusivity(ν=1/Re, κ=1/Re) : nothing

    if bottom_drag
        Δz = 2π / Nz
        z₀ = 1e-4
        κ = 0.4
        Cᴰ = (κ / log(Δz / 2z₀))^2
        u_bottom_bc = FluxBoundaryCondition(bottom_drag_func, field_dependencies=:u, parameters=Cᴰ)
        u_bcs = FieldBoundaryConditions(bottom=u_bottom_bc)
        boundary_conditions = (; u = u_bcs)
        @info string("Using a bottom drag with coefficient ", Cᴰ)
    else
        boundary_conditions = NamedTuple()
    end

    model = NonhydrostaticModel(; grid, closure, boundary_conditions,
                                advection = WENO5(),
                                timestepper = :RungeKutta3,
                                tracers = :b,
                                buoyancy = BuoyancyTracer())

    bᵢ(x, y, z) = N² * z + 1e-9 * rand()
    set!(model, b=bᵢ, u=1)

    Δx = 2π / Nx
    Δt = 0.1 * Δx
    simulation = Simulation(model; Δt, stop_time)

    u, v, w = model.velocities
    Uᵢ = mean(u)

    wall_clock = Ref(time_ns())

    function progress(sim)
        δU = mean(u) / Uᵢ
        elapsed = 1e-9 * (time_ns() - wall_clock[])
        @info @sprintf("Iter: %d, time: %.2e, δU: %.2e, wall time: %s",
                       iteration(sim), time(sim), δU, prettytime(elapsed))
        wall_clock[] = time_ns()
        return nothing
    end

    simulation.callbacks[:progress] = Callback(progress, IterationInterval(100))

    U = Average(u, dims=(1, 2, 3))
    ξ = ∂z(u) - ∂x(w)

    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, merge(model.velocities, model.tracers, (; ξ, U)),
                         schedule = TimeInterval(save_interval),
                         prefix = name,
                         force = true)

    @info "Make a simulation:"
    @show simulation

    return simulation
end

function momentum_time_series(filepath)
    U = FieldTimeSeries(filepath, "U")
    t = U.times
    δU = [U[1, 1, 1, n] / U[1, 1, 1, 1] for n=1:length(t)]
    return δU, t
end

Nx = 256
stop_time = 100.0
reference_name = "bottom_drag_reference"
#=
reference_sim = hilly_simulation(; stop_time, Nx, ϵ=0, name=reference_name, bottom_drag=true)
run!(reference_sim)
δU_ref, t_ref = momentum_time_series(name * ".jld2")
=#

experiments = []
for ϵ = [0.02, 0.05, 0.1]
    name = string("flow_over_hills_height_", ϵ)
    push!(experiments, (; ϵ, name))
    simulation = hilly_simulation(; stop_time, Nx, ϵ, name)
    run!(simulation)
end

using GLMakie

fig = Figure()
ax = Axis(fig[1, 1])
lines!(ax, t_ref, δU_ref, label="Reference with bottom drag")

for experiment in experiments
    name = experiment.name
    ϵ = experiment.ϵ
    δU, t = momentum_time_series(name * ".jld2")
    lines!(ax, t, δU, label=string("ϵ = ", ϵ))
end

axislegend(ax)

display(fig)

#=
# Animate vorticity if you like
filepath = experiments[end].name * ".jld2"
ξ = FieldTimeSeries(filepath, "ξ")
Nt = length(ξ.times)

fig = Figure()
slider = Slider(fig[2, 1], range=1:Nt, startvalue=1)
n = slider.value

title = @lift @sprintf("Vorticity in flow over hills at t = %.2e", ξ.times[$n])
ax = Axis(fig[1, 1]; title)
ξn = @lift interior(ξ[$n], :, 1, :)

#=
masked_ξn = @lift begin
    ξn = ξ[$n]
    mask_immersed_field!(ξn, NaN)
    interior(ξn, :, 1, :)
end
=#

heatmap!(ax, ξn)

display(fig)
=#


text:
# Autogenerated wrapper script for MPSolve_jll for powerpc64le-linux-gnu-cxx11
export libmps

using GMP_jll
JLLWrappers.@generate_wrapper_header("MPSolve")
JLLWrappers.@declare_library_product(libmps, "libmps.so.3")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libmps,
        "lib/libmps.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Revise
using NLSolvers
using LinearAlgebra
using SparseDiffTools
using SparseArrays
using IterativeSolvers
using DoubleFloats
using ForwardDiff
using Test

@testset "nlequations interface" begin
#
# Have Anderosn use FixedPointProblem
#
#

function theta(x)
   if x[1] > 0
       return atan(x[2] / x[1]) / (2.0 * pi)
   else
       return (pi + atan(x[2] / x[1])) / (2.0 * pi)
   end
end

function F_fletcher_powell!(Fx, x)
    if ( x[1]^2 + x[2]^2 == 0 )
        dtdx1 = 0;
        dtdx2 = 0;
    else
        dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
        dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
    end
    Fx[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
    200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
    Fx[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
    200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
    Fx[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3];
    Fx
end

function F_jacobian_fletcher_powell!(Fx, Jx, x)
    ForwardDiff.jacobian!(Jx, F_fletcher_powell!, Fx, x)
    Fx, Jx
end

jv = JacVec(F_fletcher_powell!, rand(3); autodiff=false)
function jvop(x)
    jv.u .= x
    jv
end
prob_obj = NLSolvers.NEqObjective(F_fletcher_powell!, nothing, F_jacobian_fletcher_powell!, jvop)

prob = NEqProblem(prob_obj)

x0 = [-1.0, 0.0, 0.0]
res = solve(prob, x0, LineSearch(Newton(), Backtracking()))
@test norm(res.info.best_residual, Inf) < 1e-12

#x0 = [-1.0, 0.0, 0.0]
#res = solve(prob, x0, Anderson(), NEqOptions())
#@test norm(res.info.best_residual, Inf) < 1e-12

x0 = [-1.0, 0.0, 0.0]
res = solve(prob, x0, TrustRegion(Newton()), NEqOptions())
@test norm(res.info.best_residual, Inf) < 1e-12

#x0 = [-1.0, 0.0, 0.0]
#res = solve(prob, x0, TrustRegion(DBFGS()), NEqOptions())
#@test norm(res.info.best_residual, Inf) < 1e-12

#x0 = [-1.0, 0.0, 0.0]
#res = solve(prob, x0, TrustRegion(BFGS()), NEqOptions())
#@test norm(res.info.best_residual, Inf) < 1e-12

x0 = [-1.0, 0.0, 0.0]
res = solve(prob, x0, TrustRegion(SR1()), NEqOptions())
@test norm(res.info.best_residual, Inf) < 1e-9

x0 = [-1.0, 0.0, 0.0]
state = (z=copy(x0), d=copy(x0), Fx=copy(x0), Jx=zeros(3,3))
res = solve(prob, x0, LineSearch(Newton(), Backtracking()), NEqOptions(), state)
@test norm(res.info.best_residual, Inf) < 1e-12

#x0 = [-1.0, 0.0, 0.0]
#res = solve(prob, x0, InexactNewton(FixedForceTerm(1e-3), 1e-3, 300), NEqOptions())
#@test norm(res.info.best_residual, Inf) < 1e-12

#x0 = [-1.0, 0.0, 0.0]
#@show solve(prob, x0, InexactNewton(EisenstatWalkerA(), 1e-8, 300), NEqOptions())
#@test norm(res.info.best_residual, Inf) < 1e-12

#x0 = [-1.0, 0.0, 0.0]
#@show solve(prob, x0, InexactNewton(EisenstatWalkerB(), 1e-8, 300), NEqOptions())
#@test norm(res.info.best_residual, Inf) < 1e-12

x0 = [-1.0, 0.0, 0.0]
res = solve(prob, x0, NLSolvers.DFSANE(), NLSolvers.NEqOptions(maxiter=1000))
@test norm(res.info.best_residual, Inf) < 1e-0

function dfsane_exponential(Fx, x)
    Fx[1] = exp(x[1]-1)-1
    for i = 2:length(Fx)
  	    Fx[i] = i*(exp(x[i]-1)-x[i])
    end
    Fx
end

function FJ_dfsane_exponential!(Fx, Jx, x)
  ForwardDiff.jacobian!(Jx, dfsane_exponential, Fx, x)
  Fx, Jx
end
prob_obj = NLSolvers.NEqObjective(dfsane_exponential, nothing, FJ_dfsane_exponential!, nothing)

n = 5000
x0 = fill(n/(n-1), n)
res = solve(NEqProblem(prob_obj), x0, NLSolvers.DFSANE(), NLSolvers.NEqOptions(maxiter=1000))
@test norm(res.info.best_residual, Inf)<1e-5

n = 1000
x0 = fill(n/(n-1), n)
res = solve(NEqProblem(prob_obj), x0, LineSearch(NLSolvers.Newton()), NLSolvers.NEqOptions(f_abstol=1e-6)) #converges well but doing this for time
@test norm(res.info.best_residual, Inf)<1e-6

n = 5000
x0 = fill(n/(n-1), n)
res = solve(NEqProblem(prob_obj), x0, TrustRegion(NLSolvers.Newton()), NLSolvers.NEqOptions(f_abstol=1e-6)) #converges well but doing this for time
@test norm(res.info.best_residual, Inf)<1e-6

#n = 1000
#x0 = fill(n/(n-1), n)
#@show solve(NEqProblem(prob_obj), x0, TrustRegion(NLSolvers.BFGS()), NLSolvers.NEqOptions())

function dfsane_exponential2(Fx, x)
    Fx[1] = exp(x[1])-1
    for i = 2:length(Fx)
        Fx[i] = i/10*(exp(x[i])+x[i-1]-1)
    end
    Fx
end
dfsane_prob2 = NLSolvers.NEqObjective(dfsane_exponential2, nothing, nothing,nothing)
n = 500
x0 = fill(1/(n^2), n)
res = solve(NEqProblem(dfsane_prob2), x0, NLSolvers.DFSANE(), NLSolvers.NEqOptions())
@test norm(res.info.best_residual, Inf) < 1e-5
res = solve(NEqProblem(dfsane_prob2), Double64.(x0), NLSolvers.DFSANE(), NLSolvers.NEqOptions())
@test norm(res.info.best_residual, Inf) < 1e-5
n = 2000
x0 = fill(1/(n^2), n)
res = solve(NEqProblem(dfsane_prob2), x0, NLSolvers.DFSANE(), NLSolvers.NEqOptions())
@test norm(res.info.best_residual, Inf) < 1e-5

# Take from problems
prob = NEqProblem(NLE_PROBS["rosenbrock"]["array"]["mutating"])
x0 = copy(NLE_PROBS["rosenbrock"]["array"]["x0"])
X0 = [[-0.8, 1.0], [-0.7, 1.0], [-0.7, 1.2], [-0.5, 0.7], [0.5, 0.4], [-1.2, 1.0]]
for x0 in X0
    res = solve(prob, x0, TrustRegion(Newton(), Dogleg()), NEqOptions())
    @test norm(res.info.best_residual, Inf) < 1e-12
    res = solve(prob, x0, TrustRegion(Newton(), NWI()), NEqOptions())
    @test norm(res.info.best_residual, Inf) < 1e-12
    res = solve(prob, x0, TrustRegion(Newton(), NTR()), NEqOptions())
    @test norm(res.info.best_residual, Inf) < 1e-12
    res = solve(prob, x0, LineSearch(Newton(), Static(0.6)), NEqOptions())
    @test norm(res.info.best_residual, Inf) < 1e-8
    #res = solve(prob, copy(initial), DFSANE(), NEqOptions())
    #res = solve(prob, copy(initial), InexactNewton(), NEqOptions())
end

@testset "fixedpoints" begin
# function G(x, Gx)
#     V1 = [1.0, 0.0] .+ 0.99*[0.1 0.9; 0.5 0.5]*x
#     V2 = [0.0, 2.0] .+ 0.99*[0.5 0.5; 1.0 0.0]*x
#     K = max(maximum(V1), maximum(V2))
#     Gx .= K .+ log.(exp.(V1 .- K) .+ exp.(V2 .- K))
#   end
#   fp1 = NLSolvers.fixedpoint!(G, zeros(2), Anderson(10000000,1, nothing,nothing))
#   fp2 = NLSolvers.fixedpoint!(G, zeros(2), Anderson(2, 2, 0.3, 1e2))
#   fp3 = NLSolvers.fixedpoint!(G, zeros(2), Anderson(2, 2, 0.01, 1e2))
#   @test all(fp1.x .≈ fp2.x)
#   @test all(fp1.x .≈ fp3.x)

#   fp1 = NLSolvers.solve(NEqProblem((x, F)->G(x, F).-x), zeros(2), Anderson(10000000,1, nothing,nothing))
#   fp2 = NLSolvers.solve(NEqProblem((x, F)->G(x, F).-x), zeros(2), Anderson(2, 2, 0.3, 1e2))
#   fp3 = NLSolvers.solve(NEqProblem((x, F)->G(x, F).-x), zeros(2), Anderson(2, 2, 0.01, 1e2))
#   @test all(fp1.x .≈ fp2.x)
#   @test all(fp1.x .≈ fp3.x)
end

@testset "complementarity" begin
# using NLsolve

# M = [0  0 -1 -1 ;
#      0  0  1 -2 ;
#      1 -1  2 -2 ;
#      1  2 -2  4 ]

# q = [2; 2; -2; -6]

# function f!(x, fvec)
#     fvec = M * x + q
# end

# r = mcpsolve(f!, [0., 0., 0., 0.], [Inf, Inf, Inf, Inf],
#              [1.25, 0., 0., 0.5], reformulation = :smooth, autodiff = true)

# x = r.zero  # [1.25, 0.0, 0.0, 0.5]
# @show dot( M*x + q, x )  # 0.5

# sol = [2.8, 0.0, 0.8, 1.2]
# @show dot( M*sol + q, sol )  # 0.0

end
end

@testset "lsqfit" begin
# boxbod_f(x, b) = b[1]*(1-exp(-b[2]*x))

# ydata = [109, 149, 149, 191, 213, 224]
# xdata = [1, 2, 3, 5, 7, 10]

# start1 = [1.0, 1.0]
# start2 = [100.0, 0.75]

# using Plots
# using NLSolvers

# nd = NonDiffed(t->sum(abs2, boxbod_f.(xdata, Ref(t)).-ydata))
# bounds = (fill(0.0, 2), fill(500.0, 2))
# problem = MinProblem(; obj=nd, bounds=bounds)
# @show minimize!(problem, zeros(2), APSO(), OptimizationOptions())
# @show minimize!(nd, [200.0, 1], NelderMead(), OptimizationOptions())

# @show minimize!(nd, minimize!(nd, [200.0, 10.0], NelderMead(), OptimizationOptions()).info.minimizer, NelderMead(), OptimizationOptions())
# @show minimize!(nd, minimize!(nd, [200.0, 5.0], NelderMead(), OptimizationOptions()).info.minimizer, NelderMead(), OptimizationOptions())

# function F(b, F, J=nothing)
#   @. F = b[1]*(1 - exp(-b[2]*xdata)) - ydata
#   if !isa(J, Nothing)
#   	@. J[:, 1] = 1 - exp(-b[2]*xdata)
#   	@. @views J[:, 2] = xdata*b[1]*(J[:, 1]-1)
#     return F, J
#   end
#   F
# end

# Fc = zeros(6)
# minimize!(lsqwrap, [100.0, 1.0], NelderMead(), OptimizationOptions())

# lsqwrap = NLSolvers.LsqWrapper(F, zeros(6), zeros(6,2))
# minimize!(MinProblem(;obj=lsqwrap,bounds=([0.0,0.0],[250.0,2.0])), [100.0, 1.0], APSO(), OptimizationOptions())

# function F(F, b)
#   @. F = b[1]*(1 - exp(-b[2]*xdata)) - ydata

#   F
# end
# function F(J, F, b)
#   @. F = b[1]*(1 - exp(-b[2]*xdata)) - ydata
#   if !isa(J, Nothing)
#   	@. J[:, 1] = 1 - exp(-b[2]*xdata)
#   	@. @views J[:, 2] = xdata*b[1]*(J[:, 1]-1)
#     return F, J
#   end
#   F
# end
# OnceDiffed(F)(rand(2), rand(6), rand(6,2)

# Fc = zeros(6)

# lsqwrap = NLSolvers.LsqWrapper(OnceDiffed(F), zeros(6), zeros(6,2))
# minimize!(lsqwrap, [100.0, 1.0], LineSearch(LBFGS()), OptimizationOptions())

# #using Plots
# #theme(:ggplot2)
# #gr(size=(500,500))
# #X = range(000.0, 350.0; length=420)
# #Y = range(0.0, 3.00; length=420)
# #contour(X, Y, (x, y)->nd([x, y]);
# #       fill=true,
# #       c=:turbid,levels=200, ls=:dash,
# #       xlims=(minimum(X), maximum(X)),
# #       ylims=(minimum(Y), maximum(Y)),
# #       colorbar=true)

end


text:
using MolecularConformation,PyPlot,BenchmarkTools

BenchmarkTools.DEFAULT_PARAMETERS.samples = 1000

"""
This function was created to make tests in order to perform our algorithms. Essentially,
to execute this routine is necessary to type in REPL

julia> perform(d)

where d is an integer value defined by the vector [3,4,5,10,50,100,200,300,400,500,600,700,800,900,1000,2000,3000]. Warning: some files need to be download.

"""
# ndiag is defined by [3,4,5,10,50,100,200,300,400,500,600,700,800,900,1000,2000,3000]
function perform(ndiag,allsolutions=false,LDE=false)
	#initialization
	opt_classic = ConformationSetup(0.000001,classicBP,allsolutions,LDE)
	opt_quaternion = ConformationSetup(0.000001,quaternionBP,allsolutions,LDE)
	data = preprocessing("toyinstance.nmr")
	sol = conformation(data,opt_quaternion)
	sol = conformation(data,opt_classic)
#	sol = conformation(data,opt_classicOpt)
	#
	folders = ["10","100","200","300","400","500","600","700","800","900","1000","2000","3000"]
	prob = parse.(Int,folders)
	tquat = zeros(length(folders),4)
	tclass = zeros(length(folders),4)
#	tclassopt = zeros(length(folders),4)
	k = 1
	c = 10.0^(-9)
	for foldername in folders
		cd(foldername)

		psize = parse(Int,foldername)
		if ndiag > psize
			data = preprocessing(string(foldername,"-$(psize).nmr"))
		else
			data = preprocessing(string(foldername,"-$(ndiag).nmr"))
		end
		bch  = @benchmark conformation($(data),$(opt_classic))
		tclass[k,1] = minimum(bch).time*c
		tclass[k,2] = median(bch).time*c
		tclass[k,3] = mean(bch).time*c
		tclass[k,4] = maximum(bch).time*c

#		bch  = @benchmark conformation($(data),$(opt_classicOpt))
#		tclassopt[k,1] = minimum(bch).time*c
#		tclassopt[k,2] = median(bch).time*c
#		tclassopt[k,3] = mean(bch).time*c
#		tclassopt[k,4] = maximum(bch).time*c

		bch = @benchmark conformation($(data),$(opt_quaternion))
		tquat[k,1] = minimum(bch).time*c
		tquat[k,2] = median(bch).time*c
		tquat[k,3] = mean(bch).time*c
		tquat[k,4] = maximum(bch).time*c

		k += 1
		cd("..")
	end
	#####
	improv(q,c) = (c/q -1.0)*100
	# average improvement in minimal time
	improv_min = sum(improv.(tquat[:,1],tclass[:,1]))/length(tquat[:,1])
	# average improvement in median time
	improv_med = sum(improv.(tquat[:,2],tclass[:,2]))/length(tquat[:,2])
	# average improvement in mean time
	improv_mean = sum(improv.(tquat[:,3],tclass[:,3]))/length(tquat[:,3])
	# average improvement in maximal time
	improv_max = sum(improv.(tquat[:,4],tclass[:,4]))/length(tquat[:,4])
	io = open("improvs_$(ndiag).txt", "w");
	write(io, "minimum -> $(improv_min) \n");
	write(io, "median -> $(improv_med) \n");
	write(io, "mean -> $(improv_mean) \n");
	write(io, "maximum -> $(improv_max) \n");
	close(io);
	#####
	j=1
	for information in ["minimum","median","mean","maximum"]
		ioff()
		fig = figure()
	    	xlabel("size of problems")
		ylabel("Average of $(information) processing time (s)")
    	#grid("on")
    	#title("")
		plot(prob,tquat[:,j],"--k",label="QuaternionBP");
#		plot(prob,tclassopt[:,j],"-.k",label="ClassicBPOpt");
		plot(prob,tclass[:,j],"-k",label="ClassicBP");
	    	legend(loc="upper right",fancybox="false");
		savefig("perf_$(information)_$(ndiag).png");
		close(fig)
		j+=1
	end
end
"""
	runperf

This function is used to run several tests

# Example
```
julia-repl
julia> runperf()

```
some graphs will be saved in current folder.
"""
function runperf()
	for ndiag in [3,10,100,200,400,800,1000]
		perform(ndiag)
	end
end


text:
# It's for libnode
using BinaryBuilder, Pkg

name = "libnode"
version = v"14.17.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://nodejs.org/dist/v14.17.3/node-v14.17.3.tar.gz", "dcbd156506ee79ee48439257626ca0a6db3d7eab8cb0208db6979125ae7d3a95"),
]

host_platform = HostPlatform()

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd node-*

export CC_host=$HOSTCC
export CXX_host=$HOSTCXX
patch ./node/src/inspector_agent.cc < ./patches/inspector_agent.cc.patch
cd node

if [[ $target == $MACHTYPE ]]
then
    ./configure --prefix=${prefix} --shared --no-cross-compiling
else
    DEST_CPU=x86_64
    if [[ $target == *i686* ]]; then DEST_CPU=x86; fi
    if [[ $target == *aarch64* ]]; then DEST_CPU=arm64; fi
    if [[ $target == *powerpc64* ]]; then DEST_CPU=ppc64; fi
    DEST_OS=linux
    if [[ $target == *apple* ]]; then DEST_OS=mac; fi
    if [[ $target == *freebsd* ]]; then DEST_OS=freebsd; fi
    if [[ $target == *w64* ]]; then DEST_OS=win; fi
    ./configure --prefix=${prefix} --shared --without-intl --cross-compiling --dest-cpu=$DEST_CPU --dest-os=$DEST_OS
fi

export LDFLAGS=-lrt
export CPPFLAGS=-D__STDC_FORMAT_MACROS
# make -j${nproc}
make -j12
make install
cp ./out/Release/node ${bindir}
cp ./LICENSE ${prefix}/
install_license ./LICENSE
"""

platforms = [
    x for x in supported_platforms()
    if BinaryBuilder.os(x) ∉ ("freebsd", "windows", "macos") &&
        nbits(x) ≡ 64 && !startswith(arch(x), "powerpc")
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libnode", :libnode),
    ExecutableProduct("node", :node),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

init_block = raw"""
"""

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS, name, version, sources, script,
    platforms, products, dependencies;
    preferred_gcc_version = v"10.2.0",
    init_block = init_block
)

max_stars_repo_name: sunoru/jlnode

max_issues_repo_name: sunoru/jlnode

max_forks_repo_name: sunoru/jlnode


text:
@testset "test exported executive functions" begin

    @testset "test for run_powermodels_pf: ac" begin
        result = run_powermodels_pf(case_pf_ac)

        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
        @test isapprox(result["objective"], 0.0; atol = 1e0)
        @test result["solve_time"] > 0.0
    end

    @testset "test for run_powermodels_pf: ac_native" begin
        result = run_powermodels_pf(case_pf_ac_native)
        @test result["termination_status"]
        @test string(result["optimizer"]) == "NLsolve"
        @test isapprox(result["objective"], 0.0; atol = 1e0)
        @test result["solve_time"] > 0.0
    end

    @testset "test for run_powermodels_opf: ac" begin
        result = run_powermodels_opf(case_opf_ac)

        @test isa(result, Dict{String,Any})
        @test string(result["termination_status"]) == "LOCALLY_SOLVED"

        @test isapprox(result["objective"], 8.0298; atol = 0.1)
        @test result["solve_time"] > 0.0
    end

    @testset "test for run_powermodels_opf: cl" begin
        result = run_powermodels_opf(case_opf_cl)

        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
        @test isapprox(result["objective"], 17015.5; atol = 1e0)
        @test result["solve_time"] > 0.0
    end

    @testset "test for powermodels_tnep: ac" begin
        result = run_powermodels_tnep(case_tnep_ac)

        @test string(result["termination_status"]) == "LOCALLY_SOLVED"
        @test string(result["dual_status"]) == "FEASIBLE_POINT"
        @test string(result["primal_status"]) == "FEASIBLE_POINT"

        new_branch = result["solution"]["ne_branch"]
        for idx in keys(new_branch)
            @test isapprox(new_branch[idx]["built"], 0.0, atol=1e-6, rtol=1e-6) ||
                  isapprox(new_branch[idx]["built"], 1.0, atol=1e-6, rtol=1e-6)
        end

        @test isapprox(new_branch["1"]["pt"], -14.2539, atol = 1e-1, rtol = 1e-1)
        @test isapprox(new_branch["2"]["pf"], 17.7384, atol = 1e-1, rtol = 1e-1)

        @test isapprox(result["objective_lb"], 60.0; atol = 1e-1)
        @test isapprox(result["objective"], 60.0; atol = 1e-1)

        @test result["solve_time"] > 0.0
    end

    @testset "test for powermodels_ots: dc" begin
            result = run_powermodels_ots(case_ots_dc)

            @test string(result["termination_status"]) == "LOCALLY_SOLVED"
            @test string(result["dual_status"]) == "FEASIBLE_POINT"
            @test string(result["primal_status"]) == "FEASIBLE_POINT"

            branch = result["solution"]["branch"]
            for idx in keys(branch)
                @test isapprox(branch[idx]["br_status"], 0.0, atol=1e-6, rtol=1e-6) ||
                      isapprox(branch[idx]["br_status"], 1.0, atol=1e-6, rtol=1e-6)
            end
            @test isapprox(result["objective_lb"], 14810.0; atol = 1e0)
            @test isapprox(result["objective"], 14810.0; atol = 1e0)

            @test result["solve_time"] > 0.0
    end

end


text:
"""
    centroid(simplex::AbstractArray{Float64, 2}) where {T<:Number} -> Array{Float64, 2}

Computes the centroid of a simplex given by `(dim+1)`-by-`dim` array, where
each row is a vertex. Returns the centroid as a vertex (a `1`-by-`dim`
two-dimensional array).

"""
function centroid(simplex::AbstractArray{T, 2}) where {T<:Number}
    n = size(simplex, 1) # Dimension of the space the simplex lives in
    transpose(ones(n, 1)/n) * simplex
end

"""
    radius(simplex::Array{T, 2}, centroid::Array{T, 2}) where {T<:Number}

Compute radius of a simplex (`(dim+1)`-by-`dim` sized `Array{Float64, 2}`
given its centroid (`1`-by-`dim` sized `Array{Float64, 2}`).
"""
function radius(simplex::AbstractArray{T, 2}, centroid::AbstractArray{T, 2}) where {T<:Number}

    # Express vertices with respect to origin
    dim = size(simplex, 2)
    dists_to_centroid = broadcast(-, centroid, simplex) # subtract centroid from all vertices
    maximum(sqrt.(sum(dists_to_centroid.^2, 2)))
end

"""
    radius(s::AbstractArray{T, 2}) where {T<:Number} -> Float64

Compute radius of a simplex `s`, represented by a `Array{Float64, 2}` of size
`(dim+1)`-by-`dim`.
"""
function radius(simplex::AbstractArray{T, 2}) where {T<:Number}

    # Express vertices with respect to origin
    dim = size(simplex, 2)
    dists_to_centroid = broadcast(-, centroid(simplex), simplex) # subtract centroid from all vertices
    maximum(sqrt.(sum(dists_to_centroid .^ 2, 2)))
end

"""
    orientation(simplex::AbstractArray{T, 2}) where {T<:Number} -> Float64

Compute orientation of a `simplex`, represented by a `Array{Float64, 2}` of size `(dim+1)`-by-`dim`.
"""
function orientation(simplex::AbstractArray{T, 2}) where {T<:Number}
    dim = size(simplex, 2)
    hcat(ones(dim + 1, 1), simplex) |> det
end

"""
    volume(simplex::AbstractArray{T, 2}) where {T<:Number} -> Float64

Compute the volume of a `simplex`, represented by a `Array{Float64, 2}` of size `(dim+1)`-by-`dim`.
"""
function volume(simplex::AbstractArray{T, 2}) where {T<:Number}
    orientation(simplex) |> abs
end

"""
    childsimplex(parentsimplex::AbstractArray{T, 2}) where {T<:Number} -> Array{Float64, 2}

Generates a random simplex which is entirely contained within `parentsimplex`,
which is a (dim+1)-by-dim array.
"""
function childsimplex(parentsimplex::AbstractArray{T, 2}) where {T}
    # Convex expansion coefficients of the random simplex
    dim = size(parentsimplex, 2)
    rs = rand(dim + 1, dim + 1)
    normalised_colsums = 1 ./ sum(rs, dims = 2)
    (normalised_colsums .* rs) * parentsimplex
end

"""
    issingular(simplex::AbstractArray{T, 2}) where {T<:Number}

Determines if a simplex is singular by checking if any of its vertices are
identical.
"""
function issingular(simplex::AbstractArray{T, 2}) where {T<:Number}
    size(unique(simplex, 1), 1) != size(simplex, 1)
end

"""
    insidepoints(npts::Int, parentsimplex::AbstractArray{T, 2}) where {T<:Number}

Generates `npts` points that located inside `parentsimplex`.
"""
function insidepoints(npts::Int, parentsimplex::AbstractArray{T, 2}) where {T<:Number}
    dim = size(parentsimplex, 2)
    # Random linear combination coefficients
    R = rand(Uniform(), npts, dim + 1)

    # Normalise the coefficients so that they sum to one. We can then create the new point
    # as a convex linear combination of the vertices of the parent simplex.
    normalised_coeffs = (1 ./ sum(R, dims = 2)) .* R
    normalised_coeffs * parentsimplex
end

"""
    outsidepoint(parentsimplex::AbstractArray{T, 2}) where {T<:Number}

Generate a single point that is guaranteed to lie outside `parentsimplex`.
"""
function outsidepoint(parentsimplex::AbstractArray{T, 2}) where {T<:Number}
    dim = size(parentsimplex, 2)
    # Random linear combination coefficients
    R = rand(1, dim + 1)

    # Normalise the coefficients so that they sum to one. We can then create the new point
    # as a convex linear combination of the vertices of the parent simplex.
    normalised_coeffs = (1 ./ sum(R, dims = 2)) .* R
    normalised_coeffs[1] += (1 - sum(normalised_coeffs[2:dim+1])) + rand()

    normalised_coeffs * parentsimplex
end

"""
    outsidepoints(npts::Int, parentsimplex::AbstractArray{T, 2}) where T <: Number

Generates `npts` points that located outside `parentsimplex`.
"""
function outsidepoints(npts::Int, parentsimplex::AbstractArray{T, 2}) where {T<:Number}
    vcat([outsidepoint(parentsimplex) for i in 1:npts]...)
end

"""
    nontrivially_intersecting_simplices(dim::Int) -> Array{Float64, 2}

Genereate a set of non-trivially intersecting `dim`-dimensional simplices
(i.e. they don't intersect along boundaries or vertices only).
"""
function nontrivially_intersecting_simplices(dim::Int)
    rs = rand(dim + 1, dim) # random simplex with dim+1 vertices
    n_inside = rand(1:dim, 1)[1]

    verts_inside = insidepoints(n_inside, rs)

    verts_outside = vcat([outsidepoint(rs) for i in 1:(dim+1-n_inside)]...)

    return rs, vcat(verts_inside, verts_outside)
end

"""
    simplices_sharing_vertices(dim::Int) -> Array{Float64, 2}

Genereate a set of non-trivially intersecting `dim`-dimensional simplices
(i.e. they don't intersect along boundaries or vertices only).
"""
function simplices_sharing_vertices(dim::Int)
    rs = rand(dim + 1, dim) # random simplex with dim+1 vertices

    # Determine how many vertices should be shared and how many are remaining
    n_shared_verts = rand(1:dim, 1)[1]
    n_remaining_verts = dim + 1 - n_shared_verts

    # Pick the shared vertices and generate remaining vertices
    new_simplex = zeros(Float64, 0, dim)
    shared_verts = rs[rand(1:n_shared_verts, n_shared_verts), :]

    if ndims(shared_verts) == 1
        new_simplex = vcat(new_simplex, transpose(shared_verts))
    else
        new_simplex = vcat(new_simplex, shared_verts)
    end
    for i = 1:n_remaining_verts
        if rand(Bool) == true
            new_simplex  = vcat(new_simplex, insidepoints(1, rs))
        else
            new_simplex  = vcat(new_simplex, outsidepoints(1, rs))
        end
    end

    return rs, new_simplex
end

"""
    intersecting_simplices(;dim::Int = 3, intersection_type = "nontrivial")

Generate a pair of `dim` dimensional intersecting simplices. Each resulting
simplex is a `(dim+1)`-by-`dim` two-dimensional array. The `intersection_type`
argument can be either `"nontrivial"` or `"sharingvertices"`.
"""
function intersecting_simplices(;dim::Int = 3, intersection_type = "nontrivial")
    if intersection_type == "nontrivial"
        return nontrivially_intersecting_simplices(dim)
    elseif intersection_type == "sharingvertices"
        return simplices_sharing_vertices(dim)
    end
end


text:
const AV = AbstractVector
const AM = AbstractMatrix
const AA = AbstractArray
const IN = Integer
const RN = Real
const FN = AbstractFloat


text:
using MLDatasets: MNIST
using Flux.Data: DataLoader
using Flux: onehotbatch
using Parameters: @with_kw

@with_kw mutable struct Args
    minibath_size::Int = 128  # Size of mini-batch
end

function get_data(args)

    # Load the MNIST train and test data from MLDatasets
    train_x, train_y = MNIST.traindata(Float32)
    test_x, test_y = MNIST.testdata(Float32)

    # Reshape data to 28x28x1 multi-dimensional array
    train_x = reshape(train_x, 28, 28, 1, :)
    test_x = reshape(test_x, 28, 28, 1, :)

    # Labels must be encoded as a vector with the same dimension
    # as the number of categories (unique handwritten digits) in the data set
    train_y, test_y = onehotbatch(train_y, 0:9), onehotbatch(test_y, 0:9)

    # Now, we load the train and test images and their labels onto a DataLoader object
    data_loader_train = DataLoader(train_x, train_y, batchsize=args.minibath_size, shuffle=true)
    data_loader_test = DataLoader(train_x, train_y, batchsize=args.minibath_size, shuffle=true)

    return data_loader_train, data_loader_test

end

function train(; kws...)
    args = Args(; kws...)

    @info("Loading data...")
    data_loader_train, data_loader_test = get_data(args)

    # Iterating over train data
    for (x, y) in data_loader_train
        @assert size(x) == (28, 28, 1, 128) || size(x) == (28, 28, 1, 96)
        @assert size(y) == (10, 128) || size(y) == (10, 96)
     end

end

cd(@__DIR__)
train()

max_stars_repo_name: shawn-kern/model-zoo

max_issues_repo_name: shawn-kern/model-zoo

max_forks_repo_name: shawn-kern/model-zoo


text:
"""
    Topology2D(mesh::AbstractMesh2D, elemorder=Mesh.elements(mesh))

This is a generic non-distributed topology for 2D meshes. `elemorder` is a
vector or other linear ordering of the `Mesh.elements(mesh)`.
"""
struct Topology2D{M <: Meshes.AbstractMesh{2}, EO, OI, BF} <: AbstractTopology
    "mesh on which the topology is constructed"
    mesh::M
    "`elemorder[e]` should give the `e`th element "
    elemorder::EO
    "the inverse of `elemorder`: `orderindex[elemorder[e]] == e`"
    orderindex::OI
    "a vector of all unique internal faces, (e, face, o, oface, reversed)"
    internal_faces::Vector{Tuple{Int, Int, Int, Int, Bool}}
    "the collection of all (e,vert) pairs"
    vertices::Vector{Tuple{Int, Int}}
    "the index in `vertices` of the first `(e,vert)` pair for all unique vertices"
    vertex_offset::Vector{Int}
    "a NamedTuple of vectors of `(e,face)` "
    boundaries::BF
end

function Base.show(io::IO, topology::Topology2D)
    print(io, "Topology2D on ", topology.mesh)
end

function Topology2D(
    mesh::Meshes.AbstractMesh{2},
    elemorder = Meshes.elements(mesh),
    orderindex = Meshes.linearindices(elemorder),
)
    internal_faces = Tuple{Int, Int, Int, Int, Bool}[]
    vertices = sizehint!(Tuple{Int, Int}[], length(elemorder) * 4)
    vertex_offset = sizehint!(Int[1], length(elemorder))
    boundaries = NamedTuple(
        boundary_name => Tuple{Int, Int}[] for
        boundary_name in Meshes.boundary_names(mesh)
    )
    for (e, elem) in enumerate(elemorder)
        for face in 1:4
            if Meshes.is_boundary_face(mesh, elem, face)
                boundary_name = Meshes.boundary_face_name(mesh, elem, face)
                push!(boundaries[boundary_name], (e, face))
            else
                opelem, opface, reversed =
                    Meshes.opposing_face(mesh, elem, face)
                o = orderindex[opelem]
                if (o, opface) < (e, face)
                    # new internal face
                    push!(internal_faces, (e, face, o, opface, reversed))
                end
            end
        end
        for vert in 1:4
            if !any(Meshes.SharedVertices(mesh, elem, vert)) do (velem, vvert)
                (orderindex[velem], vvert) < (e, vert)
            end
                # new vertex
                for (velem, vvert) in Meshes.SharedVertices(mesh, elem, vert)
                    o = orderindex[velem]
                    push!(vertices, (o, vvert))
                end
                push!(vertex_offset, length(vertices) + 1)
            end
        end
    end
    @assert length(vertices) == 4 * length(elemorder) == vertex_offset[end] - 1
    return Topology2D(
        mesh,
        elemorder,
        orderindex,
        internal_faces,
        vertices,
        vertex_offset,
        boundaries,
    )
end

domain(topology::Topology2D) = domain(topology.mesh)
nlocalelems(topology::Topology2D) = length(topology.elemorder)
vertex_coordinates(topology::Topology2D, e::Int) =
    ntuple(4) do vert
        Meshes.coordinates(topology.mesh, topology.elemorder[e], vert)
    end

function opposing_face(topology::Topology2D, e::Int, face::Int)
    mesh = topology.mesh
    elem = topology.elemorder[e]
    if Meshes.is_boundary_face(mesh, elem, face)
        boundary_name = Meshes.boundary_face_name(mesh, elem, face)
        b = findfirst(==(boundary_name), keys(topology.boundaries))
        return (0, b, false)
    end
    opelem, opface, reversed = Meshes.opposing_face(mesh, elem, face)
    return (topology.orderindex[opelem], opface, reversed)
end
interior_faces(topology::Topology2D) = topology.internal_faces

Base.length(vertiter::VertexIterator{<:Topology2D}) =
    length(vertiter.topology.vertex_offset) - 1
Base.eltype(::VertexIterator{T}) where {T <: Topology2D} = Vertex{T, Int}
function Base.iterate(vertiter::VertexIterator{<:Topology2D}, uvert = 1)
    topology = vertiter.topology
    if uvert >= length(topology.vertex_offset)
        return nothing
    end
    return Vertex(topology, uvert), uvert + 1
end

Base.length(vertex::Vertex{<:Topology2D}) =
    vertex.topology.vertex_offset[vertex.num + 1] -
    vertex.topology.vertex_offset[vertex.num]
Base.eltype(vertex::Vertex{<:Topology2D}) = eltype(vertex.topology.vertices)
function Base.iterate(
    vertex::Vertex{<:Topology2D},
    idx = vertex.topology.vertex_offset[vertex.num],
)
    if idx >= vertex.topology.vertex_offset[vertex.num + 1]
        return nothing
    end
    return vertex.topology.vertices[idx], idx + 1
end

boundary_names(topology::Topology2D) = keys(topology.boundaries)
boundary_tags(topology::Topology2D) = NamedTuple{boundary_names(topology)}(
    ntuple(i -> i, length(topology.boundaries)),
)
boundary_tag(topology::Topology2D, boundary_name::Symbol) =
    findfirst(==(boundary_name), boundary_names(topology))

boundary_faces(topology::Topology2D, boundary) = topology.boundaries[boundary]


text:
# this allows writing `forever() do ... end`
forever(f) = (f() for _ in Iterators.repeated(nothing))

function last(itr)
    res = iterate(itr)
    if res === nothing
        throw(ArgumentError("iterator must be non-empty"))
    end
    output, state = res
    for elem in Iterators.rest(itr, state)
        output = elem
    end
    return output
end

mutable struct Settable{T}
    v::T
end

set!(b::Settable{T}, v::T) where T = (b.v = v)

Base.IteratorSize(::Type{<:Settable}) = Base.IsInfinite()

Base.iterate(b::Settable, state=nothing) = b.v, nothing


text:
"""
    Julia Position Based Fluid

          (C) Hipgineer
             Jaybin
              2021

    dep : LightXML, Printf

"""
module  JuliaPBF
const   JULIAPBF_VERSION = v"0.0.1"
export  StructPBF,
        IO_dev,
        Solver

        ##############
        # DATA TYPES #
        ##############
        include("Struct/StructPBF.jl")  # Types are defined in the StructPBF Module
        using .StructPBF                # Load StructPBF Module

        ##############
        # In&Outputs #
        ##############
        include("IO_dev/IO.jl")
        using .IO_dev

        ##############
        # MainSolver #
        ##############
        include("Solver/Solver.jl")
        using .Solver
end

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
using OPFLearn
using Documenter

DocMeta.setdocmeta!(OPFLearn, :DocTestSetup, :(using OPFLearn); recursive=true)

makedocs(;
    modules=[OPFLearn],
    authors="Trager Joswig-Jones",
    sitename="OPFLearn",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://NREL.github.io/OPFLearn.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
		"Manual" => [
		"Getting Started" => "quickstartguide.md",
		"OPFLearn Framework" => "framework.md",
		"Result Data" => "results.md",
		"Dataset Format" => "dataset.md",
		"Distributed Processing" => "distributed.md",
		"Python Interface" => "python.md",
		],
		"Library" => [
		"Dataset Creation" => "datasetcreation.md",
		"Exporting Datasets" => "saving.md",
		"Sampling" => "sampling.md"
		]
    ],
)

deploydocs(;
    repo="github.com/NREL/OPFLearn.jl.git",
    devbranch = "main"
)


text:
export Mat
export Vec
export Point
export unit

const Mat = SMatrix

function unit(::Type{T}, i::Integer) where T <: StaticVector
    T(ntuple(Val(length(T))) do j
        ifelse(i == j, 1, 0)
    end)
end

macro fixed_vector(name, parent)
    esc(quote
        struct $(name){S, T} <: $(parent){S, T}
            data::NTuple{S, T}

            function $(name){S, T}(x::NTuple{S,T}) where {S, T}
                new{S, T}(x)
            end

            function $(name){S, T}(x::NTuple{S,Any}) where {S, T}
                new{S, T}(StaticArrays.convert_ntuple(T, x))
            end
        end
        size_or(::Type{$(name)}, or) = or
        eltype_or(::Type{$(name)}, or) = or
        eltype_or(::Type{$(name){S, T} where S}, or) where {T} = T
        eltype_or(::Type{$(name){S, T} where T}, or) where {S} = or
        eltype_or(::Type{$(name){S, T}}, or) where {S, T} = T

        size_or(::Type{$(name){S, T} where S}, or) where {T} = or
        size_or(::Type{$(name){S, T} where T}, or) where {S} = Size{(S,)}()
        size_or(::Type{$(name){S, T}}, or) where {S, T} = (S,)
        # Array constructor
        @inline function $(name){S}(x::AbstractVector{T}) where {S, T}
            @assert S <= length(x)
            $(name){S, T}(ntuple(i-> x[i], Val(S)))
        end
        @inline function $(name){S, T1}(x::AbstractVector{T2}) where {S, T1, T2}
            @assert S <= length(x)
            $(name){S, T1}(ntuple(i-> T1(x[i]), Val(S)))
        end

        @inline function $(name){S, T}(x) where {S, T}
            $(name){S, T}(ntuple(i-> T(x), Val(S)))
        end

        @inline function $(name){S}(x::T) where {S, T}
            $(name){S, T}(ntuple(i-> x, Val(S)))
        end
        @inline function $(name){1, T}(x::T) where T
            $(name){1, T}((x,))
        end
        @inline $(name)(x::NTuple{S}) where {S} = $(name){S}(x)
        @inline $(name)(x::T) where {S, T <: Tuple{Vararg{Any, S}}} = $(name){S, StaticArrays.promote_tuple_eltype(T)}(x)
        @inline function $(name){S}(x::T) where {S, T <: Tuple}
            $(name){S, StaticArrays.promote_tuple_eltype(T)}(x)
        end
        $(name){S, T}(x::StaticVector) where {S, T} = $(name){S, T}(Tuple(x))
        @generated function (::Type{$(name){S, T}})(x::$(name)) where {S, T}
            idx = [:(x[$i]) for i = 1:S]
            quote
                $($(name)){S, T}($(idx...))
            end
        end
        @generated function convert(::Type{$(name){S, T}}, x::$(name)) where {S, T}
            idx = [:(x[$i]) for i = 1:S]
            quote
                $($(name)){S, T}($(idx...))
            end
        end
        @generated function (::Type{SV})(x::StaticVector) where SV <: $(name)
            len = size_or(SV, size(x))[1]
            if length(x) == len
                :(SV(Tuple(x)))
            elseif length(x) > len
                elems = [:(x[$i]) for i = 1:len]
                :(SV($(Expr(:tuple, elems...))))
            else
                error("Static Vector too short: $x, target type: $SV")
            end
        end

        Base.@pure StaticArrays.Size(::Type{$(name){S, Any}}) where {S} = Size(S)
        Base.@pure StaticArrays.Size(::Type{$(name){S, T}}) where {S,T} = Size(S)

        Base.@propagate_inbounds function Base.getindex(v::$(name){S, T}, i::Int) where {S, T}
            v.data[i]
        end
        @inline Base.Tuple(v::$(name)) = v.data
        @inline Base.convert(::Type{$(name){S, T}}, x::NTuple{S, T}) where {S, T} = $(name){S, T}(x)
        @inline function Base.convert(::Type{$(name){S, T}}, x::Tuple) where {S, T}
            $(name){S, T}(convert(NTuple{S, T}, x))
        end

        @generated function StaticArrays.similar_type(::Type{SV}, ::Type{T}, s::Size{S}) where {SV <: $(name), T, S}
            if length(S) === 1
                $(name){S[1], T}
            else
                StaticArrays.default_similar_type(T,s(),Val{length(S)})
            end
        end

    end)
end


text:
module TestBenchSmoke

using Test
using RecordArraysBenchmarks: clear, setup_smoke

@testset "smoke test benchmarks" begin
    try
        local suite
        @test (suite = setup_smoke()) isa Any
        @test run(suite) isa Any
    finally
        clear()
    end
end

end  # module


text:
using CardiovascularModels
using Test

using OrdinaryDiffEq, ModelingToolkit

@variables t
D = Differential(t)

@testset "3-compartment model" begin
    include("test-three-compartment.jl")
end

@testset "Passive connection of two vessels" begin
    include("passive-connection.jl")
end


text:
# This file is a part of GERDAMetadata.jl, licensed under the MIT License (MIT).

__precompile__(true)

module GERDAMetadata

using Dates
using Glob
using JSON
using PropDicts

include("filekey.jl")
include("dataset.jl")
include("dataconfig.jl")
include("calibcatalog.jl")
include("gerdadata.jl")
include("calfunc.jl")

end # module


text:
# ChrBase Chr type and core functions
# Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
# Licensed under MIT License, see LICENSE.md

struct Chr{CS<:CharSet,T<:CodeUnitTypes} <: AbstractChar
    v::T
    (::Type{Chr})(::Type{CS}, v::T) where {CS<:CharSet,T<:CodeUnitTypes} = new{CS,T}(v)
end

(::Type{<:Chr{CS,T}})(v::Number) where {CS<:CharSet,T<:CodeUnitTypes} = Chr(CS, T(v))

for lst in cse_info
    nam, typ = lst
    length(lst) < 3 || continue
    chrnam = symstr(nam, "Chr")
    if String(nam)[1] != '_'
        @eval @api public $chrnam
    elseif nam == :_Latin
        @api develop _LatinChr
    else
        continue
    end
    cs = symstr(nam, "CharSet")
    @eval const $chrnam = Chr{$cs, $typ}
    @eval show(io::IO, ::Type{$chrnam}) = print(io, $(quotesym(chrnam)))
    @eval codepoint_cse(::Type{$chrnam}) = $(symstr(nam,"CSE"))
    @eval (::Type{$chrnam})(v::Number) = Chr($cs, $typ(v))
end

codepoint(ch::Chr) = ch.v
basetype(::Type{<:Chr{CS,B}}) where {CS,B} = B
charset(::Type{<:Chr{CS,B}}) where {CS,B} = CS
typemin(::Type{T}) where {CS,B,T<:Chr{CS,B}} = Chr(CS, typemin(B))
typemax(::Type{T}) where {CS,B,T<:Chr{CS,B}} = Chr(CS, typemax(B))

const LatinChars   = Union{LatinChr, _LatinChr}
const ByteChars    = Union{ASCIIChr, LatinChr, _LatinChr, Text1Chr}
const WideChars    = Union{UCS2Chr, UTF32Chr}

codepoint_cse(::Type{Char}) = RawUTF8CSE

const AbsChar = @static isdefined(Base, :AbstractChar) ? AbstractChar : Union{Char, AbstractChar}

# Promotion rules for characters

promote_rule(::Type{T}, ::Type{T}) where {T<:Chr} = T

promote_rule(::Type{Char}, ::Type{<:Chr})   = Char

Base.need_full_hex(c::Chr) = is_hex_digit(c)
Base.escape_nul(c::Chr) = ('0' <= c <= '7') ? "\\x00" : "\\0"

bytoff(::Type{UInt8},  off) = off
bytoff(::Type{UInt16}, off) = off << 1
bytoff(::Type{UInt32}, off) = off << 2
chroff(::Type{UInt8},  off) = off
chroff(::Type{UInt16}, off) = off >>> 1
chroff(::Type{UInt32}, off) = off >>> 2

chrdiff(pnt::Ptr{T}, beg::Ptr{T}) where {T<:CodeUnitTypes} = Int(chroff(T, pnt - beg))

bytoff(pnt::Ptr{T}, off) where {T<:CodeUnitTypes} = pnt + bytoff(T, off)

typemax(::Type{ASCIIChr}) = ASCIIChr(0x7f)
typemax(::Type{UTF32Chr}) = UTF32Chr(0x10ffff)

basetype(::Type{Char})        = UInt32
basetype(::Type{T}) where {T<:CodeUnitTypes} = T

convert(::Type{T}, v::S) where {T<:Integer, S<:Chr} = convert(T, codepoint(v))::T
convert(::Type{T}, v::Signed) where {T<:Chr} =
    (v >= 0 && is_valid(T, v%Unsigned)) ? convert(T, v%Unsigned) : codepoint_error(T, v)
convert(::Type{T}, v::Unsigned) where {CS,B,T<:Chr{CS,B}} =
    is_valid(T, v) ? Chr(CS, v%B) : codepoint_error(T, v)
convert(::Type{Char}, v::T) where {T<:Chr} = convert(Char, codepoint(v))
convert(::Type{T}, v::Char) where {T<:Chr} = convert(T, codepoint(v))::T

rem(x::Number, ::Type{<:Chr{CS,B}}) where {CS,B} = Chr(CS, x%B)
rem(x::Char, ::Type{T}) where {T<:Chr}   = x%UInt32%T
rem(x::Chr, ::Type{T}) where {T<:Chr}    = (x.v)%T
rem(x::Chr, ::Type{T}) where {T<:Char}   = (x.v)%T
rem(x::Chr, ::Type{T}) where {T<:Number} = (x.v)%T

(::Type{S})(v::T) where {S<:Union{UInt32, Int, UInt}, T<:Chr} = codepoint(v)%S
(::Type{Char})(v::Chr) = Char(codepoint(v))
(::Type{T})(v::Char) where {T<:Chr} = T(codepoint(v))

eltype(::Type{T}) where {T<:Chr} = T
size(cp::Chr, dim) = convert(Int, dim) < 1 ? boundserr(cp, dim) : 1
getindex(cp::Chr, i::Integer) = i == 1 ? cp : boundserr(cp, i)
getindex(cp::Chr, I::Integer...) = all(x -> x == 1, I) ? cp : boundserr(cp, I)

_uni_rng(m) = 0x00000:ifelse(m < 0xd800, m, m-0x800)
codepoint_rng(::Type{T}) where {T<:Chr} = _uni_rng(typemax(T)%UInt32)
codepoint_rng(::Type{Char}) = _uni_rng(0x10ffff)
codepoint_rng(::Type{Text2Chr}) = 0%UInt16:typemax(UInt16)
codepoint_rng(::Type{Text4Chr}) = 0%UInt32:typemax(UInt32)

codepoint_adj(::Type{T}, ch) where {T} = ifelse(ch < 0xd800, ch, ch+0x800)%T
codepoint_adj(::Type{T}, ch) where {T<:Union{Text2Chr,Text4Chr}} = ch%T

# returns a random valid Unicode scalar value in the correct range for the type of character
@static if V6_COMPAT
    import Base.Random: rand!, rand, AbstractRNG
    rand(r::AbstractRNG, ::Type{T}) where {T<:Chr} =
        codepoint_adj(T, rand(r, codepoint_rng(T)))
    rand!(rng::AbstractRNG, A::AbstractArray, r::UnitRange{<:Chr}) =
        rand!(rng, A, Base.Random.RangeGenerator(r))
else
    import Random: rand!, rand, AbstractRNG, SamplerType

    rand(r::AbstractRNG, ::SamplerType{T}) where {T<:Chr} =
        codepoint_adj(T, rand(r, codepoint_rng(T)))
end

==(x::Chr, y::AbsChar) = codepoint(x) == codepoint(y)
==(x::AbsChar, y::Chr) = codepoint(x) == codepoint(y)
==(x::Chr, y::Chr)     = codepoint(x) == codepoint(y)

isless(x::Chr, y::AbsChar) = codepoint(x) < codepoint(y)
isless(x::AbsChar, y::Chr) = codepoint(x) < codepoint(y)
isless(x::Chr, y::Chr)     = codepoint(x) < codepoint(y)

# This is so that the hash is compatible with isless, but it's very inefficient
Base.hash(x::Chr, h::UInt) = hash(Char(x), h)

# Support functions for UTF-8 handling
@inline get_utf8_2(ch) =
    (0xc0 | ((ch >>> 6)%UInt8), 0x80 | ((ch & 0x3f)%UInt8))
@inline get_utf8_3(ch) =
    (0xe0 | ((ch >>> 12)%UInt8), 0x80 | ((ch >>> 6) & 0x3f)%UInt8, 0x80 | ((ch & 0x3f)%UInt8))
@inline get_utf8_4(ch) =
    (0xf0 | ((ch >>>  18)%UInt8), 0x80 | (((ch >>> 12) & 0x3f)%UInt8),
     0x80 | ((ch >>>  6) & 0x3f)%UInt8, 0x80 | ((ch & 0x3f)%UInt8))

# Little-endian output here
@inline get_utf8_16(ch) =
    (ch >>> 6) | (((ch & 0x3f)%UInt16)<<8) | 0x80c0
@inline get_utf8_32(ch) =
    ((ch & 0xc0000) >>> 18) | ((ch & 0x3f000) >>> 4) |
    ((ch & 0xfc0) << 10) | (ch & 0x3f)<<24 | 0x808080f0

utf_trail(c::UInt8) = (0xe5000000 >>> ((c & 0xf0) >> 3)) & 0x3

is_valid_continuation(c) = ((c & 0xc0) == 0x80)

# Support functions for UTF-16 handling
@inline get_utf16(ch) = (0xd7c0 + (ch >> 10))%UInt16, (0xdc00 + (ch & 0x3ff))%UInt16

@inline get_utf16_32(ch) =
    ((0xd7c0 + (ch >>> 10))%UInt16) << 6 | (0xdc00 + (ch & 0x3ff))%UInt32

is_surrogate_lead(c::Unsigned) = ((c & ~0x003ff) == 0xd800)
is_surrogate_trail(c::Unsigned) = ((c & ~0x003ff) == 0xdc00)
is_surrogate_codeunit(c::Unsigned) = ((c & ~0x007ff) == 0xd800)

function _isvalid end
function _isvalid_chr end