text:
export conv4, deconv4, mat, pool, unpool
using NNlib: conv, DenseConvDims, maxpool, meanpool, PoolDims, ∇conv_data, ∇conv_filter, ∇maxpool, ∇meanpool
using LinearAlgebra: lmul!
using AutoGrad: AutoGrad, @primitive1

"""
    conv4(w, x; kwargs...)

Execute convolutions or cross-correlations using filters specified with `w` over tensor `x`.

If `w` has dimensions `(W1,W2,...,Cx,Cy)` and `x` has dimensions `(X1,X2,...,Cx,N)`, the
result `y` will have dimensions `(Y1,Y2,...,Cy,N)` where `Cx` is the number of input channels,
`Cy` is the number of output channels, `N` is the number of instances, and `Wi,Xi,Yi` are
spatial dimensions with `Yi` determined by:

    Yi = 1 + floor((Xi + 2*padding[i] - ((Wi-1)*dilation[i] + 1)) / stride[i])

`padding`, `stride` and `dilation` are keyword arguments that can be specified as a single
number (in which case they apply to all dimensions), or an array/tuple with entries for each
spatial dimension.

# Keywords

* `padding=0`: the number of extra zeros implicitly concatenated at the start and end of each dimension.
* `stride=1`: the number of elements to slide to reach the next filtering window.
* `dilation=1`: dilation factor for each dimension.
* `mode=0`: 0 for convolution and 1 for cross-correlation (which flips the filter).
* `alpha=1`: can be used to scale the result.
* `group=1`: can be used to perform grouped convolutions.

"""
function conv4(w, x; padding=0, stride=1, dilation=1, mode=0, alpha=1, group=1)
    @assert group == 1 "Grouped convolutions not yet implemented in NNlib, see https://github.com/JuliaGPU/CuArrays.jl/pull/523"
    N = ndims(x)
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    dilation = expand(Val(N-2), dilation)
    cdims = DenseConvDims(x, w; stride = stride, padding = padding, dilation = dilation, flipkernel = (mode!=0))
    y = conv(x, w, cdims)
    alpha == 1 ? y : lmul!(alpha, y)
end

function conv4w(w,x,dy; padding=0, stride=1, dilation=1, mode=0, alpha=1, group=1)
    @assert group == 1 "Grouped convolutions not yet implemented in NNlib, see https://github.com/JuliaGPU/CuArrays.jl/pull/523"
    N = ndims(x)
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    dilation = expand(Val(N-2), dilation)
    cdims = DenseConvDims(x, w; stride = stride, padding = padding, dilation = dilation, flipkernel = (mode!=0))
    dw = ∇conv_filter(x, dy, cdims)
    alpha == 1 ? dw : lmul!(alpha, dw)
end

function conv4x(w,x,dy; padding=0, stride=1, dilation=1, mode=0, alpha=1, group=1)
    @assert group == 1 "Grouped convolutions not yet implemented in NNlib, see https://github.com/JuliaGPU/CuArrays.jl/pull/523"
    N = ndims(x)
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    dilation = expand(Val(N-2), dilation)
    cdims = DenseConvDims(x, w; stride = stride, padding = padding, dilation = dilation, flipkernel = (mode!=0))
    dx = ∇conv_data(dy, w, cdims)
    alpha == 1 ? dx : lmul!(alpha, dx)
end

@primitive1 conv4(w,x; o...),dy,y       conv4w(w,x,dy;o...)   conv4x(w,x,dy;o...)
@primitive1 conv4w(w,x,dy;o...),ddw,dw  nothing               conv4x(ddw,x,dy;o...)  conv4(ddw,x;o...)
@primitive1 conv4x(w,x,dy;o...),ddx,dx  conv4w(w,ddx,dy;o...) nothing                conv4(w,ddx;o...)

"""
    deconv4(w, x; kwargs...)

Simulate 4-D deconvolution by using _transposed convolution_ operation. Its forward pass is
equivalent to backward pass of a convolution (gradients with respect to input
tensor). Likewise, its backward pass (gradients with respect to input tensor) is equivalent to
forward pass of a convolution. Since it swaps forward and backward passes of convolution
operation, padding and stride options belong to output tensor. See [this
report](https://arxiv.org/abs/1603.07285) for further explanation.

If `w` has dimensions `(W1,W2,...,Cy,Cx)` and `x` has dimensions `(X1,X2,...,Cx,N)`, the
result `y=deconv4(w,x)` will have dimensions `(Y1,Y2,...,Cy,N)` where

    Yi = (Xi - 1)*stride[i] + ((Wi-1)*dilation[i] + 1) - 2*padding[i]

Here Cx is the number of x channels, Cy is the number of y channels, N is the number of
instances, and Wi,Xi,Yi are spatial dimensions. Padding and stride are keyword arguments that
can be specified as a single number (in which case they apply to all dimensions), or an
array/tuple with entries for each spatial dimension.

# Keywords

* `padding=0`: the number of extra zeros implicitly concatenated at the start and at the end of each dimension.
* `stride=1`: the number of elements to slide to reach the next filtering window.
* `mode=0`: 0 for convolution and 1 for cross-correlation.
* `alpha=1`: can be used to scale the result.
* `handle`: handle to a previously created cuDNN context. Defaults to a Knet allocated handle.
* `group=1`: can be used to perform grouped convolutions.

"""
function deconv4(w,y; o...)
    x = similar(y,dcdims(w,y;o...))
    return conv4x(w,x,y;o...)
end

@primitive1 deconv4(w,x;o...),dy  conv4w(w,dy,x;o...)  conv4(w,dy;o...)

"""
    pool(x; kwargs...)

Compute pooling of input values (i.e., the maximum or average of several adjacent values) to
produce an output with smaller height and/or width.

If `x` has dimensions `(X1,X2,...,Cx,N)`, the result `y` will have dimensions
`(Y1,Y2,...,Cx,N)` where

    Yi=1+floor((Xi+2*padding[i]-window[i])/stride[i])

Here `Cx` is the number of input channels, `N` is the number of instances, and `Xi,Yi` are
spatial dimensions.  `window`, `padding` and `stride` are keyword arguments that can be
specified as a single number (in which case they apply to all dimensions), or an array/tuple
with entries for each spatial dimension.

# Keywords:

* `window=2`: the pooling window size for each dimension.
* `padding=0`: the number of extra zeros implicitly concatenated at the start and at the end of each dimension.
* `stride=window`: the number of elements to slide to reach the next pooling window.
* `mode=0`: 0 for max, 1 for average including padded values, 2 for average excluding padded values, 3 for deterministic max.
* `maxpoolingNanOpt=1`: Nan numbers are not propagated if 0, they are propagated if 1.
* `alpha=1`: can be used to scale the result.

"""
function pool(x; window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=1, alpha=1)
    mode, maxpoolingNanOpt = checkpoolopts(x, window, padding, stride, mode, maxpoolingNanOpt, alpha)
    N = ndims(x)
    window = expand(Val(N-2), window)
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    pdims = PoolDims(x, window; padding = padding, stride = stride)
    y = (mode == 0 ? maxpool(x, pdims) :
         mode == 1 ? meanpool(x, pdims) :
         mode == 2 ? error("Pool mode=2 not yet implemented in NNlib. See https://github.com/FluxML/NNlib.jl/issues/218") :
         mode == 3 ? maxpool(x, pdims) :
         error("mode=$mode is not supported for CPU pool."))
    alpha == 1 ? y : lmul!(alpha, y)
end

function poolx(x,y,dy; window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=1, alpha=1)
    mode, maxpoolingNanOpt = checkpoolopts(x, window, padding, stride, mode, maxpoolingNanOpt, alpha)
    if alpha != 1
        y = y ./ eltype(y)(alpha)
    end
    N = ndims(x)
    window = expand(Val(N-2), window)
    stride = expand(Val(N-2), stride)
    padding = expand(Val(N-2), padding)
    pdims = PoolDims(x, window; padding = padding, stride = stride)
    dx = (mode == 0 ? ∇maxpool(dy, y, x, pdims) :
          mode == 1 ? ∇meanpool(dy, y, x, pdims) :
          mode == 2 ? error("Pool mode=2 not yet implemented in NNlib. See https://github.com/FluxML/NNlib.jl/issues/218") :
          mode == 3 ? ∇maxpool(dy, y, x, pdims) :
          error("mode=$mode is not supported for CPU pool."))
    alpha == 1 ? dx : lmul!(alpha, dx)
end

@primitive1 pool(x;o...),dy,y  poolx(x,y,dy;o...)
@primitive1 poolx(x,y,dy;o...),ddx,dx  nothing  nothing  pool(ddx;o...)

function checkpoolopts(x, window, padding, stride, mode, maxpoolingNanOpt, alpha)
    @assert mode ∈ 0:3 "Bad pooling mode=$mode"
    if mode == 2
        @warn "Pool mode=2 not yet implemented in NNlib, using 1 instead. See https://github.com/FluxML/NNlib.jl/issues/218" maxlog=1
        mode = 1
    end
    @assert maxpoolingNanOpt ∈ (0,1) "Bad pooling maxpoolingNanOpt=$maxpoolingNanOpt"
    if maxpoolingNanOpt == 0
        @warn "Pool maxpoolingNanOpt=0 not yet implemented in NNlib, using 1 instead. See https://github.com/FluxML/NNlib.jl/issues/218" maxlog=1
        maxpoolingNanOpt = 1
    end
    if padding != 0 && x isa Array
        @warn "Pool padding is buggy in NNlib, use with caution. See https://github.com/FluxML/NNlib.jl/issues/229" maxlog=1
    end
    return (mode, maxpoolingNanOpt)
end

"""
    unpool(x; o...)

Perform the reverse of pooling: `x == pool(unpool(x;o...); o...)`
"""
function unpool(x; window=2, alpha=1, o...) # padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    w = prod(psize(window,x))
    y = similar(x,updims(x; window=window, o...))
    # Leave unpool as a non-primitive, it is just a poolx call
    poolx(y,x,x.*w; o..., window=window, mode=1, alpha=1/alpha)
end

"""
    mat(x; dims = ndims(x) - 1)

Reshape `x` into a two-dimensional matrix by joining the first dims dimensions, i.e.
`reshape(x, prod(size(x,i) for i in 1:dims), :)`

`dims=ndims(x)-1` (default) is typically used when turning the output of a 4-D convolution
result into a 2-D input for a fully connected layer.

`dims=1` is typically used when turning the 3-D output of an RNN layer into a 2-D input for
a fully connected layer.

`dims=0` will turn the input into a row vector, `dims=ndims(x)` will turn it into a column
vector.

"""
mat(x; dims::Int=ndims(x)-1)=reshape(x, (dims > 0 ? prod(size(x,i) for i in 1:dims) : 1), :)

## Dimension helpers:

# outputDim = 1 + ( inputDim + 2*pad - (((filterDim-1)*dilation)+1) )/convolutionStride;
# inputDim = (outputDim - 1) * convolutionStride + (((filterDim-1)*dilation)+1) - 2*pad
function dcdims(w,y; padding=0, stride=1, dilation=1, group=1, o...)
    N = ndims(y)
    @assert size(y,N-1) == size(w,N)
    ntuple(N) do i
        if i < N-1
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            di = (if isa(dilation,Number); dilation; else dilation[i]; end)
            si*(size(y,i)-1) + (((size(w,i)-1)*di)+1) - 2*pi
        elseif i == N-1
            size(w,N-1) * group
        else
            size(y,N)
        end
    end
end

# convert padding etc. size to an Int array of the right dimension
function psize(p, x)
    nd = ndims(x)-2
    if isa(p,Number)
        fill(Int(p),nd)
    elseif length(p)==nd
        collect(Int,p)
    else
        throw(DimensionMismatch("psize: $p $nd"))
    end
end

function updims(x; window=2, padding=0, stride=window, o...)
    window = psize(window,x)
    stride = psize(stride,x)
    padding = psize(padding,x)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            (size(x,i)-1)*stride[i]+window[i]-2*padding[i]
        else
            size(x,i)
        end
    end
end

# convolution padding size that preserves the input size when filter size is odd and stride=1
padsize(w)=ntuple(i->div(size(w,i)-1,2), ndims(w)-2)

expand(N, i::Tuple) = i
expand(N, i::Integer) = ntuple(_ -> i, N)

# TODO: Grouped convolutions not yet implemented in NNlib, see https://github.com/JuliaGPU/CuArrays.jl/pull/523
# TODO: maxpoolingNanOpt not yet implemented in NNlib, see https://github.com/FluxML/NNlib.jl/issues/218
# TODO: Pool mode=2 not yet implemented in NNlib, see https://github.com/FluxML/NNlib.jl/issues/218
# TODO: unpool Does not work correctly for every window, padding, mode combination?


text:
using Test
using StaticArrays
using Proj
import PROJ_jll

