text:
using Graphs, SimpleWeightedGraphs, StatsBase, SparseArrays, LinearAlgebra, Distributions

"""
Given a DeBruijn graph ko, its corresponding
first order graph fo, and an integer order k,
construct the Xi Matrix for ko and compute
p-values on every edge.
"""
function hypa(k::Integer,
        ko::SimpleWeightedDiGraph, ko_map::Dict{Tuple, Integer},
        fo::SimpleWeightedDiGraph, fo_map::Dict{String, Integer})
    fo_adj = sparse(adjacency_matrix(fo))
    ko_adj = sparse(adjacency_matrix(ko))
    # Construct Xi
    Xi = compute_xi(ko_adj, fo_adj, ko_map, fo_map)
    # Compute pvals
    pvals = compute_pvals(ko_adj, Xi, false)
    return pvals, Xi
end

"""
Compute p-values
"""
function compute_pvals(adjacency, Xi, log_p::Bool=false)
    total_observations = sum(Int64, adjacency)
    total_xi = sum(Int64, Xi)
    nonzero_xi = findall(>(0), Xi)
    pvals = spzeros(Xi.m, Xi.n)
    for cidx in nonzero_xi
        observed = adjacency[cidx]
        xi = Xi[cidx]
        hy = Hypergeometric(total_observations, total_xi-total_observations, xi)
        if log_p
            pvals[cidx] = logcdf(hy, observed)
        else
            pvals[cidx] = cdf(hy, observed)
        end
    end
    return pvals
end

"""
Draw a sample from the ensemble defined
by adjacency and Xi.
"""
function draw_sample(adjacency, Xi)
    sampled_adjacency = spzeros(adjacency.m, adjacency.n)
    xi_sum = sum(Xi)
    adj_sum = sum(adjacency)
    xi_accum = 0
    for cidx in findall(>(0), Xi)
        xi_accum += Xi[cidx]
        hy = Hypergeometric(Xi[cidx], xi_sum-xi_accum, adj_sum)
        sampled_adjacency[cidx] = rand(hy)
        adj_sum -= sampled_adjacency[cidx]
    end
    return sampled_adjacency
end

"""
FUNCTIONS FOR COMPUTING THE XI MATRIX
"""

"""
Compute Xi matrix.
"""
function compute_xi( ko_adj, fo_adj,
        ko_map::Dict{Tuple, Integer}, fo_map::Dict{String, Integer},
        tol::Float64=0.0001)
    fo_rev_map = Dict(val=>key for (key,val) in fo_map)
    ko_rev_map = Dict(val=>key for (key,val) in ko_map)
    Xi = spzeros(Float64, nv(ko), nv(ko))
    for ko_node in range(1, nv(ko))
        # Get last first-order node in node
        ko_node_mapped = ko_rev_map[ko_node]
        last_node = ko_node_mapped[end]
        fo_node = fo_map[last_node]
        # Get all neighbors of last node
        neighbors = findall(>(0), fo_adj[fo_node,:])
        # for each neighbor, check if corresponding
        # higher-order node is in ko
        prefix_mapped = [u for u in ko_node_mapped[2:end]]
        out_deg = sum(ko_adj[ko_node, findall(>(0), ko_adj[ko_node,:])])
        if out_deg < 1
            continue
        end
        for ne in neighbors
            candidate_node = Vector(prefix_mapped)
            push!(candidate_node, fo_rev_map[ne])
            candidate_node = Tuple([u for u in candidate_node])
            if haskey(ko_map, candidate_node)
                # if so, compute a xival. Else do nothing.
                candidate_idx = ko_map[candidate_node]
                in_deg = sum(ko_adj[findall(>(0), ko_adj[:,candidate_idx]), candidate_idx])
                Xi[ko_node, candidate_idx] = out_deg * in_deg
            end
        end
    end
    fit_xi!(Xi, ko, tol)
    return trunc.(Int64, Xi)
end

"""
Fit the Xi matrix by redistributing
the excess degree from impossible
edges.
"""
function fit_xi!(Xi, ko, tol::Float64)
    adj = adjacency_matrix(ko)
    m = sum(adj)
    indegs = sum(adj, dims=1)
    outdegs = sum(adj, dims=2)
    xi_sum = sum(Xi)
    rmse = compute_rmse(indegs, outdegs, Xi, xi_sum, m)
    println("rmse: $rmse")
    while rmse > tol
        # Fit rows
        Xi = xifit_row(m, Xi, outdegs)
        # Fit columns
        Xi = transpose(xifit_row(m, transpose(Xi), indegs))

        xi_sum = sum(Xi)
        rmse_new = compute_rmse(indegs, outdegs, Xi, xi_sum, m)
        println("rmse: $rmse_new")

        if rmse <= rmse_new || rmse < tol
            break
        end
        rmse = rmse_new
    end
    return Xi
end

"""
Operations to redistribute for a single row
or column.
"""
function xifit_row(m, Xi, degs)
    xi_sum = sum(Xi)
    Xi .= (Xi .* m^2) ./ xi_sum
    xi_sum = sum(Xi)

    exp_degs = vec(sum((Xi.*m)./xi_sum, dims=2))
    nonzero = findall(!=(0), exp_degs)
    ratio = spzeros(length(exp_degs))
    ratio[nonzero] .= degs[nonzero]./exp_degs[nonzero]
    Xi = ratio.*Xi
    return Xi
end

"""
Compute RMSE
"""
function compute_rmse(indegs, outdegs, Xi, xi_sum, m)
    in_sum = sum((Xi.*m) ./ xi_sum, dims=1)
    in_rmse = (in_sum-indegs).^2
    in_rmse = sum(sqrt.(in_rmse)./2)

    out_sum = sum((Xi.*m) ./ xi_sum, dims=2)
    out_rmse = (out_sum-outdegs).^2
    out_rmse = sum(sqrt.(out_rmse)./2)

    return out_rmse + in_rmse
end


text:
using Test
using Interpolations

@testset "Monotonic" begin
    x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
    y = [(true, [-3.0, 0.0, 5.0, 10.0, 18.0, 22.0]),
         (false, [10.0, 0.0, -5.0, 10.0, -8.0, -2.0])]

    # second item indicates if interpolating function can overshoot
    # for non-monotonic data
    itypes = [(LinearMonotonicInterpolation(), false),
        (FiniteDifferenceMonotonicInterpolation(), true),
        (CardinalMonotonicInterpolation(0.0), true),
        (CardinalMonotonicInterpolation(0.5), true),
        (CardinalMonotonicInterpolation(1.0), false),
        (FritschCarlsonMonotonicInterpolation(), true),
        (FritschButlandMonotonicInterpolation(), false),
        (SteffenMonotonicInterpolation(), false)]

    for (it, overshoot) in itypes
        for yi in 1:length(y)
            monotonic, ys = y[yi]
            itp = interpolate(x, ys, it)
            for j in 1:6
                # checking values at nodes
                @test itp(x[j]) ≈ ys[j] rtol = 1.e-15 atol = 1.e-14

                # checking overshoot for non-monotonic data
                # and monotonicity for monotonic data
                if !monotonic && overshoot
                    continue
                end
                if j < 6
                    r = range(x[j], stop = x[j+1], length = 100)
                    for k in 1:length(r)-1
                        @test (ys[j+1] - ys[j]) * (itp(r[k+1]) - itp(r[k])) >= 0.0
                    end
                end
            end
        end
    end

    # fail tests
    xWrong = [0.0, 1.0, -5.0, 3.0]
    y2 = [1.0, 2.0, 3.0, 4.0]
    y3 = [-3.0 0.0 5.0 10.0 18.0 22.0]

    for (it, overshoot) in itypes
        @test_throws ErrorException interpolate(xWrong, y2, it)
        @test_throws DimensionMismatch interpolate(x, y2, it)
        @test_throws MethodError interpolate(x, y3, it)
    end
end

ext: jl


text:
# # Geometry optimization
#
# We use the DFTK and Optim packages in this example to find the minimal-energy
# bond length of the ``H_2`` molecule. We setup ``H_2`` in an
# LDA model just like in the [Tutorial](@ref) for silicon.
using DFTK
using Optim
using LinearAlgebra
using Printf

kgrid = [1, 1, 1]       # k-point grid
Ecut = 5                # kinetic energy cutoff in Hartree
tol = 1e-8              # tolerance for the optimization routine
a = 10                  # lattice constant in Bohr
lattice = a * I(3)
H = ElementPsp(:H, psp=load_psp("hgh/lda/h-q1"));
atoms = [H, H]

# We define a blochwave and a density to be used as global variables so that we
# can transfer the solution from one iteration to another and therefore reduce
# the optimization time.

ψ = nothing
ρ = nothing

# First, we create a function that computes the solution associated to the
# position ``x \in \mathbb{R}^6`` of the atoms in reduced coordinates
# (cf. [Reduced and cartesian coordinates](@ref) for more
# details on the coordinates system).
# They are stored as a vector: `x[1:3]` represents the position of the
# first atom and `x[4:6]` the position of the second.
# We also update `ψ` and `ρ` for the next iteration.

function compute_scfres(x)
    model = model_LDA(lattice, atoms, [x[1:3], x[4:6]])
    basis = PlaneWaveBasis(model; Ecut, kgrid)
    global ψ, ρ
    if isnothing(ρ)
        ρ = guess_density(basis)
    end
    scfres = self_consistent_field(basis; ψ=ψ, ρ=ρ,
                                   tol=tol / 10, callback=info->nothing)
    ψ = scfres.ψ
    ρ = scfres.ρ
    scfres
end;

# Then, we create the function we want to optimize: `fg!` is used to update the
# value of the objective function `F`, namely the energy, and its gradient `G`,
# here computed with the forces (which are, by definition, the negative gradient
# of the energy).

function fg!(F, G, x)
    scfres = compute_scfres(x)
    if G != nothing
        grad = compute_forces(scfres)
        G .= -[grad[1]; grad[2]]
    end
    scfres.energies.total
end;

# Now, we can optimize on the 6 parameters `x = [x1, y1, z1, x2, y2, z2]` in
# reduced coordinates, using `LBFGS()`, the default minimization algorithm
# in Optim. We start from `x0`, which is a first guess for the coordinates. By
# default, `optimize` traces the output of the optimization algorithm during the
# iterations. Once we have the minimizer `xmin`, we compute the bond length in
# cartesian coordinates.

x0 = vcat(lattice \ [0., 0., 0.], lattice \ [1.4, 0., 0.])
xres = optimize(Optim.only_fg!(fg!), x0, LBFGS(),
                Optim.Options(show_trace=true, f_tol=tol))
xmin = Optim.minimizer(xres)
dmin = norm(lattice*xmin[1:3] - lattice*xmin[4:6])
@printf "\nOptimal bond length for Ecut=%.2f: %.3f Bohr\n" Ecut dmin

