text:
__precompile__()

module WCS

export WCSTransform,
       pix_to_world, pix_to_world!,
       world_to_pix, world_to_pix!

import Base: convert, copy, deepcopy, getindex, show, setindex!

using Compat
import Compat.ASCIIString

if isdefined(Base, :Threads)
    using Base.Threads
    const wcs_lock = SpinLock()
else
    # Before Julia 0.5 there are no threads
    const wcs_lock = 1
    lock(l) = ()
    unlock(l) = ()
end

include("../deps/deps.jl")

# -----------------------------------------------------------------------------
# Utilities

function unsafe_store_vec!{T}(p::Ptr{T}, v::Vector{T})
    for i = 1:length(v)
        unsafe_store!(p, v[i], i)
    end
end

# convert a string to a tuple of bytes
function convert_string{N}(::Type{NTuple{N,UInt8}}, s::Compat.ASCIIString)
    @assert length(s) < N
    v = zeros(UInt8, N)  # intermediate array that we can fill
    copy!(v, convert(Vector{UInt8}, s))
    (v...)
end

# load an String from a tuple of bytes, truncating at first NULL
function convert_string{N}(::Type{Compat.ASCIIString}, v::NTuple{N, UInt8})
    len = N

    # reduce length if we find a null
    i = 1
    for i=1:N
        if v[i] == 0x00
            len = i-1
            break
        end
    end
    s = Array(UInt8, len)
    copy!(s, 1, v, 1, len)
    return Compat.ASCIIString(s)  # wraps the array `s`
end

# load a String from a pointer, truncating at first NULL or maxlen
function convert_string(::Type{Compat.ASCIIString}, ptr::Ptr{UInt8},
                        maxlen::Int)
    len = maxlen

    # reduce length if we find a null
    i = 1
    for i=1:maxlen
        if unsafe_load(ptr, i) == 0x00
            len = i-1
            break
        end
    end
    return unsafe_string(ptr, len)
end

macro check_type(k, v, t)
    :(typeof($(esc(v))) <: $t || error($(esc(k))," must have type ",$t))
end

macro check_prop(k, prop, v, op, n)
    en = esc(n)
    :(($op)($prop($(esc(v))), $en) ||
      error($(esc(k))," must have ",$prop," ",$op," ",$en))
end

macro same_size(a, b)
    :(size($(esc(a))) == size($(esc(b))) ||
      error($(string(a))," must have the same dimensions as ",$(string(b))))
end

# -----------------------------------------------------------------------------
# error messaging

function get_error_message(i::Cint)
    msgptrs = cglobal((:wcs_errmsg, libwcs), Ptr{UInt8})
    msgptr = unsafe_load(msgptrs, i+1)
    unsafe_string(msgptr)
end

function assert_ok(i::Cint)
    if i != 0
        error(get_error_message(i))
    end
end

function wcslib_version()
    vers = Array(Cint, 3)
    ccall((:wcslib_version, libwcs), Ptr{UInt8}, (Ptr{Cint},), vers)
    return VersionNumber(vers[1], vers[2], vers[3])
end

# -----------------------------------------------------------------------------
# WCSTransform

# mirror of `struct pvcard`
typealias PVCard Tuple{Cint, Cint, Cdouble}  # i, m, value entries

# mirror of `struct pscard`
typealias PSCard Tuple{Cint, Cint, NTuple{72, UInt8}}  # i, m, value entries

# mirror of `struct wcserr`
immutable WCSErr
    status::Cint
    line_no::Cint
    _function::Ptr{UInt8}
    file::Ptr{UInt8}
    msg::NTuple{160, UInt8}
end

immutable linprm
    flag::Cint
    naxis::Cint
    crpix::Ptr{Cdouble}
    pc::Ptr{Cdouble}
    cdelt::Ptr{Cdouble}
    dispre::Ptr{Void}  # Ptr{disprm}
    disseq::Ptr{Void}  # Ptr{disprm}
    piximg::Ptr{Cdouble}
    imgpix::Ptr{Cdouble}
    i_naxis::Cint
    unity::Cint
    affine::Cint
    simple::Cint
    err::Ptr{WCSErr}
    tmpcrd::Ptr{Cdouble}
    m_flag::Cint
    m_naxis::Cint
    m_crpix::Ptr{Cdouble}
    m_pc::Ptr{Cdouble}
    m_cdelt::Ptr{Cdouble}
    m_dispre::Ptr{Void}  # Ptr{disprm}
    m_disseq::Ptr{Void}  # Ptr{disprm}
end

immutable prjprm
    flag::Cint
    code::NTuple{4, UInt8}
    r0::Cdouble
    pv::NTuple{30, Cdouble}
    phi0::Cdouble
    theta0::Cdouble
    bounds::Cint
    name::NTuple{40, UInt8}
    category::Cint
    pvrange::Cint
    simplezen::Cint
    equiareal::Cint
    conformal::Cint
    _global::Cint
    divergent::Cint
    x0::Cdouble
    y0::Cdouble
    err::Ptr{WCSErr}
    padding::Ptr{Void}
    w::NTuple{10, Cdouble}
    m::Cint
    n::Cint
    prjx2s::Ptr{Void}
    prjs2x::Ptr{Void}
end

immutable celprm
    flag::Cint
    offset::Cint
    phi0::Cdouble
    theta0::Cdouble
    ref::NTuple{4, Cdouble}
    prj::prjprm
    euler::NTuple{5, Cdouble}
    latpreq::Cint
    isolat::Cint
    err::Ptr{WCSErr}
    padding::Ptr{Void}
end

immutable spcprm
    flag::Cint
    _type::NTuple{8, UInt8}
    code::NTuple{4, UInt8}
    crval::Cdouble
    restfrq::Cdouble
    restwav::Cdouble
    pv::NTuple{7, Cdouble}
    w::NTuple{6, Cdouble}
    isGrism::Cint
    padding1::Cint
    err::Ptr{WCSErr}
    padding2::Ptr{Void}
    spxX2P::Ptr{Void}
    spxP2S::Ptr{Void}
    spxS2P::Ptr{Void}
    spxP2X::Ptr{Void}
end

# mirror of `wcsprm` struct in wcslib
"""
WCSTransform(naxis; kvs...)

Construct a WCS transformation with the given number of axes `naxis`.
Keyword arguments can be passed to set various attributes of the
transform. Specifying keyword arguments is equivalent to setting them
after construction:

```
julia> wcs = WCSTransform(2; crpix=[1000., 1000.])
```

is equilvalent to:

```
julia> wcs = WCSTransform(2)

julia> wcs[:crpix] = [1000., 1000.]
```
"""
type WCSTransform
    flag::Cint
    naxis::Cint
    crpix::Ptr{Cdouble}
    pc::Ptr{Cdouble}
    cdelt::Ptr{Cdouble}
    crval::Ptr{Cdouble}
    cunit::Ptr{Void}
    ctype::Ptr{Void}
    lonpole::Cdouble
    latpole::Cdouble
    restfrq::Cdouble
    restwav::Cdouble
    npv::Cint
    npvmax::Cint
    pv::Ptr{PVCard}
    nps::Cint
    npsmax::Cint
    ps::Ptr{PSCard}
    cd::Ptr{Cdouble}
    crota::Ptr{Cdouble}
    altlin::Cint
    velref::Cint
    alt::NTuple{4, UInt8}
    colnum::Cint
    colax::Ptr{Cint}
    cname::Ptr{Void}
    crder::Ptr{Cdouble}
    csyer::Ptr{Cdouble}
    dateavg::NTuple{72, UInt8}
    dateobs::NTuple{72, UInt8}
    equinox::Cdouble
    mjdavg::Cdouble
    mjdobs::Cdouble
    obsgeo::NTuple{3, Cdouble}
    radesys::NTuple{72, UInt8}
    specsys::NTuple{72, UInt8}
    ssysobs::NTuple{72, UInt8}
    velosys::Cdouble
    zsource::Cdouble
    ssyssrc::NTuple{72, UInt8}
    velangl::Cdouble
    wcsname::NTuple{72, UInt8}
    ntab::Cint
    nwtb::Cint
    tab::Ptr{Void}  # Ptr{tabprm}
    wtb::Ptr{Void}  # Ptr{wtbarr}
    lngtyp::NTuple{8, UInt8}
    lattyp::NTuple{8, UInt8}
    lng::Cint
    lat::Cint
    spec::Cint
    cubeface::Cint
    types::Ptr{Cint}
    padding::Ptr{Void}
    lin::linprm
    cel::celprm
    spc::spcprm
    err::Ptr{WCSErr}
    m_padding::Ptr{Void}
    m_flag::Cint
    m_naxis::Cint
    m_crpix::Ptr{Cdouble}
    m_pc::Ptr{Cdouble}
    m_cdelt::Ptr{Cdouble}
    m_crval::Ptr{Cdouble}
    m_cunit::Ptr{Void}
    m_ctype::Ptr{Void}
    m_pv::Ptr{PVCard}
    m_ps::Ptr{PSCard}
    m_cd::Ptr{Cdouble}
    m_crota::Ptr{Cdouble}
    m_colax::Ptr{Cint}
    m_cname::Ptr{Void}
    m_crder::Ptr{Cdouble}
    m_csyer::Ptr{Cdouble}
    m_tab::Ptr{Void}  # Ptr{tabprm}
    m_wtb::Ptr{Void}  # Ptr{wtbarr}

    function WCSTransform(naxis::Integer; kvs...)
        w = new(-1)
        status = ccall((:wcsini, libwcs), Cint,
                       (Cint, Cint, Ref{WCSTransform}),
                       1, naxis, w)
        assert_ok(status)
        finalizer(w, free!)
        for (k, v) in kvs
            w[k] = v
        end

        # wcsset() modifies the WCSTransform, so call it here so it doesn't
        # get called in wcss2p and wcsp2s. Otherwise, calling wcss2p / wcsp2s
        # on the same WCSTransform object from multiple threads would create
        # a race condition.
        status = ccall((:wcsset, libwcs), Cint, (Ref{WCSTransform},), w)
        assert_ok(status)

        return w
    end
end

function free!(w::WCSTransform)
    status = ccall((:wcsfree, libwcs), Cint, (Ref{WCSTransform},), w)
    assert_ok(status)
end

# -----------------------------------------------------------------------------
# getting attributes of a WCSTransform
# These return newly-allocated memory (not views of the WCSTransform).

function getindex(wcs::WCSTransform, k::Symbol)
    @assert wcs.flag != -1
    wcs.flag = 0
    naxis = wcs.naxis

    # double[naxis]
    if k in (:cdelt, :crder, :crota, :crpix, :crval, :csyer)
        v = Array(Float64, naxis)
        unsafe_copy!(pointer(v), getfield(wcs,k), naxis)

    # char[72,naxis]
    elseif k in (:cname, :ctype, :cunit)
        p = convert(Ptr{UInt8}, getfield(wcs, k))
        v = Array(Compat.ASCIIString, naxis)
        for i=1:naxis
            pi = p + 72*(i-1)  # Ptr{UInt8} to the i-th entry.
            v[i] = convert_string(Compat.ASCIIString, pi, 72)
        end

    # PVCard[]
    elseif k === :pv
        error("pv getter not yet implemented")

    # PSCard[]
    elseif k === :ps
        error("ps getter not yet implemented")

    # double[naxis,naxis]
    elseif k in (:cd, :pc)
        v = Array(Cdouble, naxis, naxis)
        unsafe_copy!(pointer(v), getfield(wcs,k), naxis*naxis)

    # double
    elseif k in (:equinox,:latpole,:lonpole,:mjdavg,:mjdobs,
                 :restfrq,:restwav,:velangl,:velosys,:zsource)
        v = getfield(wcs, k)

    # int
    elseif k === :colnum
        v = Int(getfield(wcs, k))

    # char[72]
    elseif k in (:dateavg,:dateobs,:radesys,:specsys,:ssysobs,:ssyssrc,
                 :wcsname)
        v = convert_string(Compat.ASCIIString, getfield(wcs, k))

    # double[3]
    elseif k === :obsgeo
        v = getfield(wcs, k)  # Tuple{Cdouble, Cdouble, Cdouble}

    # char[4], but only uses first
    elseif k === :alt
        v = Char(wcs.alt[1])

    else
        error("unrecognized keyword argument \"$k\"")
    end

    return v
