text:
module FFmpegPipe

using ImageMagick, FileIO, Images
# # alternatively: QuartsImageIO

export openvideo, readframe, writeframe

const ffmpeg = "ffmpeg" # name of executable

"""
Open a movie file using ffmpeg's image2pipe interface.
"""
openvideo(filename::String; kwargs...) = openvideo(filename, :r; kwargs...)
openvideo(filename::String, mode::Symbol; kwargs...) = openvideo(filename, Val{mode}(); kwargs...)
openvideo(filename::String, mode::Union{Char,String}; kwargs...) = openvideo(filename, Symbol(mode); kwargs...)

function addkwargs!(input_options, output_options, kwargs)
    for (namein, value) in kwargs
        m = match(r"^(.*)_(in|out)", string(namein))
        if m ≢ nothing
            name, direction = m.captures
            if direction == "in"
                input_options[Symbol(name)] = value
            else
                output_options[Symbol(name)] = value
            end
        else
            output_options[namein] = value
        end
    end
end

"""
`openvideo(file, :r)` opens movie file for reading
"""

function openvideo(filename::String, ::Val{:r}; loglevel_out = "fatal", kwargs...)
    input_options = Dict{Symbol, Any}()
    output_options = Dict{Symbol, Any}(:loglevel => loglevel_out)
    addkwargs!(input_options, output_options, kwargs)
    before = Iterators.flatten(("-$k", v) for (k,v) in pairs(input_options))
    after = Iterators.flatten(("-$k", v) for (k,v) in pairs(output_options))
    cmd = `$ffmpeg -nostats $before -i $filename $after -f image2pipe -vcodec png -compression_level 0 -`
    open(cmd)
end

"""
`openvideo(file, :w)` opens movie file for writing
"""
function openvideo(filename::String, ::Val{:w}; r_out=24, q_out=3, vcodec_out="h264", loglevel_out="fatal", kwargs...)
    input_options = Dict{Symbol, Any}()
    output_options = Dict{Symbol, Any}(:r => r_out, :q => q_out, :vcodec => vcodec_out, :loglevel => loglevel_out)
    addkwargs!(input_options, output_options, kwargs)
    before = Iterators.flatten(("-$k", v) for (k,v) in input_options)
    after = Iterators.flatten(("-$k", v) for (k,v) in output_options)
    cmd = `$ffmpeg -nostats -f image2pipe -vcodec png $before -i - $after $filename`
    open(cmd, "w")
end

function openvideo(f::Function, filename::String, mode; kwargs...)
    io = openvideo(filename, mode; kwargs...)
    try
        f(io)
    finally
        close(io)
    end
end

# Read up to the end of a .png image, one chunk at a time,
# while parsing for the "IEND" chunk.
# (Since FFmpeg sends a stream of concatenated images,
# we can't just read until eof.)
function readpngdata(io)
    blk = 65536;
    a = Array{UInt8}(undef, blk)
    readbytes!(io, a, 8)
    if view(a, 1:8) != magic(format"PNG")
        error("Bad magic.")
    end
    n = 8
    while !eof(io)
        if length(a)<n+12
            resize!(a, length(a)+blk)
        end
        readbytes!(io, view(a, n+1:n+12), 12)
        m = 0
        for i=1:4
            m = m<<8 + a[n+i]
        end
        chunktype = view(a, n+5:n+8)
        n=n+12
        if chunktype == codeunits("IEND")
            break
        end
        if length(a)<n+m
            resize!(a, max(length(a)+blk, n+m+12))
        end
        readbytes!(io, view(a, n+1:n+m), m)
        n = n+m
    end
    resize!(a,n)
    return a
end

"Read a frame from a video"
readframe(io::IO) = ImageMagick.load_(readpngdata(io))

"Write a frame to a video"
writeframe(io::IO, img) = show(io, MIME("image/png"), img)
# Bypass re-scaling:
writeframe(io::IO, img::Array{T,2} where T<:Colorant) = save(Stream(format"PNG", io), img)

end # module


text:
@testset "Update boundaries" begin for b1index in 1:16, b2index in 9:16 begin
    q1 = floor((b1index - 1)/4)
    q2 = floor((b2index - 1)/4)

    if q1 != q2
      #b1index = 1
      #b2index = 9

      print_with_color(:yellow, "------------------Testing ($(b1index), $(b2index))\n")
      # add boundaries to single element
      bb = GeometryTypes.SimpleRectangle(2.0, 3.0, 4.0, 4.0)
      mesh = QuadTreeMeshes.QuadTreeMesh(bb)
      qt = mesh.quadtree
      qtEl = qt.elements[1]

      # add boundary vertices
      push!(qt.vertices, QuadTreeMeshes.Point(get_boundary_coordinates_from_index(b1index, qtEl, qt)))
      vb1Index = length(qt.vertices)
      push!(qt.vertices, QuadTreeMeshes.Point(get_boundary_coordinates_from_index(b2index, qtEl, qt)))
      vb2Index = length(qt.vertices)

      QuadTreeMeshes.triangulate_boundary_leave(mesh, 1, vb1Index, vb2Index)
      mesh_element = get(qt.values[1])

      # subdivide new element
      QuadTreeMeshes.subdivide!(qt, 1, QuadTreeMeshes.OnChildrenCreated)

      if plot==true
        filename = "triangulate_bndy_leave.html"
        Plots.plot(qt)
        Plots.plot!(mesh, boundaries_only=true)
        Plots.savefig(filename)
      end

      b1, b2 = qt.vertices[vb1Index], qt.vertices[vb2Index]
      b1b2s = GeometryTypes.LineSegment(b1, b2)

      function check_child_element(elIndex::QuadTreeMeshes.ElIndex)
        check_leave_intersection(qt, elIndex, b1b2s, Nullable{Int64}())
        # check that there a no boundary intersections left
        @test(length(get(qt.values[elIndex]).boundaries) == 0)
      end

      # check all child elements
      check_child_element(get(qtEl.northWest))
      check_child_element(get(qtEl.northEast))
      check_child_element(get(qtEl.southWest))
      check_child_element(get(qtEl.southEast))
    end
  end
end
end
end

ext: jl


text:
using Printf

@testset "printf" begin
    for (fmt, val) in (("%7.2f", "   1.23"),
                   ("%-7.2f", "1.23   "),
                   ("%07.2f", "0001.23"),
                   ("%.0f", "1"),
                   ("%#.0f", "1."),
                   ("%.4e", "1.2345e+00"),
                   ("%.4E", "1.2345E+00"),
                   ("%.2a", "0x1.3cp+0"),
                   ("%.2A", "0X1.3CP+0")),
        num in (Float128(1.2345),)
        @test @eval(@sprintf($fmt, $num) == $val)
    end
    @test (@sprintf "%f" Float128(Inf)) == "Inf"
    @test (@sprintf "%f" Float128(NaN)) == "NaN"
    @test (@sprintf "%.0e" Float128(3e142)) == "3e+142"
    @test (@sprintf "%#.0e" Float128(3e142)) == "3.e+142"

    @test (@sprintf "%.0e" Float128(big"3e1042")) == "3e+1042"

    for (val, res) in ((Float128(12345678.), "1.23457e+07"),
                   (Float128(1234567.8), "1.23457e+06"),
                   (Float128(123456.78), "123457"),
                   (Float128(12345.678), "12345.7"))
        @test (@sprintf("%.6g", val) == res)
    end

    for (fmt, val) in (("%10.5g", "     123.4"),
                       ("%+10.5g", "    +123.4"),
                       ("% 10.5g","     123.4"),
                       ("%#10.5g", "    123.40"),
                       ("%-10.5g", "123.4     "),
                       ("%-+10.5g", "+123.4    "),
                       ("%010.5g", "00000123.4")),
        num in (Float128(123.4),)
        @test @eval(@sprintf($fmt, $num) == $val)
    end
    @test( @sprintf( "%10.5g", Float128(-123.4) ) == "    -123.4")
    @test( @sprintf( "%010.5g", Float128(-123.4) ) == "-0000123.4")
    @test( @sprintf( "%.6g", Float128(12340000.0) ) == "1.234e+07")
    @test( @sprintf( "%#.6g", Float128(12340000.0)) == "1.23400e+07")

    @test (@sprintf "%a" Float128(1.5)) == "0x1.8p+0"

end


text:
# Contains all helper functions for Tmatrix module

using LinearAlgebra

"""
    Calculate dot product for two vectors.
I made this function because the "dot" function in LinearAlgebra package doesn't work as expected for complex vectors.
Parameters
==========
A, B : two vectors

Returns
=======
scalar dot product
"""
function vector_dot_product(A, B)
    return sum(A .* B)
end

"""
    Create a meshgrid, to be used with numerical integrals.
The meshgrid created is following this rule:
x_grid[dim1,dim2]
dim1 corresponds to x index
dim2 corresponds to y index
"""
function meshgrid(x, y)
    x_grid = repeat(x, 1, length(y))
    y_grid = repeat(y', length(x), 1)
    return x_grid, y_grid
end

"""
    Create meshgrid of θ and ϕ
"""
function meshgrid_θ_ϕ(n_θ_points, n_ϕ_points; min_θ=1e-16, min_ϕ=1e-16, rotationally_symmetric=false)
    θ_1D_array = LinRange(min_θ, π, n_θ_points);
    ϕ_1D_array = LinRange(min_ϕ, 2π, n_ϕ_points);
    if rotationally_symmetric
        θ_array = collect(θ_1D_array)
        ϕ_array = zeros(size(θ_array))
    else
        θ_array, ϕ_array = meshgrid(θ_1D_array, ϕ_1D_array);
    end
    return θ_array, ϕ_array
end

"""
Create a single index from m and n
We fill the index like this:
`
idx = 0
for n = 1:n_max
    for m = -n:n
        global idx
        idx += 1
    end
end
`
"""
function single_index_from_m_n(m::Int, n::Int)
    return n * (n + 1) + m
end

"""
    Get the maximum single index, given the maximum n.
"""
function get_max_single_index_from_n_max(n_max::Int)
    return single_index_from_m_n(n_max, n_max)
end

"""
    Multiply the integrand by the `dS` element, which equals r²sin(θ)
"""
function surface_integrand(
        integrand::AbstractVecOrMat{C}, r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}
    ) where {R <: Real,C <: Complex{R}}
    return integrand .* r_array.^2 .* sin.(θ_array)
end

"""
    Multiply the integrand by the `dS` element, which equals r²sin(θ)
"""
function surface_integrand(
        integrand::AbstractVecOrMat{R}, r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}
    ) where {R <: Real}
    return integrand .* r_array.^2 .* sin.(θ_array)
end

"""
    Multiply the integrand by the `dS` element, which equals r²sin(θ)
"""
function surface_integrand(
        integrand::AbstractVecOrMat{Matrix{C}}, r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}
    ) where {R <: Real,C <: Complex{R}}
    return integrand .* r_array.^2 .* sin.(θ_array)
end

"""
    Multiply the integrand by the `dS` element, which equals r²sin(θ)
"""
function surface_integrand(
        integrand::AbstractVecOrMat{Matrix{R}}, r_array::AbstractVecOrMat{R}, θ_array::AbstractVecOrMat{R}
    ) where {R <: Real,C <: Complex{R}}
    return integrand .* r_array.^2 .* sin.(θ_array)
end

"""
    Save T-matrix to HDF5 file, with fields "Tmatrix_real_CELES_convention" and "Tmatrix_imag_CELES_convention"
"""
function save_Tmatrix_to_HDF5_file(T, HDF5_filename)
    h5write(HDF5_filename, "Tmatrix_real_CELES_convention", Float64.(real(T))) # TODO: can we store BigFloat in HDF5? I converted to `Float64` because I couldn't read the file when I stored BigFloat
    h5write(HDF5_filename, "Tmatrix_imag_CELES_convention", Float64.(imag(T)))
end

# getting indices of T-matrix
"""
    Get 2D array, first and second columns are n and m values, respectively.
m = -n : 1 : +n
n will be the repeated value of input
"""
function get_n_m_array_given_n(n)
    return hcat(repeat([n], 2 * n + 1), -n:n)
end

"""
    Get 2D array, first and second columns are n and m values, respectively.
This is looping `get_n_m_array_given_n` over all n = 1:1:n_max
"""
function get_n_m_array_given_n_max(n_max)
    return vcat(get_n_m_array_given_n.(1:n_max)...)
end

"""
    Get 2D array, first, second and third columns are n, m, and idx values, respectively.
The idx value is just the index of the row.
"""
function get_n_m_idx_array_given_n_max(n_max)
    nm_array = get_n_m_array_given_n_max(n_max)
    return hcat(nm_array, 1:size(nm_array, 1))
end