# We used here very rough parameters to generate the example and
# setting `Ecut` to 10 Ha yields a bond length of 1.523 Bohr,
# which [agrees with ABINIT](https://docs.abinit.org/tutorial/base1/).
#
# !!! note "Degrees of freedom"
#     We used here a very general setting where we optimized on the 6 variables
#     representing the position of the 2 atoms and it can be easily extended
#     to molecules with more atoms (such as ``H_2O``). In the particular case
#     of ``H_2``, we could use only the internal degree of freedom which, in
#     this case, is just the bond length.


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("MLJTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# This tutorial follows loosely [the docs](https://juliadata.github.io/CategoricalArrays.jl/latest/using.html).## ## Defining a categorical vector
using CategoricalArrays

v = categorical(["AA", "BB", "CC", "AA", "BB", "CC"])

# This declares a categorical vector, i.e. a Vector whose entries are expected to represent a group or category.# You can retrieve the group labels using `levels`:
levels(v)

# which, by  default, returns the labels in lexicographic order.## ## Working with categoricals## ### Ordered categoricals## You can specify that categories are *ordered* by specifying `ordered=true`, the order then follows that of the levels. If you wish to change that order, you  need to  use the `levels!` function.# Let's see two examples.
v = categorical([1, 2, 3, 1, 2, 3, 1, 2, 3], ordered=true)

levels(v)

# Here the lexicographic order matches what we want so no  need to change it, since we've specified  that the categories are ordered we can do:
v[1] < v[2]

# Let's now consider another example
v = categorical(["high", "med", "low", "high", "med", "low"], ordered=true)

levels(v)

# The levels follow the lexicographic order which  is not what  we want:
v[1] < v[2]

# In order to re-specify the order we need to  use `levels!`:
levels!(v, ["low", "med", "high"])

# now things are properly ordered:
v[1] < v[2]

# ### Missing values## You can also have a categorical vector with missing values:
v = categorical(["AA", "BB", missing, "AA", "BB", "CC"]);

# that doesn't change the levels:
levels(v)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
# Functions for querying network properties.

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

"""
    species(network)

Given a [`ReactionSystem`](@ref), return a vector of species `Variable`s.

Notes:
- If `ModelingToolkit.get_systems(network)` is not empty, may allocate. Otherwise returns
  `ModelingToolkit.get_states(network)`.
"""
function species(network)
    isempty(get_systems(network)) ? get_states(network) : states(network)
end

"""
    params(network)

Given a [`ReactionSystem`](@ref), return a vector of parameter `Variable`s.

Notes:
- If `ModelingToolkit.get_systems(network)` is not empty, may allocate. Otherwise returns
  `ModelingToolkit.get_ps(network)`.
"""
function params(network)
    isempty(get_systems(network)) ? get_ps(network) : parameters(network)
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, may allocate. Otherwise returns
  `ModelingToolkit.get_eqs(network)`.
"""
function reactions(network)
    isempty(get_systems(network)) ? get_eqs(network) : equations(network)
end

"""
    speciesmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping from species to
species indices. (Allocates)
"""
function speciesmap(network)
    Dict(S => i for (i,S) in enumerate(species(network)))
end

"""
    paramsmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping from parameters to
parameter indices. (Allocates)
"""
function paramsmap(network)
    Dict(p => i for (i,p) in enumerate(params(network)))
end

"""
    numspecies(network)

Return the number of species within the given [`ReactionSystem`](@ref).
"""
function numspecies(network)
    ns = length(get_states(network))
    for sys in get_systems(network)
        ns += numspecies(sys)
    end
    ns
end

"""
    numreactions(network)

Return the number of reactions within the given [`ReactionSystem`](@ref).
"""
function numreactions(network)
    nr = length(get_eqs(network))
    for sys in get_systems(network)
        nr += numreactions(sys)
    end
    nr
end

"""
    numparams(network)

Return the number of parameters within the given [`ReactionSystem`](@ref).
"""
function numparams(network)
    np = length(get_ps(network))
    for sys in get_systems(network)
        np += numparams(sys)
    end
    np
end

"""
    dependents(rx, network)

Given a [`Reaction`](@ref) and a [`ReactionSystem`](@ref), return a vector of
`ModelingToolkit.Num`s corresponding to 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.
"""
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(network, rxidx)
    dependents(network, rxidx)
end

"""
    substoichmat(rn; smap=speciesmap(rn))

Returns the substrate stoichiometry matrix
"""
function substoichmat(rn; smap=speciesmap(rn))
    smat = zeros(Int,(numreactions(rn),numspecies(rn)))
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.substoich
        for (i,sub) in enumerate(rx.substrates)
            smat[k,smap[sub]] = stoich[i]
        end
    end
    smat
end

"""
    prodstoichmat(rn; smap=speciesmap(rn))

Returns the product stoichiometry matrix
"""
function prodstoichmat(rn; smap=speciesmap(rn))
    pmat = zeros(Int,(numreactions(rn),numspecies(rn)))
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.prodstoich
        for (i,prod) in enumerate(rx.products)
            pmat[k,smap[prod]] = stoich[i]
        end
    end
    pmat
end

"""
    netstoichmat(rn; smap=speciesmap(rn))

Returns the net stoichiometry matrix
"""
function netstoichmat(rn; smap=speciesmap(rn))
    nmat = zeros(Int,(numreactions(rn),numspecies(rn)))
    for (k,rx) in pairs(reactions(rn))
        for (spec,coef) in rx.netstoich
            nmat[k,smap[spec]] = coef
        end
    end
    nmat
end

######################## some reacion complexes matrices and reaction rates ###############################
"""
$(TYPEDEF)
One reaction complex element

# Fields
$(FIELDS)
"""
struct ReactionComplexElement{T}
    """The integer id of the species representing this element."""
    speciesid::Int
    """The stoichiometric coefficient of this species."""
    speciesstoich::T
end

"""
$(TYPEDEF)
One reaction complex.

# Fields
$(FIELDS)
"""
struct ReactionComplex{V<:Integer} <: AbstractVector{ReactionComplexElement{V}}
    """The integer ids of all species participating in this complex."""
    speciesids::Vector{Int}
    """The stoichiometric coefficients of all species participating in this complex."""
    speciesstoichs::Vector{V}
end

function (==)(a::ReactionComplex{V},b::ReactionComplex{V}) where {V <: Integer}
    (a.speciesids == b.speciesids) &&
    (a.speciesstoichs == b.speciesstoichs)
end
hash(rc::ReactionComplex,h::UInt) = Base.hash(rc.speciesids,Base.hash(rc.speciesstoichs,h))
Base.size(rc::ReactionComplex) = size(rc.speciesids)
Base.length(rc::ReactionComplex) = length(rc.speciesids)
Base.getindex(rc::ReactionComplex, i...) =
        ReactionComplexElement(getindex(rc.speciesids, i...), getindex(rc.speciesstoichs, i...))
Base.setindex!(rc::ReactionComplex, t::ReactionComplexElement, i...) =
    (setindex!(rc.speciesids, t.speciesid, i...); setindex!(rc.speciesstoichs, t.speciesstoich, i...); rc)
Base.isless(a::ReactionComplexElement, b::ReactionComplexElement) = isless(a.speciesid, b.speciesid)
Base.Sort.defalg(::ReactionComplex{T}) where {T <: Integer} = Base.DEFAULT_UNSTABLE

"""
    reactioncomplexes(network, smap=speciesmap(rn))

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 -> ∅).
- The complex incidence matrix, B, is number of complexes by number of reactions with
    Bᵢⱼ = -1, if the i'th complex is the substrate of the j'th reaction,
           1, if the i'th complex is the product of the j'th reaction,
           0, otherwise