end

# -----------------------------------------------------------------------------
# modifying a WCSTransform

function setindex!(wcs::WCSTransform, v, k::Symbol)
    @assert wcs.flag != -1
    wcs.flag = 0
    naxis = wcs.naxis

    # double[naxis]
    if k in (:cdelt, :crder, :crota, :crpix, :crval, :csyer)
        @check_type k v Vector{Float64}
        @check_prop k length v (==) naxis
        unsafe_copy!(getfield(wcs,k), pointer(v), naxis)

    # char[72,naxis]
    elseif k in (:cname, :ctype, :cunit)
        @check_type k v Vector{Compat.ASCIIString}
        @check_prop k length v (==) naxis

        p = convert(Ptr{UInt8}, getfield(wcs, k))
        for i in 1:naxis
            pi = p + 72*(i-1)  # Ptr{UInt8} to i-th entry.
            n = length(v[i])
            @assert n < 72
            unsafe_copy!(pi, pointer(v[i]), n)
            unsafe_store!(pi, 0x00, n+1)  # trailing null
        end

    # PVCard[]
    elseif k === :pv
        @check_prop k length v (<=) wcs.npvmax
        npv = length(v)
        for i in 1:npv  # TODO: what if a convert fails halfway through?
            unsafe_store!(wcs.pv, convert(PVCard, v[i]), i)
        end
        wcs.npv = npv

    # PSCard[]
    elseif k === :ps
        @check_prop k length v (<=) wcs.npsmax
        nps = length(v)
        for i in 1:nps
            unsafe_store!(wcs.ps, convert(PSCard, v[i]), i)
        end
        wcs.nps = nps

    # double[naxis,naxis]
    elseif k in (:cd, :pc)
        @check_type k v Matrix{Float64}
        @check_prop k size v (==) (naxis,naxis)
        unsafe_store_vec!(getfield(wcs,k), vec(v'))

    # double
    elseif k in (:equinox,:latpole,:lonpole,:mjdavg,:mjdobs,
                 :restfrq,:restwav,:velangl,:velosys,:zsource)
        setfield!(wcs, k, convert(Float64,v))

    # int
    elseif k === :colnum
        setfield!(wcs, k, convert(Cint, v))

    # char[72]
    elseif k in (:dateavg,:dateobs,:radesys,:specsys,:ssysobs,:ssyssrc,
                 :wcsname)
        setfield!(wcs, k, convert_string(NTuple{72, UInt8}, v))

    # double[3]
    elseif k === :obsgeo
        @check_type k v Vector{Float64}
        @check_prop k length v (==) 3
        wcs.obsgeo = (v[1], v[2], v[3])

    # char[4], but only uses first
    elseif k === :alt
        @check_type k v Char
        (('A' <= v <= 'Z') || v == ' ') || error("alt must be 'A'-'Z' or ' '")
        wcs.alt = (UInt8(v), 0x00, 0x00, 0x00)

    else
        error("unrecognized keyword argument \"$k\"")
    end
end

function deepcopy(wcs::WCSTransform)
    dst = WCSTransform(wcs.naxis)
    p = convert(Ptr{Cint}, C_NULL)
    ccall((:wcssub, libwcs), Cint,
          (Cint, Ref{WCSTransform}, Ptr{Cint}, Ptr{Cint}, Ref{WCSTransform}),
          0, wcs, p, p, dst)
    return dst
end

# The default Julia `copy` would be a gotcha, as it would leave two
# WCSTransforms linked by their underlying C-allocated arrays;
# only allow deepcopy.
copy(w::WCSTransform) = deepcopy(w)

# TODO: more info here.
function show(io::IO, wcs::WCSTransform)
    print(io, "WCSTransform(naxis=$(wcs.naxis))")
end

# -----------------------------------------------------------------------------
# transforms

"""
pix_to_world(wcs, pixcoords)

Convert the array of pixel coordinates `pixcoords` to world coordinates
according to the WCSTransform `wcs`. `pixcoords` should be a 2-d array
where \"pixcoords[:, i]\" is the i-th set of coordinates, or a 1-d array
representing a single set of coordinates.

The return value is the same shape as `pixcoords`.
"""
pix_to_world(wcs::WCSTransform, pixcoords::VecOrMat{Float64}) =
    pix_to_world!(wcs, pixcoords, similar(pixcoords))

"""
pix_to_world!(wcs, pixcoords, worldcoords[; stat=, imcoords=, phi=, theta=])

Convert the array of pixel coordinates `pixcoords` to world coordinates
according to the WCSTransform `wcs`, storing the result in the
`worldcoords` and `stat` arrays. `pixcoords` should be a 2-d array where
\"pixcoords[:, i]\" is the i-th set of coordinates, or a 1-d array
representing a single set of coordinates. `worldcoords` must be
the same size and type as `pixcoords`.

If given, the arrays `stat`, `imcoords`, `phi`, `theta` will be used
to store intermediate results. Their sizes and types must all match
`pixcoords`, except for `stat` which should be the same size but of type
Cint (typically Int32).
"""
function pix_to_world!(wcs::WCSTransform, pixcoords::VecOrMat{Float64},
                       worldcoords::VecOrMat{Float64};
                       stat=similar(pixcoords, Cint),
                       imcoords=similar(pixcoords),
                       phi=similar(pixcoords),
                       theta=similar(pixcoords))
    nelem = size(pixcoords, 1)
    ncoords = size(pixcoords, 2)
    if nelem < wcs.naxis
        error("size(pixcoords, 1) must be greater than or equal to naxis")
    end
    @same_size worldcoords pixcoords
    @same_size imcoords pixcoords
    @same_size phi pixcoords
    @same_size theta pixcoords
    @same_size stat pixcoords
    ccall((:wcsp2s, libwcs), Cint,
          (Ref{WCSTransform}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}),
          wcs, ncoords, nelem, pixcoords, imcoords, phi, theta, worldcoords,
          stat)
    return worldcoords
end

"""
world_to_pix(wcs, worldcoords)

Convert the array of world coordinates `worldcoords` to pixel coordinates
according to the WCSTransform `wcs`. `worldcoords` is a 2-d array
where \"worldcoords[:, i]\" is the i-th set of coordinates, or a 1-d array
representing a single set of coordinates.

The return value is the same size as `worldcoords`.
"""
world_to_pix(wcs::WCSTransform, worldcoords::VecOrMat{Float64}) =
    world_to_pix!(wcs, worldcoords, similar(worldcoords))

"""
world_to_pix!(wcs, worldcoords, pixcoords[; stat=, phi=, theta=, imcoords=])

Convert the array of pixel coordinates `worldcoords` to pixel coordinates
according to the WCSTransform `wcs`, storing the result in the
`pixcoords` array. `worldcoords` should be a 2-d array where
\"worldcoords[:, i]\" is the i-th set of coordinates, or a 1-d array
representing a single set of coordinates. `pixcoords` must be
the same size and type as `worldcoords`.

If given, the arrays `stat`, `imcoords`, `phi`, `theta` will be used
to store intermediate results. Their sizes and types must all match
`worldcoords`, except for `stat` which should be the same size but of type
Cint (typically Int32).
"""
function world_to_pix!(wcs::WCSTransform, worldcoords::VecOrMat{Float64},
                       pixcoords::VecOrMat{Float64};
                       stat=similar(pixcoords, Cint),
                       phi=similar(pixcoords),
                       theta=similar(pixcoords),
                       imcoords=similar(pixcoords))
    nelem = size(worldcoords, 1)
    ncoords = size(worldcoords, 2)
    if nelem < wcs.naxis
        error("size(worldcoords, 1) must be greater than or equal to naxis")
    end
    @same_size pixcoords worldcoords
    @same_size imcoords worldcoords
    @same_size phi worldcoords
    @same_size theta worldcoords
    @same_size stat worldcoords
    ccall((:wcss2p, libwcs), Cint,
          (Ref{WCSTransform}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}),
          wcs, ncoords, nelem, worldcoords, phi, theta, imcoords, pixcoords,
          stat)
    return pixcoords
end

# -----------------------------------------------------------------------------
# WCSTransform <--> header

"""
from_header(header[; relax=0, ctrl=0, ignore_rejected=false, table=false])

Parse the FITS image header in the String `header`, returning a
`Vector{WCSTransform}` giving all the transforms defined in the header.
"""
function from_header(header::Compat.ASCIIString; relax::Integer=0, ctrl::Integer=0,
                     ignore_rejected::Bool=false, table::Bool=false)
    @assert ctrl >= 0  # < 0 modifies the header
    nkeyrec::Integer=div(length(header), 80)
    nreject = Ref{Cint}(0)
    nwcs = Ref{Cint}(0)
    wcsptr = Ref{Ptr{WCSTransform}}(0)
    keysel = 0
    status = Cint(0)

    # wcsbth & wcspih are not thread-safe; see
    # http://www.atnf.csiro.au/people/mcalabre/WCS/wcslib/threads.html
    lock(wcs_lock)
    if table
        colsel = convert(Ptr{Cint}, C_NULL)
        status = ccall((:wcsbth, libwcs), Cint,
                       (Ptr{UInt8}, Cint, Cint, Cint, Cint, Ptr{Cint},
                        Ref{Cint}, Ref{Cint}, Ref{Ptr{WCSTransform}}),
                       header, nkeyrec, relax, ctrl, keysel, colsel,
                       nreject, nwcs, wcsptr)
    else
        status = ccall((:wcspih, libwcs), Cint,
                       (Ptr{UInt8}, Cint, Cint, Cint, Ref{Cint}, Ref{Cint},
                        Ref{Ptr{WCSTransform}}),
                       header, nkeyrec, relax, ctrl, nreject, nwcs, wcsptr)
    end
    unlock(wcs_lock)
    assert_ok(status)
    p = wcsptr[]
    result = WCSTransform[unsafe_load(p, i) for i = 1:nwcs[]]

    # Free the array of pointers allocated by wcslib. (But not the
    # structs they point to, which we've `unsafe_load`ed into a
    # Julia-allocated array!)
    # TODO: this could call the wrong `free` function if wcslib is linked
    # against a different libc than Julia.
    Libc.free(p)

    # For each of the WCSTransforms, register a finalizer and finish
    # initialization of the struct by calling wcsset. This avoids race
    # conditions between threads using the same WCSTransform.
    for w in result
        finalizer(w, free!)
        status = ccall((:wcsset, libwcs), Cint, (Ref{WCSTransform},), w)
        assert_ok(status)
    end

    if !ignore_rejected && nreject[] != 0
        error("$(nreject[]) WCS transformations were rejected; " *
              "use ignore_rejected = true keyword to ignore")
    end
    return result
end

