text:
__precompile__()

module PyCall

export pycall, pyimport, pybuiltin, PyObject, PyReverseDims,
       PyPtr, pyincref, pydecref, pyversion, PyArray, PyArray_Info,
       pyerr_check, pyerr_clear, pytype_query, PyAny, @pyimport, PyDict,
       pyisinstance, pywrap, pytypeof, pyeval, PyVector, pystring, pystr, pyrepr,
       pyraise, pytype_mapping, pygui, pygui_start, pygui_stop,
       pygui_stop_all, @pylab, set!, PyTextIO, @pysym, PyNULL, @pydef,
       pyimport_conda, @py_str, @pywith, @pycall, pybytes

import Base: size, ndims, similar, copy, getindex, setindex!, stride,
       convert, pointer, summary, convert, show, haskey, keys, values,
       eltype, get, delete!, empty!, length, isempty, start, done,
       next, filter!, hash, splice!, pop!, ==, isequal, push!,
       unshift!, shift!, append!, insert!, prepend!, mimewritable

# Python C API is not interrupt-save.  In principle, we should
# use sigatomic for every ccall to the Python library, but this
# should really be fixed in Julia (#2622).  However, we will
# use the sigatomic_begin/end functions to protect pycall and
# similar long-running (or potentially long-running) code.
import Base: sigatomic_begin, sigatomic_end

## Compatibility import for v0.5
using Compat
import Conda
import Base.unsafe_convert
import MacroTools   # because of issue #270

if isdefined(Base, :Iterators)
    import Base.Iterators: filter
end

#########################################################################

include(joinpath(dirname(@__FILE__), "..", "deps","depsutils.jl"))
include("startup.jl")

#########################################################################

# Mirror of C PyObject struct (for non-debugging Python builds).
# We won't actually access these fields directly; we'll use the Python
# C API for everything.  However, we need to define a unique Ptr type
# for PyObject*, and we might as well define the actual struct layout
# while we're at it.
immutable PyObject_struct
    ob_refcnt::Int
    ob_type::Ptr{Void}
end

const PyPtr = Ptr{PyObject_struct} # type for PythonObject* in ccall

const PyPtr_NULL = PyPtr(C_NULL)

#########################################################################
# Wrapper around Python's C PyObject* type, with hooks to Python reference
# counting and conversion routines to/from C and Julia types.
"""
    PyObject(juliavar)

This converts a julia variable to a PyObject, which is a reference to a Python object.
You can convert back to native julia types using `convert(T, o::PyObject)`, or using `PyAny(o)`.

Given `o::PyObject`, `o[:attribute]` is equivalent to `o.attribute` in Python, with automatic type conversion.

Given `o::PyObject`, `get(o, key)` is equivalent to `o[key]` in Python, with automatic type conversion.
"""
type PyObject
    o::PyPtr # the actual PyObject*
    function PyObject(o::PyPtr)
        po = new(o)
        finalizer(po, pydecref)
        return po
    end
end

"""
    PyNULL()

Return a `PyObject` that has a `NULL` underlying pointer, i.e. it doesn't
actually refer to any Python object.

This is useful for initializing `PyObject` global variables and array elements
before an actual Python object is available.   For example, you might do `const
myglobal = PyNULL()` and later on (e.g. in a module `__init__` function),
reassign `myglobal` to point to an actual object with `copy!(myglobal,
someobject)`.   This procedure will properly handle Python's reference counting
(so that the Python object will not be freed until you are done with
`myglobal`).
"""
PyNULL() = PyObject(PyPtr_NULL)

function pydecref(o::PyObject)
    ccall(@pysym(:Py_DecRef), Void, (PyPtr,), o.o)
    o.o = PyPtr_NULL
    o
end

function pyincref_(o::PyPtr)
    ccall((@pysym :Py_IncRef), Void, (PyPtr,), o)
    return o
end

function pyincref(o::PyObject)
    pyincref_(o.o)
    o
end

# doing an incref *before* creating a PyObject may safer in the
# case of borrowed references, to ensure that no exception or interrupt
# induces a double decref.
pyincref(o::PyPtr) = PyObject(pyincref_(o))

"""
"Steal" a reference from a PyObject: return the raw PyPtr, while
setting the corresponding `o.o` field to `NULL` so that no decref
will be performed when `o` is garbage collected.  (This means that
you can no longer use `o`.)  Used for passing objects to Python.
"""
function pystealref!(o::PyObject)
    optr = o.o
    o.o = PyPtr_NULL # don't decref when o is gc'ed
    return optr
end

function Base.copy!(dest::PyObject, src::PyObject)
    pydecref(dest)
    dest.o = src.o
    return pyincref(dest)
end

pyisinstance(o::PyObject, t::PyObject) =
  t.o != C_NULL && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t.o) == 1

pyisinstance(o::PyObject, t::Union{Ptr{Void},PyPtr}) =
  t != C_NULL && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t) == 1

pyquery(q::Ptr{Void}, o::PyObject) =
  ccall(q, Cint, (PyPtr,), o) == 1

# conversion to pass PyObject as ccall arguments:
unsafe_convert(::Type{PyPtr}, po::PyObject) = po.o

# use constructor for generic conversions to PyObject
convert(::Type{PyObject}, o) = PyObject(o)
convert(::Type{PyObject}, o::PyObject) = o
PyObject(o::PyObject) = o

#########################################################################

include("pyinit.jl")
include("exception.jl")
include("gui.jl")

pytypeof(o::PyObject) = o.o == C_NULL ? throw(ArgumentError("NULL PyObjects have no Python type")) : PyObject(@pycheckn ccall(@pysym(:PyObject_Type), PyPtr, (PyPtr,), o))

#########################################################################

include("gc.jl")

# make a PyObject that embeds a reference to keep, to prevent Julia
# from garbage-collecting keep until o is finalized.
PyObject(o::PyPtr, keep::Any) = pyembed(PyObject(o), keep)

#########################################################################

@compat TypeTuple{N} = Union{Type,NTuple{N, Type}}
include("pybuffer.jl")
include("conversions.jl")
include("pytype.jl")
include("pyiterator.jl")
include("pyclass.jl")
include("callback.jl")
include("io.jl")

#########################################################################
# Pretty-printing PyObject

"""
    pystr(o::PyObject)

Return a string representation of `o` corresponding to `str(o)` in Python.
"""
pystr(o::PyObject) = convert(AbstractString,
                             PyObject(@pycheckn ccall((@pysym :PyObject_Str), PyPtr,
                                                      (PyPtr,), o)))

"""
    pyrepr(o::PyObject)

Return a string representation of `o` corresponding to `repr(o)` in Python.
"""
pyrepr(o::PyObject) = convert(AbstractString,
                              PyObject(@pycheckn ccall((@pysym :PyObject_Repr), PyPtr,
                                                       (PyPtr,), o)))

"""
    pystring(o::PyObject)

Return a string representation of `o`.  Normally, this corresponds to `repr(o)`
in Python, but unlike `repr` it should never fail (falling back to `str` or to
printing the raw pointer if necessary).
"""
function pystring(o::PyObject)
    if o.o == C_NULL
        return "NULL"
    else
        s = ccall((@pysym :PyObject_Repr), PyPtr, (PyPtr,), o)
        if (s == C_NULL)
            pyerr_clear()
            s = ccall((@pysym :PyObject_Str), PyPtr, (PyPtr,), o)
            if (s == C_NULL)
                pyerr_clear()
                return string(o.o)
            end
        end
        return convert(AbstractString, PyObject(s))
    end
end

function show(io::IO, o::PyObject)
    print(io, "PyObject $(pystring(o))")
end

function Base.Docs.doc(o::PyObject)
    if haskey(o, "__doc__")
        d = o["__doc__"]
        if d.o != pynothing[]
            return Base.Docs.Text(convert(AbstractString, d))
        end
    end
    return Base.Docs.Text("Python object (no docstring found)")
end

#########################################################################
# computing hashes of PyObjects

const pysalt = hash("PyCall.PyObject") # "salt" to mix in to PyObject hashes
hashsalt(x) = hash(x, pysalt)

function hash(o::PyObject)
    if o.o == C_NULL
        hashsalt(C_NULL)
    elseif is_pyjlwrap(o)
        # call native Julia hash directly on wrapped Julia objects,
        # since on 64-bit Windows the Python 2.x hash is only 32 bits
        hashsalt(unsafe_pyjlwrap_to_objref(o.o))
    else
        h = ccall((@pysym :PyObject_Hash), Py_hash_t, (PyPtr,), o)
        if h == -1 # error
            pyerr_clear()
            return hashsalt(o.o)
        end
        hashsalt(h)
    end
end

#########################################################################
# For o::PyObject, make o["foo"] and o[:foo] equivalent to o.foo in Python,
# with the former returning an raw PyObject and the latter giving the PyAny
# conversion.

