text:
@noinline function old_cfunction(f, r, a)
  ccall(:jl_function_ptr, Ptr{Cvoid}, (Any, Any, Any), f, r, a)
end

## Common Interface Solve Functions

mutable struct CommonFunction{F,P}
    func::F
    p::P
    neq::Cint
end

function commonfun(t::T1,y::T2,yp::T3,comfun::CommonFunction) where {T1,T2,T3}
  y_ = unsafe_wrap(Array,y,comfun.neq)
  ydot_ = unsafe_wrap(Array,yp,comfun.neq)
  comfun.func(ydot_,y_,comfun.p,t)
  return Int32(0)
end

function DiffEqBase.__solve(
    prob::DiffEqBase.AbstractODEProblem{uType,tupType,isinplace},
    alg::LSODAAlgorithm,
    timeseries=[],ts=[],ks=[];

    verbose=true,
    abstol=1/10^6,reltol=1/10^3,
    tstops=Float64[],
    saveat=Float64[], maxiter=Int(1e5),
    callback=nothing,
    timeseries_errors=true,
    save_everystep=isempty(saveat),
    save_start = save_everystep || isempty(saveat) || typeof(saveat) <: Number ? true : prob.tspan[1] in saveat,
    userdata=nothing,
    alias_u0=false,
    kwargs...) where {uType,tupType,isinplace}

    tType = eltype(tupType)

    if verbose
        warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
        if !(typeof(prob.f) <: DiffEqBase.AbstractParameterizedFunction)
            if DiffEqBase.has_tgrad(prob.f)
                @warn("Explicit t-gradient given to this stiff solver is ignored.")
                warned = true
            end
            if DiffEqBase.has_jac(prob.f)
                @warn("Explicit Jacobian given to this stiff solver is ignored.")
                warned = true
            end
        end
        warned && warn_compat()
    end

    if prob.f.mass_matrix != I
        error("This solver is not able to use mass matrices.")
    end

    if callback != nothing || :callback in keys(prob.kwargs)
      error("LSODA is not compatible with callbacks.")
    end

    tspan = prob.tspan
    t0 = tspan[1]
    T = tspan[end]

    if typeof(saveat) <: Number
        if (tspan[1]:saveat:tspan[end])[end] == tspan[end]
          saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:tspan[end]))
        else
          saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:(tspan[end]-saveat)))
        end
    else
        saveat_vec =  convert(Vector{tType},collect(saveat))
    end

    if !isempty(saveat_vec) && saveat_vec[end] == tspan[2]
        pop!(saveat_vec)
    end

    if !isempty(saveat_vec) && saveat_vec[1] == tspan[1]
        save_ts = sort(unique([saveat_vec;T]))
    else
        save_ts = sort(unique([t0;saveat_vec;T]))
    end

    if T < save_ts[end]
        error("Final saving timepoint is past the solving timespan")
    end
    if t0 > save_ts[1]
        error("First saving timepoint is before the solving timespan")
    end

    if !isempty(tstops)
        error("tstops is not supported for this solver. Please use saveat instead")
    end

    if typeof(prob.u0) <: Number
        u0 = [prob.u0]
    else
        if alias_u0
            u0 = vec(prob.u0)
        else
            u0 = vec(deepcopy(prob.u0))
        end
    end

    sizeu = size(prob.u0)

    ### Fix the more general function to Sundials allowed style
    if !isinplace && (typeof(prob.u0)<:Vector{Float64} || typeof(prob.u0)<:Number)
        f! = (du,u,p,t) -> (du[:] = prob.f(u,p,t); nothing)
    elseif !isinplace && typeof(prob.u0)<:AbstractArray
        f! = (du,u,p,t) -> (du[:] = vec(prob.f(reshape(u,sizeu),p,t)); nothing)
    elseif typeof(prob.u0)<:Vector{Float64}
        f! = prob.f
    else # Then it's an in-place function on an abstract array
        f! = (du,u,p,t) -> (prob.f(reshape(du,sizeu),reshape(u,sizeu),p,t); nothing)
    end

    ures = Vector{Float64}[]
    push!(ures,u0)
    utmp = copy(u0)
    utmp2= copy(u0)
    ttmp = [t0]
    t    = [t0]
    t2   = [t0]
    save_start ? ts = [t0] : ts = typeof(t0)[]

    neq = Int32(length(u0))
    comfun = CommonFunction(f!,prob.p,neq)
    atol = ones(Float64,neq)
    rtol = ones(Float64,neq)

    if typeof(abstol) == Float64
        atol *= abstol
    else
        atol = copy(abstol)
    end

    if typeof(reltol) == Float64
        rtol *= reltol
    else
        rtol = copy(reltol)
    end

    GC.@preserve comfun atol rtol begin

    global ___ref = comfun

    opt = lsoda_opt_t(mxstep = maxiter)
    opt.ixpr = 0
    opt.rtol = pointer(rtol)
    opt.atol = pointer(atol)
    if save_everystep
        itask_tmp = 2
    else
        itask_tmp = 1
    end
    opt.itask = itask_tmp

    function get_cfunction(comfun::T) where T
        @cfunction commonfun Cint (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ref{T})
    end

    fex_c = get_cfunction(comfun)

    ctx = lsoda_context_t()
    ctx.function_ = fex_c
    ctx.neq = neq
    ctx.state = 1
    ctx.data = pointer_from_objref(comfun)

    ch = ContextHandle(ctx)

    lsoda_prepare(ctx,opt)

    for k in 2:length(save_ts)
        ttmp[1] = save_ts[k]
        if t[1] < ttmp[1]
            while t[1] < ttmp[1]
                lsoda(ctx, utmp, t, ttmp[1])
                if t[1] > ttmp[1] # overstepd, interpolate back
                    t2[1] = t[1] # save step values
                    copyto!(utmp2,utmp) # save step values
                    opt.itask = 1 # change to interpolating
                    lsoda(ctx, utmp, t, ttmp[1])
                    opt.itask = itask_tmp
                    push!(ures, copy(utmp))
                    push!(ts, t[1])
                    # don't overstep the last timestep
                    if k != length(save_ts) && save_ts[k+1] > t2[1]
                        push!(ures, copy(utmp2))
                        push!(ts, t2[1])
                    end
                    copyto!(utmp, utmp2)
                    t[1] = t2[1]
                else
                    push!(ures, copy(utmp))
                    push!(ts,t[1])
                end
            end
      else
            t2[1] = t[1] # save step values
            copyto!(utmp2, utmp) # save step values
            opt.itask = 1 # change to interpolating
            lsoda(ctx, utmp, t, ttmp[1])
            opt.itask = itask_tmp
            push!(ures, copy(utmp))
            push!(ts, t[1])
            if k != length(save_ts) && save_ts[k+1] > t2[1] # don't overstep the last timestep
                push!(ures,copy(utmp2))
                push!(ts,t2[1])
            end
            copyto!(utmp,utmp2)
            t[1] = t2[1]
        end
    end

    ### Finishing Routine

    timeseries = uType[]
    save_start ? start_idx = 1 : start_idx = 2
    if typeof(prob.u0)<:Number
        for i=start_idx:length(ures)
            push!(timeseries,ures[i][1])
        end
    else
        for i=start_idx:length(ures)
            push!(timeseries,reshape(ures[i],sizeu))
        end
    end

    lsoda_free(ch)
    global ___ref = nothing
    end

    DiffEqBase.build_solution(prob, alg, ts, timeseries,
                   timeseries_errors = timeseries_errors,
                   retcode = :Success)
end


text:
export
    gcpool,
    gdpool,
    georec,
    getelm,
    getfat,
    getfov,
    gfdist,
    gfevnt,
    gffove!,
    gfilum,
    gfocce!,
    gfoclt,
    gfpa,
    gfposc,
    gfrefn,
    gfrepf,
    gfrepi,
    gfrepu,
    gfrfov,
    gfrr,
    gfsep,
    gfsntc,
    gfsstp,
    gfstep,
    gfstol,
    gfsubc,
    gftfov,
    gfudb!,
    gfuds!,
    gipool,
    gnpool

# Dummy routine for GF functions
function udf(et::SpiceDouble, value::Ptr{SpiceDouble})
    ccall((:udf_c, libcspice), Cvoid, (SpiceDouble, Ptr{SpiceDouble}), et, value)
end

"""
    gcpool(name; start=1, room=100, lenout=128)

Return the value of a kernel variable from the kernel pool.

### Arguments ###

- `name`: Name of the variable whose value is to be returned
- `start`: Which component to start retrieving for name (default: 1)
- `room`: The largest number of values to return (default: 100)
- `lenout`: The length of the longest string to return (default: 128)

### Output ###

Returns an array of values if the variable exists or `nothing` if not.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gcpool_c.html)
"""
function gcpool(name; start = 1, room = 100, lenout = 128)
    n = Ref{SpiceInt}()
    values = Array{SpiceChar}(undef, lenout, room)
    found = Ref{SpiceInt}()
    ccall((:gcpool_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceChar}, Ref{SpiceBoolean}),
          name, start - 1, room, lenout, n, values, found)
    handleerror()
    if Bool(found[])
        return chararray_to_string(values, n[])
    end
    nothing
end

"""
    gdpool(name; start=1, room=100)

Return the value of a kernel variable from the kernel pool.

### Arguments ###

- `name`: Name of the variable whose value is to be returned
- `start`: Which component to start retrieving for name (default: 1)
- `room`: The largest number of values to return (default: 100)

### Output ###

Returns an array of values if the variable exists or `nothing` if not.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gdpool_c.html)
"""
function gdpool(name; start = 1, room = 100)
    n = Ref{SpiceInt}()
    values = Array{SpiceDouble}(undef, room)
    found = Ref{SpiceInt}()
    ccall((:gdpool_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceDouble}, Ref{SpiceBoolean}),
          name, start - 1, room, n, values, found)
    handleerror()
    Bool(found[]) ? values[1:n[]] : nothing
end

"""
    georec(lon, lat, alt, re, f)

Convert geodetic coordinates to rectangular coordinates.

### Arguments ###

- `lon`: Geodetic longitude of point (radians)
- `lat`: Geodetic latitude  of point (radians)
- `alt`: Altitude of point above the reference spheroid
- `re`: Equatorial radius of the reference spheroid
- `f`: Flattening coefficient

### Output ###

Returns the rectangular coordinates of point.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/georec_c.html)
"""
function georec(lon, lat, alt, re, f)
    rectan = Array{SpiceDouble}(undef, 3)
    ccall((:georec_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, SpiceDouble, SpiceDouble, SpiceDouble,
           Ref{SpiceDouble}),
          lon, lat, alt, re, f, rectan)
    handleerror()
    rectan
end

"""
    getelm(frstyr, lines)

Given the "lines" of a two-line element set, parse the lines and return the elements in units
suitable for use in SPICE software.

### Arguments ###

- `frstyr`: Year of earliest representable two-line elements
- `lines`: A pair of "lines" containing two-line elements

### Output ###

- `epoch`: The epoch of the elements in seconds past J2000
- `elems`: The elements converted to SPICE units

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/getelm_c.html)
"""
function getelm(frstyr, lines)
    lines_, _, lineln = chararray(lines)
    epoch = Ref{SpiceDouble}()
    elems = Array{SpiceDouble}(undef, 10)
    ccall((:getelm_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, Ref{SpiceChar}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          frstyr, lineln, lines_, epoch, elems)
    handleerror()
    epoch[], elems
end

"""
    getfat(file, arclen=10, typlen=10)

Determine the file architecture and file type of most SPICE kernel files.

### Arguments ###

- `file`: The name of a file to be examined
- `arclen`: Maximum length of output architecture string (default: 10)
- `typlen`: Maximum length of output type string (default: 10)

### Output ###

- `arch`: The architecture of the kernel file
- `typ`: The type of the kernel file

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/getfat_c.html)
"""
function getfat(file, arclen = 10, typlen = 10)
    arch = Array{SpiceChar}(undef, arclen)
    typ = Array{SpiceChar}(undef, typlen)
    ccall((:getfat_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}, Ref{SpiceChar}),
          file, arclen, typlen, arch, typ)
    handleerror()
    chararray_to_string(arch), chararray_to_string(typ)
