text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Dierckx

using Compat

export Spline1D,
       Spline2D,
       evaluate,
       derivative,
       integrate,
       roots,
       evalgrid,
       get_knots,
       get_coeffs,
       get_residual

import Base: show

unixpath = "../deps/src/ddierckx/libddierckx"
winpath = "../deps/bin$(Sys.WORD_SIZE)/libddierckx"
const ddierckx = joinpath(dirname(@__FILE__), is_unix() ? unixpath : winpath)

function __init__()
    # Ensure library is available.
    if (Libdl.dlopen_e(ddierckx) == C_NULL)
        error("Dierckx not properly installed. Run Pkg.build(\"Dierckx\")")
    end
end

# ----------------------------------------------------------------------------
# 1-d splines

const _fit1d_messages = @compat Dict(
2=>
"""A theoretically impossible result was found during the iteration
process for finding a smoothing spline with fp = s: s too small.
There is an approximation returned but the corresponding weighted sum
of squared residuals does not satisfy the condition abs(fp-s)/s <
tol.""",
3=>
"""The maximal number of iterations maxit (set to 20 by the program)
allowed for finding a smoothing spline with fp=s has been reached: s
too small. There is an approximation returned but the corresponding
weighted sum of squared residuals does not satisfy the condition
abs(fp-s)/s < tol.""",
10=>
"""Error on entry, no approximation returned. The following conditions
must hold:
1<=k<=5
x[1] < x[2] < ... < x[end]
w[i] > 0.0 for all i

Additionally, if spline knots are given:
length(xknots) <= length(x) + k + 1
x[1] < xknots[1] < xknots[k+2] < ... < xknots[end] < x[end]
The schoenberg-whitney conditions: there must be a subset of data points
xx[j] such that t[j] < xx[j] < t[j+k+1] for j=1,2,...,n-k-1""")

const _eval1d_messages = @compat Dict(
1=>
"""Input point out of range""",
10=>
"""Invalid input data. The following conditions must hold:
length(x) != 0 and xb <= x[1] <= x[2] <= ... x[end] <= xe""")

_translate_bc(bc::AbstractString) = (bc == "extrapolate" ? 0 :
                             bc == "zero" ? 1 :
                             bc == "error" ? 2 :
                             bc == "nearest" ? 3 :
                             error("unknown boundary condition: \"$bc\""))
_translate_bc(bc::Int) = (bc == 0 ? "extrapolate" :
                          bc == 1 ? "zero" :
                          bc == 2 ? "error" :
                          bc == 3 ? "nearest" : "")

type Spline1D
    t::Vector{Float64}
    c::Vector{Float64}
    k::Int
    bc::Int
    fp::Float64
end

get_knots(spl::Spline1D) = spl.t[spl.k+1:end-spl.k]
get_coeffs(spl::Spline1D) = spl.c[1:end-spl.k+1]
get_residual(spl::Spline1D) = spl.fp

function reallycompact(a::Vector)
    io = IOBuffer()
    if length(a) <= 5
        show(io, a)
    else
        write(io, "[")
        showcompact(io, a[1])
        write(io, ",")
        showcompact(io, a[2])
        write(io, " \u2026 ")
        showcompact(io, a[end-1])
        write(io, ",")
        showcompact(io, a[end])
        write(io, "]")
        write(io, " ($(length(a)) elements)")
    end
    seekstart(io)
    readstring(io)
end

function show(io::IO, spl::Spline1D)
    print(io, """Spline1D(knots=$(reallycompact(get_knots(spl))), k=$(spl.k), extrapolation=\"$(_translate_bc(spl.bc))\", residual=$(spl.fp))""")
end

function Spline1D(x::AbstractVector, y::AbstractVector;
                  w::AbstractVector=ones(length(x)),
                  k::Int=3, s::Real=0.0, bc::AbstractString="nearest")
    m = length(x)
    length(y) == m || error("length of x and y must match")
    length(w) == m || error("length of x and w must match")
    m > k || error("length(x) must be greater than k")

    # ensure inputs are of correct type
    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)
    win = convert(Vector{Float64}, w)

    # outputs
    nest = m+k+1
    n = Array(Int32, 1)
    t = Array(Float64, nest)
    c = Array(Float64, nest)
    fp = Array(Float64, 1)
    ier = Array(Int32, 1)

    # working space
    lwrk = m * (k+1) + nest * (7 + 3*k)
    wrk = Array(Float64, lwrk)
    iwrk = Array(Int32, nest)

    ccall((:curfit_, ddierckx), Void,
          (Ptr{Int32}, Ptr{Int32},  # iopt, m
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # x, y, w
           Ptr{Float64}, Ptr{Float64},  # xb, xe
           Ptr{Int32}, Ptr{Float64},  # k, s
           Ptr{Int32}, Ptr{Int32}, # nest, n
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # t, c, fp
           Ptr{Float64}, Ptr{Int32}, Ptr{Int32},  # wrk, lwrk, iwrk
           Ptr{Int32}),  # ier
          &0, &m, xin, yin, win, &xin[1], &xin[end], &k, &@compat(Float64(s)),
          &nest, n, t, c, fp, wrk, &lwrk, iwrk, ier)

    ier[1] <= 0 || error(_fit1d_messages[ier[1]])

    # resize output arrays
    resize!(t, n[1])
    resize!(c, n[1] - k - 1)

    return Spline1D(t, c, k, _translate_bc(bc), fp[1])
end

# version with user-supplied knots
function Spline1D(x::AbstractVector, y::AbstractVector,
                  xknots::AbstractVector;
                  w::AbstractVector=ones(length(x)),
                  k::Int=3, bc::AbstractString="nearest")
    m = length(x)
    length(y) == m || error("length of x and y must match")
    length(w) == m || error("length of x and w must match")
    m > k || error("k must be less than length(x)")

    # ensure inputs are of correct type
    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)
    win = convert(Vector{Float64}, w)

    # x knots
    # (k+1) knots will be added on either end of interior knots.
    n = length(xknots) + 2 * (k+1)
    t = Array(Float64, n)  # All knots
    for i in 1:length(xknots)
        t[i+k+1] = xknots[i]
    end

    # outputs
    c = Array(Float64, n)
    fp = Array(Float64, 1)
    ier = Array(Int32, 1)

    # working space
    lwrk = m * (k+1) + n * (7 + 3*k)
    wrk = Array(Float64, lwrk)
    iwrk = Array(Int32, n)

    ccall((:curfit_, ddierckx), Void,
          (Ptr{Int32}, Ptr{Int32},  # iopt, m
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # x, y, w
           Ptr{Float64}, Ptr{Float64},  # xb, xe
           Ptr{Int32}, Ptr{Float64},  # k, s
           Ptr{Int32}, Ptr{Int32}, # nest, n
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # t, c, fp
           Ptr{Float64}, Ptr{Int32}, Ptr{Int32},  # wrk, lwrk, iwrk
           Ptr{Int32}),  # ier
          &(-1), &m, xin, yin, win, &xin[1], &xin[end], &k, &(-1.0), &n,
          &n, t, c, fp, wrk, &lwrk, iwrk, ier)

    ier[1] <= 0 || error(_fit1d_messages[ier[1]])
    resize!(c, n-k-1)

    return Spline1D(t, c, k, _translate_bc(bc), fp[1])
end

function evaluate(spline::Spline1D, x::AbstractVector)
    m = length(x)
    xin = convert(Vector{Float64}, x)
    y = Array(Float64, m)
    ier = Array(Int32, 1)
    ccall((:splev_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # t, n
           Ptr{Float64}, Ptr{Int32},  # c, k
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32},  # x, y, m
           Ptr{Int32}, Ptr{Int32}),  # e, ier
          spline.t, &length(spline.t), spline.c, &spline.k,
          xin, y, &m, &spline.bc, ier)
    ier[1] == 0 || error(_eval1d_messages[ier[1]])
    return y
end

function evaluate(spline::Spline1D, x::Real)
    y = Array(Float64, 1)
    ier = Array(Int32, 1)
    ccall((:splev_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # t, n
           Ptr{Float64}, Ptr{Int32},  # c, k
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32},  # x, y, m
           Ptr{Int32}, Ptr{Int32}),  # e, ier
          spline.t, &length(spline.t), spline.c, &spline.k,
          &x, y, &1, &spline.bc, ier)
    ier[1] == 0 || error(_eval1d_messages[ier[1]])
    return y[1]
end

# TODO: should the function name be evalder, derivative, or grad?
#       or should it be integrated with evaluate, above?
#       (problem with that: derivative doesn't accept bc="nearest")
# TODO: should `nu` be `d`?
function derivative(spline::Spline1D, x::Vector{Float64}; nu::Integer=1)

    # Unlike evaluate, this Fortran function doesn't accept bc="nearest".
    # However, "nearest" gives a derivative of 0 for all nu outside the
    # support range, so the answer will be the same as bc="zero".
    bc = spline.bc == 3 ? 1 : spline.bc

    # Check derivate requested.
    if !(1 <= nu <= spline.k)
        error("nu must be positive and less than or equal to spline order")
    end

    m = length(x)
    n = length(spline.t)
    wrk = Array(Float64, n)
    y = Array(Float64, m)
    ier = Array(Int32, 1)
    ccall((:splder_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # t, n
           Ptr{Float64}, Ptr{Int32},  # c, k
           Ptr{Int32},  # nu
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32},  # x, y, m
           Ptr{Int32}, Ptr{Float64}, Ptr{Int32}),  # e, wrk, ier
          spline.t, &n, spline.c, &spline.k, &nu,
          x, y, &m, &bc, wrk, ier)
    ier[1] == 0 || error(_eval1d_messages[ier[1]])
    return y
end

derivative(spline::Spline1D, x::Real; nu::Integer=1) =
    derivative(spline, Float64[x]; nu=nu)[1]

function integrate(spline::Spline1D, a::Real, b::Real)
    n = length(spline.t)
    wrk = Array(Float64, n)
    ccall((:splint_, ddierckx), Float64,
          (Ptr{Float64}, Ptr{Int32},  # t, n
           Ptr{Float64}, Ptr{Int32},  # c, k
           Ptr{Float64}, Ptr{Float64},  # a, b
           Ptr{Float64}),  # wrk
          spline.t, &n, spline.c, &spline.k,
          &a, &b, wrk)
end

# note: default maxn in scipy.interpolate is 3 * (length(spline.t) - 7)
function roots(spline::Spline1D; maxn::Integer=8)
    if spline.k != 3
        error("root finding only supported for cubic splines (k=3)")
    end
    n = length(spline.t)
    zeros = Array(Float64, maxn)
    m = Array(Int32, 1)
    ier = Array(Int32, 1)
    ccall((:sproot_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # t, n
           Ptr{Float64}, Ptr{Float64},  # c, zeros
           Ptr{Int32}, Ptr{Int32},  # mest, m
           Ptr{Int32}),  # ier
          spline.t, &n, spline.c, zeros,
          &maxn, m, ier)

    if ier[1] == 0
        return zeros[1:m[1]]
    elseif ier[1] == 1
        warn("number of zeros exceeded maxn; only first maxn zeros returned")
        return zeros
    elseif ier[1] == 10
        error("Invalid input data.")
    else
        error("unknown error code in sproot: $(ier[1])")
    end
end

# ----------------------------------------------------------------------------
# 2-d splines

# NOTE REGARDING ARGUMENT ORDER: In the "grid" version of the Spline2D
# constructor and evaluators, the fortran functions expects z to have
# shape (my, mx), but we'd rather have x be the fast axis in z.  So,
# in the ccall()s in these methods, all the x and y related inputs are
# swapped with regard to what the Fortran documentation says.

