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

@test integrate(x^2*cos(a+b*x+c*x^2), x) == :((1/2)*x*c^-1*sin(a+b*x+c*x^2)+-1//4*b*c^-2*sin(a+b*x+c*x^2)+-1//4*2^(1/2)*pi^(1/2)*c^-3//2*cos(a+-1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+-1//4*2^(1/2)*pi^(1/2)*c^-3//2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1)+-1//8*2^(1/2)*pi^(1/2)*b^2*c^-5//2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1)+1//8*2^(1/2)*pi^(1/2)*b^2*c^-5//2*cos(a+-1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x)))
@test integrate(x*cos(a+b*x+c*x^2), x) == :((1/2)*c^-1*sin(a+b*x+c*x^2)+-1//4*b*2^(1/2)*pi^(1/2)*c^-3//2*cos(a+-1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+1//4*b*2^(1/2)*pi^(1/2)*c^-3//2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate(cos(a+b*x+c*x^2), x) == :((1/2)*2^(1/2)*pi^(1/2)*c^-1//2*cos(a+-1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+-1//2*2^(1/2)*pi^(1/2)*c^-1//2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate(x^-2*cos(a+b*x+c*x^2)+b*x^-1*sin(a+b*x+c*x^2), x) == :(-1*x^-1*cos(a+b*x+c*x^2)+-1*2^(1/2)*pi^(1/2)*c^(1/2)*cos(a+-1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+-1*2^(1/2)*pi^(1/2)*c^(1/2)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1))
@test integrate(x^2*cos(a+b*x+-1*c*x^2), x) == :(-1//2*x*c^-1*sin(a+b*x+-1*c*x^2)+-1//4*b*c^-2*sin(a+b*x+-1*c*x^2)+-1//4*2^(1/2)*pi^(1/2)*c^-3//2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1)+1//4*2^(1/2)*pi^(1/2)*c^-3//2*cos(a+1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))+-1//8*2^(1/2)*pi^(1/2)*b^2*c^-5//2*cos(a+1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))+-1//8*2^(1/2)*pi^(1/2)*b^2*c^-5//2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1))
@test integrate(x*cos(a+b*x+-1*c*x^2), x) == :(-1//2*c^-1*sin(a+b*x+-1*c*x^2)+-1//4*b*2^(1/2)*pi^(1/2)*c^-3//2*cos(a+1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))+-1//4*b*2^(1/2)*pi^(1/2)*c^-3//2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1))
@test integrate(cos(a+b*x+-1*c*x^2), x) == :(-1//2*2^(1/2)*pi^(1/2)*c^-1//2*cos(a+1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))+-1//2*2^(1/2)*pi^(1/2)*c^-1//2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1))
@test integrate(x^-2*cos(a+b*x+-1*c*x^2)+b*x^-1*sin(a+b*x+-1*c*x^2), x) == :(-1*x^-1*cos(a+b*x+-1*c*x^2)+2^(1/2)*pi^(1/2)*c^(1/2)*cos(a+1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))+-1*2^(1/2)*pi^(1/2)*c^(1/2)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+-2*c*x))*sin(a+1//4*b^2*c^-1))
@test integrate(x^2*cos(1//4+x+x^2), x) == :(-1//4*sin(1//4+x+x^2)+(1/2)*x*sin(1//4+x+x^2)+-1//4*2^(1/2)*pi^(1/2)*fresnels((1/2)*2^(1/2)*pi^-1//2*(1+2x))+1//8*2^(1/2)*pi^(1/2)*fresnelc((1/2)*2^(1/2)*pi^-1//2*(1+2x)))
@test integrate(x*cos(1//4+x+x^2), x) == :((1/2)*sin(1//4+x+x^2)+-1//4*2^(1/2)*pi^(1/2)*fresnelc((1/2)*2^(1/2)*pi^-1//2*(1+2x)))
@test integrate(cos(1//4+x+x^2), x) == :((1/2)*2^(1/2)*pi^(1/2)*fresnelc((1/2)*2^(1/2)*pi^-1//2*(1+2x)))
@test integrate(x^2*cos(a+b*x+c*x^2)^2, x) == :(1//6*x^3+-1//16*b*c^-2*sin(2a+2*b*x+2*c*x^2)+1//8*x*c^-1*sin(2a+2*b*x+2*c*x^2)+-1//16*pi^(1/2)*c^-3//2*cos(2a+-1//2*b^2*c^-1)*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))+-1//16*pi^(1/2)*c^-3//2*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1)+-1//16*pi^(1/2)*b^2*c^-5//2*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1)+1//16*pi^(1/2)*b^2*c^-5//2*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x)))
@test integrate(x*cos(a+b*x+c*x^2)^2, x) == :(1//4*x^2+1//8*c^-1*sin(2a+2*b*x+2*c*x^2)+-1//8*b*pi^(1/2)*c^-3//2*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))+1//8*b*pi^(1/2)*c^-3//2*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1))
@test integrate(cos(a+b*x+c*x^2)^2, x) == :((1/2)*x+-1//4*pi^(1/2)*c^-1//2*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1)+1//4*pi^(1/2)*c^-1//2*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x)))
@test integrate(x^2*cos(a+b*x+-1*c*x^2)^2, x) == :(1//6*x^3+-1//8*x*c^-1*sin(2a+-2*c*x^2+2*b*x)+-1//16*b*c^-2*sin(2a+-2*c*x^2+2*b*x)+-1//16*pi^(1/2)*c^-3//2*fresnelc(pi^-1//2*c^-1//2*(b+-2*c*x))*sin(2a+(1/2)*b^2*c^-1)+1//16*pi^(1/2)*c^-3//2*cos(2a+(1/2)*b^2*c^-1)*fresnels(pi^-1//2*c^-1//2*(b+-2*c*x))+-1//16*pi^(1/2)*b^2*c^-5//2*cos(2a+(1/2)*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+-2*c*x))+-1//16*pi^(1/2)*b^2*c^-5//2*fresnels(pi^-1//2*c^-1//2*(b+-2*c*x))*sin(2a+(1/2)*b^2*c^-1))
@test integrate(x*cos(a+b*x+-1*c*x^2)^2, x) == :(1//4*x^2+-1//8*c^-1*sin(2a+-2*c*x^2+2*b*x)+-1//8*b*pi^(1/2)*c^-3//2*cos(2a+(1/2)*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+-2*c*x))+-1//8*b*pi^(1/2)*c^-3//2*fresnels(pi^-1//2*c^-1//2*(b+-2*c*x))*sin(2a+(1/2)*b^2*c^-1))
@test integrate(cos(a+b*x+-1*c*x^2)^2, x) == :((1/2)*x+-1//4*pi^(1/2)*c^-1//2*cos(2a+(1/2)*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+-2*c*x))+-1//4*pi^(1/2)*c^-1//2*fresnels(pi^-1//2*c^-1//2*(b+-2*c*x))*sin(2a+(1/2)*b^2*c^-1))
@test integrate(x^2*cos(1//4+x+x^2)^2, x) == :(-1//16*sin(1/2+2x+2*x^2)+1//6*x^3+-1//16*pi^(1/2)*fresnels(pi^-1//2*(1+2x))+1//8*x*sin(1/2+2x+2*x^2)+1//16*pi^(1/2)*fresnelc(pi^-1//2*(1+2x)))
@test integrate(x*cos(1//4+x+x^2)^2, x) == :(1//4*x^2+1//8*sin(1/2+2x+2*x^2)+-1//8*pi^(1/2)*fresnelc(pi^-1//2*(1+2x)))
@test integrate(cos(1//4+x+x^2)^2, x) == :((1/2)*x+1//4*pi^(1/2)*fresnelc(pi^-1//2*(1+2x)))
@test integrate((d+e*x)^2*cos(a+b*x+c*x^2), x) == :((1/2)*e*c^-1*(d+e*x)*sin(a+b*x+c*x^2)+1//4*e*c^-2*(-1*b*e+2*c*d)*sin(a+b*x+c*x^2)+-1//4*2^(1/2)*pi^(1/2)*c^-3//2*e^2*cos(a+-1//4*b^2*c^-1)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))+-1//4*2^(1/2)*pi^(1/2)*c^-3//2*e^2*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1)+-1//8*2^(1/2)*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1)+1//8*2^(1/2)*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*cos(a+-1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x)))
@test integrate((d+e*x)*cos(a+b*x+c*x^2), x) == :((1/2)*e*c^-1*sin(a+b*x+c*x^2)+-1//4*2^(1/2)*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*fresnels((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x))*sin(a+-1//4*b^2*c^-1)+1//4*2^(1/2)*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*cos(a+-1//4*b^2*c^-1)*fresnelc((1/2)*2^(1/2)*pi^-1//2*c^-1//2*(b+2*c*x)))
@test integrate((d+e*x)^2*cos(a+b*x+c*x^2)^2, x) == :(1//6*e^-1*(d+e*x)^3+1//8*e*c^-1*(d+e*x)*sin(2a+2*b*x+2*c*x^2)+1//16*e*c^-2*(-1*b*e+2*c*d)*sin(2a+2*b*x+2*c*x^2)+-1//16*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1)+-1//16*pi^(1/2)*c^-3//2*e^2*cos(2a+-1//2*b^2*c^-1)*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))+-1//16*pi^(1/2)*c^-3//2*e^2*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1)+1//16*pi^(1/2)*c^-5//2*(-1*b*e+2*c*d)^2*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x)))
@test integrate(cos(a+b*x+c*x^2)^2*(d+e*x), x) == :(1//4*e^-1*(d+e*x)^2+1//8*e*c^-1*sin(2a+2*b*x+2*c*x^2)+-1//8*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*fresnels(pi^-1//2*c^-1//2*(b+2*c*x))*sin(2a+-1//2*b^2*c^-1)+1//8*pi^(1/2)*c^-3//2*(-1*b*e+2*c*d)*cos(2a+-1//2*b^2*c^-1)*fresnelc(pi^-1//2*c^-1//2*(b+2*c*x)))

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.13 (d+e x)^m cos(a+b x+c x^2)^n.jl


text:
function model(model_name, model_expr, mod = Base)
    model_object = read_model(model_expr, mod)
    model_object_name = Symbol(model_name, "#OBJECT#DEFINITION#")
    quote
        struct $model_name{T <: NamedTuple}
            data::T
        end
        const $model_object_name = $model_object
        @generated function ProbabilityModels.PaddedMatrices.logdensity(
            var"###DATA###"::$model_name{var"##PROB#MODEL#NAMED#TUPLE#TYPE##"}, var"#θ#"
        ) where {var"##PROB#MODEL#NAMED#TUPLE#TYPE##"}
            var"#__m__#" = deepcopy($model_object_name)
            ProbabilityModels.preprocess!(var"#__m__#", var"##PROB#MODEL#NAMED#TUPLE#TYPE##")
            q = ProbabilityModels.ReverseDiffExpressionsBase.lower(var"#__m__#")
            pushfirst!(q.args, Expr(:(=), var"#DATA#", Expr(:(.), var"###DATA###", QuoteNode(:data))))
            q
        end
        @generated function ProbabilityModels.PaddedMatrices.∂logdensity!(
            var"#∇#", var"###DATA###"::$model_name{var"##PROB#MODEL#NAMED#TUPLE#TYPE##"}, var"#θ#"
        ) where {var"##PROB#MODEL#NAMED#TUPLE#TYPE##"}
            var"#__m__#" = deepcopy($model_object_name)
            ProbabilityModels.preprocess!(var"#__m__#", var"##PROB#MODEL#NAMED#TUPLE#TYPE##")
            var"#__∂m__#" = ProbabilityModels.ReverseDiffExpressionsBase.differentiate(var"#__m__#")
            q = ProbabilityModels.ReverseDiffExpressionsBase.lower(var"#__∂m__#")
            pushfirst!(q.args, Expr(:(=), var"#DATA#", Expr(:(.), var"###DATA###", QuoteNode(:data))))
            q
        end
    end
end

macro model(model_name, model_expr)
    model(model_name, model_expr, __module__) |> esc
end

ext: jl


text:
"Ground action definition."
struct GroundAction <: Action
    name::Symbol
    term::Compound
    preconds::Vector{Term}
    effect::Union{GenericDiff,ConditionalDiff}
end

get_name(action::GroundAction) = action.name

get_argvals(action::GroundAction) = term.args

get_precond(action::GroundAction) = Compound(:and, action.preconds)

get_effect(action::GroundAction) = as_term(action.effect)

"Group of ground actions with a shared schema."
struct GroundActionGroup <: Action
    name::Symbol
    actions::Dict{Compound,GroundAction}
end

function GroundActionGroup(name::Symbol, actions::AbstractVector{GroundAction})
    actions = Dict(act.term => act for act in actions)
    return GroundActionGroup(name, actions)
end

get_name(action::GroundActionGroup) = action.name

"Returns an iterator over all ground arguments of an `action`."
function groundargs(domain::Domain, state::State, action::Action)
    iters = (get_objects(domain, state, ty) for ty in get_argtypes(action))
    return Iterators.product(iters...)
end

"""
    groundactions(domain::Domain, state::State, action::Action)

Returns ground actions for a lifted `action` in a `domain` and initial `state`.
"""
function groundactions(domain::Domain, state::State, action::Action,
                       statics=infer_static_fluents(domain))
    ground_acts = GroundAction[]
    act_name = get_name(action)
    act_vars = get_argvars(action)
    # Dequantify and flatten preconditions and effects
    _precond = dequantify(get_precond(action), domain, state)
    _effect = dequantify(get_effect(action), domain, state)
    _conds, _effects = flatten_conditions(_effect)
    for args in groundargs(domain, state, action)
        subst = Subst(var => val for (var, val) in zip(act_vars, args))
        term = Compound(act_name, collect(args))
        # Substitute and simplify precondition
        precond = substitute(_precond, subst)
        precond = simplify_statics(precond, domain, state, statics)
        precond.name == false && continue # Skip unsatisfiable actions
        preconds = to_cnf_clauses(precond)
        # Simplify conditions of conditional effects
        conds = map(_conds) do cs
            isempty(cs) && return cs
            cond = substitute(Compound(:and, cs), subst)
            cond = simplify_statics(cond, domain, state, statics)
            return to_cnf_clauses(cond)
        end
        # Construct diffs from conditional effect terms
        diffs = map(_effects) do es
            eff = substitute(Compound(:and, es), subst)
            return effect_diff(domain, state, eff)
        end
        # Delete unsatisfiable branches
        unsat_idxs = findall(cs -> Const(false) in cs, conds)
        deleteat!(conds, unsat_idxs)
        deleteat!(diffs, unsat_idxs)
        # Construct conditional diff if necessary
        if length(diffs) > 1
            effect = ConditionalDiff(conds, diffs)
        elseif length(diffs) == 1
            effect = diffs[1]
            preconds = append!(preconds, conds[1])
            filter!(c -> c != Const(true), preconds)
        else # Skip actions with no active branches
            continue
        end
        # Construct ground action
        act = GroundAction(act_name, term, preconds, effect)
        push!(ground_acts, act)
    end
    return ground_acts
end

function groundactions(domain::Domain, state::State, action::GroundActionGroup,
                       statics=nothing)
    return values(action.actions)
end

"""
    groundactions(domain::Domain, state::State)

Returns all ground actions for a `domain` and initial `state`.
"""
function groundactions(domain::Domain, state::State)
    statics = infer_static_fluents(domain)
    iters = (groundactions(domain, state, act, statics)
             for act in values(get_actions(domain)))
    return collect(Iterators.flatten(iters))
end

"""
    ground(domain::Domain, state::State, action::Action)

Grounds a lifted `action` in a `domain` and initial `state`, returning a
group of grounded actions.
"""
function ground(domain::Domain, state::State, action::GenericAction,
                statics=infer_static_fluents(domain))
    ground_acts = groundactions(domain, state, action, statics)
    return GroundActionGroup(action.name, ground_acts)
end


text:
#
#
#
abstract type MultiTargetFilter{T,S} <: SequentialEstimator{T,S} end

"""
    NearestNeighborMTF(filter::SequentialEstimator)
    NearestNeighborMTF(filter_bank::Vector{SequentialEstimator})

Multi-target filter with nearest-neighbor data association. A multi-target
filter contains a bank of `SequentialEstimator` types representing
currently-tracked targets. To process a measurement, the `NearestNeighborMTF`
computes a Euclidean distance between the input measurement and the next
predicted measurement for all internal estimators. The `SequentialEstimator`
with the best prediction calls its own `process!` method on the measurement.

MultiTargetFilter types require at least one `SequentialEstimator`, but accept
any number in a Vector format at construction-time. Additional estimators can
be added to the MultiTargetFilter at any time with the `add!` method.

All internal `SequentialEstimator` must have the same time format. Attempting to
mix Discrete- and Continuous-time estimators will result in an error. Internal
estimators must also be unique.
"""
immutable NearestNeighborMTF{T<:AbstractFloat, S<:AbstractState{T}} <:
                                                          MultiTargetFilter{T,S}
    filter_bank::Vector{SequentialEstimator{T,S}}

    function NearestNeighborMTF(filter::SequentialEstimator{T,S}) where {T,S}
        new{T,S}(Vector{SequentialEstimator{T}}([filter]))
    end
end
function NearestNeighborMTF(filter_bank::Vector)
    if isempty(filter_bank)
        throw(ArgumentError("MTF requires at least one sequential filter."))
    end
    new_mtf = NearestNeighborMTF(filter_bank[1])
    for idx = 2:length(filter_bank)
        add!(new_mtf, filter_bank[idx])
    end
    return new_mtf
end

const DiscreteMTF{T} =
    MultiTargetFilter{T,UncertainDiscreteState{T}} where T
const ContinuousMTF{T} =
    MultiTargetFilter{T,UncertainContinuousState{T}} where T

const DiscreteSequentialEstimator{T} =
    SequentialEstimator{T,UncertainDiscreteState{T}} where T
const ContinuousSequentialEstimator{T} =
    SequentialEstimator{T,UncertainContinuousState{T}} where T

"""
    add!(mtf::MultiTargetFilter, filter::SequentialEstimator)
    add!(mtf::MultiTargetFilter, filter_bank::Array)

Add one or more `SequentialEstimator` to a `MultiTargetFilter`. All
`SequentialEstimator` in a `MultiTargetFilter` must have the same time format.
Attempting to mix Discrete- and Continuous-time estimators will result in an
error. Internal estimators must also be unique.
"""
function add!(mtf::MultiTargetFilter, filter::SequentialEstimator)
    throw(ArgumentError("Cannot mix continuous and discrete-time filters."))
end
function add!{T,S}(mtf::MultiTargetFilter{T,S},filter::SequentialEstimator{T,S})
    for idx = 1:length(mtf.filter_bank)
        if filter === mtf.filter_bank[idx]
            throw(ArgumentError("Attempted to add duplicate filter."))
        end
    end
    push!(mtf.filter_bank, filter)
    return nothing
end
function add!(mtf::MultiTargetFilter, filter_bank::Array)
    for idx = 1:length(filter_bank)
        add!(mtf, filter_bank[idx])
    end
    return nothing
end

"""
    distance(est::SequentialEstimator, z::AbstractAbsoluteState)
"""
function distance(est::DiscreteSequentialEstimator, z::DiscreteState)
    return distance(predict(predict(est.estimate, est.sys, z.t), est.obs), z)
end
function distance(est::ContinuousSequentialEstimator, z::ContinuousState)
    return distance(predict(predict(est.estimate, est.sys, z.t), est.obs), z)
end

"""
    process!(mtf::MultiTargetFilter, z::AbstractAbsoluteState)

Process a measurement with an arbitrary `MultiTargetFilter`. Uses the
appropriate data association method for the particular type of
`MultiTargetFilter`.
"""
function process!(::MultiTargetFilter) end
function process!{T<:AbstractFloat}(mtf::NearestNeighborMTF{T},
                                    zk::AbstractAbsoluteState{T})
    distances = zeros(T, length(mtf.filter_bank))
    for idx = 1:length(mtf.filter_bank)
        @inbounds distances[idx] .= distance(mtf.filter_bank[idx], zk)
    end
    process!(mtf.filter_bank[indmin(distances)], zk)
    return nothing
end


text:
module TropicalGEMM

using LinearAlgebra, TropicalNumbers, VectorizationBase, LoopVectorization

export Tropical, TropicalF64, TropicalF32, TropicalF16

include("fallbacks.jl")
include("gemm.jl")

end


text:
# condense the tree after the child was removed from the node
# accumulate the nodes that are temporary detached and need to be reinserted in `tmpdetached`
function _condense!(node::Node, tree::RTree, tmpdetached::AbstractVector{<:Node})
    min_load = floor(Int, capacity(node, tree) * tree.fill_factor)
    @debug "_condense!() lev=$(level(node)) len=$(length(node)) tmpdetached=$(length(tmpdetached)) min_load=$min_load"

    if node === tree.root
        # eliminate root if it has only one child.
        old_root = tree.root
        min_root_level = mapreduce(level, max, tmpdetached, init=1)
        while (level(node) > min_root_level) && (length(node) == 1)
            #@debug "_condense!(): shrinking single-child root (lv=$(level(node)))"
            tree.root = node[1]
            tree.root.parent = nothing
            release(tree, node)
            pop!(tree.nnodes_perlevel) # tree became 1 level shorter
            node = tree.root
        end
        if (level(node) >= min_root_level) && isempty(node)
            #@debug "_condense!(): resetting empty root (lv=$(level(node)))"
            # reset the root to a node of minimal level to accomodate the children of tmpdetached (or leaf)
            tree.root = acquire(tree, min_root_level == 1 ? Leaf : Branch, min_root_level)
            resize!(tree.nnodes_perlevel, min_root_level)
            fill!(tree.nnodes_perlevel, 0)
            tree.nnodes_perlevel[end] = 1
        elseif node === old_root
            # root didn't change, but MBR has to be updated due to data removal.
            tree.tight_mbrs && syncmbr!(node)
        end
        return node
    else
        # find the entry in the parent, that points to this node
        node_ix = pos_in_parent(node)
        @assert node_ix !== nothing

        if length(node) < min_load
            # used space less than the minimum
            # 1. eliminate node entry from the parent. deleteEntry will fix the parent's MBR.
            _detach!(parent, node_ix, tree)
            # 2. add this node to the stack in order to reinsert its entries.
            push!(tmpdetached, node)
        else
            # global recalculation necessary since the MBR can only shrink in size,
            # due to data removal.
            tree.tight_mbrs && syncmbr(parent(node))
        end

        return _condense!(parent(node), tree, tmpdetached)
    end
end

# reinsert the *children* of the `detached` nodes back to the tree
# the `detached` node themselves are released back to the pool
function _reinsert!(tree::RTree, detached::AbstractVector{<:Node})
    @debug "_reinsert!(): reinsert children of $(length(detached)) detached nodes"
    isempty(detached) && return tree
    con = SubtreeContext(tree)
    sort!(detached, by=level, rev=true) # start with the highest nodes
    for node in detached
        for child in children(node)
            _insert!(con, child, level(node))
            tree.nnode_reinsertions += 1
        end
        release(tree, node)
    end
    empty!(detached)
    return tree
end

# update the MBR of node if required due to the child_ix MBR changing from
# child_oldmbr to its current state or if forced,
# propagate the MBR update to the higher tree levels
# return true if MBR was updated
function _updatembr!(node::Branch, child_ix::Integer, child_oldmbr::Rect,
                     con::SubtreeContext; force::Bool = false)
    @debug "_updatembr!() tree_height=$(height(con.tree)) tree_len=$(length(con.tree)) node_lev=$(level(node)) node_len=$(length(node)) force=$(force) child_ix=$(child_ix) oflow=$(isoverflow(con, level(node)))"
    child = node[child_ix]
    # MBR needs recalculation if either:
    #   1. the NEW child MBR is not contained.
    #   2. the OLD child MBR is touching.
    node_oldmbr = mbr(node)
    mbr_dirty = force || !in(mbr(child), node_oldmbr) ||
                (con.tree.tight_mbrs && touches(node_oldmbr, child_oldmbr))
    mbr_dirty && syncmbr!(node)

    if mbr_dirty && hasparent(node)
        _updatembr!(parent(node), pos_in_parent(node), node_oldmbr, con, force=force)
    end
    return mbr_dirty
end

# replace the node (with oldmbr MBR) with the newnode
# return true if node's parent MBR update was necessary
function _replace!(node::Node, newnode::Node, oldmbr::Rect,
                   con::SubtreeContext)
    @debug "_replace!() lev=$(level(node)) len=$(length(node)) newlen=$(length(newnode))"
    # find an entry pointing to the old child
    @assert level(node) == level(newnode)
    par = parent(node)
    @assert par !== nothing
    node_ix = pos_in_parent(node)
    @assert node_ix !== nothing
    par[node_ix] = newnode
    return _updatembr!(par, node_ix, oldmbr, con)
end

# replace node with n1 and n2
# return true if node's parent MBR update or any other tree
# restructurings were necessary
function _replace!(node::Node, n1::Node, n2::Node,
                   oldmbr::Rect, con::SubtreeContext)
    @debug "_replace!() lev=$(level(node)) len=$(length(node)) newlens=($(length(n1)), $(length(n2)))"
    @assert level(node) == level(n1) == level(n2)
    # find entry pointing to old node
    par = parent(node)
    @assert par !== nothing
    node_ix = pos_in_parent(node)
    @assert node_ix !== nothing

    # MBR needs recalculation if either:
    #   1. MBRs of n1 and n2 are not contained in parent
    #   2. the OLD node MBR is touching parent (FIXME should it be recalced though?)
    par_oldmbr = mbr(par)
    mbr_dirty = !in(mbr(n1), par_oldmbr) || !in(mbr(n2), par_oldmbr) ||
                (con.tree.tight_mbrs && touches(par_oldmbr, oldmbr))
    # replace the node with n1
    par[node_ix] = n1
    mbr_dirty && syncmbr!(par)

    con.tree.nnode_splits += 1
    con.tree.nnodes_perlevel[level(node)] += 1

    # No registering necessary. insert!() will write the node if needed.
    #_register!(tree, node)
    adjusted = _insert!(par, n2, con)

    # if n2 is contained in the node and there was no split or reinsert,
    # we need to adjust only if recalculation took place.
    # In all other cases insertData above took care of adjustment.
    if !adjusted && mbr_dirty && hasparent(par)
        _updatembr!(parent(par), pos_in_parent(par), par_oldmbr, con)
    end
    return adjusted || mbr_dirty
end


text:
using AlgebraicPetri
using OrdinaryDiffEq
using DiffEqFlux
#using Catalyst
using LinearAlgebra: dot

"""
We will eventually be more sophisticated and use the actual projection
maps of a pullback, but for now we check whether "I" is in the name

Returns vector int
"""
get_infected_states(g::AbstractLabelledPetriNet) =
   [i for (i,s) in enumerate(g[:sname]) if occursin("I", string(s))]

function get_bounds(g::AbstractLabelledPetriNet)
  lb, ub = Float64[],Float64[]
  for (sym,_) in g[:tname]
    sym_= string(sym)
    if occursin("inf", sym_)
      push!(lb, log(1e-8)); push!(ub, log(1e-5))
    elseif occursin("rec", sym_)
      push!(lb, log(1e-4)); push!(ub, log(1e-1))
    else
      push!(lb, log(1e-10)); push!(ub, log(1))
    end
  end
  for (sym,_) in g[:sname]
    sym_ = string(sym)
    if occursin("I", sym_)
      push!(lb, log(1e2)); push!(ub, log(1e6))
    elseif occursin("R", sym_)
      push!(lb, log(1e-10)); push!(ub, log(1))
    else
      push!(lb, log(1e2)); push!(ub, log(1e8))
    end
  end
  return lb => ub
end
"""
state - Vector{Float}
"""
sum_infected_states(g::AbstractLabelledPetriNet, states) =
  sum([s for (i, s) in enumerate(states) if i ∈ get_infected_states(g)])

"""Generate a sequence of new infections"""
function generate_data(rxn::AbstractLabelledPetriNet, rxn_rates, initial_conc)
  o = ODEProblem(ReactionSystem(rxn),  initial_conc,
                 (0.0,100.0), vcat(rxn_rates, initial_conc))
  trajectory = solve(o, Tsit5())
  t_steps = [trajectory(i) for i in 0:100]
  return [sum_infected_states(rxn, t) for t in t_steps]
end

"""
Produces function from Vector{Float} -> Float.
- The input is the logged params (exponentiate to get the real params)

Outer function inputs:
prob - current problem (ODE problem)
times - each timestep (indices of timestep correspond to indices of data array)
data - Vector{Float} representing Infected at each step - target to hit
states/rates - dictionary whatever you want to be fixed
"""
function make_loss(g::AbstractLabelledPetriNet, prob, data; states=Dict(), rates=Dict())
  function loss(p)
    cur_p = exp.(p)
    u0 = exp.(p[(1:ns(g)) .+ nt(g)])
    for (k,v) in rates
      cur_p[k] = v
    end
    for (k,v) in states
        u0[k] = v
    end

    prob′ = remake(prob, p=cur_p, u0=u0, tspan=(0.0,100.0))
    sol = solve(prob′, Tsit5())
    sol_steps = [sol(i) for i in 0:100]
    prediction = [sum_infected_states(g, s) for s in sol_steps]
    log(sum(abs2, data .- prediction)), sol
  end
end

"""
Train a Petri net against data. Return the optimized parameters and minimum loss.
"""
function train(g::AbstractLabelledPetriNet,data::Vector{Float64})::Pair{Float64, Vector{Float64}}
  param_guess = vcat(repeat([1e-10], nt(g)), repeat([1e-3], ns(g)))
  #param_guess = [1e-10,1e-10,1e6,1e6,1e6] # correct [1e-8,1e-4,1e7,1e5,1e-10]
  prob = ODEProblem(ReactionSystem(g),
                    zeros(ns(g)), # initial state
                    (0.0,100.0),
                    zeros(ns(g)+nt(g))) # current variables in model
  lb, ub = get_bounds(g)
  l_func = make_loss(g, prob, data)
  p = DiffEqFlux.sciml_train(l_func,lb,ADAM(0.1), lower_bounds=lb, upper_bounds=ub, maxiters = 10000)
  return l_func(p.u)[1] => exp.(p.u)
end

function eval_petri_fn(true_net::AbstractLabelledPetriNet, true_rates::Vector{Float64}, true_initial_conc::Vector{Float64})::Function
  data = generate_data(true_net, true_rates, true_initial_conc)
  function model_search_loss_fn(g::AbstractLabelledPetriNet)#::Float64
    return first(train(g, data))
    #return train(g, data)
  end
  return model_search_loss_fn
end

#l_func = make_loss(model, prob, week_times, week_avg; calc_inf=calc_inf, rates=Dict(5=>1/14, 6=>1/14), states=Dict(2=>1e-9, 4=>1e-9));

#p = DiffEqFlux.sciml_train(l_func,log.(param_guess),ADAM(0.1),maxiters = 1000)


text:
using QuickCheck2
using Base.Test

function is_tiny(x)
    abs(x) < 10.0^(-6)
end

function is_huge(x)
    x > 10.0^(6)
end

function is_increasing(f::Function)
            eval(:(@test_formany -Inf<x::Float64<Inf, -Inf<y::Float64<Inf, x < y --> $f(x) < $f(y)))
end

@test isa( ( @test_exists 0<x::Float64<<1, is_tiny(x) ) , Base.Test.Pass )

## borrowed from julia/test/test.jl, Fail/Error results do not throw errors when tests are performed.
type NoThrowTestSet <: Base.Test.AbstractTestSet
    results::Vector
    NoThrowTestSet(desc) = new([])
end
Base.Test.record(ts::NoThrowTestSet, t::Base.Test.Result) = (push!(ts.results, t); t)
Base.Test.finish(ts::NoThrowTestSet) = ts.results

fails = @testset NoThrowTestSet begin
	#1 Fail
	@test_formany 0<x::Float64<10, x<9 ntests = 1000
	#2 Fail
	@test_exists -Inf<x::Int64<Inf, x>10000
	#3 Fail
	is_increasing(x -> x^2)
	#4 Error: types of variables must be specified
	@test_formany Inf > x > 0, x<y<10, x>0
	#5 Pass
	@test_forall i in 1:20, 1==1
	#6 Pass
	is_increasing(x -> x)
	#7 Pass
	is_increasing(x -> x^3)
end

for i in 1:3
    @test isa(fails[i], Base.Test.Fail)
end

@test isa(fails[4],Base.Test.Error)

for i in 5:7
	@test isa(fails[i],Base.Test.Pass)
end

    #Note is_increasing is a property of the function, so what I really want
    #this to mean is "not @test is_increasing(x->x^2)", i.e.
    #"this test fails."
    #This is slightly different from the @exists idea above, because
    #(a) I'm not trying to help the test find a counterexample, but just
    #declaring that when it does I'm happy, and
    #(b) I'm obscuring where the property is. In other words to express this
    #in terms of an @exists I'd have to write
    #@exists x::Float64, y::Float64, x < y, x^2 > y^2
    #But I don't want to think about that at *all*. I just want to say
    #"x^2 is not increasing" and trust that random examples will figure this
    #out.
    #Maybe think about a syntax that makes it feasible to also specify
    #a generator form with this sort of test, if I have some idea about
    #how to guide it.


text:
@testset "DNFT" begin
    @test typeof(VisuShrink(HardTH(), 8)) == VisuShrink
    @test typeof(RelErrorShrink()) == RelErrorShrink
    @test typeof(RelErrorShrink(HardTH())) == RelErrorShrink
    @test typeof(RelErrorShrink(HardTH(), 1)) == RelErrorShrink
    x = randn(8)
    tree = maketree(x, :full)
    @test typeof(SureShrink(HardTH(), 1)) == SureShrink
    @test typeof(SureShrink(x)) == SureShrink
    @test typeof(SureShrink(x, false, tree)) == SureShrink
    @test typeof(SureShrink(x, false, tree, HardTH())) == SureShrink
end

@testset "Denoise" begin
    # Single denoising
    n = 2^8
    x0 = generatesignals(:heavisine, 8)
    x = x0 + 0.5*randn(n)
    wt = wavelet(WT.haar)
    err = relativenorm(x, x0)
    dnt = VisuShrink(2, HardTH())

    ## Non-redundant Wavelet Transforms
    y = denoise(x, :sig, wt, dnt=dnt)
    @test relativenorm(y, x0) <= err
    y = denoise(dwt(x, wt, 4), :dwt, wt, L=4, dnt=dnt, smooth=:undersmooth)
    @test relativenorm(y, x0) <= 2*err
    y = denoise(dwt(x, wt), :dwt, wt, dnt=dnt, smooth=:undersmooth)
    @test relativenorm(y, x0) <= 2*err
    y = denoise(
        wpt(x, wt), :wpt, wt, tree=maketree(x, :full), dnt=dnt,
        smooth=:undersmooth
    )
    @test relativenorm(y, x0) <= 2*err

    ## Stationary Wavelet Transforms
    y = denoise(sdwt(x, wt), :sdwt, wt, dnt=dnt, smooth=:undersmooth)
    @test relativenorm(y, x0) <= 2*err
    y = denoise(swpd(x, wt), :swpd, wt, smooth=:undersmooth)
    @test relativenorm(y, x0) <= 2*err

    ## Autocorrelation wavelet transforms
    y = denoise(acdwt(x, wt), :acwt, wt, dnt=dnt, smooth=:undersmooth)
    @test relativenorm(y, x0) <= 2*err
    y = denoise(acwpd(x, wt), :acwpt, wt, smooth=:undersmooth)
    @test relativenorm(y, x0) <= 2*err

    # Group denoising
    x0 = duplicatesignals(generatesignals(:heavisine, 8), 5, 2)
    x = duplicatesignals(generatesignals(:heavisine, 8), 5, 2, true, 0.5)
    max_err = maximum([relativenorm(x[:,i],x0[:,i]) for i in 1:5])
    dnt = VisuShrink(2, HardTH())

    ## Non-redundant Transforms
    y = denoiseall(x, :sig, wt, dnt=dnt, bestTH=mean)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err
    y = denoiseall(dwtall(x, wt), :dwt, wt, dnt=dnt)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err
    dnt = RelErrorShrink(HardTH(), 0.3)
    y = denoiseall(wptall(x, wt), :wpt, wt, tree=maketree(n, 8, :full), dnt=dnt, estnoise=relerrorthreshold)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err
    y = denoiseall(wptall(x, wt), :wpt, wt, tree=maketree(n, 8, :full), dnt=dnt, estnoise=relerrorthreshold, bestTH=mean)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err

    ## Stationary Wavelet Transforms
    y = denoiseall(cat([sdwt(x[:,i], wt) for i in 1:5]..., dims=3), :sdwt, wt)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err
    y = denoiseall(cat([sdwt(x[:,i], wt) for i in 1:5]..., dims=3), :sdwt, wt,
        dnt=dnt, estnoise=relerrorthreshold, bestTH=mean)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err
    y = denoiseall(cat([swpd(x[:,i], wt) for i in 1:5]..., dims=3), :swpd, wt,
        tree=maketree(n, 7, :full), dnt=dnt, estnoise=relerrorthreshold)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err
    y = denoiseall(cat([swpd(x[:,i], wt) for i in 1:5]..., dims=3), :swpd, wt,
        tree=maketree(n, 7, :full), dnt=dnt, estnoise=relerrorthreshold,
        bestTH=mean)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err

    ## Autocorrelation Wavelet Transforms
    y = denoiseall(cat([acdwt(x[:,i], wt) for i in 1:5]..., dims=3), :acwt, wt)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err
    y = denoiseall(cat([acwpd(x[:,i], wt) for i in 1:5]..., dims=3), :acwpt, wt,
        tree=maketree(n, 7, :full), dnt=dnt, estnoise=relerrorthreshold)
    @test mean([relativenorm(y[:,i],x0[:,i]) for i in 1:5]) <= max_err
end

@testset "Threshold Determination" begin
    # non-stationary
    x = randn(32)
    tree = maketree(32, 5, :full)
    @test typeof(noisest(x, false)) <: Real
    @test typeof(Denoising.surethreshold(x, false)) <: Real
    @test typeof(relerrorthreshold(x, false)) <: Real
    # stationary - sdwt
    x = randn(32, 6)
    @test typeof(noisest(x, true, nothing)) <: Real
    @test typeof(Denoising.surethreshold(x, true, nothing)) <: Real
    @test typeof(relerrorthreshold(x, true, nothing)) <: Real
    # stationary - swpd
    x = randn(32, 63)
    @test typeof(noisest(x, true, tree)) <: Real
    @test typeof(Denoising.surethreshold(x, true, tree)) <: Real
    @test typeof(relerrorthreshold(x, true, tree)) <: Real
    # relative error plot
    @test typeof(relerrorthreshold(x, true, tree, makeplot=true)[2]) ==
        Plots.Plot{Plots.GRBackend}
end


text:
const gbm_study_uid = "1.3.6.1.4.1.14519.5.2.1.4591.4001.304604545029494418165835320551"

function download_invivo_studies(
    study = gbm_study_uid;
    destination::AbstractString,
    overwrite = false,
)
    make_folder(destination)
    (vfa_folder, dce_folder) = download_vfa_and_dce(study; destination, overwrite)
    return (vfa_folder, dce_folder)
end

function download_vfa_and_dce(study_id; destination, overwrite = false)
    vfa_folder = joinpath(destination, study_id, "vfa")
    dce_folder = joinpath(destination, study_id, "dce")

    if isdir(vfa_folder) && isdir(dce_folder)
        if !overwrite
            return (vfa_folder, dce_folder)
        end
    end

    gbm_series = tcia_series(study = study_id)
    vfa_series = find_vfa_series(gbm_series)
    dce_series = find_dce_series(gbm_series)

    if !isdir(vfa_folder) || overwrite == true
        download_series(vfa_series, destination = vfa_folder)
    end

    if !isdir(dce_folder) || overwrite == true
        download_series(dce_series, destination = dce_folder)
    end
    return (vfa_folder, dce_folder)
end

find_vfa_series(series_dataframe) = find_in_description("MAP", series_dataframe)
find_dce_series(series_dataframe) = find_in_description("DYN", series_dataframe)

function find_in_description(word_to_find::AbstractString, series_dataframe)
    descriptions = series_dataframe.SeriesDescription
    found_indices = findall(occursin.(word_to_find, descriptions))
    if length(found_indices) < 1
        errant_id = series_dataframe.PatientID[1]
        errant_study = series_dataframe.StudyInstanceUID[1]
        error("No single series in $errant_id found with $word_to_find in their description.
              The full study UID is $errant_study\n")
    elseif length(found_indices) > 1
        errant_id = series_dataframe.PatientID[1]
        @warn "Multiple series in $errant_id found containing $word_to_find in their description.
              An arbitrary series will be used among the found series"
        found_index = found_indices[end]
    else
        found_index = found_indices[1]
    end
    found_series = series_dataframe.SeriesInstanceUID[found_index]
    return found_series
end

function download_series(series_id::AbstractString; destination)
    make_folder(destination; remove_existing = true)
    zip_file = joinpath(destination, "downloaded.zip")
    tcia_images(series = series_id, file = zip_file)
    unzip_command = `unzip -o $zip_file -d $destination`
    run(unzip_command)
    rm(zip_file)
    return destination
end

function make_folder(desired_folder::AbstractString; remove_existing = false)
    if !isdir(desired_folder)
        mkpath(desired_folder)
    elseif remove_existing
        rm(desired_folder, recursive = true)
        mkpath(desired_folder)
    end
    return desired_folder
end


text:
@testset "ExternalSamplingMCCFR" begin
    env = OpenSpielEnv(
        "kuhn_poker";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = ExternalSamplingMCCFRPolicy(rng = StableRNG(123))
    run(p, env, StopAfterStep(1000))
    @test_broken RLZoo.nash_conv(p, env) < 0.05

    env = OpenSpielEnv(
        "leduc_poker";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = ExternalSamplingMCCFRPolicy(rng = StableRNG(123))
    run(p, env, StopAfterStep(1000))
    @test RLZoo.nash_conv(p, env) < 2.5

    env = OpenSpielEnv(
        "liars_dice";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = ExternalSamplingMCCFRPolicy(rng = StableRNG(123))
    run(p, env, StopAfterStep(100))
    @test RLZoo.nash_conv(p, env) < 1.6

    env = OpenSpielEnv(
        "kuhn_poker(players=3)";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = ExternalSamplingMCCFRPolicy(rng = StableRNG(123))
    run(p, env, StopAfterStep(100))
    @info "nash_conv for kuhn_poker(players=3)" RLZoo.nash_conv(p, env)
end

@testset "OutcomeSamplingMCCFR" begin
    env = OpenSpielEnv(
        "kuhn_poker";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = OutcomeSamplingMCCFRPolicy(rng = StableRNG(123))
    run(p, env, StopAfterStep(10000))
    @test_broken RLZoo.nash_conv(p, env) < 0.04

    env = OpenSpielEnv(
        "leduc_poker";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = OutcomeSamplingMCCFRPolicy(rng = StableRNG(123))
    run(p, env, StopAfterStep(10000))
    @test RLZoo.nash_conv(p, env) < 3

    env = OpenSpielEnv(
        "liars_dice";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = OutcomeSamplingMCCFRPolicy(rng = StableRNG(123))
    run(p, env, StopAfterStep(1000))
    @test RLZoo.nash_conv(p, env) < 1.7

    env = OpenSpielEnv(
        "kuhn_poker(players=3)";
        default_state_style = RLBase.Information{String}(),
        is_chance_agent_required = true,
    )
    p = OutcomeSamplingMCCFRPolicy(rng = StableRNG(123))
    run(p, env, StopAfterStep(1000))
    @info "nash_conv for kuhn_poker(players=3)" RLZoo.nash_conv(p, env)
end


text:
using CoordinateConverterGK
using Test

@testset "CoordinateConverterGK.jl" begin
    # Write your own tests here.

    # xy to lonlat
    northing = 11573.375
    easting = 22694.980
    φ, λ = yx2latlon_ja(9, northing, easting)
    @test (abs(φ - 36.10404755) < 1e-4) & (abs(λ - 140.08539843) < 1e-4)
    λ, φ = xy2lonlat(reverse(CoordinateConverterGK.Origin_LatLon_Japan[9,:])..., easting, northing)
    @test (abs(φ - 36.10404755) < 1e-4) & (abs(λ - 140.08539843) < 1e-4)

    # lonlat to xy
    lat = 36.103774791666666
    lon = 140.08785504166664
    y, x = latlon2yx_ja(9, lat, lon)
    @test (abs(y - 11543.6883) < 1e-4) & (abs(x - 22916.2436) < 1e-4)
    λ, φ = lonlat2xy(reverse(CoordinateConverterGK.Origin_LatLon_Japan[9,:])..., lon, lat)
    @test (abs(y - 11543.6883) < 1e-4) & (abs(x - 22916.2436) < 1e-4)

    # args in reverse order
    @test yx2latlon_ja(9, northing, easting) == reverse(xy2lonlat_ja(9, easting, northing), dims=2)
    @test latlon2yx_ja(9, lat, lon) == reverse(lonlat2xy_ja(9, lon, lat), dims=2)
end


text:
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIB = MOI.Bridges
const MOIU = MOI.Utilities

sympackedlen(n) = div(n*(n+1), 2)
sympackeddim(n) = div(isqrt(1+8n) - 1, 2)
function ivech!(out::AbstractMatrix{T}, v::AbstractVector{T}) where T
    n = sympackeddim(length(v))
    n1, n2 = size(out)
    @assert n == n1 == n2
    c = 0
    for j in 1:n, i in 1:j
        c += 1
        out[i,j] = v[c]
    end
    return out
end
function ivech(v::AbstractVector{T}) where T
    n = sympackeddim(length(v))
    out = zeros(n, n)
    ivech!(out, v)

    return out
end

solvers = Tuple{String, Any}[]

#MOIU.@model ProxSDPModelData () (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan) (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.PositiveSemidefiniteConeTriangle) () (MOI.SingleVariable,) (MOI.ScalarAffineFunction,) (MOI.VectorOfVariables,) (MOI.VectorAffineFunction,)

#=
    CSDP & SDPA
=#

# ]activate soi
# using CSDP
# using SDPA

# MOIU.@model(SOI_ModelData,
#             (),
#             (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
#             (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
#              MOI.PositiveSemidefiniteConeTriangle),
#             (),
#             (MOI.SingleVariable,),
#             (MOI.ScalarAffineFunction,),
#             (MOI.VectorOfVariables,),
#             (MOI.VectorAffineFunction,))

# const CSDP_optimizer =
# MOIB.full_bridge_optimizer(
#     MOIU.CachingOptimizer(
#         SOI_ModelData{Float64}(), CSDP.Optimizer()),#printlevel=0))
#     Float64
#     )

# push!(solvers, ("CSDP", CSDP_optimizer)) # eps = ???

# const SDPA_optimizer =
# MOIB.full_bridge_optimizer(
#     MOIU.CachingOptimizer(
#         SOI_ModelData{Float64}(), SDPA.Optimizer()),#printlevel=0))
#     Float64
#     )

# push!(solvers, ("SDPA", SDPA_optimizer)) # eps = ???

#=
    COSMO
=#

# using COSMO

# MOIU.@model(COSMOModelData,
#         (),
#         (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
#         (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
#          MOI.PositiveSemidefiniteConeSquare, MOI.PositiveSemidefiniteConeTriangle, MOI.ExponentialCone, MOI.DualExponentialCone),
#         (MOI.PowerCone, MOI.DualPowerCone),
#         (MOI.SingleVariable,),
#         (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
#         (MOI.VectorOfVariables,),
#         (MOI.VectorAffineFunction,),);

# const COSMO_optimizer = MOIU.CachingOptimizer(MOIU.UniversalFallback(COSMOModelData{Float64}()),
#     COSMO.Optimizer(eps_abs = 1e-4, eps_rel = 1e-4 ));

# push!(solvers, ("COSMO", COSMO_optimizer))

#=
    SeDuMi
=#

import SeDuMi
const optimizer = SeDuMi.Optimizer(fid=0)

MOIU.@model(SeDuMi_ModelData, (), (),
            (MOI.Zeros, MOI.Nonnegatives, MOI.SecondOrderCone,
             MOI.RotatedSecondOrderCone, MOI.PositiveSemidefiniteConeTriangle),
            (), (), (), (MOI.VectorOfVariables,), (MOI.VectorAffineFunction,))

const SeDuMi_optimizer = MOIB.full_bridge_optimizer(MOIU.CachingOptimizer(
    MOIU.UniversalFallback(SeDuMi_ModelData{Float64}()),
    SeDuMi.Optimizer(fid=0)),
    Float64)

push!(solvers, ("SeDuMi", SeDuMi_optimizer))


text:
"""
closer(ps::ParseState)

A magical function determining whether the parsing of an expression should continue or
stop.
"""
function closer(ps::ParseState)
    kindof(ps.nt) === Tokens.ENDMARKER ||
    (ps.closer.newline && kindof(ps.ws) == NewLineWS && !iscomma(ps.t)) ||
    (ps.closer.semicolon && kindof(ps.ws) == SemiColonWS) ||
    (isoperator(ps.nt) && precedence(ps.nt) <= ps.closer.precedence) ||
    (kindof(ps.nt) === Tokens.WHERE && ps.closer.precedence == LazyAndOp) ||
    (ps.closer.inwhere && kindof(ps.nt) === Tokens.WHERE) ||
    (ps.closer.inwhere && ps.closer.ws && kindof(ps.t) === Tokens.RPAREN && isoperator(ps.nt) && precedence(ps.nt) < DeclarationOp) ||
    (ps.closer.precedence > WhereOp && (
        (kindof(ps.nt) === Tokens.LPAREN && !(kindof(ps.t) === Tokens.EX_OR)) ||
        kindof(ps.nt) === Tokens.LBRACE ||
        kindof(ps.nt) === Tokens.LSQUARE ||
        (kindof(ps.nt) === Tokens.STRING && isemptyws(ps.ws)) ||
        ((kindof(ps.nt) === Tokens.RPAREN || kindof(ps.nt) === Tokens.RSQUARE) && isidentifier(ps.nt))
    )) ||
    (iscomma(ps.nt) && ps.closer.precedence > AssignmentOp) ||
    kindof(ps.nt) === Tokens.ENDMARKER ||
    (ps.closer.comma && iscomma(ps.nt)) ||
    (ps.closer.tuple && (iscomma(ps.nt) || isassignmentop(ps.nt))) ||
    (kindof(ps.nt) === Tokens.FOR && ps.closer.precedence > -1) ||
    (ps.closer.block && kindof(ps.nt) === Tokens.END) ||
    (ps.closer.paren && kindof(ps.nt) === Tokens.RPAREN) ||
    (ps.closer.brace && kindof(ps.nt) === Tokens.RBRACE) ||
    (ps.closer.square && kindof(ps.nt) === Tokens.RSQUARE) ||
    # tilde parsing in vect exprs needs to be special cased because `~` has assignment precedence
    (@static VERSION < v"1.4" ?
        false :
        ((ps.closer.insquare || ps.closer.inmacro) && kindof(ps.nt) === Tokens.APPROX && !isemptyws(ps.ws) && isemptyws(ps.nws))
    ) ||
    kindof(ps.nt) === Tokens.ELSEIF ||
    kindof(ps.nt) === Tokens.ELSE ||
    kindof(ps.nt) === Tokens.CATCH ||
    kindof(ps.nt) === Tokens.FINALLY ||
    (ps.closer.ifop && isoperator(ps.nt) && (precedence(ps.nt) <= 0 || kindof(ps.nt) === Tokens.COLON)) ||
    (ps.closer.range && (kindof(ps.nt) === Tokens.FOR || iscomma(ps.nt) || kindof(ps.nt) === Tokens.IF)) ||
    (ps.closer.ws && !isemptyws(ps.ws) &&
        !iscomma(ps.nt) &&
        !iscomma(ps.t) &&
        !(!ps.closer.inmacro && kindof(ps.nt) === Tokens.FOR) &&
        !(kindof(ps.nt) === Tokens.DO) &&
        !(
            (isbinaryop(ps.nt) && !(ps.closer.wsop && isemptyws(ps.nws) && isunaryop(ps.nt) && precedence(ps.nt) > 7)) ||
            (isunaryop(ps.t) && kindof(ps.ws) == WS && kindof(ps.lt) !== CSTParser.Tokens.COLON)
        )) ||
    (ps.closer.unary && (kindof(ps.t) in (Tokens.INTEGER, Tokens.FLOAT, Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE) && isidentifier(ps.nt)))
end

"""
    @closer ps rule body

Continues parsing closing on `rule`.
"""
macro closer(ps, opt, body)
    quote
        local tmp1 = getfield($(esc(ps)).closer, $opt)
        setfield!($(esc(ps)).closer, $opt, true)
        out = $(esc(body))
        setfield!($(esc(ps)).closer, $opt, tmp1)
        out
    end
end

"""
    @nocloser ps rule body

Continues parsing not closing on `rule`.
"""
macro nocloser(ps, opt, body)
    quote
        local tmp1 = getfield($(esc(ps)).closer, $opt)
        setfield!($(esc(ps)).closer, $opt, false)
        out = $(esc(body))
        setfield!($(esc(ps)).closer, $opt, tmp1)
        out
    end
end

macro closeparen(ps, body)
    quote
        local tmp1 = $(esc(ps)).closer.paren
        $(esc(ps)).closer.paren = true
        out = $(esc(body))
        $(esc(ps)).closer.paren = tmp1
        out
    end
end

macro closesquare(ps, body)
    quote
        local tmp1 = $(esc(ps)).closer.square
        $(esc(ps)).closer.square = true
        out = $(esc(body))
        $(esc(ps)).closer.square = tmp1
        out
    end
end
macro closebrace(ps, body)
    quote
        local tmp1 = $(esc(ps)).closer.brace
        $(esc(ps)).closer.brace = true
        out = $(esc(body))
        $(esc(ps)).closer.brace = tmp1
        out
    end
end

"""
    @precedence ps prec body

Continues parsing binary operators until it hits a more loosely binding
operator (with precdence lower than `prec`).
"""
macro precedence(ps, prec, body)
    quote
        local tmp1 = $(esc(ps)).closer.precedence
        $(esc(ps)).closer.precedence = $(esc(prec))
        out = $(esc(body))
        $(esc(ps)).closer.precedence = tmp1
        out
    end
end

# Closer_TMP and ancillary functions help reduce code generation
struct Closer_TMP
    newline::Bool
    semicolon::Bool
    inmacro::Bool
    tuple::Bool
    comma::Bool
    insquare::Bool
    range::Bool
    ifop::Bool
    ws::Bool
    wsop::Bool
    unary::Bool
    precedence::Int
end

@noinline function create_tmp(c::Closer)
    Closer_TMP(c.newline,
        c.semicolon,
        c.inmacro,
        c.tuple,
        c.comma,
        c.insquare,
        c.range,
        c.ifop,
        c.ws,
        c.wsop,
        c.unary,
        c.precedence)
end

@noinline function update_from_tmp!(c::Closer, tmp::Closer_TMP)
    c.newline = tmp.newline
    c.semicolon = tmp.semicolon
    c.inmacro = tmp.inmacro
    c.tuple = tmp.tuple
    c.comma = tmp.comma
    c.insquare = tmp.insquare
    c.range = tmp.range
    c.ifop = tmp.ifop
    c.ws = tmp.ws
    c.wsop = tmp.wsop
    c.unary = tmp.unary
    c.precedence = tmp.precedence
end

@noinline function update_to_default!(c::Closer)
    c.newline = true
    c.semicolon = true
    c.inmacro = false
    c.tuple = false
    c.comma = false
    c.insquare = false
    c.range = false
    c.ifop = false
    c.ws = false
    c.wsop = false
    c.unary = false
    c.precedence = -1
end

"""
    @default ps body

Parses the next expression using default closure rules.
"""
macro default(ps, body)
    quote
        TMP = create_tmp($(esc(ps)).closer)
        update_to_default!($(esc(ps)).closer)
        out = $(esc(body))
        update_from_tmp!($(esc(ps)).closer, TMP)
        out
    end
end

isajuxtaposition(ps::ParseState, ret::EXPR) = ((isnumber(ret) && (isidentifier(ps.nt) || kindof(ps.nt) === Tokens.LPAREN || kindof(ps.nt) === Tokens.CMD || kindof(ps.nt) === Tokens.STRING || kindof(ps.nt) === Tokens.TRIPLE_STRING)) ||
        ((is_prime(ret.head) && isidentifier(ps.nt)) ||
        ((kindof(ps.t) === Tokens.RPAREN || kindof(ps.t) === Tokens.RSQUARE) && (isidentifier(ps.nt) || kindof(ps.nt) === Tokens.CMD)) ||
        ((kindof(ps.t) === Tokens.STRING || kindof(ps.t) === Tokens.TRIPLE_STRING) && (kindof(ps.nt) === Tokens.STRING || kindof(ps.nt) === Tokens.TRIPLE_STRING)))) || ((kindof(ps.t) in (Tokens.INTEGER, Tokens.FLOAT) || kindof(ps.t) in (Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE)) && isidentifier(ps.nt)) ||
        (isnumber(ret) && ps.closer.inref && (ps.nt.kind === Tokens.END || ps.nt.kind === Tokens.BEGIN))

"""
    has_error(ps::ParseState)
    has_error(x::EXPR)

Determine whether a parsing error occured while processing text with the given
`ParseState`, or exists as a (sub) expression of `x`.
"""
function has_error(x::EXPR)
    return headof(x) == :errortoken || (x.args !== nothing && any(has_error, x.args)) || (x.trivia !== nothing && any(has_error, x.trivia))
end
has_error(ps::ParseState) = ps.errored

using Base.Meta

"""
    compare(x,y)

Recursively checks whether two Base.Expr are the same. Returns unequal sub-
expressions.
"""
compare(x, y) = x == y ? true : (x, y)

function compare(x::Expr, y::Expr)
    if x == y
        return true
    else
        if x.head != y.head
            return (x, y)
        end
        if length(x.args) != length(y.args)
            return (x.args, y.args)
        end
        for i = 1:length(x.args)
            t = compare(x.args[i], y.args[i])
            if t != true
                return t
            end
        end
    end
end

# code for updating CST
"""
    firstdiff(s0::AbstractString, s1::AbstractString)

Returns the last byte index, i, for which s0 and s1 are the same such that:
    `s0[1:i] == s1[1:i]`
"""
function firstdiff(s0::AbstractString, s1::AbstractString)
    minlength = min(sizeof(s0), sizeof(s1))
    @inbounds for i in 1:minlength
        if codeunits(s0)[i] !== codeunits(s1)[i]
            return i - 1 # This could return a non-commencing byte of a multi-byte unicode sequence.
        end
    end
    return minlength
end

"""
    revfirstdiff(s0::AbstractString, s1::AbstractString)

Reversed version of firstdiff but returns two indices, one for each string.
"""
function revfirstdiff(s0::AbstractString, s1::AbstractString)
    minlength = min(sizeof(s0), sizeof(s1))
    @inbounds for i in 0:minlength - 1
        if codeunits(s0)[end - i] !== codeunits(s1)[end - i]
            return sizeof(s0) - i, sizeof(s1) - i# This could return a non-commencing byte of a multi-byte unicode sequence.
        end
    end
    return 1, 1
end

"""
    find_arg_at(x, i)

Returns the index of the node of `x` within which the byte offset `i` falls.
"""
function find_arg_at(x::CSTParser.EXPR, i)
    @assert i <= x.fullspan
    offset = 0
    for (cnt, a) in enumerate(x.args)
        if i <= offset + a.fullspan
            return cnt
        end
        offset += a.fullspan
    end
    error("$(x.head) with fullspan: $(x.fullspan) at $i")
end

comp(x, y) = x == y
function comp(x::CSTParser.EXPR, y::CSTParser.EXPR)
    comp(x.head, y.head) &&
    x.span == y.span &&
    x.fullspan == y.fullspan &&
    x.val == y.val &&
    length(x) == length(y) &&
    all(comp(x[i], y[i]) for i = 1:length(x))
end

function minimal_reparse(s0, s1, x0 = CSTParser.parse(s0, true), x1 = CSTParser.parse(s1, true); inds = false)
    if has_error(x0)
        return inds ? (1:0, 1:length(x1.args), 1:0) : x1 # Error while re-parsing, so lets return the whole expression instead of patching
    end

    if sizeof(s0) !== x0.fullspan
       error("minimal reparse - original input text length doesn't match the full span of the provided CST.")
        # return inds ? (1:0, 1:length(x1.args), 1:0) : x1
    end

    if has_error(x1)
        return inds ? (1:0, 1:length(x1.args), 1:0) : x1 # Error while re-parsing, so lets return the whole expression instead of patching
    end

    if sizeof(s1) !== x1.fullspan
        error("minimal reparse - new input text length doesn't match the full span of the provided CST.")
         # return inds ? (1:0, 1:length(x1.args), 1:0) : x1
    end
    isempty(x0.args) && return inds ? (1:0, 1:length(x1.args), 1:0) : x1 # Original CST was empty
    x1.fullspan == 0 && return inds ? (1:0, 1:0, 1:0) : x1 # New CST is empty

    i0 = firstdiff(s0, s1)
    i0 > x0.fullspan && return inds ? (1:0, 1:length(x1.args), 1:0) : x1 # Should error?
    i1, i2 = revfirstdiff(s0, s1)
    (i0 > x1.fullspan || i1 > x1.fullspan || i2 > x1.fullspan) && return inds ? (1:0, 1:length(x1.args), 1:0) : x1 # Should error?
    # Find unaffected expressions at start
    # CST should be unaffected (and able to be copied across) up to this point,
    # but we need to check.
    r1 = 1:min(find_arg_at(x0, i0) - 1, length(x0.args), find_arg_at(x1, i0) - 1)
    for i = 1:min(find_arg_at(x0, i0) - 1, find_arg_at(x1, i0) - 1)
        if x0.args[i].fullspan !== x1.args[i].fullspan
            r1 = 1:(i-1)
            break
        end
        r1 = 1:i
    end
    # we can re-use x0.args[r1]

    # assume we'll just use x1.args from here on
    r2 = (last(r1) + 1):length(x1.args)
    r3 = 0:-1

    # though we now check whether there is a sequence at the end of x0.args and
    # x1.args that match
    offset = sizeof(s1)
    for i = 0:min(length(x0.args) - last(r1), length(x0.args), length(x1.args)) - 1
        if !quick_comp(x0.args[end - i], x1.args[end - i]) ||
            offset <= i1 ||
            length(x0.args) - i == last(r1) + 1 ||
            offset - x1.args[end-i].fullspan <= i2 <= offset

            r2 = first(r2):length(x1.args) - i
            r3 = length(x0.args) .+ ((-i + 1):0)
            break
        end
        offset -= x1.args[end - i].fullspan
    end
    inds && return r1, r2, r3
    x2 = CSTParser.EXPR(x0.head, CSTParser.EXPR[
        x0.args[r1]
        x1.args[r2]
        x0.args[r3]
    ], nothing)
    return x2
end

# Quick and very dirty comparison of two EXPR, makes extra effort for :errortokens
function quick_comp(a::EXPR, b::EXPR)
    a.fullspan != b.fullspan && return false
    if headof(a) === :errortoken
        headof(b) !== :errortoken && return false
        if a.args !== nothing
            b.args === nothing && return false
            return length(a.args) == length(b.args) && (length(a.args) == 0 || quick_comp(first(a.args), first(b.args)))
        end
    end
    return comp(headof(a), headof(b))
end

"""
check_span(x, neq = [])

Recursively checks whether the span of an expression equals the sum of the span
of its components. Returns a vector of failing expressions.
"""
function check_span(x::EXPR, neq = [])
    (ispunctuation(x) || isidentifier(x) || iskeyword(x) || isoperator(x) || isliteral(x) || headof(x) == :string) && return neq

    s = 0
    if x.args !== nothing
        for a in x.args
            check_span(a, neq)
            s += a.fullspan
        end
    end
    if hastrivia(x)
        for a in x.trivia
            check_span(a, neq)
            s += a.fullspan
        end
    end
    if x.head isa EXPR
        s += x.head.fullspan
    end
    if length(x) > 0 && s != x.fullspan
        push!(neq, x)
    end
    neq
end

function speed_test()
    dir = dirname(Base.find_source_file("essentials.jl"))
    println("speed test : ", @timed(for i = 1:5
        parse(read(joinpath(dir, "essentials.jl"), String), true);
        parse(read(joinpath(dir, "abstractarray.jl"), String), true);
    end)[2])
end

"""
    str_value(x)

Attempt to get a string representation of a nodeless expression.
"""
function str_value(x)
    if headof(x) === :IDENTIFIER || isliteral(x)
        return valof(x)
    elseif isidentifier(x)
        valof(x.args[2])
    elseif isoperator(x)
        return string(Expr(x))
    else
        return ""
    end
end

_unescape_string(s::AbstractString) = sprint(_unescape_string, s, sizehint=lastindex(s))
function _unescape_string(io, s::AbstractString)
    a = Iterators.Stateful(s)
    for c in a
        if !isempty(a) && c == '\\'
            c = popfirst!(a)
            if c == 'x' || c == 'u' || c == 'U'
                n = k = 0
                m = c == 'x' ? 2 :
                    c == 'u' ? 4 : 8
                while (k += 1) <= m && !isempty(a)
                    nc = Base.peek(a)
                    n = '0' <= nc <= '9' ? n << 4 + nc - '0' :
                        'a' <= nc <= 'f' ? n << 4 + nc - 'a' + 10 :
                        'A' <= nc <= 'F' ? n << 4 + nc - 'A' + 10 : break
                    popfirst!(a)
                end
                if k == 1
                    # throw(ArgumentError("invalid $(m == 2 ? "hex (\\x)" :
                    #                         "unicode (\\u)") escape sequence used in $(repr(s))"))
                    # push error to ParseState?
                    n = 0
                end
                if m == 2 # \x escape sequence
                    write(io, UInt8(n))
                else
                    print(io, Char(n))
                end
            elseif '0' <= c <= '7'
                k = 1
                n = c - '0'
                while (k += 1) <= 3 && !isempty(a)
                    c  = Base.peek(a)
                    n = ('0' <= c <= '7') ? n << 3 + c - '0' : break
                    popfirst!(a)
                end
                if n > 255
                    # throw(ArgumentError("octal escape sequence out of range"))
                    # push error to ParseState?
                    n = 255
                end
                write(io, UInt8(n))
            else
                print(io, c == 'a' ? '\a' :
                          c == 'b' ? '\b' :
                          c == 't' ? '\t' :
                          c == 'n' ? '\n' :
                          c == 'v' ? '\v' :
                          c == 'f' ? '\f' :
                          c == 'r' ? '\r' :
                          c == 'e' ? '\e' : c)
            end
        else
            print(io, c)
        end
    end
end

function valid_escaped_seq(s::AbstractString)
    l = length(s)
    l == 0 && return false # zero length chars are always invalid
    l == 1 && return true # length-one chars are always valid to Julia's parser
    a = Iterators.Stateful(s)
    if popfirst!(a) == '\\'
        c = popfirst!(a)
        if c === 'x' || c === 'u' || c === 'U'
            maxiter = c === 'x' ? 2 : c === 'u' ? 4 : 8
            0 < length(a) <= maxiter || return false
            n = 0
            while !isempty(a)
                nc = popfirst!(a)
                n = '0' <= nc <= '9' ? n << 4 + (nc - '0') :
                    'a' <= nc <= 'f' ? n << 4 + (nc - 'a' + 10) :
                    'A' <= nc <= 'F' ? n << 4 + (nc - 'A' + 10) : return false
            end
            return n <= 0x10ffff
        elseif '0' <= c <= '7'
            length(a) <= 3 || return false
            n = c - '0'
            while !isempty(a)
                nc = popfirst!(a)
                n = ('0' <= c <= '7') ? n << 3 + nc - '0' : return false
            end
            return n < 128
        else
            @static if VERSION < v"1.1.0"
                c = string(c)
            end
            return ncodeunits(c) == 1 && isempty(a)
        end
    end
    return false
end

"""
    disallowednumberjuxt(ret::EXPR)

Does this number literal end in a decimal and so cannot precede a paren for
implicit multiplication?
"""
disallowednumberjuxt(ret::EXPR) = isnumber(ret) && last(valof(ret)) == '.'

nexttokenstartsdocstring(ps::ParseState) = isidentifier(ps.nt) && val(ps.nt, ps) == "doc" && (kindof(ps.nnt) === Tokens.STRING || kindof(ps.nnt) === Tokens.TRIPLE_STRING)

"""
    _do_kw_convert(ps::ParseState, a::EXPR)

Should `a` be converted to a keyword-argument expression?
"""
_do_kw_convert(ps::ParseState, a::EXPR) = !ps.closer.brace && isassignment(a)

"""
    _kw_convert(ps::ParseState, a::EXPR)

Converted an assignment expression to a keyword-argument expression.
"""
_kw_convert(x::EXPR) = EXPR(:kw, EXPR[x.args[1], x.args[2]], EXPR[x.head], x.fullspan, x.span)

"""
    convertsigtotuple(sig::EXPR)

When parsing a function or macro signature, should it be converted to a tuple?
"""
convertsigtotuple(sig::EXPR) = isbracketed(sig) && !(istuple(sig.args[1]) || (headof(sig.args[1]) === :block) || issplat(sig.args[1]))

"""
    docable(head)

When parsing a block of expressions, can documentation be attached? Prefixed docs at the
top-level are handled within `parse(ps::ParseState, cont = false)`.
"""
docable(head) = head === :begin || head === :module || head === :baremodule || head === :quote

should_negate_number_literal(ps::ParseState, op::EXPR) = (is_plus(op) || is_minus(op)) && (kindof(ps.nt) === Tokens.INTEGER || kindof(ps.nt) === Tokens.FLOAT) && isemptyws(ps.ws) && kindof(ps.nnt) != Tokens.CIRCUMFLEX_ACCENT

"""
    can_become_comparison(x::EXPR)

Is `x` a binary comparison call (e.g. `a < b`) that can be extended to include more
arguments?
"""
can_become_comparison(x::EXPR) = (isoperator(x.head) && comp_prec(valof(x.head)) && length(x.args) > 1) || (x.head === :call && isoperator(x.args[1]) && comp_prec(valof(x.args[1])) && length(x.args) > 2 && !hastrivia(x))

"""
    can_become_chain(x::EXPR, op::EXPR)

Is `x` a binary call for `+` or `*` that can be extended to include more
arguments?
"""
can_become_chain(x::EXPR, op::EXPR) = isbinarycall(x) && (is_star(op) || is_plus(op)) && valof(op) == valof(x.args[2]) && !isdotted(x.args[2]) && x.args[2].span > 0

macro cst_str(x)
    CSTParser.parse(x)
end

function issuffixableliteral(ps::ParseState, x::EXPR)
    # prefixed string/cmd macros can be suffixed by identifiers or numeric literals
    (isidentifier(ps.nt) || isnumberliteral(ps.nt) || isbool(ps.nt)) && isemptyws(ps.ws) && ismacrocall(x) && (valof(x.args[1]) isa String && (endswith(valof(x.args[1]), "_str") || endswith(valof(x.args[1]), "_cmd")))
end

function loop_check(ps, prevpos)
    if position(ps) <= prevpos && ps.nt.kind !== Tokens.ENDMARKER
        throw(CSTInfiniteLoop("Infinite loop at $ps"))
    else
        position(ps)
    end
end


text:
"""
        SegyHeaderInfo(filename_in)

Print segy header information to screen. The input is the name of the headers file

"""
function SegyHeaderInfo(h_segy)

#requires using Printf, Statistics

println("Calculating statistics.")

#key gets all the types of the composite type SegyHeader).
key = fieldnames(SeisMain.SegyHeader)

min_h = zeros(Float64,length(key))
max_h = zeros(Float64,length(key))
mean_h = zeros(Float64,length(key))

for ikey=1:length(key)
		min_h[ikey] = convert(Float64,getfield(h_segy[1],key[ikey]))
		max_h[ikey] = convert(Float64,getfield(h_segy[1],key[ikey]))
		mean_h[ikey] += convert(Float64,getfield(h_segy[1],key[ikey]))
end

nx = size(h_segy,1)

for itrace = 2:nx
	for ikey = 1 : length(key)
		aux = convert(Float64,getfield(h_segy[itrace],ikey))
		if (aux < min_h[ikey])
			min_h[ikey] = aux
		end
		if ( aux > max_h[ikey])
			max_h[ikey] = aux
		end
		mean_h[ikey] += aux
	end
end

for ikey=1:length(key)
	mean_h[ikey] /= nx
end

	println("       Key          Minimum          Maximum             Mean");
	println("=============================================================")
	for ikey=1:length(key)
		@printf("%10s      %11.3f      %11.3f      %11.3f\n",string(key[ikey]),min_h[ikey],max_h[ikey],mean_h[ikey])
	end

end


text:
using DataFrames, PyCall
using KymatioWrapper
using Test
py"""
import torch
from kymatio.scattering1d.filter_bank import scattering_filter_factory
from kymatio import Scattering1D
from kymatio import Scattering2D
from kymatio.scattering2d.filter_bank import filter_bank

torch.cuda.is_available()
"""

@testset "scattering.jl" begin
    useGpu = py"torch.cuda.is_available()"
    J = 6; Q = 16; N = 2^13 - 2045
    T = 2^ceil.(Int, log2(N))
    Scatter = py"Scattering1D($J,$N,Q=$Q, max_order=$2, average=$true, oversampling=$0, backend='torch')"
    if useGpu
        ScatterGpu = py"Scattering1D($J,$N,Q=$Q, max_order=$2, average=$true, oversampling=$0, backend='torch')"
        ScatterGpu.cuda()
    end
    t = (range(0, 2π, length=N));
    x = zeros(Float32, N, 5, 3);
    for (ii, ω) in enumerate(floor.(Int, range(1, 150, length=5))), jj in 1:3
        x[:,ii,jj] = sin.(ω / jj .* t) + cos.(100 * ω / jj .* t.^2)
    end
    scX = Scatter.forward(py"torch.from_numpy($(permutedims(x,(2,3,1)))).contiguous()")
    res = scX.numpy()
    s = Scattering(N, Q=16, useGpu=false)
    sx = s(x)
    @test sx ≈ permutedims(res, (4, 3, 1, 2))
    if useGpu
        scXGpu = ScatterGpu.forward(py"torch.from_numpy($(permutedims(x,(2,3,1)))).cuda().contiguous()")
        resGpu = scX.cpu().numpy()
        s = Scattering(N, Q=16, useGpu=true)
        sx = s(x)
        @test sx ≈ permutedims(resGpu, (4, 3, 1, 2))
    end
    N = (26, 26); J = 2
    x = randn(N..., 4, 1); size(x)
    w = Float32.(cat(x, zeros(6, 26, 4, 1), dims=1)); size(w)
    wx = Float32.(cat(w, zeros(32, 6, 4, 1), dims=2)); size(wx)
    s = Scattering((26, 26), useGpu=false)
    Scatter = py"Scattering2D(3, shape=(32,32), L=8,backend='torch')"
    res = Scatter.forward(py"torch.from_numpy($(permutedims(wx,(3,4,1,2)))).contiguous()")
    res = res.numpy()
    @test s(x) ≈ permutedims(res, (5, 4, 3, 1:2...))
    if useGpu
        ScatterGpu = py"Scattering2D(3, shape=(32,32), L=8,backend='torch').cuda()"
        resGpu = Scatter.forward(py"torch.from_numpy($(permutedims(wx,(3,4,1,2)))).contiguous()")
        resGpu = resGpu.cpu().numpy()
        sGpu = Scattering((26, 26), useGpu=useGpu)
        @test sGpu(x) ≈ permutedims(resGpu, (5, 4, 3, 1:2...))
    end
end

@testset "filter_bank" begin
    N = 450
    nψ₁ = 23
    nψ₂ = 6
    ψ₁, ψ₂, ϕ₃, params = filter_bank(N)
    @test size(ψ₁) == (2^ceil(log2(N)), nψ₁)
    @test size(ψ₂) == (2^ceil(log2(N)), nψ₂)
    @test size(ϕ₃) == (2^ceil(log2(N)),)
    @test size(params[1]) == (nψ₁, 3)
    @test size(params[2]) == (nψ₂, 3)
    N = (25, 25); J = 2; L = 8
    JL = J * L + 1
    ψ, params, ϕₙ = filter_bank(N)
    @test size(ψ) == (N..., JL)
    @test size(params) == (JL, 2)
    @test length(ϕₙ) == J
    @test size(ϕₙ[1]) == N
    @test size(ϕₙ[2]) == N .>> 1
end


text:
using CxxAutoWrap
using Base.Test
using Cxx
using Clang
using Clang.cindex

dir = Pkg.dir("CxxAutoWrap")
folder_intcontainer = joinpath(dir, "deps", "examples", "IntContainer")
addHeaderDir(folder_intcontainer)
path = joinpath(folder_intcontainer,"library.h")

cxxinclude(path)

libpath = joinpath(folder_intcontainer, "cmake-build-debug", "libIntContainer.so")
Libdl.dlopen(libpath, Libdl.RTLD_GLOBAL)

(wrapexpr_header(path)) |> eval

@testset "IntContainer" begin

    c = IntContainer(1)
    @test getIt(c) == 1
    setIt(c, 10)
    @test getIt(c) == 10
    delete(c)

end
@testset "StringContainer" begin
    s = StringContainer("asd")
    @test getIt(s) |> unsafe_string == "asd"
    setIt(s, "42")
    @test getIt(s) |> unsafe_string == "42"

end


text:
"""percentile filter based on updating histograms while sliding over the image
implementation of Perreault and Hebert (2007) algorithm see: http://nomis80.org/ctmf.pdf

precision is limited by number of bins.
for a images with a bit depth of 8 bit (e.g. greyscale images), full precision
can be acchieved with a dual-level histogram of 16 bins
For images with a higher bit depth more bins are needed and the algorithm slows,
the precision can be set with the optional precision paramter

input
A           2d image Array
p           percentile value [1-100]
r           radius of window [n. cells]: i.e.: the window size is 2r+1 by 2r+1
is8bit      true if A has a bit-depth of 8 bits and no scaling needs to be applied (default false)
precision   number of decimals to preserve when binning
"""
function hist_filter(A::Array{T,2}, perc::Real, r::Int;
                                precision=0,
                                nodata_value=T(-9999.0),
                                hist_min=nothing,
                                hist_max=nothing) where T <: Real
    scale = T(10^precision) # default scale value
    ksize = (2*r+1)^2 # kernel size
    nrows, ncols = size(A) # image size
    mask = (A .≠ nodata_value) .& isfinite.(A)
    if T == UInt8
        im_min, im_max = 0, 255
    elseif T <: Integer
        precision > 0 && info("no scaling applied as A has integer values")
        scale = T(1)
        im_min, im_max = floor(Int, minimum(A[mask])), ceil(Int, maximum(A[mask]))
    else
        im_min, im_max = floor(Int, minimum(A[mask]))*scale, ceil(Int, maximum(A[mask]))*scale
    end
    # scale min max of hist to hist_min or hist_max if given
    if hist_min != nothing
        im_min = hist_min * scale
    end
    if hist_max != nothing
        im_max = hist_max * scale
    end

    nodata_value_scaled = nodata_value*scale
    @assert (im_max - im_min) >= 1 "the image bit depth with given 'precision', 'nodata_value', 'hist_min' and 'hist_max' params is zero"
    # create histogram bin size
    nbins = ceil(Int, sqrt(im_max-im_min))
    # fine level histograms have bins with size 1
    # course bin histogram bins are defined by:
    edges_course = Int[im_min-1:nbins:nbins*nbins+im_min;]

    # initialization of histograms
    # the counts arrays are not setup as offsetarrays because slicing doesn't work on these arrays
    # instead the arrays have an offset of r in indices
    cnts_course = zeros(Int32, (nrows+2*r, nbins))
    cnts_fine = zeros(Int32, (nrows+2*r, nbins, nbins))
    # create seperate bins for both extremes and nodata
    cnts_nodata = zeros(Int32, nrows+2*r)
    cnts_min = zeros(Int32, nrows+2*r)
    cnts_max = zeros(Int32, nrows+2*r)

    # create output array B
    B = similar(A)
    # pad A
    Apad = padarray(A, Pad(:reflect, r, r))

    function add_value!(i::Int, v::Real; w=1)
        if scale > 1
            v = round(v * scale, 0) # scaling for non 8bit images
        end
        if v ≈ nodata_value_scaled
            cnts_nodata[i+r] += w # offset r in cnts arrays
        elseif v < im_min
            cnts_min[i+r] += w
        elseif v >= im_max
            cnts_max[i+r] += w
        else
            icourse = min(Int(div(v-im_min, nbins)) + 1, nbins) # min to include right boundary
            ifine = min(Int(div(v-edges_course[icourse], 1)), nbins) # fine bins have size 1
            cnts_course[i+r, icourse] += w # offset r in cnts arrays
            cnts_fine[i+r, ifine, icourse] += w
        end

        nothing
    end

    # find bin corresponding to percentile start from FIRST bin
    function bin_count_course(i::Int, p_idx::Int, pcum::Int)
        pbin, icourse = 0, 0
        for icourse in 1:nbins
            # use of sum is faster than loop. tested with BenchmarkTools and ProfileView
            pbin = sum(cnts_course[i:i+r*2, icourse])
            pcum += pbin #[icourse]
            pcum >= p_idx && break
        end
        pcum -= pbin
        pcum, icourse
    end
    function bin_count_fine(i::Int, p_idx::Int, pcum::Int, icourse::Int)
        ifine=0
        for ifine in 1:nbins
            pcum += sum(cnts_fine[i:i+r*2, ifine, icourse])
            pcum >= p_idx && break
        end
        pcum, ifine
    end

    # find bin corresponding to percentile -> start from LAST bin
    function bin_count_course_rev(i::Int, p_idx::Int, pcum::Int)
        pbin, icourse = 0, 0
        for icourse in nbins:-1:1
            pbin = sum(cnts_course[i:i+r*2, icourse])
            pcum += pbin #[icourse]
            pcum >= (ksize - p_idx) && break
        end
        pcum -= pbin
        pcum, icourse
    end
    function bin_count_fine_rev(i::Int, p_idx::Int, pcum::Int, icourse::Int)
        ifine = 0
        for ifine in nbins:-1:1
            pcum += sum(cnts_fine[i:i+r*2, ifine, icourse])
            pcum >= (ksize - p_idx) && break
        end
        pcum, ifine
    end

    # find percentile value from multi level histogram
    function calc_perc(i::Int)
        # number of nodata_size values in kernel
        nodata_size = sum(cnts_nodata[i:i+r*2])
        if nodata_size < ksize # calculate percentile if more than half window with values
            p_idx = ceil(Int, max((ksize-nodata_size) * perc / 100., 1))
            # find percentile bin
            icourse, ifine = 0, 0
            if perc <= 100
                hist_min != nothing ? pcum = sum(cnts_min[i:i+r*2]) : pcum = 0
                # start from begin of histogram
                pcum, icourse = bin_count_course(i, p_idx, pcum)
                pcum, ifine = bin_count_fine(i, p_idx, pcum, icourse)
            else
                hist_max != nothing ? pcum = sum(cnts_max[i:i+r*2]) : pcum = 0
                # start from end in histograms
                pcum, icourse = bin_count_course_rev(i, p_idx, pcum)
                pcum, ifine = bin_count_fine_rev(i, p_idx, pcum, icourse)
            end

            # calculate the value of fine bin at
            return (edges_course[icourse] + ifine) / scale
        else
            return nodata_value
        end
    end

    ## start actual looping
    for j = 1:ncols, i in 1:nrows
        if j == 1 && i == 1 # initiate top left
            for kj = -r:r, ki = -r:r
                add_value!(i+ki, Apad[i+ki, j+kj])
            end
        elseif j == 1 && i > 1 # initiate first col
            for kj in -r:r
                add_value!(i+r, Apad[i+r, j+kj])
            end
        elseif i == 1 && j > 1 # first row
            for ki in -r:r
                add_value!(i+ki, Apad[i+ki, j-r-1]; w=-1) # remove old value
                add_value!(i+ki, Apad[i+ki, j+r])
            end
        else
            add_value!(i+r, Apad[i+r, j-r-1]; w=-1) # remove old value
            add_value!(i+r, Apad[i+r, j+r])
        end
        B[i,j] = calc_perc(i)
    end
    B
end

"""apply vosselman criterium to grid,

vosselman criterium:
z_{i} <=  min z[ki] + dist[ki] * max_slope
for ki in i-r:i+r, j-r:j+r

output
flags   true if criterium not met
B       max allow value per cell given neighbors
"""
function vosselman_filter(A::Array{T,2},
    radius,             # kernel radius [m]
    max_slope,          # mas tolorated slope [m/m]
    tolerance=0;
    res=1.,             # resolution of grid a.k.a. cell size [m]
    nodata::T=T(-9999.0)) where T <: AbstractFloat

    r = round(Int, radius/res) # length in # off cells in which to search for neighbours

    # create output array B
    nrows, ncols = size(A)
    B = ones(T, (nrows, ncols)) * Inf16
    flags = falses(A)
    # pad A
    Apad = ImageFiltering.padarray(A, Pad(:replicate, r, r))

    # loop offsets
    for kj = -r:r, ki = -r:r
        k_dist = hypot(ki,kj)
        ((k_dist > r) || ((ki==1) && (kj==1))) && continue # skip offset if out of range or zero
        k_tol = k_dist .* res .* max_slope + tolerance
        ## loop through image
        for j = 1:ncols, i = 1:nrows
            zi = A[i,j]
            zi == nodata && continue
            # replace if smaller than current (find min in kernel)
            zj = Apad[i+ki,j+kj] + k_tol
            zj == nodata && continue  # skip nodata values in kernel
            if zj < B[i, j]
                B[i, j] = zj
            end
            if zj < zi
                flags[i, j] = true
            end
        end
    end

    B, flags
end

# morphological opening & tophat = img - opening(img)
# square window, no circular kernel yet available!
function morph_tophat(img::Matrix{Float32}, boundarymask::BitMatrix, r::Int)
    # img_min = mapwindow(minimum, img, (s, s)) # erosion
    # img_open = mapwindow(maximum, img_min, (s, s)) # dilation
    img_min = morph_erosion(img, boundarymask, r)
    img_open = morph_dilation(img_min, boundarymask, r)
    img_open
end

"""Implementation of the Progressive Morphological Filter (Zhang, 2003)
works based on zmin grid where dropouts are filled with nearest neighbor value
NOTE: make sure A has no nodata values

output
flags   true if criterium not met
B       max allow value per cell given neighbors

(Zhang, 2003): http://users.cis.fiu.edu/~chens/PDF/TGRS.pdf """
function pmf_filter(A::Matrix{Float32}, boundarymask::BitMatrix,
    max_window_radius::Float64, # max radius for squared window [m]
    slope::Float64,             # terrain slope [m/m]
    dh_max::Float64,            # maximum elevation threshold between subsequent tophat operations [m]
    dh_min::Float64,            # initial elevation threshold [m]
    res::Float64)               # cellsize [m]

    @show max_window_radius, slope, dh_min, dh_max, res

    nrow, ncol = size(A)
    totcells = nrow * ncol
    # maximum window size of kernel in pixels
    rad_pix_max = round(Int, max_window_radius / res)
    # exponential increase radius; window size = 3,5,9,17,33,..,N
    kmax = ceil(Int, log2(rad_pix_max*2))
    w = [Int(exp2(k)) + 1 for k in 1:kmax]
    w[end] = rad_pix_max * 2 + 1
    # deprecate linear window growth -> does not give satisfying results
    # linear increase radius; window size = 3,5,7,9...,N
    # w = Int[2*k + 1 for k in 1:rad_pix_max]

    # give a warning when the slope will exceed dh_max
    slope * (w[end] - w[end-1]) * res + dh_min > dh_max && warn("dh_max will limit pmf filter")

    # initialize
    dh_t = dh_min
    Af = copy(A)
    # create max allowed surface (use to classify pointcloud)
    B = Af + dh_t
    # loop through kernel sizes

    @showprogress 1 "Progessive filtering..." for (k, wk) in enumerate(w)

        # threshold
        if k > 1
            dh_t = min(slope * (wk - w[k-1]) * res + dh_min, dh_max)
        end

        # calculate morph opening
        Af = morph_tophat(Af, boundarymask, fld(wk, 2))

        # update max allowed surface
        if k == 1
            B[:,:] = Af + dh_t
        else
            for I in eachindex(Af)
                # replace if smaller than current z_max
                z_max = Af[I] + dh_t
                B[I] = min(B[I], z_max)
            end
        end

    end
    # flag non ground points based on dh_t threshold
    flags = Array{Bool}(A .> B)
    # ncells = count(x -> x == 1, flags)
    # debugging print statement to be placed in loop above
    # println("  iter $(k): $(ncells)/$(totcells) filtered, window size = $(wk) cells, dh = $(@sprintf("%.2f", dh_t))")

    B, flags
end

# lee filter
# http://stackoverflow.com/questions/39785970/speckle-lee-filter-in-python
# Lee, J. S. (1980). Digital image enhancement and noise filtering by use of
# local statistics. IEEE transactions on pattern analysis and machine intelligence, (2), 165-168.
#
# enl (equivalen number of looks) -> 1 in LiDAR intensity data
# addititive noise
function lee_filter(img::Array{T,2}, size::Real; nodata=T(-9999.0)) where T <: AbstractFloat
    # local statistics
    fmean(A::AbstractArray) = mean(filter(x -> x != nodata, A))
    window_mean = mapwindow(fmean, img, (size, size))
    window_sqr_mean = mapwindow(fmean, img.^2, (size, size))
    window_variance = window_sqr_mean - window_mean.^2
    # global variance
    img_variance = var(filter(x -> x != nodata, img))

    # smoothing weights
    w = window_variance.^2 ./ (window_variance.^2 + img_variance.^2)

    # smoothed image
    img_output = (img .* w) + (window_mean .* (1 - w))
end

# enhanced lee filter
#
# Lopes, A., Touzi, R., & Nezry, E. (1990). Adaptive speckle filters and scene
# heterogeneity. IEEE transactions on Geoscience and Remote Sensing, 28(6), 992-1000.
# enl (equivalent number of looks) -> 1.0 for lidar data
function enhanced_lee_filter(img::Array{T,2}, size::Real;
        enl=1.0,
        k=1.0,
        nodata=T(-9999.0)) where T <: AbstractFloat

    max_cof_var = sqrt(1+2/enl)
    cu = 0.523/sqrt(enl)

    # local statistics
    fmean(A::AbstractArray) = mean(filter(x -> x != nodata, A))
    window_mean = mapwindow(fmean, img, (size, size))
    window_sqr_mean = mapwindow(fmean, img.^2, (size, size))
    window_variance = window_sqr_mean - window_mean.^2
    # in case of img being Float32 and having very high values,
    # such as sometimes happens with intensity
    # slightly negative values for window_variance can arise
    # these are here set to 0.0
    window_variance = clamp.(window_variance, 0.0, Inf)
    window_cof_var = sqrt.(window_variance) ./ window_mean
    window_cof_var[window_cof_var .== 0] = 0.001 # avoid zero division

    # smoothing weights
    # filter in heterogeneous area
    w = T.(exp.(-k .* (window_cof_var - cu) ./ (max_cof_var - window_cof_var)))
    # preserver original value for "point target" (e.g around edges)
    w[window_cof_var .>= max_cof_var] = 0.0
    # return local mean in homogeneous areas
    w[window_cof_var .<= cu] = 1.0

    # smoothed image
    img_output = (img .* (1-w)) + (window_mean .* w)
end

# morphological operations with circular kernel
function morph_erosion(A::Matrix{Float32}, boundarymask::BitMatrix, r::Int)
    nrows, ncols = size(A) # image size
    B = copy(A) # create output array B
    Apad = padarray(A, Pad(:replicate, r, r)) # pad A

    # loop offsets
    @inbounds for kj = -r:r, ki = -r:r
        ((hypot(ki,kj) > r) || ((ki==0) && (kj==0))) && continue # skip offset if out of range or zero
        ## loop through image
        for j = 1:ncols, i = 1:nrows
            boundarymask[i,j] || continue
            # replace if smaller than current (find min in kernel)
            zo = Apad[i+ki,j+kj]
            if zo < B[i, j]
                B[i, j] = zo
            end
        end
    end
    B
end

function morph_dilation(A::Matrix{Float32}, boundarymask::BitMatrix, r::Int)
    nrows, ncols = size(A) # image size
    B = copy(A) # create output array B
    Apad = padarray(A, Pad(:replicate, r, r)) # pad A

    # loop offsets
    @inbounds for kj = -r:r, ki = -r:r
        ((hypot(ki,kj) > r) || ((ki==0) && (kj==0))) && continue # skip offset if out of range or zero
        ## loop through image
        for j = 1:ncols, i = 1:nrows
            boundarymask[i,j] || continue
            # replace if larger than current (find max in kernel)
            zo = Apad[i+ki,j+kj]
            if zo > B[i, j]
                B[i, j] = zo
            end
        end
    end
    B
end


text:
#####
##### GameInterface Tests
#####

# Generate a set of reachable states randomly
function generate_states(gspec, n)
  traces = [play_game(gspec, RandomPlayer()) for i in 1:n]
  return Set(s for t in traces for s in t.states)
end

same_state(game1, game2) = GI.current_state(game1) == GI.current_state(game2)

# Test that the same (symmetric) actions are still available in a symmetric state
# function test_symmetry(gspec, state, (symstate, aperm))
#   mask = GI.actions_mask(GI.init(gspec, state))
#   symmask = GI.actions_mask(GI.init(gspec, symstate))
#   v = falses(length(symmask))
#   v[mask] .= true
#   v = v[aperm]
#   return all(v[symmask]) && !any(v[.~symmask])
# end

function test_symmetry(gspec, state, (symstate, aperm))
  mask = GI.actions_mask(GI.init(gspec, state))
  symmask = GI.actions_mask(GI.init(gspec, symstate))
  return symmask == mask[aperm]
end

"""

    test_game(::AbstractGameSpec)
    test_game(experiment)

Perform some sanity checks regarding the compliance of a game with the
AlphaZero.jl [Game Interface](@ref game_interface).
"""
function test_game(gspec::AbstractGameSpec; n=100)
  @assert isa(GI.two_players(gspec), Bool)

  # Testing initial states
  if all(same_state(GI.init(gspec), GI.init(gspec)) for _ in 1:n)
    @info "There appears to be a unique initial state."
  else
    @info "There appears to be a nontrivial distribution of initial states."
  end
  @assert all(GI.white_playing(GI.init(gspec)) for _ in 1:n)

  # Save the game static properties
  State = GI.state_type(gspec)
  state_dim = GI.state_dim(gspec)
  Action = GI.action_type(gspec)
  actions = GI.actions(gspec)
  @assert isa(actions, AbstractVector{Action})
  @assert !isempty(actions)
  num_actions = GI.num_actions(gspec)
  two_players = GI.two_players(gspec)

  # Check properties on random states
  for state in generate_states(gspec, n)
    game = GI.init(gspec, state)
    @assert GI.current_state(game) == state
    GI.set_state!(game, state)
    @assert GI.current_state(game) == state
    @assert isa(GI.current_state(game), State)
    @assert same_state(GI.clone(game), GI.init(gspec, state))

    @assert isa(GI.white_playing(game), Bool)
    @assert GI.white_playing(game) || two_players

    @assert isa(GI.game_terminated(game), Bool)
    @assert GI.game_terminated(game) || !isempty(GI.available_actions(game))

    @assert eltype(GI.actions_mask(game)) == Bool
    @assert length(GI.actions_mask(game)) == num_actions

    vect = GI.vectorize_state(gspec, state)
    @assert isa(vect, Array{Float32})
    @assert size(vect) == state_dim

    # Playing a move
    if !GI.game_terminated(game)
      @assert isa(GI.heuristic_value(game), Float64)
      state_copy = deepcopy(state)
      a = rand(GI.available_actions(game))
      @assert isa(GI.action_string(gspec, a), String)
      GI.play!(game, a)
      @assert state_copy == state "States must appear as persistent."
      @assert isa(GI.white_reward(game), Float64)
    end

    # Testing symmetries
    syms = GI.symmetries(gspec, state)
    @assert isa(syms, Vector)
    for sym in syms
      s, σ = sym
      @assert isa(s, State)
      @assert isa(σ, Vector{Int})
      @assert length(σ) == num_actions
      @assert test_symmetry(gspec, state, sym)
    end

    # Verifying that static properties are state invariant
    @assert GI.num_actions(GI.spec(game)) == num_actions
    @assert GI.actions(GI.spec(game)) == actions
    @assert GI.state_type(GI.spec(game)) == State
    @assert GI.action_type(GI.spec(game)) == Action
    @assert GI.two_players(GI.spec(game)) == two_players
    @assert GI.state_dim(GI.spec(game)) == state_dim
  end
end


text:
sudoku25 = Dict(
:89 =>
    [
    11 23 13 10 19 16  6  2 24  7  5  9  1 20 17 15  8 18 25  3  4 12 21 22 14;
    15 16  0 22  0 11  8  0  0  0 25  0 14  0  0  0 12 19  0  0 17  0  0  0  0;
    0  0  0  0  0  0  0  0  0  0  0 16  0  4  0 17  0 13  0 24  0 23 19 10  2;
    0  0  0  0  0 19  0 14 23  4  0 21  6 22 10  0 11  0  2  0  0  0  0  0  0;
    17 14  0  0  2  0  0 13 12  0  0  0  0  0 15  4 20 22 10  0 11  0  9 24  8;
    22  0  0  0  0  6  2  0  0  0  4  7 12  1  9  0  0  0  0  0  0 14  5  0  0;
    0 18  2  0  8 22  0 19 16 21  0  0  0 10 13 23  0  0 20  0  0  3  0 15  7;
    0  0 17  3  0  5  0  0  8  9  0  0  0  0 18  0 19  0  0  0  0  0 23 21  0;
    1 11  0  0  9  0 15 10 25  0  6  0 23  0  0  0  0  5  3  7  0 17  0  0 24;
    0  0  0  0  0  0  1  0  0 23  0  0  0 24  0  0  0 21 12  0  6  8  0 25 16;
    20 24 10  0 15 23 11 17  0  0  0  0  0  7  0 12  0  0  0  0  0 22  0  0  6;
    4  5  0 14 12 25  0 18  0  0 23  0 15  0 19  1  0  0  0 22 20  0  7  9  0;
    18  0 21  0  0  8  0 24  0  0  9  0 25  0  0  0 10  0  0  0  2  0  1 19  0;
    0  0  6  2  1  0 13  0 22  0  0  0  0  0 11  8 21 16  0  0 25  0  0 12 17;
    0 17 25  0 23  7 14  0 21  1  0  0  0  0  3  0  0 11  0  0 24  0 16  4  5;
    0  0  0  0 11 18 24  0  0  0  0  5  0 12  0 25  0  0  0 15 23  4  8 14  0;
    0  0  0 15 21  0  0  0  0  0  2  0 13 17  0  0  1  7  0  0  5  9 24  0  0;
    0  0 18  0 22 15  0  0  2 16  0 23  0  0  0 10  6 24  0 17 12  0 25 11  0;
    7  2  0  1  0  0 21  0  0  0 18 22  0  9  6 14  0  4  5 16  0  0  0  0  0;
    0  0  9  0  0  0  7 22  0  0 10  0 24  0  0  0 18  0  0  0 21  0  0  0  0;
    0 12  0 19 10  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0 14  0  4  8  0;
    24  0 11 18  0  0  0  0  0  0  0 25 17 21  0  6  0  0  1  0  0  0  0  5 12;
    16  6 22  0  0  0 23  4 15 18  8  0  0  0 20  0  0 17  0 14  0  0  0  0  0;
    0 21  0  0  4  0  9  1  7  0  0  0  0 11 14  0 16  8 15  0 22  0 18  0  0;
    8 15  0  0  0  0  0  0  5  0 24  3  0  0  4  0  0  0  9  0  0  0  0  0 20;
    ],

:90 =>
    [
        0 23 13  0 19 16  6  0 24  7  5  9  1  0  0 15  8 18 25  0  4  0 21 22  0;
       15 16  0 22  0 11  8  0  0  0 25  0 14  0  0  0 12 19  0  0 17  0  0  0  0;
        0  0  0  0  0  0  0  0  0  0  0 16  0  4  0 17  0 13  0 24  0 23 19 10  2;
        0  0  0  0  0 19  0 14 23  4  0 21  6 22 10  0 11  0  2  0  0  0  0  0  0;
       17 14  0  0  2  0  0 13 12  0  0  0  0  0 15  4 20 22 10  0 11  0  9 24  8;
       22  0  0  0  0  6  2  0  0  0  4  7 12  1  9  0  0  0  0  0  0 14  5  0  0;
        0 18  2  0  8 22  0 19 16 21  0  0  0 10 13 23  0  0 20  0  0  3  0 15  7;
        0  0 17  3  0  5  0  0  8  9  0  0  0  0 18  0 19  0  0  0  0  0 23 21  0;
        1 11  0  0  9  0 15 10 25  0  6  0 23  0  0  0  0  5  3  7  0 17  0  0 24;
        0  0  0  0  0  0  1  0  0 23  0  0  0 24  0  0  0 21 12  0  6  8  0 25 16;
       20 24 10  0 15 23 11 17  0  0  0  0  0  7  0 12  0  0  0  0  0 22  0  0  6;
        4  5  0 14 12 25  0 18  0  0 23  0 15  0 19  1  0  0  0 22 20  0  7  9  0;
       18  0 21  0  0  8  0 24  0  0  9  0 25  0  0  0 10  0  0  0  2  0  1 19  0;
        0  0  6  2  1  0 13  0 22  0  0  0  0  0 11  8 21 16  0  0 25  0  0 12 17;
        0 17 25  0 23  7 14  0 21  1  0  0  0  0  3  0  0 11  0  0 24  0 16  4  5;
        0  0  0  0 11 18 24  0  0  0  0  5  0 12  0 25  0  0  0 15 23  4  8 14  0;
        0  0  0 15 21  0  0  0  0  0  2  0 13 17  0  0  1  7  0  0  5  9 24  0  0;
        0  0 18  0 22 15  0  0  2 16  0 23  0  0  0 10  6 24  0 17 12  0 25 11  0;
        7  2  0  1  0  0 21  0  0  0 18 22  0  9  6 14  0  4  5 16  0  0  0  0  0;
        0  0  9  0  0  0  7 22  0  0 10  0 24  0  0  0 18  0  0  0 21  0  0  0  0;
        0 12  0 19 10  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0 14  0  4  8  0;
       24  0 11 18  0  0  0  0  0  0  0 25 17 21  0  6  0  0  1  0  0  0  0  5 12;
       16  6 22  0  0  0 23  4 15 18  8  0  0  0 20  0  0 17  0 14  0  0  0  0  0;
        0 21  0  0  4  0  9  1  7  0  0  0  0 11 14  0 16  8 15  0 22  0 18  0  0;
        8 15  0  0  0  0  0  0  5  0 24  3  0  0  4  0  0  0  9  0  0  0  0  0 20;
       ]
)


text:
"""
    Monte-Carlo distributions
    =========================
"""

function mc_vel(num::Int64; v_s=1.95e3, delta_v=80.0, rho=2.5e-4, z=0.15, epsilon=1.0)
    """ A random sample of velocities with a mean longitudional velocity of v_s
        and a spread of delta_v.  The transverse velocity is found by assuming
        a point source and that the beam passes through a hole at z with a radius of rho
        and with x/y eccentricity of epsilon.
    """
    # longitudional beam velocity (Gaussian)
    vz = abs.(v_s .+ randn(num) .* (delta_v / (2.0^0.5)))
    # position at the laser
    r = rho .* rand(num).^0.5
    theta = 2.0 * pi * rand(num)
    # transverse velocity (determined by skimmer)
    tof = z ./ vz
    vx = epsilon * r .* cos.(theta) ./ tof
    vy = r .* sin.(theta) ./ tof
    return Array([vx vy vz])
end

function mc_flat(num, xmin=0.0, xmax=5.0e-6)
    """ Random values, evenly sampled from xmin < x < xmax.
    """
    x = xmin .+ rand(num) .* (xmax - xmin)
    return x
end

function mc_uniform(num::Int64; tmin=0.0, tmax=5e-6, distance=0.3, z0=0.198, reverse_vz=false,
                                v_s=1.86e3, delta_v=65.0, rho=1.5e-4, epsilon=1.0)
    """ 3D He distribution with a flat, oval distribution in xy.
    """
    t =  mc_flat(num, tmin, tmax)
    vel = mc_vel(num, v_s=v_s, delta_v=delta_v, rho=rho, z=distance, epsilon=epsilon)
    tof = distance ./ vel[:, 3]
    x = tof .* vel[:, 1]
    y = tof .* vel[:, 2]
    z = zeros(num) .+ z0
    if reverse_vz
    vel[:, 3] = - vel[:, 3]
    end
    df = DataFrame(hcat(t, x, y, z, vel))
    rename!(df, [:t, :x, :y, :z, :vx, :vy, :vz]);
    return df
end

function mc_supersonic(num::Int64; tmin=0.0, tmax=5e-6, distance=0.2, z0=0., reverse_vz=false,
                                   v_s=2.0e3, delta_v=65.0, sigma_x=5e-4, sigma_y=1e-4, rho=1e-4)
    """ A random sample of velocities with a mean longitudional velocity of v_s
        and a spread of delta_v.  The transverse velocity is found by assuming
        a small source (rho) located at z=(z0 - distance) and that the beam has a
        2D Gaussian distribution in xy.
    """
    # time
    t =  mc_flat(num, tmin, tmax)
    # longitudional velocity
    vz = abs.(v_s .+ randn(num) .* (delta_v / (2.0^0.5)))
    tof = distance ./ vz
    x = randn(num) .* sigma_x
    y = randn(num) .* sigma_y
    # position at the source
    r = rho .* rand(num).^0.5
    theta = 2.0 * pi * rand(num)
    x0 = r .* sin.(theta)
    y0 = r .* cos.(theta)
    # vel
    vx = (x .- x0) ./ tof
    vy = (y .- y0) ./ tof
    z = zeros(num) .+ z0
    if reverse_vz
        vz = - vz
    end
    df = DataFrame(hcat(t, x, y, z, vx, vy, vz))
    rename!(df, [:t, :x, :y, :z, :vx, :vy, :vz]);
    return df
end


text:
function add_load!(ls::LoopSet, op::Operation, actualarray::Bool = true, broadcast::Bool = false)
    @assert isload(op)
    ref = op.ref
    id = findfirst(r -> r == ref, ls.refs_aliasing_syms)
    # try to CSE
    if id === nothing
        push!(ls.syms_aliasing_refs, name(op))
        push!(ls.refs_aliasing_syms, ref)
    else
        opp = ls.opdict[ls.syms_aliasing_refs[id]] # throw an error if not found.
        return isstore(opp) ? getop(ls, first(parents(opp))) : opp
    end
    add_vptr!(ls, op.ref.ref.array, vptr(op), actualarray, broadcast)
    pushop!(ls, op, name(op))
end

function add_load!(
    ls::LoopSet, var::Symbol, array::Symbol, rawindices, elementbytes::Int
)
    mpref = array_reference_meta!(ls, array, rawindices, elementbytes, var)
    add_load!(ls, mpref, elementbytes)
end
function add_load!(
    ls::LoopSet, mpref::ArrayReferenceMetaPosition, elementbytes::Int
)
    length(mpref.loopdependencies) == 0 && return add_constant!(ls, mpref, elementbytes)
    op = Operation( ls, varname(mpref), elementbytes, :getindex, memload, mpref )
    add_load!(ls, op, true, false)
end

# for use with broadcasting
function add_simple_load!(
    ls::LoopSet, var::Symbol, ref::ArrayReference, elementbytes::Int,
    actualarray::Bool = true, broadcast::Bool = false
)
    loopdeps = Symbol[s for s ∈ ref.indices]
    mref = ArrayReferenceMeta(
        ref, fill(true, length(loopdeps) - isdiscontiguous(ref))
    )
    add_simple_load!(ls, var, mref, loopdeps, elementbytes, actualarray, broadcast)
end
function add_simple_load!(
    ls::LoopSet, var::Symbol, mref::ArrayReferenceMeta, loopdeps::Vector{Symbol},
    elementbytes::Int, actualarray::Bool = true, broadcast::Bool = false
)
    op = Operation(
        length(operations(ls)), var, elementbytes,
        :getindex, memload, loopdeps,
        NODEPENDENCY, NOPARENTS, mref
    )
    add_vptr!(ls, op.ref.ref.array, vptr(op), actualarray, broadcast)
    pushop!(ls, op, var)
end
function add_load_ref!(ls::LoopSet, var::Symbol, ex::Expr, elementbytes::Int)
    array, rawindices = ref_from_ref!(ls, ex)
    add_load!(ls, var, array, rawindices, elementbytes)
end
function add_load_getindex!(ls::LoopSet, var::Symbol, ex::Expr, elementbytes::Int)
    array, rawindices = ref_from_getindex!(ls, ex)
    add_load!(ls, var, array, rawindices, elementbytes)
end

function add_loopvalue!(ls::LoopSet, arg::Symbol, elementbytes::Int)
    # check for CSE opportunity
    instr = Instruction(arg, arg)
    for op ∈ operations(ls)#check to CSE
        (op.variable === arg && instr == instruction(op)) && return op
    end
    op = Operation(
        length(operations(ls)), arg, elementbytes, instr, loopvalue, [arg]
    )
    pushop!(ls, op, arg)
end