end

"""
    getfov(instid, room=10, shapelen=128, framelen=128)

Return the field-of-view (FOV) parameters for a specified instrument.
The instrument is specified by its NAIF ID code.

### Arguments ###

- `instid  `: NAIF ID of an instrument
- `room    `: Maximum number of vectors that can be returned (default: 10)
- `shapelen`: Space available in the string `shape` (default: 128)
- `framelen`: Space available in the string `frame` (default: 128)

### Output ###

Returns a tuple consisting of

- `shape `: Instrument FOV shape
- `frame `: Name of the frame in which FOV vectors are defined
- `bsight`: Boresight vector
- `bounds`: FOV boundary vectors

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/getfov_c.html)
"""
function getfov(instid, room = 10, shapelen = 128, framelen = 128)
    shape = Array{SpiceChar}(undef, shapelen)
    frame = Array{SpiceChar}(undef, framelen)
    bsight = Array{SpiceDouble}(undef, 3)
    n = Ref{SpiceInt}()
    bounds = Array{SpiceDouble}(undef, 3, room)
    ccall((:getfov_c, libcspice), Cvoid,
          (SpiceInt, SpiceInt, SpiceInt, SpiceInt,
           Ref{SpiceChar}, Ref{SpiceChar}, Ref{SpiceDouble}, Ref{SpiceInt}, Ref{SpiceDouble}),
          instid, room, shapelen, framelen, shape, frame, bsight, n, bounds)
    handleerror()
    arr_bounds = cmatrix_to_array(bounds)
    shp = chararray_to_string(shape)
    frm = chararray_to_string(frame)
    shp, frm, bsight, arr_bounds[1:n[]]
end

function gfbail()
    ccall((:gfbail_c, libcspice), SpiceBoolean, ())
end

function gfclrh()
    ccall((:gfclrh_c, libcspice), Cvoid, ())
end

"""
    gfdist(target, abcorr, obsrvr, relate, refval, adjust, step, nintvls, cnfine)

Return the time window over which a specified constraint on observer-target distance is met.

### Arguments ###

- `target`: Name of the target body
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `relate`: Relational operator
- `refval`: Reference value
- `adjust`: Adjustment value for absolute extrema searches
- `step`: Step size used for locating extrema and roots
- `nintvls`: Workspace window interval count
- `cnfine`: Window to which the search is confined

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfdist_c.html)
"""
function gfdist(target, abcorr, obsrvr, relate, refval, adjust, step, nintvls, cnfine)
    result = SpiceDoubleCell(nintvls * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfdist_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, SpiceDouble, SpiceDouble, SpiceDouble,
               SpiceInt, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              target, abcorr, obsrvr, relate, refval, adjust, step, nintvls, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfevnt(udstep, udrefn, gquant, qnpars, lenvals, qpnams, qcpars, qdpars, qipars, qlpars,
           op, refval, tol, adjust, rpt, udrepi, udrepu, udrepf, nintvls, bail, udbail, cnfine)

Determine time intervals when a specified geometric quantity satisfies a specified mathematical
condition.

### Arguments ###

- `udstep`: Name of the routine that computes and returns a time step
- `udrefn`: Name of the routine that computes a refined time
- `gquant`: Type of geometric quantity
- `qpnams`: Names of quantity definition parameters
- `qcpars`: Array of character quantity definition parameters
- `qdpars`: Array of double precision quantity definition parameters
- `qipars`: Array of integer quantity definition parameters
- `qlpars`: Array of logical quantity definition parameters
- `op`: Operator that either looks for an extreme value (max, min, local, absolute) or compares the
    geometric quantity value and a number
- `refval`: Reference value
- `tol`: Convergence tolerance in second
- `adjust`: Absolute extremum adjustment value
- `rpt`: Progress reporter on (`true`) or off (`false`)
- `udrepi`: Function that initializes progress reporting
- `udrepu`: Function that updates the progress report
- `udrepf`: Function that finalizes progress reporting
- `nintvls`: Workspace window interval coun
- `cnfine`: SPICE window to which the search is restricted

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfevnt_c.html)
"""
function gfevnt(udstep, udrefn, gquant, qpnams, qcpars, qdpars, qipars, qlpars,
                op, refval, tol, adjust, rpt, udrepi, udrepu, udrepf, nintvls, cnfine)
    result = SpiceDoubleCell(nintvls * 2)
    _qpnams, qnpars, lenvals = chararray(qpnams)
    _qcpars, _, lenvals_qc = chararray(qcpars)
    lenvals != lenvals_qc &&
        throw(SpiceError("The strings in `qpnams` and `qcpars` must all have the same length"))
    @checkdims qnpars qcpars qdpars qipars qlpars
    _qipars = SpiceInt.(qipars)
    _qlpars = SpiceBoolean.(qlpars)
    function _udstep(et::SpiceDouble, step::Ptr{SpiceDouble})
        _step = unsafe_wrap(Array{SpiceDouble}, step, 1)
        _step[1] = udstep(et)
        nothing
    end
    udstep_ptr = @cfunction($_udstep, Cvoid, (SpiceDouble, Ptr{SpiceDouble}))
    function _udrefn(t1::SpiceDouble, t2::SpiceDouble, s1::SpiceBoolean, s2::SpiceBoolean,
                     t::Ptr{SpiceDouble})
        _t = unsafe_wrap(Array{SpiceDouble}, t, 1)
        _t[1] = udrefn(t1, t2, s1, s2)
        nothing
    end
    udrefn_ptr = @cfunction($_udrefn, Cvoid,
                            (SpiceDouble, SpiceDouble, SpiceBoolean, SpiceBoolean,
                             Ptr{SpiceDouble}))
    udrepf_ptr = @cfunction($udrepf, Cvoid, ())
    function _udrepi(window::Ptr{Cell{SpiceDouble}}, begmss::Cstring, endmss::Cstring)
        _window = unsafe_load(window, 1)
        udrepi(_window, begmss, endmss)
    end
    udrepi_ptr = @cfunction($_udrepi, Cvoid, (Ptr{Cell{SpiceDouble}}, Cstring, Cstring))
    udrepu_ptr = @cfunction($udrepu, Cvoid, (SpiceDouble, SpiceDouble, SpiceDouble))
    # Julia should handle SIGINT
    bail = false
    udbail = @cfunction(gfbail, SpiceBoolean, ())
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfevnt_c, libcspice), Cvoid,
              (Ptr{Cvoid}, Ptr{Cvoid}, Cstring, SpiceInt, SpiceInt, Ref{SpiceChar}, Ref{SpiceChar},
               Ref{SpiceDouble}, Ref{SpiceInt}, Ref{SpiceBoolean}, Cstring, SpiceDouble, SpiceDouble,
               SpiceDouble, SpiceBoolean, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, SpiceInt, SpiceBoolean,
               Ptr{Cvoid}, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              udstep_ptr, udrefn_ptr, gquant, qnpars, lenvals, _qpnams, _qcpars, qdpars, _qipars,
              _qlpars, op, refval, tol, adjust, rpt, udrepi_ptr, udrepu_ptr, udrepf_ptr, nintvls,
              bail, udbail, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gffove!(inst, tshape, raydir, target, tframe, abcorr, obsrvr, tol, udstep, udrefn,
                rpt, udrepi, udrepu, udrepf, cnfine, result)

Determine time intervals when a specified target body or ray intersects the space bounded by the
field-of-view (FOV) of a specified instrument.

### Arguments ###

- `inst`: Name of the instrument
- `tshape`: Type of shape model used for target body
- `raydir`: Ray's direction vector
- `target`: Name of the target body
- `tframe`: Body-fixed, body-centered frame for target body
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `tol`: Convergence tolerance in seconds
- `udstep`: Name of the routine returns a time step
- `udrefn`: Name of the routine that computes a refined time
- `rpt`: Progress report flag
- `udrepi`: Function that initializes progress reporting
- `udrepu`: Function that updates the progress report
- `udrepf`: Function that finalizes progress reporting
- `cnfine`: SPICE window to which the search is restricted
- `result`: Window containing the results

### Output ###

Returns `result`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gffove_c.html)
"""
function gffove!(inst, tshape, raydir, target, tframe, abcorr, obsrvr, tol, udstep, udrefn,
                rpt, udrepi, udrepu, udrepf, cnfine, result)
    @checkdims 3 raydir
    function _udstep(et::SpiceDouble, step::Ptr{SpiceDouble})
        _step = unsafe_wrap(Array{SpiceDouble}, step, 1)
        _step[1] = udstep(et)
        nothing
    end
    udstep_ptr = @cfunction($_udstep, Cvoid, (SpiceDouble, Ptr{SpiceDouble}))
    function _udrefn(t1::SpiceDouble, t2::SpiceDouble, s1::SpiceBoolean, s2::SpiceBoolean,
                     t::Ptr{SpiceDouble})
        _t = unsafe_wrap(Array{SpiceDouble}, t, 1)
        _t[1] = udrefn(t1, t2, s1, s2)
        nothing
    end
    udrefn_ptr = @cfunction($_udrefn, Cvoid,
                            (SpiceDouble, SpiceDouble, SpiceBoolean, SpiceBoolean,
                             Ptr{SpiceDouble}))
    udrepf_ptr = @cfunction($udrepf, Cvoid, ())
    function _udrepi(window::Ptr{Cell{SpiceDouble}}, begmss::Cstring, endmss::Cstring)
        _window = unsafe_load(window, 1)
        udrepi(_window, begmss, endmss)
    end
    udrepi_ptr = @cfunction($_udrepi, Cvoid, (Ptr{Cell{SpiceDouble}}, Cstring, Cstring))
    udrepu_ptr = @cfunction($udrepu, Cvoid, (SpiceDouble, SpiceDouble, SpiceDouble))
    # Julia should handle SIGINT
    bail = false
    udbail = @cfunction(gfbail, SpiceBoolean, ())
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gffove_c, libcspice), Cvoid,
              (Cstring, Cstring, Ref{SpiceDouble}, Cstring, Cstring, Cstring, Cstring, SpiceDouble,
               Ptr{Cvoid}, Ptr{Cvoid}, SpiceBoolean, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, SpiceBoolean,
               Ptr{Cvoid}, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              inst, tshape, raydir, target, tframe, abcorr, obsrvr, tol, udstep_ptr, udrefn_ptr, rpt,
              udrepi_ptr, udrepu_ptr, udrepf_ptr, bail, udbail, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfilum(method, angtyp, target, illmn, fixref, abcorr, obsrvr, spoint, relate, refval,
           adjust, step, nintvls, cnfine)

Return the time window over which a specified constraint on the observed phase, solar incidence,
or emission angle at a specifed target body surface point is met.

### Arguments ###

- `method`: Computation method
- `angtyp`: Type of illumination angle
- `target`: Name of the target body
- `illmn `: Name of the illumination source
- `fixref`: Body-fixed, body-centered target body frame
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `spoint`: Body-fixed coordinates of a target surface point
- `relate`: Relational operator
- `refval`: Reference value
- `adjust`: Adjustment value for absolute extrema searches
- `step`: Step size used for locating extrema and roots
- `nintvls`: Workspace window interval count
- `cnfine`: Window to which the search is confined

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfilum_c.html)
"""
function gfilum(method, angtyp, target, illmn, fixref, abcorr, obsrvr, spoint, relate, refval,
                adjust, step, nintvls, cnfine)
    @checkdims 3 spoint
    result = SpiceDoubleCell(nintvls * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfilum_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Ref{SpiceDouble},
               Cstring, SpiceDouble, SpiceDouble, SpiceDouble, SpiceInt,
               Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              method, angtyp, target, illmn, fixref, abcorr, obsrvr, spoint, relate, refval, adjust,
              step, nintvls, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    function gfocce!(occtyp, front, fshape, fframe, back, bshape, bframe, abcorr, obsrvr, tol,
        udstep, udrefn, rpt, udrepi, udrepu, udrepf, cnfine, result)

Determine time intervals when an observer sees one target occulted by another.

### Arguments ###

- `occtyp`: Type of occultation
- `front`: Name of body occulting the other
- `fshape`: Type of shape model used for front body
- `fframe`: Body-fixed, body-centered frame for front body
- `back`: Name of body occulted by the other
- `bshape`: Type of shape model used for back body
- `bframe`: Body-fixed, body-centered frame for back body
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `tol`: Convergence tolerance in seconds
- `udstep`: Name of the routine that returns a time step
- `udrefn`: Name of the routine that computes a refined time
- `rpt`: Progress report flag
- `udrepi`: Function that initializes progress reporting
- `udrepu`: Function that updates the progress report
- `udrepf`: Function that finalizes progress reporting
- `cnfine`: SPICE window to which the search is restricted
- `result`: SPICE window containing results

### Output ###

Returns `result`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfocce_c.html)
"""
function gfocce!(occtyp, front, fshape, fframe, back, bshape, bframe, abcorr, obsrvr, tol, udstep, udrefn, rpt, udrepi, udrepu, udrepf, cnfine, result)
    function _udstep(et::SpiceDouble, step::Ptr{SpiceDouble})
        _step = unsafe_wrap(Array{SpiceDouble}, step, 1)
        _step[1] = udstep(et)
        nothing
    end
    udstep_ptr = @cfunction($_udstep, Cvoid, (SpiceDouble, Ptr{SpiceDouble}))
    function _udrefn(t1::SpiceDouble, t2::SpiceDouble, s1::SpiceBoolean, s2::SpiceBoolean,
                     t::Ptr{SpiceDouble})
        _t = unsafe_wrap(Array{SpiceDouble}, t, 1)
        _t[1] = udrefn(t1, t2, s1, s2)
        nothing
    end
    udrefn_ptr = @cfunction($_udrefn, Cvoid,
                            (SpiceDouble, SpiceDouble, SpiceBoolean, SpiceBoolean,
                             Ptr{SpiceDouble}))
    udrepf_ptr = @cfunction($udrepf, Cvoid, ())
    function _udrepi(window::Ptr{Cell{SpiceDouble}}, begmss::Cstring, endmss::Cstring)
        _window = unsafe_load(window, 1)
        udrepi(_window, begmss, endmss)
    end
    udrepi_ptr = @cfunction($_udrepi, Cvoid, (Ptr{Cell{SpiceDouble}}, Cstring, Cstring))
    udrepu_ptr = @cfunction($udrepu, Cvoid, (SpiceDouble, SpiceDouble, SpiceDouble))
    # Julia should handle SIGINT
    bail = false
    udbail = @cfunction(gfbail, SpiceBoolean, ())
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfocce_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, SpiceDouble,
               Ptr{Cvoid}, Ptr{Cvoid}, SpiceBoolean, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, SpiceBoolean,
               Ptr{Cvoid}, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              occtyp, front, fshape, fframe, back, bshape, bframe, abcorr, obsrvr, tol, udstep_ptr,
              udrefn_ptr, rpt, udrepi_ptr, udrepu_ptr, udrepf_ptr, bail, udbail,
              cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfoclt(occtyp, front, fshape, fframe, back, bshape, bframe, abcorr, obsrvr, step, cnfine,
           maxwin=100)