const _fit2d_messages = @compat Dict(
-3=>

"""The coefficients of the spline returned have been computed as the
minimal norm least-squares solution of a (numerically) rank deficient
system (deficiency=%i). If deficiency is large, the results may be
inaccurate. Deficiency may strongly depend on the value of eps.""",

1=>

"""The required storage space exceeds the available storage space:
nxest or nyest too small, or s too small. Try increasing s.""",
# The weighted least-squares spline corresponds to the current set of knots.

2=>

"""A theoretically impossible result was found during the iteration
process for finding a smoothing spline with fp = s: s too small or
badly chosen eps.  Weighted sum of squared residuals does not satisfy
abs(fp-s)/s < tol.""",

3=>

"""the maximal number of iterations maxit (set to 20 by the program)
allowed for finding a smoothing spline with fp=s has been reached: s
too small.  Weighted sum of squared residuals does not satisfy
abs(fp-s)/s < tol.""",

4=>

"""No more knots can be added because the number of b-spline
coefficients (nx-kx-1)*(ny-ky-1) already exceeds the number of data
points m: either s or m too small.  The weighted least-squares spline
corresponds to the current set of knots.""",

5=>

"""No more knots can be added because the additional knot would
(quasi) coincide with an old one: s too small or too large a weight to
an inaccurate data point.  The weighted least-squares spline
corresponds to the current set of knots.""",

10=>

"""Error on entry, no approximation returned. The following conditions
must hold:
xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
If iopt==-1, then
  xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
  yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye""")

const _eval2d_message = (
"""Invalid input data. Restrictions:
length(x) != 0, length(y) != 0
x[i-1] <= x[i] for i=2,...,length(x)
y[j-1] <= y[j] for j=2,...,length(y)
""")

type Spline2D
    tx::Vector{Float64}
    ty::Vector{Float64}
    c::Vector{Float64}
    kx::Int
    ky::Int
    fp::Float64
end

get_knots(spl::Spline2D) = (spl.tx[spl.kx+1:end-spl.kx],
                            spl.ty[spl.ky+1:end-spl.ky])
get_residual(spl::Spline2D) = spl.fp

# Helper function for calculating required size of work array in surfit.
# Note that x and y here are as in the Fortran documentation.
function calc_surfit_lwrk1(m, kx, ky, nxest, nyest)
    u = nxest - kx - 1
    v = nyest - ky - 1
    km = max(kx, ky) + 1
    ne = max(nxest, nyest)
    bx = kx*v + ky + 1
    by = ky*u + kx + 1
    b1 = b2 = 0
    if (bx<=by)
        b1 = bx
        b2 = bx + v - ky
    else
        b1 = by
        b2 = by + u - kx
    end
    return u*v*(2 + b1 + b2) + 2*(u+v+km*(m+ne)+ne-kx-ky) + b2 + 1
end

# Construct spline from unstructured data
function Spline2D(x::AbstractVector, y::AbstractVector, z::AbstractVector;
                  w::AbstractVector=ones(length(x)), kx::Int=3, ky::Int=3,
                  s::Real=0.0)

    # array sizes
    m = length(x)
    (length(y) == length(z) == m) || error("lengths of x, y, z must match")
    (length(w) == m) || error("length of w must match other inputs")

    nxest = max(kx+1+ceil(Int,sqrt(m/2)), 2*(kx+1))
    nyest = max(ky+1+ceil(Int,sqrt(m/2)), 2*(ky+1))
    nmax = max(nxest, nyest)

    eps = 1.0e-16

    # bounds
    xb = minimum(x)
    xe = maximum(x)
    yb = minimum(y)
    ye = maximum(y)

    # ensure arrays are of correct type
    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)
    zin = convert(Vector{Float64}, z)
    win = convert(Vector{Float64}, w)

    # return values
    nx = Array(Int32, 1)
    tx = Array(Float64, nxest)
    ny = Array(Int32, 1)
    ty = Array(Float64, nyest)
    c = Array(Float64, (nxest-kx-1) * (nyest-ky-1))
    fp = Array(Float64, 1)
    ier = Array(Int32, 1)

    # work arrays
    # Note: in lwrk1, x and y are swapped on purpose.
    lwrk1 = calc_surfit_lwrk1(m, ky, kx, nyest, nxest)
    lwrk2 = 1
    kwrk = m + (nxest - 2*kx - 1) * (nyest - 2*ky - 1)
    wrk1 = Array(Float64, lwrk1)
    wrk2 = Array(Float64, lwrk2)
    iwrk = Array(Int32, kwrk)

    ccall((:surfit_, ddierckx), Void,
          (Ptr{Int32}, Ptr{Int32},  # iopt, m
           Ptr{Float64}, Ptr{Float64},  # y, x
           Ptr{Float64}, Ptr{Float64},  # z, w
           Ptr{Float64}, Ptr{Float64},  # yb, ye
           Ptr{Float64}, Ptr{Float64},  # xb, xe
           Ptr{Int32}, Ptr{Int32}, Ptr{Float64},  # ky, kx, s
           Ptr{Int32}, Ptr{Int32}, Ptr{Int32},  # nyest, nxest, nmax
           Ptr{Float64},  # eps
           Ptr{Int32}, Ptr{Float64},  # ny, tx
           Ptr{Int32}, Ptr{Float64},  # nx, tx
           Ptr{Float64}, Ptr{Float64},  # c, fp
           Ptr{Float64}, Ptr{Int32},  # wrk1, lwrk1
           Ptr{Float64}, Ptr{Int32},  # wrk2, lwrk2
           Ptr{Int32}, Ptr{Int32}, Ptr{Int32}),   # iwrk, kwrk, ier
          &0, &m, yin, xin, zin, win, &yb, &ye, &xb, &xe, &ky, &kx,
          &@compat(Float64(s)), &nyest, &nxest, &nmax, &eps, ny, ty, nx, tx,
          c, fp, wrk1, &lwrk1, wrk2, &lwrk2, iwrk, &kwrk, ier)

    while ier[1] > 10
        # lwrk2 is too small, i.e., there is not enough workspace
        # for computing the minimal least-squares solution of a rank
        # deficient system of linear equations. ier gives the
        # requested value for lwrk2. Rerun with that value in "continue"
        # mode, with iopt = 1.
        lwrk2 = ier[1]
        resize!(wrk2, lwrk2)
        ccall((:surfit_, ddierckx), Void,
              (Ptr{Int32}, Ptr{Int32},  # iopt, m
               Ptr{Float64}, Ptr{Float64},  # y, x
               Ptr{Float64}, Ptr{Float64},  # z, w
               Ptr{Float64}, Ptr{Float64},  # yb, ye
               Ptr{Float64}, Ptr{Float64},  # xb, xe
               Ptr{Int32}, Ptr{Int32}, Ptr{Float64},  # ky, kx, s
               Ptr{Int32}, Ptr{Int32}, Ptr{Int32},  # nyest, nxest, nmax
               Ptr{Float64},  # eps
               Ptr{Int32}, Ptr{Float64},  # ny, tx
               Ptr{Int32}, Ptr{Float64},  # nx, tx
               Ptr{Float64}, Ptr{Float64},  # c, fp
               Ptr{Float64}, Ptr{Int32},  # wrk1, lwrk1
               Ptr{Float64}, Ptr{Int32},  # wrk2, lwrk2
               Ptr{Int32}, Ptr{Int32}, Ptr{Int32}),   # iwrk, kwrk, ier
              &1, &m, yin, xin, zin, win, &yb, &ye, &xb, &xe, &ky, &kx, &s,
              &nyest, &nxest, &nmax, &eps, ny, ty, nx, tx, c, fp,
              wrk1, &lwrk1, wrk2, &lwrk2, iwrk, &kwrk, ier)
    end

    if (ier[1] == 0 || ier[1] == -1 || ier[1] == -2)
        # good values, pass.
    elseif ier[1] < -2
        warn("""
        The coefficients of the spline returned have been
        computed as the minimal norm least-squares solution of a
        (numerically) rank deficient system. The rank is $(-ier[1]).
        The rank deficiency is $((nx[1]-kx-1)*(ny[1]-ky-1)+ier[1]).
        Especially if the rank deficiency is large the results may
        be inaccurate.""")
        # "The results could also seriously depend on the value of
        # eps" (not in message because eps is currently not an input)
    else
        error(_fit2d_messages[ier[1]])
    end

    # Resize output arrays to the size actually used.
    resize!(tx, nx[1])
    resize!(ty, ny[1])
    resize!(c, (nx[1] - kx - 1) * (ny[1] - ky - 1))

    return Spline2D(tx, ty, c, kx, ky, fp[1])
end

# Construct spline from data on a grid.
function Spline2D(x::AbstractVector, y::AbstractVector, z::AbstractMatrix;
                  kx::Int=3, ky::Int=3, s::Real=0.0)
    mx = length(x)
    my = length(y)
    assert(size(z, 1) == mx && size(z, 2) == my)

    mx > kx || error("length(x) must be greater than kx")
    my > ky || error("length(y) must be greater than ky")

    # Bounds
    xb = x[1]
    xe = x[end]
    yb = y[1]
    ye = y[end]
    nxest = mx+kx+1
    nyest = my+ky+1

    # ensure arrays are of correct type
    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)
    zin = convert(Matrix{Float64}, z)

    # Return values
    nx = Array(Int32, 1)
    tx = Array(Float64, nxest)
    ny = Array(Int32, 1)
    ty = Array(Float64, nyest)
    c = Array(Float64, (nxest-kx-1) * (nyest-ky-1))
    fp = Array(Float64, 1)
    ier = Array(Int32, 1)

    # Work arrays.
    # Note that in lwrk, x and y are swapped with respect to the Fortran
    # documentation. See "NOTE REGARDING ARGUMENT ORDER" above.
    lwrk = (4 + nyest * (mx+2*ky+5) + nxest * (2*kx+5) +
            my*(ky+1) + mx*(kx+1) + max(mx, nyest))
    wrk = Array(Float64, lwrk)
    kwrk = 3 + mx + my + nxest + nyest
    iwrk = Array(Int32, kwrk)

    ccall((:regrid_, ddierckx), Void,
          (Ptr{Int32},  # iopt
           Ptr{Int32}, Ptr{Float64},  # my, y
           Ptr{Int32}, Ptr{Float64},  # mx, x
           Ptr{Float64},  # z
           Ptr{Float64}, Ptr{Float64},  # yb, ye
           Ptr{Float64}, Ptr{Float64},  # xb, xe
           Ptr{Int32}, Ptr{Int32}, Ptr{Float64},  # ky, kx, s
           Ptr{Int32}, Ptr{Int32},  # nyest, nxest
           Ptr{Int32}, Ptr{Float64},  # ny, ty
           Ptr{Int32}, Ptr{Float64},  # nx, tx
           Ptr{Float64}, Ptr{Float64},  # c, fp
           Ptr{Float64}, Ptr{Int32},  # wrk, lwrk
           Ptr{Int32}, Ptr{Int32},  # iwrk, lwrk
           Ptr{Int32}),  # ier
          &0f0, &my, yin, &mx, xin, zin, &yb, &ye, &xb, &xe, &ky, &kx,
          &@compat(Float64(s)), &nyest, &nxest, ny, ty, nx, tx, c, fp,
          wrk, &lwrk, iwrk, &kwrk, ier)

    if !(ier[1] == 0 || ier[1] == -1 || ier[1] == -2)
        error(_fit2d_messages[ier[1]])
    end

    # Resize output arrays to the size actually used.
    resize!(tx, nx[1])
    resize!(ty, ny[1])
    resize!(c, (nx[1] - kx - 1) * (ny[1] - ky - 1))

    return Spline2D(tx, ty, c, kx, ky, fp[1])
end

# Evaluate spline at individual points
function evaluate(spline::Spline2D, x::AbstractVector, y::AbstractVector)
    m = length(x)
    assert(length(y) == m)

    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)

    ier = Array(Int32, 1)
    lwrk = spline.kx + spline.ky + 2
    wrk = Array(Float64, lwrk)
    z = Array(Float64, m)

    ccall((:bispeu_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # ty, ny
           Ptr{Float64}, Ptr{Int32},  # tx, nx
           Ptr{Float64},  # c
           Ptr{Int32}, Ptr{Int32},  # ky, kx
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # y, x, z
           Ptr{Int32},  # m
           Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),  # wrk, lwrk, ier
          spline.ty, &length(spline.ty),
          spline.tx, &length(spline.tx),
          spline.c, &spline.ky, &spline.kx, yin, xin, z, &m,
          wrk, &lwrk, ier)

    ier[1] == 0 || error(_eval2d_message)

    return z
