text:
"""
Bindings indicate that an `EXPR` _may_ introduce a new name into the current scope/namespace.
Struct fields:
* `name`: the `EXPR` that defines the unqualifed name of the binding.
* `val`: what the binding points to, either a `Binding` (indicating ..), `EXPR` (this is generally the expression that defines the value) or `SymStore`.
* `type`: the type of the binding, either a `Binding`, `EXPR`, or `SymStore`.
* `refs`: a list containing all references that have been made to the binding.
"""
mutable struct Binding
    name::EXPR
    val::Union{Binding,EXPR,SymbolServer.SymStore,Nothing}
    type::Union{Binding,SymbolServer.SymStore,Nothing}
    refs::Vector{Any}
end
Binding(x::EXPR) = Binding(CSTParser.get_name(x), x, nothing, [])

function Base.show(io::IO, b::Binding)
    printstyled(io, " Binding(", to_codeobject(b.name),
        b.type === nothing ? "" : ":: ",
        b.refs isa Vector ? "($(length(b.refs)) refs))" : ")", color=:blue)
end

hasbinding(x::EXPR) = hasmeta(x) && hasbinding(x.meta)
bindingof(x) = nothing
bindingof(x::EXPR) = bindingof(x.meta)

hasref(x::EXPR) = hasmeta(x) && hasref(x.meta)
refof(x::EXPR) = hasmeta(x) ? x.meta.ref : nothing

function gotoobjectofref(x::EXPR)
    r = refof(x)
    if r isa SymbolServer.SymStore
        return r
    elseif r isa Binding

    end
end

"""
    mark_bindings!(x::EXPR, state)

Checks whether the expression `x` should introduce new names and marks them as needed. Generally this marks expressions that would introdce names to the current scope (i.e. that x sits in) but in cases marks expressions that will add names to lower scopes. This is done when it is not knowable that a child node of `x` will introduce a new name without the context of where it sits in `x` -for example the arguments of the signature of a function definition.
"""
function mark_bindings!(x::EXPR, state)
    if hasbinding(x)
        return
    end
    if !hasmeta(x)
        x.meta = Meta()
    end
    if isassignment(x)
        if CSTParser.is_func_call(x.args[1])
            name = CSTParser.get_name(x)
            mark_binding!(x)
            mark_sig_args!(x.args[1])
        elseif CSTParser.iscurly(x.args[1])
            mark_typealias_bindings!(x)
        elseif !is_getfield(x.args[1]) && state.flags & NO_NEW_BINDINGS == 0
            mark_binding!(x.args[1], x)
        end
    elseif CSTParser.defines_anon_function(x)
        mark_binding!(x.args[1], x)
    elseif CSTParser.iswhere(x)
        for i = 2:length(x.args)
            mark_binding!(x.args[i])
        end
    elseif headof(x) === :for
        markiterbinding!(x.args[2])
    elseif headof(x) === :generator || headof(x) === :filter
        for i = 2:length(x.args)
            markiterbinding!(x.args[i])
        end
    elseif headof(x) === :do
        for i in 1:length(x.args[2].args)
            mark_binding!(x.args[2].args[i])
        end
    elseif headof(x) === :function || headof(x) === :macro
        name = CSTParser.get_name(x)
        x.meta.binding = Binding(name, x, CoreTypes.Function, [])
        if isidentifier(name) && headof(x) === :macro
            setref!(name, bindingof(x))
        end
        mark_sig_args!(CSTParser.get_sig(x))
    elseif CSTParser.defines_module(x)
        x.meta.binding = Binding(x.args[2], x, CoreTypes.Module, [])
        setref!(x.args[2], bindingof(x))
    elseif headof(x) === :try && isidentifier(x.args[2])
        mark_binding!(x.args[2])
        setref!(x.args[2], bindingof(x.args[2]))
    elseif CSTParser.defines_datatype(x)
        name = CSTParser.get_name(x)
        x.meta.binding = Binding(name, x, CoreTypes.DataType, [])
        kwdef = parentof(x) isa EXPR && _points_to_Base_macro(parentof(x).args[1], Symbol("@kwdef"), state)
        if isidentifier(name)
            setref!(name, bindingof(x))
        end
        mark_parameters(CSTParser.get_sig(x))
        if CSTParser.defines_struct(x) # mark field block
            for i in 1:length(x.args[3].args)
                CSTParser.defines_function(x.args[3].args[i]) && continue
                arg = x.args[3].args[i]
                if kwdef && CSTParser.isassignment(arg)
                    arg = arg.args[1]
                end
                mark_binding!(arg)
            end
        end
    elseif headof(x) === :local
        for i = 1:length(x.args)
            if isidentifier(x.args[i])
                mark_binding!(x.args[i])
                setref!(x.args[i], bindingof(x.args[i]))
            end
        end
    end
end

function mark_binding!(x::EXPR, val=x)
    if CSTParser.iskwarg(x) || (CSTParser.isdeclaration(x) && CSTParser.istuple(x.args[1]))
        mark_binding!(x.args[1], x)
    elseif CSTParser.istuple(x) || CSTParser.isparameters(x)
        for arg in x.args
            mark_binding!(arg, val)
        end
    elseif CSTParser.isbracketed(x)
        mark_binding!(CSTParser.rem_invis(x), val)
    elseif CSTParser.issplat(x)
        mark_binding!(x.args[1], x)
    elseif !(isunarysyntax(x) && valof(headof(x)) == "::")
        if !hasmeta(x)
            x.meta = Meta()
        end
        x.meta.binding = Binding(CSTParser.get_name(x), val, nothing, [])
    end
    return x
end

function mark_parameters(sig::EXPR, params = String[])
    if CSTParser.issubtypedecl(sig)
        mark_parameters(sig.args[1], params)
    elseif iswhere(sig)
        for i = 2:length(sig.args)
            x = mark_binding!(sig.args[i])
            val = valof(bindingof(x).name)
            if val isa String
                push!(params, val)
            end
        end
        mark_parameters(sig.args[1], params)
    elseif CSTParser.iscurly(sig)
        for i = 2:length(sig.args)
            x = mark_binding!(sig.args[i])
            if bindingof(x) isa Binding && valof(bindingof(x).name) in params
                # Don't mark a new binding if a parameter has already been
                # introduced from a :where
                x.meta.binding = nothing
            end
        end
    end
    sig
end

function markiterbinding!(iter::EXPR)
    if CSTParser.isassignment(iter)
        mark_binding!(iter.args[1], iter)
    elseif CSTParser.iscall(iter) && CSTParser.isoperator(iter.args[1]) && (valof(iter.args[1]) == "in" || valof(iter.args[1]) == "∈")
        mark_binding!(iter.args[2], iter)
    elseif headof(iter) === :block
        for i = 1:length(iter.args)
            markiterbinding!(iter.args[i])
        end
    end
    return iter
end

function mark_sig_args!(x::EXPR)
    if CSTParser.iscall(x) || CSTParser.istuple(x)
        if x.args !== nothing && length(x.args) > 0
            if CSTParser.isbracketed(x.args[1]) && length(x.args[1].args) > 0 && CSTParser.isdeclaration(x.args[1].args[1])
                mark_binding!(x.args[1].args[1])
            end
            for i = (CSTParser.iscall(x) ? 2 : 1):length(x.args)
                a = x.args[i]
                if CSTParser.isparameters(a)
                    for j = 1:length(a.args)
                        aa = a.args[j]
                        mark_binding!(aa)
                    end
                elseif CSTParser.ismacrocall(a) && CSTParser.isidentifier(a.args[1]) && valofid(a.args[1]) == "@nospecialize" && length(a.args) == 3
                    mark_binding!(a.args[3])
                else
                    mark_binding!(a)
                end
            end
        end
    elseif CSTParser.iswhere(x)
        for i in 2:length(x.args)
            mark_binding!(x.args[i])
        end
        mark_sig_args!(x.args[1])
    elseif CSTParser.isbracketed(x)
        mark_sig_args!(x.args[1])
    elseif CSTParser.isdeclaration(x)
        mark_sig_args!(x.args[1])
    elseif CSTParser.isbinarycall(x)
        mark_binding!(x.args[1])
        mark_binding!(x.args[2])
    elseif CSTParser.isunarycall(x) && length(x.args) == 2 && (CSTParser.isbracketed(x.args[2]) || CSTParser.isdeclaration(x.args[2]))
        mark_binding!(x.args[2])
    end
end

function mark_typealias_bindings!(x::EXPR)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.binding = Binding(CSTParser.get_name(x.args[1]), x, CoreTypes.DataType, [])
    setscope!(x, Scope(x))
    for i = 2:length(x.args[1].args)
        arg = x.args[1].args[i]
        if isidentifier(arg)
            mark_binding!(arg)
        elseif CSTParser.issubtypedecl(arg) && isidentifier(arg.args[1])
            mark_binding!(arg.args[1])
        end
    end
    return x
end

function is_in_funcdef(x)
    if !(parentof(x) isa EXPR)
        return false
    elseif CSTParser.iswhere(parentof(x)) || CSTParser.isbracketed(parentof(x))
        return is_in_funcdef(parentof(x))
    elseif headof(parentof(x)) === :function || CSTParser.isassignment(parentof(x))
        return true
    else
        return false
    end
end

rem_wheres_subs_decls(x::EXPR) = (iswhere(x) || isdeclaration(x) || CSTParser.issubtypedecl(x)) ? rem_wheres_subs_decls(x.args[1]) : x

function _in_func_or_struct_def(x::EXPR)
    # only called in :where
    # check 1st arg contains a call (or op call)
    ex = rem_wheres_subs_decls(x.args[1])
    is_in_fexpr(x, CSTParser.defines_struct) || ((CSTParser.iscall(ex) || CSTParser.is_getfield(ex) || CSTParser.isunarycall(ex)) && is_in_funcdef(x))
end