@testset "Error handling" begin
    @test Proj.proj_errno_string(0) === nothing
    @test Proj.proj_errno_string(1) == "Unknown error (code 1)"
    @test Proj.proj_errno_string(Proj.PROJ_ERR_INVALID_OP_WRONG_SYNTAX) ==
          "Invalid PROJ string syntax"

    pj = Proj.proj_create("EPSG:4326")
    Proj.proj_errno_set(pj, 14)
    @test Proj.proj_errno(pj) == 14
    Proj.proj_errno_reset(pj)
    @test Proj.proj_errno(pj) == 0
    pj = Proj.proj_destroy(pj)

    # this throws an internal error that is handled by our log_func
    @test_throws Proj.PROJError Proj.proj_create("+proj=bobbyjoe")
    # ensure we have reset the error state
    @test Proj.proj_context_errno() == 0

end

@testset "Database" begin
    crs = Proj.proj_create_from_database("EPSG", "4326", Proj.PJ_CATEGORY_CRS, false)
    Proj.proj_errno(crs)
    Proj.proj_get_id_code(crs, 0)
end

@testset "Create" begin
    pj_latlon = Proj.proj_create("EPSG:4326")

    esriwkt = """GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]"""
    pj2 = Proj.proj_create(esriwkt)

    @test startswith(Proj.proj_as_wkt(pj2, Proj.PJ_WKT1_ESRI, C_NULL), "GEOGCS")
    @test startswith(Proj.proj_as_wkt(pj2, Proj.PJ_WKT2_2018, C_NULL), "GEOGCRS")
end

@testset "Transformation between CRS" begin
    # taken from http://osgeo-org.1560.x6.nabble.com/PROJ-PROJ-6-0-0-proj-create-operation-factory-context-behavior-td5403470.html
    src = Proj.proj_create("IGNF:REUN47GAUSSL")   # area : 55.17,-21.42,55.92,-20.76
    tgt = Proj.proj_create("IGNF:RGAF09UTM20")    # area : -63.2,14.25,-60.73,18.2
    factory = Proj.proj_create_operation_factory_context(C_NULL)
    @test factory != C_NULL
    @test Proj.proj_context_errno() == 0
    results = Proj.proj_create_operations(src, tgt, factory)
    @test results != C_NULL
    n = Proj.proj_list_get_count(results)
    @test n == 1

    operation = Proj.proj_list_get(results, 0)
    @test Bool(Proj.proj_coordoperation_has_ballpark_transformation(operation))
    @test Bool(Proj.proj_coordoperation_has_ballpark_transformation(operation))
    @test startswith(Proj.proj_get_name(operation), "Inverse of GAUSS LABORDE REUNION")

    info = Proj.proj_pj_info(operation)
    @test unsafe_string(info.id) == "pipeline"
    @test startswith(unsafe_string(info.description), "Inverse of GAUSS LABORDE REUNION")
    @test startswith(unsafe_string(info.definition), "proj=pipeline step inv proj=gstmerc")
    @test Bool(info.has_inverse)
    @test info.accuracy === -1.0
    Proj.proj_destroy(operation)

    Proj.proj_list_destroy(results)

    Proj.proj_operation_factory_context_destroy(factory)
    Proj.proj_destroy(src)
    Proj.proj_destroy(tgt)
end

@testset "Quickstart" begin
    # https://proj.org/development/quickstart.html

    pj = Proj.proj_create_crs_to_crs(
        "EPSG:4326",  # source
        "+proj=utm +zone=32 +datum=WGS84",  # target, also EPSG:32632
    )

    # This will ensure that the order of coordinates for the input CRS
    # will be longitude, latitude, whereas EPSG:4326 mandates latitude,
    # longitude
    pj_for_gis = Proj.proj_normalize_for_visualization(pj)
    Proj.proj_destroy(pj)
    pj = pj_for_gis

    # a coordinate union representing Copenhagen: 55d N, 12d E
    # Given that we have used proj_normalize_for_visualization(), the order of
    # coordinates is longitude, latitude, and values are expressed in degrees.
    a = Proj.proj_coord(12, 55)
    @test a isa Proj.Coord
    @test a isa AbstractVector
    @test eltype(a) == Float64
    @test length(a) == 4
    @test sum(a) == Inf
    @test isbits(a)
    @test a[1] === 12.0
    @test a[2] === 55.0
    @test a[3] === 0.0
    @test a[4] === Inf

    # transform to UTM zone 32
    b = Proj.proj_trans(pj, Proj.PJ_FWD, a)
    @test is_approx(b, (691875.632, 6098907.825, 0.0, Inf))

    # inverse transform, back to geographical
    b = Proj.proj_trans(pj, Proj.PJ_INV, b)
    @test is_approx(b, (12.0, 55.0, 0.0, Inf))

    # Clean up
    pj = Proj.proj_destroy(pj)
    pj = Proj.proj_destroy(pj)
    # Julia crashes if proj_destroy is called twice on the same non nullptr,
    # since the contents at that address are already wiped. To prevent this,
    # change the input binding directly after, e.g. pj = proj_destroy(pj)  # -> C_NULL
end

@testset "inverse transformation" begin
    trans = Proj.Transformation(
        "EPSG:4326",
        "+proj=utm +zone=32 +datum=WGS84",
        always_xy = true,
    )

    # for custom / proj strings, or modified axis order, no description can be looked up in
    # the database
    @test repr(trans) == """
    Transformation
        source: unknown
        target: unknown
        direction: forward
    """

    trans⁻¹ = inv(trans, always_xy = true)
    trans¹ = inv(trans⁻¹, always_xy = true)

    # inv does not flip source and target, so the WKT stays the same
    wkt_type = Proj.PJ_WKT2_2019
    @test Proj.proj_as_wkt(trans⁻¹.pj, wkt_type) == Proj.proj_as_wkt(trans.pj, wkt_type)
    @test Proj.proj_as_wkt(trans¹.pj, wkt_type) == Proj.proj_as_wkt(trans.pj, wkt_type)

    a = (12.0, 55.0)
    b = (691875.632, 6098907.825)
    @test is_approx(trans⁻¹(b), a)
    @test is_approx(trans¹(a), b)
    @test trans⁻¹.direction == PJ_INV
    @test trans¹.direction == PJ_FWD

    # we don't check if source and target are equal when constructing a identity transform
    # since Transformation is mutable, it can always be changed after construction
    trans = Proj.Transformation(
        "EPSG:4326",
        "+proj=utm +zone=32 +datum=WGS84",
        direction = PJ_IDENT,
    )
    @test trans(a) === a
    trans⁻¹ = inv(trans, always_xy = true)
    @test trans⁻¹(a) === a
    @test trans⁻¹.direction == PJ_IDENT
    trans⁻¹.direction = PJ_FWD
    @test is_approx(trans⁻¹(a), b)

    @test inv(PJ_FWD) == PJ_INV
    @test inv(PJ_IDENT) == PJ_IDENT
    @test inv(PJ_INV) == PJ_FWD
end

@testset "single transformation" begin
    source_crs = Proj.proj_create("EPSG:4326")
    target_crs = Proj.proj_create("EPSG:28992")
    @test source_crs isa Ptr{Nothing}
    @test target_crs isa Ptr{Nothing}
    trans = Proj.Transformation(source_crs, target_crs)

    # Check that altitude and time inputs are correctly forwarded from the transformation
    p = trans(SA_F64[52.16, 5.39, 5, 2020])
    @test is_approx(p, (155191.3538124342, 463537.1362732911, 5.0, 2020.0))

    a = (52.16, 5.39)
    b = Proj.proj_trans(trans.pj, Proj.PJ_FWD, a)
    @test is_approx(b, (155191.3538124342, 463537.1362732911))

    # with always_xy = true, we need to use lon/lat, and still get x/y out
    trans = Proj.Transformation(source_crs, target_crs, always_xy = true)
    a = (5.39, 52.16)
    b = Proj.proj_trans(trans.pj, Proj.PJ_FWD, a)
    @test is_approx(b, (155191.3538124342, 463537.1362732911))
end

@testset "dense 4D coord vector transformation" begin
    source_crs = Proj.proj_create("EPSG:4326")
    target_crs = Proj.proj_create("EPSG:28992")
    trans = Proj.Transformation(source_crs, target_crs, always_xy = true)
    # This array is mutated in place. Note that this array needs to have 4D elements,
    # with 2D elements it will only do every other one
    A = [Proj.Coord(5.39, 52.16) for _ = 1:5]
    err = Proj.proj_trans_array(trans.pj, Proj.PJ_FWD, length(A), A)
    @test err == 0
    B = [Proj.Coord(155191.3538124342, 463537.1362732911) for _ = 1:5]
    @test all(is_approx(a, b) for (a, b) in zip(A, B))

    # since A is not in target_crs, we will get an error on a second call
    err = Proj.PROJError("push: Invalid latitude")
    @test_throws err Proj.proj_trans_array(trans.pj, Proj.PJ_FWD, length(A), A)
    # error number is reset during PROJError construction
    @test Proj.proj_errno(trans.pj) == 0
    # test triggering finalizer manually
    finalize(trans)
    @test trans.pj === C_NULL
end

@testset "generic array transformation" begin
    source_crs = Proj.proj_create("EPSG:4326")
    target_crs = Proj.proj_create("EPSG:28992")
    trans = Proj.Transformation(source_crs, target_crs, always_xy = true)

    # inplace transformation of vector of 2D coordinates
    # using https://proj.org/development/reference/functions.html#c.proj_trans_generic
    A = [SA[5.39, 52.16] for _ = 1:5]
    st = sizeof(first(A))
    ptr = pointer(A)
    n = length(A)
    # 8 is sizeof(Float64), so ptr + 8 points to first latitude element
    n_done = Proj.proj_trans_generic(
        trans.pj,
        Proj.PJ_FWD,
        ptr,
        st,
        n,
        ptr + 8,
        st,
        n,
        C_NULL,
        C_NULL,
        C_NULL,
        C_NULL,
        C_NULL,
        C_NULL,
    )
    @test n_done == n
    B = [SA[155191.3538124342, 463537.1362732911] for _ = 1:5]
    @test A ≈ B
end

@testset "compose" begin
    trans1 = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy = true)
    trans2 = Proj.Transformation("EPSG:32632", "EPSG:2027", always_xy = true)
    trans = trans1 ∘ trans2  # same as compose(trans1, trans2)
    source_crs = Proj.proj_get_source_crs(trans.pj)
    target_crs = Proj.proj_get_target_crs(trans.pj)

    # trans1 source is the new source
    source_crs1 = Proj.proj_get_source_crs(trans1.pj)
    source_wkt = Proj.proj_as_wkt(source_crs, Proj.PJ_WKT2_2019, C_NULL)
    source_wkt1 = Proj.proj_as_wkt(source_crs1, Proj.PJ_WKT2_2019, C_NULL)
    @test source_wkt == source_wkt1

    # trans2 source is the new target
    target_crs1 = Proj.proj_get_target_crs(trans2.pj)
    target_wkt = Proj.proj_as_wkt(target_crs, Proj.PJ_WKT2_2019, C_NULL)
    target_wkt1 = Proj.proj_as_wkt(target_crs1, Proj.PJ_WKT2_2019, C_NULL)
    @test target_wkt == target_wkt1

    # which we can also see from show
    @test repr(trans) == """
    Transformation
        source: WGS 84 (with axis order normalized for visualization)
        target: NAD27(76) / UTM zone 15N
        direction: forward
    """
end

@testset "Coord constructor" begin
    coord = Proj.Coord(1.0, 2.0, 3.0, Inf)
    @test_throws MethodError Proj.Coord()
    @test_throws ErrorException Proj.Coord(0.0)
    @test Proj.Coord(1.0, 2.0) == Proj.Coord(1.0, 2.0, 0.0, Inf)
    @test Proj.Coord(1.0, 2.0, 3.0) == coord
    @test Proj.Coord(1.0, 2.0, 3.0, 4.0) == Proj.Coord(1.0, 2.0, 3.0, 4.0)

    @test Proj.Coord((1.0, 2.0, 3.0)) == coord
    @test Proj.Coord([1.0, 2.0, 3.0]) == coord
    @test Proj.Coord(1, 2, 3) == coord
    @test Proj.Coord([1, 2, 3]) == coord
    @test Proj.Coord(UInt8[1, 2, 3]) == coord
    @test Proj.Coord(1:3) == coord
    @test Proj.Coord(1.0f0:3) == coord
    @test Proj.Coord((i for i = 1:3)) == coord
    @test Proj.Coord(coord) === coord
    @test convert(Proj.Coord, coord) === coord
    @test convert(Proj.Coord, (1.0, 2.0, 3.0)) == coord

    @inferred Proj.Coord(1, 2)
    @inferred Proj.Coord((1, 2))
    @inferred Proj.Coord([1, 2])
    @inferred Proj.Coord(SA[1, 2])
end