end

function evaluate(spline::Spline2D, x::Real, y::Real)
    ier = Array(Int32, 1)
    lwrk = spline.kx + spline.ky + 2
    wrk = Array(Float64, lwrk)
    z = Array(Float64, 1)
    ccall((:bispeu_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # ty, ny
           Ptr{Float64}, Ptr{Int32},  # tx, nx
           Ptr{Float64},  # c
           Ptr{Int32}, Ptr{Int32},  # ky, kx
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # y, x, z
           Ptr{Int32},  # m
           Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),  # wrk, lwrk, ier
          spline.ty, &length(spline.ty),
          spline.tx, &length(spline.tx),
          spline.c, &spline.ky, &spline.kx, &y, &x, z, &1,
          wrk, &lwrk, ier)

    ier[1] == 0 || error(_eval2d_message)
    return z[1]
end

# Evaluate spline on the grid spanned by the input arrays.
function evalgrid(spline::Spline2D, x::AbstractVector, y::AbstractVector)
    mx = length(x)
    my = length(y)

    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)

    lwrk = mx*(spline.kx + 1) + my*(spline.ky + 1)
    wrk = Array(Float64, lwrk)
    kwrk = mx + my
    iwrk = Array(Int32, kwrk)
    ier = Array(Int32, 1)
    z = Array(Float64, mx, my)

    ccall((:bispev_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # ty, ny
           Ptr{Float64}, Ptr{Int32},  # tx, nx
           Ptr{Float64},  # c
           Ptr{Int32}, Ptr{Int32},  # ky, kx
           Ptr{Float64}, Ptr{Int32},  # y, my
           Ptr{Float64}, Ptr{Int32},  # x, mx
           Ptr{Float64},  # z
           Ptr{Float64}, Ptr{Int32},  # wrk, lwrk
           Ptr{Int32}, Ptr{Int32},  # iwrk, kwrk
           Ptr{Int32}),  # ier
          spline.ty, &length(spline.ty),
          spline.tx, &length(spline.tx),
          spline.c, &spline.ky, &spline.kx, yin, &my, xin, &mx, z,
          wrk, &lwrk, iwrk, &kwrk, ier)

    ier[1] == 0 || error(_eval2d_message)

    return z
end

# call synonyms on v0.4+
@compat if VERSION >= v"0.4.0"
    (spl::Spline1D)(x::Real) = evaluate(spl, x)
    (spl::Spline1D)(x::AbstractVector) = evaluate(spl, x)
    (spl::Spline2D)(x::Real, y::Real) = evaluate(spl, x, y)
    (spl::Spline2D)(x::AbstractVector, y::AbstractVector) =
        evaluate(spl, x, y)
end

end # module


text:
using MNIST, Nabla

""" Implementation of the Adam optimiser. """
mutable struct Adam{T<:AbstractArray}
    α::Float64
    β1::Float64
    β2::Float64
    m::T
    v::T
    β1_acc::Float64
    β2_acc::Float64
    ϵ::Float64
end
function Adam(θ0::AbstractArray, α::Float64, β1::Float64, β2::Float64, ϵ::Float64)
    return Adam(α, β1, β2, zeros(θ0), zeros(θ0), β1, β2, ϵ)
end

function iterate!(θ::AbstractArray{Float64}, ∇θ::AbstractArray{Float64}, opt::Adam)
    m, v, α, β1, β2, ϵ, β1_acc, β2_acc =
        opt.m, opt.v, opt.α, opt.β1, opt.β2, opt.ϵ, opt.β1_acc, opt.β2_acc
    @inbounds for n in eachindex(θ)
        m[n] = β1 * m[n] + (1.0 - β1) * ∇θ[n]
        v[n] = β2 * v[n] + (1.0 - β2) * ∇θ[n]^2
        m̂ = m[n] / (1 - β1_acc)
        v̂ = v[n] / (1 - β2_acc)
        θ[n] = θ[n] + α * m̂ / (sqrt(v̂) + ϵ)
    end
    opt.β1_acc *= β1
    opt.β2_acc *= β2
end

function to1hot(y_::Vector)
    y = Matrix{Int}(undef, 10, length(y_))
    for n in eachindex(y_)
        for j in 1:10
            y[j, n] = y_[n] == j - 1 ? 1 : 0
        end
    end
    return y
end
accuracy(Y, Ypr) = mean(all(Y .== Ypr, dims=1))
@inline logistic(x) = 1 ./ (1 .+ exp.(-x))

"""
    mlp_log_joint(
        X::AbstractMatrix{T} where T<:Real,
        Y::AbstractMatrix{T} where T<:Real,
        W::NTuple{N, T} where T<:AbstractMatrix{V} where V<:Real,
        b::NTuple{N, T} where T<:AbstractVector{V} where V<:Real,
        λ::AbstractFloat)

Compute the log joint probability of the data `(X, Y)` given weights `W`, biases `b` and
precision for the isotropic Gaussian prior on the weights and biases `λ` (equivalent to
quadratic weight decay) and a Bernoulli likelihood (equivalent to the cross-entropy).
"""
function mlp_log_joint(X, Y, W, b, λ)

    # Compute the log prior probability of the weights. We have assumed an isotropic
    # Gaussian prior over all of the parameters. Note that since we are not trying to learn
    # λ (it is fixed) we can neglect the normalising constant Z = 0.5 * log(2πλ).
    log_prior = 0.0
    for n in 1:length(W)
        log_prior += sum(abs2, W[n]) + sum(abs2, b[n])
    end
    log_prior *= -0.5 * λ

    # Compute the output of the MLP.
    f = logistic.(apply_transforms(X, W, b, tanh))
    f = f ./ sum(f, dims=1)

    # Compute the log likelihood of the observations given the outputs of the MLP. We assume
    ϵ = 1e-15
    log_lik = sum(Y .* log.(f .+ ϵ) .+ (1 - Y) .* log.((1 + ϵ) .- f))

    # Return the joint and the predictive distribution over the labels.
    return log_lik + log_prior, f
end

apply_transforms(X, W::Tuple{Vararg{Any, 1}}, b::Tuple{Vararg{Any, 1}}, f) =
    f.(W[1] * X .+ b[1])
apply_transforms(X, W::Tuple{Vararg{Any, N}}, b::Tuple{Vararg{Any, N}}, f) where N =
    f.(W[N] * apply_transforms(X, W[1:N-1], b[1:N-1], f) .+ b[N])

# A simple Multilayer Feedforward Neural Network (Multi-layer Perceptron (MLP)) example for
# classifying the MNIST data set.
function demo_mlp(itrs::Int, sz::Int)

    # Load the data.
    println("Loading data.")
    xtr, ytr_ = traindata()
    xte, yte_ = testdata()

    # Convert to 1-hot label encoding.
    ytr, yte = to1hot(ytr_), to1hot(yte_)

    # Initialise parameters.
    println("Initialising parameters.")
    d0, d1, d2, d3, λ = size(xtr, 1), 500, 300, size(ytr, 1), 1e-3
    W_ = (0.1 * randn(d1, d0), 0.1 * randn(d2, d1), 0.1 * randn(d3, d2))
    b_ = (0.1 * randn(d1), 0.1 * randn(d2), 0.1 * randn(d3))

    # Initialise the Adam optimiser.
    α, β1, β2, ϵ = 1e-3, 0.9, 0.999, 1e-8
    optW, optb = Adam.(W_, α, β1, β2, ϵ), Adam.(b_, α, β1, β2, ϵ)

    # Iterate to learn the parameters.
    println("Starting learning.")
    scal = size(xtr, 2) / sz
    for itr in 1:itrs

        # Pick the mini batch.
        idx = rand(eachindex(ytr_), sz)
        xtr_batch = view(xtr, :, idx)
        ytr_batch = view(ytr, :, idx)

        # Initialise computational graph.
        tape = Tape()
        W, b = Leaf.(tape, W_), Leaf.(tape, b_)

        # Compute the log marginal probability.
        logp, f = mlp_log_joint(xtr_batch, ytr_batch, W, b, λ)
        ∇logp = ∇(logp)

        # Compute most probably classification for each observation in the batch.
        ypr = zeros(d3, sz)
        for n in 1:sz
            ypr[argmax(Nabla.unbox(f)[:, n]), n] = 1.
        end
        acc = accuracy(ytr_batch, ypr)

        # Update the parameters using AdaGrad by indexing into the ∇logp tape.
        iterate!.(W_, getindex.(∇logp, W), optW)
        iterate!.(b_, getindex.(∇logp, b), optb)
        println("logp is $(Nabla.unbox(logp)) at iteration $itr. Mean loglik is ",
            "$(Nabla.unbox(logp) / size(xtr, 2)). Accuracy is $acc")
    end
end

ext: jl


text:
#Copyright Nathan Parrish 2018
#University of Colorado, Boulder
#Colorado Center for Astrodynamics Research
#Celestial and Spaceflight Mechanics Laboratory
#
#CRTBP propagator, with thrust given

function CRTBP_prop_EP_deriv(t, state, MU, DU, TU, Isp, control, time_direction)
    #pull out from state vector:
    x = state[1];
    y = state[2];
    z = state[3];
    xdot = state[4];
    ydot = state[5];
    zdot = state[6];

    if length(state) == 7
        m = state[7];
    else
        m = 1000.0; #default mass (kg)
    end

    #calculate some stuff:
    r1 = sqrt((x+MU)^2 + y^2 + z^2);
    r2 = sqrt((x+MU-1)^2 + y^2 + z^2);

    #calculate r1^3 and r2^3 once each
    r1_3 = r1.^3;
    r2_3 = r2.^3;

    #thrust acceleration magnitude (DU/TU^2):
    T_mag = norm(control)/m/1e3 * TU^2 / DU; # (kg*m/s^2) -> (km/s^2) -> (DU/TU^2) acceleration

    #thrust acceleration vector:
    if norm(control) == 0
        T = control
    else
        T = control ./ norm(control) * T_mag
    end

    g0 = 9.81; #m/s
    mdot = - time_direction * norm(control)/(Isp*g0) * TU; #mass flow rate (kg/TU)

    #allows propagation either forward or backward in the CRTBP:
    omega = time_direction; #+1 for forwards, -1 for backwards

    #calculate accelerations:
    xddot = -(1-MU)*(x+MU)/r1_3 - MU*(x-1+MU)/r2_3 + 2*omega*ydot + x + T[1];
    yddot = -(1-MU)*y/r1_3 - MU*y/r2_3 - 2*omega*xdot + y + T[2];
    zddot = -(1-MU)*z/r1_3 - MU*z/r2_3 + T[3];

    #put back into dstate vector:
    if length(state) == 7
        dstate = [xdot; ydot; zdot; xddot; yddot; zddot; mdot];
    else
        dstate = [xdot; ydot; zdot; xddot; yddot; zddot]; #constant mass
    end

    #Outputs:
    dstate
end

#In-place update of derivatives:
function CRTBP_prop_EP_deriv!(t, state, params, dstate)

    (MU, DU, TU, Isp, control, time_direction) = params

    #pull out from state vector:
    x = state[1];
    y = state[2];
    z = state[3];
    xdot = state[4];
    ydot = state[5];
    zdot = state[6];

    if length(state) == 7
        m = state[7];
    else
        m = 1000.0; #default mass (kg)
    end

    #calculate some stuff:
    r1 = sqrt((x+MU)^2 + y^2 + z^2);
    r2 = sqrt((x+MU-1)^2 + y^2 + z^2);

    #calculate r1^3 and r2^3 once each
    r1_3 = r1.^3;
    r2_3 = r2.^3;

    #thrust acceleration magnitude (DU/TU^2):
    T_mag = norm(control)/m/1e3 * TU^2 / DU; # (kg*m/s^2) -> (km/s^2) -> (DU/TU^2) acceleration

    #thrust acceleration vector:
    if norm(control) == 0
        T = control
    else
        T = control ./ norm(control) * T_mag
    end

    g0 = 9.81; #m/s
    mdot = - time_direction * norm(control)/(Isp*g0) * TU; #mass flow rate (kg/TU)

    #allows propagation either forward or backward in the CRTBP:
    omega = time_direction; #+1 for forwards, -1 for backwards

    #calculate accelerations:
    dstate[1] = xdot
    dstate[2] = ydot
    dstate[3] = zdot
    dstate[4] = -(1-MU)*(x+MU)/r1_3 - MU*(x-1+MU)/r2_3 + 2*omega*ydot + x + T[1]
    dstate[5] = -(1-MU)*y/r1_3 - MU*y/r2_3 - 2*omega*xdot + y + T[2];
    dstate[6] = -(1-MU)*z/r1_3 - MU*z/r2_3 + T[3];

    #put back into dstate vector:
    if length(state) == 7
        dstate[7] = mdot
    end

    #Outputs:
    nothing