"""
to_header(wcs[; relax=0])

Write the WCSTransform `wcs` as a FITS header.
"""
function to_header(wcs::WCSTransform; relax::Integer=0)
    nkeyrec = Ref{Cint}(0)
    hdrptr = Ref{Ptr{UInt8}}(C_NULL)
    status = ccall((:wcshdo, libwcs), Cint,
                   (Cint, Ref{WCSTransform}, Ref{Cint}, Ref{Ptr{UInt8}}),
                   relax, wcs, nkeyrec, hdrptr)
    assert_ok(status)
    header = unsafe_string(hdrptr[])
    Libc.free(hdrptr[])
    return header
end

end  # module


text:
"""
A Stochastic Gradient optimisation algorithm based on relativistic Hamiltonian dynamics

X. Lu, V. Perrone, L. Hasenclever, Y. W. Teh, Sebastian J. Vollmer relativistic Monte Carlo https://arxiv.org/abs/1609.04388

Fields:
    - x: position
    - p: momentum
    - ntiers: number of iteration
    - mass: diagonal mass matrix
    - c: speed of light
    - D: diagonal preconditioner
    - independent_momenta: true means kinetic energy is applied to each component of momenta
                           false mean kinetic energy is applied to euclidean length of momenta
Constructors:
    - construct using `RelFrozenState(x::Array{Float64};stepsize = 0.001, p=:none, mass=[1.0],c=[1.0], D=[1.0],  independent_momenta=true)`
"""

type RelFrozenState <: SamplerState
    x::Array{Float64}
    p::Array{Float64}

    stepsize::Float64
    mass::Array{Float64,1}
    c::Array{Float64,1}
    D::Array{Float64,1}

    independent_momenta::Bool
    function RelFrozenState(x::Array{Float64};stepsize = 0.001, p=:none, mass=[1.0],c=[1.0], D=[1.0],  independent_momenta=true)
        if isa(mass,Number)
          mass = mass * ones(length(x))
        end
        if p == :none
            p = sample_rel_p(mass, c, length(x))
        end
    new(x,p,stepsize,mass,c,D,independent_momenta)
    end
end

function sample!(s::RelFrozenState, sgrad)
  D = s.D
  stepsize = s.stepsize

  m = s.mass
  c = s.c
  independent_momenta=s.independent_momenta
  independent_momenta ? tmp = m .* sqrt(s.p.^2 ./ (m.^2 .* c.^2) + 1) : tmp = m .* sqrt(s.p'*s.p ./ (m.^2 .* c.^2) + 1)
  p_grad = s.p./tmp
  s.p[:] += stepsize.*(sgrad(s.x)-D[1]*p_grad)
  independent_momenta ? tmp = m .* sqrt(s.p.^2 ./ (m.^2 .* c.^2) + 1)  :  tmp = m .* sqrt(s.p'*s.p ./ (m.^2 .* c.^2) + 1)
  s.x[:] += stepsize.*s.p./tmp
end

ext: jl


text:
using BinaryBuilder

# Collection of sources required to build SAASGP4Builder
sources = [
    "https://github.com/benelsen/SAASGP4Builder.jl/raw/master/src/SGP4_small_V7.8_LINUX64.tar.gz" => "36bf3c1df0d5ed5dae59ac464548c15a888fcda256d93950a95b8bf39cc42f1a",
]

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

mkdir $prefix/lib

cp Lib/SGP4_Open_License.txt $prefix/lib/

# cp Lib/lib_external/* $prefix/lib/
# cp Lib/lib_external/libdl.so.2 $prefix/lib/
# cp Lib/lib_external/ld-linux-x86-64.so.2 $prefix/lib/
# cp Lib/lib_external/libpthread.so.0 $prefix/lib/

cp Lib/libifcore.so.5 $prefix/lib/
cp Lib/libifcoremt.so.5 $prefix/lib/
cp Lib/libifport.so.5 $prefix/lib/
cp Lib/libimf.so $prefix/lib/
cp Lib/libintlc.so.5 $prefix/lib/
# cp Lib/libiomp5.so $prefix/lib/
cp Lib/libsvml.so $prefix/lib/

mv Lib/DllMain.dll $prefix/lib/DllMain.so
mv Lib/EnvConst.dll $prefix/lib/EnvConst.so
mv Lib/AstroFunc.dll $prefix/lib/AstroFunc.so
mv Lib/TimeFunc.dll $prefix/lib/TimeFunc.so
mv Lib/Tle.dll $prefix/lib/Tle.so
mv Lib/Sgp4Prop.dll $prefix/lib/Sgp4Prop.so
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "DllMain", :dllmain),
    LibraryProduct(prefix, "EnvConst", :envconst),
    LibraryProduct(prefix, "AstroFunc", :astrofunc),
    LibraryProduct(prefix, "TimeFunc", :timefunc),
    LibraryProduct(prefix, "Tle", :tle),
    LibraryProduct(prefix, "Sgp4Prop", :sgp4prop),
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "SAASGP4", sources, script, platforms, products, dependencies)


text:
# Leaf node: contains a small hash table of size O(n^(1/7))
mutable struct Leaf{K,V}
    # Bitmap of which slots are occupied, if length(entries) <= 63.  MSB is used for
    # the mutable flag.
    slots::UInt64

    # Used if length(entries) > 63.
    big_slots::BitVector
    entries::Vector{Tuple{K,V}}

    function Leaf{K,V}(capacity::Int) where {K,V}
        entries = Vector{Tuple{K,V}}(undef, capacity)
        # MSB of slots is the 'is_mutable' bit
        slots = 0x8000000000000000
        if capacity <= 63
            big_slots = _empty_big_slots()
        else
            big_slots = BitVector(undef, capacity)
            for i=1:capacity
                @inbounds big_slots[i] = false
            end
        end
        return new(slots, big_slots, entries)
    end

    function Leaf{K,V}(slots, big_slots, entries) where {K,V}
        return new(slots, big_slots, entries)
    end
end

# Maximum number of probes we will perform in a leaf hash table.  max_probe must be a
# constant for asymptotics. The smaller the constant, the more wasted space, because leaf
# enlargement is triggered by exceeding max_probes on an insert. 16 is a reasonable number
# of probes. Because of the xor probing you get reasonable cache locality, so 16 probes
# doesn't translate into 16 cache misses, e.g. for 16-byte (k,v) pairs you're looking at
# maybe 4-5 cacheline misses as the max. But usually 1-2 cachelines for a probe sequence.
const max_probe = 16

# Sequence of hash table sizes for the leaves, where the first few table sizes are
# 1,2,3,4.. to conserve space, and each successive size is usually not more than 5/4 the
# value of the previous. The exponential sequence of sizes (increase factor ~ 5/4)
# gives Θ(1) amortization for inserts, similar to how a typical dynamic Vector doubles in
# size to give amortized Θ(1) insert cost. Using a 5/4 increase ensures that we will
# usually have no more than 20% empty slots.
const _leaf_table_sizes = [1,2,3,4,5,6,8,11,13,15,19,23,27,33,41,47,59,73,89,113,127,147,163,191,233]

# Returns the appropriate table size for storing N items
function _leaf_table_size(N::Int)
    rough_capacity = convert(Int, ceil(N*11/10))
    pos = searchsortedfirst(_leaf_table_sizes, rough_capacity)
    if pos <= length(_leaf_table_sizes)
        return _leaf_table_sizes[pos]
    else
        return convert(Int, ceil(N*5/4))
    end
end

# Returns the next largest table size, given the current capacity
function _next_leaf_table_size(capacity::Int)
    pos = searchsortedfirst(_leaf_table_sizes, capacity)
    if pos < length(_leaf_table_sizes)
        return _leaf_table_sizes[pos+1]
    else
        # This should be unreachable; would need a bad hash function or
        # N > 10^16 for this to be used.
        return convert(Int, ceil(capacity*5/4))
    end
end

@inline function _slot_occupied(leaf::Leaf{K,V}, slot::UInt64, num_slots::Int) where {K,V}
    @dassert1 (slot >= 1) && (slot <= num_slots)
    if num_slots <= 63
        return (leaf.slots & (UInt64(1) << (slot-1))) != UInt64(0)
    else
        slot2 = reinterpret(Int64, slot)
        return @inbounds leaf.big_slots[slot2]
    end
end

function is_mutable(leaf::Leaf{K,V}) where {K,V}
    return (leaf.slots & 0x8000000000000000) != 0
end

function mark_immutable!(leaf::Leaf{K,V}) where {K,V}
    leaf.slots = leaf.slots & 0x7fffffffffffffff
    @dassert1 !is_mutable(leaf)
end

function _get_mutable_version(
    leaf::Leaf{K,V},
    pathmask::UInt64,
    pathhash::UInt64
) where {K,V}
    is_mutable(leaf) && return leaf
    num_slots = length(leaf.entries)

    # Filter out entries which do not belong to this path
    count=0
    for i=1:num_slots
        slot = reinterpret(UInt64, i)
        if _slot_occupied(leaf, slot, num_slots)
            (ki,vi) = @inbounds leaf.entries[i]
            h = hash(ki)
            if (h & pathmask) == pathhash
                count += 1
            end
        end
    end

    new_capacity = _leaf_table_size(count)
    new_leaf = Leaf{K,V}(new_capacity)
    for i=1:num_slots
        slot = reinterpret(UInt64, i)
        if _slot_occupied(leaf, slot, num_slots)
            (ki,vi) = @inbounds leaf.entries[i]
            h = hash(ki)
            if (h & pathmask) == pathhash
                (new_leaf, _) = _insert!(new_leaf, h, ki, vi)
            end
        end
    end
    return new_leaf
end

function _put_entry(leaf::Leaf{K,V}, index::UInt64, k::K, v::V) where {K,V}
    @dassert1 (index >= 1) && (index <= length(leaf.entries))
    @inbounds leaf.entries[index] = (k,v)
    num_slots = length(leaf.entries)
    if num_slots <= 63
        leaf.slots |= UInt64(1) << UInt64(index-1)
    else
        @inbounds leaf.big_slots[index] = true
    end
end

# Returns the ith slot in the hash table probe sequence: i=0 is the first
# slot we check, i=1 the second etc.  We use an xor probe sequence,
# where the slot is xor(hashcode,i) modulo number of slots.  This probe
# sequence is fast to compute and gives good cache-locality, because you
# check all the slots in a cacheline before moving to a new cacheline.
# TODO: consider a probing sequence that mixes the cache-locality of xor
# with an anti-clumping measure, e.g., bottom two bits from xor, top bits
# from quadratic probe sequence.
@inline function _hashindex(hashcode::UInt64, i::Int, num_slots::Int)
    ui = reinterpret(UInt64, i)
    z = xor(hashcode, ui)

    # TODO: z = 1+((z >= num_slots) ? z-num_slots : z)
    return 1+Base.checked_urem_int(z, reinterpret(UInt64, num_slots))
end

# Unique instance of an empty BitVector
@generated function _empty_big_slots()
    return BitVector()
end


text:
using Tk
w = Toplevel("GUI enabling/disabling")
fr = Frame(w)
pack(fr, {:expand=>true, :fill => "both"})

value = Entry(fr)
increment = Button(fr, "+")
decrement = Button(fr, "-")

formlayout(value, "Value:")
formlayout(increment, " ")
formlayout(decrement, " ")

## value stores a string
set_value(value, "0") ## The field is initialized to zero.
get(value::Tk_Entry) = try parseint(get_value(value)) catch e nothing end

