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

Status:  support function.

Approximate heliocentric position and velocity of a nominated major
planet:  Mercury, Venus, EMB, Mars, Jupiter, Saturn, Uranus or
Neptune (but not the Earth itself).

Given:
   date1  double       TDB date part A (Note 1)
   date2  double       TDB date part B (Note 1)
   np     int          planet (1=Mercury, 2=Venus, 3=EMB, 4=Mars,
                           5=Jupiter, 6=Saturn, 7=Uranus, 8=Neptune)

Returned (argument):
   pv     double[2][3] planet p,v (heliocentric, J2000.0, au,au/d)

Returned (function value):
         int          status: -1 = illegal NP (outside 1-8)
                                 0 = OK
                              +1 = warning: year outside 1000-3000
                              +2 = warning: failed to converge

Notes:

   1. The date date1+date2 is in the TDB time scale (in practice TT can
   be used) and is a Julian Date, apportioned in any convenient way
   between the two arguments.  For example, JD(TDB)=2450123.7 could
   be expressed in any of these ways, among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 method is best matched to the way the
   argument is handled internally and will deliver the optimum
   resolution.  The MJD method and the date & time methods are both
   good compromises between resolution and convenience.  The limited
   accuracy of the present algorithm is such that any of the methods
   is satisfactory.

   2. If an np value outside the range 1-8 is supplied, an error status
   (function value -1) is returned and the pv vector set to zeroes.

   3. For np=3 the result is for the Earth-Moon Barycenter.  To obtain
   the heliocentric position and velocity of the Earth, use instead
   the SOFA function iauEpv00.

   4. On successful return, the array pv contains the following:

      pv[0][0]   x      }
      pv[0][1]   y      } heliocentric position, au
      pv[0][2]   z      }

      pv[1][0]   xdot   }
      pv[1][1]   ydot   } heliocentric velocity, au/d
      pv[1][2]   zdot   }

   The reference frame is equatorial and is with respect to the
   mean equator and equinox of epoch J2000.0.

   5. The algorithm is due to J.L. Simon, P. Bretagnon, J. Chapront,
   M. Chapront-Touze, G. Francou and J. Laskar (Bureau des
   Longitudes, Paris, France).  From comparisons with JPL
   ephemeris DE102, they quote the following maximum errors
   over the interval 1800-2050:

                  L (arcsec)    B (arcsec)      R (km)

      Mercury          4             1             300
      Venus            5             1             800
      EMB              6             1            1000
      Mars            17             1            7700
      Jupiter         71             5           76000
      Saturn          81            13          267000
      Uranus          86             7          712000
      Neptune         11             1          253000

   Over the interval 1000-3000, they report that the accuracy is no
   worse than 1.5 times that over 1800-2050.  Outside 1000-3000 the
   accuracy declines.

   Comparisons of the present function with the JPL DE200 ephemeris
   give the following RMS errors over the interval 1960-2025:

                     position (km)     velocity (m/s)

      Mercury            334               0.437
      Venus             1060               0.855
      EMB               2010               0.815
      Mars              7690               1.98
      Jupiter          71700               7.70
      Saturn          199000              19.4
      Uranus          564000              16.4
      Neptune         158000              14.4

   Comparisons against DE200 over the interval 1800-2100 gave the
   following maximum absolute differences.  (The results using
   DE406 were essentially the same.)

                  L (arcsec)   B (arcsec)     R (km)   Rdot (m/s)

      Mercury        7            1            500       0.7
      Venus          7            1           1100       0.9
      EMB            9            1           1300       1.0
      Mars          26            1           9000       2.5
      Jupiter       78            6          82000       8.2
      Saturn        87           14         263000      24.6
      Uranus        86            7         661000      27.4
      Neptune       11            2         248000      21.4

   6. The present SOFA re-implementation of the original Simon et al.
   Fortran code differs from the original in the following respects:

      *  C instead of Fortran.

      *  The date is supplied in two parts.

      *  The result is returned only in equatorial Cartesian form;
         the ecliptic longitude, latitude and radius vector are not
         returned.

      *  The result is in the J2000.0 equatorial frame, not ecliptic.

      *  More is done in-line: there are fewer calls to subroutines.

      *  Different error/warning status values are used.

      *  A different Kepler's-equation-solver is used (avoiding
         use of double precision complex).

      *  Polynomials in t are nested to minimize rounding errors.

      *  Explicit double constants are used to avoid mixed-mode
         expressions.

   None of the above changes affects the result significantly.

   7. The returned status indicates the most serious condition
   encountered during execution of the function.  Illegal np is
   considered the most serious, overriding failure to converge,
   which in turn takes precedence over the remote date warning.

Called:
   iauAnp       normalize angle into range 0 to 2pi

Reference:  Simon, J.L, Bretagnon, P., Chapront, J.,
            Chapront-Touze, M., Francou, G., and Laskar, J.,
            Astron.Astrophys., 282, 663 (1994).

This revision:  2017 October 12

SOFA release 2018-01-30

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

# int iauPlan94(double date1, double date2, int np, double pv[2][3])