function getindex(o::PyObject, s::AbstractString)
    if (o.o == C_NULL)
        throw(ArgumentError("ref of NULL PyObject"))
    end
    p = ccall((@pysym :PyObject_GetAttrString), PyPtr, (PyPtr, Cstring), o, s)
    if p == C_NULL
        pyerr_clear()
        throw(KeyError(s))
    end
    return PyObject(p)
end

getindex(o::PyObject, s::Symbol) = convert(PyAny, getindex(o, string(s)))

function setindex!(o::PyObject, v, s::Union{Symbol,AbstractString})
    if (o.o == C_NULL)
        throw(ArgumentError("assign of NULL PyObject"))
    end
    if -1 == ccall((@pysym :PyObject_SetAttrString), Cint,
                   (PyPtr, Cstring, PyPtr), o, s, PyObject(v))
        pyerr_clear()
        throw(KeyError(s))
    end
    o
end

function haskey(o::PyObject, s::Union{Symbol,AbstractString})
    if (o.o == C_NULL)
        throw(ArgumentError("haskey of NULL PyObject"))
    end
    return 1 == ccall((@pysym :PyObject_HasAttrString), Cint,
                      (PyPtr, Cstring), o, s)
end

#########################################################################

keys(o::PyObject) = Symbol[m[1] for m in pycall(inspect["getmembers"],
                                PyVector{Tuple{Symbol,PyObject}}, o)]

#########################################################################
# Create anonymous composite w = pywrap(o) wrapping the object o
# and providing access to o's members (converted to PyAny) as w.member.

# we skip wrapping Julia reserved words (which cannot be type members)
const reserved = Set{String}(["while", "if", "for", "try", "return", "break", "continue", "function", "macro", "quote", "let", "local", "global", "const", "abstract", "typealias", "type", "bitstype", "immutable", "ccall", "do", "module", "baremodule", "using", "import", "export", "importall", "pymember", "false", "true", "Tuple"])

"""
    pywrap(o::PyObject)

This returns a wrapper `w` that is an anonymous module which provides (read) access to converted versions of o's members as w.member.

For example, `@pyimport module as name` is equivalent to `const name = pywrap(pyimport("module"))`

If the Python module contains identifiers that are reserved words in Julia (e.g. function), they cannot be accessed as `w.member`; one must instead use `w.pymember(:member)` (for the PyAny conversion) or w.pymember("member") (for the raw PyObject).
"""
function pywrap(o::PyObject, mname::Symbol=:__anon__)
    members = convert(Vector{Tuple{AbstractString,PyObject}},
                      pycall(inspect["getmembers"], PyObject, o))
    filter!(m -> !(m[1] in reserved), members)
    m = Module(mname, false)
    consts = [Expr(:const, Expr(:(=), Symbol(x[1]), convert(PyAny, x[2]))) for x in members]
    exports = try
                  convert(Vector{Symbol}, o["__all__"])
              catch
                  [Symbol(x[1]) for x in filter(x -> x[1][1] != '_', members)]
              end
    eval(m, Expr(:toplevel, consts..., :(pymember(s) = $(getindex)($(o), s)),
                 Expr(:export, exports...)))
    m
end

#########################################################################

@static if Compat.Sys.iswindows()
    # Many python extensions are linked against a very specific version of the
    # MSVC runtime library. To load this library, libpython declares an
    # appropriate manifest, but unfortunately most extensions do not.
    # Libpython itself does extend its activation context to any extensions it
    # loads, but some python libraries (e.g. libzmq), load such extensions
    # through ctypes, which does not have this functionality. Work around
    # this by manually activating the activation context before any call to
    # PyImport_ImportModule, since extensions are most likely to be loaded
    # during import.

    immutable ACTIVATION_CONTEXT_BASIC_INFORMATION
        handle::Ptr{Void}
        dwFlags::UInt32
    end
    const ActivationContextBasicInformation = 1
    const QUERY_ACTCTX_FLAG_ACTCTX_IS_ADDRESS = 0x10
    const PyActCtx = Ref{Ptr{Void}}(C_NULL)

    function ComputePyActCtx()
        if PyActCtx[] == C_NULL
            some_address_in_libpython = @pyglobal(:PyImport_ImportModule)
            ActCtxBasicInfo = Ref{ACTIVATION_CONTEXT_BASIC_INFORMATION}()
            succeeded = ccall(:QueryActCtxW,stdcall,Bool,
                (UInt32,Ptr{Void},Ptr{Void},Culong,Ptr{Void},Csize_t,Ptr{Csize_t}),
                QUERY_ACTCTX_FLAG_ACTCTX_IS_ADDRESS, some_address_in_libpython,
                C_NULL, ActivationContextBasicInformation, ActCtxBasicInfo,
                sizeof(ActCtxBasicInfo), C_NULL)
            @assert succeeded
            PyActCtx[] = ActCtxBasicInfo[].handle
        end
        PyActCtx[]
    end

    function ActivatePyActCtx()
        cookie = Ref{Ptr{Void}}()
        succeeded = ccall(:ActivateActCtx,stdcall, Bool,
                          (Ptr{Void}, Ptr{Ptr{Void}}),
                          ComputePyActCtx(), cookie)
        @assert succeeded
        cookie[]
    end

    function DeactivatePyActCtx(cookie)
        succeeded = ccall(:DeactivateActCtx, stdcall, Bool,
                         (UInt32, Ptr{Void}), 0, cookie)
        @assert succeeded
    end
else
    ActivatePyActCtx() = nothing
    DeactivatePyActCtx(cookie) = nothing
end

"""
    pyimport(s::AbstractString)

Import the Python module `s` (a string or symbol) and return a pointer to it (a `PyObject`). Functions or other symbols in the module may then be looked up by s[name] where name is a string (for the raw PyObject) or symbol (for automatic type-conversion). Unlike the @pyimport macro, this does not define a Julia module and members cannot be accessed with `s.name`
"""
function pyimport(name::AbstractString)
    cookie = ActivatePyActCtx()
    try
        return PyObject(@pycheckn ccall((@pysym :PyImport_ImportModule), PyPtr,
                                        (Cstring,), name))
    catch e
        # on importerror, add a bit more information to the
        # exception message to help with common user confusions
        # about installing packages.
        if isa(e, PyError) && e.T.o == @pyglobalobjptr(:PyExc_ImportError)
            msg = """
The Python package $name could not be found by pyimport. Usually this means
that you did not install $name in the Python version being used by PyCall.

"""
            if conda
                msg = msg * """
PyCall is currently configured to use the Julia-specific Python distribution
installed by the Conda.jl package.  To install the $name module, you can
use `pyimport_conda("$(escape_string(name))", PKG)`, where PKG is the Anaconda
package the contains the module $name, or alternatively you can use the
Conda package directly (via `using Conda` followed by `Conda.add` etcetera).

Alternatively, if you want to use a different Python distribution on your
system, such as a system-wide Python (as opposed to the Julia-specific Python),
you can re-configure PyCall with that Python.   As explained in the PyCall
documentation, set ENV["PYTHON"] to the path/name of the python executable
you want to use, run Pkg.build("PyCall"), and re-launch Julia.
"""
            else
                msg = msg * """
PyCall is currently configured to use the Python version at:

$python

and you should use whatever mechanism you usually use (apt-get, pip, conda,
etcetera) to install the Python package containing the $name module.

One alternative is to re-configure PyCall to use a different Python
version on your system: set ENV["PYTHON"] to the path/name of the python
executable you want to use, run Pkg.build("PyCall"), and re-launch Julia.

Another alternative is to configure PyCall to use a Julia-specific Python
distribution via the Conda.jl package (which installs a private Anaconda
Python distribution), which has the advantage that packages can be installed
and kept up-to-date via Julia.  As explained in the PyCall documentation,
set ENV["PYTHON"]="", run Pkg.build("PyCall"), and re-launch Julia. Then,
To install the $name module, you can use `pyimport_conda("$(escape_string(name))", PKG)`,
where PKG is the Anaconda package the contains the module $name,
or alternatively you can use the Conda package directly (via
`using Conda` followed by `Conda.add` etcetera).
"""
            end
            e.msg *= "\n\n" * msg * "\n"
        end
        rethrow(e)
    finally
        DeactivatePyActCtx(cookie)
    end
    obj
end
pyimport(name::Symbol) = pyimport(string(name))

# convert expressions like :math or :(scipy.special) into module name strings
modulename(s::QuoteNode) = modulename(s.value)
modulename(s::Symbol) = string(s)
function modulename(e::Expr)
    if e.head == :.
        string(modulename(e.args[1]), :., modulename(e.args[2]))
    elseif e.head == :quote
        modulename(e.args...)
    else
        throw(ArgumentError("invalid module"))
    end
end

# separate this function in order to make it easier to write more
# pyimport-like functions
function pyimport_name(name, optional_varname)
    len = length(optional_varname)
    if len > 0 && (len != 2 || optional_varname[1] != :as)
        throw(ArgumentError("usage @pyimport module [as name]"))
    elseif len == 2
        optional_varname[2]
    elseif typeof(name) == Symbol
        name
    else
        mname = modulename(name)
        throw(ArgumentError("$mname is not a valid module variable name, use @pyimport $mname as <name>"))
    end