function update()
    cur_value = get(value)
    set_enabled(value,     isa(cur_value, Integer) && cur_value == 0)
    set_enabled(increment, isa(cur_value, Integer) && cur_value < 10)
    set_enabled(decrement, isa(cur_value, Integer) && cur_value > 0)
end

crement = function(step)
  set_enabled(value, true)
  set_value(value, string(get(value) + step))
  update()
end
tk_bind(increment, "command", path -> crement(1))
tk_bind(decrement, "command", path -> crement(-1))
update()

function validate_command(path, P)
    try
        if length(P) > 0 parseint(P); update() end
        tcl("expr", "TRUE")
    catch e
        tcl("expr", "FALSE")
    end
end
function invalid_command(path, W)
    println("Invalid value")
    tcl(W, "delete", "@0", "end")
end

tk_configure(value, {:validate=>"focusout", :validatecommand=>validate_command, :invalidcommand=>invalid_command })

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using FeatherFiles
using DataValues
using IteratorInterfaceExtensions
using TableTraits
using Test

@testset "FeatherFiles" begin

source = [(Name="John", Age=34., Children=2),
    (Name="Sally", Age=54., Children=1),
    (Name="Jim", Age=34., Children=0)]

output_filename = tempname() * ".feather"

source |> save(output_filename)

try
    sink = load(output_filename) |> IteratorInterfaceExtensions.getiterator |> collect

    @test source == sink

    featherfile = load(output_filename)

    @test IteratorInterfaceExtensions.isiterable(featherfile) == true
    @test TableTraits.supports_get_columns_copy_using_missing(featherfile) == true
    ff_as_cols = TableTraits.get_columns_copy_using_missing(featherfile)
    @test ff_as_cols == (Name=["John", "Sally", "Jim"], Age=[34., 54., 34.], Children=[2,1,0])
finally
    GC.gc()
    GC.gc()
    # rm(output_filename)
end

source2 = [(Name=DataValue("John"), Age=DataValue(34.), Children=DataValue{Int}()),
    (Name=DataValue("Sally"), Age=DataValue{Float64}(), Children=DataValue(1)),
    (Name=DataValue{String}(), Age=DataValue(34.), Children=DataValue(0))]

output_filename2 = tempname() * ".feather"

source2 |> save(output_filename2)

try
    sink2 = load(output_filename2) |> IteratorInterfaceExtensions.getiterator |> collect

    @test source2 == sink2

    featherfile = load(output_filename2)
    @test IteratorInterfaceExtensions.isiterable(featherfile) == true
    @test TableTraits.supports_get_columns_copy_using_missing(featherfile) == true
    ff_as_cols = TableTraits.get_columns_copy_using_missing(featherfile)
    @test isequal(ff_as_cols, (Name=["John", "Sally", missing], Age=[34., missing, 34.], Children=[missing,1,0]))
finally
    GC.gc()
    GC.gc()
    # rm(output_filename2)
end

end


text:
"""
Return Nx2 array with segment IDs and semantic labels
"""
function load_semantic_list(fn)
    return readdlm(fn, ',', Int64)
end

function load_segmentation(fn)
    return h5read(fn, "main");
end

function save_with_compression(fn, arr)
  f = h5open(fn, "w")
  f["main", "blosc", 3] = arr
  close(f)
end

"""
Subsample a 3D array
"""
function subsample3d(arr, n=(2,2,1))
    return arr[1:n[1]:end, 1:n[2]:end, 1:n[3]:end]
end

"""

    segm_overlap!{S,T}( seg1::Array{S}, seg2::Array{T}, counts=Dict{Tuple{S,T},Int}() )

  Computes an overlap matrix between two segmentations represented by an explicit dictionary
  of tuple keys. Modifies an existing counts dictionary if passed.
"""
function segm_overlap!{S,T}( seg1::Array{S}, seg2::Array{T}, counts=Dict{Tuple{S,T},Int}() )

    zS, zT = S(0), T(0)
    for i in eachindex(seg1)

      v1, v2 = seg1[i], seg2[i]

      if v1 == zS continue end
      if v2 == zT continue end

      k = (v1,v2)
      if !haskey(counts, k) counts[k] = 1
      else counts[k] += 1
      end

    end

    counts
end

segm_overlap(seg1, seg2) = segm_overlap!(seg1, seg2)

"""

    overlap_in_chunks( seg1, seg2, chunk_shape, verb=true )

  Computes an overlap matrix between two segmentations represented by an explicit dictionary
  of tuple keys. Operates over datasets which are too large to fit in RAM, building the matrix
  in chunks of size chunk_shape.
"""
function overlap_in_chunks( seg1, seg2, chunk_shape, verb=true )

    @assert size(seg1) == size(seg2)
    S, T = eltype(seg1), eltype(seg2)

    cbs = chunk_u.chunk_bounds( size(seg1), chunk_shape )
    counts = Dict{Tuple{S,T},Int}();

    for cb in cbs

      if verb println("Chunk: $cb") end

      seg1ch = chunk_u.fetch_chunk( seg1, cb )
      seg2ch = chunk_u.fetch_chunk( seg2, cb )

      segm_overlap!( seg1ch, seg2ch, counts )

    end

    counts
end

"""

    om_from_impl( impl )

  Creates an official Julia sparse matrix from the "implicit" dict of keys
  version returned by segm_overlap, overlap_in_chunks, etc.
"""
function om_from_impl( impl )

    #find max keys
    max_r, max_c = 0,0
    for (k,v) in impl
      max_r = max(k[1],max_r)
      max_c = max(k[2],max_c)
    end

    om = spzeros(Int, max_r, max_c)

    for (k,v) in impl
      om[k[1],k[2]] = v
    end

    om
end

function create_segID_map(segA, segB)
    @assert size(segA) == size(segB)
    segA_ids = unique(segA)
    return create_segID_map(segA, segB, segA_ids)
end

"""
Create map of overlapping IDs from segA to segB
"""
function create_segID_map(segA, segB, segA_ids)
    @assert size(segA) == size(segB)
    id_map = Dict{Int64, Array{UInt32,1}}()
    for id in segA_ids
        id_mask = segA .== id
        id_map[id] = unique(segB[id_mask])
    end
    return id_map
end

"""
Count number of splits per seg ID
"""
function count_splits(segID_map)
    splits = Dict()
    for (segID, overlapIDs) in segID_map
        splits[segID] = length(overlapIDs)
    end
    return splits
end

"""
Splits should be interpreted like this:
    0: seg ID not in the original segmentation
    1: there was no merger in this segment
   >1: this segment was a merger, and was split x number of times
        (the sum of these splits, minus one per segment, should equal
         the number of seg IDs not in the original segmentation)

Merges should be interpreted like this:
    This segment required x number of merges to reach its final segmentation
"""
function append_splits(seg_list, split_dict)
    splits = []
    for i in 1:size(seg_list,1)
        if haskey(split_dict, seg_list[i])
            push!(splits, split_dict[seg_list[i]])
        else
            push!(splits, 0)
        end
    end
    return hcat(seg_list, splits)
end

"""
Filter list to be only segments with splits
"""
function get_splits(seg_list, label=2)
    s = copy(seg_list)
    s[:,3] = s[:,3]-1
    splits = s[s[:,2] .== label, :]
    return splits[splits[:,3] .!= 0, :]
end


text:
# module for representing agent
# decide control order to robot

mutable struct Agent
    speed
    yaw_rate
    time_interval
    estimator
    prev_spd
    prev_yr

    # init
    function Agent(speed::Float64, yaw_rate::Float64;
                   time_interval::Float64=0.1,
                   estimator=nothing)
        self = new()
        self.speed = speed
        self.yaw_rate = yaw_rate
        self.time_interval = time_interval
        self.estimator = estimator
        self.prev_spd = 0.0
        self.prev_yr = 0.0
        return self
    end
end

function draw_decision!(self::Agent, observation)
    if self.estimator !== nothing
        motion_update(self.estimator, self.prev_spd, self.prev_yr, self.time_interval)
        self.prev_spd, self.prev_yr = self.speed, self.yaw_rate
        observation_update(self.estimator, observation)
        draw!(self.estimator)
    end
    return self.speed, self.yaw_rate
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
const _breakpoints = AbstractBreakpoint[]

"""
    breakpoints()::Vector{AbstractBreakpoint}

Return an array with all breakpoints.
"""
breakpoints() = _breakpoints

const breakpoint_update_hooks = []
"""
    on_breakpoints_updated(f)

Register a one-argument function to be called after any update to the set of all
breakpoints. This includes their creation, deletion, enabling and disabling.

The function `f` should take two inputs:
 - First argument is the function doing to update, this is provided to allow to dispatch
   on its type. It will be one:
   -  `::typeof(breakpoint)` for the creation,
   -  `::typeof(remove)` for the deletion.
   -  `::typeof(update_states)` for disable/enable/toggleing
 - Second argument is the breakpoint object that was changed.

If only desiring to handle some kinds of update, `f` should have fallback methods
to do nothing in the general case.

!!! warning
    This feature is experimental, and may be modified or removed in a minor release.
"""
on_breakpoints_updated(f) = push!(breakpoint_update_hooks, f)

"""
    firehooks(hooked_fun, bp::AbstractBreakpoint)

Trigger all hooks that were registered with [`on_breakpoints_updated`](@ref),
passing them the `hooked_fun` and the `bp`.
This should be called whenever the set of breakpoints is updated.
`hooked_fun` is the function doing the update, and `bp` is the relevent breakpoint being
updated _after_ the update is applied.

!!! warning
    This feature is experimental, and may be modified or removed in a minor release.
"""
function firehooks(hooked_fun, bp::AbstractBreakpoint)
    for hook in breakpoint_update_hooks
        try
            hook(hooked_fun, bp)
        catch err
            @warn "Hook function errored" hook hooked_fun bp exception=err
        end
    end
end

function add_to_existing_framecodes(bp::AbstractBreakpoint)
    for framecode in values(framedict)
        add_breakpoint_if_match!(framecode, bp)
    end
end

function add_breakpoint_if_match!(framecode::FrameCode, bp::AbstractBreakpoint)
    if framecode_matches_breakpoint(framecode, bp)
        stmtidx = bp.line === 0 ? 1 : statementnumber(framecode, bp.line)
        breakpoint!(framecode, stmtidx, bp.condition, bp.enabled[])
        push!(bp.instances, BreakpointRef(framecode, stmtidx))
    end
end

function framecode_matches_breakpoint(framecode::FrameCode, bp::BreakpointSignature)
    function extract_function_from_method(m::Method)
        sig = Base.unwrap_unionall(m.sig)
        ft0 = sig.parameters[1]
        ft = Base.unwrap_unionall(ft0)
        if ft <: Function && isa(ft, DataType) && isdefined(ft, :instance)
            return ft.instance
        elseif isa(ft, DataType) && ft.name === Type.body.name
            return ft.parameters[1]
        else
            return ft
        end
    end

    framecode.scope isa Method || return false
    meth = framecode.scope
    bp.f isa Method && return meth === bp.f
    bp.f === extract_function_from_method(meth) || return false
    bp.sig === nothing && return true
    return bp.sig <: meth.sig
end

