text:
export SparseMatrix

import Base.resize!, Base.insert!, Base.getindex, Base.setindex!

export from_csr, from_graph, resize!, sort!
export equals, contains, diag, transpose
export get_row,
    get_column, get_rowwise, set_rowwise!, get_columnwise, set_columnwise!

@enum Layout default_layout lower_triangle upper_triangle lower_symmetric upper_symmetric

mutable struct SparseMatrix{T}
    columns::Vector{Int64}
    columns_rows::Vector{Int64}
    values::Vector{T}

    rows::Vector{Int64}
    rows_columns::Vector{Int64}
    positions::Vector{Int64}

    layout::Layout
    nnz::Int64
    row_count::Int64
    column_count::Int64

    list::Intlist

    function SparseMatrix{T}(
        row_count::Int64,
        column_count::Int64,
        capacity::Int64,
    ) where {T}
        local sm = new()

        resize!(sm, row_count, column_count, capacity)

        return sm
    end
end

function resize!(
    sm::SparseMatrix{T},
    row_count::Int64,
    column_count::Int64,
    capacity::Int64,
) where {T}
    sm.list = Intlist(max(row_count, column_count))
    sm.layout = default_layout

    sm.nnz = 0
    sm.row_count = row_count
    sm.column_count = column_count

    sm.columns = zeros(Int64, 1 + column_count)
    sm.columns_rows = zeros(Int64, capacity)
    sm.values = zeros(T, capacity)

    sm.rows = zeros(Int64, 1 + row_count)
    sm.rows_columns = zeros(Int64, capacity)
    sm.positions = zeros(Int64, capacity)
end

function from_csr(
    row_count::Int64,
    column_count::Int64,
    rows::Vector{Int64},
    columns::Vector{Int64},
    values::Vector{T},
) where {T}
    @assert row_count == (length(rows) - 1)
    @assert length(values) == length(columns)
    @assert length(values) == rows[end] - 1

    local sm = SparseMatrix{T}(row_count, column_count, length(values))
    sm.nnz = length(values)
    sm.rows = copy(rows)

    for i = 1:row_count
        for j = rows[i]:(rows[i+1]-1)
            sm.columns[columns[j]+1] += 1
        end
    end

    sm.columns[1] = 1
    for j = 2:(1+column_count)
        sm.columns[j] += sm.columns[j-1]
    end

    local cc = copy(sm.columns)
    local rr = copy(sm.rows)
    for i = 1:row_count
        for j = rows[i]:(rows[i+1]-1)
            local c = cc[columns[j]]
            sm.columns_rows[c] = i
            sm.values[c] = values[j]

            sm.rows_columns[rr[i]] = columns[j]
            sm.positions[rr[i]] = c

            cc[columns[j]] += 1
            rr[i] += 1
        end
    end

    sort!(sm)

    return sm
end

function from_graph(
    g::DirectedGraph{T},
    row_count::Int64,
    column_count::Int64,
) where {T}
    local sm = SparseMatrix{T}(row_count, column_count, g.size)
    sm.nnz = g.size

    sm.rows[1] = 1

    local n = 1
    local values = zeros(T, sm.nnz)
    for i = 1:row_count
        local j = g.start[i]
        while NIL != j
            values[n] = g.edges[j]
            sm.rows_columns[n] = g.vertices[j]
            sm.columns[g.vertices[j]] += 1
            j = g.next[j]
            n += 1
        end
        sm.rows[i+1] = n
    end

    n = 1
    for i = 1:(column_count+1)
        local t = n
        n += sm.columns[i]
        sm.columns[i] = t
    end

    local col_index = copy(sm.columns)
    n = 1
    for i = 1:row_count
        local j = sm.rows[i]
        for j = sm.rows[i]:(sm.rows[i+1]-1)
            local col = sm.rows_columns[j]
            sm.positions[j] = col_index[col]
            sm.values[col_index[col]] = values[j]
            sm.columns_rows[col_index[col]] = i
            col_index[col] += 1
        end
    end

    sort!(sm)

    return sm
end

function equals(sm::SparseMatrix{T}, tm::SparseMatrix{T}) where {T}
    if sm.nnz != tm.nnz
        return false
    end
    if sm.row_count != tm.row_count
        return false
    end
    if sm.column_count != tm.column_count
        return false
    end
    if sm.columns != tm.columns
        return false
    end
    if sm.columns_rows != tm.columns_rows
        return false
    end
    if sm.positions != tm.positions
        return false
    end
    if sm.values != tm.values
        return false
    end

    return true
end

function Base.:transpose(sm::SparseMatrix{T}) where {T}
    local tm = SparseMatrix{T}(sm.column_count, sm.row_count, sm.nnz)

    tm.nnz = sm.nnz
    tm.columns[:] = sm.rows[:]
    tm.columns_rows[:] = sm.rows_columns[:]
    tm.rows[:] = sm.columns[:]
    tm.rows_columns[:] = sm.columns_rows[:]

    local k = 1
    for j = 1:sm.column_count
        for i = sm.columns[j]:(sm.columns[j+1]-1)
            tm.positions[sm.positions[i]] = k
            tm.values[k] = sm.values[sm.positions[i]]
            k = k + 1
        end
    end

    return tm
end

function contains(sm::SparseMatrix{T}, row::Int64, column::Int64) where {T}
    @assert column > 0 && column <= sm.column_count
    @assert row > 0 && row <= sm.row_count

    for j = sm.columns[column]:(sm.columns[column+1]-1)
        if row == sm.columns_rows[j]
            return true
        end
    end
    return false
end

function diag(sm::SparseMatrix{T}) where {T}
    @assert sm.column_count == sm.row_count

    local values = Vector{T}()
    for i = 1:sm.row_count
        push!(values, get_columnwise(sm, i, i))
    end

    return values
end

function get_row(sm::SparseMatrix{T}, row::Int64) where {T}
    @assert row > 0 && row <= sm.row_count

    local values = SparseArray{T}(sm.column_count)
    for i = sm.rows[row]:(sm.rows[row+1]-1)
        values[sm.rows_columns[i]] = sm.values[sm.positions[i]]
    end

    return values
end

function get_column(sm::SparseMatrix{T}, column::Int64) where {T}
    @assert column > 0 && column <= sm.column_count

    local values = SparseArray{T}(sm.row_count)

    for j = sm.columns[column]:(sm.columns[column+1]-1)
        values[sm.columns_rows[j]] = sm.values[j]
    end

    return values
end

function get_rowwise(sm::SparseMatrix{T}, row::Int64, column::Int64) where {T}
    @assert column > 0 && column <= sm.column_count
    @assert row > 0 && row <= sm.row_count

    for i = sm.rows[row]:(sm.rows[row+1]-1)
        if column == sm.rows_columns[i]
            return sm.values[sm.positions[i]]
        end
    end

    return missing
end

function set_rowwise!(
    sm::SparseMatrix{T},
    row::Int64,
    column::Int64,
    value::T,
) where {T}
    @assert column > 0 && column <= sm.column_count
    @assert row > 0 && row <= sm.row_count

    for i = sm.rows[row]:(sm.rows[row+1]-1)
        if column == sm.rows_columns[i]
            sm.values[sm.positions[i]] = value
            return
        end
    end

    error("out of range")
end

function get_columnwise(
    sm::SparseMatrix{T},
    row::Int64,
    column::Int64,
) where {T}
    @assert column > 0 && column <= sm.column_count
    @assert row > 0 && row <= sm.row_count

    for i = sm.columns[column]:(sm.columns[column+1]-1)
        if row == sm.columns_rows[i]
            return sm.values[i]
        end
    end

    return T(0)
end

function set_columnwise!(
    sm::SparseMatrix{T},
    row::Int64,
    column::Int64,
    value::T,
) where {T}
    @assert column > 0 && column <= sm.column_count
    @assert row > 0 && row <= sm.row_count

    for i = sm.columns[column]:(sm.columns[column+1]-1)
        if row == sm.columns_rows[i]
            sm.values[i] = value
            return
        end
    end

    error("out of range")
end

function Base.:sort!(sm::SparseMatrix{T}) where {T}
    local old_indices = 1:sm.nnz
    local new_indices = collect(old_indices)

    for j = 1:sm.column_count
        local range = sm.columns[j]:sm.columns[j+1]-1
        local perm = (sm.columns[j] - 1) .+ sortperm(sm.columns_rows[range])
        new_indices[range] = new_indices[perm]
    end

    sm.columns_rows[old_indices] = sm.columns_rows[new_indices]
    sm.values[old_indices] = sm.values[new_indices]

    for j = 1:sm.row_count
        local range = sm.rows[j]:sm.rows[j+1]-1
        local perm = (sm.rows[j] - 1) .+ sortperm(sm.rows_columns[range])
        sm.rows_columns[range] = sm.rows_columns[perm]
        sm.positions[range] = sm.positions[perm]
    end

    for i = 1:sm.nnz
        sm.positions[i] = new_indices[sm.positions[i]]
    end
end


text:
@testset "Chain" begin
    a = Chain(1, Int, Int)
    @test dim(a) == 1
    @test iszero(a)
    @test eltype(a) == (Int, Int)

    b = Chain(Int[],Int[])
    c = Chain(Int)

    push!(c, (1,0))
    @test c[1] == (1,0)

    push!(b, 1=>1)
    @test b[1] == (1,1)

    a+=c
    a+=b
    a+=(1,2)
    @test a[1] == (1,0)
    @test a[2] == (1,1)
    @test a[3] == (1,2)

    b+=(1,2)
    @test b[2] == (1,2)

    @testset "loop thrugh chain" for (chc,chel) in a
        @test chc == 1
    end

    append!(c,b)
    @testset "two chain compare" for (aa,cc) in zip(a,c)
        @test cc == aa
    end

    @testset "simplify" for (chc,chel) in simplify(a+c)
        @test chc == 2
    end

    @test (2*b)[1] == (2,1)
    @test (b*2)[1] == (2,1)

    @test length(simplify(a-c)) == 0
end

ext: jl


text:
using Interpolations: LinearInterpolation, Throw
using StaticArrays: SA

# define the hydrogenic continuum absorption

"""
    ndens_state_hydrogenic(n, nsdens_div_partition, ion_energy, T)

Calculates the LTE number density (in cm^-3) of a hydrogenic species at a given energy level.

# Arguments
- `n::Integer`: The quantum number of the state
- `nsdens_div_partition`: The total (including all states) number density of the current
ionization species (in cm^-3) divided by the species's partition function.
- `ion_energy`: The minimum energy (in eV) required to ionize the species (from the ground
state). This can be estimated as Z²*ion_energy of hydrogen or Z²*Rydberg_H.
- `T`: Temperature

# Note
This assumes that all hydrogenic species have a statistical weight of gₙ = 2*n².

This was taken from equation (5.4) of Kurucz (1970) (although this comes directly from the
boltzmann equation).
"""
function ndens_state_hydrogenic(n::Integer, nsdens_div_partition::Real, T::Real, ion_energy::Real)
    n2 = n*n
    g_n = 2.0*n2
    energy_level = ion_energy - ion_energy/n2
    nsdens_div_partition * g_n * exp(-energy_level/(kboltz_eV *T))
end

_eVtoHz(energy) = energy/hplanck_eV
_HztoeV(freq) = freq*hplanck_eV

# In the following equation, e is the elementary charge in units of statcoulombs
const _bf_σ_const = 2.815e29 # = 64*π⁴e^10 mₑ / (c h⁶ 3√3)

const _bf_σ_coef = SA[(0.9916,  2.719e13, -2.268e30),
                      (1.105,  -2.375e14,  4.077e28),
                      (1.101,  -9.863e13,  1.035e28),
                      (1.101,  -5.765e13,  4.593e27),
                      (1.102,  -3.909e13,  2.371e27),
                      (1.0986, -2.704e13,  1.229e27),
                      (1.0,     0.0,       0.0     )]

# this is helper function is the most likely part of the calculation to change.
# This uses double precision to be safe about the polynomial coefficients.
function _hydrogenic_bf_cross_section(Z::Integer, n::Integer, ν::Real, ion_freq::Real)
    # this implements equation 5.5 from Kurucz (1970)
    # - Z is the atomic number
    # - n is the energy level (remember, they start from n=1)
    # - ν should have units of Hz
    # - ion_freq is the frequency of the photon carrying the minimum energy needed to ionize the
    #   ground state configuration of the current ion (in Hz). This can be estimated as
    #   _eVtoHz(Z²*RydbergH_eV).

   if (n < 1)
        throw(DomainError(n,"n must be a positive integer"))
    end

    Z2 = Z * Z
    inv_n = 1.0/n
    inv_n2 = inv_n*inv_n
    inv_ν = 1.0/ν
    if ν < (ion_freq*inv_n2)
        0.0
    else
        # the last entry of the table should be used when n > length(_bf_σ_coef)
        A, B, C = _bf_σ_coef[min(n,length(_bf_σ_coef))]

        poly = A + (B + C*Z2*inv_ν)*Z2*inv_ν
        # equation 5.5 of Kurucz had a typo in it. In the numerator of the fraction that is
        # multiplied by the entire polynomial, Z² should be Z⁴. This was discovered during
        # comparisons with data from the Opacity Project, and can be confirmed by:
        # - looking at eqn 5.6 of Kurucz (it uses Z⁴ instead of Z²)
        # - comparison against equation 10.54 of Rybicki & Lightman
        _bf_σ_const*Z2*Z2*(inv_n2*inv_n2*inv_n)*(inv_ν^3)*poly
    end
