text:
function calculate_tgrad(sys::AbstractODESystem;
                         simplify = false)
    isempty(get_tgrad(sys)[]) || return get_tgrad(sys)[]  # use cached tgrad, if possible

    # We need to remove explicit time dependence on the state because when we
    # have `u(t) * t` we want to have the tgrad to be `u(t)` instead of `u'(t) *
    # t + u(t)`.
    rhs = [detime_dvs(eq.rhs) for eq in full_equations(sys)]
    iv = get_iv(sys)
    xs = states(sys)
    rule = Dict(map((x, xt) -> xt => x, detime_dvs.(xs), xs))
    rhs = substitute.(rhs, Ref(rule))
    tgrad = [expand_derivatives(Differential(iv)(r), simplify) for r in rhs]
    reverse_rule = Dict(map((x, xt) -> x => xt, detime_dvs.(xs), xs))
    tgrad = Num.(substitute.(tgrad, Ref(reverse_rule)))
    get_tgrad(sys)[] = tgrad
    return tgrad
end

function calculate_jacobian(sys::AbstractODESystem;
                            sparse = false, simplify = false, dvs = states(sys))
    if isequal(dvs, states(sys))
        cache = get_jac(sys)[]
        if cache isa Tuple && cache[2] == (sparse, simplify)
            return cache[1]
        end
    end

    rhs = [eq.rhs - eq.lhs for eq in full_equations(sys)] #need du terms on rhs for differentiating wrt du

    iv = get_iv(sys)

    if sparse
        jac = sparsejacobian(rhs, dvs, simplify = simplify)
    else
        jac = jacobian(rhs, dvs, simplify = simplify)
    end

    if isequal(dvs, states(sys))
        get_jac(sys)[] = jac, (sparse, simplify) # cache Jacobian
    end

    return jac
end

function calculate_control_jacobian(sys::AbstractODESystem;
                                    sparse = false, simplify = false)
    cache = get_ctrl_jac(sys)[]
    if cache isa Tuple && cache[2] == (sparse, simplify)
        return cache[1]
    end

    rhs = [eq.rhs for eq in full_equations(sys)]

    iv = get_iv(sys)
    ctrls = controls(sys)

    if sparse
        jac = sparsejacobian(rhs, ctrls, simplify = simplify)
    else
        jac = jacobian(rhs, ctrls, simplify = simplify)
    end

    get_ctrl_jac(sys)[] = jac, (sparse, simplify) # cache Jacobian
    return jac
end

function generate_tgrad(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
                        simplify = false, kwargs...)
    tgrad = calculate_tgrad(sys, simplify = simplify)
    return build_function(tgrad, dvs, ps, get_iv(sys); kwargs...)
end

function generate_jacobian(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
                           simplify = false, sparse = false, kwargs...)
    jac = calculate_jacobian(sys; simplify = simplify, sparse = sparse)
    return build_function(jac, dvs, ps, get_iv(sys); kwargs...)
end

function generate_control_jacobian(sys::AbstractODESystem, dvs = states(sys),
                                   ps = parameters(sys);
                                   simplify = false, sparse = false, kwargs...)
    jac = calculate_control_jacobian(sys; simplify = simplify, sparse = sparse)
    return build_function(jac, dvs, ps, get_iv(sys); kwargs...)
end

function generate_dae_jacobian(sys::AbstractODESystem, dvs = states(sys),
                               ps = parameters(sys); simplify = false, sparse = false,
                               kwargs...)
    jac_u = calculate_jacobian(sys; simplify = simplify, sparse = sparse)
    derivatives = Differential(get_iv(sys)).(states(sys))
    jac_du = calculate_jacobian(sys; simplify = simplify, sparse = sparse,
                                dvs = derivatives)
    dvs = states(sys)
    @variables ˍ₋gamma
    jac = ˍ₋gamma * jac_du + jac_u
    return build_function(jac, derivatives, dvs, ps, ˍ₋gamma, get_iv(sys); kwargs...)
end

function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
                           implicit_dae = false,
                           ddvs = implicit_dae ? map(Differential(get_iv(sys)), dvs) :
                                  nothing,
                           has_difference = false,
                           kwargs...)
    eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
    if !implicit_dae
        check_operator_variables(eqs, Differential)
        check_lhs(eqs, Differential, Set(dvs))
    end
    # substitute x(t) by just x
    rhss = implicit_dae ? [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] :
           [eq.rhs for eq in eqs]

    # TODO: add an optional check on the ordering of observed equations
    u = map(x -> time_varying_as_func(value(x), sys), dvs)
    p = map(x -> time_varying_as_func(value(x), sys), ps)
    t = get_iv(sys)

    pre, sol_states = get_substitutions_and_solved_states(sys,
                                                          no_postprocess = has_difference)

    if implicit_dae
        build_function(rhss, ddvs, u, p, t; postprocess_fbody = pre, states = sol_states,
                       kwargs...)
    else
        build_function(rhss, u, p, t; postprocess_fbody = pre, states = sol_states,
                       kwargs...)
    end
end

function generate_difference_cb(sys::ODESystem, dvs = states(sys), ps = parameters(sys);
                                kwargs...)
    eqs = equations(sys)
    check_operator_variables(eqs, Difference)

    var2eq = Dict(arguments(eq.lhs)[1] => eq for eq in eqs if isdifference(eq.lhs))

    u = map(x -> time_varying_as_func(value(x), sys), dvs)
    p = map(x -> time_varying_as_func(value(x), sys), ps)
    t = get_iv(sys)

    body = map(dvs) do v
        eq = get(var2eq, v, nothing)
        eq === nothing && return v
        d = operation(eq.lhs)
        d.update ? eq.rhs : eq.rhs + v
    end

    pre = get_postprocess_fbody(sys)
    f_oop, f_iip = build_function(body, u, p, t; expression = Val{false},
                                  postprocess_fbody = pre, kwargs...)

    cb_affect! = let f_oop = f_oop, f_iip = f_iip
        function cb_affect!(integ)
            if DiffEqBase.isinplace(integ.sol.prob)
                tmp, = DiffEqBase.get_tmp_cache(integ)
                f_iip(tmp, integ.u, integ.p, integ.t) # aliasing `integ.u` would be bad.
                copyto!(integ.u, tmp)
            else
                integ.u = f_oop(integ.u, integ.p, integ.t)
            end
            return nothing
        end
    end

    getdt(eq) = operation(eq.lhs).dt
    deqs = values(var2eq)
    dt = getdt(first(deqs))
    all(dt == getdt(eq) for eq in deqs) ||
        error("All difference variables should have same time steps.")

    PeriodicCallback(cb_affect!, first(dt))
end

function generate_rootfinding_callback(sys::ODESystem, dvs = states(sys),
                                       ps = parameters(sys); kwargs...)
    cbs = continuous_events(sys)
    isempty(cbs) && return nothing
    generate_rootfinding_callback(cbs, sys, dvs, ps; kwargs...)
end

function generate_rootfinding_callback(cbs, sys::ODESystem, dvs = states(sys),
                                       ps = parameters(sys); kwargs...)
    eqs = map(cb -> cb.eqs, cbs)
    num_eqs = length.(eqs)
    (isempty(eqs) || sum(num_eqs) == 0) && return nothing
    # fuse equations to create VectorContinuousCallback
    eqs = reduce(vcat, eqs)
    # rewrite all equations as 0 ~ interesting stuff
    eqs = map(eqs) do eq
        isequal(eq.lhs, 0) && return eq
        0 ~ eq.lhs - eq.rhs
    end

    rhss = map(x -> x.rhs, eqs)
    root_eq_vars = unique(collect(Iterators.flatten(map(ModelingToolkit.vars, rhss))))

    u = map(x -> time_varying_as_func(value(x), sys), dvs)
    p = map(x -> time_varying_as_func(value(x), sys), ps)
    t = get_iv(sys)
    rf_oop, rf_ip = build_function(rhss, u, p, t; expression = Val{false}, kwargs...)

    affect_functions = map(cbs) do cb # Keep affect function separate
        eq_aff = affect_equations(cb)
        affect = compile_affect(eq_aff, sys, dvs, ps; kwargs...)
    end

    if length(eqs) == 1
        cond = function (u, t, integ)
            if DiffEqBase.isinplace(integ.sol.prob)
                tmp, = DiffEqBase.get_tmp_cache(integ)
                rf_ip(tmp, u, integ.p, t)
                tmp[1]
            else
                rf_oop(u, integ.p, t)
            end
        end
        ContinuousCallback(cond, affect_functions[])
    else
        cond = function (out, u, t, integ)
            rf_ip(out, u, integ.p, t)
        end

        # since there may be different number of conditions and affects,
        # we build a map that translates the condition eq. number to the affect number
        eq_ind2affect = reduce(vcat,
                               [fill(i, num_eqs[i]) for i in eachindex(affect_functions)])
        @assert length(eq_ind2affect) == length(eqs)
        @assert maximum(eq_ind2affect) == length(affect_functions)

        affect = let affect_functions = affect_functions, eq_ind2affect = eq_ind2affect
            function (integ, eq_ind) # eq_ind refers to the equation index that triggered the event, each event has num_eqs[i] equations
                affect_functions[eq_ind2affect[eq_ind]](integ)
            end
        end
        VectorContinuousCallback(cond, affect, length(eqs))
    end
end

function compile_affect(cb::SymbolicContinuousCallback, args...; kwargs...)
    compile_affect(affect_equations(cb), args...; kwargs...)
end

"""
    compile_affect(eqs::Vector{Equation}, sys, dvs, ps; kwargs...)
    compile_affect(cb::SymbolicContinuousCallback, args...; kwargs...)

Returns a function that takes an integrator as argument and modifies the state with the affect.
"""
function compile_affect(eqs::Vector{Equation}, sys, dvs, ps; kwargs...)
    if isempty(eqs)
        return (args...) -> () # We don't do anything in the callback, we're just after the event
    else
        rhss = map(x -> x.rhs, eqs)
        lhss = map(x -> x.lhs, eqs)
        update_vars = collect(Iterators.flatten(map(ModelingToolkit.vars, lhss))) # these are the ones we're chaning
        length(update_vars) == length(unique(update_vars)) == length(eqs) ||
            error("affected variables not unique, each state can only be affected by one equation for a single `root_eqs => affects` pair.")
        vars = states(sys)

        u = map(x -> time_varying_as_func(value(x), sys), vars)
        p = map(x -> time_varying_as_func(value(x), sys), ps)
        t = get_iv(sys)
        rf_oop, rf_ip = build_function(rhss, u, p, t; expression = Val{false}, kwargs...)

        stateind(sym) = findfirst(isequal(sym), vars)

        update_inds = stateind.(update_vars)
        let update_inds = update_inds
            function (integ)
                lhs = @views integ.u[update_inds]
                rf_ip(lhs, integ.u, integ.p, integ.t)
            end
        end
    end
end

function time_varying_as_func(x, sys::AbstractTimeDependentSystem)
    # if something is not x(t) (the current state)
    # but is `x(t-1)` or something like that, pass in `x` as a callable function rather
    # than pass in a value in place of x(t).
    #
    # This is done by just making `x` the argument of the function.
    if istree(x) &&
       operation(x) isa Sym &&
       !(length(arguments(x)) == 1 && isequal(arguments(x)[1], get_iv(sys)))
        return operation(x)
    end
    return x
end

function calculate_massmatrix(sys::AbstractODESystem; simplify = false)
    eqs = [eq for eq in full_equations(sys) if !isdifferenceeq(eq)]
    dvs = states(sys)
    M = zeros(length(eqs), length(eqs))
    state2idx = Dict(s => i for (i, s) in enumerate(dvs))
    for (i, eq) in enumerate(eqs)
        if eq.lhs isa Term && operation(eq.lhs) isa Differential
            st = var_from_nested_derivative(eq.lhs)[1]
            j = state2idx[st]
            M[i, j] = 1
        else
            _iszero(eq.lhs) ||
                error("Only semi-explicit constant mass matrices are currently supported. Faulty equation: $eq.")
        end
    end
    M = simplify ? ModelingToolkit.simplify.(M) : M
    # M should only contain concrete numbers
    M == I ? I : M
end

function jacobian_sparsity(sys::AbstractODESystem)
    sparsity = torn_system_jacobian_sparsity(sys)
    sparsity === nothing || return sparsity

    jacobian_sparsity([eq.rhs for eq in full_equations(sys)],
                      [dv for dv in states(sys)])
end

function jacobian_dae_sparsity(sys::AbstractODESystem)
    J1 = jacobian_sparsity([eq.rhs for eq in full_equations(sys)],
                           [dv for dv in states(sys)])
    derivatives = Differential(get_iv(sys)).(states(sys))
    J2 = jacobian_sparsity([eq.rhs for eq in full_equations(sys)],
                           [dv for dv in derivatives])
    J1 + J2