end

macro pyimport(name, optional_varname...)
    mname = modulename(name)
    Name = pyimport_name(name, optional_varname)
    quoteName = Expr(:quote, Name)
    # VERSION 0.7
    @static if isdefined(Base, Symbol("@isdefined"))
        isdef_check = :(isdefined($__module__, $quoteName))
    else
        isdef_check = :(isdefined($quoteName))
    end
    quote
        if !$isdef_check
            const $(esc(Name)) = pywrap(pyimport($mname))
        elseif !isa($(esc(Name)), Module)
            error("@pyimport: ", $(Expr(:quote, Name)), " already defined")
        end
        nothing
    end
end

#########################################################################

"""
    @pywith

Mimics a Python 'with' statement. Usage:

@pywith EXPR[::TYPE1] [as VAR[::TYPE2]] begin
    BLOCK
end

TYPE1/TYPE2 can optionally be used to override automatic conversion to Julia
types for both the context manager and variable in cases where this is not
desired.

"""
macro pywith(EXPR,as,VAR,BLOCK)
    if isexpr(VAR,:(::))
        VAR,TYPE = VAR.args
    else
        TYPE = PyAny
    end
    if !((as==:as) && isa(VAR,Symbol))
        throw(ArgumentError("usage: @pywith EXPR[::TYPE1] [as VAR[::TYPE2]] BLOCK."))
    end
    _pywith(EXPR,VAR,TYPE,BLOCK)
end
macro pywith(EXPR,BLOCK)
    _pywith(EXPR,nothing,PyObject,BLOCK)
end

function _pywith(EXPR,VAR,TYPE,BLOCK)
    EXPR_str = string(EXPR)
    quote
        mgr = $(isexpr(EXPR,:(::)) ? :(@pycall $(esc(EXPR))) : esc(EXPR))
        if !isa(mgr,PyObject)
            if $(isexpr(EXPR,:call))
                throw(ArgumentError("@pywith: `$($EXPR_str)` did not return a PyObject. If this is a call to a Python function, try `$($EXPR_str)::PyObject` to turn off automatic conversion."))
            else
                throw(ArgumentError("@pywith: `$($EXPR_str)` should be a PyObject."))
            end
        end
        mgrT = pytypeof(mgr)
        exit = mgrT["__exit__"]
        value = @pycall mgrT["__enter__"](mgr)::$(esc(TYPE))
        exc = true
        try
            try
                $(VAR==nothing ? :() : :($(esc(VAR)) = value))
                $(esc(BLOCK))
            catch err
                exc = false
                if !(@pycall exit(mgr, pyimport(:sys)[:exc_info]()...)::Bool)
                    throw(err)
                end
            end
        finally
            if exc
                exit(mgr, nothing, nothing, nothing)
            end
        end
        nothing
    end
end

#########################################################################

"""
    anaconda_conda()

Return the path of the `conda` program if PyCall is configured to use
an Anaconda install (but *not* the Conda.jl Python), and the empty
string otherwise.
"""
function anaconda_conda()
    # Anaconda Python seems to always include "Anaconda" in the version string.
    if conda || !contains(unsafe_string(ccall(@pysym(:Py_GetVersion), Ptr{UInt8}, ())), "conda")
        return ""
    end
    aconda = joinpath(dirname(pyprogramname), "conda")
    return isfile(aconda) ? aconda : ""
end

"""
    pyimport_conda(modulename, condapkg, [channel])

Returns the result of `pyimport(modulename)` if possible.   If the module
is not found, and PyCall is configured to use the Conda Python distro
(via the Julia Conda package), then automatically install `condapkg`
via `Conda.add(condapkg)` and then re-try the `pyimport`.   Other
Anaconda-based Python installations are also supported as long as
their `conda` program is functioning.

If PyCall is not using Conda and the `pyimport` fails, throws
an exception with an error message telling the user how to configure
PyCall to use Conda for automated installation of the module.

The third argument, `channel` is an optional Anaconda "channel" to use
for installing the package; this is useful for packages that are not
included in the default Anaconda package listing.
"""
function pyimport_conda(modulename::AbstractString, condapkg::AbstractString,
                        channel::AbstractString="")
    try
        pyimport(modulename)
    catch e
        if conda
            info("Installing $modulename via the Conda $condapkg package...")
            isempty(channel) || Conda.add_channel(channel)
            Conda.add(condapkg)
            pyimport(modulename)
        else
            aconda = anaconda_conda()
            if !isempty(aconda)
                try
                    info("Installing $modulename via Anaconda's $aconda...")
                    isempty(channel) || run(`$aconda config --add channels $channel --force`)
                    run(`$aconda install -y $condapkg`)
                catch e2
                    error("""
                    Detected Anaconda Python, but \"$aconda install -y $condapkg\"  failed.

                    You should either install $condapkg manually or fix your Anaconda installation so that $aconda works.   Alternatively, you can reconfigure PyCall to use its own Miniconda installation via the Conda.jl package, by running:
                        ENV["PYTHON"]=""; Pkg.build("PyCall")
                    after re-launching Julia.

                    The exception was: $e2
                    """)
                end
                return pyimport(modulename)
            end

            error("""
                Failed to import required Python module $modulename.

                For automated $modulename installation, try configuring PyCall to use the Conda.jl package's Python "Miniconda" distribution within Julia. Relaunch Julia and run:
                    ENV["PYTHON"]=""
                    Pkg.build("PyCall")
                before trying again.

                The pyimport exception was: $e
                """)
        end
    end
end

#########################################################################

# look up a global builtin
"""
    pybuiltin(s::AbstractString)

Look up a string or symbol `s` among the global Python builtins. If `s` is a string it returns a PyObject, while if `s` is a symbol it returns the builtin converted to `PyAny`.
"""
function pybuiltin(name)
    builtin[name]
end

#########################################################################

"""
Low-level version of `pycall(o, ...)` that always returns `PyObject`.
"""
function _pycall(o::Union{PyObject,PyPtr}, args...; kwargs...)
    oargs = map(PyObject, args)
    nargs = length(args)
    sigatomic_begin()
    try
        arg = PyObject(@pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,),
                                       nargs))
        for i = 1:nargs
            @pycheckz ccall((@pysym :PyTuple_SetItem), Cint,
                             (PyPtr,Int,PyPtr), arg, i-1, oargs[i])
            pyincref(oargs[i]) # PyTuple_SetItem steals the reference
        end
        if isempty(kwargs)
            ret = PyObject(@pycheckn ccall((@pysym :PyObject_Call), PyPtr,
                                          (PyPtr,PyPtr,PyPtr), o, arg, C_NULL))
        else
            #kw = PyObject((AbstractString=>Any)[string(k) => v for (k, v) in kwargs])
            kw = PyObject(Dict{AbstractString, Any}([Pair(string(k), v) for (k, v) in kwargs]))
            ret = PyObject(@pycheckn ccall((@pysym :PyObject_Call), PyPtr,
                                            (PyPtr,PyPtr,PyPtr), o, arg, kw))
        end
        return ret::PyObject
    finally
        sigatomic_end()
    end
end

"""
    pycall(o::Union{PyObject,PyPtr}, returntype::TypeTuple, args...; kwargs...)

Call the given Python function (typically looked up from a module) with the given args... (of standard Julia types which are converted automatically to the corresponding Python types if possible), converting the return value to returntype (use a returntype of PyObject to return the unconverted Python object reference, or of PyAny to request an automated conversion)
"""
pycall(o::Union{PyObject,PyPtr}, returntype::TypeTuple, args...; kwargs...) =
    return convert(returntype, _pycall(o, args...; kwargs...))::returntype

pycall(o::Union{PyObject,PyPtr}, ::Type{PyAny}, args...; kwargs...) =
    return convert(PyAny, _pycall(o, args...; kwargs...))

(o::PyObject)(args...; kws...) = pycall(o, PyAny, args...; kws...)
(::Type{PyAny})(o::PyObject) = convert(PyAny, o)

"""
    @pycall func(args...)::T

Convenience macro which turns `func(args...)::T` into pycall(func, T, args...)
"""
macro pycall(ex)
    if !(isexpr(ex,:(::)) && isexpr(ex.args[1],:call))
        throw(ArgumentError("Usage: @pycall func(args...)::T"))
    end
    func = ex.args[1].args[1]
    args, kwargs = ex.args[1].args[2:end], []
    if isexpr(args[1],:parameters)
        kwargs, args = args[1], args[2:end]
    end
    T = ex.args[2]
    :(pycall($(map(esc,[kwargs; func; T; args])...)))
end

#########################################################################
# Once Julia lets us overload ".", we will use [] to access items, but
# for now we can define "get".