@testset "in and output types" begin
    trans = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy = true)
    trans(Proj.proj_coord(5.39, 52.16))
    b = trans(SA[5.39, 52.16, 0.0, 0.0])

    # StaticVector, Vector or Tuple
    x, y, z, t = 5.39, 52.16, 0.0, Inf

    @inferred trans(x, y)
    @inferred trans(x, y, z)
    @inferred trans(x, y, z, t)
    @inferred trans(Proj.Coord(x, y))
    @inferred trans((x, y))
    # for Vector it cannot be precisely inferred due to unknown length
    @inferred Proj.NTuple234 trans([x, y])
    @inferred trans(SA[x, y])

    # SVector{2/3/4, Float64}
    p = trans(SA[5, 52])
    @test is_approx(p, (128410.08537081012, 445806.50883314764))
    p = trans(SA[5.39, 52.16, 2.0])
    @test is_approx(p, (155191.35381147722, 463537.13624246384, 2.0))
    p = trans(SA[5.39, 52.16, 0.0, Inf])
    @test is_approx(p, (155191.3538124342, 463537.1362732911, 0.0, Inf))
end

@testset "network" begin
    as_before = Proj.network_enabled()

    # turn off network, no z transformation
    @test !Proj.enable_network!(false)
    @test !Proj.network_enabled()
    trans_z = Proj.Transformation("EPSG:4326+5773", "EPSG:7856+5711", always_xy = true)
    @test trans_z((151, -33, 5))[3] == 5
    # turn on network, z transformation
    @test Proj.enable_network!(true)
    @test Proj.network_enabled()
    trans_z = Proj.Transformation("EPSG:4326+5773", "EPSG:7856+5711", always_xy = true)
    z = trans_z((151, -33, 5))[3]
    @test z ≈ 5.280603319143665

    # 0 args turns it on as well
    Proj.enable_network!(false)
    Proj.enable_network!()
    @test Proj.network_enabled()

    # restore setting as outside the testset
    Proj.enable_network!(as_before)
end

ext: jl


text:
function build_attitude_provider(obj::SatelliteType, arg0::AbsoluteDate, arg1::AbsoluteDate, arg2::ExtendedPVCoordinatesProvider, arg3::Frame, arg4::jint)
    return jcall(obj, "buildAttitudeProvider", GNSSAttitudeProvider, (AbsoluteDate, AbsoluteDate, ExtendedPVCoordinatesProvider, Frame, jint), arg0, arg1, arg2, arg3, arg4)
end

function compare_to(obj::Enum, arg0::Enum)
    return jcall(obj, "compareTo", jint, (Enum,), arg0)
end

function describe_constable(obj::Enum)
    return jcall(obj, "describeConstable", Optional, ())
end