end

#for new version of DifferentialEquations, with control computed as a function
#of interpolated costates:
function CRTBP_prop_EP_NNControl_deriv!(dstate, state, params, t)

    (MU, DU, TU, Isp, time_direction, λv_itp, thrustLimit, p, rho) = params

    #λv_itp is an interpolation object for the control. It is prepared like this:
    # using Interpolations
    # λv_itp = scale(interpolate(λv_all', (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), t_TU, 1:3)
    # (from https://github.com/JuliaMath/Interpolations.jl)
    #where t_sec has to be an evenly-spaced time grid.

    #interpolate control from given history:
    λv = [λv_itp[t, 1], λv_itp[t, 2], λv_itp[t, 3]]

    accelLimit = thrustLimit / mass / 1e3 * TU^2 / DU #(kg*m/s^2) -> (km/s^2) -> (DU/TU^2) acceleration

    #Thrust magnitude from control law:
    if p == 0
      #Force thrust to always be on, in the optimal direction. This limits how big
      #of a search we need to make for the neural net stuff.
      umag = accelLimit #N

    elseif p == 1 #minimum mass problem
        g = norm(λv) - 1
        umag = 1/2 * (1 + tanh( g / (2*rho))) * accelLimit

    elseif p > 1 #p is in the range (1, 2]
        umag = (1/p * norm(λv))^(1 / (p-1));

        if umag > accelLimit
          umag = accelLimit
        end
    else
      error("Invalid value of p!")
    end

    #Control vector, scaled properly:
    # control_accel = -umag * λv ./ norm(λv) / mass / 1e3 * TU^2 / DU #(kg*m/s^2) -> (km/s^2) -> (DU/TU^2) acceleration
    control_accel = -umag * λv ./ norm(λv) #DU/TU^2

    if isnan(control_accel[1])
        control_accel[:] = 0.
    end

    #pull out from state vector:
    x = state[1];
    y = state[2];
    z = state[3];
    xdot = state[4];
    ydot = state[5];
    zdot = state[6];

    if length(state) == 7
        m = state[7];
    else
        m = 1000.0; #default mass (kg)
    end

    #calculate some stuff:
    r1 = sqrt((x+MU)^2 + y^2 + z^2);
    r2 = sqrt((x+MU-1)^2 + y^2 + z^2);

    #calculate r1^3 and r2^3 once each
    r1_3 = r1.^3;
    r2_3 = r2.^3;

    g0 = 9.81; #m/s
    mdot = - time_direction * norm(umag)/(Isp*g0) * TU; #mass flow rate (kg/TU)

    #allows propagation either forward or backward in the CRTBP:
    omega = time_direction; #+1 for forwards, -1 for backwards

    #calculate accelerations:
    dstate[1] = xdot
    dstate[2] = ydot
    dstate[3] = zdot
    dstate[4] = -(1-MU)*(x+MU)/r1_3 - MU*(x-1+MU)/r2_3 + 2*omega*ydot + x + control_accel[1]
    dstate[5] = -(1-MU)*y/r1_3 - MU*y/r2_3 - 2*omega*xdot + y + control_accel[2];
    dstate[6] = -(1-MU)*z/r1_3 - MU*z/r2_3 + control_accel[3];

    #if mass is variable:
    if length(state) == 7
        dstate[7] = mdot
    end

    #Outputs:
    nothing
end

max_stars_repo_name: travelingspaceman/LowThrustOpt

max_issues_repo_name: travelingspaceman/LowThrustOpt

max_forks_repo_name: travelingspaceman/LowThrustOpt


text:
@testset "Bus Constructors" begin
    tBus = Bus()
    tLoadZones = LoadZones()
end

@testset "Generation Constructors" begin
    tEconThermal = EconThermal()
    @test tEconThermal isa PowerSystems.Component
    tTechThermal = TechThermal()
    @test tTechThermal isa PowerSystems.Component
    tThermalGen = ThermalDispatch()
    @test tThermalGen isa PowerSystems.Component
    tThermalGenSeason = ThermalGenSeason()
    @test tThermalGenSeason isa PowerSystems.Component
    tTechHydro = TechHydro()
    @test tTechHydro isa PowerSystems.Component
    tEconHydro = EconHydro()
    @test tEconHydro isa PowerSystems.Component
    tHydroFix = HydroFix()
    @test tHydroFix isa PowerSystems.Component
    tHydroCurtailment = HydroCurtailment()
    @test tHydroCurtailment isa PowerSystems.Component
    tHydroStorage = HydroStorage()
    @test tHydroStorage isa PowerSystems.Component
    tTechRenewable = TechRenewable()
    @test tTechRenewable isa PowerSystems.Component
    tEconRenewable = EconRenewable()
    @test tEconRenewable isa PowerSystems.Component
    tRenewableFix = RenewableFix()
    @test tRenewableFix isa PowerSystems.Component
    tRenewableFullDispatch = RenewableFullDispatch()
    @test tRenewableFullDispatch isa PowerSystems.Component
    tRenewableCurtailment = RenewableCurtailment()
    @test tRenewableCurtailment isa PowerSystems.Component
end

@testset "Storage Constructors" begin
    tStorage = GenericBattery()
    @test tStorage isa PowerSystems.Component
end

@testset "Load Constructors" begin
    tPowerLoad = PowerLoad()
    @test tPowerLoad isa PowerSystems.Component
    tPowerLoadPF = PowerLoadPF()
    @test tPowerLoadPF isa PowerSystems.Component
    tPowerLoad = PowerLoad("init", true, Bus(), 0.0, 0.0)
    @test tPowerLoad isa PowerSystems.Component
    tPowerLoadPF = PowerLoadPF("init", true, Bus(), 0.0, 1.0)
    @test tPowerLoadPF isa PowerSystems.Component
    tLoad = InterruptibleLoad()
    @test tLoad isa PowerSystems.Component
end

@testset "Branch Constructors" begin
    tLine = Line()
    @test tLine isa PowerSystems.Component
    tMonitoredLine = MonitoredLine()
    @test tMonitoredLine isa PowerSystems.Component
    tHVDCLine = HVDCLine()
    @test tHVDCLine isa PowerSystems.Component
    tVSCDCLine = VSCDCLine()
    @test tVSCDCLine isa PowerSystems.Component
    tTransformer2W = Transformer2W()
    @test tTransformer2W isa PowerSystems.Component
    tTapTransformer = TapTransformer()
    @test tTapTransformer isa PowerSystems.Component
    tPhaseShiftingTransformer = PhaseShiftingTransformer()
    @test tPhaseShiftingTransformer isa PowerSystems.Component
end

@testset "Product Constructors" begin
    tProportionalReserve = ProportionalReserve()
    @test tProportionalReserve isa PowerSystems.Service
    tStaticReserve = StaticReserve()
    @test tStaticReserve isa PowerSystems.Service
end


text:
print("matrixops, ")
## this test for typos, mainly, not for correct functioning...

for t in (Float32, Float64)
    a = rand(t, 10, 10)
    b = rand(t, 10, 5)
    c = rand(t, 10, 10)
    apd = a' * a
    asym = a' + a
    luf = lufact(a)
    chf = cholfact(apd)
    qrf = qrfact(a)
    eif = eigfact(a)
    eigvals(a)
    hsf = hessfact(a)
    scf = schurfact(a)
    scf1 = schurfact(a, c)
    svdf = svdfact(a)
    svdvals(a)
    d = diag(a)
    diagm(d)
    cond(a)
    if VERSION ≥ v"0.4-dev"
        nullspace(a)
    else
        null(a)
    end
    kron(a, a)
    linreg(a[:,1], b[:,1])
    lyap(a, c)
    sylvester(a, c, a)
    isposdef(apd)
end


text:
using ArgParse
using Formatting
using NCDatasets
using JSON

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--POP2-monthly-profile"
            help = "POP2 ocean output monthly profile used as forcing. It is currently assume the data is centered at the middle of each month."
            arg_type = String
            required = true

        "--n-layers"
            help = "The number of vertical layers that is going to be output."
            arg_type = Int64
            required = true

        "--domain-file"
            help = "The CESM domain file that is consistent with the `--POP2-monthly-profile`."
            arg_type = String
            required = true

        "--output-dir"
            help = "Output directory"
            arg_type = String
            required = true

        "--output-forcing-file"
            help = "Output forcing file name."
            arg_type = String
            default = ""

    end

    return parse_args(s)
end

parsed = parse_commandline()

JSON.print(parsed, 4)

mkpath(parsed["output-dir"])

# First, determine output file name
if parsed["output-forcing-file"] == ""
    parsed["output-forcing-file"] = format("{:s}.{:d}-layers.nc", splitext(basename(parsed["POP2-monthly-profile"]))[1], parsed["n-layers"])
end

output_file = joinpath(parsed["output-dir"], parsed["output-forcing-file"])

println("# Output forcing file is going to be: ", output_file)

# Second, check if vairables exist
vars_check = ["z_t", "z_w_top", "z_w_bot", "TEMP", "SALT", "SHF", "SHF_QSW", "SFWF", "HBLT", "TAUX", "TAUY", "time_bound"]

println("# Check if the following variables exist: ", join(vars_check, ", "))

Dataset(parsed["POP2-monthly-profile"], "r") do ds
    local missing_vars = []
    for varname in vars_check
        if ! haskey(ds, varname)
            push!(missing_vars, varname)
        end
    end

    if length(missing_vars) != 0
        println("Error: The following variables are missing: ", join(missing_vars, ", "))
        throw(ErrorException("Please check variables in the input file."))
    else
        println("All the variables are found.")
    end

    if ds.dim["time"] != 12
        throw(ErrorException("I expect dimension of time = 12."))
    end

end

function pleaseRun(cmd)
    println(">> ", string(cmd))
    run(cmd)
end

pleaseRun(`
ncks -O -F -d z_t,1,$(parsed["n-layers"])
        -d z_w_top,1,$(parsed["n-layers"])
        -d z_w_bot,1,$(parsed["n-layers"])
       $(parsed["POP2-monthly-profile"]) $(output_file)
`)

pleaseRun(`
ncap2 -O -v -s 'time_bound=time_bound;'
            -s 'z_t=-z_t/100.0;'
            -s 'z_w_top=-z_w_top/100.0;'
            -s 'z_w_bot=-z_w_bot/100.0;'
            -s 'TEMP=TEMP;'
            -s 'SALT=SALT;'
            -s 'SWFLX=-SHF_QSW;'
            -s 'NSWFLX=-(SHF-SHF_QSW);'
            -s 'VSFLX=SFWF;'
            -s 'HMXL=HBLT/100.0;'
            -s 'TAUX=TAUX/10.0;'
            -s 'TAUY=TAUY/10.0;'
             $(output_file) $(output_file)
`)

pleaseRun(`
ncatted -a units,time,m,c,"days since 0001-01-01 00:00:00"
        -a units,time_bound,m,c,"days since 0001-01-01 00:00:00"
        $(output_file)
`)

pleaseRun(`
ncrename -d nlat,Ny -d nlon,Nx -d z_t,Nz $(output_file)
`)

