text:
for PT in (:Particles, :StaticParticles)
    @eval begin
        Base.promote_rule(::Type{Complex{S}}, ::Type{$PT{T,N}}) where {S<:Real,T<:Real,N} = Complex{$PT{promote_type(S,T),N}}

    end

    for ff in (^,)
        f = nameof(ff)
        @eval Base.$f(z::Complex{$PT{T,N}}, x::Real) where {T,N} = ℂ2ℂ_function($f, z, x)
        @eval Base.$f(z::Real, x::Complex{$PT{T,N}}) where {T,N} = ℂ2ℂ_function($f, z, x)
        @eval Base.$f(z::Complex{$PT{T,N}}, x::Complex{$PT{T,N}}) where {T,N} = ℂ2ℂ_function($f, z, x)

        @eval Base.$f(z::Complex{$PT{T,N}}, x::Int) where {T,N} = ℂ2ℂ_function($f, z, x)
        @eval Base.$f(z::Int, x::Complex{$PT{T,N}}) where {T,N} = ℂ2ℂ_function($f, z, x)
    end

end

@inline maybe_complex_particles(x,i) = x
@inline maybe_complex_particles(p::AbstractParticles,i) = complex(real(p.particles[i]), imag(p.particles[i]))
@inline maybe_complex_particles(p::Complex{<:AbstractParticles},i) = complex(p.re.particles[i], p.im.particles[i])

"""
    ℂ2ℂ_function(f::Function, z::Complex{<:AbstractParticles})

Helper function for uncertainty propagation through complex-valued functions of complex arguments.
applies `f : ℂ → ℂ ` to `z::Complex{<:AbstractParticles}`.
"""
function ℂ2ℂ_function(f::F, z::Complex{T}) where {F<:Union{Function,DataType},T<:AbstractParticles}
    s = map(1:length(z.re.particles)) do i
        @inbounds f(maybe_complex_particles(z, i))
    end
    complex(T(real.(s)), T(imag.(s)))
end

function ℂ2ℂ_function(f::F, z::Union{Complex{T},T}, a::R) where {F<:Union{Function,DataType},T<:AbstractParticles,R<:Real}
    s = map(1:length(z.re.particles)) do i
        @inbounds f(maybe_complex_particles(z, i),  vecindex(a, i))
    end
    complex(T(real.(s)), T(imag.(s)))
end

function ℂ2ℂ_function(f::F, z::R, a::Complex{S}) where {F<:Union{Function,DataType},S<:AbstractParticles,R<:Real}
    s = map(1:length(a.re.particles)) do i
        @inbounds f(vecindex(z, i), maybe_complex_particles(a, i))
    end
    complex(S(real.(s)), S(imag.(s)))
end

function ℂ2ℂ_function(f::F, z::Complex{T}, a::Complex{S}) where {F<:Union{Function,DataType},T<:AbstractParticles,S<:AbstractParticles}
    s = map(1:length(z.re.particles)) do i
        @inbounds f(maybe_complex_particles(z, i),  maybe_complex_particles(a, i))
    end
    complex(T(real.(s)), T(imag.(s)))
end

# function ℂ2ℂ_function(f::F, z::Complex{T}, a::Complex{S}) where {F<:Union{Function,DataType},T<:AbstractParticles,S<:AbstractParticles}
#     out = deepcopy(z)
#     rp, ip = out.re.particles, out.im.particles
#     @inbounds for i = 1:length(z.re.particles)
#         res = f(maybe_complex_particles(z, i),  maybe_complex_particles(a, i))
#         rp[i] = res.re
#         ip[i] = res.im
#     end
#     out
# end

function ℂ2ℂ_function!(f::F, s, z::Complex{T}) where {F,T<:AbstractParticles}
    map!(s, 1:length(z.re.particles)) do i
        @inbounds f(maybe_complex_particles(z, i))
    end
    complex(T(real.(s)), T(imag.(s)))
end

for ff in (sqrt, exp, exp10, log, log10, sin, cos, tan)
    f = nameof(ff)
    @eval Base.$f(z::Complex{<: AbstractParticles}) = ℂ2ℂ_function($f, z)
    @eval $(Symbol(f,:!))(s, z::Complex{<: AbstractParticles}) = ℂ2ℂ_function!($f, s, z)
end

Base.isinf(p::Complex{<: AbstractParticles}) = isinf(real(p)) || isinf(imag(p))
Base.isfinite(p::Complex{<: AbstractParticles}) = isfinite(real(p)) && isfinite(imag(p))

function Base.:(/)(a::Union{T, Complex{T}}, b::Complex{T}) where T<:AbstractParticles
    ℂ2ℂ_function(/, a, b)
end

function Base.FastMath.div_fast(a::Union{T, Complex{T}}, b::Complex{T}) where T<:AbstractParticles
    ℂ2ℂ_function(Base.FastMath.div_fast, a, b)
end


text:
import SDtoolbox: ct, PostShock_fr, CJspeed, zndsolve, CANTERA_MECH_FILETYPE
using PyCall
using StatsBase
using Interpolations

P₁ = 1e5
T₁ = 300

ϕ = 1

X₁ = "H2:$(42*ϕ), O2:21,N2:79"
mech = "gri30."*CANTERA_MECH_FILETYPE

gas₁ = ct.Solution(mech)
gas₁.TPX = T₁,P₁,X₁

U₁ = CJspeed(P₁,T₁,X₁,mech)

gas = PostShock_fr(U₁, P₁, T₁, X₁, mech)

out = zndsolve(gas,gas₁,U₁,advanced_output=true)

# Tplot = plot(out["time"],out["T"])
# therm_plot = plot(out["time"],out["thermicity"])

# compare with original sdtoolbox
pyimport_conda("scipy","scipy")
pyimport_conda("matplotlib","matplotlib")
pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)
zndsolve_py = pyimport("sdtoolbox.znd").zndsolve
gas = pyimport("sdtoolbox.postshock").PostShock_fr(U₁,P₁,T₁,X₁,mech)
out_py = zndsolve_py(gas,gas₁,U₁,t_end=out["time"][end])
# plot!(Tplot, out_py["time"],out_py["T"])
# plot!(therm_plot, out_py["time"],out_py["thermicity"])

# calculate cross correlation to quantify error
itp_T = LinearInterpolation(out["time"][1:(end-1)], out["T"][1:(end-1)])
sol_itp_T = itp_T(out_py["time"])
itp_therm = LinearInterpolation(out["time"][1:(end-1)], out["thermicity"][1:(end-1)])
sol_itp_therm = itp_therm(out_py["time"])

err_T = maximum(crosscor(sol_itp_T,out_py["T"]))
err_therm = maximum(crosscor(sol_itp_therm,out_py["thermicity"]))

ext: jl


text:
using DataDeps

register(DataDep(
  "shuttle-landing-control",
  "http://archive.ics.uci.edu/ml/datasets/Shuttle+Landing+Control",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/shuttle-landing-control/shuttle-landing-control.data",
  "a93ffca47434936a33e227ce3f268eaab54794dc0303d08fb8fc3856c5e7174f",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:7,
        categoric_indices=2:7,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
struct CDomain <: AbstractComponent
    dim::Int
    CDomain(dim::Int=1) = new(dim)
end

prepare!(comp::CDomain, domain::AbstractDomain{1}) = domain[:]
prepare!(comp::CDomain, domain::AbstractDomain) = domain[comp.dim]

evaluate!(buffer, comp::CDomain, domain::AbstractDomain) = nothing


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# What is included in the file `FMI3_md.jl` (model description)?
# - [Sec. 1a] the function `fmi3LoadModelDescription` to load/parse the FMI model description [exported]
# - [Sec. 1b] helper functions for the load/parse function [not exported]
# - [Sec. 2]  functions to get values from the model description in the format `fmi3Get[value](md::fmi3ModelDescription)` [exported]
# - [Sec. 3]  additional functions to get useful information from the model description in the format `fmi3Get[value](md::fmi3ModelDescription)` [exported]

using EzXML

using FMICore: fmi3ModelDescriptionModelExchange, fmi3ModelDescriptionCoSimulation, fmi3ModelDescriptionDefaultExperiment
using FMICore: fmi3ModelDescriptionFloat, fmi3ModelDescriptionBoolean, fmi3ModelDescriptionInteger, fmi3ModelDescriptionString, fmi3ModelDescriptionEnumeration
using FMICore: fmi3ModelDescriptionModelStructure
using FMICore: fmi3DependencyKind
######################################
# [Sec. 1a] fmi3LoadModelDescription #
######################################

"""
Extract the FMU variables and meta data from the ModelDescription
"""
function fmi3LoadModelDescription(pathToModellDescription::String)
    md = fmi3ModelDescription()

    md.stringValueReferences = Dict{String, fmi3ValueReference}()
    md.outputValueReferences = Array{fmi3ValueReference}(undef, 0)
    md.inputValueReferences = Array{fmi3ValueReference}(undef, 0)
    md.stateValueReferences = Array{fmi3ValueReference}(undef, 0)
    md.derivativeValueReferences = Array{fmi3ValueReference}(undef, 0)
    md.intermediateUpdateValueReferences = Array{fmi3ValueReference}(undef, 0)
    md.numberOfEventIndicators = 0

    md.enumerations = []
    typedefinitions = nothing
    modelvariables = nothing
    modelstructure = nothing

    doc = readxml(pathToModellDescription)

    root = doc.root

    # mandatory
    md.fmiVersion = root["fmiVersion"]
    md.modelName = root["modelName"]
    md.instantiationToken = root["instantiationToken"]

    # optional
    md.generationTool = parseNodeString(root, "generationTool"; onfail="[Unknown generation tool]")
    md.generationDateAndTime = parseNodeString(root, "generationDateAndTime"; onfail="[Unknown generation date and time]")
    variableNamingConventionStr = parseNodeString(root, "variableNamingConvention"; onfail= "flat")
    @assert (variableNamingConventionStr == "flat" || variableNamingConventionStr == "structured") ["fmi3ReadModelDescription(...): Unknown entry for `variableNamingConvention=$(variableNamingConventionStr)`."]
    md.variableNamingConvention = (variableNamingConventionStr == "flat" ? fmi3VariableNamingConventionFlat : fmi3VariableNamingConventionStructured)
    md.description = parseNodeString(root, "description"; onfail="[Unknown Description]")

    # defaults
    md.modelExchange = nothing
    md.coSimulation = nothing
    md.scheduledExecution = nothing

    # additionals
    md.valueReferences = []
    md.valueReferenceIndicies = Dict{UInt, UInt}()

    md.defaultStartTime = 0.0
    md.defaultStopTime = 1.0
    md.defaultTolerance = 0.0001

    for node in eachelement(root)
        if node.name == "CoSimulation" || node.name == "ModelExchange" || node.name == "ScheduledExecution"
            if node.name == "CoSimulation"
                md.coSimulation = fmi3ModelDescriptionCoSimulation()
                md.coSimulation.modelIdentifier                        = node["modelIdentifier"]
                md.coSimulation.canHandleVariableCommunicationStepSize = parseNodeBoolean(node, "canHandleVariableCommunicationStepSize"   ; onfail=false)
                md.coSimulation.canInterpolateInputs                   = parseNodeBoolean(node, "canInterpolateInputs"                     ; onfail=false)
                md.coSimulation.maxOutputDerivativeOrder               = parseNodeInteger(node, "maxOutputDerivativeOrder"                 ; onfail=nothing)
                md.coSimulation.canGetAndSetFMUstate                   = parseNodeBoolean(node, "canGetAndSetFMUState"                     ; onfail=false)
                md.coSimulation.canSerializeFMUstate                   = parseNodeBoolean(node, "canSerializeFMUState"                     ; onfail=false)
                md.coSimulation.providesDirectionalDerivatives         = parseNodeBoolean(node, "providesDirectionalDerivatives"           ; onfail=false)
                md.coSimulation.providesAdjointDerivatives             = parseNodeBoolean(node, "providesAdjointDerivatives"               ; onfail=false)
                md.coSimulation.hasEventMode                           = parseNodeBoolean(node, "hasEventMode"                             ; onfail=false)
            end

            if node.name == "ModelExchange"
                md.modelExchange = fmi3ModelDescriptionModelExchange()
                md.modelExchange.modelIdentifier                        = node["modelIdentifier"]
                md.modelExchange.canGetAndSetFMUstate                   = parseNodeBoolean(node, "canGetAndSetFMUState"                     ; onfail=false)
                md.modelExchange.canSerializeFMUstate                   = parseNodeBoolean(node, "canSerializeFMUState"                     ; onfail=false)
                md.modelExchange.providesDirectionalDerivatives         = parseNodeBoolean(node, "providesDirectionalDerivatives"           ; onfail=false)
                md.modelExchange.providesAdjointDerivatives             = parseNodeBoolean(node, "providesAdjointDerivatives"               ; onfail=false)
            end

            if node.name == "ScheduledExecution"
                md.scheduledExecution = FMICore.fmi3ModelDescriptionScheduledExecution()
                md.scheduledExecution.modelIdentifier                        = node["modelIdentifier"]
                md.scheduledExecution.needsExecutionTool                     = parseNodeBoolean(node, "needsExecutionTool"                       ; onfail=false)
                md.scheduledExecution.canBeInstantiatedOnlyOncePerProcess    = parseNodeBoolean(node, "canBeInstantiatedOnlyOncePerProcess"      ; onfail=false)
                md.scheduledExecution.canGetAndSetFMUstate                   = parseNodeBoolean(node, "canGetAndSetFMUState"                     ; onfail=false)
                md.scheduledExecution.canSerializeFMUstate                   = parseNodeBoolean(node, "canSerializeFMUState"                     ; onfail=false)
                md.scheduledExecution.providesDirectionalDerivatives         = parseNodeBoolean(node, "providesDirectionalDerivatives"           ; onfail=false)
                md.scheduledExecution.providesAdjointDerivatives             = parseNodeBoolean(node, "providesAdjointDerivatives"               ; onfail=false)
                md.scheduledExecution.providesPerElementDependencies         = parseNodeBoolean(node, "providesPerElementDependencies"           ; onfail=false)
            end
        elseif node.name == "TypeDefinitions"
            md.enumerations = createEnum(node)

        elseif node.name == "ModelVariables"
            md.modelVariables = parseModelVariables(node, md)

        elseif node.name == "ModelStructure"
            md.modelStructure = fmi3ModelDescriptionModelStructure()

            parseModelStructure(node, md)
            # for element in eachelement(node)
            #     if element.name == "ContinuousStateDerivative"
            #         parseContinuousStateDerivatives(element, md)
            #     elseif element.name == "InitialUnknowns"
            #         parseInitialUnknowns(element, md)
            #     elseif element.name == "Outputs"
            #         parseOutputs(element, md)
            #     elseif element.name == "EventIndicator"
            #         parseEventIndicator(element, md)
            #     else
            #         @warn "Unknown tag `$(element.name)` for node `ModelStructure`."
            #     end
            # end
        elseif node.name == "DefaultExperiment"
            md.defaultExperiment = fmi3ModelDescriptionDefaultExperiment()
            md.defaultExperiment.startTime  = parseNodeReal(node, "startTime")
            md.defaultExperiment.stopTime   = parseNodeReal(node, "stopTime")
            md.defaultExperiment.tolerance  = parseNodeReal(node, "tolerance"; onfail = md.defaultTolerance)
            md.defaultExperiment.stepSize   = parseNodeReal(node, "stepSize")
        end
    end

    # creating an index for value references (fast look-up for dependencies)
    for i in 1:length(md.valueReferences)
        md.valueReferenceIndicies[md.valueReferences[i]] = i
    end

    # check all intermediateUpdate variables
    for variable in md.modelVariables
        if Bool(variable.datatype.intermediateUpdate)
            push!(md.intermediateUpdateValueReferences, variable.valueReference)
        end
    end

    md
end

#######################################
# [Sec. 1b] helpers for load function #
#######################################

# Returns the indices of the state derivatives.
function fmi3getDerivativeIndices(node::EzXML.Node)
    indices = []
    for element in eachelement(node)
        if element.name == "InitialUnknown"
            ind = parse(Int, element["valueReference"])
            der = nothing
            derKind = nothing

            if haskey(element, "dependencies")
                der = split(element["dependencies"], " ")

                if der[1] == ""
                    der = fmi3Int32[]
                else
                    der = collect(parse(fmi3Int32, e) for e in der)
                end
            end

            if haskey(element, "dependenciesKind")
                derKind = split(element["dependenciesKind"], " ")
            end

            push!(indices, (ind, der, derKind))
        end
    end
    sort!(indices, rev=true)
end

# Parses the model variables of the FMU model description.
function parseModelVariables(nodes::EzXML.Node, md::fmi3ModelDescription)
    numberOfVariables = 0
    for node in eachelement(nodes)
        numberOfVariables += 1
    end
    modelVariables = Array{fmi3ModelVariable}(undef, numberOfVariables)
    index = 1
    for node in eachelement(nodes)
        name = node["name"]
        valueReference = parse(fmi3ValueReference, (node["valueReference"]))

        modelVariables[index] = fmi3ModelVariable(name, valueReference)

        if !(valueReference in md.valueReferences)
            push!(md.valueReferences, valueReference)
        end

        if haskey(node, "description")
            modelVariables[index].description = node["description"]
        end
        if haskey(node, "causality")
            modelVariables[index].causality = fmi3StringToCausality(node["causality"])

            if modelVariables[index].causality == fmi3CausalityOutput
                push!(md.outputValueReferences, valueReference)
            elseif modelVariables[index].causality == fmi3CausalityInput
                push!(md.inputValueReferences, valueReference)
            end
        end
        if haskey(node, "variability")
            modelVariables[index].variability = fmi3StringToVariability(node["variability"])
        end

        if haskey(node, "initial")
            modelVariables[index].initial = fmi3StringToInitial(node["initial"])
        end
        modelVariables[index].datatype = fmi3SetDatatypeVariables(node, md)

        # type node
        typenode = nothing
        typename = node.name

        if typename == "Float32" || typename == "Float64"
            modelVariables[index]._Float = fmi3ModelDescriptionFloat()
            typenode = modelVariables[index]._Float
            if haskey(node, "quantity")
                typenode.quantity = node["quantity"]
            end
            if haskey(node, "unit")
                typenode.unit = node["unit"]
            end
            if haskey(node, "displayUnit")
                typenode.displayUnit = node["displayUnit"]
            end
            if haskey(node, "relativeQuantity")
                typenode.relativeQuantity = fmi3parseBoolean(node["relativeQuantity"])
            end
            if haskey(node, "min")
                typenode.min = fmi3parseFloat(node["min"])
            end
            if haskey(node, "max")
                typenode.max = fmi3parseFloat(node["max"])
            end
            if haskey(node, "nominal")
                typenode.nominal = fmi3parseFloat(node["nominal"])
            end
            if haskey(node, "unbounded")
                typenode.unbounded = fmi3parseBoolean(node["unbounded"])
            end
            if haskey(node, "start")
                typenode.start = fmi3parseFloat(node["start"])
            end
            if haskey(node, "derivative")
                typenode.derivative = parse(UInt, node["derivative"])
            end
        elseif typename == "String"
            modelVariables[index]._String = fmi3ModelDescriptionString()
            typenode = modelVariables[index]._String
            if haskey(node, "start")
                modelVariables[index]._String.start = node["start"]
            end
            # ToDo: remaining attributes
        elseif typename == "Boolean"
            modelVariables[index]._Boolean = fmi3ModelDescriptionBoolean()
            typenode = modelVariables[index]._Boolean
            if haskey(node, "start")
                modelVariables[index]._Boolean.start = fmi3parseBoolean(node["start"])
            end
            # ToDo: remaining attributes
        elseif typename == "Int32"
            modelVariables[index]._Integer = fmi3ModelDescriptionInteger()
            typenode = modelVariables[index]._Integer
            if haskey(node, "start")
                modelVariables[index]._Integer.start = fmi3parseInteger(node["start"])
            end
            # ToDo: remaining attributes
        elseif typename == "Enumeration"
            modelVariables[index]._Enumeration = fmi3ModelDescriptionEnumeration()
            typenode = modelVariables[index]._Enumeration
            # ToDo: Save start value
            # ToDo: remaining attributes
        else
            @warn "Unknown data type `$(typename)`."
        end

        md.stringValueReferences[name] = valueReference

        index += 1
    end
    md.numberOfContinuousStates = length(md.stateValueReferences)
    modelVariables
end

# Parses the model variables of the FMU model description.
function parseModelStructure(nodes::EzXML.Node, md::fmi3ModelDescription)
    @assert (nodes.name == "ModelStructure") "Wrong section name."
    md.modelStructure.continuousStateDerivatives = []
    md.modelStructure.initialUnknowns = []
    md.modelStructure.eventIndicators = []
    md.modelStructure.outputs = []
    for node in eachelement(nodes)
        if haskey(node, "valueReference")
            varDep = parseDependencies(node)
            if node.name == "InitialUnknown"
                push!(md.modelStructure.initialUnknowns, varDep)
            elseif node.name == "EventIndicator"
                md.numberOfEventIndicators += 1
                push!(md.modelStructure.eventIndicators)
                # TODO parse valueReferences to another array
            elseif node.name == "ContinuousStateDerivative"

                # find states and derivatives^
                derSV = fmi3ModelVariablesForValueReference(md, fmi3ValueReference(fmi3parseInteger(node["valueReference"])))[1]
                # derSV = md.modelVariables[varDep.index]
                derVR = derSV.valueReference
                stateVR = md.modelVariables[derSV._Float.derivative].valueReference

                if stateVR ∉ md.stateValueReferences
                    push!(md.stateValueReferences, stateVR)
                end
                if derVR ∉ md.derivativeValueReferences
                    push!(md.derivativeValueReferences, derVR)
                end

                push!(md.modelStructure.continuousStateDerivatives, varDep)
            elseif node.name =="Output"
                # find outputs
                outVR = fmi3ValueReference(fmi3parseInteger(node["valueReference"]))

                if outVR ∉ md.outputValueReferences
                    push!(md.outputValueReferences, outVR)
                end

                push!(md.modelStructure.outputs, varDep)
            else
                @warn "Unknown entry in `ModelStructure` named `$(node.name)`."
            end
        else
            @warn "Invalid entry for node `$(node.name)` in `ModelStructure`, missing entry `valueReference`."
        end
    end
end

function parseDependencies(node::EzXML.Node)
    varDep = fmi3VariableDependency(fmi3ValueReference(fmi3parseInteger(node["valueReference"])))

    if haskey(node, "dependencies")
        dependencies = node["dependencies"]
        if length(dependencies) > 0
            dependenciesSplit = split(dependencies, " ")
            if length(dependenciesSplit) > 0
                varDep.dependencies = collect(parse(UInt, e) for e in dependenciesSplit)
            end
        end
    end

    if haskey(node, "dependenciesKind")
        dependenciesKind = node["dependenciesKind"]
        if length(dependenciesKind) > 0
            dependenciesKindSplit = split(dependenciesKind, " ")
            if length(dependenciesKindSplit) > 0
                varDep.dependenciesKind = collect(fmi3StringToDependencyKind(e) for e in dependenciesKindSplit)
            end
        end
    end

    if varDep.dependencies != nothing && varDep.dependenciesKind != nothing
        if length(varDep.dependencies) != length(varDep.dependenciesKind)
            @warn "Length of field dependencies ($(length(varDep.dependencies))) doesn't match length of dependenciesKind ($(length(varDep.dependenciesKind)))."
        end
    end

    return varDep
end

function parseContinuousStateDerivative(nodes::EzXML.Node, md::fmi3ModelDescription)
    @assert (nodes.name == "ContinuousStateDerivative") "Wrong element name."
    md.modelStructure.derivatives = []
    for node in eachelement(nodes)
        if node.name == "InitialUnknown"
            if haskey(node, "index")
                varDep = parseUnknwon(node)

                # find states and derivatives
                derSV = md.modelVariables[varDep.index]
                derVR = derSV.valueReference
                stateVR = md.modelVariables[derSV._Float.derivative].valueReference

                if stateVR ∉ md.stateValueReferences
                    push!(md.stateValueReferences, stateVR)
                end
                if derVR ∉ md.derivativeValueReferences
                    push!(md.derivativeValueReferences, derVR)
                end

                push!(md.modelStructure.derivatives, varDep)
            else
                @warn "Invalid entry for node `Unknown` in `ModelStructure`, missing entry `index`."
            end
        elseif node.name == "EventIndicator"
            md.numberOfEventIndicators += 1
            # TODO parse valueReferences to another array
        else
            @warn "Unknown entry in `ModelStructure.Derivatives` named `$(node.name)`."
        end
    end
end

function parseInitialUnknowns(node::EzXML.Node, md::fmi3ModelDescription)
    @assert (node.name == "InitialUnknowns") "Wrong element name."
    for node in eachelement(nodes)
        if node.name == "Unknown"
            if haskey(node, "index")
                varDep = parseUnknwon(node)

                push!(md.modelStructure.initialUnknowns, varDep)
            else
                @warn "Invalid entry for node `Unknown` in `ModelStructure`, missing entry `index`."
            end
        else
            @warn "Unknown entry in `ModelStructure.InitialUnknowns` named `$(node.name)`."
        end
    end
end

function parseOutputs(nodes::EzXML.Node, md::fmi3ModelDescription)
    @assert (nodes.name == "Outputs") "Wrong element name."
    md.modelStructure.outputs = []
    for node in eachelement(nodes)
        if node.name == "Unknown"
            if haskey(node, "index")
                varDep = parseUnknwon(node)

                # find outputs
                outVR = md.modelVariables[varDep.index].valueReference

                if outVR ∉ md.outputValueReferences
                    push!(md.outputValueReferences, outVR)
                end

                push!(md.modelStructure.outputs, varDep)
            else
                @warn "Invalid entry for node `Unknown` in `ModelStructure`, missing entry `index`."
            end
        else
            @warn "Unknown entry in `ModelStructure.Outputs` named `$(node.name)`."
        end
    end
end

# Parses a real value represented by a string.
function fmi3parseFloat(s::Union{String, SubString{String}}; onfail=nothing)
    if onfail == nothing
        return parse(fmi3Float64, s)
    else
        try
            return parse(fmi3Float64, s)
        catch
            return onfail
        end
    end
end

function fmi3parseNodeFloat(node, key; onfail=nothing)
    if haskey(node, key)
        return fmi3parseFloat(node[key]; onfail=onfail)
    else
        return onfail
    end
end

# Parses a Bool value represented by a string.
function fmi3parseBoolean(s::Union{String, SubString{String}}; onfail=nothing)
    if s == "true"
        return true
    elseif s == "false"
        return false
    else
        @assert onfail !== nothing ["parseBoolean(...) unknown boolean value '$s'."]
        return onfail
    end
end

function fmi3parseNodeBoolean(node, key; onfail=nothing)
    if haskey(node, key)
        return fmi3parseBoolean(node[key]; onfail=onfail)
    else
        return onfail
    end
end

# Parses an Integer value represented by a string.
function fmi3parseInteger(s::Union{String, SubString{String}}; onfail=nothing)
    if onfail == nothing
        return parse(Int, s)
    else
        try
            return parse(Int, s)
        catch
            return onfail
        end
    end
end

function fmi3parseNodeInteger(node, key; onfail=nothing)
    if haskey(node, key)
        return fmi3parseInteger(node[key]; onfail=onfail)
    else
        return onfail
    end
end

function fmi3parseNodeString(node, key; onfail=nothing)
    if haskey(node, key)
        return node[key]
    else
        return onfail
    end
end

# Parses a fmi3Boolean value represented by a string.
function parseFMI3Boolean(s::Union{String, SubString{String}})
    if fmi3parseBoolean(s)
        return fmi3True
    else
        return fmi3False
    end
end

# set the datatype and attributes of an model variable
function fmi3SetDatatypeVariables(node::EzXML.Node, md::fmi3ModelDescription)
    type = fmi3DatatypeVariable()
    typename = node.name
    type.canHandleMultipleSet = nothing
    type.intermediateUpdate = fmi3False
    type.previous = nothing
    type.clocks = nothing
    type.declaredType = nothing
    type.start = nothing
    type.min = nothing
    type.max = nothing
    type.initial = nothing
    type.quantity = nothing
    type.unit = nothing
    type.displayUnit = nothing
    type.relativeQuantity = nothing
    type.nominal = nothing
    type.unbounded = nothing
    type.derivative = nothing
    type.reinit = nothing
    type.mimeType = nothing
    type.maxSize = nothing
    type.datatype = nothing

    # TODO fmi3Boolean, fmi3UInt8 are the same datatype so they get recognized the same
    if typename == "Float32"
        type.datatype = fmi3Float32
    elseif typename == "Float64"
        type.datatype = fmi3Float64
    elseif typename == "Int8"
        type.datatype = fmi3Int8
    elseif typename == "UInt8"
        type.datatype = fmi3UInt8
    elseif typename == "Int16"
        type.datatype = fmi3Int16
    elseif typename == "UInt16"
        type.datatype = fmi3UInt16
    elseif typename == "Int32"
        type.datatype = fmi3Int32
    elseif typename == "UInt32"
        type.datatype = fmi3UInt32
    elseif typename == "Int64"
        type.datatype = fmi3Int64
    elseif typename == "UInt64"
        type.datatype = fmi3UInt64
    elseif typename == "Boolean"
        type.datatype = fmi3Boolean
    elseif typename == "Binary"
        type.datatype = fmi3Binary
    elseif typename == "Char"
        type.datatype = fmi3Char
    elseif typename == "String"
        type.datatype = fmi3String
    elseif typename == "Byte"
        type.datatype = fmi3Byte
    elseif typename == "Enum"
        type.datatype = fmi3Enum
    else
        @warn "Datatype for the variable $(node["name"]) is unknown!"
    end

    if haskey(node, "declaredType")
        type.declaredType = node["declaredType"]
    end

    # if haskey(node, "initial")
    #     for i in 0:(length(instances(fmi3initial))-1)
    #         if "fmi3" * node["initial"] == string(fmi3initial(i))
    #             type.initial = fmi3initial(i)
    #         end
    #     end
    # end

    if haskey(node, "start")
        if node.firstelement !== nothing && node.firstelement.name == "Dimension"
            substrings = split(node["start"], " ")
            if typename == "Float32"
                type.start = Array{fmi3Float32}(undef, 0)
                for string in substrings
                    push!(type.start, parse(fmi3Float32, string))
                end
            elseif typename == "Float64"
                type.start = Array{fmi3Float64}(undef, 0)
                for string in substrings
                    push!(type.start, parse(fmi3Float64, string))
                end
            elseif typename == "Int32"
                type.start = Array{fmi3Int32}(undef, 0)
                for string in substrings
                    push!(type.start, parse(fmi3Int32, string))
                end
            elseif typename == "UInt32"
                type.start = Array{fmi3UInt32}(undef, 0)
                for string in substrings
                    push!(type.start, parse(fmi3UInt32, string))
                end
            elseif typename == "Int64"
                type.start = Array{fmi3Int64}(undef, 0)
                for string in substrings
                    push!(type.start, parse(fmi3Int64, string))
                end
            elseif typename == "UInt64"
                type.start = Array{fmi3UInt64}(undef, 0)
                for string in substrings
                    push!(type.start, parse(fmi3UInt64, string))
                end
            else
                @warn "More array variable types not implemented yet!"
            end
        else
            if typename == "Float32"
                type.start = parse(fmi3Float32, node["start"])
            elseif typename == "Float64"
                type.start = parse(fmi3Float32, node["start"])
            elseif typename == "Int8"
                type.start = parse(fmi3Int8, node["start"])
            elseif typename == "UInt8"
                type.start = parse(fmi3UInt8, node["start"])
            elseif typename == "Int16"
                type.start = parse(fmi3Int16, node["start"])
            elseif typename == "UInt16"
                type.start = parse(fmi3UInt16, node["start"])
            elseif typename == "Int32"
                type.start = parse(fmi3Int32, node["start"])
            elseif typename == "UInt32"
                type.start = parse(fmi3UInt32, node["start"])
            elseif typename == "Int64"
                type.start = parse(fmi3Int64, node["start"])
            elseif typename == "UInt64"
                type.start = parse(fmi3UInt64, node["start"])
            elseif typename == "Boolean"
                type.start = parseFMI3Boolean(node["start"])
            elseif typename == "Binary"
                type.start = pointer(node["start"])
            elseif typename == "Char"
                type.start = parse(fmi3Char, node["start"])
            elseif typename == "String"
                type.start = parse(fmi3String, node["start"])
            elseif typename == "Byte"
                type.start = parse(fmi3Byte, node["start"])
            elseif typename == "Enum"
                for i in 1:length(md.enumerations)
                    if type.declaredType == md.enumerations[i][1] # identify the enum by the name
                        type.start = md.enumerations[i][1 + parse(Int, node["start"])] # find the enum value and set it
                    end
                end
            else
                @warn "setDatatypeVariables(...) unimplemented start value type $typename"
                type.start = node["start"]
            end
        end
    end
    if haskey(node, "intermediateUpdate")
        type.intermediateUpdate = fmi3True
    end

    if haskey(node, "min") && (type.datatype != fmi3Binary || type.datatype != fmiBoolean)
        if type.datatype == fmi3Float32 || type.datatype == fmi3Float64
            type.min = parse(fmi3Float64, node["min"])
        elseif type.datatype == fmi3Enum
            type.min = parse(fmi3Int64, node["min"])
        elseif type.datatype == fmi3Int8 || type.datatype == fmi3Int16 || type.datatype == fmi3Int32 || type.datatype == fmi3Int64
            type.min = parse(fmi3Int32, node["min"])
        else
            type.min = parse(fmi3UInt32, node["min"])
        end
    end
    if haskey(node, "max") && (type.datatype != fmi3Binary || type.datatype != fmiBoolean)
        if type.datatype == fmi3Float32 || type.datatype == fmi3Float64
            type.max = parse(fmi3Float64, node["max"])
        elseif type.datatype == fmi3Enum
            type.max = parse(fmi3Int64, node["max"])
        elseif type.datatype == fmi3Int8 || type.datatype == fmi3Int16 || type.datatype == fmi3Int32 || type.datatype == fmi3Int64
            type.max = parse(fmi3Int32, node["max"])
        else
            type.max = parse(fmi3UInt32, node["max"])
        end
    end
    if haskey(node, "quantity") && (type.datatype != Boolean || type.datatype != fmi3Binary)
        type.quantity = node["quantity"]
    end
    if haskey(node, "unit") && (type.datatype == fmi3Float32 || type.datatype == fmi3Float64)
        type.unit = node["unit"]
    end
    if haskey(node, "displayUnit") && (type.datatype == fmi3Float32 || type.datatype == fmi3Float64)
        type.displayUnit = node["displayUnit"]
    end
    if haskey(node, "relativeQuantity") && (type.datatype == fmi3Float32 || type.datatype == fmi3Float64)
        type.relativeQuantity = parseFMI3Boolean(node["relativeQuantity"])
    else
        type.relativeQuantity = fmi3False
    end
    if haskey(node, "nominal") && (type.datatype == fmi3Float32 || type.datatype == fmi3Float64)
        type.nominal = parse(fmi3Float64, node["nominal"])
    end
    if haskey(node, "unbounded") && (type.datatype == fmi3Float32 || type.datatype == fmi3Float64)
        type.unbounded = parseFMI3Boolean(node["unbounded"])
    else
        type.unbounded = fmi3False
    end
    if haskey(node, "derivative") && (type.datatype == fmi3Float32 || type.datatype == fmi3Float64)
        type.derivative = parse(fmi3UInt32, node["derivative"])
    end
    if haskey(node, "reinit") && (type.datatype == fmi3Float32 || type.datatype == fmi3Float64)
        type.reinit = parseFMI3Boolean(node["reinit"])
    end
    if haskey(node, "mimeType") && type.datatype == fmi3Binary
        type.mimeType = node["mimeType"]
    else
        type.mimeType = "application/octet"
    end
    if haskey(node, "maxSize") && type.datatype == fmi3Binary
        type.maxSize = parse(fmi3UInt32, node["maxSize"])
    end
    type
end

#=
Read all enumerations from the modeldescription and store them in a matrix. First entries are the enum names
-------------------------------------------
Example:
"enum1name" "value1"    "value2"
"enum2name" "value1"    "value2"
=#
function fmi3createEnum(node::EzXML.Node)
    enum = 1
    idx = 1
    enumerations = []
    for simpleType in eachelement(node)
        name = simpleType["name"]
        for type in eachelement(simpleType)
            if type.name == "Enumeration"
                enum = []
                push!(enum, name)
                for item in eachelement(type)
                    push!(enum, item["name"])
                end
                push!(enumerations, enum)
            end
        end
    end
    enumerations
end

################################
# [Sec. 2] get value functions #
################################

"""
Returns startTime from DefaultExperiment if defined else defaults to nothing.
"""
function fmi3GetDefaultStartTime(md::fmi3ModelDescription)
    if md.defaultExperiment == nothing
        return nothing
    end
    return md.defaultExperiment.startTime
end

"""
Returns stopTime from DefaultExperiment if defined else defaults to nothing.
"""
function fmi3GetDefaultStopTime(md::fmi3ModelDescription)
    if md.defaultExperiment == nothing
        return nothing
    end
    return md.defaultExperiment.stopTime
end

"""
Returns tolerance from DefaultExperiment if defined else defaults to nothing.
"""
function fmi3GetDefaultTolerance(md::fmi3ModelDescription)
    if md.defaultExperiment == nothing
        return nothing
    end
    return md.defaultExperiment.tolerance
end

"""
Returns stepSize from DefaultExperiment if defined else defaults to nothing.
"""
function fmi3GetDefaultStepSize(md::fmi3ModelDescription)
    if md.defaultExperiment == nothing
        return nothing
    end
    return md.defaultExperiment.stepSize
end

"""
Returns the tag 'modelName' from the model description.
"""
function fmi3GetModelName(md::fmi3ModelDescription)#, escape::Bool = true)
    md.modelName
end

"""
Returns the tag 'instantionToken' from the model description.
"""
function fmi3GetInstantiationToken(md::fmi3ModelDescription)
    md.instantiationToken
end

"""
Returns the tag 'generationtool' from the model description.
"""
function fmi3GetGenerationTool(md::fmi3ModelDescription)
    md.generationTool
end

"""
Returns the tag 'generationdateandtime' from the model description.
"""
function fmi3GetGenerationDateAndTime(md::fmi3ModelDescription)
    md.generationDateAndTime
end

"""
Returns the tag 'varaiblenamingconvention' from the model description.
"""
function fmi3GetVariableNamingConvention(md::fmi3ModelDescription)
    md.variableNamingConvention
end

"""
Returns the number of EventIndicators from the model description.
"""
function fmi3GetNumberOfEventIndicators(md::fmi3ModelDescription)
    md.numberOfEventIndicators
end

"""
Returns true, if the FMU supports co simulation
"""
function fmi3IsCoSimulation(md::fmi3ModelDescription)
    return( md.coSimulation !== nothing)
end

"""
Returns true, if the FMU supports model exchange
"""
function fmi3IsModelExchange(md::fmi3ModelDescription)
    return( md.modelExchange !== nothing)
end
"""
Returns true, if the FMU supports scheduled execution
"""
function fmi3IsScheduledExecution(md::fmi3ModelDescription)
    return( md.scheduledExecution !== nothing)
end

##################################
# [Sec. 3] information functions #
##################################

"""
Returns the tag 'modelIdentifier' from CS or ME section.
"""
function fmi3GetModelIdentifier(md::fmi3ModelDescription)
    if fmi3IsCoSimulation(md)
        return md.coSimulation.modelIdentifier
    elseif fmi3IsModelExchange(md)
        return md.modelExchange.modelIdentifier
    else
        @assert false "fmi3GetModelName(...): FMU does not support ME or CS!"
    end
end

"""
Returns true, if the FMU supports the getting/setting of states
"""
function fmi3CanGetSetState(md::fmi3ModelDescription)
    return (md.coSimulation != nothing && md.coSimulation.canGetAndSetFMUstate) || (md.modelExchange != nothing && md.modelExchange.canGetAndSetFMUstate)

end

"""
Returns true, if the FMU state can be serialized
"""
function fmi3CanSerializeFMUState(md::fmi3ModelDescription)
    return (md.coSimulation != nothing && md.coSimulation.canSerializeFMUstate) || (md.modelExchange != nothing && md.modelExchange.canSerializeFMUstate)

end

"""
Returns true, if the FMU provides directional derivatives
"""
function fmi3ProvidesDirectionalDerivatives(md::fmi3ModelDescription)
    return (md.coSimulation != nothing && md.coSimulation.providesDirectionalDerivatives) || (md.modelExchange != nothing && md.modelExchange.providesDirectionalDerivatives)
end

"""
Returns true, if the FMU provides adjoint derivatives
"""
function fmi3ProvidesAdjointDerivatives(md::fmi3ModelDescription)
    return (md.coSimulation != nothing && md.coSimulation.providesAdjointDerivatives) || (md.modelExchange != nothing && md.modelExchange.providesAdjointDerivatives)

end


text:
@testset "Diffusion" begin
    var = Variable{Float64}(1, "u")

    diffusion = Diffusion(var)

    # Generate shape function and Variable values for a dummy element
    qr = QuadratureRule{2,RefCube}(:legendre, 2)
    func_space = Lagrange{2, RefCube, 1}()
    fe_values = FECellValues(qr, func_space)
    n_qp = length(points(qr))

    # Nodal coordinates
    coords = Vec{2, Float64}[Vec{2}((-1, -1)),
                             Vec{2}((1, -1)),
                             Vec{2}((1, 1)),
                             Vec{2}((-1, 1))]

    reinit!(fe_values, coords)

    nodal_values = [0., 1., 1., 0.]

    MOOSE.reinit!(var, fe_values, [0,0,0,0], nodal_values)

    residual = zeros(Float64, var.n_dofs)

    # Test the Kernel Interface
    MOOSE.computeResidual!(residual, diffusion)

    # Verified by running the same problem in MOOSE
    moose_residual = [-0.5, 0.5, 0.5, -0.5]

    for i in 1:length(residual)
        @test abs(residual[i]-moose_residual[i]) < 1e-10
    end

    jacobian = zeros(Float64, var.n_dofs, var.n_dofs)

    MOOSE.computeJacobian!(jacobian, diffusion, var)

    moose_jac = [0.666666666667 -0.166666666667 -0.333333333333 -0.166666666667;
                 -0.166666666667 0.666666666667 -0.166666666667 -0.333333333333;
                 -0.333333333333 -0.166666666667 0.666666666667 -0.166666666667;
                 -0.166666666667 -0.333333333333 -0.166666666667 0.666666666667]

    for i in 1:length(jacobian)
        @test abs(jacobian[i] - moose_jac[i]) < 1e-9
    end
end

@testset "DiffusionFD" begin
    var = Variable{Dual{4, Float64}}(1, "u")

    diffusion = Diffusion(var)

    # Generate shape function and Variable values for a dummy element
    qr = QuadratureRule{2,RefCube}(:legendre, 2)
    func_space = Lagrange{2, RefCube, 1}()
    fe_values = FECellValues(qr, func_space)
    n_qp = length(points(qr))

    # Nodal coordinates
    coords = Vec{2, Float64}[Vec{2}((-1, -1)),
                             Vec{2}((1, -1)),
                             Vec{2}((1, 1)),
                             Vec{2}((-1, 1))]

    reinit!(fe_values, coords)

    nodal_values = [MOOSE.dualVariable(0., 1, 4),
                    MOOSE.dualVariable(1., 2, 4),
                    MOOSE.dualVariable(1., 3, 4),
                    MOOSE.dualVariable(0., 4, 4)]

    MOOSE.reinit!(var, fe_values, [0,0,0,0], nodal_values)

    residual = zeros(Dual{4, Float64}, var.n_dofs)

    # Test the Kernel Interface
    MOOSE.computeResidual!(residual, diffusion)

    # Verified by running the same problem in MOOSE
    moose_residual = [-0.5, 0.5, 0.5, -0.5]

    for i in 1:length(residual)
        @test abs(residual[i]-moose_residual[i]) < 1e-10
    end

    jacobian = zeros(Float64, var.n_dofs, var.n_dofs)

    MOOSE.computeJacobian!(jacobian, diffusion, var)

    moose_jac = [0.666666666667 -0.166666666667 -0.333333333333 -0.166666666667;
                 -0.166666666667 0.666666666667 -0.166666666667 -0.333333333333;
                 -0.333333333333 -0.166666666667 0.666666666667 -0.166666666667;
                 -0.166666666667 -0.333333333333 -0.166666666667 0.666666666667]

    for i in 1:length(jacobian)
        @test abs(jacobian[i] - moose_jac[i]) < 1e-9
    end
end


text:
include("mech.jl")

struct QueryException <: Exception
     msg::String
end

Base.showerror(io::IO, e::QueryException) = print(io, e.msg)

################### Query interface
abstract type Query end

# answers the query
answer(q::Query, data) = throw(QueryException("Not Implemented"))

# returns the query sensitivity
sens(q::Query) = throw(QueryException("Not IMplemented"))

# returns a privacy preserving answer
function protect(q::Query, data::Array, epsilon::Number)
    true_answer = answer(q, data)
    protected_answer = zeros(size(true_answer))
    for i in eachindex(true_answer)
        protected_answer[i] = laplace_mechanism(true_answer[i], sensitivity=sens(q), epsilon=epsilon)
    end
    protected_answer
end

#################### Marginal Queries

struct Marginal <: Query
    keepdims::Tuple
end

function answer(q::Marginal, data)
    numdims = length(size(data))
    sumdims = [i for i in 1:numdims if !(i in q.keepdims)]
    sum(data, dims=sumdims)
end

sens(q::Marginal) = 1

################## Range Queries

struct RangeQuery <: Query
    left::Int
    right::Int
end

answer(q::RangeQuery, data) = [sum(data[q.left:q.right])]
sens(q::RangeQuery) = 1

################## Matrix Query

struct MatrixQuery <: Query
    mat::Matrix
end

answer(q::MatrixQuery, data) = q.mat * reshape(data, length(data))
sens(q::MatrixQuery) = maximum(sum(abs.(q.mat), dims=1))
function Base.show(io::IO, q::MatrixQuery)
    println(io, "MatrixQuery")
    show(io, "text/plain", q.mat)
end

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

struct IDQuery <: Query
end

answer(q::IDQuery, data) = copy(data)
sens(q::IDQuery) = 1

##########################################
# Random range queries on 1-d dataset
#########################################
""" Returns a vector of `numqueries' random range queries on a 1-d dataset with d cells"""
function random_range_queries(numqueries, d)
    bounds = [sort([rand(1:d), rand(1:d)]) for _ in 1:numqueries]
    workload = [RangeQuery(x[1], x[2]) for x in bounds]
end

all_range_queries(d) = [RangeQuery(i,j) for i in 1:d for j in i:d]

###########################################
# Tree queries on 1-d dataset
##########################################

# tree
""" Returns a matrix query representing tree queries without the root """
function tree_queries(d, starting=1, ending=d; numsplits=2)
    mat = Vector{Int64}()
    splitsize = ceil(Int64, (ending-starting+1)/numsplits)
    for j in 1:numsplits
        left = (j-1) * splitsize + starting
        right = min(ending, left + splitsize - 1)
        tmp = zeros(Int64, 1, d)
        if left <= ending
            tmp[left:right] .= 1
            mat = cat(mat, tmp, dims=1)
            if right > left
                tmp2 = tree_queries(d, left, right, numsplits=numsplits)
                mat = cat(mat, tmp2.mat, dims=1)
            end
        end
    end
    MatrixQuery(mat)
end

max_stars_repo_name: dkifer/dpcourse

max_issues_repo_name: dkifer/dpcourse

max_forks_repo_name: dkifer/dpcourse


text:
const md_and_friends = [Symbol("@md_str"), Symbol("@html_str"), :getindex]

"""Does the cell only contain md"..." and html"..."?

This is used to run these cells first."""
function is_just_text(topology::NotebookTopology, cell::Cell)::Bool
	# https://github.com/fonsp/Pluto.jl/issues/209
	isempty(topology.nodes[cell].definitions) && isempty(topology.nodes[cell].funcdefs_with_signatures) &&
		topology.nodes[cell].references ⊆ md_and_friends &&
		no_loops(ExpressionExplorer.maybe_macroexpand(topology.codes[cell].parsedcode; recursive=true))
end

function no_loops(ex::Expr)
	if ex.head  === :while ||
		ex.head === :for ||
		ex.head === :comprehension ||
		ex.head === :generator ||
		ex.head === :try ||
		ex.head === :quote ||
		ex.head === :module
		false
	else
		all(no_loops, ex.args)
	end
end

no_loops(x) = true


text:
using Documenter, DDD
using DDD: removeNode!, removeConnection!, removeLink!

makedocs(;
    modules = [DDD],
    format = Documenter.HTML(),
    pages = ["Home" => "index.md", "API" => "api.md"],
    repo = "https://github.com/dcelisgarza/DDD.jl/blob/{commit}{path}#L{line}",
    sitename = "DDD.jl",
    authors = "Daniel Celis Garza",
    # assets=String[],
)

deploydocs(; repo = "github.com/dcelisgarza/DDD.jl.git")

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
#feedback controller (backstepping)
#inputs target is a 4 element vector containing: [x_target, v_target, phi0_target, phi0_dot_target]
function chain_feedback(chain::Chain , target_state::Array{Float64, 1}; alpha = 0.03, beta = 10, c = 10)

  #other parameters
  const g  = 9.81
  const rho = 2.3
  const L = 1.178
  const M = 6.1511

  #calculate difference between actual and desired state
  x_e = chain.q[1]-target_state[1]
  v_e = chain.u[1]-target_state[2]
  phi0_e = chain.q[2]-target_state[3]
  phi0_dot_e = chain.u[2]-target_state[4]

  #use it to generate output force
  r = -alpha*(-rho*g*L*phi0_e+c*x_e)
  r_dot = -alpha*(-rho*g*L*phi0_dot_e+c*v_e)
  ve_dot = r_dot + (1/alpha+beta)*r-beta*v_e
  F = M*ve_dot-rho*g*L*phi0_e

  return F

end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using EffectiveWaves, Test

# This case is numerically challenging, because wavenumber() has many roots close together. Make sure spacing in ω is small to help the optimisation method
@testset "strong scatterers and low frequency" begin
    medium = Acoustic(2; ρ=1.0, c=1.0)
    spatial_dim = 2

    ms = MultipleScattering # just in case Circle conflicts with a definition from another package.

    s1 = Specie(
        Acoustic(spatial_dim; ρ=5.0, c=0.002),ms.Circle(0.004);
        volume_fraction=0.2
    )
    s2 = Specie(
        Acoustic(spatial_dim; ρ=0.3, c=0.01),ms.Circle(0.002);
        volume_fraction=0.1
    )
    species = [s1, s2]
    tol=1e-7
    # ωs = [0.001,0.003]
    ωs = [0.001]
    eff_medium = effective_medium(medium, species)
    k_eff_lows = ωs./eff_medium.c

    k_eff_φs = wavenumber_low_volumefraction(ωs, medium, species, basis_order=1)

    k_effs_arr = [
        wavenumbers(ω, medium, species; tol=tol, num_wavenumbers=1, basis_order=1)
        # num_wavenumbers = 1 usually finds the wavenubmer with the smallest attenuation
    for ω in ωs]

    inds = [argmin(abs.(k_effs_arr[i] .- k_eff_φs[i])) for i in eachindex(ωs)]
    k_effs2 = [k_effs_arr[i][inds[i]] for i in eachindex(inds)]

    @test norm(k_effs2 - k_eff_lows)/norm(k_effs2) < 5e-7
    @test norm(k_effs2[1] - k_eff_lows[1])/norm(k_eff_lows[1]) < 5e-7
    @test norm(k_effs2 - k_eff_φs)/norm(k_effs2) < 0.01
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
const SLIDER_STEPS = 500
const MAX_COLUMNS = 3

"""
    InteractModel(f, model)

An [`AbstractModel`](@ref) that generates its own Interact.jl interface.
Each model [`Param`](@ref) has a slider generated for it to update the model.
After any slider updates the user-defined function `f` is passed the updated
model to generate a new output - anything that will display in a WebIO.jl
node, like a Plots.jl plot

After any slider changes, the parent model is updated, so that `parent(model)`
will return it's latest state.

## Arguments

- `f`: a function that take the model object as an argument, but with
  `Param` fields replaced with their values. Usually a `do` block.
- `model`: any object with [`Param`](@ref)s objects in some fields.

## Param fields

`Param` objects in the model need to include a `range` or `bounds`
field to define the slider range - holding an `AbsractRange` or `NTuple{2}`, respectively.

Optionally, the `Param`s can also include:
- A `label` field to use instead of field names
- A `desciption` field to use in mouse hover text

## Keyword Arguments

- `title`: `""` set a window title, if you need it.
- `submodel`: `Nothing`. `Type` or `Union` that will group sliders into labeled subsections.
- `throttle`: `0.1`. Slider throttle, in seconds. Adjust to improve performance.
- `layout`: `vbox`. This can be any three-argument function that will combine title,
  output, and sliders (all WebIO nodes) into a combined WebIO node. You can use this
  method to do any additional layout you need.

## Example

This is a simple example where the model is a `NamedTuple`, adapted from
an Interact.jl example. It will display in atoms plot pane:

```julia
using InteractModels, Interact, ColorSchemes, Colors

color(i) = Colors.hex(colors[i%length(colors)+1])
colors = ColorSchemes.viridis
width, height = 700, 300
nsamples = 256

model = (;
    sample_step=Param(val=0.05, range=0.01:0.001:0.1, label="Sample step"),
    phase=Param(val=0.0, range=0:0.1:2pi, label="Phase"),
    radii=Param(val=20,range=0:0.1:60, label="Radus")
)

ui = InteractModel(model; submodel=Nothing, throttle=0.1) do m
    cxs_unscaled = [i * m.sample_step + m.phase for i in 1:nsamples]
    cys = sin.(cxs_unscaled) .* height/3 .+ height/2
    cxs = cxs_unscaled .* width/4pi
    c = (dom"svg:circle[cx=\$(cxs[i]), cy=\$(cys[i]), r=\$(m.radii), fill=#\$(color(i))]"()
         for i in 1:nsamples)
    return dom"svg:svg[width=\$width, height=\$height]"(c...)
end
```
"""
mutable struct InteractModel{F,SM,L,Ti,Th} <: AbstractModel
    f::F
    parent::Any
    submodel::SM
    layout::L
    title::Ti
    throttle::Th
    ui::Any
    function InteractModel(
        f::F, parent; layout::L=vbox, ncolumns=nothing, submodel::SM=Nothing,
        title::Ti="", throttle::Th=0.1
    ) where {F,SM,L,Ti,Th}
        # Partially construct model
        model = new{F,SM,L,Ti,Th}(f, parent, submodel, layout, title, throttle, undef)

        # Initialise an observable output object
        output = Observable(f(parent))

        # Generate sliders and update the model and output when they change
        sliders = attach_sliders!(f, model;
            ncolumns=ncolumns, submodel=submodel, throttle=throttle, obs=output
        )

        # Define the complete user interface
        setfield!(model, :ui, layout(dom"h1"(title), output, sliders))

        return model
    end
end
InteractModel(f, model::Model; kwargs...) = InteractModel(f, parent(model); kwargs...)

ui(m::InteractModel) = getfield(m, :ui)

Base.getproperty(m::InteractModel, key::Symbol) = getindex(m, key::Symbol)
Base.setproperty!(m::InteractModel, key::Symbol, x) = setindex!(m, x, key::Symbol)
Base.display(m::InteractModel) = display(ui(m))

@WebIO.register_renderable(InteractModel) do m
    ui(m)
end

"""
    attach_sliders!(f, model::AbstractModel; submodel=false, throttle=0.1)
    attach_sliders!(model::AbstractModel; submodel=false, throttle=0.1, f=identity)

Internal method that may be useful for creating custom interfaces like `InteractModel`,
without actually using `InteracModel` directly. This interface will be less stable than
`InteractModel`.

Create sliders and attach them to the model so it will be updated when they are moved.

## Arguments

- `f`: a function that accepts the model, stripped of `Param` wrappers, with a return value
  that sets the observable `obs`. Usually this converts it to a plot or other web output.
- `model`: a [`AbstractModel`](@ref)

## Keyword Arguments

- `throttle`: `0.1` - sliders response time, in seconds.
- `submodel`: `Nothing`. `Type` or `Union` that will group sliders into labeled subsections.
- `obs`: An optional observable to be updated when sliders change

Returns a `vbox` holding the slider widgets.
"""
function attach_sliders!(f, model::AbstractModel; kwargs...)
    attach_sliders!(model; kwargs..., f=f)
end
function attach_sliders!(model::AbstractModel;
    ncolumns=nothing, submodel=Nothing, throttle=0.1, obs=nothing, f=identity
)
    length(params(model)) == 0 && return hbox()

    sliderbox = if submodel === Nothing
        objpercol = 3
        sliders, slider_obs = param_sliders(model; throttle=throttle)
        _in_columns(sliders, ncolumns, objpercol)
    else
        objpercol = 1
        sliders, slider_obs = group_sliders(f, model, submodel, obs, throttle)
        _in_columns(sliders, ncolumns, objpercol)
    end

    length(slider_obs) == 0 && return hbox()

    # `map` combining Observables is a little odd, *all* sliders here are splatted to `s`
    combined_obs = map((s...) -> s, slider_obs...)
    if length(sliders) > 0
        on(combined_obs) do values
            try
                model[:val] = stripunits(model, values)
                if obs isa Observable
                    obs[] = f(model)
                end
            catch e
                println(stdout, e)
            end
        end
    end

    return sliderbox
end

function _in_columns(objects, ncolumns, objpercol)
    nobjects = length(objects)
    nobjects == 0 && return hbox()

    if ncolumns isa Nothing
        ncolumns = max(1, min(MAX_COLUMNS, (nobjects - 1) ÷ objpercol + 1))
    end
    npercol = (nobjects - 1) ÷ ncolumns + 1
    cols = collect(objects[(npercol * (i - 1) + 1):min(nobjects, npercol * i)] for i in 1:ncolumns)
    return hbox(map(col -> vbox(col...), cols)...)
end

function param_sliders(model::AbstractModel; throttle=0.1)
    length(params(model)) == 0 && return [], []

    labels = if haskey(model, :label)
        map(model[:label], model[:fieldname]) do n, fn
            n === nothing ? fn : n
        end
    else
        model[:fieldname]
    end
    values = withunits(model)
    ranges = if haskey(model, :range)
        withunits(model, :range)
    elseif haskey(model, :bounds)
        _makerange.(withunits(model, :bounds), values)
    else
        _makerange.(Ref(nothing), values)
    end

    descriptions = if haskey(model, :description)
        model[:description]
    else
        map(x -> "", values)
    end

    # Set mouse hover text
    attributes = map(model[:component], labels, descriptions) do p, n, d
        desc = d == "" ? "" : string(": ", d)
        Dict(:title => "$p.$n $desc")
    end

    sliders = map(values, labels, ranges, attributes) do v, l, r, a
        slider(r; label=string(l), value=v, attributes=a)
    end
    # `map` combining Observables is a little odd, *all* sliders here are splatted to `s`
    slider_obs = Interact.throttle.(throttle, observe.(sliders))

    return sliders, slider_obs
end

function group_sliders(f, model::AbstractModel, submodel::Type, obs, throttle)
    # Group slider observations into a single observable
    submodels = [sm for sm in Flatten.flatten(parent(model), submodel) if length(Flatten.flatten(sm, Param)) > 0]
    length(submodels) == 0 && return [], []

    submodel_sliders = map(sm -> param_sliders(Model(sm); throttle=throttle), submodels)
    slider_groups = map(first, submodel_sliders)
    slider_obs = collect(Iterators.flatten(map(last, submodel_sliders)))

    group_doms = map(submodels, slider_groups) do sm, sl
        dom"div"(dom"h2"(string(nameof(typeof(sm)))), sl...)
    end

    return group_doms, slider_obs
end

function _makerange(bounds::Tuple, val::T) where T
    b1, b2 = map(T, bounds)
    step = (b2 - b1) / SLIDER_STEPS
    return b1:step:b2
end
function _makerange(bounds::Tuple, val::T) where T<:Integer
    b1, b2 = map(T, bounds)
    return b1:b2
end
function _makerange(bounds::Nothing, val)
    return if val == zero(val)
        LinRange(-oneunit(val), onunit(val), SLIDER_STEPS)
    else
        LinRange(zero(val), 2 * val, SLIDER_STEPS)
    end
end
function _makerange(bounds::Nothing, val::Int)
    return if val == zero(val)
        -oneunit(val):onunit(val)
    else
        zero(val):2val
    end
end
_makerange(bounds, val) = error("Can't make a range from Param bounds of $val")


text:
__precompile__()

module QuickTypes

using MacroTools: @capture, prewalk, @match, splitarg, @q, splitdef, combinedef, isdef
import ConstructionBase

export @qmutable, @qstruct
export @qmutable_fp, @qstruct_fp
export @qstruct_np, @qmutable_np
export @qfunctor
export @destruct

const special_kwargs = [:_define_show, :_concise_show]

# These are not exported for now, because they are rather specific extensions.
""" For a type X defined with `@qmutable/@qstruct` and with fields `a, b, c,
...`, `QuickTypes.construct(X, 1, 2, 3...)` is a purely-positional constructor
of `X`. This is useful for writing generic structure traversal. For any object
of a type defined by `@qstruct`, this holds:

    QuickTypes.construct(QuickTypes.roottypeof(o), QuickTypes.fieldsof(o)...) == o
"""
function construct end

""" `roottypeof(obj)` returns the type of obj with generic parametric types. Eg.
`roottypeof(a::SomeType{Int}) -> SomeType{T}`. See `QuickTypes.construct` """
@generated roottypeof(obj_type) = roottype(obj_type)
""" `roottype(typ::Type)` returns the parameterless type. Eg. `roottype(X{A}) => X` """
roottype(typ::Type) = Base.typename(typ).wrapper
type_parameters(typ) = typ.parameters
""" `fieldtypes(typ)` returns the types of the fields of a composite type. """
fieldtypes(typ::Type) = # not type-stable ATM. The generated function seemed to have
                        # fieldnames return [] for some reason. Maybe I should try
                        # (fieldtype(typ, 1), fieldtype(typ, 2), ...)
    tuple((fieldtype(typ, f) for f in fieldnames(typ))...)
""" `tuple_parameters{T<:Tuple}(::Type{T})` returns the type of each element of the
tuple, as `svec(type1, type2, ...)` """
tuple_parameters(::Type{T}) where {T<:Tuple} = type_parameters(T)
@generated tuple_parameters_stable(::Type{T}) where {T<:Tuple} =
    :(error("This is broken; it returns a type-tuples, but those cannot be fully-typed, they are instead Tuple{DataType, DataType, ...}. TODO: perhaps we can return Some{...}?"))
    #tuple(type_parameters(T)...)
@generated type_length(::Type{T}) where T = length(fieldnames(T))
is_mutable(::Type{T}) where T = T.mutable

""" `fieldsof(obj)` returns the fields of `obj` in a tuple.
See also `QuickTypes.construct` """
@generated fieldsof(obj) = :(tuple($([:(obj.$field)
                                      for field in fieldnames(obj)]...)))

type_simple_name(ty::Type)::Symbol = ty.name.name

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

""" `parse_funcall(fcall)`

Parses `fname(args; kwargs)`. Returns `(fname, args, kwargs, constraints)`
"""
function parse_funcall(fcall)
    if @capture(fcall, fun_(args__; kwargs__) do; constraints_ end)
        return (fun, args, kwargs, constraints)
    elseif @capture(fcall, fun_(args__) do; constraints_ end)
        return (fun, args, Any[], constraints)
    end

    if @capture(fcall, fun_(args__; kwargs__))
    elseif @capture(fcall, fun_(args__))
        kwargs = Any[]
    else
        error("Not a funcall: $fcall")
    end

    if !isempty(args) && @capture(args[1], ()->constraints_)
        args = args[2:end]
    else
        constraints = nothing
    end

    return (fun, args, kwargs, constraints)
end

function get_sym(e::Expr)
    @assert e.head==:(::)
    e.args[1]
end
get_sym(e::Symbol) = e

""" Build the Base.show function definition for that type """
function build_show_def(define_show::Bool, concise_show::Bool, name, fields, kwfields)
    if !define_show && !concise_show return nothing end

    quote
        function Base.show(io::IO, obj::$name)
            print(io, $(concise_show ? string(name) : @q(typeof(obj))))
            write(io, "(")
            # Positional args
            $([@q begin
               show(io, obj.$(get_sym(field)));
               $(field==last(fields) ? nothing : @q(write(io, ", ")))
               end
               for field in fields]...)
            # separating semicolon
            $(if !isempty(kwfields)
              @q(write(io, "; ")) end)
            # Keyword args
            $([@q begin
               write(io, $(string(get_sym(kwfield)))); write(io, "=");
               show(io, obj.$(get_sym(kwfield)));
               $(kwfield==last(kwfields) ? nothing : @q(write(io, ", ")))
               end
               for kwfield in kwfields]...)
            write(io, ")")
        end
    end
end

function all_type_vars_present(type_vars, args)
    # If we create an outer constructor whose parameters are not part of the args, we get
    #     > @qstruct Blah{T}()
    #     WARNING: static parameter T does not occur in signature for Type.
    #     The method will not be callable.
    # So we detect this case and don't emit an outer constructor
    s = Set(type_vars)
    for arg in args
        prewalk(arg) do x
            if isa(x, Symbol)
                delete!(s, x)
            end
            x
        end
    end
    return isempty(s)
end

narrow_typeof(t::Type{T}) where {T} = Type{T}
narrow_typeof(t::T) where {T} = T

# Helper for @qmutable/@qstruct
# narrow_types means that
function qexpansion(def, mutable, fully_parametric, narrow_types)
    if !@capture(def, typ_def_ <: parent_type_)
        typ_def = def
        parent_type = :Any
    end
    typ, args, kwargs, constraints = parse_funcall(typ_def)
    typ_def = @q($typ($(args...); $(kwargs...)))
    if fully_parametric
        typ, typ_def, args, kwargs = make_parametric(typ, typ_def, args, kwargs)
    end

    get_type_var(v::Symbol) = v
    get_type_var(e::Expr) = e.args[1]
    if @capture(typ, name_{type_params__})
        parametric = true
        type_vars = map(get_type_var, type_params)
        type_with_vars = @q($name{$(type_vars...)})
    else
        type_vars = []
        type_params = []
        parametric = false
        name = typ
        type_with_vars = name
    end
    # Parse the regular arguments
    fields = Any[]; kwfields = Any[]
    constr_args = Any[]; constr_kwargs = Any[]
    o_constr_args = Any[]; o_constr_kwargs = Any[]
    new_args = Any[]
    arg_names = Symbol[]
    reg_kwargs = Any[] # the passed kwargs, but without _concise_show et al.
    for arg in args
        arg_name, arg_type, slurp, default = splitarg(arg)
        if slurp
            @assert arg_type == :Any "Slurping with type arguments not supported"
            @assert default === nothing "Slurping with default not supported"
            arg_type = :Tuple
            push!(constr_args, arg)
            push!(o_constr_args, arg)
        else
            push!(constr_args,
                  default === nothing ? arg_name : Expr(:kw, arg_name, default))
            push!(o_constr_args, arg_name)
        end
        push!(fields, @q($arg_name::$arg_type))
        push!(new_args, arg_name)
        push!(arg_names, arg_name)
    end
    # Parse keyword-arguments
    define_show = nothing # see after the loop
    concise_show = false # default
    for kwarg in kwargs  # keyword arguments
        arg_name, arg_type, slurp, default = splitarg(kwarg)
        if arg_name == :_define_show
            define_show = default::Bool
            continue
        end
        if arg_name == :_concise_show
            concise_show = default::Bool
            continue
        end
        if slurp
            @assert arg_type == :Any "Slurping with type arguments not supported"
            @assert default === nothing "Slurping with default not supported"
            arg_type = Base.Iterators.Pairs
            push!(new_args, arg_name)
            push!(constr_kwargs, kwarg)
            push!(o_constr_kwargs, kwarg)
        else
            push!(new_args, arg_name)
            push!(constr_kwargs,
                  default === nothing ? arg_name : Expr(:kw, arg_name, default))
            push!(o_constr_kwargs, Expr(:kw, arg_name, arg_name))
        end
        push!(reg_kwargs, kwarg)
        push!(kwfields, @q($arg_name::$arg_type))
        push!(arg_names, arg_name)
    end
    # By default, only define Base.show when there are keyword arguments --- otherwise
    # the native `show` is perfectly sufficient.
    if define_show === nothing; define_show = !isempty(kwfields) end

    # -------------- end of parsing -------------

    if narrow_types
        given_types = [@q($QuickTypes.narrow_typeof($a))
                       for a in (fields..., kwfields...)]
    else
        given_types = type_vars
    end
    inner_constr = quote
        function $type_with_vars($(constr_args...);
                                 $(constr_kwargs...)) where {$(type_params...)}
            $constraints
            return new{$(type_vars...)}($(new_args...))
        end
    end
    straight_constr = :($name($(args...); $(reg_kwargs...)) where {$(type_vars...)} =
                        $name{$(given_types...)}($(o_constr_args...);
                                                 $(o_constr_kwargs...)))
    type_def =
        quote
            Base.@__doc__ $(Expr(:struct,
                                 mutable, Expr(:<:, typ, parent_type),
                                 Expr(:block, fields..., kwfields...,
                                      inner_constr,
                                      ((parametric &&
                                        all_type_vars_present(type_vars, [args; kwargs]))
                                       ? [straight_constr] : [])...)))
        end
    construct_def = quote
         function $QuickTypes.construct(::Type{$name}, $(arg_names...))
             $name($(o_constr_args...);
                   $(o_constr_kwargs...))
         end
         function $QuickTypes.ConstructionBase.constructorof(::Type{<:$name})
             (args...) -> $QuickTypes.construct($name, args...)
         end
    end
    @gensym obj
    unpack_def = quote
        macro $(Symbol(:unpack_, name))(obj_expr)
            esc(Expr(:block,
                     Expr(:(=), $(Expr(:quote, obj)), obj_expr),
                     $([Expr(:quote, :($arg = $obj.$arg)) for arg in arg_names]...)))
        end
    end
    esc(Expr(:toplevel,
             type_def,
             construct_def,
             build_show_def(define_show, concise_show, name, fields, kwfields),
             unpack_def,
             nothing))
end

""" Quick type definition.

```julia
@qstruct Car(size, nwheels::Int=4; brand::String="unnamed") <: Vehicle
```

expands into

```julia
struct Car <: Vehicle
    size
    nwheels::Int
    brand::String
    Car(size, nwheels=4; brand="unnamed") = new(size, nwheels, brand)
end
```

Also supports parametric types: `@qstruct Door{T}(size::T)`. Invariants can be
enforced using do-syntax:
```julia
@qstruct Human(name, height::Float64) do
    @assert height > 0    # arbitrary code, executed before creating the object
end
```

Note: `@qstruct` automatically defines a `Base.show` method for the new type,
unless `_define_show=false` (eg. `@qstruct(x, y; _define_show=false)`).
"""
macro qstruct(def)
    return qexpansion(def, false, false, false)
end

""" Quick mutable struct definition. See ?@qstruct """
macro qmutable(def)
    return qexpansion(def, true, false, false)
end

# -----------------------------------------------------------------------------
# Fully-parametric

# Helper for qstruct_fp
function make_parametric(typ, typ_def, args, kwargs)
    all_types = []
    type_counter = 1
    function new_type(parent)
        new_ty = Symbol(:T, type_counter)
        type_counter += 1
        push!(all_types, @q($new_ty <: $parent))
        return new_ty
    end
    #add_type(field::Symbol) = @q($field::$(new_type()))
    function add_type(field)
        name, parent_type, slurp, val = splitarg(field)
        @assert !slurp "Slurping not supported. TODO"
        if name in special_kwargs
            return field
        elseif val==nothing
            return @q($name::$(new_type(parent_type)))
        else
            return Expr(:kw, @q($name::$(new_type(parent_type))), val)
        end
    end

    typed_args = map(add_type, args)
    typed_kwargs = map(add_type, kwargs)
    new_typ = @q($typ{$(all_types...)})

    if type_counter == 1
        # Has to special-case the "no type parameters" case because of
        # https://github.com/JuliaLang/julia/issues/20878
        return (typ, typ_def, args, kwargs)
    else
        return (new_typ, @q($new_typ($(typed_args...); $(typed_kwargs...))),
                typed_args, typed_kwargs)
    end
end

""" Fully-parametric version of `@qstruct`. `@qstruct_fp Foo(a, b=2)` is like
`@qstruct Foo{T, U}(a::T, B::U=2)` """
macro qstruct_fp(def)
    return qexpansion(def, false, true, false)
end
""" Fully-parametric version of `@qmutable`. `@qmutable_fp Foo(a, b=2)` is like
`@qmutable Foo{T, U}(a::T, B::U=2)` """
macro qmutable_fp(def)
    return qexpansion(def, true, true, false)
end

""" Narrowly-parametric version of `@qstruct`. `@qstruct_np Foo(a, b=2)` is like
`@qstruct Foo{T, U}(a::T, B::U=2)`, but it will additionally specialize on types:
`Foo(Int, 2.0) ==> Foo{Type{Int64},Float64}(Int64, 2.0)` """
macro qstruct_np(def)
    return qexpansion(def, false, true, true)
end
macro qmutable_np(def)
    return qexpansion(def, true, true, true)
end

################################################################################
# @qfunctor

"""
```julia
@qfunctor function Action(verb::Symbol)(what)
     println(verb, "ing of ", what)
end <: AbstractAction
```

is equivalent to

```julia
@qstruct Action(verb::Symbol) <: AbstractAction
function (__self__::Action)(what)
    let verb = __self__.verb
        println(verb, "ing of ", what)
    end
end
```
"""
macro qfunctor(fdef0)
    if @capture(fdef0, A_ <: parenttype_)
        fdef = A
    else
        fdef = fdef0
        parenttype = :Any
    end
    di = splitdef(fdef)
    type_def = di[:name]
    if @capture(type_def, typename_(args__; kwargs__))
        all_args = map(first ∘ splitarg, vcat(args, kwargs))
    else
        @assert @capture(type_def, typename_(args__))
        all_args = map(first ∘ splitarg, args)
    end
    di[:name] = :(__self__::$typename)
    if @capture(type_def, tname_{Ts__}(__))
        di[:whereparams] = tuple(Ts..., get(di, :whereparams, ())...)
    end
    di[:body] =
        @q begin
            # I wish I could have used @unpack_Foo, but it seems we can't define a macro and use
            # it in the same top-level block.
            $(Expr(:tuple, all_args...)) =
                $(Expr(:tuple, [:(__self__.$arg) for arg in all_args]...))
            $(di[:body])
        end
    esc(quote
        $QuickTypes.@qstruct $type_def <: $parenttype
        $(combinedef(di))
        end)
end

################################################################################
# @destruct

""" Implement this function for objects that aren't meant to be destructured. Throw an error. """
check_destructurable(x) = nothing

macro destruct_assignment(ass)
    if @capture(ass, typ_(args__; kwargs__) = rhs_)
        nothing
    elseif @capture(ass, typ_(args__) = rhs_)
        kwargs = []
    elseif @capture(ass, (args__,) = rhs_)
        kwargs = []
        typ = Tuple
    else
        @assert @capture(ass, lhs_ = _)  # regular assignment
        return esc(ass)
    end
    obj = rhs isa Symbol ? rhs : gensym(:obj)  # to avoid too many gensyms
    body = []
    for (i, a) in enumerate(args)
        push!(body, :($QuickTypes.@destruct_assignment $a = $Base.getfield($obj, $i)))
    end
    for x in kwargs
        local_var, prop = @capture(x, a_ = b_) ? (a, b) : (x, x)
        prop::Symbol
        push!(body, :($local_var = $obj.$prop))
    end
    esc(@q begin
        $obj = $rhs::$typ
        $QuickTypes.check_destructurable($obj)
        $(body...)
        end)
end

macro destruct_function(fdef)
    di = splitdef(fdef)
    prologue = []
    function proc_arg(a)
        if @capture(a, f_(__))
            @gensym g
            push!(prologue, :($QuickTypes.@destruct_assignment $a = $g))
            return :($g::$f)
        else
            return a
        end
    end
    if haskey(di, :name) # anonymous functions don't have names
        di[:name] = proc_arg(di[:name])     # handle callable objects
    end
    di[:args] = map(proc_arg, di[:args])
    di[:kwargs] = map(proc_arg, get(di, :kwargs, []))
    di[:body] = @q begin
        $(prologue...)
        $(di[:body])
    end
    return esc(combinedef(di))
end

""" Destructuring for objects.

```julia
struct House
    owner
    price
    n_windows
end

@destruct function energy_cost(House(o; n_windows))
    return o == "Bob" ? 10000 : n_windows * 5
end
```

becomes

```julia
function energy_cost(temp_var::House)
    o = getfield(temp_var, 1)         # get the first field
    n_windows = temp_var.n_windows

    return o == "Bob" ? 10000 : n_windows * 5
end
```

This enables syntax like `@destruct mean_price(DataFrame(; price)) = mean(price)`. Destructuring
can also be applied to assignments with `@destruct Ref(x) := ...` and `for` loops. It can be nested:
`@destruct energy_cost(House(Landlord(name, age))) = ...`

Type annotations on fields _do not participate in dispatch_, but are instead converted to.

```julia
julia> @d foo(Ref(a::Int)) = a
foo (generic function with 1 method)

julia> foo(Ref(2.0))
2  # not 2.0
```

`@d ...` is a synonym for `@destruct`. Import it with `using QuickTypes: @d`.
"""
macro destruct(expr::Expr)
    if @capture(expr, lhs_ := rhs_)
        esc(:($QuickTypes.@destruct_assignment $lhs = $rhs))
    elseif @capture(expr, for x_ in seq_ body__ end)
        @gensym g
        esc(@q begin
            for $g in $seq
                $QuickTypes.@destruct_assignment $x = $g
                $(body...)
            end
            end)
    elseif isdef(expr)
        esc(:($QuickTypes.@destruct_function $expr))
    else
        error("@destruct does not handle expressions like $expr")
    end
end

""" Short-hand for `@destruct` """
macro d(expr)
    esc(:($QuickTypes.@destruct $expr))
end

end # module


text:
# VCF Reader
# ==========

"""
    VCFReader(input::IO)

Create a data reader of the VCF file format.

# Arguments
* `input`: data source
"""
type VCFReader <: Bio.IO.AbstractReader
    state::Bio.Ragel.State
    header::VCFHeader

    function VCFReader(input::BufferedInputStream)
        reader = new(Bio.Ragel.State(vcf_header_machine.start_state, input), VCFHeader())
        readheader!(reader)
        reader.state.cs = vcf_body_machine.start_state
        return reader
    end
end

function VCFReader(input::IO)
    return VCFReader(BufferedInputStream(input))
end

function Base.eltype(::Type{VCFReader})
    return VCFRecord
end

function Bio.IO.stream(reader::VCFReader)
    return reader.state.stream
end

function header(reader::VCFReader)
    return reader.header
end

# VCF v4.3
info("compiling VCF")
const vcf_metainfo_machine, vcf_record_machine, vcf_header_machine, vcf_body_machine = (function ()
    cat = Automa.RegExp.cat
    rep = Automa.RegExp.rep
    alt = Automa.RegExp.alt
    opt = Automa.RegExp.opt
    delim(x, sep) = opt(cat(x, rep(cat(sep, x))))

    # The 'fileformat' field is required and must be the first line.
    fileformat = let
        key = cat("fileformat")
        key.actions[:enter] = [:mark1]
        key.actions[:exit]  = [:metainfo_tag]

        version = re"[!-~]+"
        version.actions[:enter] = [:mark2]
        version.actions[:exit]  = [:metainfo_val]

        cat("##", key, '=', version)
    end
    fileformat.actions[:enter] = [:anchor]
    fileformat.actions[:exit]  = [:metainfo]

    # All kinds of meta-information line after 'fileformat' are handled here.
    metainfo = let
        tag = re"[0-9A-Za-z_]+"
        tag.actions[:enter] = [:mark1]
        tag.actions[:exit]  = [:metainfo_tag]

        str = re"[ -;=-~][ -~]*"  # does not starts with '<'
        str.actions[:enter] = [:mark2]
        str.actions[:exit]  = [:metainfo_val]

        dict = let
            dictkey = re"[0-9A-Za-z_]+"
            dictkey.actions[:enter] = [:mark1]
            dictkey.actions[:exit]  = [:metainfo_dict_key]

            dictval = let
                quoted   = cat('"', rep(alt(re"[ !#-[\]-~]", "\\\"", "\\\\")), '"')
                unquoted = rep(re"[ -~]" \ re"[\",>]")
                alt(quoted, unquoted)
            end
            dictval.actions[:enter] = [:mark1]
            dictval.actions[:exit]  = [:metainfo_dict_val]

            cat('<', delim(cat(dictkey, '=', dictval), ','), '>')
        end
        dict.actions[:enter] = [:mark2]
        dict.actions[:exit]  = [:metainfo_val]

        cat("##", tag, '=', alt(str, dict))
    end
    metainfo.actions[:enter] = [:anchor]
    metainfo.actions[:exit]  = [:metainfo]

    # The header line.
    header = let
        sampleID = re"[ -~]+"
        sampleID.actions[:enter] = [:mark1]
        sampleID.actions[:exit]  = [:header_sampleID]

        cat("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO", opt(re"\tFORMAT" * rep(re"\t" * sampleID)))
    end
    header.actions[:enter] = [:anchor]

    # Data lines (fixed fields and variable genotype fields).
    record = let
        chrom = re"[!-9;-~]+"  # no colon
        chrom.actions[:enter] = [:mark]
        chrom.actions[:exit]  = [:record_chrom]

        pos = re"[0-9]+|\."
        pos.actions[:enter] = [:mark]
        pos.actions[:exit]  = [:record_pos]

        id = let
            elm = re"[!-:<-~]+" \ cat('.')
            elm.actions[:enter] = [:mark]
            elm.actions[:exit]  = [:record_id]

            alt(delim(elm, ';'), '.')
        end

        ref = re"[!-~]+"
        ref.actions[:enter] = [:mark]
        ref.actions[:exit]  = [:record_ref]

        alt′ = let
            elm = re"[!-+--~]+" \ cat('.')
            elm.actions[:enter] = [:mark]
            elm.actions[:exit]  = [:record_alt]

            alt(delim(elm, ','), '.')
        end

        qual = re"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?|NaN|[-+]Inf|\."
        qual.actions[:enter] = [:mark]
        qual.actions[:exit]  = [:record_qual]

        filter = let
            elm = re"[!-:<-~]+" \ cat('.')
            elm.actions[:enter] = [:mark]
            elm.actions[:exit]  = [:record_filter]

            alt(delim(elm, ';'), '.')
        end

        info = let
            key = re"[A-Za-z_][0-9A-Za-z_.]*"
            key.actions[:enter] = [:mark]
            key.actions[:exit]  = [:record_info_key]

            val = opt(cat('=', re"[ -:<-~]+"))

            alt(delim(cat(key, val), ';'), '.')
        end

        format = let
            elm = re"[A-Za-z_][0-9A-Za-z_.]*"
            elm.actions[:enter] = [:mark]
            elm.actions[:exit]  = [:record_format]

            alt(delim(elm, ':'), '.')
        end

        genotype = let
            elm = re"[ -9;-~]+"  # no colon
            elm.actions[:enter] = [:mark]
            elm.actions[:exit]  = [:record_genotype_elm]

            delim(elm, ':')
        end
        genotype.actions[:enter] = [:record_genotype]

        cat(
            chrom,  '\t',
            pos,    '\t',
            id,     '\t',
            ref,    '\t',
            alt′,   '\t',
            qual,   '\t',
            filter, '\t',
            info,
            opt(cat('\t', format, rep(cat('\t', genotype)))))
    end
    record.actions[:enter] = [:anchor]
    record.actions[:exit]  = [:record]

    # A newline can be either a CR+LF or a LF.
    newline = let
        lf = re"\n"
        lf.actions[:enter] = [:countline]

        cat(opt('\r'), lf)
    end

    # The VCF file format (header and body part).
    vcfheader = cat(
        fileformat, newline,
        rep(cat(metainfo, newline)),
        header, newline)
    vcfheader.actions[:final] = [:vcfheader]

    vcfbody = rep(cat(record, newline))

    return map(Automa.compile, (metainfo, record, vcfheader, vcfbody))
end)()

const vcf_metainfo_actions = Dict(
    :metainfo_tag      => :(record.tag = (mark1:p-1) - offset),
    :metainfo_val      => :(record.val = (mark2:p-1) - offset; record.dict = data[mark2] == UInt8('<')),
    :metainfo_dict_key => :(push!(record.dictkey, (mark1:p-1) - offset)),
    :metainfo_dict_val => :(push!(record.dictval, (mark1:p-1) - offset)),
    :metainfo          => quote
        Bio.ReaderHelper.resize_and_copy!(record.data, data, offset+1:p-1)
        record.filled = (offset+1:p-1) - offset
    end,
    :anchor            => :(),
    :mark1             => :(mark1 = p),
    :mark2             => :(mark2 = p))
eval(
    Bio.ReaderHelper.generate_index_function(
        VCFMetaInfo,
        vcf_metainfo_machine,
        :(mark1 = mark2 = offset = 0),
        vcf_metainfo_actions))
eval(
    Bio.ReaderHelper.generate_readheader_function(
        VCFReader,
        VCFMetaInfo,
        vcf_header_machine,
        :(mark1 = mark2 = offset = 0),
        merge(vcf_metainfo_actions, Dict(
            :metainfo => quote
                Bio.ReaderHelper.resize_and_copy!(record.data, data, Bio.ReaderHelper.upanchor!(stream):p-1)
                record.filled = (offset+1:p-1) - offset
                @assert isfilled(record)
                push!(reader.header.metainfo, record)
                Bio.ReaderHelper.ensure_margin!(stream)
                record = VCFMetaInfo()
            end,
            :header_sampleID => :(push!(reader.header.sampleID, String(data[mark1:p-1]))),
            :vcfheader => :(finish_header = true; @escape),
            :countline => :(linenum += 1),
            :anchor    => :(Bio.ReaderHelper.anchor!(stream, p); offset = p - 1)))))

const vcf_record_actions = Dict(
    :record_chrom        => :(record.chrom = (mark:p-1) - offset),
    :record_pos          => :(record.pos = (mark:p-1) - offset),
    :record_id           => :(push!(record.id, (mark:p-1) - offset)),
    :record_ref          => :(record.ref = (mark:p-1) - offset),
    :record_alt          => :(push!(record.alt, (mark:p-1) - offset)),
    :record_qual         => :(record.qual = (mark:p-1) - offset),
    :record_filter       => :(push!(record.filter, (mark:p-1) - offset)),
    :record_info_key     => :(push!(record.infokey, (mark:p-1) - offset)),
    :record_format       => :(push!(record.format, (mark:p-1) - offset)),
    :record_genotype     => :(push!(record.genotype, UnitRange{Int}[])),
    :record_genotype_elm => :(push!(record.genotype[end], (mark:p-1) - offset)),
    :record              => quote
        Bio.ReaderHelper.resize_and_copy!(record.data, data, 1:p-1)
        record.filled = (offset+1:p-1) - offset
    end,
    :anchor              => :(),
    :mark                => :(mark = p))
eval(
    Bio.ReaderHelper.generate_index_function(
        VCFRecord,
        vcf_record_machine,
        :(mark = offset = 0),
        vcf_record_actions))
eval(
    Bio.ReaderHelper.generate_read_function(
        VCFReader,
        vcf_body_machine,
        :(mark = offset = 0),
        merge(vcf_record_actions, Dict(
            :record    => quote
                Bio.ReaderHelper.resize_and_copy!(record.data, data, Bio.ReaderHelper.upanchor!(stream):p-1)
                record.filled = (offset+1:p-1) - offset
                found_record = true
                @escape
            end,
            :countline => :(linenum += 1),
            :anchor    => :(Bio.ReaderHelper.anchor!(stream, p); offset = p - 1)))))


text:
## Scalar indexing (and related) conversion helpers

# getting
@propagate_inbounds function Base.getindex(d::AbstractDictionary{I}, i) where {I}
	return d[convert(I, i)]
end

@propagate_inbounds function Base.setindex!(d::AbstractDictionary{I}, v, i) where {I}
    # Should we do a `isequal` check here?
    return setindex!(d, v, convert(I, i))
end

# setting
@propagate_inbounds function Base.setindex!(d::AbstractDictionary{I, T}, v, i::I) where {I, T}
    return setindex!(d, convert(T, v), i)
end

## `get` helper function
@propagate_inbounds function Base.get(d::AbstractDictionary{I, T}, i, default) where {I, T}
    get(d, convert(I, i), default)
end

@propagate_inbounds function Base.get(d::AbstractDictionary{I}, i::I, default) where {I}
    (hasindex, t) = gettoken(d, i)
    if hasindex
        return gettokenvalue(d, t)
    else
        return default
    end
end

@propagate_inbounds function Base.get(f::Base.Callable, d::AbstractDictionary{I}, i) where {I}
    get(f, d, convert(I, i))
end

@propagate_inbounds function Base.get(f::Base.Callable, d::AbstractDictionary{I, T}, i::I) where {I, T}
    (hasindex, t) = gettoken(d, i)
    if hasindex
        return gettokenvalue(d, t)
    else
        return f()
    end
end

## getproperty is equivalent to indexing with a `Symbol`

# @propagate_inbounds Base.getproperty(d::AbstractDictionary, s::Symbol) = d[s]
# @propagate_inbounds function Base.setproperty!(d::AbstractDictionary, s::Symbol, x)
#     d[s] = x
#     return x
# end

## Non-scalar indexing

# Basically, getindices maps the indices over the indexee
@inline function Indexing.getindices(d, inds::AbstractDictionary)
    @boundscheck checkindices(keys(d), inds)
    return map(i -> @inbounds(d[i]), inds)
end

@inline function Indexing.setindices!(d, value, inds::AbstractDictionary)
    @boundscheck checkindices(keys(d), inds)
    foreach(i -> @inbounds(d[i] = value), inds)
    return d
end

@inline function Indexing.getindices(inds1::AbstractIndices, inds2::AbstractIndices)
    @boundscheck checkindices(inds1, inds2)
    return inds2 # TODO should this be a `copy`, perhaps? The output shouldn't alias inds1, but what about inds2?
end

## Views

struct DictionaryView{I, V, Inds <: AbstractDictionary{I}, Vals <: AbstractDictionary{<:Any, V}} <: AbstractDictionary{I, V}
    inds::Inds
    vals::Vals
end

Base.parent(d::DictionaryView) = getfield(d, :vals)
_inds(d::DictionaryView) = getfield(d, :inds)

Base.keys(d::DictionaryView) = keys(_inds(d))

function Base.isassigned(d::DictionaryView{I}, i::I) where {I}
    i2 = @inbounds _inds(d)[i]
    return isassigned(parent(d), i2)
end

@propagate_inbounds function Base.getindex(d::DictionaryView{I}, i::I) where {I}
    i2 = _inds(d)[i]
    return @inbounds parent(d)[i2]
end

@propagate_inbounds function Base.setindex!(d::DictionaryView{I, T}, value::T, i::I) where {I, T}
    i2 = _inds(d)[i]
    return @inbounds parent(d)[i2] = value
end

# `DictionaryView` shares tokens with it's `keys` (and can co-iterate quickly with other dictionaries with those keys)
@propagate_inbounds function gettoken(d::DictionaryView{I}, i::I) where {I}
    return gettoken(_inds(d), i)
end

@propagate_inbounds function istokenassigned(d::DictionaryView, t)
    i2 = gettokenvalue(_inds(d), t)
    return isassigned(parent(d), i2)
end

@propagate_inbounds function gettokenvalue(d::DictionaryView, t)
    i2 = gettokenvalue(_inds(d), t)
    return @inbounds parent(d)[i2]
end

@propagate_inbounds function settokenvalue!(d::DictionaryView{<:Any, T}, t, value::T) where {T}
    i2 = gettokenvalue(_inds(d), t)
    return @inbounds parent(d)[i2] = value
end

@inline function Base.view(vals::AbstractDictionary, inds::AbstractDictionary)
    @boundscheck checkindices(keys(vals), inds)

    return DictionaryView{keytype(inds), eltype(vals), typeof(inds), typeof(vals)}(inds, vals)
end

@inline function Base.view(inds1::AbstractIndices, inds2::AbstractIndices)
    @boundscheck checkindices(inds1, inds2)
    return inds2
end

# TODO accelerate view(::Union{Dictionary, Indices}, ::Indices) to fetch not compute the intermediate `hash`
# (similary a sort-based index could take advantage of sort-merge algorithms on iteration?)


text:
# diff_base.jl - common routines for ordinary symbolic and tensor differentiation

const IDX_NAMES = Espresso.IDX_NAMES

abstract type AbstractDiffRule end

struct DiffRule <: AbstractDiffRule
    pat::Expr        # pattern of expression to differentiate
    dpat::Any        # pattern of differentiation expression
end

const DIFF_PHS = Set([:x, :y, :z, :a, :b, :c, :m, :n])

const DIFF_RULES =
    Dict{Tuple{OpName, Vector{Type}, Int}, DiffRule}()

opname(mod, op) = canonical(mod, op)

"""
Define new differentiation rule. Arguments:

 * `ex` - original expression in a form `func(arg1::Type1, arg2::Type2, ...)`
 * `idx` - index of argument to differentiate over
 * `dex` - expression of corresponding derivative

Example:

    @scalardiff *(x::Number, y::Number) 1 y

Which means: derivative of a product of 2 numbers w.r.t. 1st argument
is a second argument.

Note that rules are always defined as if arguments were ordinary variables
and not functions of some other variables, because this case will be
automatically handled by chain rule in the differentiation engine.

"""
macro scalardiff(ex::Expr, idx::Int, dex::Any)
    @assert ex.head == :call
    op = opname(current_module(), ex.args[1])
    types = [eval(exa.args[2]) for exa in ex.args[2:end]]
    new_args = Symbol[exa.args[1] for exa in ex.args[2:end]]
    canonical_ex = Expr(:call, op, new_args...)
    # canonical_ex = canonical_calls(current_module(), ex)
    canonical_dex = canonical_calls(current_module(), dex)
    DIFF_RULES[(op, types, idx)] = DiffRule(canonical_ex, canonical_dex)
end

"""
Find differentiation rule for `op` with arguments of `types`
w.r.t. `idx`th argument. Example:

    rule = find_rule(:*, [Int, Int], 1)

Which reads as: find rule for product of 2 Ints w.r.t. 1st argument.

In addition to the types passed, rules for all combinations of all their
ansestors (as defined by `type_ansestors()`) will be checked.

Rule itself is an opaque object containing information needed for derivation
and guaranted to be compatible with `apply_rule()`.
"""
function find_rule(op::OpName, types::Vector{DataType}, idx::Int)
    type_ans = map(type_ansestors, types)
    type_products = product(type_ans...)
    ks = ((op, [tp...], idx) for tp in type_products)
    for k in ks
        if haskey(DIFF_RULES, k)
            return Nullable(DIFF_RULES[k])
        end
    end
    return Nullable()
end

"""
Apply rule retrieved using `find_rule()` to an expression.
"""
function apply_rule(rule::DiffRule, ex::Expr)
    deriv_ex = rewrite(ex, rule.pat, rule.dpat; phs=DIFF_PHS)
    return deriv_ex
end

## register rule

function without_output_tuple(ex::Expr)
    @assert(ex.head == :block && ex.args[end].head == :tuple,
            "Got unexpected derivative expression: $ex")
    return Expr(ex.head, ex.args[1:end-1]...)
end

"""
Add to all intermediate variables a prefix of a function name to avoid name conflicts
"""
function with_opname_prefix(ex::Expr, fname::OpName)
    @assert ex.head == :block
    prefix = isa(fname, Symbol) ? "__$(fname)_" : "__$(fname.args[1])_$(fname.args[2].value)_"
    st = Dict()
    for subex in ex.args
        if subex.head == :(=)
            vname = split_indexed(subex.args[1])[1]
            st[vname] = Symbol("$(prefix)_$(vname)")
        end
    end
    return subs(ex, st)
end

"""
Register new differentiation rule for function `fname` with arguments
of `types` at index `idx`, return this new rule.
"""
function register_rule(fname::OpName, types::Vector{DataType}, idx::Int)
    # TODO: check module
    f = eval(fname)
    args, ex = func_expr(f, (types...))
    # ex = sanitize(ex)
    xs = [(arg, ones(T)[1]) for (arg, T) in zip(args, types)]
    dex = xdiff(ex; xs...)
    prefixed_dex = with_opname_prefix(remove_unused(dex), fname)
    pure_dex = without_output_tuple(prefixed_dex)
    fex = Expr(:call, fname, args...)
    new_rule = DiffRule(fex, pure_dex)
    DIFF_RULES[(fname, types, idx)] = new_rule
    return new_rule
end

## derivative (for primitive expressions)

function derivative(pex::Expr, types::Vector{DataType}, idx::Int;
                    mod=current_module())
    @assert pex.head == :call
    op = canonical(mod, pex.args[1])
    maybe_rule = find_rule(op, types, idx)
    if !isnull(maybe_rule)
        rule = get(maybe_rule)
    else
        error("Primitive expression $pex with types $types at index $idx " *
              "doesn't have a registered derivative. Use `@scalardiff` to register one. " *
              "Note: currently automatic derivative inference of nested functions " *
              "is turned off as unreliable, but this may change in the future.")
        # register_rule(op, types, idx)
    end
    return apply_rule(rule, pex) |> sanitize
end

derivative(var::Symbol, types::Vector{DataType}, idx::Int) = 1


text:
module Sim6502

using MacroTools

using FunctionWrappers: FunctionWrapper

include("utils.jl")
include("cpu.jl")
include("memory.jl")
include("boilerplate.jl")
include("instructions.jl")
include("chipset.jl")
include("opcodes.jl")
include("assembler.jl")

export CPU, store!, counter!, describe, status, status!, status_string
export Memory, Π, Π8, Π16, reset!, dereference, deref, ↦, ⇾, store!
export Chipset

export AddressingMode, DirectMode, IndirectMode, Direct, DirectX, DirectY, IndirectX,
    IndirectY

export op!, tick!, @assemble

# these are all instruction pieces which probably shouldn't be exported
export ld!, lda!, ldx!, ldy!, st!, sta!, stx!, sty!, t!, tax!, tay!, txa!, tya!,
    tsx!, txs!, pha!, php!, pla!, plp!
export and!, eor!, ora!, bit!, adc!, sbc!, compare!, cmp!, cpx!, cpy!
export increment!, inc!, inx!, iny!, decrement!, dec!, dex!, dey!
export arithmetic_shiftleft!, asl!, logical_shiftright!, lsr!, rotateleft!, rol!,
    rotateright!, ror!
export jmp!, jsr!, rts!
export bcc!, bcs!, beq!, bmi!, bne!, bpl!, bvc!, bvs!
export clc!, cld!, cli!, clv!, sec!, sed!, sei!, brk!, nop!, rti!

end # module