function equals(obj::Enum, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_declaring_class(obj::Enum)
    return jcall(obj, "getDeclaringClass", Class, ())
end

function get_name(obj::SatelliteType)
    return jcall(obj, "getName", JString, ())
end

function hash_code(obj::Enum)
    return jcall(obj, "hashCode", jint, ())
end

function name(obj::Enum)
    return jcall(obj, "name", JString, ())
end

function ordinal(obj::Enum)
    return jcall(obj, "ordinal", jint, ())
end

function parse_satellite_type(::Type{SatelliteType}, arg0::JString)
    return jcall(SatelliteType, "parseSatelliteType", SatelliteType, (JString,), arg0)
end

function to_string(obj::Enum)
    return jcall(obj, "toString", JString, ())
end

function value_of(::Type{Enum}, arg0::Class, arg1::JString)
    return jcall(Enum, "valueOf", Enum, (Class, JString), arg0, arg1)
end

function value_of(::Type{SatelliteType}, arg0::JString)
    return jcall(SatelliteType, "valueOf", SatelliteType, (JString,), arg0)
end

function values(::Type{SatelliteType})
    return jcall(SatelliteType, "values", Vector{SatelliteType}, ())
end


text:
module REPLMode

using Markdown
using UUIDs

import REPL
import REPL: LineEdit, REPLCompletions

import ..devdir, ..print_first_command_header, ..API
using ..Types, ..Display, ..Operations

############
# Commands #
############
@enum(CommandKind, CMD_HELP, CMD_STATUS, CMD_SEARCH, CMD_ADD, CMD_RM, CMD_UP,
                   CMD_TEST, CMD_GC, CMD_PREVIEW, CMD_INIT, CMD_BUILD, CMD_FREE,
                   CMD_PIN, CMD_CHECKOUT, CMD_DEVELOP, CMD_GENERATE, CMD_PRECOMPILE)

struct Command
    kind::CommandKind
    val::String
end
Base.show(io::IO, cmd::Command) = print(io, cmd.val)

const cmds = Dict(
    "help"      => CMD_HELP,
    "?"         => CMD_HELP,
    "status"    => CMD_STATUS,
    "st"        => CMD_STATUS,
    "."         => CMD_STATUS,
    "search"    => CMD_SEARCH,
    "find"      => CMD_SEARCH,
    "/"         => CMD_SEARCH,
    "add"       => CMD_ADD,
    "install"   => CMD_ADD,
    "+"         => CMD_ADD,
    "rm"        => CMD_RM,
    "remove"    => CMD_RM,
    "uninstall" => CMD_RM,
    "-"         => CMD_RM,
    "up"        => CMD_UP,
    "update"    => CMD_UP,
    "upgrade"   => CMD_UP,
    "test"      => CMD_TEST,
    "gc"        => CMD_GC,
    "preview"   => CMD_PREVIEW,
    "init"      => CMD_INIT,
    "build"     => CMD_BUILD,
    "pin"       => CMD_PIN,
    "free"      => CMD_FREE,
    "checkout"  => CMD_CHECKOUT, # deprecated
    "develop"   => CMD_DEVELOP,
    "dev"       => CMD_DEVELOP,
    "generate"  => CMD_GENERATE,
    "precompile" => CMD_PRECOMPILE,
)

#################
# Git revisions #
#################
struct Rev
    rev::String
end

###########
# Options #
###########
@enum(OptionKind, OPT_ENV, OPT_PROJECT, OPT_MANIFEST, OPT_MAJOR, OPT_MINOR,
                  OPT_PATCH, OPT_FIXED, OPT_COVERAGE, OPT_NAME, OPT_PATH)

function Types.PackageMode(opt::OptionKind)
    opt == OPT_MANIFEST && return PKGMODE_MANIFEST
    opt == OPT_PROJECT  && return PKGMODE_PROJECT
    throw(ArgumentError("invalid option $opt"))
end

function Types.UpgradeLevel(opt::OptionKind)
    opt == OPT_MAJOR && return UPLEVEL_MAJOR
    opt == OPT_MINOR && return UPLEVEL_MINOR
    opt == OPT_PATCH && return UPLEVEL_PATCH
    opt == OPT_FIXED && return UPLEVEL_FIXED
    throw(ArgumentError("invalid option $opt"))
end

struct Option
    kind::OptionKind
    val::String
    argument::Union{String, Nothing}
    Option(kind::OptionKind, val::String) = new(kind, val, nothing)
    function Option(kind::OptionKind, val::String, argument::Union{String, Nothing})
        if kind in (OPT_PROJECT, OPT_MANIFEST, OPT_MAJOR,
                    OPT_MINOR, OPT_PATCH, OPT_FIXED) &&
                argument !== nothing
            cmderror("the `$val` option does not take an argument")
        elseif kind in (OPT_ENV, OPT_PATH) && argument == nothing
            cmderror("the `$val` option requires an argument")
        end
        if kind == OPT_PATH
            argument =  replace(argument, "~" => homedir())
        end
        new(kind, val, argument)
    end
end
Base.show(io::IO, opt::Option) = print(io, "--$(opt.val)", opt.argument == nothing ? "" : "=$(opt.argument)")

const opts = Dict(
    "env"      => OPT_ENV,
    "project"  => OPT_PROJECT,
    "p"        => OPT_PROJECT,
    "manifest" => OPT_MANIFEST,
    "m"        => OPT_MANIFEST,
    "major"    => OPT_MAJOR,
    "minor"    => OPT_MINOR,
    "patch"    => OPT_PATCH,
    "fixed"    => OPT_FIXED,
    "coverage" => OPT_COVERAGE,
    "name"     => OPT_NAME,
    "path"     => OPT_PATH,
)

function parse_option(word::AbstractString)::Option
    m = match(r"^(?: -([a-z]) | --([a-z]{2,})(?:\s*=\s*(\S*))? )$"ix, word)
    m == nothing && cmderror("invalid option: ", repr(word))
    k = m.captures[1] != nothing ? m.captures[1] : m.captures[2]
    haskey(opts, k) || cmderror("invalid option: ", repr(word))
    return Option(opts[k], String(k), m.captures[3] == nothing ? nothing : String(m.captures[3]))
end

###################
# Package parsing #
###################
let uuid = raw"(?i)[0-9a-z]{8}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{4}-[0-9a-z]{12}(?-i)",
    name = raw"(\w+)(?:\.jl)?"
    global const name_re = Regex("^$name\$")
    global const uuid_re = Regex("^$uuid\$")
    global const name_uuid_re = Regex("^$name\\s*=\\s*($uuid)\$")
end

function parse_package(word::AbstractString; context=nothing)::PackageSpec
    word = replace(word, "~" => homedir())
    if context in (CMD_ADD, CMD_DEVELOP) && isdir(word)
        pkg = PackageSpec()
        pkg.repo = Types.GitRepo(abspath(word))
        return pkg
    elseif occursin(uuid_re, word)
        return PackageSpec(UUID(word))
    elseif occursin(name_re, word)
        return PackageSpec(String(match(name_re, word).captures[1]))
    elseif occursin(name_uuid_re, word)
        m = match(name_uuid_re, word)
        return PackageSpec(String(m.captures[1]), UUID(m.captures[2]))
    else
        if context in (CMD_ADD, CMD_DEVELOP)
            # Guess it is a url then
            pkg = PackageSpec()
            pkg.repo = Types.GitRepo(word)
            return pkg
        else
            cmderror("`$word` cannot be parsed as a package")
        end
    end
end

################
# REPL parsing #
################
const lex_re = r"^[\?\./\+\-](?!\-) | [^@\#\s;]+\s*=\s*[^@\#\s;]+ | \#\s*[^@\#\s;]* | @\s*[^@\#\s;]* | [^@\#\s;]+|;"x

const Token = Union{Command, Option, VersionRange, String, Rev}

function tokenize(cmd::String)::Vector{Vector{Token}}
    words = map(m->m.match, eachmatch(lex_re, cmd))
    commands = Vector{Token}[]
    while !isempty(words)
        push!(commands, tokenize!(words))
    end
    return commands
end

function tokenize!(words::Vector{<:AbstractString})::Vector{Token}
    print_first_command_header()
    tokens = Token[]
    help_mode = false
    preview_mode = false
    # First parse a Command or a modifier (help / preview) + Command
    while !isempty(words)
        word = popfirst!(words)
        if word[1] == '-' && length(word) > 1
            push!(tokens, parse_option(word))
        else
            haskey(cmds, word) || cmderror("invalid command: ", repr(word))
            cmdkind = cmds[word]
            push!(tokens, Command(cmdkind, word))
            # If help / preview and not in help mode we want to eat another cmd
            if !help_mode
                cmdkind == CMD_HELP    && (help_mode    = true; continue)
                cmdkind == CMD_PREVIEW && (preview_mode = true; continue)
            end
            break
        end
    end
    if isempty(tokens) || !(tokens[end] isa Command)
        cmderror("no package command given")
    end
    # Now parse the arguments / options to the command
    while !isempty(words)
        word = popfirst!(words)
        if word == ";"
            return tokens
        elseif first(word) == '-'
            push!(tokens, parse_option(word))
        elseif first(word) == '@'
            push!(tokens, VersionRange(strip(word[2:end])))
        elseif first(word) == '#'
            push!(tokens, Rev(word[2:end]))
        else
            push!(tokens, String(word))
        end
    end
    return tokens
end

#############
# Execution #
#############

function do_cmd(repl::REPL.AbstractREPL, input::String; do_rethrow=false)
    try
        commands = tokenize(input)
        for command in commands
            do_cmd!(command, repl)
        end
    catch err
        if do_rethrow
            rethrow(err)
        end
        if err isa CommandError
            Base.display_error(repl.t.err_stream, ErrorException(err.msg), Ptr{Nothing}[])
        else
            Base.display_error(repl.t.err_stream, err, Base.catch_backtrace())
        end
    end
end

function do_cmd!(tokens::Vector{Token}, repl)
    cmd = env_opt = nothing
    while !isempty(tokens)
        token = popfirst!(tokens)
        if token isa Command
            cmd = token
            break
        elseif token isa Option
            # Only OPT_ENV is allowed before a command
            if token.kind == OPT_ENV
                env_opt = Base.parse_env(token.argument)
            else
                cmderror("unrecognized command option: `$token`")
            end
        else
            cmderror("misplaced token: ", token)
        end
    end
    ctx = Context(env = EnvCache(env_opt))
    if cmd.kind == CMD_PREVIEW
        ctx.preview = true
        isempty(tokens) && cmderror("expected a command to preview")
        cmd = popfirst!(tokens)
    end
    # Using invokelatest to hide the functions from inference.
    # Otherwise it would try to infer everything here.
    cmd.kind == CMD_INIT     ? Base.invokelatest(          do_init!, ctx, tokens) :
    cmd.kind == CMD_HELP     ? Base.invokelatest(          do_help!, ctx, tokens, repl) :
    cmd.kind == CMD_RM       ? Base.invokelatest(            do_rm!, ctx, tokens) :
    cmd.kind == CMD_ADD      ? Base.invokelatest(do_add_or_develop!, ctx, tokens, CMD_ADD) :
    cmd.kind == CMD_CHECKOUT ? Base.invokelatest(do_add_or_develop!, ctx, tokens, CMD_DEVELOP) :
    cmd.kind == CMD_DEVELOP  ? Base.invokelatest(do_add_or_develop!, ctx, tokens, CMD_DEVELOP) :
    cmd.kind == CMD_UP       ? Base.invokelatest(            do_up!, ctx, tokens) :
    cmd.kind == CMD_STATUS   ? Base.invokelatest(        do_status!, ctx, tokens) :
    cmd.kind == CMD_TEST     ? Base.invokelatest(          do_test!, ctx, tokens) :
    cmd.kind == CMD_GC       ? Base.invokelatest(            do_gc!, ctx, tokens) :
    cmd.kind == CMD_BUILD    ? Base.invokelatest(         do_build!, ctx, tokens) :
    cmd.kind == CMD_PIN      ? Base.invokelatest(           do_pin!, ctx, tokens) :
    cmd.kind == CMD_FREE     ? Base.invokelatest(          do_free!, ctx, tokens) :
    cmd.kind == CMD_GENERATE ? Base.invokelatest(      do_generate!, ctx, tokens) :
    cmd.kind == CMD_PRECOMPILE ? Base.invokelatest(  do_precompile!, ctx, tokens) :
        cmderror("`$cmd` command not yet implemented")
    return
end

const help = md"""

**Welcome to the Pkg3 REPL-mode**. To return to the `julia>` prompt, either press
backspace when the input line is empty or press Ctrl+C.

**Synopsis**

    pkg> [--env=...] cmd [opts] [args]

Multiple commands can be given on the same line by interleaving a `;` between the commands.

**Environment**

The `--env` meta option determines which project environment to manipulate. By
default, this looks for a git repo in the parents directories of the current
working directory, and if it finds one, it uses that as an environment. Otherwise,
it uses a named environment (typically found in `~/.julia/environments`) looking
for environments named `v$(VERSION.major).$(VERSION.minor).$(VERSION.patch)`,
`v$(VERSION.major).$(VERSION.minor)`,  `v$(VERSION.major)` or `default` in order.

**Commands**

What action you want the package manager to take:

`help`: show this message

`status`: summarize contents of and changes to environment

`generate`: generate files for a new project

`add`: add packages to project

`rm`: remove packages from project or manifest

`up`: update packages in manifest

`preview`: previews a subsequent command without affecting the current state

`test`: run tests for packages

`gc`: garbage collect packages not used for a significant time

`init`: initializes an environment in the current, or git base, directory

`build`: run the build script for packages

`pin`: pins the version of packages

`develop`: clone the full package repo locally for development

`free`: undos a `pin` or `develop`

`precompile`: precompile all the project dependencies
"""

const helps = Dict(
    CMD_HELP => md"""

        help

    Display this message.

        help cmd ...

    Display usage information for commands listed.

    Available commands: `help`, `status`, `add`, `rm`, `up`, `preview`, `gc`, `test`, `init`, `build`, `free`, `pin`, `develop`.
    """, CMD_STATUS => md"""

        status
        status [-p|--project]
        status [-m|--manifest]

    Show the status of the current environment. By default, the full contents of
    the project file is summarized, showing what version each package is on and
    how it has changed since the last git commit (if in a git repo), as well as
    any changes to manifest packages not already listed. In `--project` mode, the
    status of the project file is summarized. In `--project` mode, the status of
    the project file is summarized.
    """, CMD_GENERATE => md"""

        create name

    Create a project called `name` in the current folder.
    """,
    CMD_ADD => md"""

        add pkg[=uuid] [@version] [#rev] ...

    Add package `pkg` to the current project file. If `pkg` could refer to
    multiple different packages, specifying `uuid` allows you to disambiguate.
    `@version` optionally allows specifying which versions of packages. Versions
    may be specified by `@1`, `@1.2`, `@1.2.3`, allowing any version with a prefix
    that matches, or ranges thereof, such as `@1.2-3.4.5`. A git-revision can be
    specified by `#branch` or `#commit`.

    **Examples**
    ```
    pkg> add Example
    pkg> add Example@0.5
    pkg> add Example#master
    pkg> add Example#c37b675
    pkg> add https://github.com/JuliaLang/Example.jl#master
    pkg> add Example=7876af07-990d-54b4-ab0e-23690620f79a
    ```
    """, CMD_RM => md"""

        rm [-p|--project] pkg[=uuid] ...

    Remove package `pkg` from the project file. Since the name `pkg` can only
    refer to one package in a project this is unambiguous, but you can specify
    a `uuid` anyway, and the command is ignored, with a warning if package name
    and UUID do not mactch. When a package is removed from the project file, it
    may still remain in the manifest if it is required by some other package in
    the project. Project mode operation is the default, so passing `-p` or
    `--project` is optional unless it is preceded by the `-m` or `--manifest`
    options at some earlier point.

        rm [-m|--manifest] pkg[=uuid] ...

    Remove package `pkg` from the manifest file. If the name `pkg` refers to
    multiple packages in the manifest, `uuid` disambiguates it. Removing a package
    from the manifest forces the removal of all packages that depend on it, as well
    as any no-longer-necessary manifest packages due to project package removals.
    """, CMD_UP => md"""

        up [-p|project]  [opts] pkg[=uuid] [@version] ...
        up [-m|manifest] [opts] pkg[=uuid] [@version] ...

        opts: --major | --minor | --patch | --fixed

    Update the indicated package within the constraints of the indicated version
    specifications. Versions may be specified by `@1`, `@1.2`, `@1.2.3`, allowing
    any version with a prefix that matches, or ranges thereof, such as `@1.2-3.4.5`.
    In `--project` mode, package specifications only match project packages, while
    in `manifest` mode they match any manifest package. Bound level options force
    the following packages to be upgraded only within the current major, minor,
    patch version; if the `--fixed` upgrade level is given, then the following
    packages will not be upgraded at all.
    """, CMD_PREVIEW => md"""

        preview cmd

    Runs the command `cmd` in preview mode. This is defined such that no side effects
    will take place i.e. no packages are downloaded and neither the project nor manifest
    is modified.
    """, CMD_TEST => md"""

        test [opts] pkg[=uuid] ...

        opts: --coverage

    Run the tests for package `pkg`. This is done by running the file `test/runtests.jl`
    in the package directory. The option `--coverage` can be used to run the tests with
    coverage enabled.
    """, CMD_GC => md"""

    Deletes packages that are not reached from any environment used within the last 6 weeks.
    """, CMD_INIT => md"""

        init

    Creates an environment in the current directory, or the git base directory if the current directory
    is in a git repository.
    """, CMD_BUILD =>md"""

        build pkg[=uuid] ...

    Run the build script in deps/build.jl for each package in `pkg`` and all of their dependencies in depth-first recursive order.
    If no packages are given, runs the build scripts for all packages in the manifest.
    """, CMD_PIN => md"""

        pin pkg[=uuid] ...

    Pin packages to given versions, or the current version if no version is specified. A pinned package has its version fixed and will not be upgraded or downgraded.
    A pinned package has the symbol `⚲` next to its version in the status list.
    """, CMD_FREE => md"""
        free pkg[=uuid] ...

    Free a pinned package `pkg`, which allows it to be upgraded or downgraded again. If the package is checked out (see `help develop`) then this command
    makes the package no longer being checked out.
    """, CMD_DEVELOP => md"""
        develop pkg[=uuid] [#rev] ...

    Make a package available for development. If `pkg` is an existing local path that path will be recorded in
    the manifest and used. Otherwise, a full git clone of `pkg` at rev `rev` is made. The clone is stored in `devdir`,
    which defaults to `~/.julia/dev` and is set by the environment variable `JULIA_PKG_DEVDIR`.
    This operation is undone by `free`.

    *Example*
    ```jl
    pkg> develop Example
    pkg> develop Example#master
    pkg> develop Example#c37b675
    pkg> develop https://github.com/JuliaLang/Example.jl#master
    ```
    """, CMD_PRECOMPILE => md"""
        precompile

    Precompile all the dependencies of the project by running `import` on all of them in a new process.
    """
)

function do_help!(
    ctk::Context,
    tokens::Vector{Token},
    repl::REPL.AbstractREPL,
)
    disp = REPL.REPLDisplay(repl)
    if isempty(tokens)
        Base.display(disp, help)
        return
    end
    help_md = md""
    for token in tokens
        if token isa Command
            if haskey(helps, token.kind)
                isempty(help_md.content) ||
                push!(help_md.content, md"---")
                push!(help_md.content, helps[token.kind].content)
            else
                cmderror("Sorry, I don't have any help for the `$(token.val)` command.")
            end
        else
            error("invalid usage of help command")
        end
    end
    Base.display(disp, help_md)
end

function do_rm!(ctx::Context, tokens::Vector{Token})
    # tokens: package names and/or uuids
    mode = PKGMODE_PROJECT
    pkgs = PackageSpec[]
    while !isempty(tokens)
        token = popfirst!(tokens)
        if token isa String
            push!(pkgs, parse_package(token))
            pkgs[end].mode = mode
        elseif token isa VersionRange
            cmderror("`rm` does not take version specs")
        elseif token isa Option
            if token.kind in (OPT_PROJECT, OPT_MANIFEST)
                mode = PackageMode(token.kind)
            else
                cmderror("invalid option for `rm`: $token")
            end
        end
    end
    isempty(pkgs) &&
        cmderror("`rm` – list packages to remove")
    API.rm(ctx, pkgs)
end

function do_add_or_develop!(ctx::Context, tokens::Vector{Token}, cmd::CommandKind)
    @assert cmd in (CMD_ADD, CMD_DEVELOP)
    # tokens: package names and/or uuids, optionally followed by version specs
    isempty(tokens) &&
        cmderror("`add` – list packages to add")
    pkgs = PackageSpec[]
    prev_token_was_package = false
    while !isempty(tokens)
        parsed_package = false
        token = popfirst!(tokens)
        if token isa String
            push!(pkgs, parse_package(token; context=CMD_ADD))
            cmd == CMD_DEVELOP && pkgs[end].repo == nothing && (pkgs[end].repo = Types.GitRepo("", ""))
            parsed_package = true
        elseif token isa VersionRange
            prev_token_was_package ||
                cmderror("package name/uuid must precede version spec `@$token`")
            pkgs[end].version = VersionSpec(token)
        elseif token isa Rev
            prev_token_was_package ||
                cmderror("package name/uuid must precede rev spec `#$(token.rev)`")
            # WE did not get the repo from the
            pkg = pkgs[end]
            if pkg.repo == nothing
                pkg.repo = Types.GitRepo("", token.rev)
            else
                pkgs[end].repo.rev = token.rev
            end
        elseif token isa Option
            cmderror("`add` doesn't take options: $token")
        end
        prev_token_was_package = parsed_package
    end
    return API.add_or_develop(ctx, pkgs, mode=(cmd == CMD_ADD ? :add : :develop))
end

function do_up!(ctx::Context, tokens::Vector{Token})
    # tokens:
    #  - upgrade levels as options: --[fixed|patch|minor|major]
    #  - package names and/or uuids, optionally followed by version specs
    pkgs = PackageSpec[]
    mode = PKGMODE_PROJECT
    level = UPLEVEL_MAJOR
    prev_token_was_package = false
    while !isempty(tokens)
        parsed_package = false
        token = popfirst!(tokens)
        if token isa String
            push!(pkgs, parse_package(token))
            pkgs[end].version = level
            pkgs[end].mode = mode
            parsed_package = true
        elseif token isa VersionRange
            prev_token_was_package ||
                cmderror("package name/uuid must precede version spec `@$token`")
            pkgs[end].version = VersionSpec(token)
        elseif token isa Option
            if token.kind in (OPT_PROJECT, OPT_MANIFEST)
                mode = PackageMode(token.kind)
            elseif token.kind in (OPT_MAJOR, OPT_MINOR, OPT_PATCH, OPT_FIXED)
                level = UpgradeLevel(token.kind)
            else
                cmderror("invalid option for `up`: $(token)")
            end
        end
        prev_token_was_package = parsed_package
    end
    API.up(ctx, pkgs; level=level, mode=mode)
end

function do_pin!(ctx::Context, tokens::Vector{Token})
    pkgs = PackageSpec[]
    prev_token_was_package = false
    while !isempty(tokens)
        token = popfirst!(tokens)
        parsed_package = false
        if token isa String
            push!(pkgs, parse_package(token))
            parsed_package = true
        elseif token isa VersionRange
            prev_token_was_package ||
                cmderror("package name/uuid must precede version spec `@$token`")
            if token.lower != token.upper
                cmderror("pinning a package requires a single version, not a versionrange")
            end
            pkgs[end].version = VersionSpec(token)
        else
            cmderror("free only takes a list of packages ")
        end
        prev_token_was_package = parsed_package
    end
    API.pin(ctx, pkgs)
end

function do_free!(ctx::Context, tokens::Vector{Token})
    pkgs = PackageSpec[]
    while !isempty(tokens)
        token = popfirst!(tokens)
        if token isa String
            push!(pkgs, parse_package(token))
        else
            cmderror("free only takes a list of packages")
        end
    end
    API.free(ctx, pkgs)
end

function do_checkout!(ctx::Context, tokens::Vector{Token})
    Base.depwarn("`checkout`` is deprecated, use `develop`", :checkout)
    do_develop!(ctx, tokens)
end

function do_status!(ctx::Context, tokens::Vector{Token})
    mode = PKGMODE_COMBINED
    while !isempty(tokens)
        token = popfirst!(tokens)
        if token isa Option
            if token.kind in (OPT_PROJECT, OPT_MANIFEST)
                mode = PackageMode(token.kind)
            else
                cmderror("invalid option for `status`: $(token)")
            end
        else
            cmderror("`status` does not take arguments")
        end
    end
    Display.status(ctx, mode)
end

# TODO , test recursive dependencies as on option.
function do_test!(ctx::Context, tokens::Vector{Token})
    pkgs = PackageSpec[]
    coverage = false
    while !isempty(tokens)
        token = popfirst!(tokens)
        if token isa String
            pkg = parse_package(token)
            pkg.mode = PKGMODE_MANIFEST
            push!(pkgs, pkg)
        elseif token isa Option
            if token.kind == OPT_COVERAGE
                coverage = true
            else
                cmderror("invalid option for `test`: $token")
            end
        else
            # TODO: Better error message
            cmderror("invalid usage for `test`")
        end
    end
    API.test(ctx, pkgs; coverage = coverage)
end

function do_gc!(ctx::Context, tokens::Vector{Token})
    !isempty(tokens) && cmderror("`gc` does not take any arguments")
    API.gc(ctx)
end

function do_build!(ctx::Context, tokens::Vector{Token})
    pkgs = PackageSpec[]
    while !isempty(tokens)
        token = popfirst!(tokens)
        if token isa String
            push!(pkgs, parse_package(token))
        else
            cmderror("`build` only takes a list of packages")
        end
    end
    API.build(ctx, pkgs)
end

function do_init!(ctx::Context, tokens::Vector{Token})
    if !isempty(tokens)
        cmderror("`init` does currently not take any arguments")
    end
    API.init(ctx)
end

function do_generate!(ctx::Context, tokens::Vector{Token})
    local pkg
    while !isempty(tokens)
        token = popfirst!(tokens)
        if token isa String
            pkg = token
            break # TODO: error message?
        else
            cmderror("`generate` takes a name of the project to create")
        end
    end
    API.generate(pkg)
end

function do_precompile!(ctx::Context, tokens::Vector{Token})
    if !isempty(tokens)
        cmderror("`precompile` does not take any arguments")
    end
    API.precompile(ctx)
end

######################
# REPL mode creation #
######################

# Provide a string macro pkg"cmd" that can be used in the same way
# as the REPLMode `pkg> cmd`. Useful for testing and in environments
# where we do not have a REPL, e.g. IJulia.
struct MiniREPL <: REPL.AbstractREPL
    display::TextDisplay
    t::REPL.Terminals.TTYTerminal
end
function MiniREPL()
    MiniREPL(TextDisplay(stdout), REPL.Terminals.TTYTerminal(get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"), stdin, stdout, stderr))
end
REPL.REPLDisplay(repl::MiniREPL) = repl.display

__init__() = minirepl[] = MiniREPL()

const minirepl = Ref{MiniREPL}()

macro pkg_str(str::String)
    :($(do_cmd)(minirepl[], $str; do_rethrow=true))
end

pkgstr(str::String) = do_cmd(minirepl[], str; do_rethrow=true)

# handle completions
all_commands_sorted = sort!(collect(keys(cmds)))
long_commands = filter(c -> length(c) > 2, all_commands_sorted)
all_options_sorted = [length(opt) > 1 ? "--$opt" : "-$opt" for opt in sort!(collect(keys(opts)))]
long_options = filter(c -> length(c) > 2, all_options_sorted)

struct PkgCompletionProvider <: LineEdit.CompletionProvider end

function LineEdit.complete_line(c::PkgCompletionProvider, s)
    partial = REPL.beforecursor(s.input_buffer)
    full = LineEdit.input_string(s)
    ret, range, should_complete = completions(full, lastindex(partial))
    return ret, partial[range], should_complete
end

function complete_command(s, i1, i2)
    # only show short form commands when no input is given at all
    cmp = filter(cmd -> startswith(cmd, s), isempty(s) ? all_commands_sorted : long_commands)
    return cmp, i1:i2, !isempty(cmp)
end

function complete_option(s, i1, i2)
    # only show short form options if only a dash is given
    cmp = filter(cmd -> startswith(cmd, s), length(s) == 1 && first(s) == '-' ?
                                                all_options_sorted :
                                                long_options)
    return cmp, i1:i2, !isempty(cmp)
end

function complete_package(s, i1, i2, lastcommand, project_opt)
    if lastcommand in [CMD_STATUS, CMD_RM, CMD_UP, CMD_TEST, CMD_BUILD, CMD_FREE, CMD_PIN, CMD_CHECKOUT]
        return complete_installed_package(s, i1, i2, project_opt)
    elseif lastcommand in [CMD_ADD, CMD_DEVELOP]
        return complete_remote_package(s, i1, i2)
    end
    return String[], 0:-1, false
end

function complete_installed_package(s, i1, i2, project_opt)
    pkgs = project_opt ? API.installed(PKGMODE_PROJECT) : API.installed()
    pkgs = sort!(collect(keys(filter((p) -> p[2] != nothing, pkgs))))
    cmp = filter(cmd -> startswith(cmd, s), pkgs)
    return cmp, i1:i2, !isempty(cmp)
end

function complete_remote_package(s, i1, i2)
    cmp = filter(cmd -> startswith(cmd, s), collect_package_names())
    return cmp, i1:i2, !isempty(cmp)
end

function collect_package_names()
    r = r"name = \"(.*?)\""
    names = String[]
    for reg in Types.registries(;clone_default=false)
        regcontent = read(joinpath(reg, "Registry.toml"), String)
        append!(names, collect(match.captures[1] for match in eachmatch(r, regcontent)))
    end
    return sort!(names)
end

function completions(full, index)
    pre = full[1:index]

    pre_words = split(pre, ' ', keepempty=true)

    # first word should always be a command
    if isempty(pre_words)
        return complete_command("", 1:1)
    else
        to_complete = pre_words[end]
        offset = isempty(to_complete) ? index+1 : to_complete.offset+1

        if length(pre_words) == 1
            return complete_command(to_complete, offset, index)
        end

        # tokenize input, don't offer any completions for invalid commands
        tokens = try
            tokenize(join(pre_words[1:end-1], ' '))[end]
        catch
            return String[], 0:-1, false
        end

        tokens = reverse!(tokens)

        lastcommand = nothing
        project_opt = true
        for t in tokens
            if t isa Command
                lastcommand = t.kind
                break
            end
        end
        for t in tokens
            if t isa Option && t.kind in [OPT_PROJECT, OPT_MANIFEST]
                project_opt = t.kind == OPT_PROJECT
                break
            end
        end

        if lastcommand in [CMD_HELP, CMD_PREVIEW]
            return complete_command(to_complete, offset, index)
        elseif !isempty(to_complete) && first(to_complete) == '-'
            return complete_option(to_complete, offset, index)
        else
            return complete_package(to_complete, offset, index, lastcommand, project_opt)
        end
    end
end

function promptf()
    env = EnvCache()
    proj_dir = dirname(env.project_file)
    if startswith(pwd(), proj_dir) && env.pkg != nothing && !isempty(env.pkg.name)
        name = env.pkg.name
    else
        name = basename(proj_dir)
    end
    prefix = string("(", name, ") ")
    return prefix * "pkg> "
end

# Set up the repl Pkg REPLMode
function create_mode(repl, main)

    pkg_mode = LineEdit.Prompt(promptf;
        prompt_prefix = Base.text_colors[:blue],
        prompt_suffix = "",
        complete = PkgCompletionProvider(),
        sticky = true)

    pkg_mode.repl = repl
    hp = main.hist
    hp.mode_mapping[:pkg] = pkg_mode
    pkg_mode.hist = hp

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    prefix_prompt, prefix_keymap = LineEdit.setup_prefix_keymap(hp, pkg_mode)

    pkg_mode.on_done = (s, buf, ok) -> begin
        ok || return REPL.transition(s, :abort)
        input = String(take!(buf))
        REPL.reset(repl)
        do_cmd(repl, input)
        REPL.prepare_next(repl)
        REPL.reset_state(s)
        s.current_mode.sticky || REPL.transition(s, main)
    end

    mk = REPL.mode_keymap(main)

    shell_mode = nothing
    for mode in Base.active_repl.interface.modes
        if mode isa LineEdit.Prompt
            mode.prompt == "shell> " && (shell_mode = mode)
        end
    end

    repl_keymap = Dict()
    if shell_mode != nothing
        repl_keymap[';'] = function (s,o...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, shell_mode) do
                    LineEdit.state(s, shell_mode).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, ';')
            end
        end
    end

    b = Dict{Any,Any}[
        skeymap, repl_keymap, mk, prefix_keymap, LineEdit.history_keymap,
        LineEdit.default_keymap, LineEdit.escape_defaults
    ]
    pkg_mode.keymap_dict = LineEdit.keymap(b)
    return pkg_mode
end

function repl_init(repl)
    main_mode = repl.interface.modes[1]
    pkg_mode = create_mode(repl, main_mode)
    push!(repl.interface.modes, pkg_mode)
    keymap = Dict{Any,Any}(
        ']' => function (s,args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, pkg_mode) do
                    LineEdit.state(s, pkg_mode).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, ']')
            end
        end
    )
    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, keymap)
    return