"""
    add_binding(x, state, scope=state.scope)

Add the binding of `x` to the current scope. Special handling is required for:
* macros: to prefix the `@`
* functions: These are added to the top-level scope unless this syntax is used to define a closure within a function. If a function with the same name already exists in the scope then it is not replaced. This enables the `refs` list of the Binding of that 'root method' to hold a method table, the name of the new function will resolve to the binding of the root method (to get a list of actual methods -`[get_method(ref) for ref in binding.refs if get_method(ref) !== nothing]`). For example
```julia
[1] f() = 1
[2] f(x) = 2
```
[1] is the root method and the name of [2] resolves to the binding of [1]. Functions declared with qualified names require special handling, there are comments in the source.

Some simple type inference is run.
"""
function add_binding(x, state, scope=state.scope)
    if bindingof(x) isa Binding
        b = bindingof(x)
        if isidentifier(b.name)
            name = valofid(b.name)
        elseif CSTParser.ismacroname(b.name) # must be getfield
            name = string(to_codeobject(b.name))
        elseif isoperator(b.name)
            name = valof(b.name)
        else
            return
        end
        # check for global marker
        if isglobal(name, scope)
            scope = _get_global_scope(state.scope)
        end

        if CSTParser.defines_macro(x)
            scope.names[string("@", name)] = b
            mn = CSTParser.get_name(x)
            if isidentifier(mn)
                setref!(mn, b)
            end
        elseif defines_function(x)
            # TODO: Need to do check that we're not in a closure.
            tls = retrieve_toplevel_or_func_scope(scope)
            tls === nothing && return @warn "top-level scope not retrieved"
            if name_is_getfield(b.name)
                resolve_ref(parentof(parentof(b.name)).args[1], scope, state)
                lhs_ref = refof_maybe_getfield(parentof(parentof(b.name)).args[1])
                if lhs_ref isa SymbolServer.ModuleStore && haskey(lhs_ref.vals, Symbol(name))
                    # Overloading
                    if haskey(tls.names, name) && eventually_overloads(tls.names[name], lhs_ref.vals[Symbol(name)], state)
                        # Though we're explicitly naming a function for overloading, it has already been imported to the toplevel scope.
                        if !hasref(b.name)
                            setref!(b.name, tls.names[name]) # Add ref to previous overload
                            overload_method(tls, b, VarRef(lhs_ref.name, Symbol(name)))
                        end
                        # Do nothing, get_name(x) will resolve to the root method
                    elseif isexportedby(name, lhs_ref)
                        # Name is already available
                        tls.names[name] = b
                        if !hasref(b.name) # Is this an appropriate indicator that we've not marked the overload?
                            push!(b.refs, maybe_lookup(lhs_ref[Symbol(name)], state))
                            setref!(b.name, b) # we actually set the rhs of the qualified name to point to this binding
                        end
                    else
                        # Mark as overloaded so that calls to `M.f()` resolve properly.
                        overload_method(tls, b, VarRef(lhs_ref.name, Symbol(name))) # Add to overloaded list but not scope.
                    end
                elseif lhs_ref isa Binding && CoreTypes.ismodule(lhs_ref.type)
                    if hasscope(lhs_ref.val) && haskey(scopeof(lhs_ref.val).names, name)
                        # Don't need to do anything, name will resolve
                    end
                end
            else
                if scopehasbinding(tls, name)

                    existing_binding = tls.names[name]
                    if existing_binding isa Binding && (existing_binding.val isa Binding || existing_binding.val isa SymbolServer.FunctionStore || existing_binding.val isa SymbolServer.DataTypeStore)
                        # Should possibly be a while statement
                        # If the .val is as above the Binding likely won't have a proper type attached
                        # so lets use the .val instead.
                        existing_binding = existing_binding.val
                    end
                    if (existing_binding isa Binding && ((CoreTypes.isfunction(existing_binding.type) || CoreTypes.isdatatype(existing_binding.type))) || existing_binding isa SymbolServer.FunctionStore || existing_binding isa SymbolServer.DataTypeStore)
                        # do nothing name of `x` will resolve to the root method
                    else
                        seterror!(x, CannotDefineFuncAlreadyHasValue)
                    end
                else
                    scope.names[name] = b
                    if !hasref(b.name)
                        setref!(b.name, b)
                    end
                end
                if CSTParser.defines_struct(scope.expr) && parentof(scope) isa Scope
                    # hoist binding for inner constructor to parent scope
                    return add_binding(x, state, parentof(scope))
                end
            end
        elseif scopehasbinding(scope, name)
            # TODO: some checks about rebinding of consts
            check_const_decl(name, b, scope)

            scope.names[name] = b
        elseif is_soft_scope(scope) && parentof(scope) isa Scope && isidentifier(b.name) && scopehasbinding(parentof(scope), valofid(b.name)) && !enforce_hard_scope(x, scope)
            add_binding(x, state, scope.parent)
        else
            scope.names[name] = b
        end
        infer_type(b, scope, state)
    elseif bindingof(x) isa SymbolServer.SymStore
        scope.names[valofid(x)] = bindingof(x)
    end
end

function enforce_hard_scope(x::EXPR, scope)
    scope.expr.head === :for && is_in_fexpr(x, x-> x == scope.expr.args[1])
end

name_is_getfield(x) = parentof(x) isa EXPR && parentof(parentof(x)) isa EXPR && CSTParser.is_getfield_w_quotenode(parentof(parentof(x)))

"""
eventually_overloads(b, x, state)

"""
eventually_overloads(b::Binding, ss::SymbolServer.SymStore, state) = b.val == ss || (b.refs !== nothing && length(b.refs) > 0 && first(b.refs) == ss)
eventually_overloads(b::Binding, ss::SymbolServer.VarRef, state) = eventually_overloads(b, maybe_lookup(ss, state), state)
eventually_overloads(b, ss, state) = false

isglobal(name, scope) = false
isglobal(name::String, scope) = scope !== nothing && scopehasbinding(scope, "#globals") && name in scope.names["#globals"].refs

function mark_globals(x::EXPR, state)
    if headof(x) === :global
        if !scopehasbinding(state.scope, "#globals")
            state.scope.names["#globals"] = Binding(EXPR(:IDENTIFIER, EXPR[], nothing, 0, 0, "#globals", nothing, nothing), nothing, nothing, [])
        end
        for i = 2:length(x.args)
            if isidentifier(x.args[i]) && !scopehasbinding(state.scope, valofid(x.args[i]))
                push!(state.scope.names["#globals"].refs, valofid(x.args[i]))
            end
        end
    end
end

function name_extends_imported_method(b::Binding)
    if CoreTypes.isfunction(b.type) && CSTParser.hasparent(b.name) && CSTParser.is_getfield(parentof(b.name))
        if refof_maybe_getfield(parentof(b.name)[1]) !== nothing

        end
    end
end


text:
export get_renamed_quantifiers_form

"""
    get_renamed_quantifiers_form(formula, replacements=Dict{Variable, Variable}())

Return `formula` with new, unique variable symbols assigned to each quantified variable.

The returned formula is logically equivalent to `formula`.
"""
function get_renamed_quantifiers_form(formula, replacements=Dict{Variable, Variable}())
    @match formula begin
    f::Negation => Negation(
        get_renamed_quantifiers_form(f.formula, copy(replacements)),
    )
    f::Conjunction || f::Disjunction => typeof(f)(
        get_renamed_quantifiers_form(f.formula1, copy(replacements)),
        get_renamed_quantifiers_form(f.formula2, copy(replacements)),
    )
    f::AFormula || f::EFormula => begin
        replacements[f.variable] = Variable(getnextfreesymbol())
        typeof(f)(
            replacements[f.variable],
            get_renamed_quantifiers_form(f.formula, copy(replacements)),
        )
    end
    f::AtomicFormula || f::Function => begin
        typeof(f)(
            f.name,
            get_renamed_quantifiers_form.(f.arguments, copy(replacements)),
        )
    end
    f::Variable => f in keys(replacements) ? replacements[f] : f
    end
end

ext: jl


text:
"""
    AbstractSystemModel{T, N}

An abstract type that represents a set of parameters of type `T` for a dynamical system consisting
of `N` variables.

Concrete subtypes must implement the following methods:

- **[`get_dt`](@ref)**: Returns the time step.
- **[`get_ui`](@ref)**: Returns the initial condition.
- **[`get_nvars`](@ref)**: Returns the number of variables.
- **[`interaction_matrix`](@ref)**: Returns the interaction matrix for the model (respecting the
    coupling parameter for that particular model).

Concrete subtypes *may* implement the following methods:

- **[`rand`](@ref)**. If implemented for a discrete system model of type `SM`,
    returns an instance of the model with randomised parameters.
"""
abstract type AbstractSystemModel{T, N} end

""" Return the integration time step for a system model. """
get_dt(x::AbstractSystemModel) = x.dt

""" Return the initial condition for a system model. """
get_ui(x::AbstractSystemModel) = x.ui

""" Return the number of variables for a system model. """
get_nvars(x::AbstractSystemModel{T, N}) where {T, N} = N

"""
    rand(::Type{AbstractSystemModel{T, N}}) where {T, N}

Return a instance of the model with randomised parameters. Optionally implemented for
some systems.
"""
function rand(::Type{SM}) where SM <: ASM where ASM <: AbstractSystemModel
    @error "`rand` not implemented for system model of type `SM`. Consider defining `rand(::Type{$(SM)})`"
end

function display_string(x::AbstractSystemModel{T, N}) where {T, N}
    tp = typeof(x)
    fn = fieldnames(typeof(x))

    "$tp" * "(" * join(["\n  $(fd) = $(getfield(x, fd))" for fd in fn]) * ")"
end

function Base.show(io::IO, x::AbstractSystemModel{T, N}) where {T, N}
    print(io, display_string(x))
end

export get_dt, get_ui, get_nvars, randomised

"""
    rand(::Type{SM}) where {SM <: AbstractSystemModel} -> SM

Generate an instance of a model of type `SM` with randomised parameters.

Parameters are provided as keyword arguments, and must be either as scalars,
distributions or more complicated uncertain values defined as in the
[UncertainData](https://github.com/kahaaga/UncertainData.jl) package.

## See also

Concrete instructions to randomise the different systems are found in the
following docstrings.

- [`rand(::Type{RosslerLorenzUnidir}`].
"""
function rand(x::AbstractSystemModel) end

const PT = Union{Number, Distribution, AbstractUncertainValue}


text:
function cycle_parent_children(cyclic_members, parents)
    parent = -1
    for member in cyclic_members
        !isempty(parents[member]) && parents[member][1] ∈ cyclic_members && continue
        parent = member
        break
    end

    return parent, setdiff(cyclic_members, parent)
end

function lump_cycles!(cycles, cyclic_children)
    inds = Int64[]
    for (i, cycle) in enumerate(cycles)
        if cyclic_children ⊆ cycle
            return
        elseif cyclic_children ⊇ cycle
            append!(inds, i)
        end
    end
    push!(cycles, cyclic_children)
    deleteat!(cycles, inds)
    return
end


text:
"""
    shift_fields(layer, offset)

Adds `offset` to the `layer` fields.
"""
function shift_fields end

shift_fields(layer::Binary, offset::AbstractArray) = Binary(layer.θ + offset)
shift_fields(layer::Spin, offset::AbstractArray) = Spin(layer.θ + offset)
shift_fields(layer::Potts, offset::AbstractArray) = Potts(layer.θ + offset)
shift_fields(layer::Gaussian, offset::AbstractArray) = Gaussian(layer.θ + offset, layer.γ)
shift_fields(layer::ReLU, offset::AbstractArray) = ReLU(layer.θ + offset, layer.γ)

function shift_fields(layer::dReLU, offset::AbstractArray)
    return dReLU(layer.θp + offset, layer.θn + offset, layer.γp, layer.γn)
end

function shift_fields(layer::pReLU, offset::AbstractArray)
    return pReLU(layer.θ + offset, layer.γ, layer.Δ, layer.η)
end

function shift_fields(layer::xReLU, offset::AbstractArray)
    return xReLU(layer.θ + offset, layer.γ, layer.Δ, layer.ξ)
end

# In-place versions

"""
    shift_fields!(layer, offset)

In-place version of `shift_fields(layer, offset)`.
"""
function shift_fields! end

function shift_fields!(
    layer::Union{Binary, Spin, Potts, Gaussian, ReLU, pReLU, xReLU}, offset::AbstractArray)
    layer.θ .+= offset
    return layer
end

function shift_fields!(layer::dReLU, offset::AbstractArray)
    layer.θp .+= offset
    layer.θn .+= offset
    return layer
end


text:
export ModiaProblem, ModiaSolve