function iauPlan94(date1::Real, date2::Real, np::Int)
    pv = zeros(Float64, 3, 2)

    status = ccall((:iauPlan94, libsofa_c), Cint,
            (Cdouble, Cdouble, Cint, Ptr{Cdouble}),
            convert(Float64, date1),
            convert(Float64, date2),
            convert(Int32, np),
            pv)

    return status, SMatrix{2,3}(pv')
end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

"Supertype for clocks."
abstract type AbstractClock end

"Supertype for events."
abstract type AbstractEvent end

# abstract types for clock state machines
abstract type ClockState end
abstract type ClockEvent end

"""
    Timing

Enumeration type for scheduling events and timed conditions:

- `at`: schedule an event at a given time,
- `after`: schedule an event a given time after current time,
- `every`: schedule an event every given time from now on,
- `before`: a timed condition is true before a given time,
- `until`: delay until t.
"""
@enum Timing at after every before until

"""
    Action

An action is either  a `Function` or an `Expr` or a `Tuple` of them.
It can be scheduled in an event for later execution.
"""
const Action = Union{Function, Expr, Tuple{Vararg{Union{Function, Expr}}}}

"""
````
DiscreteEvent{T<:Action, X} <: AbstractEvent
```
A discrete event is an [`Action`](@ref) to be executed at an event time.

# Arguments, fields
- `ex::T`: a function or an expression or a tuple of them,
- `t::Float64`: event time,
- `Δt::X`: repeat interval (Float64, Distribution or Nothing),
- `n::Int`: number of repeats.
"""
struct DiscreteEvent{T<:Action, X} <: AbstractEvent
    ex::T
    t::Float64
    Δt::X
    n::Int
end

"""
    DiscreteCond{S<:Action, T<:Action} <: AbstractEvent

A condition to be evaluated repeatedly with expressions or functions
to be executed if conditions are met.

# Arguments, fields
- `cond::S`: a conditional function or an expression or a tuple of them
    (conditions must evaluate to `Bool`),
- `ex::T`: a function or an expression or a tuple of them to be executed
    if conditions are met,
"""
struct DiscreteCond{S<:Action, T<:Action} <: AbstractEvent
    cond::S
    ex::T
end

"""
    Sample{T<:Action} <: AbstractEvent

Sampling actions are executed at sampling time.

# Arguments, fields
- `ex<:Action`: an [`Action`](@ref) to be executed at sample time.
"""
struct Sample{T<:Action} <: AbstractEvent
    ex::T
end

"""
    PrcException(event, value)

An exception to be thrown at processes.

# Arguments, fields
- `event`: deliver an event to the interrupted task
- `value`: deliver some other value
"""
struct PrcException{U,V} <: Exception
  event::U
  value::V
end

"""
```
Prc(id, f, arg...; kw...)
Prc(f, arg...; kw...)
```

Prepare a function to run as a process and get registered to a clock.

# Arguments, Fields
- `id`: some unique identification for registration,
- `f::Function`: a function `f(clk, arg...; kw...)`, must take `clk`
    (a [`Clock`](@ref)) as its first argument,
- `arg...`: further arguments to `f`
- `kw...`: keyword arguments to `f`

# Fields identified during registration
- `id`: if it is not provided, some integer will be calculated
    for it during registration,
- `task::Union{Task,Nothing}`: a task structure used for diagnosis,
- `clk::Union{AbstractClock,Nothing}`: clock where the process is
    registered,

!!! note "The clock is identified dynamically!"

    The clock `clk` where a process runs and gets registered is
    identified during process startup and then passed as 1st
    argument to `f`.
"""
mutable struct Prc
    id::Any
    task::Union{Task,Nothing}
    clk::Union{AbstractClock,Nothing}
    f::Function
    arg::Tuple
    kw::Base.Iterators.Pairs

    Prc( id, f::Function, arg...; kw...) = new(id, nothing, nothing, f, arg, kw)
    Prc(f::Function, arg...; kw...) = new(1, nothing, nothing, f, arg, kw)
end

"""
    Schedule()

A schedule contains events, conditional events and sampling functions
to be executed or evaluated on the clock's time line.

# Fields
- `events::PriorityQueue{DiscreteEvent,Float64}`: scheduled events,
- `cevents::Array{DiscreteCond,1}`: conditional events to evaluate at each tick,
- `samples::Array{Sample,1}`: sampling expressions to evaluate at each tick,
"""
mutable struct Schedule
    events::PriorityQueue{DiscreteEvent,Float64}
    cevents::Array{DiscreteCond,1}
    samples::Array{Sample,1}

    Schedule() = new(PriorityQueue{DiscreteEvent,Float64}(), DiscreteCond[], Sample[])
end

"""
    ClockChannel{T <: ClockEvent}

Provide a message channel to an active clock or a real time clock.

# Fields
- `ref::Ref{Task}`: a reference to an active clock task, useful for diagnosis,
- `forth::Channel{T}`: a communication channel to an active clock,
- `back::Channel{T}`: response channel from the active clock,
- `thread::Int`: the thread id of the active clock,
- `done::Bool`: flag indicating if the active clock has completed its cycle,
- `load::Int`: internal flag.
"""
mutable struct ClockChannel{T <: ClockEvent}
    ref::Ref{Task}
    forth::Channel{T}
    back::Channel{T}
    thread::Int
    done::Bool
    load::Int
end

"""
    Clock{AC} <: AbstractClock

A virtual clock structure, used for global and thread local clocks.

# Fields
- `id::Int`: clock ident number 1: master clock, > 1: parallel clock,
- `ac::AC`: if id == 1: a [`Vector{ClockChannel}`](@ref ClockChannel) else: [`Ref{ActiveClock}`](@ref ActiveClock),
- `state::ClockState`: clock state,
- `time::Float64`: clock time,
- `unit::FreeUnits`: time unit,
- `end_time::Float64`: end time for simulation,
- `Δt::Float64`: sampling time, timestep between ticks,
- `sc::Schedule`: the clock [`Schedule`](@ref) (events, cond events and sampling),
- `processes::Dict{Any, Prc}`: registered `Prc`es,
- `channels::Vector{Channel}`: registered (Actor) channels,
- `tn::Float64`: next timestep,
- `tev::Float64`: next event time,
- `evcount::Int`: event counter,
- `scount::Int`: sample counter
"""
mutable struct Clock{AC} <: AbstractClock
    id::Int
    ac::AC
    state::ClockState
    time::Float64
    unit::FreeUnits
    Δt::Float64
    sc::Schedule
    processes::Dict{Any, Prc}
    channels::Vector{Channel}
    tn::Float64
    tev::Float64
    end_time::Float64
    evcount::Int
    scount::Int
end

"Global clock with vector of active clocks."
const GlobalClock = Clock{Vector{ClockChannel}}

"""
```
Clock(Δt::T=0.01; t0::U=0, unit::FreeUnits=NoUnits) where {T<:Number,U<:Number}
```
Create a new virtual clock.

# Arguments
- `Δt::T=0.01`: time increment for sampling. Δt can be set later with `sample_time!`.
- `t0::U=0`: start time for simulation
- `unit::FreeUnits=NoUnits`: clock time unit. Units can be set explicitely by
    setting e.g. `unit=minute` or implicitly by giving Δt as a time or else setting
    t0 to a time, e.g. `t0=60s`.

The created clock has id==1 (master).
"""
function Clock(Δt::T=0.01;
    t0::U=0, unit::FreeUnits=NoUnits) where {T<:Number,U<:Number}
    if 1unit isa Time
        Δt = isa(Δt, Time) ? uconvert(unit, Δt).val : Δt
        t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
    elseif Δt isa Time
        unit = Unitful.unit(Δt)
        t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
        Δt = Δt.val
    elseif t0 isa Time
        unit = Unitful.unit(t0)
        t0 = t0.val
    else
        nothing
    end
    GlobalClock(threadid(), ClockChannel[], Idle(),
        t0, unit, Δt, Schedule(), Dict{Any, Prc}(), Channel[],
        t0 + Δt, t0, t0, 0, 0)
end

"""
    ActiveClock{E <: ClockEvent} <: AbstractClock

An active clock is a wrapper around a local [`Clock`](@ref) on a
parallel thread. It is operated by a thread local task. The master
clock on thread 1 communicates with it through messages over its
channels.

# Fields
- `clock::Clock`: the thread specific local clock,
- `master::Ref{Clock}`: a pointer to the master clock (on thread 1),
- `forth::Channel{E}`: the command channel from master,
- `back::Channel{E}`: the response channel to master,
- `id::Int`: the clocks id/thread number,
- `task::Task`: the active clock`s task.

!!! note "Don't setup an active clock explicitly!"

    It is done implicitly with [`PClock`](@ref) or by [`fork!`](@ref)ing
    a [`Clock`](@ref) to other available threads.

On a parallel thread tasks can access their local clock with
[`pclock(clk)`](@ref pclock). Then they can schedule 'thread local'
events, `delay!` or `wait!` on it.

!!! note "Don't share active clocks between threads!"

    In multithreading we communicate over channels and don't want
    to share variables between threads. A user can still access
    active clocks for diagnostic purposes.
"""
mutable struct ActiveClock{E <: ClockEvent} <: AbstractClock
    clock::Clock
    master::Ref{Clock{Vector{ClockChannel}}}
    forth::Channel{E}
    back ::Channel{E}
    id::Int
    task::Task
end

"Local clock with reference to a wrapping active clock."
const LocalClock = Clock{Ref{ActiveClock{ClockEvent}}}

"""
    localClock(c::Clock)

Create a thread local clock with an undefined reference to an
active clock which inherits parameters of the master clock `c`.
"""
function localClock(c::Clock)
    LocalClock(threadid(),
        Ref{ActiveClock{ClockEvent}}(), Idle(),
        c.time, c.unit, c.Δt, Schedule(), Dict{Any, Prc}(), Channel[],
        c.time + c.Δt, c.time, c.time, 0, 0)
end

"""
    RTClock{E <: ClockEvent} <: AbstractClock

A real time clock checks every given period for scheduled events and executes them. It has a
time in seconds since its start or last reset and uses system time for updating.

Real time clocks are controlled over channels. Multiple real time clocks can be setup with
arbitrary periods (≥ 1 ms). Real time clocks should not be created directly but rather with
`CreateRTClock`.

# Fields
- `Timer::Timer`:     clock period in seconds, minimum is 0.001 (1 ms)
- `clock::Clock`:     clock work
- `cmd::Channel{T}`:  command channel to asynchronous clock
- `back::Channel{T}`: back channel from async clock
- `id::Int`:          arbitrary id number
- `thread::Int`:      thread the async clock is living in
- `time::Float64`:    clock time since start in seconds
- `t0::Float64`:      system time at clock start in seconds
- `T::Float64`:       clock period in seconds
"""
mutable struct RTClock{E <: ClockEvent} <: AbstractClock
    Timer::Timer
    clock::Clock
    cmd::Channel{E}
    back::Channel{E}
    id::Int
    thread::Int
    time::Float64
    t0::Float64
    T::Float64
    task::Task
end

ext: jl


text:
abstract PGAbstractImage

type PGCatalog <: PGAbstractImage
    name2schema::Dict{String,UInt}

    PGCatalog() = new(Dict{String,UInt}())
end

type PGSchema <: PGAbstractImage
    name::String
    name2table::Dict{String,UInt}
    name2index::Dict{String,UInt}
    name2sequence::Dict{String,UInt}
    sig2procedure::Dict{Tuple{String,Tuple{Vararg{UInt}}},UInt}
    name2type::Dict{String,UInt}
    comment::Nullable{String}

    PGSchema(name) =
        new(name,
            Dict{String,UInt}(),
            Dict{String,UInt}(),
            Dict{String,UInt}(),
            Dict{Tuple{String,Tuple{Vararg{UInt}}},UInt}(),
            Dict{String,UInt}(),
            nothing)
end

type PGTable <: PGAbstractImage
    schema::UInt
    name::String
    columns::Vector{UInt}
    name2column::Dict{String,UInt}
    name2ukey::Dict{String,UInt}
    pkey::Nullable{UInt}
    name2fkey::Dict{String,UInt}
    rfkeys::Set{UInt}
    indexes::Set{UInt}
    name2trigger::Dict{String,UInt}
    comment::Nullable{String}

    PGTable(schema, name) =
        new(schema,
            name,
            UInt[],
            Dict{String,UInt}(),
            Dict{String,UInt}(),
            nothing,
            Dict{String,UInt}(),
            Set{UInt}(),
            Set{UInt}(),
            Dict{String,UInt}(),
            nothing)
end

type PGColumn <: PGAbstractImage
    table::UInt
    name::String
    type_::UInt
    isnull::Bool
    default::Nullable{String}
    ukeys::Set{UInt}
    fkeys::Set{UInt}
    rfkeys::Set{UInt}
    indexes::Set{UInt}
    sequences::Set{UInt}
    comment::Nullable{String}

    PGColumn(table, name, type_, isnull=false, default=nothing) =
        new(table,
            name,
            type_,
            isnull,
            default,
            Set{UInt}(),
            Set{UInt}(),
            Set{UInt}(),
            Set{UInt}(),
            Set{UInt}(),
            nothing)
end

type PGUniqueKey <: PGAbstractImage
    name::String
    table::UInt
    columns::Vector{UInt}
    isprimary::Bool
    indexes::Set{UInt}
    comment::Nullable{String}

    PGUniqueKey(name, table, columns, isprimary=false) =
        new(name,
            table,
            columns,
            isprimary,
            Set{UInt}(),
            nothing)
end

type PGForeignKey <: PGAbstractImage
    name::String
    table::UInt
    columns::Vector{UInt}
    rtable::UInt
    rcolumns::Vector{UInt}
    onupdate::String
    ondelete::String
    comment::Nullable{String}

    PGForeignKey(name, table, columns, rtable, rcolumns, onupdate="NO ACTION", ondelete="NO ACTION") =
        new(name,
            table,
            columns,
            rtable,
            rcolumns,
            onupdate,
            ondelete,
            nothing)
end

type PGIndex <: PGAbstractImage
    schema::UInt
    name::String
    table::UInt
    columns::Vector{UInt}
    ukey::Nullable{UInt}
    comment::Nullable{String}

    PGIndex(schema, name, table, columns) =
        new(schema,
            name,
            table,
            columns,
            nothing,
            nothing)
end

type PGSequence <: PGAbstractImage
    schema::UInt
    name::String
    column::Nullable{UInt}
    comment::Nullable{String}

    PGSequence(schema, name) =
        new(schema,
            name,
            nothing,
            nothing)
end

type PGProcedure <: PGAbstractImage
    schema::UInt
    name::String
    argtypes::Vector{UInt}
    restype::UInt
    source::String
    triggers::Set{UInt}
    comment::Nullable{String}

    PGProcedure(schema, name, argtypes, restype, source) =
        new(schema,
            name,
            argtypes,
            restype,
            source,
            Set{UInt}(),
            nothing)
end

type PGTrigger <: PGAbstractImage
    table::UInt
    name::String
    procedure::UInt
    comment::Nullable{String}

    PGTrigger(table, name, procedure) =
        new(table,
            name,
            procedure,
            nothing)
end

type PGType <: PGAbstractImage
    schema::UInt
    name::String
    enum::Nullable{Vector{String}}
    columns::Set{UInt}
    procedures::Set{UInt}
    comment::Nullable{String}

    PGType(schema, name, enum=nothing) =
        new(schema,
            name,
            enum,
            Set{UInt}(),
            Set{UInt}(),
            nothing)
end

type PGImageSet
    catalog::PGCatalog
    schemas::Dict{UInt,PGSchema}
    tables::Dict{UInt,PGTable}
    columns::Dict{UInt,PGColumn}
    ukeys::Dict{UInt,PGUniqueKey}
    fkeys::Dict{UInt,PGForeignKey}
    indexes::Dict{UInt,PGIndex}
    sequences::Dict{UInt,PGSequence}
    procedures::Dict{UInt,PGProcedure}
    triggers::Dict{UInt,PGTrigger}
    types::Dict{UInt,PGType}

    PGImageSet() =
        new(PGCatalog(),
            Dict{UInt,PGSchema}(),
            Dict{UInt,PGTable}(),
            Dict{UInt,PGColumn}(),
            Dict{UInt,PGUniqueKey}(),
            Dict{UInt,PGForeignKey}(),
            Dict{UInt,PGIndex}(),
            Dict{UInt,PGSequence}(),
            Dict{UInt,PGProcedure}(),
            Dict{UInt,PGTrigger}(),
            Dict{UInt,PGType}())
end

immutable PGImage{T}
    set::PGImageSet
    obj::T
end

convert(::Type{UInt}, img::PGImage) =
    object_id(img.obj)

convert(::Type{UInt}, obj::PGAbstractImage) =
    object_id(obj)

PGImage(::Type{PGCatalog}, set::PGImageSet) =
    PGImage{PGCatalog}(set, set.catalog)

newcatalog() =
    PGImage(PGCatalog, PGImageSet())

function addschema(catalog::PGImage{PGCatalog}, name::String)
    obj = PGSchema(name)
    oid = object_id(obj)
    catalog.set.schemas[oid] = obj
    catalog.obj.name2schema[obj.name] = oid
    PGImage(catalog.set, obj)
end

function addtable(schema::PGImage{PGSchema}, name::String)
    obj = PGTable(schema, name)
    oid = object_id(obj)
    schema.set.tables[oid] = obj
    schema.obj.name2table[obj.name] = oid
    PGImage(schema.set, obj)
end

function addtype(schema::PGImage{PGSchema}, name::String)
    obj = PGType(schema, name)
    oid = object_id(obj)
    schema.set.types[oid] = obj
    schema.obj.name2type[obj.name] = oid
    PGImage(schema.set, obj)
end

function addsequence(schema::PGImage{PGSchema}, name::String)
    obj = PGSequence(schema, name)
    oid = object_id(obj)
    schema.set.sequences[oid] = obj
    schema.obj.name2sequence[obj.name] = oid
    PGImage(schema.set, obj)
end

function addtype(schema::PGImage{PGSchema}, name::String, labels::Vector{String})
    obj = PGType(schema, name, labels)
    oid = object_id(obj)
    schema.set.types[oid] = obj
    schema.obj.name2type[obj.name] = oid
    PGImage(schema.set, obj)
end

function addprocedure(
        schema::PGImage{PGSchema}, name::String, argtypes::Vector{PGImage{PGType}}, restype::PGImage{PGType},
        source::String)
    argtypeids = [object_id(argtype.obj) for argtype in argtypes]
    obj = PGProcedure(schema, name, argtypeids, object_id(restype.obj), source)
    oid = object_id(obj)
    schema.set.procedures[oid] = obj
    schema.obj.sig2procedure[(name, (argtypeids...))] = oid
    for type_ in argtypes
        push!(type_.obj.procedures, oid)
    end
    push!(restype.obj.procedures, oid)
    PGImage(schema.set, obj)
end

function addcolumn(table::PGImage{PGTable}, name::String, type_::PGImage{PGType}, isnull::Bool)
    obj = PGColumn(table, name, type_, isnull)
    oid = object_id(obj)
    table.set.columns[oid] = obj
    table.obj.name2column[name] = oid
    push!(table.obj.columns, oid)
    push!(type_.obj.columns, oid)
    PGImage(table.set, obj)
end

function addindex(schema::PGImage{PGSchema}, name::String, table::PGImage{PGTable}, columns::Vector{PGImage{PGColumn}})
    columnids = [object_id(column.obj) for column in columns]
    obj = PGIndex(schema, name, table, columns)
    oid = object_id(obj)
    schema.set.indexes[oid] = obj
    schema.obj.name2index[name] = oid
    push!(table.obj.indexes, oid)
    for column in columns
        push!(column.obj.indexes, oid)
    end
    PGImage(schema.set, obj)
end

function adduniquekey(table::PGImage{PGTable}, name::String, columns::Vector{PGImage{PGColumn}}, isprimary::Bool)
    columnids = [object_id(column.obj) for column in columns]
    obj = PGUniqueKey(name, table, columnids, isprimary)
    oid = object_id(obj)
    table.set.ukeys[oid] = obj
    table.obj.name2ukey[name] = oid
    if isprimary
        table.obj.pkey = oid
    end
    for column in columns
        push!(column.obj.ukeys, oid)
    end
    PGImage(table.set, obj)
end

function addforeignkey(
        table::PGImage{PGTable}, name::String, columns::Vector{PGImage{PGColumn}},
        rtable::PGImage{PGTable}, rcolumns::Vector{PGImage{PGColumn}}, onupdate::String, ondelete::String)
    columnids = [object_id(column.obj) for column in columns]
    rcolumnids = [object_id(column.obj) for column in rcolumns]
    obj = PGForeignKey(name, table, columnids, rtable, rcolumnids, onupdate, ondelete)
    oid = object_id(obj)
    table.set.fkeys[oid] = obj
    table.obj.name2fkey[name] = oid
    push!(rtable.obj.rfkeys, oid)
    for column in columns
        push!(column.obj.fkeys, oid)
    end
    for column in rcolumns
        push!(column.obj.rfkeys, oid)
    end
    PGImage(table.set, obj)
end

function addtrigger(table::PGImage{PGTable}, name::String, procedure::PGImage{PGProcedure})
    obj = PGTrigger(table, name, procedure)
    oid = object_id(obj)
    table.set.triggers[oid] = obj
    table.obj.name2trigger[name] = oid
    push!(procedure.obj.triggers, oid)
    PGImage(table.set, obj)
end

function setindex(ukey::PGImage{PGUniqueKey}, index::PGImage{PGIndex})
    index.obj.ukey = ukey
    push!(ukey.obj.indexes, object_id(index.obj))
    index
end

function setdefault(column::PGImage{PGColumn}, default::String)
    column.obj.default = default
    column
end

function setcolumn(sequence::PGImage{PGSequence}, column::PGImage{PGColumn})
    sequence.obj.column = column
    push!(column.obj.sequences, object_id(sequence.obj))
    sequence
end

function setcomment(schema::PGImage{PGSchema}, comment::String)
    schema.obj.comment = comment
    schema
end

function setcomment(type_::PGImage{PGType}, comment::String)
    type_.obj.comment = comment
    type_
end

function setcomment(table::PGImage{PGTable}, comment::String)
    table.obj.comment = comment
    table
end

function setcomment(index::PGImage{PGIndex}, comment::String)
    index.obj.comment = comment
    index
end

function setcomment(sequence::PGImage{PGSequence}, comment::String)
    sequence.obj.comment = comment
    sequence
end

function setcomment(column::PGImage{PGColumn}, comment::String)
    column.obj.comment = comment
    column
end

function setcomment(ukey::PGImage{PGUniqueKey}, comment::String)
    ukey.obj.comment = comment
    ukey
end

function setcomment(fkey::PGImage{PGForeignKey}, comment::String)
    fkey.obj.comment = comment
    fkey
end

function setcomment(trigger::PGImage{PGTrigger}, comment::String)
    trigger.obj.comment = comment
    trigger
end

function introspect(conn::PGConnection)
    catalog = newcatalog()
    # Extract schemas.
    oid2schema = Dict{UInt32,PGImage{PGSchema}}()
    execute(conn, """
        SELECT n.oid, n.nspname
        FROM pg_catalog.pg_namespace n
        ORDER BY n.nspname
    """) do row
        oid, nspname = map(get, row)
        schema = addschema(catalog, nspname)
        oid2schema[oid] = schema
    end
    # Extract ENUM labels.
    oid2labels = Dict{UInt32,Vector{String}}()
    execute(conn, """
        SELECT e.enumtypid, e.enumlabel
        FROM pg_catalog.pg_enum e
        ORDER BY e.enumtypid, e.enumsortorder, e.oid
    """) do row
        enumtypid, enumlabel = map(get, row)
        if !(enumtypid in keys(oid2labels))
            oid2labels[enumtypid] = String[]
        end
        push!(oid2labels[enumtypid], enumlabel)
    end
    # Extract data types.
    oid2type = Dict{UInt32,PGImage{PGType}}()
    execute(conn, """
        SELECT t.oid, t.typnamespace, t.typname, t.typtype,
               t.typbasetype, t.typlen, t.typtypmod
        FROM pg_catalog.pg_type t
        ORDER BY t.typnamespace, t.typname
    """) do row
        oid, typnamespace, typname, typtype, typbasetype, typlen, typtypmod = map(get, row)
        schema = oid2schema[typnamespace]
        type_ =
            if typtype == 'e'
                labels = oid2labels[oid]
                addtype(schema, typname, labels)
            else
                addtype(schema, typname)
            end
        oid2type[oid] = type_
    end
    # Extract stored procedures.
    oid2procedure = Dict{UInt32,PGImage{PGProcedure}}()
    execute(conn, """
        SELECT p.oid, p.pronamespace, p.proname,
               p.proargtypes, p.prorettype, p.prosrc
        FROM pg_catalog.pg_proc p
        ORDER BY p.pronamespace, p.proname
    """) do row
        oid, pronamespace, proname, proargtypes, prorettype, prosrc = map(get, row)
        schema = oid2schema[pronamespace]
        argtypes = PGImage{PGType}[oid2type[proargtype] for proargtype in proargtypes]
        rettype = oid2type[prorettype]
        procedure = addprocedure(schema, proname, argtypes, rettype, prosrc)
        oid2procedure[oid] = procedure
    end
    # Extract tables.
    oid2table = Dict{UInt32,PGImage{PGTable}}()
    execute(conn, """
        SELECT c.oid, c.relnamespace, c.relname
        FROM pg_catalog.pg_class c
        WHERE c.relkind IN ('r', 'v') AND
              HAS_TABLE_PRIVILEGE(c.oid, 'SELECT')
        ORDER BY c.relnamespace, c.relname
    """) do row
        oid, relnamespace, relname = map(get, row)
        schema = oid2schema[relnamespace]
        table = addtable(schema, relname)
        oid2table[oid] = table
    end
    # Extract columns.
    num2column = Dict{Tuple{UInt32,Int16},PGImage{PGColumn}}()
    execute(conn, """
        SELECT a.attrelid, a.attnum, a.attname, a.atttypid, a.atttypmod,
               a.attnotnull, a.atthasdef, a.attisdropped
        FROM pg_catalog.pg_attribute a
        ORDER BY a.attrelid, a.attnum
    """) do row
        attrelid, attnum, attname, atttypid, atttypmod, attnotnull, atthasdef, attisdropped = map(get, row)
        if attisdropped
            return
        end
        if attname in ["tableoid", "cmax", "xmax", "cmin", "xmin", "ctid"]
            return
        end
        if !(attrelid in keys(oid2table))
            return
        end
        table = oid2table[attrelid]
        type_ = oid2type[atttypid]
        isnull = !attnotnull
        column = addcolumn(table, attname, type_, isnull)
        num2column[(attrelid, attnum)] = column
    end
    # Extract default values.
    execute(conn, """
        SELECT a.adrelid, a.adnum, pg_get_expr(a.adbin, a.adrelid)
        FROM pg_catalog.pg_attrdef a
        ORDER BY a.adrelid, a.adnum
    """) do row
        adrelid, adnum, adsrc = map(get, row)
        if (adrelid, adnum) in keys(num2column)
            column = num2column[(adrelid, adnum)]
            setdefault(column, adsrc)
        end
    end
    # Extract sequences.
    oid2sequence = Dict{UInt32,PGImage{PGSequence}}()
    execute(conn, """
        SELECT c.oid, c.relnamespace, c.relname
        FROM pg_catalog.pg_class c
        WHERE c.relkind = 'S'
        ORDER BY c.relnamespace, c.relname
    """) do row
        oid, relnamespace, relname = map(get, row)
        schema = oid2schema[relnamespace]
        sequence = addsequence(schema, relname)
        oid2sequence[oid] = sequence
    end
    # Associate sequences with the columns that own them.
    execute(conn, """
        SELECT d.objid, d.refobjid, d.refobjsubid
        FROM pg_catalog.pg_depend d
        JOIN pg_catalog.pg_class c
        ON (d.classid = 'pg_class'::regclass AND d.objid = c.oid)
        WHERE c.relkind = 'S' AND
              d.refclassid = 'pg_class'::regclass AND
              d.objsubid IS NOT NULL
        ORDER BY d.objid, d.refobjid, d.objsubid
    """) do row
        objid, refobjid, refobjsubid = map(get, row)
        sequence = oid2sequence[objid]
        column = num2column[(refobjid, refobjsubid)]
        setcolumn(sequence, column)
    end
    # Extract indexes.
    oid2index = Dict{UInt32,PGImage{PGIndex}}()
    execute(conn, """
        SELECT c.oid, c.relnamespace, c.relname, i.indrelid, i.indkey
        FROM pg_catalog.pg_class c
        JOIN pg_catalog.pg_index i
        ON (c.oid = i.indexrelid)
        WHERE c.relkind = 'i'
        ORDER BY c.relnamespace, c.relname
    """) do row
        oid, relnamespace, relname, indrelid, indkey = map(get, row)
        if !(indrelid in keys(oid2table))
            return
        end
        schema = oid2schema[relnamespace]
        table = oid2table[indrelid]
        columns = PGImage{PGColumn}[]
        for num in indkey
            if (indrelid, num) in keys(num2column)
                push!(columns, num2column[(indrelid, num)])
            end
        end
        index = addindex(schema, relname, table, columns)
        oid2index[oid] = index
    end
    # Extract unique keys.
    oid2ukey = Dict{UInt32,PGImage{PGUniqueKey}}()
    execute(conn, """
        SELECT c.oid, c.conname, c.contype, c.conrelid, c.conkey
        FROM pg_catalog.pg_constraint c
        WHERE c.contype IN ('p', 'u')
        ORDER BY c.oid
    """) do row
        oid, conname, contype, conrelid, conkey = map(get, row)
        if !(conrelid in keys(oid2table))
            return
        end
        table = oid2table[conrelid]
        columns = [num2column[(conrelid, num)] for num in conkey]
        isprimary = contype == 'p'
        key = adduniquekey(table, conname, columns, isprimary)
        oid2ukey[oid] = key
        schema = table.set.schemas[table.obj.schema]
        index = PGImage(table.set, table.set.indexes[schema.name2index[conname]])
        setindex(key, index)
    end
    # Extract foreign keys.
    oid2fkey = Dict{UInt32, PGImage{PGForeignKey}}()
    execute(conn, """
        SELECT c.oid, c.conname, c.conrelid, c.conkey,
               c.confrelid, c.confkey, c.confupdtype, c.confdeltype
        FROM pg_catalog.pg_constraint c
        WHERE c.contype = 'f'
        ORDER BY c.oid
    """) do row
        oid, conname, conrelid, conkey, confrelid, confkey, confupdtype, confdeltype = map(get, row)
        if !(conrelid in keys(oid2table)) || !(confrelid in keys(oid2table))
            return
        end
        table = oid2table[conrelid]
        columns = [num2column[(conrelid, num)] for num in conkey]
        rtable = oid2table[confrelid]
        rcolumns = [num2column[(confrelid, num)] for num in confkey]
        onupdate =
            confupdtype == 'a' ? "NO ACTION" :
            confupdtype == 'r' ? "RESTRICT" :
            confupdtype == 'c' ? "CASCADE" :
            confupdtype == 'n' ? "SET NULL" :
            confupdtype == 'd' ? "SET DEFAULT" : ""
        ondelete =
            confdeltype == 'a' ? "NO ACTION" :
            confdeltype == 'r' ? "RESTRICT" :
            confdeltype == 'c' ? "CASCADE" :
            confdeltype == 'n' ? "SET NULL" :
            confdeltype == 'd' ? "SET DEFAULT" : ""
        fkey = addforeignkey(table, conname, columns, rtable, rcolumns, onupdate, ondelete)
        oid2fkey[oid] = fkey
    end
    # Extract triggers.
    oid2trigger = Dict{UInt32,PGImage{PGTrigger}}()
    execute(conn, """
        SELECT t.oid, t.tgrelid, t.tgname, t.tgfoid
        FROM pg_catalog.pg_trigger AS t
        WHERE NOT t.tgisinternal
        ORDER BY t.tgrelid, t.tgname
    """) do row
        oid, tgrelid, tgname, tgfoid = map(get, row)
        table = oid2table[tgrelid]
        procedure = oid2procedure[tgfoid]
        trigger = addtrigger(table, tgname, procedure)
        oid2trigger[oid] = trigger
    end
    # Extract comments.
    execute(conn, """
        SELECT c.relname, d.objoid, d.objsubid, d.description
        FROM pg_catalog.pg_description d
        JOIN pg_catalog.pg_class c
        ON (d.classoid = c.oid)
        WHERE d.classoid IN ('pg_catalog.pg_namespace'::regclass,
                             'pg_catalog.pg_type'::regclass,
                             'pg_catalog.pg_proc'::regclass,
                             'pg_catalog.pg_class'::regclass,
                             'pg_catalog.pg_constraint'::regclass,
                             'pg_catalog.pg_trigger'::regclass)
        ORDER BY d.objoid, d.classoid, d.objsubid
    """) do row
        relname, objoid, objsubid, description = map(get, row)
        if relname == "pg_namespace"
            schema = oid2schema[objoid]
            setcomment(schema, description)
        elseif relname == "pg_type"
            type_ = oid2type[objoid]
            setcomment(type_, description)
        elseif relname == "pg_class" && objsubid == 0 && objoid in keys(oid2table)
            table = oid2table[objoid]
            setcomment(table, description)
        elseif relname == "pg_class" && objsubid == 0 && objoid in keys(oid2index)
            index = oid2index[objoid]
            setcomment(index, description)
        elseif relname == "pg_class" && objsubid == 0 && objoid in keys(oid2sequence)
            sequence = oid2sequence[objoid]
            setcomment(sequence, description)
        elseif relname == "pg_class" && (objoid, objsubid) in keys(num2column)
            column = num2column[(objoid, objsubid)]
            setcomment(column, description)
        elseif relname == "pg_constraint" && objoid in keys(oid2ukey)
            ukey = oid2ukey[objoid]
            setcomment(ukey, description)
        elseif relname == "pg_constraint" && objoid in keys(oid2fkey)
            fkey = oid2fkey[objoid]
            setcomment(fkey, description)
        elseif relname == "pg_trigger" && objoid in keys(oid2trigger)
            trigger = oid2trigger[objoid]
            setcomment(trigger, description)
        end
    end
    return catalog
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
"""
    load_sbml_model(file_name::String)::SBMLModel

Load and return a SBML XML model in `file_name`.
"""
function load_sbml_model(file_name::String)::SBMLModel
    return SBMLModel(SBML.readSBML(file_name))
end


text:
function create_besselian_epoch(obj::TimeScales, arg0::jdouble)
    return jcall(obj, "createBesselianEpoch", AbsoluteDate, (jdouble,), arg0)
end

function create_julian_epoch(obj::TimeScales, arg0::jdouble)
    return jcall(obj, "createJulianEpoch", AbsoluteDate, (jdouble,), arg0)
end

function get_beidou_epoch(obj::TimeScales)
    return jcall(obj, "getBeidouEpoch", AbsoluteDate, ())
end

function get_ccsds_epoch(obj::TimeScales)
    return jcall(obj, "getCcsdsEpoch", AbsoluteDate, ())
end

function get_fifties_epoch(obj::TimeScales)
    return jcall(obj, "getFiftiesEpoch", AbsoluteDate, ())
end

function get_future_infinity(obj::TimeScales)
    return jcall(obj, "getFutureInfinity", AbsoluteDate, ())
end

function get_galileo_epoch(obj::TimeScales)
    return jcall(obj, "getGalileoEpoch", AbsoluteDate, ())
end

function get_glonass_epoch(obj::TimeScales)
    return jcall(obj, "getGlonassEpoch", AbsoluteDate, ())
end

function get_gmst(obj::TimeScales, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getGMST", GMSTScale, (IERSConventions, jboolean), arg0, arg1)
end

function get_gps_epoch(obj::TimeScales)
    return jcall(obj, "getGpsEpoch", AbsoluteDate, ())
end

function get_irnss_epoch(obj::TimeScales)
    return jcall(obj, "getIrnssEpoch", AbsoluteDate, ())
end

function get_j2000_epoch(obj::TimeScales)
    return jcall(obj, "getJ2000Epoch", AbsoluteDate, ())
end

function get_java_epoch(obj::TimeScales)
    return jcall(obj, "getJavaEpoch", AbsoluteDate, ())
end

function get_julian_epoch(obj::TimeScales)
    return jcall(obj, "getJulianEpoch", AbsoluteDate, ())
end

function get_modified_julian_epoch(obj::TimeScales)
    return jcall(obj, "getModifiedJulianEpoch", AbsoluteDate, ())
end

function get_past_infinity(obj::TimeScales)
    return jcall(obj, "getPastInfinity", AbsoluteDate, ())
end

function get_qzss_epoch(obj::TimeScales)
    return jcall(obj, "getQzssEpoch", AbsoluteDate, ())
end

function get_ut1(obj::TimeScales, arg0::IERSConventions, arg1::jboolean)
    return jcall(obj, "getUT1", UT1Scale, (IERSConventions, jboolean), arg0, arg1)
end

function of(::Type{TimeScales}, arg0::Collection, arg1::BiFunction)
    return jcall(TimeScales, "of", TimeScales, (Collection, BiFunction), arg0, arg1)
end


text:
function [X,J,dXdx]=JacobiAnnulusCyl(ksi,eta,F,Grid)
n1=size(ksi,1);
n2=size(eta,2);
X=zeros(n1,n2,2);
Phi1=atan2(Grid.Nodes(F.N(1)).P(2),Grid.Nodes(F.N(1)).P(1));
Rad1=sqrt(Grid.Nodes(F.N(1)).P(1)^2+Grid.Nodes(F.N(1)).P(2)^2);
Phi2=atan2(Grid.Nodes(F.N(2)).P(2),Grid.Nodes(F.N(2)).P(1));
Rad2=sqrt(Grid.Nodes(F.N(2)).P(1)^2+Grid.Nodes(F.N(2)).P(2)^2);
Phi3=atan2(Grid.Nodes(F.N(3)).P(2),Grid.Nodes(F.N(3)).P(1));
Rad3=sqrt(Grid.Nodes(F.N(3)).P(1)^2+Grid.Nodes(F.N(3)).P(2)^2);
Phi4=atan2(Grid.Nodes(F.N(4)).P(2),Grid.Nodes(F.N(4)).P(1));
Rad4=sqrt(Grid.Nodes(F.N(4)).P(1)^2+Grid.Nodes(F.N(4)).P(2)^2);
if Phi1<0
  Phi1=Phi1+2*pi;
end
if Phi2<0
  Phi2=Phi2+2*pi;
end
if Phi3<0
  Phi3=Phi3+2*pi;
end
if Phi4<0
  Phi4=Phi4+2*pi;
end
if Phi1>Phi3
  Phi3=Phi3+2*pi;
end
if Phi2>Phi4
  Phi4=Phi4+2*pi;
end
Phi=0.25*(Phi1*(1-ksi)*(1-eta)+...
  Phi2*(1+ksi)*(1-eta)+...
  Phi3*(1+ksi)*(1+eta)+...
  Phi4*(1-ksi)*(1+eta));
Rad=0.25*(Rad1*(1-ksi)*(1-eta)+...
  Rad2*(1+ksi)*(1-eta)+...
  Rad3*(1+ksi)*(1+eta)+...
  Rad4*(1-ksi)*(1+eta));
X(:,:,1)=Phi ;
X(:,:,2)=Rad;

if nargout > 2
  dPhidksi=0.25*((-Phi1+Phi2+...
    Phi3-Phi4)*ones(n1,1)*ones(1,n2)+...
    ( Phi1-Phi2+...
    Phi3-Phi4)*ones(n1,1)*eta);
  dPhideta=0.25*((-Phi1-Phi2+...
    Phi3+Phi4)*ones(n1,1)*ones(1,n2)+...
    ( Phi1-Phi2+...
    Phi3-Phi4)*ksi*ones(1,n2));
  dRaddksi=0.25*((-Rad1+Rad2+...
    Rad3-Rad4)*ones(n1,1)*ones(1,n2)+...
    ( Rad1-Rad2+...
    Rad3-Rad4)*ones(n1,1)*eta);
  dRaddeta=0.25*((-Rad1-Rad2+...
    Rad3+Rad4)*ones(n1,1)*ones(1,n2)+...
    ( Rad1-Rad2+...
    Rad3-Rad4)*ksi*ones(1,n2));

  dXdx=zeros(n1,n2,2,2);
  dXdx(:,:,1,1)=dPhidksi;
  dXdx(:,:,2,1)=dRaddksi;
  dXdx(:,:,1,2)=dPhideta;
  dXdx(:,:,2,2)=dRaddeta;
  J=dXdx(:,:,1,1).*dXdx(:,:,2,2)-dXdx(:,:,1,2).*dXdx(:,:,2,1);
end
end


text:
"""
Automatic differentiation with per-variable storage of derivative information (as a tree structure) in a new `Type` [`ADval`](@ref)` <: Real`.
This provides a way of detecting Jacobian (dy/dx) sparsity patterns for problems of form y = f(x)
where x and y are vectors.

 Notes:
  - Will fail if y = f(x) produces a tree with loops, eg by inverting a matrix to solve a linear system.
  - Very slow as each arithmetic operation generates a memory allocation (in order to form the derivative tree structure).
    Speed is approx 10 Mop/s (on a ~4Ghz laptop core) so 100 - 1000x slower than `y = f(x)` with primitive types. Intent
    is that the sparsity information is calculated once and then used to enable a separate fast AD calculator.
  - Also calculates derivative information (slowly).
  - The tree of derivative information grows without limit, so currently limited to cases `y = f(x)` where f does
    not generate too many operations. It would be straightforward to fix this by limiting tree to size `~4*length(x)`
    per element and simplifying it when this limit exceeded (as tree size = 2*length(x) per element is enough for a
    dense Jacobian)
  - Uses DiffRules.jl to autogenerate operations.
  - See ForwardDiff.Dual https://github.com/JuliaDiff/ForwardDiff.jl/blob/master/src/dual.jl
    for examples of definitions required.

# Example usage:
    x_ad = SparsityTracing.create_advec(x)           # x_ad takes values from x,  derivatives set to identity
    y_ad = f(x_ad)                              # tree of derivative information is accumulated into each element of y_ad
    jac = SparsityTracing.jacobian(y_ad, length(x))  # walk trees and generate sparse Jacobian
 """
module SparsityTracing

using SparseArrays
using DiffRules

import SpecialFunctions

"base type for tree of derivative information"
abstract type AbstractDerivNode
end

"""
    ADval

Scalar variable type T with derivatives
"""
struct ADval{T<:Real} <: Real
    val::T
    derivnode::Union{Nothing, AbstractDerivNode}
end

"""
    create_advec(x::Vector{T}) -> Vector{ADval{T, length(x)}}

Create a vector of `ADval` scalars, each with value initialised from `x`
and `length(x)` derivatives initialised to identity (each element derivative 1 wrt itself and 0 elsewhere).
Usual case for an independent vector variable `x` (eg to use as input with `y = f(x)`).

# Examples:
 ```jldoctest; setup = :(import SparsityTracing)
julia> v_ad = SparsityTracing.create_advec([1.0, 5.0])
2-element Vector{SparsityTracing.ADval{Float64}}:
 SparsityTracing.ADval{Float64}(1.0, <derivnode>)
 SparsityTracing.ADval{Float64}(5.0, <derivnode>)

julia> x_ad, y_ad = v_ad;

julia> x_ad        # deriv is sparse vector
SparsityTracing.ADval{Float64}
  val=1.0
  deriv:
  [1]  =  1.0

julia> y_ad       # deriv is sparse vector
SparsityTracing.ADval{Float64}
  val=5.0
  deriv:
  [2]  =  1.0

julia> SparsityTracing.jacobian(v_ad, 2)  # sparse 2x2 identity matrix.
2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 2 stored entries:
 1.0   ⋅
  ⋅   1.0

julia> z_ad = x_ad*y_ad^2    # d(x*y^2)/dx = y^2 = 5^2 = 25,  d(x*y^2)/dy = 2*x*y = 2*1*5 = 10
SparsityTracing.ADval{Float64}
  val=25.0
  deriv:
  [1]  =  25.0
  [2]  =  10.0
```
"""
function create_advec(x::Vector{T}) where T

    n = length(x)
    ADvec = Vector{ADval{T}}(undef, n)

    for i = 1:n
        ADvec[i] = ADval{T}(x[i], DerivLeaf(i))
    end

    return ADvec
end

###############################
# Construction and conversion #
###############################

"create from a non-AD value, set derivatives to zero"
ADval{T}(v) where {T} = ADval{T}(v, nothing)
ADval{T}(v::ADval{T}) where {T} = v

Base.convert(::Type{ADval{T}}, v) where {T} = ADval{T}(v)
Base.convert(::Type{ADval{T}}, v::Number) where {T} = ADval{T}(v)
Base.convert(::Type{ADval{T}}, d::ADval{T}) where {T} = d

Base.promote_rule(::Type{T}, ::Type{ADval{T}} ) where {T} = ADval{T}

"get value. NB: no automatic conversion to Float64 etc"
value(d::ADval) = d.val

#####################
# Generic functions #
#####################

Base.zero(d::ADval{T}) where {T}        = zero(typeof(d))
Base.zero(::Type{ADval{T}}) where {T}   = ADval{T}(zero(T), nothing)

Base.one(d::ADval{T}) where {T}         = one(typeof(d))
Base.one(::Type{ADval{T}}) where {T}    = ADval{T}(one(T), nothing)

Base.sign(d::ADval{T}) where {T} = sign(d.val)

###################
# Pretty printing #
###################

"compact form"
function Base.show(io::IO, d::ADval)
    print(io, typeof(d), "(", d.val, ", <derivnode>)")
end

"multiline form"
function Base.show(io::IO, ::MIME"text/plain", d::ADval)
    println(io, typeof(d))
    println(io, "  val=", d.val)
    println(io, "  deriv:\n", deriv(d))
end

###################
# Get derivatives #
###################

"evaluate derivnode by traversing tree, and return as a sparse vector"
function deriv(d::ADval{T}, N) where {T}
    I = Vector{Int64}()
    V = Vector{T}()
    reduce_deriv!(I, V, 1.0, d.derivnode)
    # reduce_deriv_tree!(I, V, 1.0, d.derivnode)
    return sparsevec(I, V, N)
end

function deriv(d::ADval{T}) where {T}
    I = Vector{Int64}()
    V = Vector{T}()
    reduce_deriv!(I, V, 1.0, d.derivnode)
    # reduce_deriv_tree!(I, V, 1.0, d.derivnode)
    return sparsevec(I, V)
end

"""
    jacobian(advec::Vector{ADval}, N) -> SparseMatrixCSC

Aggregate derivatives to form sparse Jacobian.
"""
function jacobian(advec::Vector{ADval{T}}, N) where {T}

    Nvar = length(advec)

    # @info "SparsityTracing.jacobian ($Nvar, $N)"

    # define arrays for sparse matrix creation
    I = Vector{Int64}()
    J = Vector{Int64}()
    V = Vector{T}()

    for i = 1:Nvar
        # @info "  i=$i"
        derivJ, derivV = findnz(deriv(advec[i], N))
        # @info "  length(derivJ) = $(length(derivJ))"
        for j = 1:length(derivJ)
            push!(I, i)
            push!(J, derivJ[j])
            push!(V, derivV[j])
        end
    end

    return sparse(I, J, V, Nvar, N)
end

"doesn't need to be mutable, but slightly faster (with Julia 1.5)"
mutable struct DerivLeaf <: AbstractDerivNode
    index::Int64
end

"doesn't need to be mutable, but slightly faster (with Julia 1.5)"
mutable struct DerivNode{T} <: AbstractDerivNode
    lmult::T
    lnode::Union{Nothing, AbstractDerivNode}
    rmult::T
    rnode::Union{Nothing, AbstractDerivNode}
end

function reduce_deriv!(I, v, mult, node::Nothing)
    return nothing
end

function reduce_deriv!(I, V, mult, node::DerivLeaf)
    push!(I, node.index)
    push!(V, mult)
    return nothing
end

function reduce_deriv!(I, V, mult, node::DerivNode)
    reduce_deriv!(I, V, mult*node.lmult, node.lnode)
    reduce_deriv!(I, V, mult*node.rmult, node.rnode)
    return nothing
end

function reduce_deriv_tree!(I, V, mult, root::Union{Nothing, DerivNode, DerivLeaf})

    stack = Vector{Tuple{Float64, Union{Nothing, DerivNode, DerivLeaf}}}()

    curr = (mult, root)
    npush = 0

    while curr[2] isa DerivNode || !isempty(stack)
        # reach the left-most node of the curr node
        while curr[2] isa DerivNode
            # put curr on stack and traverse left subtree
            push!(stack, curr)
            mult, node = curr
            curr = (mult*node.lmult, node.lnode)
            npush += 1
            if npush % 1000000 == 0
                @info "    npush $npush length(stack, I, V) $(length(stack)) $(length(I)) $(length(V))"
            end
        end
        # curr[2] must be either Nothing or a DerivLeaf
        if curr[2] isa DerivLeaf
            mult, node = curr
            push!(I, node.index)
            push!(V, mult)
            # @info "    $(node.index)   $mult"
        end
        mult, node = pop!(stack)
        # we have visited the node and its left subtree
        # visit the right subtree
        curr = (mult*node.rmult, node.rnode)
    end

    # curr[2] must be either Nothing or a DerivLeaf
    if curr[2] isa DerivLeaf
        mult, node = curr
        push!(I, node.index)
        push!(V, mult)
        # @info "    $(node.index)   $mult"
    end

    # @info " reduce_deriv_tree! length(I) = $(length(I))"
    return nothing

end

##########################################
# Autogenerate functions using DiffRules #
##########################################

"display DiffRules rules (diagnostic for development)"
function list_diffrules(arity)
    for key in DiffRules.diffrules()
        if key[3] in arity
            println(key)
        end
    end
end

"generate a unary AD rule
   ns       eg :Base
   op       eg :sin"
function gen_unary_rule(ns, op)

    deriv = DiffRules.diffrule(ns, op, :(u.val)) # eg :(cos(u.val))

    @eval $ns.$op(u::ADval{T}) where {T} = ADval{T}($ns.$op(u.val), DerivNode(convert(T, $deriv), u.derivnode, zero(T), nothing))
end

"generate a binary AD rule
   ns       eg :Base
   op       eg :+"
function gen_binary_rule(ns, op)

    u_deriv, v_deriv = DiffRules.diffrule(ns, op, :(u.val), :(v.val))
    @eval $ns.$op(u::ADval{T}, v::ADval{T}) where {T} = ADval{T}($ns.$op(u.val, v.val),
        DerivNode(convert(T, $u_deriv), u.derivnode, convert(T, $v_deriv), v.derivnode))

    u_deriv, dummy = DiffRules.diffrule(ns, op, :(u.val), :v)
    @eval $ns.$op(u::ADval{T}, v::Real) where {T} = ADval{T}($ns.$op(u.val, v),
        DerivNode(convert(T, $u_deriv), u.derivnode, zero(T), nothing))

    dummy, v_deriv = DiffRules.diffrule(ns, op, :u, :(v.val))
    @eval $ns.$op(u::Real, v::ADval{T}) where {T} = ADval{T}($ns.$op(u, v.val),
        DerivNode(convert(T, $v_deriv), v.derivnode, zero(T), nothing))
end

"generate a logical rule (ignore derivatives, just test value)"
function gen_logical_rule(op)
    @eval Base.$op(u::ADval{T}, v::ADval{T}) where {T} = $op(u.val, v.val)
    @eval Base.$op(u::ADval{T}, v::Real) where {T} = $op(u.val, v)
    @eval Base.$op(u::Real, v::ADval{T}) where {T} = $op(u, v.val)
end

function gen_unary_rules(namespaces, excludeops)
    for (ns, op, arity) in DiffRules.diffrules()
        if ns in namespaces && arity == 1 && !(op in excludeops)
            gen_unary_rule(ns, op)
        end
    end
end

function gen_binary_rules(namespaces, excludeops)
    for (ns, op, arity) in DiffRules.diffrules()
        if ns in namespaces && arity == 2 && !(op in excludeops)
            gen_binary_rule(ns, op)
        end
    end
end

function gen_logical_rules()
    for op in (:<, :>, :>=, :<=, :(==))
        gen_logical_rule(op)
    end
end

# add functions to Base
gen_unary_rules((:Base, :SpecialFunctions), ())
gen_binary_rules((:Base, :SpecialFunctions), ())
gen_logical_rules()

###############
# TODO fixups #
###############
# Errors Base.:^(d::ADval{T}, n::Integer) is ambiguous ?
# adding literal_pow apparently fixes this but might not be the right fix
Base.literal_pow(::typeof(^), d::ADval{T}, ::Val{0}) where {T} = one(d)
Base.literal_pow(::typeof(^), d::ADval{T}, ::Val{1}) where {T} = d
Base.literal_pow(::typeof(^), d::ADval{T}, ::Val{2}) where {T} = d*d
Base.literal_pow(::typeof(^), d::ADval{T}, ::Val{3}) where {T} = d*d*d
Base.literal_pow(::typeof(^), d::ADval{T}, ::Val{4}) where {T} = d*d*d*d

end


text:
# Test data container processing functionality
function t_libdata()

#############################################################
# Test data generator (test that datasets can be generated) #
#############################################################
Test.@test try j4pr.DataGenerator.iris();true catch; false end
Test.@test try j4pr.DataGenerator.quakes();true catch; false end
Test.@test try j4pr.DataGenerator.boston();true catch; false end
Test.@test try j4pr.DataGenerator.normal2d();true catch; false end
Test.@test try j4pr.DataGenerator.fish(10);true catch; false end
Test.@test try j4pr.DataGenerator.fishr();true catch; false end
Test.@test try j4pr.DataGenerator.fnoise(rand(10),x->2x);true catch; false end

###############################################################
# Test images interface (test that datasets can be generated) #
###############################################################
A = [1.0 1.1; 1.2 1.3]
img = Images.Gray.(A)
Test.@test j4pr.im2targets(A) == ([1.0 2.0 1.0 2.0; 1.0 1.0 2.0 2.0],
			    	  Images.ColorTypes.Gray{Float64}[Images.Gray{Float64}(1.0), Images.Gray{Float64}(1.2),
							   Images.Gray{Float64}(1.1), Images.Gray{Float64}(1.3)])
Test.@test j4pr.targets2im(j4pr.datacell(j4pr.im2targets(Images.Gray.(A)))) == A

########################
# Test label utilities #
########################
C = 2
N = 10
A = rand(10)
l = rand(Bool,N)
l2 = rand(Bool,C,N)

# addlabels
for Di in [j4pr.datacell(A), j4pr.datacell(A,l), j4pr.datacell(A,l2)]
	for li in [l,l2]
		Di2 = j4pr.addlabels(Di,li)
		Test.@test j4pr.getx!(Di2) == j4pr.getx!(Di)
		Test.@test j4pr.gety!(Di2) == j4pr.dcat(j4pr.gety!(Di),li)
		Test.@test j4pr.gety!(j4pr.unlabel(Di2)) == nothing
	end
end

# labelize
D = j4pr.datacell(A,l2)

wl = j4pr.labelize(l, true)
wl2 = j4pr.labelize(l, false)
Test.@test j4pr.gety!(D |> wl) == l
Test.@test j4pr.gety!(D |> wl2) == j4pr.dcat(l2,l)

wl = j4pr.labelize(x->x[1], true)
wl2 = j4pr.labelize(x->x[1], false)
Test.@test j4pr.gety!(D |> wl) == l2[1,:]
Test.@test j4pr.gety!(D |> wl2) == j4pr.dcat(l2,l2[1,:])
Test.@test try
	j4pr.unlabel(D) |> wl
	false
catch
	true # has to fail
end

A=[1 2 3; 4 5 6; 7 8 9]
wl = j4pr.labelize(1, true)
wl2 = j4pr.labelize(1, false)
wl3 = j4pr.labelize(indmax, 1:2, true)
wl4 = j4pr.labelize(indmax, 1:2, false)
D = j4pr.datacell(A)
Test.@test j4pr.getx!(D |> wl) == A[2:3,:]
Test.@test j4pr.gety!(D |> wl) == A[1,:]
Test.@test j4pr.getx!(D |> wl2) == A
Test.@test j4pr.gety!(D |> wl2) == A[1,:]
Test.@test j4pr.getx!(D |> wl3) == A[3:3,:]
Test.@test j4pr.gety!(D |> wl3) == 2*ones(3)
Test.@test j4pr.getx!(D |> wl4) == A
Test.@test j4pr.gety!(D |> wl4) == 2*ones(3)

#####################################
# Test class slicing, domain filter #
#####################################

A = [1 2 3;
     4 5 6]
l = [1 1 2]
l2 = [0 1 0;
      1 0 1]

Du = j4pr.datacell(A)
Dl = j4pr.datacell(A,l)
Dl2 = j4pr.datacell(A,l2)

# class slicing
Test.@test j4pr.cslice(Du,[1,]) isa j4pr.DataCell{<:SubArray}
Test.@test j4pr.getobs(j4pr.getx!(j4pr.cslice(Dl,[2]))) == Matrix([3;6])
Test.@test j4pr.getobs(j4pr.getx!(j4pr.cslice(Dl2,[1],1))) == Matrix([2;5])

# domain filter
A = [1.0 2 3;
     4 5 6;
     7 8 9;
     10 11 12;
     13 14 15]

l = [1 1 2]
l2 = [0 1 0;
      1 0 1]

Du = j4pr.datacell(A)
Dl = j4pr.datacell(A,l)
Dl2 = j4pr.datacell(A,l2)
(+Dl2,-Dl2) |> j4pr.filterdomain!(Dict(1=>x->x>2 ? NaN : x)) # the '3' from the first row becomes NaN
Dl2 |> j4pr.filterdomain!(Dict(2:3=>[4 9])) # the '9' from the third row becomes NaN
Dl2 |> j4pr.filterdomain!(Dict([4]=>[11 13; 14 16])) # the '10' from the fourth row becomes NaN
Dl2 |> j4pr.filterdomain!(Dict([5]=>15)) # the '13' and '14' from the fifth row become NaN

Dl2f = j4pr.getx!(Dl2)
v = [
       1.0    2.0    NaN;
       4.0    5.0    6.0;
       7.0    8.0    NaN;
       NaN     11.0   12.0;
       NaN    NaN     15.0
]
Test.@test isequal(Dl2f,v)

########################
# Test one hot encoder #
########################

A = [1.0 2 3;
     4 5 6;
    ]

# binary encoder
w = A|> j4pr.ohenc("binary")
T = [1. 6; 2 6]
Test.@test isequal(T|>w, [1.0  NaN;
			0.0  NaN;
   			0.0  NaN;
 			NaN  0.0;
 			NaN  0.0;
 			NaN  1.0])

# integer encoder
w = A|> j4pr.ohenc("integer")
T = [1. 6; 2 6]
Test.@test isequal(T|>w, [1.0  NaN; NaN 3.0])

# mixed case
A = [1.0 2 3;
     4 5 6;
     7 8 9]
w = A|> j4pr.ohenc(Dict(1=>"binary",2=>"integer"))
T = [1. 6; 2 6; -1 -2]
Test.@test isequal(T|>w, [1.0  NaN; # binary
			0.0  NaN;   # binary
   			0.0  NaN;   # binary
 			NaN  3.0;   # integer
 			-1   -2])   # as input

#################
# Test sampling #
#################
A = j4pr.datacell([1 2 3; 4 5 6],["a","a","b"])
(As, idxs) =  A |> j4pr.sample(7)
Test.@test j4pr.nobs(As) == 7
Test.@test j4pr.getx!(A)[:,idxs] == j4pr.getobs(j4pr.getx!(As))

(As, idxs) =  A |> j4pr.sample(0.5)
Test.@test j4pr.nobs(As) == 2
Test.@test j4pr.getx!(A)[:,idxs] == j4pr.getobs(j4pr.getx!(As))

(As, idxs) =  A |> j4pr.sample(Dict("a"=>10, "b"=>2))
Test.@test j4pr.nobs(As) == 12
Test.@test j4pr.getx!(A)[:,idxs] == j4pr.getobs(j4pr.getx!(As))
y = j4pr.getobs(j4pr.gety!(As))
Test.@test sum(y.=="a") == 10 && sum(y.=="b") == 2

#######################################################################
# Test generic filter (verify that works, functionality is not tested #
#######################################################################

for fo in ["mean", "c-mean", "median", "c-median", "majority", "c-majority", (x,y)->1.0, Dict(1=>"mean",2=>(x,y)->x)]
	w = j4pr.filterg(fo)
	Au = j4pr.datacell(rand(2,10))
	Al = j4pr.datacell(rand(2,10), rand([1,2],10))

	Test.@test try
       		Au |> (Au |> w)
     		Al |> (Al |> w)
		true
	catch
		false
	end
end

###############################################################
# Test scaler (verify that works, functionality is not tested #
###############################################################

for so in ["mean", "c-mean", "variance", "c-variance", "domain", Dict(1=>"mean", 2=>"variance")]
	w = j4pr.scaler!(so)
	Au = j4pr.datacell(rand(2,10))
	Al = j4pr.datacell(rand(2,10), rand([1,2],10))

	Test.@test try
       		Au |> (Au |> w)
     		Al |> (Al |> w)
		true
	catch
		false
	end
end

# Separate case for 2-sigma scaling, has to fail for unlabeled containers
w = j4pr.scaler!("2-sigma")
Au = j4pr.datacell(rand(2,10))
Al = j4pr.datacell(rand(2,10), rand([1,2],10))

Test.@test try
	Au |> (Au |>w)
	false
catch
	true
end

Test.@test try
	Al |> (Al |>w)
	true
catch
	false
end

end


text:
@deprecate pmap(fn::Function, itr) tmap(fn::Function, itr)
@deprecate pforeach(fn::Function, itr) tforeach(fn::Function, itr)
@deprecate logpmap(fn::Function, itr) logtmap(fn::Function, itr)
@deprecate logpforeach(fn::Function, itr) logtforeach(fn::Function, itr)

"""
    tmap(fn::Function, itr) -> collection

Mimics `Base.map`, but launches the function evaluations onto all available
threads, using a pre-assigned scheduling strategy appropriate for uniform
task durations.

# Example
```julia
julia> tmap(x -> begin; println((x,Threads.threadid())); x^2; end, 1:8)'
(7, 4)
(5, 3)
(8, 4)
(1, 1)
(6, 3)
(2, 1)
(3, 2)
(4, 2)
1×8 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
 1  4  9  16  25  36  49  64
```
Note that while the execution order is not guaranteed, the result order is.
Also note that the primary thread is used.
"""
function tmap(fn::Function, itr)
    result = Vector{_detect_type(fn, itr)}()
    pool = StaticPool()
    result = tmap(pool, fn, itr)
    close(pool)
    return result
end

"""
    bmap(fn::Function, itr) -> collection

Mimics `Base.map`, but launches the function evaluations onto all available
threads except the primary, using a pre-assigned scheduling strategy
appropriate for uniform task durations.

# Example
```julia
julia> bmap(x -> begin; println((x,Threads.threadid())); x^2; end, 1:8)'
(6, 4)
(1, 2)
(3, 3)
(2, 2)
(4, 3)
(7, 4)
(5, 3)
(8, 4)
1×8 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
 1  4  9  16  25  36  49  64
```
Note that while the execution order is not guaranteed, the result order is,
Also note that the primary thread is not used.
"""
function bmap(fn, itr)
    result = Vector{_detect_type(fn, itr)}()
    pool = StaticPool(2)
    result = tmap(pool, fn, itr)
    close(pool)
    return result
end

"""
    qmap(fn::Function, itr) -> collection

Mimics `Base.map`, but launches the function evaluations onto all available
threads, using a queued scheduling strategy appropriate for nonuniform
task durations.

# Example
```julia
julia> qmap(x -> begin; println((x,Threads.threadid())); x^2; end, 1:8)'
(2, 3)
(3, 2)
(4, 4)
(5, 3)
(6, 2)
(7, 4)
(8, 3)
(1, 1)
1×8 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
 1  4  9  16  25  36  49  64
```
Note that while the execution order is not guaranteed, the result order is.
Also note that the primary thread is used.
"""
function qmap(fn, itr)
    result = Vector{_detect_type(fn, itr)}()
    pool = QueuePool()
    result = tmap(pool, fn, itr)
    close(pool)
    return result
end

"""
    qbmap(fn::Function, itr) -> collection

Mimics `Base.map`, but launches the function evaluations onto all available
threads except the primary, using a queued scheduling strategy appropriate
for nonuniform task durations.

# Example
```julia
julia> qbmap(x -> begin; println((x,Threads.threadid())); x^2; end, 1:8)'
(2, 3)
(1, 2)
(3, 4)
(5, 2)
(4, 3)
(6, 4)
(7, 2)
(8, 3)
1×8 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
 1  4  9  16  25  36  49  64
```
Note that while the execution order is not guaranteed, the result order is,
Also note that the primary thread is not used.
"""
function qbmap(fn, itr)
    result = Vector{_detect_type(fn, itr)}()
    pool = QueuePool(2)
    result = tmap(pool, fn, itr)
    close(pool)
    return result
end

"""
    logtmap(fn::Function, itr) -> (pool, collection)

Mimics `Base.map`, but launches the function evaluations onto all available
threads, using a pre-assigned scheduling strategy appropriate for uniform
task durations.  Also returns a logged pool that can be analyzed with
the logging functions and `plot`ted.

# Example
```julia
julia> (pool, result) = logtmap(1:8) do x
         println((x,Threads.threadid()))
         x^2
       end;
(1, 1)
(3, 2)
(7, 4)
(5, 3)
(8, 4)
(4, 2)
(2, 1)
(6, 3)

julia> result'
1×8 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
1  4  9  16  25  36  49  64

julia> plot(pool)
```
Note that while the execution order is not guaranteed, the result order is.
Also note that the primary thread is used.
"""
function logtmap(fn::Function, itr)
    result = Vector{_detect_type(fn, itr)}()
    pool = LoggedStaticPool()
    result = tmap(pool, fn, itr)
    close(pool)
    return pool, result
end

"""
    logbmap(fn::Function, itr) -> (pool, collection)

Mimics `Base.map`, but launches the function evaluations onto all available
threads except the primary, using a pre-assigned scheduling strategy
appropriate for uniform task durations.  Also returns a logged pool that can
be analyzed with the logging functions and `plot`ted.

# Example
```julia
julia> (pool, result) = logbmap(1:8) do x
         println((x,Threads.threadid()))
         x^2
       end;
(1, 2)
(6, 4)
(3, 3)
(7, 4)
(2, 2)
(4, 3)
(8, 4)
(5, 3)

julia> result'
1×8 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
1  4  9  16  25  36  49  64

julia> plot(pool)
```
Note that while the execution order is not guaranteed, the result order is,
Also note that the primary thread is not used.
"""
function logbmap(fn, itr)
    result = Vector{_detect_type(fn, itr)}()
    pool = LoggedStaticPool(2)
    result = tmap(pool, fn, itr)
    close(pool)
    return pool, result
end

"""
    logqmap(fn::Function, itr) -> (pool, collection)

Mimics `Base.map`, but launches the function evaluations onto all available
threads, using a queued scheduling strategy appropriate for nonuniform
task durations.  Also returns a logged pool that can be analyzed with the
logging functions and `plot`ted.

# Example
```julia
julia> (pool, result) = logqmap(1:8) do x
         println((x,Threads.threadid()))
         x^2
        end;
(3, 3)
(4, 4)
(2, 2)
(5, 3)
(7, 2)
(6, 4)
(8, 3)
(1, 1)

julia> result'
1×8 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
1  4  9  16  25  36  49  64

julia> plot(pool)
```
Note that while the execution order is not guaranteed, the result order is.
Also note that the primary thread is used.
"""
function logqmap(fn, itr)
    result = Vector{_detect_type(fn, itr)}()
    pool = LoggedQueuePool()
    result = tmap(pool, fn, itr)
    close(pool)
    return pool, result
end

"""
    logqbmap(fn::Function, itr) -> (pool, collection)

Mimics `Base.map`, but launches the function evaluations onto all available
threads except the primary, using a queued scheduling strategy appropriate
for nonuniform task durations.  Also returns a logged pool that can be
analyzed with the logging functions and `plot`ted.

# Example
```julia
julia> (pool, result) = logqbmap(1:8) do x
         println((x,Threads.threadid()))
         x^2
       end;
(3, 3)
(2, 4)
(1, 2)
(4, 3)
(5, 4)
(6, 2)
(7, 3)
(8, 4)

julia> result'
1×8 LinearAlgebra.Adjoint{Int64,Array{Int64,1}}:
1  4  9  16  25  36  49  64

julia> plot(pool)
```
Note that while the execution order is not guaranteed, the result order is,
Also note that the primary thread is not used.
"""
function logqbmap(fn, itr)
    result = Vector{_detect_type(fn, itr)}()
    pool = LoggedQueuePool(2)
    result = tmap(pool, fn, itr)
    close(pool)
    return pool, result
end

"""
    tforeach(fn::Function, itr)

Mimics `Base.foreach`, but launches the function evaluations onto all available
threads, using a pre-assigned scheduling strategy appropriate for uniform
task durations.

# Example
```julia
julia> tforeach(x -> println((x,Threads.threadid())), 1:8)
(1, 1)
(3, 2)
(5, 3)
(2, 1)
(7, 4)
(4, 2)
(8, 4)
(6, 3)
```
Note that the execution order is not guaranteed, and that the primary thread
is used.
"""
function tforeach(fn::Function, itr)
    pool = StaticPool()
    tforeach(pool, fn, itr)
    close(pool)
    nothing
end

"""
    bforeach(fn::Function, itr)

Mimics `Base.foreach`, but launches the function evaluations onto all available
threads except the primary, using a pre-assigned scheduling strategy appropriate
for uniform task durations.

# Example
```julia
julia> bforeach(x -> println((x,Threads.threadid())), 1:8)
(1, 2)
(6, 4)
(2, 2)
(7, 4)
(8, 4)
(3, 3)
(4, 3)
(5, 3)
```
Note that the execution order is not guaranteed, and that the primary thread
is not used.
"""
function bforeach(fn, itr)
    pool = StaticPool(2)
    tforeach(pool, fn, itr)
    close(pool)
    nothing
end

"""
    qforeach(fn::Function, itr)

Mimics `Base.foreach`, but launches the function evaluations onto all available
threads, using a queued scheduling strategy appropriate for nonuniform
task durations.

# Example
```julia
julia> qforeach(x -> println((x,Threads.threadid())), 1:8)
(4, 3)
(2, 2)
(3, 4)
(5, 3)
(6, 2)
(7, 4)
(8, 3)
(1, 1)
```
Note that the execution order is not guaranteed, and that the primary thread
is used.
"""
function qforeach(fn, itr)
    pool = QueuePool()
    tforeach(pool, fn, itr)
    close(pool)
    nothing
end

"""
    qbforeach(fn::Function, itr)

Mimics `Base.foreach`, but launches the function evaluations onto all available
threads except the primary, using a queued scheduling strategy appropriate for
nonuniform task durations.

# Example
```julia
julia> qbforeach(x -> println((x,Threads.threadid())), 1:8)
(3, 3)
(2, 4)
(1, 2)
(4, 3)
(5, 4)
(6, 2)
(7, 3)
(8, 4)
```
Note that the execution order is not guaranteed, and that the primary thread
is not used.
"""
function qbforeach(fn, itr)
    pool = QueuePool(2)
    tforeach(pool, fn, itr)
    close(pool)
    nothing
end

"""
    logtforeach(fn::Function, itr) -> pool

Mimics `Base.foreach`, but launches the function evaluations onto all available
threads, using a pre-assigned scheduling strategy appropriate for uniform
task durations.  Returns a logged pool that can be analyzed with
the logging functions and `plot`ted.

# Example
```julia
julia> pool = logtforeach(x -> println((x,Threads.threadid())), 1:8);
(1, 1)
(3, 2)
(7, 4)
(2, 1)
(4, 2)
(5, 3)
(8, 4)
(6, 3)

julia> plot(pool)
```
Note that the execution order is not guaranteed, and that the primary thread
is used.
"""
function logtforeach(fn::Function, itr)
    pool = LoggedStaticPool()
    tforeach(pool, fn, itr)
    close(pool)
    return pool
end

"""
    logbforeach(fn::Function, itr)

Mimics `Base.foreach`, but launches the function evaluations onto all available
threads except the primary, using a pre-assigned scheduling strategy appropriate
for uniform task durations.  Returns a logged pool that can be analyzed with
the logging functions and `plot`ted.

# Example
```julia
julia> pool = logbforeach(x -> println((x,Threads.threadid())), 1:8);
(1, 2)
(3, 3)
(6, 4)
(4, 3)
(2, 2)
(7, 4)
(5, 3)
(8, 4)

julia> plot(pool)
```
Note that the execution order is not guaranteed, and that the primary thread
is not used.
"""
function logbforeach(fn, itr)
    pool = LoggedStaticPool(2)
    tforeach(pool, fn, itr)
    close(pool)
    return pool
end

"""
    logqforeach(fn::Function, itr)

Mimics `Base.foreach`, but launches the function evaluations onto all available
threads, using a queued scheduling strategy appropriate for nonuniform
task durations.

# Example
```julia
julia> pool = logqforeach(x -> println((x,Threads.threadid())), 1:8);
(2, 4)
(3, 3)
(4, 2)
(5, 4)
(6, 3)
(7, 2)
(8, 4)
(1, 1)

julia> plot(pool)
```
Note that the execution order is not guaranteed, and that the primary thread
is used.  Returns a logged pool that can be analyzed with the logging functions
and `plot`ted.
"""
function logqforeach(fn, itr)
    pool = LoggedQueuePool()
    tforeach(pool, fn, itr)
    close(pool)
    return pool
end

"""
    logqbforeach(fn::Function, itr)

Mimics `Base.foreach`, but launches the function evaluations onto all available
threads except the primary, using a queued scheduling strategy appropriate for
nonuniform task durations.  Returns a logged pool that can be analyzed with the
logging functions and `plot`ted.

# Example
```julia
julia> pool = logqbforeach(x -> println((x,Threads.threadid())), 1:8);
(2, 2)
(1, 3)
(3, 4)
(4, 2)
(5, 3)
(6, 4)
(7, 2)
(8, 3)

julia> plot(pool)
```
Note that the execution order is not guaranteed, and that the primary thread
is not used.
"""
function logqbforeach(fn, itr)
    pool = LoggedQueuePool(2)
    tforeach(pool, fn, itr)
    close(pool)
    return pool
end


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

### Route Planning for OpenStreetMap ###

### Get list of vertices (highway nodes) in specified levels of classes ###
# For all highways
function highwayVertices(highways::Dict{Int,Highway})
    vertices = Set{Int}()

    for highway in values(highways)
        union!(vertices, highway.nodes)
    end

    return vertices
end

# For classified highways
function highwayVertices(highways::Dict{Int,Highway}, classes::Dict{Int,Int})
    vertices = Set{Int}()

    for key in keys(classes)
        union!(vertices, highways[key].nodes)
    end

    return vertices
end

# For specified levels of a classifier dictionary
function highwayVertices(highways::Dict{Int,Highway}, classes::Dict{Int,Int}, levels)
    vertices = Set{Int}()

    for (key, class) in classes
        if in(class, levels)
            union!(vertices, highways[key].nodes)
        end
    end

    return vertices
end

### Form transportation network graph of map ###
function createGraph(nodes, highways::Dict{Int,Highway}, classes, levels, reverse::Bool=false)
    v = Dict{Int,Graphs.KeyVertex{Int}}()                       # Vertices
    w = Float64[]                                               # Weights
    g_classes = Int[]                                           # Road classes
    g = Graphs.inclist(Graphs.KeyVertex{Int}, is_directed=true) # Graph

    verts = highwayVertices(highways, classes, levels)
    for vert in verts
        v[vert] = Graphs.add_vertex!(g, vert)
    end

    for (key, class) in classes
        if in(class, levels)
            highway = highways[key]
            if length(highway.nodes) > 1
                # Add edges to graph and compute weights
                for n = 2:length(highway.nodes)
                    if reverse
                        node0 = highway.nodes[n]
                        node1 = highway.nodes[n-1]
                    else
                        node0 = highway.nodes[n-1]
                        node1 = highway.nodes[n]
                    end
                    edge = Graphs.make_edge(g, v[node0], v[node1])
                    Graphs.add_edge!(g, edge)
                    weight = distance(nodes, node0, node1)
                    push!(w, weight)
                    push!(g_classes, class)
                    # node_set = Set(node0, node1)

                    if !highway.oneway
                        edge = Graphs.make_edge(g, v[node1], v[node0])
                        Graphs.add_edge!(g, edge)
                        push!(w, weight)
                        push!(g_classes, class)
                    end
                end
            end
        end
    end

    return Network(g, v, w, g_classes)
end

### Form transportation network graph of map ###
function createGraph(segments::Vector{Segment}, intersections, reverse::Bool=false)
    v = Dict{Int,Graphs.KeyVertex{Int}}()                       # Vertices
    w = Float64[]                                               # Weights
    class = Int[]                                               # Road class
    g = Graphs.inclist(Graphs.KeyVertex{Int}, is_directed=true) # Graph

    for vert in keys(intersections)
        v[vert] = Graphs.add_vertex!(g, vert)
    end

    for segment in segments
        # Add edges to graph and compute weights
        if reverse
            node0 = segment.node1
            node1 = segment.node0
        else
            node0 = segment.node0
            node1 = segment.node1
        end
        edge = Graphs.make_edge(g, v[node0], v[node1])
        Graphs.add_edge!(g, edge)
        weight = segment.dist
        push!(w, weight)
        push!(class, segment.class)
        # node_set = Set(node0, node1)

        if !segment.oneway
            edge = Graphs.make_edge(g, v[node1], v[node0])
            Graphs.add_edge!(g, edge)
            push!(w, weight)
            push!(class, segment.class)
        end
    end

    return Network(g, v, w, class)
end

# Put all edges in network.g in an array, indexed by their edge index
function getEdges( network::Network )
    edges = Array(Any,Graphs.num_edges(network.g))
    vertices = Graphs.vertices(network.g)

    for v in vertices
        out_edges = Graphs.out_edges(v,network.g)
        for edge in out_edges
            edges[edge.index] = edge
        end
    end

    return edges
end

### Get distance between two nodes ###
# ENU Coordinates
function Geodesy.distance{T<:@compat(Union{ENU,ECEF})}(nodes::Dict{Int,T}, node0, node1)
    loc0 = nodes[node0]
    loc1 = nodes[node1]

    return distance(loc0, loc1)
end

### Compute the distance of a route ###
function Geodesy.distance{T<:@compat(Union{ENU,ECEF})}(nodes::Dict{Int,T}, route::Vector{Int})
    if length(route) == 0
        return Inf
    end

    dist = 0.0
    prev_point = nodes[route[1]]
    for i = 2:length(route)
        point = nodes[route[i]]
        dist += distance(prev_point, point)
        prev_point = point
    end

    return dist
end

### Shortest Paths ###
# Dijkstra's Algorithm
function dijkstra(g, w, start_vertex)
    return Graphs.dijkstra_shortest_paths(g, w, start_vertex)
end

# Bellman Ford's Algorithm
function bellmanFord(g, w, start_vertices)
    return Graphs.bellman_ford_shortest_paths(g, w, start_vertices)
end

# Extract route from Dijkstra results object
function extractRoute(dijkstra, start_index, finish_index)
    route = Int[]

    distance = dijkstra.dists[finish_index]

    if distance != Inf
        index = finish_index
        push!(route, index)
        while index != start_index
            index = dijkstra.parents[index].index
            push!(route, index)
        end
    end

    reverse!(route)

    return route, distance
end

### Generate an ordered list of edges traversed in route
function routeEdges(network::Network, route::Vector{Int})
    e = Array(Int, length(route)-1)

    # For each node pair, find matching edge
    for n = 1:length(route)-1
        s = route[n]
        t = route[n+1]

        for e_candidate in Graphs.out_edges(network.v[s],network.g)
            if t == e_candidate.target.key
                e[n] = e_candidate.index
                break
            end
        end
    end

    return e
end

### Shortest Route ###
function shortestRoute(network, node0, node1)
    start_vertex = network.v[node0]

    dijkstra_result = dijkstra(network.g, network.w, start_vertex)

    start_index = network.v[node0].index
    finish_index = network.v[node1].index
    route_indices, distance = extractRoute(dijkstra_result, start_index, finish_index)

    route_nodes = getRouteNodes(network, route_indices)

    return route_nodes, distance
end

function getRouteNodes(network, route_indices)
    route_nodes = Array(Int, length(route_indices))
    v = Graphs.vertices(network.g)
    for n = 1:length(route_indices)
        route_nodes[n] = v[route_indices[n]].key
    end

    return route_nodes
end

function networkTravelTimes(network, class_speeds)
    w = Array(Float64, length(network.w))
    for k = 1:length(w)
        w[k] = network.w[k] / class_speeds[network.class[k]]
        w[k] *= 3.6 # (3600/1000) unit conversion to seconds
    end
    return w
end

### Fastest Route ###
function fastestRoute(network, node0, node1, class_speeds=SPEED_ROADS_URBAN)
    start_vertex = network.v[node0]

    # Modify weights to be times rather than distances
    w = networkTravelTimes(network, class_speeds)
    dijkstra_result = dijkstra(network.g, w, start_vertex)

    start_index = network.v[node0].index
    finish_index = network.v[node1].index
    route_indices, route_time = extractRoute(dijkstra_result, start_index, finish_index)

    route_nodes = getRouteNodes(network, route_indices)

    return route_nodes, route_time
end

function filterVertices(vertices, weights, limit)
    if limit == Inf
        @assert length(vertices) == length(weights)
        return keys(vertices), weights
    end
    indices = Int[]
    distances = Float64[]
    for vertex in vertices
        distance = weights[vertex.index]
        if distance < limit
            push!(indices, vertex.key)
            push!(distances, distance)
        end
    end
    return indices, distances
end

# Extract nodes from BellmanFordStates object within an (optional) limit
# based on driving distance
function nodesWithinDrivingDistance(network::Network, start_indices, limit=Inf)
    start_vertices = [network.v[i] for i in start_indices]
    bellmanford = bellmanFord(network.g, network.w, start_vertices)
    return filterVertices(values(network.v), bellmanford.dists, limit)
end

function nodesWithinDrivingDistance(network::Network,
                                    loc::ENU,
                                    limit=Inf,
                                    loc_range=100.0)
    return nodesWithinDrivingDistance(network,
                                      nodesWithinRange(network.v, loc, loc_range),
                                      limit)
end

# Extract nodes from BellmanFordStates object within a (optional) limit,
# based on driving time
function nodesWithinDrivingTime(network,
                                start_indices,
                                limit=Inf,
                                class_speeds=SPEED_ROADS_URBAN)
    # Modify weights to be times rather than distances
    w = networkTravelTimes(network, class_speeds)
    start_vertices = [network.v[i] for i in start_indices]
    bellmanford = bellmanFord(network.g, w, start_vertices)
    return filterVertices(values(network.v), bellmanford.dists, limit)
end

function nodesWithinDrivingTime(network::Network,
                                loc::ENU,
                                limit=Inf,
                                class_speeds=SPEED_ROADS_URBAN,
                                loc_range=100.0)
    return nodesWithinDrivingTime(network,
                                  nodesWithinRange(network.v, loc, loc_range),
                                  limit)
end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export hamming_distance, setbit, resetbit

# we export a few bit-handling primitives
# BitArray contains too much extra functionality at the cost of O(1) extra words,
# however this could be an issue, since we represent our database as n vectors

"""
    hamming_distance(a, b)::Float64
    hamming_distance(a::AbstractVector, b::AbstractVector)::Float64 where T<:Unsigned

Computes the binary hamming distance for bit types and arrays of bit types
"""
function hamming_distance(a::T, b::T)::Float64 where {T<:Unsigned}
    count_ones(a ⊻ b)
end

function hamming_distance(a::AbstractVector{T}, b::AbstractVector{T})::Float64 where {T<:Unsigned}
    d = 0
    @inbounds @simd for i in eachindex(a)
        d += count_ones(a[i] ⊻ b[i])
    end

    d
end

"""
   setbit(a::Unsigned, i::Int)

Enables the i-th bit in `a`
"""
function setbit(a::Unsigned, i::Int)
    a | (one(a) << i)
end

"""
    resetbit(a::Unsigned, i::Int)

Disables (0 value) the i-th bit of `a`
"""
function resetbit(a::Unsigned, i::Int)
    a & ~(one(a) << i)
end


text:
export map_address, AddressInfo

struct AddressInfo
    nbits::Int
    addresses::Vector{Int}
end
AddressInfo(nbits::Int, ::AllLocs) = AddressInfo(nbits, collect(1:nbits))
AddressInfo(nbits::Int, iter) = AddressInfo(nbits, collect(iter))
Base.copy(info::AddressInfo) = AddressInfo(copy(info.addresses))
Base.:/(locs, info::AddressInfo) = map(loc -> info.addresses[loc], locs)
Base.:/(locs::AllLocs, info::AddressInfo) = info.addresses

"""
    map_address(block::AbstractBlock, info::AddressInfo) -> AbstractBlock

map the locations in `block` to target locations.

# Example

`map_address` can be used to embed a sub-circuit to a larger one.

```jldoctest; setup=:(using YaoBlocks)
julia> c = chain(5, repeat(H, 1:5), put(2=>X), kron(1=>X, 3=>Y))
nqubits: 5
chain
├─ repeat on (1, 2, 3, 4, 5)
│  └─ H
├─ put on (2)
│  └─ X
└─ kron
   ├─ 1=>X
   └─ 3=>Y

julia> map_address(c, AddressInfo(10, [6,7,8,9,10]))
nqubits: 10
chain
├─ repeat on (6, 7, 8, 9, 10)
│  └─ H
├─ put on (7)
│  └─ X
└─ kron
   ├─ 6=>X
   └─ 8=>Y
```
"""
function map_address end

function map_address(block::AbstractBlock, info::AddressInfo)
    throw(NotImplementedError(:map_address, typeof(block)))
end

function map_address(blk::Measure, info::AddressInfo)
    m = Measure{info.nbits}(blk.rng, blk.operator, (blk.locations / info...,), blk.postprocess)
    if isdefined(blk, :results)
        m.results = blk.results
    end
    return m
end

map_address(blk::PrimitiveBlock, info::AddressInfo) = blk
map_address(blk::PutBlock, info::AddressInfo) = put(info.nbits, blk.locs / info => content(blk))

function map_address(blk::ControlBlock, info::AddressInfo)
    ControlBlock{info.nbits}(blk.ctrl_locs / info, blk.ctrl_config, content(blk), blk.locs / info)
end

function map_address(blk::KronBlock, info::AddressInfo)
    kron(info.nbits, [l => G for (l, G) in zip(blk.locs / info, blk.blocks)]...)
end

function map_address(blk::RepeatedBlock, info::AddressInfo)
    repeat(info.nbits, content(blk), blk.locs / info)
end

function map_address(blk::Subroutine, info::AddressInfo)
    subroutine(info.nbits, content(blk), blk.locs / info)
end

function map_address(blk::ChainBlock, info::AddressInfo)
    chain(info.nbits, map(b -> map_address(b, info), subblocks(blk)))
end

function map_address(blk::Daggered, info::AddressInfo)
    Daggered(map_address(content(blk), info))
end

function map_address(blk::CachedBlock, info::AddressInfo)
    CachedBlock(blk.server, map_address(content(blk), info), blk.level)
end

function map_address(blk::Scale, info::AddressInfo)
    Scale(blk.alpha, map_address(content(blk), info))
end

function map_address(blk::Add, info::AddressInfo)
    Add{info.nbits}(map(b -> map_address(b, info), subblocks(blk)))
end


text:
# Called in the OrdinaryDiffEQ.__init; All `OrdinaryDiffEqAlgorithm`s have one
function OrdinaryDiffEq.initialize!(integ, cache::GaussianODEFilterCache)
    if integ.opts.dense && !integ.alg.smooth
        error("To use `dense=true` you need to set `smooth=true`!")
    elseif !integ.opts.dense && integ.alg.smooth
        @warn "If you set dense=false for efficiency, you might also want to set smooth=false."
    end
    if !integ.opts.save_everystep && integ.alg.smooth
        error("If you do not save all values, you do not need to smooth!")
    end
    @assert integ.saveiter == 1

    integ.kshortsize = 1
    resize!(integ.k, integ.kshortsize)
    integ.k[1] = integ.u

    # Update the initial state to the known (given or computed with AD) initial values
    initial_update!(integ, cache)

    # These are necessary since the solution object is not 100% initialized by default
    OrdinaryDiffEq.copyat_or_push!(integ.sol.x_filt, integ.saveiter, cache.x)
    OrdinaryDiffEq.copyat_or_push!(
        integ.sol.pu,
        integ.saveiter,
        mul!(cache.pu_tmp, cache.SolProj, cache.x),
    )
    return nothing
end

"""
    perform_step!(integ, cache::GaussianODEFilterCache, repeat_step=false)

Performs an ODE solver step.
Not necessarily successful! For that, see `step!(integ)`.

Basically consists of the following steps
- Compute the current transition and diffusion matrices
- Predict mean
- Evaluate the ODE (and Jacobian) at the predicted mean; Build measurement mean `z`
- Compute local diffusion and local error estimate
- If the step is rejected, terminate here; Else continue
- Predict the covariance and build the measurement covariance `S`
- Kalman update step
- (optional) Update the global diffusion MLE
"""
function OrdinaryDiffEq.perform_step!(
    integ,
    cache::GaussianODEFilterCache,
    repeat_step=false,
)
    @unpack t, dt = integ
    @unpack d, SolProj = integ.cache
    @unpack x, x_pred, u_pred, x_filt, u_filt, err_tmp = integ.cache
    @unpack x_tmp, x_tmp2 = integ.cache
    @unpack A, Q, Ah, Qh = integ.cache

    make_preconditioners!(cache, dt)
    @unpack P, PI = integ.cache

    tnew = t + dt

    # Build the correct matrices
    @. Ah .= PI.diag .* A .* P.diag'
    X_A_Xt!(Qh, Q, PI)

    # Predict the mean
    predict_mean!(x_pred, x, Ah)
    mul!(view(u_pred, :), SolProj, x_pred.μ)

    # Measure
    evaluate_ode!(integ, x_pred, tnew)

    # Estimate diffusion, and (if adaptive) the local error estimate; Stop here if rejected
    cache.local_diffusion = estimate_local_diffusion(cache.diffusionmodel, integ)
    if integ.opts.adaptive
        integ.EEst = compute_scaled_error_estimate!(integ, cache)
        if integ.EEst >= one(integ.EEst)
            return
        end
    end

    # Predict the covariance, using either the local or global diffusion
    extrapolation_diff =
        isdynamic(cache.diffusionmodel) ? cache.local_diffusion : cache.default_diffusion
    predict_cov!(x_pred, x, Ah, Qh, cache.C1, extrapolation_diff)

    # Compute measurement covariance only now; likelihood computation is currently broken
    compute_measurement_covariance!(cache)
    # cache.log_likelihood = logpdf(cache.measurement, zeros(d))
    # integ.sol.log_likelihood += integ.cache.log_likelihood

    # Update state and save the ODE solution value
    x_filt = update!(integ, x_pred)
    mul!(view(u_filt, :), SolProj, x_filt.μ)
    integ.u .= u_filt

    # Update the global diffusion MLE (if applicable)
    if !isdynamic(cache.diffusionmodel)
        cache.global_diffusion = estimate_global_diffusion(cache.diffusionmodel, integ)
    end

    # Advance the state
    copy!(integ.cache.x, integ.cache.x_filt)

    return nothing
end

"""
    evaluate_ode!(integ, x_pred, t, second_order::Val)

Evaluate the ODE vector field and, if necessary, the Jacobian. In addition, compute the
measurement mean (`z`) and the measurement function Jacobian (`H`). Results are saved into
`integ.cache.du`, `integ.cache.ddu`, `integ.cache.measurement.μ` and `integ.cache.H`.
Jacobians are computed either with the supplied `f.jac`, or via automatic differentiation,
following the example of OrdinaryDiffEq.jl.

For second-order ODEs and the `EK1FDB` algorithm, a specialized implementation is called.
"""
evaluate_ode!(integ, x_pred, t) =
    evaluate_ode!(integ, x_pred, t, Val(integ.f isa DynamicalODEFunction))
function evaluate_ode!(
    integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK},
    x_pred,
    t,
    second_order::Val{false},
)
    @unpack f, p, dt, alg = integ
    @unpack u_pred, du, ddu, measurement, R, H = integ.cache
    @assert iszero(R)

    @unpack E0, E1, E2 = integ.cache

    z, S = measurement.μ, measurement.Σ

    # Mean
    _eval_f!(du, u_pred, p, t, f)
    integ.destats.nf += 1
    # z .= MM*E1*x_pred.μ .- du
    if f.mass_matrix == I
        H .= E1
    elseif f.mass_matrix isa UniformScaling
        H .= f.mass_matrix.λ .* E1
    else
        _matmul!(H, f.mass_matrix, E1)
    end

    _matmul!(z, H, x_pred.μ)
    z .-= du[:]

    # If EK1, evaluate the Jacobian and adjust H
    if alg isa EK1 || alg isa IEKS
        u_lin =
            (alg isa IEKS && !isnothing(alg.linearize_at)) ? alg.linearize_at(t).μ : u_pred

        # Jacobian is computed either with the given jac, or ForwardDiff
        if !isnothing(f.jac)
            _eval_f_jac!(ddu, u_lin, p, t, f)
        else
            !isnothing(f.jac)
            @unpack du1, uf, jac_config = integ.cache
            uf.f = OrdinaryDiffEq.nlsolve_f(f, alg)
            uf.t = t
            uf.p = p
            if isinplace(f)
                OrdinaryDiffEq.jacobian!(ddu, uf, u_lin, du1, integ, jac_config)
            else
                ddu .= OrdinaryDiffEq.jacobian(uf, u_lin, integ)
            end
        end
        integ.destats.njacs += 1

        # _matmul!(H, f.mass_matrix, E1) # This is already the case (see above)
        _matmul!(H, ddu, E0, -1.0, 1.0)
    end

    return nothing
end

function evaluate_ode!(
    integ::OrdinaryDiffEq.ODEIntegrator{<:EK1FDB},
    x_pred,
    t,
    second_order::Val{false},
)
    @unpack f, p, dt, alg = integ
    @unpack d, u_pred, du, ddu, measurement, R, H = integ.cache
    @assert iszero(R)

    @unpack E0, E1, E2 = integ.cache

    z, S = measurement.μ, measurement.Σ

    (f.mass_matrix != I) && error("EK1FDB does not support mass-matrices right now")

    # Mean
    _eval_f!(du, u_pred, p, t, f)
    integ.destats.nf += 1
    # z .= MM*E1*x_pred.μ .- du
    H1, H2 = view(H, 1:d, :), view(H, d+1:2d, :)
    z1, z2 = view(z, 1:d), view(z, d+1:2d)

    H1 .= E1
    _matmul!(z1, H1, x_pred.μ)
    z1 .-= @view du[:]

    # If EK1, evaluate the Jacobian and adjust H
    u_lin = u_pred
    if !isnothing(f.jac)
        _eval_f_jac!(ddu, u_lin, p, t, f)
    elseif isinplace(f)
        ForwardDiff.jacobian!(ddu, (du, u) -> f(du, u, p, t), du, u_lin)
        integ.destats.nf += 1
    else
        ddu .= ForwardDiff.jacobian(u -> f(u, p, t), u_lin)
        integ.destats.nf += 1
    end
    integ.destats.njacs += 1
    _matmul!(H1, ddu, E0, -1.0, 1.0)

    z2 .= (E2 * x_pred.μ .- ddu * du)
    if integ.alg.jac_quality == 1
        # EK0-type approach
        H2 .= E2
    elseif integ.alg.jac_quality == 2
        H2 .= E2 - ddu * ddu * E0
    elseif integ.alg.jac_quality == 3
        _z2(m) = begin
            u_pred = E0 * m
            du = zeros(eltype(m), d)
            ddu = zeros(eltype(m), d, d)
            _eval_f!(du, u_pred, p, t, f)
            if !isnothing(f.jac)
                _eval_f_jac!(ddu, u_pred, p, t, f)
            elseif isinplace(f)
                ForwardDiff.jacobian!(ddu, (du, u) -> f(du, u, p, t), du, u_pred)
                # integ.destats.nf += 1
            else
                ddu .= ForwardDiff.jacobian(u -> f(u, p, t), u_pred)
                # integ.destats.nf += 1
            end
            # integ.destats.njacs += 1
            return (E2 * m .- ddu * du)
        end
        H2 .= ForwardDiff.jacobian(_z2, x_pred.μ)
    else
        error("EK1FDB's `jac_quality` has to be in [1,2,3]")
    end
    return nothing
end

function evaluate_ode!(
    integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK},
    x_pred,
    t,
    second_order::Val{true},
)
    @unpack f, p, dt, alg = integ
    @unpack d, u_pred, du, ddu, measurement, R, H = integ.cache
    @assert iszero(R)
    du2 = du

    @unpack E0, E1, E2 = integ.cache

    z, S = measurement.μ, measurement.Σ

    # Mean
    # _u_pred = E0 * x_pred.μ
    # _du_pred = E1 * x_pred.μ
    if isinplace(f)
        f.f1(du2, view(u_pred, 1:d), view(u_pred, d+1:2d), p, t)
    else
        du2 .= f.f1(view(u_pred, 1:d), view(u_pred, d+1:2d), p, t)
    end
    integ.destats.nf += 1
    _matmul!(z, E2, x_pred.μ)
    z .-= @view du2[:]

    # Cov
    if alg isa EK1
        (alg isa IEKS) && error("IEKS is currently not supported for SecondOrderODEProbems")

        if isinplace(f)
            H .= E2

            J = ddu
            ForwardDiff.jacobian!(
                J,
                (du2, du_u) -> f.f1(du2, view(du_u, 1:d), view(du_u, d+1:2d), p, t),
                du2,
                u_pred,
            )
            integ.destats.nf += 1
            integ.destats.njacs += 1
            _matmul!(H, J, integ.cache.SolProj, -1.0, 1.0)
        else
            J = ForwardDiff.jacobian(
                (du_u) -> f.f1(view(du_u, 1:d), view(du_u, d+1:2d), p, t),
                u_pred,
            )
            integ.destats.nf += 1
            integ.destats.njacs += 1
            H .= E2 .- J * integ.cache.SolProj
        end
    end

    return measurement