end

end

max_stars_repo_name: BeastyBlacksmith/julia

max_issues_repo_name: BeastyBlacksmith/julia

max_forks_repo_name: BeastyBlacksmith/julia


text:
####
#### Defines list of discussions
####

discussions = Any[
    "Home" => "Discussions/index.md",
    "Atmos" => Any["Discussions/Atmos/EDMFEquations.md",],
    # "Ocean" => Any[],
    # "Land" => Any[],
    # "Common" => Any[],
    # "Numerics" => Any[],
]


text:
# Computes derivatives over the timestep.
using PyPlot

include("../../../src/integrate_lightcurve.jl")

# Test it out:

include("../../../test/loglinspace.jl")

t1 = -1.5; t2 = 1.5; nt = 10000; dt = 0.3
t = zeros(nt)
t .= linearspace(t1,t2,nt)
# The following compares different tolerances and maxdepths:
r = 0.1; b0 = 0.5; u_n = [0.3,0.3]; nu = length(u_n)
favg0 = zeros(5+nu,nt)
favg1 = zeros(5+nu,nt)
neval1= zeros(Int64,nt)
depthmax1= zeros(Int64,nt)
trans = transit_init(r,b0,u_n,true)
param = [0.0,1.0,b0]   # [t_0,v,b_0]

# Now, carry out speed test:

function compute_lightcurve!(trans::Transit_Struct{T},param::Array{T,1},t::Array{T,1},favg0::Array{T,2},nt::Int64) where {T <: Real}
@inbounds for i=1:nt
  trans.b = sqrt(param[3]^2+(param[2]*(t[i]-param[1]))^2)
  favg0[1,i] =  transit_poly!(trans)-1
  favg0[2,i] =  trans.dfdrb[1]
  favg0[3,i] =  trans.dfdrb[2]/trans.b*param[2]^2*(param[1]-t[i])
  favg0[4,i] =  trans.dfdrb[2]*param[2]/trans.b*(t[i]-param[1])^2
  favg0[5,i] =  trans.dfdrb[2]*param[3]/trans.b
  favg0[6:5+trans.n,i] = trans.dfdu
end
return
end

compute_lightcurve!(trans,param,t,favg0,nt)
@time compute_lightcurve!(trans,param,t,favg0,nt)

integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-10,30,neval1,depthmax1)

