text:
export NodeInterface, IndexedNodeInterface, name, tag, messageout, messagein
export FactorNode, functionalform, interfaces, factorisation, localmarginals, localmarginalnames, metadata
export iscontain, isfactorised, getinterface
export clusters, clusterindex
export deps, connect!, activate!
export make_node, on_make_node, AutoVar
export ValidNodeFunctionalForm, UndefinedNodeFunctionalForm, as_node_functional_form
export sdtype, Deterministic, Stochastic, isdeterministic, isstochastic
export MeanField, FullFactorisation, collect_factorisation
export DefaultFunctionalDependencies, RequireInboundFunctionalDependencies, RequireEverythingFunctionalDependencies
export @node

using Rocket
using TupleTools

import Base: show, +
import Base: getindex, setindex!, firstindex, lastindex

## Node traits

struct ValidNodeFunctionalForm end
struct UndefinedNodeFunctionalForm end

as_node_functional_form(::Function) = ValidNodeFunctionalForm()
as_node_functional_form(some)       = UndefinedNodeFunctionalForm()

## Node types

"""
    Deterministic

`Deterministic` object used to parametrize factor node object with determinstic type of relationship between variables.

See also: [`Stochastic`](@ref), [`isdeterministic`](@ref), [`isstochastic`](@ref)
"""
struct Deterministic end

"""
    Stochastic

`Stochastic` object used to parametrize factor node object with stochastic type of relationship between variables.

See also: [`Deterministic`](@ref), [`isdeterministic`](@ref), [`isstochastic`](@ref)
"""
struct Stochastic end

"""
    isdeterministic(node)

Function used to check if factor node object is deterministic or not. Returns true or false.

See also: [`Deterministic`](@ref), [`Stochastic`](@ref), [`isstochastic`](@ref)
"""
function isdeterministic end

"""
    isstochastic(node)

Function used to check if factor node object is stochastic or not. Returns true or false.

See also: [`Deterministic`](@ref), [`Stochastic`](@ref), [`isdeterministic`](@ref)
"""
function isstochastic end

isdeterministic(::Deterministic)         = true
isdeterministic(::Type{ Deterministic }) = true
isdeterministic(::Stochastic)            = false
isdeterministic(::Type{ Stochastic })    = false

isstochastic(::Stochastic)            = true
isstochastic(::Type{ Stochastic })    = true
isstochastic(::Deterministic)         = false
isstochastic(::Type{ Deterministic }) = false

"""
    sdtype(object)

Returns either `Deterministic` or `Stochastic` for a given object (if defined).

See also: [`Deterministic`](@ref), [`Stochastic`](@ref), [`isdeterministic`](@ref), [`isstochastic`](@ref)
"""
function sdtype end

sdtype(::Function) = Deterministic()

## Generic factorisation constraints

"""
    MeanField

Generic factorisation constraint used to specify a mean-field factorisation for recognition distribution `q`.

See also: [`FullFactorisation`](@ref)
"""
struct MeanField end

"""
    FullFactorisation

Generic factorisation constraint used to specify a full factorisation for recognition distribution `q`.

See also: [`MeanField`](@ref)
"""
struct FullFactorisation end

"""
    collect_factorisation(nodetype, factorisation)

This function converts given factorisation to a correct internal factorisation representation for a given node.

See also: [`MeanField`](@ref), [`FullFactorisation`](@ref)
"""
function collect_factorisation end

"""
    collect_meta(nodetype, meta)

This function converts given meta object to a correct internal meta representation for a given node.
Fallbacks to `default_meta` in case if meta is `nothing`.

See also: [`default_meta`](@ref), [`FactorNode`](@ref)
"""
function collect_meta end

collect_meta(T::Any, ::Nothing) = default_meta(T)
collect_meta(T::Any, meta::Any) = meta

"""
    default_meta(nodetype)

Returns default meta object for a given node type.

See also: [`collect_meta`](@ref), [`FactorNode`](@ref)
"""
function default_meta end

default_meta(any) = nothing

## NodeInterface

"""
    NodeInterface

`NodeInterface` object represents a single node-variable connection.

See also: [`name`](@ref), [`tag`](@ref), [`messageout`](@ref), [`messagein`](@ref)
"""
mutable struct NodeInterface
    name               :: Symbol
    m_out              :: LazyObservable{AbstractMessage}
    m_in               :: LazyObservable{AbstractMessage}
    connected_variable :: Union{Nothing, AbstractVariable}
    connected_index    :: Int

    NodeInterface(name::Symbol) = new(name, lazy(AbstractMessage), lazy(AbstractMessage), nothing, 0)
end

Base.show(io::IO, interface::NodeInterface) = print(io, string("Interface(", name(interface), ")"))

"""
    name(interface)

Returns a name of the interface.

See also: [`NodeInterface`](@ref), [`tag`](@ref)
"""
name(symbol::Symbol)              = symbol
name(interface::NodeInterface)    = name(interface.name)

"""
    tag(interface)

Returns a tag of the interface in the form of `Val{ name(interface) }`.
The major difference between tag and name is that it is possible to dispath on interface's tag in message computation rule.

See also: [`NodeInterface`](@ref), [`name`](@ref)
"""
tag(interface::NodeInterface)     = Val{ name(interface) }

"""
    messageout(interface)

Returns an outbound messages stream from the given interface.

See also: [`NodeInterface`](@ref), [`messagein`](@ref)
"""
messageout(interface::NodeInterface) = interface.m_out

"""
    messagein(interface)

Returns an inbound messages stream from the given interface.

See also: [`NodeInterface`](@ref), [`messageout`](@ref)
"""
messagein(interface::NodeInterface)  = interface.m_in

"""
    connectvariable!(interface, variable, index)

Connects a variable with the interface and given index. Index is used to distinguish this connection from others in case if variable is connected to multiple interfaces.

See also: [`NodeInterface`](@ref), [`connectedvar`](@ref), [`connectedvarindex`](@ref)
"""
function connectvariable!(interface::NodeInterface, variable, index)
    interface.connected_variable = variable
    interface.connected_index    = index
end

"""
    connectedvar(interface)

Returns connected variable for the interface.

See also: [`NodeInterface`](@ref), [`connectvariable!`](@ref), [`connectedvarindex`](@ref)
"""
connectedvar(interface::NodeInterface)      = interface.connected_variable

"""
    connectedvarindex(interface)

Returns an index of connected variable for the interface.

See also: [`NodeInterface`](@ref), [`connectvariable!`](@ref), [`connectedvar`](@ref)
"""
connectedvarindex(interface::NodeInterface) = interface.connected_index

"""
    get_pipeline_stages(interface)

Returns an instance of pipeline stages of connected variable for the given interface

See also: [`NodeInterface`](@ref), [`connectvariable!`](@ref), [`connectedvar`](@ref), [`add_inbound_pipeline_stage!`](@ref)
"""
get_pipeline_stages(interface::NodeInterface) = get_pipeline_stages(connectedvar(interface))

"""
    IndexedNodeInterface

`IndexedNodeInterface` object represents a repetative node-variable connection.
Used in cases when node may connect different number of random variables with the same name, e.g. means and precisions of Gaussian Mixture node.

See also: [`name`](@ref), [`tag`](@ref), [`messageout`](@ref), [`messagein`](@ref)
"""
struct IndexedNodeInterface
    index     :: Int
    interface :: NodeInterface
end

Base.show(io::IO, interface::IndexedNodeInterface) = print(io, string("IndexedInterface(", name(interface), ",", index(interface), ")"))

name(interface::IndexedNodeInterface)  = name(interface.interface)
index(interface::IndexedNodeInterface) = interface.index
tag(interface::IndexedNodeInterface)   = (Val{ name(interface) }(), index(interface))

messageout(interface::IndexedNodeInterface) = messageout(interface.interface)
messagein(interface::IndexedNodeInterface)  = messagein(interface.interface)

connectvariable!(interface::IndexedNodeInterface, variable, index) = connectvariable!(interface.interface, variable, index)
connectedvar(interface::IndexedNodeInterface)                      = connectedvar(interface.interface)
connectedvarindex(interface::IndexedNodeInterface)                 = connectedvarindex(interface.interface)
get_pipeline_stages(interface::IndexedNodeInterface)               = get_pipeline_stages(interface.interface)

## FactorNodeLocalMarginals

"""
    FactorNodeLocalMarginal

This object represents local marginals for some specific factor node.
Local marginal can be joint in case of structured factorisation.
Local to factor node marginal also can be shared with a corresponding marginal of some random variable.

See also: [`FactorNodeLocalMarginals`](@ref)
"""
mutable struct FactorNodeLocalMarginal
    index  :: Int
    name   :: Symbol
    stream :: Union{Nothing, MarginalObservable}

    FactorNodeLocalMarginal(index::Int, name::Symbol) = new(index, name, nothing)
end

index(localmarginal::FactorNodeLocalMarginal) = localmarginal.index
name(localmarginal::FactorNodeLocalMarginal)  = localmarginal.name

getstream(localmarginal::FactorNodeLocalMarginal) = localmarginal.stream
setstream!(localmarginal::FactorNodeLocalMarginal, observable::MarginalObservable) = localmarginal.stream = observable

"""
    FactorNodeLocalMarginals

This object acts as an iterable and indexable proxy for local marginals for some node.
"""
struct FactorNodeLocalMarginals{M}
    marginals :: M
end

function FactorNodeLocalMarginals(variablenames, factorisation)
    marginal_names = map(fcluster -> clustername(map(i -> variablenames[i], fcluster)), factorisation)
    index          = 0 # its better not to use zip here to preserve tuple-like structure
    marginals      = map((mname) -> begin index += 1; FactorNodeLocalMarginal(index, mname) end, marginal_names)
    return FactorNodeLocalMarginals(marginals)
end

## AbstractFactorNode

abstract type AbstractFactorNode end

isstochastic(factornode::AbstractFactorNode)    = isstochastic(sdtype(factornode))
isdeterministic(factornode::AbstractFactorNode) = isdeterministic(sdtype(factornode))

## Generic Factor Node

struct FactorNode{F, I, C, M, A, P} <: AbstractFactorNode
    fform          :: F
    interfaces     :: I
    factorisation  :: C
    localmarginals :: M
    metadata       :: A
    pipeline       :: P
end

function FactorNode(fform::Type{F}, interfaces::I, factorisation::C, localmarginals::M, metadata::A, pipeline::P) where { F, I, C, M, A, P }
    return FactorNode{Type{F}, I, C, M, A, P}(fform, interfaces, factorisation, localmarginals, metadata, pipeline)
end

function FactorNode(fform, varnames::NTuple{N, Symbol}, factorisation, metadata, pipeline) where N
    interfaces     = map(varname -> NodeInterface(varname), varnames)
    localmarginals = FactorNodeLocalMarginals(varnames, factorisation)
    return FactorNode(fform, interfaces, factorisation, localmarginals, metadata, pipeline)
end

function Base.show(io::IO, factornode::FactorNode)
    println(io, "FactorNode:")
    println(io, string(" form            : ", functionalform(factornode)))
    println(io, string(" sdtype          : ", sdtype(factornode)))
    println(io, string(" interfaces      : ", interfaces(factornode)))
    println(io, string(" factorisation   : ", factorisation(factornode)))
    println(io, string(" local marginals : ", localmarginalnames(factornode)))
    println(io, string(" metadata        : ", metadata(factornode)))
    println(io, string(" pipeline        : ", getpipeline(factornode)))
end

functionalform(factornode::FactorNode)            = factornode.fform
sdtype(factornode::FactorNode)                    = sdtype(functionalform(factornode))
interfaces(factornode::FactorNode)                = factornode.interfaces
factorisation(factornode::FactorNode)             = factornode.factorisation
localmarginals(factornode::FactorNode)            = factornode.localmarginals.marginals
localmarginalnames(factornode::FactorNode)        = map(name, localmarginals(factornode))
metadata(factornode::FactorNode)                  = factornode.metadata
getpipeline(factornode::FactorNode)               = factornode.pipeline

clustername(cluster) = mapreduce(v -> name(v), (a, b) -> Symbol(a, :_, b), cluster)

# Cluster is reffered to a tuple of node interfaces
clusters(factornode::FactorNode) = map(factor -> map(i -> begin return @inbounds interfaces(factornode)[i] end, factor), factorisation(factornode))

clusterindex(factornode::FactorNode, v::Symbol)                              = clusterindex(factornode, (v, ))
clusterindex(factornode::FactorNode, vindex::Int)                            = clusterindex(factornode, (vindex, ))
clusterindex(factornode::FactorNode, vars::NTuple{N, NodeInterface}) where N = clusterindex(factornode, map(v -> name(v), vars))
clusterindex(factornode::FactorNode, vars::NTuple{N, Symbol})        where N = clusterindex(factornode, map(v -> interfaceindex(factornode, v), vars))

function interfaceindex(factornode::FactorNode, iname::Symbol)
    iindex = findfirst(interface -> name(interface) === iname, interfaces(factornode))
    return iindex !== nothing ? iindex : error("Unknown interface ':$(iname)' for $(functionalform(factornode)) node")
end

function clusterindex(factornode::FactorNode, vars::NTuple{N, Int}) where N
    cindex = findfirst(cluster -> all(v -> v ∈ cluster, vars), factorisation(factornode))
    return cindex !== nothing ? cindex : error("Unknown cluster '$(vars)' for $(functionalform(factornode)) node")
end

function varclusterindex(cluster, iindex::Int)
    vcindex = findfirst(index -> index === iindex, cluster)
    return vcindex !== nothing ? vcindex : error("Invalid cluster '$(vars)' for a given interface index '$(iindex)'")
end

getinterface(factornode::FactorNode, iname::Symbol)   = @inbounds interfaces(factornode)[ interfaceindex(factornode, iname) ]