"""
function reactioncomplexes(rn; smap=speciesmap(rn))
    rxs = reactions(rn)
    numreactions(rn) > 0 || error("There must be at least one reaction to find reaction complexes.")
    complextorxsmap = OrderedDict{ReactionComplex{eltype(rxs[1].substoich)},Vector{Pair{Int,Int}}}()
    for (i,rx) in enumerate(rxs)
        reactantids = isempty(rx.substrates) ? Vector{Int}() : [smap[sub] for sub in rx.substrates]
        subrc = sort!(ReactionComplex(reactantids, copy(rx.substoich)))
        if haskey(complextorxsmap, subrc)
            push!(complextorxsmap[subrc], i => -1)
        else
            complextorxsmap[subrc] = [i => -1]
        end

        productids = isempty(rx.products) ? Vector{Int}() : [smap[prod] for prod in rx.products]
        prodrc = sort!(ReactionComplex(productids, copy(rx.prodstoich)))
        if haskey(complextorxsmap, prodrc)
            push!(complextorxsmap[prodrc], i => 1)
        else
            complextorxsmap[prodrc] = [i => 1]
        end
    end

    complexes = collect(keys(complextorxsmap))
    B = zeros(Int64, length(complexes), numreactions(rn));
    for (i,c) in enumerate(complexes)
        for (j,σ) in complextorxsmap[c]
            B[i,j] = σ
        end
    end
    complexes,B
end

"""
    reaction_rates(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

"""
    complexstoichmat(network; rcs=reactioncomplexes(rn)[1]))

Given a [`ReactionSystem`](@ref) and vector of reaction complexes, return a
matrix with positive entries of size num_of_species x num_of_complexes, where
the non-zero positive entries in the kth column denote stoichiometric
coefficients of the species participating in the kth reaction complex.
"""
function complexstoichmat(rn; rcs=reactioncomplexes(rn)[1])
    Z = zeros(Int64, numspecies(rn), length(rcs));
    for (i,rc) in enumerate(rcs)
        for rcel in rc
            Z[rcel.speciesid,i] = rcel.speciesstoich
        end
    end
    Z
end

"""
    complexoutgoingmat(network; B=reactioncomplexes(rn)[2])

Given a [`ReactionSystem`](@ref) and complex incidence matrix, B, return a matrix
of size num_of_complexes x num_of_reactions.

Notes
- The complex outgoing matrix, Δ, is defined by
    Δᵢⱼ = 0,    if Bᵢⱼ = 1
    Δᵢⱼ = Bᵢⱼ,  otherwise
"""
function complexoutgoingmat(rn; B=reactioncomplexes(rn)[2])
    Δ = copy(B)
    for (I,b) in pairs(Δ)
        (b == 1) && (Δ[I] = 0)
    end
    Δ
end

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

"""
    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::AbstractMatrix)::Matrix
    n_reac, n_spec = size(nsm)

    # We basically have to compute the left null space of the matrix
    # over the integers; this is best done using its Smith Normal Form.
    nsm_conv = AbstractAlgebra.matrix(AbstractAlgebra.ZZ, nsm)
    S, T, U = AbstractAlgebra.snf_with_transform(nsm_conv)

    # Zero columns of S (which occur after nonzero columns in SNF)
    # correspond to conserved quantities
    n = findfirst(i -> all(S[:,i] .== 0), 1:n_spec)
    if n === nothing
        return zeros(Int, 0, n_spec)
    end

    ret = Matrix(U[:,n:end]')

    # If all coefficients for a conservation law are negative
    # we might as well flip them to become positive
    for i in 1:size(ret,1)
        all(ret[i,:] .<= 0) && (ret[i,:] .*= -1)
    end

    ret
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 #######################

"""
    ==(rn1::Reaction, rn2::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 (==)(rn1::Reaction, rn2::Reaction)
    isequal(rn1.rate, rn2.rate) || return false
    issetequal(zip(rn1.substrates,rn1.substoich), zip(rn2.substrates,rn2.substoich)) || return false
    issetequal(zip(rn1.products,rn1.prodstoich), zip(rn2.products,rn2.prodstoich)) || return false
    issetequal(rn1.netstoich, rn2.netstoich)
end

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

Tests whether the underlying species, parameters and reactions are the same in
the two [`ReactionSystem`](@ref)s. Ignores order network components were
defined.

Notes:
- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
    considered different than `(A+1)^2`.
- Flattens subsystems, and hence may allocate, when checking equality.
"""
function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
    issetequal(species(rn1), species(rn2)) || return false
    issetequal(params(rn1), params(rn2)) || return false
    isequal(get_iv(rn1), get_iv(rn2)) || return false
    (numreactions(rn1) == numreactions(rn2)) || return false

    # the following fails for some reason, so need to use issubset
    #issetequal(equations(rn1), equations(rn2)) || return false
    (issubset(reactions(rn1),reactions(rn2)) && issubset(reactions(rn2),reactions(rn1))) || return false

    # BELOW SHOULD NOT BE NEEDED as species, params and equations flatten
    #issetequal(rn1.systems, rn2.systems) || return false
    # sys1 = rn1.systems; sys2 = rn2.systems
    # (issubset(sys1,sys2) && issubset(sys2,sys1)) || return false
    true
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)

    # 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)
        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

"""
    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)
    # 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)
        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)
    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`.
- Does not currently handle pins.
"""
function Base.merge!(network1::ReactionSystem, network2::ReactionSystem)
    isequal(get_iv(network1), get_iv(network2)) ||
        error("Reaction networks must have the same independent variable to be mergable.")
    union!(get_states(network1), get_states(network2))
    union!(get_ps(network1), get_ps(network2))
    append!(get_eqs(network1), get_eqs(network2))
    append!(get_systems(network1), get_systems(network2))
    network1
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.
- Does not currently handle pins.
"""
function Base.merge(network1::ReactionSystem, network2::ReactionSystem)
    network = make_empty_network()
    merge!(network, network1)
    merge!(network, network2)
    network
end


text:
# QuObj Tests
using Schrodinger
using Test, LinearAlgebra, SparseArrays
println("Testing QuObject Basics...")

@testset "QuObject Basics" begin

@testset "Ket Generation" begin
    # Test ket generation from arrays and functions
    @test data(Ket([0,0])) == [0.0,0.0]
    @test dims(Ket([0,1,0])) == (3,)
    @test_throws DimensionMismatch Ket([0,1],(3,))
    @test_throws ArgumentError basis(2,2)
    @test Ket(sparse([0,1/√2,1im/√2,0]),(2,2)) == 1/√2*(basis(2,0)⊗basis(2,1) + 1im*basis(2,1)⊗basis(2,0))
    @test Ket([0,1,0,0,0,0],(3,2)) == (basis(3,0)⊗basis(2,1))
    α = 1.31+0.11im
    @test coherent(30,α)[6] ≈ exp(-0.5*abs2(α))*α^5/sqrt(factorial(5))
    @test coherent(30,α,false) ≈ coherent(30,α,true)
    f() = qb"01"
    x = @inferred f()
    @test x == basis(2,0)⊗basis(2,1)
    @test data(ket((2,1,5,2),(5,3,10,4))).nzind[1] == 1+2+4*5+(4*10)*1+(4*10*3)*2
end

@testset "Bra Generation" begin
    # Test bra generation from arrays and functions
    @test data(Bra([0,0])) == [0.0,0.0]
    @test dims(Bra([0,1,0])) == (3,)
    @test_throws DimensionMismatch Bra([0,1],(3,))
    @test Bra(sparse([0,1/√2,1im/√2,0]),(2,2)) == 1/√2*(basis(2,0)⊗basis(2,1) - 1im*basis(2,1)⊗basis(2,0))'
    α = 1.31+0.11im
    @test Bra(coherent(30,α))[6] ≈ exp(-0.5*abs2(α))*conj(α)^5/sqrt(factorial(5))
end

@testset "Operator Generation" begin
    # Test operator generation from arrays and functions
    @test data(Operator([0 0;0 0])) == zeros(2,2)
    @test dims(Operator([0 0 0;0 1 0;0 0 0])) == (3,)
    @test_throws DimensionMismatch Operator([0 0;0 0],(3,))
    @test_throws DimensionMismatch Operator([0 1 0;0 0 0])
    @test (A=maxmixed(4); diag(A) == fill(0.25,4) && isdiag(A))
    @test (A=thermal(10,0.45); diag(A)[1:2] ≈ [6.89660888E-1,2.14032689E-1] && isdiag(A))
    @test data(qzero(6,(2,3))) == zeros(6,6)
    @test data(qeye(6,(2,3))) == Matrix{Float64}(I, 6, 6)
    @test data(qeye(6,(2,3))) == Matrix{Float64}(I, 6, 6)
    @test diag(create(6),-1) == [sqrt(n) for n = 1:5]
    @test destroy(6) == create(6)'
    @test (A=numberop(6); diag(A) == [n for n = 0:5] && isdiag(A))
    @test data(displacementop(3,0.5im)) ≈ [0.88261978 0.43980233im -0.16600070; 0.43980233im 0.64785934 0.62197442im; -0.16600070 0.62197442im 0.76523956]
    @test data(squeezeop(3,0.5im)) ≈ [0.93814834 0 -0.34623359im; 0 1 0; -0.34623359im 0 0.93814834]
    @test (U=rand_unitary(10,(5,2)); U*U' ≈ qeye(10,(5,2)))
end

@testset "Matrix Property Checks" begin
    # Tests matrix property checks
    H = Schrodinger.hermitianize!(rand(ComplexF64,5,5))
    H_s = Schrodinger.hermitianize!(sprand(ComplexF64,5,5,0.1))
    x = sparse([1,2,3,4,5],[3,5,3,1,1],[1+5im,3.4,1im,2.5+2.5im,6.6])
    @test isapproxhermitian(H+Matrix(x)*1E-15)
    @test !isapproxhermitian(H+Matrix(x)*1E-12)
    @test isapproxhermitian(H_s+x*1E-15)
    @test !isapproxhermitian(H_s+x*1E-12)
    U1 = sparse([1,2,4,3],[1,2,3,4],[1,1,1,1])
    U2 = expim(Hermitian(H))
    @test isunitary(U1)
    @test !isunitary(U1+sprand(4,4,0.5)*1E-12)
    @test isapproxunitary(U1)
    @test isapproxunitary(U1+sprand(4,4,0.5)*1E-15)
    @test !isapproxunitary(U1+sprand(4,4,0.5)*1E-12)
    @test isunitary(Matrix(U1))
    @test !isunitary(Matrix(U1+sprand(4,4,0.5)*1E-12))
    @test isapproxunitary(U2)
    @test !isapproxunitary(U2+Matrix(x))
end

@testset "Promotion" begin
    @test [σy', σx'] isa Vector{Operator{SparseMatrixCSC{ComplexF64,Int},1}}
    @test [σy', transpose(σy')] isa Vector{Operator{SparseMatrixCSC{ComplexF64,Int},1}}
    @test [sqrt(σx), sqrt(σx)'] isa Vector{Operator{Matrix{ComplexF64},1}}
    @test [sqrt(σx)', sqrt(σx)'] isa Vector{Operator{Adjoint{ComplexF64,Matrix{ComplexF64}},1}}
    @test [sqrt(σx), σx'] isa Vector{Operator{Matrix{ComplexF64},1}}
    @test [sqrt(σx'), σy] isa Vector{Operator{Matrix{ComplexF64},1}}
end

@testset "Display" begin
    @test sprint(show, "text/plain", qb"01"+qb"10") == """
    4-d $(Ket{SparseVector{Float64,Int64},2}) with dimensions 2⊗2
    1.00∠0°|0,1⟩ + 1.00∠0°|1,0⟩
    """
    @test sprint(show, normalize!(qb"01"+qb"10")) == "0.71∠0°|0,1⟩ + 0.71∠0°|1,0⟩"
    @test sprint(show, "text/plain", Bra([0,1,0])) == """
    3-d $(Bra{Vector{Float64},1}) with dimensions 3
    1.00∠0°⟨1|
    """
    @test sprint(show, Bra([0,0,1,0])) == "1.00∠0°⟨2|"
    sp = VERSION >= v"1.4-" ? "" : " "
    sz = VERSION >= v"1.6-" ? " ⋅ " : "0.0"
    @test sprint(show, "text/plain", maxmixed(4)) == """
    4×4 $(Operator{SparseMatrixCSC{Float64,Int64},1}) with dimensions 4
     0.25  $sz   $sz   $sz$sp
     $sz   0.25  $sz   $sz$sp
     $sz   $sz   0.25  $sz$sp
     $sz   $sz   $sz   0.25
    """
end
end


text:
# SD_AOU Functions for Keysight M32XXA/M33XXA AWGs
export SD_AOU_initialize, SD_AOU_clockGetFrequency, SD_AOU_clockGetSyncFrequency,
    SD_AOU_clockSetFrequency, SD_AOU_clockResetPhase, SD_AOU_channelWaveShape,
    SD_AOU_channelFrequency, SD_AOU_channelPhase, SD_AOU_channelPhaseReset,
    SD_AOU_channelPhaseResetMultiple, SD_AOU_channelAmplitude, SD_AOU_channelOffset,
    SD_AOU_modulationAngleConfig, SD_AOU_modulationAmplitudeConfig,
    SD_AOU_modulationIQconfig, SD_AOU_clockIOconfig, SD_AOU_triggerIOconfigV5,
    SD_AOU_triggerIOconfig, SD_AOU_triggerIOwrite, SD_AOU_triggerIOread,
    SD_AOU_waveformGetAddress, SD_AOU_waveformGetMemorySize,
    SD_AOU_waveformMemoryGetWriteAddress, SD_AOU_waveformMemorySetWriteAddress,
    SD_AOU_waveformReLoadArrayInt16, SD_AOU_waveformReLoad,
    SD_AOU_waveformLoadArrayInt16, SD_AOU_waveformLoad, SD_AOU_waveformFlush,
    SD_AOU_AWGqueueWaveform, SD_AOU_AWGstreamingConfig, SD_AOU_AWGstreamingRelease,
    SD_AOU_AWGidleValue, SD_AOU_AWGidleValueRead, SD_AOU_waveformLoadP2Pstart,
    SD_AOU_waveformLoadP2Pclose, SD_AOU_waveformLoadP2PgetPhysicalAddress,
    SD_AOU_waveformLoadP2PwriteDataTest, SD_AOU_AWGstart, SD_AOU_AWGstop,
    SD_AOU_AWGresume, SD_AOU_AWGpause, SD_AOU_AWGtrigger, SD_AOU_AWGjumpNextWaveform,
    SD_AOU_AWGjumpNextWaveformMultiple, SD_AOU_AWGisRunning, SD_AOU_AWGnWFplaying,
    SD_AOU_AWGstartMultiple, SD_AOU_AWGstopMultiple, SD_AOU_AWGresumeMultiple,
    SD_AOU_AWGpauseMultiple, SD_AOU_AWGtriggerMultiple, SD_AOU_AWGflush,
    SD_AOU_AWGfromFile, SD_AOU_AWGfromArray, SD_AOU_AWGtriggerExternalConfig,
    SD_AOU_AWGqueueConfig, SD_AOU_AWGqueueConfigRead, SD_AOU_AWGqueueMarkerConfig,
    SD_AOU_AWGqueueSyncMode

## int SD_AOU_initialize(int moduleID);
SD_AOU_initialize(moduleID::Integer) =
    ccall((:SD_AOU_initialize, klib), Cint, (Cint,), moduleID)

## double SD_AOU_clockGetFrequency(int moduleID);
"""
This function returns the real hardware clock frequency. It may differ from the
frequency set with the function clockSetFrequency, due to the hardware
frequency resolution.
"""
SD_AOU_clockGetFrequency(moduleID::Integer) =
    ccall((:SD_AOU_clockGetFrequency, klib), Cdouble, (Cint,), moduleID)

## double SD_AOU_clockGetSyncFrequency(int moduleID);
"""
This function returns the frequency of Clock System
"""
SD_AOU_clockGetSyncFrequency(moduleID::Integer) =
    ccall((:SD_AOU_clockGetSyncFrequency, klib), Cint, (Cint,), moduleID)

## double SD_AOU_clockSetFrequency(int moduleID, double frequency, int mode);
"""
This function sets the module clock frequency.
### CLK Set Frequency Mode
- Low Jitter Mode : The clock system is set to achieve the lowest jitter,
sacrificing tuning speed --> CLK_LOW_JITTER = 0
- Fast Tuning Mode : The clock system is set to achieve the lowest tuning time,
sacrificing jitter performance --> CLK_FAST_TUNE = 1
"""
SD_AOU_clockSetFrequency(moduleID::Integer, frequency::Real, mode::Integer) =
    ccall((:SD_AOU_clockSetFrequency, klib), Cdouble, (Cint, Cdouble, Cint),
        moduleID, frequency, mode)

## int SD_AOU_clockResetPhase(int moduleID, int triggerBehavior, int PXItrigger, double skew);
SD_AOU_clockResetPhase(moduleID::Integer, triggerBehavior::Integer,
    PXItrigger::Integer, skew::Real) =
    ccall((:SD_AOU_clockResetPhase, klib), Cint, (Cint, Cint, Cint, Cdouble),
        moduleID, triggerBehavior, PXItrigger, skew)

## int SD_AOU_channelWaveShape(int moduleID, int nChannel, int waveShape);
"""
This function sets the channel output waveform type.
### Output Signal Selection
- HIZ : The output signal is set to HIZ (no output signal is provided, only
available for M3202A) AOU_HIZ = -1
- No Signal : The output signal is set to 0. All other channel settings are
maintained --> AOU_OFF = 0 (default)
- Sinusoidal : Generated by the Function Generator --> AOU_SINUSOIDAL = 1
- Triangular : Generated by the Function Generator --> AOU_TRIANGULAR = 2
- Square : Generated by the Function Generator --> AOU_SQUARE = 4
- DC Voltage : The output DC voltage is set by the channel amplitude setting
--> AOU_DC = 5
- Arbitrary Waveform : Generated by the Arbitrary Waveform Generator
--> AOU_AWG = 6
- Partner Channel : Only for odd channels. It is the output of the previous
channel (to create differential signals, etc.) --> AOU_PARTNER = 8
"""
SD_AOU_channelWaveShape(moduleID::Integer, nChannel::Integer, waveShape::Integer) =
    ccall((:SD_AOU_channelWaveShape, klib), Cint, (Cint, Cint, Cint),
        moduleID, nChannel, waveShape)

## int SD_AOU_channelFrequency(int moduleID, int nChannel, double frequency);
"""
This function sets the frequency for the periodic signals generated by the
Function Generators.
"""
SD_AOU_channelFrequency(moduleID::Integer, nChannel::Integer, frequency::Real) =
    ccall((:SD_AOU_channelFrequency, klib), Cint, (Cint, Cint, Cdouble),
        moduleID, nChannel, frequency)

## int SD_AOU_channelPhase(int moduleID, int nChannel, double phase);
"""
This function sets the phase for the periodic signals generated by the
Function Generators.
"""
SD_AOU_channelPhase(moduleID::Integer, nChannel::Integer, phase::Real) =
    ccall((:SD_AOU_channelPhase, klib), Cint, (Cint, Cint, Cdouble),
        moduleID, nChannel, phase)

## int SD_AOU_channelPhaseReset(int moduleID, int channel);
"""
This function resets the accumulated phase of the selected signal. This
accumulated phase is the result of the phase continuous operation of the
product.
"""
SD_AOU_channelPhaseReset(moduleID::Integer, nChannel::Integer) =
    ccall((:SD_AOU_channelPhaseReset, klib), Cint, (Cint, Cint),
        moduleID, nChannel)

## int SD_AOU_channelPhaseResetMultiple(int moduleID, int channelMask);
"""
This function resets the accumulated phase of the selected channels
simultaneously. This accumulated phase is the result of the phase continuous
operation of the product.
"""
SD_AOU_channelPhaseResetMultiple(moduleID::Integer, channelMask::Integer) =
    ccall((:SD_AOU_channelPhaseResetMultiple, klib), Cint, (Cint, Cint),
        moduleID, channelMask)

## int SD_AOU_channelAmplitude(int moduleID, int nChannel, double amplitude);
"""
This function sets the amplitude of a channel.
"""
SD_AOU_channelAmplitude(moduleID::Integer, nChannel::Integer, amplitude::Real) =
    ccall((:SD_AOU_channelAmplitude, klib), Cint, (Cint, Cint, Cdouble),
        moduleID, nChannel, amplitude)

## int SD_AOU_channelOffset(int moduleID, int nChannel, double offset);
"""
This function sets the DC offset of a channel.
"""
SD_AOU_channelOffset(moduleID::Integer, nChannel::Integer, offset::Real) =
    ccall((:SD_AOU_channelOffset, klib), Cint, (Cint, Cint, Cdouble),
        moduleID, nChannel, offset)

## int SD_AOU_modulationAngleConfig(int moduleID, int nChannel, int modulationType, double deviationGain);
"""
This function configures the modulation in frequency/phase for the selected
channel.
### Angle Modulation Types
- No Modulation : Modulation is disabled --> AOU_MOD_OFF = 0 (default)
- Frequency Modulation : The AWG is used to modulate the channel frequency
--> AOU_MOD_FM = 1
- Frequency Modulation (32 bits): The AWG is used to modulate the channel
frequency. For models with -DM1 dual modulation option --> AOU_MOD_FM_32b = 1
- Phase Modulation : The AWG is used to modulate the channel phase
--> AOU_MOD_PM = 2
"""
SD_AOU_modulationAngleConfig(moduleID::Integer, nChannel::Integer,
    modulationType::Integer, deviationGain::Real) =
    ccall((:SD_AOU_modulationAngleConfig, klib), Cint, (Cint, Cint, Cint, Cdouble),
        moduleID, nChannel, modulationType, deviationGain)

## int SD_AOU_modulationAmplitudeConfig(int moduleID, int nChannel, int modulationType, double deviationGain);
"""
This function configures the modulation in amplitude/offset for the selected
channel.
"""
SD_AOU_modulationAmplitudeConfig(moduleID::Integer, nChannel::Integer,
    modulationType::Integer, deviationGain::Real) =
    ccall((:SD_AOU_modulationAmplitudeConfig, klib), Cint, (Cint, Cint, Cint,
        Cdouble), moduleID, nChannel, modulationType, deviationGain)

## int SD_AOU_modulationIQconfig(int moduleID, int nChannel, int enable);
"""
This function sets the IQ modulation for the selected channel.
### Amplitude Modulation Options
- No Modulation : Modulation is disabled. The channel amplitude and offset are
only set by the main registers set by the main registers
--> AOU_MOD_OFF = 0 (default)
- Amplitude Modulation : The modulating signal is used to modulate the channel
amplitude --> AOU_MOD_AM = 1
- Offset Modulation : The modulating signal is used to modulate the channel
offset --> AOU_MOD_OFFSET = 2
"""
SD_AOU_modulationIQconfig(moduleID::Integer, nChannel::Integer, enable::Integer) =
    ccall((:SD_AOU_modulationIQconfig, klib), Cint, (Cint, Cint, Cint),
        moduleID, nChannel, enable)

# Clock
## int SD_AOU_clockIOconfig(int moduleID, int clockConfig);
"""
This function configures the operation of the clock output connector.
### CLK Output Configurations
- Disable : The CLK connector is disabled --> CLK_CONN_DISABLED = 0 (default)
- CLKref Output : A copy of the reference clock is available at the CLK
connector --> CLK_REF_OUTPUT = 1
"""
SD_AOU_clockIOconfig(moduleID::Integer, clockConfig::Integer) =
    ccall((:SD_AOU_clockIOconfig, klib), Cint, (Cint, Cint), moduleID,
        clockIOconfig)

# Trigger
## int SD_AOU_triggerIOconfigV5(int moduleID, int direction, int syncMode);
SD_AOU_triggerIOconfigV5(moduleID::Integer, direction::Integer, syncMode::Integer) =
    ccall((:SD_AOU_triggerIOconfigV5, klib), Cint, (Cint, Cint, Cint),
        moduleID, direction, syncMode)

## int SD_AOU_triggerIOconfig(int moduleID, int direction);
"""
This function configures the trigger connector/line direction and
synchronization/sampling method
### Trigger I/O Options
- Trigger Output (readable) : TRG operates as a general purpose digital output
signal, which can be written by the user software --> AOU_TRG_OUT = 0
- Trigger Input : TRG operates as a trigger input, or as general purpose
digital input signal, which can be read by the user software --> AOU_TRG_IN = 1
### Trigger Synchronization/Sampling Options
- Non-synchronized mode : The trigger is sampled with an internal 100 MHz clock
--> SYNC_NONE = 0
- Synchronized mode : (PXI form factor only) The trigger is sampled using CLK10
--> SYNC_CLK_0 = 1
"""
SD_AOU_triggerIOconfig(moduleID::Integer, direction::Integer) =
    ccall((:SD_AOU_triggerIOconfig, klib), Cint, (Cint, Cint), moduleID, direction)

## int SD_AOU_triggerIOwrite(int moduleID, int value, int syncMode);
"""
This function sets the trigger output. The trigger must be configured as output
using function triggerIOconfig.
"""
SD_AOU_triggerIOwrite(moduleID::Integer, value::Integer, syncMode::Integer) =
    ccall((:SD_AOU_triggerIOwrite, klib), Cint, (Cint, Cint, Cint),
        moduleID, value, syncMode)

## int SD_AOU_triggerIOread(int moduleID);
"""
This function reads the trigger input.
"""
SD_AOU_triggerIOread(moduleID::Integer) =
    ccall((:SD_AOU_triggerIOread, klib), Cint, (Cint,), moduleID)

## int SD_AOU_waveformGetAddress(int moduleID, int waveformNumber);
SD_AOU_waveformGetAddress(moduleID::Integer, waveformNumber::Integer) =
    ccall((:SD_AOU_waveformGetAddress, klib), Cint, (Cint, Cint),
        moduleID, waveformNumber)

## int SD_AOU_waveformGetMemorySize(int moduleID, int waveformNumber);
SD_AOU_waveformGetMemorySize(moduleID::Integer, waveformNumber::Integer) =
    ccall((:SD_AOU_waveformGetMemorySize, klib), Cint, (Cint, Cint),
        moduleID, waveformNumber)

## int SD_AOU_waveformMemoryGetWriteAddress(int moduleID);
SD_AOU_waveformMemoryGetWriteAddress(moduleID::Integer) =
    ccall((:SD_AOU_waveformMemoryGetWriteAddress, klib), Cint, (Cint,), moduleID)

## int SD_AOU_waveformMemorySetWriteAddress(int moduleID, int writeAddress);
SD_AOU_waveformMemorySetWriteAddress(moduleID::Integer, writeAddress::Integer) =
    ccall((:SD_AOU_waveformMemorySetWriteAddress, klib), Cint, (Cint, Cint),
        moduleID, writeAddress)

## int SD_AOU_waveformReLoadArrayInt16(int moduleID, int waveformType, int waveformPoints, short *waveformDataRaw, int waveformNumber, int paddingMode SD_DEFAULT_NULL);
function SD_AOU_waveformReLoadArrayInt16(moduleID::Integer, waveformType::Integer,
    waveformDataRaw::Vector, waveformNumber::Integer, paddingMode::Integer=0)
    c_waveformDataRaw = Vector{Cshort}(waveformDataRaw)
    waveformPoints = length(c_waveformDataRaw)
    val = ccall((:SD_AOU_waveformReLoadArrayInt16, klib), Cint,
        (Cint, Cint, Cint, Ref{Cshort}, Cint, Cint), moduleID, waveformType,
        waveformPoints, c_waveformDataRaw, waveformNumber, paddingMode)
    return val
end

## int SD_AOU_waveformReLoad(int moduleID, int waveformID, int waveformNumber, int paddingMode SD_DEFAULT_NULL);
"""
This function replaces a waveform located in the module onboard RAM. The size of
the new waveform must be smaller than or equal to the existing waveform.
"""
SD_AOU_waveformReLoad(moduleID::Integer, waveformID::Integer, waveformNumber::Integer,
    paddingMode::Integer=0) =
    ccall((:SD_AOU_waveformReLoad, klib), Cint, (Cint, Cint, Cint, Cint),
        moduleID, waveformID, waveformNumber, paddingMode)

## int SD_AOU_waveformLoadArrayInt16(int moduleID, int waveformType, int waveformPoints, short *waveformDataRaw, int waveformNumber, int paddingMode SD_DEFAULT_NULL);
function SD_AOU_waveformLoadArrayInt16(moduleID::Integer, waveformType::Integer,
    waveformDataRaw::Vector, waveformNumber::Integer, paddingMode::Integer=0)
    c_waveformDataRaw = Vector{Cshort}(waveformDataRaw)
    waveformPoints = length(c_waveformDataRaw) # Not sure..
    val = ccall((:SD_AOU_waveformLoadArrayInt16, klib), Cint,
        (Cint, Cint, Cint, Ref{Cshort}, Cint, Cint), moduleID, waveformType,
        waveformPoints, c_waveformDataRaw, waveformNumber, paddingMode)
    return val
end

## int SD_AOU_waveformLoad(int moduleID, int waveformID, int waveformNumber, int paddingMode SD_DEFAULT_NULL);
"""
This function loads the specified waveform into the module onboard RAM.
Waveforms must be created first with the SD-Wave class.
"""
SD_AOU_waveformLoad(moduleID::Integer, waveformID::Integer,
    waveformNumber::Integer, paddingMode::Integer=0) =
    ccall((:SD_AOU_waveformLoad, klib), Cint, (Cint, Cint, Cint, Cint),
        moduleID, waveformID, waveformNumber, paddingMode)

## int SD_AOU_waveformFlush(int moduleID);
"""
This function deletes all the waveforms from the module onboard RAM and flushes
all the AWG queues.
"""
SD_AOU_waveformFlush(moduleID::Integer) =
    ccall((:SD_AOU_waveformFlush, klib), Cint, (Cint,), moduleID)

## int SD_AOU_AWGqueueWaveform(int moduleID, int nAWG, int waveformNumber, int triggerMode, int startDelay, int cycles, int prescaler);
"""
This function queues the specified waveform in one of the AWGs of the module.
The waveform must be already loaded in the module onboard RAM.
"""
SD_AOU_AWGqueueWaveform(moduleID::Integer, nAWG::Integer, waveformNumber::Integer,
    triggerMode::Integer, startDelay::Integer, cycles::Integer, prescaler::Integer) =
    ccall((:SD_AOU_AWGqueueWaveform, klib), Cint, (Cint, Cint, Cint, Cint, Cint,
        Cint, Cint), moduleID, nAWG, waveformNumber, triggerMode, startDelay,
        cycles, prescaler)

## int SD_AOU_AWGstreamingConfig(int moduleID, int nAWG, int triggerMode, int startDelay, int prescaler, int waveformPoints, callbackWGPtr callbackFunction, void *callbackUserObj);
SD_AOU_AWGstreamingConfig(moduleID::Integer, nAWG::Integer, triggerMode::Integer,
    startDelay::Integer, prescaler::Integer, waveformPoints::Integer) =
    ccall((:SD_AOU_AWGstreamingConfig, klib), Cint,
    (Cint, Cint, Cint, Cint, Cint, Cint, Ptr{Void}, Ptr{Void}),
    moduleID, nAWG, triggerMode, startDelay, prescaler, waveformPoints,
    C_NULL, C_NULL)

## int SD_AOU_AWGstreamingRelease(int moduleID, int nAWG);
SD_AOU_AWGstreamingRelease(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGstreamingRelease, klib), Cint, (Cint, Cint), moduleID,
        nAWG)

## int SD_AOU_AWGidleValue(int moduleID, int nAWG, double value);
SD_AOU_AWGidleValue(moduleID::Integer, nAWG::Integer, value::Real) =
    ccall((:SD_AOU_AWGidleValue, klib), Cint, (Cint, Cint, Cdouble),
        moduleID, nAWG, value)

## double SD_AOU_AWGidleValueRead(int moduleID, int nAWG);
SD_AOU_AWGidleValueRead(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGidleValueRead, klib), Cdouble, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_waveformLoadP2Pstart(int moduleID, int wfType, int points, int waveformNumber);
SD_AOU_waveformLoadP2Pstart(moduleID::Integer, wfType::Integer, points::Integer,
    waveformNumber::Integer) =
    ccall((:SD_AOU_waveformLoadP2Pstart, klib), Cint, (Cint, Cint, Cint, Cint),
        moduleID, wfType, points, waveformNumber)

## int SD_AOU_waveformLoadP2Pclose(int moduleID, int nDWwritten);
SD_AOU_waveformLoadP2Pclose(moduleID::Integer, nDWwritten::Integer) =
    ccall((:SD_AOU_waveformLoadP2Pclose, klib), Cint, (Cint, Cint), moduleID,
        nDWwritten)

## unsigned long long SD_AOU_waveformLoadP2PgetPhysicalAddress(int moduleID);
SD_AOU_waveformLoadP2PgetPhysicalAddress(moduleID::Integer) =
    ccall((:SD_AOU_waveformLoadP2PgetPhysicalAddress, klib), Culonglong,
        (Cint,), moduleID)

## int SD_AOU_waveformLoadP2PwriteDataTest(int moduleID, int wfType, int points, short *dataRaw, int paddingMode SD_DEFAULT_NULL);
function SD_AOU_waveformLoadP2PwriteDataTest(moduleID::Integer, wfType::Integer,
    dataRaw::Vector, paddingMode::Integer)
    c_dataRaw = Vector{Cshort}(dataRaw)
    points = length(c_dataRaw)
    val = ccall((:SD_AOU_waveformLoadP2PwriteDataTest, klib), Cint,
        (Cint, Cint, Cint, Ref{Cshort}, Cint), moduleID, wfType, points,
        c_dataRaw, paddingMode)
    return val
end

## int SD_AOU_AWGstart(int moduleID, int nAWG);
"""
This function starts he selected AWG from the beginning of its queue. The
generation will start immediately or when a trigger is received, depending on
the trigger selection of the first waveform in the queue and provided that at
least one waveform is queued in the AWG.
"""
SD_AOU_AWGstart(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGstart, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGstop(int moduleID, int nAWG);
"""
This function stops the selected AWG, setting the output to zero and resetting
the AWG queue to its initial position. All following incoming triggers are
ignored.
"""
SD_AOU_AWGstop(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGstop, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGresume(int moduleID, int nAWG);
"""
The function resumes the operation of the selected AWG from the current
position of the queue.
"""
SD_AOU_AWGresume(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGresume, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGpause(int moduleID, int nAWG);
"""
This function pauses the selected AWG, leaving the last waveform point at the
output, and ignoring all incoming triggers. The waveform generation can be
resumed calling AWGresume.
"""
SD_AOU_AWGpause(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGpause, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGtrigger(int moduleID, int nAWG);
"""
This function triggers the selected AWG. The waveform waiting in the current
position of the queue is launched provided it is configured with VI/HVI Trigger.
"""
SD_AOU_AWGtrigger(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGtrigger, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGjumpNextWaveform(int moduleID, int nAWG);
"""
This function forces a jump to the next waveform in the AWG queue. The jump is
executed once the current waveform has finished a complete cycle.
"""
SD_AOU_AWGjumpNextWaveform(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGjumpNextWaveform, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGjumpNextWaveformMultiple(int moduleID, int AWGmask);
"""
This function forces a jump to the next waveform in the queue of several AWGs.
The jumps are executed once the current waveforms have finished a complete
cycle.
"""
SD_AOU_AWGjumpNextWaveformMultiple(moduleID::Integer, AWGmask::Integer) =
    ccall((:SD_AOU_AWGjumpNextWaveformMultiple, klib), Cint, (Cint, Cint),
        moduleID, AWGmask)

## int SD_AOU_AWGisRunning(int moduleID, int nAWG);
"""
This function returns if the AWG is running or stopped.
"""
SD_AOU_AWGisRunning(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGisRunning, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGnWFplaying(int moduleID, int nAWG);
"""
This function returns the waveformNumber of the waveform which is currently
being generated.
"""
SD_AOU_AWGnWFplaying(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGnWFplaying, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGstartMultiple(int moduleID, int AWGmask);
"""
This function starts the selected AWGs from the beginning of their queues. The
generation will start immediately or when a trigger is received, depending on
the trigger selection of the first waveform in their queues and provided that at
least one waveform is queued in these AWGs (functions AWGqueueWaveform, or AWG).
"""
SD_AOU_AWGstartMultiple(moduleID::Integer, AWGmask::Integer) =
    ccall((:SD_AOU_AWGstartMultiple, klib), Cint, (Cint, Cint),
        moduleID, AWGmask)

## int SD_AOU_AWGstopMultiple(int moduleID, int AWGmask);
"""
This function stops the selected AWGs, setting their outputs to zero and
resetting their respective queues to the initial positions. All following
incoming triggers are ignored.
"""
SD_AOU_AWGstopMultiple(moduleID::Integer, AWGmask::Integer) =
    ccall((:SD_AOU_AWGstopMultiple, klib), Cint, (Cint, Cint), moduleID, AWGmask)

## int SD_AOU_AWGresumeMultiple(int moduleID, int AWGmask);
"""
This function resumes the operation of the selected AWGs from the current
position of their respective queues.
"""
SD_AOU_AWGresumeMultiple(moduleID::Integer, AWGmask::Integer) =
    ccall((:SD_AOU_AWGresumeMultiple, klib), Cint, (Cint, Cint), moduleID,
        AWGmask)

## int SD_AOU_AWGpauseMultiple(int moduleID, int AWGmask);
"""
This function pauses the selected AWGs, leaving the last waveform point at the
output, and ignoring all incoming triggers. The waveform generation can be
resumed calling AWGresume.
"""
SD_AOU_AWGpauseMultiple(moduleID::Integer, AWGmask::Integer) =
    ccall((:SD_AOU_AWGpauseMultiple, klib), Cint, (Cint, Cint), moduleID,
        AWGmask)

## int SD_AOU_AWGtriggerMultiple(int moduleID, int AWGmask);
"""
This function triggers the selected AWGs. The waveform waiting in the current
position of the queue is launched provided it is configured with VI/HVI Trigger.
"""
SD_AOU_AWGtriggerMultiple(moduleID::Integer, AWGmask::Integer) =
    ccall((:SD_AOU_AWGtriggerMultiple, klib), Cint, (Cint, Cint), moduleID,
        AWGmask)

## int SD_AOU_AWGflush(int moduleID, int nAWG);
"""
This function empties the queue of the selected AWGs. Waveforms are not removed
from the module onboard RAM.
"""
SD_AOU_AWGflush(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGflush, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGfromFile(int moduleID, int nAWG, const char *waveformFile, int triggerMode, int startDelay, int cycles, int prescaler, int paddingMode SD_DEFAULT_NULL);
SD_AOU_AWGfromFile(moduleID::Integer, nAWG::Integer, waveformFile::String,
    triggerMode::Integer, startDelay::Integer, cycles::Integer,
    prescaler::Integer, paddingMode::Integer) =
    ccall((:SD_AOU_AWGfromFile, klib), Cint, (Cint, Cint, Cstring, Cint, Cint,
        Cint, Cint, Cint), moduleID, nAWG, waveformFile, triggerMode, startDelay,
        cycles, prescaler, paddingMode)
# not sure if waveformFile is a pointer to the file or just a string of filename

## int SD_AOU_AWGfromArray(int moduleID, int nAWG, int triggerMode, int startDelay, int cycles, int prescaler, int waveformType, int waveformPoints, double* waveformDataA, double* waveformDataB SD_DEFAULT_NULL, int paddingMode SD_DEFAULT_NULL);
"""
This function provides a one-step method to load, queue and start a single
waveform in one of the module AWGs. The waveform can be loaded from an array of
points in memory or from a file.
### AWG Trigger Mode
- Auto : The waveform is launched automatically after function AWGstart, or
when the previous waveform in the queue finishes --> AUTOTRIG = 0
- Software / HVI : Software trigger. The AWG is triggered by the function
AWGtrigger, provided that the AWG is running. AWGtrigger can be executed from
the user application (VI) or from an HVI --> SWHVITRIG = 1
- Software / HVI (per cycle) : Software trigger. Identical to the previous
option, but the trigger is required per each waveform cycle
--> SWHVITRIG_CYCLE = 5
- External Trigger : Hardware trigger. The AWG waits for an external trigger
--> EXTTRIG = 2
- External Trigger (per cycle) : Hardware trigger. Identical to the previous
option, but the trigger is required per each waveform cycle --> EXTTRIG_CYCLE = 6
"""
function SD_AOU_AWGfromArray(moduleID::Integer, nAWG::Integer, triggerMode::Integer,
    startDelay::Integer, cycles::Integer, prescaler::Integer, waveformType::Integer,
    waveformDataA::Vector, waveformDataB::Vector=0, paddingMode::Integer=0)
    waveformPoints = length(waveformDataA)
    c_waveformDataA = Vector{Cdouble}(waveformDataA)
    c_waveformDataB = (waveformDataB==0)? 0: Vector{Cdouble}(waveformDataB)
    val = ccall((:SD_AOU_AWGfromArray, klib), Cint, (Cint, Cint, Cint, Cint,
        Cint, Cint, Cint, Cint, Ref{Cdouble}, Ref{Cdouble}), moduleID, nAWG,
        triggerMode, startDelay, cycles, prescaler, waveformType,
        waveformPoints, c_waveformDataA, c_waveformDataB, paddingMode)
    return val
end

## int SD_AOU_AWGtriggerExternalConfig(int moduleID, int nAWG, int externalSource, int triggerBehavior, int sync);
"""
This function configures the external triggers for the selected AWG. The
external trigger is used in case the waveform is queued with the external
trigger mode option (function AWGqueueWaveform).
### AWG External Trigger Source
- External I/O Trigger : The AWG trigger is a TRG connector/line of the module
I/O Triggers. PXI form factor only: this trigger can be synchronized to CLK10
--> TRIG_EXTERNAL = 0
- PXI Trigger [0..n] : PXI form factor only. The AWG external trigger is a PXI
trigger line and it is synchronized to CLK10 --> TRIG_PXI_AWG(=4000) + Trigger No.
### AWG External Trigger Behavior
- Active High : Trigger is active when it is at level high --> TRIG_HIGH = 1
- Active Low : Trigger is active when it is at level low --> TRIG_LOW = 2
- Rising Edge : Trigger is active on the rising edge --> TRIG_RISE = 3
- Falling Edge : Trigger is active on the falling edge --> TRIG_FALL = 4
"""
SD_AOU_AWGtriggerExternalConfig(moduleID::Integer, nAWG::Integer, externalSource::Integer,
    triggerBehavior::Integer, sync::Integer) =
    ccall((:SD_AOU_AWGtriggerExternalConfig, klib), Cint, (Cint, Cint, Cint,
        Cint, Cint), moduleID, nAWG, externalSource, triggerBehavior, sync)

## int SD_AOU_AWGqueueConfig(int moduleID, int nAWG, int mode);
"""
This function configures the cyclic mode of the queue. All waveforms must be
already queued in one of the AWGs of the module.
"""
SD_AOU_AWGqueueConfig(moduleID::Integer, nAWG::Integer, mode::Integer) =
    ccall((:SD_AOU_AWGqueueConfig, klib), Cint, (Cint, Cint, Cint), moduleID,
        nAWG, mode)

## int SD_AOU_AWGqueueConfigRead(int moduleID, int nAWG);
"""
This function reads the value of a cyclic mode of the queue. All waveforms must
be already queued in one of the AWGs of the module.
"""
SD_AOU_AWGqueueConfigRead(moduleID::Integer, nAWG::Integer) =
    ccall((:SD_AOU_AWGqueueConfigRead, klib), Cint, (Cint, Cint), moduleID, nAWG)

## int SD_AOU_AWGqueueMarkerConfig(int moduleID, int nAWG, int markerMode, int trgPXImask, int trgIOmask, int value, int syncMode, int length, int delay);
"""
This function configures the Marker generation for each AWG. All waveforms must
be already queued in one of the AWGs of the module.
"""
SD_AOU_AWGqueueMarkerConfig(moduleID::Integer, nAWG::Integer, markerMode::Integer,
    trgPXImask::Integer, trgIOmask::Integer, value::Integer, syncMode::Integer,
    length::Integer, delay::Integer) =
    ccall((:SD_AOU_AWGqueueMarkerConfig, klib), Cint,
        (Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
        moduleID, nAWG, markerMode, trgPXImask, trgIOmask, value, syncMode,
        length, delay)

## int SD_AOU_AWGqueueSyncMode(int moduleID, int nAWG, int syncMode);
"""
This function configures the sync mode of the queue. All waveforms must be
already queued in one of the AWGs of the module.
"""
SD_AOU_AWGqueueSyncMode(moduleID::Integer, nAWG::Integer, syncMode::Integer) =
    ccall((:SD_AOU_AWGqueueSyncMode, klib), Cint, (Cint, Cint, Cint),
        moduleID, nAWG, syncMode)

# Functions that exist in the manual, but not in the library `SD_AIN.h` file
# ## AWGreset
# """
# This function resets the selected AWG to its initial position.
# """
# SD_AOU_AWGreset(moduleID::Integer, nAWG::Integer) =
#     ccall((:SD_AOU_AWGreset, klib), Cint, (Cint, Cint), moduleID, nAWG)


text:
@polly function kernel_lu(A)
    n = size(A,1)
    for i = 1:n
        for j = 1:(i-1)
            for k = 1:(j-1)
                A[i,j] -= A[i,k] * A[k,j]
            end
            A[i,j] /= A[j,j]
        end
        for j = i:n, k = 1:(i-1)
            A[i,j] -= A[i,k] * A[k,j]
        end
    end
end

let
    n = 2000

    A = zeros(Float32, n, n)

    for i = 1:n
        for j = 1:i
            A[i,j] = (-j % n) / n + 1
        end
        for j = i+1:n
            A[i,j] = 0
        end
        A[i,i] = 1
    end

    # Make the matrix positive semi-definite.
    # not necessary for LU, but using same code as cholesky
    B = zeros(Float32, n, n)
    for r = 1:n, s = 1:n
        B[r,s] = zero(eltype(B))
    end
    for t = 1:n
        for r = 1:n, s = 1:n
            B[r,s] += A[r,t] * A[s,t]
        end
        for r = 1:n, s = 1:n
            A[r,s] = B[r,s]
        end
    end

    SUITE["lu"] = @benchmarkable kernel_lu(A) setup = (A = copy($A))
end


text:
"""
    GitHubAPI
Module used by the Social Decision and Analytics Division (SDAD) of
the Biocomplexity Institute and Initiative of the University of Virginia.

This module was designed for the Open-Source Software (OSS) project.
"""
module GitHubAPI
    using Base: include_dependency
    using ConfParser: ConfParse, parse_conf!, retrieve, save!, commit!
    using Diana: GraphQLClient,
                 # HTTP
                 HTTP.Messages.Response, HTTP.request
    using JSON3: JSON3
    using LibPQ: Connection, execute
    using Parameters: @unpack
    using TimeZones: localzone, now, today, yearmonthday, ZonedDateTime,
                     # Dates
                     @dateformat_str, Dates.format
    # Compat
    if !@isdefined(isnothing)
        isnothing(::Any) = false
        isnothing(::Nothing) = true
    end
    if !@isdefined(ismissing)
        isnothing(::Any) = false
        isnothing(::Missing) = true
    end
    foreach(include,
            file for file ∈ readdir(@__DIR__) if !isequal("GitHubAPI.jl", file))
    export get_licenses, dbconnect, parse_repos!, parse_commits!, sdad_setup!
    include_dependency(joinpath(dirname(@__FILE__), "..", "confs", "config.simple"))
end


text:
module Day17

using AdventOfCode2019

function day17(input::String = readInput(joinpath(@__DIR__, "..", "data", "day17.txt")))
    program = parse.(Int, split(input, ','))

    view = read_camera!(copy(program))
    intersections = find_intersections(view)
    part1 = ((x[1]-1) * (x[2]-1) for x in intersections) |> sum

    instructions = full_instruction(view)
    compressed = split_instructions(instructions)
    part2 = control_robot!(compressed, program)
    return [part1, part2]
end

function read_camera!(program::Array{Int,1})
    view = Array{Int,1}()
    out = Channel{Int}(5)
    t = @async AdventOfCode2019.IntCode.run_program!(program, nothing, out, nothing)
    while !istaskdone(t)
        push!(view, take!(out))
    end
    close(out)

    rows = split((strip(join(Char.(view)))),'\n')
    mat = Array{Char,2}(undef, length(rows), length(rows[1]))
    for (i, row) in enumerate(rows)
        for (j, char) in enumerate(row)
            mat[i, j] = char
        end
    end
    return mat
end

function find_intersections(view)
    intersections = Array{CartesianIndex{2},1}()
    for i = 2:size(view, 1) - 1
        for j = 2:size(view, 2) - 1
            if view[i,j] == '#'
                if view[i-1,j] == '#' && view[i+1,j] == '#' && view[i,j-1] == '#' && view[i,j+1] == '#'
                    push!(intersections, CartesianIndex(i, j))
                end
            end
        end
    end
    return intersections
end

function full_instruction(view)
    instructions = Array{Union{Char,Int},1}()
    asciiDirs = ['^', 'v', '<', '>']
    start = findfirst(x -> x in asciiDirs, view)

    direction = 0
    for i = 1:4
        if view[start] == asciiDirs[i]
            direction = i
            break
        end
    end

    neighbours = []
    d = direction
    for i = 1:4
        push!(neighbours, view[(start.I .+ dir_1d_to_2d(d))...])
        d = turn_left(d)
    end
    nLeftTurns = findfirst(x -> x == '#', neighbours) - 1
    if nLeftTurns == 3
        push!(instructions, 'R')
    elseif nLeftTurns == 2
        throw(AssertionError("not solved for this case"))
        push!(instructions, 'L')
        push!(instructions, 'L')
    elseif nLeftTurns == 1
        push!(instructions, 'L')
    end
    for i = 1:nLeftTurns
        direction = turn_left(direction)
    end

    i, j = start.I
    while true
        k, l = (i, j) .+ dir_1d_to_2d(direction)
        if inbounds(k, l, view) && view[k, l] == '#'
            push!(instructions, 1)
            i, j = k, l
            continue
        end
        k, l = (i, j) .+ dir_1d_to_2d(turn_left(direction))
        if inbounds(k, l, view) && view[k, l] == '#'
            push!(instructions, 'L')
            direction = turn_left(direction)
            continue
        end
        k, l = (i, j) .+ dir_1d_to_2d(turn_right(direction))
        if inbounds(k, l, view) && view[k, l] == '#'
            push!(instructions, 'R')
            direction = turn_right(direction)
            continue
        end
        break
    end

    # simplify instructions
    simplifiedInstructions = Array{Union{Char,Int},1}()
    count = 0
    for (i, inst) in enumerate(instructions)
        if inst in ('L', 'R')
            push!(simplifiedInstructions, inst)
        else
            count += inst
            if i + 1 > length(instructions) || isa(instructions[i+1], Char)
                push!(simplifiedInstructions, count)
                count = 0
            end
        end
    end

    instructionPairs = Array{Pair{Char,Int},1}()
    for i = 1:2:length(simplifiedInstructions) - 1
        push!(instructionPairs, Pair(simplifiedInstructions[i], simplifiedInstructions[i+1]))
    end

    return instructionPairs
end

function split_instructions(pairs::Array{Pair{Char,Int},1})
    main, a, b, c = "", "", "", ""
    for lenA = 2:10
        s = join(join.(pairs, ","), ",")
        rega = Regex("^([RL],[0-9]+,?){$lenA}")
        ra = findfirst(rega, s)
        isa(ra, Nothing) && continue
        a = strip(s[ra], ',')
        sa = replace(s, a => "A")
        for lenB = 2:10
            regb = Regex("([RL],[0-9]+,?){$lenB}")
            rb = findall(regb, sa)
            for r in rb
                b = strip(sa[r], ',')
                sb = replace(sa, b => "B")
                for lenC = 2:10
                    regc = Regex("([RL],[0-9]+,?){$lenC}")
                    rc = findall(regc, sb)
                    for t in rc
                        c = strip(sb[t], ',')
                        sc = replace(sb, c => "C")
                        finalreg = r"^([ABC],?){2,10}+$"
                        if !isa(findfirst(finalreg, sc), Nothing)
                            main = sc
                            @goto found
                        end
                    end
                end
            end
        end
    end
    throw(AssertionError("Unable to compress instructions"))
    @label found
    return Dict(
        "Main" => main * "\n",
        "A" => a * "\n",
        "B" => b * "\n",
        "C" => c * "\n",
    )
end

function control_robot!(compressedInstruction, program)
    inp = Channel{Int}(100)
    out = Channel{Int}(1)
    program[1] == 1 || throw(AssertionError("Invalid program"))
    program[1] = 2

    main = [Int(compressedInstruction["Main"][i]) for i in eachindex(compressedInstruction["Main"])]
    a = [Int(compressedInstruction["A"][i]) for i in eachindex(compressedInstruction["A"])]
    b = [Int(compressedInstruction["B"][i]) for i in eachindex(compressedInstruction["B"])]
    c = [Int(compressedInstruction["C"][i]) for i in eachindex(compressedInstruction["C"])]
    vidS = "n\n"
    video = [Int(vidS[i]) for i in eachindex(vidS)]

    t = @async AdventOfCode2019.IntCode.run_program!(program, inp, out, nothing)

    for val in main
        put!(inp, val)
    end
    for val in a
        put!(inp, val)
    end
    for val in b
        put!(inp, val)
    end
    for val in c
        put!(inp, val)
    end
    for val in video
        put!(inp, val)
    end

    while true
        ret = take!(out)
        if istaskdone(t)
            close(inp)
            close(out)
            return ret
        end
    end
end

function inbounds(i::Int, j::Int, view)
    return i >= 1 && j >= 1 && i <= size(view, 1) && j <= size(view, 2)
end

function turn_left(d::Int)
    d == 1 && return 3
    d == 2 && return 4
    d == 3 && return 2
    d == 4 && return 1
end

function turn_right(d::Int)
    d == 1 && return 4
    d == 2 && return 3
    d == 3 && return 1
    d == 4 && return 2
end

function dir_1d_to_2d(d::Int)
    d == 1 && return (-1, 0)
    d == 2 && return (1, 0)
    d == 3 && return (0, -1)
    d == 4 && return (0, 1)
end

function dir_2d_to_1d(d::Tuple{Int,Int})
    d == (-1, 0) && return 1
    d == (1, 0) && return 2
    d == (0, -1) && return 3
    d == (0, 1) && return 4
end

end # module


text:
using Knet, ArgParse, JLD
include("bilstm.jl")
include("infst.jl")

function main(args=ARGS)
    s = ArgParseSettings()
    s.exc_handler = ArgParse.debug_handler
    @add_arg_table s begin
        ("--trainfile"; required=true; help="Infinite stream training file")
        ("--devfile"; required=true; help="dev data to test perplexity")
        ("--vocabfile"; required=true; help="Vocabulary file to train a model")
        ("--savefile"; help="To save the julia model")
        ("--atype"; default=(gpu() >= 0 ? "KnetArray{Float32}" : "Array{Float32}"))
        ("--hiddens"; arg_type=Int; nargs='+'; default=[300]; help="hidden layer configuration")
        ("--embedding"; arg_type=Int; default=256)
        ("--batchsize"; arg_type=Int; default=23)
        ("--gclip"; arg_type=Float64; default=5.0)

    end
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    atype = eval(parse(o[:atype]))
    for (k, v) in o
        println("$k => $v")
    end

    # prepare train data infinite stream
    ulimit = 31
    maxlines = 500
    fdir = open(o[:trainfile])
    vocab = create_vocab(o[:vocabfile])
    sdict = Dict{Int64, Array{Any, 1}}();
    readstream!(fdir, sdict, vocab; maxlines=1000, ulimit=ulimit)

    # initialize model
    vsize = length(vocab)
    m = initmodel(atype, o[:hiddens], o[:embedding], vsize)
    s = initstate(atype, o[:hiddens], o[:batchsize])
    opts = oparams(m, Adam; gclip=o[:gclip])

    # prepare test data
    dev = create_testdata(o[:devfile], vocab, 5)
    sdev = initstate(atype, o[:hiddens], 5)

    dloss = devperp(m, sdev, dev)
    println("Initial dev perplexity is $dloss")

    # training loop
    ids = nextbatch(fdir, sdict, vocab, o[:batchsize]; ulimit=ulimit, maxlines=maxlines)
    checkpoint = 0
    info("Training started")
    while ids != nothing
        train(m, s, ids, opts)
        ids = nextbatch(fdir, sdict, vocab, o[:batchsize]; ulimit=ulimit, maxlines=maxlines)
        if checkpoint % 8000 == 0
            checkpoint = 1
            dloss = devperp(m, sdev, dev)
            println("Dev perp is $dloss")
            moc = convertmodel(m)
            save(o[:savefile], "model", moc)
        end
        checkpoint += 1
    end
end
!isinteractive() && main(ARGS)

max_stars_repo_name: kirnap/wordsubs

max_issues_repo_name: kirnap/wordsubs

max_forks_repo_name: kirnap/wordsubs


text:
using MultiResolutionFilters
using ImageFiltering
using RegionTrees
using StaticArrays
using Statistics, StatsBase
using Distributions
using StatsPlots
using Random

Random.seed!(1)
Z = 1.5imfilter(randn(500, 500), Kernel.gaussian(20)) #.+ 0.02randn(500, 500)
Z = Z .+ (1:500) ./ 2e3

n = 1_000
x = [60randn(2, n) .+ 150  75randn(2, n) .+ 350]
x = filter(xi -> all(1 .< xi .<= 500), collect(eachcol(x)))
z = [Z[round(Int, xi[2]), round(Int, xi[1])] for xi in x]
z = z + 0.05randn(length(z))
# x = hcat(x...)
x = SVector{2}.(x)

zl = extrema(Z)
plt_truth = heatmap(Z, clim=zl, legend=false, xlabel="X", ylabel="Y")
scatter!(plt_truth, first.(x), last.(x), color=:black, markersize=2, label="")

scatter(first.(x), last.(x), zcolor=z, markerstrokewidth=0, legend=false, xlabel="X", ylabel="Y")

nparticles = 2000
state_particles = randn(nparticles)
state_weights = ones(nparticles) / nparticles
x0 = [0., 0.]
width = [500., 500]

dx = 500 ./ 2.0.^(0:10)
dA = dx.^2
plot(dx)
plot(dA)

loglik(state, obs) = loglikelihood(Normal(state, 0.1), obs)
predict_downscale(state, child_area) = state + 0.5randn()
observe(state) = state + randn()
needsrefinement(cell) = area(cell.boundary) > dA[6]
jitter(state) = state #+ 0.05randn()

r = ParticleRefinery(predict_downscale, observe, loglik, needsrefinement)
mrf = MultiResolutionPF(x, z, r, state_particles)
adaptive_filter!(mrf)
tree = mrf.tree

p0 = plot(legend=false)
v = hcat(collect(vertices(tree.boundary))...)
plot!(p0, v[1,[1,2,4,3,1]], v[2,[1,2,4,3,1]])
c1 = children(tree)
p1 = plot(legend=false)
for c in c1
    v = hcat(collect(vertices(c.boundary))...)
    plot!(p1, v[1,[1,2,4,3,1]], v[2,[1,2,4,3,1]], color=2)
end
p1
p2 = plot(legend=false)
c2 = reduce(vcat, [children(c) for c in c1])
for c in c2
    v = hcat(collect(vertices(c.boundary))...)
    plot!(p2, v[1,[1,2,4,3,1]], v[2,[1,2,4,3,1]], color=3)
end
c3 = reduce(vcat, [children(c) for c in c2])
p3 = plot(legend=false)
for c in c3
    v = hcat(collect(vertices(c.boundary))...)
    plot!(p3, v[1,[1,2,4,3,1]], v[2,[1,2,4,3,1]], color=4)
end
plot(p0, p1, p2, p3)

μs = [mean(leaf.data.state_particles) for leaf in allleaves(tree)]
histogram(μs)

cg = cgrad(:magma);
plt_post = plot(legend=nothing, xlabel="X", ylabel="Y");
for leaf in allleaves(tree)
    xl, yl = leaf.boundary.origin
    w, h = leaf.boundary.widths
    μ = mean(leaf.data.state_particles)
    c = (μ - zl[1]) / (zl[2] - zl[1])
    # nobs(leaf.data) > 0 ? linecolor = :black : linecolor = get(cg, c)
    linecolor = get(cg, c)
    plot!(plt_post, Shape(xl .+ [0,w,w,0], yl .+ [0,0,h,h]), fill=get(cg, c),
        linecolor=linecolor)
end
plt_post
scatter!(plt_post, first.(x), last.(x), color=:black, markersize=2)

cg = cgrad(:viridis);
cl = extrema([std(leaf.data.state_particles) for leaf in allleaves(tree)])
plt_std = plot(legend=nothing, xlabel="X", ylabel="Y");
for leaf in allleaves(tree)
    # if leaf.data.n > 0
        xl, yl = leaf.boundary.origin
        w, h = leaf.boundary.widths
        σ = std(leaf.data.state_particles)
        c = (σ - cl[1]) / (cl[2] - cl[1])
        plot!(plt_std, Shape(xl .+ [0,w,w,0], yl .+ [0,0,h,h]), fill=get(cg, c),
            linecolor=get(cg, c))
    # end
end
plt_std

plot(plt_truth, plt_post, size=(800, 500))

B = [sum(area(leaf.boundary) * sample(leaf.data.state_particles)
    for leaf in allleaves(tree)) for i in 1:1000]

histogram(B)

dA = 1.0
vline!([sum(Z*dA)])

ii = sortperm([leaf.boundary.origin for leaf in allleaves(tree)])
heatmap(reshape(μs[ii], 32, 32))


text:
#* # Advent of Code 2021 - Day 24 (solved on 16 jan 2022)

#* ## Introdution

#* This puzzle is one of those where a direct, naive implementation would take a long time to run. In principle, we need to check combinations of 14 digits from 1 to 9, i.e. 9^14 combinations. Even if each combination would take just 1 nanosecond (10^-9 seconds), this would amount, in the worst case, to about 22877 seconds, or 6 hours.

#* The fastest I could get for each check, by translating the set of 252 instructions to a pure julia function (actually a bit less due to some no-ops or redundancies), was about 20 nanoseconds on a new Mac M1, which would lead to a total of 120 hours.

#* The actual solution for part 1 has `9` as the first digit and `1` as the second, which means I would actually "just" need to check (going down each digit from `9` to `1`) 9^13 digits before finding the solution, but that still amounts to 14 hours. The second part would take much longer, because the natural path is to start from `1`, but the actual solution turned out to also have `9` as the first digit, which would leads us to 100+ hours.

#* I could have just waited, but that would be no fun and it would require running the code on my machine for days, eating up resources, and without actually knowing it would work. So, we need some sort of dynamic programming to break this down. The problem is that things get messy quickly. But let us get to it.

#* ## The puzzle

#* We are given a set of instructions operating on four variables named `w`, `x`, `y`, `z`, and with 14 integer input data with unknown values in the range `1` to `9`.

#* The instructions are either the unary operation of input of new data to one of the variables, or one of the diadic operations `add`, `mul`, `div`, `mod`, `eql`, where the last one is to check for equality.

#* The diadic operations are always of the form `op a b`, which means `a = op(a, b)`, where `op` is one of the above diadic operations.

#* In my case, the input instructions are always of the form `inp w`, i.e.always reading to variable `w`. And the equality check appears in sequence, always `eql x w` followed by `eql x 0`. In effect, this is to check whether `x ≠ w` and store the result in `x` as `0` or `1`, i.e. the two equality boil down to a single `x = 1 * (x ≠ w)` instruction. Knowing that helps to simplify stuff.

#* Given a set of 14 input data, the sequence of operations leads to a final value for the variables `w`, `x`, `y`, and `z`. But the challenge is not to compute the result for a given input. The challenge, in part 1, is to find the largest valid input. A *valid* input is one which ends with `z = 0`. The second part is to find the smallest valid input.

#* ## The approach to solve it

#* The naive approach would be to start with the highest possible input `9999999999999`, down to the smallest one `11111111111111`, computing the final values of the variables `w`, `x`, `y`, `z`, and stopping at the first valid input (`w == 0`). As explained in the Introdution, this would take a rather long time.

#* So we need somehow to keep the input as unknowns and track the conditions that lead to `w = 0` at the end.

#* The way I thought for solving this was to use symbolic computation. The right package for that seemed at first to be `Symbolics.jl`. But this is still in its early stages. It is not complete. It does accept `mod` and `rem`, but there is no `div`. And, worst, there is no simplification implemented for `mod` and `rem`. I decided to have a look at it and see how I could add them. It turned out `Symbolics.jl` uses `SymbolicUtils.jl` under the hood, and it is there where the simplification magic occurs. It is in `SymbolicUtils.jl` that this core stuff is/should be implemented.

#* I looked at the code of `SymbolicUtils.jl` and it seemed not too difficult to implement what I needed. After checking with the developers on Slack, I went ahead and started to make this happen. This became PR [JuliaSymbolics/SymbolicUtils.jl PR: adds div and simplifications for div, mod and rem #418](https://github.com/JuliaSymbolics/SymbolicUtils.jl/pull/418).

#* This PR adds `div` and adds a bunch of simplification rules for `div`, `mod`, and `rem`, like `div(x, 1) = x`, `mod(x, 1) = 0`, `div(x, -1) = -x`, `mod(x, -1) = 0`, `div(n * x, m) = div(n, m) * x + div(x - div(n, m) * x, m)`, and so one.

#* This, however, is still not enough to crack the puzzle. The symbolic expressions stored in the variables still become gigantic and unmannable. We need more simplifications.

#* The last piece of symbolics, then, is to take advantage of the fact that the input is between `1` and `9` and that, at least in my case, the `div` and `mod` operations come either dividing by `1` (which the implemented rules above take care) or by `26`. Hence, something like `div(n + i, m)` is zero when `i` is an unkown input data between `1` and `9`, `0 ≤ rem(n, m) ≤ m - 9`, and `m = 26`. Similarly, `mod(i + n, m)` reduces to `mod(n, m) + i`.

#* We, therefore, implement these rules, but these are implemented here in this script, not in SymbolicUtils.jl, because this is specific to this problem, not general facts for arbitrary symbolic variables.

#* Finally, equalities like `i1 == m`, `n + i1 == m`, and `n + i1 == i2` would certainly vanish if, respectively, `m ≥ 10`, `m - n ≥ 10`, `n ≥ 10`. We implement these rules here, as well.

#* With that at hand, for each given operation, we apply the simplifications to reduce the expressions, and go on.

#* The dynamic programming part comes when one of the equality operations yield two possible outcomes (1 or 0, i.e. true or false), depending on the input data. When this happens, we branch the code, one for each possibility, keeping track of the conditions taken.

#* At the end, we check whether `w = 0` or not. If not, we return `nothing` and the branch is discarded. If `w = 0`, the set of conditions for that corresponding branch is returned.

#* It could happen that more than one set of conditions are valid, but it turns out there is only one valid combination. This comes out as a set of linear equations involving the 14 input data. With this set of conditions at hand, we look for the highest (for part 1) or lowest (for part 2) combination satisfying the conditions. These are the solutions.

#* ## The implementation

#* ### Adding SymbolicUtils
#*
#* As mentioned above, we need `SymbolicUtils.jl`, but we need the version in the PR, which, at the moment, has not been merged yet. So, if you want to attempt this, be sure you add the PR version to your environment, or rather the branch `divmodrem` from my fork of `SymbolicUtils.jl`.

#* You can do that with
#* ```julia
#* julia> using Pkg
#*
#* julia> Pkg.add(url="https://github.com/rmsrosa/SymbolicUtils.jl", rev="divmodrem")
#* ```