Dataset(output_file, "r") do ds

    local TEMP  = permutedims(nomissing(ds["TEMP"][:, :, :, 1],  NaN), [3, 1, 2])
    local Nz, Nx, Ny = size(TEMP)

    z_w_top = nomissing(ds["z_w_top"][:], NaN)
    z_w_bot = nomissing(ds["z_w_bot"][:], NaN)

    local z_w = zeros(Float64, length(z_w_top)+1)
    z_w[1:end-1] = z_w_top
    z_w[end] = z_w_bot[end]

    Dataset(joinpath(parsed["output-dir"], "z_w.nc"), "c") do _ds

        defDim(_ds, "Nzp1", length(z_w))
        defVar(_ds, "z_w", z_w, ("Nzp1", ), ; attrib = Dict(
            "long_name" => "Vertical coordinate on W-grid",
            "units"     => "m",
        ))

    end

    mask_T = zeros(Float64, Nz, Nx, Ny)
    mask_T[isfinite.(TEMP)] .= 1.0
    valid_idx = mask_T .== 1.0
    local Nz_bot = convert(Array{Int64}, sum(mask_T, dims=1)[1, :, :])

    Dataset(joinpath(parsed["output-dir"], "Nz_bot.nc"), "c") do _ds
        defDim(_ds, "Nx", Nx)
        defDim(_ds, "Ny", Ny)
        defVar(_ds, "Nz_bot", Nz_bot, ("Nx", "Ny", ), ; attrib = Dict(
            "long_name" => "z-grid idx of the deepest cell",
            "units"     => "none",
        ))

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module Compositions

using Nemo, IterTools, Triangles

export ModuleCompositions
export T097805, L097805, V097805, M097805
# The module SetPartitionsMType contains : P097805, L097805, TL097805,

"""

* T097805, L097805, V097805, M097805
"""
const ModuleCompositions = ""

"""

Recurrence for `A097805`, the compositions of ``n`` with ``k`` parts.
"""
function R097805(n, k, prevrow::Function)
    k == 0 && return fmpz(k^n)
    prevrow(k - 1) + prevrow(k)
end

"""

Iterates over the first ``n`` rows of A097805.
"""
T097805(n) = RecTriangle(n, R097805)

"""

Lists the first ``n`` rows of A097805 by concatinating. This is the format for submissions to the OEIS.
"""
L097805(n) = vcat(T097805(n)...)

"""

Return the triangular array as a square matrix.
"""
M097805(dim) = fromΔ(T097805(dim))

"""

Return row ``n`` of A097805 based on the iteration T097805(n).
"""
V097805(n) = nth(T097805(n + 1), n + 1)

#START-TEST-########################################################

# Return row n of A097805 based on a closed formula.
VN097805(n) = n == 0 ? fmpz[1] : [Nemo.binomial(ZZ(n - 1), ZZ(k - 1)) for k ∈ 0:n]
VJ097805(n) = n == 0 ? BigInt[1] : [binomial(BigInt(n - 1), k - 1) for k ∈ 0:n]

function test() end

function demo()
    println("\nIterates over the first n rows of A097805.")
    for r ∈ T097805(9)
        println(r)
    end

    println("\nLists the first 9 rows of L097805 by concatinating.")
    println(L097805(9))

    #ShowAsMatrix(L097805(9))
    #println(M097805(9))

    println("\nReturns row 8 of A097805 based on the iteration T097805(8).")
    println(V097805(8))

    println("\nReturns row 8 of A097805 based on a closed formula.")
    println(VN097805(8))

    println("\nReturns row n of A097805 based on the iteration T097805(n).")
    for n ∈ 0:8
        println(n, ": ", V097805(n))
    end

    println("\nBenchmark the construction of the first 500 rows of A097805 based on Iteration.\n")
    #  0.111066 seconds (253.51 k allocations: 7.795 MiB)
    GC.gc()
    @time L097805(500)

    println("\nBenchmark the construction of the first 500 rows of A097805 based on closed formula.")
    # Result: 0.479466 seconds (1.72 M allocations: 21.592 MiB)
    FList(n) = vcat([VN097805(k) for k ∈ 0:n-1]...)
    GC.gc()
    @time FList(500)

    println("\nBenchmark the construction of the first 500 rows of A097805 based on closed formula.")
    # 0.476226 seconds (1.72 M allocations: 21.592 MiB)
    FList(n) = vcat([VJ097805(k) for k ∈ 0:n-1]...)
    GC.gc()
    @time FList(500)

    println("\nResume: Construction based on Iteration is about 4 times faster")
    println("than based on the closed formula and uses much less memory.")
end

"""

L097805(500) :: 0.111066 seconds (253.51 k allocations: 7.795 MiB)
"""
function perf()
    GC.gc()
    @time L097805(500)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, following here
set_line_width(cr, 6);

rectangle(cr, 12, 12, 232, 70);
new_sub_path(cr); arc(cr, 64, 64, 40, 0, 2*pi);
new_sub_path(cr); arc_negative(cr, 192, 64, 40, 0, -2*pi);

set_fill_type(cr, Cairo.CAIRO_FILL_RULE_EVEN_ODD); # should be set_fill_rule
set_source_rgb(cr, 0, 0.7, 0); fill_preserve(cr);
set_source_rgb(cr, 0, 0, 0); stroke(cr);

translate(cr, 0, 128);
rectangle(cr, 12, 12, 232, 70);
new_sub_path(cr); arc(cr, 64, 64, 40, 0, 2*pi);
new_sub_path(cr); arc_negative(cr, 192, 64, 40, 0, -2*pi);