Determine time intervals when an observer sees one target occulted by, or in transit across,
another.

The surfaces of the target bodies may be represented by triaxial ellipsoids or by topographic data
provided by DSK files.

### Arguments ###

- `occtyp`: Type of occultation
- `front`: Name of body occulting the other
- `fshape`: Type of shape model used for front body
- `fframe`: Body-fixed, body-centered frame for front body
- `back`: Name of body occulted by the other
- `bshape`: Type of shape model used for back body
- `bframe`: Body-fixed, body-centered frame for back body
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `step`: Step size in seconds for finding occultation events
- `cnfine`: Window to which the search is restricted
- `maxwin`: Maximum size of the output window (default: 100)

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfoclt_c.html)
"""
function gfoclt(occtyp, front, fshape, fframe, back, bshape, bframe, abcorr, obsrvr, step, cnfine, maxwin = 100)
    result = SpiceDoubleCell(maxwin * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfoclt_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring,
               SpiceDouble, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              occtyp, front, fshape, fframe, back, bshape, bframe, abcorr, obsrvr, step, cnfine.cell,
              result.cell)
    end
    handleerror()
    result
end

"""
    gfpa(result, target, illmn, abcorr, obsrvr, relate, refval, adjust, step, nintvls, cnfine)

Determine time intervals for which a specified constraint on the phase angle between an
illumination source, a target, and observer body centers is met.

### Arguments ###

- `target`: Name of the target body
- `illmn`: Name of the illuminating body
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `relate`: Relational operator
- `refval`: Reference value
- `adjust`: Adjustment value for absolute extrema searches
- `step`: Step size used for locating extrema and roots
- `nintvls`: Workspace window interval count
- `cnfine`: Window to which the search is confined

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfpa_c.html)
"""
function gfpa(target, illmn, abcorr, obsrvr, relate, refval, adjust, step, nintvls, cnfine)
    result = SpiceDoubleCell(nintvls * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfpa_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, Cstring,
               SpiceDouble, SpiceDouble, SpiceDouble,
               SpiceInt,
               Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              target, illmn, abcorr, obsrvr, relate, refval, adjust, step, nintvls, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfposc(target, frame, abcorr, obsrvr, crdsys, coord, relate, refval, adjust, step,
           nintvls, cnfine)

Determine time intervals for which a coordinate of an observer-target position vector satisfies a
numerical constraint.

### Arguments ###

- `target`: Name of the target body
- `frame`: Name of the reference frame for coordinate calculations
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `crdsys`: Name of the coordinate system containing `coord`
- `coord`: Name of the coordinate of interest
- `relate`: Operator that either looks for an extreme value (max, min, local, absolute) or compares
    the coordinate value and refval
- `refval`: Reference value
- `adjust`: Adjustment value for absolute extrema searches
- `step`: Step size used for locating extrema and roots
- `nintvls`: Workspace window interval count
- `cnfine`: Window to which the search is restricted

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfposc_c.html)
"""
function gfposc(target, frame, abcorr, obsrvr, crdsys, coord, relate, refval, adjust, step,
                nintvls, cnfine)
    result = SpiceDoubleCell(nintvls * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfposc_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, SpiceDouble, SpiceDouble,
               SpiceDouble, SpiceInt, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              target, frame, abcorr, obsrvr, crdsys, coord, relate, refval, adjust, step, nintvls,
              cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfrefn(t1, t2, s1, s2)

For those times when we can't do better, we use a bisection method to find the next time at which
to test for state change.

### Arguments ###

- `t1`: One of two values bracketing a state change
- `t2`: The other value that brackets a state change
- `s1`: State at `t1`
- `s2`: State at `t2`

### Output ###

Returns the new value at which to check for transition.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfrefn_c.html)
"""
function gfrefn(t1, t2, s1, s2)
    t = Ref{SpiceDouble}()
    ccall((:gfrefn_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, SpiceBoolean, SpiceBoolean, Ref{SpiceDouble}),
          t1, t2, s1, s2, t)
    t[]
end

"""
    gfrepf()

Finish a GF progress report.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfrepf_c.html)
"""
function gfrepf()
    ccall((:gfrepf_c, libcspice), Cvoid, ())
    handleerror()
end

"""
    gfrepi(window, begmss, endmss)

Initialize a search progress report.

### Arguments ###

- `window`: A window over which a job is to be performed
- `begmss`: Beginning of the text portion of the output message
- `endmss`: End of the text portion of the output message

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfrepi_c.html)
"""
function gfrepi(window::SpiceDoubleCell, begmss, endmss)
    data = window.data
    GC.@preserve data begin
        ccall((:gfrepi_c, libcspice), Cvoid,
              (Ref{Cell{SpiceDouble}}, Cstring, Cstring),
              window.cell, begmss, endmss)
    end
    handleerror()
end
function gfrepi(window::Cell{SpiceDouble}, begmss, endmss)
    ccall((:gfrepi_c, libcspice), Cvoid,
          (Ref{Cell{SpiceDouble}}, Cstring, Cstring),
          window, begmss, endmss)
    handleerror()
end

"""
    gfrepu(ivbeg, ivend, time)

Tell the progress reporting system how far a search has progressed.

### Arguments ###

- `ivbeg`: Start time of work interval
- `ivend`: End time of work interval
- `time `: Current time being examined in the search process

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfrepu_c.html)
"""
function gfrepu(ivbeg, ivend, time)
    ccall((:gfrepu_c, libcspice), Cvoid,
          (SpiceDouble, SpiceDouble, SpiceDouble),
          ivbeg, ivend, time)
    handleerror()
end

"""
    gfrfov(inst, raydir, rframe, abcorr, obsrvr, step, cnfine, maxwin=10000)

Determine time intervals when a specified ray intersects the space bounded by the
field-of-view (FOV) of a specified instrument.

### Arguments ###

- `inst`: Name of the instrument
- `raydir`: Ray's direction vector
- `rframe`: Reference frame of ray's direction vector
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `step`: Step size in seconds for finding FOV events
- `cnfine`: SPICE window to which the search is restricted
- `maxwin`: Maximum length of the output window (default: 10000)

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfrfov_c.html)
"""
function gfrfov(inst, raydir, rframe, abcorr, obsrvr, step, cnfine, maxwin = 10000)
    @checkdims 3 raydir
    result = SpiceDoubleCell(maxwin)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfrfov_c, libcspice), Cvoid,
              (Cstring, Ref{SpiceDouble}, Cstring, Cstring, Cstring, SpiceDouble,
               Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              inst, raydir, rframe, abcorr, obsrvr, step, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfrr(target, abcorr, obsrvr, relate, refval, adjust, step, nintvls, cnfine)

Determine time intervals for which a specified constraint on the observer-target range rate is met.

### Arguments ###

- `target`: Name of the target body
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `relate`: Relational operator
- `refval`: Reference value
- `adjust`: Adjustment value for absolute extrema searches
- `step`: Step size used for locating extrema and roots
- `nintvls`: Workspace window interval count
- `cnfine`: Window to which the search is confined

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfrr_c.html)
"""
function gfrr(target, abcorr, obsrvr, relate, refval, adjust, step, nintvls, cnfine)
    result = SpiceDoubleCell(nintvls * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfrr_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, SpiceDouble, SpiceDouble, SpiceDouble, SpiceInt,
               Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              target, abcorr, obsrvr, relate, refval, adjust, step, nintvls, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
Determine time intervals when the angular separation between the position vectors of two target
bodies relative to an observer satisfies a numerical relationship.

### Arguments ###

- `targ1`: Name of first body
- `shape1`: Name of shape model describing the first body
- `frame1`: The body-fixed reference frame of the first body
- `targ2`: Name of second body
- `shape2`: Name of the shape model describing the second body
- `frame2`: The body-fixed reference frame of the second body
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `relate`: Operator that either looks for an extreme value (max, min, local, absolute) or compares
    the angular separation value and refval
- `refval`: Reference value
- `adjust`: Absolute extremum adjustment value
- `step`: Step size in seconds for finding angular separation events
- `nintvls`: Workspace window interval count
- `cnfine`: Window to which the search is restricted

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfsep_c.html)
"""
function gfsep(targ1, shape1, frame1, targ2, shape2, frame2, abcorr, obsrvr, relate, refval, adjust,
               step, nintvls, cnfine)
    result = SpiceDoubleCell(nintvls * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfsep_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring,
               SpiceDouble, SpiceDouble, SpiceDouble, SpiceInt,
               Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              targ1, shape1, frame1, targ2, shape2, frame2, abcorr, obsrvr, relate, refval, adjust,
              step, nintvls, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfsntc(target, fixref, method, abcorr, obsrvr, dref, dvec, crdsys, coord, relate, refval,
           adjust, step, nintvls, cnfine)

Determine time intervals for which a coordinate of an surface intercept position vector satisfies
a numerical constraint.

### Arguments ###

- `target`: Name of the target body
- `fixref`: Body fixed frame associated with `target`
- `method`: Name of method type for surface intercept calculation
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `dref`: Reference frame of direction vector `dvec`
- `dvec`: Pointing direction vector from `obsrvr`
- `crdsys`: Name of the coordinate system containing COORD
- `coord`: Name of the coordinate of interest
- `relate`: Operator that either looks for an extreme value (max, min, local, absolute) or compares the coordinate value and `refval`
- `refval`: Reference value
- `adjust`: Adjustment value for absolute extrema searches
- `step`: Step size used for locating extrema and roots
- `nintvls`: Workspace window interval count
- `cnfine`: Window to which the search is restricted

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfsntc_c.html)
"""
function gfsntc(target, fixref, method, abcorr, obsrvr, dref, dvec, crdsys, coord, relate, refval,
                adjust, step, nintvls, cnfine)
    @checkdims 3 dvec
    result = SpiceDoubleCell(nintvls * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfsntc_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Ref{SpiceDouble}, Cstring,
               Cstring, Cstring, SpiceDouble, SpiceDouble, SpiceDouble, SpiceInt,
               Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              target, fixref, method, abcorr, obsrvr, dref, dvec, crdsys, coord, relate, refval,
              adjust, step, nintvls, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfsstp(step)

Set the step size to be returned by [`gfstep`](@ref).

### Arguments ###

- `step`: Time step to take

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfsstp_c.html)
"""
function gfsstp(step)
    ccall((:gfsstp_c, libcspice), Cvoid, (SpiceDouble,), step)
    handleerror()
end

"""
    gfstep()