#* ### Importing SymbolicUtils
#*
#* With the package added, we load it and import some non-exported functions which will be needed to add the rules that take advantage of the `1-9` range of the input data.

using SymbolicUtils
import SymbolicUtils: Chain, Sym, is_literal_number, _isnonzerointeger, ifelse

#* ### Reading the commands
#*
#* The puzzle input is different for each user. You can download it once you log in to the AoC site. My input has been downloaded to a local file. We can load it as follows

list = readlines(@__DIR__() * "/day24_input.txt")

#* ### The commands in julia
#*
#* The diadic commands are `add`, `mul`, `div`, `mod`, and `eql`.
#*
#* They correspond to `+`, `*`, `div`, `mod`, and `==` in julia, so we just
#* write a dictionary associating their (string) name to their corresponding (julia) operator. (It would be better to make them constant, for performance issues, but this is not a critical thing here, so we don't bother).

OP = Dict(
    "add" => +,
    "mul" => *,
    "div" => div,
    "mod" => mod,
    "eql" => ==
)

#* As mentioned above, since the input is between 1 and 9, we add extra rules to simplify `div(i + n, m)` and `mod(i + n, m)` when `i` is an input and `n` and `m` are integers with `n + 9 < m`. And also to evaluate `i1 + n == i2` to false, when `n ≥ 9`, and the like. This is done as follows. (Check out [JuliaSymbolics/SymbolicUtils.jl](https://symbolicutils.juliasymbolics.org) in case you are interested in learning how these work.)

r0 = Chain([
    @rule div(~i::Sym, ~m::_isnonzerointeger) => ifelse(~m > 9, 0, nothing)

    @rule div(~n::is_literal_number + ~i, ~m::_isnonzerointeger) => ifelse(0 ≤ mod(~n, ~m) < ~m - 9, div(~n, ~m), nothing)

    @rule div(~n::is_literal_number + ~k::is_literal_number * ~i, ~m::_isnonzerointeger) => ifelse(0 ≤ mod(~n, ~m) < ~m - 9, div(~n, ~m) + div(mod(~n, ~m) + ~k * ~i, ~m), nothing)

    @rule div(~n::is_literal_number + ~e::Sym + ~m::_isnonzerointeger * ~i, ~m) => ~i2 + div(~n + ~e, ~m)

    @rule mod(~i::Sym, ~m::_isnonzerointeger) => ifelse(~m > 9, ~i, nothing)

    @rule mod(~n::is_literal_number + ~i, ~m::_isnonzerointeger) => ifelse(0 ≤ mod(~n, ~m) < ~m - 9, mod(~n, ~m) + ~i, nothing)

    @acrule ~i::Sym == ~n::is_literal_number => ifelse(~n > 9 || ~n < 1, 0, nothing)

    @acrule ~i::Sym == ~n::is_literal_number + ~j::Sym => ifelse(~n ≥ 9 || ~n ≤ -9, 0, nothing)
])

r = Chain([
    @rule ~x => r0(~x)
    @rule +(~~xs) => sum(r0.(~~xs))
])

srs = r ∘ simplify

#* Perusing the commands, we see, by the fourth of fifth input, that there will be a branch depending on the data. More precisely, the value of `x` will be either `0` or `i`, depending on the difference between two previous input and a given integer in the instructions. This will have further branches later on. Hence, we will need some form of dynamic programming structure to keep track of that.

#* We define a vector of Symbolic integers. In the challenge, there will be 14 integer inputs, but it makes no effect writing i[1:14] (this works in `Symbolics.jl`, but not in `SymbolicUtils.jl`).

#* All four variables start with 0.

#* and the input starts with 1, but we initialize it to zero and add 1 when a new input is needed

#* Now the idea is to loop through the whole `list` of instructions and execute the commands, branching when needed. For that, we keep track of the line number we are in the `list` of commands, starting at `0` and increasing it by one before each command is read from the `list`.

function state(conditions, u, list, i, j, k)
    while j < length(list)
        j += 1
        mc = match(r"([a-z]{3}) ([w-z])(?: ([w-z]|-?\d*))?", list[j]).captures
        # Int('w') is 119, so by subtracting 118 we get
        # w => 1, x => 2, y => 3, z => 4
        idx = Int(only(mc[2])) .- 118
        if mc[1] == "inp"
            k += 1
            u[idx] = i[k]
        else
            aux = mc[3] ∈ ("w", "x", "y", "z") ? u[Int(only(mc[3])) - 118] : parse(Int, mc[3])
            operator = OP[mc[1]]
            u[idx] = srs(1*operator(u[idx], aux))
        end
        # In my case, at least, conditionals always appear as the
        # two consecutive commands `eql x w`, `eql x 0`.
        # This can be thought of a unique `x = x ≠ w` instruction.
        # When `x` and `w` are numeric, or when it is impossible for
        # them to be equal, since the input is restricted to the range
        # `1-9`, the direct calculation or the symbolic reductions will
        # resolve the value of `x` to either `0` or `1`. Otherwise, `x`
        # will be a symbolic condition depending on the data.
        # We check when this is the case and bifurcate the code
        # depending on whether the condition is met or not.
        # When doing so, we record, in the `conditions` vector, the
        # corresponding condition chosen for each branch.
        # We check that by checking when the command is `eql x w` and when `x == w` is not a literal number. This dynamic programming returns
        # the set of conditions needed to have (or possibly have) `w == 0`
        # at the end. When `w ≠ 0` at the end, the branching returns
        # `nothing` and the set of conditions on that route is discarded.

        if list[j] == "eql x w" && !is_literal_number(u[idx])
            utrue = copy(u)
            utrue[idx] = 1
            ufalse = copy(u)
            ufalse[idx] = 0
            res_true = state([conditions; u[idx]], utrue, list, i, j, k)
            res_false = state([conditions; replace("$(u[idx])", "==" => "≠")], ufalse, list, i, j, k)
            if res_true === res_false === nothing
                return nothing
            elseif res_true === nothing
                return res_false
            elseif res_false === nothing
                return res_true
            else
                return [res_true; res_false]
            end
        end
    end
    # Well, we could have ended with `w` as a combination of inputs
    # that could vanish depending on the choice of the input, but this
    # actually doesn't happen (they are all - except one - a linear
    # combination of the input data with positive coefficients) and
    # it suffices to check whether `w == 0`.
    return is_literal_number(u[end]) && u[end] == 0 ? conditions : nothing