getclusterinterfaces(factornode::FactorNode, cindex::Int) = @inbounds map(i -> interfaces(factornode)[i], factorisation(factornode)[ cindex ])

iscontain(factornode::FactorNode, iname::Symbol)      = findfirst(interface -> name(interface) === iname, interfaces(factornode)) !== nothing
isfactorised(factornode::FactorNode, factor)          = findfirst(f -> f == factor, factorisation(factornode)) !== nothing

function connect!(factornode::FactorNode, iname::Symbol, variable)
    return connect!(factornode::FactorNode, iname::Symbol, variable, getlastindex(variable))
end

function connect!(factornode::FactorNode, iname::Symbol, variable, index)
    vinterface = getinterface(factornode, iname)
    connectvariable!(vinterface, variable, index)
    setmessagein!(variable, index, messageout(vinterface))
end

## Node pipeline

abstract type AbstractNodeFunctionalDependenciesPipeline end

struct FactorNodePipeline{F <: AbstractNodeFunctionalDependenciesPipeline, S <: AbstractPipelineStage}
    functional_dependencies :: F
    extra_stages            :: S
end

get_pipeline_dependencies(pipeline::FactorNodePipeline) = pipeline.functional_dependencies
get_pipeline_stages(pipeline::FactorNodePipeline)       = pipeline.extra_stages

function Base.show(io::IO, pipeline::FactorNodePipeline)
    print(io, "FactorNodePipeline(functional_dependencies = $(pipeline.functional_dependencies), extra_stages = $(pipeline.extra_stages)")
end

function collect_pipeline end

collect_pipeline(T::Any, ::Nothing)                                       = FactorNodePipeline(default_functional_dependencies_pipeline(T), EmptyPipelineStage())
collect_pipeline(T::Any, stage::AbstractPipelineStage)                    = FactorNodePipeline(default_functional_dependencies_pipeline(T), stage)
collect_pipeline(T::Any, fdp::AbstractNodeFunctionalDependenciesPipeline) = FactorNodePipeline(fdp, EmptyPipelineStage())
collect_pipeline(T::Any, pipeline::FactorNodePipeline)                    = pipeline

## Functional Dependencies

function message_dependencies end
function marginal_dependencies end

Base.:+(left::AbstractNodeFunctionalDependenciesPipeline, right::AbstractPipelineStage) = FactorNodePipeline(left, right)
Base.:+(left::FactorNodePipeline, right::AbstractPipelineStage)                         = FactorNodePipeline(left.functional_dependencies, left.extra_stages + right)

### Default

"""
    DefaultFunctionalDependencies
"""
struct DefaultFunctionalDependencies <: AbstractNodeFunctionalDependenciesPipeline end

function message_dependencies(::DefaultFunctionalDependencies, nodeinterfaces, varcluster, iindex)
    # First we remove current edge index from the list of dependencies
    vdependencies = TupleTools.deleteat(varcluster, varclusterindex(varcluster, iindex))
    # Second we map interface indices to the actual interfaces
    return map(inds -> map(i -> begin return @inbounds nodeinterfaces[i] end, inds), vdependencies)
end

function marginal_dependencies(::DefaultFunctionalDependencies, nodelocalmarginals, varcluster, cindex)
    return TupleTools.deleteat(nodelocalmarginals, cindex)
end

### With inbound

struct RequireInboundFunctionalDependencies{I, S} <: AbstractNodeFunctionalDependenciesPipeline
    indices    :: I
    start_with :: S
end

struct InterfacePluginStartWithMessage{M, S}
    msg        :: M
    start_with :: S
end

name(p::InterfacePluginStartWithMessage)      = name(p.msg)
messagein(p::InterfacePluginStartWithMessage) = messagein(p.start_with, p)

messagein(::Nothing, p::InterfacePluginStartWithMessage) = messagein(p.msg)
messagein(something, p::InterfacePluginStartWithMessage) = messagein(p.msg) |> start_with(Message(something, false, true))

function message_dependencies(dependencies::RequireInboundFunctionalDependencies, nodeinterfaces, varcluster, iindex)

    # First we find dependency index in `indices`, we use it later to find `start_with` distribution
    depindex = findfirst((i) -> i === iindex, dependencies.indices)

    # If we have `depindex` in our `indices` we include it in our list of functional dependencies. It effectively forces rule to require inbound message
    if depindex !== nothing
        # `mapindex` is a lambda function here
        mapindex = let nodeinterfaces = nodeinterfaces, depindex = depindex
            (i) -> begin
                interface = @inbounds nodeinterfaces[i]
                # InterfacePluginStartWithMessage is a proxy structure for `name` and `messagein` method for an interface
                # It returns the same name but modifies `messagein` to return an observable with `start_with` operator
                return i === iindex ? InterfacePluginStartWithMessage(interface, dependencies.start_with[depindex]) : interface
            end
        end
        return map(inds -> map(mapindex, inds), varcluster)
    else
        return message_dependencies(DefaultFunctionalDependencies(), nodeinterfaces, varcluster, iindex)
    end
end

function marginal_dependencies(::RequireInboundFunctionalDependencies, nodelocalmarginals, varcluster, cindex)
    return marginal_dependencies(DefaultFunctionalDependencies(), nodelocalmarginals, varcluster, cindex)
end

### Everything

struct RequireEverythingFunctionalDependencies <: AbstractNodeFunctionalDependenciesPipeline end

function ReactiveMP.message_dependencies(::RequireEverythingFunctionalDependencies, nodeinterfaces, varcluster, iindex)
    # Return all node interfaces including the edge we are trying to compuate a message on
    return nodeinterfaces
end

function ReactiveMP.marginal_dependencies(::RequireEverythingFunctionalDependencies, nodelocalmarginals, varcluster, cindex)
    # Returns only local marginals based on local q factorisation, it does not return all possible combinations of all joint posterior marginals
    return nodelocalmarginals
end

###

default_functional_dependencies_pipeline(_) = DefaultFunctionalDependencies()

### Generic

function functional_dependencies(dependencies, factornode::FactorNode, iname::Symbol)
    return functional_dependencies(dependencies, factornode, interfaceindex(factornode, iname))
end

function functional_dependencies(dependencies, factornode::FactorNode, iindex::Int)
    cindex  = clusterindex(factornode, iindex)

    nodeinterfaces     = interfaces(factornode)
    nodeclusters       = factorisation(factornode)
    nodelocalmarginals = localmarginals(factornode)

    varcluster = @inbounds nodeclusters[ cindex ]

    messages  = message_dependencies(dependencies, nodeinterfaces, varcluster, iindex)
    marginals = marginal_dependencies(dependencies, nodelocalmarginals, varcluster, cindex)

    return tuple(messages...), tuple(marginals...)
end

function get_messages_observable(factornode, messages)
    if !isempty(messages)
        msgs_names      = Val{ map(name, messages) }
        msgs_observable = combineLatest(map(m -> messagein(m), messages), PushNew())
        return msgs_names, msgs_observable
    else
        return nothing, of(nothing)
    end
end

function get_marginals_observable(factornode, marginals)
    if !isempty(marginals)
        marginal_names       = Val{ map(name, marginals) }
        marginals_streams    = map(marginal -> getmarginal!(factornode, marginal, IncludeAll()), marginals)
        marginals_observable = combineLatestUpdates(marginals_streams, PushNew())
        return marginal_names, marginals_observable
    else
        return nothing, of(nothing)
    end
end

# Variational Message Passing
apply_mapping(msgs_observable, marginals_observable, mapping) = (dependencies) -> VariationalMessage(dependencies[1], dependencies[2], mapping)

# Fallback for Belief Propagation
apply_mapping(msgs_observable, marginals_observable::SingleObservable{Nothing}, mapping) = mapping

function activate!(model, factornode::AbstractFactorNode)
    fform = functionalform(factornode)
    meta  = metadata(factornode)
    node_pipeline = getpipeline(factornode)

    node_pipeline_dependencies = get_pipeline_dependencies(node_pipeline)
    node_pipeline_extra_stages = get_pipeline_stages(node_pipeline)

    for (iindex, interface) in enumerate(interfaces(factornode))
        message_dependencies, marginal_dependencies = functional_dependencies(node_pipeline_dependencies, factornode, iindex)

        msgs_names, msgs_observable          = get_messages_observable(factornode, message_dependencies)
        marginal_names, marginals_observable = get_marginals_observable(factornode, marginal_dependencies)

        vtag        = tag(interface)
        vconstraint = local_constraint(connectedvar(interface))

        vmessageout = combineLatest((msgs_observable, marginals_observable), PushNew()) # TODO check PushEach
        vmessageout = apply_pipeline_stage(get_pipeline_stages(interface), factornode, vtag, vmessageout)

        mapping = MessageMapping(fform, vtag, vconstraint, msgs_names, marginal_names, meta, factornode)
        mapping = apply_mapping(msgs_observable, marginals_observable, mapping)

        vmessageout = vmessageout |> map(AbstractMessage, mapping)
        vmessageout = apply_pipeline_stage(get_pipeline_stages(getoptions(model)), factornode, vtag, vmessageout)
        vmessageout = apply_pipeline_stage(node_pipeline_extra_stages, factornode, vtag, vmessageout)
        vmessageout = vmessageout |> schedule_on(global_reactive_scheduler(getoptions(model)))

        set!(messageout(interface), vmessageout |> share_recent())
        set!(messagein(interface), messageout(connectedvar(interface), connectedvarindex(interface)))
    end
end

function setmarginal!(factornode::FactorNode, cname::Symbol, marginal)
    lindex = findnext(lmarginal -> name(lmarginal) === cname, localmarginals(factornode), 1)
    @assert lindex !== nothing "Invalid local marginal id: $cname"
    lmarginal = @inbounds localmarginals(factornode)[ lindex ]
    setmarginal!(getstream(lmarginal), marginal)
end

function getmarginal!(factornode::FactorNode, localmarginal::FactorNodeLocalMarginal)
    return getmarginal!(factornode, localmarginal, IncludeAll())
end

function getmarginal!(factornode::FactorNode, localmarginal::FactorNodeLocalMarginal, skip_strategy::MarginalSkipStrategy)
    cached_stream = getstream(localmarginal)

    if cached_stream !== nothing
        return as_marginal_observable(cached_stream, skip_strategy)
    end

    clusterindex = index(localmarginal)

    marginalname = name(localmarginal)
    marginalsize = @inbounds length(factorisation(factornode)[ clusterindex ])

    if marginalsize === 1
        # Cluster contains only one variable, we can take marginal over this variable
        vmarginal = getmarginal(connectedvar(getinterface(factornode, marginalname)), IncludeAll())
        setstream!(localmarginal, vmarginal)
        return as_marginal_observable(vmarginal, skip_strategy)
    else
        cmarginal = MarginalObservable()

        message_dependencies  = tuple(getclusterinterfaces(factornode, clusterindex)...)
        marginal_dependencies = tuple(TupleTools.deleteat(localmarginals(factornode), clusterindex)...)

        msgs_names, msgs_observable          = get_messages_observable(factornode, message_dependencies)
        marginal_names, marginals_observable = get_marginals_observable(factornode, marginal_dependencies)

        fform       = functionalform(factornode)
        vtag        = Val{ name(localmarginal) }
        meta        = metadata(factornode)

        mapping = let fform = fform, vtag = vtag, msgs_names = msgs_names, marginal_names = marginal_names, meta = meta, factornode = factornode
            (dependencies) -> begin
                messages  = dependencies[1]
                marginals = getrecent(dependencies[2])

                # Marginal is clamped if all of the inputs are clamped
                is_marginal_clamped = __check_all(is_clamped, messages) && __check_all(is_clamped, marginals)

                # Marginal is initial if it is not clamped and all of the inputs are either clamped or initial
                is_marginal_initial = !is_marginal_clamped && (__check_all(m -> is_clamped(m) || is_initial(m), messages) && __check_all(m -> is_clamped(m) || is_initial(m), marginals))

                return Marginal(marginalrule(fform, vtag, msgs_names, messages, marginal_names, marginals, meta, factornode), is_marginal_clamped, is_marginal_initial)
            end
        end

        # TODO: discontinue operater is needed for loopy belief propagation? Check
        marginalout = combineLatest((msgs_observable, marginals_observable), PushNew()) |> discontinue() |> map(Marginal, mapping)

        connect!(cmarginal, marginalout) # MarginalObservable has RecentSubject by default, there is no need to share_recent() here

        setstream!(localmarginal, cmarginal)

        return as_marginal_observable(cmarginal, skip_strategy)
    end
end

## TODO
function conjugate_type end

## make_node

struct AutoVar
    name :: Symbol
end

getname(autovar::AutoVar) = autovar.name

function make_node end

function interface_get_index end
function interface_get_name end

function interface_get_index(::Type{ Val{ Node } }, ::Type{ Val{ Interface } }) where { Node, Interface }
    error("Node $Node has no interface named $Interface")
end

function interface_get_name(::Type{ Val{ Node } }, ::Type{ Val{ Interface } }) where { Node, Interface }
    error("Node $Node has no interface named $Interface")
end

make_node(fform, ::AutoVar, ::Vararg{ <: AbstractVariable }; kwargs...) = error("Unknown functional form '$(fform)' used for node specification.")
make_node(fform, args::Vararg{ <: AbstractVariable }; kwargs...)        = error("Unknown functional form '$(fform)' used for node specification.")

function make_node(fform::Function, autovar::AutoVar, args::Vararg{ <: ConstVariable }; kwargs...)
    var  = constvar(getname(autovar), fform(map((d) -> getconst(d), args)...))
    return nothing, var
end

function make_node(::Type{ T }, autovar::AutoVar, args::Vararg{ <: ConstVariable }; kwargs...) where T
    var  = constvar(getname(autovar), T(map((d) -> getconst(d), args)...))
    return nothing, var