end

"""
    _hydrogenic_bf_high_n_absorption(ν, T, Z, nmin, nsdens_div_partition, ion_energy)

Approximates the sum of absorption coefficient contributions by all bound-free transitions in which
the electron originates in an energy level of nmin or larger using an integral.

In full detail, the calculated quantity is described as
    ∑_{n'=nmin}^∞ ndens(n=n') * α_ν(n=n') * (1 - exp(-hν/k/T))/ρ,
where α(n=n') is the absorption coefficient for the bound-free atomic absorption coefficient
(uncorrected for stimulated emission). This function approximates this sum with an integral.

# Arguments
- `ν`: frequency in Hz
- `T`: temperature in K
- `Z::Integer`: Z is the atomic number of the ion (1 for HI)
- `nmin::Integer`: The lowest energy level (principle quantum number) included in the calculation
- `nsdens_div_partition` is the number density of the current species divided by the
   partition function.
- `ion_energy`: the ionization energy from the ground state (in eV).

# Notes
This is based on equation (5.6) from Kurucz (1970) (which give the formula for opacity). I think ρ
was simply omitted from that equation.
"""
function _hydrogenic_bf_high_n_absorption(ν::Real, T::Real, Z::Integer, nmin::Integer,
                                          nsdens_div_partition::Real, ion_energy::Real)

    # this function corresponds to the evaluation of a integral. We subdivide the solution into two
    # parts: (i) consts and (ii) integral. The solution is the product of both parts

    # first, compute consts (terms pulled out of the integral)
    β = 1.0/(kboltz_eV*T)
    expν_term = exp(-β*hplanck_eV*ν)
    integral_consts = nsdens_div_partition * _bf_σ_const * Z^4 * (1.0 - expν_term) * 2.0 / ν^3

    # second, compute the integral
    # need to know the ionization frequency (frequency of a photon that could ionize species
    # from ground state)
    ion_freq = _eVtoHz(ion_energy)

    ionE_times_β = ion_energy*β
    coef = 0.5/ionE_times_β

    # take 2 different approaches based on the photon's energy
    if ν >= (ion_freq/nmin^2)
        # At this frequency, photons are energetic enough to ionize electrons at n >= nmin
        # integral of exp(-ion_energy*(1-(1/n²)) /(kboltz*T))/n³ from nmin to ∞
        integral = coef * (exp(-(ionE_times_β - ionE_times_β/nmin^2)) - exp(-ionE_times_β))
    else
        # At this frequency, photon aren't energetic enough to ionize electrons at:
        # nmin ≤ n < ceil(sqrt(ion_freq/ν)). The absorption coefficients for these electrons are
        # all zero.
        # Thus, the integral's lower limit should actually be ceil(sqrt(ion_freq/ν)). Kurucz
        # approximates this by replacing ion_energy*β/nmin² in the prior solution with hplanck*ν*β
        integral = coef * (exp(-(ionE_times_β - hplanck_eV*ν*β)) - exp(-ionE_times_β))
    end

    (integral_consts*integral)
end

"""
    hydrogenic_bf_absorption(ν, T, Z, nmax_explicit_sum, nsdens_div_partition,
                             ion_energy, [integrate_high_n])

Compute the bound-free linear absorption coefficient contributed by all energy states of a
Hydrogenic species

The calculation is broken into 2 parts: (i) the contributions of the lowest energy states are
explicitly summed and (ii) the contributions of the higher energy states are estimated with an
integral.

# Arguments
- `ν`: frequency in Hz
- `T`: temperature in K
- `Z::Integer`: Z is the atomic number of the species (e.g. 1 for H I or 2 for He II)
- `nsdens_div_partition` is the total number density of the species divided by the species's
   partition function.
- `ion_energy`: the ionization energy from the ground state (in eV). This can be
   estimated as Z²*Rydberg_H (Rydberg_H is the ionization energy of Hydrogen)
- `nmax_explicit_sum::Integer`: The highest energy level whose absorption contribution is included
   in the explicit sum. The contributions from higher levels are included in the integral.
- `integrate_high_n::bool`: When this is `false`, bf absorption from higher energy states are not
   estimated at all. Default is `true`.

# Notes
This follows the approach described in section 5.1 of Kurucz (1970).
"""
function hydrogenic_bf_absorption(ν::Real, T::Real, Z::Integer, nsdens_div_partition::Real,
                                  ion_energy::Real, nmax_explicit_sum::Integer,
                                  integrate_high_n::Bool=true)
    ionization_freq = _eVtoHz(ion_energy)

    # first, directly sum the individual absorption contributions from H I atoms at each of the
    # lowest energy levels (i.e. all energy levels where n <= nmax_explicit_sum)
    partial_sum = 0.0
    for n = 1 : nmax_explicit_sum
        ndens_state = ndens_state_hydrogenic(n, nsdens_div_partition, T, ion_energy)
        hydrogenic_bf_cross_section = _hydrogenic_bf_cross_section(Z, n, ν, ionization_freq)
        cur_val = ndens_state * hydrogenic_bf_cross_section
        partial_sum += ndens_state * hydrogenic_bf_cross_section
    end
    α_low_n = partial_sum * (1.0 - exp(-hplanck_eV * ν / (kboltz_eV * T)))

    # second, estimate the absorption contributions from H I atoms at higher energy levels using an
    # integral approximation (assuming integrate_high_n is true)
    α_high_n = _hydrogenic_bf_high_n_absorption(ν, T, Z, nmax_explicit_sum+1, nsdens_div_partition,
                                                ion_energy)

    α_low_n + (α_high_n * integrate_high_n)
end

"""
    _load_gauntff_table([fname])

Returns a table of thermally-averaged free-free Gaunt factors, and the values of log₁₀(γ²) and
log₁₀(u) associated with each point.

This loads the non-relativistic free-free data published by
[van Hoof et al. (2014)](https://ui.adsabs.harvard.edu/abs/2014MNRAS.444..420V/abstract).

Note: This function code could trivially be adapted to load the relativistic free-free gaunt
factors published by
[van Hoof et al (2015)](https://ui.adsabs.harvard.edu/abs/2015MNRAS.449.2112V/abstract).
"""
function _load_gauntff_table(fname = joinpath(_data_dir, "vanHoof2014-nr-gauntff.dat"))
    parse_header_line(type, s) = (last = something(findfirst('#',s), length(s)+1) - 1;
                                  parse.(type,split(s[1:last])) )

    f = open(fname, "r")
    skipchars((c)->false, f; linecomment='#') # skip 1st section of comments

    log10_γ2, log10_u = begin # parse header values
        magic_number = parse_header_line(Int64, readline(f))[1]
        @assert magic_number == 20140210
        num_γ2, num_u = parse_header_line(Int64, readline(f))
        log10_γ2_start = parse_header_line(Float64, readline(f))[1]
        log10_u_start = parse_header_line(Float64, readline(f))[1]
        # step_size has units of dex and applies to both axes
        step_size = parse_header_line(Float64, readline(f))[1]

        (range(log10_γ2_start, length = num_γ2, step = step_size),
         range(log10_u_start, length = num_u, step = step_size))
    end
    skipchars((c)->false, f; linecomment='#') # skip 2nd section of comments

    out = Matrix{Float64}(undef, length(log10_u), length(log10_γ2))
    for i in 1:length(log10_u)
        view(out,i,:) .= parse.(Float64, split(readline(f)))
    end

    # at this point there is a number section of comments followed by a grid of uncertainties for
    # each value. We'll ignore those. That secondary section may or may not be present in versions
    # of the table that included relativistic calculations
    close(f)

    out, log10_γ2, log10_u
end