function get(o::PyObject, returntype::TypeTuple, k, default)
    r = ccall((@pysym :PyObject_GetItem), PyPtr, (PyPtr,PyPtr), o,PyObject(k))
    if r == C_NULL
        pyerr_clear()
        default
    else
        convert(returntype, PyObject(r))
    end
end

get(o::PyObject, returntype::TypeTuple, k) =
    convert(returntype, PyObject(@pycheckn ccall((@pysym :PyObject_GetItem),
                                 PyPtr, (PyPtr,PyPtr), o, PyObject(k))))

get(o::PyObject, k, default) = get(o, PyAny, k, default)
get(o::PyObject, k) = get(o, PyAny, k)

function delete!(o::PyObject, k)
    e = ccall((@pysym :PyObject_DelItem), Cint, (PyPtr, PyPtr), o, PyObject(k))
    if e == -1
        pyerr_clear() # delete! ignores errors in Julia
    end
    return o
end

function set!(o::PyObject, k, v)
    @pycheckz ccall((@pysym :PyObject_SetItem), Cint, (PyPtr, PyPtr, PyPtr),
                     o, PyObject(k), PyObject(v))
    v
end

#########################################################################
# Support [] for integer keys, and other duck-typed sequence/list operations,
# as those don't conflict with symbols/strings used for attributes.

# Index conversion: Python is zero-based.  It also has -1 based
# backwards indexing, but we don't support this, in favor of the
# Julian syntax o[end-1] etc.
function ind2py(i::Integer)
    i <= 0 && throw(BoundsError())
    return i-1
end

_getindex(o::PyObject, i::Integer, T) = convert(T, PyObject(@pycheckn ccall((@pysym :PySequence_GetItem), PyPtr, (PyPtr, Int), o, ind2py(i))))
getindex(o::PyObject, i::Integer) = _getindex(o, i, PyAny)
function setindex!(o::PyObject, v, i::Integer)
    @pycheckz ccall((@pysym :PySequence_SetItem), Cint, (PyPtr, Int, PyPtr), o, ind2py(i), PyObject(v))
    v
end
getindex(o::PyObject, i1::Integer, i2::Integer) = get(o, (ind2py(i1),ind2py(i2)))
setindex!(o::PyObject, v, i1::Integer, i2::Integer) = set!(o, (ind2py(i1),ind2py(i2)), v)
getindex(o::PyObject, I::Integer...) = get(o, map(ind2py, I))
setindex!(o::PyObject, v, I::Integer...) = set!(o, map(ind2py, I), v)
Base.endof(o::PyObject) = length(o)
length(o::PyObject) = @pycheckz ccall((@pysym :PySequence_Size), Int, (PyPtr,), o)

function splice!(a::PyObject, i::Integer)
    v = a[i]
    @pycheckz ccall((@pysym :PySequence_DelItem), Cint, (PyPtr, Int), a, i-1)
    v
end

pop!(a::PyObject) = splice!(a, length(a))
shift!(a::PyObject) = splice!(a, 1)

function empty!(a::PyObject)
    for i in length(a):-1:1
        @pycheckz ccall((@pysym :PySequence_DelItem), Cint, (PyPtr, Int), a, i-1)
    end
    a
end

# The following operations only work for the list type and subtypes thereof:
function push!(a::PyObject, item)
    @pycheckz ccall((@pysym :PyList_Append), Cint, (PyPtr, PyPtr),
                     a, PyObject(item))
    a
end

function insert!(a::PyObject, i::Integer, item)
    @pycheckz ccall((@pysym :PyList_Insert), Cint, (PyPtr, Int, PyPtr),
                     a, ind2py(i), PyObject(item))
    a
end

unshift!(a::PyObject, item) = insert!(a, 1, item)

function prepend!(a::PyObject, items)
    if isa(items,PyObject) && items.o == a.o
        # avoid infinite loop in prepending a to itself
        return prepend!(a, collect(items))
    end
    for (i,x) in enumerate(items)
        insert!(a, i, x)
    end
    a
end

function append!(a::PyObject, items)
    for item in items
        push!(a, item)
    end
    return a
end

append!(a::PyObject, items::PyObject) =
    PyObject(@pycheckn ccall((@pysym :PySequence_InPlaceConcat),
                             PyPtr, (PyPtr, PyPtr), a, items))

#########################################################################
# operators on Python objects

include("pyoperators.jl")

#########################################################################
# support IPython _repr_foo functions for MIME output of PyObjects

for (mime, method) in ((MIME"text/html", "_repr_html_"),
                       (MIME"text/markdown", "_repr_markdown_"),
                       (MIME"text/json", "_repr_json_"),
                       (MIME"application/javascript", "_repr_javascript_"),
                       (MIME"application/pdf", "_repr_pdf_"),
                       (MIME"image/jpeg", "_repr_jpeg_"),
                       (MIME"image/png", "_repr_png_"),
                       (MIME"image/svg+xml", "_repr_svg_"),
                       (MIME"text/latex", "_repr_latex_"))
    T = istextmime(mime()) ? AbstractString : Vector{UInt8}
    @eval begin
        function show(io::IO, mime::$mime, o::PyObject)
            if o.o != C_NULL && haskey(o, $method)
                r = pycall(o[$method], PyObject)
                r.o != pynothing[] && return write(io, convert($T, r))
            end
            throw(MethodError(show, (io, mime, o)))
        end
        mimewritable(::$mime, o::PyObject) =
            o.o != C_NULL && haskey(o, $method) && let meth = o[$method]
                meth.o != pynothing[] &&
                pycall(meth, PyObject).o != pynothing[]
            end
    end
end

#########################################################################
# Expose Python docstrings to the Julia doc system

if isdefined(Docs, :getdoc)  # Introduced in Julia .6
    Docs.getdoc(o::PyObject) = Text(String(o["__doc__"]))
end

#########################################################################

include("pyeval.jl")
include("serialize.jl")

#########################################################################
# Precompilation: just an optimization to speed up initialization.
# Here, we precompile functions that are passed to cfunction by __init__,
# for the reasons described in JuliaLang/julia#12256.

precompile(pyjlwrap_call, (PyPtr,PyPtr,PyPtr))
precompile(pyjlwrap_dealloc, (PyPtr,))
precompile(pyjlwrap_repr, (PyPtr,))
precompile(pyjlwrap_hash, (PyPtr,))
precompile(pyjlwrap_hash32, (PyPtr,))

# TODO: precompilation of the io.jl functions

end # module PyCall


text:
using Documenter, ReachabilityAnalysis

DocMeta.setdocmeta!(ReachabilityAnalysis, :DocTestSetup,
                   :(using ReachabilityAnalysis); recursive=true)

# pass --fast as an argument to skip rebuilding the examples and running doctests
const _FAST = findfirst(isequal("--fast"), ARGS) !== nothing

# ========================
# Generate examples
# ========================
if !_FAST
    include("generate.jl")
end

# ========================
# Tutorials
# ========================

SET_REPRESENTATIONS = ["Set representations" => [
    "Introduction"                     => "tutorials/set_representations/introduction.md",
    "Hyperrectangular sets"            => "tutorials/set_representations/hyperrectangles.md",
    "Zonotopic sets"                   => "tutorials/set_representations/zonotopes.md",
    "Polyhedral computations"          => "tutorials/set_representations/polyhedral_computations.md",
    "Lazy set representations"         => "tutorials/set_representations/lazy_sets.md",
    "Computing with support functions" => "tutorials/set_representations/support_functions.md",
    #"Other set operations"             => "tutorials/set_representations/other_set_operations.md",
    "LazySets type hierarchy"          => "tutorials/set_representations/lazysets_hierarchy.md",
    "Numerical tolerance"              => "tutorials/set_representations/tolerance.md",
    "Metric notions"                   => "tutorials/set_representations/distances.md"]]

LINEAR_METHODS = ["Linear methods" => [
    "Quickstart"                       => "tutorials/linear_methods/introduction.md",
    "Discrete time reachability"       => "tutorials/linear_methods/discrete_time.md",
    "Dense time reachability"          => "tutorials/linear_methods/dense_time.md",
#    "Propagating zonotopes"            => "tutorials/linear_methods/zonotopes_setprop.md",
#    "Propagating hyperrectangles"      => "tutorials/linear_methods/hyperrectangles_setprop.md",
#    "Propagating support functions"    => "tutorials/linear_methods/supfunc_setprop.md",
    #"Helicopter model"                 => "models/Helicopter.md",
    "Helicopter model"                 => "tutorials/linear_methods/helicopter.md",
    "Structural model"                 => "models/ISS.md"]]

UNCERTAIN_INPUTS = ["Modeling uncertain inputs" => [
    "Introduction"                     => "tutorials/uncertain_inputs/introduction.md",
    "Building model"                   => "models/Building.md",
    "Operational amplifier"            => "models/OpAmp.md"]]