set_fill_type(cr, Cairo.CAIRO_FILL_RULE_WINDING);
set_source_rgb(cr, 0, 0, 0.9); fill_preserve(cr);
set_source_rgb(cr, 0, 0, 0); stroke(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_fill_style.png");


text:
# This file was generated, do not modify it. # hide
using PyPlot
figure(figsize=(8,6))
plot(X.Volume)
xlabel("Tick number", fontsize=14)
ylabel("Volume", fontsize=14)
xticks(fontsize=12)
yticks(fontsize=12)

savefig(joinpath(@OUTPUT, "ISL-lab-4-volume.svg")) # hide

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
module TimeModels

using Base.Dates
using Distributions
using StatsBase
using TimeSeries
using Optim
using NLopt
using Compat

import Base: show

export
  # General state space model
  StateSpaceModel, simulate,

  #Kalman
  loglikelihood, kalman_filter, KalmanFiltered,
  kalman_smooth, KalmanSmoothed, KalmanSmoothedMinimal,

  # Parameter fitting
  ParametrizedMatrix,
  parametrize_full, parametrize_diag, parametrize_none,
  ParametrizedSSM, SSMParameters, fit,

  # ARIMA exports
  arima_statespace,
  arima,

  # GARCH exports
  garchFit,
  predict,

  # diagnostic tests exports
  jbtest,

  # utilities
  em_checkmodel

# Core functionality
include("statespacemodel.jl")
include("kalman_filter.jl")
include("kalman_smooth.jl")
include("parameter_estimation.jl")

# Model specifications
include("ARIMA.jl")
include("GARCH.jl")

# Tests
include("diagnostic_tests.jl")

include("utilities.jl")

end


text:
#!/usr/bin/env julia
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Sanity checking for the parsing modules.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-05-07

include("../hdr.jl")

require("conllx.jl")
require("dep.jl")
require("hybrid.jl")
require("parse.jl")

using CoNLLX
using DepGraph
using Hybrid
using Parse

# Parse using an oracle and ensure that we get the gold annotations back.
open(talpath, "r") do data_f
    coder = Coder()
    sentnum = 0
    for goldsent in conllxparse(data_f, useproj=true)
        sentnum += 1
        try
            predsent = parse(goldsent, oracle, coder)
            @test isequal(predsent, goldsent)
        catch e
            println("Sentence: $sentnum")
            println(goldsent)
            rethrow(e)
        end
    end
end

# Parse using a random oracle and ensure that we do not crash.
open(talpath, "r") do data_f
    coder = Coder()
    for goldsent in conllxparse(data_f, useproj=true)
        parse(goldsent, randoracle, coder)
    end
end

# Verify that each transition undo action is valid.
open(debugpath, "r") do debug_f
    transs = transitions()
    coder = Coder()
    for goldsent in conllxparse(debug_f, useproj=true)
        conf = config(goldsent, coder)
        # Walk some steps along the gold path to test different regions.
        for _ in rand(0:length(goldsent))
            if !isterminal(conf)
                apply!(conf, oracle(conf))
            end
        end

        for trans in transs
            if !islegal(conf, trans)
                continue
            end

            before = deepcopy(conf)
            undo = apply!(conf, trans)
            undo!(conf, undo)
            @test isequal(conf, before)
        end
    end
end

open(debugpath, "r") do debug_f
    goldsents = collect(Sentence, conllxparse(debug_f, useproj=true))

    # Can we memorise a very simple sentence?
    let
        sent = goldsents[1:1]
        scores = Any[]
        for model in train(sent, epochs=7)
            push!(scores, evaluate(model, sent))
        end
        @test_approx_eq scores[end] 1.0
    end

    # Can we memorise a few short sentences?
    let
        scores = Any[]
        for model in train(goldsents, epochs=42)
            push!(scores, evaluate(model, goldsents))
        end
        @test_approx_eq scores[end] 1.0
    end
end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
"""
Write and balance chemical equations elegantly and efficiently.
"""
module ChemEquations

using LinearAlgebraX: I, nullspacex, IntegerX
using DocStringExtensions
@template TYPES =
    """
    $TYPEDEF
    $TYPEDFIELDS
    $DOCSTRING
    """
@template (FUNCTIONS, METHODS) =
    """
    $TYPEDSIGNATURES
    $DOCSTRING
    """

export Compound, ChemEquation,
    @ce_str, @cc_str, ==, string, show,
    compounds, elements, hascharge,
    equationmatrix, balancematrix, balance

include("compound.jl")
include("chemequation.jl")
include("balance.jl")

end # module


text:
#input: model
#run WAVEFORM package for seismic full-waveform inversion for 1 or a few iterations
#output: updated model

#This package only runs in Julia 0.6
#This example used the Waveform package, which is in julia 0.6 and uses the v06 branch of JOLI (https://github.com/slimgroup/JOLI.jl)
#untill everyting is updated to current julia version, you need to got to (on linux) /.julia/v0.6/JOLI and use: git checkout v06

#this example also uses two scripts from the JUDI package (see below)

using Waveform
using JOLI
using JLD

include(string(Pkg.dir("JUDI"),"/src/Optimization/SPGSlim.jl"))
include(string(Pkg.dir("JUDI"),"/src/Optimization/OptimizationFunctions.jl"))

#load model
tmp = load("vel_model.jld")
vel_model = tmp["vel_model"]

#load other information for waveform inversion (frequencies, sources, options, bound constraints)
FWI_info = load("FWI_info.jld")
#unpack
bound_min       = FWI_info["bound_min"]
bound_max       = FWI_info["bound_max"]
nsrc            = FWI_info["nsrc"]
nfreq           = FWI_info["nfreq"]
Q               = FWI_info["Q"]
max_func_evals  = FWI_info["max_func_evals"]
v_true          = FWI_info["v_true"]
Q               = FWI_info["Q"]
n               = FWI_info["n"]
d               = FWI_info["d"]
o               = FWI_info["o"]
t0              = FWI_info["t0"]
f0              = FWI_info["f0"]
unit            = FWI_info["unit"]
freqs           = FWI_info["freqs"]
xsrc            = FWI_info["xsrc"]
ysrc            = FWI_info["ysrc"]
zsrc            = FWI_info["zsrc"]
xrec            = FWI_info["xrec"]
yrec            = FWI_info["yrec"]
zrec            = FWI_info["zrec"]
npml            = FWI_info["npml"]

model = Model{Int64,Float64}(n,d,o,t0,f0,unit,freqs,xsrc,ysrc,zsrc,xrec,yrec,zrec);
(z,x) = odn_to_grid(o,d,n);
## Size/spacing of computational domain
comp_n = n;
comp_d = d;
comp_o = o;

# Linear solver options
lsopts = LinSolveOpts(solver=:lufact);

# Misfit function for the objective
misfit = least_squares;

# PDE scheme
scheme = Waveform.helm2d_chen9p;

# Whether to remove elements in the PML when restricting to the model domain (if false, will use stacking)
cut_pml = true;

# If true, use an implicit matrix representation (false for 2D, the matrices are small enough to form explicitly)
implicit_matrix = false;

# Binary mask of the sources and frequencies to compute
srcfreqmask = trues(nsrc,nfreq);

# Computational options type
opts = PDEopts{Int64,Float64}(scheme,comp_n,comp_d,comp_o,cut_pml,implicit_matrix,npml,misfit,srcfreqmask,lsopts);

#frequency partition
size_freq_batch=length(freqs)
overlap = 1
freq_partition = partition(nfreq,size_freq_batch,overlap)

#if data was not generated previously, generate it here
if isfile("D.jld")==true
  tmp = load("D.jld")
  D = tmp["D"]
else
  D = forw_model(v_true,Q,model,opts)
  save("D.jld","D",D)#save data
end

#proj! = (xproj,x)->project_bounds!(x,bound_min,bound_max,xproj);
#proj! = in -> 1.0.*in
proj! = in -> min.(max.(in,bound_min),bound_max)

options_SPG = spg_options(verbose=3, maxIter=max_func_evals, memory=3,suffDec=1f-6)
options_SPG.progTol=1e-6
options_SPG.testOpt=false
options_SPG.interp=2

for j in 1:size(freq_partition,2)
    fbatch = freq_partition[:,j]
    srcfreqmask = falses(nsrc,nfreq)
    srcfreqmask[:,fbatch] = true
    opts.srcfreqmask = srcfreqmask
    obj! = construct_pde_misfit(vel_model,Q,D,model,opts,batch_mode=false)
    function objective_function(in,obj!)
      f=0.0
      g=zeros(length(in))
      f=obj!(in,g);
      return f,g
    end
    f_obj = in -> objective_function(in,obj!)
    #vel_model = spg(obj!,proj!,vel_model,3,maxfc=max_func_evals)
    vel_model, fsave, funEvals= minConf_SPG(f_obj, vec(vel_model), proj!, options_SPG)
end

#save new updated model
save("vel_model.jld","vel_model",vel_model)

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
#=
AbstractModel
=============
    An abstract financial optimization model, all concrete models in this package
    (including the customizable models) are children of this type.

    Methods:
    -------
    get_obj_fun(m) - retrieves the objective function of the model

Author: Fady Shoukry, Shen Wang
Date: 01/23/2016
=#

abstract AbstractModel

function getObjective(m::AbstractModel)
    """
    Return the objective function of the given Model as an Expr type
    """
    m._objective
end

function getDefaultConstraints(m::AbstractModel)
    """
    Return the default constraints of a given Model
    type
    """
    m._default_constraints
end

function getConstraints(m::AbstractModel)
    """
    Return an array of constraints as expressions
    """
    user_constraints = collect(values(m.constraints))
    all_constraints = vcat(m._default_constraints, user_constraints)
    all_constraints
end

function getSense(m::AbstractModel)
    """
    Return the Sense of the Model, Min or Max
    """
    m.sense
end

function getVariables(m::AbstractModel)
    """
    Return the list of variables in the model
    """
    m.vars
end

function optimize(m::AbstractModel, syms_dict=Dict{Symbol,Any}()::Dict; solver=JuMP.UnsetSolver())
    """
    Return the optimized weights of the model generated as an array of floats
    """
    constraints = getConstraints(m)

    vars = getVariables(m)
    sense = getSense(m)
    objective = getObjective(m)

    # Next, generate the function that will create the JuMP model
    modelGen = createJuMPModelGenFunc(vars, sense, objective, constraints, syms_dict)
    JuMPModel = modelGen(solver)
    # Once the constraints are added, and the objective is set, solve the model
    TT = STDOUT
    #redirect_stdout()
    m.status = JuMP.solve(JuMPModel)
    redirect_stdout(TT)

    # Report the status and the solution
    m.objVal = JuMP.getobjectivevalue(JuMPModel)
    m.weights = JuMP.getvalue(JuMP.getvariable(JuMPModel, :w))

    # return the weights

    return m.objVal, m.weights
end

function exportModelResultsToCSV{R<:Real}(result::Tuple{R,Array{R,1}}, path::AbstractString)

    df = DataFrames.DataFrame()
    df[:ObjValue_And_Weights] = vcat(result[1], result[2])

    DataFrames.writetable(path, df)
    return true
end

function Base.show(io::IO, m::AbstractModel)
    # print(io, "\n Sense: $(m.sense) \n")
    print(io, "\n Variables: \n")
    for vars in m.vars
        print(io, vars, "\n")
    end

    # print(io, "\n Objective Function: \n  $(m._objective) \n")

    print(io, "\n Constraints: \n")

    constraint_key_array = collect(keys(m.constraints))
    constraint_value_array = collect(values(m.constraints))
    constraint_length = length(constraint_key_array)
    constraint_df = DataFrames.DataFrame(Keys = Symbol[], Constraint = Expr[])

    def_constraint_df = DataFrames.DataFrame(Default = AbstractString[], Constraint = Expr[])

    default_constraints_value = m._default_constraints
    default_constraints_length = length(default_constraints_value)

    for i = 1:default_constraints_length
        push!(def_constraint_df, ["default", default_constraints_value[i]])
    end

    for i = 1:constraint_length
        push!(constraint_df, [constraint_key_array[i], constraint_value_array[i]])
    end

    print(io, constraint_df)

    print(io, "\n\n")

    # print(io, def_constraint_df)

    # print(io, "\n\n")

    print(io, "\n Assets: \n $(m.assets) \n")

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

const ByteArray = Union{Vector{UInt8},Vector{Int8}}

## constructors and conversions ##

# String constructor docstring from boot.jl, workaround for #16730
# and the unavailability of @doc in boot.jl context.
"""
    String(v::Vector{UInt8})

Create a new `String` from a vector `v` of bytes containing
UTF-8 encoded characters.   This function takes "ownership" of
the array, which means that you should not subsequently modify
`v` (since strings are supposed to be immutable in Julia) for
as long as the string exists.

If you need to subsequently modify `v`, use `String(copy(v))` instead.
"""
function String(v::Array{UInt8,1})
    ccall(:jl_array_to_string, Ref{String}, (Any,), v)
end

"""
    unsafe_string(p::Ptr{UInt8}, [length::Integer])

Copy a string from the address of a C-style (NUL-terminated) string encoded as UTF-8.
(The pointer can be safely freed afterwards.) If `length` is specified
(the length of the data in bytes), the string does not have to be NUL-terminated.

This function is labelled "unsafe" because it will crash if `p` is not
a valid memory address to data of the requested length.
"""
function unsafe_string(p::Union{Ptr{UInt8},Ptr{Int8}}, len::Integer)
    p == C_NULL && throw(ArgumentError("cannot convert NULL to string"))
    ccall(:jl_pchar_to_string, Ref{String}, (Ptr{UInt8}, Int), p, len)
end
function unsafe_string(p::Union{Ptr{UInt8},Ptr{Int8}})
    p == C_NULL && throw(ArgumentError("cannot convert NULL to string"))
    ccall(:jl_cstr_to_string, Ref{String}, (Ptr{UInt8},), p)
end

_string_n(n::Integer) = ccall(:jl_alloc_string, Ref{String}, (Csize_t,), n)

"""
    String(s::AbstractString)

Convert a string to a contiguous byte array representation encoded as UTF-8 bytes.
This representation is often appropriate for passing strings to C.
"""
String(s::AbstractString) = print_to_string(s)

String(s::Symbol) = unsafe_string(Cstring(s))

(::Type{Vector{UInt8}})(s::String) = ccall(:jl_string_to_array, Ref{Vector{UInt8}}, (Any,), s)

## low-level functions ##

pointer(s::String) = unsafe_convert(Ptr{UInt8}, s)
pointer(s::String, i::Integer) = pointer(s)+(i-1)

sizeof(s::String) = Core.sizeof(s)

"""
    codeunit(s::AbstractString, i::Integer)

Get the `i`th code unit of an encoded string. For example,
returns the `i`th byte of the representation of a UTF-8 string.

# Examples
```jldoctest
julia> s = "δ=γ"; [codeunit(s, i) for i in 1:sizeof(s)]
5-element Array{UInt8,1}:
 0xce
 0xb4
 0x3d
 0xce
 0xb3
```
"""
codeunit(s::AbstractString, i::Integer)

@inline function codeunit(s::String, i::Integer)
    @boundscheck if (i < 1) | (i > sizeof(s))
        throw(BoundsError(s,i))
    end
    @gc_preserve s unsafe_load(pointer(s, i))
end

write(io::IO, s::String) =
    @gc_preserve s unsafe_write(io, pointer(s), reinterpret(UInt, sizeof(s)))

## comparison ##

function cmp(a::String, b::String)
    al, bl = sizeof(a), sizeof(b)
    c = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
              a, b, min(al,bl))
    return c < 0 ? -1 : c > 0 ? +1 : cmp(al,bl)
end

function ==(a::String, b::String)
    al = sizeof(a)
    al == sizeof(b) && 0 == ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt), a, b, al)
end

## prevind and nextind ##

function prevind(s::String, i::Integer)
    j = Int(i)
    e = sizeof(s)
    if j > e
        return endof(s)
    end
    j -= 1
    @inbounds while j > 0 && is_valid_continuation(codeunit(s,j))
        j -= 1
    end
    j
end

function prevind(s::String, i::Integer, nchar::Integer)
    nchar > 0 || throw(ArgumentError("nchar must be greater than 0"))
    j = Int(i)
    e = sizeof(s)
    while nchar > 0
        if j > e
            j = endof(s)
        else
            j -= 1
            @inbounds while j > 0 && is_valid_continuation(codeunit(s,j))
                j -= 1
            end
        end
        nchar -= 1
        j <= 0 && return j - nchar
    end
    j
end

function nextind(s::String, i::Integer)
    j = Int(i)
    if j < 1
        return 1
    end
    e = sizeof(s)
    j += 1
    @inbounds while j <= e && is_valid_continuation(codeunit(s,j))
        j += 1
    end
    j
end

function nextind(s::String, i::Integer, nchar::Integer)
    nchar > 0 || throw(ArgumentError("nchar must be greater than 0"))
    j = Int(i)
    e = sizeof(s)
    while nchar > 0
        if j < 1
            j = 1
        else
            j += 1
            @inbounds while j <= e && is_valid_continuation(codeunit(s,j))
                j += 1
            end
        end
        nchar -= 1
        j > e && return j + nchar
    end
    j
end

## checking UTF-8 & ACSII validity ##

byte_string_classify(data::Vector{UInt8}) =
    ccall(:u8_isvalid, Int32, (Ptr{UInt8}, Int), data, length(data))
byte_string_classify(s::String) =
    ccall(:u8_isvalid, Int32, (Ptr{UInt8}, Int), s, sizeof(s))
    # 0: neither valid ASCII nor UTF-8
    # 1: valid ASCII
    # 2: valid UTF-8

isvalid(::Type{String}, s::Union{Vector{UInt8},String}) = byte_string_classify(s) != 0
isvalid(s::String) = isvalid(String, s)

## basic UTF-8 decoding & iteration ##

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

const utf8_offset = [
    0x00000000, 0x00003080,
    0x000e2080, 0x03c82080,
    0xfa082080, 0x82082080,
]

const utf8_trailing = [
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3,4,4,4,4,5,5,5,5,
]

## required core functionality ##

function endof(s::String)
    i = sizeof(s)
    @inbounds while i > 0 && is_valid_continuation(codeunit(s, i))
        i -= 1
    end
    i
end

function length(s::String)
    cnum = 0
    @inbounds for i = 1:sizeof(s)
        cnum += !is_valid_continuation(codeunit(s, i))
    end
    cnum
end

@noinline function slow_utf8_next(s::String, b::UInt8, i::Int, l::Int)
    @inbounds if is_valid_continuation(b)
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, codeunit(s, i)))
    end
    trailing = utf8_trailing[b + 1]
    if l < i + trailing
        return '\ufffd', i+1
    end
    c::UInt32 = 0
    @inbounds for j = 1:(trailing + 1)
        c <<= 6
        c += codeunit(s, i)
        i += 1
    end
    c -= utf8_offset[trailing + 1]
    return Char(c), i