# Now plot the results:
fig,axes = subplots(2,3,figsize=(12,8))
if VERSION >= v"0.7"
  ax = axes[1]
  ax.plot(t,favg0[1,:].+1,linestyle="--",linewidth=2,label=L"${F}$")
  ax.plot(t,favg1[1,:].+1,linewidth=2,label=L"$\overline{{F}}$")
  ax.plot([-dt/2,-dt/2,-dt/2,dt/2,dt/2,dt/2],[0.9925,0.9935,0.993,0.993,0.9925,0.9935],color="b")
  ax.annotate(L"$\Delta t$",xy=[-0.08,0.991])
  ax.legend(loc="upper center")
  ax.axis([-1.5,1.5,0.988,1.002])
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux")
  ax = axes[2]
  ax.plot(t,favg0[2,:]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial r}$")
  ax.plot(t,favg1[2,:]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial r}$")
  ax.legend(loc="upper center")
  ax.axis([-1.5,1.5,-30,2])
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [pct]")
  ax = axes[3]
  ax.plot(t,favg0[3,:]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial t_0}$")
  ax.plot(t,favg1[3,:]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial t_0}$")
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [pct]")
  ax.legend(loc="lower left")
  ax = axes[4]
  ax.plot(t,favg0[4,:]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial v}$")
  ax.plot(t,favg1[4,:]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial v}$")
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [pct]")
  ax.axis([-1.5,1.5,-1,6])
  ax.legend(loc="upper center")
  ax = axes[5]
  ax.plot(t,favg0[5,:]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial b_0}$")
  ax.plot(t,favg1[5,:]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial b_0}$")
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [pct]")
  ax.legend(loc="upper center")
  ax = axes[6]
  ax.plot(t,favg0[6,:]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial u_1}$")
  ax.plot(t,favg1[6,:]*1000.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial u_1}$")
  ax.plot(t,favg0[7,:]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial u_2}$")
  ax.plot(t,favg1[7,:]*1000.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial u_2}$")
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [ppt]")
  ax.legend(loc="upper center")
  #b = sqrt.(param[3]^2+(param[2]*(t-param[1])).^2)
  #plot(dt^2*(d2fdb2.*(param[2]*(t-param[1])./b).^2+dfdb.*param[2]^2./b)/24.0)
else
  # Now plot the results:
  fig,axes = subplots(2,3,figsize=(12,8))
  ax = axes[1]
  ax[:plot](t,favg0[1,:].+1,linestyle="--",linewidth=2,label=L"${F}$")
  ax[:plot](t,favg1[1,:].+1,linewidth=2,label=L"$\overline{{F}}$")
  ax[:plot]([-dt/2,-dt/2,-dt/2,dt/2,dt/2,dt/2],[0.9925,0.9935,0.993,0.993,0.9925,0.9935],color="b")
  ax[:annotate](L"$\Delta t$",xy=[-0.08,0.991])
  ax[:legend](loc="upper center")
  ax[:axis]([-1.5,1.5,0.988,1.002])
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux")
  ax = axes[2]
  ax[:plot](t,favg0[2,:]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial r}$")
  ax[:plot](t,favg1[2,:]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial r}$")
  ax[:legend](loc="upper center")
  ax[:axis]([-1.5,1.5,-30,2])
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [pct]")
  ax = axes[3]
  ax[:plot](t,favg0[3,:]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial t_0}$")
  ax[:plot](t,favg1[3,:]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial t_0}$")
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [pct]")
  ax[:legend](loc="lower left")
  ax = axes[4]
  ax[:plot](t,favg0[4,:]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial v}$")
  ax[:plot](t,favg1[4,:]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial v}$")
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [pct]")
  ax[:axis]([-1.5,1.5,-1,6])
  ax[:legend](loc="upper center")
  ax = axes[5]
  ax[:plot](t,favg0[5,:]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial b_0}$")
  ax[:plot](t,favg1[5,:]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial b_0}$")
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [pct]")
  ax[:legend](loc="upper center")
  ax = axes[6]
  ax[:plot](t,favg0[6,:]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial u_1}$")
  ax[:plot](t,favg1[6,:]*1000.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial u_1}$")
  ax[:plot](t,favg0[7,:]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial u_2}$")
  ax[:plot](t,favg1[7,:]*1000.,linewidth=2,label=L"$\frac{\partial\overline{F}}{\partial u_2}$")
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [ppt]")
  ax[:legend](loc="upper center")
end
savefig("integrate_transit_gradient.pdf", bbox_inches="tight")


text:
include("../src/geom_module.jl")

module btest
using BenchmarkTools
using ..geom

function test_naive_graham_scan(point_array)
    b = @benchmark geom.naive_graham_scan!($point_array)

    println("=== test_naive_graham_scan ===")
    show(
        stdout,
        MIME("text/plain"),
        b )
    println("\n\n")

end

function test_graham_scan(point_array)
    b = @benchmark geom.graham_scan!($point_array)

    println("=== test_graham_scan ===")
    show(
        stdout,
        MIME("text/plain"),
        b )
    println("\n\n")

end

function test_jarvis_march(point_array)
    b = @benchmark geom.jarvis_march!($point_array)

    println("=== test_jarvis_march ===")
    show(
        stdout,
        MIME("text/plain"),
        b )
    println("\n\n")

end

function test_jarvis_march2(point_array)
    b = @benchmark geom.jarvis2!($point_array)

    println("=== test_jarvis_march2 ===")
    show(
        stdout,
        MIME("text/plain"),
        b )
    println("\n\n")

end
function test_jarvis_march3(point_array)
    b = @benchmark geom.jarvis3!($point_array)

    println("=== test_jarvis_march3 ===")
    show(
        stdout,
        MIME("text/plain"),
        b )
    println("\n\n")

end
function test_incremental_hull(point_array)
    b = @benchmark geom.incremental_hull!($point_array)

    println("=== test_incremental_hull ===")
    show(
        stdout,
        MIME("text/plain"),
        b )
    println("\n\n")

end

function benchmark_it()
    N_point = 2_000
    # global point_array = geom.random_point(N_point)

    test_graham_scan( geom.random_point(N_point) )
    test_naive_graham_scan( geom.random_point(N_point) )
    # test_jarvis_march( geom.random_point(N_point) )
    # test_jarvis_march2( geom.random_point(N_point) )
    test_jarvis_march3( geom.random_point(N_point) )
    test_incremental_hull(geom.random_point(N_point) )

end

# using Profile

function prof_naive_graham_scan()
    N_point = 200_000
    for kk = 1:1_0
        point_array = geom.random_point(N_point)
        ch = geom.naive_graham_scan!(point_array)
    end
end

function prof_graham_scan()
    N_point = 200_000
    for kk = 1:1_0
        point_array = geom.random_point(N_point)
        ch = geom.graham_scan!(point_array)
    end
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# ## log plot
# ### semilog y

#md # ![FILE_NAME.png](images/FILE_NAME.png)

plot(exp, -5, 5, yscale=:log10, title="semilogy")


text:
module Shell

using ..Graphs: marshal, serialize

export havecommand
export minimap2, mash, mafft

havecommand(cmd) = Sys.which(cmd) !== nothing

# command line execution
function execute(cmd::Cmd; now=true)
    out = Pipe()
    err = Pipe()

    proc = run(pipeline(cmd, stdout=out, stderr=err); wait=now)

    close(out.in)
    close(err.in)

    stdout = @async String(read(out))
    stderr = @async String(read(err))

    if now
        return (
            out  = fetch(stdout),
            err  = fetch(stderr),
            code = proc.exitcode, #err,
        )
    else
        return (
            out  = stdout,
            err  = stderr,
            proc = proc,
        )
    end
end

function mash(Gs...)
    function compare(path)
        result = execute(`mash triangle $path`; now=false)
        stdout = IOBuffer(fetch(result.out))

        N     = parse(Int64,readline(stdout))
        dist  = zeros(N,N)
        names = Array{String}(undef, N)
        for (i, line) in enumerate(eachline(stdout))
            elt = split(strip(line))
            names[i] = elt[1]
            dist[i,1:(i-1)] = [parse(Float64,x) for x in elt[2:end]]
        end

        dist = dist + dist';

        return dist, names
    end

    distance, names = mktemp() do path, io
        writer = @async let
            for (i,G) ∈ enumerate(Gs)
                serialize(io, G)
            end
        end

        reader = @async compare(path)

        wait(writer)
        distance, names = fetch(reader)
    end
end

function minimap2(qry::String, ref::String)
    return execute(`minimap2 -x asm10 -m 10 -n 1 -s 30 -D -c $ref $qry`; now=false)
end

function mafft(block, case)
    out = IOBuffer()
    aln = IOBuffer()
    cmd = if case
        `mafft --auto --nuc --preservecase /dev/stdin`
    else
        `mafft --auto --nuc /dev/stdin`
    end

    names = marshal(aln, block; fmt=:fasta)
    aln = IOBuffer(String(take!(aln)))

    run(pipeline(pipeline(
        aln, cmd, out
        ); stderr=devnull ); wait=true
    )

    close(aln)

    return IOBuffer(String(take!(out))), names
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
module MPIParticleInCell

using MPI

export UniformGrid, Field, gather_field_on_root!, communicate_halo!

include("grid.jl")
include("field.jl")

end


text:
mutable struct ResettableStack{T,iip}
  data::Vector{T}
  cur::Int
  numResets::Int
  ResettableStack(ty::Type{T}) where {T} = new{T,true}(Vector{T}(),0,0)
  ResettableStack{iip}(ty::Type{T}) where {T,iip} = new{T,iip}(Vector{T}(),0,0)
end

isinplace(::ResettableStack{T,iip}) where {T,iip} = iip

isempty(S::ResettableStack) = S.cur==0
length(S::ResettableStack)  = S.cur

function push!(S::ResettableStack,x)
  if S.cur==length(S.data)
    S.cur+=1
    push!(S.data,x)
  else
    S.cur+=1
    S.data[S.cur]=x
  end
  nothing
end

safecopy(x) = copy(x)
safecopy(x::Union{Number,StaticArray}) = x
safecopy(x::Nothing) = nothing

# For DiffEqNoiseProcess S₂ fast updates
function copyat_or_push!(S::ResettableStack,x)
  if S.cur==length(S.data)
    S.cur+=1
    push!(S.data,safecopy.(x))
  else
    S.cur+=1
    curx = S.data[S.cur]
    if !isinplace(S)
      S.data[S.cur] = x
    else
      curx[2] .= x[2]
      if x[3] != nothing
        curx[3] .= x[3]
      end
      S.data[S.cur] = (x[1],curx[2],curx[3])
    end
  end
  nothing
end

function pop!(S::ResettableStack)
  S.cur-=1
  S.data[S.cur+1]
end

function iterate(S::ResettableStack, state=S.cur)
    if state == 0
        return nothing
    end

    state -= 1
    (S.data[state+1],state)
end

function reset!(S::ResettableStack,force_reset = false)
  S.numResets += 1
  S.cur = 0
  if length(S.data) > 5 && (S.numResets%FULL_RESET_COUNT==0 || force_reset)
    resize!(S.data,max(length(S.data)÷2,5))
  end
  nothing
end


text:
using Stheno, Plots, Random
using Stheno: @model

###########################  Define and inspect our model  ###########################

# Define the vanilla GP-regression generative model.
@model function model()

    # Define a smooth latent process.
    f = 1.5 * GP(EQ())

    # Define a latent noise process.
    noise = GP(Noise(1e-2))

    # Sum them to get the process of which we shall make observations.
    y = f + noise

    # Return all three processes so that we can inspect all of them.
    return f, noise, y
end

# Generate some input locations and sample from the prior.
rng, N = MersenneTwister(123546), 500;
X_prior = sort(rand(rng, N) * 10);
f, noise, y = model();

# Take a look at the mean function and kernel of y.
@show mean(y), kernel(y);

# Take a look at the cross-covariance between f and y.
@show kernel(f, y);

# Inspect the cross-covariance between noise and y.
@show kernel(noise, y);

###########################  Sample jointly from the prior  ###########################

# Sample jointly from the prior distribution over the three processes.
fX, noiseX, yX = rand(rng, [f(X_prior), noise(X_prior), y(X_prior)]);

#######################  Do posterior inference give a few samples  #######################

# Sample a few points from the prior and compute the posterior processes.
N = 15;
X = sort(rand(rng, N) * 10);
ŷ = rand(rng, y(X));
f′, noise′, y′ = (f, noise, y) | (y(X) ← ŷ);

Nplot, S = 500, 100;
Xplot = range(-2.0, stop=12.0, length=Nplot);
f′Xplot_, noise′Xplot_ = f′(Xplot), noise′(Xplot)
f′Xp, noise′Xp = rand(rng, [f′Xplot_, noise′Xplot_], S);

# Get posterior mean and marginals f′ and y′ and write them for plotting.
μ′f, σ′f = marginals(f′(Xplot));
σ′y = marginal_std(y′(Xplot));