TAYLOR_METHODS = ["Taylor methods"     => [
    "Introduction"                     => "tutorials/taylor_methods/introduction.md",
    "Taylor model reach-sets"          => "tutorials/taylor_methods/taylor_model_reachsets.md",
    "Domain splitting"                 => "tutorials/taylor_methods/domain_splitting.md",
    "Common gotchas"                   => "tutorials/taylor_methods/gotchas.md",
    "Lotka-Volterra"                   => "models/LotkaVolterra.md",
    "Van der Pol oscillator"           => "models/VanDerPol.md",
    "Lorenz system"                    => "models/Lorenz.md"]]

UNCERTAIN_PARAMETERS = ["Modeling uncertain parameters" => [
    "Introduction"                     => "tutorials/parametric_reachability/introduction.md",
    "Transmision line"                 => "models/TransmissionLine.md"]]

HYBRID_SYSTEMS = ["Hybrid systems" => [
    "Introduction"                     => "tutorials/hybrid_systems/introduction.md",
    "Square wave oscillator"           => "models/SquareWaveOscillator.md",
    "Thermostat model"                 => "tutorials/hybrid_systems/thermostat.md"]]

CLOCKED_SYSTEMS = ["Clocked systems"  => [
    "Introduction"                     => "tutorials/clocked_systems/introduction.md",
    "Platoon"                          => "models/Platoon.md",
    "Electro-mechanic break"           => "models/EMBrake.md"]]

BACKWARDS_REACHABILITY = ["Backwards reachability" => [
    "Introduction"                     => "tutorials/backwards_reachability/introduction.md",
    "Projectile motion"                => "models/Projectile.md"]]

LINEAR_PDE = ["Linear PDEs" => [
    "Introduction"                     => "tutorials/linear_pde/introduction.md",
    "Clamped beam"                     => "tutorials/linear_pde/clamped.md",
    "Heat transfer"                    => "tutorials/linear_pde/heat_transfer.md",
    "Concrete heat of hydration"       => "tutorials/linear_pde/concrete.md"]]
#   "Clamped beam"                     => "examples/Clamped.md",
#   "Heat transfer"                    => "examples/Heat1D.md",
#   "Concrete heat of hydration"       => "examples/Concrete.md"

LINEAR_SOLVERS = ["ASB07"    => "man/algorithms/ASB07.md",
                  "BFFPSV18" => "man/algorithms/BFFPSV18.md",
                  "BOX"      => "man/algorithms/BOX.md",
                  "GLGM06"   => "man/algorithms/GLGM06.md",
                  "INT"      => "man/algorithms/INT.md",
                  "LGG09"    => "man/algorithms/LGG09.md",
                  "ORBIT"    => "man/algorithms/ORBIT.md",
                  "VREP"     => "man/algorithms/VREP.md"]

LINEAR_SOLVERS_API = ["ASB07"    => "lib/algorithms/ASB07.md",
                      "BFFPSV18" => "lib/algorithms/BFFPSV18.md",
                      "BOX"      => "lib/algorithms/BOX.md",
                      "GLGM06"   => "lib/algorithms/GLGM06.md",
                      "INT"      => "lib/algorithms/INT.md",
                      "LGG09"    => "lib/algorithms/LGG09.md",
                      "ORBIT"    => "lib/algorithms/ORBIT.md",
                      "VREP"     => "lib/algorithms/VREP.md"]

NONLINEAR_SOLVERS = ["CARLIN"   => "man/algorithms/CARLIN.md",
                     "QINT"     => "man/algorithms/QINT.md",
                     "TMJets"   => "man/algorithms/TMJets.md"]

NONLINEAR_SOLVERS_API = ["CARLIN"   => "lib/algorithms/CARLIN.md",
                         "QINT"     => "lib/algorithms/QINT.md",
                         "TMJets"   => "lib/algorithms/TMJets.md"]

# ========================
# Docs contents
# ========================

makedocs(
    format = Documenter.HTML(prettyurls = haskey(ENV, "GITHUB_ACTIONS"),
                             collapselevel = 1,
                             assets = ["assets/aligned.css"]),
    sitename = "ReachabilityAnalysis.jl",
    doctest = false,
    strict = false,
    pages = [
        "Overview" => "index.md",

        "Tutorials"    => vcat(SET_REPRESENTATIONS,
                          LINEAR_METHODS,
                          UNCERTAIN_INPUTS,
                          TAYLOR_METHODS,
                          UNCERTAIN_PARAMETERS,
                          LINEAR_PDE,
                          HYBRID_SYSTEMS,
                          CLOCKED_SYSTEMS),

        "Manual"       => ["Systems"           => "man/systems.md",
                           "Reach-sets"        => "man/reachsets.md",
                           "Flowpipes"         => "man/flowpipes.md",
                           "Linear solvers"    => LINEAR_SOLVERS,
                           "Nonlinear solvers" => NONLINEAR_SOLVERS,
                           "Solutions"         => "man/solutions.md",
                           "Invariants"        => "man/invariants.md",
                           "Visualization"     => "man/visualization.md",
                           "Projections"       => "man/projections.md",
                           "Clustering"        => "man/clustering.md"],

        "Benchmarks" => ["Repeatability evaluations"   => "man/benchmarks/repeatability.md",
                         "Model library"               => "man/benchmarks/model_library.md",
                         "Comparison with other tools" => "man/benchmarks/comparison.md",
                         "Benchmark repository"        => "man/benchmarks/benchmarks.md",
                         "Filtered oscillator"         => "man/benchmarks/filtered_oscillator.md"],

        "Further examples" => ["Overview"               => "man/examples_overview.md",
                               "Duffing oscillator"     => "models/DuffingOscillator.md",
                               "Laub-Loomis"            => "models/LaubLoomis.md",
                               "Production-Destruction" => "models/ProductionDestruction.md",
                               "Brusselator"            => "models/Brusselator.md",
                               "Quadrotor"              => "models/Quadrotor.md",
                               "Epidemic model"         => "models/SEIR.md",
                               "Spacecraft"             => "models/Spacecraft.md"],

        "API Reference" => ["Systems"                       => "lib/systems.md",
                            "Reach-sets"                    => "lib/reachsets.md",
                            "Flowpipes"                     => "lib/flowpipes.md",
                            "Linear solvers"                => LINEAR_SOLVERS_API,
                            "Nonlinear solvers"             => NONLINEAR_SOLVERS_API,
                            "Solutions"                     => "lib/solutions.md",
                            "Discretization"                => "lib/discretize.md",
                            "Projections"                   => "lib/projections.md",
                            "Clustering"                    => "lib/clustering.md",
                            "Further set operations"        => "lib/operations.md",
                            "Distributed computations"      => "lib/distributed.md",
                            "Internal functions and macros" => "lib/internals.md"],

        "Frequently Asked Questions (FAQ)" => "man/faq.md",

        "Bibliography" => "references.md",

        "How to contribute" => "about.md"
    ]
)

# Deploy built documentation from Travis.
deploydocs(
    repo = "github.com/JuliaReach/ReachabilityAnalysis.jl.git",
    push_preview = true,
)

ext: jl


text:
export getVariablesLabelsWithinRange, enableSolveAllNotDRT!

"""
    $SIGNATURES

Remove all marginalization and make solvable (=1) all variables and factors that don't contain `drt` in their label.

Notes
- Dead Reckon Tether (DRT)

Related

dontMarginalizeVariablesAll!, setSolvable!, defaultFixedLagOnTree!
"""
function enableSolveAllNotDRT!(dfg::AbstractDFG; solvable::Int=1)
  dontMarginalizeVariablesAll!(dfg)
  foreach(x->setSolvable!(dfg, x, solvable), ls(dfg))
  foreach(x->setSolvable!(dfg, x, solvable), lsf(dfg))
  nothing
end

mutable struct RangeAzimuthElevation
  range::Float64
  azimuth::Float64
  elevation::Union{Nothing,Float64}
end

function convert(::Type{Rotations.UnitQuaternion}, q::TransformUtils.Quaternion)
  Rotations.UnitQuaternion(q.s, q.v...)
end
function convert(::Type{Rotations.UnitQuaternion}, x::SO3)
  q = convert(TransformUtils.Quaternion, x)
  convert(Rotations.UnitQuaternion, q)
end
function convert(::Type{T}, x::SO3) where {T <: CoordinateTransformations.AffineMap}
  LinearMap( convert(Rotations.UnitQuaternion, x) )
end

function convert(::Type{T}, x::SE3) where {T <: CoordinateTransformations.AffineMap}
  Translation(x.t...) ∘ convert(AffineMap{Rotations.UnitQuaternion{Float64}}, x.R)
end
function convert(::Type{SE3}, x::T) where {T <: CoordinateTransformations.AffineMap{Rotations.UnitQuaternion{Float64}}}
  SE3(x.translation[1:3], TransformUtils.Quaternion(x.linear.w, [x.linear.x,x.linear.y,x.linear.z]) )
end

function convert(::Type{SE3}, t::Tuple{Symbol, Vector{Float64}})
  if t[1]==:XYZqWXYZ
    return SE3(t[2][1:3],TransformUtils.Quaternion(t[2][4],t[2][5:7]))
  else
    error("Unknown conversion type $(t[1])")
  end