"""
    Get 4 matrices for m, n, m_, n_, corresponding to rank and order of incident and scattered VSWF, represented by elements of T-matrix
Arguments
---------
    n_max : Int, determine the size of the T-matrix
Returns
-------
    m_matrix, n_matrix, m__matrix, n__matrix : each is a square matrix

"""
function get_m_n_m__n__matrices_for_T_matrix(n_max)
    n_m_idx = get_n_m_idx_array_given_n_max(n_max)
    n_array = n_m_idx[:,1]
    m_array = n_m_idx[:,2]
    idx_array = n_m_idx[:,3]

    idx_max = idx_array[end]
    idx_matrix = repeat(idx_array, 1, idx_max)
    idx__matrix = repeat(idx_array', idx_max, 1)

    n_matrix = n_array[idx_matrix]
    m_matrix = m_array[idx_matrix]

    n__matrix = n_array[idx__matrix]
    m__matrix = m_array[idx__matrix]

    return m_matrix, n_matrix, m__matrix, n__matrix
end

"""
    using for loops, Get 4 matrices for m, n, m_, n_, corresponding to rank and order of incident and scattered VSWF, represented by elements of T-matrix
"""
function get_m_n_m__n__matrices_for_T_matrix_with_forloop(n_max)
    idx_max = Tmatrix.get_max_single_index_from_n_max(n_max)
    n_matrix = zeros(Int, idx_max, idx_max)
    m_matrix = zeros(Int, idx_max, idx_max)
    n__matrix = zeros(Int, idx_max, idx_max)
    m__matrix = zeros(Int, idx_max, idx_max)
    idx = 0;
    for n = 1:n_max
        for m = -n:n
            idx += 1
            idx_ = 0;
            for n_ = 1:n_max
                for m_ = -n_:n_
                    idx_ += 1
                    m_matrix[idx, idx_] = m
                    n_matrix[idx, idx_] = n
                    m__matrix[idx, idx_] = m_
                    n__matrix[idx, idx_] = n_
                end
            end
        end
    end
    return m_matrix, n_matrix, m__matrix, n__matrix
end

"""
    Create m,n,m_,n_ matrices with for loop, and compare the result with `get_m_n_m__n__matrices_for_T_matrix`
"""
function validate_get_m_n_m__n__matrices_for_T_matrix(n_max)

    m_matrix, n_matrix, m__matrix, n__matrix = get_m_n_m__n__matrices_for_T_matrix_with_forloop(n_max)
    _m_matrix, _n_matrix, _m__matrix, _n__matrix = get_m_n_m__n__matrices_for_T_matrix(n_max)

    println(_n_matrix == n_matrix)
    println(_m_matrix == m_matrix)
    println(_n__matrix == n__matrix)
    println(_m__matrix == m__matrix)
end

#######################
# TOBE moved to a new package for trapezoidal integraion, that will be compatible with autodiff
"""
    1D numerical integral using trapezoidal rule
x and y are 1D arrays
"""
function trapz_ELZOUKA(x::AbstractVector{R}, y::AbstractVector{N}) where {R <: Real,N <: Number}
    # TODO: small error if compared with Trapz.trapz
    base = x[2:end] - x[1:end - 1]
    av_height = (y[2:end] + y[1:end - 1]) / 2
    areas = base .* av_height
    total_area = sum(areas)
    return total_area
end

"""
    2D numerical integral using trapezoidal rule
x and y are 1D arrays, z is 2D array
"""
function trapz_ELZOUKA(x::AbstractVector{R}, y::AbstractVector{R}, z::AbstractMatrix{N}) where {R <: Real,N <: Number}
    integrand_wrt_x = trapz_ELZOUKA.(eachcol(repeat(x, 1, size(z, 2))), eachcol(z))
    return trapz_ELZOUKA(y, integrand_wrt_x)
end

#######################
# Calculate orientation-averaged sscattering, extinction and absorption cross sections.
"""
    calculate orientation averaged scattering cross section given a T-matrix and a wavevector.
Input T-matrix can be a complex square matrix or a concatenation of separate real and imag parts of T-matrix
"""
function get_OrentationAv_scattering_CrossSection_from_Tmatrix(T::AbstractMatrix, k1::C) where C <: Complex
    if size(T)[1] == size(T)[2] # if T-matrix is a square matrix, then T-matrix is complex
        return real(2 * pi / k1^2 * sum(T .* conj(T)))
    elseif size(T)[1] == size(T)[2] / 2 # if T-matrix has number of columns double the number of rows, then T-matrix is hcat() of real and imag parts of Tmatrix
        return get_OrentationAv_scattering_CrossSection_from_Tmatrix_SeparateRealImag(T, real(k1), imag(k1))
    end
end

"""
    calculate orientation averaged scattering cross section given a T-matrix, wavelength and optical properties of surrounding.
Input T-matrix can be a complex square matrix or a concatenation of separate real and imag parts of T-matrix
"""
function get_OrentationAv_scattering_CrossSection_from_Tmatrix(T::AbstractMatrix,  wl_or_freq_input::R, input_unit::String, Eps_r_1::C, Mu_r_1::C) where {C <: Complex,R <: Real}
    return get_OrentationAv_scattering_CrossSection_from_Tmatrix(
        T,
        get_WaveVector(
            wl_or_freq_input;
            input_unit=input_unit,
            Eps_r=Eps_r_1,
            Mu_r=Mu_r_1
        )
    )
end

"""
    calculate orientation averaged scattering cross section given a T-matrix and a wavevector. This works for "_SeparateRealImag"
"""
function get_OrentationAv_scattering_CrossSection_from_Tmatrix_SeparateRealImag(T::AbstractMatrix{R}, k1_r::R, k1_i::R) where R <: Real
    T_r, T_i = Tmatrix.separate_real_imag(T)
    T_by_conjT_sum = sum(complex_multiply.(T_r, T_i, T_r, -T_i))
    k1_squared = complex_multiply(k1_r, k1_i, k1_r, k1_i)
    two_pi_over_k1_squared = complex_divide(2 * pi, 0, k1_squared[1], k1_squared[2])
    return complex_multiply(two_pi_over_k1_squared[1], two_pi_over_k1_squared[2], T_by_conjT_sum[1], T_by_conjT_sum[2])[1]
end

"""
    calculate orientation averaged extinction cross section given a T-matrix and a wavevector.
"""
function get_OrentationAv_extinction_CrossSection_from_Tmatrix(T::AbstractMatrix, k1::C) where C <: Complex
    if size(T)[1] == size(T)[2] # if T-matrix is a square matrix, then T-matrix is complex
        return real(-2 * pi / k1^2 * tr(real(T)))
    elseif size(T)[1] == size(T)[2] / 2 # if T-matrix has number of columns double the number of rows, then T-matrix is hcat() of real and imag parts of Tmatrix
        return get_OrentationAv_extinction_CrossSection_from_Tmatrix_SeparateRealImag(T, real(k1), imag(k1))
    end
end

"""
    calculate orientation averaged extinction cross section given a T-matrix and a wavevector. This works for "_SeparateRealImag"
"""
function get_OrentationAv_extinction_CrossSection_from_Tmatrix_SeparateRealImag(T::AbstractMatrix{R}, k1_r::R, k1_i::R) where R <: Real
    T_r, T_i = Tmatrix.separate_real_imag(T)
    k1_squared = complex_multiply(k1_r, k1_i, k1_r, k1_i)
    two_pi_over_k1_squared = complex_divide(2 * pi, 0, k1_squared[1], k1_squared[2])
    return -two_pi_over_k1_squared[1] * tr(T_r)
end

"""
    calculate orientation averaged absorption cross section given a T-matrix and a wavevector.
"""
function get_OrentationAv_absorption_CrossSection_from_Tmatrix(T::AbstractMatrix, k1::C) where C <: Complex
    return get_OrentationAv_extinction_CrossSection_from_Tmatrix(T, k1) - get_OrentationAv_scattering_CrossSection_from_Tmatrix(T, k1)
end

"""
    calculate orientation averaged absorption cross section given a T-matrix and a wavevector. This works for "_SeparateRealImag"
"""
function get_OrentationAv_absorption_CrossSection_from_Tmatrix_SeparateRealImag(T::AbstractMatrix{R}, k1_r::R, k1_i::R) where R <: Real
    return get_OrentationAv_extinction_CrossSection_from_Tmatrix_SeparateRealImag(T, k1_r, k1_i) - get_OrentationAv_scattering_CrossSection_from_Tmatrix_SeparateRealImag(T, k1_r, k1_i)
end

"""
    calculate orientation averaged emissivity cross section given a T-matrix and a wavevector.
`particle_surface_area` has to be in units consistent with wavevector units. If waveevctor unit is per m, then the particle surface area has to be in m^2
"""
function get_OrentationAv_emissivity_from_Tmatrix(T::AbstractMatrix, k1::C, particle_surface_area::R) where {R <: Real,C <: Complex}
    return get_OrentationAv_absorption_CrossSection_from_Tmatrix(T, k1) * 4 / particle_surface_area
end

"""
    calculate orientation averaged emissivity cross section given a T-matrix and a wavevector. This works for "_SeparateRealImag"
`particle_surface_area` has to be in units consistent with wavevector units. If waveevctor unit is per m, then the particle surface area has to be in m^2
"""
function get_OrentationAv_emissivity_from_Tmatrix_SeparateRealImag(T::AbstractMatrix{R}, k1_r::R, k1_i::R, particle_surface_area::R) where R <: Real
    return get_OrentationAv_absorption_CrossSection_from_Tmatrix_SeparateRealImag(T, k1_r, k1_i) * 4 / particle_surface_area
end

"""
    calculate orientation averaged emissivity cross section given a T-matrix, wavelength and optical properties of surrounding. This works for "_SeparateRealImag"
`particle_surface_area` has to be in units consistent with wavevector units. If waveevctor unit is per m, then the particle surface area has to be in m^2
"""
function get_OrentationAv_emissivity_from_Tmatrix_SeparateRealImag(T::AbstractMatrix{R}, wl_or_freq_input::R, input_unit::String, Eps_r_r_1::R, Eps_r_i_1::R, Mu_r_r_1::R, Mu_r_i_1::R, particle_surface_area::R) where R <: Real
    k1_complex = get_WaveVector(
            wl_or_freq_input;
            input_unit=input_unit,
            Eps_r=Complex(Eps_r_r_1, Eps_r_i_1),
            Mu_r=Complex(Mu_r_r_1, Mu_r_i_1)
    )
    return get_OrentationAv_absorption_CrossSection_from_Tmatrix_SeparateRealImag(T, k1_r, k1_i) * 4 / particle_surface_area
end


text:
using Documenter, QuantumGlue

makedocs(;
    modules=[QuantumGlue],
    format=:html,
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jlapeyre/QuantumGlue.jl/blob/{commit}{path}#L{line}",
    sitename="QuantumGlue.jl",
    authors="John Lapeyre",
    assets=[],
)

deploydocs(;
    repo="github.com/jlapeyre/QuantumGlue.jl",
    target="build",
    julia="1.0",
    deps=nothing,
    make=nothing,
)


text:
tab(n::Int) = " "^n

Base.@kwdef struct Color
    literal::Symbol = :light_blue
    type::Symbol = :light_green
    string::Symbol = :yellow
    comment::Symbol = :light_black
    kw::Symbol = :light_magenta
    fn::Symbol = :light_blue
end

Base.@kwdef mutable struct PrintState
    line_indent::Int = 0
    content_indent::Int = line_indent
    color::Symbol = :normal
end

"""
    sprint_expr(ex; context=nothing)

Print given expression to `String`, see also [`print_expr`](@ref).
"""
function sprint_expr(ex; context=nothing)
    buf = IOBuffer()
    if context === nothing
        print_expr(buf, ex)
    else
        print_expr(IOContext(buf, context), ex)
    end
    return String(take!(buf))
end

"""
    print_expr([io::IO], ex)

Print a given expression. `ex` can be a `Expr` or a syntax type `JLExpr`.
"""
print_expr(ex) = print_expr(stdout, ex)
print_expr(io::IO, ex) = print_expr(io, ex, PrintState())
print_expr(io::IO, ex, p::PrintState) = print_expr(io, ex, p, Color())

@deprecate print_ast(ex) print_expr(ex)
@deprecate print_ast(io, ex) print_expr(io, ex)

const uni_ops = Set{Symbol}([:(+), :(-), :(!), :(¬), :(~), :(<:), :(>:), :(√), :(∛), :(∜)])
const expr_infix_wide = Set{Symbol}([
    :(=), :(+=), :(-=), :(*=), :(/=), :(\=), :(^=), :(&=), :(|=), :(÷=), :(%=), :(>>>=), :(>>=), :(<<=),
    :(.=), :(.+=), :(.-=), :(.*=), :(./=), :(.\=), :(.^=), :(.&=), :(.|=), :(.÷=), :(.%=), :(.>>>=), :(.>>=), :(.<<=),
    :(&&), :(||), :(<:), :($=), :(⊻=), :(>:), :(-->)])

Base.show(io::IO, def::JLExpr) = print_expr(io, def)

function print_expr(io::IO, ex::JLExpr, ps::PrintState, theme::Color)
    print_expr(io, codegen_ast(ex), ps, theme)
end

function print_expr(io::IO, ex, ps::PrintState, theme::Color)
    @switch ex begin
        @case Expr(:block, line::LineNumberNode, stmt)
            # stmt #= line =#
            print_expr(io, stmt, ps, theme)
            print(io, tab(2))
            print_expr(io, line, ps, theme)
        @case Expr(:block, line1::LineNumberNode, line2::LineNumberNode, stmts...)
            print_kw(io, "begin", ps, theme)
            println(io, ps)
            print_stmts_list(io, ex.args, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:block, stmt1, line::LineNumberNode, stmt2)
            printstyled(io, "("; color=ps.color)
            print_expr(io, stmt1, ps, theme)
            printstyled(io, "; "; color=ps.color)
            print_expr(io, stmt2)
            printstyled(io, ")"; color=ps.color)
        @case Expr(:block, stmts...)
            print_kw(io, "begin", ps, theme)
            println(io, ps)
            print_stmts_list(io, stmts, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:let, vars, body)
            print_kw(io, "let", ps, theme)
            if !isempty(vars.args)
                print(io, tab(1))
                print_collection(io, vars.args, ps, theme)
            end
            println(io, ps)
            print_stmts(io, body, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:if, cond, body, otherwise...) || Expr(:elseif, cond, body, otherwise...)
            within_line(io, ps) do
                print_kw(io, string(ex.head, tab(1)), ps, theme)
                print_expr(io, cond, ps, theme)
            end
            println(io, ps)
            print_stmts(io, body, ps, theme)
            if isempty(otherwise)
                print_end(io, ps, theme)
            else
                otherwise = otherwise[1]
                if Meta.isexpr(otherwise, :elseif)
                    print_expr(io, otherwise, ps, theme)
                else
                    print_kw(io, "else", ps, theme)
                    println(io, ps)
                    print_stmts(io, otherwise, ps, theme)
                    print_end(io, ps, theme)
                end
            end
        @case Expr(:for, head, body)
            within_line(io, ps) do
                print_kw(io, "for ", ps, theme)
                if Meta.isexpr(head, :block)
                    print_collection(io, head.args, ps, theme)
                else
                    print_expr(io, head, ps, theme)
                end
            end
            println(io, ps)
            print_stmts(io, body, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:function, call, body)
            within_line(io, ps) do
                print_kw(io, "function ", ps, theme)
                print_expr(io, call, ps, theme)
            end
            println(io, ps)
            print_stmts(io, body, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:macro, call, body)
            within_line(io, ps) do
                print_kw(io, "macro ", ps, theme)
                print_expr(io, call, ps, theme)
            end
            println(io, ps)
            print_stmts(io, body, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:macrocall, &(GlobalRef(Core, Symbol("@doc"))), line, doc, code)
            print_expr(io, line, ps, theme)
            println(io, ps)
            print(io, tab(ps.line_indent))
            printstyled(io, "\"\"\""; color=theme.string)
            println(io, ps)
            lines = split(doc, '\n')
            for (i, line) in enumerate(lines)
                print(io, tab(ps.line_indent))
                printstyled(io, line; color=theme.string)
                if i != length(lines)
                    println(io, ps)
                end
            end
            print(io, tab(ps.line_indent))
            printstyled(io, "\"\"\""; color=theme.string)
            println(io, ps)
            print_expr(io, code, ps, theme)
        @case Expr(:macrocall, Symbol("@__MODULE__"), line)
            with_color(theme.fn, ps) do
                print_expr(io, Symbol("@__MODULE__"), ps, theme)
            end
        @case Expr(:macrocall, name::Symbol, line, s::String)
            if endswith(string(name), "_str")
                printstyled(io, string(name)[2:end-4]; color=theme.fn)
                print_expr(s)
            else
                print_macro(io, name, line, (s, ), ps, theme)
            end
        @case Expr(:macrocall, name, line, xs...)
            print_macro(io, name, line, xs, ps, theme)
        @case Expr(:struct, ismutable, head, body)
            within_line(io, ps) do
                if ismutable
                    printstyled(io, "mutable "; color=theme.kw)
                end
                printstyled(io, "struct "; color=theme.kw)
                print_expr(io, head, ps, theme)
            end
            println(io, ps)
            print_stmts(io, body, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:try, try_body, catch_var, catch_body)
            print_try(io, try_body, ps, theme)
            print_catch(io, catch_var, catch_body, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:try, try_body, false, false, finally_body)
            print_try(io, try_body, ps, theme)
            print_finally(io, finally_body, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:try, try_body, catch_var, catch_body, finally_body)
            print_try(io, try_body, ps, theme)
            print_catch(io, catch_var, catch_body, ps, theme)
            print_finally(io, finally_body, ps, theme)
            print_end(io, ps, theme)
        @case Expr(:module, notbare, name, body)
            if notbare
                print_kw(io, "module", ps, theme)
            else
                print_kw(io, "baremodule", ps, theme)
            end
            println(io, ps)
            print_stmts(io, body, ps, theme)
            print_end(io, ps, theme)
        @case _
            print_within_line(io, ex, ps, theme)
    end
end

function print_kw(io::IO, x, ps, theme::Color)
    print(io, tab(ps.line_indent))
    printstyled(io, x; color=theme.kw)
end

function print_macro(io, name, line, xs, ps, theme)
    print_expr(io, line, ps, theme)
    println(io, ps)
    within_line(io, ps) do
        with_color(theme.fn, ps) do
            print_expr(io, name, ps, theme)
        end
        print(io, tab(1))
        print_collection(io, xs, ps, theme; delim=tab(1))
    end
end

print_end(io::IO, ps, theme) = print_kw(io, "end", ps, theme)

function Base.println(io::IO, ps::PrintState)
    ps.line_indent = ps.content_indent
    println(io)
end

function print_stmts(io, body, ps::PrintState, theme::Color)
    if body isa Expr && body.head === :block
        print_stmts_list(io, body.args, ps, theme)
    else
        within_indent(ps) do
            print_expr(io, body, ps, theme)
            println(io, ps)
        end
    end
end

function print_try(io, try_body, ps, theme)
    print_kw(io, "try", ps, theme)
    println(io, ps)
    print_stmts(io, try_body, ps, theme)
end

function print_catch(io, catch_var, catch_body, ps, theme)
    print_kw(io, "catch ", ps, theme)
    print_expr(io, catch_var, ps, theme)
    println(io, ps)
    print_stmts(io, catch_body, ps, theme)
end

function print_finally(io, finally_body, ps, theme)
    print_kw(io, "finally", ps, theme)
    println(io, ps)
    print_stmts(io, finally_body, ps, theme)
end

function print_stmts_list(io, stmts, ps::PrintState, theme::Color)
    within_indent(ps) do
        for stmt in stmts
            print_expr(io, stmt, ps, theme)
            println(io, ps)
        end
    end
    return
end

function print_collection(io, xs, ps::PrintState, theme=Color(); delim=", ")
    for i in 1:length(xs)
        print_expr(io, xs[i], ps, theme)
        if i !== length(xs)
            print(io, delim)
        end
    end
end

function with_color(f, name::Symbol, ps::PrintState)
    color = ps.color
    if color === :normal
        ps.color = name
    end
    ret = f()
    ps.color = color
    return ret
end

function within_line(f, io, ps)
    indent = ps.line_indent
    print(io, tab(indent))
    ps.line_indent = 0
    ret = f()
    ps.line_indent = indent
    return ret
end

function within_indent(f, ps)
    ps.line_indent += 4
    ps.content_indent = ps.line_indent
    ret = f()
    ps.line_indent -= 4
    ps.content_indent = ps.line_indent
    return ret
end

function print_within_line(io::IO, ex, ps::PrintState, theme::Color=Color())
    within_line(io, ps) do
        @switch ex begin
            @case ::Number
                printstyled(io, ex; color=theme.literal)
            @case ::String
                printstyled(io, "\"", ex, "\""; color=theme.string)
            @case ::Symbol
                printstyled(io, ex; color=ps.color)
            @case ::Nothing
                printstyled(io, "nothing"; color=:blue)
            @case ::QuoteNode
                if Base.isidentifier(ex.value)
                    print(io, ":", ex.value)
                else
                    print(io, ":(", ex.value, ")")
                end
            @case ::GlobalRef
                printstyled(io, ex.mod, "."; color=ps.color)
                print_expr(io, ex.name, ps, theme)
            @case ::LineNumberNode
                printstyled(io, ex; color=theme.comment)
            @case Expr(:export, xs...)
                print_kw(io, "export ", ps, theme)
                print_collection(io, xs, ps, theme)
            @case Expr(:tuple, Expr(:parameters, kwargs...), args...)
                printstyled(io, "("; color=ps.color)
                print_collection(io, args, ps, theme)
                if !isempty(kwargs)
                    printstyled(io, "; "; color=ps.color)
                    print_collection(io, kwargs, ps)
                end
                printstyled(io, ")"; color=ps.color)
            @case Expr(:tuple, xs...)
                printstyled(io, "("; color=ps.color)
                print_collection(io, xs, ps)
                printstyled(io, ")"; color=ps.color)
            @case Expr(:(::), type)
                printstyled(io, "::"; color=ps.color)
                with_color(theme.type, ps) do
                    print_expr(io, type, ps, theme)
                end
            @case Expr(:(::), name, type)
                print_expr(io, name, ps, theme)
                printstyled(io, "::"; color=ps.color)
                with_color(theme.type, ps) do
                    print_expr(io, type, ps, theme)
                end
            @case Expr(:., a, b::QuoteNode)
                print_expr(io, a, ps, theme)
                printstyled(io, "."; color=ps.color)
                print_expr(io, b.value, ps, theme)
            @case Expr(:., a, b)
                print_expr(io, a, ps, theme)
                printstyled(io, "."; color=ps.color)
                print_expr(io, b, ps, theme)
            @case Expr(:(<:), name, type)
                print_expr(io, name, ps, theme)
                printstyled(io, " <: "; color=ps.color)
                with_color(theme.type, ps) do
                    print_expr(io, type, ps, theme)
                end
            @case Expr(:kw, name, value) || Expr(:(=), name, value)
                print_expr(io, name, ps, theme)
                printstyled(io, tab(1), "=", tab(1); color=ps.color)
                print_expr(io, value, ps, theme)
            @case Expr(:..., name)
                print_expr(io, name, ps, theme)
                print(io, "...")
            @case Expr(:&, name)
                printstyled(io, "&"; color=theme.kw)
                print_expr(io, name, ps, theme)
            @case Expr(:$, name)
                printstyled(io, "\$"; color=theme.kw)
                print(io, "(")
                print_expr(io, name, ps, theme)
                print(io, ")")
            @case Expr(:curly, name, vars...)
                print_expr(io, name, ps, theme)
                print(io, "{")
                with_color(theme.type, ps) do
                    print_collection(io, vars, ps, theme)
                end
                print(io, "}")
            @case Expr(:ref, name, xs...)
                print_expr(io, name, ps, theme)
                print(io, "[")
                print_collection(io, xs, ps, theme)
                print(io, "]")
            @case Expr(:where, body, whereparams...)
                print_expr(io, body, ps, theme)
                printstyled(io, tab(1), "where", tab(1); color=theme.kw)
                print(io, "{")
                print_collection(io, whereparams, ps, theme)
                print(io, "}")
            @case Expr(:call, name, Expr(:parameters, kwargs...), args...)
                print_expr(io, name, ps, theme)
                printstyled(io, "("; color=ps.color)
                print_collection(io, args, ps)
                if !isempty(kwargs)
                    printstyled(io, "; "; color=ps.color)
                    print_collection(io, kwargs, ps)
                end
                printstyled(io, ")"; color=ps.color)
            @case Expr(:call, :(:), xs...)
                print_collection(io, xs, ps, theme; delim=":")
            @case Expr(:call, name::Symbol, x)
                if name in uni_ops
                    print_expr(io, name, ps, theme)
                    print_expr(io, x, ps, theme)
                else
                    print_call_expr(io, name, [x], ps, theme)
                end
            @case Expr(:call, :+, xs...)
                print_collection(io, xs, ps, theme; delim=" + ")
            @case Expr(:call, name, lhs, rhs)
                func_prec = Base.operator_precedence(name_only(name))
                if func_prec > 0
                    print_expr(io, lhs, ps, theme)
                    print(io, tab(1))
                    print_expr(io, name, ps, theme)
                    print(io, tab(1))
                    print_expr(io, rhs, ps, theme)
                else
                    print_call_expr(io, name, [lhs, rhs], ps, theme)
                end
            @case Expr(:call, name, args...)
                print_call_expr(io, name, args, ps, theme)
            @case Expr(:(->), call, body)
                print_expr(io, call, ps, theme)
                printstyled(io, tab(1), "->", tab(1); color=theme.kw)
                print_expr(io, body, ps, theme)
            @case Expr(:return, x)
                printstyled(io, "return", tab(1); color=theme.kw)
                @match x begin
                    Expr(:tuple, xs...) => print_collection(io, xs, ps)
                    _ => print_expr(io, x, ps, theme)
                end
            @case Expr(:string, xs...)
                printstyled(io, "\""; color=theme.string)
                for x in xs
                    if x isa String
                        printstyled(io, x, ; color=theme.string)
                    else
                        printstyled(io, "\$"; color=theme.literal)
                        with_color(theme.literal, ps) do
                            print_expr(io, x, ps, theme)
                        end
                    end
                end
                printstyled(io, "\""; color=theme.string)
            @case Expr(head, lhs, rhs)
                if head in expr_infix_wide
                    print_expr(io, lhs, ps, theme)
                    printstyled(io, tab(1), head, tab(1); color=theme.kw)
                    print_expr(io, rhs, ps, theme)
                else
                    Base.show_unquoted_quote_expr(IOContext(io, :unquote_fallback => true), ex, ps.line_indent, -1, 0)
                end
            @case ::Base.ExprNode
                Base.show_unquoted_quote_expr(IOContext(io, :unquote_fallback => true), ex, ps.line_indent, -1, 0)
            @case _
                print(io, ex)
        end
    end
    return
end

function print_call_expr(io::IO, name, args, ps, theme)
    print_expr(io, name, ps, theme)
    printstyled(io, "("; color=ps.color)
    print_collection(io, args, ps, theme)
    printstyled(io, ")"; color=ps.color)
end


text:
# ---
# title: 967. Numbers With Same Consecutive Differences
# id: problem967
# author: Indigo
# date: 2021-02-19
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/numbers-with-same-consecutive-differences/description/>
# hidden: true
# ---
#
# Return all **non-negative** integers of length `n` such that the absolute
# difference between every two consecutive digits is `k`.
#
# Note that **every** number in the answer **must not** have leading zeros
# **except** for the number `0` itself. For example, `01` has one leading zero
# and is invalid, but `0` is valid.
#
# You may return the answer in **any order**.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3, k = 7
#     Output: [181,292,707,818,929]
#     Explanation: Note that 070 is not a valid number, because it has leading zeroes.
#
#
# **Example 2:**
#
#
#
#     Input: n = 2, k = 1
#     Output: [10,12,21,23,32,34,43,45,54,56,65,67,76,78,87,89,98]
#
#
# **Example 3:**
#
#
#
#     Input: n = 2, k = 0
#     Output: [11,22,33,44,55,66,77,88,99]
#
#
# **Example 4:**
#
#
#
#     Input: n = 2, k = 1
#     Output: [10,12,21,23,32,34,43,45,54,56,65,67,76,78,87,89,98]
#
#
# **Example 5:**
#
#
#
#     Input: n = 2, k = 2
#     Output: [13,20,24,31,35,42,46,53,57,64,68,75,79,86,97]
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 9`
#   * `0 <= k <= 9`
#
#
## @lc code=start
using LeetCode

function nums_same_consec_diff(n::Int, k::Int)
    (k == 0) && return [(10 ^ n - 1) ÷ 9 * i for i in 1:9]
    function nums_same_consec_diff(n::Int, k::Int, pre::Int)
        (pre < 0 || pre > 9) && return String[]
        (n == 0) && return [""]
        n1, n2 = pre - k, pre + k
        r1 = nums_same_consec_diff(n - 1, k, n1)
        r2 = nums_same_consec_diff(n - 1, k, n2)
        [[string(n1) * r for r in r1]; [string(n2) * r for r in r2]]
    end
    res = Int[]
    for i in 1:9
        tmp = nums_same_consec_diff(n - 1, k, i)
        append!(res, parse(Int, string(i) * t) for t in tmp)
    end
    res
end
## @lc code=end


text:
"""
```
pseudo_measurement{T<:AbstractFloat}(m::Model904{T},
    TTT::Matrix{T}, RRR::Matrix{T}, CCC::Vector{T})
```

Assign pseudo-measurement equation (a linear combination of states):

```
x_t = ZZ_pseudo*s_t + DD_pseudo
```
"""
function pseudo_measurement(m::Model904{T},
                                              TTT::Matrix{T},
                                              RRR::Matrix{T},
                                              CCC::Vector{T}) where {T<:AbstractFloat}
    endo   = m.endogenous_states
    pseudo = m.pseudo_observables

    _n_states = n_states_augmented(m)
    _n_pseudo = n_pseudo_observables(m)

    # Initialize pseudo ZZ and DD matrices
    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)

    ##########################################################
    ## PSEUDO-OBSERVABLE EQUATIONS
    ##########################################################

    ## Output
    ZZ_pseudo[pseudo[:y_t],endo[:y_t]] = 1.

    ## Flexible Output
    ZZ_pseudo[pseudo[:y_f_t],endo[:y_f_t]] = 1.

    ## Natural Rate
    ZZ_pseudo[pseudo[:NaturalRate],endo[:r_f_t]] = 1.
    DD_pseudo[pseudo[:NaturalRate]]              = 100 .* (m[:rstar]-1.)

    ## π_t
    ZZ_pseudo[pseudo[:π_t],endo[:π_t]] = 1.
    DD_pseudo[pseudo[:π_t]]            = 100*(m[:π_star]-1);

    ## Output Gap
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_t]] = 1;
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_f_t]] = -1;

    ## Ex Ante Real Rate
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:R_t]]  = 1;
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:Eπ_t]] = -1;
    DD_pseudo[pseudo[:ExAnteRealRate]]             = m[:Rstarn] - 100*(m[:π_star]-1);

    ## Long Run Inflation
    ZZ_pseudo[pseudo[:LongRunInflation],endo[:π_star_t]] = 1.
    DD_pseudo[pseudo[:LongRunInflation]]                 = 100. *(m[:π_star]-1.)

    ## Marginal Cost
    ZZ_pseudo[pseudo[:MarginalCost],endo[:mc_t]] = 1.

    ## Wages
    ZZ_pseudo[pseudo[:Wages],endo[:w_t]] = 1.

    ## Flexible Wages
    ZZ_pseudo[pseudo[:FlexibleWages],endo[:w_f_t]] = 1.

    ## Hours
    ZZ_pseudo[pseudo[:Hours],endo[:L_t]] = 1.

    ## Flexible Hours
    ZZ_pseudo[pseudo[:FlexibleHours],endo[:L_f_t]] = 1.

    ## z_t
    ZZ_pseudo[pseudo[:z_t], endo[:z_t]] = 1.

    ## Z_t
    ZZ_pseudo[pseudo[:Z_t], endo[:z_t]] = 1.
    DD_pseudo[pseudo[:Z_t]]             = 100. * exp(m[:z_star]-1)

    ## ztil_t
    ZZ_pseudo[pseudo[:ztil_t], endo[:ztil_t]] = 1.

    ## Nominal FFR
    ZZ_pseudo[pseudo[:NominalFFR], endo[:R_t]] = 1.
    DD_pseudo[pseudo[:NominalFFR]] = m[:Rstarn]

    ## Real FFR
    ZZ_pseudo[pseudo[:RealFFR], endo[:R_t]] = 1.
    ZZ_pseudo[pseudo[:RealFFR], endo[:π_t]] = -1.
    DD_pseudo[pseudo[:RealFFR]] =  m[:Rstarn] - (100. * (m[:π_star] - 1.))

    ## Expected nominal natural rate
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:r_f_t]] = 1.
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:Eπ_t]]  = 1.
    DD_pseudo[pseudo[:ExpectedNominalNaturalRate]]               =  m[:Rstarn]

    ## S^Inf
    betabar     = (exp((1-m[:σ_c]) * m[:z_star]) * m[:β])
    S1 = inv(Matrix(1.0I, size(TTT, 1), size(TTT, 1)) - betabar*TTT)
    ZZ_pseudo[pseudo[:SInf], :] = S1[endo[:mc_t],:]
    ZZ_pseudo[pseudo[:SInf_t1], pseudo[:SInf]] = 1.0
    ZZ_pseudo[pseudo[:SInf_t2], pseudo[:SInf]] = 1.0
    ZZ_pseudo[pseudo[:SInf_t3], pseudo[:SInf]] = 1.0

    # π_til
    kappa =  ((1 - m[:ζ_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))*(1 - m[:ζ_p]))/(m[:ζ_p]*((m[:Φ]- 1)*m[:ϵ_p] + 1))/(1 + m[:ι_p]*m[:β]*exp((1 - m[:σ_c])*m[:z_star]))
    exp( (m[:σ_c] -1) * m[:z_star]) / m[:β]
    kappa_coef = kappa * (1+m[:ι_p]*betabar)
    ZZ_pseudo[pseudo[:Pi_Til], pseudo[:SInf]] = kappa_coef
    ZZ_pseudo[pseudo[:Pi_Til], pseudo[:SInf_t1]] = kappa_coef*m[:ι_p]
    ZZ_pseudo[pseudo[:Pi_Til], pseudo[:SInf_t2]] = kappa_coef*m[:ι_p]^2
    ZZ_pseudo[pseudo[:Pi_Til], pseudo[:SInf_t3]] = kappa_coef*m[:ι_p]^3
    DD_pseudo[pseudo[:Pi_Til]] = 100*(m[:π_star] -1)

    # Collect indices and transforms
    return PseudoMeasurement(ZZ_pseudo, DD_pseudo)
end


text:
#
# Op
#

struct Op{N, Params <: NamedTuple}
  name::String
  sites::NTuple{N,Int}
  params::Params
end

Op(name::String, sites::Tuple{Vararg{Int}}) =
  Op(name, sites, NamedTuple())

Op(name::String, site::Int,
   params::NamedTuple = NamedTuple()) =
  Op(name, (site,), params)

Op(name::String, sites::Int...) =
  Op(name, sites)

Base.convert(::Type{Op}, t::Tuple) = Op(t...)

function ITensors.op(sites::Vector{<:Index},
                     o::Op)
  return op(sites, o.name, o.sites...; o.params...)
end

function Base.show(io::IO, o::Op)
  print(io, "\"$(o.name)\"")
  if length(o.sites) == 1
    print(io, "($(only(o.sites)))")
  else
    print(io, o.sites)
  end
  if !isempty(o.params)
    print(io, " ")
    if length(o.params) == 1
      print(io, "($(only(keys(o.params))) = $(only(o.params)))")
    else
      print(io, o.params)
    end
  end
end

#
# ProductOps
#

struct ProductOps
  data::Vector{Op}
end

ProductOps() = ProductOps(Op[])

Base.copy(C::ProductOps) = ProductOps(copy(C.data))

Base.push!(C::ProductOps, O) =
  (push!(C.data, O);
   return C)

Base.pushfirst!(C::ProductOps, O) =
  (pushfirst!(C.data, O);
   return C)

Base.iterate(C::ProductOps, args...) = iterate(C.data, args...)
Base.length(C::ProductOps) = length(C.data)

import Base: *

C::ProductOps * O = push!(copy(C), O)

O * C::ProductOps = pushfirst!(copy(C), O)

C1::ProductOps * C2::ProductOps = ProductOps(vcat(C1.data, C2.data))

function Base.show(io::IO, C::ProductOps)
  println("ProductOps")
  for o in C.data
    println(io, o)
  end
end

"""
    ops(::Vector{Index}, C::ProductOps)

Return a Vector of ITensors corresponding to the input circuit.
"""
function ops(s::Vector{<:Index}, C::ProductOps)
  return [op(s, c) for c in C]
end


text:
# cellid, streak
const CountIdentifier = Tuple{UInt,UInt8}

struct PolyCounter{NUC,R,C}
    streakcounter::Accumulator{CountIdentifier,Int}
    lastcounter::Accumulator{CountIdentifier,Int}
    nuccounter::Accumulator{CountIdentifier,Int}
end

function PolyCounter(nuc::Char,read::Int,bclen::Int)
    PolyCounter{nuc,read,bclen}(
        counter(CountIdentifier),
        counter(CountIdentifier),
        counter(CountIdentifier))
end

function Base.push!(pc::PolyCounter{NUC,R},ir::InterpretedRecord) where {NUC,R}
    seq = ir.records[R].data[ir.records[R].sequence]
    streak, last, count = longeststreak(seq,UInt8(NUC))
    push!(pc.streakcounter,(cellid(ir),UInt8(streak)))
    push!(pc.lastcounter,(cellid(ir),UInt8(last)))
    push!(pc.nuccounter,(cellid(ir),UInt8(count)))
end

function Base.merge!(pc::PolyCounter{NUC,R},
                     counters::PolyCounter{NUC,R}...) where {NUC,R}
    merge!(pc.streakcounter,(c.streakcounter for c in counters)...)
    merge!(pc.lastcounter,(c.lastcounter for c in counters)...)
    merge!(pc.nuccounter,(c.nuccounter for c in counters)...)
    return pc
end

function Base.merge(counters::PolyCounter{NUC,R,C}...) where {NUC,R,C}
    pc = PolyCounter(NUC,R,C)
    return merge!(pc,counters...)
end

function Base.write(name::String,pc::PolyCounter{NUC,R,C}) where {NUC,R,C}
    if isdir(name)
        filename = joinpath(name,"poly$NUC-read$R")
    else
        filename = name
    end
    bc(id) = idtoname(id,C)
    streaks = DataFrame(
        cellid = [bc(id) for ((id,streak),count) in pc.streakcounter],
        streak = [streak for ((id,streak),count) in pc.streakcounter],
        count  = [count  for ((id,streak),count) in pc.streakcounter])
    writetable(filename*"-streaks.csv",unstack(streaks,:cellid,:streak,:count))
    last = DataFrame(
        cellid = [bc(id) for ((id,last),count)  in pc.lastcounter],
        last   = [last   for ((id,last),count)  in pc.lastcounter],
        count  = [count  for ((id,last),count)  in pc.lastcounter])
    writetable(filename*"-last.csv",unstack(last,:cellid,:last,:count))
    nuc = DataFrame(
        cellid = [bc(id) for ((id,nuc),count)   in pc.nuccounter],
        nuc    = [nuc    for ((id,nuc),count)   in pc.nuccounter],
        count  = [count  for ((id,nuc),count)   in pc.nuccounter])
    writetable(filename*"-nuc.csv",unstack(nuc,:cellid,:nuc,:count),nastring="0")
end


text:
export RandomVariable, randomvar

import Base: show

## Random variable implementation

mutable struct RandomVariable <: AbstractVariable
    name                :: Symbol
    collection_type     :: AbstractVariableCollectionType
    input_messages      :: Vector{MessageObservable{AbstractMessage}}
    output_messages     :: Vector{MessageObservable{Message}}
    output_initialised  :: Bool
    output_cache        :: Union{Nothing, EqualityChain}
    marginal            :: MarginalObservable
    pipeline            :: AbstractPipelineStage
    prod_constraint
    prod_strategy
    form_constraint
    form_check_strategy
end

Base.show(io::IO, randomvar::RandomVariable) = print(io, "RandomVariable(", indexed_name(randomvar), ")")

function randomvar(name::Symbol, collection_type::AbstractVariableCollectionType = VariableIndividual(); pipeline = EmptyPipelineStage(), prod_constraint = ProdAnalytical(), prod_strategy = FoldLeftProdStrategy(), form_constraint = UnspecifiedFormConstraint(), form_check_strategy = FormConstraintCheckPickDefault())
    return RandomVariable(name, collection_type, Vector{MessageObservable{AbstractMessage}}(), Vector{MessageObservable{Message}}(), false, nothing, MarginalObservable(), pipeline, prod_constraint, prod_strategy, form_constraint, form_check_strategy)
end

function randomvar(name::Symbol, dims::Tuple; pipeline = EmptyPipelineStage(), prod_constraint = ProdAnalytical(), prod_strategy = FoldLeftProdStrategy(), form_constraint = UnspecifiedFormConstraint(), form_check_strategy = FormConstraintCheckPickDefault())
    return randomvar(name, dims...; pipeline = pipeline, prod_constraint = prod_constraint, prod_strategy = prod_strategy, form_constraint = form_constraint, form_check_strategy = form_check_strategy)
end

function randomvar(name::Symbol, length::Int; pipeline = EmptyPipelineStage(), prod_constraint = ProdAnalytical(), prod_strategy = FoldLeftProdStrategy(), form_constraint = UnspecifiedFormConstraint(), form_check_strategy = FormConstraintCheckPickDefault())
    vars = Vector{RandomVariable}(undef, length)
    @inbounds for i in 1:length
        vars[i] = randomvar(name, VariableVector(i); pipeline = pipeline, prod_constraint = prod_constraint, prod_strategy = prod_strategy, form_constraint = form_constraint, form_check_strategy = form_check_strategy)
    end
    return vars
end

function randomvar(name::Symbol, dims::Vararg{Int}; pipeline = EmptyPipelineStage(), prod_constraint = ProdAnalytical(), prod_strategy = FoldLeftProdStrategy(), form_constraint = UnspecifiedFormConstraint(), form_check_strategy = FormConstraintCheckPickDefault())
    vars = Array{RandomVariable}(undef, dims)
    @inbounds for i in CartesianIndices(axes(vars))
        vars[i] = randomvar(name, VariableArray(i); pipeline = pipeline, prod_constraint = prod_constraint, prod_strategy = prod_strategy, form_constraint = form_constraint, form_check_strategy = form_check_strategy)
    end
    return vars
end

degree(randomvar::RandomVariable)              = length(randomvar.input_messages)
name(randomvar::RandomVariable)                = randomvar.name
collection_type(randomvar::RandomVariable)     = randomvar.collection_type
equality_chain(randomvar::RandomVariable)      = randomvar.equality_chain
prod_constraint(randomvar::RandomVariable)     = randomvar.prod_constraint
prod_strategy(randomvar::RandomVariable)       = randomvar.prod_strategy
form_constraint(randomvar::RandomVariable)     = randomvar.form_constraint
form_check_strategy(randomvar::RandomVariable) = _form_check_strategy(randomvar.form_check_strategy, randomvar)

_form_check_strategy(::FormConstraintCheckPickDefault, randomvar::RandomVariable) = default_form_check_strategy(form_constraint(randomvar))
_form_check_strategy(form_check_strategy, randomvar::RandomVariable)              = form_check_strategy

messages_prod_fn(randomvar::RandomVariable) = messages_prod_fn(prod_strategy(randomvar), prod_constraint(randomvar), UnspecifiedFormConstraint(), default_form_check_strategy(UnspecifiedFormConstraint()))
marginal_prod_fn(randomvar::RandomVariable) = marginal_prod_fn(prod_strategy(randomvar), prod_constraint(randomvar), form_constraint(randomvar), form_check_strategy(randomvar))

getlastindex(randomvar::RandomVariable) = degree(randomvar) + 1

messagein(randomvar::RandomVariable, index::Int)  = @inbounds randomvar.input_messages[index]

function messageout(randomvar::RandomVariable, index::Int)
    if randomvar.output_initialised === false
        initialize_output_messages!(randomvar)
    end
    return @inbounds randomvar.output_messages[index]
end

get_pipeline_stages(randomvar::RandomVariable)        = randomvar.pipeline
add_pipeline_stage!(randomvar::RandomVariable, stage) = randomvar.pipeline = (randomvar.pipeline + stage)

_getmarginal(randomvar::RandomVariable)              = randomvar.marginal
_setmarginal!(randomvar::RandomVariable, observable) = connect!(_getmarginal(randomvar), observable)
_makemarginal(randomvar::RandomVariable)             = collectLatest(AbstractMessage, Marginal, randomvar.input_messages, marginal_prod_fn(randomvar))

function setmessagein!(randomvar::RandomVariable, index::Int, messagein)
    if index === degree(randomvar) + 1
        push!(randomvar.input_messages, messagein)
    else
        error("Inconsistent state in setmessagein! function for random variable $(randomvar). `index` should be equal to `degree(randomvar) + 1 = $(degree(randomvar) + 1)`, $(index) is given instead")
    end
end

function activate!(model, randomvar::RandomVariable)
    # `5` here is empirical observation, maybe we can come up with better heuristic?
    # in case if number of connections is large we use cache equality nodes chain structure
    if degree(randomvar) > 5
        chain_pipeline  = schedule_on(global_reactive_scheduler(getoptions(model)))
        chain_prod_fn   = messages_prod_fn(randomvar)
        randomvar.output_cache = EqualityChain(randomvar.input_messages, chain_pipeline, chain_prod_fn)
    end

    _setmarginal!(randomvar, _makemarginal(randomvar))

    return nothing
end

initialize_output_messages!(randomvar::RandomVariable) = initialize_output_messages!(randomvar.output_cache, randomvar)

# Generic fallback for variables with small number of connected nodes, somewhere <= 5
# We do not create equality chain in this cases, but simply do eager product
function initialize_output_messages!(::Nothing, randomvar::RandomVariable)
    d          = degree(randomvar)
    outputmsgs = randomvar.output_messages
    inputmsgs  = randomvar.input_messages
    prod_fn    = messages_prod_fn(randomvar)

    resize!(outputmsgs, d)

    @inbounds for i in 1:d
        outputmsgs[i] = MessageObservable(Message)
        outputmsg     = collectLatest(AbstractMessage, Message, skipindex(inputmsgs, i), prod_fn)
        connect!(outputmsgs[i], outputmsg)
    end

    randomvar.output_initialised = true

    return nothing
end

# Equality chain initialisation for variables with large number of connected nodes, somewhere > 5
# In this cases it is more efficient to create an equality chain structure, but it does allocate way more memory
function initialize_output_messages!(chain::EqualityChain, randomvar::RandomVariable)
    d          = degree(randomvar)
    outputmsgs = randomvar.output_messages

    resize!(outputmsgs, d)

    @inbounds for i in 1:d
        outputmsgs[i] = MessageObservable(Message)
    end

    initialize!(chain, outputmsgs)
    randomvar.output_initialised = true

    return nothing
end


text:
save_model_summary(model::TierraModel) = save_model_summary(model, "final")

function save_model_summary(model::TierraModel, prefix)
    counter = Dict{String, Integer}()

    for key in model.organism_keys
        organism = model.organisms[key]

        if haskey(counter, organism.hash)
            counter[organism.hash] += 1
        else
            counter[organism.hash] = 1
        end
    end

    open(LOG_FILE * string(prefix), "w+") do io

        for (hash, num) in counter
            write(io, hash * ": \t " * string(num) * "\n")
        end

        for key in model.organism_keys
            organism = model.organisms[key]
            Base.show(io, organism, model)
            write(io, "\n")
        end

    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module ProjectEuler

include("useful_functions.jl")
include("001.jl")
include("002.jl")
include("003.jl")
include("004.jl")
include("005.jl")

include("007.jl")
include("008.jl")
include("009.jl")
include("010.jl")

include("069.jl")

end


text:
import CausalityToolsBase:
    TriangulationBinningScheme,
    TriangulationBinning,
    SimplexIntersectionType,
    ApproximateIntersection

import ..TransferOperators:
    transferoperator_triangulation_approx

import StateSpaceReconstruction:
    invariantize,
    triangulate

export TriangulationApproxInvariantMeasure, triangulationapproxinvariantmeasure

"""
TriangulationApproxInvariantMeasure

An invariant measure estimated from a triangulation using the approximate simplex
intersection estimator `transferoperator_triangulation_approx` for the
approximation of the transfer matrix.

## Fields
- **`points`**: The points from which the measure is estimated.
- **`invariantized_points`**: The points from which the measure is estimated,
possibly with the last point moved toward the center of the convex hull of
the points to ensure no information is lost when computing the transfer matrix.
- **`triangulation`**. A triangulation of `invariantized_points[1:(end - 1)]` which is used to estimate the transfer matrix.

"""
struct TriangulationApproxInvariantMeasure{PT, TT, TOT} <: AbstractTriangulationInvariantMeasure
    points::PT
    invariantized_points::PT
    triangulation::TT
    transferoperator::TOT
    measure::InvariantDistribution
end

function triangulationapproxinvariantmeasure(points::Vector{Vector{T}}; n_sample_pts = 200,
    sample_randomly::Bool = false, kwargs...) where {T}

    inv_points = invariantize(points)
    triang = triangulate(inv_points[1:(end - 1)])

    TO = transferoperator_triangulation_approx(points,
        triang.simplexindices,
        n_sample_pts = n_sample_pts, sample_randomly = sample_randomly)

    measure = invariantmeasure(TO, kwargs...)

    TriangulationApproxInvariantMeasure(
        points,
        inv_points,
        triang,
        TO,
        measure)
end

function description(μ::TriangulationApproxInvariantMeasure)
    T_pts = typeof(μ.points)
    T_invpts = typeof(μ.invariantized_points)
    triang = μ.triangulation
    to = μ.transferoperator
    measure = μ.measure

    npts =

    invmeasure_type = typeof(μ)
    return join(["TriangulationApproxInvariantMeasure", "\n",
                "  μ.points:\t\t", T_pts,   "\n",
                 "  μ.invariantized_points:\t", T_invpts, "\n",
                 "  μ.triangulation:\t", triang,
                "  μ.transferoperator:\t", to,
                "  μ.measure:\t\t", measure])
end

Base.show(io::IO, μ::TriangulationApproxInvariantMeasure) = println(io, description(μ))

export TriangulationApproxInvariantMeasure, triangulationapproxinvariantmeasure


text:
"""
    struct GreedyClique
        vertices::BitArray
        mutuals::BitArray
    end

# Arguments
- `vertices::AbstractVector{Bool}`: Logical vector; true for each vertex that is a member
    of this clique.
- `mutuals::AbstractVector{Bool}`: Vertices that neighbour the clique (potential clique members)
"""
struct GreedyClique <: AbstractClique
    mutuals::BitArray
    vertices::BitArray
    function GreedyClique(m::AbstractVector{Bool}, v::AbstractVector{Bool})
        size(v) == size(m) || error(LOGGER, DimensionMismatch("size mismatch"))
        return new(BitVector(m), BitVector(v))
    end
end

"""
    vertices(c::GreedyClique) -> BitArray

Return vertices of a clique as a logical vector.
"""
vertices(c::GreedyClique) = c.vertices

"""
    mutuals(c::GreedyClique) -> BitArray

Return mutual neighbours (potential clique members) as a logical vector.
"""
mutuals(c::GreedyClique) = c.mutuals

"""
    mergeable(c1::GreedyCliqye, c2::GreedyClique) -> Bool

Compare two cliques. They can be merged if the mutuals (neighbours) of one clique (c1)
AND the vertices of the other clique (c2) match the vertices of c2.
"""
function mergeable(c1::GreedyClique, c2::GreedyClique)
    return all(mutuals(c1) .& vertices(c2) .== vertices(c2))
end

"""
    union(c1::GreedyClique, c2::GreedyClique) -> GreedyClique

Merges two cliques.
"""
function Base.union(c1::GreedyClique, c2::GreedyClique)
    return GreedyClique(mutuals(c1) .& mutuals(c2), vertices(c1) .| vertices(c2))
end

"""
    greedycliquing(A::AbstractMatrix{Bool}, minsize::Int)
    greedycliquing!(A::AbstractMatrix{Bool}, minsize::Int)

Greedy Cliquing Algorithm based off of previous MATLAB version.

Splits an adjacency matrix into cliques[^1] and remaining non-cliquable
nodes. Compared to a previous version that achieves
maximum node set reduction through optimal cliquing, this version
finds a maximal clique, removes its nodes from the adjacency matrix,
and repeats until no cliques can be found. More discussion can be
found in this[^2] Google doc.

[^1]: <http://en.wikipedia.org/wiki/Clique_(graph_theory)>
[^2]: Private Invenia document: <https://docs.google.com/a/invenia.ca/document/d/1lLRfQT_TFJi1bTfIu_jGxDs1u9A9HYKX5s4Zta3VRNg/edit?usp=sharing>

# Arguments
- `A::AbstractMatrix{Bool}`: Adjacency Matrix
- `minsize::Int`: Min greedyclique size

# Returns
- `cliques::Vector{Clique}`: Vector array of type Clique
- `singletons::Vector{Bool}`: Vector array of Bool. These indicate nodes not in a clique
"""
function greedycliquing(A::AbstractMatrix{Bool}, minsize::Integer)
    return greedycliquing!(collect(A), minsize)
end

# Avoid promoting BitMatrix to Matrix{Bool} as it has performance implications.
function greedycliquing(A::BitMatrix, minsize::Integer)
    return greedycliquing!(copy(A), minsize)
end

function greedycliquing!(A::Union{BitMatrix,Matrix{Bool}}, minsize::Integer)
    issymmetric(A) || error(LOGGER, ArgumentError("Input matrix not symmetric: $A"))
    removediagonal!(A)
    num_nodes = size(A, 1)
    singletons = fill(true, num_nodes)
    cliques = GreedyClique[]

    debug(LOGGER, "Searching for cliques within adjacency matrix A with $num_nodes nodes.")
    for i in 1:num_nodes
        push!(cliques, anymaxclique(A))
        clique_size = count(vertices(cliques[i]))
        num_nodes -= clique_size

        if clique_size < minsize
            # we won't always have found a maximal clique of at least size `minsize` yet, so
            # just remove it from our vector of Cliques
            pop!(cliques)
            break

        elseif num_nodes == 0
            # These two lines not in the original MATLAB code but needed
            v = vertices(cliques[i])
            singletons[v] .= false
            break
        end

        v = vertices(cliques[i])
        singletons[v] .= false
        A[v, :] .= false
        A[:, v] .= false
    end
    cliques = Clique[Clique(vertices(c)) for c in cliques]
    return cliques, singletons
end

"""
    anymaxclique(A::AbstractMatrix{Bool})

Finds any maximal clique in the adjacency matrix. To find a
maximal clique, it sorts the nodes by number of neighbours then if the
node can be added to the clique, it adds it to the clique. Note that this is not the
maximum clique.

# Arguments
- `A::AbstractMatrix{Bool}`: Adjacency Matrix.

# Returns
- `clique::GreedyClique`: One of the maximal cliques in the adjacency matrix.
"""
function anymaxclique(A::AbstractMatrix{Bool})
    numNeighbours = vec(sum(A; dims=2))
    m, n = size(A)
    maxSize = maximum(numNeighbours)
    idx = sortperm(numNeighbours; rev=true)
    clique = GreedyClique(A[idx[1], :], ind2log(idx[1], m))
    vertexCount = 1
    for i in idx[2:end]
        currC = GreedyClique(A[i, :], ind2log(i, m))
        if mergeable(clique, currC)
            clique = union(clique, currC)
            vertexCount += 1
            if vertexCount > maxSize
                return clique
            end
        end
    end
    return clique
end

# Sets square logical matrix diagonal to false.
removediagonal!(A::AbstractMatrix{Bool}) = (A[diagind(A)] .= false; nothing)

# Create a Bool vector of length n that is true at the index of idx and false otherwise
ind2log(idx::Integer, n::Int) = setindex!(fill(false, (n,)), true, idx)


text:
"""
    defaultFiles(s, pids)

Make a good set of filenames for saving a dataset.
"""
function defaultFiles(s, pids)
    return [String(s) * "-$i.slice" for i in eachindex(pids)]
end

"""
    dstore(sym::Symbol, pids, files=defaultFiles(sym,pids))

Export the content of symbol `sym` by each worker specified by `pids` to a
corresponding filename in `files`.
"""
function dstore(sym::Symbol, pids, files = defaultFiles(sym, pids))
    dmap(
        files,
        (fn) -> Base.eval(Main, :(
            begin
                open(f -> $serialize(f, $sym), $fn, "w")
                nothing
            end
        )),
        pids,
    )
    nothing
end

"""
    dstore(dInfo::Dinfo, files=defaultFiles(dInfo.val, dInfo.workers))

Overloaded functionality for `Dinfo`.
"""
function dstore(dInfo::Dinfo, files = defaultFiles(dInfo.val, dInfo.workers))
    dstore(dInfo.val, dInfo.workers, files)
end

"""
    dload(sym::Symbol, pids, files=defaultFiles(sym,pids))

Import the content of symbol `sym` by each worker specified by `pids` from the
corresponding filename in `files`.
"""
function dload(sym::Symbol, pids, files = defaultFiles(sym, pids))
    dmap(files, (fn) -> Base.eval(Main, :(
        begin
            $sym = open($deserialize, $fn)
            nothing
        end
    )), pids)
    return Dinfo(sym, pids)
end

"""
    dload(dInfo::Dinfo, files=defaultFiles(dInfo.val, dInfo.workers))

Overloaded functionality for `Dinfo`.
"""
function dload(dInfo::Dinfo, files = defaultFiles(dInfo.val, dInfo.workers))
    dload(dInfo.val, dInfo.workers, files)
end

"""
    dunlink(sym::Symbol, pids, files=defaultFiles(sym,pids))

Remove the files created by `dstore` with the same parameters.
"""
function dunlink(sym::Symbol, pids, files = defaultFiles(sym, pids))
    dmap(files, (fn) -> rm(fn), pids)
    nothing
end

"""
    dunlink(dInfo::Dinfo, files=defaultFiles(dInfo.val, dInfo.workers))

Overloaded functionality for `Dinfo`.
"""
function dunlink(dInfo::Dinfo, files = defaultFiles(dInfo.val, dInfo.workers))
    dunlink(dInfo.val, dInfo.workers, files)
end


text:
###
### Sampler states
###

mutable struct HMCState{
    TV <: TypedVarInfo,
    TTraj<:AHMC.AbstractTrajectory,
    TAdapt<:AHMC.Adaptation.AbstractAdaptor,
    PhType <: AHMC.PhasePoint
} <: AbstractSamplerState
    vi       :: TV
    eval_num :: Int
    i        :: Int
    traj     :: TTraj
    h        :: AHMC.Hamiltonian
    adaptor  :: TAdapt
    z        :: PhType
end

##########################
# Hamiltonian Transition #
##########################

struct HamiltonianTransition{T, NT<:NamedTuple, F<:AbstractFloat} <: AbstractTransition
    θ    :: T
    lp   :: F
    stat :: NT
end

function HamiltonianTransition(spl::Sampler{<:Hamiltonian}, t::T) where T<:AHMC.Transition
    theta = tonamedtuple(spl.state.vi)
    lp = getlogp(spl.state.vi)
    return HamiltonianTransition(theta, lp, t.stat)
end

transition_type(spl::Sampler{<:Union{StaticHamiltonian, AdaptiveHamiltonian}}) =
    HamiltonianTransition

function additional_parameters(::Type{<:HamiltonianTransition})
    return [:lp,:stat]
end

###
### Hamiltonian Monte Carlo samplers.
###

"""
    HMC(ϵ::Float64, n_leapfrog::Int)

Hamiltonian Monte Carlo sampler with static trajectory.

Arguments:

- `ϵ::Float64` : The leapfrog step size to use.
- `n_leapfrog::Int` : The number of leapfrop steps to use.

Usage:

```julia
HMC(0.05, 10)
```

Tips:

- If you are receiving gradient errors when using `HMC`, try reducing the
`step_size` parameter, e.g.

```julia
# Original step_size
sample(gdemo([1.5, 2]), HMC(0.1, 10), 1000)

# Reduced step_size.
sample(gdemo([1.5, 2]), HMC(0.01, 10), 1000)
```
"""
mutable struct HMC{AD, space, metricT <: AHMC.AbstractMetric} <: StaticHamiltonian{AD}
    ϵ           ::  Float64   # leapfrog step size
    n_leapfrog  ::  Int       # leapfrog step number
end

transition_type(::Sampler{<:Hamiltonian}) = Transition
alg_str(::Sampler{<:Hamiltonian}) = "HMC"

HMC(args...) = HMC{ADBackend()}(args...)
function HMC{AD}(ϵ::Float64, n_leapfrog::Int, ::Type{metricT}, space::Tuple) where {AD, metricT <: AHMC.AbstractMetric}
    return HMC{AD, space, metricT}(ϵ, n_leapfrog)
end
function HMC{AD}(
    ϵ::Float64,
    n_leapfrog::Int,
    ::Tuple{};
    kwargs...
) where AD
    return HMC{AD}(ϵ, n_leapfrog; kwargs...)
end
function HMC{AD}(
    ϵ::Float64,
    n_leapfrog::Int,
    space::Symbol...;
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return HMC{AD}(ϵ, n_leapfrog, metricT, space)
end

function sample_init!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Hamiltonian},
    N::Integer;
    verbose::Bool=true,
    resume_from=nothing,
    kwargs...
) where T<:Hamiltonian

    # Resume the sampler.
    set_resume!(spl; resume_from=resume_from, kwargs...)

    # Get `init_theta`
    initialize_parameters!(spl; verbose=verbose, kwargs...)

    # Set the defualt number of adaptations, if relevant.
    if spl.alg isa AdaptiveHamiltonian
        # If there's no chain passed in, verify the n_adapts.
        if resume_from === nothing
            if spl.alg.n_adapts == 0
                n_adapts_default = Int(round(N / 2))
                spl.alg.n_adapts = n_adapts_default > 1000 ? 1000 : n_adapts_default
            else
                # Verify that n_adapts is less than the samples to draw.
                spl.alg.n_adapts < N || !ismissing(resume_from) ?
                    nothing :
                    throw(ArgumentError("n_adapt of $(spl.alg.n_adapts) is greater than total samples of $N."))
            end
        else
            spl.alg.n_adapts = 0
        end
    end

    # Convert to transformed space if we're using
    # non-Gibbs sampling.
    if !islinked(spl.state.vi, spl) && spl.selector.tag == :default
        link!(spl.state.vi, spl)
        runmodel!(model, spl.state.vi, spl)
    end
end

"""
    HMCDA(n_adapts::Int, δ::Float64, λ::Float64; ϵ::Float64=0.0)

Hamiltonian Monte Carlo sampler with Dual Averaging algorithm.

Usage:

```julia
HMCDA(200, 0.65, 0.3)
```

Arguments:

- `n_adapts::Int` : Numbers of samples to use for adaptation.
- `δ::Float64` : Target acceptance rate. 65% is often recommended.
- `λ::Float64` : Target leapfrop length.
- `ϵ::Float64=0.0` : Inital step size; 0 means automatically search by Turing.

For more information, please view the following paper ([arXiv link](https://arxiv.org/abs/1111.4246)):

- Hoffman, Matthew D., and Andrew Gelman. "The No-U-turn sampler: adaptively
  setting path lengths in Hamiltonian Monte Carlo." Journal of Machine Learning
  Research 15, no. 1 (2014): 1593-1623.
"""
mutable struct HMCDA{AD, space, metricT <: AHMC.AbstractMetric} <: AdaptiveHamiltonian{AD}
    n_adapts    ::  Int         # number of samples with adaption for ϵ
    δ           ::  Float64     # target accept rate
    λ           ::  Float64     # target leapfrog length
    ϵ           ::  Float64     # (initial) step size
end
HMCDA(args...; kwargs...) = HMCDA{ADBackend()}(args...; kwargs...)
function HMCDA{AD}(n_adapts::Int, δ::Float64, λ::Float64, ϵ::Float64, ::Type{metricT}, space::Tuple) where {AD, metricT <: AHMC.AbstractMetric}
    return HMCDA{AD, space, metricT}(n_adapts, δ, λ, ϵ)
end

function HMCDA{AD}(
    δ::Float64,
    λ::Float64;
    init_ϵ::Float64=0.0,
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return HMCDA{AD}(0, δ, λ, init_ϵ, metricT, ())
end

function HMCDA{AD}(
    n_adapts::Int,
    δ::Float64,
    λ::Float64,
    ::Tuple{};
    kwargs...
) where AD
    return HMCDA{AD}(n_adapts, δ, λ; kwargs...)
end

function HMCDA{AD}(
    n_adapts::Int,
    δ::Float64,
    λ::Float64,
    space::Symbol...;
    init_ϵ::Float64=0.0,
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return HMCDA{AD}(n_adapts, δ, λ, init_ϵ, metricT, space)
end

"""
    NUTS(n_adapts::Int, δ::Float64; max_depth::Int=5, Δ_max::Float64=1000.0, ϵ::Float64=0.0)

No-U-Turn Sampler (NUTS) sampler.

Usage:

```julia
NUTS(200, 0.6j_max)
```

Arguments:

- `n_adapts::Int` : The number of samples to use with adaptation.
- `δ::Float64` : Target acceptance rate.
- `max_depth::Float64` : Maximum doubling tree depth.
- `Δ_max::Float64` : Maximum divergence during doubling tree.
- `ϵ::Float64` : Inital step size; 0 means automatically search by Turing.

"""
mutable struct NUTS{AD, space, metricT <: AHMC.AbstractMetric} <: AdaptiveHamiltonian{AD}
    n_adapts    ::  Int         # number of samples with adaption for ϵ
    δ           ::  Float64     # target accept rate
    max_depth   ::  Int         # maximum tree depth
    Δ_max       ::  Float64
    ϵ           ::  Float64     # (initial) step size
end

NUTS(args...; kwargs...) = NUTS{ADBackend()}(args...; kwargs...)
function NUTS{AD}(
    n_adapts::Int,
    δ::Float64,
    max_depth::Int,
    Δ_max::Float64,
    ϵ::Float64,
    ::Type{metricT},
    space::Tuple
) where {AD, metricT}
    return NUTS{AD, space, metricT}(n_adapts, δ, max_depth, Δ_max, ϵ)
end

function NUTS{AD}(
    n_adapts::Int,
    δ::Float64,
    ::Tuple{};
    kwargs...
) where AD
    NUTS{AD}(n_adapts, δ; kwargs...)
end

function NUTS{AD}(
    n_adapts::Int,
    δ::Float64,
    space::Symbol...;
    max_depth::Int=10,
    Δ_max::Float64=1000.0,
    init_ϵ::Float64=0.0,
    metricT=AHMC.DiagEuclideanMetric
) where AD
    NUTS{AD}(n_adapts, δ, max_depth, Δ_max, init_ϵ, metricT, space)
end

function NUTS{AD}(
    δ::Float64;
    max_depth::Int=10,
    Δ_max::Float64=1000.0,
    init_ϵ::Float64=0.0,
    metricT=AHMC.DiagEuclideanMetric
) where AD
    NUTS{AD}(0, δ, max_depth, Δ_max, init_ϵ, metricT, ())
end

function NUTS{AD}(kwargs...) where AD
    NUTS{AD}(0, 0.65; kwargs...)
end

for alg in (:HMC, :HMCDA, :NUTS)
    @eval getmetricT(::$alg{<:Any, <:Any, metricT}) where {metricT} = metricT
end

####
#### Sampler construction
####

# Sampler(alg::Hamiltonian) =  Sampler(alg, AHMCAdaptor())
function Sampler(
    alg::SamplerType,
    model::Model,
    s::Selector=Selector()
) where {SamplerType<:Union{StaticHamiltonian, AdaptiveHamiltonian}}
    info = Dict{Symbol, Any}()
    # Create an empty sampler state that just holds a typed VarInfo.
    initial_state = SamplerState(VarInfo(model))

    # Create an initial sampler, to get all the initialization out of the way.
    initial_spl = Sampler(alg, info, s, initial_state)

    # Create the actual state based on the alg type.
    state = HMCState(model, initial_spl, GLOBAL_RNG)

    # Create a real sampler after getting all the types/running the init phase.
    return Sampler(alg, initial_spl.info, initial_spl.selector, state)
end

####
#### Transition / step functions for HMC samplers.
####

# Single step of a Hamiltonian.
function step!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{T},
    N::Integer;
    kwargs...
) where T<:Hamiltonian
    # Get step size
    ϵ = T <: AdaptiveHamiltonian ?
        AHMC.getϵ(spl.state.adaptor) :
        spl.alg.ϵ

    spl.state.i += 1
    spl.state.eval_num = 0

    Turing.DEBUG && @debug "current ϵ: $ϵ"

    # Gibbs component specified cares
    if spl.selector.tag != :default
        # Transform the space
        Turing.DEBUG && @debug "X-> R..."
        link!(spl.state.vi, spl)
        runmodel!(model, spl.state.vi, spl)
        # Update Hamiltonian
        metric = gen_metric(length(spl.state.vi[spl]), spl)
        ∂logπ∂θ = gen_∂logπ∂θ(spl.state.vi, spl, model)
        logπ = gen_logπ(spl.state.vi, spl, model)
        spl.state.h = AHMC.Hamiltonian(metric, logπ, ∂logπ∂θ)
    end

    # Get position and log density before transition
    θ_old, log_density_old = spl.state.vi[spl], spl.state.vi.logp

    # Transition
    t = AHMC.step(rng, spl.state.h, spl.state.traj, spl.state.z)
    # Update z in state
    spl.state.z = t.z

    # Adaptation
    if T <: AdaptiveHamiltonian
        spl.state.h, spl.state.traj, isadapted =
            AHMC.adapt!(spl.state.h, spl.state.traj, spl.state.adaptor,
                        spl.state.i, spl.alg.n_adapts, t.z.θ, t.stat.acceptance_rate)
    end

    Turing.DEBUG && @debug "decide whether to accept..."

    # Update `vi` based on acceptance
    if t.stat.is_accept
        spl.state.vi[spl] = t.z.θ
        setlogp!(spl.state.vi, t.stat.log_density)
    else
        spl.state.vi[spl] = θ_old
        setlogp!(spl.state.vi, log_density_old)
    end

    # Gibbs component specified cares
    # Transform the space back
    Turing.DEBUG && @debug "R -> X..."
    spl.selector.tag != :default && invlink!(spl.state.vi, spl)

    return HamiltonianTransition(spl, t)
end

#####
##### HMC core functions
#####

"""
    gen_∂logπ∂θ(vi::VarInfo, spl::Sampler, model)

Generate a function that takes a vector of reals `θ` and compute the logpdf and
gradient at `θ` for the model specified by `(vi, spl, model)`.
"""
function gen_∂logπ∂θ(vi::VarInfo, spl::Sampler, model)
    function ∂logπ∂θ(x)
        return gradient_logp(x, vi, model, spl)
    end
    return ∂logπ∂θ
end

"""
    gen_logπ(vi::VarInfo, spl::Sampler, model)

Generate a function that takes `θ` and returns logpdf at `θ` for the model specified by
`(vi, spl, model)`.
"""
function gen_logπ(vi::VarInfo, spl::Sampler, model)
    function logπ(x)::Float64
        x_old, lj_old = vi[spl], vi.logp
        vi[spl] = x
        runmodel!(model, vi, spl)
        lj = vi.logp
        vi[spl] = x_old
        setlogp!(vi, lj_old)
        return lj
    end
    return logπ
end

gen_metric(dim::Int, spl::Sampler{<:Hamiltonian}) = AHMC.UnitEuclideanMetric(dim)
gen_metric(dim::Int, spl::Sampler{<:AdaptiveHamiltonian}) = AHMC.renew(spl.state.h.metric, AHMC.getM⁻¹(spl.state.adaptor.pc))

gen_traj(alg::HMC, ϵ) = AHMC.StaticTrajectory(AHMC.Leapfrog(ϵ), alg.n_leapfrog)
gen_traj(alg::HMCDA, ϵ) = AHMC.HMCDA(AHMC.Leapfrog(ϵ), alg.λ)
gen_traj(alg::NUTS, ϵ) = AHMC.NUTS(AHMC.Leapfrog(ϵ), alg.max_depth, alg.Δ_max)

####
#### Compiler interface, i.e. tilde operators.
####
function assume(spl::Sampler{<:Hamiltonian},
    dist::Distribution,
    vn::VarName,
    vi::VarInfo
)
    Turing.DEBUG && @debug "assuming..."
    updategid!(vi, vn, spl)
    r = vi[vn]
    # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))
    # r
    Turing.DEBUG && @debug "dist = $dist"
    Turing.DEBUG && @debug "vn = $vn"
    Turing.DEBUG && @debug "r = $r" "typeof(r)=$(typeof(r))"
    return r, logpdf_with_trans(dist, r, istrans(vi, vn))
end

function assume(spl::Sampler{<:Hamiltonian},
    dists::Vector{<:Distribution},
    vn::VarName,
    var::Any,
    vi::VarInfo
)
    @assert length(dists) == 1 "[observe] Turing only support vectorizing i.i.d distribution"
    dist = dists[1]
    n = size(var)[end]

    vns = map(i -> VarName(vn, "[$i]"), 1:n)

    rs = vi[vns]  # NOTE: inside Turing the Julia conversion should be sticked to

    # acclogp!(vi, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1]))))

    if isa(dist, UnivariateDistribution) || isa(dist, MatrixDistribution)
        @assert size(var) == size(rs) "Turing.assume variable and random number dimension unmatched"
        var = rs
    elseif isa(dist, MultivariateDistribution)
        if isa(var, Vector)
            @assert length(var) == size(rs)[2] "Turing.assume variable and random number dimension unmatched"
            for i = 1:n
                var[i] = rs[:,i]
            end
        elseif isa(var, Matrix)
            @assert size(var) == size(rs) "Turing.assume variable and random number dimension unmatched"
            var = rs
        else
            error("[Turing] unsupported variable container")
        end
    end

    var, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1])))
end

observe(spl::Sampler{<:Hamiltonian},
    d::Distribution,
    value::Any,
    vi::VarInfo) = observe(nothing, d, value, vi)

observe(spl::Sampler{<:Hamiltonian},
    ds::Vector{<:Distribution},
    value::Any,
    vi::VarInfo) = observe(nothing, ds, value, vi)

####
#### Default HMC stepsize and mass matrix adaptor
####

function AHMCAdaptor(alg::AdaptiveHamiltonian, metric::AHMC.AbstractMetric; ϵ=alg.ϵ)
    pc = AHMC.Preconditioner(metric)
    da = AHMC.NesterovDualAveraging(alg.δ, ϵ)
    if metric == AHMC.UnitEuclideanMetric
        adaptor = AHMC.NaiveHMCAdaptor(pc, da)
    else
        adaptor = AHMC.StanHMCAdaptor(alg.n_adapts, pc, da)
    end
    return adaptor
end

AHMCAdaptor(::Hamiltonian, ::AHMC.AbstractMetric; kwargs...) = AHMC.Adaptation.NoAdaptation()

##########################
# HMC State Constructors #
##########################

function HMCState(
    model::Model,
    spl::Sampler{<:Hamiltonian},
    rng::AbstractRNG;
    kwargs...
)
    # Reuse the VarInfo.
    vi = spl.state.vi

    # Link everything if needed.
    !islinked(vi, spl) && link!(vi, spl)

    # Get the initial log pdf and gradient functions.
    ∂logπ∂θ = gen_∂logπ∂θ(vi, spl, model)
    logπ = gen_logπ(vi, spl, model)

    # Get the metric type.
    metricT = getmetricT(spl.alg)

    # Create a Hamiltonian.
    θ_init = Vector{Float64}(spl.state.vi[spl])
    metric = metricT(length(θ_init))
    h = AHMC.Hamiltonian(metric, logπ, ∂logπ∂θ)

    # Find good eps if not provided one
    if spl.alg.ϵ == 0.0
        ϵ = AHMC.find_good_eps(h, θ_init)
        @info "Found initial step size" ϵ
    else
        ϵ = spl.alg.ϵ
    end

    # Generate a trajectory.
    traj = gen_traj(spl.alg, ϵ)

    # Generate a phasepoint. Replaced during sample_init!
    h, t = AHMC.sample_init(rng, h, θ_init) # this also ensure AHMC has the same dim as θ.

    # Unlink everything.
    invlink!(vi, spl)

    return HMCState(vi, 0, 0, traj, h, AHMCAdaptor(spl.alg, metric; ϵ=ϵ), t.z)
end

#######################################################
# Special callback functionality for the HMC samplers #
#######################################################

mutable struct HMCCallback{
    ProgType<:ProgressMeter.AbstractProgress
} <: AbstractCallback
    p :: ProgType
end

function callback(
    rng::AbstractRNG,
    model::ModelType,
    spl::SamplerType,
    N::Integer,
    iteration::Integer,
    t::HamiltonianTransition,
    cb::HMCCallback;
    kwargs...
) where {
    ModelType<:Sampleable,
    SamplerType<:AbstractSampler
}
    AHMC.pm_next!(cb.p, t.stat, iteration, spl.state.h.metric)
end

function init_callback(
    rng::AbstractRNG,
    model::Model,
    s::Sampler{<:Union{StaticHamiltonian, AdaptiveHamiltonian}},
    N::Integer;
    dt::Real=0.25,
    kwargs...
)
    return HMCCallback(ProgressMeter.Progress(N, dt=dt, desc="Sampling ", barlen=31))
end


text:
### Extensible system for multimedia display of Python objects

const PYSHOW_RULES = Function[]

function pyshow_add_rule(rule::Function)
    push!(PYSHOW_RULES, rule)
    return
end

function pyshow(io::IO, mime::MIME, x)
    x_ = Py(x)
    for rule in PYSHOW_RULES
        rule(io, string(mime), x_) && return
    end
    throw(MethodError(show, (io, mime, x_)))
end

function pyshowable(mime::MIME, x)
    x_ = Py(x)
    for rule in PYSHOW_RULES
        rule(devnull, string(mime), x_) && return true
    end
    return false
end

### Particular rules

# x._repr_mimebundle_()
function pyshow_rule_mimebundle(io::IO, mime::String, x::Py)
    try
        ans = x._repr_mimebundle_(include=pylist((mime,)))
        if pyisinstance(ans, pybuiltins.tuple)
            data = ans[0][mime]
        else
            data = ans[mime]
        end
        data = pyconvert(Union{String,Vector{UInt8}}, data)
        if mime == "text/html"
            data = String(data)
            if occursin("altair-viz-", data) && occursin("document.currentScript.previousElementSibling", data)
                # unsatisfactory hack to get altair plots to display properly in pluto
                # TODO: fix this upstream (in altair or pluto??)
                data = replace(data, "document.currentScript.previousElementSibling" => "((document.currentScript || {}).previousElementSibling || {})")
            end
        end
        write(io, data)
        return true
    catch exc
        if exc isa PyException
            return false
        else
            rethrow()
        end
    end
end

const MIME_TO_REPR_METHOD = Dict(
    "text/plain" => "__repr__",
    "text/html" => "_repr_html_",
    "text/markdown" => "_repr_markdown_",
    "text/json" => "_repr_json_",
    "text/latex" => "_repr_latex_",
    "application/javascript" => "_repr_javascript_",
    "application/pdf" => "_repr_pdf_",
    "image/jpeg" => "_repr_jpeg_",
    "image/png" => "_repr_png_",
    "image/svg+xml" => "_repr_svg_",
)

# x._repr_FORMAT_()
function pyshow_rule_repr(io::IO, mime::String, x::Py)
    method = get(MIME_TO_REPR_METHOD, mime, "")
    isempty(method) && return false
    try
        ans = pygetattr(x, method)()
        if pyisinstance(ans, pybuiltins.tuple)
            data = ans[0]
        else
            data = ans
        end
        write(io, pyconvert(Union{String,Vector{UInt8}}, data))
        return true
    catch exc
        if exc isa PyException
            return false
        else
            rethrow()
        end
    end
end

const MIME_TO_MATPLOTLIB_FORMAT = Dict(
    "image/png" => "png",
    "image/jpeg" => "jpeg",
    "image/tiff" => "tiff",
    "image/svg+xml" => "svg",
    "application/pdf" => "pdf",
)

# x.savefig()
# Requires x to be a matplotlib.pyplot.Figure, or x.figure to be one.
# Closes the underlying figure.
function pyshow_rule_savefig(io::IO, mime::String, x::Py)
    format = get(MIME_TO_MATPLOTLIB_FORMAT, mime, "")
    isempty(format) && return false
    pyhasattr(x, "savefig") || return false
    try
        plt = pysysmodule.modules["matplotlib.pyplot"]
        Figure = plt.Figure
        fig = x
        while !pyisinstance(fig, Figure)
            fig = fig.figure
        end
        buf = pyimport("io").BytesIO()
        x.savefig(buf, format=format)
        data = pyconvert(Vector{UInt8}, buf.getvalue())
        write(io, data)
        plt.close(fig)
        return true
    catch exc
        if exc isa PyException
            return false
        else
            rethrow()
        end
    end
end

function init_pyshow()
    pyshow_add_rule(pyshow_rule_mimebundle)
    pyshow_add_rule(pyshow_rule_repr)
    pyshow_add_rule(pyshow_rule_savefig)
end


text:
abstract type AbstractFeaturedGraph end

"""
    NullGraph()

Null object for `FeaturedGraph`.
"""
struct NullGraph <: AbstractFeaturedGraph end

"""
    FeaturedGraph(graph, feature)

A feature-equipped graph structure for passing graph to layer in order to provide graph dynamically.
References to graph or feature are hold in this type. `graph` and `feature` are provided to fetch data.

# Arguments
- `graph`: should be a adjacency matrix, `SimpleGraph`, `SimpleDiGraph` (from LightGraphs) or `SimpleWeightedGraph`, `SimpleWeightedDiGraph` (from SimpleWeightedGraphs).
- `feature`: features attached to graph.
"""
struct FeaturedGraph{T,S} <: AbstractFeaturedGraph
    graph::Ref{T}
    feature::Ref{S}
    FeaturedGraph(graph::T, feature::S) where {T,S} = new{T,S}(Ref(graph), Ref(feature))
end

"""
    graph(::AbstractFeaturedGraph)

Get referenced graph.
"""
graph(::NullGraph) = nothing
graph(fg::FeaturedGraph) = fg.graph[]

"""
    feature(::AbstractFeaturedGraph)

Get feature attached to graph.
"""
feature(::NullGraph) = nothing
feature(fg::FeaturedGraph) = fg.feature[]

nv(::NullGraph) = 0
nv(fg::FeaturedGraph) = nv(fg.graph[])
nv(fg::FeaturedGraph{T}) where {T<:AbstractMatrix} = size(fg.graph[], 1)

## Linear algebra API for AbstractFeaturedGraph

adjacency_matrix(fg::FeaturedGraph, T::DataType=eltype(fg.graph[])) = adjacency_matrix(fg.graph[], T)

function degrees(fg::FeaturedGraph, T::DataType=eltype(fg.graph[]); dir::Symbol=:out)
    degrees(fg.graph[], T; dir=dir)
end

function degree_matrix(fg::FeaturedGraph, T::DataType=eltype(fg.graph[]); dir::Symbol=:out)
    degree_matrix(fg.graph[], T; dir=dir)
end

function inv_sqrt_degree_matrix(fg::FeaturedGraph, T::DataType=eltype(fg.graph[]); dir::Symbol=:out)
    inv_sqrt_degree_matrix(fg.graph[], T; dir=dir)
end

function laplacian_matrix(fg::FeaturedGraph, T::DataType=eltype(fg.graph[]); dir::Symbol=:out)
    laplacian_matrix(fg.graph[], T; dir=dir)
end

function normalized_laplacian(fg::FeaturedGraph, T::DataType=eltype(fg.graph[]); selfloop::Bool=false)
    normalized_laplacian(fg.graph[], T; selfloop=selfloop)
end

function scaled_laplacian(fg::FeaturedGraph, T::DataType=eltype(fg.graph[]))
    scaled_laplacian(fg.graph[], T)
end


text:
using Balloons
using Test

@testset "Packbits" begin
    include("packbits.jl")
end

@testset "LZW" begin
    include("lzw.jl")
end


text:
module FormFactor

using SpecialFunctions
using HypergeometricFunctions
using WignerFamilies
using FLoops

include("extra_form.jl")

export mform,mfun
@time for n1 ∈ 1:50
    for l1 ∈ 0:n1-1
               for m1 ∈ 0:l1
                   print("$n1 " , "$l1 " , "$m1 " , mform(n1, l1, m1, n1, l1, 1.0, 0.66),"\n")
               end
           end
       end
mform(10, 8, 7, 8, 7, 1.0, 0.66)
@time mform(5,4,3,4,3,1.0,0.66)
end


text:
const epsilon       = 18.016/28.966
const g             = 10u"m/s/s" #acceleration of gravity

struct Substance{T}
    cp :: Union{Nothing,T}
    cv :: Union{Nothing,T}
    R  :: Union{Nothing,T}
    Lv :: Union{Nothing,T}
    Lf :: Union{Nothing,T}
end

Substance{T}(;cp = nothing, cv = nothing, R = nothing, Lv = nothing, Lf = nothing) where T = Substance{T}(cp,cv,R,Lv,Lf)

const Dryair = Substance{Quantity}(
    cp = 1006.0u"J/kg/K", #J/kg/k at 1013 hPa
    cv = 718.0u"J/kg/K",
    R  = 287.05u"J/kg/K" # J/kg/k
)

const Liquidwater = Substance{Quantity}(
     Lv = 2.5e6u"J/kg", #J/kg
     Lf = 3.33e5u"J/kg",
     cp = 4200.0u"J/kg/K" #j/kg/k
)

const Watervapor = Substance{Quantity}(
    R = 461.52u"J/kg/K" #j/kg/K
)


text:
using ChainForth
using Test

function pushnums(forth)
    interpret(forth, "21")
    @test forth.stack[end] == 21

    interpret(forth, "10")
    @test forth.stack[end] == 10
end

@testset "Interpreter basics" begin
    forth = interpreter()
    pushnums(forth)
    interpret(forth, "+ \n")
    @test length(forth.stack) == 1
    @test forth.stack[end] == 31

    out = IOBuffer()
    forth = interpreter(stdin, out)
    interpret(forth, "5 6 +")
    @test forth.stack[end] == 11

    interpret(forth, ".")
    @test length(forth.stack) == 0
    @test String(take!(out)) == "11"
end

@testset "Defining with : and ;" begin
    forth = interpreter()
    interpret(forth, ": 5* 5 * ;")
    pushnums(forth)
    interpret(forth, "5*")
    @test length(forth.stack) == 2
    @test forth.stack[end] == 50
    interpret(forth, ": x 1 5* ; x")
    @test length(forth.stack) == 3
    @test forth.stack[end] == 5

    interpret(forth, ": double dup + ;")
    interpret(forth, ": quadruple double double ;")
    interpret(forth, "3 5* quadruple")
    @test forth.stack[end] == 60
end

@testset "Immediate words" begin
    forth = interpreter()
    interpret(forth, ": push5 5 ; immediate")
    interpret(forth, ": x 1 push5 3 ;")
    @test length(forth.stack) == 1
    @test forth.stack[1] == 5
    interpret(forth, "x")
    @test length(forth.stack) == 3
    @test forth.stack[2] == 1
    @test forth.stack[3] == 3
end

@testset "Stack manipulation" begin
    forth = interpreter()
    interpret(forth, "1 2 3 dup")
    @test length(forth.stack) == 4
    @test forth.stack[end] == 3
    interpret(forth, "4 drop")
    @test length(forth.stack) == 4
    @test forth.stack[end] == 3
    interpret(forth, "drop drop")
    @test length(forth.stack) == 2
    @test forth.stack[end] == 2
    interpret(forth, "9 rot")
    @test length(forth.stack) == 3
    @test forth.stack[end] == 1
    @test forth.stack[end-1] == 9
    @test forth.stack[end-2] == 2
end

@testset "Memory Allocation" begin
    forth = interpreter()
    interpret(forth, "5 allot here")
    @test forth.stack == [1]
    interpret(forth, "3 allot here")
    @test forth.stack == [1, 6]
    interpret(forth, "100 allot here")
    @test forth.stack == [1, 6, 9]
    @test length(forth.memory) == 108

    interpret(forth, "12 5 ! 13 8 !")
    @test forth.memory[5] == 12
    @test forth.memory[8] == 13
    @test forth.stack == [1, 6, 9]
    interpret(forth, "5 @ 8 @ 5 @")
    @test forth.stack == [1, 6, 9, 12, 13, 12]
end

@testset "Variables" begin
    forth = interpreter()
    interpret(forth, "variable x")
    interpret(forth, "variable y")
    interpret(forth, "variable zuzu1234567890")
    @test forth.stack == []
    interpret(forth, "9 zuzu1234567890 !")
    @test forth.stack == []
    interpret(forth, "zuzu1234567890 @")
    @test forth.stack == [9]
    interpret(forth, "drop 8 x ! 7 y !")
    interpret(forth, "x @ y @ zuzu1234567890 @")
    @test forth.stack == [8, 7, 9]
end

@testset "Execute" begin
    forth = interpreter()
    interpret(forth, ": GT1 123 ;")
    interpret(forth, ": GT2 ' ' , ;")
    interpret(forth, "GT2 GT1 execute")
    @test forth.stack == [123]
end

@testset "if. then, else" begin
    forth = interpreter()
    interpret(forth, ": t1 5 < if 11 else 22 then ;")
    interpret(forth, "4 t1 5 t1 6 t1")
    @test forth.stack == [11, 22, 22]
    interpret(forth, "drop drop drop")

    interpret(forth, ": t2 dup 5 < if drop 11 else 8 < if 22 else 33 then then ;")
    interpret(forth, "see t2")
    interpret(forth, "4 t2 6 t2 9 t2")
    @test forth.stack == [11, 22, 33]
end