####################################  Plot results  ####################################

plotly();

prior_plot = plot(X_prior, fX, label="f");
scatter!(prior_plot, X_prior, noiseX,
    label="noise",
    markershape=:xcross,
    markerstrokewidth=0.0,
    markersize=2);
scatter!(prior_plot, X_prior, yX,
    label="y",
    markershape=:xcross,
    markerstrokewidth=0.0,
    markersize=2);
plot!(prior_plot, X_prior, fX + noiseX - yX,
    label="f + noise - y",
    linecolor=:red);

posterior_plot = plot(Xplot, f′Xp;
    linecolor=:blue,
    linealpha=0.2,
    legend=false);
plot!(posterior_plot, Xplot, μ′f;
    linewidth=2.0,
    linecolor=:blue);
plot!(posterior_plot, Xplot, [μ′f μ′f];
    linewidth=0.0,
    fillrange=[μ′f .- 3 .* σ′y, μ′f .+ 3 * σ′y],
    fillalpha=0.3,
    fillcolor=:red);
plot!(posterior_plot, Xplot, [μ′f μ′f];
    linewidth=0.0,
    fillrange=[μ′f.- 3  .* σ′f μ′f .+ 3 .* σ′f],
    fillalpha=0.5,
    fillcolor=:blue);
scatter!(posterior_plot, X, ŷ;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    legend=false);
joint_plot = plot(prior_plot, posterior_plot, layout=(2, 1));
display(joint_plot);


text:
function generate(obj::TimeStampedGenerator, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "generate", List, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end


text:
using Glob

using Oceananigans.Fields: set!
using Oceananigans.OutputWriters: WindowedTimeAverage, checkpoint_superprefix
using Oceananigans.TimeSteppers: QuasiAdamsBashforth2TimeStepper, RungeKutta3TimeStepper, update_state!, next_time, unit_time

using Oceananigans: AbstractModel, run_diagnostic!, write_output!

import Oceananigans.OutputWriters: checkpoint_path, set!
import Oceananigans.TimeSteppers: time_step!
import Oceananigans.Utils: aligned_time_step

# Simulations are for running

#####
##### Time-step "alignment" with output and callbacks scheduled on TimeInterval
#####

function collect_scheduled_activities(sim)
    writers = values(sim.output_writers)
    callbacks = values(sim.callbacks)
    return tuple(writers..., callbacks...)
end

function schedule_aligned_Δt(sim, aligned_Δt)
    clock = sim.model.clock

    activities = collect_scheduled_activities(sim)

    for activity in activities
        aligned_Δt = aligned_time_step(activity.schedule, clock, aligned_Δt)
    end

    return aligned_Δt
end

"""
    aligned_time_step(sim, Δt)

Return a time step 'aligned' with `sim.stop_time`, output writer schedules,
and callback schedules. Alignment with `sim.stop_time` takes precedence.
"""
function aligned_time_step(sim::Simulation, Δt)
    clock = sim.model.clock

    aligned_Δt = Δt

    # Align time step with output writing and callback execution
    aligned_Δt = schedule_aligned_Δt(sim, aligned_Δt)

    # Align time step with simulation stop time
    aligned_Δt = min(aligned_Δt, unit_time(sim.stop_time - clock.time))

    # Temporary fix for https://github.com/CliMA/Oceananigans.jl/issues/1280
    aligned_Δt = aligned_Δt <= 0 ? Δt : aligned_Δt

    return aligned_Δt
end

"""
    run!(simulation; pickup=false)

Run a `simulation` until one of `simulation.stop_criteria` evaluates `true`.
The simulation will then stop.

# Picking simulations up from a checkpoint

Simulations are "picked up" from a checkpoint if `pickup` is either `true`, a `String`, or an
`Integer` greater than 0.

Picking up a simulation sets field and tendency data to the specified checkpoint,
leaving all other model properties unchanged.

Possible values for `pickup` are:

  * `pickup=true` picks a simulation up from the latest checkpoint associated with
    the `Checkpointer` in `simulation.output_writers`.

  * `pickup=iteration::Int` picks a simulation up from the checkpointed file associated
     with `iteration` and the `Checkpointer` in `simulation.output_writers`.

  * `pickup=filepath::String` picks a simulation up from checkpointer data in `filepath`.

Note that `pickup=true` and `pickup=iteration` fails if `simulation.output_writers` contains
more than one checkpointer.
"""
function run!(sim; pickup=false)

    if we_want_to_pickup(pickup)
        checkpoint_file_path = checkpoint_path(pickup, sim.output_writers)
        set!(sim.model, checkpoint_file_path)
    end

    sim.initialized = false
    sim.running = true

    while sim.running
        time_step!(sim)
    end

    return nothing
end

""" Step `sim`ulation forward by one time step. """
function time_step!(sim::Simulation)

    start_time_step = time_ns()

    if !(sim.initialized) # execute initialization step
        initialize_simulation!(sim)

        if sim.running # check that initialization didn't stop time-stepping
            @info "Executing initial time step..."

            start_time = time_ns()
            Δt = aligned_time_step(sim, sim.Δt)
            time_step!(sim.model, Δt)

            elapsed_initial_step_time = prettytime(1e-9 * (time_ns() - start_time))
            @info "    ... initial time step complete ($elapsed_initial_step_time)."
        else
            @warn "Simulation stopped during initialization."
        end
    else # business as usual...
        Δt = aligned_time_step(sim, sim.Δt)
        time_step!(sim.model, Δt)
    end

    # Callbacks and callback-like things
    [diag.schedule(sim.model)     && run_diagnostic!(diag, sim.model) for diag     in values(sim.diagnostics)]
    [callback.schedule(sim.model) && callback(sim)                    for callback in values(sim.callbacks)]
    [writer.schedule(sim.model)   && write_output!(writer, sim.model) for writer   in values(sim.output_writers)]

    end_time_step = time_ns()

    # Increment the wall clock
    sim.run_wall_time += 1e-9 * (end_time_step - start_time_step)

    return nothing
end

#####
##### Simulation initialization
#####

add_dependency!(diagnostics, output) = nothing # fallback
add_dependency!(diags, wta::WindowedTimeAverage) = wta ∈ values(diags) || push!(diags, wta)

add_dependencies!(diags, writer) = [add_dependency!(diags, out) for out in values(writer.outputs)]
add_dependencies!(sim, ::Checkpointer) = nothing # Checkpointer does not have "outputs"

we_want_to_pickup(pickup::Bool) = pickup
we_want_to_pickup(pickup::Integer) = true
we_want_to_pickup(pickup::String) = true
we_want_to_pickup(pickup) = throw(ArgumentError("Cannot run! with pickup=$pickup"))

"""
    initialize_simulation!(sim, pickup=false)

Initialize a simulation:

- Update the auxiliary state of the simulation (filling halo regions, computing auxiliary fields)
- Evaluate all diagnostics, callbacks, and output writers if sim.model.clock.iteration == 0
- Add diagnostics that "depend" on output writers

"""
function initialize_simulation!(sim)
    @info "Initializing simulation..."
    start_time = time_ns()

    model = sim.model
    clock = model.clock

    update_state!(model)

    # Output and diagnostics initialization
    [add_dependencies!(sim.diagnostics, writer) for writer in values(sim.output_writers)]

    # Evaluate all diagnostics, and then write all output at first iteration
    if clock.iteration == 0
        [run_diagnostic!(diag, model) for diag in values(sim.diagnostics)]
        [callback(sim)                for callback in values(sim.callbacks)]
        [write_output!(writer, model) for writer in values(sim.output_writers)]
    end

    sim.initialized = true

    initialization_time = prettytime(1e-9 * (time_ns() - start_time))
    @info "    ... simulation initialization complete ($initialization_time)"

    return nothing
end


text:
for f in (:log, :log2, :log10, :log1p, :exp, :exp2, :exp10, :expm1,
          :cos, :sin, :tan, :sec, :csc, :cot, :acos, :asin, :atan,
          :cosh, :sinh, :tanh, :sech, :csch, :coth, :acosh, :asinh, :atanh,
          :eint, :li2, :gamma, :lngamma, :digamma, :zeta, :erf, :erfc,
          :j0, :j1, :y0, :y1, :ai)
    fJ = Symbol(f, :!)
    @eval begin
        export $fJ
        function ($fJ)(rop::BigFloat, op::BigFloat, rnd=ROUNDING_MODE[])
            ccall(($(string("mpfr_", f)), :libmpfr), Int32,
                  (Ref{BigFloat}, Ref{BigFloat}, MPFR.MPFRRoundingMode),
                  rop, op, rnd)
            return rop
        end
    end
end

@doc """
    log!(rop, op, rnd)

Set `rop` to the natural logarithm of `op` rounded in the direction
`rnd`. Set `rop` to +0 if `op` is 1 (in all rounding modes), for
consistency with the ISOC99 and IEEE 754-2008 standards. Set `rop` to
−Inf if `op` is ±0 (i.e., the sign of the zero has no influence on the
result).
""" log!

@doc """
    log2!(rop, op, rnd)

Set `rop` to log₂ of `op` rounded in the direction `rnd`. Set `rop` to
+0 if `op` is 1 (in all rounding modes), for consistency with the
ISOC99 and IEEE 754-2008 standards. Set `rop` to −Inf if `op` is ±0
(i.e., the sign of the zero has no influence on the result).
""" log2!

@doc """
    log10!(rop, op, rnd)

Set `rop` to log\\_{10} of `op` rounded in the direction `rnd`. Set
`rop` to +0 if `op` is 1 (in all rounding modes), for consistency with
the ISOC99 and IEEE 754-2008 standards. Set `rop` to −Inf if `op` is
±0 (i.e., the sign of the zero has no influence on the result).
""" log10!

@doc """
    log1p!(rop, op, rnd)

Set `rop` to the logarithm of one plus `op`, rounded in the direction
`rnd`. Set `rop` to −Inf if `op` is −1.
""" log1p!

@doc """
    exp!(rop, op, rnd)

Set `rop` to the exponential of `op`, rounded in the direction `rnd`.
""" exp!

@doc """
    exp2!(rop, op, rnd)

Set `rop` to `2^op`, rounded in the direction `rnd`.
""" exp2!

@doc """
    exp10!(rop, op, rnd)

Set `rop` to `10^op`, rounded in the direction `rnd`.
""" exp10!

@doc """
    expm1!(rop, op, rnd)

Set `rop` to `e^op−1`, rounded in the direction `rnd`.
"""

@doc """
    cos!(rop, op, rnd)

Set `rop` to the cosine of `op`, rounded in the direction `rnd`.
""" cos!

@doc """
    sin!(rop, op, rnd)

Set `rop` to the sine of `op`, rounded in the direction `rnd`.
""" sin!

@doc """
    tan!(rop, op, rnd)

Set `rop` to the tangent of `op`, rounded in the direction `rnd`.
""" tan!

@doc """
    sec!(rop, op, rnd)

Set `rop` to the secant of `op`, rounded in the direction `rnd`.
""" sec!

@doc """
    csc!(rop, op, rnd)

Set `rop` to the cosecant of `op`, rounded in the direction `rnd`.
""" csc!

@doc """
    cot!(rop, op, rnd)

Set `rop` to the cotangent of `op`, rounded in the direction `rnd`.
""" cot!

@doc """
    acos!(rop, op, rnd)

Set `rop` to the arc-cosine of `op`, rounded in the direction
`rnd`. Note that since `acos(-1)` returns the floating-point number
closest to π according to the given rounding mode, this number might
not be in the output range `0≤rop<π` of the arc-cosine function;
still, the result lies in the image of the output range by the
rounding function.
""" acos!

@doc """
    asin!(rop, op, rnd)

Set `rop` to the arc-sine of `op`, rounded in the direction
`rnd`. Note that since `asin(1)` (resp. `asin(-1)`) returns the
floating-point number closest to π/2 (resp -π/2) according to the
given rounding mode, this number might not be in the output range
`-π/2≤rop<π/2` of the arc-sine function; still, the result lies in the
image of the output range by the rounding function.
""" asin!

@doc """
    atan!(rop, op, rnd)

Set `rop` to the arc-tangent of `op`, rounded in the direction
`rnd`. Note that since `atan(Inf)` (resp. `atan(-Inf)`) returns the
floating-point number closest to π/2 (resp -π/2) according to the
given rounding mode, this number might not be in the output range
`-π/2≤rop<π/2` of the arc-sine function; still, the result lies in the
image of the output range by the rounding function. The same holds for
`atan(op)` with large `op` and small precision of `rop`.
""" atan!

@doc """
    cosh!(rop, op, rnd)

Set `rop` to the hyperbolic cosine of `op`, rounded in the direction
`rnd`.
""" cosh!

@doc """
    sinh!(rop, op, rnd)

Set `rop` to the hyperbolic sine of `op`, rounded in the direction
`rnd`.
""" sinh!

@doc """
    tanh!(rop, op, rnd)

Set `rop` to the hyperbolic tangent of `op`, rounded in the direction
`rnd`.
""" tanh!

@doc """
    sech!(rop, op, rnd)

Set `rop` to the hyperbolic secant of `op`, rounded in the direction
`rnd`.
""" sech!

@doc """
    csch!(rop, op, rnd)

Set `rop` to the hyperbolic cosecant of `op`, rounded in the direction
`rnd`.
""" csch!

@doc """
    coth!(rop, op, rnd)

Set `rop` to the hyperbolic cotangent of `op`, rounded in the direction
`rnd`.
""" coth!

@doc """
    acosh!(rop, op, rnd)

Set `rop` to the inverse hyperbolic cosine of `op`, rounded in the
direction `rnd`.
""" acosh!

@doc """
    asinh!(rop, op, rnd)

Set `rop` to the inverse hyperbolic sine of `op`, rounded in the
direction `rnd`.
""" asinh!

@doc """
    atanh!(rop, op, rnd)

Set `rop` to the inverse hyperbolic tangent of `op`, rounded in the
direction `rnd`.
""" atanh!

@doc """
    eint!(rop, op, rnd)

Set `rop` to the exponential integral of `op`, rounded in the
direction `rnd`. This is the sum of Euler’s constant, of the logarithm
of the absolute value of `op`, and of the sum for k from 1 to infinity
of op^k/(k·k!). For positive `op`, it corresponds to the Ei function
at `op` (see formula 5.1.10 from the Handbook of Mathematical
Functions from Abramowitz and Stegun), and for negative `op`, to the
opposite of the E1 function (sometimes called eint1) at `−op` (formula
5.1.1 from the same reference).
""" eint!

@doc """
    li2!(rop, op, rnd)

Set `rop` to real part of the dilogarithm of `op`, rounded in the
direction `rnd`. MPFR defines the dilogarithm function as
−∫_{t=0}^op log(1−t)/t dt.

""" li2!

@doc """
    gamma!(rop, op, rnd)

Set `rop` to the value of the Gamma function on `op`, rounded in the
direction `rnd`. For gamma!, when `op` is a negative integer, `rop` is
set to NaN.
""" gamma!

@doc """
    lngamma!(rop, op, rnd)

Set `rop` to the value of the logarithm of the Gamma function on `op`,
rounded in the direction `rnd`. When `op` is 1 or 2, set `rop` to +0
(in all rounding modes). When `op` is an infinity or a nonpositive
integer, set `rop` to +Inf, following the general rules on special
values. When `−2k−1<op<−2k`, `k` being a nonnegative integer, set
`rop` to NaN.
""" lngamma!

@doc """
    digamma!(rop, op, rnd)

Set `rop` to the value of the Digamma (sometimes also called Psi)
function on `op`, rounded inthe direction `rnd`. When `op` is a
negative integer, set `rop` to NaN.
""" digamma!

@doc """
    zeta!(rop, op, rnd)

Set `rop` to the value of the Riemann Zeta function on `op`, rounded
in the direction `rnd`.
""" zeta!

@doc """
    erf!(rop, op, rnd)

Set `rop` to the value of the error function on `op`, rounded in the
direction `rnd`.
""" erf!

@doc """
    erfc!(rop, op, rnd)

Set `rop` to the value of the complementary error function on `op`,
rounded in the direction `rnd`.
""" erfc!

@doc """
    j0!(rop, op, rnd)

Set `rop` to the value of the first kind Bessel function of order 0 on
`op`, rounded in the direction `rnd`. When `op` is NaN, `rop` is
always set to NaN. When `op` is ±Inf, `rop` is set to +0.
""" j0!

@doc """
    j1!(rop, op, rnd)

Set `rop` to the value of the first kind Bessel function of order 1 on
`op`, rounded in the direction `rnd`. When `op` is NaN, `rop` is
always set to NaN. When `op` is ±Inf, `rop` is set to +0. When `op` is
zero, `rop` is set to +0 or −0 depending on the sign of `op`.
""" j1!

@doc """
    y0!(rop, op, rnd)

Set `rop` to the value of the second kind Bessel function of order 0
on `op`, rounded in the direction `rnd`. When `op` is NaN or negative,
`rop` is always set to NaN. When `op` is +Inf, `rop` is set to
+0. When `op` is zero, `rop` is set to -Inf.

""" y0!

@doc """
    y1!(rop, op, rnd)

Set `rop` to the value of the second kind Bessel function of order 0
on `op`, rounded in the direction `rnd`. When `op` is NaN or negative,
`rop` is always set to NaN. When `op` is +Inf, `rop` is set to
+0. When `op` is zero, `rop` is set to -Inf.
""" y1!

@doc """
    ai!(rop, x, rnd)

Set `rop` to the value of the Airy function Ai on `x`, rounded in the
direction `rnd`.  When `x` is NaN, `rop` is always set to NaN. When
`x` is +Inf or −Inf, `rop` is +0. The current implementation is not
intended to be used with large arguments. It works with `|x|`
typically smaller than 500.
""" ai!

for f in (:sin_cos, :sinh_cosh)
    fJ = Symbol(f, :!)
    @eval begin
        export $fJ
        function ($fJ)(sop::BigFloat, cop::BigFloat, op::BigFloat, rnd=ROUNDING_MODE[])
            ccall($(string("mpfr_", f), :libmpfr), Int32,
                  (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFR.MPFRRoundingMode),
                  sop, cop, op, rnd)
            return (sop, cop)
        end
    end
end

@doc """
    sin_cos!(sop, cop, op, rnd)

Set simultaneously `sop` to the sine of `op` and `cop` to the cosine
of `op`, rounded in the direction `rnd` with the corresponding
precisions of `sop` and `cop`, which must be different variables.
""" sin_cos!

@doc """
    sinh_cosh!(sop, cop, op, rnd)

Set simultaneously `sop` to the hyperbolic sine of `op` and `cop` to
the hyperbolic cosine of `op`, rounded in the direction `rnd` with the
corresponding precisions of `sop` and `cop`, which must be different
variables.
""" sinh_cosh!

for f in (:atan2, :gamma_inc, :beta, :agm, :hypot)
    fJ = Symbol(f, :!)
    @eval begin
        export $fJ
        function ($fJ)(rop::BigFloat, op1::BigFloat, op2::BigFloat, rnd=ROUNDING_MODE[])
            ccall(($(string("mpfr_", f)), :libmpfr), Int32,
                  (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFR.MPFRRoundingMode),
                  rop, op1, op2, rnd)
            return rop
        end
    end
end

@doc """
    atan2!(rop, y, x, rnd)

Set `rop` to the arc-tangent2 of `y` and `x`, rounded in the direction
`rnd`: if x > 0, `atan2(y, x) = atan(y/x)`; if x < 0, `atan2(y, x) =
sign(y)*(π - atan(|y/x|))`, thus a number from `−π` to `π`. As for
`atan`, in case the exact mathematical result is `+π` or `−π`, its
rounded result might be outside the function output range.
""" atan2!

@doc """
    gamma_inc!(rop, op, op2, rnd)

Set `rop` to the incomplete Gamma function on `op1` and `op2`, rounded
in the direction `rnd`. (In the literature, `gamma_inc` is called
upper incomplete Gamma function, or sometimes complementary incomplete
Gamma function.) For `gamma_inc`, when `op2` is zero and `op` is a
negative integer, `rop` is set to NaN.

Note: the current implementation of `gamma_inc` is slow for large
values of `rop` or `op`, in which case some internal overflow might
also occur.
""" gamma_inc!

@doc """
    beta!(rop, op1, op2, rnd)

Set `rop` to the value of the Beta function at arguments `op1` and
`op2`. Note: the current code does not try to avoid internal overflow
or underflow, and might use a huge internal precisionin some cases.
""" beta!

export factorial!
"""
    factorial!(rop, op, rnd)

Set `rop` to the factorial of `op`, rounded in the direction `rnd`.
"""
function factorial!(rop::BigFloat, op::CulongMax, rnd=ROUNDING_MODE[])
    ccall(("mpfr_fac_ui", :libmpfr), Int32, (Ref{BigFloat}, Culong, MPFR.MPFRRoundingMode),
          rop, op, rnd)
    return rop
end

for f in (:jn, :yn)
    fJ = Symbol(f, :!)
    @eval begin
        export $fJ
        function ($fJ)(rop::BigFloat, n::ClongMax, op::BigFloat, rnd=ROUNDING_MODE[])
            ccall(($(string("mpfr_", f)), :libmpfr), Int32,
                  (Ref{BigFloat}, Clong, Ref{BigFloat}, MPFR.MPFRRoundingMode),
                  rop, n, op, rnd)
            return rop
        end
    end
end

@doc """
    jn!(rop, n, op, rnd)

Set `rop` to the value of the first kind Bessel function of order `n`
on `op`, rounded in the direction `rnd`. When `op` is NaN, `rop` is
always set to NaN. When `op` is ±Inf, `rop` is set to +0. When `op` is
zero, and `n` is not zero, `rop` is set to +0 or −0 depending on the
parity and sign of `n`, and the sign of `op`.
""" jn!

@doc """
    yn!(rop, n, op, rnd)

Set `rop` to the value of the second kind Bessel function of order `n`
on `op`, rounded in the direction `rnd`. When `op` is NaN or negative,
`rop` is always set to NaN. When `op` is +Inf, `rop` is set to
+0. When `op` is zero, `rop` is set to +Inf or −Inf depending on the
parity and sign of `n`.
""" yn!

@doc """
    agm!(rop, op1, op2, rnd)

Set `rop` to the arithmetic-geometric mean of `op1` and `op2`, rounded
in the direction `rnd`. The arithmetic-geometric mean is the common
limit of the sequences `u_n` and `v_n`, where `u0=op1`, `v0=op2`,
u_{n+1} is the arithmetic mean of `u_n` and `v_n`, and `v_{n+1}` is
the geometric mean of `u_n` and `v_n`. If any operand is negative and
the other one is not zero, set `rop` to NaN. If any operand is zero
and the other one is finite (resp. infinite), set `rop` to+0
(resp. NaN).
"""

@doc """
    hypot!(rop, op1, op2, rnd)

"""

for f in (:fma, :fms,)
    fJ = Symbol(f, :!)
    @eval begin
        export $fJ
        function ($fJ)(rop::BigFloat, op1::BigFloat, op2::BigFloat, op3::BigFloat,
                       rnd=ROUNDING_MODE[])
            ccall(($(string("mpfr_", f)), :libmpfr), Int32,
                  (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat},
                   MPFR.MPFRRoundingMode),
                  rop, op1, op2, op3, rnd)
            return rop
        end
    end