"""
    breakpoint(f, [sig], [line], [condition])

Add a breakpoint to `f` with the specified argument types `sig`.¨
If `sig` is not given, the breakpoint will apply to all methods of `f`.
If `f` is a method, the breakpoint will only apply to that method.
Optionally specify an absolute line number `line` in the source file; the default
is to break upon entry at the first line of the body.
Without `condition`, the breakpoint will be triggered every time it is encountered;
the second only if `condition` evaluates to `true`.
`condition` should be written in terms of the arguments and local variables of `f`.

# Example
```julia
function radius2(x, y)
    return x^2 + y^2
end

breakpoint(radius2, Tuple{Int,Int}, :(y > x))
```
"""
function breakpoint(f::Union{Method, Function}, sig=nothing, line::Integer=0, condition::Condition=nothing)
    if sig !== nothing && f isa Function
        sig = Base.to_tuple_type(sig)
        sig = Tuple{typeof(f), sig.parameters...}
    end
    bp = BreakpointSignature(f, sig, line, condition, Ref(true), BreakpointRef[])
    add_to_existing_framecodes(bp)
    idx = findfirst(bp2 -> same_location(bp, bp2), _breakpoints)
    if idx === nothing  # creating new
        push!(_breakpoints, bp)
    else  #Replace existing breakpoint
        old_bp = _breakpoints[idx]
        _breakpoints[idx] = bp
        firehooks(remove, old_bp)
    end
    firehooks(breakpoint, bp)
    return bp
end
breakpoint(f::Union{Method, Function}, sig, condition::Condition) = breakpoint(f, sig, 0, condition)
breakpoint(f::Union{Method, Function}, line::Integer, condition::Condition=nothing) = breakpoint(f, nothing, line, condition)
breakpoint(f::Union{Method, Function}, condition::Condition) = breakpoint(f, nothing, 0, condition)

"""
    breakpoint(file, line, [condition])

Set a breakpoint in `file` at `line`. The argument `file` can be a filename, a partial path or absolute path.
For example, `file = foo.jl` will match against all files with the name `foo.jl`,
`file = src/foo.jl` will match against all paths containing `src/foo.jl`, e.g. both `Foo/src/foo.jl` and `Bar/src/foo.jl`.
Absolute paths only matches against the file with that exact absolute path.
"""
function breakpoint(file::AbstractString, line::Integer, condition::Condition=nothing)
    file = normpath(file)
    apath = CodeTracking.maybe_fix_path(abspath(file))
    ispath(apath) && (apath = realpath(apath))
    bp = BreakpointFileLocation(file, apath, line, condition, Ref(true), BreakpointRef[])
    add_to_existing_framecodes(bp)
    idx = findfirst(bp2 -> same_location(bp, bp2), _breakpoints)
    idx === nothing ? push!(_breakpoints, bp) : (_breakpoints[idx] = bp)
    firehooks(breakpoint, bp)
    return bp
end

function framecode_matches_breakpoint(framecode::FrameCode, bp::BreakpointFileLocation)
    if framecode.scope isa Method
        meth = framecode.scope::Method
        methpath = CodeTracking.maybe_fix_path(String(meth.file))
        ispath(methpath) && (methpath = realpath(methpath))
        if bp.abspath == methpath || endswith(methpath, bp.path)
            return method_contains_line(meth, bp.line)
        else
            return false
        end
    else
        path, _ = whereis(framecode, 1)
        path = CodeTracking.maybe_fix_path(path)
        ispath(path) && (path = realpath(path))

        if bp.abspath == path || endswith(path, bp.path)
            return toplevel_code_contains_line(framecode, bp.line)
        else
            return false
        end
    end
end

function shouldbreak(frame::Frame, pc::Int)
    bps = frame.framecode.breakpoints
    isassigned(bps, pc) || return false
    bp = bps[pc]
    bp.isactive || return false
    return Base.invokelatest(bp.condition, frame)::Bool
end

function prepare_slotfunction(framecode::FrameCode, body::Union{Symbol,Expr})
    ismeth = framecode.scope isa Method
    uslotnames = Set{Symbol}()
    slotnames  = Symbol[]
    for name in framecode.src.slotnames
        if name ∉ uslotnames
            push!(slotnames, name)
            push!(uslotnames, name)
        end
    end
    framename, dataname = gensym("frame"), gensym("data")
    assignments = Expr[:($dataname = $framename.framedata)]
    default = Unassigned()
    for i = 1:length(slotnames)
        slotname = framecode.src.slotnames[i]
        qslotname = QuoteNode(slotname)
        list = framecode.slotnamelists[slotname]
        if length(list) == 1
            maxexpr = :($dataname.last_reference[$(list[1])] > 0 ? $(list[1]) : 0)
        else
            maxcounter, maxidx = gensym("maxcounter"), gensym("maxidx")
            maxexpr = quote
                begin
                    $maxcounter, $maxidx = 0, 0
                    for l in $list
                        counter = $dataname.last_reference[l]
                        if counter > $maxcounter
                            $maxcounter, $maxidx = counter, l
                        end
                    end
                    $maxidx
                end
            end
        end
        maxexsym = gensym("slotid")
        push!(assignments, :($maxexsym = $maxexpr))
        push!(assignments, :($slotname = $maxexsym > 0 ? something($dataname.locals[$maxexsym]) : $default))
    end
    if ismeth
        syms = sparam_syms(framecode.scope)
        for i = 1:length(syms)
            push!(assignments, Expr(:(=), syms[i], :($dataname.sparams[$i])))
        end
    end
    funcname = ismeth ? gensym("slotfunction") : gensym(Symbol(framecode.scope, "_slotfunction"))
    return Expr(:function, Expr(:call, funcname, framename), Expr(:block, assignments..., body))
end

_unpack(condition) = isa(condition, Expr) ? (Main, condition) : condition

## The fundamental implementations of breakpoint-setting
function breakpoint!(framecode::FrameCode, pc, condition::Condition=nothing, enabled=true)
    stmtidx = pc
    if condition === nothing
        framecode.breakpoints[stmtidx] = BreakpointState(enabled)
    else
        mod, cond = _unpack(condition)
        fex = prepare_slotfunction(framecode, cond)
        framecode.breakpoints[stmtidx] = BreakpointState(enabled, Core.eval(mod, fex))
    end
end
breakpoint!(frame::Frame, pc=frame.pc, condition::Condition=nothing) =
    breakpoint!(frame.framecode, pc, condition)

function update_states!(bp::AbstractBreakpoint)
    foreach(bpref -> update_state!(bpref, bp.enabled[]), bp.instances)
    firehooks(update_states!, bp)
end
update_state!(bp::BreakpointRef, v::Bool) = bp[] = v

"""
    enable(bp::AbstractBreakpoint)

Enable breakpoint `bp`.
"""
enable(bp::AbstractBreakpoint) = (bp.enabled[] = true; update_states!(bp))
enable(bp::BreakpointRef) = update_state!(bp, true)

"""
    disable(bp::AbstractBreakpoint)

Disable breakpoint `bp`. Disabled breakpoints can be re-enabled with [`enable`](@ref).
"""
disable(bp::AbstractBreakpoint) = (bp.enabled[] = false; update_states!(bp))
disable(bp::BreakpointRef) = update_state!(bp, false)

"""
    remove(bp::AbstractBreakpoint)

Remove (delete) breakpoint `bp`. Removed breakpoints cannot be re-enabled.
"""
function remove(bp::AbstractBreakpoint)
    idx = findfirst(isequal(bp), _breakpoints)
    if idx !== nothing
        bp = _breakpoints[idx]
        deleteat!(_breakpoints, idx)
        firehooks(remove, bp)
    end
    foreach(remove, bp.instances)
end
function remove(bp::BreakpointRef)
    bp.framecode.breakpoints[bp.stmtidx] = BreakpointState(false, falsecondition)
    return nothing
end

"""
    toggle(bp::AbstractBreakpoint)

Toggle breakpoint `bp`.
"""
toggle(bp::AbstractBreakpoint) = (bp.enabled[] = !bp.enabled[]; update_states!(bp))
toggle(bp::BreakpointRef) = update_state!(bp, !bp[].isactive)

"""
    enable()

Enable all breakpoints.
"""
enable() = foreach(enable, _breakpoints)

"""
    disable()

Disable all breakpoints.
"""
disable() = foreach(disable, _breakpoints)

"""
    remove()

Remove all breakpoints.
"""
function remove()
    for bp in _breakpoints
        foreach(remove, bp.instances)
    end
    empty!(_breakpoints)
end

"""
    break_on(states...)

Turn on automatic breakpoints when any of the conditions described in `states` occurs.
The supported states are:

- `:error`: trigger a breakpoint any time an uncaught exception is thrown
- `:throw` : trigger a breakpoint any time a throw is executed (even if it will eventually be caught)
"""
function break_on(states::Vararg{Symbol})
    for state in states
        if state === :error
            break_on_error[] = true
        elseif state === :throw
            break_on_throw[] = true
        else
            throw(ArgumentError(string("unsupported state :", state)))
        end
    end
end

"""
    break_off(states...)

Turn off automatic breakpoints when any of the conditions described in `states` occurs.
See [`break_on`](@ref) for a description of valid states.
"""
function break_off(states::Vararg{Symbol})
    for state in states
        if state === :error
            break_on_error[] = false
        elseif state === :throw
            break_on_throw[] = false
        else
            throw(ArgumentError(string("unsupported state :", state)))
        end
    end
end

"""
    @breakpoint f(args...) condition=nothing
    @breakpoint f(args...) line condition=nothing

Break upon entry, or at the specified line number, in the method called by `f(args...)`.
Optionally supply a condition expressed in terms of the arguments and internal variables
of the method.
If `line` is supplied, it must be a literal integer.

# Example

Suppose a method `mysum` is defined as follows, where the numbers to the left are the line
number in the file:

```
12 function mysum(A)
13     s = zero(eltype(A))
14     for a in A
15         s += a
16     end
17     return s
18 end
```

Then

```
@breakpoint mysum(A) 15 s>10
```

would cause execution of the loop to break whenever `s>10`.
"""
macro breakpoint(call_expr, args...)
    whichexpr = InteractiveUtils.gen_call_with_extracted_types(__module__, :which, call_expr)
    haveline, line, condition = false, 0, nothing
    while !isempty(args)
        arg = first(args)
        if isa(arg, Integer)
            haveline, line = true, arg
        else
            condition = arg
        end
        args = Base.tail(args)
    end
    condexpr = condition === nothing ? nothing : esc(Expr(:quote, condition))
    if haveline
        return quote
            local method = $whichexpr
            $breakpoint(method, $line, $condexpr)
        end
    else
        return quote
            local method = $whichexpr
            $breakpoint(method, $condexpr)
        end
    end
end

const __BREAKPOINT_MARKER__ = nothing

"""
    @bp

Insert a breakpoint at a location in the source code.
"""
macro bp()
    return esc(:($(JuliaInterpreter).__BREAKPOINT_MARKER__))
end


text:
using Pkg;
Pkg.activate("knitro");
using CUTEst, Dates, NLPModels, SolverBenchmark
using NLPModelsKnitro

nmax = 10
problems = readlines("list_problems_$nmax.dat")
cutest_problems = (CUTEstModel(p) for p in problems)

max_time = 3600.0
max_ev = typemax(Int)
max_iter = typemax(Int)
atol = 1e-5
rtol = 1e-6

solvers = Dict(
    :knitro =>
        nlp -> knitro(
            nlp,
            opttol = rtol,
            opttol_abs = atol,
            maxfevals = max_ev,
            maxit = max_iter,
            maxtime_real = max_time,
            out_hints = 0,
            outlev = 0,
        ),
)
stats = bmark_solvers(solvers, cutest_problems)

using JLD2
@save "$(today())_$(prod(String.(keys(solvers)) .* "_"))cutest_$(string(length(problems)))_$(nmax).jld2" stats