end

function convert(::Type{RangeAzimuthElevation}, val::Tuple{Symbol, Vector{Float64}})
  if val[1] == :rangeazimuth
    return RangeAzimuthElevation(val[2][1],val[2][2],nothing)
  elseif val[1] == :rangeazimuthelevation
    return RangeAzimuthElevation(val[2][1],val[2][2],val[2][3])
  else
    error("Unknown conversion from $(val[1]) to RangeAzimuthElevation.")
  end
end

# should be deprecated or indicated more clearly
lsrBR(a) = [a[2,:];a[1,:]]';

function veePose3(s::SE3)
  TransformUtils.veeEuler(s)
end
function veePose(s::SE3)
  TransformUtils.veeEuler(s)
end

function \(s::SE3, wTr::CTs.Translation)
  bTr = s.R.R'*(wTr.v-s.t)
  Dtr = bTr
  range = norm(Dtr)
  azi = atan(Dtr[2], Dtr[1])
  elev = atan(Dtr[3], Dtr[1])
  RangeAzimuthElevation(range, azi, elev)
end

"""
    $SIGNATURES

Method to compare current and predicted estimate on a variable, developed for testing a new factor before adding to the factor graph.

Notes
- `fct` does not have to be in the factor graph -- likely used to test beforehand.
- function is useful for detecting if `multihypo` should be used.
- `approxConv` will project the full belief estimate through some factor but must already be in factor graph.

Example

```julia
# fg already exists containing :x7 and :l3
pp = Pose2Point2BearingRange(Normal(0,0.1),Normal(10,1.0))
# possible new measurement from :x7 to :l3
curr, pred = predictVariableByFactor(fg, :l3, pp, [:x7; :l3])
# example of naive user defined test on fit score
fitscore = minkld(curr, pred)
# `multihypo` can be used as option between existing or new variables
```

Related

approxConv
"""
function predictVariableByFactor(dfg::AbstractDFG,
                                 targetsym::Symbol,
                                 fct::Pose2Point2BearingRange,
                                 prevars::Vector{Symbol}  )
  #
  @assert targetsym in prevars
  curr = getKDE(dfg, targetsym)
  tfg = initfg()
  for var in prevars
    varnode = getVariable(dfg, var)
    addVariable!(tfg, var, getSofttype(varnode))
    if var != targetsym
      @assert isInitialized(varnode)
      initManual!(tfg,var,getKDE(varnode))
    end
  end
  addFactor!(tfg, prevars, fct, graphinit=false)
  fctsym = ls(tfg, targetsym)

  pts, infd = predictbelief(tfg, targetsym, fctsym)
  pred = manikde!(pts, getSofttype(getVariable(dfg, targetsym)))
  # return current and predicted beliefs
  return curr, pred
end

"""
    $(SIGNATURES)

Calculate the cartesian distance between two vertices in the graph using their symbol name, and by maximum belief point.
"""
function getRangeKDEMax2D(fgl::AbstractDFG, vsym1::Symbol, vsym2::Symbol)
  x1 = getKDEMax(getBelief(fgl, vsym1))
  x2 = getKDEMax(getBelief(fgl, vsym2))
  norm(x1[1:2]-x2[1:2])
end

function measureMeanDist(fg::AbstractDFG, a::AbstractString, b::AbstractString)
    #bearrang!(residual::Array{Float64,1}, Z::Array{Float64,1}, X::Array{Float64,1}, L::Array{Float64,1})
    res = zeros(2)
    A = getVal(fg,a)
    B = getVal(fg,b)
    Ax = Statistics.mean(vec(A[1,:]))
    Ay = Statistics.mean(vec(A[2,:]))
    Bx = Statistics.mean(vec(B[1,:]))
    By = Statistics.mean(vec(B[2,:]))
    dx = Bx - Ax
    dy = By - Ay
    b = atan(dy,dx)
    r = sqrt(dx^2 + dy^2)
    return r, b
end

function predictBodyBR(fg::AbstractDFG, a::Symbol, b::Symbol)
  res = zeros(2)
  A = getKDEMean(getKDE(getVariable(fg,a)))
  B = getKDEMean(getKDE(getVariable(fg,b)))
  Ax = A[1] # Statistics.mean(vec(A[1,:]))
  Ay = A[2] # Statistics.mean(vec(A[2,:]))
  Ath = getKDEMax(getKDE(getVariable(fg,a)))[3]
  Bx = B[1]
  By = B[2]
  wL = SE2([Bx;By;0.0])
  wBb = SE2([Ax;Ay;Ath])
  bL = se2vee((wBb \ Matrix{Float64}(LinearAlgebra.I, 3,3)) * wL)
  dx = bL[1] - 0.0
  dy = bL[2] - 0.0
  b = (atan(dy,dx))
  r = sqrt(dx^2 + dy^2)
  return b, r
end

"""
    $SIGNATURES

Return the last `number::Int` of poses according to `filterLabel::Regex`.

Notes
- Uses FIFO add history of variables to the distribued factor graph object as search index.
"""
function getLastPoses(dfg::AbstractDFG;
                      filterLabel::Regex=r"x\d",
                      number::Int=5)::Vector{Symbol}
  #
  # filter according to pose label
  syms = filter(l->occursin(filterLabel, string(l)), getAddHistory(dfg))

  # return the last segment of syms
  len = length(syms)
  st = number < len ? len-number+1 : 1
  return syms[st:end]
end

"""
    $SIGNATURES

Set old poses and adjacent factors to `solvable::Int=0` (default).

Notes
- `youngest::Int` and `oldest::Int` set the limits of search by count,
  - `oldest` set large enough for solver loop defintely disengage old parts in re-occuring cycle.
- `filterLabel::Regex` sets the template to search for each pose label.
- Poses are assumed to be a thread through time that connects the local exploration variables.
- Initially developed to remove old variables and factors from a solution, in combination with fix-lag marginalization.
  - `getSolverParams(fg).isfixedlag=true`

Related:

getLastPoses
"""
function setSolvableOldPoses!(dfg::AbstractDFG;
                              youngest::Int=50,
                              oldest::Int=200,
                              solvable=0,
                              filterLabel::Regex=r"x\d")
  #
  # collect old variables and factors to disable from next solve
  newPoses = getLastPoses(dfg,filterLabel=filterLabel, number=youngest)
  oldPoses = setdiff(getLastPoses(dfg,filterLabel=filterLabel, number=oldest), newPoses)
  allFcts = (oldPoses .|> x->ls(dfg,x))
  fctAdj = 0 < length(allFcts) ? union(allFcts...) : Symbol[]

  # all together
  fullList = [oldPoses; fctAdj]

  # use solvable=0 to disable variables and factors in the next solve
  map(x->setSolvable!(dfg, x, solvable), fullList)

  return fullList
end

"""
    $(SIGNATURES)

Create a new variable node and insert odometry constraint factor between
which will automatically increment latest pose symbol x<k+1> for new node new node and
constraint factor are returned as a tuple.
"""
function addOdoFG!(fg::G,
                   n::Symbol,
                   DX::Array{Float64,1},
                   cov::Array{Float64,2};
                   N::Int=0,
                   solvable::Int=1,
                   labels::Vector{<:AbstractString}=String[]  ) where G <: AbstractDFG
    #
    prev, X, nextn = getLastPose2D(fg)
    r,c = size(X)
    if N==0
      N = c
    end
    sig = diag(cov)
    XnextInit = zeros(r,c)
    # increases the number of particles based on the number of modes in the measurement Z
    for i in 1:c
        ent = [randn()*sig[1]; randn()*sig[2]; randn()*sig[3]]
        XnextInit[:,i] = addPose2Pose2(X[:,i], DX + ent)
    end

    v = addVariable!(fg, n, Pose2, N=N, solvable=solvable, tags=[labels;"POSE"])
    # v = addVariable!(fg, n, XnextInit, cov, N=N, solvable=solvable, tags=labels)
    pp = Pose2Pose2(MvNormal(DX, cov)) #[prev;v],
    f = addFactor!(fg, [prev;v], pp, solvable=solvable, graphinit=true )
    infor = inv(cov^2)
    # addOdoRemote(prev.index,v.index,DX,infor) # this is for remote factor graph ref parametric solution -- skipped internally by global flag variable
    return v, f
end

function addOdoFG!(fgl::G,
                   Z::Pose3Pose3;
                   N::Int=0,
                   solvable::Int=1,
                   labels::Vector{<:AbstractString}=String[]  ) where {G <: AbstractDFG}
  #
  vprev, X, nextn = getLastPoses(fgl)[1]
  vnext = addVariable!(fgl, nextn, Pose3, solvable=solvable, tags=labels)
  fact = addFactor!(fgl, [vprev;vnext], Z, graphinit=true)

  return vnext, fact

  # error("addOdoFG!( , ::Pose3Pose3, ) not currently usable, there were breaking changes. Work in Progress")
  # addOdoFG(fg, n, DX, cov, N=N, solvable=solvable, tags=labels)