end

function isautonomous(sys::AbstractODESystem)
    tgrad = calculate_tgrad(sys; simplify = true)
    all(iszero, tgrad)
end

for F in [:ODEFunction, :DAEFunction]
    @eval function DiffEqBase.$F(sys::AbstractODESystem, args...; kwargs...)
        $F{true}(sys, args...; kwargs...)
    end
end

"""
```julia
function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create an `ODEFunction` from the [`ODESystem`](@ref). The arguments `dvs` and `ps`
are used to set the order of the dependent variable and parameter vectors,
respectively.
"""
function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys), u0 = nothing;
                                     version = nothing, tgrad = false,
                                     jac = false,
                                     eval_expression = true,
                                     sparse = false, simplify = false,
                                     eval_module = @__MODULE__,
                                     steady_state = false,
                                     checkbounds = false,
                                     sparsity = false,
                                     kwargs...) where {iip}
    f_gen = generate_function(sys, dvs, ps; expression = Val{eval_expression},
                              expression_module = eval_module, checkbounds = checkbounds,
                              kwargs...)
    f_oop, f_iip = eval_expression ?
                   (@RuntimeGeneratedFunction(eval_module, ex) for ex in f_gen) : f_gen
    f(u, p, t) = f_oop(u, p, t)
    f(du, u, p, t) = f_iip(du, u, p, t)

    if tgrad
        tgrad_gen = generate_tgrad(sys, dvs, ps;
                                   simplify = simplify,
                                   expression = Val{eval_expression},
                                   expression_module = eval_module,
                                   checkbounds = checkbounds, kwargs...)
        tgrad_oop, tgrad_iip = eval_expression ?
                               (@RuntimeGeneratedFunction(eval_module, ex) for ex in tgrad_gen) :
                               tgrad_gen
        _tgrad(u, p, t) = tgrad_oop(u, p, t)
        _tgrad(J, u, p, t) = tgrad_iip(J, u, p, t)
    else
        _tgrad = nothing
    end

    if jac
        jac_gen = generate_jacobian(sys, dvs, ps;
                                    simplify = simplify, sparse = sparse,
                                    expression = Val{eval_expression},
                                    expression_module = eval_module,
                                    checkbounds = checkbounds, kwargs...)
        jac_oop, jac_iip = eval_expression ?
                           (@RuntimeGeneratedFunction(eval_module, ex) for ex in jac_gen) :
                           jac_gen
        _jac(u, p, t) = jac_oop(u, p, t)
        _jac(J, u, p, t) = jac_iip(J, u, p, t)
    else
        _jac = nothing
    end

    M = calculate_massmatrix(sys)

    _M = if sparse && !(u0 === nothing || M === I)
        SparseArrays.sparse(M)
    elseif u0 === nothing || M === I
        M
    else
        ArrayInterfaceCore.restructure(u0 .* u0', M)
    end

    obs = observed(sys)
    observedfun = if steady_state
        let sys = sys, dict = Dict()
            function generated_observed(obsvar, u, p, t = Inf)
                obs = get!(dict, value(obsvar)) do
                    build_explicit_observed_function(sys, obsvar)
                end
                obs(u, p, t)
            end
        end
    else
        let sys = sys, dict = Dict()
            function generated_observed(obsvar, u, p, t)
                obs = get!(dict, value(obsvar)) do
                    build_explicit_observed_function(sys, obsvar; checkbounds = checkbounds)
                end
                obs(u, p, t)
            end
        end
    end

    jac_prototype = if sparse
        uElType = u0 === nothing ? Float64 : eltype(u0)
        if jac
            similar(calculate_jacobian(sys, sparse = sparse), uElType)
        else
            similar(jacobian_sparsity(sys), uElType)
        end
    else
        nothing
    end
    ODEFunction{iip}(f,
                     jac = _jac === nothing ? nothing : _jac,
                     tgrad = _tgrad === nothing ? nothing : _tgrad,
                     mass_matrix = _M,
                     jac_prototype = jac_prototype,
                     syms = Symbol.(states(sys)),
                     indepsym = Symbol(get_iv(sys)),
                     observed = observedfun,
                     sparsity = sparsity ? jacobian_sparsity(sys) : nothing)
end

"""
```julia
function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create an `DAEFunction` from the [`ODESystem`](@ref). The arguments `dvs` and
`ps` are used to set the order of the dependent variable and parameter vectors,
respectively.
"""
function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys), u0 = nothing;
                                     ddvs = map(diff2term ∘ Differential(get_iv(sys)), dvs),
                                     version = nothing,
                                     jac = false,
                                     eval_expression = true,
                                     sparse = false, simplify = false,
                                     eval_module = @__MODULE__,
                                     checkbounds = false,
                                     kwargs...) where {iip}
    f_gen = generate_function(sys, dvs, ps; implicit_dae = true,
                              expression = Val{eval_expression},
                              expression_module = eval_module, checkbounds = checkbounds,
                              kwargs...)
    f_oop, f_iip = eval_expression ?
                   (@RuntimeGeneratedFunction(eval_module, ex) for ex in f_gen) : f_gen
    f(du, u, p, t) = f_oop(du, u, p, t)
    f(out, du, u, p, t) = f_iip(out, du, u, p, t)

    if jac
        jac_gen = generate_dae_jacobian(sys, dvs, ps;
                                        simplify = simplify, sparse = sparse,
                                        expression = Val{eval_expression},
                                        expression_module = eval_module,
                                        checkbounds = checkbounds, kwargs...)
        jac_oop, jac_iip = eval_expression ?
                           (@RuntimeGeneratedFunction(eval_module, ex) for ex in jac_gen) :
                           jac_gen
        _jac(du, u, p, ˍ₋gamma, t) = jac_oop(du, u, p, ˍ₋gamma, t)

        _jac(J, du, u, p, ˍ₋gamma, t) = jac_iip(J, du, u, p, ˍ₋gamma, t)
    else
        _jac = nothing
    end

    jac_prototype = if sparse
        uElType = u0 === nothing ? Float64 : eltype(u0)
        if jac
            J1 = calculate_jacobian(sys, sparse = sparse)
            derivatives = Differential(get_iv(sys)).(states(sys))
            J2 = calculate_jacobian(sys; sparse = sparse, dvs = derivatives)
            similar(J1 + J2, uElType)
        else
            similar(jacobian_dae_sparsity(sys), uElType)
        end
    else
        nothing
    end

    DAEFunction{iip}(f,
                     jac = _jac === nothing ? nothing : _jac,
                     syms = Symbol.(dvs),
                     jac_prototype = jac_prototype
                     # missing fields in `DAEFunction`
                     #indepsym = Symbol(get_iv(sys)),
                     #observed = observedfun,
                     )
end

"""
```julia
function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create a Julia expression for an `ODEFunction` from the [`ODESystem`](@ref).
The arguments `dvs` and `ps` are used to set the order of the dependent
variable and parameter vectors, respectively.
"""
struct ODEFunctionExpr{iip} end

struct ODEFunctionClosure{O, I} <: Function
    f_oop::O
    f_iip::I
end
(f::ODEFunctionClosure)(u, p, t) = f.f_oop(u, p, t)
(f::ODEFunctionClosure)(du, u, p, t) = f.f_iip(du, u, p, t)

function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                              ps = parameters(sys), u0 = nothing;
                              version = nothing, tgrad = false,
                              jac = false,
                              linenumbers = false,
                              sparse = false, simplify = false,
                              steady_state = false,
                              kwargs...) where {iip}
    f_oop, f_iip = generate_function(sys, dvs, ps; expression = Val{true}, kwargs...)

    dict = Dict()

    fsym = gensym(:f)
    _f = :($fsym = $ODEFunctionClosure($f_oop, $f_iip))
    tgradsym = gensym(:tgrad)
    if tgrad
        tgrad_oop, tgrad_iip = generate_tgrad(sys, dvs, ps;
                                              simplify = simplify,
                                              expression = Val{true}, kwargs...)
        _tgrad = :($tgradsym = $ODEFunctionClosure($tgrad_oop, $tgrad_iip))
    else
        _tgrad = :($tgradsym = nothing)
    end

    jacsym = gensym(:jac)
    if jac
        jac_oop, jac_iip = generate_jacobian(sys, dvs, ps;
                                             sparse = sparse, simplify = simplify,
                                             expression = Val{true}, kwargs...)
        _jac = :($jacsym = $ODEFunctionClosure($jac_oop, $jac_iip))
    else
        _jac = :($jacsym = nothing)
    end

    M = calculate_massmatrix(sys)

    _M = if sparse && !(u0 === nothing || M === I)
        SparseArrays.sparse(M)
    elseif u0 === nothing || M === I
        M
    else
        ArrayInterfaceCore.restructure(u0 .* u0', M)
    end

    jp_expr = sparse ? :(similar($(get_jac(sys)[]), Float64)) : :nothing
    ex = quote
        $_f
        $_tgrad
        $_jac
        M = $_M
        ODEFunction{$iip}($fsym,
                          jac = $jacsym,
                          tgrad = $tgradsym,
                          mass_matrix = M,
                          jac_prototype = $jp_expr,
                          syms = $(Symbol.(states(sys))),
                          indepsym = $(QuoteNode(Symbol(get_iv(sys)))))
    end
    !linenumbers ? striplines(ex) : ex
end

function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
                           implicit_dae = false, du0map = nothing,
                           version = nothing, tgrad = false,
                           jac = false,
                           checkbounds = false, sparse = false,
                           simplify = false,
                           linenumbers = true, parallel = SerialForm(),
                           eval_expression = true,
                           use_union = false,
                           kwargs...)
    eqs = equations(sys)
    dvs = states(sys)
    ps = parameters(sys)
    iv = get_iv(sys)

    defs = defaults(sys)
    defs = mergedefaults(defs, parammap, ps)
    defs = mergedefaults(defs, u0map, dvs)

    u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = true)
    p = varmap_to_vars(parammap, ps; defaults = defs, tofloat = !use_union, use_union)
    if implicit_dae && du0map !== nothing
        ddvs = map(Differential(iv), dvs)
        defs = mergedefaults(defs, du0map, ddvs)
        du0 = varmap_to_vars(du0map, ddvs; defaults = defs, toterm = identity,
                             tofloat = true)
    else
        du0 = nothing
        ddvs = nothing
    end

    check_eqs_u0(eqs, dvs, u0; kwargs...)

    f = constructor(sys, dvs, ps, u0; ddvs = ddvs, tgrad = tgrad, jac = jac,
                    checkbounds = checkbounds,
                    linenumbers = linenumbers, parallel = parallel, simplify = simplify,
                    sparse = sparse, eval_expression = eval_expression, kwargs...)
    implicit_dae ? (f, du0, u0, p) : (f, u0, p)
end

function ODEFunctionExpr(sys::AbstractODESystem, args...; kwargs...)
    ODEFunctionExpr{true}(sys, args...; kwargs...)
end

"""
```julia
function DAEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create a Julia expression for an `ODEFunction` from the [`ODESystem`](@ref).
The arguments `dvs` and `ps` are used to set the order of the dependent
variable and parameter vectors, respectively.
"""
struct DAEFunctionExpr{iip} end

struct DAEFunctionClosure{O, I} <: Function
    f_oop::O
    f_iip::I
end
(f::DAEFunctionClosure)(du, u, p, t) = f.f_oop(du, u, p, t)
(f::DAEFunctionClosure)(out, du, u, p, t) = f.f_iip(out, du, u, p, t)

function DAEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                              ps = parameters(sys), u0 = nothing;
                              version = nothing, tgrad = false,
                              jac = false,
                              linenumbers = false,
                              sparse = false, simplify = false,
                              kwargs...) where {iip}
    f_oop, f_iip = generate_function(sys, dvs, ps; expression = Val{true},
                                     implicit_dae = true, kwargs...)
    fsym = gensym(:f)
    _f = :($fsym = $DAEFunctionClosure($f_oop, $f_iip))
    ex = quote
        $_f
        ODEFunction{$iip}($fsym)
    end
    !linenumbers ? striplines(ex) : ex
end

function DAEFunctionExpr(sys::AbstractODESystem, args...; kwargs...)
    DAEFunctionExpr{true}(sys, args...; kwargs...)
end

for P in [:ODEProblem, :DAEProblem]
    @eval function DiffEqBase.$P(sys::AbstractODESystem, args...; kwargs...)
        $P{true}(sys, args...; kwargs...)
    end
end