function ModiaProblem(m1::ModiaLang.SimulationModel{FloatType1,ParType,EvaluatedParType,FloatType1},
                      m2::ModiaLang.SimulationModel{FloatType2,ParType,EvaluatedParType,FloatType2};
                      p, merge=nothing, kwargs...) where {FloatType1,ParType,EvaluatedParType,FloatType2}

    TimeType1 = FloatType1
    TimeType2 = FloatType2
    merge2 = Map(p = parameter | convert(Vector{FloatType2}, p))
    if !isnothing(merge)
        merge2 = merge | merge2
    end

    # Store states x of m1 at communication points in the DifferentialEquations result data structure solution
    m1.save_x_in_solution = true
    m2.save_x_in_solution = true

    # Save kwargs of ModiaProblem in simulation models
    empty!(m1.result)
    empty!(m2.result)
    options1 = SimulationOptions{FloatType1,TimeType1}(merge ; kwargs...)
    options2 = SimulationOptions{FloatType2,TimeType2}(merge2; kwargs...)
    if isnothing(options1) || isnothing(options2)
        return nothing
    end
    m1.options = options1
    m2.options = options2

   	tspan = (options1.startTime, options1.stopTime)
    tspan_outputs = options1.startTime:options1.interval:options1.stopTime

    # Initialize/re-initialize SimulationModel
    if options1.log || options1.logParameters || options1.logEvaluatedParameters || options1.logStates
        println("... Construct simulation problem of DifferentialEquations.jl for ", m.modelName)
    end

    success1 = ModiaLang.init!(m1)
    if !success1|| m1.eventHandler.restart == ModiaLang.Terminate
        return nothing
    end

    success2 = ModiaLang.init!(m2)
    if !success2|| m2.eventHandler.restart == ModiaLang.Terminate
        return nothing
    end

    function p_derivatives!(der_x, x, p, t::FloatType1)::Nothing
        m1.evaluatedParameters.p .= p
        Base.invokelatest(m1.getDerivatives!, der_x, x, m1, t)
        return nothing
    end

    function p_derivatives!(der_x, x, p, t::FloatType2)::Nothing
        m2.evaluatedParameters.p .= p
        Base.invokelatest(m2.getDerivatives!, der_x, x, m2, t)
        return nothing
    end

    #= This call back gives an error:
    # type TrackedAffect has no field funciter
    function outputs!(x, t, integrator)::Nothing
       if eltype(integrator.p) == FloatType1
           m1.storeResult = true
           p_derivatives!(m1.der_x, x, integrator.p, t)
           m1.storeResult = false
       end
       return nothing
    end
    callback = DifferentialEquations.FunctionCallingCallback(outputs!, funcat=tspan_outputs)
    =#

    function affect_outputs!(integrator)::Nothing
       if eltype(integrator.p) == FloatType1
           m1.storeResult = true
           p_derivatives!(m1.der_x, integrator.u, integrator.p, integrator.t)
           m1.storeResult = false
       end
       return nothing
    end
    callback = DifferentialEquations.PresetTimeCallback(tspan_outputs, affect_outputs!)

    # Define problem and callbacks based on algorithm and model type
    abstol = 0.1*options1.tolerance
    return DifferentialEquations.ODEProblem(p_derivatives!, m1.x_init, tspan, p;   # x_init
             reltol=options1.tolerance, abstol=abstol, callback = callback,
             modia_interval = options1.interval,
             modia_instantiatedModel = m1)
end

function ModiaSolve(problem, algorithm=missing; p, adaptive::Bool = true)

    startTime = problem.tspan[1]
    stopTime  = problem.tspan[2]
    interval  = problem.kwargs.data[:modia_interval]
    m         = problem.kwargs.data[:modia_instantiatedModel]
    m.algorithmType = typeof(algorithm)

    # Define problem and callbacks based on algorithm and model type
    if abs(interval) < abs(stopTime-startTime)
        tspan = startTime:interval:stopTime
    else
        tspan = [startTime, stopTime]
    end

    # Initial step size (the default of DifferentialEquations is too large) + step-size of fixed-step algorithm
    dt = adaptive ? interval/10 : interval    # initial step-size

    # Compute solution
    solution = ismissing(algorithm) ? DifferentialEquations.solve(problem, p=p,
                                        saveat = tspan, adaptive=adaptive, dt=dt) :
                                    DifferentialEquations.solve(problem, algorithm, p=p,
                                        saveat = tspan, adaptive=adaptive, dt=dt)
    m.solution = solution
    return solution
end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
struct SymEGRSSMatrix{T,UT<:AbstractArray,VT<:AbstractArray} <: AbstractMatrix{T}
	Ut::UT
	Vt::VT
	n::Int
	p::Int

	function SymEGRSSMatrix{T,UT,VT}(Ut,Vt,n,p) where
				{T,UT<:AbstractArray,VT<:AbstractArray}
		Un, Um = size(Ut)
		Vn, Vm = size(Vt)
		(Un == Vn && Um == Vm) || throw(DimensionMismatch())
		new(Ut,Vt,n,p)
	end
end

SymEGRSSMatrix(Ut::AbstractArray{T,N},Vt::AbstractArray{T,N}) where {T,N}  =
	SymEGRSSMatrix{T,typeof(Ut),typeof(Vt)}(Ut,Vt,size(Ut)[2],size(Ut)[1])

########################################################################
#### Helpful properties. Not nessecarily computionally efficient    ####
########################################################################
Matrix(K::SymEGRSSMatrix) = tril(K.Ut'*K.Vt) + triu(K.Vt'*K.Ut,1)

size(K::SymEGRSSMatrix) = (K.n,K.n)

function getindex(K::SymEGRSSMatrix{T}, i::Int, j::Int) where T
	i > j && return dot(K.Ut[:,i],K.Vt[:,j])
	return dot(K.Vt[:,i],K.Ut[:,j])
end

Base.propertynames(F::SymEGRSSMatrix, private::Bool=false) =
    (private ? fieldnames(typeof(F)) : ())

########################################################################
#### Linear Algebra routines 							   			####
########################################################################
function symegrss_mul!(Y::AbstractVecOrMat{T}, K::SymEGRSSMatrix{Q,UT,VT},
		  X::AbstractVecOrMat{S}) where
		  {T,Q,UT<:AbstractArray,VT<:AbstractArray,S}
	Ut = K.Ut;
	Vt = K.Vt;
	n = K.n;
	m = K.p;
    mx = size(X,2);
    Vbar = zeros(m,mx);
    Ubar = Ut*X;
    @inbounds for i = 1:n
        tmpV = Vt[:,i];
        tmpU = Ut[:,i];
        Ubar -= tmpU .* X[i:i,:];
        Vbar += tmpV .* X[i:i,:];
        Y[i,:] = Vbar'*tmpU + Ubar'*tmpV;
    end
	return Y
end

logdet(K::SymEGRSSMatrix) = 2.0*logdet(cholesky(K))
det(K::SymEGRSSMatrix) = det(cholesky(K))^2

#### Matrix-matrix product ####
mul!(y::AbstractVecOrMat, K::SymEGRSSMatrix, x::AbstractVecOrMat) =
		symegrss_mul!(y,K,x)
mul!(y::AbstractVecOrMat, K::Adjoint{<:Any,<:SymEGRSSMatrix}, x::AbstractVecOrMat) =
		symegrss_mul!(y,K.parent,x)
function (\)(K::SymEGRSSMatrix, x::AbstractVecOrMat)
	L = cholesky(K);
	return L'\(L\x)
end

max_stars_repo_name: mipals/SymEGRSSMatrices


text:
for (M, op, T) in (
    (:Base, :+, :KernelSum),
    (:Base, :*, :KernelProduct),
    (:TensorCore, :tensor, :KernelTensorProduct),
)
    @eval begin
        $M.$op(k1::Kernel, k2::Kernel) = $T(k1, k2)

        $M.$op(k1::$T, k2::$T) = $T(k1.kernels..., k2.kernels...)
        function $M.$op(
            k1::$T{<:AbstractVector{<:Kernel}}, k2::$T{<:AbstractVector{<:Kernel}}
        )
            return $T(vcat(k1.kernels, k2.kernels))
        end

        $M.$op(k::Kernel, ks::$T) = $T(k, ks.kernels...)
        $M.$op(k::Kernel, ks::$T{<:AbstractVector{<:Kernel}}) = $T(vcat(k, ks.kernels))

        $M.$op(ks::$T, k::Kernel) = $T(ks.kernels..., k)
        $M.$op(ks::$T{<:AbstractVector{<:Kernel}}, k::Kernel) = $T(vcat(ks.kernels, k))
    end
end


text:
@testset "base" begin
    @testset "Trig" begin
        @testset "Basics" for x = (Float64(π)-0.01, Complex(π, π/2))
            test_scalar(sec, x)
            test_scalar(csc, x)
            test_scalar(cot, x)
            test_scalar(sinpi, x)
            test_scalar(cospi, x)
        end
        @testset "Hyperbolic" for x = (Float64(π)-0.01, Complex(π-0.01, π/2))
            test_scalar(sech, x)
            test_scalar(csch, x)
            test_scalar(coth, x)
        end
        @testset "Degrees" begin
            x = 45.0
            test_scalar(sind, x)
            test_scalar(cosd, x)
            test_scalar(tand, x)
            test_scalar(secd, x)
            test_scalar(cscd, x)
            test_scalar(cotd, x)
        end
        @testset "Inverses" for x = (0.5, Complex(0.5, 0.25))
            test_scalar(asec, 1/x)
            test_scalar(acsc, 1/x)
            test_scalar(acot, 1/x)
        end
        @testset "Inverse hyperbolic" for x = (0.5, Complex(0.5, 0.25),  Complex(-2.1 -3.1im))
            test_scalar(asinh, x)
            test_scalar(acosh, x + 1)  # +1 accounts for domain for real
            test_scalar(atanh, x)
            test_scalar(asech, x)
            test_scalar(acsch, x)
            test_scalar(acoth, x + 1)
        end

        @testset "Inverse degrees" for x = (0.5, Complex(0.5, 0.25))
            test_scalar(asind, x)
            test_scalar(acosd, x)
            test_scalar(atand, x)
            test_scalar(asecd, 1/x)
            test_scalar(acscd, 1/x)
            test_scalar(acotd, 1/x)
        end

        @testset "sinc" for x = (0.0, 0.434, Complex(0.434, 0.25))
            test_scalar(sinc, x)
        end
    end  # Trig

    @testset "Angles" begin
        for x in (-0.1, 6.4, 0.5 + 0.25im)
            test_scalar(deg2rad, x)
            test_scalar(rad2deg, x)
        end
    end

    @testset "Unary complex functions" begin
        for x in (-4.1, 6.4, 0.0, 0.0 + 0.0im, 0.5 + 0.25im)
            test_scalar(real, x)
            test_scalar(imag, x)
            test_scalar(hypot, x)
            test_scalar(adjoint, x)
        end
    end

    @testset "Complex" begin
        test_scalar(Complex, randn())
        test_scalar(Complex, randn(ComplexF64))

        test_frule(Complex, randn(), randn())
        test_rrule(Complex, randn(), randn())
    end

    @testset "*(x, y) (scalar)" begin
        # This is pretty important so testing it fairly heavily
        test_points = (0.0, -2.1, 3.2, 3.7+2.12im, 14.2-7.1im)
        @testset "($x) * ($y)" for
            x in test_points, y in test_points

            # ensure all complex if any complex for FiniteDifferences
            x, y = Base.promote(x, y)

            test_frule(*, x, y)
            test_rrule(*, x, y)
        end
    end

    @testset "ldexp" begin
        for n in (0,1,20)
            test_frule(ldexp, 10rand(), n)
            test_rrule(ldexp, 10rand(), n)
        end
    end

    @testset "\\(x::$T, y::$T) (scalar)" for T in (Float64, ComplexF64)
        test_frule(*, randn(T), randn(T))
        test_rrule(*, randn(T), randn(T))
    end

    @testset "mod" begin
        test_frule(mod, 10rand(), rand())
        test_rrule(mod, 10rand(), rand())
    end

    @testset "identity" for T in (Float64, ComplexF64)
        test_frule(identity, randn(T))
        test_frule(identity, randn(T, 4))
        test_frule(
            identity,
            #Tuple(randn(T, 3)) ⊢ Composite{Tuple{T, T, T}}(randn(T, 3)...)
            Tuple(randn(T, 3))
        )

        test_rrule(identity, randn(T))
        test_rrule(identity, randn(T, 4))
        test_rrule(
            identity,
            Tuple(randn(T, 3)) ⊢ Composite{Tuple{T, T, T}}(randn(T, 3)...);
            output_tangent = Composite{Tuple{T, T, T}}(randn(T, 3)...)
        )
    end

    @testset "Constants" for x in (-0.1, 6.4, 1.0+0.5im, -10.0+0im, 0.0+200im)
        test_scalar(one, x)
        test_scalar(zero, x)
    end

    @testset "muladd(x::$T, y::$T, z::$T)" for T in (Float64, ComplexF64)
        test_frule(muladd, 10randn(), randn(), randn())
        test_rrule(muladd, 10randn(), randn(), randn())
    end

    @testset "fma" begin
        test_frule(fma, 10randn(), randn(), randn())
        test_rrule(fma, 10randn(), randn(), randn())
    end

    @testset "clamp"  begin
        # to left
        test_frule(clamp, 1., 2., 3.)
        test_rrule(clamp, 1., 2., 3.)

        # in the middle
        test_frule(clamp, 2.5, 2., 3.)
        test_rrule(clamp, 2.5, 2., 3.)

        # to right
        test_frule(clamp, 4., 2., 3.)
        test_rrule(clamp, 4., 2., 3.)
    end

    @testset "rounding" begin
        for x in (-0.6, -0.2, 0.1, 0.6)
            # thanks to RoundNearest
            if 0 > x % 1 > 0.5 || -1 < x % 1 <= 0.5
                test_scalar(round, x; fdm=backward_fdm(5,1))
            else
                test_scalar(round, x; fdm=forward_fdm(5,1))
            end
            test_scalar(floor, x; fdm=backward_fdm(5, 1))
            test_scalar(ceil, x; fdm=forward_fdm(5, 1))
        end
    end

    @testset "type" begin
        @test frule((NO_FIELDS, DoesNotExist(), DoesNotExist()), typejoin, Array{Float32,4}, Array{Float32,3}) !== nothing
        @test rrule(typejoin, Array{Float32,4}, Array{Float32,3}) !== nothing
    end

    @testset "Logging" begin
        @test frule((NO_FIELDS, DoesNotExist(), DoesNotExist()), Base.depwarn, "message", :f) !== nothing
        @test rrule(Base.depwarn, "message", :f) !== nothing
    end