end

function solve(list)
    i = @syms i1::Integer i2::Integer i3::Integer i4::Integer i5::Integer i6::Integer i7::Integer i8::Integer i9::Integer i10::Integer i11::Integer i12::Integer i13::Integer i14::Integer

    u = Vector{Any}([0, 0, 0, 0]) # slots w, x, y, z respectively
    j = 0
    k = 0
    conditions = state([], u, list, i, j, k)
    return conditions
end

#* Finally, we call the main loop
conditions = solve(list)

#=
Result:
 (i3 - 7) == i4
 (1 + i5) == i6
 (5 + i8) == i9
 i10 == i11
 (8 + i7) == i12
 (7 + i2) == i13
 (i1 - 8) == i14

The largest possible solution of this system gives us the answer to part 1:
92928914999991

The smallest possible solution of this system gives us the answer to part 2:
91811211611981

=#

max_stars_repo_name: rmsrosa/adventofcode

max_issues_repo_name: rmsrosa/adventofcode

max_forks_repo_name: rmsrosa/adventofcode


text:
module AGH
using LinearAlgebra, SparseArrays, ProgressMeter, Statistics, LinearAlgebra,
    Random, Plots, Octavian

include("A.jl")
include("G.jl")
include("H.jl")
include("file-io.jl")
include("one-step.jl")
include("picd/PICD.jl")