"""
```julia
function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    simplify=false,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```

Generates an ODEProblem from an ODESystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem, u0map, tspan,
                                    parammap = DiffEqBase.NullParameters();
                                    callback = nothing,
                                    check_length = true, kwargs...) where {iip}
    has_difference = any(isdifferenceeq, equations(sys))
    f, u0, p = process_DEProblem(ODEFunction{iip}, sys, u0map, parammap;
                                 has_difference = has_difference,
                                 check_length, kwargs...)
    if has_continuous_events(sys)
        event_cb = generate_rootfinding_callback(sys; kwargs...)
    else
        event_cb = nothing
    end
    difference_cb = has_difference ? generate_difference_cb(sys; kwargs...) : nothing
    cb = merge_cb(event_cb, difference_cb)
    cb = merge_cb(cb, callback)

    if cb === nothing
        ODEProblem{iip}(f, u0, tspan, p; kwargs...)
    else
        ODEProblem{iip}(f, u0, tspan, p; callback = cb, kwargs...)
    end
end
merge_cb(::Nothing, ::Nothing) = nothing
merge_cb(::Nothing, x) = merge_cb(x, nothing)
merge_cb(x, ::Nothing) = x
merge_cb(x, y) = CallbackSet(x, y)

"""
```julia
function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem,du0map,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    simplify=false,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```

Generates an DAEProblem from an ODESystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
                                    parammap = DiffEqBase.NullParameters();
                                    check_length = true, kwargs...) where {iip}
    has_difference = any(isdifferenceeq, equations(sys))
    f, du0, u0, p = process_DEProblem(DAEFunction{iip}, sys, u0map, parammap;
                                      implicit_dae = true, du0map = du0map,
                                      has_difference = has_difference, check_length,
                                      kwargs...)
    diffvars = collect_differential_variables(sys)
    sts = states(sys)
    differential_vars = map(Base.Fix2(in, diffvars), sts)
    if has_difference
        DAEProblem{iip}(f, du0, u0, tspan, p;
                        difference_cb = generate_difference_cb(sys; kwargs...),
                        differential_vars = differential_vars, kwargs...)
    else
        DAEProblem{iip}(f, du0, u0, tspan, p; differential_vars = differential_vars,
                        kwargs...)
    end
end

"""
```julia
function ODEProblemExpr{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, parallel=SerialForm(),
                                    skipzeros=true, fillzeros=true,
                                    simplify=false,
                                    kwargs...) where iip
```

Generates a Julia expression for constructing an ODEProblem from an
ODESystem and allows for automatically symbolically calculating
numerical enhancements.
"""
struct ODEProblemExpr{iip} end

function ODEProblemExpr{iip}(sys::AbstractODESystem, u0map, tspan,
                             parammap = DiffEqBase.NullParameters(); check_length = true,
                             kwargs...) where {iip}
    f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap; check_length,
                                 kwargs...)
    linenumbers = get(kwargs, :linenumbers, true)

    ex = quote
        f = $f
        u0 = $u0
        tspan = $tspan
        p = $p
        ODEProblem(f, u0, tspan, p; $(kwargs...))
    end
    !linenumbers ? striplines(ex) : ex
end

function ODEProblemExpr(sys::AbstractODESystem, args...; kwargs...)
    ODEProblemExpr{true}(sys, args...; kwargs...)
end

"""
```julia
function DAEProblemExpr{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, parallel=SerialForm(),
                                    skipzeros=true, fillzeros=true,
                                    simplify=false,
                                    kwargs...) where iip
```

Generates a Julia expression for constructing an ODEProblem from an
ODESystem and allows for automatically symbolically calculating
numerical enhancements.
"""
struct DAEProblemExpr{iip} end

function DAEProblemExpr{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
                             parammap = DiffEqBase.NullParameters(); check_length = true,
                             kwargs...) where {iip}
    f, du0, u0, p = process_DEProblem(DAEFunctionExpr{iip}, sys, u0map, parammap;
                                      implicit_dae = true, du0map = du0map, check_length,
                                      kwargs...)
    linenumbers = get(kwargs, :linenumbers, true)
    diffvars = collect_differential_variables(sys)
    sts = states(sys)
    differential_vars = map(Base.Fix2(in, diffvars), sts)

    ex = quote
        f = $f
        u0 = $u0
        du0 = $du0
        tspan = $tspan
        p = $p
        differential_vars = $differential_vars
        DAEProblem{$iip}(f, du0, u0, tspan, p; differential_vars = differential_vars,
                         $(kwargs...))
    end
    !linenumbers ? striplines(ex) : ex
end

function DAEProblemExpr(sys::AbstractODESystem, args...; kwargs...)
    DAEProblemExpr{true}(sys, args...; kwargs...)
end