end

@doc """
    fma(rop, op1, op2, op3, rnd)

Set `rop` to `(op1×op2)+op3` rounded in the direction
`rnd`. Concerning special values (signed zeros, infinities, NaN),
these functions behave like a multiplication followed by a separate
addition or subtraction. That is, the fused operation matters only for
rounding.
""" fma!

@doc """
    fms(rop, op1, op2, op3, rnd)

Set `rop` to `(op1×op2)-op3` rounded in the direction
`rnd`. Concerning special values (signed zeros, infinities, NaN),
these functions behave like a multiplication followed by a separate
addition or subtraction. That is, the fused operation matters only for
rounding.
""" fms!

for f in (:fmma, :fmms,)
    fJ = Symbol(f, :!)
    @eval begin
        export $fJ
        function ($fJ)(rop::BigFloat, op1::BigFloat, op2::BigFloat, op3::BigFloat,
                       op4::BigFloat, rnd=ROUNDING_MODE[])
            ccall(($(string("mpfr_", f)), :libmpfr), Int32,
                  (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat},
                   Ref{BigFloat}, MPFR.MPFRRoundingMode),
                  rop, op1, op2, op3, op4, rnd)
            return rop
        end
    end
end


text:
"""
Test costs for paths in an nxn grid lattice. Compare analysis with theory
"""

using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using JLD
using Parameters

function asymptotic_costs(gridsize::Int64)
    N = 1000
    G = GridNetwork(gridsize, gridsize)
    T = QuNet.TemporalGraph(G, 5, memory_costs = unit_costvector())
    p, dum1, dum2, dum3 = net_performance(T, N, 1, max_paths=4)
    return p["loss"], p["Z"]
end

depth = 50
num_trials = 5
num_pairs = 40
grid_size = 10

println("Collecting for time depth $depth")
G = GridNetwork(grid_size, grid_size)
# Create a Temporal Graph from G with timedepth i
T = QuNet.TemporalGraph(G, depth, memory_costs = zero_costvector())
# Get random pairs of asynchronus nodes
user_pairs = make_user_pairs(T, num_pairs)
# Get data
p, p_e, pat, pat_e = net_performance(T, num_trials, num_pairs, max_paths=4)
println(p)
println(p_e)
println(pat)
println(pat_e)

println("Getting asymptotic cost")
result = asymptotic_costs(grid_size)
println(result)

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet