text:
"""
This module provides the [`getconnection`](@ref) function with support for:
- Opening TCP and SSL connections.
- Reusing connections for multiple Request/Response Messages,
- Pipelining Request/Response Messages. i.e. allowing a new Request to be
  sent before previous Responses have been read.

This module defines a [`Connection`](@ref)
struct to manage pipelining and connection reuse and a
[`Transaction`](@ref)`<: IO` struct to manage a single
pipelined request. Methods are provided for `eof`, `readavailable`,
`unsafe_write` and `close`.
This allows the `Transaction` object to act as a proxy for the
`TCPSocket` or `SSLContext` that it wraps.

The [`pool`](@ref) is a collection of open
`Connection`s.  The `request` function calls `getconnection` to
retrieve a connection from the `pool`.  When the `request` function
has written a Request Message it calls `closewrite` to signal that
the `Connection` can be reused for writing (to send the next Request).
When the `request` function has read the Response Message it calls
`closeread` to signal that the `Connection` can be reused for
reading.
"""

module ConnectionPool

export Connection, Transaction,
       getconnection, getparser, getrawstream, inactiveseconds

using ..IOExtras

import ..@debug, ..@debugshow, ..DEBUG_LEVEL, ..taskid
import ..@require, ..precondition_error, ..@ensure, ..postcondition_error
using MbedTLS: SSLConfig, SSLContext, setup!, associate!, hostname!, handshake!

import ..Parser

const default_connection_limit = 8
const default_pipeline_limit = 16
const nolimit = typemax(Int)

const nobytes = view(UInt8[], 1:0)
const ByteView = typeof(nobytes)
byteview(bytes::ByteView) = bytes
byteview(bytes)::ByteView = view(bytes, 1:length(bytes))

"""
    Connection{T <: IO}

A `TCPSocket` or `SSLContext` connection to a HTTP `host` and `port`.

Fields:
- `host::String`
- `port::String`, exactly as specified in the URI (i.e. may be empty).
- `pipeline_limit`, number of requests to send before waiting for responses.
- `peerport`, remote TCP port number (used for debug messages).
- `localport`, local TCP port number (used for debug messages).
- `io::T`, the `TCPSocket` or `SSLContext.
- `excess::ByteView`, left over bytes read from the connection after
   the end of a response message. These bytes are probably the start of the
   next response message.
- `sequence`, number of most recent `Transaction`.
- `writecount`, number of Messages that have been written.
- `writedone`, signal that `writecount` was incremented.
- `readcount`, number of Messages that have been read.
- `readdone`, signal that `readcount` was incremented.
- `timestamp`, time data was last recieved.
- `parser::Parser`, reuse a `Parser` when this `Connection` is reused.
"""

mutable struct Connection{T <: IO}
    host::String
    port::String
    pipeline_limit::Int
    peerport::UInt16
    localport::UInt16
    io::T
    excess::ByteView
    sequence::Int
    writecount::Int
    writebusy::Bool
    writedone::Condition
    readcount::Int
    readbusy::Bool
    readdone::Condition
    timestamp::Float64
    parser::Parser
end

"""
A single pipelined HTTP Request/Response transaction`.

Fields:
 - `c`, the shared [`Connection`](@ref) used for this `Transaction`.
 - `sequence::Int`, identifies this `Transaction` among the others that share `c`.
"""

struct Transaction{T <: IO} <: IO
    c::Connection{T}
    sequence::Int
end

Connection(host::AbstractString, port::AbstractString,
           pipeline_limit::Int, io::T) where T <: IO =
    Connection{T}(host, port, pipeline_limit,
                  peerport(io), localport(io),
                  io, nobytes,
                  -1,
                  0, false, Condition(),
                  0, false, Condition(),
                  0, Parser())

Transaction(c::Connection{T}) where T <: IO =
    Transaction{T}(c, (c.sequence += 1))

function client_transaction(c)
    t = Transaction(c)
    startwrite(t)
    return t
end

getparser(t::Transaction) = t.c.parser

getrawstream(t::Transaction) = t.c.io

inactiveseconds(t::Transaction) = inactiveseconds(t.c)

function inactiveseconds(c::Connection)::Float64
    if !c.readbusy && !c.writebusy
        return Float64(0)
    end
    return time() - c.timestamp
end

Base.unsafe_write(t::Transaction, p::Ptr{UInt8}, n::UInt) =
    unsafe_write(t.c.io, p, n)

Base.isopen(c::Connection) = isopen(c.io)

Base.isopen(t::Transaction) = isopen(t.c) &&
                              t.c.readcount <= t.sequence &&
                              t.c.writecount <= t.sequence

function Base.eof(t::Transaction)
    @require isreadable(t) || !isopen(t)
    if nb_available(t) > 0
        return false
    end                 ;@debug 4 "eof(::Transaction) -> eof($typeof(c.io)): $t"
    return eof(t.c.io)
end

Base.nb_available(t::Transaction) = nb_available(t.c)
Base.nb_available(c::Connection) =
    !isempty(c.excess) ? length(c.excess) : nb_available(c.io)

Base.isreadable(t::Transaction) = t.c.readbusy && t.c.readcount == t.sequence

Base.iswritable(t::Transaction) = t.c.writebusy && t.c.writecount == t.sequence

function Base.readavailable(t::Transaction)::ByteView
    @require isreadable(t)
    if !isempty(t.c.excess)
        bytes = t.c.excess
        @debug 4 "↩  read $(length(bytes))-bytes from excess buffer."
        t.c.excess = nobytes
    else
        bytes = byteview(readavailable(t.c.io))
        @debug 4 "⬅  read $(length(bytes))-bytes from $(typeof(t.c.io))"
    end
    t.c.timestamp = time()
    return bytes
end

@static if !isdefined(Base, :copyto!)
    const copyto! = copy!
end

function Base.readbytes!(t::Transaction, a::Vector{UInt8}, nb::Int)

    if !isempty(t.c.excess)
        l = length(t.c.excess)
        copyto!(a, 1, t.c.excess, 1, min(l, nb))
        if l > nb
            t.c.excess = view(t.c.excess, nb+1:l)
            return nb
        else
            t.c.excess = nobytes
            return l + readbytes!(t.c.io, view(a, l+1:nb))
        end
    end

    return readbytes!(t.c.io, a, nb)
end

"""
    unread!(::Transaction, bytes)

Push bytes back into a connection's `excess` buffer
(to be returned by the next read).
"""

function IOExtras.unread!(t::Transaction, bytes::ByteView)
    @require isreadable(t)
    t.c.excess = bytes
    return
end

"""
    startwrite(::Transaction)

Wait for prior pending writes to complete.
"""

function IOExtras.startwrite(t::Transaction)
    @require !iswritable(t)                     ;t.c.writecount != t.sequence &&
                                                   @debug 1 "⏳  Wait write: $t"
    while t.c.writecount != t.sequence
        wait(t.c.writedone)
    end                                           ;@debug 2 "👁  Start write:$t"
    t.c.writebusy = true
    @ensure iswritable(t)
    return
end

"""
    closewrite(::Transaction)

Signal that an entire Request Message has been written to the `Transaction`.
"""

function IOExtras.closewrite(t::Transaction)
    @require iswritable(t)

    t.c.writebusy = false
    t.c.writecount += 1                           ;@debug 2 "🗣  Write done: $t"
    notify(t.c.writedone)
    notify(poolcondition)

    @ensure !iswritable(t)
    return
end

"""
    startread(::Transaction)

Wait for prior pending reads to complete.
"""

function IOExtras.startread(t::Transaction)
    @require !isreadable(t)                      ;t.c.readcount != t.sequence &&
                                                   @debug 1 "⏳  Wait read:  $t"
    t.c.timestamp = time()
    while t.c.readcount != t.sequence
        wait(t.c.readdone)
    end                                           ;@debug 2 "👁  Start read: $t"
    t.c.readbusy = true
    @ensure isreadable(t)
    return
end

"""
    closeread(::Transaction)

Signal that an entire Response Message has been read from the `Transaction`.

Increment `readcount` and wake up tasks waiting in `startread`.
"""

function IOExtras.closeread(t::Transaction)
    @require isreadable(t)

    t.c.readbusy = false
    t.c.readcount += 1
    notify(t.c.readdone)                          ;@debug 2 "✉  Read done:  $t"
    notify(poolcondition)

    @ensure !isreadable(t)
    return
end

function Base.close(t::Transaction)
    close(t.c)
    if iswritable(t)
        closewrite(t)
    end
    if isreadable(t)
        closeread(t)
    end
    return
end

function Base.close(c::Connection)
    close(c.io)
    if nb_available(c) > 0
        purge(c)
    end
    notify(poolcondition)
    return
end

"""
    purge(::Connection)

Remove unread data from a `Connection`.
"""

function purge(c::Connection)
    @require !isopen(c.io)
    while !eof(c.io)
        readavailable(c.io)
    end
    c.excess = nobytes
    @ensure nb_available(c) == 0
end

"""
The `pool` is a collection of open `Connection`s.  The `request`
function calls `getconnection` to retrieve a connection from the
`pool`.  When the `request` function has written a Request Message
it calls `closewrite` to signal that the `Connection` can be reused
for writing (to send the next Request). When the `request` function
has read the Response Message it calls `closeread` to signal that
the `Connection` can be reused for reading.
"""

const pool = Vector{Connection}()
const poollock = ReentrantLock()
const poolcondition = Condition()

"""
    closeall()

Close all connections in `pool`.
"""

function closeall()

    lock(poollock)
    for c in pool
        close(c)
    end
    empty!(pool)
    unlock(poollock)
    notify(poolcondition)
    return
end

"""
    findwritable(type, host, port) -> Vector{Connection}

Find `Connections` in the `pool` that are ready for writing.
"""

function findwritable(T::Type,
                      host::AbstractString,
                      port::AbstractString,
                      pipeline_limit::Int,
                      reuse_limit::Int)

    filter(c->(!c.writebusy &&
               typeof(c.io) == T &&
               c.host == host &&
               c.port == port &&
               c.pipeline_limit == pipeline_limit &&
               c.writecount < reuse_limit &&
               c.writecount - c.readcount < pipeline_limit + 1 &&
               isopen(c.io)), pool)
end

"""
    findoverused(type, host, port, reuse_limit) -> Vector{Connection}

Find `Connections` in the `pool` that are over the reuse limit
and have no more active readers.
"""

function findoverused(T::Type,
                      host::AbstractString,
                      port::AbstractString,
                      reuse_limit::Int)

    filter(c->(typeof(c.io) == T &&
               c.host == host &&
               c.port == port &&
               c.readcount >= reuse_limit &&
               !c.readbusy &&
               isopen(c.io)), pool)
end

"""
    findall(type, host, port) -> Vector{Connection}

Find all `Connections` in the `pool` for `host` and `port`.
"""

function findall(T::Type,
                 host::AbstractString,
                 port::AbstractString,
                 pipeline_limit::Int)

    filter(c->(typeof(c.io) == T &&
               c.host == host &&
               c.port == port &&
               c.pipeline_limit == pipeline_limit &&
               isopen(c.io)), pool)
end

"""
    purge()

Remove closed connections from `pool`.
"""
function purge()
    isdeletable(c) = !isopen(c.io) && (@debug 1 "🗑  Deleted:        $c"; true)
    deleteat!(pool, map(isdeletable, pool))
end

"""
    getconnection(type, host, port) -> Connection

Find a reusable `Connection` in the `pool`,
or create a new `Connection` if required.
"""