### Enables Steady State Problems ###
function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem, args...; kwargs...)
    SteadyStateProblem{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem,u0map,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```
Generates an SteadyStateProblem from an ODESystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.SteadyStateProblem{iip}(sys::AbstractODESystem, u0map,
                                            parammap = DiffEqBase.NullParameters();
                                            check_length = true, kwargs...) where {iip}
    f, u0, p = process_DEProblem(ODEFunction{iip}, sys, u0map, parammap;
                                 steady_state = true,
                                 check_length, kwargs...)
    SteadyStateProblem{iip}(f, u0, p; kwargs...)
end

"""
```julia
function DiffEqBase.SteadyStateProblemExpr(sys::AbstractODESystem,u0map,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    skipzeros=true, fillzeros=true,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```
Generates a Julia expression for building a SteadyStateProblem from
an ODESystem and allows for automatically symbolically calculating
numerical enhancements.
"""
struct SteadyStateProblemExpr{iip} end

function SteadyStateProblemExpr{iip}(sys::AbstractODESystem, u0map,
                                     parammap = DiffEqBase.NullParameters();
                                     check_length = true,
                                     kwargs...) where {iip}
    f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap;
                                 steady_state = true,
                                 check_length, kwargs...)
    linenumbers = get(kwargs, :linenumbers, true)
    ex = quote
        f = $f
        u0 = $u0
        p = $p
        SteadyStateProblem(f, u0, p; $(kwargs...))
    end
    !linenumbers ? striplines(ex) : ex
end

function SteadyStateProblemExpr(sys::AbstractODESystem, args...; kwargs...)
    SteadyStateProblemExpr{true}(sys, args...; kwargs...)
end

function _match_eqs(eqs1, eqs2)
    eqpairs = Pair[]
    for (i, eq) in enumerate(eqs1)
        for (j, eq2) in enumerate(eqs2)
            if isequal(eq, eq2)
                push!(eqpairs, i => j)
                break
            end
        end
    end
    eqpairs
end

function isisomorphic(sys1::AbstractODESystem, sys2::AbstractODESystem)
    sys1 = flatten(sys1)
    sys2 = flatten(sys2)

    iv2 = only(independent_variables(sys2))
    sys1 = convert_system(ODESystem, sys1, iv2)
    s1, s2 = states(sys1), states(sys2)
    p1, p2 = parameters(sys1), parameters(sys2)

    (length(s1) != length(s2)) || (length(p1) != length(p2)) && return false

    eqs1 = equations(sys1)
    eqs2 = equations(sys2)

    pps = permutations(p2)
    psts = permutations(s2)
    orig = [p1; s1]
    perms = ([x; y] for x in pps for y in psts)

    for perm in perms
        rules = Dict(orig .=> perm)
        neweqs1 = substitute(eqs1, rules)
        eqpairs = _match_eqs(neweqs1, eqs2)
        if length(eqpairs) == length(eqs1)
            return true
        end
    end
    return false
end


text:
include("/home/kloewer/Oceananigans.jl/src/Oceananigans.jl")
using .Oceananigans

# We'll set up a 2D model with an xz-slice so there's only 1 grid point in y.
Nx, Ny, Nz = 128, 1, 128    # Number of grid points in each dimension.
Lx, Ly, Lz = 1000, 1, 1000  # Domain size (meters).
Nt, Δt = 50000, 10           # Number of time steps, time step size (seconds).

# Set up the model and use an artificially high viscosity ν and diffusivity κ.
model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), Nt=Nt, Δt=Δt, ν=2e-2, κ=2e-2)

# Get location of the cell centers in x and z and reshape them to easily
# broadcast over them when calculating hot_bubble_perturbation.
xC, zC = model.grid.xC, model.grid.zC
xC, zC = reshape(xC, (Nx, 1, 1)), reshape(zC, (1, 1, Nz))

# Set a temperature perturbation with a Gaussian profile located at the center
# of the vertical slice. It roughly corresponds to a background temperature of
# T = 282.99 K and a bubble temperature of T = 283.01 K.
hot_bubble_perturbation = @. 0.01 * exp(-100 * ((xC - Lx/2)^2 + (zC + Lz/2)^2) / (Lx^2 + Lz^2))
model.tracers.T.data .= 282.99 .+ 2 .* reshape(hot_bubble_perturbation, (Nx, Ny, Nz))

# Add a NetCDF output writer that saves NetCDF files to the current directory
nc_writer = NetCDFOutputWriter(dir="data/", frequency=20)
push!(model.output_writers, nc_writer)

time_step!(model)

ext: jl


text:
#= Architecture

Core I/O
  Vital inputs and outputs for various layers
  These include reward signals, sleep/wake signals, system status signals, etc.
Main I/O
  The main inputs and outputs to the system
Gating
  Gates the connections between Main I/O and its target layers
  Gated by the Reinforcement layers
  Bistable activating, and maintains open state when communicating with target layers
  Gain-controls its activity based on activity of relays and target layers
Reinforcement
  Controls the activity of the Gating layers
  Determines its gating policy based on external and internal rewards signals, and whole-system state
Working Memory
  Stores small amounts of information in a stable manner on a short-term basis
  Storage and retrieval policy determined by associated Gating layers
Long Term Memory
  Stores high-level details and relations between them on a long-term basis
  Storage and retrieval policy determined by associated Gating layers
  High sparsity layer allows large information densities
  During storage, spatially- and temporally-close activations become bound together
  During retrieval, previously-bound fragments activate each other
  Retrieval is additionally assisted by the source layers providing reciprocal feedback based
    on newly-retrieved fragments, which allows "walking" along memories in a semi-guided manner
# TODO: Attention/Alerting?
# TODO: General Novelty?
# TODO: CPGs
# TODO: Model/Simulation?
# TODO: Error/Conflict Detection?
# TODO: Control/Monitoring

=#

### Methods ###

function BrainModel(config::Dict)
  _config = deepcopy(config)

  numCol = get!(_config, "numCol", 4)
  numTR = get!(_config, "numTR", numCol)
  szTR = get!(_config, "szTR", 64)
  sztRN = get!(_config, "sztRN", 32)
  szL23e = get!(_config, "szL23e", 32)
  szL4e = get!(_config, "szL4e", 32)
  szL5e = get!(_config, "szL5e", 32)
  szL6e = get!(_config, "szL6e", 32)

  numLanes = get!(_config, "numLanes", 4)
  szUS = get!(_config, "szUS", 8)
  szCS = get!(_config, "szCS", 8)
  szDA = get!(_config, "szDA", 8)
  szPV = get!(_config, "szPV", 8)
  szLV = get!(_config, "szLV", 8)

  numZones = get!(_config, "numZones", 4)
  szGr = get!(_config, "szGr", 32)
  szPk = get!(_config, "szPk", 8)
  szDc = get!(_config, "szDc", 4)
  szIo = get!(_config, "szIo", 4)

  agent = Agent()

  agent["Thalamus"] = ThalamusModel(_config)
  agent["Cortex"] = CortexModel(_config)
  agent["PVLV"] = PVLVModel(_config)
  agent["Cerebellum"] = CerebellumModel(_config)
  agent["Hippocampus"] = HippocampusModel(_config)
  agent["WM"] = WMModel(_config)

  # TODO: Connect inputs to thalamus[TR]
  # TODO: Connect thalamus[TR] to outputs

  # TODO: Connect inputs, cortex[L5] to cerebellum[Gr], cerebellum[Dc]
  # TODO: Connect cerebellum[Dc] to thalamus[TR], outputs

  # Connect thalamus[TR] to cortex[L4e]. cortex[L5e]
  # FIXME: MNE
  for col in keys(agent["Thalamus"]["TR"])
    addedge!(agent, GenericSynapses(szTR, szL4e), agent["Thalamus"]["TR"][col], agent["Cortex"][col]["L4e"]; name="Thalamus_TR_$(col)_Cortex_$(col)_L4e")
    addedge!(agent, GenericSynapses(szTR, szL5e), agent["Thalamus"]["TR"][col], agent["Cortex"][col]["L5e"]; name="Thalamus_TR_$(col)_Cortex_$(col)_L5e")
  end

  # Connect cortex[L6e] to thalamus[TR], thalamus[tRN]
  # FIXME: MNE
  for col in keys(agent["Thalamus"]["TR"])
    addedge!(agent, GenericSynapses(szL6e, szTR), agent["Cortex"][col]["L6e"], agent["Thalamus"]["TR"][col]; name="Cortex_$(col)_L6e_Thalamus_TR_$(col)")
    for trn in keys(agent["Thalamus"]["tRN"])
      addedge!(agent, GenericSynapses(szL6e, sztRN), agent["Cortex"][col]["L6e"], agent["Thalamus"]["tRN"][trn]; name="Cortex_$(col)_L6e_Thalamus_tRN_$(trn)")
    end
  end

  # TODO: Connect lower cortex[L5e] to upper thalamus[TR]
  # TODO: Connect inputs[US] to pvlv[US]
  # TODO: Connect cortex[L5e] to pvlv[CS]
  # TODO: Connect pvlv[???] to thalamus[TR], thalamus[tRN]

  # TODO: Connect cortex[L5e] to hippocampus[ECe]

  return agent
end


text:
## implementation of matchers using lexer methods, in particular RestrictedRegex

"lexer syntax definition for parsing RestrictedRegex definitions"
const RR_SYN = Vector{UInt32}(zero(UInt32),256)
{
    # build RestrictedLexer syntax
    l = UnsafeLexer(IOShared(),RR_SYN)
    addcategory(l,"\\",)
}


text:
module ArrayInterfaceStaticArrays

using Adapt
using ArrayInterface
using LinearAlgebra
using StaticArrays
using Static
import ArrayInterfaceStaticArraysCore

const CanonicalInt = Union{Int,StaticInt}

ArrayInterface.known_first(::Type{<:StaticArrays.SOneTo}) = 1
ArrayInterface.known_last(::Type{StaticArrays.SOneTo{N}}) where {N} = N
ArrayInterface.known_length(::Type{StaticArrays.SOneTo{N}}) where {N} = N
ArrayInterface.known_length(::Type{StaticArrays.Length{L}}) where {L} = L
function ArrayInterface.known_length(::Type{A}) where {A<:StaticArrays.StaticArray}
    ArrayInterface.known_length(StaticArrays.Length(A))
end

ArrayInterface.device(::Type{<:StaticArrays.MArray}) = ArrayInterface.CPUPointer()
ArrayInterface.device(::Type{<:StaticArrays.SArray}) = ArrayInterface.CPUTuple()
ArrayInterface.contiguous_axis(::Type{<:StaticArrays.StaticArray}) = StaticInt{1}()
ArrayInterface.contiguous_batch_size(::Type{<:StaticArrays.StaticArray}) = StaticInt{0}()
function ArrayInterface.stride_rank(::Type{T}) where {N,T<:StaticArray{<:Any,<:Any,N}}
    ntuple(static, StaticInt(N))
end
function ArrayInterface.dense_dims(::Type{<:StaticArray{S,T,N}}) where {S,T,N}
    ArrayInterface._all_dense(Val(N))
end
ArrayInterface.defines_strides(::Type{<:StaticArrays.SArray}) = true
ArrayInterface.defines_strides(::Type{<:StaticArrays.MArray}) = true

@generated function ArrayInterface.axes_types(::Type{<:StaticArrays.StaticArray{S}}) where {S}
    Tuple{[StaticArrays.SOneTo{s} for s in S.parameters]...}
end
@generated function ArrayInterface.size(A::StaticArrays.StaticArray{S}) where {S}
    t = Expr(:tuple)
    Sp = S.parameters
    for n = 1:length(Sp)
        push!(t.args, Expr(:call, Expr(:curly, :StaticInt, Sp[n])))
    end
    return t
end
@generated function ArrayInterface.strides(A::StaticArrays.StaticArray{S}) where {S}
    t = Expr(:tuple, Expr(:call, Expr(:curly, :StaticInt, 1)))
    Sp = S.parameters
    x = 1
    for n = 1:(length(Sp)-1)
        push!(t.args, Expr(:call, Expr(:curly, :StaticInt, (x *= Sp[n]))))
    end
    return t
end
if StaticArrays.SizedArray{Tuple{8,8},Float64,2,2} isa UnionAll
    @inline ArrayInterface.strides(B::StaticArrays.SizedArray{S,T,M,N,A}) where {S,T,M,N,A<:SubArray} = ArrayInterface.strides(B.data)
    ArrayInterface.parent_type(::Type{<:StaticArrays.SizedArray{S,T,M,N,A}}) where {S,T,M,N,A} = A
else
    ArrayInterface.parent_type(::Type{<:StaticArrays.SizedArray{S,T,M,N}}) where {S,T,M,N} = Array{T,N}
end

end # module


text:
# try out some ideas for sequences
#   starting with the direct sequence
using Polylogarithms
using DataFrames, CSV
using PyPlot
using Printf
include("utilities.jl")
using SpecialFunctions

# chose a "hard value"
s = -7.46526-7.56069im
z1 = -7.60722-7.41506im
z = 1/z1
Li = 0.5783+114012.0im
result1 = polylog(s, z1; min_iterations=1000 )
result2 = Polylogarithms.polylog_series_2(s, z1; min_iterations=1000 )
error1 = Li - result1[1]
error2 = Li - result2[1]
rel_error1 = abs(error1) / abs(Li)
rel_error2 = abs(error2) / abs(Li)

k = 1:30
a = z.^k ./ k.^s
k_star = Int(ceil( real(s)/real(log(z)) ))

figure()
semilogy(k,abs.(a), "-")
plot(k, abs.(real.(a)), "-.")
plot(k, abs.(imag.(a)),"--")
plot(k_star, abs.(a[k_star]), "o")

let k=0
    tmp = z
    converged = false
    b = zeros(Complex{Float64}, size(a))
    max_iterations = length(a)
    min_iterations = k_star
    accuracy= 1.0e-12
    while k<max_iterations && ~converged
        k = k+1
        b[k+1] = b[k] + tmp
        tmp *= z * ( k/(k+1.0) )^s
        # println("   total = $total")
        if k > min_iterations && abs(tmp)/abs(b[k+1]) < 0.5*accuracy
            converged = true
        end
    end
    figure()
    plot(abs.(b[1:k]))
    println(b)
end

twoπ = 2*pi
G = (twoπ*im)^s * SpecialFunctions.zeta( 1-s, 0.5 + log(complex(-z1))/(twoπ*im) ) /  SpecialFunctions.gamma(s)


text:
using JukaiNLP
using Merlin

seg = Segmenter()
trainpath = joinpath(dirname(@__FILE__),".data/webtreebank.train")
train(seg, trainpath, 10)
Merlin.save("seg_epoch10.h5", "seg"=>seg)

# load
modelpath = joinpath(dirname(@__FILE__),"seg_epoch10.h5")
seg = Merlin.load(modelpath, "seg")

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
module SubModels

# Perhaps a lot of this could be replaced with
# https://github.com/dpo/NLPModels.jl/blob/master/src/NLPModels.jl

using JuMP
using Docile

import MathProgBase
import ReverseDiffSparse
import Optim
import NLsolve

export JuMPObjective
export get_objective_val, get_objective_deriv
export get_objective_deriv!, get_objective_hess!
export optimize_subobjective

@doc """
A type for holding a particular JuMP sub-objective.
""" ->
type JuMPObjective
	name::ASCIIString
	m::JuMP.Model
	m_eval::JuMP.JuMPNLPEvaluator

	# Which variables have nonzero derivatives.
	vars::Array{Bool, 1}

	colVal::Array{Float64}

	# Stuff for the hessian
	hess_struct::Tuple{Array{Int64,1}, Array{Int64,1}}
	hess_vec::Array{Float64,1}
	numconstr::Int64
	n_params::Int64

	# An in-place gradient
	grad::Array{Float64, 1}
end

JuMPObjective(m::Model) =
	JuMPObjective(m, "model", Bool[true for i in 1:length(m.colVal)])

JuMPObjective(m::Model, name::ASCIIString, vars::Array{Bool, 1}) = begin
	colVal = m.colVal

	#m_const_mat = JuMP.prepConstrMatrix(m);
	m_eval = JuMP.JuMPNLPEvaluator(m);
	MathProgBase.initialize(m_eval, [:ExprGraph, :Grad, :Hess])

	# Structures for the hessian.
	hess_struct = MathProgBase.hesslag_structure(m_eval);
	hess_vec = zeros(length(hess_struct[1]));
	numconstr = (length(m_eval.m.linconstr) +
		           length(m_eval.m.quadconstr) +
		           length(m_eval.m.nlpdata.nlconstr))

	grad = zeros(length(m.colVal))

	n_params = sum(vars)

	JuMPObjective(name, m, m_eval, vars, colVal, hess_struct,
	              hess_vec, numconstr, n_params, grad)
end

JuMPObjective(m::Model, m_eval::JuMP.JuMPNLPEvaluator,
              name::ASCIIString, vars::Array{Bool, 1}) = begin
	colVal = m.colVal

	m_const_mat = JuMP.prepConstrMatrix(m);

	# Structures for the hessian.
	hess_struct = MathProgBase.hesslag_structure(m_eval);
	hess_vec = zeros(length(hess_struct[1]));
	numconstr = (length(m_eval.m.linconstr) +
		           length(m_eval.m.quadconstr) +
		           length(m_eval.m.nlpdata.nlconstr))

	grad = zeros(length(m.colVal))

	n_params = sum(vars)

	JuMPObjective(name, m, m_eval, vars, colVal, hess_struct, hess_vec,
	              numconstr, n_params, grad)
end

@doc """
z_par should be of length sum(jo.vars).
""" ->
function get_objective_val(
		z_par::Array{Float64, 1}, jo::JuMPObjective; verbose=false)

	@assert length(z_par) == jo.n_params
	# jo.m.colVal[jo.vars] = z_par
	jo.colVal[jo.vars] = z_par
	obj_val = jo.m_eval.eval_f_nl(jo.colVal)
	if verbose
		println("$(jo.name) elbo: $obj_val")
	end
	obj_val
end

@doc """
Calculates the derivative with respect to all variables, even excluded ones.
""" ->
function get_full_objective_deriv!(z_par::Array{Float64, 1}, jo::JuMPObjective)
	@assert length(z_par) == jo.n_params
	jo.colVal[jo.vars] = z_par
	MathProgBase.eval_grad_f(jo.m_eval, jo.grad, jo.colVal)
	jo.grad
end

@doc """
Calculates the derivative only with respect to changeable variables.
""" ->
function get_objective_deriv(z_par::Array{Float64, 1}, jo::JuMPObjective)
	get_full_objective_deriv!(z_par, jo)
	jo.grad[jo.vars]
end

@doc """
Calculates the derivative only with respect to changeable variables.
""" ->
function get_objective_deriv!(
		z_par::Array{Float64, 1}, jo::JuMPObjective, grad::Array{Float64, 1})
	grad[:] = get_objective_deriv(z_par, jo)
end

@doc """
Returns the hessian with respect to changeable variables.
""" ->
function get_objective_hess(z_par::Array{Float64, 1}, jo::JuMPObjective)
	@assert length(z_par) == jo.n_params
	jo.colVal[jo.vars] = z_par
	MathProgBase.eval_hesslag(jo.m_eval, jo.hess_vec,
		                        jo.colVal, 1.0, zeros(jo.numconstr))
	this_hess_ld = sparse(jo.hess_struct[1], jo.hess_struct[2],
		                  jo.hess_vec, length(jo.colVal), length(jo.colVal))
	this_hess =
		full(this_hess_ld + this_hess_ld' - sparse(diagm(diag(this_hess_ld))))
	this_hess[jo.vars, jo.vars]
end

@doc """
Returns the hessian only with respect to changeable variables.
""" ->
function get_objective_hess!(z_par::Array{Float64, 1}, jo::JuMPObjective, hess)
	hess[:,:] = get_objective_hess(z_par, jo)
	hess
end

@doc """
From an initial value of z_init for all the variables in m, optimize
the variables associated with the sub-objective jo.
""" ->
function optimize_subobjective(
		z_init::Array{Float64, 1}, jo::JuMPObjective;
    scale=1.0, hess_reg=0.0, show_trace=false, iterations=5, method=:Optim)

	# TODO: clean up the method names.
	@assert(method == :Optim || method == :NLsolve || method == :newton_tr,
	        "Method must be :Optim, :NLsolve, or :newton_tr")
	@assert length(z_init) == length(jo.colVal)
	z_par = z_init[jo.vars]
	jo.colVal = z_init
	function get_local_objective_val(z_par)
		val = scale * get_objective_val(z_par, jo)
		show_trace && println("Value at $(z_par): $(val)")
		val
	end
	function get_local_objective_deriv!(z_par, grad)
		get_objective_deriv!(z_par, jo, grad)
		grad[:] = grad * scale
		#show_trace && println("Gradient at $(z_par): $(grad)")
	end
	function get_local_objective_hess!(z_par, hess)
		if hess_reg < Inf
			get_objective_hess!(z_par, jo, hess)
			hess[:, :] = (hess + eye(length(z_par)) * hess_reg) * scale
		else
			hess[:, :] = scale * eye(length(z_par))
		end
		#show_trace && println("Hessian at $(z_par): $(hess)")
	end

	if method == :Optim
		optim_res = Optim.optimize(get_local_objective_val,
			                         get_local_objective_deriv!,
			                         get_local_objective_hess!,
			                         z_par, method=:newton,
															 show_trace=show_trace,
															 iterations=iterations)
		z_final = z_init
		z_final[jo.vars] = optim_res.minimum
		return z_final, optim_res.f_minimum, optim_res
	elseif method == :newton_tr
		optim_res = Optim.optimize(get_local_objective_val,
															 get_local_objective_deriv!,
															 get_local_objective_hess!,
															 z_par, method=:newton_tr,
															 show_trace=show_trace,
															 iterations=iterations)
		z_final = z_init
		z_final[jo.vars] = optim_res.minimum
		return z_final, optim_res.f_minimum, optim_res
	elseif method == :NLsolve
		nlsolve_res =
			NLsolve.nlsolve(get_local_objective_deriv!, get_local_objective_hess!,
			                z_par, method=:trust_region)
		val = get_local_objective_val(nlsolve_res.zero)
		return nlsolve_res.zero, val, nlsolve_res
	else
		error("Unknown method $method")
	end
end

end # module


text:
"""
    abstract type AbstractBlock

Abstract supertype of all blocks. You should not subtype form this,
but instead from [`Block`](#) or [`WrapperBlock`](#).
"""
abstract type AbstractBlock end

"""
    abstract type Block

A block describes the meaning of a piece of data in the context of a learning task.
For example, for supervised learning tasks, there is an input and a target
and we want to learn to predict targets from inputs. Learning to predict a
cat/dog label from 2D images is a supervised image classification task that can
be represented with the `Block`s `Image{2}()` and `Label(["cat", "dog"])`.

`Block`s are used in virtually every part of the high-level interfaces, from data
processing over model creation to visualization.

## Extending

Consider the following when subtyping `Block`. A block

- Does not hold observation data itself. Instead they are used in conjunction with
  data to annotate it with some meaning.
- If it has any fields, they should be metadata that cannot be
  derived from the data itself and is constant for every sample in
  the dataset. For example `Label` holds all possible classes which
  are constant for the learning problem.

### Interfaces

There are many interfaces that can be implemented for a `Block`. See the docstrings
of each function for more info about how to implement it.

- [`checkblock`](#)`(block, data)`: check whether a piece of data is a valid block
- [`mockblock`](#)`(block)`: randomly generate a piece of data
- [`blocklossfn`](#)`(predblock, yblock)`: loss function for comparing two blocks
- [`blockmodel`](#)`(inblock, outblock[, backbone])`: construct a task-specific model
- [`blockbackbone`](#)`(inblock)`: construct a backbone model that takes in specific data
- [`plotblock!`](#)`(block, data)`: visualize block data

"""
abstract type Block <: AbstractBlock end

"""
    checkblock(block, data)
    checkblock(blocks, datas)

Check whether `data` is compatible with `block`, returning a `Bool`.

## Examples

```julia
checkblock(Image{2}(), rand(RGB, 16, 16)) == true
```

```julia
checkblock(
    (Image{2}(),        Label(["cat", "dog"])),
    (rand(RGB, 16, 16), "cat"                ),
) == true
```

## Extending

An implementation of `checkblock` should be as specific as possible. The
default method returns `false`, so you only need to implement methods for valid types
and return `true`.
"""
checkblock(::Block, data) = false

function checkblock(blocks::Tuple, datas::Tuple)
    @assert length(blocks) == length(datas)
    return all(checkblock(block, data) for (block, data) in zip(blocks, datas))
end

"""
    mockblock(block)
    mockblock(blocks)

Randomly generate an instance of `block`.
"""
mockblock(blocks::Tuple) = map(mockblock, blocks)

"""
    setup(Block, data)

Create an instance of block type `Block` from data container `data`.

## Examples

```julia
setup(Label, ["cat", "dog", "cat"]) == Label(["cat", "dog"])
```

    setup(Encoding, block, data; kwargs...)

Create an encoding using statistics derived from a data container `data`
with observations of block `block`. Used when some arguments of the encoding
are dependent on the dataset. `data` should be the training dataset. Additional
`kwargs` are passed through to the regular constructor of `Encoding`.

## Examples

```julia
(images, labels), blocks = loaddataset("imagenette2-160", (Image, Label))
setup(ImagePreprocessing, Image{2}(), images; buffered = false)
```

```julia
data, block = loaddataset("adult_sample", TableRow)
setup(TabularPreprocessing, block, data)
```
"""
function setup end

# ## Utilities

typify(T::Type) = T
typify(t::Tuple) = Tuple{map(typify, t)...}
typify(block::FastAI.AbstractBlock) = typeof(block)

# ## Block implementations

# Image

"""
    Image{N}() <: Block

`Block` for an N-dimensional mask. `data` is valid for `Image{N}()`
if it is an N-dimensional array with color or number element type.
"""
struct Image{N} <: Block end

checkblock(::Image{N}, ::AbstractArray{T,N}) where {T <: Union{Colorant,Number},N} = true
mockblock(::Image{N}) where N = rand(RGB{N0f8}, ntuple(_ -> 16, N))

setup(::Type{Image}, data) = Image{ndims(getobs(data, 1))}()

"""
    Mask{N, T}(classes) <: Block

Block for an N-dimensional categorical mask. `data` is valid for
`Mask{N, T}(classes)`
if it is an N-dimensional array with every element in `classes`.
"""
struct Mask{N,T} <: Block
    classes::AbstractVector{T}
end
Mask{N}(classes::AbstractVector{T}) where {N,T} = Mask{N,T}(classes)

function checkblock(block::Mask{N,T}, a::AbstractArray{T,N}) where {N,T}
    return all(map(x -> x ∈ block.classes, a))
end

mockblock(mask::Mask{N}) where N = rand(mask.classes, ntuple(_ -> 16, N))

# Keypoints

"""

    Keypoints{N}(sz) <: Block

A block representing an array of size `sz` filled with keypoints of type
`SVector{N}`.
"""
struct Keypoints{N, M} <: Block
    sz::NTuple{M, Int}
end
Keypoints{N}(n::Int) where N = Keypoints{N, 1}((n,))
Keypoints{N}(t::NTuple{M, Int}) where {N, M} = Keypoints{N, M}(t)

function checkblock(
        block::Keypoints{N,M},
        a::AbstractArray{<:Union{SVector{N,T},Nothing},M}) where {M,N,T}
    return true
end

mockblock(block::Keypoints{N}) where N = rand(SVector{N, Float32}, block.sz)

# TableRow

"""
    TableRow{M, N}(catcols, contcols, categorydict) <: Block

`Block` for table rows with M categorical and N continuous columns. `data`
is valid if it satisfies the `AbstractRow` interface in Tables.jl, values
present in indices for categorical and continuous columns are consistent,
and `data` is indexable by the elements of `catcols` and `contcols`.
"""
struct TableRow{M, N, T} <: Block
    catcols::NTuple{M}
    contcols::NTuple{N}
    categorydict::T
end

function TableRow(catcols, contcols, categorydict)
    TableRow{length(catcols), length(contcols)}(catcols, contcols, categorydict)
end

function checkblock(block::TableRow, x)
    columns = Tables.columnnames(x)
    (all(col -> col ∈ columns, (block.catcols..., block.contcols...)) &&
    all(col -> haskey(block.categorydict, col) &&
        (ismissing(x[col]) || x[col] ∈ block.categorydict[col]), block.catcols) &&
    all(col -> ismissing(x[col]) || x[col] isa Number, block.contcols))
end

function mockblock(block::TableRow)
    cols = (block.catcols..., block.contcols...)
    vals = map(cols) do col
        col in block.catcols ?
            rand(block.categorydict[col]) : rand()
    end
    return NamedTuple(zip(cols, vals))
end

"""
    setup(TableRow, data[; catcols, contcols])

Create a `TableRow` block from data container `data::TableDataset`. If the
categorical and continuous columns are not specified manually, try to
guess them from the dataset's column types.
"""
function setup(::Type{TableRow}, data; catcols=nothing, contcols=nothing)
    catcols_, contcols_ = getcoltypes(data)
    catcols = isnothing(catcols) ? catcols_ : catcols
    contcols = isnothing(contcols) ? contcols_ : contcols

    return TableRow(
        catcols,
        contcols,
        gettransformdict(data, DataAugmentation.Categorify, catcols))
end

function Base.show(io::IO, block::TableRow)
    print(io, ShowCase(block, (:catcols, :contcols), show_params=false, new_lines=true))
end

# Continous

"""
    Continuous(size) <: Block

`Block` for collections of numbers. `data` is valid if it's
length is `size` and contains `Number`s.
"""

struct Continuous <: Block
    size::Int
end

function checkblock(block::Continuous, x)
    block.size == length(x) && eltype(x) <: Number
end

mockblock(block::Continuous) = rand(block.size)


text:
# =============================================================================
# rna.jl -> Funciones útiles para crear y entrenar RNAs:
#   - buildClassANN
#   - trainClassANN
# =============================================================================

using Flux: Chain, Dense, σ, softmax, Losses, params, ADAM, train!;

#
# Función que crea una rna según la topología que se le indica.
#
# @arguments
#   numInputs: número de capas de entrada.
#   topology: topología de las capas ocultas.
#   numOutputs: número de capas de salida.
#
# @return: rna
#
function buildClassANN(numInputs::Int64, topology::Array{Int64,1}, numOutputs::Int64)
    ann=Chain();
    numInputsLayer = numInputs;
    for numOutputLayers = topology
        ann = Chain(ann..., Dense(numInputsLayer, numOutputLayers, σ));
        numInputsLayer = numOutputLayers;
    end;
    #Dependiendo del nº de salidas añadimos las capas de salida adecuadas
    if (numOutputs == 1)
        ann = Chain(ann..., Dense(numInputsLayer, 1, σ));
    else
        ann = Chain(ann..., Dense(numInputsLayer, numOutputs, identity));
        ann = Chain(ann..., softmax);
    end;
    return ann;
end;

#
# Funcion para entrenar RR.NN.AA. con conjuntos de entrenamiento, validacion y
#   test. Calcular errores en el conjunto de validacion, y para el entrenamiento
#   si es necesario.
#
# @arguments
#   topology: topología de las capas ocultas.
#   trainingInputs: las entradas que representan el conjunto de entrenamiento.
#   trainingTargets: las salidas deseadas del conjunto de entrenamiento.
#   validationInputs: inputs para realizar la validación del dataset.
#   validationTargets: salidas deseadas para realizar la validación del dataset.
#   testInputs: inputs para el conjunto de test.
#   testTargets: salidas deseadas para el conjunto de test.
#   maxEpochs: número máximo de iteraciones.
#   minLoss: mínimo error de entrenamiento.
#   learningRate: tasa de aprendizaje.
#   maxEpochsVal: número máximo de iteraciones sin mejorar en el conjunto de validacion.
#   showText: saca por pantalla el texto si es cierto
#
# @return: rna entrenada. ¿¿home devolverá mais cousiñas non??? mucha mariconada esta inspiración en los javadocs
#
function trainClassANN(topology::Array{Int64,1}, trainingInputs::Array{Float64,2}, trainingTargets::Array{Bool,2},
    validationInputs::Array{Float64,2}, validationTargets::Array{Bool,2}, testInputs::Array{Float64,2}, testTargets::Array{Bool,2};
    maxEpochs::Int64=1000, minLoss::Float64=0.0, learningRate::Float64=0.1, maxEpochsVal::Int64=6, showText::Bool=false)

    # Se supone que tenemos cada patron en cada fila
    # Comprobamos que el numero de filas (numero de patrones) coincide tanto en entrenamiento como en validation como test
    @assert(size(trainingInputs,1)==size(trainingTargets,1));
    @assert(size(validationInputs,1)==size(validationTargets,1));
    @assert(size(testInputs,1)==size(testTargets,1));
    # Comprobamos que el numero de columnas coincide en los grupos de entrenamiento, validacion y test
    @assert(size(trainingInputs,2)==size(validationInputs,2)==size(testInputs,2));
    @assert(size(trainingTargets,2)==size(validationTargets,2)==size(testTargets,2));
    # Creamos la RNA
    ann = buildClassANN(size(trainingInputs,2), topology, size(trainingTargets,2));
    # Definimos la funcion de loss
    loss(x,y) = (size(y,1) == 1) ? Losses.binarycrossentropy(ann(x),y) : Losses.crossentropy(ann(x),y);
    # Creamos los vectores con los valores de loss y de precision en cada ciclo
    trainingLosses = Float64[];
    trainingAccuracies = Float64[];
    validationLosses = Float64[];
    validationAccuracies = Float64[];
    testLosses = Float64[];
    testAccuracies = Float64[];

    # Empezamos en el ciclo 0
    numEpoch = 0;

    # Una funcion util para calcular los resultados y mostrarlos por pantalla
    function calculateMetrics()
        # Calculamos el loss en entrenamiento y test. Para ello hay que pasar las matrices traspuestas (cada patron en una columna)
        trainingLoss = loss(trainingInputs', trainingTargets');
        validationLoss = loss(validationInputs', validationTargets');
        testLoss = loss(testInputs', testTargets');
        # Calculamos la salida de la RNA en entrenamiento y test. Para ello hay que pasar la matriz de entradas traspuesta (cada patron en una columna). La matriz de salidas tiene un patron en cada columna
        trainingOutputs = ann(trainingInputs');
        validationOutputs = ann(validationInputs');
        testOutputs = ann(testInputs');
        # Para calcular la precision, ponemos 2 opciones aqui equivalentes:
        #  Pasar las matrices con los datos en las columnas. La matriz de salidas ya tiene un patron en cada columna
        trainingAcc = accuracy(trainingOutputs, Array{Bool,2}(trainingTargets'); dataInRows=false);
        validationAcc = accuracy(validationOutputs, Array{Bool,2}(validationTargets'); dataInRows=false);
        testAcc = accuracy(testOutputs, Array{Bool,2}(testTargets'); dataInRows=false);
        #  Pasar las matrices con los datos en las filas. Hay que trasponer la matriz de salidas de la RNA, puesto que cada dato esta en una fila
        trainingAcc = accuracy(Array{Float64,2}(trainingOutputs'), trainingTargets;   dataInRows=true);
        validationAcc = accuracy(Array{Float64,2}(validationOutputs'), validationTargets; dataInRows=true);
        testAcc = accuracy(Array{Float64,2}(testOutputs'), testTargets; dataInRows=true);
        # Mostramos por pantalla el resultado de este ciclo de entrenamiento si nos lo han indicado
        if showText
            #println("Epoch ", numEpoch, ": Training loss: ", trainingLoss, ", accuracy: ", 100*trainingAcc, " % - Validation loss: ", validationLoss, ", accuracy: ", 100*validationAcc, " % - Test loss: ", testLoss, ", accuracy: ", 100*testAcc, " %");
            println("Epoch ", numEpoch);
            println("\t Training loss: ", trainingLoss, ", accuracy: ", 100*trainingAcc);
            println("\t Validation loss: ", validationLoss, ", accuracy: ", 100*validationAcc);
            println("\t Test loss: ", testLoss, ", accuracy: ", 100*testAcc);
        end;
        return (trainingLoss, trainingAcc, validationLoss, validationAcc, testLoss, testAcc)
    end;

    # Calculamos las metricas para el ciclo 0 (sin entrenar nada)
    (trainingLoss, trainingAccuracy, validationLoss, validationAccuracy, testLoss, testAccuracy) = calculateMetrics();
    #  y almacenamos los valores de loss y precision en este ciclo
    push!(trainingLosses, trainingLoss);
    push!(trainingAccuracies, trainingAccuracy);
    push!(validationLosses, validationLoss);
    push!(validationAccuracies, validationAccuracy);
    push!(testLosses, testLoss);
    push!(testAccuracies, testAccuracy);

    # Numero de ciclos sin mejorar el error de validacion y el mejor error de validation encontrado hasta el momento
    numEpochsValidation = 0; bestValidationLoss = validationLoss;
    # Cual es la mejor ann que se ha conseguido
    bestANN = deepcopy(ann);

    # Entrenamos hasta que se cumpla una condicion de parada
    while (numEpoch<maxEpochs) && (trainingLoss>minLoss) && (numEpochsValidation<maxEpochsVal)
        # Entrenamos 1 ciclo. Para ello hay que pasar las matrices traspuestas (cada patron en una columna)
        train!(loss, params(ann), [(trainingInputs', trainingTargets')], ADAM(learningRate));
        # Aumentamos el numero de ciclo en 1
        numEpoch += 1;
        # Calculamos las metricas en este ciclo
        (trainingLoss, trainingAccuracy, validationLoss, validationAccuracy, testLoss, testAccuracy) = calculateMetrics();
        #  y almacenamos los valores de loss y precision en este ciclo
        push!(trainingLosses, trainingLoss);
        push!(trainingAccuracies, trainingAccuracy);
        push!(validationLosses, validationLoss);
        push!(validationAccuracies, validationAccuracy);
        push!(testLosses, testLoss);
        push!(testAccuracies,testAccuracy);
        # Aplicamos la parada temprana
        if (validationLoss<bestValidationLoss)
            bestValidationLoss = validationLoss;
            numEpochsValidation = 0;
            bestANN = deepcopy(ann);
        else
            numEpochsValidation += 1;
        end;
    end;
    return (bestANN, trainingLosses, validationLosses, testLosses, trainingAccuracies, validationAccuracies, testAccuracies);
end;

#sin validacion
function trainClassANN(topology::Array{Int64,1},
    trainingInputs::Array{Float64,2}, trainingTargets::Array{Bool,2},
    testInputs::Array{Float64,2},     testTargets::Array{Bool,2};
    maxEpochs::Int64=1000, minLoss::Float64=0.0, learningRate::Float64=0.1, showText::Bool=false)

    # Se supone que tenemos cada patron en cada fila
    # Comprobamos que el numero de filas (numero de patrones) coincide tanto en entrenamiento como en test
    @assert(size(trainingInputs,1)==size(trainingTargets,1));
    @assert(size(testInputs,1)==size(testTargets,1));
    # Comprobamos que el numero de columnas coincide en los grupos de entrenamiento y test
    @assert(size(trainingInputs,2)==size(testInputs,2));
    @assert(size(trainingTargets,2)==size(testTargets,2));
    # Creamos la RNA
    ann = buildClassANN(size(trainingInputs,2), topology, size(trainingTargets,2));
    # Definimos la funcion de loss
    loss(x,y) = (size(y,1) == 1) ? Losses.binarycrossentropy(ann(x),y) : Losses.crossentropy(ann(x),y);
    # Creamos los vectores con los valores de loss y de precision en cada ciclo
    trainingLosses = Float64[];
    trainingAccuracies = Float64[];
    testLosses = Float64[];
    testAccuracies = Float64[];
    # Empezamos en el ciclo 0
    numEpoch = 0;

    # Una funcion util para calcular los resultados y mostrarlos por pantalla
    function calculateMetrics()
        # Calculamos el loss en entrenamiento y test. Para ello hay que pasar las matrices traspuestas (cada patron en una columna)
        trainingLoss = loss(trainingInputs', trainingTargets');
        testLoss     = loss(testInputs',     testTargets');
        # Calculamos la salida de la RNA en entrenamiento y test. Para ello hay que pasar la matriz de entradas traspuesta (cada patron en una columna). La matriz de salidas tiene un patron en cada columna
        trainingOutputs = ann(trainingInputs');
        testOutputs     = ann(testInputs');
        # Para calcular la precision, ponemos 2 opciones aqui equivalentes:
        #  Pasar las matrices con los datos en las columnas. La matriz de salidas ya tiene un patron en cada columna
        trainingAcc = accuracy(trainingOutputs, Array{Bool,2}(trainingTargets'); dataInRows=false);
        testAcc     = accuracy(testOutputs,     Array{Bool,2}(testTargets');     dataInRows=false);
        #  Pasar las matrices con los datos en las filas. Hay que trasponer la matriz de salidas de la RNA, puesto que cada dato esta en una fila
        trainingAcc = accuracy(Array{Float64,2}(trainingOutputs'), trainingTargets; dataInRows=true);
        testAcc     = accuracy(Array{Float64,2}(testOutputs'),     testTargets;     dataInRows=true);
        # Mostramos por pantalla el resultado de este ciclo de entrenamiento si nos lo han indicado
        if showText
            println("Epoch ", numEpoch, ": Training loss: ", trainingLoss, ", accuracy: ", 100*trainingAcc, " % - Test loss: ", testLoss, ", accuracy: ", 100*testAcc, " %");
        end;
        return (trainingLoss, trainingAcc, testLoss, testAcc)
    end;

    # Calculamos las metricas para el ciclo 0 (sin entrenar nada)
    (trainingLoss, trainingAccuracy, testLoss, testAccuracy) = calculateMetrics();
    #  y almacenamos los valores de loss y precision en este ciclo
    push!(trainingLosses,      trainingLoss);
    push!(testLosses,          testLoss);
    push!(trainingAccuracies,  trainingAccuracy);
    push!(testAccuracies,      testAccuracy);

    # Entrenamos hasta que se cumpla una condicion de parada
    while (numEpoch<maxEpochs) && (trainingLoss>minLoss)
        # Entrenamos 1 ciclo. Para ello hay que pasar las matrices traspuestas (cada patron en una columna)
        Flux.train!(loss, params(ann), [(trainingInputs', trainingTargets')], ADAM(learningRate));
        # Aumentamos el numero de ciclo en 1
        numEpoch += 1;
        # Calculamos las metricas en este ciclo
        (trainingLoss, trainingAccuracy, testLoss, testAccuracy) = calculateMetrics();
        #  y almacenamos los valores de loss y precision en este ciclo
        push!(trainingLosses,     trainingLoss);
        push!(trainingAccuracies, trainingAccuracy);
        push!(testLosses,         testLoss);
        push!(testAccuracies,     testAccuracy);
    end;
    return (ann, trainingLosses, testLosses, trainingAccuracies, testAccuracies);
end;


text:
using NiLang
using Test, LinearAlgebra

function naive_umm!(x, params)
    N = size(x, 1)
    k = 0
    for j=1:N
        for i=N-1:-1:j
            k += 1
            a, b = rot(x[i],x[i+1],params[k])
            x[i], x[i+1] = a, b
        end
    end
end

function inv_naive_umm!(x, params)
    N = size(x, 1)
    k = N*(N-1) ÷ 2
    for j=N:-1:1
        for i=j:N-1
            a, b = rot(x[i],x[i+1],-params[k])
            x[i], x[i+1] = a, b
            k -= 1
        end
    end
end

@testset "naive unitary" begin
    x = randn(200)
    params = randn(100*199).*2π

    x0 = copy(x)
    params0 = copy(params)
    naive_umm!(x, params)
    inv_naive_umm!(x, params)
    @test params ≈ params0
    @test x ≈ x0
end

@testset "unitary" begin
    x = randn(20)
    params = randn(10*19) * 2π

    x0 = copy(x)
    params0 = copy(params)
    Nx = length(x)
    @instr umm!(x, params, Nx, Nx)
    x1 = copy(x0)
    params1 = copy(params0)
    naive_umm!(x1, params1)
    @test params ≈ params1
    @test x ≈ x1
    @instr (~umm!)(x, params, Nx, Nx)
    @test params ≈ params0
    @test x ≈ x0
    @test check_inv(umm!, (x, params, Nx, Nx))
end

@testset "imapfoldl" begin
    x = randn(20)
    out = 0.0
    @test_throws InvertibilityError MapFoldl(exp, +)
    @instr MapFoldl(exp, ⊕(identity))(out, x)
    @test out ≈ sum(exp.(x))
    @test check_inv(MapFoldl(exp, ⊕(identity)), (out, x))
end


text:
# ---
# title: 109. Convert Sorted List to Binary Search Tree
# id: problem109
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List, Depth-first Search
# link: <https://leetcode.com/problems/convert-sorted-list-to-binary-search-tree/description/>
# hidden: true
# ---
#
# Given the `head` of a singly linked list where elements are **sorted in
# ascending order** , convert it to a height balanced BST.
#
# For this problem, a height-balanced binary tree is defined as a binary tree in
# which the depth of the two subtrees of _every_ node never differ by more than
# 1.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/17/linked.jpg)
#
#
#
#     Input: head = [-10,-3,0,5,9]
#     Output: [0,-3,9,-10,null,5]
#     Explanation: One possible answer is [0,-3,9,-10,null,5], which represents the shown height balanced BST.
#
#
# **Example 2:**
#
#
#
#     Input: head = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: head = [0]
#     Output: [0]
#
#
# **Example 4:**
#
#
#
#     Input: head = [1,3]
#     Output: [3,1]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in `head` is in the range `[0, 2 * 104]`.
#   * `-10^5 <= Node.val <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export KeepActor, keep, getvalues

"""
    KeepActor{D}() where D

Keep actor provides a storage actor. It saves all incoming successful `next` events in a `values` array, throws an ErrorException on `error!` event
and does nothing on completion event.

# Examples
```jldoctest
using Rocket

source = from(1:5)
actor  = keep(Int)

subscribe!(source, actor)
show(getvalues(actor))

# output
[1, 2, 3, 4, 5]
```

See also: [`Actor`](@ref), [`keep`](@ref)
"""
struct KeepActor{T} <: Actor{T}
    values :: Vector{T}

    KeepActor{T}() where T = new(Vector{T}())
end

getvalues(actor::KeepActor) = actor.values

on_next!(actor::KeepActor{T}, data::T) where T = push!(actor.values, data)
on_error!(actor::KeepActor, err)               = error(err)
on_complete!(actor::KeepActor)                 = begin end

"""
    keep(::Type{T}) where T

# Arguments
- `::Type{T}`: Type of keep data

Creation operator for the `KeepActor` actor.

# Examples

```jldoctest
using Rocket

actor = keep(Int)
actor isa KeepActor{Int}

# output
true
```

See also: [`KeepActor`](@ref), [`AbstractActor`](@ref)
"""
keep(::Type{T}) where T = KeepActor{T}()


text:
mutable struct CLCompilationDatabase
    ptr::CXCompilationDatabase
end

"""
    CLCompileCommand

Represents a command to compile the specific file.
Call [`getDirectory`](@ref), [`getArguments`](@ref) and [`getArguments`](@ref) to get its attributes.
"""
mutable struct CLCompileCommand
    # Keep a reference to its root
    root::Any
    ptr::CXCompileCommand
end

mutable struct CLCompileCommands <: AbstractVector{CLCompileCommands}
    ptr::CXCompileCommands
end

"""
    CLCompilationDatabase(path::AbstractString)

Creates a compilation database from the database found in directory `path`.
"""
function CLCompilationDatabase(path::AbstractString)
    error = Ref{CXCompilationDatabase_Error}()
    ptr = clang_CompilationDatabase_fromDirectory(path, error)
    if error[] == CXCompilationDatabase_NoError
        db = CLCompilationDatabase(ptr)
        finalizer(db) do x
            clang_CompilationDatabase_dispose(x.ptr)
        end
    else
        nothing
    end
end

"""
    getCompileCommands(db::CLCompilationDatabase)

Get all compile commands from a compilation database.
"""
function getCompileCommands(db::CLCompilationDatabase)
    ptr = clang_CompilationDatabase_getAllCompileCommands(db.ptr)
    cmds = CLCompileCommands(ptr)
    finalizer(cmds) do x
        clang_CompileCommands_dispose(x.ptr)
    end
    cmds
end

"""
    getCompileCommands(db::CLCompilationDatabase, file::AbstractString)

Get all compile commands for a file given its complete path.
"""
function getCompileCommands(db::CLCompilationDatabase, file::AbstractString)
    ptr = clang_CompilationDatabase_getCompileCommands(db.ptr, file)
    cmds = CLCompileCommands(ptr)
    finalizer(cmds) do x
        clang_CompileCommands_dispose(x.ptr)
    end
    cmds
end

Base.size(cmds::CLCompileCommands) = (clang_CompileCommands_getSize(cmds.ptr),)

function Base.getindex(cmds::CLCompileCommands, i)
    ptr = clang_CompileCommands_getCommand(cmds.ptr, i - 1)
    CLCompileCommand(cmds, ptr)
end

getDirectory(cmd::CLCompileCommand) = _cxstring_to_string(clang_CompileCommand_getDirectory(cmd.ptr))

getFileName(cmd::CLCompileCommand) = _cxstring_to_string(clang_CompileCommand_getFilename(cmd.ptr))

function getArguments(cmd::CLCompileCommand)
    num = clang_CompileCommand_getNumArgs(cmd.ptr)
    args = Vector{String}(undef, num)
    for i in 1:num
        args[i] = _cxstring_to_string(clang_CompileCommand_getArg(cmd.ptr, i - 1))
    end
    args
end

# Mapped source support is dropped by clang


text:
module PiecewiseIncreasingRanges
export PiecewiseIncreasingRange, NoNearestSampleError, findnearest, resample

using Compat, Base.Order

function combine_ranges{R<:Range}(ranges::Vector{R}, firstrg::R, firstrgidx::Int)
    newranges = R[]
    offsets = Int[1]

    step(firstrg) < 0 && throw(ArgumentError("ranges must be strictly monotonically increasing"))
    currg = firstrg

    for i = firstrgidx:length(ranges)
        newrg = ranges[i]
        isempty(newrg) && continue

        if step(newrg) == step(currg) && first(newrg) == last(currg)+step(currg)
            # Can extend current range
            currg = first(currg):step(currg):last(newrg)
        else
            if first(newrg) <= last(currg) || step(newrg) < 0
                throw(ArgumentError("ranges must be strictly monotonically increasing"))
            end

            # Need to make a new range
            push!(newranges, currg)
            push!(offsets, offsets[end]+length(currg))
            currg = newrg
        end
    end

    push!(newranges, currg)
    (newranges, offsets)
end

function combine_ranges{R<:FloatRange}(ranges::Vector{R}, firstrg::R, firstrgidx::Int)
    newranges = R[]
    offsets = Int[1]

    if signbit(firstrg.step) != signbit(firstrg.divisor)
        throw(ArgumentError("ranges must be strictly monotonically increasing"))
    end
    currg = firstrg

    for i = firstrgidx:length(ranges)
        newrg = ranges[i]
        isempty(newrg) && continue

        if step(newrg) == step(currg) && newrg.start*currg.divisor == (currg.start+currg.step*currg.len)*newrg.divisor
            currg = FloatRange(currg.start, currg.step, currg.len + newrg.len, currg.divisor)
        else
            if first(newrg) <= last(currg) || signbit(newrg.step) != signbit(newrg.divisor)
                throw(ArgumentError("ranges must be strictly monotonically increasing"))
            end

            # Need to make a new range
            push!(newranges, currg)
            push!(offsets, offsets[end]+length(currg))
            currg = newrg
        end
    end

    push!(newranges, currg)
    (newranges, offsets)
end

immutable PiecewiseIncreasingRange{T,R<:Range,S} <: AbstractVector{T}
    ranges::Vector{R}
    offsets::Vector{Int}
    divisor::S

    function PiecewiseIncreasingRange(ranges::Vector{R}, divisor)
        isempty(ranges) && return new(ranges, Int[])

        # Find first non-empty range
        firstrg = ranges[1]
        j = 0
        for j = 1:length(ranges)
            firstrg = ranges[j]
            !isempty(firstrg) && break
        end
        isempty(firstrg) && return new(R[], Int[], divisor)

        newranges, offsets = combine_ranges(ranges, firstrg, j+1)
        new(newranges, offsets, divisor)
    end
end
PiecewiseIncreasingRange{R<:Range}(ranges::Vector{R}, divisor) = PiecewiseIncreasingRange{typeof(inv(one(eltype(R)))),R,typeof(divisor)}(ranges, divisor)
PiecewiseIncreasingRange{R<:Range}(ranges::Vector{R}) = PiecewiseIncreasingRange{eltype(R),R,@compat(Void)}(ranges, nothing)

# Avoid applying the divisor if it is one, to get types right
divide_divisor{T,R}(r::PiecewiseIncreasingRange{T,R,@compat(Void)}, x) = x
multiply_divisor{T,R}(r::PiecewiseIncreasingRange{T,R,@compat(Void)}, x) = x
divide_divisor{T,R,S}(r::PiecewiseIncreasingRange{T,R,S}, x) = x/r.divisor
multiply_divisor{T,R,S}(r::PiecewiseIncreasingRange{T,R,S}, x) = x*r.divisor

function Base.size(r::PiecewiseIncreasingRange)
    isempty(r.ranges) && return 0
    return (r.offsets[end]+length(r.ranges[end])-1,)
end

function Base.getindex(r::PiecewiseIncreasingRange, i::Integer)
    rgidx = searchsortedlast(r.offsets, i, Forward)
    divide_divisor(r, r.ranges[rgidx][i-r.offsets[rgidx]+1])
end

function Base.getindex{T,R,S}(r::PiecewiseIncreasingRange{T,R,S}, x::Range{Int})
    isempty(x) && return PiecewiseIncreasingRange{T,R,S}(R[], Int[], r.divisor)
    (first(x) >= 1 && last(x) <= length(r)) || throw(BoundsError())

    # Not type-stable!
    step(x) < 0 && return invoke(getindex, (typeof(r), AbstractVector{Int}), r, x)

    firstrgidx = searchsortedlast(r.offsets, first(x), Forward)
    lastrgidx = searchsortedlast(r.offsets, last(x), Forward)
    newrgs = Array(R, lastrgidx-firstrgidx+1)
    for irange = firstrgidx:lastrgidx
        elmax = min(last(x)-r.offsets[irange]+1, length(r.ranges[irange]))
        newrg = newrgs[irange-firstrgidx+1] = r.ranges[irange][first(x)-r.offsets[irange]+1:elmax]
        x = x[length(newrg)+1:end]
    end
    PiecewiseIncreasingRange(newrgs, r.divisor)
end

function resample{T,R,S}(r::PiecewiseIncreasingRange{T,R,S}, ratio::Rational{Int})
    excess = zero(eltype(R))
    curpt = 0
    newrgs = R[]
    for i = 1:length(r.ranges)
        endpt = last(r.ranges[i])*num(ratio)
        newrg = first(r.ranges[i])*num(ratio)+excess*step(r.ranges[i]):step(r.ranges[i])*den(ratio):endpt
        curpt += length(newrg)
        push!(newrgs, newrg)
        if i != length(r.ranges)
            # Need to linearly interpolate between endpoint and next range
            # nextind is the next index in r (= 1 + curpt/ratio) times num(ratio)
            nextind = num(ratio) + curpt*den(ratio)
            if nextind < r.offsets[i+1]*num(ratio)
                # Need to interpolate between this point and the next
                weight = mod(nextind, num(ratio))
                interpt = last(r.ranges[i])*(num(ratio)-weight) + first(r.ranges[i+1])*weight
                interpstep = (first(r.ranges[i+1]) - last(r.ranges[i]))*den(ratio)
                newrg = interpt:interpstep:first(r.ranges[i+1])*num(ratio)
                push!(newrgs, newrg)
                curpt += length(newrg)
                nextind += length(newrg)*den(ratio)
            end
            excess = nextind - r.offsets[i+1]*num(ratio)
        end
    end
    PiecewiseIncreasingRange(newrgs, multiply_divisor(r, num(ratio)))
end

# searchsortedfirst, searchsortedlast
immutable PiecewiseIncreasingRangeFirstOrdering <: Ordering end
Base.Order.lt(o::PiecewiseIncreasingRangeFirstOrdering, a, b) = isless(first(a), first(b))
immutable PiecewiseIncreasingRangeLastOrdering <: Ordering end
Base.Order.lt(o::PiecewiseIncreasingRangeLastOrdering, a, b) = isless(last(a), last(b))

function Base.searchsortedfirst(r::PiecewiseIncreasingRange, x)
    isempty(r.ranges) && return 1
    xd = multiply_divisor(r, x)

    rgidx = searchsortedfirst(r.ranges, xd, PiecewiseIncreasingRangeLastOrdering())
    rgidx > length(r.ranges) && return length(r) + 1
    searchsortedfirst(r.ranges[rgidx], xd, Forward) + r.offsets[rgidx] - 1
end

function Base.searchsortedlast(r::PiecewiseIncreasingRange, x)
    isempty(r.ranges) && return 1
    xd = multiply_divisor(r, x)

    rgidx = searchsortedlast(r.ranges, xd, PiecewiseIncreasingRangeFirstOrdering())
    rgidx == 0 && return 0
    searchsortedlast(r.ranges[rgidx], xd, Forward) + r.offsets[rgidx] - 1
end

immutable NoNearestSampleError <: Exception end

function findnearest(r::PiecewiseIncreasingRange, x, within_half_step::Bool=false)
    isempty(r.ranges) && throw(NoNearestSampleError())
    xd = multiply_divisor(r, x)

    rgidx = searchsortedfirst(r.ranges, xd, PiecewiseIncreasingRangeLastOrdering())
    if rgidx > length(r.ranges)
        rgend = r.ranges[end]
        within_half_step && xd > rgend[end]+step(rgend)/2 && throw(NoNearestSampleError())
        return length(r)
    end

    rg = r.ranges[rgidx]
    idxinrg = searchsortedfirst(rg, xd, Forward)
    idx = idxinrg + r.offsets[rgidx] - 1
    d = rg[idxinrg] - xd

    if idxinrg == 1
        if rgidx == 1
            # First element of all elements
            within_half_step && d > step(rg)/2 && throw(NoNearestSampleError())
            return 1
        end

        # Could be closer to last element of preceding range
        rgprev = r.ranges[rgidx-1]
        dprev = xd - rgprev[end]
        within_half_step && d > step(rg)/2 && dprev > step(rgprev)/2 && throw(NoNearestSampleError())
        return idx - (dprev <= d)
    end

    ifelse(d >= step(rg)/2, idx-1, idx)
end

end # module


text:
using MacroTools

struct BrainFlowError <: Exception
    msg::String
    ec::Int32
end

@enum BrainflowExitCodes begin

    STATUS_OK = 0
    PORT_ALREADY_OPEN_ERROR = 1
    UNABLE_TO_OPEN_PORT_ERROR = 2
    SER_PORT_ERROR = 3
    BOARD_WRITE_ERROR = 4
    INCOMMING_MSG_ERROR = 5
    INITIAL_MSG_ERROR = 6
    BOARD_NOT_READY_ERROR = 7
    STREAM_ALREADY_RUN_ERROR = 8
    INVALID_BUFFER_SIZE_ERROR = 9
    STREAM_THREAD_ERROR = 10
    STREAM_THREAD_IS_NOT_RUNNING = 11
    EMPTY_BUFFER_ERROR = 12
    INVALID_ARGUMENTS_ERROR = 13
    UNSUPPORTED_BOARD_ERROR = 14
    BOARD_NOT_CREATED_ERROR = 15
    ANOTHER_BOARD_IS_CREATED_ERROR = 16
    GENERAL_ERROR = 17
    SYNC_TIMEOUT_ERROR = 18
    JSON_NOT_FOUND_ERROR = 19
    NO_SUCH_DATA_IN_JSON_ERROR = 20
    CLASSIFIER_IS_NOT_PREPARED_ERROR = 21
    ANOTHER_CLASSIFIER_IS_PREPARED_ERROR = 22
    UNSUPPORTED_CLASSIFIER_AND_METRIC_COMBINATION_ERROR = 23
end

"""
    @brainflow_rethrow

Macro that generates the brainflow error handling automatically after a ccall.

# Example
```
@brainflow_rethrow function prepare_session(board_shim::BoardShim)
    ccall((:prepare_session, BOARD_CONTROLLER_INTERFACE), Cint, (Cint, Ptr{UInt8}), board_shim.board_id, board_shim.input_json)
end
```

"""
macro brainflow_rethrow(defun)
    def = MacroTools.splitdef(defun)
    body = def[:body]
    name = string(def[:name])

    # find the ccall
    # assuming there is only 1 ccall
    # assuming you do NOT make an assignment, like x_ = ccall(xs__)
    args = body.args
    n_args = length(args)
    ccall_index = findfirst(x -> MacroTools.@capture(x, ccall(xs__)), args)

    # use the function name in the error message
    wrapup_ccall = quote
        if ec != Integer(STATUS_OK)
            msg = string("Error in ", $(name), " ", Integer(ec))
            throw(BrainFlowError(msg, Integer(ec)))
        end
    end

    # wrapping around the original ccall
    new_body = quote
        $(args[1:ccall_index-1]...)
        ec = Integer(STATUS_OK)
        ec = $(args[ccall_index])
        $wrapup_ccall
        $(args[ccall_index+1:end]...)
    end

    def[:body] = new_body
    esc(MacroTools.combinedef(def))
end

max_stars_repo_name: p-riebs/brainflow

max_issues_repo_name: p-riebs/brainflow

max_forks_repo_name: p-riebs/brainflow


text:
export ExtrinsicSubspaceHomotopy

"""
    ExtrinsicSubspaceHomotopy
"""
mutable struct ExtrinsicSubspaceHomotopy{T<:AbstractSystem} <: AbstractHomotopy
    F::T
    start::LinearSubspace{ComplexF64}
    target::LinearSubspace{ComplexF64}
    #cache
    A::Matrix{ComplexF64}
    b::Vector{ComplexF64}
    Ȧ::Matrix{ComplexF64}
    ḃ::Vector{ComplexF64}
    ū::Vector{ComplexDF64}
    t_cache::Base.RefValue{ComplexF64}
end

function ExtrinsicSubspaceHomotopy(
    F::System,
    start,
    target;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
)
    ExtrinsicSubspaceHomotopy(fixed(F; compile = compile), start, target)
end

function ExtrinsicSubspaceHomotopy(
    F::AbstractSystem,
    start::LinearSubspace,
    target::LinearSubspace;
    kwargs...,
)
    ExtrinsicSubspaceHomotopy(
        F,
        convert(LinearSubspace{ComplexF64}, start),
        convert(LinearSubspace{ComplexF64}, target);
        kwargs...,
    )
end

function ExtrinsicSubspaceHomotopy(
    F::AbstractSystem,
    start::LinearSubspace{ComplexF64},
    target::LinearSubspace{ComplexF64};
    kwargs...,
)
    A = copy(extrinsic(start).A)
    b = copy(extrinsic(start).b)

    Ȧ = extrinsic(start).A - extrinsic(target).A
    ḃ = extrinsic(start).b - extrinsic(target).b
    ū = zeros(ComplexDF64, size(A, 1))
    ExtrinsicSubspaceHomotopy(F, start, target, A, b, Ȧ, ḃ, ū, Ref(complex(NaN)))
end

Base.size(H::ExtrinsicSubspaceHomotopy) = (first(size(H.F)) + size(H.A, 1), last(size(H.F)))

start_parameters!(H::ExtrinsicSubspaceHomotopy, start) = parameters!(H, start, H.target)
target_parameters!(H::ExtrinsicSubspaceHomotopy, target) = parameters!(H, H.start, target)
parameters!(H::ExtrinsicSubspaceHomotopy, p, q) = set_subspaces!(H, p, q)
function set_subspaces!(
    H::ExtrinsicSubspaceHomotopy,
    start::LinearSubspace,
    target::LinearSubspace,
)
    H.start = start
    H.target = target
    H.Ȧ .= extrinsic(start).A .- extrinsic(target).A
    H.ḃ .= extrinsic(start).b .- extrinsic(target).b
    H.t_cache[] = NaN
    H
end

function tAb!(H::ExtrinsicSubspaceHomotopy, t)
    t == H.t_cache[] && return nothing

    A₁ = extrinsic(H.start).A
    A₀ = extrinsic(H.target).A
    b₁ = extrinsic(H.start).b
    b₀ = extrinsic(H.target).b
    if imag(t) == 0
        let t = real(t), t0 = (1.0 - real(t))
            @inbounds for i in eachindex(H.A)
                H.A[i] = t * A₁[i] + t0 * A₀[i]
            end
            @inbounds for i in eachindex(H.b)
                H.b[i] = t * b₁[i] + t0 * b₀[i]
            end
        end
    else
        t0 = (1.0 - t)
        @inbounds for i in eachindex(H.A)
            H.A[i] = t * A₁[i] + t0 * A₀[i]
        end
        @inbounds for i in eachindex(H.b)
            H.b[i] = t * b₁[i] + t0 * b₀[i]
        end
    end
    H.t_cache[] = t

    nothing
end

function ModelKit.evaluate!(u, H::ExtrinsicSubspaceHomotopy, x, t)
    tAb!(H, t)
    evaluate!(u, H.F, x)
    m = first(size(H.F))
    for i = 1:length(H.b)
        u[m+i] = -H.b[i]
    end
    for j = 1:size(H.A, 2)
        xj = x[j]
        for i = 1:size(H.A, 1)
            u[m+i] = muladd(H.A[i, j], xj, u[m+i])
        end
    end
    u
end
function ModelKit.evaluate!(u, H::ExtrinsicSubspaceHomotopy, x::Vector{ComplexDF64}, t)
    tAb!(H, t)
    evaluate!(u, H.F, x)
    m = first(size(H.F))
    for i = 1:length(H.b)
        H.ū[i] = -H.b[i]
    end
    for j = 1:size(H.A, 2)
        xj = x[j]
        for i = 1:size(H.A, 1)
            H.ū[i] = muladd(H.A[i, j], xj, H.ū[i])
        end
    end
    for i = 1:length(H.b)
        u[m+i] = H.ū[i]
    end

    u
end

function ModelKit.evaluate_and_jacobian!(u, U, H::ExtrinsicSubspaceHomotopy, x, t)
    tAb!(H, t)
    evaluate_and_jacobian!(u, U, H.F, x)
    m = first(size(H.F))
    for i = 1:length(H.b)
        u[m+i] = -H.b[i]
    end
    for j = 1:size(H.A, 2)
        xj = x[j]
        for i = 1:size(H.A, 1)
            u[m+i] = muladd(H.A[i, j], xj, u[m+i])
        end
    end
    for j = 1:size(H.A, 2), i = 1:size(H.A, 1)
        U[m+i, j] = H.A[i, j]
    end
    u
end

function ModelKit.taylor!(u, v::Val{1}, H::ExtrinsicSubspaceHomotopy, tx, t)
    m = first(size(H.F))
    for i = 1:m
        u[i] = zero(eltype(u))
    end
    for i = 1:length(H.b)
        u[m+i] = -H.ḃ[i]
    end
    for j = 1:size(H.Ȧ, 2)
        xⱼ, = tx[j]
        for i = 1:size(H.Ȧ, 1)
            u[m+i] = muladd(H.Ȧ[i, j], xⱼ, u[m+i])
        end
    end
    u
end

function ModelKit.taylor!(
    u,
    v::Val{K},
    H::ExtrinsicSubspaceHomotopy,
    tx::TaylorVector{K},
    t,
) where {K}
    m = first(size(H.F))
    taylor!(u, v, H.F, tx)
    for i = 1:length(H.b)
        u[m+i] = 0.0
    end
    for j = 1:size(H.Ȧ, 2)
        xⱼ = tx[j, K]
        for i = 1:size(H.Ȧ, 1)
            u[m+i] = muladd(H.Ȧ[i, j], xⱼ, u[m+i])
        end
    end
    u
end


text:
using Documenter, LogExpFunctions

# Doctest setup
DocMeta.setdocmeta!(
    LogExpFunctions, :DocTestSetup, :(using LogExpFunctions); recursive=true
)

makedocs(;
    modules=[LogExpFunctions],
    format=Documenter.HTML(),
    sitename="LogExpFunctions.jl",
    pages=Any["index.md"],
    checkdocs=:exports,
    strict=true,
)

deploydocs(; repo="github.com/JuliaStats/LogExpFunctions.jl.git")


text:
using Images
using ImageView
using LinearAlgebra

A = load("../../DATA/dog.jpg")
X = Gray.(A) # Convert RBG -> gray
ny, nx = size(X)

U, S, V = svd(X)

for r = [5, 20, 100] # Truncation value
    Xapprox = U[:, 1:r] * diagm(S)[1:r, 1:r] * V[:, 1:r]' # Approx. image
    imshow(Xapprox)
end


text:
using DataFrames
using MixedModels
using MixedModelsPermutations
using MixedModelsPermutations: resample!
using StableRNGs
using Test

isdefined(@__MODULE__, :io) || const io = IOBuffer()

@testset "LMM" begin
    sleepstudy = MixedModels.dataset(:sleepstudy)
    m1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy, progress=false)
    rm1 = resample!(StableRNG(42), deepcopy(m1))
    # test marking as not fit
    @test_logs (:warn,) show(io, rm1)

    non = nonparametricbootstrap(StableRNG(42), 1000, m1)
    @test non isa MixedModelBootstrap

    nondf = DataFrame(shortestcovint(non))

    # when scale inflation is correct, par and non should line up very closely.
    par = parametricbootstrap(StableRNG(42), 1000, m1)
    pardf = DataFrame(shortestcovint(par))

    let rho = filter(:type => ==("ρ"), nondf)
        violations = count(zip(rho.lower, rho.upper)) do (a, b)
            return !(-1 <= a <= b <= +1)
        end

        @test violations == 0
    end

    let beta = filter(:type => ==("β"), nondf)
        violations = count(zip(fixef(m1), beta.lower, beta.upper)) do (b, lower, upper)
            return !(lower <= b <= upper)
        end
        @test violations == 0
    end

    let sigma = filter(:type => ==("σ"), nondf)
        sigs = [ MixedModels.σs(m1).subj...; m1.σ ]

        violations = count(zip(sigs, sigma.lower, sigma.upper)) do (b, lower, upper)
            return !(0 <= lower <= b <= upper)
        end
        @test violations == 0
    end
end

@testset "GLMM" begin
    cbpp = MixedModels.dataset(:cbpp)
    gm1 = fit(MixedModel, @formula((incid/hsz) ~ 1 + period + (1|herd)),
              cbpp, Binomial(); wts=cbpp.hsz, fast=true, progress=false)

    @test_throws MethodError nonparametricbootstrap(1, gm1)
    @test_throws MethodError resample!(gm1)
end