end

# This implementation relies on `next` returning a value past the end of the
# String's underlying data, which is true for valid Strings
done(s::String, state) = state > sizeof(s)

@inline function next(s::String, i::Int)
    # function is split into this critical fast-path
    # for pure ascii data, such as parsing numbers,
    # and a longer function that can handle any utf8 data
    @boundscheck if (i < 1) | (i > sizeof(s))
        throw(BoundsError(s,i))
    end
    @inbounds b = codeunit(s, i)
    if b < 0x80
        return Char(b), i + 1
    end
    return slow_utf8_next(s, b, i, sizeof(s))
end

function first_utf8_byte(ch::Char)
    c = UInt32(ch)
    b = c < 0x80    ? c%UInt8 :
        c < 0x800   ? ((c>>6)  | 0xc0)%UInt8 :
        c < 0x10000 ? ((c>>12) | 0xe0)%UInt8 :
                      ((c>>18) | 0xf0)%UInt8
    return b
end

function reverseind(s::String, i::Integer)
    j = sizeof(s) + 1 - i
    @inbounds while is_valid_continuation(codeunit(s, j))
        j -= 1
    end
    return j
end

## overload methods for efficiency ##

isvalid(s::String, i::Integer) =
    (1 <= i <= sizeof(s)) && ((@inbounds b = codeunit(s, i)); !is_valid_continuation(b))

function getindex(s::String, r::UnitRange{Int})
    isempty(r) && return ""
    l = sizeof(s)
    i = first(r)
    if i < 1 || i > l
        throw(BoundsError(s, i))
    end
    @inbounds si = codeunit(s, i)
    if is_valid_continuation(si)
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, si))
    end
    j = last(r)
    if j > l
        throw(BoundsError(s, j))
    end
    @inbounds sj = codeunit(s, j)
    if is_valid_continuation(sj)
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, j, sj))
    end
    j = nextind(s,j)
    unsafe_string(pointer(s,i), j-i)
end

function search(s::String, c::Char, i::Integer = 1)
    if i < 1 || i > sizeof(s)
        i == sizeof(s) + 1 && return 0
        throw(BoundsError(s, i))
    end
    @inbounds if is_valid_continuation(codeunit(s,i))
        throw(UnicodeError(UTF_ERR_INVALID_INDEX, i, codeunit(s,i)))
    end
    c < Char(0x80) && return search(s, c%UInt8, i)
    while true
        i = search(s, first_utf8_byte(c), i)
        (i==0 || s[i] == c) && return i
        i = next(s,i)[2]
    end
end

function search(a::Union{String,ByteArray}, b::Union{Int8,UInt8}, i::Integer = 1)
    if i < 1
        throw(BoundsError(a, i))
    end
    n = sizeof(a)
    if i > n
        return i == n+1 ? 0 : throw(BoundsError(a, i))
    end
    p = pointer(a)
    q = ccall(:memchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), p+i-1, b, n-i+1)
    q == C_NULL ? 0 : Int(q-p+1)
end

function search(a::ByteArray, b::Char, i::Integer = 1)
    if isascii(b)
        search(a,UInt8(b),i)
    else
        search(a,Vector{UInt8}(string(b)),i).start
    end
end

function rsearch(s::String, c::Char, i::Integer = sizeof(s))
    c < Char(0x80) && return rsearch(s, c%UInt8, i)
    b = first_utf8_byte(c)
    while true
        i = rsearch(s, b, i)
        (i==0 || s[i] == c) && return i
        i = prevind(s,i)
    end
end

function rsearch(a::Union{String,ByteArray}, b::Union{Int8,UInt8}, i::Integer = sizeof(s))
    if i < 1
        return i == 0 ? 0 : throw(BoundsError(a, i))
    end
    n = sizeof(a)
    if i > n
        return i == n+1 ? 0 : throw(BoundsError(a, i))
    end
    p = pointer(a)
    q = ccall(:memrchr, Ptr{UInt8}, (Ptr{UInt8}, Int32, Csize_t), p, b, i)
    q == C_NULL ? 0 : Int(q-p+1)
end

function rsearch(a::ByteArray, b::Char, i::Integer = length(a))
    if isascii(b)
        rsearch(a,UInt8(b),i)
    else
        rsearch(a,Vector{UInt8}(string(b)),i).start
    end
end

## optimized concatenation, reverse, repeat ##

function string(a::String...)
    if length(a) == 1
        return a[1]::String
    end
    n = 0
    for str in a
        n += sizeof(str)
    end
    out = _string_n(n)
    offs = 1
    for str in a
        unsafe_copy!(pointer(out,offs), pointer(str), sizeof(str))
        offs += sizeof(str)
    end
    return out
end

# UTF-8 encoding length of a character
function codelen(d::Char)
    c = UInt32(d)
    if c < 0x80
        return 1
    elseif c < 0x800
        return 2
    elseif c < 0x10000
        return 3
    elseif c < 0x110000
        return 4
    end
    return 3  # '\ufffd'
end

function string(a::Union{String,Char}...)
    n = 0
    for d in a
        if isa(d,Char)
            n += codelen(d::Char)
        else
            n += sizeof(d::String)
        end
    end
    out = _string_n(n)
    offs = 1
    p = pointer(out)
    for d in a
        if isa(d,Char)
            c = UInt32(d::Char)
            if c < 0x80
                unsafe_store!(p, c%UInt8, offs); offs += 1
            elseif c < 0x800
                unsafe_store!(p, (( c >> 6          ) | 0xC0)%UInt8, offs); offs += 1
                unsafe_store!(p, (( c        & 0x3F ) | 0x80)%UInt8, offs); offs += 1
            elseif c < 0x10000
                unsafe_store!(p, (( c >> 12         ) | 0xE0)%UInt8, offs); offs += 1
                unsafe_store!(p, (((c >> 6)  & 0x3F ) | 0x80)%UInt8, offs); offs += 1
                unsafe_store!(p, (( c        & 0x3F ) | 0x80)%UInt8, offs); offs += 1
            elseif c < 0x110000
                unsafe_store!(p, (( c >> 18         ) | 0xF0)%UInt8, offs); offs += 1
                unsafe_store!(p, (((c >> 12) & 0x3F ) | 0x80)%UInt8, offs); offs += 1
                unsafe_store!(p, (((c >> 6)  & 0x3F ) | 0x80)%UInt8, offs); offs += 1
                unsafe_store!(p, (( c        & 0x3F ) | 0x80)%UInt8, offs); offs += 1
            else
                # '\ufffd'
                unsafe_store!(p, 0xef, offs); offs += 1
                unsafe_store!(p, 0xbf, offs); offs += 1
                unsafe_store!(p, 0xbd, offs); offs += 1
            end
        else
            l = sizeof(d::String)
            unsafe_copy!(pointer(out,offs), pointer(d::String), l)
            offs += l
        end
    end
    return out
end

function reverse(s::String)
    dat = Vector{UInt8}(s)
    n = length(dat)
    n <= 1 && return s
    buf = StringVector(n)
    out = n
    pos = 1
    @inbounds while out > 0
        ch = dat[pos]
        if ch > 0xdf
            if ch < 0xf0
                (out -= 3) < 0 && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
                buf[out + 1], buf[out + 2], buf[out + 3] = ch, dat[pos + 1], dat[pos + 2]
                pos += 3
            else
                (out -= 4) < 0 && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
                buf[out+1], buf[out+2], buf[out+3], buf[out+4] = ch, dat[pos+1], dat[pos+2], dat[pos+3]
                pos += 4
            end
        elseif ch > 0x7f
            (out -= 2) < 0 && throw(UnicodeError(UTF_ERR_SHORT, pos, ch))
            buf[out + 1], buf[out + 2] = ch, dat[pos + 1]
            pos += 2
        else
            buf[out] = ch
            out -= 1
            pos += 1
        end
    end
    String(buf)
end

function repeat(s::String, r::Integer)
    r < 0 && throw(ArgumentError("can't repeat a string $r times"))
    n = sizeof(s)
    out = _string_n(n*r)
    if n == 1 # common case: repeating a single ASCII char
        @inbounds ccall(:memset, Ptr{Void}, (Ptr{UInt8}, Cint, Csize_t), out, codeunit(s, 1), r)
    else
        for i=1:r
            unsafe_copy!(pointer(out, 1+(i-1)*n), pointer(s), n)
        end
    end
    return out
end

"""
    repeat(c::Char, r::Integer) -> String

Repeat a character `r` times. This can equivalently be accomplished by calling [`c^r`](@ref ^).

# Examples
```jldoctest
julia> repeat('A', 3)
"AAA"
```
"""
function repeat(c::Char, r::Integer)
    r < 0 && throw(ArgumentError("can't repeat a character $r times"))
    r == 0 && return ""
    ch = UInt(c)
    if ch < 0x80
        out = _string_n(r)
        ccall(:memset, Ptr{Void}, (Ptr{UInt8}, Cint, Csize_t), out, c, r)
    elseif ch < 0x800
        out = _string_n(2r)
        p16 = reinterpret(Ptr{UInt16}, pointer(out))
        u16 = ((ch >> 0x6) | (ch & 0x3f) << 0x8) % UInt16 | 0x80c0
        @inbounds for i = 1:r
            unsafe_store!(p16, u16, i)
        end
    elseif ch < 0x10000
        (0xd800 ≥ ch ≤ 0xdfff) || throw(ArgumentError("invalid character 0x$(hex(ch))"))
        out = _string_n(3r)
        p = pointer(out)
        b1 = (ch >> 0xc) % UInt8 | 0xe0
        b2 = ((ch >> 0x6) & 0x3f) % UInt8 | 0x80
        b3 = (ch & 0x3f) % UInt8 | 0x80
        @inbounds for i = 1:r
            unsafe_store!(p, b1)
            unsafe_store!(p, b2, 2)
            unsafe_store!(p, b3, 3)
            p += 3
        end
    elseif ch < 0x110000
        out = _string_n(4r)
        p32 = reinterpret(Ptr{UInt32}, pointer(out))
        u32 = ((ch >> 0x12) | ((ch >> 0x4) & 0x03f00) |
            ((ch << 0xa) & 0x3f0000) | ((ch & 0x3f) << 0x18)) % UInt32 | 0x808080f0
        @inbounds for i = 1:r
            unsafe_store!(p32, u32)
            p32 += 4
        end
    else
        throw(ArgumentError("invalid character 0x$(hex(ch))"))
    end
    return out
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
########################################
# Main structures
########################################
struct Locale{T}
    index::Vector{Int}
    locs::T
end

struct BlockRow{T}
    v4net::T
    geoname_id::Int
    location::Union{LLA{Float64}, Missing}
    registered_country_geoname_id::Union{Int, Missing}
    is_anonymous_proxy::Int
    is_satellite_provider::Int
    postal_code::Union{String, Missing}
    accuracy_radius::Union{Int, Missing}
end

function BlockRow(csvrow)
    net = IPNets.IPv4Net(csvrow.network)
    geoname_id = ismissing(csvrow.geoname_id) ? -1 : csvrow.geoname_id

    lat = csvrow.latitude
    lon = csvrow.longitude
    location = ismissing(lon) || ismissing(lat) ? missing : LLA(lat, lon)
    registered_country_geoname_id = csvrow.registered_country_geoname_id
    accuracy_radius = get(csvrow, :accuracy_radius, missing)
    postal_code = csvrow.postal_code

    BlockRow(
        net,
        geoname_id,
        location,
        registered_country_geoname_id,
        csvrow.is_anonymous_proxy,
        csvrow.is_satellite_provider,
        postal_code,
        accuracy_radius
    )
end

struct DB{T1, T2 <: Locale}
    index::Vector{T1}
    blocks::Vector{BlockRow{T1}}
    locs::Vector{T2}
    localeid::Int
    ldict::Dict{Symbol, Int}
end

