text:
# Functions for querying network properties.

# DEPRECIATED after v9.0
function params(network)
    Base.depwarn("`params` is depreciated, please use `ModelingToolkit.parameters` for all system and subsystem parameters, or `reactionparams` for all parameters within system and subsystem `Reaction`s.", :params, force=true)
    parameters(network)
end

function numparams(network)
    Base.depwarn("`numparams` is depreciated, please use `length(ModelingToolkit.parameters)` for the total number of parameters across all systems and subsystems, or `numreactionparams` for the number of parameters within system and subsystem `Reaction`s.", :params, force=true)
    length(parameters(network))
end

"""
    merge(network1::ReactionSystem, network2::ReactionSystem)

Create a new network merging `network1` and `network2`.

Notes:
- Duplicate reactions between the two networks are not filtered out.
- [`Reaction`](@ref)s are not deepcopied to minimize allocations, so the new
  network will share underlying data arrays.
- Subsystems are not deepcopied between the two networks and will hence be
  shared.
- Returns the merged network.
"""
function Base.merge(network1::ReactionSystem, network2::ReactionSystem)
    Base.depwarn("`merge(sys1::ReactionSystem, sys2::ReactionNetwork)` is depreciated, please use `ModelingToolkit.extend` instead.", :merge, force=true)
    network = make_empty_network()
    merge!(network, network1)
    merge!(network, network2)
    network
end

######### Accessors: #########

function filter_nonrxsys(network)
    systems   = get_systems(network)
    rxsystems = ReactionSystem[]
    for sys in systems
        (sys isa ReactionSystem) && push!(rxsystems, sys)
    end
    rxsystems
end

function species(network::ReactionSystem, sts)
    [MT.renamespace(network,st) for st in sts]
end

"""
    species(network)

Given a [`ReactionSystem`](@ref), return a vector of all species defined in the system
and any subsystems that are of type `ReactionSystem`. To get the variables in the system
and all subsystems, including non-`ReactionSystem` subsystems, uses `states(network)`.

Notes:
- If `ModelingToolkit.get_systems(network)` is non-empty will allocate.
"""
function species(network)
    sts = get_states(network)
    systems = filter_nonrxsys(network)
    isempty(systems) && return sts
    unique([sts; reduce(vcat, map(sys -> species(sys,species(sys)), systems))])
end

"""
    reactionparams(network)

Given a [`ReactionSystem`](@ref), return a vector of all parameters defined
within the system and any subsystems that are of type `ReactionSystem`. To get
the parameters in the system and all subsystems, including non-`ReactionSystem`
subsystems, use `parameters(network)`.

Notes:
- If `ModelingToolkit.get_systems(network)` is non-empty will allocate.
"""
function reactionparams(network)
    ps = get_ps(network)
    systems = filter_nonrxsys(network)
    isempty(systems) && return ps
    unique([ps; reduce(vcat, map(sys -> species(sys,reactionparams(sys)), systems))])
end

function namespace_reactions(network::ReactionSystem)
    rxs = reactions(network)
    isempty(rxs) && return Reaction[]
    map(rx -> namespace_equation(rx, network), rxs)
end

"""
    reactions(network)

Given a [`ReactionSystem`](@ref), return a vector of all `Reactions` in the system.

Notes:
- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
"""
function reactions(network)
    rxs = get_eqs(network)
    systems = filter_nonrxsys(network)
    isempty(systems) && (return rxs)
    [rxs; reduce(vcat, namespace_reactions.(systems); init=Reaction[])]
end

"""
    speciesmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping species that
participate in `Reaction`s to their index within [`species(network)`](@ref).
"""
function speciesmap(network)
    nps = get_networkproperties(network)
    if isempty(nps.speciesmap)
        nps.speciesmap = Dict(S => i for (i,S) in enumerate(species(network)))
    end
    nps.speciesmap
end

"""
    paramsmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping from all
parameters that appear within the system to their index within
`parameters(network)`.
"""
function paramsmap(network)
    Dict(p => i for (i,p) in enumerate(parameters(network)))
end

"""
    reactionparamsmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping from parameters that
appear within `Reaction`s to their index within [`reactionparams(network)`](@ref).
"""
function reactionparamsmap(network)
    Dict(p => i for (i,p) in enumerate(reactionparams(network)))
end

"""
    numspecies(network)

Return the total number of species within the given [`ReactionSystem`](@ref) and
subsystems that are `ReactionSystem`s.

Notes
- If there are no subsystems this will be fast.
- As this calls [`species`](@ref), it can be slow and will allocate if there are
  any subsystems.
"""
function numspecies(network)
    length(species(network))
end

"""
    numreactions(network)

Return the total number of reactions within the given [`ReactionSystem`](@ref)
and subsystems that are `ReactionSystem`s.
"""
function numreactions(network)
    nr = length(get_eqs(network))
    for sys in get_systems(network)
        (sys isa ReactionSystem) && (nr += numreactions(sys))
    end
    nr
end

"""
    numreactionparams(network)

Return the total number of parameters within the given [`ReactionSystem`](@ref)
and subsystems that are `ReactionSystem`s.

Notes
- If there are no subsystems this will be fast.
- As this calls [`reactionparams`](@ref), it can be slow and will allocate if
  there are any subsystems.
"""
function numreactionparams(network)
    length(reactionparams(network))
end

"""
    dependents(rx, network)

Given a [`Reaction`](@ref) and a [`ReactionSystem`](@ref), return a vector of the
*non-constant* species the reaction rate law depends on. e.g., for

`k*W, 2X + 3Y --> 5Z + W`

the returned vector would be `[W(t),X(t),Y(t)]`.

Notes:
- Allocates
- Does not check for dependents within any subsystems.
- Constant species are not considered dependents since they are internally treated as
  parameters.
"""
function dependents(rx, network)
    if rx.rate isa Number
        return rx.substrates
    else
        rvars = get_variables(rx.rate, species(network))
        return union!(rvars, rx.substrates)
    end
end

"""
    dependents(rx, network)

See documentation for [`dependents`](@ref).
"""
function dependants(rx, network)
    dependents(rx, network)
end

"""
    reactionrates(network)

Given a [`ReactionSystem`](@ref), returns a vector of the symbolic reaction
rates for each reaction.
"""
function reactionrates(rn)
    [r.rate for r in reactions(rn)]
end

"""
    substoichmat(rn; sparse=false)

Returns the substrate stoichiometry matrix, ``S``, with ``S_{i j}`` the stoichiometric
coefficient of the ith substrate within the jth reaction.

Note:
- Set sparse=true for a sparse matrix representation
- Note that constant species are not considered substrates, but just components that modify
  the associated rate law.
"""
function substoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem)
    Is=Int[];  Js=Int[];  Vs=Int[];
    smap = speciesmap(rn)
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.substoich
        for (i,sub) in enumerate(rx.substrates)
            isconstant(sub) && continue
            push!(Js, k)
            push!(Is, smap[sub])
            push!(Vs, stoich[i])
        end
    end
    sparse(Is,Js,Vs,numspecies(rn),numreactions(rn))
end
function substoichmat(::Type{Matrix{Int}},rn::ReactionSystem)
    smap = speciesmap(rn)
    smat = zeros(Int, numspecies(rn), numreactions(rn))
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.substoich
        for (i,sub) in enumerate(rx.substrates)
            isconstant(sub) && continue
            smat[smap[sub],k] = stoich[i]
        end
    end
    smat
end
function substoichmat(rn::ReactionSystem; sparse::Bool=false)
    isempty(get_systems(rn)) || error("substoichmat does not currently support subsystems.")
	sparse ? substoichmat(SparseMatrixCSC{Int,Int}, rn) : substoichmat(Matrix{Int}, rn)
end

"""
    prodstoichmat(rn; sparse=false)

Returns the product stoichiometry matrix, ``P``, with ``P_{i j}`` the stoichiometric
coefficient of the ith product within the jth reaction.

Note:
- Set sparse=true for a sparse matrix representation
- Note that constant species are not treated as products, but just components that modify
  the associated rate law.
"""
function prodstoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem)
    Is=Int[];  Js=Int[];  Vs=Int[];
    smap = speciesmap(rn)
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.prodstoich
        for (i,prod) in enumerate(rx.products)
            isconstant(prod) && continue
			push!(Js, k)
			push!(Is, smap[prod])
			push!(Vs, stoich[i])
        end
    end
    sparse(Is,Js,Vs,numspecies(rn),numreactions(rn))
end
function prodstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem)
    smap = speciesmap(rn)
    pmat = zeros(Int, numspecies(rn), numreactions(rn))
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.prodstoich
        for (i,prod) in enumerate(rx.products)
            isconstant(prod) && continue
            pmat[smap[prod],k] = stoich[i]
        end
    end
    pmat
end
function prodstoichmat(rn::ReactionSystem; sparse=false)
    isempty(get_systems(rn)) || error("prodstoichmat does not currently support subsystems.")
	sparse ? prodstoichmat(SparseMatrixCSC{Int,Int}, rn) : prodstoichmat(Matrix{Int}, rn)
end

"""
    netstoichmat(rn, sparse=false)

Returns the net stoichiometry matrix, ``N``, with ``N_{i j}`` the net stoichiometric
coefficient of the ith species within the jth reaction.

Notes:
- Set sparse=true for a sparse matrix representation
- Caches the matrix internally within `rn` so subsequent calls are fast.
- Note that constant species are not treated as reactants, but just components that modify
  the associated rate law. As such they do not contribute to the net stoichiometry matrix.
"""
function netstoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem)
    Is=Int[];  Js=Int[];  Vs=Int[];
    smap = speciesmap(rn)
    for (k,rx) in pairs(reactions(rn))
        for (spec,coef) in rx.netstoich
            isconstant(spec) && continue
			push!(Js, k)
			push!(Is, smap[spec])
			push!(Vs, coef)
        end
    end
    sparse(Is,Js,Vs,numspecies(rn),numreactions(rn))
end
function netstoichmat(::Type{Matrix{Int}},rn::ReactionSystem)
    smap = speciesmap(rn)
    nmat = zeros(Int,numspecies(rn),numreactions(rn))
    for (k,rx) in pairs(reactions(rn))
        for (spec,coef) in rx.netstoich
            isconstant(spec) && continue
            nmat[smap[spec],k] = coef
        end
    end
    nmat
end
function netstoichmat(rn::ReactionSystem; sparse=false)
    isempty(get_systems(rn)) || error("netstoichmat does not currently support subsystems, please create a flattened system before calling.")

    nps = get_networkproperties(rn)

    # if it is already calculated and has the right type
    !isempty(nps.netstoichmat) && (sparse == issparse(nps.netstoichmat)) && (return nps.netstoichmat)

	if sparse
        nps.netstoichmat = netstoichmat(SparseMatrixCSC{Int,Int}, rn)
    else
        nps.netstoichmat = netstoichmat(Matrix{Int}, rn)
    end

    nps.netstoichmat
end

"""
    setdefaults!(rn, newdefs)

Sets the default (initial) values of parameters and species in the
`ReactionSystem`, `rn`.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
setdefaults!(sir, [:S => 999.0, :I => 1.0, :R => 1.0, :β => 1e-4, :ν => .01])

# or
@parameter β ν
@variables t S(t) I(t) R(t)
setdefaults!(sir, [S => 999.0, I => 1.0, R => 0.0, β => 1e-4, ν => .01])
```
gives initial/default values to each of `S`, `I` and `β`

Notes:
- Can not be used to set default values for species, variables or parameters of
  subsystems or constraint systems. Either set defaults for those systems
  directly, or [`flatten`](@ref) to collate them into one system before setting
  defaults.
- Defaults can be specified in any iterable container of symbols to value pairs
  or symbolics to value pairs.
"""
function setdefaults!(rn, newdefs)
    defs = eltype(newdefs) <: Pair{Symbol} ? symmap_to_varmap(rn,newdefs) : newdefs
    rndefs = MT.get_defaults(rn)
    for (var,val) in defs
        rndefs[var] = value(val)
    end
    nothing
end

function __unpacksys(rn)
    ex = :(begin end)
    for key in keys(get_var_to_name(rn))
        var = MT.getproperty(rn, key, namespace=false)
        push!(ex.args, :($key = $var))
    end
    ex