end

function make_node(fform::Function, autovar::AutoVar, args::Vararg{ <: DataVariable{ <: PointMass } }; kwargs...)
    # TODO
    message_cb = let fform = fform
        (d::Tuple) -> Message(fform(d...), false, false)
    end

    subject = combineLatest(tuple(map((a) -> messageout(a, getlastindex(a)) |> map(Any, (d) -> mean(getdata(d))), args)...), PushNew()) |> map(Message, message_cb)
    var     = datavar(getname(autovar), Any, subject = subject)

    return nothing, var
end

# end

## macro helpers

import .MacroHelpers

"""
    @node(fformtype, sdtype, interfaces_list)

`@node` macro creates a node for a `fformtype` type object

# Examples
```julia

struct MyNormalDistribution
    mean :: Float64
    var  :: Float64
end

@node MyNormalDistribution Stochastic [ out, mean, var ]
```

```julia

@node typeof(+) Deterministic [ out, in1, in2 ]
```

See also: [`make_node`](@ref), [`Stochastic`](@ref), [`Deterministic`](@ref)
"""
macro node(fformtype, sdtype, interfaces_list)

    fbottomtype = MacroHelpers.bottom_type(fformtype)
    fuppertype  = MacroHelpers.upper_type(fformtype)

    @assert sdtype ∈ [ :Stochastic, :Deterministic ] "Invalid sdtype $(sdtype). Can be either Stochastic or Deterministic."

    @capture(interfaces_list, [ interfaces_args__ ]) ||
        error("Invalid interfaces specification.")

    interfaces = map(interfaces_args) do arg
        if @capture(arg, name_Symbol)
            return (name, [])
        elseif @capture(arg, (name_Symbol, aliases = [ aliases__ ]))
            @assert all(a -> a isa Symbol && !isequal(a, name), aliases)
            return (name, aliases)
        else
            error("Interface specification should have a 'name' or (name, aliases = [ alias1, alias2,... ]) signature.")
        end
    end

    @assert length(interfaces) !== 0 "Node should have at least one interface."

    names = map(d -> first(d), interfaces)

    names_quoted_tuple     = Expr(:tuple, map(name -> Expr(:quote, name), names)...)
    names_indices          = Expr(:tuple, map(i -> i, 1:length(names))...)
    names_splitted_indices = Expr(:tuple, map(i -> Expr(:tuple, i), 1:length(names))...)

    interface_args        = map(name -> :($name::AbstractVariable), names)
    interface_connections = map(name -> :(ReactiveMP.connect!(node, $(Expr(:quote, name)), $name)), names)

    # Here we create helpers function for GraphPPL.jl interfacing
    # They are used to convert interface names from `where { q = q(x, y)q(z) }` to an equivalent tuple respresentation, e.g. `((1, 2), (3, ))`
    # The general recipe to get a proper index is to call `interface_get_index(Val{ :NodeTypeName }, interface_get_name(Val{ :NodeTypeName }, Val{ :name_expression }))`
    interface_name_getters = map(enumerate(interfaces)) do (index, interface)
        name    = first(interface)
        aliases = last(interface)

        index_name_getter  = :(ReactiveMP.interface_get_index(::Type{ Val{ $(Expr(:quote, fbottomtype)) } }, ::Type{ Val{ $(Expr(:quote, name)) } }) = $(index))
        name_symbol_getter = :(ReactiveMP.interface_get_name(::Type{ Val{ $(Expr(:quote, fbottomtype)) } }, ::Type{ Val{ $(Expr(:quote, name)) } }) = $(Expr(:quote, name)))
        name_index_getter  = :(ReactiveMP.interface_get_name(::Type{ Val{ $(Expr(:quote, fbottomtype)) } }, ::Type{ Val{ $index } }) = $(Expr(:quote, name)))

        alias_getters = map(aliases) do alias
            return :(ReactiveMP.interface_get_name(::Type{ Val{ $(Expr(:quote, fbottomtype)) } }, ::Type{ Val{ $(Expr(:quote, alias)) } }) = $(Expr(:quote, name)))
        end

        return quote
            $index_name_getter
            $name_symbol_getter
            $name_index_getter
            $(alias_getters...)
        end
    end

    # By default every argument passed to a factorisation option of the node is transformed by
    # `collect_factorisation` function to have a tuple like structure.
    # The default recipe is simple: for stochastic nodes we convert `FullFactorisation` and `MeanField` objects
    # to their tuple of indices equivalents. For deterministic nodes any factorisation is replaced by a FullFactorisation equivalent
    factorisation_collectors = if sdtype === :Stochastic
        quote
            ReactiveMP.collect_factorisation(::$fuppertype, factorisation::Tuple)           = factorisation
            ReactiveMP.collect_factorisation(::$fuppertype, ::ReactiveMP.FullFactorisation) = ($names_indices, )
            ReactiveMP.collect_factorisation(::$fuppertype, ::ReactiveMP.MeanField)         = $names_splitted_indices
        end

    elseif sdtype === :Deterministic
        quote
            ReactiveMP.collect_factorisation(::$fuppertype, factorisation::Tuple)           = ($names_indices, )
            ReactiveMP.collect_factorisation(::$fuppertype, ::ReactiveMP.FullFactorisation) = ($names_indices, )
            ReactiveMP.collect_factorisation(::$fuppertype, ::ReactiveMP.MeanField)         = ($names_indices, )
        end
    else
        error("Unreachable in @node macro.")
    end

    make_node_const_mapping = if sdtype === :Stochastic
        quote
            function ReactiveMP.make_node(fform::$fuppertype, autovar::ReactiveMP.AutoVar, args::Vararg{ <: ReactiveMP.ConstVariable{ <: ReactiveMP.PointMass } }; kwargs...)
                var  = ReactiveMP.randomvar(ReactiveMP.getname(autovar))
                node = ReactiveMP.make_node(fform, var, args...; kwargs...)
                return node, var
            end
        end
    elseif sdtype === :Deterministic
        quote
            function ReactiveMP.make_node(fform::$fuppertype, autovar::ReactiveMP.AutoVar, args::Vararg{ <: ReactiveMP.ConstVariable{ <: ReactiveMP.PointMass } }; kwargs...)
                var  = ReactiveMP.constvar(ReactiveMP.getname(autovar), fform(map((d) -> ReactiveMP.getconst(d), args)...))
                return nothing, var
            end
        end
    else
        error("Unreachable in @node macro.")
    end

    res = quote

        ReactiveMP.as_node_functional_form(::$fuppertype) = ReactiveMP.ValidNodeFunctionalForm()

        ReactiveMP.sdtype(::$fuppertype) = (ReactiveMP.$sdtype)()

        function ReactiveMP.make_node(::$fuppertype; factorisation = ($names_indices, ), meta = nothing, pipeline = nothing)
            return ReactiveMP.FactorNode($fbottomtype, $names_quoted_tuple, ReactiveMP.collect_factorisation($fbottomtype, factorisation), ReactiveMP.collect_meta($fbottomtype, meta), ReactiveMP.collect_pipeline($fbottomtype, pipeline))
        end

        function ReactiveMP.make_node(::$fuppertype, $(interface_args...); factorisation = ($names_indices, ), meta = nothing, pipeline = nothing)
            node = ReactiveMP.make_node($fbottomtype, factorisation = factorisation, meta = meta, pipeline = pipeline)
            $(interface_connections...)
            return node
        end

        function ReactiveMP.make_node(fform::$fuppertype, autovar::ReactiveMP.AutoVar, args::Vararg{ <: ReactiveMP.AbstractVariable }; kwargs...)
            var  = ReactiveMP.randomvar(ReactiveMP.getname(autovar))
            node = ReactiveMP.make_node(fform, var, args...; kwargs...)
            return node, var
        end

        $(make_node_const_mapping)
        $(interface_name_getters...)

        $factorisation_collectors

    end

    return esc(res)
end


text:
@testset "Constant" begin
    # Instantiation
    N1 = 10
    A1 = rand(Float64, N1) * 100
    A2 = rand(Float64, N1, N1) * 100
    A3 = rand(Float64, N1, N1, N1) * 100

    @testset "constructor" begin
        @test Constant() == Constant{Nearest}()
        for T in (Nearest, Previous, Next)
            it = Constant{T}()
            @test it isa Constant{T}
        end
    end

    @testset "Constant{Nearest}" begin
        for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
            isinplace = constructor == interpolate!
            itp1 = @inferred(constructor(copier(A1), BSpline(Constant())))
            itp2 = @inferred(constructor(copier(A2), BSpline(Constant())))
            itp3 = @inferred(constructor(copier(A3), BSpline(Constant())))

            @test parent(itp1) === itp1.coefs
            @test Interpolations.lbounds(itp1) == (1,)
            @test Interpolations.ubounds(itp1) == (N1,)

            # Evaluation on provided data points
            for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
                check_axes(itp, A, isinplace)
                check_inbounds_values(itp, A)
                check_oob(itp)
                can_eval_near_boundaries(itp)
            end

            # Evaluation between data points (tests constancy)
            for i in 2:N1 - 1
                @test A1[i] == itp1(i + .3) == itp1(i - .3)
            end
            # 2D
            for i in 2:N1 - 1, j in 2:N1 - 1
                @test A2[i,j] == itp2(i + .4, j - .3) == itp2(i - .4, j + .3)
            end
            # 3D
            for i in 2:N1 - 1, j in 2:N1 - 1, k in 2:N1 - 1
                @test A3[i,j,k] == itp3(i + .4, j - .3, k + .1) == itp3(i + .4, j - .3, k + .2)
            end
        end
    end

    @testset "Constant{Previous}" begin
        for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
            isinplace = constructor == interpolate!
            itp1 = @inferred(constructor(copier(A1), BSpline(Constant{Previous}())))
            itp2 = @inferred(constructor(copier(A2), BSpline(Constant{Previous}())))
            itp3 = @inferred(constructor(copier(A3), BSpline(Constant{Previous}())))

            @test parent(itp1) === itp1.coefs
            @test Interpolations.lbounds(itp1) == (1,)
            @test Interpolations.ubounds(itp1) == (N1,)

            # Evaluation on provided data points
            for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
                check_axes(itp, A, isinplace)
                check_inbounds_values(itp, A)
                check_oob(itp)
                can_eval_near_boundaries(itp)
            end

            # Evaluation between data points (tests constancy)
            for i in 2:N1 - 1
                @test A1[i] == itp1(i + .3) == itp1(i + .6)
                @test A1[i - 1] == itp1(i - .3) == itp1(i - .6)
            end
            # 2D
            for i in 2:N1 - 1, j in 2:N1 - 1
                @test A2[i,j - 1] == itp2(i + .4, j - .3) == itp2(i + .4, j - .6)
            end
            # 3D
            for i in 2:N1 - 1, j in 2:N1 - 1, k in 2:N1 - 1
                @test A3[i,j - 1,k] == itp3(i + .4, j - .3, k + .1) == itp3(i + .4, j - .3, k + .2)
            end
        end
    end

    @testset "Constant{Next}" begin
        for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
            isinplace = constructor == interpolate!
            itp1 = @inferred(constructor(copier(A1), BSpline(Constant{Next}())))
            itp2 = @inferred(constructor(copier(A2), BSpline(Constant{Next}())))
            itp3 = @inferred(constructor(copier(A3), BSpline(Constant{Next}())))

            @test parent(itp1) === itp1.coefs
            @test Interpolations.lbounds(itp1) == (1,)
            @test Interpolations.ubounds(itp1) == (N1,)

            # Evaluation on provided data points
            for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
                check_axes(itp, A, isinplace)
                check_inbounds_values(itp, A)
                check_oob(itp)
                can_eval_near_boundaries(itp)
            end

            # Evaluation between data points (tests constancy)
            for i in 2:N1 - 1
                @test A1[i + 1] == itp1(i + .3) == itp1(i + .6)
                @test A1[i] == itp1(i - .3) == itp1(i - .6)
            end
            # 2D
            for i in 2:N1 - 1, j in 2:N1 - 1
                @test A2[i + 1,j] == itp2(i + .4, j - .3) == itp2(i + .4, j - .6)
            end
            # 3D
            for i in 2:N1 - 1, j in 2:N1 - 1, k in 2:N1 - 1
                @test A3[i + 1,j,k + 1] == itp3(i + .4, j - .3, k + .1) == itp3(i + .4, j - .3, k + .2)
            end
        end
    end

    @testset "Constant periodic" begin
        # Constructors
        @test Constant() === Constant(Throw(OnGrid())) === Constant{Nearest,Throw{OnGrid}}()
        @test Constant() isa Constant{Nearest,Throw{OnGrid}}
        @test Constant(Periodic()) === Constant(Periodic(OnCell())) === Constant{Nearest,Periodic{OnCell}}()
        @test Constant(Periodic()) isa Constant{Nearest,Periodic{OnCell}}
        for T in (Nearest, Previous, Next)
            it = Constant{T}()
            @test it isa Constant{T, Throw{OnGrid}}
            @test "$it" == "Constant{$T}()"
            it = Constant{T}(Periodic())
            @test it isa Constant{T, Periodic{OnCell}}
            @test "$it" == "Constant{$T}(Periodic(OnCell()))"
        end

        for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
            isinplace = constructor == interpolate!
            itp_periodic = @inferred(constructor(copier(A1), BSpline(Constant(Periodic()))))
            itp_previous = @inferred(constructor(copier(A1), BSpline(Constant{Previous}(Periodic()))))
            itp_next = @inferred(constructor(copier(A1), BSpline(Constant{Next}(Periodic()))))

            for itp in (itp_periodic, itp_previous, itp_next)
                @test parent(itp) === itp.coefs
                @test all(Interpolations.lbounds(itp) .≈ (0.5,))
                @test all(Interpolations.ubounds(itp) .≈ (N1 + 0.5,))

                check_axes(itp, A1, isinplace)
                check_inbounds_values(itp, A1)
                check_oob(itp)
                can_eval_near_boundaries(itp)
            end

            # Evaluation between data points (tests constancy)
            for i in 2:N1-1
                @test A1[i] == itp_periodic(i + .3) == itp_periodic(i - .3)
                @test A1[i] == itp_previous(i + .3) == itp_previous(i + .6)
                @test A1[i - 1] == itp_previous(i - .3) == itp_previous(i - .6)
                @test A1[i + 1] == itp_next(i + .3) == itp_next(i + .6)
                @test A1[i] == itp_next(i - .3) == itp_next(i - .6)
            end

            # Evaluation between data points in [0.5, 1.5], [N1 - 0.5, N1 + 0.5].
            @test A1[1] == itp_periodic(1 - .3) == itp_periodic(1 + .3)
            @test A1[N1] == itp_periodic(N1 - .3) == itp_periodic(N1 + .3)
            @test A1[1] == itp_previous(1 + .3)
            @test A1[N1] == itp_previous(N1 + .3) == itp_previous(1 - .3)
            @test A1[1] == itp_next(1 - .3) == itp_next(N1 + .3)
            @test A1[N1] == itp_next(N1 - .3)
        end
    end
end

ext: jl


text:
module ComposeTests

using Gridap.TensorValues
using Gridap.Arrays
using Gridap.Fields
using Gridap.Fields: MockField
using Gridap.Fields: Comp
using FillArrays

import Gridap.Fields: ∇

np = 4
p = Point(1,2)
x = fill(p,np)

fun(x,y) = 2*x
∇fun(x,y) = VectorValue(2*one(x[1]),2*one(x[1]))
∇(::typeof(fun)) = ∇fun

v = 3.0
d = 2
f = MockField{d}(v)
fx = fill(v,np)

k = Comp(fun)
test_kernel(k,(v,v),2*v)

g = compose(fun,f,f)
gx = 2*fx
∇gx = fill(∇fun(v,v),np)
test_field(g,x,gx,grad=∇gx)

l = 10
af = Fill(f,l)
ax = fill(x,l)
ag = compose(fun,af,af)
agx = fill(gx,l)
a∇gx = fill(∇gx,l)
test_array_of_fields(ag,ax,agx,grad=a∇gx)

np = 4
p = Point(1,2)
x = fill(p,np)

v = VectorValue(1,2,3)
d = 2
f = MockField{d}(v)
w = 4
d = 3
g = MockField{d}(w)

fx = evaluate(f,x)
gx = evaluate(g,fx)
∇gx = evaluate(∇(g),fx)

h = compose_fields(g,f)
test_field(h,x,gx,grad=∇gx)

h = compose(g,f)
test_field(h,x,gx,grad=∇gx)

l = 10
af = Fill(f,l)
ag = Fill(g,l)
ax = Fill(x,l)

afx = evaluate(af,ax)
agx = evaluate(ag,afx)
∇agx = evaluate(∇(ag),afx)

ah = compose_field_arrays(ag,af)
test_array_of_fields(ah,ax,agx,grad=∇agx)

ah = compose(ag,af)
test_array_of_fields(ah,ax,agx,grad=∇agx)

end # module


text:
using Documenter, JqData

makedocs(
    modules = [JqData],
    sitename="JqData.jl",
    pages=Any[
        "Home" => "index.md",
    ],
)

push!(LOAD_PATH,"../src")

deploydocs(
    repo = "github.com/hzgzh/JqData.jl",
    branch = "gh-pages",
)


text:
"""
    NamedVar{T}

Defines a named variable

# Field
 - `name::Symbol`: Variable name
 - `value::T`: Default variable value
"""
struct NamedVar{T}
    name::Symbol
    value::T
end

"""
    @var var = val

Macro which modifies the expression `var = val` to `var = NamedVar(:var, val)`
```julia-repl
julia> @var x = 5
  NamedVar{Int64}(:x, 5)

julia> x
  NamedVar{Int64}(:x, 5)
```
"""
macro var(ex)
    @assert (ex.head === :(=)) "Incorrect usage of @var macro"
    var = ex.args[1]
    name = string(var)
    val = ex.args[2]
    return esc(:($var = NamedVar(Symbol($name), $(val))))
end

"""
    name(var::NamedVar)

Return the name of a named variable.
"""
name(var::NamedVar) = var.name

"""
    value(var::NamedVar)

Return the value of a named variable.
"""
value(var::NamedVar) = var.value

"""
    value(val)

Return the value of a variable.
"""
value(val) = val

"""
    length(var::NamedVar)

Return the length of a named variable
"""
Base.length(var::NamedVar) = length(var.value)

"""
    size(var)

Return the size of a named variable
"""
Base.size(var::NamedVar) = size(var.value)

"""
    axes(var::NamedVar)

Return the axes of a named variable
"""
Base.axes(var::NamedVar) = axes(var.value)

"""
    combine(vars::Tuple)

Combine the values of all variables in `vars` into a single vector
"""
combine(vars) = vcat(vectorize.(value.(vars))...)

"""
    combine!(y, vars::Tuple)

In-place version of [`combine`](@ref)
"""
@inline function combine!(y, vars::Tuple, idx=0)
    var = first(vars)
    nval = length(var)
    y[idx+1:idx+nval] .= var
    idx += nval
    return combine!(y, tail(vars), idx)
end
@inline combine!(y, vars::Tuple{}, idx=0) = y

"""
    separate(vars::Tuple, values)

Return the values stored in `values` as separated variables rather than a single
output vector.  `vars` is a tuple of named variables corresponding to each
returned parameter.
"""
@inline function separate(vars::Tuple, values)
    # extract first variable
    n1 = length(first(vars))
    v1 = same_axes(value(first(vars)), view(values, 1:n1))
    # call recursively to extract remaining variables
    vars = tail(vars)
    values = view(values, n1+1:length(values))
    return (v1, separate(vars, values)...)
end
@inline separate(vars::Tuple{}, values) = ()

"""
    vectorize(x)

Converts `x` to a vector.
"""
vectorize
@inline vectorize(x::Number) = [x]
@inline vectorize(x::AbstractArray) = x[:]

"""
    same_axes(v1, v2)

Return `v2` with axes correponding to `v1`.  If `v1` is a scalar, return a scalar.
"""
same_axes
@inline same_axes(v1::AbstractArray, v2) = reshape(v2, axes(v1))
@inline same_axes(v1::Number, v2) = reshape(v2, axes(v1))[]

"""
    getindices(collection, indices::Tuple)

Custom version of `getindex` that allow collections to be accessed using a tuple
of indices
"""
@inline getindices(collection, indices::Tuple) = (collection[first(indices)],
    getindices(collection, tail(indices))...)
@inline getindices(collection, indices::Tuple{}) = ()

"""
    setindices(x, v::Tuple, indices::Tuple)

Creates a new collection similar to `x` with the values at the indices in `indices` set
to the values in `v`.
"""
@inline setindices(x, v::Tuple, indices::Tuple) = setindices(setindex(x, first(v), first(indices)), tail(v), tail(indices))
@inline setindices(x, v::Tuple{}, indices::Tuple{}) = x

"""
    system_component_mapping(components, argin)

Constructs the component input-output mapping for a system based on input and output
argument names.
"""
function system_component_mapping(components, argin)

    # names of input and output variables
    input_names = [name.(comp.argin) for comp in components]
    output_names = [name.(comp.argout) for comp in components]

    # construct component input mapping vector
    component_input_mapping = Vector{Vector{NTuple{2,Int}}}(undef, length(components))
    for icomp = 1:length(components)
        # initialize mapping for this component
        current_mapping = Vector{NTuple{2,Int}}(undef, length(inputs(components[icomp])))
        # input variable starting index
        idx_in = 0
        # populate mapping for this component
        for i = 1:length(input_names[icomp])
            # number of values in variable
            nval = length(components[icomp].argin[i])
            # first check system inputs for match
            idx_sys = 0
            for j = 1:length(argin)
                if input_names[icomp][i] === name(argin[j])
                    # add to mapping
                    for k = 1:nval
                        current_mapping[idx_in+k] = (0, idx_sys+k)
                    end
                    # skip to next input variable
                    @goto next_input_variable
                else
                    # check next system input
                    idx_sys += length(argin[j])
                end
            end
            # then check all components for match
            for jcomp = 1:length(components)
                # only check components other than self
                if icomp == jcomp
                    continue
                end
                # output variable starting index
                idx_out = 0
                # loop through all output variables
                for j = 1:length(output_names[jcomp])
                    # check if names match
                    if input_names[icomp][i] === output_names[jcomp][j]
                        # add to mapping
                        for k = 1:nval
                            current_mapping[idx_in+k] = (jcomp, idx_out+k)
                        end
                        # skip to next input variable
                        @goto next_input_variable
                    else
                        # check next output variable
                        idx_out += length(components[jcomp].argout[j])
                    end
                end
            end
            error("No system input or component found for input `$(input_names[icomp][i])`")
            @label next_input_variable
            idx_in += nval
        end
        component_input_mapping[icomp] = current_mapping
    end

    return component_input_mapping
end

"""
    system_output_mapping(components, argin, argout)

Constructs the output mapping for a system based on input and output argument
names of the system and components.
"""
function system_output_mapping(components, argin, argout)

    # system output vector
    ysys = combine(argout)

    # names of input and output variables
    input_names = [name.(comp.argin) for comp in components]
    output_names = [name.(comp.argout) for comp in components]

    # construct system output mapping vector
    output_mapping = Vector{NTuple{2,Int}}(undef, length(ysys))
    # output variable starting index
    idx_sys = 0
    # populate mapping for this component
    for i = 1:length(argout)
        # number of values in variable
        nval = length(argout[i])
        # first check system inputs for match
        idx_in = 0
        for j = 1:length(argin)
            if name(argout[i]) === name(argin[j])
                # add to mapping
                for k = 1:nval
                    output_mapping[idx_sys+k] = (0, idx_in+k)
                end
                # skip to next input variable
                @goto next_output_variable
            else
                # check next system input
                idx_in += length(argin[j])
            end
        end
        # then check all components for match
        for jcomp = 1:length(components)
            # output variable starting index
            idx_out = 0
            # loop through all output variables
            for j = 1:length(output_names[jcomp])
                # check if names match
                if name(argout[i]) === output_names[jcomp][j]
                    # add to mapping
                    for k = 1:nval
                        output_mapping[idx_sys+k] = (jcomp, idx_out+k)
                    end
                    # skip to next input variable
                    @goto next_output_variable
                else
                    # check next output variable
                    idx_out += length(components[jcomp].argout[j])
                end
            end
        end
        error("No system input or component output found for output `$(name(argout[i]))`")
        @label next_output_variable
        idx_sys += nval
    end

    return output_mapping
end

"""
    output_mapping_matrices(argin, argstate, argout)

Construct matrices that maps an implicit function's inputs and/or state variables
to specified output variables.
"""
function output_mapping_matrices(argin, argstate, argout)

    # length of input, state, and output vectors
    nx = sum(length.(argin))
    nu = sum(length.(argstate))
    ny = sum(length.(argout))

    # allocate matrices
    dfdx = BitArray(undef, ny, nx)
    dfdu = BitArray(undef, ny, nu)

    # initialize matrices
    dfdx .= false
    dfdu .= false

    # output variable starting index
    idx_out = 0
    # populate mapping for this component
    for i = 1:length(argout)
        # number of values in variable
        nval = length(argout[i])
        # first check system inputs for match
        idx_in = 0
        for j = 1:length(argin)
            if name(argout[i]) === name(argin[j])
                # add to mapping
                for k = 1:nval
                    dfdx[idx_out + k, idx_in + k] = 1
                end
                # skip to next input variable
                @goto next_output_variable
            else
                # check next system input
                idx_in += length(argin[j])
            end
        end
        # then check system states for match
        idx_state = 0
        for j = 1:length(argstate)
            if name(argout[i]) === name(argstate[j])
                # add to mapping
                for k = 1:nval
                    dfdu[idx_out + k, idx_state + k] = 1
                end
                # skip to next input variable
                @goto next_output_variable
            else
                # check next system input
                idx_state += length(argstate[j])
            end
        end
        error("No input or state variable found for output `$(name(argout[i]))`")
        @label next_output_variable
        idx_out += nval
    end

    return dfdx, dfdu
end


text:
using ApproxFun, Base.Test

f=Fun(exp,Jacobi(2.,.5))
@test_approx_eq f(.1) exp(.1)

f=Fun(x->cos(100x),Jacobi(2.124,.5),500)
@test_approx_eq f(.1) cos(100*.1)

sp=Jacobi(2.124,.5)
f=Fun(exp,sp)
sp2=Jacobi(2.124,1.5)
f2=Fun(exp,sp2)
sp3=Jacobi(3.124,1.5)
f3=Fun(exp,sp3)
sp4=Jacobi(4.124,2.5)
f4=Fun(exp,sp4)
@test norm((Fun(f,sp2)-f2).coefficients)<10eps()
@test norm((Fun(f,sp3)-f3).coefficients)<10eps()
@test norm((Fun(f,sp4)-f4).coefficients)<20eps()

m=20
f=Fun(x->((1-x)/2).^m.*exp(x),JacobiWeight(0.,m,Jacobi(2m+1,0.)))
@test abs(f(.1)-(x->((1-x)/2).^m.*exp(x))(.1))<10eps()

m=10
f=Fun(x->besselj(m,m*(1-x)),JacobiWeight(0.,m,Jacobi(2m+1,0.)))
@test_approx_eq f(0.) besselj(m,m)

## Conversion

@test norm(Fun(Fun(exp),Jacobi(-.5,-.5))-Fun(exp,Jacobi(-.5,-.5))) < 100eps()

x=Fun(identity)
ri=0.5./(1-x)
@test_approx_eq ((1-x)./2.*Fun(exp,JacobiWeight(0.,0.,Jacobi(1.,0.))))(.1) (1-.1)./2*exp(.1)

@test_approx_eq ((1-x)./2.*Fun(exp,JacobiWeight(0.,0.,Jacobi(1.,0.))))(.1) (1-.1)./2*exp(.1)

@test_approx_eq (ri.*Fun(exp,JacobiWeight(0.,0.,Jacobi(1.,0.))))(.1) .5/(1-.1)*exp(.1)

## Derivative

S=JacobiWeight(0.,0.,Jacobi(1.,0.,Interval(1.,0.)))
D=Derivative(S)
f=Fun(exp,domainspace(D))
@test (D*f-f).coefficients|>norm < eps(100000.)
@test (f'-f).coefficients|>norm < eps(100000.)
@test (D^2*f-f).coefficients|>norm < eps(100000000.)
@test (D*(D*f)-f).coefficients|>norm < eps(100000000.)

### Jacobi multiplication

x=Fun(identity,Jacobi(0.,0.))
f=Fun(exp,Jacobi(0.,0.))

@test_approx_eq (x*f)(.1) .1exp(.1)

x=Fun(identity,Jacobi(0.123,12.324))
f=Fun(exp,Jacobi(0.,0.))

@test_approx_eq (x*f)(.1) .1exp(.1)

x=Fun(identity,Jacobi(0.123,12.324))
f=Fun(exp,Jacobi(0.213,0.590))

@test_approx_eq (x*f)(.1) .1exp(.1)

g=Fun(cos,Jacobi(0.123,12.324))
f=Fun(exp,Jacobi(0.213,0.590))

@test_approx_eq (g*f)(.1) cos(.1)*exp(.1)

## Jacobi integrate and sum

@test_approx_eq sum(Fun(exp,Legendre([0,2]))) sum(Fun(exp,[0,2]))

a=Arc(0.,.1,0.,π/2)
g=Fun(exp,Legendre(a))

@test_approx_eq sum(g) sum(Fun(exp,a))

## Test special derivative

x=Fun()
f=exp(x)*sqrt(1-x^2)
D=Derivative(WeightedJacobi(.5,.5))
g=(D*Fun(f,domainspace(D)))
@test_approx_eq f'(0.1) g(0.1)

## Test implementation of conversion between Chebyshev and Jacobi spaces using FastTransforms

f = Fun(x->cospi(1000x))
g = Fun(f,Legendre())
h = Fun(g,Chebyshev())
@test norm(f.coefficients-h.coefficients,Inf) < 10eps()
h = Fun(h,Legendre())
@test norm(g.coefficients-h.coefficients,Inf) < 1000eps()

## Legendre conversions

@test norm(Fun(exp,Ultraspherical(1//2))-Fun(exp,Jacobi(0,0))) < 100eps()

C=Conversion(Jacobi(0,0),Chebyshev())
@test norm(C*Fun(exp,Jacobi(0,0))  - Fun(exp)) < 100eps()

C=Conversion(Ultraspherical(1//2),Chebyshev())
@test norm(C*Fun(exp,Ultraspherical(1//2))  - Fun(exp)) < 100eps()

C=Conversion(Chebyshev(),Ultraspherical(1//2))
@test norm(C*Fun(exp)-Fun(exp,Legendre())) < 100eps()

C=Conversion(Chebyshev(),Jacobi(0,0))
@test norm(C*Fun(exp)  - Fun(exp,Jacobi(0,0))) < 100eps()

C=Conversion(Chebyshev(),Jacobi(1,1))
@test norm(C*Fun(exp) - Fun(exp,Jacobi(1,1))) < 100eps()

C=Conversion(Ultraspherical(1//2),Ultraspherical(1))
@test norm(C*Fun(exp,Ultraspherical(1//2))-Fun(exp,Ultraspherical(1))) < 100eps()

C=Conversion(Jacobi(0,0),Ultraspherical(1))
@test norm(C*Fun(exp,Jacobi(0,0))-Fun(exp,Ultraspherical(1))) < 100eps()

C=Conversion(Ultraspherical(1),Jacobi(0,0))
@test norm(C*Fun(exp,Ultraspherical(1))-Fun(exp,Jacobi(0,0))) < 100eps()


text:
function ueberpruefef_nah(f_nah,A2,Messwerte,Aufloesung,Schritt)

Phantom = echteCTDaten(Aufloesung);
nrows,ncols = size(Phantom);
f_nahBild = zeros(size(Phantom,1),size(Phantom,2));

if Schritt == 1.1
    if almostequal(f_nah,(transpose(A2)*A2)\(transpose(A2)*Messwerte))
        ObjektNaiv = reshape(f_nah,(nrows,ncols));
        println("Eure Formel für das nächste Objekt am Original ist korrekt!");
        f_nahBild = convert(Array{Gray{Normed{UInt8,8}},2},ObjektNaiv);
    else
        println("Eure Formel für f_nah ist noch nicht korrekt. Überprüft erneut eure Umformungen und nehmt eine Hilfekarte zur Hilfe.")
    end
else
    if almostequal(f_nah,(transpose(A2)*A2)\(transpose(A2)*Messwerte))
        ObjektNaiv = reshape(f_nah,(nrows,ncols));
        println("Eure Formel für f_nahfehler ist korrekt!");
        f_nahBild = convert(Array{Gray{Normed{UInt8,8}},2},ObjektNaiv);
    else
        println("Eure Formel für f_nahfehler ist nicht korrekt. Überprüft, ob ihr euch beim Kopieren vertan habt.")
    end
end

return f_nah, f_nahBild;
end


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

using Markdown
using InteractiveUtils

# ╔═╡ f32b184a-ba70-4e80-8f8d-e12def369128
using Parameters

# ╔═╡ 456df07b-1be8-4adb-9302-9290314311e3
begin
	@enum  CartesianDirection XDirection = 1 YDirection ZDirection
	@enum  CartesianAxis XAxis = 1 YAxis ZAxis
	@enum  FieldComponent XComponent = 1 YComponent ZComponent
end

# ╔═╡ 4f7a7193-e75e-4756-9423-7d81e2f7e4db
begin

mutable struct StructuredGrid{nD} <: AbstractArray{Real,nD}
	nodes::NTuple{nD,Vector{Float64}}
end

function StructuredGrid(nodes...)
        StructuredGrid{length(nodes)}(nodes)
end

Base.size(grid::StructuredGrid)	= length.(grid.nodes)

nodes(sgrid::StructuredGrid) = sgrid.nodes
xnodes(grid::StructuredGrid)	= grid.nodes[1]
ynodes(grid::StructuredGrid)	= grid.nodes[2]
znodes(grid::StructuredGrid)	= grid.nodes[3]

xnodes(grid::StructuredGrid,i::Int)	= getindex(xnodes(grid),i)
ynodes(grid::StructuredGrid,j::Int)	=  getindex(ynodes(grid),j)
znodes(grid::StructuredGrid,k::Int)	= getindex(znodes(grid),k)

Base.getindex(grid::StructuredGrid{1},i::Int) = (xnodes(grid,i))
Base.getindex(grid::StructuredGrid{2},i::Int,j::Int) = (xnodes(grid,i),ynodes(grid,j))
Base.getindex(grid::StructuredGrid{3},i::Int,j::Int,k::Int) =  (xnodes(grid,i),ynodes(grid,j), znodes(grid,k))

cell_sizes(sgrid::StructuredGrid) = StructuredGrid(diff.(nodes(sgrid)))
cell_volumes(sgrid::StructuredGrid) = map(prod,cell_sizes(sgrid))
_avg(nodes::AbstractVector) = 0.5(nodes[2:end] + nodes[1:end-1])
cell_centroids(sgrid::StructuredGrid) = StructuredGrid(_avg.(nodes(sgrid)))

abstract type AbstractNodeType end
struct 	Primal <: AbstractNodeType  end
struct 	Dual <: AbstractNodeType  end

function primal_nodes(ax::FieldComponent, reference_grid::StructuredGrid{Dim}) where {Dim}
    StructuredGrid(ntuple(i -> i == Int(ax) ? _avg(nodes(reference_grid)[i]) : nodes(reference_grid)[i], Dim))
end
function dual_nodes(ax::FieldComponent, reference_grid::StructuredGrid{Dim}) where {Dim}
    StructuredGrid(ntuple(i -> i == Int(ax) ? nodes(reference_grid)[i] :  _avg(nodes(reference_grid)[i])  , Dim))
end

function nodes(::Type{Primal},ax::FieldComponent, reference_grid::StructuredGrid)
primal_nodes(ax,reference_grid)
end
function nodes(::Type{Dual},ax::FieldComponent, reference_grid::StructuredGrid)
dual_nodes(ax,reference_grid)
end

end

# ╔═╡ 1bc97871-d861-47d4-bde9-e2373b1e903a
begin
struct YeeGrid{T,nD}
	x::StructuredGrid{nD}
	y::StructuredGrid{nD}
	z::StructuredGrid{nD}
end

function YeeGrid{T}(refgrid::StructuredGrid{nD}) where {nD,T}
	YeeGrid{T,nD}(
		nodes(T,XComponent,refgrid),
		nodes(T,YComponent,refgrid),
		nodes(T,ZComponent,refgrid))
end
function YeeGrid{T}(nodes...) where T
	refgrid = StructuredGrid(nodes...)
	YeeGrid{T}(refgrid)
end
Base.size(yeegrid::YeeGrid) = (size(yeegrid.x),size(yeegrid.y),size(yeegrid.z))

end

# ╔═╡ 4b3b120f-52b0-47ae-bec2-207ea12ddd6a
begin
	grid_3d_primal = YeeGrid{Primal}(LinRange(-2,2,10),LinRange(-5,5,100),LinRange(-3,3,40))
	grid_3d_dual = YeeGrid{Dual}(LinRange(-2,2,10),LinRange(-5,5,100),LinRange(-3,3,40))
	# nonuniform grid
	grid_2d_primal = YeeGrid{Primal}([1,4,7,8,9,10,15,20],LinRange(-1,1,15))
	grid_2d_dual = YeeGrid{Primal}([1,4,7,8,9,10,15,20],LinRange(-1,1,15))

	grid_1d_dual = YeeGrid{Primal}([1,4,7,8,9,10,15,20])
end

# ╔═╡ 6ad42c82-25c5-41c3-83f1-0ff002edbafe
begin
struct VectorField{PD,T,N}
	x::Array{T,N}
	y::Array{T,N}
	z::Array{T,N}
end
function VectorField{T}(yeegrid::YeeGrid{PD,N}) where {T,PD,N}
	VectorField{PD,T,N}(zeros.(T,size(yeegrid))...)
end
struct TensorField{PD,T,N}
	xx::Array{T,N}
	xy::Array{T,N}
	xz::Array{T,N}
	yx::Array{T,N}
	yy::Array{T,N}
	yz::Array{T,N}
	zx::Array{T,N}
	zy::Array{T,N}
	zz::Array{T,N}
end
function TensorField{T}(yeegrid::YeeGrid{PD,N}) where {T,PD,N}
	TensorField{PD,T,N}(zeros.(T,size(yeegrid))...,zeros.(T,size(yeegrid))...,zeros.(T,size(yeegrid))...)
end
end

# ╔═╡ 4fd9f13e-8279-41e2-b7f5-e9865ed869be
function FaradayStep!(b,e,m,cb;ref_grid)

(nx,ny,nz) = size(ref_grid)

for k = 1:nz-1
	for j = 1:ny-1
		for i = 1:nx-1
b.x[i,j,k] =
	 b.x[i,j,k] + cb.x[i,j,k] * (e.z[i,j,k]-e.z[i,j+1,k]+e.y[i,j,k+1]-e.y[i,j,k]) -m.x[i,j,k]
b.y[i,j,k] =
	 b.y[i,j,k] + cb.y[i,j,k] * (e.x[i,j,k]-e.x[i,j,k+1]+e.z[i+1,j,k]-e.z[i,j,k]) - m.y[i,j,k]
b.z[i,j,k] =
	 b.z[i,j,k] + cb.z[i,j,k] * (e.y[i,j,k]-e.y[i+1,j,k]+e.x[i,j+1,k]-e.x[i,j,k]) -
	 m.z[i,j,k]
		end
	end
end

for k = 1:nz-1
	for j = 1:ny-1
		for i = nx:nx
b.x[i,j,k] =
	 b.x[i,j,k] + cb.x[i,j,k] * (e.z[i,j,k]-e.z[i,j+1,k]+e.y[i,j,k+1]-e.y[i,j,k]) -m.x[i,j,k]
		end
	end
end
for k = 1:nz-1
	for j = ny:ny
		for i = 1:nx-1
b.y[i,j,k] =
	 b.y[i,j,k] + cb.y[i,j,k] * (e.x[i,j,k]-e.x[i,j,k+1]+e.z[i+1,j,k]-e.z[i,j,k]) - m.y[i,j,k]
		end
	end
end
for k = nz:nz
	for j = 1:ny-1
		for i = 1:nx-1
b.z[i,j,k] =
	 b.z[i,j,k] + cb.z[i,j,k] * (e.y[i,j,k]-e.y[i+1,j,k]+e.x[i,j+1,k]-e.x[i,j,k]) -
	 m.z[i,j,k]
		end
	end
end

end

# ╔═╡ fcb5f651-3e20-4181-9e72-9310f823adf0
function AmpereStep!(d,h,J,cb;ref_grid)

(nx,ny,nz) = size(ref_grid)

for k = 2:nz-1
	for j = 2:ny-1
		for i = 2:nx-1
d.x[i,j,k] = d.x[i,j,k] +
cb.x[i,j,k] * (h.z[i,j,k]-h.z[i,j-1,k]-h.y[i,j,k]+h.y[i,j,k-1]) -J.x[i,j,k]
d.y[i,j,k] = d.y[i,j,k] +
cb.y[i,j,k] * (h.x[i,j,k]-h.x[i,j,k-1]-h.z[i,j,k]+h.z[i-1,j,k]) -J.y[i,j,k]
d.z[i,j,k] = d.z[i,j,k] +
cb.z[i,j,k] * (h.y[i,j,k]-h.y[i-1,j,k]-h.x[i,j,k]+h.z[i,j-1,k]) - J.z[i,j,k]
		end
	end
end

for k = 2:nz-1
	for j = 2:ny-1
		for i = 1:1
d.x[i,j,k] = d.x[i,j,k] +
cb.x[i,j,k] * (h.z[i,j,k]-h.z[i,j-1,k]-h.y[i,j,k]+h.y[i,j,k-1]) -J.x[i,j,k]		end
	end
end
for k = 2:nz-1
	for j = 1:1
		for i = 2:nx-1
d.y[i,j,k] = d.y[i,j,k] +
cb.x[i,j,k] * (h.x[i,j,k]-h.x[i,j,k-1]-h.z[i,j,k]+h.z[i-1,j,k]) -J.y[i,j,k]		end
	end
end
for k = 1:1
	for j = 2:ny-1
		for i = 2:nx-1
d.z[i,j,k] = d.z[i,j,k] +
cb.x[i,j,k] * (h.y[i,j,k]-h.y[i-1,j,k]-h.x[i,j,k]+h.z[i,j-1,k]) - J.z[i,j,k]
		end
	end
end
end

# ╔═╡ 7f1e505f-7223-4ae8-a6d4-9ff455553256
function _update_primal_constitutive!(e,d,C,nx,ny,nz)
for k = 2:nz-1
	for j = 2:ny-1
		for i = 2:nx-1
dy =  0.25 * (d.y[i,j,k] + d.y[i+1,j,k] + d.y[i+1,j-1,k] + d.y[i,j-1,k] )
dz =  0.25 * (d.z[i,j,k] + d.z[i+1,j,k] + d.z[i+1,j,k-1] + d.z[i,j,k-1] )
e.x[i,j,k]  = C.xx[i,j,k] * d.x[i,j,k]  + C.xy[i,j,k] * dy  + C.xz[i,j,k] * dz

dx = 0.25*(d.x[i,j,k]+d.x[i,j+1,k]+d.x[i-1,j,k]+d.x[i-1,j+1,k])
dz = 0.25*(d.z[i,j,k] + d.z[i+1,j,k]+ d.z[i+1,j,k-1] + d.z[i,j,k-1])

e.y[i,j,k]  = C.yx[i,j,k] * dx  + C.yy[i,j,k] * d.y[i,j,k]  + C.yz[i,j,k] * dz

dx = 0.25 * (d.x[i,j,k] + d.x[i,j,k+1] + d.x[i-1,j,k+1]+d.x[i-1,j,k])
dy = 0.25 * (d.y[i,j,k] + d.y[i,j,k+1] + d.y[i,j-1,k+1]+d.y[i,j-1,k])

e.z[i,j,k]  = C.zx[i,j,k] * dx  + C.zy[i,j,k] * dy  + C.yz[i,j,k] *   d.z[i,j,k]

		end
	end
end

for k = 2:nz-1
	for j = 2:ny-1
		for i = 1:1

dy =  0.25 * (d.y[i,j,k] + d.y[i+1,j,k] + d.y[i+1,j-1,k] + d.y[i,j-1,k] )
dz =  0.25 * (d.z[i,j,k] + d.z[i+1,j,k] + d.z[i+1,j,k-1] + d.z[i,j,k-1] )
e.x[i,j,k]  = C.xx[i,j,k] * d.x[i,j,k]  + C.xy[i,j,k] * dy  + C.xz[i,j,k] * dz

	end
	end
end
for k = 2:nz-1
	for j = 1:1
		for i = 2:nx-1

dx = 0.25*(d.x[i,j,k]+d.x[i,j+1,k]+d.x[i-1,j,k]+d.x[i-1,j+1,k])
dz = 0.25*(d.z[i,j,k] + d.z[i,j+1,k]+ d.z[i,j+1,k-1] + d.z[i,j,k-1])

e.y[i,j,k]  = C.yx[i,j,k] * dx  + C.yy[i,j,k] * d.y[i,j,k]  + C.yz[i,j,k] * dz
		end
	end
end
for k = 1:1
	for j = 2:ny-1
		for i = 2:nx-1
dx = 0.25 * (d.x[i,j,k] + d.x[i,j,k+1] + d.x[i-1,j,k+1]+d.x[i-1,j,k])
dy = 0.25 * (d.y[i,j,k] + d.y[i,j,k+1] + d.y[i,j-1,k+1]+d.y[i,j-1,k])

e.z[i,j,k]  = C.zx[i,j,k] * dx  + C.zy[i,j,k] * dy  + C.yz[i,j,k] *   d.z[i,j,k]
		end
	end
end

end

# ╔═╡ 8176e6dd-0768-4d9b-a2b7-957236e4c55b
function _update_dual_constitutive!(h,b,C,nx,ny,nz)
	for k = 2:nz-1
	for j = 2:ny-1
		for i = 2:nx-1
by = 0.25 * (b.y[i,j,k] + b.y[i-1,j,k] + b.y[i-1,j+1,k] + b.y[i,j+1,k])
bz = 0.25 * (b.z[i,j,k] + b.z[i,j,k+1] + b.z[i-1,j,k+1] + b.z[i-1,j,k])
h.x[i,j,k]  = C.xx[i,j,k] * b.x[i,j,k]  + C.xy[i,j,k] * by  + C.xz[i,j,k] * bz

bx = 0.25*(b.x[i,j,k]+b.x[i,j-1,k]+b.x[i+1,j-1,k]+b.x[i+1,j,k])
bz = 0.25*(b.z[i,j,k] + b.z[i,j-1,k]+ b.z[i,j-1,k+1] + b.z[i,j,k+1])

h.y[i,j,k]  = C.yx[i,j,k] * bx  + C.yy[i,j,k] * b.y[i,j,k]  + C.yz[i,j,k] * bz

bx = 0.25*(b.x[i,j,k]+b.x[i+1,j,k]+b.x[i+1,j,k-1]+b.x[i,j,k-1])
by = 0.25*(b.y[i,j,k]+b.y[i,j+1,k]+b.y[i,j+1,k-1]+b.y[i,j,k-1])

h.z[i,j,k]  = C.zx[i,j,k] * bx  + C.zy[i,j,k] * by  + C.zz[i,j,k] * b.z[i,j,k]

		end
	end
end

for k = 1:1
	for j = 1:1
		for i = 2:nx-1
by = 0.25 * (b.y[i,j,k] + b.y[i-1,j,k] + b.y[i-1,j+1,k] + b.y[i,j+1,k])
bz = 0.25 * (b.z[i,j,k] + b.z[i,j,k+1] + b.z[i-1,j,k+1] + b.z[i-1,j,k])
h.x[i,j,k]  = C.xx[i,j,k] * b.x[i,j,k]  + C.xy[i,j,k] * by  + C.xz[i,j,k] * bz
		end
	end
end
for k = 1:1
	for j = 2:ny-1
		for i = 1:1
bx = 0.25*(b.x[i,j,k]+b.x[i,j-1,k]+b.x[i+1,j-1,k]+b.x[i+1,j,k])
bz = 0.25*(b.z[i,j,k] + b.z[i,j-1,k]+ b.z[i,j-1,k+1] + b.z[i,j,k+1])
h.y[i,j,k]  = C.yx[i,j,k] * bx  + C.yy[i,j,k] * b.y[i,j,k]  + C.yz[i,j,k] * bz
		end
	end
end
for k = 2:ny-1
	for j = 1:1
		for i = 1:1
bx = 0.25*(b.x[i,j,k]+b.x[i+1,j,k]+b.x[i+1,j,k-1]+b.x[i,j,k-1])
by = 0.25*(b.y[i,j,k]+b.y[i,j+1,k]+b.y[i,j+1,k-1]+b.y[i,j,k-1])

h.z[i,j,k]  = C.zx[i,j,k] * bx  + C.zy[i,j,k] * by  + C.zz[i,j,k] * b.z[i,j,k]
		end
	end
end

end

# ╔═╡ d87ef4f9-e4e9-4b1d-907d-054505d3ff54
begin
	ref_grid = StructuredGrid(LinRange(-5,15,30),LinRange(-5,17,10),LinRange(-5,1,101))

	e = VectorField{Float64}(YeeGrid{Primal}(ref_grid))
	d = VectorField{Float64}(YeeGrid{Primal}(ref_grid))
	ϵ = TensorField{Float64}(YeeGrid{Primal}(ref_grid))

	b = VectorField{Float64}(YeeGrid{Dual}(ref_grid))
	h = VectorField{Float64}(YeeGrid{Dual}(ref_grid))
	μ = TensorField{Float64}(YeeGrid{Dual}(ref_grid))

	j =  VectorField{Float64}(YeeGrid{Primal}(ref_grid))
	σₑ = TensorField{Float64}(YeeGrid{Primal}(ref_grid))

	m =  VectorField{Float64}(YeeGrid{Dual}(ref_grid))
	σₘ = TensorField{Float64}(YeeGrid{Dual}(ref_grid))

	cf = VectorField{Float64}(YeeGrid{Dual}(ref_grid))
	ca = VectorField{Float64}(YeeGrid{Primal}(ref_grid))

end

# ╔═╡ 03ac2bca-403e-478e-b947-208427fdf579
begin
	FaradayStep!(b,e,m,cf, ref_grid = ref_grid)
	_update_dual_constitutive!(h,b,μ,size(ref_grid)...)
	_update_dual_constitutive!(m,h,σₘ,size(ref_grid)...)

	AmpereStep!(d,h,j,ca, ref_grid = ref_grid)
	_update_primal_constitutive!(e,d,ϵ,size(ref_grid)...)
	_update_primal_constitutive!(j,e,σₑ,size(ref_grid)...)

end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"

[compat]
Parameters = "~0.12.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"
"""

# ╔═╡ Cell order:
# ╠═456df07b-1be8-4adb-9302-9290314311e3
# ╠═4f7a7193-e75e-4756-9423-7d81e2f7e4db
# ╠═1bc97871-d861-47d4-bde9-e2373b1e903a
# ╠═4b3b120f-52b0-47ae-bec2-207ea12ddd6a
# ╠═6ad42c82-25c5-41c3-83f1-0ff002edbafe
# ╠═f32b184a-ba70-4e80-8f8d-e12def369128
# ╠═4fd9f13e-8279-41e2-b7f5-e9865ed869be
# ╠═fcb5f651-3e20-4181-9e72-9310f823adf0
# ╠═7f1e505f-7223-4ae8-a6d4-9ff455553256
# ╠═8176e6dd-0768-4d9b-a2b7-957236e4c55b
# ╠═d87ef4f9-e4e9-4b1d-907d-054505d3ff54
# ╠═03ac2bca-403e-478e-b947-208427fdf579
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
Pkg.add(["TimeSeries", "MarketData", "Plots", "Pipe"])
using Pipe: @pipe
using TimeSeries, MarketData, Plots
# include("src/Portfolio.jl")
using Portfolio

function strategy_tearsheet(weights::Weights, returns::Returns)
    strategy_returns = run_weights(returns, weights, 0.0)
    equity_curves = plot(TimeSeries.merge(equity_curve(returns), equity_curve(strategy_returns)))
    plot(
        equity_curves,
        plot(weights),
        size = (1200, 1200),
        layout = grid(2, 1, heights = [0.8, 0.2])
    )
end

@time dell = DELL[:AdjClose] |> percentchange

# Target mean volatility, 30 day lookback
weight_vol_target_mean = @pipe strategy_vol_target_mean(dell, 30) |> lowpass(1.0).(_)
strategy_tearsheet(weight_vol_target_mean, dell)

# Ensemble lookback 20:40, still target the mean volatility
ensembled_weight_vol_target = @pipe ensemble_parameters(strategy_vol_target_mean, dell, 20:40) |> lowpass(1.0).(_)
strategy_tearsheet(ensembled_weight_vol_target, dell)

# If 200 day momentum is positive, 100% long, otherwise zero allocation
weight_momentum_above_0 = strategy_positive_mom_binary(dell, 200)
strategy_tearsheet(weight_momentum_above_0, dell)

# Ensemble 50:200 day momentum
ensembled_weight_momentum_above_0 = ensemble_parameters(strategy_positive_mom_binary, dell, 50:200)
strategy_tearsheet(ensembled_weight_momentum_above_0, dell)

# Merge the ensembled mom / vol targeting weights
all_weights_ensembled = merge_weights([ensembled_weight_vol_target, ensembled_weight_momentum_above_0])
strategy_tearsheet(all_weights_ensembled, dell)


text:
function accuracy(y_pred, y)
    return mean(onecold(y_pred) .== onecold(y))
end


text:
# Run this file using julia
# Ex: julia runme.jl
import GitWorkers
GitWorkers._runme(@__DIR__, ARGS)

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
mvs = []

function reverse(Side::HorizonSide)
    return HorizonSide((mod(Int(Side)+2, 4)))
end

function main(r::Robot)
    for i in [0, 1, 2, 3]
        while !isborder(r, HorizonSide(i)) && !isborder(r, HorizonSide(mod(i+1,4)))
            move!(r, HorizonSide(i))
            push!(mvs, HorizonSide(i))
            move!(r, HorizonSide(mod(i+1,4)))
            putmarker!(r)
            push!(mvs, HorizonSide(mod(i+1,4)))
        end
        while length(mvs) != 0
            move!(r, reverse(pop!(mvs)))
        end
    end
    putmarker!(r)
end


text:
function relabelSitesWithIntIndex(
            unitcell :: U
        ) :: Unitcell{Site{Int64,D},B} where {LS,LB,D,N,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},U<:AbstractUnitcell{S,B}}

    # create a list of new sites
    sitelist = Site{Int64,D}[
        newSite(Site{Int64,D}, point(site(unitcell,i)), i)
        for i in 1:numSites(unitcell)
    ]

    # return a new lattice object
    return newUnitcell(
        Unitcell{Site{Int64,D},B},
        deepcopy(latticeVectors(unitcell)),
        sitelist,
        deepcopy(bonds(unitcell))
    )
end

# export the relabeling function
export relabelSitesWithIntIndex


text:
function get_default_env(; version::VersionNumber=GCC_MIN_VER, is_cxx=false)
    p = HostPlatform()
    if os(p) == "macos"
        return MacEnv(p, version, is_cxx)
    elseif os(p) == "windows"
        return WindowsEnv(p, version, is_cxx)
    elseif arch(p) == "armv7l"
        return ArmEnv(p, version, is_cxx)
    elseif libc(p) == "musl"
        return MuslEnv(p, version, is_cxx)
    elseif libc(p) == "gnu"
        return GnuEnv(p, version, is_cxx)
    else
        @warn "unknown platform! using default GNU environment."
        return GnuEnv(p, version, is_cxx)
    end
end

function get_system_includes(env::AbstractJLLEnv=get_default_env())
    info = get_environment_info(env.platform, env.gcc_version)
    gcc_info = info[1]
    sys_info = info[2]

    # download shards
    download_artifact(Base.SHA1(gcc_info.id), gcc_info.url, gcc_info.chk)
    # download_artifact(Base.SHA1(sys_info.id), sys_info.url, sys_info.chk)

    # -isystem paths
    gcc_prefix = artifact_path(Base.SHA1(gcc_info.id))
    # sys_prefix = artifact_path(Base.SHA1(sys_info.id))

    isys = String[]
    get_system_includes!(env, gcc_prefix, isys)

    for dir in isys
        @assert isdir(dir) "failed to setup environment due to missing dir: $dir, please file an issue."
    end

    return normpath.(isys)
end

function get_system_includes!(env::MacEnv, prefix::String, isys::Vector{String})
    p = env.platform
    triple = __triplet(env.platform)
    version = env.gcc_version
    if os(p) == "macos" && arch(p) == "x86_64"
        if env.is_cxx
            push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include", "c++", "v1"))
            push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include"))
            push!(isys, joinpath(prefix, triple, "include", "c++", string(version)))
            push!(isys, joinpath(prefix, triple, "include", "c++", string(version), triple))
            push!(isys, joinpath(prefix, triple, "include", "c++", string(version), "backward"))
            push!(isys, joinpath(prefix, triple, "include"))
            push!(isys, joinpath(prefix, triple, "sys-root", "System", "Library", "Frameworks"))
        else
            push!(isys, joinpath(prefix, "lib", "gcc", triple, string(version), "include"))
            push!(isys, joinpath(prefix, "lib", "gcc", triple, string(version), "include-fixed"))
            push!(isys, joinpath(prefix, triple, "include"))
            push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include"))
            push!(isys, joinpath(prefix, triple, "sys-root", "System", "Library", "Frameworks"))
        end
    else  # "aarch64-apple-darwin20"
        if env.is_cxx
            push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include", "c++", "v1"))
            push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include"))
            push!(isys, joinpath(prefix, triple, "include", "c++", "11.0.0"))
            push!(isys, joinpath(prefix, triple, "include", "c++", "11.0.0", triple))
            push!(isys, joinpath(prefix, triple, "include", "c++", "11.0.0", "backward"))
        else
            push!(isys, joinpath(prefix, "lib", "gcc", triple, "11.0.0", "include"))
            push!(isys, joinpath(prefix, "lib", "gcc", triple, "11.0.0", "include-fixed"))
            push!(isys, joinpath(prefix, triple, "include"))
            push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include"))
            push!(isys, joinpath(prefix, triple, "sys-root", "System", "Library", "Frameworks"))
        end
    end
end

function get_system_includes!(env::WindowsEnv, prefix::String, isys::Vector{String})
    triple = __triplet(env.platform)
    version = env.gcc_version
    if env.is_cxx
        push!(isys, joinpath(prefix, triple, "include", "c++", string(version)))
        push!(isys, joinpath(prefix, triple, "include", "c++", string(version), triple))
        push!(isys, joinpath(prefix, triple, "include", "c++", string(version), "backward"))
        push!(isys, joinpath(prefix, triple, "include"))
        push!(isys, joinpath(prefix, triple, "sys-root", "include"))
    else
        push!(isys, joinpath(prefix, "lib", "gcc", triple, string(version), "include"))
        push!(isys, joinpath(prefix, "lib", "gcc", triple, string(version), "include-fixed"))
        push!(isys, joinpath(prefix, triple, "include"))
        push!(isys, joinpath(prefix, triple, "sys-root", "include"))
    end
end

function get_system_includes!(env::GnuEnv, prefix::String, isys::Vector{String})
    triple = __triplet(env.platform)
    version = env.gcc_version
    if env.is_cxx
        push!(isys, joinpath(prefix, triple, "include", "c++", string(version)))
        push!(isys, joinpath(prefix, triple, "include", "c++", string(version), triple))
        push!(isys, joinpath(prefix, triple, "include", "c++", string(version), "backward"))
        push!(isys, joinpath(prefix, triple, "include"))
        push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include"))
    else
        push!(isys, joinpath(prefix, "lib", "gcc", triple, string(version), "include"))
        push!(isys, joinpath(prefix, "lib", "gcc", triple, string(version), "include-fixed"))
        push!(isys, joinpath(prefix, triple, "include"))
        push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include"))
    end
end

function get_system_includes!(env::MuslEnv, prefix::String, isys::Vector{String})
    triple = __triplet(env.platform)
    version = env.gcc_version
    if env.is_cxx
        push!(isys, joinpath(prefix, triple, "include", "c++", string(version)))
        push!(isys, joinpath(prefix, triple, "include", "c++", string(version), triple))
        push!(isys, joinpath(prefix, triple, "include", "c++", string(version), "backward"))
        push!(isys, joinpath(prefix, triple, "include"))
        push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include"))
    else
        push!(isys, joinpath(prefix, "lib", "gcc", triple, string(version), "include"))
        push!(isys, joinpath(prefix, triple, "include"))
        push!(isys, joinpath(prefix, triple, "sys-root", "usr", "include"))
    end
end

function get_system_includes!(env::ArmEnv, prefix::String, isys::Vector{String})
    triple = __triplet(env.platform)
    version = env.gcc_version
    if env.platform == Platform("armv7l", "linux")
        if env.is_cxx
            push!(isys, joinpath(prefix, "arm-linux-gnueabihf", "include", "c++", string(version)))
            push!(isys, joinpath(prefix, "arm-linux-gnueabihf", "include", "c++", string(version), "arm-linux-gnueabihf"))
            push!(isys, joinpath(prefix, "arm-linux-gnueabihf", "include", "c++", string(version), "backward"))
            push!(isys, joinpath(prefix, "arm-linux-gnueabihf", "include"))
            push!(isys, joinpath(prefix, "arm-linux-gnueabihf", "sys-root", "usr", "include"))
        else
            push!(isys, joinpath(prefix, "lib", "gcc", "arm-linux-gnueabihf", string(version), "include"))
            push!(isys, joinpath(prefix, "lib", "gcc", "arm-linux-gnueabihf", string(version), "include-fixed"))
            push!(isys, joinpath(prefix, "arm-linux-gnueabihf", "include"))
            push!(isys, joinpath(prefix, "arm-linux-gnueabihf", "sys-root", "usr", "include"))
        end
    else  # "armv7l-linux-musleabihf"
        if env.is_cxx
            push!(isys, joinpath(prefix, "arm-linux-musleabihf", "include", "c++", string(version)))
            push!(isys, joinpath(prefix, "arm-linux-musleabihf", "include", "c++", string(version), "arm-linux-musleabihf"))
            push!(isys, joinpath(prefix, "arm-linux-musleabihf", "include", "c++", string(version), "backward"))
            push!(isys, joinpath(prefix, "arm-linux-musleabihf", "include"))
            push!(isys, joinpath(prefix, "arm-linux-musleabihf", "sys-root", "usr", "include"))
        else
            push!(isys, joinpath(prefix, "lib", "gcc", "arm-linux-musleabihf", string(version), "include"))
            push!(isys, joinpath(prefix, "lib", "gcc", "arm-linux-musleabihf", string(version), "include-fixed"))
            push!(isys, joinpath(prefix, "arm-linux-musleabihf", "include"))
            push!(isys, joinpath(prefix, "arm-linux-musleabihf", "sys-root", "usr", "include"))
        end
    end
end


text:
using UUIDs

mutable struct Family
    id::UUID
    size::Int
    applications::Vector{Application}
    Family() = new(UUIDs.uuid4(), 0, Application[])
    Family(id, s, apps) = new(id, s, apps)
end

id(f::Family) = f.id
applications(f::Family) = f.applications
Base.size(f::Family) = f.size

jurisdiction(f::Family) = unique(jurisdiction(a) for a in applications(f))
dates(f::Family) = [date(a) for a in applications(f)]
earliest_filing(f::Family) = minimum(dates(f))
latest_filing(f::Family) = maximum(dates(f))

title(f::Family) = [title(a) for a in applications(f)]
abstract(f::Family) = [abstract(a) for a in applications(f)]

function applicants(f::Family)
    apps = applications(f)
    appl = reduce(vcat, applicants.(apps))
    return unique(appl)
end

function inventors(f::Family)
    apps = applications(f)
    inv = reduce(vcat, inventors.(apps))
    return unique(inv)
end

function classification(fam::Family)
	apps = applications(fam)
	class = reduce(vcat, classification.(apps))
	return unique(class)
end

cites(f::Family) = reduce(vcat, cites(a) for a in applications(f)) |> unique
cites_count(f::Family) = length(cites(f))
cites_npl(f::Family) = reduce(vcat, cites_npl(a) for a in applications(f)) |> unique
citedby(f::Family) = reduce(vcat, citedby(a) for a in applications(f)) |> unique
citedby_count(f::Family) = length(citedby(f))

# function citedby_count(f::Family)
#     res = Set{String}()
#     count = 0
#     for a in applications(f)
#         for c in citedby(a)
#             incr = ifelse(id(c) in res, 0, 1)
#             push!(res, id(c))
#             count += incr
#         end
#     end
#     count
# end

function Base.show(io::IO, f::Family)
    print(io, "Family with $(f.size) members" )
end

function Base.show(io::IO, ::MIME"text/plain", f::Family)
    tit = reduce(vcat, title.(applications(f)))
    if length(tit) == 0
        tit = ""
    else
        idx_tit = findfirst(lang.(tit) .== "en")
        tit = isnothing(idx_tit) ? text(first(tit)) : text(tit[idx_tit])
    end

    abs = reduce(vcat, abstract.(applications(f)))
    if length(abs) == 0
        abs = ""
    else
        idx_abs = findfirst(lang.(abs) .== "en")
        abs = isnothing(idx_abs) ? text(first(abs)) : text(abs[idx_abs])
    end

    println(io, "----------------------")
    println(io, "Size: $(f.size)")
    println(io, "Earliest filing: $(earliest_filing(f))")
    println(io, "Latest filing: $(maximum(dates(f)))")
    println(io, "Cited by: $(citedby(f) |> length)")
    println(io, "Filed in: $(join(jurisdiction(f), ", "))")
    println(io, "Applicants: $(join(applicants(f), ", "))")
    println(io, "Subgroups: $(join(subgroup.(classification(f)), ", "))")
    println(io, "----------------------")
    println(io, "Title: $tit\n")
    println(io, "Abstract: $abs")
end

# function Base.show(io::IO, ::MIME"text/plain", f::Family)

# end

# function Base.show(io::IO, ::MIME"text/markdown", f::Family)
#     tit = reduce(vcat, title.(applications(f)))
#     idx_tit = findfirst(lang.(tit) .== "en")
#     tit = isnothing(idx_tit) ? text(first(tit)) : text(tit[idx_tit])
#     abs = reduce(vcat, abstract.(applications(f)))
#     idx_abs = findfirst(lang.(abs) .== "en")
#     abs = isnothing(idx_abs) ? text(first(abs)) : text(abs[idx_abs])

#     md"""
#     ##### $tit

#     **Size:** $(f.size)
#     **Earliest filing: $(earliest)
#     **Earliest filing:** $(earliest_filing(f))
#     **Latest filing:** $(maximum(dates(f)))
#     **Cited by:** $(citedby(f) |> length)
#     **Filed in:** $(join(jurisdiction(f), ", "))
#     **Applicants:** $(join(applicants(f), ", "))
#     ---
#     **Abstract:** $abs
#     """
# end

function aggregate_families(apps::Vector{Application})
    visited = Dict(id(a) => false for a in apps)
    idx = Dict(id(a) => i for (i, a) in enumerate(apps))
    families = Family[]
    for a in apps
        visited[id(a)] && continue
        applications = Application[]
        push!(applications, a)
        for s in siblings(a)
            haskey(idx, id(s)) || continue
            push!(applications, apps[idx[id(s)]])
            visited[id(s)] = true
        end
        push!(families, Family(UUIDs.uuid4(), familysize(a), applications))
    end
    return families
end


text:
using Documenter, AuraLighting

makedocs(sitename = "AuraLighting Module Documentation", format = Documenter.HTML(prettyurls = false))


text:
__precompile__(false)
module IBKRCore

using CxxWrap
using CxxWrap.CxxWrapCore
using CxxWrap.CxxWrapCore: CxxInt, CxxLong
using MacroTools

using ...Emporos: PACKAGE_ROOT, @debugcb

const TWS_DEFAULT_PORT = 7497

include("types.jl")

# NOTE: Functions called after @wrapmodule extend those defined in
# @wrapmodule while those defined before can be referenced in C
@wrapmodule(abspath(joinpath(PACKAGE_ROOT, "../result/lib", "libIBApiWrapper")))
Base.convert(::Type{Any}, x::CxxWrap.CxxWrapCore.SmartPointer{DerivedT}) where {BaseT, DerivedT<:BaseT} = x

function __init__()
    @initcxx
end

# TODO See: https://github.com/JuliaInterop/CxxWrap.jl/issues/301
Base.convert(::Type{Any}, x::CxxWrap.CxxWrapCore.SmartPointer{DerivedT}) where {BaseT, DerivedT<:BaseT} = x

include("jlclient.jl")
include("ibclient.jl")

# Generated
include("generated/ComboLeg.jl")
include("generated/Contract.jl")

end


text:
function value(obj::MultivariateVectorFunction, arg0::Vector{jdouble})
    return jcall(obj, "value", Vector{jdouble}, (Vector{jdouble},), arg0)
end


text:
# The functions add the calculated values to the existing values of the output array.
# Therefore, if the derivative values themselves are desired, pass the output array
# initialized with zeros.

# To-dos
# - Test if using a separate index vector (which is an identity map) to iterate over a 3D
# array makes iteration significantly slower.  This is to simulate the case of FEM.

export apply_curl!

# Wrapper to apply the discrete curl by default
apply_curl!(G::AbsArrNumber{K₊₁},  # output field; G[i,j,k,w] is w-component of G at (i,j,k)
            F::AbsArrNumber{K₊₁},  # input field; F[i,j,k,w] is w-component of F at (i,j,k)
            ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
            isfwd::AbsVecBool,  # isfwd[w] = true|false: create ∂w by forward|backward difference
            ∆l⁻¹::NTuple{K,Number}=ntuple(k->1.0,length(isfwd)),  # ∆l⁻¹[w]: inverse of uniform distance between grid planes in x-direction
            isbloch::AbsVecBool=fill(true,K),  # boundary conditions in x, y, z
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(K);  # Bloch phase factor in x, y, z
            cmp_shp::AbsVecInteger=1:K,
            cmp_out::AbsVecInteger=1:size(G,K₊₁),
            cmp_in::AbsVecInteger=1:size(F,K₊₁),
            α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α ∇×F
            ) where {K,K₊₁,OP} =
    (N = size(G)[1:K]; apply_curl!(G, F, Val(OP), isfwd, fill.(∆l⁻¹,N), isbloch, e⁻ⁱᵏᴸ; cmp_shp, cmp_out, cmp_in, α))

# Wrapper for converting AbstractVector's to SVec's
apply_curl!(G::AbsArrNumber{K₊₁},  # output field; G[i,j,k,w] is w-component of G at (i,j,k)
            F::AbsArrNumber{K₊₁},  # input field; F[i,j,k,w] is w-component of F at (i,j,k)
            ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
            isfwd::AbsVecBool,  # isfwd[w] = true|false: create ∂w by forward|backward difference
            ∆l⁻¹::NTuple{K,AbsVecNumber},  # ∆l⁻¹[w]: inverse of distances between grid planes in x-direction
            isbloch::AbsVecBool=fill(true,K),  # boundary conditions in x, y, z
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(K);  # Bloch phase factor in x, y, z
            cmp_shp::AbsVecInteger=1:K,
            cmp_out::AbsVecInteger=1:size(G,K₊₁),
            cmp_in::AbsVecInteger=1:size(F,K₊₁),
            α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α ∇×F
            ) where {K,K₊₁,OP} =
    # I should not cast e⁻ⁱᵏᴸ into a complex vector, because then the entire curl matrix
    # becomes a complex matrix.  Sometimes I want to keep it real (e.g., when no PML and
    # Bloch phase factors are used).  In fact, this is the reason why I accept e⁻ⁱᵏᴸ instead
    # of constructing it from k and L as exp.(-im .* k .* L), which is always complex even
    # if k = 0.
    #
    # I should not cast ∆l⁻¹ to a vector of any specific type (e.g., Float, ComplexF), either,
    # because sometimes I would want to even create an integral curl operator.
    (Kout = length(cmp_out); Kin = length(cmp_in);
     apply_curl!(G, F, Val(OP), SVec{K}(isfwd), ∆l⁻¹, SVec{K}(isbloch), SVec{K}(e⁻ⁱᵏᴸ);
                 cmp_shp=SInt{K}(cmp_shp), cmp_out=SInt{Kout}(cmp_out), cmp_in=SInt{Kin}(cmp_in), α))

# Concrete implementation
function apply_curl!(G::AbsArrNumber{K₊₁},  # output field; G[i,j,k,w] is w-component of G at (i,j,k)
                     F::AbsArrNumber{K₊₁},  # input field; F[i,j,k,w] is w-component of F at (i,j,k)
                     ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                     isfwd::SBool{K},  # isfwd[w] = true|false: create ∂w by forward|backward difference
                     ∆l⁻¹::NTuple{K,AbsVecNumber},  # ∆l⁻¹[w]: inverse of distances between grid planes in x-direction
                     isbloch::SBool{K},  # boundary conditions in x, y, z
                     e⁻ⁱᵏᴸ::SNumber{K};  # Bloch phase factor in x, y, z
                     cmp_shp::SInt{K}=SVec(ntuple(identity, Val(K))),
                     cmp_out::SInt{Kout}=SVec(ntuple(identity, size(G,K₊₁))),
                     cmp_in::SInt{Kin}=SVec(ntuple(identity, size(G,K₊₁))),
                     α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α ∇×F
                     ) where {K,K₊₁,Kout,Kin,OP}
    @assert K₊₁ == K+1

    for ind_nv = 1:Kout
        nv = cmp_out[ind_nv]  # Cartesian compotent of output field
        Gv = selectdim(G, K₊₁, ind_nv)  # v-component of output field

        op = Val(OP)
        for ind_nu = 1:Kin
            nu = cmp_in[ind_nu]  # Cartesian compotent of input field
            parity = CURL_BLK[nv,nu]  # CURL_BLK defined in matrix/curl.jl

            if iszero(parity)
                if op == Val(:(=))
                    Gv .= 0
                    op = Val(:(+=))
                end
                continue
            end

            nw = 6 - nv - nu  # direction of differentiation; 6 = 1 + 2 + 3
            is_nw = cmp_shp.==nw
            sum(is_nw)==1 || @error "cmp_shp = $cmp_shp does not have one and only one nw = $nw"
            ind_nw = findfirst(is_nw)

            Fu = selectdim(F, K₊₁, ind_nu)  # u-component of input field
            apply_∂!(Gv, Fu, op, ind_nw, isfwd[ind_nw], ∆l⁻¹[ind_nw], isbloch[ind_nw], e⁻ⁱᵏᴸ[ind_nw], α=parity*α)  # Gv += α (±∂Fu/∂w)

            op = Val(:(+=))

            # # Equivalent to the above code:
            # if op == Val(:(=))
            #     if iszero(parity)
            #         Gv .= 0
            #     else
            #         nw = 6 - nv - nu  # direction of differentiation; 6 = 1 + 2 + 3
            #         is_nw = cmp_shp.==nw
            #         @assert sum(is_nw)==1  # in cmp_shp, one and only one entry is nw
            #         ind_nw = findfirst(is_nw)
            #
            #         Fu = @view F[:,:,:,ind_nu]  # u-component of input field
            #         apply_∂!(Gv, Fu, op, ind_nw, isfwd[ind_nw], ∆l⁻¹[ind_nw], isbloch[ind_nw], e⁻ⁱᵏᴸ[ind_nw], α=parity*α)  # Gv += α (±∂Fu/∂w)
            #     end
            #
            #     op = Val(:(+=))
            # else  # op == Val(:(+=))
            #     if iszero(parity)
            #         # Do nothing.
            #     else
            #         nw = 6 - nv - nu  # direction of differentiation; 6 = 1 + 2 + 3
            #         is_nw = cmp_shp.==nw
            #         @assert sum(is_nw)==1  # in cmp_shp, one and only one entry is nw
            #         ind_nw = findfirst(is_nw)
            #
            #         Fu = @view F[:,:,:,ind_nu]  # u-component of input field
            #         apply_∂!(Gv, Fu, op, ind_nw, isfwd[ind_nw], ∆l⁻¹[ind_nw], isbloch[ind_nw], e⁻ⁱᵏᴸ[ind_nw], α=parity*α)  # Gv += α (±∂Fu/∂w)
            #     end
            # end
        end
    end

    # @sync begin
    #     @spawn apply_∂!(@view(G[:,:,:,1]), @view(F[:,:,:,3]), Val(OP), 2, isfwd[2], ∆l⁻¹[2], isbloch[2], e⁻ⁱᵏᴸ[2], α=α)  # Gv += α (±∂Fu/∂w)
    #     @spawn apply_∂!(@view(G[:,:,:,2]), @view(F[:,:,:,1]), Val(OP), 3, isfwd[3], ∆l⁻¹[3], isbloch[3], e⁻ⁱᵏᴸ[3], α=α)  # Gv += α (±∂Fu/∂w)
    #     @spawn apply_∂!(@view(G[:,:,:,3]), @view(F[:,:,:,2]), Val(OP), 1, isfwd[1], ∆l⁻¹[1], isbloch[1], e⁻ⁱᵏᴸ[1], α=α)  # Gv += α (±∂Fu/∂w)
    # end
    #
    # @sync begin
    #     @spawn apply_∂!(@view(G[:,:,:,1]), @view(F[:,:,:,2]), Val(:(+=)), 3, isfwd[3], ∆l⁻¹[3], isbloch[3], e⁻ⁱᵏᴸ[3], α=-α)  # Gv += α (±∂Fu/∂w)
    #     @spawn apply_∂!(@view(G[:,:,:,2]), @view(F[:,:,:,3]), Val(:(+=)), 1, isfwd[1], ∆l⁻¹[1], isbloch[1], e⁻ⁱᵏᴸ[1], α=-α)  # Gv += α (±∂Fu/∂w)
    #     @spawn apply_∂!(@view(G[:,:,:,3]), @view(F[:,:,:,1]), Val(:(+=)), 2, isfwd[2], ∆l⁻¹[2], isbloch[2], e⁻ⁱᵏᴸ[2], α=-α)  # Gv += α (±∂Fu/∂w)
    # end

    return nothing
end


text:
@testset "DB operations" begin
    db = setup_db("test/data/TEST_PowerFlowNetworks_SQLite.sqlite"; delete_if_exists=true)
    load_instance_in_db!(db, "case6ww", 0, "test/data/case6ww.m", "MATPOWERM", Dates.now())
    load_instance_in_db!(db, "C2FEN02312", 3, "test/data/C2FEN02312/scenario_3/case.raw", "RAWGO", Dates.now())
    load_instance_in_db!(db, "C2S6N02045", 1, "test/data/C2S6N02045/scenario_1/", "RAWGO", Dates.now())
    @test table_count(db, "instances") == 3

    save_basic_features_instances!(db)
    save_features_instances!(db)
    serialize_instances!(db;
                         serialize_path="test/data/networks_serialized",
                         graphs_path="test/data/graphs")
    generate_decompositions!(db;
                             cliques_path="test/data/cliques",
                             cliquetrees_path="test/data/cliquetrees",
                             graphs_path="test/data/graphs",
                             extension_alg="cholesky",
                             preprocess_path="test/data/configs/preprocess_default.json",
                             subset=[1, 2])
    @test table_count(db, "decompositions") == 2
    generate_decompositions!(db;
                             cliques_path="test/data/cliques",
                             cliquetrees_path="test/data/cliquetrees",
                             graphs_path="test/data/graphs",
                             extension_alg="minimum_degree",
                             preprocess_path="test/data/configs/preprocess_add_random_50.json",
                             subset=[3])
    @test table_count(db, "decompositions") == 3
    merge_decompositions!(db;
                          heuristic=["molzahn"],
                          heuristic_switch=[0],
                          treshold_name="clique_nv_up",
                          merge_kwargs=Dict("treshold_percent" => 0.5),
                          subset=[1, 2])
    @test table_count(db, "decompositions") == 5
    @test table_count(db, "merges") == 2
    generate_decompositions!(db;
                             cliques_path="test/data/cliques",
                             cliquetrees_path="test/data/cliquetrees",
                             graphs_path="test/data/graphs",
                             extension_alg="minimum_degree",
                             preprocess_path="test/data/configs/preprocess_add_random_50.json",
                             subset=[1])
    combine_decompositions!(db; how="vertices_intersect", extension_alg="cholesky", subset=[1])
    combine_decompositions!(db; how="vertices_union", extension_alg="minimum_degree", subset=[1])
    combine_decompositions!(db; how="clique_intersect", extension_alg="cholesky", subset=[1])
    combine_decompositions!(db; how="clique_union", extension_alg="minimum_degree", subset=[1])
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
type Instance
    nb::Int
    rect::Vector{Vector{Float64}}
    Instance() = new(0, [])
end

inst = Instance()

open("../instance.txt") do f

    lines = readlines(f)

    inst.nb = parse(Int, lines[1])

    for i in 2:(inst.nb+1)
        l = split(lines[i], " ")
        push!(inst.rect, [parse(Float64, l[1]), parse(Float64, l[2])])
    end

end