end

"""
    $(SIGNATURES)

Create a new variable node and insert odometry constraint factor between
which will automatically increment latest pose symbol x<k+1> for new node new node and
constraint factor are returned as a tuple.

"""
function addOdoFG!(
        fgl::AbstractDFG,
        odo::Pose2Pose2;
        N::Int=0,
        solvable::Int=1,
        labels::Vector{<:AbstractString}=String[] ) # where {PP <: RoME.BetweenPoses}
    #
    vprev, X, nextn = getLastPose(fgl)
    if N==0
      N = size(X,2)
    end
    # vnext = addVariable!(fgl, nextn, X⊕odo, ones(1,1), N=N, solvable=solvable, tags=labels)
    vnext = addVariable!(fgl, nextn, Pose2, N=N, solvable=solvable, tags=labels)
    fact = addFactor!(fgl, [vprev;vnext], odo, graphinit=true)

    return vnext, fact
end

"""
    $(SIGNATURES)

Initialize a factor graph object as Pose2, Pose3, or neither and returns variable and factor symbols as array.
"""
function initFactorGraph!(fg::AbstractDFG;
                          P0::Union{Array{Float64,2},Nothing}=nothing,
                          init::Union{Vector{Float64},Nothing}=nothing,
                          N::Int=100,
                          lbl::Symbol=:x0,
                          solvable::Int=1,
                          firstPoseType=Pose2,
                          labels::Vector{Symbol}=Symbol[])
  #
  nodesymbols = Symbol[]
  if firstPoseType == Pose2
      init = init!=nothing ? init : zeros(3)
      P0 = P0!=nothing ? P0 : Matrix(Diagonal([0.03;0.03;0.001]))
      # init = vectoarr2(init)
      addVariable!(fg,lbl,Pose2,N=N, solvable=solvable, tags=labels )
      push!(nodesymbols, lbl)
      # v1 = addVariable!(fg, lbl, init, P0, N=N, solvable=solvable, tags=labels)
      fctVert = addFactor!(fg, [lbl;], PriorPose2(MvNormal(init, P0)), solvable=solvable, tags=labels) #[v1],
      push!(nodesymbols, Symbol(fctVert.label))
  end
  if firstPoseType == Pose3
      init = init!=nothing ? init : zeros(6)
      P0 = P0!=nothing ? P0 : Matrix(Diagonal([0.03;0.03;0.03;0.001;0.001;0.001]))
      addVariable!(fg,lbl,Pose2,N=N,solvable=solvable,tags=labels )
      push!(nodesymbols, lbl)
      # v1 = addVariable!(fg, lbl, init, P0, N=N, solvable=solvable, tags=labels)
      fctVert = addFactor!(fg, [lbl;], PriorPose3(MvNormal(init, P0)), solvable=solvable, tags=labels) #[v1],
      push!(nodesymbols, Symbol(fctVert.label))
  end
  return nodesymbols
end

function malahanobisBR(measA, preA, cov::Array{Float64,2})
    # measure landmark with noise
    res = measA - preA
    mala2 = Union{}
    #Malahanobis distance
    if false
      lambda = cov \ Matrix{Float64}(LinearAlgebra.I, 2,2)
      mala2 = res' * lambda * res
    else
      mala2 = res' * (cov \ res)
    end

    mala = sqrt(mala2)
    return mala
end

# ------------------------------------
# Transfered from IncrementalInference

function get2DSamples(fg::G; #::Union{Symbol, S};
                      from::Int=0, to::Int=999999999,
                      minnei::Int=0,
                      regexKey::Regex=r"x") where {G <: AbstractDFG, S <: AbstractString}
  #
  X = Array{Float64,1}()
  Y = Array{Float64,1}()

  # if sym = 'l', ignore single measurement landmarks
  allids = listVariables(fg, regexKey)  # fg.IDs
  saids = DFG.sortDFG(allids)
  for id in saids
    # vertlbl = string(id[1])
    vertlbl = string(id)
    # if vertlbl[1] == sym
      val = parse(Int,split(vertlbl[2:end],'_')[1])
      if from <= val && val <= to
        # if length( getOutNeighbors(fg, vertlbl[2] , needdata=true ) ) >= minnei
        if length( DFG.getNeighbors(fg, id ) ) >= minnei
          # if length(out_neighbors(fg.v[id[2]],fg.g)) >= minnei
          X=[X; vec(getVal(fg,id)[1,:]) ]
          Y=[Y; vec(getVal(fg,id)[2,:]) ]
        end
      end
    # end
  end
  return X,Y
end

"""
    $SIGNATURES

List all variables that fall in numerical range `from`, `to`, and with prefix key as specified.

Related

DFG.getVariableLabelNumber, DFT.findFactorsBetweenNaive
"""
function listVariablesLabelsWithinRange(fg::AbstractDFG,
                                        regexKey::Regex=r"x";
                                        from::Int=0, to::Int=9999999999,
                                        minnei::Int=0)
  #

  # if sym = 'l', ignore single measurement landmarks
  allids = listVariables(fg, regexKey)  # fg.IDs
  saids = DFG.sortDFG(allids)
  mask = Array{Bool,1}(undef, length(saids))
  fill!(mask, false)
  count = 0
  for id in saids
    count += 1
    if length( DFG.getNeighbors(fg, id) ) >= minnei
      mask[count] = true
    end
    if occursin(regexKey, string(id)) && (from != 0 || to != 9999999999)
      vertlbl = string(id)
        # TODO won't work with nested labels
        m = match(r"\d+", string(id))
        val = parse(Int, m.match)
        # val_ = split(vertlbl[2:end],'_')[1]
        # val = parse(Int,val_)
        if !(from <= val <= to)
          mask[count] = false
        end
    end
  end

  saids[mask]
end

@deprecate getVariablesLabelsWithinRange(fg::AbstractDFG,regexKey::Regex=r"x";from::Int=0, to::Int=9999999999,minnei::Int=0) listVariablesLabelsWithinRange(fg,regexKey,from=from,to=to,minnei=minnei)

function get2DSampleMeans(fg::AbstractDFG,
                          regexKey::Regex=r"x";
                          from::Int=0, to::Int=9999999999,
                          minnei::Int=0)
  #
  X = Array{Float64,1}()
  Y = Array{Float64,1}()
  Th = Array{Float64,1}()
  LB = String[]

  vsyms = getVariablesLabelsWithinRange(fg, regexKey, from=from, to=to, minnei=minnei)

  for id in vsyms
    X=[X; Statistics.mean( vec( getVal(fg, id )[1,:] ) )]
    Y=[Y; Statistics.mean( vec( getVal(fg, id )[2,:] ) )]
    # crude test for pose TODO probably not going to always work right
    if string(id)[1] == 'x'
      Th=[Th; Statistics.mean( vec( getVal(fg, id )[3,:] ) )]
    end
    push!(LB, string(id))
  end
  return X,Y,Th,LB
end

#draw landmark positions
function getAll2DMeans(fg, sym::Regex)
  return get2DSampleMeans(fg, sym )
end

function getAll2DPoses(fg::G; regexKey=r"x") where G <: AbstractDFG
    return getAll2DSamples(fg, regexKey=regexKey )
end

function get2DPoseSamples(fg::G; from::Int=0, to::Int=999999999, regexKey=r"x") where G <: AbstractDFG
  return get2DSamples(fg, regexKey=regexKey, from=from, to=to )
end

function get2DPoseMeans(fg::G; from::Int=0, to::Int=999999999, regexKey=r"x") where G <: AbstractDFG
  return get2DSampleMeans(fg, regexKey, from=from, to=to )
end

function get2DPoseMax(fgl::G;
					  regexKey::Regex=r"x",
                      from::Int=-99999999999, to::Int=9999999999 ) where G <: AbstractDFG
  #
  # xLB,ll = ls(fgl) # TODO add: from, to, special option 'x'
  xLB = DFG.getVariableIds(fgl, regexKey)
  saids = DFG.sortDFG(xLB)
  X = Array{Float64,1}()
  Y = Array{Float64,1}()
  Th = Array{Float64,1}()
  LB = String[]
  for slbl in saids
    lbl = string(slbl)
    if from <= parse(Int,split(lbl[2:end],'_')[1]) <=to
      mv = getKDEMax(getBelief(fgl,slbl))
      push!(X,mv[1])
      push!(Y,mv[2])
      push!(Th,mv[3])
      push!(LB, string(lbl))
    end
  end
  return X, Y, Th, LB
end

function get2DLandmSamples(fg::G;
                           from::Int=0,
                           to::Int=999999999,
                           minnei::Int=0 ) where G <: AbstractDFG
  #
  return get2DSamples(fg, regexKey=r"l", from=from, to=to, minnei=minnei )
end