end

"""
    @unpacksys sys::ModelingToolkit.AbstractSystem

Loads all species, variables, parameters, and observables defined in `sys` as
variables within the calling module.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
@unpacksys sir
```
will load the symbolic variables, `S`, `I`, `R`, `ν` and `β`.

Notes:
- Can not be used to load species, variables, or parameters of subsystems or
  constraints. Either call `@unpacksys` on those systems directly, or
  [`flatten`](@ref) to collate them into one system before calling.
- Note that this places symbolic variables within the calling module's scope, so
  calling from a function defined in a script or the REPL will still result in
  the symbolic variables being defined in the `Main` module.
"""
macro unpacksys(rn)
    quote
    ex = Catalyst.__unpacksys($(esc(rn)))
    Base.eval($(__module__), ex)
    end
end

# convert symbol of the form :sys.a.b.c to a symbolic a.b.c
function _symbol_to_var(sys, sym)
    if hasproperty(sys, sym)
        var = getproperty(sys, sym, namespace=false)
    else
        strs = split(String(sym), "₊")   # need to check if this should be split of not!!!
        if length(strs) > 1
            var = getproperty(sys, Symbol(strs[1]), namespace=false)
            for str in view(strs, 2:length(strs))
                var = getproperty(var, Symbol(str), namespace=true)
            end
        else
            throw(ArgumentError("System $(nameof(sys)): variable $sym does not exist"))
        end
    end
    var
end

"""
    symmap_to_varmap(sys, symmap)

Given a system and map of `Symbol`s to values, generates a map from
corresponding symbolic variables/parameters to the values that can be used to
pass initial conditions and parameter mappings.

For example,
```julia
sir = @reaction_network sir begin
    β, S + I --> 2I
    ν, I --> R
end β ν
subsys = @reaction_network subsys begin
    k, A --> B
end k
@named sys = compose(sir, [subsys])
```
gives
```
Model sys with 3 equations
States (5):
  S(t)
  I(t)
  R(t)
  subsys₊A(t)
  subsys₊B(t)
Parameters (3):
  β
  ν
  subsys₊k
```
to specify initial condition and parameter mappings from *symbols* we can use
```julia
symmap = [:S => 1.0, :I => 1.0, :R => 1.0, :subsys₊A => 1.0, :subsys₊B => 1.0]
u0map  = symmap_to_varmap(sys, symmap)
pmap   = symmap_to_varmap(sys, [:β => 1.0, :ν => 1.0, :subsys₊k => 1.0])
```
`u0map` and `pmap` can then be used as input to various problem types.

Notes:
- Any `Symbol`, `sym`, within `symmap` must be a valid field of `sys`. i.e.
  `sys.sym` must be defined.
"""
function symmap_to_varmap(sys, symmap::Tuple)
    if all(p -> p isa Pair{Symbol}, symmap)
        return ((_symbol_to_var(sys,sym) => val for (sym,val) in symmap)...,)
    else  # if not all entries map a symbol to value pass through
        return symmap
    end
end

symmap_to_varmap(sys, symmap::AbstractArray{Pair{Symbol,T}}) where {T} =
    [_symbol_to_var(sys,sym) => val for (sym,val) in symmap]

symmap_to_varmap(sys, symmap::Dict{Symbol,T}) where {T} =
    Dict(_symbol_to_var(sys,sym) => val for (sym,val) in symmap)

# don't permute any other types and let varmap_to_vars handle erroring
symmap_to_varmap(sys, symmap) = symmap
#error("symmap_to_varmap requires a Dict, AbstractArray or Tuple to map Symbols to values.")

######################## reaction complexes and reaction rates ###############################

"""
    reset_networkproperties!(rn::ReactionSystem)

Clears the cache of various properties (like the netstoichiometry matrix). Use if such
properties need to be recalculated for some reason.
"""
function reset_networkproperties!(rn::ReactionSystem)
    reset!(get_networkproperties(rn))
    nothing
end

# get the species indices and stoichiometry while filtering out constant species.
function filter_constspecs(specs, stoich::AbstractVector{V}, smap) where {V <: Integer}
    isempty(specs) && (return Vector{Int}(),Vector{V}())

    if any(isconstant, specs)
        ids = Vector{Int}()
        filtered_stoich = Vector{V}()
        for (i,s) in enumerate(specs)
            if !isconstant(s)
                push!(ids, smap[s])
                push!(filtered_stoich, stoich[i])
            end
        end
    else
        ids = map(Base.Fix1(getindex, smap), specs)
        filtered_stoich = copy(stoich)
    end
    ids,filtered_stoich
end

"""
    reactioncomplexmap(rn::ReactionSystem)

Find each [`ReactionComplex`](@ref) within the specified system, constructing a mapping from
the complex to vectors that indicate which reactions it appears in as substrates and
products.

Notes:
- Each [`ReactionComplex`](@ref) is mapped to a vector of pairs, with each pair having the
  form `reactionidx => ± 1`, where `-1` indicates the complex appears as a substrate and
  `+1` as a product in the reaction with integer label `reactionidx`.
- Constant species are ignored as part of a complex. i.e. if species `A` is constant then
  the reaction `A + B --> C + D` is considered to consist of the complexes `B` and `C + D`.
  Likewise `A --> B` would be treated as the same as `0 --> B`.
"""
function reactioncomplexmap(rn::ReactionSystem)
    isempty(get_systems(rn)) || error("reactioncomplexmap does not currently support subsystems.")

    # check if previously calculated and hence cached
    nps = get_networkproperties(rn)
    !isempty(nps.complextorxsmap) && return nps.complextorxsmap
    complextorxsmap = nps.complextorxsmap

    rxs = reactions(rn)
    smap = speciesmap(rn)
    numreactions(rn) > 0 || error("There must be at least one reaction to find reaction complexes.")
    for (i,rx) in enumerate(rxs)
        subids,substoich = filter_constspecs(rx.substrates, rx.substoich, smap)
        subrc = sort!(ReactionComplex(subids, substoich))
        if haskey(complextorxsmap, subrc)
            push!(complextorxsmap[subrc], i => -1)
        else
            complextorxsmap[subrc] = [i => -1]
        end

        prodids,prodstoich = filter_constspecs(rx.products, rx.prodstoich, smap)
        prodrc = sort!(ReactionComplex(prodids, prodstoich))
        if haskey(complextorxsmap, prodrc)
            push!(complextorxsmap[prodrc], i => 1)
        else
            complextorxsmap[prodrc] = [i => 1]
        end
    end
    complextorxsmap
end

function reactioncomplexes(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem, complextorxsmap)
    complexes = collect(keys(complextorxsmap))
    Is=Int[];  Js=Int[];  Vs=Int[];
	for (i,c) in enumerate(complexes)
        for (j,σ) in complextorxsmap[c]
			push!(Is, i)
			push!(Js, j)
			push!(Vs, σ)
        end
    end
	B = sparse(Is,Js,Vs,length(complexes),numreactions(rn))
    complexes,B
end
function reactioncomplexes(::Type{Matrix{Int}}, rn::ReactionSystem, complextorxsmap)
    complexes = collect(keys(complextorxsmap))
    B = zeros(Int, length(complexes), numreactions(rn));
    for (i,c) in enumerate(complexes)
        for (j,σ) in complextorxsmap[c]
            B[i,j] = σ
        end
    end
    complexes,B
end

@doc raw"""
    reactioncomplexes(network::ReactionSystem; sparse=false)

Calculate the reaction complexes and complex incidence matrix for the given
[`ReactionSystem`](@ref).

Notes:
- returns a pair of a vector of [`ReactionComplex`](@ref)s and the complex incidence matrix.
- An empty [`ReactionComplex`](@ref) denotes the null (∅) state (from reactions like ∅ -> A
  or A -> ∅).
- Constant species are ignored in generating a reaction complex. i.e. if A is constant then
  A --> B consists of the complexes ∅ and B.
- The complex incidence matrix, ``B``, is number of complexes by number of reactions with
```math
B_{i j} = \begin{cases}
-1, &\text{if the i'th complex is the substrate of the j'th reaction},\\
1, &\text{if the i'th complex is the product of the j'th reaction},\\
0, &\text{otherwise.}
\end{cases}
```
- Set sparse=true for a sparse matrix representation of the incidence matrix
"""
function reactioncomplexes(rn::ReactionSystem; sparse=false)
    isempty(get_systems(rn)) || error("reactioncomplexes does not currently support subsystems.")
    nps = get_networkproperties(rn)
    if isempty(nps.complexes) || (sparse != issparse(nps.complexes))
        complextorxsmap = reactioncomplexmap(rn)
        nps.complexes,nps.incidencemat = if sparse
            reactioncomplexes(SparseMatrixCSC{Int,Int}, rn, complextorxsmap)
        else
            reactioncomplexes(Matrix{Int}, rn, complextorxsmap)
        end
    end
    nps.complexes,nps.incidencemat
end

"""
    incidencemat(rn::ReactionSystem; sparse=false)

Calculate the incidence matrix of `rn`, see [`reactioncomplexes`](@ref).

Notes:
- Is cached in `rn` so that future calls, assuming the same sparsity, will also be fast.
"""
incidencemat(rn::ReactionSystem; sparse=false) = reactioncomplexes(rn; sparse)[2]

function complexstoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem, rcs)
    Is=Int[];  Js=Int[];  Vs=Int[];
    for (i,rc) in enumerate(rcs)
        for rcel in rc
			push!(Is, rcel.speciesid)
			push!(Js, i)
			push!(Vs, rcel.speciesstoich)
        end
    end
    Z = sparse(Is,Js,Vs, numspecies(rn), length(rcs))
end
function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem, rcs)
    Z=zeros(Int, numspecies(rn), length(rcs))
    for (i,rc) in enumerate(rcs)
        for rcel in rc
            Z[rcel.speciesid,i] = rcel.speciesstoich
        end
    end
    Z
end

"""
    complexstoichmat(network::ReactionSystem; sparse=false)

Given a [`ReactionSystem`](@ref) and vector of reaction complexes, return a
matrix with positive entries of size number of species by number of complexes,
where the non-zero positive entries in the kth column denote stoichiometric
coefficients of the species participating in the kth reaction complex.

Notes:
- Set sparse=true for a sparse matrix representation
"""
function complexstoichmat(rn::ReactionSystem; sparse=false)
    isempty(get_systems(rn)) || error("complexstoichmat does not currently support subsystems.")
    nps = get_networkproperties(rn)
    if isempty(nps.complexstoichmat) || (sparse != issparse(nps.complexstoichmat))
        nps.complexstoichmat = if sparse
            complexstoichmat(SparseMatrixCSC{Int,Int}, rn, keys(reactioncomplexmap(rn)))
        else
            complexstoichmat(Matrix{Int}, rn, keys(reactioncomplexmap(rn)))
        end
    end
    nps.complexstoichmat
end

function complexoutgoingmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem, B)
    n = size(B,2)
	rows = rowvals(B)
	vals = nonzeros(B)
    Is = Int[]; Js = Int[]; Vs = Int[]
    sizehint!(Is, div(length(vals),2))
    sizehint!(Js, div(length(vals),2))
    sizehint!(Vs, div(length(vals),2))
	for j = 1:n
	   for i in nzrange(B, j)
	      if vals[i] != one(eltype(vals))
            push!(Is, rows[i])
            push!(Js, j)
            push!(Vs, vals[i])
          end
	   end
	end
    sparse(Is,Js,Vs,size(B,1),size(B,2))
end
function complexoutgoingmat(::Type{Matrix{Int}}, rn::ReactionSystem, B)
    Δ = copy(B)
    for (I,b) in pairs(Δ)
        (b == 1) && (Δ[I] = 0)
    end
    Δ
end