Base.broadcastable(db::DB) = Ref(db)
"""
    setlocale(db, localename)

Set new locale which should be used in return results. If locale is not found, then current locale is going to be used.
"""
function setlocale(db::DB, localename)
    if localename in keys(db.ldict)
        return DB(db.index, db.blocks, db.locs, db.ldict[localename], db.ldict)
    else
        @warn "Unable to find locale $localename"
        return db
    end
end

# Path to directory with data, can define GEOIP_DATADIR to override
# the default (useful for testing with a smaller test set)
function getdatadir(datadir)
    isempty(datadir) || return datadir
    haskey(ENV, "GEOIP_DATADIR") ?  ENV["GEOIP_DATADIR"] : datadir
end

function getzipfile(zipfile)
    isempty(zipfile) || return zipfile
    haskey(ENV, "GEOIP_ZIPFILE") ? ENV["GEOIP_ZIPFILE"] : zipfile
end

getlocale(x::Pair) = x
function getlocale(x::Symbol)
    if x == :en
        return :en => r"Locations-en.csv$"
    elseif x == :de
        return :de => r"Locations-de.csv$"
    elseif x == :ru
        return :ru => r"Locations-ru.csv$"
    elseif x == :ja
        return :ja => r"Locations-ja.csv$"
    elseif x == :es
        return :es => r"Locations-es.csv$"
    elseif x == :fr
        return :fr => r"Locations-fr.csv$"
    elseif x == :pt_br
        return :pt_br => r"Locations-pt-BR.csv$"
    elseif x == :zh_cn
        return :zh_cn => r"Locations-zh_cn.csv$"
    end
end

function loadzip(datadir, zipfile, locales)
    zipfile = joinpath(datadir, zipfile)
    isfile(zipfile) || throw(ArgumentError("Unable to find data file in $(zipfile)"))

    r = ZipFile.Reader(zipfile)
    ldict = Dict{Symbol, Int}()
    locid = 1
    local blocks
    locs = []
    try
        for f in r.files
            for (l, s) in locales
                if occursin(s, f.name)
                    v = Vector{UInt8}(undef, f.uncompressedsize)
                    ls = read!(f, v) |> CSV.File
                    push!(locs, ls)
                    ldict[l] = locid
                    locid += 1
                end
            end
            if occursin(r"Blocks-IPv4.csv$", f.name)
                v = Vector{UInt8}(undef, f.uncompressedsize)
                blocks = read!(f, v) |> x -> CSV.File(x; types = Dict(:postal_code => String))
            end
        end
    catch
        @error "Geolocation data cannot be read. Data directory may be corrupt..."
        rethrow()
    finally
        close(r)
    end

    return blocks, locs, ldict
end

"""
    load(; datadir, zipfile, locales, deflocale)

Load GeoIP database from compressed CSV file. The argument `zipfile` should be provided, otherwise `load` function error. `datadir` defines where data files are located and can be either set as an argument or read from the `ENV` variable `GEOIP_DATADIR`. In the same way if `ENV` variable `GEOIP_ZIPFILE` is set, then it is used for determining `zipfile` argument.

Argument `locales` determine locale files which should be loaded. Locales can be given as `Symbol`s or `Pair`s of locale name and filename which contains corresponding locale, e.g. `locales = [:en, :fr]` or `locales = [:en => r"-en.csv"]`. Following locales are supported in `Symbol` version `:en, :de, :ru, :ja, :es, :fr, :pt_br, :zh_cn`. To set default locale use `deflocale` argument, e.g. `deflocale = :en`.
"""
function load(; zipfile = "",
                datadir = "",
                locales = [:en],
                deflocale = :en)
    datadir = getdatadir(datadir)
    zipfile = getzipfile(zipfile)
    locales = getlocale.(locales)
    blocks, locs, ldict = loadzip(datadir, zipfile, locales)

    blockdb = BlockRow.(blocks)
    sort!(blockdb, by = x -> x.v4net)
    index = map(x -> x.v4net, blockdb)
    locsdb = map(locs) do loc
        ldb = collect(loc)
        sort!(ldb, by = x -> x.geoname_id)
        lindex = map(x -> x.geoname_id, ldb)
        Locale(lindex, ldb)
    end

    localeid = if deflocale in keys(ldict)
        ldict[deflocale]
    else
        cd = collect(d)
        idx = findfirst(x -> x[2] == 1, cd)
        locname = cd[idx][1]
        @warn "Default locale $deflocale was not found, using locale $locname"
        1
    end

    return DB(index, blockdb, locsdb, localeid, ldict)
end


text:
using Test
using HTTP
using HTTP.IOExtras, HTTP.Sockets
using Sockets

@testset "websockets.jl" begin
    p = 8085 # rand(8000:8999)
    socket_type = ["wss", "ws"]

    function listen_localhost()
        @async HTTP.listen(Sockets.localhost, p) do http
            if HTTP.WebSockets.is_upgrade(http.message)
                HTTP.WebSockets.upgrade(http) do ws
                    while !eof(ws)
                        data = readavailable(ws)
                        write(ws, data)
                    end
                end
            end
        end
    end

    @testset "External Host - $s" for s in socket_type
        HTTP.WebSockets.open("$s://echo.websocket.org") do io
            write(io, "Foo")
            @test !eof(io)
            @test String(readavailable(io)) == "Foo"

            write(io, "Hello")
            write(io, " There")
            write(io, " World", "!")
            closewrite(io)

            buf = IOBuffer()
            write(buf, io)
            @test String(take!(buf)) == "Hello There World!"
        end
    end

    @testset "Localhost" begin
       listen_localhost()

        HTTP.WebSockets.open("ws://127.0.0.1:$(p)") do ws
            write(ws, "Foo")
            @test String(readavailable(ws)) == "Foo"

            write(ws, "Bar")
            @test String(readavailable(ws)) == "Bar"
        end
    end

    @testset "extened feautre support for listen" begin
        port=UInt16(8086)
        tcpserver = listen(port)
        target = "/query?k1=v1&k2=v2"

        servertask =  @async HTTP.WebSockets.listen("127.0.0.1", port; server=tcpserver) do ws
            @testset "request access" begin
                @test ws.request isa HTTP.Request
                write(ws, ws.request.target)
                while !eof(ws)
                    write(ws, readavailable(ws))
                end
                close(ws)
            end
        end

        HTTP.WebSockets.open("ws://127.0.0.1:$(port)$(target)") do ws
            @test String(readavailable(ws)) == target
            @test write(ws, "Bye!") == 4
            @test String(readavailable(ws)) == "Bye!"
            close(ws)
        end

        close(tcpserver)
        @test timedwait(()->servertask.state === :failed, 5.0) === :ok
        @test_throws Exception wait(servertask)
    end
end


text:
using Documenter, DoubleExponentialFormulas

makedocs(;
    modules=[DoubleExponentialFormulas],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Reference" => "reference.md",
    ],
    repo="https://github.com/machakann/DoubleExponentialFormulas.jl/blob/{commit}{path}#L{line}",
    sitename="DoubleExponentialFormulas.jl",
    authors="Masaaki Nakamura",
    assets=String[],
)

deploydocs(;
    repo="github.com/machakann/DoubleExponentialFormulas.jl",
    deploy_config=Documenter.GitHubActions(),
)


text:
"""
	gmt2kml(cmd0::String="", arg1=[], kwargs...)

Convert GMT data tables to KML files for Google Earth

Full option list at [`gmt2kml`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html)

Parameters
----------

- **A** : **altitude_mode** : -- Str --       Flags = a|g|s[alt|xscale]

    Select one of three altitude modes recognized by Google Earth that determines the altitude (in m)
    of the feature: ``a`` absolute altitude, ``g`` altitude relative to sea surface or ground,
    ``s`` altitude relative to seafloor or ground.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#A)
- $(GMT.opt_C)
- **D** : **descript** : -- Str --   Flags = descriptfile

    File with HTML snippets that will be included as part of the main description content for the KML file.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#d)
- **E** : **extrude** : -- Str or [] --   Flags = [altitude]

    Extrude feature down to ground level.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#e)
- **F** : **feature_type** : -- Str --   Flags = e|s|t|l|p|w

    Sets the feature type. Choose from points (event, symbol, or timespan), line, polygon, or wiggle.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#f)
- **G** : **fill_color** : -- Str --  Flags = f|nfill

    Sets color fill (G=:f) or label font color (G=:n).
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#g)
- **I** : **icon** : -- Str --      Flags = icon

    Specify the URL to an alternative icon that should be used for the symbol
    [Default is a Google Earth circle].
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#i)
- **K** : **not_finished** : -- Bool or [] --

    Allow more KML code to be appended to the output later [finalize the KML file].
    [`-K`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#i)
- **L** : **extended_data** : -- Str --      Flags = name1,name2,…

    Extended data given. Append one or more column names separated by commas.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#l)
- **N** : **feature_name** : -- Str or Number --      Flags = [t|col |name_template|name]

    By default, if segment headers contain a -L”label string” then we use that for the name of the KML feature.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#n)
- **O** : **overlay** : -- Bool or [] --

    Append KML code to an existing KML file [initialize a new KML file].
    [`-O`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#n)
- **Qa** : **wiggles** : -- Str --      Flags =  azimuth

    Option in support of wiggle plots (requires F=:w).
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#q)
- **Qs** : **wiggle_scale** : -- Number or Str --      Flags =  scale[unit]

    Required setting for wiggle plots (i.e., it requires F=:w). Sets a wiggle scale
    in z-data units per the user’s units
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#q)
- $(GMT.opt_R)
- **S** : **scale** : -- Str --      Flags =  c|nscale

    Scale icons or labels. Here, S=:c sets a scale for the symbol icon, whereas S=:n sets
    a scale for the name labels
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#s)
- **T** : **title ** : -- Str --    Flags = title[/foldername]

    Sets the document title [default is unset]. Optionally, append /FolderName;
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#t)
- **W** : **pen** : -- Str or [] --      Flags =  [pen][attr]

    Set pen attributes for lines, wiggles or polygon outlines.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#w)
- **Z** : **attrib** : -- Str --      Flags =  args

    Set one or more attributes of the Document and Region tags.
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/gmt2kml.html#z)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_swap_xy)
"""
function gmt2kml(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("gmt2kml", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd, = parse_R("", d)
	cmd  = parse_V_params(cmd, d)
	cmd, = parse_bi(cmd, d)
	cmd, = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'A', d, [:A :altitude_mode])
	cmd = add_opt(cmd, 'D', d, [:D :descript])
	cmd = add_opt(cmd, 'E', d, [:E :extrude])
	cmd = add_opt(cmd, 'F', d, [:F :feature_type])
	cmd = add_opt(cmd, 'G', d, [:G :fill_color])
	cmd = add_opt(cmd, 'I', d, [:I :icon])
	cmd = add_opt(cmd, 'K', d, [:K :not_finished])
	cmd = add_opt(cmd, 'L', d, [:L :extended_data])
	cmd = add_opt(cmd, 'N', d, [:N :feature_name])
	cmd = add_opt(cmd, 'O', d, [:O :overlay])
	cmd = add_opt(cmd, "Qa", d, [:Qa :wiggles])
	cmd = add_opt(cmd, "Qs", d, [:Qs :wiggle_scale])
	cmd = add_opt(cmd, 'S', d, [:S :scale])
	cmd = add_opt(cmd, 'T', d, [:T :title])
	cmd = add_opt(cmd, 'Z', d, [:Z :attrib])

	for symb in [:W :pen]
		if (haskey(d, symb))
			if (isa(d[symb], Tuple))
				cmd = cmd * parse_pen(d[symb])
			else
				cmd = cmd * arg2str(d[symb]);
				lc = parse_pen_color(d)
				if (!isempty(lc))
					cmd = cmd * "," * lc
					ls = parse_pen_style(d)
					if (!isempty(ls))		cmd = cmd * "," * ls	end
				end
			end
			break
		end
	end

    cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)

	N_used = got_fname == 0 ? 1 : 0				# To know whether a cpt will go to arg1 or arg2
	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_used, arg1, arg2)

	return common_grd(d, cmd, got_fname, 1, "gmt2kml", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmt2kml(arg1=[], cmd0::String=""; kw...) = gmt2kml(cmd0, arg1; kw...)