Return the time step set by the most recent call to [`gfsstp`](@ref).

### Arguments ###

- `step`: Time step to take

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfstep_c.html)
"""
function gfstep()
    step = Ref{SpiceDouble}()
    ccall((:gfstep_c, libcspice), Cvoid,
          (SpiceDouble, Ref{SpiceDouble}),
          0.0, step)
    handleerror()
    step[]
end
gfstep(_) = gfstep()

"""
    gfstol(value)

Override the default GF convergence value used in the high level GF routines.

### Arguments ###

- `value`: Double precision value returned or to store

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfstol_c.html)
"""
function gfstol(value)
    ccall((:gfstol_c, libcspice), Cvoid, (SpiceDouble,), value)
    handleerror()
end

"""
    gfsubc(target, fixref, method, abcorr, obsrvr, crdsys, coord, relate, refval, adjust, step,
           nintvls, cnfine)

Determine time intervals for which a coordinate of an subpoint position vector satisfies a
numerical constraint.

### Arguments ###

- `target`: Name of the target body
- `fixref`: Body fixed frame associated with `target`
- `method`: Name of method type for subpoint calculation
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `crdsys`: Name of the coordinate system containing `coord`
- `coord`: Name of the coordinate of interest
- `relate`: Operator that either looks for an extreme value (max, min, local, absolute) or compares
    the coordinate value and refval
- `refval`: Reference value
- `adjust`: Adjustment value for absolute extrema searches
- `step`: Step size used for locating extrema and roots
- `nintvls`: Workspace window interval count
- `cnfine`: Window to which the search is restricted

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfsubc_c.html)
"""
function gfsubc(target, fixref, method, abcorr, obsrvr, crdsys, coord, relate, refval, adjust, step,
                nintvls, cnfine)
    result = SpiceDoubleCell(nintvls * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfsubc_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, SpiceDouble,
               SpiceDouble, SpiceDouble, SpiceInt, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              target, fixref, method, abcorr, obsrvr, crdsys, coord, relate, refval, adjust, step,
              nintvls, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gftfov(inst, target, tshape, tframe, abcorr, obsrvr, step, nintvls, cnfine)

Determine time intervals when a specified ephemeris object intersects the space bounded by the
field-of-view (FOV) of a specified instrument.

### Arguments ###

- `inst`: Name of the instrument
- `target`: Name of the target body
- `tshape`: Type of shape model used for target body
- `tframe`: Body-fixed, body-centered frame for target body
- `abcorr`: Aberration correction flag
- `obsrvr`: Name of the observing body
- `step`: Step size in seconds for finding FOV events
- `nintvls`: Workspace window interval count
- `cnfine`: Window to which the search is restricted

### Output ###

Returns a window containing the results.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gftfov_c.html)
"""
function gftfov(inst, target, tshape, tframe, abcorr, obsrvr, step, nintvls, cnfine)
    result = SpiceDoubleCell(nintvls * 2)
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gftfov_c, libcspice), Cvoid,
              (Cstring, Cstring, Cstring, Cstring, Cstring, Cstring, SpiceDouble,
               Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              inst, target, tshape, tframe, abcorr, obsrvr, step, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfudb!(udfuns, udfunb, step, cnfine, result)

Perform a GF search on a user defined boolean quantity.

### Arguments ###

- `udfuns`: Name of the routine that computes a scalar quantity of interest corresponding to an `et`, e.g. `f(et) = ...`
- `udfunb`: Name of the routine returning the boolean value corresponding to an `et`, e.g. `g(f, et) = ...`
- `step  `: Step size used for locating extrema and roots
- `cnfine`: Window to which the search is restricted
- `result`: Window containing results

### Output ###

Returns `result`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfudb_c.html)
"""
function gfudb!(udfuns, udfunb, step, cnfine, result)
    function _udfuns(et::SpiceDouble, value::Ptr{SpiceDouble})
        _value = unsafe_wrap(Array{SpiceDouble}, value, 1)
        _value[1] = udfuns(et)
        nothing
    end
    udfuns_ptr = @cfunction($_udfuns, Cvoid, (SpiceDouble, Ptr{SpiceDouble}))
    function _udfunb(_::Ptr{Cvoid}, et::SpiceDouble, value::Ptr{SpiceBoolean})
        _value = unsafe_wrap(Array{SpiceBoolean}, value, 1)
        _value[1] = udfunb(udfuns, et)
        nothing
    end
    udfunb_ptr = @cfunction($_udfunb, Cvoid, (Ptr{Cvoid}, SpiceDouble, Ptr{SpiceBoolean}))
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfudb_c, libcspice), Cvoid,
              (Ptr{Cvoid}, Ptr{Cvoid}, SpiceDouble, Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              udfuns_ptr, udfunb_ptr, step, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gfuds!(udfuns, udqdec, relate, refval, adjust, step, nintvls, cnfine, result)

Perform a GF search on a user defined scalar quantity.

### Arguments ###

- `udfuns`: Name of the routine that computes the scalar quantity of interest at some time, e.g. `f(et) = ...`
- `udqdec`: Name of the routine that computes whether the scalar quantity is decreasing, e.g. `g(f, et) = ...`
- `relate`: Operator that either looks for an extreme value (max, min, local, absolute) or compares the geometric quantity value and a number
- `refval`: Value used as reference for scalar quantity condition
- `adjust`: Allowed variation for absolute extremal geometric conditions
- `step`: Step size used for locating extrema and roots
- `nintvls`: Workspace window interval count
- `cnfine`: SPICE window to which the search is restricted
- `result`: SPICE window containing results

### Output ###

Returns `result`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gfuds_c.html)
"""
function gfuds!(udfuns, udqdec, relate, refval, adjust, step, nintvls, cnfine, result)
    function _udfuns(et::SpiceDouble, value::Ptr{SpiceDouble})
        _value = unsafe_wrap(Array{SpiceDouble}, value, 1)
        _value[1] = udfuns(et)
        nothing
    end
    udfuns_ptr = @cfunction($_udfuns, Cvoid, (SpiceDouble, Ptr{SpiceDouble}))
    function _udqdec(_::Ptr{Cvoid}, et::SpiceDouble, value::Ptr{SpiceBoolean})
        _value = unsafe_wrap(Array{SpiceBoolean}, value, 1)
        _value[1] = udqdec(udfuns, et)
        nothing
    end
    udqdec_ptr = @cfunction($_udqdec, Cvoid, (Ptr{Cvoid}, SpiceDouble, Ptr{SpiceBoolean}))
    cnfine_data = cnfine.data
    result_data = result.data
    GC.@preserve cnfine_data result_data begin
        ccall((:gfuds_c, libcspice), Cvoid,
              (Ptr{Cvoid}, Ptr{Cvoid}, Cstring, SpiceDouble, SpiceDouble, SpiceDouble, SpiceInt,
               Ref{Cell{SpiceDouble}}, Ref{Cell{SpiceDouble}}),
              udfuns_ptr, udqdec_ptr, relate, refval, adjust, step, nintvls, cnfine.cell, result.cell)
    end
    handleerror()
    result
end

"""
    gipool(name; start=1, room=100)

Return the value of a kernel variable from the kernel pool.

### Arguments ###

- `name`: Name of the variable whose value is to be returned
- `start`: Which component to start retrieving for name (default: 1)
- `room`: The largest number of values to return (default: 100)

### Output ###

Returns an array of values if the variable exists or `nothing` if not.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gipool_c.html)
"""
function gipool(name; start = 1, room = 100)
    n = Ref{SpiceInt}()
    values = Array{SpiceInt}(undef, room)
    found = Ref{SpiceBoolean}()
    ccall((:gipool_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceInt, Ref{SpiceInt}, Ref{SpiceInt}, Ref{SpiceBoolean}),
          name, start - 1, room, n, values, found)
    handleerror()
    Bool(found[]) ? Int.(values[1:n[]]) : nothing
end

"""
    gnpool(name, start, room, lenout=128)

Return names of kernel variables matching a specified template.

### Arguments ###

- `name`: Template that names should match
- `start`: Index of first matching name to retrieve
- `room`: The largest number of values to return
- `lenout`: Length of strings in output array `kvars` (default: 128)

### Output ###

Returns lernel pool variables whose names match `name`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/gnpool_c.html)
"""
function gnpool(name, start, room, lenout = 128)
    n = Ref{SpiceInt}()
    kvars = Array{SpiceChar}(undef, lenout, room)
    found = Ref{SpiceBoolean}()
    ccall((:gnpool_c, libcspice), Cvoid,
          (Cstring, SpiceInt, SpiceInt, SpiceInt,
           Ref{SpiceInt}, Ref{SpiceChar}, Ref{SpiceBoolean}),
          name, start - 1, room, lenout, n, kvars, found)
    handleerror()
    Bool(found[]) || return nothing
    chararray_to_string(kvars)[1:n[]]
end

ext: jl


text:
export Newlbfgs

"""
A globalized L-BFGS algorithm with Line Search
Inputs:
    - An AbstractNLPModel (our problem)
    - An NLPStopping (stopping fonctions/criterion of our problem)
    - (opt) verbose = true will show function value, norm of gradient and the
        inner product of the gradient and the descent direction
    - (opt) hessian_rep, different way to implement the Hessian. The default
        way is a LinearOperator
Outputs:
    - An NLPAtX, so all the information of the last iterate
    - A boolean true if we reahced an optimal solution, false otherwise
"""
function Newlbfgs(nlp         :: AbstractNLPModel,
                  nlp_stop    :: NLPStopping;
                  linesearch  :: Function = armijo_ls,
                  verbose     :: Bool=false,
                  mem         :: Int=5,
                  kwargs...)

    nlp_at_x = nlp_stop.current_state
    x = copy(nlp.meta.x0)
    n = nlp.meta.nvar

    xt = copy(nlp.meta.x0)
    ∇f = grad(nlp, nlp.meta.x0)  # à changer pour 0.7

    f = obj(nlp, x)
    # g = grad(nlp, x) # Nécessaire??
    H = InverseLBFGSOperator(n, mem, scaling=true)
    update!(nlp_at_x, x = x, fx = f, gx = ∇f, Hx = H)
    OK = start!(nlp_stop)
    iter = 0
    ∇fNorm = BLAS.nrm2(n, nlp_at_x.gx, 1)

    verbose && @printf("%4s  %8s  %7s  %8s\n", "iter", "f", "‖∇f‖", "∇f'd")
    verbose && @printf("%4d  %8.1e  %7.1e", iter, f, ∇fNorm)

    # d = zeros(nlp_at_x.gx)
    d = zero(nlp_at_x.gx)
    h = LineModel(nlp, x, d)

    while !OK
        d = - H * ∇f
        slope = BLAS.dot(n, d, 1, nlp_at_x.gx, 1)
        if slope > 0.0 # qu'est ce qu'on fait maintenant?
          verbose && @printf("  %8.1e", slope)
        else
          # Perform linesearch.
          h = redirect!(h, xt, d)

          ls_at_t = LSAtT(0.0, h₀ = nlp_at_x.fx, g₀ = slope)
          stop_ls = LS_Stopping(h, (x, y) -> armijo(x, y, τ₀ = 0.01), ls_at_t)
          verbose && println(" ")
          ls_at_t, good_step_size = linesearch(h, stop_ls; kwargs...)

          good_step_size || (nlp_stop.meta.stalled_linesearch = true)

          xt = nlp_at_x.x + ls_at_t.x * d
          ft = obj(nlp, xt); ∇ft = grad(nlp, xt)
          BLAS.blascopy!(n, nlp_at_x.x, 1, xt, 1) # BLAS.blascopy!(n, x, 1, xt, 1)
          BLAS.axpy!(n, ls_at_t.x, d, 1, xt, 1) # BLAS.axpy!(n, t, d, 1, xt, 1)
          ∇ft = grad!(nlp, xt, ∇ft)

          # Update L-BFGS approximation.
          push!(H, ls_at_t.x * d, ∇ft - ∇f)

          # Move on.
          x = xt
          f = ft

          update!(nlp_at_x, x = xt, fx = ft, Hx = H)

          BLAS.blascopy!(n, ∇ft, 1, ∇f, 1)
          # norm(∇f) bug: https://github.com/JuliaLang/julia/issues/11788
          ∇fNorm = BLAS.nrm2(n, ∇f, 1)
          iter = iter + 1

          verbose && @printf("%4d  %8.1e  %7.1e", iter, f, ∇fNorm)
        end
        OK = stop!(nlp_stop)
    end
    verbose && @printf("\n")

    return nlp_stop, nlp_stop.meta.optimal
end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ e9af4550-b632-11eb-1d0d-e11b45468706
html"<button onclick=present()> Presente </button>"

# ╔═╡ 4fb54a13-3db2-4475-bd49-fd0f9d5cefa4
md"""
## Plano Unidade 05 - Matrizes, Tabelas e Dataframes

### Tema relacionado na ementa:

3. Comandos de uma linguagem procedimental: atribuição, entrada e saída, condicionais, repetitivos e seletivos.

- Variáveis escalares e estruturadas homogêneas e heterogêneas. Vetores e encadeamento de estruturas.

### Subitens

**Capítulo 7 - Matriz**

- 7.1 Matriz em algoritmos
- 7.3 Matriz em C/C++
- 7.4 Matriz em JAVA

**Fonte:** ASCENCIO, Ana Fernanda Gomes; DE CAMPOS, Edilene Aparecida Veneruchi. **Fundamentos da programação de computadores**. Pearson Educación, 2008.

**Working with Data Frames**
[Working with Data Frames.](https://dataframes.juliadata.org/stable/man/working_with_dataframes/)

## Resultados de aprendizagem esperados

✔ Discutir estruturas de arrays, tabelas e Dataframes.

✔ Desenvolver estratégias para solução de problemas.

✔ Estabelecer relações interpessoais e trabalhar em equipe, comunicar-se oralmente e ter responsabilidade coletiva.

# Roteiro da Unidade

#### Momento Pré

🕐 1. * Atividade a ser executada por cada grupo (momento aula) **(16/07/2021)**

- _**Apresentar a utilização de arrays mutidimensionais (Matrizes, tabelas e Matrizes) na linguagem julia**_.

🕒 2. (30 min) Apresentação do tema pelo professor **(16/07/2021)**.

🕑 3. (90 min) Ler os subitens 7.1, 7.3 e 7.4 e o material indicado das páginas do [Working with Data Frames.](https://dataframes.juliadata.org/stable/man/working_with_dataframes/) e postar no Classroom resumo de até uma página do que estudou (Resumo pré) **(até 20/07/2021)**;

#### Bibliografia.

* ASCENCIO, Ana Fernanda Gomes; DE CAMPOS, Edilene Aparecida Veneruchi. Fundamentos da programação de computadores. Pearson Educación, 2008.

* Pesquisa a rede computadores. [As fontes devem ser formalmente citadas na norma da ABNT](http://www.biblioteca.fsp.usp.br/~biblioteca/guia/a_modelos.htm).

* [Working with Data Frames.](https://dataframes.juliadata.org/stable/man/working_with_dataframes/)

# Roteiro da aula

#### Momento Aula

🕔 4. (30 min) Seleção dos temas e subtemas com justificativa para que o que será aprendendido e experienciado durante a unidade **(20/07/2021)**.

🕕 5. (60 min) Cada grupo divide os subtemas entre os estudantes do grupo; elaboram os estudos dos subtemas (modelos, algoritmos, códigos, simulações e análises); como um painel de especialistas, apresentam ao grupo os conceitos e questões chaves e apresentações formais dos subtemas em cada grupo (reunião do grupo) **(até 23/07/2021)**.

🕙 6. (10 min/subtema) Apresentações de alguns subtemas pelos membros responsávies em cada grupo (seleção por sorteio a critério do professsor) (sala) **(23/07/2021)**.

## Roteiro

#### Momento Pós

🕚 7. (60 min) Preparação das apresentações de cada um dos grupos e envio do relatório da atividade via classroom (atividade do grupo) **(27/07/2021)**.

🕐 8. (10 min/grupo) Apresentações de cada um dos grupos (sala) **(30/07/2021)**.

🕑 9. (20 min/grupo) Avaliações (sala ou extra sala) **(30/07/2021)**.

🕒 10. (10 min/grupo) Feedback (sala) **(30/07/2021)**.

"""

# ╔═╡ Cell order:
# ╟─e9af4550-b632-11eb-1d0d-e11b45468706
# ╟─4fb54a13-3db2-4475-bd49-fd0f9d5cefa4

max_stars_repo_path: Plano Unidade 05 - Matrizes, Tabelas e Dataframes.jl

max_stars_repo_name: jefferds/Programming-languages

max_issues_repo_name: jefferds/Programming-languages

max_forks_repo_name: jefferds/Programming-languages


text:
export getinstruments, getmiditracks

using MIDI: NoteOnEvent, NoteOffEvent, ProgramChangeEvent, ControlChangeEvent, PitchBendEvent, channelnumber
import MIDI: toabsolutetime!

"""
    Instrument <: Any

Data structure representing a musical instrument.

## Fields:
* `program::Int`: MIDI program of the instrument.
* `notes::Vector{Note}`: Musical notes for the instrument.
* `pitchbends::Vector{PitchBendEvent}`: Pitch bend events for the instrument.
* `controlchanges::Vector{ControlChangeEvent}`: Control change events for the instrument.
"""
Base.@kwdef struct Instrument
    program::Int = 0
    notes::Vector{Note} = []
    pitchbends::Vector{PitchBendEvent} = []
    controlchanges::Vector{ControlChangeEvent} = []
end

function Base.show(io::IO, instrument::Instrument)
    N = length(instrument.notes)
    C = length(instrument.controlchanges)
    P = length(instrument.pitchbends)
    pn = instrument.program
    print(io, "Instrument(program = $pn) with $N Notes, $C Control Changes, $P Pitch Bends")
end

"""
    getinstruments(midi::MIDIFile, time::Symbol=:relative)

Extract [`Instrument`](@ref)s from each [`MIDITrack`](@ref) of a [`MIDIFile`](@ref).
The argument `time` informs whether the `MIDIFile` is in `:absolute` or `:relative` time.
Internally, `getinstruments` converts relative time `MIDIFile`s to absolute time
for extracting the instruments.

If there are notes in a track that do not have a [`MIDI.ProgramChangeEvent`](@ref)
preceding them, an `Instrument` with a default program number of zero is created.

See also: [`MIDI.getnotes`](@ref)
"""
function getinstruments(midi::MIDIFile, time::Symbol=:relative)
    midi = deepcopy(midi)

    if time === :relative
        toabsolutetime!(midi)
    elseif time != :absolute
        throw(ArgumentError("`time` should be either :absolute or :relative"))
    end

    # Create an instrument map which maps:
    # (program, channel, track) => Instrument
    instrument_map = OrderedDict{Tuple{Int64, Int64, Int64}, Instrument}()

    for (track_num, track) in enumerate(midi.tracks)
        # Create a map of channel => program to map the channels to the current playing instrument
        # Initially all channels map to the default program (0)
        current_instrument = zeros(Int, 16)

        # Create a map of (channel, note) => (NoteOn position, velocity)
        note_map = Dict{Tuple{Int, Int}, Vector{Tuple{Int, Int}}}() # TODO: Use default dict instead

        for event in track.events
            if event isa NoteOnEvent && event.velocity > 0
                key = (channel(event), event.note)
                if haskey(note_map, key)
                    push!(note_map[key], (event.dT, event.velocity))
                else
                    note_map[key] = [(event.dT, event.velocity)]
                end
            elseif event isa NoteOffEvent || (event isa NoteOnEvent && event.velocity == 0)
                k = (channel(event), event.note)
                if haskey(note_map, k)
                    end_pos = event.dT
                    notes_playing = note_map[k]
                    turnoff_notes = [(pos, vel) for (pos, vel) in notes_playing if pos != end_pos]
                    continue_notes = [(pos, vel) for (pos, vel) in notes_playing if pos == end_pos]

                    program = current_instrument[channel(event) + 1]
                    key = (program, channel(event), track_num)
                    if !haskey(instrument_map, key)
                        # Create an instrument with the current program number
                        instrument = Instrument(program=program)
                        instrument_map[key] = instrument
                    else
                        instrument = instrument_map[key]
                    end

                    for (position, velocity) in turnoff_notes
                        note = Note(event.note, velocity, position, end_pos - position, channel(event))
                        push!(instrument.notes, note)
                    end

                    if length(turnoff_notes) > 0 && length(continue_notes) > 0
                        note_map[k] = continue_notes
                    else
                        delete!(note_map, k)
                    end
                end
            elseif event isa ProgramChangeEvent
                current_instrument[channel(event) + 1] = event.program
                # If the intrument map has an instrument at the default program with no notes,
                # Change it's program to the current program
                key = (0, channel(event), track_num)
                if haskey(instrument_map, key) && event.program != 0
                    instrument = instrument_map[key]
                    if isempty(instrument.notes)
                        newinstrument = Instrument(event.program, instrument.notes, instrument.pitchbends, instrument.controlchanges)
                        instrument_map[(event.program, channel(event), track_num)] = newinstrument
                        delete!(instrument_map, key)
                    end
                end
            elseif typeof(event) in [ControlChangeEvent, PitchBendEvent]
                program = current_instrument[channel(event) + 1]
                key = (program, channel(event), track_num)
                if !haskey(instrument_map, key)
                    instrument = Instrument(program=program)
                    instrument_map[key] = instrument
                else
                    instrument = instrument_map[key]
                end

                # Store a copy of the events
                if event isa ControlChangeEvent
                    push!(instrument.controlchanges, ControlChangeEvent(event.dT, event.status, event.controller, event.value))
                elseif event isa PitchBendEvent
                    push!(instrument.pitchbends, PitchBendEvent(event.dT, event.status, event.pitch))
                end
            end
        end
    end

    instruments = Vector{Instrument}()
    for instrument in values(instrument_map)
        push!(instruments, instrument)
    end
    return instruments
end

"""
    getmiditracks(instruments::Vector{Instrument})

Return [`MIDITrack`](@ref)s from the [`Instrument`](@ref)s.
"""
function getmiditracks(instruments::Vector{Instrument})
    tracks = MIDITrack[]
    for (ins_num, ins) in enumerate(instruments)
        track = MIDITrack()
        # Get the channel number for the current instrument
        ins_channel = (ins_num - 1) % 16

        # Add a program change event at the beginning
        push!(track.events, ProgramChangeEvent(0, ins.program, channel = ins_channel))

        # Add control change events
        for event in ins.controlchanges
            push!(track.events, ControlChangeEvent(event.dT, event.controller, event.value, channel = ins_channel))
        end

        # Add pitch bend events
        for event in ins.pitchbends
            push!(track.events, PitchBendEvent(event.dT, event.pitch, channel = ins_channel))
        end

        # Sort the events by time and convert them to relative time
        sort!(track.events, by=event->event.dT)
        torelativetime!(track)

        # Alter the channel for each note
        for note in ins.notes
            note.channel = ins_channel
        end
        addnotes!(track, Notes(ins.notes))

        push!(tracks, track)
    end

    tracks
end

# TODO: Send a PR with the below 3 functions to MIDI.jl

"""
    toabsolutetime!(midi::MIDIFile)

Convert a midi file from relative time to absolute time.

See also: [`MIDI.toabsolutetime!`](@ref), [`torelativetime!`](@ref)
"""
function toabsolutetime!(midi::MIDIFile)
    for track in midi.tracks
        toabsolutetime!(track)
    end
    midi
end

"""
    torelativetime!(midi::MIDITrack)

Convert a midi track from absolute time to relative time.

See also: [`MIDI.toabsolutetime!`](@ref)
"""
function torelativetime!(track::MIDITrack)
    time = 0
    for event in track.events
        event.dT -= time
        time += event.dT
    end
end

"""
    torelativetime!(midi::MIDIFile)

Convert a midi file from absolute time to relative time.

See also: [`MIDI.toabsolutetime!`](@ref)
"""
function torelativetime!(midi::MIDIFile)
    for track in midi.tracks
        torelativetime!(track)
    end
    midi
end

channel(event::MIDIEvent) = Int(channelnumber(event))


text:
module GmmHmcTest
import Pkg; Pkg.activate("../../")  # Cytof5

include("HMC.jl")
using Distributions, Random
using Flux, Flux.Tracker
using Flux.Tracker: @grad, track, lgamma
include("../../src/VB/ADVI/ADVI.jl")
include("../../src/VB/ADVI/custom_grads.jl")
include("../../src/VB/ADVI/custom_functions.jl")
include("../../src/VB/ADVI/StickBreak.jl")

using RCall
@rimport graphics as rgraphics
@rimport rcommon

Random.seed!(0)

struct State
  mu
  log_sig
  stickbreak_w
end

pretty(s::State) = ("State(μ => $(s.mu), " *
                    "σ => $(exp.(s.log_sig)), " *
                    "w => $(SB_transform(s.stickbreak_w)))")

State(K::Int) = State(param(randn(K)),  # mu
                      param(randn(K)),  # sig
                      param(randn(K - 1)))  # stickbreak_w

function sim_data(; N=100, mu=[-2, 3], sig=[.2, .4], w=[.2, .8], seed=missing)
  if !ismissing(seed)
    Random.seed!(seed)
  end

  K = length(mu)
  @assert K == length(sig) == length(w)
  @assert sum(w) == 1

  cluster_labels = wsample(1:K, w, N)
  y = randn(N) .* sig[cluster_labels] + mu[cluster_labels]

  return Dict(:cluster_labels => cluster_labels,
              :mu => mu, :sig => sig, :w => w, :y => y)
end

### MAIN ###

# Simulate data
simdat = sim_data(N=300)
N = length(simdat[:y])
K = 2

# Log posterior
function logpost(s::State)
  K = length(s.mu)
  mu = reshape(s.mu, 1, K)
  sig = reshape(exp.(s.log_sig), 1, K)
  # w = SB_transform(s.stickbreak_w)  # FIXME: I expected this would work.
  real_w = reshape(s.stickbreak_w, 1, K - 1)
  w = SB_transform(real_w)

  ll = (ADVI.lpdf_normal.(reshape(simdat[:y], N, 1), mu, sig) .+
        log.(w))
  ll = sum(ADVI.logsumexp(ll, dims=2))

  lp = (sum(ADVI.compute_lpdf(Normal(0, 3), mu)) +
        sum(ADVI.compute_lpdf(LogNormal(0, 1), sig)) +
        sum(ADVI.compute_lpdf(Dirichlet(K, K), w)))

  lpabsjacobian = (sum(s.log_sig) +
                   sum(SB_logabsdetJ(real_w, w)))

  return ll + lp + lpabsjacobian
end

state = State(K)
_ = logpost(state)

function simulate(init; nburn, nsamps, eps, num_leapfrog_steps,
                  kappa=nothing, momentum_sd=1.0)
  state = deepcopy(init)

  samps = [state for i in 1:nsamps]
  log_prob_hist = zeros(nsamps)

  for i in 1:(nburn + nsamps)
    print("\rProgress: $i / $(nburn + nsamps)")
    eta = (kappa == nothing ? eps : eps * i ^ -kappa)
    state, curr_log_prob = HMC.hmc_update(state, logpost,
                                          num_leapfrog_steps, eta,
                                          momentum_sd=momentum_sd)
    if i > nburn
      samps[i - nburn] = state
      log_prob_hist[i - nburn] = curr_log_prob
    end
  end
  println()

  return Dict(:samps => samps, :log_prob => log_prob_hist)
end

# Compile
_ = simulate(state, nburn=1, nsamps=1, num_leapfrog_steps=1, eps=.1, kappa=.7)

# Simulate
@time out = simulate(state, nburn=500, nsamps=300,
                     num_leapfrog_steps=2^5, eps=1/N)

# FIXME: runs. but can't get correct answer...
println(simdat[:mu], simdat[:sig], simdat[:w])
println()

println("mu post mean: $(mean([o.mu for o in out[:samps]]))")
println("sig post mean: $(mean([exp.(o.log_sig) for o in out[:samps]]))")
println("w post mean: $(mean([SB_transform(o.stickbreak_w) for o in out[:samps]]))")

R"plot($(out[:log_prob]), type='l')"
end # module GmmHmcTest


text:
s = open("input.txt") do file
    read(file, String);
end

function ParseInput(s)
    startStateFlat = map(x->x=='#', replace(s,"\n"=>"") |> collect);
    dim = length(split(s,"\n")[1])
    startState = reshape(startStateFlat, (dim,dim));

    simSpace = falses(21, 21, 21, 21);

    halfDim = dim ÷ 2;

    simSpace[11-halfDim:11+halfDim - ifelse(dim%2==0,1,0),
        11-halfDim:11+halfDim - ifelse(dim%2==0,1,0),
        11, 11] .= startState;

    return simSpace;
end

function numNeighbors(space, loc)
    count = 0;
    for i = -1:1
        for j = -1:1
            for k = -1:1
                for l = -1:1
                    tempLoc = loc .+ [i,j,k,l];
                    if all(tempLoc .== loc)
                        continue;
                    end

                    if all(tempLoc.>0) && all(tempLoc.<22)
                        if space[tempLoc...]
                            count += 1;
                        end
                    end
                end
            end
        end
    end

    return count;
end

function BootUp(simSpace)
    nextSpace = falses(21,21,21,21);
    for step = 1:6
        for k = 1:21
            for i = 1:21
                for j = 1:21
                    for l = 1:21
                        curLoc = [i, j, k, l];
                        nCount = numNeighbors(simSpace, curLoc);

                        active = simSpace[curLoc...];
                        if active
                            # println(nCount);
                            if 1 < nCount < 4
                                nextSpace[curLoc...] = true;
                            else
                                nextSpace[curLoc...] = false;
                            end
                        else
                            if nCount == 3
                                nextSpace[curLoc...] = true;
                            else
                                nextSpace[curLoc...] = false;
                            end
                        end
                    end
                end
            end

            # println(k);
            # println(simSpace[:,:,k])
        end
        println("--------");
        simSpace .= nextSpace;
    end

    totalActive = sum(simSpace);
    println(totalActive);

    return totalActive;
end

simSpace = ParseInput(s);
# simSpace = ParseInput(".#.\n..#\n###");
answer = BootUp(simSpace)


text:
# Julia wrapper for header: include/blockmat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# TODO detect size and use bitstype because if the DLL changes and gets
#      compiled with NOSHORTS we are screwed with the following code...
@static if is_windows()
    typealias csdpshort Cushort
  else
    typealias csdpshort Cint
end

# begin enum blockcat
typealias blockcat Cuint
const DIAG = (blockcat)(0)
const MATRIX = (blockcat)(1)
const PACKEDMATRIX = (blockcat)(2)
# end enum blockcat

immutable blockdatarec
    _blockdatarec::Ptr{Cdouble}
end

immutable blockrec
    data::blockdatarec
    blockcategory::blockcat
    blocksize::csdpshort
end

type blockmatrix
    nblocks::Cint
    blocks::Ptr{blockrec}
end

type sparseblock
    next::Ptr{sparseblock}
    nextbyblock::Ptr{sparseblock}
    entries::Ptr{Cdouble}
    iindices::Ptr{csdpshort}
    jindices::Ptr{csdpshort}
    numentries::Cint
    blocknum::csdpshort
    blocksize::csdpshort
    constraintnum::csdpshort
    issparse::csdpshort
end

immutable constraintmatrix
    blocks::Ptr{sparseblock}
end

# Skipping MacroDefinition: ijtok ( iiii , jjjj , lda ) ( ( jjjj - 1 ) * lda + iiii - 1 )
# Skipping MacroDefinition: ijtokp ( iii , jjj , lda ) ( ( iii + jjj * ( jjj - 1 ) / 2 ) - 1 )
# Skipping MacroDefinition: ktoi ( k , lda ) ( ( k % lda ) + 1 )
# Skipping MacroDefinition: ktoj ( k , lda ) ( ( k / lda ) + 1 )

type paramstruc
    axtol::Cdouble
    atytol::Cdouble
    objtol::Cdouble
    pinftol::Cdouble
    dinftol::Cdouble
    maxiter::Cint
    minstepfrac::Cdouble
    maxstepfrac::Cdouble
    minstepp::Cdouble
    minstepd::Cdouble
    usexzgap::Cint
    tweakgap::Cint
    affine::Cint
    perturbobj::Cdouble
    fastmode::Cint
end


text:
using ITensors
using ITensorMakie
using GLMakie
using Graphs: Graph
using Random

using ITensors.ContractionSequenceOptimization: optimal_contraction_sequence

Random.seed!(1234)

N = 5
g = Graph(N, N)
A = itensornetwork(g; linkspaces=5)
sequence = optimal_contraction_sequence(A)
edge_labels = (; tags=true)
R = @visualize_sequence fig ITensors.contract(A; sequence=sequence) edge_labels =
  edge_labels

fig


text:
@inline Base.map!(f::F, A::AbstractStrideArray, arg1::AbstractArray, args::Vararg{AbstractArray,K}) where {F,K} = vmap!(f, A, arg1, args...)
# these two definitions are to avoid ambiguities
@inline Base.map!(f::F, A::AbstractStrideArray, arg::AbstractArray) where {F} = vmap!(f, A, arg)
@inline Base.map!(f::F, A::AbstractStrideArray, arg1::AbstractArray, arg2::AbstractArray) where {F} = vmap!(f, A, arg1, arg2)
@inline Base.map(f::F, A::AbstractStrideArray, args::Vararg{Any,K}) where {F,K} = vmap(f, A, args...)
@inline Base.reduce(op::O, A::AbstractStrideArray) where {O} = vreduce(op, A)
@inline Base.mapreduce(f::F, op::O, A::AbstractStrideArray, args::Vararg{AbstractArray,K}) where {F, O, K} = vmapreduce(f, op, A, args...)
@inline function Base.mapreduce(f::F, op::O, A::AbstractStrideArray) where {F, O}
    @gc_preserve vmapreduce(f, op, A)
end

for (op,r) ∈ ((:+,:sum), (:max,:maximum), (:min,:minimum))
    @eval begin
        @inline Base.reduce(::typeof($op), A::AbstractStrideArray; dims = nothing) = @gc_preserve vreduce($op, A, dims = dims)
        @inline Base.$r(A::AbstractStrideArray; dims = nothing) = @gc_preserve vreduce($op, A, dims = dims)
    end
end

function maximum(::typeof(abs), A::AbstractStrideArray{S,T}) where {S,T}
    s = typemin(T)
    @avx for i ∈ eachindex(A)
        s = max(s, abs(A[i]))
    end
    s
end

function Base.vcat(A::AbstractStrideMatrix, B::AbstractStrideMatrix)
    MA, NA = size(A)
    MB, NB = size(B)
    @assert NA == NB
    TC = promote_type(eltype(A), eltype(B))
    C = StrideArray{TC}(undef, (MA + MB, NA))
    # TODO: Actually handle offsets
    @assert offsets(A) === offsets(B)
    @assert offsets(A) === offsets(C)
    @avx for j ∈ axes(A,2), i ∈ axes(A,1)
        C[i,j] = A[i,j]
    end
    @avx for j ∈ axes(B,2), i ∈ axes(B,1)
        C[i + MA,j] = B[i,j]
    end
    C
end


text:
using ReachabilityBenchmarks, MathematicalSystems, Reachability, Plots

include(@relpath "projectile_model.jl")
include(@relpath "projectile_specifications.jl")

S = projectile_model()
X0, options = projectile_specification()
options = Options(options)

# initial value problem
problem = InitialValueProblem(S, X0)

# partition with 2D blocks
partition = [2*i-1:2*i for i in 1:2]

# reachability algorithms
algorithm_reach = BFFPSV18(:δ => 0.5, :vars => [1, 3], :partition => partition)

# compute flowpipe
options[:mode] = "reach"
solution = solve(problem, options; op=algorithm_reach)

# project flowpipe to x₁ and x₃
solution.options[:plot_vars] = [1, 3]
solution_proj = project(solution)

# plot projection
plot(solution_proj)
savefig("projectile")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
function SelfEnergy(domains::NTuple{L,Domain{1}},surface,domain_index,α_half,a,nnm,nnp,indices,interior,domain_wall) where L
    li = findfirst(surface)
    ri = findlast(surface)
    # surface[1] || throw("first site should be an endpoint, something is amiss")
    # surface[end] || throw("last site should be an endpoint, something is amiss")
    sum(surface)>2 && throw("no point except first and last should be an endpoint, something is amiss")

    N = length(surface)
    α_half⁻¹ = 1 ./α_half

    dx1 = domains[domain_index[li]].lattice.dx
    bc1 = domains[domain_index[li]].boundary.bcs[1]

    if typeof(bc1)<:DirichletBC{1}
        Σ0L = sparse([li],[li],[-α_half⁻¹[li]/dx1^2],N,N)
        Σ1L = sparse([li],[li],[+1/dx1],N,N)
        fL = FM(dx1,0)
    elseif typeof(bc1)<:NeumannBC{1}
        Σ0L = sparse([li],[li],[+α_half⁻¹[li]/dx1^2],N,N)
        Σ1L = sparse([li],[li],[-1/dx1],N,N)
        fL = FM(dx1,0)
    elseif typeof(bc1)<:FloquetBC{1}
        # Σ0L = sparse([li],[ri],[α_half⁻¹[ri+1]/dx1^2],N,N)
        # Σ1L = sparse([li],[ri],[-1/dx1],N,N)
        # fL = (ω,ka,_...) -> ffl(ka,a)
        throw("not yet")
    elseif typeof(bc1)<:MatchedBC{1}
        Σ0L = sparse([li],[li],[+α_half⁻¹[li]/dx1^2],N,N)
        Σ1L = sparse([li],[li],[-1/dx1],N,N)
        if isempty(bc1.in)
            fL = FM(dx1,+1)
        else
            fL = FM(dx1,-1)
        end
    else
        Σ0L = spzeros(N,N)
        Σ1L = spzeros(N,N)
        fL = FM(dx1,0)
    end

    dx2 = domains[domain_index[ri]].lattice.dx
    bc2 = domains[domain_index[ri]].boundary.bcs[2]
    if typeof(bc2)<:DirichletBC{2}
        Σ0R = sparse([ri],[ri],[-α_half⁻¹[ri+1]/dx2^2],N,N)
        Σ1R = sparse([ri],[ri],[-1/dx2],N,N)
        fR = FM(dx2,0)
    elseif typeof(bc2)<:NeumannBC{2}
        Σ0R = sparse([ri],[ri],[+α_half⁻¹[ri+1]/dx2^2],N,N)
        Σ1R = sparse([ri],[ri],[+1/dx2],N,N)
        fR = FM(dx2,0)
    elseif typeof(bc2)<:FloquetBC{2}
        # Σ0R = sparse([ri],[1],[-α_half⁻¹[li]/dx2^2],N,N)
        # Σ1R = sparse([ri],[1],[+1/dx2],N,N)
        # fR = (ω,ka,_...) -> ffr(ka,a)
        throw("not yet")
    elseif typeof(bc2)<:MatchedBC{2}
        Σ0R = sparse([ri],[ri],[+α_half⁻¹[ri+1]/dx2^2],N,N)
        Σ1R = sparse([ri],[ri],[+1/dx2],N,N)
        if isempty(bc2.in)
            fR = FM(dx2,+1)
        else
            fR = FM(dx2,-1)
        end
    else
        Σ0R = spzeros(N,N)
        Σ1R = spzeros(N,N)
        fR = FM(dx1,0)
    end
    Σ2 = spzeros(N,N)

    fs = (fL,fR)

    return SelfEnergy{1,2,typeof(fs)}((-Σ0L,-Σ0R),(Σ1L,Σ1R),Σ2,(fL,fR))
end

# fdl(ω,_...) = one(ω)
# fnl(ω,_...) = one(ω)
# ffl(ka,a) = exp(-1im*ka*a)

# fdr(ω,_...) = one(ω)
# fnr(ω,_...) = one(ω)
# ffr(ka,a) = exp(+1im*ka*a)

# fmin(ω,ka,ky,kz,dx)  = exp(-1im*sqrt(ω^2-one(ω)*ky^2-one(ω)*kz^2)*dx)
# fmout(ω,ka,ky,kz,dx) = exp(+1im*sqrt(ω^2-one(ω)*ky^2-one(ω)*kz^2)*dx)

struct FM
    dx::Float64
    sign::Int
end
function (fm::FM)(ω,ka,ky,kz)
    if iszero(fm.sign)
        complex(1,0)
    else
        exp(fm.sign*1im*sqrt(ω^2-one(ω)*ky^2-one(ω)*kz^2)*fm.dx)
    end
end


text:
type SequencePerLine; file; dict; unk;
    function SequencePerLine(file; dict=nothing, o...)
        if dict == nothing
            unk = 0             # we will construct dict from data
            dict = Dict{Any,Int}()
        else
            isa(dict,Dict) || (dict = readvocab(dict))
            unk = length(dict)+1
        end
        new(file, dict, unk)
    end
end

unk(s::SequencePerLine)=s.unk
maxtoken(s::SequencePerLine)=(s.unk > 0 ? s.unk : error("Please specify a dict for maxtoken"))

start(s::SequencePerLine)=open(s.file)
done(s::SequencePerLine,io)=(isa(io,Tuple) && (io=io[1]); eof(io) && (close(io); true))

function next(s::SequencePerLine, io)
    line = readline(isa(io,Tuple) ? io[1] : io)
    sent = Int32[]
    if s.unk == 0
        for token in split(line)
            push!(sent, get!(s.dict, token, 1+length(s.dict)))
        end
    else
        for token in split(line)
            push!(sent, get(s.dict, token, s.unk))
        end
    end
    return (sent, io)
end

### DEAD CODE

# elseif OLD_S2S==1
#     n = length(dict)
#     unk = get(dict, "<unk>", n+1)
#     eos = get(dict, "<s>", unk==n+1 ? n+2 : n+1)
#     maxtoken = length(dict)
# elseif OLD_S2S==2
#     unk = get(dict, "<unk>", length(dict)+1)
#     eos = length(dict)+1
#     maxtoken = length(dict)+1


text:
using HJBSolver, FactCheck

type TestProblem2D{T<:Real}
    gamma::Vector{T}
    c::Vector{T}
    model::HJBTwoDim{T}
    truevaluefun::Function
    truecontrolfun::Function
end

function createmodel(;T::Float64=1.)
    γ = [5e-2, 0.1]
    c = [-1., -0.9]

    xmin = zeros(2); xmax = [0.5,0.5]

    truesol(t,x) = (t-T)*(x[1]+x[2])^2+vecdot(c,x)
    truepol(t,x) = (1+c+2*(t-T)*(x[1]+x[2]-γ.^2))./(2-2*(t-T)*γ.^2)
    hatf(t,x) = (x[1]+x[2])^2+sum((1-c-2*(t-T)*(x[1]+x[2])).^2./(4-4*(t-T)*γ.^2))

    amin = truepol(0., xmax) - 1e-1
    amax = truepol(0., xmin) + 1e-1

    b(t,x,a) = a-1.
    σ(t,x,a) = (γ.*b(t,x,a)).^2
    f(t,x,a) = -vecdot(a,b(t,x,a)) - hatf(t,x)
    g(x) = truesol(T,x)

    model = HJBTwoDim(b, σ, f, g, T, amin, amax, xmin, xmax, truesol)
    return TestProblem2D(γ, c, model, truesol, truepol)
end

function calculateerror(prob::TestProblem2D, K::Vector{Int}, v,pol)
    model = prob.model

    x1 = linspace(model.xmin[1], model.xmax[1], K[1])
    x2 = linspace(model.xmin[2], model.xmax[2], K[2])
    x = (collect(x1), collect(x2))

    i,j = ceil(Int, K/2)
    idxi = K[2]*(i-1) + j
    xij = [x1[i], x2[j]]
    w = prob.truevaluefun(0., xij)
    α = prob.truecontrolfun(0., xij)

    return abs(w-v[idxi,1])/abs(w), abs(α[1]-pol[1][idxi,1])/abs(α[1]),
    abs(α[2]-pol[2][idxi,1])/abs(α[2])
end

facts("2D, Policy iteration") do
    K = [26, 26]; N = 10
    Δt = 5e-2
    prob = createmodel(T=Δt*N)
    model = prob.model

    @time v, pol = solve(model, K, N)
    errv, erra1, erra2 = calculateerror(prob, K, v, pol)

    @fact errv --> roughly(0., 5e-3)
    @fact erra1 --> roughly(0., 5e-2)
    @fact erra2 --> roughly(0., 5e-2)
end

facts("2D, Policy timestepping") do
    K = [21, 21]; N = 5; M = (25,25)
    Δt = 1e-1
    prob = createmodel(T=Δt*N)
    model = prob.model

    avals1 = linspace(model.amin[1], model.amax[1], M[1])
    avals2 = linspace(model.amin[2], model.amax[2], M[2])

    @time v, pol = solve(model, K, N, (avals1,avals2))

    errv, erra1, erra2 = calculateerror(prob, K, v, pol)
    @fact errv --> roughly(0., 1e-2)
    @fact erra1 --> roughly(0., 0.3)
    @fact erra2 --> roughly(0., 0.3)
end

#==

function calctruesols(t, x, prob, flat=false)
    K = [length(xi) for xi in x]
    v = zeros(reverse(K)...)
    pol1 = zeros(v)
    pol2 = zeros(v)

    for i = 1:K[1], j = 1:K[2]
        xij = [x[1][i],x[2][j]]
        v[j,i] = prob.truevaluefun(t, xij)
        pol1[j,i], pol2[j,i] = prob.truecontrolfun(t, xij)
    end
    if flat == true
        v = v[:]
        pol1 = pol1[:]
        pol2 = pol2[:]
    end
    return v, (pol1, pol2)
end

K = [51, 51]; N = 50
Δt = 2e-2
prob = createmodel(T=Δt*N)
model = prob.model
x1 = linspace(model.xmin[1], model.xmax[1], K[1])
x2 = linspace(model.xmin[2], model.xmax[2], K[2])
x = (collect(x1), collect(x2))
v, pol = solve(prob.model, K, N)

tv = zeros(v)
tpol1 = zeros(pol[1])
tpol2 = zeros(pol[2])
for i = 1:N
tv[:,i], tmppol = calctruesols((i-1)*Δt, x, prob, true)
tpol1[:,i] = tmppol[1]
tpol2[:,i] = tmppol[2]
end
tv[:,end], tmppol = calctruesols(model.T,x,prob, true)
==#

max_stars_repo_name: anriseth/HJBSolver

max_issues_repo_name: anriseth/HJBSolver

max_forks_repo_name: anriseth/HJBSolver


text:
using DataDeps

register(DataDep(
    "lymphography",
    """
    Dataset: lymphography
    Author: Igor Kononenko, Bojan Cestnik
    Website: https://archive.ics.uci.edu/ml/datasets/Lymphography

    Observations: 148
    Features:     18
    Classes:      4

    This is one of three domains provided by the Oncology Institute that has repeatedly appeared in the machine learning literature. (See also breast-cancer and primary-tumor.)
    """,
    "https://archive.ics.uci.edu/ml/machine-learning-databases/lymphography/lymphography.data",
    "e17d22d071341af08a505a10eb558afc835e302b02c8b381c1fa3dd305731935",
    post_fetch_method = (path -> begin
        MVDA.process_dataset(path,
            header=false,
            target_index=1,
            feature_indices=2:19,
            ext=".csv",
        )
    end),
))


text:
module TestMarkovApprox

using QuantEcon
using Base.Test
using FactCheck

rough_kwargs = Dict(:atol => 1e-8, :rtol => 1e-8)

# set up
srand(42)
ρ, σ_u = rand(2)
μ = 0.0

facts("Testing markov_approx.jl") do
    for n=3:25
        x, P = rouwenhorst(n, ρ, σ_u, μ)

        @fact size(P, 1) --> size(P, 2)
        @fact ndims(x) --> 1
        @fact ndims(P) --> 2
        @fact sum(P, 2) --> roughly(ones(n, 1); rough_kwargs...)
        @fact all(P .>= 0.0) --> true
        @fact sum(x) --> roughly(0.0; rough_kwargs...)

        for m=1:4
            x, P = tauchen(n, ρ, σ_u, μ, m)

            @fact size(P, 1) --> size(P, 2)
            @fact ndims(x) --> 1
            @fact ndims(P) --> 2
            @fact sum(P, 2) --> roughly(ones(n, 1); rough_kwargs...)
            @fact all(P .>= 0.0) --> true
            @fact sum(x) --> roughly(0.0; rough_kwargs...)
        end
    end
end  # facts
end  # module


text:
function BilinearInterpolatingFunction(arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{Vector{jdouble}})
    return BilinearInterpolatingFunction((Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function get_x_inf(obj::BilinearInterpolatingFunction)
    return jcall(obj, "getXInf", jdouble, ())
end

function get_x_sup(obj::BilinearInterpolatingFunction)
    return jcall(obj, "getXSup", jdouble, ())
end

function get_y_inf(obj::BilinearInterpolatingFunction)
    return jcall(obj, "getYInf", jdouble, ())
end

function get_y_sup(obj::BilinearInterpolatingFunction)
    return jcall(obj, "getYSup", jdouble, ())
end

function to_real_field_bivariate_function(obj::FieldBivariateFunction, arg0::Field)
    return jcall(obj, "toRealFieldBivariateFunction", RealFieldBivariateFunction, (Field,), arg0)
end

function value(obj::BilinearInterpolatingFunction, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(obj, "value", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function value(obj::BilinearInterpolatingFunction, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "value", jdouble, (jdouble, jdouble), arg0, arg1)
end


text:
# See @ianshmean for name

"""
Returns a list of dependencies required by a package that aren't actually used.
"""
function trim_dependencies(s)
    named_deps = get_dependencies(s)
    test_path = joinpath(Pkg.dir(s), "test", "runtests.jl")
    @info "Snooping $test_path"

    # snoop the tests
    SnoopCompile.@snoopc "/tmp/$s.log" begin
        using ColorTypes, Pkg
        include(joinpath(dirname(dirname(pathof(ColorTypes))), "test", "runtests.jl"))
    end

    # parse the test coverage
    snoopfile = read("/tmp/$s.log")
    used_deps = filter(dep -> occursin(dep, snoopfile), named_deps)
    return setdiff(named_deps, used_deps)
end


text:
using OrdinaryDiffEq
using DiffEqUncertainty
using DiffEqCallbacks
using Statistics
using Random
using Plots
using BenchmarkTools

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[2]
end;

function affect!(integrator)
  integrator.u[2] = 0.0
end;

positive_cb = ContinuousCallback(condition,affect!);

tmax = 40.0
tspan = (0.0,tmax);

u0 = [990.0,10.0,0.0]; # S,I,R

p = [0.05,10.0,0.25]; # β,c,γ

Random.seed!(1234);

prob_ode = ODEProblem(sir_ode!, u0, tspan, p);

probint_cb_4th = AdaptiveProbIntsUncertainty(4);

num_samples = 100
ensemble_prob_ode = EnsembleProblem(prob_ode)

samples_ode = solve(ensemble_prob_ode,
                             ROS34PW3(),
                             trajectories=num_samples,
                             callback=CallbackSet(positive_cb,probint_cb_4th));

s20 = [s(20.0) for s in samples_ode]
[[mean([s[i] for s in s20]) for i in 1:3] [std([s[i] for s in s20]) for i in 1:3]]

p = plot(samples_ode[1],
     label=["S" "I" "R"],
     color=[:blue :red :green],
     xlabel="Time",
     ylabel="Number")
for i in 2:num_samples
    plot!(p,
          samples_ode[i],
          label="",
          color=[:blue :red :green])
end;

plot(p,yaxis=:log10,xlim=(15,20),ylim=(100,1000))

plot(p,yaxis=:log10,xlim=(35,40),ylim=(10,1000))

@benchmark solve(ensemble_prob_ode,
                 ROS34PW3(),
                 trajectories=100,
                 callback=CallbackSet(positive_cb,probint_cb_4th))


text:
# Script to create the first three panels of Figure 3

using Plots
include("./optimizers_abstract.jl")
#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(4.5)
using Plots.PlotMeasures

γ = 1.00
η = 0.20
ε = 1e-6

nIter = 100
xHistGDA = zeros(1, nIter + 1)
yHistGDA = zeros(1, nIter + 1)
xHistLCGD = zeros(1, nIter + 1)
yHistLCGD = zeros(1, nIter + 1)
xHistConOpt = zeros(1, nIter + 1)
yHistConOpt = zeros(1, nIter + 1)
xHistOGDA = zeros(1, nIter + 1)
yHistOGDA = zeros(1, nIter + 1)
xHistCGDA = zeros(1, nIter + 1)
yHistCGDA = zeros(1, nIter + 1)

α = 6
f(x,y) = -α * (x'*x - y'*y)
g(x,y) = + α * (x'*x - y'*y)

x0 = [0.5]
y0 = [0.5]

xHistGDA[:,1] .= x0
yHistGDA[:,1] .= y0
xHistLCGD[:,1] .= x0
yHistLCGD[:,1] .= y0
xHistConOpt[:,1] .= x0
yHistConOpt[:,1] .= y0
xHistOGDA[:,1] .= x0
yHistOGDA[:,1] .= y0
xHistCGDA[:,1] .= x0
yHistCGDA[:,1] .= y0

∇x_f, ∇y_f, ∇x∇x_f, ∇x∇y_f, ∇y∇y_f, ∇y∇x_f = get_closures(f)
∇x_g, ∇y_g, ∇x∇x_g, ∇x∇y_g, ∇y∇y_g, ∇y∇x_g = get_closures(g)

ΔxC = zeros(1)
ΔyC = zeros(1)

ΔxO = zeros(1)
ΔyO = zeros(1)

for k = 1 : nIter
  global ΔxC
  global ΔyC
  global ΔxO
  global ΔyO

  x, y = xHistGDA[:, k], yHistGDA[:, k]
  GDA!(x,
       y,
       ∇x_f, ∇y_g, η)
  xHistGDA[:, k + 1], yHistGDA[:, k + 1] = x,y

    x, y = xHistOGDA[:, k], yHistOGDA[:, k]
  OGDA!(x,
        y,
        ΔxO, ΔyO, ∇x_f, ∇y_g, η)
  xHistOGDA[:, k + 1], yHistOGDA[:, k + 1] = x, y

  x, y = xHistCGDA[:, k], yHistCGDA[:, k]
  CG_CGDA!(x,
           y,
           ΔxC, ΔyC , ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, η, ε)
  xHistCGDA[:, k + 1], yHistCGDA[:, k + 1] = x, y

  x, y = xHistLCGD[:, k], yHistLCGD[:, k]
  LCGDA!(x,
         y,
         ∇x_f, ∇x∇y_f, ∇y_g, ∇y∇x_g, η,)
  xHistLCGD[:, k + 1], yHistLCGD[:, k + 1] = x, y

  x, y = xHistConOpt[:, k], yHistConOpt[:, k]
  conOpt!(x,
          y,
          ∇x_f, ∇y_g, ∇x∇x_f, ∇x∇y_f, ∇y∇y_g, ∇y∇x_g, η, γ)
  xHistConOpt[:, k + 1], yHistConOpt[:, k + 1] = x, y
end

l10DistGDA = log10.(sqrt.(xHistGDA.^2 + yHistGDA.^2))[:]./2
l10DistOGDA = log10.(sqrt.(xHistOGDA.^2 + yHistOGDA.^2))[:]./2
l10DistCGDA = log10.(sqrt.(xHistCGDA.^2 + yHistCGDA.^2))[:]./2
l10DistLCGD = log10.(sqrt.(xHistLCGD.^2 + yHistLCGD.^2))[:]./2
l10DistConOpt = log10.(sqrt.(xHistConOpt.^2 + yHistConOpt.^2))[:]./2

pl = plot(# xlabel="Iterations",
          # ylabel="\$\\log_{10}(r)\$",
          linewidth = 5.0,
          legend=:false,
          # bottom_margin= 10 * mm,
          # left_margin= 7 * mm,
          right_margin= 5 * mm,
          #top_margin= 15 * mm
          )

plot!(pl, l10DistGDA, label="GDA, SGA, LCGD, CGD", linewidth = 5.0)
plot!(pl, l10DistConOpt, label="ConOpt", linewidth = 5.0)
plot!(pl, l10DistOGDA, label="OGDA", linewidth = 5.0)

 savefig(pl, "./out/quadratic_noequi_alpha$α")

max_stars_repo_name: f-t-s/CGD

max_issues_repo_name: f-t-s/CGD

max_forks_repo_name: f-t-s/CGD