@doc raw"""
    complexoutgoingmat(network::ReactionSystem; sparse=false)

Given a [`ReactionSystem`](@ref) and complex incidence matrix, ``B``, return a
matrix of size num of complexes by num of reactions that identifies substrate
complexes.

Notes:
- The complex outgoing matrix, ``\Delta``, is defined by
```math
\Delta_{i j} = \begin{cases}
    = 0,    &\text{if } B_{i j} = 1, \\
    = B_{i j}, &\text{otherwise.}
\end{cases}
```
- Set sparse=true for a sparse matrix representation
"""
function complexoutgoingmat(rn::ReactionSystem; sparse=false)
    isempty(get_systems(rn)) || error("complexoutgoingmat does not currently support subsystems.")
    nps = get_networkproperties(rn)
    if isempty(nps.complexoutgoingmat) || (sparse != issparse(nps.complexoutgoingmat))
        B = reactioncomplexes(rn, sparse=sparse)[2]
        nps.complexoutgoingmat = if sparse
            complexoutgoingmat(SparseMatrixCSC{Int,Int}, rn, B)
        else
            complexoutgoingmat(Matrix{Int}, rn, B)
        end
    end
    nps.complexoutgoingmat
end

function incidencematgraph(incidencemat::Matrix{Int})
    @assert all(∈([-1,0,1]) ,incidencemat)
    n = size(incidencemat,1)  # no. of nodes/complexes
    graph = Graphs.DiGraph(n)
    for col in eachcol(incidencemat)
        src = 0; dst = 0;
        for i in eachindex(col)
                (col[i] == -1) && (src = i)
                (col[i] == 1) && (dst = i)
                (src != 0) && (dst != 0) && break
        end
        Graphs.add_edge!(graph, src, dst)
    end
    return graph
end
function incidencematgraph(incidencemat::SparseMatrixCSC{Int,Int})
    @assert all(∈([-1,0,1]) ,incidencemat)
    m,n = size(incidencemat)
    graph = Graphs.DiGraph(m)
    rows = rowvals(incidencemat)
    vals = nonzeros(incidencemat)
    for j = 1:n
        inds=nzrange(incidencemat, j)
        row = rows[inds];
        val = vals[inds];
        if val[1] == -1
            Graphs.add_edge!(graph, row[1], row[2])
        else
            Graphs.add_edge!(graph, row[2], row[1])
        end
    end
    return graph
end

"""
    incidencematgraph(rn::ReactionSystem)

Construct a directed simple graph where nodes correspond to reaction complexes and directed
edges to reactions converting between two complexes.

Notes:
- Requires the `incidencemat` to already be cached in `rn` by a previous call to
  `reactioncomplexes`.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
complexes,incidencemat = reactioncomplexes(sir)
incidencematgraph(sir)
```
"""
function incidencematgraph(rn::ReactionSystem)
    nps = get_networkproperties(rn)
    if Graphs.nv(nps.incidencegraph) == 0
        isempty(nps.incidencemat) && error("Please call reactioncomplexes(rn) first to "
                                           * "construct the incidence matrix.")
        nps.incidencegraph = incidencematgraph(nps.incidencemat)
    end
    nps.incidencegraph
end

linkageclasses(incidencegraph) = Graphs.connected_components(incidencegraph)

"""
    linkageclasses(rn::ReactionSystem)

Given the incidence graph of a reaction network, return a vector of the
connected components of the graph (i.e. sub-groups of reaction complexes that
are connected in the incidence graph).

Notes:
- Requires the `incidencemat` to already be cached in `rn` by a previous call to
  `reactioncomplexes`.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
complexes,incidencemat = reactioncomplexes(sir)
linkageclasses(sir)
```
gives
```julia
2-element Vector{Vector{Int64}}:
 [1, 2]
 [3, 4]
```
"""
function linkageclasses(rn::ReactionSystem)
    nps = get_networkproperties(rn)
    if isempty(nps.linkageclasses)
        nps.linkageclasses = linkageclasses(incidencematgraph(rn))
    end
    nps.linkageclasses
end

@doc raw"""
    deficiency(rn::ReactionSystem)

Calculate the deficiency of a reaction network.

Here the deficiency, ``\delta``, of a network with ``n`` reaction complexes,
``\ell`` linkage classes and a rank ``s`` stoichiometric matrix is

```math
\delta = n - \ell - s
```

Notes:
- Requires the `incidencemat` to already be cached in `rn` by a previous call to
  `reactioncomplexes`.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
δ = deficiency(sir)
```
"""
function deficiency(rn::ReactionSystem)
    nps = get_networkproperties(rn)
    conservationlaws(rn)
    r = nps.rank
    ig = incidencematgraph(rn)
    lc = linkageclasses(rn)
    nps.deficiency = Graphs.nv(ig) - length(lc) - r
    nps.deficiency
end

function subnetworkmapping(linkageclass, allrxs, complextorxsmap, p)
    rxinds  = sort!(collect(Set(rxidx for rcidx in linkageclass for rxidx in complextorxsmap[rcidx])))
    rxs     = allrxs[rxinds]
    specset = Set(s for rx in rxs for s in rx.substrates if !isconstant(s))
    for rx in rxs
        for product in rx.products
            !isconstant(product) && push!(specset, product)
        end
    end
    specs = collect(specset)
    newps = Vector{eltype(p)}()
    for rx in rxs
        Symbolics.get_variables!(newps, rx.rate, p)
    end
    rxs, specs, newps   # reactions and species involved in reactions of subnetwork
end

"""
    subnetworks(rn::ReactionSystem)

Find subnetworks corresponding to each linkage class of the reaction network.

Notes:
- Requires the `incidencemat` to already be cached in `rn` by a previous call to
  `reactioncomplexes`.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
complexes,incidencemat = reactioncomplexes(sir)
subnetworks(sir)
```
"""
function subnetworks(rs::ReactionSystem)
    isempty(get_systems(rs)) || error("subnetworks does not currently support subsystems.")
    lcs = linkageclasses(rs)
    rxs = reactions(rs)
    p = parameters(rs)
    t = get_iv(rs)
    complextorxsmap = [map(first,rcmap) for rcmap in values(reactioncomplexmap(rs))]
    subnetworks = Vector{ReactionSystem}()
    for i in 1:length(lcs)
        reacs,specs,newps = subnetworkmapping(lcs[i], rxs, complextorxsmap, p)
        newname = Symbol(nameof(rs), "_", i)
        push!(subnetworks, ReactionSystem(reacs, t, specs, newps; name=newname))
    end
    subnetworks
end

"""
    linkagedeficiencies(network::ReactionSystem)

Calculates the deficiency of each sub-reaction network within `network`.

Notes:
- Requires the `incidencemat` to already be cached in `rn` by a previous call to
  `reactioncomplexes`.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
linkage_deficiencies = linkagedeficiencies(sir)
```
"""
function linkagedeficiencies(rs::ReactionSystem)
    lcs = linkageclasses(rs)
    subnets = subnetworks(rs)
    δ = zeros(Int,length(lcs))
    for (i,subnet) in enumerate(subnets)
        conservationlaws(subnet)
        nps = get_networkproperties(subnet)
        δ[i] = length(lcs[i]) - 1 - nps.rank
    end
    δ
end

"""
    isreversible(rn::ReactionSystem)

Given a reaction network, returns if the network is reversible or not.

Notes:
- Requires the `incidencemat` to already be cached in `rn` by a previous call to
  `reactioncomplexes`.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
isreversible(sir)
```
"""
function isreversible(rn::ReactionSystem)
    ig = incidencematgraph(rn)
    Graphs.reverse(ig) == ig
end

"""
    isweaklyreversible(rn::ReactionSystem, subnetworks)

Determine if the reaction network with the given subnetworks is weakly reversible or not.

Notes:
- Requires the `incidencemat` to already be cached in `rn` by a previous call to
  `reactioncomplexes`.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
subnets = subnetworks(rn)
isweaklyreversible(rn, subnets)
```
"""
function isweaklyreversible(rn::ReactionSystem, subnets)
    im = get_networkproperties(rn).incidencemat
    isempty(im) && error("Error, please call reactioncomplexes(rn::ReactionSystem) to "
                         * "ensure the incidence matrix has been cached.")
    sparseig = issparse(im)
    for subnet in subnets
        nps = get_networkproperties(subnet)
        isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse=sparseig)
    end
    all(Graphs.is_strongly_connected ∘ incidencematgraph, subnets)
end

############################################################################################
######################## conservation laws ###############################

"""
    conservedequations(rn::ReactionSystem)

Calculate symbolic equations from conservation laws, writing dependent variables as
functions of independent variables and the conservation law constants.

Notes:
- Caches the resulting equations in `rn`, so will be fast on subsequent calls.

Examples:
```@repl
rn = @reaction_network begin
    k, A + B --> C
    k2, C --> A + B
    end k k2
conservedequations(rn)
```
gives
```
2-element Vector{Equation}:
 B(t) ~ A(t) + _ConLaw[1]
 C(t) ~ _ConLaw[2] - A(t)
```
"""
function conservedequations(rn::ReactionSystem)
    conservationlaws(rn)
    nps = get_networkproperties(rn)
    nps.conservedeqs
end

"""
    conservationlaw_constants(rn::ReactionSystem)

Calculate symbolic equations from conservation laws, writing the conservation law constants
in terms of the dependent and independent variables.

Notes:
- Caches the resulting equations in `rn`, so will be fast on subsequent calls.

Examples:
```@julia
rn = @reaction_network begin
    k, A + B --> C
    k2, C --> A + B
    end k k2
conservationlaw_constants(rn)
```
gives
```
2-element Vector{Equation}:
 _ConLaw[1] ~ B(t) - A(t)
 _ConLaw[2] ~ A(t) + C(t)
```
"""
function conservationlaw_constants(rn::ReactionSystem)
    conservationlaws(rn)
    nps = get_networkproperties(rn)
    nps.constantdefs
end

"""
    conservationlaws(netstoichmat::AbstractMatrix)::Matrix

Given the net stoichiometry matrix of a reaction system, computes a matrix of
conservation laws, each represented as a row in the output.
"""
function conservationlaws(nsm::T; col_order=nothing) where {T <: AbstractMatrix}

    # compute the left nullspace over the integers
    N = MT.nullspace(nsm'; col_order)

    # if all coefficients for a conservation law are negative, make positive
    for Nrow in eachcol(N)
        all(r -> r <= 0, Nrow) && (Nrow .*= -1)
    end

    # check we haven't overflowed
    iszero(N' * nsm) || error("Calculation of the conservation law matrix was inaccurate, "
                            * "likely due to numerical overflow. Please use a larger integer "
                            * "type like Int128 or BigInt for the net stoichiometry matrix.")

    T(N')
end

function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_order)
    nullity = size(N,1)
    r = numspecies(rn) - nullity     # rank of the netstoichmat
    sts = species(rn)
    indepidxs = col_order[begin:r]
    indepspecs = sts[indepidxs]
    depidxs = col_order[(r+1):end]
    depspecs = sts[depidxs]
    constants = MT.unwrap.(MT.scalarize((@parameters _ConLaw[1:nullity])[1]))

    conservedeqs = Equation[]
    constantdefs = Equation[]
    for (i,depidx) in enumerate(depidxs)
        scaleby = (N[i,depidx] != 1) ? N[i,depidx] : one(eltype(N))
        (scaleby != 0) || error("Error, found a zero in the conservation law matrix where "
                                 * "one was not expected.")
        coefs = @view N[i,indepidxs]
        terms = sum(p -> p[1]/scaleby * p[2], zip(coefs,indepspecs))
        eq = depspecs[i] ~ constants[i] - terms
        push!(conservedeqs, eq)
        eq = constants[i] ~ depspecs[i] + terms
        push!(constantdefs,eq)
    end

    # cache in the system
    nps = get_networkproperties(rn)
    nps.rank = r
    nps.nullity = nullity
    nps.indepspecs = Set(indepspecs)
    nps.depspecs = Set(depspecs)
    nps.conservedeqs = conservedeqs
    nps.constantdefs = constantdefs

    nothing
end

"""
    conservationlaws(rs::ReactionSystem)

Return the conservation law matrix of the given `ReactionSystem`, calculating it if it is
not already stored within the system, or returning an alias to it.

Notes:
- The first time being called it is calculated and cached in `rn`, subsequent calls should
  be fast.
"""
function conservationlaws(rs::ReactionSystem)
    nps = get_networkproperties(rs)
    !isempty(nps.conservationmat) && (return nps.conservationmat)
    nsm = netstoichmat(rs)
    nps.conservationmat = conservationlaws(nsm; col_order=nps.col_order)
    cache_conservationlaw_eqs!(rs, nps.conservationmat, nps.col_order)
    nps.conservationmat
end

"""
    conservedquantities(state, cons_laws)

Compute conserved quantities for a system with the given conservation laws.
"""
conservedquantities(state, cons_laws) = cons_laws * state

######################## reaction network operators #######################

"""
    ==(rx1::Reaction, rx2::Reaction)

Tests whether two [`Reaction`](@ref)s are identical.

Notes:
- Ignores the order in which stoichiometry components are listed.
- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
    considered different than `(A+1)^2`.
"""
function (==)(rx1::Reaction, rx2::Reaction)
    isequal(rx1.rate, rx2.rate) || return false
    issetequal(zip(rx1.substrates,rx1.substoich), zip(rx2.substrates,rx2.substoich)) || return false
    issetequal(zip(rx1.products,rx1.prodstoich), zip(rx2.products,rx2.prodstoich)) || return false
    issetequal(rx1.netstoich, rx2.netstoich) || return false
    rx1.only_use_rate == rx2.only_use_rate
end

function hash(rx::Reaction, h::UInt)
    h = Base.hash(rx.rate, h)
    h = Base.hash(rx.substrates, h)
    h = Base.hash(rx.products, h)
    h = Base.hash(rx.prodstoich, h)
    h = Base.hash(rx.substoich, h)
    h = Base.hash(rx.netstoich, h)
    Base.hash(rx.only_use_rate, h)
end

"""
    isequal_ignore_names(rn1::ReactionSystem, rn2::ReactionSystem)

Tests whether the underlying species, parameters and reactions are the same in
the two [`ReactionSystem`](@ref)s. Ignores the names of the systems in testing
equality.

Notes:
- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
    considered different than `(A+1)^2`.
- Does not include `defaults` in determining equality.
"""
function isequal_ignore_names(rn1::ReactionSystem, rn2::ReactionSystem)
    (get_combinatoric_ratelaws(rn1) == get_combinatoric_ratelaws(rn2)) || return false
    isequal(get_iv(rn1), get_iv(rn2)) || return false
    issetequal(get_states(rn1), get_states(rn2)) || return false
    issetequal(get_ps(rn1), get_ps(rn2)) || return false
    issetequal(MT.get_observed(rn1), MT.get_observed(rn2)) || return false
    issetequal(get_eqs(rn1), get_eqs(rn2)) || return false
    (get_constraints(rn1) == get_constraints(rn2)) || return false

    # subsystems
    (length(get_systems(rn1)) == length(get_systems(rn2))) || return false
    issetequal(get_systems(rn1), get_systems(rn2)) || return false

    true
end

"""
    ==(rn1::ReactionSystem, rn2::ReactionSystem)

Tests whether the underlying species, parameters and reactions are the same in
the two [`ReactionSystem`](@ref)s. Requires the systems to have the same names
too.

Notes:
- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
    considered different than `(A+1)^2`.
- Does not include `defaults` in determining equality.
"""
function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
    (nameof(rn1) == nameof(rn2)) || return false
    isequal_ignore_names(rn1,rn2)
end

######################## functions to extend a network ####################

"""
    make_empty_network(; iv=DEFAULT_IV, name=gensym(:ReactionSystem))

Construct an empty [`ReactionSystem`](@ref). `iv` is the independent variable,
usually time, and `name` is the name to give the `ReactionSystem`.
"""
function make_empty_network(; iv=DEFAULT_IV, name=gensym(:ReactionSystem))
    ReactionSystem(Reaction[], iv, [], []; name=name)
end

"""
    addspecies!(network::ReactionSystem, s::Symbolic; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the species corresponding to the variable
`s` to the network (if it is not already defined). Returns the integer id of the
species within the system.

Notes:
- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addspecies!(network::ReactionSystem, s::Symbolic; disablechecks=false)
    reset_networkproperties!(network)

    isconstant(s) && error("Constant species should be added via addparams!.")

    # we don't check subsystems since we will add it to the top-level system...
    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, s), get_states(network))
    if curidx === nothing
        push!(get_states(network), s)
        MT.process_variables!(get_var_to_name(network), get_defaults(network), [s])
        return length(get_states(network))
    else
        return curidx
    end
end

"""
    addspecies!(network::ReactionSystem, s::Num; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the species corresponding to the
variable `s` to the network (if it is not already defined). Returns the
integer id of the species within the system.

- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addspecies!(network::ReactionSystem, s::Num; disablechecks=false)
    addspecies!(network, value(s), disablechecks=disablechecks)
end

"""
    reorder_states!(rn, neworder)

Given a [`ReactionSystem`](@ref) and a vector `neworder`, orders the states of `rn` accordingly to `neworder`.

Notes:
- Currently only supports `ReactionSystem`s without constraints or subsystems.
"""
function reorder_states!(rn, neworder)
   (get_constraints(rn) === nothing) && isempty(get_systems(rn)) ||
           error("Reordering of states is only supported for systems without constraints or subsystems.")
    permute!(get_states(rn), neworder)
    reset_networkproperties!(rn)
    nothing
end

"""
    addparam!(network::ReactionSystem, p::Symbolic; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the parameter corresponding to the
variable `p` to the network (if it is not already defined). Returns the integer
id of the parameter within the system.

- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addparam!(network::ReactionSystem, p::Symbolic; disablechecks=false)
    reset_networkproperties!(network)

    # we don't check subsystems since we will add it to the top-level system...
    if istree(p) && !(operation(p) isa Sym)
        error("If the passed in parameter is an expression, it must correspond to an underlying Variable.")
    end
    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, p), get_ps(network))
    if curidx === nothing
        push!(get_ps(network), p)
        MT.process_variables!(get_var_to_name(network), get_defaults(network), [p])
        return length(get_ps(network))
    else
        return curidx
    end
end

"""
    addparam!(network::ReactionSystem, p::Num; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the parameter corresponding to the
variable `p` to the network (if it is not already defined). Returns the
integer id of the parameter within the system.

- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addparam!(network::ReactionSystem, p::Num; disablechecks=false)
    addparam!(network, value(p); disablechecks=disablechecks)
end

"""
    addreaction!(network::ReactionSystem, rx::Reaction)

Add the passed in reaction to the [`ReactionSystem`](@ref). Returns the
integer id of `rx` in the list of `Reaction`s within `network`.

Notes:
- Any new species or parameters used in `rx` should be separately added to
    `network` using [`addspecies!`](@ref) and [`addparam!`](@ref).
"""
function addreaction!(network::ReactionSystem, rx::Reaction)
    reset_networkproperties!(network)
    push!(get_eqs(network), rx)
    length(get_eqs(network))
end

"""
    merge!(network1::ReactionSystem, network2::ReactionSystem)

Merge `network2` into `network1`.

Notes:
- Duplicate reactions between the two networks are not filtered out.
- [`Reaction`](@ref)s are not deepcopied to minimize allocations, so both
  networks will share underlying data arrays.
- Subsystems are not deepcopied between the two networks and will hence be
  shared.
- Returns `network1`.
- `combinatoric_ratelaws` is the value of `network1`.
"""
function Base.merge!(network1::ReactionSystem, network2::ReactionSystem)
    ((get_constraints(network1) === nothing) && (get_constraints(network2) === nothing)) ||
        error("merge! does not currently support ReactionSystems with constraints, consider ModelingToolkit.extend instead.")
    isequal(get_iv(network1), get_iv(network2)) ||
        error("Reaction networks must have the same independent variable to be mergable.")
    append!(get_eqs(network1), get_eqs(network2))
    union!(get_states(network1), get_states(network2))
    union!(get_ps(network1), get_ps(network2))
    append!(get_observed(network1), get_observed(network2))
    append!(get_systems(network1), get_systems(network2))
    merge!(get_defaults(network1), get_defaults(network2))
    reset_networkproperties!(network1)
    network1
end

###############################   units   #####################################

"""
    validate(rx::Reaction; info::String = "")

Check that all substrates and products within the given [`Reaction`](@ref) have
the same units, and that the units of the reaction's rate expression are
internally consistent (i.e. if the rate involves sums, each term in the sum has
the same units).

"""
function validate(rx::Reaction; info::String = "")
    validated = MT._validate([rx.rate], [string(rx, ": rate")], info = info)

    subunits = isempty(rx.substrates) ? nothing : get_unit(rx.substrates[1])
    for i in 2:length(rx.substrates)
        if get_unit(rx.substrates[i]) != subunits
            validated = false
            @warn(string("In ", rx, " the substrates have differing units."))
        end
    end

    produnits = isempty(rx.products) ? nothing : get_unit(rx.products[1])
    for i in 2:length(rx.products)
        if get_unit(rx.products[i]) != produnits
            validated = false
            @warn(string("In ", rx, " the products have differing units."))
        end
    end

    if (subunits !== nothing) && (produnits !== nothing) && (subunits != produnits)
        validated = false
        @warn(string("in ", rx, " the substrate units are not consistent with the product units."))
    end

    validated
end

"""
    validate(rs::ReactionSystem, info::String="")

Check that all species in the [`ReactionSystem`](@ref) have the same units, and
that the rate laws of all reactions reduce to units of (species units) / (time
units).

Notes:
- Does not check subsystems too.
"""
function validate(rs::ReactionSystem, info::String="")
    specs = get_states(rs)

    # if there are no species we don't check units on the system
    isempty(specs) && return true

    specunits = get_unit(specs[1])
    validated = true
    for spec in specs
        if get_unit(spec) != specunits
            validated = false
            @warn(string("Species are expected to have units of ", specunits, " however, species ", spec, " has units ", get_unit(spec), "."))
        end
    end
    timeunits = get_unit(get_iv(rs))

    # no units for species, time or parameters then assume validated
    (specunits in (MT.unitless,nothing)) && (timeunits in (MT.unitless,nothing)) &&
        MT.all_dimensionless(get_ps(rs)) && return true

    rateunits = specunits / timeunits
    for rx in get_eqs(rs)
        rxunits = get_unit(rx.rate)
        for (i,sub) in enumerate(rx.substrates)
            rxunits *= get_unit(sub) ^ rx.substoich[i]
        end

        if rxunits != rateunits
            validated = false
            @warn(string("Reaction rate laws are expected to have units of ", rateunits, " however, ", rx, " has units of ", rxunits, "."))
        end
    end

    validated
end


text:
# Julia wrapper for header: /usr/include/scip/scipdefplugins.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeDefaultPlugins(scip)
    ccall((:SCIPincludeDefaultPlugins, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

ext: jl


text:
push!(LOAD_PATH, joinpath(pwd(), "src"))
using JuML
using Test
using MLJ
import MLJBase
import DataFrames
import CSV

# train_df = DataFrame(joinpath("data", "airlinetrain"))
# test_df = DataFrame(joinpath("data", "airlinetest"))

# distance = train_df[:Distance]
# deptime = train_df[:DepTime]

# deptime = factor(train_df[:DepTime])
# distance = factor(train_df[:Distance])

# factors = [train_df.factors; [deptime, distance]]

# label = train_df[:dep_delayed_15min]

# model = XGBClassifier(maxdepth = 2, nrounds = 5)

# mach = machine(model, factors, label)
# fit!(mach; force = true)

# yhat = MLJ.predict(mach, train_df)

df = CSV.read(abspath(joinpath("data", "airlinetrain1M.csv")))

DataFrames.names(df)
DataFrames.categorical!(df, [1, 2, 3, 5, 6, 7, 9])

train_df = JuML.DataFrame(df)

distance = train_df["Distance"]
deptime = train_df["DepTime"]
label = covariate(train_df["dep_delayed_15min"], level -> level == "Y" ? 1.0 : 0.0)
deptime = factor(train_df["DepTime"])
distance = factor(train_df["Distance"])

factors = [train_df.factors; [deptime, distance]]

@time model7 = xgblogit(label, factors; η = 0.1, λ = 1.0, γ = 0.0, minchildweight = 1.0, nrounds = 100, maxdepth = 8, ordstumps = false, pruning = false, caching = true, usefloat64 = false, singlethread = true, slicelength = 0);
auc = getauc(model7.pred, label)
model7.pred[1:5]

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
"""
Convenience code for plotting result of DCPCA

Author: Eric Ford
Created: August 2020
"""

function plot_basis_vectors(λ::VR1, f_mean::V1, deriv::V2, proj::A3;
                              num_basis::Integer = 4, idx_plt::AR = 1:length(f_mean),
                              label = 0 ) where {
                              R1<:AbstractRange, VR1<:AbstractVector{R1}, T1<:Real, V1<:AbstractVector{T1},
                              T2<:Real, V2<:AbstractVector{T2}, T3<:Real, A3<:AbstractArray{T3,2}, AR<:AbstractUnitRange }
  @assert num_basis >=1
  @assert sum(length.(λ)) == length(f_mean)
  @assert length(f_mean) == length(deriv)
  @assert size(proj,1) == length(f_mean)
  @assert minimum(idx_plt) >= 1
  @assert maximum(idx_plt) <= length(f_mean)
  if num_basis > 4   @warn "plot_basis_vectors hasn't been generalized for num_basis!=4."   end
  λ_plt = get_λs(λ,idx_plt)
  λc = mean(λ_plt)
  plt0 = plot(λ_plt, (f_mean[idx_plt].-1.0)./maximum(abs.(f_mean[idx_plt].-1.0)),linecolor=:black, label=:none) #"Std Mean")
  plt0 = plot!(plt0,λ_plt, deriv[idx_plt]./maximum(abs.(deriv[idx_plt])),linecolor=:red,label=:none) #"Std Deriv")
  plt0 = plot!(plt0,[λc,λc],[-1.1,1.0], linecolor=:black, label=:none)
  ylims!(-1.1,1.1)
  pltall = [plt0]
  if num_basis >= 1
      plt0 = plot!(plt0,λ_plt, proj[idx_plt,1]./maximum(abs.(proj[idx_plt,1])),linecolor=:blue,label=:none)
      plt1 = plot(λ_plt, proj[idx_plt,1]./maximum(abs.(proj[idx_plt,1])),linecolor=:blue,label=:none) #"PC 1")
      plt1 = plot!(plt1,[λc,λc],[-1.1,1.0], linecolor=:black, label=:none)
      ylabel!(plt1,"PC1")
      ylims!(-1.1,1.1)
      push!(pltall,plt1)
  end
  if num_basis >= 2
      plt0 = plot!(plt0,λ_plt, proj[idx_plt,2]./maximum(abs.(proj[idx_plt,2])),linecolor=:green,label=:none)
      plt2 = plot(λ_plt, proj[idx_plt,2]./maximum(abs.(proj[idx_plt,2])),linecolor=:green,label=:none) #"PC 2")
      plt2 = plot!(plt2,[λc,λc],[-1.1,1.0], linecolor=:black, label=:none)
      ylabel!(plt2,"PC2")
      ylims!(-1.1,1.1)
      push!(pltall,plt2)
  end
  if num_basis >= 3
      plt3 = plot(λ_plt, proj[idx_plt,3]./maximum(abs.(proj[idx_plt,3])),linecolor=:cyan,label=:none) #"PC 3")
      plt3 = plot!(plt3,[λc,λc],[-1.1,1.0], linecolor=:black, label=:none)
      ylabel!(plt3,"PC3")
      ylims!(-1.1,1.1)
      push!(pltall,plt3)
  end
  if num_basis >= 4
      plt4 = plot(λ_plt, proj[idx_plt,4]./maximum(abs.(proj[idx_plt,4])),linecolor=:magenta,label=:none) #"PC 4")
      plt4 = plot!(plt4,[λc,λc],[-1.1,1.0], linecolor=:black, label=:none)
      ylabel!(plt4,"PC4")
      ylims!(-1.1,1.1)
      push!(pltall,plt4)
  end
  if num_basis == 1
      xlabel!(plt1,"λ (Å)  line# "  * string(label))
      pltall = plot(plt0,plt1, layout = (2,1) )
  elseif num_basis == 2
      xlabel!(plt2,"λ (Å)  line# "  * string(label))
      pltall = plot(plt0,plt1,plt2, layout = (3,1) )
  elseif num_basis == 3
      xlabel!(plt3,"λ (Å)  line# "  * string(label))
      pltall = plot(plt0,plt1,plt2,plt3, layout = (4,1) )
  elseif num_basis == 4
      xlabel!(plt4,"λ (Å)  line# "  * string(label))
      pltall = plot(plt0,plt1,plt2,plt3,plt4, layout = (5,1) )
  elseif num_basis == 5
      xlabel!(plt5,"λ (Å)  line# "  * string(label))
      pltall = plot(plt0,plt1,plt2,plt3,plt4,plt5, layout = (6,1) )
  end

  pltall
end

function plot_basis_scores(times::V1, rvs::V2, scores::A3;
                              num_basis::Integer = 4, idx_plt::AR = 1:length(times) ) where {
                              T1<:Real, V1<:AbstractVector{T1}, T2<:Real, V2<:AbstractVector{T2}, T3<:Real, A3<:AbstractArray{T3,2}, AR<:AbstractUnitRange }
  @assert length(times) >=1
  @assert length(times) == length(rvs)
  @assert length(times) == size(scores,2)
  @assert size(scores,1) >= num_basis
  @assert minimum(idx_plt) >= 1
  @assert maximum(idx_plt) <= length(times)
  if num_basis > 4   @warn "plot_basis_scores hasn't been generalized for num_basis>4."   end
  local plt0 = scatter(times[idx_plt],rvs[idx_plt],label=:none,color=:black)
  ylabel!(plt0,"RV (m/s)")
  if num_basis >= 1
      plt1 = scatter(times[idx_plt],scores[1,idx_plt],label=:none,color=:blue)
      pltall = plot(plt0,plt1, layout = (2,1) )
  end
  if num_basis >= 2
      plt2 = scatter(times[idx_plt],scores[2,idx_plt],label=:none,color=:green)
      pltall = plot(plt0,plt1,plt2, layout = (3,1) )
  end
  if num_basis >= 3
      plt3 = scatter(times[idx_plt],scores[3,idx_plt],label=:none,color=:cyan)
      pltall = plot(plt0,plt1,plt2,plt3, layout = (4,1) )
  end
  if num_basis >= 4
      plt4 = scatter(times[idx_plt],scores[4,idx_plt],label=:none,color=:magenta)
      pltall = plot(plt0,plt1,plt2,plt3,plt4, layout = (5,1) )
  end
  xlabel!(pltall,"Time")
  pltall
end

function plot_basis_scores_cor(rvs::V2, scores::A3;
                              num_basis::Integer = 4, idx_plt::AR = 1:length(rvs) ) where {
                              T2<:Real, V2<:AbstractVector{T2}, T3<:Real, A3<:AbstractArray{T3,2}, AR<:AbstractUnitRange }
  @assert length(rvs) >=1
  @assert length(rvs) == size(scores,2)
  @assert size(scores,1) >= num_basis
  @assert minimum(idx_plt) >= 1
  @assert maximum(idx_plt) <= length(rvs)
  if num_basis > 4   @warn "plot_basis_scores_cor hasn't been generalized for num_basis>4."   end
  if num_basis >= 1
      plt1 = scatter(rvs[idx_plt],vec(scores[1,idx_plt]),xlabel="RV",ylabel="PC1",legend=:none)
      pltall = plt1
  end
  if num_basis >= 2
      plt2 = scatter(rvs[idx_plt],vec(scores[2,idx_plt]),xlabel="RV",ylabel="PC2",legend=:none)
      plt4 = scatter(vec(scores[1,idx_plt]),vec(scores[2,idx_plt]),xlabel="PC1",ylabel="PC2",legend=:none)
      pltall = plot(plt1,plt2,plt4,layout=3)
  end
  if num_basis >= 3
      plt3 = scatter(rvs[idx_plt],vec(scores[3,idx_plt]),xlabel="RV",ylabel="PC3",legend=:none)
      plt5 = scatter(vec(scores[1,idx_plt]),vec(scores[3,idx_plt]),xlabel="PC1",ylabel="PC3",legend=:none)
      plt6 = scatter(vec(scores[2,idx_plt]),vec(scores[3,idx_plt]),xlabel="PC2",ylabel="PC3",legend=:none)
      pltall = plot(plt1,plt2,plt3,plt4,plt5,plt6,layout=(2,3))
  end
  if num_basis >= 4
      plt7 = scatter(vec(scores[1,idx_plt]),vec(scores[4,idx_plt]),xlabel="PC2",ylabel="PC4",legend=:none)
      plt8 = scatter(vec(scores[2,idx_plt]),vec(scores[4,idx_plt]),xlabel="PC2",ylabel="PC4",legend=:none)
      plt9 = scatter(vec(scores[3,idx_plt]),vec(scores[4,idx_plt]),xlabel="PC3",ylabel="PC4",legend=:none)
      pltall = plot(plt1,plt2,plt3,plt4,plt5,plt6,plt7,plt8,plt9,layout=(3,3))
  end
  pltall
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct LoadBalancerBackendAddress2 <: SwaggerModel
    properties::Any # spec type: Union{ Nothing, LoadBalancerBackendAddressPropertiesFormat2 } # spec name: properties
    name::Any # spec type: Union{ Nothing, String } # spec name: name

    function LoadBalancerBackendAddress2(;properties=nothing, name=nothing)
        o = new()
        validate_property(LoadBalancerBackendAddress2, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        validate_property(LoadBalancerBackendAddress2, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        o
    end
end # type LoadBalancerBackendAddress2

const _property_map_LoadBalancerBackendAddress2 = Dict{Symbol,Symbol}(Symbol("properties")=>Symbol("properties"), Symbol("name")=>Symbol("name"))
const _property_types_LoadBalancerBackendAddress2 = Dict{Symbol,String}(Symbol("properties")=>"LoadBalancerBackendAddressPropertiesFormat2", Symbol("name")=>"String")
Base.propertynames(::Type{ LoadBalancerBackendAddress2 }) = collect(keys(_property_map_LoadBalancerBackendAddress2))
Swagger.property_type(::Type{ LoadBalancerBackendAddress2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LoadBalancerBackendAddress2[name]))}
Swagger.field_name(::Type{ LoadBalancerBackendAddress2 }, property_name::Symbol) =  _property_map_LoadBalancerBackendAddress2[property_name]

function check_required(o::LoadBalancerBackendAddress2)
    true
end

function validate_property(::Type{ LoadBalancerBackendAddress2 }, name::Symbol, val)
end


text:
export add_timestamps_to_logger

# Add timestamps to logger
const logger_time_format = "HH:MM:SS.sss"
timestamp_logger(logger) = TransformerLogger(logger) do log
  merge(log, (; message = "$(format(now(), logger_time_format)) $(log.message)"))
end

logger_timestamps_on::Bool = false
function add_timestamps_to_logger()
    if !logger_timestamps_on
        logger = global_logger()
        logger |> timestamp_logger |> global_logger
        global logger_timestamps_on = true
    end
    return logger_timestamps_on
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
const turing_key_map = Dict(
    "hamiltonian_energy" => "energy",
    "hamiltonian_energy_error" => "energy_error",
    "is_adapt" => "tune",
    "max_hamiltonian_energy_error" => "max_energy_error",
    "nom_step_size" => "step_size_nom",
    "numerical_error" => "diverging",
)
const stan_key_map = Dict(
    "accept_stat__" => "acceptance_rate",
    "divergent__" => "diverging",
    "energy__" => "energy",
    "lp__" => "lp",
    "n_leapfrog__" => "n_steps",
    "stepsize__" => "step_size",
    "treedepth__" => "tree_depth",
)
const stats_key_map = merge(turing_key_map, stan_key_map)

"""
    reshape_values(x::AbstractArray) -> AbstractArray

Convert from `MCMCChains` variable values with dimensions `(ndraw, nchain, size...)` to
ArviZ's expected `(nchain, ndraw, size...)`.
"""
reshape_values(x::AbstractArray{T,N}) where {T,N} = permutedims(x, (2, 1, (3:N)...))

headtail(x) = x[1], x[2:end]

function split_locname(name)
    name = replace(name, r"[\[,]" => '.')
    name = replace(name, ']' => "")
    name, loc = headtail(split(name, '.'))
    isempty(loc) && return name, ()
    loc = tryparse.(Int, loc)
    Nothing <: eltype(loc) && return name, ()
    return name, tuple(loc...)
end

function varnames_locs_dict(loc_names, loc_str_to_old)
    vars_to_locs = Dict()
    for loc_name in loc_names
        var_name, loc = split_locname(loc_name)
        if var_name ∉ keys(vars_to_locs)
            vars_to_locs[var_name] = ([loc_str_to_old[loc_name]], [loc])
        else
            push!(vars_to_locs[var_name][1], loc_str_to_old[loc_name])
            push!(vars_to_locs[var_name][2], loc)
        end
    end
    # ensure that elements are ordered in the same order as they would be iterated
    for loc_name_locs in values(vars_to_locs)
        perm = sortperm(loc_name_locs[2]; by=CartesianIndex)
        permute!(loc_name_locs[1], perm)
        permute!(loc_name_locs[2], perm)
    end
    return vars_to_locs
end

function attributes_dict(chns::Chains)
    info = delete(chns.info, :hashedsummary)
    return Dict{String,Any}((string(k), v) for (k, v) in pairs(info))
end

function section_dict(chns::Chains, section)
    ndraws, _, nchains = size(chns)
    loc_names_old = getfield(chns.name_map, section) # old may be Symbol or String
    loc_names = string.(loc_names_old)
    loc_str_to_old = Dict(
        name_str => name_old for (name_str, name_old) in zip(loc_names, loc_names_old)
    )
    vars_to_locs = varnames_locs_dict(loc_names, loc_str_to_old)
    vars_to_arrays = Dict{String,Array}()
    for (var_name, names_locs) in vars_to_locs
        loc_names, locs = names_locs
        sizes = reduce((a, b) -> max.(a, b), locs)
        ndim = length(sizes)
        var_views = (@view(chns.value[:, n, :]) for n in loc_names)
        oldarr = reshape_values(replacemissing(convert(Array, cat(var_views...; dims=3))))
        if iszero(ndim)
            arr = dropdims(oldarr; dims=3)
        else
            arr = Array{Union{typeof(NaN),eltype(oldarr)}}(undef, nchains, ndraws, sizes...)
            fill!(arr, NaN)
            for i in eachindex(locs)
                arr[:, :, locs[i]...] = oldarr[:, :, i]
            end
        end
        vars_to_arrays[var_name] = arr
    end
    return vars_to_arrays
end

function chains_to_dict(
    chns::Chains; ignore=String[], section=:parameters, rekey_fun=identity
)
    section in sections(chns) || return Dict()
    chns_dict = section_dict(chns, section)
    removekeys!(chns_dict, ignore)
    return rekey_fun(chns_dict)
end

"""
    convert_to_inference_data(obj::Chains; group = :posterior, kwargs...) -> InferenceData

Convert the chains `obj` to an [`InferenceData`](@ref) with the specified `group`.

Remaining `kwargs` are forwarded to [`from_mcmcchains`](@ref).
"""
function convert_to_inference_data(chns::Chains; group=:posterior, kwargs...)
    group = Symbol(group)
    group === :posterior && return from_mcmcchains(chns; kwargs...)
    return from_mcmcchains(; group => chns)
end

@doc doc"""
    from_mcmcchains(posterior::MCMCChains.Chains; kwargs...) -> InferenceData
    from_mcmcchains(; kwargs...) -> InferenceData
    from_mcmcchains(
        posterior::MCMCChains.Chains,
        posterior_predictive::Any,
        predictions::Any,
        log_likelihood::Any;
        kwargs...
    ) -> InferenceData

Convert data in an `MCMCChains.Chains` format into an [`InferenceData`](@ref).

Any keyword argument below without an an explicitly annotated type above is allowed, so long
as it can be passed to [`convert_to_inference_data`](@ref).

# Arguments

- `posterior::MCMCChains.Chains`: Draws from the posterior

# Keywords

- `posterior_predictive::Any=nothing`: Draws from the posterior predictive distribution or
    name(s) of predictive variables in `posterior`
- `predictions::Any=nothing`: Out-of-sample predictions for the posterior.
- `prior::Any=nothing`: Draws from the prior
- `prior_predictive::Any=nothing`: Draws from the prior predictive distribution or name(s)
    of predictive variables in `prior`
- `observed_data::Dict{String,Array}=nothing`: Observed data on which the `posterior` is
    conditional. It should only contain data which is modeled as a random variable. Keys are
    parameter names and values.
- `constant_data::Dict{String,Array}=nothing`: Model constants, data included in the model
    which is not modeled as a random variable. Keys are parameter names and values.
- `predictions_constant_data::Dict{String,Array}=nothing`: Constants relevant to the model
     predictions (i.e. new `x` values in a linear regression).
- `log_likelihood::Any=nothing`: Pointwise log-likelihood for the data. It is recommended
     to use this argument as a dictionary whose keys are observed variable names and whose
     values are log likelihood arrays.
- `log_likelihood::String=nothing`: Name of variable in `posterior` with log likelihoods
- `library=MCMCChains`: Name of library that generated the chains
- `coords::Dict{String,Vector}=Dict()`: Map from named dimension to named indices
- `dims::Dict{String,Vector{String}}=Dict()`: Map from variable name to names of its
    dimensions
- `eltypes::Dict{String,DataType}=Dict()`: Apply eltypes to specific variables. This is used
    to assign discrete eltypes to discrete variables.

# Returns

- `InferenceData`: The data with groups corresponding to the provided data
"""
from_mcmcchains

function from_mcmcchains(
    posterior,
    posterior_predictive,
    predictions,
    log_likelihood;
    library=MCMCChains,
    eltypes=Dict(),
    kwargs...,
)
    kwargs = convert(Dict, merge((; dims=Dict()), kwargs))
    library = string(library)
    rekey_fun = d -> rekey(d, stats_key_map)

    # Convert chains to dicts
    if posterior === nothing
        post_dict = nothing
        stats_dict = nothing
    else
        post_dict = convert_to_eltypes(chains_to_dict(posterior), eltypes)
        stats_dict = chains_to_dict(posterior; section=:internals, rekey_fun=rekey_fun)
        stats_dict = enforce_stat_eltypes(stats_dict)
        stats_dict = convert_to_eltypes(stats_dict, Dict("is_accept" => Bool))
    end

    all_idata = InferenceData()
    for (group, group_data) in [
        :posterior_predictive => posterior_predictive,
        :predictions => predictions,
        :log_likelihood => log_likelihood,
    ]
        group_data === nothing && continue
        if group_data isa Union{Symbol,String}
            group_data = [string(group_data)]
        end
        if group_data isa Union{AbstractVector{Symbol},NTuple{N,Symbol} where {N}}
            group_data = map(string, group_data)
        end
        if group_data isa Union{AbstractVector{String},NTuple{N,String} where {N}}
            group_data = popsubdict!(post_dict, group_data)
        end
        group_dataset = if group_data isa Chains
            convert_to_dataset(group_data; library=library, eltypes=eltypes, kwargs...)
        else
            convert_to_dataset(group_data; library=library, kwargs...)
        end
        setattribute!(group_dataset, "inference_library", library)
        concat!(all_idata, InferenceData(; group => group_dataset))
    end

    attrs_library = Dict("inference_library" => library)
    if posterior === nothing
        attrs = attrs_library
    else
        attrs = merge(attributes_dict(posterior), attrs_library)
    end
    kwargs = convert(Dict, merge((; attrs=attrs, dims=Dict()), kwargs))
    post_idata = _from_dict(post_dict; sample_stats=stats_dict, kwargs...)
    concat!(all_idata, post_idata)
    return all_idata
end
function from_mcmcchains(
    posterior=nothing;
    posterior_predictive=nothing,
    predictions=nothing,
    prior=nothing,
    prior_predictive=nothing,
    observed_data=nothing,
    constant_data=nothing,
    predictions_constant_data=nothing,
    log_likelihood=nothing,
    library=MCMCChains,
    eltypes=Dict(),
    kwargs...,
)
    kwargs = convert(Dict, merge((; dims=Dict(), coords=Dict()), kwargs))

    all_idata = from_mcmcchains(
        posterior,
        posterior_predictive,
        predictions,
        log_likelihood;
        library=library,
        eltypes=eltypes,
        kwargs...,
    )

    if prior !== nothing
        pre_prior_idata = convert_to_inference_data(
            prior;
            posterior_predictive=prior_predictive,
            library=library,
            eltypes=eltypes,
            kwargs...,
        )
        prior_idata = rekey(
            pre_prior_idata,
            Dict(
                :posterior => :prior,
                :posterior_predictive => :prior_predictive,
                :sample_stats => :sample_stats_prior,
            ),
        )
        concat!(all_idata, prior_idata)
    end

    for (group, group_data) in [
        :observed_data => observed_data,
        :constant_data => constant_data,
        :predictions_constant_data => predictions_constant_data,
    ]
        group_data === nothing && continue
        group_data = convert_to_eltypes(group_data, eltypes)
        group_dataset = convert_to_constant_dataset(group_data; library=library, kwargs...)
        concat!(all_idata, InferenceData(; group => group_dataset))
    end

    return all_idata
end

"""
    from_cmdstan(posterior::Chains; kwargs...) -> InferenceData

Call [`from_mcmcchains`](@ref) on output of `CmdStan`.
"""
function from_cmdstan(posterior::Chains; kwargs...)
    return from_mcmcchains(posterior; library="CmdStan", kwargs...)
end


text: using ExpressCommands; ExpressCommands.comonicon_install()


text:
"""
A ProjMPO computes and stores the projection of an
MPO into a basis defined by an MPS, leaving a
certain number of site indices of the MPO unprojected.
Which sites are unprojected can be shifted by calling
the `position!` method.

Drawing of the network represented by a ProjMPO `P(H)`,
showing the case of `nsite(P)==2` and `position!(P,psi,4)`
for an MPS `psi`:

```
o--o--o-      -o--o--o--o--o--o <psi|
|  |  |  |  |  |  |  |  |  |  |
o--o--o--o--o--o--o--o--o--o--o H
|  |  |  |  |  |  |  |  |  |  |
o--o--o-      -o--o--o--o--o--o |psi>
```
"""
mutable struct ProjMPO <: AbstractProjMPO
  lpos::Int
  rpos::Int
  nsite::Int
  H::MPO
  LR::Vector{ITensor}
  ProjMPO(H::MPO) = new(0, length(H) + 1, 2, H, Vector{ITensor}(undef, length(H)))
end


text:
@enum ERROR_CODE begin
    ERROR_NONE=0x0
    ERROR_PROTOCOL=0x1
    ERROR_INTERNAL=0x2
    ERROR_FLOW_CONTROL=0x3
    ERROR_SETTINGS_TIMEOUT=0x4
    ERROR_STREAM_CLOSED=0x5
    ERROR_FRAME_SIZE=0x6
    ERROR_REFUSED_STREAM=0x7
    ERROR_CANCEL=0x8
    ERROR_COMPRESSION=0x9
    ERROR_CONNECT=0xa
    ERROR_ENHANCE_YOUR_CALM=0xb
    ERROR_INADEQUATE_SECURITY=0xc
    ERROR_HTTP_1_1_REQUIRED=0xd
end

struct ProtocolError
    message::String
end

struct InternalError
    message::String
end

struct NullError end

function goaway!(connection::HTTPConnection, error)
    error_code = if typeof(error) == NullError
        ERROR_NONE
    elseif typeof(error) == InternalError
        ERROR_INTERNAL
    else
        ERROR_PROTOCOL
    end

    frame = GoawayFrame(0x0, UInt32(error_code), Array{UInt8, 1}())

    put!(connection.channel_act_raw, frame)
end

Base.close(connection::HTTPConnection) = goaway!(connection, NullError())


text:
using ContactAlgorithms
using Test

@testset "ContactAlgorithms.jl" begin
    @test 1==1
end

max_stars_repo_name: SteffenPL/ContactAlgorithms

max_issues_repo_name: SteffenPL/ContactAlgorithms


text:
"""
    wcs_to_celestial_frame(wcs::WCSTransform)

Returns the reference frame of a WCSTransform.
The reference frame supported in Julia are FK5, ICRS and Galactic.
"""
function wcs_to_celestial_frame(wcs::WCSTransform)
    radesys = wcs.radesys

    xcoord = wcs.ctype[1][1:4]
    ycoord = wcs.ctype[2][1:4]

    if radesys == ""
        if xcoord == "GLON" && ycoord == "GLAT"
            radesys = "Gal"
        elseif xcoord == "TLON" && ycoord == "TLAT"
            radesys = "ITRS"
        end
    end

    return radesys
end

"""
    pad_edges(array_in::Matrix{T}) where {T}

Pads a given array and creates a border with edge elements.
"""
function pad_edges(array_in::Matrix{T}) where {T}
    image = Matrix{T}(undef, size(array_in)[1] + 2, size(array_in)[2] + 2)
    image[2:end-1,2:end-1] = array_in
    image[2:end-1,1] = array_in[:,1]
    image[2:end-1,end] = array_in[:,end]
    image[1,:] = image[2,:]
    image[end,:] = image[end-1,:]
    return image
end


text:
using Statistics
using Plots
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using JLD2
using Random
using Distributions
using FITSIO
using SparseArrays

push!(LOAD_PATH, pwd()*"/../../../main")
using DHC_2DUtils
push!(LOAD_PATH, pwd()*"/../../../scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using ReconFuncs

println("Log of Coefficients Check: Coeff_mask added in get_dbn line in invvar2")
numfile = Base.parse(Int, ENV["SLURM_ARRAY_TASK_ID"])
println(numfile, ARGS[1], ARGS[2])

if ARGS[1]=="log"
    logbool=true
else
    if ARGS[1]!="reg" error("Invalid log arg") end
    logbool=false
end

if ARGS[2]=="apd"
    apdbool = true
else
    if ARGS[2]!="nonapd" error("Invalid apd arg") end
    apdbool=false
end

if ARGS[3]=="iso"
    isobool = true
else
    if ARGS[3]!="noiso" error("Invalid iso arg") end
    isobool=false
end

#read in data
direc = ARGS[4]
datfile = direc * "Data_" * string(numfile) * ".jld2" #Replace w SLURM array
loaddf = load(datfile)
true_img = loaddf["true_img"]
init = loaddf["init"]

fname_save = direc * "SFDTargSFDCov/" * ARGS[1] * "_" * ARGS[2] * "_" * ARGS[3] * "/LogCoeff/" * string(numfile) * ARGS[5] #Change
Nx=size(true_img)[1]
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(S1iso, Nf) = size(filter_hash["S1_iso_mat"])
(S2iso, Nfsq) = size(filter_hash["S2_iso_mat"])
dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool) #Iso #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense
if dhc_args[:iso]
    coeff_mask = falses(2+S1iso+S2iso)
    coeff_mask[S1iso+3:end] .= true
else #Not iso
    coeff_mask = falses(2+Nf+Nf^2)
    coeff_mask[Nf+3:end] .= triu(trues(Nf, Nf))[:]
end

recon_settings = Dict([("log", logbool), ("covar_type", "sfd_dbn"), ("target_type", "sfd_dbn"), ("Invcov_matrix", ARGS[6]), ("epsvalue", 1e-6)])
println("Apdbool", apdbool)
println("Isobool", isobool)
println(recon_settings)

sfdall = readsfd(Nx, logbool=logbool)
dbnocffs = log.(get_dbn_coeffs(sfdall, filter_hash, dhc_args, coeff_mask = coeff_mask))
fmean = mean(dbnocffs, dims=1)
fcov = (dbnocffs .- fmean)' * (dbnocffs .- fmean) ./(size(dbnocffs)[1] -1)
sig2 = diag(fcov)
if recon_settings["Invcov_matrix"]=="Diagonal"
    fcovinv = invert_covmat(sig2)
elseif recon_settings["Invcov_matrix"]=="Diagonal+Eps"
    fcovinv = invert_covmat(sig2, recon_settings["epsvalue"]) #DEBUG
elseif recon_settings["Invcov_matrix"]=="Full+Eps"
    fcovinv = invert_covmat(fcov, recon_settings["epsvalue"])
else#s2icov
    fcovinv = invert_covmat(fcov)
end

invvar1 = diag(fcovinv)

#Check wrt calculation here
dbn = log.(get_dbn_coeffs(sfdall, filter_hash, dhc_args, coeff_mask = coeff_mask))
invvar2 = std(dbn, dims=1)[:]
println(size(invvar2))
invvar2 = (invvar2.^(-2))
println("InvVar1 ", invvar1[1:10])
println("InvVar2 ", invvar2[1:10])
println("Mean Abs Error = ", mean(abs.(invvar1 .- invvar2)))

var1 = diag(fcov)
var2 = std(dbn, dims=1)[:]
var2 = var2.^2
println("Var1 ", var1[1:10])
println("Var2 ", var2[1:10])
println("Mean Abs Error = ", mean(abs.(var1 .- var2)))
println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
@testset "VF Inteprolation" begin
    let prim = dcdp_primitives(flowfuncname, β, wp, zspace, Πp, ψspace),
        tmpv = dcdp_tmpvars(prim),
        royalty_rates = [0.25,],
        geology_types = 2:2

        sev = SharedEV([1,], prim, geology_types, royalty_rates)
        EVcopy = similar(sev.EV)
        EVcopy2 = similar(sev.EV)

        serial_solve_vf_all!(sev, tmpv, prim, θfull, false)
        EVcopy .= sev.EV
        EVcopy2 .= sev.EV
        @test !(EVcopy === sev.EV)

        EVcopy .= sev.EV
        @test all(EVcopy .== sev.EV)
        spec = (BSpline(Quadratic(InPlace())), BSpline(Quadratic(InPlace())), BSpline(Quadratic(InPlace())), NoInterp(), NoInterp(), NoInterp(),)
        sitp_test = interpolate!(EVcopy, spec, OnCell() )
        @test sitp_test.coefs === EVcopy
        @test !all(sitp_test.coefs .== sev.EV)

        sitp = ItpSharedEV(sev, prim)
        sitp.itypes

        serial_prefilterByView!(sev, sitp)

        @test all(sitp_test.coefs .== sev.EV)
        @test !all(sev.EV .== EVcopy2)
    end
end


text:
# Define kernel interface

"""
Abstract type representing the operations to be used in the [`apply`](@ref) function.

Derived types must implement the following method:

- [`apply_kernel!(cache,k,x...)`](@ref)

and optionally these ones:

- [`kernel_cache(k,x...)`](@ref)
- [`kernel_return_type(k,x...)`](@ref)

The kernel interface can be tested with the [`test_kernel`](@ref) function.

Note that most of the functionality implemented in terms of this interface
relies in duck typing. That is, it is not strictly needed to work with types
that inherit from `Kernel`. This is specially useful in order to accommodate
existing types into this framework without the need to implement a wrapper type
that inherits from `Kernel`. For instance, a default implementation is available
for `Function` objects.  However, we recommend that new types inherit from `Kernel`.

"""
abstract type Kernel <: GridapType end

"""
    kernel_return_type(f,x...)

Returns the type of the result of calling kernel `f` with
arguments of the types of the objects `x`.

It defaults to `typeof(kernel_testitem(f,x...))`
"""
function kernel_return_type(f,x...)
  typeof(kernel_testitem(f,x...))
end

"""
    kernel_cache(f,x...)

Returns the `cache` needed to apply kernel `f` with arguments
of the same type as the objects in `x`.
This function returns `nothing` by default.
"""
kernel_cache(f,x...) = nothing

"""
    apply_kernel!(cache,f,x...)

Applies the kernel `f` at the arguments `x...` using
the scratch data provided in the given `cache` object. The `cache` object
is built with the [`kernel_cache`](@ref) function using arguments of the same type as in `x`.
In general, the returned value `y` can share some part of its state with the `cache` object.
If the result of two or more invocations of this function need to be accessed simultaneously
(e.g., in multi-threading), create and use various `cache` objects (e.g., one cache
per thread).
"""
function apply_kernel!(cache,f,x...)
  @abstractmethod
end

# Testing the interface

"""
    test_kernel(f,x::Tuple,y,cmp=(==))

Function used to test if the kernel `f` has been
implemented correctly. `f` is a kernel object, `x` is a tuple containing the arguments
of the kernel, and `y` is the expected result. Function `cmp` is used to compare
the computed result with the expected one. The checks are done with the `@test`
macro.
"""
function test_kernel(f,x::Tuple,y,cmp=(==))
  z = apply_kernel(f,x...)
  @test cmp(z,y)
  @test typeof(z) == kernel_return_type(f,x...)
  cache = kernel_cache(f,x...)
  z = apply_kernel!(cache,f,x...)
  @test cmp(z,y)
  z = apply_kernel!(cache,f,x...)
  @test cmp(z,y)
  z = kernel_testitem!(cache,f,x...)
  #@test cmp(typeof(z),typeof(y))
end

# Functions on kernel objects

"""
    apply_kernel(f,x...)

apply the kernel `f` at the arguments in `x` by creating a temporary cache
internally. This functions is equivalent to
```jl
cache = kernel_cache(f,x...)
apply_kernel!(cache,f,x...)
```
"""
function apply_kernel(f,x...)
  cache = kernel_cache(f,x...)
  y = apply_kernel!(cache,f,x...)
  y
end

# Work with several kernels at once

"""
    kernel_caches(fs::Tuple,x...) -> Tuple

Returns a tuple with the cache corresponding to each kernel in `fs`
for the arguments `x...`.
"""
function kernel_caches(fs::Tuple,x...)
  _kernel_caches(x,fs...)
end

function _kernel_caches(x::Tuple,a,b...)
  ca = kernel_cache(a,x...)
  cb = kernel_caches(b,x...)
  (ca,cb...)
end

function _kernel_caches(x::Tuple,a)
  ca = kernel_cache(a,x...)
  (ca,)
end

"""
    apply_kernels!(caches::Tuple,fs::Tuple,x...) -> Tuple

Applies the kernels in the tuple `fs` at the arguments `x...`
by using the corresponding cache objects in the tuple `caches`.
The result is also a tuple containing the result for each kernel in `fs`.
"""
@inline function apply_kernels!(cfs::Tuple,f::Tuple,x...)
  _apply_kernels!(cfs,x,f...)
end

@inline function _apply_kernels!(cfs,x,f1,f...)
  cf1, cf = _split(cfs...)
  f1x = apply_kernel!(cf1,f1,x...)
  fx = apply_kernels!(cf,f,x...)
  (f1x,fx...)
end

@inline function _apply_kernels!(cfs,x,f1)
  cf1, = cfs
  f1x = apply_kernel!(cf1,f1,x...)
  (f1x,)
end

@inline function _split(a,b...)
  (a,b)
end

"""
    kernel_return_types(f::Tuple,x...) -> Tuple

Computes the return types of the kernels in `f` when called
with arguments `x`.
"""
function kernel_return_types(f::Tuple,x...)
  _kernel_return_types(x,f...)
end

function _kernel_return_types(x::Tuple,a,b...)
  Ta = kernel_return_type(a,x...)
  Tb = kernel_return_types(b,x...)
  (Ta,Tb...)
end

function _kernel_return_types(x::Tuple,a)
  Ta = kernel_return_type(a,x...)
  (Ta,)
end

function kernel_testitem(k,x...)
  cache = kernel_cache(k,x...)
  kernel_testitem!(cache,k,x...)
end

@inline function kernel_testitem!(cache,k,x...)
  apply_kernel!(cache,k,x...)
end

# Include some well-known types in this interface

function kernel_return_type(f::Function,x...)
  Ts = map(typeof,x)
  return_type(f,Ts...)
end

@inline apply_kernel!(::Nothing,f::Function,args...) = f(args...)

# Some particular cases

const NumberOrArray = Union{Number,AbstractArray}

"""
    bcast(f::Function)

Returns a kernel object that represents the "broadcasted" version of the given
function `f`.
"""
bcast(f::Function) = BCasted(f)

struct BCasted{F<:Function} <: Kernel
  f::F
end

function kernel_return_type(f::BCasted,x::Number...)
  Ts = map(typeof,x)
  return_type(f.f,Ts...)
end

function kernel_cache(f::BCasted,x::Number...)
  nothing
end

@inline function apply_kernel!(::Nothing,f::BCasted,x::Number...)
  f.f(x...)
end

function kernel_return_type(f::BCasted,x::NumberOrArray...)
  typeof(kernel_cache(f,x...).array)
end

function kernel_cache(f::BCasted,x::NumberOrArray...)
  s = _sizes(x...)
  bs = Base.Broadcast.broadcast_shape(s...)
  Te = map(numbertype,x)
  T = return_type(f.f,Te...)
  N = length(bs)
  r = Array{T,N}(undef,bs)
  ri = testvalue(T)
  fill!(r,ri)
  cache = CachedArray(r)
   _prepare_cache(cache,x...)
end

numbertype(a::AbstractArray) = eltype(a)

numbertype(a::Number) = typeof(a)

@inline function apply_kernel!(cache,f::BCasted,x::NumberOrArray...)
  r = _prepare_cache(cache,x...)
  a = r.array
  broadcast!(f.f,a,x...)
  a
end

@inline function _prepare_cache(c,x...)
  s = _sizes(x...)
  bs = Base.Broadcast.broadcast_shape(s...)
  if bs != size(c)
    setsize!(c,bs)
  end
  c
end

# TODO use map
@inline function _sizes(a,x...)
  (_size(a), _sizes(x...)...)
end

@inline function _sizes(a)
  (_size(a),)
end

@inline _size(a) = size(a)
@inline _size(a::Number) = (1,)

"""
    elem(f::Function)

Returns a kernel that represents the element-wise
version of the operation `f`
It does not broadcast in singleton axes. Thus, allows some
performance optimizations with respect to broadcast.
"""
elem(f::Function) = Elem(f)

struct Elem{F} <: Kernel
  f::F
  Elem(f::Function) = new{typeof(f)}(f)
end

# TODO more tests here

# It defaults to bcast (TODO test these ones)

@inline function apply_kernel!(cache,k::Elem,x::NumberOrArray...)
  b = bcast(k.f)
  apply_kernel!(cache,b,x...)
end

function kernel_cache(k::Elem,x::NumberOrArray...)
  b = bcast(k.f)
  kernel_cache(b,x...)
end

function kernel_return_type(k::Elem,x::NumberOrArray...)
  b = bcast(k.f)
  kernel_return_type(b,x...)
end

# More Efficient implementations

# Number

function kernel_return_type(k::Elem,a::Number)
  return_type(k.f,typeof(a))
end

function kernel_cache(k::Elem,a::Number)
  nothing
end

@inline function apply_kernel!(c,k::Elem,a::Number)
  k.f(a)
end

# Array

function kernel_return_type(k::Elem,a::AbstractArray)
  typeof(kernel_cache(k,a).array)
end

function kernel_cache(k::Elem,a::AbstractArray)
  T = return_type(k.f,eltype(a))
  CachedArray(similar(a,T))
end

@inline function apply_kernel!(c,f::Elem,a::AbstractArray)
  setsize!(c,size(a))
  r = c.array
  for i in eachindex(a)
    @inbounds r[i] = f.f(a[i])
  end
  r
end

# Array vs Array

function kernel_return_type(k::Elem,a::AbstractArray,b::AbstractArray)
  typeof(kernel_cache(k,a,b).array)
end

function kernel_cache(k::Elem,a::AbstractArray,b::AbstractArray)
  _checks(a,b)
  T = return_type(k.f,eltype(a),eltype(b))
  CachedArray(similar(a,T))
end

@inline function apply_kernel!(c,f::Elem,a::AbstractArray,b::AbstractArray)
  _checks(a,b)
  setsize!(c,size(a))
  r = c.array
  for i in eachindex(a)
    @inbounds r[i] = f.f(a[i],b[i])
  end
  r
end

# Number vs Number

function kernel_return_type(k::Elem,a::Number,b::Number)
  return_type(k.f,typeof(a),typeof(b))
end

function kernel_cache(k::Elem,a::Number,b::Number)
  nothing
end

@inline function apply_kernel!(c,k::Elem,a::Number,b::Number)
  k.f(a,b)
end

# Array vs Number

function kernel_return_type(k::Elem,a::AbstractArray,b::Number)
  typeof(kernel_cache(k,a,b).array)
end

function kernel_cache(k::Elem,a::AbstractArray,b::Number)
  T = return_type(k.f,eltype(a),typeof(b))
  CachedArray(similar(a,T))
end

@inline function apply_kernel!(c,k::Elem,a::AbstractArray,b::Number)
  setsize!(c,size(a))
  r = c.array
  for i in eachindex(a)
    @inbounds r[i] = k.f(a[i],b)
  end
  r
end

# Number vs Array

function kernel_return_type(k::Elem,a::Number,b::AbstractArray)
  typeof(kernel_cache(k,a,b).array)
end

function kernel_cache(k::Elem,a::Number,b::AbstractArray)
  T = return_type(k.f,typeof(a),eltype(b))
  CachedArray(similar(b,T))
end

@inline function apply_kernel!(c,k::Elem,a::Number,b::AbstractArray)
  setsize!(c,size(b))
  r = c.array
  for i in eachindex(b)
    @inbounds r[i] = k.f(a,b[i])
  end
  r
end

function _checks(a,b)
  @assert size(a) == size(b) "Sizes must agree."
  nothing
end

"""
    contract(f::Function)

Like the dot product between to vectors, but using operation `f` instead
of `*` between components.

!!! warning
    not needed any more, to be deleted

# Examples

```jldoctests
using Gridap.Arrays
k = contract(-)
apply_kernel(k,[1,2],[2,4]) # Equivalent to (1-2) + (2-4)
# output
-3
```
"""
contract(f::Function) = Contracted(f)

struct Contracted{F} <: Kernel
  f::F
  Contracted(f::Function) = new{typeof(f)}(f)
end

function kernel_return_type(k::Contracted,a::AbstractArray,b::AbstractArray)
  return_type(k.f,eltype(a),eltype(b))
end

function kernel_cache(k::Contracted,a::AbstractArray,b::AbstractArray)
  kernel_return_type(k,a,b)
end

@inline function apply_kernel!(T,f::Contracted,a::AbstractArray,b::AbstractArray)
  _checks(a,b)
  c = zero(T)
  for i in eachindex(a)
    c += f.f(a[i],b[i])
  end
  c
end

struct MulKernel <: Kernel end

function kernel_cache(k::MulKernel,a,b)
  c = a*b
  CachedArray(c)
end

@inline function apply_kernel!(cache,k::MulKernel,a::AbstractMatrix,b::AbstractVector)
  m = axes(a,1)
  setaxes!(cache,(m,))
  c = cache.array
  mul!(c,a,b)
  c
end

@inline function apply_kernel!(cache,k::MulKernel,a::AbstractMatrix,b::AbstractMatrix)
  m = axes(a,1)
  n = axes(b,2)
  setaxes!(cache,(m,n))
  c = cache.array
  mul!(c,a,b)
  c
end

struct MulAddKernel{T} <: Kernel
  α::T
  β::T
end

function kernel_cache(k::MulAddKernel,a,b,c)
  d = a*b+c
  CachedArray(d)
end

@inline function apply_kernel!(cache,k::MulAddKernel,a,b,c)
  setaxes!(cache,axes(c))
  d = cache.array
  copyto!(d,c)
  mul!(d,a,b,k.α,k.β)
  d
end


text:
module TestUtilities

# using Revise
using Test
using MLJ
using ScientificTypes
using Tables
using OrderedCollections

d = Dict(
    :x => Set([1, 2]),
    :y => Set([2, 3, 5]),
    :z => Set([4, 7]),
    :a => Set([8, 1]),
    :b => Set([4,]),
    :w => Set([3, 1, 2]),
    :t => Set([0,]))

dinv = Dict(
    0 => Set([:t,]),
    1 => Set([:x, :a, :w]),
    2 => Set([:x, :y, :w]),
    3 => Set([:y, :w]),
    4 => Set([:z, :b]),
    5 => Set([:y,]),
    7 => Set([:z,]),
    8 => Set([:a,]))

@test MLJ.inverse(d) == dinv

d = LittleDict(
    :x => Set([1, 2]),
    :y => Set([2, 3, 5]),
    :z => Set([4, 7]),
    :a => Set([8, 1]),
    :b => Set([4,]),
    :w => Set([3, 1, 2]),
    :t => Set([0,]))

dinv = LittleDict(
    0 => Set([:t,]),
    1 => Set([:x, :a, :w]),
    2 => Set([:x, :y, :w]),
    3 => Set([:y, :w]),
    4 => Set([:z, :b]),
    5 => Set([:y,]),
    7 => Set([:z,]),
    8 => Set([:a,]))

@test MLJ.inverse(d) == dinv

d = LittleDict()
d[:test] = Tuple{Union{Continuous,Missing}, Finite}
d["junk"] = LittleDict{Any,Any}("H" => Missing, :c => Table(Finite,Continuous),
                                :cross => "lemon", :t => :w, "r" => "r")
d["a"] = "b"
d[:f] = true
d["j"] = :post

@test MLJ.decode_dic(MLJ.encode_dic(d)) == d

mutable struct M
    a1
    a2
end
mutable struct A1
    a11
    a12
end
mutable struct A2
    a21
end
mutable struct A21
    a211
    a212
end

@testset "recursive getproperty, setproperty!" begin

    m = (a1 = (a11 = 10, a12 = 20), a2 = (a21 = (a211 = 30, a212 = 40),))

    @test MLJ.recursive_getproperty(m, :(a1.a12)) == 20
    @test MLJ.recursive_getproperty(m, :a1) == (a11 = 10, a12 = 20)
    @test MLJ.recursive_getproperty(m, :(a2.a21.a212)) == 40

    m = M(A1(10, 20), A2(A21(30, 40)))
    MLJ.recursive_setproperty!(m, :(a2.a21.a212), 42)
    @test MLJ.recursive_getproperty(m, :(a1.a11)) == 10
    @test MLJ.recursive_getproperty(m, :(a1.a12)) == 20
    @test MLJ.recursive_getproperty(m, :(a2.a21.a211)) == 30
    @test MLJ.recursive_getproperty(m, :(a2.a21.a212)) == 42
    @test MLJ.recursive_getproperty(
        MLJ.recursive_getproperty(m, :(a2.a21)), :a212) == 42

end

end # module
true


text:
for dim in (2,3)
    @eval begin
        function (::Type{QuadratureRule{$dim,RefTetrahedron}})(quad_type::GrundmannMoeller, s::Int)
            points, weigths = _compute_p_w(s, $dim)
            return QuadratureRule{$dim,RefTetrahedron,Float64}(weigths, points)
        end
    end
end

function _compute_p_w(s::Int, n::Int, T = Float64)
    d = 2 * s + 1
    k = 0
    nnodes = binomial(n + s + 1, s)
    weights = Vector{T}(undef, nnodes)
    kk = nnodes
    points = Vector{Vec{n,T}}(undef, nnodes)
    exponentials = _get_all_exponentials(s, n+1)
    for i = 0:s
        w::T = ((-1)^i * 2.0^(-2*s) * (d+n-2*i)^d)/(factorial(i) * factorial(d+n-i))
        for p in _get_all_exponentials(s-i, n+1)
            k += 1
            points[k] = Vec{n}(j -> (2*p[j+1] + 1)/(d+n-2*i))
            weights[k] = w
        end
    end
    return points, weights./sum(2*weights)
end

function _get_all_exponentials(n::Int,k::Int)
    a = Vector{Int}(undef, k)
    exponentials = Vector{typeof(a)}(undef,binomial(n+2,k-1))
    t::Int = n
    h = 0
    a[1] = n
    a[2:k] .= 0
    exponentials[1] = copy(a)
    idx = 2
    while a[k] != n
        if ( 1 < t ); h = 0; end
        h = h + 1
        t = a[h]
        a[h] = 0
        a[1] = t - 1
        a[h+1] = a[h+1] + 1
        exponentials[idx] = copy(a)
        idx +=1
    end
    return exponentials
end