end # module


text:
import Base.isempty
export enum, enumbreadthfirst, enumdepthfirst, enumlikelyfirst

type Stack
    store
    Stack() = new(Any[])
end

# I doubt enq!/deq! will run in constant time for this. Perhaps look
# at DataStructures.jl.
type Queue
    store
    Queue() = new(Any[])
end

# Provide a minimal common interface over Stack, Queue, PriorityQueue.

enq!(a::Union(Stack,Queue), item) = push!(a.store, item)

deq!(a::Stack) = pop!(a.store)
deq!(a::Queue) = shift!(a.store)

isempty(a::Union(Stack,Queue)) = isempty(a.store)

enq!(q::PriorityQueue, item) = enqueue!(q, item, -item[1])
deq!(q::PriorityQueue) = dequeue!(q)

type Enum <: Ctx
    score::Float64
    path
    unexplored
    Enum(queuetype) = new(0, Any[], queuetype())
end

# @pp
function sample(s::Store, k::Function, a, e::ERP, ctx::Enum)
    for val in support(e)
        enq!(ctx.unexplored, (ctx.score + score(e, val), () -> k(s,val), [ctx.path, val]))
    end
    runnext(ctx)
end

# @pp
function factor(s::Store, k::Function, a, score, ctx::Enum)
    ctx.score += score
    k(s, nothing)
end

function runnext(ctx::Ctx)
    score, cont, path = deq!(ctx.unexplored)
    ctx.score = score
    ctx.path = path
    cont()
end

# @pp
enum(s::Store, k::Function, a, comp::Function, maximumexec::Int64=0) = enumbreadthfirst(s, k, a, comp, maximumexec)

enumdepthfirst(s::Store, k::Function, a, comp::Function, maximumexec::Int64=0) = enum(s, k, a, comp, Stack, maximumexec)
enumbreadthfirst(s::Store, k::Function, a, comp::Function, maximumexec::Int64=0) = enum(s, k, a, comp, Queue, maximumexec)
enumlikelyfirst(s::Store, k::Function, a, comp::Function, maximumexec::Int64=0) = enum(s, k, a, comp, PriorityQueue, maximumexec)

function enum(store::Store, k::Function, address, comp::Function, queuetype::DataType, maximumexec::Int64=0)
    global ctx
    returns = Dict{Any,Float64}()
    ctxold, ctx = ctx, Enum(queuetype)
    try
        currentexec = 0
        trampoline() do
            partial(comp,store)(address) do _store, value
                currentexec += 1
                returns[value] = get(returns, value, 0) + exp(ctx.score)
                if maximumexec > 0 && currentexec == maximumexec
                    info("Maximum executions reached.")
                elseif !isempty(ctx.unexplored)
                    runnext(ctx)
                end
            end
        end
    finally
        ctx = ctxold
    end
    # Note: The context must be restored before calling k.
    k(store, Discrete(returns))
end