end


text:
@bitmask_flag FlagBits::UInt32 begin
    BIT_A = 1
    BIT_B = 2
    BIT_C = 3
end

@bitmask_flag FlagBits2::UInt32 begin
    BIT_A_2 = 1
    BIT_B_2 = 2
    BIT_C_2 = 3
end

@testset "Bitmasks" begin
    @testset "Bitmask creation & operations" begin
        @test BIT_A & BIT_B == FlagBits(0) == zero(FlagBits)
        @test iszero(BIT_A & BIT_B)
        @test BIT_A | BIT_B == BIT_C
        @test !iszero(BIT_C)
        @test xor(BIT_A, BIT_B) == BIT_C
        @test BIT_A < BIT_B && BIT_B < BIT_C && BIT_A < BIT_C
        @test BIT_A & 1 == 1 & BIT_A == BIT_A
        @test BIT_A | 2 == 2 | BIT_A == BIT_C
        @test xor(BIT_A, 2) == xor(2, BIT_A) == BIT_C
        @test 1 < BIT_B && BIT_B < 3
        @test BIT_A <= BIT_B <= BIT_C
        @test BIT_A == BIT_A
        @test BIT_B in BIT_C
        @test BIT_A ∉ BIT_B
        @test 2 in BIT_C
        @test BIT_B in 3
        @test_throws ErrorException BIT_A & BIT_A_2
        @test_throws ErrorException BIT_A == BIT_A_2
        @test_throws ErrorException BIT_A in BIT_A_2
    end

    all_graphics = SHADER_STAGE_VERTEX_BIT | SHADER_STAGE_FRAGMENT_BIT | SHADER_STAGE_TESSELLATION_CONTROL_BIT | SHADER_STAGE_TESSELLATION_EVALUATION_BIT | SHADER_STAGE_GEOMETRY_BIT

    @testset "Generation of Vulkan bitmasks" begin
        @test all_graphics == SHADER_STAGE_ALL_GRAPHICS
    end

    @testset "Correct printing of bitmasks" begin
        @test string(BIT_C) == "FlagBits(BIT_C)"
        @test string(all_graphics) == "ShaderStageFlag(SHADER_STAGE_ALL_GRAPHICS)"
        @test string(SHADER_STAGE_VERTEX_BIT | SHADER_STAGE_FRAGMENT_BIT) == "ShaderStageFlag(SHADER_STAGE_VERTEX_BIT | SHADER_STAGE_FRAGMENT_BIT)"
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
include("square.jl")
include("chimera.jl")
include("second_neighbor.jl")
include("cubic.jl")


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

using Markdown
using InteractiveUtils

# ╔═╡ dac66550-6a42-11eb-18fd-1dbb51ebb356
using DataFrames, PlutoUI

# ╔═╡ b9438ea0-6a42-11eb-3203-bdef1a59bb19
md"""
## Introduction to DataFrames / 2

ref. [Introduction to DataFrames, Bogumił Kamiński](https://github.com/bkamins/Julia-DataFrames-Tutorial/blob/master/02_basicinfo.ipynb)

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ be11851f-b9a5-4140-9d56-fa0e1ec896f9
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ 0c9061f5-92fc-4861-8920-acb7f16aea51
md"""
### Getting basic information about a DataFrame
"""

# ╔═╡ e8266538-6a42-11eb-0919-f7a19fd021ba
xdf = DataFrame(A = [1, 2], B = [1.0, missing], C = ["a", "b"])

# ╔═╡ f28d6da0-6a42-11eb-1bfd-edda225bd80a
size(xdf), size(xdf, 1), size(xdf, 2)

# ╔═╡ 02630c26-6a43-11eb-375b-d3bf7b722863
nrow(xdf), ncol(xdf)

# ╔═╡ 09ef1732-6a43-11eb-3ce5-6904ca52514f
describe(xdf)

# ╔═╡ 1436e9b8-6a43-11eb-1dde-d97c58bf5218
describe(xdf, cols=1:2)

# ╔═╡ 6ca3a755-9f7b-4b49-b15c-1989144e4e44
## A custom describe
describe(xdf, :eltype, :nmissing, :first => first, :last => last)

# ╔═╡ 20d87024-6a43-11eb-0a14-d5c00e2aafde
names(xdf)

# ╔═╡ 2c9fb4c6-6a43-11eb-3e4c-5d12a291df4a
propertynames(xdf)

# ╔═╡ 3534cd7e-6a43-11eb-2eb7-0fa46c82b12d
eltype.(eachcol(xdf))

# ╔═╡ 7f69dd67-3fa6-4c74-bf8c-4fecb087b74f
md"""
Let's create a large dataframe and have peek into it using `first` and `last`.
"""

# ╔═╡ 3f6d9596-6a43-11eb-3469-5dcd7f0332e0
begin
	n = 10
	ydf = DataFrame(rand(1:n, 1000, n))
end

# ╔═╡ 546cdee8-6a43-11eb-38a2-35fe89ae2dde
first(ydf, 3)

# ╔═╡ 5bc30e06-6a43-11eb-1625-015e8cb5bce1
last(ydf, 3)

# ╔═╡ cc4ffb73-ada8-4446-8914-9e1d5cf669d4
md"""
### DataFrames get and set operations

Re-using `xdf` DataFrame as defined a few cells above.
"""

# ╔═╡ c6786444-6a43-11eb-21b9-3d4a2722cdfd
# all get the vector stored in our DataFrame without copying it
xdf.A, xdf[!, 1], xdf[!, :A]

# ╔═╡ 04e69443-e80d-45ea-9d59-fb5dd8bf11e5
xdf."A", xdf[!, "A"]   # the same using string indexing

# ╔═╡ d105dcca-6a43-11eb-07fe-cd24a5efcf2a
begin
	xdf[:, 1]                                  # note that this creates a copy...
	xdf[:, 1] === xdf[:, 1]                    # ...here is the proof

	objectid(xdf[:, 1]), objectid(xdf[:, 1])   # ...and another way to show this
end

# ╔═╡ 49729e83-ebb2-4514-bdaf-68db8847800a
md"""
To grab one row as a DataFrame, we can index as follows:
"""

# ╔═╡ 0ea18045-6d7e-4202-ab94-709ab6150201
xdf[1:1, :]

# ╔═╡ 288c1879-77db-4246-b087-5f4e7ee600c5
xdf[1, :] # this produces a DataFrameRow which is treated as 1-dim object similar to a NamedTuple

# ╔═╡ 2140fccb-bbbb-4bae-b3e0-40a850ca40b9
typeof(xdf[1:1, :]), typeof(xdf[1, :])

# ╔═╡ b9566a8c-748b-4696-8aed-d9add8a085d0
md"""
We can grab a single cell or element with the same syntax to grab an element of an array.
"""

# ╔═╡ 0ea6ef6d-d908-4888-9fa3-4bd5edefeea6
xdf[1, 2]

# ╔═╡ 2aae7a6a-21db-45ee-bcdf-619ee2ea5240
md"""
or a new DataFrame that is a subset of rows and columns:
"""

# ╔═╡ 9b8d271c-5def-46d4-a38d-a98a63e4d62c
xdf[1:2, 1:2]

# ╔═╡ 2e1155ca-0cfc-4113-b2ce-2ff770e51e52
md"""
We can also use a `Regexp` to select columns and `Not` from `InvertedIndices.jl` both to select rows and columns>
"""

# ╔═╡ af6641c1-d1eb-4ae0-9e31-63aed4793b14
xdf[Not(1), r"A"]

# ╔═╡ bd1f96af-21fb-4e9e-933b-c4b99a05d13a
md"""
**Reference (`!`) and copy (`:`)**
"""

# ╔═╡ 6f4ab986-81af-4fbd-93b2-546b3ad7c425
xdf[!, Not(1)]     # ! indicates that underlying columns are not copied

# ╔═╡ da813272-f45a-4fcd-a4ff-539bb4c924c1
xdf[:, Not(1)]     # (however) : means that the columns will get copied

# ╔═╡ 17f8d99c-0bf0-4442-b234-d8d4955482e0
md"""
Assignment of a scalar to a dataframe can be done in ranges using broadcasting:
"""

# ╔═╡ 3be6e5f8-777b-4f9a-9430-362880e842ad
begin
	xdf[!, :B] = convert.(Float64, xdf[!, :B])   ## Type conversion for column :B
	xdf[1:2, :B] .= π;                           ## broadcaset π over all rows
	xdf
end

# ╔═╡ 137dfe1f-0a3d-4d76-978a-70da3c178479
begin
	xdf[1:2, 1:2] .= 3;   ## Assignment over both columns :A(1) and :B(1), over all rows
	xdf
end

# ╔═╡ 5f890d28-43ef-40f3-b1ff-b32c9cd5773d
md"""
Assignment of a vector of length equal to the number of assigned rows using broadcasting:
"""

# ╔═╡ 20ae267b-c26b-49c5-9096-96190f868c18
begin
	xdf[1:2, 1:2] .= [2, 4]
	xdf
end

# ╔═╡ 49daec8a-0263-4f69-a27b-dc0efdf89c71
md"""
Assignment or of another data frame of matching size and column names, again using broadcasting:
"""

# ╔═╡ ed5ee6d9-ec3e-417c-8206-f62666d88c4f
begin
	xdf[1:2, 1:2] .= DataFrame([5 6; 7 8], [:A, :B])
	xdf
end

# ╔═╡ 05d0a958-6071-41ba-8219-98f040ac90ae
md"""
#### Caution

  - With `df[!, :col]` and `df.col` syntax we get a direct (non copying) access to a column of a data frame. This is *potentially unsafe* as we can easily corrupt data in the `df` dataframe if we resize, sort, etc. the column obtained in this way. *Therefore such access should be used with caution*.

  - Similarly `df[!, cols]` when `cols` is a collection of columns produces a new data frame that holds the same (not copied) columns as the source `df` dataframe. *Similarly, modifying the data frame obtained via `df[!, cols]` might cause problems with the consistency of `df`*.

  - The `df[:, :col]` and `df[:, cols]` syntaxes always copy columns so they are safe to use (and should generally be preferred except for performance or memory critical use cases).
"""

# ╔═╡ c0960534-f6a1-40f8-8b35-b19aa5b78ac2
md"""
Here are examples of how `Cols` and `Between` can be used to select columns of a data frame:
"""

# ╔═╡ 0bfe2061-cddb-4f9f-9b7e-33a3e703d084
xdf₁ = DataFrame(rand(4, 5), :auto)

# ╔═╡ 5055355b-a405-4d27-b2aa-9477bd5ee2c4
ydf₀ = xdf₁[:, Between(:x2, :x4)]     ## make a copy

# ╔═╡ 2623b999-d8b1-4cba-9f2b-4d644c672354
ydf₁ = xdf₁[:, Cols("x1", Between("x3", "x4"))]   ## make a copy

# ╔═╡ 3c8b682a-5387-40b6-b3cd-d8beaa09383c
begin
	ydf₁[!, :x3] .= 3
	ydf₁
end

# ╔═╡ b1d6b7e3-8c67-44b3-ae30-87302295e976
xdf₁   ## was not changed

# ╔═╡ 61ac92e7-f181-4554-809a-dba4b786acec
ydf₀   ## nor this dataframe

# ╔═╡ e5e723bd-aa8e-415c-b7c6-ffd49bbf1841
md"""
#### Views

We can simply create a view of a DataFrame (it is more efficient than creating a materialized selection). Here are the possible return value options.
"""

# ╔═╡ fe37621d-c903-4055-bdb9-0dbc98f9d524
@view xdf[1:2, 1]

# ╔═╡ 8624b878-775b-4e37-a5c5-9b1cf38ef614
@view xdf[1, 1]

# ╔═╡ 9d40b847-5ebc-4656-aa03-f918ccee3974
@view xdf[1, 1:2]   ## a DataFrameRow, the same as for xdf[1, 1:2] without a view

# ╔═╡ 12a63294-6a44-11eb-0f29-7dcf33acface
md"""
### Adding new columns to a DataFrame
"""

# ╔═╡ 8c59bad5-6a7c-4923-b088-07a225c57a96
xdf₂ = DataFrame();    ## Create an empty DataFrame

# ╔═╡ 91e94361-3b46-4e0f-b8c3-14ced0b2e71b
md"""
Using `setproperty`!
"""

# ╔═╡ 48c59474-ffb8-4c92-96d1-39dddcd39076
begin
	v = collect(1:5)
	xdf₂.a = v
	xdf₂
end

# ╔═╡ 2180b2ba-3201-49a6-a0d5-5549c00c7b37
xdf₂.a === v   ## no copy is performed

# ╔═╡ 4c3feedd-f3e0-4121-91b6-9b43b0840c9b
md"""
using `setindex`!
"""

# ╔═╡ dcf94057-a263-4b27-b086-db0c74953405
begin
	xdf₂[!, :b] = v
	xdf₂[:, :c] = v

	@assert xdf₂[!, :b] === v  ## no copy is performed
	@assert xdf₂[!, :c] !== v  ## copy performed

	xdf₂
end

# ╔═╡ 56ba4b95-5a94-46be-ab6c-c80a40daccec
begin
	xdf₂[!, :d] .= v
	xdf₂[:, :e] .= v

	## using broadcasting implies copy:
	@assert xdf₂[!, :d] !== v  ## copy performed
	@assert xdf₂[!, :e] !== v  ## copy performed

	xdf₂
end

# ╔═╡ f0c302ae-23c1-4019-80d9-d727f83b9815
md"""
Remember that columns `:a` and `:b` are not copies of `v`.
This can lead to silent errors.

For example the following code leads to a bug (note that calling `pairs` on `eachcol(df)` creates an iterator of (column name, column) pairs):
"""

# ╔═╡ e0f1d7c5-039a-4be6-9b84-e7434faafad6
with_terminal() do
	for (n, c) in pairs(eachcol(xdf₂))
    	println("$n: ", pop!(c))           ## side-effect !
	end
end

# ╔═╡ 73b293e6-4609-419f-afe6-5f1185ab58f7
md"""
note (above cell) that for column `:b` we printed 4 as 5 was removed from it when we used `pop!` on column `:a`.

Such mistakes sometimes happen. Because of this `DataFrames.jl` performs consistency checks before doing an expensive operation (most notably before showing a DataFrame).
"""

# ╔═╡ 09290340-e74c-4f2c-9962-ccac5dd334a8
collect(pairs(eachcol(xdf₂)))  ## The output confirms that the DataFrame got corrupted.

# ╔═╡ 5a53e3f7-1585-4d16-ad4a-1e88ec08ef83
md"""
`DataFrames.jl` supports a complete set of `getindex`, `getproperty`, `setindex!`, `setproperty!`, view, broadcasting, and broadcasting assignment operations. The details are explained here: [http://juliadata.github.io/DataFrames.jl/latest/lib/indexing/](http://juliadata.github.io/DataFrames.jl/latest/lib/indexing/).
"""

# ╔═╡ 96f4c0ff-8bd7-4622-9d32-d382f3674485
md"""
### Comparisons
"""

# ╔═╡ 0d0af2ac-a896-4f70-bc8f-8715342037f0
xdf₃ = DataFrame(rand(2,3), :auto)

# ╔═╡ f2715e53-7a02-4eb8-ac0e-28902540e3b9
xdf₄ = copy(xdf₃)

# ╔═╡ 2471303f-8618-4c09-a511-2decf60a1c6c
@assert xdf₃ == xdf₄   ## compares column names and contents

# ╔═╡ 1aa0c67e-c1b0-4f0d-958f-84562eb39d76
md"""
Create a minimally different DataFrame and use `isapprox` for comparison:
"""

# ╔═╡ bef79b4f-d498-4e4f-b635-d769d0e191fd
xdf₅ = xdf₃ .+ eps()

# ╔═╡ 14df0d21-dc58-4385-a29c-d339d54a0719
begin
	@assert xdf₅ ≠ xdf₃                          ## different
	@assert isapprox(xdf₅, xdf₅, atol=eps()/2)   ## approx. equal
end

# ╔═╡ 1bbbb4ca-5d4a-466a-8541-8cdf7fe6c92c
md"""
`missings` are handled as in `Julia Base`:
"""

# ╔═╡ 3dc2357e-0cc7-4259-95e3-bc9b41243ebe
xdfₘ = DataFrame(a=missing)

# ╔═╡ 2ea87bb2-2d60-40f3-8728-ed5460c6c0c8
xdfₘ == xdfₘ

# ╔═╡ 8d06d5d7-2fb9-4b6a-be3c-1f6fcc4fd7ba
begin
	@assert xdfₘ === xdfₘ
	isequal(xdfₘ, xdfₘ)
end

# ╔═╡ 56edf702-2972-482a-b371-5d9d5de38fcd
html"""

"""

# ╔═╡ Cell order:
# ╟─b9438ea0-6a42-11eb-3203-bdef1a59bb19
# ╠═dac66550-6a42-11eb-18fd-1dbb51ebb356
# ╟─be11851f-b9a5-4140-9d56-fa0e1ec896f9
# ╟─0c9061f5-92fc-4861-8920-acb7f16aea51
# ╠═e8266538-6a42-11eb-0919-f7a19fd021ba
# ╠═f28d6da0-6a42-11eb-1bfd-edda225bd80a
# ╠═02630c26-6a43-11eb-375b-d3bf7b722863
# ╠═09ef1732-6a43-11eb-3ce5-6904ca52514f
# ╠═1436e9b8-6a43-11eb-1dde-d97c58bf5218
# ╠═6ca3a755-9f7b-4b49-b15c-1989144e4e44
# ╠═20d87024-6a43-11eb-0a14-d5c00e2aafde
# ╠═2c9fb4c6-6a43-11eb-3e4c-5d12a291df4a
# ╠═3534cd7e-6a43-11eb-2eb7-0fa46c82b12d
# ╟─7f69dd67-3fa6-4c74-bf8c-4fecb087b74f
# ╠═3f6d9596-6a43-11eb-3469-5dcd7f0332e0
# ╠═546cdee8-6a43-11eb-38a2-35fe89ae2dde
# ╠═5bc30e06-6a43-11eb-1625-015e8cb5bce1
# ╟─cc4ffb73-ada8-4446-8914-9e1d5cf669d4
# ╠═c6786444-6a43-11eb-21b9-3d4a2722cdfd
# ╠═04e69443-e80d-45ea-9d59-fb5dd8bf11e5
# ╠═d105dcca-6a43-11eb-07fe-cd24a5efcf2a
# ╟─49729e83-ebb2-4514-bdaf-68db8847800a
# ╠═0ea18045-6d7e-4202-ab94-709ab6150201
# ╠═288c1879-77db-4246-b087-5f4e7ee600c5
# ╠═2140fccb-bbbb-4bae-b3e0-40a850ca40b9
# ╟─b9566a8c-748b-4696-8aed-d9add8a085d0
# ╠═0ea6ef6d-d908-4888-9fa3-4bd5edefeea6
# ╟─2aae7a6a-21db-45ee-bcdf-619ee2ea5240
# ╠═9b8d271c-5def-46d4-a38d-a98a63e4d62c
# ╟─2e1155ca-0cfc-4113-b2ce-2ff770e51e52
# ╠═af6641c1-d1eb-4ae0-9e31-63aed4793b14
# ╟─bd1f96af-21fb-4e9e-933b-c4b99a05d13a
# ╠═6f4ab986-81af-4fbd-93b2-546b3ad7c425
# ╠═da813272-f45a-4fcd-a4ff-539bb4c924c1
# ╟─17f8d99c-0bf0-4442-b234-d8d4955482e0
# ╠═3be6e5f8-777b-4f9a-9430-362880e842ad
# ╠═137dfe1f-0a3d-4d76-978a-70da3c178479
# ╟─5f890d28-43ef-40f3-b1ff-b32c9cd5773d
# ╠═20ae267b-c26b-49c5-9096-96190f868c18
# ╟─49daec8a-0263-4f69-a27b-dc0efdf89c71
# ╠═ed5ee6d9-ec3e-417c-8206-f62666d88c4f
# ╟─05d0a958-6071-41ba-8219-98f040ac90ae
# ╟─c0960534-f6a1-40f8-8b35-b19aa5b78ac2
# ╠═0bfe2061-cddb-4f9f-9b7e-33a3e703d084
# ╠═5055355b-a405-4d27-b2aa-9477bd5ee2c4
# ╠═2623b999-d8b1-4cba-9f2b-4d644c672354
# ╠═3c8b682a-5387-40b6-b3cd-d8beaa09383c
# ╠═b1d6b7e3-8c67-44b3-ae30-87302295e976
# ╠═61ac92e7-f181-4554-809a-dba4b786acec
# ╟─e5e723bd-aa8e-415c-b7c6-ffd49bbf1841
# ╠═fe37621d-c903-4055-bdb9-0dbc98f9d524
# ╠═8624b878-775b-4e37-a5c5-9b1cf38ef614
# ╠═9d40b847-5ebc-4656-aa03-f918ccee3974
# ╟─12a63294-6a44-11eb-0f29-7dcf33acface
# ╠═8c59bad5-6a7c-4923-b088-07a225c57a96
# ╟─91e94361-3b46-4e0f-b8c3-14ced0b2e71b
# ╠═48c59474-ffb8-4c92-96d1-39dddcd39076
# ╠═2180b2ba-3201-49a6-a0d5-5549c00c7b37
# ╟─4c3feedd-f3e0-4121-91b6-9b43b0840c9b
# ╠═dcf94057-a263-4b27-b086-db0c74953405
# ╠═56ba4b95-5a94-46be-ab6c-c80a40daccec
# ╟─f0c302ae-23c1-4019-80d9-d727f83b9815
# ╠═e0f1d7c5-039a-4be6-9b84-e7434faafad6
# ╟─73b293e6-4609-419f-afe6-5f1185ab58f7
# ╠═09290340-e74c-4f2c-9962-ccac5dd334a8
# ╟─5a53e3f7-1585-4d16-ad4a-1e88ec08ef83
# ╟─96f4c0ff-8bd7-4622-9d32-d382f3674485
# ╠═0d0af2ac-a896-4f70-bc8f-8715342037f0
# ╠═f2715e53-7a02-4eb8-ac0e-28902540e3b9
# ╠═2471303f-8618-4c09-a511-2decf60a1c6c
# ╟─1aa0c67e-c1b0-4f0d-958f-84562eb39d76
# ╠═bef79b4f-d498-4e4f-b635-d769d0e191fd
# ╠═14df0d21-dc58-4385-a29c-d339d54a0719
# ╟─1bbbb4ca-5d4a-466a-8541-8cdf7fe6c92c
# ╠═3dc2357e-0cc7-4259-95e3-bc9b41243ebe
# ╠═2ea87bb2-2d60-40f3-8728-ed5460c6c0c8
# ╠═8d06d5d7-2fb9-4b6a-be3c-1f6fcc4fd7ba
# ╟─56edf702-2972-482a-b371-5d9d5de38fcd

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
using NumbersFromText
using CSV

data = rand(Int, 32*1024*1024)
buf = IOBuffer()
foreach(x -> println(buf, x), data)
seek(buf, 0); @time a, = readarrays!(buf, zeros(Int, 0))
seek(buf, 0); @time a, = readarrays!(Val{true}, buf, zeros(Int, 0))


text:
module CoordinateSystems

export CoordinateSystem,
Polar,
Cartesian,
isPolar,
isCartesian

abstract type CoordinateSystem end
struct Polar<:CoordinateSystem end
struct Cartesian<:CoordinateSystem end

"""
    bool = isCartesian(coordinate_system)
"""
isCartesian(coordinate_system::Cartesian) = true
isCartesian(coordinate_system) = false

"""
    bool = isPolar(coordinate_system)
"""
isPolar(coordinate_system::Polar) = true
isPolar(coordinate_system) = false

end # module


text:
include("euler/euler.jl")
using .Primes: prime_numbers
using .Numbers: sum_proper_divisors

function compute()::Int
    n = 28_123
    primes = prime_numbers(trunc(Int, √n))
    abundant_sieve, not_expressible = repeat([false], n), repeat([true], n)
    result = 0
    for i ∈ 1:n
        abundant_sieve[i] = i < sum_proper_divisors(i, primes)
        for j ∈ 1:i ÷ 2
            if abundant_sieve[j] && abundant_sieve[i - j]
                not_expressible[i] = false
                break
            end
        end
        if not_expressible[i] result += i end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq

using LinearAlgebra
using Plots
gr()

# Create a test problem
function lorenz(u,p,t)
    x, y, z = u

    ẋ = 10.0*(y - x)
    ẏ = x*(28.0-z) - y
    ż = x*y - (8/3)*z
    return [ẋ, ẏ, ż]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
dt = 0.005
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob, Tsit5(), saveat = dt)

plot(sol,vars=(1,2,3))

# Differential data from equations
X = Array(sol)
DX = similar(X)
for (i, xi) in enumerate(eachcol(X))
    DX[:,i] = lorenz(xi, [], 0.0)
end

# Estimate differential data from state variables via the Savitzky-Golay filter
# Test on a single variable
windowSize, polyOrder = 9, 4
DX1_cropped, DX1_sg = savitzky_golay(X[1,:], windowSize, polyOrder, deriv=1, dt=dt)

# By default, the savitzky_golay function crops the borders, where the estimation
# is less accurate. Optionally, this can be turned off by passing the
# argument `crop = false` to savitzky_golay function.

# Check if the estimated derivatives are approximate to the "ground truth"
halfWindow = Int(ceil((windowSize+1)/2))
DX = DX[:,halfWindow+1:end-halfWindow]
isapprox(DX1_sg, DX[1,:], rtol=1e-2)

plot(DX1_sg, label = "Estimated with Savitzky-Golay filter")
plot!(DX[1,:],label="Ground truth")

# Now let's estimate the derivatives for all variables and use them to infer the equations
X_cropped, DX_sg = savitzky_golay(X, windowSize, polyOrder, deriv=1, dt=dt)

# Create a basis
@variables u[1:3]

# Lots of polynomials
polys = [u[1]^0]
for i ∈ 0:3
    for j ∈ 0:3
        for  k ∈ 0:3
            push!(polys, u[1]^i * u[2]^j * u[3]^k)
        end
    end
end

# And some other stuff
h = [1u[1];1u[2]; cos(u[1]); sin(u[1]); u[1]*u[2]; u[1]*sin(u[2]); u[2]*cos(u[2]); polys...]
basis = Basis(h, u)

# Get the reduced basis via the sparse regression
opt = STRRidge(0.1)
Ψ = SInDy(X_cropped, DX_sg, basis, maxiter = 100, opt = opt, normalize = false)
print(Ψ)
print_equations(Ψ)
print_equations(Ψ, show_parameter = true)

# Let's try adding some noise
using Random
seed = MersenneTwister(3)
X_noisy = X + 0.01*randn(seed,size(X))

X_noisy_cropped, DX_sg = savitzky_golay(X_noisy, windowSize, polyOrder, deriv=1, dt=dt)
X_noisy_cropped, X_smoothed = savitzky_golay(X_noisy, windowSize, polyOrder, deriv=0, dt=dt)

Ψ = SInDy(X_smoothed, DX_sg, basis, maxiter = 100, opt = opt, normalize = false)
print_equations(Ψ)


text:
"""
    symmetries
A Dictionary that contains the function to apply for each symmetry to avoid searching a whole space.
"""
const usual_symmetries = Dict(
    :permutable => sort,
)

"""
    Constraint
Parametric stucture with the following fields.
- `concept`: a Boolean function that, given an assignment `x`, outputs `true` if `x` satisfies the constraint, and `false` otherwise.
- `error`: a positive function that works as preferences over invalid assignements. Return `0.0` if the constraint is satisfied, and a strictly positive real otherwise.
"""
mutable struct Constraint{FConcept <: Function,FError <: Function}
    args_length::Union{Nothing,Int}
    concept::FConcept
    error::FError
    params_length::Union{Nothing,Int}
    symmetries::Set{Symbol}

    function Constraint(;
        args_length=nothing,
        concept=x -> true,
        error=(x; param=nothing, dom_size=0) -> Float64(!concept(x)),
        param=0,
        syms=Set{Symbol}(),
    )
        new{typeof(concept),typeof(error)}(args_length, concept, error, param, syms)
    end
end

"""
    concept(c::Constraint)
Return the concept (function) of constraint `c`.
    concept(c::Constraint, x...; param = nothing)
Apply the concept of `c` to values `x` and optionally `param`.
"""
concept(c::Constraint) = c.concept
function concept(c::Constraint, x; param=nothing)
    return isnothing(param) ? concept(c)(x) : concept(c)(x; param)
end

"""
    error_f(c::Constraint)
Return the error function of constraint `c`.
    error_f(c::Constraint, x; param = nothing)
Apply the error function of `c` to values `x` and optionally `param`.
"""
error_f(c::Constraint) = c.error
function error_f(c::Constraint, x; param=nothing, dom_size=0)
    return isnothing(param) ? error_f(c)(x; dom_size) : error_f(c)(x; param, dom_size)
end

"""
    args_length(c::Constraint)
Return the expected length restriction of the arguments in a constraint `c`. The value `nothing` indicates that any strictly positive number of value is accepted.
"""
args_length(c::Constraint) = c.args_length

"""
    params_length(c::Constraint)
Return the expected length restriction of the arguments in a constraint `c`. The value `nothing` indicates that any strictly positive number of parameters is accepted.
"""
params_length(c::Constraint) = c.params_length

"""
    symmetries(c::Constraint)
Return the list of symmetries of `c`.
"""
symmetries(c::Constraint) = c.symmetries

function make_error(symb::Symbol)
    return begin
        isdefined(Constraints, Symbol("icn_$symb")) ? eval(Symbol("icn_$symb")) :
        isdefined(Constraints, Symbol("error_$symb")) ? eval(Symbol("error_$symb")) :
        ((x; param=nothing, dom_size=0) -> Float64(!eval(Symbol("concept_$symb"))(x)))
    end
end


text:
# GPU implementation of generator update of ADMM.
@kernel function generator_kernel(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u, v,
    l, rho,
    pgmin, pgmax,
    qgmin, qgmax,
    c2, c1
)
    I_ = @index(Group, Linear)
    J_ = @index(Local, Linear)
    I = J_ + (@groupsize()[1] * (I_ - 1))
    if (I <= ngen)
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1
        u[pg_idx] = max(pgmin[I],
                        min(pgmax[I],
                            (-(c1[I]*baseMVA + l[pg_idx] -
                            rho[pg_idx]*v[pg_idx])) / (2*c2[I]*(baseMVA^2) + rho[pg_idx])))
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] - rho[qg_idx]*v[qg_idx])) / rho[qg_idx]))
    end
end

function generator_kernel(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u, v, l, rho,
    pgmin::Array{Float64,1}, pgmax::Array{Float64,1},
    qgmin::Array{Float64,1}, qgmax::Array{Float64,1},
    c2::Array{Float64,1}, c1::Array{Float64,1}
)
    for I=1:ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1
        u[pg_idx] = max(pgmin[I],
                        min(pgmax[I],
                            (-(c1[I]*baseMVA + l[pg_idx] -
                               rho[pg_idx]*v[pg_idx])) / (2*c2[I]*(baseMVA^2) + rho[pg_idx])))
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] - rho[qg_idx]*v[qg_idx])) / rho[qg_idx]))
    end

    return
end

function generator_kernel(
    gen_mod::GeneratorModel,
    baseMVA::Float64, u, v, l, rho,
    device::KA.GPU
)
    nblk = div(gen_mod.ngen, 32, RoundUp)
    nblk
    wait(generator_kernel(device, 32, gen_mod.ngen)(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                u, v, l, rho, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1,
                dependencies=Event(device)
                )
    )
    return 0.0
end

function generator_kernel(
    gen_mod::GeneratorModel,
    baseMVA::Float64, u, v, l, rho,
    device::KA.CPU
)
    tcpu = @timed generator_kernel(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                u, v, l, rho, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1)
    return tcpu
end


text:
module TravelingSalesmanHeuristics

using Random, LinearAlgebra

include("helpers.jl")
include("simulated_annealing.jl")
include("lowerbounds.jl")

export solve_tsp, lowerbound, repetitive_heuristic, two_opt,
       nearest_neighbor, cheapest_insertion, simulated_annealing, farthest_insertion
"""
    solve_tsp(distmat; quality_factor = 40)

Approximately solve a TSP by specifying a distance matrix. Return a tuple `(path, cost)`.

The optional keyword `quality_factor` (real number in [0,100]; defaults to 40) specifies the
tradeoff between computation time and quality of solution returned. Higher values
tend to lead to better solutions found at the cost of more computation time.

!!! note
    It is not guaranteed that a call with a high `quality_factor` will always run slower or
    return a better solution than a call with a lower `quality_factor`, though this is
    typically the case. A `quality_factor` of 100 neither guarantees an optimal solution
    nor the best solution  that can be found via extensive use of the methods in this package.

See also: individual heuristic methods such as [`farthest_insertion`](@ref)
"""
function solve_tsp(distmat::AbstractMatrix{T}; quality_factor::Real = 40.0) where {T<:Real}
    if quality_factor < 0 || quality_factor > 100
        @warn "quality_factor keyword passed to solve_tsp must be in [0,100]"
        quality_factor = clamp(quality_factor, 0, 100)
    end

    lowest_threshold = 5

    # begin adding heuristics as dictated by the quality_factor,
    # starting with the very fastest
    answers = Vector{Tuple{Vector{Int}, T}}()
    push!(answers, farthest_insertion(distmat; do2opt = false))
    if quality_factor < lowest_threshold # fastest heuristic and out
        return answers[1]
    end
    # otherwise, we'll try several heuristics and return the best

    # add any nearest neighbor heuristics as dictated by quality_factor
    if quality_factor >= 60 # repetitive-NN, 2-opt on each iter
        push!(answers, repetitive_heuristic(distmat, nearest_neighbor; do2opt = true))
    elseif quality_factor >= 25 # repetitive-NN, 2-opt on final
        rnnpath, _ = repetitive_heuristic(distmat, nearest_neighbor; do2opt = false)
        push!(answers, two_opt(distmat, rnnpath))
    elseif quality_factor >= 15 # NN w/ 2-opt
        push!(answers, nearest_neighbor(distmat; do2opt = true))
    end

    # farthest insertions as needed
    if quality_factor >= 70 # repetitive-FI, 2-opt on each
        push!(answers, repetitive_heuristic(distmat, farthest_insertion; do2opt = true))
    elseif quality_factor >= 5 # FI w/ 2-opt
        push!(answers, farthest_insertion(distmat; do2opt = true))
    end

    # cheapest insertions
    if quality_factor >= 90 # repetitive-CI w/ 2-opt on each
        push!(answers, repetitive_heuristic(distmat, cheapest_insertion; do2opt = true))
    elseif quality_factor >= 35
        push!(answers, cheapest_insertion(distmat; do2opt = true))
    end

    # simulated annealing refinement, seeded with best so far
    if quality_factor >= 80
        _, bestind = findmin([pc[2] for pc in answers])
        bestpath, bestcost = answers[bestind]
        nstart = ceil(Int, (quality_factor - 79)/5)
        push!(answers,
              simulated_annealing(distmat; num_starts = nstart, init_path = bestpath)
        )
    end

    # pick best
    _, bestind = findmin([pc[2] for pc in answers])
    return answers[bestind]
end

###
# path generation heuristics
###

"""
    nearest_neighbor(distmat)

Approximately solve a TSP using the nearest neighbor heuristic. `distmat` is a square real matrix
where distmat[i,j] represents the distance from city `i` to city `j`. The matrix needn't be
symmetric and can contain negative values. Return a tuple `(path, pathcost)`.

# Optional keyword arguments:
- `firstcity::Union{Int, Nothing}`: specifies the city to begin the path on. Passing `nothing` or
    not specifying a value corresponds to random selection.
- `closepath::Bool = true`: whether to include the arc from the last city visited back to the
    first city in cost calculations. If true, the first city appears first and last in the path.
- `do2opt::Bool = true`: whether to refine the path found by 2-opt switches (corresponds to
    removing path crossings in the planar Euclidean case).
"""
function nearest_neighbor(distmat::AbstractMatrix{T} where {T<:Real};
                          firstcity::Union{Int, Nothing} = rand(1:size(distmat, 1)),
                          repetitive::Bool = false,
                          closepath::Bool = true,
                          do2opt::Bool = true)
    # must have a square matrix
    num_cities = check_square(distmat, "Must pass a square distance matrix to nearest_neighbor")

    # extract a safe int value for firstcity
    if firstcity == nothing # random first city
        firstcityint = rand(1:num_cities)
    else # an int was passed
        firstcityint = firstcity
        if !(1 <= firstcityint <= num_cities)
            error("firstcity of $(firstcity) passed to `nearest_neighbor`" *
                  " out of range [1, $(num_cities)]")
        end
    end

    # calling with KW repetitive is deprecated; pass the call to repetitive_heuristic
    if repetitive
        @warn ("Calling `nearest_neighbor` with keyword `repetitive` is deprecated;'" *
             " instead call `repetitive_heuristic(distmat, nearest_neighbor; kwargs...)`")
        return repetitive_heuristic(distmat, nearest_neighbor;
                                    closepath = closepath, do2opt = do2opt)
    end

    # put first city on path
    path = Vector{Int}()
    push!(path, firstcityint)

    # cities to visit
    citiesToVisit = collect(1:(firstcityint - 1))
    append!(citiesToVisit, collect((firstcityint + 1):num_cities))

    # nearest neighbor loop
    while !isempty(citiesToVisit)
        curCity = path[end]
        dists = distmat[curCity, citiesToVisit]
        _, nextInd = findmin(dists)
        nextCity = citiesToVisit[nextInd]
        push!(path, nextCity)
        deleteat!(citiesToVisit, nextInd)
    end

    # complete cycle? (duplicates first city)
    if closepath
        push!(path, firstcityint)
    end

    # do swaps?
    if do2opt
        path, _ = two_opt(distmat, path)
    end

    return (path, pathcost(distmat, path))
end

"""
    cheapest_insertion(distmat::Matrix, initpath::AbstractArray{Int})

Given a distance matrix and an initial path, complete the tour by repeatedly doing the cheapest
insertion. Return a tuple `(path, cost)`. The initial path must have length at least 2, but can
be simply `[i, i]` for some city index `i` which corresponds to starting with a loop at city `i`.

!!! note
    Insertions are always in the interior of the current path so this heuristic can also be used for
    non-closed TSP paths.

Currently the implementation is a naive ``n^3`` algorithm.
"""
function cheapest_insertion(distmat::AbstractMatrix{T}, initpath::AbstractVector{S}) where {T<:Real, S<:Integer}
    check_square(distmat, "Distance matrix passed to cheapest_insertion must be square.")

    n = size(distmat, 1)
    path = Vector{Int}(initpath)

    # collect cities to visited
    visitus = setdiff(collect(1:n), initpath)

    while !isempty(visitus)
        bestCost = Inf
        bestInsertion = (-1, -1)
        for k in visitus
            for after in 1:(length(path) - 1) # can't insert after end of path
                c = inscost(k, after, path, distmat)
                if c < bestCost
                    bestCost = c
                    bestInsertion = (k, after)
                end
            end
        end
        # bestInsertion now holds (k, after)
        # insert into path, remove from to-do list
        k, after = bestInsertion
        insert!(path, after + 1, k)
        visitus = setdiff(visitus, k)
    end

    return (path, pathcost(distmat, path))
end
"""
    cheapest_insertion(distmat; ...)

Complete a tour using cheapest insertion with a single-city loop as the initial path. Return a
tuple `(path, cost)`.

### Optional keyword arguments:
- `firstcity::Union{Int, Nothing}`: specifies the city to begin the path on. Passing `nothing` or
    not specifying a value corresponds to random selection.
- `do2opt::Bool = true`: whether to improve the path found by 2-opt swaps.
"""
function cheapest_insertion(distmat::AbstractMatrix{T} where{T<:Real};
                            firstcity::Union{Int, Nothing} = rand(1:size(distmat, 1)),
                            repetitive::Bool = false,
                            do2opt::Bool = true)
    #infer size
    num_cities = size(distmat, 1)

    # calling with KW repetitive is deprecated; pass the call to repetitive_heuristic
    if repetitive
        @warn "Calling `cheapest_insertionr` with keyword `repetitive` is deprecated;'" *
             " instead call `repetitive_heuristic(distmat, cheapest_insertion; kwargs...)`"
        return repetitive_heuristic(distmat, cheapest_insertion; do2opt = do2opt)
    end

    # extract a safe int value for firstcity
    if firstcity == nothing # random first city
        firstcityint = rand(1:num_cities)
    else # an int was passed
        firstcityint = firstcity
        if !(1 <= firstcityint <= num_cities)
            error("firstcity of $(firstcity) passed to `cheapest_insertion`" *
                  " out of range [1, $(num_cities)]")
        end
    end

    # okay, we're not repetitive, put a loop on the first city
    path, cost = cheapest_insertion(distmat, [firstcityint, firstcityint])

    # user may have asked for 2-opt refinement
    if do2opt
        path, cost = two_opt(distmat, path)
    end

    return path, cost
end

"""
    farthest_insertion(distmat; ...)

Generate a TSP path using the farthest insertion strategy. `distmat` must be a square real matrix.
Return a tuple `(path, cost)`.

### Optional arguments:
- `firstCity::Int`: specifies the city to begin the path on. Not specifying a value corresponds
    to random selection.
- `do2opt::Bool = true`: whether to improve the path by 2-opt swaps.
"""
function farthest_insertion(distmat::AbstractMatrix{T};
                            firstcity::Int = rand(1:size(distmat, 1)),
                            do2opt::Bool = true) where {T<:Real}
    n = check_square(distmat, "Must pass square distance matrix to farthest_insertion.")
    if firstcity < 1 || firstcity > n
        error("First city for farthest_insertion must be in [1,...,n]")
    end
    smallval = minimum(distmat) - one(T) # will never be the max

    path = Int[firstcity, firstcity]
    sizehint!(path, n + 1)
    dists_to_tour = (distmat[firstcity, :] + distmat[:, firstcity]) / 2
    dists_to_tour[firstcity] = smallval

    while length(path) < n + 1
        # city farthest from tour
        _, nextcity = findmax(dists_to_tour)

        # find where to add it
        bestcost = inscost(nextcity, 1, path, distmat)
        bestind = 1
        for ind in 2:(length(path) - 1)
            altcost = inscost(nextcity, ind, path, distmat)
            if altcost < bestcost
                bestcost = altcost
                bestind = ind
            end
        end
        # and add the next city at the best spot found
        insert!(path, bestind + 1, nextcity) # +1 since arg to insert! is where nextcity ends up

        # update distances to tour
        dists_to_tour[nextcity] = smallval
        for i in 1:n
            c = dists_to_tour[i]
            if c == zero(T) # i already in tour
                continue
            end
            altc = (distmat[i, nextcity] + distmat[nextcity, i]) / 2 # cost i--nextcity
            if altc < c
                @inbounds dists_to_tour[i] = altc
            end
        end
    end

    if do2opt
        path, _ = two_opt(distmat, path)
    end

    return path, pathcost(distmat, path)
end

###
# path improvement heuristics
###

"""
    two_opt(distmat, path)

Improve `path` by doing 2-opt switches (i.e. reversing part of the path) until doing so no
longer reduces the cost. Return a tuple `(improved_path, improved_cost)`.

On large problem instances this heuristic can be slow, but it is highly recommended on small and
medium problem instances.

See also [`simulated_annealing`](@ref) for another path generation heuristic.
"""
function two_opt(distmat::AbstractMatrix{T}, path::AbstractVector{S}) where {T<:Real, S<:Integer}
    # size checks
    n = length(path)
    if size(distmat, 1) != size(distmat, 2)
        error("Distance matrix passed to two_opt must be square.")
    end

    # don't modify input
    path = copy(path) # Int is a bitstype

    # main loop
    # check every possible switch until no 2-swaps reduce objective
    # if the path passed in is a loop (first/last nodes are the same)
    # then we must keep these the endpoints of the path the same
    # ie just keep it a loop, and therefore it doesn't matter which node is at the end
    # if the path is not a cycle, we should respect the endpoints
    switchLow = 2
    switchHigh = n - 1
    prevCost = Inf
    curCost = pathcost(distmat, path)
    while prevCost > pathcost(distmat, path)
        prevCost = curCost
        # we can't change the first
        for i in switchLow:(switchHigh-1)
            for j in switchHigh:-1:(i+1)
                altCost = pathcost_rev(distmat, path, i, j)
                if altCost < curCost
                    curCost = altCost
                    reverse!(path, i, j)
                end
            end
        end
    end
    return path, pathcost(distmat, path)
end

end # module


text:
strip_line_numbers!(arg) = arg
function strip_line_numbers!(ex::Expr)
    keep = trues(length(ex.args))
    for i = 1:length(ex.args)
        if isa(ex.args[i], LineNumberNode)
            keep[i] = false
        end
        isa(ex.args[i], Expr) || continue
        if ex.args[i].head == :line
            keep[i] = false
        else
            ex.args[i] = strip_line_numbers!(ex.args[i])
        end
    end
    ex.args = ex.args[keep]
    ex
end


text:
"""
    MarkovModulatedPoissonProcess{M,Tr<:ContinuousMarkovChain,Em<:TemporalPointProcess{M}}

Markov-Modulated Poisson Process with mark type `M`.

# Fields
- `transitions::Tr`: state evolution process.
- `emissions::Vector{Em}`: one emission distribution per state.
"""
struct MarkovModulatedPoissonProcess{
    M,
    Tr<:ContinuousMarkovChain,
    Em<:TemporalPointProcess{M},
} <: AbstractMeasure
    transitions::Tr
    emissions::Vector{Em}
end

const MMPP = MarkovModulatedPoissonProcess

## Access

transitions(mmpp::MMPP) = mmpp.transitions
initial_distribution(mmpp::MMPP) = initial_distribution(transitions(mmpp))
transition_matrix(mmpp::MMPP) = transition_matrix(transitions(mmpp))
rate_matrix(mmpp::MMPP) = rate_matrix(transitions(mmpp))

emissions(mmpp::MMPP) = mmpp.emissions
emission(mmpp::MMPP, s::Int) = mmpp.emissions[s]
nb_states(mmpp::MMPP) = length(emissions(mmpp))

## Simulation

function Base.rand(rng::AbstractRNG, mmpp::MMPP{M}, tmin, tmax) where {M}
    state_history = rand(rng, transitions(mmpp), tmin, tmax)
    transition_times, states = event_times(state_history), event_marks(state_history)
    observations = TemporalHistory(Float64[], M[], tmin, tmax)
    for k = 1:length(transition_times)
        local_s = states[k]
        local_tmin = transition_times[k]
        local_tmax = k < length(transition_times) ? transition_times[k+1] : tmax
        local_observations = rand(rng, emission(mmpp, local_s), local_tmin, local_tmax)
        append!(observations, local_observations)
    end
    return state_history, observations
end

Base.rand(mmpp::MMPP, args...) = rand(GLOBAL_RNG, mmpp, args...)


text:
@testset "one_component_3d" begin

	for T in (Int, Bool, Gray{Bool}, Gray{N0f8}, Gray{N0f16}, Gray{N0f32}, Gray{Float64})
		for i = 1:20
			test_volume = T.(eval(Symbol("test_volume_$(i)")))

			for n = 3:2:max(size(test_volume)...)
				# Call one component 3D algorithm and compare the number of components.
				labels = ImageComponentAnalysis.label_components(OneComponent3D(), test_volume, (n, n, n))
				num_components = maximum(labels)

				@test num_components == test_volume_results[i]
			end
		end
	end

end


text:
# Particularly interesting graph coloring instance derived from an AC power flow model.
using HDF5, JLD
using ReverseDiffSparse
#using ProfileView

d = load("benchmark2.jld");
I = d["I"];
J = d["J"];
n = d["n"];

g = gen_adjlist(I,J, n)
color, num_colors = acyclic_coloring(g)
@time color, num_colors = acyclic_coloring(g)
#@profile color, num_colors = acyclic_coloring(g)

rinfo = ReverseDiffSparse.recovery_preprocess(g, color, num_colors)
Profile.clear_malloc_data()
@time rinfo = ReverseDiffSparse.recovery_preprocess(g, color, num_colors)

#ProfileView.view()


text:
using Test
using NeXLSpectrum
using Random

between(x, a, b) = x >= min(a, b) && x <= max(a, b)

N = 100
mt = MersenneTwister(0xBADF00D)
ci1s = CartesianIndex.(zip(rand(mt, -1000:1000, N), rand(mt, -1000:1000, N)))
ci2s = CartesianIndex.(zip(rand(mt, -1000:1000, N), rand(mt, -1000:1000, N)))

for (ci1, ci2) in zip(ci1s, ci2s)
    # println(ci1, " to ", ci2)
    drawline(ci1, ci2, false) do pt
        @test between(pt[1], ci1[1], ci2[1])
        @test between(pt[2], ci1[2], ci2[2])
        t1 = (pt[1] - ci1[1]) / (ci2[1] - ci1[1])
        t2 = (pt[2] - ci1[2]) / (ci2[2] - ci1[2])
        @test abs(pt[1] - (ci1[1] + t2 * (ci2[1] - ci1[1]))) <= 0.5 || #
              abs(pt[2] - (ci1[2] + t1 * (ci2[2] - ci1[2]))) <= 0.5
    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
export maxima,
       minima

for (funcname, comp) in ((:maxima, :<),
                         (:minima, :>))

    @eval begin
        function ($funcname)(x::AbstractVector,
                        w::Integer=1,
                        includebounds::Bool=false)
            w > 0 || throw(ArgumentError("window cannot be negative"))
            xlen = length(x)
            idxs = Int[]

            # There can't be more than one peak every `w` elements
            maxN = includebounds ? cld(xlen,2*w) : fld(xlen,2*w)
            N = 0

            ww = [ collect(w:-1:1); collect(-w:-1) ]

            sizehint!(idxs,maxN)

            i = 1

            if includebounds
                @inbounds while i <= w
                    peak = true
                    for j in max(i-w, 1):(i+w)
                        if ($comp)(x[i], x[j])
                            peak &= false
                            break # No reason to continue checking if the rest of the elements qualify
                        end
                    end

                    if peak
                        push!(idxs,i)
                        N += 1
                        i += w # There can't be another peak for at least `w` more elements
                    end
                    i += 1
                end
            end

            i = 1+w
            maxI = xlen-w

            @inbounds while i <= maxI
                peak = true
                for j in ww # For all elements within the window
                    if ($comp)(x[i], x[i+j])
                        peak &= false
                        break # No reason to continue checking if the rest of the elements qualify
                    end
                end

                if peak
                    push!(idxs,i)
                    N += 1
                    i += w # There can't be another peak for at least `w` more elements
                end
                i += 1
            end

            if includebounds
                @inbounds while i <= xlen
                    peak = true
                    for j in (i-w):min((i+w),xlen)
                        if ($comp)(x[i], x[j])
                            peak &= false
                            break # No reason to continue checking if the rest of the elements qualify
                        end
                    end

                    if peak
                        push!(idxs,i)
                        N += 1
                        i += w # There can't be another peak for at least `w` more elements
                    end
                    i += 1
                end
            end

            resize!(idxs,N)
            sizehint!(idxs,N)

            return idxs
        end
    end
end

@doc """
    maxima(x[, w=1, includebounds=false])

Find the local maxima of `x`.

`w` sets the minimum allowed distance between maxima. If `includebounds` is `true`, maxima
are allowed to be less than `w` away from the bounds of `x`, otherwise, maxima may not be
located any closer to the ends of `x` than `w`.
"""
maxima

@doc """
    minima(x[, w=1, includebounds=false])

Find the local minima of `x`.

`w` sets the minimum allowed distance between minima. If `includebounds` is `true`, minima
are allowed to be less than `w` away from the bounds of `x`, otherwise, minima may not be
located any closer to the ends of `x` than `w`.
"""
minima


text:
# Reduction along specific dimensions
#
# This is going to be obsolete, the purpose of this remains here is for supporting foldl along dims
#

## auxiliary functions

import ArrayViews: parent, offset
offset_view(a::Number, ::Int, ::Int, ::Int) = a
offset_view(a::ContiguousArray, o::Int, m::Int) = ContiguousView(parent(a), offset(a) + o, (m,))
offset_view(a::ContiguousArray, o::Int, m::Int, n::Int) = ContiguousView(parent(a), offset(a) + o, (m, n))

#################################################
#
#    codegen facilities
#
#################################################

# reduction type dependent codes

type FoldlReduc end

extra_params(::FoldlReduc) = [:(op::Functor{2}), :(s::Number)]
extra_args(::FoldlReduc) = [:op, :s]

update_code(R::FoldlReduc, s, x) = :( @inbounds $s = evaluate(op, $s, $x) )

function emptyreduc_code(R::FoldlReduc, dst::Symbol, T::Symbol, n::Symbol)
    quote
        for i = 1 : $n
            @inbounds ($dst)[i] = s
        end
    end
end

reduce_result(R::FoldlReduc, ty) = ty

# core skeleton

function generate_reducedim_codes(AN::Int, accum::Symbol, reducty)

    # function names
    _accum_eachcol! = symbol("_$(accum)_eachcol!")
    _accum_eachrow! = symbol("_$(accum)_eachrow!")
    _accum = symbol("_$(accum)")
    _accum! = symbol("_$(accum)!")
    accum! = symbol("$(accum)!")

    # code preparation
    h = codegen_helper(AN)
    aparams = h.contiguous_aparams
    exparams = extra_params(reducty)
    exargs = extra_args(reducty)

    if AN == 0
        offset_args = [:(offset_view(a, ao, m, n))]
    else
        offset_args = [:fun, [:(offset_view($a, ao, m, n)) for a in h.args[2:end]]...]
    end

    quote
        global $(_accum_eachcol!)
        function $(_accum_eachcol!){R<:Number}(m::Int, n::Int, r::ContiguousArray{R}, $(exparams...), $(aparams...))
            offset = 0
            if m > 0
                for j = 1 : n
                    rj = ($_accum)(offset+1, offset+m, $(exargs...), $(h.args...))
                    @inbounds r[j] = rj
                    offset += m
                end
            else
                $(emptyreduc_code(reducty, :r, :R, :n))
            end
        end

        global $(_accum_eachrow!)
        function $(_accum_eachrow!){R<:Number}(m::Int, n::Int, r::ContiguousArray{R}, $(exparams...), $(aparams...))
            if n > 0
                for i = 1 : m
                    @inbounds vi = $(h.term(:i))
                    @inbounds r[i] = vi
                end

                offset = m
                for j = 2 : n
                    for i = 1 : m
                        idx = offset + i
                        @inbounds vi = $(h.term(:idx))
                        $(update_code(reducty, :(r[i]), :vi))
                    end
                    offset += m
                end
            else
                $(emptyreduc_code(reducty, :r, :R, :m))
            end
        end

        global $(_accum!)
        function $(_accum!)(r::ContiguousArray, $(exparams...), $(aparams...), dim::Int)
            shp = $(h.inputsize)

            if dim == 1
                m = shp[1]
                n = succ_length(shp, 1)
                $(_accum_eachcol!)(m, n, r, $(exargs...), $(h.args...))

            else
                m = prec_length(shp, dim)
                n = shp[dim]
                k = succ_length(shp, dim)

                if k == 1
                    $(_accum_eachrow!)(m, n, r, $(exargs...), $(h.args...))
                else
                    mn = m * n
                    ro = 0
                    ao = 0
                    for l = 1 : k
                        $(_accum_eachrow!)(m, n, offset_view(r, ro, m), $(exargs...), $(offset_args...))
                        ro += m
                        ao += mn
                    end
                end
            end
            return r
        end

        global $(accum!)
        function $(accum!)(r::ContiguousArray, $(exparams...), $(aparams...), dim::Int)
            rshp = Base.reduced_dims($(h.inputsize), dim)
            length(r) == prod(rshp) || error("Invalid argument dimensions.")
            $(_accum!)(r, $(exargs...), $(h.args...), dim)
        end

        global $(accum)
        function $(accum)($(exparams...), $(aparams...), dim::Int)
            rshp = Base.reduced_dims($(h.inputsize), dim)
            $(_accum!)(Array($(h.termtype), rshp), $(exargs...), $(h.args...), dim)
        end
    end
end

macro code_reducedim(AN, fname, reducty)
    R = eval(reducty)
    esc(generate_reducedim_codes(AN, fname, R()))
end

#################################################
#
#   folding along dims
#
#################################################

macro code_foldldim(AN, fname)
    esc(generate_foldldim_codes(AN, fname))
end

@code_reducedim 0 foldl FoldlReduc
@code_reducedim 1 foldl FoldlReduc
@code_reducedim 2 foldl FoldlReduc
@code_reducedim 3 foldl FoldlReduc
@code_reducedim (-2) foldl_fdiff FoldlReduc


text:
module AffineArithmetic

using IntervalArithmetic

import Base: +, -, *, /, ^, ==,
            zero, one, range,
            show

include("aff.jl")
include("affine.jl")
include("full_affine.jl")

export Affine, affine
export FullAffine, reset_affine_index

end

max_stars_repo_name: dpsanders/AffineArithmetic

max_issues_repo_name: dpsanders/AffineArithmetic

max_forks_repo_name: dpsanders/AffineArithmetic