# maybe also initialize λ_bounds and T_bounds
const _gauntff_interpolator, _gauntff_T_bounds, _gauntff_λ_bounds = begin

    # load in the tabulated data
    table_val, log10_γ2, log10_u = _load_gauntff_table()

    # the table is bigger than it needs to be. We can easily cut its size by a factor of ~10
    T_extrema = [100.0, 1e6] # units of K
    λ_extrema = [1.0e-6, 1.0e-2] # units of cm (equivalent to 100 Å and 100 μm)
    Z_extrema = [1, 2]

    function _find_bound_inds(range, min, max) # this is very similar to some duplicated code
        lb,ub = searchsortedlast(range, min), searchsortedfirst(range, max)
        @assert (range[lb] ≤ min) && (range[ub] ≥ max)
        lb, ub
    end
    calc_log10_γ2(Z, Tₑ) = log10(Rydberg_eV * Z^2 / (kboltz_eV * Tₑ))
    calc_log10_u(λ, Tₑ) = log10(hplanck_cgs * c_cgs / (λ * kboltz_cgs * Tₑ))

    γ2_lb, γ2_ub = _find_bound_inds(log10_γ2, extrema(calc_log10_γ2.(Z_extrema, T_extrema'))...)
    u_lb, u_ub   = _find_bound_inds(log10_u,  extrema( calc_log10_u.(λ_extrema, T_extrema'))...)

    # make a copy of the selected data so that the unneeded data can be garbage collected
    revised_table = copy(table_val[u_lb:u_ub, γ2_lb:γ2_ub])

    (LinearInterpolation((log10_u[u_lb:u_ub], log10_γ2[γ2_lb:γ2_ub]), revised_table,
                         extrapolation_bc = Throw()),
     closed_interval(T_extrema...),
     closed_interval(λ_extrema...))
end

"""
    gaunt_ff_vanHoof(log_u, log_γ2)

computes the thermally averaged, non-relativistic free-free gaunt factor by interpolating the table
provided by [van Hoof et al. (2014)](https://ui.adsabs.harvard.edu/abs/2014MNRAS.444..420V).

# Arguments
- `log_u`: Equal to log₁₀(u) = log₁₀(h*ν/(k*Tₑ))
- `log_γ2`: Equal to log₁₀(γ²) = log₁₀(Rydberg*Z²/(k*Tₑ))

`Rydberg` is the "infinite mass unit of energy" and `Tₑ` is the temperature of free electrons (for
our purposes, we assume that free electrons are in thermal equilibrium with ions and neutral
species).

# Notes

van Hoof et al. (2014) computed the associated data table with a non-relativistic approach, which
is invalid at very high temperatures. They conclude (from comparisons with a different paper) that
their "results should be accurate up to electron temperatures of roughly 100 MK". This is more than
adequate for stellar atmospheres.
In [van Hoof et al. (2015)](https://ui.adsabs.harvard.edu/abs/2015MNRAS.449.2112V), they find that
relativistic effects introduce a ∼0.75% at 100MK, for Z = 1 (when Z > 1, the change is smaller).

This function currently uses linear interpolation. However, van Hoof et al. (2014) provides an
implementation of a third-order Lagrange scheme, which "reaches a relative precision better than
1.5e-4 everywhere." The C and Fortran implementations of this scheme can be found
[here](http://data.nublado.org/gauntff/), and are copyrighted by a BSD-style license.

Earlier variants of this function used less-accurate data from section 5.1 of Kurucz (1970) that
extended over a smaller interval of data. That table was originally derived from a figure in Karsas
and Latter (1961) and it's now used for testing purposes.
"""
gaunt_ff_vanHoof(log_u, log_γ2) = _gauntff_interpolator(log_u, log_γ2)

"""
    hydrogenic_ff_absorption(ν, T, Z, ni, ne)

computes the free-free linear absorption coefficient for a hydrogenic species

The naming convention for free-free absorption is counter-intuitive. A free-free interaction is
named as though the species interacting with the free electron had one more bound electron (in
other words it's named as though the free-electron and ion were bound together). In practice, this
means that `ni` should refer to:
- the number density of H II if computing the H I free-free absorption
- the number density of He III if computing the He II free-free absorption
- the number density of Li IV if computing the Li III free-free absorption

# Arguments
- `Z::Integer`: the charge of the ion. For example, this is 1 for ionized H.
- `ni`: the number density of the ion species in cm⁻³.
- `ne`: the number density of free electrons.
- `ν`: frequency in Hz
- `T`: temperature in K

# Note
This approach was adopted from equation 5.8 from section 5.1 of Kurucz (1970). Comparison against
equation 5.18b of Rybicki & Lightman (2004), reveals that the equation in Kurucz (1970) omits the
dependence on ρ. According to Rybicki & Lightman (2004) the free-free absorption coefficient
(corrected for stimulated emission) is:
```
    α = coef * Z² *ne * ni * (1 - exp(-hplanck*ν/(kboltz*T))) * g_ff / (sqrt(T) * ν³)
```
Note that the g_ff is the free-free gaunt factor and coef is ∼3.7e8 (a more exact coefficient can
be computed from eqn 5.18a).

With this in mind, equation 5.8 of Kurucz (1970) should actually read
```
    κ = ne * n(H II) * F_ν(T) * (1 - exp(-hplanck*ν/(kboltz*T))) / ρ
```
where F_ν(T) = coef * Z² * g_ff / (sqrt(T) * ν³).

See `gaunt_ff_vanHoof` for details about where our gaunt factor data comes from. For simplicity, we
enforce temperature and ν bounds constraints that don't include all of the available gaunt factor
data. In practice, this should never be a concern for stellar spectroscopy.
"""
function hydrogenic_ff_absorption(ν::Real, T::Real, Z::Integer, ni::Real, ne::Real)
    inv_T = 1.0/T
    Z2 = Z*Z

    hν_div_kT = (hplanck_eV/kboltz_eV) * ν * inv_T
    log_u = log10(hν_div_kT)
    log_γ2 = log10((Rydberg_eV/kboltz_eV) * Z2 * inv_T)
    gaunt_ff = gaunt_ff_vanHoof(log_u, log_γ2)

    F_ν = 3.6919e8*gaunt_ff*Z2*sqrt(inv_T)/(ν*ν*ν)

    ni*ne*F_ν*(1-exp(-hν_div_kT))
end


text:
export soundspeed, halfwavelength, wavenumber
export flutelength, toneholelength, closedholecorrection

"""
  𝑐 = soundspeed(ϑ=25.0)

Calculate the speed of sound in air of temperature ϑ in m/s
"""
function soundspeed(ϑ=25.0)
  𝛾 = 1.400            # heat capacity ratio of air
  𝑅 = 8.31446261815324 # molar gas constant (J/mol/K)
  𝑀 = 0.028965369      # mean molar mass of air (kg/mol)
  𝑐 = √(𝛾*𝑅/𝑀 *273.15)*√(1.0 + ϑ/273.15)
  return 𝑐
end

"""
  𝜔 = angularfreq(𝑓=440.0)

calculate angular frequency of 𝑓 in rad/s
"""
function angularfreq(𝑓=440.0)
  𝜔 = 2π*𝑓
  return 𝜔
end

"""
  𝑘 = wavenumber(𝑓=440.0; 𝜗=25.0)

calculate wavenumber of given frequency 𝑓 in air of temperature 𝜗 in m⁻¹
"""
function wavenumber(𝑓=440.0; ϑ=25.0)
  𝑐 = soundspeed(ϑ)
  𝜔 = angularfreq(𝑓)
  𝑘 = 𝜔/𝑐
  return 𝑘
end

"""
  𝜆 = halfwavelength(𝑓=440.0; ϑ=25.0)

calculate half-wavelength of given frequency 𝑓 in air of temperature ϑ
"""
function halfwavelength(𝑓=440.0; ϑ=25.0)
  𝑐 = soundspeed(ϑ)
  𝜆 = 𝑐/2𝑓
  return 1000𝜆 # (to mm)
end

"""
  ℓₜ = flutelength(𝑓=440.0; ϑ=25.0, ⌀=19.0, ℎ=3.5, 𝛥ℓₑ=57.0, 𝛥ℓᵥ=0.0)

Calculate flute length from embouchure-hole to open-end
  for fundamental frequency 𝑓, temperature ϑ, open-end bore diameter ⌀,
  wall thickness ℎ, embouchure correction 𝛥ℓₑ and closed-hole correction 𝛥ℓᵥ
"""
function flutelength(𝑓=440.0; ϑ=25.0, ⌀=19.0, ℎ=3.5, 𝛥ℓₑ=57.0, 𝛥ℓᵥ=0.0)
  𝜆 = halfwavelength(𝑓; ϑ=ϑ)
  ⌀₊ = ⌀ + 2ℎ
  ϖ = ⌀/⌀₊
  𝜍 = 0.6133⌀/2
  𝜚 = 0.8216⌀/2
  𝛥ℓₜ = 𝜚 + ϖ * (𝜍-𝜚) + 0.057ϖ * (1-ϖ ^5)
  ℓₜ = 𝜆 - 𝛥ℓₑ - 𝛥ℓᵥ - 𝛥ℓₜ
  return ℓₜ
end

"""
  ℓₕ = toneholelength(𝑓=440.0, 𝑓ₜ=415.305, 𝑑=9.0; ϑ=25.0, ⌀=19.0, ℎ=3.5, 𝛥ℓₑ=57.0, 𝛥ℓᵥ=0.0)

Calculate distance from embouchure hole center to tone hole center
  for open frequency 𝑓, closed frequency 𝑓ₜ, tone-hole diameter 𝑑,
  temperature ϑ, tone-hole bore diameter ⌀, tone-hole height ℎ,
  embouchure correction 𝛥ℓₑ and closed-hole correction 𝛥ℓᵥ
"""
function toneholelength(𝑓=440.0, 𝑓ₜ=415.305, 𝑑=9.0; ϑ=25.0, ⌀=19.0, ℎ=3.5, 𝛥ℓₑ=57.0, 𝛥ℓᵥ=0.0)
  𝜆 = halfwavelength(𝑓; ϑ=ϑ)
  𝑔 = 𝑓/𝑓ₜ - 1
  𝜙 = (ℎ+𝑑) * (⌀/𝑑)^2 - 0.45⌀
  𝑧 = 𝑔/2 * √(1 + 4𝜙/(𝑔*𝜆)) - 𝑔/2
  𝛥ℓₕ = 𝑧 * 𝜆
  ℓₕ = 𝜆 - 𝛥ℓₑ - 𝛥ℓᵥ - 𝛥ℓₕ
  return ℓₕ
end

"""
  𝛥𝜆ₕ = closedholecorrection(𝑓=440.0, 𝑓ₜ=415.305, 𝑑=9.0; ϑ=25.0, ⌀=19.0, ℎ=3.5, 𝛥ℓₑ=57.0, 𝛥ℓᵥ=0.0)

Calculate correction due to closed hole
  for open frequency 𝑓, closed frequency 𝑓ₜ, tone-hole diameter 𝑑,
  temperature ϑ, tone-hole bore diameter ⌀, tone-hole height ℎ,
  embouchure correction 𝛥ℓₑ and closed-hole correction 𝛥ℓᵥ
"""
function closedholecorrection(𝑓=440.0, 𝑓ₜ=415.305, 𝑑=9.0; ϑ=25.0, ⌀=19.0, ℎ=3.5, 𝛥ℓₑ=57.0, 𝛥ℓᵥ=0.0)
  𝑘 = wavenumber(𝑓; ϑ=ϑ)
  𝜆ₜ= halfwavelength(𝑓ₜ; ϑ=ϑ)
  ℓₕ = toneholelength(𝑓, 𝑓ₜ, 𝑑; ϑ=ϑ, ⌀=⌀, ℎ=ℎ, 𝛥ℓₑ=𝛥ℓₑ, 𝛥ℓᵥ=𝛥ℓᵥ)
  𝜙 = 𝑑^2*ℎ / ⌀^2
  ϵ = 2/π * atan(2𝑑/13ℎ)
  𝛿ₕ = 𝜆ₜ - 𝛥ℓₑ - ℓₕ
  𝛥𝜆ₕ = (sin(𝑘*𝛿ₕ)^2 - ϵ*cos(𝑘*𝛿ₕ)^2) * 𝜙
  return 𝛥𝜆ₕ
end

max_stars_repo_name: isaacwhanson/scad-flute

max_issues_repo_name: isaacwhanson/scad-flute

max_forks_repo_name: isaacwhanson/scad-flute


text:
struct MinDomainVariableSelection{TakeObjective} <: AbstractVariableSelection{TakeObjective} end

MinDomainVariableSelection(;take_objective=true) = MinDomainVariableSelection{take_objective}()

function (::MinDomainVariableSelection{false})(cpmodel::CPModel)
    selectedVar = nothing
    minSize = typemax(Int)
    for (k, x) in branchable_variables(cpmodel)
        if x !== cpmodel.objective && !isbound(x) && length(x.domain) < minSize
            selectedVar = x
            minSize = length(x.domain)
        end
    end
    if isnothing(selectedVar) && !isbound(cpmodel.objective)
        return cpmodel.objective
    end
    return selectedVar
end

function (::MinDomainVariableSelection{true})(cpmodel::CPModel)
    selectedVar = nothing
    minSize = typemax(Int)
    for (k, x) in branchable_variables(cpmodel)
        if !isbound(x) && length(x.domain) < minSize
            selectedVar = x
            minSize = length(x.domain)
        end
    end
    return selectedVar
end


text:
push!(LOAD_PATH, "$(homedir())/Development/")
using JSort

path = "$(homedir())/master/sortering/"
parameters = Parameters(path*"si28.yaml")
coefficients = featurealign2d("/home/erdos/master/sortering/sirius/")
savecoefficients_eΔe(parameters, coefficients...)
println(coefficients)

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
function plot_DrDtNr(Dr::Array{Float64,1},Dt::Array{Float64,1},Nr::Array{Float64,1},qNum::Int,tol::Array{Float64,1})
    fig = figure()
    p1=plot(tol,Dr,color="red",linewidth=2.0,linestyle="--")
    p2=plot(tol,Dt,color="blue",linewidth=2.0,linestyle="--")
    p3=plot(tol,Nr,color="green",linewidth=2.0,linestyle="--")
    axis("tight") # Fit the axis tightly to the plot
    ax = gca() # Get the handle of the current axis
    Title=string("Query",qNum)
    title(Title)
    xlabel("Tolerance")
    ylabel("No. of Documents")
    grid("on")
    #legend(loc="upper right",fancybox="true")
    legend([p1,p2,p3],["Dr","Dt","Nr"])
end

Colors=Array{String,1}
Colors=["red","blue","green","black","cyan","magenta","yellow"]
ColorCount=0
PlotCount=0

function plotNew_RecPrec(Rec::Array{Float64,1},Prec::Array{Float64,1},strMethod::String)
    fig = figure()
    global ColorCount+=1
    global PlotCount+=1
    p1=plot(Rec,Prec,color=Colors[ColorCount],linewidth=2.0,linestyle="--",label=strMethod)
    axis("tight") # Fit the axis tightly to the plot
    ax = gca() # Get the handle of the current axis
    Title=string("Precision Vs Recall")
    title(Title)
    xlabel("Recall")
    ylabel("Precision")
    grid("on")
    #push!(PlotIDs,string(PlotCount))
    #push!(MethodIDs,strMethod)
    legend(loc="upper right",fancybox="true")
end

function plotAdd_RecPrec(Rec::Array{Float64,1},Prec::Array{Float64,1},strMethod::String)
    global ColorCount+=1
    global PlotCount+=1
    PlotID=string("p",PlotCount)
    PlotID=plot(Rec,Prec,color=Colors[ColorCount],linewidth=2.0,linestyle="--",label=strMethod)
    #push!(PlotIDS,PlotID)
    #push!(MethodIDs,strMethod)
    #legend(PlotIDs,MethodIDs)
    legend(loc="upper right",fancybox="true")
end


text:
module FrostHelperCoreBerry

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/CoreBerry" FTCoreBerry(x::Integer, y::Integer, order::Integer=-1, checkpointID::Integer=-1, isIce::Bool=false)

const placements = Ahorn.PlacementDict(
	"Core Berry (Frost Helper)" => Ahorn.EntityPlacement(
		FTCoreBerry
	)
)

function Ahorn.selection(entity::FTCoreBerry)
	x, y = Ahorn.position(entity)

	return Ahorn.getSpriteRectangle("collectables/strawberry/normal00", x, y)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::FTCoreBerry, room::Maple.Room)
	x, y = Ahorn.position(entity)

	Ahorn.drawSprite(ctx, "collectables/strawberry/normal00", x, y)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
using Distributions, SpecialFunctions, Roots, LinearAlgebra
"""
GSmileOriginal : SMiLe rule (Faraji, et al., 2018)
    for the Gaussian task (section 2.3.1)
    - No knowledge of environment's prior (i.e. μ0 and σ0) for the calculation
        of Scc (it is there for the calculation of the χ and ν though)
    - Optimization to find γ (There is a closed formula for this task)
"""
struct GSmileOriginal <: GLearner
    m::Float64
    χ_0::Float64
    ν_0::Float64
    sigma::Float64 # std of observations
    χ_n::Array{Float64, 1}
    ν_n::Array{Float64, 1}
    mu_n::Array{Float64, 1} #parameter we are interested in
    Scc::Array{Float64, 1}
    γ_n::Array{Float64, 1}
end
function GSmileOriginal(;m = .1, mu_0 = 0., sigma_0 = 1., sigma = 0.1)
    χ_n = Array{Float64, 1}(undef, 1)
    ν_n = Array{Float64, 1}(undef, 1)
    mu_n = Array{Float64, 1}(undef, 1)
    ν_0 = calcν(sigma, sigma_0)
    χ_0 = calcχ(mu_0, ν_0, sigma)
    χ_n[1] = copy(χ_0)
    ν_n[1] = copy(ν_0)
    mu_n[1] = sigma * χ_n[1] / ν_n[1]
    GSmileOriginal(m, χ_0, ν_0, sigma, χ_n, ν_n, mu_n, [1.], [m/(1. + m)])
end
export GSmileOriginal
function update!(learnerG::GSmileOriginal, y)

    σ_n_squared = learnerG.sigma^2 / learnerG.ν_n[1]

    # The scaled likelihood is ~N(y, σ^2), assuming uniform prior (σ0 -> Inf)
    learnerG.Scc[1] = KL(learnerG.mu_n[1], σ_n_squared, y, learnerG.sigma^2)
    Bmax = KL(y, learnerG.sigma^2, learnerG.mu_n[1], σ_n_squared)
    B = learnerG.m * learnerG.Scc[1]/(1. + learnerG.m * learnerG.Scc[1]) * Bmax

    learnerG.γ_n[1] = find_γ0(y, learnerG.sigma^2, learnerG.mu_n[1], σ_n_squared, B)

    # Update using scaled likelihood: x0 = 0, and v0 = 0 (cause σ0 -> Inf)
    update_χ_ν!(learnerG, calcφ(learnerG, y), 1.)

    computemu_n!(learnerG)
end
export update!
function KL(μ1, σ1_squared, μ2, σ2_squared)
    σ_ratio = σ1_squared / σ2_squared
    max((μ1 - μ2)^2/(2. * σ2_squared) + 0.5 * (σ_ratio - 1. - log(σ_ratio)), 0.)
end
function update_χ_ν!(learnerG::GSmileOriginal, χ_0_y, ν_0_y)
    learnerG.χ_n[1] = (1. - learnerG.γ_n[1]) * learnerG.χ_n[1] + learnerG.γ_n[1] * χ_0_y
    learnerG.ν_n[1] = (1. - learnerG.γ_n[1]) * learnerG.ν_n[1] + learnerG.γ_n[1] * ν_0_y
end
function find_γ0(y, σ_squared, μ_n, σ_n_squared, B)
    # f = γ -> KL(γ .* betas .+ (1. - γ) .* alphas, alphas) - B
    ρ = σ_squared / σ_n_squared
    # w = 1. / (1. + (1. - γ)*ρ/γ)
    # μ_{n+1} = y*w + μ_n*(1. - w)
    #σ_{n+1} = 1. / (γ/σ^2 + (1. - γ)/σ_{n}^2)
    f = γ -> KL(y*get_w_smile(γ, ρ) + μ_n*(1. - get_w_smile(γ, ρ)), #μ_{n+1}
                get_σ_squared_smile(γ, σ_squared, σ_n_squared), #σ_{n+1}
                μ_n, σ_n_squared) - B
    if abs(f(0.)) < 5*eps()
        γ0 = 0.
    elseif abs(f(1.)) < 5*eps()
        γ0 = 1.
    else
        γ0 = find_zero(f, (0., 1.))
    end
    γ0
end
function get_w_smile(γ, ρ)
    1. / (1. + (1. - γ)*ρ/γ)
end
function get_σ_squared_smile(γ, σ_squared, σ_n_squared)
    1. /(γ/σ_squared + (1. - γ)/σ_n_squared)
end


text:
export learn_parameters, to_onehot

using CUDA
using LoopVectorization: @avx, vifelse

"""
LogisticCircuit Parameter learning through gradient descents
Note: data need to be DataFrame and Labels need to be in one-hot form.
"""
function learn_parameters(lc::LogisticCircuit, nc::Int, data, labels; num_epochs=25, step_size=0.01)
    bc = ParamBitCircuit(lc, nc, data)
    if isgpu(data)
        @assert isgpu(labels) "Data and labels must be both stored in either GPU or CPU."
        for _ = 1:num_epochs
            cl = class_likelihood_per_instance(bc, data)
            update_parameters_gpu(to_gpu(bc), data, labels, cl, step_size)
        end
    else
        @assert !isgpu(labels) "Data and labels must be both stored in either GPU or CPU."
        for _ = 1:num_epochs
            cl = class_likelihood_per_instance(bc, data)
            update_parameters_cpu(bc, data, labels, cl, step_size)
        end
    end
end

function update_parameters_cpu(bc, data, labels, cl, step_size)
    ne::Int = num_examples(data)
    nc::Int = size(bc.params, 2)
    params_lock::Threads.ReentrantLock = Threads.ReentrantLock()

    @inline function on_edge_binary(flows, values, prime, sub, element, grandpa, single_child, weights)
        lock(params_lock) do # TODO: move lock to inner loop?
            for i = 1:size(flows, 1)
                @inbounds edge_flow = values[i, prime] & values[i, sub] & flows[i, grandpa]
                first_true_bit = trailing_zeros(edge_flow) + 1
                last_true_bit = 64 - leading_zeros(edge_flow)
                @simd for j = first_true_bit:last_true_bit
                    ex_id = ((i - 1) << 6) + j
                    if get_bit(edge_flow, j)
                        for class = 1:nc
                            @inbounds bc.params[element, class] -= (cl[ex_id, class] - labels[ex_id, class]) * step_size
                        end
                    end
                end
            end
        end
    end

    @inline function on_edge_float(flows, values, prime, sub, element, grandpa, single_child, weights)
        lock(params_lock) do # TODO: move lock to inner loop?
            @avx for i = 1:size(flows, 1)
                @inbounds edge_flow = values[i, prime] * values[i, sub] / values[i, grandpa] * flows[i, grandpa]
                edge_flow = vifelse(isfinite(edge_flow), edge_flow, zero(eltype(flows)))
                for class = 1:nc
                    @inbounds bc.parames[element, class] -= (cl[i, class] - labels[i, class]) * edge_flow * step_size
                end
            end
        end
        nothing
    end

    if isbinarydata(data)
        v,f = satisfies_flows(bc.bitcircuit, data; on_edge = on_edge_binary)
    else
        @assert isfpdata(data) "Only floating point and binary data are supported"
        v,f = satisfies_flows(bc.bitcircuit, data; on_edge = on_edge_float)
    end

    nothing
end

function update_parameters_gpu(bc, data, labels, cl, step_size)
    ne::Int = num_examples(data)
    nc::Int = size(bc.params, 2)
    cl_device = CUDA.cudaconvert(cl)
    label_device = CUDA.cudaconvert(labels)
    params_device = CUDA.cudaconvert(bc.params)

    @inline function on_edge_binary(flows, values, prime, sub, element, grandpa, chunk_id, edge_flow, single_child, weights)
        first_true_bit = 1 + trailing_zeros(edge_flow)
        last_true_bit = 64 - leading_zeros(edge_flow)
        for j = first_true_bit:last_true_bit
            if get_bit(edge_flow, j)
                ex_id = ((chunk_id - 1) << 6) + j
                for class = 1:nc
                    CUDA.@atomic params_device[element, class] -= (cl_device[ex_id, class] - label_device[ex_id, class]) * step_size
                end
            end
        end
        nothing
    end

    @inline function on_edge_float(flows, values, prime, sub, element, grandpa, ex_id, edge_flow, single_child, weights)
        for class = 1:nc
            CUDA.@atomic params_device[element, class] -= (cl_device[ex_id, class] - label_device[ex_id, class]) * edge_flow * step_size
        end
        nothing
    end

    if isbinarydata(data)
        v,f = satisfies_flows(bc.bitcircuit, data; on_edge = on_edge_binary)
    else
        @assert isfpdata(data) "Only floating point and binary data are supported"
        v,f = satisfies_flows(bc.bitcircuit, data; on_edge = on_edge_float)
    end
    CUDA.unsafe_free!(v) # save the GC some effort
    CUDA.unsafe_free!(f) # save the GC some effort

    nothing
end

function to_onehot(labels::Vector, nc::Integer)
    ne = length(labels)
    one_hot_labels = zeros(Float32, ne, nc)
    for (i, j) in enumerate(labels)
        one_hot_labels[i, j + 1] = 1.0
    end
    one_hot_labels
end


text:
# Compute the longest streak of `nuc` in vector `seq`
function longeststreak(seq::Vector{UInt8},nuc::UInt8)
    len=0
    current=0
    last=0
    count=0
    for i in linearindices(seq)
        if seq[i]==nuc
            count+=1
            current+=1
            if current > len
                len = current
                last = i
            end
        else
            current=0
        end
    end
    return len, last, count
end

# Compute a Hamming distance between `id1` and `id2`, assuming that
# both are of length `len`
function hamming(id1::UInt,id2::UInt,len)
    id1-=1
    id2-=1
    dist::UInt8 = 0
    for n in 1:len
        if (id1&0b11) $ (id2&0b11) != 0
            dist+=1
        end
        id1>>=2
        id2>>=2
    end
    return dist
end

function getminhammingdistance(codes::AbstractVector{UInt},N)
    nearest = zeros(UInt8,length(codes))
    @inbounds for i1 in 1:length(codes)
        nearest[i1] = N       # start with a maximal distance
        for i2 in 1:(i1-1)
            dist = hamming(codes[i1],codes[i2],N)
            nearest[i1] = min(nearest[i1],dist)
            # the minimum has been reached
            if dist == 1
                continue
            end
        end
    end
    return nearest
end


text:
VERSION >= v"0.4.0" && __precompile__()

module Metamath

using Compat
using CSV

import Base: empty!, show

export mmverify!

type MetamathException <: Exception
  text::AbstractString
end
metamath_error(x...) = throw(MetamathException(string(x...)))
metamath_warn(x...) = warn("(Metamath) ",x...)
macro warn_and_ret(rv,x...)
  :(metamath_warn($x...) ; return $rv)
end

typealias Expression Vector{Symbol}
typealias Hypothesis Pair{Expression,Bool}

# An axiom or theorem
immutable Assertion
  hypotheses::Vector{Symbol}
  disjvars::Set{Pair{Symbol,Symbol}}
  expression::Expression
  Assertion(expression) = new(Symbol[],Set{Pair{Symbol,Symbol}}(),expression)
end

type Scope
  activevariables::Set{Symbol}
  activehyp::Vector{Symbol}
  disjvars::Vector{Set{Symbol}}
  floatinghyp::Dict{Symbol,Symbol}
  Scope() = new(Set{Symbol}(),Symbol[],Set{Symbol}[],Dict{Symbol,Symbol}())
end

# An environment holding truths and definitions accumulated by verifier
immutable Environment
  filenames::Set{ASCIIString}
  tokens::Vector{Symbol}
  compressedproofs::Vector{ASCIIString}
  constants::Set{Symbol}
  hypotheses::Dict{Symbol,Hypothesis}
  variables::Set{Symbol}
  assertions::Dict{Symbol,Assertion}
  scopes::Vector{Scope}
  substitutionssrc::Vector{Symbol}
  substitutionsdst::Vector{Expression}
  dirty::Vector{Bool}
  Environment() = new(Set{ASCIIString}(),Symbol[],ASCIIString[],Set{Symbol}(),
    Dict{Symbol,Hypothesis}(), Set{Symbol}(), Dict{Symbol,Assertion}(),
    Scope[Scope()],Symbol[],Expression[],Bool[false])
end

"""a global Metamath.Environment variable for defualt use by mmverify!"""
const globalenv = Environment()

# true if label already used in hypotheses or assertions in env
labelused(env,label) = haskey(env.hypotheses,label) || haskey(env.assertions,label)

# return floating hypothesis in which variable appears
function getfloatinghyp(env,variable)
  for scope in env.scopes
    if haskey(scope.floatinghyp,variable)
      return scope.floatinghyp[variable]
    end
  end
  return :dummytoken
end

# Character set not allowed in Metamath labels
const tokenspecials = Set(['.','-','_'])

# true if ch is a Metamath whitespace
ismmws(ch) = return ch == ' ' || ch == '\n' || ch == '\t' || ch == '\f' || ch == '\r'
# true if sym represent a string which is legal Metamath label
islabeltoken(str) = findfirst(c->!(c in tokenspecials || isalnum(c)),str)==0
# true if str is a legal Metamath symbol
ismathsymboltoken(str::ASCIIString) = findfirst(str,'$')==0
# true is sym represents a string which is legal Metamth symbol
ismathsymboltoken(sym::Symbol) = ismathsymboltoken(string(sym))
# true if c is not uppercase or question mark
isntupperorq(c) = !(isupper(c) || c=='?')
# true is str has only upper-case letters or question marks
containsonlyupperorq(str) = findfirst(isntupperorq,str)==0

# true if str is an active variable according to env
isactivevariable(env,str) = any(scope->str in scope.activevariables,env.scopes)
# true if str is an active hypothesis according to env
isactivehyp(env,str) = any(scope->str in scope.activehyp,env.scopes)

# true if str1 and str2 should be distinct according to env
function isdvr(env,str1,str2)
  if str1==str2
    return false
  end
  for scope in env.scopes
    for dvset in scope.disjvars
      if str1 in dvset && str2 in dvset
        return true
      end
    end
  end
  return false
end

# Output a string summary representation of an Environment type
function show(io::IO,env::Environment)
  print(io,"<Metamath environment: ")
  print(io,join([length(getfield(env,f)) for f in fieldnames(Environment)],","))
  print(io,">")
end

# Return a single token from buffer discarding surrounding whitespace.
function nexttoken(stream,buf::IOBuffer = IOBuffer())
  c = ' '
  while !eof(stream) && ismmws(c)
    c = read(stream,Char)
  end
  while !ismmws(c)
    if c<'!' || c>'~'
      @warn_and_ret("","Invalid character read with code $(UInt8(c))")
    end
    write(buf,c)
    if eof(stream)
      break
    end
    c = read(stream,Char)
  end
  return takebuf_string(buf)
end

# Parse filename into tokens and compressed proofs and store
# parsed output into environment. Discard comments and read
# other included files.
function readtokens!(env, filename; use_mmap::Bool=true)
  if filename in env.filenames
    return nothing
  else
    push!(env.filenames,filename)
  end

  isfile(filename) || throw(ArgumentError("\"$filename\" is not a valid or existing filename"))
  stream = CSV.UnsafeBuffer(use_mmap ? Mmap.mmap(filename) : open(readbytes,filename))

  incomment = false
  infileinclusion = false
  newfilename = ""
  state = 0

  strtoken = ""
  buf = IOBuffer()
  compressedbuf = IOBuffer()
  while (strtoken = nexttoken(stream,buf))!=""
    if incomment
      if strtoken == "\$)"
        incomment = false
        continue
      end
      if !isempty(search(strtoken,"\$("))
        metamath_error("Character \$( found in comment")
      end
      if !isempty(search(strtoken,"\$)"))
        metamath_error("Character \$) found in comment")
      end
      continue
    end

    if strtoken=="\$("
      incomment = true
      continue
    end

    if infileinclusion
      if newfilename == ""
        if findfirst(strtoken,'$')!=0
          metamath_error("Filename contains a \$")
        end
        newfilename = strtoken
        continue
      else
        if strtoken != "\$]"
          metamath_error("Didn't find closing file inclusion delimeter")
        end
        readtokens(env,newfilename,use_mmap=user_mmap)
        newfilename = ""
        continue
      end
    end

    if strtoken == "\$["
      infileinclusion = true
      continue
    end

    if state == 0 && strtoken == "\$p"
      state = 1
    elseif state == 1 && strtoken == "\$="
      state = 2
    elseif state == 2 && strtoken == "("
      state = 3
    elseif state == 3 && strtoken == ")"
      state = 4
    elseif state == 4
      if strtoken != "\$."
        if !containsonlyupperorq(strtoken)
          metamath_error("Compressed proof $strtoken contains bogus character")
        end
        write(compressedbuf, strtoken)
      else
        push!(env.compressedproofs,takebuf_string(compressedbuf))
        state = 0
      end
      continue
    elseif state > 0 && strtoken == "\$."
      state = 0
    end
    push!(env.tokens,Symbol(strtoken))
  end
  if !eof(stream)
    metamath_error("Error reading from $filename")
  end
  if incomment
    metamath_error("Unclosed comment in $filename")
  end
  if infileinclusion
    metamath_error("Unfinished file inclusion command in $filename")
  end
  if state > 0 || nb_available(compressedbuf) != 0
    metamath_error("Compressed proof parsing confused in $filename")
  end
  return nothing
end

# Construct Assertion from expression by determining mandatory
# hypotheses and disjoint variable restrictions. Assertion is inserted
# into environment and returned
function constructassertion(env,label, expression)
  assertion = Assertion(expression)
  varsused = Set{Symbol}()
  for token in expression
    if token in env.variables
      push!(varsused, token)
    end
  end
  for scope in reverse(env.scopes)
    for hyp in reverse(scope.activehyp)
      ind = Base.ht_keyindex(env.hypotheses,hyp)
      if ind != -1
        if last(env.hypotheses.vals[ind]) && (first(env.hypotheses.vals[ind])[2] in varsused)
          unshift!(assertion.hypotheses,hyp)
        elseif !last(env.hypotheses.vals[ind])
          unshift!(assertion.hypotheses,hyp)
          for token in first(env.hypotheses.vals[ind])
            if token in env.variables
              push!(varsused, token)
            end
          end
        end
      end
    end
  end
  varusedvec = collect(varsused)
  for scope in env.scopes
    for dvset in scope.disjvars
      dvec = intersect(varusedvec,dvset)
      if length(dvec)<2 continue ; end
      for i=1:length(dvec)
        for j=(i+1):length(dvec)
          push!(assertion.disjvars,dvec[i]=>dvec[j])
        end
      end
    end
  end
  return assertion
end

# Return token expressions appearing in various statements
function readexpression!(env,stattype::Char, label, terminator)
  expression = Expression()
  if isempty(env.tokens)
    metamath_error("Unfinished \$$stattype statement $label")
  end
  token = shift!(env.tokens)
  if !(token in env.constants)
    metamath_error("First symbol in \$$stattype statement $label"*
      " is $token which is not a constant")
  end
  push!(expression, token)
  while !isempty(env.tokens)
    token = shift!(env.tokens)
    if token == terminator
      break
    end
    if !(token in env.constants) && getfloatinghyp(env,token)==:dummytoken
      metamath_error("In \$$stattype statement $label"*
        " $token found which is not a constant or a variable"*
        " in an active \$f statement")
    end
    push!(expression, token)
  end
  if isempty(env.tokens) && token != terminator
    metamath_error("Unfinished \$$stattype statement $label")
  end
  return expression
end

# Return number vector encoded alphabetically in compressed proof
function getproofnumbers(label, proof)
  proofnumbers = Vector{Int}()
  sizehint!(proofnumbers,length(proof))
  justgotnum = false
  num = zero(UInt)
  for c in proof
    if c <= 'T'
      addval = c-('A'-1)
      if num > typemax(UInt)/20 || 20*num > typemax(UInt)-addval
        metamath_error("Overflow computing numbers in compressed"*
          " proof of $label")
      end
      push!(proofnumbers, 20*num+addval)
      num = zero(UInt)
      justgotnum = true
    elseif c <= 'Y'
      addval = c-'T'
      if num > typemax(UInt)/5 || 5*num > typemax(UInt)-addval
        metamath_error("Overflow computing numbers in compressed"*
          " proof of $label")
      end
      num = 5*num+addval
      justgotnum = false
    else # it must be 'Z'
      if !justgotnum
        metamath_error("Stray Z found in compressed proof of $label")
      end
      push!(proofnumbers,0)
      justgotnum = false
    end
  end
  if num != 0
    metamath_error("Compressed proof of theorem $label ends"*
      " in unfinished number")
  end
  return proofnumbers
end

# Perform substitutions of variables specified in substmap and return
# resulting Expression
function makesubstitution(original, substmapsrc, substmapdst)
  destination = Expression()
  for token in original
    ind = findfirst(substmapsrc,token)
    if ind != 0
      append!(destination,substmapdst[ind])
    else
      push!(destination,token)
    end
  end
  return destination
end

# Verify a proof step referencing assertion (i.e. not hypothesis)
function verifyassertionref(env, thlabel, reflabel, stack)
  assertion = env.assertions[reflabel]
  if length(stack)<length(assertion.hypotheses)
    metamath_error("In proof of theorem $thlabel not enough"*
      " items found on stack")
  end
  base = length(stack)-length(assertion.hypotheses)
  # empty!(substitutions)
  empty!(env.substitutionssrc)
  empty!(env.substitutionsdst)
  for i=1:length(assertion.hypotheses)
    hypothesis = env.hypotheses[assertion.hypotheses[i]]
    if last(hypothesis)
      # Floating hypothesis
      if first(hypothesis)[1] != stack[base+i][1]
        metamath_error("In proof of theorem $thlabel the"*
          " unification failed")
      end
      subst = copy(stack[base+i])
      shift!(subst)
      # substitutions[first(hypothesis)[2]] = subst
      push!(env.substitutionssrc,first(hypothesis)[2])
      push!(env.substitutionsdst,subst)
    else
      # Essential hypothesis
      dest = makesubstitution(first(hypothesis),env.substitutionssrc,env.substitutionsdst)
      if dest != stack[base+i]
        metamath_error("In proof of theorem $thlabel unification failed")
      end
    end
  end
  resize!(stack,base)
  # Verify disjoint variable conditions
  for vpair in assertion.disjvars
    exp1 = env.substitutionsdst[findfirst(env.substitutionssrc,first(vpair))]
    exp2 = env.substitutionsdst[findfirst(env.substitutionssrc,last(vpair))]
    # exp2 = substitutions[last(vpair)]
    exp1vars = Set{Symbol}()
    for token in exp1
      if token in env.variables
        push!(exp1vars,token)
      end
    end
    exp2vars = Set{Symbol}()
    for token in exp2
      if token in env.variables
        push!(exp2vars,token)
      end
    end
    for exp1var in exp1vars
      for exp2var in exp2vars
        if !isdvr(env,exp1var,exp2var)
          metamath_error("In proof of theorem $thlabel disjoint"*
            " variable restriction violated")
        end
      end
    end
  end
  dest = makesubstitution(assertion.expression, env.substitutionssrc, env.substitutionsdst)
  push!(stack,dest)
  return nothing
end

# Verify regular proof
function verifyregularproof(env, label, theorem, proof)
  stack = Vector{Expression}()
  for proofstep in proof
    if haskey(env.hypotheses,proofstep)
      push!(stack, first(env.hypotheses[proofstep]))
      continue
    end
    verifyassertionref(env, label, proofstep, stack)
  end
  if length(stack) != 1
    metamath_error("Proof of theorem $label does not end"*
      " with only one item on the stack")
  end
  if stack[1] != theorem.expression
    metamath_error("Proof of theorem $label proves wrong statement")
  end
  return nothing
end

# Verify compressed proof
function verifycompressedproof(env,label,theorem,labels,proofnumbers)
  stack = Vector{Expression}()
  mandhypt = length(theorem.hypotheses)
  labelt = mandhypt+length(labels)
  savedsteps = Vector{Expression}()
  for pn in proofnumbers
    if pn==0
      push!(savedsteps,last(stack))
      continue
    end
    if pn<=mandhypt
      push!(stack,first(env.hypotheses[theorem.hypotheses[pn]]))
    elseif pn<=labelt
      proofstep = labels[pn-mandhypt]
      if haskey(env.hypotheses,proofstep)
        push!(stack,first(env.hypotheses[proofstep]))
        continue
      end
      verifyassertionref(env,label, proofstep, stack)
    else
      if pn > labelt+length(savedsteps)
        metamath_error("Number in compressed proof of $label"*
          " is too high")
      end
      push!(stack, savedsteps[pn-labelt])
    end
  end
  if length(stack) != 1
    metamath_error("Proof of theorem $label does not end with"*
      " only one item on the stack")
  end
  if first(stack) != theorem.expression
    metamath_error("Proof of theorem $label proves the wrong statement")
  end
  return nothing
end

# Parse $p statement
function parsep!(env,label)
  newtheorem = readexpression!(env,'p', label, Symbol("\$="))
  assertion = constructassertion(env,label, newtheorem)
  env.assertions[label] = assertion
  if isempty(env.tokens)
    metamath_error("Unfinished \$p statement")
  end
  token = shift!(env.tokens)
  if token == Symbol("(")
    # Compressed proof
    labels = Vector{Symbol}()
    while !isempty(env.tokens)
      token = shift!(env.tokens)
      if token == Symbol(")")
        break
      end
      push!(labels, token)
      if token == label
        metamath_error("Proof of theorem $label refers to itself")
      elseif token in assertion.hypotheses
        metamath_error("Compressed proof of theorem $label has"*
          " mandatory hypothesis $token in label list")
      elseif !haskey(env.assertions,token) && !isactivehyp(env,token)
        metamath_error("Proof of theorem $label refers to $token"*
          " which is not an active statement")
      end
    end
    if isempty(env.tokens)
      metamath_error("Unfinished \$p statement $label")
    end
    if isempty(env.compressedproofs)
      metamath_error("Missing compressed proof for $label")
    end
    proof = shift!(env.compressedproofs)
    if isempty(proof)
      metamath_error("Theorem $label has no proof")
    end
    if findfirst(proof,'?')!=0
      metamath_error("Proof of theorem $label is incomplete")
    end
    proofnumbers = getproofnumbers(label, proof)
    verifycompressedproof(env,label, assertion, labels, proofnumbers)
  else
    # Regular (uncompressed) proof
    proof = Vector{Symbol}()
    push!(proof,token)
    incomplete = false
    token = :dummytoken
    while !isempty(env.tokens)
      token = shift!(env.tokens)
      if token == Symbol("\$.")
        break
      end
      push!(proof, token)
      if token == Symbol("?")
        incomplete = true
      elseif token == label
        metamath_error("Proof of theorem $label refers to itself")
      elseif !haskey(env.assertions,token) && !isactivehyp(env,token)
        metamath_error("Proof of theorem $label refers to $token"*
          " which is not an active statement")
      end
    end
    if isempty(env.tokens) && token != Symbol("\$.")
      metamath_error("Unfinished \$p statement")
    end
    if isempty(proof)
      metamath_error("Theorem $label has no proof")
    end
    if incomplete
      @warn_and_ret(nothing, "Warning: Proof of theorem $label is incomplete")
    end
    verifyregularproof(env,label, assertion, proof)
  end
  return nothing
end

# Parse $e statement
function parsee!(env,label)
  newhyp = readexpression!(env,'e', label, Symbol("\$."))
  env.hypotheses[label] = newhyp=>false
  push!(last(env.scopes).activehyp,label)
  return nothing
end

# Parse $a statement
function parsea!(env,label)
  newaxiom = readexpression!(env,'a',label,Symbol("\$."))
  assertion = constructassertion(env,label, newaxiom)
  env.assertions[label] = assertion
  return nothing
end

# Parse $f statement
function parsef!(env,label)
  if isempty(env.tokens)
    metamath_error("Unfinished \$f statement $label")
  end
  typetoken = shift!(env.tokens)
  if !(typetoken in env.constants)
    metamath_error("First symbol in \$f statement $label is $typetoken"*
      " which is not a constant")
  end
  if isempty(env.tokens)
    metamath_error("Unfinished \$f statement $label")
  end
  variable = shift!(env.tokens)
  if !isactivevariable(env,variable)
    metamath_error("Second symbol in \$f statement $label is $variable"*
      " which is not an active variable")
  end
  if !(getfloatinghyp(env,variable)==:dummytoken)
    metamath_error("The variable $variable appears in a second"*
      " \$f statement $label")
  end
  if isempty(env.tokens)
    metamath_error("Unfinished \$f statement $label")
  end
  token = shift!(env.tokens)
  if token != Symbol("\$.")
    metamath_error("Expected end of \$f statement but found $token")
  end
  newhyp = Expression()
  push!(newhyp, typetoken)
  push!(newhyp, variable)
  env.hypotheses[label] = newhyp=>true
  push!(last(env.scopes).activehyp,label)
  last(env.scopes).floatinghyp[variable] = label
  return nothing
end

# Parse labeled statement
function parselabel!(env,label)
  if label in env.constants
    metamath_error("Attempt to reuse constant $label as a label")
  end
  if label in env.variables
    metamath_error("Attempt to reuse variable $label as a label")
  end
  if labelused(env,label)
    metamath_error("Attempt to reuse label $label")
  end
  if isempty(env.tokens)
    metamath_error("Unfinished labeled statement")
  end
  labeltype = shift!(env.tokens)
  if labeltype == Symbol("\$p")
    parsep!(env,label)
  elseif labeltype == Symbol("\$e")
    parsee!(env,label)
  elseif labeltype == Symbol("\$a")
    parsea!(env,label)
  elseif labeltype == Symbol("\$f")
    parsef!(env,label)
  else
    metamath_error("Unexpected token $labeltype encountered")
  end
  return nothing
end

# Parse $d statement
function parsed!(env)
  dvars = Set{Symbol}()
  token = :dummytoken
  while !isempty(env.tokens)
    token = shift!(env.tokens)
    if token == Symbol("\$.")
      break
    end
    if !isactivevariable(env,token)
      metamath_error("Token $token is not an active variable but was"*
        " found in a \$d statement")
    end
    if token in dvars
      metamath_error("\$d statement mentions $token twice")
    end
    push!(dvars, token)
  end
  if isempty(env.tokens) && token != Symbol("\$.")
    metamath_error("Unfinished \$d statement")
  end
  if length(dvars)<2
    metamath_error("Not enough variables in \$d statement")
  end
  push!(last(env.scopes).disjvars, dvars)
  return nothing
end

# Parse $c statement
function parsec!(env)
  if length(env.scopes)>1
    metamath_error("\$c statement occurs in inner block")
  end
  token = ""
  listempty = true
  while !isempty(env.tokens)
    token = shift!(env.tokens)
    if token == Symbol("\$.")
      break
    end
    listempty = false
    if !ismathsymboltoken(token)
      metamath_error("Attempt to declare $token as a constant")
    end
    if token in env.variables
      metamath_error("Attempt to redeclare variable $token as a constant")
    end
    if labelused(env,token)
      metamath_error("Attempt to reuse label $token as a constant")
    end
    if token in env.constants
      metamath_error("Attempt to redeclare constant $token")
    end
    push!(env.constants, token)
  end
  if isempty(env.tokens) && token != Symbol("\$.")
    metamath_error("Unterminated \$c declaration")
  end
  if listempty
    metamath_error("Empty \$c statement")
  end
  return nothing
end

# Parse $v statement
function parsev!(env)
  token = ""
  listempty = true
  while !isempty(env.tokens)
    token = shift!(env.tokens)
    if token == Symbol("\$.")
      break
    end
    listempty = false
    if !ismathsymboltoken(token)
      metamath_error("Attempt to declare $token as a variable")
    end
    if token in env.constants
      metamath_error("Attempt to redeclare constant $token as a variable")
    end
    if labelused(env,token)
      metamath_error("Attempt to reuse label $token as a variable")
    end
    if isactivevariable(env,token)
      metamath_error("Attempt to redeclare active variable $token")
    end
    push!(env.variables,token)
    push!(last(env.scopes).activevariables,token)
  end
  if isempty(env.tokens) && token != "\$."
    metamath_error("Unterminated \$v statement")
  end
  if listempty
    metamath_error("Empty \$v statement")
  end
  return nothing
end

"""`mmverify!(env,filename)` verifies database in filename using and then return env"""
function mmverify!(env::Environment,filename::AbstractString;
             use_mmap::Bool=true)
  if env.dirty[1]
    metamath_warn("Starting verification of $filename with dirty"*
      " environment")
  else
    env.dirty[1] = true
  end
  readtokens!(env,filename)

  while length(env.tokens)>0
    token = shift!(env.tokens)
    if islabeltoken(string(token))
      parselabel!(env,token)
    elseif token==Symbol("\$d")
      parsed!(env)
    elseif token==Symbol("\${")
      push!(env.scopes,Scope())
    elseif token==Symbol("\$}")
      pop!(env.scopes)
      if isempty(env.scopes)
        metamath_error("\$} without corresponding \${")
      end
    elseif token==Symbol("\$c")
      parsec!(env)
    elseif token==Symbol("\$v")
      parsev!(env)
    else
      metamath_error("Unexpected token $token encountered")
    end
  end
  if length(env.scopes)>1
    metamath_error("\${ without corresponding \$}")
  end
  empty!(env.substitutionssrc)
  empty!(env.substitutionsdst)
  env.dirty[1] = false
  return nothing
end

"""`main(filename ; use_mmap=true, use_globalenv=true)`

Verify contents of filename with a new environment or `globalenv` by default.
`use_mmap`, true by default, makes file processing faster, may be impossible
with other data sources.

Return value: the Metamath environment after processing file.

Example: `Metamath.main("set.mm")`
"""
function main(filename::AbstractString;
             use_mmap::Bool=true,use_globalenv::Bool=true)
  if use_globalenv
    env = globalenv
  else
    env = Environment()
  end
  mmverify!(env,filename)
  return env
end

# Restore scope to initial state without reallocating scope
function empty!(scope::Scope)
  empty!(scope.activevariables)
  empty!(scope.activehyp)
  empty!(scope.disjvars)
  empty!(scope.floatinghyp)
  return scope
end

"""`empty!(env)` returns `env` to initial empty state and returns `env`"""
function empty!(env::Environment)
  empty!(env.filenames)
  empty!(env.tokens)
  empty!(env.compressedproofs)
  empty!(env.constants)
  empty!(env.hypotheses)
  empty!(env.variables)
  empty!(env.assertions)
  resize!(env.scopes,1)
  empty!(env.scopes[1])
  empty!(env.substitutionssrc)
  empty!(env.substitutionsdst)
  env.dirty[1] = false
  return env
end

end # module


text:
module Destructuring
using  ..Helper
export destructure, @letds, @macrods, @anonds, @funds

include("Function.jl")
include("Applications.jl")

import .Function:     destructure
import .Applications: @letds, @macrods, @anonds, @funds

end


text:
# This file includes the Model object

"""
    Node(component, idx, label)

Constructs a model `Node` with `component`. `idx` is the index and `label` is label of `Node`.
"""
struct Node{CP, L}
    component::CP
    idx::Int
    label::L
end

show(io::IO, node::Node) = print(io, "Node(component:$(node.component), idx:$(node.idx), label:$(node.label))")

"""
    Branch(nodepair, indexpair, links)

Constructs a `Branch` connecting the first and second element of `nodepair` with `links`. `indexpair` determines the subindices by which the elements of `nodepair` are connected.
"""
struct Branch{NP, IP, LN<:AbstractVector{<:Link}}
    nodepair::NP
    indexpair::IP
    links::LN
end

show(io::IO, branch::Branch) = print(io, "Branch(nodepair:$(branch.nodepair), indexpair:$(branch.indexpair), ",
    "links:$(branch.links))")

"""
    Model(components::AbstractVector)

Constructs a `Model` whose with components `components` which are of type `AbstractComponent`.

    Model()

Constructs a `Model` with empty components. After the construction, components can be added to `Model`.

!!! warning
    `Model`s are units that can be simulated. As the data flows through the branches i.e. input output busses of the components, its is important that the components must be connected to each other. See also: [`simulate!`](@ref)
"""
struct Model{GR, ND, BR, CK, TM, CB}
    graph::GR
    nodes::ND
    branches::BR
    clock::CK
    taskmanager::TM
    callbacks::CB
    name::Symbol
    id::UUID
    function Model(nodes::AbstractVector=[], branches::AbstractVector=[];
        clock=Clock(0, 0.01, 1.), callbacks=nothing, name=Symbol())
        graph = SimpleDiGraph()
        taskmanager = TaskManager()
        new{typeof(graph), typeof(nodes), typeof(branches), typeof(clock), typeof(taskmanager),
            typeof(callbacks)}(graph, nodes, branches, clock, taskmanager, callbacks, name, uuid4())
    end
end

show(io::IO, model::Model) = print(io, "Model(numnodes:$(length(model.nodes)), ",
    "numedges:$(length(model.branches)), timesettings=($(model.clock.t), $(model.clock.dt), $(model.clock.tf)))")

##### Addinng nodes and branches.
"""
    addnode!(model, component; label=nothing)

Adds a node to `model`. Component is `component` and `label` is `label` the label of node. Returns added node.

# Example
```jldoctest
julia> model = Model()
Model(numnodes:0, numedges:0, timesettings=(0.0, 0.01, 1.0))

julia> addnode!(model, SinewaveGenerator(), label=:gen)
Node(component:SinewaveGenerator(amp:1.0, freq:1.0, phase:0.0, offset:0.0, delay:0.0), idx:1, label:gen)
```
"""
function addnode!(model::Model, component::AbstractComponent; label=nothing)
    label === nothing || label in [node.label for node in model.nodes] && error(label," is already assigned.")
    node = Node(component, length(model.nodes) + 1, label)
    push!(model.nodes, node)
    register(model.taskmanager, component)
    add_vertex!(model.graph)
    node
end

"""
    getnode(model, idx::Int)

Returns node of `model` whose index is `idx`.

    getnode(model, label)

Returns node of `model` whose label is `label`.

# Example
```jldoctest
julia> model = Model()
Model(numnodes:0, numedges:0, timesettings=(0.0, 0.01, 1.0))

julia> addnode!(model, SinewaveGenerator(), label=:gen)
Node(component:SinewaveGenerator(amp:1.0, freq:1.0, phase:0.0, offset:0.0, delay:0.0), idx:1, label:gen)

julia> addnode!(model, Gain(), label=:gain)
Node(component:Gain(gain:1.0, input:Inport(numpins:1, eltype:Inpin{Float64}), output:Outport(numpins:1, eltype:Outpin{Float64})), idx:2, label:gain)

julia> getnode(model, :gen)
Node(component:SinewaveGenerator(amp:1.0, freq:1.0, phase:0.0, offset:0.0, delay:0.0), idx:1, label:gen)

julia> getnode(model, 2)
Node(component:Gain(gain:1.0, input:Inport(numpins:1, eltype:Inpin{Float64}), output:Outport(numpins:1, eltype:Outpin{Float64})), idx:2, label:gain)
```
"""
getnode(model::Model, idx::Int) = model.nodes[idx]
getnode(model::Model, label) = filter(node -> node.label === label, model.nodes)[1]

"""
    $(SIGNATURES)

Returns the component of `model` corresponding to `specifier` that can be either index or label of the component.
"""
getcomponent(model::Model, specifier) = getnode(model, specifier).component

function register(taskmanager, component)
    triggerport, handshakeport = taskmanager.triggerport, taskmanager.handshakeport
    triggerpin, handshakepin = Outpin(), Inpin{Bool}()
    connect!(triggerpin, component.trigger)
    connect!(component.handshake, handshakepin)
    push!(triggerport.pins, triggerpin)
    push!(handshakeport.pins, handshakepin)
    taskmanager.pairs[component] = nothing
end

"""
    addbranch!(model::Model, branch::Branch)

Adds `branch` to branched of `model`.
"""
function addbranch!(model::Model, nodepair::Pair, indexpair::Pair=(:)=>(:))
    srcnode, dstnode = getnode(model, nodepair.first), getnode(model, nodepair.second)
    links = connect!(srcnode.component.output[indexpair.first], dstnode.component.input[indexpair.second])
    typeof(links) <: AbstractVector{<:Link} || (links = [links])
    srcidx, dstidx = srcnode.idx, dstnode.idx
    branch =  Branch(srcidx => dstidx, indexpair, links)
    push!(model.branches, branch)
    add_edge!(model.graph, srcidx, dstidx)
    branch
end

getbranch(model::Model, nodepair::Pair{Int, Int}) = filter(branch -> branch.nodepair == nodepair, model.branches)[1]
getbranch(model::Model, nodepair::Pair{Symbol, Symbol}) =
    getbranch(model, getnode(model, nodepair.first).idx => getnode(model, nodepair.second).idx)

"""
    $(SIGNATURES)

Returns the links of `model` corresponding to the `pair` which can be a pair of integers or symbols to specify the source and destination nodes of the branch.
"""
getlinks(model::Model, pair) = getbranch(model, nodepair).links

"""
    deletebranch!(model::Model, branch::Branch)

Deletes `branch` from branched of `model`.

    deletebranch!(model::Model, srcnode::Node, dstnode::Node)

Deletes branch between `srcnode` and `dstnode` of the `model`.
"""
function deletebranch!(model::Model, nodepair::Pair{Int, Int})
    srcnode, dstnode = getnode(model, nodepair.first), getnode(model, nodepair.second)
    branch = getbranch(model, nodepair)
    srcidx, dstidx = branch.indexpair
    disconnect!(srcnode.component.output[srcidx], dstnode.component.input[dstidx])
    deleteat!(model.branches, findall(br -> br == branch, model.branches))
    rem_edge!(model.graph, srcnode.idx, dstnode.idx)
    branch
end
deletebranch!(model::Model, nodepair::Pair{Symbol, Symbol}) =
    deletebranch!(model, getnode(model, nodepair.first).idx, getnode(model, nodepair.second).idx)

##### Model inspection.
"""
    inspect!(model::Model)

Inspects the `model`. If `model` has some inconsistencies such as including algebraic loops or unterminated busses and
error is thrown.
"""
function inspect!(model, breakpoints::Vector{Int}=Int[])
    # Check unbound pins in ports of componensts
    checknodeports(model)

    # Check links of the model
    checkchannels(model)

    # Break algebraic loops if there exits.
    loops = getloops(model)
    if !isempty(loops)
        msg = "\tThe model has algrebraic loops:$(loops)"
        msg *= "\n\t\tTrying to break these loops..."
        @info msg
        while !isempty(loops)
            loop = popfirst!(loops)
            if hasmemory(model, loop)
                @info "\tLoop $loop has a Memory component.  The loops is broken"
                continue
            end
            breakpoint = isempty(breakpoints) ? length(loop) : popfirst!(breakpoints)
            breakloop!(model, loop, breakpoint)
            @info "\tLoop $loop is broken"
            loops = getloops(model)
        end
    end

    # Return model
    model
end

hasmemory(model, loop) = any([getnode(model, idx).component isa Memory for idx in loop])

"""
    getloops(model)

Returns idx of nodes that constructs algrebraic loops.
"""
getloops(model::Model) = simplecycles(model.graph)

# LoopBreaker to break the loop
@def_static_system struct LoopBreaker{OP, RO} <: AbstractStaticSystem
    input::Nothing = nothing
    output::OP
    readout::RO
end

"""
    breakloop!(model, loop, breakpoint=length(loop))

Breaks the algebraic `loop` of `model`. The `loop` of the `model` is broken by inserting a `Memory` at the `breakpoint`
of loop.
"""
function breakloop!(model::Model, loop, breakpoint=length(loop))
    nftidx = findfirst(idx -> !isfeedthrough(getnode(model, idx).component), loop)
    nftidx === nothing || (breakpoint = nftidx)

    # Delete the branch at the breakpoint.
    srcnode = getnode(model, loop[breakpoint])
    if breakpoint == length(loop)
        dstnode = getnode(model, loop[1])
    else
        dstnode = getnode(model, loop[(breakpoint + 1)])
    end
    branch = getbranch(model, srcnode.idx => dstnode.idx)

    # Construct the loopbreaker.
    if nftidx === nothing
        nodefuncs = wrap(model, loop)
        ff = feedforward(nodefuncs, breakpoint)
        n = length(srcnode.component.output)
        breaker = LoopBreaker(readout = (u,t) -> findroot(ff, n, t), output=Outport(n))
    else
        component = srcnode.component
        n = length(component.output)
        breaker = LoopBreaker(readout = (u,t)->component.readout(component.state, nothing, t), output=Outport(n))
    end
    # newidx = length(model.nodes) + 1
    newnode = addnode!(model, breaker)

    # Delete the branch at the breakpoint
    deletebranch!(model, branch.nodepair)

    # Connect the loopbreker to the loop at the breakpoint.
    addbranch!(model, newnode.idx => dstnode.idx, branch.indexpair)
    return newnode
end

function wrap(model, loop)
    graph = model.graph
    map(loop) do idx
        node = getnode(model, idx)
        innbrs = filter(i -> i ∉ loop, inneighbors(graph, idx))
        outnbrs = filter(i -> i ∉ loop, outneighbors(graph, idx))
        if isempty(innbrs) && isempty(outnbrs)
            zero_in_zero_out(node)
        elseif isempty(innbrs) && !isempty(outnbrs)
            zero_in_nonzero_out(node, getoutmask(model, node, loop))
        elseif !isempty(innbrs) && isempty(outnbrs)
            nonzero_in_zero_out(node, getinmask(model, node, loop))
        else
            nonzero_in_nonzero_out(node, getinmask(model, node, loop), getoutmask(model, node, loop))
        end
    end
end

function zero_in_zero_out(node)
    component = node.component
    function func(ut)
        u, t = ut
        out = [_computeoutput(component, u, t)...]
        out, t
    end
end

function zero_in_nonzero_out(node, outmask)
    component = node.component
    function func(ut)
        u, t = ut
        out = [_computeoutput(component, u, t)...]
        out[outmask], t
    end
end

function nonzero_in_zero_out(node, inmask)
    component = node.component
    nin = length(inmask)
    function func(ut)
        u, t = ut
        uu = zeros(nin)
        uu[inmask] .= readbuffer(component.input, inmask)
        uu[.!inmask] .= u
        out = [_computeoutput(component, uu, t)...]
        out, t
    end
end

function nonzero_in_nonzero_out(node, inmask, outmask)
    component = node.component
    nin = length(inmask)
    function func(ut)
        u, t = ut
        uu = zeros(nin)
        uu[inmask] .= readbuffer(component.input, inmask)
        uu[.!inmask] .= u
        out = [_computeoutput(component, uu, t)...]
        out[outmask]
        out, t
    end
end

function getinmask(model, node, loop)
    idx = node.idx
    inmask = falses(length(node.component.input))
    for nidx in filter(n -> n ∉ loop, inneighbors(model.graph, idx)) # Not-in-loop inneighbors
        k = getbranch(model, nidx => idx).indexpair.second
        if length(k) == 1
            inmask[k] = true
        else
            inmask[k] .= trues(length(k))
        end
    end
    inmask
end

function getoutmask(model, node, loop)
    idx = node.idx
    outmask = falses(length(node.component.output))
    for nidx in filter(n -> n ∈ loop, outneighbors(model.graph, idx)) # In-loop outneighbors
        k = getbranch(model, idx => nidx).indexpair.first
        if length(k) == 1
            outmask[k] = true
        else
            outmask[k] .= trues(length(k))
        end
    end
    outmask
end

readbuffer(input, inmask) = map(pin -> read(pin.link.buffer), input[inmask])
_computeoutput(comp::AbstractStaticSystem, u, t) = comp.readout(u, t)
_computeoutput(comp::AbstractDynamicSystem, u, t) = comp.readout(comp.state, map(uu -> t -> uu, u), t)

function feedforward(nodefuncs, breakpoint=length(nodefuncs))
    (u, t) -> ∘(reverse(circshift(nodefuncs, -breakpoint))...)((u, t))[1] - u
end

function findroot(ff, n, t)
    sol = nlsolve((dx, x) -> (dx .= ff(x, t)), rand(n))
    sol.zero
end

function isfeedthrough(component)
    try
        out = typeof(component) <: AbstractStaticSystem ?
            component.readout(nothing, 0.) : component.readout(component.state, nothing, 0.)
        return false
    catch ex
        return true
    end
end

# Check if components of nodes of the models has unbound pins. In case there are any unbound pins,
# the simulation is got stuck since the data flow through an unbound pin is not possible.
checknodeports(model) = foreach(node -> checkports(node.component), model.nodes)
function checkports(comp::T) where T
    if hasfield(T, :input)
        idx = unboundpins(comp.input)
        isempty(idx) || error("Input port of $comp has unbound pins at index $idx")
    end
    if hasfield(T, :output)
        idx = unboundpins(comp.output)
        isempty(idx) || error("Output port of $comp has unbound pins at index $idx")
    end
end
unboundpins(port::AbstractPort) = findall(.!isbound.(port))
unboundpins(port::Nothing) = Int[]

# Checks if all the channels the links in the model is open. If a link is not open, than
# it is not possible to bind a task that reads and writes data from the channel.
function checkchannels(model)
    # Check branch links
    for branch in model.branches
        for link in branch.links
            isopen(link) || refresh!(link)
        end
    end
    # Check taskmanager links
    for pin in model.taskmanager.triggerport
        link = only(pin.links)
        isopen(link) || refresh!(link)
    end
    for pin in model.taskmanager.handshakeport
        link = pin.link
        isopen(link) || refresh!(link)
    end
end

##### Model initialization

"""
    $(SIGNATURES)

Initializes `model` by launching component task for each of the component of `model`. The pairs component and component tasks are recordedin the task manager of the `model`. The `model` clock is [`set!`](@ref) and the files of [`Writer`](@ref) are openned.
"""
function initialize!(model::Model)
    # NOTE: Tasks to make the components be triggerable are launched here.
    # The important point here is that the simulation should be cancelled if an error is thrown in any of the tasks
    # launched here. This is done by binding the task to the chnnel of the trigger link of the component. Hrence the
    # lifetime of the channel of the link connecting the component to the taskmanger is determined by the lifetime of
    # the task launched for the component. To cancel the simulation and report the stacktrace the task is `fetch`ed.
    bind!(model)

    # Turn on clock model clock if it is running.
    initclock!(model)

    # Clean the model
    clean!(model)

    # Open the files, GUI's for sink components.
    opensinks!(model)
end

# Find the link connecting `component` to `taskmanager`.
function whichlink(taskmanager, component)
    tpin = component.trigger
    tport = taskmanager.triggerport
    # NOTE: `component` must be connected to `taskmanager` by a single link which is checked by `only`
    # `outpin.links` must have just a single link which checked by `only`
    outpin = filter(pin -> isconnected(pin, tpin), tport) |> only
    outpin.links |> only
end

# Bind tasks to connections
function bind!(model::Model)
    taskmanager = model.taskmanager
    pairs = taskmanager.pairs
    nodes = model.nodes
    for node in nodes
        component = node.component
        link = whichlink(taskmanager, component)  # Link connecting the component to taskmanager.
        task = launch(component)    # Task launched to make `componnent` be triggerable.
        bind(link.channel, task)    # Bind the task to the channel of the link.
        pairs[component] = task
    end
    model
end

# Open AbstractSink components
function opensinks!(model::Model)
    foreach(node -> open(node.component), filter(node->isa(node.component, AbstractSink), model.nodes))
    model
end

# Initialze model clock
function initclock!(model::Model)
    if isoutoftime(model.clock)
        msg = "Model clock is out of time. Its current time $(model.clock.t) should be less than its final time "
        msg *= "$(model.clock.tf). Resettting the model clock to its defaults."
        @warn msg
        set!(model.clock)
    end
    isrunning(model.clock) || set!(model.clock)
    model
end

"""
    $(SIGNATURES)

Cleans the buffers of the links of the connections, internal buffers of components, current time of dynamical systems.
"""
function clean!(model::Model)
    cleanbranches!(model)
    cleannodes!(model)
    resetdynamicalsystems!(model)
    model
end

function cleanbranches!(model)
    foreach(branch -> foreach(link -> clean!(link), branch.links), model.branches)
    model
end

function cleannodes!(model)
    foreach(node -> cleancomponent!(node.component), model.nodes)
    model
end

function cleancomponent!(comp::T) where T
    # Clean all internal buffers
    foreach(name -> (field = getfield(comp, name); field isa Buffer && clean!(field)), fieldnames(T))
    # If comp is a Writer, renew the file.
    comp isa Writer && (close(comp); open(comp, "w"); close(comp))
    comp
end

function resetdynamicalsystems!(model)
    t = model.clock.ti
    for comp in filter(comp -> comp isa AbstractDynamicSystem, getfield.(model.nodes, :component))
        comp.t != t && (comp.t = t)
        reinit!(comp.integrator)
        comp.input === nothing || (interp = comp.interpolant; clean!(interp.timebuf); clean!(interp.databuf))
    end
    model
end

##### Model running
# Copy-paste loop body. See `run!(model, withbar)`.
# NOTE: We first trigger the component, Then the tasks fo the `taskmanager` is checked. If an error is thrown in one
# of the tasks, the simulation is cancelled and stacktrace is printed reporting the error. In order to ensure the
# time synchronization between the components of the model, `handshakeport` of the taskmanger is read. When all the
# components take step succesfully, then the simulation goes with the next step after calling the callbacks of the
# components.
# Note we first check the tasks of the taskmanager and then read the `handshakeport` of the taskmanager. Otherwise,
# the simulation gets stuck without printing the stacktrace if an error occurs in one of the tasks of the taskmanager.
@def loopbody begin
    put!(triggerport, fill(t, ncomponents))
    checktaskmanager(taskmanager)
    all(take!(handshakeport)) || @warn "Taking step could not be approved."
    applycallbacks(model)
end

"""
    run!(model::Model, withbar::Bool=true)

Runs the `model` by triggering the components of the `model`. This triggering is done by generating clock tick using the model clock `model.clock`. Triggering starts with initial time of model clock, goes on with a step size of the sampling period of the model clock, and finishes at the finishing time of the model clock. If `withbar` is `true`, a progress bar indicating the simulation status is displayed on the console.

!!! warning
    The `model` must first be initialized to be run. See also: [`initialize!`](@ref).
"""
function run!(model::Model, withbar::Bool=true)
    taskmanager = model.taskmanager
    triggerport, handshakeport = taskmanager.triggerport, taskmanager.handshakeport
    ncomponents = length(model.nodes)
    clock = model.clock
    withbar ? (@showprogress clock.dt for t in clock @loopbody end) : (for t in clock @loopbody end)
    model
end

# ##### Model termination
# """
#     release(model::Model)

# Releaes the each component of `model`, i.e., the input and output bus of each component is released.
# """
# release(model::Model) = foreach(release, model.nodes)

"""
    terminate!(model::Model)

Terminates `model` by terminating all the components of the `model`, i.e., the components tasks in the task manager of the `model` is terminated.
"""
function terminate!(model::Model)
    taskmanager = model.taskmanager
    tasks = unwrap(collect(values(taskmanager.pairs)), Task, depth=length(taskmanager.pairs))
    any(istaskstarted.(tasks)) && put!(taskmanager.triggerport, fill(NaN, length(model.nodes)))
    isrunning(model.clock) && stop!(model.clock)
    model
end

##### Model simulation.
function _simulate(sim::Simulation, reportsim::Bool, withbar::Bool, breakpoints::Vector{Int})
    model = sim.model
    @siminfo "Started simulation..."
    sim.state = :running

    @siminfo "Inspecting model..."
    inspect!(model, breakpoints)
    @siminfo "Done."

    @siminfo "Initializing the model..."
    initialize!(model)
    @siminfo "Done..."

    @siminfo "Running the simulation..."
    run!(model, withbar)
    sim.state = :done
    sim.retcode = :success
    @siminfo "Done..."

    @siminfo "Terminating the simulation..."
    terminate!(model)
    @siminfo "Done."

    reportsim && report(sim)
    return sim
end

"""
    simulate!(model::Model; simdir::String=tempdir(), simprefix::String="Simulation-", simname=string(uuid4()),
        logtofile::Bool=false, loglevel::LogLevel=Logging.Info, reportsim::Bool=false, withbar::Bool=true)

Simulates `model`. `simdir` is the path of the directory into which simulation files are saved. `simprefix` is the prefix of the simulation name `simname`. If `logtofile` is `true`, a log file for the simulation is constructed. `loglevel` determines the logging level. If `reportsim` is `true`, model components are saved into files. If `withbar` is `true`, a progress bar indicating the simualation status is displayed on the console.
"""
function simulate!(model::Model; simdir::String=tempdir(), simprefix::String="Simulation-", simname=string(uuid4()),
    logtofile::Bool=false, loglevel::LogLevel=Logging.Info, reportsim::Bool=false, withbar::Bool=true,
    breakpoints::Vector{Int}=Int[])

    # Construct a Simulation
    sim = Simulation(model, simdir=simdir, simprefix=simprefix, simname=simname)
    sim.logger = logtofile ? SimpleLogger(open(joinpath(sim.path, "simlog.log"), "w+"), loglevel) : ConsoleLogger(stderr, loglevel)

    # Simualate the modoel
    with_logger(sim.logger) do
        _simulate(sim, reportsim, withbar, breakpoints)
    end
    logtofile && flush(sim.logger.stream)  # Close logger file stream.
    return sim
end

"""
    simulate!(model::Model, t0::Real, dt::Real, tf::Real; kwargs...)

Simulates the `model` starting from the initial time `t0` until the final time `tf` with the sampling interval of `tf`. For `kwargs` are

* `logtofile::Bool`: If `true`, a log file is contructed logging each step of the simulation.
* `reportsim::Bool`: If `true`, `model` components are written files after the simulation. When this file is read back, the model components can be consructed back with their status at the end of the simulation.
* `simdir::String`: The path of the directory in which simulation file are recorded.
"""
function simulate!(model::Model, t0::Real, dt::Real, tf::Real; kwargs...)
    set!(model.clock, t0, dt, tf)
    simulate!(model; kwargs...)
end

#### Troubleshooting
"""
    troubleshoot(model)

Prints the exceptions of the tasks that are failed during the simulation of `model`.
"""
function troubleshoot(model::Model)
    fails = filter(pair -> istaskfailed(pair.second), model.taskmanager.pairs)
    if isempty(fails)
        @info "No failed tasks in $model."
    else
        for (comp, task) in fails
            println("", comp)
            @error task.exception
        end
    end
end

##### Plotting signal flow of the model
"""
    signalflow(model, args...; kwargs...)

Plots the signal flow of `model`. `args` and `kwargs` are passed into [`gplot`](https://github.com/JuliaGraphs/GraphPlot.jl) function.
"""
signalflow(model::Model, args...; kwargs...) =
    gplot(model.graph, args...; nodelabel=[node.label for node in model.nodes], kwargs...)

##### @model macro

function check_macro_syntax(name, ex)
    name isa Symbol || error("Invalid usage of @model")
    ex isa Expr && ex.head == :block || error("Invalid usage of @model")
end

function check_block_syntax(node_expr, branch_expr)
    #-------------------  Node expression check ---------------
    # Check syntax the following syntax
    # @nodes begin
    #   label1 = Component1()
    #   label2 = Component2()
    #       ⋮
    # end
    (
        node_expr isa Expr &&
        node_expr.head === :(macrocall) &&
        node_expr.args[1] === Symbol("@nodes")
    ) ||  error("Invalid usage of @nodes")
    node_block = node_expr.args[3]
    (
        node_block.head === :block &&
        all([ex.head === :(=) for ex in filter(arg -> isa(arg, Expr), node_block.args)])
    ) || error("Invalid usage of @nodes")

    #---------------------  Branch expression check --------------
    # Check syntax the following syntax
    # @branches begin
    #   src1[srcidx1] => dst1[dstidx1]
    #   src2[srcidx2] => dst2[dstidx2]
    #       ⋮
    # end
    (
        branch_expr isa Expr &&
        branch_expr.head === :(macrocall) &&
        branch_expr.args[1] === Symbol("@branches")
    ) || error("Invalid usage of @branches")
    branch_block = branch_expr.args[3]
    (
        branch_block.head === :block &&
        all([ex.head === :call && ex.args[1] == :(=>) for ex in filter(arg -> isa(arg, Expr), branch_block.args)])
    ) || error("Invalid usage of @branches")
end

"""
    @defmodel name ex

Construts a model. The expected syntax is.
```
    @defmodel mymodel begin
        @nodes begin
            label1 = Component1()
            label2 = Component1()
                ⋮
        end
        @branches begin
            src1 => dst1
            src2 => dst2
                ⋮
        end
    end
```
Here `@nodes` and `@branches` blocks adefine the nodes and branches of the model, respectively.
"""
macro defmodel(name, ex)
    # Check syntax
    check_macro_syntax(name, ex)
    node_expr = ex.args[2]
    branch_expr = ex.args[4]
    check_block_syntax(node_expr, branch_expr)

    # Extract nodes info
    node_block = node_expr.args[3]
    node_labels = [expr.args[1] for expr in node_block.args if expr isa Expr]
    node_components = [expr.args[2] for expr in node_block.args if expr isa Expr]

    # Extract branches info
    branch_block = branch_expr.args[3]
    lhs = [expr.args[2] for expr in filter(ex -> isa(ex, Expr), branch_block.args)]
    rhs = [expr.args[3] for expr in filter(ex -> isa(ex, Expr), branch_block.args)]
    quote
        # Construct model
        $name = Model()

        # Add nodes to model
        for (node_label, node_component) in zip($node_labels, $node_components)
            addnode!($name, eval(node_component), label=node_label)
        end

        # Add braches to model
        for (src, dst) in zip($lhs, $rhs)
            if src isa Symbol && dst isa Symbol
                addbranch!($name, src => dst)
            elseif src isa Expr && dst isa Expr # src and dst has index.
                if src.args[2] isa Expr && dst.args[2] isa Expr
                    # array or range index.
                    addbranch!($name, src.args[1] => dst.args[1], eval(src.args[2]) => eval(dst.args[2]))
                else
                    # integer index
                    addbranch!($name, src.args[1] => dst.args[1], src.args[2] => dst.args[2])
                end
            else
                error("Ambbiuos connection. Specify the indexes explicitely.")
            end
        end
    end |> esc
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 2cb37ccd-f298-404a-8dea-0815bcd7fda9
using PlutoTest

# ╔═╡ 4f66ae50-b0e9-11eb-0350-3535845809d1
md"""
# [Day 19: Medicine for Rudolph](https://adventofcode.com/2015/day/19)

Rudolph the Red-Nosed Reindeer is sick! His nose isn't shining very brightly, and he needs medicine.

Red-Nosed Reindeer biology isn't similar to regular reindeer biology; Rudolph is going to need custom-made medicine. Unfortunately, Red-Nosed Reindeer chemistry isn't similar to regular reindeer chemistry, either.

The North Pole is equipped with a Red-Nosed Reindeer nuclear fusion/fission plant, capable of constructing any Red-Nosed Reindeer molecule you need. It works by starting with some input molecule and then doing a series of replacements, one per step, until it has the right molecule.

However, the machine has to be calibrated before it can be used. Calibration involves determining the number of molecules that can be generated in one step from a given starting point.

For example, imagine a simpler machine that supports only the following replacements:

```
H => HO
H => OH
O => HH

```

Given the replacements above and starting with `HOH`, the following molecules could be generated:

-   `HOOH` (via `H => HO` on the first `H`).
-   `HOHO` (via `H => HO` on the second `H`).
-   `OHOH` (via `H => OH` on the first `H`).
-   `HOOH` (via `H => OH` on the second `H`).
-   `HHHH` (via `O => HH`).

So, in the example above, there are `4` distinct molecules (not five, because `HOOH` appears twice) after one replacement from `HOH`. Santa's favorite molecule, `HOHOHO`, can become `7` distinct molecules (over nine replacements: six from `H`, and three from `O`).

The machine replaces without regard for the surrounding characters. For example, given the string `H2O`, the transition `H => OO` would result in `OO2O`.

Your puzzle input describes all of the possible replacements and, at the bottom, the medicine molecule for which you need to calibrate the machine. **How many distinct molecules can be created after all the different ways you can do one replacement on the medicine molecule?**
"""

# ╔═╡ 1a84a9cf-60bb-423b-b802-b35763e96032
puzzle_input = "Al => ThF
Al => ThRnFAr
B => BCa
B => TiB
B => TiRnFAr
Ca => CaCa
Ca => PB
Ca => PRnFAr
Ca => SiRnFYFAr
Ca => SiRnMgAr
Ca => SiTh
F => CaF
F => PMg
F => SiAl
H => CRnAlAr
H => CRnFYFYFAr
H => CRnFYMgAr
H => CRnMgYFAr
H => HCa
H => NRnFYFAr
H => NRnMgAr
H => NTh
H => OB
H => ORnFAr
Mg => BF
Mg => TiMg
N => CRnFAr
N => HSi
O => CRnFYFAr
O => CRnMgAr
O => HP
O => NRnFAr
O => OTi
P => CaP
P => PTi
P => SiRnFAr
Si => CaSi
Th => ThCa
Ti => BP
Ti => TiTi
e => HF
e => NAl
e => OMg

ORnPBPMgArCaCaCaSiThCaCaSiThCaCaPBSiRnFArRnFArCaCaSiThCaCaSiThCaCaCaCaCaCaSiRnFYFArSiRnMgArCaSiRnPTiTiBFYPBFArSiRnCaSiRnTiRnFArSiAlArPTiBPTiRnCaSiAlArCaPTiTiBPMgYFArPTiRnFArSiRnCaCaFArRnCaFArCaSiRnSiRnMgArFYCaSiRnMgArCaCaSiThPRnFArPBCaSiRnMgArCaCaSiThCaSiRnTiMgArFArSiThSiThCaCaSiRnMgArCaCaSiRnFArTiBPTiRnCaSiAlArCaPTiRnFArPBPBCaCaSiThCaPBSiThPRnFArSiThCaSiThCaSiThCaPTiBSiRnFYFArCaCaPRnFArPBCaCaPBSiRnTiRnFArCaPRnFArSiRnCaCaCaSiThCaRnCaFArYCaSiRnFArBCaCaCaSiThFArPBFArCaSiRnFArRnCaCaCaFArSiRnFArTiRnPMgArF"

# ╔═╡ 03e179fb-d380-4d8d-919c-7dc940533876
example_input = """
H => HO
H => OH
O => HH
""" |> strip

# ╔═╡ 43daa5a3-6c7f-4afd-9bba-bb514bf41e83
function parse1(input)
	d = Dict()
	for line in split(input, "\n")
		parts = split(line, " => ")

		key = parts[1]
		if haskey(d, key)
			push!(d[key], parts[2])
		else
			d[key] = [parts[2]]
		end
	end
	d
end

# ╔═╡ fcad707b-3039-4753-bfe1-7db9a5846c48
parse1(example_input)

# ╔═╡ d393d954-9dc7-4aa3-a219-7a79458babe1
map(f -> (f.match, f.offset), collect(eachmatch(Regex("H"), "HOH")))

# ╔═╡ 19be711d-ca4c-478e-8422-0ece0efefc7e
function molcombinations(mol, c)
	vals = []
	for key in keys(c)

		for m in eachmatch(Regex(key), mol)

			for val in c[key]

				n = collect(mol)
				n = collect(val)
				push!(vals,
					(
					mol[m.offset:m.offset+length(key)-1],

					)
				)

			end
		end

		#push!(n, c[key] * length(collect(eachmatch(Regex(key), mol))))
	end
	vals
end

# ╔═╡ 0662d133-08c2-4f14-857e-a77a8cb2ba43
molcombinations("HOH", parse1(example_input))

# ╔═╡ b30427b8-b806-4197-9879-b351587b5fb6
@assert length(molcombinations("HOH", example_input)) == 4

# ╔═╡ d29f7880-d4c0-40a7-a92a-c45883c3ae6f
@assert length(molcombinations("HOHOHO", example_input)) == 7

# ╔═╡ Cell order:
# ╠═2cb37ccd-f298-404a-8dea-0815bcd7fda9
# ╟─4f66ae50-b0e9-11eb-0350-3535845809d1
# ╟─1a84a9cf-60bb-423b-b802-b35763e96032
# ╟─03e179fb-d380-4d8d-919c-7dc940533876
# ╠═43daa5a3-6c7f-4afd-9bba-bb514bf41e83
# ╠═fcad707b-3039-4753-bfe1-7db9a5846c48
# ╠═d393d954-9dc7-4aa3-a219-7a79458babe1
# ╠═19be711d-ca4c-478e-8422-0ece0efefc7e
# ╠═0662d133-08c2-4f14-857e-a77a8cb2ba43
# ╠═b30427b8-b806-4197-9879-b351587b5fb6
# ╠═d29f7880-d4c0-40a7-a92a-c45883c3ae6f