function getconnection(::Type{Transaction{T}},
                       host::AbstractString,
                       port::AbstractString;
                       connection_limit::Int=default_connection_limit,
                       pipeline_limit::Int=default_pipeline_limit,
                       reuse_limit::Int=nolimit,
                       kw...)::Transaction{T} where T <: IO

    while true

        lock(poollock)
        try

            # Close connections that have reached the reuse limit...
            if reuse_limit != nolimit
                for c in findoverused(T, host, port, reuse_limit)
                    close(c)
                end
            end

            # Remove closed connections from `pool`...
            purge()

            # Try to find a connection with no active readers or writers...
            writable = findwritable(T, host, port, pipeline_limit, reuse_limit)
            idle = filter(c->!c.readbusy, writable)
            if !isempty(idle)
                c = rand(idle)                 ;@debug 2 "♻  Idle:           $c"
                return client_transaction(c)
            end

            # If there are not too many connections to this host:port,
            # create a new connection...
            busy = findall(T, host, port, pipeline_limit)
            if length(busy) < connection_limit
                io = getconnection(T, host, port; kw...)
                c = Connection(host, port, pipeline_limit, io)
                push!(pool, c)                ;@debug 1 "🔗  New:            $c"
                return client_transaction(c)
            end

            # Share a connection that has active readers...
            if !isempty(writable)
                c = rand(writable)             ;@debug 2 "⇆  Shared:         $c"
                return client_transaction(c)
            end

        finally
            unlock(poollock)
        end

        # Wait for `closewrite` or `close` to signal that a connection is ready.
        wait(poolcondition)
    end
end

function getconnection(::Type{TCPSocket},
                       host::AbstractString,
                       port::AbstractString;
                       kw...)::TCPSocket

    p::UInt = isempty(port) ? UInt(80) : parse(UInt, port)
    @debug 2 "TCP connect: $host:$p..."
    connect(getaddrinfo(host), p)
end

const nosslconfig = SSLConfig()
default_sslconfig = nothing
noverify_sslconfig = nothing

function global_sslconfig(require_ssl_verification::Bool)::SSLConfig
    global default_sslconfig
    global noverify_sslconfig
    if default_sslconfig == nothing
        default_sslconfig = SSLConfig(true)
        noverify_sslconfig = SSLConfig(false)
    end
    return require_ssl_verification ? default_sslconfig : noverify_sslconfig
end

function getconnection(::Type{SSLContext},
                       host::AbstractString,
                       port::AbstractString;
                       require_ssl_verification::Bool=true,
                       sslconfig::SSLConfig=nosslconfig,
                       kw...)::SSLContext

    if sslconfig === nosslconfig
        sslconfig = global_sslconfig(require_ssl_verification)
    end

    port = isempty(port) ? "443" : port
    @debug 2 "SSL connect: $host:$port..."
    io = SSLContext()
    setup!(io, sslconfig)
    associate!(io, getconnection(TCPSocket, host, port))
    hostname!(io, host)
    handshake!(io)
    return io
end

function Base.show(io::IO, c::Connection)
    nwaiting = nb_available(tcpsocket(c.io))
    print(
        io,
        tcpstatus(c), " ",
        lpad(c.writecount,3),"↑", c.writebusy ? "🔒  " : "   ",
        lpad(c.readcount,3), "↓", c.readbusy ? "🔒   " : "    ",
        c.host, ":",
        c.port != "" ? c.port : Int(c.peerport), ":", Int(c.localport),
        " ≣", c.pipeline_limit,
        length(c.excess) > 0 ? " $(length(c.excess))-byte excess" : "",
        inactiveseconds(c) > 5 ?
            " inactive $(round(inactiveseconds(c),1))s" : "",
        nwaiting > 0 ? " $nwaiting bytes waiting" : "",
        DEBUG_LEVEL > 1 ? " $(Base._fd(tcpsocket(c.io)))" : "")
end

Base.show(io::IO, t::Transaction) = print(io, "T$(rpad(t.sequence,2)) ", t.c)

function tcpstatus(c::Connection)
    s = Base.uv_status_string(tcpsocket(c.io))
        if s == "connecting" return "🔜🔗"
    elseif s == "open"       return "🔗 "
    elseif s == "active"     return "🔁 "
    elseif s == "paused"     return "⏸ "
    elseif s == "closing"    return "🔜💀"
    elseif s == "closed"     return "💀 "
    else
        return s
    end
end

function showpool(io::IO)
    lock(poollock)
    println(io, "ConnectionPool[")
    for c in pool
        println(io, "   $c")
    end
    println(io, "]\n")
    unlock(poollock)
end