function get2DLandmMeans(fg::G;
                         from::Int=0, to::Int=999999999,
                         minnei::Int=0,
                         landmarkRegex::Regex=r"l" ) where G <: AbstractDFG
  #
  return get2DSampleMeans(fg, landmarkRegex, from=from, to=to, minnei=minnei )
end

function removeKeysFromArr(fgl::G,
                           torm::Array{Int,1},
                           lbl::Array{String,1}) where G <: AbstractDFG
  #
  retlbs = String[]
  for i in 1:length(lbl)
    id = parse(Int,split(lbl[i][2:end],'_')[1])
    if something(findfirst(isequal(id), torm), 0) == 0 #findfirst(torm,id) == 0
      push!(retlbs, lbl[i])
    else
      println("removeKeysFromArr -- skipping $(lbl[i]), id=$(id)")
    end
  end
  return retlbs
end
function removeKeysFromArr(fgl::G,
                           torm::Array{Int,1},
                           lbl::Array{Symbol,1} ) where G
  #
  removeKeysFromArr(fgl, torm, string.(lbl))
end

function get2DLandmMax(fgl::G;
                       from::Int=-99999999999,
                       to::Int=9999999999,
                       showmm=false, MM::Dict{Int,T}=Dict{Int,Int}(),
                       regexLandmark::Regex=r"l"  ) where {G <: AbstractDFG, T}
  #
  # xLB,lLB = ls(fgl) # TODO add: from, to, special option 'x'
  lLB = DFG.getVariableIds(fgl, regexLandmark)
  if !showmm lLB = removeKeysFromArr(fgl, collect(keys(MM)), lLB); end
  X = Array{Float64,1}()
  Y = Array{Float64,1}()
  Th = Array{Float64,1}()
  LB = String[]
  for lb in lLB
    @show lb
    lbl = string(lb)
    if from <= parse(Int,split(lbl[2:end],'_')[1]) <=to
      mv = getKDEMax(getBelief(fgl, Symbol(lb)))
      push!(X,mv[1])
      push!(Y,mv[2])
      push!(LB, string(lbl))
    end
  end
  return X, Y, Th, LB
end

# convenience function to add DIDSON sonar constraints to graph
function addLinearArrayConstraint(fgl::G,
                                  rangebearing::Union{Tuple{Float64, Float64}, Vector{Float64}},
                                  pose::Symbol,
                                  landm::Symbol ;
                                  rangecov::Float64=3e-4,
                                  bearingcov::Float64=3e-4 ) where G <: AbstractDFG

  #

  cl = LinearRangeBearingElevation((rangebearing[1],rangecov),(rangebearing[2],bearingcov))
  if !haskey(fgl.IDs, landm)
    pts = getVal(fgl, pose) + cl
    N = size(pts,2)
    vl1 = addVariable!(fgl, landm,  pts,  N=N)
    println("Automatically added $(landm) to the factor graph")
  end
  addFactor!(fgl, [getVert(fgl, pose); getVert(fgl, landm)], cl)
  nothing
end

function addSoftEqualityPoint2D(fgl::G,
                                l1::Symbol,
                                l2::Symbol;
                                dist=MvNormal([0.0;0.0],Matrix{Float64}(LinearAlgebra.I, 2,2)),
                                solvable::Int=1  )  where G <: AbstractDFG
  #
  pp = Point2DPoint2D(dist)
  addFactor!(fgl, [l1,l2], pp, solvable=solvable)
end


text:
let
    mpsim = MPsim()

    attrib = Attribute("Attr A")
    addPossibleAttributeValue!( attrib, "Val AA", "Val AB" )
    addSimulationAttribute!( mpsim, attrib )

    node = BaseNode("Test A")
    addNodeRequirement!( node, ("Attr A", "Val AA") )
    addSimulationBaseNode!( mpsim, node )

    recr = Recruitment("EW")
    setRecruitmentTarget!( recr, "Test A" )
    setRecruitmentSchedule!( recr, 12 )
    setRecruitmentFixed!( recr, 10 )
    setRecruitmentAgeFixed!( recr, 18*12 )
    addSimulationRecruitment!( mpsim, recr )

    attrit = Attrition()
    setAttritionRate!( attrit, .5 )
    setAttritionPeriod!( attrit, 12 )
    addSimulationAttrition!( mpsim, attrit )

    setSimulationLength!( mpsim, 24 )
    run( mpsim )
end


text:
function odd(n::Integer)
  n & 1 == 1
end

function odd(n::Real)
  int(n) & 1 == 1
end

# Given integers a, n, and m with n ≥ 0 and 0 ≤ a < m, compute mod(a^n, m)
# using the fast exponentiation method.
# Based on
#   http://homepages.math.uic.edu/~leon/cs-mcs401-s08/handouts/fastexp.pdf
#
function fast_exponentiation(a, n, m)

  x = a
  y = odd(n) ? a : 1
  nprim = floor(n / 2)

  while nprim > 0
    x = mod(x .* x, m)

    if odd(nprim)
      y = (y == 1) ? x : mod(x .* y, m)
    end

    nprim = floor(nprim / 2)
  end

  return y

end


text:
digits_prec(prec::Integer) = floor(Int, prec * log(2) / log(10))

Base.show(io::IO, x::Union{Arb, Acb}) = print(io, string_nice(x))
Base.show(io::IO, x::Arf) = print(io, string_decimal(x))

for ArbT in (Arf, Arb, Acb, Mag)
    arbf = Symbol(cprefix(ArbT), :_, :print)
    @eval begin
        function _string(x::$ArbT)
            Libc.flush_cstdio()
            Base.flush(stdout)
            original_stdout = stdout
            out_rd, out_wr = redirect_stdout()
            try
                ccall(@libarb($arbf), Cvoid, (Ref{$ArbT},), x)
                Libc.flush_cstdio()
            finally
                redirect_stdout(original_stdout)
                close(out_wr)
            end
            return read(out_rd, String)
        end
    end

    ArbT == Mag && continue # no mag_printd and mag_printn

    arbf = Symbol(cprefix(ArbT), :_printd)
    @eval begin
        function string_decimal(x::$ArbT, digits::Integer=digits_prec(precision(x)))
            Libc.flush_cstdio()
            Base.flush(stdout)
            original_stdout = stdout
            out_rd, out_wr = redirect_stdout()
            try
                ccall(@libarb($arbf), Cvoid, (Ref{$ArbT}, Clong), x, digits)
                Libc.flush_cstdio()
            finally
                redirect_stdout(original_stdout)
                close(out_wr)
            end
            return read(out_rd, String)
        end
    end

    ArbT == Arf && continue #no arf_printn

    arbf = Symbol(cprefix(ArbT), :_, :printn)
    @eval begin
        function string_nice(x::$ArbT, digits::Integer=digits_prec(precision(x)), flags::UInt=UInt(0))
            Libc.flush_cstdio()
            Base.flush(stdout)
            original_stdout = stdout
            out_rd, out_wr = redirect_stdout()
            try
                ccall(@libarb($arbf), Cvoid, (Ref{$ArbT}, Clong, Culong), x, digits, flags)
                Libc.flush_cstdio()
            finally
                redirect_stdout(original_stdout)
                close(out_wr)
            end
            return read(out_rd, String)
        end
    end
end

for ArbT in (Arf, Arb, Mag)
    jlf = :_load_string!
    arbf = Symbol(cprefix(ArbT), :_load_str)
    @eval begin
        function $(jlf)(x::$ArbT, str::AbstractString)
            res = ccall(@libarb($arbf), Cint, (Ref{$ArbT}, Cstring), x, str)
            iszero(res) || throw(ArgumentError("arblib could not load_str $str as "*$(string(ArbT))))
            return x
        end
    end

    jlf = :_dump_string
    arbf = Symbol(cprefix(ArbT), :_dump_str)
    @eval begin
        function $(jlf)(x::$ArbT)
            char_ptr = ccall(@libarb($arbf), Cstring, (Ref{$ArbT},), x)
            str = unsafe_string(char_ptr)
            ccall((:flint_free, libflint), Cvoid, (Cstring,), char_ptr)
            return str
        end
    end
end


text:
module TestExamples1DHypDiff

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "1d")

@testset "Hyperbolic diffusion" begin

  @testset "elixir_hypdiff_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic.jl"),
      l2   = [1.3655114954641076e-7, 1.0200345025539218e-6],
      linf = [7.173285538342178e-7, 4.507116681651269e-6])
  end

  @testset "elixir_hypdiff_harmonic_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_harmonic_nonperiodic.jl"),
      l2   = [3.0130941075207524e-12, 2.6240829677090014e-12],
      linf = [4.054534485931072e-12, 3.8826719617190975e-12])
  end
end

end # module


text:
## preprocess.jl - preprocess an expression for better parsing

const PREPROCESS_RULES = [
    :(mean(abs2, _)) => :(mean(abs2.(_))),
]

function preprocess(ex)
    return rewrite_all(ex, PREPROCESS_RULES)
end