text:
using KernelAbstractions, Test
using CUDAapi
if CUDAapi.has_cuda_gpu()
    using CuArrays
    CuArrays.allowscalar(false)
end

@kernel function kernel_print()
  I = @index(Global)
  @print("Hello from thread ", I, "!\n")
end

function test_print(backend)
  kernel = kernel_print(backend, 4)
  kernel(ndrange=(4,))
end

@testset "print test" begin
    if CUDAapi.has_cuda_gpu()
        wait(test_print(CUDA()))
        @test true
    end

    wait(test_print(CPU()))
    @test true

    @print("Why this should work")
    @test true
end


text:
struct VGrid{T <: Real, N <: Int}
                   nxv :: N                                     # Number of frid cells in x-direction in UGrid
                   nyv :: N                                     # Number of grid cells in y-direction in UGrid
                 mask :: Array{Bool,2}                         # Mask specifying model domain wrt U grid
                    n :: N                                     # Total number of cells in model domain
                 crop :: Diagonal{T,Array{T,1}}                # Crop matrix: diagonal matrix with mask entries on diag
                 samp :: SparseMatrixCSC{T,N}                  # Sampling matrix: take full domain to model domain
               spread :: SparseMatrixCSC{T,N}                  # Spread matrix: take model domain to full domain
                 cent :: KronType{T,N}                         # Map from U grid to H grid
                   ∂x :: KronType{T,N}                         # Matrix representation of differentiation wrt x
                   ∂y :: KronType{T,N}                         # Matrix representation of differentiation wrt y
               levels :: N                                     # Number of levels in the preconditioner
                  dwt :: KronType{T,N}                         # Wavelet matrix product on u grid
                    s :: Array{T,2}                            # Ice surface elevation
                    h :: Array{T,2}                            # Ice thickness
    grounded_fraction :: Array{T,2}                            # Grid cell grounded fraction
                 βeff :: Array{T,2}                            # Effective β value on u grid(eqn 12 in Arthern 2015 JGeophysRes)
             dnegβeff :: Base.RefValue{Diagonal{T,Array{T,1}}} # Rheological operator (diagonal of βeff)
                    v :: Array{T,2}                            # Ice velocity in the y-direction
end