function showpoolhtml(io::IO)
    lock(poollock)
    println(io, "
")
    for c in pool
        print(io, "
")
        for x in split("$c")
            print(io, "$x	")
        end
        println(io, "
")
    end
    println(io, "
")
    unlock(poollock)
end

end # module ConnectionPool


text:
function sphcesaro1(A::AbstractMatrix{T}) where T
    M, N = size(A)
    m, n = M, N÷2
    C = zeros(T, m, 2n+1)
    @inbounds for i = 1:m
        C[i,1] = (m+one(T)-i)/m
    end
    @inbounds for j = 1:n
        for i = 1:m-j
            C[i,2j] = C[i+j,1]
            C[i,2j+1] = C[i+j,1]
        end
    end
    C
end

function sphcesaro2(A::AbstractMatrix{T}) where T
    M, N = size(A)
    m, n = M, N÷2
    C = zeros(T, m, 2n+1)
    @inbounds for i = 1:m
        C[i,1] = (m+2one(T)-i)*(m+one(T)-i)/((m+one(T))*m)
    end
    @inbounds for j = 1:n
        for i = 1:m-j
            C[i,2j] = C[i+j,1]
            C[i,2j+1] = C[i+j,1]
        end
    end
    C
end

function sphcesaro(A::AbstractMatrix{T}, κ::Int) where T
    M, N = size(A)
    m, n = M, N÷2
    C = zeros(T, m, 2n+1)
    @inbounds for i = 1:m
        C[i,1] = binomial(m+κ-i, m-i)/binomial(m+κ-1, m-1)
    end
    @inbounds for j = 1:n
        for i = 1:m-j
            C[i,2j] = C[i+j,1]
            C[i,2j+1] = C[i+j,1]
        end
    end
    C
end

ext: jl


text:
function init(f, xsplits, lower::AbstractVector, upper::AbstractVector)
    # Validate the inputs
    n = length(lower)
    length(upper) == length(xsplits) == n || throw(DimensionMismatch("inconsistent dimensionality among lower, upper, and splits"))
    for i = 1:n
        xsi = xsplits[i]
        length(xsi) == 3 || throw(DimensionMismatch("3 values along each dimension required, got $(xsi) along dimension $i"))
        lower[i] <= xsi[1] < xsi[2] < xsi[3] <= upper[i] || error("splits must be in strictly increasing order and lie between lower and upper, got $xsi along dimension $i")
        all(isfinite, xsi) || error("not all splits along dimension $i are finite")
    end
    # Compute a consistent eltype
    xs1 = xsplits[1]
    Tx = promote_type(typeof(xs1[1]), typeof(xs1[2]), typeof(xs1[3]))
    for i = 2:n
        xsi = xsplits[i]
        Tx = promote_type(Tx, typeof(xsi[1]), typeof(xsi[2]), typeof(xsi[3]))
    end
    Tx = boxeltype(Tx)
    x0 = Tx[x[2] for x in xsplits]
    box, xstar = _init(f, copy(x0), xsplits, lower, upper)
    box, x0, xstar
end

@noinline function _init(f, xstar, xsplits, lower, upper)
    T = boxeltype(promote_type(eltype(xstar), eltype(lower), eltype(upper)))
    n = length(xstar)
    root = box = Box{T,n}()
    xtmp = copy(xstar)
    fcur = dummyvalue(T)
    for i = 1:n
        box = split!(box, f, xtmp, i, xsplits[i], lower[i], upper[i], xstar[i], fcur)
        xcur, fcur = box.parent.xvalues[box.parent_cindex], box.parent.fvalues[box.parent_cindex]
        xstar = replacecoordinate!(xstar, i, xcur)
    end
    box, xstar
end

function split!(box::Box{T}, f, xtmp, splitdim, xsplit, lower::Real, upper::Real, xcur, fcur) where T
    # Evaluate f along the splits, keeping track of the best
    fsplit = MVector3{T}(typemax(T), typemax(T), typemax(T))
    fmin, idxmin = typemax(T), 0
    for l = 1:3
        @assert(isfinite(xsplit[l]))
        if xsplit[l] == xcur && !isdummy(fcur)
            ftmp = fcur
        else
            xtmp = replacecoordinate!(xtmp, splitdim, xsplit[l])
            ftmp = f(xtmp)
        end
        if ftmp < fmin
            fmin = ftmp
            idxmin = l
        end
        fsplit[l] = ftmp
    end
    idxmin == 0 && error("function was not finite at any evaluation point $xsplit")
    if idxmin == 1 && fsplit[1] == fsplit[2]
        idxmin = 2  # prefer the middle in case of ties
    end
    add_children!(box, splitdim, xsplit, fsplit, lower, upper)
    xtmp = replacecoordinate!(xtmp, splitdim, xsplit[idxmin])
    return box.children[idxmin]
end

# Splits a box. If the "tilt" over the domain of the
# box suggests that one of its neighbors might be even lower,
# recursively calls itself to split that box, too.
# Returns `box, minval`.
function autosplit!(box::Box{T}, mes::Vector{<:MELink}, f::WrappedFunction, x0, xtmp, splitdim, xsplitdefaults, lower, upper, minwidth, visited::Set) where T
    box ∈ visited && error("already visited box")
    box.qnconverged && return (box, value(box))
    if !isleaf(box)
        # This box already got split along a different dimension
        box = find_smallest_child_leaf(box)
        xtmp = position(box, x0)
    end
    if splitdim == 0
        # If we entered this box from a neighbor, just choose an axis that has been split the least
        nsplits = count_splits(box)
        splitdim = argmin(nsplits)
        if nsplits[splitdim] > 0
            # If all dimensions have been split, prefer splitting unbounded dimensions
            bbs = boxbounds(box, lower, upper)
            if all(isfinite, bbs[splitdim])
                for (i, bb) in enumerate(bbs)
                    if !all(isfinite, bb)
                        splitdim = i
                        break
                    end
                end
            end
        end
    end
    xsplitdefault = xsplitdefaults[splitdim]
    lwr, upr = lower[splitdim], upper[splitdim]
    p = find_parent_with_splitdim(box, splitdim)
    if isroot(p) && p.splitdim != splitdim
        xcur, fcur = x0[splitdim], box.parent.fvalues[box.parent_cindex]
        split!(box, f, xtmp, splitdim, xsplitdefault, lwr, upr, xcur, fcur)
        trimschedule!(mes, box, splitdim, x0, lower, upper)
        return box, minimum(box.fvalues)
    end
    bb = boxbounds(p)
    bb[2]-bb[1] > max(minwidth[splitdim], epswidth(bb)) || return (box, value(box))
    xp, fp = p.parent.xvalues, p.parent.fvalues
    Δx = max(xp[2]-xp[1], xp[3]-xp[2])  # a measure of the "pragmatic" box scale even when the box is infinite
    xcur = xp[p.parent_cindex]
    @assert(xcur == xtmp[splitdim])
    fcur = box.parent.fvalues[box.parent_cindex]
    # Do a quadratic fit along splitdim
    xvert, fvert, qcoef = qfit(xp[1]=>fp[1], xp[2]=>fp[2], xp[3]=>fp[3])
    if qcoef > 0
        # xvert is a minimum. However, in general we can't trust it: for example, for the
        # "canyon" function
        #     f(x,y) = (x+y)^2/k + k*(x-y)^2,
        # at fixed `x` the minimum `y` is at
        #     ystar = (k^2-1)/(k^2+1)*x.
        # If the current box has a different `x` than the one used for box `p`,
        # then this estimate will be wrong.
        # (If `p == box.parent` it would be safe, but given that we'd have to pick a third
        # point anyway, this case does not seem worth treating separately.)
        # Instead of trusting `xvert`, we take `qcoef` (the quadratic coefficient) at face value. We then
        # choose a second (arbitrary) point along splitdim (e.g., a new `y` in the "canyon" example)
        # to evaluate `f`, and then use these two position/value pairs to estimate the
        # remaining two parameters of the quadratic.
        if isinf(bb[1])
            xnew = xcur - 8*(bb[2]-xcur)
        elseif isinf(bb[2])
            xnew = xcur + 8*(xcur-bb[1])
        else
            xnew = bb[2]-xcur > xcur-bb[1] ? (xcur+2*bb[2])/3 : (2*bb[1]+xcur)/3
        end
        @assert(isfinite(xnew))
        xtmp = replacecoordinate!(xtmp, splitdim, xnew)
        fnew = f(xtmp)
        # Solve for the minimum of the quadratic given the two pairs xcur=>fcur, xnew=>fnew
        xvert = (xcur+xnew)/2 - (fnew-fcur)/(2*qcoef*(xnew-xcur))
        if bb[1] <= xvert <= bb[2]
            # xvert is in the box, use it---but first, make sure it's not "degenerate" with
            # xcur or xnew
            xvert = ensure_distinct(xvert, xcur, xnew, bb)
            xtmp = replacecoordinate!(xtmp, splitdim, xvert)
            fvert = f(xtmp)
            xf1, xf2, xf3 = order_pairs(xcur=>fcur, xnew=>fnew, xvert=>fvert)
            add_children!(box, splitdim, MVector3{T}(xf1[1], xf2[1], xf3[1]),
                          MVector3{T}(xf1[2], xf2[2], xf3[2]), lwr, upr)
            trimschedule!(mes, box, splitdim, x0, lower, upper)
            return box, minimum(box.fvalues)
        end
        # xvert is not in the box. Prepare to split the neighbor, but for this box
        # just bisect xcur and xnew
        if xvert < bb[1]
            xtmp, dir = replacecoordinate!(xtmp, splitdim, bb[1]), -1
        else
            xtmp, dir = replacecoordinate!(xtmp, splitdim, bb[2]), +1
        end
        nbr, success = find_leaf_at_edge(get_root(box), xtmp, splitdim, dir)
        xmid = (xcur+xnew)/2
        xtmp = replacecoordinate!(xtmp, splitdim, xmid)
        fmid = f(xtmp)
        xf1, xf2, xf3 = order_pairs(xcur=>fcur, xnew=>fnew, xmid=>fmid)
        add_children!(box, splitdim, MVector3{T}(xf1[1], xf2[1], xf3[1]),
                      MVector3{T}(xf1[2], xf2[2], xf3[2]), lwr, upr)
        trimschedule!(mes, box, splitdim, x0, lower, upper)
        minbox = minimum(box.fvalues)
        if !success || nbr ∈ visited || !isfinite(value(nbr)) || length(visited) > ndims(box)
            return (box, minbox) # check nbr ∈ visited to avoid cycles
        end
        nbr, mv = autosplit!(nbr, mes, f, x0, position(nbr, x0), 0, xsplitdefaults, lower, upper, minwidth, push!(visited, box))
        if mv < minbox
            return nbr, mv
        end
        return box, minbox
    end
    trisect!(box, f, xtmp, splitdim, bb, Δx, xcur, fcur)
    trimschedule!(mes, box, splitdim, x0, lower, upper)
    return box, isleaf(box) ? value(box) : minimum(box.fvalues)
end

function trisect!(box::Box{T}, f, xtmp, splitdim, bb, Δx, xcur, fcur) where T
    l, r = bb
    if isinf(l)
        l = r - 6*Δx  # gap is usually 1/3 of the box size, so undo this at the level of "box size"
    elseif isinf(r)
        r = l + 6*Δx
    end
    a, b, c = 5l/6+r/6, (l+r)/2, l/6+5r/6
    # Ensure xcur is one of the three
    if xcur < (a+b)/2
        a = xcur
    elseif xcur < (b+c)/2
        b = xcur
    else
        c = xcur
    end
    if a < b < c
        return split!(box, f, xtmp, splitdim, MVector3{T}(a, b, c), bb..., xcur, fcur)
    end
    return box
end

# This requires there to be a (non-root) parent split along splitdim
function trisect!(box::Box, f, xtmp, splitdim, p = find_parent_with_splitdim(box, splitdim))
    bb = boxbounds(p)
    xp, fp = p.parent.xvalues, p.parent.fvalues
    Δx = max(xp[2]-xp[1], xp[3]-xp[2])  # a measure of the "pragmatic" box scale even when the box is infinite
    xcur = xp[p.parent_cindex]
    @assert(xcur == xtmp[splitdim])
    fcur = box.parent.fvalues[box.parent_cindex]
    trisect!(box, f, xtmp, splitdim, bb, Δx, xcur, fcur)
end

# Returns target box, value, and α. Uses α=0 to signal failure.
function quasinewton!(box::Box{T}, B, g, c, scale, f::Function, x0, lower, upper, itermax = 20) where T
    x = position(box, x0)
    fbox = fmin = value(box)
    isfinite(fbox) || return box, fmin, T(0)
    # Solve the bounded linear least-squares problem forcing B to be positive-definite
    cB = cholesky(Positive, B)
    Δx = -(cB \ g)   # unconstrained solution
    lsc, usc = (lower.-x)./scale, (upper.-x)./scale
    Δx[:] .= clamp.(Δx, lsc, usc)
    Δx = lls_bounded!(Δx, Matrix(cB), g, lsc, usc)
    # Test whether the update moves the point appreciably. If not, mark this point
    # as converged.
    issame = true
    for i = 1:length(Δx)
        issame &= abs(Δx[i]) < sqrt(eps(scale[i]))
    end
    if issame
        box.qnconverged = true
        return box, fmin, T(0)
    end
    Δx[:] .= Δx .* scale

    α = T(1.0)
    iter = 0
    while !isinside(x + α*Δx, lower, upper) && iter < itermax
        α /= 2
        iter += 1
    end
    iter == itermax && return box, fmin, T(0)
    iter = 0 # the above weren't "real" iterations, so reset
    root = get_root(box)
    # Do a backtracking linesearch
    local xtarget, ftarget
    while iter < itermax
        iter += 1
        xtarget = x + α*Δx
        # Sadly, the following function evaluations get discarded. But recording them
        # while preserving the box structure would take ndims(box)-1 additional
        # evaluations, which is not worth it unless xtarget itself is an improvement.
        ftarget = f(xtarget)
        ftarget < fbox && break
        α /= 2
    end
    ftarget >= fbox && return box, fmin, T(0)

    leaf = find_leaf_at(root, xtarget)
    isfinite(value(leaf)) || return box, fmin, T(0)
    # The new point should improve on what was already obtained in `leaf`
    ftarget >= value(leaf) && return box, fmin, T(0)

    # For consistency with the tree structure, we can change only one coordinate of
    # xleaf per split. Cycle through all coordinates, picking the one at each stage
    # that yields the smallest value as predicted by the quadratic model among the
    # not-yet-selected coordinates.
    dims_targeted = falses(ndims(leaf))
    q(x) = (x'*B*x)/2 + g'*x + c
    for j = 1:ndims(leaf)
        xleaf = position(leaf, x0)
        xtest = copy(xleaf)
        imin, qmin = 0, typemax(T)
        for i = 1:ndims(leaf)
            dims_targeted[i] && continue
            xtest = replacecoordinate!(xtest, i, xtarget[i])
            qx = q(xtest)
            if qx < qmin
                qmin = qx
                imin = i
            end
        end
        imin == 0 && break
        xcur, xt = xleaf[imin], xtarget[imin]
        if abs(xcur - xt) < sqrt(eps(scale[imin]))
            # We're already at this point, continue
            dims_targeted[imin] = true
            continue
        end
        bb = boxbounds(leaf, imin, lower, upper)
        a, b, c = pick3(xcur, xt, bb)
        minsep = eps(T)*(c-a)
        if a + minsep >= b || b + minsep >= c
            # The box might be so narrow that there are not enough distinct floating-point
            # numbers that lie between the bounds.
            dims_targeted[imin] = true
            continue
        end
        split!(leaf, f, xleaf, imin, MVector3{T}(a, b, c), bb..., xleaf[imin], value(leaf))
        fmin = min(fmin, minimum(leaf.fvalues))
        childindex = findfirst(isequal(xt), leaf.xvalues)
        leaf = leaf.children[childindex]
        isfinite(value(leaf)) || return leaf, fmin, T(0)
        dims_targeted[imin] = true
    end
    return leaf, fmin, α
end

# Linear least squares with box-bounds. From
#    Sequential Coordinate-wise Algorithm for the
#    Non-negative Least Squares Problem
#    Vojtˇch e Franc, V ́clav aHlav ́ acˇ, and Mirko Navara
# Solves min(x'*H*x/2 + f'*x) over the domain bounded by lower and upper
function lls_bounded(H::AbstractMatrix{T}, f::VT, lower::V, upper::V, tol = sqrt(eps(T))) where {T,VT<:AbstractVector{T},V<:AbstractVector}
    x = clamp.(-H \ f, lower, upper)
    lls_bounded!(x, H, f, lower, upper, zeros(T, length(f)), tol)
    return x
end

lls_bounded!(x::VT, H::AbstractMatrix{T}, f::VT, lower::V, upper::V, tol = sqrt(eps(T))) where {T,VT<:AbstractVector{T},V<:AbstractVector} =
    lls_bounded!(x, H, f, lower, upper, zeros(T, length(f)), tol)

function lls_bounded!(x::VT, H::AbstractMatrix{T}, f::VT, lower::V, upper::V, g::VT, tol = sqrt(eps(T))) where {T,VT<:AbstractVector{T},V<:AbstractVector}
    assert_oneindexed(A) = @assert(all(r->first(r)==1, axes(A)))
    assert_oneindexed(x); assert_oneindexed(H); assert_oneindexed(f); assert_oneindexed(lower)
    assert_oneindexed(upper); assert_oneindexed(g)
    m, n = size(H)
    m == n || error("H must be square")
    length(f) == n || error("Mismatch between H and f")
    length(g) == n || error("Mismatch between H and g")
    niter = 0
    while niter <= 1000
        # Once per iter, calculate gradient freshly to avoid accumulation of roundoff
        mul!(g, H, x)
        for i = 1:n
            g[i] += f[i]
        end
        xchange = zero(T)
        for k = 1:n
            xold = x[k]
            xnew = clamp(xold - g[k]/H[k,k], lower[k], upper[k])
            x[k] = xnew
            dx = xnew-xold
            xchange = max(xchange, abs(dx))
            if dx != 0
                for i = 1:n
                    g[i] += dx*H[i,k]
                end
            end
        end
        if xchange < tol
            break
        end
        niter += 1
    end
    x
end

# A dumb O(N) algorithm for building the minimum-edge structures
# For large trees, this is the bottleneck
function minimum_edges(root::Box{T,N}, x0, lower, upper, minwidth=zeros(eltype(x0), ndims(root)); extrapolate::Bool=true) where {T,N}
    mes = [MELink{T,T}(root) for i = 1:N]
    # Make it a priority to have enough splits along each coordinate to provide adequate
    # information for the quasi-Newton method. To compute an estimate of the full quadratic
    # model, we'll need (N+1)*(N+2)÷2 independent points, so prioritze splitting
    # boxes so that each coordinate has at least 1/Nth of the requisite number of splits.
    nthresh = ceil(Int, (((N+1)*(N+2))÷2)/N)
    nsplits = Vector{Int}(undef, N)
    for box in leaves(root)
        box.qnconverged && continue
        fval = value(box)
        isfinite(fval) || continue
        count_splits!(nsplits, box)
        nmin = minimum(nsplits)
        for i = 1:N
            if nmin < nthresh
                nsplits[i] == nmin || continue
            end
            bb = boxbounds(find_parent_with_splitdim(box, i), lower[i], upper[i])
            bb[2]-bb[1] < minwidth[i] && continue
            insert!(mes[i], width(box, i, x0, lower, upper), box=>extrapolate ? fval+qdelta(box, i) : fval)
        end
    end
    return mes
end

function trimschedule!(mes::Vector{<:MELink}, box::Box, splitdim, x0, lower, upper)
    isleaf(box) && return mes
    for child in box.children
        fval = child.parent.fvalues[child.parent_cindex]
        for i = 1:ndims(box)
            trim!(mes[i], width(child, i, x0, lower, upper), child=>fval)
        end
    end
    return mes
end

function mesprint(mes)
    for i = 1:length(mes)
        print(i, ": ")
        for item in mes[i]
            print(", (", item.w, ',', item.f, ')')
        end
        println()
    end
    nothing
end

function sweep!(root::Box, f, x0, splits, lower, upper; extrapolate::Bool = true, fvalue=-Inf, minwidth=zeros(eltype(x0), ndims(root)))
    mes = minimum_edges(root, x0, lower, upper, minwidth; extrapolate=extrapolate)
    sweep!(root, mes, f, x0, splits, lower, upper; fvalue=fvalue, minwidth=minwidth)
end
function sweep!(root::Box{T}, mes::Vector{<:MELink}, f, x0, splits, lower, upper; fvalue=-Inf, minwidth=zeros(eltype(x0), ndims(root))) where T
    xtmp = similar(x0)
    flag = similar(x0, Bool)
    nsplits = similar(x0, Int)
    visited = Set{typeof(root)}()
    splitboxes = Set{typeof(root)}()
    dimorder = sortperm(length.(mes))  # process the dimensions with the shortest queues first
    fvalueT = T(fvalue)
    for i in dimorder
        me = mes[i]
        while !isempty(me)
            item = popfirst!(me)
            box = item.l
            box.qnconverged && continue
            bb = boxbounds(find_parent_with_splitdim(box, i), lower[i], upper[i])
            bb[2]-bb[1] < minwidth[i] && continue
            position!(xtmp, flag, box)
            default_position!(xtmp, flag, x0)
            count_splits!(nsplits, box)
            if nsplits[i] > 0 && any(iszero, nsplits)
                # println("discarding ", box, " along dimension ", i)
                continue
            end
            empty!(visited)
            finalbox, minval = autosplit!(box, mes, f, x0, xtmp, i, splits, lower, upper, minwidth, visited)
            push!(splitboxes, finalbox)
            minval <= fvalueT && return root, unique_smallest_leaves(splitboxes)
        end
    end
    root, unique_smallest_leaves(splitboxes)
end

"""
    root, x0 = analyze(f, splits, lower, upper; rtol=1e-3, atol=0.0, fvalue=-Inf, maxevals=2500, nquasinewton=<auto>, minwidth=zeros(length(splits)), print_interval=typemax(Int))

Analyze the behavior of `f`, searching for minima, over the rectangular box specified by
`lower` and `upper` (`lower[i] <= x[i] <= upper[i]`). The bounds may be infinite.
`splits` is a list of 3-vectors containing the initial values along each coordinate
axis at which to potentially evaluate `f`; the values must be in increasing order.

`root` is a tree structure that stores information about the behavior of `f`.
`x0` contains the initial evaluation point, the position constructed from the middle value
of `splits` in each dimension.

There are several keywords which can be important:

- `rtol` and `atol` represent relative and absolute, respectively, changes in minimum
  function value required for the exploration to terminate. (These limits must be hit on
  `ndims` successive sweeps.)
  Alternatively, the analysis is terminated if the function
  value is ever reduced below `fvalue`. The combination `rtol=0, fvalue=<user-supplied-value>`
  can be especially useful in benchmarking situations in which you know the value of the
  global minimum; it will force the algorithm to keep trying until it reaches `fvalue`,
  or until it exceeds the maximum number of allowed function evaluations.
- `maxevals` allows you to control the maximum number of function evaluations that the
  algorithm is permitted to use. The actual number may be a bit higher, as the criterion
  is not checked while, e.g., in the middle of a splitting-sweep or a quasi-Newton step.
- `nquasinewton` is the number of "explore" function evaluations required before attempting
  the first quasi-Newton step. The default value is the minimum number of evaluations
  required to estimate the parameters of a quadratic model, `(N+1)*(N+2)÷2` in `N` dimensions.
  Using a higher number may encourage the algorithm to spend more time "exploring" before
  thoroughly "exploiting" the most promising leads so far.
- `minwidth[i]` sets the minimum box size along dimension `i`. This is not strictly enforced,
  but boxes that become smaller will not be split during "explore" (sweep) phases of the algorithm.
  You can use this to improve coverage of space when you know that some parameter-value
  changes are too small to lead to measurable improvements.
- setting `print_interval` displays a small amount of information about progress.

# Example:

    # A function that has a long but skinny valley aligned at 45 degrees (minimum at [0,0])
    function canyon(x)
        x1, x2 = x[1], x[2]
        return 0.1*(x1+x2)^2 + 10*(x1-x2)^2
    end

    lower, upper = [-Inf, -Inf], [Inf, Inf]  # unbounded domain

    # We'll initially consider a grid that covers -11->-9 along the first dimension
    # and -7->-5 along the second. This isn't a great initial guess, but that's OK.
    splits = ([-11,-10,-9], [-7,-6,-5])

    # Since this problem has a minimum value of 0, relying on `rtol` is not ideal
    # (it's possible to make infinitesimal relative improvements nearly forever),
    # so specify an absolute tolerance.
    julia> root, x0 = analyze(canyon, splits, lower, upper; atol=0.01)
    (BoxRoot@[NaN, NaN], [-10.0, -6.0])

    julia> box = minimum(root)
    Box1.9407745316864587e-25@[-6.94556e-13, -6.98108e-13]

    julia> value(box)
    1.9407745316864587e-25

    julia> position(box, x0)
    2-element Array{Float64,1}:
     -6.94556e-13
     -6.98108e-13

See also [`minimize`](@ref).
"""
function analyze(f, splits, lower, upper; kwargs...)
    box, x0 = init(f, splits, lower, upper)
    root = get_root(box)
    analyze!(root, f, x0, splits, lower, upper; kwargs...)
    return root, x0
end

"""
    root = analyze!(root, f, x0, splits, lower, upper; kwargs...)

Further refinement of `root`. See [`analyze`](@ref) for details.
"""
function analyze!(root::Box, f::Function, x0, splits, lower, upper; rtol=1e-3, atol=0.0, fvalue=-Inf, maxevals=2500, print_interval=typemax(Int), kwargs...)
    analyze!(root, CountedFunction(f), x0, splits, lower, upper; rtol=rtol, atol=atol, fvalue=fvalue, maxevals=maxevals, print_interval=print_interval, kwargs...)
end
function analyze!(root::Box{T}, f::WrappedFunction, x0, splits, lower, upper; rtol=1e-3, atol=0.0, fvalue=-Inf, maxevals=2500, print_interval=typemax(Int), nquasinewton=qnthresh(ndims(root)), kwargs...) where T
    box = minimum(root)
    boxval = value(box)
    lastval = typemax(boxval)
    tol_counter = 0
    lastprint = 0
    baseline_evals = len = lenold = length(leaves(root))
    if baseline_evals == numevals(f)
        baseline_evals = 0   # already counted
    end
    if print_interval < typemax(Int)
        println("Initial minimum ($len evaluations): ", box)
    end
    extrapolate = true
    # The quasi-Newton step can reduce the function value so significantly, insist on
    # using it at least once (unless the threshold to kick in not reachable)
    used_quasinewton = nquasinewton > maxevals
    nsweep, nqn = baseline_evals+numevals(f), 0
    thresh = sqrt(eps(T))
    nsplits = zeros(Int, ndims(root))
    while boxval > fvalue && (tol_counter <= ndims(root) || !used_quasinewton) && len < maxevals
        lastval = boxval
        numeval0 = numevals(f)
        _, splitboxes = sweep!(root, f, x0, splits, lower, upper; extrapolate=extrapolate, fvalue=fvalue, kwargs...)
        box = minimum(root)
        numeval1 = numevals(f)
        nsweep += numeval1 - numeval0
        if nquasinewton < baseline_evals + numevals(f) < maxevals && boxval > fvalue && !isempty(splitboxes)
            # Quasi-Newton step for boxes split in the sweep
            spbxs = sort(collect(splitboxes); by=value) # if unsorted, set hashing would introduce randomness
            # determine which are in separate "basins" using a convexity test
            # basinboxes = spbxs
            # basinboxes = different_basins(spbxs, x0, lower, upper)
            basinboxes = filter(is_diag_convex, spbxs)
            badbox = Dict(box=>false for box in basinboxes)
            for (ibox, box) in enumerate(basinboxes)
                box.qnconverged && continue
                isleaf(box) || continue   # "outdated"
                badbox[box] && continue   # in same basin as a previous successful quadratic fit
                baseline_evals + numevals(f) >= maxevals && break
                # Ensure at least one split per coordinate before trying to build the quadratic model
                count_splits!(nsplits, box)
                xtmp = position(box, x0)
                for i = 1:ndims(box)
                    if nsplits[i] == 0
                        box = split!(box, f, xtmp, i, splits[i], lower[i], upper[i], xtmp[i], value(box))
                        xtmp[i] = box.parent.xvalues[box.parent_cindex]
                    end
                end
                scale = boxscale(box, splits)
                Q, xbase, c, success_build = build_quadratic_model(box, x0, scale, thresh)
                if success_build && all(x->x>0, Q.dimpiv)
                    # Add more points, as needed, to fill in any missing rows from the regression
                    @assert isleaf(box)
                    complete_quadratic_model!(Q, c, box, f, x0, xbase, scale, splits, lower, upper, thresh)
                    # If the model is complete, try Quasi-Newton algorithm
                    if Q.nzrows[] == size(Q.coefs, 1) && minimum(abs.(diag(Q.coefs))) >= thresh
                        g, B = solve(Q)
                        # Bscale = (1./scale) .* B .* (1./scale)'
                        # display(Bscale)
                        # error("stop")
                        leaf, minvalue, α = quasinewton!(box, B, g, c, scale, f, x0, lower, upper)
                        success_qn = α > 0
                        used_quasinewton |= success_qn
                        minvalue < fvalue && break
                        if success_qn
                            @assert(isleaf(leaf))
                            boxvalue = value(box)
                            # If the improvement was tiny, mark the solution to block
                            # future attempts at quasi-Newton in this region
                            if boxvalue - minvalue <= max(atol, rtol*(abs(boxvalue) + abs(minvalue)))
                                # Check that leaf has roughly the same coordinates, too.
                                xbox = position(box, x0)
                                xleaf = position(leaf, x0)
                                if issame(xbox, xleaf, scale)
                                    leaf.qnconverged = true
                                end
                            end
                            # check future boxes and eliminate them if they appear to lie in the same basin
                            if α == 1
                                for j = ibox+1:length(basinboxes)
                                    boxj = basinboxes[j]
                                    vj, xj = value(boxj), position(boxj, x0)
                                    Δx = (xj - xbase) ./ scale
                                    pred = c + g'*Δx + (Δx'*B*Δx)/2
                                    if abs(pred - vj) <= 0.1*abs(pred - minvalue)   # FIXME hardcoded
                                        badbox[boxj] = true
                                    end
                                end
                            end
                        end
                    end
                end
            end
            # if !used_quasinewton
            #     nquasinewton *= 2
            # end
            box = minimum(root)
            boxval = value(box)
        end
        nqn += numevals(f) - numeval1
        len = baseline_evals + numevals(f)
        len == lenold && break  # couldn't split any boxes
        lenold = len
        if len-lastprint > print_interval
            println("minimum ($len evaluations): ", box)
            lastprint = len
        end
        @assert(boxval <= lastval)
        if lastval - boxval < atol || lastval - boxval < rtol*(abs(lastval) + abs(boxval))
            tol_counter += 1
        else
            tol_counter = 0
        end
    end
    if print_interval < typemax(Int)
        println("Final minimum ($len evaluations): ", minimum(root))
        println("$nsweep evaluations for initialization+sweeps, $nqn for Quasi-Newton")
    end
    root
end

"""
    xmin, fmin = minimize(f, splits, lower, upper; kwargs...)

Return the position `xmin` and value `fmin` of the minimum of `f` over the specified domain.
See [`analyze`](@ref) for information about the input arguments.
"""
function minimize(f, splits, lower, upper; kwargs...)
    root, x0 = analyze(f, splits, lower, upper; kwargs...)
    box = minimum(root)
    return position(box, x0), value(box)
end


text:
using ArtifactUtils
using Artifacts

add_artifact!(
           "Artifacts.toml",
           "JuliaMono",
           "https://github.com/cormullion/juliamono/releases/download/v0.043/JuliaMono-ttf.tar.gz",
           force=true,
       )

import Pkg; Pkg.instantiate()


text:
using Pkg
current_env = Base.load_path()[1]
Pkg.activate(@__DIR__)

using Literate

example_dir = joinpath(@__DIR__, "..", "examples")
src_dir = joinpath(@__DIR__, "..", "src")

function replace_comments( content )
    content = replace(content, r"^(\h*)#~(.*)$"m => s"\1## \2")
    return content
end

#%% Functions to convert notebook
using PlutoSliderServer
function make_notebook_html( notebook_jl_path; execute = true )
    out_dir = joinpath(@__DIR__, "src", "custom_assets")
    if execute
        PlutoSliderServer.export_notebook( notebook_jl_path; Export_output_dir = out_dir )
    end
    base_name = splitpath(splitext( notebook_jl_path )[1])[end]
    return joinpath( out_dir, base_name * ".html" )
end

function modify_notebook_html( html_name )
    content = open( html_name, "r" ) do html_file
        read( html_file, String )
    end
    content = replace(content,
        r"(<script\b[^>]*>[\s\S]*?<\/script>)" =>
        s"\1\n\t<script src='./iframeResizer.contentWindow.min.js'></script>\n";
        count = 1
    )
    open( html_name, "w" ) do html_file
        write( html_file, content )
    end
    return nothing
end

const md_nb_template = """
```@raw html
<iframe id="fdnotebook" src="../custom_assets/HTML_NAME" width="100%"></iframe>
<!--<script src="../custom_assets/iframeResizer.min.js"></srcipt>-->
<script>
const iFrameResizerPath = '../custom_assets/iframeResizer.min.js';

if (require) {
  require([iFrameResizerPath], (iFrameResize) => iFrameResize())
} else {
  const script = document.createElement('script')
  script.onload = () => iFrameResize()
  script.src = iFrameResizerPath
}
</script>
<script>
document.addEventListener('DOMContentLoaded', function(){
	var myIframe = document.getElementById("fdnotebook");
	iFrameResize({log:true}, myIframe);
});
</script>
```
"""

function create_md_nb_file( html_name )
    global md_nb_template;
    html_filename = splitpath(html_name)[end]
    content = replace( md_nb_template, "HTML_NAME" => html_filename)
    md_path = joinpath(@__DIR__, "src", splitext( html_filename )[1] * ".md" )
    open(md_path, "w") do md_file
        write(md_file, content)
    end
    return md_path
end

function make_notebook_md( notebook_jl_path; execute = true )
    html_path = make_notebook_html( notebook_jl_path; execute )
    modify_notebook_html( html_path )
    return create_md_nb_file( html_path )
end

#%%
Literate.markdown(
    joinpath( example_dir, "example_two_parabolas.jl"),
    joinpath( @__DIR__, "src" );
    preprocess = replace_comments,
#    codefence = "````julia" => "````",  # disables execution; REMOVE when it works again
    )

Literate.markdown(
    joinpath( example_dir, "example_zdt.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````", # disables execution; REMOVE when it works again
    preprocess = replace_comments
    )

Literate.markdown(
    joinpath( src_dir, "TaylorModel.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````",
)

Literate.markdown(
    joinpath( src_dir, "ExactModel.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````",
)

Literate.markdown(
    joinpath( src_dir, "LagrangeModel.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````",
)

Literate.markdown(
    joinpath( src_dir, "RbfModel.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````",
)

Literate.markdown(
    joinpath( src_dir, "custom_logging.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````",
)

Literate.markdown(
    joinpath( src_dir, "Interfaces.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````",
)

Literate.markdown(
    joinpath( src_dir, "AbstractResultInterface.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````",
)

Literate.markdown(
    joinpath( src_dir, "AbstractIterDataInterface.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````",
)

Literate.markdown(
    joinpath( src_dir, "AbstractDBInterface.jl"),
    joinpath( @__DIR__, "src" );
    codefence = "````julia" => "````",
)
#%%

#make_notebook_md( joinpath( example_dir, "notebook_finite_differences.jl" ); execute = true )
#make_notebook_md( joinpath( example_dir, "notebook_polynomial_interpolation.jl" ); execute = true )

#%%
Pkg.activate(current_env)


text:
DIR=pwd()*"/cuda_kernels/"

println("loading CUDA kernels from $DIR")

md=CuModule(DIR*"vcopyshift.ptx",false)
vcopyshift_kernel=CuFunction(md,"vcopyshift")
function copyinto!(out::CudaArray{Float64},A::CudaArray{Float64},location::Int)
    nblocks=round(Int,ceil(length(A)/1024))
    launch(vcopyshift_kernel,nblocks,1024,(length(A),location-1,A,out))
end

md=CuModule(DIR*"vcopyshift_32.ptx",false)
vcopyshift_kernel_32=CuFunction(md,"vcopyshift_32")
function copyinto!(out::CudaArray{Float32},A::CudaArray{Float32},location::Int)
    nblocks=round(Int,ceil(length(A)/1024))
    launch(vcopyshift_kernel_32,nblocks,1024,(length(A),location-1,A,out))
end
export copyinto!

md=CuModule(DIR*"vcopyfrom.ptx",false)
vcopyfrom_kernel=CuFunction(md,"vcopyfrom")
function copyfrom!(out::CudaArray{Float64},A::CudaArray{Float64},location::Int)
    nblocks=round(Int,ceil(length(out)/1024))
    launch(vcopyfrom_kernel,nblocks,1024,(length(out),location-1,A,out))
end

md=CuModule(DIR*"vcopyfrom_32.ptx",false)
vcopyfrom_kernel_32=CuFunction(md,"vcopyfrom_32")
function copyfrom!(out::CudaArray{Float32},A::CudaArray{Float32},location::Int)
    nblocks=round(Int,ceil(length(out)/1024))
    launch(vcopyfrom_kernel_32,nblocks,1024,(length(out),location-1,A,out))
end
export copyfrom!

md=CuModule(DIR*"vcopyfrom_update.ptx",false)
vcopyfrom_update_kernel=CuFunction(md,"vcopyfrom_update")
function copyfrom_update!(out::CudaArray{Float64},A::CudaArray{Float64},location::Int)
    nblocks=round(Int,ceil(length(out)/1024))
    launch(vcopyfrom_update_kernel,nblocks,1024,(length(out),location-1,A,out))
end

md=CuModule(DIR*"vcopyfrom_update_32.ptx",false)
vcopyfrom_update_kernel_32=CuFunction(md,"vcopyfrom_update_32")
function copyfrom_update!(out::CudaArray{Float32},A::CudaArray{Float32},location::Int)
    nblocks=round(Int,ceil(length(out)/1024))
    launch(vcopyfrom_update_kernel_32,nblocks,1024,(length(out),location-1,A,out))
end
export copyfrom_update!

md=CuModule(DIR*"vsign.ptx",false)
vsign_kernel=CuFunction(md,"vsign")
function vsign!(A::CudaArray{Float64},out::CudaArray{Float64})
    nblocks=round(Int,ceil(length(out)/1024))
    launch(vsign_kernel,nblocks,1024,(length(out),A,out))
end

md=CuModule(DIR*"vsign_32.ptx",false)
vsign_kernel_32=CuFunction(md,"vsign_32")
function vsign!(A::CudaArray{Float32},out::CudaArray{Float32})
    nblocks=round(Int,ceil(length(out)/1024))
    launch(vsign_kernel_32,nblocks,1024,(length(out),A,out))
end
export vsign!

md=CuModule(DIR*"xsigny_update.ptx",false)
xsigny_update_kernel=CuFunction(md,"xsigny_update")
function xsigny_update!(X::CudaArray{Float64},Y::CudaArray{Float64},out::CudaArray{Float64})
    nblocks=round(Int,ceil(length(out)/1024))
    launch(xsigny_update_kernel,nblocks,1024,(length(out),X,Y,out))
end

md=CuModule(DIR*"xsigny_update_32.ptx",false)
xsigny_update_kernel_32=CuFunction(md,"xsigny_update_32")
function xsigny_update!(X::CudaArray{Float32},Y::CudaArray{Float32},out::CudaArray{Float32})
    nblocks=round(Int,ceil(length(out)/1024))
    launch(xsigny_update_kernel_32,nblocks,1024,(length(out),X,Y,out))
end
export xsigny_update!

md=CuModule(DIR*"vabs.ptx",false)
vabs_kernel=CuFunction(md,"vabs")
function abs!(A::CudaArray{Float64},out::CudaArray{Float64})
    nblocks=round(Int,ceil(length(out)/1024))
    launch(vabs_kernel,nblocks,1024,(length(out),A,out))
end
export abs!

md=CuModule(DIR*"gfill.ptx",false)
gfill_kernel=CuFunction(md,"gfill")
function gfill!(out::CudaArray{Float64},fillval::CudaArray{Float64})
    nblocks=round(Int,ceil(length(out)/1024))
    launch(gfill_kernel,nblocks,1024,(length(out),fillval,out))
end

md=CuModule(DIR*"gfill_32.ptx",false)
gfill_kernel_32=CuFunction(md,"gfill_32")
function gfill!(out::CudaArray{Float32},fillval::CudaArray{Float32})
    nblocks=round(Int,ceil(length(out)/1024))
    launch(gfill_kernel_32,nblocks,1024,(length(out),fillval,out))
end
export gfill!

md=CuModule(DIR*"tx1mx_32.ptx",false)
tx1mx_kernel_32=CuFunction(md,"tx1mx_32")
function tx1mx!(t::CudaArray{Float32},x::CudaArray{Float32},out::CudaArray{Float32})
    nblocks=round(Int,ceil(length(x)/1024))
    launch(tx1mx_kernel_32,nblocks,1024,(length(x),t,x,out))
end

md=CuModule(DIR*"tx1mx.ptx",false)
tx1mx_kernel=CuFunction(md,"tx1mx")
function tx1mx!(t::CudaArray{Float64},x::CudaArray{Float64},out::CudaArray{Float64})
    nblocks=round(Int,ceil(length(x)/1024))
    launch(tx1mx_kernel,nblocks,1024,(length(x),t,x,out))
end
export tx1mx!

md=CuModule(DIR*"Dstanh.ptx",false)
Dstanh_kernel=CuFunction(md,"Dstanh")
function Dstanh!(sf,grad_c::CudaArray{Float64},f_c::CudaArray{Float64},grad_n::CudaArray{Float64})
    nblocks=round(Int,ceil(length(f_c)/1024))
    launch(Dstanh_kernel,nblocks,1024,(length(f_c),sf,grad_c,f_c,grad_n))
end

md=CuModule(DIR*"Dstanh_32.ptx",false)
Dstanh_kernel_32=CuFunction(md,"Dstanh_32")
function Dstanh!(sf,grad_c::CudaArray{Float32},f_c::CudaArray{Float32},grad_n::CudaArray{Float32})
    nblocks=round(Int,ceil(length(f_c)/1024))
    launch(Dstanh_kernel_32,nblocks,1024,(length(f_c),sf,grad_c,f_c,grad_n))
end
export Dstanh!

md=CuModule(DIR*"DmeanSquareLoss.ptx",false)
DmeanSquareLoss_kernel=CuFunction(md,"DmeanSquareLoss")
function DmeanSquareLoss!(grad_c::CudaArray{Float64},x::CudaArray{Float64},y::CudaArray{Float64},grad_n::CudaArray{Float64})
    nblocks=round(Int,ceil(length(x)/1024))
    launch(DmeanSquareLoss_kernel,nblocks,1024,(length(x),2.0/length(x),grad_c,x,y,grad_n))
end

md=CuModule(DIR*"DmeanSquareLoss_32.ptx",false)
DmeanSquareLoss_kernel_32=CuFunction(md,"DmeanSquareLoss_32")
function DmeanSquareLoss!(grad_c::CudaArray{Float32},x::CudaArray{Float32},y::CudaArray{Float32},grad_n::CudaArray{Float32})
    nblocks=round(Int,ceil(length(x)/1024))
    launch(DmeanSquareLoss_kernel_32,nblocks,1024,(length(x),Float32(2.0/length(x)),grad_c,x,y,grad_n))
end
export DmeanSquareLoss!

md=CuModule(DIR*"stanh.ptx",false)
stanh_kernel=CuFunction(md,"stanh")
function stanh!(sf,A::CudaArray{Float64},B::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(stanh_kernel,nblocks,1024,(length(A),sf,A,B))
end

md=CuModule(DIR*"stanh_32.ptx",false)
stanh_kernel_32=CuFunction(md,"stanh_32")
function stanh!(sf,A::CudaArray{Float32},B::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(stanh_kernel_32,nblocks,1024,(length(A),sf,A,B))
end
export stanh!

md=CuModule(DIR*"binaryentropy.ptx",false)
binaryentropy_kernel=CuFunction(md,"binaryentropy")
function binaryentropy(X::CudaArray{Float64},Y::CudaArray{Float64})
    tmp=CudaArray(Float64,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(binaryentropy_kernel,nblocks,1024,(length(X),X,Y,tmp))
    out=mean(tmp)
    free(tmp)
    return out
end

md=CuModule(DIR*"binaryentropy_32.ptx",false)
binaryentropy_kernel_32=CuFunction(md,"binaryentropy_32")
function binaryentropy(X::CudaArray{Float32},Y::CudaArray{Float32})
    tmp=CudaArray(Float64,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(binaryentropy_kernel_32,nblocks,1024,(length(X),X,Y,tmp))
    out=mean(tmp)
    free(tmp)
    return out
end
export binaryentropy

md=CuModule(DIR*"binaryentropyXsigmoidY.ptx",false)
binaryentropyXsigmoidY_kernel=CuFunction(md,"binaryentropyXsigmoidY")
function binaryentropyXsigmoidY(X::CudaArray{Float64},Y::CudaArray{Float64})
    tmp=CudaArray(Float64,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(binaryentropyXsigmoidY_kernel,nblocks,1024,(length(X),X,Y,tmp))
    out=mean(tmp)
    free(tmp)
    return out
end

md=CuModule(DIR*"binaryentropyXsigmoidY_32.ptx",false)
binaryentropyXsigmoidY_kernel_32=CuFunction(md,"binaryentropyXsigmoidY_32")
function binaryentropyXsigmoidY(X::CudaArray{Float32},Y::CudaArray{Float32})
    tmp=CudaArray(Float32,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(binaryentropyXsigmoidY_kernel_32,nblocks,1024,(length(X),X,Y,tmp))
    out=mean(tmp)
    free(tmp)
    return out
end
export binaryentropyXsigmoidY

function binaryentropyXsigmoidY!(X::CudaArray{Float64},Y::CudaArray{Float64},out)
    tmp=CudaArray(Float64,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(binaryentropyXsigmoidY_kernel,nblocks,1024,(length(X),X,Y,tmp))
    copy!(out,mean(tmp))
    free(tmp)
end

function binaryentropyXsigmoidY!(X::CudaArray{Float32},Y::CudaArray{Float32},out::CudaArray{Float32})
    tmp=CudaArray(Float32,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(binaryentropyXsigmoidY_kernel_32,nblocks,1024,(length(X),X,Y,tmp))
    copy!(out,mean(tmp))
    free(tmp)
end
export binaryentropyXsigmoidY!

md=CuModule(DIR*"DXbinaryentropy.ptx",false)
DXbinaryentropy_kernel=CuFunction(md,"DXbinaryentropy")
function DXbinaryentropy!(X::CudaArray{Float64},Y::CudaArray{Float64},T::CudaArray{Float64},Out::CudaArray{Float64})
    tmp=CudaArray(Float64,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(DXbinaryentropy_kernel,nblocks,1024,(length(X),X,Y,T,Out))
end

md=CuModule(DIR*"DXbinaryentropy_32.ptx",false)
DXbinaryentropy_kernel_32=CuFunction(md,"DXbinaryentropy_32")
function DXbinaryentropy!(X::CudaArray{Float32},Y::CudaArray{Float32},T::CudaArray{Float32},Out::CudaArray{Float32})
    tmp=CudaArray(Float32,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(DXbinaryentropy_kernel_32,nblocks,1024,(length(X),X,Y,T,Out))
end
export DXbinaryentropy!

md=CuModule(DIR*"DYbinaryentropy.ptx",false)
DYbinaryentropy_kernel=CuFunction(md,"DYbinaryentropy")
function DYbinaryentropy!(X::CudaArray{Float64},Y::CudaArray{Float64},T::CudaArray{Float64},Out::CudaArray{Float64})
    tmp=CudaArray(Float64,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(DYbinaryentropy_kernel,nblocks,1024,(length(X),X,Y,T,Out))
end

md=CuModule(DIR*"DYbinaryentropy_32.ptx",false)
DYbinaryentropy_kernel_32=CuFunction(md,"DYbinaryentropy_32")
function DYbinaryentropy!(X::CudaArray{Float32},Y::CudaArray{Float32},T::CudaArray{Float32},Out::CudaArray{Float32})
    tmp=CudaArray(Float32,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(DYbinaryentropy_kernel_32,nblocks,1024,(length(X),X,Y,T,Out))
end
export DYbinaryentropy!

md=CuModule(DIR*"DXbinaryentropyXsigmoidY.ptx",false)
DXbinaryentropyXsigmoidY_kernel=CuFunction(md,"DXbinaryentropyXsigmoidY")
function DXbinaryentropyXsigmoidY!(X::CudaArray{Float64},Y::CudaArray{Float64},T::CudaArray{Float64},Out::CudaArray{Float64})
    tmp=CudaArray(Float64,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(DXbinaryentropyXsigmoidY_kernel,nblocks,1024,(length(X),X,Y,T,Out))
end

md=CuModule(DIR*"DXbinaryentropyXsigmoidY_32.ptx",false)
DXbinaryentropyXsigmoidY_kernel32=CuFunction(md,"DXbinaryentropyXsigmoidY_32")
function DXbinaryentropyXsigmoidY!(X::CudaArray{Float32},Y::CudaArray{Float32},T::CudaArray{Float32},Out::CudaArray{Float32})
    tmp=CudaArray(Float32,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(DXbinaryentropyXsigmoidY_kernel32,nblocks,1024,(length(X),X,Y,T,Out))
end
export DXbinaryentropyXsigmoidY!

md=CuModule(DIR*"DYbinaryentropyXsigmoidY.ptx",false)
DYbinaryentropyXsigmoidY_kernel=CuFunction(md,"DYbinaryentropyXsigmoidY")
function DYbinaryentropyXsigmoidY!(X::CudaArray,Y::CudaArray,T::CudaArray,Out::CudaArray)
    tmp=CudaArray(Float64,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(DYbinaryentropyXsigmoidY_kernel,nblocks,1024,(length(X),X,Y,T,Out))
end

md=CuModule(DIR*"DYbinaryentropyXsigmoidY_32.ptx",false)
DYbinaryentropyXsigmoidY_kernel32=CuFunction(md,"DYbinaryentropyXsigmoidY_32")
function DYbinaryentropyXsigmoidY!(X::CudaArray{Float32},Y::CudaArray{Float32},T::CudaArray{Float32},Out::CudaArray{Float32})
    tmp=CudaArray(Float32,size(X))
    nblocks=round(Int,ceil(length(X)/1024))
    launch(DYbinaryentropyXsigmoidY_kernel32,nblocks,1024,(length(X),X,Y,T,Out))
end
export DYbinaryentropyXsigmoidY!

md=CuModule(DIR*"sigmoid.ptx",false)
sigmoid_kernel=CuFunction(md,"sigmoid")
function sigmoid(A::CudaArray{Float64})
    out=CudaArray(Float64,size(A))
    nblocks=round(Int,ceil(length(A)/1024))
    launch(sigmoid_kernel,nblocks,1024,(length(A),A,out))
    return out
end

md=CuModule(DIR*"sigmoid32.ptx",false)
sigmoid_kernel_32=CuFunction(md,"sigmoid32")
function sigmoid(A::CudaArray{Float32})
    out=CudaArray(Float32,size(A))
    nblocks=round(Int,ceil(length(A)/1024))
    launch(sigmoid_kernel_32,nblocks,1024,(length(A),A,out))
    return out
end
export sigmoid

function sigmoid!(A::CudaArray{Float64},out::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(sigmoid_kernel,nblocks,1024,(length(A),A,out))
end

function sigmoid!(A::CudaArray{Float32},out::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(sigmoid_kernel_32,nblocks,1024,(length(A),A,out))
end
export sigmoid!

md=CuModule(DIR*"vsquare.ptx",false)
vsquarekernel=CuFunction(md,"vsquare")
function vsquare(A::CudaArray{Float64})
    out=zeros(Float64,A)
    launch(vsquarekernel,size(A,1),size(A,2),(A,out))
    return out
end

md=CuModule(DIR*"vsquare_32.ptx",false)
vsquarekernel_32=CuFunction(md,"vsquare_32")
function vsquare(A::CudaArray{Float32})
    out=zeros(Float32,A)
    launch(vsquarekernel_32,size(A,1),size(A,2),(A,out))
    return out
end

#function vsquareout!(A::CudaArray{Float64},Out::CudaArray{Float64})
#    launch(vsquarekernel,size(A,1),size(A,2),(A,Out))
#end

#function vsquareout!(A::CudaArray{Float32},Out::CudaArray{Float32})
#    launch(vsquarekernel_32,size(A,1),size(A,2),(A,Out))
#end

import Base.exp
md=CuModule(DIR*"exp.ptx",false)
kernel_exp=CuFunction(md,"expkernel")
function exp!(A::CudaArray{Float64},B::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_exp,nblocks,1024,(length(A),A,B))
end

md=CuModule(DIR*"exp_32.ptx",false)
kernel_exp_32=CuFunction(md,"expkernel_32")
function exp!(A::CudaArray{Float32},B::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_exp_32,nblocks,1024,(length(A),A,B))
end

function exp(A::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    B=CudaArray(Float64,size(A))
    launch(kernel_exp,nblocks,1024,(length(A),A,B))
    return B # MEMORY LEAK -- HOW TO REMOVE TEMPORARY B?
end

function exp(A::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    B=CudaArray(Float32,size(A))
    launch(kernel_exp_32,nblocks,1024,(length(A),A,B))
    return B # MEMORY LEAK -- HOW TO REMOVE TEMPORARY B?
end

export exp,exp!

import Base.log
md=CuModule(DIR*"log.ptx",false)
kernel_log=CuFunction(md,"logkernel")
function log!(A::CudaArray{Float64},B::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_log,nblocks,1024,(length(A),A,B))
end

md=CuModule(DIR*"log_32.ptx",false)
kernel_log_32=CuFunction(md,"logkernel_32")
function log!(A::CudaArray{Float32},B::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_log_32,nblocks,1024,(length(A),A,B))
end

function log(A::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    B=CudaArray(Float64,size(A))
    launch(kernel_log,nblocks,1024,(length(A),A,B))
    return B # MEMORY LEAK -- HOW TO REMOVE TEMPORARY B?
end

function log(A::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    B=CudaArray(Float32,size(A))
    launch(kernel_log_32,nblocks,1024,(length(A),A,B))
    return B # MEMORY LEAK -- HOW TO REMOVE TEMPORARY B?
end

export log,log!

md=CuModule(DIR*"rectlin.ptx",false)
kernel_rectlin=CuFunction(md,"rectlin")
function rectlin!(A::CudaArray{Float64},B::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_rectlin,nblocks,1024,(length(A),A,B))
end

md=CuModule(DIR*"rectlin_32.ptx",false)
kernel_rectlin_32=CuFunction(md,"rectlin_32")
function rectlin!(A::CudaArray{Float32},B::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_rectlin_32,nblocks,1024,(length(A),A,B))
end
export rectlin!

md=CuModule(DIR*"alex.ptx",false)
kernel_alex=CuFunction(md,"alex")
function alex!(A::CudaArray{Float64},B::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_alex,nblocks,1024,(length(A),A,B))
end

md=CuModule(DIR*"alex_32.ptx",false)
kernel_alex_32=CuFunction(md,"alex_32")
function alex!(A::CudaArray{Float32},B::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_alex_32,nblocks,1024,(length(A),A,B))
end
export alex!

md=CuModule(DIR*"gradalex.ptx",false)
kernel_gradalex=CuFunction(md,"gradalex")
function gradalex!(A::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_alex,nblocks,1024,(length(A),A,B,C))
end

md=CuModule(DIR*"gradalex_32.ptx",false)
kernel_gradalex_32=CuFunction(md,"gradalex_32")
function gradalex!(A::CudaArray{Float32},B::CudaArray{Float32},C::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_gradalex_32,nblocks,1024,(length(A),A,B,C))
end
export gradalex!

md=CuModule(DIR*"kinklin.ptx",false)
kernel_kinklin=CuFunction(md,"kinklin")
function kinklin!(A::CudaArray{Float64},B::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_kinklin,nblocks,1024,(length(A),0.25,A,B))
end

md=CuModule(DIR*"kinklin_32.ptx",false)
kernel_kinklin_32=CuFunction(md,"kinklin_32")
function kinklin!(A::CudaArray{Float32},B::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_kinklin_32,nblocks,1024,(length(A),0.25,A,B))
end
export kinklin!

function rectlin(A::CudaArray)
    out=CudaArray(eltype(A),size(A))
    rectlin!(A,out)
    return out
end
export rectlin

function alex(A::CudaArray)
    out=CudaArray(eltype(A),size(A))
    alex!(A,out)
    return out
end
export alex

# A.*(B.>0):
md=CuModule(DIR*"A_emult_Bg0.ptx",false)
kernel_A_emult_Bg0=CuFunction(md,"A_emult_Bg0")
function A_emult_Bg0!(A::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_A_emult_Bg0,nblocks,1024,(length(A),A,B,C))
end

md=CuModule(DIR*"A_emult_Bg0_32.ptx",false)
kernel_A_emult_Bg0_32=CuFunction(md,"A_emult_Bg0_32")
function A_emult_Bg0!(A::CudaArray{Float32},B::CudaArray{Float32},C::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    launch(kernel_A_emult_Bg0_32,nblocks,1024,(length(A),A,B,C))
end

export A_emult_Bg0!

md=CuModule(DIR*"alphaaxpy.ptx",false)
alphaaxpy=CuFunction(md,"alphaaxpy")

function alphaaxpy!(alpha::Float64,A::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(alphaaxpy,nblocks,1024,(length(C),alpha,A,B,C))
end

md=CuModule(DIR*"alphaaxpy_32.ptx",false)
alphaaxpy_32=CuFunction(md,"alphaaxpy_32")

function alphaaxpy!(alpha::Real,A::CudaArray{Float32},B::CudaArray{Float32},C::CudaArray{Float32})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(alphaaxpy_32,nblocks,1024,(length(C),Float32(alpha),A,B,C))
end
export alphaaxpy!

md=CuModule(DIR*"alphaax.ptx",false)
alphaax=CuFunction(md,"alphaax")
function alphaax!(alpha::Float64,A::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(alphaax,nblocks,1024,(length(C),alpha,A,B,C))
end

md=CuModule(DIR*"alphaax_32.ptx",false)
alphaax_32=CuFunction(md,"alphaax_32")
function alphaax!(alpha::Float32,A::CudaArray{Float32},B::CudaArray{Float32},C::CudaArray{Float32})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(alphaax_32,nblocks,1024,(length(C),alpha,A,B,C))
end
export alphaax!

md=CuModule(DIR*"gaxpy.ptx",false)
gaxpy=CuFunction(md,"gaxpy")
function gaxpy!(alpha::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(gaxpy,nblocks,1024,(length(C),alpha,B,C))
end

md=CuModule(DIR*"gaxpy_32.ptx",false)
gaxpy_32=CuFunction(md,"gaxpy_32")
function gaxpy!(alpha::CudaArray{Float32},B::CudaArray{Float32},C::CudaArray{Float32})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(gaxpy_32,nblocks,1024,(length(C),alpha,B,C))
end
export gaxpy!

axpy!(alpha::CudaArray,B::CudaArray,C::CudaArray)=gaxpy!(alpha::CudaArray,B::CudaArray,C::CudaArray)

export axpy!

# NB: I've defined this differently from the Base.scale! -- here scale!(a,X) scales X by factor a and overwrites X. This is unfortunate, but I defined this according to the convention that overwrite functions overwrite their last argument.
import Base.scale!
md=CuModule(DIR*"gscale.ptx",false)
kernel_gscale=CuFunction(md,"gscale")
function scale!(alpha::CudaArray{Float64},B::CudaArray{Float64})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(kernel_gscale,nblocks,1024,(length(B),alpha,B))
end

md=CuModule(DIR*"gscale_32.ptx",false)
kernel_gscale_32=CuFunction(md,"gscale_32")
function scale!(alpha::CudaArray{Float32},B::CudaArray{Float32})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(kernel_gscale_32,nblocks,1024,(length(B),alpha,B))
end
export scale!

md=CuModule(DIR*"diag.ptx",false)
diag_kernel=CuFunction(md,"diag_kernel")
function diag!(A::CudaArray{Float64},B::CudaArray{Float64})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(diag_kernel,nblocks,1024,(length(B),A,B))
end

md=CuModule(DIR*"diag_32.ptx",false)
diag_kernel_32=CuFunction(md,"diag_kernel_32")
function diag!(A::CudaArray{Float32},B::CudaArray{Float32})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(diag_kernel_32,nblocks,1024,(length(B),A,B))
end
export diag!

md=CuModule(DIR*"diagm.ptx",false)
diagm_kernel=CuFunction(md,"diagm_kernel")
function diagm!(A::CudaArray{Float64},B::CudaArray{Float64})
    nblocks=round(Int,ceil(length(A)/1024))
    fill!(B,0.0)
    launch(diagm_kernel,nblocks,1024,(length(A),A,B))
end

md=CuModule(DIR*"diagm_32.ptx",false)
diagm_kernel_32=CuFunction(md,"diagm_kernel_32")
function diagm!(A::CudaArray{Float32},B::CudaArray{Float32})
    nblocks=round(Int,ceil(length(A)/1024))
    fill!(B,Float32(0.0))
    launch(diagm_kernel_32,nblocks,1024,(length(A),A,B))
end
export diagm!

md=CuModule(DIR*"gax.ptx",false)
gax=CuFunction(md,"gax")
function gax!(alpha::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(gax,nblocks,1024,(length(C),alpha,B,C))
end

md=CuModule(DIR*"gax_32.ptx",false)
gax_32=CuFunction(md,"gax_32")
function gax!(alpha::CudaArray{Float32},B::CudaArray{Float32},C::CudaArray{Float32})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(gax_32,nblocks,1024,(length(C),alpha,B,C))
end
export gax!

md=CuModule(DIR*"ax.ptx",false)
ax=CuFunction(md,"ax")
function ax!(alpha::Float64,B::CudaArray{Float64},C::CudaArray{Float64})
       nblocks=round(Int,ceil(length(B)/1024))
    launch(ax,nblocks,1024,(length(C),alpha,B,C))
end

md=CuModule(DIR*"ax_32.ptx",false)
ax_32=CuFunction(md,"ax_32")
function ax!(alpha::Float32,B::CudaArray{Float32},C::CudaArray{Float32})
       nblocks=round(Int,ceil(length(B)/1024))
    launch(ax_32,nblocks,1024,(length(C),alpha,B,C))
end
export ax!

md=CuModule(DIR*"vmultbang.ptx",false)
vmultbangkernel=CuFunction(md,"vmultbang")
function vmult!(alpha::Float64,A::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(vmultbangkernel,nblocks,1024,(length(A),alpha,A,B,C))
end

md=CuModule(DIR*"vmultbang_32.ptx",false)
vmultbangkernel_32=CuFunction(md,"vmultbang_32")
function vmult!(alpha::Real,A::CudaArray{Float32},B::CudaArray{Float32},C::CudaArray{Float32})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(vmultbangkernel_32,nblocks,1024,(length(A),Float32(alpha),A,B,C))
end
export vmult!

md=CuModule(DIR*"vdivbang.ptx",false)
vdivbangkernel=CuFunction(md,"vdivbang")
function vdiv!(alpha::Float64,A::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(vdivbangkernel,nblocks,1024,(length(A),alpha,A,B,C))
end

md=CuModule(DIR*"vdivbang_32.ptx",false)
vdivbangkernel_32=CuFunction(md,"vdivbang_32")
function vdiv!(alpha::Real,A::CudaArray{Float32},B::CudaArray{Float32},C::CudaArray{Float32})
    nblocks=round(Int,ceil(length(B)/1024))
    launch(vdivbangkernel,nblocks,1024,(length(A),Float32(alpha),A,B,C))
end
export vdiv!

md=CuModule(DIR*"vmultbangupdate.ptx",false)
vmultbangupdate_kernel=CuFunction(md,"vmultbangupdate")
function vmultupdate!(alpha::Float64,A::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
        nblocks=round(Int,ceil(length(B)/1024))
    launch(vmultbangupdate_kernel,nblocks,1024,(length(A),alpha,A,B,C))
end

md=CuModule(DIR*"vmultbangupdate_32.ptx",false)
vmultbangupdate_kernel_32=CuFunction(md,"vmultbangupdate_32")
function vmultupdate!(alpha::Real,A::CudaArray{Float32},B::CudaArray{Float32},C::CudaArray{Float32})
        nblocks=round(Int,ceil(length(B)/1024))
    launch(vmultbangupdate_kernel_32,nblocks,1024,(length(A),alpha,A,B,C))
end

#md=CuModule(DIR*"vdivbangupdate.ptx",false)
#vdivbangupdate_kernel=CuFunction(md,"vdivbangupdate")
#function vdivupdate!(alpha::Float64,A::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64})
#    nblocks=round(Int,ceil(length(B)/1024))
#    launch(vdivbangupdate_kernel,nblocks,1024,(length(A),alpha,A,B,C))
#end

#md=CuModule(DIR*"vAoverBupdate.ptx",false)
#vAoverBupdate_kernel=CuFunction(md,"vAoverBupdate")
#function vAoverBupdate!(alpha::Float64,A::CudaArray{Float64},B::CudaArray{Float64},C::CudaArray{Float64},D::CudaArray{Float64})
#    nblocks=round(Int,ceil(length(B)/1024))
#    launch(vAoverBupdate_kernel,nblocks,1024,(length(A),alpha,A,B,C,D))
#end

function axpy(a::Float64,x::CudaArray{Float64},y::CudaArray{Float64})
    out=copy(y)
    CUBLAS.axpy!(a,x,out)
    return out
end

function axpy(a::Real,x::CudaArray{Float32},y::CudaArray{Float32})
    out=copy(y)
    CUBLAS.axpy!(Float32(a),x,out)
    return out
end
export axpy

flatten(T::Type,g::CudaArray)=reinterpret(T,g,(length(g),1)) # make a (prod(dims),1) CudaArray from a CudaArray
# This is useful since we can then call CUBLAS routines that operate on vectors

flatten(g::CudaArray)=reinterpret(eltype(g),g,(length(g),1)) # make a (prod(dims),1) CudaArray from a CudaArray
# This is useful since we can then call CUBLAS routines that operate on vectors
export flatten

export flatten

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using GNSSTools
using Statistics
using Random
using StatsBase
using LinearAlgebra
using Base.Threads
using JLD
using PyPlot
pygui(true)
include("ch4_data_plot_functions.jl")
include("ch4_constellation_hist_structs.jl")
include("ch4_simulation_function.jl")

directory = "/home/sjbilardi/Dropbox/Apps/Overleaf/MastersThesis/"
folder = "figures/l1ca_sim/"

# doppler_file_name = string(directory, "ch1_doppler_curves.jld")
# doppler_raw_gps, doppler_rate_raw_gps, sig_freq = load(doppler_file_name,
#                                                        "doppler_raw_gps",
#                                                        "doppler_rate_raw_gps",
#                                                        "freq")
# gps_hist = initiate_hist(doppler_raw_gps, doppler_rate_raw_gps; nbins=1000)
# iridium_hist = initiate_hist(doppler_raw_iridium, doppler_rate_raw_iridium;
#                              nbins=1000)

# t_length = 2
# f_s = 25e6
# f_if = 0
# Tsys = 535
# nADC = 8
# N = floor(Int, f_s*1e-3)
# B = 25e6
# prn_num = 9
# upsample_factor = 4
# prns = []
# for i in 1:prn_num
#     prn = rand(1:32)
#     while prn in prns
#         prn = rand(1:32)
#     end
#     push!(prns, prn)
# end
# prns = sort(prns)
# phis = rand(0:0.0001:2π, prn_num)
# cn0s = rand(35:0.01:45, prn_num)
# f_ds = Array{Float64}(undef, prn_num)
# fd_rates = Array{Float64}(undef, prn_num)
# for i in 1:prn_num
#     f_d, fd_rate = sample_distribution(gps_hist)
#     f_ds[i] = f_d
#     fd_rates[i] = fd_rate
# end
# code_start_idxs = rand(1:N, prn_num) + Rational.(1, (rand(1:upsample_factor, prn_num)))
# h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15] ./ 16.3676e6^2  # Rakon IT5300B TCXO

# save(string(directory, folder, "ch4_l1ca_sim_parms.jld"),
#             "t_length", t_length,
#             "f_s", f_s,
#             "f_if", f_if,
#             "Tsys", Tsys,
#             "nADC", nADC,
#             "N", N,
#             "B", B,
#             "prn_num", prn_num,
#             "upsample_factor", upsample_factor,
#             "prns", prns,
#             "phis", phis,
#             "cn0s", cn0s,
#             "f_ds", f_ds,
#             "fd_rates", fd_rates,
#             "code_start_idxs", code_start_idxs,
#             "h_parms", h_parms)

parms = load(string(directory, folder, "ch4_l1ca_sim_parms.jld"),
             "t_length",
             "f_s",
             "f_if",
             "Tsys",
             "nADC",
             "N",
             "B",
             "prn_num",
             "upsample_factor",
             "prns",
             "phis",
             "cn0s",
             "f_ds",
             "fd_rates",
             "code_start_idxs",
             "h_parms")

t_length, f_s, f_if, Tsys, nADC, N, B, prn_num, upsample_factor, prns,
phis, cn0s, f_ds, fd_rates, code_start_idxs, h_parms = parms

channel = "I"
sigtype = define_l1ca_code_type(t_length; B=B)

data = simulate_signal(sigtype, f_s, t_length, prns, cn0s, f_ds, fd_rates,
                       phis, code_start_idxs, h_parms;
                       nADC=nADC, f_if=f_if, Tsys=Tsys)

replica = definereplica(sigtype, 25e6, 4e-3)
detectable_prns = []
for i in 1:32
    fd_est, n0_est, SNR_est,
    P_d, above_threshold = courseacquisition(data, replica, i; M=20,
                                             Δfd=1/(1*replica.t_length))
    if above_threshold && (P_d >= 0.9)
        detectable = true
        push!(detectable_prns, i)
    else
        detectable = false
    end
    f_code_d, f_code_dd = GNSSTools.calc_doppler_code_rate(l1ca_chipping_rate,
                                                           replica.signal_type.sig_freq,
                                                           fd_est, 0)
    n0 = calcinitcodephase(l1ca_code_length, f_code_d, f_code_dd, data.f_s,
                           n0_est)
    println("PRN $i\t$SNR_est\t$fd_est\t$n0_est\t$n0\t$(Float64(round(P_d, digits=2)))\t$detectable")
end
println("$(length(detectable_prns)) Detected PRNs: $detectable_prns")

fig = figure(figsize=(7.5, 5))
acquisition_T = 4e-3
fine_acq_T = 10e-3
fine_acq_sub_T = 1e-3
tracking_T = 1e-3
M = 20
linewidth = 0.5
tracking_σᵩ = Array{Float64}(undef, length(detectable_prns))
tracking_snr = Array{Float64}(undef, length(detectable_prns))
tracking_code_phase_err = Array{Float64}(undef, length(detectable_prns))
for i in 1:length(detectable_prns)
    prn = detectable_prns[i]
    ax = fig.add_subplot(3, 3, i)
    tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
    title("PRN $prn")
    acqresults, trackresults = process(data, define_l1ca_code_type(), prn;
                                       show_plot=false,
                                       fine_acq_method=:fft,
                                       state_num=3,
                                       dynamickf=true,
                                    #    h₀=h_parms[3], h₋₂=h_parms[1],
                                       dll_b=0.01,
                                       acquisition_T=acquisition_T,
                                       fine_acq_T=fine_acq_T,
                                       fine_acq_sub_T=fine_acq_sub_T,
                                       tracking_T=tracking_T,
                                       use_fine_acq=true,
                                       M=M)
    ax.plot(trackresults.t, real.(trackresults.ZP), linewidth=linewidth)
    ax.plot(trackresults.t, imag.(trackresults.ZP), linewidth=linewidth)
    N = floor(Int, length(trackresults.t)/2)
    tracking_σᵩ[i] = std(trackresults.dphi_meas[N:end]) * 180 / π
    tracking_snr[i] = mean(trackresults.SNR[N:end])
    tracking_code_phase_err[i] = std(trackresults.code_err_meas[N:end])
end
subplots_adjust(hspace=0.4, wspace=0.4,  top=0.88, left=0.02, right=0.98,
                bottom=0.02)
plot_title = string("Course Acq. T=$(floor(Int, acquisition_T*1000))ms    ",
                    "Fine Acq. T=$(floor(Int, fine_acq_T*1000))ms    ",
                    "Tracking T=$(floor(Int, tracking_T*1000))ms    ",
                    "M=$(M)")
suptitle(plot_title)
savefig(string(directory, folder, "ch4_l1ca_acquired_prns.svg"), dpi=300)

for i in 1:length(detectable_prns)
    detected_prn = detectable_prns[i]
    j = argmax(detected_prn .== prns)
    acqresults, trackresults, corr_result, SNR_est,
    P_d, above_threshold = process(data,
                                sigtype,
                                detected_prn, channel;
                                fine_acq_method=:fft,
                                return_corrresult=true,
                                return_Pd=true,
                                state_num=3,
                                #    h₀=h_parms[3],
                                #    h₋₂=h_parms[1],
                                q_a=1,
                                dll_b=1,
                                acquisition_T=4e-3,
                                fine_acq_T=10e-3,
                                fine_acq_sub_T=1e-3,
                                tracking_T=1e-3,
                                show_plot=false,
                                σω=1000,
                                M=20,
                                use_fine_acq=true);
    # Calculate truth info
    dt = 0.001
    code_start_idx = code_start_idxs[j]
    f_d = f_ds[j]
    fd_rate = fd_rates[j]
    phi = phis[j]
    cn0 = cn0s[j]
    t, doppler, doppler_rate, phi, code_phase, snr = get_signal_truth_parms(t_length,
                                                                        f_s,
                                                                        code_start_idx,
                                                                        f_d, fd_rate,
                                                                        phi, cn0, Tsys,
                                                                        sigtype, channel,
                                                                        1e-3; dt=dt, B=B)
    name = "pdfs/ch4_l1ca_prn_$(trackresults.prn).pdf"
    file_name = string(directory, folder, name)
    plot_example_track_results(trackresults; f_if=0, code_length=l1ca_code_length,
                            save_to_file=file_name,
                            bottom=0.075, top=0.95, left=0.125, right=0.875,
                            hspace=0.75, wspace=0.35, marker_size=1,
                            line_width=0.8,
                            truth_t=t, truth_doppler=doppler,
                            truth_doppler_rate=doppler_rate,
                            truth_code_phase=code_phase,
                            truth_snr=snr, truth_phi=phi)
end


text:
using MeshArrays, NCDatasets, NCTiles
grid = GridSpec("LLC90","grids/")

exampledir = joinpath("data","ex4")
datadir = joinpath(exampledir,"model_output")
saveloc = joinpath(exampledir,"tiles")
fnames = joinpath.(Ref(datadir),filter(x->occursin("data",x),readdir(datadir)))

metafile = parsemeta(joinpath(datadir,"ptr_3d_set1.0000000732.meta"))
README = readlines(joinpath(datadir,"README"))

timeunits = "days since 1992-1-1 0:0:0"
nsteps = 12
timeinterval = 30
time_steps = timeinterval/2:timeinterval:timeinterval*nsteps
dep = -1 .* collect([-5 -15 -25 -35 -45 -55 -65 -75.0049972534180 -85.0250015258789 -95.0950012207031 -105.309997558594 -115.870002746582 -127.150001525879 -139.740005493164 -154.470001220703 -172.399993896484 -194.735000610352 -222.710006713867 -257.470001220703 -299.929992675781 -350.679992675781 -409.929992675781 -477.470001220703 -552.710021972656 -634.734985351563 -722.400024414063 -814.469970703125 -909.739990234375 -1007.15502929688 -1105.90502929688 -1205.53503417969 -1306.20495605469 -1409.15002441406  -1517.09497070313 -1634.17504882813 -1765.13500976563 -1914.15002441406 -2084.03491210938 -2276.22509765625 -2491.25000000000 -2729.25000000000 -2990.25000000000 -3274.25000000000 -3581.25000000000 -3911.25000000000 -4264.25000000000 -4640.25000000000 -5039.25000000000 -5461.25000000000 -5906.25000000000]')
prec = Float32
iosize = (90,1170,50)
tilesize = (90,90)
dims = [
    NCvar("i_c","1",tilesize[1],1:tilesize[1],Dict("long_name" => "Cartesian coordinate 1"),NCDatasets),
    NCvar("j_c","1",tilesize[2],1:tilesize[2],Dict("long_name" => "Cartesian coordinate 2"),NCDatasets),
    NCvar("dep_c","m",size(dep),dep,Dict("long_name" => "depth","standard_name" => "depth","positive" => "down"),NCDatasets),
    NCvar("tim",timeunits,Inf,time_steps,Dict(("long_name" => "time","standard_name" => "time")),NCDatasets)
]

gridvars = GridLoad(grid)

land = gridvars["hFacC"]
for f in land.fIndex
    for d in 1:size(land,2)
        land[f,d][land[f,d].==0] .= NaN
        land[f,d][land[f,d].>0] .= 1
    end
end

tilarea = TileData(gridvars["RAC"],tilesize,grid)
tilland = TileData(land,tilesize,grid)
thic = gridvars["DRC"]

fldidx = 1:106
for fidx in fldidx
    @time begin
        fldname = metafile["fldList"][fidx]
        println("Processing "*fldname)
        diaginfo = readAvailDiagnosticsLog(joinpath(exampledir,"available_diagnostics.log"),fldname)
        fielddata = BinData(fnames,prec,iosize,fidx)
        tilfld = TileData(fielddata,tilesize,grid)
        tillat = TileData(gridvars["YC"],tilfld.tileinfo,tilfld.tilesize,tilfld.precision,tilfld.numtiles)
        tillon = TileData(gridvars["XC"],tilfld.tileinfo,tilfld.tilesize,tilfld.precision,tilfld.numtiles)
        flds = Dict([fldname => NCvar(fldname,diaginfo["units"],dims,tilfld,Dict(),NCDatasets),
                    "lon" => NCvar("lon","degrees_east",dims[1:2],tillon,Dict("long_name" => "longitude"),NCDatasets),
                    "lat" => NCvar("lat","degrees_north",dims[1:2],tillat,Dict("long_name" => "latitude"),NCDatasets),
                    "area" => NCvar("area","m^2",dims[1:2],tilarea,Dict(["long_name" => "grid cell area", "standard_name" => "cell_area"]),NCDatasets),
                    "land" => NCvar("land","1",dims[1:3],tilland,Dict(["long_name" => "land mask", "standard_name" => "land_binary_mask"]),NCDatasets),
                    "thic" => NCvar("thic","m",dims[3],thic,Dict("standard_name" => "cell_thickness"),NCDatasets)
        ])
        savepath = joinpath(saveloc,fldname)
        if !isdir(savepath); mkpath(savepath); end

        numtiles = flds[fldname].values.numtiles
        savenames = joinpath.(Ref(savepath),fldname*".".*lpad.(string.(1:numtiles),4,"0").*".nc")
        rm.(savenames, force=true)

        datasets = [createfile(savenames[tidx],flds,README, itile = tidx, ntile = length(savenames)) for tidx in 1:length(savenames)]

        ds = [x[1] for x in datasets]
        fldvars = [x[2] for x in datasets]
        #dims = [x[3] for x in datasets]

        for k in keys(flds)
            if isa(flds[k].values,TileData)
                addData(fldvars,flds[k])
            else
                tmpfldvars = [fv[findfirst(isequal(k),name.(fv))] for fv in fldvars]
                addData.(tmpfldvars,Ref(flds[k]))
            end
        end

        for dim in dims
            addDimData.(ds,Ref(dim))
        end

        close.(ds)
    end
end