end

# Helper functions to handle both IIP and OOP easily
_eval_f!(du, u, p, t, f::AbstractODEFunction{true}) = f(du, u, p, t)
_eval_f!(du, u, p, t, f::AbstractODEFunction{false}) = (du .= f(u, p, t))
_eval_f_jac!(ddu, u, p, t, f::AbstractODEFunction{true}) = f.jac(ddu, u, p, t)
_eval_f_jac!(ddu, u, p, t, f::AbstractODEFunction{false}) = (ddu .= f.jac(u, p, t))

compute_measurement_covariance!(cache) =
    X_A_Xt!(cache.measurement.Σ, cache.x_pred.Σ, cache.H)

function update!(integ, prediction)
    @unpack measurement, H, R, x_filt = integ.cache
    @unpack K1, K2, x_tmp2, m_tmp = integ.cache
    update!(x_filt, prediction, measurement, H, K1, x_tmp2.Σ.mat, m_tmp)
    return x_filt
end

"""
    compute_scaled_error_estimate!(integ, cache)

Computes the scaled, local error estimate, that should satisfy `Eest < 1`. The actual local
error is computed with [`estimate_errors!`](@ref). Then, `DiffEqBase.calculate_residuals!`
handles the scaling with adaptive and relative tolerances, and `integ.opts.internalnorm`
provides the norm that should be used to return only a scalar.
"""
function compute_scaled_error_estimate!(integ, cache)
    @unpack u_pred, err_tmp = cache
    t = integ.t + integ.dt
    err_est_unscaled = estimate_errors!(cache)
    if integ.f isa DynamicalODEFunction # second-order ODE
        DiffEqBase.calculate_residuals!(
            err_tmp,
            integ.dt * err_est_unscaled,
            integ.u[1, :],
            u_pred[1, :],
            integ.opts.abstol,
            integ.opts.reltol,
            integ.opts.internalnorm,
            t,
        )
    else # regular first-order ODE
        DiffEqBase.calculate_residuals!(
            err_tmp,
            integ.dt * err_est_unscaled,
            integ.u,
            u_pred,
            integ.opts.abstol,
            integ.opts.reltol,
            integ.opts.internalnorm,
            t,
        )
    end
    return integ.opts.internalnorm(err_tmp, t) # scalar
end

"""
    estimate_errors!(cache)

Computes a local error estimate, as
```math
E_i = ( σ_{loc}^2 ⋅ (H Q(h) H^T)_{ii} )^(1/2)
```
To save allocations, the function modifies the given `cache` and writes into
`cache.m_tmp.Σ.squareroot` during some computations.
"""
function estimate_errors!(cache::GaussianODEFilterCache)
    @unpack local_diffusion, Qh, H, d = cache

    if local_diffusion isa Real && isinf(local_diffusion)
        return Inf
    end

    L = cache.m_tmp.Σ.squareroot

    if local_diffusion isa Diagonal
        _matmul!(L, H, sqrt.(local_diffusion) * Qh.squareroot)
        error_estimate = sqrt.(diag(L * L'))
        return view(error_estimate, 1:d)

    elseif local_diffusion isa Number
        _matmul!(L, H, Qh.squareroot)
        # error_estimate = local_diffusion .* diag(L*L')
        @tullio error_estimate[i] := L[i, j] * L[i, j]
        error_estimate .*= local_diffusion

        # @info "it's small anyways I guess?" error_estimate cache.measurement.μ .^ 2
        # error_estimate .+= cache.measurement.μ .^ 2

        error_estimate .= sqrt.(error_estimate)
        return view(error_estimate, 1:d)
    end
end