"""
    VGrid(;
            nxv,
            nyv,
            mask = trues(nxv,nyv),
            levels,
            dx,
            dy)

Construct a WAVI.jl VGrid with size (nxv,nyv)
VGrid stores fields that are defined on the problem's V grid.
(Co-ordinates of HGrid stored in a Grid under xxv and yyv fields)

Keyword arguments
=================
    - 'nxv': (required) Number of grid cells in x-direction in VGrid (should be same as grid.nx)
            Note that we store the grid size here, even though it can be easily inferred from grid, to increase transparency in velocity solve.
    - 'nyv': (required) Number of grid cells in y-direction in VGrid (should be same as grid.ny + 1)
    - 'mask': Mask specifying the model domain with respect to V grid
    - levels: (required) Number of levels in the preconditioner
    - dx: (required) Grid spacing in the x direction
    - dy: (required) Grid spacing in the y direction
"""
function VGrid(;
        nxv,
        nyv,
        mask = trues(nxv,nyv),
        levels,
        dx,
        dy)

    #check the sizes of inputs
    (size(mask) == (nxv,nyv)) || throw(DimensionMismatch("Sizes of inputs to UGrid must all be equal to nxv x nyv (i.e. $nxv x $nyv)"))

    #construct operators
    n = count(mask)
    crop = Diagonal(float(mask[:]))
    samp = crop[mask[:],:]
    spread = sparse(samp')
    cent = c(nyv-1) ⊗ spI(nxv)
    ∂x = χ(nyv-2) ⊗ ∂1d(nxv-1,dx)
    ∂y = ∂1d(nyv-1,dy) ⊗ spI(nxv)
    dωt = wavelet_matrix(nyv,levels,"forward" ) ⊗ wavelet_matrix(nxv,levels,"forward")

    #fields stored on UGrid
    s = zeros(nxv,nyv)
    h = zeros(nxv,nyv)
    grounded_fraction = ones(nxv,nyv)
    βeff = zeros(nxv,nyv)
    dnegβeff = Ref(crop*Diagonal(-βeff[:])*crop)
    v = zeros(nxv,nyv)

    #size assertions
    @assert size(mask)==(nxv,nyv)
    n == count(mask)
    @assert crop == Diagonal(float(mask[:]))
    @assert samp == crop[mask[:],:]
    @assert spread == sparse(samp')
    @assert size(s)==(nxv,nyv)
    @assert size(h)==(nxv,nyv)
    @assert size(grounded_fraction)==(nxv,nyv)
    @assert size(βeff)==(nxv,nyv)
    @assert size(v)==(nxv,nyv)

    #make sure boolean type rather than bitarray
    mask = convert(Array{Bool,2}, mask)

    return VGrid(
                nxv,
                nyv,
                mask,
                n,
                crop,
                samp,
                spread,
                cent,
                ∂x,
                ∂y,
                levels,
                dωt,
                s,
                h,
                grounded_fraction,
                βeff,
                dnegβeff,
                v)
end


text:
module SimpleTropical

import Base: isinf, show, (+), (*), inv, (==), isequal, (^)
import Base: convert, zero, one, string

export Tropical, TropicalInf, ⊕, ⊗, long_tropical_show

_long_show = true

"""
`Tropical(x::T) where T<:Real` defines a new `Tropical`
number.
"""
struct Tropical{T<:Real} <: Number
    val::T
    inf_flag::Bool

    function Tropical{T}(xx::Real, ii::Bool = false) where {T}
        TT = typeof(xx)
        if isinf(xx) || ii
            return new(zero(TT), true)
        end
        return new(xx, false)
    end
end

Tropical(x::T) where {T<:Real} = Tropical{T}(x)
function Tropical(x::T, i::Bool) where {T<:Real}
    if i
        return Tropical{T}(zero(T), true)
    end
    return Tropical(x)
end

"""
`TropicalInf` is a constant that represents infinity in the tropical
semiring.
"""
const TropicalInf = Tropical{Bool}(0, true)

isinf(X::Tropical) = X.inf_flag

Base.promote_rule(::Type{Tropical{T}}, ::Type{S}) where {T<:Real,S<:Real} =
    Tropical{promote_type(T, S)}
Base.promote_rule(::Type{Tropical{T}}, ::Type{Tropical{S}}) where {T<:Real,S<:Real} =
    Tropical{promote_type(T, S)}

convert(::Type{Tropical}, x::T) where {T<:Real} = Tropical{T}(x)
convert(::Type{Tropical{T}}, x::S) where {T<:Real,S<:Tropical} =
    Tropical(convert(T, x.val), x.inf_flag)

function long_tropical_show(t::Bool)::Bool
    global _long_show = t
end

"""
`long_tropical_show(t::Bool)` determines the display style for
`Tropical` numbers.
+ When `t` is `true`: display numbers like this: `Tropical(5)` or `Tropical(∞)`.
+ When `t` is `false`: display numbers like this: `5` or `∞`.

`long_tropical_show()` returns the current state (`true` or `false`).
"""
long_tropical_show()::Bool = _long_show

function string(x::Tropical{T})::String where T
    val = string(x.val)
    if isinf(x)
        val = "∞"
    end
    if _long_show
        return "Tropical($val)"
    else
        return val
    end
end

function show(io::IO, t::Tropical)
    print(io, string(t))
end

# Calling zero(Tropical) returns Tropical(∞) because that's the identity
# element of addition. Likewise, one(Tropical) returns Tropical(0) because
# that's the identity element of multiplication.

zero(::Tropical{T}) where {T} = TropicalInf
zero(::Type{Tropical}) = TropicalInf
zero(::Type{Tropical{T}}) where {T} = Tropical(zero(T), true)

one(::Tropical{T}) where {T} = Tropical{T}(0)
one(::Type{Tropical}) = Tropical(0)
one(::Type{Tropical{T}}) where {T} = Tropical{T}(0)

function (+)(x::Tropical{T}, y::Tropical{T}) where {T}
    if isinf(x)
        if isinf(y)   # when X,Y both are infinite
            return Tropical(zero(T), true)  # create common infinite
        else
            return Tropical(y)
        end
    end

    if isinf(y)
        return Tropical(x)
    end

    return Tropical(min(x.val, y.val))
end
(+)(x::Tropical{T}, y::Tropical{S}) where {T,S} = +(promote(x, y)...)
(+)(x::Tropical{T}, y::Real) where {T} = +(promote(x, y)...)
(+)(x::Real, y::Tropical{T}) where {T} = +(promote(x, y)...)

function (*)(x::Tropical{T}, y::Tropical{T}) where {T}
    if isinf(x) || isinf(y)
        return Tropical(zero(T), true)
    end

    return Tropical(x.val + y.val)
end
(*)(x::Tropical{T}, y::Tropical{S}) where {T,S} = *(promote(x, y)...)
(*)(x::Tropical{T}, y::Real) where {T} = *(promote(x, y)...)
(*)(x::Real, y::Tropical{T}) where {T} = *(promote(x, y)...)

(⊕)(x::S, y::T) where {S<:Union{Tropical,Real},T<:Union{Tropical,Real}} =
    Tropical(x) + Tropical(y)
(⊗)(x::S, y::T) where {S<:Union{Tropical,Real},T<:Union{Tropical,Real}} =
    Tropical(x) * Tropical(y)

function inv(X::Tropical)
    @assert !isinf(X) "TropicalInf is not invertible"
    return Tropical(-X.val)
end

function (^)(X::Tropical, p::Integer)
    if isinf(X)
        @assert p > 0 "Cannot raise tropical infinity to a nonpositive power."
        return X
    end

    return Tropical(X.val * p)
end

function isequal(X::Tropical, Y::Tropical)
    if !isinf(X) && !isinf(Y)
        return isequal(X.val, Y.val)
    else
        return isinf(X) && isinf(Y)
    end
end

function ==(X::Tropical, Y::Tropical)
    if !isinf(X) && !isinf(Y)
        return X.val == Y.val
    else
        return isinf(X) && isinf(Y)
    end
end

end # end of module


text:
export iauFane03
"""
Fundamental argument, IERS Conventions (2003):
mean longitude of Neptune.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   t     double    TDB, Julian centuries since J2000.0 (Note 1)

Returned (function value):
         double    mean longitude of Neptune, radians (Note 2)

Notes:

   1. Though t is strictly TDB, it is usually more convenient to use
   TT, which makes no significant difference.

   2. The expression used is as adopted in IERS Conventions (2003) and
   is adapted from Simon et al. (1994).

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
   Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauFane03(double t)
function iauFane03(t::Real)
   return ccall((:iauFane03, libsofa_c), Cdouble, (Cdouble,), convert(Float64, t))
end


text:
#!/usr/bin/julia

__precompile__(true)

module FunctionWrappers

# Used to bypass NULL check
@inline function assume(v::Bool)
    Base.llvmcall(("declare void @llvm.assume(i1)",
                   """
                   %v = trunc i8 %0 to i1
                   call void @llvm.assume(i1 %v)
                   ret void
                   """), Void, Tuple{Bool}, v)
end

# Convert return type and generates cfunction signatures
Base.@pure map_argtype(T) = (isbits(T) || T === Any) ? T : Ref{T}
Base.@pure get_cfunc_argtype(Obj, Args) =
    Tuple{Ref{Obj}, (map_argtype(Arg) for Arg in Args.parameters)...}

# Call wrapper since `cfunction` does not support non-function
# or closures
immutable CallWrapper{Ret} <: Function
end
(::CallWrapper{Ret}){Ret}(f, args...)::Ret = f(args...)

# Specialized wrapper for
for nargs in 0:128
    @eval (::CallWrapper{Ret}){Ret}(f, $((Symbol("arg", i) for i in 1:nargs)...))::Ret =
        f($((Symbol("arg", i) for i in 1:nargs)...))
end

type FunctionWrapper{Ret,Args<:Tuple}
    ptr::Ptr{Void}
    objptr::Ptr{Void}
    obj
    objT
    function FunctionWrapper{objT}(obj::objT)
        objref = Base.cconvert(Ref{objT}, obj)
        new(cfunction(CallWrapper{Ret}(), map_argtype(Ret),
                      get_cfunc_argtype(objT, Args)),
            Base.unsafe_convert(Ref{objT}, objref), objref, objT)
    end
    FunctionWrapper(obj::FunctionWrapper{Ret,Args}) = obj
end

Base.convert{T<:FunctionWrapper}(::Type{T}, obj) = T(obj)

@noinline function reinit_wrapper{Ret,Args}(f::FunctionWrapper{Ret,Args})
    objref = f.obj
    objT = f.objT
    ptr = cfunction(CallWrapper{Ret}(), map_argtype(Ret),
                    get_cfunc_argtype(objT, Args))
    f.ptr = ptr
    f.objptr = Base.unsafe_convert(Ref{objT}, objref)
    return ptr
end

@generated function do_ccall{Ret,Args}(f::FunctionWrapper{Ret,Args}, args::Args)
    # Has to be generated since the arguments type of `ccall` does not allow
    # anything other than tuple (i.e. `@pure` function doesn't work).
    quote
        $(Expr(:meta, :inline))
        ptr = f.ptr
        if ptr == C_NULL
            # For precompile support
            ptr = reinit_wrapper(f)
        end
        assume(ptr != C_NULL)
        objptr = f.objptr
        ccall(ptr, $(map_argtype(Ret)),
              (Ptr{Void}, $((map_argtype(Arg) for Arg in Args.parameters)...)),
              objptr, $((:(args[$i]) for i in 1:length(Args.parameters))...))
    end
end

@inline (f::FunctionWrapper)(args...) = do_ccall(f, args)

# Testing only
const identityAnyAny = FunctionWrapper{Any,Tuple{Any}}(identity)

end


text:
@testset "PermutationTracker" begin
    @testset for n = 0
        v = []
        k2i = Int[]
        i2k = Int[]
        heapify!(
            PermutationTracker(v,k2i,i2k),
            By(((k,vk),) -> vk, Forward)
        )
        @test v == []
        @test k2i == []
        @test i2k == []
    end
    @testset for n = 1
        v = [0]
        k2i = Int[1]
        i2k = Int[1]
        heapify!(
            PermutationTracker(v,k2i,i2k),
            By(((k,vk),) -> vk, Forward)
        )
        @test v == [0]
        @test k2i == [1]
        @test i2k == [1]
    end
    @testset for n = 2:10
        Random.seed!(42)
        for _ = 1:10
            v̂ = rand(1:100,n)
            v = copy(v̂)
            k2i = collect(1:n)
            i2k = collect(1:n)
            heapify!(
                PermutationTracker(v,k2i,i2k),
                By(((k,vk),) -> vk, Forward)
            )
            @test isheap(v,Forward)
            @test v[k2i] == v̂
            @test k2i[i2k] == 1:n
        end
    end
end


text:
module Land

using DiffEqOperators
using DocStringExtensions
using OrdinaryDiffEq
using Parameters
using PkgUtility: T_25
using WaterPhysics: TraceGasAir, TraceGasCO₂, diffusive_coefficient

# export public types
export VerticalLayers

# export public functions
export vertical_layers!

# include sub modules
include("CanopyLayers/CanopyLayers.jl"                  )
include("Photosynthesis/Photosynthesis.jl"              )
include("PlantHydraulics/PlantHydraulics.jl"            )
include("StomataModels/StomataModels.jl"                )
include("SoilPlantAirContinuum/SoilPlantAirContinuum.jl")

# include types
include("Land/Types.jl"    )
include("Land/AirLayers.jl")

end # module

max_stars_repo_name: climate-machine/Land

max_issues_repo_name: climate-machine/Land

max_forks_repo_name: climate-machine/Land


text:
import RData
using RNGPool
using CoupledConditionalSMC

# data generated with A = 0.95 (not 0.9)
dataJLS = RData.load("demo/ar1data.RData")
ys = vec(dataJLS["observations"])

include("../demo/lgModel.jl")

# they start at time 0 with observations from time 1,2,...
# so we have an initial distribution of N(0,1+0.9^2)
theta = LGTheta(0.9, 1.0, 1.0, 1.0, 0.0, 1.81)
model = makeLGModel(theta, ys)
lM = LinearGaussian.makelM(theta)
ko = kalman(theta, ys)

function checkBoundary(ref1::Vector{Float64Particle}, ref2::Vector{Float64Particle}) :: Int64
    b = 0
    for i in 1:length(ref1)
        (ref1[i].x != ref2[i].x) && break
        b = i
    end
    b
end

function couplingBoundary(model, lM, algorithm=:AS, N=2000, n=500, maxit=500,
    rngCouple = true, independentInitialisation = true)
    io = CCSMCIO{model.particle, model.pScratch}(N, n)
    initializeCCSMC(model, lM, io, algorithm, independentInitialisation)
    bd = fill(n, maxit)
    for i in 1:maxit
        ccXpf!(model, lM, io, algorithm, rngCouple)
        bd_ = checkBoundary(io.ref1, io.ref2)
        bd[i] = bd_
        (bd_ == n) && break
    end
    bd
end

setRNGs(12345)
N = 64    # Number of particles
n = 800     # Length of time series
maxit = 500 # Maximum number of iterations
bdAT = couplingBoundary(model, lM, :AT, N, n, maxit)
bdAS = couplingBoundary(model, lM, :AS, N, n, maxit)
bdBS = couplingBoundary(model, lM, :BS, N, n, maxit)

using GR
figure(figsize=(15,5))
beginprint("paper/output/lgDemoPaths.pdf")
GR.legend("BS", "AS", "ATT"; loc=2)
GR.xlabel("Iteration")
GR.ylabel("Coupling boundary")
GR.ylim([-30, n+30])
GR.xlim([0, maxit])
GR.plot(1:maxit, bdBS, "k-.", 1:maxit, bdAS, "k--", 1:maxit, bdAT, "k-")
endprint()


text:
using eFEMpart

function brinkmanMP_example()
  x0 = 2.0
  mesh = squareMeshFluid([-x0,x0,-1,1],32)

  μ0 = 1.0
  a10 = 0.0
  a20 = 1.0

  μ  = μ0*ones(length(mesh.xy))
  a1  = a10*ones(length(mesh.xy))
  a2  = a20*ones(length(mesh.xy))
  G  = 1e0
  param = BrinkmanMPParam(μ,a1,a2)

  OperatorType = :BrinkmanMP2D

  dNodes = Dirichlet(:left,:top,:bottom)
  nNodes = Neumann(:right)
  Nodes = [dNodes,nNodes]

  dUBCf = Dirichlet((x,y) -> (x==-x0 ? 1.0-y^2 : 0.0))
  dVBCf = Dirichlet((x,y) -> 0.0)
  bcfun = [dUBCf,dVBCf]

  prob = Problem(mesh,Nodes,bcfun,OperatorType)
  sol = solve(prob,mesh,param)

  vtkname = Path("solutions","brinkmanMP_example")
  sd = ScalarData(sol.p)
  sn = ScalarNames("pressure")
  vd = VectorData([sol.u,sol.v])
  vn = VectorNames("velocity")
  vtksave(mesh,sd,sn,vd,vn,vtkname)

  return nothing
end

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((1+x^4+x^6+-1*x^3+-1*x^5)*exp(x), x) == :(871*exp(x)+x^6*exp(x)+-870*x*exp(x)+-145*x^3*exp(x)+-7*x^5*exp(x)+36*x^4*exp(x)+435*x^2*exp(x))
@test integrate((x^3+2x)^-1*(2+-1*x^2)*exp(x*(2+x^2)^-1), x) == :(expintegralei(x*(2+x^2)^-1))
@test integrate((x^3+2x)^-1*(2+-1*x^3+2x+2*x^4+3*x^2)*exp(x*(2+x^2)^-1), x) == :((2+x^2)*exp(x*(2+x^2)^-1)+expintegralei(x*(2+x^2)^-1))
@test integrate((x+exp(x))^-1*(1+exp(x))*exp(x+exp(x)), x) == :(expintegralei(x+exp(x)))
@test integrate((1+x^3+-1x+-1*x^2)^-1*(1+x^3+-1*x^2+-3x)*exp((-1+x^2)^-1), x) == :((1+x)*exp((-1+x^2)^-1))
@test integrate(log(x)^-2*(-1+log(x)^2)*exp(1+log(x)^-1), x) == :(x*exp(1+log(x)^-1))
@test integrate(log(x)^-2*(-1+log(x)^2*(1+x))*exp(x+log(x)^-1), x) == :(x*exp(x+log(x)^-1))

max_stars_repo_path: test/0 Independent test suites/Hebisch Problems.jl


text:
"""
    NonadiabaticDistributions

Module containing types for representing distributions of combined nuclear and electronic configurations.
"""
module NonadiabaticDistributions

"""
    NonadiabaticDistribution

Top level type for all nonadiabatic distributions.

These distributions contain a both electronic and nuclear degrees of freedom.
"""
abstract type NonadiabaticDistribution end

include("nuclear_distributions.jl")

include("electronic_distributions.jl")
export Diabatic
export Adiabatic
export SingleState
export ElectronicPopulation

struct CombinedDistribution{N<:NuclearDistribution, E<:ElectronicDistribution} <: NonadiabaticDistribution
    nuclear::N
    electronic::E
end

Base.:*(N::NuclearDistribution, E::ElectronicDistribution) = CombinedDistribution(N, E)
Base.:*(E::ElectronicDistribution, N::NuclearDistribution) = CombinedDistribution(N, E)

Base.lastindex(dist::CombinedDistribution) = lastindex(dist.nuclear)

end # module


text:
#!/usr/bin/env julia

using ImageMagick
using Images
using ArgParse
using FileIO
using Statistics
import APTDecoder
using PyPlot

s = ArgParseSettings()
@add_arg_table s begin
    "--satellite", "-s"
        help = "the satellite name ( \"NOAA 15\",  \"NOAA 18\",  \"NOAA 19\")"
        required = true
    "--starttime", "-t"
        help = "start time of the capture in UTC. It can be deduced from the file name, if it has the structure like gqrx_20190811_075102_137620000.png."
    "--output-prefix"
        help = "an option without argument, i.e. a flag"
    "wav"
        help = "a positional argument"
        required = true
end

parsed_args = parse_args(ARGS, s)

@show parsed_args
wavname = parsed_args["wav"]
satellite_name = parsed_args["satellite"]

# missing paramter have a value nothing
kwargs = filter(kv -> (kv[2] != nothing) && !(kv[1] ∈ ["wav","satellite"]),parsed_args)

APTDecoder.makeplots(wavname,satellite_name,kwargs...)

max_stars_repo_name: Alexander-Barth/APTDecoder

max_issues_repo_name: Alexander-Barth/APTDecoder

max_forks_repo_name: Alexander-Barth/APTDecoder


text:
module Types

export VersionInterval, VersionSet, Requires, Available, Fixed, merge_requires!, satisfies
import Base: show, isempty, in, intersect, hash, deepcopy_internal

immutable VersionInterval
    lower::VersionNumber
    upper::VersionNumber
end
VersionInterval(lower::VersionNumber) = VersionInterval(lower,typemax(VersionNumber))
VersionInterval() = VersionInterval(typemin(VersionNumber))

show(io::IO, i::VersionInterval) = print(io, "[$(i.lower),$(i.upper))")
isempty(i::VersionInterval) = i.upper <= i.lower
in(v::VersionNumber, i::VersionInterval) = i.lower <= v < i.upper
intersect(a::VersionInterval, b::VersionInterval) = VersionInterval(max(a.lower,b.lower), min(a.upper,b.upper))
==(a::VersionInterval, b::VersionInterval) = a.lower == b.lower && a.upper == b.upper

immutable VersionSet
    intervals::Vector{VersionInterval}
end
function VersionSet(versions::Vector{VersionNumber})
    intervals = VersionInterval[]
    if isempty(versions)
        push!(intervals, VersionInterval())
    else
        isodd(length(versions)) && push!(versions, typemax(VersionNumber))
        while !isempty(versions)
            push!(intervals, VersionInterval(shift!(versions), shift!(versions)))
        end
    end
    VersionSet(intervals)
end
VersionSet(versions::VersionNumber...) = VersionSet(VersionNumber[versions...])

show(io::IO, s::VersionSet) = print_joined(io, s.intervals, " ∪ ")
isempty(s::VersionSet) = all(i->isempty(i), s.intervals)
in(v::VersionNumber, s::VersionSet) = any(i->in(v,i), s.intervals)
function intersect(A::VersionSet, B::VersionSet)
    ivals = vec([ intersect(a,b) for a in A.intervals, b in B.intervals ])
    ivals = copy(ivals) # temporary bandaid for issue #4592
    filter!(i->!isempty(i), ivals)
    sort!(ivals, by=i->i.lower)
    VersionSet(ivals)
end
==(A::VersionSet, B::VersionSet) = A.intervals == B.intervals
hash(s::VersionSet, h::Uint) = hash(s.intervals, h + uint(0x2fd2ca6efa023f44))
deepcopy_internal(vs::VersionSet, ::ObjectIdDict) = VersionSet(copy(vs.intervals))

typealias Requires Dict{ByteString,VersionSet}

function merge_requires!(A::Requires, B::Requires)
    for (pkg,vers) in B
        A[pkg] = haskey(A,pkg) ? intersect(A[pkg],vers) : vers
    end
    return A
end

satisfies(pkg::String, ver::VersionNumber, reqs::Requires) =
    !haskey(reqs, pkg) || in(ver, reqs[pkg])

immutable Available
    sha1::ASCIIString
    requires::Requires
end

==(a::Available, b::Available) = a.sha1 == b.sha1 && a.requires == b.requires

show(io::IO, a::Available) = isempty(a.requires) ?
    print(io, "Available(", repr(a.sha1), ")") :
    print(io, "Available(", repr(a.sha1), ",", a.requires, ")")

immutable Fixed
    version::VersionNumber
    requires::Requires
end
Fixed(v::VersionNumber) = Fixed(v,Requires())

==(a::Fixed, b::Fixed) = a.version == b.version && a.requires == b.requires

show(io::IO, f::Fixed) = isempty(f.requires) ?
    print(io, "Fixed(", repr(f.version), ")") :
    print(io, "Fixed(", repr(f.version), ",", f.requires, ")")

# TODO: Available & Fixed are almost the same – merge them?
# Free could include the same information too, it just isn't
# required by anything that processes these things.

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# productdomain.jl
# Routines for a cartesian product of domains

###################
# Helper functions
###################

"""
Create an eltype that is suitable for a product domain. The result is a tuple
type, where each of the elements is the eltype of the corresponding element
of the product domain.
"""
_product_eltype(t) = t
_product_eltype(::Type{Tuple{T,V}}) where {T<:Real,V<:Real} =
	NTuple{2,promote_type(T,V)}
_product_eltype(::Type{Tuple{T,V}}) where {T<:Complex,V<:Complex} =
	NTuple{2,promote_type(T,V)}
product_eltype(domains...) = _product_eltype(Tuple{map(eltype, domains)...})

"""
Try to simplify the type of a product domain to a type to which it is isomorphic.
The goal is to automatically embed the product domain in ℝ^N if possible.

Examples of simplifications:
`Tuple{Float64,Float64} -> SVector{2,Float64}`
`Tuple{Tuple{Float64,Float64},Float64} -> SVector{3,Float64}`
"""
simplify_product_eltype(::Type{T}) where {T} = T
simplify_product_eltype(::Type{NTuple{N,T}}) where {N,T} = SVector{N,T}

@generated function simplify_product_eltype(::Type{Tuple{SVector{N,T},T}}) where {N,T}
    M = N+1
    quote SVector{$M,T} end
end
@generated function simplify_product_eltype(::Type{Tuple{T,SVector{N,T}}}) where {N,T}
    M = N+1
    quote SVector{$M,T} end
end
@generated function simplify_product_eltype(::Type{Tuple{SVector{N,T},SVector{K,T}}}) where {N,K,T}
    M = N+K
    quote SVector{$M,T} end
end

@generated function simplify_product_eltype(::Type{Tuple{NTuple{N,T},T}}) where {N,T}
    M = N+1
    quote SVector{$M,T} end
end
@generated function simplify_product_eltype(::Type{Tuple{T,NTuple{N,T}}}) where {N,T}
    M = N+1
    quote SVector{$M,T} end
end
@generated function simplify_product_eltype(::Type{Tuple{NTuple{N,T},NTuple{K,T}}}) where {N,K,T}
    M = N+K
    quote SVector{$M,T} end
end

#######################
# Main type definition
#######################

"""
```using DomainSets```

A `ProductDomain` represents the cartesian product of other domains.

A product domain has two eltypes, an internal type `S` and an external type `T`.
The internal type `S` is a tuple containing the eltypes of the elements of the
product domain. The external eltype `T` is a type whose associated space is
isomorphic to that of `S`, but which has been simplified. (See also
`simplify_product_eltype`).

For example, if `S` is `Tuple{Float64,Float64}`, then `T` is `SVector{2,Float64}`.
"""
struct ProductDomain{DD,S,T} <: Domain{T}
	# D is the type of an indexable list of domains, such as a tuple
	domains	::	DD

	# Inner constructor to verify that S and T are correct
	function ProductDomain{DD,S,T}(domains) where {DD,S,T}
		@assert S == product_eltype(domains...)
		@assert isomorphic(spacetype(S),spacetype(T))
		new{DD,S,T}(domains)
	end
end

function ProductDomain(domains...)
    DD = typeof(domains)
    S = product_eltype(domains...)
	T = simplify_product_eltype(S)
    ProductDomain{DD,S,T}(domains)
end

convert(::Type{Domain{SVector{N,T}}}, d::ProductDomain{<:Tuple{Vararg{<:Any,N}}}) where {N,T} =
	ProductDomain(convert.(Domain{T}, d.domains)...)

elements(d::ProductDomain) = d.domains

internal_eltype(::Type{ProductDomain{DD,S,T}}) where {DD,S,T} = S
internal_eltype(::Type{P}) where {P <: ProductDomain} = S
internal_eltype(d::ProductDomain) = internal_eltype(typeof(d))

cross(x::Domain...) = cartesianproduct(x...)

# One can use the cartesianproduct routine to create product domains
cartesianproduct(domains::Domain...) = ProductDomain(domains...)

# We try to avoid creating nested cartesian products
cartesianproduct(d1::Domain, d2::ProductDomain) = cartesianproduct(d1, elements(d2)...)
cartesianproduct(d1::ProductDomain, d2::Domain) = cartesianproduct(elements(d1)..., d2)
cartesianproduct(d1::ProductDomain, d2::ProductDomain) = cartesianproduct(elements(d1)..., elements(d2)...)

^(d::Domain, n::Int) = cartesianproduct(d, n)

^(d::Domain{T}, ::Type{Val{N}}) where {N,T} = ProductDomain{NTuple{N,typeof(d)},NTuple{N,T},SVector{N,T}}(ntuple(i->d, Val{N}))

indomain(x, d::ProductDomain) = _indomain(convert_space(spacetype(internal_eltype(d)), x), d, elements(d))

# The line below allocates a little bit of memory...
_indomain(x, d::ProductDomain, el) = reduce(&, map(in, x, el))

# ...hence these special cases:
_indomain(x::Tuple{A,B}, d::ProductDomain, el) where {A,B} = in(x[1], el[1]) &&
    in(x[2], el[2])
_indomain(x::Tuple{A,B,C}, d::ProductDomain, el) where {A,B,C} = in(x[1], el[1]) &&
    in(x[2], el[2]) && in(x[3], el[3])
_indomain(x::Tuple{A,B,C,D}, d::ProductDomain, el) where {A,B,C,D} = in(x[1], el[1]) &&
    in(x[2], el[2]) && in(x[3], el[3]) && in(x[4], el[4])

approx_indomain(x, d::ProductDomain, tolerance) =
    _approx_indomain(convert_space(spacetype(internal_eltype(d)), x), d, elements(d), tolerance)

_approx_indomain(x, d::ProductDomain, el, tolerance) = reduce(&, map((u,v) -> approx_indomain(u,v,tolerance), x, el))

point_in_domain(d::ProductDomain) = convert_space(spaceof(d), map(point_in_domain, elements(d)))

infimum(d::ProductDomain) = convert_space(spaceof(d), map(infimum, elements(d)))
supremum(d::ProductDomain) = convert_space(spaceof(d), map(supremum, elements(d)))

isempty(d::ProductDomain) = any(isempty, d.domains)

function show(io::IO, t::ProductDomain)
    L = numelements(t)
    for i in 1:L-1
        show(io, element(t, i))
        print(io, " x ")
    end
    show(io, element(t, L))
end


text:
struct Gisette <: Tabular end
function url(::Gisette)
    [
        "https://archive.ics.uci.edu/ml/machine-learning-databases/gisette/gisette_valid.labels",
        "https://archive.ics.uci.edu/ml/machine-learning-databases/gisette/GISETTE/gisette_valid.data",
        "https://archive.ics.uci.edu/ml/machine-learning-databases/gisette/GISETTE/gisette_train.labels",
        "https://archive.ics.uci.edu/ml/machine-learning-databases/gisette/GISETTE/gisette_train.data",
    ]
end

function checksum(::Gisette)
    [
        "a6b857a0448023f033c4dda2ef848714b4be2ae45ce598d088fb3efb406e08c5",
        "5cea897956dd172a006132738254a27a8f61ecc1ceb6f5b20639c281d2942254",
        "42bd681fe51b161f033df773df14a0116e492676555ab14616c1b72edc054075",
        "6d4c5e998afe67937b9e77a3334e03c85e545ebc65a6eb1333ffc14125cfc389"
    ]
end
function prep(::Gisette)
    [
        path -> preprocess(path, Labels),
        path -> preprocess(path, Gisette()),
        path -> preprocess(path, Labels),
        path -> preprocess(path, Gisette())
    ]
end
target(::Gisette) = Labels
size(::Gisette) = (6000, 1000, 0)
problem(::Gisette) = Classification

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData